text:
# This file is part of EndemicModel project, which is licensed under BDS-3-Clause.
# See file LICENSE.md for more information.

# Parameter limits

const MIN_μ = 0.001
const MAX_μ = 0.300

# Source: https://journals.lww.com/cmj/Fulltext/2020/05050/Persistence_and_clearance_of_viral_RNA_in_2019.6.aspx
# Average infeccious period of 9.5 (6.0 to 11.0) days
# Source: https://www.sciencedirect.com/science/article/pii/S0163445320301195
# Average infeccious period of 11.0 (10.0 to 12.0) days
const MIN_α = 1 / 12.0
const MAX_α = 1 / 6.0

# Source: https://www.acpjournals.org/doi/10.7326/M20-0504
# Incubation period 5.1 (4.5 to 5.8) days
const MIN_γ = 1 / 5.8
const MAX_γ = 1 / 4.5

const MIN_β = 0.001
const MAX_β = 0.5

const ε = 1e-2

phuber_loss(a::Real, δ::Real) = δ^2 * (sqrt(1 + (a / δ)^2) - 1)
diff_phuber_loss(a::Real, δ::Real) = (sqrt(1 + (a / δ)^2)^-1) * a
phuber_loss(a::Vector{<:Real}, δ::Real) = δ^2 * (sqrt(1 + sum(a .* a) / δ^2) - 1)
diff_phuber_loss(a::Vector{<:Real}, δ::Real) = (sqrt(1 + (sum(a .* a) / δ^2))^-1) * a

l2_loss(a::Real) = a^2
diff_l2_loss(a::Real) = 2 * a
l2_loss(a::Vector{<:Real}) = sum(a .* a)
diff_l2_loss(a::Vector{<:Real}) = 2 * a

function phuber_loss(δ::Real)
    function _phuber_loss(a)
        phuber_loss(a, δ)
    end
    _phuber_loss
end

function diff_phuber_loss(δ::Real)
    function _diff_phuber_loss(a)
        diff_phuber_loss(a, δ)
    end
    _diff_phuber_loss
end

function model_loss(
    model::SEIRModel;
    loss_func::Function = phuber_loss(1.0),
    model_days = 15,
    kwargs...,
    # diff_loss_func::Function = diff_phuber_loss(1.0),
)
    params = parameters(model)
    ng = model.ngroups
    data = realdata(model)
    initial_date = default_kwarg(model, :initial_date, today() - Day(model_days))
    ini = findfirst(isequal(initial_date), data.date)
    if isnothing(ini)
        throw(ErrorException("Initial date not found in model data", model, initial_date))
    end
    iend = nrow(data)
    data = data[ini:iend, :]
    ndays = nrow(data)
    mdata = modeldata(model)
    nrow(mdata) < ndays && return Inf

    @debug(
        "Computing maximum loss for model",
        model,
        recovered = Tuple(data.recovered),
        deaths = Tuple(data.deaths),
        active = Tuple(data.active),
    )
    #rec_tot = loss_func(data.recovered .- mean(data.recovered))
    dth_tot = loss_func(data.deaths .- mean(data.deaths))
    act_tot = loss_func(data.active .- mean(data.active))
    max_loss = act_tot + dth_tot # + rec_tot

    #rec_loss = Float[]
    act_loss_vec = data.active .- mdata.active[1:ndays]
    dth_loss_vec = data.deaths .- mdata.deaths[1:ndays]
    #rec_loss_vec = data.recovered .- mdata.recovered[1:ndays]

    act_loss, dth_loss = loss_func(act_loss_vec), loss_func(dth_loss_vec)
    @debug(
        "Computing loss",
        Tuple(act_loss_vec),
        Tuple(dth_loss_vec),
        act_loss,
        dth_loss,
        max_loss
    )
    loss_value = act_loss + dth_loss
    loss_value / max_loss
end

function optimize_parameters!(
    model::SEIRModel,
    params = (:M, :β, :E);
    model_days = 30,
    kwargs...,
)
    packed_params = intersect(SEIR_PARAMETERS, params)
    n = model.ngroups
    lastparams = nothing
    model_params = parameters(model)
    model_vars = variables(model)

    function _calc_loss(arg)
        if :E ∈ params
            newparams = arg[1:(end-n)]
            newE = arg[(end-n+1):end]
        else
            newparams = arg
        end
        newparams = unpack_params(
            newparams,
            model.ngroups;
            packed_params = packed_params,
            default_params = model_params,
        )
        parameters!(model, newparams)
        if :E ∈ params || :M ∈ params
            (S, E, I, R) = variables(model)
            if :E ∈ params
                E = newE
            end
            M = newparams[:M]
            if :M ∈ params
                S = M .- E .- I .- R
            end
            variables!(model, SEIRVariables(S, E, I, R))
        end
        modeldata!(
            model,
            to_dataframe(model; diff_columns = false, maxtime = model_days + 1, kwargs...),
        )
        model_loss(model; kwargs...)
    end

    maxM = fill((1 + ε) * sum(model_params[:M]), n)
    minM = (1 + ε) * (model_vars[:I] .+ model_vars[:R])
    maxp = (; M = maxM, β = MAX_β, γ = MAX_γ, α = MAX_α, μ = MAX_μ)
    minp = (; M = minM, β = MIN_β, γ = MIN_γ, α = MIN_α, μ = MIN_μ)

    upper = pack_params(maxp; packed_params = packed_params)
    lower = pack_params(minp; packed_params = packed_params)
    opt_params = pack_params(model_params; packed_params = packed_params)

    if :E ∈ params
        maxE, minE = maxM ./ 2, fill(0.0, n)
        append!(upper, maxE)
        append!(lower, minE)
        E0 = model_vars[:E]
        append!(opt_params, Float.(E0))
    end
    @debug "Initial parameter values:" opt_params, lower, upper

    opt = optimize(
        _calc_loss,
        lower,
        upper,
        opt_params,
        Fminbox(NelderMead()),
        Optim.Options(x_tol = 1 / 256),
    )
    modeldata!(model, to_dataframe(model))
    opt
end

optimize_parameters(model::SEIRModel; kwargs...) =
    optimize_parameters!(copy(model); kwargs...)

function optimize_parameters!(data::AbstractDict; kwargs...)
    result = DataDict()
    for (key, subdata) ∈ data
        @debug "Optimizing parameters." key _debuginfo(subdata)
        result[key] = optimize_parameters!(subdata; kwargs...)
    end
    result
end

function optimize_parameters(data::AbstractDict; kwargs...)
    result = DataDict()
    for (key, subdata) ∈ data
        @debug "Optimizing parameters." key _debuginfo(subdata)
        result[key] = optimize_parameters(subdata; kwargs...)
    end
    result
end

optimize_parameters(data::AbstractDatabase; kwargs...) =
    optimize_parameters(modeldict(data); kwargs...)

function optimize_parameters!(data::AbstractDatabase; kwargs...)
    result = optimize_parameters!(modeldict(data); kwargs...)
    parameters!(data)
    result
end

optimize_parameters(data; kwargs...) = data
optimize_parameters!(data; kwargs...) = data

function optimize_variables!(model::SEIRModel; kwargs...)
    n = model.ngroups
    lastparams = nothing
    model_vars, model_params = variables(model), parameters(model)
    function _calc_loss(vars)
        newvars = unpack_vars(vars)

        variables!(model, newvars)
        model_loss(model; kwargs...)
    end

    maxS = (1 + ε) .* model_params[:M]
    maxv = (; S = maxS, E = maxS, I = maxS, R = maxS)
    minS = fill(-0.01, n)
    minv = (; S = minS, E = minS, I = minS, R = minS)
    upper = pack_vars(maxv)
    lower = pack_vars(minv)
    minbox = Fminbox(NelderMead())
    opt_vars = pack_vars(model_vars)
    optimize(_calc_loss, lower, upper, opt_vars, minbox)
end

optimize_variables(model::SEIRModel; kwargs...) =
    optimize_variables!(copy(model); kwargs...)

function search_parameters(model::AbstractEndemicModel) end

function estimate_μ(data::AbstractDataFrame; ndays = 7, kwargs...)
    deaths, recovered = data.deaths, data.recovered
    deaths[end] == 0 && return [NaN]
    ind1 = findfirst(.!ismissing.(deaths) .& (deaths .> 0))
    ind2 = findfirst(.!ismissing.(recovered) .& (recovered .>= 100))
    (isnothing(ind1) || isnothing(ind2)) && return [NaN]

    ind = max(ind1, ind2, length(deaths) - ndays)
    dth, rec = deaths[ind:end], recovered[ind:end]
    vals = filter(x -> !ismissing(x) && isfinite(x) && x > 0, dth ./ (dth .+ rec))
    @debug(
        "Estimating μ = deaths / (deaths + recovered)",
        deaths = Tuple(dth),
        recovered = Tuple(rec),
        values = Tuple(vals)
    )
    isempty(vals) && return (NaN, Inf)
    val = mean(vals)
    if val >= 0.3 || val <= MIN_μ
        return NaN
    else
        return val
    end
end

estimate_μ(model::AbstractEndemicModel; kwargs...) = estimate_μ(realdata(model); kwargs...)

function estimate_α(data::AbstractDataFrame; ndays = 7, kwargs...)
    # Get numbers from last `ndays` days
    # ini = max(1, nrow(data) - ndays)
    # data = data[ini:end, :]
    # # Filter valid entries
    # data = @where(data, (.!ismissing.(:diff_closed)) .& (:active .> 0))
    # dR = data.diff_closed
    # I = data.active
    # isempty(I) && return (NaN, Inf)
    # @debug "Estimating α = dR / I." dR = Tuple(dR) I = Tuple(I)
    # vals = dR ./ I
    #
    # if any(val -> val >= 1 - ε, vals)
    #     return (MAX_α, StatsBase.std(vals; mean = MAX_α))
    # elseif any(val -> val <= ε, vals)
    #     return (MIN_α, StatsBase.std(vals; mean = MIN_α))
    # else
    #     val = min(MAX_α, max(MIN_α, geomean(vals .+ 1) - 1))
    #     return (val, StatsBase.std(vals; mean = val))
    # end

    # Source: https://journals.lww.com/cmj/Fulltext/2020/05050/Persistence_and_clearance_of_viral_RNA_in_2019.6.aspx
    # Average infeccious period of 9.5 (6.0 to 11.0) days
    # Source: https://www.sciencedirect.com/science/article/pii/S0163445320301195
    # Average infeccious period of 11.0 (10.0 to 12.0) days
    1 / 11.0
end

estimate_α(model::AbstractEndemicModel; kwargs...) = estimate_α(realdata(model); kwargs...)

# function _γ_root(d1, d2, d3, I, α)
#     a = -I .* d2 .+ d1 .^ 2 .- I .* α .* d1
#     b = d1 .* d2 .- I .* α .* d2
#     c = -I .* d3
#     Δ = b .^ 2 - 4 .* a .* c
#
#     vals = [
#         ai < 0 ? (-sqrt(abs(Δi)) - bi) / (2ai) : (sqrt(abs(Δi)) - bi) / (2ai)
#         for (ai, bi, Δi) ∈ zip(a, b, Δ)
#     ]
#
#     @debug(
#         "Estimating γ as the mean of max roots of equations `a * γ^2 + b * γ + c = 0`.",
#         a = Tuple(a),
#         b = Tuple(b),
#         c = Tuple(c),
#         Δ = Tuple(Δ),
#         values = Tuple(vals),
#     )
#
#     vals
# end

function estimate_γ(
    data::AbstractDataFrame;
    # ndays = 7,
    # α = estimate_α(data; kwargs...),
    kwargs...,
)
    # # Try to find at least 4 valid entries
    # ini = something(findlast(!ismissing, data.diff3_infected), 8) - 3
    # # Get numbers from last `ndays` days
    # ini = min(ini, max(1, nrow(data) - ndays))
    # data = data[ini:end, :]
    # # Filter valid entries
    # data = @where(data, .!ismissing.(:diff3_infected) .& .!ismissing.(:active))
    #
    # d1 = data.diff_infected
    # d2 = data.diff2_infected
    # d3 = data.diff3_infected
    #
    # I = data.active
    # isempty(I) && return (NaN, Inf)
    # vals = _γ_root(d1, d2, d3, I, α)
    #
    # if any(val -> val >= 1 - ε, vals)
    #     return (MAX_γ, StatsBase.std(vals; mean = MAX_γ))
    # elseif any(val -> val <= ε, vals)
    #     return (MIN_γ, StatsBase.std(vals; mean = MIN_γ))
    # else
    #     val = min(MAX_γ, max(MIN_γ, geomean(vals .+ 1) - 1))
    #     return (val, StatsBase.std(vals; mean = val))
    # end

    # Source: https://www.acpjournals.org/doi/10.7326/M20-0504
    # Incubation period 5.1 (4.5 to 5.8) days
    1 / 5.1
end

estimate_γ(model::AbstractEndemicModel; kwargs...) = estimate_γ(realdata(model); kwargs...)

function estimate_β(
    data::AbstractDataFrame;
    ndays = 7,
    α = estimate_α(data; kwargs...),
    γ = estimate_γ(data; α = α, kwargs...),
    kwargs...,
)
    # Get numbers from last `ndays` days
    ini = max(1, nrow(data) - ndays)
    data = data[ini:end, :]
    # Filter valid entries
    data = @where(data, .!ismissing.(:diff2_infected))

    if nrow(data) == 0
        return (mean([MIN_β, MAX_β]), (MAX_β - MIN_β) / 2)
    end

    d1 = data.diff_infected
    d2 = data.diff2_infected
    I = data.active
    M = data.estimated_population[1]
    R = data.closed

    @debug(
        "Estimating β = (d²(I + R) + γ * d(I + R)) / (γ * I + d(I + R)).",
        d1 = Tuple(d1),
        d2 = Tuple(d2),
        I = Tuple(I),
        γ = γ,
    )
    s_inv = M ./ (M .- I .- R)
    vals = s_inv .* (d2 .+ γ .* d1) ./ (γ .* I .+ d1)
    isempty(vals) && return (NaN, Inf)

    if any(val -> val >= 1 - ε, vals)
        return (1 - ε) * MAX_β
    elseif any(val -> val <= ε, vals)
        return (1 + ε) * MIN_β
    else
        val = (1 - ε) * MAX_β
        return val
    end

end

function estimate_exposed(
    data::AbstractDataFrame;
    μ = estimate_μ(data; kwargs...),
    α = estimate_α(data; kwargs...),
    γ = estimate_γ(data; α = α, kwargs...),
    kwargs...,
)
    d1 = data.diff_infected

    @debug "Estimating E = d(I + R) / γ." d1 = Tuple(d1) γ = Tuple(γ)
    E = d1 ./ γ
    E = [ismissing(elt) ? 0 : isfinite(elt) ? round(Int, elt) : 0 for elt ∈ E]
    max.(E, 0)
end

function estimate_exposed!(data::AbstractDataFrame; kwargs...)
    E = estimate_exposed(data; kwargs...)
    if "exposed" ∉ names(data)
        insertcols!(data, ncol(data) + 1, :exposed => E)
    else
        data.exposed = E
    end
    E
end

estimate_exposed(model::AbstractEndemicModel; kwargs...) =
    estimate_exposed(realdata(model); kwargs...)
estimate_exposed!(model::AbstractEndemicModel; kwargs...) =
    estimate_exposed!(realdata(model); kwargs...)

function SEIRModel(
    data::AbstractDataFrame;
    modeldata = nothing,
    minimum_infected = option(:minimum_infected),
    model_days = 30,
    optimize::Bool = true,
    kwargs...,
)
    if !isnothing(modeldata)
        return SEIRModel(
            data,
            modeldata;
            minimum_infected = minimum_infected,
            model_days = model_days,
            kwargs...,
        )
    end
    !hasproperty(data, :date) && return missing
    nrows = nrow(data)
    if any(data.date[i] + Day(1) != data.date[i+1] for i ∈ 1:nrows-1)
        # If dates are not sequential, give up
        return missing
    end

    @debug "Computing SEIR model for data" _debuginfo(data)

    numpeople = data.estimated_population[1]
    istart = findfirst(data.infected .>= minimum_infected)
    iend = findlast(!ismissing, data.recovered)
    (isnothing(istart) || isnothing(iend) || (iend <= istart)) && return missing
    data = data[istart:iend, :]

    idx = nrow(data) - model_days + 1
    idx <= 0 && return missing

    initial_date = data.date[idx]
    M = Float[numpeople]
    μ = estimate_μ(data; kwargs...)
    α = estimate_α(data; μ = μ, kwargs...)
    γ = estimate_γ(data; μ = μ, α = α, kwargs...)
    β = estimate_β(data; μ = μ, α = α, γ = γ, kwargs...)

    E_vec = estimate_exposed!(data; μ = μ, α = α, γ = γ, kwargs...)
    E = Float[E_vec[idx]]
    I = Float[data.active[idx]]
    R = Float[data.closed[idx]]
    S = M - E - I - R
    vars, params = SEIRVariables(S, E, I, R), SEIRParameters(M, [β], [γ], [α], [μ])

    model = SEIRModel(vars, params; realdata = data, initial_date = initial_date, kwargs...)

    if optimize
        @debug "Finding optimal parameters for model" model
        result = optimize_parameters!(model; model_days = model_days, kwargs...)
        @debug "Optimal parameters for model computed" result
    end

    model
end

function SEIRModel(
    data::AbstractDataFrame,
    modeldata::AbstractDataFrame;
    minimum_infected = option(:minimum_infected),
    kwargs...,
)
    initial_date = modeldata.date[1]
    M = data.estimated_population[1]

    E, I = modeldata.exposed, modeldata.active
    R = modeldata.recovered .+ modeldata.deaths
    S = M .- (E .+ I .+ R)

    dS = S[end] - S[end-1]
    dE = E[end] - E[end-1]
    dI = I[end] - I[end-1]
    dR = R[end] - R[end-1]

    β = mean(-dS / (I[end-1] * S[end-1] / M))
    γ = mean((dI + dR) / E[end-1])
    α = mean(dR / I[end-1])
    μ = modeldata.deaths[end] / (modeldata.recovered[end] + modeldata.deaths[end])

    SEIRModel(
        [S[1]],
        [E[1]],
        [I[1]],
        [R[1]],
        Float[M],
        [β],
        [γ],
        [α],
        [μ];
        initial_date = initial_date,
        realdata = data,
        modeldata = modeldata,
    )
end

_SEIRModel_exceptions = []

function SEIRModel!(
    destiny::AbstractDataDict,
    data::AbstractDict;
    modeldata = nothing,
    kwargs...,
)
    function _try(key, data)
        try
            submodel = isnothing(modeldata) ? nothing : get(modeldata, key, nothing)
            val = SEIRModel(data; modeldata = submodel, kwargs...)
            val isa AbstractDataDict && isempty(val) ? missing : val
        catch exception
            exception isa InterruptException && rethrow(exception)
            @debug "Exception during computation of SEIR model." exception
            push!(_SEIRModel_exceptions, (data, exception))
            missing
        end
    end
    for (key, subdata) ∈ data
        val = _try(key, subdata)
        !ismissing(val) && (destiny[key] = val)
    end
    destiny
end

SEIRModel(data::AbstractDict; kwargs...) =
    SEIRModel!(empty(data, Symbol, Any), data; kwargs...)

SEIRModel(database::AbstractDatabase; kwargs...) = SEIRModel(datadict(database); kwargs...)

function SEIRModel!(database::AbstractDatabase; kwargs...)
    modeldict!(database, SEIRModel(database; kwargs...))
    parameters!(database; kwargs...)
end

function parameters!(destiny::AbstractDict, data::AbstractDict; kwargs...)
    columns = [param => OptFloat[] for param ∈ SEIR_PARAMETERS]
    paramdf = DataFrame(
        :key => String[],
        :group_name => String[],
        columns...,
        :R0 => OptFloat[],
        :Rt => OptFloat[],
        :loss_7_days => OptFloat[],
        :loss_14_days => OptFloat[],
    )
    # Auxiliar value
    for (key, subdata) ∈ data
        @debug "Computing parameters." key _debuginfo(subdata)
        if subdata isa AbstractDict
            subdestiny = get!(destiny, key, empty(destiny))
            try
                parameters!(subdestiny, subdata; kwargs...)
            catch
            end
            continue
        end
        !(subdata isa AbstractEndemicModel) && continue

        n = ngroups(subdata)
        gnames = groupnames(subdata)
        params = pairs(parameters(subdata))
        vars = pairs(variables(subdata))
        loss7 = model_loss(subdata; ndays = 7, kwargs...)
        loss7 = isfinite(loss7) ? loss7 : missing
        loss14 = model_loss(subdata; ndays = 14, kwargs...)
        loss14 = isfinite(loss14) ? loss14 : missing

        E_last = subdata.realdata.exposed[end]
        I_last = subdata.realdata.infected[end]
        R_last = subdata.realdata.recovered[end]
        S_last = params[:M] .- E_last .- I_last .- R_last

        for (i, gname) ∈ enumerate(gnames)
            params = [k => (isfinite(v[i]) ? v[i] : missing) for (k, v) ∈ params]
            R0 = params[:β][i] / params[:α][i]
            R0 = isfinite(R0) ? R0 : missing
            Rt = ismissing(R0) ? R0 : R0 * S_last[i] / params[:M][i]
            moreparams = (; R0 = R0, Rt = Rt, loss_7_days = loss7, loss_14_days = loss14)
            row = (;
                key = string(key),
                group_name = string(prettify(gname)),
                params...,
                moreparams...,
            )
            @debug "Computed params for SEIR model" row
            push!(paramdf, row)
        end
    end
    sort!(paramdf, :key)
    !isempty(paramdf) && (destiny[:MODEL_PARAMETERS] = paramdf)
    destiny
end

parameters(data::AbstractDict; kwargs...) =
    parameters!(empty(data, Symbol, Any), data; kwargs...)
parameters!(data::AbstractDataDict; kwargs...) = parameters!(data, data; kwargs...)

parameters(data::AbstractDatabase; kwargs...) = parameters(modeldict(data); kwargs...)
parameters!(data::AbstractDatabase; kwargs...) = parameters!(modeldict(data); kwargs...)


text:
include("entity.jl")

function generatePopulation(n, genesLength)
   data = []
   for i in 1:n
       entity1 = generateEntity(genesLength)
       push!(data, entity1)
   end
   return data
end

function printPopulation(data)
   for i in 1:length(data)
      printEntity(data[i])
   end
end

function calculatePopulationFitness!(data)
   for i in 1:length(data)
      calculateFitness!(data[i])
   end
   sort!(data, by = d -> d.fitness, rev=false)
end

function crossoverPopulation(data, crossoverPoint)
    newData = []
    for i in 1:2:length(data)
       entity1 = deepcopy(data[i])
       entity2 = deepcopy(data[i+1])
       crossover!(entity1, entity2, crossoverPoint)
       push!(newData, entity1)
       push!(newData, entity2)
    end
    return newData
end

function mutatePopulation!(data, mutationPercentage)
   for i in 1:length(data)
      mutate!(data[i], mutationPercentage)
   end;
end

function selectPopulation(data, n)
   return copy(data[1:n])
end

ext: jl

max_stars_repo_name: pikacxe/Primenjeni-Algoritmi

max_issues_repo_name: pikacxe/Primenjeni-Algoritmi

max_forks_repo_name: pikacxe/Primenjeni-Algoritmi


text:
"constraint: `c^2 + d^2 <= a*b`"
function relaxation_complex_product_conic(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef)
    a_lb, a_ub = InfrastructureModels.variable_domain(a)
    b_lb, b_ub = InfrastructureModels.variable_domain(b)

    @assert (a_lb >= 0 && b_lb >= 0) || (a_ub <= 0 && b_ub <= 0)

    JuMP.@constraint(m, [a/sqrt(2), b/sqrt(2), c, 0] in JuMP.RotatedSecondOrderCone())
end


text:
task = @task begin
                 sleep(5)
                 println("Done")
             end

println(task)

# exec
schedule(task)
wait(task)

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
module MinimumViableTables

using Base: @pure, @propagate_inbounds, @_propagate_inbounds_meta, @_inline_meta,
            getproperty

import Base: size, axes, getindex, setindex!, show, similar, copy, filter, map, findall,
             getproperty, in, broadcast, map

import LinearAlgebra: cross, ×

using Unicode
using Indexing
#using DataStreams

export colnames, columns, accelerators, project, Project, accelerate, rename, Rename

export AbstractIndex, AbstractUniqueIndex, NoIndex, UniqueIndex, HashIndex, UniqueHashIndex,
       SortIndex, UniqueSortIndex

export Table, ProductArray, ProductTable, ×, cross

export Predicate
export IsEqual, IsLess, IsLessEqual, IsGreater, IsGreaterEqual, In
export Equals, LessThan, LessEqualThan, GreaterThan, GreaterEqualThan, Within

include("interval.jl")

include("util.jl")
include("project.jl")
include("rename.jl")
include("index.jl")
include("indices.jl")
include("table.jl")
include("rowoperations.jl")
include("columnoperations.jl")
include("accelerate.jl")
include("predicates.jl")

include("custom_arrays.jl")
include("producttable.jl")

include("show.jl")
#include("datastreams.jl")

end # module

# Fundamental relational algebra operations:
#
# Set and bag operations: union, setdiff, intersect, isunique, vcat (from Julia AbstractVector interface)
# Projections: project, Project{names} (also map)
# Selections: filter (accelerated when using filter predicates like IsEqual)
# Renames: rename, Rename{oldnames, newnames} (also map)
# Cartesian product: ProductTable - a matrix of named tuples?
#
# Add a new column based on existing columns: ?? (also map)
#
# The above also relate to joins with (abstract) relations... we can think about this...
#
# Some performance improvements:
#
# * The unique searching cases tend to return 0 or 1 elements. A "Nullable" like vector
#   (an immutable AbstractVector of length 0 or 1) could be a performance help in these
#   cases.
#
# * Perhaps use faster sorting (radix, etc)
#
# * Perhaps investigate performance of Dict
#
# * Primary keys (particulary lexicographical sort) and the ability to reorder a table
#   to suite one (or more) acceleration indexes (i.e. improve memory locality).
#
# * Think about memory locality and hash-based acceleration indexes.
#
# * Consider tables with row storage.
#
# * More views.


text:
using ECCodes

withmessages(ARGS[1]) do msg
    for k in sort(keys(msg))
        v=msg[k]
        try
            println("$(k): $(v)")
        catch
        end
    end
    println()
end


text:
# function posterior_2pl_quick(
#     parameters::Vector{Vector{Float64}},
#     responses::Vector{Float64},
#     X::Vector{Float64},
#     W::Vector{Float64}
#     )
#     return map( (x, w) ->
#             mapreduce( (pars, r) ->
#             begin
#                 if (r > 0)
#                     _sig_c(pars[1]*(x - pars[2]))
#                 else
#                     _sig_cplus(pars[1]*(x - pars[2]))
#                 end
#             end
#             ,
#             *,
#             parameters,
#             responses,
#             )*w,
#         X,
#         W
#         )
# end

# function posterior_2pl_quick!(
#     posteriors::Vector{Vector{Float64}},
#     i_index::Vector{Vector{Int64}},
#     parameters_vectors::Vector{Vector{Float64}},
#     responses::Vector{Vector{Float64}},
#     X::Vector{Float64},
#     W::Vector{Float64}
#     )
#     likelihood = zero(Float64)
#     for n = 1 : size(posteriors, 1)
#         p = posterior_2pl_quick(parameters_vectors[i_index[n]], responses[n], X, W)
#         normalizer = sum(p)
#         if normalizer > typemin(Float64)
#             posteriors[n] = p ./ normalizer
#             likelihood += _log_c(normalizer)
#         else
#             posteriors[n] = copy(p)
#         end
#     end
#     return likelihood
# end

# # function posterior_quick(
# #     responses::Vector{Union{Missing, Float64}},
# #     sig_phi::Matrix{Float64},
# #     W::Vector{Float64}
# #     )
# #     K = size(X, 1)
# #     post = ones(Float64, K)
# #     for k = 1:K
# #             post_k = 1.0
# #             for i in 1:size(responses, 1)
# #                 if (responses[i] > 0)
# #                     post_k *= sig_phi[i, k]
# #                 else
# #                     post_k *= 1 - sig_phi[i, k]
# #                 end
# #             end
# #             post[k] = post_k * W[k]
# #         end
# #     return post::Vector{Float64}
# # end


text:
function hash_store(A, hash_key)

    H = zeros(hash_key+1)
    for x in A
        k = Int(x % hash_key)
        while H[k+1] != 0
          k = (k + 1) % hash_key
        end
        H[k+1] = x
    end

    return H

end

function hash_search(ls, x, hash_key)

    k = Int(x % hash_key)
    while ls[k+1] != 0
        if ls[k+1] == x
            return k + 1
        else
            k = (k + 1) % hash_key
        end
    end

    return nothing

end

function hashSearch(A, hash_key, t)

    hash_array = hash_store(A, hash_key)
    hash_index = hash_search(hash_array, t, hash_key)

    return hash_array, hash_index

end

function main()

  # list to be sorted
  rng = MersenneTwister(1234);
  list_size = 10
  A = shuffle(rng, Vector(1:list_size))

  # Search
  search_num = 10
  hash_array, index_found = @time hashSearch(A, 11, search_num)

  print("Search for : ")
  println(search_num)

  print("Found index: ")
  print(index_found)
  println(" (in hash array)")

  print("Found value is ... ")
  println(hash_array[index_found])

end

if contains(@__FILE__, PROGRAM_FILE)
    main()
end

max_stars_repo_name: szmlb/algorithmsInJulia

max_issues_repo_name: szmlb/algorithmsInJulia

max_forks_repo_name: szmlb/algorithmsInJulia


text:
import Base.exp, Base./, Base.*, Base.-, Base.sum, Base.mod, Base.sin
#import NumericExtensions.fma

immutable Number8{T <: FloatingPoint}
    array :: Array{T,1}
    s0 :: T
    s1 :: T
    s2 :: T
    s3 :: T
    s4 :: T
    s5 :: T
    s6 :: T
    s7 :: T
end

immutable Number4 {T <: FloatingPoint}
    s0 :: T
    s1 :: T
    s2 :: T
    s3 :: T
end

function Number8{T <: FloatingPoint}(elems :: Array{T, 1})
    @assert length(elems) == 8
    Number8(elems, elems[1], elems[2], elems[3], elems[4], elems[5], elems[6], elems[7], elems[8])
end

function Number8{T <: FloatingPoint}(s0 :: T, s1 :: T, s2 :: T, s3 :: T, s4 :: T, s5 :: T, s6 :: T, s7 :: T)
    Number8([s0, s1, s2, s3, s4, s5, s6, s7])
end

exp(x :: Number8) = Number8(exp(x.array))
mod(n :: Number8, r :: Real) = Number8(mod(n.array, r))
sin(n :: Number8) = Number8(sin(n.array))

*(a :: Real, n :: Number8) = Number8(a * n.array)
*(n1 :: Number8, n2 :: Number8) = Number8(n1.array .* n2.array)
-(a :: Real, n :: Number8) = Number8(a - n.array)
-(n :: Number8) = Number8(-n.array)
/(n1 :: Number8, n2 :: Number8) = Number8(n1.array ./ n2.array)
/(n :: Number8, a :: Real) = Number8(n.array ./ a)

sum(n:: Number8) = sum(n.array)

fmod(a,b) = mod(a,b)

const PI = pi


text:
### A Pluto.jl notebook ###
# v0.14.4

using Markdown
using InteractiveUtils

# ╔═╡ f0a762cf-fd95-4350-8078-bab36104831c
using Pkg, DrWatson

# ╔═╡ 41d572f9-8e02-4751-a3a3-3a93f55eb418
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Symbolics, SymbolicUtils, ModelingToolkit, LaTeXStrings
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ b2c90da5-9f4b-4580-a4e4-7802fb8e827d
function lotka_volterra!(du1, u1, p1, t1)
  x, y = u1
  α, β, δ, γ = p1
  du1[1] = dx = α*x - β*x*y
  du1[2] = dy = -δ*y + γ*x*y
end

# ╔═╡ 8f0d514c-8d3e-4810-97ac-264317186c6c
begin
	@variables t du[1:2] u[1:2] p[1:4]
	lotka_volterra!(du, u, p, t)
	du = Num[
		p[1] * u[1] - (p[2] * u[1]) * u[2],
		-p[3] * u[2] + (p[4] * u[1]) * u[2]
	]
end;

# ╔═╡ f3513df2-c301-4581-ace2-b250eec62b72
D = Differential(t);

# ╔═╡ 7ff6f706-c7a7-4cbe-a8f9-12b1460e3bb3
eqs = @. D(u) ~ du

# ╔═╡ 1cdf03b5-fbac-44ff-b128-0771d9f9e1f7
build_function(eqs, u, p, t, target=Symbolics.CTarget());

# ╔═╡ 3a58c1e6-2bf7-494d-bae8-836d2482aba7
f = build_function(eqs, u, p, t, target=Symbolics.CTarget(),
	expression=Val{false});

# ╔═╡ 8c4074bc-8d24-4d4b-88f6-f753542abd05
du1 = rand(2)

# ╔═╡ dacbad6c-5fe5-4959-9ebe-b29a93f4237c
u1 = rand(2)

# ╔═╡ 183abe62-0e93-4ea5-9ff4-13afc6236e03
p1 = rand(4)

# ╔═╡ 79799876-9444-44d5-939b-4b6149add315
t1 = rand()

# ╔═╡ 99e6e2ca-0a7a-4cf9-b587-312486e1fc45
f(du1, u1, p1, t1)

# ╔═╡ 4e19a394-cd7b-41ee-91c5-357c7f0bf117
du1

# ╔═╡ aad3f488-7185-473c-948a-ff155bc81c2f
du3 = zeros(2)

# ╔═╡ f25c0f0f-6c38-4bf4-b18e-c1f6dd129d86
lotka_volterra!(du3, u1, p1, t1)

# ╔═╡ 3fde3bc3-94cd-4bb7-953d-1126d498d2ef
du3

# ╔═╡ c0b776a6-2e37-4b08-a9a7-e66576e10254
@test du1 == du3

# ╔═╡ Cell order:
# ╠═f0a762cf-fd95-4350-8078-bab36104831c
# ╠═41d572f9-8e02-4751-a3a3-3a93f55eb418
# ╠═b2c90da5-9f4b-4580-a4e4-7802fb8e827d
# ╠═8f0d514c-8d3e-4810-97ac-264317186c6c
# ╠═f3513df2-c301-4581-ace2-b250eec62b72
# ╠═7ff6f706-c7a7-4cbe-a8f9-12b1460e3bb3
# ╠═1cdf03b5-fbac-44ff-b128-0771d9f9e1f7
# ╠═3a58c1e6-2bf7-494d-bae8-836d2482aba7
# ╠═8c4074bc-8d24-4d4b-88f6-f753542abd05
# ╠═dacbad6c-5fe5-4959-9ebe-b29a93f4237c
# ╠═183abe62-0e93-4ea5-9ff4-13afc6236e03
# ╠═79799876-9444-44d5-939b-4b6149add315
# ╠═99e6e2ca-0a7a-4cf9-b587-312486e1fc45
# ╠═4e19a394-cd7b-41ee-91c5-357c7f0bf117
# ╠═aad3f488-7185-473c-948a-ff155bc81c2f
# ╠═f25c0f0f-6c38-4bf4-b18e-c1f6dd129d86
# ╠═3fde3bc3-94cd-4bb7-953d-1126d498d2ef
# ╠═c0b776a6-2e37-4b08-a9a7-e66576e10254


text:
# # Chapter 2

# ## Timing Julia
# `@time` is the simplest way to measure execution times
@time sqrt.(rand(1000));

# `@time` can also be used for any expression
s=0
@time for i=1:1000
   global s=s+sqrt(i)
end

# `@timev` provides additoinal information
@timev sqrt.(rand(1000));

# `@time` returns the value of the expression
sum(@time sqrt.(rand(1000)))

# `@elapsed` returns the elapsed execution time as a value ...
@elapsed sqrt.(rand(1000))

# ... and hence can be used in automated tests
using Test
@test @elapsed(sqrt.(rand(1000))) <= 10e-4

# ## Profiling

# Load the `Profile` stdlib module
using Profile

# We profile a simple function that creates a random matrix, and then
# computes the mean of the square of each row.
using Statistics
function randmsq()
    x = rand(10000, 1000)
    y = mean(x.^2, dims = 1)
    return y
end

# Run the function once before profiling, to compile it
randmsq()

# Execute the function while profiling it
@profile randmsq()

# `print` displays a tree view of the execution traces
Profile.print()

# Clear saved profiles
Profile.clear()

# Use the `init` function to configure the profiler
Profile.init(delay=.01)

# Run the profile using new setting
@profile randmsq()

# Profileview can be used to visualise the
# `] add ProfileView` to add the package

using Pkg
Pkg.add("ProfileView")

using ProfileView
ProfileView.view()
ProfileView.svgwrite("profile_results.svg")

# ## Statistically significant benchmarking
# `]add BenchmarkTools` to add the package
Pkg.add("BenchmarkTools")
using BenchmarkTools
@benchmark sqrt.(rand(1000))

@btime mean(rand(1000));


text:
using ChaosTools, DelayEmbeddings
using Test, OrdinaryDiffEq
using DelayEmbeddings: localextrema, exponential_decay_fit

println("\nTesting delay estimation...")

testval = (val, vmin, vmax) -> @test vmin ≤ val ≤ vmax

@testset "Estimate Delay" begin
    # test exponential decay fit
    x = 0:200
    y = @. exp(-x/5)
    @test exponential_decay_fit(x, y, :equal) ≈ 5
    @test exponential_decay_fit(x, y, :small) ≈ 5

    ds = Systems.henon()
    data = trajectory(ds,1000;dt=1)
    x = data[:,1]
    @test estimate_delay(x,"ac_zero", 0:10) ≤ 2
    @test estimate_delay(x,"ac_min", 0:10)  ≤ 2
    @test estimate_delay(x,"exp_extrema", 0:10)  ≤ 4
    @test estimate_delay(x,"mi_min", 0:10) == 10

    ds = Systems.roessler(ones(3))
    dt = 0.02
    data = trajectory(ds,500.0;dt=dt,diffeq...)
    x = data[:,1]
    @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:2:500)*dt ≤ 1.7
    @test 2.6 ≤ estimate_delay(x,"ac_min", 1:2:500)*dt  ≤ 3.4
    @test 1.0 ≤ estimate_delay(x,"mi_min", 1:2:500)*dt ≤ 1.6

    dt = 0.1
    data = trajectory(ds,2000.0;dt=dt,diffeq...)
    x = data[:,1]
    @test 1.3 ≤ estimate_delay(x,"ac_zero", 1:1:50)*dt ≤ 1.7
    @test 2.6 ≤ estimate_delay(x,"ac_min", 1:1:50)*dt  ≤ 3.4
    @test 1.15 ≤ estimate_delay(x,"mi_min", 1:1:50)*dt ≤ 1.6

    ds = Systems.lorenz()
    dt = 0.1
    data = trajectory(ds,5000;dt=dt,diffeq...)
    x = data[500:end,1]
    @test 0 < estimate_delay(x,"exp_extrema", 0:2:200)  < 200

    # Issue #18
    ds = Systems.gissinger(ones(3)) # 3D continuous chaotic system, also shown in orbit diagrams tutorial
    dt = 0.05
    data = trajectory(ds, 1000.0, dt = dt)
    s = data[:, 1]

    τ = estimate_delay(s, "mi_min", 0:1:400) # this was the erroring line
    @test τ > 0
end


text:
# References
#
# 18.8.30 numFmt (Number Format)
# lists a table with predefined numFmtId.
# In this case, a numFmtId value is written on the xf record,
# but no corresponding numFmt element is written.
# Some of these Ids can be interpreted differently,
# depending on the UI language of the implementing application.
#
# General formatCodes ranges from 0 to 81.
#
# 18.8.10 cellXfs (Cell Formats)

const STYLES_NAMESPACE_XPATH_ARG = [ "xpath" => "http://schemas.openxmlformats.org/spreadsheetml/2006/main" ]

_styles_root(wb::Workbook) = EzXML.root(wb.styles)

"""
Returns the xf XML node element for style `index`.
`index` is 0-based.
"""
function styles_cell_xf(wb::Workbook, index::Int) :: EzXML.Node
    xroot = _styles_root(wb)
    xf_elements = find(xroot, "/xpath:styleSheet/xpath:cellXfs/xpath:xf", STYLES_NAMESPACE_XPATH_ARG)
    return xf_elements[index+1]
end

styles_cell_xf(wb::Workbook, index::AbstractString) = styles_cell_xf(wb, parse(Int, index))
styles_cell_xf(ws::Worksheet, cell::Cell) = styles_cell_xf(ws.package.workbook, cell)

function styles_cell_xf(wb::Workbook, cell::Cell)
    @assert !isempty(cell.style) "Cell $(cell.ref.name) has empty style."
    styles_cell_xf(wb, cell.style)
end

"Queries numFmtId from cellXfs -> xf nodes."
function styles_cell_xf_numFmtId(wb::Workbook, index::Int) :: Int
    el = styles_cell_xf(wb, index)
    return parse(Int, el["numFmtId"])
end

"""
Queries numFmt formatCode field by numFmtId.
"""
function styles_numFmt_formatCode(wb::Workbook, numFmtId::AbstractString) :: String
    xroot = _styles_root(wb)
    elements_found = find(xroot, "/xpath:styleSheet/xpath:numFmts/xpath:numFmt[@numFmtId='$(numFmtId)']", STYLES_NAMESPACE_XPATH_ARG)
    @assert length(elements_found) == 1 "numFmtId $numFmtId not found."
    return elements_found[1]["formatCode"]
end

styles_numFmt_formatCode(wb::Workbook, numFmtId::Int) = styles_numFmt_formatCode(wb, string(numFmtId))

function styles_is_datetime(wb::Workbook, index::Int) :: Bool
    if !haskey(wb.buffer_styles_is_datetime, index)
        isdatetime = false

        numFmtId = styles_cell_xf_numFmtId(wb, index)

        if (14 <= numFmtId && numFmtId <= 22) || (45 <= numFmtId && numFmtId <= 47)
            isdatetime = true
        elseif numFmtId > 81
            code = styles_numFmt_formatCode(wb, numFmtId)
            if contains(code, "dd") || contains(code, "mm") || contains(code, "yy") || contains(code, "hh") || contains(code, "ss")
                isdatetime = true
            end
        end

        wb.buffer_styles_is_datetime[index] = isdatetime
    end

    return wb.buffer_styles_is_datetime[index]
end

function styles_is_datetime(wb::Workbook, index::AbstractString)
    @assert !isempty(index)
    styles_is_datetime(wb, parse(Int, index))
end

styles_is_datetime(ws::Worksheet, index) = styles_is_datetime(ws.package.workbook, index)

function styles_is_float(wb::Workbook, index::Int) :: Bool
    if !haskey(wb.buffer_styles_is_float, index)
        isfloat = false

        numFmtId = styles_cell_xf_numFmtId(wb, index)

        if numFmtId == 2 || numFmtId == 4 || (7 <= numFmtId && numFmtId <= 11) || numFmtId == 39 || numFmtId == 40 || numFmtId == 44 || numFmtId == 48
            isfloat = true
        elseif numFmtId > 81
            code = styles_numFmt_formatCode(wb, numFmtId)
            if contains(code, ".0")
                isfloat = true
            end
        end

        wb.buffer_styles_is_float[index] = isfloat
    end

    return wb.buffer_styles_is_float[index]
end

function styles_is_float(wb::Workbook, index::AbstractString)
    @assert !isempty(index)
    styles_is_float(wb, parse(Int, index))
end

styles_is_float(ws::Worksheet, index) = styles_is_float(ws.package.workbook, index)


text:
"""
    Region4

    Returns the vapour-liquid phase boundary (saturation line).
    Valid from triple point to critical point
    273.15K ≤ T ≤ 647.096K
    InputType is either :T or :P to indicate that InputValue is temperature [K]
    or pressure [MPa]. The complimentary value is returned.
"""
function Region4 end

const Region4_tup = (;n = [0.116_705_214_527_67E4,
    -0.724_213_167_032_06E6,
    -0.170_738_469_400_92E2,
    0.120_208_247_024_70E5,
    -0.323_255_503_223_33E7,
    0.149_151_086_135_30E2,
    -0.482_326_573_615_91E4,
    0.405_113_405_420_57E6,
    -0.238_555_575_678_49,
    0.650_175_348_447_98E3]
)

function pressure_impl(mt::SingleSatT,model::IF97Region{:r4},t)
    n = Region4_tup.n
    Θ = t + n[9]/(t - n[10])
    A =      Θ^2 + n[1]*Θ + n[2]
    B = n[3]*Θ^2 + n[4]*Θ + n[5]
    C = n[6]*Θ^2 + n[7]*Θ + n[8]
    P = (2C / (-B + √(B^2 - 4*A*C)))^4
    return convert_unit(u"MPa",u"Pa",P)
end

function temperature_impl(mt::SingleSatP,model::IF97Region{:r4},p)
    n = Region4_tup.n
    P = convert_unit(u"Pa",u"MPa",p)
    β = P^0.25
    E =      β^2 + n[3]*β + n[6]
    F = n[1]*β^2 + n[4]*β + n[7]
    G = n[2]*β^2 + n[5]*β + n[8]
    D = 2G / (-F - √(F^2 - 4*E*G))
    T = (n[10]+D-√((n[10]+D)^2 - 4(n[9]+n[10]*D)))/2.0
    return T
end

pressure_impl(mt::SingleΦT,model::IF97Region{:r4},t) = pressure_impl(QuickStates.sat_p(),model,t)
temperature_impl(mt::SingleΦP,model::IF97Region{:r4},p) = temperature_impl(QuickStates.sat_t(),model,p)


text:
# Submodule for basic wavelets operations

###################################################
# Function that remove empty trailing dimensions
function trim_dim(x)
	if size(x,3)==1
		x = x[:,:,1];
		if size(x,2)==1
			x = x[:,1];
		end
	end
	return x;
end # trim_dim

###################################################
# Function that make a subsampling of an array
# along dimension d
"""
...
subsampling(x,d=1,p=2) subsampling of an array x, along dimension d, by a
factor of p, for d≤3
...
"""
function subsampling(x,d=1,p=2)
	if d==1
	        y = x[1:p:end,:,:];
	elseif d==2
	        y = x[:,1:p:end,:];
	elseif d==3
	        y = x[:,:,1:p:end];
	else
	        error("Not implemented");
	end
	y = trim_dim(y);
	return y;
end # subsampling

###################################################
# Function that make a upsampling of an array
# along dimension d
"""
...
upsampling(x,d=1,p=2) upsampling of an array x, along dimension d, by a
factor of p, for d≤3
...
"""
function upsampling(x,d=1,p=2)
	if d==1
	        y = zeros(p*size(x,1),size(x,2),size(x,3));
	        y[1:p:end,:,:] = x;
	elseif d==2
	        y = zeros(size(x,1),p*size(x,2),size(x,3));
	        y[:,1:p:end,:] = x;
	elseif d==3
	        y = zeros(size(x,1),size(x,2),p*size(x,3));
	        y[:,:,1:p:end] = x;
	else
	        error("Not implemented");
	end
	y = trim_dim(y);

	return y;
end # upsampling

###################################################
# Function to compute the spatial domain circular
# convolution of two signals
"""
...
cconv(x,h,d=1) circular convolution along dimension d, of
signals x and h
...
"""
function cconvol(x,h,d=1)
	#   If p=length(h) is odd then h((p+1)/2) corresponds to the zero position.
	#   If p is even, h(p/2) corresponds to the zero position.
	if d==2
	    return ( cconvol( x',h, 1 )' );
	elseif d==3
		# permute do not exists in python
	    return permute( cconvol( permute(x,[3, 2, 1]),h), [3 2 1]);
	end
	p = length(h);
	if mod(p,2)==0
	    pc = p/2;
	else
	    pc = (p+1)/2;
	end
	y = zeros(size(x));
	for i=1:length(h)
	    y = y + h[i]*circshift(x,i-pc);
	end
	return y;
end # cconvol

#####################################
# Function that flips a vector
"""
...
reverse(x::Array) flips a vector
...
"""
function reverse(x)
	return x[length(x):-1:1];
end  # function x = reverse(x)

#####################################
# Function that computes the mirror
# filter
"""
...
mirror(h::Array) computes the mirror filter of h[n], namely
g[n] = (-1)^(1-n) h[1-n]
...
"""
function mirror(h)
	cat(1, 0, h[length(h):-1:2]) .* ( (-1).^(1:length(h)) )
end #mirror

#####################################
# Function that clamps a value
"""
...
clamp(x::Array, a=0,b=1) clamps a value in [a,b](=[0,1] by default)
...
"""
function clamp(x,a=0,b=1)
	y = max(x,a);
	y = min(y,b);
	return y;

end # clamp

#########################################
# Function to subselect dimensions in an
# array
function subselectdim(f,sel,d)
	g = [];
	if d==1
	    g = f[sel,:,:];
	elseif d==2
	    g = f[:,sel,:];
	elseif d==3
	    g = f[:,:,sel];
	end
	return trim_dim(g);
end

#########################################
# Function to assign value to elements
# of an array
function subassign(f,sel,g)
	d = ndims(f);
	if d==1
	    f[sel] = g;
	elseif d==2
	    f[sel,sel] = g;
	elseif d==3
	    f[sel,sel,sel] = g;
	end
	return f;
end # subassign

#########################################
# Function to select elements of an array
function subselect(f,sel)
	d = ndims(f)
    if d==1
        return f[sel];
    elseif d==2
        return f[sel,sel];
    elseif d==3
        return f[sel,sel,sel];
	end
	return [];
end

##############################################################
# Function to perform wavelet orthognal transform of a signal
# f0, with Jmin minimum scale, dir dirección, and h a filter
"""
...
perform_wavortho_transform(f0:::Array, Jmin:::int,dir,h) computes
the wavelet transform of a signal f0, with Jmin as the minimum
scale in the direction dir∈{-1,1}
...
"""
function perform_wavortho_transf(f0,Jmin,dir,h)

	n = size(f0,1);
	Jmax = round(Int64,log2(n))-1;

	g = [0; h[length(h):-1:2]] .* (-1).^(1:length(h));

	f = copy(f0)
	if dir==1
	    ### FORWARD ###
	    for j=Jmax:-1:Jmin
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            a = cat(d, subsampling(cconvol(a,h,d),d), subsampling(cconvol(a,g,d),d) );
	        end
	        f = subassign(f,sel,a);
	    end
	else
	    ### FORWARD ###
	    for j=Jmin:Jmax
	        sel = 1:2^(j+1);
	        a = subselect(f,sel);
	        for d=1:ndims(f)
	            w = subselectdim(a,2^j+1:2^(j+1),d);
	            a = subselectdim(a,1:2^j,d);
	            a = cconvol(upsampling(a,d),reverse(h),d) + cconvol(upsampling(w,d),reverse(g),d);
	        end
	        f = subassign(f,sel,a);
	    end
	end
	return f;
end # perform_wavortho_transf

#####################################################
# Function that computes the signal to noise radio
"""
...
snr(x::Array,y:::Array) computes the siganl to noise radio
of a original clean signal x and a denoised signal y
...
"""
function snr(x,y)
	return 20*log10(norm(x[:])/norm(x[:]-y[:]));
end # snr


text:
function IBClient()
    return IBClient(JLClient())
end


text:
const operatorList = Dict{String, Function}( "IS" => Base.:(==),
    "IS NOT" => Base.:(!=), "IN" => Base.:∈, "NOT IN" => Base.:∉,
    ">" => Base.:>, ">=" => Base.:(>=), "<" => Base.:<, "<=" => Base.:(<=) )

function readTransitions( mpSim::MPsim, xf::XF )

    if !XLSX.hassheet( xf, "Transition types" )
        @warn "Configuration file is missing a 'Base Nodes' sheet, cannot set preferred transition order."
    else
        ws = xf[ "Transition types" ]
        nTrans = ws[ "B4" ] isa Integer ? ws[ "B4" ] : 0

        # Set the preferred transition order.
        transitions = ws[ CR( 7, 1, 6 + nTrans, 1 ) ]
        order = ws[ CR( 7, 2, 6 + nTrans, 2 ) ]
        isValid = .!isa.( transitions, Missing ) .& isa.( order, Integer )
        transitions = string.( transitions[ isValid ] )
        order = Vector{Int}( order[ isValid ] )
        setSimulationTransitionTypeOrder!( mpSim, transitions, order )
    end  # if !XLSX.hassheet( xf, "Transition types" )

    readRecruitment( mpSim, xf )
    readTransitionsThrough( mpSim, xf )
    readTransitionsOut( mpSim, xf )

end  # readTransitions( mpSim, xf )

function readTransitionsThrough( mpSim::MPsim, xf::XF )

    if !XLSX.hassheet( xf, "THROUGH Transitions" )
        @warn "Configuration file is missing a 'THROUGH Transitions' sheet, cannot read transitions within the organisation."
        return
    end  # if !XLSX.hassheet( xf, "THROUGH Transitions" )

    ws = xf[ "THROUGH Transitions" ]
    nTrans = ws[ "B4" ] isa Integer ? max( ws[ "B4" ], 0 ) : 0

    for sLine in (1:nTrans) * 4 .+ 3
        sLine1 = sLine + 1
        sLine2 = sLine + 2

        # Read transition base information.
        name = ws[ sLine, 1 ]
        sourceNode, targetNode = ws[ sLine, 2 ], ws[ sLine, 3 ]

        if ( name isa Missing ) || ( sourceNode isa Missing ) ||
            ( targetNode isa Missing )
            continue
        end  # if ( name isa Missing ) || ...

        transition = Transition( name, sourceNode, targetNode )
        freq, offset = ws[ sLine, 4 ], ws[ sLine, 5 ]
        setTransitionSchedule!( transition, freq, offset )

        # Read time conditions of transition.
        sCol = 7
        nTimeConds = ws[ sLine, sCol ]

        for ii in (1:nTimeConds) .+ sCol
            newCond, isCondOkay = processCondition( ws[ sLine, ii ],
                ws[ sLine1, ii ], ws[ sLine2, ii ] )

            if isCondOkay
                addTransitionCondition!( transition, newCond )
            end  # if isCondOkay
        end  # for ii in (1:nTimeConds) .+ sCol

        # Read other conditions of transitions.
        sCol += nTimeConds + 2
        nOtherConds = ws[ sLine, sCol ]

        for ii in (1:nOtherConds) .+ sCol
            newCond, isCondOkay = processCondition( ws[ sLine, ii ],
                ws[ sLine1, ii ], ws[ sLine2, ii ] )

            if isCondOkay
                addTransitionCondition!( transition, newCond )
            end  # if isCondOkay
        end  # for ii in (1:nOtherConds) .+ sCol

        # Set transition flux limits.
        sCol += nOtherConds + 2
        minFlux = ws[ sLine, sCol ] isa Integer ? ws[ sLine, sCol ] : 0
        maxFlux = ws[ sLine, sCol + 1 ] isa Integer ? ws[ sLine, sCol + 1 ] : 0
        setTransitionFluxLimits!( transition, minFlux, maxFlux )
        setTransitionHasPriority!( transition,
            strip( string( ws[ sLine, sCol + 2 ] ) ) == "NO" )

        # Set transition attempts and probabilities.
        nAttempts = ws[ sLine, sCol + 3 ] isa Integer ? ws[ sLine, sCol + 3 ] :
            -1
        setTransitionMaxAttempts!( transition, nAttempts )

        sCol += 4
        nProbs = ws[ sLine, sCol ]
        transitionProbs = Float64.( ws[ CR( sLine, sCol + 1, sLine,
            sCol + nProbs ) ] )[ : ]
        setTransitionProbabilities!( transition, transitionProbs )

        # Set extra transition changes.
        sCol += nProbs + 2
        nChanges = ws[ sLine, sCol ]
        attributes = Vector{String}()
        values = Vector{String}()

        for ii in (1:nChanges) .+ sCol
            if !any( isa.( [ ws[ sLine, ii ], ws[ sLine + 1, ii ] ], Missing ) )
                push!( attributes, strip( string( ws[ sLine, ii ] ) ) )
                push!( values, strip( string( ws[ sLine1, ii ] ) ) )
            end  # if !any( isa.( ... ) )
        end  # for ii in 1:nChanges

        setTransitionAttributeChanges!( transition, attributes, values )
        addSimulationTransition!( mpSim, transition )
    end  # for sLine (1:nTrans) * 4 .+ 3

end  # readTransitionsThrough( mpSim, xf )

function readTransitionsOut( mpSim::MPsim, xf::XF )

    if !XLSX.hassheet( xf, "OUT Transitions" )
        @warn "Configuration file is missing a 'OUT Transitions' sheet, cannot read transitions within the organisation."
        return
    end  # if !XLSX.hassheet( xf, "OUT Transitions" )

    ws = xf[ "OUT Transitions" ]
    nTrans = ws[ "B4" ] isa Integer ? max( ws[ "B4" ], 0 ) : 0

    for sLine in (1:nTrans) * 4 .+ 3
        sLine1 = sLine + 1
        sLine2 = sLine + 2

        # Read transition base information.
        name = ws[ sLine, 1 ]
        sourceNode = ws[ sLine, 2 ]

        if ( name isa Missing ) || ( sourceNode isa Missing )
            continue
        end  # if ( name isa Missing ) || ...

        transition = Transition( name, sourceNode )
        freq, offset = ws[ sLine, 3 ], ws[ sLine, 4 ]
        setTransitionSchedule!( transition, freq, offset )

        # Read time conditions of transition.
        sCol = 6
        nTimeConds = ws[ sLine, sCol ]

        for ii in (1:nTimeConds) .+ sCol
            newCond, isCondOkay = processCondition( ws[ sLine, ii ],
                ws[ sLine1, ii ], ws[ sLine2, ii ] )

            if isCondOkay
                addTransitionCondition!( transition, newCond )
            end  # if isCondOkay
        end  # for ii in (1:nTimeConds) .+ sCol

        # Read other conditions of transitions.
        sCol += nTimeConds + 2
        nOtherConds = ws[ sLine, sCol ]

        for ii in (1:nOtherConds) .+ sCol
            newCond, isCondOkay = processCondition( ws[ sLine, ii ],
                ws[ sLine1, ii ], ws[ sLine2, ii ] )

            if isCondOkay
                addTransitionCondition!( transition, newCond )
            end  # if isCondOkay
        end  # for ii in (1:nOtherConds) .+ sCol

        # Set transition flux limits.
        sCol += nOtherConds + 2
        minFlux = ws[ sLine, sCol ] isa Integer ? ws[ sLine, sCol ] : 0
        maxFlux = ws[ sLine, sCol + 1 ] isa Integer ? ws[ sLine, sCol + 1 ] : 0
        setTransitionFluxLimits!( transition, minFlux, maxFlux )

        # Set transition attempts and probabilities.
        nAttempts = ws[ sLine, sCol + 2 ] isa Integer ? ws[ sLine, sCol + 2 ] :
            -1
        setTransitionMaxAttempts!( transition, nAttempts )

        sCol += 3
        nProbs = ws[ sLine, sCol ]
        transitionProbs = Float64.( ws[ CR( sLine, sCol + 1, sLine,
            sCol + nProbs ) ] )[ : ]
        setTransitionProbabilities!( transition, transitionProbs )
        addSimulationTransition!( mpSim, transition )
    end  # for sLine (1:nTrans) * 4 .+ 3

end  # readTransitionsOut( mpSim, xf )

function processCondition( attribute::Union{String, Missing},
    operator::Union{String, Missing}, value::Union{Real, String, Missing} )

    # Check if any of the fields are missing.
    if any( isa.( [ attribute, operator, value ], Missing ) )
        return MPcondition( "condition", Base.:!=, "okay" ), false
    end  # if any( isa.( [ attribute, operator, value ], Missing ) )

    attribute = string( strip( attribute ) )
    operator = uppercase( string( strip( operator ) ) )

    # Check if the operator is valid.
    if !haskey( operatorList, operator )
        return MPcondition( "condition", Base.:!=, "okay" ), false
    end  # if !haskey( operatorList, operator )

    operator = operatorList[ operator ]

    # Time conditions must have a number value.
    if attribute ∈ MP.timeAttributes
        if ( value isa String ) || ( operator ∈ [ ∈ , ∉ ] )
            return MPcondition( "condition", Base.:!=, "okay" ), false
        end  # if ( value isa String ) || ...

        return MPcondition( attribute, operator, value * 12.0 ), true
    end  # if attribute ∈ MP.timeAttributes

    value = strip( string( value ) )

    if operator ∈ [ "IN", "NOT IN" ]
        value = Vector{String}( strip.( split( value, "," ) ) )
    end  # if operator ∈ [ "IN", "NOT IN" ]

    return MPcondition( attribute, operator, value ), true

end  # processCondition( attribute, operator, value )

function readTransitions( mpSim::MPsim, configDB::SQLite.DB,
    configName::String )

    transitionPars = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Transition'" ) ) )

    if !isempty( transitionPars )
        transitions = map( eachindex( transitionPars[ :, :parName ] ) ) do ii
            pars = string.( split( transitionPars[ ii, :parValue ], ";" ) )
            transition = uppercase( pars[ 2 ] ) == "OUT" ?
                Transition( transitionPars[ ii, :parName ], pars[ 1 ] ) :
                Transition( transitionPars[ ii, :parName ], pars[ 1 ],
                    pars[ 2 ] )

            if ( tryparse( Float64, pars[ 3 ] ) isa Float64 ) &&
                ( tryparse( Float64, pars[ 4 ] ) isa Float64 )
                setTransitionSchedule!( transition, parse( Float64, pars[ 3 ] ),
                    parse( Float64, pars[ 4 ] ) )
            end  # if ( tryparse( Float64, pars[ 3 ] ) isa Float64 ) && ...

            if tryparse( Int, pars[ 5 ] ) isa Int
                setTransitionMaxAttempts!( transition, parse( Int, pars[ 5 ] ) )
            end  # if tryparse( Int, pars[ 5 ] ) isa Int

            if ( tryparse( Int, pars[ 6 ] ) isa Int ) &&
                ( tryparse( Int, pars[ 7 ] ) isa Int )
                setTransitionFluxLimits!( transition, parse( Int, pars[ 6 ] ),
                    parse( Int, pars[ 7 ] ) )
            end  # if ( tryparse( Int, pars[ 6 ] ) isa Int ) && ...

            if tryparse( Bool, pars[ 8 ] ) isa Bool
                setTransitionHasPriority!( transition,
                    parse( Bool, pars[ 8 ] ) )
            end  # if tryparse( Bool, pars[ 8 ] ) isa Bool

            if length( pars[ 9 ] ) > 2
                conditionPars = split.( split( pars[ 9 ][ 2:(end-1) ], "," ),
                    ":" )
                conditions = map( conditionPars ) do condPars
                    if length( condPars ) == 3
                        attrVal = tryparse( Float64, condPars[ 3 ] )
                        attrVal = attrVal isa Nothing ?
                            string( condPars[ 3 ] ) : attrVal / 12.0
                            # To ensure proper reading of conditions.
                    else
                        attrVal = join( condPars[ 3:end ], "," )
                    end  # if length( condPars ) == 3

                    condition, isOkay = processCondition(
                        string( condPars[ 1 ] ), string( condPars[ 2 ] ),
                        attrVal )
                    return condition
                end  # map( conditionPars ) do condPars

                setTransitionConditions!( transition, conditions )
            end  # if length( pars[ 9 ] ) > 2

            if length( pars[ 10 ] ) > 2
                changes = split.( split( pars[ 10 ][ 2:(end-1) ], "," ), ":" )
                changesDict = Dict{String, String}()

                for change in changes
                    changesDict[ change[ 1 ] ] = change[ 2 ]
                end  # for nodeReq in nodeReqs

                setTransitionAttributeChanges!( transition, changesDict )
            end  # if length( pars[ 10 ] ) > 2

            if length( pars[ 11 ] ) > 2
                setTransitionProbabilities!( transition,
                    parse.( Float64, split( pars[ 11 ][ 2:(end-1) ] ) ) )
            end  # if length( pars[ 11 ] ) > 2

            return transition
        end  # map( eachindex( transitionPars[ :parName ] ) ) do ii

        setSimulationTransitions!( mpSim, transitions )
    end  # if !isempty( transitionPars )

    transOrder = DataFrame( DBInterface.execute( configDB,
        string( "SELECT * FROM `", configName,
        "` WHERE parType IS 'Transition Order'" ) ) )

    if !isempty( transOrder ) && ( length( transOrder[ 1, :parValue ] ) > 2 )
        order = split.( split( transOrder[ 1, :parValue ], ";" ), ":" )
        orderDict = Dict{String, Int}()

        for orderPair in order
            orderDict[ orderPair[ 1 ] ] = parse( Int, orderPair[ 2 ] )
        end  # for orderPair in order

        setSimulationTransitionTypeOrder!( mpSim, orderDict )
    end  # if !isempty( nodeOrder ) && ...

end  # readTransitions( mpSim, configDB, configName )


text:
using Yao
using BitBasis
using QuDiffEq
using Test, LinearAlgebra

function crot(n_reg::Int, C_value::Real)
    n_rot = n_reg + 1
    rot   = chain(n_rot)
    θ = zeros(1<<n_reg - 1)
    for i = 1:(1<<n_reg - 1)
        c_bit = Vector(2:n_rot)
        λ = 0.0
        for j = 1:n_reg
            if (readbit(i,j) == 0)
                c_bit[j] = -c_bit[j]
            end
        end
        λ = i/(1<<n_reg)
        # println("\nλ($i) = $λ")
        # println("c_bit($i) = $c_bit\n")
        sin_value = C_value / λ
        if (sin_value) > 1
            return println("C_value = $C_value, λ = $λ, sinθ = $sin_value > 1, please lower C_value.\n")
        end
        θ[(i)] = 2.0*asin(C_value / λ)
        push!(rot, control(c_bit, 1=>Ry(θ[(i)])))
    end
    rot
end

@testset "HHLCRot" begin
    hr = HHLCRot{4}([4,3,2], 1, 0.01)
    reg = rand_state(4)
    @test reg |> copy |> hr |> isnormalized

    hr2 = crot(3, 0.01)
    reg1 = reg |> copy |> hr
    reg2 = reg |> copy |> hr2
    @test fidelity(reg1, reg2)[] ≈ 1
end

"""
    hhl_problem(nbit::Int) -> Tuple

Returns (A, b), where
    * `A` is positive definite, hermitian, with its maximum eigenvalue λ_max < 1.
    * `b` is normalized.
"""
function hhl_problem(nbit::Int)
    siz = 1<<nbit
    base_space = qr(randn(ComplexF64, siz, siz)).Q
    phases = rand(siz)
    signs = Diagonal(phases)
    A = base_space*signs*base_space'

    # reinforce hermitian, see issue: https://github.com/JuliaLang/julia/issues/28885
    A = (A+A')/2
    b = normalize(rand(ComplexF64, siz))
    A, b
end

@testset "HHLtest" begin
    # Set up initial conditions.
    ## A: Matrix in linear equation A|x> = |b>.
    ## signs: Diagonal Matrix of eigen values of A.
    ## base_space: the eigen space of A.
    ## x: |x>.
    using Random
    Random.seed!(2)
    N = 3
    A, b = hhl_problem(N)
    x = A^(-1)*b # base_i = base_space[:,i] ϕ1 = (A*base_i./base_i)[1]

    ## n_b  : number of bits for |b>.
    ## n_reg: number of PE register.
    ## n_all: number of all bits.
    n_reg = 12

    ## C_value: value of constant C in control rotation.
    ## It should be samller than the minimum eigen value of A.
    C_value = minimum(eigvals(A) .|> abs)*0.25
    #C_value = 1.0/(1<<n_reg) * 0.9
    res = hhlsolve(A, b, n_reg, C_value)

    # Test whether HHL circuit returns correct coefficient of |1>|00>|u>.
    @test isapprox.(x, res, atol=0.5) |> all
end


text:
#MAIN ENTRY POINT
#final columns for the decompressed file
#TODO: Fix the assertion checks

#merges crsp and compustat
# note refresh merge will be automatically run if any of the other two are run
#pre-op allows for a function to be applied before univ is written
function mergecrsp(crsp::Union{DataFrame,Nothing}=nothing,
  comp::Union{DataFrame,Nothing}=nothing;
  crsptype::Symbol = :monthly,
  refreshcomp::Bool = true,
  refreshmerge::Bool = true,
  datapath::String=DATA_PATH,
  prefix::String = "",
  univname::String = "$(prefix)$(crsptype==:monthly ? "CRSP-M" : "CRSP-D")_$(COMP_NAME)_$UNIV_NAME",
  inbinstream::Function = IN_BIN_STREAM,
  outbinstream::Function = OUT_BIN_STREAM,
  binextension::String = BIN_EXTENSION,
  validatemerged::Bool = VALIDATE_MERGED,
  retainedcolsuniv::Vector{Symbol} = [propertynames(comp); propertynames(crsp); :compid],
  cachedata::Bool = true
  )::DataFrame

  local univ::DataFrame

  #as of now, refreshmerged must be true if any of the other flags are true
  refreshmerge = refreshmerge || refreshcrsp || refreshcomp

  #makes a serialization object so we don't have to keep parsing the CSV
  if refreshmerge
    univ = mergecrspcomp!(crsp, comp, validatemerged=validatemerged)
    select!(univ, retainedcolsuniv)

    if TEST_OUTPUT
      d = Date(2015,12,31)
      CSV.write("output\\$(univname)_$d.csv", univ[(cd->cd==d).(univ.date), :])
    end

    if cachedata
      outbinstream("$datapath\\$univname.$binextension", univ)
    end

  else
    univ = inbinstream("$datapath\\$univname.$binextension")
  end

  @info "CRSP and Compustat merged in file $(univname).$binextension."

  return univ
end

#this is core code to prepare the merge
function matchcompid!(crsp::DataFrame, comp::DataFrame;
    mergetype::Symbol = :full,
    Fcrspdate=error("Fcrspdate is required"))

  (issorted(crsp, [:permno, Fcrspdate])) || error(
    "crsp comp merge must be sorted by :permno, $Fcrspdate")

  scrsps::GroupedDataFrame = groupby(crsp, :permno)
  scomps::Dict = Dict([scomp.lpermno[1]=>scomp for scomp ∈ groupby(comp, :lpermno)])

  #join(crsp, comp, kind)

@mpar for i ∈ 1:length(scrsps)
    local scrsp::SubDataFrame = scrsps[i]
    local permno::Int = scrsp.permno[1]

    if haskey(scomps, permno)
      local scomp = scomps[permno]
      for rscomp ∈ eachrow(scomp)
        begindate::Date =  rscomp[:begindate]
        enddate::Date = rscomp[:enddate]
        sscrsp::SubDataFrame = view(scrsp, (d->d>begindate && d≤enddate).(scrsp[!, Fcrspdate]),:)
        if size(sscrsp,1) ≥ 1
          (sum((!ismissing).(sscrsp.compid)) > 0) && error("multiple valid compids for one permno")
          sscrsp.compid .= rscomp[:compid]
        end
      end
    end
  end

end

function mergecrspcomp!(crsp::DataFrame, comp::DataFrame;
    validatemerged::Bool=true,
    joinkind::Symbol = :inner,
    Fcrspdate::Symbol = :date,
    allowdatecollisions::Bool=false)::DataFrame

  local univ::DataFrame
  local scomp::SubDataFrame

  local Ncrsp::Int = size(crsp,1)
  #create a record id to form the dictionary
  comp.compid = Vector{MInt}(1:(size(comp,1)))
  comp.validinterval = ((begindate,enddate)->
    (begindate+Day(1)):Day(1):(enddate)).(comp.begindate, comp.enddate)

  crsp.compid = Vector{MInt}(undef, Ncrsp)

  #make the necessary allocations ahead of time
  print("matching...time to match: ")
  @time matchcompid!(crsp,comp, Fcrspdate=Fcrspdate)

  #println("sum crsp compid: $(sum(skipmissing(crsp.compid)))")
  println("Matching complete: $(100*(Ncrsp - sum((ismissing).(crsp.compid)))/Ncrsp)% matched.")
  #println("crsp records: $(Ncrsp)")
  if joinkind==:inner
    delete!(crsp, (ismissing).(crsp.compid))
    univ = innerjoin(crsp, comp, on=:compid)
  elseif joinkind == :left
    univ = leftjoin(crsp, comp, on=:compid)
  else
    error("Unrecognized joinkind: $joinkind")
  end

  validatemerged && testmerged(univ, Fdate=Fcrspdate, allowdatecollisions=allowdatecollisions)
  println("Total records: $(size(univ))")

  if TEST_OUTPUT
    d = Date(2015,12,31)
    CSV.write("output\\univ_mid_$d.csv", univ[(cd->cd==d).(univ[!, Fcrspdate]), :])
  end

  return univ
end

#do a bunch of integrity checks on the merged dataset
function testmerged(univ::DataFrame;
  months2stale::Int = MONTHS_2_STALE_LONG,
  Fdate::Symbol = error("Fdate is required"),
  allowdatecollisions::Bool = false)

  @mpar for r ∈ eachrow(univ)
    ismissing(r.gvkey) && continue
    #check date is valid
    d::Date = r[Fdate]
    (d ≥ r.linkeffdate) || error("failed check: (d ≥ r.linkeffdate)")
    (d ≤ r.linkenddate) || error("failed check: (d ≤ r.linkenddate)")
    (d ∈ r.validinterval) || error("failed check: (d ∈ r.validinterval)")
    (d > r.adate) || error("failed check: (d > r.adate)")
    (d ≤ r.Nadate) || error("failed check: (d ≤ r.Nadate)")
    period2stale = Month(months2stale)

    (!(d ≤ lastdayofmonth(firstdayofmonth(r.adate) + period2stale))) && error(
      "lastdayofmonth(firstdayofmonth(r.adate) + period2stale)))
      $r")

    (r.lpermno == r.permno) || error("failed check: (r.lpermno == r.permno)")
  end

  univ = sort(univ, [:permno, :adate, Fdate])

  sunivs::GroupedDataFrame = groupby(univ, :permno)
  @mpar for i ∈ 1:length(sunivs)
    local suniv::SubDataFrame = sunivs[i]
    maxfiscal::Int= 1900
    maxadate::Date = Date(1900,1,1)

    for r ∈ eachrow(suniv)
      ismissing(r.gvkey) && continue
      (!allowdatecollisions) && ((r[Fdate] > maxadate) || error("failed check: r[Fdate] > maxadate
        permno=$(r[:permno]), r[Fdate]=$(r[Fdate]), maxadate=$maxadate\n\n$r"))

      (allowdatecollisions) && ((r[Fdate] ≥ maxadate) || error("failed check: r[Fdate] ≥ maxadate
        permno=$(r[:permno]), r[Fdate]=$(r[Fdate]), maxadate=$maxadate"))

      maxadate = r[Fdate]

      (maxfiscal ≤ r.fyear) ||  error("failed check: maxfiscal ≤ r.fyear")
      (maxfiscal ≠ r.fyear) && (maxfiscal=r.fyear)
    end

  end
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# using Revise
using mujoco_sim
import Distributions: Uniform
using MatrixEquations

keypath = joinpath(@__DIR__, "../key/mjkey.txt")
activate_mujoco(keypath)

function lqr(A, B, Q, R)
    G = B*inv(R)*B'
    P = arec(A, G, Q)[1]  #continuous ricatti eqn solver
    K = inv(R)*B'*P
    return K
end

T = 100
mk = 1.5
mp = 0.5
g = 9.81
lp = 0.3
noise = 0.5

a = g/(lp*(4.0/3 - mp/(mp+mk)))
A =[[0. 1. 0. 0.];
    [0. 0. a  0.];
    [0. 0. 0. 1.];
    [0. 0. a  0.]]
b = -1/(lp*(4.0/3 - mp/(mp+mk)))
B = [[0.]; [1.0/(mp+mk)]; [0.]; [b]]

Q = [[5. 0. 0. 0.];
     [0. 5. 0. 0.];
     [0. 0. 5. 0.];
     [0. 0. 0. 5.]]
R = 1.

K = lqr(A, B, Q, R)

function apply_ctrl(K,x)
    u = -K*x
    if u > 1.0
        return [1.0]
    elseif u < -1.0
        return [-1.0]
    else
        return [u]
    end
end

function ctrler!(s)
    x = [s.d.qpos[1], s.d.qvel[1],s.d.qpos[2], s.d.qvel[2]]
    u = apply_ctrl(K,x)

    #randomize init pole deviation to make things interesting
    if s.d.time < 0.25
        u = [rand(Uniform(-noise, noise))]
    end
    setaction!(s, u)
    # println("state: ",s.d.qpos[2], " action: ",s.d.ctrl)
end

m = get_model("models/cartpole.xml")
d = get_data(m)
sim = MJSim(m,d)
simulate( sim, controller = ctrler!; mode="active")


text:
# Parameters:
#	N::Int
#	- the problem size
#	SNF::Array{Foat64, 1}
#	- SNF[i] is the value associated with the Permutation that permutation_index() maps to i
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1 from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2 from yor()
# Return Values:
#	FFT::Array{Float64, 2}
#	- FFT is the Fast Fourier Transform of SNF
function sn_fft(N::Int, SNF::Array{Float64, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1})
	np = nprocs()
	if np == 1 || N < 10
		C = Counter(1)
		return compute_fft(N, SNF, YOR, PT, C)
	else #This is a modified version of Julia's pmap()
		sFFT = Array(Array{Array{Float64, 2}, 1}, N)
		RR_YOR = Array(RemoteRef, np)
		RR_PT = Array(RemoteRef, np)
		for p = 1:np
			if p != myid()
				RR_YOR[p] = RemoteRef(p)
				put!(RR_YOR[p], YOR)
				RR_PT[p] = RemoteRef(p)
				put!(RR_PT[p], PT)
			end
		end
		BS = factorial(N - 1)
		i = 1
		nextidx() = (idx = i; i += 1; idx)
		@sync begin
			for p = 1:np
				if p != myid()
					@async begin
						while true
							idx = nextidx()
							if idx > N
								break
							end
							lb = (idx - 1) * BS + 1
							ub = lb + BS - 1
							SNF_subgroup = SNF[lb:ub]
							sFFT[idx] = remotecall_fetch(p, compute_fft_remote, N - 1, SNF_subgroup, RR_YOR[p], RR_PT[p], Counter(1))
						end
					end
				end
			end
		end
		YORn = YOR[N]
		NP = length(YORn)
		PTn = PT[N]
		RR_sFFT = Array(RemoteRef, np)
		for p = 1:np
			if p != myid()
				RR_sFFT[p] = RemoteRef(p)
				put!(RR_sFFT[p], sFFT)
			end
		end
		FFT = Array(Array{Float64, 2}, NP)
		i = 1
		nextidx() = (idx = i; i += 1; idx)
		@sync begin
			for p = 1:np
				if p != myid()
					@async begin
						while true
							idx = nextidx()
							if idx > NP
								break
							end
							FFT[idx] = remotecall_fetch(p, fc_remote, N, idx, RR_YOR[p], RR_PT[p], RR_sFFT[p])
						end
					end
				end
			end
		end
		return FFT
	end
end

# Parameters:
#	N::Int
#	- the size of the FFT being calculated is N
#	SNF::Array{Foat64, 1}
#	- SNF[i] is the value associated with the Permutation that permutation_index() maps to i
#	YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- output1() from yor()
#	PT::Array{Array{Array{Int, 1}, 1}, 1}
#	- output2() from yor()
#	C::Counter
#	- a wrapper for an Int that is used to coordinate the recursive access to the elements of SNF
# Return Values:
#	- FFT is the Fast Fourier Transfrom of SNF
function compute_fft(N::Int, SNF::Array{Float64, 1}, YOR::Array{Array{Array{SparseMatrixCSC, 1}, 1}, 1}, PT::Array{Array{Array{Int, 1}, 1}, 1}, C::Counter)
	if N == 1
		sFFT = Array(Array{Float64, 2}, 1)
		sFFTi = Array(Float64, 1, 1)
		sFFTi[1, 1] = SNF[C.N]
		sFFT[1] = sFFTi
		C.N += 1
		return sFFT
	end
	sFFT = Array(Array{Array{Float64, 2}, 1}, N)
	for n = 1:N
		sFFT[n] = compute_fft(N - 1, SNF, YOR, PT, C)
	end
	FFT = combine_sfft(N, YOR[N], PT[N], sFFT)
	return FFT
end

function compute_fft_remote(N::Int, SNF::Array{Float64, 1}, RR_YOR::RemoteRef, RR_PT::RemoteRef, C::Counter)
	YOR = fetch(RR_YOR)
	PT = fetch(RR_PT)
	FFT =  compute_fft(N, SNF, YOR, PT, C)
	return FFT
end

# Parameters:
#	N::Int
#	- the size of the FFT being calculated is N
#	YORn::Array{Array{SparseMatrixCSC, 1}, 1}, 1}
#	- Young's Orthogonal Representations for the Partitions of N
#	PTn::Array{Array{Int, 1}, 1}
#	- the decomposition indices for the Partitions of N
#	sFFT::Array{Array{Array{Float64, 2}, 1}, 1}
#	- the array of N FFT's of size N-1 that is used to compute this FFT of size N
# Return Values:
#	FFT::Array{Array{Float64, 2}, 1}
#	- FFT is the Fast Fourier Transform for this group
function combine_sfft(N::Int, YORn::Array{Array{SparseMatrixCSC, 1}, 1}, PTn::Array{Array{Int, 1}, 1}, sFFT::Array{Array{Array{Float64, 2}, 1}, 1})
	NP = length(YORn)
	FFT = Array(Array{Float64, 2}, NP)
	for p = 1:NP
		YORnp = YORn[p]
		PTnp = PTn[p]
		FC = fc(N, YORnp, PTnp, sFFT)
		FFT[p] = FC
	end
	return FFT
end

# Parameters:
#	N::Int
#	- the Fourier Component that is being calculated corresponds to a Partition, P, of N
#	YORnp::Array{SparseMatrixCSC, 1}
#	- Young's Orthogonal Representations for P
#	PTnp::Array{Int, 1}
#	- the indices for the Partitions that P decomposes into
#	sFFT::Array{Array{Array{Float64, 2}, 1}, 1}
#	- the array of N FFT's of size N-1 that is used to compute this FFT of size N
# Return Values:
#	FC::Array{Float64, 2}
#	- FC is the Fourier Coefficient corresponding to the Partition P
function fc(N::Int, YORnp::Array{SparseMatrixCSC, 1}, PTnp::Array{Int, 1}, sFFT::Array{Array{Array{Float64, 2}, 1}, 1})
	Dim = size(YORnp[1], 1)
	FC = dsm(Dim, sFFT[N], PTnp)
	CCM = eye(Dim)
	for n = (N - 1):-1:1
		CCM = YORnp[n] * CCM
		DSM = dsm(Dim, sFFT[n], PTnp)
		FC_n = CCM * DSM
		FC += FC_n
	end
	return FC
end

function fc_remote(N::Int, p::Int, RR_YOR::RemoteRef, RR_PT::RemoteRef, RR_sFFT::RemoteRef)
	YOR = fetch(RR_YOR)
	PT = fetch(RR_PT)
	sFFT = 	fetch(RR_sFFT)
	FC = fc(N, YOR[N][p], PT[N][p], sFFT)
	return FC
end

# Parameters:
#	Dim::Int
#	- the size of the DSM that is being caculated
#	sFFTn::Array{Array{Float64, 2}, 1}
#	- one of the elements of sFFT from fc()
#	PTnp::Array{Int, 1}
#	- same as in fc()
# Return Values:
#	DSM::Array{Float64, 2} (Direct Sum Matrix)
#	- the direct sum of the coefficients of sFFTn that correspond to the Partitions indicated by PTnp
function dsm(Dim::Int, sFFTn::Array{Array{Float64, 2}, 1}, PTnp::Array{Int, 1})
	DSM = zeros(Float64, Dim, Dim)
	offset = 0
	for i = 1:length(PTnp)
		FC = sFFTn[PTnp[i]]
		sDim = size(FC, 1)
		for r = 1:sDim
			for c = 1:sDim
				v = FC[r, c]
				if v != 0.0
					DSM[offset + r, offset + c] = v
				end
			end
		end
		offset += sDim
	end
	return DSM
end


text:
module GPMTD

## distributed packages
using Distances
using Distributions
using PDMats
using Random, Random123
using LinearAlgebra
using Distributed
using Dates
# import StatsBase: counts, sample, Weights
using StatsBase
using SpecialFunctions

## personal packages
using SparseProbVec

## files
include("Cor.jl")
include("General.jl")
include("Dist.jl")
include("Update_int.jl")
include("Update_Cov.jl")
include("Update_mixcomps.jl")
include("Update_alloc.jl")
include("Update_weights.jl")
include("Update_CovHypers.jl")
include("Lik.jl")
include("Sim.jl")
include("mcmc.jl")
include("DensityEstimation.jl")

end # module


text:
@testset "Test CIR simulation functions" begin
    @test_nowarn DualOptimalFiltering.rCIR(2, 0.5, 0.5, 0.5, 0.5, 0.5)
    @test_nowarn DualOptimalFiltering.generate_CIR_trajectory(range(0, stop = 2, length = 20), 3, 3., 0.5, 1);
    @test length(DualOptimalFiltering.generate_CIR_trajectory(range(0, stop = 2, length = 20), 3, 3., 0.5, 1)) == 20
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
using Beep
using Test

@testset "Beep.jl" begin
    @testset "resources" begin
    end

    @testset "functions" begin
        @test beep(:boing) === nothing
        @test beep(:back) === nothing
    end
end


text:
# stuct Som holds all important parameters of a trained som:
#
"""
    struct Som

Structure to hold all data of a trained SOM.

# Fields:
- `codes::Array{Float64,2}`: 2D-array of codebook vectors. One vector per row
- `colNames::Array{String,1}`: names of the attribute with which the SOM is trained
- `normParams::DataFrame`: normalisation parameters for each column
                of training data. Column headers corresponds with
                colNames.
- `norm::Symbol`: normalisation type; one of `:none, :minmax, :zscore`
- `xdim::Int`: number of neurons in x-direction
- `ydim::Int`: number of neurons in y-direction
- `nCodes::Int`: total number of neurons
- `grid::Array{Float64,2}`: 2D-array of coordinates of neurons on the map
          (2 columns (x,y)] for rectangular and hexagonal maps
           3 columns (x,y,z) for spherical maps)
- `indices::DataFrame`: X-, Y-indices of the neurons
- `topol::Symbol`: topology of the SOM; one of `:rectangular, :hexagonal, :spherical`
- `toroidal::Bool`: if `true`, the SOM is toroidal (has no edges)
- `population::Array{Int,1}`: 1D-array of numbers of training samples mapped to
                each neuron.
"""
struct Som
    codes::Array{Float64,2}
    colNames::Array{String}
    normParams::DataFrame
    norm::Symbol        # one of :none :minmax :zscore
    xdim::Int
    ydim::Int
    nCodes::Int
    grid::Array{Float64,2}
    indices::DataFrame
    topol::Symbol       # one of :rectangular :hexagonal :spherical
    toroidal::Bool
    population::Array{Int,1}

    Som(;codes::Array{Float64} = Array{Float64}(0),
        colNames::Array{String,1} = Array{String}(0),
        normParams::DataFrame = DataFrame(),
        norm::Symbol = :none,
        xdim::Int = 1,
        ydim::Int = 1,
        nCodes::Int = 1,
        grid::Array{Float64,2} = zeros(1,1),
        indices::DataFrame = DataFrame(),
        topol::Symbol = :hexagonal,
        toroidal::Bool = false,
        population::Array{Int,1} = [1]) = new(codes,
                                              colNames,
                                              normParams,
                                              norm,
                                              xdim,
                                              ydim,
                                              nCodes,
                                              grid,
                                              indices,
                                              topol,
                                              toroidal,
                                              population)
end


text:
export EmptyGridWorld

mutable struct EmptyGridWorld <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
end

function EmptyGridWorld(;n=8, agent_start_pos=CartesianIndex(2,2), agent_start_dir=RIGHT)
    objects = (EMPTY, WALL, GOAL)
    w = GridWorldBase(objects, n, n)
    w[EMPTY, 2:n-1, 2:n-1] .= true
    w[WALL, [1,n], 1:n] .= true
    w[WALL, 1:n, [1,n]] .= true
    w[GOAL, n-1, n-1] = true
    w[EMPTY, n-1, n-1] = false
    EmptyGridWorld(w, agent_start_pos, Agent(dir=agent_start_dir))
end

function (w::EmptyGridWorld)(::MoveForward)
    dir = get_dir(w.agent)
    dest = dir(w.agent_pos)
    if !w.world[WALL, dest]
        w.agent_pos = dest
    end
    w
end


text:
using Test, JLD2, UnicodePlots
using LinearAlgebra, ForwardDiff, Arpack
using SparseArrays, SparseDiffTools, Symbolics

# load mesh
meshfile = "$(@__DIR__)/data/squaremesh.jld2"
JLD2.@load(meshfile, p, t, ps, ts)
#p, t =  ps, ts

np, nt = size(p, 1), size(t, 1)
ε = 1e-4
x, dp = p, 2 * rand(size(p)...) .- 1
Ib = unique(FemAdjoint.btri(t)[:])
dp[Ib, :] .= 0.# fix Dirichlet points
xp, xm = p .+ ε * dp, p .- ε * dp

# finite differences
f = x -> FemAdjoint.costnormU(x, t)
fp, fm = f(xp), f(xm)
@time g = FemAdjoint.∇costnormU_direct(x, t, plotsolutions = true)
@time g = FemAdjoint.∇costnormU_direct(x, t, plotsolutions = false)
println("")
println(fp)
println(fm)
println(" ")
println((fp - fm) / (2 * ε))
println(dot(g, dp))
println("")
@testset "test eval direct gradient dot function for costnormU" begin
    @test abs(dot(g, dp) - (fp - fm) / (2 * ε)) < 1e-4
end


text:
mutable struct RINEXHeader
    format_version::Float64
    file_type::String
    satellite_system::String
    pgm::String
    run_by::String
    date::String
    marker_name::String
    marker_number::String
    observer::String
    agency::String
    ionospheric_corrections::String
    time_system_corrections::String
    leap_seconds::String
    comments::String
end

# empty constructor
function RINEXHeader()
    RINEXHeader(
        0.0,
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        "",
        ""
    )
end

function set_format_version(RINEXHeader, set_format_version)
end


text:
module GraphJuMP

using JuMP
using LightGraphs
using MetaGraphs

export setgraph!, removegraph!, addpath!

const extKey = :GraphJuMP

"""
Several formulations of the flows are possible:

* `EdgeFlow`: each edge of the graph is explicitly modelled. This model tends to scale poorly to large number of nodes
  and edges in the graph
* `Path`: only some paths in the graph are modelled. This formulation is typical for column-generation schemes
  (set `delayed=true` for this)

To create a new formulation, subtype `GraphFormulation` and `GraphFormulationExt`; implement the required functions:

* `supportspath(::GraphFormulation)`: whether this formulation uses the notion of path. This implies that functions like
  `addpath!` must be supported.
* `buildformulation(m::Model, gf::GraphFormulation)`: builds the formulation, i.e. creates the required JuMP variables and
  constraints for the formulation.
"""
abstract type GraphFormulation end
struct EdgeFlow <: GraphFormulation end
struct Path <: GraphFormulation end

supportspath(::GraphFormulation) = false
supportspath(::Path) = true

"""
Each formulation should have its subtype of `GraphFormulationExt` to store its internal parameters.
"""
abstract type GraphFormulationExt end

mutable struct EdgeFlowFormulationExt <: GraphFormulationExt
  flows::Matrix{Variable} # First index: edge; second index: commodity.
  id_to_edge::Dict{Int, Edge} # Mapping an integer index (used in the model) to a graph edge.
  edge_to_id::Dict{Edge, Int} # Reverse mapping: from a graph edge to an integer index.
end

"""
Stores the internal state of the graph model.
"""
mutable struct GraphModelExt{T <: Int, U <: Real}
  graph::MetaDiGraph{T, U} # The graph below this model. Must be a "meta" graph, as many properties might be useful for each edge or vertex (such as max capacity).
  formulation::GraphFormulation # A graph formulation, see GraphFormulation enumeration.
  delayed::Bool # For a path formulation, controls when paths are computed: false, before solving; true, after solving.
  flowcategory::Symbol # The JuMP category for flow variables. For now, these are: :Cont, :Int, :Bin, :SemiCont, :SemiInt, :SDP
  commodities::Vector{String} # Whether several commodities are present.
  formulationext::Union{GraphFormulationExt, Void} # The formulation object, once it is built (after calling `buildformulation`).

  function GraphModelExt{T, U}(graph::MetaDiGraph{T, U}, formulation::GraphFormulation=EdgeFlow(),
                               delayed::Bool=false, flowcategory::Symbol=:Cont, commodities::Vector{String}=String["default"]
                              ) where{T, U}
    return new(graph, formulation, delayed, flowcategory, commodities, nothing)
  end
end

# Link between GraphModelExt and JuMP's model.
hasgraph(m::Model) = haskey(m.ext, extKey)
checkgraph(m::Model) = hasgraph(m) || error("The model has no associated graph model.")
getgraphext(m::Model) = hasgraph(m) && m.ext[extKey]
getgraph(m::Model) = getgraphext(m).graph

hasformulation(m::Model) = getgraphext(m).formulationext !== nothing
checkformulation(m::Model) = hasformulation(m) || error("The model has no associated graph formulation. " *
  "Did you forget to call `buildformulation` or to link the graph model to a JuMP model with `setgraph!`?")

function setgraph!(m::Model, graph::GraphModelExt)
  if hasgraph(m)
    warn("Model already has an associated graph model! It will be replaced.")
  end
  m.ext[extKey] = graph
  m.ext[extKey].formulationext = buildformulation(m)
end

setgraph!{T, U}(m::Model, graph::MetaDiGraph{T, U}, args...) = setgraph!(m, GraphModelExt{T, U}(graph, args...))
setgraph!{T}(m::Model, graph::AbstractMetaGraph{T}, args...) = error("GraphJuMP only uses graphs with metadata from the MetaGraphs.jl package. " *
  "You can convert your `graph` to a MetaGraph using `MetaGraph(graph)`. ")

function removegraph!(m::Model)
  if hasgraph(m)
    delete!(m.ext, extKey)
  end
  # No error if there is no graph.
end

# Build the formulations.
function buildformulation(m::Model) # Makes a dispatch on the formulation.
  checkgraph(m)
  return buildformulation(m, getgraphext(m).formulation)
end

function buildformulation(m::Model, gf::EdgeFlow)
  # Create the flow variables, one per edge and per commodity.
  n_edges = ne(getgraph(m))
  n_commodities = length(getgraphext(m).commodities)
  flows = @variable(m, [e=1:n_edges, c=1:n_commodities], category=getgraphext(m).flowcategory)

  # Map edge and their IDs in the model.
  id_to_edge = Dict([idx => edge for (idx, edge) in enumerate(edges(getgraph(m)))])
  edge_to_id = map(reverse, id_to_edge)

  # Give names to the flow variables.
  for e in 1:n_edges
    for c in 1:n_commodities
      edge = id_to_edge[e]
      setname(flows[e, c], "flow_" * string(edge.src) * "_to_" * string(edge.dst) * "_commodity_" * string(getgraphext(m).commodities[c]))
    end
  end

  # Done!
  return EdgeFlowFormulationExt(flows, id_to_edge, edge_to_id)
end

# For path formulations.
function addpath!(m::Model, path)
  checkgraph(m)

  if ! supportspath(getgraphext(m).formulation)
    error("addpath! can only be called on graph models using a path-based formulation, such as `Path`.")
  end

  # TODO!
end

end


text:
export expmv

function expmv(t, A, b; M = [], prec = "double", shift = false, full_term = false, prnt = false)
               # bal = false,

    #EXPMV   Matrix exponential times vector or matrix.
    #   [F,S,M,MV,MVD] = EXPMV(t,A,B,[],PREC) computes EXPM(t*A)*B without
    #   explicitly forming EXPM(t*A). PREC is the required accuracy, 'double',
    #   'single' or 'half', and defaults to CLASS(A).
    #
    #   A total of MV products with A or A^* are used, of which MVD are
    #   for norm estimation.
    #
    #   The full syntax is
    #
    #     [f,s,m,mv,mvd,unA] = expmv(t,A,b,M,prec,shift,bal,full_term,prnt).
    #
    #   unA = 1 if the alpha_p were used instead of norm(A).
    #
    #   If repeated invocation of EXPMV is required for several values of t
    #   or B, it is recommended to provide M as an external parameter as
    #   M = SELECT_TAYLOR_DEGREE(A,m_max,p_max,prec,shift,bal,true).
    #   This also allows choosing different m_max and p_max.

    #   Reference: A. H. Al-Mohy and N. J. Higham, Computing the action of
    #   the matrix exponential, with an application to exponential
    #   integrators. MIMS EPrint 2010.30, The University of Manchester, 2010.

    #   Awad H. Al-Mohy and Nicholas J. Higham, October 26, 2010.

    # if bal
    #     [D,B] = balance(A)
    #     if norm(B,1) < norm(A,1)
    #         A = B
    #         b = D\b
    #     else
    #         bal = false
    #     end
    # end

    n = size(A,1)

    if shift
        mu = trace(A)/n
        #mu = full(mu) # Much slower without the full!
        A = A-mu*speye(n)
    end

    if isempty(M)
        tt = 1
        (M,mvd,alpha,unA) = select_taylor_degree(t*A,b)
        mv = mvd
    else
        tt = t
        mv = 0
        mvd = 0
    end

    tol =
      if prec == "double"
          2.0^(-53)
      elseif prec == "single"
          2.0^(-24)
      elseif prec == "half"
          2.0^(-10)
      end

    s = 1;

    if t == 0
        m = 0;
    else
        (m_max,p) = size(M);
        U = diagm(1:m_max);
        C = ( (ceil(abs(tt)*M))'*U );
        zero_els = find(x->x==0, C)
        for el in zero_els
            C[el] = Inf
        end
        if p > 1
            cost,m = findmin(minimum(C,1)); # cost is the overall cost.
        else
            cost,m = findmin(C);  # when C is one column. Happens if p_max = 2.
        end
        if cost == Inf
            cost = 0
        end
        s = max(cost/m,1);
    end

    eta = 1;

    if shift
        eta = exp(t*mu/s)
    end

    f = b;

    # if prnt
    #     fprintf("m = %2.0f, s = %g, m_actual = ", m, s)
    # end

    for i = 1:s
        c1 = norm(b,Inf);
        for k = 1:m
            b = (t/(s*k))*(A*b);
            mv = mv + 1;
            f =  f + b;
            c2 = norm(b,Inf);
            if !full_term
                if c1 + c2 <= tol*norm(f,Inf)
                    # if prnt
                    #     fprintf(" %2.0f, ", k)
                    # end
                    break
                end
                c1 = c2;
            end

        end
        f = eta*f
        b = f
    end

    # if prnt
    #     fprintf("\n")
    # end

    #if bal
    #    f = D*f
    #end

    #return (f,s,m,mv,mvd,unA)
    return f
end


text:
# readme example

#addprocs(2)
using Caesar, RoME, TransformUtils, Distributions

# load scene and ROV model (might experience UDP packet loss LCM buffer not set)
vc = startdefaultvisualization()
sc1 = loadmodel(:scene01); sc1(vc)
rovt = loadmodel(:rov); rovt(vc)

initCov = 0.001*eye(6); [initCov[i,i] = 0.00001 for i in 4:6];
odoCov = 0.0001*eye(6); [odoCov[i,i] = 0.00001 for i in 4:6];
rangecov, bearingcov = 3e-4, 2e-3

# start and add to a factor graph
fg = identitypose6fg(initCov=initCov)
tf1 = SE3([0.0;0.0;0.0], Euler(0,0.0,0.0) )
tf2 = SE3([0.0;0.0;0.0], Euler(pi/2,0.0,0.0) )
tf3 = SE3([0.0;0.0;0.0], Euler(-pi/2,0.0,0.0) )

tf4 = SE3([0.0;5.0;0.0], Euler(0,0.0,0.0) )
tf5 = SE3([4.0;4.0;0.0], Euler(0,0.0,-pi/2) )

tf6 =  SE3([0.0;5.0;0.0], Euler(0,0.0,0.0) )
tf7 = SE3([4.0;4.0;0.0], Euler(0,0.0,-pi/2) )

tf8 =  SE3([0.0;5.0;0.0], Euler(0,0.0,0.0) )
tf9 = SE3([4.0;4.0;0.0], Euler(0,0.0,-pi/2) )

addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf1), odoCov) ) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf2), odoCov) ) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf3), odoCov) ) )

addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf4), odoCov) ) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf5), odoCov) ) )

addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf6), odoCov) ) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf7), odoCov) ) )

addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf8), odoCov) ) )
addOdoFG!(fg, Pose3Pose3(MvNormal(veeEuler(tf9), odoCov) ) )

visualizeallposes!(vc, fg, drawlandms=false)

# addLinearArrayConstraint(fg, (4.0, 0.0), :x0, :l1, rangecov=rangecov,bearingcov=bearingcov)
# visualizeDensityMesh!(vc, fg, :l1)
# addLinearArrayConstraint(fg, (4.0, 0.0), :x1, :l1, rangecov=rangecov,bearingcov=bearingcov)
# addLinearArrayConstraint(fg, (4.0, 0.0), :x2, :l1, rangecov=rangecov,bearingcov=bearingcov)
# visualizeDensityMesh!(vc, fg, :l1)

addLinearArrayConstraint(fg, (4.0, 0.0), :x4, :l2, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (4.0, 0.0), :x5, :l2, rangecov=rangecov,bearingcov=bearingcov)

addLinearArrayConstraint(fg, (4.0, 0.0), :x6, :l3, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (4.0, 0.0), :x7, :l3, rangecov=rangecov,bearingcov=bearingcov)

addLinearArrayConstraint(fg, (4.0, 0.0), :x8, :l4, rangecov=rangecov,bearingcov=bearingcov)
addLinearArrayConstraint(fg, (4.0, 0.0), :x9, :l4, rangecov=rangecov,bearingcov=bearingcov)

solveTree!(fg)

# visualize(fg, vc, drawlandms=true, densitymeshes=[:l1;:x2])
visualize(fg, vc, drawlandms=true, densitymeshes=[:l2;:x4])
visualize(fg, vc, drawlandms=true, densitymeshes=[:l3;:x7])
visualize(fg, vc, drawlandms=true, densitymeshes=[:l4;:x9])

Graphs.plot(fg.g)

#


text:
using ControlSystems, RobustAndOptimalControl
plotly(size=(1850, 900), show=false)
##
w = 2pi*exp10.(LinRange(-2, 2, 400))
function flexible_servo_loop(; kp=1, kv=1, ki=1, Tf=0.1, Jm=1, Jl=10, c=1, k=10000)

    A = [
        0.0 1 0 0
        -k/Jm -c/Jm k/Jm c/Jm
        0 0 0 1
        k/Jl c/Jl -k/Jl -c/Jl
    ]
    B = [0, 1/Jm, 0, 0]
    C = I(4)
    D11 = D12 = D21 = D22 = 0
    nx,nu,ny = 4,1,1
    P = named_ss(ss(A,B,C,0), u=:u, x=[:qm, :qdm, :qa, :qda], y = [:qm, :qdm, :qa, :qda])

    K = let
        K = cascade_controller(; kp, kv, ki, Tf) # TODO: extend to take in several measurements instead of differentiating
        K = named_ss(K, u=:e, y = :Cu)
        # extended_controller(K)
    end

    sumE = sumblock("e = r - y")

    u1 = [:u,  :y, :e]
    y1 = [:Cu, :qm, :e]
    w1 = [:r]
    G = connect([P, K, sumE]; y1, u1, w1)
    #F = ss([tf(0, [1, 0]); tf(1)])

    (; P, K, G)
end

function cascade_controller(; kp=1, kv=1, ki=1, Tf=0.1)
    ss(tf([kv, kp*kv + ki, kp*ki], [Tf, 1, 0]))
end

#
P,K,G = flexible_servo_loop(kp=0.1, kv=250, ki=15, Tf=0.001)

plot(
    # plot(step(c2d(G.sys, 0.001), 3)),
    plot(step(G, 3)),
    bodeplot(P, w, hz=true, plotphase=false),
    bodeplot(G, w, hz=true, plotphase=false),
    layout=(1,3),
    title="",
)
display(current())