text:
"""
    fit!(m::WSVarLmmModel,
    solver=IpoptSolver(print_level=0, mehrotra_algorithm="yes", max_iter=100);
    init=init_ls!(m), runs = 2)

Fit a `WSVarLmmModel` object using a weighted NLS method.

# Positional arguments
- `m::WSVarLmmModel`: Model to fit.
- `solver`: Nonlinear programming solver to use. Common choices include:
    - `Ipopt.IpoptSolver(print_level=0, mehrotra_algorithm="yes", warm_start_init_point="yes", max_iter=100)`.
    - `Ipopt.IpoptSolver(print_level=0, watchdog_shortened_iter_trigger=3, max_iter=100)`.
    - `Ipopt.IpoptSolver(print_level=0, max_iter=100)`.
    - `KNITRO.KnitroSolver(outlev=3)`. (Knitro is commercial software)
    - `NLopt.NLoptSolver(algorithm=:LD_MMA, maxeval=4000)`.
    - `NLopt.NLoptSolver(algorithm=:LD_LBFGS, maxeval=4000)`.

# Keyword arguments
- `init`: Initialization strategy. `fit!` will use `m.τ` and `m.Lγ` to set the
    weight matrices `Vi` and solve the weighted NLS to obtain an
    estimate for `m.β`, `m.τ`, and `m.Lγ`.  Choices for `init` include
    - `init_ls!(m)` (default): initialize by the least squares analytical solution.
    - `init_mom!(m)`: initialize by the unweighted NLS (MoM).
    - `m`: initilize from user supplied values in `m.τ` and `m.Lγ`.
- `runs::Integer`: Number of weighted NLS runs; default is 2. Each run will use the
    newest `m.τ` and `m.Lγ` to update the weight matrices `Vi` and solve the
    new weighted NLS.
- `parallel::Bool`: Multi-threading or not. Default is `false`.
- `verbose::Bool`: Verbose display or not, Default is `true`.
"""
function fit!(
    m        :: WSVarLmmModel{T},
    solver   :: MOI.AbstractOptimizer = Ipopt.Optimizer();
    solver_config::Dict =
        Dict("print_level"           => 0,
             "mehrotra_algorithm"    => "yes",
             "warm_start_init_point" => "yes",
             "max_iter"              => 100),
    init     :: WSVarLmmModel{T} = init_ls!(m),
    runs     :: Integer = 2,
    parallel :: Bool = false,
    verbose  :: Bool = true
    ) where {T<:BlasReal}
    solvertype = typeof(solver)
    solvertype <: Ipopt.Optimizer ||
        @warn("Optimizer object is $solvertype, `solver_config` may need to be defined.")
    # Pass options to solver
    config_solver(solver, solver_config)
    # set up NLP optimization problem
    npar = m.l + ◺(m.q)
    MOI.empty!(solver)
    lb   = T[]
    ub   = T[]

    NLPBlock = MOI.NLPBlockData(
        MOI.NLPBoundsPair.(lb, ub), m, true
    )

    par0 = Vector{T}(undef, npar)
    modelpar_to_optimpar!(par0, m)

    solver_pars = MOI.add_variables(solver, npar)
    for i in 1:npar
        MOI.set(solver, MOI.VariablePrimalStart(), solver_pars[i], par0[i])
    end

    MOI.set(solver, MOI.NLPBlock(), NLPBlock)
    MOI.set(solver, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    # optimize weighted NLS
    m.iswtnls[1] = true
    m.ismthrd[1] = parallel
    xsol = similar(par0)
    βprev, τprev, Lγprev = similar(m.β), similar(m.τ), similar(m.Lγ)
    for run in 1:runs
        copyto!(βprev, m.β); copyto!(τprev, m.τ), copyto!(Lγprev, m.Lγ)
        tic = time() # start timing
        # update Vi, then β and residuals with WLS
        update_wtmat!(m)
        # update τ and Lγ by WNLS
        modelpar_to_optimpar!(par0, m)
        for i in 1:npar
            MOI.set(solver, MOI.VariablePrimalStart(), solver_pars[i], par0[i])
        end
        MOI.optimize!(solver)
        toc = time()
        optstat = MOI.get(solver, MOI.TerminationStatus())
        optstat in (MOI.LOCALLY_SOLVED, MOI.ALMOST_LOCALLY_SOLVED) ||
            @warn("Optimization unsuccessful; got $optstat")
        # Get solver solution values
        fill!(xsol, zero(T))
        for i in eachindex(xsol)
            xsol[i] = MOI.get(solver, MOI.VariablePrimal(), MOI.VariableIndex(i))
        end
        optimpar_to_modelpar!(m, xsol)
        verbose && @printf(
            "run = %d, ‖Δβ‖ = %f, ‖Δτ‖ = %f, ‖ΔL‖ = %f, status = %s, time(s) = %f\n",
            run,
            norm(m.β  - βprev ),
            norm(m.τ  - τprev ),
            norm(m.Lγ - Lγprev),
            optstat,
            toc - tic)
    end
    m.isfitted[1] = true
    # refresh objective, gradient, and Hessian
    mul!(m.Σγ, m.Lγ, transpose(m.Lγ))
    nlsv_obj!(m, true, true, false)
    # sandwich estimator
    sandwich!(m)
    m
end

"""
    modelpar_to_optimpar!(par, m)

Translate model parameters in `m` to optimization variables in `par`.
"""
function modelpar_to_optimpar!(
    par :: AbstractVector{T},
    m   :: WSVarLmmModel{T}
    ) where {T<:BlasReal}
    q, l = m.q, m.l
    # τ
    copyto!(par, m.τ)
    # Lγ
    offset = l + 1
    @inbounds for j in 1:q, i in j:q
        par[offset] = i == j ? log(max(m.Lγ[i, j], floatmin(T))) : m.Lγ[i, j]
        offset += 1
    end
    par
end

"""
    optimpar_to_modelpar!(m, par)

Translate optimization variables in `par` to the model parameters in `m`.
"""
function optimpar_to_modelpar!(
    m   :: WSVarLmmModel{T},
    par :: Vector{T}
    ) where {T<:BlasReal}
    q, l = m.q, m.l
    # τ
    copyto!(m.τ, 1, par, 1, l)
    # Lγ
    fill!(m.Lγ, 0)
    offset = l + 1
    @inbounds for j in 1:q, i in j:q
        m.Lγ[i, j] = i == j ? exp(par[offset]) : par[offset]
        offset += 1
    end
    m
end

function MOI.initialize(
    m                  :: WSVarLmmModel,
    requested_features :: Vector{Symbol}
    )
    for feat in requested_features
        if !(feat in MOI.features_available(m))
            error("Unsupported feature $feat")
        end
    end
end

MOI.features_available(m::WSVarLmmModel) = [:Grad, :Hess]

function MOI.eval_objective(
    m   :: WSVarLmmModel{T},
    par :: AbstractVector{T}
    ) where {T<:BlasReal}
    optimpar_to_modelpar!(m, par)
    nlsv_obj!(m, false, false, false)
end

function MOI.eval_objective_gradient(
    m    :: WSVarLmmModel{T},
    grad :: Vector{T},
    par  :: AbstractVector{T}
    ) where {T<:BlasReal}
    q, l = m.q, m.l
    optimpar_to_modelpar!(m, par)
    obj = nlsv_obj!(m, true, false, false)
    # gradient wrt τ
    copyto!(grad, m.∇τ)
    # gradient wrt Lγ
    offset = l + 1
    @inbounds for j in 1:q, i in j:q # traverse lower triangular of Lγ[i, j]
        grad[offset] = m.∇Lγ[i, j]
        i == j && (grad[offset] *= m.Lγ[i, j])
        offset += 1
    end
end

function MOI.eval_constraint(
    m   :: WSVarLmmModel,
    g   :: Vector{T},
    par :: AbstractVector{T}
    ) where {T<:BlasReal}
    return nothing
end

function MOI.hessian_lagrangian_structure(m::WSVarLmmModel)
    # our Hessian is a dense matrix, work on the upper triangular part
    npar = m.l + ◺(m.q)
    arr1 = Vector{Int}(undef, ◺(npar))
    arr2 = Vector{Int}(undef, ◺(npar))
    idx  = 1
    @inbounds for j in 1:npar, i in 1:j
        arr1[idx] = i
        arr2[idx] = j
        idx      += 1
    end
    return zip(arr1, arr2)
end

function MOI.eval_hessian_lagrangian(
    m   :: WSVarLmmModel,
    H   :: AbstractVector{T},
    par :: AbstractVector{T},
    σ   :: T,
    μ   :: AbstractVector{T}
    ) where {T<:BlasReal}
    q, l = m.q, m.l
    # refresh obj, gradient, and hessian
    optimpar_to_modelpar!(m, par)
    nlsv_obj!(m, true, true, false)
    # Hττ
    idx = 1
    @inbounds for j in 1:l, i in 1:j
        H[idx] = m.Hττ[i, j]
        idx   += 1
    end
    j = 1 # index columns of HτLγ and HLγLγ
    @inbounds for j2 in 1:q, j1 in j2:q # traverse lower triangular of Lγ[j1, j2]
        # HτLγ
        for i in 1:l # i index rows of HτLγ
            H[idx] = m.HτLγ[i, j]
            j1 == j2 && (H[idx] *= m.Lγ[j1, j2])
            idx += 1
        end
        # HLγLγ
        i = 1 # index rows of HLγLγ
        for i2 in 1:q, i1 in i2:q # traverse lower triangular of Lγ[i1, i2]
            i > j && break # skip lower triangular of HLγLγ
            H[idx] = m.HLγLγ[i, j]
            # different diagonal entries of Lγ
            i1 == i2 && (H[idx] *= m.Lγ[i1, i2])
            j1 == j2 && (H[idx] *= m.Lγ[j1, j2])
            # same diagonal entry of Lγ
            i1 == i2 == j1 == j2 && (H[idx] += m.∇Lγ[j1, j2] * m.Lγ[j1, j2])
            idx += 1
            i   += 1
        end
        j += 1
    end
    lmul!(σ, H)
end


text:
mutable struct BinaryNode{T}
    left::Union{Nothing, T}
    right::Union{Nothing, T}

    BinaryNode{T}() where {T} = new{T}(nothing, nothing)
    BinaryNode(left::T, right::T) where {T} = new{T}(nothing, nothing)
    BinaryNode{T}(left::T, right::T) where {T} = new{T}(nothing, nothing)
end

mutable struct RectanglePacker{T}
    children::BinaryNode{RectanglePacker{T}}
    area::Rect2D{T}
end

left(a::RectanglePacker) = a.children.left
left(a::RectanglePacker{T}, r::RectanglePacker{T}) where {T} = (a.children.left = r)
right(a::RectanglePacker) = a.children.right
right(a::RectanglePacker{T}, r::RectanglePacker{T}) where {T} = (a.children.right = r)
RectanglePacker(area::Rect2D{T}) where {T} = RectanglePacker{T}(BinaryNode{RectanglePacker{T}}(), area)
isleaf(a::RectanglePacker) = (a.children.left) == nothing && (a.children.right == nothing)
# This is rather append, but it seems odd to use another function here.
# Maybe its a bad idea, to call it push regardless!?
function Base.push!(node::RectanglePacker{T}, areas::Vector{Rect2D{T}}) where T
    sort!(areas, by=GeometryBasics.norm ∘ widths)
    return RectanglePacker{T}[push!(node, area) for area in areas]
end

function Base.push!(node::RectanglePacker{T}, area::Rect2D{T}) where T
    if !isleaf(node)
        l = push!(left(node), area)
        l !== nothing && return l
        # if left does not have space, try right
        return push!(right(node), area)
    end
    newarea = RectanglePacker(area).area
    if all(widths(newarea) .<= widths(node.area))
        neww, newh = widths(newarea)
        xmin, ymin = minimum(node.area)
        xmax, ymax = maximum(node.area)
        w, h = widths(node.area)
        oax,oay,oaxw,oayh = xmin + neww, ymin, xmax, ymin + newh
        nax,nay,naxw,nayh = xmin, ymin + newh, xmax, ymax
        rax,ray,raxw,rayh = xmin, ymin, xmin + neww, ymin + newh
        left(node, RectanglePacker(Rect2D(oax, oay, oaxw - oax, oayh - oay)))
        right(node, RectanglePacker(Rect2D(nax, nay, naxw - nax, nayh - nay)))
        return RectanglePacker(Rect2D(rax, ray, raxw - rax, rayh - ray))
    end
    return nothing
end

ext: jl


text:
using Test
using Quon
using Yao

function check_inv(a, b, c)
    theta = [a,b,c] * im

    phi = yang_baxter_param_inv(yang_baxter_param(theta[1], theta[2], theta[3])...)
    phi = imag.([phi...])

    circ1 = chain(1, Rz(a), Rx(b), Rz(c))
    circ2 = chain(1, Rz(phi[1]), Rx(phi[2]), Rz(phi[3]))
    return operator_fidelity(circ1, circ2) ≈ 1
end
check_inv() = check_inv(rand(3)...)
check_inv()

@test check_inv()


text:
# ------------------------------------------------------------------------------------------#
# constants (from bionumbers)       units
# ------------------------------------------------------------------------------------------#
cell_diameter = 1.1                 # mum
number_of_rnapII = 4600            	# copies/cells
number_of_ribosome = 50000         	# copies/cells
mRNA_half_life_TF = 0.083           # hrs
protein_half_life = 70              # hrs
doubling_time_cell = 0.33           # hrs
max_translation_rate = 16.5         # aa/sec
max_transcription_rate = 60.0       # nt/sec
average_transcript_length = 1200   	# nt
average_protein_length = 400       	# aa
fraction_nucleus = 0.0             	# dimensionless
av_number = 6.02e23                 # number/mol
avg_gene_number = 2                 # number of copies of a gene
polysome_number = 4									# number of ribsomoses per transcript
# ------------------------------------------------------------------------------------------#
#
# ------------------------------------------------------------------------------------------#
# Calculate constants using bionumber values
# ------------------------------------------------------------------------------------------#
# Calculate the volume (convert to L)
V = ((1-fraction_nucleus)*(1/6)*(3.14159)*(cell_diameter)^3)*(1e-15)     # volume

# Calculate the rnapII_concentration and ribosome_concentration
rnapII_concentration = number_of_rnapII*(1/av_number)*(1/V)*1e9                   # nM
ribosome_concentration = number_of_ribosome*(1/av_number)*(1/V)*1e9               # nM

# degrdation rate constants -
degradation_constant_mRNA = -(1/mRNA_half_life_TF)*log(0.5)                       # hr^-1
degradation_constant_protein = -(1/protein_half_life)*log(0.5)                    # hr^-1

# kcats for transcription and translation -
kcat_transcription = max_transcription_rate*(3600/average_transcript_length)      # hr^-1
kcat_translation = polysome_number*max_translation_rate*(3600/average_protein_length)             # hr^-1

# Maximum specific growth rate -
maximum_specific_growth_rate = (1/doubling_time_cell)*log(2)                      # hr^-1

# What is the average gene concentration -
avg_gene_concentration = avg_gene_number*(1/av_number)*(1/V)*1e9                  # nM

# How fast do my cells die?
death_rate_constant = 0.05*maximum_specific_growth_rate                            # hr^-1

# Saturation constants for translation and trascription -
saturation_transcription = 4600*(1/av_number)*(1/V)*1e9                           # nM
saturation_translation = 150000*(1/av_number)*(1/V)*1e9                           # nM
# -------------------------------------------------------------------------------------------#


max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test

@testset "App" begin
    include("HealthHandler.jl")
end

max_stars_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_issues_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia

max_forks_repo_name: DevOpsDemoTF/DevOpsDemo-template-Julia


text:
## FitnessFun.m states the current design objective
# We want to minimize the score; rate by closeness to zero
Position=1:Brood;
Fitness=[GenerationScores,Position']
Ranked=sortrows(Fitness,1)
clear Position Fitness;

max_stars_repo_name: DesignEngrLab/ColorGraph

max_issues_repo_name: DesignEngrLab/ColorGraph

max_forks_repo_name: DesignEngrLab/ColorGraph


text:
using RoME
using Test

##

@testset "test DynPose2 and velocity..." begin

##

N = 100
fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

# initialize the first pose
IncrementalInference.doautoinit!(fg, [getVariable(fg,:x0);])

addVariable!(fg, :x1, DynPose2;  nanosecondtime=1000_000_000)

# conditional likelihood between Dynamic Point2
dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

pts = approxConv(fg, :x0x1f1, :x1)

# ensureAllInitialized!(fg)

tree, smt, hist = solveTree!(fg);

X1 = getVal(fg, :x1)

@test 0.9*N <= sum(abs.(X1[1,:] .- 10.0) .< 0.75)
@test 0.9*N <= sum(abs.(X1[2,:] .- 0.0) .< 0.75)
# @show TU.wrapRad.(X1[3,:])
@test 0.8*N <= sum(abs.(TU.wrapRad.(X1[3,:]) .- 0.0) .< 0.25)
# @warn "wrapRad issue, accepting 80% as good enough until issue JuliaRobotics/RoME.jl#90 is fixed."
@test 0.9*N <= sum(abs.(X1[4,:] .- 10.0) .< 0.5)
@test 0.9*N <= sum(abs.(X1[5,:] .- 0.0) .< 0.5)

# using RoMEPlotting
# # plotLocalProduct(fg, :x10, dims=[1;2])
# # plotSLAM2DPoses(fg)
# xx1 = marginal(getKDE(fg, :x1),[1;2;3])
# plotPose(Pose2(), [xx1])
# plotPose(fg, [:x1], levels=1, show=false)
#
# plotKDE(marginal(getKDE(fg, :x1),[4;5]), levels=5)

##

end

@testset "test distribution compare functions..." begin

##

mu = randn(6)
mv1 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
mv2 = MvNormal(deepcopy(mu), Matrix{Float64}(LinearAlgebra.I, 6,6))
mv3 = MvNormal(randn(6), Matrix{Float64}(LinearAlgebra.I, 6,6))
@test RoME.compareDensity(mv1, mv2)
@test !RoME.compareDensity(mv1, mv3)
@test !RoME.compareDensity(mv2, mv3)

##

end

@testset "test DynPose2 packing converters..." begin

##

pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([10.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

pp = convert(PackedDynPose2VelocityPrior, pp0)
ppu = convert(DynPose2VelocityPrior, pp)

@test RoME.compare(pp0, ppu)

dp2dp2 = VelPose2VelPose2(MvNormal([10.0;0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))

pp = convert(PackedVelPose2VelPose2, dp2dp2)
ppu = convert(VelPose2VelPose2, pp)

@test RoME.compare(dp2dp2, ppu)

##

end

@testset "test many DynPose2 chain stationary and 'pulled'..." begin

##

N = 100
fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
pp0 = DynPose2VelocityPrior(MvNormal(zeros(3), Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;], pp0)

sym = :x0
k = 0
for sy in Symbol[Symbol("x$i") for i in 1:10]

k+=1
addVariable!(fg, sy, DynPose2; nanosecondtime=1000_000_000*k)

# conditional likelihood between Dynamic Point2
dp2dp2 = VelPose2VelPose2(MvNormal([0.0;0;0], Matrix(Diagonal([1.0;0.1;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [sym;sy], dp2dp2)
sym =sy

end # for

x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.25
@test abs(x5[2]) < 1.25
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

ensureAllInitialized!(fg)

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 1.25
@test abs(x10[2]) < 1.25
@test abs(TU.wrapRad(x10[3])) < 0.4
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# drawGraph(fg, show=true)
# tree = buildTreeReset!(fg)
# drawTree(tree, show=true)

# using RoMEPlotting
# Gadfly.set_default_plot_size(35cm, 25cm)
# plotSLAM2DPoses(fg)
# plotPose(fg, [:x10])

##

@error ".useMsgLikelihoods = false required until IIF #1010 completed."
getSolverParams(fg).useMsgLikelihoods = false

# solve
smtasks = Task[]
tree, smt, hist = solveTree!(fg, smtasks=smtasks); #, recordcliqs=ls(fg));

##

x5 = getPPE(getVariable(fg, :x5)).suggested
# x5 = KDE.getKDEMean(getKDE(getVariable(fg, :x5)))

@test abs(x5[1]) < 1.5
@test abs(x5[2]) < 1.5
@test abs(TU.wrapRad(x5[3])) < 0.4
@test abs(x5[4]) < 0.5
@test abs(x5[5]) < 0.5

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test abs(x10[1]) < 2.75
@test abs(x10[2]) < 2.75
@test abs(TU.wrapRad(x10[3])) < 0.5
@test abs(x10[4]) < 0.5
@test abs(x10[5]) < 0.5

# pull the tail end out with position
pp10 = DynPose2VelocityPrior(MvNormal([10.0;0;0], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x10;], pp10)

fg2 = deepcopy(fg)

tree, mst, hist = solveTree!(fg) # N=N

x10 = KDE.getKDEMean(getKDE(getVariable(fg, :x10)))

@test 5.0 < x10[1]
@test abs(x10[2]) < 1.0
@test abs(TU.wrapRad(x10[3])) < 0.6
@test -0.1 < x10[4] < 1.0
@test abs(x10[5]) < 0.5

for sym in [Symbol("x$i") for i in 2:9]

XX = KDE.getKDEMean(getKDE(getVariable(fg, sym)))

@show sym, round.(XX,digits=5)
@test -2.0 < XX[1] < 10.0
@test abs(XX[2]) < 1.0
@test abs(TU.wrapRad(XX[3])) < 1.3
@test -0.5 < XX[4] < 2.0
@test abs(XX[5]) < 0.5

end

##

end

# using RoMEPlotting
# plotLocalProduct(fg, :x10, dims=[1;2])
# drawPoses(fg)
# plotPose(fg, [:x9],levels=1);
#
# plotPose(fg, [:x1;:x2;:x3;:x4;:x5;:x6;:x7;:x8;:x9;:x10],levels=1);

# savejld(fg) # tempfg.jld

@testset "test many DynPose2 sideways velocity..." begin

##

global N = 100
global fg = initfg()

# add first pose locations
addVariable!(fg, :x0, DynPose2; nanosecondtime=0)

# Prior factor as boundary condition
global pp0 = DynPose2VelocityPrior(MvNormal([0.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x0;], pp0)

addVariable!(fg, :x1, DynPose2; nanosecondtime=1000_000_000)

global pp0 = DynPose2VelocityPrior(MvNormal([1.0;0.0;pi/2], Matrix(Diagonal([0.01; 0.01; 0.001].^2))),
                            MvNormal([0.0;0], Matrix(Diagonal([0.5; 0.5].^2))))
addFactor!(fg, [:x1;], pp0)

# conditional likelihood between Dynamic Point2
global dp2dp2 = VelPose2VelPose2(MvNormal([0.0;-1.0;0], Matrix(Diagonal([0.01;0.01;0.001].^2))),
                          MvNormal([0.0;0], Matrix(Diagonal([0.1; 0.1].^2))))
addFactor!(fg, [:x0;:x1], dp2dp2)

getSolverParams(fg).N = N
solveTree!(fg)

# test for velocity in the body frame
global x0 = KDE.getKDEMean(getKDE(getVariable(fg, :x0)))

@test -0.4 < x0[1] < 2.0
@test abs(x0[2]) < 0.5
@test abs(x0[3] - pi/2) < 0.1
@test abs(x0[4]) < 0.4
@test -1.5 < x0[5] < -0.5

global x1 = KDE.getKDEMean(getKDE(getVariable(fg, :x1)))

@test -0.1 < x1[1] < 2.0
@test abs(x1[2]) < 0.5
@test abs(x1[3] - pi/2) < 0.1
@test abs(x1[4]) < 0.4
@test -1.5 < x1[5] < -0.5

##

end

#


text:
"""
direct product of two or more sets
"""
struct DirectProduct{T} <: AbstractProduct
    itr :: T
end
const DProd = DirectProduct

# Direct product INITIALIZERS

"""
return direct product of finite sets xs
"""
DirectProduct(xs...) = DirectProduct(product(xs...))

# from vector/tuple of finite sets
function DirectProduct(v :: Union{Vector{T}, NTuple{N, T}}) where T <: FinSet where N
    return DirectProduct(v...)
end

# from vector/tuple of Integers
function DirectProduct(v :: Union{Vector{T}, NTuple{N, T}}) where T <: Integer where N
    dom = map(Segment, v)
    return DirectProduct(dom)
end

# INTERFACE implementation
iterator(dp :: DirectProduct) = dp.itr
factors(dp :: DirectProduct) = dp.itr.iterators

max_stars_repo_name: kirtsar/DiscreteFunctions

max_issues_repo_name: kirtsar/DiscreteFunctions


text:
# # Lorenz equations
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/models/Lorenz.ipynb)
#
#md # !!! note "Overview"
#md #     System type: Continuous blackbox system\
#md #     State dimension: 3\
#md #     Application domain: Atmospheric convection
#
# ## Model description
#
# The model is a system of three ordinary differential equations now known as
# the Lorenz equations:
#
# ```math
#   \begin{array}{lcl}
#   \dfrac{dx}{dt} & = & \sigma (y-x), \\ \\
#   \dfrac{dy}{dt} & = & x(\rho - z) - y, \\ \\
#   \dfrac{dz}{dt} & = & xy-\beta z
#   \end{array}
# ```
#
# The equations relate the properties of a two-dimensional fluid layer uniformly
# warmed from below and cooled from above. In particular, the equations describe
# the rate of change of three quantities with respect to time:``x`` is
# proportional to the rate of convection, ``y`` to the horizontal temperature
# variation, and ``z`` to the vertical temperature variation. The constants
# ``\sigma`` , ``\rho`` , and ``\beta``  are system parameters proportional to
# the Prandtl number, Rayleigh number, and certain physical dimensions of the
# layer itself.

using ReachabilityAnalysis, Plots

@taylorize function lorenz!(dx, x, params, t)
    local σ = 10.0
    local β = 8.0 / 3.0
    local ρ = 28.0
    dx[1] = σ * (x[2] - x[1])
    dx[2] = x[1] * (ρ - x[3]) - x[2]
    dx[3] = x[1] * x[2] - β * x[3]
    return dx
end

# ## Reachability settings
#
# The initial values considered are ``X_0 \in [0.9, 1.1] \times [0, 0] \times [0, 0]``,
# for a time span of `10`.

X0 = Hyperrectangle(low=[0.9, 0.0, 0.0], high=[1.1, 0.0, 0.0])
prob = @ivp(x' = lorenz!(x), dim=3, x(0) ∈ X0);

# ## Results

# We compute the flowpipe using the TMJets algorithm with ``n_T=10`` and ``n_Q=2``.

alg = TMJets(abstol=1e-15, orderT=10, orderQ=2, maxsteps=50_000);

sol = solve(prob, T=10.0, alg=alg);

solz = overapproximate(sol, Zonotope);

#-

plot(solz, vars=(0, 1), xlab="t", ylab="x")

# It is apparent by inspection that variable $x(t)$ does not exceed 20.0 in the computed
# time span:

plot(solz(0.0 .. 1.5), vars=(0, 1), xlab="t", ylab="x", lw=0.0)
plot!(x -> 20.0, c=:red, xlims=(0.0, 1.5), lab="")

# We can prove that it is the case by computing the support function of the flowpipe
# along direction `[1.0, 0.0, 0.0]`:

ρ([1.0, 0.0, 0.0], solz(0.0 .. 1.5))

# In a similar fashion, we can compute extremal values of variable $y(t)$:

plot(solz, vars=(0, 2), xlab="t", ylab="y")

# Since we have computed overapproximations of the exact flowipe, the following
# quantities are a lower bound on the exact minimum (resp. an uppper bound on the
# exact maximum):

-ρ([0.0, -1.0, 0.0], solz)

ρ([0.0, 1.0, 0.0], solz)

#-

plot(solz, vars=(0, 3), xlab="t", ylab="z")

#-

# Below we plot the flowpipe projected on the `(1, 3)` plane.

plot(solz, vars=(1, 3))


text:
# Vectors defining a plane source: the plane normal and the direction of dipole current.

# PlaneSrc supports the polarization normal to the plane.  In that case, the source strength
# is interpreted as J∆n rather than sheet current density K.

# For the TE equations, the electric plane source has K = 2 and Kf = 2, whereas the magnetic
# plane source has K = 2 and Kf = 1.  This might be nonintuitive to the users, but it
# becomes obvious when the equations are written down: in the TE equations on the xy-plane,
# the electric current source has both the x- and y-components (so Kf = 2), whereas the
# magnetic current source has only the z-component (so Kf = 1).

export PlaneSrc

mutable struct PlaneSrc{K,Kf} <: Source{K,Kf}
    # Specify geometry first, current next (because the value of current is optional).
    n::SFloat{K}  # direction normal to plane; currently supports only Cartesian directions
    c::Float  # intercept (c means center)
    p::SFloat{Kf}  # polarization direction in field dimension
    J∆n::ComplexF # sheet current density (current per unit in-plane length normal to flow direction); complex in order to represent phase difference between plane sources
    isfield˔shp::Bool  # true if spaces where field and shapes are defined are orthogonal complement of each other; false if they are the same

    PlaneSrc{K,Kf}(n::AbsVecReal, c::Real, p::AbsVecReal, J∆n::Number, isfield˔shp::Bool) where {K,Kf} = new(normalize(n), c, normalize(p), J∆n, isfield˔shp)
end

function PlaneSrc(n::AbsVecReal, c::Real, p::AbsVecReal, J∆n::Number=1.0, isfield˔shp::Bool=isfield_ortho_shape(length(p),length(n)))
    count(n.!=0) == 1 || throw(ArgumentError("n = $n must be along Cartesian direction."))

    K = length(n)
    Kf = length(p)

    return PlaneSrc{K,Kf}(n, c, p, J∆n, isfield˔shp)
end

function add_src!(jKd::AbsArrNumber{K₊₁},  # (K+1)-dimensional array of Je (electric) or Jm (magnetic); first K dimensions specify location; last 1 dimension specify field component
                  gt₀::SVec{K,GridType},  # grid type of voxel corners; generated by ft2gt.(ft, boundft)
                  bounds::Tuple2{SFloat{K}},  # bounds[NEG][k] = boundary of domain at negative end in k-direction
                  l::Tuple2{NTuple{K,VecFloat}},  # l[PRIM][k] = primal vertex locations in k-direction
                  ∆l::Tuple2{NTuple{K,VecFloat}},  # ∆l[PRIM][k] = (∆l at primal vertices in w) == diff(l[DUAL][k] including ghost point)
                  isbloch::SBool{K},  # Bloch boundary conditions
                  src::PlaneSrc{K,Kf}  # plane source to add
                  ) where {K,K₊₁,Kf}
    nn = findfirst(src.n.==1)  # index of Cartesian direction corresponding to plane normal
    for nw = 1:Kf  # set Jw (w-component of J)
        gt_cmp = src.isfield˔shp ? gt₀ : gt_w(nw, gt₀)  # grid type of Jw

        lcmp = t_ind(l, gt_cmp)
        ∆lcmp = t_ind(∆l, gt_cmp)

        gtn = gt_cmp[nn]
        ln = lcmp[nn]
        ∆ln = ∆lcmp[nn]
        # Above, ln is the field locations in the plane normal direction.
        # - If src.isfield˔shp = true, then gt_cmp = gt₀ for all nw, so ln does not change
        # with nw.
        # - If src.isfield˔shp = false, then gt_cmp changes with nw, but ln changes only
        # when nw = nn, because lcmp[nn], which is determined by gt_cmp[nn], remains gt₀[nn]
        # for nw ≠ nn and changes to alter(gt₀[nn]) only for nw = nn.
        # Using the above observation, it would be possible to simplify the above code, but
        # I leave it the way it is in order to maintain the similarity of the implementation
        # between concrete Source types.

        ind, wt = distweights(src.c, gtn, t_ind(bounds,nn,nn), ln, ∆ln, isbloch[nn])

        Jw∆n = src.J∆n * src.p[nw]  # nw is @inbounds, because src.p is length-Kf and 1 ≤ nw ≤ Kf
        Jw = Jw∆n * wt  # Jw has two elements like wt

        # The code below mimics the implementation of slicedim.
        # For example, V[Base.setindex(axes(V), ir:ir, nr)...] means V[:,ir:ir,:] for r = y.
        # (axes(V) = (1:Nx,1:Ny,1:Nz) and setindex(t,v,i) does t[i] = v.  Therefore,
        # Base.setindex(axes(V), ir:ir, nr) creates (1:Nx, ir:ir, 1:Nz) for r = y.)
        # The use of ir:ir instead of ir is to support 1D arrays.  If V is 1D, then V[ir] is
        # a scalar and the dot equal on V[ir] fails.
        for k = 1:2  # 2 == length(ind)
            # Below, [1:K] is used to take the first K entries of the length-K₊₁ tuple.
            jKd[Base.setindex(axes(jKd), ind[k]:ind[k], nn)[1:K]..., nw] .+= Jw[k]
        end
    end

    return nothing
end


text:
module TestRuns

using Test
using RLEVectors

@testset begin

# numruns single
@test numruns([1,2,2,3,3]) == 3
@test numruns([1,1,2,2,3]) == 3

# numruns pair
@test numruns(["a","b","c","d","d","e","f","g"],[3,6,9,12,15,18,21,24]) == 7 # Run of 2 runvalues
@test numruns(["a","b","c","d","d","d","e","f"],[3,6,9,12,15,18,21,24]) == 6 # Run of 3 runvalues
@test numruns(["a","b","c","d","e","f","g"],[3,6,9,12,12,15,17]) == 6 # run of 2 runends
@test numruns(["a","b","c","d","e","f"],[3,6,12,12,12,15]) == 4 # run of 3 runends
@test numruns([1,2,3],[1,2,2]) == 2
@test_throws ArgumentError numruns([1,2,3],[3,2,4])

# one vector
@test ree( [1,1,2,2,3,3,4,4,4] ) == ([1,2,3,4],[2,4,6,9])
@test ree( [1,1,2,2,3,3,4,4,7] ) == ([1,2,3,4,7],[2,4,6,8,9])
@test ree( [0,1,2,2,3,3,4,4,7] ) ==([0,1,2,3,4,7],[1,2,4,6,8,9])

# Clean up an RLEVector
@test ree(["a","b","c","d","d","e","f","g"],[3,6,9,12,15,18,21,24]) == (["a","b","c","d","e","f","g"],[3,6,9,15,18,21,24]) # Run of 2 runvalues
@test ree(["a","b","c","d","d","d","e","f"],[3,6,9,12,15,18,21,24]) == (["a","b","c","d","e","f"],[3,6,9,18,21,24]) # Run of 3 runvalues
@test ree(["a","b","c","d","e","f","g"],[3,6,9,12,12,15,17]) == (["a","b","c","d","f","g"],[3,6,9,12,15,17]) # run of 2 runends
@test ree(["a","b","c","d","e","f"],[3,6,12,12,12,15]) == (["a","b","c","f"],[3,6,12,15]) # run of 3 runends
@test ree([1,2,3],[1,2,2]) == ([1,2],[1,2])

# two vectors
@test ree([9],[1]) == ([9],[1])
@test ree([4,3],[2,4]) == ([4,3],[2,4])
@test ree([4,3],[0,4]) == ([3],[4])

# scalar
@test ree(5) == ([5],[1])
@test ree([5]) == ([5],[1])

# inverse_ree
@test inverse_ree([4,5,6],[2,5,7]) == [4,4,5,5,5,6,6]
@test inverse_ree(Real[],Int64[]) == Real[]
@test_throws ArgumentError inverse_ree([],[1])

end # testset

end # module


text:
"Abstract type to define residue alphabet types."
abstract type ResidueAlphabet end

"""
This type defines the usual alphabet of the 20 natural residues and a gap character.

```jldoctest
julia> using MIToS.MSA

julia> GappedAlphabet()
GappedAlphabet of length 21. Residues : res"ARNDCQEGHILKMFPSTWYV-"

```
"""
struct GappedAlphabet <: ResidueAlphabet end

"""
This type defines the usual alphabet of the 20 natural residues, without the gap character.

```jldoctest
julia> using MIToS.MSA

julia> UngappedAlphabet()
UngappedAlphabet of length 20. Residues : res"ARNDCQEGHILKMFPSTWYV"

```
"""
struct UngappedAlphabet <: ResidueAlphabet end

"""
`ReducedAlphabet` allows the construction of reduced residue alphabets, where residues
inside parenthesis belong to the same group.

```jldoctest
julia> using MIToS.MSA

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> ab[Residue('K')]
2

```
"""
struct ReducedAlphabet <: ResidueAlphabet
    mapping::Array{Int,1} # Residue -> Int
    named::NamedArray{Int,1,Array{Int,1},Tuple{OrderedDict{String,Int}}}
    len::Int
end

# Creation
# --------

function ReducedAlphabet(str::AbstractString)
    N = Int(XAA)
    mapping = fill!(Array{Int}(undef, N), N)
    ingroup = false
    pos = 0
    group_names = fill!(Array{String}(undef, N), "")
    for char in str
        if char == '('
            ingroup = true
            pos += 1
            continue
        end
        if char == ')'
            ingroup = false
            continue
        end
        if !ingroup
            pos += 1
        end
        int_residue = Int(Residue(char))
        @assert int_residue != 22 "$char isn't valid for a residue alphabet." # N == 22
        group_names[pos] = string(group_names[pos], char)
        mapping[int_residue] = pos
    end
    named = NamedArray(collect(1:pos), (
        OrderedDict{String,Int}(group_names[i] => i for i in 1:pos),), ("Groups",))
    ReducedAlphabet(mapping, named, pos)
end

# Iteration Interface
# -------------------

Base.length(ab::UngappedAlphabet)  = 20
Base.length(ab::GappedAlphabet)    = 21
Base.length(ab::ReducedAlphabet) = ab.len

Base.eltype(::ResidueAlphabet) = Int

function Base.iterate(ab::ResidueAlphabet, state=0)
    if state ≥ length(ab)
        return nothing
    else
        next_state = state + 1
        return (next_state, next_state)
    end
end

# Show
# ----

function Base.show(io::IO, ab::ResidueAlphabet)
    print(io, typeof(ab), " of length ", length(ab), ". Residues : res\"",
          join(_to_char[1:length(ab)]), "\"")
end

function Base.show(io::IO, ab::ReducedAlphabet)
    groups = names(ab.named, 1)
    print(io, "ReducedAlphabet of length ", length(ab), " : \"")
    for i in ab
        chars = groups[i]
        if length(chars) == 1
            print(io, chars[1])
        elseif length(chars) > 1
            print(io, "(", join(chars),")")
        end
    end
    print(io, "\"")
end

# getindex
# --------

@inline function Base.getindex(ab::ResidueAlphabet, res::Int)::Int
    ifelse(res <= length(ab), res, 22)
end

@inline function Base.getindex(ab::ReducedAlphabet, res::Int)::Int
    @inbounds value = ab.mapping[res]
    value
end

@inline function Base.getindex(ab::ResidueAlphabet, res::Residue)::Int
    ab[Int(res)]
end

@inline function Base.getindex(ab::ResidueAlphabet, res::String)::Int
    @assert length(res) == 1 "The string with the residue should have only one character."
    i = Int(Residue(res[1]))
    ifelse(i <= length(ab), i, 22)
end

@inline function Base.getindex(ab::ReducedAlphabet, res::String)::Int
    @inbounds value = ab.named[res]
    value
end

# In Alphabet
# -----------

Base.in(res::Residue, alphabet::ResidueAlphabet) = Int(res) <= length(alphabet)
Base.in(res::Residue, alphabet::ReducedAlphabet) = alphabet[res] != 22

# Names
# -----

"""
It returns the name of each group. The name is a string with the one letter code of each
residue that belong to the group.

```jldoctest
julia> using MIToS.MSA

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> names(ab)
8-element Array{String,1}:
 "AILMV"
 "RHK"
 "NQST"
 "DE"
 "FWY"
 "C"
 "G"
 "P"

```
"""
Base.names(alphabet::ReducedAlphabet) = names(alphabet.named, 1)
Base.names(alphabet::ResidueAlphabet) = String[ string(Residue(i)) for i in alphabet ]

# Dict of names to indexes
# ------------------------

@inline _getdict(n::NamedArray{Int,1,Array{Int,1},Tuple{OrderedDict{String,Int}}}) = n.dicts[1]

const _UngappedAlphabet_Names = OrderedDict{String,Int}(string(Residue(i))=>i for i in 1:20)
const _GappedAlphabet_Names = OrderedDict{String,Int}(string(Residue(i))=>i for i in 1:21)

@inline getnamedict(alphabet::UngappedAlphabet) = _UngappedAlphabet_Names
@inline getnamedict(alphabet::GappedAlphabet) = _GappedAlphabet_Names

"""
It takes a `ResidueAlphabet` and returns a dictionary from group name to group position.

```jldoctest
julia> using MIToS.MSA

julia> ab = ReducedAlphabet("(AILMV)(RHK)(NQST)(DE)(FWY)CGP")
ReducedAlphabet of length 8 : "(AILMV)(RHK)(NQST)(DE)(FWY)CGP"

julia> getnamedict(ab)
OrderedCollections.OrderedDict{String,Int64} with 8 entries:
  "AILMV" => 1
  "RHK"   => 2
  "NQST"  => 3
  "DE"    => 4
  "FWY"   => 5
  "C"     => 6
  "G"     => 7
  "P"     => 8

```
"""
@inline getnamedict(alphabet::ReducedAlphabet) = _getdict(copy(alphabet.named))


text:
# Test creation of elastic tensors from parameters

using CIJ, Test

@testset "Synthesis" begin
    @testset "Iso" begin
        let vp=8000, vs=4500, λ=vp^2-2vs^2, μ=vs^2, K=vp^2-4/3*vs^2, G=μ
            @test_throws ArgumentError iso(vp=64e6, lam=23.5e6)
            # Vp, Vs
            c = iso(vp=vp, vs=vs)
            @test c isa EC{CIJ.DEFAULT_FLOAT}
            @test c[1,1] == c[2,2] == c[3,3] == 64_000_000
            @test c[4,4] == c[5,5] == c[6,6] == 20_250_000
            @test c[1,2] == c[1,3] == c[2,3] == 23_500_000
            # λ, μ
            @test iso(lam=λ, mu=μ) == c
            # K, G
            @test iso(K=K, G=G) == c
        end
    end
end


text:
struct Baggr
    fm
    s_row  # in-bag
    s_col
    s_oob  # out-of-bag
end

"""
    baggr(X, Y, weights = nothing; rep = 50,
        withr = false, rowsamp = .7, colsamp = 1, fun, kwargs...)
Bagging of regression models.
* `X` : X-data (n obs., p variables).
* `Y` : Y-data (n obs., q variables).
* `weights` : Weights of the observations.
* `rep` : Nb. of bagging repetitions.
* `withr`: Type of sampling of the observations
    (`true` => with replacement).
* `rowsamp` : Proportion of rows sampled in `X`
    at each repetition.
* `colsamp` : Proportion of columns sampled (without replacement) in `X`
    at each repetition.
* `fun` : Name of the function computing the model to bagg.
* `kwargs` : Optional named arguments to pass in 'fun`.

## References

Breiman, L., 1996. Bagging predictors. Mach Learn 24, 123–140.
https://doi.org/10.1007/BF00058655

Breiman, L., 2001. Random Forests. Machine Learning 45, 5–32.
https://doi.org/10.1023/A:1010933404324

Genuer, R., 2010. Forêts aléatoires : aspects théoriques,
sélection de variables et applications. PhD Thesis. Université Paris Sud - Paris XI.

Gey, S., 2002. Bornes de risque, détection de ruptures, boosting :
trois thèmes statistiques autour de CART en régression (These de doctorat).
Paris 11. http://www.theses.fr/2002PA112245
"""
function baggr(X, Y, weights = nothing, wcol = nothing; rep = 50,
        withr = false, rowsamp = .7, colsamp = 1, fun, kwargs...)
    X = ensure_mat(X)
    Y = ensure_mat(Y)
    n, p = size(X)
    q = size(Y, 2)
    fm = list(rep)
    nrow = Int64(round(rowsamp * n))
    ncol = max(Int64(round(colsamp * p)), 1)
    s_row = fill(1, (nrow, rep))        # (nrow, rep)
    s_col = similar(s_row, ncol, rep)   # (ncol, rep)
    s_oob = list(rep, Vector{Int64})
    srow = similar(s_row, nrow)
    scol = similar(s_row, ncol)
    w = similar(X, nrow)
    zncol = collect(1:ncol)
    zX = similar(X, nrow, ncol)
    zY = similar(Y, nrow, q)
    #Threads.@threads for i = 1:rep
    #@sync @distributed for i = 1:rep
    #@distributed for i = 1:rep
    @inbounds for i = 1:rep
        srow .= sample(1:n, nrow; replace = withr)
        s_oob[i] = findall(in(srow).(1:n) .== 0)
        if colsamp == 1
            scol .= zncol
        else
            if isnothing(wcol)
                scol .= sample(1:p, ncol; replace = false)
            else
                scol .= sample(1:p, StatsBase.weights(wcol), ncol; replace = false)
            end
        end
        zX .= X[srow, scol]
        zY .= Y[srow, :]
        if(isnothing(weights))
            fm[i] = fun(zX, zY; kwargs...)
        else
            w .= mweights(weights[srow])
            fm[i] = fun(zX, zY, w; kwargs...)
        end
        s_row[:, i] .= srow
        s_col[:, i] .= scol
    end
    Baggr(fm, s_row, s_col, s_oob)
end

function predict(object::Baggr, X)
    rep = length(object.fm)
    scol = vcol(object.s_col, 1)
    # @view is not accepted by XGBoost.predict
    # @view(X[:, scol])
    acc = predict(object.fm[1], X[:, scol]).pred
    @inbounds for i = 2:rep
        scol = vcol(object.s_col, i)
        acc .+= predict(object.fm[i], X[:, scol]).pred
    end
    pred = acc ./ rep
    (pred = pred,)
end


text:
struct NExp{S}
    base::Var
    exp::Union{Poly,Var}
end

function (x::NExp{S})(n::Int) where {S}
    exp = MultivariatePolynomials.subs(x.exp, mkvar(S)=>n)
    @assert isconstant(exp)
    ex = convert(Number, exp)
    x.base^ex
end

_exp_map = Dict{NExp,Var}()

function replacement_pair(exp::NExp)
    var = get(_exp_map, exp, mkvar(gensym_unhashed(:r)))
    _exp_map[exp] = var
    var=>exp
end

struct CFiniteExpr{S}
    poly::Poly
    subs::Dict{Var, NExp{S}}
end

function _merge(x::Dict{Var, NExp{S}}, y::Dict{Var, NExp{S}}) where {S}
    @assert isempty(intersect(keys(x), keys(y)))
    merge(x, y)
end

Base.:+(x::CFiniteExpr{S}, y::CFiniteExpr{S}) where {S} = CFiniteExpr{S}(x.poly + y.poly, merge(x.subs, y.subs))
Base.:-(x::CFiniteExpr{S}, y::CFiniteExpr{S}) where {S} = CFiniteExpr{S}(x.poly - y.poly, merge(x.subs, y.subs))
Base.:*(x::CFiniteExpr{S}, y::CFiniteExpr{S}) where {S} = CFiniteExpr{S}(x.poly * y.poly, merge(x.subs, y.subs))
Base.:^(x::CFiniteExpr{S}, y::Int) where {S} = CFiniteExpr{S}(x.poly^y, x.subs)

Base.:+(x, y) = Base.:+(promote(x, y)...)
Base.:-(x, y) = Base.:-(promote(x, y)...)
Base.:*(x, y) = Base.:*(promote(x, y)...)

Base.convert(::Type{CFiniteExpr{S}}, x::Number) where {S} = CFiniteExpr{S}(mkpoly(x), Dict{Var, NExp{S}}())

Base.promote_rule(::Type{CFiniteExpr{S}}, ::Type{<:Number}) where {S} = CFiniteExpr{S}
Base.promote_rule(::Type{<:Number}, ::Type{CFiniteExpr{S}}) where {S} = CFiniteExpr{S}

function (expr::CFiniteExpr{S})(n::Int) where {S}
    s = mkvar(S)
    newsubs = Dict(k=>v(n) for (k,v) in expr.subs)
    MultivariatePolynomials.subs(expr.poly, newsubs...)
end

function order(expr::CFiniteExpr)
    ms, us = destructterm(expr.poly, collect(keys(expr.subs)))
    @assert sum(m*u for (m,u) in zip(ms,us)) == expr.poly "Factorization bug"
    all(isone, ms) ? 1 : length(ms)
end

function constraints(expr::CFiniteExpr{S}; split_vars::Vector{Symbol}=Symbol[]) where {S}
    cs = ClauseSet()
    # qs = map(x->CFiniteExpr{S}(x, expr.subs), destructpoly(expr.poly, split_vars))
    qs = destructpoly([expr.poly], split_vars)
    for (i, q) in enumerate(qs)
        cfin = CFiniteExpr{S}(q, expr.subs)
        l = order(cfin)
        for i in 0:l-1
            ex = Meta.parse(string(cfin(i)))
            cs &= Constraint{EQ}(ex)
        end
    end
    cs
end

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

coeffs(p::Poly, v::Var) = [coefficient(p, v^i, [v]) for i in 0:maxdegree(p, v)]

destructpoly(p::Poly, var::Var) = coeffs(p, var)
destructpoly(p::Poly, var::Var, left::Var...) =
    reduce(union, map(x->destructpoly(x, left...), coeffs(p, var)))

function destructpoly(ps, vars)
    xvars = map(mkvar, vars)
    isempty(xvars) ? ps : reduce(union, map(x->destructpoly(x, xvars...), ps))
end

function destructterm(p::Poly, rs::Vector{<:Var})
    ms = Poly[]
    cs = Poly[]
    for term in terms(p)
        ds = [maxdegree(term, r) for r in rs]
        m = prod(r^d for (r,d) in zip(rs,ds))
        c = div(term, m)
        push!(ms, m)
        push!(cs, c)
    end
    factor(ms, cs)
end

function factor(ms::Vector{<:Poly}, us::Vector{<:Poly})
    map = Dict{Poly,Poly}()
    for (m,u) in zip(ms,us)
        if haskey(map, m)
            map[m] += u
        else
            map[m] = u
        end
    end
    keys(map), Base.values(map)
end


text:
function sort(M::AbstractMatrix{Int})
    X = deepcopy(M)
    for i in 1:size(M, 2)
        for k in 1:size(M, 1)
            for j in 1:(size(M, 1)-1)
                if(X[j, i] > X[j+1, i])
                    X[j, i],X[j+1, i] = X[j+1, i],X[j, i]
                end
            end
        end
    end
    return X
end

function sort!(M::AbstractMatrix{Int})
    for i in 1:size(M, 2)
        for j in 1:size(M[i], 1)
            for k in 1:size(M[i], 1)

            end
        end
    end
end

print(sort([[6, 3, 0] [7, 4, 1] [8, 5, 2]]))


text:
#the calling function is responsible for allocating the proper number of array slots to do this correctly.
function flatten!{F, i, o, tf}(v::Vector{F}, fcl::FullyConnectedLayer{F, i, o, tf}, offset::Integer = 0)

  #first slot in the biases.
  @inbounds for idx = 1:o
    v[idx + offset] = fcl.bias[idx]
  end

  #next slot in the data.
  @inbounds for idx = 1:(i * o)
    v[idx + offset + o] = fcl.transition[idx]
  end
end

function unflatten!{F, i, o, tf}(fcl::FullyConnectedLayer{F, i, o, tf}, v::Vector{F}, offset::Integer = 0)
  #reverse the process
  @inbounds for idx = 1:o
    fcl.bias[idx] = v[idx + offset]
  end
  #next slot in the data.
  @inbounds for idx = 1:(i * o)
    fcl.transition[idx] = v[idx + offset + o]
  end
end

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
struct ConvexTotalChunker{F}
    f::F
end

struct ConvexTotalSplitter{F}
    f::F
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::ConvexTotalChunker, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f, A, args...)

        ftr = CircularDeque{Tuple{Ti, Ti}}(n + 1)
        spl = zeros(Ti, n + 1)
        cst = fill(typemax(cost_type(f)), n + 1)
        cst[1] = zero(cost_type(f))
        f′(j, j′) = cst[j] + f(j, j′)
        chunk_convex!(cst, spl, f′, 1, n + 1, ftr)

        return unravel_chunks!(spl, n)
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ConvexTotalSplitter, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f, A, args...)

        if K == 1
            return SplitPartition{Ti}(1, [Ti(1), Ti(n + 1)])
        end

        ftr = CircularDeque{Tuple{Ti, Ti}}(n + 1)

        ptr = zeros(Ti, n + 1, K)
        cst = fill(typemax(cost_type(f)), n + 1, K)

        for j′ = 1:n + 1
            cst[j′, 1] = f(1, j′, 1)
            ptr[j′, 1] = 1
        end
        for k = 2:K
            f′(j, j′) = cst[j, k - 1] + f(j, j′, k)
            for j′ = 1:n + 1
                cst[j′, k] = f′(j′, j′)
                ptr[j′, k] = j′
            end
            chunk_convex!((@view cst[:, k]), (@view ptr[:, k]), f′, 1, n + 1, ftr)
        end
        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end

function chunk_convex!(cst, ptr, f::F, j₀, j′₁, ftr) where {F}
    @inbounds begin
        empty!(ftr)
        push!(ftr, (j₀, j′₁ + 1))
        for j′ = j₀ + 1:j′₁
            (j, h) = last(ftr)
            c = f(j, j′)
            c′ = f(j′ - 1, j′)
            if c <= c′
                if c <= cst[j′]
                    cst[j′] = c
                    ptr[j′] = j
                end
                if h == j′ + 1
                    pop!(ftr)
                end
            else
                if c′ <= cst[j′]
                    cst[j′] = c′
                    ptr[j′] = j′ - 1
                end
                while !isempty(ftr) && ((j, h) = last(ftr); (f(j′ - 1, h - 1) < f(j, h - 1)))
                    pop!(ftr)
                end
                if isempty(ftr)
                    push!(ftr, (j′ - 1, j′₁ + 1))
                else
                    (j, h) = last(ftr)
                    #=
                    h_ref = j′ + 1
                    while h_ref < h - 1 && f(j′ - 1, h_ref) < f(j, h_ref)
                        h_ref += 1
                    end
                    h = h_ref
                    =#

                    h_lo = j′ + 1
                    h_hi = h - 1
                    while h_lo <= h_hi
                        h = fld2(h_lo + h_hi)
                        if f(j′ - 1, h - 1) < f(j, h - 1)
                            h_lo = h + 1
                        else
                            h_hi = h - 1
                        end
                    end
                    h = h_hi

                    if j′ + 1 != h
                        push!(ftr, (j′ - 1, h))
                    end
                end
            end
        end
    end
end

#=
function chunk_convex!(cst, ptr, f, j₀, j′₁, ftr)
    @inbounds begin
        for a = j₀:j′₁
            for b = a:j′₁
                for c = b + 1:j′₁
                    for d = c:j′₁
                        @assert f(a, c) + f(b, d) >= f(a, d) + f(b, c)
                    end
                end
            end
        end
    end

    @inbounds begin
        for j′ = j₀ + 1:j′₁
            for j = j₀:j′-1
                if f(j, j′) <= cst[j′]
                    cst[j′] = f(j, j′)
                    ptr[j′] = j
                end
            end
        end
    end
end
=#

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::ConvexTotalChunker{<:ConstrainedCost}, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f.f, A, args...)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max

        ftr = CircularDeque{Tuple{Ti, Ti}}(2n + 1)
        σ_j = undefs(Ti, 2n + 1)
        σ_j′ = undefs(Ti, 2n + 1)
        σ_ptr = undefs(Ti, 2n + 1)
        σ_cst = undefs(cost_type(f), 2n + 1)

        spl = zeros(Ti, n + 1)
        cst = fill(typemax(cost_type(f)), n + 1)
        cst[1] = zero(cost_type(f))
        f′ = let cst=cst, f=f
            @inline f′(j, j′) = @inbounds cst[j] + f(j, j′)
        end
        chunk_convex_constrained!(cst, spl, f′, w, w_max, 1, n + 1, ftr, σ_j, σ_j′, σ_cst, σ_ptr)

        return unravel_chunks!(spl, n)
    end
end

function partition_stripe(A::SparseMatrixCSC{Tv, Ti}, K, method::ConvexTotalSplitter{<:ConstrainedCost}, args...) where {Tv, Ti}
    @inbounds begin
        (m, n) = size(A)

        f = oracle_stripe(RandomHint(), method.f.f, A, args...)
        w = oracle_stripe(StepHint(), method.f.w, A, args...)
        w_max = method.f.w_max

        ftr = CircularDeque{Tuple{Ti, Ti}}(2n + 1)
        σ_j = undefs(Ti, 2n + 1)
        σ_j′ = undefs(Ti, 2n + 1)
        σ_ptr = undefs(Ti, 2n + 1)
        σ_cst = undefs(extend(cost_type(f)), 2n + 1)

        (j′_lo, j′_hi) = column_constraints(A, K, w, w_max)

        if j′_hi[K] < n + 1
            spl = ones(Ti, K + 1)
            spl[end] = n + 1
            #TODO throw(ArgumentError("infeasible"))
            return SplitPartition(K, spl)
        end

        ptr = WindowConstrainedMatrix{Ti, Ti}(zero(Ti), n + 1, K, j′_lo, j′_hi)
        cst = WindowConstrainedMatrix{extend(cost_type(f)), Ti}(infinity(cost_type(f)), n + 1, K, j′_lo, j′_hi, ptr.pos)

        for j′ = j′_lo[1] : j′_hi[1]
            cst[j′, 1] = extend(f(1, j′, 1))
            ptr[j′, 1] = 1
        end

        for k = 2:K
            f′(j, j′) = cst[j, k - 1] + extend(f(j, j′, k))
            for j′ = j′_lo[k]:j′_hi[k]
                cst[j′, k] = f′(j′, j′)
                ptr[j′, k] = j′
            end
            chunk_convex_constrained!((@view cst[:, k]), (@view ptr[:, k]), f′, w, w_max, j′_lo[k - 1], j′_hi[k], ftr, σ_j, σ_j′, σ_cst, σ_ptr)
        end

        return unravel_splits(K, n, PermutedDimsArray(ptr, (2, 1)))
    end
end

function chunk_convex_constrained!(cst, ptr, f::F, w, w_max, J₀, J′₁, ftr, σ_j, σ_j′, σ_cst, σ_ptr) where {F}
    @inbounds begin
        j′₁ = J₀ + 1
        while j′₁ < J′₁ && w(J₀, j′₁ + 1) <= w_max
            j′₁ += 1
        end
        j₀ = J₀

        while true
            chunk_convex!(cst, ptr, f, j₀, j′₁, ftr)

            if j′₁ == J′₁
                break
            end
            j′ = j′₁
            I = 1
            for j = j₀ + 1 : j′₁ # j = j₀ is redundant, as is j′ = j′₁.
                if j′ > j′₁
                    σ_j′[I] = j′
                    I += 1
                    σ_j[I] = j
                end
                while j′ < J′₁ && w(j, j′ + 1) <= w_max
                    j′ += 1
                    #if j′ > j′₁
                        σ_j′[I] = j′
                        I += 1
                        σ_j[I] = j
                    #end
                end
            end
            I += 1

            for i = 2:I-1
                σ_cst[i] = typemax(eltype(cst))
            end
            #have i < i′
            #need if j = σ_j[I], then j′ < σ_j′[I] (can do this with i < i′)
            #also need j < j′ but this is satisfied for any i, i′ since rand(σ_j) < rand(σ_j′)

            f′ = let I=I, σ_j=σ_j, σ_j′=σ_j′, f=f
                @inline f′(i, i′) = @inbounds f(σ_j[I - i], σ_j′[I - i′])
            end
            chunk_convex!(σ_cst, σ_ptr, f′, 1, I - 1, ftr)

            for i′ = 2:I - 1
                cst[σ_j′[I - i′]] = σ_cst[i′]
                ptr[σ_j′[I - i′]] = σ_j[I - σ_ptr[i′]]
            end

            j₀ = j′₁
            j′₁ = σ_j′[I - 2]
        end
    end
end


text:
### Instructions for adding a new version of the Rust toolchain
#
# * update the `version` variable
# * To deploy the shard and automatically update your BinaryBuilderBase's
#   `Artifacts.toml`, use the `--deploy` flag to the `build_tarballs.jl` script.
#   You can build & deploy by running:
#
#      julia build_tarballs.jl --debug --verbose --deploy
#

using BinaryBuilderBase, BinaryBuilder, Pkg.Artifacts
using BinaryBuilderBase: map_rust_target

include("../common.jl")

# We first download Rustup, and use that to install rust
rustup_name = "RustStage1"
rustup_version = v"1.24.3"

# This is the version of the Rust toolchain we install
version = v"1.57.0"

sources = [
    # We'll use rustup to install rust
    FileSource("https://static.rust-lang.org/rustup/archive/$(rustup_version)/x86_64-unknown-linux-musl/rustup-init",
               "bdf022eb7cba403d0285bb62cbc47211f610caec24589a72af70e1e900663be9"),
]

# The first thing we're going to do is to install Rust for all targets into a single prefix
script = "version=$(version)\n" * raw"""
cd ${WORKSPACE}/srcdir

# We reinstall RustBase because that's easier than re-using it
export CARGO_HOME=${prefix}
export RUSTUP_HOME=${prefix}
chmod +x rustup-init
./rustup-init -y --no-modify-path --default-host=${rust_host} --default-toolchain ${version}

# Collection of all rust targets we will download toolchains for:
RUST_TARGETS=(
    aarch64-apple-darwin
    aarch64-unknown-linux-gnu
    aarch64-unknown-linux-musl
    arm-unknown-linux-gnueabihf
    arm-unknown-linux-musleabihf
    armv7-unknown-linux-gnueabihf
    armv7-unknown-linux-musleabihf
    i686-unknown-linux-gnu
    i686-unknown-linux-musl
    x86_64-unknown-linux-gnu
    x86_64-unknown-linux-musl
    powerpc64le-unknown-linux-gnu
    i686-pc-windows-gnu
    x86_64-pc-windows-gnu
    x86_64-apple-darwin
    x86_64-unknown-freebsd
)

for rust_target in "${RUST_TARGETS[@]}"; do
    # Install our target-specific stuffs for the toolchain we're requesting
    ${CARGO_HOME}/bin/rustup target add --toolchain ${version} ${rust_target}
done

# We're going to bundle cargo-edit since it's a useful dep
export OPENSSL_STATIC=yes
${CARGO_HOME}/bin/cargo install cargo-edit
"""

# We assemble this giant tarball, then will split it up immediately after this:
platforms = [Platform("x86_64", "linux"; libc="musl")]
products = [
    ExecutableProduct("cargo", :cargo),
]
dependencies = [
    Dependency("OpenSSL_jll"),
]
ndARGS = filter(a -> !occursin("--deploy", a), ARGS)
build_info = build_tarballs(ndARGS, rustup_name, rustup_version, sources, script, platforms, products, dependencies; skip_audit=true)

# We don't actually need the .tar.gz it creates, so delete that to save space
rm(joinpath("products", first(values(build_info))[1]))

# Take the hash of the unpacked Rustup artifact, then split it into a bunch of smaller ones
mega_rust_path = artifact_path(first(values(build_info))[3])
rust_host = Platform("x86_64", "linux"; libc="musl")
rust_host_triplet = map_rust_target(rust_host)

for target_platform in supported_platforms(; experimental=true)
    rust_target_triplet = map_rust_target(target_platform)
    @info("Generating artifacts for $(rust_target_triplet)...")
    unpacked_hash = create_artifact() do dir
        srcpath = joinpath(mega_rust_path, "toolchains", "$(version)-$(rust_host_triplet)", "lib", "rustlib", rust_target_triplet)
        dstpath = joinpath(dir, "toolchains", "$(version)-$(rust_host_triplet)", "lib", "rustlib")
        mkpath(dstpath)
        cp(srcpath, joinpath(dstpath, rust_target_triplet))

        # Our mingw rust shards need to have their crt2.o updated
        # https://github.com/rust-lang/rust/issues/48272#issuecomment-429596397
        if Sys.iswindows(target_platform)
            # Find the corresponding mingw toolchain within a GCC shard
            all_cs = BinaryBuilderBase.all_compiler_shards()
            cs = first(filter(cs -> cs.name == "GCCBootstrap" && cs.target == target_platform && cs.archive_type == :unpacked, all_cs))
            cs_path = BinaryBuilderBase.mount(cs, "")

            # Locate our GCC's mingw's crt2.o
            crt_src = joinpath(cs_path, triplet(target_platform), "sys-root", "lib", "crt2.o")
            # Overwrite the one that rust ships with
            crt_dst = joinpath(dstpath, rust_target_triplet, "lib", "crt2.o")
            cp(crt_src, crt_dst; force=true)
        end
    end
    squashfs_hash = unpacked_to_squashfs(unpacked_hash, "RustToolchain", version; platform=rust_host, target=target_platform)

    # Upload them both to GH releases on Yggdrasil
    upload_and_insert_shards("JuliaPackaging/Yggdrasil", "RustToolchain", version, unpacked_hash, squashfs_hash, rust_host; target=target_platform)
end

# Finally, we do RustBase:
unpacked_hash = create_artifact() do dir
    cp(joinpath(mega_rust_path, "bin"), joinpath(dir, "bin"))
    cp(joinpath(mega_rust_path, "toolchains"), joinpath(dir, "toolchains"))
    rm(joinpath(dir, "toolchains", "$(version)-$(rust_host_triplet)", "share"); recursive=true)
    rm(joinpath(dir, "toolchains", "$(version)-$(rust_host_triplet)", "etc"); recursive=true)
    for rust_target_triplet in map_rust_target.(supported_platforms(; experimental=true))
        rm(joinpath(dir, "toolchains", "$(version)-$(rust_host_triplet)", "lib", "rustlib", rust_target_triplet); recursive=true)
    end
end

squashfs_hash = unpacked_to_squashfs(unpacked_hash, "RustBase", version; platform=rust_host)
upload_and_insert_shards("JuliaPackaging/Yggdrasil", "RustBase", version, unpacked_hash, squashfs_hash, rust_host)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
@testset "Core" begin

  @test sizeof(Jib.Client.Core.HEADTYPE) == 4

  @test Jib.Client.Core.MAX_LEN < typemax(Jib.Client.Core.HEADTYPE)

  @test Jib.Client.Core.isascii([0x45, 0x79])

  @test !Jib.Client.Core.isascii([0x45, 0x80])

end


text:
function orthopoly_evaluate_all(T, x)
    ##  Evaluate 0th through nth orthogonal polynomial at x
    #   Parameters
    #     T: Tridiagonal matrix encoding the three term recurrence
    #        size is (n+1)x(n+1)
    #     x: evaluation points
    #
    #  Returns
    #    phi: Polynomials 0 through n at x
    #         Method:  Spectrally solve (T-XI)u = constant * (last column of I)
    #         WARNING: potential for dividing by 0
    #         1st row of phi is 1 -- for w not a probability measure
    #                divide by 1/sqrt(c), where c is the total weight
    #                of the real line
    n = size(T,1)
    Lambda, H = eig(T)
    Hn = H[end, :]'
    phi = Array(Float64, length(Lambda), length(x))
    for m = 1:length(x)
        xi = x[m]
        # generate the u vector as in (2.5)
        u = H * (Hn./(Lambda - xi))
        phi[:,m] = u/u[1]
    end
    phi
end

function orthopoly_evaluate_n(T, x)
    ##  Evaluate nth orthogonal polynomial at x
    #   Parameters
    #     T: Tridiagonal matrix encoding the three term recurrence
    #        size is (n+1)x(n+1)
    #     x: evaluation points
    #
    #  Returns
    #    phi: nth at x
    #         normalized so that 0'th is 1 (w has unit weight)
    #         Method: characteristic polynomial of n x n T
    #                 with leading coefficient taken from the
    #                 product of the superdiagonal of the
    #                 (n+1)x(n+1) T

    n = size(T,1)-1
    ## compute the characteristic polynomial of T[1:n,1:n]
    Lambda = eigvals(T[1:n, 1:n])
    phi = ones(size(x))
    for i = 1:n
        phi = phi.*(x - Lambda[i])
    end
    ## normalization
    if n > 0
        phi = phi / prod(diag(T,1))
    end
    phi
end


text:
using Dash, DashHtmlComponents, DashCoreComponents, DataFrames, CSV

states = ["Baden-Wuerttemberg",
            "Bayern",
            "Berlin",
            "Brandenburg",
            "Bremen",
            "Hamburg",
            "Hessen",
            "Mecklenburg-Vorpommern",
            "Niedersachsen",
            "Nordrhein-Westfalen",
            "Rheinland-Pfalz",
            "Saarland",
            "Sachsen-Anhalt",
            "Sachsen",
            "Schleswig-Holstein",
            "Thueringen"]
country = "Germany"
dropdown_opts = [Dict("label"=>state, "value"=>state) for state in [country; states]]

# when running on heroku
download = true
github_url = "https://raw.githubusercontent.com/timueh/COVID-19/master/example/"

# when running locally
# download = false
# github_url = "../example/"

function get_data(path_to_file::String, ylabel::String, title::String)
    df = DataFrame(CSV.File(path_to_file))
    index, values = names(df)[1], names(df)[2:end]

    data = [ Dict( :x => df[!, index],
                :y => df[!, value],
                :type => "scatter",
                :name => value, )
                for value in values ]
    layout = Dict("yaxis" => Dict("title" => ylabel),
                  "title" => title)

    data, layout
end

function get_R_data_germany()
    file_ger = "results-R-reported-Germany.csv"
    data, layout = get_data(file_ger, "R", "")
    data_proj = data[end]
    data_proj[:name] = "Germany"
    data_proj
end

function download_data(github_url, suffix, states; download=true::Bool)
    values = ["R", "N"]
    for state in states
        for value in values
            for suff in suffix
                file_name = "results-"*value*"-"*suff*"-"*state*".csv"
                url = github_url*file_name
                display(url)
                if download
                    run(`curl $url --output $file_name`)
                else
                    run(`cp $url $file_name`)
                end
            end
        end
    end
end

download_state_data(url; download=true) = download_data(url, ["reported"], states, download=download)
download_germany_data(url; download=true) = download_data(url, ["reported", "nowcasting"], [country], download=download)

download_germany_data(github_url; download=download)
download_state_data(github_url; download=download)

## APP
app = dash(external_stylesheets = ["https://codepen.io/chriddyp/pen/bWLwgP.css"])

app.layout = html_div() do
    dcc_markdown("### Reproduction number `R` for COVID-19 pandemic in Germany"),
    html_h1(id="header"),
    dcc_dropdown(
            id = "location-source",
            options= dropdown_opts,
            value="Germany",
            ),
    dcc_dropdown(id="data-source", value="RKI-reported"),
    dcc_graph(id = "R-values"),
    dcc_graph(id = "N-values"),
    html_table(children = [
                    html_tr(children = [
                        html_th("Method"),
                        html_th("Explanation")
                    ]),
                    html_tr(children = [
                        html_td("RKI Nowcast 4 days"),
                        html_td(children = [dcc_markdown("[Estimator long used by Robert Koch Institut (RKI)](https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Projekte_RKI/R-Wert-Erlaeuterung.pdf?__blob=publicationFile); effectively a 4-day moving average.")]),
                        ]),
                    html_tr(children = [
                        html_td("RKI Nowcast 7 days"),
                        html_td(children = [dcc_markdown("[Estimator used by Robert Koch Institut (RKI)](https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Projekte_RKI/R-Wert-Erlaeuterung.pdf?__blob=publicationFile); effectively a 7-day moving average, taking 5 days from the past and one day from the future.")]),
                        ]),
                    html_tr(children = [
                        html_td("Projected 7 days"),
                        html_td(children = [dcc_markdown("[Acausal estimator that accounts for three days of the past, the current day, and three days of the future](https://www.medrxiv.org/content/10.1101/2020.11.27.20238618v1); future values are based on the respective values from the previous week.")]),
                        ]),
        ]
        ),
        html_table(children = [
            html_tr(children = [
                html_th("Data source"),
                html_th("Explanation")
            ]),
            html_tr(children = [
                html_td("RKI Reported cases"),
                html_td(children=[dcc_markdown("[Full data set of reported cases](https://npgeo-corona-npgeo-de.hub.arcgis.com/datasets/dd4580c810204019a7b8eb3e0b329dd6_0), curated and updated daily by the Robert Koch Institut (RKI)")]),
                ]),
            html_tr(children = [
                html_td("RKI Nowcasting"),
                html_td(children = [dcc_markdown("[Nowcasting data](https://www.rki.de/DE/Content/InfAZ/N/Neuartiges_Coronavirus/Projekte_RKI/Nowcasting.html;jsessionid=C42873168F44ED8B13EA88FDCEF2DE7A.internet062?nn=13490888) on total number of cases provided by Robert Koch Institut (RKI)")]),
                ]),
            ]
            )
end

callback!(app, [Output("R-values", "figure"), Output("N-values", "figure"), Output("header", "children")], [Input("location-source", "value"), Input("data-source", "value")]) do input_location, input_source

    suffix = "reported"
    if input_location == country && input_source == "RKI-nowcasting"
        suffix = "nowcasting"
    end

    file_R = "results-R-"*suffix*"-"*input_location*".csv"
    file_N = "results-N-"*suffix*"-"*input_location*".csv"

    R, l_R = get_data(file_R, "R", input_location)
    R_value = round(R[3][:y][end], digits=2)
    N, l_N = get_data(file_N, "Cases", "")

    if input_location != country
        # add R values for Germany to compare
        push!(R, get_R_data_germany())
    end

    Dict(:data=>R, :layout=>l_R), Dict(:data=>N, :layout=>l_N), "R = $R_value"
end

callback!(app, Output("data-source", "options"), Input("location-source", "value")) do input_value
    options = [Dict("label" => "Robert Koch Institut (RKI) - reported cases", "value" => "RKI-reported")]
    if input_value == country
        push!(options, Dict("label" => "Robert Koch Institut (RKI) - nowcasting cases", "value" => "RKI-nowcasting"))
    end
    options
end

run_server(app, "0.0.0.0", 8080)


text:
@testset "End of Week Crossover" begin
    t = 0
    @test t == (PositionVelocityTime.check_crossover(deepcopy(t)))

    t = 350000
    @test t - 604800 == PositionVelocityTime.check_crossover(deepcopy(t))

    t = -350000
    @test t + 604800 == PositionVelocityTime.check_crossover(deepcopy(t))

end

correct_raw_times =  [
    208863.71771887812,
    208863.72169284965,
    208863.72290867585,
    208863.7173908643,
    208863.7304948417]

@testset "Calc uncorrected Time" begin
    for i in 1:length(satellite_states)
        out = PositionVelocityTime.calc_uncorrected_time(satellite_states[i])
        @test out == correct_raw_times[i]
    end
end


text:
@inline function Flux(CellData,s,z,Def)
    """
    Flux Transfer Function
    # Add Ins and Outs
        # Cell Data
        # Frequency Vector
        # Discretisation Locations
        # Electrode Definition
    """

 if Def == "Pos"
    Electrode = CellData.Pos #Electrode Length
 else
    Electrode = CellData.Neg #Electrode Length
 end

CC_A = CellData.Const.CC_A   # Current-collector area [m^2]
as = 3*Electrode.ϵ_s/Electrode.Rs # Specific interfacial surf. area
κ_eff = CellData.Const.κ*Electrode.ϵ_e^Electrode.κ_brug #Effective Electrolyte Conductivity
σ_eff = Electrode.σ*Electrode.ϵ_s^Electrode.σ_brug #Effective Electrode Conductivity

#Defining SOC
θ = CellData.Const.SOC * (Electrode.θ_100-Electrode.θ_0) + Electrode.θ_0

#Beta's
β = @. Electrode.Rs*sqrt(s/Electrode.Ds )

#Prepare for j0
ce0 = CellData.Const.ce0
cs_max = Electrode.cs_max
cs0 = cs_max * θ
α = Electrode.α

#Current Flux Density
if CellData.Const.CellTyp == "Doyle_94"
    κ = Electrode.k_norm/Electrode.cs_max/ce0^(1-α)
    j0 = κ*(ce0*(cs_max-cs0))^(1-α)*cs0^α
else
   j0 = Electrode.k_norm*(ce0*(cs0/cs_max*(1-cs0/cs_max)))^(1-α)
end

#Resistances
Rtot = R*CellData.Const.T /(j0*F^2) + Electrode.RFilm

#∂Uocp_Def
∂Uocp_elc = CellData.Const.∂Uocp(Def,θ)/cs_max

#Condensing Variable
ν = @. Electrode.L*sqrt((as/σ_eff+as/κ_eff)/(Rtot+∂Uocp_elc*(Electrode.Rs/(F*Electrode.Ds))*(tanh(β)/(tanh(β)-β))))
ν_∞ = @. Electrode.L*sqrt(as*((1/κ_eff)+(1/σ_eff))/(Rtot))

#Transfer Function
j_tf = @. ν*(σ_eff*cosh(ν*z)+κ_eff*cosh(ν*(z-1)))/(as*F*Electrode.L*CC_A*(κ_eff+σ_eff)*sinh(ν))
D = @. ν_∞*(σ_eff*cosh(ν_∞*z)+κ_eff*cosh(ν_∞*(z-1)))/(as*F*Electrode.L*CC_A*(κ_eff+σ_eff)*sinh(ν_∞))
D_term = "@. $ν_∞*($σ_eff*cosh($ν_∞*$z)+$κ_eff*cosh($ν_∞*($z-1)))/($as*$F*$(Electrode.L)*$CC_A*($κ_eff+$σ_eff)*sinh($ν_∞))"
zero_tf =ones(size(z,1))*1/(CellData.Const.CC_A*as*F*Electrode.L)
j_tf[:,findall(s.==0)] .= zero_tf[:,findall(s.==0)]
res0 = zeros(length(z))

if Def == "Pos"
    j_tf = -j_tf
    D = -D
    D_term = "@. -$ν_∞*($σ_eff*cosh($ν_∞*$z)+$κ_eff*cosh($ν_∞*($z-1)))/($as*$F*$(Electrode.L)*$CC_A*($κ_eff+$σ_eff)*sinh($ν_∞))"
    if Debug == 1
        println("D:Flux:Pos",D)
        println("z:Flux:Pos",z)
        println("ν_∞:Flux:Pos",ν_∞)
        println("j0:Flux:Pos",j0)
        println("Rtot:Flux:Pos",Rtot)
        println("σ_eff:Flux:Pos",σ_eff)
        println("ν_∞:Flux:Pos",ν_∞)
        println("L:Flux:Pos",Electrode.L)
    end
else
    if Debug == 1
        println("D:Flux:Neg",D)
        println("z:Flux:Neg",z)
        println("ν_∞:Flux:Neg",ν_∞)
        println("D:Flux:Neg",D)
    end
end

return j_tf, D, res0, D_term

end

max_stars_repo_name: BradyPlanden/RealisationAlgorithms

max_issues_repo_name: BradyPlanden/RealisationAlgorithms

max_forks_repo_name: BradyPlanden/RealisationAlgorithms


text:
immutable Binomial <: DiscreteUnivariateDistribution
    size::Int
    prob::Float64
    function Binomial(n::Real, p::Real)
        n >= zero(n) || error("size must be positive")
        zero(p) <= p <= one(p) || error("prob must be in [0, 1]")
        new(int(n), float64(p))
    end
end

Binomial(size::Integer) = Binomial(size, 0.5)
Binomial() = Binomial(1, 0.5)

minimum(d::Binomial) = 0
maximum(d::Binomial) = d.size

@_jl_dist_2p Binomial binom

function _probs(d::Binomial, f::Int, l::Int)
    n = d.size
    p = d.prob
    b = f - 1
    r = Array(Float64, l - b)
    r[1] = v = pdf(d, f)
    if l > f
        c = p / (1.0 - p)
        for k = f+1:l
            v *= ((n - k + 1) / k * c)
            r[k-b] = v
        end
    end
    return r
end

probs(d::Binomial) = _probs(d, 0, d.size)

function probs(d::Binomial, rgn::UnitRange)
    f, l = rgn[1], rgn[end]
    0 <= f <= l <= d.size || throw(BoundsError())
    _probs(d, f, l)
end

function entropy(d::Binomial; approx::Bool=false)
    n = d.size
    p1 = d.prob

    (p1 == 0.0 || p1 == 1.0 || n == 0) && return 0.0
    p0 = 1.0 - p1

    if approx
        return 0.5 * (log(2.0pi * n * p0 * p1) + 1.0)
    else
        lg = log(p1 / p0)
        lp = n * log(p0)
        s = exp(lp) * lp
        for k = 1:n
           lp += log((n - k + 1) / k) + lg
           s += exp(lp) * lp
        end
        return -s
    end
end

kurtosis(d::Binomial) = (1.0 - 6.0 * d.prob * (1.0 - d.prob)) / var(d)

mean(d::Binomial) = d.size * d.prob

var(d::Binomial) = d.size * d.prob * (1.0 - d.prob)

skewness(d::Binomial) = (1.0 - 2.0 * d.prob) / std(d)

median(d::Binomial) = iround(d.size * d.prob)

# TODO: May need to subtract 1 sometimes
# two modes possible e.g. size odd, p = 0.5
mode(d::Binomial) = d.size > 0 ? iround((d.size + 1.0) * d.prob) : 0

function mgf(d::Binomial, t::Real)
    p = d.prob
    (1.0 - p + p * exp(t))^d.size
end

function cf(d::Binomial, t::Real)
    p = d.prob
    (1.0 - p + p * exp(im * t))^d.size
end

modes(d::Binomial) = iround([d.size * d.prob])

### handling support

isupperbounded(d::Union(Binomial, Type{Binomial})) = true
islowerbounded(d::Union(Binomial, Type{Binomial})) = true
isbounded(d::Union(Binomial, Type{Binomial})) = true

minimum(d::Union(Binomial, Type{Binomial})) = 0
maximum(d::Binomial) = d.size
support(d::Binomial) = 0:d.size

insupport(d::Binomial, x::Real) = isinteger(x) && 0 <= x <= d.size

## Fit model

immutable BinomialStats <: SufficientStats
    ns::Float64   # the total number of successes
    ne::Float64   # the number of experiments
    n::Int        # the number of trials in each experiment

    function BinomialStats(ns::Real, ne::Real, n::Integer)
        new(float64(ns), float64(ne), int(n))
    end
end

function suffstats{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T})
    ns = zero(T)
    for xi in x
        0 <= xi <= n || throw(DomainError())
        ns += xi
    end
    BinomialStats(ns, length(x), n)
end

function suffstats{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64})
    ns = 0.
    ne = 0.
    for i = 1:length(x)
        @inbounds xi = x[i]
        0 <= xi <= n || throw(DomainError())
        @inbounds wi = w[i]
        ns += xi * wi
        ne += wi
    end
    BinomialStats(ns, ne, n)
end

suffstats{T<:Integer}(::Type{Binomial}, data::(Int, Array{T})) = suffstats(Binomial, data...)
suffstats{T<:Integer}(::Type{Binomial}, data::(Int, Array{T}), w::Array{Float64}) = suffstats(Binomial, data..., w)

fit_mle(::Type{Binomial}, ss::BinomialStats) = Binomial(ss.n, ss.ns / (ss.ne * ss.n))

fit_mle{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T}) = fit_mle(Binomial, suffstats(Binomial, n, x))
fit_mle{T<:Integer}(::Type{Binomial}, n::Integer, x::Array{T}, w::Array{Float64}) = fit_mle(Binomial, suffstats(Binomial, n, x, w))
fit_mle{T<:Integer}(::Type{Binomial}, data::(Int, Array{T})) = fit_mle(Binomial, suffstats(Binomial, data))
fit_mle{T<:Integer}(::Type{Binomial}, data::(Int, Array{T}), w::Array{Float64}) = fit_mle(Binomial, suffstats(Binomial, data, w))


text:
function create_initial_tour(data, graph)
    """Create an initial tour with a very high cost and wich cover all
    required arcs/edges."""
    initial_tour = Dict()
    initial_tour["depot"] = Nothing
    initial_tour["reduced_cost"] = 10^3*(sum(graph["req_arcs"][:serv_cost]) +
        sum(graph["req_arcs"][:trav_cost]) + sum(graph["req_edges"][:serv_cost]) +
        sum(graph["req_edges"][:trav_cost]))
    initial_tour["cost"] = initial_tour["reduced_cost"]
    initial_tour["X"] = Dict(
        "req_arcs" => [1 for arc in 1:data[1,:req_arcs]],
        "req_edges" => [1 for arc in 1:data[1,:req_edges]],
        "req_edges_swap" => [1 for arc in 1:data[1,:req_edges]]
        )
    return(initial_tour)
end

function cut_branch(node, global_upper_bound, data, verbose)
    """Cut the branch if the current node can't improve the upper bound."""
    if (node["lower_bound"]>global_upper_bound)
        if verbose; println("Cut branch because of bound"); end
    elseif (isnan(node["lower_bound"]))
        if verbose; println("Cut branch because of NaN"); end
    end
    return((node["lower_bound"]>global_upper_bound) | (isnan(node["lower_bound"])))
end

function create_child_node(parent, branch_index)
    """From a parent node and a branch, create two child nodes."""

    even_node = Dict(
        "type" => "even",
        "branch" => push!(copy(parent["branch"]), branch_index),
        "tours" => copy(parent["tours"]),
        "dual_bound" => 0,
        "alpha" => [],
    )

    odd_node = Dict(
        "type" => "odd",
        "branch" => push!(copy(parent["branch"]), branch_index),
        "tours" => copy(parent["tours"]),
        "dual_bound" => 0,
        "alpha" => [],
    )

    return(even_node, odd_node)
end

function possible_branch(node, depot, arc_type, arc, branch_pool, data)
    """Checks if the branch doesn't already exists for the node and if it
    doesn't make the problem infeable by setting to 0 the number of passages
    in a required arc/edge."""

    if [depot, arc_type, arc] in [branch_pool[i] for i in node["branch"]]
        return(false)
    end

    depot_avoided_for_arc = 0
    for branch in [branch_pool[i] for i in node["branch"]]
        if (branch[2]==arc_type) & (branch[3]==arc)
            depot_avoided_for_arc += 1
        end
    end
    if depot_avoided_for_arc >= (length(data[1, :depots])-1)
        return(false)
    else
        return(true)
    end
end

function create_branch(node, tour_pool, branch_pool, data)
    """Create a branch from node depending of the coefficients of the
    master problem. Also checking if the branch hasen't already been created."""

    # Constructing delta
    delta = Dict()

    for depot in data[1, :depots]
        delta[depot] = Dict("req_arcs"=>[], "req_edges"=>[], "req_edges_swap"=>[])

        for arc in 1:data[1,:req_arcs]
            if possible_branch(node, depot, "req_arcs", arc, branch_pool, data)
                sum_arc = 0
                for (i,j) in enumerate(node["tours"])
                    if tour_pool[j]["depot"] == depot
                        sum_arc += node["alpha"][i]*tour_pool[j]["X"]["req_arcs"][arc]
                    end
                end
                push!(delta[depot]["req_arcs"], sum_arc)
            else
                push!(delta[depot]["req_arcs"], typemax(Int64))
            end
        end

        for arc in 1:data[1,:req_edges]
            if possible_branch(node, depot, "req_edges", arc, branch_pool, data)
                sum_arc = 0
                for (i,j) in enumerate(node["tours"])
                    if tour_pool[j]["depot"] == depot
                        sum_arc += node["alpha"][i]*tour_pool[j]["X"]["req_edges"][arc]
                    end
                end
                push!(delta[depot]["req_edges"], sum_arc)
            else
                push!(delta[depot]["req_edges"], typemax(Int64))
            end
        end

        for arc in 1:data[1,:req_edges]
            if possible_branch(node, depot, "req_edges_swap", arc, branch_pool, data)
                sum_arc = 0
                for (i,j) in enumerate(node["tours"])
                    if tour_pool[j]["depot"] == depot
                        sum_arc += node["alpha"][i]*tour_pool[j]["X"]["req_edges_swap"][arc]
                    end
                end
                push!(delta[depot]["req_edges_swap"], sum_arc)
            else
                push!(delta[depot]["req_edges_swap"], typemax(Int64))
            end
        end
    end

    min_depot = Dict()
    for depot in data[1, :depots]
        min_arcs = try
                findmin(abs.(delta[depot]["req_arcs"].-0.5))
            catch
                (1000000,1)
            end
        min_edges = try
                findmin(abs.(delta[depot]["req_edges"].-0.5))
            catch
                (1000000,1)
            end
        min_edges_swap = try
                findmin(abs.(delta[depot]["req_edges_swap"].-0.5))
            catch
                (1000000,1)
            end
        mins = [min_arcs, min_edges, min_edges_swap]
        type_min = ["req_arcs", "req_edges", "req_edges_swap"]
        min_tot = findmin([min_arcs[1], min_edges[1], min_edges_swap[1]])
        # Type d'arc/arête | Numéro de l'arc/arête | Valeur du min
        min_depot[depot] = [type_min[min_tot[2]], mins[min_tot[2]][2], mins[min_tot[2]][1]]
    end

    min_tot = findmin([val[3] for val in values(min_depot)])

    # Dépot | Type d'arc/arête | Numéro de l'arc/arête
    branch = [collect(keys(min_depot))[min_tot[2]], min_depot[collect(keys(min_depot))[min_tot[2]]][1], min_depot[collect(keys(min_depot))[min_tot[2]]][2]]

    branch_index = length(branch_pool)+1
    push!(branch_pool, branch)

    return(branch_index, branch_pool)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem