text:
using Random
using Plots
using DifferentialEquations,DiffEqFlux,LinearAlgebra
using Optim
include("../src/GLOQ.jl")
pyplot()
#import GLOQ
############################################################
# Step 1: Generate Synthetic Data
############################################################
H = [ 0 im 1.0;
     -im 0 0;
     1.0   0 0]
L1 = [0.0 0.25 0.0;
      0.0 0.0  0.5;
      0.0 0.0  0.0]
L2 = [0.0 0.0 0.0;
      0.0 0.5 0.0;
      0.0 0.0 1.0]
L_list = (L1,L2)
LH,LD = GLOQ.make_lindblad_operator_complex(H,L_list)
L = LH+LD
println("Lindblad operator assembled")
# define the time interval
tspan_length = 21
t_span = range(0.0,pi,length=tspan_length)
# forward solves with different initial conditions,
# use exponential_solver
u1 = [1.0;0;0]
rho1_history = GLOQ.exponential_solver_complex(u1,L,t_span;initial_type = "states")
u2 = [0;1.0;0]
rho2 = (u2*u2')[:]
rho2_history = GLOQ.exponential_solver_complex(u2,L,t_span;initial_type = "states")
u3 = [1.0/sqrt(2.0);-im/sqrt(2.0);0]
rho3_history = GLOQ.exponential_solver_complex(u3,L,t_span;initial_type = "states")
u4 = [0.0;1.0/sqrt(2.0);-im/sqrt(2.0)]
rho4_history = GLOQ.exponential_solver_complex(u4,L,t_span;initial_type = "states")
u5 = [1.0/sqrt(2.0);0.0;-im/sqrt(2.0)]
rho5_history = GLOQ.exponential_solver_complex(u5,L,t_span;initial_type = "states")

population1_history = GLOQ.get_population(rho1_history)
population2_history = GLOQ.get_population(rho2_history)
population3_history = GLOQ.get_population(rho3_history)
population4_history = GLOQ.get_population(rho4_history)
population5_history = GLOQ.get_population(rho5_history)

fig = plot(t_span,population1_history,lab=:false)
plot!(fig,t_span,population2_history,lab=:false)
plot!(fig,t_span,population3_history,lab=:false)
plot!(fig,t_span,population4_history,lab=:false)
plot!(fig,t_span,population5_history,lab=:false)
display(fig)

# add noise to the synthetic data
noise_amp = 0.025
noisy_data1 = population1_history+noise_amp.*rand(Float64,size(population1_history))
noisy_data2 = population2_history+noise_amp.*rand(Float64,size(population2_history))
noisy_data3 = population3_history+noise_amp.*rand(Float64,size(population3_history))
noisy_data4 = population4_history+noise_amp.*rand(Float64,size(population4_history))
noisy_data5 = population5_history+noise_amp.*rand(Float64,size(population5_history))
# present the noisy data
fig = plot(t_span,noisy_data1,lab=:false)
plot!(fig,t_span,noisy_data2,lab=:false)
plot!(fig,t_span,noisy_data3,lab=:false)
plot!(fig,t_span,noisy_data4,lab=:false)
plot!(fig,t_span,noisy_data5,lab=:false)
display(fig)
############################################################
# Step 2: Learn the operators
############################################################
rho1 = convert(Vector{ComplexF64},(u1*u1')[:])
rho2 = convert(Vector{ComplexF64},(u2*u2')[:])
rho3 = convert(Vector{ComplexF64},(u3*u3')[:])
rho4 = convert(Vector{ComplexF64},(u4*u4')[:])
rho5 = convert(Vector{ComplexF64},(u5*u5')[:])
# p[1]-p[5] upper diagonal of the Hamiltonian part
# p[6]-p[9] parameters in the Lindblad term
function loss(p)
      _X = [ 0 p[1]*im p[2];
             0   0  0;
             0   0  0]
      _H = 0.5.*(_X+adjoint(_X))
      _L1 = [0.0 p[3] 0.0;
            0.0 0.0  p[4];
            0.0 0.0  0.0]
      _L2 = [0.0 0.0 0.0;
            0.0 p[5] 0.0;
            0.0 0.0 2.0*p[5]]
      _L_list = (_L1,_L2)
      _LH,_LD = GLOQ.make_lindblad_operator_complex(_H,_L_list)
      _L = _LH+_LD

      _problem = GLOQ.LindbladODEProblemComplex(rho1,_L,t_span[end];initial_type="density")
      rho1_data = solve(_problem,saveat=t_span)
      rho2_data = solve(remake(_problem, u0=rho2), saveat=t_span)
      rho3_data = solve(remake(_problem, u0=rho3), saveat=t_span)
      rho4_data = solve(remake(_problem, u0=rho4), saveat=t_span)
      rho5_data = solve(remake(_problem, u0=rho5), saveat=t_span)

      pop_data1 = GLOQ.get_population(Array(rho1_data))
      pop_data2 = GLOQ.get_population(Array(rho2_data))
      pop_data3 = GLOQ.get_population(Array(rho3_data))
      pop_data4 = GLOQ.get_population(Array(rho4_data))
      pop_data5 = GLOQ.get_population(Array(rho5_data))

      _loss = sum(abs2, pop_data1-noisy_data1)
      _loss += sum(abs2, pop_data2-noisy_data2)
      _loss += sum(abs2, pop_data3-noisy_data3)
      _loss += sum(abs2, pop_data4-noisy_data4)
      _loss += sum(abs2, pop_data5-noisy_data5)
      return _loss
end

p_initial = [1.0;1.0;0.5;0.5;0.5]
#loss_initial = loss(p_initial)
#println(loss_initial)
p_true = [2.0 2.0 0.25 0.5 0.5]
@time p_DEFlux = DiffEqFlux.sciml_train(loss, p_initial)
@time p_Optim = Optim.optimize(loss, p_initial)
println("True p: ",p_true)
println("DEFlux p: ",p_DEFlux.u)
println("Optim p: ",p_Optim.minimizer)
println(p_DEFlux)
println(p_Optim)

################################
function loss3(p)
      _X = [ 0 p[1]*im p[2];
             0   p[3]  p[4];
             0   0     p[5]]
      _H = 0.5.*(_X+adjoint(_X))
      _L1 = [0.0 p[6] 0.0;
            0.0 0.0  p[7];
            0.0 0.0  0.0]
      _L2 = [0.0 0.0 0.0;
            0.0 p[8] 0.0;
            0.0 0.0 2.0*p[8]]
      _L_list = (_L1,_L2)
      _LH,_LD = GLOQ.make_lindblad_operator_complex(_H,_L_list)
      _L = _LH+_LD

      _problem = GLOQ.LindbladODEProblemComplex(rho1,_L,t_span[end];initial_type="density")
      rho1_data = solve(_problem,saveat=t_span)
      rho2_data = solve(remake(_problem, u0=rho2), saveat=t_span)
      rho3_data = solve(remake(_problem, u0=rho3), saveat=t_span)
      rho4_data = solve(remake(_problem, u0=rho4), saveat=t_span)
      rho5_data = solve(remake(_problem, u0=rho5), saveat=t_span)

      pop_data1 = GLOQ.get_population(Array(rho1_data))
      pop_data2 = GLOQ.get_population(Array(rho2_data))
      pop_data3 = GLOQ.get_population(Array(rho3_data))
      pop_data4 = GLOQ.get_population(Array(rho4_data))
      pop_data5 = GLOQ.get_population(Array(rho5_data))

      _loss = sum(abs2, pop_data1-noisy_data1)
      _loss += sum(abs2, pop_data2-noisy_data2)
      _loss += sum(abs2, pop_data3-noisy_data3)
      _loss += sum(abs2, pop_data4-noisy_data4)
      _loss += sum(abs2, pop_data5-noisy_data5)
      return _loss
end

p3_initial = [1.0;1.0;0.25;0.1;0.1;0.5;0.5;0.5]
@time p3_DEFlux = DiffEqFlux.sciml_train(loss3, p3_initial)
@time p3_Optim = Optim.optimize(loss3, p3_initial)
println("True p: ",p_true)
println("DEFlux p: ",p3_DEFlux.u)
println("Optim p: ",p3_Optim.minimizer)
println(p3_DEFlux)
println(p3_Optim)
################################################
#=
function loss2(p)
      _H = p[1]+1.0im.*p[2]
      _L1 = p[3]
      _L2 = p[4]
      _L_list = (_L1,_L2)
      _LH,_LD = GLOQ.make_lindblad_operator(_H,_L_list)
      _L = _LH+_LD

      _problem = GLOQ.LindbladODEProblem(rho1,_L,t_span[end];initial_type="density")
      rho1_data = solve(_problem,saveat=t_span)
      rho2_data = solve(remake(_problem, u0=rho2), saveat=t_span)
      rho3_data = solve(remake(_problem, u0=rho3), saveat=t_span)
      rho4_data = solve(remake(_problem, u0=rho4), saveat=t_span)
      rho5_data = solve(remake(_problem, u0=rho5), saveat=t_span)

      pop_data1 = GLOQ.get_population(Array(rho1_data))
      pop_data2 = GLOQ.get_population(Array(rho2_data))
      pop_data3 = GLOQ.get_population(Array(rho3_data))
      pop_data4 = GLOQ.get_population(Array(rho4_data))
      pop_data5 = GLOQ.get_population(Array(rho5_data))

      _loss = sum(abs2, pop_data1-noisy_data1)
      _loss += sum(abs2, pop_data2-noisy_data2)
      _loss += sum(abs2, pop_data3-noisy_data3)
      _loss += sum(abs2, pop_data4-noisy_data4)
      _loss += sum(abs2, pop_data5-noisy_data5)
      return _loss
end

Hr_initial = [ 0 0 1.2;
               0 0 0;
            1.2  0 0]
Hi_initial = [ 0 1.5 0;
              -1.5 0 0;
               1  0 0]
L1_initial = [0.0 0.5 0.0;
              0.0 0.0 0.3;
              0.0 0.0 0.0]
L2_initial = [0.0 0.0 0.0;
              0.0 0.3 0.0;
              0.0 0.0 0.6]
p2_initial = (Hr_initial,Hi_initial,L1_initial,L2_initial)
loss2_initial = loss2(p2_initial)

@time p_opt2 = DiffEqFlux.sciml_train(loss, p2_initial,BFGS(),maxiters=100)
#println("True p: ",p_true)
println("Optimized Operators: ",p_opt2.u)
println(H," ",L1," ",L2)
=#


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule tippecanoe_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("tippecanoe")
JLLWrappers.@generate_main_file("tippecanoe", UUID("6e4bcf07-11bd-542f-a056-5857bcc30079"))
end  # module tippecanoe_jll

ext: jl


text:
function incident_infections(params::InfParams, S::Int, I::Vector{Int}, R::Int)
    iszero(S) && return 0
    infSum = 0
    N = S + sum(I) + R
    for (i, inf) in enumerate(I)
        d = params.infectiousness[i]
        k,θ = Distributions.params(d)
        d′ = Gamma(k*S/N, θ) # conversion from R₀ to Rₜ
        infSum += rand(RVsum(d′, inf))
    end

    return min(S, floor(Int, infSum)) # Can't infect more people than are susceptible
end

function symptomatic_isolation(params::InfParams, I::Vector{Int})
    isolating = zero(I)
    for (i, inf) in enumerate(I)
        symptomatic_prob = params.symptom_probs[i]
        isolation_prob = symptomatic_prob*(1-params.asymptomatic_prob)*params.symptomatic_isolation_prob
        isolating[i] = rand(Binomial(inf,isolation_prob))
    end
    return isolating
end

function positive_tests(params::InfParams, I::Vector{Int}, T::Matrix{Int}, a::CovidAction)
    pos_tests = zeros(Int, length(params.pos_test_probs))

    for (i, inf) in enumerate(I)
        num_already_tested = sum(@view T[:,i])
        num_tested = floor(Int,(inf-num_already_tested)*a.testing_prop)
        pos_tests[i] = rand(Binomial(num_tested,params.pos_test_probs[i]))
    end
    return pos_tests
end

# Only record number that have taken the test, the number that return postive is
# Binomial dist, such that s' is stochastic on s.
function update_isolations(params::InfParams, I, R, T, a::CovidAction)

    sympt = symptomatic_isolation(params, I) # Number of people isolating due to symptoms
    pos_tests = positive_tests(params, I, T, a)

    sympt_prop = sympt ./ I # Symptomatic Isolation Proportion
    replace!(sympt_prop, NaN=>0.0)

    R += sum(sympt)
    I .-= sympt

    T[end,:] .= pos_tests

    @. T = floor(Int, (1 - sympt_prop)' * T)

    R += sum(@view T[1,:])
    @views I .-= T[1,:]

    @assert all(≥(0), I)

    CovidPOMDPs.shift_test!(T)

    return I, R, T, pos_tests
end

function CovidPOMDPs.sim_step(pomdp::SingleCovidPOMDP, state::SingleCovidState, a::CovidAction)
    (;S, I, R, T, params, prev_action) = state
    I = copy(I)
    T = copy(T)

    # Update symptomatic and testing-based isolations
    I, R, T, pos_tests = update_isolations(params, I, R, T, a)

    # Incident Infections
    R += I[end]
    I = circshift(I, 1) # shift_inf!(I)
    new_infections = incident_infections(params, S, I, R)
    I[1] = new_infections
    S -= new_infections
    sp = SingleCovidState(S, I, R, T, params, a)

    return sp, new_infections, pos_tests
end


text:
function __init__()
    init_acache()
    init_bcache()
    nt = init_num_tasks()
    if nt < num_cores() && ("SUPPRESS_OCTAVIAN_WARNING" ∉ keys(ENV))
        msg = string(
            "Your system has $(num_cores()) physical cores, but `Octavian.jl` only has ",
            "$(nt > 1 ? "$(nt) threads" : "$(nt) thread") available. ",
            "For the best performance, you should start Julia with at least $(num_cores()) threads.",
        )
        @warn msg
    end
    reseet_bcache_lock!()
end

function init_bcache()
    BCACHEPTR[] = VectorizationBase.valloc(second_cache_size() * bcache_count(), Cvoid, ccall(:jl_getpagesize, Int, ()))
    nothing
end

function init_acache()
    Sys.WORD_SIZE ≤ 32 || return
    ACACHEPTR[] = VectorizationBase.valloc(first_cache_size() * init_num_tasks(), Cvoid, ccall(:jl_getpagesize, Int, ()))
    nothing
end

function init_num_tasks()
    num_tasks = _read_environment_num_tasks()::Int
    OCTAVIAN_NUM_TASKS[] = num_tasks
end

function _read_environment_num_tasks()
    environment_variable = get(ENV, "OCTAVIAN_NUM_TASKS", "")::String
    nt = min(Threads.nthreads(), VectorizationBase.num_cores())::Int
    if isempty(environment_variable)
        return nt
    else
        return min(parse(Int, environment_variable)::Int, nt)
    end
end


text:
@enum VersionSpec begin
    patch
    minor
    major
    current
end

"""
    release(version_spec[, path=pwd()]; kw...)

Release a new version for a package, e.g

```julia-repl
julia> pwd()
"/home/roger/julia/Comonicon"

julia> Ion.release(patch) # will release a new patch version on github
```

# Arguments

- `version_spec`: required, semantic version bump (`patch`, `minor`, `major`, `current`)
    or a specific version number to release (`VersionNumber` or valid version string).
- `path`: optional, path to the pacakge to release, default is `pwd()`.

# Keyword Arguments

- `registry`: optional, name of the registry to register, if not specified, will lookup
    if the package has been registered in local registries.
- `branch`: optional, branch to register, default is current branch.
- `note`: optional, release note.
"""
function release(
        version_spec::Union{VersionSpec, VersionNumber, String},
        path::String=pwd();
        registry::String="",
        branch::String=Internal.current_branch(path),
        note::String="", debug::Bool=false,
    )

    spec = if version_spec == patch
        "patch"
    elseif version_spec == minor
        "minor"
    elseif version_spec == major
        "major"
    else
        "current"
    end

    try
        Internal.release(spec, path; registry, branch, note, debug)
    catch e
        if e isa Comonicon.CommandExit
            return e.exitcode
        else
            rethrow(e)
        end
    end
end

for api in [:clone, :compat]
    md = @eval @doc(Internal.$api)
    @eval @doc $md function $api(args...; kw...)
        try
            Internal.$api(args...; kw...)
        catch e
            if e isa Comonicon.CommandExit
                return e.exitcode
            else
                rethrow(e)
            end
        end
    end
end


text:
"""
We define a `MultiRandomWalk` struct herein. This is a discrete-time, discrete-space
stochastic process. A `MultiRandomWalk` is a `RandomWalk` with multiple agents walking on
the same domain.
"""

"""
    MultiRandomWalk{S,T,U<:Sampleable} <: AbstractMarkovChain

##### Fields
- `domain::S` : the domain within which walkers move (e.g. an integer lattice or a graph).
- `initialstate::T<:AbstractVector` : the state at which the random walk begins (at
time-step 1) -- contains each walker's state.
- `stepsampler::U<:MultiSampler` : sampler that allows i.i.d. steps to be drawn.
"""
struct MultiRandomWalk{S,T<:AbstractVector,U<:MultiSampler} <: AbstractMarkovChain
    domain::S
    initialstate::T
    stepsampler::U
end

"""
    Base.show(io::IO, mrw::MultiRandomWalk) -> nothing

Show a `MultiRandomWalk` in a human-friendly manner.
"""
function Base.show(io::IO, mrw::MultiRandomWalk)
    println(io, typeof(mrw))
    println(io, " Domain:")
    println(io, "  ", mrw.domain)
    println(io, " Number of walkers:")
    println(io, "  ", length(mrw.initialstate))
    println(io, " Initial state:")
    println(io, "  ", mrw.initialstate)
    println(io, "")
    println(io, " MultiSampler that allows i.i.d. steps to be drawn:")
    println(io, "  ", mrw.stepsampler)
    return nothing
end

"""
    Base.eltype(mrw::MultiRandomWalk) -> DataType

Return the type associated with `mrw`'s domain.
"""
function Base.eltype(mrw::MultiRandomWalk)
    return typeof(mrw.initialstate[1])
end

"""
    MultiRandomWalk(nwalkers::Integer, latticedim::Integer = 1,
                    initialstate::AbstractVector{<:AbstractVector{<:Integer}} = (
                                                  repeat(zeros(Int, latticedim), nwalkers)),
                    stepsampler::MultiSampler = MultiSampler(
                            repeat([EVecSampler(latticedim)], nwalkers))) -> MultiRandomWalk

Construct and return a `MultiRandomWalk` on the integer lattice with dimension `latticedim`.

##### Complexity
- Time complexity: `O(1)`.
"""
function MultiRandomWalk(nwalkers::Integer, latticedim::Integer = 1,
                         initialstate::AbstractVector = (
                                              [zeros(Int, latticedim) for _ in 1:nwalkers]),
                         stepsampler::MultiSampler = MultiSampler(
                                               repeat([EVecSampler(latticedim)], nwalkers)))
    # check arguments
    latticedim < 1 && throw(ArgumentError("latticedim must be positive."))
    nwalkers < 1 && throw(ArgumentError("nwalkers must be positive."))
    for i in 1:length(initialstate)
        length(initialstate[i]) != latticedim &&
            throw(DimensionMismatch("initialstate's elements' lengths should equal
                                    latticedim"))
    end

    if latticedim == 1
        return MultiRandomWalk(Int, initialstate, stepsampler)
    else
        return MultiRandomWalk(repeat([Int], latticedim), initialstate, stepsampler)
    end
end

"""
    MultiRandomWalk(graph::AbstractGraph, nwalkers::Integer,
                    initialvertex::AbstractVector = repeat([rand(1:nv(graph))], nwalkers);
                    checkargs::Bool = true) -> MultiRandomWalk

Construct and return a `MultiRandomWalk` on a graph.

##### Complexity
- Time complexity: `O(n)` when `checkargs == true`; otherwise, `O(1)`.
"""
function MultiRandomWalk(graph::AbstractGraph, nwalkers::Integer,
                        initialvertex::AbstractVector=repeat([rand(1:nv(graph))], nwalkers);
                        checkargs::Bool = true)

    if checkargs && !(initialvertex in vertices(graph))
        throw(ArgumentError("initialvertex must be a vertex in graph)."))
    else
        return RandomWalk(graph, initialvertex, NeighborSampler(graph, initialvertex))
    end
end

"""
    nwalkers(mrw::MultiRandomWalk) -> Int

Return the number of walkers in `mrw`.

##### Complexity
- Time complexity: `O(1)`.
"""
function nwalkers(mrw::MultiRandomWalk)
    return length(initialstate)
end

max_stars_repo_name: LoganDGraham/StochasticProcesses

max_issues_repo_name: LoganDGraham/StochasticProcesses

max_forks_repo_name: LoganDGraham/StochasticProcesses


text:
function rotation_matrix(roll::Real, pitch::Real, yaw::Real)
    R_roll =  [1.0      0.0          0.0;
               0.0      cos(roll)   -sin(roll);
               0.0      sin(roll)   cos(roll)]

    R_pitch = [cos(pitch)   0.0      sin(pitch);
               0.0          1.0      0.0;
               -sin(pitch)  0.0      cos(pitch)]

    R_yaw = [cos(yaw)     -sin(yaw)  0.0;
             sin(yaw)     cos(yaw)   0.0;
             0.0          0.0        1.0]

    return R_roll * R_pitch * R_yaw
end

# TODO accept camera intrinsics
function simple_projection(znear, zfar)
    [
        1             0              0                                  0;
        0             1              0                                  0;
        0             0              -(zfar + znear)/(zfar - znear)     -2*zfar*znear/(zfar - znear);
        0             0              -1                                 0
    ]
end

export rotation_matrix

max_stars_repo_name: probcomp/DepthRenderer

max_issues_repo_name: probcomp/DepthRenderer

max_forks_repo_name: probcomp/DepthRenderer


text:
## Optimization Model for dtu42380a1
## Strategy I: New production lines
## Mar 11, 2020

using JuMP, Gurobi, CSV

include("./get_data.jl")

model = Model(with_optimizer(Gurobi.Optimizer, Presolve=0, OutputFlag=0));

## PRODUCTION VARIABLES
@variable(model, x[1:L,1:P,1:T] >= 0)
## SOURCE VARIABLES (PREBLENDING AND PACKAGING)
@variable(model, y[1:L,1:S,1:P,1:T] >= 0);
## MACHINE FLOW VARIABLES (ALLOCATION OF REGULAR TIME AND OVERTIME FOR BLENDING AND PACKAGING MACHINES)
@variable(model, mBR[1:L, 1:P, 1:T] >= 0);
@variable(model, mBO[1:L, 1:P, 1:T] >= 0);
@variable(model, mPR[1:L, 1:P, 1:T] >= 0);
@variable(model, mPO[1:L, 1:P, 1:T] >= 0);
## INVENTORY VARIABLES
@variable(model, I[1:L,1:P,1:T] >= 0);
## DEMAND NOT MET VARIABLES
@variable(model, delta[1:L,1:P,1:T] >= 0);

@variable(model, s_b[1:L] >= 1)
@variable(model, s_p[1:L] >= 1)

## New decision variables regarding promotion or not
@variable(model, whe_pro[1:L, 1:2], Bin);

## OBJECTIVE FUNCTION
@objective(model, Min, sum( #=
=#           sum(y[i,s,p,t]*c[s] for s = 1:S)    #=  material cost
=#         + mBR[i,p,t]*cBR                      #= regular time blending machine cost
=#         + mBO[i,p,t]*cBO                      #= overtime blending machine cost
=#         + mPR[i,p,t]*cPR                      #= regular time packaging machine cost
=#         + mPO[i,p,t]*cPO                      #= overtime time packaging machine cost
=#         + I[i,p,t]*h[p]                       #= inventory cost
=#         + PI*delta[i,p,t]                     #= unsatisfied demand penalty
=#         for i = 1:L, p = 1:P, t = 1:T)
           + sum(s_b[i] + s_p[i] for i = 1:L) * 1000000
);

## PROCESS CONSTRAINTS - material and machine inflow must match (according to
##   the requirements) product outflow
@constraint(model, [i = 1:L, p = 1:P, t = 1:T,s = 1:S],
  y[i,s,p,t] == pi[p,s] * x[i,p,t]);

@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  mBR[i,p,t] + mBO[i,p,t] >=  muB[p] * x[i,p,t] );
@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  mPR[i,p,t] + mPO[i,p,t] >= muP[p] * x[i,p,t]);

## STORAGE CONSTRAINTS - inflow (production and start-of-period inventory) must
##   equal outflow (demand and end-of-period inventory) at each DC

## Constraints regarding promotion
@constraint(model, [i = 1:L, p = 3, t = 1],
  x[i,p,t] + I0[i,p] == I[i,p,t] + D[i,p,t] +
  0.3 * whe_pro[i, 1] * D[i,p,t] - delta[i,p,t]);
@constraint(model, [i = 1:L, p = 3, t = 2],
  x[i,p,t] + I[i,p,t-1] == I[i,p,t] + D[i,p,t] +
  0.3 * whe_pro[i, 2] * D[i,p,t] - delta[i,p,t]);
@constraint(model, [i = 1:L], whe_pro[i, 1] + whe_pro[i, 2] >= 1)

@constraint(model, [i = 1:L, p = 1:2, t = 1],
  x[i,p,t] + I0[i,p] == I[i,p,t] + D[i,p,t] - delta[i,p,t]);
@constraint(model, [i = 1:L, p = 1:2, t = 2],
  x[i,p,t] + I[i,p,t-1] == I[i,p,t] + D[i,p,t] - delta[i,p,t]);
@constraint(model, [i = 1:L, p = 1:P, t = 3:4],
  x[i,p,t] + I[i,p,t-1] == I[i,p,t] + D[i,p,t] - delta[i,p,t]);

@constraint(model, [i = 1:L, p = 1:P, t = 1:T],
  I[i,p,t] >= ss[i,p,t]);

## Modifications
## CAPACITY CONSTRAINTS
@constraint(model, [i = 1:L, t = 1:T],
  sum(mBR[i,p,t] for p = 1:P) <= 24*5*60*60 * s_b[i]);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mPR[i,p,t] for p = 1:P) <= 24*5*60*60 * s_p[i]);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mBO[i,p,t] for p = 1:P) <= 24*1*60*60 * s_b[i]);
@constraint(model, [i = 1:L, t = 1:T],
  sum(mPO[i,p,t] for p = 1:P) <= 24*1*60*60 * s_p[i]);
@constraint(model, [i = 1:L, p = 1:P, t = 1:T], delta[i,p,t] <= 0.00001);

optimize!(model)
println(value(whe_pro[1, 1]))
println(value(whe_pro[1, 2]))
println(value(whe_pro[2, 1]))
println(value(whe_pro[2, 2]))
println(value(whe_pro[3, 1]))
println(value(whe_pro[3, 2]))

for i = 1:L
  println(value(s_b[i]))
  println(value(s_p[i]))
end
extra = value(s_b[1]) + value(s_p[1]) + value(s_b[2]) +
  value(s_p[2]) + value(s_b[3]) + value(s_p[3])

println(objective_value(model) - extra * 1000000)

include("./export_result.jl")


text:
"""
    Configuration(configfile::String)

Parses configuration file `configfile`.
"""
struct Configuration
    equation_name::String
    scenario_name::String
    equation_type::String
    source::String
    end_time::Float64
    order::Int64
    timeintegrator_name::String
    riemann::String
    courant::Float64

    filter_name::String
    filter_order::Int64
    plot_start::Float64
    plot_step::Float64
    grid_elements::Int64
    physicalsize::Array{Float64,1}
    cellsize::Array{Float64,1}

    function Configuration(configfile::String)
        config = YAML.load(open(configfile))
        equation_name = config["equation"]["equation"]
        scenario_name = config["equation"]["scenario"]
        equation_type = config["equation"]["type"]
        source = config["equation"]["source"]
        end_time = config["simulation"]["end_time"]
        order = config["solver"]["order"]
        timeintegrator_name = config["solver"]["timeintegrator"]
        riemann = config["solver"]["riemann"]
        courant = config["solver"]["courant"]

        filter_name = get(config["solver"], "filter", "identity")
        filter_order = get(config["solver"], "filter_order", order)
        plot_start = config["output"]["start"]
        plot_step = config["output"]["step"]
        grid_elements = config["simulation"]["grid_elements"]
        physicalsize_1d = config["simulation"]["grid_size"]
        physicalsize = [physicalsize_1d, physicalsize_1d]
        cellsize = physicalsize ./ grid_elements

        new(
            equation_name,
            scenario_name,
            equation_type,
            source,
            end_time,
            order,
            timeintegrator_name,
            riemann,
            courant,
            filter_name,
            filter_order,
            plot_start,
            plot_step,
            grid_elements,
            physicalsize,
            cellsize
        )
    end
end


text:
using GingerGUI
using Test

@testset "GingerGUI.jl" begin
    # Write your own tests here.
end

max_stars_repo_name: pacific-zitter/GingerGUI

max_issues_repo_name: pacific-zitter/GingerGUI

max_forks_repo_name: pacific-zitter/GingerGUI


text:
# run 2 processor tests

push!(LOAD_PATH, abspath(joinpath(pwd(), "..")))

using PDESolver
#using Base.Test
using Base.Test
using ODLCommonTools
using PdePumiInterface  # common mesh interface - pumi
using SummationByParts  # SBP operators
using AdvectionEquationMod
using LinearSolvers
using NonlinearSolvers   # non-linear solvers
using OptimizationInterface
using ArrayViews
import ArrayViews.view
using Utils
using Input

#------------------------------------------------------------------------------
# define tests and tags

#include("../TestSystem.jl")
using TestSystem
# define tags that will be used
include("../tags.jl")

# test list
global const AdvectionTests = TestList()

"""
  Run parallel tests and compare to serial results calculated as part of
  serial tests.
"""
function runtests_parallel()
  @testset "----- Testing Parallel -----" begin

    start_dir = pwd()
    cd("./rk4/parallel")
    fname = "input_vals_parallel_runp.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")

    @test isapprox( datas[1], datap[1]) atol=1e-13
    @test isapprox( datas[2], datap[2]) atol=1e-13
    cd("../../")

    cd("./newton/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("./error_calc.dat")
    @test isapprox( datas[1], datap[1]) atol=1e-13

    cd(start_dir)

    cd("./rk4_3d/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")
    @test isapprox( datas[1], datap[1]) atol=1e-13

    cd(start_dir)

    cd("./newton_3d/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    datas = readdlm("../serial/error_calc.dat")
    datap = readdlm("error_calc.dat")
    @test isapprox( datas[1], datap[1]) atol=1e-13

    cd(start_dir)
  end  # end facts block

  return nothing
end

#runtests_parallel()
add_func1!(AdvectionTests, runtests_parallel, [TAG_SHORTTEST])

function test_precompute()
  @testset "----- testing non-precompute functions -----" begin
    start_dir = pwd()

    # test rk4
    cd("./rk4/parallel")
    fname = "input_vals_parallel_runp.jl"
    #TODO: set opts["solve"] = false before doing this
    mesh, sbp, eqn, opts = solvePDE(fname)

    fill!(eqn.res, 0.0)
    evalResidual(mesh, sbp, eqn, opts)

    res_orig = copy(eqn.res)

    opts["precompute_face_flux"] = false
    evalResidual(mesh, sbp, eqn, opts)

    @test isapprox( norm(vec(eqn.res - res_orig)), 0.0) atol=1e-13

    # test newton
    cd(start_dir)
    cd("./newton/parallel")
    fname = "input_vals_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    fill!(eqn.res, 0.0)
    evalResidual(mesh, sbp, eqn, opts)
    res_orig = copy(eqn.res)

    opts["precompute_face_flux"] = false
    evalResidual(mesh, sbp, eqn, opts)

    @test isapprox( norm(vec(eqn.res - res_orig)), 0.0) atol=1e-13

    cd(start_dir)
  end

  return nothing
end

add_func1!(AdvectionTests, test_precompute, [TAG_SHORTTEST])

function test_adjoint_parallel()

  @testset "--- Testing Adjoint Computation on a Geometric Boundary ---" begin

    fname = "input_vals_functional_DG_parallel.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    objective = createFunctional(mesh, sbp, eqn, opts, 1)
    evalFunctional(mesh, sbp, eqn, opts, objective)
    pc, lo = getNewtonPCandLO(mesh, sbp, eqn, opts)
    ls = StandardLinearSolver(pc, lo, eqn.comm, opts)

    adjoint_vec = zeros(Complex{Float64}, mesh.numDof)
    calcAdjoint(mesh, sbp, eqn, opts, ls, objective, adjoint_vec, recalc_jac=true, recalc_pc=true)

    for i = 1:length(adjoint_vec)
      @test isapprox( real(adjoint_vec[i]), 1.0) atol=1e-10
    end

  end # End testset("--- Testing Functional Computation on a Geometric Boundary ---")

  return nothing
end

add_func1!(AdvectionTests, test_adjoint_parallel, [TAG_ADJOINT, TAG_LONGTEST])

#------------------------------------------------------------------------------
# run tests
@testset "----- Running Advection 2 processor tests -----" begin
  runTestSystem(AdvectionTests, solvePDE, ARGS)
end

#------------------------------------------------------------------------------
# cleanup


text:
module Transports

include("./transports/common.jl")
include("./transports/memory.jl")
include("./transports/tcp.jl")

end


text:
using CSV, DataFrames, Test, Plots, Binscatters
df = DataFrame(CSV.File(joinpath(dirname(pathof(Binscatters)), "../dataset/Cigar.csv")))

# bin
out = Binscatters.bin(df, @formula(Sales ~ Price))
@test issorted(out.Price)
@test minimum(out.Price) >= minimum(df.Price)
@test maximum(out.Price) <= maximum(df.Price)
Binscatters.bin(df, @formula(Sales ~ Price), weights = :Pop)
@test minimum(out.Price) >= minimum(df.Price)
@test maximum(out.Price) <= maximum(df.Price)
Binscatters.bin(df, @formula(Sales ~ Price), 10)
Binscatters.bin(df, @formula(Sales ~ Price + NDI))
Binscatters.bin(df, @formula(Sales ~ Price + fe(State)))
Binscatters.bin(df, @formula(Sales + NDI ~ Price))
df.Price_missing = ifelse.(df.Price .>= 30, df.Price, missing)
Binscatters.bin(df, @formula(Sales ~ Price_missing))
Binscatters.bin(groupby(df, :State), @formula(Sales ~ Price))

binscatter(df, @formula(Sales ~ Price))
binscatter(df, @formula(Sales ~ Price), seriestype = :scatterpath)
binscatter(df, @formula(Sales ~ Price), seriestype = :linearfit)

binscatter(df, @formula(Sales + NDI ~ Price))
binscatter(df, @formula(Sales + NDI ~ Price), seriestype = :scatterpath)
binscatter(df, @formula(Sales + NDI ~ Price), seriestype = :linearfit)

binscatter(df, @formula(Sales ~ Price + fe(State)))
binscatter(df, @formula(Sales ~ Price), weights = :Pop)
binscatter(df, @formula(Sales + NDI~ Price))
binscatter(df, @formula(Sales ~ Price + NDI))

df.dummy = df.State .>= 25
binscatter(groupby(df, :dummy), @formula(Sales ~ Price))
binscatter(groupby(df, :dummy), @formula(Sales ~ Price), seriestype = :scatterpath)
binscatter(groupby(df, :dummy), @formula(Sales ~ Price), seriestype = :linearfit)

binscatter(groupby(df, :dummy), @formula(Sales + NDI ~ Price))
binscatter(groupby(df, :dummy), @formula(Sales + NDI ~ Price), seriestype = :scatterpath)
binscatter(groupby(df, :dummy), @formula(Sales + NDI ~ Price), seriestype = :linearfit)


text:
"""
    add_community_objective!(
        community::CoreModel,
        objective_mets_weights::Dict{String, Float64};
        objective_id = "community_biomass",
    )

Add an objective column to the `community` model with optional id `objective_id`. Supply a
dictionary mapping the string names of the objective metabolites to their weights in
`objective_mets_weights`. Note, the weights are negated inside the function so that positive
weights are seen as reagents/substrates, NOT products in the reaction equation.

# Example
```
add_community_objective!(model, Dict("met1"=>1.0, "met2"=>2.0))
```

See also: [`update_community_objective!`](@ref)
"""
function add_community_objective!(
    community::CoreModel,
    objective_mets_weights::Dict{String,Float64};
    objective_id = "community_biomass",
)
    obj_inds = indexin(keys(objective_mets_weights), metabolites(community))
    nothing in obj_inds && throw(ArgumentError("Some objective metabolite(s) not found."))

    nr, _ = size(community.S)
    objcol = spzeros(nr)
    objcol[obj_inds] .= -collect(values(objective_mets_weights))

    # extend model by one reaction
    community.S = hcat(community.S, objcol)
    community.xl = [community.xl; 0.0]
    community.xu = [community.xu; _constants.default_reaction_bound]
    community.rxns = [community.rxns; objective_id]
    community.c = spzeros(size(community.S, 2))
    community.c[end] = 1.0

    return nothing # stop annoying return value
end

"""
    add_community_objective!(
        community::StandardModel,
        objective_mets_weights::Dict{String, Float64};
        objective_id = "community_biomass"
    )

Variant of [`add_community_objective!`] that takes a `StandardModel` community model as input.
"""
function add_community_objective!(
    community::StandardModel,
    objective_mets_weights::Dict{String,Float64};
    objective_id = "community_biomass",
)
    nothing in indexin(keys(objective_mets_weights), metabolites(community)) &&
        throw(ArgumentError("Some objective metabolite(s) not found."))

    rdict = Dict(k => -float(v) for (k, v) in objective_mets_weights)
    rxn = Reaction(objective_id)
    rxn.metabolites = rdict
    rxn.lb = 0.0
    rxn.ub = _constants.default_reaction_bound
    rxn.objective_coefficient = 1.0
    community.reactions[rxn.id] = rxn

    return nothing # stop annoying return value
end

"""
    update_community_objective!(
        community::CoreModel,
        objective_id::String,
        objective_mets_weights::Dict{String, Float64}
    )

Update the weights for the objective column with id `objective_id` in `community` using
`objective_mets_weights`, which maps metabolite ids to weights. The current weights are
reset to 0 before being updated to the supplied weights. Note, the weights are negated
inside the function so that the objective metabolites are seen as reagents/substrates, NOT
products in the reaction equation.

# Example
```
update_community_objective!(model, "community_biomass", Dict("met1"=>1.0, "met2"=>2.0))
```

See also: [`add_community_objective!`](@ref)
"""
function update_community_objective!(
    community::CoreModel,
    objective_id::String,
    objective_mets_weights::Dict{String,Float64},
)
    obj_inds = indexin(keys(objective_mets_weights), metabolites(community))
    nothing in obj_inds && throw(ArgumentError("Some objective metabolite(s) not found."))

    objective_column_index = first(indexin([objective_id], reactions(community)))
    community.S[:, objective_column_index] .= 0.0 # reset
    community.S[obj_inds, objective_column_index] .=
        -collect(values(objective_mets_weights))
    dropzeros!(community.S)
    community.c = spzeros(size(community.S, 2))
    community.c[objective_column_index] = 1.0
    community.xl[objective_column_index] = 0.0
    community.xu[objective_column_index] = _constants.default_reaction_bound

    return nothing # stop annoying return value
end

"""
    update_community_objective!(
        community::StandardModel,
        objective_id::String,
        objective_mets_weights::Dict{String, Float64}
    )

Variant of [`update_community_objective!`] that takes a `StandardModel` community model as input.
"""
function update_community_objective!(
    community::StandardModel,
    objective_id::String,
    objective_mets_weights::Dict{String,Float64},
)
    delete!(community.reactions, objective_id)
    add_community_objective!(community, objective_mets_weights; objective_id = objective_id)

    return nothing # stop annoying return value
end

"""
    join_with_exchanges(
        ::Type{CoreModel},
        models::Vector{M},
        exchange_rxn_mets::Dict{String,String};
        biomass_ids = String[],
        model_names = String[],
    ) where {M<:MetabolicModel}

Return a `CoreModel` representing the community model of `models` joined through their
exchange reactions and metabolites in the dictionary `exchange_rxn_mets`, which maps
exchange reactions to their associated metabolite. These exchange reactions and metabolites
link model metabolites to environmental metabolites and reactions. Optionally specify
`model_names` to append a specific name to each reaction and metabolite of an organism for
easier reference (default is `species_i` for each model index i in `models`). Note, the
bounds of the environmental variables are all set to zero. Thus, to run a simulation you
need to constrain them appropriately. All the other bounds are inherited from the models
used to construct the community model.

If `biomass_ids` is supplied, then a community model is returned that has an extra reaction
added to the end of the stoichiometric matrix (last column) that can be assigned as the
objective reaction. It also creates biomass "metabolites" that can be used in this objective
reaction. In the returned mode, these biomass metabolites are produced by the reaction
corresponding to `biomass_ids` in each model respectively. Note, this reaction is
unspecified, further action needs to be taken to specify it, e.g. assign weights to the last
column of the stoichiometric matrix in the rows corresponding to the biomass metabolites.

To further clarify how this `join` works. Suppose you have 2 organisms with stoichiometric
matrices S₁ and S₂ and you want to link them with `exchange_rxn_mets = Dict(er₁ => em₁, er₂
=> em₂, er₃ => em₃, ...)`. Then a new community stoichiometric matrix is constructed that
looks like this:
```
            _      er₁  er₂  er₃  ...  b_
           |S₁                           |
           |   S₂                        |
        em₁|                             |
S   =   em₂|                             |
        em₃|                             |
        ...|                             |
        bm₁|                             |
        bm₂|_                           _|

```
The exchange reactions in each model get linked to environmental metabolites, `emᵢ`, and
these get linked to environmental exchanges, `erᵢ`. These `erᵢ` behave like normal single
organism exchange reactions. When `biomass_ids` are supplied, each model's biomass reaction
produces a pseudo-metabolite (`bmᵢ`). These can be weighted in column `b`, called the
`community_biomass` reaction in the community model, if desired. Refer to the tutorial if
this is unclear.

# Example
```
m1 = load_model(core_model_path)
m2 = load_model(CoreModel, core_model_path)

# need to list ALL the exchanges that will form part of the entire model
exchange_rxn_mets = Dict(k => first(keys(reaction_stoichiometry(m1, ex_rxn)))
    for filter(looks_like_exchange_reaction, reactions(m1)))

biomass_ids = ["BIOMASS_Ecoli_core_w_GAM", "BIOMASS_Ecoli_core_w_GAM"]

community = join_with_exchanges(
    CoreModel,
    [m1, m2],
    exchange_rxn_mets;
    biomass_ids = biomass_ids,
)
```
"""
function join_with_exchanges(
    ::Type{CoreModel},
    models::Vector{M},
    exchange_rxn_mets::Dict{String,String};
    biomass_ids = String[],
    model_names = String[],
) where {M<:MetabolicModel}

    exchange_rxn_ids = keys(exchange_rxn_mets)
    exchange_met_ids = values(exchange_rxn_mets)
    add_biomass_objective = isempty(biomass_ids) ? false : true

    # get offsets to construct community S
    reaction_lengths = [n_reactions(model) for model in models]
    metabolite_lengths = [n_metabolites(model) for model in models]
    reaction_offset = [0; cumsum(reaction_lengths[1:end-1])]
    metabolite_offset = [0; cumsum(metabolite_lengths[1:end-1])]

    # get each model's S matrix (needed for the size calculations)
    stoichs = [stoichiometry(model) for model in models]
    nnzs = [findnz(stoich) for stoich in stoichs] # nonzero indices per model

    # size calculations
    column_add = add_biomass_objective ? 1 : 0 # objective rxn
    row_add = add_biomass_objective ? length(models) : 0 # biomass as metabolites
    nnz_add = add_biomass_objective ? (1 + length(models)) : 0
    nnz_total =
        sum(length(first(nnz)) for nnz in nnzs) +
        length(models) * length(exchange_rxn_ids) +
        length(exchange_met_ids) +
        nnz_add
    n_reactions_metabolic = sum(reaction_lengths)
    n_reactions_total = n_reactions_metabolic + length(exchange_rxn_ids) + column_add
    n_metabolites_metabolic = sum(metabolite_lengths)
    n_metabolites_total = n_metabolites_metabolic + length(exchange_met_ids) + row_add

    # Create empty storage vectors
    lbs = spzeros(n_reactions_total)
    ubs = spzeros(n_reactions_total)
    rxns = Array{String,1}(undef, n_reactions_total)
    mets = Array{String,1}(undef, n_metabolites_total)
    I = Array{Int,1}(undef, nnz_total)
    J = Array{Int,1}(undef, nnz_total)
    V = Array{Float64,1}(undef, nnz_total)

    # build metabolic components, block diagonals
    kstart = 1
    for i = 1:length(models)
        kend = kstart + length(nnzs[i][3]) - 1
        rng = kstart:kend
        I[rng] .= nnzs[i][1] .+ metabolite_offset[i]
        J[rng] .= nnzs[i][2] .+ reaction_offset[i]
        V[rng] .= nnzs[i][3]
        kstart += length(nnzs[i][3])
    end

    # build environmental - exchange links
    for i = 1:length(models)
        exchange_rxn_inds = indexin(exchange_rxn_ids, reactions(models[i]))
        exchange_met_inds = indexin(exchange_met_ids, metabolites(models[i]))
        for (n, (ex_rxn, ex_met)) in enumerate(zip(exchange_rxn_inds, exchange_met_inds)) # each exchange rxn has one exchange met
            isnothing(ex_rxn) && continue
            isnothing(ex_met) && continue
            # connect environmental metabolite with exchange metabolite
            I[kstart] = n_metabolites_metabolic + n
            J[kstart] = ex_rxn + reaction_offset[i]
            V[kstart] = -stoichs[i][ex_met, ex_rxn] # ex is normally negative, make positive
            kstart += 1
        end
    end

    # # build diagonal environmental exchanges
    for i = 1:length(exchange_rxn_ids)
        I[kstart] = n_metabolites_metabolic + i
        J[kstart] = n_reactions_metabolic + i
        V[kstart] = -1.0
        kstart += 1
    end

    if add_biomass_objective
        n_before_biomass_row = n_metabolites_metabolic + length(exchange_met_ids)
        for i = 1:length(models)
            biomass_ind = first(indexin([biomass_ids[i]], reactions(models[i])))
            I[kstart] = i + n_before_biomass_row
            J[kstart] = biomass_ind + reaction_offset[i]
            V[kstart] = 1.0
            kstart += 1
        end
    end

    S = sparse(
        I[1:kstart-1],
        J[1:kstart-1],
        V[1:kstart-1],
        n_metabolites_total,
        n_reactions_total,
    ) # could be that some microbes don't have all the exchanges, hence kstart-1

    _reaction_offsets = cumsum(reaction_lengths)
    _metabolite_offsets = cumsum(metabolite_lengths)
    for i = 1:length(models)
        species = isempty(model_names) ? "species_$(i)_" : model_names[i] * "_"
        tlbs, tubs = bounds(models[i])
        lbs[reaction_offset[i]+1:_reaction_offsets[i]] .= tlbs
        ubs[reaction_offset[i]+1:_reaction_offsets[i]] .= tubs
        rxns[reaction_offset[i]+1:_reaction_offsets[i]] = species .* reactions(models[i])
        mets[metabolite_offset[i]+1:_metabolite_offsets[i]] =
            species .* metabolites(models[i])
    end
    mets[_metabolite_offsets[end]+1:_metabolite_offsets[end]+length(exchange_met_ids)] .=
        exchange_met_ids
    rxns[_reaction_offsets[end]+1:_reaction_offsets[end]+length(exchange_rxn_ids)] .=
        exchange_rxn_ids

    if add_biomass_objective
        rxns[end] = "community_biomass"
        for i = 1:length(models)
            species = isempty(model_names) ? "species_$(i)_" : model_names[i] * "_"
            mets[end-length(biomass_ids)+i] = species .* biomass_ids[i]
        end
    end

    return CoreModel(S, spzeros(size(S, 1)), spzeros(size(S, 2)), lbs, ubs, rxns, mets)
end

"""
    join_with_exchanges(
        ::Type{StandardModel},
        models::Vector{M},
        exchange_rxn_mets::Dict{String,String};
        biomass_ids = [],
        model_names = [],
    )::StandardModel where {M<:MetabolicModel}

A variant of [`join_with_exchanges`](@ref) that returns a `StandardModel`.
"""
function join_with_exchanges(
    ::Type{StandardModel},
    models::Vector{M},
    exchange_rxn_mets::Dict{String,String};
    biomass_ids = [],
    model_names = [],
)::StandardModel where {M<:MetabolicModel}

    community = StandardModel()
    rxns = OrderedDict{String,Reaction}()
    mets = OrderedDict{String,Metabolite}()
    genes = OrderedDict{String,Gene}()
    sizehint!(rxns, sum(n_reactions(m) for m in models) + length(keys(exchange_rxn_mets)))
    sizehint!(
        mets,
        sum(n_metabolites(m) for m in models) + length(values(exchange_rxn_mets)),
    )
    sizehint!(genes, sum(n_genes(m) for m in models))

    for (i, model) in enumerate(models)
        species = isempty(model_names) ? "species_$(i)" : model_names[i] # underscore gets added in add_model_with_exchanges!
        biomass_id = isempty(biomass_ids) ? nothing : biomass_ids[i]
        add_model_with_exchanges!(
            community,
            model,
            exchange_rxn_mets;
            model_name = species,
            biomass_id = biomass_id,
        )
    end

    if !isempty(biomass_ids)
        bm = Dict{String,Float64}() # unassigned
        community.reactions["community_biomass"] =
            Reaction("community_biomass", bm, :forward)
    end

    # Add environmental exchange reactions and metabolites.TODO: add annotation details
    for (rid, mid) in exchange_rxn_mets
        community.reactions[rid] = Reaction(rid)
        community.reactions[rid].metabolites = Dict{String,Float64}(mid => -1.0)
        community.reactions[rid].lb = 0.0
        community.reactions[rid].ub = 0.0
        community.metabolites[mid] = Metabolite(mid)
        community.metabolites[mid].id = mid
    end

    return community
end

"""
    add_model_with_exchanges(
        community::CoreModel,
        model::MetabolicModel,
        exchange_rxn_mets::Dict{String,String};
        model_name = "unknown_species",
        biomass_id = nothing,
    )

Add `model` to `community`, which is a pre-existing community model with exchange reactions
and metabolites in the dictionary `exchange_rxn_mets`. The `model_name` is appended to each
reaction and metabolite, see [`join_with_exchanges`](@ref). If `biomass_id` is specified
then a biomass metabolite for `model` is also added to the resulting model. The column
corresponding to the `biomass_id` reaction then produces this new biomass metabolite with
unit coefficient. The exchange reactions and metabolites in `exchange_rxn_mets` must already
exist in `community`. Always returns a new community model because it is more efficient than
resizing all the matrices.

No in-place variant for `CoreModel`s exists yet.

# Example
```
community = add_model_with_exchanges(community,
    model,
    exchange_rxn_mets;
    model_name="species_2",
    biomass_id="BIOMASS_Ecoli_core_w_GAM")
```
"""
function add_model_with_exchanges(
    community::CoreModel,
    model::MetabolicModel,
    exchange_rxn_mets::Dict{String,String};
    model_name = "unknown_species",
    biomass_id = nothing,
)

    exchange_rxn_ids = keys(exchange_rxn_mets)
    exchange_met_ids = values(exchange_rxn_mets)
    exchange_met_community_inds = indexin(exchange_met_ids, metabolites(community))
    exchange_rxn_community_inds = indexin(exchange_rxn_ids, reactions(community))
    if any(isnothing.(exchange_met_community_inds)) ||
       any(isnothing.(exchange_rxn_community_inds))
        throw(
            DomainError(
                "exchange metabolite/reaction not found.",
                "Exchange metabolites/reactions must already be contained in the community model.",
            ),
        )
    end

    n_cmodel_rows, n_cmodel_cols = size(stoichiometry(community))
    n_model_rows, n_model_cols = size(stoichiometry(model))
    # A note on the variable names here.Suppose M is some sparse matrix, then I
    # = row indices, J = column indices and V = values at the associated
    # indices. So I[a] = i, J[a]=j and then M[i,j] = V[a]
    Iadd, Jadd, Vadd = findnz(stoichiometry(model))

    # shift to fit into community
    Iadd .+= n_cmodel_rows
    Jadd .+= n_cmodel_cols

    # when adding a single model not that many reactions, push! okay?
    exchange_rxn_model_idxs = indexin(exchange_rxn_ids, reactions(model))
    for i = 1:length(exchange_rxn_ids)
        isnothing(exchange_rxn_model_idxs[i]) && continue
        push!(Iadd, exchange_met_community_inds[i]) # already present ex met in community model
        push!(Jadd, n_cmodel_cols + exchange_rxn_model_idxs[i]) # new column hence the offset
        push!(Vadd, 1.0)
    end

    biomass_met = 0.0
    if biomass_id != "" # add biomass metabolite
        biomass_rxn = first(indexin([biomass_id], reactions(model)))
        push!(Iadd, n_model_rows + n_cmodel_rows + 1)
        push!(Jadd, biomass_rxn + n_cmodel_cols)
        push!(Vadd, 1.0)
        biomass_met = 1
    end

    n_metabolites_total = n_model_rows + n_cmodel_rows + biomass_met
    n_reactions_total = n_cmodel_cols + n_model_cols

    I, J, V = findnz(stoichiometry(community))
    I = [I; Iadd]
    J = [J; Jadd]
    V = [V; Vadd]
    S = sparse(I, J, V, n_metabolites_total, n_reactions_total)

    # A note on the variables here. The bounds are vectors of upper and lower
    # bounds for each reaction. So lbs = [lb_1, lb_2, lb_i, ...], ubs = [ub_1,
    # ub_2, ub_i, ...] for reaction i. See the bounds function for more
    # information
    lbsadd, ubsadd = bounds(model)
    lbs, ubs = bounds(community)
    lbs = [lbs; lbsadd]
    ubs = [ubs; ubsadd]

    rxnsadd = "$(model_name)_" .* reactions(model)
    if !isnothing(biomass_id)
        metsadd = ["$(model_name)_" .* metabolites(model); "$(model_name)_" * biomass_id]
    else
        metsadd = "$(model_name)_" .* metabolites(model)
    end
    rxns = [reactions(community); rxnsadd]
    mets = [metabolites(community); metsadd]

    # adds to the original community data, could possibly reset?
    I, V = findnz(balance(community))
    b = sparsevec(I, V, n_metabolites_total)
    I, V = findnz(objective(community))
    c = sparsevec(I, V, n_reactions_total)

    return CoreModel(S, b, c, lbs, ubs, rxns, mets)
end

"""
    add_model_with_exchanges!(
        community::StandardModel,
        model::MetabolicModel,
        exchange_rxn_mets::Dict{String,String};
        model_name = "unknown_species",
        biomass_id = nothing,
    )

The `StandardModel` variant of [`add_model_with_exchanges`](@ref), but is in-place.
"""
function add_model_with_exchanges!(
    community::StandardModel,
    model::MetabolicModel,
    exchange_rxn_mets::Dict{String,String};
    model_name = "unknown_species",
    biomass_id = nothing,
)
    stdm = model isa StandardModel ? deepcopy(model) : convert(StandardModel, model)
    model_name = model_name * "_"

    for met in values(stdm.metabolites)
        met.id = model_name * met.id
        community.metabolites[met.id] = met
    end

    for rxn in values(stdm.reactions)
        # rename reaction string
        rxn.metabolites = Dict(model_name * k => v for (k, v) in rxn.metabolites)
        # change direction of exchange
        if rxn.id in keys(exchange_rxn_mets)
            rxn.metabolites[model_name*exchange_rxn_mets[rxn.id]] = 1.0 # exchanges should be negative originally. TODO: test if they are?
        end
        # add biomass metabolite if applicable
        if rxn.id == biomass_id
            rxn.metabolites[model_name*rxn.id] = 1.0 # produces one biomass
            community.metabolites[model_name*rxn.id] = Metabolite(model_name * rxn.id)
        end
        # add environmental connection
        if rxn.id in keys(exchange_rxn_mets)
            rxn.metabolites[exchange_rxn_mets[rxn.id]] = -1.0
        end
        # reset objective
        rxn.objective_coefficient = 0.0
        rxn.id = model_name * rxn.id
        # add to community model
        community.reactions[rxn.id] = rxn
    end

    for gene in values(stdm.genes)
        gene.id = model_name * gene.id
        community.genes[gene.id] = gene
    end

    return nothing
end

"""
    add_model_with_exchanges(
        community::StandardModel,
        model::MetabolicModel,
        exchange_rxn_mets::Dict{String,String};
        model_name = "unknown_species",
        biomass_id = nothing,
    )

The `StandardModel` variant of [`add_model_with_exchanges`](@ref). Makes a deepcopy of
`community` and calls the inplace variant of this function on that copy.
"""
function add_model_with_exchanges(
    community::StandardModel,
    model::MetabolicModel,
    exchange_rxn_mets::Dict{String,String};
    model_name = "unknown_species",
    biomass_id = nothing,
)
    new_comm = deepcopy(community)
    add_model_with_exchanges!(
        new_comm,
        model,
        exchange_rxn_mets;
        model_name = model_name,
        biomass_id = biomass_id,
    )
    return new_comm
end


text:
command_pool = CommandPool(device, 0)
buffer = first(unwrap(allocate_command_buffers(device, CommandBufferAllocateInfo(command_pool, COMMAND_BUFFER_LEVEL_PRIMARY, 1))))


text:
struct DevTools
    ui::Bool
    props_check::Bool
    serve_dev_bundles::Bool
    hot_reload::Bool
    hot_reload_interval::Float64
    hot_reload_watch_interval::Float64
    hot_reload_max_retry::Int
    silence_routes_logging::Bool
    prune_errors::Bool
    function DevTools(debug = false; ui = nothing,
            props_check = nothing,
            serve_dev_bundles = nothing,
            hot_reload = nothing,
            hot_reload_interval = nothing,
            hot_reload_watch_interval = nothing,
            hot_reload_max_retry = nothing,
            silence_routes_logging = nothing,
            prune_errors = nothing)

        @env_default!(ui, Bool, debug)
        @env_default!(props_check, Bool, debug)
        @env_default!(serve_dev_bundles, Bool, debug)
        @env_default!(hot_reload, Bool, debug)
        @env_default!(silence_routes_logging, Bool, debug)
        @env_default!(prune_errors, Bool, debug)

        @env_default!(hot_reload_interval, Float64, 3.)
        @env_default!(hot_reload_watch_interval, Float64, 0.5)
        @env_default!(hot_reload_max_retry, Int, 8)

        return new(
            ui,
            props_check,
            serve_dev_bundles,
            hot_reload,
            hot_reload_interval,
            hot_reload_watch_interval,
            hot_reload_max_retry,
            silence_routes_logging,
            prune_errors
        )
    end
end


text:
# --------------------------------------------------------
# Cheng Cao, Berkeley 2019
# --------------------------------------------------------
# Translated and extended from C
# Original Author: Andy Zeng, Princeton University, 2016
# --------------------------------------------------------

using CUDAdrv, CUDAnative, CuArrays
using DelimitedFiles, Printf, FileIO, Images
using MeshCat
using CoordinateTransformations
using GeometryTypes: GeometryTypes, Vec, Point, Point3f0
using Colors: RGBA, RGB

@enum occupancy::UInt8 begin
    unknown = 0
    occupied = 1
    free = 255
end

function Integrate(
    cam_K, cam2base, cam2world,
    depth_im, im_height, im_width,
    voxel_grid_dim_x, voxel_grid_dim_y, voxel_grid_dim_z,
    voxel_grid_origin_x, voxel_grid_origin_y, voxel_grid_origin_z,
    voxel_size, trunc_margin,
    voxel_grid_TSDF, voxel_grid_weight, voxel_grid_occupancy
)
    pt_grid_z = blockIdx().x
    pt_grid_y = threadIdx().x

    for pt_grid_x::Int = 1:voxel_grid_dim_x
        pt_base_x = voxel_grid_origin_x + pt_grid_x * voxel_size
        pt_base_y = voxel_grid_origin_y + pt_grid_y * voxel_size
        pt_base_z = voxel_grid_origin_z + pt_grid_z * voxel_size

        tmp_pt_1 = pt_base_x - cam2base[1, 4]
        tmp_pt_2 = pt_base_x - cam2base[2, 4]
        tmp_pt_3 = pt_base_x - cam2base[3, 4]

        pt_cam_x = cam2base[1, 1] * tmp_pt_1 + cam2base[2, 1] * tmp_pt_2 + cam2base[3, 1] * tmp_pt_3
        pt_cam_y = cam2base[1, 2] * tmp_pt_1 + cam2base[2, 2] * tmp_pt_2 + cam2base[3, 2] * tmp_pt_3
        pt_cam_z = cam2base[1, 3] * tmp_pt_1 + cam2base[2, 3] * tmp_pt_2 + cam2base[3, 3] * tmp_pt_3

        if pt_cam_z <= 0
            continue
        end

        pt_pix_x::Int = round(cam_K[1, 1] * (pt_cam_x / pt_cam_z) + cam_K[1, 3])
        pt_pix_y::Int = round(cam_K[2, 2] * (pt_cam_x / pt_cam_z) + cam_K[2, 3])
        if pt_pix_x < 0 || pt_pix_x >= im_width || pt_pix_y < 0 || pt_pix_y >= im_height
            continue
        end

        depth_val = depth_im[pt_pix_y, pt_pix_x]

        diff = depth_val - pt_cam_z

        if diff > 0 && voxel_grid_occupancy[pt_grid_x, pt_grid_y, pt_grid_z] != occupied
            voxel_grid_occupancy[pt_grid_x, pt_grid_y, pt_grid_z] = free
        end

        if depth_val <= 0 || depth_val > 6
            continue
        end

        if diff <= -trunc_margin
            continue
        end

        dist = min(1.0, diff / trunc_margin)
        weight_old = voxel_grid_weight[pt_grid_x, pt_grid_y, pt_grid_z]
        weight_new = weight_old + 1.0
        voxel_grid_weight[pt_grid_x, pt_grid_y, pt_grid_z] = weight_new
        tsdf = (voxel_grid_TSDF[pt_grid_x, pt_grid_y, pt_grid_z] * weight_old + dist) / weight_new
        voxel_grid_TSDF[pt_grid_x, pt_grid_y, pt_grid_z] = tsdf

        if abs(tsdf) < voxel_size / trunc_margin && weight_new > 0.0
            voxel_grid_occupancy[pt_grid_x, pt_grid_y, pt_grid_z] = occupied
        end
    end

    return
end

function main(args)
    # Location of camera intrinsic file
    cam_K_file = "data/camera-intrinsics.txt";

    # Location of folder containing RGB-D frames and camera pose files
    data_path = "data/rgbd-frames";
    base_frame_idx::Int = 150;
    first_frame_idx::Int = 150;
    num_frames::Int = 50;

    cam_K = Array{Float32, 2}(undef, 3, 3)
    base2world = Array{Float32, 2}(undef, 4, 4)
    cam2base = Array{Float32, 2}(undef, 4, 4)
    cam2world = Array{Float32, 2}(undef, 4, 4)
    cam2world_inv = Array{Float32, 2}(undef, 4, 4)
    im_width::Int = 640;
    im_height::Int = 480;
    depth_im = Array{Float32, 2}(undef, im_width, im_height);

    # Voxel grid parameters (change these to change voxel grid resolution, etc.)
    voxel_grid_origin = Float32[-1.5 -1.5 0.5] # Location of voxel grid origin in base frame camera coordinates
    voxel_size::Float32 = 0.006
    trunc_margin::Float32 = voxel_size * 5.0
    voxel_grid_dim = [500 500 500]

    # Manual parameters
    if length(args) > 1
        cam_K = readdlm(args[1], ' ', Float32)
        data_path = args[2];
        base_frame_idx = parse(Int, args[3]);
        first_frame_idx = parse(Int, args[4]);
        num_frames = parse(Int, args[5]);
        voxel_grid_origin = [parse(Float32, args[6]) parse(Float32, args[7]) parse(Float32, args[8])]
        voxel_size = parse(Float32, args[9]);
        trunc_margin = parse(Float32, args[10]);
    end

    @show cam_K data_path base_frame_idx first_frame_idx num_frames voxel_grid_origin voxel_size trunc_margin

    base2world = Float32[1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 1]
    base2world_inv = inv(base2world)

    voxel_grid_TSDF = fill(1.0, (voxel_grid_dim[1], voxel_grid_dim[2], voxel_grid_dim[3]))
    voxel_grid_weight = fill(0.0, (voxel_grid_dim[1], voxel_grid_dim[2], voxel_grid_dim[3]))
    voxel_grid_occupancy = fill(unknown, (voxel_grid_dim[1], voxel_grid_dim[2], voxel_grid_dim[3]))

    gpu_voxel_grid_TSDF = CuArray(voxel_grid_TSDF)
    gpu_voxel_grid_weight = CuArray(voxel_grid_weight)
    gpu_voxel_grid_occupancy = CuArray(voxel_grid_occupancy)

    gpu_cam_K = CuArray(cam_K)
    gpu_cam2base = CuArray(cam2base)
    gpu_cam2world = CuArray(cam2world)
    gpu_voxel_grid_dim = CuArray(voxel_grid_dim)
    gpu_voxel_grid_origin = CuArray(voxel_grid_origin)

    for frame_idx::Int = first_frame_idx:first_frame_idx + num_frames
        depth_im_file = @sprintf "%s/frame-%06d.depth.png" data_path frame_idx
        depth_im = convert(Array{Float32, 2}, Gray.(load(depth_im_file)))

        cam2world_file = @sprintf "%s/frame-%06d.pose.txt" data_path frame_idx
        cam2world = readdlm(cam2world_file, ' ', Float32)

        base2world_inv = cam2world * cam2base;

        cam2world_inv = inv(cam2world)

        gpu_cam2base = CuArray(cam2base)
        gpu_cam2world = CuArray(cam2world)
        gpu_depth_im = CuArray(depth_im)

        println("Fusing: ", depth_im_file)

        @cuda blocks=voxel_grid_dim[3] threads=voxel_grid_dim[2] Integrate(
            gpu_cam_K, gpu_cam2base, gpu_cam2world,
            gpu_depth_im, im_height, im_width,
            gpu_voxel_grid_dim[1], gpu_voxel_grid_dim[2], gpu_voxel_grid_dim[3],
            gpu_voxel_grid_origin[1], gpu_voxel_grid_origin[2], gpu_voxel_grid_origin[3],
            voxel_size, trunc_margin,
            gpu_voxel_grid_TSDF, gpu_voxel_grid_weight, gpu_voxel_grid_occupancy
        )
    end

    voxel_grid_TSDF = Array(gpu_voxel_grid_TSDF)
    voxel_grid_occupancy = Array(gpu_voxel_grid_occupancy)

    write("tsdf.bin", voxel_grid_occupancy)

    vis = Visualizer()
    open(vis)

    verts = Point3f0[]
    colors = RGBA[]

    for x = 1:voxel_grid_dim[1]
        for y = 1:voxel_grid_dim[2]
            for z = 1:voxel_grid_dim[3]
                #if voxel_grid_occupancy[x, y, z] == occupied
                    push!(verts, Point3f0(
                        x * voxel_size + voxel_grid_origin[1],
                        y * voxel_size + voxel_grid_origin[2],
                        z * voxel_size + voxel_grid_origin[3]
                    ))
                    push!(colors, RGBA(1.0, 1.0, 1.0, voxel_grid_TSDF[x, y, z]))
                #end
            end
        end
        println(x)
    end

    setobject!(vis, PointCloud(verts, colors))

    return readline()
end

main(ARGS)


text:
# -----------------------------------------------------
# This script gives the methods for maniupulating the
# Cartesian indices corresponding to multilevels that
# are either kept or cut in various schemes of
# "sparsification"
# -----------------------------------------------------

# Efficiency Criticality: LOW
# This is not evaluated very often

# Accuracy Criticality: N/A
# No float manipulation

# Gives the appropriate boolean cutoff corresponding
# to a given scheme, e.g. sparse basis, full basis,
# and in the future possibly the energy basis of
# Bungartz and Griebel

function cutoff(::Val{:sparse}, x::CartesianIndex{D}, n::Int) where {D}
    sum(x.I) > n+D
end

function cutoff(::Val{:full}, x::CartesianIndex{D}, n::Int) where {D}
    false
end

# Same functions but for anisotropic grids
# all grids where the multi-indices are bigger than the current level get cutoff
# index1 <=! n(1)+1, index2 <=! n(2)+1
function cutoff(::Val{:full}, x::CartesianIndex{D}, n::Vector{Int}) where {D}
    test = 0
    if sum(x.I) <= sum(n) + length(n) # eigentlich <=sum(n+1) aber n Vector
        for d in 1:D
            if x.I[d] > n[d]+1 # n+1 weil CartesianIndex bei 1 startet obwohl unsere theoretischen Levels bei 0 starten
                test += 1
            end
        end
        if test == 0
            return false
        end
    end
    return true
end

# cutoff for anisotropic grids (n is vector) only needed for full grids
# function cutoff(::Val{:sparse}, x::CartesianIndex{D}, n::Vector{Int}) where D
#     false
# end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate((a+a*csc(x))^-1*csc(x)^5, x) == :(-4*a^-1*cot(x)+-4//3*a^-1*cot(x)^3+3//2*a^-1*arctanh(cos(x))+(a+a*csc(x))^-1*csc(x)^3*cot(x)+3//2*a^-1*cot(x)*csc(x))
@test integrate((a+a*csc(x))^-1*csc(x)^4, x) == :(2*a^-1*cot(x)+-3//2*a^-1*arctanh(cos(x))+(a+a*csc(x))^-1*csc(x)^2*cot(x)+-3//2*a^-1*cot(x)*csc(x))
@test integrate((a+a*csc(x))^-1*csc(x)^3, x) == :(a^-1*arctanh(cos(x))+-1*a^-1*cot(x)+-1*(a+a*csc(x))^-1*cot(x))
@test integrate((a+a*csc(x))^-1*csc(x)^2, x) == :((a+a*csc(x))^-1*cot(x)+-1*a^-1*arctanh(cos(x)))
@test integrate((a+a*csc(x))^-1*csc(x), x) == :(-1*(a+a*csc(x))^-1*cot(x))
@test integrate((a+a*csc(c+d*x))^-1, x) == :(x*a^-1+d^-1*(a+a*csc(c+d*x))^-1*cot(c+d*x))
@test integrate((a+a*csc(x))^-1*sin(x), x) == :((a+a*csc(x))^-1*cos(x)+-1*x*a^-1+-2*a^-1*cos(x))
@test integrate((a+a*csc(x))^-1*sin(x)^2, x) == :(2*a^-1*cos(x)+3//2*x*a^-1+(a+a*csc(x))^-1*cos(x)*sin(x)+-3//2*a^-1*cos(x)*sin(x))
@test integrate((a+a*csc(x))^-1*sin(x)^3, x) == :(-4*a^-1*cos(x)+-3//2*x*a^-1+4//3*a^-1*cos(x)^3+(a+a*csc(x))^-1*sin(x)^2*cos(x)+3//2*a^-1*cos(x)*sin(x))
@test integrate((a+a*csc(x))^-1*sin(x)^4, x) == :(4*a^-1*cos(x)+-4//3*a^-1*cos(x)^3+15//8*x*a^-1+(a+a*csc(x))^-1*sin(x)^3*cos(x)+-15//8*a^-1*cos(x)*sin(x)+-5//4*a^-1*sin(x)^3*cos(x))
@test integrate((a+a*csc(c+d*x))^-2, x) == :(x*a^-2+1//3*d^-1*(a+a*csc(c+d*x))^-2*cot(c+d*x)+4//3*a^-2*d^-1*(1+csc(c+d*x))^-1*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^-3, x) == :(x*a^-3+1//5*d^-1*(a+a*csc(c+d*x))^-3*cot(c+d*x)+22//15*d^-1*(a^3+a^3*csc(c+d*x))^-1*cot(c+d*x)+7//15*a^-1*d^-1*(a+a*csc(c+d*x))^-2*cot(c+d*x))
@test integrate((a+a*csc(x))^5//2, x) == :(-2*a^5//2*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x))+-14//3*a^3*(a+a*csc(x))^-1//2*cot(x)+-2//3*a^2*(a+a*csc(x))^(1/2)*cot(x))
@test integrate((a+a*csc(x))^3//2, x) == :(-2*a^3//2*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x))+-2*a^2*(a+a*csc(x))^-1//2*cot(x))
@test integrate((a+a*csc(x))^(1/2), x) == :(-2*a^(1/2)*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x)))
@test integrate((a+a*csc(x))^-1//2, x) == :(-2*a^-1//2*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x))+2^(1/2)*a^-1//2*arctan((1/2)*2^(1/2)*a^(1/2)*(a+a*csc(x))^-1//2*cot(x)))
@test integrate((a+a*csc(x))^-3//2, x) == :((1/2)*(a+a*csc(x))^-3//2*cot(x)+-2*a^-3//2*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x))+5//4*2^(1/2)*a^-3//2*arctan((1/2)*2^(1/2)*a^(1/2)*(a+a*csc(x))^-1//2*cot(x)))
@test integrate((a+a*csc(x))^-5//2, x) == :(-2*a^-5//2*arctan(a^(1/2)*(a+a*csc(x))^-1//2*cot(x))+1//4*(a+a*csc(x))^-5//2*cot(x)+11//16*a^-1*(a+a*csc(x))^-3//2*cot(x)+43//32*2^(1/2)*a^-5//2*arctan((1/2)*2^(1/2)*a^(1/2)*(a+a*csc(x))^-1//2*cot(x)))
@test integrate((a+a*csc(e+f*x))^(1/2)*csc(e+f*x)^(1/2), x) == :(-2*a^(1/2)*f^-1*arcsinh(a^(1/2)*(a+a*csc(e+f*x))^-1//2*cot(e+f*x)))
@test integrate((-1*csc(e+f*x))^(1/2)*(a+-1*a*csc(e+f*x))^(1/2), x) == :(-2*a^(1/2)*f^-1*arcsinh(a^(1/2)*(a+-1*a*csc(e+f*x))^-1//2*cot(e+f*x)))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^4//3, x) == :(-6//5*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^4//3*cos(c+d*x)+-4//5*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^1//3, x) == :(-2*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^-2//3, x) == :(-3//2*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^1//3*cos(c+d*x)+-1//2*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^5//3, x) == :(-6//7*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^5//3*cos(c+d*x)+24//7*a*d^-1*(a+a*csc(c+d*x))^-1//2*(1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*cot(c+d*x)+-12//7*3^1//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+-1*3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.E(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x)+8//7*2^(1/2)*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^2//3, x) == :(6*a*d^-1*(a+a*csc(c+d*x))^-1//2*(1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*cot(c+d*x)+-3*3^1//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+-1*3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.E(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x)+2*2^(1/2)*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^-1//3, x) == :(-3*a*d^-1*(a+a*csc(c+d*x))^-1//2*(1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*cot(c+d*x)+-3*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^2//3*cos(c+d*x)+-1*2^(1/2)*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x)+3//2*3^1//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+-1*3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.E(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+a*csc(c+d*x))^(1/2)*csc(c+d*x)^-4//3, x) == :(-15//8*a*d^-1*(a+a*csc(c+d*x))^-1//2*(1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*cot(c+d*x)+-15//8*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^2//3*cos(c+d*x)+-3//4*a*d^-1*(a+a*csc(c+d*x))^-1//2*csc(c+d*x)^-1//3*cos(c+d*x)+-5//8*2^(1/2)*3^3//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.F(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x)+15//16*3^1//4*a^2*d^-1*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+-1*csc(c+d*x)^1//3))^-1//2*((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-2*(1+csc(c+d*x)^1//3+csc(c+d*x)^2//3))^(1/2)*(2+-1*3^(1/2))^(1/2)*(a+a*csc(c+d*x))^-1//2*(a+-1*a*csc(c+d*x))^-1*(1+-1*csc(c+d*x)^1//3)*Elliptic.E(arcsin((1+3^(1/2)+-1*csc(c+d*x)^1//3)^-1*(1+-1*3^(1/2)+-1*csc(c+d*x)^1//3)),-7+-4*3^(1/2))*cot(c+d*x))
@test integrate((a+b*csc(c+d*x))^4, x) == :(x*a^4+-1//3*b^2*d^-1*(a+b*csc(c+d*x))^2*cot(c+d*x)+-1//3*b^2*d^-1*(2*b^2+17*a^2)*cot(c+d*x)+-2*a*b*d^-1*(b^2+2*a^2)*arctanh(cos(c+d*x))+-4//3*a*b^3*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((a+b*csc(c+d*x))^3, x) == :(x*a^3+-5//2*a*b^2*d^-1*cot(c+d*x)+-1//2*b*d^-1*(b^2+6*a^2)*arctanh(cos(c+d*x))+-1//2*b^2*d^-1*(a+b*csc(c+d*x))*cot(c+d*x))
@test integrate((a+b*csc(c+d*x))^2, x) == :(x*a^2+-1*b^2*d^-1*cot(c+d*x)+-2*a*b*d^-1*arctanh(cos(c+d*x)))
@test integrate((a+b*csc(x))^-1*csc(x)^5, x) == :(-1//3*b^-1*csc(x)^2*cot(x)+-1//3*b^-3*(2*b^2+3*a^2)*cot(x)+(1/2)*a*b^-4*(b^2+2*a^2)*arctanh(cos(x))+(1/2)*a*b^-2*cot(x)*csc(x)+-2*a^4*b^-4*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*csc(x)^4, x) == :(a*b^-2*cot(x)+-1//2*b^-1*cot(x)*csc(x)+-1//2*b^-3*(b^2+2*a^2)*arctanh(cos(x))+2*a^3*b^-3*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*csc(x)^3, x) == :(-1*b^-1*cot(x)+a*b^-2*arctanh(cos(x))+-2*a^2*b^-2*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*csc(x)^2, x) == :(-1*b^-1*arctanh(cos(x))+2*a*b^-1*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*csc(x), x) == :(-2*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(c+d*x))^-1, x) == :(x*a^-1+2*b*a^-1*d^-1*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*csc(x))^-1*sin(x), x) == :(-1*a^-1*cos(x)+-1*b*x*a^-2+-2*a^-2*b^2*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*sin(x)^2, x) == :(b*a^-2*cos(x)+(1/2)*x*a^-3*(a^2+2*b^2)+-1//2*a^-1*cos(x)*sin(x)+2*a^-3*b^3*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x))))
@test integrate((a+b*csc(x))^-1*sin(x)^3, x) == :(-1//3*a^-1*sin(x)^2*cos(x)+-1//3*a^-3*(2*a^2+3*b^2)*cos(x)+(1/2)*b*a^-2*cos(x)*sin(x)+-2*a^-4*b^4*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x)))+-1//2*b*x*a^-4*(a^2+2*b^2))
@test integrate((a+b*csc(x))^-1*sin(x)^4, x) == :(-1//4*a^-1*sin(x)^3*cos(x)+1//8*x*a^-5*(3*a^4+8*b^4+4*a^2*b^2)+2*a^-5*b^5*(a^2+-1*b^2)^-1//2*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*x)))+-1//8*a^-3*(3*a^2+4*b^2)*cos(x)*sin(x)+1//3*b*a^-4*(2*a^2+3*b^2)*cos(x)+1//3*b*a^-2*sin(x)^2*cos(x))
@test integrate((a+b*csc(c+d*x))^-2, x) == :(x*a^-2+-1*a^-1*b^2*d^-1*(a+b*csc(c+d*x))^-1*(a^2+-1*b^2)^-1*cot(c+d*x)+2*b*a^-2*d^-1*(a^2+-1*b^2)^-3//2*(-1*b^2+2*a^2)*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*c+(1/2)*d*x))))
@test integrate((a+b*csc(c+d*x))^-3, x) == :(x*a^-3+b*a^-3*d^-1*(a^2+-1*b^2)^-5//2*(2*b^4+6*a^4+-5*a^2*b^2)*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*c+(1/2)*d*x)))+-1//2*a^-1*b^2*d^-1*(a+b*csc(c+d*x))^-2*(a^2+-1*b^2)^-1*cot(c+d*x)+-1//2*a^-2*b^2*d^-1*(a+b*csc(c+d*x))^-1*(a^2+-1*b^2)^-2*(-2*b^2+5*a^2)*cot(c+d*x))
@test integrate((a+b*csc(c+d*x))^-4, x) == :(x*a^-4+b*a^-4*d^-1*(a^2+-1*b^2)^-7//2*(-2*b^6+8*a^6+-8*a^4*b^2+7*a^2*b^4)*arctanh((a^2+-1*b^2)^-1//2*(a+b*tan((1/2)*c+(1/2)*d*x)))+-1//3*a^-1*b^2*d^-1*(a+b*csc(c+d*x))^-3*(a^2+-1*b^2)^-1*cot(c+d*x)+-1//6*a^-3*b^2*d^-1*(a+b*csc(c+d*x))^-1*(a^2+-1*b^2)^-3*(6*b^4+26*a^4+-17*a^2*b^2)*cot(c+d*x)+-1//6*a^-2*b^2*d^-1*(a+b*csc(c+d*x))^-2*(a^2+-1*b^2)^-2*(-3*b^2+8*a^2)*cot(c+d*x))
@test integrate((3+5*csc(c+d*x))^-1, x) == :(-1//12*x+-5//6*d^-1*arctan((3+sin(c+d*x))^-1*cos(c+d*x)))
@test integrate((5+3*csc(c+d*x))^-1, x) == :(1//5*x+-3//20*d^-1*log(3*sin((1/2)*c+(1/2)*d*x)+cos((1/2)*c+(1/2)*d*x))+3//20*d^-1*log(3*cos((1/2)*c+(1/2)*d*x)+sin((1/2)*c+(1/2)*d*x)))

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.1.2 (d csc)^n (a+b csc)^m.jl


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Structures of SatelliteToolbox.jl.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

include("./atmospheric_models.jl")
include("./gravity_models.jl")
include("./iers_eop.jl")
include("./orbit.jl")
include("./propagators.jl")
include("./reference_frames.jl")
include("./ellipsoid.jl")


text:
module Custom

using Flux
using BSON:@load
using DataFrames
using Pandas
using Pkg

export Flux, load_model, transform, score
# export loadModel, score

function load_model(code_dir)
    @load "$(code_dir)/mymodel.bson" model
    return model
end

function transform(df, model)
    try
        df = DataFrames.select(df, Not(:MEDV))
    catch
        nothing
    end
    return df
end

function score(df, model; kwargs...)
    # X = convert(Matrix, df)
    X = Matrix(df)
    yhat = model(X')
    yhat = yhat'
    return Pandas.DataFrame(yhat, columns = ["Predictions"])
end

end

max_stars_repo_name: amperie/user-models

max_issues_repo_name: amperie/user-models

max_forks_repo_name: amperie/user-models


text:
using FLoops
using FoldsPreview
using UnicodePlots

function online_histogram(;
    N = 100_000_000_000,
    interval = 10,
    ex = ThreadedEx(),
    nbins = 1000,
)
    function on_preview(((hist,), (nsamples,)))
    # function on_preview(((hist,), (nsamples,), (tids,)))
        # display(heatmap(log10.(hist .+ 1)))
        display(heatmap(hist))
        @show nsamples
        # @show sort(collect(tids))
    end
    with_preview(on_preview, ex, interval) do previewer
        @floop ex for n in previewer(1:N)
            x = randn() * 0.3 + rand(-1:1)
            y = randn() * 0.3 + rand(-1:1)
            indexfor(v) = floor(Int, (v + 1.5) / 3 * nbins)
            i = indexfor(x)
            j = indexfor(y)
            0 < i <= nbins || continue
            0 < j <= nbins || continue
            sample = (i, j)
            @reduce() do (hist = zeros(Int, nbins, nbins); sample)
                if sample isa Tuple  # sequential case
                    hist[i, j] += 1
                else
                    sample::Matrix{Int}  # otherwise, we are merging
                    hist .+= sample
                end
            end
            @reduce(nsamples = 0 + 1)
            # @reduce  tids = union!(Set{Int}(), Threads.threadid())
        end
        return hist
    end
end


text:
import Base: convert, ==

function convert(::Type{DFGFactorSummary}, f::DFGFactor)
    return DFGFactorSummary(f.label, deepcopy(f.tags), f._internalId, deepcopy(f._variableOrderSymbols))
end

function packFactor(dfg::G, f::DFGFactor)::Dict{String, Any} where G <: AbstractDFG
    # Construct the properties to save
    props = Dict{String, Any}()
    props["label"] = string(f.label)
    props["tags"] = JSON2.write(f.tags)
    # Pack the node data
    fnctype = f.data.fnc.usrfnc!
    try
        packtype = getfield(_getmodule(fnctype), Symbol("Packed$(_getname(fnctype))"))
        packed = convert(PackedFunctionNodeData{packtype}, f.data)
        props["data"] = JSON2.write(packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while packing '$(f.label)' as '$fnctype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end
    # Include the type
    props["fnctype"] = String(_getname(fnctype))
    props["_variableOrderSymbols"] = JSON2.write(f._variableOrderSymbols)
    props["backendset"] = f.backendset
    props["ready"] = f.ready

    return props
end

function unpackFactor(dfg::G, packedProps::Dict{String, Any}, iifModule)::DFGFactor where G <: AbstractDFG
    label = packedProps["label"]
    tags = JSON2.read(packedProps["tags"], Vector{Symbol})

    data = packedProps["data"]
    @debug "Decoding $label..."
    datatype = packedProps["fnctype"]
    packtype = getfield(Main, Symbol("Packed"*datatype))
    packed = nothing
    fullFactor = nothing
    try
        packed = JSON2.read(data, GenericFunctionNodeData{packtype,String})
        fullFactor = iifModule.decodePackedType(dfg, packed)
    catch ex
        io = IOBuffer()
        showerror(io, ex, catch_backtrace())
        err = String(take!(io))
        msg = "Error while unpacking '$label' as '$datatype', please check the unpacking/packing converters for this factor - \r\n$err"
        error(msg)
    end

    # Include the type
    _variableOrderSymbols = JSON2.read(packedProps["_variableOrderSymbols"], Vector{Symbol})
    backendset = packedProps["backendset"]
    ready = packedProps["ready"]

    # Rebuild DFGVariable
    factor = DFGFactor{typeof(fullFactor.fnc), Symbol}(Symbol(label))
    factor.tags = tags
    factor.data = fullFactor
    factor._variableOrderSymbols = _variableOrderSymbols
    factor.ready = ready
    factor.backendset = backendset

    # GUARANTEED never to bite us in the ass in the future...
    # ... TODO: refactor if changed: https://github.com/JuliaRobotics/IncrementalInference.jl/issues/350
    factor.data.fncargvID = deepcopy(_variableOrderSymbols)

    # Note, once inserted, you still need to call IIF.rebuildFactorMetadata!
    return factor
end

# Compare FunctionNodeData
function compare(a::GenericFunctionNodeData{T1,S},b::GenericFunctionNodeData{T2,S}) where {T1, T2, S}
  # TODO -- beef up this comparison to include the gwp
  TP = true
  TP = TP && a.fncargvID == b.fncargvID
  TP = TP && a.eliminated == b.eliminated
  TP = TP && a.potentialused == b.potentialused
  TP = TP && a.edgeIDs == b.edgeIDs
  TP = TP && a.frommodule == b.frommodule
  # TP = TP && typeof(a.fnc) == typeof(b.fnc)
  return TP
end

"""
    $(SIGNATURES)
Equality check for DFGFactor.
"""
function ==(a::DFGFactor, b::DFGFactor)::Bool
    return compareFactor(a, b)
end


text:
using MPIReco

@testset "single- and multi-channel in-memory reconstruction" begin
  bSF = MPIFile("./data/SF_MP")
  b = MPIFile("./data/MP01")
  names = (:color, :x, :y, :z, :time)
  values = (1:1,
	    -27.375u"mm":1.25u"mm":11.375u"mm",
	    -11.375u"mm":1.25u"mm":27.375u"mm",
	    0.0u"mm":1.0u"mm":0.0u"mm",
	    0.0u"ms":0.6528u"ms":0.0u"ms")

  # standard reconstruction
  c1 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, spectralLeakageCorrection=true)
  @test axisnames(c1) == names
  @test axisvalues(c1) == values
  exportImage("./img/Reconstruction1.png", arraydata(data(c1[1,:,:,1,1])))

  # fused lasso
  c2 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=100, solver="fusedlasso",
		      loadasreal=true, λ=[0.01,0.01])
  @test axisnames(c2) == names
  @test axisvalues(c2) == values
  exportImage("./img/Reconstruction2.png", arraydata(data(c2[1,:,:,1,1])))

  # with interpolation
  c3 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, gridsize=[100,100,1],iterations=1)
  @test axisnames(c3) == names
  @test axisvalues(c3) == (values[1], -27.8u"mm":0.4u"mm":11.8u"mm", -11.8u"mm":0.4u"mm":27.8u"mm", values[4:5]...)
  exportImage("./img/Reconstruction3.png", arraydata(data(c3[1,:,:,1,1])))

  # with fov adpation and center shift
  c4 = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, gridsize=[50,50,1],fov=[0.04,0.04,1],
		      center=[0.0,-0.01,0], iterations=1)
  @test axisnames(c4) == names
  # TODO Tobi: does this make sense?
  @test axisvalues(c4) == (values[1], -27.6u"mm":0.8u"mm":11.6u"mm", -11.6u"mm":0.8u"mm":27.6u"mm", 499.5u"mm":1000.0u"mm":499.5u"mm", values[5])
  exportImage("./img/Reconstruction4.png", arraydata(data(c4[1,:,:,1,1])))

  # multi colored reconstruction (deprecated interface)
  c5 = reconstruction([bSF,bSF], b;
		      SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1)
  @test axisnames(c5) == names
  @test axisvalues(c5) == (1:2,values[2:end]...)
  exportImage("./img/Reconstruction5a.png", arraydata(data(c5[1,:,:,1,1])))
  exportImage("./img/Reconstruction5b.png", arraydata(data(c5[2,:,:,1,1])))

  # multi colored reconstruction (new interface)
  c5 = reconstruction(MultiContrastFile([bSF,bSF]), b;
          SNRThresh=5, frames=1, minFreq=80e3,
          recChannels=1:2, iterations=1)
  @test axisnames(c5) == names
  @test axisvalues(c5) == (1:2,values[2:end]...)
  exportImage("./img/Reconstruction5c.png", arraydata(data(c5[1,:,:,1,1])))
  exportImage("./img/Reconstruction5d.png", arraydata(data(c5[2,:,:,1,1])))

  # dict based reco
  r = defaultRecoParams()
  r[:measPath] = filepath(b)
  r[:SFPath] = filepath(bSF)
  r[:frames] =  1
  r[:minFreq] = 80e3
  r[:SNRThresh] = 4
  r[:lambd] = 0.001
  r[:iterations] = 1
  c6 = reconstruction(r)
  @test axisnames(c6) == names
  @test axisvalues(c6) == values
  exportImage("./img/Reconstruction6.png", arraydata(data(c6[1,:,:,1,1])))

  # channel weighting
  c7a = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[1.0,1.0,1.0])

  exportImage("./img/Reconstruction7a.png", arraydata(data(c7a[1,:,:,1,1])))

  c7b = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[1.0,0.001,1.0])

  exportImage("./img/Reconstruction7b.png", arraydata(data(c7b[1,:,:,1,1])))

  c7c = reconstruction(bSF, b; SNRThresh=5, frames=1, minFreq=80e3,
		      recChannels=1:2, iterations=1, λ=0.1,
          weightType=WeightingType.Channel, channelWeights=[0.001,1.0,1.0])

  exportImage("./img/Reconstruction7c.png", arraydata(data(c7c[1,:,:,1,1])))

end


text:
@testset "Single-line Strings can be diffed" begin
    @testset "bπiz -> bπaz" begin
        s1 = "bπiz"
        s2 = "bπaz"
        diff = deepdiff(s1, s2)
        @test before(diff) == s1
        @test after(diff) == s2
        # the indices are assuming a Vector of chars, like you'd get from `collect`
        @test removed(diff) == [3]
        @test added(diff) == [3]
        @test changed(diff) == []
    end
end

@testset "Multi-line Strings can be diffed" begin
    s1 = """differences can
           be hard to find
           in
           multiline
           output"""
    s2 = """differences can
           be hurd to find
           multiline
           output"""
    diff = deepdiff(s1, s2)
    @test before(diff) == s1
    @test after(diff) == s2
    @test removed(diff) == [2, 3]
    @test added(diff) == [2]
    @test changed(diff) == []
end


text:
using JuMP
using CPLEX
using Plots

function Rolling(Tprim::Int64)
  T, M, Emax = 36, 4, 3
  f = [10,30,60,90, 10, 30]
  e = [8,6,4,2, 8, 6]
  h = Array{Int64,1}(ones(T))
  p = Array{Int64,2}(zeros(T,M))
  d = Array{Int64,1}(ones(T))
  for t in 1:T
    d[t] = rand(20:70)
  end

  # Create the model
  model = Model(CPLEX.Optimizer)

  ## Variables
  @variable(model, x[1:T,1:M] >= 0)
  @variable(model, y[1:T,1:M], Bin)
  @variable(model, s[1:T+1] >= 0)

  ## Constraints
  @constraint(model, s[1] == 0)
  @constraint(model, [t in 2:T+1], sum(x[t-1,m] - s[t] + s[t-1] for m in 1:M) == d[t-1])
  @constraint(model, [t in 1:T, m in 1:M], x[t,m] <= y[t,m]*sum(d[tprim] for tprim in t:T))
  @constraint(model, [t in 1:(T-Tprim)], sum((e[m] - Emax)*x[tprim,m] for tprim in t:(t+Tprim), m in 1:M) <= 0)

  ## Objective
  @objective(model, Min, sum(f[m]*y[t,m] for m in 1:M, t in 1:T) + sum(s[t] for t in  1:T+1))

  #resolution
  optimize!(model)

  x_sol = JuMP.getvalue.(model[:x])
  cout_carbone = sum(e[m]*x_sol[t,m] for m in 1:M, t in 1:T)
  println(fout, "cout carbone = ", cout_carbone)
  println(fout, "objective value = ", JuMP.getobjectivevalue.(model))
  return cout_carbone, JuMP.getobjectivevalue.(model)
end

fout = open("resultatsEmax3M4T36.txt", "w")

T = 36
cout_carbone_history = Array{Float64,1}(zeros(T))
cout_history = Array{Float64,1}(zeros(T))

for Tprim in 1:T
  start = time()
  println(fout, "Tprim = ", Tprim)
  cout_carbone, cout = Rolling(Tprim)
  cout_carbone_history[Tprim] = cout_carbone
  cout_history[Tprim] = cout
  println(fout, "execution time = ", time() - start, "s")
end
close(fout)

plotlyjs()
# plot(1:T,cout_carbone_history)
# xlabel!("Tprim")
# ylabel!("cout carbone")

plot(1:T,cout_history)
xlabel!("Tprim")
ylabel!("cout")

max_stars_repo_name: boumpteryx/RODD

max_issues_repo_name: boumpteryx/RODD

max_forks_repo_name: boumpteryx/RODD


text:
"""
    plotld!(ax::Axis, LD::AbstractMatrix; color::AbstractString)

Visualize a symmetric correlation matrix `LD` with the diagonal elements on the x axis.
`color` can be either "red", "green", "blue", or "black".
"""
function plotld!(ax::Axis, LD::AbstractMatrix; color::AbstractString = "blue")
    if any(isnan.(LD))
        throw(ErrorException("LD matrix contains NaNs"))
    end
    n = size(LD, 1)
    m = count(>(1 / 9), LD)
    addx1 = [0, 1, 0, 1]
    addx2 = [1, 1, 0, 1]
    addy = [0, -1, 0, 1]
    ps = Vector{Polygon}(undef, Int((m - n) / 2 + n))
    LDvech = Vector{Float32}(undef, Int((m - n) / 2 + n))
    counter = 1
    for i in 1:n
        for j in i:n
            if LD[i, j] <= 1 / 9
                continue
            else
                ps[counter] = Polygon([Point2f(
                    abs(i - j) * 5 / n + 10 / n * (i - addx2[k]) + 5 / n * addx1[k], # x coords
                    5 - abs(i - j) * 5 / n + 5 / n * addy[k] # y coords
                    ) for k in 1:4])
                LDvech[counter] = LD[i, j]
                counter += 1
            end
        end
    end
    if color == "black"
        poly!(ax, ps, color = LDvech, colorrange = (0, 1),
            colormap = cgrad(:Greys_9, 9, categorical = true), strokewidth = 0)
    elseif color == "green"
        poly!(ax, ps, color = LDvech, colorrange = (0, 1),
            colormap = cgrad(:Greens_9, 9, categorical = true), strokewidth = 0)
    elseif color == "red"
        poly!(ax, ps, color = LDvech, colorrange = (0, 1),
            colormap = cgrad(:Reds_9, 9, categorical = true), strokewidth = 0)
    else
        poly!(ax, ps, color = LDvech, colorrange = (0, 1),
            colormap = cgrad(:Blues_9, 9, categorical = true), strokewidth = 0)
    end
    ax.aspect = DataAspect()
    ax.spinewidth = 0.75
    xlims!(0, 10)
    ylims!(0, 5)
    hidedecorations!(ax)
    hidespines!(ax, :r, :l, :b)
end


text:
# Determine the best parameters to set for FlowSeed and PPR Push

BestParams = zeros(Int64,4,17)
include("Plots/Plot_Helper.jl")

epsis = [0.05; .1; .25; .5]
Sizes = []
# There are 4 seed types:
# 1: 100 nodes from target region
# 2: 1% of target region
# 3: 2% of target region
# 4: 3% of target region
for seedtype = 1:4

    ## Get the F1 scores for different runs of FlowSeed, with both soft and strict penalties
    p1, r1, f1, Sizes, seednodes =  LoadFlowSeedTrain(.05,seedtype)
    p2, r2, f2, Sizes, seednodes =  LoadFlowSeedTrain(.1,seedtype)
    p3, r3, f3, Sizes, seednodes =  LoadFlowSeedTrain(.25,seedtype)
    p4, r4, f4, Sizes, seednodes =  LoadFlowSeedTrain(.5,seedtype)

    FS_f1 = [f1[4,:]'; f2[4,:]'; f3[4,:]'; f4[4,:]']

    # For each region, get the best F1 score, and which tolerance lead to it
    for i = 1:17
        f1b, loc = findmax(FS_f1[:,i],1)
        BestParams[seedtype,i] = loc[1]
    end

end
## Get the best of PPR, with various tolerances
BestParamsPPR = zeros(Int64,4,17)

alpha = .6
for seedtype = 1:4

    prPPR, rePPR, f1PPR, Sizes = LoadPPRTrain(alpha,seedtype)
    PR = f1PPR
    for i = 1:17
        f1b, loc = findmax(PR[:,i],1)
        BestParamsPPR[seedtype,i] = loc[1]
    end

end

Tols = [1e-7; 1e-8; 1e-9; 1e-10; 1e-11]
epsis = [0.05; .1; .25; .5]

using MAT

#matwrite("Training17_Results.mat", Dict( "Tols" => Tols, "epsis" => epsis,
#"BestParamsPPR" => BestParamsPPR, "BestParamsFlowSeed" => BestParams, "Sizes" => Sizes))

max_stars_repo_name: nveldt/FlowSeed

max_issues_repo_name: nveldt/FlowSeed

max_forks_repo_name: nveldt/FlowSeed


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "SDL2_net"
version = v"2.0.1"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://libsdl.org/projects/SDL_net/release/SDL2_net-$(version).tar.gz",
                  "15ce8a7e5a23dafe8177c8df6e6c79b6749a03fff1e8196742d3571657609d21"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SDL*/

if [[ "${target}" == powerpc64le-* ]] || [[ "${target}" == *-freebsd* ]]; then
    # We need to regenerate `configure` for this platforms

    # Include directory with M4 macros
    atomic_patch -p1 ../patches/configure-add-macro-dir.patch
    # Create absolutely important files
    touch AUTHORS NEWS README ChangeLog
    # Regnerate `configure` script
    autoreconf -fiv
fi

mkdir build && cd build
../configure --prefix=${prefix} --build=${MACHTYPE} --host=${target}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["libSDL2_net", "SDL2_net"], :libSDL2_net)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="SDL2_jll", uuid="ab825dc5-c88e-5901-9575-1e5e20358fcf"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
@testset "ad.jl" begin
    @testset "logp" begin
        # Hand-written log probabilities for vector `x = [s, m]`.
        function logp_gdemo_default(x)
            s = x[1]
            m = x[2]
            dist = Normal(m, sqrt(s))

            return logpdf(InverseGamma(2, 3), s) + logpdf(Normal(0, sqrt(s)), m) +
                logpdf(dist, 1.5) + logpdf(dist, 2.0)
        end

        test_model_ad(gdemo_default, logp_gdemo_default)

        @model function wishart_ad()
            v ~ Wishart(7, [1 0.5; 0.5 1])
        end

        # Hand-written log probabilities for `x = [v]`.
        function logp_wishart_ad(x)
            dist = Wishart(7, [1 0.5; 0.5 1])
            return logpdf(dist, reshape(x, 2, 2))
        end

        test_model_ad(wishart_ad(), logp_wishart_ad)
    end
end


text:
struct AnalyticalProblem{uType,tType,isinplace,P,F,C} <: AbstractAnalyticalProblem{uType,tType,isinplace}
  f::F
  u0::uType
  tspan::tType
  p::P
  callback::C
  @add_kwonly function AnalyticalProblem{iip}(f,u0,tspan,p=nothing;
           callback = nothing) where {iip}
    _tspan = promote_tspan(tspan)
    new{typeof(u0),typeof(_tspan),iip,typeof(p),
        typeof(f),typeof(callback)}(f,u0,_tspan,p,callback)
  end
end

function AnalyticalProblem(f,u0,tspan,p=nothing;kwargs...)
  iip = DiffEqBase.isinplace(f,4)
  AnalyticalProblem{iip}(f,u0,tspan,p;kwargs...)
end

export AnalyticalProblem, AbstractAnalyticalProblem


text:
using Test, DynamicPPL
dir = splitdir(splitdir(pathof(DynamicPPL))[1])[1]
include(dir*"/test/Turing/Turing.jl")
using .Turing

turnprogress(false)

@testset "DynamicPPL.jl" begin
    include("utils.jl")
    include("compiler.jl")
    include("varinfo.jl")
    include("prob_macro.jl")
    include("independence.jl")
end


text:
function read_hypergraph_label_names(dataname::String)
    names = String[]
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    # open("data/$dataname/label-names-$dataname.txt") do f
    open("$pathname/$dataname/label-names-$dataname.txt") do f
        for line in eachline(f)
            push!(names, line)
        end
    end
    return names
end

function read_hypergraph_labels(dataname::String)
    labels = Int64[]
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    # open("data/$dataname/node-labels-$dataname.txt") do f
    open("$pathname/$dataname/node-labels-$dataname.txt") do f
        for line in eachline(f)
            push!(labels, parse(Int64, line))
        end
    end
    return labels
end

function read_hypergraph_edges(dataname::String, maxsize::Int64=25, minsize::Int64=2)
    E = Dict{Integer, Dict}()
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    # open("data/$dataname/hyperedges-$dataname.txt") do f
    open("$pathname/$dataname/hyperedges-$dataname.txt") do f
        for line in eachline(f)
            edge = [parse(Int64, v) for v in split(line, ',')]
            sort!(edge)
            if minsize <= length(edge) <= maxsize
                sz = length(edge)
                if !haskey(E, sz)
                    E[sz] = Dict{}()
                end
                E[sz][edge] = 1
            end
        end
    end
    return E
end

function read_hypergraph_data(dataname::String, maxsize::Int64=25, minsize::Int64=2, return_labels=true)
    E = read_hypergraph_edges(dataname, maxsize, minsize)

    n = maximum([maximum(e) for k in keys(E) for e in keys(E[k])])
    D = zeros(Int64, n)
    for (sz, edges) in E
        for (e, _) in edges
            D[e] .+= 1
        end
    end

    maxedges = maximum(keys(E))
    for k in 1:maxedges
        if !haskey(E, k)
            E[k] = Dict{}()
        end
    end

    N = 1:n

    if return_labels
        labels = read_hypergraph_labels(dataname)
        return hypergraph(N, E, D), labels[N]
    end

    return hypergraph(N, E, D)
end

function hypermodularity_datasets()
    println("HyperModularity Package Datasets \n")
    pathname = joinpath(dirname(dirname(@__FILE__)),"data")
    foreach(readdir(pathname)) do f
       if isdir("$pathname/"*f)
           println("\t $f")
       end
    end
end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg