text:
@testset "getsubfield" begin
    type_field = Dict("type" => "type_field")
    @test GraphQLClient.getsubfield(type_field) == "type_field"
    ofType_field = Dict("type" => "ofType_field")
    @test GraphQLClient.getsubfield(ofType_field) == "ofType_field"
    not_a_field = Dict()
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.getsubfield(not_a_field)
end

@testset "is tests" begin
    # istype
    generic_field = Dict("type" => Dict("kind" => "CUSTOM_TYPE"))
    @test GraphQLClient.istype(generic_field, "CUSTOM_TYPE")
    @test !GraphQLClient.istype(generic_field, "NOT_CUSTOM_TYPE")
    generic_field = Dict("ofType" => Dict("kind" => "CUSTOM_TYPE"))
    @test GraphQLClient.istype(generic_field, "CUSTOM_TYPE")
    @test !GraphQLClient.istype(generic_field, "NOT_CUSTOM_TYPE")

    # specific, not nested
    object = Dict("type" => Dict("kind" => "OBJECT"))
    @test GraphQLClient.isobject(object)
    non_null = Dict("type" => Dict("kind" => "NON_NULL"))
    @test GraphQLClient.isnonnull(non_null)
    list = Dict("type" => Dict("kind" => "LIST"))
    @test GraphQLClient.islist(list)
    scalar = Dict("type" => Dict("kind" => "SCALAR"))
    @test GraphQLClient.isscalar(scalar)
    input_object = Dict("type" => Dict("kind" => "INPUT_OBJECT"))
    @test GraphQLClient.isinputobject(input_object)
    enum = Dict("type" => Dict("kind" => "ENUM"))
    @test GraphQLClient.isenum(enum)

    # specific and nested
    nonnull_input_object = Dict(
        "type" => Dict(
            "kind" => "NON_NULL",
            "ofType" => Dict("kind" => "INPUT_OBJECT")
        )
    )
    @test GraphQLClient.is_nonnull_input_object(nonnull_input_object)
    not_a_nonnull_input_object = Dict(
        "type" => Dict(
            "kind" => "NON_NULL",
            "ofType" => Dict("kind" => "SCALAR")
        )
    )
    @test !GraphQLClient.is_nonnull_input_object(not_a_nonnull_input_object)
    @test GraphQLClient.isroottypescalar(not_a_nonnull_input_object)
    listofobjects = Dict(
        "type" => Dict(
            "kind" => "LIST",
            "ofType" => Dict("kind" => "OBJECT")
        )
    )
    @test GraphQLClient.isroottypeobject(listofobjects)
    @test !GraphQLClient.isroottypeenum(listofobjects)
    listofenums = Dict(
        "type" => Dict(
            "kind" => "LIST",
            "ofType" => Dict("kind" => "ENUM")
        )
    )
    @test GraphQLClient.isroottypeenum(listofenums)
    @test !GraphQLClient.isroottypeobject(listofenums)
    nonnull_enum_list = build_type("kind", "Fieldname",
        type=build_type("NON_NULL", "",
            ofType=build_type("LIST", "",
                ofType=build_type("ENUM",""))))
    @test GraphQLClient.isroottypeenum(nonnull_enum_list)
    nonnull_object_list = build_type("kind", "Fieldname",
        type=build_type("NON_NULL", "",
            ofType=build_type("LIST", "",
                ofType=build_type("OBJECT",""))))
    @test GraphQLClient.isroottypeobject(nonnull_object_list)
end

@testset "get_field_type_string and getjuliatype" begin
    # SCALAR
    arg = build_arg("SCALAR", nothing, "String", nothing)
    @test GraphQLClient.get_field_type_string(arg) == "String"
    @test GraphQLClient.getroottype(arg) == "String"
    @test GraphQLClient.getjuliatype(arg) == String

    # Custom scalar
    scalar_dict = merge(
        GraphQLClient.GQL_DEFAULT_SCALAR_TO_JULIA_TYPE,
        Dict("CUSTOM_FLOAT_32" => Float32),
    )
    arg = build_arg("SCALAR", nothing, "CUSTOM_FLOAT_32", nothing)
    @test GraphQLClient.get_field_type_string(arg) == "CUSTOM_FLOAT_32"
    @test GraphQLClient.getroottype(arg) == "CUSTOM_FLOAT_32"
    @test GraphQLClient.getjuliatype(arg; scalar_types=scalar_dict) == Float32

    # INPUT_OBJECT
    arg = build_arg("INPUT_OBJECT", nothing, "MyObject", nothing)
    @test GraphQLClient.get_field_type_string(arg) == "MyObject"
    @test GraphQLClient.getroottype(arg) == "MyObject"
    @test GraphQLClient.getjuliatype(arg) == GraphQLClient.InputObject

    # ENUM
    arg = build_arg("ENUM", nothing, "MyEnum", nothing)
    @test GraphQLClient.get_field_type_string(arg) == "MyEnum"
    @test GraphQLClient.getroottype(arg) == "MyEnum"
    @test GraphQLClient.getjuliatype(arg) == String

    # NON_NULL
    arg = Dict(
        "name" => "MyArg",
        "defaultValue" =>  nothing,
        "description" =>  nothing,
        "type" => build_type("NON_NULL", nothing; ofType=build_type("SCALAR", "String"))
    )
    @test GraphQLClient.get_field_type_string(arg) == "String!"
    @test GraphQLClient.getroottype(arg) == "String"
    @test GraphQLClient.getjuliatype(arg) == String

    # LIST
    arg = Dict(
        "name" => "MyArg",
        "defaultValue" =>  nothing,
        "description" =>  nothing,
        "type" => build_type("LIST", nothing; ofType=build_type("SCALAR", "Float"))
    )
    @test GraphQLClient.get_field_type_string(arg) == "[Float]"
    @test GraphQLClient.getroottype(arg) == "Float"
    @test GraphQLClient.getjuliatype(arg) == Vector{Float64}

    # NESTED NON-NULL LIST
    arg = Dict(
        "name" => "MyArg",
        "defaultValue" =>  nothing,
        "description" =>  nothing,
        "type" => build_type("NON_NULL", nothing;
            ofType=build_type("LIST", nothing,
                ofType=build_type("LIST", nothing,
                    ofType=build_type("ENUM", "MyEnum")))))
    @test GraphQLClient.get_field_type_string(arg) == "[[MyEnum]]!"
    @test GraphQLClient.getroottype(arg) == "MyEnum"
    @test GraphQLClient.getjuliatype(arg) == Vector{Vector{String}}

    # OBJECT
    arg = build_arg("OBJECT", nothing, "MyObject", nothing)
    @test GraphQLClient.get_field_type_string(arg) == "MyObject"
    @test GraphQLClient.getroottype(arg) == "MyObject"
    @test GraphQLClient.getjuliatype(arg) == GraphQLClient.Object

    # Errors - not handled types
    arg = build_arg("UNION", nothing, "MyEnum", nothing)
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.get_field_type_string(arg)
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.getjuliatype(arg)
    arg = build_arg("INTERFACE", nothing, "MyEnum", nothing)
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.get_field_type_string(arg)
    @test_throws GraphQLClient.GraphQLClientException GraphQLClient.getjuliatype(arg)

    # Errors - unrecognised scalars
    arg = build_arg("SCALAR", nothing, "NOT_A_DEFAULT", nothing)
    @test_throws ArgumentError GraphQLClient.getjuliatype(arg)
    arg = build_arg("SCALAR", nothing, "NOT_IN_CUSTOM_DICT", nothing)
    @test_throws ArgumentError GraphQLClient.getjuliatype(arg; scalar_types=scalar_dict)
end

@testset "getroottypefield" begin
    type = build_type("NON_NULL", "name",
        type=build_type("LIST", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("SCALAR", "Boolean"))))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
    type = build_type("LIST", "name",
        type=build_type("LIST", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("SCALAR", "Boolean"))))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
    type = build_type("LIST", "name", type=build_type("SCALAR", "Boolean"))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
end

@testset "getroottypefield" begin
    type = build_type("NON_NULL", "name",
        type=build_type("LIST", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("SCALAR", "Boolean"))))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
    type = build_type("LIST", "name",
        type=build_type("LIST", nothing,
            ofType=build_type("LIST", nothing,
                ofType=build_type("SCALAR", "Boolean"))))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
    type = build_type("LIST", "name", type=build_type("SCALAR", "Boolean"))
    @test GraphQLClient.getroottypefield(type) == build_type("SCALAR", "Boolean")
end


text:
#
# Extracting a slot from the environment.
#

SlotQuery(tag::Symbol, oty::Output) =
    Query(
        SlotSig(tag),
        Input(Any) |> setslots([InputSlot(tag, oty)]),
        oty)

SlotQuery(tag::Symbol, oty) =
    SlotQuery(tag, convert(Output, oty))

SlotQuery(ity::Input, tag::Symbol) =
    SlotQuery(tag, slot(ity, tag))

immutable SlotSig <: AbstractPrimitive
    tag::Symbol
end

function ev(sig::SlotSig, ::Input, oty::Output, iflow::InputFlow)
    for (name, flow) in slotflows(iflow)
        if name == sig.tag
            return flow
        end
    end
end

ext: jl

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
export SUNIrrep

# is this type piracy?
Base.getindex(::TensorKit.IrrepTable, ::Type{TensorKit.SU{N}}) where {N} = SUNIrrep{N}
Base.convert(::Type{SUNIrrep{N}}, I::NTuple{N,Int}) where N = SUNIrrep{N}(I)
Base.IteratorSize(::Type{TensorKit.SectorValues{T}}) where T<:SUNIrrep = Base.IsInfinite()

function Base.iterate(::TensorKit.SectorValues{SUNIrrep{N}}, I = ntuple(zero, N)) where N
    s = SUNIrrep(I)
    k = N-1
    while k > 1 && I[k] == I[k-1]
        k -= 1
    end
    I = Base.setindex(I, I[k]+1, k)
    for l = k+1:N-1
        I = Base.setindex(I, 0, l)
    end
    return s, I
end

Base.:(==)(s::SUNIrrep, t::SUNIrrep) = ==(s.I, t.I)
Base.hash(s::SUNIrrep, h::UInt) = hash(s.I, h)
Base.conj(s::SUNIrrep) = SUNIrrep(s.I[1] .- reverse(s.I))
Base.one(::Type{SUNIrrep{N}}) where N = SUNIrrep(ntuple(n->0, N))

TensorKit.FusionStyle(::Type{SUNIrrep{N}}) where N = TensorKit.GenericFusion()
Base.isreal(::Type{<:SUNIrrep}) = true
TensorKit.BraidingStyle(::Type{<:SUNIrrep}) = TensorKit.Bosonic();

TensorKit.:⊗(s1::SUNIrrep{N}, s2::SUNIrrep{N}) where N =
    TensorKit.SectorSet{SUNIrrep{N}}( keys(directproduct(s1,s2)) )
TensorKit.Nsymbol(s1::SUNIrrep{N}, s2::SUNIrrep{N}, s3::SUNIrrep{N}) where N =
    get(directproduct(s1, s2), s3, 0)

TensorKit.fusiontensor(s1::SUNIrrep{N}, s2::SUNIrrep{N}, s3::SUNIrrep{N}) where N =
    CGC(Float64, s1, s2, s3)

const FCACHE = Vector{Any}()
function TensorKit.Fsymbol(a::SUNIrrep{N}, b::SUNIrrep{N}, c::SUNIrrep{N},
                            d::SUNIrrep{N}, e::SUNIrrep{N}, f::SUNIrrep{N}) where N
    key = (a, b, c, d, e, f)
    K = typeof(key)
    V = Array{Float64,4}
    if length(FCACHE) < N || !isassigned(FCACHE, N)
        resize!(FCACHE, max(length(FCACHE), N))
        cache = Dict{K,V}()
        FCACHE[N] = cache
    else
        cache::Dict{K,V} = FCACHE[N]
    end
    return get!(cache, key) do
        _Fsymbol(a,b,c,d,e,f)
    end
end
function _Fsymbol(a::SUNIrrep{N}, b::SUNIrrep{N}, c::SUNIrrep{N},
                            d::SUNIrrep{N}, e::SUNIrrep{N}, f::SUNIrrep{N}) where N
    N1 = Nsymbol(a,b,e)
    N2 = Nsymbol(e,c,d)
    N3 = Nsymbol(b,c,f)
    N4 = Nsymbol(a,f,d)

    (N1 == 0 || N2 == 0 || N3 == 0 || N4 == 0) && return fill(0.0, N1, N2, N3, N4)

    # computing first diagonal element
    A = fusiontensor(a,b,e)
    B = fusiontensor(e,c,d)[:, :, 1, :]
    C = fusiontensor(b,c,f)
    D = fusiontensor(a,f,d)[:, :, 1, :]

    @tensor F[-1,-2,-3,-4] := conj(D[1,5,-4]) * conj(C[2,4,5,-3]) *
                                A[1,2,3,-1] * B[3,4,-2]
    return Array(F)
end

const RCACHE = Vector{Any}()
function TensorKit.Rsymbol(a::SUNIrrep{N}, b::SUNIrrep{N}, c::SUNIrrep{N}) where N
    key = (a, b, c)
    K = typeof(key)
    V = Array{Float64,2}
    if length(RCACHE) < N || !isassigned(RCACHE, N)
        resize!(RCACHE, max(length(RCACHE), N))
        cache = Dict{K,V}()
        RCACHE[N] = cache
    else
        cache::Dict{K,V} = RCACHE[N]
    end
    return get!(cache, key) do
        _Rsymbol(a,b,c)
    end
end
function _Rsymbol(a::SUNIrrep{N}, b::SUNIrrep{N}, c::SUNIrrep{N}) where N
    N1 = Nsymbol(a, b, c);
    N2 = Nsymbol(b, a, c);

    (N1 == 0 || N2 == 0) && return fill(0.0, N1, N2)

    A = fusiontensor(a, b, c)[:, :, 1, :]
    B = fusiontensor(b, a, c)[:, :, 1, :]

    @tensor R[-1;-2] := conj(B[1, 2, -2]) * A[2, 1, -1]
    Array(R)
end


text:
## A simple token-based interface

"""
    istokenizable(indices::AbstractIndices)

Return `true` if the indices `indices` obeys the token interface, or `false` otherwise.

A token is a more efficient way of refering to an element of `indices`. Using tokens may
help avoid multiple index lookups for a single operation.

A tokenizable indices must implement:

 * `tokentype(indices) --> T::Type`
 * `iteratetoken(indices, s...)` iterates the tokens of `indices`, like `iterate`
 * `gettoken(indices, i) --> (hasindex::Bool, token)`
 * `gettokenvalue(indices, token)` returning the value of the index at `token`

An `isinsertable` tokenizable indices must implement

 * `gettoken!(indices, i) --> (hadtoken::Bool, token)`
 * `deletetoken!(indices, token) --> indices`

    istokenizable(dict::AbstractDictionary)

Return `true` if the dictionary `dict` obeys the token interface, or `false` otherwise.

A token is a more efficient way of refering to an element of `dict`. Using tokens may
help avoid multiple index lookups for a single operation.

An tokenizable dictionary must implement:

 * `keys(dict)` must be `istokenizable` and share tokens with `dict`
 * `gettokenvalue(dict, token)` returning the dictionary value at `token`
 * `istokenassigned(dict, token) --> Bool`

An `issettable` tokenizable dictionary must implement:

 * `settokenvalue!(dict, token)`

An `isinsertable` tokenizable dictionary must implement:

 * `gettoken!(dict, i) --> (hadtoken::Bool, token)`
 * `deletetoken!(dict, token) --> dict`
"""
istokenizable(d::AbstractIndices) = false
istokenizable(d::AbstractDictionary) = istokenizable(keys(d))

"""
    tokentype(dict)

Return the type of the tokens retrieved via `gettoken` for dictionary `dict`, for
dictionaries that are `istokenizable`.
"""
tokentype(d::AbstractDictionary) = tokentype(keys(d))

function tokentype(inds::AbstractIndices)
    if istokenizable(inds)
        error("`tokentype` needs to be defined for tokenizable indices: $(typeof(inds))")
    else
        error("Indices are not tokenizable: $(typeof(inds))")
    end
end

"""
    IndicesTokens(indices)

Return a dictionary mapping from `i ∈ indices` to a valid token.
"""
struct IndicesTokens{I, T, Parent <: AbstractIndices{I}} <: AbstractDictionary{I, T}
    parent::Parent
end

IndicesTokens(indices::AbstractIndices{I}) where {I} = IndicesTokens{I}(indices)
IndicesTokens{I}(indices::AbstractIndices{I}) where {I} = IndicesTokens{I, tokentype(indices)}(indices)

Base.parent(ts::IndicesTokens) = getfield(ts, :parent)

Base.keys(ts::IndicesTokens) = parent(ts)

# These containers are obviously tokenizable (but never settable or insertable)
istokenizable(::IndicesTokens) = true
tokentype(::IndicesTokens{<:Any, T}) where {T} = T
tokens(ts::IndicesTokens) = ts
istokenassigned(ts::IndicesTokens, t) = istokenassigned(parent(ts), t)
@propagate_inbounds iteratetoken(ts::IndicesTokens, s...) = iteratetoken(parent(ts), s...)
@propagate_inbounds gettoken(ts::IndicesTokens{I}, i::I) where {I} = gettoken(parent(ts), i)
gettokenvalue(ts::IndicesTokens, t) = t

Base.IteratorSize(ts::IndicesTokens) = Base.IteratorSize(parent(ts))
Base.length(ts::IndicesTokens) = length(parent(ts))

"""
    tokens(dict::AbstractDictionary)

Return a new dictionary mapping from the `keys` of `dict` to a "token". The token can be
used to fetch a value with `gettokenvalue`. For mutable containers, the token can be used to
overwrite a value with `settokenvalue!`.
"""
@inline tokens(d::AbstractDictionary) = tokens(keys(d))
@inline function tokens(i::AbstractIndices{I}) where {I}
    if istokenizable(i)
        return IndicesTokens{I, tokentype(i), typeof(i)}(i)
    else
        return i
    end
end

"""
    Tokenized(dict)

Return an indexable container mapping from the `token ∈ tokens(dict)` to the values of
`dict`.
"""
struct Tokenized{Parent <: AbstractDictionary}
    parent::Parent
end

Base.parent(dict::Tokenized) = getfield(dict, :parent)

@propagate_inbounds function Base.getindex(d::Tokenized, t)
    return gettoken(parent(d), t)
end

function Base.isassigned(d::Tokenized, t)
    return istokenassigned(parent(d), t)
end

@propagate_inbounds function Base.setindex!(d::Tokenized, value, t)
    return settokenvalue!(parent(d), t, value)
end

"""
    tokenized(dict::AbstractDictionary)

For `istokenizable` dictionary `dict`, return a container that supports
`getindex(dict, token)` and `isassigned(dict, token)`, where `token ∈ tokens(dict)`.

If `issettable(dict)` then the result should also support `setindex!(dict, value, token)`.

Note: the output is not necessarily an `AbstractDictionary`.
"""
function tokenized(d::AbstractDictionary)
    if istokenizable(d)
        return Tokenized{tokentype(d), eltype(d), typeof(d)}(d)
    else
        return d
    end
end

"""
    gettoken(dict, i)

Return the tuple `(hasindex, token)`, where `hasindex` is `true` if `i` exists in `dict`.
The `token` can be used to retrieve a value using the `gettokenvalue` function. You can
check if a token is assigned to a valid Julia object (i.e. not `#undef`) via
`istokenassigned`.

Settable (i.e. mutable) dictionaries allow you to set a corresponding value via the
`settokenvalue!` function (see the `issettable` trait).

Insertable dictionaries provide the `gettoken!` function (see the `isinsertable` trait).
"""
@propagate_inbounds function gettoken(inds::AbstractIndices{I}, i::I) where I
    if istokenizable(inds)
        error("gettoken needs to be defined for tokenizable indices: $(typeof(inds))")
    end

    @boundscheck if i ∉ inds
        return (false, i)
    end
    return (true, i)
end

@propagate_inbounds function gettoken(d::AbstractDictionary{I}, i::I) where I
    return gettoken(keys(d), i)
end

@propagate_inbounds function gettoken(d::AbstractDictionary{I}, i) where I
    return gettoken(d, convert(I, i))
end

@propagate_inbounds function gettokenvalue(d::AbstractDictionary, token)
    if istokenizable(d)
        error("gettokenvalue needs to be defined for tokenizable $(d isa AbstractIndices ? "indices" : "dictionary"): $(typeof(d))")
    end

    return d[token]
end

@propagate_inbounds function istokenassigned(d::AbstractDictionary, token)
    if istokenizable(d)
        error("istokenassigned needs to be defined for tokenizable $(d isa AbstractIndices ? "indices" : "dictionary"): $(typeof(d))")
    end

    return isassigned(d, token)
end

function istokenassigned(::AbstractIndices, token)
    return true  # We are checking that given a valid token a value exists
end

@propagate_inbounds function settokenvalue!(d::AbstractDictionary{<:Any,T}, t, value) where {T}
    return settokenvalue!(d, t, convert(T, value))
end

@propagate_inbounds function settokenvalue!(d::AbstractDictionary{<:Any,T}, i, value::T) where {T}
    if !issettable(d)
        error("Cannot mutate values of dictionary: $(typeof(d))")
    end
    if istokenizable(d)
        error("settoken! needs to be defined for settable, tokenizable dictionary: $(typeof(d))")
    end

    return d[i] = value
end

## Check if we can do fast mutual iteration over multiple containers

# I guess someone could overload this method if they can check this faster than O(N),
# but more precise than ===

"""
    sharetokens(dict1, dict2)

Return `true` if `dict1` and `dict2` obviously share tokens, using a test which can be
performed quickly (e.g. O(1) rather than O(N)). Return `false` otherwise.

Note: the test may not be precise, this defaults to `tokens(dict1) === tokens(dict2)`.
"""
sharetokens(i1::AbstractIndices, i2::AbstractIndices) = istokenizable(i1) && i1 === i2
sharetokens(d1, d2) = sharetokens(keys(d1), keys(d2))
sharetokens(d1, d2, ds...) = sharetokens(d1, d2) && sharetokens(d1, ds...)


text:
"""

"""
struct InfiniteOptMappingData
    ref_to_num::Dict{InfiniteOpt.GeneralVariableRef, _Symbolics.Num}
    refs_to_arr::Dict{Vector{InfiniteOpt.GeneralVariableRef}, _Symbolics.Arr{Num, 1}}
    param_to_diff::Dict{InfiniteOpt.GeneralVariableRef, _Symbolics.Differential}

    # default constructor
    function InfiniteOptMappingData()
        return new(
            Dict{InfiniteOpt.GeneralVariableRef, _Symbolics.Num}(),
            Dict{Vector{InfiniteOpt.GeneralVariableRef}, _Symbolics.Arr{Num, 1}}(),
            Dict{InfiniteOpt.GeneralVariableRef, _Symbolics.Differential}()
        )
    end
end

# Extract the root name of a variable (removes the bracketed container indices)
function _remove_name_index(name::String)
    first_bracket = findfirst(isequal('['), name)
    if isnothing(first_bracket)
        return name
    else
        return name[1:prevind(name, first_bracket, 1)] # can handle unicode
    end
end

# convert the brackets into a compatible format
function _format_name(name::String)
    first_bracket = findfirst(isequal('['), name)
    if first_bracket == 1
        error("Cannot create a `Symbolics.jl` representation using ",
              "anonymous parameters/variables.")
    elseif isnothing(first_bracket)
        return name
    else
        last_root_idx = prevind(name, first_bracket, 1)
        root = name[1:last_root_idx]
        inds_name = replace(name[first_bracket+1:end-1], "," => "_")
        return string(root, "_", inds_name)
    end
end

"""

"""
function create_symbolic_params(
    model::InfiniteOpt.InfiniteModel,
    data::InfiniteOptMappingData
    )
    # handle the scalar parameters
    for p in InfiniteOpt.all_parameters(model, ScalarParameter)
        name = _formate_name(_JuMP.name(p))
        data.ref_to_num[p] = first(ModelingToolkit.@parameters($(Symbol(name))))
    end
    # handle the multi-variate infinite parameters
    for (idx, data) in model.dependent_parameters
        name = _remove_name_index(first(data.names))
        if isempty(name) || !_all_equal(_remove_name_index.(data.names))
            error("Cannot create a `Symbolics.jl` representation using ",
                  "anonymous infinite parameters or ones with heterogeneous ",
                  "naming.")
        end
        num_ps = length(data.names)
        ps = [InfiniteOpt.GeneralVariableRef(model, idx.value, InfiniteOpt.DependentParameterIndex, i)
              for i in 1:num_ps]
        data.refs_to_arr[ps] = first(ModelingToolkit.@parameters($(Symbol(name))[1:num_ps]))
        for (i, p) in enumerate(ps)
            data.ref_to_num[p] = data.refs_to_arr[ps][i]
        end
    end
    return
end

"""

"""
function create_derivative_operators(
    model::InfiniteOpt.InfiniteModel,
    data::InfiniteOptMappingData
    )
    deriv_params = unique!([p[2] for p in keys(model.deriv_lookup)])
    for p in deriv_params
        data.param_to_diff[p] = _Symbolics.Differential(data.ref_to_num[p])
    end
    return
end

## Helper functions to retrieve infinite parameter mappings
# Single parameter
function _map_param_group(
    p::InfiniteOpt.GeneralVariableRef,
    data::InfiniteOptMappingData
    )
    return data.ref_to_num[p]
end

# Multi-variate parameter
function _map_param_group(
    ps::Vector{InfiniteOpt.GeneralVariableRef},
    data::InfiniteOptMappingData
    )
    return data.refs_to_arr[ps]
end

"""

"""
function create_symbolic_vars(
    model::InfiniteModel,
    data::InfiniteOptMappingData
    )
    # finite variables
    for v in _JuMP.all_variables(model, FiniteVariable)
        name = _format_name(_JuMP.name(v))
        data.ref_to_num[v] = first(_Symbolics.@variables($(Symbol(name))))
    end
    # infinite variables
    for v in _JuMP.all_variables(model, InfiniteVariable)
        name = _format_name(_JuMP.name(v))
        param_tuple = Tuple(InfiniteOpt.raw_parameter_refs(v), use_indices = false)
        num_tuple = _map_param_group.(param_tuple)
        data.ref_to_num[v] = first(_Symbolics.@variables($(Symbol(name))($(num_tuple...))))
    end
    # derivatives
    for d in InfiniteOpt.all_derivatives(model)
        name = _format_name(_JuMP.name(d))

    end
    # point variables

    # semi-infinite variables

end


text:
### functions for posterior odds and Bayes p-value
#
"""
    todds(theta_hat,theta_hat_se,v)

Input:

        v = degrees of freedom

        theta_hat, theta_hat_se
        = estimate and standard error for theta

Optional argument:

        h0 = value in hull hypothesis (default = 0).

Returns:

        t-posterior odds ratio of evidence against h0 (default = 0).
"""
function todds(theta_hat, theta_se, v; h0 = 0.0)
  # compute posterior odds assuming a Student-t distribution
  # mcs = MC sample
  # v = degrees of freedom, n - k
  # h0 = value under null hypothesis
  that = abs.(theta_hat - h0)/theta_se
  odds = (1.0 + (that^2)/v)^((v+1)/2.0)
end

"""
    mcodds(theta_draws)

Input:

        theta_draws = MC sample from posterior for theta

Optional argument:

        h0 = value in hull hypothesis (default = 0).

Returns:

        posterior odds ratio of evidence against h0.
"""
function mcodds(mcs; h0=0.0)
  # dependency: KernelDensity
  # compute posterior odds from nonparametric density
  # mcs = MC sample

### NEED catches for when odds too small/large - out of bounds of sample
  postden = KernelDensity.kde(mcs)
  ax = abs.(postden.x .- h0)
  indx = findall(ax .== minimum(ax))  # find all x that = 4 in x
  numodds = maximum(postden.density)
  denodds = postden.density[indx[1]]
  odds = numodds/denodds[1]
  if odds > 5000000.0
    odds = 5000000.0
  end
  return odds
end

# Example: if bs[:,2] is vector of MCMC draws for beta
# to test if beta = 1.00:
#mcodds(bs[:,2],h0=1.00)

"""
    pdr_pval(theta_draws)

Input:

        theta_draws = MC sample from posterior for theta

Optional argument:

        h0 = value in hull hypothesis (default = 0).

Returns:

        posterior density ratio, one-tailed and two-tailed 'p-value' of evidence against h0.
        one-tailed 'p-value' = minimum[Prob(theta < 0), Prob(theta>0)]
"""
# x = MC sample, h0 = null hypothesis value
# usage: odds, p_val, pval_2 = post_odds_pval(x, h0 = -2.3)
function pdr_pval(x; h0 = 0.0)
    d = kde(x)
    p,ind = findmax(d.density)
    ph0 = ifelse(pdf(d,h0) == 0.0,0.0000001,pdf(d,h0))
    odds = pdf(d,d.x[ind])/ph0
    p_val = length(x[x .<= h0])
    p_val2 = length(x[x .>= h0])
    if p_val <= p_val2
        p_value = p_val/length(x)
    else
        p_value = p_val2/length(x)
    end
    p_value_2tail = 2*p_value
    return odds, p_value, p_value_2tail
end

"""
    bayespval(theta_draws)

Input:

        theta_draws = MC sample from posterior for theta

Optional argument:

        h0 = value in hull hypothesis (default = 0).

Returns:

        posterior 'p-value' of evidence against h0.
        pval = one tailed value minimum[(Prob(theta <0), Prob(theta>0)]
        pval2 = two-tailed value (2*pval)
"""
function bayespval(mcs; h0=0.0)
  p_val = length(mcs[mcs .<= h0])/length(mcs)
  if p_val == 0.0
    p_val = 0.00001
  end
  if p_val > 0.5
    p_val = (1.0 - p_val)
  end
  p_val2 = 2.0*p_val
  return p_val, p_val2
end

#= old deprecated todds assuming t
function todds(mcs,n,k; h0 = 0.0)
  # compute posterior odds assuming a Student-t distribution
  # mcs = MC sample
  # v = degrees of freedom, n - k
  # h0 = value under null hypothesis
  v = n-k
  unbiasedseb = std(mcs)*sqrt(v)/sqrt(n)
  that = abs.(mean(mcs) - h0)/unbiasedseb
  odds = (1 + (that^2)/v)^((v+1)/2)
end
=#


text:
@testset "Interface" begin
    @test TPCA(5) == TPCA(ncomps=5)
    @test isnothing(TPCA().ncomps)
end

@testset "Decomposition" begin
    data = 4 .* randn(rng, 30, 101, 101) .+ 10
    angles = sort!(90randn(rng, 30)) |> normalize_par_angles

    # get sizes correct for ncomps
    for N in [1, 3, 5]
        # A, w = @inferred decompose(TPCA(N), data, angles) TODO type stability
        A, w = decompose(TPCA(N), data, angles)
        @test size(A) == (N, 101 * 101)
        @test size(w) == (30, N)
    end
    @test_throws ErrorException decompose(TPCA(40), data, angles)

    # default is to use whole cube
    S = reconstruct(TPCA(), data, angles)
    @test size(S) == (30, 101, 101)
end

# @testset "ADI Trivial" begin
#     cube = ones(10, 100, 100)
#     angles = zeros(10)

#     reduced_5 = TPCA(5)(cube, angles)
#     reduced_10 = TPCA(10)(cube, angles)

#     @test size(reduced_5) == size(reduced_10) == (100, 100)
#     @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_5)
#     @test all(x -> isapprox(x, 0, atol = 1e-9), reduced_10)
# end

# @testset "RDI Trivial" begin
#     data = randn(rng, 30, 101, 101)
#     angles = sort!(90rand(rng, 30)) |> normalize_par_angles

#     S = data .- reconstruct(TPCA(2), data, angles, zeros(10, 101, 101))
#     @test S ≈ data rtol=2e-2
# end


text:
#Comparing Gini Index to L0 as a measure of sparsity

using PyPlot
using CompressedSensing
srand(0)

#Create 2 data sets, one without noise, and one with 5% gaussian noise
example_data_no_noise=zeros(100)
example_data_noise=randn(100)*.05

#Create arrays that will contain the measures of sparsity for the two data sets
calculated_GI_no_noise = []
calculated_GI_noise = []
calculated_L0_no_noise = []
calculated_L0_noise = []

#for both data sets, set an increasing number of points to 1, IE: reduce sparsity on element at a time
for i=1:length(example_data_no_noise)

    example_data_no_noise[i]+=1
    example_data_noise[i]+=1

    #calculate the Gini Index for both data sets and save the value
    calculated_GI_no_noise = [calculated_GI_no_noise,GI(example_data_no_noise)]
    calculated_GI_noise = [calculated_GI_noise,GI(example_data_noise)]

    #calculate the L0 sparsity measure for each data set and save the value
    calculated_L0_no_noise = [calculated_L0_no_noise,sum(abs(example_data_no_noise).>0)]
    calculated_L0_noise = [calculated_L0_noise,sum(abs(example_data_noise).>0.)]
end

#Plot the results
figure(figsize=(6,5))
title("Comparison of Gini Index to L0 as measures of Sparsity")
plot([1:length(example_data_no_noise)],calculated_GI_no_noise,label="Gini Index",color="red")
plot([1:length(example_data_noise)],calculated_GI_noise,label="Gini Index w/ Noise",color="purple")
xlabel("Number of Non-Zero Elements in Data")
ylabel("Gini Index")
ylim([0,1.05])
grid()
#put the L0 on a different axes, as it goes from 0-100 not 0-1
twinx()
plot([1:length(example_data_no_noise)],calculated_L0_no_noise,label="L0 Sparsity")
plot([1:length(example_data_noise)],calculated_L0_noise,label="L0 w/ Noise")
ylabel("L0 Sparsity Measure")
ylim([0,105])
annotate("L0 w/ Noise", xy=[70, 100], xytext=(65, 89),arrowprops=["width"=>.1,"headwidth"=>5]);
annotate("L0", xy=[80, 80], xytext=(85, 69),arrowprops=["width"=>.1,"headwidth"=>5]);
annotate("GI", xy=[20, 80], xytext=(25, 90),arrowprops=["width"=>.1,"headwidth"=>5]);
annotate("GI w/ Noise", xy=[20, 71], xytext=(10, 50),arrowprops=["width"=>.1,"headwidth"=>5]);
savefig("Example_3_Fig_1.png")

max_stars_repo_path: examples/Example 3.jl

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_path: examples/Example 3.jl

max_forks_repo_path: examples/Example 3.jl


text:
abstract type AbstractDomain end

abstract type AbstractSphericalDomain <: AbstractDomain end

abstract type AbstractPlanarDomain <: AbstractDomain end

"""
    SphericalTriangle((α::arb, β::arb, γ::arb)[, RR::ArbField = RealField(64)])
    SphericalTriangle((α::fmpq, β::fmpq, γ::fmpq)[, RR::ArbField = RealField(64)])
> Create a spherical triangle with angles α, β, γ based on the field
  RR.

  If the angles are of type arb then α, β and γ represents the angles
  directly. If they are of type fmpq then they represent the angle as
  a rational multiple of π. The vertex with angle α is places on the
  north pole and the vertex with angle β is taken to have y = 0.
  Vertex 1, 2 and 3 correspond to angles α, β and γ respectively. The
  boundaries are ordered by which vertex they are opposite of, so
  boundary 1 is opposite of vertex 1.
"""
struct SphericalTriangle{T<:Union{fmpq,arb}} <: AbstractSphericalDomain
    angles::Tuple{T,T,T}
    parent::ArbField
end

"""
    LShape()
> Create the classical L-shaped domain.

  The only relevant vertex is the re-entrant one with angle 3π/2 and
  the boundaries are enumerated from 1 to 4 in positive direction
  skipping the two which are directly adjacent to the re-entrant
  vertex.
"""
struct LShape <: AbstractPlanarDomain
    parent::ArbField
end

"""
    Triangle(α, β[, parent::ArbField = RealField(64)])

Create a (planar) triangle with two of the angles being `α` and `β`.
The length of the side in between is defined to be one and the last
angle is called `γ`.

If the angles are of type `arb` then `α` and `β` represents the angles
directly. If they are of type `fmpq` then they represent the angle as
a rational multiple of `π`. The vertex with angle `α` is placed at the
origin the vertex with angle `β` is taken to have `y = 0`.

Vertex 1, 2 and 3 are opposite of angles α, β and γ respectively.
"""
struct Triangle{T<:Union{fmpq,arb}} <: AbstractPlanarDomain
    angles::NTuple{3,T}
    parent::ArbField

    function Triangle(α::fmpq, β::fmpq, parent::ArbField = RealField(64))
        α > 0 || throw(DomainError(α, "angle must be positive"))
        β > 0 || throw(DomainError(β, "angle must be positive"))
        α + β < 1 || throw(ArgumentError("α + β must be less than 1"))
        γ = 1 - α - β
        return new{fmpq}((α, β, γ), parent)
    end

    function Triangle(α::arb, β::arb, parent::ArbField = parent(α))
        α > 0 || throw(DomainError(α, "angle must be positive"))
        β > 0 || throw(DomainError(β, "angle must be positive"))
        !(α + β > parent(π)) || throw(ArgumentError("α + β must be less than π"))
        α + β < parent(π) || @warn "α + β might not be less than π"
        γ = parent(π) - α - β
        return new{arb}((α, β, γ), parent)
    end
end

"""
    Polygon(angles, vertices)

Create a (planar) polygon with the given angles. To simplify the
implementation we also require that the location. of the vertices are
given.

The boundaries are enumerated by which vertex they are next to, in
positive order.
"""
struct Polygon{T<:Union{fmpq,arb}} <: AbstractPlanarDomain
    angles::Vector{T}
    vertices::Vector{SVector{2,arb}}
    parent::ArbField
end

"""
    TransformedDomain{T<:AbstractPlanarDomain}(domain::T, rotation, scaling, translation)

Represents a transformation of `domain` corresponding of a rotation,
a uniform scaling and a translation (in that order).
"""
struct TransformedDomain{T<:Union{fmpq,arb},S<:AbstractPlanarDomain} <: AbstractPlanarDomain
    original::S
    rotation::T
    scaling::arb
    translation::SVector{2,arb}
end

"""
    IntersectedDomain{T}(exterior::T, interiors::Vector{AbstractPlanarDomain})

Represents a domain given by removing from `exterior` the domains in
`interiors`. It's assumed that all domains in `interiors` are
contained in `exterior` and that they are pairwise disjoint.

The boundaries are enumerated in a way such that first comes the
exterior domains boundaries and then comes the interior domains.
"""
struct IntersectedDomain{T<:AbstractPlanarDomain} <: AbstractPlanarDomain
    exterior::T
    interiors::Vector{AbstractPlanarDomain}

    function IntersectedDomain(
        exterior::T,
        interiors::Vector{<:AbstractPlanarDomain},
    ) where {T<:AbstractPlanarDomain}
        domain = new{T}(exterior, convert(Vector{AbstractPlanarDomain}, interiors))
    end
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Decision #
# ========================== #
@enum DecisionState NotTaken Taken Known

mutable struct Decision{T}
    state::DecisionState
    value::T
end
state(decision::Decision) = decision.state
decision_value(decision::Decision) = state(decision) == NotTaken ? NaN : decision.value

function Decision(value, ::Type{T}) where T
    return Decision(Taken, T(value))
end

function KnownDecision(value, ::Type{T}) where T
    return Decision(Known, T(value))
end

function Decision(info::JuMP.VariableInfo, ::Type{T}) where T
    if info.has_fix
        return Decision(Taken, T(info.fixed_value))
    end
    return Decision(NotTaken, T(NaN))
end

function KnownDecision(info::JuMP.VariableInfo, ::Type{T}) where T
    value = if info.has_fix
        # Value fixed at construction
        T(info.fixed_value)
    else
        # Value is known, but not yet set. Pick a feasible value for now
        value = if info.has_lb
            T(info.lower_bound)
        elseif info.has_ub
            T(info.upper_bound)
        else
            zero(T)
        end
    end
    return Decision(Known, value)
end

Base.copy(decision::Decision) = Decision(decision.state, decision.value)

# Decisions #
# ========================== #
struct IgnoreDecisions end

const DecisionMap = OrderedDict{MOI.VariableIndex, Decision{Float64}}
const StageMap = Dict{MOI.VariableIndex, Int}
const ConstraintMap = MOIU.DoubleDicts.MainIndexDoubleDict

struct Decisions{N}
    decisions::NTuple{N, DecisionMap}
    stage_objectives::NTuple{N, Vector{Tuple{MOI.OptimizationSense, MOI.AbstractScalarFunction}}}
    stage_map::StageMap
    constraint_map::ConstraintMap
    is_node::Bool

    function Decisions(::Val{N}; is_node::Bool = false) where N
        decisions = ntuple(Val(N)) do i
            DecisionMap()
        end
        stage_objectives = ntuple(Val(N)) do i
            Vector{Tuple{MOI.ObjectiveSense, MOI.AbstractScalarFunction}}()
        end
        return new{N}(decisions, stage_objectives, StageMap(), MOIU.DoubleDicts.IndexDoubleDict(), is_node)
    end

    function Decisions(decisions::NTuple{N, DecisionMap}; is_node::Bool = false) where N
        stage_objectives = ntuple(Val(N)) do i
            Vector{Tuple{MOI.ObjectiveSense, MOI.AbstractScalarFunction}}()
        end
        return new{N}(decisions, stage_objectives, StageMap(), MOIU.DoubleDicts.IndexDoubleDict(), is_node)
    end
end

function num_stages(::IgnoreDecisions)
    return 0
end
function num_stages(decisions::Decisions{N}) where N
    return N
end

function Base.getindex(decisions::Decisions, stage::Integer)
    return decisions.decisions[stage]
end

function has_decision(decisions::Decisions{N}, stage::Integer, index::MOI.VariableIndex) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    decisions.stage_map[index] == stage || error("Decision $index is mapped to stage $(decisions.stage_map[index])")
    return haskey(decisions.decisions[stage], index)
end
function has_decision(decisions::Decisions, index::MOI.VariableIndex)
    haskey(decisions.stage_map, index) || error("Decision $index not properly stage mapped.")
    return has_decision(decisions, decisions.stage_map[index], index)
end

function stage(decisions::Decisions, index::MOI.VariableIndex)
    haskey(decisions.stage_map, index) || error("Decision $index not properly stage mapped.")
    return decisions.stage_map[index]
end
function set_stage!(decisions::Decisions{N}, index::MOI.VariableIndex, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    decisions.stage_map[index] = stage
    return nothing
end

function decision(decisions::DecisionMap, index::MOI.VariableIndex)
    return decisions[index]
end
function decision(decisions::Decisions{N}, stage::Integer, index::MOI.VariableIndex) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return decision(decisions.decisions[stage], index)
end
function decision(decisions::Decisions, index::MOI.VariableIndex)
    haskey(decisions.stage_map, index) || error("Decision $index not properly stage mapped.")
    return decision(decisions, decisions.stage_map[index], index)
end

function decision_value(decisions::Decisions{N}, stage::Integer, index::MOI.VariableIndex) where N
    return decision_value(decision(decisions, stage, index))
end
function decision_value(decisions::Decisions, index::MOI.VariableIndex)
    return decision_value(decision(decisions, index))
end

function set_decision!(decisions::DecisionMap, index::MOI.VariableIndex, decision::Decision)
    decisions[index] = decision
    return nothing
end
function set_decision!(decisions::Decisions{N}, stage::Integer, index::MOI.VariableIndex, decision::Decision) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    set_decision!(decisions.decisions[stage], index, decision)
    return nothing
end
function set_decision!(decisions::Decisions{N}, index::MOI.VariableIndex, decision::Decision) where N
    haskey(decisions.stage_map, index) || error("Decision $index not properly stage mapped.")
    set_decision!(decisions, decisions.stage_map[index], index, decision)
    return nothing
end

function remove_decision!(::IgnoreDecisions, ::MOI.VariableIndex)
    return nothing
end
function remove_decision!(decisions::DecisionMap, index::MOI.VariableIndex)
    delete!(decisions, index)
    return nothing
end
function remove_decision!(decisions::Decisions{N}, stage::Integer, index::MOI.VariableIndex) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    if !has_decision(decisions, stage, index)
        return nothing
    end
    delete!(decisions.decisions[stage], index)
    delete!(decisions.stage_map, index)
    return nothing
end
function remove_decision!(decisions::Decisions{N}, index::MOI.VariableIndex) where N
    haskey(decisions.stage_map, index) || error("Decision $index not properly stage mapped.")
    remove_decision!(decisions, decisions.stage_map[index], index)
    return nothing
end

function get_stage_objective(decisions::Decisions{N}, stage::Integer, scenario_index::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = length(decisions.stage_objectives[stage])
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    return decisions.stage_objectives[stage][scenario_index]
end
function set_stage_objective!(decisions::Decisions{N},
                              stage::Integer,
                              scenario_index::Integer,
                              sense::MOI.OptimizationSense,
                              objective::MOI.AbstractScalarFunction) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    n = length(decisions.stage_objectives[stage])
    1 <= scenario_index <= n || error("Scenario index $scenario_index not in range 1 to $n.")
    decisions.stage_objectives[stage][scenario_index] = (sense, objective)
    return nothing
end

function add_stage_objective!(decisions::Decisions{N},
                              stage::Integer,
                              sense::MOI.OptimizationSense,
                              objective::MOI.AbstractScalarFunction) where N
    if decisions.is_node
        # No need to cache objective if model is a node problem
        return nothing
    end
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    push!(decisions.stage_objectives[stage], (sense, objective))
    return nothing
end

function mapped_constraint(decisions::Decisions, ci::MOI.ConstraintIndex)
    if haskey(decisions.constraint_map, ci)
        return decisions.constraint_map[ci]
    else
        return typeof(ci)(0)
    end
end
function map_constraint!(decisions::Decisions, ci::MOI.ConstraintIndex, inner::MOI.ConstraintIndex)
    decisions.constraint_map[ci] = inner
    return nothing
end
function remove_mapped_constraint!(decisions::Decisions, ci::MOI.ConstraintIndex)
    haskey(decisions.constraint_map, ci) || error("Constraint $ci not properly mapped.")
    delete!(decisions.constraint_map, ci)
    return nothing
end

function clear!(decisions::Decisions)
    map(decisions.decisions) do decisions
        empty!(decisions)
    end
    map(empty!, decisions.stage_objectives)
    empty!(decisions.stage_map)
    empty!(decisions.constraint_map)
    return nothing
end

function all_decisions(decisions::DecisionMap)
    return filter(decisions) do (index, decision)
        state(decision) != Known
    end |> keys |> collect
end
function all_decisions(decisions::Decisions{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return all_decisions(decisions[stage])
end

function all_known_decisions(decisions::DecisionMap)
    return filter(decisions) do (index, decision)
        state(decision) == Known
    end |> keys |> collect
end
function all_known_decisions(decisions::Decisions{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return all_known_decisions(decisions[stage])
end

function num_decisions(decisions::DecisionMap)
    return count(decisions) do (index, decision)
        return state(decision) != Known
    end
end
function num_decisions(decisions::Decisions{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return num_decisions(decisions[stage])
end

function num_known_decisions(decisions::DecisionMap)
    return count(decisions) do (index, decision)
        state(decision) == Known
    end
end
function num_known_decisions(decisions::Decisions{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return num_known_decisions(decisions[stage])
end

function untake_decisions!(decisions::DecisionMap)
    need_update = false
    # Update states and values
    for key in all_decisions(decisions)
        if decisions[key].state != NotTaken
            decisions[key].state = NotTaken
            need_update |= true
        end
    end
    return need_update
end
function untake_decisions!(decisions::Decisions{N}, stage::Integer) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    return untake_decisions!(decisions[stage])
end

function take_decisions!(decisions::DecisionMap, x::AbstractVector)
    # Check decision length
    num_decisions(decisions) == length(x) || error("Given decision of length $(length(x)) not compatible with number of defined decision variables $(num_decisions(decisions)).")
    # Update states and values (assume x given in sorted order)
    for (key, val) in zip(all_decisions(decisions), x)
        decisions[key].state = Taken
        decisions[key].value = val
    end
    return nothing
end
function take_decisions!(decisions::Decisions{N}, stage::Integer, x::AbstractVector) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    take_decisions!(decisions[stage], x)
    return nothing
end

function update_known_decisions!(decisions::DecisionMap, x::AbstractVector)
    # Check decision length
    num_known_decisions(decisions) == length(x) || error("Given decision of length $(length(x)) not compatible with number of defined known decision variables $(num_known_decisions(decisions)).")
    # Update values (assume x given in sorted order)
    for (key, val) in zip(all_known_decisions(decisions), x)
        decisions[key].value = val
    end
    return nothing
end
function update_known_decisions!(decisions::Decisions{N}, stage::Integer, x::AbstractVector) where N
    1 <= stage <= N || error("Stage $stage not in range 1 to $N.")
    update_known_decisions!(decisions[stage], x)
end

is_decision_type(::DataType) = false

include("variable_interface.jl")
include("expressions/expressions.jl")
include("functions/functions.jl")
include("sets.jl")
include("modifications.jl")
include("bridges/bridges.jl")
include("macros.jl")
include("updates.jl")
include("moi_overrides.jl")


text:
"""
    PerformanceOpt

Abstract supertype for performance optimisation flags.
"""
abstract type PerformanceOpt end

"""
    SparseOpt <: PerformanceOpt

    SparseOpt()

An optimisation flag that ignores all zero values in the grid.

For low-density simulations performance may improve by
orders of magnitude, as only used cells are run.

This is complicated for optimising neighborhoods - they
must run if they contain just one non-zero cell.

Specifiy with:

```julia
ruleset = Ruleset(rule; opt=SparseOpt())
# or
output = sim!(output, rule; opt=SparseOpt())
```

`SparseOpt` is best demonstrated with this simulation, where the grey areas do not
run except where the neighborhood partially hangs over an area that is not grey:

![SparseOpt demonstration](https://raw.githubusercontent.com/cesaraustralia/DynamicGrids.jl/media/complexlife_spareseopt.gif)
"""
struct SparseOpt{F<:Function} <: PerformanceOpt
    f::F
end
SparseOpt() = SparseOpt(==(0))

@inline _isactive(val, opt::SparseOpt{<:Function}) = !opt.f(val)

"""
    NoOpt <: PerformanceOpt

    NoOpt()

Flag to run a simulation without performance optimisations besides basic high performance
programming. Still fast, but not intelligent about the work that it does: all cells are run
for all rules.

`NoOpt` is the default `opt` method.
"""
struct NoOpt <: PerformanceOpt end

"""
    Processor

Abstract supertype for selecting a hardware processor, such as ia CPU or GPU.
"""
abstract type Processor end

"""
    CPU <: Processor

Abstract supertype for CPU processors.
"""
abstract type CPU <: Processor end

"""
    SingleCPU <: CPU

    SingleCPU()

[`Processor`](@ref) flag that specifies to use a single thread on a single CPU.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=SingleCPU())
# or
output = sim!(output, rule; proc=SingleCPU())
```
"""
struct SingleCPU <: CPU end

"""
    ThreadedCPU <: CPU

    ThreadedCPU()

[`Processor`](@ref) flag that specifies to use a `Threads.nthreads()` CPUs.

Specifiy with:

```julia
ruleset = Ruleset(rule; proc=ThreadedCPU())
# or
output = sim!(output, rule; proc=ThreadedCPU())
```
"""
struct ThreadedCPU{L} <: CPU
    spinlock::L
end
ThreadedCPU() = ThreadedCPU(Base.Threads.SpinLock())
Base.Threads.lock(opt::ThreadedCPU) = lock(opt.spinlock)
Base.Threads.unlock(opt::ThreadedCPU) = unlock(opt.spinlock)

"""
    BoundaryCondition

Abstract supertype for flags that specify the boundary conditions used in the simulation,
used in [`inbounds`](@ref) and to update [`NeighborhoodRule`](@ref) grid padding.
These determine what happens when a neighborhood or jump extends outside of the grid.
"""
abstract type BoundaryCondition end

"""
    Wrap <: BoundaryCondition

    Wrap()

[`BoundaryCondition`](@ref) flag to wrap cordinates that boundary boundaries back to the
opposite side of the grid.

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Wrap())
# or
output = sim!(output, rule; boundary=Wrap())
```
"""
struct Wrap <: BoundaryCondition end

"""
    Remove <: BoundaryCondition

    Remove()

[`BoundaryCondition`](@ref) flag that specifies to assign `padval` to cells that overflow
grid boundaries. `padval` defaults to `zero(eltype(grid))` but can be assigned as a keyword
argument to an [`Output`](@ref).

Specifiy with:

```julia
ruleset = Ruleset(rule; boundary=Remove())
# or
output = sim!(output, rule; boundary=Remove())
```
"""
struct Remove <: BoundaryCondition end


text:
####### SIMULATION OF MARGINAL ###########################################

### +++ INITIALIZATION +++

### import packages
using Random, LinearAlgebra, RDatasets, MixedModels, Distributions, Optim, CSV, Dates, Statistics, CategoricalArrays

### define functions
vcov2 = function(x, xbar, σ, τ, n, γ)
    m = length(xbar)
    a = γ / τ^2 * m # 1-1 entry
    d = (sum(x.^2) - n * γ * sum(xbar.^2) ) / σ^2 # 2-2 entry
    b = γ / τ^2 / n * sum(x) # 2-1 entry
    [d -b; -b a] / (a * d - b^2) # return inverse (X^t V^-1 X)^-1
end

beta = function(vcov1, x, xbar, y, ybar, σ, τ, n, γ)
    c1 = γ / τ^2 * sum(ybar) # = 1^t V^-1 y
    c2 = (x' * y - n * γ * sum(ybar .* xbar))/σ^2 # = x^t V^-1 y
    vcov1 * [c1 c2]'  # (X^tV^-1X)^-1 X^tV^-1y
end

Q = function(data, y, ybar, x, xbar, n, N, μ0, quantile)

    ### fit model
    model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)

    ### make sure the variance parameters are positive
    τ = maximum([VarCorr(model).σρ.v.σ[1], 0.001])
    σ = maximum([VarCorr(model).s, 0.001])

    ###     MARGINAL

    ### helpful parameters
    γ = τ^2 / (τ^2 + σ^2 / n)
    vcov1 = vcov2(x, xbar, σ, τ, n, γ) # = (X^tV^-1X)^-1
    d = (1-γ) .* hcat(fill(1, m), xbar) # d' = l' - b'X
    dd = d' * d
    a = n * τ^2 + σ^2

    ### assemble diagonal matrix parts
    g1 = γ * σ^2 / n
    g3 = 2 * σ^4 / τ^2 * γ / N/ (n - 1)

    ### assemble inverse of covariance matrix with Woodbury
    Σdiag = g1 + 2 * g3
    #Σ = Σdiag * I + d * vcov1 * d'
    Σi = (I - d * inv(Σdiag * inv(vcov1) + dd) * d' ) / Σdiag

    ###     CONDITIONAL

    ### assemble diagonal matrix parts
    l1 = γ^2 / n * σ^2

    l3_1 = 2 * σ^2 * τ^2 / N/ a^3
    l3_2 = 8 / N/ (n - 1) / a^4 / n^2 * τ^2 * σ^2 * (
        + (a^2 - σ^2) / n * (
            + σ^6 / a^2 * n
            - (σ^4  / a^2) * (τ + a / 2 * n)
            - (τ^2 - a / 2 + σ^2 * (1 + n^2)) * (n - 1))
        + σ^10 / a^2 * (n^3 - (n - 1)^2 + 1 / n)
        + a * σ^4 * (n - 1 + 1 / n / 2)  +
        - σ^6 * (2 * n - 1 + 1 / n) +
        + σ^8 / a * (1 + (n-1) * (n^2 + 1) / n) / 2)
    l3_3 = - 8 / N/ n * τ^2 * σ^4 / a^4 * (n - 1) * (a^2 - σ^4)
    l3_4 = - 8 / N/ n * τ^2 * σ^4 / a^3 / (n - 1) * (
        + a * ((n - 1) * n - 1)
        + σ^2 * (n + 1))
    l3 = 2 * (l3_1 + l3_2 + l3_3 + l3_4)

    l4 = 2 * σ^4 / N / a^2 * ( 1 / n + ((n - 1) * τ^2 +
            (σ^2 + τ^2)^2) / ((n - 1) * a^2))

    l5 = 2 * n / (n - 1) / N * σ^4 * (σ^2 - 2 * n * τ^2) / a^2

    ### assemble inverse of covariance matrix with Woodbury
    cΣdiag = maximum([l1 + l3 + l4 - l5, 0.001])
    cΣi = (I - d * inv(cΣdiag * σ^2 *
                inv(vcov1 * (2 * γ * I + n^2 * dd * vcov1)) + dd) * d' ) / cΣdiag

    ###     COMPUTE MEAN
    β = beta(vcov1, x, xbar, y, ybar, σ, τ, n, γ) # = fixef(model)
    μ = (1 - γ) * (β[1] .+ β[2] .* xbar) + γ * ybar

    ###     SAVE PIVOTS
    MarginalPivot = sum((μ - μ0)' .* ((μ - μ0)' * Σi), dims = 2)[1]
    ConditionalPivot = sum((μ - μ0)' .* ((μ - μ0)' * cΣi), dims = 2)[1]

    if quantile == true

        ###     NONCENTRALITY PARAMETER
        ay = d * β - (1 - γ) * ybar
        λ = maximum([(ay' * cΣi * ay)[1] - tr(cΣi) * σ^2 / a^2 / n, 0])

        ###     QUANTILES
        qc = cquantile(NoncentralChisq(m, λ), 0.05) # quantile of alpha = 0.05

        return [MarginalPivot ConditionalPivot qc β[1] β[2] τ σ Σdiag cΣdiag]
    else
        return [MarginalPivot ConditionalPivot]
    end
end

### fix parameters
p = 2
nsim = 5000
R = 300
Random.seed!(2021)

### +++ PER SPECIFICATION SETTING +++

### to compile
τ0 = 8^0.5
σ0 = 4^0.5
m = 50
n = 10
REdistr = "normal"
Edistr = "normal"

### variable parameters
N = m * n

β0 = [-4.9, 0.03]
x = 10 * randn(N)
xbar = vec(sum(reshape(x, n, :), dims=1) / n)

qm = cquantile(Chisq(m), 0.05) # quantile of alpha = 0.05

### conditional evaluation
if REdistr == "normal"
    v = τ0 * randn(m)
elseif REdistr == "tdist"
    v = rand(TDist(2), m)
elseif REdistr == "chi2"
    v = rand(Chisq(τ0^2 / 2), m) .- τ0^2 / 2
end
u = β0[1] .+ β0[2] * x + repeat(v, inner = n)
μ0 = β0[1] .+ β0[2] .* xbar + v

### obtain conditional variance
σ1 = zeros(10000, 1)
τ1 = zeros(10000, 1)
for i in 1:10000
    if Edistr == "normal"
        e = σ0 * randn(N)
    elseif Edistr == "tdist"
        e = rand(TDist(2), N)
    elseif Edistr == "chi2"
        e = rand(Chisq(σ0^2 / 2), N) .- σ0^2 / 2
    end
    y = u + e # conditional
    data = DataFrame(x = x, v = repeat(categorical(v), inner = n), y = y)
    model = fit(MixedModel, @formula(y ~ x + (1 | v)), data)
    τ1[i] = VarCorr(model).σρ.v.σ[1]
    σ1[i] = VarCorr(model).s
end
σ1 = mean(σ1) # should be (close to) σ0, but theory for conditional convergence is missing
τ1 = mean(τ1) # should be (close to) τ0, but theory for conditional convergence is missing

### simulation
Coverage = zeros(6)
λ = zeros(nsim)
Volume = ones(nsim, 6)
QM1 = zeros(nsim)
QC1 = zeros(nsim)
QM2 = zeros(nsim)
QC2 = zeros(nsim)
for i in 1:nsim

    ### +++ EMPIRICAL +++

    ### variable data
    if Edistr == "normal"
        e = σ0 * randn(N)
    elseif Edistr == "tdist"
        e = rand(TDist(2), N)
    elseif Edistr == "chi2"
        e = rand(Chisq(σ0^2 / 2), N) .- σ0^2 / 2
    end
    y = u + e # conditional
    ybar = vec(sum(reshape(y, n, :), dims=1) / n)
    data = DataFrame(x = x, v = repeat(categorical(v), inner = n), y = y)

    MP_bt_Alg1 = zeros(R-1)
    CP_bt_Alg1 = zeros(R-1)

    β1 = 1;  β2 = 1; τ = 1;  σ = 1

    MP, CP, qc, β1, β2, τ, σ, Σdiag, cΣdiag= Q(data, y, ybar, x, xbar, n, N, μ0, true)

    for r in 2:R

        v_bt = τ * randn(m)
        u_bt = β1 .+ β2 * x + repeat(v_bt, inner = n)
        μ0_bt = β1 .+ β2 .* xbar + v_bt

            ### generate new bootstrap data
            y_bt = u_bt + σ * randn(N) # conditional
            ybar_bt = vec(sum(reshape(y_bt, n, :), dims = 1) / n)
            data_bt = DataFrame(x = x, v = repeat(categorical(v_bt), inner = n), y = y_bt)

            MP_bt_Alg1[r-1], CP_bt_Alg1[r-1] = Q(data_bt, y_bt, ybar_bt, x, xbar, n, N, μ0_bt, false)
    end

    MP_bt_Alg2 = zeros(R-1)
    CP_bt_Alg2 = zeros(R-1)

    v_bt = τ * randn(m)
    u_bt = β1 .+ β2 * x + repeat(v_bt, inner = n)
    μ0_bt = β1 .+ β2 .* xbar + v_bt

    for r in 2:R

            ### generate new bootstrap data
            y_bt = u_bt + σ * randn(N) # conditional
            ybar_bt = vec(sum(reshape(y_bt, n, :), dims = 1) / n)
            data_bt = DataFrame(x = x, v = repeat(categorical(v_bt), inner = n), y = y_bt)

            MP_bt_Alg2[r-1], CP_bt_Alg2[r-1] = Q(data_bt, y_bt, ybar_bt, x, xbar, n, N, μ0_bt, false)
    end

    ###     COMPUTE BOOTSTRAP QUANTILES
    QM1[i] = quantile!(MP_bt_Alg1, 0.95)
    QC1[i] = quantile!(CP_bt_Alg1, 0.95)

    ###     COMPUTE BOOTSTRAP QUANTILES
    QM2[i] = quantile!(MP_bt_Alg2, 0.95)
    QC2[i] = quantile!(CP_bt_Alg2, 0.95)

    ###     COMPUTE COVERAGE
    Coverage[1] += (MP .< qm) / nsim
    Coverage[2] += (CP .< qc) / nsim
    Coverage[3] += (MP .< QM1[i]) / nsim
    Coverage[4] += (CP .< QC1[i]) / nsim
    Coverage[5] += (MP .< QM2[i]) / nsim
    Coverage[6] += (CP .< QC2[i]) / nsim

    ###     COMPUTE RELATIVE VOLUME
    # (via approximations as the determinant of large matrices is numerically unstable)
    Volume[i, 2] = (Σdiag / cΣdiag * qm / qc)^(-m / 2)
    Volume[i, 3] = (qm / QM1[i])^(-m / 2)
    Volume[i, 4] = (Σdiag / cΣdiag * qc / QC1[i])^(-m / 2)
    Volume[i, 5] = (qm / QM2[i])^(-m / 2)
    Volume[i, 6] = (Σdiag / cΣdiag * qc / QC2[i])^(-m / 2)

end


text:
## target

export Target,
       name, description,
       hasjit, hastargetmachine, hasasmparser

@checked struct Target
    ref::API.LLVMTargetRef
end

Base.unsafe_convert(::Type{API.LLVMTargetRef}, target::Target) = target.ref

function Target(; name=nothing, triple=nothing)
    (name !== nothing) ⊻ (triple !== nothing) ||
        throw(ArgumentError("Specify either name or triple."))

    if triple !== nothing
        target_ref = Ref{API.LLVMTargetRef}(0)
        error_ref = Ref{Cstring}(C_NULL)
        status = convert(Core.Bool, API.LLVMGetTargetFromTriple(triple, target_ref, error_ref))
        if status && error_ref[] !== C_NULL
            error = unsafe_message(error_ref[])
            throw(ArgumentError(error))
        elseif status
            throw(ArgumentError("Cannot find a target for triple '$triple'"))
        end
        @assert target_ref[] != C_NULL
        return Target(target_ref[])
    elseif name !== nothing
        target = API.LLVMGetTargetFromName(name)
        if target == C_NULL
            throw(ArgumentError("Cannot find target '$triple'"))
        end
        return Target(target)
    end
end

name(t::Target) = unsafe_string(API.LLVMGetTargetName(t))

description(t::Target) = unsafe_string(API.LLVMGetTargetDescription(t))

hasjit(t::Target) = convert(Core.Bool, API.LLVMTargetHasJIT(t))
hastargetmachine(t::Target) = convert(Core.Bool, API.LLVMTargetHasTargetMachine(t))
hasasmparser(t::Target) = convert(Core.Bool, API.LLVMTargetHasAsmBackend(t))

function Base.show(io::IO, ::MIME"text/plain", target::Target)
  print(io, "LLVM.Target($(name(target))): $(description(target))")
end

## target iteration

export targets

struct TargetIterator end

targets() = TargetIterator()

Base.eltype(::TargetIterator) = Target

function Base.iterate(iter::TargetIterator, state=API.LLVMGetFirstTarget())
    state == C_NULL ? nothing : (Target(state), API.LLVMGetNextTarget(state))
end

Base.IteratorSize(::TargetIterator) = Base.SizeUnknown()


text:
## ExoplanetsSysSim/src/setup.jl
## (c) 2015 Eric B. Ford

# Since bitbucket messed up capitalization of package name
if ! isdir( joinpath(Pkg.dir(),"ExoplanetsSysSim") )
     symlink( joinpath(Pkg.dir(),"exoplanetssyssim"), joinpath(Pkg.dir(),"ExoplanetsSysSim") )
end

startdir = pwd()
# Code to be run once to install unregistered package dependencies
# Code to be run once to install non-standard dependencies
Pkg.clone("git@github.com:eford/ABC.jl.git")
# switch to dragozzine's fork for any minor SysSim specific ABC changes

Pkg.clone("git@github.com:jbrakensiek/CORBITS.git")

# Compile CORBITS library and put it somewhere we can find
cd(joinpath(Pkg.dir(),"CORBITS"))
run(`make lib`)
#cd(homedir())
cd(workdir)

if !is_windows()
   symlink( joinpath(Pkg.dir("CORBITS"),"libcorbits.so"), joinpath(Pkg.dir("ExoplanetsSysSim"),"libcorbits.so") )
else
   cp( joinpath(Pkg.dir("CORBITS"),"libcorbits.so"), joinpath(Pkg.dir("ExoplanetsSysSim"),"libcorbits.so") )
end


text:
#!/usr/bin/env julia

import Celeste.DeterministicVI
import Celeste.DeterministicVI: ElboIntermediateVariables, elbo_likelihood

include(joinpath(Pkg.dir("Celeste"), "test", "SampleData.jl"))

import SampleData

function benchmark_elbo_likelihood()
    ea, vp, catalog = SampleData.gen_sample_star_dataset()

    # Warm up---this compiles the code
    DeterministicVI.elbo_likelihood(ea, vp)

    # clear allocations in case julia is running with --track-allocations=user
    Profile.clear_malloc_data()

    if isempty(ARGS)
        @time DeterministicVI.elbo_likelihood(ea, vp)
    elseif ARGS[1] == "--profile"
        Profile.init(delay=0.0001)
        @profile DeterministicVI.elbo_likelihood(ea, vp)
        Profile.print(format=:flat, sortedby=:count)
    end
end

benchmark_elbo_likelihood()


text:
require("ArgParse")
require("Iterators")

using ArgParse, Iterators
import JSON
push!(LOAD_PATH, dirname(@__FILE__()))

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--E-ONE", "-o"
            help = "converter for e-one energy term"
            action = :store_true
            #required = true
        "--input-file"
            help = "file with E-one values"
            required = true
        "--output-file"
            help = ".json file name form saving results"
            required = true
    end

    return parse_args(s)
end
aaNames = ["GLY", "ALA", "SER", "CYS",
  "VAL", "THR", "ILE", "PRO",
  "MET", "ASP", "ASN", "LEU",
  "LYS", "GLU", "GLN", "ARG",
  "HIS", "PHE", "TYR", "TRP",
  "CYX"]
aaIndices = Dict(collect(zip(aaNames, 1:21)))

immutable E12Data
  aa1 :: String
  aa2 :: String
  n14 :: Int
  e214 :: Array{Number, 1}
end

function process_r14aa12(result :: Dict{String, Dict{String, E12Data}})
  {"data"=>
  [
  aa1 =>
    [
      aa2 => result[aa1][aa2]
      for aa2 in intersect(aaNames, keys(result[aa1]))
    ]
    for aa1 in intersect(aaNames, keys(result))
    ],
    "startIndex" => -86,
    "endIndex" => 91,
    "binRightBorders" => [-79, -55, -35, -25, -12, 14, 26, 46, 65, 72, 81]
    }
end

function readE12part(input_file)
  (aa1, aa2, n14) = split(rstrip(readline(input_file)))
  numbers = Number[]
  while(!eof(input_file))
    s = readline(input_file)
    if beginswith(s, "=")
      break
    end
    append!(numbers, map(x->float(x), split(rstrip(s))))
  end
  return (aa1, aa2, int(n14), numbers)
end

function main()
    parsed_args = parse_commandline()
    input_file = open(parsed_args["input-file"], "r")
    result = Dict{String, Dict{String, E12Data}}()
    while !eof(input_file)
      (aa1, aa2, n14, numbers) = readE12part(input_file)
      if !haskey(result, aa1)
        result[aa1] = Dict{String, E12Data}()
      end
      result[aa1][aa2] = E12Data(aa1, aa2, n14, numbers)
      #push!(result, rstrip(readline(input_file)))
    end
    close(input_file)
    r1 = process_r14aa12(result)
    output_file = open(parsed_args["output-file"], "w")
    println(output_file, JSON.json(r1, 1))
    close(output_file)
end

main()


text:
using PATHSolver

# define functions to calculate residual
function index(m, n, M,N)
    return (n-1)*M + m
end

function alln(m, M,N)
    return [index(m, n, M,N) for n = 1:N]
end

function allm(n, M,N)
    return [index(m, n, M,N) for m = 1:M]
end

function f(c,h,γ)
    return (c.^γ + h.^γ).^(1./γ)
end

function clearing(q, v, g, M, N)
    ret = zeros(eltype(q),N)
    for n = 1:N
        ret[n] = g[n] - sum(q[allm(n, M, N)].*v)
    end
    return ret
end

function wages(q, ch, ρ, γ, r, pc, M, N)
    w = zeros(eltype(q),M*N)
    Q = Qvec(q, M, N)
    yj = matchoutput(ch[:,2],ch[:,1],ρ, γ,r)
    kj = zeros(M*N)

     for m = 1:M
        for n = 1:N
            i = index(m, n, M,N)
            ip1 = index(m, n+1, M,N)
            if q[i] > 0.
                w[i] = pc[m]*yj[i] - r*kj[i]
                for nprime in 1:n-1
                    iprime = index(m, nprime, M,N)
                    iprimep1 = index(m, nprime+1, M,N)
                    w[i] -= exp(-(Q[iprimep1]-Q[i]))* (1-exp(-q[iprime])) * max(0., pc[m]*yj[iprime] - r*kj[iprime])
                end
                w[i] *= q[i]*exp(-q[i])/(1-exp(-q[i]))
            end
        end
    end
    return w
end

function expected_u(q, w, M, N)
    Q = Qvec(q, M, N)
    u = zeros(eltype(q),M*N)
    for m = 1:M
        for n = 1:N
            i = index(m, n, M,N)
            if n < N
                ip1 = index(m, n+1, M,N)
                u[i] = exp(-Q[ip1])  * w[i]
            else
                u[i] =  w[i]
            end
            if q[i]>1e-9
                u[i] = u[i]*(1-exp(-q[i]))/q[i]
            end
        end
    end
    return u
end

function allocation(q, ch, ρ, γ, v, g, r, σ, M, N)
    e = employment(q, v, M, N)
    yj = e.*matchoutput(ch[:,2],ch[:,1],ρ, γ,r)
    yc = [sum(yj[alln(m, M, N)]) for m = 1:M]
    yf = sum(yc.^((σ-1)/σ))^((σ/(σ-1)))
    pc = (yc./yf).^(-1/σ)

    w = wages(q, ch, ρ, γ, r, pc, M, N)
    u = expected_u(q, w, M, N)
    return w, u, pc, yc, yf, e
end

function matchoutput(h,c,ρ, γ,r)
    fj = f(c,h,γ)
    k = 0.
    return (fj.^ρ + k.^ρ).^(1./ρ)
end

function Qvec(q, M, N)
    Q = zeros(eltype(q),M*N)
    for m = 1:M
        for n=N:-1:1
            i = index(m, n, M,N)
            if n == N
                Q[i] = q[i]
            else
                ip1 = index(m, n+1, M,N)
                Q[i] = Q[ip1] + q[i]
            end
        end
    end
    return Q
end

function employment(q, v, M, N)
    # employment
    Q = Qvec(q, M, N)
    e = zeros(eltype(q),M*N)
    qi = 0.
    Qp1 = 0.
    i = 0
    for m = 1:M
        for n=1:N
            i = index(m, n, M, N)
            if n < N
                Qp1 = Q[index(m, n+1, M, N)]
            else
                Qp1 = 0.
            end
            if q[i] > 0.
                e[i] = v[m]*exp(-Qp1) * (1-exp(-q[i]))
            end
        end
    end
    return e
end

function resid(qu, N, M, ch, ρ, γ, v, g, r, σ)
    res = zeros(eltype(qu), size(qu))
    q = qu[1:end-N]
    umax = qu[end-N+1:end]
    w, u, pc, yc, yf, e = allocation(q, ch, ρ, γ, v, g, r, σ, M, N)
    res[1:end-N] = -[u[index(m,n, M, N)] - umax[n] for n=1:N, m=1:M][:]
    res[end-N+1:end] =  clearing(q, v, g, M, N)
    return res
end

# Set parameters
# workers
N = 5
# jobs
M = 5
c = linspace(1.,2.,M)
if N > 1
    h = linspace(0.0,2.,N)
else
    h = [1.]
end

ch = hcat(repeat(c, outer=N), repeat(h, inner=M))
γ = -.5
r = .01
σ = 2.

qguess = ones(M*N)
uguess = 1e-6*ones(N)
quguess = vcat(qguess, uguess)
ρ = ones(qguess)*0.5

v = ones(M)
g = ones(N)

# wrapper to objective function obj(x) to be set to zero
obj(x) = resid(x, N, M, ch, ρ, γ, v, g, r, σ)

# set bounds for x
lower = vcat(zeros(M*N), zeros(N))
upper = 10000.*ones(M*N+N)

status, zero, f_zero = solveMCP(obj, lower, upper, quguess)


text:
const JLDMODE=Val(0)
const GPUMODE=Val(1)
const CPUMODE=Val(2)

serialize(x) = _ser(x,IdDict(),JLDMODE)
gpucopy(x)   = _ser(x,IdDict(),GPUMODE)
cpucopy(x)   = _ser(x,IdDict(),CPUMODE)

function _ser(x::KnetPtr,s::IdDict,::typeof(JLDMODE))
    if !haskey(s,x)
        if isa(x.ptr, Cptr) && (x.dev >= 0)
            a = Array{UInt8}(undef,x.len)
            @cudart(cudaMemcpy,(Cptr,Cptr,Csize_t,UInt32),pointer(a),x.ptr,x.len,2)
            s[x] = KnetPtr(a,x.len)
        elseif isa(x.ptr, Array{UInt8,1})
            if gpu() >= 0
                g = knetMalloc(x.len); if g === nothing; error("No space left on GPU."); end
                @cudart(cudaMemcpy,(Cptr,Cptr,Csize_t,UInt32),g,pointer(x.ptr),x.len,1)
                s[x] = KnetPtr(g,x.len,gpu())
            else
                s[x] = x  # Leave conversion to array to KnetArray
            end
        else
            error("Unrecognized KnetPtr")
        end
    end
    return s[x]
end

function _ser(x::KnetArray{T,N},s::IdDict,m::typeof(JLDMODE)) where {T,N}
    if !haskey(s,x)
        if isa(x.ptr.ptr, Array) && gpu() < 0
            s[x] = copy(reshape(reinterpret(eltype(x),view(x.ptr.ptr,1:sizeof(T)*length(x))),size(x)))
        else
            s[x] = KnetArray{T,N}(_ser(x.ptr,s,m),x.dims)
        end
    end
    return s[x]
end

function _ser(x::RNN, s::IdDict, m::typeof(JLDMODE))
    if !haskey(s,x)
        dropoutDesc = (x.dropoutDesc != nothing || gpu() < 0 ? nothing : DD(handle=gethandle(),dropout=x.dropout,seed=x.seed))
        rnnDesc = (x.rnnDesc != nothing || gpu() < 0 ? nothing : RD(x.hiddenSize,x.numLayers,dropoutDesc,x.inputMode,x.direction,x.mode,x.algo,x.dataType))
        s[x] = RNN(_ser(x.w,s,m), _ser(x.h,s,m), _ser(x.c,s,m), x.inputSize, x.hiddenSize, x.numLayers, x.dropout, x.seed, x.inputMode, x.direction, x.mode, x.algo, x.dataType, rnnDesc, dropoutDesc, _ser(x.dx,s,m), _ser(x.dhx,s,m), _ser(x.dcx,s,m))
    end
    return s[x]
end
# Partially fixes the issue: when KA converts to A because no gpu, surrounding parametric types remain Param{KA}.
# However other container types that include KnetArray may still have an inconsistent parametric type problem.
_ser(x::Param, s::IdDict, m::Val)=(haskey(s,x) ? s[x] : s[x]=Param(_ser(x.value,s,m),_ser(x.opt,s,m)))

_ser(x::KnetArray,s::IdDict,::typeof(GPUMODE))=x
_ser(x::KnetArray,s::IdDict,::typeof(CPUMODE))=(haskey(s,x) ? s[x] : s[x]=Array(x))
_ser(x::Array, s::IdDict, m::Val) = (haskey(s, x) ? s[x] : s[x] = _ser_array_t(x, eltype(x), s, m))

function _ser_array_t(@nospecialize(x), T, s::IdDict, m::Val)
    if !isbitstype(T)
        map(xi->_ser(xi,s,m), x)
    elseif m === GPUMODE
        KnetArray(x)
    else
        x
    end
end

# Generic serialization rules from deepcopy.jl
_ser(x::Union{Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,UnionAll,Task},::IdDict,::Val) = x
_ser(x::Tuple, s::IdDict, m::Val) = ntuple(i->_ser(x[i], s, m), length(x))
_ser(x::Module, ::IdDict, ::Val) = error("serialize of Modules not supported")
_ser(x::Core.SimpleVector, s::IdDict,m::Val) = (haskey(s, x) ? s[x] : s[x] = Core.svec(Any[_ser(x[i], s, m) for i = 1:length(x)]...))
_ser(x::String, s::IdDict,::Val) = (haskey(s, x) ? s[x] : s[x] = (GC.@preserve x unsafe_string(pointer(x), sizeof(x))))

function _ser(@nospecialize(x), s::IdDict, m::Val)
    T = typeof(x)::DataType
    nf = nfields(x)
    (isbitstype(T) || nf == 0) && return x
    if haskey(s, x)
        return s[x]
    end
    y = ccall(:jl_new_struct_uninit, Any, (Any,), T)
    if T.mutable
        s[x] = y
    end
    for i in 1:nf
        if isdefined(x,i)
            ccall(:jl_set_nth_field, Cvoid, (Any, Csize_t, Any), y, i-1,
                  _ser(getfield(x,i), s, m))
        end
    end
    return y::T
end

function _ser(x::Union{Dict,IdDict}, s::IdDict,m::Val)
    if haskey(s, x)
        return s[x]
    end
    if isbitstype(eltype(x))
        return (s[x] = x)
    end
    dest = empty(x)
    s[x] = dest
    for (k, v) in x
        dest[_ser(k, s, m)] = _ser(v, s, m)
    end
    dest
end


text:
## Julia 1.1.0

# Initialize
using Plots
plotlyjs() # Backend PlotlyJS

# Archimedean Spiral
# r = a + bΘ

a = 1.0
b = 2.0
angles = 8π
Θ = range(0, stop = angles, length = 1000)
r = a .+ b .* Θ

plot(Θ, r, proj=:polar)
savefig("./png/archimedean-spiral.png")

steps = 100
anim = @animate for i = 1:steps
    Θ = range(0, stop = angles * i / steps, length = 10 * i)
    r = a .+ b .* Θ
    plot(Θ, r, proj=:polar)
end

gif(anim, "./gif/archimedean-spiral.gif", fps = 30)

max_stars_repo_name: fengyfei/sapphire

max_issues_repo_name: fengyfei/sapphire

max_forks_repo_name: fengyfei/sapphire


text:
module Util

# import MacroTools: postwalk, prewalk
using MLStyle
using ..AExpressions

export postwalkstate,
       postwalkpos

"""
Wraps MacroToolls.postwalk:

postwalk(f, expr)
Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *after* they have been transformed by the walk. See also
`prewalk`.
"""

walk(x, inner, outer, mapf) = outer(x)
walk(x::AExpr, inner, outer, mapf = map) =
  outer(AExpr(x.head, mapf(inner, x.args)...))

postwalk(f, x, mapf = map) = walk(x, x -> postwalk(f, x, mapf), f, mapf)

@inline mapenumerate(f, xs) = map(f, enumerate(xs))
@inline mapid(f, xs) = map(f, 1:length(xs), xs)

# ## With state

"""
Postwalk with state.

`postwalkstate(x::Expr, f, state, statef, mapf = map)`

postwalk(f, expr)
Applies `f(expr, state)` to each node in the given expression tree, returning the result.
`f` sees expressions *after* they have been transformed by the walk.
`state` seen by each node is inductively defined:
  `state` for `x` is `state`
  `state` for each child of is `stataf(state, parent, i)` where
    `parent` is the parent of
    `i` is position of child in `parent`
```
s0 = Int[]
statef(state, arg, i) = @show [i; state]
expr = AExpr(:(let x = 4, y = 3
  x + y
end))
f(expr, state) = expr
postwalkstate(expr, f, s0, statef)
```
"""
postwalkstate(f, x::AExpr, state, statef) =
  let g(i, x_) = postwalkstate(f, x_,  statef(state, x, i), statef)
    # @show x.args
    f(AExpr(x.head, mapid(g, x.args)...), state)
  end

postwalkstate(f, x, state, statef) = f(x, state)

idappend(state, arg, i) = [state; i]
postwalkpos(f, x, p0 = Int[]) = postwalkstate(f, x, p0, idappend)
"""

"""
function statewalk()
end

end


text:
# degree matrix
degree_matrix(x::T) where {T<:AbstractHyperGraph} = Diagonal(degrees(x))

# incidence matrix (oriented case)
@traitfn function incidence_matrix(x::T::IsOriented) where {T<:AbstractHyperGraph}
    idx = vertices_to_indices(x)
    I = zeros(Int, nv(x), nhe(x))
    @inbounds for (i, e) in enumerate(hyperedges(x))
        [I[idx[v], i] -= s for (v, s) in zip(objects(e.src), src_multiplicities(e))]
        [I[idx[v], i] += s for (v, s) in zip(objects(e.tgt), tgt_multiplicities(e))]
    end
    I
end

# incidence matrix (unoriented case)
@traitfn function incidence_matrix(x::T::(!IsOriented)) where {T<:AbstractHyperGraph}
    idx = HyperGraphs.vertices_to_indices(x)
    I = zeros(Int, nv(x), nhe(x))
    @inbounds for (i, e) in enumerate(hyperedges(x))
        [I[idx[v], i] += s for (v, s) in zip(vertices(e), fill(1, cardinality(e)))]
    end
    I
end

# cardinality matrix
cardinality_matrix(x::T) where {T<:AbstractHyperGraph} = Diagonal(cardinalities(x))

# weight matrix
weight_matrix(x::T) where {T<:AbstractHyperGraph} = Diagonal(weights(x))


text:
@testset "Document" begin
    @testset "count unicode literals in bytes" begin
        s = """
        \"""
        𝔽𝔽

        \"""
        struct A end
        """
        d = JuliaFormatter.Document(s)
        ranges = Dict(1 => 1:4, 2 => 5:7, 3 => 8:8, 4 => 9:12, 5 => 13:25, 6 => 26:25)
        @test ranges == d.line_to_range
    end

    @testset "count unicode whitespace in bytes" begin
        s0 = """a = b || c ;
               f("A")"""
        d = JuliaFormatter.Document(s0)
        ranges = Dict(1 => 1:13, 2 => 14:19)
        @test ranges == d.line_to_range

        # this string has a nbsp after 'c'
        # so it should have an additional byte because
        # it's unicode
        s = """a = b || c ;
               f("A")"""
        d = JuliaFormatter.Document(s)
        ranges = Dict(1 => 1:14, 2 => 15:20)
        @test ranges == d.line_to_range
    end
end


text:
# tempoptions() = Dict(
#     :gisregion => "Europe8",            # "Europe8", "Eurasia38", "Scand3"
#
#     :scenarioyear => "ssp2_2050",       # default scenario and year for population and grid access datasets
#     :era_year => 2018,                  # which year of the ERA5 time series to use
#
#     :numcenters => 3                    # number of population centers in each region (for average hourly temperature)
#     :mindist => 3.3                     # minimum distance between population centers [in ERA5 pixels]
# )

function GIStemp(gisregion::String, scenarioyear::String, era_year::Int, numcenters::Int, mindist::Float64)
    println("\nReading temperature data for $gisregion...")

    regions, offshoreregions, regionlist, lonrange, latrange, pop, meantemp, temp, eralonranges, eralatrange =
                read_temperature_datasets(gisregion, scenarioyear, era_year)

    erapop = rescale_population_to_ERA5_res(era_year, regions, offshoreregions, regionlist, lonrange, latrange, pop, temp)
    tempmask = dropdims(sum(abs.(temp[1000:1100,:,:]), dims=1), dims=1) .> 0    # detect cells with any temp data (coasts?)

    numreg = length(regionlist)
    popcenters = [findpopcenters(erapop[i,:,:].*tempmask, numcenters, mindist) for i = 1:numreg]
    # if any region lacks enough population centers, just repeat the ones found (only happens for tiny regions, e.g. Malta)
    for (i,p) in enumerate(popcenters)
        length(p) == 0 && error("No population centers found for region $(regionlist[i]).")
        if length(p) < numcenters
            popcenters[i] = repeat(p, outer=numcenters)[1:numcenters]
        end
    end

    eralonrange = vcat(eralonranges...)
    erares = 0.28125
    eralons = (-180+erares/2:erares:180-erares/2)[eralonrange]     # longitude values (pixel center)
    eralats = (90-erares/2:-erares:-90+erares/2)[eralatrange]      # latitude values (pixel center)

    hours = DateTime(era_year, 1, 1, 0) : Hour(1) : DateTime(era_year, 12, 31, 23)
    numhours = length(hours)
    numhours != size(temp,1) && error("Inconsistent number of hours.")

    temp_popcenters = [temp[h, popcenters[r][i]...] for h = 1:numhours, r = 1:numreg, i = 1:numcenters]

    return hours, temp_popcenters
end

function read_temperature_datasets(gisregion, scenarioyear, era_year)
    res = 0.01          # resolution of auxiliary datasets [degrees per pixel]
    erares = 0.28125    # resolution of ERA5 datasets [degrees per pixel]

    println("\nReading auxiliary datasets...")
    regions, offshoreregions, regionlist, lonrange, latrange = loadregions(gisregion)
    pop = JLD.load(in_datafolder("population_$scenarioyear.jld"), "population")[lonrange,latrange]

    println("Reading ERA5 temperature dataset...")
    eralonranges, eralatrange = eraranges(lonrange, latrange, res, erares)

    @time meantemp, temp = h5open(in_datafolder("era5temp$era_year.h5"), "r") do file
        if length(eralonranges) == 1
            file["meantemp"][eralonranges[1], eralatrange],
                file["temp"][:,eralonranges[1], eralatrange]
        else
            [file["meantemp"][eralonranges[1], eralatrange]; file["meantemp"][eralonranges[2], eralatrange]],
                [file["temp"][:, eralonranges[1], eralatrange] file["temp"][:, eralonranges[2], eralatrange]]
        end
    end
    return regions, offshoreregions, regionlist, lonrange, latrange, pop, meantemp, temp, eralonranges, eralatrange
end

function rescale_population_to_ERA5_res(era_year, regions, offshoreregions, regionlist, lonrange, latrange, pop, temp)
    res = 0.01          # resolution of auxiliary datasets [degrees per pixel]
    erares = 0.28125    # resolution of ERA5 datasets [degrees per pixel]

    println("Calculating population centers at ERA5 resolution...")
    eralons, eralats, lonmap, latmap, cellarea = eralonlat(Dict(:res=>res, :erares=>erares), lonrange, latrange)

    numreg = length(regionlist)
    yearlength, nlons, nlats = size(temp)
    firsttime = DateTime(era_year, 1, 1)

    erapop = zeros(numreg, nlons, nlats)

    # Run times vary wildly depending on geographical area (because of far offshore regions with mostly zero wind speeds).
    # To improve the estimated time of completing the progress bar, iterate over latitudes in random order.
    Random.seed!(1)
    updateprogress = Progress(nlats, 1)
    for j in randperm(nlats)
        eralat = eralats[j]
        colrange = latmap[lat2col(eralat+erares/2, res):lat2col(eralat-erares/2, res)-1]
        for i = 1:nlons
            eralon = eralons[i]
            # get all high resolution row and column indexes within this ERA5 cell
            rowrange = lonmap[lon2row(eralon-erares/2, res):lon2row(eralon+erares/2, res)-1]

            for c in colrange, r in rowrange
                (c == 0 || r == 0) && continue
                reg = regions[r,c]
                if reg > 0 && reg != NOREGION
                    erapop[reg,i,j] += pop[r,c]
                end
            end
        end
        next!(updateprogress)
    end

    return erapop
end

# Find the (row,col) coordinates of the n largest population cells which are at least mindist cells apart from each other.
function findpopcenters(pop, n, mindist)
    n = min(n, sum(pop .> 0))   # make sure there are enough non-zero cells
    p = copy(pop)
    coords = CartesianIndex{2}[]
    for i = 1:n
        val, index = findmax(p)
        if val == 0   # if the largest cell has zero pop then we were too ambitious - redo with lower distance
            # println("lowering distance")
            return findpopcenters(pop, n, mindist/2)
        end
        push!(coords, index)
        fillcircle!(p, Tuple(index), mindist, 0)
    end
    # display([dist(coords[i], coords[j]) for i = 1:length(coords)-1, j = 2:length(coords) if j > i])
    return Tuple.(coords)
end

dist(a::CartesianIndex, b::CartesianIndex) = sqrt(sum((Tuple(a) .- Tuple(b)).^2))

function fillcircle!(a::AbstractMatrix, center::Tuple{Int,Int}, radius::Real, fillvalue)
    nrows, ncols = size(a)
    rowrange = max(1, ceil(Int, center[1] - radius)) : min(nrows, floor(Int, center[1] + radius))
    colrange = max(1, ceil(Int, center[2] - radius)) : min(ncols, floor(Int, center[2] + radius))
    for c in colrange, r in rowrange
        if (r - center[1])^2 + (c - center[2])^2 <= radius^2
            a[r,c] = fillvalue
        end
    end
end

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
"""
    struct PolynomialV <: AbstractVApproximator{Int}
        weights::Vector{Float64}
    end

See more details at Section (9.5.1) on Page 210 of the book *Sutton, Richard S., and Andrew G. Barto. Reinforcement learning: An introduction. MIT press, 2018.*
"""
struct PolynomialV <: AbstractVApproximator{Int}
    weights::Vector{Float64}
end

PolynomialV(order::Int) = PolynomialV(zeros(Float64, order+1))

function (p::PolynomialV)(s)
    sum(w * s^(i-1) for (i, w) in enumerate(p.weights))
end

function update!(p::PolynomialV, s, e)
    for i in 1:length(p.weights)
        p.weights[i] += e * s^(i-1)
    end
end


text:
module ConstantTests

using Interpolations
using Base.Test

# Instantiation
N1 = 10
A1 = rand(Float64, N1) * 100
A2 = rand(Float64, N1, N1) * 100
A3 = rand(Float64, N1, N1, N1) * 100

for (constructor, copier) in ((interpolate, x->x), (interpolate!, copy))
    itp1c = @inferred(constructor(copier(A1), BSpline(Constant()), OnCell()))
    itp1g = @inferred(constructor(copier(A1), BSpline(Constant()), OnGrid()))
    itp2c = @inferred(constructor(copier(A2), BSpline(Constant()), OnCell()))
    itp2g = @inferred(constructor(copier(A2), BSpline(Constant()), OnGrid()))
    itp3c = @inferred(constructor(copier(A3), BSpline(Constant()), OnCell()))
    itp3g = @inferred(constructor(copier(A3), BSpline(Constant()), OnGrid()))

    # Evaluation on provided data points
    # 1D
    for i in 1:length(A1)
        @test A1[i] == itp1c[i] == itp1g[i]
        @test A1[i] == itp1c[convert(Float64,i)] == itp1g[convert(Float64,i)]
    end
    @test @inferred(size(itp1c)) == size(A1)
    @test @inferred(size(itp1g)) == size(A1)
    # 2D
    for i in 1:N1, j in 1:N1
        @test A2[i,j] == itp2c[i,j] == itp2g[i,j]
        @test A2[i,j] == itp2c[convert(Float64,i),convert(Float64,j)] == itp2g[convert(Float64,i),convert(Float64,j)]
    end
    @test @inferred(size(itp2c)) == size(A2)
    @test @inferred(size(itp2g)) == size(A2)
    # 3D
    for i in 1:N1, j in 1:N1, k in 1:N1
        @test A3[i,j,k] == itp3c[i,j,k] == itp3g[i,j,k]
        @test A3[i,j,k] == itp3c[convert(Float64,i),convert(Float64,j),convert(Float64,k)] == itp3g[convert(Float64,i),convert(Float64,j),convert(Float64,k)]
    end
    @test @inferred(size(itp3c)) == size(A3)
    @test @inferred(size(itp3g)) == size(A3)

    # Evaluation between data points
    for i in 2:N1-1
        @test A1[i] == itp1c[i+.3] == itp1g[i+.3] == itp1c[i-.3] == itp1g[i-.3]
    end
    # 2D
    for i in 2:N1-1, j in 2:N1-1
        @test A2[i,j] == itp2c[i+.4,j-.3] == itp2g[i+.4,j-.3]
    end
    # 3D
    for i in 2:N1-1, j in 2:N1-1, k in 2:N1-1
        @test A3[i,j,k] == itp3c[i+.4,j-.3,k+.1] == itp3g[i+.4,j-.3,k+.2]
    end

    # Edge behavior
    @test A1[1] == itp1c[.7]
    @test A1[N1] == itp1c[N1+.3]
end

end


text:
module FileMgrTest
  using Base.Test
  using HelloDB.FileMgrs
  using HelloDB

  function test_init()
    dbname = "some_test_db"
    dbdir = homedir() * "/$dbname"
    @test !isdir(dbdir)
    FileMgr(dbname)
    @test isdir(dbdir)
    @test is_dbdir(dbdir)
    if is_dbdir(dbdir)
      rm(dbdir; recursive=true)
    end
  end

  test_init()
end


text:
export Constraint, <<, >>, ==

"""
    Constraint(p::AbstractPolynomial, lb::Number, ub::Number)

Define a mathematical complex constraint, from a polynomial body and two complex
bounds. Equivalent to the two real inequality produced from real and imag part
of `p`, `lb` and `ub`.

### Attributes
- `p::Polynomial` : the polynomial body of the constraint.
- `lb::Number` : complex lower bound.
- `ub::Number` : complex upper bound.
- `precond::Symbol` : a symbol describing the preconditioning method to be
applied for this constraint. Default value is `:none`.

### Exemple
```julia
julia > a, b = Variable("a", Complex), Variable("b", Complex)
julia > cstr = abs2(a^2) + abs2(3*b) << 25
julia > cstr.precond = :sqrt
```
"""
mutable struct Constraint
  p::Polynomial
  lb::Number
  ub::Number
  precond::Symbol

  Constraint(p::AbstractPolynomial, lb::Number, ub::Number) = new(p, lb, ub, :none)
end

#############################
## Constraint constructors
#############################
function <<(p::T, bound::Number) where T<: AbstractPolynomial
  return Constraint(p, -Inf-im*Inf, bound)
end

function <<(bound::Number, p::T) where T<: AbstractPolynomial
  return Constraint(p, bound, Inf+im*Inf)
end

function >>(p::T, bound::Number) where T<:AbstractPolynomial
  return bound << p
end

function >>(bound::Number, p::T) where T<:AbstractPolynomial
  return p << bound
end

function <<(bound::Number, cstr::Constraint)
  if real(bound) > real(cstr.ub) || imag(bound) > imag(cstr.ub)
    warn(LOGGER, "<<(bnd, cstr): Creating a constraint with lower bound ", bound, " and upper bound ", cstr.ub)
  end
  return Constraint(cstr.p, bound, cstr.ub)
end
>>(bound::Number, cstr::Constraint) = cstr << bound

function <<(cstr::Constraint, bound::Number)
  if real(cstr.lb) > real(bound) || imag(cstr.lb) > imag(bound)
    warn(LOGGER, "<<(cstr, bnd): Creating a constraint with lower bound ", cstr.lb, " and upper bound ", bound)
  end
  return Constraint(cstr.p, cstr.lb, bound)
end
>>(cstr::Constraint, bound::Number) = bound << cstr

function ==(p::T, x::Number) where T<:AbstractPolynomial
  return Constraint(p, x, x)
end

#############################
## Equality, hash
#############################
==(c1::Constraint, c2::Constraint) = (c1.p == c2.p) &&
                                     (c1.ub == c2.ub) &&
                                     (c1.lb == c2.lb) &&
                                     (c1.precond == c2.precond)

!=(c1::Constraint, c2::Constraint) = !(c1 == c2)

#############################
## Print
#############################
function Base.print(io::IO, cstr::Constraint)
  if cstr.lb == cstr.ub
    print(io, cstr.p, " = ", cstr.ub)
  else
    if cstr.lb != (-Inf-im*Inf)
      print(io, cstr.lb, " < ")
    end
    print(io, cstr.p)
    if cstr.ub !=(Inf+im*Inf)
      print(io, " < ", cstr.ub)
    end
  end
end


text:
#Seed

Random.seed!(1234);

#iesh or electronic friction, set 0 for iesh and 1 for electronic friction
global const latticeopt = 0

#conversion factors to simulation units.

global const conv1 = 6.02214076 # Newton/meter * Angstrom^2 in kilojoule/mole
global const conv2 = 96.48 # eV in kilojoule per mole
global const conv3 = 1e-4 #kilojoule per mole in amu*angstrom^2/femtosecond^2
#constants in SI units
global const Å = 1.0e-10 # 1 Angström in meter
global const ev = 1.60217646e-19 # 1 eV in Joule
global const amu = 1.66053886e-27	# 1 amu in kg
global const hbar = 6.582119569e-16*conv2*conv3*1e15	# hbar in amu*angstrom^2/femtosecond^2*fs
global const fs = 1.0e-15  	# fs in s

global const kjpermol = 1.660538863e-21 # 1 kJ/mol in Joules
global const N_a = 602214076000000000000000

global const kb = 0.08617*1e-3 * conv2 * conv3 	# boltzman constant [mev/Kelvin] in simulation units []

#Simulation parameters
const global Ne = 20     #number of electrons
const global Ms = 40     #number of orbitals
const global numtraj = 1   #number of trajectories
const global tsteps = 1000
const global thop = 1   #number of timesteps between surface hops
const global twrite = 5     #number of timesteps between data writing
const global dt = 0.5  #time stepsize in Femtoseconds
#these energy parameters are used in subsequent calculations with other variables defined
# in the simulation units Angstrom, amu, femtoseconds. Therefore we need to convert
# kilojoule/mole to amu*angstrom^2/femtosecond^2
const global r_no = 1.15077		# Initial bond length of NO
const global e_trans_i = 0.05*conv2*conv3 #translational energy in eV -> kilojoule/mole
const global e_vib_i = 3.1*conv2*conv3 #vibrational energy in eV-> kilojoule/mole
const global e_rot_i = 0.05*conv2*conv3 #translational energy in eV-> kilojoule/mole
const global phi_inc = 0.0 #incident angle in radians
const global tsurf = 300.0 #temperature of surface in Kelvin
const global sim_time = tsteps * dt

#parameters of surface etc.
const global delta_E = 7.0 *conv2*conv3
const global z_end = 11.0

#Parameters of Neutral Diabatic Matrix Element
global const A_n = 457000.052788*conv3 #Au--O: exponential repulsion: A
global const α_n = 3.75257871821 #Au--O: exponential repulsion: alpha
global const B_n = 30788.8486039*conv3 #Au--N: exponential repulsion: B
global const β_n = 2.9728905911 #Au--N: exponential repulsion: beta
global const cutoff = 10.0 #Au--N: cutoff
global const exp_beta_n_cutoff = exp(-β_n * cutoff)
global const exp_alpha_n_cutoff = exp(-α_n * cutoff)
global const F_n = 638.5000000000*conv3 #N--O: Morse: F
global const δ_n=2.7430000000 #N--O: Morse: delta
global const r_0_NO = 1.1507700000#N--O: Morse: cutoff

#Parameters of IONIC Diabatic Matrix Element

global const C_i = 1.25581276843  # Image potential: C
global const D_i = 347.2225355*conv3  #Image potential: D
global const z0_i = 1.153606314 #Image potential: z0
global const A_i = 457000.052788*conv3  #Au--O: exponential repulsion: A
global const α_i = 3.75257871821  #Au--O: exponential repulsion: alpha
global const B_i = 23.8597594272*conv3  #Au--N: Morse: B
global const β_i=1.91014033785  #Au--N: Morse: beta
global const rN_sur_e_i = 2.38958832878 #Au--N: Morse: rN_sur_e
global const F_i = 495.9809807*conv3   #N--O: Morse: F
global const δ_i = 2.47093477934  #N--O: Morse: delta
global const rNO_e_i = 1.29289837288  #N--O: Morse: cutoff
global const KI = 512.06425722*conv3  #N--O: Morse: cutoff
global const exp_alpha_i_cutoff = exp(-α_i * cutoff)
global const exp_beta_i_cutoff1 = exp(-2*β_i * (cutoff - rN_sur_e_i))
global const exp_beta_i_cutoff2 = exp(-β_i * (cutoff - rN_sur_e_i))
# DEFINE PARAMETERS FOR COUPLING FUNCTION

global const coup_a_N=-70.5259924491*conv3 # Au--N: Exponential decay: a
global const coup_b_N=0.00470023958504 # Au--N: Exponential decay: b
global const coup_β_N=1.95982478112 # Au--N: Exponential decay: beta
global const coup_a_O=-16.7488672932*conv3 # Au--O: Exponential decay: a
global const coup_b_O=0.00617151653727 # Au--O: Exponential decay: b
global const coup_β_O=1.35353579356  # Au--O: Exponential decay: beta
global const Au_O_coupling_cutoff=10.0 # Au--O: Exponential decay: cutoff
global const coup_cutoff_O = coup_a_O/(1+coup_b_O*exp(coup_β_O * Au_O_coupling_cutoff))
global const coup_cutoff_N = coup_a_N/(1+coup_b_N*exp(coup_β_N * Au_O_coupling_cutoff))
#Define parameters for Au-Au interaction potential

# in Newton/meters, multiply with conv1 so that resulting units of energy will
# be in kilojoule per mole
global const α = -4.94*conv1
global const β = 17.15*conv1
global const γ = 19.40*conv1

function def_d_matrices(α, β, γ)
    #see paper from 1947 for definitions
    d17 = [[α 0.0 0.0]; [0.0 β γ]; [0.0 γ β]]
    d28 = [[α 0.0 0.0]; [0.0 β -γ]; [0.0 -γ β]]
    d39 = [[β 0.0 γ]; [0.0 α 0.0]; [γ 0.0 β]]
    d410 = [[β 0.0 -γ]; [0.0 α 0.0]; [-γ 0.0 β]]
    d511 = [[β γ 0.0]; [γ β 0.0]; [0.0 0.0 α]]
    d612 = [[β -γ 0.0]; [-γ β 0.0]; [0.0 0.0 α]]

    return d17, d28, d39, d410, d511, d612
end
const global d17, d28, d39, d410, d511, d612 = def_d_matrices(α, β, γ)

function compute_d_new_basis(d17, d28, d39, d410, d511, d612)
    U = permutedims([[-1.0 0.0 1.0]/sqrt(2.0);[1.0 -2.0 1.0]/sqrt(6.0);[-1.0 -1.0 -1.0]/sqrt(3.0)])
    d1_new = U' * d17
    d2_new = U' * d28
    d3_new = U' * d39
    d4_new = U' * d410
    d5_new = U' * d511
    d6_new = U' * d612
    return d1_new, d2_new, d3_new, d4_new, d5_new, d6_new
end
const global d1_new, d2_new, d3_new, d4_new, d5_new, d6_new = compute_d_new_basis(d17, d28, d39, d410, d511, d612)

#set masses
const global mass_arr =  [14.00307440, 15.99491502, 196.966548]
const global m_N = mass_arr[1]
const global m_O = mass_arr[2]
const global m_au = mass_arr[3]
const global μ = m_N*m_O/(m_N + m_O)


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

mutable struct PublicIPAddressSku <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name

    function PublicIPAddressSku(;name=nothing)
        o = new()
        validate_property(PublicIPAddressSku, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        o
    end
end # type PublicIPAddressSku

const _property_map_PublicIPAddressSku = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"))
const _property_types_PublicIPAddressSku = Dict{Symbol,String}(Symbol("name")=>"String")
Base.propertynames(::Type{ PublicIPAddressSku }) = collect(keys(_property_map_PublicIPAddressSku))
Swagger.property_type(::Type{ PublicIPAddressSku }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_PublicIPAddressSku[name]))}
Swagger.field_name(::Type{ PublicIPAddressSku }, property_name::Symbol) =  _property_map_PublicIPAddressSku[property_name]

const _allowed_PublicIPAddressSku_name = ["Basic", "Standard"]

function check_required(o::PublicIPAddressSku)
    true
end

function validate_property(::Type{ PublicIPAddressSku }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "PublicIPAddressSku", :enum, val, _allowed_PublicIPAddressSku_name)
    end
end


text:
module MyInterpolations

export lin_int

function lin_int(grid,vals)
    function func(x)
        if x < grid[1]
            return "error"
        end
        if x > grid[end]
            return "error"
        end
        if x == grid[end]
            return vals[end]
        end
        if grid[1] <= x < length(grid)
            z = searchsortedlast(grid,x)
            slope = (vals[z+1]-vals[z])/(grid[z+1]-grid[z])
            y = slope*(x-grid[z])+vals[z]
        return y
        end
    end
    return func
end

end


text:
using DSGE
using BenchmarkTools, JLD
import DSGE: jacobian

path = dirname(@__FILE__)
filepath = "$path/../../reference/KrusellSmith/jacobian.jld"

# Benchmark current function
m = KrusellSmith()
jacobian(m)
trial = @benchmark jacobian($m) gcsample = true

# Optionally over-write the existing reference trial
# write_ref_trial(trial, "jacobian", filepath = filepath)

# Compute time differential
print_all_benchmarks(trial, filepath, "jacobian")