text:
function varstruct_construct(modul::Module, T, args_field, args_defaultvalue, args_type, shared_var = false)
    return quote
        $(get_struct_definition(modul, T, args_field, args_defaultvalue, args_type, shared_var))
    end
end

function get_struct_definition(modul::Module, T, args_field, args_defaultvalue, args_type, shared_var = false)
    esc_T = esc(T)
    T_declaration = Symbol("___", T , "_declaration___")

    if isdefined(modul, T) && isa(getfield(modul, T), Type)
        # if already defined don't define it again
        out = quote
            # Struct Initialization (returns an instance)
            $(initialize_struct(modul, T, args_field, args_defaultvalue, args_type, T_declaration))
        end
    else
        strdef = :(mutable struct $T
            fieldtable::FieldTable
        end)

        out = quote

            # Struct Declaration
            # we need to eval manually to support definition inside non-toplevel scope (like functions). Otherwise this is not necessary
            Base.eval($modul, $strdef)

            # Struct Interface
            $(get_struct_interface(esc_T))

            # Struct Constructor
            $(get_struct_constructor(modul, T, T_declaration, shared_var))

            # Show
            $(show_struct(esc_T))

            # Struct Initialization (returns an instance)
            # @eval $modul $strinit
            $(initialize_struct(modul, T, args_field, args_defaultvalue, args_type, T_declaration))
        end

    end
    return out
end

function get_struct_interface(esc_T)
    return quote
        Base.getproperty(vs::$esc_T, fieldname::Symbol)= getfield(vs, :fieldtable)[fieldname].value
        Base.getproperty(vs::$esc_T, fieldname::String)= getfield(vs, :fieldtable)[fieldname].value

        function Base.setproperty!(vs::$esc_T, fieldname::Symbol, value)
            fieldtable = getfield(vs, :fieldtable)

            if haskey(fieldtable, fieldname)
                # Type conversion - checks for the type errors too
                value_converted = convert(fieldtable[fieldname].type, value)

                # Setting the fields
                fieldtable[fieldname].value = value_converted
            else
                # new field without type
                fieldtable[fieldname] = Props(value, typeof(value))
            end

            return vs
        end
        function Base.setproperty!(vs::$esc_T, fieldname::String, value)

            if haskey(fieldtable, fieldname)
                # Type conversion - checks for the type errors too
                value_converted = convert(fieldtable[fieldname].type, value)

                # Setting the fields
                fieldtable[fieldname].value = value_converted
            else
                # new field without type
                fieldtable[fieldname] = Props(value, typeof(value))
            end

            return vs
        end

        Base.fieldnames(vs::$esc_T) = collect(keys(getfield(vs, :fieldtable)))

    end # end quote
end

function get_struct_constructor(modul, T, T_declaration, shared_var = false)
    T_declaration = QuoteNode(T_declaration)
    esc_T = esc(T)
    return quote
        # kw method
        function $esc_T(; args...)
            # Type Checking for already defined fields
            fieldtableDeclaration = getfield(getfield($modul, $T_declaration), :fieldtable)

            # shared instance
            @static if ($shared_var)
                fieldtable = fieldtableDeclaration
            else
                fieldtable = deepcopy(fieldtableDeclaration)
            end

            for (fieldname, value) in args
                if haskey(fieldtable, fieldname)
                    # Type conversion - checks for the type errors too
                    value_converted = convert(fieldtable[fieldname].type, value)

                    # Setting the fields
                    fieldtable[fieldname].value = value_converted
                else
                    # new field without type
                    fieldtable[fieldname] = Props(value, typeof(value))
                end
            end
            return $esc_T(fieldtable)
        end
    end
end

function initialize_struct(modul, T, esc_args_field, args_defaultvalue, args_type, T_declaration)
    esc_T_declaration = esc(T_declaration)

    # we need to store declaration in the memory. This is needed for type checking, conversion, etc.
    # If someday Julia decides to include this in the Base, we can hide this from the memory.

    # also, we can't just call the constrctor directly because some default values might be missing in the declaration. We should use FieldTable

    return quote
        # initialize the struct
        $esc_T_declaration = $modul.$T( VarStructs.FieldTable( n => VarStructs.Props(v, t) for (n, v, t) in zip($esc_args_field, [$(esc.(args_defaultvalue)...)],  [$(esc.(args_type)...)]) ) )
    end
end

function show_struct(esc_T)
    return quote
        function Base.show(io::IO, vs::$esc_T)
            fieldtable = getfield(vs, :fieldtable)
            # arg_show = Vector{String}(undef, length(fieldtable))
            # iArg = 1
            arg_show_all = ""
            for (name, props) in fieldtable
                arg_show_all = arg_show_all * "    $name::$(props.type) = $(props.value), \n"
                # iArg = iArg+1
            end

            out = "$($esc_T)(\n$arg_show_all)"
            return print(io, out)
        end
    end
end


text:
# TODO: modularize code, write simulate_unit_cell independent of the geometry, all geometry related code should be in geometry_code

"""
```x, y, Ez, dpml, dsource, resolution = simulate_unit_cell(ps; refractive_indexes=zeros(3), frequency=1, interstice = 0.5, hole = 0.75, Ly = 17)```
This function simulates Helmholtz equation for a unit-cell consisting of holes in a substrate."""
function simulation_hole_layers_unit_cell(ps; refractive_indexes=zeros(3), frequency=1, interstice = 0.5, hole = 0.75, Ly = 17)
    #parameters for Maxwell solver
    # specific to the unit-cell
    Lx = 0.95
    # Ly = 17
    ω = 2pi*frequency
    dpml = 2
    resolution = 40
    dsource = 1

    # specific to the unit-cell
    if refractive_indexes == zeros(3)
        refractive_indexes = Float64[1.0, 1.0, 1.45]
    end

    # ϵ_hole_layers is specific to the unit cell, all the rest is not
    A, nx, ny, x, y = Maxwell_2d(Lx, Ly,
    (x,y)-> ϵ_hole_layers(x, y, ps, refractive_indexes=refractive_indexes, interstice = interstice, hole = hole),
    ω, dpml, resolution)

    J = zeros(ComplexF64, (ny, nx))
    J[end-(dpml + dsource) * resolution, :]  = 1im  * ω * ones(nx)*resolution
    Ez = reshape(A \ J[:], (ny, nx))

    return x, y, Ez, dpml, dsource, resolution
end

"""
```function simulation_pillar_unit_cell(p; kwargs=Dict())```
This function simulated Helmholtz equation for a unit-cell consisting of a pillar on top of a substrate.

@ parameters:

- p: width of the pillar in adimensional units
- kwargs: dictionary with parameters of the unit cell we would want to change

@return (everything needed to compute the local field)

- x, y: coordinates
- Ez: values of the electric field out of plane
- dpml
- dsource is the distance of the source from the pml (in adimensional unit), also use for the distance between the monitor and the pml
- resolution
"""
function simulation_pillar_unit_cell(p; kwargs=Dict())
    @assert length(p) == 1

    # parameter unit cell
    height= get(kwargs, "height", 600/443)

    #parameters for Maxwell solver
    refractive_indexes= get(kwargs, "refractive_indexes", [1., 2., 1.4])
    frequency= get(kwargs, "frequency", 1.)
    Ly = get(kwargs, "Ly", 7)
    resolution = get(kwargs, "resolution", 40)
    dpml = get(kwargs, "dpml", 2)
    Lx = get(kwargs, "Lx", 1)
    dsource = get(kwargs, "dsource", 1)

    ω = 2pi*frequency

    # specific to the unit-cell
    A, nx, ny, x, y = Maxwell_2d(Lx, Ly,
    (x,y)-> ϵ_pillar_function(x, y, [p, height], refractive_indexes=refractive_indexes),
    ω, dpml, resolution)

    J = zeros(ComplexF64, (ny, nx))
    J[end-(dpml + dsource) * resolution, :]  = 1im  * ω * ones(nx)*resolution
    Ez = reshape(A \ J[:], (ny, nx))

    return x, y, Ez, dpml, dsource, resolution
end

"""
```x, y, zeroth_order = get_local_field(ps; refractive_indexes=zeros(3), frequency=1, interstice = 0.5, hole = 0.75, Ly = 17)```
    this function simulates the structure with holes using Maxwell's equations,
    and computes the zeroth-order Fourier coefficient.

input arguments:

- ps: array with parametrization of the unit-cell corresponding the widths of holes in the substrate
- refractive_indexes: optional argument with refractive indexes of background, hole and substrate. For reference simulation: set refractive indexes to ones(3)*eps_substrate
- frequency (optional): frequency of simulation
- interstice (optional): number of wavelength between holes
- hole (optional): number of wavelength inside holes
- Ly (optional): length of the computational domain

returns:

- return the complex zeroth order Fourier coefficient of the transmitted field
"""
function get_local_field(ps; simulate_unit_cell=simulation_hole_layers_unit_cell, refractive_indexes=zeros(3), frequency=1, interstice = 0.5, hole = 0.75, Ly = 17)
    x, y, Ez, dpml, dsource, resolution = simulate_unit_cell(ps,
    refractive_indexes=refractive_indexes, frequency=frequency,
    interstice = interstice, hole = hole, Ly = Ly)
    return x, y, mean(Ez[(dpml + dsource) * resolution, :])
end

function get_local_field(ps, simulate_unit_cell; kwargs=Dict())
    x, y, Ez, dpml, dsource, resolution = simulate_unit_cell(ps, kwargs=kwargs)
    return x, y, mean(Ez[(dpml + dsource) * resolution, :])
end

ext: jl


text:
using GPMLj, Test, Random

include("kernels.jl")
include("likelihoods.jl")
include("mean_functions.jl")
include("models.jl")
include("parameter_priors.jl")
include("train.jl")


text:
using BenchmarkTools
using MutualInformationImageRegistration
using MutualInformationImageRegistration: FastHistograms.SingleThreadFixedWidth2DHistogram

function bench_mi()
    fixed = rand(UInt8, 70, 100)
    moving = rand(UInt8, 70, 100)

    mi = MutualInformationContainer(SingleThreadFixedWidth2DHistogram())
    bm = @benchmark mutual_information!($mi, $fixed, $moving)
    display(bm)
    println()

    # Current benchmark:
    # Minimum time should be ~32 μs
    #
    # BenchmarkTools.Trial:
    #   memory estimate:  8.75 KiB
    #   allocs estimate:  5
    #   --------------
    #   minimum time:     32.375 μs (0.00% GC)
    #   median time:      34.322 μs (0.00% GC)
    #   mean time:        37.622 μs (0.36% GC)
    #   maximum time:     1.386 ms (96.61% GC)
    #   --------------
    #   samples:          10000
    #   evals/sample:     1
end

function bench_register()
    MAX_SHIFT = 11
    padding = [-10, -10, 10, 10]
    mi = MutualInformationContainer(SingleThreadFixedWidth2DHistogram())
    full_image = rand(UInt8, 500, 300)
    # Create an asymmetric pattern of black pixels to register against
    view(full_image, 300:320, 200:210) .= 0xff

    fixed = full_image[(300-10):(330+10), (200-10):(220+10)]

    buffer = Array{UInt8}(undef, (size(fixed) .+ (MAX_SHIFT * 2))...)

    expected_x = 5
    expected_y = -5
    moving =
        [300, 200, 330, 220] .+ padding .+ [expected_x, expected_y, expected_x, expected_y]

    bm = @benchmark register!(
        $mi,
        $full_image,
        $fixed,
        $moving,
        $MAX_SHIFT,
        $MAX_SHIFT,
        $buffer,
    )
    display(bm)
    println()

    # Current benchmark:
    # Minimum time should be ~900 μs
    #
    # BenchmarkTools.Trial:
    #   memory estimate:  558.05 KiB
    #   allocs estimate:  332
    #   --------------
    #   minimum time:     1.122 ms (0.00% GC)
    #   median time:      1.229 ms (0.00% GC)
    #   mean time:        1.289 ms (1.36% GC)
    #   maximum time:     4.280 ms (59.71% GC)
    #   --------------
    #   samples:          3877
    #   evals/sample:     1
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMSparse.jl/blob/master/LICENSE

import Base: size, resize!, convert, getindex, get!, isempty, empty!, sparsevec, full
import Base.SparseArrays: nonzeroinds, nonzeros, nnz

global const SPARSEVECTORCOO_DEFAULT_BLOCK_SIZE = 1024

type SparseVectorCOO{Tv,Ti<:Integer} <: AbstractSparseArray{Tv,Ti<:Integer,1}
    I :: Vector{Ti}
    V :: Vector{Tv}
    cnt :: Int64
    max_cnt :: Int64
    blk_size :: Int64
end

function SparseVectorCOO()
    return SparseVectorCOO(Int64[], Float64[], 0, 0, SPARSEVECTORCOO_DEFAULT_BLOCK_SIZE)
end

function SparseVectorCOO{Tv,Ti<:Integer}(I::Vector{Ti}, V::Vector{Tv})
    @assert length(I) == length(V)
    return SparseVectorCOO(I, V, length(V), length(V), SPARSEVECTORCOO_DEFAULT_BLOCK_SIZE)
end

function SparseVectorCOO{T}(b::Vector{T})
    return SparseVectorCOO(sparsevec(b))
end

function SparseVectorCOO{Tv,Ti<:Integer}(b::SparseVector{Tv,Ti})
    I, V = findnz(b)
    return SparseVectorCOO(I, V)
end

function convert{Tv,Ti<:Integer}(::Type{SparseVectorCOO{Tv,Ti}}, b::SparseVector{Tv,Ti})
    I, V = findnz(b)
    return SparseVectorCOO(I, V)
end

function size(b::SparseVectorCOO)
    if b.cnt == 0
        return (0,)
    end
    return (maximum(b.I),)
end

function sparsevec(b::SparseVectorCOO, args...)
    b_new = Base.sparsevec(b.I[1:b.cnt], b.V[1:b.cnt], args...)
    SparseArrays.dropzeros!(b_new)
    return b_new
end

function empty!(b::SparseVectorCOO)
    empty!(b.I)
    empty!(b.V)
    b.cnt = 0
    b.max_cnt = 0
    return b
end

function nnz(b::SparseVectorCOO)
    return nnz(sparsevec(b))
end

function nonzeroinds(b::SparseVectorCOO)
    return nonzeroinds(sparsevec(b))
end

function nonzeros(b::SparseVectorCOO)
    return nonzeros(sparsevec(b))
end

function isempty(b::SparseVectorCOO)
    return b.cnt == 0
end

function full(b::SparseVectorCOO, args...)
    return full(sparsevec(b, args...))
end

#=

function get!{Tv,Ti<:Integer}(b::SparseVectorCOO{Tv,Ti}, dofs::Vector{Ti}, data::Vector{Tv})
    @assert length(dofs) == length(data)
    for i=1:length(data)
        data[i] = get(b, i)
    end
    return nothing
end

function get{Tv,Ti<:Integer}(b::SparseVectorCOO{Tv,Ti}, idx::Ti)
    value = 0.0
    for i=1:b.cnt
        if b.I[i] == idx
            value += b.V[i]
        end
    end
    return value
end

=#

""" Allocate memory for I,V vectors of SparseVectorCOO. """
function allocate!(b::SparseVectorCOO, N::Int)
    if N < b.max_cnt
        warn("Down-sizing SparseVectorCOO vectors")
    end
    b.max_cnt = N
    Base.resize!(b.I, b.max_cnt)
    Base.resize!(b.V, b.max_cnt)
    return b
end

""" Add new value to sparse matrix. """
function add!{T}(b::SparseVectorCOO{T}, I::Int, V::T)
    if b.cnt == b.max_cnt
        allocate!(b, b.cnt+b.blk_size)
    end
    b.cnt += 1
    @inbounds begin
        b.I[b.cnt] = I
        b.V[b.cnt] = V
    end
    return b
end

""" Add new data to COO Sparse vector. """
function add!(A::SparseVectorCOO, dofs::Vector{Int}, data::Array{Float64})
    @assert length(dofs) == length(data)
    for i=1:length(data)
        @inbounds add!(A, dofs[i], data[i])
    end
    return nothing
end

""" Add SparseVector to SparseVectorCOO. """
function add!(a::SparseVectorCOO, b::SparseVector)
    I, V = findnz(b)
    for (i, v) in zip(I, V)
        add!(a, i, v)
    end
    return nothing
end

""" Add one or more SparseVectorCOO to SparseVectorCOO. """
function add!(a::SparseVectorCOO, rest::SparseVectorCOO...)
    for b in rest
        for (i, v) in zip(b.I[1:b.cnt], b.V[1:b.cnt])
            add!(a, i, v)
        end
    end
    return nothing
end


text:
@inline function issquare(m::AbstractArray{T,2}) where {T<:Number}
    r, c = size(m)
    return r === c
end

@inline function issquare(m::Array{DoubleFloat{T},2}) where {T<:IEEEFloat}
    r, c = size(m)
    return r === c
end

function norm(v::Array{DoubleFloat{T},N}, p::Real=2.0) where {N, T<:IEEEFloat}
    if isinf(p)
        return signbit(p) ? minimum(abs.(v)) : maximum(abs.(v))
    elseif p==2
        return vp = sqrt(sum(v .* v))
    else
        vp = sum((v).^(p))
        r = inv(DoubleFloat{T}(p))
        return vp^r
    end
end

function norm(v::Array{Complex{DoubleFloat{T}},N}, p::Real=2.0) where {N, T<:IEEEFloat}
    if isinf(p)
        return signbit(p) ? minimum(abs.(real.(v))) : maximum(abs.(real.(v)))
    elseif p==2
        return vp = sqrt(sum(conj.(v) .* v))
    else
        vp = sum(abs.(v).^(p))
        r = inv(DoubleFloat{T}(p))
        return vp^r
    end
end

function LinearAlgebra.normalize(v::Array{DoubleFloat{T},N}, p::Real=2.0) where {N, T<:IEEEFloat}
    return v ./ norm(v, p)
end

function LinearAlgebra.normalize(v::Array{Complex{DoubleFloat{T}},N}, p::Real=2.0) where {N, T<:IEEEFloat}
    return v ./ norm(v, p)
end


text:
module MITgcmTools

using Dates, Printf, SparseArrays, Artifacts, LazyArtifacts, UUIDs, Suppressor
using OrderedCollections, DataFrames, NetCDF, MeshArrays, ClimateModels

include("Types.jl")
include("ReadFiles.jl")
include("ModelSteps.jl")
include("FormatConversions.jl")
include("PhysicalOceanography.jl")
include("verification_experiments.jl")

export MITgcm_path, MITgcm_download
export PICKUP_hs94_path, PICKUP_hs94_download
export MITgcm_config, MITgcm_namelist, MITgcm_launch
export testreport, build, compile, setup, clean
#export pause, stop, clock, monitor, train, help
export verification_experiments, read_namelist, write_namelist
export read_mdsio, read_meta, read_available_diagnostics
export scan_rundir, scan_stdout
export read_bin, read_flt, read_mnc, read_nctiles, findtiles
export GridLoad_mnc, GridLoad_mdsio
export cube2compact, compact2cube, convert2array, convert2gcmfaces
export SeaWaterDensity, MixedLayerDepth

p=dirname(pathof(MITgcmTools))
artifact_toml = joinpath(p, "../Artifacts.toml")
MITgcm_hash = artifact_hash("MITgcm", artifact_toml)

"""
    MITgcm_path

Path to a MITgcm folder. `MITgcm_path[1]` should generally be used. `MITgcm_path[2]` is mostly
meant to facilitate comparisons between e.g. MITgcm releases when needed.
"""
MITgcm_path = [ joinpath(artifact_path(MITgcm_hash)*"/","MITgcm_test/"),
                joinpath(artifact_path(MITgcm_hash)*"/","MITgcm_test/")]

PICKUP_hs94_path = artifact_path(artifact_hash("PICKUP_hs94", artifact_toml))

MITgcm_download() = artifact"MITgcm"
PICKUP_hs94_download() = artifact"PICKUP_hs94"

#more:
#
#using Plots; include("recipes_plots.jl"); export qwckplot

#deprecated:
#
#include("deprecated.jl"); export prep_MTRX, read_SPM

end # module


text:
@doc raw"""
Axisymmetric tokamak equilibrium in (r,θ,ϕ) coordinates with covariant
components of the vector potential given by
```math
A (r, \theta, \phi) = B_0 \, \bigg( 0 , \, \frac{r R_0}{\cos (\theta)} - \bigg( \frac{R_0}{\cos (\theta)} \bigg)^2 \, \ln \bigg( \frac{R}{R_0} \bigg) , \, - \frac{r^2}{2 q_0} \bigg)^T ,
```
resulting in the magnetic field with covariant components
```math
B (r, \theta, \phi) = \frac{B_0}{q_0} \, \bigg( 0 , \, \frac{r^2}{R}, \, q_0 R_0 \bigg)^T ,
```
where $R = R_0 + r \cos \theta$.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `q₀`: safety factor at magnetic axis
"""
struct AxisymmetricTokamakToroidal{T <: Number} <: AnalyticEquilibrium
    name::String
    R₀::T
    B₀::T
    q₀::T

    function AxisymmetricTokamakToroidal{T}(R₀::T, B₀::T, q₀::T) where T <: Number
        new("AxisymmetricTokamakEquilibriumToroidal", R₀, B₀, q₀)
    end
end

AxisymmetricTokamakToroidal(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakToroidal{T}(R₀, B₀, q₀)

function Base.show(io::IO, equ::AxisymmetricTokamakToroidal)
    print(io, "Axisymmetric Tokamak Equilibrium in Toroidal Coordinates with\n")
    print(io, "  R₀ = ", equ.R₀, "\n")
    print(io, "  B₀ = ", equ.B₀, "\n")
    print(io, "  q₀ = ", equ.q₀)
end

function X(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    R(x,equ) * cos(ϕ(x,equ))
end

function Y(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    R(x,equ) * sin(ϕ(x,equ))
end

function Z(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    r(x,equ) * sin(θ(x,equ))
end

function R(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    equ.R₀ + r(x,equ) * cos(θ(x,equ))
end

function R²(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    R(x,equ)^2
end

function r(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    x[1]
end

function r²(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    r(x,equ)^2
end

function θ(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    x[2]
end

function ϕ(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    x[3]
end

function J(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    r(x,equ) * R(x,equ)
end

function periodicity(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    p = zero(x)
    p[2] = 2π
    p[3] = 2π
    return p
end

function A₁(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    zero(T)
end

function A₂(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    equ.B₀ * equ.R₀ / cos(θ(x,equ))^2 * ( r(x,equ) * cos(θ(x,equ)) - equ.R₀ * log(R(x,equ) / equ.R₀) )
end

function A₃(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    - equ.B₀ * r²(x,equ) / equ.q₀ / 2
end

function g₁₁(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    one(T)
end

function g₂₂(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    r²(x, equ)
end

function g₃₃(x::AbstractArray{T,1}, equ::AxisymmetricTokamakToroidal) where {T <: Number}
    R²(x, equ)
end

macro axisymmetric_tokamak_equilibrium_toroidal(R₀, B₀, q₀)
    generate_equilibrium_code(AxisymmetricTokamakToroidal(R₀, B₀, q₀); output=false)
end


text:
@testset "Test ideal lens" begin

    arr = randn((10, 10))
    @test all(abs.(lens(10.0, 12.0, size(arr), radius=1000)) .≈ 1)

    arr2 = zeros(size(arr))
    arr2[6, 6] = 1
    @test arr2 ≈ abs.(lens(10.0, 12.0, size(arr), radius=0))
end

@testset "Test microlens" begin
    mla = micro_lens_array(1000e-3, 10e-3, 10e-3, (64, 64))
    l = lens(1000e-3, 10e-3, (64, 64), radius=10e-3/2)
    @test mla ≈ l

    mla = micro_lens_array(1000e-3, 10e-3, 10e-3, (64, 64), aperture=false)
    l = lens(1000e-3, 10e-3, (64, 64), radius=Inf)
    @test mla ≈ l

    mla = micro_lens_array(1000e-3, 20e-3, 20e-3 / 2, (128, 128), aperture=false, centered=false)
    l = lens(1000e-3, 10e-3, (64, 64), radius=Inf)
    @test mla[1:64, 1:64] ≈ l

    mla = micro_lens_array(1000e-3, 20e-3, 20e-3 / 2, (128, 128), aperture=false, centered=true)
    l = lens(1000e-3, 10e-3, (64, 64), radius=Inf)
    @test center_extract(mla, (64, 64)) ≈ l

end

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
using Test
using LinearAlgebra

import ClimaCore:
    ClimaCore,
    slab,
    Spaces,
    Domains,
    Meshes,
    Geometry,
    Topologies,
    Spaces,
    Fields,
    Operators
import ClimaCore.Utilities: half

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

const n_vert = 10
const n_horz = 4
const p_horz = 4

const R = 6.4e6 # radius
const Ω = 7.2921e-5 # Earth rotation (radians / sec)
const z_top = 3.0e4 # height position of the model top
const grav = 9.8 # gravitational constant
const p_0 = 1e5 # mean sea level pressure

const R_d = 287.058 # R dry (gas constant / mol mass dry air)
const T_tri = 273.16 # triple point temperature
const γ = 1.4 # heat capacity ratio
const cv_d = R_d / (γ - 1)
const cp_d = R_d + cv_d

const T_0 = 300 # isothermal atmospheric temperature
const H = R_d * T_0 / grav # scale height

# set up function space
function sphere_3D(
    R = 6.4e6,
    zlim = (0, 30.0e3),
    helem = 4,
    zelem = 15,
    npoly = 5,
)
    FT = Float64
    vertdomain = Domains.IntervalDomain(
        Geometry.ZPoint{FT}(zlim[1]),
        Geometry.ZPoint{FT}(zlim[2]);
        boundary_names = (:bottom, :top),
    )
    vertmesh = Meshes.IntervalMesh(vertdomain, nelems = zelem)
    vert_center_space = Spaces.CenterFiniteDifferenceSpace(vertmesh)

    horzdomain = Domains.SphereDomain(R)
    horzmesh = Meshes.EquiangularCubedSphere(horzdomain, helem)
    horztopology = Topologies.Topology2D(horzmesh)
    quad = Spaces.Quadratures.GLL{npoly + 1}()
    horzspace = Spaces.SpectralElementSpace2D(horztopology, quad)

    hv_center_space =
        Spaces.ExtrudedFiniteDifferenceSpace(horzspace, vert_center_space)
    hv_face_space = Spaces.FaceExtrudedFiniteDifferenceSpace(hv_center_space)
    return (hv_center_space, hv_face_space)
end

# geopotential
Φ(z) = grav * z

# pressure
function pressure(ρ, e, normuvw, z)
    I = e - Φ(z) - normuvw^2 / 2
    T = I / cv_d + T_tri
    return ρ * R_d * T
end

# initial conditions for density and total energy
function init_sbr_thermo(z)

    p = p_0 * exp(-z / H)
    ρ = 1 / R_d / T_0 * p

    e = cv_d * (T_0 - T_tri) + Φ(z)
    ρe = ρ * e

    return (ρ = ρ, ρe = ρe)
end

function rhs!(dY, Y, _, t)
    cρ = Y.Yc.ρ # density on centers
    fw = Y.w # Covariant3Vector on faces
    cuₕ = Y.uₕ # Covariant12Vector on centers
    cρe = Y.Yc.ρe # total energy on centers

    dρ = dY.Yc.ρ
    dw = dY.w
    duₕ = dY.uₕ
    dρe = dY.Yc.ρe
    z = c_coords.z

    # # 0) update w at the bottom
    # fw = -g^31 cuₕ/ g^33 ????????

    hdiv = Operators.Divergence()
    hgrad = Operators.Gradient()
    hcurl = Operators.Curl()

    dρ .= 0 .* cρ
    dw .= 0 .* fw
    duₕ .= 0 .* cuₕ
    dρe .= 0 .* cρe

    # hyperdiffusion not needed in SBR

    # 1) Mass conservation
    If2c = Operators.InterpolateF2C()
    Ic2f = Operators.InterpolateC2F(
        bottom = Operators.Extrapolate(),
        top = Operators.Extrapolate(),
    )
    cw = If2c.(fw)
    cuvw = Geometry.Covariant123Vector.(cuₕ) .+ Geometry.Covariant123Vector.(cw)

    # 1.a) horizontal divergence
    dρ .-= hdiv.(cρ .* (cuvw))

    # 1.b) vertical divergence
    vdivf2c = Operators.DivergenceF2C(
        top = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
        bottom = Operators.SetValue(Geometry.Contravariant3Vector(0.0)),
    )
    # we want the total u³ at the boundary to be zero: we can either constrain
    # both to be zero, or allow one to be non-zero and set the other to be its
    # negation

    # explicit part
    dρ .-= vdivf2c.(Ic2f.(cρ .* cuₕ))
    # implicit part
    dρ .-= vdivf2c.(Ic2f.(cρ) .* fw)

    # 2) Momentum equation

    # curl term
    # effectively a homogeneous Dirichlet condition on u₁ at the boundary
    vcurlc2f = Operators.CurlC2F(
        bottom = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
        top = Operators.SetCurl(Geometry.Contravariant12Vector(0.0, 0.0)),
    )
    cω³ = hcurl.(cuₕ) # Contravariant3Vector
    fω¹² = hcurl.(fw) # Contravariant12Vector
    fω¹² .+= vcurlc2f.(cuₕ) # Contravariant12Vector

    # cross product
    # convert to contravariant
    # these will need to be modified with topography
    fu¹² =
        Geometry.Contravariant12Vector.(
            Geometry.Covariant123Vector.(Ic2f.(cuₕ)),
        ) # Contravariant12Vector in 3D
    fu³ = Geometry.Contravariant3Vector.(Geometry.Covariant123Vector.(fw))
    @. dw -= fω¹² × fu¹² # Covariant3Vector on faces
    @. duₕ -= If2c(fω¹² × fu³)

    # Needed for 3D:
    @. duₕ -=
        (f + cω³) ×
        Geometry.Contravariant12Vector(Geometry.Covariant123Vector(cuₕ))

    ce = @. cρe / cρ
    cp = @. pressure(cρ, ce, norm(cuvw), z)

    @. duₕ -= hgrad(cp) / cρ
    vgradc2f = Operators.GradientC2F(
        bottom = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
        top = Operators.SetGradient(Geometry.Covariant3Vector(0.0)),
    )
    @. dw -= vgradc2f(cp) / Ic2f(cρ)

    cE = @. (norm(cuvw)^2) / 2 + Φ(z)
    @. duₕ -= hgrad(cE)
    @. dw -= vgradc2f(cE)

    # 3) total energy

    @. dρe -= hdiv(cuvw * (cρe + cp))
    @. dρe -= vdivf2c(fw * Ic2f(cρe + cp))
    @. dρe -= vdivf2c(Ic2f(cuₕ * (cρe + cp)))

    Spaces.weighted_dss!(dY.Yc)
    Spaces.weighted_dss!(dY.uₕ)
    Spaces.weighted_dss!(dY.w)

    return dY

end

# set up 3D spherical domain and coords
hv_center_space, hv_face_space =
    sphere_3D(R, (0, z_top), n_horz, n_vert, p_horz)
c_coords = Fields.coordinate_field(hv_center_space)
f_coords = Fields.coordinate_field(hv_face_space)

# Coriolis
const f = @. Geometry.Contravariant3Vector(
    Geometry.WVector(2 * Ω * sind(c_coords.lat)),
)

# discrete hydrostatic profile
zc_vec = parent(c_coords.z) |> unique

N = length(zc_vec)
ρ = zeros(Float64, N)
p = zeros(Float64, N)
ρe = zeros(Float64, N)

# compute ρ, ρe, p from analytical formula; not discretely balanced
for i in 1:N
    var = init_sbr_thermo(zc_vec[i])
    ρ[i] = var.ρ
    ρe[i] = var.ρe
    p[i] = pressure(ρ[i], ρe[i] / ρ[i], 0.0, zc_vec[i])
end

ρ_ana = copy(ρ) # keep a copy for analytical ρ which will be used in correction ρe

function discrete_hydrostatic_balance!(ρ, p, dz, grav)
    for i in 1:(length(ρ) - 1)
        ρ[i + 1] = -ρ[i] - 2 * (p[i + 1] - p[i]) / dz / grav
    end
end

discrete_hydrostatic_balance!(ρ, p, z_top / n_vert, grav)
# now ρ (after correction) and p (computed from analytical relation) are in discrete hydrostatic balance
# only need to correct ρe without changing ρ and p, i.e., keep ρT unchanged before vs after the correction on ρ
ρe = @. ρe + (ρ - ρ_ana) * Φ(zc_vec) - (ρ - ρ_ana) * cv_d * T_tri

# Note: In princile, ρe = @. cv_d * p /R_d - ρ * cv_d * T_tri + ρ * Φ(zc_vec) should work,
#       however, it is not as accurate as the above correction

# set up initial condition: not discretely balanced; only create a Field as a place holder
Yc = map(coord -> init_sbr_thermo(coord.z), c_coords)
# put the dicretely balanced ρ and ρe into Yc
parent(Yc.ρ) .= ρ  # Yc.ρ is a VIJFH layout
parent(Yc.ρe) .= ρe

# initialize velocity: at rest
uₕ = map(_ -> Geometry.Covariant12Vector(0.0, 0.0), c_coords)
w = map(_ -> Geometry.Covariant3Vector(0.0), f_coords)
Y = Fields.FieldVector(Yc = Yc, uₕ = uₕ, w = w)

# initialize tendency
dYdt = similar(Y)
# set up rhs
rhs!(dYdt, Y, nothing, 0.0)

# run!
using OrdinaryDiffEq
# Solve the ODE
T = 3600
dt = 5
prob = ODEProblem(rhs!, Y, (0.0, T))

haskey(ENV, "CI_PERF_SKIP_RUN") && exit() # for performance analysis

# solve ode
sol = solve(
    prob,
    SSPRK33(),
    dt = dt,
    saveat = dt,
    progress = true,
    adaptive = false,
    progress_message = (dt, u, p, t) -> t,
)

uₕ_phy = Geometry.transform.(Ref(Geometry.UVAxis()), sol.u[end].uₕ)
w_phy = Geometry.transform.(Ref(Geometry.WAxis()), sol.u[end].w)

@test maximum(abs.(uₕ_phy.components.data.:1)) ≤ 1e-11
@test maximum(abs.(uₕ_phy.components.data.:2)) ≤ 1e-11
@test maximum(abs.(w_phy |> parent)) ≤ 1e-11
@test norm(sol.u[end].Yc.ρ) ≈ norm(sol.u[1].Yc.ρ) rtol = 1e-2
@test norm(sol.u[end].Yc.ρe) ≈ norm(sol.u[1].Yc.ρe) rtol = 1e-2


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    data = readdlm(fname, String)
    return data
end

# useful constants
const openers = ['(', '[', '{', '<']
const closers = Dict([(')','('), (']','['), ('}','{'), ('>','<')])
const errorScores  = Dict([(')',3), (']',57), ('}',1197), ('>',25137)])
const incomScores  = Dict([('(',1), ('[',2),  ('{',3),    ('<',4)])

function Part1(data)
    score = 0
    for line in data
        stack = []
        for c in line
            if c in openers
                pushfirst!(stack, c)
            else
                o = popfirst!(stack)
                if closers[c] != o
                    score += errorScores[c]
                    break
                end
            end
        end
    end

    return score
end

function Part2(data)
    scores = []
    for line in data
        stack = []
        skipline = false
        for c in line
            if c in openers
                pushfirst!(stack, c)
            else
                o = popfirst!(stack)
                if closers[c] != o
                    skipline = true
                    break
                end
            end
        end

        if !skipline && length(stack) > 0
            score = 0
            while length(stack) > 0
                o = popfirst!(stack)
                score = 5*score + incomScores[o]
            end
            push!(scores, score)
        end
    end

    sort!(scores)
    return scores[div(length(scores),2) + 1]
end

#data = ReadData("../data/test10.txt")
data = ReadData("../data/day10.txt")

println(Part1(data))
println(Part2(data))


text:
################
# Content Type #
################

type Content <: GitLabType
    file_name::Nullable{GitLabString}
    file_path::Nullable{GitLabString}
    size::Nullable{Int}
    encoding::Nullable{GitLabString}
    content::Nullable{GitLabString}
    ref::Nullable{GitLabString}
    blob_id::Nullable{GitLabString}
    commit_id::Nullable{GitLabString}
    last_commit_id::Nullable{GitLabString}

#=
    typ::Nullable{GitLabString}
    name::Nullable{GitLabString}
    target::Nullable{GitLabString}
    url::Nullable{HttpCommon.URI}
    git_url::Nullable{HttpCommon.URI}
    html_url::Nullable{HttpCommon.URI}
    download_url::Nullable{HttpCommon.URI}
=#
end

Content(data::Dict) = json2gitlab(Content, data)
Content(file_path::AbstractString) = Content(Dict("file_path" => file_path))

namefield(content::Content) = content.file_path

###############
# API Methods #
###############

function file(repo, path, ref; options...)
    result = gh_get_json(content_uri(repo, path, ref); options...)
    return Content(result)
end

#= TODO No equivalent API
function directory(repo, path; options...)
    results, page_data = gh_get_paged_json(content_uri(repo, path); options...)
    return map(Content, results), page_data
end
=#

function create_file(repo, path; options...)
    result = gh_put_json(content_uri(repo, path); options...)
    return build_content_response(result)
end

function update_file(repo, path; options...)
    result = gh_put_json(content_uri(repo, path); options...)
    return build_content_response(result)
end

function delete_file(repo, path; options...)
    result = gh_delete_json(content_uri(repo, path); options...)
    return build_content_response(result)
end

function readme(repo; options...)
    ## result = gh_get_json("/api/v3/projects/$(get(repo.id))/readme"; options...)
    result = gh_get_json(content_uri(repo, "README.md"); options...)
    return Content(result)
end

function permalink(content::Content, commit)
    url = string(get(content.html_url))
    prefix = get(content.typ) == "file" ? "blob" : "tree"
    rgx = Regex(string("\/", prefix, "\/.*?\/"))
    replacement = string("/", prefix, "/", name(commit), "/")
    return HttpCommon.URI(replace(url, rgx, replacement))
end

###########################
# Content Utility Methods #
###########################

## content_uri(repo, path) = "/api/v3/projects/$(get(repo.id))/contents/$(name(path))"
## content_uri(repo, path) = "/api/v3/projects/$(get(repo.id))/files"
content_uri(repo, path, ref) = "/api/v3/projects/$(get(repo.id))/repository/files?file_path=$(name(path))&ref=$(name(ref))"
content_uri(repo, path) = "/api/v3/projects/$(get(repo.id))/repository/files?file_path=$(name(path))&ref=master"

function build_content_response(json::Dict)
    results = Dict()
    haskey(json, "commit") && setindex!(results, Commit(json["commit"]), "commit")
    haskey(json, "content") && setindex!(results, Content(json["content"]), "content")
    return results
end


text:
struct DBReference
    database::SQLite.DB
    pkey::Int
    campaign::Int
    spectrum::DBSpectrum
    elements::Vector{Element} # The elements evident in the reference spectrum
end

Base.show(io::IO, dbr::DBReference) =
    print(io, "Reference[$(dbr.spectrum.name) with $(join(symbol.(dbr.elements),','))]")

NeXLUncertainties.asa(::Type{Spectrum}, dbr::DBReference) = NeXLUncertainties.asa(Spectrum, dbr.spectrum)

function Base.write(db::SQLite.DB, ::Type{DBReference}, campaign::Int, spectrum::Int, elements::Vector{Element})::Int
    stmt1 = SQLite.Stmt(db, "INSERT INTO REFERENCESPECTRUM(CAMPAIGN, SPECTRUM, ELEMENTS) VALUES ( ?, ?, ? );")
    q = DBInterface.execute(stmt1, (campaign, spectrum, _elmstostr(elements)))
    return DBInterface.lastrowid(q)
end

function NeXLUncertainties.asa(::Type{DataFrame}, refs::Vector{DBReference})
    fm, sm, cm = Union{Float64, Missing}, Union{String, Missing}, Union{Material,Missing}
    elms, spname, e0, lt, pc, comp = String[], String[], fm[],  fm[], fm[], cm[]
    for ref in refs
        spec = asa(Spectrum, ref)
        push!(elms, _elmstostr(ref.elements))
        push!(spname, spec[:Name])
        push!(e0, get(spec, :BeamEnergy, missing))
        push!(lt, get(spec, :LiveTime, missing))
        push!(pc, get(spec, :ProbeCurrent, missing))
        push!(comp, get(spec, :Composition, missing))
    end
    return DataFrame(Elements=elms, Name=spname, BeamEnergy=e0, LiveTime=lt, ProbeCurrent=pc, Composition=comp)
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: marketplace_metering
using AWS.Compat
using AWS.UUIDs

"""
    BatchMeterUsage()

BatchMeterUsage is called from a SaaS application listed on the AWS Marketplace to post metering records for a set of customers. For identical requests, the API is idempotent; requests can be retried with the same records or a subset of the input records. Every request to BatchMeterUsage is for one product. If you need to meter usage for multiple products, you must make multiple calls to BatchMeterUsage. BatchMeterUsage can process up to 25 UsageRecords at a time. A UsageRecord can optionally include multiple usage allocations, to provide customers with usagedata split into buckets by tags that you define (or allow the customer to define). BatchMeterUsage requests must be less than 1MB in size.

# Required Parameters
- `ProductCode`: Product code is used to uniquely identify a product in AWS Marketplace. The product code should be the same as the one used during the publishing of a new product.
- `UsageRecords`: The set of UsageRecords to submit. BatchMeterUsage accepts up to 25 UsageRecords at a time.

"""
batch_meter_usage(ProductCode, UsageRecords; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("BatchMeterUsage", Dict{String, Any}("ProductCode"=>ProductCode, "UsageRecords"=>UsageRecords); aws_config=aws_config)
batch_meter_usage(ProductCode, UsageRecords, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("BatchMeterUsage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ProductCode"=>ProductCode, "UsageRecords"=>UsageRecords), args)); aws_config=aws_config)

"""
    MeterUsage()

API to emit metering records. For identical requests, the API is idempotent. It simply returns the metering record ID. MeterUsage is authenticated on the buyer's AWS account using credentials from the EC2 instance, ECS task, or EKS pod. MeterUsage can optionally include multiple usage allocations, to provide customers with usage data split into buckets by tags that you define (or allow the customer to define).

# Required Parameters
- `ProductCode`: Product code is used to uniquely identify a product in AWS Marketplace. The product code should be the same as the one used during the publishing of a new product.
- `Timestamp`: Timestamp, in UTC, for which the usage is being reported. Your application can meter usage for up to one hour in the past. Make sure the timestamp value is not before the start of the software usage.
- `UsageDimension`: It will be one of the fcp dimension name provided during the publishing of the product.

# Optional Parameters
- `DryRun`: Checks whether you have the permissions required for the action, but does not make the request. If you have the permissions, the request returns DryRunOperation; otherwise, it returns UnauthorizedException. Defaults to false if not specified.
- `UsageAllocations`: The set of UsageAllocations to submit. The sum of all UsageAllocation quantities must equal the UsageQuantity of the MeterUsage request, and each UsageAllocation must have a unique set of tags (include no tags).
- `UsageQuantity`: Consumption value for the hour. Defaults to 0 if not specified.
"""
meter_usage(ProductCode, Timestamp, UsageDimension; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("MeterUsage", Dict{String, Any}("ProductCode"=>ProductCode, "Timestamp"=>Timestamp, "UsageDimension"=>UsageDimension); aws_config=aws_config)
meter_usage(ProductCode, Timestamp, UsageDimension, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("MeterUsage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ProductCode"=>ProductCode, "Timestamp"=>Timestamp, "UsageDimension"=>UsageDimension), args)); aws_config=aws_config)

"""
    RegisterUsage()

Paid container software products sold through AWS Marketplace must integrate with the AWS Marketplace Metering Service and call the RegisterUsage operation for software entitlement and metering. Free and BYOL products for Amazon ECS or Amazon EKS aren't required to call RegisterUsage, but you may choose to do so if you would like to receive usage data in your seller reports. The sections below explain the behavior of RegisterUsage. RegisterUsage performs two primary functions: metering and entitlement.    Entitlement: RegisterUsage allows you to verify that the customer running your paid software is subscribed to your product on AWS Marketplace, enabling you to guard against unauthorized use. Your container image that integrates with RegisterUsage is only required to guard against unauthorized use at container startup, as such a CustomerNotSubscribedException/PlatformNotSupportedException will only be thrown on the initial call to RegisterUsage. Subsequent calls from the same Amazon ECS task instance (e.g. task-id) or Amazon EKS pod will not throw a CustomerNotSubscribedException, even if the customer unsubscribes while the Amazon ECS task or Amazon EKS pod is still running.    Metering: RegisterUsage meters software use per ECS task, per hour, or per pod for Amazon EKS with usage prorated to the second. A minimum of 1 minute of usage applies to tasks that are short lived. For example, if a customer has a 10 node Amazon ECS or Amazon EKS cluster and a service configured as a Daemon Set, then Amazon ECS or Amazon EKS will launch a task on all 10 cluster nodes and the customer will be charged: (10 * hourly_rate). Metering for software use is automatically handled by the AWS Marketplace Metering Control Plane -- your software is not required to perform any metering specific actions, other than call RegisterUsage once for metering of software use to commence. The AWS Marketplace Metering Control Plane will also continue to bill customers for running ECS tasks and Amazon EKS pods, regardless of the customers subscription state, removing the need for your software to perform entitlement checks at runtime.

# Required Parameters
- `ProductCode`: Product code is used to uniquely identify a product in AWS Marketplace. The product code should be the same as the one used during the publishing of a new product.
- `PublicKeyVersion`: Public Key Version provided by AWS Marketplace

# Optional Parameters
- `Nonce`: (Optional) To scope down the registration to a specific running software instance and guard against replay attacks.
"""
register_usage(ProductCode, PublicKeyVersion; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("RegisterUsage", Dict{String, Any}("ProductCode"=>ProductCode, "PublicKeyVersion"=>PublicKeyVersion); aws_config=aws_config)
register_usage(ProductCode, PublicKeyVersion, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("RegisterUsage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ProductCode"=>ProductCode, "PublicKeyVersion"=>PublicKeyVersion), args)); aws_config=aws_config)

"""
    ResolveCustomer()

ResolveCustomer is called by a SaaS application during the registration process. When a buyer visits your website during the registration process, the buyer submits a registration token through their browser. The registration token is resolved through this API to obtain a CustomerIdentifier and product code.

# Required Parameters
- `RegistrationToken`: When a buyer visits your website during the registration process, the buyer submits a registration token through the browser. The registration token is resolved to obtain a CustomerIdentifier and product code.

"""
resolve_customer(RegistrationToken; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("ResolveCustomer", Dict{String, Any}("RegistrationToken"=>RegistrationToken); aws_config=aws_config)
resolve_customer(RegistrationToken, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = marketplace_metering("ResolveCustomer", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("RegistrationToken"=>RegistrationToken), args)); aws_config=aws_config)


text:
include("Include.jl")

# Estimates model parameters for HL60 model -
function main()

  number_of_subdivisions = 10
  number_of_parameters = 153
  number_of_objectives = 12

  # Setup initial parameter guess -
  initial_parameter_array = zeros(number_of_parameters);

  # Load the data dictionary -
  data_dictionary = DataDictionary(0.0,0.0,0.0);

  # Update data dictionary to match new parameters before calculating obj
  parameter_mapping_array = data_dictionary["parameter_name_mapping_array"]
  for index = 1:length(parameter_mapping_array)
    if index <= data_dictionary["number_of_binding"]
      value = data_dictionary["binding_parameter_dictionary"][parameter_mapping_array[index]];
      initial_parameter_array[index] = value;
    else
      value = data_dictionary["control_parameter_dictionary"][parameter_mapping_array[index]];
      initial_parameter_array[index] = value;
    end
  end

  # Setup the error archives -
  ec_array = zeros(number_of_objectives)
  pc_array = zeros(number_of_parameters)
  for index in collect(1:number_of_subdivisions)

    # Run JuPOETs -
    (EC,PC,RA) = estimate_ensemble(objective_function,neighbor_function,acceptance_probability_function,cooling_function,initial_parameter_array;rank_cutoff=4,maximum_number_of_iterations=10,show_trace=true)

    # Package -
    ec_array = [ec_array EC]
    pc_array = [pc_array PC]

    # Take the *best* value from the current ensemble, refine it and go around again -
    total_error = sum(ec_array[:,2:end],1)

    # Which col is the min error?
    min_index = indmin(total_error)
    @show index,total_error[min_index]

    # dump to disk (in case we need to restart)
    # error cache -
    path_name = "./ensemble/ec_array.dat.$(index)"
    writedlm(path_name,ec_array)

    # parameter cache -
    path_name = "./ensemble/pc_array.dat.$(index)"
    writedlm(path_name,pc_array)

    # Refine the best solution -
    initial_parameter_array = pc_array[:,min_index].*(1+0.15*randn(number_of_parameters))
  end

  #return (ec_array,pc_array)
end

# Run -
main()


text:
module MatrixReduction

using BenchmarkTools
using Ripserer
using Ripserer: AbstractRipsFiltration, simplex_type, BoundaryMatrix, CoboundaryMatrix, compute_death_simplices!, compute_intervals!
using IterTools

include(joinpath(@__DIR__, "utils.jl"))

"""
    collect_simplices(flt::AbstractRipsFiltration, dim)

Collect all simplices in `dim`-th dimension of Rips-like filtration.
"""
function collect_simplices(flt::AbstractRipsFiltration, dim)
    res = simplex_type(flt, dim)[]
    for vs in IterTools.subsets(vertices(flt), Val(dim + 1))
        σ = simplex(flt, Val(dim), vs)
        !isnothing(σ) && push!(res, σ)
    end
    return res
end

function coboundary_matrix(flt, dim_max; field=Mod{2}, implicit=true)
    _, to_reduce, to_skip = Ripserer.zeroth_intervals(
        flt, 0, false, field, false
    )
    matrix = CoboundaryMatrix{implicit}(field, flt, to_reduce, to_skip)
    for dim in 1:(dim_max - 1)
        compute_intervals!(matrix, 0, false, false)
        matrix = Ripserer.next_matrix(matrix, false)
    end
    return matrix
end

function boundary_matrix(flt, dim; field=Mod{2}, implicit=false)
    columns = collect_simplices(flt, dim+1)
    sort!(columns)
    return BoundaryMatrix{implicit}(field, flt, columns)
end

function involuted_boundary_matrix(flt, dim; field=Mod{2}, implicit=false)
    d = coboundary_matrix(flt, dim; field)
    columns, _ = compute_death_simplices!(d, false, 0)
    return BoundaryMatrix{implicit}(field, flt, columns)
end

"""
    reduce!(matrix; verbose=false, reps=true)

Reduce a matrix. Note: this functions is mutating! Don't run it twice on the same matrix.
"""
function reduce!(matrix; verbose=false, reps=true)
    compute_intervals!(matrix, 0, verbose, reps)
end

"""
    add_suite!(suite, data, name, dim; seconds=100, samples=10)

Add a suite to the `BenchmarkGroup`.
"""
function add_suite!(suite, data, name, dim, threshold, hom; seconds=100, samples=10)
    msg = " preparing $name ($dim). "
    println("="^length(msg))
    println(msg)
    suite[name] = BenchmarkGroup()

    if isnothing(threshold)
        ℜ = Rips(data)
    else
        ℜ = Rips(data; threshold)
    end
    d = coboundary_matrix(ℜ, dim)
    println("  d has ", length(d.columns_to_reduce), " columns.")
    if hom
        ∂ = boundary_matrix(ℜ, dim)
        println("  ∂ has ", length(∂.columns_to_reduce), " columns.")
    end
    ∂′= involuted_boundary_matrix(ℜ, dim)
    println("  ∂′ has ", length(∂′.columns_to_reduce), " columns.")

    suite[name]["cohomology"] = @benchmarkable reduce!(M) setup=(
        print("  reducing d for ", $name, "                    "); M=deepcopy($d)
    ) samples=samples seconds=seconds gcsample=true

    if hom
        suite[name]["homology"] = @benchmarkable reduce!(M) setup=(
            print("  reducing ∂ for ", $name, "                    "); M=deepcopy($∂)
        ) samples=samples seconds=seconds gcsample=true
    end

    suite[name]["involuted"] = @benchmarkable reduce!(M) setup=(
        print("  reducing ∂ for ", $name, "                    "); M=deepcopy($∂′)
    ) samples=samples seconds=seconds gcsample=true

    return suite
end

"""
    MatrixReduction.suite()

The suite will benchmark the reduction of a matrix in the `dim`-th dimension for a `Rips`
filtration built on each dataset in `datasets`. Notes:

* All (co)boundary matrices are built on the fly.
* The coboundary matrix uses an implicit representation of reduced columns, while the
  representation of boundary matrices is more explicit (they're reduced faster that way).
* The coboundary matrix uses the twist algorithm to select which columns need to be
  reduced. Without it, coboundary matrix reduction would be orders of magnitude slower (and
  the timings would not reflect real-world scenarios)
* The boundary matrix is the full matrix. The assumption here is that `dim` is the highest
  dimension we are computing in, so the twist algorithm would not help.
* The involuted matrix sometimes has fewer columns than the coboundary matrix because it is
  only collected up to the time of the last non-zero interval.
* The time it takes to collect the simplices is not accounted for. Sorting the simplices
  is. Including this time would work in favor of cohomology.
* Due to how Ripserer is implemented, 1-dimensional columns in the coboundary matrix are not
  sorted.
"""
function suite(
    datasets=[
        # filename            name           dim threshold homology?
        ("gcycle.dist",       "gcycle (1)",   1,  nothing,  true),
        ("gcycle.dist",       "gcycle (3)",   3,  40,       true),
        ("2_sphere_100.pts",  "2 sphere",     2,  nothing,  true),
        ("rand_4_50.pts",     "rand_4",       3,  nothing,  true),
        ("random16.dist",     "random16 (3)", 3,  nothing,  true),
        ("celegans.dist",     "celegans",     2,  nothing,  true),
        ("dragon1000.dist",   "dragon1000",   1,  nothing,  true),
        ("o3_1024.pts",       "o3_1024",      3,  1.8,      false),
        ("o3_4096.pts",       "o3_4096",      3,  1.4,      false),
        ("random16.dist",     "random16 (7)", 7,  nothing,  false),
        ("hiv.dist",          "hiv",          1,  nothing,  true),
    ];
    seconds=1000, samples=10,
)
    suite = BenchmarkGroup()
    for (file, name, dim, threshold, hom) in datasets
        data = load_data(joinpath(@__DIR__, "..", "datasets", file))
        add_suite!(suite, data, name, dim, threshold, hom; seconds, samples)
    end
    return suite
end

end


text:
# ---
# title: 1383. Maximum Performance of a Team
# id: problem1383
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Greedy, Sort
# link: <https://leetcode.com/problems/maximum-performance-of-a-team/description/>
# hidden: true
# ---
#
# There are `n` engineers numbered from 1 to `n` and two arrays: `speed` and
# `efficiency`, where `speed[i]` and `efficiency[i]` represent the speed and
# efficiency for the i-th engineer respectively. _Return the maximum
# **performance** of a team composed of at most `k` engineers, since the answer
# can be a huge number, return this modulo 10^9 + 7._
#
# The **performance** of a team is the sum of their engineers' speeds multiplied
# by the minimum efficiency among their engineers.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 6, speed = [2,10,3,1,5,8], efficiency = [5,4,3,9,7,2], k = 2
#     Output: 60
#     Explanation:
#     We have the maximum performance of the team by selecting engineer 2 (with speed=10 and efficiency=4) and engineer 5 (with speed=5 and efficiency=7). That is, performance = (10 + 5) * min(4, 7) = 60.
#
#
# **Example 2:**
#
#
#
#     Input: n = 6, speed = [2,10,3,1,5,8], efficiency = [5,4,3,9,7,2], k = 3
#     Output: 68
#     Explanation: This is the same example as the first but k = 3. We can select engineer 1, engineer 2 and engineer 5 to get the maximum performance of the team. That is, performance = (2 + 10 + 5) * min(5, 4, 7) = 68.
#
#
# **Example 3:**
#
#
#
#     Input: n = 6, speed = [2,10,3,1,5,8], efficiency = [5,4,3,9,7,2], k = 4
#     Output: 72
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 10^5`
#   * `speed.length == n`
#   * `efficiency.length == n`
#   * `1 <= speed[i] <= 10^5`
#   * `1 <= efficiency[i] <= 10^8`
#   * `1 <= k <= n`
#
#
## @lc code=start
using LeetCode

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


text:
VERSION < v"0.6" && warn("The package has not been tested on Julia versions < 0.6")
using SimplexSplitting
using BenchmarkTools
using InvariantDistribution

# Define a parent BenchmarkGroup to contain our suite
const suite = BenchmarkGroup()

#########
# Setup
#########

# Create a triangulation to work with
t = SimplexSplitting.triang_from_embedding(SimplexSplitting.Embedding(InvariantDistribution.invariant_gaussian_embedding(npts = 30, covariance = 0.3, tau = 1)))

# Add some child groups to our benchmark suite.
suite["markovmatrix"] = BenchmarkGroup(["discrete", "exact"])

#########
# Benchmark
#########
#suite["markovmatrix"]["exact"] = @benchmarkable mm_p(t)
suite["markovmatrix"]["discrete"] = @benchmarkable mm_discrete_dense(t, n_randpts = 100)

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `suite` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(suite, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(suite)
    BenchmarkTools.save(paramspath, params(suite));
end

max_stars_repo_name: kahaaga/InvariantDistribution

max_issues_repo_name: kahaaga/InvariantDistribution

max_forks_repo_name: kahaaga/InvariantDistribution


text:
# https://oxfordstrat.com/coasdfASD32/uploads/2016/03/How-Sharp-Is-the-Sharpe-Ratio.pdf

"""
    sharpe_ratio(returns; multiplier=1.0, risk_free=0.0)

Calculates the Sharpe Ratio (SR) according to the original definition by William F. Sharpe in 1966. For calculating the Sharpe Ratio according to Sharpe's revision in 1994, please see function `information_ratio` (IR).

# Formula

    SR = E[returns - risk_free] / std(returns) * sqrt(multiplier)

    IR = E[asset_returns - benchmark_returns] / std(asset_returns - benchmark_returns) * sqrt(multiplier)

# Arguments
- `returns`:        Vector of asset returns.
- `multiplier`:     Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.
- `risk_free`:      Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).

# Sources
- Sharpe, W. F. (1966). Mutual Fund Performance. Journal of Business.
- Sharpe, William F. (1994). The Sharpe Ratio. The Journal of Portfolio Management.
"""
function sharpe_ratio(returns; multiplier=1.0, risk_free=0.0)
    mean(returns .- risk_free) / std(returns) * sqrt(multiplier)
end

"""
    adjusted_sharpe_ratio(returns; multiplier=1.0, risk_free=0.0)

Calculates the adjusted Sharpe Ratio introduced by Pezier and White (2006) by penalizing negative skewness and excess kurtosis.

# Formula

    ASR = SR*[1 + (S/6)SR - (K-3)/24*SR^2] * sqrt(multiplier)

# Arguments
- `returns`:        Vector of asset returns.
- `multiplier`:     Optional scalar multiplier, i.e. use `12` to annualize monthly returns, and use `252` to annualize daily returns.
- `risk_free`:      Optional vector or scalar value denoting the risk-free return (must have same frequency as the provided returns, e.g. daily).

# Sources
- Pezier, Jaques and White, Anthony (2006). The Relative Merits of Investable Hedge Fund Indices and of Funds of Hedge Funds in Optimal Passive Portfolios. ICMA Centre Discussion Papers in Finance.
"""
function adjusted_sharpe_ratio(returns; multiplier=1.0, risk_free=0.0)
    excess = returns .- risk_free
    SR = mean(excess) / std(excess)
    S = skewness(excess)
    K = kurtosis(excess; method=:excess)
    SR*(1 + (S/6)SR - K/24*SR^2) * sqrt(multiplier)
end


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################
# src/aff_expr.jl
# Defines all types relating to affine expressions
# - GenericAffExpr              ∑ aᵢ xᵢ  +  c
#   - AffExpr                   Alias for (Float64, VariableRef)
#   - AffExprConstraint         AffExpr-in-set constraint
# Operator overloads in src/operators.jl
#############################################################################

# Utilities for OrderedDict
function _add_or_set!(dict::OrderedDict{K,V}, k::K, v::V) where {K,V}
    # Adding zero terms to this dictionary leads to unacceptable performance
    # degradations. See, e.g., https://github.com/jump-dev/JuMP.jl/issues/1946.
    if iszero(v)
        return dict  # No-op.
    end
    # TODO: This unnecessarily requires two lookups for k.
    dict[k] = get!(dict, k, zero(V)) + v
    return dict
end

function _new_ordered_dict(
    ::Type{K},
    ::Type{V},
    kv::AbstractArray{<:Pair},
) where {K,V}
    dict = OrderedDict{K,V}()
    sizehint!(dict, length(kv))
    for pair in kv
        _add_or_set!(dict, convert(K, pair.first), convert(V, pair.second))
    end
    return dict
end

function _new_ordered_dict(
    ::Type{K},
    ::Type{V},
    kv::Vararg{Pair,N},
) where {K,V,N}
    dict = OrderedDict{K,V}()
    sizehint!(dict, length(kv))
    for pair in kv
        _add_or_set!(dict, convert(K, pair.first), convert(V, pair.second))
    end
    return dict
end
# Shortcut for one and two arguments to avoid creating an empty dict and add
# elements one by one with `JuMP._add_or_set!`
function _new_ordered_dict(::Type{K}, ::Type{V}, kv::Pair) where {K,V}
    return OrderedDict{K,V}(kv)
end
function _new_ordered_dict(
    ::Type{K},
    ::Type{V},
    kv1::Pair,
    kv2::Pair,
) where {K,V}
    if isequal(kv1.first, kv2.first)
        return OrderedDict{K,V}(kv1.first => kv1.second + kv2.second)
    else
        return OrderedDict{K,V}(kv1, kv2)
    end
end

# As `!isbits(VariableRef)`, creating a pair allocates, with this API, we avoid
# this allocation.
function _build_aff_expr(constant::V, coef::V, var::K) where {V,K}
    terms = OrderedDict{K,V}()
    terms[var] = coef
    return GenericAffExpr{V,K}(constant, terms)
end
function _build_aff_expr(
    constant::V,
    coef1::V,
    var1::K,
    coef2::V,
    var2::K,
) where {V,K}
    if isequal(var1, var2)
        return _build_aff_expr(constant, coef1 + coef2, var1)
    end
    terms = OrderedDict{K,V}()
    terms[var1] = coef1
    terms[var2] = coef2
    return GenericAffExpr{V,K}(constant, terms)
end

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

"""
    mutable struct GenericAffExpr{CoefType,VarType} <: AbstractJuMPScalar
        constant::CoefType
        terms::OrderedDict{VarType,CoefType}
    end

An expression type representing an affine expression of the form:
``\\sum a_i x_i + c``.

## Fields

 * `.constant`: the constant `c` in the expression.
 * `.terms`: an `OrderedDict`, with keys of `VarType` and values of `CoefType`
   describing the sparse vector `a`.
"""
mutable struct GenericAffExpr{CoefType,VarType} <: AbstractJuMPScalar
    constant::CoefType
    terms::OrderedDict{VarType,CoefType}
end

"""
    variable_ref_type(::GenericAffExpr{C, V}) where {C, V}

A helper function used internally by JuMP and some JuMP extensions. Returns the
variable type `V` from a [`GenericAffExpr`](@ref)
"""
variable_ref_type(::GenericAffExpr{C,V}) where {C,V} = V

"""
    GenericAffExpr(constant::V, kv::AbstractArray{Pair{K,V}}) where {K,V}

Create a [`GenericAffExpr`](@ref) by passing a constant and a vector of pairs.

## Examples

```jldoctest; setup=:(using JuMP; model = Model(); @variable(model, x))
julia> GenericAffExpr(1.0, [x => 1.0])
x + 1
"""
function GenericAffExpr(constant::V, kv::AbstractArray{Pair{K,V}}) where {K,V}
    return GenericAffExpr{V,K}(constant, _new_ordered_dict(K, V, kv))
end

"""
    GenericAffExpr(constant::V, kv::Vararg{Pair{K,V},N}) where {K,V,N}

Create a [`GenericAffExpr`](@Ref) by passing a constant and pairs of additional
arguments.

## Examples

```jldoctest; setup=:(using JuMP; model = Model(); @variable(model, x))
julia> GenericAffExpr(1.0, x => 1.0)
x + 1
"""
function GenericAffExpr(constant::V, kv::Vararg{Pair{K,V},N}) where {K,V,N}
    return GenericAffExpr{V,K}(constant, _new_ordered_dict(K, V, kv...))
end

function GenericAffExpr{V,K}(constant, kv::AbstractArray{<:Pair}) where {K,V}
    return GenericAffExpr{V,K}(
        convert(V, constant),
        _new_ordered_dict(K, V, kv),
    )
end

function GenericAffExpr{V,K}(constant, kv::Vararg{Pair,N}) where {K,V,N}
    return GenericAffExpr{V,K}(
        convert(V, constant),
        _new_ordered_dict(K, V, kv...),
    )
end

function Base.iszero(expr::GenericAffExpr)
    return iszero(expr.constant) && all(iszero, values(expr.terms))
end
function Base.zero(::Type{GenericAffExpr{C,V}}) where {C,V}
    return GenericAffExpr{C,V}(zero(C), OrderedDict{V,C}())
end
function Base.one(::Type{GenericAffExpr{C,V}}) where {C,V}
    return GenericAffExpr{C,V}(one(C), OrderedDict{V,C}())
end
Base.zero(a::GenericAffExpr) = zero(typeof(a))
Base.one(a::GenericAffExpr) = one(typeof(a))
Base.copy(a::GenericAffExpr) = GenericAffExpr(copy(a.constant), copy(a.terms))
Base.broadcastable(a::GenericAffExpr) = Ref(a)

"""
    coefficient(a::GenericAffExpr{C,V}, v::V) where {C,V}

Return the coefficient associated with variable `v` in the affine expression `a`.
"""
coefficient(a::GenericAffExpr{C,V}, v::V) where {C,V} = get(a.terms, v, zero(C))
coefficient(a::GenericAffExpr{C,V}, v1::V, v2::V) where {C,V} = zero(C)

"""
    drop_zeros!(expr::GenericAffExpr)

Remove terms in the affine expression with `0` coefficients.
"""
function drop_zeros!(expr::GenericAffExpr)
    for (key, coef) in expr.terms
        if iszero(coef)
            delete!(expr.terms, key)
        end
    end
    return
end

GenericAffExpr{C,V}() where {C,V} = zero(GenericAffExpr{C,V})

function _affine_coefficient(f::GenericAffExpr{C,V}, variable::V) where {C,V}
    return get(f.terms, variable, zero(C))
end

"""
    map_coefficients_inplace!(f::Function, a::GenericAffExpr)

Apply `f` to the coefficients and constant term of an [`GenericAffExpr`](@ref)
`a` and update them in-place.

See also: [`map_coefficients`](@ref)

## Example

```jldoctest; setup=:(using JuMP; model = Model(); @variable(model, x))
julia> a = GenericAffExpr(1.0, x => 1.0)
x + 1

julia> map_coefficients_inplace!(c -> 2 * c, a)
2 x + 2

julia> a
2 x + 2
```
"""
function map_coefficients_inplace!(f::Function, a::GenericAffExpr)
    # The iterator remains valid if existing elements are updated.
    for (coef, var) in linear_terms(a)
        a.terms[var] = f(coef)
    end
    a.constant = f(a.constant)
    return a
end

"""
    map_coefficients(f::Function, a::GenericAffExpr)

Apply `f` to the coefficients and constant term of an [`GenericAffExpr`](@ref)
`a` and return a new expression.

See also: [`map_coefficients_inplace!`](@ref)

## Example

```jldoctest; setup=:(using JuMP; model = Model(); @variable(model, x))
julia> a = GenericAffExpr(1.0, x => 1.0)
x + 1

julia> map_coefficients(c -> 2 * c, a)
2 x + 2

julia> a
x + 1
```
"""
function map_coefficients(f::Function, a::GenericAffExpr)
    return map_coefficients_inplace!(f, copy(a))
end

Base.sizehint!(a::GenericAffExpr, n::Int) = sizehint!(a.terms, n)

"""
    value(ex::GenericAffExpr, var_value::Function)

Evaluate `ex` using `var_value(v)` as the value for each variable `v`.
"""
function value(ex::GenericAffExpr{T,V}, var_value::Function) where {T,V}
    S = Base.promote_op(var_value, V)
    U = Base.promote_op(*, T, S)
    ret = convert(U, ex.constant)
    for (var, coef) in ex.terms
        ret += coef * var_value(var)
    end
    return ret
end

"""
    constant(aff::GenericAffExpr{C, V})::C

Return the constant of the affine expression.
"""
constant(aff::GenericAffExpr) = aff.constant

# Iterator protocol - iterates over tuples (aᵢ,xᵢ)
struct LinearTermIterator{GAE<:GenericAffExpr}
    aff::GAE
end

"""
    linear_terms(aff::GenericAffExpr{C, V})

Provides an iterator over coefficient-variable tuples `(a_i::C, x_i::V)` in the
linear part of the affine expression.
"""
linear_terms(aff::GenericAffExpr) = LinearTermIterator(aff)

_reverse_pair_to_tuple(p::Pair) = (p.second, p.first)
function Base.iterate(lti::LinearTermIterator)
    ret = iterate(lti.aff.terms)
    if ret === nothing
        return nothing
    else
        return _reverse_pair_to_tuple(ret[1]), ret[2]
    end
end
function Base.iterate(lti::LinearTermIterator, state)
    ret = iterate(lti.aff.terms, state)
    if ret === nothing
        return nothing
    else
        return _reverse_pair_to_tuple(ret[1]), ret[2]
    end
end
Base.length(lti::LinearTermIterator) = length(lti.aff.terms)
function Base.eltype(lti::LinearTermIterator{GenericAffExpr{C,V}}) where {C,V}
    return Tuple{C,V}
end

"""
    add_to_expression!(expression, terms...)

Updates `expression` *in place* to `expression + (*)(terms...)`. This is
typically much more efficient than `expression += (*)(terms...)`. For example,
`add_to_expression!(expression, a, b)` produces the same result as `expression
+= a*b`, and `add_to_expression!(expression, a)` produces the same result as
`expression += a`.

Only a few methods are defined, mostly for internal use, and only for the cases
when (1) they can be implemented efficiently and (2) `expression` is capable of
storing the result. For example, `add_to_expression!(::AffExpr, ::VariableRef,
::VariableRef)` is not defined because a `GenericAffExpr` cannot store the
product of two variables.
"""
function add_to_expression! end

# With one factor.

function add_to_expression!(aff::GenericAffExpr, other::_Constant)
    aff.constant += _constant_to_number(other)
    return aff
end

function add_to_expression!(aff::GenericAffExpr{C,V}, new_var::V) where {C,V}
    _add_or_set!(aff.terms, new_var, one(C))
    return aff
end

function add_to_expression!(
    aff::GenericAffExpr{C,V},
    other::GenericAffExpr{C,V},
) where {C,V}
    # Note: merge!() doesn't appear to call sizehint!(). Is this important?
    merge!(+, aff.terms, other.terms)
    aff.constant += other.constant
    return aff
end

# With two factors.

function add_to_expression!(
    aff::GenericAffExpr{C,V},
    new_coef::_Constant,
    new_var::V,
) where {C,V}
    _add_or_set!(aff.terms, new_var, convert(C, _constant_to_number(new_coef)))
    return aff
end

function add_to_expression!(
    aff::GenericAffExpr{C,V},
    new_var::V,
    new_coef::_Constant,
) where {C,V}
    return add_to_expression!(aff, new_coef, new_var)
end

function add_to_expression!(
    aff::GenericAffExpr{C,V},
    coef::_Constant,
    other::GenericAffExpr{C,V},
) where {C,V}
    sizehint!(aff, length(linear_terms(aff)) + length(linear_terms(other)))
    for (term_coef, var) in linear_terms(other)
        _add_or_set!(aff.terms, var, coef * term_coef)
    end
    aff.constant += coef * other.constant
    return aff
end

function add_to_expression!(
    aff::GenericAffExpr{C,V},
    other::GenericAffExpr{C,V},
    coef::_Constant,
) where {C,V}
    return add_to_expression!(aff, coef, other)
end

function Base.isequal(
    aff::GenericAffExpr{C,V},
    other::GenericAffExpr{C,V},
) where {C,V}
    return isequal(aff.constant, other.constant) &&
           isequal(aff.terms, other.terms)
end

Base.hash(aff::GenericAffExpr, h::UInt) = hash(aff.constant, hash(aff.terms, h))

function SparseArrays.dropzeros(aff::GenericAffExpr)
    result = copy(aff)
    for (coef, var) in linear_terms(aff)
        if iszero(coef)
            delete!(result.terms, var)
        end
    end
    if iszero(result.constant)
        # This is to work around isequal(0.0, -0.0) == false.
        result.constant = zero(typeof(result.constant))
    end
    return result
end

"""
    isequal_canonical(
        aff::GenericAffExpr{C,V},
        other::GenericAffExpr{C,V}
    ) where {C,V}

Return `true` if `aff` is equal to `other` after dropping zeros and disregarding
the order. Mainly useful for testing.
"""
function isequal_canonical(
    aff::GenericAffExpr{C,V},
    other::GenericAffExpr{C,V},
) where {C,V}
    aff_nozeros = dropzeros(aff)
    other_nozeros = dropzeros(other)
    # Note: This depends on equality of OrderedDicts ignoring order.
    # This is the current behavior, but it seems questionable.
    return isequal(aff_nozeros, other_nozeros)
end

function Base.convert(::Type{GenericAffExpr{T,V}}, v::V) where {T,V}
    return GenericAffExpr(zero(T), v => one(T))
end
function Base.convert(::Type{GenericAffExpr{T,V}}, v::_Constant) where {T,V}
    return GenericAffExpr{T,V}(convert(T, _constant_to_number(v)))
end

"""
    AffExpr

Alias for `GenericAffExpr{Float64,VariableRef}`, the specific
[`GenericAffExpr`](@ref) used by JuMP.
"""
const AffExpr = GenericAffExpr{Float64,VariableRef}

# Check all coefficients are finite, i.e. not NaN, not Inf, not -Inf
function _assert_isfinite(a::AffExpr)
    for (coef, var) in linear_terms(a)
        isfinite(coef) || error("Invalid coefficient $coef on variable $var.")
    end
    if isnan(a.constant)
        error(
            "Expression contains an invalid NaN constant. This could be " *
            "produced by `Inf - Inf`.",
        )
    end
end

"""
    value(v::GenericAffExpr; result::Int = 1)

Return the value of the `GenericAffExpr` `v` associated with result index
`result` of the most-recent solution returned by the solver.

Replaces `getvalue` for most use cases.

See also: [`result_count`](@ref).
"""
function value(a::GenericAffExpr; result::Int = 1)
    return value(a, (x) -> value(x; result = result))
end

function check_belongs_to_model(a::GenericAffExpr, model::AbstractModel)
    for variable in keys(a.terms)
        check_belongs_to_model(variable, model)
    end
end

# Note: No validation is performed that the variables in the AffExpr belong to
# the same model. The verification is done in `check_belongs_to_model` which
# should be called before calling `MOI.ScalarAffineFunction`.
function MOI.ScalarAffineFunction(a::AffExpr)
    _assert_isfinite(a)
    terms = MOI.ScalarAffineTerm{Float64}[
        MOI.ScalarAffineTerm(t[1], index(t[2])) for t in linear_terms(a)
    ]
    return MOI.ScalarAffineFunction(terms, a.constant)
end

"""
    moi_function(x)

Given a JuMP object `x`, return the MathOptInterface equivalent.

See also: [`jump_function`](@ref).
"""
function moi_function end

"""
    moi_function_type(::Type{T}) where {T}

Given a JuMP object type `T`, return the MathOptInterface equivalent.

See also: [`jump_function_type`](@ref).
"""
function moi_function_type end

"""
    jump_function(x)

Given an MathOptInterface object `x`, return the JuMP equivalent.

See also: [`moi_function`](@ref).
"""
function jump_function end

"""
    jump_function_type(::Type{T}) where {T}

Given an MathOptInterface object type `T`, return the JuMP equivalent.

See also: [`moi_function_type`](@ref).
"""
function jump_function_type end

moi_function(a::GenericAffExpr) = MOI.ScalarAffineFunction(a)
function moi_function_type(::Type{<:GenericAffExpr{T}}) where {T}
    return MOI.ScalarAffineFunction{T}
end

function AffExpr(m::Model, f::MOI.ScalarAffineFunction)
    aff = AffExpr()
    for t in f.terms
        add_to_expression!(aff, t.coefficient, VariableRef(m, t.variable_index))
    end
    aff.constant = f.constant
    return aff
end
function jump_function_type(
    ::Model,
    ::Type{MOI.ScalarAffineFunction{T}},
) where {T}
    return GenericAffExpr{T,VariableRef}
end
function jump_function(model::Model, f::MOI.ScalarAffineFunction{T}) where {T}
    return GenericAffExpr{T,VariableRef}(model, f)
end
function jump_function_type(
    ::Model,
    ::Type{MOI.VectorAffineFunction{T}},
) where {T}
    return Vector{GenericAffExpr{T,VariableRef}}
end
function jump_function(model::Model, f::MOI.VectorAffineFunction{T}) where {T}
    return GenericAffExpr{T,VariableRef}[
        GenericAffExpr{T,VariableRef}(model, f) for f in MOIU.eachscalar(f)
    ]
end

"""
    _fill_vaf!(terms::Vector{<:MOI.VectorAffineTerm}, offset::Int, oi::Int,
               aff::AbstractJuMPScalar)

Fills the vectors terms at indices starting at `offset+1` with the affine terms
of `aff`. The output index for all terms is `oi`. Return the index of the last
term added.
"""
function _fill_vaf!(
    terms::Vector{<:MOI.VectorAffineTerm},
    offset::Int,
    oi::Int,
    aff::AbstractJuMPScalar,
)
    i = 1
    for (coef, var) in linear_terms(aff)
        terms[offset+i] = MOI.VectorAffineTerm(
            Int64(oi),
            MOI.ScalarAffineTerm(coef, index(var)),
        )
        i += 1
    end
    return offset + length(linear_terms(aff))
end

function MOI.VectorAffineFunction(affs::Vector{AffExpr})
    len = sum(aff -> length(linear_terms(aff)), affs)
    terms = Vector{MOI.VectorAffineTerm{Float64}}(undef, len)
    constant = Vector{Float64}(undef, length(affs))
    offset = 0
    for (i, aff) in enumerate(affs)
        constant[i] = aff.constant
        offset = _fill_vaf!(terms, offset, i, aff)
    end
    return MOI.VectorAffineFunction(terms, constant)
end
moi_function(a::Vector{<:GenericAffExpr}) = MOI.VectorAffineFunction(a)
function moi_function_type(::Type{<:Vector{<:GenericAffExpr{T}}}) where {T}
    return MOI.VectorAffineFunction{T}
end

# TODO: Find somewhere to put this error message.
#add_constraint(m::Model, c::Array{AffExprConstraint}) =
#    error("The operators <=, >=, and == can only be used to specify scalar constraints. If you are trying to add a vectorized constraint, use the element-wise dot comparison operators (.<=, .>=, or .==) instead")

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using EvolutionaryStrategies
using Cambrian
using Test

cfg = get_config("test.yaml")

function test_es(constructor, init)
    center = ones(cfg.n_genes)
    e = constructor(cfg, i->sphere(i; center=ones(cfg.n_genes)))
    @test ~any(isnan.(e.state.μ))

    population = init(e.config, e.state)
    @test length(population) == e.config.n_population
    @test all([p.fitness[1] .== -Inf for p in population])

    step!(e)
    @test all([p.fitness[1] .!= -Inf for p in e.population])

    @test ~any(isnan.(e.state.μ))
    best = e.elites[end]
    @test best.fitness[1] <= 0.0
    @test e.gen == 1

    for i in 1:10
        step!(e)
        new_best = e.elites[end]
        @test new_best.fitness[1] >= best.fitness[1]
        @test ~any(isnan.(e.state.μ))
        @test all([p.fitness[1] .!= -Inf for p in e.population])
    end
    @timev step!(e)
    @test e.elites[end].fitness[1] < 0.1
end

@testset "Exponential Natural Evolution Strategy" begin
    test_es(xNES, EvolutionaryStrategies.xnes_init)
end

@testset "Separable Natural Evolution Strategy" begin
    test_es(sNES, EvolutionaryStrategies.snes_init)
end


text:
"""
## module Potentials

### Summary

This module implements some basic interatomic potentials in pure Julia, as well
as provides building blocks and prototypes for further implementations
The implementation is done in such a way that they can be used either in "raw"
form or withinabstractframeworks.

### Types

### `evaluate`, `evaluate_d`, `evaluate_dd`, `grad`

### The `@D`, `@DD`, `@GRAD` macros

TODO: write documentation

"""
module Potentials

using JuLIP: AbstractAtoms, AbstractCalculator,
      JVec, mat, vec, JMat, SVec, vecs, SMat,
      positions, set_positions!

using StaticArrays: @SMatrix

using NeighbourLists

using LinearAlgebra: norm

using SparseArrays: sparse

import JuLIP: energy, forces, cutoff, virial, hessian_pos, hessian,
              site_energies, r_sum,
              site_energy, site_energy_d,
              energy!, forces!, virial!,
              alloc_temp, alloc_temp_d, alloc_temp_dd

export PairPotential, SitePotential, ZeroSitePotential

# the following are prototypes for internal functions around which IPs are
# defined

function evaluate end
function evaluate_d end
function evaluate_ed end
function evaluate_dd end
function evaluate! end
function evaluate_d! end
function evaluate_dd! end
function grad end
function precond end

include("potentials_base.jl")
# * @pot, @D, @DD, @GRAD and related things

"""
`SitePotential`:abstractsupertype for generic site potentials
"""
abstract type SitePotential <: AbstractCalculator end

"""
`PairPotential`:abstractsupertype for pair potentials

Can also be used as a constructor for analytic pair potentials, e.g.,
```julia
lj = @analytic r -> r^(-12) - 2 * r^(-6)
```
"""
abstract type PairPotential <: SitePotential end

evaluate(V::SitePotential, R) =
      evaluate!(alloc_temp(V, length(R)), V, R)

evaluate_d(V::SitePotential, R::AbstractVector{JVec{T}}) where {T} =
      evaluate_d!(zeros(JVec{T}, length(R)),
                  alloc_temp_d(V, length(R)),
                  V, R)

evaluate_dd(V::SitePotential, R::AbstractVector{JVec{T}}) where {T} =
      evaluate_dd!(zeros(JMat{T}, length(R), length(R)),
                   alloc_temp_dd(V, length(R)),
                   V, R)

evaluate(V::PairPotential, r::Number) = evaluate!(alloc_temp(V, 1), V, r)
evaluate_d(V::PairPotential, r::Number) = evaluate_d!(alloc_temp_d(V, 1), V, r)
evaluate_dd(V::PairPotential, r::Number) = evaluate_dd!(alloc_temp_d(V, 1), V, r)

NeighbourLists.sites(at::AbstractAtoms, rcut::AbstractFloat) =
      sites(neighbourlist(at, rcut))

NeighbourLists.pairs(at::AbstractAtoms, rcut::AbstractFloat) =
      pairs(neighbourlist(at, rcut))

"a site potential that just returns zero"
mutable struct ZeroSitePotential <: SitePotential
end

@pot ZeroSitePotential

cutoff(::ZeroSitePotential) = Bool(0)
energy(V::ZeroSitePotential, at::AbstractAtoms{T}; kwargs...) where T = zero(T)
forces(V::ZeroSitePotential, at::AbstractAtoms{T}; kwargs...) where T = zeros(JVec{T}, length(at))
evaluate!(tmp, p::ZeroSitePotential, args...) = Bool(0)
evaluate_d!(dEs, tmp, V::ZeroSitePotential, args...) = fill!(dEs, zero(eltype(dEs)))
evaluate_dd!(hEs, tmp, V::ZeroSitePotential, args...) = fill!(hEs, zero(eltype(hEs)))

# Implementation of a generic site potential
# ================================================

alloc_temp(V::SitePotential, at::AbstractAtoms) =
      alloc_temp(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp(V::SitePotential, N::Integer) = ( R = zeros(JVecF, N), )

alloc_temp_d(V::SitePotential, at::AbstractAtoms) =
      alloc_temp_d(V, maxneigs(neighbourlist(at, cutoff(V))))

alloc_temp_d(V::SitePotential, N::Integer) = ( R = zeros(JVecF, N),
                                              dV = zeros(JVecF, N) )

alloc_temp_dd(V::SitePotential, N::Integer) = nothing

energy(V::SitePotential, at::AbstractAtoms; kwargs...) =
      energy!(alloc_temp(V, at), V, at; kwargs...)

virial(V::SitePotential, at::AbstractAtoms; kwargs...) =
      virial!(alloc_temp_d(V, at), V, at; kwargs...)

forces(V::SitePotential, at::AbstractAtoms{T}; kwargs...) where {T} =
      forces!(zeros(JVec{T}, length(at)), alloc_temp_d(V, at), V, at; kwargs...)

function energy!(tmp, V::SitePotential, at::AbstractAtoms{T};
                 domain=1:length(at)) where {T}
   E = zero(T)
   nlist = neighbourlist(at, cutoff(V))
   for i in domain
      _j, R = neigs!(tmp.R, nlist, i)
      E += evaluate!(tmp, V, R)
   end
   return E
end

function forces!(frc, tmp, V::SitePotential, at::AbstractAtoms{T};
                 domain=1:length(at), reset=true) where {T}
   if reset; fill!(frc, zero(JVec{T})); end
   nlist = neighbourlist(at, cutoff(V))
   for i in domain
      j, R = neigs!(tmp.R, nlist, i)
      evaluate_d!(tmp.dV, tmp, V, R)
      for a = 1:length(j)
         frc[j[a]] -= tmp.dV[a]
         frc[i]    += tmp.dV[a]
      end
   end
   return frc
end

site_virial(dV, R::AbstractVector{JVec{T}}) where {T} =  (
      length(R) > 0 ? (- sum( dVi * Ri' for (dVi, Ri) in zip(dV, R) ))
                    : zero(JMat{T}) )

function virial!(tmp, V::SitePotential, at::AbstractAtoms{T};
                 domain=1:length(at)) where {T}
   vir = zero(JMat{T})
   nlist = neighbourlist(at, cutoff(V))
   for i in domain
      _j, R = neigs!(tmp.R, nlist, i)
      evaluate_d!(tmp.dV, tmp, V, R)
      vir += site_virial(tmp.dV, R)
   end
   return vir
end

site_energies(V::SitePotential, at::AbstractAtoms{T}; kwargs...) where {T} =
      site_energies!(zeros(T, length(at)), alloc_temp(V, at), V, at)

function site_energies!(Es, tmp, V::SitePotential, at::AbstractAtoms{T};
         domain = 1:length(at)) where {T}
   nlist = neighbourlist(at, cutoff(V))
   for i in domain
      _j, R = neigs!(tmp.R, nlist, i)
      Es[i] = evaluate!(tmp, V, R)
   end
   return Es
end

site_energy(V::SitePotential, at::AbstractAtoms, i0::Integer) =
      energy(V, at; domain = (i0,))

site_energy_d(V::SitePotential, at::AbstractAtoms, i0::Integer) =
      rmul!(forces(V, at; domain = (i0,)), -one(eltype(at)))

# ------------------------------------------------
#  specialisation for Pair potentials

include("analyticpotential.jl")
# * AnalyticFunction
# * F64fun, WrappedAnalyticFunction

include("cutoffs.jl")
#   * SWCutoff
#   * ShiftCutoff
#   * SplineCutoff

include("pairpotentials.jl")
# * PairCalculator
# * LennardJonesPotential
# * MorsePotential
# * SimpleExponential
# * WrappedPairPotential

include("adsite.jl")
# * ADPotential : Site potential using ForwardDiff

include("stillingerweber.jl")
# * type StillingerWeber

include("splines.jl")
include("eam.jl")
# EAM, FinnisSinclair

include("onebody.jl")
# code for 1-body functions ; TODO: move into `nbody`

include("multi.jl")
# experimental multi-species code
# -> eventually this is to be integrated into all the main codebase

include("hessians.jl")
# code for hessians of site potentials

include("emt.jl")
# a simple analytic EAM potential. (EMT -> embedded medium theory)

end


text:
using Random, JuMP, JuDGE, Test

include("solvers/setup_gurobi.jl")
#include("solvers/setup_cplex.jl")
#include("solvers/setup_coin.jl")
#include("solvers/setup_glpk.jl")

function VRP(;seed::Int=648,num_packages=20,num_drivers=4,max_deliveries=6,use_heuristic=true,use_branch_callback=false,reltol=0.001,rlx_reltol=0.0001)
	Random.seed!(seed)
	package_values=rand(num_packages+1)
	positions=rand(num_packages+1,2)
	package_values[1]=0.0
	positions[1,1]=0.5
	positions[1,2]=0.5

	dist=Dict{Any,Float64}()

	for i in 1:num_packages+1
	    for j in 1:num_packages+1
	        dist[i,j]=((positions[i,1]-positions[j,1])^2+(positions[i,2]-positions[j,2])^2)^0.5
	    end
	end
	dist

	mytree = narytree(0,0)

	mytree.ext[:sum_max]=num_drivers
	mytree.ext[:sum_min]=1

	num_packages+=1
	max_deliveries+=1

	function initialise_columns(judge::JuDGEModel)
		sequences=[]
		package=2
		for i in 1:num_drivers
			sequence=[1]
			for j in 1:(num_packages-1) ÷ num_drivers
				push!(sequence,package)
			end
			push!(sequences,sequence)
			package+=1
		end

		for s in sequences
			sp=sub_problems(judge.tree,deliveries=s,mipgap=0.0)
			sp.ext[:form]=:binarycolumns
			sp.ext[:objective]=@variable(sp, obj)
			sp.ext[:objective_con]=@constraint(sp, sp.ext[:objective]-objective_function(sp) == 0)
			optimize!(sp)
			col=JuDGE.add_column(judge.master_problem, sp, judge.tree,branches=judge.ext[:branches])
		end
	end

	function sub_problems(node;deliveries=nothing,mipgap=0.00)
		sp=JuMP.Model(JuDGE_SP_Solver)
	    @shutdown(sp, delivered[1:num_packages], Bin)

		@enforced(sp, 1<=numdrivers<=1, Int, lb=1, ub=num_drivers)

		@capitalcosts(sp, -sum(package_values[i]*delivered[i]*10 for i in 2:num_packages)+2*numdrivers)

		@variable(sp,trip[1:num_packages,1:num_packages],Bin)
		@variable(sp,0<=counter[1:num_packages,1:num_packages]<=max_deliveries,Int)
		@constraint(sp,delivered[1]<=0)
	    @constraint(sp,sum(trip[1,j] for j in 1:num_packages)==1)
		@constraint(sp,sum(counter[1,j] for j in 1:num_packages)==1)

	    @constraint(sp,sum(trip[i,1] for i in 1:num_packages)==1)

	    for i in 2:num_packages
	        @constraint(sp,sum(trip[j,i] - trip[i,j] for j in 1:num_packages)==0)
			@constraint(sp,sum(counter[i,j] - trip[j,i] - counter[j,i] for j in 1:num_packages)==0)
	    end

	    for i in 2:num_packages
	        @constraint(sp,sum(trip[i,j] for j in 1:num_packages)>=delivered[i])
	    end
		@constraint(sp,sum(trip[i,i] for i in 1:num_packages)==0)

		if typeof(deliveries)==Nothing
			@constraint(sp,sum(trip[i,j] for i in 1:num_packages for j in 1:num_packages)<=max_deliveries)
			@constraint(sp,sum(counter[i,j] for i in 1:num_packages for j in 1:num_packages)<=max_deliveries*(max_deliveries+1)/2)
			for i in 1:num_packages
				for j in 1:num_packages
					@constraint(sp,counter[i,j]<=max_deliveries*trip[i,j])
				end
			end

		else
			@constraint(sp,sum(trip[i,j] for i in 1:num_packages for j in 1:num_packages)<=length(deliveries))
			@constraint(sp,sum(counter[i,j] for i in 1:num_packages for j in 1:num_packages)<=length(deliveries)*(length(deliveries)+1)/2)
			for i in 1:num_packages
				for j in 1:num_packages
					@constraint(sp,counter[i,j]<=length(deliveries)*trip[i,j])
				end
			end
			for i in 2:num_packages
				if i in deliveries
					@constraint(sp,delivered[i]==1)
				else
					@constraint(sp,delivered[i]==0)
				end
			end
		end

		@objective(sp, Min, sum(dist[i,j]*trip[i,j] for i in 1:num_packages for j in 1:num_packages))

	    sp
	end

	function VRP_branch(jmodel::JuDGEModel, inttol::Float64)
		master=jmodel.master_problem
		subproblems=jmodel.sub_problems
		tree=jmodel.tree

		branches=JuDGE.variable_branch(jmodel, inttol)
		if branches!=nothing
			return branches
		end

		onebranch=zeros(length(subproblems[tree][:delivered]),length(subproblems[tree][:delivered]))
		map(Main.eval,[:(jmodel=$jmodel)])
		for i in 1:length(subproblems[tree][:delivered])
			for j in i+1:length(subproblems[tree][:delivered])
				for col in master.ext[:columns][tree]
					if i in keys(col.coeffs[:delivered]) && j in keys(col.coeffs[:delivered])
						onebranch[i,j]+=value(col.var)
					end
				end
				onebranch[i,j] = onebranch[i,j] > 1.0-10^-6 ? 0 : onebranch[i,j]
			end
		end

		index=0
		#onebranch=abs.(onebranch-0.5*ones(length(subproblems[tree][:delivered]),length(subproblems[tree][:delivered])))

		# if minimum(onebranch)<0.4999
		# 	index=argmin(onebranch)
		# end

		if maximum(onebranch)>0.1
			index=argmax(onebranch)
		end

		if index!=0
			println("Branching on: "*string(index))

			filter1(col::JuDGE.Column)=((index[1] in keys(col.coeffs[:delivered])) ⊻ (index[2] in keys(col.coeffs[:delivered]))) ? :ban : nothing
			filter2(col::JuDGE.Column)=((index[1] in keys(col.coeffs[:delivered])) & (index[2] in keys(col.coeffs[:delivered]))) ? :ban : nothing

			constraint1=JuDGE.BranchConstraint(@expression(subproblems[tree],subproblems[tree][:delivered][index[1]]-subproblems[tree][:delivered][index[2]]),:eq,0.0,tree)
			constraint2=JuDGE.BranchConstraint(@expression(subproblems[tree],subproblems[tree][:delivered][index[1]]+subproblems[tree][:delivered][index[2]]),:le,1.0,tree)

			return [JuDGE.Branch(constraint1,filter1),JuDGE.Branch(constraint2,filter2)]
		end

	end

	function bp_callback(termination::Termination,search::Symbol,log::Array{Array{JuDGE.ConvergenceState,1},1})
		termination.rlx_reltol=0.0002
		return termination,search
		# if (termination.allow_frac==:first_fractional && log[end][end].num_frac==0)
		# 	termination.allow_frac=:binary_solve_return_relaxation
		# 	return (termination,:lowestLB)
		# else
		# 	termination.allow_frac=:first_fractional
		#  	return (termination,:depth_first_dive)
		# end
	end

	function heuristic(incumbent::JuDGEModel)
		JuDGE.solve_binary(incumbent,10^-9,0.01,true,nothing)
		sequences=get_tours(incumbent,:nodes)
		JuDGE.remove_binary(incumbent)
		changed=true
		improve=0.0
		while changed
			changed=false
			for n=0:3
				for i=1:length(sequences)-1
					for j=i+1:length(sequences)
						for k in 1:length(sequences[i])
							skip=false
							for p=0:n
								if sequences[i][(k+p-1) % length(sequences[i]) + 1]==1
									skip=true
									break
								end
							end
							if skip
								continue
							end
							for l in 1:length(sequences[j])
								skip=false
								for p=0:n
									if sequences[j][(l+p-1) % length(sequences[j]) + 1]==1
										skip=true
										break
									end
								end
								if skip
									continue
								end
								prev_k = (k-2+length(sequences[i])) % length(sequences[i]) + 1
								k2=(k+n-1) % length(sequences[i]) + 1
								next_k = (k+n) % length(sequences[i]) + 1
								prev_l = (l-2+length(sequences[j])) % length(sequences[j]) + 1
								l2=(l+n-1) % length(sequences[j]) + 1
								next_l = (l+n) % length(sequences[j]) + 1
								temp=0.0
								temp-=dist[sequences[i][k2],sequences[i][next_k]]
								temp-=dist[sequences[i][prev_k],sequences[i][k]]
								temp-=dist[sequences[j][l2],sequences[j][next_l]]
								temp-=dist[sequences[j][prev_l],sequences[j][l]]

								temp2=0.0
								temp2+=dist[sequences[i][k2],sequences[j][next_l]]
								temp2+=dist[sequences[j][l2],sequences[i][next_k]]
								temp2+=dist[sequences[i][prev_k],sequences[j][l]]
								temp2+=dist[sequences[j][prev_l],sequences[i][k]]

								if temp+temp2>=-10^-6
									temp2=0.0
									temp2+=dist[sequences[i][k2],sequences[j][prev_l]]
									temp2+=dist[sequences[j][l2],sequences[i][next_k]]
									temp2+=dist[sequences[i][prev_k],sequences[j][l]]
									temp2+=dist[sequences[j][next_l],sequences[i][k]]
								end

								if temp+temp2>=-10^-6
									temp2=0.0
									temp2+=dist[sequences[i][k2],sequences[j][prev_l]]
									temp2+=dist[sequences[j][l2],sequences[i][prev_k]]
									temp2+=dist[sequences[i][next_k],sequences[j][l]]
									temp2+=dist[sequences[j][next_l],sequences[i][k]]
								end

								if temp+temp2>=-10^-6
									temp2=0.0
									temp2+=dist[sequences[i][k2],sequences[j][next_l]]
									temp2+=dist[sequences[j][l2],sequences[i][prev_k]]
									temp2+=dist[sequences[i][next_k],sequences[j][l]]
									temp2+=dist[sequences[j][prev_l],sequences[i][k]]
								end

								if temp+temp2<-10^-6
									changed=true
									improve+=temp+temp2
									for p=0:n
										store=sequences[i][(k+p-1)%length(sequences[i])+1]
										sequences[i][(k+p-1)%length(sequences[i])+1]=sequences[j][(l+p-1)%length(sequences[j])+1]
										sequences[j][(l+p-1)%length(sequences[j])+1]=store
									end
									break
								end
							end
							if changed
								break
							end
						end
						if changed
							break
						end
					end
					if changed
						break
					end
				end
				if changed
					break
				end
			end
			if changed
				break
			end
		end
		if improve<0.0
			println("Heuristic improved solution by "*string(improve))
			for s in sequences
				sp=sub_problems(incumbent.tree,deliveries=s,mipgap=0.0)
				sp.ext[:form]=:binary
				sp.ext[:objective]=@variable(sp, obj)
				sp.ext[:objective_con]=@constraint(sp, sp.ext[:objective]-objective_function(sp) == 0)
				optimize!(sp)
				col=JuDGE.add_column(incumbent.master_problem, sp, incumbent.tree,branches=incumbent.ext[:branches])
			end
		end
		return improve
	end

	function get_tours(jmodel::JuDGEModel,type::Symbol)
		active = JuDGE.get_active_columns(jmodel,inttol=10^-5)

		sequences=[]
		for a in active[mytree]
			deliveries=[]
			for d in keys(a[1].coeffs[:delivered])
				push!(deliveries,d)
			end
			if 1 ∉ deliveries
				push!(deliveries,1)
			end
			sp=sub_problems(mytree,deliveries=deliveries,mipgap=0.0)
			optimize!(sp)
			if type==:arcs
				arcs=Tuple{Int,Int}[]
				count=length(deliveries)
				i=1
				while count>0
					for j in 1:num_packages
						if value(sp[:trip][i,j])>10^-7
							push!(arcs,(i,j))
							i=j
							count-=1
							break
						end
					end
				end
				push!(sequences,arcs)
			elseif type==:nodes
				nodes=Int[]
				count=length(deliveries)
				i=1
				while count>0
					for j in 1:num_packages
						if value(sp[:trip][i,j])>10^-7
							push!(nodes,i)
							i=j
							count-=1
							break
						end
					end
				end
				push!(sequences,nodes)
			end
		end
		sequences
	end

	function show_vrp(jmodel::JuDGEModel)

		sequences = get_tours(jmodel,:arcs)

		nodetemp = "{id:LABEL,label:\"LABEL\",posX:x,posY:y},"
		temp1="var nodes = ["

		for i in 1:num_packages
			temp1*=replace(replace(replace(nodetemp,"LABEL"=>string(i)),"x"=>positions[i,1]*800),"y"=>positions[i,2]*800)
		end
		temp1*="]"

		temp2="var edges = ["
		edgetemp="{from:I,to:J,color:D},"
		d=0
		colors=["\'#ff0000\'","\'#0000ff\'","\'#00ee00\'","\'#eeee00\'","\'#ee00ee\'","\'#00eeee\'","\'#ee8800\'","\'#88ee00\'","\'#8800ee\'","\'#0088ee\'"]

		for s in sequences
			d+=1
			for e in s
				temp2*=replace(replace(replace(edgetemp,"I"=>e[1]),"J"=>e[2]),"D"=>colors[d])
			end
		end
		temp2*="]"

		s = read(joinpath(dirname(@__DIR__),"visualise","vrp.html"), String)
		filename=joinpath("vrp"*string(Int(round((time()*100)%100000)))*".html")
		file=open(filename,"w")
		println(file,replace(s,"SOLUTION" => temp1*"\n"*temp2))
		close(file)

		if Sys.iswindows()
			run(`$(ENV["COMSPEC"]) /c start $(filename)`)
		elseif Sys.isapple()
			run(`open $(filename)`)
		elseif Sys.islinux() || Sys.isbsd()
			run(`xdg-open $(filename)`)
		else
			error("Unable to show plot. Try opening the file $(filename) manually.")
		end
		nothing
	end

	judy = JuDGEModel(mytree, ConditionallyUniformProbabilities, sub_problems, JuDGE_MP_Solver)
	initialise_columns(judy)

	h=nothing
	cb=nothing
	if use_heuristic
		h=heuristic
	end
	if use_branch_callback
		cb=bp_callback
	end
	judy=JuDGE.branch_and_price(judy,termination=Termination(reltol=reltol,rlx_reltol=rlx_reltol,inttol=10^-6),bp_callback=cb,branch_method=VRP_branch,max_no_int=-10,heuristic=h)
	show_vrp(judy)

	JuDGE.get_objval(judy)
end

@time @test VRP(seed=648,num_packages=20,num_drivers=3,max_deliveries=8,use_heuristic=false) ≈ -103.43628 atol = 1e-3
@time @test VRP(seed=648,num_packages=20,num_drivers=3,max_deliveries=8,use_heuristic=true) ≈ -103.43628 atol = 1e-3
@time @test VRP(seed=648,num_packages=20,num_drivers=3,max_deliveries=8,use_heuristic=true,use_branch_callback=true,rlx_reltol=0.15) ≈ -103.43628 atol = 1e-4

@test VRP(seed=120,num_packages=30,num_drivers=4,max_deliveries=8,use_heuristic=true,reltol=0.005,rlx_reltol=0.002) ≈ -130.8 atol = 1e-1