text:
import Base: display

Tp, Tm, Ts, Td = typeof(+), typeof(*), typeof(-), typeof(÷)

struct CageConstraint{T} <: Constraint where T <: Union{Tp, Tm, Ts, Td}
    total::Int
    coords::Array{Tuple{Int, Int}}
    op::T
end

"""
Creates a KenKen board given constraints. It inserts values for any `CellConstraints`
and discards them. It also creates row an column constraints, so those don't need
to be part of the input. As per the traditional formulation of constraint
satisfaction problems, each instance contains `domains` (valid remaining values for
each square), `variables` (values that have been completely determined), and
`constraints`.
"""
function KenKen(constraints, k::Int)
    variables = Array{Int, 2}(undef, k, k)
    domains = Array{BitArray, 2}(undef, k, k)
    for i = 1:k
        for j = 1:k
            domains[i, j] = BitArray{1}(undef, k) .+ 1
            variables[i, j] = -1
        end
    end
    clean_constraints = []
    for constr ∈ constraints
        # if it's a cell constraint, implement and discard. otherwise, keep
        if constr isa CellConstraint
            i, j = constr.coord
            v = constr.value
            set_cell!(i, j, v, variables)
            set_domain!(i, j, v, domains)
        else
            push!(clean_constraints, constr)
        end
    end
    row_constraints = [RowConstraint(i) for i = 1:k]
    col_constraints = [ColumnConstraint(i) for i = 1:k]
    constrs = vcat(clean_constraints, row_constraints, col_constraints)
    return GridCSP(variables, constrs, domains)
end

inverse(_::Tp) = -
inverse(_::Tm) = ÷

id(_::Tp) = 0
id(_::Tm) = 1

function _constraint_value_and_undetermined(
    constr::CageConstraint{T},
    variables::Array{Int, 2}
) where T <: Union{Tp, Tm}
    acc = id(constr.op)
    undetermined = []
    for (i, j) ∈ constr.coords
        v = variables[i, j]
        if v > 0
            acc = constr.op(acc, v)
        else
            push!(undetermined, (i, j))
        end
    end
    return acc, undetermined
end

"""
Eliminates candidate values for an addition cage constraint.
"""
function eliminate!(constr::CageConstraint{Tp}, csp::GridCSP)
    acc, undetermined = _constraint_value_and_undetermined(constr, csp.variables)
    if length(undetermined) == 0
        return
    end
    k = size(csp)
    # if there is one cell left, anything not equal to the total - acc is invalid.
    # if there more than one undetermined, any single cell ≥ to the total - acc
    # is invalid
    comp = length(undetermined) == 1 ? (≠) : (≥)
    for (i, j) ∈ undetermined
        for cand ∈ 1:k
            if comp(cand, constr.total - acc)
                csp.domains[i, j][cand] = 0
            end
        end
    end
end

"""
Eliminates candidate values for a multiplication cage constraint.
"""
function eliminate!(constr::CageConstraint{Tm}, csp::GridCSP)
    acc, undetermined = _constraint_value_and_undetermined(constr, csp.variables)
    if length(undetermined) == 0
        return
    end
    k = size(csp)
    for (i, j) ∈ undetermined
        for cand ∈ 1:k
            # to remain valid, candidates must divide the remaining amount
            if (constr.total ÷ acc) % cand != 0
                csp.domains[i, j][cand] = 0
            end
        end
    end
end

"""
Eliminates candidate values for a subtraction cage constraint.
"""
function eliminate!(constr::CageConstraint{Ts}, csp::GridCSP)
    # We can assume there are only two cells in this cage.
    (i1, j1), (i2, j2) = constr.coords
    n1, n2 = csp.variables[i1, j1], csp.variables[i2, j2]
    if (n1 > 0 && n2 > 0) || (n1 < 0 && n2 < 0)
        return
    end
    k = size(csp)
    i, j = n1 > 0 ? (i2, j2) : (i1, j1)
    n = n1 > 0 ? n1 : n2
    csp.domains[i, j] .= 0
    # because the difference between n and the undetermined value must be
    # constr.total, once n has been set, the only other valid values are
    # n + constr.total (provided that it's ≤ k) and n - constr.total
    # (provided that it's ≥ 1)
    if n + constr.total ≤ k
        csp.domains[i, j][n + constr.total] = 1
    end
    if n - constr.total ≥ 1
        csp.domains[i, j][n - constr.total] = 1
    end
end

"""
Eliminates candidate values for a division cage constraint.
"""
function eliminate!(constr::CageConstraint{Td}, csp::GridCSP)
    # We can assume that there are only two cells in this cage.
    (i1, j1), (i2, j2) = constr.coords
    n1, n2 = csp.variables[i1, j1], csp.variables[i2, j2]
    if (n1 > 0 && n2 > 0) || (n1 < 0 && n2 < 0)
        return
    end
    # We can now assume that exactly one value has been set
    k = size(csp)
    i, j = n1 > 0 ? (i2, j2) : (i1, j1)
    n = n1 > 0 ? n1 : n2
    csp.domains[i, j] .= 0
    # Similar logic as the subtraction case
    if n * constr.total ≤ k
        csp.domains[i, j][n * constr.total] = 1
    end
    if n ÷ constr.total ≥ 1
        csp.domains[i, j][n ÷ constr.total] = 1
    end
end

function neighbors(coord::Tuple{Int, Int}, k::Int)
    i, j = coord
    neighbors = Set()
    for d = [1, -1]
        if 1 ≤ i + d ≤ k
            push!(neighbors, (i + d, j))
        end
        if 1 ≤ j + d ≤ k
            push!(neighbors, (i, j + d))
        end
    end
    return neighbors
end

function generate_cage!(
    cell::Tuple{Int, Int},
    used::Set{Tuple{Int, Int}},
    range::UnitRange,
    k::Int
)
    cage = Set{Tuple{Int, Int}}()
    generate_cage!(cell, cage, used, rand(range), k)
    return cage
end

function generate_cage!(
    cell::Tuple{Int, Int},
    acc::Set{Tuple{Int, Int}},
    used::Set{Tuple{Int, Int}},
    sz::Int,
    k::Int
)
    # if the cage has reached the right size
    if length(acc) == sz
        return
    end
    push!(acc, cell)
    push!(used, cell)
    nbrs = setdiff(neighbors(cell, k), used)
    # with some probability, recurse
    denom = 1.
    for (i, j) ∈ nbrs
        if rand() < 1 / denom
            generate_cage!((i, j), acc, used, sz, k)
            denom += 1
        end
    end
    return acc
end

function next_unused_cell_with_max_neighbors(used::Set{Tuple{Int, Int}}, k::Int)
    @assert length(used) < k^2
    max_neighbors = 0
    imax, jmax = nothing, nothing
    for i = 1:k
        for j = 1:k
            if (i, j) ∉ used
                n_neighbors = length(neighbors((i, j), k))
                if isnothing(imax) && isnothing(jmax) || n_neighbors ≥ max_neighbors
                    imax, jmax = i, j
                end
                max_neighbors = max(max_neighbors, n_neighbors)
            end
        end
    end
    return imax, jmax
end

function cage_value(op::Union{Tp, Tm}, cage::Set{Tuple{Int, Int}}, board::Array{Int, 2})
    # the order of operands doesnt matter for multiplication and addition. we can
    # also have more than 2 cells
    t = id(op)
    for (i, j) ∈ cage
        t = op(t, board[i, j])
    end
    return t
end

function cage_value(op::Union{Td, Ts}, cage::Set{Tuple{Int, Int}}, board::Array{Int, 2})
    (i1, j1), (i2, j2) = cage
    n1, n2 = board[i1, j1], board[i2, j2]
    M, m = max(n1, n2), min(n1, n2)
    return op(M, m)
end

function get_operation_options(cage::Set{Tuple{Int, Int}}, board::Array{Int, 2})
    options = []
    # in order to have ÷ and - as options, the cage must be of size 2
    if length(cage) == 2
        (i1, j1), (i2, j2) = cage
        n1, n2 = board[i1, j1], board[i2, j2]
        M, m = max(n1, n2), min(n1, n2)
        # in order to have division, one of the cell values must divide the other
        if M % m == 0
            push!(options, ÷)
        end
        push!(options, -)
    end
    return vcat(options, [+, *])
end

"""
Generate a random KenKen of a certain size. The solution is guaranteed to be
unique.
"""
function generate_random_kenken(k::Int)
    row_constrs = [RowConstraint(i) for i = 1:k]
    col_constrs = [ColumnConstraint(i) for i = 1:k]
    constrs = vcat(row_constrs, col_constrs)
    kk, solution = nothing, nothing
    solution_is_unique = false
    # we need to keep trying until the solution is unique
    while !solution_is_unique
        # create a grid csp that will provide a grid whose rows and columns
        # respectively have no duplicates
        g = backtrack(GridCSP(constrs, k), 1; randomize = true)
        solution = g.variables
        used = Set{Tuple{Int, Int}}()
        kk_constrs = Constraint[]
        # tile the grid with cages
        while length(used) < k^2
            cell = next_unused_cell_with_max_neighbors(used, k)
            cage = generate_cage!(cell, used, 2:4, k)
            # if the cage has one cell, must be a cell constraint
            if length(cage) == 1
                i, j = collect(cage)[1]
                next_constr = CellConstraint(g.variables[i, j], (i, j))
            else
                # otherwise, choose a valid operation at random collect the cell
                # value and add a cage constraint
                options = get_operation_options(cage, g.variables)
                op = rand(options)
                val = cage_value(op, cage, g.variables)
                next_constr = CageConstraint(val, collect(cage), op)
            end
            push!(kk_constrs, next_constr)
        end
        # solve the resulting problem and if the solution is unique, we're done
        # otherwise, try again
        kk = KenKen(kk_constrs, k)
        solutions = backtrack(kk, 2)
        solution_is_unique = length(solutions) == 1
    end
    return kk, solution
end

strop(_::Tp) = '+'
strop(_::Tm) = 'x'
strop(_::Ts) = '-'
strop(_::Td) = '÷'

"""
Prints out a kenken board. Example output:
a=5+ b=9x c=2x d=3
3×3 Array{Char,2}:
 'a'  'a'  'd'
 'b'  'a'  'c'
 'b'  'b'  'c'
"""
function print_kenken(kk::GridCSP)
    k = size(kk)
    legend = []
    grid = -ones(Int, k, k)
    c = 'a'
    # cage constraints
    for constr ∈ kk.constraints
        if constr isa CageConstraint
            for (i, j) ∈ constr.coords
                grid[i, j] = c
            end
            push!(legend, "$c=$(constr.total)$(strop(constr.op))")
            c += 1
        end
    end
    # cell constraints won't have been addressed yet
    for i = 1:k
        for j = 1:k
            if grid[i, j] < 0
                grid[i, j] = c
                push!(legend, "$c=$(kk.variables[i, j])")
                c += 1
            end
        end
    end
    println(join(legend, " "))
    display(convert(Array{Char, 2}, grid))
    println()
end

max_stars_repo_name: gindij/CSPSolver

max_issues_repo_name: gindij/CSPSolver

max_forks_repo_name: gindij/CSPSolver


text:
TestParam(
  skip_tex_files = ["preamble.tex"]
)

ext: jl

max_stars_repo_name: jgoldfar/paper-build-demo

max_issues_repo_name: jgoldfar/paper-build-demo

max_forks_repo_name: jgoldfar/paper-build-demo


text:
struct Vertice
    valor_alocado::Float64
    valor_mam::Float64
    function Vertice(valor_alocado::Float64, valor_mam::Float64)
        @assert valor_alocado > -eps() "valor alocado deve ser maior ou igual a zero"
        @assert valor_mam > -eps()  "valor mam deve ser maior ou igual a zero"
        new(valor_alocado, valor_mam)
    end
end

function Base.:+(x::Vertice, y::Vertice)
    return Vertice(x.valor_alocado + y.valor_alocado, x.valor_mam + y.valor_mam)
end

struct Fluxos
    vertices::Dict{Symbol, Vertice}
    Fluxos() = new(Dict{Symbol, Vertice}())
end

function Base.getindex(X::Fluxos, i::Symbol)::Vertice
    return X.vertices[i]
end
Base.getindex(X::Fluxos, i::Int64)::Vertice = X[get_codigo_vertice(i)]

function Base.setindex!(X::Fluxos, v::Vertice, i::Symbol)
    if i != get_codigo_vertice(12)
        @assert v.valor_mam < eps() "apenas o vertice :12 pode ter valor MaM maior que zero"
    end
    X.vertices[i] = v
end
Base.setindex!(X::Fluxos, v::Vertice, i::Int64) = X[get_codigo_vertice(i)] = v

function add_vertice!(fluxos::Fluxos, codigo_vertice::Symbol, vertice::Vertice)
    if haskey(fluxos.vertices, codigo_vertice)
        fluxos[codigo_vertice] += vertice
    else
        fluxos[codigo_vertice] = vertice
    end
end
add_vertice!(fluxos::Fluxos, i::Int64, vertice::Vertice) = add_vertice!(fluxos, get_codigo_vertice(i), vertice)

sorted_keys(fluxos::Fluxos) = sort(collect(keys(fluxos.vertices)))

function get_codigo_vertice(v::Int64)::Symbol
    @assert v >= 1
    @assert v <= 12
    if v < 10
        # add leading zero
        return Symbol("0$v")
    else
        return Symbol("$v")
    end
end

SymbolOrNothing = Union{Symbol, Nothing}

struct ItemCarteira
    item::Symbol
    id_posicao::SymbolOrNothing
    fator_risco::Symbol
    local_registro::SymbolOrNothing
    carteira_negoc::Symbol
    is_abertura_fundo::Bool
    function ItemCarteira(
        item::Symbol,
        id_posicao::SymbolOrNothing,
        fator_risco::Symbol,
        local_registro::SymbolOrNothing,
        carteira_negoc::Symbol,
        is_abertura_fundo::Bool=false)

        @assert item in CONTAS "item invalido: $item"
        id_posicao != nothing && @assert id_posicao in POSICOES "posicao invalida: $id_posicao"
        @assert fator_risco in FATORES_RISCO "fator de risco invalido: $fator_risco"
        local_registro != nothing && @assert local_registro in LOCAIS_REGISTRO "local_registro invalido: $local_registro"
        @assert carteira_negoc in CARTEIRAS "carteira invalida: $carteira_negoc"
        is_abertura_fundo && @assert item in CONTAS "item fundo invalido: $item"

        return new(item, id_posicao, fator_risco, local_registro, carteira_negoc, is_abertura_fundo)
    end
end

is_ativo(item::ItemCarteira)::Bool                = !item.is_abertura_fundo && (item.item in CONTAS_ATIVO)
is_ativo_fundo(item::ItemCarteira)::Bool          = item.is_abertura_fundo && (item.item in CONTAS_ATIVO)

is_passivo(item::ItemCarteira)::Bool              = !item.is_abertura_fundo && (item.item in CONTAS_PASSIVO)
is_passivo_fundo(item::ItemCarteira)::Bool        = item.is_abertura_fundo && (item.item in CONTAS_PASSIVO)

is_derivativo(item::ItemCarteira)::Bool           = !item.is_abertura_fundo && (item.item in CONTAS_DERIVATIVO)
is_derivativo_fundo(item::ItemCarteira)::Bool     = item.is_abertura_fundo && (item.item in CONTAS_DERIVATIVO)

is_atividade_financeira(item::ItemCarteira)::Bool = item.item in CONTAS_ATIVIDADE_FINANCEIRA

struct Documento
    id_docto::String
    id_docto_versao::String
    data_base::String
    id_inst_financ::Int64
    tipo_arq::Symbol
    nome_contato::String
    fone_contato::String
    ativo::Dict{ItemCarteira, Fluxos}                   # internal
    passivo::Dict{ItemCarteira, Fluxos}                 # internal
    derivativo::Dict{ItemCarteira, Fluxos}              # internal
    ativo_fundo::Dict{ItemCarteira, Fluxos}             # internal
    passivo_fundo::Dict{ItemCarteira, Fluxos}           # internal
    derivativo_fundo::Dict{ItemCarteira, Fluxos}        # internal
    atividade_financeira::Dict{ItemCarteira, Fluxos}    # internal

    function Documento(
        id_docto::String,
        id_docto_versao::String,
        data_base::String,
        id_inst_financ::Int64,
        tipo_arq::Symbol,
        nome_contato::String,
        fone_contato::String
    )
        @assert id_docto in ["2060"] "id documento invalido: $id_docto"
        @assert occursin(r"^v\d+$", id_docto_versao) "versao documento invalida: $id_docto_versao. exemplo: v1"
        @assert occursin(r"^\d{4}-\d{2}", data_base) "data base invalida: $data_base. exemplo: 2020-06"
        @assert tipo_arq in TIPOS_ARQUIVO "tipo arquivo invalido: $tipo_arq"
        @assert nome_contato != "" "nome contato obrigatorio"
        @assert fone_contato != "" "telefone contato obrigatorio"
        return new(
            id_docto, id_docto_versao, data_base, id_inst_financ,
            tipo_arq, nome_contato, fone_contato,
            # blank dicts
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}(),
            Dict{ItemCarteira, Fluxos}()
        )
    end
end

"""
Adds a vertice data into a document.
Document section is infered. Please check `_infer_document_section()`
"""
function Base.push!(doc::Documento, item::ItemCarteira, codigo_vertice::Union{Symbol, Int64}, vertice::Vertice)
    section = _infer_document_section(doc, item)
    if !haskey(section, item)
        section[item] = Fluxos() # init with blank
    end
    add_vertice!(section[item], codigo_vertice, vertice)
end

"Returns the `Document` section where `ItemCarteira` should be included. Depends on `ItemCarteira.item` code."
function _infer_document_section(doc::Documento, item_carteira::ItemCarteira)
    if is_ativo(item_carteira)
        return doc.ativo
    elseif is_ativo_fundo(item_carteira)
        return doc.ativo_fundo
    elseif is_passivo(item_carteira)
        return doc.passivo
    elseif is_passivo_fundo(item_carteira)
        return doc.passivo_fundo
    elseif is_derivativo(item_carteira)
        return doc.derivativo
    elseif is_derivativo_fundo(item_carteira)
        return doc.derivativo_fundo
    elseif is_atividade_financeira(item_carteira)
        return doc.atividade_financeira
    else
        error("unicorn!")
    end
end


text:
function instantiate_problem(dev, stepper)
    problem = Problem(nx=4, dev=dev, stepper=stepper)
    return typeof(problem) <: FourierFlows.Problem
end

function instantiate_problem_with_filter_kwargs(dev, stepper)
    dummy_problem = Problem(nx=16, dev=dev, stepper=stepper)

    real_problem = FourierFlows.Problem(dummy_problem.eqn,
                                        stepper,
                                        1.0,
                                        dummy_problem.grid,
                                        dummy_problem.vars,
                                        dummy_problem.params,
                                        dev,
                                        innerK=0.0,
                                        outerK=1/16)

    stepper = real_problem.timestepper

    return stepper.filter[3] < 1e-16
end


text:
"""
`init_subspec!(m::RealBondMkup)`

Initializes a model subspecification by overwriting parameters from
the original model object with new parameter objects. This function is
called from within the model constructor.
"""
function init_subspec!(m::RealBondMkup)
    if subspec(m) == "ss0"
        return
    elseif subspec(m) == "ss1"
        return ss1!(m)
    elseif subspec(m) == "ss2"
        return ss2!(m)
    elseif subspec(m) == "ss3"
        return ss3!(m)
    elseif subspec(m) == "ss4"
        return ss4!(m)
    elseif subspec(m) == "ss5"
        return ss5!(m)
    elseif subspec(m) == "ss6"
        return ss6!(m)
    else
        error("This subspec should be a 0")
    end
end

"""
```
ss1!(m::RealBondMkup)
```

Initializes subspec 1 of `RealBondMkup`.
This shuts down all shocks except the technology process (z shock).
"""
function ss1!(m::RealBondMkup)
    m <= parameter(:ρ_z, 0.0, fixed = true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mon, 0.0, fixed = true,
                   description="ρ_mon: Persistence of the monetary policy shock",
                   tex_label="\\rho_{mon}")
    m <= parameter(:ρ_mkp, 0.0, fixed = true,
                   description="ρ_mkp: Persistence of the markup shock",
                   tex_label="\\rho_{mkp}")

    m <= parameter(:σ_mon, 0.0, fixed = true,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
    m <= parameter(:σ_mkp, 0.0, fixed = true,
                   description="σ_mkp: The standard deviation of the markup shock.",
                   tex_label="\\sigma_{mkp}")
end

"""
```
ss2!(m::RealBondMkup)
```

Initializes subspec 2 of `RealBondMkup`.
This shuts down all shocks except the markup shock process (mkp shock).
"""
function ss2!(m::RealBondMkup)
    m <= parameter(:ρ_z, 0.0, fixed = true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mon, 0.0, fixed = true,
                   description="ρ_mon: Persistence of the monetary policy shock",
                   tex_label="\\rho_{mon}")
    m <= parameter(:ρ_mkp, 0.0, fixed = true,
                   description="ρ_mkp: Persistence of the markup shock",
                   tex_label="\\rho_{mkp}")

    m <= parameter(:σ_z, 0.0, fixed = true,
                   description="σ_z: The standard deviation of the technology process.",
                   tex_label="\\sigma_z")
    m <= parameter(:σ_mon, 0.0, fixed = true,
                   description="σ_mon: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_{mon}")
end

"""
```
ss3!(m::RealBondMkup)
```

Initializes subspec 3 of `RealBondMkup`.
This shuts down all shocks except the monetary policy process (mon shock).
"""
function ss3!(m::RealBondMkup)
    m <= parameter(:ρ_z, 0.0, fixed = true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mon, 0.0, fixed = true,
                   description="ρ_mon: Persistence of the monetary policy shock",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mkp, 0.0, fixed = true,
                   description="ρ_mkp: Persistence of the markup shock",
                   tex_label="\\rho_z")

    m <= parameter(:σ_z, 0.0, fixed = true,
                   description="σ_z: The standard deviation of the monetary policy shock.",
                   tex_label="\\sigma_z")
    m <= parameter(:σ_mkp, 0.0, fixed = true,
                   description="σ_mkp: The standard deviation of the markup shock.",
                   tex_label="\\sigma_{mkp}")
end

"""
```
ss4!(m::RealBondMkup)
```

Initializes subspec 4 of `RealBondMkup`.
This fixes all parameters except for the sigmas. Used for estimating only the sigmas
"""
function ss4!(m::RealBondMkup)
    # Steady State Parameters
    m <= parameter(:R, 1.04, fixed = true,
                   description = "R: Steady-state gross real interest rate.", tex_label = "R")
    m <= parameter(:γ, 1.0, fixed = true,
                   description = "γ: CRRA Parameter.", tex_label = "\\gamma")
    m <= parameter(:ν, 1.0, fixed = true,
                   description = "Inverse Frisch elasticity of labor supply.", tex_label = "\\nu")
    m <= parameter(:abar, -0.5, fixed = true,
                   description = "Borrowing floor.", tex_label = "\\bar{a}")

    # Ar for observables (non-steadystate)
    m <= parameter(:ρ_z, 0.03081565350294113, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed=true,
                   description="ρ_z: AR(1) coefficient in the technology process.",
                   tex_label="\\rho_z")
    m <= parameter(:ρ_mon, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = true,
                   description = "ρ_mon: Persistence of monetary policy shock",
                   tex_label = "\\rho_{mon}")
    m <= parameter(:ρ_mkp, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = true,
                   tex_label = "\\rho_{mkp}", description = "ρ_mkp: Persistence of the markup shock")

    # More non-steadystate parameters
    m <= parameter(:ρ_tay, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = true,
                   description = "ρ_tay: Persistence in the taylor rule",
                   tex_label = "rho_{tay}")
    m <= parameter(:κ, 1.0, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),
                   fixed = true,
                   description = "κ: The slope of the Phillips curve",
                   tex_label = "\\kappa")
    m <= parameter(:phipi, 1.5, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.25),
                   fixed = true,
                   description = "phipi: The slope of the taylor rule",
                   tex_label = "\\phi_\\pi")

    # More steadystate parameters
    m <= parameter(:μ_s, 0., fixed = true, description = "μ_s: Mu of log normal in income")
    m <= parameter(:σ_s, 0.1, fixed = true,
                   description = "σ_s: Sigma of log normal in income")
    m <= parameter(:e_y, 1e-3, fixed = true, description = "e_y: Measurement error on GDP",
                   tex_label = "e_y")
end

"""
```
ss5!(m::RealBondMkup)
```

Initializes subspec 5 of `RealBondMkup`.
This fixes all parameters except for the sigmas and rhos. Used for estimating only the sigmas and rhos.
"""
function ss5!(m::RealBondMkup)
    # Steady State Parameters
    m <= parameter(:R, 1.04, fixed = true,
                   description = "R: Steady-state gross real interest rate.", tex_label = "R")
    m <= parameter(:γ, 1.0, fixed = true,
                   description = "γ: CRRA Parameter.", tex_label = "\\gamma")
    m <= parameter(:ν, 1.0, fixed = true,
                   description = "Inverse Frisch elasticity of labor supply.", tex_label = "\\nu")
    m <= parameter(:abar, -0.5, fixed = true,
                   description = "Borrowing floor.", tex_label = "\\bar{a}")

    # More non-steadystate parameters
    m <= parameter(:ρ_tay, 0.5, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), BetaAlt(0.5, 0.2),
                   fixed = true,
                   description = "ρ_tay: Persistence in the taylor rule",
                   tex_label = "rho_{tay}")
    m <= parameter(:κ, 1.0, (1e-3, 0.999), (1e-3, 0.999), SquareRoot(), Distributions.Uniform(0.0, 1.0),
                   fixed = true,
                   description = "κ: The slope of the Phillips curve",
                   tex_label = "\\kappa")
    m <= parameter(:phipi, 1.5, (1e-5, 10.), (1e-5, 10.00), ModelConstructors.Exponential(), Normal(1.5, 0.25),
                   fixed = true,
                   description = "phipi: The slope of the taylor rule",
                   tex_label = "\\phi_\\pi")

    # More steadystate parameters
    m <= parameter(:μ_s, 0., fixed = true, description = "μ_s: Mu of log normal in income")
    m <= parameter(:σ_s, 0.1, fixed = true,
                   description = "σ_s: Sigma of log normal in income")
    m <= parameter(:e_y, 1e-3, fixed = true, description = "e_y: Measurement error on GDP",
                   tex_label = "e_y")
end

"""
```
ss6!(m::RealBondMkup)
```

Initializes subspec 6 of `RealBondMkup`.
This fixes only the steady-state parameters. Used for estimating all non-steady state parameters. Right not, it's the same as ss0 since by default (when constructing the model), we do this.
"""
function ss6!(m::RealBondMkup)
    # Steady State Parameters
    m <= parameter(:R, 1.04, fixed = true,
                   description = "R: Steady-state gross real interest rate.", tex_label = "R")
    m <= parameter(:γ, 1.0, fixed = true,
                   description = "γ: CRRA Parameter.", tex_label = "\\gamma")
    m <= parameter(:ν, 1.0, fixed = true,
                   description = "Inverse Frisch elasticity of labor supply.", tex_label = "\\nu")
    m <= parameter(:abar, -0.5, fixed = true,
                   description = "Borrowing floor.", tex_label = "\\bar{a}")

    # More steadystate parameters
    m <= parameter(:μ_s, 0., fixed = true, description = "μ_s: Mu of log normal in income")
    m <= parameter(:σ_s, 0.1, fixed = true,
                   description = "σ_s: Sigma of log normal in income")
    m <= parameter(:e_y, 1e-3, fixed = true, description = "e_y: Measurement error on GDP",
                   tex_label = "e_y")
end


text:
struct ValueTest <: TrainingTest
    idx::Integer
    state::AbstractState
    image::AbstractImage
    value::Float32
end

mutable struct ValueTests <: TrainingTests
    tests::Vector{ValueTest}
    images::AbstractArray{Float32, 4}
    actual_value::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function ValueTests(tests; period=1, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            Flux.unsqueeze(get_value.(tests), 1),
            period,
            warmup,
        )
    end
end

function gettestfn(tests::ValueTests)
    function fn(n::Network)
        getvaluemse(
            tests,
            initialinferencefn(tests)(n)
        )
    end
end

struct ValueTestsSVG
    tests::ValueTests
end
svg(x::ValueTests) = ValueTestsSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::ValueTestsSVG)
    num_tests = length(x.tests)
    range = 1:(min(10, num_tests))
    subset = 1:num_tests
    subset_indices = subset[range]
    svgbatch = [gridsvg(x.tests[i]) for i in subset_indices]
    Base.show(io, m, GridSVG(svgbatch, main_title="Value Tests"))
end


text:
using Test
using KernelDensityEstimatePlotting
using IncrementalInference
using RoME
using RoMEPlotting

##

include("testPose2Point2Plotting.jl")
include("testPlotKDE.jl")
include("testPlotSLAMandPose.jl")

# @warn "plotMCMC needs ImageMagick on osx, not running test yet."
# plotMCMC(tree, :x1, show=false)
# println("Success")

# test ellipse function
# Z = MvNormal([0.0;0],[9.0 5; 5 4])
# eX = covEllipseParameterized(Z)

## legacy

# println("[TEST] with local Graphs.jl dictionary and arrays only (multicore)...")
# include(joinpath(dirname(@__FILE__),"..","..","IncrementalInference","test","fourdoortest.jl"))
# println("Success")
#
# println("[TEST] plot functions...")
# using Gadfly
# # draw all beliefs
# DOYTICKS = false
# xx,ll = ls(fg)
# msgPlots = drawHorBeliefsList(fg, xx, gt=gt,nhor=2);
# pl = vstack(msgPlots...);
# # Gadfly.draw(PDF("/tmp/test.pdf",15cm,30cm),pl)
# println("Success")


text:
include("Stubs.jl")
include("Stubs_test.jl")


text:
using OSQP
using Test, SparseArrays, LinearAlgebra
using Random: seed!

tests = [
    "basic.jl",
    "dual_infeasibility.jl",
    "feasibility.jl",
    "non_convex.jl",
    "polishing.jl",
    "primal_infeasibility.jl",
    "unconstrained.jl",
    "warm_start.jl",
    "update_matrices.jl",
    "MPB_wrapper.jl",
    "MOI_wrapper.jl",
    "interface.jl",
    ]

println("Running tests:")

for curtest in tests
    println(" Test: $(curtest)")
    include(curtest)
end


text:
@testset "Integration" begin
    res = ArbTools.integrate(AcbField(333), x -> sin(x + exp(x)), 0, 8)
    @test string(res) == "[0.34740017265724780787951215911989312465745625486618018388549271361674821398878532052968510434660 +/- 5.31e-96] + i*0"
end

@testset "Integration tolerance" begin
    C = AcbField(64)
    f(x) = exp(-1000 + x)*sin(10x)
    res = ArbTools.integrate(C, f, 0, 1)
    @test string(res) == "[+/- 1.38e-434] + i*0"
    res = ArbTools.integrate(C, f, 0, 1, abs_tol = 0)
    @test string(res) == "[1.574528586972758e-435 +/- 7.34e-451] + i*0"
    res = ArbTools.integrate(C, f, 0, 1, abs_tol = exp(-1000)/(2.0)^64)
    @test string(res) == "[1.574528586972758e-435 +/- 7.34e-451] + i*0"
    g(x) = exp(1000 + x)*sin(10x)
    res = ArbTools.integrate(C, g, 0, 1, rel_tol = 0)
    @test string(res) == "[+/- 6.94e+433] + i*0"
    res = ArbTools.integrate(C, g, 0, 1, rel_tol = 1e-10)
    @test string(res) == "[6.111029167e+433 +/- 2.66e+423] + i*0"
end

@testset "Integration branch cut" begin
    h(x, analytic = false) = begin
        onemx2 = 1 - x^2
        if analytic && contains_nonpositive(real(onemx2))
            return parent(x)(NaN, NaN)
        else
            return sqrt(onemx2)
        end
    end
    res = ArbTools.integrate(AcbField(64), h, 0, 1)
    @test string(res) == "[0.796113019377273535 +/- 7.17e-19] + i*[+/- 2.32e-31]"
end


text:
## --- Weighted means

    # Plain vanilla weighted mean
    @test awmean([0,1,2,3],[1,1,1,1]) == (1.5, 0.5, 5/3)
    @test awmean(1:10,ones(10)) == (5.5, 0.31622776601683794, 9.166666666666666)

    # Weighted mean with geochronogists' MSWD-correction to uncertainty
    @test gwmean([0,1,2,3],[1,1,1,1]) == (1.5, 0.6454972243679028, 5/3)
    @test gwmean(1:10,ones(10)) == (5.5, 0.9574271077563381, 9.166666666666666)

    # Mean Square of Weighted Deviates (aka reduced chi-squared)
    @test MSWD(0:10, ones(11)) == 11
    @test MSWD([0,1,2],[1,1,1]) == 1.0

## --- York fit

    x = [0.9304, 2.2969, 2.8047, 3.7933, 5.3853, 6.1995, 6.7479, 8.1856, 8.7423, 10.2588]
    y = [0.8742, 2.1626, 3.042, 3.829, 5.0116, 5.5614, 6.7675, 7.8856, 9.6414, 10.4955]
    fobj = yorkfit(x, ones(10)/4, y, ones(10)/4)
    @test fobj.intercept ≈  -0.23498964673701916
    @test fobj.intercept_sigma ≈ 0.02250863813481163
    @test fobj.slope ≈ 1.041124018512526
    @test fobj.slope_sigma ≈ 0.0035683808205783673
    @test fobj.mswd ≈ 1.1419901440278089

    @test display(fobj) != NaN

## ---

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
# A chain should behave exactly like an array of `ChainElements`. Internally, it may store
# its elements in a more efficient way.

# It can be used as a heap with `heapify!`, `heappush!`, and `heapop!`. When used this way,
# it will only return unique elements, summing elements with the same simplex together.

# `clean!` can be used to sum duplicates together and sort the chain.
"""
    Chain{F,S,E} <: AbstractVector{ChainElement{S,F}}

An internal representation of a chain. Behaves like an array of pairs `S => F`, where `S` is
the simplex type, and `F` is the coefficient type (probably a subtype of [`Mod`](@ref)).

Most functions that can be called on `AbstractCell`s can also be called on the elements of a
`Chain`.

# Examples

```jldoctest; setup=(using Random; Random.seed!(1337))
julia> data = [(rand(), rand(), rand()) for _ in 1:100];

julia> result = ripserer(data; reps=true, modulus=7);

julia> chain = result[end][end].representative
6-element Chain{Mod{7},Simplex{1, Float64, Int64}}:
 +Simplex{1}((87, 59), 0.23148225999797645) => 6 mod 7
 +Simplex{1}((59, 46), 0.3054281021426286) => 1 mod 7
 +Simplex{1}((87, 14), 0.32453294355760326) => 6 mod 7
 +Simplex{1}((87, 1), 0.34642558062390577) => 6 mod 7
 +Simplex{1}((100, 87), 0.3480194268484163) => 1 mod 7
 +Simplex{1}((79, 1), 0.36519064466525686) => 6 mod 7

julia> simplex.(chain)
6-element Vector{Simplex{1, Float64, Int64}}:
 +Simplex{1}((87, 59), 0.23148225999797645)
 +Simplex{1}((59, 46), 0.3054281021426286)
 +Simplex{1}((87, 14), 0.32453294355760326)
 +Simplex{1}((87, 1), 0.34642558062390577)
 +Simplex{1}((100, 87), 0.3480194268484163)
 +Simplex{1}((79, 1), 0.36519064466525686)

julia> vertices.(chain)
6-element Vector{Tuple{Int64, Int64}}:
 (87, 59)
 (59, 46)
 (87, 14)
 (87, 1)
 (100, 87)
 (79, 1)

julia> coefficient.(chain)
6-element Vector{Mod{7}}:
 6 mod 7
 1 mod 7
 6 mod 7
 6 mod 7
 1 mod 7
 6 mod 7

```
"""
struct Chain{F,S,E} <: AbstractVector{ChainElement{S,F}}
    elements::Vector{E}
end

_internal_eltype(::Type{Mod{2}}, ::Type{S}) where {S} = S
_internal_eltype(::Type{F}, ::Type{S}) where {F,S} = chain_element_type(S, F)
_internal_eltype(chain::Chain) = eltype(chain.elements)

function Chain{F,S}() where {F,S}
    E = _internal_eltype(F, S)
    return Chain{F,S,E}(_internal_eltype(F, S)[])
end
function Chain{F,S}(elements) where {F,S}
    E = _internal_eltype(F, S)
    return Chain{F,S,E}(convert.(E, elements))
end
function Chain(elements::AbstractVector{ChainElement{S,F}}) where {S,F}
    return Chain{F,S}(collect(elements))
end

function Base.summary(io::IO, chain::Chain{F,S}) where {F,S}
    return print(io, length(chain), "-element Chain{$F,$S}")
end

simplex_type(chain::Chain) = simplex_type(eltype(chain))
coefficient_type(chain::Chain) = coefficient_type(eltype(chain))

# Array overloads
function Base.size(chain::Chain)
    return size(chain.elements)
end
@propagate_inbounds function Base.getindex(chain::Chain{Mod{2}}, i::Int)
    return eltype(chain)(chain.elements[i])
end
@propagate_inbounds function Base.getindex(chain::Chain, i::Int)
    el = chain.elements[i]
    return eltype(chain)(simplex(el), coefficient(el))
end
@propagate_inbounds function Base.getindex(chain::Chain{F,S,E}, is) where {F,S,E}
    return Chain{F,S,E}(chain.elements[is])
end
function Base.setindex!(chain::Chain, val, i::Int)
    return chain.elements[i] = convert(_internal_eltype(chain), val)
end
function Base.setindex!(chain::Chain{F,S}, (s, c)::Tuple{S,F}, i::Int) where {F,S}
    element = ChainElement{S,F}(s, c)
    return chain[i] = element
end
function Base.IndexStyle(::Chain)
    return Base.IndexLinear()
end
function Base.resize!(chain::Chain, n)
    resize!(chain.elements, n)
    return chain
end
function Base.empty!(chain::Chain)
    empty!(chain.elements)
    return chain
end
function Base.pop!(chain::Chain{F,S}) where {F,S}
    return convert(ChainElement{S,F}, pop!(chain.elements))
end
function Base.copy(chain::Chain{F,S,E}) where {F,S,E}
    return Chain{F,S,E}(copy(chain.elements))
end
function Base.sizehint!(chain::Chain, n)
    sizehint!(chain.elements, n)
    return chain
end
# Compat: these methods are here for compatibility reasons.
# TODO: Remove when 1.6 becomes LTS
function Base.append!(chain::Chain, elements)
    n = length(chain)
    resize!(chain, n + length(elements))
    @inbounds for (i, e) in enumerate(elements)
        chain[i + n] = e
    end
    return chain
end
function Base.append!(chain1::C, chain2::C) where {C<:Chain}
    return append!(chain1.elements, chain2.elements)
end
function Base.push!(chain::Chain, element)
    resize!(chain, length(chain) + 1)
    return @inbounds chain[end] = element
end

# Heap stuff
function DataStructures.heapify!(chain::Chain, ordering::Base.Ordering)
    heapify!(chain.elements, ordering)
    return chain
end
function DataStructures.heappush!(chain::Chain, val, ordering::Base.Ordering)
    return heappush!(chain.elements, convert(_internal_eltype(chain), val), ordering)
end
function DataStructures.heappush!(
    chain::Chain{F,S}, (s, c)::Tuple{S,F}, ordering::Base.Ordering
) where {F,S}
    element = ChainElement{S,F}(s, c)
    return heappush!(chain, element, ordering)
end
function DataStructures.heappop!(chain::Chain, ordering::Base.Ordering)
    if isempty(chain)
        return nothing
    else
        heap = chain.elements
        E = eltype(chain)

        top = convert(E, heappop!(heap, ordering))
        while !isempty(heap)
            if iszero(top)
                top = convert(E, heappop!(heap, ordering))
            elseif convert(E, first(heap)) == top
                top += convert(E, heappop!(heap, ordering))
            else
                break
            end
        end
        return iszero(top) ? nothing : top
    end
end
function heapmove!(dst::Chain{F,S}, chain::Chain{F,S}, ordering::Base.Ordering) where {F,S}
    while (pivot = heappop!(chain, ordering)) ≢ nothing
        push!(dst, pivot)
    end
    return dst
end
function heapmove!(chain::Chain{F,S}, ordering::Base.Ordering) where {F,S}
    return heapmove!(Chain{F,S}(), chain, ordering)
end

_no_ordering_error() = error("no ordering given")
DataStructures.heapify!(::Chain) = _no_ordering_error()
DataStructures.heappop!(::Chain) = _no_ordering_error()
DataStructures.heappush!(::Chain, _) = _no_ordering_error()
heapmove!(::Chain) = _no_ordering_error()
heapmove!(::Chain, ::Chain) = _no_ordering_error()

# Other stuff
function clean!(chain::Chain{F}, ordering::Base.Ordering, factor=one(F)) where {F}
    @inbounds if !isempty(chain)
        sort!(chain.elements; alg=QuickSort, order=ordering)
        i = 0
        current = chain[1]
        for j in 2:length(chain)
            top = chain[j]
            if top == current
                current += top
            else
                if !iszero(current)
                    i += 1
                    chain[i] = current * factor
                end
                current = top
            end
        end
        if !iszero(current)
            i += 1
            chain[i] = current * factor
        end
        resize!(chain, i)
    end
    return chain
end

"""
    is_cocycle(filtration, chain::Chain, time)

Check whether a `chain` is a cocycle in `filtration` at specified `time`.
"""
function is_cocycle(filtration, chain::Chain{F,S}, time) where {F,S}
    buffer = Chain{F,simplex_type(filtration, dim(S) + 1)}()
    for (simplex, coefficient) in chain
        for cofacet in coboundary(filtration, simplex)
            if birth(cofacet) < time
                heappush!(buffer, (cofacet, coefficient), Base.Order.Forward)
            end
        end
    end
    return isnothing(heappop!(buffer, Base.Order.Forward))
end


text:
# This file is a part of SimilaritySearch.jl

import Base: push!
export ParallelExhaustiveSearch, search

struct ParallelExhaustiveSearch{DistanceType<:SemiMetric,DataType<:AbstractDatabase} <: AbstractSearchContext
    dist::DistanceType
    db::DataType
    lock::Threads.SpinLock
end

ParallelExhaustiveSearch(dist::SemiMetric, db::AbstractVecOrMat) = ParallelExhaustiveSearch(dist, convert(AbstractDatabase, db))
ParallelExhaustiveSearch(dist::SemiMetric, db::AbstractDatabase) = ParallelExhaustiveSearch(dist, db, Threads.SpinLock())

"""
    ParallelExhaustiveSearch(; dist=SqL2Distance(), db=VectorDatabase{Float32}())

Solves queries evaluating `dist` in parallel for the query and all elements in the dataset.
Note that this should not be used in conjunction with `searchbatch(...; parallel=true)` since they will compete for resources.
"""
function ParallelExhaustiveSearch(; dist=SqL2Distance(), db=VectorDatabase{Float32}())
    ParallelExhaustiveSearch(dist, db, Threads.SpinLock())
end

getpools(index::ParallelExhaustiveSearch) = nothing
Base.copy(ex::ParallelExhaustiveSearch; dist=ex.dist, db=ex.db) = ParallelExhaustiveSearch(dist, db, Threads.SpinLock())

"""
    search(ex::ParallelExhaustiveSearch, q, res::KnnResult; pools=nothing)

Solves the query evaluating all items in the given query.
"""
function search(ex::ParallelExhaustiveSearch, q, res::KnnResult; pools=nothing)
    dist = ex.dist
    elock = ex.lock
    Threads.@threads for i in eachindex(ex)
        d = evaluate(dist, ex[i], q)
        try
            lock(elock)
            push!(res, i, d)
        finally
            unlock(elock)
        end
    end

    (res=res, cost=length(ex))
end

function Base.push!(ex::ParallelExhaustiveSearch, u)
    push!(ex.db, u)
end

max_stars_repo_name: sadit/SimilaritySearch

max_issues_repo_name: sadit/SimilaritySearch

max_forks_repo_name: sadit/SimilaritySearch


text:
################################################################################
# Roessler SRK for second order weak approx

struct DRI1ConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T
  a031::T
  a032::T

  a121::T
  a131::T
  #a132::T

  #a221::T
  #a231::T
  #a232::T

  b021::T
  b031::T
  #b032::T

  b121::T
  b131::T
  #b132::T

  b221::T
  b222::T
  b223::T
  b231::T
  b232::T
  b233::T

  α1::T
  α2::T
  α3::T

  c02::T2
  c03::T2

  #c11::T2
  c12::T2
  c13::T2

  #c21::T2
  #c22::T2
  #c23::T2

  beta11::T
  beta12::T
  beta13::T

  #beta21::T
  beta22::T
  beta23::T

  beta31::T
  beta32::T
  beta33::T

  #beta41::T
  beta42::T
  beta43::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function DRI1ConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  a021 = convert(T, 1//2)
  a031 = convert(T, -1)
  a032 = convert(T, 2)

  a121 = convert(T, 342//491)
  a131 = convert(T, 342//491)

  b021 = convert(T, (6-sqrt(6))/10)
  b031 = convert(T, (3+2*sqrt(6))/5)

  b121 = convert(T, 3*sqrt(38//491))
  b131 = convert(T, -3*sqrt(38/491))

  b221 = convert(T, -214//513*sqrt(1105//991))
  b222 = convert(T, -491//513*sqrt(221//4955))
  b223 = convert(T, -491//513*sqrt(221//4955))
  b231 = convert(T, 214//513*sqrt(1105//991))
  b232 = convert(T, 491//513*sqrt(221//4955))
  b233 = convert(T, 491//513*sqrt(221//4955))

  α1 = convert(T, 1//6)
  α2 = convert(T, 2//3)
  α3 = convert(T, 1//6)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 1)

  c12 = convert(T2, 342//491)
  c13 = convert(T2, 342//491)

  beta11 = convert(T, 193//684)
  beta12 = convert(T, 491//1368)
  beta13 = convert(T, 491//1368)

  beta22 = convert(T, 1//6*sqrt(491//38))
  beta23 = convert(T, -1//6*sqrt(491//38))

  beta31 = convert(T, -4955//7072)
  beta32 = convert(T, 4955//14144)
  beta33 = convert(T, 4955//14144)

  beta42 = convert(T, -1//8*sqrt(4955//221))
  beta43 = convert(T, 1//8*sqrt(4955//221))

  NORMAL_ONESIX_QUANTILE = convert(T, -0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::DRI1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  DRI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI1ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 2//3) # convert(T, 2//3)
  a031 = convert(T, -1//3)
  a032 = convert(T, 1)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//4)
  α2 = convert(T, 1//2)
  α3 = convert(T, 1//4)

  c02 = convert(T2, 2//3)
  c03 = convert(T2, 2//3)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI3ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 1//4)
  a032 = convert(T, 1//4)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, (3-2*sqrt(6))/5)
  b031 = convert(T, (6+sqrt(6))/10)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//6)
  α2 = convert(T, 1//6)
  α3 = convert(T, 2//3)

  c02 = convert(T2, 1)
  c03 = convert(T2, 1//2)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI3,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI5ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 25//144)
  a032 = convert(T, 35//144)

  a121 = convert(T, 1//4)
  a131 = convert(T, 1//4)

  b021 = convert(T, 1//3)
  b031 = convert(T, -5//6)

  b121 = convert(T, 1//2)
  b131 = convert(T, -1//2)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//10)
  α2 = convert(T, 3//14)
  α3 = convert(T, 24//35)

  c02 = convert(T2, 1)
  c03 = convert(T2, 5//12)

  c12 = convert(T2, 1//4)
  c13 = convert(T2, 1//4)

  beta11 = convert(T, 1)
  beta12 = convert(T, -1)
  beta13 = convert(T, -1)

  beta22 = convert(T, 1)
  beta23 = convert(T, -1)

  beta31 = convert(T, 1//2)
  beta32 = convert(T, -1//4)
  beta33 = convert(T, -1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI5,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI5ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RI6ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 0)
  a032 = convert(T, 0)

  a121 = convert(T, 1)
  a131 = convert(T, 1)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, 1)
  b131 = convert(T, -1)

  b221 = convert(T, 1)
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -1)
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)
  α3 = convert(T, 0)

  c02 = convert(T2, 1)
  c03 = convert(T2, 0)

  c12 = convert(T2, 1)
  c13 = convert(T2, 1)

  beta11 = convert(T, 1//2)
  beta12 = convert(T, 1//4)
  beta13 = convert(T, 1//4)

  beta22 = convert(T, 1//2)
  beta23 = convert(T, -1//2)

  beta31 = convert(T, -1//2)
  beta32 = convert(T, 1//4)
  beta33 = convert(T, 1//4)

  beta42 = convert(T, 1//2)
  beta43 = convert(T, -1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RI6,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RI6ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI2WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1)
  a031 = convert(T, 0)
  a032 = convert(T, 0)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, 1)
  b031 = convert(T, 0)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)
  α3 = convert(T, 0)

  c02 = convert(T2, 1)
  c03 = convert(T2, 0)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI2WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI2WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI3WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1//2)
  a031 = convert(T, 0)
  a032 = convert(T, 3//4)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, (9-2*sqrt(15))/14)
  b031 = convert(T, (18+3*sqrt(15))/28)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 2//9)
  α2 = convert(T, 1//3)
  α3 = convert(T, 4//9)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 3//4)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI3WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI3WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RDI4WMConstantCache(T::Type, T2::Type)
  a021 = convert(T, 1//2)
  a031 = convert(T, -1)
  a032 = convert(T, 2)

  a121 = convert(T, 2//3)
  a131 = convert(T, 2//3)

  b021 = convert(T, (6-sqrt(6))/10)
  b031 = convert(T, (3+2*sqrt(6))/5)

  b121 = convert(T, sqrt(2//3))
  b131 = convert(T, -sqrt(2//3))

  b221 = convert(T, sqrt(2))
  b222 = convert(T, 0)
  b223 = convert(T, 0)
  b231 = convert(T, -sqrt(2))
  b232 = convert(T, 0)
  b233 = convert(T, 0)

  α1 = convert(T, 1//6)
  α2 = convert(T, 2//3)
  α3 = convert(T, 1//6)

  c02 = convert(T2, 1//2)
  c03 = convert(T2, 1)

  c12 = convert(T2, 2//3)
  c13 = convert(T2, 2//3)

  beta11 = convert(T, 1//4)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)

  beta22 = convert(T, sqrt(6)/4)
  beta23 = convert(T, -sqrt(6)/4)

  beta31 = convert(T, -1//4)
  beta32 = convert(T, 1//8)
  beta33 = convert(T, 1//8)

  beta42 = convert(T, sqrt(2)/4)
  beta43 = convert(T, -sqrt(2)/4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  DRI1ConstantCache(a021,a031,a032,a121,a131,b021,b031,b121,b131,b221,b222,b223,b231,b232,b233,α1,α2,α3,c02,c03,c12,c13,beta11,beta12,beta13,beta22,beta23,beta31,beta32,beta33,beta42,beta43,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI4WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI4WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct DRI1Cache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uhat::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType
  g2::Vector{rateNoiseType}
  g3::Vector{rateNoiseType}

  k1::rateType
  k2::rateType
  k3::rateType

  H02::uType
  H03::uType
  H12::Vector{uType}
  H13::Vector{uType}
  H22::Vector{uType}
  H23::Vector{uType}

  tmp1::possibleRateType
  tmpg::rateNoiseType

  tmp::uType
  resids::uType

end

@cache struct DRI1NMCache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType
  uhat::uType

  _dW::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType
  g2::rateNoiseType
  g3::rateNoiseType

  k1::rateType
  k2::rateType
  k3::rateType

  H02::uType
  H03::uType
  H12::rateNoiseType
  H13::rateNoiseType

  tmp1::possibleRateType
  tmpg::rateNoiseType

  tmp::uType
  resids::uType

end

function alg_cache(alg::DRI1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = DRI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::DRI1NM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = DRI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = zero(noise_rate_prototype)
  g3 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = zero(noise_rate_prototype)
  H13 = zero(noise_rate_prototype)
  #H22 = zero(noise_rate_prototype)
  #H23 = zero(noise_rate_prototype)

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1NMCache(u,uprev,uhat,_dW,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RI1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RI3,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI3ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RI5,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI5ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RI6,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RI6ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RDI2WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI2WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RDI3WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI3WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

function alg_cache(alg::RDI4WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI4WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  uhat = copy(uprev)
  tmp = zero(u)
  resids = zero(u)

  DRI1Cache(u,uprev,uhat,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,k1,k2,k3,H02,H03,H12,H13,H22,H23,tmp1,tmpg,tmp,resids)
end

# Roessler SRK for first order weak approx
struct RDI1WMConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T

  #a121::T

  #a221::T

  b021::T

  #b121::T

  #b221::T

  α1::T
  α2::T

  c02::T2

  #c11::T2
  #c12::T2

  #c21::T2
  #c22::T2

  beta11::T
  #beta12::T

  #beta21::T
  #beta22::T

  #beta31::T
  #beta32::T

  #beta41::T
  #beta42::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function RDI1WMConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  a021 = convert(T, 2//3)

  b021 = convert(T, 2//3)

  α1 = convert(T, 1//4)
  α2 = convert(T, 3//4)

  c02 = convert(T2, 2//3)

  beta11 = convert(T, 1)

  NORMAL_ONESIX_QUANTILE = convert(T, -0.9674215661017014)

  RDI1WMConstantCache(a021,b021,α1,α2,c02,beta11,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RDI1WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RDI1WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct RDI1WMCache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType

  k1::rateType
  k2::rateType

  H02::uType

  tmp1::possibleRateType
end

function alg_cache(alg::RDI1WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RDI1WMConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  k1 = zero(rate_prototype); k2 = zero(rate_prototype)

  H02 = zero(u)

  tmp1 = zero(rate_prototype)

  RDI1WMCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,H02,tmp1)
end

# Stratonovich sense

struct RSConstantCache{T,T2} <: StochasticDiffEqConstantCache
  # hard-coded version
  a021::T
  a031::T
  a032::T

  a131::T
  a141::T

  b031::T
  b032::T

  b121::T
  b131::T
  b132::T
  b141::T
  b142::T
  b143::T

  b221::T
  b231::T

  b331::T
  b332::T
  b341::T
  b342::T

  α1::T
  α2::T
  α3::T
  α4::T

  c02::T2
  c03::T2

  c13::T2
  c14::T2

  beta11::T
  beta12::T
  beta13::T
  beta14::T

  beta22::T
  beta23::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function RS1ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 0)
  a031 = convert(T, 1)
  a032 = convert(T, 0)

  a131 = convert(T, 1)
  a141 = convert(T, 1)

  b031 = convert(T, 1//4)
  b032 = convert(T, 3//4)

  b121 = convert(T, 2//3)
  b131 = convert(T, 1//12)
  b132 = convert(T, 1//4)
  b141 = convert(T, -5//4)
  b142 = convert(T, 1//4)
  b143 = convert(T, 2)

  b221 = convert(T, 1)
  b231 = convert(T, -1)

  b331 =  convert(T, 1//4)
  b332 =  convert(T, 3//4)
  b341 =  convert(T, 1//4)
  b342 =  convert(T, 3//4)

  α1 = convert(T, 0)
  α2 = convert(T, 0)
  α3 = convert(T, 1//2)
  α4 = convert(T, 1//2)

  c02 = convert(T2, 0)
  c03 = convert(T2, 1)

  c13 = convert(T2, 1)
  c14 = convert(T2, 1)

  beta11 = convert(T, 1//8)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)
  beta14 = convert(T, 1//8)

  beta22 = convert(T, -1//4)
  beta23 = convert(T, 1//4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  RSConstantCache(a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RS1,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RS1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function RS2ConstantCache(T::Type, T2::Type)
  a021 = convert(T, 2//3)
  a031 = convert(T, 1//6)
  a032 = convert(T, 1//2)

  a131 = convert(T, 1)
  a141 = convert(T, 1)

  b031 = convert(T, 1//4)
  b032 = convert(T, 3//4)

  b121 = convert(T, 2//3)
  b131 = convert(T, 1//12)
  b132 = convert(T, 1//4)
  b141 = convert(T, -5//4)
  b142 = convert(T, 1//4)
  b143 = convert(T, 2)

  b221 = convert(T, 1)
  b231 = convert(T, -1)

  b331 =  convert(T, 1//4)
  b332 =  convert(T, 3//4)
  b341 =  convert(T, 1//4)
  b342 =  convert(T, 3//4)

  α1 = convert(T, 1//4)
  α2 = convert(T, 1//4)
  α3 = convert(T, 1//2)
  α4 = convert(T, 0)

  c02 = convert(T2, 2//3)
  c03 = convert(T2, 2//3)

  c13 = convert(T2, 1)
  c14 = convert(T2, 1)

  beta11 = convert(T, 1//8)
  beta12 = convert(T, 3//8)
  beta13 = convert(T, 3//8)
  beta14 = convert(T, 1//8)

  beta22 = convert(T, -1//4)
  beta23 = convert(T, 1//4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  RSConstantCache(a021,a031,a032,a131,a141,b031,b032,b121,b131,b132,b141,b142,b143,b221,b231,b331,b332,b341,b342,α1,α2,α3,α4,c02,c03,c13,c14,beta11,beta12,beta13,beta14,beta22,beta23,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::RS2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  RS2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct RSCache{uType,randType,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType
  g2::Vector{rateNoiseType}
  g3::Vector{rateNoiseType}
  g4::Vector{rateNoiseType}

  k1::rateType
  k2::rateType
  k3::rateType

  H02::uType
  H03::uType
  H12::Vector{uType}
  H13::Vector{uType}
  H14::Vector{uType}
  H22::Vector{uType}
  H23::Vector{uType}

  tmp1::possibleRateType
  tmpg::rateNoiseType

end

function alg_cache(alg::RS1,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RS1ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  g4 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H14 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H14,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  RSCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg)
end

function alg_cache(alg::RS2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = RS2ConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
  g1 = zero(noise_rate_prototype)
  g2 = [zero(noise_rate_prototype) for k=1:m]
  g3 = [zero(noise_rate_prototype) for k=1:m]
  g4 = [zero(noise_rate_prototype) for k=1:m]
  k1 = zero(rate_prototype); k2 = zero(rate_prototype); k3 = zero(rate_prototype)

  H02 = zero(u)
  H03 = zero(u)
  H12 = Vector{typeof(u)}()
  H13 = Vector{typeof(u)}()
  H14 = Vector{typeof(u)}()
  H22 = Vector{typeof(u)}()
  H23 = Vector{typeof(u)}()

  for k=1:m
    push!(H12,zero(u))
    push!(H13,zero(u))
    push!(H14,zero(u))
    push!(H22,zero(u))
    push!(H23,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg = zero(noise_rate_prototype)

  RSCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,g2,g3,g4,k1,k2,k3,H02,H03,H12,H13,H14,H22,H23,tmp1,tmpg)
end

# PL1WM
struct PL1WMConstantCache{T} <: StochasticDiffEqConstantCache
  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function PL1WMConstantCache(T::Type)
  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  PL1WMConstantCache(NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::PL1WM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  PL1WMConstantCache(real(uBottomEltypeNoUnits))
end

struct PL1WMAConstantCache{T} <: StochasticDiffEqConstantCache
  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function PL1WMAConstantCache(T::Type)
  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  PL1WMAConstantCache(NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::PL1WMA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  PL1WMAConstantCache(real(uBottomEltypeNoUnits))
end

@cache struct PL1WMCache{uType,randType,rand2Type,MType1,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::rand2Type
  chi1::randType
  Ihat2::MType1

  tab::tabType

  g1::rateNoiseType

  k1::rateType
  k2::rateType

  Y::uType
  Yp::Vector{uType}
  Ym::Vector{uType}

  tmp1::possibleRateType
  tmpg1::rateNoiseType
  tmpg2::rateNoiseType
  Ulp::uType
  Ulm::uType
end

function alg_cache(alg::PL1WM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔZ)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔZ)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = PL1WMConstantCache(real(uBottomEltypeNoUnits))
  g1 = zero(noise_rate_prototype)

  k1 = zero(rate_prototype); k2 = zero(rate_prototype);

  Y = zero(u)
  Yp = Vector{typeof(u)}()
  Ym = Vector{typeof(u)}()

  for k=1:m
    push!(Yp,zero(u))
    push!(Ym,zero(u))
  end

  tmp1 = zero(rate_prototype)
  tmpg1 = zero(noise_rate_prototype)
  tmpg2 = zero(noise_rate_prototype)

  Ulp = zero(u)
  Ulm = zero(u)

  PL1WMCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,g1,k1,k2,Y,Yp,Ym,tmp1,tmpg1,tmpg2,Ulp,Ulm)
end

# additive noise
@cache struct PL1WMACache{uType,randType,tabType,rateNoiseType,rateType,possibleRateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  chi1::randType
  tab::tabType

  g1::rateNoiseType
  k1::rateType
  k2::rateType

  Y::uType
  tmp1::possibleRateType
end

function alg_cache(alg::PL1WMA,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    chi1 = zero(ΔW)
  end

  tab = PL1WMConstantCache(real(uBottomEltypeNoUnits))
  g1 = zero(noise_rate_prototype)

  k1 = zero(rate_prototype); k2 = zero(rate_prototype);

  Y = zero(u)

  tmp1 = zero(rate_prototype)

  PL1WMACache(u,uprev,_dW,chi1,tab,g1,k1,k2,Y,tmp1)
end

#NON
struct NONConstantCache{T} <: StochasticDiffEqConstantCache
  c01::T
  c02::T
  c03::T
  c04::T

  cj1::T
  cj2::T
  cj3::T
  cj4::T

  cjl2::T
  cjl3::T

  clj2::T
  clj3::T

  a0021::T
  a0032::T
  a0043::T

  aj021::T
  aj041::T

  a0j21::T
  a0j31::T
  a0j32::T
  a0j41::T

  ajj21::T
  ajj31::T
  ajj32::T
  ajj41::T
  ajj42::T
  ajj43::T

  ajl31::T
  ajl32::T
  ajl41::T
  ajl42::T

  ajljj31::T

  aljjl21::T
  aljjl31::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function NONConstantCache(T::Type)
  c01 = convert(T, 1//6)
  c02 = convert(T, 1//3)
  c03 = convert(T, 1//3)
  c04 = convert(T, 1//6)

  cj1 = convert(T, 1//8)
  cj2 = convert(T, 3//8)
  cj3 = convert(T, 3//8)
  cj4 = convert(T, 1//8)

  cjl2 = convert(T, 1//2)
  cjl3 = convert(T, -1//2)

  clj2 = convert(T, -1//2)
  clj3 = convert(T, 1//2)

  a0021 = convert(T, 1//2)
  a0032 = convert(T, 1//2)
  a0043 = convert(T, 1)

  aj021 = convert(T, 2)
  aj041 = convert(T, -2)

  a0j21 = convert(T, 1)
  a0j31 = convert(T, -9//8)
  a0j32 = convert(T, 9//8)
  a0j41 = convert(T, 1)

  ajj21 = convert(T, 2//3)
  ajj31 = convert(T, 1//12)
  ajj32 = convert(T, 1//4)
  ajj41 = convert(T, -5//4)
  ajj42 = convert(T, 1//4)
  ajj43 = convert(T, 2)

  ajl31 = convert(T, 1//4)
  ajl32 = convert(T, 3//4)
  ajl41 = convert(T, 1//4)
  ajl42 = convert(T, 3//4)

  ajljj31 = convert(T, 1)

  aljjl21 = convert(T, -1//2)
  aljjl31 = convert(T, 1//2)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  NONConstantCache(c01,c02,c03,c04,cj1,cj2,cj3,cj4,cjl2,cjl3,clj2,clj3,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,ajljj31,aljjl21,aljjl31,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::NON,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  NONConstantCache(real(uBottomEltypeNoUnits))
end

@cache struct NONCache{uType,randType,MType1,tabType,rateNoiseType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  gtmp::rateNoiseType
  ktmp::rateType

  Y100::uType
  Y200::uType
  Y300::uType
  Y400::uType
  Y1jajb::Array{uType}
  Y2jajb::Array{uType}
  Y3jajb::Array{uType}
  Y4jajb::Array{uType}

  tmpu::uType
  #tmpg::rateNoiseType

end

function alg_cache(alg::NON,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = NONConstantCache(real(uBottomEltypeNoUnits))

  gtmp = zero(noise_rate_prototype)
  ktmp = zero(rate_prototype)

  Y100 = zero(u)
  Y200 = zero(u)
  Y300 = zero(u)
  Y400 = zero(u)
  Y1jajb = [zero(u) for ja=1:m, jb=1:m]
  Y2jajb = [zero(u) for ja=1:m, jb=1:m]
  Y3jajb = [zero(u) for ja=1:m, jb=1:m]
  Y4jajb = [zero(u) for ja=1:m, jb=1:m]

  tmpu = zero(u)

  NONCache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,gtmp,ktmp,Y100,Y200,Y300,Y400,Y1jajb,Y2jajb,Y3jajb,Y4jajb,tmpu)

end

# COM
struct COMConstantCache{T} <: StochasticDiffEqConstantCache
  c01::T
  c02::T
  c03::T
  c04::T

  cj1::T
  cj2::T
  cj3::T
  cj4::T

  a0021::T
  a0032::T
  a0043::T

  aj021::T
  aj041::T

  a0j21::T
  a0j31::T
  a0j32::T
  a0j41::T

  ajj21::T
  ajj31::T
  ajj32::T
  ajj41::T
  ajj42::T
  ajj43::T

  ajl31::T
  ajl32::T
  ajl41::T
  ajl42::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function COMConstantCache(T::Type)
  c01 = convert(T, 1//6)
  c02 = convert(T, 1//3)
  c03 = convert(T, 1//3)
  c04 = convert(T, 1//6)

  cj1 = convert(T, 1//8)
  cj2 = convert(T, 3//8)
  cj3 = convert(T, 3//8)
  cj4 = convert(T, 1//8)

  a0021 = convert(T, 1//2)
  a0032 = convert(T, 1//2)
  a0043 = convert(T, 1)

  aj021 = convert(T, 2)
  aj041 = convert(T, -2)

  a0j21 = convert(T, 1)
  a0j31 = convert(T, -9//8)
  a0j32 = convert(T, 9//8)
  a0j41 = convert(T, 1)

  ajj21 = convert(T, 2//3)
  ajj31 = convert(T, 1//12)
  ajj32 = convert(T, 1//4)
  ajj41 = convert(T, -5//4)
  ajj42 = convert(T, 1//4)
  ajj43 = convert(T, 2)

  ajl31 = convert(T, 1//4)
  ajl32 = convert(T, 3//4)
  ajl41 = convert(T, 1//4)
  ajl42 = convert(T, 3//4)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  COMConstantCache(c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::COM,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  COMConstantCache(real(uBottomEltypeNoUnits))
end

@cache struct COMCache{uType,randType,tabType,rateNoiseType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType

  tab::tabType

  gtmp::rateNoiseType
  ktmp::rateType

  Y10::uType
  Y20::uType
  Y30::uType
  Y40::uType
  Y1j::rateNoiseType
  Y2j::rateNoiseType
  Y3j::rateNoiseType
  Y4j::rateNoiseType

  tmpu::uType
  tmpu2::uType
end

function alg_cache(alg::COM,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
  else
    _dW = zero(ΔW)
  end
  m = length(ΔW)
  tab = COMConstantCache(real(uBottomEltypeNoUnits))

  gtmp = zero(noise_rate_prototype)
  ktmp = zero(rate_prototype)

  Y10 = zero(u)
  Y20 = zero(u)
  Y30 = zero(u)
  Y40 = zero(u)
  Y1j = zero(noise_rate_prototype)
  Y2j = zero(noise_rate_prototype)
  Y3j = zero(noise_rate_prototype)
  Y4j = zero(noise_rate_prototype)

  tmpu = zero(u)
  tmpu2 = zero(u)

  COMCache(u,uprev,_dW,tab,gtmp,ktmp,Y10,Y20,Y30,Y40,Y1j,Y2j,Y3j,Y4j,tmpu,tmpu2)
end

# NON2
struct NON2ConstantCache{T} <: StochasticDiffEqConstantCache
  c01::T
  c02::T
  c03::T
  c04::T

  cj1::T
  cj2::T
  cj3::T
  cj4::T

  a0021::T
  a0032::T
  a0043::T

  aj021::T
  aj041::T

  a0j21::T
  a0j31::T
  a0j32::T
  a0j41::T

  ajj21::T
  ajj31::T
  ajj32::T
  ajj41::T
  ajj42::T
  ajj43::T

  ajl31::T
  ajl32::T
  ajl41::T
  ajl42::T

  # for non-commuting terms
  #ckj3::T
  #ckj4::T
  #akjjl32::T
  #akjjl42::T
  # are all expressed in terms of γ
  γ::T

  #quantile(Normal(),1/6)
  NORMAL_ONESIX_QUANTILE::T
end

function NON2ConstantCache(T::Type)
  c01 = convert(T, 1//6)
  c02 = convert(T, 1//3)
  c03 = convert(T, 1//3)
  c04 = convert(T, 1//6)

  cj1 = convert(T, 1//8)
  cj2 = convert(T, 3//8)
  cj3 = convert(T, 3//8)
  cj4 = convert(T, 1//8)

  a0021 = convert(T, 1//2)
  a0032 = convert(T, 1//2)
  a0043 = convert(T, 1)

  aj021 = convert(T, 2)
  aj041 = convert(T, -2)

  a0j21 = convert(T, 1)
  a0j31 = convert(T, -9//8)
  a0j32 = convert(T, 9//8)
  a0j41 = convert(T, 1)

  ajj21 = convert(T, 2//3)
  ajj31 = convert(T, 1//12)
  ajj32 = convert(T, 1//4)
  ajj41 = convert(T, -5//4)
  ajj42 = convert(T, 1//4)
  ajj43 = convert(T, 2)

  ajl31 = convert(T, 1//4)
  ajl32 = convert(T, 3//4)
  ajl41 = convert(T, 1//4)
  ajl42 = convert(T, 3//4)

  γ = convert(T, 1)

  NORMAL_ONESIX_QUANTILE = convert(T,-0.9674215661017014)

  NON2ConstantCache(c01,c02,c03,c04,cj1,cj2,cj3,cj4,a0021,a0032,a0043,aj021,aj041,a0j21,a0j31,a0j32,a0j41,ajj21,ajj31,ajj32,ajj41,ajj42,ajj43,ajl31,ajl32,ajl41,ajl42,γ,NORMAL_ONESIX_QUANTILE)
end

function alg_cache(alg::NON2,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  NON2ConstantCache(real(uBottomEltypeNoUnits))
end

@cache struct NON2Cache{uType,randType,MType1,tabType,rateNoiseType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  _dW::randType
  _dZ::randType
  chi1::randType
  Ihat2::MType1

  tab::tabType

  gtmp::rateNoiseType
  gtmp1::rateNoiseType
  ktmp::rateType

  Y10::uType
  Y20::uType
  Y30::uType
  Y40::uType
  Y1j::rateNoiseType
  Y2j::rateNoiseType
  Y3j::rateNoiseType
  Y4j::rateNoiseType

  tmpu::uType
  tmpu2::uType
end

function alg_cache(alg::NON2,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})
  if typeof(ΔW) <: Union{SArray,Number}
    _dW = copy(ΔW)
    _dZ = copy(ΔW)
    chi1 = copy(ΔW)
  else
    _dW = zero(ΔW)
    _dZ = zero(ΔW)
    chi1 = zero(ΔW)
  end
  m = length(ΔW)
  Ihat2 = zeros(eltype(ΔW), m, m)
  tab = NON2ConstantCache(real(uBottomEltypeNoUnits))

  gtmp = zero(noise_rate_prototype)
  gtmp1 = zero(noise_rate_prototype)
  ktmp = zero(rate_prototype)

  Y10 = zero(u)
  Y20 = zero(u)
  Y30 = zero(u)
  Y40 = zero(u)
  Y1j = zero(noise_rate_prototype)
  Y2j = zero(noise_rate_prototype)
  Y3j = zero(noise_rate_prototype)
  Y4j = zero(noise_rate_prototype)

  tmpu = zero(u)
  tmpu2 = zero(u)

  NON2Cache(u,uprev,_dW,_dZ,chi1,Ihat2,tab,gtmp,gtmp1,ktmp,Y10,Y20,Y30,Y40,Y1j,Y2j,Y3j,Y4j,tmpu,tmpu2)
end

# SIE / SME methods

struct SIESMEConstantCache{T,T2} <: StochasticDiffEqConstantCache
  α1::T
  α2::T

  γ1::T

  λ1::T
  λ2::T
  λ3::T

  µ1::T
  µ2::T
  µ3::T

  µ0::T2
  µbar0::T2

  λ0::T
  λbar0::T

  ν1::T
  ν2::T

  β2::T
  β3::T

  δ2::T
  δ3::T
end

function SIEAConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)

  γ1 = convert(T, 1//2)

  λ1 = convert(T, 1//4)
  λ2 = convert(T, -1//4)
  λ3 = convert(T, 1//4)

  µ1 = convert(T, 1//4)
  µ2 = convert(T, 1//4)
  µ3 = convert(T, -1//4)

  µ0 = convert(T2, 1//1)
  µbar0 = convert(T2, 1//1)

  λ0 = convert(T, 1//1)
  λbar0 = convert(T, 1//1)

  ν1 = convert(T, 1//1)
  ν2 = convert(T, 0)

  β2 = convert(T, 1//1)
  β3 = convert(T, 0)

  δ2 = convert(T, -1//1)
  δ3 = convert(T, 0)

  SIESMEConstantCache(α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3)
end

function alg_cache(alg::SIEA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SIEAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function SMEAConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  α1 = convert(T, 0)
  α2 = convert(T, 1//1)

  γ1 = convert(T, 1//2)

  λ1 = convert(T, 1//4)
  λ2 = convert(T, -1//4)
  λ3 = convert(T, 1//4)

  µ1 = convert(T, 1//4)
  µ2 = convert(T, 1//4)
  µ3 = convert(T, -1//4)

  µ0 = convert(T2, 1//2)
  µbar0 = convert(T2, 1//1)

  λ0 = convert(T, 1//2)
  λbar0 = convert(T, 1//1)

  ν1 = convert(T, (2-sqrt(6))/4)
  ν2 = convert(T, sqrt(6)/12)

  β2 = convert(T, 1//1)
  β3 = convert(T, 0)

  δ2 = convert(T, -1//1)
  δ3 = convert(T, 0)

  SIESMEConstantCache(α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3)
end

function alg_cache(alg::SMEA,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SMEAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function SIEBConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  α1 = convert(T, 1//2)
  α2 = convert(T, 1//2)

  γ1 = convert(T, -1//5)

  λ1 = convert(T, 3//5)
  λ2 = convert(T, 3//2)
  λ3 = convert(T, -1//2)

  µ1 = convert(T, 3//5)
  µ2 = convert(T, -3//2)
  µ3 = convert(T, 1//2)

  µ0 = convert(T2, 1//1)
  µbar0 = convert(T2, 5//12)

  λ0 = convert(T, 1//1)
  λbar0 = convert(T, 5//12)

  ν1 = convert(T, 1//1)
  ν2 = convert(T, 0)

  β2 = convert(T, 0)
  β3 = convert(T, -1//6)

  δ2 = convert(T, 0)
  δ3 = convert(T, 1//6)

  SIESMEConstantCache(α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3)
end

function alg_cache(alg::SIEB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SIEBConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

function SMEBConstantCache(::Type{T}, ::Type{T2}) where {T,T2}

  α1 = convert(T, 0)
  α2 = convert(T, 1//1)

  γ1 = convert(T, -1//5)

  λ1 = convert(T, 3//5)
  λ2 = convert(T, 3//2)
  λ3 = convert(T, -1//2)

  µ1 = convert(T, 3//5)
  µ2 = convert(T, -3//2)
  µ3 = convert(T, 1//2)

  µ0 = convert(T2, 1//2)
  µbar0 = convert(T2, 5//12)

  λ0 = convert(T, 1//2)
  λbar0 = convert(T, 5//12)

  ν1 = convert(T, (2-sqrt(6))/4)
  ν2 = convert(T, sqrt(6)/12)

  β2 = convert(T, 0)
  β3 = convert(T, -1//6)

  δ2 = convert(T, 0)
  δ3 = convert(T, 1//6)

  SIESMEConstantCache(α1,α2,γ1,λ1,λ2,λ3,µ1,µ2,µ3,µ0,µbar0,λ0,λbar0,ν1,ν2,β2,β3,δ2,δ3)
end

function alg_cache(alg::SMEB,prob,u,ΔW,ΔZ,p,rate_prototype,noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{false}})
  SMEBConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))
end

@cache struct SIESMECache{uType,randType,tabType,rateNoiseType,rateType} <: StochasticDiffEqMutableCache
  u::uType
  uprev::uType

  W2::randType
  W3::randType

  tab::tabType

  k0::rateType
  k1::rateType

  g0::rateNoiseType
  g1::rateNoiseType
  g2::rateNoiseType

  #tmp1::possibleRateType
  tmpu::uType
end

function alg_cache(alg::SIEA,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  if typeof(ΔW) <: Union{SArray,Number}
    W2 = copy(ΔW)
    W3 = copy(ΔW)
  else
    W2 = zero(ΔW)
    W3 = zero(ΔW)
  end

  k0 = zero(u)
  k1 = zero(u)

  g0 = zero(noise_rate_prototype)
  g1 = zero(noise_rate_prototype)
  g2 = zero(noise_rate_prototype)

  tmpu = zero(u)

  tab = SIEAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))

  SIESMECache(u,uprev,W2,W3,tab,k0,k1,g0,g1,g2,tmpu)
end

function alg_cache(alg::SMEA,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  if typeof(ΔW) <: Union{SArray,Number}
    W2 = copy(ΔW)
    W3 = copy(ΔW)
  else
    W2 = zero(ΔW)
    W3 = zero(ΔW)
  end

  k0 = zero(u)
  k1 = zero(u)

  g0 = zero(noise_rate_prototype)
  g1 = zero(noise_rate_prototype)
  g2 = zero(noise_rate_prototype)

  tmpu = zero(u)

  tab = SMEAConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))

  SIESMECache(u,uprev,W2,W3,tab,k0,k1,g0,g1,g2,tmpu)
end

function alg_cache(alg::SIEB,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  if typeof(ΔW) <: Union{SArray,Number}
    W2 = copy(ΔW)
    W3 = copy(ΔW)
  else
    W2 = zero(ΔW)
    W3 = zero(ΔW)
  end

  k0 = zero(u)
  k1 = zero(u)

  g0 = zero(noise_rate_prototype)
  g1 = zero(noise_rate_prototype)
  g2 = zero(noise_rate_prototype)

  tmpu = zero(u)

  tab = SIEBConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))

  SIESMECache(u,uprev,W2,W3,tab,k0,k1,g0,g1,g2,tmpu)
end

function alg_cache(alg::SMEB,prob,u,ΔW,ΔZ,p,rate_prototype,
                   noise_rate_prototype,jump_rate_prototype,uEltypeNoUnits,
                   uBottomEltypeNoUnits,tTypeNoUnits,uprev,f,t,dt,::Type{Val{true}})

  if typeof(ΔW) <: Union{SArray,Number}
    W2 = copy(ΔW)
    W3 = copy(ΔW)
  else
    W2 = zero(ΔW)
    W3 = zero(ΔW)
  end

  k0 = zero(u)
  k1 = zero(u)

  g0 = zero(noise_rate_prototype)
  g1 = zero(noise_rate_prototype)
  g2 = zero(noise_rate_prototype)

  tmpu = zero(u)

  tab = SMEBConstantCache(real(uBottomEltypeNoUnits), real(tTypeNoUnits))

  SIESMECache(u,uprev,W2,W3,tab,k0,k1,g0,g1,g2,tmpu)
end


text:
#!/usr/bin/env julia

## =========================== ## ## =========================== ##
## 	                  Minkowski Dimension                      	 ##
##	                Ollin Demian Langle Chimal					 ##
##	                         						             ##
##                         10 / 08 / 2018	    		         ##
## =========================== ## ## =========================== ##

using LibGEOS

# An AbstractType with an empty property
#  mutable struct box_2d
#      squares::Polygon
#      empty::Bool
#  end

"""
    gen_box(start,Δx,Δy,empty)

Generates a rectangle with the given start and the Δ in every axis.
Can be either empty or not.
"""
function gen_box(start::Array{Float64,1},Δx::Float64,Δy::Float64,empty::Bool)
    x = start[1]
    y = start[2]
    nx = x + Δx
    ny = y + Δy
    polyg = readgeom("POLYGON(($x $y,$nx $y,$nx $ny,$x $ny,$x $y))")
    box = box_2d(polyg,empty)
end

"""
    make_grid(min_x,max_y,min_y,max_y,Nx,Ny,empty)

Generates a grid of boxes given by the number of divisions in every axis.
Needs a start and end point. Currently you can only make them all empty,
not empty or randomize them.
"""
function make_grid(min_x,max_x,min_y,max_y,Nx,Ny,empty)
    Δx = (max_x-min_x)/Nx
    Δy = (max_y-min_y)/Ny
    tN = Nx*Ny
    boxes = Array{box_2d}(undef,tN)
    x = min_x
    y = min_y
    box = 1
    for j ∈ 1:Ny
        for i ∈ 1:Nx
            boxes[box] = gen_box([x,y],Δx,0.,empty)
            x = x + Δx
            box += 1
        end
        y = y + Δy
        x = min_x
    end
    return boxes,Δx,Δy
end

"""
    make_square_grid(min_x,max_x,min_y,max_y,Δ,empty)

Generates a grid of square boxes given by the Δ which is the box size length.
Needs a start and end point. Currently you can only make them all empty,
not empty or randomize them.
"""
function make_square_grid(min_x,max_x,min_y,max_y,Δ,empty)
    Nx = ceil((max_x-min_x)/Δ)
    Ny = ceil((max_y-min_y)/Δ)
    tN = Int(Nx*Ny)
    boxes = Array{box_2d}(undef,tN)
    x = min_x
    y = min_y
    box = 1
    for j ∈ 1:Ny
        for i ∈ 1:Nx
            boxes[box] = gen_box([x,y],Δ,Δ,empty)
            x = x + Δ
            box += 1
        end
        y = y + Δ
        x = min_x
    end
    return boxes
end

"""
    inside_box(point,box)

Checks if there's a point inside a box and if there is change the state of the box
to not empty.
"""
function inside_box(point,box)
    if within(point,box.squares)
        box.empty = false
    end
    return box.empty
end


text:
using Documenter
using Spectra

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

makedocs(sitename = "Spectra.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true",
    ),
    authors = "Miles Lucas and contributors.",
    linkcheck = !("skiplinks" in ARGS),
    modules = [Spectra],
    pages = [
        "Home" => "index.md",
        "spectrum.md",
        "transforms.md",
        "fitting.md",
        "analysis.md",
        "contrib.md",
    ],
    # strict = true,
)

deploydocs(repo = "github.com/JuliaAstro/Spectra.jl.git")


text:
if OS_NAME == :Windows
#	println("Compiling sigma.c...")
#	wd = pwd()
#	cd(dirname(Base.source_path()))
#	run(`gcc -fPIC -O3 -c sigma.c`)
#	println("Linking libsigma...")
#	run(`gcc -shared sigma.o -o libsigma.dll`)
#	cd(wd)
	println("Windows: skipping binary deps, please compile manually.")
elseif OS_NAME == :Darwin
	println("Darwin: skipping binary deps, using julia fallback.")
else
	println("Compiling sigma.c...")
	wd = pwd()
	cd(dirname(Base.source_path()))
	run(`gcc -fPIC -O3 -c sigma.c`)
	println("Linking libsigma...")
	run(`gcc -shared sigma.o -o libsigma.so`)
	cd(wd)
end