text:
using Test

using EquationsOfState.FiniteStrains
using EquationsOfState.LinearFitting

@testset "Test strain-volume-derivative" begin
    volumes = [
        25.987454833,
        26.9045702104,
        27.8430241908,
        28.8029649591,
        29.7848370694,
        30.7887887064,
        31.814968055,
        32.8638196693,
        33.9353435494,
        35.0299842495,
        36.1477417695,
        37.2892088485,
        38.4543854865,
        39.6437162376,
        40.857201102,
        42.095136449,
        43.3579668329,
        44.6456922537,
        45.9587572656,
        47.2973100535,
        48.6614988019,
        50.0517680652,
        51.4682660281,
        52.9112890601,
        54.3808371612,
        55.8775030703,
        57.4014349722,
        58.9526328669,
    ]
    v0 = 40.98926572870838
    @test [strain_volume_derivative(EulerianStrain(), v0, v, 3) for v in volumes] ≈ [
        -0.00025163164656036914,
        -0.00022158317701314234,
        -0.0001954064296869249,
        -0.00017256893337511754,
        -0.00015261005948744114,
        -0.00013514276136071226,
        -0.00011983495896811362,
        -0.00010639785139417004,
        -9.458902516425933e-5,
        -8.419491823560175e-5,
        -7.503571748024374e-5,
        -6.695156194310254e-5,
        -5.980866702052554e-5,
        -5.348840953851233e-5,
        -4.789028506457302e-5,
        -4.2925492613058147e-5,
        -3.85164573279167e-5,
        -3.459717265567486e-5,
        -3.110870879523787e-5,
        -2.8000473828556818e-5,
        -2.5228159545528874e-5,
        -2.275242969867225e-5,
        -2.0539381634855355e-5,
        -1.855902166896191e-5,
        -1.6785363112484185e-5,
        -1.5194859884358581e-5,
        -1.3767304534537318e-5,
        -1.2484962604556655e-5,
    ]
    @test [strain_volume_derivative(LagrangianStrain(), v0, v, 3) for v in volumes] ≈ [
        -4.360792585446927e-5,
        -4.021795998862133e-5,
        -3.7125806529670504e-5,
        -3.430262724447532e-5,
        -3.172185978064063e-5,
        -2.9360589247189284e-5,
        -2.719823329453535e-5,
        -2.5215761270023045e-5,
        -2.3396932620829496e-5,
        -2.172639361674028e-5,
        -2.0191013020817355e-5,
        -1.8778177794469247e-5,
        -1.7477278166117428e-5,
        -1.6278235883744374e-5,
        -1.5172392867363944e-5,
        -1.4151626497343523e-5,
        -1.3208496564809248e-5,
        -1.2336591903572137e-5,
        -1.1529793498552908e-5,
        -1.0782740265865238e-5,
        -1.009054720444801e-5,
        -9.448630234274912e-6,
        -8.852961760583228e-6,
        -8.299807942932397e-6,
        -7.78586914649821e-6,
        -7.307932527332951e-6,
        -6.863217457745042e-6,
        -6.449212613211634e-6,
    ]
    @test [strain_volume_derivative(NaturalStrain(), v0, v, 3) for v in volumes] ≈ [
        0.00011395640842831919,
        0.00010269559592471744,
        9.265754335710983e-5,
        8.369864119123017e-5,
        7.569101157907677e-5,
        6.852549573283301e-5,
        6.210628060176279e-5,
        5.634765917870933e-5,
        5.117681395143497e-5,
        4.6527548395465794e-5,
        4.234347362622052e-5,
        3.857273983853629e-5,
        3.517161972472347e-5,
        3.2100139446683775e-5,
        2.9324064018820016e-5,
        2.681230993011352e-5,
        2.4537102053587615e-5,
        2.2474565886168472e-5,
        2.0602748703699993e-5,
        1.890256196526239e-5,
        1.7356955567635e-5,
        1.5950405389513548e-5,
        1.466936735556791e-5,
        1.3501591655278468e-5,
        1.2436332606671682e-5,
        1.1463547171440715e-5,
        1.0574546405738849e-5,
        9.761585855581402e-6,
    ]
    @test [strain_volume_derivative(InfinitesimalStrain(), v0, v, 3) for v in volumes] ≈ [
        0.00017195382081596378,
        0.00015318072693119346,
        0.0001366374157858279,
        0.00012203947625819057,
        0.00010913738901171863,
        9.77197419230021e-5,
        8.76030749059566e-5,
        7.862563981068296e-5,
        7.065075730232455e-5,
        6.355618910218268e-5,
        5.723834723338108e-5,
        5.160365044520063e-5,
        4.657341353283599e-5,
        4.207684089417924e-5,
        3.80535882418961e-5,
        3.444963171076881e-5,
        3.1217248819476426e-5,
        2.8315601024208774e-5,
        2.5707710435800747e-5,
        2.3361614304347095e-5,
        2.1249043437552616e-5,
        1.9344592130116945e-5,
        1.7626220489559587e-5,
        1.6074217776270604e-5,
        1.4671393271712412e-5,
        1.3401940487277108e-5,
        1.2252230465391096e-5,
        1.121020710562245e-5,
    ]
end

ext: jl


text:
export isjwt

function isjwt(str::AbstractString)::Bool
    splitted = split(str, ".")

    if length(splitted) > 3 || length(splitted) < 2
        return false
    end

    return all(Validator.isbase64.(splitted, true))
end


text:
# Basic Julia implementation of adaptive shrinkage with normal
# likelihood and mixture-of-centered-normals prior. The outputs are
# the fitted mixture weights (x) and the runtime of the likelihood
# computations (t_likelihood), model fitting (t_fit) and posterior
# computations (t_posterior).
function ash(x, s; method = "mixSQP", gridmult = 1.3, lowrank = "svd",
             pqrtol = 1e-8)

  # Select the variances of the normal mixture components, then
  # compute the n x m conditional likelihood matrix.
  out, t_likelihood, bytes, gctime, memallocs = @timed begin
    sd = autoselectmixsd(x,s,gridmult = gridmult);
    L  = normlikmatrix(x,s,sd = sd);
  end

  # Fit the model using either the SQP or interior point solvers.
  if method == "mixSQP"
    fit, t_fit, bytes, gctime, memallocs = @timed mixSQP(L,pqrtol = pqrtol,
      lowrank = lowrank,eps = 1e-6,verbose = false);
    w = fit["x"];
  elseif method == "REBayes"
    w, t_fit = REBayes(L);
  end

  # Posterior calculations.
  out, t_posterior, bytes, gctime, memallocs = @timed begin
    ind = findall(w .> 0);
    ps2 = sd[ind].^2;
    s2  = s.^2;
    t   = s2 .+ ps2';
    comp_post_mean  = (x * ps2') ./ t;
    comp_post_sd2   = (s2 * ps2') ./ t;
    comp_post_mean2 = comp_post_sd2 .+ comp_post_mean.^2;
    comp_post_prob  = L[:,ind] .* w[ind]';
    comp_post_prob  = comp_post_prob ./ sum(comp_post_prob,dims = 2);
    post_mean       = sum(comp_post_prob .* comp_post_mean,dims = 2);
    post_mean2      = sum(comp_post_prob .* comp_post_mean2,dims = 2);
  end

  return Dict([("x",x),
               ("timing-likelihood",t_likelihood),
               ("timing-fit",t_fit),
               ("timing-posterior",t_posterior)])
end

max_stars_repo_name: stephenslab/mixsqp-paper

max_issues_repo_name: stephenslab/mixsqp-paper

max_forks_repo_name: stephenslab/mixsqp-paper


text:
"Holds the tableau of an variational partitioned additive Runge-Kutta method."
struct TableauVPARK{T} <: AbstractTableau{T}
    name::Symbol
    o::Int
    s::Int
    r::Int

    q::CoefficientsARK{T}
    p::CoefficientsARK{T}

    q̃::CoefficientsPRK{T}
    p̃::CoefficientsPRK{T}

    λ::CoefficientsMRK{T}

    d::Vector{T}

    function TableauVPARK{T}(name, o, s, r, q, p, q̃, p̃, λ, d) where {T}
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)

        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s==length(d)
        @assert r==q.r==p.r==q̃.r==p̃.r==λ.r

        new(name, o, s, r, q, p, q̃, p̃, λ, d)
    end

    function TableauVPARK{T}(name, o, s, r, q, p, q̃, p̃, λ) where {T}
        @assert isa(name, Symbol)
        @assert isa(s, Integer)
        @assert isa(r, Integer)
        @assert isa(o, Integer)

        @assert s > 0 "Number of stages s must be > 0"
        @assert r > 0 "Number of stages r must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s
        @assert r==q.r==p.r==q̃.r==p̃.r==λ.r

        new(name, o, s, r, q, p, q̃, p̃, λ)
    end
end

function TableauVPARK(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        α_q::Matrix{T}, α_p::Matrix{T},
                        a_q̃::Matrix{T}, a_p̃::Matrix{T},
                        α_q̃::Matrix{T}, α_p̃::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        β_q::Vector{T}, β_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T},
                        c_λ::Vector{T}, d_λ::Vector{T},
                        d::Vector{T}) where {T <: Real}

    s = length(c_q)
    r = length(c_λ)

    @assert s > 0 "Number of stages s must be > 0"
    @assert r > 0 "Number of stages r must be > 0"

    @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
    @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
    @assert s==size(α_q,1)==size(α_p,1)
    @assert r==size(α_q,2)==size(α_p,2)
    @assert s==length(d)
    @assert r==length(c_λ)==length(d_λ)
    @assert r==size(a_q̃,1)==size(a_p̃,1)
    @assert s==size(a_q̃,2)==size(a_p̃,2)
    @assert r==size(α_q̃,1)==size(α_q̃,2)==length(β_q)
    @assert r==size(α_p̃,1)==size(α_p̃,2)==length(β_p)

    q = CoefficientsARK{T}(name, order, s, r, a_q, b_q, c_q, α_q, β_q)
    p = CoefficientsARK{T}(name, order, s, r, a_p, b_p, c_p, α_p, β_p)
    q̃ = CoefficientsPRK{T}(name, order, s, r, a_q̃, c_λ, α_q̃)
    p̃ = CoefficientsPRK{T}(name, order, s, r, a_p̃, c_λ, α_p̃)
    λ = CoefficientsMRK{T}(name, r, d_λ, c_λ)

    TableauVPARK{T}(name, order, s, r, q, p, q̃, p̃, λ, d)
end

function TableauVPARK(name::Symbol, order::Int,
                        a_q::Matrix{T}, a_p::Matrix{T},
                        α_q::Matrix{T}, α_p::Matrix{T},
                        a_q̃::Matrix{T}, a_p̃::Matrix{T},
                        α_q̃::Matrix{T}, α_p̃::Matrix{T},
                        b_q::Vector{T}, b_p::Vector{T},
                        β_q::Vector{T}, β_p::Vector{T},
                        c_q::Vector{T}, c_p::Vector{T},
                        c_λ::Vector{T}, d_λ::Vector{T}) where {T <: Real}

    s = length(c_q)
    r = length(c_λ)

    @assert s > 0 "Number of stages s must be > 0"
    @assert r > 0 "Number of stages r must be > 0"

    @assert s==size(a_q,1)==size(a_q,2)==length(b_q)==length(c_q)
    @assert s==size(a_p,1)==size(a_p,2)==length(b_p)==length(c_p)
    @assert s==size(α_q,1)==size(α_p,1)
    @assert r==size(α_q,2)==size(α_p,2)
    @assert r==length(c_λ)==length(d_λ)
    @assert r==size(a_q̃,1)==size(a_p̃,1)
    @assert s==size(a_q̃,2)==size(a_p̃,2)
    @assert r==size(α_q̃,1)==size(α_q̃,2)==length(β_q)
    @assert r==size(α_p̃,1)==size(α_p̃,2)==length(β_p)

    q = CoefficientsARK{T}(name, order, s, r, a_q, b_q, c_q, α_q, β_q)
    p = CoefficientsARK{T}(name, order, s, r, a_p, b_p, c_p, α_p, β_p)
    q̃ = CoefficientsPRK{T}(name, order, s, r, a_q̃, c_λ, α_q̃)
    p̃ = CoefficientsPRK{T}(name, order, s, r, a_p̃, c_λ, α_p̃)
    λ = CoefficientsMRK{T}(name, r, d_λ, c_λ)

    TableauVPARK{T}(name, order, s, r, q, p, q̃, p̃, λ)
end

# TODO function readTableauVPARKFromFile(dir::AbstractString, name::AbstractString)

"Parameters for right-hand side function of variational partitioned additive Runge-Kutta methods."
mutable struct ParametersVPARK{DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT} <: Parameters{DT,TT}
    f_ϑ::ϑT
    f_f::FT
    f_u::UT
    f_g::GT
    f_ϕ::ϕT

    Δt::TT

    t_q::CoefficientsARK{TT}
    t_p::CoefficientsARK{TT}
    t_q̃::CoefficientsPRK{TT}
    t_p̃::CoefficientsPRK{TT}
    t_λ::CoefficientsMRK{TT}
    d_v::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function ParametersVPARK{DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT}(f_ϑ, f_f, f_u, f_g, f_ϕ, Δt, t_q, t_p, t_q̃, t_p̃, t_λ, d_v) where {DT,TT,D,S,R,ϑT,FT,UT,GT,ϕT}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new(f_ϑ, f_f, f_u, f_g, f_ϕ, Δt,
            t_q, t_p, t_q̃, t_p̃, t_λ, d_v,
            zero(TT), q, p, λ)
    end
end

"""
Variational partitioned additive Runge-Kutta integrator cache.

### Fields

* `n`: time step number
* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `p`: current solution of p
* `p̅`: previous solution of p
* `v`: vector field of q
* `v̅`: vector field of q̅
* `f`: vector field of p
* `f̅`: vector field of p̅
* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: vector field of internal stages of q
* `Z`: vector field of internal stages of p
"""
mutable struct IntegratorCacheVPARK{ST,TT,D,S,R} <: IDAEIntegratorCache{ST,D}
    n::Int
    t::TT
    t̅::TT

    q::Vector{ST}
    q̅::Vector{ST}
    p::Vector{ST}
    p̅::Vector{ST}
    λ::Vector{ST}
    λ̅::Vector{ST}
    μ::Vector{ST}
    μ̅::Vector{ST}

    qₑᵣᵣ::Vector{ST}
    pₑᵣᵣ::Vector{ST}

    v::Vector{ST}
    v̅::Vector{ST}
    f::Vector{ST}
    f̅::Vector{ST}
    u::Vector{ST}
    u̅::Vector{ST}
    g::Vector{ST}
    g̅::Vector{ST}

    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    s̃::Vector{ST}

    Qi::Vector{Vector{ST}}
    Pi::Vector{Vector{ST}}
    Vi::Vector{Vector{ST}}
    Fi::Vector{Vector{ST}}
    Yi::Vector{Vector{ST}}
    Zi::Vector{Vector{ST}}
    Φi::Vector{Vector{ST}}

    Qp::Vector{Vector{ST}}
    Pp::Vector{Vector{ST}}
    Λp::Vector{Vector{ST}}
    Up::Vector{Vector{ST}}
    Gp::Vector{Vector{ST}}
    Yp::Vector{Vector{ST}}
    Zp::Vector{Vector{ST}}
    Φp::Vector{Vector{ST}}

    function IntegratorCacheVPARK{ST,TT,D,S,R}() where {ST,TT,D,S,R}
        q = zeros(ST,D)
        q̅ = zeros(ST,D)
        p = zeros(ST,D)
        p̅ = zeros(ST,D)
        λ = zeros(ST,D)
        λ̅ = zeros(ST,D)
        μ = zeros(ST,D)
        μ̅ = zeros(ST,D)

        # create error vectors
        qₑᵣᵣ = zeros(ST,D)
        pₑᵣᵣ = zeros(ST,D)

        # create update vectors
        v = zeros(ST,D)
        v̅ = zeros(ST,D)
        f = zeros(ST,D)
        f̅ = zeros(ST,D)
        u = zeros(ST,D)
        u̅ = zeros(ST,D)
        g = zeros(ST,D)
        g̅ = zeros(ST,D)

        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Qi = create_internal_stage_vector(ST, D, S)
        Pi = create_internal_stage_vector(ST, D, S)
        Vi = create_internal_stage_vector(ST, D, S)
        Fi = create_internal_stage_vector(ST, D, S)
        Yi = create_internal_stage_vector(ST, D, S)
        Zi = create_internal_stage_vector(ST, D, S)
        Φi = create_internal_stage_vector(ST, D, S)

        Qp = create_internal_stage_vector(ST, D, R)
        Pp = create_internal_stage_vector(ST, D, R)
        Λp = create_internal_stage_vector(ST, D, R)
        Up = create_internal_stage_vector(ST, D, R)
        Gp = create_internal_stage_vector(ST, D, R)
        Yp = create_internal_stage_vector(ST, D, R)
        Zp = create_internal_stage_vector(ST, D, R)
        Φp = create_internal_stage_vector(ST, D, R)

        new(0, zero(TT), zero(TT), q, q̅, p, p̅, λ, λ̅, μ, μ̅,
                                   qₑᵣᵣ, pₑᵣᵣ,
                                   v, v̅, f, f̅, u, u̅, g, g̅,
                                   q̃, p̃, ṽ, f̃, s̃,
                                   Qi, Pi, Vi, Fi, Yi, Zi, Φi,
                                   Qp, Pp, Λp, Up, Gp, Yp, Zp, Φp)
    end
end

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheVPARK{ST,TT,D,S,R},
                                        params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    local tpᵢ::TT
    local tλᵢ::TT

    for i in 1:S
        for k in 1:D
            # copy x to Y, Z
            cache.Yi[i][k] = x[3*(D*(i-1)+k-1)+1]
            cache.Zi[i][k] = x[3*(D*(i-1)+k-1)+2]
            cache.Vi[i][k] = x[3*(D*(i-1)+k-1)+3]

            # compute Q and P
            cache.Qi[i][k] = params.q[k] + params.Δt * cache.Yi[i][k]
            cache.Pi[i][k] = params.p[k] + params.Δt * cache.Zi[i][k]
        end

        # compute f(X)
        tpᵢ = params.t + params.Δt * params.t_p.c[i]
        params.f_ϑ(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Φi[i])
        params.f_f(tpᵢ, cache.Qi[i], cache.Vi[i], cache.Fi[i])

        cache.Φi[i] .-= cache.Pi[i]
    end

    for i in 1:R
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+2]
            cache.Λp[i][k] = x[3*D*S+3*(D*(i-1)+k-1)+3]

            # compute Q and V
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        tλᵢ = params.t + params.Δt * params.t_λ.c[i]
        params.f_u(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Up[i])
        params.f_g(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Λp[i], cache.Gp[i])
        params.f_ϕ(tλᵢ, cache.Qp[i], cache.Pp[i], cache.Φp[i])
    end

    if length(params.d_v) > 0
        for k in 1:D
            cache.μ[k] = x[3*D*S+3*D*R+k]
        end
    end
end

"Compute stages of variational partitioned additive Runge-Kutta methods."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersVPARK{DT,TT,D,S,R}) where {ST,DT,TT,D,S,R}
    cache = IntegratorCacheVPARK{ST,TT,D,S,R}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:S
            for k in 1:D
                b[3*(D*(i-1)+k-1)+1] = - $cache.Yi[i][k]
                b[3*(D*(i-1)+k-1)+2] = - $cache.Zi[i][k]
                b[3*(D*(i-1)+k-1)+3] = - $cache.Φi[i][k]
                for j in 1:S
                    b[3*(D*(i-1)+k-1)+1] += params.t_q.a[i,j] * $cache.Vi[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.t_p.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*(D*(i-1)+k-1)+1] += params.t_q.α[i,j] * $cache.Up[j][k]
                    b[3*(D*(i-1)+k-1)+2] += params.t_p.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ]
        for i in 1:R
            for k in 1:D
                b[3*D*S+3*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[3*D*S+3*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                b[3*D*S+3*(D*(i-1)+k-1)+3] = - $cache.Φp[i][k]
                for j in 1:S
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.a[i,j] * $cache.Vi[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.a[i,j] * $cache.Fi[j][k]
                end
                for j in 1:R
                    b[3*D*S+3*(D*(i-1)+k-1)+1] += params.t_q̃.α[i,j] * $cache.Up[j][k]
                    b[3*D*S+3*(D*(i-1)+k-1)+2] += params.t_p̃.α[i,j] * $cache.Gp[j][k]
                end
            end
        end

        # compute b = - [Λ₁-λ]
        if params.t_λ.c[1] == 0
            for k in 1:D
                b[3*D*S+3*(k-1)+3] = - $cache.Λp[1][k] + params.λ[k]
            end
        end

        if length(params.d_v) > 0
            for i in 1:S
                for k in 1:D
                    b[3*(D*(i-1)+k-1)+3] -= $cache.μ[k] * params.d_v[i]
                end
            end

            for k in 1:D
                b[3*D*S+3*D*R+k] = 0
                for i in 1:S
                    b[3*D*S+3*D*R+k] -= $cache.Vi[i][k] * params.d_v[i]
                end
            end
        end
    end
end

@doc raw"""
Variational partitioned additive Runge-Kutta integrator.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} U_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} a_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \alpha_{ij} G_{n,j} , & i &= 1, ..., s , \\
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} V_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} U_{n,j} , & i &= 1, ..., r , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \tilde{a}_{ij} F_{n,j} + h \sum \limits_{j=1}^{r} \tilde{\alpha}_{ij} G_{n,j} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= 0 , & i &= 1, ..., r ,
\end{align}
```
with definitions
```math
\begin{align}
P_{n,i} &= \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
F_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
U_{n,i} &= \hphantom{-} \frac{\partial \phi}{\partial p} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
G_{n,i} &=           -  \frac{\partial \phi}{\partial q} (\tilde{Q}_{n,i}, \tilde{P}_{n,i})^{T} \Lambda_{n,i} , & i &= 1, ..., r , \\
\tilde{\Phi}_{n,i} &= \phi(\tilde{Q}_{n,i}, \tilde{P}_{n,i}) , & i &= 1, ..., r ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} V_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} U_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b_{i} F_{n,i} + h \sum \limits_{i=1}^{r} \beta_{i} G_{n,i} .
\end{align}
```
"""
struct IntegratorVPARK{DT, TT, ET <: IDAE{DT,TT},
                               PT <: ParametersVPARK{DT,TT},
                               ST <: NonlinearSolver{DT},
                               IT <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorSPARK{DT, TT}
    equation::ET
    tableau::TableauVPARK{TT}

    params::PT
    solver::ST
    iguess::IT
end

function IntegratorVPARK(equation::IDAE{DT,TT,FT,PT,UT,GT,ϕT,VT},
                         tableau::TableauVPARK{TT}, Δt::TT) where {DT,TT,FT,PT,UT,GT,ϕT,VT}
    D = equation.d
    S = tableau.s
    R = tableau.r

    N = 3*D*S + 3*D*R

    if isdefined(tableau, :d)
        N += D
        d_v = tableau.d
    else
        d_v = DT[]
    end

    # create params
    params = ParametersVPARK{DT,TT,D,S,R,FT,PT,UT,GT,ϕT}(
                                equation.ϑ, equation.f, equation.u, equation.g, equation.ϕ, Δt,
                                tableau.q, tableau.p, tableau.q̃, tableau.p̃, tableau.λ, d_v)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create integrator
    IntegratorVPARK{DT, TT, typeof(equation), typeof(params), typeof(solver), typeof(iguess)}(
                                        equation, tableau, params, solver, iguess)
end

equation(int::IntegratorVPARK) = int.equation
timestep(int::IntegratorVPARK) = int.params.Δt
tableau(int::IntegratorVPARK) = int.tableau
nstages(int::IntegratorVPARK) = int.tableau.s
pstages(int::IntegratorVPARK) = int.tableau.r

function create_integrator_cache(int::IntegratorVPARK{DT,TT}) where {DT,TT}
    IntegratorCacheVPARK{DT, TT, ndims(int), nstages(int), pstages(int)}()
end

function update_params!(params::ParametersVPARK, cache::IntegratorCacheVPARK)
    # set time for nonlinear solver and copy previous solution
    params.t  = cache.t
    params.q .= cache.q
    params.p .= cache.p
    params.λ .= cache.λ
end

function initialize!(int::IntegratorVPARK, cache::IntegratorCacheVPARK)
    cache.t̅ = cache.t - timestep(int)

    equation(int).v(cache.t, cache.q, cache.p, cache.v)
    equation(int).f(cache.t, cache.q, cache.p, cache.f)

    initialize!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f,
                            cache.t̅, cache.q̅, cache.p̅, cache.v̅, cache.f̅)
end

function initial_guess!(int::IntegratorVPARK, cache::IntegratorCacheVPARK)
    for i in 1:nstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q.c[i], tableau(int).p.c[i])

        for k in 1:ndims(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[3*(ndims(int)*(i-1)+k-1)+3] = cache.ṽ[k]
        end
    end

    for i in 1:pstages(int)
        evaluate!(int.iguess, cache.q, cache.p, cache.v, cache.f,
                              cache.q̅, cache.p̅, cache.v̅, cache.f̅,
                              cache.q̃, cache.p̃, cache.ṽ, cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+1] = (cache.q̃[k] - cache.q[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+2] = (cache.p̃[k] - cache.p[k])/timestep(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(ndims(int)*(i-1)+k-1)+3] = 0
        end
    end

    if int.params.t_λ.c[1] == 0
        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*(k-1)+3] = cache.λ[k]
        end
    end

    if isdefined(tableau(int), :d)
        for k in 1:ndims(int)
            int.solver.x[3*ndims(int)*nstages(int)+3*ndims(int)*pstages(int)+k] = 0
        end
    end
end

"Integrate DAE with variational partitioned additive Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPARK{DT,TT}, cache::IntegratorCacheVPARK{DT,TT}) where {DT,TT}
    # update nonlinear solver parameters from cache
    update_params!(int.params, cache)

    # compute initial guess
    initial_guess!(int, cache)

    # reset cache
    reset!(cache, timestep(int))

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params, cache.n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params, cache.n)

    # compute vector fields at internal stages
    compute_stages!(int.solver.x, cache, int.params)

    # compute final update
    update_solution!(cache.q, cache.qₑᵣᵣ, cache.Vi, int.params.t_q.b, timestep(int))
    update_solution!(cache.p, cache.pₑᵣᵣ, cache.Fi, int.params.t_p.b, timestep(int))

    # compute projection
    update_solution!(cache.q, cache.qₑᵣᵣ, cache.Up, int.params.t_q.β, timestep(int))
    update_solution!(cache.p, cache.pₑᵣᵣ, cache.Gp, int.params.t_p.β, timestep(int))
    update_multiplier!(cache.λ, cache.Λp, int.params.t_λ.b)

    # copy solution to initial guess
    update!(int.iguess, cache.t, cache.q, cache.p, cache.v, cache.f)

    # take care of periodic solutions
    cut_periodic_solution!(cache, equation(int).periodicity)
end


text:
using NPZ
using Plots
include("plot_utils.jl")
include("voronoi_utils.jl")
include("lambda_continuum.jl")
include("irregular_ray_tracing.jl")

global my_seed = 1998
Random.seed!(my_seed)

pyplot()

function plot_sites(sites::VoronoiSites)
    z = ustrip.(sites.positions[1,:])
    x = ustrip.(sites.positions[2,:])
    y = ustrip.(sites.positions[3,:])
    scatter(x/1e6,
            y/1e6,
            z/1e6,
            title="Sites",
            dpi=300)
    savefig("../img/sites.png")
end

function compare(DATA, quadrature)
    maxiter = 50
    ϵ = 1e-4

    θ = 170.0
    ϕ = 5.0

    n_skip = 1

    function regular()

        atmos = Atmosphere(get_atmos(DATA; periodic=true, skip=n_skip)...)

        @time J_mean, S_λ, α_tot = Λ_regular(ϵ, maxiter, atmos, quadrature)

        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        I_top = short_characteristics_up(k, S_λ, α_tot, atmos, I_0=S_λ[1,:,:])

        I_top = ustrip(uconvert.(u"kW*nm^-1*m^-2", I_top[end, 2:end-1, 2:end-1]))

        global min_lim, max_lim
        min_lim = minimum(I_top)
        max_lim = maximum(I_top)

        heatmap(ustrip.(atmos.x[2:end-1]),
                ustrip.(atmos.y[2:end-1]),
                transpose.(I_top),
                xaxis="x",
                yaxis="y",
                dpi=300,
                rightmargin=10Plots.mm,
                title="Regular Grid",
                aspect_ratio=:equal,
                clim=(min_lim,max_lim))

        μ = abs(k[1])
        savefig("../img/compare_continuum/regular_cont_$(floor(Int, 100*μ))")

        return 0
    end

    function voronoi()

        atmos = Atmosphere(get_atmos(DATA; periodic=false, skip=n_skip)...)

        nx = length(atmos.x)
        nz = length(atmos.z)
        ny = length(atmos.y)

        x_min = ustrip(atmos.x[1])
        x_max = ustrip(atmos.x[end])
        y_min = ustrip(atmos.y[1])
        y_max = ustrip(atmos.y[end])
        z_min = ustrip(atmos.z[1])
        z_max = ustrip(atmos.z[end])

        n_sites = floor(Int, nz*nx*ny)
        positions = rejection_sampling(n_sites, atmos, log10.(ustrip.(atmos.hydrogen_populations)))

        # rand(3, n_sites)

        # positions[1, :] = positions[1, :].*(z_max - z_min) .+ z_min
        # positions[2, :] = positions[2, :].*(x_max - x_min) .+ x_min
        # positions[3, :] = positions[3, :].*(y_max - y_min) .+ y_min

        # positions = positions*1u"m"

        # rejection_sampling(n_sites, atmos, log10.(ustrip.(atmos.hydrogen_populations)))
        # sample_from_extinction(atmos, 500.0u"nm", n_sites)

        sites_file = "../data/sites_continuum_half.txt"
        neighbours_file = "../data/neighbours_continuum_half.txt"
        # write sites to file
        write_arrays(ustrip.(positions[2, :]),
                     ustrip.(positions[3, :]),
                     ustrip.(positions[1, :]),
                     sites_file)

        # export sites to voro++, and compute grid information
        println("---Preprocessing grid---")

        # compute neigbours
        # run(`./voro.sh $sites_file $neighbours_file
                       # $(x_min) $(x_max)
                       # $(y_min) $(y_max)
                       # $(z_min) $(z_max)`)

        # Voronoi grid
        sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                             _initialise(positions, atmos)...,
                             z_min*1u"m", z_max*1u"m",
                             x_min*1u"m", x_max*1u"m",
                             y_min*1u"m", y_max*1u"m",
                             n_sites)

        # plot_sites(sites)

        @time J_mean, S_λ, α_tot = Λ_voronoi(ϵ, maxiter, sites, quadrature)

        atmos_from_voronoi, S_λ_grid, α_grid = Voronoi_to_Raster(sites, atmos,
                                                                 S_λ, α_tot, 2;
                                                                 periodic=true)

        k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
        I_top = short_characteristics_up(k, S_λ_grid, α_grid,
                                        atmos_from_voronoi, I_0=S_λ_grid[1,:,:])

        I_top = ustrip(uconvert.(u"kW*nm^-1*m^-2", I_top[end, 2:end-1, 2:end-1]))

        heatmap(ustrip.(atmos_from_voronoi.x[2:end-1]),
                ustrip.(atmos_from_voronoi.y[2:end-1]),
                transpose.(I_top),
                xaxis="x",
                yaxis="y",
                dpi=300,
                rightmargin=10Plots.mm,
                title="Irregular Grid",
                aspect_ratio=:equal,
                clim=(min_lim,max_lim))

        μ = abs(k[1])
        savefig("../img/compare_continuum/irregular_cont_$(floor(Int, 100*μ))")

        return 0
    end

    regular();
    # voronoi();

end

function LTE_compare(DATA::String, n_sites::Int)

    # atmos = Atmosphere(get_atmos(DATA; periodic=true, skip=1)...)

    # choose a wavelength
    λ = 500u"nm"  # nm

    """
    # Lte populations
    LTE_pops = LTE_populations(atmos)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])

    # Planck function
    S_λ = blackbody_λ.(λ, atmos.temperature)

    θ = 180.0
    ϕ = 0.0

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    intensity = short_characteristics_up(k, S_λ, S_λ[1,:,:], α_cont, atmos)

    I_regular = transpose.(ustrip.(uconvert.(u"kW*nm^-1*m^-2", intensity[end, 2:end-1, 2:end-1])))
    x_regular = atmos.x[2:end-1]
    y_regular = atmos.y[2:end-1]
    """

    atmos = Atmosphere(get_atmos(DATA; periodic=false, skip=1)...)

    nx = length(atmos.x)
    nz = length(atmos.z)
    ny = length(atmos.y)

    x_min = ustrip(atmos.x[1])
    x_max = ustrip(atmos.x[end])
    y_min = ustrip(atmos.y[1])
    y_max = ustrip(atmos.y[end])
    z_min = ustrip(atmos.z[1])
    z_max = ustrip(atmos.z[end])

    #positions = rand(3, n_sites)

    #positions[1, :] = positions[1, :].*(z_max - z_min) .+ z_min
    #positions[2, :] = positions[2, :].*(x_max - x_min) .+ x_min
    #positions[3, :] = positions[3, :].*(y_max - y_min) .+ y_min

    #positions = positions*1u"m"

    positions = sample_from_ionised_hydrogen(atmos, n_sites)
    # sample_from_extinction(atmos, λ, n_sites)
    # positions = sample_from_destruction(atmos, n_sites)
    # positions = sample_from_temp_gradient(atmos, n_sites)
    # positions = rejection_sampling(n_sites, atmos, ustrip.(atmos.temperature))

    sites_file = "../data/sites_continuum_$n_sites.txt"
    neighbours_file = "../data/neighbours_continuum_$n_sites.txt"
    # write sites to file
    write_arrays(ustrip.(positions[2, :]),
                 ustrip.(positions[3, :]),
                 ustrip.(positions[1, :]),
                 sites_file)

    # export sites to voro++, and compute grid information
    println("---Preprocessing grid---")

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
                   $(x_min) $(x_max)
                   $(y_min) $(y_max)
                   $(z_min) $(z_max)`)

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         _initialise(positions, atmos)...,
                         z_min*1u"m", z_max*1u"m",
                         x_min*1u"m", x_max*1u"m",
                         y_min*1u"m", y_max*1u"m",
                         n_sites)

    run(`rm ../data/$sites_file ../data/$neighbours_file`)

    # Lte populations
    LTE_pops = LTE_populations(sites)

    S_λ = Matrix{Float64}(undef, (1, n_sites))u"kW*m^-2*nm^-1"
    S_λ[1,:] = blackbody_λ.(λ, sites.temperature)

    # println(typeof(S_λ))
    # println(typeof(LTE_pops))

    atmos_size = (nz, nx, ny).*2
    atmos_size = floor.(Int, atmos_size)

    atmos, S_λ_grid, populations_grid = Voronoi_to_Raster(sites, atmos_size,
                                                          S_λ, LTE_pops;
                                                          periodic=true)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           populations_grid[:,:,:,1].+populations_grid[:,:,:,2],
                           populations_grid[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           populations_grid[:,:,:,1])

    θ = 180.0
    ϕ = 0.0

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]

    I_0 = S_λ_grid[1,1,:,:]
    intensity = short_characteristics_up(k, S_λ_grid[1,:,:,:], I_0, α_cont, atmos; n_sweeps=3)

    I_irregular = transpose.(ustrip.(uconvert.(u"kW*nm^-1*m^-2", intensity[end, 2:end-1, 2:end-1])))
    x_irregular = atmos.x[2:end-1]
    y_irregular = atmos.y[2:end-1]

    npzwrite("../data/LTE/I_irregular_$(n_sites)_ionised_hydrogen.npy", I_irregular)

end

function LTE_regular(DATA::String, n_skip::Integer)

    if n_skip == 4
        RES = "quarter"
    elseif n_skip == 3
        RES = "third"
    elseif n_skip == 2
        RES = "half"
    elseif n_skip == 1
        RES = "full"
    else
        RES = string(n_skip)
    end

    println("Continuum at $(RES) resolution, regular grid")

    atmos = Atmosphere(get_atmos(DATA; periodic=true, skip=n_skip)...)

    # choose a wavelength
    λ = 500u"nm"  # nm

    # Lte populations
    LTE_pops = LTE_populations(atmos)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])

    # Planck function
    S_λ = blackbody_λ.(λ, atmos.temperature)

    θ = 180.0
    ϕ = 0.0

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    intensity = short_characteristics_up(k, S_λ, S_λ[1,:,:], α_cont, atmos)

    I_regular = transpose.(ustrip.(uconvert.(u"kW*nm^-1*m^-2", intensity[end, 2:end-1, 2:end-1])))
    x_regular = atmos.x[2:end-1]
    y_regular = atmos.y[2:end-1]

    npzwrite("../data/LTE/I_regular_$(RES).npy", I_regular)
    npzwrite("../data/LTE/y_regular_$(RES).npy", ustrip.(x_regular))
    npzwrite("../data/LTE/x_regular_$(RES).npy", ustrip.(y_regular))

end

function test_interpolation(DATA)
    atmos = Atmosphere(get_atmos(DATA; periodic=false, skip=1)...)

    nx = length(atmos.x)
    nz = length(atmos.z)
    ny = length(atmos.y)

    x_min = ustrip(atmos.x[1])
    x_max = ustrip(atmos.x[end])
    y_min = ustrip(atmos.y[1])
    y_max = ustrip(atmos.y[end])
    z_min = ustrip(atmos.z[1])
    z_max = ustrip(atmos.z[end])

    n_sites = floor(Int, nz*nx*ny)
    positions = rand(3, n_sites)

    positions[1, :] = positions[1, :].*(z_max - z_min) .+ z_min
    positions[2, :] = positions[2, :].*(x_max - x_min) .+ x_min
    positions[3, :] = positions[3, :].*(y_max - y_min) .+ y_min

    positions = positions*1u"m"

    # rejection_sampling(n_sites, atmos, log10.(ustrip.(atmos.hydrogen_populations)))
    # sample_from_extinction(atmos, 500.0u"nm", n_sites)

    sites_file = "../data/sites_continuum.txt"
    neighbours_file = "../data/neighbours_continuum.txt"
    # write sites to file
    write_arrays(ustrip.(positions[2, :]),
                 ustrip.(positions[3, :]),
                 ustrip.(positions[1, :]),
                 sites_file)

    # export sites to voro++, and compute grid information
    println("---Preprocessing grid---")

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
                   $(x_min - 0.1) $(x_max + 0.1)
                   $(y_min - 0.1) $(y_max + 0.1)
                   $(z_min - 0.1) $(z_max + 0.1)`)

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         _initialise(positions, atmos)...,
                         z_min*1u"m", z_max*1u"m",
                         x_min*1u"m", x_max*1u"m",
                         y_min*1u"m", y_max*1u"m",
                         n_sites)

    # Lte populations
    LTE_pops = LTE_populations(sites)
    λ = 500u"nm"
    # Find continuum extinction
    α_cont = α_absorption.(λ,
                           sites.temperature,
                           sites.electron_density*1.0,
                           LTE_pops[:,1].+LTE_pops[:,2],
                           LTE_pops[:,3]) .+
             α_scattering.(λ,
                           sites.electron_density,
                           LTE_pops[:,1])

    S_λ = blackbody_λ.(λ, sites.temperature)

    θ = 175.0
    ϕ = 0.1

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    bottom_layer = sites.layers_up[2] - 1
    bottom_layer_idx = sites.perm_up[1:bottom_layer]
    println("---Ray tracing---")
    I_0 = blackbody_λ.(500u"nm", sites.temperature[bottom_layer_idx])
    intensity = Delaunay_upII(k, S_λ, α_cont, sites, I_0, 3)

    x = collect(LinRange(sites.x_min, sites.x_max, 10*nx))
    y = collect(LinRange(sites.y_min, sites.y_max, 10*ny))
    top_z = atmos.z[end]

    tree = KDTree(ustrip(sites.positions))

    I_top = Matrix{Float64}(undef, (length(x), length(y)))u"kW*m^-2*nm^-1"
    for i in 1:length(x)
        for j in 1:length(y)
            position = [top_z, x[i], y[j]]
            idx, dist = nn(tree, ustrip(position))
            I_top[i, j] = intensity[idx]
        end
    end

    heatmap(ustrip.(x),
            ustrip.(y),
            ustrip.(transpose.(I_top)),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="LTE Intensity (500 nm), Irregular Grid",
            aspect_ratio=:equal)

    savefig("../img/compare_continuum/irregular_top_I")

    atmos_from_voronoi = Voronoi_to_Raster(sites, atmos, 1; periodic=false)

    T_diff = abs.(1 .- atmos_from_voronoi.temperature ./ atmos.temperature)
    N_e_diff = abs.(1 .- atmos_from_voronoi.electron_density ./ atmos.electron_density)
    N_H_diff = abs.(1 .- atmos_from_voronoi.hydrogen_populations ./ atmos.hydrogen_populations)

    heatmap(ustrip.(atmos.x),
            ustrip.(atmos.y),
            ustrip.(transpose.(T_diff[end-20,:,:])),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="Temperature Difference",
            aspect_ratio=:equal)

    savefig("../img/compare_continuum/T_diff")

    heatmap(ustrip.(atmos.x),
            ustrip.(atmos.y),
            ustrip.(transpose.(N_e_diff[end-20,:,:])),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="Electron Density Difference",
            aspect_ratio=:equal)

    savefig("../img/compare_continuum/N_e_diff")

    heatmap(ustrip.(atmos.x),
            ustrip.(atmos.y),
            ustrip.(transpose.(N_H_diff[end-20,:,:])),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="Hydrogen Density Difference",
            aspect_ratio=:equal)

    savefig("../img/compare_continuum/N_H_diff")

    return 0
end

function test_with_regular_grid(DATA, quadrature)
    atmos = Atmosphere(get_atmos(DATA; periodic=false, skip=1)...)

    nx = length(atmos.x)
    nz = length(atmos.z)
    ny = length(atmos.y)

    x_min = ustrip(atmos.x[1])
    x_max = ustrip(atmos.x[end])
    y_min = ustrip(atmos.y[1])
    y_max = ustrip(atmos.y[end])
    z_min = ustrip(atmos.z[1])
    z_max = ustrip(atmos.z[end])

    n_sites = floor(Int, nz*nx*ny)
    positions = Matrix{Float64}(undef, (3, n_sites))u"m"
    temperature = Vector{Float64}(undef, n_sites)u"K"
    hydrogen_populations = Vector{Float64}(undef, n_sites)u"m^-3"
    electron_density = Vector{Float64}(undef, n_sites)u"m^-3"
    velocity_z = Vector{Float64}(undef, n_sites)u"m/s"
    velocity_x = Vector{Float64}(undef, n_sites)u"m/s"
    velocity_y = Vector{Float64}(undef, n_sites)u"m/s"

    for k in eachindex(atmos.z)
        for i in eachindex(atmos.x)
            for j in eachindex(atmos.y)
                p = [atmos.z[k], atmos.x[i], atmos.y[j]]
                m = k + (i-1)*nz + (j-1)*nz*nx
                positions[:,m] = p
                temperature[m] = atmos.temperature[k, i, j]
                hydrogen_populations[m] = atmos.hydrogen_populations[k, i, j]
                electron_density[m] = atmos.electron_density[k, i, j]
                velocity_z[m] = atmos.velocity_z[k, i, j]
                velocity_x[m] = atmos.velocity_x[k, i, j]
                velocity_y[m] = atmos.velocity_y[k, i, j]
            end
        end
    end

    sites_file = "../data/sites_continuum.txt"
    neighbours_file = "../data/neighbours_continuum.txt"
    # write sites to file
    write_arrays(ustrip.(positions[2, :]),
                 ustrip.(positions[3, :]),
                 ustrip.(positions[1, :]),
                 sites_file)

    # export sites to voro++, and compute grid information
    println("---Preprocessing grid---")

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
                   $(x_min-0.1) $(x_max+0.1)
                   $(y_min-0.1) $(y_max+0.1)
                   $(z_min-0.1) $(z_max+0.1)`)

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         temperature, electron_density, hydrogen_populations,
                         velocity_z, velocity_x, velocity_y,
                         z_min*1u"m", z_max*1u"m",
                         x_min*1u"m", x_max*1u"m",
                         y_min*1u"m", y_max*1u"m",
                         n_sites)

    # Lte populations
    LTE_pops = LTE_populations(sites)
    λ = 500u"nm"
    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           LTE_pops[:,:,:,1].+LTE_pops[:,:,:,2],
                           LTE_pops[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           LTE_pops[:,:,:,1])
    S_λ = blackbody_λ.(λ, sites.temperature)

    θ = 170.0
    ϕ = 30.0
    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]
    bottom_layer = sites.layers_up[2] - 1
    bottom_layer_idx = sites.perm_up[1:bottom_layer]
    I_0 = blackbody_λ.(500u"nm", sites.temperature[bottom_layer_idx])
    intensity = Delaunay_upII(k, S_λ, α_cont, sites, I_0, 3)

    x = collect(LinRange(sites.x_min, sites.x_max, 5*nx))
    y = collect(LinRange(sites.y_min, sites.y_max, 5*ny))
    top_z = sites.z_max

    tree = KDTree(ustrip(sites.positions))

    I_top = Matrix{Float64}(undef, (length(x), length(y)))u"kW*m^-2*nm^-1"
    for i in 1:length(x)
        for j in 1:length(y)
            position = [top_z, x[i], y[j]]
            idx, dist = nn(tree, ustrip(position))
            I_top[i, j] = intensity[idx]
        end
    end

    heatmap(ustrip.(x),
            ustrip.(y),
            ustrip.(transpose.(I_top)),
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="intensity",
            aspect_ratio=:equal,
            clim=(19.0, maximum(ustrip.(I_top))))

    savefig("../img/compare_continuum/irregular_top_I")

    return 0
end

function compare_interpolations(DATA::String, n_sites::Int)

    atmos = Atmosphere(get_atmos(DATA; periodic=false, skip=1)...)

    # choose a wavelength
    λ = 500u"nm"  # nm

    nx = length(atmos.x)
    nz = length(atmos.z)
    ny = length(atmos.y)

    x_min = ustrip(atmos.x[1])
    x_max = ustrip(atmos.x[end])
    y_min = ustrip(atmos.y[1])
    y_max = ustrip(atmos.y[end])
    z_min = ustrip(atmos.z[1])
    z_max = ustrip(atmos.z[end])

    positions = sample_from_extinction(atmos, λ, n_sites)
    # positions = sample_from_temp_gradient(atmos, n_sites)
    # positions = rejection_sampling(n_sites, atmos, ustrip.(atmos.temperature))

    sites_file = "../data/sites_write.txt"
    neighbours_file = "../data/neighbours_write.txt"
    # write sites to file
    write_arrays(ustrip.(positions[2, :]),
                 ustrip.(positions[3, :]),
                 ustrip.(positions[1, :]),
                 sites_file)

    # export sites to voro++, and compute grid information
    println("---Preprocessing grid---")

    # compute neigbours
    run(`./voro.sh $sites_file $neighbours_file
                   $(x_min) $(x_max)
                   $(y_min) $(y_max)
                   $(z_min) $(z_max)`)

    # Voronoi grid
    sites = VoronoiSites(read_cell(neighbours_file, n_sites, positions)...,
                         _initialise(positions, atmos)...,
                         z_min*1u"m", z_max*1u"m",
                         x_min*1u"m", x_max*1u"m",
                         y_min*1u"m", y_max*1u"m",
                         n_sites)

    run(`rm ../data/$sites_file ../data/$neighbours_file`)

    # Lte populations
    populations = LTE_populations(sites)

    # LTE source function --> Planck function
    S_λ = Matrix{Float64}(undef, (1, n_sites))u"kW*m^-2*nm^-1"
    S_λ[1,:] = blackbody_λ.(λ, sites.temperature)

    CMAX = 60.52755753644262
    CMIN = 20.169713390042126

    atmos_size = (430, 256, 256)

    θ = 180.0
    ϕ = 0.0

    k = [cos(θ*π/180), cos(ϕ*π/180)*sin(θ*π/180), sin(ϕ*π/180)*sin(θ*π/180)]

    # Inverse distance interpolation
    atmos, S_λ_grid, populations_grid = Voronoi_to_Raster_inv_dist(sites, atmos_size, S_λ, populations)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           populations_grid[:,:,:,1].+populations_grid[:,:,:,2],
                           populations_grid[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           populations_grid[:,:,:,1])

    I_0 = S_λ_grid[1,1,:,:]
    intensity = short_characteristics_up(k, S_λ_grid[1,:,:,:], I_0, α_cont, atmos; n_sweeps=3)

    I_top = transpose.(ustrip.(uconvert.(u"kW*nm^-1*m^-2", intensity[end, 2:end-1, 2:end-1])))
    x = ustrip.(atmos.x[2:end-1])
    y = ustrip.(atmos.y[2:end-1])

    heatmap(x,
            y,
            I_top,
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="Inverse Distance interpolation",
            aspect_ratio=:equal,
            clim=(CMIN,CMAX))

    savefig("../img/compare_continuum/inv_dist")

    # Nearest neighbour interpolation
    atmos, S_λ_grid, populations_grid = Voronoi_to_Raster(sites, atmos_size, S_λ, populations)

    # Find continuum extinction (only with Thomson and Rayleigh)
    α_cont = α_absorption.(λ,
                           atmos.temperature,
                           atmos.electron_density*1.0,
                           populations_grid[:,:,:,1].+populations_grid[:,:,:,2],
                           populations_grid[:,:,:,3]) .+
             α_scattering.(λ,
                           atmos.electron_density,
                           populations_grid[:,:,:,1])

    I_0 = S_λ_grid[1,1,:,:]
    intensity = short_characteristics_up(k, S_λ_grid[1,:,:,:], I_0, α_cont, atmos; n_sweeps=3)

    I_top = transpose.(ustrip.(uconvert.(u"kW*nm^-1*m^-2", intensity[end, 2:end-1, 2:end-1])))
    x = ustrip.(atmos.x[2:end-1])
    y = ustrip.(atmos.y[2:end-1])

    heatmap(x,
            y,
            I_top,
            xaxis="x",
            yaxis="y",
            dpi=300,
            rightmargin=10Plots.mm,
            title="Nearest Neighbour interpolation",
            aspect_ratio=:equal,
            clim=(CMIN,CMAX))

    savefig("../img/compare_continuum/nearest_neighbour")

end

# compare("../data/bifrost_qs006023_s525_half.hdf5", "../quadratures/ul2n3.dat");
n_list = [100_000, 250_000, 500_000, 1_000_000, 2_500_000, 5_000_000, 10_000_000, 15_000_000]

for num_sites in n_list
    LTE_compare("../data/bifrost_qs006023_s525.hdf5", num_sites)
end
# LTE_regular("../data/bifrost_qs006023_s525_quarter.hdf5")
# LTE_regular("../data/bifrost_qs006023_s525_half.hdf5")
# LTE_regular("../data/bifrost_qs006023_s525.hdf5", 1)
# LTE_regular("../data/bifrost_qs006023_s525.hdf5", 2)
# LTE_regular("../data/bifrost_qs006023_s525.hdf5", 3)
# LTE_regular("../data/bifrost_qs006023_s525.hdf5", 4)
# test_interpolation("../data/bifrost_qs006023_s525_half.hdf5", "../quadratures/n1.dat")
# test_with_regular("../data/bifrost_qs006023_s525_quarter.hdf5", "../quadratures/n1.dat")
# write_grid("../data/bifrost_qs006023_s525.hdf5", 10_000_000)
# compare_interpolations("../data/bifrost_qs006023_s525.hdf5", 1_000_000)
print("")

max_stars_repo_name: meudnaes/VoronoiRT

max_issues_repo_name: meudnaes/VoronoiRT

max_forks_repo_name: meudnaes/VoronoiRT


text:
using OmegaCore
using Distributions
using OmegaDistributions
using Test

@testset "Pointwie" begin
  x = 1 ~ Normal(0, 1)
  y = abs.(x)
  y_, x_ = randsample((y, x))
  @test y_ == abs(x_)
  @test isinferred(randsample, x)
end

@testset "Pointwise class" begin
  μ = @~ Normal(0, 1)
  ϵs = ~ Normal(0, 1)
  ys = μ .+ ϵs
  y1 = 1 ~ ys
  y2 = 2 ~ ys
  randsample((y1, y2))
end

@testset "Pointwise func" begin
  issin = @~ Bernoulli(0.5)
  f = ifelse.(issin, sin, sqrt)

  μ = @~ Normal(0, 1)
  ϵs = ~ Normal(0, 1)
  ys = μ .+ ϵs
  y1 = 1 ~ ys
  y2 = 2 ~ ys
  randsample((y1, y2))
end

function f()
  y = @~ Normal(0, 1)
  z = y .+ y
  randsample(z)
  ω = defω()
  z(ω)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using Franklin

serve(single=true, clear=true)
publish(minify=false)

max_stars_repo_name: jbmuir/JackBMuir

max_issues_repo_name: jbmuir/JackBMuir

max_forks_repo_name: jbmuir/JackBMuir


text:
module IVPTestSuite
# Differential equation initial value problem test suite

import Base: +, -, *, /, .+, .-, .*, ./,==,>,>=

# to dispose of the @doc macro:
if VERSION<v"0.4"
    macro doc(a...)
        if length(a)>1
            esc(a[end])
        else
            esc(a[1].args[end])
        end
    end
end

using Requires

# modules
export Solvers
# variables
export nonstiff, mildlystiff, stiff, explicit_eq, explicit_mass_eq, implicit_eq #, imex_eq
# types
export Stiff, Solver, TestCase, TestCaseExplicit, TestRun, TestSuite, TestResults, @doc,
       TestRunAdapt, TestRunFixedStep
# functions
export isadaptive, isimplicit, isapplicable, isconsistent
export hasjacobian, hasmass, solvertype, name

# some enumerations
abstract Enum
immutable Stiff <: Enum
    val::Int
    function Stiff(i)
        @assert isa(i, Integer)
        @assert 0<=i<=2
        new(i)
    end
end
const nonstiff = Stiff(0)
const mildlystiff = Stiff(1)
const stiff = Stiff(2)

immutable EqType <: Enum # type of equation a solver can handle
    val::Int
    function EqType(i)
        @assert isa(i, Integer)
        @assert 0<=i<=3
        new(i)
    end
end
const explicit_eq = EqType(0)            # dy/dt = F(t, y)
const explicit_mass_eq = EqType(1)       # M(t,y)*dy/dt = F(t, y)
const implicit_eq = EqType(2)            # 0 = F(t, y, dy/dt)
# const imex_eq = EqType(3)                # G(t,y) = F(t, y, dy/dt)

# common enumeration functions
==(s1::Enum, s2::Enum) = s1.val==s2.val
>(s1::Enum, s2::Enum) = s1.val>s2.val
>=(s1::Enum, s2::Enum) = s1.val>=s2.val

# Holds all the constants for a certain test
# TODO: think about also implementing implicitly defined and
# PETSc-style problems
abstract TestCase{Name, T<:Real, Tarr<:AbstractMatrix}
immutable TestCaseExplicit{Name, T, Tarr} <: TestCase{Name,T,Tarr}  # for explicit ODE/DAE equation
    stiff::Stiff          # how stiff the problem is: nonstiff, mildlystiff, stiff
    dae::Int              # index of DAE, ==0 for ODE
    dof::Int              # degrees of freedom
    fn!::Function          # objective function
    jac!::Union{Function,Void}  # Jacobian of fn, can also initialize an array for the Jacobian
    mass!::Union{Function,Void} # mass matrix, can also initialize an array for the mass matrix
    ic::Vector{T}         # Initial condition
    tspan::Vector{T}      # [t_start, t_end]
    refsol::Vector{T}     # reference solution at t=t_end
    refsolinds::Union{BitVector,Vector{Int}}  # which indices of the solution should compared to refsol
    scd_absinds::Union{BitVector,Vector{Int}}   # components for which to use the absolute
                                                # error instead of relative, c.f. calc_error_scd
end
# TODO make:
# immutable TestCaseImplicit{Name, T} <: TestCase  # for implicit ODE/DAE equation

name{Name}(tc::TestCase{Name}) = Name
hasjacobian(tc::TestCase) = (tc.jac!)==nothing ? false : true
Base.issparse{N,T,Tarr}(tc::TestCase{N,T,Tarr}) = Tarr<:AbstractSparseMatrix
hasmass(tc::TestCaseExplicit) = (tc.mass!)==nothing ? false : true
function Base.show(io::IO, tc::TestCase)
    dae = tc.dae>0 ? "DAE" : "ODE"
    st = tc.stiff> nonstiff ? "stiff" : "non-stiff"
    print(io, "Test case \"$(name(tc))\" ($dae, DOF=$(tc.dof), $st)")
end

@doc "Holds information about a solver" ->
immutable Solver{Typ}    # Typ: :ex (explicit method), :im (implicit method), :imex (IMEX method)
    solverfn::Function   # the actual solver
    package::Module      # the package-module where it is defined
    wrapper::Function    # a function (tr::TestRun) ->
    # capabilities:
    stiff::Stiff         # stiff, mildly stiff or non-stiff solver
    adaptive::Bool
    dae::Int             # maximum index of DAE it can solve, ==0 for ODE only
    eq_type::EqType      # dy/dt=F(t,y); M(t,y)*dy/dt=F(t,y); F(t,y,dy/dt)=0
end
isadaptive(s::Solver) = s.adaptive
solvertype{Typ}(s::Solver{Typ}) = Typ
isimplicit(s::Solver) = (solvertype(s)==:im || solvertype(s)==:imex) ? true : false # i.e. uses an implicit method
name(s::Solver) = s.solverfn
function Base.show(io::IO, s::Solver)
    ad = isadaptive(s) ? "adaptive" : "fixed step"
    im = isimplicit(s) ? "implicit" : "explicit"
    dae = s.dae>0 ? "DAE" : "ODE"
    if s.stiff==nonstiff
        st = "non-stiff"
    elseif s.stiff==mildlystiff
        st = "mildly-stiff"
    else
        st = "stiff"
    end
    print(io, "Solver $(name(s)) ($ad, $dae, $st, $im)")
end

# Holds a few extra bits for a specific test run.
abstract TestRun{Name, T<:Real}
immutable TestRunAdapt{Name, T} <: TestRun{Name, T}
    tc::TestCase{Name,T}
    solver::Solver          # which solver to use
    solverpara::Dict{Symbol,Any}  # extra parameters a solver needs for
                                 # a particular TestCase
    abstol::T               # absolute error tolerance
    reltol::T               # relative error tolerance
    h0::T                   # first time step size
end
isadaptive(::TestRunAdapt) = true
_show(tr::TestRunAdapt) = "an adaptive solver using abstol = $(tr.abstol) and reltol = $(tr.reltol)"
immutable TestRunFixedStep{Name, T} <: TestRun{Name, T}
    tc::TestCase{Name,T}
    solver::Solver          # which solver to use
    solverpara::Dict{Symbol,Any} # extra parameters a solver needs for
                                 # a particular TestCase
    tsteps::Vector{T}       # time steps vector
end
isadaptive(::TestRunFixedStep) = false
_show(tr::TestRunFixedStep) = "a fixed step solver using $(length(tr.tsteps)) time steps"

immutable TestSuite{Name, T<:Real}
    tc::TestCase{Name, T}
    solver::Solver                   # which solver to use
    solverpara::Dict{Symbol,Any}     # extra parameters a solver needs for
                                     # a particular TestCase
    abstols::Vector{T}               # absolute error tolerance
    reltols::Vector{T}               # relative error tolerance
    h0s::Vector{T}                   # first time step size
    tstepss::Vector{Vector{T}}       # for non-adaptive steppers set to the steps
    nruns::Int
    adapt::Bool
    function TestSuite(tc, solver, abstols, reltols, h0s; solverpara=Dict{Symbol,Any}())
        # constructor for adaptive solvers
        if !isadaptive(solver)
            error("Constructor `TestSuite{Name,T}(tc::TestCase{Name,T}, solver, abstols, reltols, h0s)` is for adaptive solvers only")
        end
        len = 1
        for v in Any[abstols, reltols, h0s] # then need to be length==1 or same length
            if length(v)>1
                if len>1 && length(v)!=len
                    error("abstols, reltols, and h0s need to be length==1 or same length")
                else
                    len = length(v)
                end
            end
        end
        new(tc, solver, solverpara, abstols, reltols, h0s, Vector{T}[], len, true)
    end
    function TestSuite(tc, solver, tstepss; solverpara=Dict{Symbol,Any}())
        # constructor for fixed step solvers
        if isadaptive(solver)
            error("Constructor `TestSuite{Name,T}(tc::TestCase{Name,T}, solver, tstepss)` is for fixed step solvers only.")
        end
        len = length(tstepss)
        new(tc, solver, solverpara, T[], T[], T[], tstepss, len, false)
    end

end
TestSuite{Name,T}(tc::TestCase{Name,T}, solver, abstols, reltols, h0s; solverpara=Dict{Symbol,Any}() ) = TestSuite{Name,T}(tc, solver, abstols, reltols, h0s; solverpara=solverpara )
TestSuite{Name,T}(tc::TestCase{Name,T}, solver, tstepss; solverpara=Dict{Symbol,Any}() ) = TestSuite{Name,T}(tc, solver, tstepss; solverpara=solverpara )
isadaptive(ts::TestSuite) = ts.adapt

Base.length(ts::TestSuite) = ts.nruns
# make iterator which returns a TestRun for each member in TestSuite
Base.start(ts::TestSuite) = 1
Base.done(ts::TestSuite, state) = length(ts)+1==state ? true : false
function Base.next(ts::TestSuite, state)
    if isadaptive(ts)
        abstol = length(ts.abstols)>1 ? ts.abstols[state] : ts.abstols[1]
        reltol = length(ts.reltols)>1 ? ts.reltols[state] : ts.reltols[1]
        h0 = length(ts.h0s)>1 ? ts.h0s[state] : ts.h0s[1]
        return TestRunAdapt(ts.tc, ts.solver, ts.solverpara, abstol, reltol, h0), state+1
    else
        tsteps = ts.tstepss[state]
        if tsteps[1]!=ts.tc.tspan[1] ||  tsteps[end]!=ts.tc.tspan[2]
            error("Time interval of tsteps not equal to tspan.")
        end
        return TestRunFixedStep(ts.tc, ts.solver, ts.solverpara, tsteps), state+1
    end
end

@doc "Hold the results from a test run" ->
immutable TestResults{Name,T<:Real}
    testrun::TestRun{Name}
    tend::T   # the returned end-time
    yend::Vector{T}     # the calculated solution.  Needs to be initialized to Array(T,dof)
    steps_total::Int
    steps_accepted::Int
    fn_evals::Int
    jac_evals::Int
    linear_solves::Int
    walltime::Float64 # (s)
    mem::Int # memory allocated (bytes)
    gc_time::Float64 # time used by garbage collode23sector
    scd::Float64  # Error estimate: significant digits, see Mazzia & Magherini p.II-ii
    mescd::Float64  # a variation on scd, see Mazzia & Magherini p.II-ii
    error::Union{Void, Exception} # holds the exception if one occured
end
function Base.show( io::IO, res::TestResults)
    tr = res.testrun
    tc = res.testrun.tc
    so = res.testrun.solver
    scd = round(res.scd, 3)
    wt = round(res.walltime*1000, 3)
    mb = round(res.mem/1e6, 3)
    trdetails = _show(tr)
    print(io,
          """
          Results of TestCase \"$(name(tc))\" solved with $(so.solverfn),
          $trdetails

          Significant digits: $(scd)
          Walltime:           $wt ms
          Memory allocated:   $mb MB
          """)

end

# misc helper functions
name{Name}(t::Union{TestRun{Name}, TestResults{Name}, TestCase{Name}}) = Name

## misc helper funs
include("helpers.jl")

## solvers functions
include("solvers/solvers.jl")

## Running tests
include("running_tests.jl")

## Evaluating tests
include("eval_tests.jl")

## The test cases
include("testcases/testcases.jl")

## Plotting
@require PyPlot include(Pkg.dir()*"/IVPTestSuite/src/plotting.jl") # bug requires full path: https://github.com/one-more-minute/Requires.jl/issues/2

end # module


text:
using DashBootstrapComponents

spinners = html_div([
    dbc_spinner(size = "sm"),
    html_hr(),
    dbc_spinner(spinner_style = Dict("width" => "3rem", "height" => "3rem")),
]);

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
"""
    struct LangOptions <: Any
Hold a pointer to a `clang::LangOptions` object.
"""
struct LangOptions
    ptr::CXLangOptions
end

Base.unsafe_convert(::Type{CXLangOptions}, x::LangOptions) = x.ptr
Base.cconvert(::Type{CXLangOptions}, x::LangOptions) = x


text:
lm = LakeModel()
e_0 = 0.92     # Initial employment rate
u_0 = 1 - e_0  # Initial unemployment rate
T = 50         # Simulation length

xbar = rate_steady_state(lm)
x_0 = [u_0; e_0]
x_path = simulate_rate_path(lm, x_0, T)

titles = ["Unmployment rate" "Employment rate"]

fig, axes = subplots(2, 1, figsize=(10, 8))

for (i, ax) in enumerate(axes)
    ax[:plot](1:T, x_path[i, :], c="blue", lw=2, alpha=0.5)
    ax[:hlines](xbar[i], 0, T, "r", "--")
    ax[:set](title=titles[i])
    ax[:grid]("on")
end


text:
# ---
# title: 191. Number of 1 Bits
# id: problem191
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/number-of-1-bits/description/>
# hidden: true
# ---
#
# Write a function that takes an unsigned integer and returns the number of '1'
# bits it has (also known as the [Hamming
# weight](http://en.wikipedia.org/wiki/Hamming_weight)).
#
# **Note:**
#
#   * Note that in some languages such as Java, there is no unsigned integer type. In this case, the input will be given as a signed integer type. It should not affect your implementation, as the integer's internal binary representation is the same, whether it is signed or unsigned.
#   * In Java, the compiler represents the signed integers using [2's complement notation](https://en.wikipedia.org/wiki/Two%27s_complement). Therefore, in **Example 3**  above, the input represents the signed integer. `-3`.
#
# **Follow up** : If this function is called many times, how would you optimize
# it?
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 00000000000000000000000000001011
#     Output: 3
#     Explanation: The input binary string **00000000000000000000000000001011** has a total of three '1' bits.
#
#
# **Example 2:**
#
#
#
#     Input: n = 00000000000000000000000010000000
#     Output: 1
#     Explanation: The input binary string **00000000000000000000000010000000** has a total of one '1' bit.
#
#
# **Example 3:**
#
#
#
#     Input: n = 11111111111111111111111111111101
#     Output: 31
#     Explanation: The input binary string **11111111111111111111111111111101** has a total of thirty one '1' bits.
#
#
#
#
# **Constraints:**
#
#   * The input must be a **binary string** of length `32`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
const mac_gdrive = "/Volumes/GoogleDrive/My Drive"
const linux_gdrive = "/home/beatriz/mainen.flipping.5ht@gmail.com"
const files_dir = "Flipping/Datasets/FMRIOpto"
const figs_dir = "Flipping/Datasets/FMRIOpto/figures/"
## load
if ispath(linux_gdrive)
    ongoing_dir = linux_gdrive
else
    ongoing_dir = mac_gdrive
end

files_loc = joinpath(ongoing_dir,files_dir)
figs_loc = joinpath(ongoing_dir,figs_dir)

const columns_types = Dict(
    :Mouse => String,
    :Genotype => String,
    :Area => String,
    :Concentration => Float64,
    :Stim => Bool,
    :Repetition_nr => Int,
    :Trial_nr => Int,
    :Block => Int,
    :BOLD_90 => Float64,
    :BOLD_85 => Float64,
    :BOLD_0 => Float64,
    :Frame => Int,
    :Valid_frs => Bool,
    :Outlier_count_per_rep => Int
    )
using FMRIOpto
fulldata =  CSV.read(joinpath(files_loc,"odour_dose_data.csv"), DataFrame;
    types = columns_types,
    missingstring = "NaN")
mapcols!(Array, fulldata)
fulldata.Time = fulldata.Frame .+ 3
sort(fulldata.Area)
lvA = union(fulldata.Stim)
# fulldata.Area = categorical(fulldata.Area; levels = lvA)
# lvC = union(fulldata.Concentration)
# fulldata.Concentration = categorical(fulldata.Concentration; levels = lvC)
########
des = describe(fulldata)
open_html_table(des)
des2 = combine(groupby(fulldata, :Mouse),
    [:Genotype, :Area, :Concentration, :Stim, :Repetition_nr, :Trial_nr, :Block] .=> (c-> [union(c)]), nrow)
open_html_table(des2)
## average response per concentraition and stim over all presentation
fdata = filter(r ->
    r.Area != "PIR" &&
    # r.Area != "Ent" &&
    r.Valid_frs,
    # r.Mouse != "FM14",
    fulldata)
gd0 = groupby(fdata,[:Genotype, :Mouse, :Stim,:Area, :Concentration, :Time])
df0 = combine(gd0, [:BOLD_90, :BOLD_85,:BOLD_0] .=> mean .=> [:BOLD_90, :BOLD_85,:BOLD_0])
##
fdata2 = filter(r -> r.Genotype == "sert" &&
    r.Concentration != 0.0001,
    # r.Mouse != "FM14",
    df0)

gd1 = groupby(fdata2, [:Mouse, :Area, :Stim, :Concentration])

df1 = combine(gd1, [:BOLD_90, :BOLD_85,:BOLD_0] .=> max_bold .=> [:BOLD90_max6, :BOLD85_max6,:BOLD0_max6],
    [:BOLD_90, :Time] => ((b,t) -> average_boldtime(b,t)) => :BOLD90_4to10,
    [:BOLD_85, :Time] => ((b,t) -> average_boldtime(b,t)) => :BOLD85_4to10,
    [:BOLD_0, :Time] => ((b,t) -> average_boldtime(b,t)) => :BOLD0_4to10,
    )
## Areas that respond to concentration
df2 = filter(r -> !r.Stim, df1)
vaform2c = @formula(BOLD90_4to10 ~ 1 + Concentration + (1|Mouse))
mm2c = fit(MixedModel, vaform2c,df2)
vaform2 = @formula(BOLD90_4to10 ~ 1 + Concentration*Area + (1|Mouse))
mm2 = fit(MixedModel, vaform2,df2)
MixedModels.likelihoodratiotest(mm2c,mm2)

res = DataFrame(Area = coefnames(mm2), P = mm2.pvalues, Z = mm2.beta ./ mm2.stderror)
res2 = filter!(r -> occursin("Concentration & Area: ", r.Area), res)
res2.Area = replace.(res2.Area, "Concentration & Area: " => "")
sort!(res2,:P)
res3 = filter(r -> r.P <= 0.049, res2)
##
df3 = filter(r -> r.Area in res3.Area, df1)
vaform3c = @formula(BOLD90_4to10 ~ 1 + Concentration*Stim + Area*Stim + (1|Mouse))
mm3c = fit(MixedModel, vaform3c,df3)
vaform3 = @formula(BOLD90_4to10 ~ 1 + Concentration*Area*Stim + (1|Mouse))
mm3 = fit(MixedModel, vaform3,df3)
MixedModels.likelihoodratiotest(mm3c,mm3)

# show(MIME("text/latex"), mm1)
##
Area = "OB"
df4 = filter(r -> r.Area == Area, df1)
vaform4c = @formula(BOLD90_4to10 ~ 1 + Concentration + Stim + (1|Mouse))
mm4c = fit(MixedModel, vaform4c,df4)
vaform4 = @formula(BOLD90_4to10 ~ 1 + Concentration*Stim + (1|Mouse))
mm4 = fit(MixedModel, vaform4,df4)
MixedModels.likelihoodratiotest(mm4c,mm4)
##
Area = "OFC"
df5 = filter(r -> r.Area == Area, df1)
vaform5c = @formula(BOLD85_max6 ~ 1 + Concentration + Stim + (1|Mouse))
mm5c = fit(MixedModel, vaform5c,df5)
vaform5 = @formula(BOLD85_max6 ~ 1 + Concentration*Stim + (1|Mouse))
mm5 = fit(MixedModel, vaform5,df5)
MixedModels.likelihoodratiotest(mm5c,mm5)
##
fdata = filter(r ->
    r.Area != "PIR" &&
    r.Area != "OT" &&
    # r.Area != "Ent" &&
    r.Valid_frs,
    # r.Mouse != "FM14",
    fulldata)
    # lvA = union(fdata.Area)
    # fdata.Area = categorical(fdata.Area; levels = lvA)
gd0 = groupby(fdata,[:Genotype, :Mouse, :Stim, :Area, :Concentration])
df0 = combine(gd0) do dd
    any(ismissing,dd.BOLD_0) && println("got NaN ", dd.Mouse[1], dd.Area[1], dd.Stim[1], dd.Concentration[1])
    if all(ismissing.(dd.BOLD_0))
        return filter(r -> r.Area == "NoData",dd)
    else
        i = findmax(skipmissing(dd.BOLD_0))[2]
        return dd[i:end,:]
    end
end

df1 = filter(r -> r.Genotype == "sert" &&
    r.Concentration != 0.0001,
    # r.Mouse != "FM14",
    df0)
form_1c = @formula(BOLD_0 ~ 1 + Time + Concentration*Area + Stim*Area + (1+Time|Mouse+Area))
model_1c = fit(MixedModel, form_1c,df1)
form_1 = @formula(BOLD_0 ~ 1 + Time + Concentration*Area*Stim + (1+Time|Mouse+Area))
model_1 = fit(MixedModel, form_1,df1)
MixedModels.likelihoodratiotest(model_1c,model_1)
## Areas that respond to concentration
df2 = filter(r -> !r.Stim, df0)
form_2c = @formula(BOLD_90 ~ 1 + Time + Concentration + (1+Time|Mouse))
model_2c = fit(MixedModel, form_2c,df2)
form_2 = @formula(BOLD_90 ~ 1 + Time + Concentration*Area + (1+Time|Mouse))
model_2 = fit(MixedModel, form_2,df2)
MixedModels.likelihoodratiotest(model_2c,model_2)
res = DataFrame(Area = coefnames(model_2), P = model_2.pvalues, Z = model_2.beta ./ model_2.stderror)
res2 = filter!(r -> occursin("Concentration & Area: ", r.Area), res)
res2.Area = replace.(res2.Area, "Concentration & Area: " => "")
sort!(res2,:P)
res3 = filter(r -> r.P <= 0.049, res2)
##
df3 = filter(r -> r.Area in res3.Area, df1)
union(df3.Area)
df3.Area = categorical(df3.Area; levels = union(df3.Area))

form_3c = @formula(BOLD_90 ~ 1 + Time + Concentration*Stim + Area*Stim + (1+Time|Mouse))
model_3c = fit(MixedModel, form_3c,df3)

form_3 = @formula(BOLD_90 ~ 1 + Time + Concentration*Area*Stim + (1+Time|Mouse))
model_3 = fit(MixedModel, form_3,df3)

form_3b = @formula(BOLD_90 ~ 1 + Time + Concentration + Area + Stim + Concentration & Area + Concentration & Area & Stim + (1+Time|Mouse))
model_3b = fit(MixedModel, form_3b,df3)
MixedModels.likelihoodratiotest(model_3b,model_3)
##
Area = "OB"
df4 = filter(r -> r.Area == Area, df1)
form_4c = @formula(BOLD_90 ~ 1 + Time + Concentration + Stim + (1+Time|Mouse))
model_4c = fit(MixedModel, form_4c,df4)
form_4 = @formula(BOLD_90 ~ 1 + Time + Concentration*Stim + (1+Time|Mouse))
model_4 = fit(MixedModel, form_4,df4)
MixedModels.likelihoodratiotest(model_4c,model_4)
##
Area = "OFC"
df5 = filter(r -> r.Area == Area, df1)
form_5c = @formula(BOLD_90 ~ 1 + Time + Concentration + Stim + (1+Time|Mouse))
model_5c = fit(MixedModel, form_5c,df5)
form_5 = @formula(BOLD_90 ~ 1 + Time + Concentration*Stim + (1+Time|Mouse))
model_5 = fit(MixedModel, form_5,df5)
MixedModels.likelihoodratiotest(model_5c,model_5)
##
Area = "aPIR"
df5 = filter(r -> r.Area == Area, df1)
form_5c = @formula(BOLD_90 ~ 1 + Time + Concentration + Stim + (1+Time|Mouse))
model_5c = fit(MixedModel, form_5c,df5)
form_5 = @formula(BOLD_90 ~ 1 + Time + Concentration*Stim + (1+Time|Mouse))
model_5 = fit(MixedModel, form_5,df5)
MixedModels.likelihoodratiotest(model_5c,model_5)

model_5

##
vaform2 = @formula(BOLD_90 ~ 1 + Concentration*Area + (1|Mouse));
df2 = filter(r -> !r.Stim, df1)
mm2 = fit(MixedModel, vaform2,df2)

cnames = ["z(Intercept)",
    "zConcentration",
    "zStim",
    "zConcentration & Stim",
    "p(Intercept)",
    "pConcentration",
    "pStim",
    "pConcentration & Stim"]
allmodels = combine(groupby(fdata, :Area)) do dd
    try
        mm = fit(MixedModel, vaform2, dd)
        DataFrame(
            vcat("z".*coefnames(mm),"p".*coefnames(mm)) .=>
            vcat(mm.beta ./ mm.stderror,
            mm.pvalues)
            )
    catch
        # (mm = missing,)
        DataFrame(
            cnames .=>
            missing
            )
    end
end
open_html_table(allmodels)
##
coef(allmodels.mm[6])
fixef(allmodels.mm[6])
propertynames(allmodels.mm[6])

coefnames(allmodels.mm[6])
allmodels.mm[6].beta
allmodels.mm[6].stderror
allmodels.mm[6].beta ./ allmodels.mm[6].stderror
allmodels.mm[6].pvalues
DataFrame(
    vcat("z".*coefnames(allmodels.mm[6]),"p".*coefnames(allmodels.mm[6])) .=>
    vcat(allmodels.mm[6].beta ./ allmodels.mm[6].stderror,
    allmodels.mm[6].pvalues)
    )
DataFrame(
    vcat("z".*coefnames(allmodels.mm[6]),"p".*coefnames(allmodels.mm[6])) .=>
    missing
    )

vcat("z".*coefnames(allmodels.mm[6]),"p".*coefnames(allmodels.mm[6]))


text:
using Makie

"""
convert_to_cartesian(p)
Converts a point in the geographic coordinate system to a point in the
cartesian one with the given point in radians.
"""
function convert_to_cartesian(p)
    [cos(p[2]) * cos(p[1]),
     cos(p[2]) * sin(p[1]),
     sin(p[2])]
end

"""
convert_to_geographic(p)
Converts a point in the cartesian coordinate system to a point in the
geographic one with the given point.
"""
function convert_to_geographic(p)
    r = sqrt(p[1]^2 + p[2]^2 + p[3]^2)
    if p[1] > 0
          ϕ = atan(p[2] / p[1])
    elseif p[2] > 0
          ϕ = atan(p[2] / p[1]) + pi
    else
          ϕ = atan(p[2] / p[1]) - pi
    end
    θ = asin(p[3] / r)
    [ϕ, θ]
end

"""
λmap(p)
Sends a point on a 3-sphere to a point in the plane x₄=0 with the given point
in the form of 2 complex numbers representing a unit quaternion. This is the
stereographic projection of a 3-sphere. S³ ↦ R³
"""
function λmap(p)
    [real(p[1]), imag(p[1]), real(p[2])] ./ (1 - imag(p[2]))
end

"""
λ⁻¹map(p)
Sends a point on the plane back to a point on a unit sphere with the given
point. This is the inverse stereographic projection of a 3-sphere.
"""
function λ⁻¹map(p)
    x₁ = 2p[1] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₂ = 2p[2] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₃ = 2p[3] / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    x₄ = (-1 + p[1]^2 + p[2]^2 + p[3]^2) / (1 + p[1]^2 + p[2]^2 + p[3]^2)
    [Complex(x₁, x₂), Complex(x₃, x₄)]
end

"""
S¹action(α, p)
Performs a group action corresponding to moving along the circumference of a
circle with the given angle and the point in the form of 2 complex numbers
representing a unit quaternion.
"""
function S¹action(α, p)
    [exp(im * α) * p[1], exp(im * α) * p[2]]
end

"""
get_manifold(points, segments, start, finish)

Calculates a grid of points in R³ for constructing a surface in a specific way
with the given points in the base space, the number of segments, and the start
and finish angles that determine where to cut.
"""
function get_manifold(points, segments, start, finish)
    samples = size(points, 1)
    manifold = Array{Float64}(undef, segments, samples, 3)
    for i in 1:samples
            ϕ, θ = points[i, :]
        f = θ / (pi / 2)
        α = -(finish - start) / (segments - 1) / f
        p = S¹action(-start / f - pi/2, λ⁻¹map(convert_to_cartesian([ϕ, θ])))
        for j in 1:segments
            x₁, x₂, x₃ = λmap(p)
            manifold[j, i, :] = [x₁, x₂, x₃]
            p = S¹action(α, p)
        end
    end
    manifold
end

"""
build_surface(scene, points, color; transparency, shading)

Builds a surface with the given scene, points, color, transparency and shading.
"""
function build_surface(scene,
                       points,
                       color;
                       transparency = false,
                       shading = true)
    surface!(scene,
             @lift($points[:, :, 1]),
             @lift($points[:, :, 2]),
             @lift($points[:, :, 3]),
             color = color,
             transparency = transparency,
             shading = shading)
end

# The scene object that contains other visual objects
universe = Scene(backgroundcolor = :orange, show_axis=false, resolution = (360, 360))
# Use a slider for rotating the base space in an interactive way
sg, og = textslider(0:0.05:2pi, "g", start = 0)
segments = 100
samples = 100
# Construct a manifold for each country in the dictionary
latitudes = [pi/12, pi/6, pi/4, pi/3, 5pi/12]
for l in 1:length(latitudes)
    latitude = latitudes[l]
    # Sample a circle of points
    points = Array{Float64}(undef, samples, 2)
    for i in 1:samples
        longitude = i * 2pi / (samples - 1)
        points[i, :] = [longitude, latitude]
    end
    if (-1) ^ l > 0
        color = fill(RGBAf0(0.0, 0.0, 1.0, 1.0), segments, samples)
    else
        color = fill(RGBAf0(0.0, 0.0, 0.75, 1.0), segments, samples)
    end
    for i in range(1, stop = samples, step = 10)
        color[:, i] = fill(RGBAf0(rand(), rand(), rand(), 1.0), segments)
    end
    # Rotate the points by adding to the longitudes
    rotated = @lift begin
        R = similar(points)
        for i in 1:samples
                  ϕ, θ = points[i, :]
            R[i, :] = [ϕ + $og, θ]
        end
        R
    end
    start = -latitude
    finish = latitude
    manifold = @lift(get_manifold($rotated,
                                  segments,
                                  start,
                                  finish))
    build_surface(universe, manifold, color)
end

# Instantiate a horizontal box for holding the visuals and the controls
#scene = hbox(universe,
#             vbox(sg),
#             parent = Scene(resolution = (360, 360)))

# update eye position
eye_position, lookat, upvector = Vec3f0(-2, 0, 2), Vec3f0(0), Vec3f0(0, 0, 1.0)
update_cam!(universe, eye_position, lookat)
universe.center = false # prevent scene from recentering on display
rotate_cam!(universe, 0.0, 0.0, pi/2)

record(universe, "gallery/londontsai.gif") do io
    frames = 90
    for i in 1:frames
        og[] = i*2pi/frames # animate scene
        rotate_cam!(universe, 0.0, 2pi/frames, 0.0)
        recordframe!(io) # record a new frame
    end
end

max_issues_repo_name: iamazadi/Porta

max_forks_repo_name: iamazadi/Porta


text:
function file_header_V(file, params::Video_Params)
    #set up file header
    print(file, "order_run&")
    print(file, "video_number&")
    for i = 1:params.n_possible_objects
    	print(file, "fa_", string(i), "&")
    	print(file, "m_", string(i), "&")
    end

    for j = 1:num_particles
        for i = 1:params.n_possible_objects
        	print(file, "fa_", string(i), "_", string(j), "&")
        	print(file, "m_", string(i), "_", string(j), "&")
        end
        print(file, "weight_", string(j), "&")
    end
    print(file, "inferred_best_world_state")
    print(file, "\n")
end

function file_header_ws(file, params::Video_Params, num_particles::Int64)
    #set up file header
    print(file, "order_run&")
    print(file, "video_number&")
    print(file, "inferred_best_world_state&")
    for i = 1:num_particles-1
        print(file, "world_state_", string(i), "&")
        print(file, "weight_", string(i), "&")
    end
    print(file, "world_state_", num_particles, "&")
    print(file, "weight_", num_particles)
    print(file, "\n")
end

function file_header(file, params::Video_Params, num_particles::Int64)
    #set up file header
    print(file, "video_number&")
    for i = 1:params.n_possible_objects
    	print(file, "fa_", string(i), "&")
    	print(file, "m_", string(i), "&")
    end
    print(file, "inferred_world_states&inferred_best_world_state&")
    for i = 1:num_particles-1
        print(file, "world_state_", string(i), "&")
        print(file, "weight_", string(i), "&")
    end
    print(file, "world_state_", num_particles, "&")
    print(file, "weight_", num_particles)
    print(file, "\n")
end

function make_observations_full_COCO(dict::Array{Any,1}, receptive_fields::Vector{Receptive_Field}, num_videos::Int64, num_frames::Int64)
    objects_observed = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)
    #getting undefined reference when I change to Array{Array{}} instead of matrix

    camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)

    #temporary for getting min and max values of positions
    x_min = 0
    x_max = 0
    y_min = 3
    y_max = 0
    z_min = 0
    z_max = 0

    f_x_min = 0
    f_x_max = 0
    f_y_min = 3
    f_y_max = 0
    f_z_min = 0
    f_z_max = 0

    labels_max = 0

    for v=1:num_videos
        for f=1:num_frames
            #indices is where confidence was > 0.5
            indices = dict[v]["views"][f]["detections"]["scores"] .> 0.5
            labels = round.(Int64, dict[v]["views"][f]["detections"]["labels"][indices])
            center = dict[v]["views"][f]["detections"]["center"][indices]
            temp = Array{Detection2D}(undef, length(labels))
            if length(labels) > 5
                println("uh oh. too many obs ", length(labels))
                println("video ", v)
                println("frame ", f)
            end
            for i = 1:length(labels)
                label = labels[i]
                x = center[i][1]
                y = center[i][2]
                if x < 0 || x > 320
                    println("uh oh. x ", x)
                    println("video ", v)
                    println("frame ", f)
                elseif y < 0 || y > 240
                    println("uh oh. y ", y)
                    println("video ", v)
                    println("frame ", f)
                end
                temp[i] = (x, y, label) #so now I have an array of detections
            end
            #turn that array of detections into an array of an array of detections sorted by receptive_field
            #temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
            objects_observed[v, f] = temp

            #camera trajectory
            x = dict[v]["views"][f]["camera"][1]
            y = dict[v]["views"][f]["camera"][2]
            z = dict[v]["views"][f]["camera"][3]
            #focus
            f_x = dict[v]["views"][f]["lookat"][1]
            f_y = dict[v]["views"][f]["lookat"][2]
            f_z = dict[v]["views"][f]["lookat"][3]
            c = Camera_Params(camera_location = Coordinate(x,y,z), camera_focus = Coordinate(f_x,f_y,f_z))
            camera_trajectories[v, f] = c

            length(labels) > labels_max ? labels_max = length(labels) : labels_max = labels_max

            x > x_max ? x_max = x : x_max = x_max
            x < x_min ? x_min = x : x_min = x_min
            y > y_max ? y_max = y : y_max = y_max
            y < y_min ? y_min = y : y_min = y_min
            z > z_max ? z_max = z : z_max = z_max
            z < z_min ? z_min = z : z_min = z_min

            f_x > f_x_max ? f_x_max = f_x : f_x_max = f_x_max
            f_x < f_x_min ? f_x_min = f_x : f_x_min = f_x_min
            f_y > f_y_max ? f_y_max = f_y : f_y_max = f_y_max
            f_y < f_y_min ? f_y_min = f_y : f_y_min = f_y_min
            f_z > f_z_max ? f_z_max = f_z : f_z_max = f_z_max
            f_z < f_z_min ? f_z_min = f_z : f_z_min = f_z_min

        end
    end

    println("labels_max ", labels_max)

    println("x_max ", x_max)
    println("x_min ", x_min)
    println("y_max ", y_max)
    println("y_min ", y_min)
    println("z_max ", z_max)
    println("z_min ", z_min)

    println("f_x_max ", f_x_max)
    println("f_x_min ", f_x_min)
    println("f_y_max ", f_y_max)
    println("f_y_min ", f_y_min)
    println("f_z_max ", f_z_max)
    println("f_z_min ", f_z_min)

    return objects_observed, camera_trajectories
end

function make_observations_office(dict::Array{Any,1}, receptive_fields::Vector{Receptive_Field}, office_subset::Vector{String}, num_videos::Int64, num_frames::Int64, threshold=0.5, top_n=Inf)

    COCO_CLASSES = ["person", "bicycle", "car", "motorcycle",
    			"airplane", "bus", "train", "truck", "boat", "traffic light", "fire hydrant",
    			"N/A", "stop sign","parking meter", "bench", "bird", "cat", "dog", "horse",
    			"sheep", "cow", "elephant", "bear", "zebra", "giraffe", "N/A", "backpack",
    			"umbrella", "N/A", "N/A", "handbag", "tie", "suitcase", "frisbee", "skis", "snowboard",
    			"sports ball", "kite", "baseball bat", "baseball glove", "skateboard",
    			"surfboard", "tennis racket","bottle", "N/A", "wine glass", "cup", "fork", "knife",
    			"spoon","bowl", "banana", "apple", "sandwich", "orange", "broccoli", "carrot",
    			"hot dog", "pizza","donut", "cake", "chair", "couch", "potted plant", "bed",
    			"N/A", "dining table","N/A", "N/A", "toilet", "N/A", "tv", "laptop", "mouse",
    			"remote", "keyboard", "cell phone","microwave", "oven", "toaster", "sink", "refrigerator", "N/A",
    			"book","clock", "vase", "scissors", "teddy bear", "hair drier",
    			"toothbrush"]

    #office_subset = ["chair", "keyboard", "laptop", "dining table", "potted plant", "cell phone", "bottle"]
    #office_subset = ["book", "chair", "keyboard", "laptop", "table", "potted plant", "cell phone", "wine bottle"]
    #office_subset = ["dining table", "microwave", "backpack", "bowl", "couch"]

    objects_observed = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)
    #getting undefined reference when I change to Array{Array{}} instead of matrix

    camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)

    #temporary for getting min and max values of positions
    x_min = 0
    x_max = 0
    y_min = 0
    y_max = 0
    z_min = 0
    z_max = 0

    f_x_min = 0
    f_x_max = 0
    f_y_min = 3
    f_y_max = 0
    f_z_min = 0
    f_z_max = 0

    labels_max = 0

    for v=1:num_videos
        for f=1:num_frames
            #first sort by confidence level
            idx = sortperm(-1 .* dict[v]["views"][f]["detections"]["scores"]) #descending order
            dict[v]["views"][f]["detections"]["scores"] = dict[v]["views"][f]["detections"]["scores"][idx]
            dict[v]["views"][f]["detections"]["center"] = dict[v]["views"][f]["detections"]["center"][idx]
            dict[v]["views"][f]["detections"]["labels"] = dict[v]["views"][f]["detections"]["labels"][idx]

            indices = dict[v]["views"][f]["detections"]["scores"] .> threshold #indices should be sorted by confidence
            arr = round.(Int64, dict[v]["views"][f]["detections"]["labels"][indices]) #for turning to ints
            center = dict[v]["views"][f]["detections"]["center"][indices]
            #temp = Array{Detection2D}(undef, length(labels))
            # if length(labels) > 5
            #     println("uh oh. too many obs ", length(labels))
            #     println("video ", v)
            #     println("frame ", f)
            # end
            temp = []
            for i = 1:length(arr)
                label = findfirst(COCO_CLASSES[arr[i]] .== office_subset)
                if !isnothing(label)
                    x = center[i][1]
                    y = center[i][2]
                    push!(temp, (x, y, label))
                end
            end

            if length(temp) > 5
                println("v ", v, " f ", f)
                println(length(temp))
            end

            if length(temp) > top_n #temp should already be sorted by confidence. if it's not, will have to add code to sort it first
                #println(length(temp))
                #println("v ", v, " f ", f)
                println("v gets trimmed ", v)
                temp = temp[1:top_n]
            end

            #turn that array of detections into an array of an array of detections sorted by receptive_field
            #temp_sorted_into_rfs = map(rf -> filter(p -> within(p, rf), temp), receptive_fields)
            objects_observed[v, f] = convert(Array{Detection2D}, temp)

            #camera trajectory
            # x = dict[v]["views"][f]["camera"][1]
            # y = dict[v]["views"][f]["camera"][2]
            # z = dict[v]["views"][f]["camera"][3]
            x = dict[v]["views"][f]["camera"]["x"]
            y = dict[v]["views"][f]["camera"]["y"]
            z = dict[v]["views"][f]["camera"]["z"]
            #focus
            #f_x = dict[v]["views"][f]["lookat"][1]
            #f_y = dict[v]["views"][f]["lookat"][2]
            #f_z = dict[v]["views"][f]["lookat"][3]
            f_x = dict[v]["views"][f]["lookat"]["x"]
            f_y = dict[v]["views"][f]["lookat"]["y"]
            f_z = dict[v]["views"][f]["lookat"]["z"]
            c = Camera_Params(camera_location = Coordinate(x,y,z), camera_focus = Coordinate(f_x,f_y,f_z))
            camera_trajectories[v, f] = c

            length(temp) > labels_max ? labels_max = length(temp) : labels_max = labels_max

            x > x_max ? x_max = x : x_max = x_max
            x < x_min ? x_min = x : x_min = x_min
            y > y_max ? y_max = y : y_max = y_max
            y < y_min ? y_min = y : y_min = y_min
            z > z_max ? z_max = z : z_max = z_max
            z < z_min ? z_min = z : z_min = z_min

            f_x > f_x_max ? f_x_max = f_x : f_x_max = f_x_max
            f_x < f_x_min ? f_x_min = f_x : f_x_min = f_x_min
            f_y > f_y_max ? f_y_max = f_y : f_y_max = f_y_max
            f_y < f_y_min ? f_y_min = f_y : f_y_min = f_y_min
            f_z > f_z_max ? f_z_max = f_z : f_z_max = f_z_max
            f_z < f_z_min ? f_z_min = f_z : f_z_min = f_z_min

        end
    end

    println("labels_max ", labels_max)
    #
    println("x_max ", x_max)
    println("x_min ", x_min)
    println("y_max ", y_max)
    println("y_min ", y_min)
    println("z_max ", z_max)
    println("z_min ", z_min)

    println("f_x_max ", f_x_max)
    println("f_x_min ", f_x_min)
    println("f_y_max ", f_y_max)
    println("f_y_min ", f_y_min)
    println("f_z_max ", f_z_max)
    println("f_z_min ", f_z_min)

    return objects_observed, camera_trajectories
end

#for testing: feed ground-truth 2D detections
function make_observations_office_from_gt(dict::Array{Any,1}, receptive_fields::Vector{Receptive_Field}, office_subset::Vector{String},
    num_videos::Int64, num_frames::Int64, threshold=0.5)

    COCO_CLASSES = ["person", "bicycle", "car", "motorcycle",
    			"airplane", "bus", "train", "truck", "boat", "traffic light", "fire hydrant",
    			"N/A", "stop sign","parking meter", "bench", "bird", "cat", "dog", "horse",
    			"sheep", "cow", "elephant", "bear", "zebra", "giraffe", "N/A", "backpack",
    			"umbrella", "N/A", "N/A", "handbag", "tie", "suitcase", "frisbee", "skis", "snowboard",
    			"sports ball", "kite", "baseball bat", "baseball glove", "skateboard",
    			"surfboard", "tennis racket","bottle", "N/A", "wine glass", "cup", "fork", "knife",
    			"spoon","bowl", "banana", "apple", "sandwich", "orange", "broccoli", "carrot",
    			"hot dog", "pizza","donut", "cake", "chair", "couch", "potted plant", "bed",
    			"N/A", "dining table","N/A", "N/A", "toilet", "N/A", "tv", "laptop", "mouse",
    			"remote", "keyboard", "cell phone","microwave", "oven", "toaster", "sink", "refrigerator", "N/A",
    			"book","clock", "vase", "scissors", "teddy bear", "hair drier",
    			"toothbrush"]

    #office_subset = ["chair", "keyboard", "laptop", "dining table", "potted plant", "cell phone", "bottle"]
    #office_subset = ["book", "chair", "keyboard", "laptop", "table", "potted plant", "cell phone", "wine bottle"]
    #office_subset = ["dining table", "microwave", "backpack", "bowl", "couch"]

    objects_observed = Matrix{Array{Detection2D}}(undef, num_videos, num_frames)
    #getting undefined reference when I change to Array{Array{}} instead of matrix

    camera_trajectories = Matrix{Camera_Params}(undef, num_videos, num_frames)

    for v=1:num_videos
        for f=1:num_frames
            arr = dict[v]["views"][f]["ground_truth"]["labels"]
            center = dict[v]["views"][f]["ground_truth"]["centers"]

            temp = []
            for i = 1:length(arr)
                x = center[i][1]
                y = center[i][2]
                push!(temp, (x, y, arr[i]))
            end
            objects_observed[v, f] = convert(Array{Detection2D}, temp)

            #camera postion
            x = dict[v]["views"][f]["camera"]["x"]
            y = dict[v]["views"][f]["camera"]["y"]
            z = dict[v]["views"][f]["camera"]["z"]
            #focus
            f_x = dict[v]["views"][f]["lookat"]["x"]
            f_y = dict[v]["views"][f]["lookat"]["y"]
            f_z = dict[v]["views"][f]["lookat"]["z"]
            c = Camera_Params(camera_location = Coordinate(x,y,z), camera_focus = Coordinate(f_x,f_y,f_z))
            camera_trajectories[v, f] = c
        end
    end
    return objects_observed, camera_trajectories
end

#dict is actually an array of dictionaries
function write_to_dict(dict::Array{Any,1}, camera_trajectories::Matrix{Camera_Params}, inferred_realities, num_videos::Int64, num_frames::Int64)
    params = Video_Params()

    #add inferences about the objects in the scenes
    for v = 1:num_videos
        #for each object
        #parsed = eval(Meta.parse(inferred_realities[v])) #or eval(Meta.parse(inferred_realities[v][1]))
        a = Array{Any}(undef, length(inferred_realities[v]))
        for i in 1:length(inferred_realities[v])
            r = inferred_realities[v][i]
            a[i] = Dict("label" => r[4], "position" => r[1:3])
        end
        d = dict[v]
        merge!(d, Dict("inferences" => a))

        #add per frame
        for f = 1:num_frames
            labels = []
            centers = []
            for i in 1:length(inferred_realities[v])
                r = inferred_realities[v][i]
                x, y = get_image_xy(camera_trajectories[v,f], params, Coordinate(r[1],r[2],r[3]))
                if within((x, y), Receptive_Field((0,0), (params.image_dim_x, params.image_dim_y)))
                    push!(labels, r[4])
                    push!(centers, (x, y))
                end
            end
            d_f = d["views"][f]
            a = Dict("labels" => labels, "centers" => centers)
            merge!(d_f, Dict("inferences" => a))
        end
    end

    return dict
end

# function count_observations(objects_observed::Matrix{Array{Array{Detection2D}}})
#     num_videos = 100
#     num_frames = 75
#
#     zeros(91)
#
#     arr = []
#
#     for v=1:num_videos
#         arr_per_vid = []
#         for f=1:num_frames
#             n_obj = length(objects_observed[v, f])
#             for i = 1:n
#                 push!(arr_per_vid, objects_observed[v, f][])
#             end
#         end
#     end
#
# end

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen