text:
"""
    add_reactions(
        m::CoreModelCoupled,
        s::V1,
        b::V2,
        c::AbstractFloat,
        xl::AbstractFloat,
        xu::AbstractFloat;
        check_consistency = false,
    ) where {V1<:VecType,V2<:VecType}

Add reaction(s) to a `CoreModelCoupled` model `m`.

"""
function add_reactions(
    m::CoreModelCoupled,
    s::V1,
    b::V2,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat;
    check_consistency = false,
) where {V1<:VecType,V2<:VecType}
    new_lm = add_reactions(m.lm, s, b, c, xl, xu, check_consistency = check_consistency)
    return CoreModelCoupled(
        new_lm,
        hcat(m.C, spzeros(size(m.C, 1), n_reactions(new_lm) - n_reactions(m.lm))),
        m.cl,
        m.cu,
    )
end

"""
    add_reactions(
        m::CoreModelCoupled,
        s::V1,
        b::V2,
        c::AbstractFloat,
        xl::AbstractFloat,
        xu::AbstractFloat,
        rxn::String,
        mets::K;
        check_consistency = false,
    ) where {V1<:VecType,V2<:VecType,K<:StringVecType}

"""
function add_reactions(
    m::CoreModelCoupled,
    s::V1,
    b::V2,
    c::AbstractFloat,
    xl::AbstractFloat,
    xu::AbstractFloat,
    rxn::String,
    mets::K;
    check_consistency = false,
) where {V1<:VecType,V2<:VecType,K<:StringVecType}
    new_lm = add_reactions(
        m.lm,
        s,
        b,
        c,
        xl,
        xu,
        rxn,
        mets,
        check_consistency = check_consistency,
    )
    return CoreModelCoupled(
        new_lm,
        hcat(m.C, spzeros(size(m.C, 1), n_reactions(new_lm) - n_reactions(m.lm))),
        m.cl,
        m.cu,
    )
end

"""
    add_reactions(
        m::CoreModelCoupled,
        Sp::M,
        b::V,
        c::V,
        xl::V,
        xu::V;
        check_consistency = false,
    ) where {M<:MatType,V<:VecType}

"""
function add_reactions(
    m::CoreModelCoupled,
    Sp::M,
    b::V,
    c::V,
    xl::V,
    xu::V;
    check_consistency = false,
) where {M<:MatType,V<:VecType}
    new_lm = add_reactions(m.lm, Sp, b, c, xl, xu, check_consistency = check_consistency)
    return CoreModelCoupled(
        new_lm,
        hcat(m.C, spzeros(size(m.C, 1), n_reactions(new_lm) - n_reactions(m.lm))),
        m.cl,
        m.cu,
    )
end

"""
    add_reactions(m1::CoreModelCoupled, m2::CoreModel; check_consistency = false)

Add all reactions from `m2` to `m1`.

"""
function add_reactions(m1::CoreModelCoupled, m2::CoreModel; check_consistency = false)
    new_lm = add_reactions(m1.lm, m2, check_consistency = check_consistency)
    return CoreModelCoupled(
        new_lm,
        hcat(m1.C, spzeros(size(m1.C, 1), n_reactions(new_lm) - n_reactions(m1.lm))),
        m1.cl,
        m1.cu,
    )
end

"""
    add_reactions(
        m::CoreModelCoupled,
        Sp::M,
        b::V,
        c::V,
        xl::V,
        xu::V,
        rxns::K,
        mets::K;
        check_consistency = false,
    ) where {M<:MatType,V<:VecType,K<:StringVecType}

"""
function add_reactions(
    m::CoreModelCoupled,
    Sp::M,
    b::V,
    c::V,
    xl::V,
    xu::V,
    rxns::K,
    mets::K;
    check_consistency = false,
) where {M<:MatType,V<:VecType,K<:StringVecType}
    new_lm = add_reactions(
        m.lm,
        Sp,
        b,
        c,
        xl,
        xu,
        rxns,
        mets,
        check_consistency = check_consistency,
    )
    return CoreModelCoupled(
        new_lm,
        hcat(m.C, spzeros(size(m.C, 1), n_reactions(new_lm) - n_reactions(m.lm))),
        m.cl,
        m.cu,
    )
end

"""
    remove_reactions(m::CoreModelCoupled, rxns::Vector{Int})

Remove reaction(s) from a `CoreModelCoupled`.

Also removes any metabolites not involved in any reaction after the deletion.

"""
function remove_reactions(m::CoreModelCoupled, rxns::Vector{Int})
    return CoreModelCoupled(
        remove_reactions(m.lm, rxns),
        m.C[:, filter(e -> e ∉ rxns, 1:n_reactions(m))],
        m.cl,
        m.cu,
    )
end

"""
    remove_reactions(m::CoreModelCoupled, rxn::Integer)

"""
function remove_reactions(m::CoreModelCoupled, rxn::Integer)
    return remove_reactions(m, [rxn])
end

"""
    remove_reactions(m::CoreModelCoupled, rxn::String)

"""
function remove_reactions(m::CoreModelCoupled, rxn::String)
    return remove_reactions(m, [rxn])
end

"""
    remove_reactions(m::CoreModelCoupled, rxns::Vector{String})

"""
function remove_reactions(m::CoreModelCoupled, rxns::Vector{String})
    rxn_indices =
        [findfirst(isequal(name), m.lm.rxns) for name in intersect(rxns, m.lm.rxns)]
    if isempty(rxn_indices)
        return m
    else
        return remove_reactions(m, rxn_indices)
    end
end

"""
Add constraints of the following form to a CoreModelCoupled and return a modified one.

Add constraints to a [`CoreModelCoupled`](@ref) and return a modified one.

The arguments are same as for in-place [`add_coupling_constraints!`](@ref).
"""
function add_coupling_constraints(m::CoreModelCoupled, args...)
    new_lp = deepcopy(m)
    add_coupling_constraints!(new_lp, args...)
    return new_lp
end

"""
    add_coupling_constraints(m::CoreModel, args...)

Add coupling constraints to a plain [`CoreModel`](@ref) (returns a
[`CoreModelCoupled`](@ref)).
"""
add_coupling_constraints(m::CoreModel, args...) =
    add_coupling_constraints(convert(CoreModelCoupled, m), args...)

"""
    add_coupling_constraints!(
        m::CoreModelCoupled,
        c::VecType,
        cl::AbstractFloat,
        cu::AbstractFloat,
    )

Overload for adding a single coupling constraint.
"""
function add_coupling_constraints!(
    m::CoreModelCoupled,
    c::VecType,
    cl::AbstractFloat,
    cu::AbstractFloat,
)
    return add_coupling_constraints!(
        m,
        sparse(reshape(c, (1, length(c)))),
        sparse([cl]),
        sparse([cu]),
    )
end

"""
    add_coupling_constraints!(
        m::CoreModelCoupled,
        C::MatType,
        cl::V,
        cu::V,
    ) where {V<:VecType}

In-place add a single coupling constraint in form
```
    cₗ ≤ C x ≤ cᵤ
```
"""
function add_coupling_constraints!(
    m::CoreModelCoupled,
    C::MatType,
    cl::V,
    cu::V,
) where {V<:VecType}

    all([length(cu), length(cl)] .== size(C, 1)) ||
        throw(DimensionMismatch("mismatched numbers of constraints"))
    size(C, 2) == n_reactions(m) ||
        throw(DimensionMismatch("mismatched number of reactions"))

    m.C = vcat(m.C, sparse(C))
    m.cl = vcat(m.cl, sparse(cl))
    m.cu = vcat(m.cu, sparse(cu))
end

"""
    remove_coupling_constraints(m::CoreModelCoupled, args...)

Remove coupling constraints from the linear model, and return the modified
model. Arguments are the same as for in-place version
[`remove_coupling_constraints!`](@ref).
"""
function remove_coupling_constraints(m::CoreModelCoupled, args...)
    new_model = deepcopy(m)
    remove_coupling_constraints!(new_model, args...)
    return new_model
end

"""
    remove_coupling_constraints!(m::CoreModelCoupled, constraint::Int)

Removes a single coupling constraints from a [`CoreModelCoupled`](@ref)
in-place.
"""
remove_coupling_constraints!(m::CoreModelCoupled, constraint::Int) =
    remove_coupling_constraints!(m, [constraint])

"""
    remove_coupling_constraints!(m::CoreModelCoupled, constraints::Vector{Int})

Removes a set of coupling constraints from a [`CoreModelCoupled`](@ref)
in-place.
"""
function remove_coupling_constraints!(m::CoreModelCoupled, constraints::Vector{Int})
    to_be_kept = filter(e -> e ∉ constraints, 1:n_coupling_constraints(m))
    m.C = m.C[to_be_kept, :]
    m.cl = m.cl[to_be_kept]
    m.cu = m.cu[to_be_kept]
end

"""
    change_coupling_bounds!(
        model::CoreModelCoupled,
        constraints::Vector{Int};
        cl::V = Float64[],
        cu::V = Float64[],
    ) where {V<:VecType}

Change the lower and/or upper bounds (`cl` and `cu`) for the given list of
coupling constraints.
"""
function change_coupling_bounds!(
    model::CoreModelCoupled,
    constraints::Vector{Int};
    cl::V = Float64[],
    cu::V = Float64[],
) where {V<:VecType}
    found = [index ∈ 1:n_coupling_constraints(model) for index in constraints]
    red_constraints = constraints[found]

    length(red_constraints) == length(unique(red_constraints)) ||
        error("`constraints` appears to contain duplicates")
    if !isempty(cl)
        length(constraints) == length(cl) ||
            throw(DimensionMismatch("`constraints` size doesn't match with `cl`"))
        model.cl[red_constraints] = cl[found]
    end

    if !isempty(cu)
        length(constraints) == length(cu) ||
            throw(DimensionMismatch("`constraints` size doesn't match with `cu`"))
        model.cu[red_constraints] = cu[found]
    end
end

"""
    change_bounds!(
        model::CoreModelCoupled,
        rxns::Vector{Int};
        xl::V = Float64[],
        xu::V = Float64[],
    )

Change the lower and/or upper bounds ('xl' and 'xu') for given reactions.

"""
function change_bounds!(
    model::CoreModelCoupled,
    rxns::Vector{Int};
    xl::V = Float64[],
    xu::V = Float64[],
) where {V<:VecType}
    change_bounds!(model.lm, rxns, xl = xl, xu = xu)
end

"""
    change_bounds!(
        model::CoreModelCoupled,
        rxns::Vector{String};
        xl::V = Float64[],
        xu::V = Float64[],
    ) where {V<:VecType}

"""
function change_bounds!(
    model::CoreModelCoupled,
    rxns::Vector{String};
    xl::V = Float64[],
    xu::V = Float64[],
) where {V<:VecType}
    change_bounds!(model.lm, rxns, xl = xl, xu = xu)
end


text:
using StructuralDDP

# load models
include("models/CapitalAdjustModel.jl")
include("models/CapitalAdjustModel2.jl")
createmodel(model::Symbol; kwargs...) = eval(model)(; kwargs...)

function warn_time(t, border, preprint)
	t < border || @warn preprint * " took longer than usual: $t vs $border seconds"
end
function warn_memory(bytes, border, preprint)
	bytes < border || @warn preprint * " allocated more memory than usual: " *
		Base.format_bytes(bytes) * " vs. " * Base.format_bytes(border)
end

# Single choice variable
model = :CapitalAdjustModel
dipar = Dict(:nK => 40, :nz => 5, :γ => 0.5, :β => 0.9)
p_1_All = createmodel(model; dipar..., intdim = :All)
p_1_Separable = createmodel(model; dipar..., intdim = :Separable)
p_1_Separable_States = createmodel(model; dipar..., intdim = :Separable_States)
p_1_Separable_ExogStates = createmodel(model; dipar..., intdim = :Separable_ExogStates)

diopt = Dict(:concavity => false, :monotonicity => false, :rewardcall=>:pre)

sol_1_All = solve(p_1_All; diopt...)
val, t, bytes = @timed solve(p_1_All; diopt...)
warn_time(t, 0.12, "CapitalAdjustModel All: Solver")
warn_memory(bytes, 75555555, "CapitalAdjustModel All: Solver")

sol_1_Separable = solve(p_1_Separable; diopt...)
val, t, bytes = @timed solve(p_1_Separable; diopt...)
warn_time(t, 0.045, "CapitalAdjustModel Separable: Solver")
warn_memory(bytes, 15000000, "CapitalAdjustModel Separable: Solver")

sol_1_Separable_States = solve(p_1_Separable_States; diopt...)
val, t, bytes = @timed solve(p_1_Separable_States; diopt...)
warn_time(t, 0.007, "CapitalAdjustModel Separable_States: Solver")
warn_memory(bytes, 580000, "CapitalAdjustModel Separable_States: Solver")

sol_1_Separable_ExogStates = solve(p_1_Separable_ExogStates; diopt...)
val, t, bytes = @timed solve(p_1_Separable_ExogStates; diopt...)
warn_time(t, 0.006, "CapitalAdjustModel Separable_ExogStates: Solver")
warn_memory(bytes, 230000, "CapitalAdjustModel Separable_ExogStates: Solver")

# SIMULATE
shocks = drawshocks(p_1_All, nFirms=10^4, nPeriods=50)

simulate(sol_1_All, shocks)
val, t, bytes = @timed simulate(sol_1_All, shocks)
warn_time(t, 0.25, "CapitalAdjustModel All: Simulator")
warn_memory(bytes, 231000000, "CapitalAdjustModel All: Simulator")

simulate(sol_1_Separable, shocks)
val, t, bytes = @timed simulate(sol_1_Separable, shocks)
warn_time(t, 0.25, "CapitalAdjustModel Separable: Simulator")
warn_memory(bytes, 260000000, "CapitalAdjustModel Separable: Simulator")

simulate(sol_1_Separable_States, shocks)
val, t, bytes = @timed simulate(sol_1_Separable_States, shocks)
warn_time(t, 0.24, "CapitalAdjustModel Separable_States: Simulator")
warn_memory(bytes, 260000000, "CapitalAdjustModel Separable_States: Simulator")

simulate(sol_1_Separable_ExogStates, shocks)
val, t, bytes = @timed simulate(sol_1_Separable_ExogStates, shocks)
warn_time(t, 0.23, "CapitalAdjustModel Separable_ExogStates: Simulator")
warn_memory(bytes, 191000000, "CapitalAdjustModel Separable_ExogStates: Simulator")

# Two choice variables
model = :CapitalAdjustModel2
dipar = Dict(:nK => 10, :nN => 10, :nz => 5, :β => 0.9)
p_2_All = createmodel(model; dipar..., intdim = :All)
p_2_Sep = createmodel(model; dipar..., intdim = :Separable)
p_2_Sep_States = createmodel(model; dipar..., intdim = :Separable_States)
p_2_Sep_ExogStates = createmodel(model; dipar..., intdim = :Separable_ExogStates)

diopt = Dict(:monotonicity=>[false,false], :concavity=>[false,false],
	:rewardcall=>:pre)

sol_2_All = solve(p_2_All; diopt...)
val, t, bytes = @timed solve(p_2_All; diopt...)
warn_time(t, 2.1, "CapitalAdjustModel2 All: Solver")
warn_memory(bytes, 870000000, "CapitalAdjustModel2 All: Solver")

sol_2_Sep = solve(p_2_Sep; diopt...)
val, t, bytes = @timed solve(p_2_Sep; diopt...)
warn_time(t, 0.45, "CapitalAdjustModel2 Separable: Solver")
warn_memory(bytes, 120000000, "CapitalAdjustModel2 Separable: Solver")

sol_2_Sep_States = solve(p_2_Sep_States; diopt...)
val, t, bytes = @timed solve(p_2_Sep_States; diopt...)
warn_time(t, 0.031, "CapitalAdjustModel2 Separable_States: Solver")
warn_memory(bytes, 1700000, "CapitalAdjustModel2 Separable_States: Solver")

sol_2_Sep_ExogStates = solve(p_2_Sep_ExogStates; diopt...)
val, t, bytes = @timed solve(p_2_Sep_ExogStates; diopt...)
warn_time(t, 0.03, "CapitalAdjustModel2 Separable_ExogStates: Solver")
warn_memory(bytes, 570000, "CapitalAdjustModel2 Separable_ExogStates: Solver")

# SIMULATE
shocks = drawshocks(p_2_All, nFirms=10^3, nPeriods=50)

simulate(sol_2_All, shocks)
val, t, bytes = @timed simulate(sol_2_All, shocks)
warn_time(t, 0.12, "CapitalAdjustModel2 All: Simulator")
warn_memory(bytes, 62100000, "CapitalAdjustModel2 All: Simulator")

simulate(sol_2_Sep, shocks)
val, t, bytes = @timed simulate(sol_2_Sep, shocks)
warn_time(t, 0.07, "CapitalAdjustModel2 Separable: Simulator")
warn_memory(bytes, 48000000, "CapitalAdjustModel2 Separable: Simulator")

simulate(sol_2_Sep_States, shocks)
val, t, bytes = @timed simulate(sol_2_Sep_States, shocks)
warn_time(t, 0.07, "CapitalAdjustModel2 Separable_States: Simulator")
warn_memory(bytes, 48000000, "CapitalAdjustModel2 Separable_States: Simulator")

simulate(sol_2_Sep_ExogStates, shocks)
val, t, bytes = @timed simulate(sol_2_Sep_ExogStates, shocks)
warn_time(t, 0.07, "CapitalAdjustModel2 Separable_ExogStates: Simulator")
warn_memory(bytes, 48000000, "CapitalAdjustModel2 Separable_ExogStates: Simulator")

ext: jl


text:
using Thebes, Luxor

include(dirname(pathof(Thebes)) * "/../data/moreobjects.jl")

function main()
    eyepoint(0, 100, 300)
    perspective(300)
    Drawing(500, 500, "docs/src/assets/logo.svg")
    origin()
    squircle(O, 240, 240, :clip, rt=0.3)
    background("midnightblue")

    sethue(0.0, 0.3, 0.5)
    circle(O + (0, 850), 800, :fill)

    object = make(rhombitruncated_icosidodecahedron, "solid")
    scaleby!(object, 180, 180, 180)
    rotateby!(object, π/3, 0, 0)
    sortfaces!(object, eyepoint=eyepoint())
    pin(object, gfunction = (args...) -> begin
        vertices, faces, labels = args
        cols = (Luxor.julia_green, Luxor.julia_red,Luxor.julia_purple, Luxor.julia_blue)
        setlinejoin("bevel")
        if !isempty(faces)
            @layer begin
                for (n, p) in enumerate(faces)
                    setopacity(0.8)
                    sethue(cols[mod1(n, end)])
                    poly(p, :fillpreserve, close=true)
                    setline(6)
                    sethue("gold")
                    strokepreserve()
                    setline(2)
                    sethue("black")
                    strokepreserve()
                    setline(1)
                    sethue("white")
                    strokepath()
                end
            end
        end
    end)
    finish()
    preview()
end

main()


text:
"""
    Code for a half-cosine mask shape for use with CCF
Author: Eric Ford
Created: August 2020
"""

"""   CosCCFMask
Cosine mask with one parameter, it's quarter period, i.e., where to truncate it, as a velocity in m/s.
Mask weights are stored separately in a line list.
"""
struct CosCCFMask <: AbstractCCFMaskShape
    half_width::Float64

    """ CosCCFMask( full_width ) """
    function CosCCFMask( w::Real )
        @assert 100 <= w <= 300000   # 0.1 and 300km/s are arbitrary choices for an upper limit
        new(w/2)
    end

end

""" `CosCCFMask( inst  )` """
function CosCCFMask(inst::AbstractInstrument; scale_factor::Real = 1)
    w = scale_factor * default_ccf_mask_v_width(inst)
    CosCCFMask(w)
end

λ_min(m::CosCCFMask,λ::Real) = λ/calc_doppler_factor(m.half_width)
λ_max(m::CosCCFMask,λ::Real) = λ*calc_doppler_factor(m.half_width)

function integrate(m::CosCCFMask, v_lo::Real,v_hi::Real)
    #quadgk(m, v_lo, v_hi)[1]
    0.5*(sin(0.5*π*v_hi/m.half_width)-sin(0.5*π*v_lo/m.half_width))
end

""" Functor for returning PSF for Δv <= half_width.  """
function (m::CosCCFMask)(Δv::Real)
    if abs2(Δv) > abs2(m.half_width)
        return zero(Δv)
    else
        return cos(0.5*π*Δv/m.half_width)*π/(4*m.half_width)
    end
end

function mask_with_increased_fwhm(m::CosCCFMask, Δfwhm::Real )
    width = 2*m.half_width+Δfwhm
    return CosCCFMask(width)
end

max_stars_repo_name: RvSpectML/EchelleCCFs

max_issues_repo_name: RvSpectML/EchelleCCFs

max_forks_repo_name: RvSpectML/EchelleCCFs


text:
OVERRIDE_MSG = "This should be overridden in concrete implemention of AbstractState"

last_reward(::AbstractState)            = @assert false OVERRIDE_MSG
# @todo change to get_next_state
next_state(::AbstractState, ::Integer)  = @assert false OVERRIDE_MSG
states(::AbstractState)                 = @assert false OVERRIDE_MSG
dynamics(::AbstractState, ::Integer)    = @assert false OVERRIDE_MSG

# @todo make this into a generic string
hex_str(::AbstractState)                = @assert false OVERRIDE_MSG
key(::AbstractState)                    = @assert false OVERRIDE_MSG
to_play(::AbstractState)                = @assert false OVERRIDE_MSG
make_image(::AbstractState)             = @assert false OVERRIDE_MSG
gridsvg(::AbstractState)                = @assert false OVERRIDE_MSG
targetsvg(::AbstractState, ::Target)    = @assert false OVERRIDE_MSG

# @todo should be in AbstractBoardGame
function cartesian_coords_to_action(
    ::Type{<:AbstractState},
    x::Int,
    y::Int)
    @assert false OVERRIDE_MSG
end

# @todo should be in AbstractBoardGame
# @todo rename to get_cell_states
calc_cell_states(::AbstractState)       = @assert false OVERRIDE_MSG

# @todo using the decode{T} pattern
encode(s::AbstractEnvState) = base64encode(serialize, s)
function decode(encoded::String)
    encoded |> base64decode |> IOBuffer |> deserialize
end

function init(T::Type{<:AbstractState}, actions::Vector{<:Integer})
    reduce(next_state, actions; init=T())
end

Base.copy(T::Type{<:AbstractState}, state::AbstractState) = T(key(state))

function Base.show(io::IO, s::AbstractState)
    encoded = encode(s)
    if length(encoded) >= 40
        encoded = join([first(encoded, 20), "...", last(encoded, 20)])
    end
    print(io, "$(string(AbstractState)){\"$encoded\"}")
end


text:
function hdrswap!(BUF::SeisIOBuf)
  u16 = getfield(BUF, :u16)
  @inbounds for i = 1:5
    u16[i] = bswap(u16[i])
  end
  setfield!(BUF, :n, bswap(getfield(BUF, :n)))
  setfield!(BUF, :r1, bswap(getfield(BUF, :r1)))
  setfield!(BUF, :r2, bswap(getfield(BUF, :r2)))
  setfield!(BUF, :tc, bswap(getfield(BUF, :tc)))
  return nothing
end

function update_dt!(BUF::SeisIOBuf)
  r1 = getfield(BUF, :r1)
  r2 = getfield(BUF, :r2)
  dt = 0.0
  if r1 > 0 && r2 > 0
    dt = 1.0/Float64(r1*r2)
  elseif r1 > 0
    dt = -1.0*Float64(r2/r1)
  elseif r2 > 0
    dt = -1.0*Float64(r1/r2)
  else
    dt = Float64(r1*r2)
  end
  setfield!(BUF, :dt, dt)
  setfield!(BUF, :Δ, round(Int64, sμ*dt))
  setfield!(BUF, :r1_old, r1)
  setfield!(BUF, :r2_old, r2)
  return nothing
end

function update_hdr!(BUF::SeisIOBuf)
  id_j = 0
  id = getfield(BUF, :id)
  hdr = getfield(BUF, :hdr)
  for p in id_positions
    if hdr[p] != 0x20
      id_j += 1
      id[id_j] = hdr[p]
    end
    if p == 12 || p == 5 || p == 7
      id_j += 1
      id[id_j] = id_spacer
    end
  end
  unsafe_copyto!(getfield(BUF, :hdr_old), 1, getfield(BUF, :hdr), 1, 12)
  setfield!(BUF, :id_str, unsafe_string(pointer(getfield(BUF, :id)), id_j))
  return nothing
end

###############################################################################
function parserec!(S::SeisData, BUF::SeisIOBuf, sid::IO, v::Int64, nx_new::Int64, nx_add::Int64)
  # =========================================================================
  u16 = getfield(BUF, :u16)
  flags = getfield(BUF, :flags)
  xi = 0
  te = 0

  # Fixed section of data header (48 bytes)
  pos = position(sid)
  read!(sid, BUF.seq)
  read!(sid, BUF.hdr)
  if v > 2
      println(stdout, join(map(Char,BUF.seq), map(Char,BUF.hdr)))
  end
  u16[1]          = read(sid, UInt16)
  u16[2]          = read(sid, UInt16)
  hh              = read(sid, UInt8)
  mm              = read(sid, UInt8)
  ss              = read(sid, UInt8)
  skip(sid, 1)
  u16[3]          = read(sid, UInt16)
  BUF.n          = read(sid, UInt16)
  BUF.r1         = read(sid, Int16)
  BUF.r2         = read(sid, Int16)
  read!(sid, flags)
  BUF.tc         = read(sid, Int32)
  u16[4]          = read(sid, UInt16)
  u16[5]          = read(sid, UInt16)

  if getfield(BUF, :swap) == true
    hdrswap!(BUF)
  end

  # =========================================================================
  # Post-read header processing

  # This is the standard check for correct byte order...?
  yy = u16[1]
  jj = u16[2]
  if (jj > 0x0200 || ((jj == 0x0000 || jj == 0x0100) &&
      (yy > 0x0907 || yy < 0x707)) || yy>0x0bb8)
	  setfield!(BUF, :swap, !BUF.swap)
    if ((BUF.swap == true) && (BUF.wo == 0x01))
      BUF.xs = true
    end
    hdrswap!(BUF)
  end

  if BUF.r1 != BUF.r1_old || BUF.r2 != BUF.r2_old
    update_dt!(BUF)
  end

  n = getfield(BUF, :n)

  # =========================================================================
  # Channel handling for S

  # Check this SEED id and whether or not it exists in S
  if BUF.hdr != BUF.hdr_old
    update_hdr!(BUF)
  end
  id = getfield(BUF, :id_str)
  c = findid(id, S)

  if c == 0
    if v > 2
      println(stdout, "New channel; ID = ", id, ", S.id = ", S.id)
    end
    x = Array{Float32, 1}(undef, nx_new)
    L = nx_new
    nt = 2

    C = SeisChannel()
    setfield!(C, :id, id)
    setfield!(C, :name, identity(id))
    setfield!(C, :fs, 1.0/getfield(BUF, :dt))
    setfield!(C, :x, x)
    push!(S, C)
    c = S.n

    (v > 1) && println(stdout, "Added channel: ", id)
  else
    # assumes fs doesn't change within a SeisData structure
    t = getindex(getfield(S, :t), c)
    x = getindex(getfield(S, :x), c)
    nt = div(lastindex(t), 2)
    L = lastindex(x)

    if nt > 0
      xi = getindex(t, nt)
      te = endtime(t, getindex(getfield(S, :fs), c))
    end
    if xi + n > L
      resize!(x, xi + max(n, nx_add))
      (v > 1) && println(stdout, id, ": resized from length ", L,
                                      " to length ", length(x))
    end
  end

  # =========================================================================
  # Parse blockettes

  nsk = u16[4] - 0x0030
  u16[6] = u16[5] - 0x0030
  nblk = flags[4]
  v > 2 && println(string("Blockettes to read: ", nblk))
  @inbounds for i = 0x01:0x01:nblk

    # DND DND DND
    skip(sid, u16[6])
    nsk = nsk - u16[6]
    u16[5] = UInt16(position(sid) - pos)
    # DND DND DND

    bt            = read(sid, UInt16)
    u16[6]        = read(sid, UInt16)
    if getfield(BUF, :swap) == true
      bt = bswap(bt)
      setindex!(u16, bswap(u16[6]), 6)
    end

    # debug
    if v > 2
      printstyled(string("Position = ", position(sid), "\n"), color=:light_green)
      printstyled(string("Blockette type to read: ", bt, "\n"), color=:light_yellow)
      println(stdout, "Relative position u16[5] = ", u16[5], " bytes from record begin")
      println(stdout, "We are nsk = ", nsk, " bytes to data begin")
    end

    # Blockette parsing moved to individual functions
    if bt == 0x03e8
      blk_len = blk_1000(S, sid, c)
    elseif bt == 0x03e9
      blk_len = blk_1001(S, sid, c)
    elseif bt == 0x0064
      blk_len = blk_100(S, sid, c)
    elseif bt == 0x00c9
      blk_len = blk_201(S, sid, c)
    elseif bt == 0x01f4
      blk_len = blk_500(S, sid, c)
    elseif bt == 0x07d0
      blk_len = blk_2000(S, sid, c)
    elseif bt in BUF.calibs
      blk_len = blk_calib(S, sid, c, bt)
    else
      v > 1 && println(stdout, id, ": no support for Blockette Type ", bt, "; skipped.")
      blk_len = (u16[6] == 0x0000 ? nsk : u16[6])
      skip(sid, blk_len - 0x0004)
    end
    nsk = nsk - blk_len
    if u16[6] != 0x0000
      u16[6] = u16[6] - blk_len - u16[5]
    end
  end

  # =========================================================================
  # Data parsing: originally adapted from rdmseed.m by Francois Beauducel
  # (not very similar anymore)
  if nsk > 0x0000
    skip(sid, Int(nsk))
  end

  # Get data format
  fmt = getfield(BUF, :fmt)
  nb = getfield(BUF, :nx) - u16[4]

  # debug output
  if v > 2
    printstyled(string("Position = ", position(sid), "\n"), color=:light_green)
    println(stdout, "To parse: nx = ", n, " sample blockette, ",
    "compressed size = ", nb, " bytes, fmt = ", fmt)
  end

  if fmt == 0x0a || fmt == 0x0b
    SEED_Steim!(sid, BUF, nb)
  elseif fmt == 0x00
    # ASCII is a special case as it's typically not data
    D = getindex(getfield(S, :misc), c)
    if !haskey(D, "seed_ascii")
      D["seed_ascii"] = Array{String,1}(undef,0)
    end
    push!(D["seed_ascii"], SEED_Char(sid, BUF, nb))
  elseif fmt in UInt8[0x01, 0x03, 0x04]
    SEED_Unenc!(sid, S, c, xi, nb, n)
  elseif fmt == 0x05
    SEED_Float64!(sid, S, c, xi, nb)
  elseif fmt == 0x0d || fmt == 0x0e
    SEED_Geoscope!(sid, BUF)
  elseif fmt == 0x10
    SEED_CDSN!(sid, BUF)
  elseif fmt == 0x1e
    SEED_SRO!(sid, BUF)
  elseif fmt == 0x20
    SEED_DWWSSN!(sid, BUF)
  else
    warn_str = string("readmseed, unsupported format = ", fmt, ", ", nb, " bytes skipped.")
    @warn(warn_str); note!(S, c, warn_str)
    skip(sid, nb)
    return nothing
  end

  if fmt > 0x00
    # Update S.x[c]
    if fmt > 0x05
      unsafe_copyto!(x, xi+1, getfield(BUF, :x), 1, getfield(BUF, :k))
    end

    # Update S.t[c]

    # Check for time correction
    is_tc = flags[2] >> 1 & 0x01
    tc = getfield(BUF, :tc)
    if is_tc == false && tc != zero(Int32)
      δt = Int64(tc)*100
    else
      δt = zero(Int64)
    end

    # Sample rate in μs
    Δ = getfield(BUF, :Δ)

    # Elapsed time since S.t[c] ended
    τ = y2μs(u16[1]) +
        Int64(u16[2]-one(UInt16))*86400000000 +
        Int64(hh)*3600000000 +
        Int64(mm)*60000000 +
        Int64(ss)*1000000 +
        Int64(u16[3])*100 +
        δt# -
        # te -
        # Δ

    # New channel
    if te == 0
      t = Array{Int64, 2}(undef, 2, 2)
      setindex!(t, one(Int64), 1)
      setindex!(t, n, 2)
      setindex!(t, τ, 3)
      setindex!(t, zero(Int64), 4)
      setindex!(getfield(S, :t), t, c)

    # Existing channel
    else
      check_for_gap!(S, c, τ, n, v)
    end

    v > 2 && printstyled(string("Position = ", position(sid), "\n"), color=:light_green)
  end
  return nothing
end


text:
using DataDeps

register(DataDep(
    "spiral",
    """
    Dataset: spiral
    Credit: https://smorbieu.gitlab.io/generate-datasets-to-understand-some-clustering-algorithms-behavior/

    A simulated dataset of three noisy spirals.

    Observations: 1000
    Features:     2
    Classes:      3
    """,
    "script", # nothing to download
    "51d3cf8d223edcacb65db752ff3b0487722951b6e00cb566c3e41a574189c20e";
    fetch_method=function(unused, localdir)
        # Parameters
        N           = 1000
        max_A       = 600
        max_B       = 300
        max_C       = N - max_A - max_B
        max_radius  = 7.0
        x0          = -3.5
        y0          = 3.5
        angle_start = π / 8
        seed        = 1903

        # Simulate the data.
        rng = StableRNG(seed)
        target, x, y = Vector{Char}(undef, N), zeros(N), zeros(N)
        for i in 1:N
        if i ≤ max_A
            (class, k, n, θ) = ('A', i, max_A, angle_start)
            noise = 0.1
        elseif i ≤ max_A + max_B
            (class, k, n, θ) = ('B', i-max_A+1, max_B, angle_start + 2π/3)
            noise = 0.2
        else
            (class, k, n, θ) = ('C', i-max_A-max_B+1, max_C, angle_start + 4π/3)
            noise = 0.3
        end
        # Compute coordinates.
        angle = θ + π * k / n
        radius = max_radius * (1 - k / (n + n / 5))

        target[i] = class
        x[i] = x0 + radius*cos(angle) + noise*randn(rng)
        y[i] = y0 + radius*sin(angle) + noise*randn(rng)
        end
        local_file = joinpath(localdir, "data.csv")
        df = DataFrame(target=target, x1=x, x2=y)
        perm = Random.randperm(rng, size(df, 1))
        foreach(col -> permute!(col, perm), eachcol(df))
        CSV.write(local_file, df)
        return local_file
    end
))


text:
using Test
using Science_ABM
using Agents

include("model.jl")
include("basicfunctions.jl")


text:
@testset "706.design-hashmap.jl" begin
    hm = MyHashMap()
    @test put_706!(hm, 1, 1) == nothing
    @test put_706!(hm, 2, 2) == nothing
    @test get_706(hm, 1) == 1
    @test get_706(hm, 3) == -1
    @test put_706!(hm, 2, 1) == nothing
    @test get_706(hm, 2) == 1
    @test remove_706!(hm, 2) == nothing
    @test get_706(hm, 2) == -1
end


text:
#testfile.jl
#= Jacob Child
    Jan 26, 2022
        I learned how to comment and commit!!!
        I will be using test code from Medium
=#

function sphere_vol(r) #near as i can tell this creates a function with input r
    return 4/3*pi*r^3
end

print("volume: ", sphere_vol(5)) #looks like you can print and call all in one line. Go Julia!
#it ran! It also outputs 13 decimal places! answer= 523.599 fyi


text:
  using Base.Test
  using Sigma
  import Sigma: LazyOmega, flip
  import AbstractDomains: tf, f, t
  # When cond is a random variable, @If and ifelse return a random variable
  begin
    #FIXME adding local before x and a caused errors
    x = flip(0.6)
    a = ifelse(x,false,true)
    @show typeof(a)
    @test isa(a,Sigma.RandVar)
    @test call(a,LazyOmega()) === tf
    @test call(a,[0.3]) == false
    @test call(a,[0.7]) == true

    a = ifelse(x,false,true)
    @test isa(a,Sigma.RandVar)
    @test call(a,LazyOmega()) === tf
    @test call(a,[0.3]) == false
    @test call(a,[0.7]) == true
  end

# When cond is random variable, and a/b also random variable
# @If and ifelse should return random variable which "pipes"
# input to consequent and/or alternative.
begin
  x = flip(0.6)
  a = ifelse(x,flip(0.9),flip(0.2))
  @test isa(a,Sigma.RandVar)
  # LazyOmega() should be piped through to a or b, and hence:
  @test !isa(call(a,LazyOmega()),Sigma.RandVar)
  @test call(a,[0.3]) == true
  @test call(a,[0.7]) == false

  a = ifelse(x,flip(0.9), flip(0.2))
  @test isa(a,Sigma.RandVar)
  # LazyOmega() should be piped through to a or b, and hence:
  @test !isa(call(a,LazyOmega()),Sigma.RandVar)
  @test call(a,[0.3]) == true
  @test call(a,[0.7]) == false
end


text:
"""
$(TYPEDSIGNATURES)

Two panel visualization of gridfactory with input and resulting grid
See [`default_options`](@ref) for available `kwargs`.
"""
builderplot(gb::SimplexGridBuilder; Plotter=nothing,kwargs...)= builderplot(plottertype(Plotter), gb,Plotter;kwargs...)

builderplot(::Type{Nothing}, builder::SimplexGridBuilder, PyPlot ;kwargs...)=nothing

function builderplot(::Type{PyPlotType}, builder::SimplexGridBuilder, PyPlot ;kwargs...)

    p=GridVisualizer(Plotter=PyPlot, layout=(1,2); kwargs...)

    opts=blendoptions!(copy(builder.options);kwargs...)

    Triangulate=builder.Generator
    @assert(istriangulate(Triangulate))

    flags=makeflags(opts,:triangle)

    if opts[:verbose]
        @show flags
    end

    triin=nothing
    try
        triin=triangulateio(builder)
    catch err
        @error "Incomplete geometry description"
        rethrow(err)
    end

    if !isnothing(opts[:unsuitable])
        Triangulate.triunsuitable(opts[:unsuitable])
    end

    triout,vorout=Triangulate.triangulate(flags,triin)
    PyPlot.subplot(121)
    PyPlot.title("In")
    Triangulate.plot(PyPlot,triin)
    PyPlot.subplot(122)
    PyPlot.title("Out")
    Triangulate.plot(PyPlot,triout)
    PyPlot.tight_layout()
    PyPlot.gcf()
end


text:
using Test

@testset "Test Dijkstra.jl" begin
    @testset "test_graph" begin
        include("test_graph.jl")
    end
end


text:
function ALcheck(userinput; verbose=true)
    try
      inputpass = string(replace(userinput, r"^(check|prove|all|any|✓)[\\:\\-\\ ]*"=>""))
      occursin(r"^(check|✓)[\\:\\-\\ ]+", userinput) &&
        checkfeasible(inputpass, replset, verbose = verbose)
      occursin(r"^(prove|all)[\\:\\-\\ ]+", userinput) &&
        checkfeasible(inputpass, replset, countall=true, verbose = verbose)
      occursin(r"^any[\\:\\-\\ ]+", userinput) &&
        checkfeasible(inputpass, replset, countany=true, verbose = verbose)
      # push!(logicset, replset)
  catch er
      replthrow("\nWarning! Check Fail: $er")
      (length(userinput) == 5) && println("Nothing to check")
    end
end


text:
# ---
# title: 804. Unique Morse Code Words
# id: problem804
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/unique-morse-code-words/description/>
# hidden: true
# ---
#
# International Morse Code defines a standard encoding where each letter is
# mapped to a series of dots and dashes, as follows: `"a"` maps to `".-"`, `"b"`
# maps to `"-..."`, `"c"` maps to `"-.-."`, and so on.
#
# For convenience, the full table for the 26 letters of the English alphabet is
# given below:
#
#
#
#     [".-","-...","-.-.","-..",".","..-.","--.","....","..",".---","-.-",".-..","--","-.","---",".--.","--.-",".-.","...","-","..-","...-",".--","-..-","-.--","--.."]
#
# Now, given a list of words, each word can be written as a concatenation of the
# Morse code of each letter. For example, "cab" can be written as "-.-..--...",
# (which is the concatenation "-.-." \+ ".-" \+ "`-...`"). We'll call such a
# concatenation, the transformation of a word.
#
# Return the number of different transformations among all words we have.
#
#
#
#     **Example:**
#     Input: words = ["gin", "zen", "gig", "msg"]
#     Output: 2
#     Explanation:
#     The transformation of each word is:
#     "gin" -> "--...-."
#     "zen" -> "--...-."
#     "gig" -> "--...--."
#     "msg" -> "--...--."
#
#     There are 2 different transformations, "--...-." and "--...--.".
#
#
# **Note:**
#
#   * The length of `words` will be at most `100`.
#   * Each `words[i]` will have length in range `[1, 12]`.
#   * `words[i]` will only consist of lowercase letters.
#
#
## @lc code=start
using LeetCode

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


text:
"""
    stepforward!(prob)

Step forward `prob` one time step.
"""
function stepforward!(prob::Problem)
  stepforward!(prob.sol, prob.clock, prob.timestepper, prob.eqn, prob.vars, prob.params, prob.grid)
end

"""
    stepforward!(prob, nsteps)

Step forward `prob` for `nsteps`.
"""
function stepforward!(prob::Problem, nsteps::Int)
  for step = 1:nsteps
    stepforward!(prob)
  end
  nothing
end

"""
    stepforward!(prob, diags, nsteps)

Step forward `prob` for `nsteps`, incrementing `diags` along the way. `diags` may be a single `Diagnostic` or
a `Vector` of `Diagnostic`s.
"""
function stepforward!(prob::Problem, diags, nsteps::Int)
  for step = 1:nsteps
    stepforward!(prob)
    increment!(diags)
  end
  nothing
end

const fullyexplicitsteppers= [
  :ForwardEuler,
  :RK4,
  :AB3,
  :FilteredForwardEuler,
  :FilteredRK4,
  :FilteredAB3
]

isexplicit(stepper) = any(Symbol(stepper) .== fullyexplicitsteppers)

"""
    TimeStepper(stepper, eq, dt=nothing)

Generalized timestepper constructor. If `stepper` is explicit, `dt` is not used.
"""
function TimeStepper(stepper, eq, dt=nothing, dev::Device=CPU())
  fullsteppername = Symbol(stepper, :TimeStepper)
  if isexplicit(stepper)
    return eval(Expr(:call, fullsteppername, eq, dev))
  else
    return eval(Expr(:call, fullsteppername, eq, dt, dev))
  end
end

# The following time-steppers are implemented below
#
#   * Forward Euler
#   * Filtered Forward Euler
#   * RK4
#   * Filtered RK4
#   * ETDRK4
#   * Filtered ETDRK4
#   * AB3
#   * Filtered AB3
#
# Explicit time-steppers are constricted with the signature
#   ts = ExplicitTimeStepper(eq::Equation)
#
# Implicit time-steppers are constricted with the signature
#   ts = ImplicitTimeStepper(eq::Equation, dt)

# --
# Forward Euler
# --

"""
    ForwardEulerTimeStepper(eq, dev)

Initialize a forward Euler timestepper.
"""
struct ForwardEulerTimeStepper{T} <: AbstractTimeStepper{T}
  N::T # Explicit linear and nonlinear terms
  ForwardEulerTimeStepper(N::T) where T = new{T}(0N)
end

ForwardEulerTimeStepper(eq::Equation, dev::Device=CPU()) = ForwardEulerTimeStepper(devzeros(dev, eq.T, eq.dims))

function stepforward!(sol, cl, ts::ForwardEulerTimeStepper, eq, v, p, g)
  eq.calcN!(ts.N, sol, cl.t, cl, v, p, g)
  @. sol += cl.dt*(eq.L*sol + ts.N)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

function stepforward!(sol::AbstractArray{T}, cl, ts::ForwardEulerTimeStepper, eq, v, p, g) where T<:AbstractArray
  eq.calcN!(ts.N, sol, cl.t, cl, v, p, g)
  for i = 1:length(sol)
    @views @. sol[i] += cl.dt*(eq.L[i]*sol[i] + ts.N[i])
  end
  cl.t += cl.dt
  cl.step += 1
  nothing
end

"""
    FilteredForwardEulerTimeStepper(eq, dev; filterkwargs...)

Construct a forward Euler timestepper with spectral filtering.
"""
struct FilteredForwardEulerTimeStepper{T,Tf} <: AbstractTimeStepper{T}
  N::T
  filter::Tf
end

function FilteredForwardEulerTimeStepper(eq::Equation, dev::Device=CPU(); filterkwargs...)
  filter = makefilter(eq; filterkwargs...)
  FilteredForwardEulerTimeStepper(devzeros(dev, eq.T, eq.dims), filter)
end

function stepforward!(sol, cl, ts::FilteredForwardEulerTimeStepper, eq, v, p, g)
  eq.calcN!(ts.N, sol, cl.t, cl, v, p, g)
  @. sol = ts.filter*(sol + cl.dt*(ts.N + eq.L*sol))
  cl.t += cl.dt
  cl.step += 1
  nothing
end

# --
# RK4
# --

"""
    RK4TimeStepper(eq, dev)

Construct a 4th-order Runge-Kutta time stepper.
"""
struct RK4TimeStepper{T} <: AbstractTimeStepper{T}
  sol₁::T
  RHS₁::T
  RHS₂::T
  RHS₃::T
  RHS₄::T
end

"""
    FilteredRK4TimeStepper(eq, dev; filterkwargs...)

Construct a 4th-order Runge-Kutta time stepper with spectral filtering for the equation `eq`.
"""
struct FilteredRK4TimeStepper{T,Tf} <: AbstractTimeStepper{T}
  sol₁::T
  RHS₁::T
  RHS₂::T
  RHS₃::T
  RHS₄::T
  filter::Tf
end

function RK4TimeStepper(eq::Equation, dev::Device=CPU())
  @devzeros typeof(dev) eq.T eq.dims N sol₁ RHS₁ RHS₂ RHS₃ RHS₄
  RK4TimeStepper(sol₁, RHS₁, RHS₂, RHS₃, RHS₄)
end

function FilteredRK4TimeStepper(eq::Equation, dev::Device=CPU(); filterkwargs...)
  ts = RK4TimeStepper(eq, dev)
  filter = makefilter(eq; filterkwargs...)
  FilteredRK4TimeStepper(getfield.(Ref(ts), fieldnames(typeof(ts)))..., filter)
end

function addlinearterm!(RHS, L, sol)
  @. RHS += L*sol
  nothing
end

function substepsol!(newsol, sol, RHS, dt)
  @. newsol = sol + dt*RHS
  nothing
end

function RK4substeps!(sol, cl, ts, eq, v, p, g, t, dt)
  # Substep 1
  eq.calcN!(ts.RHS₁, sol, t, cl, v, p, g)
  addlinearterm!(ts.RHS₁, eq.L, sol)
  # Substep 2
  substepsol!(ts.sol₁, sol, ts.RHS₁, dt/2)
  eq.calcN!(ts.RHS₂, ts.sol₁, t+dt/2, cl, v, p, g)
  addlinearterm!(ts.RHS₂, eq.L, ts.sol₁)
  # Substep 3
  substepsol!(ts.sol₁, sol, ts.RHS₂, dt/2)
  eq.calcN!(ts.RHS₃, ts.sol₁, t+dt/2, cl, v, p, g)
  addlinearterm!(ts.RHS₃, eq.L, ts.sol₁)
  # Substep 4
  substepsol!(ts.sol₁, sol, ts.RHS₃, dt)
  eq.calcN!(ts.RHS₄, ts.sol₁, t+dt, cl, v, p, g)
  addlinearterm!(ts.RHS₄, eq.L, ts.sol₁)
  nothing
end

function RK4substeps!(sol::AbstractArray{T}, cl, ts, eq, v, p, g, t, dt) where T<:AbstractArray
  # Substep 1
  eq.calcN!(ts.RHS₁, sol, t, cl, v, p, g)
  addlinearterm!.(ts.RHS₁, eq.L, sol)
  # Substep 2
  substepsol!.(ts.sol₁, sol, ts.RHS₁, dt/2)
  eq.calcN!(ts.RHS₂, ts.sol₁, t+dt/2, cl, v, p, g)
  addlinearterm!.(ts.RHS₂, eq.L, ts.sol₁)
  # Substep 3
  substepsol!.(ts.sol₁, sol, ts.RHS₂, dt/2)
  eq.calcN!(ts.RHS₃, ts.sol₁, t+dt/2, cl, v, p, g)
  addlinearterm!.(ts.RHS₃, eq.L, ts.sol₁)
  # Substep 4
  substepsol!.(ts.sol₁, sol, ts.RHS₃, dt)
  eq.calcN!(ts.RHS₄, ts.sol₁, t+dt, cl, v, p, g)
  addlinearterm!.(ts.RHS₄, eq.L, ts.sol₁)
  nothing
end

function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, dt)
  @. sol += dt*(RHS₁ / 6 + RHS₂ / 3  + RHS₃ / 3 + RHS₄ / 6)
  nothing
end

function RK4update!(sol, RHS₁, RHS₂, RHS₃, RHS₄, filter, dt)
  @. sol = filter * (sol + dt*(RHS₁ / 6 + RHS₂ / 3  + RHS₃ / 3 + RHS₄ / 6))
  nothing
end

function stepforward!(sol, cl, ts::RK4TimeStepper, eq, v, p, g)
  RK4substeps!(sol, cl, ts, eq, v, p, g, cl.t, cl.dt)
  RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, cl.dt)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

function stepforward!(sol::AbstractArray{T}, cl, ts::RK4TimeStepper, eq, v, p, g) where T<:AbstractArray
  RK4substeps!(sol, cl, ts, eq, v, p, g, cl.t, cl.dt)
  RK4update!.(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, cl.dt)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

function stepforward!(sol, cl, ts::FilteredRK4TimeStepper, eq, v, p, g)
  RK4substeps!(sol, cl, ts, eq, v, p, g, cl.t, cl.dt)
  RK4update!(sol, ts.RHS₁, ts.RHS₂, ts.RHS₃, ts.RHS₄, ts.filter, cl.dt)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

# ------
# ETDRK4
# ------

"""
    ETDRK4TimeStepper(eq, dt)

Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper. The Rolls Royce of timestepping.
"""
struct ETDRK4TimeStepper{T,TL} <: AbstractTimeStepper{T}
  # ETDRK4 coefficents
  ζ::TL
  α::TL
  β::TL
  Γ::TL
  expLdt::TL
  expLdt2::TL
  sol₁::T
  sol₂::T
  N₁::T
  N₂::T
  N₃::T
  N₄::T
end

"""
    FilteredETDRK4TimeStepper(eq, dt; filterkwargs...)

Construct a 4th-order exponential-time-differencing Runge-Kutta time stepper with spectral filtering.
"""
struct FilteredETDRK4TimeStepper{T,TL,Tf} <: AbstractTimeStepper{T}
  # ETDRK4 coefficents:
  ζ::TL
  α::TL
  β::TL
  Γ::TL
  expLdt::TL
  expLdt2::TL
  sol₁::T
  sol₂::T
  N₁::T
  N₂::T
  N₃::T
  N₄::T
  filter::Tf
end

function ETDRK4TimeStepper(eq::Equation, dt, dev::Device=CPU())
  dt = fltype(eq.T)(dt) # ensure dt is correct type.
  expLdt, expLdt2 = getexpLs(dt, eq)
  ζ, α, β, Γ = getetdcoeffs(dt, eq.L)
  @devzeros typeof(dev) eq.T eq.dims sol₁ sol₂ N₁ N₂ N₃ N₄
  ETDRK4TimeStepper(ζ, α, β, Γ, expLdt, expLdt2, sol₁, sol₂, N₁, N₂, N₃, N₄)
end

function FilteredETDRK4TimeStepper(eq::Equation, dt, dev::Device=CPU(); filterkwargs...)
  ts = ETDRK4TimeStepper(eq, dt, dev)
  filter = makefilter(eq; filterkwargs...)
  FilteredETDRK4TimeStepper(getfield.(Ref(ts), fieldnames(typeof(ts)))..., filter)
end

function ETDRK4update!(sol, expLdt, α, β, Γ, N₁, N₂, N₃, N₄)
  @. sol = (expLdt*sol +  α * N₁
                       + 2β * (N₂ + N₃)
                       +  Γ * N₄ )
  nothing
end

function ETDRK4update!(sol, ts, filter)
  @. sol = filter*(ts.expLdt*sol +   ts.α * ts.N₁
                                 + 2*ts.β * (ts.N₂ + ts.N₃)
                                 +   ts.Γ * ts.N₄ )
  nothing
end

function ETDRK4substep12!(sol₁, expLdt2, sol, ζ, N)
  @. sol₁ = expLdt2*sol + ζ*N
  nothing
end

function ETDRK4substep3!(sol₂, expLdt2, sol₁, ζ, N₁, N₃)
  @. sol₂ = expLdt2*sol₁ + ζ*(2N₃ - N₁)
  nothing
end

function ETDRK4substeps!(sol, cl, ts, eq, v, p, g)
  # Substep 1
  eq.calcN!(ts.N₁, sol, cl.t, cl, v, p, g)
  ETDRK4substep12!(ts.sol₁, ts.expLdt2, sol, ts.ζ, ts.N₁)
  @. ts.sol₁ = ts.expLdt2*sol + ts.ζ*ts.N₁
  # Substep 2
  t2 = cl.t + cl.dt/2
  eq.calcN!(ts.N₂, ts.sol₁, t2, cl, v, p, g)
  ETDRK4substep12!(ts.sol₂, ts.expLdt2, sol, ts.ζ, ts.N₂)
  # Substep 3
  eq.calcN!(ts.N₃, ts.sol₂, t2, cl, v, p, g)
  ETDRK4substep3!(ts.sol₂, ts.expLdt2, ts.sol₁, ts.ζ, ts.N₁, ts.N₃)
  # Substep 4
  t3 = cl.t + cl.dt
  eq.calcN!(ts.N₄, ts.sol₂, t3, cl, v, p, g)
  nothing
end

function ETDRK4substeps!(sol::AbstractArray{T}, cl, ts, eq, v, p, g) where T<:AbstractArray
  # Substep 1
  eq.calcN!(ts.N₁, sol, cl.t, cl, v, p, g)
  ETDRK4substep12!.(ts.sol₁, ts.expLdt2, sol, ts.ζ, ts.N₁)
  @. ts.sol₁ = ts.expLdt2*sol + ts.ζ*ts.N₁
  # Substep 2
  t2 = cl.t + cl.dt/2
  eq.calcN!(ts.N₂, ts.sol₁, t2, cl, v, p, g)
  ETDRK4substep12!.(ts.sol₂, ts.expLdt2, sol, ts.ζ, ts.N₂)
  # Substep 3
  eq.calcN!(ts.N₃, ts.sol₂, t2, cl, v, p, g)
  ETDRK4substep3!.(ts.sol₂, ts.expLdt2, ts.sol₁, ts.ζ, ts.N₁, ts.N₃)
  # Substep 4
  t3 = cl.t + cl.dt
  eq.calcN!(ts.N₄, ts.sol₂, t3, cl, v, p, g)
  nothing
end

function stepforward!(sol, cl, ts::ETDRK4TimeStepper, eq, v, p, g)
  ETDRK4substeps!(sol, cl, ts, eq, v, p, g)
  ETDRK4update!(sol, ts.expLdt, ts.α, ts.β, ts.Γ, ts.N₁, ts.N₂, ts.N₃, ts.N₄)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

function stepforward!(sol::AbstractArray{T}, cl, ts::ETDRK4TimeStepper, eq, v, p, g) where T<:AbstractArray
  ETDRK4substeps!(sol, cl, ts, eq, v, p, g)
  ETDRK4update!.(sol, ts.expLdt, ts.α, ts.β, ts.Γ, ts.N₁, ts.N₂, ts.N₃, ts.N₄)
  cl.t += cl.dt
  cl.step += 1
  nothing
end

function stepforward!(sol, cl, ts::FilteredETDRK4TimeStepper, eq, v, p, g)
  ETDRK4substeps!(sol, cl, ts, eq, v, p, g)
  ETDRK4update!(sol, ts, ts.filter) # update
  cl.t += cl.dt
  cl.step += 1
  nothing
end

# ---
# AB3
# ---

"""
    AB3TimeStepper(eq, dev)

Construct a 3rd order Adams-Bashforth time stepper.
"""
const ab3h1 = 23/12
const ab3h2 = 16/12
const ab3h3 = 5/12

struct AB3TimeStepper{T} <: AbstractTimeStepper{T}
  RHS::T
  RHS₋₁::T
  RHS₋₂::T
end

function AB3TimeStepper(eq::Equation, dev::Device=CPU())
  @devzeros typeof(dev) eq.T eq.dims RHS RHS₋₁ RHS₋₂
  AB3TimeStepper(RHS, RHS₋₁, RHS₋₂)
end

"""
    FilteredAB3TimeStepper(eq, dev; filterkwargs...)_

Construct a 3rd order Adams-Bashforth time stepper with spectral filtering.
"""
struct FilteredAB3TimeStepper{T,Tf} <: AbstractTimeStepper{T}
  RHS::T
  RHS₋₁::T
  RHS₋₂::T
  filter::Tf
end

function FilteredAB3TimeStepper(eq::Equation, dev::Device=CPU(); filterkwargs...)
  ts = AB3TimeStepper(eq, dev)
  filter = makefilter(eq; filterkwargs...)
  FilteredAB3TimeStepper(getfield.(Ref(ts), fieldnames(typeof(ts)))..., filter)
end

function AB3update!(sol, ts, cl)
  if cl.step < 3  # forward Euler steps to initialize AB3
    @. sol += cl.dt*ts.RHS    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. sol += cl.dt*(ab3h1*ts.RHS - ab3h2*ts.RHS₋₁ + ab3h3*ts.RHS₋₂)
  end
  nothing
end

function AB3update!(sol, ts::FilteredAB3TimeStepper, cl)
  if cl.step < 3  # forward Euler steps to initialize AB3
    @. sol = ts.filter*(sol + cl.dt*ts.RHS)    # Update
  else   # Otherwise, stepforward with 3rd order Adams Bashforth:
    @. sol = ts.filter*(sol + cl.dt*(ab3h1*ts.RHS - ab3h2*ts.RHS₋₁ + ab3h3*ts.RHS₋₂))
  end
  nothing
end

function stepforward!(sol, cl, ts::AB3TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS, sol, cl.t, cl, v, p, g)
  addlinearterm!(ts.RHS, eq.L, sol)
  AB3update!(sol, ts, cl)
  cl.t += cl.dt
  cl.step += 1
  @. ts.RHS₋₂ = ts.RHS₋₁          # Store
  @. ts.RHS₋₁ = ts.RHS            # ... previous values of RHS
  nothing
end

function stepforward!(sol::AbstractArray{T}, cl, ts::AB3TimeStepper, eq, v, p, g) where T<:AbstractArray
  eq.calcN!(ts.RHS, sol, cl.t, cl, v, p, g)
  addlinearterm!.(ts.RHS, eq.L, sol)
  AB3update!(sol, ts, cl)
  cl.t += cl.dt
  cl.step += 1
  @. ts.RHS₋₂ = ts.RHS₋₁          # Store
  @. ts.RHS₋₁ = ts.RHS            # ... previous values of RHS
  nothing
end

function stepforward!(sol, cl, ts::FilteredAB3TimeStepper, eq, v, p, g)
  eq.calcN!(ts.RHS, sol, cl.t, cl, v, p, g)
  addlinearterm!(ts.RHS, eq.L, sol)
  AB3update!(sol, ts, cl)
  cl.t += cl.dt
  cl.step += 1
  @. ts.RHS₋₂ = ts.RHS₋₁          # Store
  @. ts.RHS₋₁ = ts.RHS            # ... previous values of RHS
  nothing
end

# --
# Timestepper utils
# --

function getexpLs(dt, eq::Equation{T,TL,Tg}) where {T,TL<:Array{TLi},Tg} where TLi<:AbstractArray
  expLdt = [ @.(exp(dt*L)) for L in eq.L ]
  expLdt2 = [ @.(exp(dt*L/2)) for L in eq.L ]
  expLdt, expLdt2
end

function getexpLs(dt, eq)
  expLdt  = @. exp(dt*eq.L)
  expLdt2 = @. exp(dt*eq.L/2)
  expLdt, expLdt2
end

function getetdcoeffs(dt, L::AbstractArray{T}; kwargs...) where T<:AbstractArray
  ζαβΓ = [ getetdcoeffs(dt, Li) for Li in L ]
  ζ = map(x->x[1], ζαβΓ)
  α = map(x->x[2], ζαβΓ)
  β = map(x->x[3], ζαβΓ)
  Γ = map(x->x[4], ζαβΓ)
  ζ, α, β, Γ
end

"""
    getetdcoeffs(dt, L; ncirc=32, rcirc=1)

Calculate ETDRK4 coefficients associated with the (diagonal) linear coefficient
L by integrating over a small circle in complex space.
"""
function getetdcoeffs(dt, L; ncirc=32, rcirc=1)

  shape = Tuple(cat(ncirc, ones(Int, ndims(L)), dims=1))

  circ = zeros(Complex{Float64}, shape) # use double precision for this calculation
  circ .= rcirc * exp.(2π*im/ncirc*(0.5:1:(ncirc-0.5)))
  circ = permutedims(circ, ndims(circ):-1:1)

  zc = dt*L .+ circ
  M = ndims(L)+1

  # Four coefficients: ζ, α, β, Γ
  ζc = @.             ( exp(zc/2)-1 ) / zc
  αc = @. ( -4 - zc + exp(zc)*(4 - 3zc + zc^2) ) / zc^3
  βc = @.    ( 2  + zc + exp(zc)*(-2 + zc) ) / zc^3
  Γc = @. ( -4 - 3zc - zc^2 + exp(zc)*(4 - zc) ) / zc^3

  ζ = dt*dropdims(mean(ζc, dims=M), dims=M)
  α = dt*dropdims(mean(αc, dims=M), dims=M)
  β = dt*dropdims(mean(βc, dims=M), dims=M)
  Γ = dt*dropdims(mean(Γc, dims=M), dims=M)

  if eltype(L) <: Real # this is conservative, but unclear if necessary
    ζ = real.(ζ)
    α = real.(α)
    β = real.(β)
    Γ = real.(Γ)
  end

  ζ, α, β, Γ
end


text:
using Test

using NeXLCore

@testset "Stoichiometry" begin

    # Check against DTSA-II
    @test isapprox(
        obystoichiometry(
            n"Mg" => 0.116567,
            n"Al" => 0.0490615,
            n"Si" => 0.211982,
            n"Ca" => 0.10899,
            n"Fe" => 0.0774196,
        ),
        0.4276,
        atol = 0.0001,
    )

    sio2 = asoxide(n"Si")
    @test name(sio2) == "SiO₂"
    asio2 = atomicfraction("SiO2", n"Si" => 1, n"O" => 2)
    asio2 = atomicfraction("SiO₂", n"Si" => 1, n"O" => 2)
    @test isapprox(sio2[n"Si"], asio2[n"Si"], atol = 1.0e-6)
    @test isapprox(sio2[n"O"], asio2[n"O"], atol = 1.0e-6)

    mix = sum(asoxide(
        n"Mg" => 0.116567,
        n"Al" => 0.0490615,
        n"Si" => 0.211982,
        n"Ca" => 0.10899,
        n"Fe" => 0.0774196))

    @test isapprox(mix[n"O"], 0.4276, atol = 0.0001)
    @test isapprox(mix[n"Al"], 0.0490615, atol = 0.0001)
    @test isapprox(mix[n"Fe"], 0.0774196, atol = 0.0001)

    ox = sum(asoxide(massfraction(delete(srm470_k412, n"O"))), name="K412")
    @test all(isapprox(ox[elm], srm470_k412[elm], atol=1.0e-8) for elm in keys(srm470_k412))

    ox = sum(asoxide(massfraction(delete(srm470_k411, [ n"O" ]))), name="K411")
    @test all(isapprox(ox[elm], srm470_k411[elm], atol=1.0e-8) for elm in keys(srm470_k411))

    m=sum(asoxide(n"Al"=>0.1029, n"Na"=>0.0877, n"Si"=>0.3213), name="Albite", density=5.0, description="a common mineral", pedigree="Natural")
    @test isapprox(m[n"O"], 0.4881, atol=0.0001)
    @test isapprox(m[n"Si"], 0.3212, atol=0.0001)
    @test isapprox(m[n"Na"], 0.0877, atol=0.0001)
    @test m[:Density] == 5.0
    @test m[:Description] == "a common mineral"
    @test m[:Pedigree] == "Natural"
    @test name(m) == "Albite"

    amix =
        0.4 * atomicfraction("SiO2", n"Si" => 1, n"O" => 2) +
        0.6 * atomicfraction("Al2O3", n"Al" => 2, n"O" => 3)
    @test isapprox(amix[n"O"], 0.4955, atol = 0.0001)
    @test isapprox(amix[n"Al"], 0.3176, atol = 0.0001)
    @test isapprox(amix[n"Si"], 0.1870, atol = 0.0001)

    # Check against DTSA-II
    @test isapprox(amix[n"Si"], 0.1870, atol = 0.0001)
    @test isapprox(amix[n"Al"], 0.3176, atol = 0.0001)
    @test isapprox(amix[n"O"], 0.4955, atol = 0.0001)
end


text:
"""
    dict_list(c::Dict)
Expand the dictionary `c` into a vector of dictionaries.
Each entry has a unique combination from the product of the `Vector`
values of the dictionary while the non-`Vector` values are kept constant
for all possibilities. The keys of the entries are the same.

Whether the values of `c` are iterable or not is of no concern;
the function considers as "iterable" only subtypes of `Vector`.

To restrict some values in the dictionary so that they only appear in the
resulting dictionaries, if a certain condition is met, the macro
[`@onlyif`](@ref) can be used on those values.

Use the function [`dict_list_count`](@ref) to get the number of
dictionaries that `dict_list` will produce.

## Examples
```julia
julia> c = Dict(:a => [1, 2], :b => 4);

julia> dict_list(c)
2-element Array{Dict{Symbol,Int64},1}:
 Dict(:a=>1,:b=>4)
 Dict(:a=>2,:b=>4)

julia> c[:model] = "linear"; c[:run] = ["bi", "tri"];

julia> dict_list(c)
4-element Array{Dict{Symbol,Any},1}:
 Dict(:a=>1,:b=>4,:run=>"bi",:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:model=>"linear")

julia> c[:e] = [[1, 2], [3, 5]];

julia> dict_list(c)
8-element Array{Dict{Symbol,Any},1}:
 Dict(:a=>1,:b=>4,:run=>"bi",:e=>[1, 2],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:e=>[1, 2],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:e=>[1, 2],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:e=>[1, 2],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"bi",:e=>[3, 5],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"bi",:e=>[3, 5],:model=>"linear")
 Dict(:a=>1,:b=>4,:run=>"tri",:e=>[3, 5],:model=>"linear")
 Dict(:a=>2,:b=>4,:run=>"tri",:e=>[3, 5],:model=>"linear")
```
"""
function dict_list(c::Dict)
    if contains_partially_restricted(c)
        # The method for generating the restricted parameter set is as follows:
        # 1. Create an array of trial combinations containing all possible
        # combinations of all parameters.
        # 2. For each solution, remove all parameters where the conditions
        # functions returns false
        # 3. Do (2) until the length of the obtained parameter set stops
        # changing.
        # 4. Replace all `DependentParameter` types by their respective values
        # 5. This gives a parameter dict with a valid combination
        # 6. From the resulting list of valid combinations remove the
        # duplicates.
        # 7. Remove solutions which are only a subset of others.
        parameter_sets = Set(map(_dict_list(c)) do trial
                                 n = length(trial)
                                 for i in 1:100_000
                                     for key in keys(trial)
                                         val = trial[key]
                                         if val isa DependentParameter && !val.condition(c,trial)
                                             delete!(trial,key)
                                         end
                                     end
                                     length(trial) == n && break
                                     n = length(trial)
                                     i == 100_000 && error("There are too many parameters with a serial dependency. The limit is set to 100000.")
                                 end
                                 Dict([k=>lookup_candidate(c,trial,k) for k in keys(trial)])
                             end)
        return collect(filter(parameter_sets) do trial
            !is_solution_subset_of_existing(trial, parameter_sets)
        end)
    end
    return _dict_list(c)
end

function is_solution_subset_of_existing(trial, trial_solutions)
    ks = Set(keys(trial))
    for _trial in trial_solutions
        trial == _trial && continue
        ks ⊆ Set(keys(_trial)) || continue
        all(trial[k] == _trial[k] for k in ks) && return true
    end
    return false
end

function _dict_list(c::Dict)
    iterable_fields = filter(k -> typeof(c[k]) <: Vector, keys(c))
    non_iterables = setdiff(keys(c), iterable_fields)

    iterable_dict = Dict(iterable_fields .=> getindex.(Ref(c), iterable_fields))
    non_iterable_dict = Dict(non_iterables .=> getindex.(Ref(c), non_iterables))

    vec(
        map(Iterators.product(values(iterable_dict)...)) do vals
            dd = Dict(keys(iterable_dict) .=> vals)
            if isempty(non_iterable_dict)
                dd
            elseif isempty(iterable_dict)
                non_iterable_dict
            else
                merge(non_iterable_dict, dd)
            end
        end
    )
end

"""
    dict_list_count(c) -> N
Return the number of dictionaries that will be created by
calling `dict_list(c)`.
"""
function dict_list_count(c)
    if DependentParameter in eltype.(values(c))
        return length(dict_list(c))
    end
    iterable_fields = filter(k -> typeof(c[k]) <: Vector, keys(c))
    prod(length(c[i]) for i in iterable_fields)
end

# Basis taken from https://github.com/FluxML/MacroTools.jl
walk(x, inner, outer, ex) = outer(x, ex)
walk(x::Expr, inner, outer, ex) = outer(Expr(x.head, map(y->inner(y,x), x.args)...), ex)
postwalk(f, x, ex=:()) = walk(x, (x,y) -> postwalk(f, x, y), f, ex)

export @onlyif

struct DependentParameter{T}
    value::T
    condition::Function
end

contains_partially_restricted(d::Dict) = any(contains_partially_restricted,values(d))
contains_partially_restricted(d::Vector) = any(contains_partially_restricted,d)
contains_partially_restricted(::DependentParameter) = true
contains_partially_restricted(::Any) = false

function toDependentParameter(value::T,condition) where T
    if T <: Vector
        return DependentParameter.(value,Ref(condition))
    end
    return DependentParameter(value,condition)
end

struct KeyDeletedFromDictError <: Exception end

"""
    lookup_candidate(original_dict, d, name)
If `name` is a key name from `original_dict` either return it's value from `d`
or throw a `KeyDeletedFromDictError` error if it's not in `d`.
"""
function lookup_candidate(original_dict,d, name)
    if name in keys(original_dict)
        name in keys(d) || throw(KeyDeletedFromDictError())
        if d[name] isa DependentParameter
            return d[name].value
        end
        return d[name]
    end
    return name
end

"""
    @onlyif(ex, value)
Tag `value` to only appear in a dictionary created with [`dict_list`](@ref) if
the Julia expression `ex` (see below) is evaluated as true.  If `value` is a subtype of
`Vector`, `@onlyif` is applied to each entry.
Since `@onlyif` is applied to a value and not to a dictionary key, it is
possible to restrict only some of the values of a vector. This means that based
on `ex` the number of options for a particular key varies.

Within `ex` it is possible to extract values of the dictionary passed to
[`dict_list`](@ref) by a shorthand notation where only the key must be
provided.  For example `ex = :(:N == 1)` is tranformed in the call
`dict_list(d)` to an expression analogous to `:(d[:N] == 1)` by using the
function `lookup_candidate`.  This is supported for `Symbol` and `String` keys.

## Examples
```julia
julia> d = Dict(:a => [1, 2], :b => 4, :c => @onlyif(:a == 1, [10, 11]));

julia> dict_list(d) # only in case `:a` is `1` the dictionary will get key `:c`
3-element Array{Dict{Symbol,Int64},1}:
 Dict(:a => 1,:b => 4,:c => 10)
 Dict(:a => 1,:b => 4,:c => 11)
 Dict(:a => 2,:b => 4)

 julia> d = Dict(:a => [1, 2], :b => 4, :c => [10, @onlyif(:a == 1, 11)]);

 julia> dict_list(d) # only in case `:a` is `1` the dictionary will get extra value `11` for key `:c`
 3-element Array{Dict{Symbol,Int64},1}:
 Dict(:a => 1,:b => 4,:c => 10)
 Dict(:a => 1,:b => 4,:c => 11)
 Dict(:a => 2,:b => 4,:c => 10)
```
 See the [Defining parameter sets with restrictions](@ref) section for more examples.
 """
 macro onlyif(ex, value)
     pd = gensym()
     original_dict = gensym()
     anon_function = quote
         ($(original_dict),$pd)->try
             $(ex)
         catch e
             e isa DrWatson.KeyDeletedFromDictError && return false
             rethrow(e)
         end
     end
     condition = postwalk(anon_function) do x, parent
         # Check if the parent expression in the postwalk is a dot expression
         # like Foo.Bar. In this case Bar is a QuoteNode, however it should not
         # be wrapped by the lookup_candidate function, because without Foo it's
         # not valid.
         if (x isa QuoteNode && !(parent isa Expr && parent.head == :.)) || x isa String
             return :(DrWatson.lookup_candidate($(original_dict),$pd, $x))
         end
         return x
     end
     :(toDependentParameter($(esc(value)),$(esc(condition))))
 end


text:
"""
"""
module StateEstimation

export
    # Top-level abstract types
    Estimator, SequentialEstimator, BatchEstimator,
    # States
    AbstractState, AbstractAbsoluteState, AbstractUncertainState,
    DiscreteState, UncertainDiscreteState, ContinuousState,
    UncertainContinuousState, make_uncertain, make_absolute,
    sample, distance, mahalanobis, expand_state, reduce_state,
    # Systems
    AbstractSystem, LinearSystem, NonlinearSystem, predict, predict!,
    state_transition_matrix,
    # Observers
    AbstractObserver, LinearObserver, NonlinearObserver, observable,
    # Data archiving
    EstimatorHistory, plot_state_history, plot_residuals,
    # Kalman filters
    AbstractKalmanFilter, KalmanFilter, ExtendedKalmanFilter, simulate,
    process!, kalman_filter, kalman_filter!, extended_kalman_filter,
    extended_kalman_filter!,
    # Unscented Kalman filter
    UnscentedTransform, compute_weights, compute_sigma_points, transform,
    trasnform!, augment, UnscentedKalmanFilter, unscented_kalman_filter,
    unscented_kalman_filter!,
    # Least squares estimation
    LeastSquaresEstimator, NonlinearLeastSquaresEstimator, add!, solve, solve!,
    # Multi-target filtering
    MultiTargetFilter, NearestNeighborMTF,
    # Simulation
    SingleStateSimulator, MultiStateSimulator, make_simulator

abstract type Estimator{T,S} end
abstract type SequentialEstimator{T,S} <: Estimator{T,S} end
abstract type BatchEstimator{T,S} <: Estimator{T,S} end

include("states.jl")
include("systems.jl")
include("observers.jl")
include("archive.jl")

include("kalman_filter.jl")
include("extended_kalman_filter.jl")
include("unscented_transform.jl")
#include("unscented_kalman_filter.jl")
include("multi_target_filter.jl")

include("least_squares.jl")

include("simulator.jl")

end


text:
module VectorCalculus

using ForwardDiff
using LinearAlgebra

import LinearAlgebra: dot, cross, norm, div
import Base: div
import Base: div, getindex, setindex!, convert, size, similar, parent, axes, IndexStyle

using Base: @propagate_inbounds

using Combinatorics: levicivita

include("coordinates.jl")
include("coordinate_vectors.jl")
include("algebraic_operations.jl")
include("differential_operations.jl")

export Cartesian, Cylindrical, Spherical
export CartesianCoordinate, CylindricalCoordinate, SphericalCoordinate
export Coordinate, CoordinateVector
export CovariantVector, ContravariantVector
export metric
export covariant_basis, covariant_metric
export contravariant_basis, contravariant_metric
export unit_basis_components
export grad, laplacian, curl, vector_laplacian

end


text:
#!/usr/bin/julia

# Date: 21 August 2016
# https://github.com/trizen

# An efficient algorithm for computing large Fibonacci numbers, modulus some n.
# Algorithm from: http://codeforces.com/blog/entry/14516

function fibmod(n::Int64, mod::Int64)

    n <= 1 && return n
    cache = Dict{Int64, Int64}()

    function f(n::Int64)

        n <= 1           && return 1
        haskey(cache, n) && return cache[n]

        k = div(n, 2)

        cache[n] =
            if n % 2 == 0
                (f(k) * f(k  ) + f(k-1) * f(k-1)) % mod
            else
                (f(k) * f(k+1) + f(k-1) * f(k  )) % mod
            end
    end

    f(n-1)
end

println(fibmod(1000, 10^4))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
# Analyze convergence of Monte Carlo tree search for various horizon lenghts

using SubmodularMaximization
using PyPlot
using Statistics
using Base.Threads
using Base.Iterators
using JLD2

close()

data_file = "./data/mcvs_convergence_data.jld2"
reprocess = false

num_trials = 100
num_mcts_samples = 2 .^ (4:16)
horizons = 1:6

grid_size = 10
sensor = RangingSensor()

# place the target in the center of the grid
grid = Grid(grid_size, grid_size)
target_state = (5,5)
robot_state = (1,1)

# Create some uncertainty. Two steps is probably a reasonable approximation of
# the kind of target uncertainty that we will encounter
histogram_filter = Filter(grid, target_state)
process_update!(histogram_filter, transition_matrix(grid))

function get_data()
  data = Array{Float64}(undef, num_trials,
                               length(num_mcts_samples),
                               length(horizons))

  if !isfile(data_file) || reprocess
    for (kk, horizon) in Iterators.reverse(enumerate(horizons)),
      (jj, num_mcts_samples) in Iterators.reverse(enumerate(num_mcts_samples))
      println("Horizon: ", horizon, " Num. samples: ", num_mcts_samples)

      problem = SingleRobotTargetTrackingProblem(grid, sensor, horizon,
                                                 [histogram_filter])

      @time for (block_num, block) in enumerate(partition(1:num_trials,
                                                          Threads.nthreads()))
        println("Block num: ", block_num)
        # Ensure that thread execution does not overlap
        @threads for ii in block
          # solve via mcts
          trajectory = solve_single_robot(problem, robot_state,
                                          n_iterations = num_mcts_samples).trajectory

          # evaluate mutual information using the default number of samples (which
          # should be quite accurate)
          reward = finite_horizon_information(grid, histogram_filter, sensor,
                                              trajectory).reward

          data[ii, jj, kk] = reward
        end
      end
    end

    @save data_file data
  else
    @load data_file data
  end

  data
end

data = get_data()

mean_rewards = reshape(mean(data, dims=1), length.((num_mcts_samples,
                                                    horizons)))

for (ii, horizon) in enumerate(horizons)
  horizon_data = mean_rewards[:,ii]
  maximum_reward = maximum(horizon_data)
  fraction_rewards =  horizon_data ./ maximum_reward

  println("Horizon: ", horizon, " Max reward: ", maximum_reward)

  label = string("Horz.: ", horizon)
  loglog(num_mcts_samples, fraction_rewards, label=label)
end

title("MCTS convergence")
xlabel("Num rollouts")
ylabel("Frac. of max. reward")
legend(loc="upper left")

save_fig("fig", "mcts_convergence")

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
"Construction of the little-d Wigner function, and the counting modes for rotations"

using WignerFunctions

export tp, ld, λ, μ

"`tp(l, m₁, m₂)` is the redefinition of the tupled of indexes for the Wigner little-_d_ matrix"
tp(l, m₁, m₂) = WignerFunctions.Index((l, m₁, m₂))

"`ld(β, l, m₁, m₂) is the redefinition of the Wigner little-_d_ matrix`"
ld(β, l, m₁, m₂) = WignerFunctions.naive(β, tp(l, m₁, m₂))

"`λ(n, j₁, j₂)` is a romboid counter for the little-_d_ matrix"
function λ(n, j₁, j₂)
    if 0 ≤ n ≤ 2 * j₂
        return (1/2) * n
    elseif 2 * j₂ < n ≤ 2 * j₁
        return j₂
    elseif 2 * j₁ ≤ n ≤ 2 * (j₁ + j₂)
        return j₁ + j₂ - (1/2) * n
    end
end

"`μ(n, m, j₁, j₂)` is a romboid counter for the little-_d_ matrix. We have the equivalence m = n₁ - n₂, and n = n₁ + n₂"
function μ(n, m, j₁, j₂)
    if 0 ≤ n ≤ 2 * j₂
        return (1/2) * m
    elseif 2 * j₂ < n ≤ 2 * j₁
        return j₂ + (1/2) * (m - n)
    elseif 2 * j₁ ≤ n ≤ 2 * (j₁ + j₂)
        return (1/2) * m - j₁ + j₂
    end
end

"`expf(σ, m)` is a index-dependent phase"
expf(σ, n₁, n₂) = exp(σ * 1im * (1/4) * π * m)

function Ld(j₁, j₂)
    v₀ = zeros((2 * j₁ + 1, 2 * j₂ + 1))
    for n in 0:2 * (j₁ + j₂)
        for m in -n:2:n
            v₀[n + 1, m + n + 1] = ld(π/2, λ(n, j₁, j₂), μ(n, m, j₁, j₂), μ(n, m, j₁, j₂))
        end
    end
    return v₀
end


text:
using ImageGather, Test

@test true


text: Dict{String,Any}("people" => Any[Dict{String,Any}("first_name" => Dict{String,Any}("value" => "Bruce","type" => "string"),"last_name" => Dict{String,Any}("value" => "Springsteen","type" => "string")), Dict{String,Any}("first_name" => Dict{String,Any}("value" => "Eric","type" => "string"),"last_name" => Dict{String,Any}("value" => "Clapton","type" => "string")), Dict{String,Any}("first_name" => Dict{String,Any}("value" => "Bob","type" => "string"),"last_name" => Dict{String,Any}("value" => "Seger","type" => "string"))])

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
using Access
using Test

fn = cp(
    joinpath(@__DIR__, "TestDatabase.accdb"),
    joinpath(mktempdir(; cleanup=false), "Test.accdb"),
)
chmod(fn, 0o777) # Make sure file is writable for tests

@testset "Read data from db" begin
    @test isfile(fn)
    df = Access.query(fn, "SELECT * FROM Person WHERE BirthYear > 1980")

    @test size(df, 1) == 1 # One Row
    @test df.FirstName[1] == "Slim"
    @test df.LastName[1] == "Shady"
    @test df.ID[1] == 3
end

@testset "Modify db data" begin

    # Read table content to DataFrame
    StoredValues = Access.query(fn, "SELECT * FROM StoredValues;")

    # Delete everything in table
    Access.query!(fn, "DELETE * FROM StoredValues;")

    # Check table empty?
    shouldbeemtpy = Access.query(fn, "SELECT * FROM StoredValues;")
    @test size(shouldbeemtpy, 1) == 0

    # Re-insert data from DataFrame
    Access.@insertdf(fn, StoredValues)
    shouldnotbeemtpy = Access.query(fn, "SELECT * FROM StoredValues;")
    @test size(shouldnotbeemtpy, 1) > 0

    # Check table content
    @test shouldnotbeemtpy[6, 4] == "Mashall Mathers"
    @test shouldnotbeemtpy[2, 3] ≈ 3.3

    # Alternatively, specify table explicitly:
    Access.query!(fn, "DELETE * FROM StoredValues;")
    Access.insertdf(fn, StoredValues, "StoredValues")
    shouldnotbeemtpy = Access.query(fn, "SELECT * FROM StoredValues;")
    @test shouldnotbeemtpy[2, 3] ≈ 3.3
end

# TODO: Figure out why this (and cleanup) fails
# rm(dirname(fn);recursive=true)

max_stars_repo_name: hellemo/Access

max_issues_repo_name: hellemo/Access

max_forks_repo_name: hellemo/Access


text:
# HTML Output Utilities

using Base64: Base64EncodePipe
using Colors
using ImageTransformations
using DomTerm

make_newline() = make_div("clear:both", "")
function make_div(st, content::String; id="", class="")
    if !isempty(id) && !isempty(class)
        return "
"*content*"
"
    end
    if !isempty(id) && isempty(class)
        return "
"*content*"
"
    end
    if isempty(id) && !isempty(class)
        return "
"*content*"
"
    end
    if isempty(id) && isempty(class)
        return "
"*content*"
"
    end
end

function make_img(st, content::String; id="", class="")
    if !isempty(id) && !isempty(class)
        return ""
    end
    if !isempty(id) && isempty(class)
        return ""
    end
    if isempty(id) && !isempty(class)
        return ""
    end
    if isempty(id) && isempty(class)
        return ""
    end
end

DTPIOContext(io, KVs::Pair...) = IOContext(
    io,
    :limit=>true, :color=>true, :domterm=>true,
    KVs...
)

function cmatrix_to_pngstr(img::AbstractMatrix{<:Colorant}; ratio=0.5, maxlen=10^6)
    buf = IOBuffer()
    b64 = Base64EncodePipe(buf)
    while applicable(length, img) && length(img) > maxlen
        img = imresize(img; ratio)
    end
    show(DTPIOContext(b64), MIME("image/png"), img)
    png_str = String(take!(buf))
    close(b64)

    return png_str
end

function make_img(st, im::AbstractMatrix{<:Colorant}; id="", class="")
    return make_img(st, cmatrix_to_pngstr(im); id, class)
end

to_html_ele(x::AbstractMatrix{<:Colorant}) = make_img("", cmatrix_to_pngstr(im))
to_html_ele(x::String) = make_div("", x)  # FIXME

function hprint(args...)
    content = ""
    for arg in args
        content *= to_html_ele(arg)
    end

    if !isempty(content)
        hput(content)
    end
end

function hprintln(args...)
    !isempty(args) && hprint(args...)
    hput("

")
end

# `x`: an HTML element
# dt_addst(".centerizer * {margin: 0 auto}")
function centerize(x::String)
    return make_div("text-align:center", x; class="centerizer")
end


text:
module MadNLPGPU

import LinearAlgebra
# CUDA
import CUDA: CUBLAS, CUSOLVER, CuVector, CuMatrix, CuArray, toolkit_version, R_64F, has_cuda, @allowscalar, runtime_version
# Kernels
import KernelAbstractions: @kernel, @index, wait, Event
import CUDAKernels: CUDADevice

import MadNLP

import MadNLP:
    @kwdef, Logger, @debug, @warn, @error,
    AbstractOptions, AbstractLinearSolver, AbstractNLPModel, set_options!, MadNLPLapackCPU,
    SymbolicException,FactorizationException,SolveException,InertiaException,
    introduce, factorize!, solve!, improve!, is_inertia, inertia, tril_to_full!

include("kernels.jl")

if has_cuda()
    include("lapackgpu.jl")
    export MadNLPLapackGPU
end

end # module


text:
struct MockDofBasis{P} <: AbstractVector{PointValue}
  nodes::Vector{P}
end

function return_cache(b::MockDofBasis,field)
  return_cache(field,b.x)
end

@inline function evaluate!(cache,b::MockDofBasis,field)
  evaluate!(cache,field,b.x)
end

@inline Base.size(a::MockDofBasis) = (length(a.nodes),)
@inline Base.axes(a::MockDofBasis) = (axes(a.nodes,1),)
# @santiagobadia : Not sure we want to create the moment dofs
@inline Base.getindex(a::MockDofBasis,i::Integer) = PointValue(a.x[i])
@inline Base.IndexStyle(::MockDofBasis) = IndexLinear()


text:
##############################################################################
##
## Hamming
##
##############################################################################
function evaluate(dist::Hamming, s1::AbstractString, s2::AbstractString)
    current = abs(length(s2) - length(s1))
    for (ch1, ch2) in zip(s1, s2)
        current += ch1 != ch2
    end
    return current
end

##############################################################################
##
## Jaro
##
##############################################################################
"""
    Jaro()

Creates the Jaro metric

The Jaro distance is defined as

``1 - (m / |s1| + m / |s2| + (m - t) / m) / 3``

where ``m`` is the number of matching characters and
``t`` is half the number of transpositions.
"""
struct Jaro <: SemiMetric end

## http://alias-i.com/lingpipe/docs/api/com/aliasi/spell/JaroWinklerDistance.html
function evaluate(dist::Jaro, s1::AbstractString, s2::AbstractString)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    # if both are empty, m = 0 so should be 1.0 according to wikipedia. Add this line so that not the case
    len2 == 0 && return 0.0
    maxdist = max(0, div(len2, 2) - 1)
    flag = fill(false, len2)
    prevstate1 = firstindex(s1)
    i1_match = fill(prevstate1, len1)
    #  m counts number matching characters
    m = 0
    i1 = 1
    i2 = 1
    x1 = iterate(s1)
    x2 = iterate(s2)
    while x1 !== nothing
        ch1, state1 = x1
        if i2 <= i1 - maxdist - 1
            ch2, state2 = x2
            i2 += 1
            x2 = iterate(s2, state2)
        end
        i2curr = i2
        x2curr = x2
        while x2curr !== nothing
            i2curr > i1 + maxdist && break
            ch2, state2 = x2curr
            if (ch1 == ch2) && !flag[i2curr]
                m += 1
                flag[i2curr] = true
                i1_match[m] = prevstate1
                break
            end
            x2curr = iterate(s2, state2)
            i2curr += 1
        end
        x1 = iterate(s1, state1)
        i1 += 1
        prevstate1 = state1
    end
    m == 0 && return 1.0
    # t counts number of transpositions
    t = 0
    i1 = 0
    i2 = 0
    for ch2 in s2
        i2 += 1
        if flag[i2]
            i1 += 1
            t += ch2 != iterate(s1, i1_match[i1])[1]
        end
    end
    current = 1.0 - (m / len1 + m / len2 + (m - t/2) / m) / 3.0
    return current
end

##############################################################################
##
## Levenshtein
##
## Return max_dist +1 if distance higher than max_dist
## This makes it possible to differentiate distance equalt to max_dist vs strictly higher
## This is important for find_all
##
##############################################################################
"""
    Levenshtein()

Creates the Levenshtein metric

The Levenshtein distance is the minimum number of operations (consisting of insertions, deletions, substitutions of a single character) required to change one string into the other.
"""
struct Levenshtein <: SemiMetric end

## Source: http://blog.softwx.net/2014/12/optimizing-levenshtein-algorithm-in-c.html
function evaluate(dist::Levenshtein, s1::AbstractString, s2::AbstractString;
    max_dist = nothing)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    max_dist !== nothing && len2 - len1 > max_dist && return max_dist + 1
    # prefix common to both strings can be ignored
    k, x1, x2start = remove_prefix(s1, s2)
    x1 == nothing && return len2 - k
    # distance initialized to first row of matrix
    # => distance between "" and s2[1:i}
    v0 = collect(1:(len2 - k))
    current = 0
    i1 = 1
    while x1 !== nothing
        ch1, state1 = x1
        left = i1 - 1
        current = i1 - 1
        min_dist = i1 - 2
        i2 = 1
        x2 = x2start
        while x2 !== nothing
            ch2, state2 = x2
            #  update
            above, current, left = current, left, v0[i2]
            if ch1 != ch2
                current = min(current + 1, above + 1, left + 1)
            end
            min_dist = min(min_dist, left)
            v0[i2] = current
            x2 = iterate(s2, state2)
            i2 += 1
        end
        max_dist !== nothing && min_dist > max_dist && return max_dist + 1
        x1 = iterate(s1, state1)
        i1 += 1
    end
    max_dist !== nothing && current > max_dist && return max_dist + 1
    return current
end

##############################################################################
##
## Damerau Levenshtein
##
##############################################################################
"""
    DamerauLevenshtein()

Creates the DamerauLevenshtein metric

The DamerauLevenshtein distance is the minimum number of operations (consisting of insertions, deletions or substitutions of a single character, or transposition of two adjacent characters) required to change one string into the other.
"""
struct DamerauLevenshtein <: SemiMetric end

## http://blog.softwx.net/2015/01/optimizing-damerau-levenshtein_15.html
function evaluate(dist::DamerauLevenshtein, s1::AbstractString, s2::AbstractString;
    max_dist = nothing)
    s1, s2 = reorder(s1, s2)
    len1, len2 = length(s1), length(s2)
    max_dist !== nothing && len2 - len1 > max_dist && return max_dist + 1
    # prefix common to both strings can be ignored
    k, x1, x2start = remove_prefix(s1, s2)
    (x1 == nothing) && return len2 - k
    v0 = collect(1:(len2 - k))
    v2 = similar(v0)
    if max_dist !== nothing
        offset = 1 + max_dist - (len2 - len1)
        i2_start = 1
        i2_end = max_dist
    end
    i1 = 1
    current = i1
    prevch1, = x1
    while x1 !== nothing
        ch1, state1 = x1
        left = (i1 - 1)
        current = i1
        nextTransCost = 0
        prevch2, = x2start
        if max_dist !== nothing
            i2_start += (i1 > offset) ? 1 : 0
            i2_end = min(i2_end + 1, len2)
        end
        x2 = x2start
        i2 = 1
        while x2 !== nothing
            ch2, state2 = x2
            if max_dist == nothing || (i2_start <= i2 <= i2_end)
                above = current
                thisTransCost = nextTransCost
                nextTransCost = v2[i2]
                # cost of diagonal (substitution)
                v2[i2] = current = left
                # left now equals current cost (which will be diagonal at next iteration)
                left = v0[i2]
                if ch1 != ch2
                    # insertion
                    if left < current
                        current = left
                    end
                    # deletion
                    if above < current
                        current = above
                    end
                    current += 1
                    if (i1 != 1) & (i2 != 1) & (ch1 == prevch2) & (prevch1 == ch2)
                        thisTransCost += 1
                        if thisTransCost < current
                            current = thisTransCost
                        end
                    end
                end
                v0[i2] = current
            end
            x2 = iterate(s2, state2)
            i2 += 1
            prevch2 = ch2
        end
        max_dist !== nothing && v0[i1 + len2 - len1] > max_dist && return max_dist + 1
        x1 = iterate(s1, state1)
        i1 += 1
        prevch1 = ch1
    end
    max_dist !== nothing && current > max_dist && return max_dist + 1
    return current
end

##############################################################################
##
## Ratcliff/Obershelp
##
##############################################################################
"""
    RatcliffObershelp()

Creates the RatcliffObershelp metric

The distance between two strings is defined as one minus  the number of matching characters divided by the total number of characters in the two strings. Matching characters are those in the longest common subsequence plus, recursively, matching characters in the unmatched region on either side of the longest common subsequence.
"""
struct RatcliffObershelp <: PreMetric end

function evaluate(dist::RatcliffObershelp, s1::AbstractString, s2::AbstractString)
    n_matched = sum(last.(matching_blocks(s1, s2)))
    len1, len2 = length(s1), length(s2)
    len1 + len2 == 0 ? 0 : 1.0 - 2 *  n_matched / (len1 + len2)
end

function matching_blocks(s1::AbstractString, s2::AbstractString)
    matching_blocks!(Set{Tuple{Int, Int, Int}}(), s1, s2, length(s1), length(s2), 1, 1)
end

function matching_blocks!(x::Set{Tuple{Int, Int, Int}}, s1::AbstractString, s2::AbstractString, len1::Integer, len2::Integer, start1::Integer, start2::Integer)
    a = longest_common_substring(s1, s2, len1 , len2)
    # exit if there is no common substring
    a[3] == 0 && return x
    # add the info of the common to the existing set
    push!(x, (a[1] + start1 - 1, a[2] + start2 - 1, a[3]))
    # add the longest common substring that happens before
    s1before = SubString(s1, 1, nextind(s1, 0, a[1] - 1))
    s2before = SubString(s2, 1, nextind(s2, 0, a[2] - 1))
    matching_blocks!(x, s1before, s2before, a[1] - 1, a[2] - 1, start1, start2)
    # add the longest common substring that happens after
    s1after = SubString(s1, nextind(s1, 0, a[1] + a[3]), lastindex(s1))
    s2after = SubString(s2, nextind(s2, 0, a[2] + a[3]), lastindex(s2))
    matching_blocks!(x, s1after, s2after, len1 - (a[1] + a[3]) + 1, len2 - (a[2] + a[3]) + 1, start1 + a[1] + a[3] - 1, start2 + a[2] + a[3] - 1)
    return x
end


text:
using Bio.Seq
using BioBenchmarks

function searchall(seq, pat)
    query = ExactSearchQuery(pat)
    count = 0
    from = 0
    while (from = searchindex(seq, query, from + 1)) > 0
        count += 1
    end
    return count
end

chrom = first(open(datafile("chr1.fa"), FASTA)).seq
@benchmark searchall(chrom, dna"ATTACGTGACGT")


text:
using Revise
using JuMP, EAGO

m = Model(optimizer_with_attributes(EAGO.Optimizer, "verbosity" => 1,
                                                    "output_iterations" => 1000,
                                                    "iteration_limit" => 100000,

                                                    "cp_depth" => -1,
                                                    "cp_repetitions" => -1,
                                                    "cp_forward_reverse_limit" => 2,

                                                    "cut_min_iterations" => 3,
                                                    "cut_max_iterations" => 3,
                                                    "objective_cut_on" => true,
                                                    "subgrad_tighten" => true,

                                                    "obbt_depth" => 4,
                                                    "obbt_repetitions" => 4,
                                                    "obbt_aggressive_on" => true,

                                                    "upper_bounding_depth" => 4,

                                                    "fbbt_lp_depth" => 100000,
                                                    "fbbt_lp_repetitions" => 3))

# cp_repetitions
# obbt_repetitions

# ----- Variables ----- #
x_Idx = Any[1, 2]
@variable(m, x[x_Idx])
JuMP.set_lower_bound(x[1], -2.0)
JuMP.set_upper_bound(x[1], 2.0)
JuMP.set_lower_bound(x[2], -2.0)
JuMP.set_upper_bound(x[2], 2.0)

# FALSE FATHOM.... sqr, mult, +, -

# ----- Objective ----- #
#@NLobjective(m, Min, 3*x[1]^2 + 6*x[1]*x[2] + 3*x[2]^2 - 14*x[1] - 14*x[2])
#@NLobjective(m, Min, 19+3*(x[1])^2-14*x[1]+6*x[1]*x[2]-14*x[2]+3* (x[2])^2)
@NLobjective(m, Min, (1+ (1+x[1]+x[2])^2*(19+3* (x[1])^2-14*x[1]+6*x[1]*x[2]-14*x[2]+3* (x[2])^2))*
(30+ (2*x[1]-3*x[2])^2*(18+12* (x[1])^2-32*x[1]-36*x[1]*x[2]+48*x[2]+27* (x[2])^2)))

JuMP.optimize!(m)
run_time = backend(m).optimizer.model.optimizer._run_time
println("run time: $run_time")


text:
module Schedulers

include("../../docstrings.jl")

import ..Rule
import ..EGraph

include("./abstractscheduler.jl")
include("./backoffscheduler.jl")
include("./simplescheduler.jl")

export AbstractScheduler
export BackoffScheduler
export SimpleScheduler
export cansaturate
export shouldskip
export readstep!
export writestep!

end


text:
"""
    cumulants_to_raw_moments(N::Int, max_order::Int)

Express all `N`-variate cumulants ``κ`` up to order given by `max_order`
in terms of raw moments ``μ``. Return a Dictionary mapping from
vector ``\\mathbf{i}`` (that indicates the cumulant ``κ_{\\mathbf{i}}``)
to the corresponding raw moment expressions.

# Example

```julia
cumulants_to_raw_moments(2, 2)

Dict{Any,Any} with 5 entries:
  (1, 0) => μ₁₀(t)
  (2, 0) => μ₂₀(t) - ((μ₁₀(t))^2)
  (0, 1) => μ₀₁(t)
  (0, 2) => μ₀₂(t) - ((μ₀₁(t))^2)
  (1, 1) => μ₁₁(t) - ((μ₀₁(t))*(μ₁₀(t)))
```
"""
function cumulants_to_raw_moments(N::Int, max_order::Int, μ=nothing)

    # following Smith (1995)

    iter_all = construct_iter_all(N, max_order)
    iter_1 = filter(x -> sum(x) == 1, iter_all)

    if μ == nothing
        μ = define_μ(N, max_order, iter_all)
    else
        @assert (μ isa Dict) && length(μ) == length(iter_all) "passed arguments are inconsistent (μ vs N & order)"
    end

    K = Dict()
    μ_star = Dict()

    μ_star[Tuple(fill(0, N))] = 1.0
    for i in 1:N
        eᵢ = iter_1[i]
        K[eᵢ] = μ[eᵢ]
        μ_star[eᵢ] = -μ[eᵢ]
    end

    for order in 2:max_order

        iter_order = filter(x -> sum(x) == order, iter_all)
        for r in iter_order

            ind = findall(x -> x!= 0, r)[end]
            r_sub = r .- iter_1[ind]
            iter_i = filter(x -> all(x .<= r_sub), iter_all)

            suma = 0.0
            for i in iter_i
                factor = 1.0
                for j in 1:N
                    factor *= binomial(r_sub[j], i[j])
                end
                suma += factor*μ[r.-i]*μ_star[i]
            end
            K[r] = simplify(suma)

            suma = 0.0
            for i in iter_i
                factor = 1.0
                for j in 1:N
                    factor *= binomial(r_sub[j], i[j])
                end
                suma += factor*(-K[r.-i])*μ_star[i]
            end
            μ_star[r] = simplify(suma)
        end

    end

    K

end

"""
    cumulants_to_central_moments(N::Int, max_order::Int)

Express all `N`-variate cumulants ``κ`` up to order given by `max_order`
in terms of raw moments ``M``. Return a Dictionary mapping from
vector ``\\mathbf{i}`` (that indicates the cumulant ``κ_{\\mathbf{i}}``)
to the corresponding central moment expressions.
"""
function cumulants_to_central_moments(N::Int, max_order::Int)

    # obtain cumulants up to (m_order)^th order in terms of
    # central moments using formula from Balakrishan et al. (1998)

    K = Dict()
    M_star = Dict()

    iter_all = construct_iter_all(N, max_order)
    iter_1 = filter(x -> sum(x) == 1, iter_all)
    μ = define_μ(N, 1, iter_1)
    M = define_M(N, max_order, iter_all)

    M_star[Tuple(zeros(N))] = 1.0
    for i in 1:N
        eᵢ = iter_1[i]
        K[eᵢ] = μ[eᵢ]
        M_star[eᵢ] = 0.0
    end

    for order in 2:max_order
        # going up through the orders iteratively to build up the expressions (could do recursively?)

        iter_order = filter(x -> sum(x) == order, iter_all)
        for r in iter_order

            ind = findall(x -> x!= 0, r)[end]
            r_sub = r .- iter_1[ind]
            iter_i = filter(x -> all(x .<= r_sub), iter_all)
            # find the cumulant \kappa_{\bm{r}}}
            suma = 0.0
            for i in iter_i
                factor = 1.0
                for j in 1:N
                    factor *= binomial(r_sub[j], i[j])
                end
                suma += factor*M[r.-i]*M_star[i]
            end
            K[r] = simplify(suma)

            # Find the central moment \M^*_{\bm{r}}
            suma = 0.0
            for i in iter_i
                factor = 1.0
                for j in 1:N
                    factor *= binomial(r_sub[j], i[j])
                end
                suma += factor*(-K[r.-i])*M_star[i]
            end
            suma -= -μ[iter_1[ind]]*M_star[r_sub]
            M_star[r] = simplify(suma)

        end

    end

    K

end

function raw_to_central_moments(N::Int, order::Int, μ=nothing; bernoulli=false)

    # Return a dictionary of central moments expressed in terms of raw moments
    # example use:
    # 1 raw_to_central = raw_to_central_moments(2, 3)
    # 2 M₁₂ = raw_to_central[(1,2)] = 2 μ₁₀ μ₀₁² - 2μ₁₁μ₀₁ - μ₀₂μ₁₀ + μ₁₂
    # note that μ is an optional argument which can be used to pass
    # arbitrary values/symbols for each raw moment (different from default μᵢ)

    iter_all = construct_iter_all(N, order)
    iter_μ = filter(x -> sum(x) == 1, iter_all)
    M = define_M(N, order, iter_all)
    if μ == nothing
        μ = define_μ(N, order, iter_all)
    elseif !(μ isa Dict) || length(μ) != length(iter_all)
        if bernoulli
            iter_all = keys(μ)
        else
            error("passed arguments are inconsistent (μ vs N & order)")
        end
    else
        iter_all = construct_iter_all(N, order)
    end
    raw_to_central = Dict()

    for i in iter_all

        iter_j = Iterators.filter(x -> all(x .<= i), iter_all)
        suma = 0.0
        for j in iter_j
            term = μ[i.-j]
            for (k, e_k) in zip(1:N, iter_μ)
                term *= (-1)^(j[k])*binomial(i[k], j[k])*μ[e_k]^j[k]
            end
            suma += term
        end
        raw_to_central[i] = simplify(suma)
    end

    raw_to_central

end

function central_to_raw_moments(N::Int, order::Int)

    # Return a dictionary of raw moments expressed in terms of central moments
    # example use:
    # 1 central_to_raw = central_to_raw_moments(2, 3)
    # 2 μ₁₂ = central_to_raw[(1,2)] = 2 M₁₁ μ₀₁ + M₀₂μ₁₀ + M₁₂ + μ₁₀ μ₀₁²

    iter_all = construct_iter_all(N, order)
    iter_μ = filter(x -> sum(x) == 1, iter_all)

    M = define_M(N, order, iter_all)
    μ = define_μ(N, order, iter_μ)

    central_to_raw = Dict()

    central_to_raw[iter_all[1]] = 1.0
    for i in iter_all[2:end]

        iter_j = Iterators.filter(x -> all(x .<= i), iter_all)
        suma = 0.0
        for j in iter_j
            term = M[i.-j]
            for (k, e_k) in zip(1:N, iter_μ)
                term *= binomial(i[k], j[k])*μ[e_k]^j[k]
            end
            suma += term
        end
        central_to_raw[i] = simplify(suma)
    end

    central_to_raw

end


text:
## Convex Polyhedra
## ================
convert( ::Type{Ptr{Ptr{UInt8}}}, s::Array{ASCIIString,1} ) = map(pointer ,s)

pplisinitailsed = false

function isppl_initialzied()
  global pplisinitailsed
  pplisinitailsed
end

function ppl_initialize()
  @assert !isppl_initialzied()
  ccall( (:ppl_initialize, "libppl_c"), Int32, ())
  global pplisinitailsed = true
end

function ppl_finalize()
  @assert isppl_initialzied()
  ccall( (:ppl_finalize, "libppl_c"), Int32, ())
  global pplisinitailsed = false
end

function ppl_version_vec()
  @assert isppl_initialzied()
  Int[ccall((:ppl_version_major, "libppl_c"), Int32, ()),
      ccall((:ppl_version_minor, "libppl_c"), Int32, ()),
      ccall((:ppl_version_revision, "libppl_c"), Int32, ()),
      ccall((:ppl_version_beta, "libppl_c"), Int32, ())]
end

function ppl_version()
  @assert isppl_initialzied()
#   versionnum = Array(UInt8, 20)
  versionnum = convert(Ptr{Ptr{Uint8}},["H", "E","L","L","O"])
  res = ccall((:ppl_version, "libppl_c"), Int, (Ptr{Ptr{Uint8}},), versionnum)
  bytestring(convert(Ptr{Uint8}, versionnum[1]))
end


text:
# Input data
data = map(
    x -> (toks = split(x, " ");
    (name, v) = (toks[1], parse(Int64, toks[2]));
    ((name == "up") * v, (name == "down") * v, (name == "forward") * v)),
    readlines("Day02.input"),
);

# Part 1
ans = prod(reduce(.+, map(((u, d, f),) -> (d - u, f), data)));
println("Part 1: $ans")

# Part 2
ans = prod(
    reduce(
        .+,
        map(
            (aim, (_, _, f)) -> (f * aim, f),
            cumsum(map(((u, d, _),) -> d - u, data)),
            data,
        ),
    ),
);
println("Part 2: $ans")


text:
using BenchmarkTools, ARCH2018_RE

SUITE = BenchmarkGroup()

include("../src/SpaceStation/iss.jl")
include("../src/Building/build.jl")

#=
# From BenchmarkTools.jl
# ----------------------
paramspath = joinpath(dirname(@__FILE__), "params.json")

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