text:
import StatsFuns.chisqccdf

## generic goodness of fit testing

function discreteProbabilities(pmf::F, low::Int64,
  high::Int64) where F <: Function
  n::Int64 = high - low + 2
  probs::Vector{Float64} = Vector{Float64}(undef, n)
  for k in low:high
    probs[k-low+1] = pmf(k)
  end
  probs[n] = 1 - sum(probs[1:n-1])
  return probs
end

function continuousProbabilities(cdf::F, low::Float64, high::Float64,
  h::Float64) where F <: Function
  divisions::Vector{Float64} = Vector(low:h:high)
  n::Int64 = length(divisions)
  probs::Vector{Float64} = Vector{Float64}(undef, n)
  for i in 1:n-1
    probs[i] = cdf(divisions[i+1]) - cdf(divisions[i])
  end
  probs[n] = 1 - sum(probs[1:n-1])
  return probs
end

function testGOFMultinomial(probs::Vector{Float64}, counts::Vector{Int64})
  n::Int64 = sum(counts)
  es::Vector{Float64} = n .* probs
  if minimum(es) < 10.0
    println("warning: minimum expected value of ", minimum(es),
      " in bin ", indmin(es), " of ", n)
  end
  testStatistic::Float64 = sum((counts .- es).^2 ./ es)
  pValue::Float64 = chisqccdf(length(probs)-1, testStatistic)
  return pValue
end

function testGOFDiscrete(pmf::F1, sampler::F2, low::Int64, high::Int64,
  N::Int64) where {F1 <: Function, F2 <: Function}
  n = high - low + 2
  probs::Vector{Float64} = discreteProbabilities(pmf, low, high)
  counts::Vector{Int64} = Vector{Int64}(undef, n)
  fill!(counts, 0)
  for i in 1:N
    v::Int64 = sampler()
    if v < low || v > high
      counts[n] += 1
    else
      counts[v-low+1] += 1
    end
  end
  return testGOFMultinomial(probs, counts)
end

function testGOFContinuous(cdf::F1, sampler::F2, low::Float64, high::Float64,
  h::Float64, N::Int64) where {F1 <: Function, F2 <: Function}
  divisions::Vector{Float64} = Vector(low:h:high)
  n::Int64 = length(divisions)
  probs::Vector{Float64} = continuousProbabilities(cdf, low, high, h)
  counts::Vector{Int64} = Vector{Int64}(undef, n)
  fill!(counts, 0)
  for i in 1:N
    v::Float64 = sampler()
    if v < divisions[1] || v > divisions[n]
      counts[n] += 1
    else
      bin::Int64 = floor(Int64, (v - divisions[1])/h) + 1
      counts[bin] += 1
    end
  end
  return testGOFMultinomial(probs, counts)
end


text:
givensids(P) = [ (idx1, idx2) for idx1 = 1:P-1 for idx2 = (idx1+1):P ]

ngivensangles(n::Integer) = fld(n*(n-1),2)

mat2rotations(args...) = mat2rotations_lowmemory(args...)

function mat2rotations_reference(mtx::AbstractMatrix{T}) where T <: Real
    P = size(mtx, 1)

    ids = givensids(P)
    res = similar(mtx, length(ids))

    for (nr, (idx1, idx2)) in enumerate(ids)
        a = givens(mtx, idx1, idx2, idx1)
        g = a[1]

        res[nr] = atan(g.s, g.c)

        R = Matrix(I,P,P)
        R[g.i1, g.i1] =  g.c
        R[g.i1, g.i2] =  g.s
        R[g.i2, g.i1] = -g.s
        R[g.i2, g.i2] =  g.c

        mtx .= R*mtx
    end
    (res, round.(diag(mtx)))
end

function mat2rotations_lowmemory(mtx::AbstractMatrix{T}) where {T<:Real}
    ids = givensids(size(mtx, 1))
    res = similar(mtx, length(ids))

    mtx = Array(mtx)
    for (nr, (idx1, idx2)) in enumerate(ids)
        a = givens(mtx, idx1, idx2, idx1)
        g = a[1]

        res[nr] = atan(g.s, g.c)

        row1 = mtx[g.i1, :]
        vrw2 = @view mtx[g.i2, :]
        mtx[g.i1, :] =  g.c * row1 + g.s * vrw2
        mtx[g.i2, :] = -g.s * row1 + g.c * vrw2
    end
    (res, round.(diag(mtx)))
end

function rotations2mat(θs::AbstractArray, sig::AbstractArray)
    rotations2mat(θs, sig, round(Integer, (1 + sqrt(1+8*length(θs))) / 2))
end

rotations2mat(θs, sig, P) = rotations2mat_lowmemory(θs, sig, P)

function rotations2mat_reference(θs::AbstractArray{TA}, sig::AbstractArray{TS}, P::Integer) where {TA<:Real,TS<:Number}
    mtx = Matrix{TA}(I,P,P)

    for (nr, (idx1, idx2)) in enumerate(givensids(P))
        c, s = cos(θs[nr]), sin(θs[nr])

        R = Matrix(I,P,P)
        R[idx1, idx1] =  c
        R[idx1, idx2] = -s
        R[idx2, idx1] =  s
        R[idx2, idx2] =  c

        mtx .= mtx*R
    end
    mtx * diagm(0 => sig)
end

function rotations2mat_lowmemory(θs::AbstractArray{TA}, sig::AbstractArray{TS}, P::Integer) where {TA<:Real,TS<:Number}
    mtx = Matrix{TA}(I,P,P)

    for (nr, (idx1, idx2)) in enumerate(givensids(P))
        s, c = sincos(θs[nr])

        col1 = mtx[:, idx1]
        vcl2 = @view mtx[:, idx2]

        mtx[:, idx1] =  c * col1 + s * vcl2
        mtx[:, idx2] = -s * col1 + c * vcl2
    end
    mtx * diagm(0 => sig)
end

graph_one_factorization(P::Integer) = graph_one_factorization(Val(isodd(P)), P)

# When P is odd.
function graph_one_factorization(::Val{true}, P::Integer)
    map(0:P-1) do pivot_idx
        map(1:fld(P,2)) do offset
            ids = [ pivot_idx - offset, pivot_idx + offset ]
            (sort(mod.(ids, P) .+ 1)...,)
        end
    end
end

# When P is even.
function graph_one_factorization(::Val{false}, P::Integer)
    oddgf = graph_one_factorization(P-1)
    [ [ elem..., (idx, P) ] for (idx, elem) in enumerate(oddgf) ]
end

# function mat2rotations_permutated(mtx::AbstractMatrix{T}) where T <: Real
#     P = size(mtx, 1)
#
#     res = similar(mtx, fld(P*(P-1),2))
#     mtx = Array(mtx)
#
#     for (nr, (idx1, idx2)) in enumerate(vcat(graph_one_factorization(P)...))
#         a = givens(mtx, idx1, idx2, idx1)
#         g = a[1]
#
#         res[nr] = atan(g.s, g.c)
#
#         R = Matrix{T}(I,P,P)
#         R[g.i1, g.i1] =  g.c
#         R[g.i1, g.i2] =  g.s
#         R[g.i2, g.i1] = -g.s
#         R[g.i2, g.i2] =  g.c
#         mtx .= R*mtx
#     end
#     (res, round.(diag(mtx)))
# end
#
# function rotations2mat_permutated(θs::AbstractArray{TA}, sig::AbstractArray{TS}, P::Integer) where {TA<:Real,TS<:Number}
#     θrsp = reshape(θs, 2*cld(P,2)-1, fld(P,2))
#     mtx = mapreduce(*, enumerate(graph_one_factorization(P))) do (blk, ids)
#         R = Matrix{TA}(I,P,P)
#         for (nr, (idx1, idx2)) in enumerate(ids)
#             s = sin(θrsp[blk, nr])
#             c = cos(θrsp[blk, nr])
#             R[idx1, idx1] =  c
#             R[idx1, idx2] = -s
#             R[idx2, idx1] =  s
#             R[idx2, idx2] =  c
#         end
#         R
#     end
#     mtx * diagm(0 => sig)
# end

ext: jl

max_stars_repo_name: nixir/MDCDL

max_issues_repo_name: nixir/MDCDL

max_forks_repo_name: nixir/MDCDL


text:
export readplist_string, readplist

function parse_obj(p::Parser)
    token = next_token(p)

    if     NUMBER == token.kind
        Meta.parse(token.lexeme)
    elseif STRING == token.kind || IDENT == token.kind
        token.lexeme
    elseif HEXBINARY == token.kind
        parse_hexbinary(token.lexeme)
    elseif LPAREN == token.kind
        backup_token(p)
        parse_array(p)
    elseif LBRACE == token.kind
        backup_token(p)
        parse_dict(p)
    end
end

parse_hexbinary(s::AbstractString) = hex2bytes(filter(isxdigit, s))

function parse_array(p::Parser)
    array = []

    expect(p, LPAREN)
    token = peek_token(p)
    if token.kind == RPAREN
        return array
    end
    push!(array, parse_obj(p))
    token = peek_token(p)
    while token.kind != RPAREN
        expect(p, COMMA)
        push!(array, parse_obj(p))
        token = peek_token(p)
    end
    expect(p, RPAREN)
    array
end

function parse_dict(p::Parser)
    dict =  Dict{Any, Any}()

    expect(p, LBRACE)
    token = peek_token(p)
    while token.kind != RBRACE
        key = parse_obj(p)::String
        expect(p, EQUAL)
        dict[key] = parse_obj(p)
        expect(p, SEMICOLON)
        token = peek_token(p)
    end
    expect(p, RBRACE)
    dict
end

"""
    readplist_string(text::AbstractString) -> Dict

Parse a NeXTSTEP style property list stored in the string `text` and return result
as a dictionary object.
"""
function readplist_string(text::AbstractString)
    l = lex_plist(text)
    p = Parser(l)
    parse_obj(p)
end

"""
    readplist(stream::IO) -> Dict
    readplist(filename::AbstractString) -> Dict

Read a NeXTSTEP style PList from file named `filename`, or from an `IO` stream object
you have opened previously. This could be from a file, a string or a socket. In each case
you get a dictionary back.
"""
function readplist(stream::IO)
    text = read(stream, String)
    readplist_string(text)
end

function readplist(filename::AbstractString)
    open(readplist, filename)
end


text:
using PriceIndex
using Test

@testset "PriceIndex.jl" begin

@test foo() == 0.11388071406436832
@test foo(1, 1.5) == 0.2731856314283442
@test_broken foo(1, 0) # tells us this is broken Write your own tests here.

end


text:
"""
    function dl_pfdsj()
---
## Prepare a template
```bash
pandoc -D latex >template.tex
```
Insert
```latex
\\XeTeXlinebreaklocale "zh"
```
after
```latex
\$if(mathfont)\$
% blah blah
\$endif\$

## Compile to pdf
```bash
pandoc pf.md
    -V 'mainfont:FangSong'
    -V papersize:a5
    --pdf-engine=xelatex
    --template=template.tex
    -o pf.pdf
```
"""
function dl_pfdsj()
    url = "http://www.dushu369.com/zhongguomingzhu/HTML"
    pf = open("pf.md", "w")
    for page in 12349:-1:12188
        page == 12213 && continue
        download("$url/$page.html", "tmp.txt")
        run(`iconv -cf gbk -t utf8 tmp.txt -o page.txt`)
        chapter = parsehtml(read("page.txt", String))
        for t in eachmatch(sel"title", chapter.root)
            println(pf, "## ", text(t), "\n\n")
        end

        for c in eachmatch(Selector(".content"), chapter.root)
            println(pf, text(c), "\n")
        end
    end
    close(pf)
end


text:
export Tropical, TropicalF64, TropicalF32, TropicalF16, content

neginf(::Type{T}) where T = typemin(T)
neginf(::Type{T}) where T<:AbstractFloat = typemin(T)
neginf(::Type{T}) where T<:Rational = typemin(T)
neginf(::Type{T}) where T<:Integer = T(-999999)
neginf(::Type{Int16}) = Int16(-16384)
neginf(::Type{Int8}) = Int8(-64)

struct Tropical{T} <: Number
    n::T
    Tropical{T}(x) where T = new{T}(T(x))
    function Tropical(x::T) where T
        new{T}(x)
    end
    function Tropical{T}(x::Tropical{T}) where T
        x
    end
    function Tropical{T1}(x::Tropical{T2}) where {T1,T2}
        new{T1}(T2(x.n))
    end
end

Base.show(io::IO, t::Tropical) = Base.print(io, "$(t.n)ₜ")

Base.:*(a::Tropical, b::Tropical) = Tropical(a.n + b.n)
function Base.:*(a::Tropical{<:Rational}, b::Tropical{<:Rational})
    if a.n.den == 0
        a
    elseif b.n.den == 0
        b
    else
        Tropical(a.n + b.n)
    end
end
Base.:+(a::Tropical, b::Tropical) = Tropical(max(a.n, b.n))
Base.zero(::Type{Tropical{T}}) where T = Tropical(neginf(T))
Base.zero(::Tropical{T}) where T = zero(Tropical{T})

Base.one(::Type{Tropical{T}}) where T = Tropical(zero(T))
Base.one(::Tropical{T}) where T = one(Tropical{T})

# inverse and division
Base.inv(x::Tropical) = Tropical(-x.n)
Base.:/(x::Tropical, y::Tropical) = Tropical(x.n - y.n)
Base.div(x::Tropical, y::Tropical) = Tropical(x.n - y.n)

Base.isapprox(x::Tropical, y::Tropical; kwargs...) = isapprox(x.n, y.n; kwargs...)

# promotion rules
Base.promote_type(::Type{Tropical{T1}}, b::Type{Tropical{T2}}) where {T1, T2} = Tropical{promote_type(T1,T2)}


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    BinomialProcess(n)

A Binomial point process with `n` points.
"""
struct BinomialProcess <: PointProcess
  n::Int
end

ishomogeneous(p::BinomialProcess) = true

default_sampling_algorithm(::BinomialProcess, ::Geometry) = ProductSampling()

function rand_single(rng::Random.AbstractRNG,
                     p::BinomialProcess, b::Box{Dim,T},
                     ::ProductSampling) where {Dim,T}
  # region configuration
  lo, up = coordinates.(extrema(b))

  # product of uniform distributions
  U = product_distribution([Uniform(lo[i], up[i]) for i in 1:Dim])

  # return point pattern
  PointSet(rand(rng, U, p.n))
end


text:
using Omega
using Test
using Statistics: mean

function testrid()
  θ = betarv(2.0, 2.0)
  x = bernoulli(θ)
  ridxθ = Omega.rid(x, θ)
  mean1, mean2 = rand((θ, lmean(ridxθ)))
  @test mean1 == mean2
  ω1 = Omega.rand(defΩ())
  ω2 = Omega.rand(defΩ())
  x_ = ridxθ(ω1)
  θ_ = θ(ω1)
  nsamples = 10000
  θ_approx = mean([x_(Omega.rand(defΩ())) for i = 1:nsamples])
  @test isapprox(θ_approx, θ_; atol = 0.01)
end

testrid()

function testrid2()
  t1 = normal(0.0, 1.0)
  t2 = normal(0.0, 1.0)
  x = normal(t1 + t2, 1.0)
  xrid = rid(x, t1, t2)
  samples = rand(t1, lmean(xrid) ==ₛ 0.0, 1000; alg = SSMH)
end

testrid2()


text:
function create_quadratic_spline(x::Union{Array{DateTime,1},Array{Date,1},Array{Union{Date,DateTime},1}},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve), left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing)
    x_as_Floats = years_from_global_base.(x)
    return create_quadratic_spline(x_as_Floats, y; gradients = gradients, extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
end

function create_quadratic_spline(x::Array{<:Real,1},y::Array{<:Real,1} ; gradients::Union{Missing,Array{<:Real,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve),
                                 left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing)
    schum = Schumaker(x, y; gradients = gradients, extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
    return create_quadratic_spline(schum)
end

function create_quadratic_spline(schum::Schumaker)
    starts_ = schum.IntStarts_
    coefficients = schum.coefficient_matrix_
    number_of_intervals = size(coefficients)[1]
    funcs_ = Array{Sum_Of_Functions,1}(undef, number_of_intervals)
    for i in 1:number_of_intervals
        quadratic = PE_Function(coefficients[i,1], 0.0, starts_[i], 2)
        linear    = PE_Function(coefficients[i,2], 0.0, starts_[i], 1)
        constant  = PE_Function(coefficients[i,3], 0.0, 0.0       , 0)
        polynomial = Sum_Of_Functions([quadratic, linear, constant])
        funcs_[i] = polynomial
    end
    return Piecewise_Function(starts_, funcs_)
end

function create_quadratic_spline(x::Union{Array{D,1},Array{<:DatePeriod,1}},y::Array{<:Real,1}; gradients::Union{Missing,Array{<:Real,1}} = missing, extrapolation::Tuple{Schumaker_ExtrapolationSchemes,Schumaker_ExtrapolationSchemes} = (Curve,Curve),
                                 left_gradient::Union{Missing,Real} = missing, right_gradient::Union{Missing,Real} = missing) where D<:DatePeriod
    x_as_Floats = period_length.(x)
    return create_quadratic_spline(x_as_Floats, y; gradients = gradients, extrapolation = extrapolation, left_gradient = left_gradient, right_gradient = right_gradient)
end

function create_constant_interpolation_to_right(x::Array{Date,1},y::Array{<:Real,1})
    x_Float = years_from_global_base.(x)
    return create_constant_interpolation_to_right(x_Float,y)
end

function create_constant_interpolation_to_right(x::Array{<:Real,1},y::Array{<:Real,1})
    x_ = vcat(-Inf,x)
    y = vcat(y[1], y)
    funcs_ = PE_Function.(y,0.0,0.0,0)
    return Piecewise_Function(x_, funcs_)
end

function create_constant_interpolation_to_right(x::Union{Array{D,1},Array{<:DatePeriod,1}},y::Array{<:Real,1}) where D<:DatePeriod
    x_as_Floats = period_length.(x)
    return create_constant_interpolation_to_right(x_as_Floats,y)
end

function create_constant_interpolation_to_left(x::Array{Date,1},y::Array{<:Real,1})
    x_Float = years_from_global_base.(x)
    return create_constant_interpolation_to_left(x_Float,y)
end

function create_constant_interpolation_to_left(x::Array{<:Real,1},y::Array{<:Real,1})
    x_ = vcat(-Inf,x[1:(length(x)-1)])
    funcs_ = PE_Function.(y,0.0,0.0,0)
    return Piecewise_Function(x_, funcs_)
end

function create_constant_interpolation_to_left(x::Union{Array{D,1},Array{<:DatePeriod,1}},y::Array{<:Real,1}) where D<:DatePeriod
    x_as_Floats = period_length.(x)
    return create_constant_interpolation_to_left(x_as_Floats,y)
end

function create_linear_interpolation(x::Array{Date,1},y::Array{<:Real,1})
    x_Float = years_from_global_base.(x)
    return create_linear_interpolation(x_Float,y)
end

function create_linear_interpolation(x::Union{Array{D,1},Array{<:DatePeriod,1}},y::Array{<:Real,1}) where D<:DatePeriod
    x_as_Floats = period_length.(x)
    return create_linear_interpolation(x_as_Floats,y)
end

function create_linear_interpolation(x::Array{R,1},y::Array{<:Real,1}) where R<:Real
    len = length(x)
    if len < 2
        error("Insufficient data to linearly interpolate")
    end
    starts_ = Array{R,1}(undef, len-1)
    funcs_  = Array{UnivariateFunction}(undef, len-1)
    coefficient = (y[2] - y[1])/(x[2] - x[1])
    starts_[1] = -Inf
    funcs_[1]  = Sum_Of_Functions([PE_Function(y[1],0.0,0.0,0), PE_Function(coefficient,0.0,x[1],1)])
    if len > 2
        for i in 2:(len-1)
            starts_[i] = x[i]
            coefficient = (y[i+1] - y[i])/(x[i+1] - x[i])
            funcs_[i]  = Sum_Of_Functions([PE_Function(y[i],0.0,0.0,0), PE_Function(coefficient,0.0,x[i],1)])
        end
    end
    return Piecewise_Function(starts_, funcs_)
end


text:
using Test
using QuantumLattices.Prerequisites.NamedVectors
using QuantumLattices.Prerequisites: Float

@namedvector mutable struct NHNV
    scope::String
    site::Int
end

@testset "NHNV" begin
    @test fieldnames(NHNV) == (:scope, :site)
    @test length(NHNV) == 2

    pid = NHNV("A", 0)
    @test pid == convert(NHNV, ("A", 0))
    @test string(pid) == "NHNV(\"A\", 0)"
    @test length(pid) == 2
    @test collect(Iterators.reverse(pid)) == [0, "A"]
    @test pid[1] == pid.scope == "A"
    @test pid[2] == pid.site == 0
    @test replace(pid, scope="B") == NHNV("B", 0)
    @test replace(pid, site=1) == NHNV("A", 1)
    @test replace(pid, scope="B", site=1) == NHNV("B", 1)
    @test (pid[1] = "B"; pid[1] == "B")
    @test (pid.site = 2; pid.site == 2)

    @test NHNV("A", 2.0) < NHNV("B", 0.0)
    @test NHNV("A", 2.0) < NHNV("A", 3.0)
    @test isless(NHNV("A", 2.0), NHNV("B", 0.0))
    @test isless(NHNV("A", 2.0), NHNV("A", 3.0))

    dict = Dict(NHNV("A", 0)=>1, NHNV("A", 1)=>2)
    @test dict[NHNV("A", 0)] == 1
    @test dict[NHNV("A", 1)] == 2
end

@homonamedvector "FHNV" (:scope, :site) Float64 mutable = true
@homonamedvector "RHNV" (:scope, :site) (<:Real) mutable = true

@testset "FHNV" begin
    @test fieldnames(FHNV) == (:scope, :site)
    @test length(FHNV) == 2
    @test eltype(FHNV) == Float64
    @test zero(FHNV) == FHNV(0.0, 0.0)
    @test isequal(zero(FHNV), FHNV(0.0, 0.0))

    pid = FHNV(1.0, 0.0)
    @test pid == convert(FHNV, (1.0, 0.0))
    @test string(pid) == "FHNV(1.0, 0.0)"
    @test length(pid) == 2
    @test eltype(pid) == Float64
    @test zero(pid) == FHNV(0.0, 0.0)
    @test pid[1] == pid.scope == 1.0
    @test pid[2] == pid.site == 0.0
    @test replace(pid, scope=2.0) == FHNV(2.0, 0.0)
    @test replace(pid, site=1.0) == FHNV(1.0, 1.0)
    @test replace(pid, scope=2.0, site=1.0) == FHNV(2.0, 1.0)
    @test (pid[1] = 2.0; pid[1] == 2.0)
    @test (pid.site = 3.0; pid.site == 3.0)

    @test FHNV(1.0, 2.0) < FHNV(2.0, 0.0)
    @test FHNV(1.0, 2.0) < FHNV(1.0, 3.0)
    @test isless(FHNV(1.0, 2.0), FHNV(2.0, 0.0))
    @test isless(FHNV(1.0, 2.0), FHNV(1.0, 3.0))

    dict = Dict(FHNV(0.0, 0.0)=>1, FHNV(0.0, 1.0)=>2)
    @test dict[FHNV(0.0, 0.0)] == 1
    @test dict[FHNV(0.0, 1.0)] == 2
end

@testset "RHNV" begin
    @test fieldnames(RHNV) == (:scope, :site)
    @test length(RHNV) == 2
    @test eltype(RHNV{Int}) == Int
    @test zero(RHNV{Int}) == RHNV(0, 0)
    @test eltype(RHNV{Float64}) == Float64
    @test zero(RHNV{Float64}) == RHNV(0.0, 0.0)

    pid = RHNV(1, 0)
    @test pid == convert(RHNV{Int}, (1, 0))
    @test string(pid) == "RHNV(1, 0)"
    @test length(pid) == 2
    @test eltype(pid) == Int
    @test zero(pid) == RHNV(0, 0)
    @test pid[1] == pid.scope == 1
    @test pid[2] == pid.site == 0
    @test replace(pid, scope=2) == FHNV(2, 0)
    @test replace(pid, site=1) == FHNV(1, 1)
    @test replace(pid, scope=2, site=1) == FHNV(2, 1)
    @test (pid[1] = 2; pid[1] == 2)
    @test (pid.site = 3; pid.site == 3)

    @test RHNV(1.0, 2.0) < RHNV(2.0, 0.0)
    @test RHNV(1.0, 2.0) < RHNV(1.0, 3.0)

    dict = Dict(RHNV(0, 0)=>1, RHNV(0, 1)=>2)
    @test dict[RHNV(0, 0)] == 1
    @test dict[RHNV(0, 1)] == 2
end


text:
@testset "test renderer" begin

    iob = IOBuffer()
    Presentation.render(iob, joinpath(@__DIR__, "sample.md"))

    @test String(take!(iob)) == "\e[1mPresentation.jl\e[22m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━┃┃┃┃┏┓┗┛"

end


text:
struct GeoTiff
    raster::AbstractRaster
end
GeoTiff(path::AbstractString) = GeoTiff(_to_rgb(path))
Base.parent(geotiff::GeoTiff) = geotiff.raster

function _to_rgb(fn)
    rast = Raster(fn)
    rast = reverse(rast, dims = 2)
    red = view(rast, Rasters.Band(1)) |> read
    green = view(rast, Rasters.Band(2)) |> read
    blue = view(rast, Rasters.Band(3)) |> read

    replace!.(x -> isnan(x) ? 0.0 : x, [red, green, blue])

    rgb = fill(RGBf(0,0,0), dims(rast, :X), dims(rast, :Y));

    for i in eachindex(red)
        rgb[i] = RGBf(red[i], green[i], blue[i])
    end
    eltype(rast) == UInt8 ? rgb./255 : rgb
end


text:
# The following is taken from StatsBase
# https://github.com/JuliaStats/StatsBase.jl/blob/master/src/sampling.jl
# which is available under the following MIT license:

# Copyright (c) 2012-2016: Dahua Lin, Simon Byrne, Andreas Noack, Douglas Bates, John Myles White, Simon Kornblith, and other contributors.
# Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

function fisher_yates_sample!(rng::AbstractRNG, a::AbstractArray, x::AbstractArray)
    n = length(a)
    k = length(x)
    k <= n || error("length(x) should not exceed length(a)")

    inds = Vector{Int}(undef, n)
    for i = 1:n
        @inbounds inds[i] = i
    end

    @inbounds for i = 1:k
        j = rand(rng, i:n)
        t = inds[j]
        inds[j] = inds[i]
        inds[i] = t
        x[i] = a[t]
    end
    return x
end
fisher_yates_sample!(a::AbstractArray, x::AbstractArray) =
    fisher_yates_sample!(Random.GLOBAL_RNG, a, x)


text:
function layout_choice(n)
  return (n,1)
end

@recipe function plot(sim::Simulation; vars = observables(sim))
  @assert !isempty(vals(sim)) "Results don't contain output. You should probably add output observables to your model"

  time = times(sim)
  simvals = [sim[id,:] for id in vars]

  xguide --> "time"
  yguide --> "output"
  legend --> :right
  label --> permutedims(string.(vars))
  xlims --> (time[1],time[end])
  linewidth --> 3
  (time, simvals)
end

@recipe function plot(sr::SimResult; vars=observables(sr), show_measurements=true)

  # this code should be replaced with
  #   @series begin
  #     sr.sim
  #   end
  # when we separate plots for sim and measurements

  @assert !isempty(vals(sr)) "Results don't contain output. You should probably add output observables to your model"

  #color_choice = theme_palette(:auto)

  time = times(sr)
  for (i,v) in enumerate(vars)
    @series begin
      xguide --> "time"
      yguide --> "output"
      label --> String(v)
      seriescolor --> i
      #seriescolor --> color_choice[i]
      xlims --> (time[1],time[end])
      linewidth --> 3
      (time, sr[v,:])
    end
  end

  # todo: separate plots for measurement tables
  if show_measurements && !isempty(measurements(sr))
    t_meas = NamedTuple{Tuple(vars)}([Float64[] for i in eachindex(vars)])
    vals_meas = NamedTuple{Tuple(vars)}([Float64[] for i in eachindex(vars)])
    for meas in measurements(sr)
      μ = meas.μ
      if isa(μ,Symbol) && μ ∈ vars
        push!(t_meas[μ], meas.t)
        push!(vals_meas[μ], meas.val)
      end
    end
    for (i,v) in enumerate(vars)
      if !isempty(t_meas[v])
        @series begin
          seriestype --> :scatter
          xguide --> "time"
          yguide --> "output"
          label --> false
          seriescolor --> i
          #seriescolor --> color_choice[i]
          (t_meas[v], vals_meas[v])
        end
      end
    end
  end
  nothing
end
#= XXX: do we need it?
@recipe function plot(sim::Vector{S}) where S <:AbstractResult
  [Symbol("_$i")=>s for (i,s) in enumerate(sim)]
end
=#
@recipe function plot(s::Pair{Symbol,S}) where S<:AbstractResult
  @series begin
    title := "$(first(s))"
    last(s)
  end
end

@recipe function plot(sim::Vector{Pair{Symbol,S}}) where S<:AbstractResult
  (m,n) = layout_choice(length(sim))
  layout := (m,n)
  size := (400*n,300*m)

  for (i, s) in enumerate(sim)
    @series begin
      subplot := i
      s
    end
  end
end

#https://github.com/SciML/SciMLBase.jl/blob/7151bbe784df70cc572073d76d3a818aa8d1f4d0/src/ensemble/ensemble_solutions.jl#L102
@recipe function plot(sol::MCResult)
  for i in 1:length(sol)
    @series begin
      legend := false
      sol[i]
    end
  end
end


text:
# Script for further analysis of simulation results.
# Diamagnetic velocity, plasma beta
#
# Hongyang Zhou, hyzhou@umich.edu 10/01/2019

using Batsrus, VisAna, PyPlot, PyCall, Glob, Printf

include("constants.jl")

"Series of contour plots from box output."
function plot_contour_from_box(filename::String, dir=".")

   data = readdata(filename, dir=dir, verbose=false);
   npict = data.list.npictinfiles

   saveDir = "ycutsFromBox"
   mkdir(saveDir)

   for ipict = 1:npict
      data = readdata(filename, npict=ipict, dir=dir)

      plotdata(data, "by bx;bz",plotmode="contbar streamover",
         cut="y",cutPlaneIndex=3, density=1.5)

      fig = matplotlib.pyplot.gcf()
      fig.set_size_inches(3.5,6)
      plt.axis("scaled")

      plt.savefig(saveDir*"/y_"*string(ipict)*".png")
      close()
   end

end

"Series of contour plots from cuts."
function plot_contour_from_cuts(filename::String, dir=".")

   data = readdata(filename, dir=dir)
   npict = data.list.npictinfiles
   saveDir = "ycuts"
   mkdir(saveDir)

   for ipict = 1:npict
      data = readdata(filename,npict=ipict,dir=dir)

      plotdata(data, "by bx;bz", plotmode="contbar streamover",
         plotrange=[1.0 1.8 -2.0 2.0], density=1.5)

      fig = matplotlib.pyplot.gcf()
      fig.set_size_inches(3.5,6)
      plt.axis("scaled")

      plt.savefig(saveDir*"/y_"*string(ipict)*".png")
      close()
   end

end

"""
   get_diamagnetic_velocity(filename, filedir, nType)

Calculate diamagnetic velocity from outputs.

nType = 1 is Hall MHD, nType = 2 is PIC.
"""
function get_diamagnetic_velocity(filename::String, filedir=".", nType=1)

   saveDir = "ycuts"
   mkdir(saveDir)

   if filename[end-3:end] == "outs"
      data = readdata(filename, dir=filedir, verbose=false)
      npict = data.list.npictinfiles
      for i = 1:npict
         data = readdata(filename, dir=filedir, npict=i)
         processing(data, saveDir, nType)
      end
   else
      filenames = glob(filename, filedir)
      npict = length(filenames)
      for i = 1:npict
         data = readdata(filenames[i], dir=filedir)
         processing(data, saveDir, nType)
      end
   end

end

function processing(data::Data, saveDir::String, nType=1)

   np = pyimport("numpy")

   @views x = data.x[:,:,:,1]
   @views y = data.x[:,:,:,2]
   @views z = data.x[:,:,:,3]
   w = data.w

   if nType == 1 # Hall MHD
      ni_ = findfirst(isequal("Rho"), data.head.wnames)
      pi_ = findfirst(isequal("P"),   data.head.wnames)
   elseif nType == 2 # PIC
      ni_ = findfirst(isequal("rhos1"), data.head.wnames)
      pi_ = findfirst(isequal("ps1"),   data.head.wnames)
   end

   bx_ = findfirst(isequal("Bx"),  data.head.wnames)
   by_ = findfirst(isequal("By"),  data.head.wnames)
   bz_ = findfirst(isequal("Bz"),  data.head.wnames)

   ni = w[:,:,:,ni_] # [/cm^3]
   pi = w[:,:,:,pi_] # [nPa]
   Bx = w[:,:,:,bx_] # [nT]
   By = w[:,:,:,by_] # [nT]
   Bz = w[:,:,:,bz_] # [nT]

   # This assumes uniform grid resolution!
   px, py, pz = np.gradient(pi, x[2]-x[1])

   nX, nY, nZ = size(w)

   vDiaMag = Array{Float64,4}(undef,3,nX,nY,nZ)
   vMag    = Array{Float64,3}(undef,nX,nY,nZ)

   cUnit = 1.0e7/(14*1.602*2.634)

   # v = B×▽p/neB^2
   @inbounds for k = 1:nZ, j = 1:nY, i = 1:nX
      Bmag2 = Bx[i,j,k]^2 + Bx[i,j,k]^2 + Bz[i,j,k]^2
      vDiaMag[1,i,j,k] = (By[i,j,k]*pz[i,j,k] - Bz[i,j,k]*py[i,j,k]) /
         (ni[i,j,k]*Bmag2) * cUnit
      vDiaMag[2,i,j,k] = (Bz[i,j,k]*px[i,j,k] - Bx[i,j,k]*pz[i,j,k]) /
         (ni[i,j,k]*Bmag2) * cUnit
      vDiaMag[3,i,j,k] = (Bx[i,j,k]*py[i,j,k] - By[i,j,k]*px[i,j,k]) /
         (ni[i,j,k]*Bmag2) * cUnit
      vMag[i,j,k] = hypot(vDiaMag[:,i,j,k]...)
   end

   yMid = floor(Int,nY/2)
   xrange = floor(Int,nX/4):floor(Int,nX/4*3)
   zrange = floor(Int,nZ/4):floor(Int,nZ/4*3)

   clf()

   fig, ax = subplots()
   cont_levels = range(0.0, 10000.0, length=50)
   c = ax.contourf(x[xrange,yMid,zrange], z[xrange,yMid,zrange],
      vMag[xrange,yMid,zrange], levels=cont_levels)
   fig.colorbar(c, ax=ax, ticks=range(0.0, 10000.0, step=1000.0))
   xlabel("x"); ylabel("z")
   title(L"$\mathbf{B}\times \nabla P/(neB^2)$"*", t=$(i)")

   fig.set_size_inches(3.5,6)
   plt.axis("scaled")

   plt.savefig(saveDir*"/diamagnetic_$(i).png")
end

#=
filename = "box*outs";
dir = "/Users/hyzhou/Ganymede/run_mercury_80s_newbox/GM";

filename = "y*outs";
dir = "/Users/hyzhou/Ganymede/run_mercury_80s_newbox/GM";
=#

#filename = "cut*.dat";
#dir = "GM/IO2";
#vtk_auto_conversion(filename, dir)

function y_by_cuts()
   np = pyimport("numpy")
   filename = "y*out"
   data = readdata(filename)

   fig, ax = plt.subplots(1, 1)
   var = "by"; plotinterval = 0.1
   VarIndex_ = findfirst(x->x==lowercase(var), lowercase.(data.head.wnames))
   X = vec(data.x[:,:,1])
   Y = vec(data.x[:,:,2])
   W = vec(data.w[:,:,VarIndex_])

   # Create grid values first.
   xi = range(minimum(X), stop=maximum(X), step=plotinterval)
   yi = range(minimum(Y), stop=maximum(Y), step=plotinterval)
   # Perform linear interpolation of the data (x,y) on grid(xi,yi)
   triang = matplotlib.tri.Triangulation(X,Y)
   interpolator = matplotlib.tri.LinearTriInterpolator(triang, W)
   Xi, Yi = np.meshgrid(xi, yi)
   wi = interpolator(Xi, Yi)

   c = contourf(xi,yi,wi,50)
   fig.set_size_inches(6,6)
   plt.axis("scaled")
   colorbar(c, ax=ax,ticks=np.arange(-100.0, 0.0, 12.5))
   xlabel("x"); ylabel("z")
   title(L"By")
   plt.savefig("test.png")
end

function plot_beta(filename::String)

   data = readdata(filename)

   cutPlaneIndex = 65
   VarIndex_ = 18

   X = @view data.x[:,:,:,1]
   Y = @view data.x[:,:,:,2]
   Z = @view data.x[:,:,:,3]

   fig, ax = plt.subplots(1, 1)
   fig.set_size_inches(3.3,6);

   W = @view data.w[:,:,:,VarIndex_]
   Bx = @view data.w[:,:,:,5]
   By = @view data.w[:,:,:,6]
   Bz = @view data.w[:,:,:,7]

   cut1 = @view X[:,cutPlaneIndex,:]
   cut2 = @view Z[:,cutPlaneIndex,:]
   W    = @view W[:,cutPlaneIndex,:]
   Bx   = @view Bx[:,cutPlaneIndex,:]
   By   = @view By[:,cutPlaneIndex,:]
   Bz   = @view Bz[:,cutPlaneIndex,:]
   PB   = hypot.(Bx, By, Bz)

   c = ax.contourf(cut1, cut2, W./PB)
   fig.colorbar(c,ax=ax)
   ax.axis("scaled")
   title(data.head.wnames[VarIndex_])

   xlabel("x"); ylabel("z")

   dim = [0.125, 0.013, 0.2, 0.045]
   str = @sprintf "it=%d, time=%4.2f" filehead[:it] filehead[:time]
   at = matplotlib.offsetbox.AnchoredText(str,
      loc="lower left", prop=Dict("size"=>8), frameon=true,
      bbox_to_anchor=(0., 1.),
      bbox_transform=ax.transAxes)
      at.patch.set_boxstyle("round,pad=0.,rounding_size=0.2")
      ax.add_artist(at)
   #plt.savefig("beta.png")
end


text:
"""
    NondeterministicThreading(; basesize, ntasks = nthreads())

Parallelize inner reducing function using `ntasks`.

Given

```julia
#                   ,-- Not parallelized
#            ______/__
  foldxl(rf, xs |> xfo |> NondeterministicThreading() |> xfi)
#        ==                                              ===
#        Parallelized                                    Parallelized
```

the inner transducer `xfi` and the reducing function `rf` are
parallelized but the outer transducer `xfo` and the iteration over
data collection `xs` are not parallelized.

The scheduling of the tasks (hence the call tree of the inner reducing
function) is non-deterministic.  It means that the result is
deterministic only if the inner reducing function is exactly
associative.  If the inner reducing function is only approximately
associative (e.g., addition of floating point numbers), the result of
reduction is likely different for each time.

The outer transducers and `iterate` are processed sequentially.  In
particular, the data collection does not have to implement
`SplittablesBase.halve`.

!!! note
    Currently, the default `basesize` is 1.  However, it may be
    changed in the future (e.g. it may be automatically tuned at
    run-time).

!!! note
    `NondeterministicThreading` does not work with Julia < 1.3.

# Keyword Arguments
- `basesize::Integer`: The number of input elements to be accumulated
  in a buffer before sent to a task.
- `ntasks::Integer`: The number of tasks `@spawn`ed.  The default
  value is `Threads.nthreads()`.  A number larger than
  `Threads.nthreads()` may be useful if the inner reducing function
  contains I/O and does not consume too much resource (e.g., memory).

# Examples
```jldoctest
julia> using Transducers

julia> collect(1:4 |> Filter(isodd))
2-element Vector{Int64}:
 1
 3

julia> collect(1:4 |> NondeterministicThreading() |> Filter(isodd))
2-element Vector{Int64}:
 1
 3
```
"""
struct NondeterministicThreading <: Transducer
    basesize::Int
    ntasks::Int
end

NondeterministicThreading(; basesize = 1, ntasks = Threads.nthreads()) =
    NondeterministicThreading(basesize, ntasks)

function start(rf::R_{NondeterministicThreading}, init)
    buffer = UndefVector(xform(rf).basesize)
    n = xform(rf).basesize - 1
    ntasks = xform(rf).ntasks
    irf = inner(rf)
    taskref = Ref{Task}()
    err_promise = Promise()
    ichan = Channel(0; taskref = taskref) do ichan
        async_foreach(1:ntasks) do _
            sync_spawn_foreach(ichan) do (lbridge, buffer, rbridge)
                try
                    acc = foldl_nocomplete(irf, start(irf, init), buffer)
                    while true
                        x = tryfetch(lbridge::Union{Nothing,Promise})
                        x isa Some{<:Ok} || break
                        lacc, lbridge = something(x).value
                        acc = combine_reduced(irf, lacc, acc)
                    end
                    while true
                        x = tryfetch(rbridge::Promise)
                        x isa Some{<:Ok} || break
                        racc, rbridge = something(x).value
                        acc = combine_reduced(irf, acc, racc)
                    end
                    while true
                        if lbridge isa Promise && isodd(time_ns())
                            x = tryput!(lbridge::Promise, Ok((acc, rbridge)))
                            x isa Some{<:Ok} || break
                            lacc, lbridge = something(x).value
                            acc = combine_reduced(irf, lacc, acc)
                        else
                            x = tryput!(rbridge::Promise, Ok((acc, lbridge)))
                            x isa Some{<:Ok} || break
                            racc, rbridge = something(x).value
                            acc = combine_reduced(irf, acc, racc)
                        end
                    end
                catch err
                    @debug(
                        "start(::R_{NondeterministicThreading}, _)",
                        exception = (err, catch_backtrace())
                    )
                    tryput!(err_promise, err)
                    close(ichan)
                    rethrow()
                end
            end
        end
    end
    return (n, buffer, nothing, ichan, taskref[], init, err_promise)
end

function next(
    rf::R_{NondeterministicThreading},
    (n, buffer, lbridge, ichan, task, init, err_promise),
    x,
)
    buffer = setindex!!(buffer, x, length(buffer) - n)
    if n == 0
        rbridge = Promise()
        try
            put!(ichan, (lbridge, buffer, rbridge))
        catch err
            if err isa InvalidStateException
                istaskdone(task) && fetch(task)  # let it throw
                let x = tryfetch(err_promise)
                    x === nothing || throw(something(x))
                end
            end
            rethrow()
        end
        lbridge = rbridge
        buffer = typeof(buffer)(undef, xform(rf).basesize)
        n = xform(rf).basesize - 1
    else
        n -= 1
    end
    return (n, buffer, lbridge, ichan, task, init, err_promise)
end

function complete(
    rf::R_{NondeterministicThreading},
    (n, buffer, lbridge, ichan, task, init, err_promise),
)
    if !(buffer isa UndefVector)
        ilast = length(buffer) - n - 1
        if ilast > 0
            resize!(buffer, ilast)
            rbridge = Promise()
            put!(ichan, (lbridge, buffer, rbridge))
            lbridge = rbridge
        end
    end
    close(ichan)
    if lbridge === nothing
        acc = start(inner(rf), init)
    else
        # Two tasks below share `lbridge`
        local errtask
        @sync begin
            errtask = @async begin
                local err = fetch(err_promise)
                err === nothing && return
                # Chase `lbridge` and propagate the error in the other
                # task below:
                while lbridge isa Promise
                    local x = tryput!(lbridge, Err(err))
                    x isa Some{<:Ok} || break
                    _, lbridge = something(x).value
                end
                Err(err)
            end
            acc = fetch(@async begin
                local acc
                local x = fetch(lbridge)
                x isa Ok || return
                acc, lbridge = x.value
                while lbridge isa Promise
                    local lacc
                    x = fetch(lbridge)
                    x isa Ok || return
                    lacc, lbridge = x.value
                    acc = combine_reduced(inner(rf), lacc, acc)
                end
                tryput!(err_promise, nothing)  # cancel `errtask`
                return acc
            end)
        end
        let x = fetch(errtask)
            x isa Err && throw(x.value)
        end
    end
    wait(task)
    return complete(inner(rf), acc)
end


text:
struct TimestampArray{D, T} <: AbstractVector{Timestamp{D, T}}
    tsps::Vector{Timestamp{D, T}}
    names::Vector{Symbol}
    TimestampArray(v::AbstractVector{Timestamp{D, T}}, names) where {D, T} = new{D, T}(sort!(v), names)
end

size(tsparr::TimestampArray, i) = size(tsparr)[i]
size(tsparr::TimestampArray) = (length(tsparr.tsps), length(tsparr.names))

@Base.propagate_inbounds getindex(tsparr::TimestampArray, i::Integer) = getindex(tsparr.tsps, i)
@Base.propagate_inbounds getindex(tsparr::TimestampArray, r::AbstractUnitRange) = TimestampArray(getindex(tsparr.tsps, r), tsparr.names)

Base.isempty(v::TimestampArray) = isempty(v.tsps)

function Base.push!(v::TimestampArray{D, T}, tsp::Timestamp{D, T}) where {D, T}
    if isempty(v)
        push!(v.tsps, tsp)
    elseif v[end] <= tsp
        push!(v.tsps, tsp)
    elseif v[1] >= tsp
        pushfirst!(v.tsps, tsp)
    else
        push!(v.tsps, tsp)
        sort!(v.tsps)
    end

    return v
end


text:
# A palindromic number reads the same both ways. The largest palindrome made
# from the product of two 2-digit numbers is 9009 = 91 × 99.
#
# Find the largest palindrome made from the product of two 3-digit numbers.

using ProjectEulerSolutions

# Treat numbers as strings and test by reversing them.
function p004solution_string(n_digits::Integer=1)::Integer
    testnums = sort(unique_products_of_ndigits(n_digits), rev=true)
    for i in testnums
        if ispalindrome_string(i)
            return i
        end
    end
    return 0
end

# Treat numbers as integers and test by reversing them using modulus and floor
# to pull out individual digits.
function p004solution_integer(n_digits::Integer=1)::Integer
    testnums = sort(unique_products_of_ndigits(n_digits), rev=true)
    for i in testnums
        if ispalindrome(i)
            return i
        end
    end
    return 0
end

# Same as previous, but no initial list construction using list comprehension.
function p004solution_integer_fast(n_digits::Integer=1)::Integer
    maxnum = 10^n_digits - 1
    minnum = 10^(n_digits - 1) # Can probably assume this can be larger
    maxval = 0
    for i in minnum:maxnum
        for j in minnum:maxnum
            if j > i && ispalindrome_integer(i*j)
                maxval = max(maxval, i*j)
            end
        end
    end
    return maxval
end

p004 = Problems.Problem(Dict("String reversal" => p004solution_string,
                             "Integer math" => p004solution_integer,
                             "Integer math fast" => p004solution_integer_fast))

Problems.benchmark(p004, 3)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# This file contains methods to solve an instance (heuristically or with CPLEX)
using CPLEX
using Printf
include("generation.jl")

TOL = 0.00001
directions=[1 0;-1 0;0 1;0 -1]

"""
Solve an instance with CPL

Argument
-t: array which stores the puzzle data

Return

"""

function cplexSolve(t::Array{Int,2})

    rowNum=size(t,1)
    colNum=size(t,2)
    if rowNum!=colNum
        println("not square matrix")
        return
    end
    @printf("%d-%d\n",rowNum,colNum)
    # Create the model
    m = Model(CPLEX.Optimizer)

    # x[i,j]=1 if the cell is white, 0 if the cell is black
    @variable(m,x[1:rowNum,1:colNum],Bin)

    ss=0
    # set the cell white if it has a number
    for l in 1:rowNum
        for c in 1:colNum
            if t[l,c]!=0
                @constraint(m,x[l,c]==1)
                ss=ss+1
            end
        end
    end

    # No two black squares are orthogonally adjacent
    @constraint(m,[r in 1:rowNum-1, c in 1:colNum],x[r,c]+x[r+1,c]>=1)
    @constraint(m,[r in 1:rowNum, c in 1:colNum-1],x[r,c]+x[r,c+1]>=1)

    # For any two white squares, there is a path between them using only white squares.
    @constraint(m,[c in 2:colNum-1],sum(x[1,c-1:c+1])+x[2,c]>=2) # the first line except the corners
    @constraint(m,[c in 2:colNum-1],sum(x[rowNum,c-1:c+1])+x[rowNum-1,c]>=2) # the last line except the corners
    @constraint(m,[r in 2:rowNum-1],sum(x[r-1:r+1,1])+x[r,2]>=2) # the first column except the corners
    @constraint(m,[r in 2:rowNum-1],sum(x[r-1:r+1,colNum])+x[r,colNum-1]>=2) # the last colNum except the corners
    @constraint(m,sum(x[1:2,1])+x[1,2]>=2) # the left top corner
    @constraint(m,sum(x[1:2,colNum])+x[1,colNum-1]>=2) # the right top corner
    @constraint(m,sum(x[rowNum-1:rowNum,1])+x[rowNum,2]>=2) # the left bottom corner
    @constraint(m,sum(x[rowNum-1:rowNum,colNum])+x[rowNum,colNum-1]>=2) # the right bottom corner

    @constraint(m,[r in 2:rowNum-1, c in 2:colNum-1],sum(x[r-1:r+1,c-1:c+1])>=2)

    #case when the black is in diagonal line
    @constraint(m,sum(x[r,r] for r=1:rowNum)>=1)
    @constraint(m, sum(x[rowNum-r+1,r] for r=1:rowNum)>=1)

    function my_callback_function(cb_data)
        x_vals = callback_value.(Ref(cb_data), x)
        zz=0
        para=Array{Int64,2}(undef,0,6)
        for l in 1:rowNum
            for c in 1:colNum
                if t[l,c]!=0
                    up=0 # the number of white cells in top direction until meet black cell
                    right=0 # the number of white cells in right direction until meet black cell
                    down=0 # the number of white cells in down direction until meet black cell
                    left=0 # the number of white cells in left direction until meet black cell
                    for i1 in 1:l-1
                        if  x_vals[l-i1,c]== 0  #meet the black cell
                            down=i1-1
                            break
                        else
                            down=l-1
                        end
                    end
                    for j1 in 1:rowNum-l
                        if  x_vals[l+j1,c]== 0
                            up=j1-1
                            break
                        else
                            up=rowNum-l
                        end
                    end
                    for i2 in 1:c-1
                        if  x_vals[l,c-i2]== 0
                            left=i2-1
                            break
                        else
                            left=c-1
                        end
                    end
                    for j2 in 1:colNum-c
                        if  x_vals[l,c+j2]== 0
                            right=j2-1
                            break
                        else
                            right=colNum-c
                        end
                    end
                    if up+down+left+right+1>t[l,c]
                        # println(up+down+left+right+1)
                        con1 = @build_constraint((sum(x[l:l+up,c])+sum(x[l-down:l,c])+sum(x[l,c-left:c])+sum(x[l,c:c+right])-3)>=t[l,c]) #counting the white cells from four orthogonal directions
                        # # # con1=@build_constraint(up+down+left+right+1==t[l,c])
                        MOI.submit(m, MOI.LazyConstraint(cb_data), con1)
                    elseif up+down+left+right+1<t[l,c]
                        con = @build_constraint((sum(x[l:l+up,c])+sum(x[l-down:l,c])+sum(x[l,c-left:c])+sum(x[l,c:c+right])-3)<=t[l,c])
                        MOI.submit(m, MOI.LazyConstraint(cb_data), con)
                    elseif up+down+left+right+1==t[l,c]
                        zz=zz+1
                        para=vcat(para,[l c up down left right])
                    end
                end
            end
        end
        # println(zz,"/",ss)
        if(zz>=ss)

            for i=1:ss
                l,c,up,down,left,right=para[i,:]
                # @printf("func:%d+%d+%d+%d+1=%d\n",up,down,left,right,t[l,c]) #up,down,left,right,target
                # con = @build_constraint((sum(x[l:l+up,c])+sum(x[l-down:l,c])+sum(x[l,c-left:c])+sum(x[l,c:c+right])-3)==t[l,c])

                # MOI.submit(m, MOI.LazyConstraint(cb_data), con)

            end
        end
    end

    MOI.set(m, MOI.LazyConstraintCallback(), my_callback_function)
    # Maximize the number of white cells
    @objective(m,Max, sum(x))

    # Start a chronometer
    start = time()

    # Solve the model
    optimize!(m)

    # println(num_constraints(m, VariableRef, MOI.ZeroOne))
    # println("-----",termination_status(m))
    # Return:
    # 1 - true if an optimum is found(type: Bool)
    # 2 - the resolution time(type Float64)
    # 3 - the value of each cell(type:Array{undef,2})
    return JuMP.primal_status(m) == JuMP.MathOptInterface.FEASIBLE_POINT, time() - start,x

end

"""
Heuristically solve an instance
"""
function visiable(x_vals::Array{Bool, 2},l::Int,c::Int)
    rowNum=size(x_vals,1)
    colNum=size(x_vals,2)

    up=0 # the number of white cells in top direction until meet black cell
    right=0 # the number of white cells in right direction until meet black cell
    down=0 # the number of white cells in down direction until meet black cell
    left=0 # the number of white cells in left direction until meet black cell
    for i1 in 1:l-1
        if  x_vals[l-i1,c]== 0  #meet the black cell
            down=i1-1
            break
        else
            down=l-1
        end
    end
    for j1 in 1:rowNum-l
        if  x_vals[l+j1,c]== 0
            up=j1-1
            break
        else
            up=rowNum-l
        end
    end
    for i2 in 1:c-1
        if  x_vals[l,c-i2]== 0
            left=i2-1
            break
        else
            left=c-1
        end
    end
    for j2 in 1:colNum-c
        if  x_vals[l,c+j2]== 0
            right=j2-1
            break
        else
            right=colNum-c
        end
    end
    return up+down+left+right+x_vals[l,c]
end

function dfs(l::Int,c::Int,x_vals::Array{Bool, 2},walked::Array{Bool, 2})

    rowNum=size(x_vals,1)
    colNum=size(x_vals,2)
    ans=1
    for i in 1:4
        r1=l+directions[i,1]
        c1=c+directions[i,2]
        if r1>=1 && r1<=rowNum && c1>=1 && c1<=colNum && walked[r1,c1]==0 && x_vals[r1,c1]!=0
            walked[r1,c1]=1
            ans+=dfs(r1,c1,x_vals,walked)
        end
    end
    walked=zeros(Bool,rowNum,colNum)
    return ans
end

"""
check if it is possible to place a black cell in position(r,c) and if possible, place black cell
return true if possible, false if impossible
"""
function placeBlack(r::Int, c::Int,x_vals::Array{Bool,2})
    rowNum=size(x_vals,1)
    colNum=size(x_vals,2)
    # if the cell is already a black one
    if x_vals[r,c]!=1
        return false
    end

    # there is black cells in its 4 orthogonal neighbors
    for i=1:size(directions,1)
        if r+directions[i,1]>=1 && r+directions[i,1]<=rowNum && c+directions[i,2]>=1 &&c+directions[i,2]<=colNum && x_vals[r+directions[i,1],c+directions[i,2]]==0
            return false
        end
    end
    x_vals[r,c]=0
    return true
end

"""
unplace the black cell, change it to white
"""
function unplaceBlack(r::Int, c::Int,x_vals::Array{Bool,2})
    x_vals[r,c]=1
end

"""
iterativly solve the puzzle
"""
function solve(r::Int,c::Int,rw::Int,cw::Int,x_vals::Array{Bool,2},walked::Array{Bool,2},t::Array{Int,2})
    rowNum=size(x_vals,1)
    colNum=size(x_vals,2)

    if c > colNum
        c = 1
        r+=1
    end
    # check if white cells are continuous
    reach=0
    walked[rw,cw]=1
    if x_vals[rw,cw]!=1
        x_vals[rw,cw]=1
        reach+=dfs(rw,cw,x_vals,walked)
        x_vals=0
    else
        reach+=dfs(rw,cw,x_vals,walked)
    end
    walked=zeros(Bool,rowNum,colNum)
    if reach!=sum(x_vals)
        # println("white cells are not contiguous,reach=",reach,"sum=",sum(x_vals))
        return 0
    end
    # check after all iteration whether the puzzle is finally solve
    if r>rowNum
        solved=true
        for i in 1:rowNum
            for j in 1:colNum
                if t[i,j]>0 # for the cell who has a number
                    if visiable(x_vals,i,j)!=t[i,j] # not solved
                        # println("no solution")
                        solved=false
                        return 0
                    end
                end
            end
        end
        if solved==true
            # finally solved
            # println("solve successfully")
            # println(x_vals)
            return x_vals
        end
        return 0
    end
    # check in each iteration
    for i in 1:rowNum
        for j in 1:colNum
            if t[i,j]>0 # for the cell who has a number
                if visiable(x_vals,i,j)<t[i,j] # not solved
                    # println(visiable(x_vals,i,j)," less than ",t[i,j])
                    return 0
                end
                if i<r && visiable(x_vals[1:r-1,:],i,j)>t[i,j]
                    # println("before ",r)
                    return 0
                end
                # println("nothing:",visiable(x_vals,i,j),"->",t[i,j])
            end
        end
    end
    if placeBlack(r,c,x_vals)==true
        rw1=rw
        cw1=cw
        # update the position of white cell
        if rw1==r && cw1==c
            cw1+=1 # move right
            if cw1>colNum
                cw1=1
                rw1+=1
            end
        end
        # println("new->",r,c+1,rw1,cw1,x_vals,walked,t)
        ret_val=solve(r,c+1,rw1,cw1,x_vals,walked,t)
        # println(ret_val)
        if ret_val!=0 && ret_val!=nothing
            return ret_val
        end
        unplaceBlack(r,c,x_vals)
    end
    walked=zeros(Bool,rowNum,colNum)
    # println("new start",r,c+1,rw,cw,x_vals,walked,t)
    ret_val=solve(r,c+1,rw,cw,x_vals,walked,t)
    if ret_val!=0 && ret_val!=nothing
        return ret_val
    end
    # println(ret_val)
end

function heuristicSolve(t::Array{Int, 2})
    rowNum=size(t,1)
    colNum=size(t,2)
    walked=zeros(Bool,rowNum,colNum)
    x_vals=ones(Bool,rowNum,colNum)
    ret_val=solve(1,1,1,1,x_vals,walked,t)
    isOptimal=false
    if ret_val!=0 && ret_val!=nothing
        isOptimal=true
        println("solve successfully")
        println(ret_val)
    else
        isOptimal=false
        println("no solution")
        println(ret_val)
    end
    return isOptimal,ret_val
end

"""
Solve all the instances contained in "../data" through CPLEX and heuristics

The results are written in "../res/cplex" and "../res/heuristic"

Remark: If an instance has previously been solved (either by cplex or the heuristic) it will not be solved again
"""
function solveDataSet()

    dataFolder = "../data/"
    resFolder = "../res/"

    # Array which contains the name of the resolution methods
    # resolutionMethod = ["cplex"]
    resolutionMethod = ["cplex", "heuristique"]

    # Array which contains the result folder of each resolution method
    resolutionFolder = resFolder .* resolutionMethod

    # Create each result folder if it does not exist
    for folder in resolutionFolder
        if !isdir(folder)
            mkdir(folder)
        end
    end

    global isOptimal = false
    global solveTime = -1

    # For each instance
    # (for each file in folder dataFolder which ends by ".txt")
    for file in filter(x->occursin(".txt", x), readdir(dataFolder))

        println("-- Resolution of ", file)
        t=readInputFile(dataFolder * file)

        # For each resolution method
        for methodId = 1:size(resolutionMethod,1)

            outputFile = resolutionFolder[methodId] * "/" * file

            # If the instance has not already been solved by this method
            if !isfile(outputFile)

                fout = open(outputFile, "w")

                resolutionTime = -1
                isOptimal = false

                # If the method is cplex
                if resolutionMethod[methodId] == "cplex"

                    # Solve it and get the results
                    isOptimal, resolutionTime,x = cplexSolve(t)

                    # If a solution is found, write it
                    if isOptimal
                        writeSolution(fout, x)
                    end

                # If the method is one of the heuristics
                else

                    isSolved = false

                    # Start a chronometer
                    startingTime = time()

                    solution=nothing

                    # While the grid is not solved and less than 100 seconds are elapsed

                    # Solve it and get the results
                    isOptimal, solution = heuristicSolve(t)

                    # Stop the chronometer
                    resolutionTime = time() - startingTime

                    # Write the solution (if any)
                    if isOptimal

                        writeSolution(fout, solution)

                    end
                end

                println(fout, "solveTime = ", resolutionTime)
                println(fout, "isOptimal = ", isOptimal)

                close(fout)
            end

            # Display the results obtained with the method on the current instance
            include(outputFile)
            println(resolutionMethod[methodId], " optimal: ", isOptimal)
            println(resolutionMethod[methodId], " time: " * string(round(solveTime, sigdigits=2)) * "s\n")
        end
    end
end


text:
function _rand_continuous(p::D1, q::D2) where {D1 <: Distribution, D2 <: Distribution}
    X = rand(p)
    log_pX = logpdf(p, X)
    log_qX = logpdf(q, X)
    if isnan(log_qX)
        # logpdf can return NaN if X is outside of
        # the support of q, in which case we should
        # change it to -Inf
        log_qX = -Inf
    end
    logU = -randexp()
    local Y
    if logU <= ((log_qX - log_pX) + eps()) # eps gives a bit of margin
        Y = X
    else
        found_uncoupled = false
        for _ in 1:1_000_000 # give up after a while
            Ystar = rand(q)
            log_pYstar = logpdf(p, Ystar)
            if isnan(log_pYstar)
                # logpdf can return NaN if Ystar is outside of
                # the support of p, in which case we should
                # change it to -Inf
                log_pYstar = -Inf
            end
            log_qYstar = logpdf(q, Ystar)
            if log_qYstar < log_pYstar
                continue
            end
            logUstar = -randexp()
            if logUstar >= ((log_pYstar - log_qYstar) - eps())
                Y = Ystar
                found_uncoupled = true
                break
            end
        end
        if !found_uncoupled
            @warn("failed to find uncoupled sample between ", p, "\n and \n", q)
            # just return a coupled sample after all
            Y = X
        end
    end
    return X, Y
end
function _rand_discrete(p::D1, q::D2) where {D1 <: DiscreteDistribution, D2 <: DiscreteDistribution}
    q_support = support(q)
    # overlap = p_support ∩ q_support
    # prob_min = [min(pdf(p, x), pdf(q, x)) for x in overlap]
    # prob_couple = sum(prob_min)

    X = rand(p)
    log_pX = logpdf(p, X)
    log_qX = logpdf(q, X)
    if isnan(log_qX)
        # logpdf can return NaN if X is outside of
        # the support of q, in which case we should
        # change it to -Inf
        log_qX = -Inf
    end
    logU = -randexp()
    local Y
    if logU <= ((log_qX - log_pX) + eps()) # eps gives a bit of margin
        Y = X
    else
        q_minus_p = [max(pdf(q, x)-pdf(p,x), 0.0) for x in q_support]
        Y = sample(q_support, Weights(q_minus_p))

    end

    return X, Y
end

get_Σ_PD(d::MvNormal) = d.Σ
standard(d::MvNormal) = MvNormal(length(d), 1.0)
samecov(d1::MvNormal, d2::MvNormal) = d1.Σ.chol == d2.Σ.chol

"""
    Obtain coupled random numbers using the algorithm of Bou-Rabee et al. [2018].
    See also Jacob et al. 2018, section 4.1 on Sampling from maximal couplings.
    I follow the procedure and notation used in the Jacob paper below.
"""
function _rand_BouRabee(p::D, q::D) where {D <: ContinuousMultivariateDistribution}
    @assert samecov(p, q)
    Σp, Σq = get_Σ_PD(p), get_Σ_PD(q)
    Σ = Σp
    s_distr = standard(p)
    μp, μq = mean(p), mean(q)
    z = whiten(Σ, μp-μq)
    e = z/norm(z)

    Ẋ = rand(s_distr)
    X = μp + unwhiten(Σ, Ẋ)
    logU = -randexp()

    Ẏ = Ẋ + z
    @assert μq+unwhiten(Σ, Ẏ) ≈ X
    logsẊ = logpdf(s_distr, Ẋ)
    logsẎ = logpdf(s_distr, Ẏ)
    local Y
    if logU <= ((logsẎ - logsẊ) + eps())
        Y = X
    else
        Ẏ = Ẋ - 2*dot(e, Ẋ)*e
        Y = μq + unwhiten(Σ, Ẏ)
    end
    return X, Y
end

function rand(coup::MaximalCoupling{D1,D2}) where {D1, D2}
    _rand_continuous(coup.p, coup.q)
end
function rand(coup::MaximalCoupling{D1,D2}) where {D1<:DiscreteDistribution, D2<:DiscreteDistribution}
    p_finite_support = isfinite(minimum(coup.p)) && isfinite(maximum(coup.p))
    q_finite_support = isfinite(minimum(coup.q)) && isfinite(maximum(coup.q))
    if q_finite_support
        return _rand_discrete(coup.p, coup.q)
    elseif p_finite_support
        nx = length(coup.p)
        ny = length(coup.q)
        Y, X = _rand_discrete(coup.q, coup.p)
        return X, Y
    else
        return _rand_continuous(coup.p, coup.q)
    end
end
function rand(coup::MaximalCoupling{D1,D2}) where {D1<:MvNormal, D2<:MvNormal}
    if samecov(coup.p, coup.q)
        return _rand_BouRabee(coup.p, coup.q)
    else
        return _rand_continuous(coup.p, coup.q)
    end
end

function rand(coup::MaximalCoupling{D1, D2}, n::Int64) where {D1, D2}
    X = Vector{eltype(coup)}(undef, n)
    Y = Vector{eltype(coup)}(undef, n)
    for i in 1:n
        xi, yi = rand(coup)
        X[i] = xi
        Y[i] = yi
    end
    return X, Y
end
function rand(coup::MaximalCoupling{D1, D2}, n::Int64) where {D1 <: MultivariateDistribution, D2 <: MultivariateDistribution}
    X = Matrix{eltype(coup)}(undef, length(coup), n)
    Y = Matrix{eltype(coup)}(undef, length(coup), n)
    for i in 1:n
        xi, yi = rand(coup)
        X[:, i] = xi
        Y[:, i] = yi
    end
    return X, Y
end


text:
function isless(a::Float120, b::Float120)
    (a.hi < b.hi) || (a.hi==b.hi && a.lo<b.lo)
end
function isless(a::Float120, b::Float64)
    (a.hi < b) || (a.hi==b && a.lo<zero(Float64))
end
function isless(a::Float64, b::Float120)
    (a < b.hi) || (a==b.hi && b.lo<zero(Float64))
end
function isequal(a::Float120, b::Float120)
    (a.hi == b.hi) && (a.lo == b.lo)
end
function isequal(a::Float120, b::Float64)
    (a.hi == b) && (a.lo == zero(Float64))
end
function isequal(a::Float64, b::Float120)
    (a == b.hi) && (b.lo == zero(Float64))
end

@inline (==)(a::Float120,b::Float120) = (a.hi == b.hi) && (a.lo == b.lo)
@inline (< )(a::Float120,b::Float120) = (a.hi < b.hi) || (a.hi==b.hi && a.lo<b.lo)
@inline (<=)(a::Float120,b::Float120) = (a.hi < b.hi) || (a.hi==b.hi && a.lo<=b.lo)
@inline (> )(a::Float120,b::Float120) = (a.hi > b.hi) || (a.hi==b.hi && a.lo>b.lo)
@inline (>=)(a::Float120,b::Float120) = (a.hi > b.hi) || (a.hi==b.hi && a.lo>=b.lo)

min{T<:Float120}(a::T, b::T) = (a<b) ? a : b
max{T<:Float120}(a::T, b::T) = (b<a) ? a : b
minmax{T<:Float120}(a::T, b::T) = (a<b) ? (a,b) : (b,a)
maxmin{T<:Float120}(a::T, b::T) = (b<a) ? (a,b) : (b,a)


text:
push!(LOAD_PATH,"./src/")

using sail_route

using Base.Test
@time include("aerodynamics.jl")
@time include("route.jl")
@time include("comparison.jl")
@time include("route.jl")
@time include("performance.jl")


text:
using  NeuronBuilder, ModelingToolkit, OrdinaryDiffEq, Plots

## xolotl versions divided by 10
AB1_channels = [NaV(100.), CaS(6.), CaT(2.5), H(0.01), Ka(50.), KCa(5.), Kdr(100.), Leak(0.)]

LP1_channels = [NaV(100.), CaS(4.), CaT(0.), H(0.05), Ka(20.), KCa(0.), Kdr(25.), Leak(0.03)]

PY1_channels = [NaV(100.), CaS(2.), CaT(2.4), H(0.05), Ka(50.), KCa(0.), Kdr(125.), Leak(0.01)]

## parameters from O'Leary
AB1_channels = [NaV(1831.), CaS(27.), CaT(23.), H(10.1), Ka(246.), KCa(980.), Kdr(610.), Leak(0.99)]

PY1 = build_neuron(PY1_channels; name = :PY1)
LP1 = build_neuron(PY1_channels; name = :LP1)
AB1 = build_neuron(AB1_channels; name = :AB1)

final_PY1 = structural_simplify(PY1)
probPY1 = ODEProblem(final_PY1, [], (0., 3000.), [])

final_LP1 = structural_simplify(LP1)
probLP1 = ODEProblem(final_LP1,[],
                (0.,3000.), [])

final_AB1 = structural_simplify(AB1)

probAB1 = ODEProblem(final_AB1, [], (0., 3000.), [])

sol = solve(probAB1, ImplicitEuler(), dt = 0.025)

plot(sol)

max_stars_repo_name: sg-s/NeuronBuilder

max_issues_repo_name: sg-s/NeuronBuilder

max_forks_repo_name: sg-s/NeuronBuilder


text:
module TwoVarianceComponentTest

using VarianceComponentModels, MathProgBase, Ipopt, Base.Test

srand(123)

# generate data from a d-variate response variane component model
n = 100   # no. observations
d = 2     # no. categories
m = 2     # no. variance components
p = 2     # no. covariates
X = randn(n, p)
B = ones(p, d)
Σ = ntuple(x -> zeros(d, d), m)
for i in 1:m
  Σi = randn(d, d)
  copy!(Σ[i], Σi' * Σi)
end
## make the first variance component 0 matrix
#fill!(Σ[1], 0.0)
V = ntuple(x -> zeros(n, n), m)
for i = 1:m-1
  Vi = randn(n, 50)
  copy!(V[i], Vi * Vi')
end
copy!(V[m], eye(n))
# form Ω
Ω = zeros(n*d, n*d)
for i = 1:m
  Ω += kron(Σ[i], V[i])
end
Ωchol = cholfact(Ω)
Y = X * B + reshape(Ωchol[:L] * randn(n*d), n, d)

info("Forming VarianceComponentModel from data")
@inferred VarianceComponentVariate(Y, X, V)
vcdata  = VarianceComponentVariate(Y, X, V)
vcmodel = VarianceComponentModel(vcdata)

info("Pre-compute eigen-decomposition and rotate data")
vcdatarot  = TwoVarCompVariateRotate(vcdata)
vcmodelrot = TwoVarCompModelRotate(vcmodel)

info("Evaluate log-pdf")
#@code_warntype logpdf(vcmodelrot, vcdatarot)
@inferred logpdf(vcmodelrot, vcdatarot)
@test logpdf(vcmodel, vcdata) == logpdf(vcmodelrot, vcdatarot)
@test (logpdf(vcmodelrot, [vcdatarot vcdatarot; vcdatarot vcdatarot]) -
  logpdf(vcmodel, [vcdata vcdata; vcdata vcdata])) ≈ 0.0

info("Evaluate gradient")
∇ = zeros(2d^2)
#@code_warntype gradient!(∇, vcmodelrot, vcdatarot)
@inferred gradient!(∇, vcmodelrot, vcdatarot)
@test vecnorm(gradient(vcmodel, vcdata) - gradient(vcmodelrot, vcdatarot)) ≈ 0.0
@test vecnorm(gradient(vcmodel, vcdata) - gradient(vcmodel, vcdatarot)) ≈ 0.0
@test vecnorm(gradient(vcmodel, [vcdata vcdata]) -
  2.0gradient(vcmodel, vcdata)) ≈ 0.0
@test vecnorm(gradient(vcmodel, [vcdata vcdata]) -
  gradient(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

info("Evaluate Fisher information matrix of Σ")
H = zeros(2d^2, 2d^2)
#@code_warntype fisher!(H, vcmodelrot, vcdatarot)
@inferred fisher_Σ!(H, vcmodelrot, vcdatarot)
@test vecnorm(fisher_Σ(vcmodel, vcdata) - fisher_Σ(vcmodelrot, vcdatarot)) ≈ 0.0
@test vecnorm(fisher_Σ(vcmodel, vcdata) - fisher_Σ(vcmodel, vcdatarot)) ≈ 0.0
@test vecnorm(fisher_Σ(vcmodel, [vcdata vcdata]) -
  2fisher_Σ(vcmodel, vcdata)) ≈ 0.0
@test vecnorm(fisher_Σ(vcmodel, [vcdata vcdata]) -
  fisher_Σ(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

info("Evaluate Fisher information matrix of B")
H = zeros(p * d, p * d)
#@code_warntype fisher_B!(H, vcmodelrot, vcdatarot)
@inferred fisher_B!(H, vcmodelrot, vcdatarot)
@test vecnorm(fisher_B(vcmodel, vcdata) - fisher_B(vcmodelrot, vcdatarot)) ≈ 0.0
@test vecnorm(fisher_B(vcmodel, vcdata) - fisher_B(vcmodel, vcdatarot)) ≈ 0.0
@test vecnorm(fisher_B(vcmodel, [vcdata vcdata]) -
  2.0fisher_B(vcmodel, vcdata)) ≈ 0.0
@test vecnorm(fisher_B(vcmodel, [vcdata vcdata]) -
  fisher_B(vcmodelrot, [vcdatarot vcdatarot])) ≈ 0.0

info("Find MLE using Fisher scoring")
vcmfs = deepcopy(vcmodel)
#@code_warntype mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
#@inferred mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
logl_fs, _, _, Σcov_fs, Bse_fs, = mle_fs!(vcmfs, vcdatarot; solver = :Ipopt)
logl_fs_array, = mle_fs!(vcmfs, [vcdatarot vcdatarot]; solver = :Ipopt)
@show vcmfs.B
@show Bse_fs
@show B

info("Find MLE using MM algorithm")
vcmmm = deepcopy(vcmodel)
#@code_warntype mle_mm!(vcmm, vcdatarot)
#@inferred mle_mm!(vcmmm, vcdatarot)
logl_mm, _, _, _, Bse_mm, = mle_mm!(vcmmm, vcdatarot)
@test abs(logl_fs - logl_mm) / (abs(logl_fs) + 1.0) < 1.0e-4
@show vcmmm.B
@show Bse_mm
@show B
vcmmm = deepcopy(vcmodel)
logl_mm_array, = mle_mm!(vcmmm, [vcdatarot vcdatarot])
@test abs(logl_fs_array - logl_mm_array) / (abs(logl_fs_array) + 1.0) < 1.0e-4

info("Find MLE using Fisher scoring (linear equality + box constraints)")
vcmfs       = deepcopy(vcmodel)
vcmfs.A     = [1.0 -1.0 zeros(1, p*d-2)]
vcmfs.sense = '='
vcmfs.b     = 0.0
vcmfs.lb    = 0.0
vcmfs.ub    = 1.0
logl_fs, _, _, Σcov_fs, Bse_fs, = mle_fs!(vcmfs, vcdatarot; solver = :Ipopt, qpsolver = :Ipopt)
@show vcmfs.B
@test vcmfs.B[1] ≈ vcmfs.B[2]
@test all(vcmfs.B .≥ 0.0)
@test all(vcmfs.B .≤ 1.0)

info("Find MLE using MM algorithm (linear equality + box constraints)")
vcmm = deepcopy(vcmodel)
vcmm.A = [1.0 -1.0 zeros(1, p*d-2)]
vcmm.sense = '='
vcmm.b = 0.0
vcmm.lb = 0.0
vcmm.ub = 1.0
logl_mm, _, _, Σcov_mm = mle_mm!(vcmm, vcdatarot; qpsolver = :Ipopt)
@show vcmm.B
@test vcmm.B[1] ≈ vcmm.B[2]
@test all(vcmm.B .≥ 0.0)
@test all(vcmm.B .≤ 1.0)
@test abs(logl_fs - logl_mm) / (abs(logl_fs) + 1.0) < 1.0e-4

info("Heritability estimation")
h, h_se = heritability(vcmfs.Σ, Σcov_fs)
@show h, h_se
@test all(0.0 .≤ h .≤ 1.0)

info("test fit_mle (FS)")
vcmmle = deepcopy(vcmodel)
logl_mle, _, _, Σcov_mle, Bse_mle, = fit_mle!(vcmmle, vcdata; algo = :FS)
@show vcmmle.B, Bse_mle, B
vcmle = deepcopy(vcmodel)
fit_mle!(vcmmle, [vcdata vcdata]; algo = :FS)

info("test fit_mle (MM)")
vcmmle = deepcopy(vcmodel)
logl_mle, _, _, Σcov_mle, Bse_mle, = fit_mle!(vcmmle, vcdata; algo = :MM)
@show vcmmle.B, Bse_mle, B

info("test fit_reml (FS)")
vcmreml = deepcopy(vcmodel)
logl_reml, _, _, Σcov_reml, Bse_reml, = fit_reml!(vcmreml, vcdata; algo = :FS)
@show vcmreml.B, Bse_reml, B

info("test fit_reml (MM)")
vcmreml = deepcopy(vcmodel)
logl_reml, _, _, Σcov_reml, Bse_reml, = fit_reml!(vcmreml, vcdata; algo = :MM)
@show vcmreml.B, Bse_reml, B

end # module VarianceComponentTypeTest


text:
using Pkg, DrWatson, Random, Test

Random.seed!(123)

@quickactivate "StatisticalRethinkingTuring"
using Turing
using StatisticalRethinking

#=
using Turing, DataFrames
=#

Turing.setprogress!(false);

@model function linear_reg(x, y)
    β ~ Normal(0, 1)
    α ~ Normal(0, 5)
    σ ~ Exponential(1)

    for i ∈ eachindex(y)
        y[i] ~ Normal(α + β * x[i], σ)
    end
    return (α, β, σ)
end

f(x) = 5 + 2 * x + 0.5 * randn()

Δ = 1.0; xs_train = 0:Δ:20; ys_train = f.(xs_train);
xs_test = [10 + i * Δ for i in 1:4]; ys_test = f.(xs_test);

# Infer
m_train = linear_reg(xs_train, ys_train);
#precis(m_train) |> display  # From StatisticalRethinking.precis

chain_lin_reg = sample(m_train, NUTS(100, 0.65), 200);
chain_lin_reg |> display

#m_test = linear_reg(xs_test, fill(missing, length(ys_test)));
m_test = linear_reg(xs_test, similar(xs_test, Missing));
#m_test = linear_reg(xs_test, missing);

# Use the new predict function!
predictions = predict(m_test, chain_lin_reg)
predictions |> display

# Get the mean predicted values.
ys_mean_pred = vec(mean(Array(group(predictions, :y)); dims = 1))
ys_std_pred = vec(std(Array(group(predictions, :y)); dims = 1))

# Get the prediction error:
ys_test - ys_mean_pred |> display
ys_test - ys_mean_pred2 |> display

#@test sum(abs2, ys_test - ys_pred) ≤ σ


text:
"""
The Node.jl module defines the Node type and methods

------------------------------------------------------------------------------
COPYRIGHT/LICENSE. This file is part of the XYZ package. It is subject to
the license terms in the LICENSE file found in the top-level directory of
this distribution. No part of the XYZ package, including this file, may be
copied, modified, propagated, or distributed except according to the terms
contained in the LICENSE file.
------------------------------------------------------------------------------
"""
# --- Exports

# ------ Types

export Node

# ------ Functions

import Base.run

# --- Type definitions

using ZMQ

"""
    struct Node

TODO
"""
struct Node
    #=
      Fields
      ------
      * __Attributes__

        * `id`: id

      * __Components__

        * `processing_core`: processing core

        * `control_unit`: control unit

        * `output_channel`: channel for publishing output

        * `input_channels`: channels for listening for inputs

        * `input_data`: TODO
    =#

    # Attributes
    id::String

    # Components
    processing_core::AbstractProcessingCore
    control_unit::ControlUnit
    output_channel::OutputChannel
    input_channels::Vector{InputChannel}
    input_data::Vector
end

# --- Constructors

"""
    Node(id::String,
         processing_core::AbstractProcessingCore,
         data_type::Type{<:AbstractNodeData},
         control_unit_params::Dict,
         output_channel_params::Dict,
         input_channel_params::Vector{Dict})

Construct a distributable, autonomous compute node identified by `id` and that
uses `processing_core` to process data that it receives. `control_unit_params`,
`output_channel_params`, and `input_channel_params` define the Node's control
unit, output channel, and input channels, respectively.

### ControlUnit Parameters

* `logic_core`<:AbstractControlLogicCore - control logic unit for Node

* `url`::String - URL that the ControlUnit listens for control signals on

* `copy_logic_core`::Bool - set to true if the ControlUnit should be
  initialized with a copy of `logic_core`; set to false if the ControlUnit
  is set to point to the instance referenced by `logic_core`.

* `use_bind`::Bool - set to true if the ZMQ Socket that the ControlUnit
  listens for control signals on should be connected to `url` using the
  `bind()` method; set to false if the ZMQ Socket should be connected using
  the `connect()` method. For the ControlUnit, defaults to `true`.

### OutputChannel Parameters

* `data_type`::Type - subtype of AbstractChannelData that output channel
  publishes

* `url`::String - URL that the OutputChannel publishes data to

* `use_bind`::Bool - set to true if the ZMQ Socket that the OutputChannel
  publishes data on should be connected to `url` using the `bind()` method;
  set to false if the ZMQ Socket should be connected using the `connect()`
  method. For the OutputChannel, defaults to `true`.

### Input Channel Parameters

* `data_type`::Type - subtype of AbstractChannelData that input channel accepts

* `url`::String - URL that the InputChannel listens for data on

* `use_bind`::Bool - set to true if the ZMQ Socket that the InputChannel
  listens for data on should be connected to `url` using the `bind()` method;
  set to false if the ZMQ Socket should be connected using the `connect()`
  method. For the InputChannel, defaults to `false`.
"""
function Node(id::String,
              processing_core::AbstractProcessingCore,
              control_unit_params::Dict,
              output_channel_params::Dict,
              input_channel_params::Vector{Dict})

    # --- Check arguments

    # Check for required parameters for control unit
    required_params = ["logic_core", "url"]
    for param in required_params
        if !haskey(control_unit_params, param)
            message = "`control_unit_params` is missing required parameter: " *
                      "`$param`"
            throw(ArgumentError(message))
        end
    end

    # Check for required parameters for output channel
    required_params = ["data_type", "url"]
    for param in required_params
        if !haskey(output_channel_params, param)
            message =
                "`output_channel_params` is missing required parameter: " *
                "`$param`"
            throw(ArgumentError(message))
        end
    end

    # Check for required parameters for input channels
    required_params = ["data_type", "url"]
    for i in 1:length(input_channel_params)
        for param in required_params
            if !haskey(input_channel_params[i], param)
                message =
                    "`input_channel_params[$i]` is missing required " *
                    "parameter: `$param`"
                throw(ArgumentError(message))
            end
        end
    end

    # --- Construct components

    # Construct control unit
    control_unit =
        ControlUnit(control_unit_params["logic_core"],
                    control_unit_params["url"],
                    copy_logic_core=get(control_unit_params,
                                        "copy_logic_core", true),
                    use_bind=get(control_unit_params, "use_bind", true))

    # Construct output channel
    output_channel =
        OutputChannel(output_channel_params["data_type"],
                      output_channel_params["url"],
                      use_bind=get(output_channel_params, "use_bind", true))

    # Construct input channels
    input_channels = Vector{InputChannel}()
    for params in input_channel_params
        input_channel =
            InputChannel(params["data_type"],
                         params["url"],
                         use_bind=get(params, "use_bind", true))

        push!(input_channels, input_channel)
    end

    # Construct storage for input data
    # TODO: initialize to default value
    input_data = Vector(undef, length(input_channels))

    # --- Construct new Node

    Node(id, processing_core,
         control_unit, output_channel, input_channels, input_data)
end

# --- Method definitions

# Constants
INIT_SLEEP_TIME = 0.001

"""
    run(node::Node)

Activate all of the operational components of `node` and start running the
main processing loop.
"""
function run(node::Node)
    # --- Start listening on all input channels

    for channel in node.input_channels
        loop_listen(channel)
    end

    # --- Start control unit

    # TODO

    # --- TODO: Do we need to verify input connections?

    # TODO: Consider this alternative as simplification to async implementation
    #
    # --- Set control unit to be in 'running' state
    #
    # set_running(node.control_unit)
    #
    # --- Main processing loop
    #
    # set_running!(node.control_unit)
    #
    # while is_running(node)
    #     channel = Channel()
    #     @async process_control_signal(); put!(channel, :received)
    #     sleep(CONTROL_CHANNEL_TIMEOUT)
    #     close(channel)
    #
    #     process_inputs()
    # end

    # --- Start control unit

    # Wait for control unit to be in 'running' state
    while !is_running(node)
        sleep(INIT_SLEEP_TIME)
    end

    # --- Run main processing loop

    run_processing_loop!(node)
end

"""
    run_processing_loop!(node::Node)

Run main processing loop.
"""
function run_processing_loop!(node::Node)
    while is_running(node)
        # Wait for input to be ready for processing
        wait_for_input_ready(node.control_unit, node.input_channels)

        # Gather input
        fetch_input!(node)

        # Process input
        result = process_data!(node.processing_core, node.input_data)

        # Publish output
        publish(node.output_channel, result)

        # Sleep to give other tasks a chance to run
        # TODO: improve design and implementation. One possibility would
        #       be to use round-robin between control and input processing
        sleep(0.001)
    end
end

"""
    is_running(node::Node)

Return true if the Node is running; return false otherwise.
"""
is_running(node::Node) = is_running(node.control_unit)

"""
    is_terminated(node::Node)

Return true if the Node is terminated; return false otherwise.
"""
is_terminated(node::Node) = is_terminated(node.control_unit)

"""
    fetch_input(node::Node)

Gather input from the input channels of `node`.
"""
function fetch_input!(node::Node)
    for (i, channel) in enumerate(node.input_channels)
        node.input_data[i] = get_last_value(channel)
    end
end


text:
export Chain

"""
    Chain{B} <: BCLattice{B, 1}

general chain lattice with boundary condition `B`.

    Chain([Boundary], length)

Construct a chain lattice with boundary (optional,
default to be `Fixed`)

```julia-repl
julia> Chain(10)
QMTK.Chain{QMTK.Fixed}(10)

```
"""
struct Chain{B<:Boundary} <: BCLattice{B, 1}
    length::Int
end

Chain(length::Integer) = Chain{Fixed}(length)
Chain(::Type{B}, length::Integer) where {B<:Boundary} = Chain{B}(length)

shape(chain::Chain) = (chain.length, )
length(chain::Chain) = chain.length
sites(chain::Chain) = ChainSiteIter(chain)
bonds(chain::Chain, bond::Int) = ChainBondIter(chain, bond)

struct ChainSiteIter{B<:Boundary} <: SiteIterator{Chain{B}}
    length::Int
end

ChainSiteIter(chain::Chain{B}) where B = ChainSiteIter{B}(chain.length)

struct ChainBondIter{B<:Boundary, K} <: BondIterator{Chain{B}}
    length::Int
end

ChainBondIter(chain::Chain{B}, bond::Int) where {B<:Boundary} =
    ChainBondIter{B, bond}(chain.length)

lattice(itr::LatticeIterator{Chain{B}}) where B = Chain(B, itr.length)

############
# Iterators
############

# import iterator interface
import Base: start, next, done, length, eltype

start(itr::ChainSiteIter) = 1
next(itr::ChainSiteIter, state::Int) = (state, state+1)
done(itr::ChainSiteIter, state::Int) = state > itr.length
length(itr::ChainSiteIter) = itr.length
eltype(itr::ChainSiteIter) = Int

# Bond Iterators
eltype(itr::ChainBondIter) = Tuple{Int, Int}
start(itr::ChainBondIter) = 1

next(itr::ChainBondIter{Fixed, K}, state::Int) where K = ((state, state + K), state+1)
done(itr::ChainBondIter{Fixed, K}, state::Int) where K = state > (itr.length - K)
length(itr::ChainBondIter{Fixed, K}) where K = itr.length - K

next(itr::ChainBondIter{Periodic, K}, state::Int) where K = ((state, (state+K-1)%itr.length+1), state+1)
done(itr::ChainBondIter{Periodic, K}, state::Int) where K = state > itr.length
length(itr::ChainBondIter{Periodic, K}) where K = itr.length

getid(chain::Chain, pos::Integer) = pos


text:
module BracketInserter

using OhMyREPL

# If a user writes an opening bracket
# automatically complete it with a closing bracket
# unless the next character is that closing bracket
mutable struct BracketInserterSettings
    complete_brackets::Bool
end

import REPL
import REPL.LineEdit
import REPL.LineEdit: edit_insert, edit_move_left, edit_move_right, buffer, char_move_left,
                      edit_backspace, terminal, transition, state

import REPL.Terminals.beep
import OhMyREPL.Prompt.rewrite_with_ANSI

const BRACKET_INSERTER = BracketInserterSettings(true)

function leftpeek(b::IOBuffer)
    p = position(b)
    c = char_move_left(b)
    seek(b, p)
    return c
end

function peek(b::IOBuffer)
    p = position(b)
    c = read(b, Char)
    seek(b, p)
    return c
end

const AUTOMATIC_BRACKET_MATCH = Ref(!Sys.iswindows())
enable_autocomplete_brackets(v::Bool) = AUTOMATIC_BRACKET_MATCH[] = v

const pkgmode = Ref{Any}()
import Pkg
function insert_into_keymap!(D::Dict)
    left_brackets = ['(', '{', '[']
    right_brackets = [')', '}', ']']
    right_brackets_ws = vcat(right_brackets, [' ', '\t'])
    left_brackets_ws = vcat(left_brackets, [' ', '\t'])
    all_brackets_ws = vcat(left_brackets, right_brackets_ws)
    for (l, r) in zip(left_brackets, right_brackets)
        # If we enter a left bracket automatically complete it if the next
        # char is a whitespace or a right bracket
        D[l] = (s, o...) -> begin
            edit_insert(buffer(s), l)
            if AUTOMATIC_BRACKET_MATCH[] && (eof(buffer(s)) || peek(buffer(s)) in right_brackets_ws)
                edit_insert(buffer(s), r)
                edit_move_left(buffer(s))
            end
            rewrite_with_ANSI(s)
        end
        # If we enter a right bracket and the next char is that right bracket just move right
        D[r] = (s, o...) -> begin
            if AUTOMATIC_BRACKET_MATCH[] && !eof(buffer(s)) && peek(buffer(s)) == r
                edit_move_right(buffer(s))
            else
                edit_insert(buffer(s), r)
            end
            rewrite_with_ANSI(s)
        end
    end

    f = D[']']
    D[']'] = (s, o...) -> begin
        if isempty(s) || position(LineEdit.buffer(s)) == 0
            if !isassigned(pkgmode)
                found_pkg = false
                for mode in Base.active_repl.interface.modes
                    if mode isa LineEdit.Prompt
                        if mode.prompt == Pkg.REPLMode.promptf
                            found_pkg = true
                            pkgmode[] = mode
                        end
                    end
                end
                if !found_pkg
                    pkgmode[] = nothing
                end
            end
            if pkgmode[] !== nothing
                buf = copy(LineEdit.buffer(s))
                transition(s, pkgmode[]) do
                    LineEdit.state(s, pkgmode[]).input_buffer = buf
                end
            else
                f(s, o...)
            end
        else
            f(s, o...)
        end
    end

    # Similar to above but with quotation marks that need to be handled a bit differently
    for v in ['\"', '\'']
        D[v] = (s, o...) -> begin
            b = buffer(s)
            # Next char is the quote symbol so just move right
            if AUTOMATIC_BRACKET_MATCH[] && !eof(b) && peek(b) == v
                edit_move_right(buffer(s))
            elseif AUTOMATIC_BRACKET_MATCH[] &&
                    ((position(b) > 0 && leftpeek(b) in left_brackets_ws && (eof(b) || peek(b) in right_brackets_ws)) ||
                     ((!eof(b) && peek(b) in right_brackets_ws) || b.size == 0) && (position(b) == 0 || leftpeek(b) in left_brackets_ws))
                edit_insert(buffer(s), v)
                edit_insert(buffer(s), v)
                edit_move_left(buffer(s))
            else
                edit_insert(buffer(s), v)
            end
            rewrite_with_ANSI(s)
        end
    end

    # On backspace, also remove a corresponding right bracket
    # to the right if we remove a left bracket
    left_brackets2 = ['(', '{', '[', '\"', '\'']
    right_brackets2 = [')', '}', ']', '\"', '\'']
    D['\b'] = (s, data, c) -> begin
        repl = Base.active_repl
        mirepl = isdefined(repl,:mi) ? repl.mi : repl
        main_mode = mirepl.interface.modes[1]
        if isempty(s) || position(buffer(s)) == 0
            buf = copy(buffer(s))
            transition(s, main_mode) do
                state(s, main_mode).input_buffer = buf
            end
        else
            b = buffer(s)
            str = String(take!(copy(b)))
            if AUTOMATIC_BRACKET_MATCH[] && !eof(buffer(s)) && position(buffer(s)) != nothing
                i = findfirst(isequal(str[prevind(str, position(b) + 1)]), left_brackets2)
                if i != nothing && peek(b) == right_brackets2[i]
                    edit_move_right(buffer(s))
                    edit_backspace(buffer(s))
                    edit_backspace(buffer(s))
                    rewrite_with_ANSI(s)
                    return
                end
            end
            edit_backspace(buffer(s))
        end
        rewrite_with_ANSI(s)
    end
end

insert_into_keymap!(OhMyREPL.Prompt.NEW_KEYBINDINGS)

end # module


text:
using CompScienceMeshes, BEAST
o, x, y, z = euclidianbasis(3)

trc = X->ntrace(X,γ)
dvg = divergence

h1 = h2 = h3 = 1/15; h = max(h1, h2, h3)
width, height = 1.0, 0.5
Γ1 = meshrectangle(width, height, h1)
Γ2 = meshrectangle(width, height, h2); CompScienceMeshes.rotate!(Γ2, 0.5π * [1,0,0])
Γ3 = meshrectangle(width, height, h3); CompScienceMeshes.rotate!(Γ3, 1.0π * [1,0,0])
γ = meshsegment(width, width, 3)

κ = 1.0
S, T, I = SingleLayerTrace(κ), MWSingleLayer3D(κ), Identity()
E = planewavemw3d(direction=z, polarization=x, wavenumber=κ)
e = (n×E)×n

X1, X2, X3 = raviartthomas(Γ1, γ), raviartthomas(Γ2, γ), raviartthomas(Γ3, γ)
X = X1 × X2 × X3

@hilbertspace j
@hilbertspace k
α, β = 1/(im*κ), log(abs(h))/(im*κ)
Eq = @varform T[k,j] + α*S[trc(k), dvg(j)] + α*S.'[dvg(k), trc(j)] + β*I[trc(k), trc(j)] == e[k]
eq = @discretise Eq  j∈X k∈X
u = solve(eq)

fcr, geo = facecurrents(u, X)
include(Pkg.dir("CompScienceMeshes","examples","plotlyjs_patches.jl"))
p = patch(geo, real.(norm.(fcr)))
PlotlyJS.plot(p)


text:
function value_of(::Type{SignalCode}, arg0::JString)
    return jcall(SignalCode, "valueOf", SignalCode, (JString,), arg0)
end

function values(::Type{SignalCode})
    return jcall(SignalCode, "values", Vector{SignalCode}, ())
end


text:
"""
    QDHT{p, n}(R, N; dim=1)
    QDHT([p, ] R, N; dim=1)
    QDHT(p, [n, ] R, N; dim=1)

`p`-th order quasi-discrete Hankel transform over aperture radius `R` with `N` samples
which transforms along dimension `dim`. If not given, `p` defaults to 0, and `n` the
spherical dimension defaults to 1 (cylindrical).

After:

[1] L. Yu, M. Huang, M. Chen, W. Chen, W. Huang, and Z. Zhu, Optics Letters 23 (1998)

[2] M. Guizar-Sicairos and J. C. Gutiérrez-Vega, JOSA A 21, 53 (2004)

[3] H. F. Johnson, Comput. Phys. Commun., 43, 2 (1987)

but with some alterations:

The transform matrix T is not the same as C/T defined in [1, 2] but is more like the form
used in equation 14 of [3].
Instead of dividing by ``J_1(α_{pn}) J_1(α_{pm})`` we divide by ``J_1(α_{pn})^2``.
This cancels out the factor between ``f`` and ``F`` so we do not have to mutltiply (divide)
by ``J_1(α_{pn})`` (``J_1(α_{pm})``) before and after applying the transform matrix.

Follows [`AbstractFFT`](https://github.com/JuliaMath/AbstractFFTs.jl) approach of applying
fwd and inv transform with `mul` and `ldiv`.

To calculate radial integrals of functions sampled using `QDHT`, use [`integrateR`](@ref)
and [`integrateK`](@ref).

The type of the coefficients is inferred from the type of `R` (but is promoted to be at
least `Float`), so for arbitrary precision use `QDHT([p, ] BigFloat(R), ...)`.
"""
struct QDHT{p, n, T<:Real}
    N::Int # Number of samples
    T::Array{T, 2} # Transform matrix
    j1sq::Array{T, 1} # (j₁)² factors
    K::T # Highest spatial frequency
    k::Vector{T} # Spatial frequency grid
    R::T # Aperture size (largest real-space coordinate)
    r::Vector{T} # Real-space grid
    scaleR::Vector{T} # Scale factor for real-space integration
    scaleK::Vector{T} # Scale factor for frequency-space integration
    scaleRK::T # Scalar factor multiplied/divided during fwd/rev transform
    dim::Int # Dimension along which to transform
end

function QDHT{p, n}(R, N; dim=1) where {p, n}
    R = float(R)
    pf = oftype(R, p)
    cn = oftype(R, sphbesselj_scale(n))
    roots = sphbesselj_zero.(pf, n, 1:N) # type of sphbesselj_zero is inferred from first argument
    S = sphbesselj_zero(pf, n, N+1)
    r = roots .* R/S # real-space vector
    K = S/R # Highest spatial frequency
    k = roots .* K/S # Spatial frequency vector
    j₁sq = abs2.(sphbesselj.(p+1, n, roots))
    T = 2 * cn * S^(-(n + 1) / 2) * sphbesselj.(p, n, (roots .* roots') ./ S) ./ j₁sq' # Transform matrix

    scaleR = 2 * cn^2 / K^(n + 1) ./ j₁sq # scale factor for real-space integration
    scaleK = 2 * cn^2 / R^(n + 1) ./ j₁sq # scale factor for reciprocal-space integration
    scaleRK = (R / K) ^ ((n + 1) / 2)
    QDHT{p, n, eltype(T)}(N, T, j₁sq, K, k, R, r, scaleR, scaleK, scaleRK, dim)
end

QDHT{p}(R, N; dim=1) where {p} = QDHT{p, 1}(R, N; dim=dim)
QDHT(R, N; dim=1) = QDHT{0}(R, N; dim=dim)
QDHT(p, R, N; dim=1) = QDHT{p}(R, N; dim=dim)
QDHT(p, n, R, N; dim=1) = QDHT{p, n}(R, N; dim=dim)

"""
    order(q::QDHT)

Order of transform `q`.
"""
order(::QDHT{p}) where {p} = p

"""
    sphericaldim(q::QDHT)

Spherical dimension ``n`` of transform `q`.

Dimension of the ``n``-sphere over whose radial coordinate the bases of the transform `q`
are orthogonal polynomials. For the circle, ``n=1``, and for the sphere, ``n=2``.

See [`sphbesselj`](@ref).
"""
sphericaldim(::QDHT{p,n}) where {p,n} = n

"""
    inv(Q::QDHT) -> QDHT

Compute the `QDHT` for the inverse of the Hankel transform `Q`, that is, the transform such
that `Q * A == inv(Q) \\ A`
"""
function Base.inv(Q::T) where {T<:QDHT}
    return T(
        Q.N,
        Q.T,
        Q.j1sq,
        Q.R,
        Q.r,
        Q.K,
        Q.k,
        Q.scaleK,
        Q.scaleR,
        inv(Q.scaleRK),
        Q.dim,
    )
end

"""
    reshape(Q::QDHT, N::Int) -> QDHT

Compute a new `QDHT` with the same order, spherical dimension, and radius but with `N`
points.
"""
function Base.reshape(Q::QDHT{p, n}, N::Int) where {p, n}
    N == Q.N && return Q
    return QDHT{p, n}(Q.R, N, dim=Q.dim)
end

"
    mul!(Y, Q::QDHT, A)

Calculate the forward quasi-discrete Hankel transform of array `A` using the QDHT `Q`
and store the result in `Y`.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R)); Y = similar(A);
julia> mul!(Y, q, A)
8-element Array{Float64,1}:
  4.326937831591551e-6
  2.3341589529175126e-6
  7.689558743828849e-7
  1.546419420523699e-7
  1.8999259906096856e-8
  1.4159642663129888e-9
  7.013670190083954e-11
 -6.07681871673291e-13
```
"
function mul!(Y, Q::QDHT, A)
    dot!(Y, Q.T, A, dim=Q.dim)
    Y .*= Q.scaleRK
end

"
    ldiv!(Y, Q::QDHT, A)

Calculate the inverse quasi-discrete Hankel transform of array `A` using the QDHT `Q`
and store the result in `Y`.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R)); Y = similar(A);
julia> mul!(Y, q, A);
julia> YY = similar(Y); ldiv!(YY, q, Y);
julia> YY ≈ A
true
```
"
function ldiv!(Y, Q::QDHT, A)
    dot!(Y, Q.T, A, dim=Q.dim)
    Y ./= Q.scaleRK
end

"""
    *(Q::QDHT, A)

Calculate the forward quasi-discrete Hankel transform of array `A` using the QDHT `Q`.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R));
julia> q*A
8-element Array{Float64,1}:
  4.326937831591551e-6
  2.3341589529175126e-6
  7.689558743828849e-7
  1.546419420523699e-7
  1.8999259906096856e-8
  1.4159642663129888e-9
  7.013670190083954e-11
 -6.07681871673291e-13
```
"""
function *(Q::QDHT, A)
    out = similar(A)
    mul!(out, Q, A)
end

"""
    \\(Q::QDHT, A)

Calculate the inverse quasi-discrete Hankel transform of array `A` using the QDHT `Q`.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(1e-2, 8); A = exp.(-q.r.^2/(1e-3*q.R));
julia> Ak = q*A;
julia> q \\ Ak ≈ A
true
```
"""
function \(Q::QDHT, A)
    out = similar(A)
    ldiv!(out, Q, A)
end

"""
    integrateR(A, Q::QDHT; dim=1)

Radial integral of `A`, over the aperture of `Q` in real space.

Assuming `A` contains samples of a function `f(r)` at sample points `Q.r`, then
`integrateR(A, Q)` approximates ∫f(r)r dr from r=0 to r=∞.

!!! note
    `integrateR` and `integrateK` fulfill Parseval's theorem, i.e. for some array `A`,
    `integrateR(abs2.(A), q)` and `integrateK(abs2.(q*A), q)` are equal, **but**
    `integrateR(A, q)` and `integrateK(q*A, q)` are **not** equal.

!!! warning
    using `integrateR` to integrate a function (i.e. `A` rather than `abs2(A)`) is only
    supported for the 0th-order QDHT. For more details see [Derivations](@ref).

# Examples
```jldoctest
julia> q = QDHT{0, 1}(10, 128); A = exp.(-q.r.^2/2);
julia> integrateR(abs2.(A), q) ≈ 0.5 # analytical solution of ∫exp(-r²)r dr from 0 to ∞
true
```
"""
integrateR(A, Q::QDHT; dim=1) = dimdot(Q.scaleR, A; dim=dim)

"""
    integrateK(Ak, Q::QDHT; dim=1)

Radial integral of `A`, over the aperture of `Q` in reciprocal space.

Assuming `A` contains samples of a function `f(k)` at sample points `Q.k`, then
`integrateR(A, Q)` approximates ∫f(k)k dk from k=0 to k=∞.

!!! note
    `integrateR` and `integrateK` fulfill Parseval's theorem, i.e. for some array `A`,
    `integrateR(abs2.(A), q)` and `integrateK(abs2.(q*A), q)` are equal, **but**
    `integrateR(A, q)` and `integrateK(q*A, q)` are **not** equal.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(10, 128); A = exp.(-q.r.^2/2);
julia> integrateR(abs2.(A), q) ≈ 0.5 # analytical solution of ∫exp(-r²)r dr from 0 to ∞
true
julia> Ak = q*A;
julia> integrateK(abs2.(Ak), q) ≈ 0.5 # Same result
true

```
"""
integrateK(Ak, Q::QDHT; dim=1) = dimdot(Q.scaleK, Ak; dim=dim)

"""
    onaxis(Ak, Q::QDHT; dim=Q.dim)

Calculate on-axis sample in space (i.e. at r=0) from transformed array `Ak`.

!!! note
    `onaxis` is currently only supported for 0-order transforms

# Examples
```jldoctest
julia> q = QDHT{0, 1}(10, 128); A = exp.(-q.r.^2/2);
julia> onaxis(q*A, q) ≈ 1 # should be exp(0) = 1
true
```
"""
function onaxis(Ak, Q::QDHT{0,n}; dim = Q.dim) where {n}
    j00ocn = inv(gamma((n + 1) / 2) * 2^((n - 1) / 2))
    return j00ocn .* integrateK(Ak, Q; dim = dim)
end
onaxis(Ak, Q::QDHT{0,1}; dim=Q.dim) = J₀₀ .* integrateK(Ak, Q; dim=dim)

"""
    symmetric(A, Q::QDHT)

Create symmetric array from samples in `A`, including on-axis sample.

Given `A`, sampled at `[r₁, r₂, r₃, ...]`, generates array sampled at
`[...-r₃, -r₂, -r₁, 0, r₁, r₂, r₃...]`

# Examples
```jldoctest
julia> q = QDHT{0, 1}(10, 128); A = exp.(-q.r.^2);
julia> As = symmetric(A, q);
julia> size(As)
(257,)
julia> As[1:128] == A[128:-1:1]
true
julia> As[129] ≈ 1 # should be exp(0) = 1
true
julia> As[130:end] == A
true
```
"""
function symmetric(A, Q::QDHT; dim=Q.dim)
    order(Q) == 0 || throw(DomainError("`symmetric` is only supported for 0-order transforms"))
    s = size(A)
    N = s[dim]
    idxlo = CartesianIndices(s[1:(dim - 1)])
    idxhi = CartesianIndices(s[(dim + 1):end])
    A0 = asarray(onaxis(Q * A, Q))
    out = cat(view(A, idxlo, N:-1:1, idxhi), A0, view(A, idxlo, :, idxhi); dims=dim)
    return out
end

"""
    Rsymmetric(Q::QDHT)

Create radial coordinate array to go along with `symmetric(A, Q::QDHT)`.

# Examples
```jldoctest
julia> q = QDHT{0, 1}(10, 4);
julia> q.r
4-element Array{Float64,1}:
 1.6106347946239767
 3.697078919099734
 5.795844623798052
 7.8973942990196395
julia> Rsymmetric(q)
9-element Array{Float64,1}:
 -7.8973942990196395
 -5.795844623798052
 -3.697078919099734
 -1.6106347946239767
  0.0
  1.6106347946239767
  3.697078919099734
  5.795844623798052
  7.8973942990196395
```
"""
Rsymmetric(Q::QDHT) = vcat(-view(Q.r, lastindex(Q.r):-1:firstindex(Q.r)), 0, Q.r)

"""
    oversample(A, Q::QDHT; factor::Int=4) -> Tuple{AbstractArray,QDHT}
    oversample(A, Q::QDHT, QNew::QDHT) -> AbstractArray

Oversample (smooth) the array `A`, which is sampled with the `QDHT` `Q`, by a `factor`.

If calling the function many times, it is more efficient to precompute `QNew` with
[`reshape`](@ref) and then provide it to the function.

This works like Fourier-domain zero-padding: a new `QDHT` is created with the same radius,
but `factor` times more points. The existing array is transformed and placed onto this
new spatial frequency grid, and the rest filled with zeros. Transforming back yields the
same shape in space but with more samples.

!!! note
    Unlike in zero-padding using FFTs, the old and oversampled **spatial** grids do not
    have any sampling points in common.
"""
function oversample(A, Q::QDHT; factor::Int=4)
    QNew = reshape(Q, Q.N * factor)
    return oversample(A, Q, QNew), QNew
end
function oversample(A, Q::QDHT, QNew::QDHT)
    Ak = Q * A
    Ako = padzeros(Ak, QNew.N; dim = QNew.dim)
    Ao = QNew \ Ako
    return Ao
end

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
import Pluto
Pluto.run(1234)