text:
"""
update_θ_jensen!(gcm)

Update θ using the MM algorithm and Jensens inequality, given β.
"""
function update_θ_jensen!(
    gcm::Union{GLMCopulaVCModel{T, D, Link}, NBCopulaVCModel{T, D, Link}, GLMCopulaARModel{T, D, Link}, GLMCopulaCSModel{T, D, Link}, NBCopulaARModel{T, D, Link}, NBCopulaCSModel{T, D, Link}, Poisson_Bernoulli_VCModel{T, VD, VL}},
    maxiter::Integer=50000,
    reltol::Number=1e-6,
    verbose::Bool=false) where {T <: BlasReal, D<:Union{Poisson, Bernoulli, NegativeBinomial}, Link, VD, VL}
    rsstotal = zero(T)
    @inbounds for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
        rsstotal += abs2(norm(gcm.data[i].res))  # needed for updating τ in normal case
        standardize_res!(gcm.data[i])            # standardize the residuals GLM variance(μ)
        QuasiCopula.update_quadform!(gcm.data[i]) # with standardized residuals
        gcm.QF[i, :] = gcm.data[i].q
    end
    # MM iteration
    for iter in 1:maxiter
        # store previous iterate
        copyto!(gcm.storage_θ, gcm.θ)

        # update τ if necessary
        mul!(gcm.storage_n, gcm.QF, gcm.θ) # gcm.storage_n[i] = sum_k^m qi[k] sigmai_[k] # denom of numerator
        fill!(gcm.τ, 1.0)

        ##### Numerator in the multiplicative update ##########
        # when its logistic gcm.τ = 1.0
        gcm.storage_n .= inv.(inv.(gcm.τ) .+ gcm.storage_n) # 1/ (1 + sum_k^m qi[k] sigmai_[k]) # denom of numerator
        mul!(gcm.storage_m, transpose(gcm.QF), gcm.storage_n) # store numerator of numerator
        gcm.θ .*= gcm.storage_m # multiply by the numerator

        ######## Denominator in the multiplicative update  ###########
        mul!(gcm.storage_n, gcm.TR, gcm.storage_θ)
        gcm.storage_n .= inv.(1 .+ gcm.storage_n) # denominator of the denominator

        mul!(gcm.storage_m, transpose(gcm.TR), gcm.storage_n) # numerator of the denominator
        gcm.θ ./= gcm.storage_m # divide by the denominator

        # monotonicity diagnosis (for the normal distribution there is an extra term)
        verbose && println(sum(log, 1 .+ (gcm.QF * gcm.θ)) -
            sum(log, 1 .+ gcm.TR * gcm.θ))
        # 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

"""
    update_res!(gc, β)
Update the residual vector according to `β` given link function and distribution.
"""
function update_res!(
   gc::Union{GLMCopulaVCObs, NBCopulaVCObs, GLMCopulaCSObs, GLMCopulaARObs, NBCopulaARObs, NBCopulaCSObs},
   β::Vector)
   mul!(gc.η, gc.X, β)
   @turbo for i in 1:gc.n
       gc.μ[i] = GLM.linkinv(gc.link, gc.η[i])
       gc.varμ[i] = GLM.glmvar(gc.d, gc.μ[i]) # Note: for negative binomial, d.r is used
       gc.dμ[i] = GLM.mueta(gc.link, gc.η[i])
       gc.w1[i] = gc.dμ[i] / gc.varμ[i]
       gc.w2[i] = gc.w1[i] * gc.dμ[i]
       gc.res[i] = gc.y[i] - gc.μ[i]
   end
   return gc.res
end

"""
    update_res!(gc, β)
Update the residual vector according to `β` given link function and distribution.
"""
function update_res!(
   gc::Poisson_Bernoulli_VCObs,
   β::Vector)
   mul!(gc.η, gc.X, β)
   @turbo for i in 1:gc.n
       gc.μ[i] = GLM.linkinv(gc.veclink[i], gc.η[i])
       gc.varμ[i] = GLM.glmvar(gc.vecd[i], gc.μ[i]) # Note: for negative binomial, d.r is used
       gc.dμ[i] = GLM.mueta(gc.veclink[i], gc.η[i])
       gc.w1[i] = gc.dμ[i] / gc.varμ[i]
       gc.w2[i] = gc.w1[i] * gc.dμ[i]
       gc.res[i] = gc.y[i] - gc.μ[i]
   end
   return gc.res
end

function update_res!(
    gcm::Union{GLMCopulaVCModel, GLMCopulaARModel, GLMCopulaCSModel, NBCopulaVCModel, NBCopulaARModel, NBCopulaCSModel, Poisson_Bernoulli_VCModel}
    )
    @inbounds for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
    end
    nothing
end

function standardize_res!(
    gc::Union{GLMCopulaVCObs, NBCopulaVCObs, GLMCopulaARObs, GLMCopulaCSObs, NBCopulaARObs, NBCopulaCSObs, Poisson_Bernoulli_VCObs}
    )
    @turbo for j in eachindex(gc.y)
        gc.res[j] /= sqrt(gc.varμ[j])
    end
end

function standardize_res!(
    gcm::Union{GLMCopulaVCModel, NBCopulaVCModel, GLMCopulaARModel, GLMCopulaCSModel}
    )
    # standardize residual
    @inbounds for i in eachindex(gcm.data)
        standardize_res!(gcm.data[i])
    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::Union{GLMCopulaVCObs{T, D, Link}, NBCopulaVCObs{T, D, Link}, Poisson_Bernoulli_VCObs{T, VD, VL}}) where {T<:Real, D, Link, VD, VL}
    @inbounds for k in 1:length(gc.V)
        mul!(gc.storage_n, gc.V[k], gc.res)
        gc.q[k] = dot(gc.res, gc.storage_n) / 2
    end
    gc.q
end

"""
    update_quadform!(gc)

Update the quadratic forms `(r^T V[k] r) / 2` according to the current residual `r`.
"""
function update_quadform!(gc::Union{GLMCopulaARObs{T, D, Link}, GLMCopulaCSObs{T, D, Link}, NBCopulaARObs{T, D, Link}, NBCopulaCSObs{T, D, Link}}) where {T<:Real, D, Link}
    mul!(gc.storage_n, gc.V, gc.res)
    gc.q .= dot(gc.res, gc.storage_n) / 2
    gc.q
end

"""
update_res!(gc, β)
Update the residual vector according to `β`.
"""
function update_res!(
    gc::GaussianCopulaVCObs{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::GaussianCopulaVCModel{T}
    ) where T <: BlasReal
    @inbounds for i in eachindex(gcm.data)
        update_res!(gcm.data[i], gcm.β)
    end
    nothing
end

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

function standardize_res!(
    gcm::GaussianCopulaVCModel{T}
    ) where T <: BlasReal
    σinv = sqrt(gcm.τ[1])
    # standardize residual
    @inbounds 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)
    @inbounds 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)
    @inbounds 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)
    @inbounds 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!

"""
    update_r!(gc)

Performs maximum loglikelihood estimation of the nuisance paramter for negative
binomial model using Newton's algorithm. Will run a maximum of `maxIter` and
convergence is defaulted to `convTol`.
"""
function update_r_newton!(gcm::Union{NBCopulaVCModel, NBCopulaARModel, NBCopulaCSModel};
    maxIter=100, convTol=1e-6)

    T = eltype(gcm.β)
    r = gcm.r[1] # estimated r in previous iteration
    new_r = one(T)
    stepsize = one(T)
    for i in 1:maxIter
        old_logl = negbin_component_loglikelihood(gcm, r)

        # run 1 iteration of Newton's algorithm
        increment = newton_increment(gcm, r)
        new_r = r - stepsize * increment

        if new_r <= 0
            new_r = 1.0
            @warn "New r estimated to be negative, restarting at r = 1..."
        end

        # linesearch
        # for j in 1:20
        #     if new_r <= 0
        #         stepsize = stepsize / 2
        #         new_r = r - stepsize * increment
        #     else
        #         new_logl = negbin_component_loglikelihood(gcm, new_r)
        #         if old_logl >= new_logl
        #             stepsize = stepsize / 2
        #             new_r = r - stepsize * increment
        #         else
        #             break
        #         end
        #     end
        # end

        #check convergence
        if abs(r - new_r) ≤ convTol
            break
        else
            r = new_r
        end
    end
    return new_r
end

function update_r!(gcm::Union{NBCopulaVCModel, NBCopulaARModel, NBCopulaCSModel})
    new_r = update_r_newton!(gcm, maxIter=10)
    gcm.r[1] = new_r
    for gc in gcm.data
        gc.d = NegativeBinomial(new_r)
    end
    return nothing
end

function first_derivative(gcm::NBCopulaVCModel, r::T) where T <: AbstractFloat
    s = zero(T)
    @inbounds for i in eachindex(gcm.data)
        # 2nd term of logl
        y = gcm.data[i].y
        μ = gcm.data[i].μ
        for j in eachindex(y)
            s += dLdr_2ndterm(y[j], μ[j], r)
        end
        # 3rd term of logl
        resid = gcm.data[i].res
        Γ = gcm.data[i].storage_nn # Γ = a1*V1 + ... + am*Vm
        η = gcm.data[i].η
        # Γ = gcm.θ' * gcm.data[i].V
        # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
        # dD = Diagonal([-exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) for j in 1:length(η)])
        # dresid = -inv(D)*dD*resid
        # s += resid'*Γ*dresid / (1 + 0.5resid'*Γ*resid)
        for j in 1:length(η)
            gcm.data[i].storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n[i] = 1 / Di = 1 / sqrt(var(yi))
        end
        for j in 1:length(η)
            gcm.data[i].storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) # storage_n = inv(D) * dD
        end
        gcm.data[i].storage_n .*= -1.0 .* resid # storage_n = dr(β) (derivative of residuals)
        mul!(gcm.data[i].storage_n2, Γ, gcm.data[i].storage_n) # storage_n2 = Γ * dresid
        numer = dot(resid, gcm.data[i].storage_n2) # numer = r' * Γ * dr
        mul!(gcm.data[i].storage_n2, Γ, resid) # storage_n = Γ * resid
        denom = 1 + 0.5 * dot(resid, gcm.data[i].storage_n2) # denom = 1 + 0.5(r * Γ * r)
        s += numer / denom
    end
    return s
end

function first_derivative(gcm::Union{NBCopulaARModel, NBCopulaCSModel}, r::T) where T <: AbstractFloat
    s = zero(T)
    @inbounds for i in eachindex(gcm.data)
        # 2nd term of logl
        y = gcm.data[i].y
        μ = gcm.data[i].μ
        for j in eachindex(y)
            s += dLdr_2ndterm(y[j], μ[j], r)
        end
        # 3rd term of logl
        get_V!(gcm.ρ[1], gcm.data[i])
        Γ = gcm.σ2[1] * gcm.data[i].V
        η = gcm.data[i].η
        resid = gcm.data[i].res
        # Γ = gcm.Σ' * gcm.data[i].V
        # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
        # dD = Diagonal([-exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) for j in 1:length(η)])
        # dresid = -inv(D)*dD*resid
        # s += resid'*Γ*dresid / (1 + 0.5resid'*Γ*resid)
        for j in 1:length(η)
            gcm.data[i].storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n[i] = 1 / Di = 1 / sqrt(var(yi))
        end
        for j in 1:length(η)
            gcm.data[i].storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) # storage_n = inv(D) * dD
        end
        gcm.data[i].storage_n .*= -1.0 .* resid # storage_n = dr(β) (derivative of residuals)
        mul!(gcm.data[i].storage_n2, Γ, gcm.data[i].storage_n) # storage_n2 = Γ * dresid
        numer = dot(resid, gcm.data[i].storage_n2) # numer = r' * Γ * dr
        mul!(gcm.data[i].storage_n2, Γ, resid) # storage_n = Γ * resid
        denom = 1 + 0.5 * dot(resid, gcm.data[i].storage_n2) # denom = 1 + 0.5(r * Γ * r)
        s += numer / denom
    end
    return s
end

function second_derivative(gcm::NBCopulaVCModel, r::T) where T <: AbstractFloat
    s = zero(T)
    @inbounds for i in eachindex(gcm.data)
        # 2nd term of logl
        y = gcm.data[i].y
        μ = gcm.data[i].μ
        for j in eachindex(y)
            s += dLdr2_2ndterm(y[j], μ[j], r)
        end
        # 3rd term of logl
        Γ = gcm.data[i].storage_nn # Γ = a1*V1 + ... + am*Vm
        η = gcm.data[i].η
        resid = gcm.data[i].res
        # Γ = gcm.θ' * gcm.data[i].V
        # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
        # dD = Diagonal([-exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) for j in 1:length(η)])
        # d2D = Diagonal([(exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
        #     (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))) for j in 1:length(η)])
        # resid = gcm.data[i].res
        # dresid = -inv(D)*dD*resid
        # d2resid = (2inv(D)*dD*inv(D)*dD - inv(D)*d2D)*resid
        # denom = 1 + 0.5resid'*Γ*resid
        # term1 = (resid'*Γ*dresid / denom)^2
        # term2 = dresid'*Γ*dresid / denom
        # term3 = resid'*Γ*d2resid / denom
        # s += -term1 + term2 + term3
        for j in 1:length(η)
            gcm.data[i].storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n = inv(Di) = 1 / sqrt(var(yi))
        end
        for j in 1:length(η)
            # storage_n2 = -inv(D) * d2D
            gcm.data[i].storage_n2[j] = -gcm.data[i].storage_n[j] *
                ((exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
                (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))))
        end
        for j in 1:length(η)
            # storage_n = inv(D) * dD
            gcm.data[i].storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r)))
        end
        for j in 1:length(η)
            # storage_n2 = 2inv(D)*dD*inv(D)*dD -inv(D)*d2D
            gcm.data[i].storage_n2[j] += 2 * abs2(gcm.data[i].storage_n[j])
        end
        gcm.data[i].storage_n .*= -1.0 .* resid # storage_n = dr(β) = derivative of residuals
        gcm.data[i].storage_n2 .*= resid # storage_n2 = dr²(β) = 2nd derivative of residuals
        mul!(gcm.data[i].storage_n3, Γ, resid) # storage_n3 = Γ * resid
        denom = 1 + 0.5 * dot(resid, gcm.data[i].storage_n3)
        mul!(gcm.data[i].storage_n3, Γ, gcm.data[i].storage_n) # storage_n3 = Γ * dresid
        term1 = (dot(resid, gcm.data[i].storage_n3) / denom)^2 # (resid'*Γ*dresid / denom)^2
        term2 = dot(gcm.data[i].storage_n, gcm.data[i].storage_n3) / denom # term2 = dresid'*Γ*dresid / denom
        mul!(gcm.data[i].storage_n3, Γ, gcm.data[i].storage_n2) # storage_n3 = Γ * d2resid
        term3 = dot(resid, gcm.data[i].storage_n3) / denom # term3 = resid'*Γ*d2resid / denom
        s += -term1 + term2 + term3
    end
    return s
end

function second_derivative(gcm::Union{NBCopulaARModel, NBCopulaCSModel}, r::T) where T <: AbstractFloat
    s = zero(T)
    @inbounds for i in eachindex(gcm.data)
        # 2nd term of logl
        y = gcm.data[i].y
        μ = gcm.data[i].μ
        for j in eachindex(y)
            s += dLdr2_2ndterm(y[j], μ[j], r)
        end
        # 3rd term of logl
        get_V!(gcm.ρ[1], gcm.data[i])
        Γ = gcm.σ2[1] * gcm.data[i].V
        η = gcm.data[i].η
        resid = gcm.data[i].res
        # Γ = gcm.Σ' * gcm.data[i].V
        # D = Diagonal([sqrt(exp(η[j])*(exp(η[j])+r) / r) for j in 1:length(η)])
        # dD = Diagonal([-exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r))) for j in 1:length(η)])
        # d2D = Diagonal([(exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
        #     (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))) for j in 1:length(η)])
        # resid = gcm.data[i].res
        # dresid = -inv(D)*dD*resid
        # d2resid = (2inv(D)*dD*inv(D)*dD - inv(D)*d2D)*resid
        # denom = 1 + 0.5resid'*Γ*resid
        # term1 = (resid'*Γ*dresid / denom)^2
        # term2 = dresid'*Γ*dresid / denom
        # term3 = resid'*Γ*d2resid / denom
        # s += -term1 + term2 + term3
        for j in 1:length(η)
            gcm.data[i].storage_n[j] = inv(sqrt(exp(η[j])*(exp(η[j])+r) / r)) # storage_n = inv(Di) = 1 / sqrt(var(yi))
        end
        for j in 1:length(η)
            # storage_n2 = -inv(D) * d2D
            gcm.data[i].storage_n2[j] = -gcm.data[i].storage_n[j] *
                ((exp(3η[j]) / (4r^1.5 * (exp(η[j])*(exp(η[j])+r))^(1.5))) +
                (3exp(2η[j]) / (4r^(2.5)*sqrt(exp(η[j])*(exp(η[j])+r)))))
        end
        for j in 1:length(η)
            # storage_n = inv(D) * dD
            gcm.data[i].storage_n[j] *= -exp(2η[j]) / (2r^1.5 * sqrt(exp(η[j])*(exp(η[j])+r)))
        end
        for j in 1:length(η)
            # storage_n2 = 2inv(D)*dD*inv(D)*dD -inv(D)*d2D
            gcm.data[i].storage_n2[j] += 2 * abs2(gcm.data[i].storage_n[j])
        end
        gcm.data[i].storage_n .*= -1.0 .* resid # storage_n = dr(β) = derivative of residuals
        gcm.data[i].storage_n2 .*= resid # storage_n2 = dr²(β) = 2nd derivative of residuals
        mul!(gcm.data[i].storage_n3, Γ, resid) # storage_n3 = Γ * resid
        denom = 1 + 0.5 * dot(resid, gcm.data[i].storage_n3)
        mul!(gcm.data[i].storage_n3, Γ, gcm.data[i].storage_n) # storage_n3 = Γ * dresid
        term1 = (dot(resid, gcm.data[i].storage_n3) / denom)^2 # (resid'*Γ*dresid / denom)^2
        term2 = dot(gcm.data[i].storage_n, gcm.data[i].storage_n3) / denom # term2 = dresid'*Γ*dresid / denom
        mul!(gcm.data[i].storage_n3, Γ, gcm.data[i].storage_n2) # storage_n3 = Γ * d2resid
        term3 = dot(resid, gcm.data[i].storage_n3) / denom # term3 = resid'*Γ*d2resid / denom
        s += -term1 + term2 + term3
    end
    return s
end

function negbin_component_loglikelihood(gcm::Union{NBCopulaARModel, NBCopulaCSModel}, r::T) where T <: AbstractFloat
    return loglikelihood!(gcm)
end

function negbin_component_loglikelihood(gcm::NBCopulaVCModel, r::T) where T <: AbstractFloat
    logl = zero(T)
    for gc in gcm.data
        n, p, m = size(gc.X, 1), size(gc.X, 2), length(gc.V)
        # fill!(gc.∇β, 0.0)
        update_res!(gc, gcm.β)
        standardize_res!(gc)
        # fill!(gc.∇resβ, 0.0) # fill gradient of residual vector with 0
        std_res_differential!(gc) # this will compute ∇resβ
        fill!(gc.storage_nn, 0)
        @inbounds for k in 1:m
            mul!(gc.storage_n, gc.V[k], gc.res) # storage_n = V[k] * res
            gc.q[k] = dot(gc.res, gc.storage_n) / 2 # gc.q[k] = 0.5res' * V[k] * res
            # BLAS.gemv!('T', gcm.θ[k], gc.∇resβ, gc.storage_n, 1.0, gc.∇β) # gc.∇β += ∇resβ*Γ*res (standardized residual)
            gc.storage_nn .+= gcm.θ[k] .* gc.V[k] # compute Γ = a1*V1 + ... am*Vm
        end
        # 2nd term of logl
        logl += component_loglikelihood(gc, r)
        # 3rd term of logl
        qsum  = dot(gcm.θ, gc.q) # q[k] = res_i' * V_i[k] * res_i / 2, so qsum = 0.5r(β)*Γ*r(β)
        logl += log(1 + qsum)
    end
    return logl
end

function newton_increment(gcm, r)
    dx = first_derivative(gcm, r)
    dx2 = second_derivative(gcm, r)
    increment = dx / dx2
    # use gradient ascent if hessian not negative definite
    if dx2 < 0
        increment = dx / dx2
    else
        increment = dx
    end
    return increment
end

function dLdr_2ndterm(yi::T, μi::T, r::T) where T <: AbstractFloat
    return -(yi+r)/(μi+r) - log(μi+r) + 1 + log(r) + digamma(r+yi) - digamma(r) :: T
end

function dLdr2_2ndterm(yi::T, μi::T, r::T) where T <: AbstractFloat
    return (yi+r)/(μi+r)^2 - 2/(μi+r) + 1/r + trigamma(r+yi) - trigamma(r) :: T
end


text:
module Juno

using Media, Base64, Profile

import Media: render

export Media, media, @render

_active = false

isprecompiling() = ccall(:jl_generating_output, Cint, ()) == 1

activate() = return nothing #!isprecompiling() && @eval using Atom

setactive!(active) = (global _active = active)

"""
    isactive()

Will return `true` when the current Julia process is connected to a running Juno
frontend.
"""
isactive() = _active

include("types.jl")
include("frontend.jl")
include("progress.jl")
include("user.jl")
include("utils.jl")

# precompile
# include("../deps/SnoopCompile/precompile/precompile_Juno.jl")
# _precompile_()

end # module

ext: jl


text:
#
# User-facing Model types providing a simplified API to model definitions and instances.
#

abstract type AbstractModel <: MimiStruct end

"""
    Model

A user-facing API containing a `ModelInstance` (`mi`) and a `ModelDef` (`md`).
This `Model` can be created with the optional keyword argument `number_type` indicating
the default type of number used for the `ModelDef`.  If not specified the `Model` assumes
a `number_type` of `Float64`.
"""
mutable struct Model <: AbstractModel
    md::ModelDef
    mi::Union{Nothing, ModelInstance}

    function Model(number_type::DataType=Float64)
        return new(ModelDef(number_type), nothing)
    end

    # Create a copy of a model, e.g., to create marginal models
    function Model(m::Model)
        return new(deepcopy(m.md), nothing)
    end

    # Create a model from a ModelInstance (temporary for explore call)
    function Model(mi::ModelInstance)
        return new(deepcopy(mi.md), deepcopy(mi))
    end
end

"""
    MarginalModel

A Mimi `Model` whose results are obtained by subtracting results of one `base` Model
from those of another `marginal` Model that has a difference of `delta`.
"""
struct MarginalModel <: AbstractModel
    base::Model
    modified::Model
    delta::Float64

    function MarginalModel(base::Model, delta::Float64=1.0)
        return new(base, Model(base), delta)
    end
end

function Base.getindex(mm::MarginalModel, comp_name::Symbol, name::Symbol)
    return (mm.modified[comp_name, name] .- mm.base[comp_name, name]) ./ mm.delta
end

function Base.getindex(mm::MarginalModel, comp_path::ComponentPath, name::Symbol)
    return (mm.modified.mi[comp_path, name] .- mm.base.mi[comp_path, name]) ./ mm.delta
end

# DEPRECATION - EVENTUALLY REMOVE (and go back to default getproperty behavior)
function Base.getproperty(base::MarginalModel, s::Symbol)
    if (s == :marginal)
        error("Use of `MarginalModel.marginal` is deprecated in favor of `MarginalModel.modified`.")
    end
    return getfield(base, s);
end


text:
# =============================================================
# Experimenting with Stefan Karpinski's "Shelling out sucks"
# https://julialang.org/blog/2012/03/shelling-out-sucks
# =============================================================

# NOTE (quirky things that may be ok one way or the other)
# - skipping a line in a list --> in jekyll stays part of the same block; not in JuDoc (but that's the julia markdown parser) ==> would require processing lists yourself
#
# =============================================================

# t1 = process shortcuts for links

st = """
    [Perl]:     http://www.perl.org/
    [Python]:   http://python.org/

    (...) like [Perl] and [Ruby].
    """

# t2 = quadruple spaces at the start of a non-empty line indicate code (like ```) provided the
# line is preceded by an empty line

st = """
    Blah

        this is code

    end
    """

st = """
    This will fail:

        irb(main):001:0> dir="src"
        => "src"
        irb(main):002:0> `find #{dir} -type f -print0 | xargs -0 grep foo | wc -l`.to_i
        => 5

    The simple
    """


text:
const _breakpoints = BreakpointRef[]
breakpoints() = copy(_breakpoints)

Base.getindex(bp::BreakpointRef) = bp.framecode.breakpoints[bp.stmtidx]
function Base.setindex!(bp::BreakpointRef, isactive::Bool)
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(isactive, bp[].condition)
end
function toggle!(bp::BreakpointRef)
    state = bp[]
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(!state.isactive, state.condition)
end

function add_breakpoint(framecode, stmtidx)
    bp = BreakpointRef(framecode, stmtidx)
    # Since there can be only one BreakpointState for a given framecode/stmtidx,
    # check whether _breakpoints is already storing a reference to that location
    idx = findfirst(isequal(bp), _breakpoints)
    if idx === nothing
        push!(_breakpoints, bp)
    end
    return bp
end

function shouldbreak(frame::Frame, pc::Int)
    bps = frame.framecode.breakpoints
    isassigned(bps, pc) || return false
    bp = bps[pc]
    bp.isactive || return false
    return Base.invokelatest(bp.condition, frame)::Bool
end

function prepare_slotfunction(framecode::FrameCode, body::Union{Symbol,Expr})
    ismeth = framecode.scope isa Method
    uslotnames = Set{Symbol}()
    slotnames  = Symbol[]
    for name in framecode.src.slotnames
        if name ∉ uslotnames
            push!(slotnames, name)
            push!(uslotnames, name)
        end
    end
    framename, dataname = gensym("frame"), gensym("data")
    assignments = Expr[:($dataname = $framename.framedata)]
    default = Unassigned()
    for i = 1:length(slotnames)
        slotname = framecode.src.slotnames[i]
        qslotname = QuoteNode(slotname)
        getexpr = :(something($dataname.locals[$dataname.last_reference[$qslotname]]))
        push!(assignments, Expr(:(=), slotname, :(haskey($dataname.last_reference, $qslotname) ? $getexpr : $default)))
    end
    if ismeth
        syms = sparam_syms(framecode.scope)
        for i = 1:length(syms)
            push!(assignments, Expr(:(=), syms[i], :($dataname.sparams[$i])))
        end
    end
    funcname = ismeth ? gensym("slotfunction") : gensym(Symbol(framecode.scope.name, "_slotfunction"))
    return Expr(:function, Expr(:call, funcname, framename), Expr(:block, assignments..., body))
end

const Condition = Union{Nothing,Expr,Tuple{Module,Expr}}
_unpack(condition) = isa(condition, Expr) ? (Main, condition) : condition

## The fundamental implementations of breakpoint-setting
function breakpoint!(framecode::FrameCode, pc, condition::Condition=nothing)
    stmtidx = pc
    if condition === nothing
        framecode.breakpoints[stmtidx] = BreakpointState()
    else
        mod, cond = _unpack(condition)
        fex = prepare_slotfunction(framecode, cond)
        framecode.breakpoints[stmtidx] = BreakpointState(true, Core.eval(mod, fex))
    end
    return add_breakpoint(framecode, stmtidx)
end
breakpoint!(frame::Frame, pc=frame.pc, condition::Condition=nothing) =
    breakpoint!(frame.framecode, pc, condition)

"""
    enable(bp::BreakpointRef)

Enable breakpoint `bp`.
"""
enable(bp::BreakpointRef)  = bp[] = true

"""
    disable(bp::BreakpointRef)

Disable breakpoint `bp`. Disabled breakpoints can be re-enabled with [`enable`](@ref).
"""
disable(bp::BreakpointRef) = bp[] = false

"""
    remove(bp::BreakpointRef)

Remove (delete) breakpoint `bp`. Removed breakpoints cannot be re-enabled.
"""
function remove(bp::BreakpointRef)
    idx = findfirst(isequal(bp), _breakpoints)
    deleteat!(_breakpoints, idx)
    bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(false, falsecondition)
    return nothing
end

"""
    enable()

Enable all breakpoints.
"""
enable() = for bp in _breakpoints enable(bp) end

"""
    disable()

Disable all breakpoints.
"""
disable() = for bp in _breakpoints disable(bp) end

"""
    remove()

Remove all breakpoints.
"""
function remove()
    for bp in _breakpoints
        bp.framecode.breakpoints[bp.stmtidx] = BreakpointState(false, falsecondition)
    end
    empty!(_breakpoints)
    return nothing
end

"""
    break_on(states...)

Turn on automatic breakpoints when any of the conditions described in `states` occurs.
The supported states are:

- `:error`: trigger a breakpoint any time an uncaught exception is thrown
"""
function break_on(states::Vararg{Symbol})
    for state in states
        if state == :error
            break_on_error[] = true
        else
            throw(ArgumentError(string("unsupported state :", state)))
        end
    end
end

"""
    break_off(states...)

Turn off automatic breakpoints when any of the conditions described in `states` occurs.
See [`break_on`](@ref) for a description of valid states.
"""
function break_off(states::Vararg{Symbol})
    for state in states
        if state == :error
            break_on_error[] = false
        else
            throw(ArgumentError(string("unsupported state :", state)))
        end
    end
end

"""
    breakpoint(f, sig)
    breakpoint(f, sig, line)
    breakpoint(f, sig, condition)
    breakpoint(f, sig, line, condition)
    breakpoint(...; enter_generated=false)

Add a breakpoint to `f` with the specified argument types `sig`.
Optionally specify an absolute line number `line` in the source file; the default
is to break upon entry at the first line of the body.
Without `condition`, the breakpoint will be triggered every time it is encountered;
the second only if `condition` evaluates to `true`.
`condition` should be written in terms of the arguments and local variables of `f`.

# Example
```julia
function radius2(x, y)
    return x^2 + y^2
end

breakpoint(radius2, Tuple{Int,Int}, :(y > x))
```
"""
function breakpoint(f, sig::Type, line::Integer, condition::Condition=nothing; enter_generated=false)
    method = which(f, sig)
    framecode, _ = prepare_framecode(method, sig; enter_generated=enter_generated)
    # Don't use statementnumber(method, line) in case it's enter_generated
    linec = line - whereis(method)[2] + method.line
    stmtidx = statementnumber(framecode, linec)
    breakpoint!(framecode, stmtidx, condition)
end
function breakpoint(f, sig::Type, condition::Condition=nothing; enter_generated=false)
    method = which(f, sig)
    framecode, _ = prepare_framecode(method, sig; enter_generated=enter_generated)
    breakpoint!(framecode, 1, condition)
end

"""
    breakpoint(method::Method)
    breakpoint(method::Method, line)
    breakpoint(method::Method, condition::Expr)
    breakpoint(method::Method, line, condition::Expr)

Add a breakpoint to `method`.
"""
function breakpoint(method::Method, line::Integer, condition::Condition=nothing)
    framecode, stmtidx = statementnumber(method, line)
    breakpoint!(framecode, stmtidx, condition)
end
function breakpoint(method::Method, condition::Condition=nothing)
    framecode = get_framecode(method)
    breakpoint!(framecode, 1, condition)
end

"""
    breakpoint(f)
    breakpoint(f, condition)

Break-on-entry to all methods of `f`.
"""
function breakpoint(f, condition::Condition=nothing)
    bps = BreakpointRef[]
    for method in methods(f)
        push!(bps, breakpoint(method, condition))
    end
    return bps
end

"""
    breakpoint(filename, line)

Set a breakpoint at the specified file and line number.
"""
function breakpoint(filename::AbstractString, line::Integer, args...)
    local sigs
    try
        sigs = signatures_at(filename, line)
    catch
        sigs = nothing
    end
    if sigs === nothing
        # TODO: build a Revise-free fallback. Note this won't work well for methods with keywords.
        error("no signatures found at $filename, $line.\nRestarting and `using Revise` and the relevant package may fix this problem.")
    end
    for sig in sigs
        method = JuliaInterpreter.whichtt(sig)
        method === nothing && continue
        # Check to see if this method really contains that line. Methods that fill in a default positional argument,
        # keyword arguments, and @generated sections may not contain the line.
        _, line1 = whereis(method)
        offset = line1 - method.line
        src = JuliaInterpreter.get_source(method)
        lastline = src.linetable[end]
        if getline(lastline) + offset >= line
            return breakpoint(method, line, args...)
        end
    end
    error("no signatures found at $filename, $line among the signatures $sigs")
end

"""
    @breakpoint f(args...) condition=nothing
    @breakpoint f(args...) line condition=nothing

Break upon entry, or at the specified line number, in the method called by `f(args...)`.
Optionally supply a condition expressed in terms of the arguments and internal variables
of the method.
If `line` is supplied, it must be a literal integer.

# Example

Suppose a method `mysum` is defined as follows, where the numbers to the left are the line
number in the file:

```
12 function mysum(A)
13     s = zero(eltype(A))
14     for a in A
15         s += a
16     end
17     return s
18 end
```

Then

```
@breakpoint mysum(A) 15 s>10
```

would cause execution of the loop to break whenever `s>10`.
"""
macro breakpoint(call_expr, args...)
    whichexpr = InteractiveUtils.gen_call_with_extracted_types(__module__, :which, call_expr)
    haveline, line, condition = false, 0, nothing
    while !isempty(args)
        arg = first(args)
        if isa(arg, Integer)
            haveline, line = true, arg
        else
            condition = arg
        end
        args = Base.tail(args)
    end
    condexpr = condition === nothing ? nothing : Expr(:quote, condition)
    if haveline
        return quote
            local method = $whichexpr
            $breakpoint(method, $line, $condexpr)
        end
    else
        return quote
            local method = $whichexpr
            $breakpoint(method, $condexpr)
        end
    end
end

const __BREAKPOINT_MARKER__ = nothing

"""
    @bp

Insert a breakpoint at a location in the source code.
"""
macro bp()
    return esc(:($(JuliaInterpreter).__BREAKPOINT_MARKER__))
end


text:
function OilInflow(pr,J;pb=0,n=10)

    pres=zeros(n,size(pr,1))
    Qrate=zeros(n,size(pr,1)*size(J,1))
    Sat=pr.<=pb

    leg=Array{String,1}(undef,size(pr,1)*size(J,1))
    count=0
    for p=1:size(pr,1)

        pres[:,p]=range(0,stop=pr[p]-50,length=n)

        for j=1:size(J,1)
             count=count+1
             leg[count]="Res. Pressure=$(pr[p]) psi; J=$(J[j]) bbl/d*psi"
            if Sat[p]==1
               Qrate[:,count]=map(pwf->(pr[p]*J[j]/1.8).*(1-0.2*(pwf./pr[p])-0.8*(pwf./pr[p]).^2),pres[:,p])
            else
                Qrate[pres[:,p].>=pb,count]=map(pwf->J[j]*(pr[p]-pwf),pres[:,p][pres[:,p].>=pb])
                Qrate[pres[:,p].<pb,count]=map(pwf->(J[j]*(pr[p]-pb))+(J[j]*pb/1.8)*(1-0.2*(pwf./pb)-0.8*(pwf./pb).^2),pres[:,p][pres[:,p].<pb])

            end

        end

    end
    return pres, Qrate, leg
end


text:
using TrajectoryOptimization
using Test
using BenchmarkTools
using LinearAlgebra
using Random
using SparseArrays
using ForwardDiff
using Logging

disable_logging(Logging.Info)

@testset "Simple Pendulum" begin
    include("simple_pendulum.jl")
end
@testset "Dubins Car" begin
    include("dubins_car_tests.jl")
end
@testset "Quadrotor" begin
    include("quadrotor_tests.jl")
end
@testset "Constrained Objective" begin
    include("objective_tests.jl")
end
@testset "Results" begin
    include("results_tests.jl")
end
@testset "Jacobians" begin
    include("jacobian_tests.jl")
end
@testset "Square Root Method" begin
    include("sqrt_method_tests.jl")
end
@testset "Infeasible Start" begin
    include("infeasible_start_tests.jl")
end
@testset "Minimum Time" begin
    include("minimum_time_tests.jl")
end
# @testset "Direct Collocation" begin
#     include("dircol_test.jl")
#     include("ipopt_test.jl")
# end

disable_logging(Logging.Debug)

"""
# NEEDED TESTS:
- constraint generator

- Attempt undefined integration scheme
- All dynamics
- Custom terminal constraints in objective
- Try/catch in solve_sqrt
- UnconstrainedResults constructor
- more advanced infeasible start
"""


text:
using Iris
using SpecialFunctions

c = Circle(1,2,3)
lat = Lattices.Cartesian(1/21,1/21)
bnd = Boundary(c,DirichletBC())
sim = Simulation(Cavity(bnd,lat))
k,ψ = eig_kl(sim,eps(),nev=20)
@testset "Zeros of Bessel J" begin
    @test abs(besselj(0,k[1])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(1,k[2])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(1,k[3])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(2,k[4])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(2,k[5])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(0,k[6])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(3,k[7])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(3,k[8])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(1,k[9])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(1,k[10])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(4,k[11])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(4,k[12])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(2,k[13])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(2,k[14])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(0,k[15])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(5,k[16])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(5,k[17])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(3,k[18])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(3,k[19])) ≤ max(lat.dx,lat.dy)
    @test abs(besselj(6,k[20])) ≤ max(lat.dx,lat.dy)
end

r = Rectangle(1,2,3,4,.1)
bnd = Boundary(r,NeumannBC())
sim = Simulation(Cavity(bnd,lat))
k,ψ = eig_kl(sim,eps(),nev=5)
@testset "1:2 Rectangle" begin
    @test abs((k[1])^2 - (0π/1)^2 - (0π/2)^2) ≤ 100*2*max(lat.dx,lat.dy)^2
    @test abs((k[2])^2 - (0π/1)^2 - (1π/2)^2) ≤ 100*2*max(lat.dx,lat.dy)^2
    @test abs((k[3])^2 - (0π/1)^2 - (2π/2)^2) ≤ 100*2*max(lat.dx,lat.dy)^2
    @test abs((k[4])^2 - (1π/1)^2 - (0π/2)^2) ≤ 100*2*max(lat.dx,lat.dy)^2
    @test abs((k[5])^2 - (1π/1)^2 - (1π/2)^2) ≤ 100*2*max(lat.dx,lat.dy)^2
end


text:
module Registrator

using UUIDs, LibGit2

import Base: PkgId

include("slack.jl")
include("regedit/RegEdit.jl")
include("Server.jl")

end # module


text:
import UUIDs

abstract type UnitsData end

@scoped_enum(UnitSystem, SYSTEM_BASE = 0, DEVICE_BASE = 1, NATURAL_UNITS = 2,)

const _UNIT_SYSTEM_MAP = Dict(string(x) => x for x in instances(UnitSystem))

mutable struct SystemUnitsSettings <: UnitsData
    base_value::Float64
    unit_system::UnitSystem
end

function serialize(val::SystemUnitsSettings)
    return Dict("base_value" => val.base_value, "unit_system" => string(val.unit_system))
end

function deserialize(::Type{SystemUnitsSettings}, data::Dict)
    return SystemUnitsSettings(data["base_value"], _UNIT_SYSTEM_MAP[data["unit_system"]])
end

"""
Internal storage common to InfrastructureSystems types.
"""
mutable struct InfrastructureSystemsInternal <: InfrastructureSystemsType
    uuid::Base.UUID
    units_info::Union{Nothing, UnitsData}
    ext::Union{Nothing, Dict{String, Any}}
end

"""
Creates InfrastructureSystemsInternal with a new UUID.
"""
InfrastructureSystemsInternal(; uuid=make_uuid(), units_info=nothing, ext=nothing) =
    InfrastructureSystemsInternal(uuid, units_info, ext)

"""
Creates InfrastructureSystemsInternal with an existing UUID.
"""
InfrastructureSystemsInternal(u::UUIDs.UUID) =
    InfrastructureSystemsInternal(u, nothing, nothing)

"""
Return a user-modifiable dictionary to store extra information.
"""
function get_ext(obj::InfrastructureSystemsInternal)
    if isnothing(obj.ext)
        obj.ext = Dict{String, Any}()
    end

    return obj.ext
end

"""
Clear any value stored in ext.
"""
function clear_ext!(obj::InfrastructureSystemsInternal)
    obj.ext = nothing
    return
end

get_uuid(internal::InfrastructureSystemsInternal) = internal.uuid
set_uuid!(internal::InfrastructureSystemsInternal, uuid) = internal.uuid = uuid

get_units_info(internal::InfrastructureSystemsInternal) = internal.units_info
set_units_info!(internal::InfrastructureSystemsInternal, value) =
    internal.units_info = value

"""
Gets the UUID for any InfrastructureSystemsType.
"""
function get_uuid(obj::InfrastructureSystemsType)
    return get_internal(obj).uuid
end

"""
Assign a new UUID.
"""
function assign_new_uuid!(obj::InfrastructureSystemsType)
    get_internal(obj).uuid = make_uuid()
    return
end

function serialize(internal::InfrastructureSystemsInternal)
    data = Dict{String, Any}()

    for field in fieldnames(InfrastructureSystemsInternal)
        val = getfield(internal, field)
        # reset the units data since this is a struct related to the system the components is
        # added which is resolved later in the de-serialization.
        if val isa UnitsData
            val = nothing
        else
            val = serialize(val)
        end
        data[string(field)] = val
    end

    return data
end

function compare_values(
    x::InfrastructureSystemsInternal,
    y::InfrastructureSystemsInternal;
    compare_uuids=false,
)
    match = true
    for name in fieldnames(InfrastructureSystemsInternal)
        if name == :uuid && !compare_uuids
            continue
        end
        if name == :ext
            val1 = getfield(x, name)
            if val1 isa Dict && isempty(val1)
                val1 = nothing
            end
            val2 = getfield(y, name)
            if val2 isa Dict && isempty(val2)
                val2 = nothing
            end
            if !compare_values(val1, val2, compare_uuids=compare_uuids)
                @error "ext does not match" val1 val2
                match = false
            end
        elseif !compare_values(
            getfield(x, name),
            getfield(y, name),
            compare_uuids=compare_uuids,
        )
            @error "InfrastructureSystemsInternal field=$name does not match"
            match = false
        end
    end

    return match
end

make_uuid() = UUIDs.uuid4()


text:
@safetestset "isiso31661alpha3" begin
    using Validator

    @test Validator.isiso31661alpha3("ABW") === true
    @test Validator.isiso31661alpha3("HND") === true
    @test Validator.isiso31661alpha3("KHM") === true
    @test Validator.isiso31661alpha3("RWA") === true
    @test Validator.isiso31661alpha3("rOU") === true

    @test Validator.isiso31661alpha3("") === false
    @test Validator.isiso31661alpha3("FR") === false
    @test Validator.isiso31661alpha3("fR") === false
    @test Validator.isiso31661alpha3("GB") === false
    @test Validator.isiso31661alpha3("PT") === false
    @test Validator.isiso31661alpha3("CM") === false
    @test Validator.isiso31661alpha3("JP") === false
    @test Validator.isiso31661alpha3("PM") === false
    @test Validator.isiso31661alpha3("ZW") === false
end


text:
using Tokenize
using BenchmarkTools
using Printf

function speed_test(::Type{T}=Tokenize.Tokens.Token) where T <: Tokenize.Tokens.AbstractToken
    tot_files = 0
    tot_tokens = 0
    tot_errors = 0
    basedir = dirname(Base.find_source_file("int.jl"))
    for dir in (basedir, Sys.STDLIB)
        for (root, dirs, files) in walkdir(dir)
            for file in files
                if endswith(file, ".jl")
                    tot_files += 1
                    file = joinpath(root, file)
                    str = read(file, String)::String
                    l = tokenize(str, T)
                    while !Tokenize.Lexers.eof(l)
                        t = Tokenize.Lexers.next_token(l)
                        tot_tokens += 1
                        if t.kind == Tokens.ERROR
                            tot_errors += 1
                        end
                    end
                end
            end
        end
    end
    tot_files, tot_tokens, tot_errors
end

tot_files, tot_tokens, tot_errors = speed_test()
tot_time_token = @belapsed speed_test()
tot_time_rawtoken = @belapsed speed_test(Tokenize.Tokens.RawToken)
println("Lexed ", tot_files, " files, with a total of ", tot_tokens,
        " tokens with ", tot_errors, " errors")
println("Time Token: ", @sprintf("%3.4f", tot_time_token), " seconds")
println("Time RawToken: ", @sprintf("%3.4f", tot_time_rawtoken), " seconds")


text:
using AlgebraOfGraphics
using Documenter
using Literate, Glob
using CairoMakie

CairoMakie.activate!(type="svg")

ENV["DATADEPS_ALWAYS_ACCEPT"] = true

# generate examples
GENERATED = joinpath(@__DIR__, "src", "generated")
SOURCE_FILES = Glob.glob("*.jl", GENERATED)
foreach(fn -> Literate.markdown(fn, GENERATED), SOURCE_FILES)

DocMeta.setdocmeta!(AlgebraOfGraphics, :DocTestSetup, :(using AlgebraOfGraphics); recursive=true)

makedocs(;
    modules=[AlgebraOfGraphics],
    authors="Pietro Vertechi <pietro.vertechi@veos.digital>",
    repo="https://github.com/JuliaPlots/AlgebraOfGraphics.jl/blob/{commit}{path}#{line}",
    sitename="Algebra of Graphics",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=["assets/favicon.ico"],
    ),
    pages=Any[
        "Home" => "index.md",
        "Getting Started" => [
            "generated/penguins.md",
            "generated/gallery.md",
        ],
        "Algebra of Layers" => [
            "layers/introduction.md",
            "layers/data.md",
            "layers/mappings.md",
            "Transformations" => [
                "generated/visualtransformations.md",
                "generated/datatransformations.md",
            ],
            "layers/operations.md",
            "layers/draw.md",
        ],
        "Internals" => [
            "generated/entries.md",
        ],
        "API.md",
        "philosophy.md",
    ],
    strict = true,
)

deploydocs(;
    repo="github.com/JuliaPlots/AlgebraOfGraphics.jl",
    push_preview = true,
)


text:
using Images, Luxor, BenchmarkTools

function calc(dims)
    img = fill(RGB(0,0,0), (dims, dims))
    a = 2.24; b = 0.43; c = -0.65; d = -2.43
    x = y = z = .0
    for _ in 1:dims^2
        x, y, z = sin(a * y) - z * cos(b * x), z * sin(c * x) - cos(d * y), sin(x)
        xx = trunc(Int, rescale(x, -2., 2., 0, dims))
        yy = trunc(Int, rescale(y, -2., 2., 0, dims))
        zz = rescale(z, -1, 1, .3, .8)
        img[yy,xx] = RGB(.7xx / dims, .3yy / dims + .2, zz)
    end
    return img
end

@time save("strange.png",calc(1000))

max_stars_repo_name: BoardTM/fractals

max_issues_repo_name: BoardTM/fractals

max_forks_repo_name: BoardTM/fractals


text:
"""
    SoftSphere(; cutoff, nl_only, lorentz_mixing, force_unit, energy_unit, skip_shortcut)

The soft-sphere potential.
"""
struct SoftSphere{S, C, F, E} <: GeneralInteraction
    cutoff::C
    nl_only::Bool
    lorentz_mixing::Bool
    force_unit::F
    energy_unit::E
end

function SoftSphere(;
                    cutoff=NoCutoff(),
                    nl_only=false,
                    lorentz_mixing=true,
                    force_unit=u"kJ * mol^-1 * nm^-1",
                    energy_unit=u"kJ * mol^-1",
                    skip_shortcut=false)
    return SoftSphere{skip_shortcut, typeof(cutoff), typeof(force_unit), typeof(energy_unit)}(
        cutoff, nl_only, lorentz_mixing, force_unit, energy_unit)
end

@inline @inbounds function force(inter::SoftSphere{S, C},
                                    coord_i,
                                    coord_j,
                                    atom_i,
                                    atom_j,
                                    box_size) where {S, C}
    dr = vector(coord_i, coord_j, box_size)
    r2 = sum(abs2, dr)

    if !S && iszero(atom_i.σ) || iszero(atom_j.σ)
        return ustrip.(zero(coord_i)) * inter.force_unit
    end

    # Lorentz-Berthelot mixing rules use the arithmetic average for σ
    # Otherwise use the geometric average
    σ = inter.lorentz_mixing ? (atom_i.σ + atom_j.σ) / 2 : sqrt(atom_i.σ * atom_j.σ)
    ϵ = sqrt(atom_i.ϵ * atom_j.ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        f = force_divr_nocutoff(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_unit

        f = force_divr_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip.(zero(coord_i)) * inter.force_unit

        if r2 < cutoff.activation_dist
            f = force_divr_nocutoff(inter, r2, inv(r2), params)
        else
            f = force_divr_cutoff(cutoff, r2, inter, params)
        end
    end

    return f * dr
end

@fastmath function force_divr_nocutoff(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return (24ϵ * invr2) * 2 * six_term ^ 2
end

@inbounds function potential_energy(inter::SoftSphere{S, C},
                                    s::System,
                                    i::Integer,
                                    j::Integer) where {S, C}
    dr = vector(s.coords[i], s.coords[j], s.box_size)
    r2 = sum(abs2, dr)

    if !S && iszero(s.atoms[i].σ) || iszero(s.atoms[j].σ)
        return ustrip(zero(s.box_size[1])) * inter.energy_unit
    end

    σ = inter.lorentz_mixing ? (s.atoms[i].σ + s.atoms[j].σ) / 2 : sqrt(s.atoms[i].σ * s.atoms[j].σ)
    ϵ = sqrt(s.atoms[i].ϵ * s.atoms[j].ϵ)

    cutoff = inter.cutoff
    σ2 = σ^2
    params = (σ2, ϵ)

    if cutoff_points(C) == 0
        potential(inter, r2, inv(r2), params)
    elseif cutoff_points(C) == 1
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(s.box_size[1])) * inter.energy_unit

        potential_cutoff(cutoff, r2, inter, params)
    elseif cutoff_points(C) == 2
        r2 > cutoff.sqdist_cutoff && return ustrip(zero(s.box_size[1])) * inter.energy_unit

        if r2 < cutoff.activation_dist
            potential(inter, r2, inv(r2), params)
        else
            potential_cutoff(cutoff, r2, inter, params)
        end
    end
end

@fastmath function potential(::SoftSphere, r2, invr2, (σ2, ϵ))
    six_term = (σ2 * invr2) ^ 3

    return 4ϵ * (six_term ^ 2)
end


text:
module StingerGraphs

# package code goes here
include("stinger_lib.jl")
include("stinger_graph.jl")
include("stinger_config.jl")
include("stinger_core.jl")
include("fields.jl")
include("traversal.jl")

include("algorithms/bfs.jl")
include("algorithms/parallelbfs.jl")
include("algorithms/kcore.jl")
include("generators/kronecker.jl")

end # module


text:
### Resample Interface ###
"""
    resample(resampler, bp::AbstractParticleBelief, rng::AbstractRNG)

Sample a new ParticleCollection from `bp`.

Generic domain-independent resamplers should implement this version.

    resample(resampler, bp::WeightedParticleBelief, predict_model, reweight_model, b, a, o, rng)

Sample a new particle collection from bp with additional information from the arguments to the update function.

This version defaults to `resample(resampler, bp, rng)`. Domain-specific resamplers that wish to add noise to particles, etc. should implement this version.
"""
function resample end

resample(resampler, bp::WeightedParticleBelief, pm, rm, b, a, o, rng) = resample(resampler, bp, rng)

function resample(resampler, bp::WeightedParticleBelief, pm::Union{POMDP,MDP}, rm, b, a, o, rng)
    if weight_sum(bp) == 0.0 && all(isterminal(pm, s) for s in particles(b))
        error("Particle filter update error: all states in the particle collection were terminal.")
    end
    resample(resampler, bp, rng)
end

### Resamplers ###

abstract type Resampler end

"""
    ImportanceResampler(n)

Simple resampler. Uses alias sampling to attain O(n log(n)) performance with uncorrelated samples.
"""
struct ImportanceResampler <: Resampler
    n::Int
end

function resample(r::ImportanceResampler, b::AbstractParticleBelief{S}, rng::AbstractRNG) where {S}
    ps = Array{S}(undef, r.n)
    if weight_sum(b) <= 0
        warn("Invalid weights in particle filter: weight_sum = $(weight_sum(b))")
    end
    #XXX this may break if StatsBase changes
    StatsBase.alias_sample!(rng, particles(b), Weights(weights(b), weight_sum(b)), ps)
    return ParticleCollection(ps)
end

"""
    LowVarianceResampler(n)

Low variance sampling algorithm on page 110 of Probabilistic Robotics by Thrun Burgard and Fox. O(n) runtime, correlated samples, but produces a useful low-variance set.
"""
struct LowVarianceResampler <: Resampler
    n::Int
end

function resample(re::LowVarianceResampler, b::AbstractParticleBelief{S}, rng::AbstractRNG) where {S}
    ps = Array{S}(undef, re.n)
    r = rand(rng)*weight_sum(b)/re.n
    c = weight(b,1)
    i = 1
    U = r
    for m in 1:re.n
        while U > c && i < n_particles(b)
            i += 1
            c += weight(b, i)
        end
        U += weight_sum(b)/re.n
        ps[m] = particles(b)[i]
    end
    return ParticleCollection(ps)
end

function resample(re::LowVarianceResampler, b::ParticleCollection{S}, rng::AbstractRNG) where {S}
    r = rand(rng)*n_particles(b)/re.n
    chunk = n_particles(b)/re.n
    inds = ceil.(Int, chunk*(0:re.n-1).+r)
    ps = particles(b)[inds]
    return ParticleCollection(ps)
end

n_init_samples(r::Union{LowVarianceResampler, ImportanceResampler}) = r.n

resample(f::Function, d::Any, rng::AbstractRNG) = f(d, rng)

resample(re, d::Any, rng::AbstractRNG) = ParticleCollection([rand(rng, d) for i=1:re.n])


text:
module VisualEFM

using Plots
gr();
using Images
using ImageFiltering
using ImageBinarization
using ImageMorphology

include("utils.jl")
include("sandErosion.jl")
include("smoke.jl")

export
    binarymask,
    read_image,
    erosion_one,
    erosion_anim,
    erosion_colormap,
    smoke_anim,
    smoke_statsmap

end


text:
using Revise
using Plots
using Basins
using DifferentialEquations
using DynamicalSystems

# Blowout bifurcation with non-normal parameters in population dynamicsBernard Cazelle
function franke_yabuku(dz,z, p, n)
    r=2.825; s=20.25; c1=1.2; c2=0.1;
    x = z[1]; y = z[2];
    dx = x*exp(r-s*(x+y))
    dy = c1*y/(c2+x+y)
    dz = [dx,dy]
    return
end

# dummy function to keep the initializator happy
function franke_yabuku_J(J,z0, p, n)
   return
end

ds = DiscreteDynamicalSystem(franke_yabuku,[0.1, 0.2], [] , franke_yabuku_J)
integ  = integrator(ds)

xg=range(8.,16.,length=100)
yg=range(0.5,1.,length=100)

# compute basin

@time basin=basins_map2D(xg, yg, integ)
plot(xg,yg,basin',seriestype=:heatmap)


text:
module V

const NAMES = [
    "ppMEKc"
    "CREBn"
    "pCREBn"
    "ERKc"
    "ERKn"
    "pERKc"
    "pERKn"
    "ppERKc"
    "ppERKn"
    "Elk1n"
    "pElk1n"
    "cFOSc"
    "cFOSn"
    "pcFOSc"
    "pcFOSn"
    "DUSPc"
    "DUSPn"
    "pDUSPc"
    "pDUSPn"
    "DUSPn_ERKn"
    "DUSPn_pERKn"
    "DUSPn_ppERKn"
    "pDUSPn_ERKn"
    "pDUSPn_pERKn"
    "pDUSPn_ppERKn"
    "RSKc"
    "pRSKc"
    "pRSKn"
    "PrecfosmRNAn"
    "PreduspmRNAn"
    "cfosmRNAc"
    "duspmRNAc"
    "Fc"
    "Fn"
    "FmRNAc"
    "PreFmRNAn"
]

for (idx,name) in enumerate(NAMES)
    eval(Meta.parse("const $name = $idx"))
end

const NUM = length(NAMES)

end  # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim


text:
## test_circulant_embedding.jl : test circulant embedding for GRF generation

@testset "circulant embedding      " begin

    ## 1d Exponential ##
    cov = CovarianceFunction(1,Exponential(0.1))
    pts = 0:1//1001:1
    grf = GaussianRandomField(cov,CirculantEmbedding(),pts)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Exponential)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{CirculantEmbedding})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## 1d Matern with padding ##
    m = CovarianceFunction(1,Matern(1.0,1.0))
    grf = GaussianRandomField(m,CirculantEmbedding(),pts,minpadding=8000)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,Matern)
    @test ndims(grf.cov) == 1
    @test isa(grf,GaussianRandomField{CirculantEmbedding})
    @test length(grf.pts[1]) == length(pts)
    @test length(sample(grf)) == length(pts)

    ## 2d Anisotropic ##
    A = [1000 0; 0 1000]
    a = AnisotropicExponential(A)
    c = CovarianceFunction(2,a)
    pts1 = range(0,stop = 1,length = 128)
    pts2 = range(0,stop = 1,length = 256)
    @suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2,minpadding=(128,256))
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,AnisotropicExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{CirculantEmbedding})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## 3d AnisotropicExponential ##
    A = [100 0 0; 0 100 0; 0 0 100]
    a = AnisotropicExponential(A)
    c = CovarianceFunction(3,a)
    pts = 0:0.025:1
    @suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts,pts,pts,measure=false)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,AnisotropicExponential)
    @test ndims(grf.cov) == 3
    @test isa(grf,GaussianRandomField{CirculantEmbedding})
    @test length(grf.pts) == 3
    @test length(grf.pts[1]) == length(pts)
    @test length(grf.pts[2]) == length(pts)
    @test length(grf.pts[3]) == length(pts)
    @test size(sample(grf),1) == length(pts)
    @test size(sample(grf),2) == length(pts)
    @test size(sample(grf),3) == length(pts)

    ## Larger domain ##
    A = [1 0.8; 0.8 1]
    m = AnisotropicExponential(A)
    c = CovarianceFunction(2,m)
    pts1 = range(-5,stop = 5,length = 128)
    pts2 = range(10,stop = 0,length = 128)
    @suppress grf = GaussianRandomField(c,CirculantEmbedding(),pts1,pts2,minpadding=128,primes=false)
    @test isa(grf,GaussianRandomField)
    @test isa(grf.cov,CovarianceFunction)
    @test isa(grf.cov.cov,AnisotropicExponential)
    @test ndims(grf.cov) == 2
    @test isa(grf,GaussianRandomField{CirculantEmbedding})
    @test length(grf.pts) == 2
    @test length(grf.pts[1]) == length(pts1)
    @test length(grf.pts[2]) == length(pts2)
    @test size(sample(grf),1) == length(pts1)
    @test size(sample(grf),2) == length(pts2)

    ## Dedicated plotting commands ##
    m = Matern(0.5,2)
    c = CovarianceFunction(3,m)
    pts1 = range(0,stop = 1,length = 16)
    pts2 = range(0,stop = 1,length = 8)
    pts3 = range(0,stop = 1,length = 4)
    #plot_covariance_matrix(c,pts1,pts2,pts3); close()
    m = Matern(0.1,1)
    c = CovarianceFunction(2,m)
    pts = range(0,stop = 1,length = 512)
    g = GaussianRandomField(c,CirculantEmbedding(),pts,pts,minpadding=512,measure=false)
    #contourf(g); close()
    #plot_eigenvalues(g); close()

end


text:
abstract type FSNode end
abstract type File <: FSNode end

struct CommonNode
    name::String
    created::String
    modified::String
    thumb::Bool
    ismine::Bool
    isshared::Bool
    id::String
    parentfolderid::Int
    comments::String
    metadata::NamedTuple
end

struct CommonFile
    fileid::Int
    hash::String
    category::Int
    size::Int
    contenttype::String
    icon::String
end

struct Folder <: FSNode
    folderid::Int
    g::CommonNode
    icon::String
    contents::Vector{FSNode}
end

struct Plain <: File
    f::CommonFile
    g::CommonNode
end

struct Document <: File
    f::CommonFile
    g::CommonNode
end

struct Image <: File
    f::CommonFile
    g::CommonNode
    width::Int
    height::Int
end

# TODO: Add necessary fields
struct Video <: File
    f::CommonFile
    g::CommonNode
end

# TODO: Add necessary fields
struct Audio <: File
    f::CommonFile
    g::CommonNode
end

# TODO: Add necessary fields
struct Archive <: File
    f::CommonFile
    g::CommonNode
end

name(x::FSNode) = x.g.name
created(x::FSNode) = x.g.created
modified(x::FSNode) = x.g.modified
thumb(x::FSNode) = x.g.thumb
ismine(x::FSNode) = x.g.ismine
isshared(x::FSNode) = x.g.isshared
id(x::FSNode) = x.g.id
parentfolderid(x::FSNode) = x.g.parentfolderid
comments(x::FSNode) = x.g.comments
metadata(x::FSNode) = x.g.metadata

folderid(x::Folder) = x.folderid

fileid(x::File) = x.f.fileid
size(x::File) = x.f.size
hash(x::File) = x.f.hash
category(x::File) = x.f.category
contenttype(x::File) = x.f.contenttype
icon(x::File) = x.f.icon

isfolder(x::Folder) = true
isfolder(x::File) = false

catname(x::Plain) = "Uncategorized"
catname(x::Image) = "Image"
catname(x::Video) = "Video"
catname(x::Audio) = "Audio"
catname(x::Document) = "Document"
catname(x::Archive) = "Archive"

Base.:length(x::Folder) = length(x.contents)
Base.:iterate(x::Folder) = iterate(x.contents)
Base.:iterate(x::Folder, i::Int) = iterate(x.contents, i)

function CommonNode(x)
    properties = propertynames(x)
    parentfolderid = :parentfolfolderid in properties ? x.parentfolderid : -1
    comments = :comments in properties ? string(x.comments) : ""
    y = []
    for property in propertynames(x)
        property == :contents && continue
        push!(y, property => getproperty(x, property))
    end
    metadata = (; y...)
    CommonNode(x.name, x.created, x.modified, x.thumb, x.ismine, x.isshared, x.id, parentfolderid, comments, metadata)
end

function CommonFile(x)
    CommonFile(x.fileid, string(x.hash), x.category, x.size, x.contenttype, x.icon)
end

function Folder(x, g::CommonNode)
    properties = propertynames(x)
    contents = :contents in properties ? FSNode.(x.contents) : FSNode[]
    Folder(x.folderid, g, x.icon, contents)
end

function Plain(x, g::CommonNode)
    Plain(CommonFile(x), g)
end

function Document(x, g::CommonNode)
    Document(CommonFile(x), g)
end

function Image(x, g::CommonNode)
    Image(CommonFile(x), g, x.width, x.height)
end

function Video(x, g::CommonNode)
    Video(CommonFile(x), g)
end

function Audio(x, g::CommonNode)
    Audio(CommonFile(x), g)
end

function Archive(x, g::CommonNode)
    Archive(CommonFile(x), g)
end

function FSNode(x)
    g = CommonNode(x)
    x.isfolder && return Folder(x, g)
    x.category == 0 && return Plain(x, g)
    x.category == 1 && return Image(x, g)
    x.category == 2 && return Video(x, g)
    x.category == 3 && return Audio(x, g)
    x.category == 4 && return Document(x, g)
    x.category == 5 && return Archive(x, g)

    @warn "Unknown category " * string(x.category)
    return Plain(x, g)
end

children(x::Folder) = x.contents
children(x::File) = ()

printnode(io::IO, x::Folder) = print(io, folderid(x), ": " * name(x) * "/")
printnode(io::IO, x::File) = print(io, fileid(x), " (", catname(x), ")", ": ", name(x))

function gettree(client::PCloudClient = OPTS.client; params...)
    if !haskey(params, :folderid)
        params = (; params..., :folderid => 0)
    end

    dir = API.listfolder(client; params...)
    return FSNode(dir.metadata)
end

function gettree!(client::PCloudClient, node::Folder; params...)
    dir = API.listfolder(client; :folderid => folderid(node), params...)
    empty!(node.contents)
    if haskey(dir.metadata, :contents)
        for x in dir.metadata.contents
            push!(node.contents, FSNode(x))
        end
    end

    return
end

gettree!(client::PCloudClient, node::File; params...) = nothing

gettree!(nodes::AbstractVector; params...) = gettree!(OPTS.client, nodes; params...)
function gettree!(client::PCloudClient, nodes::AbstractVector; params...)
    for node in nodes
        gettree!(client, node; params...)
    end

    return
end

function search(f, x::FSNode)
    res = FSNode[]
    for node in PreOrderDFS(x)
        f(node) && push!(res, node)
    end

    return res
end


text:
using DelimitedFiles
using UnicodePlots

function analyze_flow(xs0, logp0, xs, logp; xlim=[-5,5], ylim=[0,1])
    order = sortperm(xs0)[1:10:end]
    lineplot(xs0[order], exp.(logp0[order]); xlim=xlim, ylim=ylim,title="p(x), before transformation") |> println

    fy = fit(Histogram, xs, range(xlim..., length=200), closed=:left)
    xticks = (fy.edges[1][1:end-1] .+ fy.edges[1][2:end]) ./ 2
    lineplot(xticks, fy.weights; xlim=xlim, title="data, after transformation") |> println

    order = sortperm(xs)[1:10:end]
    lineplot(xs[order], exp.(logp[order]); xlim=xlim, ylim=ylim,title="p(x), after transformation") |> println
end

function forward_sample(μ, σ, field, θ; nsample=100, Nt=100, dt=0.01, xlim=[-5,5], ylim=[0,1])
    # source distribution
    source_distri = Normal(μ, σ)
    xs0 = rand(source_distri, nsample)
    logp0 = max.(logpdf.(source_distri, xs0), -10000)

    # solving the ode, and obtain the probability change
    @newvar v_xs = copy(xs0)
    @newvar v_ks = copy(xs0)
    @newvar v_logp = copy(logp0)
    @newvar v_θ = θ
    @newvar field_out = 0.0
    tape = ode!(v_xs, v_ks, v_logp, field, field_out, v_θ, Nt, dt)

    play!(tape)
    analyze_flow(xs0, logp0, v_xs[], v_logp[]; xlim=xlim)
end

function get_loss_gradient!(forward, loss_out, θ)
    reset_grad!(forward)
    resetreg(forward)
    play!(forward)

    ll = loss_out[]

    grad(loss_out)[] = 1
    play!(forward')
    return ll, grad(θ)
end

function train(xs_target, θ; niter=100, Nt=100, dt=0.01, lr=0.1)
    μ = 0.0
    σ = 1.0
    @newvar loss_out = 0.0
    @newvar logp_out = zeros(Float64, length(xs_target))
    @newvar v_xst = xs_target
    @newvar v_kst = copy(xs_target)
    @newvar field_out = 0.0
    @newvar jacobian_out = zero(logp_out[])
    forward = ode_loss!(μ, σ, logp_out, v_xst, v_kst, field, ode!,
        inv_normal_logpdf!, jacobian_out, field_out, loss_out, θ, Nt, dt)

    init_grad!(forward)
    local ll
    for i=1:niter
        logp_out[] .= 0.0
        loss_out[] = 0.0
        ll, θδ = get_loss_gradient!(forward, loss_out, θ)
        println("Step $i, log-likelihood = $ll")
        θ[] += θδ[] * lr
        @show θ[]
    end
    return ll, forward
end

=#


text:
struct FDSwell <: FDSpectralModel
    θ::Vector{Float64}
end

parameter(F::FDSwell) = F.θ

function sdf(F::FDSwell, ω)
    θ = F.θ
    α = θ[1]
    ωₚ = θ[2]
    κ = θ[3]
    sigω = sign(ω)
    ω = abs(ω)
    logbit = log(ω) - log(ωₚ)
    h = logbit/κ - κ
    sdf = α * exp(-0.5 * (h)^2)
    sdf *= (!isnan(sdf))
end

spreading(F::FDSwell, ω, ϕ) = WrappedGaussian(ϕ, F.θ[4], F.θ[5])

####

struct FDSwellJS <: FDSpectralModel
    θ::Vector{Float64}
end

parameter(F::FDSwellJS) = F.θ

spreading(F::FDSwellJS, ω, ϕ) = WrappedGaussian(ϕ, F.θ[5], F.θ[6])


text:
# common single qubit gate defns
export Id, X, Y, Z, H, S, T, Rn, Rx, Ry, Rz

"""
    Id

Single qubit Id (identity) := ComplexF64[1. 0.; 0. 1.]
"""
Id = ComplexF64[1. 0.; 0. 1.]

"""
    X

Single qubit X (not) := ComplexF64[0 1; 1 0]
"""
X = ComplexF64[0 1; 1 0]

"""
    Z

Single qubit Z (phase-flip) := ComplexF64[1 0; 0 -1]
"""
Z = ComplexF64[1 0; 0 -1]

"""
    Y

Single qubit Y (π rotation about Y-axis) := ComplexF64[1 0; 0 -1]
"""
Y = ComplexF64[0 -1im; 1im 0]

"""
    H

Single qubit H `(Hadamard - single qubit QFT)` := ComplexF64[1 1; 1 -1] ./ sqrt(2)
"""
H = ComplexF64[1 1; 1 -1] ./ sqrt(2)

"""
    S

Single qubit S `(== Z^2)` := ComplexF64[1 0; 0 im]
"""
S = ComplexF64[1 0; 0 im]

"""
    T

Single qubit T `(π/8 gate, == Z^4)` := ComplexF64[1 0; 0 exp(im * π/4)]
"""
T = ComplexF64[1 0; 0 exp(im * π/4)]

"""
    Rn(θ, n̂)

Rn returns an operator that rotates a qubit about the n̂ axis.
Given an angle [0, 2π) in radians, and a normed vector n̂,
"""
Rn(θ, n̂) = cos(θ/2) .* Id - 1im*sin(θ/2) .* (n̂[1] * X + n̂[2] * Y + n̂[3] * Z)

# TODO: generate instead of manual
Rx(θ) = cos(θ/2) .* Id - 1im*sin(θ/2) .* X
Ry(θ) = cos(θ/2) .* Id - 1im*sin(θ/2) .* Y
Rz(θ) = cos(θ/2) .* Id - 1im*sin(θ/2) .* Z

max_stars_repo_name: anandijain/MyQuantum

max_issues_repo_name: anandijain/MyQuantum

max_forks_repo_name: anandijain/MyQuantum


text:
# The Orthopolys packages is not in the package repo
hafta_installed = false
try
    import Hafta
    println("Hafta.jl is already available.")
    hafta_installed = true
catch end

if ! hafta_installed
    Pkg.clone("https://github.com/mortenpi/Hafta.jl.git")
    Pkg.build("Hafta")
end


text:
module ReflectionFg

using ..Ahorn, Maple

placements = Maple.ReflectionFg

function Ahorn.canFgBg(effect::Maple.ReflectionFg)
    return true, false
end

end


text:
#Heatmap plots

import .LibCImPlot: PlotHeatmap

function PlotHeatmap(x::AbstractArray{T}, rows, cols, scale_min = 0.0, scale_max = 1.0;
                 label_id::String = "", label_fmt="%.1f", bounds_min = LibCImPlot.ImPlotPoint(0.0,0.0),
                 bounds_max = LibCImPlot.ImPlotPoint(1.0,1.0)) where {T<:ImPlotData}

    LibCImPlot.PlotHeatmap(label_id, x, rows, cols, scale_min, scale_max, label_fmt, bounds_min, bounds_max)
end


text:
## Binomial Example ####
## Note: Adapted from the Rate_4 example in Bayesian Cognitive Modeling
##  https://github.com/stan-dev/example-models/tree/master/Bayesian_Cognitive_Modeling

using Mamba, Stan

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  binomialstanmodel = "
  // Inferring a Rate
  data {
    int<lower=1> n;
    int<lower=0> k;
  }
  parameters {
    real<lower=0,upper=1> theta;
    real<lower=0,upper=1> thetaprior;
  }
  model {
    // Prior Distribution for Rate Theta
    theta ~ beta(1, 1);
    thetaprior ~ beta(1, 1);

    // Observed Counts
    k ~ binomial(n, theta);
  }
  generated quantities {
    int<lower=0> postpredk;
    int<lower=0> priorpredk;

    postpredk <- binomial_rng(n, theta);
    priorpredk <- binomial_rng(n, thetaprior);
  }
  "

  global stanmodel, rc, sim1
  stanmodel = Stanmodel(name="binomial", model=binomialstanmodel)

  binomialdata = [
    Dict("n" => 10, "k" => 5)
  ]

  rc, sim1 = stan(stanmodel, binomialdata, ProjDir, diagnostics=false,
    CmdStanDir=CMDSTAN_HOME)

  if rc == 0
    theta_sim = sim1[1:size(sim1,1), ["theta", "thetaprior"], :]
    predk_sim = sim1[1:size(sim1,1), ["postpredk", "priorpredk"], :]

    p = plot(theta_sim, [:trace, :mean, :density, :autocor], legend=false);
    draw(p, ncol=4, nrow=2, filename="$(stanmodel.name)-thetas", fmt=:svg)

    p = plot(predk_sim, [:bar], legend=false)
    draw(p, ncol=2, nrow=2, filename="$(stanmodel.name)-predictiv", fmt=:svg)
  end # rc == 0
end # cd


text:
using CommonLH, CollegeStratCollege
using Random, Test

cs = CollegeStratCollege;

function hprod_main_test(hS ::cs.AbstractHcProdFct)
	@testset "Basics $hS" begin
		println(hS);

		st = study_time_per_course(hS, 0.3, 3);
		@test isa(st, Float64)
		stV = study_time_per_course(hS, [0.3, 0.4], 3);
		@test isa(stV, Vector{Float64})
		stV = study_time_per_course(hS, 0.3, [3,4]);
		@test isa(stV, Vector{Float64})
		@test all(stV .> 0.0)
	end
end

function tfp_test(learnRelativeToH0 :: Bool, tfpSpec :: AbstractTfpSpec)
	@testset "TFP $learnRelativeToH0, $tfpSpec" begin
		hS = cs.make_test_hprod_bounded(
			learnRelativeToH0 = learnRelativeToH0, tfpSpec = tfpSpec);

		# Set h values so that learning is either 0 or maxLearn
		h0V = [2.1, 3.4];
		maxLearn = cs.max_learn(hS);
		if learnRelativeToH0
			h2 = h0V[2] * (1.0 + maxLearn);
		else
			h2 = h0V[2] + maxLearn;
		end
		hV = [h0V[1], h2];

		learnV = cs.learned_h(hS, hV, h0V);
		@test learnV[1] ≈ 0.0;
		@test learnV[2] ≈ maxLearn;

		tfpV = cs.base_tfp(hS, hV, h0V);
		tfpRange = cs.tfp_range(hS);
		@test all(tfpV .> tfpRange[1] .- 1e-6);
		@test all(tfpV .<= tfpRange[2] .+ 1e-6);
		@test all(tfpV .>= 0.0);
	end
end

function hprod_dh_test(hS ::cs.AbstractHcProdFct)
	@testset "dh $hS" begin
		nc = 4;

		n = 3;
		hV = collect(range(2, 3, length = n));
		h0V = hV .* range(0.5, 1.5, length = n);
		timeV = collect(range(0.4, 0.2, length = n));
		nBarV = round.(UInt8, collect(range(1, 4, length = n)));
		abilV = collect(range(-0.3, 0.4, length = n));
		dhV = cs.dh(hS, abilV,  hV,  h0V,  timeV,  nBarV);
		@test all(dhV .>= 0.0)
		@test size(dhV) == size(hV)
		@test eltype(dhV) == Float64

		dhV = cs.dh(hS, 0.83,  hV,  h0V,  timeV,  nBarV);
		@test all(dhV .> 0.0)
		@test isa(dhV, Vector{Float64})

		dh = cs.dh(hS, 0.83, 1.2, 0.9, 0.3, 3);
		@test isa(dh, Float64)
	end
end

function dh_ben_porath_test()
	@testset "dh Ben Porath" begin
		hS = cs.make_test_hprod();
		# TFP increases dh
		set_tfp(hS, 1.2);
		dh1 = cs.dh(hS, 0.8, 1.2, 0.9, 0.3, 3);
		set_tfp(hS, 1.3);
		dh2 = cs.dh(hS, 0.8, 1.2, 0.9, 0.3, 3);
		@test dh2 - dh1 > 0.001
    end
end

# For a set of colleges
function hprod_set_test(hs :: cs.AbstractHcProdSet)
	@testset "Set $hs" begin
		nc = 3;
		for ic = 1 : nc
			hp = cs.make_h_prod(hs, ic);
			@test cs.validate_hprod(hp)
			# @test CollegeStrat.tfp(hs, ic) ≈ hp.tfp
		end
	end
end

function h_path_test(hS :: cs.AbstractHcProdFct)
	@testset "H path $hS" begin
		# hS = CollegeStrat.make_test_hprod();
		T = 6;
		abil = -0.8;
		h0 = 1.2;
		sTimeV = collect(range(0.4, 0.1, length = T));
		nTriedV = round.(Int, range(5, 15, length = T));
		hV = h_path(hS, abil, h0, sTimeV, nTriedV);
		@test check_float_array(hV, h0, 100.0 * h0)
		@test length(hV) == T+1

		# Increasing TFP or max learning should increase h
		if cs.has_tfp(hS)
			tfp0 = tfp(hS);
			set_tfp(hS, tfp0 + 0.1);
			@test tfp(hS) ≈ tfp0 + 0.1
			h2V = h_path(hS, abil, h0, sTimeV, nTriedV);
			@test h2V[1] ≈ hV[1]
			@test all(h2V[2 : (T+1)] .- hV[2 : (T+1)] .> 0.001)
		end
	end
end

@testset "H production" begin
	for tfpSpec in cs.tfp_spec_list()
		for learnRelativeToH0 in (true, false)
			tfp_test(learnRelativeToH0, tfpSpec);
		end
	end

	for h ∈ [
		cs.make_test_hprod(),
		cs.make_test_hprod_bounded(),
		cs.make_test_hprod_ces(cs.hCesAggrAhl),
		cs.make_test_hprod_ces(cs.hCesAggrA)
		]

		hprod_main_test(h);
		hprod_dh_test(h);
		h_path_test(h);
	end
	dh_ben_porath_test();
	for hS ∈ [
		cs.make_test_hc_prod_set(),
		cs.make_test_hc_bounded_set(learnRelativeToH0 = false),
		cs.make_test_hc_bounded_set(learnRelativeToH0 = true),
		cs.make_test_hc_bounded_set(learnRelativeToH0 = true,
			tfpSpec = cs.TfpOneMinusLearnOverMaxLearn()),
		cs.make_test_hc_ces_set(cs.hCesAggrAhl),
		cs.make_test_hc_ces_set(cs.hCesAggrA)
		]

		hprod_set_test(hS);
	end
end

# ----------------


text:
using PushGP: Stack

@testset "PushStack" begin

@testset "Stack{Int}" begin
    ist = Stack{Int}()
    @test length(ist) == 0

    @test push!(ist, 1) == Int[1]
    @test length(ist) == 1
    @test last(ist) == 1

    push!(ist, 2)
    push!(ist, 3)
    @test length(ist) == 3

    a, b = pop2!(ist)
    @test a == 3
    @test b == 2

    # Can't push a Float on an Int stack
    @test_throws ErrorException push!(ist, 1.0)
end

@testset "Stack{Float64}" begin
    fst = Stack{Float64}()
    @test length(fst) == 0

    @test push!(fst, 3.14) == Float64[3.14]
    @test length(fst) == 1
    @test last(fst) == 3.14

    push!(fst, 2e-3)
    push!(fst, 3.1e2)
    @test length(fst) == 3

    a, b = pop2!(fst)
    @test a == 3.1e2
    @test b == 2e-3

    # Can't push a Float on an Int stack
    @test_throws ErrorException push!(fst, 1)
end # @testset "Stack{Float64}" begin

@testset "Stack{Bool}" begin
    bst = Stack{Bool}()
    @test length(bst) == 0

    push!(bst, false)
    push!(bst, false)
    push!(bst, true)
    @test length(bst) == 3

    a, b = pop2!(bst)
    @test a === true
    @test b === false

    # Can't push a Float on an Int stack
    @test_throws ErrorException push!(bst, 1)
    @test_throws ErrorException push!(bst, 3.14)
end # @testset "Stack{Bool}" begin

@testset "Stack{String}" begin
    sst = Stack{String}()
    @test length(sst) == 0

    push!(sst, "ab")
    push!(sst, "cde")
    push!(sst, "f")
    @test length(sst) == 3

    a, b = pop2!(sst)
    @test a === "f"
    @test b === "cde"

    @test length(sst) == 1
    empty!(sst)
    @test length(sst) == 0

    # Can't push a Float on an Int stack
    @test_throws ErrorException push!(sst, 1)
    @test_throws ErrorException push!(sst, 3.14)
    @test_throws ErrorException push!(sst, true)
end # @testset "Stack{String}" begin

end # @testset "PushStack" begin


text:
# direct integration algorithms
export Bathe,
       CentralDifference,
       Houbolt,
       Linear,
       Newmark,
       BackwardEuler,
       Trapezoidal

# solution API
export solve,
       displacements,
       velocities,
       accelerations,
       times,
       dim


text:
# r[3,Npoints]
# atpos[3,Natoms]
# A[Natoms]
# R0[Natoms]
function to_curvilinear( r::Array{Float64,2}, atpos, A, R0)
    #
    Npoints = size(r)[2]
    Natoms = size(atpos)[2]
    #
    ξ = zeros(Float64,3,Npoints)
    #
    for ip = 1:Npoints
        ξ[:,ip] = r[:,ip]
        for ia = 1:Natoms
            xx = norm( (r[:,ip] - atpos[:,ia])/R0[ia] )
            ξ[:,ip] = ξ[:,ip] + A[ia]*( r[:,ip] - atpos[:,ia] )/R0[ia] * sech(xx)
        end
    end
end


text:
"""
    create_geometry(setup)

Create cells, surfaces and finite element mesh.
This function does the following:
- Check geometry setup consistency
- Create or load cell configuration
- Create or load surface triangulation
- Call TetGen
- Deform domain
- Split mesh into compartments
For custom geometries with more than one compartment, call `split_mesh`
directly instead. This requires facet and element labels.
"""
function create_geometry(setup::Setup)

    @unpack ncell, ecs_shape, refinement = setup

    # File name for saving or loading geometry
    filename = "meshfiles" * "/" * setup.name

    # Make sure that folder exists
    dir = dirname(filename)
    isdir(dir) || mkpath(dir)

    # Check correct input format
    if isa(setup, NeuronSetup) && ncell != 1
        error("Neuron cell type is only available for ncell=1.")
    end
    !(ecs_shape == "no_ecs" && ncell > 1) ||
        @error "Geometry must include ECS if more than one cell"
    @assert ecs_shape ∈ ["no_ecs", "box", "convex_hull", "tight_wrap"]

    # Check if cell description file is already available
    cellfilename = filename * "_cells"
    if isfile(cellfilename)
        cells = read_cells(cellfilename)
    elseif isa(setup, SphereSetup) || isa(setup, CylinderSetup)
        # Create cells
        cells = create_cells(setup)

        # Save cell configuration
        save_cells(cells, cellfilename)
    else
        # No cells
        cells = nothing
    end

    # Make directory for storing finite elements mesh
    is_stl = endswith(filename, ".stl")
    if is_stl
        # ECS is currently only available for surface meshes
        @assert ecs_shape == "no_ecs"
        filename = filename[1:end-4]
    end
    save_meshdir_path = filename * "_dir"
    if !isdir(save_meshdir_path)
        mkpath(save_meshdir_path)
    end

    # Use an existing finite elements mesh or create a new finite
    # elements mesh. The name of the finite elements mesh is stored in the string
    # fname_tetgen_femesh
    refinement_str = isnothing(refinement) ? "" : "_refinement$refinement"
    fname_tetgen =
        save_meshdir_path * "/" * split(filename, "/")[end] * refinement_str * "_mesh"

    # Read or create surface triangulation
    if isfile(fname_tetgen * ".node") && isfile(fname_tetgen * ".poly")
        surfaces = read_surfaces(fname_tetgen)
    else
        if isa(setup, SphereSetup)
            # Create surface geometry of spheres
            surfaces = create_surfaces_sphere(cells, setup)
        elseif isa(setup, CylinderSetup)
            # Create surface geometry of cylinders
            surfaces = create_surfaces_cylinder(cells, setup)
        elseif isa(setup, NeuronSetup)
            if is_stl
                surfaces = nothing
            else
                surfaces = create_surfaces_neuron(filename, setup)
            end
        end

        if !is_stl
            save_surfaces(fname_tetgen, surfaces)
        end
    end

    # Add ".1" suffix to output file name, since this is what Tetgen does
    fname_tetgen_femesh = fname_tetgen * ".1"

    if isfile(fname_tetgen_femesh * ".node")
        # Read global mesh from Tetgen output
        mesh_all = read_tetgen(fname_tetgen_femesh)
    elseif is_stl
        error("Not implemented")
    else
        mesh_all = call_tetgen(surfaces, refinement)
        save_tetgen(mesh_all, fname_tetgen_femesh)
    end

    # # Check that at correct number of compartments and boundaries has been found
    # compartments_new = unique(mesh_all.elementmarkers)
    # boundaries_new = unique(mesh_all.facetmarkers)

    # solution = "use smaller refinement or change surface triangulation."
    # length(compartments_new) == ncompartment ||
    #     @error "Incorrect number of compartments, " * solution
    # length(boundaries_new) == nboundary ||
    #     @error "Incorrect number of boundaries, " * solution

    # Deform domain
    if isa(setup, CylinderSetup) && (setup.bend > 1e-16 || setup.twist > 1e-16)
        println("Deforming domain with bend $(setup.bend) and twist $(setup.twist)")
        deform_domain!(mesh_all.points, setup.bend, setup.twist)
    end

    # Split mesh into compartments
    @time mesh = split_mesh(mesh_all)

    mesh, surfaces, cells
end


text:
next_variable_index!(model::ParametricOptimizer) = model.last_variable_index_added += 1
next_parameter_index!(model::ParametricOptimizer) = model.last_parameter_index_added += 1
update_number_of_parameters!(model::ParametricOptimizer) = model.number_of_parameters_in_model += 1

function is_parameter_in_model(model::ParametricOptimizer, v::MOI.VariableIndex)
    return PARAMETER_INDEX_THRESHOLD < v.value <= model.last_parameter_index_added
end

function is_variable_in_model(model::ParametricOptimizer, v::MOI.VariableIndex)
    return 0 < v.value <= model.last_variable_index_added
end

function function_has_parameters(model::ParametricOptimizer, f::MOI.ScalarAffineFunction{T}) where T
    for term in f.terms
        if is_parameter_in_model(model, term.variable_index)
            return true
        end
    end
    return false
end
function function_has_parameters(model::ParametricOptimizer, f::MOI.VectorOfVariables)
    for variable in f.variables
        if is_parameter_in_model(model, variable)
            return true
        end
    end
    return false
end
function function_has_parameters(model::ParametricOptimizer, f::MOI.VectorAffineFunction{T}) where T
    for term in f.terms
        if is_parameter_in_model(model, term.variable_index)
            return true
        end
    end
    return false
end
function function_has_parameters(model::ParametricOptimizer, f::MOI.ScalarQuadraticFunction{T}) where T
    return function_affine_terms_has_parameters(model, f.affine_terms) ||
           function_quadratic_terms_has_parameters(model, f.quadratic_terms)
end
function function_affine_terms_has_parameters(model::ParametricOptimizer, affine_terms::Vector{MOI.ScalarAffineTerm{T}}) where T
    for term in affine_terms
        if is_parameter_in_model(model, term.variable_index)
            return true
        end
    end
    return false
end
function function_quadratic_terms_has_parameters(model::ParametricOptimizer, quadratic_terms::Vector{MOI.ScalarQuadraticTerm{T}}) where T
    for term in quadratic_terms
        if is_parameter_in_model(model, term.variable_index_1) || is_parameter_in_model(model, term.variable_index_2)
            return true
        end
    end
    return false
end

function separate_possible_terms_and_calculate_parameter_constant(model::ParametricOptimizer, terms::Vector{MOI.ScalarAffineTerm{T}}) where T
    vars = MOI.ScalarAffineTerm{T}[]
    params = MOI.ScalarAffineTerm{T}[]
    param_constant = zero(T)

    for term in terms
        if is_variable_in_model(model, term.variable_index)
            push!(vars, term)
        elseif is_parameter_in_model(model, term.variable_index)
            push!(params, term)
            param_constant += term.coefficient * model.parameters[term.variable_index]
        else
            error("Constraint uses a variable that is not in the model")
        end
    end
    return vars, params, param_constant
end
# This version is used on SQFs
function separate_possible_terms_and_calculate_parameter_constant(
                                model::ParametricOptimizer,
                                terms::Vector{MOI.ScalarAffineTerm{T}},
                                variables_associated_to_parameters::Vector{MOI.VariableIndex}
                            ) where T
    vars = MOI.ScalarAffineTerm{T}[]
    params = MOI.ScalarAffineTerm{T}[]
    terms_with_variables_associated_to_parameters = MOI.ScalarAffineTerm{T}[]
    param_constant = zero(T)

    for term in terms
        if is_variable_in_model(model, term.variable_index)
            push!(vars, term)
            if term.variable_index in variables_associated_to_parameters
                push!(terms_with_variables_associated_to_parameters, term)
            end
        elseif is_parameter_in_model(model, term.variable_index)
            push!(params, term)
            param_constant += term.coefficient * model.parameters[term.variable_index]
        else
            error("Constraint uses a variable that is not in the model")
        end
    end
    return vars, params, terms_with_variables_associated_to_parameters, param_constant
end
function separate_possible_terms_and_calculate_parameter_constant(model::ParametricOptimizer, terms::Vector{MOI.ScalarQuadraticTerm{T}}) where T
    quad_params = MOI.ScalarQuadraticTerm{T}[] # parameter x parameter
    quad_aff_vars = MOI.ScalarQuadraticTerm{T}[] # parameter x variable
    quad_vars = MOI.ScalarQuadraticTerm{T}[] # variable x variable
    aff_terms = MOI.ScalarAffineTerm{T}[]
    variables_associated_to_parameters = MOI.VariableIndex[]
    quad_param_constant = zero(T)

    # TODO (guilherme) do we really need the convention ?
    # When we have a parameter x variable or a variable x parameter the convention is to rewrite
    # the SQT with parameter as variable_index_1 and variable as variable_index_2
    for term in terms
        if (is_variable_in_model(model, term.variable_index_1) &&
            is_variable_in_model(model, term.variable_index_2))
            push!(quad_vars, term) # if there are only variables, it remains a quadratic term

        elseif (is_parameter_in_model(model, term.variable_index_1) &&
                is_variable_in_model(model, term.variable_index_2))
            push!(quad_aff_vars, term)
            push!(variables_associated_to_parameters, term.variable_index_2)
            push!(model.multiplicative_parameters, term.variable_index_1.value)
            push!(aff_terms,
                  MOI.ScalarAffineTerm(
                      term.coefficient * model.parameters[term.variable_index_1], term.variable_index_2)
            )

        elseif (is_variable_in_model(model, term.variable_index_1) &&
                is_parameter_in_model(model, term.variable_index_2))
            # Check convention defined above
            push!(quad_aff_vars, MOI.ScalarQuadraticTerm(term.coefficient, term.variable_index_2, term.variable_index_1))
            push!(variables_associated_to_parameters, term.variable_index_2)
            push!(aff_terms,
                  MOI.ScalarAffineTerm(
                      term.coefficient * model.parameters[term.variable_index_2], term.variable_index_1)
            )
            push!(model.multiplicative_parameters, term.variable_index_2.value)

        elseif (is_parameter_in_model(model, term.variable_index_1) &&
                is_parameter_in_model(model, term.variable_index_2))
            push!(quad_params, term)
            push!(model.multiplicative_parameters, term.variable_index_1.value)
            push!(model.multiplicative_parameters, term.variable_index_2.value)
            quad_param_constant += term.coefficient * model.parameters[term.variable_index_1] * model.parameters[term.variable_index_2]
        else
            throw(ErrorException("Constraint uses a variable or parameter that is not in the model"))
        end
    end
    return quad_vars, quad_aff_vars, quad_params, aff_terms, variables_associated_to_parameters, quad_param_constant
end

function fill_quadratic_constraint_caches!(
                                    model::ParametricOptimizer,
                                    new_ci::MOI.ConstraintIndex,
                                    quad_aff_vars::Vector{MOI.ScalarQuadraticTerm{T}},
                                    quad_params::Vector{MOI.ScalarQuadraticTerm{T}},
                                    aff_params::Vector{MOI.ScalarAffineTerm{T}},
                                    terms_with_variables_associated_to_parameters::Vector{MOI.ScalarAffineTerm{T}},
                                ) where T
    if !isempty(quad_aff_vars)
        model.quadratic_constraint_cache_pv[new_ci] = quad_aff_vars
    end
    if !isempty(quad_params)
        model.quadratic_constraint_cache_pp[new_ci] = quad_params
    end
    if !isempty(aff_params)
        model.quadratic_constraint_cache_pc[new_ci] = aff_params
    end
    if !isempty(terms_with_variables_associated_to_parameters)
        model.quadratic_constraint_variables_associated_to_parameters_cache[new_ci] = terms_with_variables_associated_to_parameters
    end
    return nothing
end


text:
using PtFEM

data = Dict(
  # Plane(ndim, nst, nxe, nye, nip, direction, finite_element(nod, nodof), axisymmetric)
  :struc_el => Plane(2, 3, 2, 2, 6, :x, Triangle(3, 2), false),
  :properties => [1.0e6 0.3;],
  :x_coords => [0.0,  0.5,  1.0],
  :y_coords => [0.0,  -0.5,  -1.0],
  :support => [
    (1, [0 1]),
    (4, [0 1]),
    (7, [0 0]),
    (8, [1 0]),
    (9, [1 0])
    ],
  :loaded_nodes => [
    (1, [0.0 -0.25]),
    (2, [0.0 -0.50]),
    (3, [0.0 -0.25])
    ]
)

println()
data |> display
println()

@time m = p51_skyline(data)
println()


text:
module TestLocalCounterExample

using StaticStorages
using Test

const COUNTER_BUCKET = StaticStorages.BucketKey()

struct Counter
    tag::Symbol
    source::LineNumberNode
    _module::Module
    count::Ref{UInt}
end

macro count(tag::Symbol)
    counter = Counter(tag, __source__, __module__, Ref(UInt(0)))
    key = StaticStorages.put!(__module__, COUNTER_BUCKET, counter)
    quote
        $(QuoteNode(counter.count))[] += 1
        $(QuoteNode(key))
    end
end

count_user_a() = @count a
count_user_b() = @count b

counters() = values(StaticStorages.getbucket(COUNTER_BUCKET))

function test_count()
    for c in counters()
        c.count[] = 0
    end

    count_user_a()
    count_user_b()
    count_user_a()
    count_user_a()
    counts = Dict(c.tag => c.count[] for c in counters())
    @test counts[:a] == 3
    @test counts[:b] == 1

    ka = count_user_a()
    ca = StaticStorages.get(COUNTER_BUCKET, ka)
    @test ca.count[] == 4
    @test ca.tag == :a
end

end  # module


text:
#-----------------------------------------------------------------------
#   MINKOWSKI METRIC
#-----------------------------------------------------------------------
η   = x->ημν()

#-----------------------------------------------------------------------
#   INITIAL DATA GENERATION
#-----------------------------------------------------------------------

function idg( ne::Int , Xc::RealVec=rand(4)/100 , spacelike::Bool=true )
    V = rand(4,ne)

    if spacelike
        V[1,:] = -ones(ne)

        for i=1:ne
            V[2:4,i] = V[2:4,i]/norm(V[2:4,i])
        end
        for i=1:ne
            @test ηdot(V[:,i],V[:,i]) ≈ 0 atol=1e-15
        end

        Y   = zeros(4,ne)
        Zi   = zeros(8,ne)

        for i=1:ne
            Y[:,i] = V[:,i] + Xc
        end
        for i=1:ne
        Zi[:,i] = vcat( Y[:,i] , Xc - Y[:,i] )
        end
    else
        V[4,:] = ones(ne)

        for i=1:ne
            V[1,i] = √( 1 + V[2,i]^2 + V[3,i]^2 )
        end
        for i=1:ne
            @test ηdot(V[:,i],V[:,i]) ≈ 0 atol=1e-15
        end

        Y   = zeros(4,ne)
        Zi   = zeros(8,ne)

        for i=1:ne
            Y[:,i] = V[:,i] + Xc
        end
        for i=1:ne
        Zi[:,i] = vcat( Y[:,i] , Xc - Y[:,i] )
        end
    end

    return (Y,Xc,Zi,V)
end