text:
"""
```
Setting{T}
```

The `Setting` type is an interface for computational settings that affect how the code runs
without affecting the mathematical definition of the model. It also provides support for
non-conflicting file names for output of 2 models that differ only in the values of their
computational settings.

### Fields
- `key::Symbol`: Name of setting
- `value::T`: Value of setting
- `print::Bool`: Indicates whether to append this setting's code and value to output file
  names. If true, output file names will include a suffix of the form _code1=val1_code2=val2
  etc. where codes are listed in alphabetical order.
- `code::AbstractString`: string of <=4 characters to print to output file suffixes when
  `print=true`.
- `description::AbstractString`: Short description of what the setting is used for.
"""
type Setting{T}
    key::Symbol                  # name of setting
    value::T                     # whatever the setting is
    print::Bool                  # whether or not to add this setting to the print
    code::AbstractString         # what gets printed to the print
    description::AbstractString  # description of what the setting is for
end

# for printing codes to filename string
Base.convert{T<:Number, U<:Number}(::Type{T}, s::Setting{U}) = convert(T, s.value)
Base.convert{T<:AbstractString, U<:AbstractString}(::Type{T}, s::Setting{U}) = convert(T, s.value)

Base.promote_rule{T<:Number,U<:Number}(::Type{Setting{T}}, ::Type{U}) = promote_rule(T,U)
Base.promote_rule{T<:AbstractString,U<:AbstractString}(::Type{Setting{T}}, ::Type{U}) = promote_rule(T,U)
Base.promote_rule(::Type{Setting{Bool}}, ::Type{Bool}) = promote_rule(Bool, Bool)

Base.string(s::Setting{AbstractString}) = string(s.value)

to_filestring(s::Setting) = "$(s.code)=$(s.value)"

# key, value constructor
Setting(key, value) = Setting(key, value, false, "", "")
Setting(key, value, description) = Setting(key, value, false, "", description)

function Base.show(io::IO, s::Setting)
    @printf io "%s: %s" s.key s.value
end

"""
```
(<=)(m::AbstractModel, s::Setting)
```
Syntax for adding a setting to a model/overwriting a setting via `m <= Setting(...)`
"""
function (<=)(m::AbstractModel, s::Setting)
    if !m.testing
        setting_field_name = :settings
    else
        setting_field_name = :test_settings
    end

    if !haskey(getfield(m, setting_field_name), s.key)
        getfield(m, setting_field_name)[s.key] = s
    else
        update!(getfield(m, setting_field_name)[s.key], s)
    end
end

function update!(a::Setting, b::Setting)
    # Make sure Setting a can appropriately be updated by b.
    if a.key ≠ b.key
        return
    end

    a.value = b.value

    # b overrides the print boolean of a if:
    # - a.print is false and b.print is true.
    # - a.print is true, b.print is false, and b.code is non-empty.
    # We must be able to tell if b was created via a constructor like `Setting(:key,
    # value)`, in which case the print, code, and description values are set to defaults. We
    # do not overwrite if we can't determine whether or not those fields are just the
    # defaults.
    if (a.print == false && b.print == true) ||
       (a.print == true && b.print == false && !isempty(b.code))
       a.print = b.print
   end
    if !isempty(b.code) && b.code ≠ a.code
        a.code = b.code
    end
    if !isempty(b.description) && b.description ≠ a.description
        a.description = b.description
    end
end

"""
```
get_setting(m::AbstractModel, setting::Symbol)
```

Returns the value of the setting
"""
function get_setting(m::AbstractModel, s::Symbol)

    if m.testing && in(s, keys(m.test_settings))
        return m.test_settings[s].value
    end

    return m.settings[s].value
end

"""
```
default_test_settings!(m::AbstractModel)
```

The following Settings are constructed, initialized and added to
`m.test_settings`. Their purposes are identical to those in
`m.settings`, but these values are used to test DSGE.jl.

### I/O Locations and identifiers
- `saveroot::Setting{ASCIIString}`: A temporary directory in /tmp/
- `dataroot::Setting{ASCIIString}`: dsgeroot/test/reference/
- `data_vintage::Setting{ASCIIString}`: "_REF"

### Metropolis-Hastings
- `n_mh_simulations::Setting{Int}`: 100
- `n_mh_blocks::Setting{Int}`: 1
- `n_mh_burn::Setting{Int}`: 0
- `mh_thin::Setting{Int}`: 1
"""
function default_test_settings!(m::AbstractModel)

    test = m.test_settings

    # I/O
    dataroot = normpath(joinpath(dirname(@__FILE__), "..","test","reference"))
    saveroot = mktempdir()

    #General
    test[:saveroot] = Setting(:saveroot, saveroot,
        "Where to write files when in test mode")
    test[:dataroot] = Setting(:dataroot, dataroot,
        "Location of input files when in test mode" )
    test[:data_vintage] = Setting(:data_vintage, "REF", true, "vint",
        "Reference data identifier")
    test[:date_mainsample_end] = Setting(:date_mainsample_end, quartertodate("2015-Q3"),
        "End date of main sample")
    test[:use_parallel_workers] = Setting(:use_parallel_workers, false, false, "parw",
        "Use available parallel workers in computations")
    test[:n_hessian_test_params] = Setting(:n_hessian_test_params, 3, false, "mhfp",
        "Max number of free params for which to calculate Hessian")

    # Metropolis-Hastings
    test[:n_mh_simulations] = Setting(:n_mh_simulations, 100, false, "nsim",
        "Number of parameter draws per block for testing Metropolis-Hastings")
    test[:n_mh_blocks] = Setting(:n_mh_blocks, 1, false, "nblc",
        "Number of blocks to draw parameters for testing Metropolis-Hastings")
    test[:n_mh_burn] = Setting(:n_mh_burn, 0, false, "nbrn",
        "Number of burn-in blocks for testing Metropolis-Hastings")
    test[:mh_thin] = Setting(:mh_thin, 1, false, "thin",
        "Thinning step for testing Metropolis-Hastings")

    return test
end


text:
using HyperbolicPlane, Plots

"""
`dual_tess()` draws dual tesselations of the hyperbolic plane.
The primal is a tiling by pentagons (drawn with solid lines)
and the dual is a tiling by quadrilaterals (drawn dotted).
"""
function dual_tess()
    C = tesselation(5,4,5)  # primal
    D = tesselation(4,5,5,true)  # dual

    # rotate all the pentagons 36 degrees
    r = rotation(pi/5)
    CC = [r(P) for P in C.objs]

    # make the dual quadrilaterals dotted
    for X in D.objs
        set_line_style(X,:dot)
    end

    plot()   # clear the screen
    draw(CC) # draw primal
    draw(D)  # draw dual
    finish()
end

ext: jl


text:
"""
    struct CompressedArray{T,N,A,P} <: AbstractArray{T,N}
      values::A
      ptrs::P
    end
Type representing an array with a reduced set of values.
The array is represented by a short array of values, namely
the field `values`, and a large array of indices, namely the
field `ptrs`. The `i`-th component of the resulting array is
defined as `values[ptrs[i]]`. The type parameters `A`, and `P`
are restricted to be array types by the inner constructor of this `struct`.
"""
struct CompressedArray{T,N,A,P} <: AbstractArray{T,N}
  values::A
  ptrs::P
  @doc """
      CompressedArray(values::AbstractArray,ptrs::AbstractArray)

  Creates a `CompressedArray` object by the given arrays of `values` and
  `ptrs`.
  """
  function CompressedArray(values::AbstractArray,ptrs::AbstractArray)
    A = typeof(values)
    P = typeof(ptrs)
    T = eltype(values)
    N = ndims(ptrs)
    new{T,N,A,P}(values,ptrs)
  end
end

function testitem(a::CompressedArray)
  if length(a.ptrs) == 0
    testitem(a.values)
  else
    a.values[first(a.ptrs)]
  end
end

size(a::CompressedArray) = size(a.ptrs)

@propagate_inbounds function getindex(a::CompressedArray,i::Integer)
  j = a.ptrs[i]
  a.values[j]
end

@propagate_inbounds function getindex(a::CompressedArray{T,N},i::Vararg{Integer,N}) where {T,N}
  j = a.ptrs[i...]
  a.values[j]
end

function IndexStyle(a::Type{CompressedArray{T,N,A,P}}) where {T,N,A,P}
  IndexStyle(P)
end

function lazy_map(::typeof(evaluate),::Type{T},g::CompressedArray...) where T
  if _have_same_ptrs(g)
    _lazy_map_compressed(g...)
  else
    LazyArray(T,g...)
  end
end

function lazy_map(::typeof(evaluate),::Type{T},g::Union{CompressedArray,Fill}...) where T
  g_compressed = _find_compressed_ones(g)
  if _have_same_ptrs(g_compressed)
    g1 = first(g_compressed)
    g_all_compressed = map(gi->_compress(gi,g1),g)
    _lazy_map_compressed(g_all_compressed...)
  else
    LazyArray(T,g...)
  end
end

function _find_compressed_ones(g)
  g_compressed = ( gi for gi in g if isa(gi,CompressedArray) )
  g_compressed
end

function _lazy_map_compressed(g::CompressedArray...)
  vals = map(evaluate, map(gi->gi.values,g)...)
  ptrs = first(g).ptrs
  CompressedArray(vals,ptrs)
end

function _have_same_ptrs(g)
  g1 = first(g)
  all(map( gi -> gi.ptrs === g1.ptrs ||  gi.ptrs == g1.ptrs, g))
end

function _compress(a::CompressedArray,b::CompressedArray)
  @check _have_same_ptrs((a,b))
  a
end

function _compress(a::Fill,b::CompressedArray)
  vals = fill(a.value,length(b.values))
  CompressedArray(vals,b.ptrs)
end


text:
function add_dual_equality_constraints(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_dual_map::PrimalDualMap,
    dual_names::DualNames,
    primal_objective::PrimalObjective{T},
    con_types::Vector{Tuple{Type,Type}},
    variable_parameters::Vector{VI},
) where {T}
    sense_change = ifelse(
        MOI.get(primal_model, MOI.ObjectiveSense()) == MOI.MIN_SENSE,
        one(T),
        -one(T),
    )

    all_variables = MOI.get(primal_model, MOI.ListOfVariableIndices())
    non_parameter_variables = setdiff(all_variables, variable_parameters)

    # Loop at every constraint to collect the scalar affine terms in the
    # `scalar_affine_terms` list.
    # TODO: flip these signs a priorie instead of require post processing later
    scalar_affine_terms = get_scalar_affine_terms(
        primal_model,
        primal_dual_map.primal_con_dual_var,
        all_variables,
        con_types,
        T,
    )

    # get constants (rhs) of dual constraints from primal objective coefficients
    scalar_terms = get_scalar_terms(primal_objective)

    # Collect affine terms of dual constraints that come from the quadratic
    # part of the primal objective function, and add them into
    # `scalar_affine_terms`.
    # These terms are added with flipped signs (because the sign will be flipped again).
    # TODO: unflip these signs
    add_scalar_affine_terms_from_quad_obj(
        scalar_affine_terms,
        primal_dual_map.primal_var_dual_quad_slack,
        primal_objective,
    )

    # terms from mixing variables and parameters
    # These terms are added with flipped signs (because the sign will be flipped again).
    # TODO: unflip these signs
    add_scalar_affine_terms_from_quad_params(
        scalar_affine_terms,
        primal_dual_map.primal_parameter,
        primal_objective,
    )

    # Constrained variables
    for ci in keys(primal_dual_map.constrained_var_dual)
        # Add constraints associated with constrained variables
        # These are constraints that will not be regular equality constraints
        # they will be  function-in-set, where set is the dual set of the
        # primal constrained variable.
        # TODO: flip these signs a priori instead of requiring post-processing later
        _add_constrained_variable_constraint(
            dual_model,
            primal_model,
            primal_dual_map.constrained_var_zero,
            primal_dual_map.constrained_var_dual,
            ci,
            scalar_affine_terms,
            scalar_terms,
            sense_change,
            T,
            dual_names,
        )
    end

    # Free variables
    for primal_vi in non_parameter_variables
        if primal_vi in keys(primal_dual_map.constrained_var_idx)
            continue # constrained variable
        end
        # Add equality constraint
        # these are constraints associated to primal variables that# are not
        # treated as constrained variables, that is "free variables" (x \in R)
        # therefore their associated dual constraints are equalities.
        dual_ci = MOIU.normalize_and_add_constraint(
            dual_model,
            MOI.ScalarAffineFunction(
                # TODO: flip these two signs bellow to mathc _add_constrained_variable_constraint
                # MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
                # sense_change * get(scalar_terms, primal_vi, zero(T))),
                MOIU.operate_terms(+, scalar_affine_terms[primal_vi]),
                -sense_change * get(scalar_terms, primal_vi, zero(T)),
            ),
            MOI.EqualTo(zero(T)),
        )
        # Set constraint name with the name of the associated priaml variable
        if !is_empty(dual_names)
            set_dual_constraint_name(
                dual_model,
                primal_model,
                primal_vi,
                dual_ci,
                dual_names.dual_constraint_name_prefix,
            )
        end
        # Add primal variable to dual contraint to the link dictionary
        push!(primal_dual_map.primal_var_dual_con, primal_vi => dual_ci)
    end
    return scalar_affine_terms
end

# `_add_constrained_variable_constraint` is a function barrier to avoid type
# unstability.

function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables,MOI.Zeros},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
) where {T}
    # The dual is `Reals`, adding a constraint `func`-in-`Reals` is equivalent
    # to not adding any constraint.
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    zero_map[ci] = MOIU.vectorize([
        MOI.ScalarAffineFunction(
            MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
            sense_change * get(scalar_terms, primal_vi, zero(T)),
        ) for primal_vi in func_primal.variables
    ])
    return
end
function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
    dual_names::DualNames,
) where {T}
    set_primal = MOI.get(primal_model, MOI.ConstraintSet(), ci)
    set_dual = MOI.dual_set(set_primal)
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    func_dual = MOIU.vectorize([
        MOI.ScalarAffineFunction(
            MOIU.operate_term.(
                *,
                -inv(set_dot(i, set_primal, T)),
                scalar_affine_terms[primal_vi],
            ),
            sense_change *
            inv(set_dot(i, set_primal, T)) *
            get(scalar_terms, primal_vi, zero(T)),
        ) for (i, primal_vi) in enumerate(func_primal.variables)
    ])
    ci_map[ci] = MOI.add_constraint(dual_model, func_dual, set_dual)
    if !is_empty(dual_names)
        @warn(
            "dual names for constrained vector of variables not supported yet."
        )
    end
    return
end

function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VariableIndex,<:MOI.EqualTo},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
    dual_names::DualNames,
) where {T}
    # Nothing to add as the set is `EqualTo`.
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    primal_vi = func_primal
    zero_map[ci] = MOI.ScalarAffineFunction(
        MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
        sense_change * get(scalar_terms, primal_vi, zero(T)),
    )
    return
end
function _add_constrained_variable_constraint(
    dual_model,
    primal_model,
    zero_map,
    ci_map,
    ci::MOI.ConstraintIndex{MOI.VariableIndex},
    scalar_affine_terms,
    scalar_terms,
    sense_change,
    ::Type{T},
    dual_names::DualNames,
) where {T}
    func_primal = MOI.get(primal_model, MOI.ConstraintFunction(), ci)
    primal_vi = func_primal
    func_dual = MOI.ScalarAffineFunction(
        MOIU.operate_terms(-, scalar_affine_terms[primal_vi]),
        sense_change * get(scalar_terms, primal_vi, zero(T)),
    )
    set_primal = MOI.get(primal_model, MOI.ConstraintSet(), ci)
    set_dual = _dual_set(set_primal)
    ci_map[ci] =
        MOIU.normalize_and_add_constraint(dual_model, func_dual, set_dual)
    if !is_empty(dual_names)
        set_dual_constraint_name(
            dual_model,
            primal_model,
            primal_vi,
            ci_map[ci],
            dual_names.dual_constraint_name_prefix,
        )
    end
    return
end

function add_scalar_affine_terms_from_quad_obj(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_var_dual_quad_slack::Dict{VI,VI},
    primal_objective::PrimalObjective{T},
) where {T}
    for term in primal_objective.obj.quadratic_terms
        if term.variable_1 == term.variable_2
            dual_vi = primal_var_dual_quad_slack[term.variable_1]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_1],
                -MOI.coefficient(term),
                dual_vi,
            )
        else
            dual_vi_1 = primal_var_dual_quad_slack[term.variable_1]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_2],
                -MOI.coefficient(term),
                dual_vi_1,
            )
            dual_vi_2 = primal_var_dual_quad_slack[term.variable_2]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[term.variable_1],
                -MOI.coefficient(term),
                dual_vi_2,
            )
        end
    end
    return
end

function add_scalar_affine_terms_from_quad_params(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_parameter::Dict{VI,VI},
    primal_objective::PrimalObjective{T},
) where {T}
    for (key, val) in primal_objective.quad_cross_parameters
        for term in val
            dual_vi = primal_parameter[term.variable]
            push_to_scalar_affine_terms!(
                scalar_affine_terms[key],
                -MOI.coefficient(term),
                dual_vi,
            )
        end
    end
end

function set_dual_constraint_name(
    dual_model::MOI.ModelLike,
    primal_model::MOI.ModelLike,
    primal_vi::VI,
    dual_ci::CI,
    prefix::String,
)
    MOI.set(
        dual_model,
        MOI.ConstraintName(),
        dual_ci,
        prefix * MOI.get(primal_model, MOI.VariableName(), primal_vi),
    )
    return
end

function get_scalar_terms(primal_objective::PrimalObjective{T}) where {T}
    scalar_terms = Dict{VI,T}()
    for term in get_affine_terms(primal_objective)
        if haskey(scalar_terms, term.variable)
            scalar_terms[term.variable] += MOI.coefficient(term)
        else
            scalar_terms[term.variable] = MOI.coefficient(term)
        end
    end
    return scalar_terms
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ::Type{F},
    ::Type{S},
) where {T,F,S}
    for ci in MOI.get(primal_model, MOI.ListOfConstraintIndices{F,S}())
        fill_scalar_affine_terms!(
            scalar_affine_terms,
            primal_con_dual_var,
            primal_model,
            ci,
        )
    end
    return
end

function get_scalar_affine_terms(
    primal_model::MOI.ModelLike,
    primal_con_dual_var::Dict{CI,Vector{VI}},
    variables::Vector{VI},
    con_types::Vector{Tuple{Type,Type}},
    ::Type{T},
) where {T}
    scalar_affine_terms = Dict{VI,Vector{MOI.ScalarAffineTerm{T}}}(
        vi => MOI.ScalarAffineTerm{T}[] for vi in variables
    )
    for (F, S) in con_types
        fill_scalar_affine_terms!(
            scalar_affine_terms,
            primal_con_dual_var,
            primal_model,
            F,
            S,
        )
    end
    return scalar_affine_terms
end

function push_to_scalar_affine_terms!(
    scalar_affine_terms::Vector{MOI.ScalarAffineTerm{T}},
    affine_term::T,
    vi::VI,
) where {T}
    if !iszero(affine_term) # if term is different than 0 add to the scalar affine terms vector
        push!(scalar_affine_terms, MOI.ScalarAffineTerm(affine_term, vi))
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{SAF{T},S},
) where {T,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T}}}
    moi_function = get_function(primal_model, ci)
    for term in moi_function.terms
        dual_vi = primal_con_dual_var[ci][1] # In this case we only have one vi
        push_to_scalar_affine_terms!(
            scalar_affine_terms[term.variable],
            MOI.coefficient(term),
            dual_vi,
        )
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{VI,S},
) where {T,S<:Union{MOI.GreaterThan{T},MOI.LessThan{T},MOI.EqualTo{T}}}
    dual_var = get(primal_con_dual_var, ci, nothing)
    if dual_var === nothing
        # No variables created as the primal constraint is the constraint
        # of a constrained variable. Hence, its duality information goes to
        # the dual constraint associated to that primal variable.
        return
    end
    moi_function = get_function(primal_model, ci)
    dual_vi = dual_var[1] # In this case we only have one vi
    push_to_scalar_affine_terms!(
        scalar_affine_terms[moi_function],
        one(T),
        dual_vi,
    )
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{VAF{T},S},
) where {T,S<:MOI.AbstractVectorSet}
    moi_function = get_function(primal_model, ci)
    set = get_set(primal_model, ci)
    for term in moi_function.terms
        dual_vi = primal_con_dual_var[ci][term.output_index]
        # term.output_index is the row of the VAF,
        # it corresponds to the dual variable associated with this constraint
        push_to_scalar_affine_terms!(
            scalar_affine_terms[term.scalar_term.variable],
            set_dot(term.output_index, set, T) * MOI.coefficient(term),
            dual_vi,
        )
    end
    return
end

function fill_scalar_affine_terms!(
    scalar_affine_terms::Dict{VI,Vector{MOI.ScalarAffineTerm{T}}},
    primal_con_dual_var::Dict{CI,Vector{VI}},
    primal_model::MOI.ModelLike,
    ci::CI{VVF,S},
) where {T,S<:MOI.AbstractVectorSet}
    dual_vars = get(primal_con_dual_var, ci, nothing)
    if dual_vars === nothing
        # No variables created as the primal constraint is the constraint
        # of a constrained variable. Hence, its duality information goes to
        # the dual constraint associated to that primal variable.
        return
    end
    moi_function = get_function(primal_model, ci)
    set = get_set(primal_model, ci)
    for (i, variable) in enumerate(moi_function.variables)
        dual_vi = dual_vars[i]
        push_to_scalar_affine_terms!(
            scalar_affine_terms[variable],
            set_dot(i, set, T) * one(T),
            dual_vi,
        )
    end
    return
end

struct CanonicalVector{T} <: AbstractVector{T}
    index::Int
    n::Int
end
Base.eltype(::Type{CanonicalVector{T}}) where {T} = T
Base.length(v::CanonicalVector) = v.n
Base.size(v::CanonicalVector) = (v.n,)
function Base.getindex(v::CanonicalVector{T}, i::Integer) where {T}
    return convert(T, i == v.index)
end

# This is much faster than the default implementation that goes
# through all entries even if only one is nonzero.
function LinearAlgebra.dot(
    x::CanonicalVector{T},
    y::CanonicalVector{T},
) where {T}
    return convert(T, x.index == y.index)
end
function MOI.Utilities.triangle_dot(
    x::CanonicalVector{T},
    y::CanonicalVector{T},
    dim::Int,
    offset::Int,
) where {T}
    if x.index != y.index || x.index <= offset
        return zero(T)
    elseif MOI.Utilities.is_diagonal_vectorized_index(x.index - offset)
        return one(T)
    else
        return 2one(T)
    end
end

function set_dot(i::Integer, s::MOI.AbstractVectorSet, T::Type)
    vec = CanonicalVector{T}(i, MOI.dimension(s))
    return MOIU.set_dot(vec, vec, s)
end
function set_dot(::Integer, ::MOI.AbstractScalarSet, T::Type)
    return one(T)
end


text:
push!(LOAD_PATH, "../src/")
using Documenter, Neuroimaging

makedocs(
    modules = [Neuroimaging],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "Neuroimaging.jl",
    authors = "Robert Luke",
    pages = [
        "Home" => "index.md",
        "Types" => "types.md",
        "Steady State Responses" => Any[
            "Overview"=>"assr/assr.md",
            "Example"=>"assr/examples.md",
            "API"=>"assr/functions.md",
        ],
        "API" => "api.md",
    ],
)

deploydocs(repo = "github.com/rob-luke/EEG.jl.git", push_preview = true, devbranch = "main")


text:
using Flux

mutable struct robust_rnn_stable
    nu
    nx
    nv
    ny
    ϕ
    E
    F
    B1
    B2
    C1
    D12
    C2
    D21
    D22
    λ
    P
    bx
    bv
    by
end

function robust_rnn_stable(nu, nx, nv, ny; nl=Flux.relu, T=Float64)
    T = Float64

    B1 = randn(nx, nv) / sqrt(nx + nv)
    B2 = randn(nx, nu) / sqrt(2 * (nx + nu))

    C1 = randn(nv, nx) / sqrt(nv + nx)
    D12 = zeros(nv, nu)

    λ = ones(nv, 1) # Choose sum abs values to make diagonally dominant
    Λ = 2diagm(λ[:, 1])

    C2 = zeros(ny, nx) / sqrt(ny)
    D21 = zeros(ny, nv) / sqrt(ny)
    D22 = zeros(ny, nu) / sqrt(ny)
    V = randn(2nx, 2nx) / sqrt(2nx)

    H = V' * V + [-C1'; B1] * [-C1 B1'] / 2
    H = (H + H') / 2
    S = randn(nx, nx) / sqrt(nx)

    P = H[nx + 1:end, nx + 1:end]
    E = (H[1:nx, 1:nx] + P) / 2 + S - S'
    F = H[nx + 1:end, 1:nx]

    # Check definiteness of LMI
    # lmi = [(E + E' - P) -C1' F';
    #         -C1 2Λ B1';
    #         F B1 P]

    # eigvals(lmi)

    bx = zeros(nx)
    bv = randn(nv) / sqrt(nv)
    by = zeros(ny)

    return robust_rnn_stable(nu, nx, nv, ny, nl, E, F, B1, B2, C1, D12, C2, D21, D22, λ, P, bx, bv, by)

end

Flux.trainable(M::robust_rnn_stable) = [M.E, M.F, M.B1, M.B2, M.C1, M.D12, M.C2,
                                    M.D21, M.D22, M.λ, M.P, M.bx, M.bv, M.by]

init_state(model::robust_rnn_stable, batches) = zeros(model.nx, batches)

function (M::robust_rnn_stable)(x0, u)

    E = M.E
    A = E \ M.F
    B1 = E \ M.B1
    B2 = E \ M.B2

    # Λ = Diagonal(M.λ[:, 1])
    Λ = I
    C1 = Λ \ M.C1
    D12 = Λ \ M.D12

    C2 = M.C2
    D21 = M.D21
    D22 = M.D22

    function cell(xt, ut)
        wt = M.ϕ.(C1 * xt + D12 * ut)
        xn = A * xt + B1 * wt + B2 * ut
        yt = C2 * xt + D21 * wt + D22 * ut
        return xn, (yt, xt, wt)
    end
    recurrent_cell  = Flux.Recur(cell, x0)

    return collect(unzip(recurrent_cell.(u)))
end

function eval_LMIs(model::robust_rnn_stable; ϵ=1E-6)

    E = model.E

    # return [ (E + E' - 2ϵ * I) ]

    F = model.F
    B1 = model.B1
    C1 = model.C1
    λ = model.λ
    P = model.P

    Λ = Diagonal(λ[:, 1])

    # Apprently assembly mat all at once mutates some memory
    row1 = [(E + E' - P) -C1' F']
    row2 = [-C1 2Λ B1']
    row3 = [F B1 P]
    mat = [row1; row2; row3]

    return [mat + mat' - 2ϵ * I] / 2
end

function eval_inequalities(model::robust_rnn_stable; ϵ=1E-6)
    return model.λ
end

# ------------------------------------------------------------------------------------------------

mutable struct robust_rnn_lipschitz
    nu
    nx
    nv
    ny
    γ
    ϕ
    E
    F
    B1
    B2
    C1
    D12
    C2
    D21
    D22
    λ
    P
    bx
    bv
    by
end

init_state(model::robust_rnn_lipschitz, batches) = zeros(model.nx, batches)

function robust_rnn_lipschitz(nu, nx, nv, ny, γ; nl=tanh)

    F = randn(nx, nx) / sqrt(nx)
    B1 = randn(nx, nv) / sqrt(nx)
    B2 = randn(nx, nu) / sqrt(nx)

    C1 = randn(nv, nx) / sqrt(nv)
    D12 = zeros(nv, nu)

    C2 = randn(ny, nx) / sqrt(ny)
    D21 = randn(ny, nv) / sqrt(ny)
    D22 = randn(ny, nu) / sqrt(ny)

    # Rescale input and output to ensure feasability
    Γy = sum(abs.([C2 D21 D22]), dims=2)
    sfy = γ ./ Γy
    C2 = C2 .* sfy
    D21 = D21 .* sfy
    D22 = D22 .* sfy

    Γu = sum(abs.([D12' B2' D22']), dims=2)
    sfu = γ ./ Γu
    D12 = D12 .* sfu'
    B2 = B2 .* sfu'
    D22 = D22 .* sfu'

    # Ensures feasability via diagonal dominance
    p = sum(abs.([F B1 B2]), dims=2)
    P = diagm(p[:, 1])
    λ = sum(abs.([C1 D12 B1' D21']), dims=2) / 2
    Λ = diagm(λ[:, 1])
    e = p + sum(abs.([C1' F' C2']), dims=2)
    E = diagm(e[:, 1]) / 2

    bx = zeros(nx)
    bv = randn(nv) / sqrt(nv)
    by = zeros(ny)

    return robust_rnn_lipschitz(nu, nx, nv, ny, γ,  nl, E, F, B1, B2, C1, D12, C2, D21, D22, λ, P, bx, bv, by)

end

Flux.trainable(M::robust_rnn_lipschitz) = [M.E, M.F, M.B1, M.B2, M.C1, M.D12, M.C2,
                                    M.D21, M.D22, M.λ, M.P, M.bx, M.bv, M.by]

function (M::robust_rnn_lipschitz)(x0, u)
    E = M.E
    A = E \ M.F
    B1 = E \ M.B1
    B2 = E \ M.B2

    Λ = Diagonal(M.λ[:, 1])
    C1 = Λ \ M.C1
    D12 = Λ \ M.D12

    C2 = M.C2
    D21 = M.D21
    D22 = M.D22

    function cell(xt, ut)
        wt = M.ϕ.(C1 * xt + D12 * ut)
        xn = A * xt + B1 * wt + B2 * ut
        yt = C2 * xt + D21 * wt + D22 * ut
        return xn, (yt, xt, wt)
    end
    recurrent_cell  = Flux.Recur(cell, x0)

    return collect(unzip(recurrent_cell.(u)))
end

function eval_inequalities(model::robust_rnn_lipschitz; ϵ=1E-6)
    return model.λ
end

function eval_LMIs(model::robust_rnn_lipschitz; ϵ=1E-6)
    E = model.E

    return [E + E' - 2ϵ * I] /  2

    F = model.F
    B1 = model.B1
    B2 = model.B2
    C1 = model.C1
    D12 = model.D12
    C2 = model.C2
    D21 = model.D21
    D22 = model.D22
    λ = model.λ
    P = model.P

    Λ = Diagonal(λ[:, 1])
    γ = model.γ

    Iu = Matrix(I, model.nu, model.nu)
    Iy = Matrix(I, model.ny, model.ny)

    row1 = [(E + E' - P) -C1' zeros(model.nx, model.nu) F' C2']
    row2 = [-C1 2Λ -D12 B1' D21']
    row3 = [zeros(model.nu, model.nx) -D12' (γ * Iu) B2' D22']
    row4 = [F B1 B2 P zeros(model.nx, model.ny)]
    row5 = [C2 D21 D22 zeros(model.ny, model.nx) (γ * Iy)]
    mat = [row1; row2; row3; row4; row5]

    return [mat + mat' - 2ϵ * I] /  2
end

max_stars_repo_name: imanchester/REN

max_issues_repo_name: imanchester/REN

max_forks_repo_name: imanchester/REN


text:
"""
    fractionaluncertainty(::Type{MonatomicGas}, z::Integer, energy)

Determines from the element and energy, the approximate range of fractional uncertainties to
associate with the total and photoelectric components of the mass attenuation coefficients
for monatomic gas samples.
Based on [this table](https://physics.nist.gov/PhysRefData/FFast/Text2000/sec06.html#tab2).
"""
function fractionaluncertainty(::Type{MonatomicGas}, z::Integer, energy)
    low, high = 0.01, 0.01
    if energy < 200.0
        low, high = 0.5, 1.0
    elseif energy < 500.0
        low, high = 0.20, 0.30
    elseif energy < 1.0
        low, high = 0.03, 0.10
    end
    distance = ( (energy - edgeenergy(z, sh)) / energy for sh in eachedge(z) )
    if minimum(abs.(distance)) < 0.001
        low, high = max(low, 0.2), max(high, 0.3)
    end
    u = [ energy / edgeenergy(z, sh) for sh in eachedge(z) ]
    if (u[1] > 1.0) && (u[1] < 1.1)
        low, high = max(low, 0.1), max(high, 0.1)
    elseif (u[1] >= 1.1) && (u[1] < 1.2)
        low, high = max(low, 0.03), max(high, 0.03)
    end
    # L1, M1, M2, M3
    for sh in filter(sh->get(u, sh, 0.0) >= 1.0, [ 2, 5, 6, 7 ])
        if u[sh] < 1.15
            low, high = max(low, 0.15), max(high, 0.15)
        elseif u[sh] < 1.4
            low, high = max(low, 0.04), max(high, 0.04)
        end
    end
    # L2, L3, M4, M5
    for sh in filter(sh->get(u, sh, 0.0) >= 1.0, [ 3, 4, 8, 9 ])
        if u[sh] < 1.15
            low, high = max(low, 0.20), max(high, 0.20)
        elseif u[sh] < 1.4
            low, high = max(low, 0.04), max(high, 0.04)
        end
    end
    if energy > 200.0e3
        low, high = max(low, 0.02), max(high, 0.03)
    end
    return (low, high)
end

"""
    fractionaluncertainty(::Type{SolidLiquid}, z::Integer, energy)

Determines from the element and energy, the approximate range of fractional uncertainties to
associate with the total and photoelectric components of the mass attenuation coefficients
for solids and liquids.
Based on [this table](https://physics.nist.gov/PhysRefData/FFast/Text2000/sec06.html#tab2).
"""
function fractionaluncertainty(::Type{SolidLiquid}, z::Integer, energy)
    low, high = 0.01, 0.01
    if energy < 200.0
        low, high = 1.0, 2.0
    elseif energy < 500.0
        low, high = 0.50, 1.0
    elseif energy < 1.0
        low, high = 0.05, 0.20
    end
    distance = ( (energy - edgeenergy(z, sh)) / energy for sh in eachedge(z) )
    if minimum(abs.(distance)) < 0.001
        low, high = max(low, 0.5), max(high, 0.5)
    end
    u = [ energy / edgeenergy(z, sh) for sh in eachedge(z) ]
    if (u[1] > 1.0) && (u[1] < 1.1)
        low, high = max(low, 0.1), max(high, 0.2)
    elseif (u[1] >= 1.1) && (u[1] < 1.2)
        low, high = max(low, 0.03), max(high, 0.03)
    end
    # L1, M1, M2, M3
    for sh in filter(sh->get(u, sh, 0.0) >= 1.0, [ 2, 5, 6, 7 ])
        if u[sh] < 1.15
            low, high = max(low, 0.15), max(high, 0.30)
        elseif u[sh] < 1.4
            low, high = max(low, 0.04), max(high, 0.04)
        end
    end
    # L2, L3, M4, M5
    for sh in filter(sh->get(u, sh, 0.0) >= 1.0, [ 3, 4, 8, 9 ] )
        if u[sh] < 1.15
            low, high = max(low, 0.20), max(high, 0.40)
        elseif u[sh] < 1.4
            low, high = max(low, 0.04), max(high, 0.04)
        end
    end
    if energy > 200.0e3
        low, high = max(low, 0.02), max(high, 0.03)
    end
    return (low, high)
end


text:
module SolidState

using ColorTypes
using Images
using Interpolations: LinearInterpolation
using LinearAlgebra
using Printf
using StaticArrays

include("types.jl")
include("transformations.jl")
include("cameras.jl")
include("materials.jl")
include("shapes.jl")
include("render.jl")

end # module


text:
function minimal_EGM(p::Param)
    nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
    yvec          = sqrt(2.0) * p.sigma .* nodes
    ywgt          = weights .* pi^(-0.5)
    avec          = collect(range(0.0,p.a_high,length = p.na))
    m             = Vector{Float64}[Float64[] for i in 1:p.nT]   # endogenous grid
    c             = Vector{Float64}[Float64[] for i in 1:p.nT]   # consumption function on m
    m[p.nT]       = [0.0,p.a_high]
    c[p.nT]       = [0.0,p.a_high]

    cg = cgrad(:viridis)
    cols = cg[range(0.0,stop=1.0,length=p.nT)]

    pl = plot(m[p.nT],c[p.nT],label="$(p.nT)",leg=:topright,title="Consumption Function",
              xlims = (0,p.a_high),ylims = (0,p.a_high), color = cols[p.nT],
              xlab = "Cash on Hand", ylab = "Consumption")
    # cycle back in time
    for it in p.nT-1:-1:1
        w1 = 0.0 .+ exp.(yvec) .+ p.R.*avec'   # w1 = y + yshock*R*savings:  next period wealth at all states. (p.ny,p.na)
        # get next period consumption on that wealth w1
        # interpolate on next period's endogenous grid m[it+1].
        # notice that the `interpolate` object needs to be able to extrapolate
        c1 = reshape(extrapolate(interpolate((m[it+1],),c[it+1],Gridded(Linear())),Line())(w1[:]) ,p.ny,p.na)
        c1[c1.<0] .= p.cfloor     # don't allow negative consumption
        Emu   = ywgt' * (1 ./ c1)   # Expected marginal utility (with log utility!). (p.na,1)
        rhs   = p.beta * p.R * Emu[:]   # RHS of euler equation
        c[it] = 1.0 ./ rhs   # inverse marginal utility function gives current consumption
        m[it] = avec .+ c[it]   # current period endogenous cash on hand grid. (p.na+1,1)

        # add credit constraint region
        # the lowest asset grid point on tomorrow's cash on hand (agrid[1]) corresponds to
        # saving exactly zero (or location on the lower bound of the asset grid)
        # well, we allow that people save zero (and consume all they have currently on hand).
        c[it] = vcat(0.0, c[it])   # prepend with 0
        m[it] = vcat(0.0, m[it])   #

        plot!(pl,m[it],c[it],label= it == 1 ? "$it" : "", color = cols[it])
    end
    pl = lens!(pl, [0, 2], [0, 2], inset = (1, bbox(0.2, 0.1, 0.25, 0.25)))
    pl
end

# artifial example: lowest m state corresponds to zero consumption,
# even if that means negative cash on hand
# simple EGM that allows borrowing also in final period of life
# just shifts everything left
function minimal_EGM_borrow(p::Param)
    nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
    yvec          = sqrt(2.0) * p.sigma .* nodes
    ywgt          = weights .* pi^(-0.5)
    avec          = collect(range(p.a_low,p.a_high,length = p.na))
    m             = Vector{Float64}[Float64[] for i in 1:p.nT]   # endogenous grid
    c             = Vector{Float64}[Float64[] for i in 1:p.nT]   # consumption function on m
    m[p.nT]       = [p.a_low,p.a_high]
    c[p.nT]       = [0.0,p.a_high]

    cg = cgrad(:viridis)
    cols = cg[range(0.0,stop=1.0,length=p.nT)]

    pl = plot(m[p.nT],c[p.nT],label="$(p.nT)",leg=:topright,title="Consumption Function",
              xlims = (p.a_low,p.a_high),ylims = (0.0,p.a_high), color = cols[p.nT],
              xlab = "Cash on Hand", ylab = "Consumption")
    # cycle back in time
    for it in p.nT-1:-1:1
        w1 = 0.0 .+ exp.(yvec) .+ p.R.*avec'   # w1 = y + yshock*R*savings:  next period wealth at all states. (p.ny,p.na)
        # get next period consumption on that wealth w1
        # interpolate on next period's endogenous grid m[it+1].
        # notice that the `interpolate` object needs to be able to extrapolate
        c1 = reshape(extrapolate(interpolate((m[it+1],),c[it+1],Gridded(Linear())),Line())(w1[:]) ,p.ny,p.na)
        c1[c1.<0] .= p.cfloor     # don't allow negative consumption
        Emu   = ywgt' * (1 ./ c1)   # Expected marginal utility (with log utility!). (p.na,1)
        rhs   = p.beta * p.R * Emu[:]   # RHS of euler equation
        c[it] = 1.0 ./ rhs   # inverse marginal utility function gives current consumption
        m[it] = avec .+ c[it]   # current period endogenous cash on hand grid. (p.na+1,1)

        # add credit constraint region
        # the lowest asset grid point on tomorrow's cash on hand (agrid[1]) corresponds to
        # saving exactly zero (or location on the lower bound of the asset grid)
        # well, we allow that people save zero (and consume all they have currently on hand).
        c[it] = vcat(0.0, c[it])   # prepend with 0
        m[it] = vcat(p.a_low, m[it])   #

        plot!(pl,m[it],c[it],label= it == 1 ? "$it" : "", color = cols[it])
    end
    pl = lens!(pl, [p.a_low, 2], [0.0, 2], inset = (1, bbox(0.2, 0.1, 0.25, 0.25)))
    m,c,pl
end

# artifial example: lowest m state corresponds to zero consumption,
# even if that means negative cash on hand
# simple EGM that allows borrowing also in final period of life
# just shifts everything left
# function minimal_EGM_borrow(p::Param)
#     nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
#     yvec          = sqrt(2.0) * p.sigma .* nodes
#     ywgt          = weights .* pi^(-0.5)
#     avec          = collect(range(p.a_low,p.a_high,length = p.na))
#     m             = Vector{Float64}[Float64[] for i in 1:p.nT]   # endogenous grid
#     c             = Vector{Float64}[Float64[] for i in 1:p.nT]   # consumption function on m
#     m[p.nT]       = [0.0,p.a_high]
#     c[p.nT]       = [0.0,p.a_high]

#     cg = cgrad(:viridis)
#     cols = cg[range(0.0,stop=1.0,length=p.nT)]

#     pl = plot(m[p.nT],c[p.nT],label="$(p.nT)",leg=:topright,title="Consumption Function",
#               xlims = (p.a_low,p.a_high),ylims = (0.0,p.a_high), color = cols[p.nT],
#               xlab = "Cash on Hand", ylab = "Consumption")
#     # cycle back in time
#     for it in p.nT-1:-1:1
#         w1 = 0.0 .+ exp.(yvec) .+ p.R.*avec'   # w1 = y + yshock*R*savings:  next period wealth at all states. (p.ny,p.na)
#         # get next period consumption on that wealth w1
#         # interpolate on next period's endogenous grid m[it+1].
#         # notice that the `interpolate` object needs to be able to extrapolate
#         c1 = reshape(extrapolate(interpolate((m[it+1],),c[it+1],Gridded(Linear())),Line())(w1[:]) ,p.ny,p.na)
#         c1[c1.<0] .= p.cfloor     # don't allow negative consumption
#         Emu   = ywgt' * (1 ./ c1)   # Expected marginal utility (with log utility!). (p.na,1)
#         rhs   = p.beta * p.R * Emu[:]   # RHS of euler equation
#         c[it] = 1.0 ./ rhs   # inverse marginal utility function gives current consumption
#         m[it] = avec .+ c[it]   # current period endogenous cash on hand grid. (p.na+1,1)

#         # add credit constraint region
#         # the lowest asset grid point on tomorrow's cash on hand (agrid[1]) corresponds to
#         # saving exactly zero (or location on the lower bound of the asset grid)
#         # well, we allow that people save zero (and consume all they have currently on hand).
#         c[it] = vcat(0.0, c[it])   # prepend with 0
#         m[it] = vcat(p.a_low, m[it])   #

#         plot!(pl,m[it],c[it],label= it == 1 ? "$it" : "", color = cols[it])
#     end
#     pl = lens!(pl, [p.a_low, 2], [0.0, 2], inset = (1, bbox(0.2, 0.1, 0.25, 0.25)))
#     m,c,pl
# end

function minimal_EGM_bequest(p::Param)
    nodes,weights = gausshermite(p.ny)  # from FastGaussQuadrature
    yvec          = sqrt(2.0) * p.sigma .* nodes
    ywgt          = weights .* pi^(-0.5)
    # avec          = collect(range(p.a_low,stop = p.a_high,length = p.na))
    avec          = scaleGrid(0.0,p.a_high,p.na,logorder = 1)
    m             = Vector{Float64}[Float64[] for i in 1:p.nT]   # endogenous grid
    c             = Vector{Float64}[Float64[] for i in 1:p.nT]   # consumption function on m
    v             = Vector{Float64}[Float64[] for i in 1:p.nT]   # value function on m

    m[p.nT]       = avec    # no debt in last period possible
    c[p.nT]       = avec    # consume all you have left
    v[p.nT]       = p.ν > 0 ? bequest(avec,p) : zeros(p.na)

    # plot setup
    cg = cgrad(:plasma)
    cols = cg[range(0.0,stop=1.0,length=p.nT)]
    pl = plot(m[p.nT],c[p.nT],label="$(p.nT)",leg=false,title="Consumption",color = cols[p.nT], xlims = (0,p.a_high),ylims = (0,p.a_high))
    pv = plot(m[p.nT],v[p.nT],label="$(p.nT)",leg=:bottomright,title="Values",color = cols[p.nT], xlims = (0,p.a_high), ylims = (-20,2))

    for it in (p.nT-1):-1:1
        w1 = 0.0 .+ exp.(yvec) .+ p.R.*avec'   # w1 = y + yshock*R*savings:  next period wealth at all states. (p.ny,p.na)
        # get next period consumption on that wealth w1
        # interpolate on next period's endogenous grid m[it+1].
        # notice that the `interpolate` object needs to be able to extrapolate
        c1 = reshape(extrapolate(interpolate((m[it+1],),c[it+1],Gridded(Linear())),Line())(w1[:]) ,p.ny,p.na)
        c1[c1.<0] .= p.cfloor     # don't allow negative consumption

        # take care of bequest: if this is the penultimate period,
        # need to decide how much to leave for bequest
        if it == p.nT-1
            if p.ν > 0
                # if there is a bequest motive at all
                # this is like any other savings decision, except for that
                # next (i.e. final) period utility is the bequest function.
                # this looks like u(b + bbar) where bbar > 0 governs says that people can leave
                # zero bequests and survive (i.e. do not get u(0) = -Inf).
                mu1 = ywgt' * (1 ./ (c1 .+ p.bbar))  # expected marginal utility of consumption next period
                rhs = p.beta * p.R * p.ν * mu1[:]   # RHS of Euler equation
                c[it] = vcat(0.0, (1.0 ./ rhs[:])...)   # current period consumption vector. (p.na+1,1)
            else
                # if there is no bequest motive, this is like a standard period
                mu1 = ywgt' * (1 ./ c1)   # rhs of euler equation (with log utility!). (p.na,1)
                rhs = p.beta * p.R * mu1[:]
                c[it] = vcat(0.0, (1.0 ./ rhs[:])...)   # current period consumption vector. (p.na+1,1)
            end
        else
            mu1 = ywgt' * (1 ./ c1)   # rhs of euler equation (with log utility!). (p.na,1)
            rhs = p.beta * p.R * mu1[:]
            c[it] = vcat(0.0, (1.0 ./ rhs[:])...)   # current period consumption vector. (p.na+1,1)
        end

        # we prepend c and m with (0) here to create the borrowing constrained region
        m[it] = vcat(0.0, avec .+ c[it][2:end]...)   # current period endogenous cash on hand grid. (p.na+1,1)
        plot!(pl,m[it],c[it],label = it == 1 ? "$it" : "",color = cols[it])

        # compute value function
        # next period values at all states w1
        v1 = reshape(extrapolate(interpolate((m[it+1],),v[it+1],Gridded(Linear())),Line())(w1[:]) ,p.ny,p.na)
        ev = ywgt' * v1  # integrate out wage shocks
        v[it] = vcat(ev[1], u(c[it][2:end],p) .+ p.beta * ev[:])  # ev[1] is expected value of saving exactly zero.
        plot!(pv, m[it][2:end], v[it][2:end],label = it == 1 ? "$it" : "",color = cols[it])

    end

    pl = lens!(pl, [0, 2], [0, 2], inset = (1, bbox(0.2, 0.05, 0.3, 0.25)))
    return plot(pl,pv,layout = (1,2))
end


text:
module MatchV4 # LOLTools

using ..LOLTools: AbstractDTO, lol_api_server, http_action, call_api, nothing_in_event
using HTTP

const ParticipantIdentityDTO = NamedTuple{(:participantId, :player), Tuple{Int, NamedTuple{(:platformId, :accountId, :summonerName, :summonerId, :currentPlatformId, :currentAccountId, :matchHistoryUri, :profileIcon)}}}
const TeamStatsDTO = NamedTuple{(:teamId, :win, :firstBlood, :firstTower, :firstInhibitor, :firstBaron, :firstDragon, :firstRiftHerald, :towerKills, :inhibitorKills, :baronKills, :dragonKills, :vilemawKills, :riftHeraldKills, :dominionVictoryScore, :bans)}
const ParticipantDTO = NamedTuple{(:participantId, :teamId, :championId, :spell1Id, :spell2Id, :stats, :timeline)}

struct MatchDTO <: AbstractDTO
    seasonId
    queueId
    gameId
    participantIdentities::Vector{ParticipantIdentityDTO}
    gameVersion
    platformId
    gameMode
    mapId
    gameType
    teams::Vector{TeamStatsDTO}
    participants::Vector{ParticipantDTO}
    gameDuration
    gameCreation
    MatchDTO(seasonId, queueId, gameId, participantIdentities, gameVersion, platformId, gameMode, mapId, gameType, teams, participants, gameDuration, gameCreation) = new(seasonId, queueId, gameId, participantIdentities, gameVersion, platformId, gameMode, mapId, gameType, teams, participants, gameDuration, gameCreation)
end

const MatchReferenceDTO = NamedTuple{(:platformId, :gameId, :champion, :queue, :season, :timestamp, :role, :lane)}

struct MatchlistDTO <: AbstractDTO
    matches::Vector{MatchReferenceDTO}
    totalGames
    startIndex
    endIndex
    MatchlistDTO(matches, totalGames, startIndex, endIndex) = new(matches, totalGames, startIndex, endIndex)
end

const MatchFrameDTO = NamedTuple{(:participantFrames, :events, :timestamp)}

struct MatchTimelineDTO <: AbstractDTO
    frames::Vector{MatchFrameDTO}
    frameInterval
    MatchTimelineDTO(frames, frameInterval) = new(frames, frameInterval)
end

"""
         match_by_tournament_code(api_key::String,
                                  platform::String,
                                  tournamentCode::String ;
                                  endpoint::HTTP.URI = lol_api_server(platform),
                                  action::Function = http_action,
                                  event::Function = nothing_in_event)::MatchDTO
"""
function match_by_tournament_code(api_key::String,
                                  platform::String,
                                  tournamentCode::String ;
                                  endpoint::HTTP.URI = lol_api_server(platform),
                                  action::Function = http_action,
                                  event::Function = nothing_in_event)::MatchDTO
    call_api(MatchDTO, api_key, action, endpoint, "/lol/match/v4/matches/by-tournament-code/$tournamentCode/ids", event, match_by_tournament_code)
end

"""
         match_by_id(api_key::String,
                     platform::String,
                     matchId::Int64 ;
                     endpoint::HTTP.URI = lol_api_server(platform),
                     action::Function = http_action,
                     event::Function = nothing_in_event)::MatchDTO
"""
function match_by_id(api_key::String,
                     platform::String,
                     matchId::Int64 ;
                     endpoint::HTTP.URI = lol_api_server(platform),
                     action::Function = http_action,
                     event::Function = nothing_in_event)::MatchDTO
    call_api(MatchDTO, api_key, action, endpoint, "/lol/match/v4/matches/$matchId", event, match_by_id)
end

"""
         matchlists(api_key::String,
                    platform::String,
                    encryptedAccountId::String ;
                    endpoint::HTTP.URI = lol_api_server(platform),
                    action::Function = http_action,
                    event::Function = nothing_in_event,
                    query::Union{Nothing,Dict{String,String}} = nothing)::MatchlistDTO
"""
function matchlists(api_key::String,
                    platform::String,
                    encryptedAccountId::String ;
                    endpoint::HTTP.URI = lol_api_server(platform),
                    action::Function = http_action,
                    event::Function = nothing_in_event,
                    query::Union{Nothing,Dict{String,String}} = nothing)::MatchlistDTO
    call_api(MatchlistDTO, api_key, action, endpoint, "/lol/match/v4/matchlists/by-account/$encryptedAccountId", event, matchlists; query=query)
end

"""
         timelines(api_key::String,
                   platform::String,
                   matchId::Int64 ;
                   endpoint::HTTP.URI = lol_api_server(platform),
                   action::Function = http_action,
                   event::Function = nothing_in_event)::MatchTimelineDTO
"""
function timelines(api_key::String,
                   platform::String,
                   matchId::Int64 ;
                   endpoint::HTTP.URI = lol_api_server(platform),
                   action::Function = http_action,
                   event::Function = nothing_in_event)::MatchTimelineDTO
    call_api(MatchTimelineDTO, api_key, action, endpoint, "/lol/match/v4/timelines/by-match/$matchId", event, timelines)
end

end # module LOLTools.MatchV4


text:
struct Layers
    layers::Vector{Layer}
end

Base.convert(::Type{Layers}, s::Layer) = Layers([s])
Base.convert(::Type{Layers}, s::Layers) = s

Base.getindex(v::Layers, i::Int) = v.layers[i]
Base.length(v::Layers) = length(v.layers)
Base.eltype(::Type{Layers}) = Layer
Base.iterate(v::Layers, args...) = iterate(v.layers, args...)

const OneOrMoreLayers = Union{Layers, Layer}

function Base.:+(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers(vcat(l1.layers, l2.layers))
end

function Base.:*(s1::OneOrMoreLayers, s2::OneOrMoreLayers)
    l1::Layers, l2::Layers = s1, s2
    return Layers([el1 * el2 for el1 in l1 for el2 in l2])
end

summary(v) = iscontinuous(v) ? extrema(v) : collect(uniquesorted(vec(v)))
mergesummaries(s1::AbstractVector, s2::AbstractVector) = mergesorted(s1, s2)
mergesummaries(s1::Tuple, s2::Tuple) = extend_extrema(s1, s2)

mergelabels(a, b) = a

function Entries(s::OneOrMoreLayers, palettes=NamedTuple())
    layers::Layers = s
    summaries = arguments()
    labels = arguments()
    entries = Entry[]
    for labeledentries in process_transformations(layers)
        for le in labeledentries
            entry = Entry(le.plottype, map(getvalue, le.mappings), le.attributes)
            push!(entries, entry)
            mergewith!(mergesummaries, summaries, map(summary, entry.mappings))
            mergewith!(mergelabels, labels, map(getlabel, le.mappings))
        end
    end
    palettes = merge!(default_palettes(), arguments(; palettes...))
    scales = default_scales(summaries, palettes)
    return Entries(entries, scales, labels)
end

function AbstractPlotting.plot!(fig, s::OneOrMoreLayers;
                                axis=NamedTuple(), palettes=NamedTuple())
    return plot!(fig, Entries(s, palettes); axis)
end

function AbstractPlotting.plot(s::OneOrMoreLayers;
                               axis=NamedTuple(), figure=NamedTuple(), palettes=NamedTuple())
    return plot(Entries(s, palettes); axis, figure)
end

# Convenience function, which may become superfluous if `plot` also calls `facet!`
function draw(s::OneOrMoreLayers;
              axis = NamedTuple(), figure=NamedTuple(), palettes=NamedTuple())
    fg = plot(s; axis, figure, palettes)
    facet!(fg)
    Legend(fg, Entries(s))
    resizetocontent!(fg)
    return fg
end

function draw!(fig, s::OneOrMoreLayers;
               axis=NamedTuple(), palettes=NamedTuple())
    ag = plot!(fig, s; axis, palettes)
    facet!(fig, ag)
    return ag
end


text:
# Integrations with NBodySimulator.jl

struct ElementMassBody{cType<:Real,mType<:Real} <: Body
    r::SVector{3,cType}
    v::SVector{3,cType}
    m::mType
    s::Symbol
end
function ElementMassBody(r::SVector{3,<:Unitful.Length}, v::SVector{3,<:Unitful.Velocity}, e::Element)
    ElementMassBody{Float64,Float64}(austrip.(r), austrip.(v), austrip(e.atomic_mass), Symbol(e.symbol))
end

function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, L::Unitful.Length, reference_temp::Unitful.Temperature; rng = MersenneTwister(n))
    e = elements[symbol]
    average_velocity = √(u"k" * reference_temp / e.atomic_mass)
    generate_bodies_in_cell_nodes(n, symbol, average_velocity, L, rng = rng)
end
function NBodySimulator.generate_bodies_in_cell_nodes(n::Integer, symbol::Union{Integer,AbstractString,Symbol}, average_velocity::Unitful.Velocity, L::Unitful.Length; rng = MersenneTwister(n))
    velocities = average_velocity * randn(rng, Float64, (3, n))
    e = elements[symbol]
    bodies = ElementMassBody[]

    count = 1
    dL = L / (ceil(n^(1 / 3)))
    for x ∈ dL/2:dL:L, y ∈ dL/2:dL:L, z ∈ dL/2:dL:L
        if count > n
            break
        end
        push!(bodies, ElementMassBody(SVector(x, y, z), SVector{3}(velocities[:, count]), e))
        count += 1
    end
    return bodies
end

function get_symbols(system::PotentialNBodySystem{ElementMassBody})
    [system.bodies[i].s for i ∈ 1:length(system.bodies)]
end

function bodies(system::AbstractSystem)
    [ElementMassBody(position(a), velocity(a), species(a)) for a ∈ system]
end

function nbody_boundary_conditions(system::AbstractSystem)
    # TODO: support more boundary conditions
    box = bounding_box(system)
    @assert box[1][1] == box[2][2] == box[3][3]
    @assert boundary_conditions(system) == [Periodic(), Periodic(), Periodic()]
    CubicPeriodicBoundaryConditions(austrip(box[1][1]))
end

# TODO: support more boundary conditions
function DynamicAtom(b::ElementMassBody, boundary_conditions::CubicPeriodicBoundaryConditions)
    DynamicAtom(b, boundary_conditions.L * u"bohr")
end
function DynamicAtom(b::ElementMassBody, box_size::Unitful.Length)
    DynamicAtom(mod.(b.r .* u"bohr", box_size), b.v .* u"bohr * hartree / ħ_au", elements[b.s])
end

# TODO: support more boundary conditions
function DynamicSystem(bodies::Vector{<:ElementMassBody}, boundary_conditions::CubicPeriodicBoundaryConditions, time::Real = 0.0)
    DynamicSystem(bodies, boundary_conditions.L * u"bohr", time * u"ħ_au / hartree")
end
function DynamicSystem(bodies::Vector{<:ElementMassBody}, box_size::Unitful.Length, time::Unitful.Time = 0.0u"s")
    DynamicSystem(DynamicAtom.(bodies, box_size), box_size, time)
end

@kwdef struct InteratomicPotentialParameters <: PotentialParameters
    potential::ArbitraryPotential
    timestep_cache::RefValue{Real} = Ref{Real}()
    force_cache::RefValue{Vector{SVector{3,Real}}} = Ref{Vector{SVector{3,Real}}}()
end

function NBodySimulator.get_accelerating_function(parameters::InteratomicPotentialParameters, simulation::NBodySimulation)
    masses = get_masses(simulation.system)
    symbols = get_symbols(simulation.system)
    (dv, u, v, t, i) -> begin
        if !isassigned(parameters.timestep_cache) || t != parameters.timestep_cache[]
            particles = [ElementMassBody(SVector{3}(u[:, j]), SVector{3}(v[:, j]), masses[j], symbols[j]) for j ∈ 1:length(symbols)]
            system = DynamicSystem(particles, simulation.boundary_conditions, t)
            parameters.timestep_cache[] = t
            parameters.force_cache[] = force(system, parameters.potential)
        end
        dv .+= parameters.force_cache[][i] / masses[i]
    end
end

@kwdef struct LJPotential
    ϵ::Real
    σ::Real
    R::Real
end
function LJPotential(ϵ::Unitful.Energy, σ::Unitful.Length, R::Unitful.Length)
    LJPotential(austrip(ϵ), austrip(σ), austrip(R))
end

"""
    NBSimulator <: MolecularDynamicsSimulator

A wrapper around NBodySimulator to implement the Atomistic API.

**Field descriptions**
- `Δt::Real` the time between timesteps, assumed to be in atomic units
- `steps::Integer` the number of timesteps for the simulation
- `t₀::Real` the starting time of the simulation, assumed to be in atomic units;
    defaults to 0
- `thermostat::Thermostat` the thermostat for the simulation;
    many options are defined by `NBodySimulator`, but a user could also define a custom thermostat;
    defaults to the `NullThermostat`
- `simulator::OrdinaryDiffEqAlgorithm` the algorithm to be used for the ODE;
    defaults to VelocityVerlet
- `potentials::Dict{Symbol,PotentialParameters}` dictionary of potentials;
    shouldn't be manipulated directly by the user
"""
@kwdef struct NBSimulator <: MolecularDynamicsSimulator
    Δt::Real
    steps::Integer
    t₀::Real = 0.0
    thermostat::Thermostat = NullThermostat()
    simulator::OrdinaryDiffEqAlgorithm = VelocityVerlet()
    potentials::Dict{Symbol,PotentialParameters} = Dict{Symbol,PotentialParameters}()
end
function NBSimulator(Δt::Unitful.Time,
    steps::Integer,
    t₀::Unitful.Time = 0.0u"s",
    thermostat::Thermostat = NullThermostat(),
    simulator::OrdinaryDiffEqAlgorithm = VelocityVerlet(),
    potentials::Dict{Symbol,PotentialParameters} = Dict{Symbol,PotentialParameters}())
    NBSimulator(austrip(Δt), steps, austrip(t₀), thermostat, simulator, potentials)
end

function simulate(system::AbstractSystem, simulator::NBSimulator, potential::ArbitraryPotential)
    simulator.potentials[:custom] = InteratomicPotentialParameters(potential = potential)
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator, potential::LJPotential)
    simulator.potentials[:lennard_jones] = LennardJonesParameters(potential.ϵ, potential.σ, potential.R)
    simulate(system, simulator)
end
function simulate(system::AbstractSystem, simulator::NBSimulator)
    nb_system = PotentialNBodySystem{ElementMassBody}(bodies(system), simulator.potentials)
    simulation = NBodySimulation(nb_system, (simulator.t₀, simulator.t₀ + simulator.steps * simulator.Δt), nbody_boundary_conditions(system), simulator.thermostat, 1.0)
    NBSResult(run_simulation(simulation, simulator.simulator, dt = simulator.Δt))
end

"""
    NBSResult <: MolecularDynamicsResult

The result generating from running a `MolecularDynamicsSimulator`.

**Field descriptions**
- `result::SimulationResult` the standard simulation result from `NBodySimulator`
"""
struct NBSResult <: MolecularDynamicsResult
    result::SimulationResult
end

function get_system(result::NBSResult, t::Integer = 0)
    sr = result.result
    time = get_time_range(result)[t > 0 ? t : end]
    positions = get_position(sr, time)
    velocities = get_velocity(sr, time)
    masses = get_masses(sr.simulation.system)
    symbols = get_symbols(sr.simulation.system)
    particles = [ElementMassBody(SVector{3}(positions[:, i]), SVector{3}(velocities[:, i]), masses[i], symbols[i]) for i ∈ 1:length(symbols)]
    DynamicSystem(particles, sr.simulation.boundary_conditions, time)
end

function get_time_range(result::NBSResult)
    result.result.solution.t
end

function temperature(result::NBSResult, t::Integer = 0)
    time = get_time_range(result)[t > 0 ? t : end]
    NBodySimulator.temperature(result.result, time)
end

function reference_temperature(result::NBSResult)
    result.result.simulation.thermostat isa NullThermostat ? missing : result.result.simulation.thermostat.T
end

function kinetic_energy(result::NBSResult, t::Integer = 0)
    time = get_time_range(result)[t > 0 ? t : end]
    NBodySimulator.kinetic_energy(result.result, time)
end

function potential_energy(result::NBSResult, t::Integer = 0)
    potentials = result.result.simulation.system.potentials
    # https://github.com/SciML/NBodySimulator.jl/issues/44
    if :custom ∈ keys(potentials)
        return InteratomicPotentials.potential_energy(get_system(result, t), potentials[:custom].potential)
    end
    time = get_time_range(result)[t > 0 ? t : end]
    NBodySimulator.potential_energy(result.result, time)
end

function rdf(result::NBSResult, sample_fraction::Float64 = 1.0)
    @assert 0 < sample_fraction ≤ 1
    sr = result.result
    n = length(sr.simulation.system.bodies)
    pbc = sr.simulation.boundary_conditions
    trange = get_time_range(result)[end-floor(Int, length(sr.solution.t) * sample_fraction)+1:end]

    maxbin = 1000
    dr = pbc.L / 2 / maxbin
    hist = zeros(maxbin)
    for t ∈ trange
        cc = get_position(sr, t)
        for i ∈ 1:n
            ri = SVector{3}(cc[:, i])
            for j ∈ i+1:n
                rj = SVector{3}(cc[:, j])
                (rij, r, r2) = NBodySimulator.get_interparticle_distance(ri, rj, pbc)
                if r2 < (0.5 * pbc.L)^2
                    bin = ceil(Int, r / dr)
                    if bin > 1 && bin <= maxbin
                        hist[bin] += 2
                    end
                end
            end
        end
    end

    c = 4 / 3 * π * n / pbc.L^3
    gr = zeros(maxbin)
    rs = zeros(maxbin)
    tlen = length(trange)
    for bin ∈ 1:maxbin
        rlower = (bin - 1) * dr
        rupper = rlower + dr
        nideal = c * (rupper^3 - rlower^3)
        gr[bin] = (hist[bin] / (tlen * n)) / nideal
        rs[bin] = rlower + dr / 2
    end

    rs, gr
end

max_stars_repo_name: jrdegreeff/cesmix-julia

max_issues_repo_name: jrdegreeff/cesmix-julia

max_forks_repo_name: jrdegreeff/cesmix-julia


text:
## number-theoretic functions ##

function gcd{T<:Integer}(a::T, b::T)
    while b != 0
        t = b
        b = rem(a, b)
        a = t
    end
    abs(a)
end

# explicit a==0 test is to handle case of lcm(0,0) correctly
lcm{T<:Integer}(a::T, b::T) = a == 0 ? a : abs(a * div(b, gcd(b,a)))

gcd(a::Integer) = a
lcm(a::Integer) = a
gcd(a::Integer, b::Integer) = gcd(promote(a,b)...)
lcm(a::Integer, b::Integer) = lcm(promote(a,b)...)
gcd(a::Integer, b::Integer...) = gcd(a, gcd(b...))
lcm(a::Integer, b::Integer...) = lcm(a, lcm(b...))

gcd{T<:Integer}(abc::AbstractArray{T}) = reduce(gcd,abc)
lcm{T<:Integer}(abc::AbstractArray{T}) = reduce(lcm,abc)

# return (gcd(a,b),x,y) such that ax+by == gcd(a,b)
function gcdx{T<:Integer}(a::T, b::T)
    s0, s1 = one(T), zero(T)
    t0, t1 = s1, s0
    while b != 0
        q = div(a, b)
        a, b = b, rem(a, b)
        s0, s1 = s1, s0 - q*s1
        t0, t1 = t1, t0 - q*t1
    end
    a < 0 ? (-a, -s0, -t0) : (a, s0, t0)
end
gcdx(a::Integer, b::Integer) = gcdx(promote(a,b)...)

# multiplicative inverse of n mod m, error if none
function invmod(n, m)
    g, x, y = gcdx(n, m)
    if g != 1 || m == 0
        error("no inverse exists")
    end
    x < 0 ? abs(m) + x : x
end

# ^ for any x supporting *
to_power_type(x::Number) = oftype(x*x, x)
to_power_type(x) = x
function power_by_squaring(x, p::Integer)
    x = to_power_type(x)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        throw(DomainError())
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end
power_by_squaring(x::Bool, p::Unsigned) = ((p==0) | x)
function power_by_squaring(x::Bool, p::Integer)
    p < 0 && throw(DomainError())
    return (p==0) | x
end

^{T<:Integer}(x::T, p::T) = power_by_squaring(x,p)
^(x::Number, p::Integer)  = power_by_squaring(x,p)
^(x, p::Integer)          = power_by_squaring(x,p)

# b^p mod m
function powermod{T}(b::Integer, p::Integer, m::T)
    p < 0 && throw(DomainError())
    b = oftype(m,mod(b,m))  # this also checks for divide by zero
    p == 0 && return mod(one(b),m)
    (m == 1 || m == -1) && return zero(m)

    t = prevpow2(p)
    local r::T
    r = 1
    while true
        if p >= t
            r = mod(widemul(r,b),m)
            p -= t
        end
        t >>>= 1
        t <= 0 && break
        r = mod(widemul(r,r),m)
    end
    return r
end

# smallest power of 2 >= x
nextpow2(x::Unsigned) = one(x)<<((sizeof(x)<<3)-leading_zeros(x-1))
nextpow2(x::Integer) = oftype(x,x < 0 ? -nextpow2(unsigned(-x)) : nextpow2(unsigned(x)))

prevpow2(x::Unsigned) = (one(x)>>(x==0)) << ((sizeof(x)<<3)-leading_zeros(x)-1)
prevpow2(x::Integer) = oftype(x,x < 0 ? -prevpow2(unsigned(-x)) : prevpow2(unsigned(x)))

ispow2(x::Integer) = count_ones(x)==1

# smallest a^n >= x, with integer n
function nextpow(a::Real, x::Real)
    (a <= 1 || x <= 0) && throw(DomainError())
    x <= 1 && return one(a)
    n = iceil(log(a, x))
    p = a^(n-1)
    # guard against roundoff error, e.g., with a=5 and x=125
    p >= x ? p : a^n
end
# largest a^n <= x, with integer n
function prevpow(a::Real, x::Real)
    (a <= 1 || x < 1) && throw(DomainError())
    n = ifloor(log(a, x))
    p = a^(n+1)
    p <= x ? p : a^n
end

# decimal digits in an unsigned integer
const powers_of_ten = [
    0x0000000000000001, 0x000000000000000a, 0x0000000000000064, 0x00000000000003e8,
    0x0000000000002710, 0x00000000000186a0, 0x00000000000f4240, 0x0000000000989680,
    0x0000000005f5e100, 0x000000003b9aca00, 0x00000002540be400, 0x000000174876e800,
    0x000000e8d4a51000, 0x000009184e72a000, 0x00005af3107a4000, 0x00038d7ea4c68000,
    0x002386f26fc10000, 0x016345785d8a0000, 0x0de0b6b3a7640000, 0x8ac7230489e80000,
]
function ndigits0z(x::Union(Uint8,Uint16,Uint32,Uint64))
    lz = (sizeof(x)<<3)-leading_zeros(x)
    nd = (1233*lz)>>12+1
    nd -= x < powers_of_ten[nd]
end
function ndigits0z(x::Uint128)
    n = 0
    while x > 0x8ac7230489e80000
        x = div(x,0x8ac7230489e80000)
        n += 19
    end
    return n + ndigits0z(uint64(x))
end
ndigits0z(x::Integer) = ndigits0z(unsigned(abs(x)))

const ndigits_max_mul = WORD_SIZE==32 ? 69000000 : 290000000000000000

function ndigits0z(n::Unsigned, b::Int)
    b == 2  && return (sizeof(n)<<3-leading_zeros(n))
    b == 8  && return div((sizeof(n)<<3)-leading_zeros(n)+2,3)
    b == 16 && return (sizeof(n)<<1)-(leading_zeros(n)>>2)
    b == 10 && return ndigits0z(n)
    d = 0
    while ndigits_max_mul < n
        n = div(n,b)
        d += 1
    end
    m = 1
    while m <= n
        m *= b
        d += 1
    end
    return d
end
ndigits0z(x::Integer, b::Integer) = ndigits0z(unsigned(abs(x)),int(b))

ndigits(x::Unsigned, b::Integer) = x==0 ? 1 : ndigits0z(x,int(b))
ndigits(x::Unsigned)             = x==0 ? 1 : ndigits0z(x)

ndigits(x::Integer, b::Integer) = ndigits(unsigned(abs(x)),int(b))
ndigits(x::Integer) = ndigits(unsigned(abs(x)))

## integer to string functions ##

function bin(x::Unsigned, pad::Int, neg::Bool)
    i = neg + max(pad,sizeof(x)<<3-leading_zeros(x))
    a = Array(Uint8,i)
    while i > neg
        a[i] = '0'+(x&0x1)
        x >>= 1
        i -= 1
    end
    if neg; a[1]='-'; end
    ASCIIString(a)
end

function oct(x::Unsigned, pad::Int, neg::Bool)
    i = neg + max(pad,div((sizeof(x)<<3)-leading_zeros(x)+2,3))
    a = Array(Uint8,i)
    while i > neg
        a[i] = '0'+(x&0x7)
        x >>= 3
        i -= 1
    end
    if neg; a[1]='-'; end
    ASCIIString(a)
end

function dec(x::Unsigned, pad::Int, neg::Bool)
    i = neg + max(pad,ndigits0z(x))
    a = Array(Uint8,i)
    while i > neg
        a[i] = '0'+rem(x,10)
        x = oftype(x,div(x,10))
        i -= 1
    end
    if neg; a[1]='-'; end
    ASCIIString(a)
end

function hex(x::Unsigned, pad::Int, neg::Bool)
    i = neg + max(pad,(sizeof(x)<<1)-(leading_zeros(x)>>2))
    a = Array(Uint8,i)
    while i > neg
        d = x & 0xf
        a[i] = '0'+d+39*(d>9)
        x >>= 4
        i -= 1
    end
    if neg; a[1]='-'; end
    ASCIIString(a)
end

num2hex(n::Integer) = hex(n, sizeof(n)*2)

const base36digits = ['0':'9','a':'z']
const base62digits = ['0':'9','A':'Z','a':'z']

function base(b::Int, x::Unsigned, pad::Int, neg::Bool)
    if !(2 <= b <= 62) error("invalid base: $b") end
    digits = b <= 36 ? base36digits : base62digits
    i = neg + max(pad,ndigits0z(x,b))
    a = Array(Uint8,i)
    while i > neg
        a[i] = digits[1+rem(x,b)]
        x = div(x,b)
        i -= 1
    end
    if neg; a[1]='-'; end
    ASCIIString(a)
end
base(b::Integer, n::Integer, pad::Integer=1) = base(int(b), unsigned(abs(n)), pad, n<0)

for sym in (:bin, :oct, :dec, :hex)
    @eval begin
        ($sym)(x::Unsigned, p::Int) = ($sym)(x,p,false)
        ($sym)(x::Unsigned)         = ($sym)(x,1,false)
        ($sym)(x::Char, p::Int)     = ($sym)(unsigned(x),p,false)
        ($sym)(x::Char)             = ($sym)(unsigned(x),1,false)
        ($sym)(x::Integer, p::Int)  = ($sym)(unsigned(abs(x)),p,x<0)
        ($sym)(x::Integer)          = ($sym)(unsigned(abs(x)),1,x<0)
    end
end

bits(x::Union(Bool,Int8,Uint8))           = bin(reinterpret(Uint8,x),8)
bits(x::Union(Int16,Uint16,Float16))      = bin(reinterpret(Uint16,x),16)
bits(x::Union(Char,Int32,Uint32,Float32)) = bin(reinterpret(Uint32,x),32)
bits(x::Union(Int64,Uint64,Float64))      = bin(reinterpret(Uint64,x),64)
bits(x::Union(Int128,Uint128))            = bin(reinterpret(Uint128,x),128)

function digits{T<:Integer}(n::Integer, base::T=10, pad::Int=1)
    2 <= base || error("invalid base: $base")
    m = max(pad,ndigits0z(n,base))
    a = zeros(T,m)
    for i = 1:m
        a[i] = rem(n,base)
        n = div(n,base)
    end
    return a
end

isqrt(x::Integer) = oftype(x, trunc(sqrt(x)))

function isqrt(x::Union(Int64,Uint64,Int128,Uint128))
    x==0 && return x
    s = oftype(x, trunc(sqrt(x)))
    # fix with a Newton iteration, since conversion to float discards
    # too many bits.
    s = (s + div(x,s)) >> 1
    s*s > x ? s-1 : s
end

max_stars_repo_name: jakebolewski/julia

max_issues_repo_name: jakebolewski/julia

max_forks_repo_name: jakebolewski/julia


text:
"""
    generate_realizations(n::Int,
                          domdist::DomainDistribution)

Generate `n` realizations from the domain distribution `domdist`, return as a
vector of matrices. Currently restricted to `Float64` elements. Useful for
generating habitats and spatially varying catchabilities that are realization within
replicates.
"""
function generate_realizations(domdist::DomainDistribution, n::Int)
    realizations = Vector{Matrix{Float64}}()
    for rlz in 1:n
        push!(realizations, rand(domdist))
    end
    realizations
end

"""
    save_realizations(fn::AbstractString,
                      habitats::Vector{Matrix{Float64}},
                      name::AbstractString)

Save a vector of `DomainDistribution` realizations as a 3-dimensional array
using the HDF5 package as a data set `name` in the file `fn`.
"""
function save_realizations(fn::AbstractString,
                           name::AbstractString,
                           realizations::Vector{Matrix{Float64}})
    fid = h5open(fn, "cw")
    try
        _dims = size(realizations[1])
        _n = length(realizations)
        dset = create_dataset(fid, name, datatype(Float64),
                              dataspace(_dims..., _n),
                              chunk = (_dims..., 1))
        for (idx, rlzn) in enumerate(realizations)
            dset[:, :, idx] = rlzn
        end
        return nothing
    catch
        @error "Problem encountered; realizations probably not saved."
    finally
        close(fid)
    end
end

"""
    load_realization(fn::AbstractString,
                     name::AbstractString,
                     n::Integer)

Load the `n`th realization of the dataset `name` in the HDF5 file `fn`.
"""
function load_realization(fn::AbstractString,
                          name::AbstractString,
                          n::Integer)
    h5open(fn, "r")[name][:, :, n]
end

"""
    save_movements(fn::AbstractString,
                   moves::Vector{M}) where M<:MovementModel

Save vector of `MovementModel` objects as 3-dimensional HDF5 array.
"""
function save_realizations(fn::AbstractString, name::AbstractString,
                           moves::Vector{M}) where M<:MovementModel
    move_ops = getfield.(moves, :M)
    save_realizations(fn, name, move_ops)
end

function load_movement(fn::AbstractString, name::AbstractString, n::Integer)
    M = load_realization(fn, name, n)
    MovementModel(M, isqrt.(size(M)))
end

function save_realizations(fn::AbstractString, name::AbstractString,
                           popstates::Vector{P}) where P<:PopState
    pops = getfield.(popstates, :P)
    save_realizations(fn, name, pops)
end

function load_spatial_eq(fn::AbstractString, name::AbstractString, n::Integer)
    PopState(load_realization(fn, "spatial_eq", n))
end

function load_popstate(fn::AbstractString, name::AbstractString, n::Integer)
    PopState(load_realization(fn, name * "/popstate", n))
end

function save_realizations(fn::AbstractString, name::AbstractString,
                           spatqs::Vector{C}) where C<:Catchability
    qs = getfield.(spatqs, :catchability)
    save_realizations(fn, name, qs)
end

function load_catchability(fn::AbstractString, name::AbstractString, n::Integer)
    Catchability(load_realization(fn, name, n))
end


text:
module Day11

using AdventOfCode2019

struct Panel
    position::CartesianIndex
    color::UInt8
end

const directions = Dict(
    0 => (1, 0),   # right
    1 => (0, -1),  # down
    2 => (-1, 0),  # left
    3 => (0, 1)    # up
)

function day11(input::String = readInput(joinpath(@__DIR__, "..", "data", "day11.txt")))
    program = parse.(Int, split(input, ","))

    path = run_painting_robot(program, 0)
    nPaintedPanels = length(unique([p.position for p in path]))

    path = run_painting_robot(program, 1)
    minX = minimum(x->x.position[1], path)
    maxX = maximum(x->x.position[1], path)
    minY = minimum(x->x.position[2], path)
    maxY = maximum(x->x.position[2], path)
    panelMap = zeros(Int, (maxX-minX+1, maxY-minY+1))
    for panel in path
        panelMap[panel.position[1]-minX+1, (maxY-minY+1)-(panel.position[2]-minY)] = (panel.color == 1) ? 1 : 0
    end

    return [nPaintedPanels, panelMap]
end

function run_painting_robot(program::Array{Int,1}, initialColor::Int)
    path = Array{Panel,1}()
    position = CartesianIndex(0, 0)
    direction = 3

    status = Channel{Bool}(1)
    input = Channel{Int}(1)
    put!(input, initialColor)
    output = Channel{Int}(2)

    t = @task AdventOfCode2019.IntCode.run_program!(program, input, output, nothing)
    schedule(t)

    while !istaskdone(t)
        color = take!(output)
        turn = take!(output)
        push!(path, Panel(position, color))
        if turn == 0
            direction = turn_left(direction)
        elseif turn == 1
            direction = turn_right(direction)
        end
        position = move(position, direction)
        colorIn = 0
        for i = length(path):-1:1
            if path[i].position == position
                colorIn = path[i].color
                break
            end
        end
        put!(input, colorIn)
    end
    close(output)
    close(input)
    return path
end

function turn_right(dir::Int)
    return mod(dir + 1, 4)
end

function turn_left(dir::Int)
    return mod(dir - 1, 4)
end

function move(pos::CartesianIndex, dir::Int)
    return CartesianIndex(pos.I .+ directions[dir])
end

end # module


text:
function init_constraint_struct(com, ::CPE.AllEqual, internals)
    EqualConstraint(internals, ones(Int, length(internals.indices)))
end

implements_activate(::Type{EqualConstraint}, ::Type{MOI.VectorOfVariables}, ::Type{CS.CPE.AllEqual}) = true

function activate_constraint!(
    com::CS.CoM,
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CS.CPE.AllEqual
)
    indices = constraint.indices
    search_space = com.search_space
    intersect_vals = Set(intersect(CS.values.(search_space[indices])...))
    if isempty(intersect_vals)
        return false
    end
    for vidx in indices
        for val in CS.values(search_space[vidx])
            if !(val in intersect_vals)
                !rm!(com, search_space[vidx], val) && return false
            end
        end
    end

    return true
end

function apply_changes!(
    com::CS.CoM,
    v::Variable,
    changes,
    first_ptr::Int,
)
    for i in first_ptr:length(changes)
        change = changes[i]
        if change[1] == :remove_below
            !remove_below!(com, v, change[2]) && return false
        elseif change[1] == :remove_above
            !remove_above!(com, v, change[2]) && return false
        elseif change[1] == :rm && has(v, change[2])
            !rm!(com, v, change[2]) && return false
        end
    end
    return true
end

"""
    prune_constraint!(com::CS.CoM, constraint::EqualConstraint, fct::MOI.VectorOfVariables, set::CPE.AllEqual; logs = true)

Reduce the number of possibilities given the equality constraint which sets all variables in `MOI.VectorOfVariables` to the same value.
Return if still feasible and throw a warning if infeasible and `logs` is set to `true`
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CPE.AllEqual;
    logs = true,
)
    indices = constraint.indices

    search_space = com.search_space
    # is only needed if we want to set more
    if length(indices) > 2
        fixed_vals, unfixed_indices = fixed_vs_unfixed(search_space, indices)

        fixed_vals_set = Set(fixed_vals)
        # check if one value is used more than once
        if length(fixed_vals_set) > 1
            logs && @warn "The problem is infeasible"
            return false
        elseif length(fixed_vals_set) == 0
            # sync the changes in each variable
            for i in 1:length(indices)
                v1 = search_space[indices[i]]
                v1_changes = view_changes(v1, com.c_step_nr)
                isnothing(v1_changes) && continue
                for j in 1:length(indices)
                    i == j && continue
                    v2 = search_space[indices[j]]
                    !apply_changes!(com, v2, v1_changes, constraint.first_ptrs[i]) &&
                        return false
                end
                constraint.first_ptrs[i] = length(v1_changes) + 1
            end
            return true
        end

        # otherwise prune => set all variables to fixed value
        for i in unfixed_indices
            vidx = indices[i]
            feasible = fix!(com, search_space[vidx], fixed_vals[1])
            if !feasible
                return false
            end
        end
    else # faster for two variables
        v1 = search_space[indices[1]]
        v2 = search_space[indices[2]]
        fixed_v1 = isfixed(v1)
        fixed_v2 = isfixed(v2)
        if !fixed_v1 && !fixed_v2
            num_changes_v1 = num_changes(v1, com.c_step_nr)
            num_changes_v2 = num_changes(v2, com.c_step_nr)
            if num_changes_v1 == 0 && num_changes_v2 == 0
                return true
            end
            !apply_changes!(com, v2, view_changes(v1, com.c_step_nr), constraint.first_ptrs[1]) && return false
            !apply_changes!(com, v1, view_changes(v2, com.c_step_nr), constraint.first_ptrs[2]) && return false
            constraint.first_ptrs[1] = num_changes(v1, com.c_step_nr) + 1
            constraint.first_ptrs[2] = num_changes(v2, com.c_step_nr) + 1
            return true
        elseif fixed_v1 && fixed_v2
            if CS.value(v1) != CS.value(v2)
                return false
            end
            return true
        end
        # one is fixed and one isn't
        if fixed_v1
            fix_v = 2
            feasible = fix!(com, v2, CS.value(v1))
            if !feasible
                return false
            end
        else
            feasible = fix!(com, v1, CS.value(v2))
            if !feasible
                return false
            end
            fix_v = 1
        end
    end
    return true
end

"""
    finished_pruning_constraint!(com::CS.CoM,
        constraint::EqualConstraint,
        fct::MOI.VectorOfVariables,
        set::CPE.AllEqual)

Reset the first_ptrs to one for the next pruning step
"""
function finished_pruning_constraint!(
    com::CS.CoM,
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CPE.AllEqual,
)

    constraint.first_ptrs .= 1
end

"""
    still_feasible(com::CoM, constraint::EqualConstraint, fct::MOI.VectorOfVariables, set::CPE.AllEqual, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled.
"""
function still_feasible(
    com::CoM,
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CPE.AllEqual,
    vidx::Int,
    value::Int,
)
    variables = com.search_space
    for cvidx in constraint.indices
        if cvidx == vidx
            continue
        end
        v = variables[cvidx]
        !has(v, value) && return false
    end
    return true
end

function is_constraint_solved(
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CPE.AllEqual,
    values::Vector{Int},
)
    return all(v -> v == values[1], values)
end

"""
    is_constraint_violated(
        com::CoM,
        constraint::EqualConstraint,
        fct::MOI.VectorOfVariables,
        set::CPE.AllEqual
    )

Checks if the constraint is violated as it is currently set. This can happen inside an
inactive reified or indicator constraint.
"""
function is_constraint_violated(
    com::CoM,
    constraint::EqualConstraint,
    fct::MOI.VectorOfVariables,
    set::CPE.AllEqual,
)
    found_first_val = false
    first_val = 0
    for var in com.search_space[constraint.indices]
        if isfixed(var)
            found_first_val = true
            first_val = CS.value(var)
            break
        end
    end
    !found_first_val && return false
    return !all(
        CS.value(var) == first_val
        for var in com.search_space[constraint.indices] if isfixed(var)
    )
end


text:
module AtomsBase
using Unitful
using UnitfulAtomic
import PeriodicTable: elements
using StaticArrays

include("interface.jl")
include("properties.jl")
include("flexible_system.jl")
include("atom.jl")
include("atomview.jl")
include("fast_system.jl")

end