text:
# This package redefines Complex operations so that there is a single
# complex infinity and a single complex NaN.

import Base.inv, Base.Complex, Base.show, Base.showcompact
import Base.+, Base.-, Base.*, Base./, Base.==

function Complex(x::Real, y::Real)
    if isnan(x) || isnan(y)
        return ComplexNaN
    end

    if isinf(x) || isinf(y)
        return ComplexInf
    end

    if x==0
        x = zero(x)
    end
    if y==0
        y = zero(y)
    end

    return Complex(promote(x,y)...)
end

module RiemannComplexNumbers
"""
Complex not-a-number in the `RiemannComplexNumbers` module.
"""
const ComplexNaN = Complex(NaN,NaN)

"""
Single complex infinity in the `RiemannComplexNumbers` module.
"""
const ComplexInf = Complex(Inf,Inf)

export ComplexNaN, ComplexInf

function my_inv(z::Complex)
    if isnan(z)
        return ComplexNaN
    end
    if isinf(z)
        return zero(z)
    end
    if z == 0
        return ComplexInf
    end
    return my_inv_work(z)
end

my_inv_work(z::Complex)  = conj(z)/abs2(z)
my_inv_work{T<:Integer}(z::Complex{T}) = inv(float(z))

function my_inv_work(w::Complex128)
    c, d = reim(w)
    half = 0.5
    two = 2.0
    cd = max(abs(c), abs(d))
    ov = realmax(c)
    un = realmin(c)
    ep = eps(Float64)
    bs = two/(ep*ep)
    s = 1.0
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    cd <= un*two/ep && (c=c*bs; d=d*bs; s=s*bs     ) # scale up c,d
    if abs(d)<=abs(c)
        r = d/c
        t = 1.0/(c+d*r)
        p = t
        q = -r * t
    else
        c, d = d, c
        r = d/c
        t = 1.0/(c+d*r)
        p = r * t
        q = -t
    end
    return Complex128(p*s,q*s) # undo scaling
end

# z-->b w-->a
function my_div(a::Complex, b::Complex)
    a,b = promote(a,b)
    if isnan(a) || isnan(b)
        return ComplexNaN
    end
    if isinf(b)
        if isinf(a)
            return ComplexNaN
        end
        return zero(typeof(z=a))
    end
    if b==0
        if a==0
            return ComplexNaN
        end
        return ComplexInf
    end
    return my_div_work(a,b)
end

# Special cases are already handled before this is called. Code taken
# from Julia's complex.jl
function my_div_work(a::Complex, b::Complex)
    are = real(a); aim = imag(a); bre = real(b); bim = imag(b)
    if abs(bre) <= abs(bim)
        r = bre / bim
        den = bim + r*bre
        Complex((are*r + aim)/den, (aim*r - are)/den)
    else
        r = bim / bre
        den = bre + r*bim
        Complex((are + aim*r)/den, (aim - are*r)/den)
    end
end

# Also taken from complex.jl
function my_div_work(z::Complex128, w::Complex128)
    a, b = reim(z); c, d = reim(w)
    half = 0.5
    two = 2.0
    ab = max(abs(a), abs(b))
    cd = max(abs(c), abs(d))
    ov = realmax(a)
    un = realmin(a)
    ep = eps(Float64)
    bs = two/(ep*ep)
    s = 1.0
    ab >= half*ov  && (a=half*a; b=half*b; s=two*s ) # scale down a,b
    cd >= half*ov  && (c=half*c; d=half*d; s=s*half) # scale down c,d
    ab <= un*two/ep && (a=a*bs; b=b*bs; s=s/bs     ) # scale up a,b
    cd <= un*two/ep && (c=c*bs; d=d*bs; s=s*bs     ) # scale up c,d
    abs(d)<=abs(c) ? ((p,q)=robust_cdiv1(a,b,c,d)  ) : ((p,q)=robust_cdiv1(b,a,d,c); q=-q)
    return Complex128(p*s,q*s) # undo scaling
end

end #end of module

show(io::IO, z::Complex)        = print(io, string(z,false))
showcompact(io::IO, z::Complex) = print(io, string(z,true))

#### BASIC FOUR OPERATIONS ####

# Addition

function +(w::Complex, z::Complex)
    w,z = promote(w,z)
    if isnan(w) || isnan(z)
        return ComplexNaN
    end
    if isinf(w)
        if isinf(z)
            return ComplexNaN
        end
        return ComplexInf
    end
    if isinf(z)
        return ComplexInf
    end
    return Complex(w.re+z.re, w.im+z.im)
end

+(w::Complex, x::Real) = w + Complex(x)
+(x::Real, w::Complex) = Complex(x) + w

# Binary minus

function -(w::Complex, z::Complex)
    w,z = promote(w,z)
    if isnan(w) || isnan(z)
        return ComplexNaN
    end
    if isinf(w)
        if isinf(z)
            return ComplexNaN
        end
        return ComplexInf
    end
    if isinf(z)
        return ComplexInf
    end
    return Complex(w.re-z.re, w.im-z.im)
end

-(w::Complex, x::Real) = w - Complex(x)
-(x::Real, w::Complex) = Complex(x) - w

# Unary minus

function -(w::Complex)
    if isnan(w)
        return ComplexNan
    end
    if isinf(w)
        return ComplexInf
    end
    return Complex(-w.re, -w.im)
end

# Multiplication

function *(w::Complex, z::Complex)
    w,z = promote(w,z)
    if isnan(w) || isnan(z)
        return ComplexNaN
    end
    if isinf(w)
        if z==0
            return ComplexNaN
        end
        return ComplexInf
    end
    if isinf(z)
        if w==0
            return ComplexNaN
        end
        return ComplexInf
    end
    return Complex(w.re*z.re - w.im*z.im , w.re*z.im + w.im*z.re)
end

*(w::Complex, x::Real) = w * Complex(x)
*(x::Real, w::Complex) = Complex(x) * w

# Division

/(w::Complex, x::Real) = w * RiemannComplexNumbers.my_inv(Complex(x))
/(x::Real, z::Complex) = Complex(x)*RiemannComplexNumbers.my_inv(z)

/(w::Complex, z::Complex) = RiemannComplexNumbers.my_div(w,z)
/(w::Complex128, z::Complex128) = RiemannComplexNumbers.my_div(w,z)

# These cover the cases in complex.jl (Julia 0.3.10)
inv(w::Complex{Float64}) = RiemannComplexNumbers.my_inv(w)
inv{T<:Integer}(w::Complex{T}) = RiemannComplexNumbers.my_inv(w)
inv{T<:AbstractFloat}(w::Complex{T}) = RiemannComplexNumbers.my_inv(w)
inv(w::Complex) = RiemannComplexNumbers.my_inv(w)

# Equality

function ==(w::Complex, z::Complex)
    if isnan(w) || isnan(z)
        return false
    end
    if isinf(w) && isinf(z)
        return true
    end
    return w.re == z.re && w.im == z.im
end

==(w::Complex, x::Real) = w == Complex(x)
==(x::Real, z::Complex) = Complex(x) == z

import Base.string

function string(z::Complex,compact::Bool=true)
    if isnan(z)
        if compact
            return "C_NaN"
        else
            return "ComplexNaN"
        end
    end

    if isinf(z)
        if compact
            return "C_Inf"
        else
            return "ComplexInf"
        end
    end

    a,b = reim(z)

    # This is to hide -0.0
    if a==0
        a = zero(a)
    end

    if b==0
        b = zero(b)
    end

    sp = " "
    if compact
        sp = ""
    end

    op = "+"
    if b<0
        op = "-"
    end

    return string(a) * sp * op * sp * string(abs(b)) * "im"
end


text:
module FiberNlse
using FFTW
using ProgressBars

include("datatypes.jl")
# Physical units & constants

nm = ns = 1e-9
ps = pm = 1e-12
km = 1e3
mW = mm = 1e-3
GHz = 1e9
Thz = 1e12
m = 1
W = 1
c = 299792458

export nm, ns, ps, pm, km, mW, mm, GHz, THz, m, W, c

function configure(Nₜ::Int,
    Nₗ::Int,
    fib::Fiber,
    T,
    λ, conf::SimulationConfig=default_config)

    return configure(Nₜ::Int,
        Nₗ::Int,
        fib.D,
        fib.γ,
        fib.α,
        fib.L,
        T,
        λ, conf)

end

function configure(
    Nₜ::Int,
    Nₗ::Int,
    D,
    γ,
    α,
    L,
    T,
    λ, conf::SimulationConfig=default_config)
    t = T * 0.5 * range(-1, stop=1, length=Nₜ) # Time vector
    l = L * range(0,stop=1,length=Nₗ)
    Ψ = Matrix{ComplexF32}(zeros((Nₗ, Nₜ)))
    return Simulation(Ψ, L / Nₗ, T / Nₜ, Nₜ, Nₗ, D, γ, α, L, T, λ, -D * λ^2 / (2 * pi * c), t, conf), t, l
end

function inputSignal(sim, ψₒ)
    sim.Ψ[1, :] = ψₒ
end

function simulate(sim::Simulation, progress::Bool=false)

    D(sim, ν) = @. -sim.β2 * 0.5im * (2 * pi * ν)^2
    N̂(sim, u) = @. abs(u)^2 * sim.γ * 1im

    # Check wether to show progressbar or not
    if progress
        iter_z = ProgressBar(range(2, sim.Nₗ))
    else
        iter_z = range(2, sim.Nₗ)
    end

    ν = FFTW.fftfreq(sim.Nₜ, 1.0 / (sim.dt))

    dl = sim.dz
    α = sim.α

    D̂ = D(sim, ν)

    for i in iter_z

        ψₗ = sim.Ψ[i-1, :]

        ψₗ = ifft(exp.(0.5 * dl .* (D̂ .- 0.5α)) .* fft(ψₗ))
        ψₗ = exp.(dl * N̂(sim, ψₗ)) .* ψₗ
        ψₗ = ifft(exp.(0.5 * dl .* (D̂ .- 0.5α)) .* fft(ψₗ))

        sim.Ψ[i, :] = ψₗ

    end

end

function transition(sim1::Simulation, sim2::Simulation) # Use the final state of sim1 as initial state of sim2
    inputSignal(sim2, sim1.Ψ[end, :])
end

function rhs()
    # ! Lire le code de gnlse-python (les mettre dans les crédits)
end

function simulate2(sim::Simulation)
    # ! Integration RK4 de l'ode rhs
end
end

ext: jl


text:
Fk(w, x, u2k, c2k) = apply_quad(w, x, u2k) - c2k

Gk(w, x, u2kp1, c2kp1) = apply_quad(w, x, u2kp1) - c2kp1

function leastsquares_approximate(dict::Dictionary, points, values)
    @assert length(points) == length(values)

    M = length(points)
    N = min(length(values), max(6,round(Int, sqrt(M))))
    basis = rescale(resize(dict, N), minimum(points), maximum(points))
    A = BasisFunctions.evaluation_matrix(basis, points)
    Expansion(basis, A\values)
end

default_threshold(dict::Dictionary) = default_threshold(codomaintype(dict))

default_threshold(::Type{T}) where {T <: AbstractFloat} = sqrt(eps(T))/100
default_threshold(::Type{BigFloat}) = big(1e-28)

supportleft(dict) = leftendpoint(support(dict))
supportright(dict) = rightendpoint(support(dict))

function compute_one_point_rule(dict, moments; verbose = false)
    @assert length(dict) == 2
    @assert length(moments) == 2

    x0 = 1/2 * (supportleft(dict) + supportright(dict))
    w0 = moments[1] / eval_element(dict, 1, x0)
    w, x = compute_principal_representation_lower(dict, moments, w0, x0)
    # w0 = supportright(dict) - supportleft(dict)
    # w, x, converged = compute_one_point_rule(dict, moments, x0, w0,
    #     maxiter=maxiter, verbose=verbose)
    # if !converged
    #     x0 = 1/2 * (supportleft(dict) + supportright(dict))
    #     w0 = supportright(dict) - supportleft(dict)
    # end
end

# function compute_one_point_rule(dict, moments, x0, w0; maxiter, verbose)
#     sys = LowerPrincipalOdd(dict, moments)
#     result = newton(sys, quad_to_newton(sys, w0, x0),
#         maxiter=maxiter, verbose=verbose)
#     iterations = result[2]
#     w, x = newton_to_quad(sys, result[1])
#     w, x, iterations<maxiter
# end

function compute_canonical_representation_upper(dict, moments, xi, w0, x0)
    @assert iseven(length(dict))
    @assert length(moments) == length(dict)
    @assert length(w0) == (length(dict)>>1)+1

    rule = CanonicalRepresentationOdd_K1(dict, xi, moments)
    nx, iter, normx = newton_with_restart(rule, quad_to_newton(rule, w0, x0),
        threshold = default_threshold(dict), verbose = false)
    w, x = newton_to_quad(rule, nx)
end

function compute_many_canonical_representation_upper(dict, moments, gw, gx, n,
    pts = linspace(supportleft(dict), gx[1], n))

    l = length(dict) >> 1
    T = eltype(gw)
    w = zeros(T,l+1,n)
    x = zeros(T,l+1,n)
    for i = n:-1:1
        xi = pts[i]
        if i == n
            w0 = [gw; zero(T)]
            x0 = [gx; supportright(dict)]
        else
            w0 = w[:,i+1]
            x0 = x[:,i+1]
        end
        w1,x1 = compute_canonical_representation_upper(dict, moments, xi, w0, x0)
        w[:,i] = w1
        x[:,i] = x1
    end
    w,x
end

function chebyshev_approximation(basis, vals)
    A = approximation(basis)
    Expansion(basis, A*vals)
end

function approximate_canonical_representation_upper(dict, moments, gw, gx, n)
    @assert isodd(length(dict))
    @assert length(dict) == length(moments)
    l = (length(dict)-1) >> 1
    @assert length(gw) == l
    @assert length(gx) == l

    a = supportleft(dict)
    b = gx[1]
    a1 = 9a/10+b/10
    basis = ChebyshevT(n) → a1..b
    # basis = ChebyshevT(n) → supportleft(dict)..gx[1]
    pts = interpolation_grid(basis)
    w,x = compute_many_canonical_representation_upper(dict[1:2*l], moments[1:2*l], gw, gx, n, pts)
    T = eltype(w)
    Fvals = T[Fk(w[:,i],x[:,i],dict[2*l+1], moments[2*l+1]) for i in 1:size(w,2)]
    F_fun = chebyshev_approximation(basis, Fvals)
    w_funs = [chebyshev_approximation(basis, w[i,:]') for i in 1:size(w,1)]
    x_funs = [chebyshev_approximation(basis, x[i,:]') for i in 1:size(x,1)]
    F_fun, w_funs, x_funs
end

function compute_principal_representation_upper(dict, moments, w0, x0)
    rule = UpperPrincipalEven(dict, moments)
    nx, iter, normx = newton_with_restart(rule, quad_to_newton(rule, w0, x0),
        threshold = default_threshold(dict), verbose = false)
    w, x = newton_to_quad(rule, nx)
end

function compute_canonical_representation_lower(dict, moments, xi, w0, x0)
    @assert isodd(length(dict))
    @assert length(moments) == length(dict)
    @assert length(w0) == (length(dict)+1)>>1

    rule = CanonicalRepresentationEven_J1(dict, xi, moments)
    nx, iter, normx = newton_with_restart(rule, quad_to_newton(rule, w0, x0),
        threshold = default_threshold(dict), verbose = false, maxstep = 10)
    w, x = newton_to_quad(rule, nx)
end

function compute_many_canonical_representation_lower(dict, moments, gw, gx, n,
    pts = linspace(supportleft(dict), gx[1], n))

    l = (length(dict)+1) >> 1
    T = eltype(gw)
    w = zeros(T,l,n)
    x = zeros(T,l,n)
    w[:,n] = gw
    x[:,n] = gx
    for i = n:-1:1
        xi = pts[i]
        if i == n
            w0 = gw
            x0 = gx
        else
            w0 = w[:,i+1]
            x0 = x[:,i+1]
        end
        w1,x1 = compute_canonical_representation_lower(dict, moments, xi, w0, x0)
        w[:,i] = w1
        x[:,i] = x1
    end
    w,x
end

function approximate_canonical_representation_lower(dict, moments, gw, gx, n)
    @assert iseven(length(dict))
    @assert length(dict) == length(moments)
    l = length(dict) >> 1
    @assert length(gw) == l
    @assert length(gx) == l

    a = supportleft(dict)
    b = gx[1]
    a1 = 9a/10+b/10
    basis = ChebyshevT(n) → a1..b
    # basis = ChebyshevT(n) → supportleft(dict)..gx[1]
    pts = interpolation_grid(basis)
    w,x = compute_many_canonical_representation_lower(dict[1:2*l-1], moments[1:2*l-1], gw, gx, n, pts)
    T = eltype(w)
    Fvals = T[Fk(w[:,i],x[:,i],dict[2*l], moments[2*l]) for i in 1:size(w,2)]
    F_fun = chebyshev_approximation(basis, Fvals)
    w_funs = [chebyshev_approximation(basis, w[i,:]') for i in 1:size(w,1)]
    x_funs = [chebyshev_approximation(basis, x[i,:]') for i in 1:size(x,1)]
    F_fun, w_funs, x_funs
end

function compute_principal_representation_lower(dict, moments, w0, x0)
    rule = LowerPrincipalOdd(dict, moments)
    nx, iter, normx = newton_with_restart(rule, quad_to_newton(rule, w0, x0),
        threshold = default_threshold(dict), verbose = false)
    w, x = newton_to_quad(rule, nx)
end

function compute_gauss_rules(dict::Dictionary, moments = compute_moments(dict); verbose = false)
    @assert iseven(length(dict))

    n = length(dict)
    l = n >> 1

    T = codomaintype(dict)

    xk_upper = fill(T[],0)
    wk_upper = fill(T[],0)
    xk_lower = fill(T[],0)
    wk_lower = fill(T[],0)
    f_funs_upper = Array{Any}(undef,0)
    w_funs_upper = Array{Any}(undef,0)
    x_funs_upper = Array{Any}(undef,0)
    f_funs_lower = Array{Any}(undef,0)
    w_funs_lower = Array{Any}(undef,0)
    x_funs_lower = Array{Any}(undef,0)

    verbose && println("Computing initial one point rule")
    wk,xk = compute_one_point_rule(dict[1:2], moments[1:2], verbose=verbose)
    verbose && println("One point quadrature rule is: ", xk, ", ", wk)
    xi_lower = zeros(T,l)
    xi_upper = zeros(T,l-1)
    xi_lower[1] = xk[1]

    push!(wk_lower, wk)
    push!(xk_lower, xk)

    if l == 1
        return wk, xk, xi_upper, xi_lower, wk_lower, wk_upper, xk_lower, xk_upper, f_funs_lower, f_funs_upper, w_funs_lower, w_funs_upper, x_funs_lower, x_funs_upper
    end

    bn = 32
    for k = 1:l-1
        f_fun, w_funs, x_funs = approximate_canonical_representation_upper(
            dict[1:2*k+1], moments[1:2*k+1], wk, xk, bn)
        xi0 = bisection(f_fun, supportleft(f_fun), supportright(f_fun); threshold = default_threshold(T))
        w0 = T[w_funs[k](xi0) for k in 1:length(x_funs)]
        x0 = T[x_funs[k](xi0) for k in 1:length(x_funs)]
        wk, xk = compute_principal_representation_upper(dict[1:2*k+1], moments[1:2*k+1], w0, x0)
        xi = xk[1]
        xi_upper[k] = xi
        verbose && println("Upper principal representation ", k, " : xi is ", xi)

        # Store the results for posterity
        push!(wk_upper, wk)
        push!(xk_upper, xk)
        push!(f_funs_upper, f_fun)
        push!(w_funs_upper, w_funs)
        push!(x_funs_upper, x_funs)

        f_fun, w_funs, x_funs = approximate_canonical_representation_lower(
            dict[1:2*k+2], moments[1:2*k+2], wk, xk, bn)
        xi0 = bisection(f_fun, supportleft(f_fun), supportright(f_fun); threshold = default_threshold(T))
        w0 = T[w_funs[k](xi0) for k in 1:length(w_funs)]
        x0 = T[x_funs[k](xi0) for k in 1:length(x_funs)]
        wk, xk = compute_principal_representation_lower(dict[1:2*k+2], moments[1:2*k+2], w0, x0)
        xi = xk[1]
        xi_lower[k+1] = xi
        if verbose
            println("Lower principal representation ", k+1, " : xi is ", xi)
        end

        push!(wk_lower, wk)
        push!(xk_lower, xk)
        push!(f_funs_lower, f_fun)
        push!(w_funs_lower, w_funs)
        push!(x_funs_lower, x_funs)

        residual = abs(sum([wk[i]*dict[2*k+2](xk[i]) for i=1:length(wk)]) - moments[2*k+2])
        # verbose && @show residual

        if f_fun.coefficients[end] > 1e-10
            bn *= 2
            verbose && println("Increasing bn parameter to $(bn)")
        end
    end
    wk, xk, xi_upper, xi_lower, wk_lower, wk_upper, xk_lower, xk_upper, f_funs_lower, f_funs_upper, w_funs_lower, w_funs_upper, x_funs_lower, x_funs_upper
end

function compute_gauss_rule(dict::Dictionary, moments = compute_moments(dict); options...)
    wk, xk, xi_upper, xi_lower, wk_lower, wk_upper, xk_lower, xk_upper, f_funs_lower, f_funs_upper, w_funs_lower, w_funs_upper, x_funs_lower, x_funs_upper =
        compute_gauss_rules(dict, moments; options...)
    wk, xk
end


text:
using PrettyTables

function prettyprint(a, cnames="", rnames="")
if rnames !=""
    rnames = rnames[:]
    a = [rnames a]
    if cnames != ""
        cnames = cnames[:]
        cnames = vcat("", cnames)
    end
end
if cnames !=""
    pretty_table(a, cnames; formatters=ft_printf("%12.5f"))
else
    pretty_table(a; formatters=ft_printf("%12.5f"))
end
end

#=
#using Printf
# formatted print of array, with column names
function prettyprint(a, cnames="", rnames="",digits=12, decimals=5)
    # TBD: try to use this to allow using specified digits and decimals
    #fmt = @sprintf("%d",digits)"."@sprintf("%d",decimals)"%f"
    #@eval dofmt(x) = @sprintf($fmt, x)

    # print column names
    if cnames != ""
        for i = 1:size(a,2)
            pad = digits
            if rnames != "" && i==1
                pad = 2*digits
            end
            @printf("%s", lpad(cnames[i],pad," "))
        end
        @printf("\n")
    end
    # print the rows
    for i = 1:size(a,1)
        if rnames != ""
            @printf("%s", lpad(rnames[i],digits," "))
        end
        for j = 1:size(a,2)
            # TBD: use fmt defined above to print array contents
            @printf("%12.5f",(a[i,j]))
        end
        @printf("\n")
    end
    return
end
=#

max_issues_repo_name: DawievLill/Econometrics

max_forks_repo_name: DawievLill/Econometrics


text:
## Uniform

abstract type Uniform{T} <: Distribution{T} end

struct UniformType{T} <: Uniform{T} end

Uniform(::Type{T}) where {T} = UniformType{T}()

Base.getindex(::UniformType{T}) where {T} = T

struct UniformWrap{T,E} <: Uniform{E}
    val::T
end

Uniform(x::T) where {T} = Uniform(x, isuniform(T))

Uniform(x::T, ::Val{true}) where {T} = UniformWrap{T,gentype(T)}(x)

Base.getindex(x::UniformWrap) = x.val

const ImplicitUniform = Union{AbstractArray,AbstractDict,AbstractSet,
                              AbstractString,Tuple}

isuniform(::Type{T}) where {T} = Val(false)
isuniform(::Type{T}) where {T<:ImplicitUniform} = Val(true)

Base.show(io::IO, x::Uniform) = print(io, "Uniform($(x[]))")

## Wrap

struct Wrap{T,E} <: Distribution{E}
    val::T
end

Wrap(x::T) where {T} = Wrap{T,gentype(T)}(x)

Base.getindex(x::Wrap) = x.val

## wrap

wrap(::Type{T}) where {T} = Uniform(T)

wrap(x::T) where {T} = wrap(x, isuniform(T))

wrap(x::T, ::Val{true}) where {T} = Uniform(x)
wrap(x::T, ::Val{false}) where {T} = Wrap(x)

wrap(x::Distribution) = x
wrap(x::Sampler) = x

## sampling

Sampler(RNG::Type{<:AbstractRNG},
        d::Union{UniformWrap,UniformType,Wrap}, n::Repetition) =
            Sampler(RNG, d[], n)

## support / pmf

support(d::Union{Wrap,UniformWrap}) = support(d.val)
pmf(d::Union{Wrap,UniformWrap}, x) = pmf(d.val, x)


text:
using DrWatson
@quickactivate "Medical segmentation evaluation"
2+2
datadir()


text:
# Parameters:
#	N::Int
#	- the problem size
#	P::Array{Array{Array{Int, 1}, 1}, 1}
#	- output1 from partitions()
#	WI::Array{Int, 2}
#	- output2 from partitions()
# Return Values:
#	PT::Array{Array{Array{Int, 1}, 1}, 1} (Partition Tree)
#	- for each value, i, in PT[n][p], P[n][p] decomposes into P[n-1][i]
#	- length(PT[1]) = 0
function partition_tree(N::Int, P::Array{Array{Array{Int, 1}, 1}, 1}, WI::Array{Int, 2})
	PT = Array(Array{Array{Int, 1}, 1}, N)
	PT[1] = Array(Array{Int, 1}, 0)
	for n = 2:N
		Pd = P[n - 1]
		Pn = P[n]
		PTn = Array(Array{Int, 1} , length(Pn))
		for p = 1:length(Pn)
			PDA = pda(Pn[p])
			if length(PDA) == 1 || PDA[1][1] == PDA[2][1]
				lb = 1
				if PDA[1][1] != 1
					lb = WI[n - 1, PDA[1][1] - 1] + 1
				end
				PTn[p] = dia(PDA, Pd, lb)
			else
				lb1 = 1
				if PDA[1][1] != 1
					lb1 = WI[n - 1, PDA[1][1] - 1] + 1
				end
				lb2 = WI[n - 1, PDA[2][1] - 1]
				PTn[p] = dia(PDA, Pd, lb1, lb2)
			end
		end
		PT[n] = PTn
	end
	return PT
end

# Parameters:
#	P::Array{Int, 1}
#	- P is a Partition
# Return Values:
#	PDA::Array{Array{Int, 1}, 1} (Partition Decomposition Array)
#	- PDA[i] is the the ith Partition that P decomposes into
function pda(P::Array{Int, 1})
	P_L = length(P)
	num = 1
	for i = 1:(P_L - 1)
		if P[i] > P[i + 1]
			num += 1
		end
	end
	PDA = Array(Array{Int, 1}, num)
	c = 1
	for i = 1:(P_L - 1)
		if P[i] > P[i + 1]
			D = copy(P)
			D[i] -= 1
			PDA[c] = D
			c += 1
		end
	end
	if P[P_L] == 1
		D = copy(P[1:(P_L - 1)])
		PDA[c] = D
	else
		D = copy(P)
		D[P_L] -= 1
		PDA[c] = D
	end
	return PDA
end

# Parameters:
#	PDA::Array{Array{Int, 1}, 1}
#	- output1 from pda()
#	Pd::Array{Array{Int, 1}, 1}
#	- Pd is array containing all the partitions of the same size as those in PDA
#	lb::Int
#	- the lower bound for the index where the first element of PDA will be found
# Return Values:
#	DIA::Array{Int, 1} (Decomposition Index Array)
#	- PDA[i] = Pd[DIA[i]]
function dia(PDA::Array{Array{Int, 1}, 1}, Pd::Array{Array{Int, 1}, 1}, lb::Int)
	PDA_L = length(PDA)
	DIA = Array(Int, PDA_L)
	c = 1
	p = lb
	while true
		if Pd[p] == PDA[c]
			DIA[c] = p
			c += 1
			if c > PDA_L
				return DIA
			end
		end
		p += 1
	end
end

# Parameters:
#	PDA::Array{Array{Int, 1}, 1}
#	- output1 from pda()
#	Pd::Array{Array{Int, 1}, 1}
#	- Pd is array containing all the partitions of the same size as those in PDA
#	lb1::Int
#	- the lower bound for the index where the first element of PDA will be found
#	lb2::Int
#	- the lower bound for the index where the second element of PDA will be found
# Return Values:
#	DIA::Array{Int, 1}
#	- PDA[i] = Pd[DIA[i]]
function dia(PDA::Array{Array{Int, 1}, 1}, Pd::Array{Array{Int, 1}, 1}, lb1::Int, lb2::Int)
	PDA_L = length(PDA)
	DIA = Array(Int, PDA_L)
	p = lb1
	while true
		if Pd[p] == PDA[1]
 			DIA[1] = p
			break
		end
		p += 1
	end
	p = lb2
	c = 2
	while true
		if Pd[p] == PDA[c]
			DIA[c] = p
			c += 1
			if c > PDA_L
				return DIA
			end
		end
		p += 1
	end
end


text:
@testset "Mutual Information" begin
    @test_throws MethodError MutualInfo(0.5, 1.0)
    @test_throws MethodError MutualInfo(tuple(), 3)
    @test_throws MethodError MutualInfo(3, tuple())
    @test_throws MethodError MutualInfo(tuple(), tuple())
    @test_throws MethodError MutualInfo(3, (1.3,))
    for b in -2:1
        @test_throws ArgumentError MutualInfo(b, 2)
        @test_throws ArgumentError MutualInfo(2, b)
    end

    let mi = MutualInfo(2, 3)
        @test mi.joint == zeros(Int, 2, 3)
        @test mi.m1 == zeros(Int, 2)
        @test mi.m2 == zeros(Int, 3)
        @test mi.b1 == (2,)
        @test mi.b2 == (3,)
        @test mi.N == zero(Int)
    end

    let mi = MutualInfo((2, 2), 3)
        @test mi.joint == zeros(Int, 4, 3)
        @test mi.m1 == zeros(Int, 4)
        @test mi.m2 == zeros(Int, 3)
        @test mi.b1 == (2, 2)
        @test mi.b2 == (3,)
        @test mi.N == 0
    end

    let mi = MutualInfo(3, (2, 2))
        @test mi.joint == zeros(Int, 3, 4)
        @test mi.m1 == zeros(Int, 3)
        @test mi.m2 == zeros(Int, 4)
        @test mi.b1 == (3,)
        @test mi.b2 == (2, 2)
        @test mi.N == 0
    end

    let mi = MutualInfo(2, 2)
        observe!(mi, [2,1,2,1], [1,1,2,2])
        @test mi.joint == [1 1; 1 1]
        @test mi.m1 == [2, 2]
        @test mi.m2 == [2, 2]
        @test mi.b1 == (2,)
        @test mi.b2 == (2,)
        @test mi.N == 4
    end

    let mi = MutualInfo(3, 2)
        observe!(mi, [1,2,3,1,2,3], [1,1,2,2,1,2])
        @test mi.joint == [1 1; 2 0; 0 2]
        @test mi.m1 == [2, 2, 2]
        @test mi.m2 == [3, 3]
        @test mi.b1 == (3,)
        @test mi.b2 == (2,)
        @test mi.N == 6
    end

    let mi = MutualInfo(2, 2)
        observe!(mi, [1,1,1,1], [2,2,2,2])
        observe!(mi, [1,2,2,2], [1,1,2,2])
        @test mi.joint == [1 4; 1 2]
        @test mi.m1 == [5, 3]
        @test mi.m2 == [2, 6]
        @test mi.N == 8
    end

    let mi = MutualInfo(2, 2)

        xs = zeros(Int, 1, 4, 2)
        xs[:,:,1] = [1,1,1,1]
        xs[:,:,2] = [1,2,2,2]

        ys = zeros(Int, 1, 4, 2)
        ys[:,:,1] = [2,2,2,2]
        ys[:,:,2] = [1,1,2,2]

        observe!(mi, xs, ys)

        @test mi.joint == [1 4; 1 2]
        @test mi.m1 == [5, 3]
        @test mi.m2 == [2, 6]
        @test mi.N == 8
    end

    let mi = MutualInfo((2,2), (2,2))
        observe!(mi, [1 1 2 2; 1 2 1 2], [2 1 2 1; 1 1 1 1])
        @test mi.joint == [0 1 0 0; 0 1 0 0; 1 0 0 0; 1 0 0 0]
        @test mi.m1 == [1,1,1,1]
        @test mi.m2 == [2,2,0,0]
        @test mi.b1 == (2,2)
        @test mi.b2 == (2,2)
        @test mi.N == 4
    end

    let mi = MutualInfo((2,2), (2,2))
        observe!(mi, [1 2 2 1; 1 1 2 2], [1 1 2 1; 2 1 1 2])
        @test mi.joint == [0 0 1 0; 1 0 0 0; 0 0 1 0; 0 1 0 0]
        @test mi.m1 == [1,1,1,1]
        @test mi.m2 == [1,1,2,0]
        @test mi.b1 == (2,2)
        @test mi.b2 == (2,2)
        @test mi.N == 4
    end

    let mi = MutualInfo((2,2), (2,2))
        xs = zeros(Int, 2, 4, 2)
        xs[:,:,1] = [1 1 2 2; 1 2 1 2]
        xs[:,:,2] = [1 2 2 1; 1 1 2 2]

        ys = zeros(Int, 2, 4, 2)
        ys[:,:,1] = [2 1 2 1; 1 1 1 1]
        ys[:,:,2] = [1 1 2 1; 2 1 1 2]

        observe!(mi, xs, ys)

        @test mi.joint == [0 1 1 0; 1 1 0 0; 1 0 1 0; 1 1 0 0]
        @test mi.m1 == [2,2,2,2]
        @test mi.m2 == [3,3,2,0]
        @test mi.b1 == (2,2)
        @test mi.b2 == (2,2)
        @test mi.N == 8
    end

    @test_throws ArgumentError MutualInfo([0,1,1], [1,1,2])
    @test_throws ArgumentError MutualInfo([2,1,1], [0,1,2])
    @test_throws ArgumentError MutualInfo(Int[], Int[])
    @test_throws ArgumentError MutualInfo(Int[], [1,2])
    @test_throws ArgumentError MutualInfo([1,2], Int[])
    @test_throws ArgumentError MutualInfo([1,2,2], [1,2])
    @test_throws ArgumentError MutualInfo([1,2], [1,2,2])
    @test_throws ArgumentError MutualInfo(ones(Int, 3, 5), ones(Int, 3, 2))
    @test_throws ArgumentError MutualInfo(ones(Int, 3, 2), ones(Int, 3, 5))
    @test_throws ArgumentError MutualInfo(ones(Int, 3, 5, 7), ones(Int, 3, 5, 6))
    @test_throws ArgumentError MutualInfo(ones(Int, 3, 5, 7), ones(Int, 3, 4, 7))

    let mi = MutualInfo([1,1,1,1,1,2,2,2], [2,2,2,2,1,1,2,2])
        @test mi.joint == [1 4; 1 2]
        @test mi.m1 == [5, 3]
        @test mi.m2 == [2, 6]
        @test mi.b1 == (2,)
        @test mi.b2 == (2,)
        @test mi.N == 8
    end

    let mi = MutualInfo([1,1,1,1,1,1,1,1], [2,2,2,2,1,1,2,2])
        @test mi.joint == [2 6; 0 0]
        @test mi.m1 == [8, 0]
        @test mi.m2 == [2, 6]
        @test mi.b1 == (2,)
        @test mi.b2 == (2,)
        @test mi.N == 8
    end

    let xs = zeros(Int, 1, 4, 2)
        xs[:,:,1] = [1,1,1,1]
        xs[:,:,2] = [1,2,2,2]

        ys = zeros(Int, 1, 4, 2)
        ys[:,:,1] = [2,2,2,2]
        ys[:,:,2] = [1,1,2,2]

        mi = MutualInfo(xs, ys)

        @test mi.joint == [1 4; 1 2]
        @test mi.m1 == [5, 3]
        @test mi.m2 == [2, 6]
        @test mi.N == 8
    end

    let xs = [1 1 2 2; 1 2 1 2]
        ys = [2 1 2 1; 1 1 1 1]

        mi = MutualInfo(xs, ys)

        @test mi.joint == [0 1 0 0; 0 1 0 0; 1 0 0 0; 1 0 0 0]
        @test mi.m1 == [1,1,1,1]
        @test mi.m2 == [2,2,0,0]
        @test mi.N == 4
    end

    let xs = [1 2 2 1; 1 1 2 2]
        ys = [1 1 2 1; 2 1 1 2]

        mi = MutualInfo(xs, ys)

        @test mi.joint == [0 0 1 0; 1 0 0 0; 0 0 1 0; 0 1 0 0]
        @test mi.m1 == [1,1,1,1]
        @test mi.m2 == [1,1,2,0]
        @test mi.N == 4
    end

    let xs = zeros(Int, 2, 4, 2)
        xs[:,:,1] = [1 1 2 2; 1 2 1 2]
        xs[:,:,2] = [1 2 2 1; 1 1 2 2]

        ys = zeros(Int, 2, 4, 2)
        ys[:,:,1] = [2 1 2 1; 1 1 1 1]
        ys[:,:,2] = [1 1 2 1; 2 1 1 2]

        mi = MutualInfo(xs, ys)

        @test mi.joint == [0 1 1 0; 1 1 0 0; 1 0 1 0; 1 1 0 0]
        @test mi.m1 == [2,2,2,2]
        @test mi.m2 == [3,3,2,0]
        @test mi.N == 8
    end

    @test_throws ArgumentError observe!(MutualInfo(2,2), [1,2,3], [1,2])
    @test_throws ArgumentError observe!(MutualInfo(2, 2), [0,1,0,1], [0,0,1,1])
    @test_throws BoundsError observe!(MutualInfo(2, 2), [2,3,2,3], [2,2,3,3])

    let mi = observe!(MutualInfo(2, 2), [1,1,2,2], [1,2,1,2])
        clear!(mi)
        @test mi.joint == zeros(Int, 2, 2)
        @test mi.m1 == zeros(Int, 2)
        @test mi.m2 == zeros(Int, 2)
        @test mi.b1 == (2,)
        @test mi.b2 == (2,)
        @test mi.N == zero(Int)
    end

    let mi = observe!(MutualInfo((2,2), (2,3)),
                      [1 1 2 2 1 1; 1 2 1 2 2 1],
                      [1 1 1 1 2 2; 2 3 2 1 1 3])
        clear!(mi)
        @test mi.joint == zeros(Int, 4, 6)
        @test mi.m1 == zeros(Int, 4)
        @test mi.m2 == zeros(Int, 6)
        @test mi.b1 == (2,2)
        @test mi.b2 == (2,3)
        @test mi.N == 0
    end

    let xs = zeros(Int, 2, 4, 2)
        xs[:,:,1] = [1 1 2 2; 1 2 1 2]
        xs[:,:,2] = [1 2 2 1; 1 1 2 2]

        ys = zeros(Int, 2, 4, 2)
        ys[:,:,1] = [2 1 2 1; 1 1 1 1]
        ys[:,:,2] = [1 1 2 1; 2 1 1 2]

        mi = observe!(MutualInfo((2,2), (2,2)), xs, ys)
        clear!(mi)

        @test mi.joint == zeros(Int, 4, 4)
        @test mi.m1 == zeros(Int, 4)
        @test mi.m2 == zeros(Int, 4)
        @test mi.N == 0
    end

    let as = [1,1,1,1,2,2,2,2]
        bs = [2,2,2,2,1,1,1,1]
        @test mutualinfo(as, bs) ≈ 1.0 atol=1e-6
        @test mutualinfo(bs, as) ≈ 1.0 atol=1e-6
    end

    let as = [1,1,2,2,2,2,1,1,1]
        bs = [2,2,1,1,1,1,2,2,2]
        @test mutualinfo(as, bs) ≈ 0.991076 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.991076 atol=1e-6
    end

    let as = [2,2,1,2,1,2,2,2,1]
        bs = [2,2,1,1,1,2,1,2,2]
        @test mutualinfo(as, bs) ≈ 0.072780 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.072780 atol=1e-6
    end

    let as = [1,1,1,1,1,1,1,1,1]
        bs = [2,2,2,1,1,1,2,2,2]
        @test mutualinfo(as, bs) ≈ 0.0 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.0 atol=1e-6
    end

    let as = [2,2,2,2,1,1,1,1,2]
        bs = [2,2,2,1,1,1,2,2,2]
        @test mutualinfo(as, bs) ≈ 0.072780 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.072780 atol=1e-6
    end

    let as = [2,2,1,1,2,2,1,1,2]
        bs = [2,2,2,1,1,1,2,2,2]
        @test mutualinfo(as, bs) ≈ 0.018311 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.018311 atol=1e-6
    end

    let as = [2 2 2 2 1 1 1 1 2;
              2 2 1 1 2 2 1 1 2]
        bs = [2 2 2 1 1 1 2 2 2;
              2 2 2 1 1 1 2 2 2]
        mi = MutualInfo(2, 2)
        @test mutualinfo!(mi, as[1,:], bs[1,:]) ≈ 0.072780 atol=1e-6
        @test mutualinfo!(mi, as[2,:], bs[2,:]) ≈ 0.004497 atol=1e-6
    end

    let as = zeros(Int, 1, 9, 2)
        as[:,:,1] = [2,2,2,2,1,1,1,1,2]
        as[:,:,2] = [2,2,1,1,2,2,1,1,2]

        bs = zeros(Int, 1, 9, 2)
        bs[:,:,1] = [2,2,2,1,1,1,2,2,2]
        bs[:,:,2] = [2,2,2,1,1,1,2,2,2]

        @test mutualinfo(as, bs) ≈ 0.004497 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.004497 atol=1e-6
    end

    let as = [2 2 2 2 1 1 1 1 2;
              2 2 1 1 2 2 1 1 2]
        bs = [2 2 2 1 1 1 2 2 2;
              2 2 2 1 1 1 2 2 2]

        @test mutualinfo(as, bs) ≈ 0.696074 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.696074 atol=1e-6
    end

    let as = [2 1 2 2 1 2 2 1 2;
              2 2 1 1 2 2 1 1 2]
        bs = [2 2 2 1 1 1 2 2 2;
              2 1 2 1 2 1 2 1 2]

        @test mutualinfo(as, bs) ≈ 1.002172 atol=1e-6
        @test mutualinfo(bs, as) ≈ 1.002172 atol=1e-6
    end

    let as = zeros(Int, 2, 9, 2)
        as[:,:,1] = [2 1 2 2 1 2 2 1 2;
                     2 2 1 1 2 2 1 1 2]
        as[:,:,2] = [2 2 2 2 1 1 1 1 2;
                     2 2 1 1 2 2 1 1 2]

        bs = zeros(Int, 2, 9, 2)
        bs[:,:,1] = [2 2 2 1 1 1 2 2 2;
                     2 1 2 1 2 1 2 1 2]
        bs[:,:,2] = [2 2 2 1 1 1 2 2 2;
                     2 2 2 1 1 1 2 2 2]

        @test mutualinfo(as, bs) ≈ 0.595553 atol=1e-6
        @test mutualinfo(bs, as) ≈ 0.595553 atol=1e-6
    end

    @testset "Lagged Mutual Information" begin
        let xs = [1,2,1,1,2,2]
            ys = [3,1,2,1,1,2]
            @test mutualinfo(xs, ys; lag=1) == mutualinfo(xs[1:end-1], ys[2:end])
            @test mutualinfo(xs, ys; lag=2) == mutualinfo(xs[1:end-2], ys[3:end])
            @test mutualinfo(xs, ys; lag=-1) == mutualinfo(ys[1:end-1], xs[2:end])
            @test mutualinfo(xs, ys; lag=-2) == mutualinfo(ys[1:end-2], xs[3:end])
        end

        let xs = [2 1 2 2 1 2 2 1 2;
                  2 2 1 1 2 2 1 1 2]
            ys = [2 2 2 1 1 1 2 2 2;
                  2 1 2 1 2 1 2 1 2]

            @test mutualinfo(xs, ys; lag=1) == mutualinfo(xs[:,1:end-1], ys[:,2:end])
            @test mutualinfo(xs, ys; lag=2) == mutualinfo(xs[:,1:end-2], ys[:,3:end])
            @test mutualinfo(xs, ys; lag=-1) == mutualinfo(ys[:,1:end-1], xs[:,2:end])
            @test mutualinfo(xs, ys; lag=-2) == mutualinfo(ys[:,1:end-2], xs[:,3:end])
        end

        let xs = zeros(Int, 2, 9, 2)
            xs[:,:,1] = [2 1 2 2 1 2 2 1 2;
                         2 2 1 1 2 2 1 1 2]
            xs[:,:,2] = [2 2 2 2 1 1 1 1 2;
                         2 2 1 1 2 2 1 1 2]

            ys = zeros(Int, 2, 9, 2)
            ys[:,:,1] = [2 2 2 1 1 1 2 2 2;
                         2 1 2 1 2 1 2 1 2]
            ys[:,:,2] = [2 2 2 1 1 1 2 2 2;
                         2 2 2 1 1 1 2 2 2]

            @test mutualinfo(xs, ys; lag=1) == mutualinfo(xs[:,1:end-1,:], ys[:,2:end,:])
            @test mutualinfo(xs, ys; lag=2) == mutualinfo(xs[:,1:end-2,:], ys[:,3:end,:])
            @test mutualinfo(xs, ys; lag=-1) == mutualinfo(ys[:,1:end-1,:], xs[:,2:end,:])
            @test mutualinfo(xs, ys; lag=-2) == mutualinfo(ys[:,1:end-2,:], xs[:,3:end,:])
        end
    end
end


text:
## https://developer.spotify.com/documentation/web-api/reference/#/operations/get-available-markets

"""
    markets_get()

**Summary**: Get the list of markets where Spotify is available.

# Example
```julia-repl
julia> Spotify.markets_get()[1]
JSON3.Object{Base.CodeUnits{UInt8, String}, Vector{UInt64}} with 1 entry:
  :markets => ["AD", "AE", "AG", "AL", "AM", "AO", "AR", "AT", "AU", "AZ"  …
```
"""
function markets_get()

    return Spotify.spotify_request("markets")

end


text:
# ellipse to multiple forward simulated shapes
# only final configurations observed

using BridgeLandmarks
using Random
using Distributions
using StaticArrays
using LinearAlgebra
using JLD2

workdir = @__DIR__
cd(workdir)

Random.seed!(9)

σobs = 0.01

n = 15
nshapes = 10

a0 = 2.0
c0 = 0.1
γ0 = 0.7

q0 = [PointF(2.0cos(t), sin(t))  for t in collect(0:(2pi/n):2pi)[2:end]]
xobsT = Vector{PointF}[]
Ptrue = MarslandShardlow(a0, c0, γ0, 0.0, n)

for k in 1:nshapes
    # Ptrue = MarslandShardlow(a0 * exp(0.5*randn()), c0 * exp(0.5*randn()),
    #                 γ0 *  exp(0.5*randn()), 0.0, n)
    x0 = State(q0, randn(PointF,n))
    T = 1.0; dt = 0.01; t = 0.0:dt:T
    Wf, Xf = landmarksforward(t, x0, Ptrue)
    push!(xobsT, [Xf.yy[end].q[i] for i in 1:n ] + σobs * randn(PointF,n))
end
xobs0 = []

JLD2.@save "data_exp2.jld2" xobs0 xobsT n nshapes q0


text:
using LinearAlgebra

function to_rpr_camera(context::RPR.Context, cam_controls, cam)
    # Create camera
    camera = RPR.Camera(context)

    map(cam_controls.eyeposition) do position
        l, u = cam_controls.lookat[], cam_controls.upvector[]
        RPR.rprCameraSetFocusDistance(camera, norm(l - position))
        lookat!(camera, position, l, u)
        return
    end

    map(cam_controls.far) do far
        RPR.rprCameraSetFarPlane(camera, far)
        return
    end

    map(cam_controls.near) do near
        RPR.rprCameraSetNearPlane(camera, near)
        return
    end

    map(cam_controls.fov, cam.resolution) do fov, res
        l, p = cam_controls.lookat[], cam_controls.eyeposition[]
        wd = norm(l-p)
        h = norm(res)
        return RPR.rprCameraSetFocalLength(camera, (h*wd)/fov)
    end
    # TODO:
    # RPR_CAMERA_FSTOP
    # RPR_CAMERA_FOCAL_LENGTH
    # RPR_CAMERA_SENSOR_SIZE
    # RPR_CAMERA_MODE
    # RPR_CAMERA_FOCUS_DISTANCE
    return camera
end

function to_rpr_object(context, matsys, scene, plot)
    @warn("Not supported $(typeof(plot))")
    return nothing
end

function insert_plots!(context, matsys, scene, mscene::Makie.Scene, @nospecialize(plot::Combined))
    if isempty(plot.plots) # if no plots inserted, this truly is an atomic
        object = to_rpr_object(context, matsys, mscene, plot)
        if !isnothing(object)
            if object isa AbstractVector
                foreach(x-> push!(scene, x), object)
            else
                push!(scene, object)
            end
        end
    else
        for plot in plot.plots
            insert_plots!(context, matsys, scene, mscene, plot)
        end
    end
end

function to_rpr_scene(context::RPR.Context, matsys, mscene::Makie.Scene)
    scene = RPR.Scene(context)
    set!(context, scene)
    cam_controls = Makie.cameracontrols(mscene)
    cam = Makie.camera(mscene)
    camera = to_rpr_camera(context, cam_controls, cam)
    set!(scene, camera)

    env_light = RPR.EnvironmentLight(context)
    image_path = RPR.assetpath("studio032.exr")
    # env_img = fill(to_color(mscene.backgroundcolor[]), 1, 1)
    img = RPR.Image(context, image_path)
    set!(env_light, img)
    setintensityscale!(env_light, 0.5)
    push!(scene, env_light)

    light = RPR.PointLight(context)
    transform!(light, Makie.translationmatrix(Vec3f0(10, 10, 10)))
    RPR.setradiantpower!(light, 1000, 1000, 1000)
    push!(scene, light)

    for plot in mscene.plots
        insert_plots!(context, matsys, scene, mscene, plot)
    end
    return scene, camera
end

function replace_scene_rpr!(scene,
        context=RPR.Context(resource=RPR.RPR_CREATION_FLAGS_ENABLE_GPU0),
        matsys = RPR.MaterialSystem(context, 0); refresh=Observable(nothing))
    set_standard_tonemapping!(context)
    rpr_scene, rpr_camera = RPRMakie.to_rpr_scene(context, matsys, scene)
    # hide Makie scene
    scene.visible[] = false
    # foreach(p-> delete!(scene, p), copy(scene.plots))
    sub = campixel(scene)
    fb_size = size(scene)
    im = image!(sub, zeros(RGBAf, fb_size))
    framebuffer1 = RPR.FrameBuffer(context, RGBA, fb_size)
    framebuffer2 = RPR.FrameBuffer(context, RGBA, fb_size)
    RPR.rprCameraSetSensorSize(rpr_camera, fb_size...)
    onany(camera(scene).projection, refresh) do proj, refresh
        clear!(framebuffer1)
    end
    set!(context, RPR.RPR_AOV_COLOR, framebuffer1)
    RPR.rprContextSetParameterByKey1u(context, RPR.RPR_CONTEXT_ITERATIONS, 1)
    task = @async while isopen(scene)
        if fb_size != size(scene)
            @info("resizing scene")
            fb_size = size(scene)
            RPR.release(framebuffer1)
            RPR.release(framebuffer2)
            framebuffer1 = RPR.FrameBuffer(context, RGBA, fb_size)
            framebuffer2 = RPR.FrameBuffer(context, RGBA, fb_size)
            set!(context, RPR.RPR_AOV_COLOR, framebuffer1)
            RPR.rprCameraSetSensorSize(rpr_camera, fb_size...)
        end
        RPR.render(context)
        RPR.rprContextResolveFrameBuffer(context, framebuffer1, framebuffer2, false)
        data = RPR.get_data(framebuffer2)
        im[1] = reverse(reshape(data, fb_size), dims=2)
        sleep(0.01)
    end
    return context, task
end


text:
"""
Operation{rank} <: AbstractTensor{rank}

A supertype for all tensor operations.
"""
abstract type Operation{rank} <: AbstractTensor{rank} end

"""
    IndexSumOperation(A::AbstractTensor, indices::Indices, freeindices::Vararg{FreeIndex})

Symbolises summing over given free index of a tensor.
"""
struct IndexSumOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function IndexSumOperation(A::AbstractTensor, indices::Indices, freeindices::Vararg{FreeIndex})
        new{length(A.shape)}(A.shape, (A, indices), freeindices)
    end
end

"""
    contractioncheck(A::AbstractTensor)
Checks if we have have an upper and lower index. If there is, we create
an [`IndexSumOperation`](@ref) node.
"""
function contractioncheck(A::AbstractTensor)
    contractionindices = []
    for (i, index) in enumerate(A.freeindices)
        if index' in A.freeindices[(i+1):end]
            push!(contractionindices, index)
        end
    end
    if isempty(contractionindices)
        return A
    end
    freeindices = tuple(setdiff(A.freeindices, [i for i in contractionindices],
        [i' for i in contractionindices])...)
    return IndexSumOperation(A, Indices(contractionindices...), freeindices...)
end

"""
    AddOperation(children::Tuple{Vararg{AbstractTensor{rank}}},
        freeindices::Tuple{Vararg{FreeIndex}}) where rank

A node symbolising addition of tensors.
"""
struct AddOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{Vararg{AbstractTensor}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function AddOperation(children::Tuple{Vararg{AbstractTensor{rank}}},
            freeindices::Tuple{Vararg{FreeIndex}}) where rank
        if length(children) > 1
            for child in children[2:length(children)]
                if child.shape != children[1].shape
                    throw(DimensionMismatch(string(
                        "Shapes ", children[1].shape, " and ", child.shape,
                        " don't match")))
                end
            end
        end
        newchildren = tuple(filter!(x -> !(x isa ZeroTensor),
            [children...])...)
        if length(newchildren) == 1
            return newchildren[1]
        end
        if length(newchildren) == 0
            return children[1]
        end
        return contractioncheck(new{rank}(children[1].shape, newchildren,
            freeindices))
    end
end

"""
    Base.:+(nodes::Vararg{Node})

Creates an [`AddOperation`](@ref) whose children are the nodes. The shapes of
the nodes must match.
"""
function Base.:+(nodes::Vararg{Node})
    return AddOperation(nodes, tuple(union([node.freeindices for node=nodes]...)...))
end

"""
    OuterProductOperation(shape::Tuple{Vararg{AbstractVectorSpace}},
        children::Tuple{AbstractTensor, AbstractTensor},

A node symbolising tensor outer product."""
struct OuterProductOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, AbstractTensor}
    freeindices::Tuple{Vararg{FreeIndex}}
    function OuterProductOperation(shape::Tuple{Vararg{AbstractVectorSpace}},
            children::Tuple{AbstractTensor, AbstractTensor},
            freeindices::Tuple{Vararg{FreeIndex}})
        return contractioncheck(new{length(shape)}(shape, children,
            freeindices))
    end
end

"""
    ⊗(A::AbstractTensor, B::AbstractTensor)

Returns an [`OuterProductOperation`](@ref) with A and B as its children.
"""
function ⊗(A::AbstractTensor, B::AbstractTensor)
    return OuterProductOperation((A.shape..., B.shape...), (A, B),  (A.freeindices..., B.freeindices...))
end

"""
    Base.:*(x::Scalar, A::AbstractTensor)

Shorthand for multiplying a tensor by a scalar.
If multiplying a number by a variable, the number must be first.
"""
function Base.:*(x::Scalar, A::AbstractTensor)
    if x == 1 || x == ConstantTensor(1)
        return A
    elseif A == ConstantTensor(1)
        return Tensor(x)
    elseif x == 0 || x isa ZeroTensor || A isa ZeroTensor
        return ZeroTensor(A.shape...)
    else
        return Tensor(x) ⊗ A
    end
end

"""
    Base.:-(A::AbstractTensor)

Unary minus. Creates an [`OuterProductOperation`](@ref) between ConstantTensor(-1) and A.
"""
function Base.:-(A::AbstractTensor)
    return -1*A
end

"""
    Base.:-(A::AbstractTensor{rank}, B::AbstractTensor{rank}) where rank

Binary minus. Creates an [`AddOperation`](@ref) of A and -1*B.
"""
function Base.:-(A::AbstractTensor{rank}, B::AbstractTensor{rank}) where rank
    return A + (-1*B)
end

"""
    PowerOperation(x::AbstractTensor{0}, y::AbstractTensor{0})

Symbolises raising a scalar-shaped tensor x to a scalar-shaped power y.
"""
struct PowerOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor{0}, AbstractTensor{0}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function PowerOperation(x::AbstractTensor{0}, y::AbstractTensor{0})
        new{0}((), (x, y), (x.freeindices..., y.freeindices...))
    end
end

"""
    Base.:^(x::AbstractTensor{0}, y::AbstractTensor{0})

Raises the scalar-shaped tensor x to the scalar-shaped power y. Returns a
[`PowerOperation`](@ref).
"""
function Base.:^(x::AbstractTensor{0}, y::AbstractTensor{0})
    return PowerOperation(x, y)
end

"""
    Base.:^(x::AbstractTensor{0}, y::Number)

Raises the scalar-shaped tensor x to the power y. Returns a
[`PowerOperation`](@ref).
"""
function Base.:^(x::AbstractTensor{0}, y::Number)
    return x^Tensor(y)
end

"""
    Base.sqrt(x::AbstractTensor{0})

Returns the square root of x as a [`PowerOperation`](@ref).
"""
function Base.sqrt(x::AbstractTensor{0})
    return PowerOperation(x, Tensor(1//2))
end

"""
    DivisionOperation(A::AbstractTensor, y::AbstractTensor{0})

Symbolises dividing tensor A by a scalar-shaped tensor y.
"""
struct DivisionOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, AbstractTensor{0}}
    freeindices::Tuple{Vararg{FreeIndex}}
    function DivisionOperation(A::AbstractTensor, y::AbstractTensor{0})
        if y isa ZeroTensor
            throw(DivideError())
        elseif y == ConstantTensor(1)
            return A
        end
        return contractioncheck(new{length(A.shape)}((A.shape), (A, y), (A.freeindices..., y.freeindices...)))
    end
end

"""
    Base.:/(A::AbstractTensor, y::AbstractTensor{0})

Creates a DivisionOperation where A is divived by y. y must be scalar-shaped.
"""
function Base.:/(A::AbstractTensor, y::AbstractTensor{0})
    return DivisionOperation(A, y)
end

"""
    Base.:/(A::AbstractTensor, y::Number)

Creates a DivisionOperation where A is divived by y.
"""
function Base.:/(A::AbstractTensor, y::Number)
    return DivisionOperation(A, Tensor(y))
end

"""
    Base.:/(y::Number, A::AbstractTensor{0})

Creates a DivisionOperation where y is divided by A, where A must be a scalar.
"""
function Base.:/(y::Number, A::AbstractTensor{0})
    return DivisionOperation(Tensor(y), A)
end

"""
    IndexingOperation(x::AbstractTensor, indices::Indices)

A node symbolising indexing a tensor by its every dimension.
"""
struct IndexingOperation{rank} <: Operation{rank}
    shape::Tuple{}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function IndexingOperation(x::AbstractTensor, indices::Indices)
        return contractioncheck(new{0}((),(x, indices), tuple(x.freeindices..., [i for i=indices.indices if i isa FreeIndex]...)))
    end
end

"""
    ComponentTensorOperation(shape::Tuple{Vararg{AbstractVectorSpace}},
        children::Tuple{AbstractTensor, Indices},
        freeindices::Tuple{Vararg{FreeIndex}})

A node symbolising component tensor.
"""
struct ComponentTensorOperation{rank} <: Operation{rank}
    shape::Tuple{Vararg{AbstractVectorSpace}}
    children::Tuple{AbstractTensor, Indices}
    freeindices::Tuple{Vararg{FreeIndex}}
    function ComponentTensorOperation(
            shape::Tuple{Vararg{AbstractVectorSpace}},
            children::Tuple{AbstractTensor, Indices},
            freeindices::Tuple{Vararg{FreeIndex}})
        new{length(shape)}(shape, children, freeindices)
    end
end

"""
    componenttensor(A::AbstractTensor, indices::Vararg{Index})

Creates a [`ComponentTensorOperation`](@ref) of A over indices. Indices must be a
subset of the free indices of A.
"""
function componenttensor(A::AbstractTensor, indices::Vararg{Index})
    if length(indices) == 0
        return A
    end
    if !(indices ⊆ A.freeindices)
        throw(DomainError(string("The free indices ", indices, " are not a subset of ", A.freeindices)))
    end
    shape = tuple(A.shape..., [i.V for i in indices]...)
    freeindices = tuple(setdiff(A.freeindices, indices, [i' for i in indices])...)
    return ComponentTensorOperation(shape, (A, Indices(indices...)), freeindices)
end

idcounter = 0

"""
    Base.getindex(A::AbstractTensor, ys::Vararg{Index})

Creates an [`IndexingOperation`](@ref) of A indexed by ys. If every dimension of A is
not indexed, creates a ComponentTensorOperation over the unindexed dimensions."""
function Base.getindex(A::AbstractTensor, ys::Vararg{Index})
    if length(ys) == 0
        return A
    end
    if length(A.shape) < length(ys)
        throw(BoundsError(A.shape, ys))
    end
    for i in 1:length(ys)
        if A.shape[i] != ys[i].V
            throw(DimensionMismatch(string(ys[i], " is not in the vector space ", A.shape[i])))
        end
    end
    addedindices = []
    global idcounter
    for i in 1:(length(A.shape) - length(ys))
        push!(addedindices, FreeIndex(A.shape[length(ys) + i], "", idcounter))
        idcounter += 1
    end
    return componenttensor(IndexingOperation(A, Indices(ys..., addedindices...)), addedindices...)
end

"""
    Base.getindex(A::AbstractTensor, ys::Vararg{Union{String, Int, Index}})

A convenience function that allows indexing a tensor by an integer. The
function creates a corresponding [`FixedIndex`](@ref) object in the appropiate
vector space.

Technically the function can also create a [`FreeIndex`](@ref) from a string,
but this is not recommended, because the syntax makes no distinction between
upper and lower indices.
"""
function Base.getindex(A::AbstractTensor, ys::Vararg{Union{String, Int, Index}})
    indexarray = []
    for i in 1:length(ys)
        push!(indexarray, toindex(A.shape[i], ys[i]))
    end
    return Base.getindex(A, indexarray...)
end

"""
    Base.transpose(A::AbstractTensor)

Reverses the shape of A. Does this by indexing A and then creating a
[`ComponentTensorOperation`](@ref) with the indices reversed.
"""
function Base.transpose(A::AbstractTensor)
    indices = []
    for V in A.shape
        global idcounter
        push!(indices, FreeIndex(V, "tranpose", idcounter))
        idcounter += 1
    end
    return componenttensor(A[indices...], reverse(indices)...)
end

"""
    trace(A::AbstractTensor{2})

Returns the trace of a matrix of the shape (V, V').
"""
function trace(A::AbstractTensor{2})
    if A.shape[1] != dual(A.shape[2])
        throw(DomainError(A, "Cannot get a trace from a tensor whose shape is not of the form (V, V')"))
    end
    global counter
    counter += 1
    return A[FreeIndex(A.shape[1], "trace", counter), FreeIndex(A.shape[2], "trace", counter)]
end

function Base.show(io::IO, op::Operation, depth::Int)
    println(io, ["    " for i in 1:depth]..., typeof(op))
    for child in op.children
        if child isa Operation
            Base.show(io, child, depth + 1)
        else
            println(io, ["    " for i in 1:(depth+1)]..., child)
        end
    end
end

Base.show(io::IO, op::Operation) = Base.show(io, op, 0)

max_stars_repo_name: Sohva/oolong

max_issues_repo_name: Sohva/oolong

max_forks_repo_name: oolongtensor/oolong


text:
export data
export children
export update!
# export insert!
export root

"""
    data(tree, node)

Retrieve the data aka payload associated with the given node.
"""
function data end

"""
    root(tree)

Return a proxy for the root of the tree.
"""
function root end

"""
    insert!(tree, parent, data)

Insert a node carrying 'data' as a new child of 'parent'
"""
function insert! end

"""
    The expression `children(tree,node)` returns an iterator that will produce
    a sequence of node iterators. These values do not have a lot of meaning by
    themselves, but can be used in conjunction with the tree object. E.g:

        data(tree, node_itr)
        children(tree, node_itr)

    In fact, the node iterators should be regarded as lightweight proxies for
    the underlying node and their attached data payload. The node objects
    themselves are of limited use for the client programmer as they are an
    implementation detail of the specific tree being used.
"""
function children end

function haschildren end

"""
Traverse the tree depth first, executing the function `f(tree, node, level)`
at every node. If `f` returns `false`, recursion is halted and the next node
on the current level is visited.
"""
function depthfirst(f, tree, node=root(tree), level=1)
    f(tree, node, level) || return
    for c in FMMTrees.children(tree, node)
        depthfirst(f, tree, c, level+1)
    end
end

function print_tree(tree, node=root(tree); maxdepth=0)
    depthfirst(tree) do tree, node, level
        print("-"^(level-1))
        println(FMMTrees.data(tree, node))
        return level == maxdepth ? false : true
    end
end

function route! end

"""
    update!(f, tree, node, data, router!)

Algorithm to update or add nodes of the tree. `router!` and `updater!` are
user supplied functions:

    router!(tree, node)

Returns the next candidate `node` until the node for insertion is reaches. Note
that this function potentially created new nodes. Arrival at the destination is
indicated by returning the same node that was passed as the second argument.

    f(tree, node, data)

Update the destination node `node`. Typically, `data` is added in some sense
to the data residing at the desitination node.
"""
function update!(f, tree, state, data, target)
    while true
        next_state = route!(tree, state, target)
        next_state == state && break
        state = next_state
    end
    node = state[1]
    f(tree, node, data)
    return node
end

struct DepthFirstIterator{T,N}
    tree::T
    node::N
end

Base.IteratorSize(::DepthFirstIterator) = Base.SizeUnknown()

function Base.iterate(itr::DepthFirstIterator)
    chitr = children(itr.tree,itr.node)
    stack = Any[(chitr, iterate(chitr))]
    iterate(itr, stack)
end

function Base.iterate(itr::DepthFirstIterator, stack)
    isempty(stack) && return nothing
    while true
        chditr, next = last(stack)
        if next != nothing
            node, state = next
            chitr = children(itr.tree, node)
            push!(stack, (chitr, iterate(chitr)))
        else
            pop!(stack)
            isempty(stack) && return (itr.node, stack)
            chitr, (node, state) = last(stack)
            stack[end] = (chitr, iterate(chitr, state))
            return node, stack
        end
    end
end

function leaves(tree, node=root(tree))
    Iterators.filter(n->!haschildren(tree,n), DepthFirstIterator(tree,node))
end


text:
function romanize(cur_int::Integer)
  @assert cur_int < 19

  if cur_int < 9
    ( cur_int < 4 ) && return "I" ^ cur_int
    ( cur_int == 4 ) && return "IV"

    cur_string = "V"
    cur_string *= "I" ^ ( cur_int - 5 )
    cur_string
  else
    ( cur_int == 9 ) && return "IX"

    cur_string = "X"
    ( cur_int < 14 ) && return cur_string * ( "I" ^ ( cur_int - 10 ) )
    ( cur_int == 14 ) && return "XIV"

    cur_string *= "V"
    cur_string *= "I" ^ ( cur_int - 5 - 10 )
    cur_string
  end
end


text:
using LinearAlgebra
"""
        (V,H,hist)=iar(prob,startvec,N)

### Infinite Arnoldi method
Infinite Arnoldi method, as described in Algorithm 2 in  "A linear eigenvalue algorithm for the nonlinear eigenvalue problem",
    by Jarlebring, Elias and Michiels, Wim and Meerbergen, Karl.

Runs `N` steps of the infinite Arnoldi method started with `startvec` with target `σ=0`.
The argument `prob` represents a problem `M(λ)` and has associated functions
`compute_Mder`, `compute_Mlincomb`:

* `M=compute_Mder(prob,s)` computes the matrix function value at `s=λ`.
* `z=compute_Mlincomb(prob,s,X)` computes a linear combination of derivatives. The linear combination coefficients are given in the matrix `X`:  `z=M(s)*X[:,1]+M'(s)*X[:,2]+...+M^{(k)}(s)*X[:,k+1]`.

"""
function iar(prob,startvec,N,precomputed_factorization=nothing)

    T=ComplexF64; # Problem type

    n=length(startvec);

    sigma=0;  # Only target implemented

    x0=zeros(T,n);
    x0[:]=normalize(startvec);
    H=zeros(T,N+1,N);
    print("IAR Iteration:");
    V=zeros(T,n*(N+1),N+1);
    V[1:n,1]=x0;

    start_time=time();
    hist=Dict(:time_count => NaN*zeros(N))

    if (isnothing(precomputed_factorization))
        Msolve=factorize(compute_Mder(prob,sigma));
    else
        Msolve=precomputed_factorization;
    end

    for k=1:N
        print("$k ")
        # Create L matrix: Equation (17) in nummath paper
        L=Diagonal(1 ./(1:k));

        vecx=V[1:(k*n),k];
        X=reshape(vecx,n,k);

        # Compute y1,..yn: Equation (13) in nummath paper
        Y=zeros(T,size(X,1),size(X,2)+1);
        Y[:,2:end]=X*L;

        # Compute y0: Equation (18) in nummath paper
        # y0=nep.M_lin_comb(X,Y,sigma);
        if (true)
            # Note: First col of Y is zero
            z0=compute_Mlincomb(prob,sigma,Y);
            #z0=Mlinprob.M_lin_comb(Y(:,2:end),1);
        else
            z0=compute_Mlincomb(prob,sigma,X[:,1:k]);
        end

        y0=-(Msolve\z0);
        Y[:,1]=y0;
        vecY=reshape(Y,n*(k+1),1);

        if (true)
            kk=k+1;
        else
            kk=1; # short vector orthogonalization
        end
        V0=V[1:(kk*n),1:k];
        h=V0'*vecY[1:(kk*n)];
        vecY=vecY-V[1:((k+1)*n),1:k]*h;

        g=V0'*vecY[1:(kk*n)];   # twice to be sure
        vecY=vecY-V[1:((k+1)*n),1:k]*g;
        h=h+g;

        beta=norm(vecY);  # normalize
        vecY=vecY/beta;
        H[1:(k+1),k]=[h;beta];
        V[1:((k+1)*n),k+1]=vecY;

        hist[:time_count][k]=time()-start_time;

    end

    println("done");

    return (V,H,hist);

end

max_stars_repo_name: siobhanie/InfGMRES

max_issues_repo_name: siobhanie/InfGMRES

max_forks_repo_name: siobhanie/InfGMRES


text:
@test super(BasicKalmanFilter) == LinearKalmanFilter

s = State([1.,0],[0.5 0; 0 1e-10])

@test typeof(s) == State{Float64}
@test s.x == [1,0]
@test s.p == [0.5 0; 0 1e-10]

o = convert(Observation,2.0)

@test typeof(o) == Observation{Float64}
@test o.y == [2]

x0 = State([1.,0],[0.5 0; 0 1e-10])
a = [1 0.01; 0.01*(-1) 1]
g = [0 0; 0 1]
q = [1e-10 0; 0 1e-10]
f = LinearModel(a,g,q)
h = [1 0]
r = [0.5]'
z = LinearObservationModel(h,r)

dt = 0.01
t = [0:dt:10]
x = zeros(length(t))
m = cos(t)

kf0 = BasicKalmanFilter(x0,f,z)
kf = kf0
x[1] = kf.x.x[1]

y = map(t->Observation(m[t]+kf.z.r[1]*randn(1)),1:length(t))

kf1 = predictupdate(kf,y[1])
predictupdate!(kf,y[1])
@test kf1.x == kf.x

for i in 2:length(t)-1
    kf2 = predict(kf)
    kf = update(kf2,y[i])
    x[i+1] = kf.x.x[1]
end

@test predictupdate(kf0,y[1]).x == update(predict(kf0),y[1]).x


text:
"""

## fsparv!

Function fsparv! assembles fin_el matrices into a symmetric skyline
global matrix. The Skyline vector kv is updated.

### Method
```julia
fsparv!(kv, km, g, km)
```

### Arguments
```julia
* kv::Vector{Float64}        : Sparse stiffnes matrix (Updated)
* km::Matrix{Float64}        : Symmetric element stiffnes matrix
* g::Vector{Int}           : Global steering vector.
* kdiag::Vector{Int}       : Location of diagoinal terms
```
"""
function fsparv!(kv::Vector{Float64}, km::Matrix{Float64},
  g::Vector{Int}, kdiag::Vector{Int})
  ndof = size(g, 1)
  for i in 1:ndof
    k = g[i]
    if k !== 0
      for j in 1:ndof
        if g[j] !== 0
          iw = k - g[j]
          if iw >= 0
            ival = kdiag[k] - iw
            kv[ival] += km[i, j]
          end
        end
      end
    end
  end
end


text:
using MarcuciusUtility
using Test

@testset "MarcuciusUtility.jl" begin
    # Write your tests here.
end


text:
# Animate DAMM in the REPL
# Not included in DAMMmodel module for now

using DAMMmodel
using Plots; unicodeplots()
df = DAMMfdata(20)
d = DAMMmat(df.Tₛ, df.θ, df.Rₛ, 20)

main() = begin
  anim = @animate for i = -180:10:180
    surface(d.x, d.y.*100, d.DAMM_Matrix, colormap = :jet, camera = (i, 30))
  end
  #gif(anim, "anim_fps5.gif", fps=5)
end

main()

# other method, similar result
for i = -180:1:180
  println(i)
  surface(d.x, d.y.*100, d.DAMM_Matrix, colormap = :jet, border = :dotted, camera = (i, 30)) |> display
  sleep(0.01)
  print("\e[2J")
end


text:
using BandedMatrices, LinearAlgebra, FillArrays, Test
import BandedMatrices: _BandedMatrix, DefaultBandedMatrix

@testset "misc tests" begin
    @testset "Diagonal of banded" begin
        D = Diagonal([BandedMatrix(Eye(5),(2,3)), BandedMatrix(Eye(6),(1,1))])
        @test @inferred(D[1,1]) == I
        @test @inferred(D[1,2]) == zeros(5,6)
        @test bandwidths(D[1,2]) == (2,1)
    end

    @time @testset "BandedMatrix of BandedMatrix" begin
        A = BandedMatrix{DefaultBandedMatrix{Float64}}(undef, (1, 2), (0, 1))
        A[1,1] = BandedMatrix(Eye(1),(0,1))
        A[1,2] = BandedMatrix(Zeros(1,2),(0,1))
        A[1,2][1,1] = -1/3
        A[1,2][1,2] = 1/3
        B = BandedMatrix{DefaultBandedMatrix{Float64}}(undef, (2, 1), (1, 1))
        B[1,1] = 0.2BandedMatrix(Eye(1),(0,1))
        B[2,1] = BandedMatrix(Zeros(2,1), (1,0))
        B[2,1][1,1] = -2/30
        B[2,1][2,1] = 1/3

        @test (A*B)[1,1][1,1] ≈ 1/3

        A = BandedMatrix(Diagonal([randn(2,3), randn(3,4)]), (0,0))
        @test @inferred(A[1,2]) == zeros(2,4)

        A = BandedMatrix(Diagonal([brand(5,6,2,1), brand(7,8,1,1)]), (0,0))
        @test @inferred(A[1,2]) == zeros(5,8)
    end

    @time @testset "dense overrides" begin
        A = rand(10,11)
        @test bandwidths(A) == (9,10)
        A = rand(10)
        @test bandwidths(A) == (9,0)
        @test bandwidths(A') == (0,9)
    end

    @time @testset "trivial convert routines" begin
        A = brand(3,4,1,2)
        @test isa(BandedMatrix{Float64}(A), BandedMatrix{Float64})
        @test isa(AbstractMatrix{Float64}(A), BandedMatrix{Float64})
        @test isa(AbstractArray{Float64}(A), BandedMatrix{Float64})
        @test isa(BandedMatrix(A), BandedMatrix{Float64})
        @test isa(AbstractMatrix(A), BandedMatrix{Float64})
        @test isa(AbstractArray(A), BandedMatrix{Float64})
        @test isa(BandedMatrix{ComplexF16}(A), BandedMatrix{ComplexF16})
        @test isa(AbstractMatrix{ComplexF16}(A), BandedMatrix{ComplexF16})
        @test isa(AbstractArray{ComplexF16}(A), BandedMatrix{ComplexF16})
    end
    @time @testset "show" begin
      @test occursin("10×10 BandedMatrix{Float64,Array{Float64,2},"*string(Base.OneTo{Int})*"}",
         sprint() do io
            show(io, MIME"text/plain"(), brand(10, 10, 3, 3))
         end)
        needle = "1.0  0.0   ⋅ \n 0.0  1.0  0.0\n  ⋅   0.0  1.0"
        @test occursin(needle, sprint() do io
             show(io, MIME"text/plain"(), BandedMatrix(Eye(3),(1,1)))
          end)
    end
    @time @testset "Issue #27" begin
        A=brand(1,10,0,9)
        B=brand(10,10,255,255)
        @test Matrix(A*B)  ≈ Matrix(A)*Matrix(B)
    end

    @testset "defaultdot" begin
        A = randn(5)
        @test BandedMatrices.dot(A,A) ≡ LinearAlgebra.dot(A,A)
    end

    @testset "Offset axes" begin
        A = _BandedMatrix(Ones((Base.OneTo(3),Base.Slice(-5:5),)), Base.Slice(-4:2), 1,1)
        @test size(A) == (7,11)
        @test A[-4,-5] == 1
        @test A[-4,-5+2] == 1
        @test A[-4,-5+3] == 0
        @test_throws BoundsError A[-5,-5+3]
    end
end


text:
using MagicMarshalling
using Test

is_simple_type(o) = false
is_simple_type(o::Int) = true
is_simple_type(o::Float64) = true
is_simple_type(o::String) = true

function is_simple_type(o::Vector)
    for x in o
        is_simple_type(x) || return false
    end
    return true
end
function is_simple_type(o::Dict)
    for (k,v) in pairs(o)
        k isa String || return false
        is_simple_type(v) || return false
    end
    return true
end

@testset "Basic bitstypes" begin
    for o in [1,1.0,-2,NaN,-Inf,]
        m = marshal(o)
        @test is_simple_type(m)
        o2 = unmarshal(m)
        @test o === o2
    end
end

@testset "NamedTuple" begin
    o = (a=1, b=2)
    m = marshal(o)
    @test is_simple_type(m)
    o2 = unmarshal(m)
    @test o === o2
end

@testset "Tuple $o" for o in [(), (1,), (1,2), ((1,2), (3,4))]
    m = marshal(o)
    @test is_simple_type(m)
    o2 = unmarshal(m)
    @test o === o2
end

@testset "Array" begin
    for o in [
                [1,2],
                [1 2; 3 4],
                ["a", 2],
                [],
                [[1]],
                [[1], [1,[3]]],
        ]
        m = marshal(o)
        @test is_simple_type(m)
        o2 = unmarshal(m)
        @test o == o2
    end
end

@testset "Structs" begin
    struct ParametricStruct{A,B}
        a::A
        b::B
    end

    for o in [
        ParametricStruct(1,2),
        ParametricStruct(ParametricStruct(1,2), ParametricStruct(1.0, 3)),
       ]
        m = marshal(o)
        @test is_simple_type(m)
        o2 = unmarshal(m)
        @test o === o2
    end
end


text:
module API # module H3

# types
export H3Index, GeoCoord, CoordIJ, Vec2d, Vec3d, CoordIJK, FaceIJK

# Indexing functions
export geoToH3, h3ToGeo, h3ToGeoBoundary

# Index inspection functions
export h3GetResolution, h3GetBaseCell, stringToH3, h3ToString, h3IsValid, h3IsResClassIII, h3IsPentagon

# Grid traversal functions
export kRing, maxKringSize, kRingDistances, hexRange, hexRangeDistances, hexRanges, hexRing, h3Line, h3LineSize, h3Distance, experimentalH3ToLocalIj, experimentalLocalIjToH3

# Hierarchical grid functions
export h3ToParent, h3ToChildren, maxH3ToChildrenSize, compact, uncompact, maxUncompactSize

# Unidirectional edge functions
export h3IndexesAreNeighbors, getH3UnidirectionalEdge, h3UnidirectionalEdgeIsValid, getOriginH3IndexFromUnidirectionalEdge, getDestinationH3IndexFromUnidirectionalEdge, getH3IndexesFromUnidirectionalEdge, getH3UnidirectionalEdgesFromHexagon, getH3UnidirectionalEdgeBoundary

# Miscellaneous H3 functions
export hexAreaKm2, hexAreaM2, edgeLengthKm, edgeLengthM, numHexagons, getRes0Indexes, res0IndexCount

# Coordinate Systems
export ijToIjk, ijkToHex2d, ijkToIj, ijkDistance, ijkNormalize, h3ToLocalIjk, h3ToFaceIjk, localIjkToH3, faceIjkToH3, hex2dToCoordIJK, geoToVec3d, geoToFaceIjk

using ..Lib
using .Lib: H3Index, GeoCoord, GeoBoundary, CoordIJ
using .Lib: Vec2d, Vec3d, CoordIJK, FaceIJK
using .Lib: H3_NULL

###
#
# * descriptions are taken from
#   - https://github.com/uber/h3/tree/master/docs/api
#   - https://github.com/uber/h3/blob/master/src/h3lib/include/h3api.h.in
#
###

# Indexing functions

"""
    geoToH3(location::GeoCoord, resolution::Int)::H3Index

Indexes the location at the specified resolution.
"""
function geoToH3(location::GeoCoord, resolution::Int)::H3Index
    Lib.geoToH3(Ref(location), resolution)
end

"""
    h3ToGeo(h::H3Index)::GeoCoord

Finds the centroid of the index.
"""
function h3ToGeo(h::H3Index)::GeoCoord
    refcenter = Ref{GeoCoord}()
    Lib.h3ToGeo(h, refcenter)
    refcenter[]
end

"""
    h3ToGeoBoundary(h::H3Index)::Vector{GeoCoord}

Finds the boundary of the index.
"""
function h3ToGeoBoundary(h::H3Index)::Vector{GeoCoord}
    refboundary = Ref{GeoBoundary}()
    Lib.h3ToGeoBoundary(h, refboundary)
    numVerts = refboundary[].numVerts
    verts = refboundary[].verts[1:numVerts]
    collect(verts)
end

# Index inspection functions

"""
    h3GetResolution(h::H3Index)::Int

Returns the resolution of the index.
"""
function h3GetResolution(h::H3Index)::Int
    Lib.h3GetResolution(h)
end

"""
    h3GetBaseCell(h::H3Index)::Int

Returns the base cell number of the index.
"""
function h3GetBaseCell(h::H3Index)::Int
    Lib.h3GetBaseCell(h)
end

"""
    stringToH3(str::String)::H3Index

Converts the string representation to H3Index (UInt64) representation.
"""
function stringToH3(str::String)::H3Index
    Lib.stringToH3(str)
end

"""
    h3ToString(h::H3Index)::String

Converts the H3Index representation of the index to the string representation.
"""
function h3ToString(h::H3Index)::String
    string(h, base=16)
end

"""
    h3IsValid(h::H3Index)::Bool

Returns `true` if this is a valid H3 index.
"""
function h3IsValid(h::H3Index)::Bool
    Bool(Lib.h3IsValid(h))
end

"""
    h3IsResClassIII(h::H3Index)::Bool

Returns `true` if this index has a resolution with Class III orientation.
"""
function h3IsResClassIII(h::H3Index)::Bool
    Bool(Lib.h3IsResClassIII(h))
end

"""
    h3IsPentagon(h::H3Index)::Bool

Returns `true` if this index represents a pentagonal cell.
"""
function h3IsPentagon(h::H3Index)::Bool
    Bool(Lib.h3IsPentagon(h))
end

# Grid traversal functions

"""
    kRing(origin::H3Index, k::Int)::Vector{H3Index}

k-rings produces indices within `k` distance of the origin index.
"""
function kRing(origin::H3Index, k::Int)::Vector{H3Index}
    array_len = Lib.maxKringSize(k)
    krings = Vector{H3Index}(undef, array_len)
    Lib.kRing(origin, k, krings)
    krings
end

"""
    maxKringSize(k::Int)::Int

Maximum number of indices that result from the kRing algorithm with the given `k`.
"""
function maxKringSize(k::Int)::Int
    Lib.maxKringSize(k)
end

"""
    kRingDistances(origin::H3Index, k::Int)::NamedTuple{(:out, :distances)}

k-rings produces indices within `k` distance of the origin index.
"""
function kRingDistances(origin::H3Index, k::Int)::NamedTuple{(:out, :distances)}
    array_len = Lib.maxKringSize(k)
    out = Vector{H3Index}(undef, array_len)
    distances = Vector{Cint}(undef, array_len)
    Lib.kRingDistances(origin, k, out, distances)
    (out = out, distances = distances)
end

"""
    hexRange(origin::H3Index, k::Int)::Vector{H3Index}

`hexRange` produces indexes within `k` distance of the origin index.
"""
function hexRange(origin::H3Index, k::Int)::Vector{H3Index}
    array_len = Lib.maxKringSize(k)
    out = Vector{H3Index}(undef, array_len)
    Lib.hexRange(origin, k, out)
    out
end

"""
    hexRangeDistances(origin::H3Index, k::Int)::NamedTuple{(:out, :distances)}

`hexRange` produces indexes within `k` distance of the origin index.
"""
function hexRangeDistances(origin::H3Index, k::Int)::NamedTuple{(:out, :distances)}
    array_len = Lib.maxKringSize(k)
    out = Vector{H3Index}(undef, array_len)
    distances = Vector{Cint}(undef, array_len)
    Lib.hexRangeDistances(origin, k, out, distances)
    (out = out, distances = distances)
end

"""
    hexRanges(h3Set::Vector{H3Index}, k::Int)::Vector{H3Index}

`hexRanges` takes an array of input hex IDs and a max k-ring and returns an array of hexagon IDs sorted first by the original hex IDs and then by the k-ring (0 to max), with no guaranteed sorting within each k-ring group.
"""
function hexRanges(h3Set::Vector{H3Index}, k::Int)::Vector{H3Index}
    array_len = Lib.maxKringSize(k)
    out = Vector{H3Index}(undef, array_len)
    Lib.hexRanges(h3Set, length(h3Set), k, out)
    out
end

"""
    hexRing(origin::H3Index, k::Int)::Vector{H3Index}

Produces the hollow hexagonal ring centered at `origin` with sides of length `k`.
"""
function hexRing(origin::H3Index, k::Int)::Vector{H3Index}
    out = Vector{H3Index}(undef, 6k)
    Lib.hexRing(origin, k, out)
    out
end

"""
    h3Line(origin::H3Index, destination::H3Index)::Vector{H3Index}

Given two H3 indexes, return the line of indexes between them (inclusive).
"""
function h3Line(origin::H3Index, destination::H3Index)::Vector{H3Index}
    line_size = Lib.h3LineSize(origin, destination)
    out = Vector{H3Index}(undef, line_size)
    Lib.h3Line(origin, destination, out)
    out
end

"""
    h3LineSize(origin::H3Index, destination::H3Index)::Int

Number of indexes in a line from the start index to the end index, to be used for allocating memory.
"""
function h3LineSize(origin::H3Index, destination::H3Index)::Int
    Lib.h3LineSize(origin, destination)
end

"""
    h3Distance(origin::H3Index, h::H3Index)::Int

Returns the distance in grid cells between the two indexes.
"""
function h3Distance(origin::H3Index, h::H3Index)::Int
    Lib.h3Distance(origin, h)
end

"""
    experimentalH3ToLocalIj(origin::H3Index, h::H3Index)::CoordIJ

Produces local IJ coordinates for an H3 index anchored by an origin.
"""
function experimentalH3ToLocalIj(origin::H3Index, h::H3Index)::CoordIJ
    refij = Ref{CoordIJ}()
    Lib.experimentalH3ToLocalIj(origin, h, refij)
    refij[]
end

"""
    experimentalLocalIjToH3(origin::H3Index, ij::CoordIJ)::H3Index

Produces an H3 index from local IJ coordinates anchored by an `origin`.
"""
function experimentalLocalIjToH3(origin::H3Index, ij::CoordIJ)::H3Index
    refh = Ref{H3Index}()
    Lib.experimentalLocalIjToH3(origin, Ref(ij), refh)
    refh[]
end

# Hierarchical grid functions

"""
    h3ToParent(h::H3Index, parentRes::Int)::H3Index

Returns the parent (coarser) index containing h.
"""
function h3ToParent(h::H3Index, parentRes::Int)::H3Index
    Lib.h3ToParent(h, parentRes)
end

"""
    h3ToChildren(h::H3Index, childRes::Int)::Vector{H3Index}

Populates children with the indexes contained by h at resolution childRes. children must be an array of at least size maxH3ToChildrenSize(h, childRes).
"""
function h3ToChildren(h::H3Index, childRes::Int)::Vector{H3Index}
    children_size = Lib.maxH3ToChildrenSize(h, childRes)
    children = Vector{H3Index}(undef, children_size)
    Lib.h3ToChildren(h, childRes, children)
    children
end

"""
    maxH3ToChildrenSize(h::H3Index, childRes::Int)::Int

Returns the size of the array needed by h3ToChildren for these inputs.
"""
function maxH3ToChildrenSize(h::H3Index, childRes::Int)::Int
    Lib.maxH3ToChildrenSize(h, childRes)
end

"""
    compact(h3Set::Vector{H3Index})::Vector{H3Index}

Compacts the set h3Set of indexes as best as possible, into the array compactedSet. compactedSet must be at least the size of h3Set in case the set cannot be compacted.
"""
function compact(h3Set::Vector{H3Index})::Vector{H3Index}
    numHexes = length(h3Set)
    compactedSet = fill(H3_NULL, numHexes)
    Lib.compact(h3Set, compactedSet, numHexes)
    compactedSet
end

"""
    uncompact(compactedSet::Vector{H3Index}, res::Int)::Vector{H3Index}

Uncompacts the set compactedSet of indexes to the resolution res.
"""
function uncompact(compactedSet::Vector{H3Index}, res::Int)::Vector{H3Index}
    hexCount = length(compactedSet)
    maxHexes = Lib.maxUncompactSize(compactedSet, hexCount, res)
    h3Set = Vector{H3Index}(undef, maxHexes)
    Lib.uncompact(compactedSet, hexCount, h3Set, maxHexes, res)
    h3Set
end

"""
    maxUncompactSize(compactedSet::Vector{H3Index}, res::Int)::Int

Returns the size of the array needed by uncompact.
"""
function maxUncompactSize(compactedSet::Vector{H3Index}, res::Int)::Int
    hexCount = length(compactedSet)
    Lib.maxUncompactSize(compactedSet, hexCount, res)
end

# Unidirectional edge functions

"""
    h3IndexesAreNeighbors(origin::H3Index, destination::H3Index)::Bool

Returns whether or not the provided H3Indexes are neighbors.
Returns `true` if the indexes are neighbors, `false` otherwise.
"""
function h3IndexesAreNeighbors(origin::H3Index, destination::H3Index)::Bool
    Bool(Lib.h3IndexesAreNeighbors(origin, destination))
end

"""
    getH3UnidirectionalEdge(origin::H3Index, destination::H3Index)::H3Index

Returns a unidirectional edge H3 index based on the provided origin and destination.
"""
function getH3UnidirectionalEdge(origin::H3Index, destination::H3Index)::H3Index
    Lib.getH3UnidirectionalEdge(origin, destination)
end

"""
    h3UnidirectionalEdgeIsValid(edge::H3Index)::Bool

Determines if the provided H3Index is a valid unidirectional edge index.
Returns `true` if it is a unidirectional edge H3Index, otherwise `false`.
"""
function h3UnidirectionalEdgeIsValid(edge::H3Index)::Bool
    Bool(Lib.h3UnidirectionalEdgeIsValid(edge))
end

"""
    getOriginH3IndexFromUnidirectionalEdge(edge::H3Index)::H3Index

Returns the origin hexagon from the unidirectional edge H3Index.
"""
function getOriginH3IndexFromUnidirectionalEdge(edge::H3Index)::H3Index
    Lib.getOriginH3IndexFromUnidirectionalEdge(edge)
end

"""
    getDestinationH3IndexFromUnidirectionalEdge(edge::H3Index)::H3Index

Returns the destination hexagon from the unidirectional edge H3Index.
"""
function getDestinationH3IndexFromUnidirectionalEdge(edge::H3Index)::H3Index
    Lib.getDestinationH3IndexFromUnidirectionalEdge(edge)
end

"""
    getH3IndexesFromUnidirectionalEdge(edge::H3Index)::Tuple{H3Index, H3Index}

Returns the origin, destination pair of hexagon IDs for the given edge ID, which are placed at originDestination[0] and originDestination[1] respectively.
"""
function getH3IndexesFromUnidirectionalEdge(edge::H3Index)::Tuple{H3Index, H3Index}
    originDestination = Vector{H3Index}(undef, 2)
    Lib.getH3IndexesFromUnidirectionalEdge(edge, originDestination)
    tuple(originDestination...)
end

"""
    getH3UnidirectionalEdgesFromHexagon(origin::H3Index)::Vector{H3Index}

Provides all of the unidirectional edges from the current H3Index. edges must be of length 6, and the number of undirectional edges placed in the array may be less than 6.
"""
function getH3UnidirectionalEdgesFromHexagon(origin::H3Index)::Vector{H3Index}
    edges = Vector{H3Index}(undef, 6)
    Lib.getH3UnidirectionalEdgesFromHexagon(origin, edges)
    edges
end

"""
    getH3UnidirectionalEdgeBoundary(edge::H3Index)::Vector{GeoCoord}

Provides the coordinates defining the unidirectional edge.
"""
function getH3UnidirectionalEdgeBoundary(edge::H3Index)::Vector{GeoCoord}
    refboundary = Ref{GeoBoundary}()
    Lib.getH3UnidirectionalEdgeBoundary(edge, refboundary)
    numVerts = refboundary[].numVerts
    verts = refboundary[].verts[1:numVerts]
    collect(verts)
end

# Miscellaneous H3 functions

"""
    hexAreaKm2(res::Int)::Cdouble

Average hexagon area in square kilometers at the given resolution.
"""
function hexAreaKm2(res::Int)::Cdouble
    Lib.hexAreaKm2(res)
end

"""
    hexAreaM2(res::Int)::Cdouble

Average hexagon area in square meters at the given resolution.
"""
function hexAreaM2(res::Int)::Cdouble
    Lib.hexAreaM2(res)
end

"""
    edgeLengthKm(res::Int)::Cdouble

Average hexagon edge length in kilometers at the given resolution.
"""
function edgeLengthKm(res::Int)::Cdouble
    Lib.edgeLengthKm(res)
end

"""
    edgeLengthM(res::Int)::Cdouble

Average hexagon edge length in meters at the given resolution.
"""
function edgeLengthM(res::Int)::Cdouble
    Lib.edgeLengthM(res)
end

"""
    numHexagons(res::Int)::Int

Number of unique H3 indexes at the given resolution.
"""
function numHexagons(res::Int)::Int
    Lib.numHexagons(res)
end

"""
    getRes0Indexes()::Vector{H3Index}

All the resolution 0 H3 indexes.
"""
function getRes0Indexes()::Vector{H3Index}
    out = Vector{H3Index}(undef, Lib.res0IndexCount())
    Lib.getRes0Indexes(out)
    out
end

"""
    res0IndexCount()::Cint

Number of resolution 0 H3 indexes.
"""
function res0IndexCount()::Int
    Lib.res0IndexCount()
end

### Coordinate Systems

"""
    ijToIjk(c::CoordIJ)::CoordIJK

Transforms coordinates from the IJ coordinate system to the IJK+ coordinate system.
"""
function ijToIjk(c::CoordIJ)::CoordIJK
    ref = Ref{CoordIJK}()
    ccall((:ijToIjk, Lib.libh3), Cvoid, (Ptr{CoordIJ}, Ptr{CoordIJK}), Ref(c), ref)
    ref[]
end

"""
    ijkToHex2d(c::CoordIJK)::Vec2d

Find the center point in 2D cartesian coordinates of a hex.
"""
function ijkToHex2d(c::CoordIJK)::Vec2d
    ref = Ref{Vec2d}()
    ccall((:_ijkToHex2d, Lib.libh3), Cvoid, (Ptr{CoordIJK}, Ptr{Vec2d}), Ref(c), ref)
    ref[]
end

"""
    ijkToIj(c::CoordIJK)::CoordIJ

Transforms coordinates from the IJK+ coordinate system to the IJ coordinate system.
"""
function ijkToIj(c::CoordIJK)::CoordIJ
    ref = Ref{CoordIJ}()
    ccall((:ijkToIj, Lib.libh3), Cvoid, (Ptr{CoordIJK}, Ptr{CoordIJ}), Ref(c), ref)
    ref[]
end

"""
    ijkDistance(c1::CoordIJK, c2::CoordIJK)::Int

Finds the distance between the two coordinates. Returns result.
"""
function ijkDistance(c1::CoordIJK, c2::CoordIJK)::Int
    ccall((:ijkDistance, Lib.libh3), Cint, (Ptr{CoordIJK}, Ptr{CoordIJK}), Ref(c1), Ref(c2))
end

"""
    ijkNormalize(c::CoordIJK)::CoordIJK

Normalizes ijk coordinates by setting the components to the smallest possible values. Works in place.
"""
function ijkNormalize(c::CoordIJK)::CoordIJK
    ref = Ref(c)
    ccall((:_ijkNormalize, Lib.libh3), Cvoid, (Ptr{CoordIJK},), ref)
    ref[]
end

"""
    h3ToLocalIjk(origin::H3Index, h3::H3Index)::CoordIJK

Produces ijk+ coordinates for an index anchored by an origin.
The coordinate space used by this function may have deleted regions or warping due to pentagonal distortion.
Coordinates are only comparable if they come from the same origin index.
Failure may occur if the index is too far away from the origin or if the index is on the other side of a pentagon.
"""
function h3ToLocalIjk(origin::H3Index, h3::H3Index)::CoordIJK
    ref = Ref{CoordIJK}()
    ccall((:h3ToLocalIjk, Lib.libh3), Cint, (H3Index, H3Index, Ptr{CoordIJK}), origin, h3, ref)
    ref[]
end

"""
    h3ToFaceIjk(h::H3Index)::FaceIJK

Convert an H3Index to a FaceIJK address.
"""
function h3ToFaceIjk(h::H3Index)::FaceIJK
    ref = Ref{FaceIJK}()
    ccall((:_h3ToFaceIjk, Lib.libh3), Cvoid, (H3Index, Ptr{FaceIJK}), h, ref)
    ref[]
end

"""
    localIjkToH3(origin::H3Index, ijk::CoordIJK)::H3Index

Produces an index for ijk+ coordinates anchored by an origin.
The coordinate space used by this function may have deleted regions or warping due to pentagonal distortion.
Failure may occur if the coordinates are too far away from the origin or if the index is on the other side of a pentagon.
"""
function localIjkToH3(origin::H3Index, ijk::CoordIJK)::H3Index
    ref = Ref{H3Index}()
    ccall((:localIjkToH3, Lib.libh3), Cint, (H3Index, Ptr{CoordIJK}, Ptr{H3Index}), origin, Ref(ijk), ref)
    ref[]
end

"""
    faceIjkToH3(faceijk::FaceIJK, res::Int)::H3Index

Convert an FaceIJK address to the corresponding H3Index.
"""
function faceIjkToH3(faceijk::FaceIJK, res::Int)::H3Index
    ccall((:_faceIjkToH3, Lib.libh3), H3Index, (Ptr{FaceIJK}, Cint), Ref(faceijk), res)
end

"""
    hex2dToCoordIJK(v::Vec2d)::CoordIJK

Determine the containing hex in ijk+ coordinates for a 2D cartesian coordinate vector (from DGGRID).
"""
function hex2dToCoordIJK(v::Vec2d)::CoordIJK
    ref = Ref{CoordIJK}()
    ccall((:_hex2dToCoordIJK, Lib.libh3), Cvoid, (Ptr{Vec2d}, Ptr{CoordIJK}), Ref(v), ref)
    ref[]
end

"""
    geoToVec3d(geo::GeoCoord)::Vec3d

Calculate the 3D coordinate on unit sphere from the latitude and longitude.
"""
function geoToVec3d(geo::GeoCoord)::Vec3d
    ref = Ref{Vec3d}()
    ccall((:_geoToVec3d, Lib.libh3), Cvoid, (Ptr{GeoCoord}, Ptr{Vec3d}), Ref(geo), ref)
    ref[]
end

"""
    geoToFaceIjk(geo::GeoCoord, res::Int)::FaceIJK

Encodes a coordinate on the sphere to the FaceIJK address of the containing cell at the specified resolution.
"""
function geoToFaceIjk(geo::GeoCoord, res::Int)::FaceIJK
    ref = Ref{FaceIJK}()
    ccall((:_geoToFaceIjk, Lib.libh3), Cvoid, (Ptr{GeoCoord}, Cint, Ptr{FaceIJK}), Ref(geo), res, ref)
    ref[]
end

### ≈

function Base.isapprox(A::Vector{GeoCoord}, B::Vector{GeoCoord})
    A === B && return true
    axes(A) != axes(B) && return false
    for (a, b) in zip(A, B)
        if !(isapprox(a.lat, b.lat) && isapprox(a.lon, b.lon))
            return false
        end
    end
    return true
end

function Base.isapprox(A::Vector{Tuple{Float64,Float64}}, B::Vector{Tuple{Float64,Float64}})
    A === B && return true
    axes(A) != axes(B) && return false
    for (a, b) in zip(A, B)
        if !(isapprox(a[1], b[1]) && isapprox(a[2], b[2]))
            return false
        end
    end
    return true
end

end # module H3.API


text:
using ChebyshevPicardIteration
using Test

@testset "ChebyshevPicardIteration.jl" begin
    # Write your tests here.
end


text:
@testset "$(format("Loss: Gamma[construction]"))" begin
    @test typeof(Loss{Gamma}(Gamma())) == Loss{Gamma}
    @test typeof(Loss(Gamma())) == Loss{Gamma}
    @test typeof(Loss(:Gamma)) == Loss{Gamma}
end

@testset "$(format("Optimizer: Gamma Loss [Small][Forgiving][Native][]"))"  begin
    @timeit to "Optimizer: Gamma Loss [Small][Forgiving][Native]"  begin
        let
            is_admissible(result)  = begin
                if result["relative-error[L2]"] < 1000
                    return true
                end
                @warn @sprintf("expected %f, got %f",0.1,result["relative-error[L2]"])
                return false
            end
            @test !isnothing(Loss(Gamma()))
            @test !isnothing(provide(Loss(:Gamma)))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 1000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 3000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 5000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 1000,
                                                        ρ                     = 0.1,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 3000,
                                                        ρ                     = 0.1,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(Gamma()) ,
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 5000,
                                                        ρ                     = 0.2,
                                                        step_size             = 0.005,
                                                        max_iter              = 100))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 1000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 3000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 5000,
                                                        ρ                     = 0,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 1000,
                                                        ρ                     = 0.1,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 3000,
                                                        ρ                     = 0.1,
                                                        step_size             =0.005,
                                                        max_iter              = 200))
            @test is_admissible(unit_test_train_subloss(Gamma(),gradient_eval = Loss(:Gamma),
                                                        input_distribution    = Distributions.Gamma(5,0.5),
                                                        input_size            = 5000,
                                                        ρ                     = 0.2,
                                                        step_size             = 0.005,
                                                        max_iter              = 200))
        end
    end
end

@warn "Autograd version of the gamma case is yet to be implemented due to numerical instability (forward pass goes into complex domain...)"


text:
# tutorial on conventional 2D SLAM example

# 0. Load required packages
using Caesar, RoME

# 1. Get the authentication and session information
include(joinpath(ENV["HOME"],"Documents","blandauthlocal.jl"))
backend_config, user_config = standardcloudgraphsetup(addrdict=addrdict)
user_config["session"] = "SLAM2D_TUTORIAL"
addrdict["robotId"] = "TestRobot"

# Start with an empty graph (local dictionary version) # fg = initfg(sessionname="SLAM2D_TUTORIAL")
## TODO -- ISSUE Julia 0.6.0-0.6.2 dives into some StackOverflow problem using the functions, but fine when called separately.
fg = Caesar.initfg(sessionname=user_config["session"], robotname=addrdict["robotId"], cloudgraph=backend_config)

# fg = RoME.initfg(sessionname=addrdict["session"])
# fg = IncrementalInference.emptyFactorGraph()
# user_config["session"]
# fg.sessionname
# fg.sessionname = user_config["session"]
# fg.cg = backend_config

# also add a PriorPose2 to pin the first pose at a fixed location
addNode!(fg, :x0, Pose2) # , labels=["POSE"]
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), 0.01*eye(3))) )

# ls(fg, :x0)
# isInitialized(fg, :x0)
# # to debug the zero lengths issue
# ensureAllInitialized!(fg)
# getVal(fg, :x0)
# Prior(MvNormal([0.0;0.0;0], diagm([1.0;1.0;0.01].^2)))

# Drive around in a hexagon
for i in 0:5
  psym = Symbol("x$i")
  nsym = Symbol("x$(i+1)")
  addNode!(fg, nsym, Pose2) # , labels=["VARIABLE";"POSE"]
  addFactor!(fg, [psym;nsym], Pose2Pose2(reshape([10.0;0;pi/3],3,1), 0.01*eye(3), [1.0]), autoinit=true )
  # Pose2Pose2_NEW(MvNormal([10.0;0;pi/3], diagm([0.1;0.1;0.1].^2)))
end

# Graphs.plot(fg.g)
isInitialized(fg, :x5)
isInitialized(fg, :x6)

ensureAllInitialized!(fg)

using RoMEPlotting

drawPoses(fg)

getVal(fg, :x0)
v = getVert(fg, :x0)
getVal(v)

importall CloudGraphs

exvid = fg.IDs[:x0]
neoID = fg.cgIDs[exvid]
cvr = CloudGraphs.get_vertex(fg.cg, neoID, false)
exv = CloudGraphs.cloudVertex2ExVertex(cvr)
getVal(exv)
# getData(exv)

getData(fg.g.vertices[1])

tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree)
inferOverTreeR!(fg, tree)

using RoMEPlotting, Gadfly

pl = plotKDE(fg, [:x0; :x1; :x2; :x3; :x4; :x5; :x6]);

Gadfly.draw(PDF("tmpX0123456.pdf", 15cm, 15cm), pl)

@async run(`evince tmpX0123456.pdf`)

# pl = drawPoses(fg)
pl = drawPosesLandms(fg)
Gadfly.draw(PDF("tmpPosesFg.pdf", 15cm, 15cm), pl)
@async run(`evince tmpPosesFg.pdf`)

tree = wipeBuildNewTree!(fg)

@async Graphs.plot(tree.bt)

@time inferOverTree!(fg, tree)

# Graphs.plot(tree.bt)

@async Graphs.plot(fg.g)


text:
function show_manifold_measure(io, mime, μ)
    name = nameof(typeof(μ))
    fields = (getfield(μ, n) for n in fieldnames(typeof(μ)) if n !== :par)
    field_string = sprint(show, mime, Tuple(fields))[2:(end - 1)]
    if endswith(field_string, ",")
        field_string = field_string[1:(end - 1)]
    end
    print(io, name, "(", field_string)
    par = params(μ)
    if !isempty(par)
        par_string = sprint(show, mime, par)[2:(end - 1)]
        if endswith(par_string, ",")
            par_string = par_string[1:(end - 1)]
        end
        print(io, "; ", par_string)
    end
    return print(io, ")")
end


text:
using PICDataStructures
using Test, SafeTestsets

@testset "PICDataStructures.jl" begin
    @safetestset "Code quality" begin include("quality.jl") end
    @safetestset "Grids" begin include("grids.jl") end
    @safetestset "Scalar quantities" begin include("scalar.jl") end
    @safetestset "Vector quantities" begin include("vector.jl") end
    @safetestset "Plotting" begin include("plots.jl") end
    @safetestset "Integrations" begin include("lasertypes.jl") end
end


text:
module Estimation

using Docile
@docstrings

export
    apes,
    quin_fernandes

include("frequency.jl")
include("APES.jl")

end # module


text:
__precompile__(true)

module KenyaCoVaccines

using DynamicHMC,Dates,OrdinaryDiffEq,DiffEqCallbacks,JLD2,Parameters,Distributions
using NamedArrays,TransformVariables,LogDensityProblems,LinearAlgebra,BlackBoxOptim,Random
using DynamicHMC.Diagnostics,MCMCDiagnostics,MCMCChains, Interpolations, FileIO, CSV, DataFrames
using RecursiveArrayTools
using DataInterpolations
using StatsPlots
using Plots.PlotMeasures
using ForwardDiff

relative_testing_rate = vcat([0.0 for d in Date(2020,2,21):Day(1):Date(2020,3,14)],
        [t/120 for (t,d) in enumerate(Date(2020,3,15):Day(1):Date(2020,7,12))],
        [1.0 for k = 1:500])

include("transmissionmodel.jl")
include("types.jl")
include("likelihood.jl")
include("priors.jl")
include("inference.jl")
include("dataflow.jl")
include("vacc_scenarios.jl")

#Load defaults for the PCR and serological sensitivity after infection
@load("data/rel_sero_detection_after_infection.jld2")
@load("data/smoothedPCR_detection_after_infection.jld2")

#Load the prior for the p_test and chi parameters
@load("data/non_urban_detection_fits.jld2")

#@load("data/semi_urban_rural_detection_fits.jld2")
#@load("data/rural_fits.jld2")

end # module


text:
@testset "377.combination-sum-iv.jl" begin
    @test combination_sum4([1, 2, 3], 4) == 7
    @test combination_sum4([9], 3) == 0
    @test combination_sum4([10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140,
                            150, 160, 170, 180, 190, 200, 210, 220, 230, 240, 250, 260, 270,
                            280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380, 390, 400,
                            410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530,
                            540, 550, 560, 570, 580, 590, 600, 610, 620, 630, 640, 650, 660,
                            670, 680, 690, 700, 710, 720, 730, 740, 750, 760, 770, 780, 790,
                            800, 810, 820, 830, 840, 850, 860, 870, 880, 890, 900, 910, 920,
                            930, 940, 950, 960, 970, 980, 990, 111], 999) == 1
end


text:
using MetacommunityDynamics
using DynamicGrids
using DynamicGridsGtk
using Dispersal
using Distributions
using ColorSchemes
using Plots

consumermodel =
    Eating{Tuple{:C,:R}}(
        functionalresponse=LotkaVolterra(0.3),
        dt=0.1) +
    LinearMortality{:C}(0.1);

resourcemodel =
    MetacommunityDynamics.LogisticGrowth{:R}(λ=1.4, K=10., dt=0.1);

model = resourcemodel + consumermodel

latticesize = 1
initconsumer = zeros(Biomass, latticesize,latticesize)
initresource = zeros(Biomass, latticesize,latticesize)

initresource[1] = 2.
initconsumer[1] = 0.5

arrayout = ArrayOutput((C=initconsumer, R=initresource ), tspan=1:500)
@time sim!(arrayout, model)

Cs = Float64[]
Rs = Float64[]
for t in arrayout
    push!(Cs, sum(t[:C]))
    push!(Rs, sum(t[:R]))
end

plot(1:length(Cs),Cs, label="consumer")
plot!(1:length(Rs),Rs, label="resource")
xlabel!("time")
ylabel!("biomass")