text:
"""
    DecisionVariable <: AbstractVariableRef

Holds a reference to the stochastic program and the corresponding MOI.VariableIndex.
"""
struct DecisionVariable <: JuMP.AbstractVariableRef
    stochasticprogram::StochasticProgram
    index::MOI.VariableIndex
end

# Getters #
# ========================== #
function get_decisions(dvar::DecisionVariable)
    return decisions(structure(owner_model(dvar)))
end
"""
    decision(dvar::DecisionVariable)

Return the internal `Decision` associated with `dvar`.
"""
function decision(dvar::DecisionVariable)
    decisions = get_decisions(dvar)::Decisions
    return decision(decisions, index(dvar))
end
"""
    state(dvar::DecisionVariable)

Return the `DecisionState` of `dvar`.
"""
function state(dvar::DecisionVariable)
    return decision(dvar).state
end

# Setters #
# ========================== #
function take_decisions!(stochasticprogram::StochasticProgram, dvars::Vector{DecisionVariable}, vals::AbstractVector)
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, stochasticprogram), dvars)
    # Check decision length
    length(dvars) == length(vals) || error("Given decision of length $(length(vals)) not compatible with number of decision variables $(length(dvars)).")
    # Update decisions
    for (dvar, val) in zip(dvars, vals)
        d = decision(dvar)
        # Update state
        d.state = Taken
        # Update value
        d.value = val
    end
    # Update objective and constraints in model
    update_decisions!(stochasticprogram, DecisionsStateChange())
    return nothing
end

function untake_decisions!(model::JuMP.Model, dvars::Vector{DecisionVariable})
    # Check that all given decisions are in model
    map(dvar -> check_belongs_to_model(dvar, model), dvars)
    # Update decisions
    need_update = false
    for dvar in dvars
        d = decision(dvar)
        if state(d) == Taken
            need_update |= true
            # Update state
            d.state = NotTaken
        end
    end
    # Update objective and constraints in model (if needed)
    need_update && update_decisions!(stochasticprogram, DecisionsStateChange())
    return nothing
end

function update_decisions!(stochasticprogram::StochasticProgram, change::DecisionModification)
    update_decisions!(structure(stochasticprogram), change)
end

# JuMP overloads #
# ========================== #
function JuMP.delete(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    if stochasticprogram !== owner_model(dvar)
        error("The decision variable you are trying to delete does not " *
              "belong to the stochastic program.")
    end
    MOI.delete(structure(stochasticprogram), index(dvar))
end

function delete(stochasticprogram::StochasticProgram, con_ref::ConstraintRef{<:StochasticProgram})
    if stochasticprogram !== con_ref.model
        error("The constraint reference you are trying to delete does not " *
              "belong to the stochasticprogram.")
    end
    MOI.delete(structure(sto), index(con_ref))
end

# MOI #
# ========================== #
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable)
    check_belongs_to_model(dvar, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        if haskey(cache, :solution)
            # Returned cached solution
            return MOI.get(cache[:solution], attr, index(dvar))
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        return MOI.get(optimizer(stochasticprogram), attr, index(dvar))
    end
    return MOI.get(structure(stochasticprogram), attr, index(dvar))
end
function MOI.get(stochasticprogram::StochasticProgram, attr::MOI.AbstractConstraintAttribute,
                 cr::ConstraintRef{<:StochasticProgram})
    check_belongs_to_model(cr, stochasticprogram)
    if MOI.is_set_by_optimize(attr)
        # Check if there is a cached solution
        cache = solutioncache(stochasticprogram)
        if haskey(cache, :solution)
            # Returned cached solution
            return MOI.get(cache[:solution], attr, index(cr))
        end
        check_provided_optimizer(stochasticprogram.optimizer)
        if MOI.get(stochasticprogram, MOI.TerminationStatus()) == MOI.OPTIMIZE_NOT_CALLED
            throw(OptimizeNotCalled())
        end
        return MOI.get(optimizer(stochasticprogram), attr, index(cr))
    else
        return MOI.get(structure(stochasticprogram), attr, index(cr))
    end
end

function MOI.set(stochasticprogram::StochasticProgram, attr::MOI.AbstractVariableAttribute,
                 dvar::DecisionVariable, value)
    check_belongs_to_model(dvar, stochasticprogram)
    MOI.set(structure(stochasticprogram), attr, index(dvar), value)
    return nothing
end
function MOI.set(stochasticprogram::StochasticProgram, attr::MOI.AbstractConstraintAttribute,
                 cr::ConstraintRef, value)
    check_belongs_to_model(cr, stochasticprogram)
    MOI.set(structure(stochasticprogram), attr, index(cr), value)
end

# JuMP variable interface #
# ========================== #
JuMP.name(dvar::DecisionVariable) = MOI.get(owner_model(dvar), MOI.VariableName(), dvar)::String
function JuMP.set_name(dvar::DecisionVariable, name::String)
    return MOI.set(owner_model(dvar), MOI.VariableName(), dvar, name)
end

function decision_by_name(stochasticprogram::StochasticProgram, name::String)
    index = MOI.get(structure(stochasticprogram), MOI.VariableIndex, name)
    if index isa Nothing
        return nothing
    else
        return DecisionVariable(stochasticprogram, index)
    end
end

JuMP.index(dvar::DecisionVariable) = dvar.index

function JuMP.value(dvar::DecisionVariable; result::Int = 1)::Float64
    d = decision(dvar)
    if d.state == Taken
        # If decision has been fixed the value can be fetched
        # directly
        return d.value
    end
    return MOI.get(owner_model(dvar), MOI.VariablePrimal(result), dvar)
end

function JuMP.is_fixed(dvar::DecisionVariable)
    if state(dvar) == Taken
        return true
    end
    return false
end
"""
    unfix(dvar::DecisionVariable)

Unfix the decision associated with `dvar`. If the decision is already in a `NotTaken` state, this does nothing.

See also [`fix`](@ref).
"""
function JuMP.unfix(dvar::DecisionVariable)
    if state(dvar) == NotTaken
        # Nothing to do, just return
        return nothing
    end
    d = decision(dvar)
    # Update state
    d.state = NotTaken
    # Prepare modification
    change = DecisionStateChange(index(dvar), NotTaken, -value(dvar))
    # Update objective and constraints
    update_decisions!(JuMP.owner_model(dvar), change)
    return nothing
end
"""
    fix(dvar::DecisionVariable, val::Number)

Fix the decision associated with `dvar` to `val`. In contexts where `dvar` is a variable, the variable is fixed to the value. In contexts where `dvar` is a known parameter value, the value is updated.

See also [`unfix`](@ref).
"""
function JuMP.fix(dvar::DecisionVariable, val::Number)
    d = decision(dvar)
    if state(dvar) == NotTaken
        # Prepare modification
        change = DecisionStateChange(index(dvar), Taken, val)
        # Update state
        d.state = Taken
        # Update value
        d.value = val
    else
        # Prepare modification
        change = DecisionStateChange(index(dvar), Taken, val - d.value)
        # Just update value
        d.value = val
    end
    # Update objective and constraints
    update_decisions!(JuMP.owner_model(dvar), change)
    return nothing
end

JuMP.owner_model(dvar::DecisionVariable) = dvar.stochasticprogram

struct DecisionVariableNotOwned <: Exception
    dvar::DecisionVariable
end

function JuMP.check_belongs_to_model(dvar::DecisionVariable, stochasticprogram::StochasticProgram)
    if owner_model(dvar) !== stochasticprogram
        throw(DecisionNotOwned(dvar))
    end
end

function JuMP.is_valid(stochasticprogram::StochasticProgram, dvar::DecisionVariable)
    return stochasticprogram === owner_model(dvar)
end

function JuMP.has_lower_bound(dvar::DecisionVariable)
    index = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), index)
end
function JuMP.LowerBoundRef(dvar::DecisionVariable)
    moi_lb =  MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}
    index = moi_lb(index(dvar).value)
    sp = owner_model(dvar)
    SP = typeof(sp)
    return ConstraintRef{SP, moi_lb, ScalarShape}(sp,
                                                  index,
                                                  ScalarShape())
end
function JuMP.set_lower_bound(dvar::DecisionVariable, lower::Number)
    new_set = MOI.GreaterThan(convert(Float64, lower))
    if has_lower_bound(dvar)
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.GreaterThan{Float64}}(index(dvar).value)
        MOI.set(structure(owner_model(dvar)), MOI.ConstraintSet(), cindex, new_set)
    else
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set)
    end
    return nothing
end
function JuMP.delete_lower_bound(dvar::DecisionVariable)
    JuMP.delete(owner_model(dvar), LowerBoundRef(dvar))
end
function JuMP.lower_bound(dvar::DecisionVariable)
    if !has_lower_bound(dvar)
        error("Decision variable $(dvar) does not have a lower bound.")
    end
    cset = MOI.get(structure(owner_model(dvar)), MOI.ConstraintSet(),
                   LowerBoundRef(dvar))::MOI.GreaterThan{Float64}
    return cset.lower
end

function JuMP.has_upper_bound(dvar::DecisionVariable)
    index = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), index)
end
function JuMP.UpperBoundRef(dvar::DecisionVariable)
    moi_ub =  MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}
    index = moi_ub(index(dvar).value)
    sp = owner_model(dvar)
    SP = typeof(sp)
    return ConstraintRef{SP, moi_ub, ScalarShape}(sp,
                                                  index,
                                                  ScalarShape())
end
function JuMP.set_upper_bound(dvar::DecisionVariable, lower::Number)
    new_set = MOI.LessThan(convert(Float64, lower))
    if has_upper_bound(dvar)
        cindex = MOI.ConstraintIndex{SingleDecision, MOI.LessThan{Float64}}(index(dvar).value)
        MOI.set(structure(owner_model(dvar)), MOI.ConstraintSet(), cindex, new_set)
    else
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), new_set)
    end
    return nothing
end
function JuMP.delete_upper_bound(dvar::DecisionVariable)
    JuMP.delete(owner_model(dvar), UpperBoundRef(dvar))
end
function JuMP.upper_bound(dvar::DecisionVariable)
    if !has_upper_bound(dvar)
        error("Decision $(dvar) does not have a upper bound.")
    end
    cset = MOI.get(structure(owner_model(dvar)), MOI.ConstraintSet(),
                   UpperBoundRef(dvar))::MOI.LessThan{Float64}
    return cset.upper
end

function JuMP.is_integer(dvar::DecisionVariable)
    index = MOI.ConstraintIndex{SingleDecision, MOI.Integer}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), index)
end
function JuMP.set_integer(dvar::DecisionVariable)
    if is_integer(dvar)
        return nothing
    elseif is_binary(dvar)
        error("Cannot set the decision $(dvar) to integer as it is already binary.")
    else
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.Integer())
    end
end
function unset_integer(dvar::DecisionVariable)
    JuMP.delete(owner_model(dvar), IntegerRef(dvar))
    return nothing
end
function JuMP.IntegerRef(dvar::DecisionVariable)
    moi_int =  MOI.ConstraintIndex{SingleDecision, MOI.Integer}
    index = moi_int(index(dvar).value)
    sp = owner_model(dvar)
    SP = typeof(sp)
    return ConstraintRef{SP, moi_int, ScalarShape}(sp,
                                                   index,
                                                   ScalarShape())
end

function JuMP.is_binary(dvar::DecisionVariable)
    index = MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}(index(dvar).value)
    return MOI.is_valid(structure(owner_model(dvar)), index)
end
function JuMP.set_binary(dvar::DecisionVariable)
    if is_binary(dvar)
        return nothing
    elseif is_integer(dvar)
        error("Cannot set the decision $(dvar) to binary as it is already integer.")
    else
        MOI.add_constraint(structure(owner_model(dvar)), SingleDecision(index(dvar)), MOI.ZeroOne())
    end
end
function unset_binary(dvar::DecisionVariable)
    JuMP.delete(owner_model(dvar), BinaryRef(dvar))
    return nothing
end
function JuMP.BinaryRef(dvar::DecisionVariable)
    moi_bin =  MOI.ConstraintIndex{SingleDecision, MOI.ZeroOne}
    index = moi_bin(index(dvar).value)
    sp = owner_model(dvar)
    SP = typeof(sp)
    return ConstraintRef{SP, moi_bin, ScalarShape}(sp,
                                                   index,
                                                   ScalarShape())
end

function JuMP.start_value(dvar::DecisionVariable)
    return MOI.get(owner_model(dvar), MOI.VariablePrimalStart(), dvar)
end
function JuMP.set_start_value(dvar::DecisionVariable, value::Number)
    MOI.set(owner_model(dvar), MOI.VariablePrimalStart(), dvar, Float64(value))
end

function Base.hash(dvar::DecisionVariable, h::UInt)
    return hash(objectid(owner_model(dvar)), hash(dvar.index, h))
end

JuMP.isequal_canonical(d::DecisionVariable, other::DecisionVariable) = isequal(d, other)
function Base.isequal(dvar::DecisionVariable, other::DecisionVariable)
    return owner_model(dvar) === owner_model(other) && dvar.index == other.index
end

Base.iszero(::DecisionVariable) = false
Base.copy(dvar::DecisionVariable) = DecisionVariable(dvar.stochasticprogram, dvar.index)
Base.broadcastable(dvar::DecisionVariable) = Ref(dvar)


text:
module new_logistic_regression

using Distributions

struct data
    y::Vector
    x::Array{Float64}
    m::Vector
end

## Logistic regression nobs
function nobs(data::data)
    nx = size(data.x)[1]
    ny = length(data.y)
    nm = length(data.m)
    if (nx != ny)
        error("number of rows in of x is not equal to the length of y")
    elseif (nx != nm)
        error("number of rows in of x is not equal to the length of m")
    elseif (ny != nm)
        error("length of y is not equal to the length of m")
    end
    return nx
end

function observation(data::data, i::Int64)
    vcat(data.y[i], data.m[i], data.x[i,:])
end

function loglik(theta::AbstractArray,
                obs::AbstractArray)
    y = obs[1]
    m = obs[2]
    x = obs[3:end]
    eta = sum(x .* theta)
    mu = exp.(eta) ./ (1 .+ exp.(eta))
    return y .* log.(mu) + (m - y) .* log.(1 .- mu)
end

function logistic_ef(theta::AbstractArray,
                     obs::AbstractArray)
    y = obs[1]
    m = obs[2]
    x = obs[3:end]
    eta = sum(x .* theta)
    mu = exp.(eta) ./ (1 .+ exp.(eta))
    return x * (y - m * mu)
end

function simulate(theta::Vector,
                  x::Matrix{Float64},
                  m::Vector{Int64})
    n = size(x)[2]
    y = rand.(Binomial.(m, cdf.(Logistic(), x * theta)))
    data(y, x, m);
end

end

ext: jl


text:
"""
    ObjectAgentModule
Module implementing all the operations on agents.
"""
module ObjectAgentModule

include("SquaredMatrixPlusIModule.jl")

using ..ActuatorArrayModule
using ..LinSolvers

using .SquaredMatrixPlusIModule

using LinearAlgebra

export ObjectAgent, updatex!, updatex_dirFixed!, updateu!, broadcastx!,
    actuatorsInCommon, resolveNeighbrRelations!, costFun, costFun!, saturate_control_action

    abstract type ObjectAgent{T<:Real,U<:Unsigned} end

    function resolveNeighbrRelations!(agents::Array{OA, 1}) where {OA<:ObjectAgent}
        Threads.@threads for k in 1:length(agents)
            for l in (k+1):length(agents)
                al = agents[l]
                ak = agents[k]
                il1, il2 = actuatorsInCommon(ak.actList, al.actList)

                push!(ak.neighbors, (al, il1, il2));
                push!(al.neighbors, (ak, il2, il1));
            end
        end
    end

    function updatex!(oa::ObjectAgent{T,U}, λ::T) where {T<:Real, U<:Unsigned}
        oa.xk .= oa.zk .- oa.uk;

        # Compute the jacobian
        jac!(oa);

        costFun!(oa); # update fvk,  < 20us

        #lsd.xsol sits in the same spot in the memory as oa.decDit (they are the same)
        kernel_solve!(oa.lsd, oa.J, λ, oa.fvk);

        for i in 1:length(oa.act_used)
            oa.xk[oa.act_used[i]] -= oa.decDir[i];
        end

        # return norm(oa.fvk);
        nothing;
    end

    function updatex_dirFixed!(oa::ObjectAgent{T,U}, λ::T) where {T<:Real, U<:Unsigned}
        error(":fixdir has not been supported yet")

        n_red = length(act_used);

        # Compute the jacobian
        jac!(oa); # ACU: ~10 us, DEP ~40 us

        c = 1e-3;
        A_kkt = vcat([(1/λ + c)*I zeros(n_red,1) -oa.J],
                     [zeros(1,n_red) 1 [oa.Fdes[1] oa.Fdes[2]]],
                     [-oa.J' [oa.Fdes[1] oa.Fdes[2]]' zeros(2,2)]);
        kktb = vcat(zeros(T, n_red, 1), [1.0; 0.0; 0.0]);

        # return norm(oa.fvk);
        nothing;
    end

    function saturate_control_action(oa::ObjectAgent{T,U}, x::T) where {T<:Real, U<:Unsigned}
        """
        This method is used to saturate the control aciton. This default implementation
        assumes that there are no saturations.
        """
        x
    end

    function broadcastx!(oa::ObjectAgent)
        for (neighbor, il1, il2) in oa.neighbors
            receivex!(neighbor, oa.xk, oa.uk, il1, il2);
        end
    end

    function receivex!(oa::ObjectAgent{T,U}, xk_nghbr::AbstractVector{T}, uk_nghbr::AbstractVector{T}, il1::Array{U, 1}, il2::Array{U, 1}) where {T<:Real,U<:Unsigned}
        for k in 1:length(il1)
            oa.rcvBuffer_x[il2[k]] += xk_nghbr[il1[k]];
            oa.rcvBuffer_u[il2[k]] += uk_nghbr[il1[k]];
            oa.rcvBuffer_N[il2[k]] += 1;
        end
    end

    function updateu!(oa::ObjectAgent, ρ::Real)
        for i in 1:length(oa.xk)
            xbari = (oa.xk[i] + oa.rcvBuffer_x[i])/(oa.rcvBuffer_N[i]+1);
            ubari = (oa.uk[i] + oa.rcvBuffer_u[i])/(oa.rcvBuffer_N[i]+1);
            oa.zk[i] = saturate_control_action(oa, xbari + ubari);
            oa.uk[i] += ρ*(oa.xk[i] - oa.zk[i]);

            oa.rcvBuffer_x[i] = 0;
            oa.rcvBuffer_u[i] = 0;
            oa.rcvBuffer_N[i] = 0;
        end
    end

    function Base.show(io::IO, oa::ObjectAgent)
        xp, yp = oa.pos;
        println(io, "--Agent: $(oa.name) --")
        println(io, "Position: ($xp, $yp)")

        println(io, "List of actuators:")
        for (i, j) in oa.actList
            print(io, "($i, $j), ")
        end
    end

    include("MagMan.jl")
    include("AcouMan.jl")
    include("DEPman.jl")
end


text:
using JDF
using Test
using DataFrames
using Random: randstring
using WeakRefStrings

include("test-categorical-ararys.jl")
include("test-stringarray.jl")
include("test-symbol.jl")
include("test-jdf-str.jl")
include("test-eachcol.jl")
include("test-ZonedDateTime.jl")
include("test-substring.jl")

@testset "JDF.jl parallel" begin
    df = DataFrame([collect(1:100) for i = 1:3000])
    df[!, :int_missing] =
        rand([rand(rand([UInt, Int, Float64, Float32, Bool])), missing], nrow(df))

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i = 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i = 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i = 1:nrow(df)]
    df[!, :stringarray_missing] =
        StringVector([rand([missing, randstring(8)]) for i = 1:nrow(df)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i = 1:nrow(df)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:nrow(df), 10), :char_missing] .= missing

    @time JDF.save("a.jdf", df)
    @time df2 = DataFrame(JDF.load("a.jdf"), copycols=false)

    isequal(df, df2)

    df2 = DataFrame(JDF.load("a.jdf"), copycols=false)
    @test ncol(df2) == 3009
    @test nrow(df2) == 100

    df2[!, :stringarray_missing]

    @test all(all(isequal(df[!, n], df2[!, n])) for n in names(df))

    # clean up
    rm("a.jdf", force = true, recursive = true)
end

@testset "JDF.jl serial" begin
    df = DataFrame([collect(1:100) for i = 1:3000])
    df[!, :int_missing] =
        rand([rand(rand([UInt, Int, Float64, Float32, Bool])), missing], nrow(df))

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i = 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i = 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i = 1:nrow(df)]
    df[!, :stringarray_missing] =
        StringVector([rand([missing, randstring(8)]) for i = 1:nrow(df)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i = 1:nrow(df)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:nrow(df), 10), :char_missing] .= missing

    ssavejdf("a.jdf", df)

    df2 = DataFrame(sloadjdf("a.jdf"), copycols=false)
    @test ncol(df2) == 3009
    @test nrow(df2) == 100

    @test all(all(isequal(df[!, n], df2[!, n])) for n in names(df))

    # clean up
    rm("a.jdf", force = true, recursive = true)
end


text:
function update_A!(
    A::AbstractArray{T,3} where {T},
    θ_Q::AbstractArray{T,3} where {T},
    ξ::AbstractArray,
    α::AbstractMatrix,
    β::AbstractMatrix,
    LL::AbstractMatrix,
    n2t::AbstractArray{Int},
    model_A::Model;
    warm_start = true
)
    @argcheck size(α, 1) == size(β, 1) == size(LL, 1) == size(ξ, 1)
    @argcheck size(α, 2) ==
              size(β, 2) ==
              size(LL, 2) ==
              size(A, 1) ==
              size(A, 2) ==
              size(ξ, 2) ==
              size(ξ, 3)

    N, K = size(LL)
    T = size(A, 3)
    @inbounds for n in OneTo(N - 1)
        t = n2t[n] # periodic t
        m = SHHMM.vec_maximum(view(LL, n + 1, :)) #TODO Should be imported from HMMBase maxmouchet
        c = 0.0

        for i in OneTo(K), j in OneTo(K)
            ξ[n, i, j] = α[n, i] * A[i, j, t] * exp(LL[n+1, j] - m) * β[n+1, j]
            c += ξ[n, i, j]
        end

        for i in OneTo(K), j in OneTo(K)
            ξ[n, i, j] /= c
        end
    end
    # ξ are the filtering probablies
    pklj_jump = model_A[:pklj_jump]
    πkl = model_A[:πkl]
    πk = model_A[:πk]

    for k = 1:K
        ## Update the smoothing parameters in the JuMP model
        [set_value(πkl[n, l], ξ[n, k, l]) for n = 1:N-1, l = 1:K-1]
        [set_value(πk[n], sum(ξ[n, k, l] for l = 1:K)) for n = 1:N-1]
        warm_start && set_start_value.(pklj_jump, θ_Q[k, :, :])
        # Optimize the updated model
        optimize!(model_A)
        # Obtained the new parameters
        θ_Q[k, :, :] = value.(pklj_jump)
    end

    [A[k, l, t] = exp(polynomial_trigo(t, θ_Q[k, l, :], T = T)) for k = 1:K, l = 1:K-1, t = 1:T]
    [A[k, K, t] = 1 for k = 1:K, t = 1:T] # last colum is 1/normalization (one could do otherwise)
    normalization_polynomial = [1 + sum(A[k, l, t] for l = 1:K-1) for k = 1:K, t = 1:T]
    [A[k, l, t] /= normalization_polynomial[k, t] for k = 1:K, l = 1:K, t = 1:T]
end

function fit_mle_B_trig1!(θ_Y, γ, k, s, h, model_B; warm_start = true)
    θ_jump = model_B[:θ_jump]
    mle = model_B[:mle][k, s, h]
    warm_start && set_start_value.(θ_jump, θ_Y[k, s, h, :])
    # set_silent(model_B)
    @NLobjective(
        model_B, Max,
        mle
    )
    optimize!(model_B)
    θ_Y[k, s, h, :] = value.(θ_jump)
end

function fit_mle_B_trig01!(θ_Y, model_B, mle; warm_start = true)
    θ_jump = model_B[:θ_jump]
    warm_start && set_start_value.(θ_jump, θ_Y[:])
    # set_silent(model_B)
    @NLobjective(
        model_B, Max,
        mle
    )
    optimize!(model_B)
    θ_Y[:] = value.(θ_jump)
end

function fit_mle_B_trig!(θ_Y::AbstractVector, model::Model, γ::AbstractVector; warm_start = true)
    N = size(γ, 1)
    θ_jump = model[:θ_jump]
    πk = model[:πk]
    ## Update the smoothing parameters in the JuMP model
    for n = 1:N
        set_value(πk[n], γ[n])
    end
    warm_start && set_start_value.(θ_jump, θ_Y[:])
    #set_silent(model)

    optimize!(model)
    θ_Y[:] = value.(θ_jump)
end

function update_B!(B::AbstractArray{T,4} where {T}, θ_Y::AbstractArray{N,4} where {N}, γ::AbstractMatrix, observations, model_B::Model, mles; warm_start = true)
    @argcheck size(γ, 1) == size(observations, 1)
    @argcheck size(γ, 2) == size(B, 1)
    N = size(γ, 1)
    K = size(B, 1)
    T = size(B, 2)
    D = size(B, 3)
    size_memory = size(B, 4)
    ## For periodicHMM only the n observations corresponding to B(t) are used to update B(t)
    θ_jump = model_B[:θ_jump]
    πk = model_B[:πk]
    ## Update the smoothing parameters in the JuMP model
    for n = 1:N, k = 1:K
        set_value(πk[n, k], γ[n, k])
    end
    for k = 1:K, s = 1:D, h = 1:size_memory
        warm_start && set_start_value.(θ_jump, θ_Y[k, s, h, :])
        # set_silent(model_B)
        @NLobjective(
            model_B, Max,
            mles[k, s, h]
        )
        optimize!(model_B)
        θ_Y[k, s, h, :] = value.(θ_jump)
    end
    p = [1 / (1 + exp(polynomial_trigo(t, θ_Y[k, s, h, :], T = T))) for k = 1:K, t = 1:T, s = 1:D, h = 1:size_memory]
    B[:, :, :, :] = Bernoulli.(p)
end

# JuMP model use to increase R(θ,θ^i) for the Q(t) matrix
function model_for_A(ξ::AbstractArray, n2t::AbstractArray{Int}, d::Int, T::Int; silence = true)
    K = size(ξ, 2)
    N = length(n2t)
    model = Model(Ipopt.Optimizer)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * n2t[n]) for n = 1:(N-1), j = 1:d]
    sin_nj = [sin(f * j * n2t[n]) for n = 1:(N-1), j = 1:d]

    trig_nj = [[1; interleave2(cos_nj[n, :], sin_nj[n, :])] for n = 1:N-1]

    @variable(model, pklj_jump[l = 1:(K-1), j = 1:(2d+1)], start = 0.01)
    # Polynomial P_kl

    @NLexpression(model, Pkl[n = 1:N-1, l = 1:K-1], sum(trig_nj[n][j] * pklj_jump[l, j] for j = 1:length(trig_nj[n])))

    @NLparameter(model, πkl[n = 1:N-1, l = 1:K-1] == ξ[n, l])
    @NLparameter(model, πk[n = 1:N-1] == sum(ξ[n, l] for l = 1:K))

    @NLobjective(
        model,
        Max,
        sum(sum(πkl[n, l] * Pkl[n, l] for l = 1:K-1) - πk[n] * log1p(sum(exp(Pkl[n, l]) for l = 1:K-1)) for n = 1:N-1)
    )
    # To add NL parameters to the model for later use https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:πkl] = πkl
    model[:πk] = πk
    return model
end

function model_for_B(γ::AbstractMatrix, n2t, n_occurence_history, d::Int, T::Int, D::Int, size_memory::Int; silence = true)
    K = size(γ, 2)
    N = length(n2t)
    model = Model(Ipopt.Optimizer)
    set_optimizer_attribute(model, "max_cpu_time", 60.0)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * n2t[n]) for n = 1:N, j = 1:d]
    sin_nj = [sin(f * j * n2t[n]) for n = 1:N, j = 1:d]

    trig = [[1; interleave2(cos_nj[n, :], sin_nj[n, :])] for n = 1:N]
    @variable(model, θ_jump[j = 1:(2d+1)])
    # Polynomial P
    @NLexpression(model, Pn[n = 1:N], sum(trig[n][j] * θ_jump[j] for j = 1:length(trig[n])))

    @NLparameter(model, πk[n = 1:N, k = 1:K] == γ[n, k])

    @NLexpression(model, mle[k = 1:K, s = 1:D, h = 1:size_memory],
        -sum(πk[n, k] * log1p(exp(-Pn[n])) for n in n_occurence_history[s, h, 1]) - sum(πk[n, k] * log1p(exp(Pn[n])) for n in n_occurence_history[s, h, 2])
    ) # 1 is where it did not rain # 2 where it rained

    # I don't know if it is the best but https://discourse.julialang.org/t/jump-updating-nlparameter-of-a-model-in-a-loop/35081/3
    model[:πk] = πk
    model[:mle] = mle
    return model
end

function fit_mle2(hmm::HierarchicalPeriodicHMM, observations, n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}; all_iters = false, kwargs...)
    hmm = copy(hmm)
    θ_Q = copy(θ_Q)
    θ_Y = copy(θ_Y)

    history, all_θ_Q, all_θ_Y = fit_mle2!(hmm, observations, n2t, θ_Q, θ_Y; kwargs...)
    if all_iters == true
        return hmm, θ_Q, θ_Y, history, all_θ_Q, all_θ_Y
    else
        return hmm, θ_Q, θ_Y, history
    end
end

function fit_mle3(hmm::HierarchicalPeriodicHMM, observations, n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}; all_iters = false, kwargs...)
    hmm = copy(hmm)
    θ_Q = copy(θ_Q)
    θ_Y = copy(θ_Y)

    history, all_θ_Q, all_θ_Y = fit_mle3!(hmm, observations, n2t, θ_Q, θ_Y; kwargs...)
    if all_iters == true
        return hmm, θ_Q, θ_Y, history, all_θ_Q, all_θ_Y
    else
        return hmm, θ_Q, θ_Y, history
    end
end

function fit_mle(hmm::HierarchicalPeriodicHMM, observations, n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}; all_iters = false, kwargs...)
    hmm = copy(hmm)
    θ_Q = copy(θ_Q)
    θ_Y = copy(θ_Y)

    history, all_θ_Q, all_θ_Y = fit_mle!(hmm, observations, n2t, θ_Q, θ_Y; kwargs...)
    if all_iters == true
        return hmm, θ_Q, θ_Y, history, all_θ_Q, all_θ_Y
    else
        return hmm, θ_Q, θ_Y, history
    end
end

function fit_mle1(hmm::HierarchicalPeriodicHMM, observations, n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}; all_iters = false, kwargs...)
    hmm = copy(hmm)
    θ_Q = copy(θ_Q)
    θ_Y = copy(θ_Y)

    history, all_θ_Q, all_θ_Y = fit_mle1!(hmm, observations, n2t, θ_Q, θ_Y; kwargs...)
    if all_iters == true
        return hmm, θ_Q, θ_Y, history, all_θ_Q, all_θ_Y
    else
        return hmm, θ_Q, θ_Y, history
    end
end

function fit_mle2!(
    hmm::HierarchicalPeriodicHMM,
    observations,
    n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}
    ;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    silence = true,
    warm_start = true
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T, size_memory, D = size(observations, 1), size(hmm, 1), size(hmm, 3), size(hmm, 4), size(hmm, 2)

    deg_Q = (size(θ_Q, 3) - 1) ÷ 2
    deg_Y = (size(θ_Y, 4) - 1) ÷ 2

    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    all_θ_Q = [copy(θ_Q)]
    all_θ_Y = [copy(θ_Y)]
    # Allocate memory for in-place updates
    c = zeros(N)
    α = zeros(N, K)
    β = zeros(N, K)
    γ = zeros(N, K)
    ξ = zeros(N, K, K)
    LL = zeros(N, K)

    # assign category for observation depending in the past observations
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(observations, memory)
    n_occurence_history = [findall(.&(observations[:, j] .== y, lag_cat[:, j] .== h)) for j = 1:D, h = 1:size_memory, y = 0:1]

    model_A = model_for_A(ξ[:, 1, :], n2t, deg_Q, T, silence = silence) # JuMP Model for transition matrix
    model_B = [model_for_B2(γ[:, 1], n2t, n_occurence_history[j, h, :], deg_Y, T, D, size_memory, silence = silence) for j = 1:D, h = 1:size_memory] # JuMP Model for transition matrix

    loglikelihoods!(LL, hmm, observations, n2t, lag_cat)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
    backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, θ_Q, ξ, α, β, LL, n2t, model_A; warm_start = warm_start)
        update_B2!(hmm.B, θ_Y, γ, observations, model_B; warm_start = warm_start)
        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        push!(all_θ_Q, copy(θ_Q))
        push!(all_θ_Y, copy(θ_Y))

        # loglikelihoods!(LL, hmm, observations, n2t)
        loglikelihoods!(LL, hmm, observations, n2t, lag_cat)

        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
        backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history, all_θ_Q, all_θ_Y
end

function fit_mle3!(
    hmm::HierarchicalPeriodicHMM,
    observations,
    n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}
    ;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    silence = true,
    warm_start = true
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T, size_memory, D = size(observations, 1), size(hmm, 1), size(hmm, 3), size(hmm, 4), size(hmm, 2)

    deg_Q = (size(θ_Q, 3) - 1) ÷ 2
    deg_Y = (size(θ_Y, 4) - 1) ÷ 2

    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    all_θ_Q = [copy(θ_Q)]
    all_θ_Y = [copy(θ_Y)]
    # Allocate memory for in-place updates
    c = zeros(N)
    α = zeros(N, K)
    β = zeros(N, K)
    γ = zeros(N, K)
    ξ = zeros(N, K, K)
    LL = zeros(N, K)

    # assign category for observation depending in the past observations
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(observations, memory)
    n_occurence_history = [findall(.&(observations[:, j] .== y, lag_cat[:, j] .== h)) for j = 1:D, h = 1:size_memory, y = 0:1]

    model_A = model_for_A(ξ[:, 1, :], n2t, deg_Q, T, silence = silence) # JuMP Model for transition matrix
    model_B = [model_for_B3(γ[:, 1], n2t, n_occurence_history[j, h, :], deg_Y, T, D, size_memory, silence = silence) for j = 1:D, h = 1:size_memory] # JuMP Model for transition matrix

    loglikelihoods!(LL, hmm, observations, n2t, lag_cat)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
    backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, θ_Q, ξ, α, β, LL, n2t, model_A; warm_start = warm_start)
        update_B2!(hmm.B, θ_Y, γ, observations, model_B; warm_start = warm_start)
        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        push!(all_θ_Q, copy(θ_Q))
        push!(all_θ_Y, copy(θ_Y))

        # loglikelihoods!(LL, hmm, observations, n2t)
        loglikelihoods!(LL, hmm, observations, n2t, lag_cat)

        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
        backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history, all_θ_Q, all_θ_Y
end
function fit_mle!(
    hmm::HierarchicalPeriodicHMM,
    observations,
    n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}
    ;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    silence = true,
    warm_start = true
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T, size_memory, D = size(observations, 1), size(hmm, 1), size(hmm, 3), size(hmm, 4), size(hmm, 2)

    deg_Q = (size(θ_Q, 3) - 1) ÷ 2
    deg_Y = (size(θ_Y, 4) - 1) ÷ 2

    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    all_θ_Q = [copy(θ_Q)]
    all_θ_Y = [copy(θ_Y)]
    # Allocate memory for in-place updates
    c = zeros(N)
    α = zeros(N, K)
    β = zeros(N, K)
    γ = zeros(N, K)
    ξ = zeros(N, K, K)
    LL = zeros(N, K)

    # assign category for observation depending in the past observations
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(observations, memory)
    n_occurence_history = [findall(.&(observations[:, j] .== y, lag_cat[:, j] .== h)) for j = 1:D, h = 1:size_memory, y = 0:1]

    model_A = model_for_A(ξ[:, 1, :], n2t, deg_Q, T, silence = silence) # JuMP Model for transition matrix
    model_B = model_for_B(γ, n2t, n_occurence_history, deg_Y, T, D, size_memory, silence = silence) # JuMP Model for transition matrix
    mles = model_B[:mle]

    loglikelihoods!(LL, hmm, observations, n2t, lag_cat)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
    backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, θ_Q, ξ, α, β, LL, n2t, model_A; warm_start = warm_start)
        update_B!(hmm.B, θ_Y, γ, observations, model_B, mles; warm_start = warm_start)
        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        push!(all_θ_Q, copy(θ_Q))
        push!(all_θ_Y, copy(θ_Y))

        # loglikelihoods!(LL, hmm, observations, n2t)
        loglikelihoods!(LL, hmm, observations, n2t, lag_cat)

        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
        backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history, all_θ_Q, all_θ_Y
end

function fit_mle1!(
    hmm::HierarchicalPeriodicHMM,
    observations,
    n2t::AbstractArray{Int},
    θ_Q::AbstractArray{TQ,3} where {TQ},
    θ_Y::AbstractArray{TY,4} where {TY}
    ;
    display = :none,
    maxiter = 100,
    tol = 1e-3,
    robust = false,
    silence = true,
    warm_start = true
)
    @argcheck display in [:none, :iter, :final]
    @argcheck maxiter >= 0

    N, K, T, size_memory, D = size(observations, 1), size(hmm, 1), size(hmm, 3), size(hmm, 4), size(hmm, 2)

    deg_Q = (size(θ_Q, 3) - 1) ÷ 2
    deg_Y = (size(θ_Y, 4) - 1) ÷ 2

    @argcheck T == size(hmm.B, 2)
    history = EMHistory(false, 0, [])

    all_θ_Q = [copy(θ_Q)]
    all_θ_Y = [copy(θ_Y)]
    # Allocate memory for in-place updates
    c = zeros(N)
    α = zeros(N, K)
    β = zeros(N, K)
    γ = zeros(N, K)
    ξ = zeros(N, K, K)
    LL = zeros(N, K)

    # assign category for observation depending in the past observations
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(observations, memory)
    n_occurence_history = [findall(.&(observations[:, j] .== y, lag_cat[:, j] .== h)) for j = 1:D, h = 1:size_memory, y = 0:1]

    model_A = model_for_A(ξ[:, 1, :], n2t, deg_Q, T, silence = silence) # JuMP Model for transition matrix
    model_B = model_for_B(γ, n2t, n_occurence_history, deg_Y, T, D, size_memory, silence = silence) # JuMP Model for transition matrix

    loglikelihoods!(LL, hmm, observations, n2t, lag_cat)
    robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

    forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
    backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
    posteriors!(γ, α, β)

    logtot = sum(c)
    (display == :iter) && println("Iteration 0: logtot = $logtot")

    for it = 1:maxiter
        update_a!(hmm.a, α, β)
        update_A!(hmm.A, θ_Q, ξ, α, β, LL, n2t, model_A; warm_start = warm_start)
        update_B1!(hmm.B, θ_Y, γ, observations, model_B; warm_start = warm_start)
        # Ensure the "connected-ness" of the states,
        # this prevents case where there is no transitions
        # between two extremely likely observations.
        robust && (hmm.A .+= eps())

        @check isprobvec(hmm.a)
        @check istransmats(hmm.A)

        push!(all_θ_Q, copy(θ_Q))
        push!(all_θ_Y, copy(θ_Y))

        # loglikelihoods!(LL, hmm, observations, n2t)
        loglikelihoods!(LL, hmm, observations, n2t, lag_cat)

        robust && replace!(LL, -Inf => nextfloat(-Inf), Inf => log(prevfloat(Inf)))

        forwardlog!(α, c, hmm.a, hmm.A, LL, n2t)
        backwardlog!(β, c, hmm.a, hmm.A, LL, n2t)
        posteriors!(γ, α, β)

        logtotp = sum(c)
        (display == :iter) && println("Iteration $it: logtot = $logtotp")

        push!(history.logtots, logtotp)
        history.iterations += 1

        if abs(logtotp - logtot) < tol
            (display in [:iter, :final]) &&
                println("EM converged in $it iterations, logtot = $logtotp")
            history.converged = true
            break
        end

        logtot = logtotp
    end

    if !history.converged
        if display in [:iter, :final]
            println("EM has not converged after $(history.iterations) iterations, logtot = $logtot")
        end
    end

    history, all_θ_Q, all_θ_Y
end


text:
### A Pluto.jl notebook ###
# v0.14.7

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 37080ff2-a3b7-11eb-3634-af524d492479
using Plots, SSpline, SpectralVis

# ╔═╡ accd501b-c865-485b-9c8c-4d7b074b22a3
using PlutoUI

# ╔═╡ 201148e9-e528-435e-afc8-9e4390e2729d
#plotlyjs()
plotly()

# ╔═╡ 22394b49-4e75-444e-b19f-1bfc637dbe98
md"Cell width:"

# ╔═╡ e62085f0-b38a-43cc-a277-0504da926958
html"""<style>main {max-width: 950px;}</style>"""

# ╔═╡ 9e05057a-216a-44eb-95cc-b15b46fc2d3d
md"Color matching function defined in the spectral environment SPECENV:"

# ╔═╡ 54e4792f-0294-4c58-931b-4c3a5ac624ed
begin
	#SPECENV.cmf= SpectralVis.colormatch(CMF2012_2)
	cmfunc = SPECENV.cmf
end

# ╔═╡ 28c0b9e7-9faa-4814-8e1d-5c67a58dfe14
md"Cubic spline interpolation and extrapolations to mach the CMF range:"

# ╔═╡ 54dfc832-dcc2-463c-9be9-ba785f0fec30
md"colorpatch:"

# ╔═╡ 455c4526-1dd0-461d-a40d-d4c9be7b9c2c
md"Definitions of source and destination white points and illuminant spectra:"

# ╔═╡ d96a11d2-6321-4b88-9798-9809533a734f
md"Definitions of Macbeth Chart ground truth, source, and destination colors according to user choice of CCT and CAT:"

# ╔═╡ 07df29f9-b7aa-4877-bda6-8e6c32c6d891
md"MER Caltarget patches"

# ╔═╡ 27c3f7c1-afc6-4095-86a2-85c4651ab75d
md"Color Checker patch spectrum:"

# ╔═╡ d128b785-6a9b-48c0-973a-4acfd24fa97c
md"""## Macbeth Color Checker
All colors are shown in sRGB color space. **Warning**: color patch **Cyan** (3rd row, rightmost column) is just outside of sRGB color gamut. Even ground truth color in the outer ring will display a slight error."""

# ╔═╡ bab11da2-681f-41a8-ae66-c73b6acaf38b
md"""* **Outer ring**: spectral ground truth color at whitepoint D65 (sRGB white point)

* **Inner circle**: color transformed from source CCT to D65 using selected CAT
"""

# ╔═╡ 04857ffb-aa2e-4dd8-af49-10bcc22c4356
charttype=:Macbeth_chart2

# ╔═╡ b69148f5-789e-4e95-8f3b-c78ed71de534
md"Color Adaptation Transform:"

# ╔═╡ 798d5450-00fd-4f18-a962-56374e3e0fcc
@bind CAT Select(["xyzscaling" => "XYZ scaling", "vonkries" => "vonKries scaling", "sharp" =>"Sharp transform", "bradford" => "Bradford transform", "fairchild" => "Fairchild transform", "cat02" => "CAT02 transform", "hpe" => "Hunt-Pointer-Estevez", "bestof20" => "Süsstrunk"], default="vonkries")

# ╔═╡ e59483cc-5366-40cb-8086-381915eb3531
@bind sourcetemp Select(["500" => "*D5", "600" => "*D6", "700" => "*D7", "750" => "*D7.5", "800" => "*D8", "900" => "*D9", "1000" => "*D10", "1500" => "*D15", "2000" => "*D20", "2200" => "*D22", "2700" => "*D27", "3400" => "*D34", "4000" => "D40", "5000" => "D50", "5500" => "D55", "6500" => "D65", "7500" => "D75", "9300" => "D93", "10000" => "D100", "15000" => "D150", "17500" => "D175", "20000" => "D200", "22500" => "D225", "25000" => "D250", "30000" => "*D300", "40000" => "*D400", "50000" => "*D500", "75000" => "*D750"], default = "4000")

#@bind sourcetemp Select(["1000"=> 1000,"1500"=> 1500,"1750"=>1750,"2000"=>2000,"2250"=>2250,"2500"=>2500,"2700"=>2700,"3000"=>3000,"3400"=>3400,"3750"=>3750,"4000"=>4000,"4500"=>4500,"5000"=>5000,"5500"=>5500,"6000"=>6000,"6500"=>6500,"7000"=>7000,"9300"=>9300,"10000"=>10000,"15000"=>15000,"20000"=>20000,"25000"=>25000], default = "4000")

# ╔═╡ c71b58d4-a969-4153-ade5-3ee27044b482
begin
	stemp = parse(Int,sourcetemp)
	cels = round(stemp-273.15,digits=2)
	adapt = Symbol(CAT)
	sWP = D_series_whitepoint(SPECENV,stemp)
	dWP = D_series_whitepoint(SPECENV,6504.0)
	source = normalized_D_series(SPECENV,stemp)
	dest = normalized_D_series(SPECENV,6504.0)
	#sWP = blackbody_whitepoint(SPECENV,stemp)
	#dWP = blackbody_whitepoint(SPECENV,6500)
	#source = normalize_blackbody(blackbody_illuminant(SPECENV,stemp))
	#dest = normalize_blackbody(blackbody_illuminant(SPECENV,6500))
	#source = normalize_spec(blackbody_illuminant(SPECENV,stemp),1.0)
	#dest = normalize_spec(blackbody_illuminant(SPECENV,6500),1.0)
end

# ╔═╡ e94743af-b474-4801-ac3e-53703c4ab56b
begin
	groundtruth_patches=Array{Any}(undef,24)
	source_patches=Array{Any}(undef,24)
	dest_patches=Array{Any}(undef,24)
	@inbounds for i in 1:24
		macbeth_colorpatch= macbeth(charttype,i)
		cpspline2=cubicspline(macbeth_colorpatch.λ,macbeth_colorpatch.r,:natural)
		cpextr2=extrap(cpspline2,SPECENV.λmin:SPECENV.λmax,:linear)
    	λ2,s2 = interp(cpextr2,(390.0:1.0:830.0))
    	colorpatch2=RSpec(λ2,s2)
		groundtruth_patches[i]= ×(dest, colorpatch2, cmfunc)
		source_patches[i]= ×(source, colorpatch2, cmfunc)
		dest_patches[i]= chromatic_adaptation(×(source, colorpatch2, cmfunc), sWP, dWP, Symbol(CAT))
	end
end

# ╔═╡ ad9c6ea9-f6c3-415e-bc65-ca817809db8b
begin
	merpatches = Array{Any}(undef,7)
	mersource = Array{Any}(undef,7)
	merdest = Array{Any}(undef,7)

	@inbounds for i in 1:7
		merpatch=mer_caltarget(i)
		merspline=spline3(merpatch.λ,merpatch.r,:natural)
		merextr=extrap(merspline,SPECENV.λmin:SPECENV.λmax,:linear)
		λmer,rmer = interp(merextr,(SPECENV.λmin:SPECENV.Δλ:SPECENV.λmax))
		merspec= reflectance_spec(λmer,rmer)
		merpatches[i] = ×(dest, merspec, cmfunc)
		mersource[i] = ×(source, merspec, cmfunc)
		merdest[i] = chromatic_adaptation(×(source, merspec, cmfunc), sWP, dWP, Symbol(CAT))
	end
end

# ╔═╡ e99970e1-5b00-42d0-a464-51cded141acf
md"Source CCT: $stemp K / $cels °C"

# ╔═╡ 261ebe67-9f43-455e-ab37-51f254d26e67
begin
	xcoord=zeros(24)
	ycoord=zeros(24)

	plot(xrange=(0.2,6.8),yrange=(-0.4,3.6),background=Lab(30,0,0),foreground=Lab(26,0,0),ticks=false,legend=false,size=(820,560))
	for i=1:24
		row,col=divrem(i-1,6).+(1,1)
    	xcoord[i]=col
		ycoord[i]=4.1-row
	end
	scatter!(xcoord,ycoord,color=groundtruth_patches,markersize=52,markerstrokewidth=0,markershape= :circle, aspect_ratio=:equal)
	scatter!(xcoord,ycoord,color=dest_patches,markersize=30,markerstrokewidth=0,markershape= :circle, aspect_ratio=:equal)
end

# ╔═╡ fa59297b-e992-4285-8295-b5c8273d8a9d
md"Patch Array initialization"

# ╔═╡ c3370bb8-b574-4c82-a1a2-5cbe04e3da29
md"Color shift from pre-adaptation to adaptation with applied CAT vs. ground truth"

# ╔═╡ 5e8c03cc-c67a-445e-b775-b339b0beba03
md"Source WP | destination WP | source color | destin. color | ground truth color:"

# ╔═╡ e22a766c-1ab4-4f6d-9b38-98bc579f6de9
md"rel. to D65 | D65 whitepoint | relat. to D65 | adapted |"

# ╔═╡ 57ac1b78-5d10-4ff4-b30b-430cc111f925
md"Patch color:"

# ╔═╡ 6eafc330-c225-4f26-8a17-561467a8be13
begin
	colordict=["1" => "dark skin", "2" => "light skin", "3" => "blue sky", "4" => "foliage", "5" => "blue flower", "6" => "bluish green", "7" => "orange", "8" => "purplish blue", "9" => "moderate red", "10" => "purple", "11" => "yellow green", "12" => "orange yellow", "13" => "blue", "14" => "green", "15" => "red", "16" => "yellow", "17" => "magenta", "18" => "cyan", "19" => "white 9.5", "20" => "neutral 8", "21" => "neutral 6.5", "22" => "neutral 5", "23" => "neutral 3.5", "24" => "black 2"]
@bind patch Select(colordict, default = "1")
end

# ╔═╡ f9c57d5c-b947-445f-b479-91515cd473b8
mbcolorpatch= macbeth(:Macbeth_chart2,parse(Int,patch))

# ╔═╡ e9951abb-915a-4799-8ed7-5c97490e8c8e
begin
	cpspline=cubicspline(mbcolorpatch.λ,mbcolorpatch.r,:natural)
	cpextr=extrap(cpspline,390.0:830.0,:linear)
    l,s = interp(cpextr,(390.0:1.0:830.0))
    colorpatch=RSpec(l,s)
end

# ╔═╡ 4d5584ab-2ed2-4831-8b5e-9c766a59c4a5
begin
	index=parse(Int,patch)
	groundtruth=groundtruth_patches[index]
	sourcecolor=source_patches[index]
	destcolor=dest_patches[index]
end

# ╔═╡ 0e8135ab-8592-4e9a-b0db-41b83c5fccaf
[sWP,dWP,sourcecolor,destcolor,groundtruth]

# ╔═╡ 4c35d989-97fe-47cf-8db2-1218e4055072
begin
	gt = convert(DIN99o,groundtruth)
	sc = convert(DIN99o,sourcecolor)
	dc = convert(DIN99o,destcolor)
	plot(legend = :outertopright, ratio = 1.0, xaxis="a*99",yaxis="b*99", title="Color coordinates in DIN99o color space", size = (800,450))
	#scatter(zeros(10),zeros(10), color=Lab(50,0,0),markersize=2, shape = :hexagon,  label="neutral axis",legend=:outertopleft, aspect_ratio= :equal,background=Lab(20,0,0))
	scatter!([sc.a],[sc.b], color = sourcecolor,markersize=9, shape = :diamond, label="source color")
	scatter!([gt.a],[gt.b], color=groundtruth,markersize=10, markershape = :hexagon, label="ground truth")
	scatter!([dc.a],[dc.b], color = destcolor,markersize=8, shape = :circle, label="destination color")
	quiver!([sc.a],[sc.b],quiver=([dc.a-sc.a],[dc.b-sc.b]),arrowsize=1,marker=:none,color=:black)
end

# ╔═╡ df224fbf-e32f-47b0-824b-3ada60992734
begin
	ΔE = sqrt((gt.a - dc.a)^2 + (gt.b - dc.b)^2)
	nothing
end

# ╔═╡ 91454a50-f4ab-48b8-b9b0-51d0c73fab3b
md"ΔE $(string(typeof(gt))[12:end-9]) = $ΔE"

# ╔═╡ 28f92850-856b-4539-8ef8-111b3ca06313
md"Spectrum of color patch and light source:"

# ╔═╡ 0855f0cc-2e79-4916-a785-87c45aef0f4f
begin
	#scatter(mbf.λ,mbf.r,label="original")
	plot(colorpatch.λ,colorpatch.r, ylimits = (0.0,1.5), color=groundtruth, linewidth = 2,label="Color Checker $patch",background=DIN99o(100-gt.l,0,0),legend=:outertopright)
	#scatter!(colorpatch.λ,colorpatch.r,color=groundtruth,markersize=1, label=false)
	plot!(source.λ,source.l,color=sWP,linewidth = 4,label="light source $sourcetemp K")
	#scatter!(source.λ,source.l,color=sWP,markersize = 1, label=false)
end

# ╔═╡ 25ddb7d2-584b-4fcc-bae4-2cdf68878751
source

# ╔═╡ Cell order:
# ╠═37080ff2-a3b7-11eb-3634-af524d492479
# ╠═accd501b-c865-485b-9c8c-4d7b074b22a3
# ╠═201148e9-e528-435e-afc8-9e4390e2729d
# ╟─22394b49-4e75-444e-b19f-1bfc637dbe98
# ╠═e62085f0-b38a-43cc-a277-0504da926958
# ╟─9e05057a-216a-44eb-95cc-b15b46fc2d3d
# ╠═54e4792f-0294-4c58-931b-4c3a5ac624ed
# ╟─28c0b9e7-9faa-4814-8e1d-5c67a58dfe14
# ╟─54dfc832-dcc2-463c-9be9-ba785f0fec30
# ╟─e9951abb-915a-4799-8ed7-5c97490e8c8e
# ╟─455c4526-1dd0-461d-a40d-d4c9be7b9c2c
# ╟─c71b58d4-a969-4153-ade5-3ee27044b482
# ╟─d96a11d2-6321-4b88-9798-9809533a734f
# ╟─e94743af-b474-4801-ac3e-53703c4ab56b
# ╟─07df29f9-b7aa-4877-bda6-8e6c32c6d891
# ╟─ad9c6ea9-f6c3-415e-bc65-ca817809db8b
# ╟─27c3f7c1-afc6-4095-86a2-85c4651ab75d
# ╟─f9c57d5c-b947-445f-b479-91515cd473b8
# ╟─d128b785-6a9b-48c0-973a-4acfd24fa97c
# ╟─bab11da2-681f-41a8-ae66-c73b6acaf38b
# ╟─04857ffb-aa2e-4dd8-af49-10bcc22c4356
# ╟─b69148f5-789e-4e95-8f3b-c78ed71de534
# ╟─798d5450-00fd-4f18-a962-56374e3e0fcc
# ╟─e99970e1-5b00-42d0-a464-51cded141acf
# ╟─e59483cc-5366-40cb-8086-381915eb3531
# ╟─261ebe67-9f43-455e-ab37-51f254d26e67
# ╟─fa59297b-e992-4285-8295-b5c8273d8a9d
# ╟─c3370bb8-b574-4c82-a1a2-5cbe04e3da29
# ╟─5e8c03cc-c67a-445e-b775-b339b0beba03
# ╟─e22a766c-1ab4-4f6d-9b38-98bc579f6de9
# ╟─0e8135ab-8592-4e9a-b0db-41b83c5fccaf
# ╟─57ac1b78-5d10-4ff4-b30b-430cc111f925
# ╟─6eafc330-c225-4f26-8a17-561467a8be13
# ╟─4d5584ab-2ed2-4831-8b5e-9c766a59c4a5
# ╟─91454a50-f4ab-48b8-b9b0-51d0c73fab3b
# ╟─df224fbf-e32f-47b0-824b-3ada60992734
# ╟─4c35d989-97fe-47cf-8db2-1218e4055072
# ╟─28f92850-856b-4539-8ef8-111b3ca06313
# ╠═0855f0cc-2e79-4916-a785-87c45aef0f4f
# ╟─25ddb7d2-584b-4fcc-bae4-2cdf68878751


text:
using Plots, LaTeXStrings
import DustExtinction: ccm89_ca, ccm89_cb, od94_ca, od94_cb, cal00_invum, ccm89_invum, vcg04_invum, gcc09_invum, f99_invum, f04_invum, f19_invum, m14_invum, FM90, G16

dir = joinpath(@__DIR__, "src", "assets")

#--------------------------------------------------------------------------------
# ccm89

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = ccm89_invum.(w, rv, Ref(ccm89_ca), Ref(ccm89_cb))
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "ccm89_plot.svg"))

#--------------------------------------------------------------------------------
# od94

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = ccm89_invum.(w, rv, Ref(od94_ca), Ref(od94_cb))
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "od94_plot.svg"))

#--------------------------------------------------------------------------------
# cal00

w = range(0.46, 8.3, length=1000)
plot()
for rv in [2.0, 3.0, 4.05, 5.0, 6.0]
  m = cal00_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "cal00_plot.svg"))

#--------------------------------------------------------------------------------
# sfd98

l = range(-pi, pi, length=400)
b = range(-pi/64, pi/64, length=300)
dustmap = SFD98Map()
m = [dustmap(li, bj) for li in l, bj in b]
heatmap(l, b, m, label="", transpose=true, colorbar_title="E(B-V)")
xlabel!("l (rad)")
ylabel!("b (rad)")
savefig(joinpath(dir, "sfd98_plot.svg"))

#--------------------------------------------------------------------------------
# gcc09

w = range(3.3, 11.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = gcc09_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "gcc09_plot.svg"))

#--------------------------------------------------------------------------------
# vcg04

w = range(3.3, 8.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = vcg04_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"\mu m ^{-1}")
ylabel!("E(B-V)")
savefig(joinpath(dir, "vcg04_plot.svg"))

#--------------------------------------------------------------------------------
# F99

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = f99_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "F99_plot.svg"))

#--------------------------------------------------------------------------------
# F04

w = range(0.3, 10.0, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = f04_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "F04_plot.svg"))

#--------------------------------------------------------------------------------
# F19

w = range(0.3, 8.7, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = f19_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "F19_plot.svg"))

#--------------------------------------------------------------------------------
# M14

w = range(0.3, 3.3, length=1000)
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = m14_invum.(w, rv)
  plot!(w, m, label="Rv=$rv")
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "M14_plot.svg"))

#--------------------------------------------------------------------------------
# FM90

w = range(3.8, 8.6, step = 0.001)
x = 1e4 ./ w
plot()

m1 = FM90().(x)
plot!(w, m1, label = "total")

m2 = FM90(c3=0.0, c4=0.0).(x)
plot!(w, m2, label = "linear term")

m3 = FM90(c1=0.0, c2=0.0, c4=0.0).(x)
plot!(w, m3, label = "bump term")

m4 = FM90(c1=0.0, c2=0.0, c3=0.0).(x)
plot!(w, m4, label = "FUV rise term")

xlabel!(L"\mu m ^{-1}")
ylabel!(L"E(\lambda - V)/E(B - V)")
savefig(joinpath(dir, "FM90_plot.svg"))

#--------------------------------------------------------------------------------
# G16

# Fixed f_A = 1.0, variable Rv
w = range(0.3, 10.0, length=1000)
x = 1e4 ./ w
plot()
for rv in [2.0, 3.1, 4.0, 5.0, 6.0]
  m = G16(Rv=rv, f_A=1.0).(x)
  plot!(w, m, label="Rv=$rv", legendtitle="f_A = 1.0", legend=:topleft)
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "G16_fixed_f_A_plot.svg"))

# Fixed Rv = 3.1, variable f_A
plot()
for f_A in [0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
  m = G16(Rv=3.1, f_A=f_A).(x)
  plot!(w, m, label="f_A=$f_A", legendtitle="Rv = 3.1", legend=:topleft)
end
xlabel!(L"x\ \left[\mu m ^{-1}\right]")
ylabel!(L"A(x)/A(V)")
savefig(joinpath(dir, "G16_fixed_Rv_plot.svg"))


text:
using DelimitedFiles

datafiles = readdir("/Users/carsonc/github/covid-sicr/data")
datafiles = readdir("/Users/carsonc/github/VaxBRAModelEPI/USStates/")

datafiles = datafiles[occursin.("covidtimeseries_US_",datafiles)]

regions = Array{String,1}(undef,0)
for datafile in datafiles
    push!(regions,replace(replace(datafile,".csv"=>""),"covidtimeseries_"=>""))
end

# models = ["fulllinearmodel","reducedlinearmodelq0","reducedlinearmodelq0ctime","reducedlinearmodelNegBinom","nonlinearmodelq0ctime","nonlinearmodel"]
# models = ["nonlinearmodelq0ctime"]
# models = ["nonlinearmodel"]
# models = ["linearm2R","linearmqc2R"]
models = ["SICRM","SICRLM","SICRMQC","SICRLMQC","SICRM2R","SICRMQC2R","SICRLM2R","SICRLMQC2R"]
models = ["SICRMQ","SICRLMQ","SICRMQ2R","SICRLMQ2R"]
models = ["SICRMQC2R2DX2"]
models = ["SICR2HV"]

for modelname in models

    iter = "-it=10000"
    warm = "-wm=6000"
    format = "-f=1"
    fitpath = "-fp='/data/carsonc/covid-19/fitUS06-27-20-2HV/'"
    modelpath = "-mp='/home/carsonc/covid-sicr/models/'"
    # datapath = "-dp='/home/carsonc/covid-sicr/data/'"
    datapath = "-dp='/home/carsonc/VaxBRAModelEPI/USStates/'"
    ad = "-ad=.85"
    ft = "-ft=1"
    ld = "-ld='06/27/20'"
    ch = "-ch=4"
    runfile = "&& python /home/carsonc/covid-sicr/scripts/runV.py"

    swarmname = modelname * ".swarm"

    f = open(swarmname,"w")

    for region in regions
        # init = "-i='/data/carsonc/covid-19/fit10/" * modelname * "_" * region * ".pkl'"
        init = ""
        regionarg = "-r='" * region * "'"
        writedlm(f,["source /data/carsonc/conda/etc/profile.d/conda.sh \\" ])
        writedlm(f,[ "&& conda activate covid-19 \\" ])
        writedlm(f,[ runfile modelname regionarg modelpath datapath fitpath iter warm ch format ad ft ld])
        # writedlm(f,[ "/data/carsonc/conda/envs/covid-19/bin/python run.py" modelname region fitpath iter warm format ])
    end

    close(f)
end

max_stars_repo_name: nih-niddk-mbs/covid-sicr

max_issues_repo_name: nih-niddk-mbs/covid-sicr

max_forks_repo_name: nih-niddk-mbs/covid-sicr


text:
"""
    MetaGraph{T,Label,Graph,VertexMeta,EdgeMeta,GraphMeta,WeightFunction,U} <: AbstractGraph{T}

A graph type with custom vertex labels containing vertex-, edge- and graph-level metadata.

Vertex labels have type `Label`, while vertex (resp. edge, resp. graph) metadata has type `VertexMeta` (resp. `EdgeMeta`, resp. `GraphMeta`).
It is recommended not to set `Label` to `Int` to avoid confusion between vertex labels and vertex codes (which have type `T <: Integer`).

# Fields
- `g::Graph`: underlying, data-less graph with vertex indices of type `T`
- `vprops::Dict{Label,Tuple{T,VertexMeta}}`: dictionary mapping vertex labels to vertex codes and metadata
- `eprops::Dict{Tuple{Label,Label},EdgeMeta}`: dictionary mapping edge labels such as `(label_u, label_v)` to edge metadata
- `gprops::GraphMeta`: graph metadata
- `weightfunction::WeightFunction`: function defining edge weight from edge metadata
- `defaultweight::U`: default weight for the edges
- `metaindex::Dict{T,Label}`: dictionary mapping vertex codes to vertex labels
"""
struct MetaGraph{
    T<:Integer,
    Label,
    Graph,
    VertexMeta,
    EdgeMeta,
    GraphMeta,
    WeightFunction,
    U<:Real,
} <: AbstractGraph{T}
    graph::Graph
    vprops::Dict{Label,Tuple{T,VertexMeta}}
    eprops::Dict{Tuple{Label,Label},EdgeMeta}
    gprops::GraphMeta
    weightfunction::WeightFunction
    defaultweight::U
    metaindex::Dict{T,Label}
end

"""
    MetaGraph(g;
        Label = Symbol,
        VertexMeta = nothing,
        EdgeMeta = nothing,
        gprops = nothing,
        weightfunction = eprops -> 1.0,
        defaultweight = 1.0
    )
"""
function MetaGraph(
    g::AbstractGraph{T};
    Label = Symbol,
    VertexMeta = Nothing,
    EdgeMeta = Nothing,
    gprops = nothing,
    weightfunction = eprops -> 1.0,
    defaultweight = 1.0,
) where {Vertex,T}
    MetaGraph(
        g,
        Dict{Label,Tuple{T,VertexMeta}}(),
        Dict{Tuple{Label,Label},EdgeMeta}(),
        gprops,
        weightfunction,
        defaultweight,
        Dict{T,Label}(),
    )
end

zero(
    g::MetaGraph{T,Label,Graph,VertexMeta,EdgeMeta,GraphMeta},
) where {T,Label,Graph,VertexMeta,EdgeMeta,GraphMeta} = MetaGraph(
    Graph();
    Label = Label,
    VertexMeta = VertexMeta,
    EdgeMeta = EdgeMeta,
    gprops = g.gprops,
    weightfunction = g.weightfunction,
    defaultweight = g.defaultweight,
)

==(x::MetaGraph, y::MetaGraph) = x.graph == y.graph

copy(g::T) where {T<:MetaGraph} = deepcopy(g)

function show(
    io::IO,
    g::MetaGraph{<:Any,Label,<:Any,VertexMeta,EdgeMeta},
) where {Label,VertexMeta,EdgeMeta}
    print(
        io,
        "Meta graph based on a $(g.graph) with vertices indexed by $Label(s), $VertexMeta(s) vertex metadata, $EdgeMeta(s) edge metadata, $(repr(g.gprops)) as graph metadata, and default weight $(g.defaultweight)",
    )
end


text:
# Much of this is copied from Keno's CXX.jl

import Base: LineEdit, REPL

function create_mode(prompt, name;
    on_enter=LineEdit.default_enter_cb,
    on_done=()->nothing,
    prefix_color=Base.text_colors[:green],
    suffix_color=Base.text_colors[:red],
    repl=Base.active_repl,
    main_mode = repl.interface.modes[1])

    mirepl = isdefined(repl,:mi) ? repl.mi : repl

    done_func = REPL.respond(on_done, repl, main_mode)

    # Setup mode
    mode = LineEdit.Prompt(prompt;
        # Copy colors from the prompt object
        prompt_prefix=prefix_color,
        prompt_suffix=suffix_color,
        on_enter = on_enter,
        on_done = done_func)

    main_mode == mirepl.interface.modes[1] &&
            push!(mirepl.interface.modes,mode)

    hp = main_mode.hist
    hp.mode_mapping[name] = mode
    mode.hist = hp

    search_prompt, skeymap = LineEdit.setup_search_keymap(hp)
    mk = REPL.mode_keymap(main_mode)

    b = Dict{Any,Any}[skeymap, mk, LineEdit.history_keymap, LineEdit.default_keymap, LineEdit.escape_defaults]
    mode.keymap_dict = LineEdit.keymap(b)

    return mode
end

function add_mode_to_repl(key, mode)
    repl = Base.active_repl
    mirepl = isdefined(repl,:mi) ? repl.mi : repl
    main_mode = mirepl.interface.modes[1]

    const new_keymap = Dict{Any,Any}(
            key => function (s,args...)
                if isempty(s) || position(LineEdit.buffer(s)) == 0
                    buf = copy(LineEdit.buffer(s))
                    LineEdit.transition(s, mode) do
                        LineEdit.state(s, mode).input_buffer = buf
                    end
                else
                    LineEdit.edit_insert(s,key)
                end
            end
        )

    main_mode.keymap_dict = LineEdit.keymap_merge(main_mode.keymap_dict, new_keymap);
    nothing
end


text:
# ---
# title: 1171. Remove Zero Sum Consecutive Nodes from Linked List
# id: problem1171
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Linked List
# link: <https://leetcode.com/problems/remove-zero-sum-consecutive-nodes-from-linked-list/description/>
# hidden: true
# ---
#
# Given the `head` of a linked list, we repeatedly delete consecutive sequences
# of nodes that sum to `0` until there are no such sequences.
#
# After doing so, return the head of the final linked list.  You may return any
# such answer.
#
#
#
# (Note that in the examples below, all sequences are serializations of
# `ListNode` objects.)
#
# **Example 1:**
#
#
#
#     Input: head = [1,2,-3,3,1]
#     Output: [3,1]
#     **Note:** The answer [1,2,1] would also be accepted.
#
#
# **Example 2:**
#
#
#
#     Input: head = [1,2,3,-3,4]
#     Output: [1,2,4]
#
#
# **Example 3:**
#
#
#
#     Input: head = [1,2,3,-3,-2]
#     Output: [1]
#
#
#
#
# **Constraints:**
#
#   * The given linked list will contain between `1` and `1000` nodes.
#   * Each node in the linked list has `-1000 <= node.val <= 1000`.
#
#
## @lc code=start
using LeetCode

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


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/display.jl
# Functions used to print information about solution routine to console.
# Printing is done with reference to the input problem is there is any
# ambiguity.
#############################################################################

const PRINTING_IOFORMAT = :SCI
const PRINTING_CHARSET = :ASCII

"""
$(FUNCTIONNAME)

Prints solution information for the B&B problem. Displays first node found, solution value,
solution, and time spent solving subproblems.
"""
function print_solution!(m::Optimizer)
    if m._parameters.verbosity > 0
        println(" ")
        println("First Solution Found at Node $(m._first_solution_node)")
        lower_bound = m._min_converged_value
        if (m._input_problem._optimization_sense !== MOI.MIN_SENSE)
            lower_bound *= -1.0
        end
        println("LBD = $(lower_bound)")
        println("UBD = $(MOI.get(m, MOI.ObjectiveValue()))")
        println("Solution is :")
        if m._feasible_solution_found
            for i = 1:m._input_problem._variable_count
                println("    X[$i] = $(m._continuous_solution[i])")
            end
        end
        println(" ")
     end
     return
end

"""
$(FUNCTIONNAME)

Prints node information for the B&B problem. Node id, bound, and interval box.
"""
function print_node!(m::Optimizer)
    n = m._current_node
    bound = (m._input_problem._optimization_sense === MOI.MIN_SENSE) ? n.lower_bound : -n.lower_bound
    println(" ")
    println("Node ID: $(n.id), Lower Bound: $(bound), Lower Variable Bounds:
             $(n.lower_variable_bounds), Upper Variable Bounds: $(n.upper_variable_bounds)")
    println(" ")
    return
end

"""
$(FUNCTIONNAME)

Prints the iteration information based on verbosity. The header is displayed
every `header_interval`, the iteration info is displayed every `iteration_interval`.
"""
function print_iteration!(m::Optimizer)

    if m._parameters.verbosity > 0

        # prints header line every B.hdr_intv times
        if mod(m._iteration_count, m._parameters.header_iterations) === 0 || m._iteration_count === 1
            println("-----------------------------------------------------------------------------------------------------------------------------")
            println("|  Iteration #  |     Nodes    | Lower Bound  |  Upper Bound  |      Gap     |     Ratio    |     Time     |    Time Left   |")
            println("-----------------------------------------------------------------------------------------------------------------------------")
        end

        # prints iteration summary every B.itr_intv times
        if mod(m._iteration_count, m._parameters.output_iterations) === 0

            print_str = "| "

            max_len = 12
            temp_str = string(m._iteration_count)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  | "

            max_len = 12
            temp_str = string(m._node_count)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 12
            if m._input_problem._optimization_sense === MOI.MIN_SENSE
                lower = m._global_lower_bound
                upper = m._global_upper_bound
            else
                lower = -m._global_upper_bound
                upper = -m._global_lower_bound
            end
            #temp_str = string(round(lower, sigdigits = 5))
            #temp_str = string(lower, sigdigits = 3))
            temp_str = formatted(lower, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            #temp_str = formatted(upper, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            #temp_str = string(upper, sigdigits = 3))
            temp_str = formatted(upper, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  |"

            max_len = 12
            #temp_str = string(round(abs(x._global_upper_bound - x._global_lower_bound), sigdigits = 3))
            temp_str = formatted(abs(m._global_upper_bound - m._global_lower_bound), PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  | "

            max_len = 12
            #temp_str = string(round(relative_gap(x._global_lower_bound, x._global_upper_bound), sigdigits = 3))
            temp_str = formatted(relative_gap(m._global_lower_bound, m._global_upper_bound), PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" | "

            max_len = 12
            #temp_str = string(round(x._run_time, sigdigits = 3))
            temp_str = formatted(m._run_time, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*" |  "

            max_len = 12
            #temp_str = string(round(x._time_left, sigdigits = 4))
            temp_str = formatted(m._time_left, PRINTING_IOFORMAT, ndigits=4, charset=PRINTING_CHARSET)
            len_str = length(temp_str)
            print_str *= (" "^(max_len - len_str))*temp_str*"  |"

            println(print_str)
        end
    end

    return
end

"""
$(FUNCTIONNAME)

Prints the results of a single bounding problem.
"""
function print_results!(m::Optimizer, flag::Bool)
    if m._parameters.verbosity > 1
        println(" ")
        if flag
            obj_val = m._lower_objective_value
            if m._input_problem._optimization_sense === MOI.MIN_SENSE
                print("Lower Bound (First Iteration): $(obj_val),")
            else
                print("Upper Bound (First Iteration): $(-obj_val),")
            end
            print(" Solution: $(m._lower_solution), Feasibility: $(m._lower_feasibility)\n")
            println("Termination Status Code: $(m._lower_termination_status)")
            println("Result Code: $(m._lower_result_status)")
        else
            obj_val = m._upper_objective_value
            if m._input_problem._optimization_sense === MOI.MIN_SENSE
                print("Upper Bound: $(obj_val), ")
            else
                print("Lower Bound: $(-obj_val), ")
            end
            print(" Solution: $(m._upper_solution), Feasibility: $(m._upper_feasibility)\n")
            println("Termination Status Code: $(m._upper_termination_status)")
            println("Result Code: $(m._upper_result_status)")
        end
        println(" ")
    end
    return
end

"""
$(FUNCTIONNAME)

Prints the results after performing various cuts.
"""
function print_results_post_cut!(m::Optimizer)
    if m._parameters.verbosity > 1
        println(" ")
        if m._input_problem._optimization_sense === MOI.MIN_SENSE
            print("Lower Bound (Last Iteration): $(m._lower_objective_value)")
        else
            print("Upper Bound (Last Iteration): $(-m._lower_objective_value)")
        end
        print(", Solution: $(m._lower_solution), Feasibility: $(m._lower_feasibility)\n")
        println(" ")
    end
    return
end


text:
logpdf0(x, P) = logdensity(Gaussian{(:Σ,)}(P), x)

function backward(::BF, k::Union{AffineGaussianKernel,LinearGaussianKernel}, q::Gaussian{(:μ,:Σ)})
    ν, Σ = q.μ, q.Σ
    B, β, Q = params(k)
    B⁻¹ = inv(B)
    νp = B⁻¹*(ν - β)
    Σp = B⁻¹*(Σ + Q)*B⁻¹'
    q, Gaussian{(:μ,:Σ)}(νp, Σp)
end

function backward(::BF, k::ConstantGaussianKernel, q::Gaussian{(:F,:Γ)})
    message = q
    message, nothing
end

function backward(::BF, k::Union{AffineGaussianKernel,LinearGaussianKernel}, q::Gaussian{(:F,:Γ)})
    @unpack F, Γ = q
    # Theorem 7.1 [Automatic BFFG]
    B, β, Q = params(k)
    Σ = inv(Γ) # requires invertibility of Γ
    K = B'*inv(Σ + Q)
    ν̃ = Σ*F - β
    Fp = K*ν̃
    Γp = K*B
    message = q
    message, Gaussian{(:F,:Γ)}(Fp, Γp)
end

function backward(::BF, k::Union{AffineGaussianKernel,LinearGaussianKernel}, y)
    # Theorem 7.1 [Automatic BFFG]
    B, β, Q = params(k)
    K = B'/Q
    Fp = K*(y - β)
    Γp = K*B
    message = Leaf(y)
    message, Gaussian{(:F,:Γ)}(Fp, Γp)
end

function backward(::BFFG, k::Union{AffineGaussianKernel,LinearGaussianKernel}, q::WGaussian{(:F,:Γ,:c)}; unfused=false)
    @unpack F, Γ, c = q
    # Theorem 7.1 [Automatic BFFG]
    B, β, Q = params(k)
    Σ = inv(Γ) # requires invertibility of Γ
    K = B'/(Σ + Q)
    ν̃ = Σ*F - β
    Fp = K*ν̃
    Γp = K*B
    # Corollary 7.2 [Automatic BFFG]
    if !unfused
        cp = c - logdet(B)
    else
        cp = c + logpdf0(ν̃, Σ + Q)
    end
    message = q
    message, WGaussian{(:F,:Γ,:c)}(Fp, Γp, cp)
end

function backward(::BFFG, k::Union{AffineGaussianKernel,LinearGaussianKernel}, y; unfused=false)
    # Theorem 7.1 [Automatic BFFG]
    B, β, Q = params(k)
    K = B'/Q
    Fp = K*(y - β)
    Γp = K*B
    # Corollary 7.2 [Automatic BFFG]
    if !unfused
        cp = -logdet(B)
    else
        cp = +logpdf0(y - β, Q)
    end
    message = Leaf(y)
    message, WGaussian{(:F,:Γ,:c)}(Fp, Γp, cp)
end

function forward(::BF, k::Union{AffineGaussianKernel,LinearGaussianKernel}, m::Gaussian{(:F,:Γ)})
    @unpack F, Γ = m
    B, β, Q = params(k)

    Q⁻ = inv(Q)
    Qᵒ = inv(Q⁻ + Γ)
    Bᵒ = Qᵒ*Q⁻*B
    βᵒ = Qᵒ*(Q⁻*β + F)

    kernel(Gaussian; μ=AffineMap(Bᵒ, βᵒ), Σ=ConstantMap(Qᵒ))
end
function forward(::BF, k::ConstantGaussianKernel, m::Gaussian{(:F,:Γ)})
    @unpack F, Γ = m
    β, Q = params(k)

    Q⁻ = inv(Q)
    Qᵒ = inv(Q⁻ + Γ)
    βᵒ = Qᵒ*(Q⁻*β + F)

    kernel(Gaussian; μ=ConstantMap(βᵒ), Σ=ConstantMap(Qᵒ))
end

function forward(::BFFG, k::Union{AffineGaussianKernel,LinearGaussianKernel}, m::WGaussian{(:F,:Γ,:c)})
    @unpack F, Γ, c = m
    B, β, Q = params(k)

    Q⁻ = inv(Q)
    Qᵒ = inv(Q⁻ + Γ)
    #μᵒ = Qᵒ*(Q⁻*(B*x + β) + F )
    Bᵒ = Qᵒ*Q⁻*B
    βᵒ = Qᵒ*(Q⁻*β + F)

    kernel(WGaussian; μ=AffineMap(Bᵒ, βᵒ), Σ=ConstantMap(Qᵒ), c=ConstantMap(0.0))
end

function forward(::BFFG, k::GaussKernel, k̃::Union{AffineGaussianKernel,LinearGaussianKernel}, m::WGaussian{(:F,:Γ,:c)}, x)
    @unpack F, Γ, c = m
    μ, Q = k.ops
    μ̃, Q̃ = k̃.ops

    # Proposition 7.3.
    Q⁻ = inv(Q(x))
    Qᵒ = inv(Q⁻ + Γ)
    μᵒ = Qᵒ*(Q⁻*(μ(x)) + F)

    Q̃⁻ = inv(Q̃(x))
    Q̃ᵒ = inv(Q̃⁻ + Γ)
    μ̃ᵒ = Q̃ᵒ*(Q̃⁻*(μ̃(x)) + F)

    c = logpdf0(μ(x), Q(x)) - logpdf0(μ̃(x), Q̃(x))
    c += logpdf0(μ̃ᵒ, Q̃ᵒ) - logpdf0(μᵒ, Qᵒ)

    WGaussian{(:μ,:Σ,:c)}(μᵒ, Qᵒ, c)
end

function backward(::BFFG, ::Copy, args::WGaussian{(:μ,:Σ,:c)}...; unfused=true)
    F, H, c = params(convert(WGaussian{(:F,:Γ,:c)}, args[1]))
    unfused || (c += logdensity(Gaussian{(:F,:Γ)}(F, H), 0F))
    for b in args[2:end]
        F2, H2, c2 = params(convert(WGaussian{(:F,:Γ,:c)}, b))
        F += F2
        H += H2
        c += c2
        unfused || (c += logdensity(Gaussian{(:F,:Γ)}(F2, H2), 0F2))
    end
    Δ = -logdensity(Gaussian{(:F,:Γ)}(F, H), 0F)
    m = ()
    m, convert(WGaussian{(:μ,:Σ,:c)}, WGaussian{(:F,:Γ,:c)}(F, H, Δ + c))
end

function backward(::Union{BFFG,BF}, ::Copy, a::Gaussian{(:F,:Γ)}, args...)
    F, H = params(a)
    for b in args
        F2, H2 = params(b::Gaussian{(:F,:Γ)})
        F += F2
        H += H2
    end
    m = ()
    m, Gaussian{(:F,:Γ)}(F, H)
end

function backward(::BFFG, ::Copy, a::WGaussian{(:F,:Γ,:c)}, args...; unfused=true)
    F, H, c = params(a)
    unfused || (c += logdensity(Gaussian{(:F,:Γ)}(F, H), 0F))
    for b in args
        F2, H2, c2 = params(b::WGaussian{(:F,:Γ,:c)})
        F += F2
        H += H2
        c += c2
        unfused || (c += logdensity(Gaussian{(:F,:Γ)}(F2, H2), 0F2))
    end
    Δ = -logdensity(Gaussian{(:F,:Γ)}(F, H), 0F)
    m = ()
    m, WGaussian{(:F,:Γ,:c)}(F, H, Δ + c)
end

function forward(::BFFG, k::Union{AffineGaussianKernel,LinearGaussianKernel}, y::Leaf, x::Weighted)
    Dirac(weighted(y[], x.ll))
end
function forward(::BFFG, ::Copy{2}, _, x::Weighted)
    MeasureTheory.Dirac((x, weighted(x[])))
end


text:
function toric_ed(A)
    d, n = size(A)
    @polyvar t[1:d] y[1:n] u[1:n]

    φ = map(j -> prod(i -> t[i]^A[i,j], 1:d), 1:n)
    Dφ = [differentiate(φ[i], t[j]) for i=1:n, j=1:d]

    F = [φ + y - u; Dφ' * y]

    # We sample a random t, find a normal vector on this point and then assemble
    # an u where we know the solution to.
    t_rand = randn(Float64, size(A, 1))
    null = nullspace(map(fi -> fi(t => t_rand), Dφ'))
    y_rand = null * randn(Float64, size(null, 2))
    u₀ = map(fi -> fi(t => t_rand), φ) + y_rand

    x₀ = complex.([t_rand; y_rand])

    F, u, u₀, x₀
end

@testset "Monodromy" begin

    @testset "monodromy_solve" begin
        F, p, p₀, x₀ = toric_ed([3 2 1 0; 0 1 2 3])

        # test that timeout works
        Random.seed!(51232)
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21, timeout=1e-12)
        @test length(result.solutions) < 21

        result = monodromy_solve(F, x₀, p₀, parameters=p,
                target_solutions_count=21,
                maximal_number_of_iterations_without_progress=20)
        @test result.returncode == :success
        @test length(solutions(result)) == 21
        @test length(solutions(result, onlyreal = true)) >= 1
        @test result.statistics.ntrackedpaths ≥ 21
        @test result.statistics.nparametergenerations ≥ 1
        @test length(HC.UniquePoints(result.solutions).points) == 21
        @test isempty(sprint(show, result)) == false

        # test input of length > 1
        result = monodromy_solve(F, [x₀ for _ in 1:30], p₀, parameters=p)
        @test length(solutions(result)) == 21

        # test with false input
        result = monodromy_solve(F, [rand(6) for _ in 1:10], p₀, parameters=p)
        @test result.returncode == :invalid_startvalue
        result = monodromy_solve(F, [x₀, rand(6)], p₀, parameters=p)
        @test length(solutions(result)) == 21
        result = monodromy_solve(F, [x₀, rand(6)], p₀, parameters=p, check_startsolutions = false)
        @test result.returncode == :invalid_startvalue

        # different distance function
        result = monodromy_solve(F, x₀, p₀, parameters=p, distance = (x,y) -> 0.0)
        @test length(solutions(result)) == 1

        # Test stop heuristic using too high target_solutions_count
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=25)
        @test result.returncode == :heuristic_stop
        # Test stop heuristic with no target solutions count
        result = monodromy_solve(F, x₀, p₀, parameters=p)
        @test result.returncode == :heuristic_stop
        # Test stop heuristic with no target solutions count
        result = monodromy_solve(F, x₀, p₀, parameters=p, strategy=Triangle(useweights=true))
        @test result.returncode == :heuristic_stop

        # By group_actions=nothing we force that complex conjugation is not used.
        result2 = monodromy_solve(F, x₀, p₀, parameters=p,
                        target_solutions_count=21, complex_conjugation=false,
                        maximal_number_of_iterations_without_progress=100)
        @test result2.returncode == :success

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
                done_callback=((_, _) -> true))
        @test length(result.solutions) == 2

        roots_of_unity(s) = begin
            t = cis(π*2/3)
            t² = t * t
            (vcat(t * s[1], t * s[2], s[3:end]),
             vcat(t² * s[1], t² * s[2], s[3:end]))
        end

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            maximal_number_of_iterations_without_progress=100,
            equivalence_classes=false,
            group_action=roots_of_unity)
        @test length(result.solutions) == 21

        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            maximal_number_of_iterations_without_progress=100,
            complex_conjugation=false, # disable complex conjugation to test it as a group action.
            equivalence_classes=false,
            group_actions=(roots_of_unity, s -> (conj.(s),)))
        @test length(result.solutions) == 21
        @test length(solutions(result)) == 21
        @test length(realsolutions(result)) < 21
        test_treeviews(result)

        # group_actions as a vector
        result = monodromy_solve(F, x₀, p₀, parameters=p, target_solutions_count=21,
            maximal_number_of_iterations_without_progress=100,
            complex_conjugation=false,
            equivalence_classes=false,
            group_actions=[roots_of_unity, s -> (conj.(s),)])
        @test length(result.solutions) == 21

        # equivalence classes
        result = monodromy_solve(F, x₀, p₀, parameters=p,
            equivalence_classes=true,
            target_solutions_count=7,
            maximal_number_of_iterations_without_progress=100,
            group_actions=roots_of_unity)
        @test length(result.solutions) == 7
        # Test that equivalence classes are on by default if we supply a group action
        result = monodromy_solve(F, x₀, p₀, parameters=p,
                            group_action=roots_of_unity,
                            maximal_number_of_iterations_without_progress=20)
        @test length(result.solutions) == 7

        # Test affine tracking
        result = monodromy_solve(F, x₀, p₀, parameters=p, affine_tracking=true,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        maximal_number_of_iterations_without_progress=200)
        @test length(result.solutions) == 7

        # AbstractSystem as input
        F_p = SPSystem(F; parameters=p)
        result = monodromy_solve(F_p, x₀, p₀, affine_tracking=true,
                        group_action=roots_of_unity,
                        target_solutions_count=7,
                        maximal_number_of_iterations_without_progress=200)
        @test length(result.solutions) == 7
    end

    @testset "Method of Moments" begin
        Random.seed!(130793)
        @polyvar a[1:3] x[1:3] s[1:3]

        f0 = a[1]+a[2]+a[3];
        f1 = a[1]*x[1]+a[2]*x[2]+a[3]*x[3];
        f2 = a[1]*(x[1]^2+s[1])+a[2]*(x[2]^2+s[2])+a[3]*(x[3]^2+s[3]);
        f3 = a[1]*(x[1]^3+3*s[1]*x[1])+a[2]*(x[2]^3+3*s[2]*x[2])+a[3]*(x[3]^3+3*s[3]*x[3]);
        f4 = a[1]*(x[1]^4+6*s[1]*x[1]^2+3*s[1]^2)+a[2]*(x[2]^4+6*s[2]*x[2]^2+3*s[2]^2)+
        	a[3]*(x[3]^4+6*s[3]*x[3]^2+3*s[3]^2);
        f5= a[1]*(x[1]^5+10*s[1]*x[1]^3+15*x[1]*s[1]^2)+a[2]*(x[2]^5+10*s[2]*x[2]^3+15*x[2]*s[2]^2)+a[3]*(x[3]^5+10*s[3]*x[3]^3+15*x[3]*s[3]^2);
        f6= a[1]*(x[1]^6+15*s[1]*x[1]^4+45*x[1]^2*s[1]^2+15*s[1]^3)+
        	a[2]*(x[2]^6+15*s[2]*x[2]^4+45*x[2]^2*s[2]^2+15*s[2]^3)+
        	a[3]*(x[3]^6+15*s[3]*x[3]^4+45*x[3]^2*s[3]^2+15*s[3]^3);
        f7= a[1]*(x[1]^7+21*s[1]*x[1]^5+105*x[1]^3*s[1]^2+105*x[1]*s[1]^3)+
        	a[2]*(x[2]^7+21*s[2]*x[2]^5+105*x[2]^3*s[2]^2+105*x[2]*s[2]^3)+
        	a[3]*(x[3]^7+21*s[3]*x[3]^5+105*x[3]^3*s[3]^2+105*x[3]*s[3]^3);
        f8= a[1]*(x[1]^8+28*s[1]*x[1]^6+210*x[1]^4*s[1]^2+420*x[1]^2*s[1]^3+105*s[1]^4)+
        	a[2]*(x[2]^8+28*s[2]*x[2]^6+210*x[2]^4*s[2]^2+420*x[2]^2*s[2]^3+105*s[2]^4)+
        	a[3]*(x[3]^8+28*s[3]*x[3]^6+210*x[3]^4*s[3]^2+420*x[3]^2*s[3]^3+105*s[3]^4)

        f = [f0, f1, f2, f3, f4, f5, f6, f7, f8]

        S₃ = SymmetricGroup(3)
        relabeling(v) = map(p -> (v[1:3][p]..., v[4:6][p]..., v[7:9][p]...), S₃)

        @polyvar p[1:9]

        function sample_moments(_=nothing)
        	a₀ = rand(3)
        	a₀ ./= sum(a₀)
        	x₀ = randn(3)
        	s₀ = rand(3) .* 2

        	y₀ = [a₀; x₀; s₀]
        	p₀ = [fᵢ(a => a₀, x => x₀, s => s₀) for fᵢ in f]
        	y₀, p₀
        end

        y₀, p₀ = sample_moments()

        R = monodromy_solve(f - p, y₀, p₀;
        			parameters=p, group_action=relabeling,
        			parameter_sampler=last ∘ sample_moments,
                    show_progress=false,
                    maximal_number_of_iterations_without_progress=5)
        @test length(solutions(R)) ≤ 225

        R = monodromy_solve(f - p, y₀, p₀;
        			parameters=p, group_action=relabeling,
        			target_solutions_count=225,
                    affine_tracking=true,
        			parameter_sampler=last ∘ sample_moments,
                    show_progress=false,
                    maximal_number_of_iterations_without_progress=20)
        @test length(solutions(R)) == 225
    end
end


text:
@testset "Frankot" begin
    p, q = synthetic_gradient(SynthSphere(), radius = 50)
    Z = convert_gradient(Frankot(), p, q)
    @test Z[38,38] ≈ 10.2715249411254
    @test Z[1,1] ≈ 10.647589591760878
    @test Z[76,76] ≈ 36.6976320484175
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitorSuccessThreshold <: SwaggerModel
    checksFailedPercent::Any # spec type: Union{ Nothing, Int32 } # spec name: checksFailedPercent
    roundTripTimeMs::Any # spec type: Union{ Nothing, Float32 } # spec name: roundTripTimeMs

    function ConnectionMonitorSuccessThreshold(;checksFailedPercent=nothing, roundTripTimeMs=nothing)
        o = new()
        validate_property(ConnectionMonitorSuccessThreshold, Symbol("checksFailedPercent"), checksFailedPercent)
        setfield!(o, Symbol("checksFailedPercent"), checksFailedPercent)
        validate_property(ConnectionMonitorSuccessThreshold, Symbol("roundTripTimeMs"), roundTripTimeMs)
        setfield!(o, Symbol("roundTripTimeMs"), roundTripTimeMs)
        o
    end
end # type ConnectionMonitorSuccessThreshold

const _property_map_ConnectionMonitorSuccessThreshold = Dict{Symbol,Symbol}(Symbol("checksFailedPercent")=>Symbol("checksFailedPercent"), Symbol("roundTripTimeMs")=>Symbol("roundTripTimeMs"))
const _property_types_ConnectionMonitorSuccessThreshold = Dict{Symbol,String}(Symbol("checksFailedPercent")=>"Int32", Symbol("roundTripTimeMs")=>"Float32")
Base.propertynames(::Type{ ConnectionMonitorSuccessThreshold }) = collect(keys(_property_map_ConnectionMonitorSuccessThreshold))
Swagger.property_type(::Type{ ConnectionMonitorSuccessThreshold }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorSuccessThreshold[name]))}
Swagger.field_name(::Type{ ConnectionMonitorSuccessThreshold }, property_name::Symbol) =  _property_map_ConnectionMonitorSuccessThreshold[property_name]

function check_required(o::ConnectionMonitorSuccessThreshold)
    true
end

function validate_property(::Type{ ConnectionMonitorSuccessThreshold }, name::Symbol, val)
end