text:
module Probs
include("./samplesGenerator.jl")
export generateChanges
# export ifDetected, notDetected, extractFps, extractTps, belowH, aboveH

function ifDetected(cde::Int64, changes::Vector{Int64}, delay::Int64)
    # 1) If change detection event is TP
    for chp in changes
        if abs(chp - cde) <= delay
            return true
        end
    end
    return false
end

function notDetected(cde::Int64, changes::Vector{Int64}, delay::Int64)
    # 2) If change detection event is not TP
    return ! ifDetected(cde, changes, delay)
end

function extractFps(detections::Vector{Int64},
                    trueChanges::Vector{Int64},
                    delay::Int64)
    # 3) Extract False Positives from a sequence of CDEs
    filter(x->notDetected(x, trueChanges, delay), detections)
end

function extractTps(detections::Vector{Int64},
                    trueChanges::Vector{Int64},
                    delay::Int64)
    # 4) Extract True Positives from a sequence of CDEs
    filter(x->ifDetected(x, trueChanges, delay), detections)
end

function ifPredictedTp(cde::Int64, probs::Vector{Float64}, h::Float64)
    # 5) Check if an individual TP is correctly predicted
    probs[cde] >= h
end

function ifPredictedFp(cde::Int64, probs::Vector{Float64}, h::Float64)
    # 6) Check if an individual FP is correctly ignored
    probs[cde] < h
end

function correctlyPredictedTps(probs::Vector{Float64},
                               trueChanges::Vector{Int64},
                               h::Float64)
    # 7) Given probs vector count those provided true changes for
    # which corresponding probs >= h. Time locations must be integers
    filter(x->ifPredictedTp(x, probs, h), trueChanges)
end

function correctlyPredictedFps(probs::Vector{Float64},
                               fPositives::Vector{Int64},
                               h::Float64)
    # 8) Given probs vector count those provided false alarms for
    # which corresponding probs < h. Time locations must be integers
    filter(x->ifPredictedFp(x, probs, h), fPositives)
end

function belowH(h::Float64, probs::Vector{Float64}, delay::Int64)
    # 9) Returns _indices_ given a vector of probs and threshold
    # function. Returs indices of those elements in probs vector which
    # are below threshold value
    ret :: Vector{Int64}
    subprobs = probs[1 : min(delay, length(probs))]
    ret = find( x -> x <= h, subprobs)
end

function aboveH(h, probs, delayConfirm)
    # 10)
    subprobs = probs[1 : min(delayConfirm, length(probs))]
    find(x->(x>h), subprobs)
end

#function initProbs(mu, sigma)
#    # 11)
#    0
#end
#
#function whichBin(num, edges)
#    # 12) Edges:[1,2,3] Bins:[(-inf,1) (1,2) (2,3) (3, inf)] Bins nums: 1,2,3,4
#    n = length(edges)
#    if num <= edges[1]
#        return 1
#    elseif num > edges[end]
#        return n+1
#    else
#        for i in 2:n
#            if edges[i-1] < num <= edges[i]
#                return i
#            end
#        end
#    end
#end

function gaussPdf(x, mu, sigma)
    k = 1/(sigma * sqrt(2*pi))
    p = k * exp(-((x-mu)^2/(2*sigma^2)))
end

# 13, 14 - don't need
function generateChanges(n, mu, sigma, width)
    # 15) Generate sequence of changes
    changes::Array{Int,1}
    values = 1:width
    weights = map(x->gaussPdf(x, mu, sigma), 1:width)
    samples = genSamples(values, weights, n)
    changes = cumsum(samples)
end

function dynamicH(probs,threshold,hmin,hmax,tmax)
    dynamich::Vector{Float64}
    dynamich=map(probs) do p
        if p>=threshold
            hmin
        else
            hmax
        end
    end
    dynamich[1:tmax]
end

end
#------------------ End module Probs-------------------------#

#------------------------------------------------------------#
# module Pmf
# export gaussPmf
#
# using Distributions
# using Bisect
#
# type GaussPmf
#     # Prototype is in thinkbayes.py/def MakeGaussianPmf(mu, sigma, num_sigmas, n=201):
#     mu         :: Float64
#     sigma      :: Float64
#     num_sigmas :: Int
#     n          :: Int
#     edges      :: Vector{Float64}
#     probs      :: Vector{Float64}
#     function GaussPmf(mu::Float64, sigma::Float64)
#         n_default          = 251
#         num_sigmas_default = 3
#         low                = mu - num_sigmas_default * sigma
#         high               = mu + num_sigmas_default * sigma
#         d                  = Normal(mu, sigma)
#         edges              = linspace(low, high, n_default)
#         probs              = pdf(d, edges)
#         new(mu, sigma, num_sigmas_default, n_default, edges, probs)
#     end
# end
#
# function gaussPmf(gp::GaussPmf, val)
#      Bisect.idx_of_closest(gp.probs, float(val))
# end
#
#end
#----------------- End module Pmf ------------------------------#

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
function setref!(x::EXPR, binding::Binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
    binding.refs !== nothing && push!(binding.refs, x)
end

function setref!(x::EXPR, binding)
    if !hasmeta(x)
        x.meta = Meta()
    end
    x.meta.ref = binding
end

# Main function to be called. Given the `state` tries to determine what `x`
# refers to. If it remains unresolved and is in a delayed evaluation scope
# (i.e. a function) it gets pushed to list (.urefs) to be resolved after we've
# run over the entire top-level scope.
function resolve_ref(x, state)
    if !(parentof(x) isa EXPR && typof(parentof(x)) == CSTParser.Quotenode)
        resolved = resolve_ref(x, state.scope, state)
    end
end

# The first method that is tried. Searches the current scope for local bindings
# that match `x`. Steps:
# 1. Check whether we've already checked this scope (inifinite loops are
# possible when traversing nested modules.)
# 2. Check what sort of EXPR we're dealing with, separate name from EXPR that
# binds.
# 3. Look in the scope's variable list for a binding matching the name.
# 4. If 3. is unsuccessful, check whether the scope imports any modules then check them.
# 5. If no match is found within this scope check the parent scope.
# The return value is a boolean that is false if x should point to something but
# can't be resolved.

function resolve_ref(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = false

    if is_getfield(x)
        return resolve_getfield(x, scope, state)
    elseif is_kwarg(x)
        # Note to self: this seems wronge - Binding should be attached to entire Kw EXPR.
        if isidentifier(x[1])
            setref!(x[1], Binding(noname, nothing, nothing, [], nothing, nothing))
        elseif is_declaration(x[1]) && isidentifier(x[1][1])
            setref!(x[1][1], Binding(noname, nothing, nothing, [], nothing, nothing))
        end
        return true
    elseif is_special_macro_term(x) || new_within_struct(x)
        setref!(x, Binding(noname, nothing, nothing, [], nothing, nothing))
        return true
    end

    x1, mn = nameof_expr_to_resolve(x)
    mn == true && return true

    if scopehasbinding(scope, mn)
        setref!(x1, scope.names[mn])
        resolved = true
    elseif scope.modules isa Dict && length(scope.modules) > 0
        for m in values(scope.modules)
            resolved = resolve_ref_from_module(x, m, state)
            resolved && return true
        end
    end
    if !resolved && !CSTParser.defines_module(scope.expr) && parentof(scope) isa Scope
        return resolve_ref(x, parentof(scope), state)
    end
    return resolved
end

# Searches a module store for a binding/variable that matches the reference `x1`.
function resolve_ref_from_module(x1::EXPR, m::SymbolServer.ModuleStore, state::State)::Bool
    hasref(x1) && return true
    if isidentifier(x1)
        x = x1
        if Symbol(valof(x)) == m.name.name
            setref!(x, m)
            return true
        elseif isexportedby(x, m)
            setref!(x, maybe_lookup(m[Symbol(valof(x))], state.server))
            return true
        end
    elseif is_macroname(x1)
        x = x1[2]
        if valof(x) == "." && m.name == VarRef(nothing, :Base)
            # @. gets converted to @__dot__, probably during lowering.
            setref!(x, m[:Broadcast][Symbol("@__dot__")])
            return true
        end

        mn = Symbol("@", valof(x))
        if isexportedby(mn, m)
            setref!(x, maybe_lookup(m[mn], state.server))
            return true
        end
    elseif typof(x1) === x_Str
        mac = x1
        mn = Symbol("@", valof(mac[1]), "_str")
        if isexportedby(mn, m)
            setref!(mac[1], maybe_lookup(m[mn], state.server))
            return true
        end
    end
    return false
end

function resolve_ref_from_module(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = false

    x1, mn = nameof_expr_to_resolve(x)
    mn == true && return true

    if scope_exports(scope, mn, state)
        setref!(x1, scope.names[mn])
        resolved = true
    end
    return resolved
end

"""
    scope_exports(scope::Scope, name::String)

Does the scope export a variable called `name`?
"""
function scope_exports(scope::Scope, name::String, state)
    if scopehasbinding(scope, name) && (b = scope.names[name]) isa Binding
        initial_pass_on_exports(scope.expr, state.server)
        for ref in b.refs
            if ref isa EXPR && parentof(ref) isa EXPR && typof(parentof(ref)) === CSTParser.Export
                return true
            end
        end
    end
    return false
end

"""
    initial_pass_on_exports(x::EXPR, server)

Export statements need to be (pseudo) evaluated each time we consider
whether a variable is made available by an import statement.
"""
function initial_pass_on_exports(x::EXPR, server)
    # @assert CSTParser.defines_module(x)
    for a in x[3] # module block expressions
        if typof(a) === CSTParser.Export
            traverse(a, Delayed(retrieve_scope(a), server))
        end
    end
end

# Fallback method
function resolve_ref(x::EXPR, m, state::State)::Bool
    return hasref(x)::Bool
end

rhs_of_getfield(x::EXPR) = is_getfield_w_quotenode(x) ? isquoted(x[3]) ? x[3][2] : x[3][1] : x
lhs_of_getfield(x::EXPR) = rhs_of_getfield(x[1])

"""
    resolve_getfield(x::EXPR, parent::Union{EXPR,Scope,ModuleStore,Binding}, state::State)::Bool

Given an expression of the form `parent.x` try to resolve `x`. The method
called with `parent::EXPR` resolves the reference for `parent`, other methods
then check whether the Binding/Scope/ModuleStore to which `parent` points has
a field matching `x`.
"""
function resolve_getfield(x::EXPR, scope::Scope, state::State)::Bool
    hasref(x) && return true
    resolved = resolve_ref(x[1], scope, state)
    if isidentifier(x[1])
        lhs = x[1]
    elseif is_getfield_w_quotenode(x[1])
        lhs = lhs_of_getfield(x)
    else
        return resolved
    end
    if resolved && (rhs = rhs_of_getfield(x)) !== nothing
        resolved = resolve_getfield(rhs, refof(lhs), state)
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type::EXPR, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x)
        if CSTParser.defines_module(parent_type) && scopeof(parent_type) isa Scope
            resolved = resolve_ref(x, scopeof(parent_type), state)
        elseif CSTParser.defines_struct(parent_type)
            if scopehasbinding(scopeof(parent_type), valof(x))
                setref!(x, scopeof(parent_type).names[valof(x)])
                resolved = true
            end
        end
    end
    return resolved
end

function resolve_getfield(x::EXPR, b::Binding, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if b.val isa Binding
        resolved = resolve_getfield(x, b.val, state)
    elseif b.val isa SymbolServer.ModuleStore || (b.val isa EXPR && CSTParser.defines_module(b.val))
        resolved = resolve_getfield(x, b.val, state)
    elseif b.type isa Binding
        resolved = resolve_getfield(x, b.type.val, state)
    elseif b.type isa SymbolServer.DataTypeStore
        resolved = resolve_getfield(x, b.type, state)
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent_type, state::State)::Bool
    hasref(x) && return true
    return false
end

function is_overloaded(val::SymbolServer.SymStore, scope::Scope)
    (vr = val.name isa SymbolServer.FakeTypeName ? val.name.name : val.name)
    haskey(scope.overloaded, vr)
end

function resolve_getfield(x::EXPR, m::SymbolServer.ModuleStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x) && (val = maybe_lookup(SymbolServer.maybe_getfield(Symbol(CSTParser.str_value(x)), m, getsymbolserver(state.server)), state.server)) !== nothing
        # Check whether variable is overloaded in top-level scope
        tls = retrieve_toplevel_scope(state.scope)
        if tls.overloaded !== nothing  && (vr = val.name isa SymbolServer.FakeTypeName ? val.name.name : val.name; haskey(tls.overloaded, vr))
            setref!(x, tls.overloaded[vr])
            return true
        end
        setref!(x, val)
        resolved = true
    elseif is_macroname(x) && (val = maybe_lookup(SymbolServer.maybe_getfield(Symbol("@", CSTParser.str_value(x[2])), m, getsymbolserver(state.server)), state.server)) !== nothing
        setref!(x[2], val)
        resolved = true
    end
    return resolved
end

function resolve_getfield(x::EXPR, parent::SymbolServer.DataTypeStore, state::State)::Bool
    hasref(x) && return true
    resolved = false
    if isidentifier(x) && Symbol(valof(x)) in parent.fieldnames
        fi = findfirst(f -> Symbol(valof(x)) == f, parent.fieldnames)
        ft = parent.types[fi]
        val = SymbolServer._lookup(ft, getsymbolserver(state.server), true)
        # TODO: Need to handle the case where we get back a FakeUnion, etc.
        setref!(x, Binding(noname, nothing, val, [], nothing, nothing))
        resolved = true
    end
    return resolved
end

resolvable_macroname(x::EXPR) = is_macroname(x) && isidentifier(x[2]) && refof(x[2]) === nothing

"""
    module_safety_trip(scope::Scope,  visited_scopes)

Checks whether the scope is a module and we've visited it before,
otherwise adds the module to the list.
"""
function module_safety_trip(scope::Scope,  visited_scopes)
    if CSTParser.defines_module(scope.expr) && length(scope.expr) > 1 && isidentifier(scope.expr[2])
        s_m_name = valofid(scope.expr[2])
        if s_m_name in visited_scopes
            return true
        else
            push!(visited_scopes, s_m_name)
        end
    end
    return false
end

function nameof_expr_to_resolve(x)
    if isidentifier(x)
        x1 = x
        mn = valofid(x)
    elseif resolvable_macroname(x)
        x1 = x[2]
        mn = string("@", valofid(x1))
    elseif typof(x) === x_Str && isidentifier(x[1])
        x1 = x[1]
        mn = string("@", valofid(x1), "_str")
    else
        return x, true
    end
    x1, mn
end

"""
    valofid(x)

Returns the string value of an expression for which `isidentifier` is true,
i.e. handles NONSTDIDENTIFIERs.
"""
valofid(x::EXPR) = typof(x) === CSTParser.IDENTIFIER ? valof(x) : valof(x[2])

"""
new_within_struct(x::EXPR)

Checks whether x is a reference to `new` within a datatype constructor.
"""
new_within_struct(x::EXPR) = isidentifier(x) && valofid(x) == "new" && is_in_fexpr(x, CSTParser.defines_struct)
is_special_macro_term(x::EXPR) = isidentifier(x) && (valofid(x) == "__source__" || valofid(x) == "__module__") && is_in_fexpr(x, CSTParser.defines_macro)

ext: jl


text:
using Test, ChainRulesCore, ChainRulesTestUtils

@nospecialize

using Base.Broadcast: broadcastable
using ChainRules
using ChainRulesCore
using ChainRulesTestUtils
using ChainRulesTestUtils: rand_tangent, _fdm
using FiniteDifferences
using LinearAlgebra
using LinearAlgebra.BLAS
using LinearAlgebra: dot
using Random
using SparseArrays
using StaticArrays
using Statistics
using Test
using JuliaInterpreter

union!(JuliaInterpreter.compiled_modules, Any[Base, Base.Broadcast, LinearAlgebra, Random, StaticArrays, Statistics])

Random.seed!(1) # Set seed that all testsets should reset to.

function include_test(path)
    if isempty(ARGS) || any(occursin(a, path) for a in ARGS)
        println("Testing $path:")  # print so TravisCI doesn't timeout due to no output
        @time Base.include(@__MODULE__(), path) do ex
            Meta.isexpr(ex, :macrocall) && ex.args[1] == Symbol("@testset") || return ex
            return :(@interpret (() -> $ex)())  # interpret testsets using JuliaInterpreter
        end
    else
        # If you provide ARGS like so, then it runs only matching testsets:
        # Pkg.test("ChainRules", test_args = ["index", "LinearAlgebra"])
        println("(Not testing $path)")
    end
end

if isempty(ARGS)
    println("Testing ChainRules.jl")
else
    println("Testing ChainRules.jl with test_args = ", ARGS)
end

@testset "ChainRules" begin  # One overall @testset ensures it keeps going after failures
    include("test_helpers.jl")
    println()

    test_method_tables()  # Check the global method tables are consistent

    # Each file puts all tests inside one or more @testset blocks
    include_test("rulesets/Base/base.jl")
    include_test("rulesets/Base/fastmath_able.jl")
    include_test("rulesets/Base/evalpoly.jl")
    include_test("rulesets/Base/array.jl")
    include_test("rulesets/Base/arraymath.jl")
    include_test("rulesets/Base/indexing.jl")
    include_test("rulesets/Base/mapreduce.jl")
    include_test("rulesets/Base/sort.jl")

    println()

    include_test("rulesets/Statistics/statistics.jl")

    println()

    include_test("rulesets/LinearAlgebra/dense.jl")
    include_test("rulesets/LinearAlgebra/norm.jl")
    include_test("rulesets/LinearAlgebra/matfun.jl")
    include_test("rulesets/LinearAlgebra/structured.jl")
    include_test("rulesets/LinearAlgebra/symmetric.jl")
    include_test("rulesets/LinearAlgebra/factorization.jl")
    include_test("rulesets/LinearAlgebra/blas.jl")
    include_test("rulesets/LinearAlgebra/lapack.jl")
    include_test("rulesets/LinearAlgebra/uniformscaling.jl")

    println()

    include_test("rulesets/SparseArrays/sparsematrix.jl")

    println()

    include_test("rulesets/Random/random.jl")
    println()
end


text:
struct IterativeInitializationCollisionModel{T <: Real, UT <: AbstractArray{T, 3}} <:
       CollisionModel
    τ::T

    # The incompressible equilibrium funciton can be split in a linear and a nonlinear component
    # the linear component depends on the local density while the nonlinear component depends
    # on the local velocity.
    # Since this collision operator keeps a constant velocity for each lattice node we can
    # precompute the nonlinear term for some performance improvements
    nonlinear_term::UT
end
function IterativeInitializationCollisionModel(
    q::Quadrature,
    τ::T,
    problem,
) where {T <: Real}
    nx, ny, nf = problem.NX, problem.NY, length(q.weights)
    nonlinear_term = Array{T}(undef, nx, ny, nf)
    ρ_0 = one(T)

    # Initialize f with the given velocity field
    x_range, y_range = range(problem)
    for x_idx in 1:nx, y_idx in 1:ny
        u = lattice_velocity(q, problem, x_range[x_idx], y_range[y_idx])
        u_squared = u[1]^2 + u[2]^2

        for f_idx in 1:nf
            u_dot_xi = dot(q.abscissae[:, f_idx], u)

            cs = q.speed_of_sound_squared
            a_H_1 = cs * u_dot_xi
            a_H_2 = cs^2 * (u_dot_xi * u_dot_xi) - cs * u_squared

            nonlinear_term[x_idx, y_idx, f_idx] =
                q.weights[f_idx] * ρ_0 * (a_H_1 + a_H_2 / 2)
        end
    end

    return IterativeInitializationCollisionModel(τ, nonlinear_term)
end
function collide!(
    collision_model::IterativeInitializationCollisionModel,
    q::Quadrature,
    f_in,
    f_out;
    time = 0.0,
)
    τ = collision_model.τ
    nx, ny, nf = size(f_out)
    for x_idx in 1:nx, y_idx in 1:ny
        ρ = density(q, f_in[x_idx, y_idx, :])
        @inbounds for f_idx in 1:nf
            feq = q.weights[f_idx] * ρ + collision_model.nonlinear_term[x_idx, y_idx, f_idx]

            f_out[x_idx, y_idx, f_idx] =
                (1 - 1 / τ) * f_in[x_idx, y_idx, f_idx] + (1 / τ) * feq
        end
    end
    return
end


text:
__precompile__()

module SimpleStructs

export @defstruct, @defimmutable

"""A convenient macro copied from Mocha.jl that could be used to define structs
with default values and type checks. For example
```julia
@defstruct MyStruct (
  field1 :: Int = 0,
  (field2 :: AbstractString = "", !isempty(field2))
)
```
where each field could be either
```julia
field_name :: field_type = default_value
```
or put within a tuple, with the second element
specifying a validation check on the field value.
In the example above, the default value for
field2 does not satisfy the assertion, this
could be used to force user to provide a
valid value when no meaningful default value
is available.
The macro will define a constructor that could accept
the keyword arguments.
"""
macro defstruct(name, fields)
  _defstruct_impl(false, name, fields)
end

"""A convenient macro to define immutable structs. The same as
`@defstruct` except that the defined type is immutable.
"""
macro defimmutable(name, fields)
  _defstruct_impl(true, name, fields)
end

"""Internal use only, this value is used to indicate a required value
is not specified.
"""
immutable __Undefined
end

function _contains_type_param(expr, tp::AbstractVector)
  for t in tp
    if _contains_type_param(expr, t)
      return true
    end
  end
  return false
end

function _contains_type_param(expr, tp::Symbol)
  if isa(expr, Symbol)
    return expr == tp
  else
    for t in expr.args
      if _contains_type_param(t, tp)
        return true
      end
    end
  end
  return false
end

function _type_param_name(expr)
  if isa(expr, Symbol)
    expr
  else
    @assert isa(expr, Expr) && expr.head == :<:
    expr.args[1]
  end
end

function _name_to_call(name; func_name=nothing)
    if isa(name, Expr) && name.head == :curly
        type_param_names = map(_type_param_name, name.args[2:end])
        Expr(:curly, isa(func_name, Void)? name.args[1] : func_name, type_param_names...)
    else
        isa(func_name, Void)? name : func_name
    end
end

# from Foo{T,T2<:Real} to Foo{T,T2}
function _name_to_type(name)
  if isa(name, Expr) && name.head == :curly
    Expr(:curly, name.args[1], map(x -> isa(x,Symbol)?x:x.args[1], name.args[2:end])...)
  else
    name
  end
end

function _defstruct_impl(is_immutable, name, fields)
  if isa(fields, Expr) && fields.head == :tuple
    fields = fields.args
  else
    fields = [fields]
  end
  @assert length(fields) > 0

  type_param_names = Array(Symbol, 0)

  if isa(name, Symbol)
    super_name = :Any
    name       = name
  elseif isa(name, Expr) && name.head == :curly
    type_param_names = map(_type_param_name, name.args[2:end]) # :T
    super_name = :Any
    name       = name
  else
    @assert(isa(name, Expr) && name.head == :comparison && length(name.args) == 3 && name.args[2] == :(<:),
            "name must be of form 'Name <: SuperType'")
    #@assert(isa(name.args[1], Symbol) && isa(name.args[3], Symbol))
    if isa(name.args[1], Expr) && name.args[1].head == :curly
        type_param_names = map(_type_param_name, name.args[1].args[2:end]) # :T
    end
    super_name = name.args[3]
    name       = name.args[1]
  end

  field_defs     = Array(Expr, length(fields))        # :(field2 :: Int)
  field_names    = Array(Symbol, length(fields))      # :field2
  field_defaults = Array(Expr, length(fields))        # :(field2 = 0)
  field_types    = Array(Any, length(fields))         # Int
  field_asserts  = Array(Expr, length(fields))        # :(field2 >= 0)
  required_field = Symbol[]

  for i = 1:length(fields)
    field = fields[i]
    if field.head == :tuple
      field_asserts[i] = field.args[2]
      field = field.args[1]
    end
    if field.head == :(=)
      fname             = field.args[1].args[1]
      field_defs[i]     = field.args[1]
      field_names[i]    = fname
      field_types[i]    = field.args[1].args[2]
      if _contains_type_param(field.args[1].args[2], type_param_names)
        field_defaults[i] = Expr(:kw, Expr(:(::), fname, field.args[1].args[2]), field.args[2])
      else
        field_defaults[i] = Expr(:kw, fname, field.args[2])
      end
    else
      # no default value provided, required field
      fname             = field.args[1]
      field_defs[i]     = field
      field_names[i]    = fname
      field_types[i]    = field.args[2]
      if _contains_type_param(field.args[2], type_param_names)
        field_defaults[i] = Expr(:kw, Expr(:(::), fname, field.args[2]), __Undefined())
      else
        field_defaults[i] = Expr(:kw, fname, __Undefined())
      end
      push!(required_field, fname)
    end
  end

  # body of layer type, defining fields
  type_body = Expr(:block, field_defs...)

  # constructor
  requires = map(required_field) do fname
    :(@assert(!isa($fname, SimpleStructs.__Undefined), "value for " * $(string(fname)) * " is required"))
  end
  converts = map(zip(field_names, field_types)) do param
    f_name, f_type = param
    :($f_name = convert($f_type, $f_name))
  end
  asserts = map(filter(i -> isdefined(field_asserts,i), 1:length(fields))) do i
    :(@assert($(field_asserts[i])))
  end
  construct = Expr(:call, _name_to_call(name), field_names...)
  ctor_body = Expr(:block, requires..., converts..., asserts..., construct)
  ctor_def = Expr(:call, name, Expr(:parameters, field_defaults...))
  ctor = Expr(:(=), ctor_def, ctor_body)

  # Base.show function
  show_fields = map(enumerate(field_names)) do i_fname
    i, fname = i_fname
    sep = i == 1 ? :() : :(print(io, ", "))
    quote
      $sep
      print(io, $(string(fname)), "=")
      show(io, obj.$fname)
    end
  end
  show_fields = Expr(:block, show_fields...)
  defshow = quote
    function $(_name_to_call(name, func_name=:show))(io::IO, obj::$(_name_to_type(name)))
      print(io, typeof(obj))
      print(io, "(")
      $show_fields
      print(io, ")")
    end
  end

  if is_immutable
    s_def =  quote
      immutable $(name) <: $(super_name)
        $type_body
      end
    end
  else
    s_def = quote
      type $(name) <: $(super_name)
        $type_body
      end
    end
  end

  return esc(quote
    $s_def
    $ctor

    import Base.show
    $defshow
  end)
end

end # module


text:
using Documenter

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

using StochasticOptimalTransport

DocMeta.setdocmeta!(
    StochasticOptimalTransport,
    :DocTestSetup,
    :(using StochasticOptimalTransport);
    recursive=true,
)

makedocs(;
    modules=[StochasticOptimalTransport],
    authors="David Widmann <david.widmann@it.uu.se>",
    repo="https://github.com/JuliaOptimalTransport/StochasticOptimalTransport.jl/blob/{commit}{path}#L{line}",
    sitename="StochasticOptimalTransport.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://juliaoptimaltransport.github.io/StochasticOptimalTransport.jl",
        assets=String[],
    ),
    pages=["Home" => "index.md"],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/JuliaOptimalTransport/StochasticOptimalTransport.jl",
    push_preview=true,
    devbranch="main",
)


text:
@testset "basic predicate tests" begin
    pred1 = Predicate((x,y) -> x == y)
    pred2 = Predicate((x,y) -> x != y)
    pred3 = Predicate((x,y) -> all(x .== -y))
    @test pred1(1,1)
    @test !pred1(1,2)
    @test !pred2(1,1)
    @test pred2(1,2)
    @test pred3([1,2],[-1,-2])
    @test !pred3([1,2],[2,3])
end

@testset "is predicate equivalence relation" begin
    pred1 = Predicate((x,y) -> x == y)
    pred2 = Predicate((x,y) -> x != y)
    pred3 = Predicate((x,y) -> all(x .== -y))
    pred4 = Predicate((x,y) -> x == y+3)
    @test check_for_equivalence_relation(pred1, -11:10)
    @test check_for_equivalence_relation(pred2, -7:10)
    @test check_for_equivalence_relation(pred3, [[i,j] for i in -7:10 for j in -2:8])
    @test !check_for_equivalence_relation(pred4, -7:10)
end

@testset "mappings" begin
    m1 = Mapping(x -> [x])
    m2 = Mapping(x -> [-x])
    sol2 = Dict(-4 => 1, -3 => 2, 3 => 2, -2 => 3, 2 => 3, -1 => 4, 1 => 4, 0 => 5)
    m3 = Mapping(x -> [x+3,x-3])
    m5 = Mapping(x -> [x-3])
    sol3 = Dict(-4 => 1, -1 => 1, 2 => 1, -3 => 2, 0 => 2, 3 => 2, -2 => 3, 1 => 3)
    @test find_classes(m1, 1:5) == Dict(zip(1:5,collect(1:5)))
    @test find_classes(m2, -4:3) == sol2
    @test find_classes(m3, -4:3) == sol3
    @test !(find_classes(m5, -4:3) == sol3)
end


text:
module TestBuilderRCFile

using Test
using Comonicon
using Comonicon.Builder:
    contains_path, contains_fpath, write_path, write_fpath, detect_rcfile, install_env_path

home_dir = mktempdir()
test_dir = pkgdir(Comonicon, "test")
usr_dir = joinpath(test_dir, "usr")

options = Comonicon.Options.Comonicon(
    name = "test",
    install = Comonicon.Options.Install(path = usr_dir, completion = false),
    sysimg = nothing,
)

@testset "contains path/write path" begin
    withenv("PATH" => nothing, "FPATH" => nothing) do
        dir = mktempdir()
        rcfile = joinpath(dir, ".bashrc")
        touch(rcfile)
        @test contains_path(rcfile, usr_dir, Base.EnvDict()) == false

        write_path(rcfile, usr_dir)
        @test contains_path(rcfile, usr_dir) == true

        @test contains_fpath(rcfile, usr_dir, Base.EnvDict()) == false

        write_fpath(rcfile, usr_dir)
        @test contains_fpath(rcfile, usr_dir) == true
    end

    withenv("PATH" => joinpath(usr_dir, "bin"), "FPATH" => joinpath(usr_dir, "completions")) do
        @test contains_path("rcfile", usr_dir) == true
        @test contains_fpath("rcfile", usr_dir) == true
    end
end

@testset "detect rcfile" begin
    @test detect_rcfile("zsh", home_dir) == joinpath(home_dir, ".zshrc")
    @test withenv("ZDOTDIR" => "test") do
        detect_rcfile("zsh", home_dir)
    end == joinpath("test", ".zshrc")

    @test detect_rcfile("bash", home_dir) == joinpath(home_dir, ".bash_profile")
    touch(joinpath(home_dir, ".bashrc"))
    @test detect_rcfile("bash", home_dir) == joinpath(home_dir, ".bashrc")
end

@testset "install env path" begin
    install_env_path(Main, options; shell = "zsh", home_dir, env = ENV, yes = true)

    withenv("PATH" => nothing, "FPATH" => nothing) do
        @test contains_path(joinpath(home_dir, ".zshrc"), usr_dir)
        @test contains_fpath(joinpath(home_dir, ".zshrc"), usr_dir)
    end
end

end # TestBuilderRCFile


text:
function update_da(da, BA, τ; δt=1)

	δda = (-da ./ τ) .+ BA

	da = da + δda .* δt
end

function update_ACh!(t, l, m, params)

	ACh, tpre = get_matrices(update_ACh!, m, l)
	Φ, synt, δt = get_parameters(update_ACh!, params, l)

	δACh = -ACh ./ synt .+ Φ .* Δ(t .- tpre)

	ACh .= ACh .+ δACh .* δt
end


text:
function jldopen(f::Function, args...; kws...)
    jld = jldopen(args...; kws...)
    try
        return f(jld)
    finally
        close(jld)
    end
end

"""
    @save filename var1 [var2 ...]
    @save filename {compress=true} var1 name2=var2

Write one or more variables `var1,...` from the current scope to a JLD2 file
`filename`.

For interactive use you can save all variables in the current module's global
scope using `@save filename`. More permanent code should prefer the explicit
form to avoid saving unwanted variables.

# Example

To save the string `hello` and array `xs` to the JLD2 file example.jld2:

    hello = "world"
    xs = [1,2,3]
    @save "example.jld2" hello xs

For passing options to the saving command use {}

    @save "example.jld2" {compress=true} hello xs

For saving variables under a different name use regular assignment syntax

    @save "example.jld2" greeting=hello xarray = xs
"""
macro save(filename, vars...)
    fields = []
    options = []
    for var in vars
        # Capture options of the form {compress = true, mmaparrays=false}
        if @capture(var, {opts__})
            for opt in opts
                if @capture(opt, key_ = val_)
                    push!(options, :($key = $(esc(val))))
                else
                    return :(throw(ArgumentError("Invalid option syntax")))
                end
            end
        # Allow assignment syntax a = b
        elseif @capture(var, a_ = b_)
            push!(fields, :(write(f, $(string(a)), $(esc(b)), wsession)))
        # Allow single arg syntax a   → "a" = a
        elseif @capture(var, a_Symbol)
            push!(fields, :(write(f, $(string(a)), $(esc(a)), wsession)))
        else
            return :(throw(ArgumentError("Invalid field syntax")))
        end
    end
    if !isempty(fields)
        return quote
            let
                f = jldopen($(esc(filename)), "w"; $(Expr(:tuple,options...))...)
                wsession = JLDWriteSession()
                try
                    $(Expr(:block, fields...))
                catch e
                    rethrow(e)
                finally
                    close(f)
                end
            end
        end
    else
        # The next part is old code that handles saving the whole workspace
        quote
            let
                m = $(__module__)
                f = jldopen($(esc(filename)), "w"; $(Expr(:tuple,options...))...)
                wsession = JLDWriteSession()
                try
                    for vname in names(m; all=true)
                        s = string(vname)
                        if !occursin(r"^_+[0-9]*$", s) # skip IJulia history vars
                            v = getfield(m, vname)
                            if !isa(v, Module)
                                try
                                    write(f, s, v, wsession)
                                catch e
                                    if isa(e, PointerException)
                                        @warn("skipping $vname because it contains a pointer")
                                    else
                                        rethrow(e)
                                    end
                                end
                            end
                        end
                    end
                finally
                    close(f)
                end
            end
        end
    end
end

"""
    @load filename var1 [var2 ...]

Load one or more variables `var1,...` from JLD2 file `filename` into the
current scope and return a vector of the loaded variable names.

For interactive use, the form `@load "somefile.jld2"` will load all variables
from `"somefile.jld2"` into the current scope. This form only supports literal
file names and should be avoided in more permanent code so that it's clear
where the variables come from.

# Example

To load the variables `hello` and `foo` from the file example.jld2, use

    @load "example.jld2" hello foo
"""
macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            throw(ArgumentError("filename argument must be a string literal unless variable names are specified"))
        end
        # Load all variables in the top level of the file
        readexprs = Expr[]
        vars = Symbol[]
        f = jldopen(filename)
        try
            for n in keys(f)
                if !isgroup(f, lookup_offset(f.root_group, n))
                    push!(vars, Symbol(n))
                end
            end
        finally
            close(f)
        end
    end
    return quote
        ($([esc(x) for x in vars]...),) = jldopen($(esc(filename))) do f
            ($([:(read(f, $(string(x)))) for x in vars]...),)
        end
        $(Symbol[v for v in vars]) # convert to Array
    end
end

function loadtodict!(d::Dict, g::Union{JLDFile,Group}, prefix::String="")
    for k in keys(g)
        v = g[k]
        if v isa Group
            loadtodict!(d, v, prefix*k*"/")
        else
            d[prefix*k] = v
        end
    end
    return d
end


text:
abstract type ProcessingMethod end

include("processing_methods/stopping_criteria/stopping_criteria.jl")
include("processing_methods/stopping_criteria/density_convergence.jl")
include("processing_methods/track_hydrodynamic_errors.jl")
include("processing_methods/visualize.jl")
include("processing_methods/take_snapshots.jl")
include("processing_methods/process_iterative_initialization.jl")

ProcessingMethod(problem, should_process, n_steps, stop_criteria = StopCriteria(problem)) =
    CompareWithAnalyticalSolution(problem, should_process, n_steps, stop_criteria)
ProcessingMethod(
    problem::TaylorGreenVortex,
    should_process,
    n_steps,
    stop_criteria = StopCriteria(problem),
) = TrackHydrodynamicErrors(problem, should_process, n_steps, stop_criteria)
ProcessingMethod(
    problem::DecayingShearFlow,
    should_process,
    n_steps,
    stop_criteria = StopCriteria(problem),
) = TrackHydrodynamicErrors(problem, should_process, n_steps, stop_criteria)
ProcessingMethod(
    problem::TGV,
    should_process,
    n_steps,
    stop_criteria = StopCriteria(problem),
) = TrackHydrodynamicErrors(problem, should_process, n_steps, stop_criteria)
# ProcessingMethod(problem::TGV, should_process, n_steps) =
#     TrackHydrodynamicErrors(problem, should_process, n_steps)

struct CompareWithAnalyticalSolution{T} <: ProcessingMethod
    problem::FluidFlowProblem
    should_process::Bool
    n_steps::Int64
    stop_criteria::StopCriteria
    df::T
end
CompareWithAnalyticalSolution(
    problem,
    should_process,
    n_steps,
    stop_criteria = StopCriteria(problem),
) = CompareWithAnalyticalSolution(
    problem,
    should_process,
    n_steps,
    stop_criteria,
    Vector{
        NamedTuple{
            (
                :density,
                :momentum,
                :total_energy,
                :kinetic_energy,
                :internal_energy,
                :density_a,
                :momentum_a,
                :total_energy_a,
                :kinetic_energy_a,
                :internal_energy_a,
                :error_u,
                :error_p,
                :error_σ_xx,
                :error_σ_xy,
                :error_σ_yy,
                :error_σ_yx,
            ),
            Tuple{
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
                Float64,
            },
        },
    }(),
)

function next!(process_method::CompareWithAnalyticalSolution, q, f_in, t::Int64)
    if mod(t, 100) == 0
        if (should_stop!(process_method.stop_criteria, q, f_in))
            @info "Stopping after $t steps "

            Δt = delta_t(process_method.problem)
            process!(
                process_method.problem,
                q,
                f_in,
                t * Δt,
                process_method.df;
                should_visualize = false,
            )
            return true
        end
    end

    if (!process_method.should_process)
        if (t != process_method.n_steps)
            return false
        end
    end

    if mod(t, 1) == 0
        should_visualize = false
        if (process_method.should_process)
            if t == process_method.n_steps
                should_visualize = true
            end

            if mod(t, max(10, round(Int, process_method.n_steps / 5))) == 0
                should_visualize = true
            end
        end

        Δt = delta_t(process_method.problem)
        process!(
            process_method.problem,
            q,
            f_in,
            t * Δt,
            process_method.df;
            should_visualize = should_visualize,
        )
    end

    return false
end

function process!(
    problem::FluidFlowProblem,
    q::Quadrature,
    f_in,
    time,
    stats;
    should_visualize = false,
)
    f = Array{Float64}(undef, size(f_in, 3))
    u = zeros(dimension(q))
    expected_u = zeros(dimension(q))

    Nx = size(f_in, 1)
    Ny = size(f_in, 2)

    x_range, y_range = range(problem)

    total_density = 0.0
    total_momentum = 0.0
    total_energy = 0.0
    total_kinetic_energy = 0.0
    total_internal_energy = 0.0

    expected_total_density = 0.0
    expected_total_momentum = 0.0
    expected_total_energy = 0.0
    expected_total_kinetic_energy = 0.0
    expected_total_internal_energy = 0.0

    rho_error_squared = 0.0
    ux_error_squared = 0.0
    uy_error_squared = 0.0
    error_u = 0.0
    error_p = 0.0

    @inbounds for x_idx in 1:Nx, y_idx in 1:Ny
        # Calculated
        @inbounds for f_idx in 1:size(f_in, 3)
            f[f_idx] = f_in[x_idx, y_idx, f_idx]
        end
        ρ = density(q, f)
        velocity!(q, f, ρ, u)

        # Adding the forcing term moves the optimal tau for poiseuille flows
        # F = cm.force(x_idx, y_idx, 0.0)
        # u += cm.τ * F

        T = temperature(q, f, ρ, u)
        p = pressure(q, f, ρ, u)

        ρ = dimensionless_density(problem, ρ)
        u = dimensionless_velocity(problem, u)
        T = dimensionless_temperature(q, problem, T)
        p = dimensionless_pressure(q, problem, p)

        opp = 1.0 / (Nx * Ny)
        # CHECK!
        opp = 1.0
        total_density += opp * ρ
        total_momentum += opp * (u[1] + u[2]) * ρ
        kinetic_energy = opp * (u[1]^2 + u[2]^2) * ρ

        internal_energy = T

        total_kinetic_energy += opp * kinetic_energy
        total_internal_energy += opp * internal_energy
        total_energy += opp * (kinetic_energy + internal_energy)

        # Analytical
        x = x_range[x_idx]
        y = y_range[y_idx]

        # Compute statistics
        expected_ρ = LatticeBoltzmann.density(q, problem, x, y, time)
        expected_p = LatticeBoltzmann.pressure(q, problem, x, y, time)
        expected_ϵ = (dimension(q) / 2) * expected_p / expected_ρ
        expected_v = LatticeBoltzmann.velocity(problem, x, y, time)
        expected_T = expected_p / expected_ρ

        expected_kinetic_energy = (expected_v[1]^2 + expected_v[2]^2)
        expected_internal_energy = expected_T

        expected_total_density += expected_ρ
        expected_total_momentum += expected_ρ * (expected_v[1] + expected_v[2])
        expected_total_kinetic_energy += expected_kinetic_energy
        expected_total_internal_energy += expected_internal_energy
        expected_total_energy += expected_kinetic_energy + expected_internal_energy

        rho_error_squared += (ρ - expected_ρ)^2
        ux_error_squared += (u[1] - expected_v[1])^2
        uy_error_squared += (u[2] - expected_v[2])^2
        # error_u += sqrt((u[1] - expected_v[1])^2 + (u[2] - expected_v[2])^2)
        opp = Float64(y_range.step) * Float64(x_range.step)
        # opp = 1.0
        # opp = Float64(y_range.step)
        error_u += (opp * (((u[1] - expected_v[1]))^2 + ((u[2] - expected_v[2]))^2))
        error_p += opp * ((p - expected_p)^2)
    end

    # Compare with analytical results?
    push!(
        stats,
        (
            density = total_density,
            momentum = total_momentum,
            total_energy = total_energy,
            kinetic_energy = total_kinetic_energy,
            internal_energy = total_internal_energy,
            density_a = expected_total_density,
            momentum_a = expected_total_momentum,
            total_energy_a = expected_total_energy,
            kinetic_energy_a = expected_total_kinetic_energy,
            internal_energy_a = expected_total_internal_energy,
            error_u = sqrt(error_u),
            error_p = sqrt(error_p),
            error_σ_xx = 0.0,
            error_σ_xy = 0.0,
            error_σ_yy = 0.0,
            error_σ_yx = 0.0,
        ),
    )

    if should_visualize
        visualize(problem, q, f_in, time, stats)
    end

    return false
end

struct ShowVelocityError{V} <: ProcessingMethod
    problem::FluidFlowProblem
    plot_every::Int
    # TODO: change to vector of vector?
    l2_norms::Vector{V}
    stop_criteria::StopCriteria
end
ShowVelocityError(problem, plot_every = 1, l2_norms = Float64[]) =
    ShowVelocityError(problem, plot_every, l2_norms, NoStoppingCriteria())

function next!(process_method::ShowVelocityError, q, f, t::Int64)
    should_stop = false
    if mod(t, 100) == 0
        if (should_stop!(process_method.stop_criteria, q, f))
            @info "Stopping after $t iterations"
            should_stop = true
        end
    end
    # if mod(t, process_method.plot_every) != 1
    #     return false
    # end

    problem = process_method.problem
    analytical_velocity = (y) -> velocity(problem, 1.0, y)
    x_domain = (0.0, problem.domain_size[1])
    y_domain = (0.0, problem.domain_size[2])

    x_idx = 1

    v_y = zeros(problem.NY)
    v_e = zeros(problem.NY)
    v_a = zeros(problem.NY)
    u = zeros(dimension(q))
    x_range, y_range = range(problem)
    total_expected_momentum = 0.0
    for y_idx in 1:(problem.NY)
        ρ = density(q, f[x_idx, y_idx, :])
        velocity!(q, f[x_idx, y_idx, :], ρ, u)
        u = dimensionless_velocity(problem, u)
        # v_e[y_idx] = norm(u - analytical_velocity(y_range[y_idx]))
        v_e[y_idx] = (
            (u[1] - analytical_velocity(y_range[y_idx])[1])^2 +
                (u[1] - analytical_velocity(y_range[y_idx])[1])^2
        )
        v_y[y_idx] = u[1]
        v_a[y_idx] = analytical_velocity(y_range[y_idx])[1]

        total_expected_momentum +=
            analytical_velocity(y_range[y_idx])[1]^2 +
            analytical_velocity(y_range[y_idx])[2]^2
    end

    push!(process_method.l2_norms, sqrt(sum(v_e) / total_expected_momentum))
    # push!(process_method.l2_norms, v_e)

    # return false
    if mod(t, process_method.plot_every) != 1 && !should_stop
        return should_stop
    end

    y_idx = problem.NY
    ρ = density(q, f[x_idx, y_idx, :])
    velocity!(q, f[x_idx, y_idx, :], ρ, u)
    # @show v_y[problem.NY], v_a[problem.NY], v_e[problem.NY], u[1], ρ

    velocity_plot = plot(
        y_range,
        v_y,
        label = "Numerical solution",
        legend = :topleft,
        title = string(q),
    )
    plot!(velocity_plot, y_range, v_a, label = "Exact solution")

    plot(scatter(y_range, v_e, title = "Abs. Error", legend = nothing), velocity_plot)
    gui()

    return should_stop
end


text:
using ConstrainedDynamics
using ConstrainedDynamicsVis
using Rotations
import MeshCat

boxtemplate = Box(rand(4)...; color = RGBA(rand(4)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
cylindertemplate = Cylinder(rand(3)...; color = RGBA(rand(4)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
spheretemplate = Sphere(rand(2)...; color = RGBA(rand(4)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
pyramidtemplate = Pyramid(rand(3)...; color = RGBA(rand(4)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))
meshtemplate = Mesh("test_files/test.obj", rand(), rand(3,3); color = RGBA(rand(4)...), xoffset = rand(3), qoffset = rand(UnitQuaternion))

originbox = Origin(deepcopy(boxtemplate))
linkbox = deepcopy(boxtemplate)
nonmechbox = deepcopy(boxtemplate)

origincylinder = Origin(deepcopy(cylindertemplate))
linkcylinder = deepcopy(cylindertemplate)
nonmechcylinder = deepcopy(cylindertemplate)

originsphere = Origin(deepcopy(spheretemplate))
linksphere = deepcopy(spheretemplate)
nonmechsphere = deepcopy(spheretemplate)

originpyramid = Origin(deepcopy(pyramidtemplate))
linkpyramid = deepcopy(pyramidtemplate)
nonmechpyramid = deepcopy(pyramidtemplate)

originmesh = Origin(deepcopy(meshtemplate))
linkmesh = deepcopy(meshtemplate)
nonmechmesh = deepcopy(meshtemplate)

mechbox = Mechanism(originbox, [linkbox])
mechcylinder = Mechanism(origincylinder, [linkcylinder])
mechsphere = Mechanism(originsphere, [linksphere])
mechpyramid = Mechanism(originpyramid, [linkpyramid])
mechmesh = Mechanism(originmesh, [linkmesh])

steps = Base.OneTo(10)
storagebox = simulate!(mechbox, 0.1, record = true)
storagecylinder = simulate!(mechcylinder, 0.1, record = true)
storagesphere = simulate!(mechsphere, 0.1, record = true)
storagepyramid = simulate!(mechpyramid, 0.1, record = true)
storagemesh = simulate!(mechmesh, 0.1, record = true)


text:
__precompile__()

module AutoDiff
using DataStructures
export Constant,Variable, Call, CoTangent, Lift, val
export rect_lin,sigmoid,convolve
export minimize

import Base: (+), (-), (*), (/), (^),
(.+), (.-), (.*), (./), (.^), exp, log, sum

begin
	#my_zero(x) returns an array of zeros of the same shape as x
	#x can be a number, an array, an array of arrays, etc.
	function my_zero(A::Array)
		[my_zero(x) for x in A]
	end
	function my_zero{T<:Number}(A::Array{T})
		zeros(A)
	end
	function my_zero(x)
		zero(x)
	end
end

#Type CoTangent{T} stores a value of type T and a derivative
abstract CoTangent{T}
type Variable{T} <: CoTangent{T}
	val::T
	delta::T
end
type Constant{T} <: CoTangent{T}
	val::T
	delta::T
end
Variable{T}(val::T)=Variable{T}(val,my_zero(val))
Constant{T}(val::T)=Constant{T}(val,my_zero(val))

is_CoTangent(x::CoTangent)=true
is_CoTangent(x)=false

type Call{F,T} <: CoTangent{T}
	f::Function
	args::Tuple
	val::T
	delta::T
	function Call(f,args,v,delta)
		new(f,args,v,delta)
	end
end

type Lift{T}<:CoTangent{T}
	args
	val::T
	delta::T
end
function Lift(args)
	tmp=map(val,args)
	Lift(args,tmp,my_zero(tmp))
end

val(x)=x
val(x::CoTangent)=x.val

Call(f,args)=Call(f,args,f(map(val,args)...))
Call{T}(f::Function,args,val::T)=Call{Base.function_name(f),T}(f,args,val,my_zero(val))

begin
	function execute(x::Call)
		x.val=x.f(map(val,x.args)...)
	end
	function execute(x::Lift)
		x.val=map(val,x.args)
	end
	function execute(x) end
end

begin
	function clear_array{T<:Real}(A::Array{T})
		fill!(A,zero(T))
	end
	function clear_array{T<:Array}(A::Array{T})
		map(clear_array,A)
	end
	function reset_delta{T<:Real}(t::CoTangent{T})
		t.delta=zero(T)
	end
	function reset_delta{T<:Array}(t::CoTangent{T})
		clear_array(t.delta)
	end
end

function backprop end

#returns a function which computes f'
#expects a function f which takes an array of vars
#and an array of constants
function derivative(f, vars, constants)
	wrapped_vars=map(x->Variable(x),vars)
	wrapped_constants=map(x->Constant(x),constants)
	output=f(wrapped_vars,wrapped_constants)
	l=topological_order(output)
	function update_params(new_params)
		for (x,y) in zip(wrapped_vars,new_params)
			x.val=y
		end
	end
	function gradient(constants)
		map(reset_delta,l)
		for (x,y) in zip(wrapped_constants,constants)
			x.val=y
		end
		map(execute,l)
		output.delta=1
		map(backprop,reverse(l))
		return map(x->x.delta,wrapped_vars)
	end

	return (update_params,gradient)
end

function compute_grad(expr,l)
	map(reset_delta,l)
	expr.delta=1
	map(backprop,reverse(l))
end
begin
	function descend(d::Variable,delta)
		d.val -= delta.*d.delta
	end
	function descend(d::CoTangent,delta) end
end
function minimize(c::CoTangent;delta=0.01,its=1000,f=(it->nothing))
	l=topological_order(c)
	for i in 1:its
		map(execute,l)
		f(i)
		compute_grad(c,l)
		map(x->descend(x,delta), l)
	end
end

function depth_first_map(f,t::CoTangent)
	visited=Set()
	function traverse(c::Union{Call,Lift})
		if !(c in visited)
			for child in c.args
				traverse(child)
			end
			f(c)
			push!(visited,c)
		end
	end
	function traverse(c::Variable)
		if !(c in visited)
			f(c)
			push!(visited,c)
		end
	end
	function traverse(c) end
	traverse(t)
end

function topological_order(c::CoTangent)
	q=Any[]
	depth_first_map(x->push!(q,x),c)
	q
end

include("AutoDiff/printing.jl")
include("AutoDiff/backprop.jl")
include("AutoDiff/call_defs.jl")

function sigmoid(x)
	1 ./(1 .+exp(-x))
end
function rect_lin(x)
	max(x,0)
end
end


text:
function fib(i::Int64)::Int64
    if i == 0 || i == 1
        return 1
    end
    return (fib((i - 1)) + fib((i - 2)))
end

function main()
    println(join([fib(5)], " "))
end

main()


text:
#  MIPLearn: Extensible Framework for Learning-Enhanced Mixed-Integer Optimization
#  Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
#  Released under the modified BSD license. See COPYING.md for more details.

using MIPLearn

@testset "Parse" begin
    @test MIPLearn.parse_name("x") == ["x"]
    @test MIPLearn.parse_name("x[3]") == ["x", "3"]
    @test MIPLearn.parse_name("test_eq[x]") == ["test_eq", "x"]
    @test MIPLearn.parse_name("test_eq[x,y,z]") == ["test_eq", "x", "y", "z"]
end


text:
module Session

import HPack
import HPack: DynamicTable
import HTTP2: bytearr, Frame, Headers, Maybe
import HTTP2.Frame: ContinuationFrame, DataFrame, GoawayFrame, HeadersFrame, PingFrame,
                    PriorityFrame, PushPromiseFrame, RstStreamFrame, SettingsFrame,
                    WindowUpdateFrame

@enum STREAM_STATE IDLE=1 RESERVED_LOCAL RESERVED_REMOTE OPEN HALF_CLOSED_REMOTE HALF_CLOSED_LOCAL CLOSED

mutable struct Priority
    dependent_stream_identifier::UInt32
    weight::UInt8
end

## Actions, which should be feeded in to `in` channel.

struct ActPromise
    stream_identifier::UInt32
    promised_stream_identifier::UInt32
    headers::Headers
end

struct ActSendHeaders
    stream_identifier::UInt32
    headers::Headers
    is_end_stream::Bool
end

struct ActSendData
    stream_identifier::UInt32
    data::Vector{UInt8}
    is_end_stream::Bool
end

## Events, which should be fetched from `out` channel.

struct EvtPromise
    stream_identifier::UInt32
    promised_stream_identifier::UInt32
    headers::Headers
end

struct EvtRecvHeaders
    stream_identifier::UInt32
    headers::Headers
    is_end_stream::Bool
end

struct EvtRecvData
    stream_identifier::UInt32
    data::Vector{UInt8}
    is_end_stream::Bool
end

struct EvtGoaway end

mutable struct HTTPStream
    stream_identifier::UInt32
    state::STREAM_STATE
    window_size::UInt32
    priority::Maybe{Priority}
end

mutable struct HTTPSettings
    push_enabled::Bool
    max_concurrent_streams::Maybe{UInt}
    initial_window_size::UInt
    max_frame_size::UInt
    max_header_list_size::Maybe{UInt}
end

HTTPSettings() = HTTPSettings(true, nothing, 65535, 16384, nothing)

const DEFAULT_CHANNEL_SZ = 1024
mutable struct HTTPConnection
    dynamic_table::DynamicTable
    streams::Vector{HTTPStream}
    window_size::UInt32
    isclient::Bool
    last_stream_identifier::UInt32
    settings::HTTPSettings
    closed::Bool

    channel_act::Channel{Any} # Process actions
    channel_act_raw::Channel{Any} # Process raw frames
    channel_evt::Channel{Any} # Output events
    channel_evt_raw::Channel{Any} # Output raw frames

    ## actions -> channel_act -> channel_act_raw -> io
    ## io -> channel_evt_raw -> channel_evt -> events
    act_processor::Maybe{Task}
    act_raw_processor::Maybe{Task}
    evt_raw_processor::Maybe{Task}
    evt_processor::Maybe{Task}

    function HTTPConnection(isclient)
        new(HPack.DynamicTable(),
            Vector{HTTPStream}(),
            1024,
            isclient,
            isclient ? 1 : 2,
            HTTPSettings(),
            false,
            Channel(DEFAULT_CHANNEL_SZ), Channel(DEFAULT_CHANNEL_SZ),
            Channel(DEFAULT_CHANNEL_SZ), Channel(DEFAULT_CHANNEL_SZ),
            nothing, nothing, nothing, nothing)
    end
end

function next_free_stream_identifier(connection::HTTPConnection)
    return connection.last_stream_identifier + 2
end

include("Session/utils.jl")
include("Session/settings.jl")
include("Session/errors.jl")
include("Session/states.jl")
include("Session/handlers.jl")
include("Session/channels.jl")

function new_connection(buffer; isclient::Bool=true, skip_preface=false)
    connection = HTTPConnection(isclient)
    initialize_loop_async(connection, buffer; skip_preface=skip_preface)
    return connection
end

function put_act!(connection::HTTPConnection, act)
    put!(connection.channel_act, act)
end

function take_evt!(connection::HTTPConnection)
    return take!(connection.channel_evt)
end

end # module Session


text:
using Docker
using Test

# Make sure "hildebrandmw/hello-world" is not in the local list of images. If
# so, remove it.
images = Docker.list_images()
test_image() = "hildebrandmw/hello-world:latest"

# Need to handle the case where there are no repo tags and `image["RepoTags"] === nothing`
function contains_test_image(image)
    tag = image["RepoTags"]
    tag === nothing && return false
    return any(x -> occursin(test_image(), x), tag)
end

if any(contains_test_image, images)
    @info """
    Found local copy of the test repo: $(test_image())
    Removing for test validation
    """
    Docker.remove_image(test_image(); force = true)
end

getids(x) = Docker.getid.(x)
@testset "Testing Functions" begin
    Docker.pull_image(test_image())

    # Make sure that the image we just pulled show up in the list of images.
    images = Docker.list_images()
    @test any(contains_test_image, images) == true

    # Create a container based on this image.
    container = Docker.create_container(test_image())
    id = Docker.getid(container)
    @info "Created Container: $container"
    Docker.start(container)
    # Container should quite fairly quickly after running
    sleep(5)

    # Get the expected response for the log.
    # Wrap the string returned by "get_log" in an IOBuffer and eachline to
    # avoid issues with newlines and carriage returns
    expected_log = eachline("log.txt")
    log = eachline(IOBuffer(Docker.log(container)))
    for (e,l) in zip(expected_log, log)
        @test e == l
    end

    # Test that the contain is in the list of containers
    containers = Docker.list_containers(all = true)
    @test length(containers) > 0
    @test in(id, getids(containers))

    # Delete the container
    Docker.remove(container)
    containers = Docker.list_containers(all = true)
    @test !in(id, getids(containers))

    # Remove the test image
    Docker.remove_image(test_image(); force = true)
end


text:
using Combinatorics
using SymEngine

abstract type AnalyticField <: ElectromagneticField end
abstract type AnalyticEquilibrium <: AnalyticField end
abstract type AnalyticPerturbation <: AnalyticField end

function get_functions end
function get_parameters end

x¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x¹() not implemented for ", ET)
x²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x²() not implemented for ", ET)
x³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("x³() not implemented for ", ET)

ξ¹(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ¹() not implemented for ", ET)
ξ²(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ²() not implemented for ", ET)
ξ³(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("ξ³() not implemented for ", ET)

J(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("J() not implemented for ", ET)

g₁₁(::AbstractVector{T}, ::AnalyticField) where {T} = one(T)
g₁₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₁₃(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₂₁(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₂₂(::AbstractVector{T}, ::AnalyticField) where {T} = one(T)
g₂₃(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₃₁(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₃₂(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)
g₃₃(::AbstractVector{T}, ::AnalyticField) where {T} = one(T)

periodicity(x::AbstractVector, ::AnalyticField) = zero(x)

from_cartesian(x::AbstractVector, equ::AnalyticField) = [ξ¹(x,equ), ξ²(x,equ), ξ³(x,equ)]
to_cartesian(ξ::AbstractVector, equ::AnalyticField) = [x¹(ξ,equ), x²(ξ,equ), x³(ξ,equ)]

A₁(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₁() not implemented for ", ET)
A₂(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₂() not implemented for ", ET)
A₃(::AbstractVector, ::ET) where {ET <: AnalyticField} = error("A₃() not implemented for ", ET)

φ(::AbstractVector{T}, ::AnalyticField) where {T} = zero(T)

function A(x, equ)
    [A₁(x,equ), A₂(x,equ), A₃(x,equ)]
end

function g(x, equ)
    [g₁₁(x,equ)  g₁₂(x,equ)  g₁₃(x,equ);
     g₂₁(x,equ)  g₂₂(x,equ)  g₂₃(x,equ);
     g₃₁(x,equ)  g₃₂(x,equ)  g₃₃(x,equ);]
end

struct ZeroPerturbation <: AnalyticPerturbation
    name::String
    ZeroPerturbation() = new("ZeroPerturbation")
end

A₁(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)
A₂(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)
A₃(::AbstractVector{T}, ::AnalyticPerturbation) where {T} = zero(T)

"Returns the i-th component of the vector corresponding to the one-form α"
function covariant_to_contravariant(α, g̅, i)
    g̅[i,1] * α[1] + g̅[i,2] * α[2] + g̅[i,3] * α[3]
end

"Returns the i-th component of the one-form corresponding to the vector v"
function contravariant_to_covariant(v, g, i)
    g[i,1] * v[1] + g[i,2] * v[2] + g[i,3] * v[3]
end

"Returns the i-th component of the physical coordinate representation of the one-form α"
function covariant_to_physical(α, DF̄, i)
    DF̄[1,i] * α[1] + DF̄[2,i] * α[2] + DF̄[3,i] * α[3]
end

"Returns the i-th component of the physical coordinate representation of the one-form α"
function contravariant_to_physical(v, DF, i)
    DF[i,1] * v[1] + DF[i,2] * v[2] + DF[i,3] * v[3]
end

"Returns the m-th component of the one-form corresponding to the two-form β"
function hodge²¹(β, g̅, J, m)
    α = 0

    for i in 1:3
        for j in 1:3
            for k in 1:3
                for l in 1:3
                    α += β[i,j] * g̅[i,k] * g̅[j,l] * levicivita([k,l,m])
                end
            end
        end
    end

    return J*α
end

"Returns the Christoffel symbol Γⱼₖˡ"
function Γ(g, g̅, x, j, k, l)
    γ = Basic(0)
    for r in 1:3
        γ += g̅[l,r] * diff(g[r,j], x[k])
        γ += g̅[l,r] * diff(g[r,k], x[j])
        γ -= g̅[l,r] * diff(g[j,k], x[r])
    end
    return 1 // 2 * γ
end

"Returns the l-th component of the Levi-Civita connection"
function connection(u, v, x, g, g̅, l)
    w = Basic(0)
    for j in 1:3
        for k in 1:3
            w += u[j] * ( diff(v[l], x[j]) + v[k] * Γ(g, g̅, x, j, k, l) )
        end
    end
    return w
end

"Returns the m-th component of the cross-product between the vectors v and w"
function crossproduct(v, w, g̅, J, l)
    u = zero(J)

    for i in 1:3
        for j in 1:3
            for k in 1:3
                u += v[i] * w[j] * g̅[k,l] * levicivita([i,j,k])
            end
        end
    end

    return J*u
end

"Returns the length of the vector v"
function magnitude(v, g)
    l = Basic(0)

    for i in 1:3
        for j in 1:3
            l += v[i] * g[i,j] * v[j]
        end
    end

    return sqrt(l)
end

"Normalises the vector v in the metric g"
function normalize(v, g)
    return v ./ magnitude(v, g)
end

"Normalises the vector v in the metric g"
function normalize!(v, g)
    v ./= magnitude(v, g)
end

"""
Generate functions for evaluating analytic equilibria.
"""
function generate_equilibrium_functions(equ::AnalyticEquilibrium, pert::AnalyticPerturbation; output=0)
    # define symbols for time t and coordinates x = (x₁, x₂, x₃),
    # positive=true is set so that sqrt(x^2) does not become |x^2|
    t, x₁, x₂, x₃, ξ₁, ξ₂, ξ₃ = symbols("t, x₁, x₂, x₃, ξ₁, ξ₂, ξ₃")#, real=true, positive=true)
    x = [x₁, x₂, x₃]
    ξ = [ξ₁, ξ₂, ξ₃]
    symprint("x", x, output, 2)
    symprint("ξ", ξ, output, 2)

    # check for compatible metric
    if typeof(pert) != ZeroPerturbation
        @assert J(x,equ) == J(x,pert)
        @assert g(x,equ) == g(x,pert)
    end

    # cartesian coordinates
    x̂ = [x¹(ξ, equ), x²(ξ, equ), x³(ξ, equ)]

    # curvilinear coordinates
    ξ̂ = [ξ¹(x, equ), ξ²(x, equ), ξ³(x, equ)]

    # Jacobian
    DF = [diff(x̂[i], ξ[j]) for i in 1:3, j in 1:3]
    symprint("DF", DF, output, 2)

    DF̄ = [subs(diff(ξ̂[i], x[j]), x₁=>x¹(ξ, equ), x₂=>x²(ξ, equ), x₃=>x³(ξ, equ)) for i in 1:3, j in 1:3]
    symprint("DF̄", DF̄, output, 2)

    # obtain metric
    gmat = g(ξ, equ)
    symprint("g", gmat, output, 2)

    # invert metric
    ginv = inv(gmat)
    symprint("g⁻¹", ginv, output, 2)

    # derivatives of metric coefficients
    Dg = [diff(gmat[i,j], ξ[k]) for i in 1:3, j in 1:3, k in 1:3]
    symprint("Dg", Dg, output, 3)

    Dḡ = [diff(ginv[i,j], ξ[k]) for i in 1:3, j in 1:3, k in 1:3]
    symprint("Dḡ", Dḡ, output, 3)

    # compute Jacobian determinant
    # Jdet² = expand(det(gmat))
    # symprint("J²", Jdet², output, 2)

    # Jdet = sqrt(Jdet²)
    Jdet = J(ξ, equ)
    symprint("J", Jdet, output, 2)

    # obtain vector potential
    A¹ = A(ξ, equ) .+ A(ξ, pert)
    symprint("A¹", A¹, output, 2)

    # compute vector potential in contravariant coordinates
    Avec = [covariant_to_contravariant(A¹, ginv, i) for i in 1:3]
    symprint("Avec", Avec, output, 2)

    # compute Jacobian of vector potential A
    DA = [(diff(A¹[i], ξ[j])) for i in 1:3, j in 1:3]
    symprint("DA", DA, output, 3)

    # compute second derivative of vector potential A
    DDA = [diff(DA[i,j], ξ[k]) for i in 1:3, j in 1:3, k in 1:3]
    symprint("DDA", DDA, output, 3)

    # compute components of magnetic field B
    Bᶜ = [DA[3,2] - DA[2,3],
          DA[1,3] - DA[3,1],
          DA[2,1] - DA[1,2]]
    symprint("Bᶜ", Bᶜ, output, 2)

    # compute magnetic field two-form B²
    B² = [ 0      +Bᶜ[3]    -Bᶜ[2];
          -Bᶜ[3]   0        +Bᶜ[1];
          +Bᶜ[2]   -Bᶜ[1]    0    ] .* 1//2
    symprint("B²", B², output, 2)

    # compute magnetic field one-form B¹ = ⋆B²
    B¹ = [hodge²¹(B², ginv, Jdet, i) for i in 1:3]
    symprint("B¹", B¹, output, 2)

    # compute magnetic field in physical coordinates
    Bphys = [covariant_to_physical(B¹, DF̄, i) for i in 1:3]
    symprint("B̂", Bphys, output, 2)

    # compute magnetic field in contravariant coordinates
    Bvec = [covariant_to_contravariant(B¹, ginv, i) for i in 1:3]
    symprint("B⃗", Bvec, output, 2)

    # compute absolute value |B| of B
    Babs = sqrt(transpose(Bvec) * B¹)
    symprint("|B|", Babs, output, 2)

    # compute magnetic unit one-form
    b¹ = [B¹[i] / Babs for i in 1:3]
    symprint("b¹", b¹, output, 2)

    # compute unit magnetic field in physical coordinates
    bphys = [Bphys[i] / Babs for i in 1:3]
    symprint("b̂", bphys, output, 2)

    # compute unit magnetic field in contravariant coordinates
    bvec = [Bvec[i] / Babs for i in 1:3]
    symprint("b⃗", bvec, output, 2)

    # compute Jacobian of magnetic field B
    DB = [diff(B¹[i], ξ[j]) for i in 1:3, j in 1:3]
    symprint("DB", DB, output, 3)

    # compute Jacobian of magnetic unit oneform b
    Db = [diff(b¹[i], ξ[j]) for i in 1:3, j in 1:3]
    symprint("Db", Db, output, 3)

    # compute second derivative of magnetic unit vector b
    DDb = [diff(diff(b¹[i], ξ[j]), ξ[k]) for i in 1:3, j in 1:3, k in 1:3]
    symprint("DDb", DDb, output, 3)

    # compute first derivatives of absolute value of magnetic field
    DBabs = [diff(Babs, ξ[j]) for j in 1:3]
    symprint("D|B|", DBabs, output, 3)

    # compute second derivatives of absolute value of magnetic field
    DDBabs = [diff(diff(Babs, ξ[i]), ξ[j]) for i in 1:3, j in 1:3]
    symprint("DD|B|", DDBabs, output, 3)

    # compute unit vectors perpendicular to magnetic field
    avec = [Basic(0), Basic(0), Basic(0)]
    for tvec ∈ ([Basic(1), Basic(0), Basic(0)],
                [Basic(0), Basic(1), Basic(0)],
                [Basic(0), Basic(0), Basic(1)])
        avec .= [crossproduct(tvec, bvec, ginv, Jdet, i) for i in 1:3]
        if avec != [Basic(0), Basic(0), Basic(0)]
            break
        end
    end
    cvec = [crossproduct(bvec, avec, ginv, Jdet, i) for i in 1:3]

    normalize!(avec, gmat)
    normalize!(cvec, gmat)

    # compute components of magnetic unit vectors in physical coordinates
    aphys = DF * avec
    cphys = DF * cvec

    # compute components of magnetic unit vectors in covariant coordinates
    a¹ = gmat * avec
    c¹ = gmat * cvec

    # obtain scalar potential
    φ⁰ = φ(ξ, equ) .+ φ(ξ, pert)

    # compute components of electric field E
    E¹ = [diff(-φ⁰, ξ[i]) for i in 1:3]
    symprint("E¹", E¹, output, 2)

    # compute electric field in contravariant coordinates
    Evec = [covariant_to_contravariant(E¹, ginv, i) for i in 1:3]
    symprint("Evec", Evec, output, 2)

    # collect all functions to generate code for
    functions = Dict{String,Any}()
    indices   = ["₁", "₂", "₃"]
    indicesup = ["¹", "²", "³"]
    indicesph = ["₍₁₎", "₍₂₎", "₍₃₎"]

    try
        for f in pairs(get_functions(equ))
            functions[string(f[1])] = f[2](ξ, equ)
        end
    catch
    end

    # cartesian coordinates
    functions["x¹"] = x̂[1]
    functions["x²"] = x̂[2]
    functions["x³"] = x̂[3]

    # curvilinear coordinates (replacing x with ξ is needed for the code generator)
    functions["ξ¹"] = subs(ξ̂[1], x₁=>ξ₁, x₂=>ξ₂, x₃=>ξ₃)
    functions["ξ²"] = subs(ξ̂[2], x₁=>ξ₁, x₂=>ξ₂, x₃=>ξ₃)
    functions["ξ³"] = subs(ξ̂[3], x₁=>ξ₁, x₂=>ξ₂, x₃=>ξ₃)

    functions["periodicity"] = periodicity([ξ₁, ξ₂, ξ₃], equ)

    functions["J"] = Jdet
    functions["B"] = Babs
    functions["φ"] = φ⁰

    for i in 1:3
        functions["A" * indices[i]] = A¹[i]
        functions["B" * indices[i]] = B¹[i]
        functions["a" * indices[i]] = a¹[i]
        functions["b" * indices[i]] = b¹[i]
        functions["c" * indices[i]] = c¹[i]
        functions["E" * indices[i]] = E¹[i]

        functions["B" * indicesph[i]] = Bphys[i]
        functions["a" * indicesph[i]] = aphys[i]
        functions["b" * indicesph[i]] = bphys[i]
        functions["c" * indicesph[i]] = cphys[i]

        functions["A" * indicesup[i]] = Avec[i]
        functions["B" * indicesup[i]] = Bvec[i]
        functions["a" * indicesup[i]] = avec[i]
        functions["b" * indicesup[i]] = bvec[i]
        functions["c" * indicesup[i]] = cvec[i]
        functions["E" * indicesup[i]] = Evec[i]

        functions["dBdx" * indices[i]] = DBabs[i]
    end

    for i in 1:3
        for j in 1:3
            functions["g"  * indices[i]   * indices[j]]   = gmat[i,j]
            functions["g"  * indicesup[i] * indicesup[j]] = ginv[i,j]

            functions["DF" * indices[i]   * indices[j]]   = DF[i,j]
            functions["DF̄" * indices[i]   * indices[j]]   = DF̄[i,j]

            functions["B"  * indices[i]   * indices[j]]   = B²[i,j]

            functions["dA" * indices[i]   * "dx" * indices[j]] = DA[i,j]
            functions["dB" * indices[i]   * "dx" * indices[j]] = DB[i,j]
            functions["db" * indices[i]   * "dx" * indices[j]] = Db[i,j]

            functions["d²B"  * "dx" * indices[i] * "dx" * indices[j]] = DDBabs[i,j]
        end
    end

    for i in 1:3
        for j in 1:3
            for k in 1:3
                functions["d²A" * indices[i] * "dx" * indices[j] * "dx" * indices[k]] = DDA[i,j,k]
                functions["d²b" * indices[i] * "dx" * indices[j] * "dx" * indices[k]] = DDb[i,j,k]
                functions["dg" * indices[i]   * indices[j]   * "dx" * indices[k]] = Dg[i,j,k]
                functions["dg" * indicesup[i] * indicesup[j] * "dx" * indices[k]] = Dḡ[i,j,k]
            end
        end
    end

    functions
end

function replace_expr!(e, old, new)
    for (i,a) in enumerate(e.args)
        if a==old
            e.args[i] = new
        elseif a isa Expr
            replace_expr!(a, old, new)
        end
    end
    e
end

fnesc(name, escape) = escape ? esc(name) : name

"""
Generate code for evaluating analytic equilibria.
"""
function code(equ, pert=ZeroPerturbation(); export_parameters=true, escape=false, output=0)

    if output ≥ 1
        println("Generating code for ")
        println(equ)
        if typeof(pert) != ZeroPerturbation
            println("   and ")
            println(pert)
        end
        println()
    end

    functions = generate_equilibrium_functions(equ, pert; output=output)

    equ_code = quote end

    # generate Julia code and export parameters
    if export_parameters
        try
            global parameters = get_parameters(equ)
        catch
            global parameters = fieldnames(typeof(equ))
        end

        for param in parameters
            if param != :name
                value = getfield(equ, param)

                p_code = quote
                    export $param
                    $(fnesc(param, escape)) = $value
                end

                # append p_code to equ_code
                append!(equ_code.args, p_code.args)
            end
        end
    end

    # add wrapper functions
    functions["a"] = quote
       [$(fnesc(:a₁, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:a₂, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:a₃, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["b"] = quote
       [$(fnesc(:b₁, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:b₂, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:b₃, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["c"] = quote
       [$(fnesc(:c₁, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:c₂, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:c₃, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["aₚ"] = quote
       [$(fnesc(:a₍₁₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:a₍₂₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:a₍₃₎, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["bₚ"] = quote
       [$(fnesc(:b₍₁₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:b₍₂₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:b₍₃₎, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["cₚ"] = quote
       [$(fnesc(:c₍₁₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:c₍₂₎, escape))(t, ξ₁, ξ₂, ξ₃),
        $(fnesc(:c₍₃₎, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["a⃗"] = quote
        [$(fnesc(:a¹, escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:a², escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:a³, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["b⃗"] = quote
        [$(fnesc(:b¹, escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:b², escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:b³, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["c⃗"] = quote
        [$(fnesc(:c¹, escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:c², escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:c³, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["from_cartesian"] = quote
        [$(fnesc(:ξ¹, escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:ξ², escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:ξ³, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["to_cartesian"] = quote
        [$(fnesc(:x¹, escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:x², escape))(t, ξ₁, ξ₂, ξ₃),
         $(fnesc(:x³, escape))(t, ξ₁, ξ₂, ξ₃)]
    end

    functions["DF"] = quote
        [$(fnesc(:DF₁₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₁₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₁₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:DF₂₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₂₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₂₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:DF₃₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₃₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF₃₃, escape))(t, ξ₁, ξ₂, ξ₃);]
    end

    functions["DF̄"] = quote
        [$(fnesc(:DF̄₁₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₁₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₁₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:DF̄₂₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₂₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₂₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:DF̄₃₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₃₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:DF̄₃₃, escape))(t, ξ₁, ξ₂, ξ₃);]
    end

    functions["g"] = quote
        [$(fnesc(:g₁₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₁₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₁₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:g₂₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₂₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₂₃, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:g₃₁, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₃₂, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g₃₃, escape))(t, ξ₁, ξ₂, ξ₃);]
    end

    functions["ḡ"] = quote
        [$(fnesc(:g¹¹, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g¹², escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g¹³, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:g²¹, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g²², escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g²³, escape))(t, ξ₁, ξ₂, ξ₃);
         $(fnesc(:g³¹, escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g³², escape))(t, ξ₁, ξ₂, ξ₃)  $(fnesc(:g³³, escape))(t, ξ₁, ξ₂, ξ₃);]
    end

    # generate Julia code and export functions
    for (key, value) in functions
        f_symb = fnesc(Symbol(key), escape)
        f_expr = value

        output ≥ 1 ? println("Generating function ", key) : nothing

        f_body = convert(Expr, f_expr)
        replace_expr!(f_body, :atan2, :atan)
        output ≥ 2 ? println("   ", f_body) : nothing

        f_code = quote
            export $f_symb
            function $f_symb(t, ξ₁, ξ₂, ξ₃)
                $f_body
            end
            function $f_symb(t, ξ)
                $f_symb(t,ξ[1],ξ[2],ξ[3])
            end
        end

        # append f_code to equ_code
        append!(equ_code.args, f_code.args)
    end

    # generate Julia code and export wrapper functions
    f_code = quote
        $(fnesc(:periodicity, escape))(ξ) = $(fnesc(:periodicity, escape))(0, ξ)
        $(fnesc(:periodicity, escape))(ξ₁, ξ₂, ξ₃) = $(fnesc(:periodicity, escape))(0, ξ₁, ξ₂, ξ₃)
    end

    # append f_code to equ_code
    append!(equ_code.args, f_code.args)

    output ≥ 1 ? println() : nothing

    return equ_code
end

"""
Evaluate functions for evaluating analytic equilibria.
"""
function load_equilibrium(equ, pert=ZeroPerturbation(); target_module=Main, output=0)
    equ_code = code(equ, pert; output=output)
    Core.eval(target_module, equ_code)
end

function symprint(name, symexpr, output=1, detail_level=0)
    if output ≥ detail_level
        println(name, " = ", symexpr, "\n")
    end
end


text:
export GaussCopula, ProfitLoss, PLInsurance, PLInvestments,
PLTotal, BusinessUnit, BuInsurance, BuInvestments, Total

"Gaussian copula"
type GaussCopula
  "`Int`: Number of marginal distributions"
  n::Int
  "`Array{Float64,2} `: Correlation matrix"
  Σ::Array{Float64,2}
  "`Vector{ContinuousUnivariateDistribution}`:
  Marginal distributions"
  marginals::Vector{ContinuousUnivariateDistribution}
end

"Profit and loss account of an insurance company component"
abstract ProfitLoss

"Profit loss account for the insurance result from a line
of business (lob)"
type PLInsurance <: ProfitLoss
  "`Real`: Premium for lob"
  premium::Real
  "`Real`: Administration costs of lob"
  costs::Real
  "`Real`: Fraction of premium which is quota reinsured"
  ceded::Real
  "`Vector{Real}`: Insurance profit from lob"
  profit::Vector{Real}
  "`Real`: Expected insurance profit from lob"
  profit_mean::Real
  "`Real`: Economic capital for lob"
  eco_cap::Real
  "`Real`: Return on risk adjusted capital for lob"
  rorac::Real
end

"Profit loss account for the investment result result"
type PLInvestments <: ProfitLoss
  "`Real`: Invested funds beginning of period"
  invest_bop::Real
  "`Real`: Administration costs"
  costs::Real
  "`Vector{Real}`: Investment profit"
  profit::Vector{Real}
  "`Real`: Expected investment profit"
  profit_mean::Real
  "`Real`: Economic capital for asset management"
  eco_cap::Real
  "`Real`: Return on risk adjusted capital for asset management"
  rorac::Real
end

"Profit loss account for the aggregated insurance company"
type PLTotal <: ProfitLoss
  "`Vector{Real}`: Profit of company"
  profit::Vector{Real}
  "Real`: Expected profit of company"
  profit_mean::Real
  "Economic capital for the company"
  eco_cap::Real
  "`Real`: Return on risk adjusted capital for the company"
  rorac::Real
end

" A business unit is a line of business or investments"
abstract BusinessUnit

"Business unit: Line of business"
type BuInsurance <: BusinessUnit
  "`AbstractString`:  id of business unit"
  id::Symbol
  "`AbstractString`:  Name of business unit"
  name::AbstractString
  "`PLInsurance`: profit loss account for lob
  (gross of reinsurance)"
  gross::PLInsurance
  "`PLInsurance`: profit loss account for lob
  (net of reinsurance)"
  net::PLInsurance
end

"Business unit: Investments"
type BuInvestments <: BusinessUnit
  "`AbstractString`:  id of business unit"
  id::Symbol
  "`AbstractString`:  Name of business unit"
  name::AbstractString
  "`Real`:  Initial capital"
  init::Real
  "`PLInvestments`: profit loss account for investments
  (gross of reinsurance)"
  gross::PLInvestments
  "`PLInvestments`: profit loss account for investments
  (net of reinsurance)"
  net::PLInvestments
end

"Total company result"
type Total
  "`PLTotal`: profit loss account for company
  (gross of reinsurance)"
  gross::PLTotal
  "`PLTotal`: profit loss account for company
  (net of reinsurance)"
  net::PLTotal
end


text:
function matching_brackets(input)
    input = split(input, "")
    open_brackets,close_brackets = ["(", "[", "{"], [")", "]", "}"]
    input = filter(i -> (i in union(open_brackets, close_brackets)), input)
    stack = Array{String}(undef, 0)
    brackets = Dict([("(", ")"),(")", "("), ("[", "]"), ("]", "["), ("{", "}"), ("}", "{")]) #global
    count, valid = 0,1 #global
    length(input) % 2 != 0 && return false
    for i in input
        if i in open_brackets
            push!(stack, i)
            count += 1 #global
        elseif count == 0 # & downward: assumes i not in open_brackets
            valid = 0
            break
        elseif stack[end] != brackets[i] #& downward: assumes count != 0
            valid = 0
            break
        else # assumes stack[end] == d[i]
            pop!(stack)
            count -= 1
        end
    end
    return !(valid == 0 || count != 0)
end


text:
using Gadfly, Distributions

function plot2D(G::Array{ConjugatePostDistribution};
    clim::(Float64, Float64, Float64, Float64) = (-10.0, 10.0, -10.0, 10.0), res::Int=50)

    sx = (clim[2]-clim[1])/(res-1)
    sy = (clim[4]-clim[3])/(res-1)
    A = Array(Float64,2,res^2)
    for (i,a) in enumerate(clim[1]:sx:clim[2])
        for (j,b) in enumerate(clim[3]:sy:clim[4])
            A[:,(i-1)*res+j] = [a,b]
        end
    end

    p = zeros(res * res)
    for d in G

        μ = vec(d.sums / d.n)
        if μ != μ
            μ = vec(zeros(size(d.sums)))
        end

        Σ = (d.ssums / d.n) + (μ * μ')

        p += pdf(MvNormal(μ, Σ), A)
    end

    z= reshape(p,res,res)
    return Gadfly.plot(z=z,x=[clim[1]:sx:clim[2]],y=[clim[3]:sy:clim[4]],Gadfly.Geom.contour)
end


text:
using DrWatson
@quickactivate "StatReth"

# %%
using DataFrames
using CSV
using Distributions
using Turing
using StatsPlots

include(srcdir("quap.jl"))
include(srcdir("tools.jl"))

# %% 4.26
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
d2 = d[d.age .>= 18, :]

# %% 4.37
scatter(d2.weight, d2.height)

# %% 4.38
N = 100
a = rand(Normal(178.0, 20.0), N)
b = rand(Normal(0.0, 10.0), N)

# %% 4.39
plot(xlims = extrema(d2.weight), ylims = (-100, 400), xlable = "weight", ylabel = "height")
hline!([0.0, 272.0])
x̄ = mean(d2.weight)
x = range(minimum(d2.weight), maximum(d2.weight), length = 100)     # either
x = range(extrema(d2.weight)..., length = 100)                      # or
for (a, b) in zip(a, b)
    plot!(x, a .+ b .* (x .- x̄), color = "black", alpha = 0.2)
end
plot!(legend = false)

# %% 4.40
b = rand(LogNormal(0.0, 1.0), 10_000)
density(b, xlims = (0, 5))

# %% 4.41
N = 100
a = rand(Normal(178.0, 20.0), N)
b = rand(LogNormal(0.0, 1.0), N)

plot(xlims = extrema(d2.weight), ylims = (-100, 400), xlable = "weight", ylabel = "height")
hline!([0.0, 272.0])
x̄ = mean(d2.weight)
x = range(minimum(d2.weight), maximum(d2.weight), length = 100)
foreach(zip(a, b)) do (a, b)
    plot!(x, a .+ b .* (x .- x̄), color = "black", alpha = 0.2)
end
plot!(legend = false)

# %% 4.42
d = DataFrame(CSV.File(datadir("exp_raw/Howell_1.csv")))
d2 = d[d.age .>= 18, :]
x̄ = mean(d2.weight)

@model function height(weights, heights)
    a ~ Normal(178, 20)
    b ~ LogNormal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ b .* (weights .- x̄)

    # This can't do the predictive posterior:
    # heights .~ Normal.(μ, σ)

    # This is pretty slow but works:
    # for i ∈ eachindex(heights)
    #     heights[i] ~ Normal(μ[i], σ)
    # end

    # This seems to work, fast:
    heights ~ MvNormal(μ, σ)
end

m4_3 = height(d2.weight, d2.height)
q4_3 = quap(m4_3, NelderMead())

# %% 4.43
@model function height_log(weights, heights)
    a ~ Normal(178, 20)
    log_b ~ Normal(0, 1)
    σ ~ Uniform(0, 50)
    μ = a .+ exp(log_b) .* (weights .- mean(weights))
    heights .~ Normal.(μ, σ)
end

m4_3b = height(d2.weight, d2.height)
q4_3b = quap(m4_3b, NelderMead())

# %% 4.44, 4.45
# precis(m4_3)

round.(q4_3.vcov, digits = 3)

# %% 4.46
scatter(d2.weight, d2.height)
post = DataFrame(rand(q4_3.distr, 10_000)', q4_3.params)
a_map = mean(post.a)
b_map = mean(post.b)
plot!(x, a_map .+ b_map .* (x .- x̄))

# %% 4.47
post = DataFrame(rand(q4_3.distr, 10_000)', q4_3.params)
post[1:5, :]

# %% 4.48, 4.49
N = 10  # rerun with 50, 150, 352
dN = d2[1:N, :]
mN = quap(height(dN.weight, dN.height))

post = rand(mN.distr, 20)
post = DataFrame(post', mN.params)

scatter(dN.weight, dN.height)
for p in eachrow(post)
    plot!(x, p.a .+ p.b .* (x .- mean(dN.weight)), color = "black", alpha = 0.3)
end
plot!(legend = false, xlabel = "weight", ylabel = "height")

# %% 4.50 - 4.52
post = rand(q4_3.distr, 1_000)
postdf = DataFrame(post', q4_3.params)
mu_at_50 = postdf.a + postdf.b * (50 - x̄)

density(mu_at_50)

quantile(mu_at_50, (0.1, 0.9))

# %% 4.53 - 4.55
weight_seq = 25:70

# It's a little unfortunate that you have to write out the formula you have already put
# into the model. I don't have a better way at the moment though.
mu = postdf.a' .+ postdf.b' .* (weight_seq .- x̄)

scatter(weight_seq, mu[:, 1:100], legend = false, c = 1, alpha = 0.1)

# %% 4.56, 4.57
mu_m = mean.(eachrow(mu))
mu_lower = quantile.(eachrow(mu), 0.055)
mu_upper = quantile.(eachrow(mu), 0.945)

scatter(d2.weight, d2.height, ms = 3)
plot!(weight_seq, mu_m, ribbon = (mu_m .- mu_lower, mu_upper .- mu_m))

# or

mu = meanlowerupper(mu)

scatter(d2.weight, d2.height, ms = 3)
plot!(weight_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))

# %% 4.59 - 4.61
# This isn't really pretty either. I'm not sure how to put this into a function since
# there isn't a way to know how to create `predict_model` from `height`.
chn = Chains(post', ["a", "b", "σ"])
predict_model = height(weight_seq, missing)
sim = predict(predict_model, chn) |> Array
sim = meanlowerupper(sim')

scatter(d2.weight, d2.height, ms = 3, legend = false)
plot!(weight_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(weight_seq, sim.lower, fillrange = sim.upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %% 4.62
sim = predict(predict_model, Chains(rand(q4_3.distr, 10_000)', ["a", "b", "σ"])) |> Array
sim = meanlowerupper(sim')

scatter(d2.weight, d2.height, ms = 3, legend = false)
plot!(weight_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(weight_seq, sim.lower, fillrange = sim.upper, alpha = 0.3, linealpha = 0.0, c = 2)

# %% 4.63
post = rand(q4_3.distr, 1_000)
postdf = DataFrame(post', q4_3.params)
weight_seq = 25:70
normals = Normal.(postdf.a' .+ postdf.b' .* (weight_seq .- x̄), postdf.σ')
sim = rand.(normals)
sim = meanlowerupper(sim)

scatter(d2.weight, d2.height, ms = 3, legend = false)
plot!(weight_seq, mu.mean, ribbon = (mu.mean .- mu.lower, mu.upper .- mu.mean))
plot!(weight_seq, sim.lower, fillrange = sim.upper, alpha = 0.3, linealpha = 0.0, c = 2)


text:
export   SparsePolynomial

"""
    SparsePolynomial(coeffs::Dict, var)

Polynomials in the standard basis backed by a dictionary holding the
non-zero coefficients. For polynomials of high degree, this might be
advantageous. Addition and multiplication with constant polynomials
are treated as having no symbol.

Examples:

```jldoctest
julia> using Polynomials

julia> P  = SparsePolynomial
SparsePolynomial

julia> p, q = P([1,2,3]), P([4,3,2,1])
(SparsePolynomial(1 + 2*x + 3*x^2), SparsePolynomial(4 + 3*x + 2*x^2 + x^3))

julia> p + q
SparsePolynomial(5 + 5*x + 5*x^2 + x^3)

julia> p * q
SparsePolynomial(4 + 11*x + 20*x^2 + 14*x^3 + 8*x^4 + 3*x^5)

julia> p + 1
SparsePolynomial(2 + 2*x + 3*x^2)

julia> q * 2
SparsePolynomial(8 + 6*x + 4*x^2 + 2*x^3)

julia> p = Polynomials.basis(P, 10^9) - Polynomials.basis(P,0) # also P(Dict(0=>-1, 10^9=>1))
SparsePolynomial(-1.0 + 1.0*x^1000000000)

julia> p(1)
0.0
```

"""
struct SparsePolynomial{T <: Number} <: StandardBasisPolynomial{T}
    coeffs::Dict{Int, T}
    var::Symbol
    function SparsePolynomial{T}(coeffs::Dict{Int, T}, var::SymbolLike) where {T <: Number}
        for (k,v)  in coeffs
            iszero(v) && pop!(coeffs,  k)
        end
        new{T}(coeffs, var)
    end
end

@register SparsePolynomial

function SparsePolynomial{T}(coeffs::AbstractVector{T}, var::SymbolLike=:x) where {T <: Number}
    D = Dict{Int,T}()
    for (i,val) in enumerate(coeffs)
        if !iszero(val)
            D[i-1] = val
        end
    end
    return SparsePolynomial{T}(D, var)
end

function SparsePolynomial{T}(coeffs::OffsetArray{T,1, Array{T, 1}}, var::SymbolLike=:x) where {T <: Number}
    firstindex(coeffs) >= 0 || throw(ArgumentError("Use the `LaurentPolynomial` type for offset arrays with negative first index"))
    D = Dict{Int, T}()
    for i in eachindex(coeffs)
        D[i] = coeffs[i]
    end
    SparsePolynomial{T}(D, var)
end

function SparsePolynomial(coeffs::Dict{Int, T}, var::SymbolLike=:x) where {T <: Number}
    SparsePolynomial{T}(coeffs, Symbol(var))
end

# Interface through `Polynomial`. As with ImmutablePolynomial, this may not be  good idea...
# Deprecated
function Polynomial{T}(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T}
    Base.depwarn("Use of `Polynomial(Dict, var)` is deprecated. Use the `SparsePolynomial` constructor",
            :Polynomial)
    SparsePolynomial{T}(coeffs, var)
end

function Polynomial(coeffs::Dict{Int,T}, var::SymbolLike = :x) where {T}
    Base.depwarn("Use of `Polynomial(Dict, var)` is deprecated. Use the `SparsePolynomial` constructor",
                 :Polynomial)
    SparsePolynomial{T}(coeffs, var)
end

# conversion
function Base.convert(P::Type{<:Polynomial}, q::SparsePolynomial)
    ⟒(P)(coeffs(q), q.var)
end

function Base.convert(P::Type{<:SparsePolynomial}, q::StandardBasisPolynomial{T}) where {T}
    R = promote(eltype(P), T)
    ⟒(P){R}(coeffs(q), q.var)
end

## changes to common
degree(p::SparsePolynomial) = isempty(p.coeffs) ? -1 : maximum(keys(p.coeffs))
function isconstant(p::SparsePolynomial)
    n = length(keys(p.coeffs))
    (n > 1 || (n==1 && iszero(p[0]))) && return false
    return true
end

basis(P::Type{<:SparsePolynomial}, n::Int, var::SymbolLike=:x) =
    SparsePolynomial(Dict(n=>one(eltype(one(P)))), var)

# return coeffs as  a vector
# use p.coeffs to get Dictionary
function  coeffs(p::SparsePolynomial{T})  where {T}

    n = degree(p)
    cs = zeros(T, n+1)
    for (k,v) in p.coeffs
        cs[k+1]=v
    end
    cs

end

# get/set index
function Base.getindex(p::SparsePolynomial{T}, idx::Int) where {T <: Number}
    get(p.coeffs, idx, zero(T))
end

function Base.setindex!(p::SparsePolynomial, value::Number, idx::Int)
    idx < 0  && return p
    if iszero(value)
        haskey(p.coeffs, idx) && pop!(p.coeffs, idx)
    else
        p.coeffs[idx]  = value
    end
    return p
end

Base.firstindex(p::SparsePolynomial) = sort(collect(keys(p.coeffs)), by=x->x[1])[1]
Base.lastindex(p::SparsePolynomial) = sort(collect(keys(p.coeffs)), by=x->x[1])[end]
Base.eachindex(p::SparsePolynomial) = sort(collect(keys(p.coeffs)), by=x->x[1])

# only from tail
function chop!(p::SparsePolynomial{T};
               rtol::Real = Base.rtoldefault(real(T)),
               atol::Real = 0,) where {T}

    for k in sort(collect(keys(p.coeffs)), by=x->x[1], rev=true)
        if isapprox(p[k], zero(T); rtol = rtol, atol = atol)
            pop!(p.coeffs, k)
        else
            return p
        end
    end

    return p

end

function truncate!(p::SparsePolynomial{T};
                   rtol::Real = Base.rtoldefault(real(T)),
                   atol::Real = 0,) where {T}

    max_coeff = maximum(abs, coeffs(p))
    thresh = max_coeff * rtol + atol

    for (k,val) in  p.coeffs
        if abs(val) <= thresh
            pop!(p.coeffs,k)
        end
    end

    return p

end

##
## ----
##

function (p::SparsePolynomial{T})(x::S) where {T,S}

    tot = zero(T) * _one(x)
    for (k,v) in p.coeffs
        tot = _muladd(x^k, v, tot)
    end

    return tot

end

#  map: over values -- not keys
function Base.map(fn, p::P, args...) where {P <: SparsePolynomial}
    ks, vs = keys(p.coeffs), values(p.coeffs)
    vs′ = map(fn, vs, args...)
    _convert(p, Dict(Pair.(ks, vs′)))
end

function Base.:+(p1::SparsePolynomial{T}, p2::SparsePolynomial{S}) where {T, S}

    isconstant(p1) && return p2 + p1[0]
    isconstant(p2) && return p1 + p2[0]

    p1.var != p2.var && error("SparsePolynomials must have same variable")

    R = promote_type(T,S)
    P = SparsePolynomial

    p = zero(P{R}, p1.var)

    # this allocates in the union
#    for i in union(eachindex(p1), eachindex(p2))
#        p[i] = p1[i] + p2[i]
#    end

    # this seems faster
    for i in eachindex(p1)
        p[i] = p1[i] + p2[i]
    end
    for i in eachindex(p2)
        if iszero(p[i])
            @inbounds p[i] = p1[i] + p2[i]
        end
    end

    return  p

end

function Base.:+(p::SparsePolynomial{T}, c::S) where {T, S <: Number}

    R = promote_type(T,S)
    P = SparsePolynomial

    q = zero(P{R}, p.var)
    for k in eachindex(p)
        @inbounds q[k] = R(p[k])
    end
    q[0] = q[0] + c

    return q
end

function Base.:*(p1::SparsePolynomial{T}, p2::SparsePolynomial{S}) where {T,S}

    isconstant(p1) && return p2 * p1[0]
    isconstant(p2) && return p1 * p2[0]
    p1.var != p2.var && error("SparsePolynomials must have same variable")

    R = promote_type(T,S)
    P = SparsePolynomial

    p  = zero(P{R},  p1.var)
    for i in eachindex(p1)
        p1ᵢ = p1[i]
        for j in eachindex(p2)
            @inbounds p[i+j] = muladd(p1ᵢ, p2[j], p[i+j])
        end
    end

    return p

end

function Base.:*(p::P, c::S) where {T, P <: SparsePolynomial{T}, S <: Number}

    R = promote_type(T,S)
    q  = zero(⟒(P){R},  p.var)
    for k in eachindex(p)
        q[k] = p[k] * c
    end

    return q
end

function derivative(p::SparsePolynomial{T}, order::Integer = 1) where {T}

    order < 0 && error("Order of derivative must be non-negative")
    order == 0 && return p

    R = eltype(one(T)*1)
    P = SparsePolynomial
    hasnan(p) && return P(Dict(0 => R(NaN)), p.var)

    n = degree(p)
    order > n && return zero(P{R}, p.var)

    dpn = zero(P{R}, p.var)
    @inbounds for k in eachindex(p)
        dpn[k-order] =  reduce(*, (k - order + 1):k, init = p[k])
    end

    return dpn

end

function integrate(p::SparsePolynomial{T}, k::S) where {T, S<:Number}

    R = eltype((one(T)+one(S))/1)
    P = SparsePolynomial

    if hasnan(p) || isnan(k)
        return P(Dict(0 => NaN), p.var) # not R(NaN)!! don't like XXX
    end

    ∫p = P{R}(R(k), p.var)
    for k in eachindex(p)
        ∫p[k + 1] = p[k] / (k+1)
    end

    return ∫p

end


text:
include("../common/PWGrid_v02.jl")
include("../common/ortho_gram_schmidt.jl")
include("../common/wrappers_fft.jl")

include("EnergiesT.jl")
include("PotentialsT.jl")
include("gen_dr.jl")
include("op_K.jl")
include("op_V_loc.jl")
include("op_H.jl")
include("calc_rho.jl")
include("calc_grad.jl")
include("calc_Energies.jl")
include("KS_solve_Emin_sd.jl")
include("KS_solve_Emin_cg.jl")
include("Poisson_solve.jl")
include("LDA_VWN.jl")
include("Kprec.jl")

include("../common/gen_dr_center.jl")
include("../common/calc_strfact_v2.jl")
include("../common/calc_ewald_v2.jl")

include("diag_lobpcg.jl")
include("KS_solve_SCF.jl")

function test_main( Ns )

    const a = 5.66/0.52917721 # Lattice constant (converted from angstroms to bohrs)
    const LatVecs = a*diagm(ones(3))

    pw = PWGrid( Ns, LatVecs )

    const Ω  = pw.Ω
    const r  = pw.r
    const G  = pw.gvec.G
    const G2 = pw.gvec.G2
    const Npoints = prod(Ns)
    const Ngwx = pw.gvecw.Ngwx

    @printf("Ns   = (%d,%d,%d)\n", Ns[1], Ns[2], Ns[3])
    @printf("Ngwx = %d\n", Ngwx)

    const actual = Npoints/Ngwx
    const theor = 1/(4*pi*0.25^3/3)
    @printf("Compression: actual, theor: %f , %f\n", actual, theor)

    # diamond lattice in cubic cell
    Xpos = zeros(3,8)
    Xpos[:,1] = [0.00 0.00 0.00]
    Xpos[:,2] = [0.25 0.25 0.25]
    Xpos[:,3] = [0.00 0.50 0.50]
    Xpos[:,4] = [0.25 0.75 0.75]
    Xpos[:,5] = [0.50 0.00 0.50]
    Xpos[:,6] = [0.75 0.25 0.75]
    Xpos[:,7] = [0.50 0.50 0.00]
    Xpos[:,8] = [0.75 0.75 0.25]
    Xpos = a*Xpos[:,:]

    for ia = 1:size(Xpos)[2]
        @printf("Ge %18.10f %18.10f %18.10f\n", Xpos[1,ia], Xpos[2,ia], Xpos[3,ia] )
    end

    println(LatVecs)

    Nspecies = 1
    atmsymb = ["Ge"] # unique list of atomic symbols
    atm2species = [1, 1, 1, 1, 1, 1, 1, 1]  # mapping from atom to species
    Zv = [4.0]  # only valence ?

    Sf = calc_strfact( Xpos, Nspecies, atm2species, pw.gvec.G )

    E_nn = calc_ewald( pw, Sf, Xpos, Nspecies, atm2species, Zv )
    @printf("E_nn = %18.10f\n", E_nn)

    Vps = zeros(ComplexF64, Npoints)

    # Ge pseudopotential
    Z = Zv[1]
    λ = 18.5
    rc = 1.052
    Gm = sqrt(G2)

    Vps[:] = -2*pi*exp(-pi*Gm/λ).*cos(Gm*rc).*(Gm/λ)./(1-exp(-2*pi*Gm/λ))
    for n = 0:4
        Vps[:] = Vps[:] + (-1)^n*exp(-λ*rc*n)./(1+(n*λ./Gm).^2)
    end

    Vps[:] = Vps.*4*pi*Z./Gm.^2*(1+exp(-λ*rc)) - 4*pi*Z./Gm.^2

    n = collect(1:4)
    Vps[1] = 4*pi*Z*(1+exp(-λ*rc))*(rc^2/2+1/λ^2 *
             (pi^2/6 + sum((-1).^n.*exp(-λ*rc*n)./n.^2)))

    Vps = Vps[:]/Ω

    V_ionic = real( G_to_R(Ns, Vps .* Sf) ) * Npoints

    # Need to sum up over Nspecies for more than one species type
    # We simply need reshape because we only have one species type here.
    V_ionic = reshape( V_ionic, (Npoints) )

    const Nstates = 16
    Focc = 2.0*ones(Nstates)

    psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates, NiterMax=1000 )
    #psi, Energies, Potentials = KS_solve_Emin_cg( pw, V_ionic, Focc, Nstates,
    #                            NiterMax=1000, Potentials0=Potentials, psi0=psi )

    Y = ortho_gram_schmidt(psi)
    mu = Y' * op_H( pw, Potentials, Y )
    evals, evecs = eig(mu)
    psi = Y*evecs

    #Energies, Potentials, psi, evals = KS_solve_SCF( pw, V_ionic, Focc, Nstates )

    for st = 1:Nstates
        @printf("State # %4d, Energy = %18.10f\n", st, real(evals[st]))
    end

    @printf("E_nn    = %18.10f\n", E_nn)
    @printf("E total = %18.10f\n", E_nn + Energies.Total)

    @printf("\nE total (eV per atom)= %18.10f\n", (E_nn + Energies.Total)/8*27.21 )

end

#@time test_main( [48,48,48] )
test_main( [46, 46, 46] )