text:
##### temperature function
@inline function tempFunc(temp, p)
    k = exp(-p.Ea/(8.3145*(temp+273.15)))*(1.0-exp(temp+273.15 - p.T⁺))
    k = max(0.0, k)
    OGT_rate = exp(-p.Ea/(8.3145*(p.T⁺-2)))
    return k/OGT_rate
end

##### calculate photosynthesis rate (mmolC/individual/second)
@inline function calc_PS(par, temp, Chl, PRO, p)
    αI  = par * p.α * p.Φ
    PCm = p.PCmax * tempFunc(temp, p)
    PS  = PCm * (1.0 - exp(-αI / max(1.0e-30, PCm) * Chl / max(1.0e-30, PRO))) * PRO
    return PS
end

##### calculate nutrient uptake rate (mmolN/individual/second)
@inline function calc_NP_uptake(NH4, NO3, PO4, temp, NST, PST, PRO, DNA, RNA, Chl, p, ac)
    N_tot = total_N_biomass(PRO, DNA, RNA, NST, Chl, p)
    P_tot = total_P_biomass(DNA, RNA, PST, p)
    regQN = max(0.0, min(1.0, (p.NSTmax - NST / max(1.0e-30, N_tot)) / (p.NSTmax - p.NSTmin)))
    regQP = max(0.0, min(1.0, (p.PSTmax - PST / max(1.0e-30, P_tot)) / (p.PSTmax - p.PSTmin)))
    VNH4 = p.VNH4max * regQN * NH4/max(1.0e-30, NH4+p.KsatNH4) * tempFunc(temp, p) * PRO * ac
    VNO3 = p.VNO3max * regQN * NO3/max(1.0e-30, NO3+p.KsatNO3) * tempFunc(temp, p) * PRO * ac
    VPO4 = p.VPO4max * regQP * PO4/max(1.0e-30, PO4+p.KsatPO4) * tempFunc(temp, p) * PRO * ac
    return VNH4, VNO3, VPO4
end

@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.PRO[i], p) * plank.ac[i]

    @inbounds plank.VNH4[i], plank.VNO3[i], plank.VPO4[i] =
                            calc_NP_uptake(nuts.NH4[i], nuts.NO3[i], nuts.PO4[i], nuts.T[i],
                                        plank.NST[i], plank.PST[i], plank.PRO[i],
                                        plank.DNA[i], plank.RNA[i], plank.Chl[i], p, plank.ac[i])
end
function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P reserves for the first time of each time step
@kernel function update_quotas_1_kernel!(plank, ΔT)
    i = @index(Global)
    @inbounds plank.CH[i]  += plank.PS[i]   * ΔT
    @inbounds plank.PST[i] += plank.VPO4[i] * ΔT
    @inbounds plank.NST[i] +=(plank.VNH4[i] + plank.VNO3[i]) * ΔT
end
function update_quotas_1!(plank, ΔT, arch)
    kernel! = update_quotas_1_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT)
    wait(device(arch), event)
    return nothing
end

##### calculate DOC uptake rate (mmolC/individual/second)
##### DOC uptake needs support of photosynthesis for at least 5% of total C acquisition.
@inline function calc_DOC_uptake(DOC, temp, CH, PRO, DNA, RNA, Chl, p)
    C_tot = total_C_biomass(PRO, DNA, RNA, CH, Chl)
    regQ = max(0.0, min(1.0, (p.CHmax - CH / max(1.0e-30, C_tot)) / (p.CHmax - p.CHmin)))
    VN = p.VDOCmax * regQ * DOC/max(1.0e-30, DOC+p.KsatDOC) * tempFunc(temp, p) * PRO
    return VN
end
@kernel function calc_organic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.VDOC[i] = calc_DOC_uptake(nuts.DOC[i], nuts.T[i], plank.CH[i], plank.PRO[i], plank.DNA[i], plank.RNA[i], plank.Chl[i], p) * plank.ac[i]

    @inbounds plank.VDOC[i] = plank.VDOC[i] * isless(0.05, plank.PS[i]/(plank.VDOC[i]+plank.PS[i]))
end
function calc_organic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_organic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### calculate ρChl
@kernel function calc_ρChl_kernel!(plank, par, p)
    i = @index(Global)
    @inbounds plank.ρChl[i] = plank.PS[i]/max(1.0e-30, plank.PRO[i]) / max(1.0e-30, par[i] * p.α * p.Φ * plank.Chl[i]/max(1.0e-30, plank.PRO[i])) *
                             isless(1.0e-1, par[i]) * plank.ac[i]
end
function calc_ρChl!(plank, par, p, arch)
    kernel! = calc_ρChl_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, par, p)
    wait(device(arch), event)
    return nothing
end

##### calculate respiration (mmolC/individual/second)
@kernel function calc_respir_kernel!(plank, T, p)
    i = @index(Global)
    @inbounds plank.resp[i] = p.respir * plank.PRO[i] * tempFunc(T[i], p) * plank.ac[i]
end
function calc_respir!(plank, T, p, arch)
    kernel! = calc_respir_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, T, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P reserves for the second time of each time step
##### respiration first use Carbohydrate, if it's not enough, use protein then.
@kernel function update_quotas_2_kernel!(plank, ΔT, p)
    i = @index(Global)
    @inbounds plank.CH[i]  = plank.CH[i]  + (plank.VDOC[i] - plank.resp[i]) * ΔT
    @inbounds plank.CH[i]  = plank.CH[i]  + max(0.0, (0.0 - plank.CH[i]))
    @inbounds plank.PRO[i] = plank.PRO[i] - max(0.0, (0.0 - plank.CH[i]))
    @inbounds plank.NST[i] = plank.NST[i] + max(0.0, (0.0 - plank.CH[i])) * p.R_NC_PRO
end
function update_quotas_2!(plank, ΔT, p, arch)
    kernel! = update_quotas_2_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT, p)
    wait(device(arch), event)
    return nothing
end

##### calculate protein, DNA, RNA synthesis (mmol C /individual/second)
@kernel function calc_BS_kernel!(plank, T, p)
    i  = @index(Global)
    @inbounds limit_PRO = min(plank.CH[i]/(plank.CH[i] + p.k_sat_pro * p.Nsuper), plank.NST[i]/(plank.NST[i] + p.k_sat_pro * p.Nsuper * p.R_NC_PRO))
    @inbounds limit_DNA = min(plank.CH[i]/(plank.CH[i] + p.k_sat_dna * p.Nsuper), plank.NST[i]/(plank.NST[i] + p.k_sat_dna * p.Nsuper * p.R_NC_DNA),
                                plank.PST[i]/(plank.PST[i] + p.k_sat_dna * p.Nsuper * p.R_PC_DNA))
    @inbounds limit_RNA = min(plank.CH[i]/(plank.CH[i] + p.k_sat_rna * p.Nsuper), plank.NST[i]/(plank.NST[i] + p.k_sat_rna * p.Nsuper * p.R_NC_RNA),
                                plank.PST[i]/(plank.PST[i] + p.k_sat_rna * p.Nsuper * p.R_PC_RNA))

    @inbounds plank.S_PRO[i] = p.k_pro * plank.RNA[i] * limit_PRO * tempFunc(T[i], p)
    @inbounds plank.S_DNA[i] = p.k_dna * plank.PRO[i] * limit_DNA * tempFunc(T[i], p) * isless(plank.DNA[i]/(p.C_DNA * p.Nsuper), 2.0)
    @inbounds plank.S_RNA[i] = p.k_rna * plank.PRO[i] * limit_RNA * tempFunc(T[i], p)
end
function calc_BS!(plank, T, p, arch)
    kernel! = calc_BS_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, T, p)
    wait(device(arch), event)
    return nothing
end

##### update C, N, P reserves, protein, DNA, RNA, Chla
@kernel function update_biomass_kernel!(plank, p, ΔT)
    i = @index(Global)
    @inbounds S_Chl = plank.S_PRO[i] * plank.ρChl[i]
    @inbounds plank.PRO[i] += ΔT * plank.S_PRO[i]
    @inbounds plank.DNA[i] += ΔT * plank.S_DNA[i]
    @inbounds plank.RNA[i] += ΔT * plank.S_RNA[i]
    @inbounds plank.CH[i]  -= ΔT *(plank.S_PRO[i] + plank.S_DNA[i] + plank.S_RNA[i] + S_Chl)
    @inbounds plank.NST[i] -= ΔT *(plank.S_PRO[i] * p.R_NC_PRO + plank.S_DNA[i] * p.R_NC_DNA + plank.S_RNA[i] * p.R_NC_RNA + S_Chl*4.0/55.0)
    @inbounds plank.PST[i] -= ΔT *(plank.S_DNA[i] * p.R_PC_DNA + plank.S_RNA[i] * p.R_PC_RNA)
    @inbounds plank.Chl[i] += ΔT * S_Chl * 893.49 / 55.0 # chl unit is mgChl/cell
    @inbounds plank.age[i] += ΔT / 3600.0 * plank.ac[i]
end
function update_biomass!(plank, p, ΔT, arch)
    kernel! = update_biomass_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p, ΔT)
    wait(device(arch), event)
    return nothing
end

##### calculate exudation of carbon. Nitrogen and phosphorus will not be exuded for now
@kernel function calc_exudation_kernel!(plank, p)
    i = @index(Global)
    @inbounds tot_C = total_C_biomass(plank.PRO[i], plank.DNA[i], plank.RNA[i], plank.CH[i], plank.Chl[i])
    @inbounds plank.exu[i] = max(0.0, plank.CH[i] - p.CHmax * tot_C)
end
function calc_exudation!(plank, p, arch)
    kernel! = calc_exudation_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end

@kernel function update_CH_kernel!(plank)
    i = @index(Global)
    @inbounds plank.CH[i] -= plank.exu[i]
end
function update_CH!(plank, arch)
    kernel! = update_CH_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank)
    wait(device(arch), event)
    return nothing
end


text:
module Jlrs
struct LocalTask
    func::Function
    args::Tuple
    kwargs
    wakeptr::Ptr{Cvoid}
end

const color = Ref{Bool}(false)
const wakerust = Ref{Ptr{Cvoid}}(C_NULL)

function asynccall(func::Function, wakeptr::Ptr{Cvoid}, args...; kwargs...)::Task
    @nospecialize func wakeptr args kwargs
    Base.Threads.@spawn begin
        try
            func(args...; kwargs...)
        finally
            if wakeptr != C_NULL
                ccall(wakerust[], Cvoid, (Ptr{Cvoid},), wakeptr)
            end
        end
    end
end

function asynccall(func::Function, args...; kwargs...)::Task
    @nospecialize func args kwargs
    Base.Threads.@spawn func(args...; kwargs...)
end

const inchannel = Channel{LocalTask}(1)
const outchannel = Channel{Task}(1)
Base.Threads.@spawn begin
    while true
        local_task::LocalTask = take!(inchannel)
        task::Task = @async begin
            try
                local_task.func(local_task.args...; local_task.kwargs...)
            finally
                if local_task.wakeptr != C_NULL
                    ccall(wakerust[], Cvoid, (Ptr{Cvoid},), local_task.wakeptr)
                end
            end
        end
        put!(outchannel, task)
    end
end

function scheduleasync(func::Function, wakeptr::Ptr{Cvoid}, args...; kwargs...)::Task
    @nospecialize func wakeptr args kwargs
    task::LocalTask = LocalTask(func, args, kwargs, wakeptr)
    put!(inchannel, task)
    take!(outchannel)
end

function scheduleasync(func::Function, args...; kwargs...)::Task
    @nospecialize func args kwargs
    task::LocalTask = LocalTask(func, args, kwargs, C_NULL)
    put!(inchannel, task)
    take!(outchannel)
end

function valuestring(@nospecialize(value::Any))::String
    io = IOBuffer()
    show(io, "text/plain", value)
    String(take!(io))
end

function errorstring(@nospecialize(value::Any))::String
    io = IOBuffer()
    showerror(IOContext(io, :color => color[]), value)
    String(take!(io))
end
end

ext: jl

max_stars_repo_name: marhkb/jlrs

max_issues_repo_name: marhkb/jlrs

max_forks_repo_name: marhkb/jlrs


text:
# Check if file is a data file
function is_solution_restart_file(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # If attribute "mesh_file" exists, this must be a data file
    return haskey(attributes(file), "mesh_file")
  end
end

# Use data file to extract mesh filename from attributes
function extract_mesh_filename(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract filename relative to data file
    mesh_file = read(attributes(file)["mesh_file"])

    return joinpath(dirname(filename), mesh_file)
  end
end

function extract_mesh_information(mesh::Trixi.TreeMesh)
  center_level_0 = mesh.tree.center_level_0
  length_level_0 = mesh.tree.length_level_0

  coordinates = mesh.tree.coordinates[:, Trixi.leaf_cells(mesh.tree)]
  levels = mesh.tree.levels[Trixi.leaf_cells(mesh.tree)]

  return coordinates, levels, center_level_0, length_level_0
end

# Read in data file and return all relevant information
function read_datafile(filename::String)
  # Open file for reading
  h5open(filename, "r") do file
    # Extract basic information
    ndims_ = read(attributes(file)["ndims"])
    polydeg = read(attributes(file)["polydeg"])
    n_elements = read(attributes(file)["n_elements"])
    n_variables = read(attributes(file)["n_vars"])
    time = read(attributes(file)["time"])

    # Extract labels for legend
    labels = Array{String}(undef, 1, n_variables)
    for v = 1:n_variables
      labels[1, v] = read(attributes(file["variables_$v"])["name"])
    end

    # Extract data arrays
    n_nodes = polydeg + 1

    if ndims_ == 2
      data = Array{Float64}(undef, n_nodes, n_nodes, n_elements, n_variables)
      for v = 1:n_variables
        vardata = read(file["variables_$v"])
        @views data[:, :, :, v][:] .= vardata
      end
    elseif ndims_ == 3
      data = Array{Float64}(undef, n_nodes, n_nodes, n_nodes, n_elements, n_variables)
      for v = 1:n_variables
        vardata = read(file["variables_$v"])
        @views data[:, :, :, :, v][:] .= vardata
      end
    else
      error("Unsupported number of spatial dimensions: ", ndims_)
    end

    # Extract element variable arrays
    element_variables = Dict{String, Union{Vector{Float64}, Vector{Int}}}()
    index = 1
    while haskey(file, "element_variables_$index")
      varname = read(attributes(file["element_variables_$index"])["name"])
      element_variables[varname] = read(file["element_variables_$index"])
      index +=1
    end

    return labels, data, n_elements, n_nodes, element_variables, time
  end
end


text:
using PolyChaos, QuadGK

N = 10
gauss = OrthoPoly("gaussian",N)
pipj(t,i,j) = gauss.w(t)*evaluate(i,t,gauss)*evaluate(j,t,gauss)
[ quadgk(t->pipj(t,i,i),-Inf,Inf) - factorial(i) for i=1:N ]

op = OrthoPoly("uniform01",N)
pipj(t,i,j) = op.w(t)*evaluate(i,t,op)*evaluate(j,t,op)
for i = 0:N, j=0:N
    display(quadqk(t->op.w(t)*evaluate(i,t,op)*evaluate(j,t,op),0,1)[1])
end


text:
### Initial Value Problem: INPUT is state_matrix, p_drive, p_dest, C, T and OUTPUT initial_state

function solveinitialvalueproblem(state_matrix, transition_matrix, p_drive, p_dest, C, number_zones)

    initial_state = zeros(C,1);

    for t=1:(T-1)

        # Sampling if a car drives
        for i=1:C
            origin = state_matrix[i,t];
            RndVar = rand();
            driving_probability = p_drive[origin,t];
            if(RndVar <= driving_probability)
                drive = 1;
            else
                drive = 0;
                transition_matrix[i,t,2]= origin;
            end
            transition_matrix[i,t,1] = drive;
        end

        # Sampling travel destinations
        for i=1:C
            drive = transition_matrix[i,t,1];
            if (drive == 1)
                RndVar = rand();
                range_up = 0;
                range_low = 0;
                destination = 0;
                origin = state_matrix[i,t];
                distribution = p_dest[origin, :, t];
                if(sum(distribution) == 0)
                    destination = origin;
                else
                    for j=1:number_zones
                        range_up = range_up + distribution[j];
                        if(range_low < RndVar <= range_up)
                            destination = j;
                            break;
                        end
                        range_low = range_up;
                    end
                end
                transition_matrix[i,t,2] = destination;
            end
        end

        # state update
        state_matrix[:,t+1] = transition_matrix[:,t,2];
        println(trunc(Int,t/0.23),"% of the Initial Value Problem is solved")
    end

    initial_state = state_matrix[:,24];
    initial_state = round.(Int, initial_state)

    return initial_state

end

max_stars_repo_name: ArsamAryandoust/CarParkingMaps

max_issues_repo_name: ArsamAryandoust/CarParkingMaps

max_forks_repo_name: ArsamAryandoust/CarParkingMaps


text:
function get_U0(γ, β, rdisp, Nₚ)
    psamples_raw = get_psamples(rdisp, Nₚ)
    s = 0.0
    for psample in psamples_raw
        εₚ⁺, εₚ⁻, wₚ = psample
        κₚ = sqrt(εₚ⁻^2 + γ^2)
        s += wₚ*(tanh(β*0.5*(εₚ⁺+κₚ)) + tanh(β*0.5*(κₚ-εₚ⁺)))/κₚ
    end
    return 4.0/s
end

check_gamma(γ::Real, params::AbstractParams, rdisp::ReducedDispersion, Nₚ) = check_gamma(γ, params, get_samples(rdisp, Nₚ))
check_gamma(γ::Real, params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}}) = check_gamma(γ, params.β, params.u, psamples_raw)

function check_gamma(γ::Real, β::Real, u::Real, psamples_raw::Vector{NTuple{3, Float64}})
    s = 0.0
    for psample in psamples_raw
        εₚ⁺, εₚ⁻, wₚ = psample
        κₚ = sqrt(εₚ⁻^2 + γ^2)
        s += wₚ*(tanh(β*0.5*(εₚ⁺+κₚ)) + tanh(β*0.5*(κₚ-εₚ⁺)))/κₚ
    end
    return 0.25*s-u
end

function get_gamma(params::AbstractParams, rdisp::ReducedDispersion, Nₚ, range)
    psamples_raw = get_psamples(rdisp, Nₚ)
    return find_zero((γ)->check_gamma(γ, params, psamples_raw), range, Bisection())
end

function get_gamma(params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}}, range)
    return find_zero((γ)->check_gamma(γ, params.β, params.u, psamples_raw), range, Bisection())
end

function get_gamma(β::Real, u::Real, psamples_raw::Vector{NTuple{3, Float64}}, range)
    return find_zero((γ)->check_gamma(γ, β, u, psamples_raw), range, Bisection())
end

calculate_sce_rhs(k, params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64) = calculate_sce_rhs(k, params, get_psamples(rdisp, Nₚ))

function calculate_sce_rhs(k::T, params::AbstractParams, psamples_raw::AbstractVector{NTuple{3, Float64}}) where {T<:Number}
    β, m, U0 = params.β, params.m, 1/params.u
    k² = k^2
    γ = 4*K(k²)*m/β
    kk = 4k/(1+k)^2
    s = zero(T)
    @fastmath for psample in psamples_raw
        εₚ⁺, εₚ⁻, wₚ = psample
        nₚ = γ^2*k/(εₚ⁻^2+0.25*γ^2*(1+k)^2)
		κₚ = sqrt((εₚ⁻^2+0.25*γ^2*(1-k)^2)/(1.0+0.25*(γ/εₚ⁻)^2*(1+k)^2))*Pi(nₚ, pi/2, kk)/K(kk)
        s += wₚ*(tanh(β*0.5*(εₚ⁺+κₚ)) + tanh(β*0.5*(κₚ-εₚ⁺)))*abs(εₚ⁻)/sqrt((εₚ⁻^2+0.25*γ^2*(1+k)^2)*(εₚ⁻^2+0.25*γ^2*(1-k)^2))
    end
    return 1.0 - 0.25*U0*s
end

function solve_sce_nob1(params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}}, range = (0.5, 0.99999))
    k = find_zero((x)->calculate_sce_rhs(x, params, psamples_raw), range, Bisection())
    γ = 4*K(k^2)*params.m/params.β
    return k, γ
end

function _seed_sn(params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}}, range = (0.5, 0.99999))
    β, m, N, N′ = params.β, params.m, params.N, get_length(params)
    k, γ = solve_sce_nob1(params, psamples_raw, range)
    dτ = β/(N*m)
    τs = [dτ*(i-1) for i = 1:N]
    sns = zeros(N′)
    sns[1:N] .= k*γ*Jacobi.sn.(γ*τs, k^2)
    pin_bs!(sns, params)
    return sns
end

function _seed_sn_asymp(params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}}, range = (0.001, 1000.0))
    β, m, N, N′ = params.β, params.m, params.N, get_length(params)
    γ = get_gamma(params, psamples_raw, range)
    dτ = β/(N*m)
    τs = [dτ*(i-1) for i = 1:N]
    sns = zeros(N′)
    N1 = div(N, 4)
    N2 = div(N, 2)
    sns[1:N1] .= γ.*tanh.(γ.*τs[1:N1])
    sns[(N1+1):N2] .= γ.*tanh.(γ.*(τs[N2+1].-τs[(N1+1):N2]))
    sns[(N2+1):N] .= -sns[1:N2]
    pin_bs!(sns, params)
    return sns
end

function get_static_fmin(β::Real, u::Real, psamples_raw::Vector{NTuple{3, Float64}}, range = (0.001, 1000.0))
	γ = get_gamma(β, u, psamples_raw, range)
	γ′ = Float128(γ)
	β′ = Float128(β)
	s = β*u*γ^2
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
		εₚ⁻′ = Float128(εₚ⁻)
		εₚ⁺′ = Float128(εₚ⁺)
		κₚ′ = sqrt(γ′^2 + εₚ⁻′^2)
		s -= wₚ*Float64(log((cosh(β′*κₚ′) + cosh(β′*εₚ⁺′))/(cosh(β′*εₚ⁻′) + cosh(β′*εₚ⁺′))))
	end
	return γ, s
end

seed_sn(params::AbstractParams, rdisp::ReducedDispersion, Nₚ::Int64) = seed_sn(params, get_psamples(rdisp, Nₚ))
function seed_sn(params::AbstractParams, psamples_raw::Vector{NTuple{3, Float64}})
    try
        _seed_sn(params, psamples_raw)
    catch e
        _seed_sn_asymp(params, psamples_raw)
    end
end

function get_τs(params::AbstractParams)
    β, m, N = params.β, params.m, params.N
    dτ = β/(N*m)
    return [dτ*(i-1) for i = 1:N]
end
function get_τs(β::Real, m::Int64, N::Int64)
    dτ = β/(N*m)
    return [dτ*(i-1) for i = 1:N]
end

show_sn(params::AbstractParams, rdisp::ReducedDispersion, Nₚ, range = (0.5, 0.99999)) = get_τs(params), seed_sn(params, rdisp, Nₚ, range)

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

seed_const(c, params::TTCrystal.AbstractParamsNoB1) = fill!(Vector{Float64}(undef, params.N),c)
function seed_const(c1, c2, params::AbstractParamsB1)
    N = params.N
    bs = Vector{Float64}(undef, 2N)
    bs[1:N] .= c1
    bs[(N+1):2N] .= c2
    return bs
end

function seed_const(c1, c2, c3, params::Union{ParamsB1B3, ParamsB3})
    N = params.N
	bs = Vector{Float64}(undef, get_length(params))
    bs[1:N] .= c1
    bs[(N+1):2N] .= c2
	bs[(2N+1):3N] .= c3
    return bs
end

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

isB1(params::AbstractParamsB1) = Val(true)
isB1(params::AbstractParamsNoB1) = Val(false)

function show_bs(res::Optim.MultivariateOptimizationResults, ::Val{true})
    bs = Optim.minimizer(res)
    N = div(length(bs), 2)
    b0 = bs[1:N]
    b1 = bs[(N+1):2N]
    plot(b0)
    plot!(b1)
end
function show_bs(res::Optim.MultivariateOptimizationResults, ::Val{false})
    bs = Optim.minimizer(res)
    plot(bs)
end

show_bs(res, params::AbstractParams) = show_bs(res, isB1(params))

function check_derivative(bs)
    N = div(length(bs), 2)
    b0 = bs[1:N]
    b1 = bs[(N+1):2N]
    b_der = similar(b0)
    b_der[1] = b0[2]-b0[end]
    b_der[end] = b0[1]-b0[end-1]
    for i=2:(N-1)
        b_der[i] = b0[i+1] - b0[i-1]
    end
    #b_der .*= -sum(abs2, b1)/sum(abs2, b_der)
	b_der .*= b1[div(N,2)+1]/b_der[div(N,2)+1]
    plot(b_der)
    plot!(b1)
end
check_derivative(res::Optim.MultivariateOptimizationResults) = check_derivative(Optim.minimizer(res))

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

function el_part(γ, P, β, psamples_raw)
	β′ = Float128(β)
    γ′ = Float128(γ)
    s = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
		εₚ⁻′ = Float128(εₚ⁻)
		εₚ⁺′ = Float128(εₚ⁺)
        κₚ′ = sqrt(γ′^2 + (εₚ⁻′ - P)^2)
		s += wₚ*Float64(log((cosh(β′*κₚ′) + cosh(β′*εₚ⁺′))/(cosh(β′*εₚ⁻′) + cosh(β′*εₚ⁺′))))
	end
    return s
end

function static_f_energy(γ, P, β, u, u₃, psamples_raw)
    β*(u*γ^2 + u₃*P^2) - el_part(γ, P, β, psamples_raw)
end

function energy_grad!(g, v, β, u, u₃, psamples_raw)
    γ, P = v[1], v[2]
    g[1] = 2*β*u*γ
    g[2] = 2*β*u₃*P
    s = 0.0
    s1 = 0.0
	for (εₚ⁺, εₚ⁻, wₚ) in psamples_raw
        κₚ = sqrt(γ^2 + (εₚ⁻-P)^2)
        s += wₚ*(tanh(0.5*β*(κₚ+εₚ⁺)) + tanh(0.5*β*(κₚ-εₚ⁺)))/κₚ
        s1 += wₚ*(εₚ⁻-P)*(tanh(0.5*β*(κₚ+εₚ⁺)) + tanh(0.5*β*(κₚ-εₚ⁺)))/κₚ
    end
    g[1] -= β*0.5*γ*s
    g[2] += β*0.5*s1
end

function optimize_γ_P(γ₀, P₀, β, u, u₃, psamples_raw)
    f(v) = static_f_energy(v[1], v[2], β, u, u₃, psamples_raw)
    g!(g, v) = energy_grad!(g, v, β, u, u₃, psamples_raw)
    x0 = [γ₀, P₀]
    return optimize(f, g!, x0, Optim.Options(g_tol = 1e-10, show_trace = true))
end


text:
using Base.Order

abstract AbstractKDTree <: AbstractNearestNeighborTree

type EmptyKDTree <: AbstractKDTree
end

type KDTreeNode{T <: Real} <: AbstractKDTree
    k::     Vector{T}       # multidimensional point
    i::		Int             # index
    d::     Int             # split dimension
    s::     Float64         # split value
    leaf::  Bool            # is leaf?
    left::  AbstractKDTree  # left node
    right:: AbstractKDTree  # right node

    KDTreeNode(k::Vector{T}, i::Int) = new(k, i, 1, NaN, true, EmptyKDTree(), EmptyKDTree())
    KDTreeNode(k::Vector{T}, i::Int, d) = new(k, i, d, NaN, true, EmptyKDTree(), EmptyKDTree())
end

type KDTree <: NearestNeighborTree
	metric::Metric
	root:: AbstractKDTree

	KDTree() = new (Euclidean(), EmptyKDTree())
	KDTree(metric::Metric) = new(metric, EmptyKDTree())
	KDTree(k, i, metric::Metric) = new(metric, setindex!(EmptyKDTree(), i, k))

	function KDTree{T <: Real}(X::Matrix{T}, metric::Metric = Euclidean())
		n = size(X, 2)
		t = KDTree(X[:,1], 1, metric)
		for i = 2 : n
			setindex!(t.root, i, X[:,i])
		end
		new(t.metric, t.root)
	end
end

function Base.show(io::IO, t::KDTree)
	println(io, "KDTree: ", t.metric)
end

Base.setindex!{T <: Real}(t::EmptyKDTree, i::Int, k::Vector{T}) = KDTreeNode{T}(k, i)
Base.setindex!(t::KDTree, i::Int, k) = (t.root = setindex!(t.root, i, k); t)

function Base.setindex!{T <: Real}(n::KDTreeNode{T}, i::Int, k::Vector{T})
	if n.leaf
		left_k = k
		left_i = i
		right_k = n.k
		right_i = n.i
		if right_k[n.d] < left_k[n.d]
			left_k = n.k
			left_i = n.i
			right_k = k
			right_i = i
		end
		n.s = 0.5*(right_k[n.d] + left_k[n.d])
		next_sd = ((n.d) % length(k))+1

		n.left = KDTreeNode{T}(left_k, left_i, next_sd)
		n.right = KDTreeNode{T}(right_k, right_i, next_sd)
		n.leaf = false
	else
		setindex!(k[n.d] <= n.s ? n.left : n.right, i, k)
	end
	return n
end

disp(t::KDTree) = disp(t.root)
function disp(n::KDTreeNode, l::Int = 0)
	print(repeat("\t", l), "L: $(l), SD: $(n.d), SV: $(n.s)")

	if typeof(n.left) <: EmptyKDTree && typeof(n.right) <: EmptyKDTree
		println(" = K: $(n.k) ($(n.i))")
	else
		println("")
		disp(n.right, l+1)
		disp(n.left, l+1)
	end
end

function search_leaf{T <: Real}(n::KDTreeNode{T},
						x::Vector{T},
						h::Array{Float64,1},
						index::Array{Int},
						m::Metric)

	if n.leaf
		d = evaluate(m, x, n.k)
		if (d < h[1])
			h[end]=d
			index[end]=n.i
			percolate_down!(h, index, 1, d, n.i, Base.Order.Reverse)
		end
	else
		# Determine nearest and furtherest branch
		if (x[n.d] > n.s)
			near = n.right
			far = n.left
		else
			near = n.left
			far = n.right
		end

		# Search the nearest branch
		search_leaf(near, x, h, index, m)

		# Only search far tree if do not have enough neighbors
		d = evaluate(m, float64(x[n.d]), n.s)
		if d <= h[1]
			search_leaf(far, x, h, index, m)
		end
	end
end

function nearest{T <: Real}(t::KDTree,
							x::Vector{T},
							k::Int)
	h = fill(Inf, k+1)
	index = fill(0, k+1)

	search_leaf(t.root, x, h, index, t.metric)

	return index[1:end-1], h[1:end-1]
end

function inball{T <: Real}(t::KDTree,
						   v::Vector{T},
	                       r::Real,
	                       exclude::Integer = -1)
	error("Not yet implemented")
end

# Binary heap indexing
heapleft(i::Integer) = 2i
heapright(i::Integer) = 2i + 1

# Binary min-heap percolate down.
function percolate_down!(xs::AbstractArray, xis::AbstractArray, i::Integer, x=xs[i], xi=xis[i], o::Ordering=Forward, len::Integer=length(xs))
    @inbounds while (l = heapleft(i)) <= len
        r = heapright(i)
        j = r > len || lt(o, xs[l], xs[r]) ? l : r
        if lt(o, xs[j], x)
            xs[i] = xs[j]
            xis[i] = xis[j]
            i = j
        else
            break
        end
    end
    xs[i] = x
    xis[i] = xi
end


text:
using CSV, DataFrames
using Plots; Plots.plotly()
using StatsBase

realdata = CSV.read("korea.csv", DataFrame)
구분 = unique(realdata.gubun)

temp = replace.(realdata.stdDay, "2021년" => "21년")
temp = replace.(temp, "2020년" => "20년")
for i in 1:9
    temp = replace.(temp, "년 $(i)월" => "년 0$(i)월")
    temp = replace.(temp, "월 $(i)일" => "월 0$(i)일")
end
for i in lpad.(0:24, 2, "0")
    temp = replace.(temp, "일 $(i)시" => "")
end
temp = replace.(temp, "년 " => "")
temp = replace.(temp, "월 " => "")
realdata.stdDay = temp
sort!(realdata, :stdDay)

realdata

Plots.plotly()

time_evolution = plot()
for 지역 in 구분
    plot!(time_evolution, realdata[realdata.gubun .== 지역,:defCnt], label = 지역)
end
time_evolution

time_evolution = plot()
for 지역 in 구분
    if 지역 == "합계" continue end
    plot!(time_evolution, realdata[realdata.gubun .== 지역,:defCnt], label = 지역)
end
time_evolution

Plots.gr()

누적확진자 = DataFrame()
누적확진자["stdDay"] = unique(realdata.stdDay)
for 지역 in 구분
    누적확진자 = rightjoin(누적확진자, unique(realdata[realdata.gubun .== 지역,[1,3]],1), on = :stdDay)
    rename!(누적확진자, :defCnt => 지역)
end

일일확진자 = DataFrame(Matrix(누적확진자[2:end,2:end]) - Matrix(누적확진자[1:(end-1),2:end]))
일일확진자 = coalesce.(일일확진자,0)
일일확진자 = max.(일일확진자, 0)
rename!(일일확진자, 구분)
일일확진자 = hcat(DataFrame(stdDay = unique(누적확진자.stdDay)[2:end]),일일확진자)
일일확진자["합계"] = vec(sum(Matrix(일일확진자[2:(end-1)]), dims = 2))

time_evolution = plot()
for 지역 in 구분
    if 지역 == "합계" continue end
    plot!(time_evolution, 일일확진자[:,지역], label = 지역)
end
time_evolution

savefig(time_evolution,"time_evolution.html")

entropy_ = zeros(size(일일확진자)[1])
for i in 1:length(entropy_)
    entropy_[i] = entropy(collect(일일확진자[i,2:(end-1)]) ./ 일일확진자.합계[i],18)
end
plot(entropy_, ylims = (-0.1,1.1))
savefig("entropy.html")

CSV.write("일일확진자.csv", 일일확진자)
CSV.write("누적확진자.csv", 누적확진자)

max_stars_repo_name: rmsmsgood/aespa

max_issues_repo_name: rmsmsgood/aespa

max_forks_repo_name: rmsmsgood/aespa


text:
# Generate 3 Gaussian peak at position : 5, 10, 20 with difference
# shape factors σ 1, 1.5, 2
#
using DelimitedFiles,Random
using NLS_Fit

Random.seed!(1234)

model = Gaussian_Peak() + Gaussian_Peak() + Gaussian_Peak()
θ1 = Float64[1,5,1]
θ2 = Float64[1,10,1.5]
θ3 = Float64[1,20,2]
θ = vcat(θ1,θ2,θ3)

X=Float64[1:0.25:30;]
n=length(X)
Y=eval_y(model,X,θ) + 0.1*(rand(n) .- 0.5)

rootDir = joinpath(dirname(pathof(NLS_Fit)), "..")
dataDir = joinpath(rootDir,"data")
dataFile = joinpath(dataDir,"mapped_parameters.txt")
writedlm(dataFile,hcat(X,Y))


text:
using Test
using MiraiBots, MiraiBots.Messages, MiraiBots.Events, MiraiBots.MessageChains, MiraiBots.Commands
using MiraiBots: send, Commands
using Dates
using JSON3

(; qq, key, server, image_url) = open(JSON3.read, "config.json")

# ENV["JULIA_DEBUG"] = MiraiBots

@testset "No Broadcaster" begin
    using MiraiBots, MiraiBots.Messages, MiraiBots.Events, MiraiBots.MessageChains
    using MiraiBots: Commands, send
    bot = MiraiBots.WebSocketAdapter()
    @sync begin
        @async MiraiBots.loop(bot, server, qq, key)
        # 需等待首条消息才算连接到HTTP API
        MiraiBots.receive_or_throw(bot)

        @async for msg in MiraiBots.get_output_channel(bot)
            if msg isa MiraiBots.ExceptionAndBacktrace
                @error "Error occurred in adapter" exception = (msg.exception, msg.backtrace)
            elseif msg isa FriendMessage
                chain = msg.messageChain
                send(bot, Commands.sendFriendMessage(
                    target = msg.sender.id, quoteId = chain[1].id,
                    messageChain = chain[2:end]))
            elseif msg isa NewFriendRequestEvent
                send(bot, Commands.resp_newFriendRequestEvent(
                    eventId = msg.eventId, fromId = msg.fromId, groupId = msg.groupId,
                    operate = Commands.NewFriendOperations.REFUSE, message = "sorry~"))
            else
                @show msg
            end
        end

        sleep(10)
        close(bot)
    end
end

@testset "Broadcaster" begin
    broadcaster = MiraiBots.Broadcaster(MiraiBots.HTTPAdapter())
    # 匿名函数
    MiraiBots.register(broadcaster) do bot, msg::FriendMessage
        chain = msg.messageChain
        send(bot, Commands.sendFriendMessage(target = msg.sender.id, quoteId = chain[1].id, messageChain = chain[2:end]))
        throw(MiraiBots.ShutdownBroadcaster())
    end
    # 命名函数
    register(f) = MiraiBots.register(f, broadcaster)
    function logger(_, msg)
        @info "$(now()): $msg"
    end |> register

    function respond_image(bot, msg::NudgeEvent)
        if msg.subject.kind == Commands.NudgeKinds.FRIEND
            @info 1
            send(bot, Commands.sendFriendMessage(
                target = msg.subject.id,
                messageChain = [Image(url = image_url)]
            ))
        end
    end |> register
    # 启动
    MiraiBots.launch(broadcaster, server, qq, key)
end


text:
# # Maximum-minimum driving force analysis

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# In this notebook we use the max-min driving force analysis (MMDFA) to find
# optimal concentrations for the metabolites in glycolysis to ensure that the
# smallest driving force across all the reactions in the model is as large as
# possible. For more information, see Flamholz, Avi, et al.  "Glycolytic
# strategy as a tradeoff between energy yield and protein cost.", Proceedings
# of the National Academy of Sciences 110.24 (2013): 10039-10044.

using COBREXA, GLPK, Tulip

# Let's load the core E. coli model

model = load_model("e_coli_core.json")

# We need some thermodynamic data. You can get reaction Gibbs free energies (ΔG⁰) from
# e.g. [eQuilibrator](https://equilibrator.weizmann.ac.il/), possibly using the
# [Julia wrapper](https://github.com/stelmo/eQuilibrator.jl) that allows you to
# automate this step. Here, we make a dictionary that maps the reaction IDs to
# calculated Gibbs free energy of reaction for each reaction (including the transporters).

reaction_standard_gibbs_free_energies = Dict( # kJ/mol
    "ACALD" => -21.26,
    "PTAr" => 8.65,
    "ALCD2x" => 17.47,
    "PDH" => -34.24,
    "PYK" => -24.48,
    "CO2t" => 0.00,
    "MALt2_2" => -6.83,
    "CS" => -39.33,
    "PGM" => -4.47,
    "TKT1" => -1.49,
    "ACONTa" => 8.46,
    "GLNS" => -15.77,
    "ICL" => 9.53,
    "FBA" => 23.37,
    "SUCCt3" => -43.97,
    "FORt2" => -3.42,
    "G6PDH2r" => -7.39,
    "AKGDH" => -28.23,
    "TKT2" => -10.31,
    "FRD7" => 73.61,
    "SUCOAS" => -1.15,
    "FBP" => -11.60,
    "ICDHyr" => 5.39,
    "AKGt2r" => 10.08,
    "GLUSy" => -47.21,
    "TPI" => 5.62,
    "FORt" => 13.50,
    "ACONTb" => -1.62,
    "GLNabc" => -30.19,
    "RPE" => -3.38,
    "ACKr" => 14.02,
    "THD2" => -33.84,
    "PFL" => -19.81,
    "RPI" => 4.47,
    "D_LACt2" => -3.42,
    "TALA" => -0.94,
    "PPCK" => 10.65,
    "ACt2r" => -3.41,
    "NH4t" => -13.60,
    "PGL" => -25.94,
    "NADTRHD" => -0.01,
    "PGK" => 19.57,
    "LDH_D" => 20.04,
    "ME1" => 12.08,
    "PIt2r" => 10.41,
    "ATPS4r" => -37.57,
    "PYRt2" => -3.42,
    "GLCpts" => -45.42,
    "GLUDy" => 32.83,
    "CYTBD" => -59.70,
    "FUMt2_2" => -6.84,
    "FRUpts2" => -42.67,
    "GAPD" => 0.53,
    "H2Ot" => 0.00,
    "PPC" => -40.81,
    "NADH16" => -80.37,
    "PFK" => -18.54,
    "MDH" => 25.91,
    "PGI" => 2.63,
    "O2t" => 0.00,
    "ME2" => 12.09,
    "GND" => 10.31,
    "SUCCt2_2" => -6.82,
    "GLUN" => -14.38,
    "ETOHt2r" => -16.93,
    "ADK1" => 0.38,
    "ACALDt" => 0.00,
    "SUCDi" => -73.61,
    "ENO" => -3.81,
    "MALS" => -39.22,
    "GLUt2r" => -3.49,
    "PPS" => -6.05,
    "FUM" => -3.42,
);

# In general we cannot be certain that all fluxes will be positive for a given flux
# solution. This poses problems for systematically enforcing that ΔᵣG ≤ 0 for each reaction,
# because it implicitly assumes that all fluxes are positive, as done in the original
# formulation of MMDF. In `max_min_driving_force` we instead enforce ΔᵣG ⋅ vᵢ ≤ 0, where vᵢ
# is the flux of reaction i. By default all fluxes are assumed to be positive, but by
# supplying thermodynamically consistent flux solution it is possible to drop this implicit
# assumption and makes it easier to directly incorporate the max min driving force into
# non-customized models. Here, customized model means a model written such that a negative
# ΔᵣG is associated with each positive flux in the model, and only positive fluxes are used
# by the model.

flux_solution = flux_balance_analysis_dict( # find a thermodynamically consistent solution
    model,
    GLPK.Optimizer;
    modifications = [add_loopless_constraints()],
)

# Run max min driving force analysis with some reasonable constraints on metabolite
# concentration bounds. To remove protons and water from the concentration calculations, we
# explicitly specify their IDs. Note, protons and water need to be removed from the
# concentration calculation of the optimization problem, because the Gibbs free energies of
# biochemical reactions are measured at constant pH, so proton concentration is fixed, and
# reactions occur in aqueous environments, hence water concentration does not change.

sol = max_min_driving_force(
    model,
    reaction_standard_gibbs_free_energies,
    Tulip.Optimizer;
    flux_solution = flux_solution,
    proton_ids = ["h_c", "h_e"],
    water_ids = ["h2o_c", "h2o_e"],
    concentration_ratios = Dict(
        ("atp_c", "adp_c") => 10.0,
        ("nadh_c", "nad_c") => 0.13,
        ("nadph_c", "nadp_c") => 1.3,
    ),
    concentration_lb = 1e-6, # M
    concentration_ub = 100e-3, # M
    ignore_reaction_ids = [
        "H2Ot", # ignore water transporter
    ],
)

sol.mmdf

#md # !!! note "Note: transporters"
#md #       Transporters can be included in MMDF analysis, however water and proton
#md #       transporters must be excluded explicitly in `ignore_reaction_ids`. Due to
#md #       the way the method is implemented, the ΔᵣG for these transport reactions
#md #       will always be 0. If not excluded, the MMDF will only have a zero solution (if
#md #       these reactions are used in the flux solution).

# Next, we plot the results to show how the concentrations can be used to ensure that
# each reach proceeds "down hill" (ΔᵣG < 0) and that the driving force is as
# large as possible across all the reactions in the model. Compare this to the
# driving forces at standard conditions. Note, we only plot glycolysis for simplicity.

# We additionally scale the fluxes according to their stoichiometry in the
# pathway. From the output, we can clearly see that that metabolite concentrations
# play a large role in ensuring the thermodynamic consistency of in vivo reactions.

rids = ["GLCpts", "PGI", "PFK", "FBA", "TPI", "GAPD", "PGK", "PGM", "ENO", "PYK"] # glycolysis
rid_rf = [flux_solution[rid] for rid in rids]
dg_standard = cumsum([
    reaction_standard_gibbs_free_energies[rid] * flux_solution[rid] for rid in rids
])
dg_standard .-= first(dg_standard)
dg_opt = cumsum([sol.dg_reactions[rid] * flux_solution[rid] for rid in rids])
dg_opt .-= first(dg_opt)

using CairoMakie

fig = Figure();
ax = Axis(
    fig[1, 1],
    xticklabelrotation = -pi / 2,
    xlabel = "Reaction",
    ylabel = "Cumulative ΔG [kJ/mol]",
);

lines!(ax, 1:length(rids), dg_standard; color = :red, label = "Standard")
lines!(ax, 1:length(rids), dg_opt, color = :blue, label = "Optimized")
ax.xticks = (1:length(rids), rids)
fig[1, 2] = Legend(fig, ax, "ΔG'", framevisible = false)
fig

#md # !!! tip "Tip: Thermodynamic variability"
#md #       Much like flux variability, thermodynamic constraints in a model are also
#md #       degenerate. Check out [`max_min_driving_force_variability`](@ref) for ways
#md #       to explore the thermodynamic solution space.


text:
#aliases
get_cell_lines(experiment_id::String) = to_json(collect(values(get_experiment(experiment_id).cell_lines)))
get_drugs(experiment_id::String) = to_json(collect(values(get_experiment(experiment_id).drugs)))
get_genes(experiment_id::String) = to_json(union(collect(values(get_experiment(experiment_id).genes)), collect(values(get_experiment(experiment_id).genes_by_hgnc)), collect(values(get_experiment(experiment_id).genes_by_ensembl))))
get_proteins(experiment_id::String) = to_json(get_experiment(experiment_id).proteins)
get_pathways(experiment_id::String) = to_json(collect(values(get_experiment(experiment_id).pathway_information)))
get_outcomes(experiment_id::String) = to_json(collect(values(get_experiment(experiment_id).results)))

get_cell_line(experiment_id::String, cell_line_id::String) = to_json(get_experiment(experiment_id).cell_lines[cell_line_id])
get_drug(experiment_id::String, drug_id::String) = to_json(get_experiment(experiment_id).drugs[drug_id])
get_gene(experiment_id::String, gene_id::String) = to_json(get_experiment(experiment_id).genes[parse(Int64, gene_id)])
get_protein(experiment_id::String, protein_id::String) = to_json(get_experiment(experiment_id).proteins[protein_id])
get_pathway(experiment_id::String, pw_id::String) = to_json(get_experiment(experiment_id).pathway_information[pw_id])

function get_outcome(experiment_id::String, drug_id::String)
    experiment = get_experiment(experiment_id)
    d = experiment.drugs[drug_id]
    to_json(experiment.results[d])
end

function to_json(data)
    response = Dict{String, Any}()
    response["message"] = "JSON serialization not yet implemented for $(typeof(data))"
    response
end

function to_json(experiment_id::String)
    data = get_experiment(experiment_id)
    response = Dict{String, Any}()
    response["results"] = to_json(collect(values(data.results)))
    response["cell_lines"] = to_json(collect(values(data.cell_lines)))
    response["drugs"] = to_json(collect(values(data.drugs)))
    response["genes"] = to_json(union(collect(values(data.genes)), collect(values(data.genes_by_hgnc)), collect(values(data.genes_by_ensembl))))
    response["proteins"] = to_json(data.proteins)
    response["views"] = map(v -> string(v), data.views)
    response["pathways"] = to_json(collect(values(data.pathway_information)))
    response
end

function to_json(data::Vector)
    ret = Vector{Dict{String, Any}}()
    for obj in data
        push!(ret, to_json(obj))
    end
    ret
end

function to_json(data::CellLine)
    response = Dict{String, Any}()
    response["id"] = data.id
    response["cancer_type"] = data.cancer_type
    response["views"] = Dict{String, Any}()
    for (view_type, view_data) in data.views
        response[string(view_type)] = to_json(view_data)
    end
    response["in_test_set"] = data.in_test_set
    response
end

function to_json(data::Gene)
    response = Dict{String, Any}()
    if data.entrez_id != -1 response["entrez_id"] = data.entrez_id end
    if data.hgnc_id != "" response["hgnc_id"] = data.hgnc_id end
    if data.ensembl_id != "" response["ensembl_id"] = data.ensembl_id end
    response["cancer_gene"] = data.cancer_gene
    response
end

function to_json(data::Protein)
    response = Dict{String, Any}()
    response["hgnc_id"] = data.hgnc_id
    response["antibody_validated"] = data.antibody_validated
end

function to_json(data::Drug)
    response = Dict{String, Any}()
    response["id"] = data.id
    if length(data.affected_genes) > 0
        response["affected_genes"] = Vector{Dict{String, Any}}()
        for gene in data.affected_genes
            push!(response["affected_genes"], to_json(gene))
        end
    end
    if data.chemical_structure != ""
        response["chemical_structure"] = data.chemical_structure
    end
    response
end

function to_json(data::Outcome)
    response = Dict{String, Any}()
    response["drug_id"] = data.drug_id
    response["outcome_values"] = Dict{String, Any}()
    for (cl, val) in data.outcome_values
        response["outcome_values"][cl.id] = val
    end
    response["outcome_type"] = data.outcome_type
    response
end

function to_json(data::Pathway)
    response = Dict{String, Any}()
    response["id"] = data.id
    response["name"] = data.name
    response["genes"] = Vector{Dict{String, Any}}()
    for gene in data.genes
        push!(response["genes"], to_json(gene))
    end
    response
end


text:
module RootedTrees

const IDInt = UInt16    # Maximum of 65535

#======================#
#         CORE
#======================#

abstract type AbstractTree end
abstract type AbstractRTree <: AbstractTree end
abstract type AbstractNode end
abstract type AbstractRNode <: AbstractNode end

export AbstractTree
export AbstractRTree
export AbstractNode
export AbstractRNode

include("Branch.jl")

include("RNode.jl")
export
    RNode,
    getid,
    brlength, brlength!,
    getlabel, setlabel!,
    parent, children,
    indegree, outdegree,
    isouter, isinner

include("RTree.jl")
export
    RTree,
    getroot,
    getnode, getnodes,
    nnode, ninner, nouter,
    createtree

include("link.jl")
export link!, cut!, swap!

#------end of CORE------#

include("tree_builders.jl")
export symmetric_tree

include("traverse.jl")
export
    pretraverse, posttraverse,
    preorder, postorder

include("reroot.jl")
export reroot!

include("ancestors.jl")
export
    path_to_root,
    findmrca,
    mrca_matrix

include("read_newick.jl")
export parsenewick

include("write_newick.jl")
export newick, print_newick

include("preplot.jl")

include("plot_recipe.jl")
export plot, plot!

end # module

max_stars_repo_name: eascarrunz/RootedTrees

max_issues_repo_name: eascarrunz/RootedTrees

max_forks_repo_name: eascarrunz/RootedTrees


text:
using Test
using PolynomialBasis
using IntervalArithmetic
# using Revise
using ImplicitDomainQuadrature

IDQ = ImplicitDomainQuadrature
PB = PolynomialBasis

function allapprox(v1, v2)
    return all(v1 .≈ v2)
end

f(x) = (x - 1.0) * (x - 1.5)
g(x) = (x - 0.5) * (x - 1.25)
r = IDQ.roots_and_ends([f, g], 0.0, 2.0)
testr = [0.0, 0.5, 1.0, 1.25, 1.5, 2.0]
@test allapprox(r, testr)

@test IDQ.extend(1.0, 1, 2.0) ≈ [2.0, 1.0]
@test IDQ.extend(1.0, 2, 2.0) ≈ [1.0, 2.0]
@test_throws ArgumentError IDQ.extend(1.0, 3, 2.0)

@test allapprox(IDQ.extend([1.0], 1, 2.0), [2.0, 1.0])
@test_throws ArgumentError IDQ.extend([1.0, 2.0], 1, 2.0)

x0 = [1.0]
x = [
    2.0 3.0 4.0 5.0
    6.0 7.0 8.0 9.0
]
@test_throws ArgumentError IDQ.extend(x0, 1, x)

x0 = [1.0]
x = [2.0 3.0 4.0 5.0]
testx = [
    2.0 3.0 4.0 5.0
    1.0 1.0 1.0 1.0
]
@test allapprox(IDQ.extend(x0, 1, x), testx)

x0 = [1.0]
x = [2.0 3.0 4.0 5.0]
testx = [
    1.0 1.0 1.0 1.0
    2.0 3.0 4.0 5.0
]
@test allapprox(IDQ.extend(x0, 2, x), testx)

x0 = [1.0]
x = [2.0 3.0 4.0 5.0]
@test_throws ArgumentError IDQ.extend(x0, 3, x)

f2(x) = x[2]
P = InterpolatingPolynomial(1, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
@test IDQ.height_direction(x -> gradient(P, x), [0.0, 0.0]) == 2
box = IntervalBox(-1..1, 2)
@test IDQ.height_direction(x -> gradient(P, x), box) == 2

flag, s = IDQ.is_suitable(2, x -> gradient(P, x), box)
@test flag == true
@test s == 1

f2(x) = (x[2] + 0.5) * (x[2] - 0.5)
P = InterpolatingPolynomial(1, 2, 2)
coeffs = [f2(P.basis.points[:, i]) for i = 1:size(P.basis.points)[2]]
update!(P, coeffs)
flag, s = IDQ.is_suitable(2, x -> gradient(P, x), box)
@test flag == false

@test IDQ.sign(1, 1, true, -1) == -1
@test IDQ.sign(1, -1, false, -1) == -1
@test IDQ.sign(1, -1, false, 1) == 0

function testcubic(coeffs, v)
    a, b, c, d, e, f, g, h, i, j = coeffs
    x, y = v
    return a * x^3 +
           b * x^2 * y +
           c * x * y^2 +
           d * y^3 +
           e * x^2 +
           f * x * y +
           g * y^2 +
           h * x +
           i * y +
           j
end

function testcubicgrad(coeffs, v)
    a, b, c, d, e, f, g, h, i, j = coeffs
    x, y = v

    fx = 3 * a * x^2 + 2 * b * x * y + c * y^2 + 2 * e * x + f * y + h
    fy = b * x^2 + 2 * c * x * y + 3 * d * y^2 + f * x + 2 * g * y + i
    return [fx, fy]
end

coeffs = [10,12,6,17,9,5,2,1,8,11.0]
poly = InterpolatingPolynomial(1,2,3)
points = PolynomialBasis.interpolation_points(PolynomialBasis.basis(poly))
interpcoeffs = vec(mapslices(x->testcubic(coeffs,x),points,dims=1))
update!(poly,interpcoeffs)

testp = rand(2,10)
vals = mapslices(poly,testp,dims=1)
testvals = mapslices(x->testcubic(coeffs,x),testp,dims=1)
@test allapprox(vals,testvals)

interpgrad = InterpolatingPolynomial(2,2,3)
IDQ.update_interpolating_gradient!(interpgrad,poly)
vals = mapslices(interpgrad,testp,dims=1)
testvals = mapslices(x->testcubicgrad(coeffs,x),testp,dims=1)
@test allapprox(vals,testvals)

max_stars_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_issues_repo_name: ArjunNarayanan/ImplicitDomainQuadrature

max_forks_repo_name: ArjunNarayanan/ImplicitDomainQuadrature


text:
export ft_rescan

function ft_rescan(ft_handle::Culong)
  ft_status = ccall((:FT_Rescan, d2xx),Cuint,(Culong,),ft_handle)
  checkstatus(ft_status)
  return nothing
end


text:
function H2basis(z::Number, k::Int64)
    return 1/(√π*(z+1.0im)) * h(z)^k
end

function hardy_expand(z::Number, H::Int64, param::Matrix)
    res = 0.
    for k = 1:H
        Bk = H2basis(z, k)
        res += param[k,1] * Bk + param[k,2] * Bk'
    end
    return res
end

function loss(Aω::Vector)

end

# 2nd order derivative


text:
"""tests if point p1 and p2 are on the same side of line ab"""
function same_side(p1::Vector{U},p2::Vector{U}, a::Vector{U}, b::Vector{U}) where U <: AbstractFloat
    cp1 = cross(b-a, p1-a)
    cp2 = cross(b-a, p2-a)
    dot(cp1, cp2) >= 0
end

""" tests if a point p is Left|On|Right of an infinite line ab
- >0 for p left of the line ab
- =0 for p on the line
- <0 for p right of the line
"""
function is_left(p::Vector{U}, a::Vector{U}, b::Vector{U}) where U <: Real
    (b[1]-a[1])*(p[2]-a[2]) - (p[1]-a[1])*(b[2]-a[2])
end

"""calculates 2D convec hull for set of points
returns points which determine the hull and its indices in the input list

based on Andrew's Monotone Chain Algorithm (Andrew, 1979), explained here:
http://geomalgorithms.com/a10-_hull-1.html"""
function get_hull(P::Vector{Vector{T}}) where T <: Real
    np = length(P)
    # sort indices, first by x then y
    pm = zeros(T,(np, 3))
    [pm[i,:] = [P[i][1], P[i][2], i] for i in 1:np] # to 2d array
    pm = sortrows(pm, by=x->(x[1],x[2]))
    P = [pm[i,1:2] for i in 1:np]
    idx = Int[pm[i,3] for i in 1:np]
    # get indices of points with 1st x min or max and 2nd y min or max
    min_min = 1 # min x, min y
    i = 0 ## initiate so it is known outside loop, arghh...
    for i in 2:np
        P[i][1] != P[1][1] && break
    end
    min_max = i-1
    for i in np-1:-1:1
        P[i][1] != P[np][1] && break
    end
    max_min = i+1
    max_max = np
    # initiate convex hull stack with point indices
    stack = Int[]

    ## LOWER HULL
    push!(stack, min_min)
    top = 1 # no. of points in stack
    for i in 2:max_min-1 # loop through points with increasing x sequence
        is_left(P[i], P[min_min], P[max_min]) >= 0 && continue # ignore P[i] above or on the lower line
        while top >= 2 # at least two points in stack
            is_left(P[i], P[stack[top-1]], P[stack[top]]) > 0 && break # P[i] in hull
            pop!(stack)
            top -= 1
        end
        push!(stack, i)
        top += 1
    end
    push!(stack, max_min)
    top += 1

    ## UPPER HULL
    # if  distinct xmax points push max_min point onto stack
    if max_max != max_min
        push!(stack, max_max)
        top += 1
    end
    bot = top # bottom point of upper hull
    for i in np-1:-1:2 # loop through points with decreasing x sequence
        is_left(P[i], P[max_max], P[min_max]) >= 0 && continue # ignore P[i] below or on the upper line
        while bot > top # at least two points in upper stack
            is_left(P[i], P[stack[top-1]], P[stack[top]]) > 0 && break # P[i] in hull
            pop!(stack)
            top -= 1
        end
        push!(stack, i)
        top += 1
    end
    push!(stack, min_max)
    top += 1

    # return points of hull and indices
    P[stack], idx[stack]
end

"""calculates bbox coordinates for set of points
returns the list with coordinates the indices of the points which determine the bbox"""
function get_bbox(P::Vector{Vector{T}}) where T <: Real
    @assert length(P) >= 2
    xmin, xmax, ymin, ymax = Inf, -Inf, Inf, -Inf
    xmini, xmaxi, ymini, ymaxi = 0, 0, 0, 0  # declare indices of points that make up bbox
    i = 0
    @inbounds for p in P
        addp = false
        i += i  # index of point in vector P
        if p[1] < xmin
            xmin = p[1]
            xmini = i
        end
        if p[1] > xmax
            xmax = p[1]
            xmaxi = i
        end
        if p[2] < ymin
            ymin = p[2]
            ymini = i
        end
        if p[2] > ymax
            ymax = p[2]
            ymaxi = p
        end
    end

    # return bbox and indices that make up bbox
    [xmin,xmax,ymin,ymax], [xmini, xmaxi, ymini, ymaxi]
end

function triangle_area(a::Vector{U}, b::Vector{U}, c::Vector{U}) where U <: AbstractFloat
    abs(a[1]*(b[2]-c[2]) + b[1]*(c[2]-a[2]) + c[1]*(a[2]-b[2])) / 2.
end

"""check if point p [x,y] in triangle made out of points a,b,c"""
function in_triangle(p::Vector{U}, a::Vector{U},b::Vector{U},c::Vector{U}) where U <: AbstractFloat
    p,a,b,c = [[x[1],x[2],0.] for x in [p,a,b,c]]
    same_side(p,a, b,c) && same_side(p,b, a,c) && same_side(p,c, a,b)
end

"""check if point (x,y) in traingle made out of points vector"""
function in_triangle(x::AbstractFloat, y::AbstractFloat, points::Vector{Vector{U}}) where U <: AbstractFloat
    length(points)
    npoints != 3 && warn("more than three points in list, check based on first three points only")
    a,b,c = points[1], points[2], points[3]
    xypoint_in_triangle([x, y], a, b,c)
end

"""check if point (x,y) in bbox"""
function in_bbox(x::AbstractFloat, y::AbstractFloat, points::Vector{Vector{U}}) where U <: AbstractFloat
    xmin,xmax,ymin,ymax = get_bbox(points)[1]
    xmin <= x <= xmax && ymin <= y <= ymax
end

"""check if point in hull of all points based on hull computation"""
function in_hull(x::Real, y::Real, points::Vector{Vector{U}}) where U <: Real
    npoints = length(points)
    @assert npoints >= 3 ["point vector should contain at least three points"]
    push!(points, [x,y])
    p_hull, idx_hull = get_hull(points)
    !in(npoints+1, idx_hull)
end


text:
using Vehicles;

malcolm = Contact("Malcolm","malcolm@abc.net","07777555999");
myCar = Ford(malcolm, "Model T", 1000, 50.0);
myBike = Scooter(malcolm, "Vespa", 125, 35.0);

james = Contact("James","james@abc.net","07777666888");
jmCar = BMW(james,"Series 500", 3200, 125.0);
jmBoat = Yacht(james,"Oceanis 44",14.6);
jmBike = MotorBike(james, "Harley", 850, 120.0);

david = Contact("David","dave@abc.net","07777222444");
dvCar = VW(david,"Golf", "diesel", "red", 1800, 85.0);
dvBoat = Speedboat(david,"Sealine","petrol", 600, 45.0, 8.2);

-----------------------

dsBoat.owner

msCar
msCar.owner = david;
msCar

s = [myCar, jmCar, dvCar]
3-element Array{Car,1}:
 Ford(Contact("Malcolm","malcolm@abc.net","07777555999"),"Model T","Petrol","Black",1000,50.0)
 BMW(Contact("James","james@abc.net","07777666888"),"Series 500","Petrol","Blue",3200,125.0)
 VW(Contact("David","dave@abc.net","07777222444"),"Golf","diesel","red",1800,85.0)

for c in s
    who = c.owner.name
    model = c.make
    println("$who has a $model car")
end

s =[jmCar, jmBike, jmBoat]
for c in s
    who = c.owner.name
    model = c.make
    println("$who has a $model vehicle")
end

print "James owns these vehicles: "
for c in s
    print("$(s.make) ")
end
println()

Series 500
Oceanis 44

max_issues_repo_path: Module 3/Chapter03/vehicles-ex.jl

max_forks_repo_path: Module 3/Chapter03/vehicles-ex.jl


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    Remainder([total])

The transform that takes a table with columns `x₁, x₂, …, xₙ`
and returns a new table with an additional column containing the
remainder value `xₙ₊₁ = total .- (x₁ + x₂ + ⋯ + xₙ)` If the `total`
value is not specified, then default to the maximum sum across rows.
"""
struct Remainder <: Transform end

isrevertible(::Type{Remainder}) = true

function apply(::Remainder, table)
  # TODO
end

function revert(::Remainder, newtable, cache)
  # TODO
end

function reapply(::Remainder, table, cache)
  # TODO
end


text:
@info "Testing common.jl"

# Dimensionality Checks

@testset "check_dims(X; dims)" begin
    n = 20
    p = 5
    for T in (Float32, Float64)
        @test_throws ArgumentError DA.check_dims(zeros(T, p, p), dims=0)
        @test_throws ArgumentError DA.check_dims(zeros(T, p, p), dims=3)

        @test (n, p) == DA.check_dims(zeros(T, n, p), dims=1)
        @test (n, p) == DA.check_dims(zeros(T, p, n), dims=2)
        @test (n, p) == DA.check_dims(transpose(zeros(T, n, p)), dims=2)
    end
end

# Data Validation

@testset "validate_priors(π)" begin
    k = 10
    for T in (Float32, Float64)
        # Check probability domain
        @test_throws DomainError DA.validate_priors(T[0.5; 0.5; 0.0])
        @test_throws DomainError DA.validate_priors(T[1.0; 0.0; 0.0])

        # Check totals
        @test_throws ArgumentError DA.validate_priors(T[0.3; 0.3; 0.3])
        @test_throws ArgumentError DA.validate_priors(T[0.4; 0.4; 0.4])

        # Test valid π
        π = rand(T, k)
        π ./= sum(π)
        #broadcast!(/, π, π, sum(π))

        @test k == DA.validate_priors(π)
    end
end

# Class operations

@testset "class_counts(nₘ, y)" begin
    m = 3
    nₖ = 300

    y = repeat(vec(1:m), outer=nₖ)

    # Test bounds
    y_tst = copy(y)

    y_tst[1] = 0
    @test_throws BoundsError DA.class_counts!(Vector{Int}(undef, m), y_tst)

    y_tst[1] = m + 1
    @test_throws BoundsError DA.class_counts!(Vector{Int}(undef, m), y_tst)

    # Test computation
    @test DA.class_counts!(Vector{Int}(undef, m), y) == [nₖ for i = 1:m]
end

@testset "_class_statistics!(M, nₘ, X, y)" begin
    nₘ = [45, 55, 100]
    m = length(nₘ)
    n = sum(nₘ)
    p = 5

    for T in (Float32, Float64)
        X, y, M = random_data(T, nₘ, p)
        Xtt = transpose(copy(transpose(X)))
        Mtt = transpose(copy(transpose(M)))

        for X_test in (X, Xtt)
            nₘ_test = copy(nₘ)
            # test predictor dimensionality
            @test_throws DimensionMismatch DA._class_statistics!(zeros(T,2,p+1), nₘ_test, X_test,  y)
            @test_throws DimensionMismatch DA._class_statistics!(zeros(T,2,p-1), nₘ_test, X_test,  y)

            # test observation dimensionality
            @test_throws DimensionMismatch DA._class_statistics!(similar(M), nₘ_test, X_test, zeros(Int,n+1))
            @test_throws DimensionMismatch DA._class_statistics!(similar(M), nₘ_test, X_test, zeros(Int,n-1))

            # test count dimensionality
            @test_throws DimensionMismatch DA._class_statistics!(similar(M), zeros(Int,m+1), X_test, y)
            @test_throws DimensionMismatch DA._class_statistics!(similar(M), zeros(Int,m-1), X_test, y)

            # test indexing of class_statistics
            y_test = copy(y)
            y_test[1] = 0
            @test_throws BoundsError DA._class_statistics!(similar(M), nₘ_test, X_test, y_test)

            y_test[1] = m+1
            @test_throws BoundsError DA._class_statistics!(similar(M), nₘ_test, X_test, y_test)
        end

        # test mean computation
        for (M_test, X_test) in ((copy(M), X), (deepcopy(Mtt), Xtt))
            nₘ_test = zeros(Int, m)

            M_res, nₘ_res = DA._class_statistics!(M_test, nₘ_test, X, y)

            @test nₘ_res === nₘ_test
            @test nₘ_res == nₘ

            @test M_res === M_test
            @test isapprox(M, M_test)
        end
    end
end

### Regularization

@testset "regularize!(Σ₁, Σ₂, λ)" begin
    for T in (Float32, Float64)
        S = zeros(T, 3, 3)
        @test_throws DimensionMismatch DA.regularize!(zeros(T, 2, 3), S, T(0))
        @test_throws DimensionMismatch DA.regularize!(S, zeros(T, 2, 3), T(0))
        @test_throws DimensionMismatch DA.regularize!(S, zeros(T, 2, 2), T(0))

        @test_throws DomainError DA.regularize!(S, S, T(0) - eps(T(0)))
        @test_throws DomainError DA.regularize!(S, S, T(1) + eps(T(1)))

        S1 = cov(rand(T, 10, 3))
        S2 = cov(rand(T, 10, 3))
        for λ in (T(0), T(0.25), T(0.5), T(0.75), T(1))
            S1_test = copy(S1)
            S2_test = copy(S2)

            DA.regularize!(S1_test, S2_test, λ)

            @test isapprox(S1_test, (1-λ)S1 + λ*S2)
            @test S2_test == S2
        end
    end
end

@testset "regularize!(Σ, γ)" begin
    n = 10
    p = 3
    for T in (Float32, Float64)
        @test_throws DimensionMismatch DA.regularize!(zeros(T, 2, 3), T(0))
        @test_throws DimensionMismatch DA.regularize!(zeros(T, 3, 2), T(0))

        S = zeros(T, p, p)
        @test_throws DomainError DA.regularize!(S, T(0) - eps(T(0)))
        @test_throws DomainError DA.regularize!(S, T(1) + eps(T(1)))

        S = cov(rand(T, n, p))
        for γ in (T(0), T(0.25), T(0.5), T(0.75), T(1))
            S_test = copy(S)

            DA.regularize!(S_test, γ)

            @test isapprox(S_test, (1-γ)S + γ*(tr(S)/p)I)
        end
    end
end


text:
struct CentroidalTrajectoryProblem
    model::JuMP.Model

    # Axes
    c_coeffs
    f_coeffs
    contacts
    regions
    pieces
    coords
    coords2d

    # Variables
    c_vars
    f_vars
    f̄_vars
    r_vars
    p_vars
    p̄_vars
    Δts
    z_vars

    # Other data
    normals
    friction_cone_quadratic_constraints
end

module ObjectiveTypes
export ObjectiveType
@enum ObjectiveType begin
    FEASIBILITY
    MIN_TIME
    MIN_EXCURSION
    MAX_HEIGHT
end
end
using .ObjectiveTypes

function CentroidalTrajectoryProblem(optimizer_factory::JuMP.OptimizerFactory,
        region_data::AbstractVector{<:ContactRegion},
        c0::AbstractVector{<:Number},
        ċ0::AbstractVector{<:Number},
        contacts0::AbstractVector{<:Union{Pair{<:ContactRegion, <:AbstractVector}, Nothing}};
        cf = nothing,
        com_degree = 3,
        cop_degree = com_degree,
        num_pieces = 2,
        g = SVector(0.0, 0.0, -9.81),
        max_cop_distance = 0.1,
        min_com_to_contact_distance,
        max_com_to_contact_distance,
        min_inter_contact_distance,
        objective_type::ObjectiveType = FEASIBILITY,
        min_Δt = 0.6,
        max_Δt = 1.5,
        c_margin_xy = 0.5,
        c_margin_z_min = 0.5,
        c_margin_z_max = 1.2,
        max_force = 3 * norm(g),
        Δrmax = 0.7 # TODO
    )

    c_num_coeffs = com_degree + 1
    f_num_coeffs = c_num_coeffs - 2
    p_num_coeffs = cop_degree + 1
    num_regions = length(region_data)
    num_contacts = length(contacts0)

    # Axes
    c_coeffs = Axis{:c_coeff}(1 : c_num_coeffs)
    f_coeffs = Axis{:f_coeff}(1 : f_num_coeffs)
    p_coeffs = Axis{:r_coeff}(1 : p_num_coeffs)
    contacts = Axis{:contact}(1 : num_contacts)
    regions = Axis{:region}(1 : num_regions)
    pieces = Axis{:piece}(1 : num_pieces)
    coords = Axis{:coord}(SVector(:x, :y, :z))
    coords2d = Axis{:coord}(SVector(:x, :y))

    # Time stuff
    Δt_tolerance = 0.05 # TODO: constructor arg

    model = Model(optimizer_factory)
    optimizer = backend(model)#.optimizer.model

    optimizer_does_sos1 = MOI.supports_constraint(optimizer, MOI.VectorOfVariables, MOI.SOS1{Float64})
    optimizer_does_sos1 || @info "Optimizer does not support SOS1 constraints; using sum <= 1 instead."

    optimizer_does_soc = MOI.supports_constraint(optimizer, MOI.VectorOfVariables, MOI.SecondOrderCone)
    optimizer_does_soc || @info "Optimizer does not support second order cone constraints; using quadratic constraints instead."

    optimizer_does_indicator_constraints = false#optimizer.optimizer.model isa SCIP.Optimizer # TODO: generalize

    # Indexing convention:
    # i: piece index
    # j: contact index
    # k: coordinate index
    # l: coefficient index
    # m: region index

    if min_Δt == max_Δt
        Δts = AxisArray(fill(min_Δt, length(pieces)), pieces)
        Δtsqs = map(x -> x^2, Δts)
    else
        Δts = axis_array_vars(model, i -> "Δt[$i]", pieces)
        Δtsqs = axis_array_vars(model, i -> "Δtsq[$i]", pieces)
        for (Δt, Δtsq) in zip(Δts, Δtsqs)
            set_lower_bound(Δt, min_Δt)
            set_upper_bound(Δt, max_Δt)
            set_lower_bound(Δtsq, min_Δt^2)
            set_upper_bound(Δtsq, max_Δt^2)
            @constraint model Δtsq == Δt^2
        end
        if objective_type == ObjectiveTypes.MIN_TIME
            @objective model Min sum(Δts)
        end
    end
    # Δts = AxisArray([ 1.0910477104335914
    #     0.6545247110231425
    #     0.6014271582441553
    #     0.848420003941298 ], pieces)
    # Δtsqs = map(x -> x^2, Δts)

    # Contact indicators. z[pieces(i), contacts(j), regions(m)] == 1 implies that
    # during piece i, contact j is in contact with region m.
    z_vars = axis_array_vars(model, (i, j, m) -> "z[p$i, c$j, r$m]", pieces, contacts, regions)
    foreach(set_binary, z_vars)

    # Contact normals
    Rs = AxisArray(map(data -> data.transform.linear, region_data), regions)
    ns = map(R -> R[:, 3], Rs)

    # Contact location extrema (local coordinates)
    r̄_extrema = AxisArray(map(region -> polyhedron_extrema(polyhedron(hrep(region.A, region.b))), region_data), regions)
    p̄_extrema = map(((p̄_min, p̄_max),) -> (p̄_min .- max_cop_distance, p̄_max .+ max_cop_distance), r̄_extrema)
    r_vertices = mapreduce(vcat, 1 : num_regions) do m
        let region = region_data[m]
            map(Polyhedra.points(polyhedron(hrep(region.A, region.b)))) do p̄_vertex
                region.transform([p̄_vertex; 0])
            end
        end
    end
    r_min, r_max = polyhedron_extrema(polyhedron(vrep(r_vertices)))
    p_min, p_max = r_min .- max_cop_distance, r_max .+ max_cop_distance
    c_min, c_max = p_min - SVector(c_margin_xy, c_margin_xy, -c_margin_z_min), p_max + SVector(c_margin_xy, c_margin_xy, c_margin_z_max)

    # Continuous variables
    c_vars = axis_array_vars(model, (i, k, l) -> "C[p$i, $k, $l]",pieces, coords, c_coeffs)
    f_vars = axis_array_vars(model, (i, j, k, l) -> "F[p$i, c$j, $k, $l]", pieces, contacts, coords, f_coeffs)
    f̄_vars = axis_array_vars(model, (i, j, k, l, m) -> "F̄[p$i, c$j, $k, $l, r$m]", pieces, contacts, coords, f_coeffs, regions)
    r_vars = axis_array_vars(model, (i, j, k) -> "P[p$i, c$j, $k]", pieces, contacts, coords)
    r̄_vars = axis_array_vars(model, (i, j, k, m) -> "P̄[p$i, c$j, $k, r$m]", pieces, contacts, coords2d, regions)
    p_vars = axis_array_vars(model, (i, j, k, l) -> "R[p$i, c$j, $k, $l]", pieces, contacts, coords, p_coeffs)
    p̄_vars = axis_array_vars(model, (i, j, k, l, m) -> "R̄[p$i, c$j, $k, $l, r$m]", pieces, contacts, coords2d, p_coeffs, regions)

    #@variable model objval
    if objective_type == ObjectiveTypes.MIN_EXCURSION
        # @constraint model objval >= sum((c_vars[pieces(i), c_coeffs(l)] - c0) ⋅ (c_vars[pieces(i), c_coeffs(l)] - c0) for i in 1 : num_pieces, l in 1 : c_num_coeffs)
        @objective model Min sum(x -> x ⋅ x, (c_vars[pieces(i), c_coeffs(l)] - c0) for i in 1 : num_pieces, l in 1 : c_num_coeffs)
    elseif objective_type == ObjectiveTypes.MAX_HEIGHT
        @objective model Max sum(c_vars[c_coeffs(3)])
    else
        @objective model Min 0
    end
    # @objective model Max sum(z_vars)

    cprev = nothing
    c′prev = nothing
    friction_cone_quadratic_constraints = []

    for i in 1 : num_pieces
        piece = pieces(i)

        # CoM and derivatives
        c = [BezierCurve(c_vars[piece, coords(k)]...) for k in 1 : length(coords)]
        c′ = derivative.(c)
        c′′ = derivative.(c′)

        # CoM bounds
        for l in 1 : c_num_coeffs
            cpoint = map(x -> x.coeffs[l], c)
            set_lower_bound.(cpoint, c_min)
            set_upper_bound.(cpoint, c_max)
        end

        # Contact/region assignment constraints
        for j in 1 : num_contacts
            contact = contacts(j)

            # Each contact can be assigned to at most one region
            if optimizer_does_sos1
                @constraint model z_vars[piece, contact] in MOI.SOS1(collect(Float64, 1 : num_regions))
            else
                @constraint model sum(z_vars[piece, contact]) <= 1
            end

            r = r_vars[piece, contact]
            set_lower_bound.(r, r_min)
            set_upper_bound.(r, r_max)
            if i == 1
                # Initial contact assignment.
                if contacts0[j] === nothing
                    @constraint model sum(z_vars[piece, contact]) == 0
                else
                    region_data0, r0 = contacts0[j]
                    m = findfirst(isequal(region_data0), region_data)
                    z_var = z_vars[piece, contact, regions(m)]
                    fix(z_var, 1, force=true)
                    # unset_binary(z_var) # to make Alpine happpy
                    fix.(r, r0, force=true)
                end
            else
                # Each contact must be unassigned for one piece before it can be reassigned to a region.
                # Let Δzᵢ,ⱼ,ₘ = zᵢ,ⱼ,ₘ - zᵢ₋₁,ⱼ,ₘ. Then there are three cases for ∑ₘ |Δzᵢ,ⱼ,ₘ|:
                # 1) ∑ₘ |Δzᵢ,ⱼ,ₘ| = 0: the assignment during piece i is the same as the assignment during piece i - 1
                # 2) ∑ₘ |Δzᵢ,ⱼ,ₘ| = 1: unassigned during piece i and assigned during piece i - 1 or vice versa
                # 3) ∑ₘ |Δzᵢ,ⱼ,ₘ| = 2: different assignment during piece i and piece i - 1
                # The third case must be prevented, so we want ∑ₘ |Δzᵢ,ⱼ,ₘ| ≤ 1.
                # This is an ℓ₁-norm constraint.
                Δz = z_vars[piece, contact] - z_vars[pieces(i - 1), contact]
                constrain_l1_norm(model, Δz, 1; add_bounds=true)

                # Two swing segments in a row is not allowed.
                @constraint model sum(z_vars[piece, contact]) + sum(z_vars[pieces(i - 1), contact]) >= 1

                # Contact reference position r may only change when the contact is not assigned to a region,
                # i.e., when ∑ₘ zᵢ,ⱼ,ₘ = 0.
                Δr = r - r_vars[pieces(i - 1), contact]
                # @constraint model sum(x -> x^2, Δp) <= (1 - sum(z_vars[piece, contact])) * Δpmax^2
                z_no_contact = @variable model
                set_binary(z_no_contact)
                @constraint model z_no_contact == 1 - sum(z_vars[piece, contact])
                @constraint model  Δr .<= z_no_contact * Δrmax
                @constraint model -Δr .<= z_no_contact * Δrmax
            end
        end

        # Initial conditions
        if i == 1
            fix.(map(x -> x.coeffs[1], c), c0, force=true)
            @constraint model map(x -> x(0), c′) .== ċ0 .* Δts[i]
        end

        # Final conditions
        if i == num_pieces
            if cf !== nothing
                fix.(map(x -> x.coeffs[end], c), cf, force=true)
            end
            @constraint model map(x -> x(1), c′) .== 0
            @constraint model map(x -> x(1), c′′) .== 0
        end

        # CoM continuity
        if i > 1
            @constraint model map(x -> x(1), cprev) .== map(x -> x(0), c)
            @constraint model map(x -> x(1), c′prev) .* Δts[i] .== map(x -> x(0), c′) .* Δts[i - 1]
        end

        # Forces (global)
        fs = [[BezierCurve(f_vars[piece, contacts(j), coords(k)]...) for k in coords.val] for j in contacts.val]

        # CoPs (global)
        ps = [[BezierCurve(p_vars[piece, contacts(j), coords(k)]...) for k in coords.val] for j in contacts.val]

        # Dynamics
        ftot = sum(fs)
        constrain_poly_equal.(model, c′′, identity((g + ftot) .* Δtsqs[i]))

        # Torque balance
        τ = sum(Polynomial.(ps[j]) × Polynomial.(fs[j]) for j = 1 : num_contacts) -
            Polynomial.(c) × Polynomial.(ftot)
        τ = map(x -> Polynomial(simplify.(model, x.coeffs)), τ)
        constrain_poly_equal.(model, τ, 0)

        # CoP and contact position constraints
        # TODO: could derive tighter big-M from constraints on p and p̄, but it's unclear whether that will actually help performance,
        # because changing Mr to
        Mr = 2.0
        for j in 1 : num_contacts
            contact = contacts(j)
            r = r_vars[piece, contact]
            for m in 1 : num_regions
                region = regions(m)
                r̄ = r̄_vars[piece, contact, region]
                r̄_min, r̄_max = r̄_extrema[region]
                set_lower_bound.(r̄, r̄_min)
                set_upper_bound.(r̄, r̄_max)
                A = region_data[m].A
                b = region_data[m].b
                transform = region_data[m].transform
                z = z_vars[piece, contact, region]

                # Contact position constraints
                @constraint model A * r̄ .<= b
                r_aux = @variable model [1 : 3]
                @constraint model r_aux .== r - transform([r̄; 0])
                if optimizer_does_indicator_constraints
                    for i in eachindex(r_aux)
                        MOI.add_constraint(optimizer, MOI.VectorOfVariables([z; r_aux]), SCIP.IndicatorSet(setindex!(zeros(3), +1, i), 0.0))
                        MOI.add_constraint(optimizer, MOI.VectorOfVariables([z; r_aux]), SCIP.IndicatorSet(setindex!(zeros(3), -1, i), 0.0))
                    end
                else
                    @constraint model  r_aux .<= Mr * (1 - z)
                    @constraint model -r_aux .<= Mr * (1 - z)
                end

                # CoP constraints
                p̄_min, p̄_max = p̄_extrema[region]
                for l = 1 : p_num_coeffs
                    coeff = p_coeffs(l)
                    p = p_vars[piece, contact, coeff]
                    set_lower_bound.(p, p_min)
                    set_upper_bound.(p, p_max)
                    p̄ = p̄_vars[piece, contact, coeff, region]
                    set_lower_bound.(p̄, p̄_min)
                    set_upper_bound.(p̄, p̄_max)
                    # @constraint model A * p̄ .<= b
                    p_aux = @variable model [1 : 3]
                    @constraint model p_aux .== p - transform([p̄; 0])
                    if optimizer_does_indicator_constraints
                        for i in eachindex(p_aux)
                            MOI.add_constraint(optimizer, MOI.VectorOfVariables([z; p_aux]), SCIP.IndicatorSet(setindex!(zeros(3), +1, i), 0.0))
                            MOI.add_constraint(optimizer, MOI.VectorOfVariables([z; p_aux]), SCIP.IndicatorSet(setindex!(zeros(3), -1, i), 0.0))
                        end
                    else
                        @constraint model  p_aux .<= Mr * (1 - z)
                        @constraint model -p_aux .<= Mr * (1 - z)
                    end
                    if max_cop_distance == 0
                        @constraint model p .== r
                    else
                        @constraint model sum(x -> x^2, p - r) <= max_cop_distance^2
                        # constrain_l1_norm(model, r - p, √2 / 2 * max_cop_distance; add_bounds=true)
                    end
                end
            end
        end

        # Contact force/torque constraints
        for j in 1 : num_contacts
            contact = contacts(j)
            μmax = 0.0
            for m in 1 : num_regions
                region = regions(m)
                μ = region_data[m].μ
                μmax = max(μ, μmax)
                μrot = region_data[m].μrot
                R = Rs[region]
                n = ns[region]
                z = z_vars[piece, contact, region]
                for l in 1 : f_num_coeffs
                    coeff = f_coeffs(l)
                    f̄ = f̄_vars[piece, contact, coeff, region]
                    f̄xy = f̄[1 : 2]
                    f̄z = f̄[3]
                    set_lower_bound(f̄z, 0)
                    set_upper_bound(f̄z, max_force)
                    set_lower_bound.(f̄xy, -μ * max_force)
                    set_upper_bound.(f̄xy, +μ * max_force)
                    # auxiliary variable needed for SecondOrderCone constraint:
                    μf̄z = @variable model lower_bound=0 upper_bound=μ * upper_bound(f̄z)
                    @constraint model μf̄z == μ * f̄z
                    if optimizer_does_indicator_constraints
                        z_opposite = @variable model
                        set_binary(z_opposite)
                        @constraint model z_opposite == 1 - z
                        MOI.add_constraint(optimizer, MOI.VectorOfVariables([z_opposite; f̄z]), SCIP.IndicatorSet([1.0], 0.0))
                    else
                        # @constraint model μf̄z <= max_force * μ * z
                        @constraint model f̄z <= max_force * z
                    end
                    if optimizer_does_soc
                        @constraint model [μf̄z; f̄xy] in SecondOrderCone()
                    else
                        friction_cone_quadratic_constraint = @constraint model f̄xy ⋅ f̄xy <= μf̄z^2
                        push!(friction_cone_quadratic_constraints, friction_cone_quadratic_constraint)
                    end
                end
            end
            for l in 1 : f_num_coeffs
                coeff = f_coeffs(l)
                f = f_vars[piece, contact, coeff]
                # set_lower_bound.(f, -max(μmax, 1) * max_force)
                # set_upper_bound.(f, +max(μmax, 1) * max_force)
                @constraint model f .== sum(Rs[regions(m)] * f̄_vars[piece, contact, coeff, regions(m)] for m in 1 : num_regions)
            end
        end

        # CoM kinematic constraints
        for j in 1 : num_contacts
            contact = contacts(j)
            r = r_vars[piece, contact]
            for l in 1 : c_num_coeffs
                cpoint = map(x -> x.coeffs[l], c)
                # @constraint model cpoint[3] - r[3] >= 0.75 # TODO
                @constraint model sum(x -> x^2, cpoint - r) >= min_com_to_contact_distance^2
                # set_lower_bound(cpoint[3], 0.7)
                # @constraint model cpoint - r .<= 1.5
                # @constraint model r - cpoint .<= 1.5
                # @constraint model sum(x -> x^2, cpoint - r) >= 0.5^2 # TODO
                @constraint model sum(x -> x^2, cpoint - r) <= max_com_to_contact_distance^2
            end
        end

        # Inter-contact kinematic constraints
        if i > 1
            for j1 in 1 : num_contacts
                contact1 = contacts(j1)
                r1 = r_vars[piece, contact1]
                for j2 in j1 + 1 : num_contacts
                    contact2 = contacts(j2)
                    r2 = r_vars[piece, contact2]
                    @constraint model sum(x -> x^2, (r1 - r2)) >= min_inter_contact_distance^2
                end
            end
            # FIXME: hack:
            @constraint model r_vars[piece, contacts(1), coords(2)] - r_vars[piece, contacts(2), coords(2)] >= -max_cop_distance
        end

        cprev = c
        c′prev = c′
    end

    CentroidalTrajectoryProblem(model,
        c_coeffs, f_coeffs, contacts, regions, pieces, coords, coords2d,
        c_vars, f_vars, f̄_vars, r_vars, p_vars, p̄_vars, Δts, z_vars,
        ns, friction_cone_quadratic_constraints
    )
end

function disallow_jumping!(problem::CentroidalTrajectoryProblem)
    model = problem.model
    pieces = problem.pieces
    for i in 1 : length(pieces)
        @constraint model sum(problem.z_vars[pieces(i)]) >= 1
    end
end

normals(problem) = problem.normals


text:
using Test

@testset "unique" begin

source_1 = [1,2,2,3,3,3,4]
enum = QueryOperators.query(source_1)

@test collect(QueryOperators.@unique(enum, i->i)) == [1,2,3,4]

source_1 = [1,2,3,4]
enum = QueryOperators.query(source_1)

@test collect(QueryOperators.@unique(enum, i->i)) == [1,2,3,4]

source_1 = [1,-1,2,3,4]
enum = QueryOperators.query(source_1)

@test collect(QueryOperators.@unique(enum, i->abs(i))) == [1,2,3,4]

end


text:
export fibonacci_bound
export fibonacci_count
export unique_products_of_ndigits
export multiset_combination!
export permutation_matrix!

"""
Returns an array of the Fibonacci numbers up to a bound n.
"""
function fibonacci_bound(n::BigInt=1000)::Array{BigInt}
    a = BigInt(1)
    b = BigInt(2)
    fibs = zeros(BigInt, 1)
    fibs[1] = 1
    while a < n
        (a, b) = (b, a + b)
        push!(fibs, a)
    end
    return fibs
end

"""
Returns an array of the first n Fibonacci numbers.
"""
function fibonacci_count(n::Integer=20)::Array{BigInt}
    a = BigInt(1)
    b = BigInt(2)
    fibs = zeros(BigInt, n)
    fibs[1] = 1
    for i in 2:n
        (a, b) = (b, a + b)
        @inbounds fibs[i] = a
    end
    return fibs
end

"""
Return a list of products of numbers with "ndigit" digits.
"""
function unique_products_of_ndigits(n_digits::Integer)::Array{Integer}
    maxnum = 10^n_digits - 1
    minnum = 10^(n_digits - 1)
    return [i * j for i in minnum:maxnum for j in minnum:maxnum if i < j]
end

"""
Return a multiset (combinations with replacement, order doesn't matter).  Uses
a recursive solution to build up a full set of indexes.
"""
function multiset_combination!(storage::Array{Int8,2},
                               digit_index::Integer,
                               index::Integer,
                               rstart::Integer,
                               rend::Integer)
    if digit_index == 0
        storage[:, index + 1] = storage[:, index]
        return index + 1
    end

    for i in rstart:rend
        storage[digit_index, index] = i
        index = multiset_combination!(storage, digit_index - 1, index, i, rend)
    end
    return index
end

"""
Create array of permutations for which valid_func is true.
"""
function permutation_matrix!(elementlist::Array, matrix::Array, valid_func, k::Integer=1)
    if k == length(elementlist)
        if valid_func(elementlist)
            append!(matrix, elementlist)
        end
    else
        for i in k:length(elementlist)
            elementlist[k], elementlist[i] = elementlist[i], elementlist[k]
            permutation_matrix!(elementlist, matrix, valid_func, k + 1)
            elementlist[k], elementlist[i] = elementlist[i], elementlist[k]
        end
    end
end

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# Cheke, L. G. and Clayton, N. S., Eurasian jays (Garrulus glandarius) overcome
# their current desires to anticipate two distinct future needs and plan for
# them appropriately, Biology Letters 8:171-175, 2011,
# http://dx.doi.org/10.1098/rsbl.2011.0909
#
# id-bird mapping: 1: Hoy, 2: Wiggins, 3: Hunter, 4: Ainsley
# Wiggins and Ainsley were fed with raisins instead of suet pellets. For
# simplicity I use suet pellets for all birds in this protocol.

function _summarize(::Experiment{:Cheke11_specsat}, data)
	combine(groupby(data, [:prefed, :action, :foodtype]),
				df -> DataFrame(μ = mean(df.items),
                                sem = sem(float.(df.items)),
								n = length(df.items)))
end
function statistical_tests(exp::Experiment{:Cheke11_specsat}, data)
    tests = Test[]
    push!(tests, Test("overall",
                      splitplotanova(data, exp.tests["overall"].locals,
                                     betweenfactors = [],
                                     withinfactors = [:prefed, :foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    for action in ("cache", "eat")
        push!(tests, Test(action,
                          splitplotanova(@where(data, :action .== action),
                                         exp.tests[action].locals,
                                         betweenfactors = [],
                                         withinfactors = [:prefed, :foodtype],
                                         id = :id, data = :items),
                          exp.tests))
    end
    TestCollection(tests)
end
function run!(exp::Experiment{:Cheke11_specsat}, models; N = 4)
	results = DataFrame(prefed = String[], action = String[],
					    foodtype = String[], id = Int[], items = Float64[])
	for id in 1:N
		m = models[id]
		add!(m, MaintenanceDiet)
		for food in (PowderedPeanut, PowderedSuetPellet)
			remove!(m, MaintenanceDiet)
			wait!(m, 2u"hr")
			add!(m, food, 20)
			wait!(m, 15u"minute")
			remove!(m, Any)
			add!(m, Peanut, 40)
			add!(m, SuetPellet, 40)
			tray = Tray(0, rand(1:100))
			add!(m, tray)
			wait!(m, 15u"minute")
			prefed = lowercase(string(food))
			for typ in (Peanut, SuetPellet)
                push!(results, [prefed, "eat", lowercase(string(typ)), id,
								40 - countfooditems(m, typ)])
                push!(results, [prefed, "cache", lowercase(string(typ)), id,
								countcache(tray, typ)])
			end
			remove!(m, Any)
			add!(m, MaintenanceDiet)
			wait!(m, 21.5u"hr")
		end
	end
    results
end

function getfraction_Cheke11_planning(data)
	groups = groupby(data, :foodtype)
	df = join(groups..., on = [:trial, :tray, :id])
	total = df.items .+ df.items_1
	df.fractionprefed = zeros(length(df.items))
	for (i, t) in enumerate(total)
		if t == 0
			df.fractionprefed[i] = NA
		else
			df.fractionprefed[i] = df.items_1[i] / t
		end
	end
	df
end
function _summarize(::Experiment{:Cheke11_planning}, data)
	combine(groupby(data, [:trial, :foodtype, :tray, :action]),
				df -> DataFrame(μ = mean(df.items),
                                sem = sem(float.(df.items)),
								n = length(df.items)))
end
function statistical_tests(exp::Experiment{:Cheke11_planning}, data)
    tests = Test[]
    push!(tests, Test("foodtype & tray",
                      splitplotanova(data,
                                     exp.tests["foodtype & tray"].locals,
                                     betweenfactors = [],
                                     withinfactors = [:tray, :foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    for trial in 1:3
        id = "Tray × foodtype interaction in trial $trial"
        push!(tests, Test(id,
                          splitplotanova(@where(data, :trial .== trial),
                                         exp.tests[id].locals,
                                         betweenfactors = [],
                                         withinfactors = [:tray, :foodtype],
                                         id = :id, data = :items),
                          exp.tests))
    end
    id = "Tray × foodtype interaction in trial 3 without hunter"
    push!(tests, Test(id,
                      splitplotanova(@where(data, :trial .== 3, :id .!= 3),
                                     exp.tests[id].locals,
                                     betweenfactors = [],
                                     withinfactors = [:tray, :foodtype],
                                     id = :id, data = :items),
                      exp.tests))
    TestCollection(tests)
end
function run!(exp::Experiment{:Cheke11_planning}, models; N = 4)
    results = DataFrame(trial = Int[], tray = String[], action = String[],
                        foodtype = String[],
						id = Int[], items = Float64[])
	for id in 1:N
		m = models[id]
		add!(m, MaintenanceDiet)
		for trial in 1:3
			# stage 1
		    trays = [Tray(traypos) for traypos in 1:2]
			remove!(m, MaintenanceDiet)
			wait!(m, 1u"hr")
			if trial == 1
				add!(m, MaintenanceDiet)
			elseif id == 2
				add!(m, PowderedSuetPellet, 20)
			else
				add!(m, PowderedPeanut, 20)
			end
            samefirst = true
			if id <= 2 # Hoy and Wiggins
				samefirst = false
			end
			foodtypeorder = [Peanut, SuetPellet]
			if id == 2
				reverse!(foodtypeorder)
			end
			wait!(m, 15u"minute")
            remove!(m, Any)
			add!(m, Peanut, 40)
			add!(m, SuetPellet, 40)
			add!(m, trays)
			wait!(m, 15u"minute")
			remove!(m, Any)
			add!(m, MaintenanceDiet)
			for traypos in 1:2
				for foodid in 1:2
					push!(results, [trial,
                                    ((samefirst && traypos == 1) || (!samefirst && traypos == 2)) ? "same" : "different",
                                    "cache",
                                    foodid == 1 ? "prefed" : "nonprefed",
									id, countcache(trays[traypos],
													foodtypeorder[foodid])])
				end
			end
			wait!(m, 3u"hr")
			# stage 2 & 3
			for stage in 2:3
				remove!(m, MaintenanceDiet)
				wait!(m, 1u"hr")
				if id == 1 && stage == 2 || id != 1 && stage == 3
					add!(m, PowderedSuetPellet, 20)
				else
					add!(m, PowderedPeanut, 20)
				end
				wait!(m, 15u"minute")
                remove!(m, Any)
                # cover tray 2 in stage 2 and tray 1 in stage 3
				cover!(trays[stage == 2 ? 2 : 1])
				add!(m, trays)
                # retrieve from tray 1 in stage 2 and tray 2 in stage 3
				wait!(m, 15u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				uncover!(trays[stage == 2 ? 2 : 1])
				wait!(m, 22.5u"hr")
			end
            pilfer!.(trays)
			wait!(m, 19.5u"hr")
		end
	end
    results
end


text:
using Test

include("util.jl")

using Random
Random.seed!(1)

using CuArrays

using GPUArrays
import GPUArrays: allowscalar, @allowscalar

testf(f, xs...; kwargs...) = GPUArrays.TestSuite.compare(f, CuArray, xs...; kwargs...)

allowscalar(false)

@testset "CuArrays" begin

include("base.jl")
include("dnn.jl")
include("blas.jl")
include("sparse.jl")
include("solver.jl")
include("fft.jl")
include("rand.jl")
include("sparse_solver.jl")

end


text:
using RMInsurance
using Test
using DataFrames
using Random
using LinearAlgebra
import LinearAlgebra.⋅

println("Testing ChainLadder")
include("ChainLadder_Test.jl")
println("Testing S2 Life")
include("S2Life_Test.jl")
println("Testing SST Life")
include("SSTLife_Test.jl")
println("Testing S2 Non-Life")
include("S2NonLife_Test.jl")
println("Testing ECModel")
include("ECModel_Test.jl")

println("tests completed")


text:
using ConvDiffMIPDECO
using jInv.Mesh
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.LinearSolvers
using MUMPSjInv
using MAT
using LinearAlgebra

# read results from l-curve
dataset    = "Sources3D"
reg        = "wTVReg"
m          = [96 48 48]
noiseLevel = 0.1

filename= joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples","$(dataset).mat")
data = matread(filename)
domain = data["domain"]
v      = data["v"]
sig    = data["sig"]
dtrue  = data["dtrue"]
rec    = data["rec"];
mf     = round.(Int64,data["m"])

M = getRegularMesh(domain,m)

resFile = joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples",
                 "$(dataset)-$(reg)-noise-$(noiseLevel)-$(m[1])x$(m[2])x$(m[3]).mat")
res = matread(resFile)
dobs = res["dobs"]
srcRelaxed = res["SourcesRelaxed"]
alphaRelaxed = res["alphaRelaxed"]

resFilePenalty = joinpath(dirname(pathof(ConvDiffMIPDECO)),"..","examples",
                 "$(dataset)-$(reg)-noise-$(noiseLevel)-$(m[1])x$(m[2])x$(m[3])-penalty.mat")

println("solve using alpha = $(alphaRelaxed)")
noiseLevel = norm(dobs-dtrue)/norm(dtrue)
println("estimated noise level: $noiseLevel")
if noiseLevel > 0.5
    error("noise level too large, double check that data was loaded correctly. ")
end

# build inverse problem
Mfine = getRegularMesh(domain,vec(data["m"]));
x1c,x2c,x3c = getCellCenteredAxes(Mfine)
rec3D = [kron(ones(length(x3c)),rec) kron(x3c,ones(size(rec,1)))]
# build linear interpolation matrix from nodes to receiver locations
x1,x2,x3 = getNodalAxes(M)
P = interpmat(x1,x2,x3, rec3D);
pFor = getConvDiffFEMParam(M,v=v,sig=sig,P=P,Ainv=getMUMPSsolver());
## configure misfit
Wt         = ones(size(dobs))/sqrt(mf[3])
sigback    = 0.0
pMis       = getMisfitParam(pFor,Wt,dobs,SSDFun)

## Configure regularization
mref       = zeros(M.nc,2)
regTV      = (m,mr,M,I=1.0) -> wTVReg(m,mr,M,eps=1e-8)
regP       = (m,mr,M,I=1.0) -> mipReg(m,mr,M)
reg        = [regTV;regP]

## Configure optimization
maxIter    = 50
minUpdate  = 1e-3
HesPrec    = getSSORRegularizationPreconditioner(1.0,1e-15,50)
cgit       = 5
pcgTol     = 1e-1
modFun     = identityMod
boundsLow  = 0*ones(M.nc)
boundsHigh = 1*ones(M.nc)
maxStep	   = 0.1*maximum(boundsHigh)

## store the configuration
pInv       = getInverseParam(M,modFun,reg,[alphaRelaxed;0.0],mref,
                             boundsLow,boundsHigh,maxStep=maxStep,
                            pcgMaxIter=cgit,pcgTol=pcgTol,minUpdate=minUpdate,maxIter=maxIter,
                            HesPrec=HesPrec);

mc = mref[:,1].+0.1
# viewImage2D(mc,M)
println("alpha:\t\t\t\t$(pInv.alpha)")
println("no. cells:\t\t\t$(length(mc))")
println("no. nodes:\t\t\t$(size(P,2))")
println("no. measurements:\t$(size(P,1))")
runtimeRelaxed = @elapsed begin
    mc,Dc,flag,HisRelaxed = projGNCG(mc,pInv,pMis)
    uc = pFor.Fields;          # u variables
    wr = pInv.modelfun(mc)[1]; # w variables
end

println("\n\n--Penalty method--\n")
Imax = 30
epsilon = 1e-4
betaMin = 1e-6
pInv.alpha[2] = betaMin
pInv.maxIter = 10
Ms = zeros(length(mc),Imax+1)
iter=1;
Ms[:,1] = mc;
HIS = Array{Any}(undef,Imax)
runtimePenalty = @elapsed begin
for iter=1:Imax
    mround,Dc,flag,HIS[iter] = projGNCG(Ms[:,iter],pInv,pMis)
    uc = pFor.Fields;          # u variables
    wr = pInv.modelfun(mround)[1]; # w variables

    Ms[:,iter+1] = mround;
    pInv.alpha[2] *= 2
end
end
println("runtime:\t$runtimePenalty")
println("no. PDE solves:\t$(pFor.Ainv.nSolve)")
println("time PDE factor:\t$(pFor.Ainv.facTime)")
println("time PDE solves:\t$(pFor.Ainv.solveTime)")

res["HIS"] = HIS
res["Ms"] = Ms
res["runtimeRelaxed"] = runtimeRelaxed
res["runtimePenalty"] = runtimeRelaxed
res["alphaPenalty"] = pInv.alpha
res["HisRelaxed"] = HisRelaxed

matwrite(resFilePenalty,res)

max_stars_repo_name: JuliaInv/ConvDiffMIPDECO

max_issues_repo_name: JuliaInv/ConvDiffMIPDECO

max_forks_repo_name: JuliaInv/ConvDiffMIPDECO


text:
const SimpleGraphEdge = SimpleEdge

"""
    SimpleGraph{T}

A type representing an undirected graph.
"""
mutable struct SimpleGraph{T<:Integer} <: AbstractSimpleGraph
    ne::Int
    fadjlist::Vector{Vector{T}} # [src]: (dst, dst, dst)
end

eltype(x::SimpleGraph{T}) where T<:Integer = T

# Graph{UInt8}(6), Graph{Int16}(7), Graph{UInt8}()
function (::Type{SimpleGraph{T}})(n::Integer = 0) where T<:Integer
    fadjlist = Vector{Vector{T}}()
    sizehint!(fadjlist, n)
    for _ = one(T):n
        push!(fadjlist, Vector{T}())
    end
    vertices = one(T):T(n)
    return SimpleGraph{T}(0, fadjlist)
end

# Graph()
SimpleGraph() = SimpleGraph{Int}()

# Graph(6), Graph(0x5)
SimpleGraph(n::T) where T<:Integer = SimpleGraph{T}(n)

# SimpleGraph(UInt8)
SimpleGraph(::Type{T}) where T<:Integer = SimpleGraph{T}(zero(T))

# Graph{UInt8}(adjmx)
function (::Type{SimpleGraph{T}})(adjmx::AbstractMatrix) where T<:Integer
    dima,dimb = size(adjmx)
    isequal(dima,dimb) || error("Adjacency / distance matrices must be square")
    issymmetric(adjmx) || error("Adjacency / distance matrices must be symmetric")

    g = SimpleGraph(T(dima))
    for i in find(triu(adjmx))
        ind = ind2sub((dima,dimb),i)
        add_edge!(g,ind...)
    end
    return g
end

# converts Graph{Int} to Graph{Int32}
function (::Type{SimpleGraph{T}})(g::SimpleGraph) where T<:Integer
    h_fadj = [Vector{T}(x) for x in fadj(g)]
    return SimpleGraph(ne(g), h_fadj)
end

# Graph(adjmx)
SimpleGraph(adjmx::AbstractMatrix) = SimpleGraph{Int}(adjmx)

# Graph(digraph)
function SimpleGraph(g::SimpleDiGraph)
    gnv = nv(g)
    edgect = 0
    newfadj = deepcopy(g.fadjlist)
    for i in vertices(g)
        for j in badj(g,i)
            if (_insert_and_dedup!(newfadj[i], j))
                edgect += 2     # this is a new edge only in badjlist
            else
                edgect += 1     # this is an existing edge - we already have it
                if i == j
                    edgect += 1 # need to count self loops
                end
            end
        end
    end
    iseven(edgect) || throw(AssertionError("invalid edgect in graph creation - please file bug report"))
    return SimpleGraph(edgect ÷ 2, newfadj)
end

edgetype(::SimpleGraph{T}) where T<:Integer = SimpleGraphEdge{T}

"""
    badj(g::SimpleGraph[, v::Integer])

Return the backwards adjacency list of a graph. If `v` is specified,
return only the adjacency list for that vertex.

###Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
badj(g::SimpleGraph) = fadj(g)
badj(g::SimpleGraph, v::Integer) = fadj(g, v)

"""
    adj(g[, v])

Return the adjacency list of a graph. If `v` is specified, return only the
adjacency list for that vertex.

### Implementation Notes
Returns a reference, not a copy. Do not modify result.
"""
adj(g::SimpleGraph) = fadj(g)
adj(g::SimpleGraph, v::Integer) = fadj(g, v)

copy(g::SimpleGraph) =  SimpleGraph(g.ne, deepcopy(g.fadjlist))

==(g::SimpleGraph, h::SimpleGraph) =
vertices(g) == vertices(h) &&
ne(g) == ne(h) &&
fadj(g) == fadj(h)

"""
    is_directed(g)

Return `true` if `g` is a directed graph.
"""
is_directed(::Type{SimpleGraph}) = false
is_directed(::Type{SimpleGraph{T}}) where T = false
is_directed(g::SimpleGraph) = false

function has_edge(g::SimpleGraph, e::SimpleGraphEdge)
    u, v = Tuple(e)
    u > nv(g) || v > nv(g) && return false
    if degree(g,u) > degree(g,v)
        u, v = v, u
    end
    return length(searchsorted(fadj(g,u), v)) > 0
end

function add_edge!(g::SimpleGraph, e::SimpleGraphEdge)
    T = eltype(g)
    s, d = T.(Tuple(e))
    (s in vertices(g) && d in vertices(g)) || return false
    inserted = _insert_and_dedup!(g.fadjlist[s], d)
    if inserted
        g.ne += 1
    end
    if s != d
        inserted = _insert_and_dedup!(g.fadjlist[d], s)
    end
    return inserted
end

function rem_edge!(g::SimpleGraph, e::SimpleGraphEdge)
    i = searchsorted(g.fadjlist[src(e)], dst(e))
    length(i) > 0 || return false   # edge not in graph
    i = i[1]
    deleteat!(g.fadjlist[src(e)], i)
    if src(e) != dst(e)     # not a self loop
        i = searchsorted(g.fadjlist[dst(e)], src(e))[1]
        deleteat!(g.fadjlist[dst(e)], i)
    end
    g.ne -= 1
    return true # edge successfully removed
end

"""
    add_vertex!(g)

Add a new vertex to the graph `g`. Return `true` if addition was successful.
"""
function add_vertex!(g::SimpleGraph)
    T = eltype(g)
    (nv(g) + one(T) <= nv(g)) && return false       # test for overflow
    push!(g.fadjlist, Vector{T}())
    return true
end

max_stars_repo_name: AsileBcd/JohnsonCycles

max_issues_repo_name: AsileBcd/JohnsonCycles

max_forks_repo_name: AsileBcd/JohnsonCycles