text:
# This file contains funtions for analysing results from chains.

"""
    projectpositions(chains.planets[1], mjd("2020-02-02"))

Given the posterior for a particular planet in the model and a modified julian date(s),
return `ra` and `dec` offsets in mas for each sampling in the posterior.
"""
function projectpositions(chains, planet_key, times)

    ras = zeros(size(chains,1),size(chains,3),length(times))
    decs = zeros(size(chains,1),size(chains,3),length(times))

    for is in collect(Iterators.partition(1:size(chains,1), 5000))
        for j in 1:size(chains,3)
            els = DirectDetections.construct_elements(chains, planet_key, is .* j)
            for (i,el) in zip(is,els)
                for (k, t) in enumerate(times)
                    ra, dec, _ = kep2cart(el, t)
                    ras[i,j,k] = ra
                    decs[i,j,k] = dec
                end
            end
        end
    end
    return ras, decs
end
export projectpositions

# """
#     sampleorbits(chains, planet_key, 100)

# Given the posterior for a particular planet in the model, and the number of orbits to sample,
# return a random subset of length N.
# """
# function sampleorbits(chains, planetnum, N)
#     return map(rand(eachindex(chains.planets[planetnum].a), N)) do i
#         return KeplerianElements(;
#             a = chains.planets[planetnum].a[i],
#             i = chains.planets[planetnum].i[i],
#             e = chains.planets[planetnum].e[i],
#             ω = chains.planets[planetnum].ω[i],
#             Ω = chains.planets[planetnum].Ω[i],
#             μ = chains.μ[i],
#             plx = chains.plx[i],
#             τ = chains.planets[planetnum].τ[i],
#         )
#     end
# end
# export sampleorbits

@recipe function f(astrom::Astrometry) where T
    xerror := astrom.σ_ra
    yerror := astrom.σ_dec
    xflip --> true
    xguide --> "Δ right ascension (mas)"
    yguide --> "Δ declination (mas)"

    return astrom.ra, astrom.dec
end

function plotposterior end
export plotposterior
function plotposterior! end
export plotposterior!
function plotmodel end
export plotmodel
function plotmodel! end
export plotmodel!

# Optionally depend on Plots. If the user imports it, this code will be run to set up
# our `imshow` function.
using Requires
function init_plots()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin

        function plotposterior(args...;kwargs...)
            Plots.plot()
            plotposterior!(args...;kwargs...)
        end
        function plotposterior!(
            chain,
            planet_key,
            keys::Nothing,
            N=1500;
            alpha=0.02,
            lw=0.3,
            color=1,
            kwargs...
        )
            ii = rand(1:size(chain,1)*size(chain,3), N)
            elements = construct_elements(chain, planet_key, ii)

            Plots.plot!(;
                size=(700,700),
                dpi=200,
                fmt=:png,
                framestyle=:box,
                grid=:none,
                minorticks=true,
                aspectratio=1,
                fontfamily="Arial",
                margin=8Plots.mm,
                kwargs...
            )
            Plots.xlabel!("Δ right ascension (mas)")
            Plots.ylabel!("Δ declination (mas)")

            for i in eachindex(elements)
                Plots.plot!(elements[i], label="",color=color, lw=lw, alpha=alpha,)
            end
            Plots.scatter!([0],[0], marker=(:star,:black, 5), markerstrokewidth=1, markerstrokecolor=:white, label="")
        end
        function plotposterior!(
            chain,
            planet_key,
            property,
            N=1500;
            alpha=0.02,
            cmap=:turbo,
            rev=true,
            colorbartitle= property isa Symbol ? "$planet_key[$property]" : "",#"semi-major axis (au)",
            clims=nothing,
            lw=0.3,
            kwargs...
        )
            ii = rand(1:size(chain,1)*size(chain,3), N)
            elements = construct_elements(chain, planet_key, ii)

            Plots.plot!(;
                size=(550,550),
                dpi=200,
                fmt=:png,
                framestyle=:box,
                grid=:none,
                minorticks=true,
                aspectratio=1,
                fontfamily="Arial",
                margin=8Plots.mm,
                kwargs...
            )
            Plots.xlabel!("Δ right ascension (mas)")
            Plots.ylabel!("Δ declination (mas)")

            if property isa Symbol
                k = string(planet_key)*"[$property]"
                colours = chain[k][ii]
            else # Assume they passed in the values directly
                colours = property[ii]
            end

            if isnothing(clims)
                clims = extrema(colours)
            end

            minc = first(clims)
            dc = last(clims) - minc
            clims = (minc, dc+minc)

            cmap= Plots.cgrad(cmap; rev)
            for i in eachindex(colours)#sortperm(colours, rev=false)
                c = colours[i]
                Plots.plot!(elements[i]; label="",color=cmap[(c-minc)/dc], lw, alpha)
            end
            Plots.scatter!([0],[0], marker=(:star,:black, 5), markerstrokewidth=1, markerstrokecolor=:white,  label="")
            Plots.scatter!([0],[0]; marker_z=[0], ms=0, clims=clims, color=cmap, label="", colorbartitle)

        end

        function plotmodel(args...; kwargs...)
            Plots.plot()
            plotmodel!(args...;kwargs...)
        end
        function plotmodel!(images::Images,i; kwargs...)
            img = DirectImages.DirectImage(images.image[i])
            img.PLATESCALE = images.platescale[i]
            imshow!(img;skyconvention=true,kwargs...)
            # plot!(sampled,color=:white,alpha=5/N,label="")
            Plots.xlims!(img.PLATESCALE.*extrema(axes(img,1)))
            Plots.ylims!(img.PLATESCALE.*extrema(axes(img,2)))

        end

        function plotmodel!(
            chain,
            N=1500,
            system=chain.info.model;
            alpha=(N <= 15 ? 1 : 30/N),
            color= isnothing(system.propermotionanom) ? :a : :mass,
            plotpma=!isnothing(system.propermotionanom),
            # TODO: ideally this is based on if there is a mass variable
            plotmass=!isnothing(system.propermotionanom),
            cmap=:plasma,
            imagecmap=:Greys,
            pma_scatter=nothing,
            clims=quantile(Iterators.flatten(
                # extrema(planet.a)
                # extrema([sample.planets[key].a for sample in chain])
                # for key in keys(chain[1].planets)
                extrema(chain["$pk[$color]"] for pk in keys(system.planets))
            ), (0.01,0.99)),
            lims=nothing,
            kwargs...
        )

            # Plot orbits
            p_orbits = Plots.plot(
                xlims=:symmetric,
                ylims=:symmetric,
            )

            # plot images?
            if !isnothing(system.images)
                if length(unique(system.images.platescale)) != 1
                    @warn "Plotmodel does not yet support images with multiple platescales. Taking first image only."
                    img = DirectImages.DirectImage(first(system.images.image))
                else
                    img = DirectImages.DirectImage(
                        DirectImages.stack(maximum,system.images.image)
                    )
                end
                img.PLATESCALE = system.images.platescale[1]
                # To get the colour scales to work out
                img ./= only(quantile(filter(isfinite, arraydata(img)), [0.98]))#[0.9995]))
                img .*= maximum(clims) - minimum(clims)
                img .+= minimum(clims)
                imshow!(img; color=imagecmap, skyconvention=true, lims)
            end

            for planet_key in eachindex(system.planets)
                plotposterior!(
                    chain, planet_key, color, N; lw=1, alpha=alpha,
                    cmap=cmap, rev=false,
                    # cmap=:turbo,
                    clims=clims,
                    kwargs...
                )

                # Planet astrometry?
                astrom = DirectDetections.astrometry(system.planets[planet_key])
                if !isnothing(astrom)
                    # Put black error bars over thicker white error bars to make them easier to see
                    # when plotted over e.g. images, orbits.
                    Plots.scatter!(astrom,marker=(:white,:circle,0),label="",linewidth=0,color=:white,markerstrokewidth=3, markerstrokecolor=:white)
                    Plots.scatter!(astrom,marker=(:black,:circle,0),label="",linewidth=0,color=:black,markerstrokewidth=1, markerstrokecolor=:black)
                end
            end
            # We will override this if we have more information available further down
            final_plot = p_orbits

            # astrometric acceleration?
            if plotpma

                lpma =length(system.propermotionanom.ra_epoch)
                pma_ii = sortperm(system.propermotionanom.ra_epoch)
                if lpma == 2 &&
                    isapprox(minimum(system.propermotionanom.ra_epoch), mjd("1991"), atol=365) &&
                    isapprox(maximum(system.propermotionanom.ra_epoch), mjd("2016"), atol=365)

                    titles=["Hipparcos", "GAIA EDR3"]
                else
                    titles = string.(system.propermotionanom.ra_epoch)
                end
                system_pma = system.propermotionanom

                pma_extrema_x = [-0.0,0.0]
                pma_extrema_y = [-0.0,0.0]
                pma_plots = map(sortperm(system_pma.ra_epoch)) do i
                    vx = zeros(prod(size(chain,[1,3])))
                    vy = zeros(prod(size(chain,[1,3])))
                    for j in keys(system.planets)
                        elements = construct_elements(chain, j, 1:prod(size(chain,[1,3])))
                        # mass = [sample.planets[j].mass for sample in chain]
                        mass = reshape(chain["$j[mass]"],:)
                        vx .+= getindex.(propmotionanom.(elements, system_pma.ra_epoch[i], mass.*mjup2msol),1)
                        vy .+= getindex.(propmotionanom.(elements, system_pma.dec_epoch[i], mass.*mjup2msol),2)
                    end
                    Plots.plot(
                        framestyle=:box,
                        minorticks=true,
                        aspectratio=1,
                        grid=false,
                        xlims=:symmetric,
                        ylims=:symmetric,
                    )
                    if !isnothing(pma_scatter)
                        if length(system.planets) > 1
                            @warn "Multiplanet PMA scatter plots not yet implemented"
                        end
                        ii = rand(1:size(chain,1), N)
                        planet_key = first(system.planets).name
                        prop = chain["$planet_key[$pma_scatter]"][ii]
                        # if pma_scatter == :mass
                        #     prop ./=mjup2msol
                        # end
                        if color == pma_scatter
                            clims_pma = clims
                        else
                            clims_pma = quantile(prop, (0.01, 0.99))
                        end
                        Plots.scatter!(vx[ii], vy[ii], marker_z=prop, alpha=1, color=:plasma, legend=false, colorbar=false, label="", markerstrokewidth=0, ms=1, clims=clims_pma)
                    else
                        h = fit(Histogram, (vx, vy))#, (-1:0.05:1, -1:0.05:1))
                        Plots.plot!(h, color=Plots.cgrad([Plots.RGBA(0,0,0,0), Plots.RGBA(0,0,0,1)]), colorbar=false)
                    end
                    pma_plot = Plots.scatter!(
                        [system_pma.pm_ra[i]],
                        [system_pma.pm_dec[i]],
                        xerror=[system_pma.σ_pm_ra[i]],
                        yerror=[system_pma.σ_pm_dec[i]],
                        label="",
                        color=:red,
                        markercolor=:red,
                        markerstrokecolor=:red,
                    )
                    # Plots.scatter!([0], [0], marker=(5, :circle, :red),label="")
                    Plots.hline!(pma_plot, [0], color=:black, label="")
                    Plots.vline!(pma_plot, [0], color=:black, label="")
                    Plots.title!(pma_plot, titles[i])
                    # Plots.xlims!(-1,1)
                    # Plots.ylims!(-1,1)
                    Plots.xlabel!(pma_plot, "Δμ ra - mas/yr")
                    Plots.ylabel!(pma_plot, "Δμ dec - mas/yr")
                    pma_extrema_x[1] = min(pma_extrema_x[1], minimum(vx), system_pma.pm_ra[i] - system_pma.σ_pm_ra[i], system_pma.pm_ra[i] + system_pma.σ_pm_ra[i])
                    pma_extrema_x[2] = max(pma_extrema_x[2], maximum(vx), system_pma.pm_ra[i] - system_pma.σ_pm_ra[i], system_pma.pm_ra[i] + system_pma.σ_pm_ra[i])
                    pma_extrema_y[1] = min(pma_extrema_y[1], minimum(vy), system_pma.pm_dec[i] - system_pma.σ_pm_dec[i], system_pma.pm_dec[i] + system_pma.σ_pm_dec[i])
                    pma_extrema_y[2] = max(pma_extrema_y[2], maximum(vy), system_pma.pm_dec[i] - system_pma.σ_pm_dec[i], system_pma.pm_dec[i] + system_pma.σ_pm_dec[i])
                    return pma_plot
                end
                for plot in pma_plots
                    # Plots.xlims!(plot, (pma_extrema_x...,))
                    # Plots.ylims!(plot, (pma_extrema_y...,))
                    lim = maximum(abs, ([pma_extrema_x; pma_extrema_y]))
                    Plots.xlims!(plot, -lim, lim)
                    Plots.ylims!(plot, -lim, lim)
                end
                # pma_plot = plot(pma_plots..., size=(700,300),margin=5Plots.mm, guidefontsize=9, titlefontsize=9, top_margin=0Plots.mm)
                # pma_plot = Plots.plot(pma_plots, link=:both)

                # Crete new final plot
                l = eval(:(Plots.@layout [
                    A{0.65h}
                    Plots.grid(1, $lpma)
                ]))
                final_plot = Plots.plot(p_orbits, pma_plots..., layout=l, size=(550,750),
                    margin=5Plots.mm, guidefontsize=9, titlefontsize=9, top_margin=0Plots.mm)
            end

            if plotmass
                p = Plots.plot(;xlabel="", linewidth=3, label="")
                for p in keys(system.planets)
                    m = vec(chain["$p[mass]"])
                    h = fit(Histogram, m, nbins=round(Int,sqrt(length(m))))#, nbins=100)
                    Plots.plot!(h.edges[1][begin:end-1],seriestype=:step, h.weights, label="$p[mass]")
                end

                layout = eval(:(Plots.@layout [
                    A{0.8h}
                    B
                ]))
                final_plot = Plots.plot(final_plot, p, layout=layout, size=(550,850))
            end

            return final_plot
        end
    end
end


text:
# важный костыль - кто угадает зачем?
# подсказка: туториал по generated function
foo(::Type{Type{T}}) where {T} = T

## функция для чтения произвольного тупла типов из бинарного потока

@generated function read_t(io::IOBuffer, types::Type...)
    Expr(:tuple, [Expr(:call, :read, :io, foo(i)) for i in types]...)
end
# альтернативная форма
@generated function read_t(io::IOBuffer, types...)
  args = Expr[]
  for t in types
    ex = :(read(io, (foo($t))))
    push!(args, ex)
  end
  return :(tuple($(args...)))
end

## чтение из большого бинарного массива - в вектор туплов
function read_tvector(data::Vector{UInt8}, types::Type...)
    elsize = mapreduce(sizeof, +, types)
    len = sizeof(data) ÷ elsize
    result = Vector{Tuple{types...}}(undef, len)
    io = IOBuffer(data)
    @inbounds for i in 1:len
        result[i] = read_t(io, types...)
    end
    result
end
## генерит тупл векторов заданной длины из тупла типов
@generated function tvectors(len, types...)
  Expr(:tuple, [Expr(:call, Expr(:curly, :Vector, foo(T)), :undef, :len) for T in types]...)
end
# альтернативная форма
@generated function tvectors(len, types...)
  args = Expr[]
  for T in types
    ex = :(Vector{foo($T)}(undef, len))
    push!(args, ex)
  end
  return :(tuple($(args...)))
end

## чтение из большого бинарного массива - в тупл векторов
function read_tvectors(data::Vector{UInt8}, types::Type...)
    elsize = mapreduce(sizeof, +, types)
    len = sizeof(data) ÷ elsize
    vectors = tvectors(len, types...)
    io = IOBuffer(data)
    eltypes = length(types)
    @inbounds for i in 1:len
        for j = 1:eltypes
            vectors[j][i] = read(io, types[j])
        end
    end
    vectors
end

##
bytevec = UInt8[0x01, 0x02, 0x00, 0x03, 0x00, 0x00, 0x00]
types = (Int8, Int16, Int32)
io = IOBuffer(bytevec)
x = read_t(io, types...)

types = (Int8, Int16, Int32)
elsize = mapreduce(sizeof, +, types)
bytevec = rand(UInt8, 1000 * elsize)
vec = read_tvector(bytevec, types...)

types = (Int8, Int16, Int32)
vectors = tvectors(1000, types...)

elsize = mapreduce(sizeof, +, types)
bytevec = rand(UInt8, 1000 * elsize)
vec = read_tvectors(bytevec,types...)

ext: jl


text:
export
        LaneChangeAction,
        LaneChangeChoice,
        LaneChangeModel,
        get_lane_offset,

        DIR_RIGHT,
        DIR_MIDDLE,
        DIR_LEFT,

        TimLaneChanger

abstract type LaneChangeAction <: DriveAction end

"""
    LaneChangeChoice
A choice of whether to change lanes, and what direction to do it in
"""
const DIR_RIGHT = -1
const DIR_MIDDLE =  0
const DIR_LEFT =  1
immutable LaneChangeChoice <: LaneChangeAction
    dir::Int # -1, 0, 1
end
Base.show(io::IO, a::LaneChangeChoice) = @printf(io, "LaneChangeChoice(%d)", dir)
Base.length(::Type{LaneChangeChoice}) = 1
Base.convert(::Type{LaneChangeChoice}, v::Vector{Float64}) = LaneChangeChoice(convert(Int, v[1]))
function Base.copy!(v::Vector{Float64}, a::LaneChangeChoice)
    v[1] = a.dir
    v
end
function get_lane_offset(a::LaneChangeChoice, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    if a.dir == DIR_MIDDLE
        rec[vehicle_index, pastframe].state.posF.t
    elseif a.dir == DIR_LEFT
        convert(Float64, get(LANEOFFSETLEFT, rec, roadway, vehicle_index, pastframe))
    else
        @assert(a.dir == DIR_RIGHT)
        convert(Float64, get(LANEOFFSETRIGHT, rec, roadway, vehicle_index, pastframe))
    end
end

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

abstract type LaneChangeModel{LaneChangeAction} end
get_name(::LaneChangeModel) = "???"
set_desired_speed!(::LaneChangeModel, v_des::Float64) = model # # do nothing by default
reset_hidden_state!(model::LaneChangeModel) = model # do nothing by default
observe!(model::LaneChangeModel, scene::Scene, roadway::Roadway, egoid::Int) = model  # do nothing by default
Base.rand(model::LaneChangeModel) = error("rand not implemented for model $model")

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

type TimLaneChanger <: LaneChangeModel{LaneChangeAction}
    dir::Int
    rec::SceneRecord

    v_des::Float64
    threshold_fore::Float64
    threshold_lane_change_gap_fore::Float64
    threshold_lane_change_gap_rear::Float64

    function TimLaneChanger(
        action_context::IntegratedContinuous;
        v_des::Float64=29.0,
        rec::SceneRecord=SceneRecord(2,action_context.Δt),
        threshold_fore::Float64 = 50.0,
        threshold_lane_change_gap_fore::Float64 = 10.0,
        threshold_lane_change_gap_rear::Float64 = 10.0,
        dir::Int=DIR_MIDDLE,
        )

        retval = new()
        retval.dir = dir
        retval.rec = rec
        retval.v_des = v_des
        retval.threshold_fore = threshold_fore
        retval.threshold_lane_change_gap_fore = threshold_lane_change_gap_fore
        retval.threshold_lane_change_gap_rear = threshold_lane_change_gap_rear
        retval
    end
end
get_name(::TimLaneChanger) = "TimLaneChanger"
function set_desired_speed!(model::TimLaneChanger, v_des::Float64)
    model.v_des = v_des
    model
end
function observe!(model::TimLaneChanger, scene::Scene, roadway::Roadway, egoid::Int)
    #debug = open("debug.log", "a"); println(debug, "observe!(model::TimLaneChanger, scene::Scene, roadway::Roadway, egoid::Int)");  close(debug)

    rec = model.rec
    update!(rec, scene)
    vehicle_index = get_index_of_first_vehicle_with_id(rec, egoid)

    veh_ego = scene[vehicle_index]
    v = veh_ego.state.v

    left_lane_exists = convert(Float64, get(N_LANE_LEFT, rec, roadway, vehicle_index)) > 0
    right_lane_exists = convert(Float64, get(N_LANE_RIGHT, rec, roadway, vehicle_index)) > 0
    fore_M = get_neighbor_fore_along_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointRear(), VehicleTargetPointFront(), max_distance_fore=model.threshold_fore)
    fore_L = get_neighbor_fore_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    fore_R = get_neighbor_fore_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointRear(), VehicleTargetPointRear(), VehicleTargetPointFront())
    rear_L = get_neighbor_rear_along_left_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())
    rear_R = get_neighbor_rear_along_right_lane(scene, vehicle_index, roadway, VehicleTargetPointFront(), VehicleTargetPointFront(), VehicleTargetPointRear())

    model.dir = DIR_MIDDLE
    if fore_M.Δs < model.threshold_fore # there is a lead vehicle
        veh_M = scene[fore_M.ind]
        speed_M = veh_M.state.v
        if speed_M ≤ min(model.v_des, v) # they are driving slower than we want

            speed_ahead = speed_M

            # consider changing to a different lane
            if right_lane_exists &&
               fore_R.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_R.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_R.ind == 0 || scene[rear_R.ind].state.v ≤ v) && # we are faster than any follower
               (fore_R.ind == 0 || scene[fore_R.ind].state.v > speed_ahead) # we are faster than any leader

                speed_ahead = fore_R.ind != 0 ? scene[fore_R.ind].state.v : Inf
                model.dir = DIR_RIGHT
            end
            if left_lane_exists &&
               fore_L.Δs > model.threshold_lane_change_gap_rear && # there is space rear
               rear_L.Δs > model.threshold_lane_change_gap_fore && # there is space fore
               (rear_L.ind == 0 || scene[rear_L.ind].state.v ≤ v) && # we are faster than any follower
               (fore_L.ind == 0 || scene[fore_L.ind].state.v > speed_ahead) # we are faster than any leader

                speed_ahead = fore_L.ind != 0 ? scene[fore_L.ind].state.v : Inf
                model.dir = DIR_LEFT
            end
        end
    end

    model
end
Base.rand(model::TimLaneChanger) = LaneChangeChoice(model.dir)


text:
type Vehiculo
    velocidad::Int8
    posicion::Int64
    tipo::Int8
    cambio::Int8
    num::Int64
end

type Carretera1D
	cuenta::Int64
	N::Int64
	carretera::Vector{Vehiculo}
end

function Carretera1D(N::Int64)
    carretera = Vehiculo[]
    for i in 1:N
    # Aquí no es necesario poner la función Celda_Vacia
        celda = Vehiculo(round(Int8, -1), i, round(Int8, -1), zero(Int8), 0)
        push!(carretera, celda)
    end
    return Carretera1D(0, N, carretera)
end

type Carretera2D
  carriles::Int8
  carretera::Vector{Carretera1D}
end

function Carretera2D(carriles::Int64, N::Int64)
    C = Carretera1D[]
    if carriles > 1
        for i = 1:carriles
            push!(C, Carretera1D(N))
        end
        return Carretera2D(round(Int8,carriles), C)
    end
    if carriles == 1
        return Carretera1D(N)
    end
end


text:
using RecipesBase

@recipe function f(s::T) where {T <: AbstractSimplex}
    seriestype := :path
    splitaxes(connectvertices(s))
end

@recipe function f(simplices::Vararg{T,N}) where {T <: AbstractSimplex, N}
    seriestype := :path
    legend --> false

    for simplex in simplices
        @series begin
            label --> ""
            splitaxes(connectvertices(simplex))
        end
    end

end

@recipe function f(simplices::AbstractVector{T}) where {T <: AbstractSimplex}
    seriestype := :path
    legend --> false

    for simplex in simplices
        @series begin
            label --> ""
            splitaxes(connectvertices(simplex))
        end
    end

end


text:
using RAFF
using Distributed

#
# This example is the "Parallel running" example from the
# documentation.
#

# Add 3 worker processes
addprocs(3)

# Distribute RAFF, the model and its gradien (the gradient is not
# mandatory, just an example)
@everywhere using RAFF

@everywhere function model(x, θ)
   θ[1] * x[1]^2 + θ[2]
end

@everywhere function gmodel!(g, x, θ)
    g[1] = x[1]^2
    g[2] = 1.0
end

# Define the data matrix
A=[-2.0  5.0;
  -1.5  3.25;
  -1.0  2.0 ;
  -0.5  1.25;
   0.0  1.0 ;
   0.5  2.55;
   1.0  2.0 ;
   1.5  3.25;
   2.0  5.0 ;];

# Number of parameters of the model
n = 2

# Run Parallel RAFF
output = praff(model, gmodel!, A, n)

# Print output
println(output)


text:
include("SourceParallelism.jl");
include("modelUtils.jl");
include("SrcRcvUtils.jl");
include("BoundModel.jl");


text:
@testset "B-splines" begin
    include("knot_sets.jl")
    include("splines.jl")
    include("operators.jl")
    include("derivatives.jl")
end


text:
using Base.Sys
"""
Convert DynamcHMC samples to a chain
* `posterior`: an array of NamedTuple consisting of mcmcm samples
"""
function nptochain(results, posterior)
    Np = length(vcat(posterior[1]...)) + 2 #include lf_eps
    Ns = length(posterior)
    a3d = Array{Float64, 3}(undef, Ns, Np, 1)
    depth = map(x->x.depth, results.tree_statistics)
    ϵ = results.ϵ
    i = 0
    for (post,ch) in zip(posterior, depth)
        i += 1
        temp = Float64[]
        for p in post
            push!(temp, values(p)...)
        end
        push!(temp, ϵ, ch)
        a3d[i,:,1] = temp'
    end
    parameter_names = getnames(posterior)
    push!(parameter_names, "lf_eps","tree_depth")
    chns = MCMCChains.Chains(a3d, parameter_names,
        Dict(:internals => ["lf_eps", "tree_depth"]); sorted=true)
    return chns
end

function getnames(post)
    nt = post[1]
    Np = length(vcat(nt...))
    parm_names = fill("", Np)
    cnt = 0
    for (k,v) in pairs(nt)
        N = length(v)
        if isa(v, Array)
            for i in 1:N
                cnt += 1
                parm_names[cnt] = string(k, "[",i,"]")
            end
        else
            cnt+=1
            parm_names[cnt] = string(k)
        end
    end
    return parm_names
end

function setprocs(n)
    np = nprocs()-1
    m = max(n-np, 0)
    addprocs(m)
end

function save(results,ProjDir)
    str = string(round(now(), Dates.Minute))
    str = replace(str, "-"=>"_")
    str = replace(str, ":"=>"_")
    dir = string(ProjDir, "/results")
    !isdir(dir) ? mkdir(dir) : nothing
    newdir = dir*"/"*str
    mkdir(newdir)
    CSV.write(newdir*"/results.csv", results)
    metadata = getMetadata()
    CSV.write(newdir*"/metadata.csv", metadata)
end

function getMetadata()
    path = getpath()
    dict = TOML.parsefile(path)
    df = DataFrame()
    pkgs = [:CmdStan, :DynamicHMC,
        :Turing, :AdvancedHMC]
    map(p->df[!,p]=[dict[string(p)][1]["version"]], pkgs)
    df[!,:julia] = [VERSION]
    df[!,:os] = [MACHINE]
    cpu = cpu_info()
    df[!,:cpu] = [cpu[1].model]
    return df
end

function getpath()
    userdir = expanduser("~")
    str = split(string(VERSION), ".")[1:end-1]
    str1 = [str[s]*"." for s in 1:length(str)-1]
    push!(str1, str[end])
    version = string(str1...)
    path = joinpath(userdir, string(".julia/environments/v", version,"/Manifest.toml"))
    return path
end

"""
Iterates over all permutations of factors in the benchmark simulation.
"""
mutable struct Permutation{T1,T2,T3,T4}
    a::T1
    names::T2
    idx::Vector{Int64}
    N::T3
    start::T4
end

function Permutation(kwargs)
    names = kwargs.itr
    a = values(kwargs.data)
    idx = fill(1,length(a))
    N = length.(a)
    v = [i[1] for i in a]
    start = NamedTuple{names}(v)
    return Permutation(a, names, idx, N, start)
end

Base.length(P::Permutation) = prod(P.N)

Base.eltype(::Type{Permutation{T}}) where {T} = Vector{eltype(T)}

Base.collect(P::Permutation) = [p for p in P]

function Base.iterate(p::Permutation,state=(p.start,0))
    @unpack a,idx,names,N=p
    element,count = state
    count == length(p) ? (return nothing) : nothing
    K = length(p.idx)
    for j in K:-1:1
        if idx[j] == N[j]
            idx[j] = 1
        else
            idx[j] += 1
            break
        end
    end
    count += 1
    v = [a[i][j] for (i,j) in enumerate(idx)]
    newVal = NamedTuple{names}(v)
    return (element, (newVal, count))
end

"""
remove compile time from benchmarks
"""
function compile(samplers, fun::F; kwargs...) where {F<:Function}
    p = collect(Permutation(kwargs))[1]
    data = fun(;Nd=1, p...); N=nprocs()
    pmap(s->compile(samplers, data; p...), 1:N)
end

function compile(samplers, data; kwargs...)
    for s in samplers
        println("Compiling $(typeof(s))")
        println()
        modifyConfig!(s; kwargs...)
        runSampler(s, data; kwargs...)
    end
end

function initStan(s)
    base = string(s.dir,"/tmp/",s.model.name)
    for p in procs()
        stream = open(base*".stan","r")
        str = read(stream, String)
        close(stream)
        stream = open(string(base, p,".stan"),"w")
        write(stream, str)
        close(stream)
    end
end

function setreps(Nreps)
    p = nprocs()-1
    v = Int(floor(Nreps/p))
    reps = fill(v,p)
    r = mod(Nreps,p)
    reps[1:r] .+= 1
    return reps
end

"""
change name parm.x to parm[x]
"""
function format_names(chain)
    parm_names = names(chain)
    new_names = map(format_name, parm_names)
    d = Dict(o=>n for (o,n) in zip(parm_names,new_names))
    chain = set_names(chain, d)
    return chain
end

function format_name(str)
    s = split(str, ".")
    if length(s) > 1
        return string(s[1], "[", s[2],"]")
    end
    return str
end


text:
"""
    traintensor([indices]; [dir], [decimal=true]) -> Array{Float64}

Returns the MNIST **training** images corresponding to the given
`indices` as a multi-dimensional array.

The corresponding source file of the dataset is expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

```julia
julia> MNIST.traintensor(dir="/home/user/MNIST")
WARNING: The MNIST file "train-images-idx3-ubyte.gz" was not found in "/home/user/MNIST". You can download [...]
```

The image(s) is/are returned in the native horizontal-major
memory layout as a single floating point array. If `decimal=true`
all values are scaled to be between `0.0` and `1.0`, otherwise
the values will be between `0.0` and `255.0`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 3D array (i.e. a `Array{Float64,3}`), in
which the first dimension corresponds to the pixel *rows* (x) of
the image, the second dimension to the pixel *columns* (y) of the
image, and the third dimension denotes the index of the image.

```julia
julia> MNIST.traintensor() # load all training images
28×28×60000 Array{Float64,3}:
[...]

julia> MNIST.traintensor(1:3) # load first three training images
28×28×3 Array{Float64,3}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Matrix{Float64}` in horizontal-major layout, which means that
the first dimension denotes the pixel *rows* (x), and the second
dimension denotes the pixel *columns* (y) of the image.

```julia
julia> MNIST.traintensor(1) # load first training image
28×28 Array{Float64,2}:
[...]
```

As mentioned above, the images are returned in the native
horizontal-major layout to preserve the original feature
ordering. You can use the utility function
[`convert2image`](@ref) to convert an MNIST array into a
vertical-major Julia image with the corrected color values.

```
julia> MNIST.convert2image(MNIST.traintensor(1)) # convert to column-major colorant array
28×28 Array{Gray{Float64},2}:
[...]
```
"""
function traintensor(args...; dir=DEFAULT_DIR, decimal=true)
    if decimal
        Reader.readtrainimages(dir, args...) ./ 255
    else
        convert(Array{Float64}, Reader.readtrainimages(dir, args...))
    end
end

"""
    testtensor([indices]; [dir], [decimal=true]) -> Array{Float64}

Returns the MNIST **test** images corresponding to the given
`indices` as a multi-dimensional array.

The corresponding source file of the dataset is expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

```julia
julia> MNIST.testtensor(dir="/home/user/MNIST")
WARNING: The MNIST file "t10k-images-idx3-ubyte.gz" was not found in "/home/user/MNIST". You can download [...]
```

The image(s) is/are returned in the native horizontal-major
memory layout as a single floating point array. If `decimal=true`
all values are scaled to be between `0.0` and `1.0`, otherwise
the values will be between `0.0` and `255.0`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 3D array (i.e. a `Array{Float64,3}`), in
which the first dimension corresponds to the pixel *rows* (x) of
the image, the second dimension to the pixel *columns* (y) of the
image, and the third dimension denotes the index of the image.

```julia
julia> MNIST.testtensor() # load all test images
28×28×10000 Array{Float64,3}:
[...]

julia> MNIST.testtensor(1:3) # load first three test images
28×28×3 Array{Float64,3}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Matrix{Float64}` in horizontal-major layout, which means that
the first dimension denotes the pixel *rows* (x), and the second
dimension denotes the pixel *columns* (y) of the image.

```julia
julia> MNIST.testtensor(1) # load first test image
28×28 Array{Float64,2}:
[...]
```

As mentioned above, the images are returned in the native
horizontal-major layout to preserve the original feature
ordering. You can use the utility function
[`convert2image`](@ref) to convert an MNIST array into a
vertical-major Julia image with the corrected color values.

```
julia> MNIST.convert2image(MNIST.testtensor(1)) # convert to column-major colorant array
28×28 Array{Gray{Float64},2}:
[...]
```
"""
function testtensor(args...; dir=DEFAULT_DIR, decimal=true)
    if decimal
        Reader.readtestimages(dir, args...) ./ 255
    else
        convert(Array{Float64}, Reader.readtestimages(dir, args...))
    end
end

"""
    trainlabels([indices]; [dir])

Returns the MNIST **trainset** labels corresponding to the given
`indices` as an `Int` or `Vector{Int}`. The values of the labels
denote the digit that they represent. If `indices` is omitted,
all labels are returned.

```julia
julia> MNIST.trainlabels() # full training set
60000-element Array{Int64,1}:
 5
 0
 ⋮
 6
 8

julia> MNIST.trainlabels(1:3) # first three labels
3-element Array{Int64,1}:
 5
 0
 4

julia> MNIST.trainlabels(1) # first label
5
```

The corresponding source file of the dataset is expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

```julia
julia> MNIST.trainlabels(dir="/home/user/MNIST")
WARNING: The MNIST file "train-labels-idx1-ubyte.gz" was not found in "/home/user/MNIST". You can download [...]
```
"""
trainlabels(args...; dir=DEFAULT_DIR) = Vector{Int}(Reader.readtrainlabels(dir, args...))
trainlabels(index::Integer; dir=DEFAULT_DIR) = Int(Reader.readtrainlabels(dir, index))

"""
    testlabels([indices]; [dir])

Returns the MNIST **testset** labels corresponding to the given
`indices` as an `Int` or `Vector{Int}`. The values of the labels
denote the digit that they represent. If `indices` is omitted,
all labels are returned.

```julia
julia> MNIST.testlabels() # full test set
10000-element Array{Int64,1}:
 7
 2
 ⋮
 5
 6

julia> MNIST.testlabels(1:3) # first three labels
3-element Array{Int64,1}:
 7
 2
 1

julia> MNIST.testlabels(1) # first label
7
```

The corresponding source file of the dataset is expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

```julia
julia> MNIST.testlabels(dir="/home/user/MNIST")
WARNING: The MNIST file "t10k-labels-idx1-ubyte.gz" was not found in "/home/user/MNIST". You can download [...]
```
"""
testlabels(args...; dir=DEFAULT_DIR) = Vector{Int}(Reader.readtestlabels(dir, args...))
testlabels(index::Integer; dir=DEFAULT_DIR) = Int(Reader.readtestlabels(dir, index))

"""
    traindata([indices]; [dir], [decimal=true]) -> Tuple

Returns the MNIST **trainingset** corresponding to the given
`indices` as a two-element tuple. If `indices` is omitted the
full trainingset is returned. The first element of thre return
value will be the images as a multi-dimensional array, and the
second element the corresponding labels as integers.

The images are returned in the native horizontal-major memory
layout as a single floating point array. If `decimal=true` all
values are scaled to be between `0.0` and `1.0`, otherwise the
values will be between `0.0` and `255.0`. The integer values of
the labels correspond 1-to-1 the digit that they represent.

```julia
train_x, train_y = MNIST.traindata() # full datatset
train_x, train_y = MNIST.traindata(2) # only second observation
train_x, train_y = MNIST.traindata(dir="./MNIST") # custom folder
```

The corresponding source files of the dataset are expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

Take a look at [`traintensor`](@ref) and [`trainlabels`](@ref)
for more information.
"""
function traindata(args...; dir=DEFAULT_DIR, decimal=true)
    (traintensor(args...; dir=dir, decimal=decimal),
     trainlabels(args...; dir=dir))
end

"""
    testdata([indices]; [dir], [decimal=true]) -> Tuple

Returns the MNIST **testset** corresponding to the given
`indices` as a two-element tuple. If `indices` is omitted the
full testset is returned. The first element of thre return value
will be the images as a multi-dimensional array, and the second
element the corresponding labels as integers.

The images are returned in the native horizontal-major memory
layout as a single floating point array. If `decimal=true` all
values are scaled to be between `0.0` and `1.0`, otherwise the
values will be between `0.0` and `255.0`. The integer values of
the labels correspond 1-to-1 the digit that they represent.

```julia
test_x, test_y = MNIST.testdata() # full datatset
test_x, test_y = MNIST.testdata(2) # only second observation
test_x, test_y = MNIST.testdata(dir="./MNIST") # custom folder
```

The corresponding source files of the dataset are expected to be
located in the specified directory `dir`. If `dir` is omitted it
will default to `MLDatasets/datasets/mnist`. In the case the
source files have not been downloaded yet, you can use
`MNIST.download_helper(dir)` to assist in the process.

Take a look at [`testtensor`](@ref) and [`testlabels`](@ref)
for more information.
"""
function testdata(args...; dir=DEFAULT_DIR, decimal=true)
    (testtensor(args...; dir=dir, decimal=decimal),
     testlabels(args...; dir=dir))
end


text:
module MetabolicSim

    import DiffEqBase, OrdinaryDiffEq, DiffEqCallbacks

    include("parameters/parameters.jl")
    include("simulation/assembly//assembly.jl")
    include("simulation/dxdt.jl")
    include("simulation/simulate.jl")
    include("analysis/analysis.jl")

    #parameter functions
    export make_parameters

    #callback functions
    # export add_at_equilibrium, add_at_t

    #simulate functions
    export dx!

    #analysis
    export get_timeseries, get_Nsp
    export get_resource_a, get_overlap_competition, get_overlap_facilitation
    export get_R,get_ω
end


text:
"""
```
sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == sum(pg[g] for g in bus_gens) + sum(pconvac[c] for c in bus_convs) - pd - gs*v^2
sum(q[a] for a in bus_arcs) + sum(q_dc[a_dc] for a_dc in bus_arcs_dc) == sum(qg[g] for g in bus_gens) + sum(qconvac[c] for c in bus_convs) - qd + bs*v^2
```
"""
function constraint_kcl_shunt(pm::AbstractACPModel, n::Int, cnd::Int, i::Int, bus_arcs, bus_arcs_dc, bus_gens, bus_convs_ac, bus_loads, bus_shunts, pd, qd, gs, bs)
    vm = PowerModels.var(pm, n, cnd, :vm, i)
    p = PowerModels.var(pm, n, cnd, :p)
    q = PowerModels.var(pm, n, cnd, :q)
    pg = PowerModels.var(pm, n, cnd, :pg)
    qg = PowerModels.var(pm, n, cnd, :qg)
    pconv_grid_ac = PowerModels.var(pm, n, cnd, :pconv_tf_fr)
    qconv_grid_ac = PowerModels.var(pm, n, cnd, :qconv_tf_fr)

    PowerModels.con(pm, n, cnd, :kcl_p)[i] = @NLconstraint(pm.model, sum(p[a] for a in bus_arcs) + sum(pconv_grid_ac[c] for c in bus_convs_ac)  == sum(pg[g] for g in bus_gens)   - sum(pd[d] for d in bus_loads) - sum(gs[s] for s in bus_shunts)*vm^2)
    PowerModels.con(pm, n, cnd, :kcl_q)[i] = @NLconstraint(pm.model, sum(q[a] for a in bus_arcs) + sum(qconv_grid_ac[c] for c in bus_convs_ac)  == sum(qg[g] for g in bus_gens)  - sum(qd[d] for d in bus_loads) + sum(bs[s] for s in bus_shunts)*vm^2)
end
"""
Creates Ohms constraints for DC branches

```
p[f_idx] == p * g[l] * vmdc[f_bus] * (vmdc[f_bus] - vmdc[t_bus])
```
"""
function constraint_ohms_dc_branch(pm::AbstractACPModel, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, r, p)
    p_dc_fr = PowerModels.var(pm, n, cnd, :p_dcgrid, f_idx)
    p_dc_to = PowerModels.var(pm, n, cnd, :p_dcgrid, t_idx)
    vmdc_fr = PowerModels.var(pm, n, cnd, :vdcm, f_bus)
    vmdc_to = PowerModels.var(pm, n, cnd, :vdcm, t_bus)

    if r == 0
        @constraint(pm.model, p_dc_fr + p_dc_to == 0)
    else
        g = 1 / r
        @NLconstraint(pm.model, p_dc_fr == p * g * vmdc_fr * (vmdc_fr - vmdc_to))
        @NLconstraint(pm.model, p_dc_to == p * g * vmdc_to * (vmdc_to - vmdc_fr))
    end
end
"`vdc[i] == vdcm`"
function constraint_dc_voltage_magnitude_setpoint(pm::AbstractACPModel, n::Int, cnd::Int, i, vdcm)
    v = PowerModels.var(pm, n, cnd, :vdcm, i)
    PowerModels.con(pm, n, cnd, :v_dc)[i] = @constraint(pm.model, v == vdcm)
end

function constraint_dc_branch_current(pm::AbstractACPModel, n::Int, cnd::Int, f_bus, f_idx, ccm_max, p)
# do nothing
end


text:
using IntegralArrays, LinearAlgebra, IntervalSets
using ColorTypes, ColorVectorSpace
using ColorTypes.FixedPointNumbers
using OffsetArrays
using Test
using Documenter

@testset "meta-quality" begin
    if VERSION >= v"1.6"
        doctest(IntegralArrays; manual = false)
    end
end

@testset "IntegralArrays" begin
    a = zeros(10, 10)
    int_array = IntegralArray(a)
    @test all(int_array == a)

    a = ones(10,10)
    int_array = IntegralArray(a)
    chk = Array(1:10)
    @test all([vec(int_array[i, :]) == chk * i for i in 1:10])

    int_sum = int_array[1..5, 1..2]
    @test int_sum == 10.0
    int_sum = int_array[3 ± 2, 1..2]
    @test int_sum == 10.0
    int_sum = int_array[3 ± 2, 3 ± 2]
    @test int_sum == 25.0
    int_sum = int_array[1..2, 1..5]
    @test int_sum == 10.0
    int_sum = int_array[4..8, 4..8]
    @test int_sum == 25.0
    int_sum = int_array[6 ± 2, 6 ± 2]
    @test int_sum == 25.0
    int_sum = int_array[6 ± 2, 6 ± 2]
    @test int_sum == 25.0

    a = Array(reshape(1:100, 10, 10))
    int_array = IntegralArray(a)
    @test diag(int_array) == Array([1, 26,  108,  280,  575, 1026, 1666, 2528, 3645, 5050])

    int_sum = int_array[1..3, 1..3]
    @test int_sum == 108
    int_sum = int_array[2 ± 1, 2 ± 1]
    @test int_sum == 108
    int_sum = int_array[2 ± 1, 2 ± 1]
    @test int_sum == 108
    int_sum = int_array[1..5, 1..2]
    @test int_sum == 80
    int_sum = int_array[3 ± 2, 1..2]
    @test int_sum == 80
    int_sum = int_array[4..8, 4..8]
    @test int_sum == 1400
    int_sum = int_array[6 ± 2, 6 ± 2]
    @test int_sum == 1400
    int_sum = int_array[6 ± 2, 6 ± 2]
    @test int_sum == 1400

    I, Δ = CartesianIndex(6, 6), CartesianIndex(2, 2)
    @test int_array[I-Δ..I+Δ] == 1400

    @testset "OffsetArray" begin
        X = rand(1:5, 5, 5)
        iX = IntegralArray(X)
        @test axes(iX) == (1:5, 1:5)
        Xo = OffsetArray(X, -1, -1)
        iXo = IntegralArray(Xo)
        @test axes(iXo) == (0:4, 0:4)

        @test_throws DimensionMismatch IntegralArray(X, Xo)
    end

    @testset "in-place with buffer" begin
        X = collect(reshape(1:25, 5, 5))
        iX = similar(X)
        IntegralArray(iX, X)
        IntegralArray(X, X)
        @test iX == X == IntegralArray(reshape(1:25, 5, 5))
    end
end

@testset "Color integral arrays" begin
    A = [RGB{N0f8}(0.2, 0.8, 0), RGB{N0f8}(0.5, 0.3, 0)]
    intA = IntegralArray(A)
    a1, a2 = A
    @test intA[2] == RGB(Float32(red(a1))+Float32(red(a2)),
                         Float32(green(a1))+Float32(green(a2)),
                         Float32(blue(a1))+Float32(blue(a2)))
    @test intA[1..2] == intA[2]
    @test intA[2..2] ≈ A[2]
end


text:
#=
reverser.jl
based on
https://stackoverflow.com/questions/27411401/julia-reverse-n-dimensional-arrays
=#

export reverser

using Test: @test

"""
    y = reverser(x, dims)
reverse array along specified dimensions (or all if unspecified)
"""
function reverser(x::AbstractArray, dims::AbstractVector{<:Int})
    y = copy(x)
    for d in dims
        y = reverse(y, dims=d)
    end
    return y
end

reverser(x::AbstractArray) = reverser(x, 1:ndims(x)) # all dimensions
reverser(x::AbstractArray, d::Int) = reverser(x, [d])

"""
    reverser(:test)
self test
"""
function reverser(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	@test reverser(1:3) == 3:-1:1
	@test reverser(1:3, 1) == 3:-1:1
	@test reverser((1:3)', 1) == (1:3)'
	@test reverser((1:3)', 2) == (3:-1:1)'
	true
end


text:
@testset "other-fhir-versions.jl" begin
    @test FHIRClient.DSTU2() isa FHIRClient.DSTU2
    @test FHIRClient.STU3() isa FHIRClient.STU3
    @test FHIRClient.UnknownFHIRVersion() isa FHIRClient.UnknownFHIRVersion
    @test FHIRClient.OtherFHIRVersion{:DSTU2}() isa FHIRClient.OtherFHIRVersion
    @test FHIRClient.OtherFHIRVersion{:DSTU2}() isa FHIRClient.OtherFHIRVersion{:DSTU2}
end


text:
import Pkg; Pkg.add("HDF5")
import Pkg; Pkg.add("ArgParse")
import Pkg; Pkg.add("Plots")
import Pkg; Pkg.add("LinearAlgebra")
import Pkg; Pkg.add("StaticArrays")
import Pkg; Pkg.add("SphericalHarmonics")

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
"""
    Airtable.request(method::AbstractString, cred::Credential, path[, headers = []]; query_kwargs...)

Make a request to the Airtable API.

Required arguments:

- `method`: one of "GET", "PUT", "POST", or "PATCH",
- `cred`: an [`Airtable.Credential`](@ref) containing your API key
- `path`: This often takes the form `<base id>/<table name>` or `<base id>/<table name>/<record id>`

Optional arguments:

- `headers`: a vector of `Pair`s representing header arguments (the things noted in the airtable API as -H arguments).
  Eg. `["Content-Type" => "application/json"]`

Query parameters are in the form of keyword arguments,
eg `filterByFormla = "NOT({Name} = '')", maxRecords=2`.
See [Airtable API](https://airtable.com/api) reference for more information.
"""
function request(method::AbstractString, cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...)
    method in ("GET", "PUT", "POST", "PATCH", "DELETE") || error("Invalid API method: $method")

    headers = append!(["Authorization"=> "Bearer $(cred.api_key)"], headers)
    query = []
    for (key, value) in query_kwargs
        isempty(value) && continue
        push!(query, string(key) => string(value))
    end
    path = joinpath("/", HTTP.escapeuri.(splitpath(lstrip(path, '/')))...)
    uri = HTTP.URI(; host="api.airtable.com", scheme="https", path, query)

    _ratelimit!(_ratelimiter)
    resp = isnothing(body) ? HTTP.request(method, uri, headers) :
                             HTTP.request(method, uri, headers, body)
    return JSON3.read(String(resp.body))
end

"""
    Airtable.get(cred::Credential, path[, headers=[], body=nothing]; query_kwargs...)

Shorthand for [`Airtable.request("GET", cred, path[, headers, body]; query_kwargs)`](@ref Airtable.request)
"""
get(cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...) = request("GET", cred, path, headers, body; query_kwargs...)

"""
    Airtable.put!(cred::Credential, path[, headers=[], body=nothing]; query_kwargs...)

Shorthand for [`Airtable.request("PUT", cred, path[, headers, body]; query_kwargs)`](@ref Airtable.request)
"""
put!(cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...) = request("PUT", cred, path, headers, body; query_kwargs...)

"""
    Airtable.post!(cred::Credential, path[, headers=[], body=nothing]; query_kwargs...)

Shorthand for [`Airtable.request("POST", cred, path[, headers, body]; query_kwargs)`](@ref Airtable.request)
"""
post!(cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...) = request("POST", cred, path, headers, body; query_kwargs...)

"""
    Airtable.patch!(cred::Credential, path[, headers=[], body=nothing]; query_kwargs...)

Shorthand for [`Airtable.request("PATCH", cred, path[, headers, body]; query_kwargs)`](@ref Airtable.request)
"""
patch!(cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...) = request("PATCH", cred, path, headers, body; query_kwargs...)

"""
    Airtable.delete!(cred::Credential, path[, headers=[], body=nothing]; query_kwargs...)

Shorthand for [`Airtable.request("DELETE", cred, path[, headers, body]; query_kwargs)`](@ref Airtable.request)
"""
delete!(cred::Credential, path::AbstractString, headers=[], body=nothing; query_kwargs...) = request("DELETE", cred, path, headers, body; query_kwargs...)


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

mutable struct ApplicationGatewayProtocol <: SwaggerModel

    function ApplicationGatewayProtocol(;)
        o = new()
        o
    end
end # type ApplicationGatewayProtocol

const _property_map_ApplicationGatewayProtocol = Dict{Symbol,Symbol}()
const _property_types_ApplicationGatewayProtocol = Dict{Symbol,String}()
Base.propertynames(::Type{ ApplicationGatewayProtocol }) = collect(keys(_property_map_ApplicationGatewayProtocol))
Swagger.property_type(::Type{ ApplicationGatewayProtocol }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ApplicationGatewayProtocol[name]))}
Swagger.field_name(::Type{ ApplicationGatewayProtocol }, property_name::Symbol) =  _property_map_ApplicationGatewayProtocol[property_name]

function check_required(o::ApplicationGatewayProtocol)
    true
end

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


text:
readintegers() = [parse(Int, x) for x in split(readline())]

function solve()
    n, k = readintegers()
    n = n % k
    ret = n
    n = abs(n - k)
    min(ret, n)
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
include("_models.jl")
include("ensembles.jl")


text:
################################################################################
#
#   METHODS FOR CONSTRUCTION OF PATHS INSIDE THE BZ OF A UNITCELL
#
#   STRUCTURE OF THE FILE
#
#   1) TYPE PATH
#       - type definition
#       - printInfo function
#       - path String function (e.g. X-G-M-X)
#
#   2) CONSTRUCTION FUNCTIONS PATH
#       - function to add points
#       - TODO function to remove points
#       - function to scale the total resolution
#       - function to set the total resolution
#
#   3) DEFAULT PATHS
#       - getDefaultPathTriangular
#       - getDefaultPathSquare
#       - getDefaultPathSquareOctagon
#       - getDefaultPathFCC
#
#   4) TODO FUNCTION TO CREATE A DEFAULT PATH BASED ON A UNITCELL OBJECT
#
#   5) FUNCTION TO PLOT A PATH (in both 2D and 3D)
#
################################################################################

################################################################################
#
#   The type Path
#
################################################################################
"""
    mutable struct Path

The type that contains information on a path (in momentum space). Fields are

    points             :: Array{Array{Float64, 1}, 1}
    point_names        :: Array{String, 1}
    segment_resolution :: Array{Int64, 1}

Note that the length of the `segment_resolution` array should be smaller than
the length of the `points` array by exactly 1.

New `Path` objects can be created with

    Path(points::Array{Array{Float64,1},1}, point_names::Array{String,1}, segment_resolution::Array{Int64,1})
    Path(points::Array{Array{Float64,1},1}, point_names::Array{String,1})
    Path()

or by one of the several default functions to create a default path.

# Examples

```julia-repl
julia> path = Path()
```
"""
mutable struct Path

    # Array of Point coordinates
    points::Array{Array{Float64,1}, 1}

    # Array of Point names
    point_names::Array{String, 1}

    # Array of resolutions (for later calculations)
    segment_resolution::Array{Int64, 1}

    # The Default constructor
    function Path(points::Array{Array{Float64,1},1}, point_names::Array{String,1}, segment_resolution::Array{Int64,1})
        return new(points, point_names, segment_resolution)
    end

    # The constructor for a new path without segment resolution information
    function Path(points::Array{Array{Float64,1},1}, point_names::Array{String,1})
        return new(points, point_names, ones(Int64, length(point_names)-1).*100)
    end

    # The constructor for a new path without information
    function Path()
        return new(Array{Float64,1}[], String[], Int64[])
    end

end

# export the type
export Path

################################################################################
#
#   Some functions to print information on the path
#
################################################################################

# Function to print some information on a path
"""
    printInfo(path::Path [; detailed::Bool])

prints (detailed) information on a `Path` object `path`. If detailed output is desired, the complete path will be printed.

# Examples

```julia-repl
julia> printInfo(path)
...

julia> printInfo(path, detailed=true)
...
```
"""
function printInfo(path::Path; detailed::Bool=false)
    # distinguish detailed vs. non-detailed
    if detailed
        # print the complete path
        println("Path overview:")
        # maybe already abort if no points or only one point in path
        if length(path.points) == 0
            println("   (no points defined)")
            return nothing
        elseif length(path.points) == 1
            println("  ($(1)) \"$(path.point_names[1])\" at $(path.points[1])  (only point in path)")
            return nothing
        end
        # alternate between points and segments
        for i in 1:length(path.points)-1
            # print the point
            println("  ($(i)) \"$(path.point_names[i])\" at $(path.points[i])")
            # print the outgoing segment
            println("         |  (resolution: $(path.segment_resolution[i]))")
        end
        # print the last point
        println("  ($(length(path.points))) \"$(path.point_names[length(path.points)])\" at $(path.points[length(path.points)])")
    else
        # not detailed, just give the number of segments and the total resolution
        println("Path contains $(length(path.points)) points ($(length(path.segment_resolution)) segments) with a total resolution of $(sum(path.segment_resolution)).")
        println("($(getPathString(path)))")
    end
end

# export the function
export printInfo

"""
    getPathString(path::Path)

compiles a string that represents the given path by chaining together all point names.

# Examples

```julia-repl
julia> getPathString(path)
"X--Gamma--K--M--X"

```
"""
function getPathString(path::Path)
    # create a new string
    path_string = ""
    # distinguish by point number
    if length(path.points) == 0
        # no points in path
        path_string = "(no points defined)"
    elseif length(path.points) == 1
        # one point in path
        path_string = "$(path.point_names[1])"
    else
        # alternate between points and segments
        for i in 1:length(path.points)-1
            # add the point and a segment
            path_string = "$(path_string)$(path.point_names[i])--"
        end
        # add the last point
        path_string = "$(path_string)$(path.point_names[end])"
    end
    # return the string
    return path_string
end

# export the function
export getPathString

################################################################################
#
#   CONSTRUCTION OF PATHS
#
################################################################################

# Add a point
"""
    addPointToPath!(path::Path, point::Array{Float64,1}, point_name::String [, resolution::Int64=100])

adds a new point to an existing `Path` object. The new point has to be specified by a name and a location.
Optionally, a resolution of the segment to the preceeding point can be given as well.
Note, that only a resolution will be added if there are already points in the path.
In this function, the path object will be changed and no new object will be created.

# Examples

```julia-repl
julia> addPointToPath!(path, [0.0, 0.0], "Gamma")

julia> addPointToPath!(path, [0.0, pi], "M", 150)
```
"""
function addPointToPath!(path::Path, point::Array{Float64,1}, point_name::String, resolution::Int64=100)
    # push the values into the lists
    push!(path.points, point)
    push!(path.point_names, point_name)
    # maybe push resolution, if there were already some points
    if length(path.points) > 1
        push!(path.segment_resolution, resolution)
    end
    # return nothing
    return nothing
end
# Function for when points including pi are added (not of type Float64)
function addPointToPath!(path::Path, point::Array, point_name::String, resolution::Int64=100)
    # push the values into the lists
    push!(path.points, point)
    push!(path.point_names, point_name)
    # maybe push resolution, if there were already some points
    if length(path.points) > 1
        push!(path.segment_resolution, resolution)
    end
    # return nothing
    return nothing
end

# export the function
export addPointToPath!

# scale the resolution by some factor
"""
    scaleResolution!(path::Path, factor::Float64)

scales all segment resolutions of the path by a factor and converts them back to `Int64`.
The path object will be changed and no new object will be created.

# Examples

```julia-repl
julia> scaleResolution!(path, 1.5)
```
"""
function scaleResolution!(path::Path, factor::Float64)
    # multiply all segments
    for s in 1:length(path.segment_resolution)
        path.segment_resolution[s] = round(Int64, path.segment_resolution[s]*factor)
    end
end

# export the function
export scaleResolution!

# set the total resolution
"""
    setTotalResolution!(path::Path, resolution::Int64)

scales all segment resolutions of the path by a factor to match the total resolution `resolution`
and converts them back to `Int64`. The new sum over all segments will give approximately `resolution` (up to float/int conversion).
The path object will be changed and no new object will be created.

# Examples

```julia-repl
julia> setTotalResolution!(path, 1500)
```
"""
function setTotalResolution!(path::Path, resolution::Int64)
    # determine the factor
    factor = resolution / sum(path.segment_resolution)
    # apply the factor
    scaleResolution!(path, factor)
end

# export the function
export setTotalResolution!

################################################################################
#
#   DEFAULT PATHS
#
################################################################################

"""
    getDefaultPathTriangular( [; resolution::Int64=900])

creates the default path for the triangular lattice. Points in this path are

    [0.0, 0.0]               (Gamma)
    [2*pi/sqrt(3.0), 2*pi/3] (K)
    [2*pi/sqrt(3.0), 0.0]    (M)
    [0.0, 0.0]               (Gamma)

Additionally, a resolution can be set so that the entire path in total has this resolution.

# Examples
```julia-repl
julia> path = getDefaultPathTriangular()
LatticePhysics.Path(...)

julia> path = getDefaultPathTriangular(resolution=1200)
LatticePhysics.Path(...)
```
"""
function getDefaultPathTriangular( ; resolution::Int64=900)
    # create a new path object
    path = Path()
    # insert points
    addPointToPath!(path, [0.0, 0.0], "Gamma")
    addPointToPath!(path, [2*pi/sqrt(3.0), 2*pi/3], "K")
    addPointToPath!(path, [2*pi/sqrt(3.0), 0.0], "M")
    addPointToPath!(path, [0.0, 0.0], "Gamma")
    # set the total resolution
    setTotalResolution!(path, resolution)
    # return the path
    return path
end

# export the function
export getDefaultPathTriangular

"""
    getDefaultPathSquare(version::Int64=1 [; resolution::Int64])

creates the default path for the square lattice. Points in this path are dependent on the version.

Additionally, a resolution can be set so that the entire path in total has this resolution.

# Versions

#### Version 1 (DEFAULT) - short

Points are given by

    [0.0, 0.0]  (Gamma)
    [ pi, 0.0]  (M)
    [ pi,  pi]  (K)
    [0.0, 0.0]  (Gamma)

#### Version 2 - long / extended

Points are given by

    [ pi, 0.0]  (M)
    [0.0, 0.0]  (Gamma)
    [ pi, -pi]  (K')
    [ pi, 0.0]  (M)
    [0.0,  pi]  (M')
    [ pi,  pi]  (K)
    [0.0, 0.0]  (Gamma)

# Examples
```julia-repl
julia> path = getDefaultPathSquare()
LatticePhysics.Path(...)

julia> path = getDefaultPathSquare(2)
LatticePhysics.Path(...)

julia> path = getDefaultPathSquare(resolution=1200)
LatticePhysics.Path(...)
```
"""
function getDefaultPathSquare(version::Int64=1; resolution::Int64=1000)
    # create a new path object
    path = Path()
    # distinguish between version
    if version == 1
        # insert points
        addPointToPath!(path, [0.0, 0.0], "Gamma")
        addPointToPath!(path, [ pi, 0.0], "M")
        addPointToPath!(path, [ pi,  pi], "K")
        addPointToPath!(path, [0.0, 0.0], "Gamma")
    elseif version == 2
        # insert points
        addPointToPath!(path, [ pi, 0.0], "M")
        addPointToPath!(path, [0.0, 0.0], "Gamma")
        addPointToPath!(path, [ pi, -pi], "K'")
        addPointToPath!(path, [ pi, 0.0], "M")
        addPointToPath!(path, [0.0,  pi], "M'")
        addPointToPath!(path, [ pi,  pi], "K")
        addPointToPath!(path, [0.0, 0.0], "Gamma")
    else
        println("version $(version) unknown")
    end
    # set the total resolution
    setTotalResolution!(path, resolution)
    # return the path
    return path
end

# export the function
export getDefaultPathSquare

"""
    getDefaultPathSquareOctagon( [; resolution::Int64=900])

creates the default path for the square octagon lattice (version 2). Points in this path are

    [ 0.0,  0.0]                                 (Gamma)
    [ 2.0,  0.0] * (pi / (1.0 + 1.0/sqrt(2.0)))  (K)
    [-1.0, -1.0] * (pi / (1.0 + 1.0/sqrt(2.0)))  (M)
    [ 0.0,  0.0]                                 (Gamma)

Additionally, a resolution can be set so that the entire path in total has this resolution.

# Examples
```julia-repl
julia> path = getDefaultPathSquareOctagon()
LatticePhysics.Path(...)

julia> path = getDefaultPathSquareOctagon(resolution=1200)
LatticePhysics.Path(...)
```
"""
function getDefaultPathSquareOctagon( ; resolution::Int64=900)
    # create a new path object
    path = Path()
    # insert points
    addPointToPath!(path, [0.0, 0.0], "Gamma")
    addPointToPath!(path, [2,  0].*(pi / (1.0 + 1.0/sqrt(2.0))), "K")
    addPointToPath!(path, [1, -1].*(pi / (1.0 + 1.0/sqrt(2.0))), "M")
    addPointToPath!(path, [0.0, 0.0], "Gamma")
    # set the total resolution
    setTotalResolution!(path, resolution)
    # return the path
    return path
end

# export the function
export getDefaultPathSquareOctagon

"""
    getDefaultPathFCC( [; resolution::Int64=1200])

creates the default path for the FCC lattice (3D). Points in this path are

    [   0.0,    0.0,  0.0]   (Gamma)
    [  2*pi,    0.0,  0.0]   (X)
    [  2*pi,     pi,  0.0]   (W)
    [    pi,     pi,   pi]   (L)
    [   0.0,    0.0,  0.0]   (Gamma)
    [3*pi/2, 3*pi/2,  0.0]   (K)
    [  2*pi,     pi,  0.0]   (W)

Additionally, a resolution can be set so that the entire path in total has this resolution.

# Examples
```julia-repl
julia> path = getDefaultPathFCC()
LatticePhysics.Path(...)

julia> path = getDefaultPathFCC(resolution=600)
LatticePhysics.Path(...)
```
"""
function getDefaultPathFCC( ; resolution::Int64=1200)
    # create a new path object
    path = Path()
    # insert points
    addPointToPath!(path, [   0.0,    0.0,  0.0], "Gamma")
    addPointToPath!(path, [  2*pi,    0.0,  0.0], "X")
    addPointToPath!(path, [  2*pi,     pi,  0.0], "W")
    addPointToPath!(path, [    pi,     pi,   pi], "L")
    addPointToPath!(path, [   0.0,    0.0,  0.0], "Gamma")
    addPointToPath!(path, [3*pi/2, 3*pi/2,  0.0], "K")
    addPointToPath!(path, [  2*pi,     pi,  0.0], "W")
    # set the total resolution
    setTotalResolution!(path, resolution)
    # return the path
    return path
end

# export the function
export getDefaultPathFCC

# PLOTTING IN 2D (not exported)
function plotPath2D(
            path::Path;
            plot_color="r",
            new_figure::Bool=true,
            showPlot::Bool=true
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # only if new figure is desired
    if new_figure

        # configure plot environment
        rc("font", family="serif")

        # create a new figure
        fig = figure()

    end

    ###########################
    #   PLOT BRILLOUIN ZONE
    ###########################

    # compile lists of x and y values
    x_values = Float64[p[1] for p in path.points]
    y_values = Float64[p[2] for p in path.points]

    # STEP 1 - scatter the points
    scatter(x_values, y_values, color=plot_color)

    # STEP 2 - draw all lines of edges
    plot(
            x_values, y_values, "-",
            color=plot_color
        )

    # STEP 3 - set all labels
    for i in 1:length(path.points)
        text(
            x_values[i], y_values[i], path.point_names[i],
            ha="left", va="top", size=15, color=plot_color
        )
    end

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return gcf()

end

# PLOTTING IN 3D (not exported)
function plotPath3D(
            path::Path;
            plot_color="r",
            new_figure::Bool=true,
            showPlot::Bool=true
        )

    ###########################
    #   INITIAL SETTINGS
    ###########################

    # only if new figure is desired
    if new_figure

        # configure plot environment
        rc("font", family="serif")

        # create a new figure
        fig = figure()

    end

    ###########################
    #   PLOT BRILLOUIN ZONE
    ###########################

    # compile lists of x and y values
    x_values = Float64[p[1] for p in path.points]
    y_values = Float64[p[2] for p in path.points]
    z_values = Float64[p[3] for p in path.points]

    # STEP 1 - scatter the points
    scatter3D(x_values, y_values, z_values, color=plot_color)

    # STEP 2 - draw all lines of edges
    plot3D(
            x_values, y_values, z_values, "-",
            color=plot_color
        )

    # STEP 3 - set all labels
    for i in 1:length(path.points)
        text3D(
            x_values[i], y_values[i], z_values[i], path.point_names[i],
            ha="left", va="top", size=15, color=plot_color
        )
    end

    ###########################
    #   CONFIGURE AXIS & TITLE
    ###########################

    # equal axis
    gca()[:set_aspect]("equal")

    ###########################
    #   FINISH THE PLOT
    ###########################

    # tighten the layout
    tight_layout()

    # maybe show the plot
    if showPlot
        show()
    end

    # return the figure object
    return gcf()

end

# GENERAL PLOTTING (exported)

"""
    plotPath(
                path::Path
             [; plot_color="r",
                new_figure::Bool=true,
                showPlot::Bool=true ]
            )

Plots the Path using `PyPlot`. First, all corner points are scattered using `scatter3d`.
Second, all vertices are plotted using `plot3d`.
Third, all labels are put at the location of the points.

For specifying details while plotting, optional keywords can be passed:
- `plot_color` is the color of the BZ in the plot.
- `new_figure` specifies if a new figure should be created before plotting (or if the current open figure should be used)
- `showPlot` specifies if the plot should be opened after plotting.

# Examples

```julia-repl
julia> plotPath(path)
PyPlot.Figure(...)

julia> plotPath(path, plot_color="m")
PyPlot.Figure(...)

```
"""
function plotPath(
            path::Path;
            plot_color="r",
            new_figure::Bool=true,
            showPlot::Bool=true
        )
    # check if points are 3D or 2D
    if length(path.points[1]) == 2
        return plotPath2D(
            path,
            plot_color=plot_color,
            new_figure=new_figure,
            showPlot=showPlot
        )
    elseif length(path.points[1]) == 3
        return plotPath3D(
            path,
            plot_color=plot_color,
            new_figure=new_figure,
            showPlot=showPlot
        )
    else
        println("Bad dimension of k space: $(length(path.points[1]))")
    end

end
export plotPath


text:
InputFolder = './Images/CrazyCells/';
OutputFolder = './Results/Images/CrazyCells/';
Lbl = 1;

@iA = '*_C00*.tif';

@fxg_sLoGLocMax [iA] > [S];
params.GRad = 15;
params.LocalMaxBox = [9 9];
params.ThrLocalMax = 0.025;
params.MinLocThr = 0;
/endf

@fxgs_lSeededPolarGeod [iA, S] > [L];
params.ThDiv = 180;
params.ObjRad = 35;
params.MinArea = 175;
params.Method = 'simple';
/endf

/show iA > S;
/show iA > L;
/keep L > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
parsetree(org::Org, maxdepth=-1) = parsetree(stdout, org, maxdepth)
parsetree(org::OrgComponent, maxdepth=-1) = parsetree(stdout, org, maxdepth)

function parsetree(io::IO, org::Org, maxdepth=-1, depth=0)
    printstyled(io, "Org Parse Tree", color=:yellow)
    if length(org.contents) == 0
        printstyled(io, " (empty)\n", color=:light_black)
    else
        print(io, '\n')
        parsetree.(Ref(io), org.contents, maxdepth, depth+1)
    end
    nothing
end

function printstructure(io::IO, component::AbstractString, description,
                        contents, maxdepth, depth, color=:yellow, bold=false)
    printstyled(io, ' '^4depth, component; bold, color)
    if !isnothing(description)
        printstyled(io, ' ', description, color=:green)
    end
    if isnothing(contents)
        print(io, '\n')
    elseif contents isa Vector && length(contents) == 0
        printstyled(io, " (empty)\n", color=:light_black)
    elseif ! (contents isa Vector{<:OrgComponent})
        print(io, '\n')
    elseif depth == maxdepth
        printstyled(io, " [hidden]\n", color=:light_black)
    else
        print(io, '\n')
    end
    if !isnothing(contents)
        parsetree.(Ref(io), contents, maxdepth, depth+1)
    end
    nothing
end

parsetree(io::IO, h::Heading, maxdepth::Integer, depth::Integer) =
    printstructure(io, "Heading",
                   '(' * '*'^h.level * ") " * string(Paragraph(h.title)),
                   filter(!isnothing, [h.planning, h.properties, h.section]),
                   maxdepth, depth, :yellow, true)

parsetree(io::IO, s::Section, maxdepth::Integer, depth::Integer) =
    printstructure(io, "Section", nothing, s.contents, maxdepth, depth, :yellow, true)

parsetree(io::IO, l::L, maxdepth::Integer, depth::Integer) where {L <: List} =
    printstructure(io, string(L), l.items[1].bullet, l.items, maxdepth, depth, :magenta, true)

function parsetree(io::IO, component::T, maxdepth::Integer, depth::Integer=0) where {T <: OrgComponent}
    contents = if applicable(iterate, component)
        collect(component)
    end
    bold, color = if T <: GreaterElement
        (false, :magenta)
    elseif T <: Element
        (false, :blue)
    else
        (false, :default)
    end
    printstructure(io, structurename(component), structuredesc(component),
                   contents, maxdepth, depth, color, bold)
end

structurename(::C) where {C<:OrgComponent} = string(nameof(C))
structuredesc(_::OrgComponent) = nothing

structuredesc(d::Drawer) = d.name
structuredesc(k::CitationReference) = string('@', k.key)
structuredesc(c::Citation) =
    string(something(c.style[1], "⋅"), " ", something(c.style[2], "⋅"))
structuredesc(f::FootnoteDefinition) = f.label
structuredesc(k::Keyword) = k.key
structuredesc(a::AffiliatedKeyword) = a.key
structuredesc(f::FootnoteReference) = something(f.label, "")
structuredesc(n::NodeProperty) = string(n.name, if n.additive "+" else "" end)
structuredesc(s::SourceBlock) = s.lang
structuredesc(b::SpecialBlock) = b.name
structuredesc(l::Link) = l.path.protocol
structuredesc(m::Macro) = m.name
structurename(t::TextMarkup) = "Text" * uppercasefirst(string(t.formatting))
structuredesc(t::TextMarkup{<:AbstractString}) =
    string(structuredesc(TextPlain(t.contents)))
structuredesc(t::TextPlain) = if length(t.text) <= 40
    sprint(show, t.text)
else
    string(sprint(show, t.text[1:20])[1:end-1], " … ",
           sprint(show, t.text[end-20:end])[2:end])
end


text:
using DSGE
using HDF5, Base.Test
using DataFrames

include("../util.jl")
path = dirname(@__FILE__)

# Set up model for testing
m = Model990()
m.testing=true

# Read in the covariance matrix for Metropolis-Hastings and reference parameter draws
file = h5open("$path/../reference/metropolis_hastings.h5","r")
propdist_cov  = read(file, "propcov")
ref_draws     = read(file, "ref_draws")
ref_cov       = read(file, "ref_cov")
close(file)

# Read in the data
df = load_data(m; try_disk=true, verbose=:none)

# Set up and run metropolis-hastings
specify_mode!(m, inpath(m, "user", "paramsmode.h5"), verbose=:none)
specify_hessian(m, inpath(m, "user", "hessian.h5"), verbose=:none)

# Note Distributions (used in earlier test) also exports `estimate`.
DSGE.estimate(m, df; verbose=:none, proposal_covariance = propdist_cov)

# Read in the parameter draws and covariance just generated from estimate.
test_draws = h5open(rawpath(m, "estimate", "mhsave.h5"), "r") do file
    read(file, "mhparams")
end
test_cov = h5open(workpath(m, "estimate", "parameter_covariance.h5"), "r") do file
    read(file, "mhcov")
end

return

# Test that the parameter draws and covariance matrices are equal
@test_matrix_approx_eq ref_draws test_draws
@test_matrix_approx_eq ref_cov test_cov

# Test that the fixed parameters are all fixed
for fixed_param in [:δ, :λ_w, :ϵ_w, :ϵ_p, :g_star]
    @test test_draws[:,m.keys[fixed_param]] == fill(Float32(m[fixed_param].value), 100)
end

# Make sure that compute_moments runs appropriately
compute_moments(m, verbose=:none)

nothing


text:
const ORDER_URL = "https://api-fxpractice.oanda.com/v3/accounts/"*ACCOUNT*"/orders"

function order(pair::AbstractString, lotsize::Int8, stoploss::Int8, takeprofit::Int8, volatility::Int8, price::Float16)
	order = Dict{String,Dict}(
		order => Dict{String,String}(
			type => "MARKET_IF_TOUCHED", #["MARKET","LIMIT","STOP","MARKET_IF_TOUCHED","TAKE_PROFIT","STOP_LOSS","TRAILING_STOP_LOSS"]
			instrument => pair,
			# positive = long, negative = short
			units => string(lotsize),
			# The time-in-force requested for the Market Order. Restricted to FOK or IOC for a MarketOrder.
			# GTC The Order is “Good unTil Cancelled”
			# GTD The Order is “Good unTil Date” and will be cancelled at the provided time
			# GFD The Order is “Good For Day” and will be cancelled at 5pm New York time
			# FOK The Order must be immediately “Filled Or Killed”
			# IOC The Order must be “Immediatedly paritally filled Or Cancelled”
			timeInForce => "FOK",
			# How Positions are modified when Order is filled.
			# OPEN_ONLY only allow Positions to be opened or extended.
			# REDUCE_FIRST always fully reduce an existing Position before opening a new Position.
			# REDUCE_ONLY only reduce an existing Position.
			# DEFAULT use REDUCE_FIRST behaviour for non-client hedging Accounts, and OPEN_ONLY behaviour for client hedging Accounts.
			positionFill => "REDUCE_FIRST",
			# worst price acceptable
			# priceBound => string(price),
		)
	)

	Requests.post(ORDER_URL; headers = HEADERS; json = order)
end

max_stars_repo_name: cadwallader/Oanda

max_issues_repo_name: cadwallader/Oanda

max_forks_repo_name: cadwallader/Oanda


text:
using ProximalOperators
using Base.Test

TOL_ASSERT = 1e-12

# measures time and returns the result of the call method
function call_test(f, x)
  print("* call                        : "); @time fx = f(x)
  @test typeof(fx) == eltype(real(x))
  return fx
end

# measures time of the three possible calls to prox and prox!
# then tests equality of the results and returns them if they agree
function prox_test(f, x, gamma::Real=1.0)
  print("* prox                        : "); @time yf, fy = prox(f, x, gamma)
  print("* prox! (preallocated output) : "); yf_prealloc = copy(x); @time fy_prealloc = prox!(f, x, yf_prealloc, gamma)
  print("* prox! (in place)            : "); yf_inplace = copy(x); @time fy_inplace = prox!(f, yf_inplace, gamma)
  @test typeof(fy) == eltype(real(x))
  @test vecnorm(yf_prealloc - yf, Inf)/(1+vecnorm(yf, Inf)) <= TOL_ASSERT
  @test vecnorm(yf_inplace - yf, Inf)/(1+vecnorm(yf, Inf)) <= TOL_ASSERT
  if ProximalOperators.is_prox_accurate(f)
    @test fy_prealloc == fy || abs(fy_prealloc - fy)/(1+abs(fy)) <= TOL_ASSERT
    @test fy_inplace == fy || abs(fy_inplace - fy)/(1+abs(fy_inplace)) <= TOL_ASSERT
  end
  return yf, fy
end

println("*********************************************************************")
include("test_utilities.jl")
println("*********************************************************************")
include("test_calls.jl")
println("*********************************************************************")
include("test_calculus.jl")
println("*********************************************************************")
include("test_SeparableSum.jl")
include("test_SlicedSeparableSum.jl")
println("*********************************************************************")
include("test_equivalences.jl")
println("*********************************************************************")
include("test_condition.jl")
println("*********************************************************************")
include("test_results.jl")
println("*********************************************************************")


text:
module TestModule

end


text:
@testset "Properties" begin

    # Fake atoms
    atoms = [
        Atom(1, 1.008, [5.0, 5.0, 0.0])
        Atom(2, 1.008, [5.0, -5.0, 0.0])
        Atom(3, 1.008, [-5.0, -5.0, 0.0])
        Atom(4, 1.008, [-5.0, 5.0, 0.0])
    ]

    @test Molecules.center_of_mass(atoms) == [0.0, 0.0, 0.0]
    @test Molecules.nuclear_repulsion(atoms) ≈ 1.6816285798739845
end


text:
# This file contains basic structures and method dispatched when the network
# is evaluated with respect to a state with an associated lookup table.
# At the moment, those have a bad performance and should not be used

# !!!!! Ignore this file! Lookup tables are not really needed....

abstract type NNLookUp
    #valid
end

"""
    lookup(net)

Constructs the `NNCache{typeof(net)}` object that holds the cache for this network.
If it has not been implemented returns nothing.
"""
lookup(net::NeuralNetwork) = nothing
lookup(net::CachedNet) = lookup(net.net)


text:
mutable struct TspLibProblem
    name::String
    numcities::Int
    weights::Matrix{Float64}
end

size(t::TspLibProblem) = t.numcities

function resetweights!(t::TspLibProblem, nc::Int)
    t.numcities = nc
    t.weights = zeros(Float64, nc, nc)
end

function cost(p::TspLibProblem, visitorder::Vector{Int})
    startcity = city = first(visitorder)
    cost = 0.0
    for i in 2:length(visitorder)
        nextcity = visitorder[i]
        cost += p.weights[city, nextcity]
        city = nextcity
    end
    cost += p.weights[city, startcity]
    return cost
end

function whenmatch(fn, re, s)
    m = match(re, s)
    if !isnothing(m)
        fn(m)
        return true
    end
    return false
end

function parse_edge_weights!(p, lines, lineindex)
    row = col = 1
    while !isnothing(match(r"^(\s*(\d+))+\s*$", lines[lineindex]))
        map(split(lines[lineindex], r"\s+")) do dstr
            if length(dstr) > 0
                p.weights[row, col] = p.weights[col, row] = parse(Int, strip(dstr))
                col += 1
                if col > row
                    row += 1
                    col = 1
                end
            end
        end
        lineindex += 1
    end
    return lineindex - 1
end

filename = "dantzig42.tsp"

function read_tsplib_file(filename)
    p = TspLibProblem("", 1, zeros(Float64, 1, 1))
    lines = readlines(filename)
    lineindex = 0
    while lineindex < length(lines)
        lineindex += 1
        l = lines[lineindex]
        whenmatch(m -> p.name = strip(m[1]),
            r"\s*NAME\s*:\s*(.+)$", l) && continue
        whenmatch(m -> resetweights!(p, parse(Int, m[1])),
            r"\s*DIMENSION\s*:\s*(\d+)$", l) && continue
        whenmatch(m -> lineindex = parse_edge_weights!(p, lines, lineindex+1),
            r"\s*EDGE_WEIGHT_SECTION\s*$", l)
    end
    return p
end


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

p = [0.0, 0.0]

mergemapper = CutCell.MergeMapper()

south = CutCell.cell_map_to_south()
@test allapprox(south(p), [0.0, -2.0])
@test allapprox(mergemapper[1](p), [0.0, -2.0])

southeast = CutCell.cell_map_to_south_east()
@test allapprox(southeast(p), [2.0, -2.0])
@test allapprox(mergemapper[5](p), [2.0, -2.0])

east = CutCell.cell_map_to_east()
@test allapprox(east(p), [2.0, 0.0])
@test allapprox(mergemapper[2](p), [2.0, 0.0])

northeast = CutCell.cell_map_to_north_east()
@test allapprox(northeast(p), [2.0, 2.0])
@test allapprox(mergemapper[6](p), [2.0, 2.0])

north = CutCell.cell_map_to_north()
@test allapprox(north(p), [0.0, 2.0])
@test allapprox(mergemapper[3](p), [0.0, 2.0])

northwest = CutCell.cell_map_to_north_west()
@test allapprox(northwest(p), [-2.0, 2.0])
@test allapprox(mergemapper[7](p), [-2.0, 2.0])

west = CutCell.cell_map_to_west()
@test allapprox(west(p), [-2.0, 0.0])
@test allapprox(mergemapper[4](p), [-2.0, 0.0])

southwest = CutCell.cell_map_to_south_west()
@test allapprox(southwest(p), [-2.0, -2.0])
@test allapprox(mergemapper[8](p), [-2.0, -2.0])

function test_linear_merge_cell_assembly()

    polyorder = 1
    numqp = 2

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)
    nf = CutCell.number_of_basis_functions(basis)
    mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [2, 1], nf)

    normal = [1.0, 0.0]
    x0 = [1.1, 0.0]
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> plane_distance_function(x, normal, x0), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    mergecutmesh = CutCell.MergeCutMesh(cutmesh)
    mergemapper = CutCell.MergeMapper()
    @test CutCell.number_of_nodes(mergecutmesh) == 10

    CutCell.merge_cells_in_mesh!(mergecutmesh,cellquads,interfacequads,mergemapper)
    mergedmesh = CutCell.MergedMesh(mergecutmesh)
    @test CutCell.number_of_nodes(mergedmesh) == 8

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    dx = 0.1
    e11 = dx / 2.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22
    s11 = (lambda + 2mu) * e11 + lambda * e22

    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, mergedmesh)
    interfacecondition =
        CutCell.coherent_interface_condition(basis, interfacequads, stiffness, mergedmesh, penalty)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, mergedmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, mergedmesh)

    matrix = CutCell.make_sparse(sysmatrix, mergedmesh)
    rhs = CutCell.rhs(sysrhs, mergedmesh)

    CutCell.apply_dirichlet_bc!(matrix, rhs, 9, 0.0)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 10, 0.0)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 11, 0.0)

    CutCell.apply_dirichlet_bc!(matrix, rhs, 5, dx)
    CutCell.apply_dirichlet_bc!(matrix, rhs, 7, dx)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    testdisp = [dx/2  dx/2  dx  dx  0.  0.  dx/2  dx/2
                0.    dy    0.  dy  0.  dy  0.    dy]
    @test allapprox(disp,testdisp,1e2eps())
end

function test_curved_merge_assembly()
    polyorder = 2
    numqp = 4

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)
    nf = CutCell.number_of_basis_functions(basis)
    mesh = CutCell.Mesh([0.0, 0.0], [2.0, 1.0], [2, 1], nf)

    xc = [2.0, 0.5]
    rad = 0.9
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> circle_distance_function(x, xc, rad), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    mergemapper = CutCell.MergeMapper()
    mergecutmesh = CutCell.MergeCutMesh(cutmesh)
    @test CutCell.number_of_nodes(mergecutmesh) == 24

    CutCell.merge_cells_in_mesh!(mergecutmesh,cellquads,interfacequads,mergemapper)
    mergedmesh = CutCell.MergedMesh(mergecutmesh)
    @test CutCell.number_of_nodes(mergedmesh) == 18

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    dx = 0.1
    e11 = dx / 2.0
    e22 = -lambda / (lambda + 2mu) * e11
    dy = e22
    s11 = (lambda + 2mu) * e11 + lambda * e22

    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, mergedmesh)
    interfacecondition =
        CutCell.coherent_interface_condition(basis, interfacequads, stiffness, mergedmesh, penalty)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, mergedmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, mergedmesh)

    matrix = CutCell.make_sparse(sysmatrix, mergedmesh)
    rhs = CutCell.rhs(sysrhs, mergedmesh)

    CutCell.apply_dirichlet_bc!(matrix, rhs, [10, 11, 12], 1, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [10], 2, 0.0, 2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [7, 8, 9], 1, dx, 2)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    testdisp = [dx/2  dx/2  dx/2  3dx/4  3dx/4  3dx/4  dx  dx    dx  0.  0.    0.   dx/4  dx/4  dx/4  dx/2  dx/2  dx/2
                0.    dy/2  dy    0.     dy/2   dy     0.  dy/2  dy  0.  dy/2  dy   0.    dy/2  dy    0.    dy/2  dy]
    @test allapprox(disp,testdisp,1e2eps())
end

function test_four_cell_vertical_tension()
    polyorder = 1
    numqp = 3

    basis = TensorProductBasis(2, polyorder)
    levelset = InterpolatingPolynomial(1, basis)
    nf = CutCell.number_of_basis_functions(basis)
    mesh = CutCell.Mesh([0.0, 0.0], [2.0, 2.0], [2, 2], nf)

    normal = [0.0, 1.0]
    x0 = [0.0, 1.1]
    levelsetcoeffs =
        CutCell.levelset_coefficients(x -> plane_distance_function(x, normal, x0), mesh)

    cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
    cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
    interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

    mergemapper = CutCell.MergeMapper()
    mergecutmesh = CutCell.MergeCutMesh(cutmesh)

    CutCell.merge_cells_in_mesh!(mergecutmesh,cellquads,interfacequads,mergemapper)
    mergedmesh = CutCell.MergedMesh(mergecutmesh)

    lambda, mu = (1.0, 2.0)
    penalty = 1.0
    dy = 0.1
    e22 = dy / 2.0
    e11 = -lambda / (lambda + 2mu) * e22
    dx = e11*2.0

    stiffness = CutCell.HookeStiffness(lambda, mu, lambda, mu)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition =
        CutCell.coherent_interface_condition(basis, interfacequads, stiffness, cutmesh, penalty)

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, mergedmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, mergedmesh)

    matrix = CutCell.make_sparse(sysmatrix, mergedmesh)
    rhs = CutCell.rhs(sysrhs, mergedmesh)

    CutCell.apply_dirichlet_bc!(matrix, rhs, [7,9,11],2,0.0,2)
    CutCell.apply_dirichlet_bc!(matrix, rhs, [7],1,0.0,2)

    CutCell.apply_dirichlet_bc!(matrix, rhs, [2,4,6],2,dy,2)

    sol = matrix \ rhs
    disp = reshape(sol, 2, :)

    testdisp = [0.   0.  dx/2  dx/2  dx    dx  0.  0.   dx/2  dx/2  dx  dx
                dy/2 dy  dy/2  dy    dy/2  dy  0.  dy/2 0.    dy/2  0.  dy/2]
    @test allapprox(disp,testdisp,1e2eps())
end

test_linear_merge_cell_assembly()
test_curved_merge_assembly()
test_four_cell_vertical_tension()


text:
module CharibdeOptim

export constraint, charibde_min, charibde_max, ibc_maximise, diffevol_maximise
export ibc_minimise, diffevol_minimise

using IntervalArithmetic
using Distributed
using IntervalOptimisation: SortedVector, filter_elements!
using IntervalConstraintProgramming
using ModelingToolkit
using MathOptInterface
using StaticArrays
using MacroTools

import Base: invokelatest, push!
import ForwardDiff: gradient

struct Constraint{T}
   bound::Interval{T}
   C::Contractor
end

function constraint(vars, constraint_expr::Operation, bound::Interval{T}; epsilon = 1e-4) where{T}
   C = Contractor(vars, constraint_expr)
   if diam(bound) == 0.0
       bound = Interval(bound.lo - epsilon, bound.hi + epsilon)
   end
   Constraint{T}(bound, C)
end

mutable struct Information
    de_to_ibc::Int
    ibc_to_de::Int
    iterations::Int
end

include("IBC.jl")
include("DifferentialEvolution.jl")
include("ConstrainedDifferentialEvolution.jl")
include("ConstrainedIBC.jl")
include("BoundEnsure.jl")
include("GenerateRandom.jl")
include("SecondOrder.jl")

function create_channels(X::IntervalBox{N, T}, workers::Int64) where{N, T}
    if workers > 1
        return (RemoteChannel(()->Channel{Tuple{IntervalBox{N, T}, T}}(1)), RemoteChannel(()->Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(1)))
    else
        return (Channel{Tuple{IntervalBox{N, T}, T}}(1), Channel{Tuple{SVector{N, T}, T, Union{Nothing, IntervalBox{N, T}}}}(1))
    end
end

"""Usage:
```
f = X->((x,y)=X;x^3 + 2y + 5)
A = IntervalBox(2..4, 2..3)
(global_min, minimisers, info) = charibde_min(f, A)
(global_max, maximisers, info) = charibde_max(f, A)

charibde_min/charibde_max find the global minimum/maximum value of the function in given search space by using two algorithms Interval Bound & Contract(IBC) and Differential Evolution inparallel
TODO: Fix the parallel implementation of charibde_min/max for Constrained Optimisation
```
"""
function charibde_min(f::Function, X::IntervalBox{N,T}; workers = 2, tol = 1e-6, np = N*10, debug = false) where{N,T}

    worker_ids = Distributed.workers()
    if workers > 1
        if worker_ids[1] == 1
            error("Not enough workers available: Add one worker and load the package on it")
        end

        (chnl1, chnl2) = create_channels(X, workers)                     #IBC recieve element from chnl1 and DiffEvolution from chnl2

        r1 = remotecall(diffevol_minimise, worker_ids[1], f, X, chnl1, chnl2, np = np)
        r2 = remotecall(ibc_minimise, worker_ids[2], f, X, ibc_chnl = chnl1, diffevol_chnl = chnl2, tol = tol, debug = debug)
        return fetch(r2)
    else
        (chnl1, chnl2) = create_channels(X, workers)

        r1 = @async diffevol_minimise(f, X, chnl1, chnl2, np = np)
        r2 = @async ibc_minimise(f, X, ibc_chnl = chnl1, diffevol_chnl = chnl2, tol = tol, debug = debug)
        return fetch(r2)
    end
end

function charibde_min(f::Function, X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}; workers = 2, tol = 1e-6, np = N*10, debug = false) where{N,T}

    worker_ids = Distributed.workers()
    if workers > 1
        if nprocs() < 3
            error("Not enough workers available: Session should have atleast 2 workers more")
        end

        (chnl1, chnl2) = create_channels(X, workers)

        r1 = remotecall(diffevol_minimise, worker_ids[1], f, X, constraints, chnl1, chnl2, np = np)
        r2 = remotecall(ibc_minimise, worker_ids[2], f, X, constraints, ibc_chnl = chnl1, diffevol_chnl = chnl2, tol = tol, debug = debug)
        return fetch(r2)
    else
        (chnl1, chnl2) = create_channels(X, workers)

        r1 = @async diffevol_minimise(f, X, constraints, chnl1, chnl2, np = np)
        r2 = @async ibc_minimise(f, X, constraints, ibc_chnl = chnl1, diffevol_chnl = chnl2, tol= tol, debug = debug)
        return fetch(r2)
    end
end

function charibde_max(f::Function, X::IntervalBox{N,T}; workers = 2, tol = 1e-6, np = N*10, debug = false) where{N,T}
    bound, minimizers, info = charibde_min(x -> -f(x), X, workers = workers, tol = tol, np = np, debug = debug)
    return -bound, minimizers, info
end

function charibde_max(f::Function, X::IntervalBox{N,T}, constraints::Vector{Constraint{T}}; workers = 2, tol = 1e-6, np = N*10, debug = false) where{N,T}
    bound, minimizers, info = charibde_min(x -> -f(x), X, constraints, workers = workers, tol = tol, np = np, debug = debug)
    return -bound, minimizers, info
end

include("MOI_wrapper.jl")

end


text:
using Documenter, OpticsPolynomials

makedocs(modules = [OpticsPolynomials], doctest = false, sitename = "OpticsPolynomials")
deploydocs(repo = "github.com/JuliaOptics/OpticsPolynomials.jl.git", devbranch="main")


text:
_annotations(m::Metabolite) = m.annotations
_annotations(r::Reaction) = r.annotations
_annotations(g::Gene) = g.annotations

"""
    annotation_index(
        xs::AbstractDict{String};
        annotations = _annotations,
    )::Dict{String,Dict{String,[String]}}

Extract annotations from a dictionary of items `xs` and build an index that
maps annotation "kinds" (e.g. `"PubChem"`) to the mapping from the annotations
(e.g.  `"COMPOUND_12345"`) to item IDs that carry the annotations.

Function `annotations` is used to access the `Annotations` object in the
dictionary values.

This is extremely useful for finding items by annotation data.
"""
function annotation_index(
    xs::AbstractDict{String};
    annotations = _annotations,
)::Dict{String,Dict{String,Set{String}}}
    res = Dict{String,Dict{String,Set{String}}}()
    for (n, ax) in xs
        a = annotations(ax)
        for (k, anns) in a
            if !haskey(res, k)
                res[k] = Dict{String,Set{String}}()
            end
            for v in anns
                if !haskey(res[k], v)
                    res[k][v] = Set([n])
                else
                    push!(res[k][v], n)
                end
            end
        end
    end
    res
end

"""
    ambiguously_identified_items(
        index::Dict{String,Dict{String,[String]}},
    )::Vector{String}

Find items (genes, metabolites, ...) from the annotation index that are
identified non-uniquely by at least one of their annotations.

This often indicates that the items are duplicate or miscategorized.
"""
function ambiguously_identified_items(
    index::Dict{String,Dict{String,Set{String}}},
)::Set{String}
    res = Set{String}()
    for (_, idents) in index
        for (_, items) in idents
            if length(items) > 1
                push!(res, items...)
            end
        end
    end
    res
end