text:
struct BoundUpdate{S<:Union{LT,GT}}
    cvi::CVI
    s::S
end

mutable struct Node
    lt_bounds::Vector{BoundUpdate{LT}}
    gt_bounds::Vector{BoundUpdate{GT}}
    lt_general_constrs::Vector{AffineConstraint{LT}}
    gt_general_constrs::Vector{AffineConstraint{GT}}
    depth::Int
    dual_bound::Float64
    bound_update::BoundUpdate
    fractional_value::Float64

    function Node()
        node = new()
        node.lt_bounds = BoundUpdate{LT}[]
        node.gt_bounds = BoundUpdate{GT}[]
        node.lt_general_constrs = AffineConstraint{LT}[]
        node.gt_general_constrs = AffineConstraint{LT}[]
        node.depth = 0
        node.dual_bound = -Inf
        # bound_update & fractional_value are left unitialized.
        return node
    end
end

function Node(
    lt_bounds::Vector{BoundUpdate{LT}},
    gt_bounds::Vector{BoundUpdate{GT}},
    depth::Int,
    dual_bound::Float64 = -Inf,
)
    node = Node()
    node.lt_bounds = lt_bounds
    node.gt_bounds = gt_bounds
    node.lt_general_constrs = []
    node.gt_general_constrs = []
    node.depth = depth
    node.dual_bound = dual_bound
    return node
end

function Node(
    lt_bounds::Vector{BoundUpdate{LT}},
    gt_bounds::Vector{BoundUpdate{GT}},
    lt_general_constrs::Vector{AffineConstraint{LT}},
    gt_general_constrs::Vector{AffineConstraint{GT}},
    depth::Int,
)
    node = Node()
    node.lt_bounds = lt_bounds
    node.gt_bounds = gt_bounds
    node.lt_general_constrs = lt_general_constrs
    node.gt_general_constrs = gt_general_constrs
    node.depth = depth
    node.dual_bound = -Inf
    return node
end

function Base.copy(node::Node)
    new_node =  Node(
        copy(node.lt_bounds),
        copy(node.gt_bounds),
        copy(node.lt_general_constrs),
        copy(node.gt_general_constrs),
        node.depth)
        new_node.dual_bound = node.dual_bound
        return new_node
end

mutable struct Tree
    open_nodes::DataStructures.Stack{Node}

    Tree() = new(DataStructures.Stack{Node}())
end

Base.isempty(tree::Tree) = Base.isempty(tree.open_nodes)
Base.length(tree::Tree) = Base.length(tree.open_nodes)
push_node!(tree::Tree, node::Node) = push!(tree.open_nodes, node)
pop_node!(tree::Tree) = pop!(tree.open_nodes)
num_open_nodes(tree::Tree) = length(tree.open_nodes)


text:
######### Stan program example  ###########

using StanSample

ProjDir = @__DIR__

diagnose_1 = "
parameters {
  real y;
  vector[9] x;
}
model {
  y ~ normal(0, 3);
  x ~ normal(0, exp(y/2));
}
";

diagnose_2 = "
parameters {
  real y_raw;
  vector[9] x_raw;
}
transformed parameters {
  real y;
  vector[9] x;

  y = 3.0 * y_raw;
  x = exp(y/2) * x_raw;
}
model {
  y_raw ~ std_normal(); // implies y ~ normal(0, 3)
  x_raw ~ std_normal(); // implies x ~ normal(0, exp(y/2))
}
";

sm1 = SampleModel("diagnose.1", diagnose_1);
rc1 = stan_sample(sm1);
if success(rc1)
  read_summary(sm1, true)
  diagnose(sm1)
end

sm2 = SampleModel("diagnose.2", diagnose_2);
rc2 = stan_sample(sm2);
if success(rc2)
  read_summary(sm2, true)
  diagnose(sm2)
end

ext: jl


text:
include("../src/SFFM.jl")
T =[-2.0 1.0 1; 1.0 -2.0 1; 4 3 -7]# [-2.0 2.0; 1.0 -1.0]#
C = [3.0; -2.0; 0]
fn(x) = [ones(size(x)) -ones(size(x)) -ones(size(x))]
model = SFFM.Model(;T=T,C=C,r=(r=fn,R=fn),Bounds=[0 10;-Inf Inf])

order = 2
mesh = SFFM.MakeMesh(model=model,Nodes = collect(0:5:10),NBases=order,Basis="lagrange")

# me = SFFM.MakeME(SFFM.CMEParams[order])
# erlang = SFFM.MakeErlang(order)

# me = SFFM.MakeME(SFFM.CMEParams[order], mean=mesh.Δ[1])
# B = SFFM.MakeBFRAP(model=model,mesh=mesh,me=me)
# R = SFFM.MakeR(model = model, mesh = mesh, approxType = "interpolation", probTransform = false)
# D = SFFM.MakeD(R=R,B=B,model=model,mesh=mesh)
# Ψme = SFFM.PsiFun(D=D)

mat = SFFM.MakeMatrices(model=model,mesh=mesh)
B2 = SFFM.MakeB(model=model,mesh=mesh,Matrices=mat)
All = SFFM.MakeAll(model = model, mesh = mesh, approxType = "interpolation")
Ψ = SFFM.PsiFun(D=All.D)

max_stars_repo_name: angus-lewis/SFFM

max_issues_repo_name: angus-lewis/SFFM

max_forks_repo_name: angus-lewis/SFFM


text:
"""
init_β(gcm)

Initialize the linear regression parameters `β` and `τ=σ0^{-2}` by the least
squares solution.
"""
function init_β!(
    gcm::Union{GaussianCopulaVCModel{T},GaussianCopulaLMMModel{T}}
    ) where T <: BlasReal
    # accumulate sufficient statistics X'y
    xty = zeros(T, gcm.p)
    for i in eachindex(gcm.data)
        BLAS.gemv!('T', one(T), gcm.data[i].X, gcm.data[i].y, one(T), xty)
    end
    # least square solution for β
    ldiv!(gcm.β, cholesky(Symmetric(gcm.XtX)), xty)
    # accumulate residual sum of squares
    rss = zero(T)
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rss += abs2(norm(gcm.data[i].res))
    end
    gcm.τ[1] = gcm.ntotal / rss
    gcm.β
end

"""
update_res!(gc, β)

Update the residual vector according to `β`.
"""
function update_res!(
    gc::Union{GaussianCopulaVCObs{T}, GaussianCopulaLMMObs{T}},
    β::Vector{T}
    ) where T <: BlasReal
    copyto!(gc.res, gc.y)
    BLAS.gemv!('N', -one(T), gc.X, β, one(T), gc.res)
    gc.res
end

function update_res!(
    gcm::Union{GaussianCopulaVCModel{T}, GaussianCopulaLMMModel{T}}
    ) where T <: BlasReal
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
    end
    nothing
end

function standardize_res!(
    gc::Union{GaussianCopulaVCObs{T}, GaussianCopulaLMMObs{T}},
    σinv::T
    ) where T <: BlasReal
    gc.res .*= σinv
end

function standardize_res!(
    gcm::Union{GaussianCopulaVCModel{T}, GaussianCopulaLMMModel{T}}
    ) where T <: BlasReal
    σinv = sqrt(gcm.τ[1])
    # standardize residual
    for i in eachindex(gcm.data)
        standardize_res!(gcm.data[i], σinv)
    end
    nothing
end

"""
update_quadform!(gc)
Update the quadratic forms `(r^T V[k] r) / 2` according to the current residual `r`.
"""
function update_quadform!(gc::GaussianCopulaVCObs)
    for k in 1:length(gc.V)
        gc.q[k] = dot(gc.res, mul!(gc.storage_n, gc.V[k], gc.res)) / 2
    end
    gc.q
end

"""

MM update to minimize ``n \\log (\\tau) - rss / 2 \\ln (\\tau) +
\\sum_i \\log (1 + \\tau * q_i)``.
"""
function update_τ(
    τ0::T,
    q::Vector{T},
    n::Integer,
    rss::T,
    maxiter::Integer=50000,
    reltol::Number=1e-6,
    ) where T <: BlasReal
    @assert τ0 ≥ 0 "τ0 has to be nonnegative"
    τ = τ0
    for τiter in 1:maxiter
        τold = τ
        tmp = zero(T)
        for i in eachindex(q)
            tmp += q[i] / (1 + τ * q[i])
        end
        τ = (n + 2τ * tmp) / rss
        abs(τ - τold) < reltol * (abs(τold) + 1) && break
    end
    τ
end

function update_Σ_jensen!(
    gcm::GaussianCopulaVCModel{T},
    maxiter::Integer=50000,
    reltol::Number=1e-6,
    verbose::Bool=false) where T <: BlasReal
    rsstotal = zero(T)
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rsstotal += abs2(norm(gcm.data[i].res))
        update_quadform!(gcm.data[i])
        gcm.QF[i, :] = gcm.data[i].q
    end
    # MM iteration
    for iter in 1:maxiter
        # store previous iterate
        copyto!(gcm.storage_Σ, gcm.Σ)
        # update τ
        mul!(gcm.storage_n, gcm.QF, gcm.Σ) # gcm.storage_n[i] = q[i]
        gcm.τ[1] = update_τ(gcm.τ[1], gcm.storage_n, gcm.ntotal, rsstotal, 1)
        # numerator in the multiplicative update
        gcm.storage_n .= inv.(inv(gcm.τ[1]) .+ gcm.storage_n) # use newest τ to update Σ
        mul!(gcm.storage_m, transpose(gcm.QF), gcm.storage_n)
        gcm.Σ .*= gcm.storage_m
        # denominator in the multiplicative update
        mul!(gcm.storage_n, gcm.TR, gcm.storage_Σ)
        gcm.storage_n .= inv.(1 .+ gcm.storage_n)
        mul!(gcm.storage_m, transpose(gcm.TR), gcm.storage_n)
        gcm.Σ ./= gcm.storage_m
        # monotonicity diagnosis
        verbose && println(sum(log, 1 .+ gcm.τ[1] .* (gcm.QF * gcm.Σ)) -
            sum(log, 1 .+ gcm.TR * gcm.Σ) +
            gcm.ntotal / 2 * (log(gcm.τ[1]) - log(2π)) -
            rsstotal / 2 * gcm.τ[1])
        # convergence check
        gcm.storage_m .= gcm.Σ .- gcm.storage_Σ
        # norm(gcm.storage_m) < reltol * (norm(gcm.storage_Σ) + 1) && break
        if norm(gcm.storage_m) < reltol * (norm(gcm.storage_Σ) + 1)
            verbose && println("iters=$iter")
            break
        end
        verbose && iter == maxiter && @warn "maximum iterations $maxiter reached"
    end
    gcm.Σ
end

function update_Σ_quadratic!(
    gcm::GaussianCopulaVCModel{T},
    maxiter::Integer=50000,
    reltol::Number=1e-6,
    qpsolver=Ipopt.IpoptSolver(print_level=0),
    verbose::Bool=false) where T <: BlasReal
    n, m = length(gcm.data), length(gcm.data[1].V)
    # pre-compute quadratic forms and RSS
    rsstotal = zero(T)
    for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rsstotal += abs2(norm(gcm.data[i].res))
        update_quadform!(gcm.data[i])
        gcm.QF[i, :] = gcm.data[i].q
    end
    qcolsum = sum(gcm.QF, dims=1)[:]
    # define NNLS optimization problem
    H = Matrix{T}(undef, m, m)  # quadratic coefficient in QP
    c = Vector{T}(undef, m)     # linear coefficient in QP
    w = Vector{T}(undef, n)
    # MM iteration
    for iter in 1:maxiter
        # store previous iterate
        copyto!(gcm.storage_Σ, gcm.Σ)
        # update τ
        mul!(gcm.storage_n, gcm.QF, gcm.Σ) # gcm.storage_n[i] = q[i]
        # a, b = zero(T), - rsstotal / 2
        # for i in eachindex(gcm.data)
        #     a += abs2(gcm.storage_n[i]) / (1 + gcm.τ[1] * gcm.storage_n[i])
        #     b += gcm.storage_n[i]
        # end
        # gcm.τ[1] = (b + sqrt(abs2(b) + 2a * gcm.ntotal)) / 2a
        tmp = zero(T)
        for i in eachindex(gcm.data)
            tmp += gcm.storage_n[i] / (1 + gcm.τ[1] * gcm.storage_n[i])
        end
        gcm.τ[1] = (gcm.ntotal + 2gcm.τ[1] * tmp) / rsstotal  # update τ
        # update variance components
        for i in eachindex(gcm.data)
            w[i] = abs2(gcm.τ[1]) / (1 + gcm.τ[1] * gcm.storage_n[i])
        end
        mul!(H, transpose(gcm.QF) * Diagonal(w), gcm.QF)
        mul!(gcm.storage_n, gcm.TR, gcm.storage_Σ)
        gcm.storage_n .= inv.(1 .+ gcm.storage_n)
        mul!(gcm.storage_m, transpose(gcm.TR), gcm.storage_n)
        c .= gcm.τ[1] .* qcolsum .- gcm.storage_m
        # try unconstrained solution first
        ldiv!(gcm.Σ, cholesky(Symmetric(H)), c)
        # if violate nonnegativity constraint, resort to quadratic programming
        if any(x -> x < 0, gcm.Σ)
            @show "use QP"
            qpsol = quadprog(-c, H, Matrix{T}(undef, 0, m),
                Vector{Char}(undef, 0), Vector{T}(undef, 0),
                fill(T(0), m), fill(T(Inf), m), qpsolver)
            gcm.Σ .= qpsol.sol
        end
        # monotonicity diagnosis
        verbose && println(sum(log, 1 .+ gcm.τ[1] .* (gcm.QF * gcm.Σ)) -
            sum(log, 1 .+ gcm.TR * gcm.Σ) +
            gcm.ntotal / 2 * (log(gcm.τ[1]) - log(2π)) -
            rsstotal / 2 * gcm.τ[1])
        # convergence check
        gcm.storage_m .= gcm.Σ .- gcm.storage_Σ
        if norm(gcm.storage_m) < reltol * (norm(gcm.storage_Σ) + 1)
            println("iters=$iter")
            break
        end
        verbose && iter == maxiter && @warn "maximum iterations $maxiter reached"
    end
    gcm.Σ
end

"""
update_Σ!(gc)

Update `τ` and variance components `Σ` according to the current value of
`β` by an MM algorithm. `gcm.QF` needs to hold qudratic forms calculated from
un-standardized residuals.
"""
update_Σ! = update_Σ_jensen!

function fitted(
    gc::GaussianCopulaVCObs{T},
    β::Vector{T},
    τ::T,
    Σ::Vector{T}) where T <: BlasReal
    n, m = length(gc.y), length(gc.V)
    μ̂ = gc.X * β
    Ω = Matrix{T}(undef, n, n)
    for k in 1:m
        Ω .+= Σ[k] .* gc.V[k]
    end
    σ02 = inv(τ)
    c = inv(1 + dot(Σ, gc.t)) # normalizing constant
    V̂ = Matrix{T}(undef, n, n)
    for j in 1:n
        for i in 1:j-1
            V̂[i, j] = c * σ02 * Ω[i, j]
        end
        V̂[j, j] = c * σ02 * (1 + Ω[j, j] + tr(Ω) / 2)
    end
    LinearAlgebra.copytri!(V̂, 'U')
    μ̂, V̂
end

function loglikelihood!(
    gc::GaussianCopulaVCObs{T},
    β::Vector{T},
    τ::T, # inverse of linear regression variance
    Σ::Vector{T},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: BlasReal
    n, p, m = size(gc.X, 1), size(gc.X, 2), length(gc.V)
    needgrad = needgrad || needhess
    if needgrad
        fill!(gc.∇β, 0)
        fill!(gc.∇τ, 0)
        fill!(gc.∇Σ, 0)
    end
    needhess && fill!(gc.Hβ, 0)
    # evaluate copula loglikelihood
    sqrtτ = sqrt(τ)
    update_res!(gc, β)
    standardize_res!(gc, sqrtτ)
    rss  = abs2(norm(gc.res)) # RSS of standardized residual
    tsum = dot(Σ, gc.t)
    logl = - log(1 + tsum) - (n * log(2π) -  n * log(τ) + rss) / 2
    for k in 1:m
        mul!(gc.storage_n, gc.V[k], gc.res) # storage_n = V[k] * res
        if needgrad # ∇β stores X'*Γ*res (standardized residual)
            BLAS.gemv!('T', Σ[k], gc.X, gc.storage_n, one(T), gc.∇β)
        end
        gc.q[k] = dot(gc.res, gc.storage_n) / 2
    end
    qsum  = dot(Σ, gc.q)
    logl += log(1 + qsum)
    # gradient
    if needgrad
        inv1pq = inv(1 + qsum)
        if needhess
            BLAS.syrk!('L', 'N', -abs2(inv1pq), gc.∇β, one(T), gc.Hβ) # only lower triangular
            gc.Hτ[1, 1] = - abs2(qsum * inv1pq / τ)
        end
        BLAS.gemv!('T', one(T), gc.X, gc.res, -inv1pq, gc.∇β)
        gc.∇β .*= sqrtτ
        gc.∇τ  .= (n - rss + 2qsum * inv1pq) / 2τ
        gc.∇Σ  .= inv1pq .* gc.q .- inv(1 + tsum) .* gc.t
    end
    # output
    logl
end

function loglikelihood!(
    gcm::Union{GaussianCopulaVCModel{T}, GaussianCopulaLMMModel{T}},
    needgrad::Bool = false,
    needhess::Bool = false
    ) where T <: BlasReal
    logl = zero(T)
    if needgrad
        fill!(gcm.∇β, 0)
        fill!(gcm.∇τ, 0)
        fill!(gcm.∇Σ, 0)
    end
    if needhess
        gcm.Hβ .= - gcm.XtX
        gcm.Hτ .= - gcm.ntotal / 2abs2(gcm.τ[1])
    end
    for i in eachindex(gcm.data)
        logl += loglikelihood!(gcm.data[i], gcm.β, gcm.τ[1], gcm.Σ, needgrad, needhess)
        if needgrad
            gcm.∇β .+= gcm.data[i].∇β
            gcm.∇τ .+= gcm.data[i].∇τ
            gcm.∇Σ .+= gcm.data[i].∇Σ
        end
        if needhess
            gcm.Hβ .+= gcm.data[i].Hβ
            gcm.Hτ .+= gcm.data[i].Hτ
        end
    end
    needhess && (gcm.Hβ .*= gcm.τ[1])
    logl
end

function fit!(
    gcm::GaussianCopulaVCModel,
    solver=Ipopt.IpoptSolver(print_level=0)
    )
    optm = MathProgBase.NonlinearModel(solver)
    lb = fill(-Inf, gcm.p)
    ub = fill( Inf, gcm.p)
    MathProgBase.loadproblem!(optm, gcm.p, 0, lb, ub, Float64[], Float64[], :Max, gcm)
    MathProgBase.setwarmstart!(optm, gcm.β)
    MathProgBase.optimize!(optm)
    optstat = MathProgBase.status(optm)
    optstat == :Optimal || @warn("Optimization unsuccesful; got $optstat")
    copy_par!(gcm, MathProgBase.getsolution(optm))
    loglikelihood!(gcm)
    gcm
end

function MathProgBase.initialize(
    gcm::GaussianCopulaVCModel,
    requested_features::Vector{Symbol})
    for feat in requested_features
        if !(feat in [:Grad])
            error("Unsupported feature $feat")
        end
    end
end

MathProgBase.features_available(gcm::GaussianCopulaVCModel) = [:Grad]

function MathProgBase.eval_f(
    gcm::GaussianCopulaVCModel,
    par::Vector)
    copy_par!(gcm, par)
    # maximize σ2 and τ at current β using MM
    update_Σ!(gcm)
    # evaluate loglikelihood
    loglikelihood!(gcm, false, false)
end

function MathProgBase.eval_grad_f(
    gcm::GaussianCopulaVCModel,
    grad::Vector,
    par::Vector)
    copy_par!(gcm, par)
    # maximize σ2 and τ at current β using MM
    update_Σ!(gcm)
    # evaluate gradient
    logl = loglikelihood!(gcm, true, false)
    copyto!(grad, gcm.∇β)
    nothing
end

function copy_par!(
    gcm::GaussianCopulaVCModel,
    par::Vector)
    copyto!(gcm.β, par)
    par
end

function MathProgBase.hesslag_structure(gcm::GaussianCopulaVCModel)
    Iidx = Vector{Int}(undef, (gcm.p * (gcm.p + 1)) >> 1)
    Jidx = similar(Iidx)
    ct = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            Iidx[ct] = i
            Jidx[ct] = j
            ct += 1
        end
    end
    Iidx, Jidx
end

function MathProgBase.eval_hesslag(
    gcm::GaussianCopulaVCModel{T},
    H::Vector{T},
    par::Vector{T},
    σ::T,
    μ::Vector{T}) where T <: BlasReal
    copy_par!(gcm, par)
    # maximize σ2 and τ at current β using MM
    update_Σ!(gcm)
    # evaluate Hessian
    loglikelihood!(gcm, true, true)
    # copy Hessian elements into H
    ct = 1
    for j in 1:gcm.p
        for i in j:gcm.p
            H[ct] = gcm.Hβ[i, j]
            ct += 1
        end
    end
    H .*= σ
end


text:
"""
    NodePosition

Represents the position and velocity of a node during force directed layouts. The
fields are `x` and `y` for position, and `vx` and `vy` for their relative
velocity.
"""
mutable struct NodePosition
    x::Float64
    y::Float64
    vx::Float64
    vy::Float64
    r::Number
end

function NodePosition()
    n1, n2 = (rand(2).*2.0).-1.0
    while ((n1^2.0)+(n2^2.0))≥1.0
        n1, n2 = (rand(2).*2.0).-1.0
    end
    NodePosition(n1, n2, 0.0, 0.0, 0.0)
end
NodePosition(x::Float64, y::Float64) = NodePosition(x, y, 0.0, 0.0, 0.0)
NodePosition(x::Float64, y::Float64, vx::Float64, vy::Float64) = NodePosition(x, y, vx, vy, 0.0)

function NodePosition(r::T) where {T <: Number}
    return NodePosition(0.0, 0.0, 0.0, 0.0, r)
end

function NodePosition(x::Float64, y::Float64, r::T) where {T <: Number}
    return NodePosition(x, y, 0.0, 0.0, r)
end

"""
    BipartiteInitialLayout

This type is used to generate an initial bipartite layout, where the nodes are
placed on two levels, but their horizontal position is random.
"""
struct BipartiteInitialLayout end

"""
    FoodwebInitialLayout

This type is used to generate an initial layout, where the nodes are
placed on their trophic levels, but their horizontal position is random.
"""
struct FoodwebInitialLayout end

"""
    RandomInitialLayout

This type is used to generate an initial layout, where the nodes are
placed at random within the unit circle.
"""
struct RandomInitialLayout end

"""
    CircularInitialLayout

This type is used to generate an initial layout, where the nodes are
placed at random along a circle.
"""
struct CircularInitialLayout end

"""
    UnravelledInitialLayout

This type is used to generate an initial unravelled layout, where the nodes are
sorted vertically by trophic level, and horizontally by omnivory index. Credit
for this approach goes to @gvdr -- https://github.com/gvdr/unravel#unravel
"""
struct UnravelledInitialLayout end


text:
module TestLoadAllPackages

using LoadAllPackages
using Test

@testset """@loadall("../Project.toml")""" begin
    @test_logs(
        (:info, r"Loading packages from .*"),
        (:info, r"packages loaded in .* seconds"),
        match_mode = :any,
        LoadAllPackages.@loadall("../Project.toml")
    )
end

@testset """@loadall("..")""" begin
    @test_logs(
        (:info, r"Loading packages from .*"),
        (:info, r"packages loaded in .* seconds"),
        match_mode = :any,
        LoadAllPackages.@loadall("..")
    )
end

@testset "error handling" begin
    mktempdir() do emptydir
        @test_throws(
            ErrorException("`$emptydir` does not have `JuliaProject.toml` or `Project.toml`"),
            LoadAllPackages.@loadall(emptydir)
        )

        nonexisting = joinpath(emptydir, "nonexisting")
        @test_throws(
            ErrorException("`$nonexisting` is not a file or directory"),
            LoadAllPackages.@loadall(nonexisting)
        )
    end
end

end  # module


text:
versions = 56:-1:36
apt_versions = [52, 48, 42]


text:
using Phylodendron2, Test

@testset "All nodes unlabelled , no branch lengths." begin
    nwk = "(,,(,));"
    tree_origin = parse_newick(nwk)

    @testset "Tree structure" begin
        @test n_neighbour(tree_origin) == 3

        global a, b, c = neighbours(tree_origin)

        @test neighbours(a) == [tree_origin]
        @test neighbours(b) == [tree_origin]
        @test n_neighbour(c) == 3
        @test neighbours(c)[1] == tree_origin

        global d, e = neighbours(c)[2:3]
        @test neighbours(d) == [c]
        @test neighbours(e) == [c]
    end
    @testset "Node and branch properties" begin
        @test label(tree_origin) == ""
        @test label(c) == ""
        @test isnothing(brlength(tree_origin, a))
        @test isnothing(brlength(c, d))
        @test brlabel(tree_origin, a) == ""
        @test brlabel(c, d) == ""
    end
    @testset "Writing Newick" begin
        @test newick_string(tree_origin) == nwk
    end
end

# Clear variables
tree_origin, a, b, c, d, e = nothing, nothing, nothing, nothing, nothing, nothing

@testset "Labelled tips, no branch lengths" begin
    nwk = "(A,B,(C,D));"
    tree_origin = parse_newick(nwk)
    @testset "Tree structure" begin
        @test n_neighbour(tree_origin) == 3

        global a, b, e = neighbours(tree_origin)

        @test neighbours(a) == [tree_origin]
        @test neighbours(b) == [tree_origin]
        @test n_neighbour(e) == 3
        @test neighbours(e)[1] == tree_origin
        global c, d = neighbours(e)[2:3]
        @test neighbours(c) == [e]
        @test neighbours(d) == [e]
    end
    @testset "Node and branch properties" begin
        @test label(tree_origin) == ""
        @test label(a) == "A"
        @test label(b) == "B"
        @test label(c) == "C"
        @test label(d) == "D"
        @test label(e) == ""
        @test isnothing(brlength(tree_origin, a))
        @test isnothing(brlength(c, e))
        @test brlabel(tree_origin, a) == ""
        @test brlabel(c, e) == ""
    end
    @testset "Writing Newick" begin
        @test newick_string(tree_origin) == nwk
    end
end

# Clear variables
tree_origin, a, b, c, d, e = nothing, nothing, nothing, nothing, nothing, nothing

@testset "All nodes labelled, no branch lengths" begin
    nwk = "(A,B,(C,D)E)F;"
    f = parse_newick(nwk)
    @testset "Tree structure" begin
        @test n_neighbour(f) == 3

        global a, b, e = neighbours(f)

        @test neighbours(a) == [f]
        @test neighbours(b) == [f]
        @test n_neighbour(e) == 3
        @test neighbours(e)[1] == f
        global c, d = neighbours(e)[2:3]
        @test neighbours(c) == [e]
        @test neighbours(d) == [e]
    end
    @testset "Node and branch properties" begin
        @test label(a) == "A"
        @test label(b) == "B"
        @test label(c) == "C"
        @test label(d) == "D"
        @test label(e) == "E"
        @test label(f) == "F"
        @test isnothing(brlength(f, a))
        @test isnothing(brlength(e, d))
        @test brlabel(f, a) == ""
        @test brlabel(e, d) == ""
    end
    @testset "Writing Newick" begin
        @test newick_string(f) == nwk
    end
end

# Clear variables
tree_origin, a, b, c, d, e = nothing, nothing, nothing, nothing, nothing, nothing

@testset "All nodes unlabelled, with branch lengths." begin
    nwk = "(:0.1,:0.2,(:0.3,:0.4):0.5);"
    tree_origin = parse_newick(nwk)

    @testset "Tree structure" begin
        @test n_neighbour(tree_origin) == 3

        global a, b, c = neighbours(tree_origin)

        @test neighbours(a) == [tree_origin]
        @test neighbours(b) == [tree_origin]
        @test n_neighbour(c) == 3
        @test neighbours(c)[1] == tree_origin

        global d, e = neighbours(c)[2:3]
        @test neighbours(d) == [c]
        @test neighbours(e) == [c]
    end
    @testset "Node and branch properties" begin
        @test label(tree_origin) == ""
        @test label(c) == ""
        @test brlength(tree_origin, a) == 0.1
        @test brlength(tree_origin, b) == 0.2
        @test brlength(tree_origin, c) == 0.5
        @test brlength(c, d) == 0.3
        @test brlength(c, e) == 0.4
        @test brlabel(tree_origin, a) == ""
        @test brlabel(c, d) == ""
    end
    @testset "Writing Newick" begin
        @test newick_string(tree_origin) == nwk
    end
end

# Clear variables
tree_origin, a, b, c, d, e = nothing, nothing, nothing, nothing, nothing, nothing

@testset "Labelled tips, with branch lengths" begin
    nwk = "(A:0.1,B:0.2,(C:0.3,D:0.4):0.5);"
    tree_origin = parse_newick(nwk)
    @testset "Tree structure" begin
        @test n_neighbour(tree_origin) == 3

        global a, b, e = neighbours(tree_origin)

        @test neighbours(a) == [tree_origin]
        @test neighbours(b) == [tree_origin]
        @test n_neighbour(e) == 3
        @test neighbours(e)[1] == tree_origin
        global c, d = neighbours(e)[2:3]
        @test neighbours(c) == [e]
        @test neighbours(d) == [e]
    end
    @testset "Node and branch properties" begin
        @test label(tree_origin) == ""
        @test label(a) == "A"
        @test label(b) == "B"
        @test label(c) == "C"
        @test label(d) == "D"
        @test label(e) == ""
        @test brlength(tree_origin, a) == 0.1
        @test brlength(tree_origin, b) == 0.2
        @test brlength(tree_origin, e) == 0.5
        @test brlength(e, c) == 0.3
        @test brlength(e, d) == 0.4
        @test brlabel(tree_origin, a) == ""
        @test brlabel(c, e) == ""
    end
    @testset "Writing Newick" begin
        @test newick_string(tree_origin) == nwk
    end
end

# Clear variables
tree_origin, a, b, c, d, e = nothing, nothing, nothing, nothing, nothing, nothing

@testset "All nodes labelled, with branch lengths" begin
    nwk = "(A:0.1,B:0.2,(C:0.3,D:0.4)E:0.5)F;"
    f = parse_newick(nwk)
    @testset "Tree structure" begin
        @test n_neighbour(f) == 3

        global a, b, e = neighbours(f)

        @test neighbours(a) == [f]
        @test neighbours(b) == [f]
        @test n_neighbour(e) == 3
        @test neighbours(e)[1] == f
        global c, d = neighbours(e)[2:3]
        @test neighbours(c) == [e]
        @test neighbours(d) == [e]
    end
    @testset "Node and branch properties" begin
        @test label(a) == "A"
        @test label(b) == "B"
        @test label(c) == "C"
        @test label(d) == "D"
        @test label(e) == "E"
        @test label(f) == "F"
        @test brlength(f, a) == 0.1
        @test brlength(f, b) == 0.2
        @test brlength(f, e) == 0.5
        @test brlength(e, c) == 0.3
        @test brlength(e, d) == 0.4
        @test brlabel(f, a) == ""
        @test brlabel(e, d) == ""
    end
    @testset "Making Newick strings" begin
        @test newick_string(f) == nwk
    end
end

@testset "Reading from file" begin
    @testset "Default settings" begin
        trees = read_newick("../data/snouters.nwk")
        @test typeof(trees) <: TreeVector
        @test length(trees) == 1
        @test n_tip(trees[1]) == 10
        trees = read_newick("../data/wikipedia.nwk")
        @test typeof(trees) <: TreeVector
        @test length(trees) == 14
        @test n_tip(trees[1]) == 4
    end
    @testset "n hint too low" begin
        trees = read_newick("../data/wikipedia.nwk"; nhint=1)
        @test typeof(trees) <: TreeVector
        @test length(trees) == 14
        @test n_tip(trees[1]) == 4
    end
end

@testset "Reading trees from string" begin
    text = [randtree(12) for _ in 1:20] .|> newick_string |> join
    trees = read_newick(text=text)

    @test length(trees) == 20
    @test text == newick_string.(trees) |> join
end

@testset "Newick trees with emoji" begin
    text_tree = "((((🐇,🐳),(🐍:12.0,🐢)),🐠),🐌);"
    tree = parse_newick(text_tree)

    @test text_tree == newick_string(tree)
end

# tree = read_newick("../../data/snouters.nwk")
# trees = read_newick("trees/wikipedia.nwk")

# @testset "Writing to file" begin
#     write_newick("trees/test1.tmp", tree[1])
#     test1 = read_newick("trees/test1.tmp")
#     @test typeof(test1) == TreeVector
#     @test length(test1) == 1

#     write_newick("trees/test2.tmp", tree)
#     test2 = read_newick("trees/test2.tmp")
#     @test typeof(test2) == TreeVector
#     @test length(test2) == 1

#     write_newick("trees/test3.tmp", trees)
#     test3 = read_newick("trees/test3.tmp")
#     @test typeof(test3) == TreeVector
#     @test length(test3) == 14
# end

# rm("trees/test1.tmp")
# rm("trees/test2.tmp")
# rm("trees/test3.tmp")


text:
@testset "acrobot_env" begin

    env = AcrobotEnv(; rng = MersenneTwister(123))
    RLBase.test_interfaces!(env)
    RLBase.test_runnable!(env)

end


text:
module FileIO

using DelimitedFiles, MatrixMarket, MAT, CSV, DataFrames, HDF5
export readtxt,   # Text file
	   readmtx,   # MatrixMarket file
	   readmat,   # MATLAB mat file
	   readhdf,   # HDF5 file
       readgenelist

"""
    readmtx(filename::AbstractString) -> AbstractMatrix

Load and return the content of the MatrixMarket file.

# Arguments

- `filename::AbstractString`: the filename of the file to load,

# Keywords

# Returns
- `AbstractMatrix`: the UMI matrix of the file content,

# Throws
- `Error`: in the case of file does not exist.
"""

function readmtx(filename::AbstractString)
    if isfile(filename)
	# read MatrixMarket file
	X=MatrixMarket.mmread(filename);
    else
        error("File $(filename) does not exist.")
    end
    return X
end

function readcsv(filename::AbstractString)
    # read CSV file with header and row name
    df = CSV.File(filename; datarow=2) |> DataFrame!
    X=convert(Matrix, df[:,2:end])
    genelist=df[:,1]
    return X,genelist
end

function readmat(filename::AbstractString)
    # read Matlab Mat file
    file=matopen(filename)
    X=read(file,"X")
    genelist=read(file,"genelist")
    close(file)
    X = convert(Array{Float64,2}, X)
    return X,genelist
end

function readtxt(filename::AbstractString)
    # read DLM text file
    X=readdlm(filename,',',Int16)
end

function readhdf(filename::AbstractString)
    # https://anndata.readthedocs.io/en/latest/anndata.AnnData.html
end

function readgenelist(filename::AbstractString,colidx::Integer=1)
    # read genelist
    genelist=readdlm(filename,'\t',String)
    genelist=vec(genelist[:,colidx])
end

# using UnicodePlots
# function showx(X)
#    spy(X)
# end

end


text:
using Oceananigans.Utils: prettysummary

struct ZeroField{T, N} <: AbstractField{Nothing, Nothing, Nothing, Nothing, T, N} end
struct OneField{T, N} <: AbstractField{Nothing, Nothing, Nothing, Nothing, T, N} end

ZeroField(T=Int) = ZeroField{T, 3}() # default 3D, integer 0
OneField(T=Int) = OneField{T, 3}() # default 3D, integer 1

@inline Base.getindex(::ZeroField{T, N}, ind...) where {N, T} = zero(T)
@inline Base.getindex(::OneField{T, N}, ind...) where {N, T} = one(T)

struct ConstantField{T, N} <: AbstractField{Nothing, Nothing, Nothing, Nothing, T, N}
    constant :: T
    ConstantField{N}(constant::T) where {T, N} = new{T, N}(constant)
end

# Default 3-dimensional field
ConstantField(constant) = ConstantField{3}(constant)

@inline Base.getindex(f::ConstantField, ind...) = f.constant

const CF = Union{ConstantField, ZeroField, OneField}

@inline Base.summary(f::CF) = string("ConstantField(", prettysummary(f.constant), ")")
Base.show(io::IO, f::CF) = print(io, summary(f))


text:
function reduce_array_for_convex_hull!(point_array::Vector{<:AbstractPoint})::Vector{<:AbstractPoint}
    if length(point_array) <= 8
        return point_array
    end

    (min_x, max_x, min_y, max_y) = extrema(point_array)
    (skew_min_x, skew_max_x, skew_min_y, skew_max_y) = skew_extrema(point_array)

    approx_hull_array = [
        min_x,
        skew_min_x,
        min_y,
        skew_min_y,
        max_x,
        skew_max_x,
        max_y,
        skew_max_y]

    # approx_hull_array = [
    #     min_x,
    #     min_y,
    #     max_x,
    #     max_y]

    u_approx_hull_array = unique(p->p.x, approx_hull_array)
    uu_approx_hull_array = unique(p->p.y, u_approx_hull_array)
    approx_hull = naive_graham_scan!(uu_approx_hull_array)

    N_triangles = length(approx_hull)-2
    triangles = Vector{Triangle{Point{Float64}}}(undef, N_triangles)
    for kk = 1:N_triangles
        triangles[kk] = Triangle(approx_hull[1], approx_hull[kk+1], approx_hull[kk+2])
    end

    # bool = falses( length(point_array) )
    # for kk = 1:N_triangles
    #     bool_tmp = is_in_triangle.([triangles[kk]], point_array)
    #     bool = bool .| bool_tmp
    # end

    # reduced_array = [
    #     approx_hull_array...,
    #     point_array[.!bool]...
    #     ]

    for kk = 1:N_triangles
        bool_tmp = is_in_triangle.([triangles[kk]], point_array)
        point_array = point_array[.!bool_tmp]
    end

    uu_reduced_array = [
        approx_hull...,
        point_array...
        ]

    # u_reduced_array = unique(p->p.x, reduced_array)
    # uu_reduced_array = unique(p->p.y, u_reduced_array)
    return uu_reduced_array
end

function graham_scan!(point_array::Vector{<:AbstractPoint})
    reduced_array = reduce_array_for_convex_hull!(point_array)
    return naive_graham_scan!(reduced_array)
end

function naive_graham_scan!(point_array::Vector{<:AbstractPoint})
    # opti possibility: dont calculate the angle, but only the ordering is relevant
    L = length(point_array)
    if L <= 3
        return point_array
    end

    y_ = y(point_array);
    (miny, idx) = findmin(y_)
    point_array[idx], point_array[end] = point_array[end], point_array[idx]
    p0 = point_array[end]

    phi = Vector{Float64}(undef, L)
    for kk = 1:(L-1)
        # v = point_array[kk] - p0
        # phi[kk] = atan(v[2], v[1])
        phi[kk] = atan(point_array[kk].y - p0.y, point_array[kk].x - p0.x)
    end
    phi[end] = 1e308 # missing # missing is greater than anything else
    iU2S = sortperm(phi)
    point_array = point_array[iU2S]

    ElementType = eltype(point_array)
    stack = datastructs.Stack{ElementType}(L)

    push!(stack, p0)
    push!(stack, point_array[1])
    push!(stack, point_array[2])

    for kk = 3:(L-1)
        while nonleft_turn(
                datastructs.peek(stack, 1),
                datastructs.peek(stack, 0),
                point_array[kk] )
            pop!(stack)
        end
        push!(stack, point_array[kk])
    end
    ch = stack.data[1:stack.ptr]
    return ch
end

function nonleft_turn(p1::T, p2::T, p3::T)::Bool where {T<:AbstractPoint}
    # a = p2 - p1
    # b = p3 - p2
    # return cross2(a, b) <= 0

    a_x = p2.x - p1.x
    a_y = p2.y - p1.y

    b_x = p3.x - p2.x
    b_y = p3.y - p2.y
    return cross2(a_x, a_y, b_x, b_y) <= 0
end

function jarvis3!(point_array::Vector{<:AbstractPoint})
	L = length(point_array)
    if L <= 3
        return point_array
    end

    reduced_array = reduce_array_for_convex_hull!(point_array)
    L = length(reduced_array)
    (min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)

	stack = datastructs.Stack{Point{Float64}}(L)
	push!(stack, min_y)

	# util.swappop!(reduced_array, idx_min)
	util.swap!(reduced_array, idx_min, L)
	L = L - 1

    phi = Vector{Float64}(undef, L)
	ref_point = peek(stack, 0)

	a_x = zero(Float64)
	a_y = zero(Float64)
	b_x = zero(Float64)
	b_y = zero(Float64)

	for kk = 1:L
		# a_x = 1
		# a_y = 0
		b_x = reduced_array[kk].x - ref_point.x
		b_y = reduced_array[kk].y - ref_point.y
		sin_ = cross2(1.0, 0.0, b_x, b_y)
		cos_ = dot2(1.0, 0.0, b_x, b_y)
		phi[kk] = cos_ / sin_ # cotangent, the angle is irrelevant, only the ordering is important
		# phi[kk] = atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
	end
	(minphi, idx) = findmax(phi) # for cotangent the min angle is the max value

	push!(stack, reduced_array[idx])
	util.swappop!(reduced_array, idx)
	# pop!(phi) # must be popped, in order to not search in the last not used element in findmin
	# L = L - 1
	first = true
	while true
		prev_point = peek(stack, 1)
		ref_point = peek(stack, 0)
		for kk = 1:L
			a_x = ref_point.x - prev_point.x
			a_y = ref_point.y - prev_point.y
			b_x = reduced_array[kk].x - ref_point.x
			b_y = reduced_array[kk].y - ref_point.y
			sin_ = cross2(a_x, a_y, b_x, b_y)
			cos_ = dot2(a_x, a_y, b_x, b_y)
			phi[kk] = cos_ / sin_ # cotangent, the angle is irrelevant, only the ordering is important
			# phi[kk] =  sin_ / cos_
			# phi[kk] = atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
		end
		if first
			phi[isnan.(phi)] .= -Inf
			first = false
		end
		phi_view = view(phi, 1:L)
		(minphi, idx) = findmax(phi_view) # for cotangent the min angle is the max value
		if min_y == reduced_array[idx]
			break
		end
		push!(stack, reduced_array[idx])
		util.swap!(reduced_array, idx, L)
		# util.swappop!(reduced_array, idx)
		# pop!(phi) # must be popped, in order to not search in the last not used element in findmin
		L = L - 1
	end

	ch = [stack.data[1:stack.ptr]...]
	return ch

end

function jarvis2!(point_array::Vector{<:AbstractPoint})
	L = length(point_array)
    if L <= 3
        return point_array
    end

    reduced_array = reduce_array_for_convex_hull!(point_array)
    L = length(reduced_array)
    (min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)

	stack = datastructs.Stack{Point{Float64}}(div(L,2)+1)
	push!(stack, min_y)

	# util.swappop!(reduced_array, idx_min)
	util.swap!(reduced_array, idx_min, L)
	L = L - 1

    phi = Vector{Float64}(undef, L)
	ref_point = peek(stack, 0)

	a_x = zero(Float64)
	a_y = zero(Float64)
	b_x = zero(Float64)
	b_y = zero(Float64)

	for kk = 1:L
		# a_x = 1
		# a_y = 0
		b_x = reduced_array[kk].x - ref_point.x
		b_y = reduced_array[kk].y - ref_point.y
		sin_ = cross2(1.0, 0.0, b_x, b_y)
		cos_ = dot2(1.0, 0.0, b_x, b_y)
		phi[kk] = cos_ / sin_ # cotangent, the angle is irrelevant, only the ordering is important
		# phi[kk] = atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
	end
	(minphi, idx) = findmax(phi) # for cotangent the min angle is the max value

	push!(stack, reduced_array[idx])
	util.swappop!(reduced_array, idx)
	# pop!(phi) # must be popped, in order to not search in the last not used element in findmin
	# L = L - 1
	first = true
	while true
		prev_point = peek(stack, 1)
		ref_point = peek(stack, 0)
		for kk = 1:L
			a_x = ref_point.x - prev_point.x
			a_y = ref_point.y - prev_point.y
			b_x = reduced_array[kk].x - ref_point.x
			b_y = reduced_array[kk].y - ref_point.y
			sin_ = cross2(a_x, a_y, b_x, b_y)
			cos_ = dot2(a_x, a_y, b_x, b_y)
			phi[kk] = cos_ / sin_ # cotangent, the angle is irrelevant, only the ordering is important
			# phi[kk] =  sin_ / cos_
			# phi[kk] = atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
		end
		if first
			phi[isnan.(phi)] .= -Inf
			first = false
		end
		(minphi, idx) = findmax(phi) # for cotangent the min angle is the max value
		if min_y == reduced_array[idx]
			break
		end
		push!(stack, reduced_array[idx])
		util.swappop!(reduced_array, idx)
		pop!(phi) # must be popped, in order to not search in the last not used element in findmin
		L = L - 1
	end

	ch = [stack.data[1:stack.ptr]...]
	return ch

end

function jarvis_march!(point_array::Vector{<:AbstractPoint})

    L = length(point_array)
    if L <= 3
        return point_array
    end

    reduced_array = reduce_array_for_convex_hull!(point_array)
    L = length(reduced_array)
    (min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)
    # p_1 = reduced_array[1]
    # p_L = reduced_array[end]
    # reduced_array[1], reduced_array[idx_min] = reduced_array[idx_min], p_1
    # reduced_array[end], reduced_array[idx_max] = reduced_array[idx_max], p_L
    # at the index 1 is the point with min y coord, at L with max y

    stack_left = datastructs.Stack{Point{Float64}}(div(L,2)+1)
    stack_right = datastructs.Stack{Point{Float64}}(div(L,2)+1)

    push!(stack_right, min_y)
    push!(stack_left, max_y)

	# pop the min element, the angles shall not be computed for this point
	util.swappop!(reduced_array, idx_min)
	L = L - 1

    phi = Vector{Float64}(undef, L)

	idx = 0 # new local declaration
    # right segment
    while true
		ref_point = peek(stack_right, 0)
        for kk = 1:L
            phi[kk] = atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
        end

        (minphi, idx) = findmin(phi)
        if max_y == reduced_array[idx]
            break
        end
        push!(stack_right, reduced_array[idx])
		util.swappop!(reduced_array, idx)
		pop!(phi) # must be popped, in order to not search in the last not used element in findmin
		L = L - 1
    end

	# at idx it is known that tha max_y is there, but min_y is missing now in the array
	# put min_y at idx
	reduced_array[idx] = min_y
    # left segment
	while true
		ref_point = peek(stack_left, 0)
        for kk = 1:L
            phi[kk] = π + atan(reduced_array[kk].y - ref_point.y, reduced_array[kk].x - ref_point.x)
        end

        (minphi, idx) = findmin(phi)
        if min_y == reduced_array[idx]
            break
        end
        push!(stack_left, reduced_array[idx])
		util.swappop!(reduced_array, idx)
		pop!(phi) # must be popped, in order to not search in the last not used element in findmin
		L = L - 1
    end

	ch = [
	stack_right.data[1:stack_right.ptr]...
	stack_left.data[1:stack_left.ptr]...
	]
	return ch
end

# function quick_hull!(point_array::Vector{<:AbstractPoint})
#
# 	L = length(point_array)
# 	if L <= 3
# 		return point_array
# 	end
#
# 	reduced_array = reduce_array_for_convex_hull!(point_array)
# 	L = length(reduced_array)
# 	(min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)
#
#
#
# end

function incremental_hull!(point_array::Vector{<:AbstractPoint})

	L = length(point_array)
    if L <= 3
        return point_array
    end

    reduced_array = reduce_array_for_convex_hull!(point_array)

    L = length(reduced_array)
    (min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)

	aL = reduced_array[end]
	aLm1 = reduced_array[end-1]

	reduced_array[idx_min] = aL
	reduced_array[idx_max] = aLm1
	pop!(reduced_array)
	pop!(reduced_array)

	list = initial_hull!(reduced_array, max_y, min_y)
	L = length(reduced_array)-1 # p3 is in the array

	for kk = L:-1:1
		new_point = reduced_array[kk]
		bools = is_lefts(list, new_point)
		if all(bools)
			continue
		else
			(t1, t2) = tangent_points_to_hull(list, bools)
			datastructs.remove_between!(list, t1, t2)
			new_node = datastructs.new_node(list)
			new_node.data = new_point
			datastructs.insert_after!(list, t1, new_node)
		end
	end

	return datastructs.linked_list_to_array(list)
end

function tangent_points_to_hull(list, bools)
	node = list.node
	kk = 1
	L = length(list)
	if bools[1] == true
		while kk <= L && bools[kk] == true
			node = node.next
			kk += 1
		end
		t1 = node
		while kk <= L && bools[kk] == false
			node = node.next
			kk += 1
		end
		t2 = node
	else
		while kk <= L && bools[kk] == false
			node = node.next
			kk += 1
		end
		t2 = node
		while kk <= L && bools[kk] == true
			node = node.next
			kk += 1
		end
		t1 = node
	end
	return (t1, t2)
end

function initial_hull!(reduced_array, p1, p2)

	list = datastructs.LinkedList{Point{Float64}}()
	# n1 = datastructs.new_node(list)
	n2 = datastructs.new_node(list)
	n3 = datastructs.new_node(list)

	list.node.data = p1

	idx = length(reduced_array)
	while idx >= 1
		p3 = reduced_array[idx]
		if !is_on_line(p1, p2, p3)

			if is_left(p1, p2, p3)
				# n1.data = p1
				n2.data = p2
				n3.data = p3
			else
				# n1.data = p1
				n2.data = p3
				n3.data = p2
			end
			break
		end
		pop!(reduced_array)
		idx -= 1
	end

	datastructs.insert_after!(list, list.node, n2)
	datastructs.insert_after!(list, n2, n3)

	return list
end

function is_lefts(polygon_list::datastructs.LinkedList, new_point::AbstractPoint)

	bools = falses(polygon_list.count)
	node2 = polygon_list.node.next

	p1 = polygon_list.node.data
	p2 = node2.data

	for kk = 1:length(polygon_list)
		b = is_left_on(p1, p2, new_point)
		bools[kk] = b
		node2 = node2.next
		p1 = p2
		p2 = node2.data
	end

	return bools
end

#=
function incremental_hull!(point_array::Vector{<:AbstractPoint})

	L = length(point_array)
    if L <= 3
        return point_array
    end

    reduced_array = reduce_array_for_convex_hull!(point_array)

    L = length(reduced_array)
    (min_y, max_y, idx_min, idx_max) = extrema_y(reduced_array)

	list = datastructs.LinkedList{Point{Float64}}()
	list.node.data = min_y

	node = datastructs.new_node(list)
	node.data = max_y
	datastructs.insert_after!(list, list.node, node)

	aL = reduced_array[end]
	aLm1 = reduced_array[end-1]

	reduced_array[idx_min] = aL
	reduced_array[idx_max] = aLm1
	pop!(reduced_array)
	pop!(reduced_array)
	L = L - 2

	node = datastructs.new_node(list)

	# ensure that the 3. point is not on the line
	p1 = list.node.data
	p2 = list.node.next.data
	line_12 = Line(list.node.data, list.node.next.data)
	while true
		p = reduced_array[end]
		if is_on_line(line_12, p)
			pop!(reduced_array)
		else
			break
		end
	end
	node.data = reduced_array[end]
	pop!(reduced_array)
	datastructs.insert_after!(list, list.node.next, node)

	# enforce that the third point is om the left of 1 and 2
	p3 = list.node.next.next.data
	if is_left(line_12, p3)
		# fine, counter clockwise order
	else
		# swap 2 and 3
		p2 = list.node.next.data
		p3 = list.node.next.next.data
		list.node.next.data = p3
		list.node.next.next.data = p2
	end

	idx = length(reduced_array)
	while idx >= 1
		new_point = reduced_array[idx]
		bools = is_lefts(list, new_point)

		idx -= 1
		if all(bools)
			continue
		else
			(n1, n2) = find_transition_points(list, bools)
			new_node = datastructs.new_node(list)
			new_node.data = new_point
			if n1.data == n2.data
				datastructs.insert_after!(list, n2, new_node)
				continue
			end
			n2 = n2.next
			next_node = n1.next
			while next_node.data != n2.data
				nn = next_node.next
				datastructs.remove_element!(list, next_node)
				next_node = nn
			end
			datastructs.insert_after!(list, n1, new_node)
		end
	end

	return datastructs.linked_list_to_array(list)

end

function find_transition_points(list, bools)
	L = length(list)
	first_zero = 0
	last_zero = 0

	if bools[1] == false
		# backwards search for first zero
		kk = 1
		while bools[kk] == false
			kk = kk - 1
			if kk == 0
				kk = L
			end
		end
		first_zero = util.modulo_index(kk + 1, L)

		kk = 1
		while bools[kk] == false
			kk = kk + 1
		end
		last_zero = kk - 1

	else # bools[1] == true
		# forward search
		kk = 1
		while bools[kk] == true
			kk = kk + 1
		end
		first_zero = kk

		# kk = kk + 1
		while bools[kk] == false
			if kk == L
				kk = 0
			end
			kk = kk + 1
		end
		if kk == 1
			kk = L + 1
		end
		last_zero = util.modulo_index(kk - 1, L)
	end

	node = list.node
	kk = 1

	if first_zero > last_zero
		# search backwards

		### TO DO ###
		# what happens if first_zero == L
		kk = L
		while kk >= first_zero
			kk = kk - 1
			node = node.prev
		end
		n1 = node

		kk = 1
		node = list.node
		while kk < last_zero
			kk = kk + 1
			node = node.next
		end
		n2 = node
	else
		# search forward
		while kk < first_zero
			kk = kk + 1
			node = node.next
		end
		n1 = node
		while kk < last_zero
			kk = kk + 1
			node = node.next
		end
		n2 = node
	end

	return (n1, n2)
end

function is_lefts(polygon_list::datastructs.LinkedList, new_point::AbstractPoint)

	bools = falses(polygon_list.count)

	node2 = polygon_list.node.next

	p1 = polygon_list.node.data
	p2 = node2.data

	for kk = 1:length(polygon_list)
		b = is_left_on(p1, p2, new_point)
		bools[kk] = b
		node2 = node2.next
		p1 = p2
		p2 = node2.data
	end

	return bools
end

# function is_in_polygon_list(list::datastructs.LinkedList, new_point::AbstractPoint)
# 	#stupid linear algo
#
# 	bool = true
# 	p1 = list.node.data
# 	node2 = list.node.next
# 	p2 = node2.data
# 	line_12 = Line(p1, p2)
#
# 	while p2 != list.node
# 		if !is_left_on(line_12, new_point)
# 			bool = false
# 			break
# 		end
# 		p1 = p2
# 		node2 = node2.next
# 		p2 = node2.data
# 	end
#
# 	return bool
# end
=#

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Differential geometry.

# This struct represents the geometry on a point on a surface where a ray has hit.

struct DifferentialGeometry
    # Location and normal vector on the surface
    p::Point3
    n::Normal3
    T::Transformation
end

function DifferentialGeometry(p::Point3, n::Normal3, s::Vector3)
    n = normalize(n)
    s = normalize(s)
    return DifferentialGeometry(p, n, local_transformation(n, s))
end

"""
    flip_normal(n::Normal3, d::Vector3)

Inverts the normal vector `n` if the normal and the direction vector `d` are in the same
hemisphere. This is to always get normals pointing "out" of a surface regardless of whether
the ray comes from "inside" or "outside" the shape.

"""
flip_normal(n::Normal3, d::Vector3) = dot(n,d) < 0.0 ? n : -n

"""
    local_transformation(n::Normal3, s::Vector3)

Get the transformation which transforms a vector from the world coordinates to
surface-local coordinates, where the Z-axis is along the surface normal.

Makes an orthonormal basis of the normal vector `n`, a vector `s` tangent to the surface,
and their cross product. Both vectors are assumed to be normalized.

"""
function local_transformation(n::Normal3, s::Vector3)
    t = cross(n, s)
    M = zeros(4,4)
    for i = 1:3
        M[1,i] = s[i]
        M[2,i] = t[i]
        M[3,i] = n[i]
    end
    M[4,4] = 1.0

    return Transformation(M, M')
end

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using SIAN

println("Setting up the problem")

ode = @ODEmodel(
  s'(t) = mu - mu * s(t) - b0 * (1 + b1 * x1(t)) * i(t) * s(t) + g * r(t),
  i'(t) = b0 * (1 + b1 * x1(t)) * i(t) * s(t) - (nu + mu) * i(t),
  r'(t) = nu * i(t) - (mu + g) * r(t),
  x1'(t) = -M * x2(t),
  x2'(t) = M * x1(t),
  y1(t) = i(t),
  y2(t) = r(t)
)

res = identifiability_ode(ode, get_parameters(ode); p=0.99, p_mod=2^29 - 3, nthrds=1)

println(res)

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using Base.Broadcast: check_broadcast_indices, broadcast_indices

Base.@pure nullable_eltypestuple(a) = Tuple{eltype(eltype(a))}
Base.@pure nullable_eltypestuple(T::Type) = Tuple{Type{eltype(T)}}
Base.@pure nullable_eltypestuple(a, b...) =
    Tuple{nullable_eltypestuple(a).types..., nullable_eltypestuple(b...).types...}

Base.@pure function nullable_broadcast_eltype(f, As...)
    T = Core.Inference.return_type(f, nullable_eltypestuple(As...))
    T === Union{} ? Any : T
end

invoke_broadcast!{F, N}(f::F, dest, As::Vararg{DataValueArray, N}) =
    invoke(broadcast!, Tuple{F, AbstractArray, Vararg{AbstractArray, N}}, f, dest, As...)

function Base.broadcast{F}(f::F, As::DataValueArray...)
    T = nullable_broadcast_eltype(f, As...)
    dest = similar(DataValueArray{T}, broadcast_indices(As...))
    invoke_broadcast!(f, dest, As...)
end

function Base.broadcast!{F}(f::F, dest::DataValueArray, As::DataValueArray...)
    invoke_broadcast!(f, dest, As...)
end

# To fix ambiguity
function Base.broadcast!{F}(f::F, dest::DataValueArray)
    invoke_broadcast!(f, dest)
end


text:
using Plots

function gaussian(x; μ, σ)
    y= @. exp(-(x-μ)^2/2)/sqrt(2π)/sigma
end

function p(y;x,a,b)
    return gaussian(y, μ=a*sin(bx), σ=1)
end

function main()
    a=1.
    b=1.

    #function body
end

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function wilson_estimate(c::MolecularProperty, p, T)
    return wilson_estimate(p, T, c.ω, c.p_c, c.T_c)
end

"""
    wilson_estimate(p, T, ω, p_c, T_c)

Estimate K-values for a given acentric factor ω and pressure and temperature at current and critical conditions.

Reference [Vapor-Liquid Equilibrium. XI. A New Expression for the Excess Free Energy of Mixing by GM Wilson](https://doi.org/10.1021/ja01056a002)
"""
function wilson_estimate(p::R, T::R, ω::R, p_c::R, T_c::R) where R<:Real
    K = exp(5.37*(1.0 + ω)*(1.0 - T_c/T))*(p_c/p)
    return K::R
end

"""
    wilson_estimate!(K, properties, p, T)

Update a vector K in-place with K-values from `wilson_estimate`.
"""
function wilson_estimate!(K::AbstractVector{R}, props, p::R, T::R) where R<:Real
    @inbounds for i in eachindex(K)
        K[i] = wilson_estimate(props[i], p, T)::R
    end
end

"""
    wilson_estimate(properties, p, T)

Create vector of K-values that holds the `wilson_estimate` for each species.
"""
function wilson_estimate(mixture::MultiComponentMixture, p, T)
    K = zeros(number_of_components(mixture))
    wilson_estimate!(K, mixture.properties, p, T)
    return K
end

"""
    initial_guess_K(eos, cond)

Produce a plausible initial guess for K values for `eos` under current `cond`.
"""
function initial_guess_K(eos, cond)
    K = zeros(number_of_components(eos))
    initial_guess_K!(K, eos, cond)
    return K
end

"""
    initial_guess_K(eos, cond)

In-place version of `initial_guess_K`.
"""
initial_guess_K!(K, eos, cond) = wilson_estimate!(K, eos.mixture.properties, cond.p, cond.T)


text:
# This file is a part of RRRMC.jl. License is MIT: http://github.com/carlobaldassi/RRRMC.jl/LICENCE.md

module TLE

using ExtractMacro
using ..Interface
using ..Common

using ...RRRMC # this is silly but it's required for correct cross-linking in docstrings, apparently

export GraphTLE, GraphTopologicalLocalEntropy, TLEenergies

import ..Interface: energy, delta_energy, neighbors, allΔE,
                    update_cache!, delta_energy_residual,
                    cenergy, distances

mutable struct GraphTLE{M,γT,λT} <: DiscrGraph{Float64}
    N::Int
    Nk::Int
    neighb::Vector{Vector{Int}}
    neighb_jcs::Vector{Vector{Int}}
    neighb_Us::Vector{Vector{UnitRange{Int}}}
    cache1::LocalFields{Int} # replica-center
    cache2::LocalFields{Int} # topological
    ncache::Vector{Int}

    function GraphTLE(N::Integer, neighb::Vector{Vector{Int}}) where {M,γT,λT}
        isa(M, Int) || throw(ArgumentError("invalid parameter M, expected Int, given: $(typeof(M))"))
        M > 2 || throw(ArgumentError("M must be greater than 2, given: $M"))
        isa(γT, Float64) || throw(ArgumentError("invalid parameter γT, expected Float64, given: $(typeof(γT))"))
        N % (M + 1) == 0 || throw(ArgumentError("N must be divisible by M+1, given: N=$N M=$M"))
        Nk = N ÷ (M+1)
        length(neighb) == Nk || throw(ArgumentError("invalid neighb argument, expected length $Nk, given: $(length(neighb))"))
        for (i,ni) in enumerate(neighb)
            all(j->(1 ≤ j ≤ Nk), ni) || throw(ArgumentError("invalid neighb[$i] argument, all elements must be between 1 and $Nk"))
            i ∈ ni && throw(ArgumentError("invalid neighb[$i], contains itself"))
            lni = length(ni)
            for i1 in 1:lni, i2 in (i1+1):lni
                ni[i1] ≠ ni[i2] || throw(ArgumentError("duplicated element $i1 in neighb[$i]"))
            end
        end
        for (i,ni) in enumerate(neighb)
            for i1 in ni
                i ∈ neighb[i1] || throw(ArgumentError("neighb is not symmetrical, $i1 ∈ ∂$i but $i ∉ ∂$i1"))
            end
        end

        neighb_jcs = [[1 + (i1-1) * (M+1) for i1 in ni] for ni in neighb]
        neighb_Us = [[jinterval(i1, M) for i1 in ni] for ni in neighb]

        cache1 = LocalFields{Int}(N)
        cache2 = LocalFields{Int}(N)
        ncache = empty!(Array{Int}(undef, (M+1) * maximum(length.(neighb))))
        return new{M,γT,λT}(N, Nk, neighb, neighb_jcs, neighb_Us, cache1, cache2, ncache)
    end
end

@doc """
    GraphTLE{M,γT,λT}(N::Integer) <: DiscrGraph

An auxiliary `DiscrGraph` used to implement the interactions in the
Local Entropy Ensemble.

It is only useful when implementing other graph types; see [`GraphTopologicalLocalEntropy`](@ref).
""" -> GraphTLE{M,γT,λT}(N::Integer)

GraphTLE(X::GraphTLE{M,oldγ,oldλ}, newγ::Float64, newλ::Float64) where {M,oldγ,oldλ} = GraphTLE{M,newγ,newλ}(X.N, X.neighb)

# @generated function ΔElist{M,γT}(::Type{GraphTLE{M,γT}})
#     Expr(:tuple, ntuple(d->fk(2*(d - 1 - (M-1) >>> 0x1) - iseven(M), γT), M)...)
# end
# lstind(μ̄::Int, M::Int) = (μ̄ + M-1) >>> 0x1 + 1
#
# function getk{M,γT}(X::GraphTLE{M,γT}, μ̄::Int)
#     @inbounds k = ΔElist(GraphTLE{M,γT})[lstind(μ̄, M)]
#     return k
# end

function energy(X::GraphTLE{M,γT,λT}, C::Config) where {M,γT,λT}
    # @assert X.N == C.N
    @extract X : Nk neighb cache1 cache2
    @extract cache1 : lfields1=lfields lfields_last1=lfields_last
    @extract cache2 : lfields2=lfields lfields_last2=lfields_last
    @extract C : s

    n = 0
    j = 0
    for i = 1:Nk
        j += 1
        # @assert j == 1 + (i-1) * (M+1)
        jc = j
        σc = 2s[jc] - 1

        μ = 0
        for k = 2:(M+1)
            j += 1
            # @assert j == k + (i-1) * (M+1)
            σj = 2s[j] - 1
            lfields1[j] = σc * σj
            μ += σj
        end
        f = σc * μ
        lfields1[jc] = f
        n -= f
    end
    cache1.move_last = 0
    fill!(lfields_last1, 0.0)

    r = 0
    fill!(lfields2, 0.0)
    for (i1, ni) in enumerate(neighb)
        jc1 = 1 + (i1-1) * (M+1)
        σc1 = 2s[jc1] - 1
        μ = 0
        for i2 in ni
            jc2 = 1 + (i2-1) * (M+1)
            σc2 = 2s[jc2] - 1
            σc = σc1 * σc2

            for k = 2:(M+1)
                j1 = k + (i1-1) * (M+1)
                j2 = k + (i2-1) * (M+1)
                σj = (2s[j1] - 1) * (2s[j2] - 1)
                lfields2[j1] += σc * σj
                lfields2[j2] += σc * σj
                μ += σc2 * σj
            end
        end
        f = σc1 * μ
        lfields2[jc1] = 2f
        r -= f
    end
    @assert r % 2 == 0
    r ÷= 2
    lfields2 .÷= 2
    # @show n,r
    # @show n * γT, r * λT
    cache2.move_last = 0
    fill!(lfields_last2, 0.0)

    return n * γT + r * λT
end

function kinterval(move::Integer, M::Integer)
    j0 = move - ((move-1) % (M+1))
    j1 = j0 + M
    return j0:j1
end

function jinterval(i::Integer, M::Integer)
    a = (M+1) * (i-1) + 1
    b = (M+1) * i
    return a:b
end

function update_cache!(X::GraphTLE{M,γT,λT}, C::Config, move::Int) where {M,γT,λT}
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract C : N s
    @extract X : Nk neighb neighb_jcs neighb_Us cache1 cache2

    k = mod1(move, M+1)         # replica index
    i = (move-1) ÷ (M+1) + 1    # variable index
    jc = 1 + (i-1) * (M+1)      # corresponding reference
    if k ≠ 1
        jcs = neighb_jcs[i]
        # jcs = [1 + (i1-1) * (M+1) for i1 in neighb[i]]  # other neighboring references
        # js = [k + (i1-1) * (M+1) for i1 in neighb[i]]   # other neighboring variables
    else
        Ux = kinterval(move, M)                         # the whole slice
        # Us = [jinterval(i1, M) for i1 in neighb[i]]     # other neighboring slices
        Us = neighb_Us[i]
    end

    @extract cache1 : lfields1=lfields lfields_last1=lfields_last move_last1=move_last
    @extract cache2 : lfields2=lfields lfields_last2=lfields_last move_last2=move_last
    @assert move_last1 == move_last2
    if move_last1 == move
        @inbounds begin
            #for y in neighbors(X, move)
            if k ≠ 1
                lfields1[jc], lfields_last1[jc] = lfields_last1[jc], lfields1[jc]
                lfields1[move], lfields_last1[move] = lfields_last1[move], lfields1[move]
            else
                for y in Ux
                    lfields1[y], lfields_last1[y] = lfields_last1[y], lfields1[y]
                end
            end
            #lfields1[move] = -lfields1[move]
            #lfields_last1[move] = -lfields_last1[move]
        end

        @inbounds begin
            if k ≠ 1
                lfields2[jc], lfields_last2[jc] = lfields_last2[jc], lfields2[jc]
                lfields2[move], lfields_last2[move] = lfields_last2[move], lfields2[move]
                for jc1 in jcs
                    lfields2[jc1], lfields_last2[jc1] = lfields_last2[jc1], lfields2[jc1]
                    j1 = jc1 - 1 + k
                    lfields2[j1], lfields_last2[j1] = lfields_last2[j1], lfields2[j1]
                end
                # for j in js
                #     lfields2[j], lfields_last2[j] = lfields_last2[j], lfields2[j]
                # end
            else
                for y in Ux
                    lfields2[y], lfields_last2[y] = lfields_last2[y], lfields2[y]
                end
                for u in Us, y in u
                    lfields2[y], lfields_last2[y] = lfields_last2[y], lfields2[y]
                end
            end
        end
        return
    end

    @inbounds begin
        if k ≠ 1
            #jc = 1 + (i-1) * (M+1)
            σx = 2s[move] - 1
            σc = 2s[jc] - 1
            lfc = lfields1[jc]
            lfm = lfields1[move]
            lfields_last1[jc] = lfc
            lfields_last1[move] = lfm
            lfields1[jc] = lfc + 2 * (σc * σx)
            lfields1[move] = -lfm
        else
            for y in Ux
                lfy = lfields1[y]
                lfields_last1[y] = lfy
                lfields1[y] = -lfy
            end
        end
        #lfm = lfields1[move]
        #lfields_last1[move] = lfm
        #lfields1[move] = -lfm
    end
    cache1.move_last = move

    # lfields_bk1 = copy(lfields1)
    # energy(X, C)
    # lfields_bk1 ≠ lfields1 && @show move hcat(lfields1,lfields_bk1) find(lfields1 .≠ lfields_bk1)
    # @assert lfields_bk1 == lfields1

    @inbounds begin
        if k ≠ 1
            σx = 2s[move] - 1
            lfm = lfields2[move]
            lfields_last2[move] = lfm
            lfields2[move] = -lfm

            σc = 2s[jc] - 1
            lfc = lfields2[jc]
            lfields_last2[jc] = lfc
            # for (jc1,j1) in zip(jcs, js)
            for jc1 in jcs
                j1 = jc1 - 1 + k
                lfj1 = lfields2[j1]
                lfc1 = lfields2[jc1]
                lfields_last2[j1] = lfj1
                lfields_last2[jc1] = lfc1
                σj1 = 2s[j1] - 1
                σc1 = 2s[jc1] - 1
                σ = σx * σc * σj1 * σc1
                lfields2[j1] = lfj1 + 2 * σ
                lfields2[jc1] = lfc1 + 2 * σ
                lfc += 2σ
            end
            lfields2[jc] = lfc
        else
            for y in Ux
                lfy = lfields2[y]
                lfields_last2[y] = lfy
                lfields2[y] = -lfy
            end
            σc = 2s[jc] - 1
            for Ux1 in Us
                jc1 = first(Ux1)
                lfc1 = lfields2[jc1]
                lfields_last2[jc1] = lfc1
                σc1 = 2s[jc1] - 1
                for j1 = (jc1+1):last(Ux1)
                    lfj1 = lfields2[j1]
                    lfields_last2[j1] = lfj1
                    σj1 = 2s[j1] - 1
                    σx = 2s[jc+(j1-jc1)] - 1
                    σ = σx * σc * σj1 * σc1
                    lfields2[j1] = lfj1 + 2 * σ
                    lfc1 += 2 * σ
                end
                lfields2[jc1] = lfc1
            end
        end
    end
    cache2.move_last = move

    return
end

@inline function delta_energy(X::GraphTLE{M,γT,λT}, C::Config, move::Int) where {M,γT,λT}
    # @assert X.N == C.N
    # @assert 1 ≤ move ≤ C.N
    @extract X : cache1 cache2
    @extract cache1 : lfields1=lfields
    @extract cache2 : lfields2=lfields

    # @show 2γT * lfields1[move], 2λT * lfields2[move]
    @inbounds Δ = 2γT * lfields1[move] + 2λT * lfields2[move]
    return Δ
end

@inline function neighbors(X::GraphTLE{M}, j::Int) where {M}
    @extract X : neighb ncache
    k = mod1(j, M+1)         # replica index
    i = (j-1) ÷ (M+1) + 1    # variable index
    if k == 1
        resize!(ncache, M + (M+1) * length(neighb[i]))
        ncache[1:M] = j + (1:M)
        for (n,i1) in enumerate(neighb[i])
            ncache[M + jinterval(n, M)] = jinterval(i1, M)
        end
        # return (j+1):(j+M)
    else
        resize!(ncache, 1 + 2 * length(neighb[i]))
        ncache[1] = j - k + 1
        for (n,i1) in enumerate(neighb[i])
            ncache[2n] = 1 + (i1-1) * (M+1)
            ncache[2n+1] = k + (i1-1) * (M+1)
        end
    end
    return ncache
end

function allΔE(X::GraphTLE{M,γT,λT}) where {M,γT,λT}
    @extract X : neighb
    ΔE1 = iseven(M) ? insert!([4*d*γT for d = 0:(M÷2)], 2, 2γT) :
                      [2*(2d-1)*γT for d = 1:((M+1)÷2)]

    mn = M * maximum(length.(neighb))
    ΔE2 = [2*d*λT for d = -mn:mn]
    # @show ΔE1
    # @show ΔE2

    ΔE = sort!(unique(abs.(ΔE1 .+ ΔE2')))
    return tuple(ΔE...)
end

# Replicate an existsing graph

mutable struct GraphTopologicalLocalEntropy{M,γT,λT,G<:AbstractGraph} <: DoubleGraph{DiscrGraph{Float64},Float64}
    N::Int
    Nk::Int
    X0::GraphTLE{M,γT,λT}
    Xc::G
    X1::Vector{G}
    Cc::Config
    C1::Vector{Config}
    function GraphTopologicalLocalEntropy(N::Integer, neighb::Vector{Vector{Int}}, g0::G, Gconstr, args...) where {M,γT,λT,G}
        X0 = GraphTLE{M,γT,λT}(N, neighb)
        Nk = X0.Nk
        X1 = Array{G}(undef, M)
        Xc = g0
        for k = 1:M
            X1[k] = Gconstr(args...)
        end
        Cc = Config(Nk, init=false)
        C1 = [Config(Nk, init=false) for k = 1:M]
        return new{M,γT,λT,G}(N, Nk, X0, Xc, X1, Cc, C1)
    end
end

"""
    GraphTopologicalLocalEntropy(N::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64, Gconstr, args...) <: DoubleGraph

A `DoubleGraph` which implements a "Local Entropy" model, given any other Ising model previously defined.
This simulates a replicated system in which each replica interacts with an extra "reference"
configuration.
This kind of graph can be simulated efficiently with [`rrrMC`](@ref).

`N` is the number of spins, `M` the number of replicas, `γ` the interaction strength of the replicas with the reference,
`λ` is the topological interaction strenght,
`β` the inverse temperature. `GConstr` is the (original) graph constructor, and `args` the
arguments to the contructor.

This is similar to [`GraphLocalEntropy`](@ref RRRMC.GraphLocalEntropy), but with the additional topological interaction.

See also [`TLEenergies`](@ref) and [`cenergy`](@ref).
"""
function GraphTopologicalLocalEntropy(Nk::Integer, M::Integer, γ::Float64, λ::Float64, β::Float64, Gconstr, args...)
    g0 = Gconstr(args...)
    getN(g0) == Nk || throw(ArgumentError("expected $Nk spins, got $(getN(g0))"))
    G = typeof(g0)
    neighb = [Int[neighbors(g0, i)...] for i = 1:Nk]
    return GraphTopologicalLocalEntropy{M,γ/β,λ/β,G}(Nk * (M+1), neighb, g0, Gconstr, args...)
end

function update_cache!(X::GraphTopologicalLocalEntropy{M}, C::Config, move::Int) where {M}
    @extract X : X0 Xc X1 Cc C1
    k = mod1(move, M+1)
    i = (move - 1) ÷ (M+1) + 1

    if k == 1
        #spinflip!(Xc, Cc, i)
        spinflip!(Cc, i) # Xc cache not updated!
    else
        spinflip!(X1[k-1], C1[k-1], i)
    end

    update_cache!(X0, C, move)
end

function energy(X::GraphTopologicalLocalEntropy{M}, C::Config) where {M}
    # @assert X.N == C.N
    @extract X : Nk X0 Xc X1 Cc C1
    @extract C : s

    E = energy(X0, C)

    s1 = Cc.s
    for (i,j) = enumerate(1:(M+1):(1 + (M+1) * (Nk-1)))
        s1[i] = s[j]
    end
    energy(Xc, Cc)

    for k = 1:M
        s1 = C1[k].s
        for (i,j) = enumerate((k+1):(M+1):(k+1 + (M+1) * (Nk-1)))
            s1[i] = s[j]
        end
        E += energy(X1[k], C1[k])
    end

    return E
end

"""
    TLEenergies(X::GraphTopologicalLocalEntropy)

Returns a Vector with the individual energy (as defined by the original model)
of each replica in a [`GraphTopologicalLocalEntropy`](@ref) graph.
"""
function TLEenergies(X::GraphTopologicalLocalEntropy{M}) where {M}
    @extract X : X1 C1

    Es = zeros(M)

    for k = 1:M
        Es[k] = energy(X1[k], C1[k])
    end

    return Es
end

function cenergy(X::GraphTopologicalLocalEntropy{M}) where {M}
    @extract X : Xc Cc
    return energy(Xc, Cc)
end

function delta_energy_residual(X::GraphTopologicalLocalEntropy{M}, C::Config, move::Int) where {M}
    @extract X : X1 C1

    k = mod1(move, M+1)
    k == 1 && return 0.0

    i = (move - 1) ÷ (M+1) + 1

    return delta_energy(X1[k-1], C1[k-1], i)
end

function delta_energy(X::GraphTopologicalLocalEntropy, C::Config, move::Int)
    return delta_energy(X.X0, C, move) +
           delta_energy_residual(X, C, move)
end

# This is type-unstable and inefficient. On the other hand, it is
# basically only written for testing purposes...
function neighbors(X::GraphTopologicalLocalEntropy{M}, i::Int) where {M}
    @extract X : X0 X1
    jts = neighbors(X0, i)

    k = mod1(i, M+1)

    k == 1 && return jts

    i1 = (i - 1) ÷ (M+1) + 1

    is = neighbors(X1[k-1], i1)

    return tuple(jts..., map(i->((i-1)*(M+1)+k), is)...)
end

function distances(X::GraphTopologicalLocalEntropy{M}) where {M}
    @extract X : C1
    # @extract C : s
    # @assert Nk == C.N ÷ (M+1)

    # rs = BitVector[s[(k-1)*(M+1) + (2:(M+1):end)] for k = 1:M]

    # return [sum(s1 .⊻ s2) for s1 in rs, s2 in rs]
    return [sum(C1[k1].s .⊻ C1[k2].s) for k1 = 1:M, k2 = 1:M]
end

end # module


text:
"""Returns preallocated variables of different size that derive from u."""
function preallocate_u_vars()
    # with full halo
    du = zeros(Numtype,nux+2*halo,nuy+2*halo)
    u0 = zero(du)
    u1 = zero(du)

    # derivative: one less in x or y direction
    dudx = zeros(Numtype,nux+2*halo-1,nuy+2*halo)
    dudy = zeros(Numtype,nux+2*halo,nuy+2*halo-1)

    return du,u0,u1,dudx,dudy
end

"""Returns preallocated variables of different size that derive from v."""
function preallocate_v_vars()
    # with full halo
    dv = zeros(Numtype,nvx+2*halo,nvy+2*halo)
    v0 = zero(dv)
    v1 = zero(dv)

    # derivative: one less in x or y direction
    dvdx = zeros(Numtype,nvx+2*halo-1,nvy+2*halo)
    dvdy = zeros(Numtype,nvx+2*halo,nvy+2*halo-1)

    return dv,v0,v1,dvdx,dvdy
end

"""Returns preallocated variables of different size that derive from η."""
function preallocate_η_vars()
    # with full halo
    dη = zeros(Numtype,nx+2*haloη,ny+2*haloη)
    η0 = zero(dη)
    η1 = zero(dη)
    h = zero(dη)

    return dη,η0,η1,h
end

"""Returns preallocated variables of different size that will be used in the contuinity equation."""
function preallocate_continuity()
    # interpolation: one less in x-direction
    h_u = zeros(Numtype,nx+2*haloη-1,ny+2*haloη)
    U = zero(h_u)

    # interpolation: one less in y-direction
    h_v = zeros(Numtype,nx+2*haloη,ny+2*haloη-1)
    V = zero(h_v)

    # Derivatives: two less in x- or y-direction
    dUdx = zeros(Numtype,nx+2*haloη-2,ny+2*haloη)
    dVdy = zeros(Numtype,nx+2*haloη,ny+2*haloη-2)

    return h_u,U,h_v,V,dUdx,dVdy
end

"""Returns preallocated variables of different size that will be used for PV and the Sadourny adv scheme."""
function preallocate_Sadourny()

    # interpolation from h: ones less in both directions
    h_q = zeros(Numtype,nx+2*haloη-1,ny+2*haloη-1)
    q = zero(h_q)

    # two less in x direction, one less in y
    q_v = zeros(Numtype,nx+2*haloη-2,ny+2*haloη-1)
    U_v = zero(q_v)

    # two less in y direction, one less in x
    q_u = zeros(Numtype,nx+2*haloη-1,ny+2*haloη-2)
    V_u = zero(q_u)

    #
    qhu = zeros(Numtype,nvx,nvy)
    qhv = zeros(Numtype,nux,nuy)

    return h_q,q,q_v,qhu,U_v,q_u,qhv,V_u
end

"""Returns preallocated variables of different size for the PV linear combinations in the Arakawa&Hsu advection scheme."""
function preallocate_ArakawaHsu()
    qα = zeros(Numtype,nx+2*haloη-2,ny+2*haloη-2)
    qβ = zeros(Numtype,nx+2*haloη-1,ny+2*haloη-2)
    qγ = zeros(Numtype,nx+2*haloη-1,ny+2*haloη-2)
    qδ = zeros(Numtype,nx+2*haloη-2,ny+2*haloη-2)
    return qα,qβ,qγ,qδ
end

"""Returns preallocated variables of different size that will be used for the Bernoulli potential."""
function preallocate_Bernoulli()
    # with full halo
    u² = zeros(Numtype,nux+2*halo,nuy+2*halo)
    v² = zeros(Numtype,nvx+2*halo,nvy+2*halo)

    KEu = zeros(Numtype,nux+2*halo-1,nuy+2*halo)
    KEv = zeros(Numtype,nvx+2*halo,nvy+2*halo-1)

    p = zeros(Numtype,nx+2*haloη,ny+2*haloη)
    dpdx = zeros(Numtype,nx+2*haloη-1,ny+2*haloη)
    dpdy = zeros(Numtype,nx+2*haloη,ny+2*haloη-1)

    return u²,v²,KEu,KEv,p,dpdx,dpdy
end

"""Returns preallocated variables of different size that will be used for the bottom drag term."""
function preallocate_bottomdrag()
    sqrtKE = zeros(Numtype,nx+2*haloη,ny+2*haloη)
    sqrtKE_u = zeros(Numtype,nx+2*haloη-1,ny+2*haloη)
    sqrtKE_v = zeros(Numtype,nx+2*haloη,ny+2*haloη-1)

    Bu = zeros(Numtype,nx+2*haloη-1,ny+2*haloη)
    Bv = zeros(Numtype,nx+2*haloη,ny+2*haloη-1)

    return sqrtKE,sqrtKE_u,sqrtKE_v,Bu,Bv
end

"""Returns preallocated variables of different size that will be used for the momentum diffusion term."""
function preallocate_Laplace()
    # two less in both directions
    Lu = zeros(Numtype,nux+2*halo-2,nuy+2*halo-2)
    Lv = zeros(Numtype,nvx+2*halo-2,nvy+2*halo-2)

    # Derivatives of Lu,Lv
    dLudx = zeros(Numtype,nux+2*halo-3,nuy+2*halo-2)
    dLudy = zeros(Numtype,nux+2*halo-2,nuy+2*halo-3)
    dLvdx = zeros(Numtype,nvx+2*halo-3,nvy+2*halo-2)
    dLvdy = zeros(Numtype,nvx+2*halo-2,nvy+2*halo-3)

    return Lu,Lv,dLudx,dLudy,dLvdx,dLvdy
end

"""Returns preallocated variables of different size that will be used in the contuinity equation."""
function preallocate_Smagorinsky()
    # on the η-grid including halo
    DT = zeros(Numtype,nx+2*haloη,ny+2*haloη)
    DS = zero(DT)
    νSmag = zero(DT)

    # DS_q has same size as dvdx
    DS_q = zeros(Numtype,nvx+2*halo-1,nvy+2*halo)

    # one less in both directions, the q-grid
    νSmag_q = zeros(Numtype,nx+2*haloη-1,ny+2*haloη-1)
    S12 = zero(νSmag_q)
    S21 = zero(νSmag_q)

    S11 = zeros(Numtype,nux+2*halo-3,nuy+2*halo-2)
    S22 = zeros(Numtype,nvx+2*halo-2,nvy+2*halo-3)

    LLu1 = zeros(Numtype,nux+2*halo-4,nuy+2*halo-2)
    LLu2 = zeros(Numtype,nx+1,ny)

    LLv1 = zeros(Numtype,nx,ny+1)
    LLv2 = zeros(Numtype,nvx+2*halo-2,nvy+2*halo-4)

    return DT,DS,DS_q,νSmag,νSmag_q,S11,S12,S21,S22,LLu1,LLu2,LLv1,LLv2
end

"""Returns preallocated variables of different size that will be used for the semi-Lagrangian tracer advection."""
function preallocate_semiLagrange()
    xd = zeros(Numtype,nx,ny)
    yd = zero(xd)

    um = zeros(Numtype,nux+2*halo,nuy+2*halo)
    vm = zeros(Numtype,nvx+2*halo,nvy+2*halo)

    # u on T-grid one less in x-direction
    u_T = zeros(Numtype,nux+2*halo-1,nuy+2*halo)
    um_T = zero(u_T)

    # v on T-grid one less in y-direction
    v_T = zeros(Numtype,nvx+2*halo,nvy+2*halo-1)
    vm_T = zero(v_T)

    uinterp = zeros(Numtype,nx,ny)
    vinterp = zero(uinterp)

    ssti = zeros(Numtype,nx+2*halosstx,ny+2*halossty)

    return xd,yd,um,vm,u_T,um_T,v_T,vm_T,uinterp,vinterp,ssti
end

max_stars_repo_name: eschnett/juls

max_issues_repo_name: eschnett/juls

max_forks_repo_name: eschnett/juls


text:
function sortrows(a)

n = size(a,1);
b = sortslices([a collect(1:n)], dims=1);
ind = Int.(b[:,end]);
b = b[:,1:end-1];

return b,ind

end

max_stars_repo_name: glwagner/Exasim

max_issues_repo_name: glwagner/Exasim

max_forks_repo_name: glwagner/Exasim


text:
"""
    χ²(sem_fit::SemFit)

Return the χ² value.
"""
function χ² end

############################################################################################
# Single Models
############################################################################################

# SemFit splices loss functions ------------------------------------------------------------
χ²(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: AbstractSemSingle, O}) =
    χ²(
        sem_fit,
        sem_fit.model.observed,
        sem_fit.model.imply,
        sem_fit.model.optimizer,
        sem_fit.model.loss.functions...
        )

# RAM + SemML
χ²(sem_fit::SemFit, observed, imp::Union{RAM, RAMSymbolic}, optimizer, loss_ml::SemML) =
    (n_obs(sem_fit)-1)*(sem_fit.minimum - logdet(observed.obs_cov) - observed.n_man)

# bollen, p. 115, only correct for GLS weight matrix
χ²(sem_fit::SemFit, observed, imp::Union{RAM, RAMSymbolic}, optimizer, loss_ml::SemWLS) =
    (n_obs(sem_fit)-1)*sem_fit.minimum

# FIML
function χ²(sem_fit::SemFit, observed::SemObservedMissing, imp, optimizer, loss_ml::SemFIML)
    ll_H0 = minus2ll(sem_fit)
    ll_H1 = minus2ll(observed)
    chi2 = ll_H0 - ll_H1
    return chi2
end

############################################################################################
# Collections
############################################################################################

# SemFit splices loss functions ------------------------------------------------------------
χ²(sem_fit::SemFit{Mi, So, St, Mo, O} where {Mi, So, St, Mo <: SemEnsemble, O}) =
    χ²(
        sem_fit,
        sem_fit.model,
        sem_fit.model.sems[1].loss.functions[1]
        )

function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemWLS}
    check_ensemble_length(model)
    check_lossfun_types(model, L)
    return (sum(n_obs.(model.sems))-1)*sem_fit.minimum
end

function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemML}
    check_ensemble_length(model)
    check_lossfun_types(model, L)
    F_G = sem_fit.minimum
    F_G -= sum([w*(logdet(m.observed.obs_cov) + m.observed.n_man) for (w, m) in zip(model.weights, model.sems)])
    return (sum(n_obs.(model.sems))-1)*F_G
end

function χ²(sem_fit::SemFit, model::SemEnsemble, lossfun::L) where {L <: SemFIML}
    check_ensemble_length(model)
    check_lossfun_types(model, L)

    ll_H0 = minus2ll(sem_fit)
    ll_H1 = sum(minus2ll.(observed.(models(model))))
    chi2 = ll_H0 - ll_H1

    return chi2
end

function check_ensemble_length(model)
    for sem in model.sems
        if length(sem.loss.functions) > 1
            @error "A model for one of the groups contains multiple loss functions."
        end
    end
end

function check_lossfun_types(model, type)
    for sem in model.sems
        for lossfun in sem.loss.functions
            if !isa(lossfun, type)
                @error "Your model(s) contain multiple lossfunctions with differing types."
            end
        end
    end
end


text:
using Printf
using Statistics

const LOGDIR = joinpath(@__DIR__, "../../logs/rmp")
const RMP_INSTANCES = readdir(joinpath(@__DIR__,  "../../dat/rmp"))
const SOLVERS = ["CPLEX", "Gurobi", "Mosek", "Tulip", "Tulip_UBA"]
const SUFFIX = Dict{String, String}(
    "CPLEX" => ".cpx",
    "Gurobi" => ".grb",
    "Mosek" => ".msk",
    "Tulip" => ".tlp",
    "Tulip_UBA" => ".tlpu"
)
const SOLVERS_SHORT = Dict{String, String}(
    "CPLEX" => "CPX",
    "Gurobi" => "GRB",
    "Mosek" => "MSK",
    "Tulip" => "TLP",
    "Tulip_UBA" => "TLP*"
)

function parse_logname(finst)
    if finst[1:3] == "DER"
        fname, T, R, cg = split(finst, r"_|\.")[1:4]
        fname = "DER_$T"
    else
        fname, R, cg = split(finst, r"_|\.")[1:3]
    end
    R = parse(Int, R)
    cg = parse(Int, cg)
    return fname, R, cg
end

# Get number of CG iterations for each instance
function extract_cg_iter(instances)
    CG_ITER = Dict{Tuple{String, Int}, Int}()

    for finst in instances
        fname, R, cg = parse_logname(finst)
        cg_ = get(CG_ITER, (fname, R), 0)
        if cg > cg_
            CG_ITER[fname, R] = cg
        end
    end
    return CG_ITER
end

function parse_res(flog::String, res=Dict{String, Any}())
    open(flog, "r") do f
        for ln in eachline(f)
            if occursin("Solve time:", ln)
                res["time"] = parse(Float64, split(ln)[end])
                continue
            end
            if occursin("Barrier iterations:", ln)
                res["iter"] = parse(Int, split(ln)[end])
                continue
            end
        end
    end
    return res
end

function parse_all(instances, logdir)
    res_all = Dict()

    for finst in instances

        fname, R, cg = parse_logname(finst)

        res_all[fname, R, cg] = Dict()
        for solver in SOLVERS
            res_all[fname, R, cg][solver] = Dict()
            # Parse log file
            flog = joinpath(logdir, lowercase(solver), finst*SUFFIX[solver])
            r = parse_res(flog)

            # Check result status
            t = get(r, "time", Inf)
            iter = get(r, "iter", Inf)

            res_all[fname, R, cg][solver]["time"] = t
            res_all[fname, R, cg][solver]["iter"] = iter
        end
    end

    return res_all
end

function print_table(res, last_only::Bool=false)
    # Header line
    @printf "\\toprule\n"
    @printf "%-16s &   \$R\$" "Problem"
    last_only || @printf " & CG"
    for solver in SOLVERS
        @printf " & \\multicolumn{2}{c}{%s}" solver
    end
    @printf "\\\\\n"

    CG_ITER = extract_cg_iter(RMP_INSTANCES)

    # Individual times
    @printf "\\midrule\n"
    KEYS = sort(collect(keys(res)))
    for (fname, R, cg) in KEYS

        if last_only
            cg == CG_ITER[fname, R] || continue
        end

        @printf "%-16s & %5d" replace(fname, "_" => "-") R
        last_only || @printf " & %2d" cg
        finst = "$(fname)_$(R)_$(cg).mps"

        # Time
        tmin = minimum([res[fname, R, cg][solver]["time"] for solver in SOLVERS])
        for solver in SOLVERS
            # Parse log file
            flog = joinpath(LOGDIR, lowercase(solver), finst*SUFFIX[solver])
            r = parse_res(flog)

            # Check result status
            t = res[fname, R, cg][solver]["time"]
            iter = res[fname, R, cg][solver]["iter"]
            titer = t / iter
            if t == tmin

                @printf " & \\textbf{%6.1f} & %3d & %6.3f" t iter titer
            else
                @printf " & %6.1f & %3d & %6.3f" t iter titer
            end
        end

        @printf "\\\\\n"
    end
    @printf "\\bottomrule\n"
    return nothing
end

gmean(u, δ) = exp(mean(log.(u .+ δ))) - δ

if abspath(PROGRAM_FILE) == @__FILE__
    res = parse_all(RMP_INSTANCES, LOGDIR)
    print_table(res, true)
    return nothing
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
# Several strategies for missing data handling were attempted.
# 1. Use `missing`s as expected. This turned out to be problematic for type-stability.
# 2. Sentinel values (NaNs). Also problematic for type-stability because Zygote.
# 3. (The adopted strategy) - replace missings with arbitrary observations and _large_
#   observation noises. While not optimal, type-stability is preserved inside the
#   performance-sensitive code.
#
# In an ideal world, strategy 1 would work. Unfortunately Zygote isn't up to it yet.

function AbstractGPs.logpdf(
    model::LGSSM, y::AbstractVector{Union{Missing, T}},
) where {T<:Union{<:AbstractVector, <:Real}}
    model_with_missings, y_filled_in = transform_model_and_obs(model, y)
    return logpdf(model_with_missings, y_filled_in) + _logpdf_volume_compensation(y, model)
end

function _filter(model::LGSSM, y::AbstractVector{Union{Missing, T}}) where {T}
    model_with_missings, y_filled_in = transform_model_and_obs(model, y)
    return _filter(model_with_missings, y_filled_in)
end

function posterior(model::LGSSM, y::AbstractVector{Union{Missing, T}}) where {T}
    model_with_missings, y_filled_in = transform_model_and_obs(model, y)
    return posterior(model_with_missings, y_filled_in)
end

function transform_model_and_obs(
    model::LGSSM, y::AbstractVector{<:Union{Missing, T}},
) where {T<:Union{<:AbstractVector, <:Real}}
    Σs_filled_in, y_filled_in = fill_in_missings(
        zygote_friendly_map(noise_cov, emissions(model)), y,
    )
    model_with_missings = replace_observation_noise_cov(model, Σs_filled_in)
    return model_with_missings, y_filled_in
end

function replace_observation_noise_cov(model::LGSSM, Σs_new::AbstractVector)
    return LGSSM(transitions(model), replace_noise_cov(emissions(model), Σs_new))
end

function replace_noise_cov(emissions::StructArray{T}, Qs_new) where {T<:AbstractLGC}
    return StructArray{T}((emissions.A, emissions.a, Qs_new))
end

_large_var_const() = 1e15

function _logpdf_volume_compensation(y, model)
    emissions = model.emissions
    y_obs_count = sum(n -> y[n] === missing ? dim_out(emissions[n]) : 0, eachindex(y))
    return y_obs_count * log(2π * _large_var_const()) / 2
end

function _logpdf_volume_compensation(y::AbstractVector{<:Union{Missing, <:Real}})
    return count(ismissing, y) * log(2π * _large_var_const()) / 2
end

function Zygote._pullback(::AContext, ::typeof(_logpdf_volume_compensation), y)
    return _logpdf_volume_compensation(y), nograd_pullback
end

function fill_in_missings(Σs::Vector, y::AbstractVector{Union{Missing, T}}) where {T}
    return _fill_in_missings(Σs, y)
end

function _fill_in_missings(Σs::Vector, y::AbstractVector{Union{Missing, T}}) where {T}

    # Fill in observation covariance matrices with very large values.
    Σs_filled_in = map(eachindex(y)) do n
        (y[n] === missing) ? build_large_var(Σs[n]) : Σs[n]
    end

    # Fill in missing y's with zero vectors (any other choice would also have been fine).
    y_filled_in = Vector{T}(undef, length(y))
    map!(
        n -> y[n] === missing ? get_zero(size(Σs[n], 1), T) : y[n],
        y_filled_in,
        eachindex(y),
    )
    return Σs_filled_in, y_filled_in
end

function fill_in_missings(Σ::Diagonal, y::AbstractVector{<:Union{Missing, <:Real}})
    Σ_diag_filled, y_filled = fill_in_missings(Zygote.literal_getfield(Σ, Val(:diag)), y)
    return Diagonal(Σ_diag_filled), y_filled
end

# We need to densify anyway, might as well do it here and save having to implement the
# pullback twice.
function fill_in_missings(Σs::Fill, y::AbstractVector{Union{Missing, T}}) where {T}
    return fill_in_missings(collect(Σs), y)
end

fill_in_missings(Σ::Diagonal, y::AbstractVector{<:Real}) = (Σ, y)

function Zygote._pullback(
    ::AContext,
    ::typeof(_fill_in_missings),
    Σs::Vector,
    y::AbstractVector{Union{T, Missing}},
) where {T}
    pullback_fill_in_missings(Δ::Nothing) = nothing
    function pullback_fill_in_missings(Δ)
        ΔΣs_filled_in = Δ[1]
        Δy_filled_in = Δ[2]

        # The cotangent of a `Missing` doesn't make sense, so should be a `NoTangent`.
        Δy = if Δy_filled_in === nothing
            nothing
        else
            Δy = Vector{Union{eltype(Δy_filled_in), NoTangent}}(undef, length(y))
            map!(
                n -> y[n] === missing ? NoTangent() : Δy_filled_in[n],
                Δy, eachindex(y),
            )
            Δy
        end

        # Fill in missing locations with zeros. Opting for type-stability to keep things
        # simple.
        ΔΣs = map(
            n -> y[n] === missing ? zero(Σs[n]) : ΔΣs_filled_in[n],
            eachindex(y),
        )

        # return nothing, ΔΣs, Δy
        return nothing, ΔΣs, Δy
    end
    return fill_in_missings(Σs, y), pullback_fill_in_missings
end

get_zero(D::Int, ::Type{Vector{T}}) where {T} = zeros(T, D)

get_zero(::Int, ::Type{T}) where {T<:SVector} = zeros(T)

get_zero(::Int, ::Type{T}) where {T<:Real} = zero(T)

build_large_var(S::T) where {T<:Matrix} = T(_large_var_const() * I, size(S))

build_large_var(::T) where {T<:SMatrix} = T(_large_var_const() * I)

build_large_var(S::T) where {T<:Diagonal} = T(fill(_large_var_const(), length(diag(S))))

build_large_var(::T) where {T<:Real} = T(_large_var_const())

ChainRulesCore.@non_differentiable build_large_var(::Any)


text:
"""
```
measurement{T<:AbstractFloat}(m::RealBondMkup{T}, TTT::Matrix{T}, RRR::Matrix{T},
                              CCC::Vector{T})
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::RealBondMkup{T}, TTT::Matrix{T},
                     TTT_jump::Matrix{T},
                     RRR::Matrix{T}, CCC::Vector{T}, GDPeqn::AbstractArray{T}) where {T<:AbstractFloat}

    #@info "enter"
    endo        = m.endogenous_states
    endo_unnorm = m.endogenous_states_unnormalized
    endo_new    = m.model_states_augmented
    exo         = m.exogenous_shocks
    obs         = m.observables

    _n_model_states = n_model_states(m)
    _n_model_states_aug = _n_model_states + length(endo_new)
    _n_states       = n_backward_looking_states(m)
    _n_jumps        = n_jumps(m)

    _n_observables = n_observables(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    # Load in parameters, steady-state parameters, and grids
    R::Float64     = m[:R].value

    ZZ = zeros(_n_observables, _n_model_states_aug)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    # GDP
    ZZ[obs[:obs_gdp], 1:_n_states]     = GDPeqn
    ZZ[obs[:obs_gdp], _n_states+1:_n_states + _n_jumps] = zeros(_n_jumps)
    #@show size(ZZ[obs[:obs_gdp], _n_states + _n_jumps + 1:end])
    ZZ[obs[:obs_gdp], first(_n_states + _n_jumps + 1:end)] = -1

    # Inflation
    ZZ[obs[:obs_corepce], first(endo[:π′_t])] = 1.

    # Nominal FFR
    ZZ[obs[:obs_nominalrate], first(endo[:i′_t])] = 1.0 / R

    # Measurement error
    EE[obs[:obs_gdp], obs[:obs_gdp]] = m[:e_y]

    # Variance of innovations
    QQ[exo[:z_sh], exo[:z_sh]] = m[:σ_z]^2
    QQ[exo[:mon_sh], exo[:mon_sh]] = m[:σ_mon]^2
    QQ[exo[:mkp_sh], exo[:mkp_sh]] = m[:σ_mkp]^2

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        @show CCC[CCC.!=0]
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    return Measurement(ZZ, DD, QQ, EE)
end

function construct_GDPfn_realbondmkp(nx::Int, ns::Int,
                                  xgrid_total::Vector{Float64}, sgrid::Vector{Float64},
                                  xwts::Vector{Float64}, swts::Vector{Float64},
                                  γ::Float64, ν::Float64, abar::Float64,
                                  R::Float64, aborrow::Float64,
                                  μ::Vector{Float64}, c::Vector{Float64}, η::Vector{Float64},
                                  ell::Vector{Float64}, χss::Vector{Float64})
    dGDP_dMU  = zeros(1,nx*ns)
    dGDP_dELL = zeros(1,nx*ns)
    dGDP_dWW  = 0.0
    dGDP_dRR  = 0.0
    dGDP_dTT  = 0.0
    dGDP_dZ   = 0.0

    unc       = zeros(nx*ns)
    chipW, chipR, chipX = construct_chip_realbond(xgrid_total, γ, ν, aborrow,
                                                  abar, R, χss)
    for ix =1:nx
        for is=1:ns
            i = ix + nx*(is-1)
            unc[i] = ((ell[i]^(-1/γ))<=χss[i]) # =1 if unconstrained
            dGDP_dMU[i]  = xwts[ix]*swts[is]*sgrid[is]*η[i]
            dGDP_dELL[i] = xwts[ix]*swts[is]*μ[ix]*(sgrid[is]*γ*η[i]/(ν*c[i]))*(1/γ)*unc[i]*ell[i]^(-(1/γ)-1)
            dGDP_dWW    += swts[is]*xwts[ix]*μ[i]*((sgrid[is]*η[i])/ν - (sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipW[i])
            dGDP_dRR    += -swts[is]*xwts[ix]*μ[i]*(sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipR[i]
            dGDP_dTT    += -swts[is]*xwts[ix]*μ[i]*(sgrid[is]*γ*η[i]/(ν*c[i]))*(1-unc[i])*chipX[i]
            dGDP_dZ     += swts[is]*xwts[ix]*sgrid[is]*η[i]*μ[i]
        end
    end
    return dGDP_dMU, dGDP_dZ, dGDP_dELL, dGDP_dRR, dGDP_dWW, dGDP_dTT
end

function update_measurement_covariance_matrices!(m::RealBondMkup, system::System{T}) where {T<:AbstractFloat}
    exo       = m.exogenous_shocks

    # Variance of innovations
    system[:QQ][exo[:z_sh], exo[:z_sh]] = m[:σ_z]^2
    system[:QQ][exo[:mon_sh], exo[:mon_sh]] = m[:σ_mon]^2
    system[:QQ][exo[:mkp_sh], exo[:mkp_sh]] = m[:σ_mkp]^2
end

function construct_GDPeqn(m::RealBondMkup, TTT_jump::Matrix{T}) where {T<:AbstractFloat}
    endo_unnorm = m.endogenous_states_unnormalized

    # Load in parameters, steady-state parameters, and grids
    γ::Float64     = m[:γ].value
    ν::Float64     = m[:ν].value
    R::Float64     = m[:R].value
    abar::Float64  = m[:abar].value
    aborrow        = abar/R

    ell::Vector{Float64}  = m[:lstar].value
    c::Vector{Float64}    = m[:cstar].value
    μ::Vector{Float64}    = m[:μstar].value
    η::Vector{Float64}    = m[:ηstar].value
    χss::Vector{Float64}  = m[:χstar].value

    xwts::Vector{Float64}  = m.grids[:xgrid].weights
    sgrid::Vector{Float64} = m.grids[:sgrid].points
    swts::Vector{Float64}  = m.grids[:sgrid].weights
    xgrid_total::Vector{Float64}  = m.grids[:xgrid_total]

    nx::Int = get_setting(m, :nx)
    ns::Int = get_setting(m, :ns)

    # Construct GDP
    dGDP_dMU, dGDP_dZ, dGDP_dELL, dGDP_dRR, dGDP_dWW, dGDP_dTT =
        construct_GDPfn_realbondmkp(nx, ns, xgrid_total, sgrid, xwts, swts,
                                 γ, ν, abar, R, aborrow, μ, c, η, ell, χss)

    GDPfn = zeros(n_model_states_unnormalized(m))
    # GDP as function of un-normalized MU Z MON ELL RR II WW PI TT
    # note: GDP is only a function of contemporaneous variables
    # we are using the indices (MUP,ZP, etc.) corresponding to date t+1 variables
    # simply because these indices happen to work here also
    # this does not mean that GDP is a function of date t+1 variables
    GDPfn[endo_unnorm[:μ′_t]]  = vec(dGDP_dMU)
    GDPfn[first(endo_unnorm[:z′_t])]  = dGDP_dZ      #first because the index is a range of single number
    GDPfn[endo_unnorm[:l′_t]]  = vec(dGDP_dELL)
    GDPfn[first(endo_unnorm[:R′_t])]  = dGDP_dRR
    GDPfn[first(endo_unnorm[:w′_t])]  = dGDP_dWW
    GDPfn[first(endo_unnorm[:t′_t])]  = dGDP_dTT

    ########################################

    Qx, Qy, _, _ = compose_normalization_matrices(m)

    gx2  = Qy'*TTT_jump*Qx

    # now we need to create GDP as a function of the normalized states
    n_backward_looking_states_unnorm = n_backward_looking_states_unnormalized(m)
    GDPeqn = GDPfn'*[Matrix{Float64}(I, n_backward_looking_states_unnorm, n_backward_looking_states_unnorm); gx2]*Qx'
             # this is for level of GDP
             # to use the log of gdp, front multiply by (1/dGDP_dZ)
    return GDPeqn
end


text:
using Dates
using DataFrames
using JSON
using StatsPlots
pyplot(size = (1600,1200))

dateformat() = "yyyy-mm-dd-T-HH-MM-SS"
getnow() = Dates.format(Dates.now(), dateformat())
parsetime(s) = DateTime(s, dateformat())
parsedataset(s) = occursin("48echo", s) ? "240x240x48x48" : occursin("56echo", s) ? "240x240x113x56" : error("Unknown data source: $s")

function load_results()
    df = DataFrame()
    for dir in readdir(joinpath(@__DIR__, "results"); join = true)
        for res in JSON.parsefile(joinpath(dir, "results.json"))["results"]
            push!(df, (
                date = parsetime(split(basename(dir), '_')[1]),
                name = split(basename(dir), '_')[2] |> String,
                time = res["min"] |> Float64,
                dataset = parsedataset(res["parameters"]["input"]),
                julia = res["parameters"]["julia"],
                optlevel = parse(Int, res["parameters"]["optimize"]),
                threads = parse(Int, res["parameters"]["threads"]),
                version = res["parameters"]["version"],
            ))
        end
    end
    return df
end

function plot_results(df; compare_version = "v0.3")
    df = sort(df, [:dataset, :threads, :julia, :version, :optlevel])
    gd = groupby(df, [:dataset, :threads])
    ptimes = []
    pspeedups = []
    for (k,g) in zip(keys(gd), gd)
        groupkeystr = join(["$k = $v" for (k,v) in pairs(k)], ", ")

        ptime = @df g scatter(
            :version, :time, group = :julia,
            title = groupkeystr, xlabel = "version", ylabel = "time [s]",
            marker = (10, [:circle :square :utriangle :dtriangle]),
            legend = :topright,
        )
        push!(ptimes, ptime)

        if false
            # Times relative to `compare_version` runs with exactly the same parameters
            grel = DataFrame()
            for gsub in groupby(deepcopy(g), [:julia, :threads, :optlevel])
                if compare_version ∈ gsub.version
                    gsub.time .= only(gsub.time[gsub.version .== compare_version]) ./ gsub.time
                    append!(grel, gsub)
                end
            end
        else
            # Times relative to minimum `compare_version` time
            grel = deepcopy(g)
            grel.time .= minimum(grel[grel.version .== compare_version, :time]) ./ grel.time
        end
        pspeedup = @df grel scatter(
            :version, :time, group = :julia,
            title = groupkeystr, xlabel = "version", ylabel = "speedup w.r.t. fastest $compare_version",
            marker = (10, [:circle :square :utriangle :dtriangle]),
            legend = :topleft,
        )
        push!(pspeedups, pspeedup)
    end
    plot(ptimes...) |> display
    plot(pspeedups...) |> display
end

df = load_results()
plot_results(df)


text:
"""
    IterationPrintPs(i,N,msgOld,freq=1,msgExtra="")

Prints iteration progress, overwrites old printing

# Input
- `i::Integer`:       iteration number
- `N::Integer`:       total number of iterations
- `msgOld::String`:   old string, exactly this space is overwritten
- `freq::Integer`:    (optional) indicating when to print
- `msgExtra::String`: (optional) extra string to print

# Output
- `msg::String`:msg, s tring printed to screen. Use as msgOld in the next iteration

# Uses
- Printf

Paul.Soderlind@unisg.ch

"""
function IterationPrintPs(i,N,msgOld,freq=1,msgExtra="")
  if mod(i,freq) == 0                      #print only if i/freq is an integer
    reverseStr = "\b" ^ length(msgOld)
    msg        = string(@sprintf("Processed %d/%d",i,N),msgExtra)
    msg        = rpad(msg,60," ")                         #fill at least 60 columns
    @printf("%s",reverseStr * msg)
  else
    msg = msgOld
  end
  return msg
end
#------------------------------------------------------------------------------

max_stars_repo_name: PaulSoderlind/PaulSoderlindCode

max_issues_repo_name: PaulSoderlind/PaulSoderlindCode

max_forks_repo_name: PaulSoderlind/PaulSoderlindCode