text:
module MathOptInterfaceToMathProgBase

using MathOptInterface
const MOI = MathOptInterface
const MOIU = MOI.Utilities

import MathProgBase
const MPB = MathProgBase.SolverInterface

MOIU.@model(InnerModel,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (),
    (MOI.SingleVariable,),
    (),
    (),
    ()
)

# Make `Model` a constant for use in the rest of the file.
const Model = MOIU.UniversalFallback{InnerModel{Float64}}

"Attribute for the MathProgBase solver."
struct MPBSolver <: MOI.AbstractOptimizerAttribute end

"Attribute for the MathProgBase status."
struct MPBSolutionAttribute <: MOI.AbstractModelAttribute end

"Struct to contain the MPB solution."
struct MPBSolution
    status::Symbol
    is_minimization::Bool
    objective_value::Float64
    primal_solution::Dict{MOI.VariableIndex, Float64}
end

"""
    Optimizer(; solver)

# Example

    model = Optimizer(solver = AmplNLWriter("/usr/bin/bonmin"))
"""
function Optimizer(; solver)
    model = MOIU.UniversalFallback(InnerModel{Float64}())
    MOI.set(model, MPBSolver(), solver)
    return model
end

Base.show(io::IO, ::Model) = println(io, "A MathProgBase model")

# We re-define is_empty and empty! to prevent the universal fallback from
# deleting the solver that we are caching in it.
function MOI.is_empty(model::Model)
    return MOI.is_empty(model.model) &&
        isempty(model.constraints) &&
        isempty(model.modattr) &&
        isempty(model.varattr) &&
        isempty(model.conattr) &&
        length(model.optattr) == 1 &&
        haskey(model.optattr, MPBSolver())
end

function MOI.empty!(model::Model)
    MOI.empty!(model.model)
    empty!(model.constraints)
    model.nextconstraintid = 0
    empty!(model.con_to_name)
    model.name_to_con = nothing
    empty!(model.modattr)
    empty!(model.varattr)
    empty!(model.conattr)
    mpb_solver = model.optattr[MPBSolver()]
    empty!(model.optattr)
    model.optattr[MPBSolver()] = mpb_solver
    return
end

set_to_bounds(set::MOI.LessThan) = (-Inf, set.upper)
set_to_bounds(set::MOI.GreaterThan) = (set.lower, Inf)
set_to_bounds(set::MOI.EqualTo) = (set.value, set.value)
set_to_bounds(set::MOI.Interval) = (set.lower, set.upper)
set_to_cat(set::MOI.ZeroOne) = :Bin
set_to_cat(set::MOI.Integer) = :Int

struct NLPEvaluator{T <: MOI.AbstractNLPEvaluator} <: MPB.AbstractNLPEvaluator
    inner::T
    variable_map::Dict{MOI.VariableIndex, Int}
end

"""
MathProgBase expects expressions with variables denoted by `x[i]` for contiguous
`i`. However, JuMP 0.19 creates expressions with `x[MOI.VariableIndex(i)]`. So
we have to recursively walk the expression replacing instances of
MOI.VariableIndex by a corresponding integer.
"""
function replace_variableindex_by_int(variable_map, expr::Expr)
    for (i, arg) in enumerate(expr.args)
        expr.args[i] = replace_variableindex_by_int(variable_map, arg)
    end
    return expr
end
function replace_variableindex_by_int(variable_map, expr::MOI.VariableIndex)
    return variable_map[expr]
end
replace_variableindex_by_int(variable_map, expr) = expr

function MPB.initialize(d::NLPEvaluator, requested_features::Vector{Symbol})
    MOI.initialize(d.inner, requested_features)
    return
end

function MPB.features_available(d::NLPEvaluator)
    return MOI.features_available(d.inner)
end

function MPB.eval_f(d::NLPEvaluator, x)
    return MOI.eval_objective(d.inner, x)
end

function MPB.eval_g(d::NLPEvaluator, g, x)
    return MOI.eval_constraint(d.inner, g, x)
end

function MPB.eval_grad_f(d::NLPEvaluator, g, x)
    return MOI.eval_objective_gradient(d.inner, g, x)
end

function MPB.jac_structure(d::NLPEvaluator)
    return MOI.jacobian_structure(d.inner)
end

function MPB.obj_expr(d::NLPEvaluator)
    expr = MOI.objective_expr(d.inner)
    return replace_variableindex_by_int(d.variable_map, expr)
end

function MPB.constr_expr(d::NLPEvaluator, i)
    expr = MOI.constraint_expr(d.inner, i)
    return replace_variableindex_by_int(d.variable_map, expr)
end

function MOI.optimize!(model::Model)
    mpb_solver = MOI.get(model, MPBSolver())

    # Get the variable data.
    variables = MOI.get(model, MOI.ListOfVariableIndices())
    num_var = length(variables)
    x_l = fill(-Inf, num_var)
    x_u = fill(Inf, num_var)
    x_cat = fill(:Cont, num_var)

    variable_map = Dict{MOI.VariableIndex, Int}()
    for (i, variable) in enumerate(variables)
        variable_map[variable] = i
    end

    for set_type in (MOI.LessThan{Float64}, MOI.GreaterThan{Float64})
        for c_ref in MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())
            c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref)
            c_set = MOI.get(model, MOI.ConstraintSet(), c_ref)
            v_index = variable_map[c_func.variable]
            lower, upper = set_to_bounds(c_set)
            x_l[v_index] = lower
            x_u[v_index] = upper
        end
    end

    for set_type in (MOI.ZeroOne, MOI.Integer)
        for c_ref in MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())
            c_func = MOI.get(model, MOI.ConstraintFunction(), c_ref)
            c_set = MOI.get(model, MOI.ConstraintSet(), c_ref)
            v_index = variable_map[c_func.variable]
            x_cat[v_index] = set_to_cat(c_set)
        end
    end

    # Get the optimzation sense.
    opt_sense = MOI.get(model, MOI.ObjectiveSense())
    sense = opt_sense == MOI.MAX_SENSE ? :Max : :Min

    nlp_block = try
        MOI.get(model, MOI.NLPBlock())
    catch ex
        error("Expected a NLPBLock.")
    end

    # Extract constraint bounds.
    num_con = length(nlp_block.constraint_bounds)
    g_l = fill(-Inf, num_con)
    g_u = fill(Inf, num_con)
    for (i, bound) in enumerate(nlp_block.constraint_bounds)
        g_l[i] = bound.lower
        g_u[i] = bound.upper
    end

    mpb_model = MPB.NonlinearModel(mpb_solver)
    MPB.loadproblem!(
        mpb_model, num_var, num_con, x_l, x_u, g_l, g_u, sense,
        NLPEvaluator(nlp_block.evaluator, variable_map)
    )

    MPB.setvartype!(mpb_model, x_cat)

    MPB.optimize!(mpb_model)

    # MathProbBase solution
    primal_solution = Dict{MOI.VariableIndex, Float64}()
    for (variable, sol) in zip(variables, MPB.getsolution(mpb_model))
        primal_solution[variable] = sol
    end
    MOI.set(model, MPBSolutionAttribute(), MPBSolution(
        MPB.status(mpb_model),
        sense == :Min,
        MPB.getobjval(mpb_model),
        primal_solution
    ))
    return
end

function MOI.get(model::Model, ::MOI.VariablePrimal, var::MOI.VariableIndex)
    mpb_solution = MOI.get(model, MPBSolutionAttribute())
    if mpb_solution === nothing
        return nothing
    end
    return mpb_solution.primal_solution[var]
end

function MOI.get(model::Model, ::MOI.ObjectiveValue)
    mpb_solution = MOI.get(model, MPBSolutionAttribute())
    if mpb_solution === nothing
        return nothing
    end
    if mpb_solution.is_minimization
        return -mpb_solution.objective_value
    else
        return mpb_solution.objective_value
    end
end

function MOI.get(model::Model, ::MOI.TerminationStatus)
    mpb_solution = MOI.get(model, MPBSolutionAttribute())
    if mpb_solution === nothing
        return MOI.OPTIMIZE_NOT_CALLED
    end
    status = mpb_solution.status
    if status == :Optimal
        return MOI.OPTIMAL
    elseif status == :Infeasible
        return MOI.INFEASIBLE
    elseif status == :Unbounded
        return MOI.DUAL_INFEASIBLE
    elseif status == :UserLimit
        return MOI.OTHER_LIMIT
    elseif status == :Error
        return MOI.OTHER_ERROR
    end
    return MOI.OTHER_ERROR
end

function MOI.get(model::Model, ::MOI.PrimalStatus)
    mpb_solution = MOI.get(model, MPBSolutionAttribute())
    if mpb_solution === nothing
        return MOI.NO_SOLUTION
    end
    status = mpb_solution.status
    if status == :Optimal
        return MOI.FEASIBLE_POINT
    end
    return MOI.NO_SOLUTION
end

function MOI.get(model::Model, ::MOI.DualStatus)
    return MOI.NO_SOLUTION
end

end


text:
import DataFrames: DataFrame, DataFrameRow
import Flatten: flatten
import Dates

using Mixers
using Agtor

abstract type AgParameter end

AgUnion = Union{Date, Int64, Float64, AgParameter}

@with_kw mutable struct ConstantParameter <: AgParameter
    name::String
    default_val::Any
    value::Any

    ConstantParameter(name, default_val) = new(name, default_val, default_val)
end

@with_kw mutable struct RealParameter <: AgParameter
    name::String
    min_val::Float64
    max_val::Float64
    default_val::Float64
    value::Float64

    RealParameter(name, min_val, max_val, d_val, value) = new(name, min_val, max_val, d_val, value)
    RealParameter(name, min_val, max_val, value) = new(name, min_val, max_val, value, value)
    RealParameter(name, range_vals::Array, value) = new(name, range_vals..., value, value)
end

@doc """
    Categorical parameters.

    Min and max values will map to (integer) element position in an CategoricalArray.

    Sampling between the min and max values will be mapped to their categorical value in the given array
    using `floor()` of the Float value `x`.

    Valid values for the CategoricalParameter will therefore be:
    `1 <= x < (n+1)`, where n is number of options.
"""
@with_kw mutable struct CategoricalParameter <: AgParameter
    name::String
    cat_val::CategoricalArray
    min_val::Int64
    max_val::Int64
    default_val::Any
    value::Any

    function CategoricalParameter(name::String, cat_val::CategoricalArray, value::Any)::CategoricalValue
        min_val = 1
        max_val = length(cat_val) + 1

        return new(name, cat_val, min_val, max_val, value, value)
    end

    function CategoricalParameter(name::String, cat_val::CategoricalArray, default_value::Any, value::Any)::CategoricalValue
        min_val = 1
        max_val = length(cat_val) + 1

        return new(name, cat_val, min_val, max_val, default_value, value)
    end
end

"""Setter for CategoricalParameter to handle float to categorical element position."""
function Base.setproperty!(cat::CategoricalParameter, v::Symbol, value)::Nothing
    if v == :value
        if value isa AbstractFloat || value isa Integer
            cat_pos = floor(Int, value)
            cat.value = cat.cat_val[cat_pos]
        elseif value isa String
            pos = cat.cat_val(findfirst(x-> x == value, cat.cat_val))
            if isnothing(pos)
                throw(ArgumentError("$(pos) is not a valid option for $(cat.name)"))
            end

            cat.value = cat.cat_val[pos]
        else
            throw(ArgumentError("Type $(typeof(value)) is not a valid type option for $(cat.name), has to be Integer, Float or String"))
        end
    else
        setfield!(f, Symbol(v), value)
    end

    return nothing
end

"""Ag specific in function.

Purposeful decision not to override Base.in() method as this
performs the check on any Array expecting NamedTuples.
"""
function agin(p::AgParameter, store::Array)::Bool
    for st in store
        if st.name == p.name
            return true
        end
    end

    return false
end

# Below is equivalent to defining methods for each operation, e.g:
#    function Base.:+(x::AgParameter, y::AgParameter) x.value + y.value end
for op = (:+, :-, :/, :*, :^, :<, :>)
    @eval Base.$op(x::AgParameter, y::AgParameter) = Base.$op(x.value, y.value)
end

for op = (:+, :-, :/, :*, :^, :<, :>)
    @eval Base.$op(x::Number, y::AgParameter) = Base.$op(x, y.value)
    @eval Base.$op(x::AgParameter, y::Number) = Base.$op(x.value, y)
end

for op = (:isless, :isgreater)
    @eval Base.$op(x::AgParameter, y::Number) = Base.$op(x.value, y)
end

Base.:*(x::String, y::Union{ConstantParameter, CategoricalParameter}) = x * y.value

Base.convert(::Type{T}, x::AgParameter) where {T<:Number} = T(x.value)
Base.convert(x::Type{Any}, y::Agtor.ConstantParameter) = convert(x, y.value)

for op = (:Year, :Month, :Day)
    @eval Dates.$op(x::AgParameter) = Dates.$op(x.value)
end

"""Returns min/max values"""
function min_max(p::AgParameter)
    if is_const(p)
        return p.value
    end

    return p.min_val, p.max_val
end

"""Returns min/max values"""
function min_max(dataset::Dict)::Dict
    mm::Dict{Union{Dict, Symbol}, Union{Dict, Any}} = Dict()
    for (k, v) in dataset
        if v isa AgParameter
            mm[k] = min_max(v)
        else
            mm[k] = v
        end
    end

    return mm
end

"""
Generate AgParameter definitions.

Parameters
----------
prefix : str

dataset : Dict, of parameters for given component

Returns
----------
* Dict matching structure of dataset
"""
function generate_agparams(prefix::Union{String, Symbol}, dataset::Dict)::Dict{Symbol, Any}
    comp_sep = Agtor.component_sep
    name_sep = Agtor.component_name_sep
    field_sep = Agtor.component_field_sep

    if "component" in keys(dataset)
        comp = dataset["component"]
        comp_name = dataset["name"]

        if prefix === ""
            prefix *= "$(comp)$(name_sep)$(comp_name)"
        else
            prefix *= "$(comp_sep)$(comp)$(name_sep)$(comp_name)"
        end
    end

    created::Dict{Any, Any} = deepcopy(dataset)
    for (n, vals) in dataset
        var_id = "$(prefix)$(field_sep)$n"

        s = Symbol(n)
        pop!(created, n)

        if isa(vals, Dict)
            created[s] = generate_agparams(prefix, vals)
            continue
        elseif endswith(String(n), "_spec")
            # Recurse into other defined specifications
            tmp = load_yaml(vals)
            created[s] = generate_agparams(prefix, tmp)
            continue
        end

        if !in(vals[1], ["CategoricalParameter", "RealParameter", "ConstantParameter"])
            created[s] = vals
            continue
        end

        if isa(vals, Array)
            val_type, param_vals = vals[1], vals[2:end]

            def_val = param_vals[1]

            if length(unique(param_vals)) == 1
                created[s] = ConstantParameter(var_id, def_val)
                continue
            end

            valid_vals = param_vals[2:end]
            if val_type == "CategoricalParameter"
                valid_vals = categorical(valid_vals, ordered=true)
            end

            created[s] = eval(Symbol(val_type))(var_id, valid_vals, def_val)
        else
            created[s] = vals
        end
    end

    return created
end

"""Extract parameter values from AgParameter"""
function extract_values(p::AgParameter; prefix::Union{String, Nothing}=nothing)::NamedTuple
    if !isnothing(prefix)
        name = prefix * p.name
    else
        name = p.name
    end

    if is_const(p)
        return (name=name, ptype=typeof(p), default=p.value, min_val=p.value, max_val=p.value)
    end

    return (name=name, ptype=typeof(p), default=p.default_val, min_val=p.min_val, max_val=p.max_val)
end

function collect_agparams!(dataset::Dict, store::Array; ignore::Union{Array, Nothing}=nothing)::DataFrame
    for (k, v) in dataset
        if !isnothing(ignore) && (k in ignore)
            continue
        end

        if v isa AgParameter && !is_const(v)
            if !agin(v, store)
                push!(store, extract_values(v))
            end
        elseif v isa Dict || v isa Array
            collect_agparams!(v, store; ignore=ignore)
        end
    end

    return DataFrame(store)
end

"""Extract parameter values from Dict specification."""
function collect_agparams(dataset::Dict; ignore::Union{Array, Nothing}=nothing)::DataFrame
    return collect_agparams!(dataset, []; ignore=ignore)
end

function collect_agparams!(dataset::Union{Array, Tuple}, store::Array; ignore::Union{Array, Nothing}=nothing)::Nothing
    for v in dataset
        if v isa AgParameter && !is_const(v)
            if !agin(v, store)
                push!(store, extract_values(v))
            end
        elseif v isa Dict || v isa Array || v isa AgComponent || v isa Tuple
            collect_agparams!(v, store; ignore=ignore)
        end
    end

    return nothing
end

"""Extract parameter values from AgComponent."""
function collect_agparams!(dataset::Union{AgComponent, AgParameter}, store::Array; ignore=nothing)::DataFrame
    match = (Array, Tuple, Dict, AgComponent, AgParameter)

    for fn in fieldnames(typeof(dataset))

        fv = getfield(dataset, fn)
        if fv isa AgParameter && !is_const(fv)
            if !agin(fv, store)
                push!(store, extract_values(fv))
            end
        elseif any(map(x -> fv isa x, match))
            if fv isa Array{Date}
                continue
            end
            collect_agparams!(fv, store; ignore=ignore)
        end
    end

    return DataFrame(store)
end

function collect_agparams!(dataset::Union{AgComponent, AgParameter}; ignore=nothing)::DataFrame
    return collect_agparams!(dataset, []; ignore=nothing)
end

"""Extract parameter values from Dict specification and store in a common Dict."""
function collect_agparams!(dataset::Dict, mainset::Dict)::Nothing
    for (k, v) in dataset
        if Symbol(v.name) in mainset
            throw(ArgumentError("$(v.name) is already set!"))
        end

        if v isa AgParameter && !is_const(v)
            mainset[Symbol(v.name)] = extract_values(v)
        end
    end

    return
end

"""Checks if parameter is constant."""
function is_const(p::AgParameter)::Bool
    if (p isa ConstantParameter) || (value_range(p) == 0.0)
        return true
    end

    return false
end

"""Returns max - min, or 0.0 if no min value defined"""
function value_range(p::AgParameter)::Float64
    if hasproperty(p, :min_val) == false
        return 0.0
    end

    return p.max_val - p.min_val
end

"""Modify AgParameter name in place by adding a prefix."""
function add_prefix!(prefix::String, component)::Nothing
    params = flatten(component, AgParameter)
    for pr in params
        setfield!(pr, :name, prefix*pr.name)
    end

    return
end

"""
Usage:
    zone_dir = "data_dir/zones/"
    zone_specs = load_yaml(zone_dir)
    zone_params = generate_agparams("", zone_specs["Zone_1"])

    collated_specs = []
    collect_agparams!(zone_params, collated_specs; ignore=["crop_spec"])

    # Expect only CSV for now...
    climate_fn::String = "data/climate/farm_climate_data.csv"
    climate::Climate = load_climate(climate_fn)
    z1 = create(zone_params, climate)

    param_info = DataFrame(collated_specs)

    # Generate dataframe of samples
    samples = sample(param_info, 1000, sampler)  # where sampler is some function

    # Update z1 components with values found in first row
    set_params!(z1, samples[1])
"""
function set_params!(comp, sample)::Nothing
    match = (Array, Tuple, Dict, AgComponent, AgParameter)
    for f_name in fieldnames(typeof(comp))
        tmp_f = getfield(comp, f_name)

        if any(map(x -> isa(tmp_f, x), match))
            set_params!(tmp_f, sample)
        end
    end

    return
end

function set_params!(comp::Array, sample)::Nothing
    arr_type = eltype(comp)
    tmp_flat = reduce(vcat, Flatten.flatten(comp, Array{arr_type}))
    for i in tmp_flat
        set_params!(i, sample)
    end

    return
end

function set_params!(comp::Dict, sample)::Nothing
    for (k,v) in comp
        set_params!(v, sample)
    end

    return
end

function set_params!(p::AgParameter, sample::Union{DataFrame, DataFrameRow})::Nothing
    p_name::Symbol = Symbol(p.name)
    if p_name in names(sample)
        setfield!(p, :value, sample[p_name])
    end

    return
end

function set_params!(p::AgParameter, sample::Union{Number, String})::Nothing
    setfield!(p, :value, sample)

    return
end

function set_params!(p::AgParameter, sample::Union{Dict, NamedTuple})::Nothing
    p_name::Symbol = Symbol(p.name)
    if p_name in keys(sample)
        setfield!(p, :value, sample)
    end

    return
end

function update_model!(comp, sample)::Nothing
    set_params!(comp, sample)

    return
end

"""
Create relational mapping between components and values.
"""
function component_relationship(agparams::DataFrame)::Dict
    comp_sep = Agtor.component_sep
    name_sep = Agtor.component_name_sep
    field_sep = Agtor.component_field_sep

    relation = Dict()
    names = agparams.name

    for n in names
        nstr = String(n)
        components = split(nstr, comp_sep)
        top_level = components[1]

        if !haskey(relation, top_level)
            relation[top_level] = Dict()
        end

        upper = relation[top_level]

        for c in components[2:end]
            if contains(c, field_sep)
                comp, fld_name = Tuple(split(c, field_sep))

                field_values = agparams[in([n]).(names), [:default, :min_val, :max_val]]
                field_values = copy(field_values[1, :])

                if haskey(upper, comp)
                    upper[comp][fld_name] = field_values
                else
                    upper[comp] = Dict(fld_name=>field_values)
                end

                continue
            end

            if !haskey(upper, c)
                upper[c] = Dict()
            end

            upper = upper[c]
        end
    end

    return relation
end

ext: jl

max_stars_repo_name: ConnectedSystems/agtor-jl

max_issues_repo_name: ConnectedSystems/agtor-jl

max_forks_repo_name: ConnectedSystems/agtor-jl


text:
using CompositeConvex
using LinearAlgebra

τ0 = 2

@testset "Testing 1D norms" begin
    n = 100
    x = randn(n)
    y = similar(x)
    for p in [0,1,2]
        compute_norm_projection!(y, x, τ0, pNorm{p}())
        @test abs(norm(y, p) - τ0) < 1e-2
    end
end

@testset "Testing 2D Norms" begin
    m,n = 100,100
    x = randn(m,n)
    y = similar(x)
    row_norms = t->reshape(sum(t.^2,dims=2).^(0.5),size(t,2))
    col_norms = t->reshape(sum(t.^2,dims=1).^(0.5),size(t,1))

    compute_norm_projection!(y,x,τ0, pqNorm{0,2}())
    @test abs(norm(col_norms(y),0) - τ0) < 1e-2

    compute_norm_projection!(y,x,τ0, pqNorm{1,2}())
    @test abs(norm(col_norms(y),1) - τ0) < 1e-2

    compute_norm_projection!(y,x,τ0, pqNorm{2,0}())
    @test abs(norm(row_norms(y),0) - τ0) < 1e-2

    compute_norm_projection!(y,x,τ0, pqNorm{2,1}())
    @test abs(norm(row_norms(y),1) - τ0) < 1e-2
end


text:
using Test
using MLE
using Distributions
using Random
using Optim
using RCall

@testset "Test MLE" begin
    Random.seed!(1234)
    data = rand(Gumbel(1, 1), 1_000)

    @testset "Test Normal" begin
	f = fit_numerical_mle(Normal, data)
	@test f.converged

	f = fit_numerical_mle(Normal, data, (μ=10, σ=2))
	@test f.converged

	f = fit_numerical_mle(Normal, data, (μ=10, σ=2), bounds=(μ=[5,Inf], ))
	@test f.converged
    end

    f = fit_numerical_mle(Gumbel, data)
    @test f.converged
    @test isapprox(f.result, [0.9843, 0.9863]; atol=1e-3)

    f = fit_numerical_mle(Gumbel, data, holding=(; θ=1))
    @test f.converged

    pois_data = [1,2,1,3]

    f = fit_numerical_mle(Poisson, pois_data)
    @test f.converged
    @test isapprox(f.result, [1.75])

    R"""
    library(MASS)
    r_mle <- fitdistr($pois_data, "poisson")
    """
    @rget r_mle
    @test isapprox(f.result[1], r_mle[:estimate])
    @test isapprox(f.sd[1], r_mle[:sd])
    @test isapprox(f.loglikelihood, r_mle[:loglik])

    f = fit_numerical_mle(NegativeBinomial, pois_data)
    @test f.converged
    @test isapprox(mean(f.fit), 1.75; atol=1e-4)

    R"""
    library(MASS)
    r_mle <- fitdistr($pois_data, "negative binomial")
    """
    @rget r_mle
    @test isapprox(mean(f.fit), r_mle[:estimate][2]; atol=1e-4)

    f = fit_numerical_mle(NegativeBinomial, pois_data, holding=(;r=0.5))
    @test f.converged
    @test isapprox(mean(f.fit), 1.75; atol=1e-4)

    binom_data = [0, 1, 2, 0, 1, 2]
    f = fit_numerical_mle(Binomial, binom_data, holding=(;n=4))
    @test f.converged
    @test isapprox(f.result[1], 0.25; atol=1e-4)
    @test isapprox(f.result[1], mean(binom_data) / 4)

    f = fit_numerical_mle(Geometric, pois_data)
    @test f.converged
    @test isapprox(mean(f.fit), mean(pois_data); atol=1e-4)

    beta_data = [0.1, 0.2, 0.1, 0.2, 0.9]
    f = fit_numerical_mle(Beta, beta_data)
    @test f.converged
    @test isapprox(mean(f.fit), 0.35409; atol=1e-4)

    R"""
    library(MASS)
    r_mle <- fitdistr($beta_data, "beta", start=list(shape1=1, shape2=1))
    """
    @rget r_mle
    @test isapprox(f.result, r_mle[:estimate]; atol=1e-4)
    @test isapprox(f.vcov, r_mle[:vcov]; atol=1e-4)
    @test isapprox(f.sd, r_mle[:sd]; atol=1e-4)
    @test isapprox(f.loglikelihood, r_mle[:loglik])

    expon_data = [0.1, 0.2, 1.0, .3, .05]
    f = fit_numerical_mle(Exponential, expon_data)
    @test f.converged
    @test isapprox(mean(f.fit), mean(expon_data); atol=1e-2)

    R"""
    library(MASS)
    r_mle <- fitdistr($expon_data, "exponential")
    """
    @rget r_mle
    @test isapprox(1 ./ f.result[1], r_mle[:estimate]; atol=1e-2)
    @test isapprox(f.loglikelihood, r_mle[:loglik])

    f = fit_numerical_mle(Weibull, expon_data)
    @test f.converged

    R"""
    library(MASS)
    r_mle <- fitdistr($expon_data, "weibull")
    """
    @rget r_mle
    @test isapprox(f.result, r_mle[:estimate]; atol=1e-4)
    @test isapprox(f.sd, r_mle[:sd]; atol=1e-4)
    @test isapprox(f.loglikelihood, r_mle[:loglik])
end


text:
using StructuralCausalModels, Test

ProjDir = @__DIR__

a= [
  0    0    0    0;
  1    0    0    1;
  0    1    0    0;
  0    0    0    0
]

vars = [Symbol("n$i") for i in 1:size(a, 1)]
amat = NamedArray(Int.(a), (vars, vars), ("Rows", "Cols"));

dag = DAG("m_separation", amat)

m = maximize(ribbon_graph(amat, c=[:n1], m=[:n4]))

f = [:n1]
s = [:n4]
c = [:n3]

ms1 = m_separation(dag, f, s; debug=false)
ms2 = m_separation(dag, f, s; c=c)

@testset "m_separation" begin

  @test ms1 == true
  @test ms2 == false

end


text:
"""
    s2c(theta, phi)

Convert spherical coordinates to Cartesian.

### Given ###

- `theta`: Longitude angle (radians)
- `phi`: Latitude angle (radians)

### Returned ###

- `c`: Direction cosines

"""
function s2c(theta, phi)
    c = zeros(3)
    ccall((:eraS2c, liberfa), Cvoid,
          (Cdouble, Cdouble, Ptr{Cdouble}),
          theta, phi, c)
    c
end

"""
    s2p(theta, phi, r)

Convert spherical polar coordinates to p-vector.

### Given ###

- `theta`: Longitude angle (radians)
- `phi`: Latitude angle (radians)
- `r`: Radial distance

### Returned ###

- `p`: Cartesian coordinates

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraSxp`: multiply p-vector by scalar

"""
function s2p(theta, phi, r)
    p = zeros(3)
    ccall((:eraS2p, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
          theta, phi, r, p)
    p
end

"""
    s2pv(theta, phi, r, td, pd, rd)

Convert position/velocity from spherical to Cartesian coordinates.

### Given ###

- `theta`: Longitude angle (radians)
- `phi`: Latitude angle (radians)
- `r`: Radial distance
- `td`: Rate of change of theta
- `pd`: Rate of change of phi
- `rd`: Rate of change of r

### Returned ###

- `pv`: Pv-vector

"""
function s2pv(theta, phi, r, td, pd, rd)
    pv = zeros((2, 3))
    ccall((:eraS2pv, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
          theta, phi, r, td, pd, rd, pv)
    pv
end

"""
    s2xpv(s1, s2, pv)

Multiply a pv-vector by two scalars.

### Given ###

- `s1`: Scalar to multiply position component by
- `s2`: Scalar to multiply velocity component by
- `pv`: Pv-vector

### Returned ###

- `spv`: Pv-vector: p scaled by s1, v scaled by s2

### Note ###

   It is permissible for pv and spv to be the same array.

### Called ###

- `eraSxp`: multiply p-vector by scalar

"""
function s2xpv(s1, s2, pv)
    spv = zeros((2, 3))
    ccall((:eraS2xpv, liberfa), Cvoid,
          (Cdouble, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
          s1, s2, pv, spv)
    spv
end

"""
    starpm(ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b)

Star proper motion:  update star catalog data for space motion.

### Given ###

- `ra1`: Right ascension (radians), before
- `dec1`: Declination (radians), before
- `pmr1`: RA proper motion (radians/year), before
- `pmd1`: Dec proper motion (radians/year), before
- `px1`: Parallax (arcseconds), before
- `rv1`: Radial velocity (km/s, +ve = receding), before
- `ep1a`: "before" epoch, part A (Note 1)
- `ep1b`: "before" epoch, part B (Note 1)
- `ep2a`: "after" epoch, part A (Note 1)
- `ep2b`: "after" epoch, part B (Note 1)

### Returned ###

- `ra2`: Right ascension (radians), after
- `dec2`: Declination (radians), after
- `pmr2`: RA proper motion (radians/year), after
- `pmd2`: Dec proper motion (radians/year), after
- `px2`: Parallax (arcseconds), after
- `rv2`: Radial velocity (km/s, +ve = receding), after

### Notes ###

1. The starting and ending TDB dates ep1a+ep1b and ep2a+ep2b are
   Julian Dates, apportioned in any convenient way between the two
   parts (A and B).  For example, JD(TDB)=2450123.7 could be
   expressed in any of these ways, among others:

           epna          epnb

       2450123.7           0.0       (JD method)
       2451545.0       -1421.3       (J2000 method)
       2400000.5       50123.2       (MJD method)
       2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. In accordance with normal star-catalog conventions, the object's
   right ascension and declination are freed from the effects of
   secular aberration.  The frame, which is aligned to the catalog
   equator and equinox, is Lorentzian and centered on the SSB.

   The proper motions are the rate of change of the right ascension
   and declination at the catalog epoch and are in radians per TDB
   Julian year.

   The parallax and radial velocity are in the same frame.

3. Care is needed with units.  The star coordinates are in radians
   and the proper motions in radians per Julian year, but the
   parallax is in arcseconds.

4. The RA proper motion is in terms of coordinate angle, not true
   angle.  If the catalog uses arcseconds for both RA and Dec proper
   motions, the RA proper motion will need to be divided by cos(Dec)
   before use.

5. Straight-line motion at constant speed, in the inertial frame,
   is assumed.

6. An extremely small (or zero or negative) parallax is interpreted
   to mean that the object is on the "celestial sphere", the radius
   of which is an arbitrary (large) value (see the eraStarpv
   function for the value used).  When the distance is overridden in
   this way, the status, initially zero, has 1 added to it.

7. If the space velocity is a significant fraction of c (see the
   constant VMAX in the function eraStarpv), it is arbitrarily set
   to zero.  When this action occurs, 2 is added to the status.

8. The relativistic adjustment carried out in the eraStarpv function
   involves an iterative calculation.  If the process fails to
   converge within a set number of iterations, 4 is added to the
   status.

### Called ###

- `eraStarpv`: star catalog data to space motion pv-vector
- `eraPvu`: update a pv-vector
- `eraPdp`: scalar product of two p-vectors
- `eraPvstar`: space motion pv-vector to star catalog data

"""
function starpm(ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b)
    ra2 = Ref(0.0)
    dec2 = Ref(0.0)
    pmr2 = Ref(0.0)
    pmd2 = Ref(0.0)
    px2 = Ref(0.0)
    rv2 = Ref(0.0)
    i = ccall((:eraStarpm, liberfa), Cint,
              (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble},
              Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
              ra1, dec1, pmr1, pmd1, px1, rv1, ep1a, ep1b, ep2a, ep2b, ra2, dec2, pmr2, pmd2, px2, rv2)
    if i == -1
        throw(ERFAException("system error"))
    elseif i == 1
        @warn "distance overridden"
        return ra2[], dec2[], pmr2[], pmd2[], px2[], rv2[]
    elseif i == 2
        @warn "excessive velocity"
        return ra2[], dec2[], pmr2[], pmd2[], px2[], rv2[]
    elseif i == 4
        throw(ERFAException("solution didn't converge"))
    end
    ra2[], dec2[], pmr2[], pmd2[], px2[], rv2[]
end

"""
    starpv(ra, dec, pmr, pmd, px, rv)

Convert star catalog coordinates to position+velocity vector.

### Given (Note 1) ###

- `ra`: Right ascension (radians)
- `dec`: Declination (radians)
- `pmr`: RA proper motion (radians/year)
- `pmd`: Dec proper motion (radians/year)
- `px`: Parallax (arcseconds)
- `rv`: Radial velocity (km/s, positive = receding)

### Returned (Note 2) ###

- `pv`: pv-vector (au, au/day)

### Notes ###

1. The star data accepted by this function are "observables" for an
   imaginary observer at the solar-system barycenter.  Proper motion
   and radial velocity are, strictly, in terms of barycentric
   coordinate time, TCB.  For most practical applications, it is
   permissible to neglect the distinction between TCB and ordinary
   "proper" time on Earth (TT/TAI).  The result will, as a rule, be
   limited by the intrinsic accuracy of the proper-motion and
   radial-velocity data;  moreover, the pv-vector is likely to be
   merely an intermediate result, so that a change of time unit
   would cancel out overall.

   In accordance with normal star-catalog conventions, the object's
   right ascension and declination are freed from the effects of
   secular aberration.  The frame, which is aligned to the catalog
   equator and equinox, is Lorentzian and centered on the SSB.

2. The resulting position and velocity pv-vector is with respect to
   the same frame and, like the catalog coordinates, is freed from
   the effects of secular aberration.  Should the "coordinate
   direction", where the object was located at the catalog epoch, be
   required, it may be obtained by calculating the magnitude of the
   position vector pv[0][0-2] dividing by the speed of light in
   au/day to give the light-time, and then multiplying the space
   velocity pv[1][0-2] by this light-time and adding the result to
   pv[0][0-2].

   Summarizing, the pv-vector returned is for most stars almost
   identical to the result of applying the standard geometrical
   "space motion" transformation.  The differences, which are the
   subject of the Stumpff paper referenced below, are:

   (i) In stars with significant radial velocity and proper motion,
   the constantly changing light-time distorts the apparent proper
   motion.  Note that this is a classical, not a relativistic,
   effect.

   (ii) The transformation complies with special relativity.

3. Care is needed with units.  The star coordinates are in radians
   and the proper motions in radians per Julian year, but the
   parallax is in arcseconds; the radial velocity is in km/s, but
   the pv-vector result is in au and au/day.

4. The RA proper motion is in terms of coordinate angle, not true
   angle.  If the catalog uses arcseconds for both RA and Dec proper
   motions, the RA proper motion will need to be divided by cos(Dec)
   before use.

5. Straight-line motion at constant speed, in the inertial frame,
   is assumed.

6. An extremely small (or zero or negative) parallax is interpreted
   to mean that the object is on the "celestial sphere", the radius
   of which is an arbitrary (large) value (see the constant PXMIN).
   When the distance is overridden in this way, the status,
   initially zero, has 1 added to it.

7. If the space velocity is a significant fraction of c (see the
   constant VMAX), it is arbitrarily set to zero.  When this action
   occurs, 2 is added to the status.

8. The relativistic adjustment involves an iterative calculation.
   If the process fails to converge within a set number (IMAX) of
   iterations, 4 is added to the status.

9. The inverse transformation is performed by the function
   eraPvstar.

### Called ###

- `eraS2pv`: spherical coordinates to pv-vector
- `eraPm`: modulus of p-vector
- `eraZp`: zero p-vector
- `eraPn`: decompose p-vector into modulus and direction
- `eraPdp`: scalar product of two p-vectors
- `eraSxp`: multiply p-vector by scalar
- `eraPmp`: p-vector minus p-vector
- `eraPpp`: p-vector plus p-vector

### Reference ###

- Stumpff, P., 1985, Astron.Astrophys. 144, 232-240.

"""
function starpv(ra, dec, pmr, pmd, px, rv)
    pv = zeros((2, 3))
    i = ccall((:eraStarpv, liberfa), Cint,
              (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ptr{Cdouble}),
              ra, dec, pmr, pmd, px, rv, pv)
    if i == 1
        @warn "distance overridden"
        return pv
    elseif i == 2
        @warn "excessive speed "
        return pv
    elseif i == 4
        throw(ERFAException("solution didn't converge"))
    end
    pv
end

"""
    sxp(s, p)

Multiply a p-vector by a scalar.

### Given ###

- `s`: Scalar
- `p`: P-vector

### Returned ###

- `sp`: S * p

### Note ###

   It is permissible for p and sp to be the same array.

"""
function sxp(s, p)
    sp = zeros(3)
    ccall((:eraSxp, liberfa), Cvoid,
          (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
          s, p, sp)
    sp
end

"""
    sxpv(s, pv)

Multiply a pv-vector by a scalar.

### Given ###

- `s`: Scalar
- `pv`: Pv-vector

### Returned ###

- `spv`: S * pv

### Note ###

   It is permissible for pv and spv to be the same array

### Called ###

- `eraS2xpv`: multiply pv-vector by two scalars

"""
function sxpv(s, pv)
    spv = zeros((2, 3))
    ccall((:eraSxpv, liberfa), Cvoid,
          (Cdouble, Ptr{Cdouble}, Ptr{Cdouble}),
          s, pv, spv)
    spv
end

"""
    seps(al, ap, bl, bp)

Angular separation between two sets of spherical coordinates.

### Given ###

- `al`: First longitude (radians)
- `ap`: First latitude (radians)
- `bl`: Second longitude (radians)
- `bp`: Second latitude (radians)

### Returned ###

- Angular separation (radians)

### Called ###

- `eraS2c`: spherical coordinates to unit vector
- `eraSepp`: angular separation between two p-vectors

"""
function seps(al, ap, bl, bp)
    ccall((:eraSeps, liberfa), Cdouble, (Cdouble, Cdouble, Cdouble, Cdouble), al, ap, bl, bp)
end

"""
    sepp(a, b)

Angular separation between two p-vectors.

### Given ###

- `a`: First p-vector (not necessarily unit length)
- `b`: Second p-vector (not necessarily unit length)

### Returned ###

- Angular separation (radians, always positive)

### Notes ###

1. If either vector is null, a zero result is returned.

2. The angular separation is most simply formulated in terms of
   scalar product.  However, this gives poor accuracy for angles
   near zero and pi.  The present algorithm uses both cross product
   and dot product, to deliver full accuracy whatever the size of
   the angle.

### Called ###

- `eraPxp`: vector product of two p-vectors
- `eraPm`: modulus of p-vector
- `eraPdp`: scalar product of two p-vectors

"""
function sepp(a, b)
    ccall((:eraSepp, liberfa), Cdouble, (Ptr{Cdouble}, Ptr{Cdouble}), a, b)
end

"""
    s00a(date1, date2)

The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, using the IAU 2000A
precession-nutation model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- The CIO locator s in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems.  The two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The CIO locator s remains a small fraction of
   1 arcsecond throughout 1900-2100.

3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series
   is more compact than a direct series for s would be.  The present
   function uses the full IAU 2000A nutation model when predicting
   the CIP position.  Faster results, with no significant loss of
   accuracy, can be obtained via the function eraS00b, which uses
   instead the IAU 2000B truncated model.

### Called ###

- `eraPnm00a`: classical NPB matrix, IAU 2000A
- `eraBnp2xy`: extract CIP X,Y from the BPN matrix
- `eraS00`: the CIO locator s, given X,Y, IAU 2000A

### References ###

- Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

- n.b. The celestial ephemeris origin (CEO) was renamed "celestial
    intermediate origin" (CIO) by IAU 2006 Resolution 2.

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
s00a

"""
    s00b(date1, date2)

The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, using the IAU 2000B
precession-nutation model.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- The CIO locator s in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems.  The two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The CIO locator s remains a small fraction of
   1 arcsecond throughout 1900-2100.

3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series
   is more compact than a direct series for s would be.  The present
   function uses the IAU 2000B truncated nutation model when
   predicting the CIP position.  The function eraS00a uses instead
   the full IAU 2000A model, but with no significant increase in
   accuracy and at some cost in speed.

### Called ###

- `eraPnm00b`: classical NPB matrix, IAU 2000B
- `eraBnp2xy`: extract CIP X,Y from the BPN matrix
- `eraS00`: the CIO locator s, given X,Y, IAU 2000A

### References ###

- Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

- n.b. The celestial ephemeris origin (CEO) was renamed "celestial
    intermediate origin" (CIO) by IAU 2006 Resolution 2.

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
s00b

"""
    s06a(date1, date2)

The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, using the IAU 2006
precession and IAU 2000A nutation models.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- The CIO locator s in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems.  The two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The CIO locator s remains a small fraction of
   1 arcsecond throughout 1900-2100.

3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series is
   more compact than a direct series for s would be.  The present
   function uses the full IAU 2000A nutation model when predicting
   the CIP position.

### Called ###

- `eraPnm06a`: classical NPB matrix, IAU 2006/2000A
- `eraBpn2xy`: extract CIP X,Y coordinates from NPB matrix
- `eraS06`: the CIO locator s, given X,Y, IAU 2006

### References ###

- Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

- n.b. The celestial ephemeris origin (CEO) was renamed "celestial
    intermediate origin" (CIO) by IAU 2006 Resolution 2.

- Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

- McCarthy, D. D., Petit, G. (eds.), 2004, IERS Conventions (2003),
    IERS Technical Note No. 32, BKG

- Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

"""
s06a

"""
    sp00(date1, date2)

The TIO locator s', positioning the Terrestrial Intermediate Origin
on the equator of the Celestial Intermediate Pole.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)

### Returned ###

- The TIO locator s' in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The TIO locator s' is obtained from polar motion observations by
   numerical integration, and so is in essence unpredictable.
   However, it is dominated by a secular drift of about
   47 microarcseconds per century, which is the approximation
   evaluated by the present function.

### Reference ###

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
sp00

for name in ("s00a",
             "s00b",
             "s06a",
             "sp00")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d1, d2) = ccall(($fc, liberfa), Cdouble, (Cdouble, Cdouble), d1, d2)
end

"""
    s00(date1, date2, x, y)

The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, given the CIP's X,Y
coordinates.  Compatible with IAU 2000A precession-nutation.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)
- `x`, `y`: CIP coordinates (Note 3)

### Returned ###

- The CIO locator s in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems:  the two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The quantity s remains below 0.1 arcsecond
   throughout 1900-2100.

3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series
   is more compact than a direct series for s would be.  This
   function requires X,Y to be supplied by the caller, who is
   responsible for providing values that are consistent with the
   supplied date.

4. The model is consistent with the IAU 2000A precession-nutation.

### Called ###

- `eraFal03`: mean anomaly of the Moon
- `eraFalp03`: mean anomaly of the Sun
- `eraFaf03`: mean argument of the latitude of the Moon
- `eraFad03`: mean elongation of the Moon from the Sun
- `eraFaom03`: mean longitude of the Moon's ascending node
- `eraFave03`: mean longitude of Venus
- `eraFae03`: mean longitude of Earth
- `eraFapa03`: general accumulated precession in longitude

### References ###

- Capitaine, N., Chapront, J., Lambert, S. and Wallace, P.,
    "Expressions for the Celestial Intermediate Pole and Celestial
    Ephemeris Origin consistent with the IAU 2000A precession-
    nutation model", Astron.Astrophys. 400, 1145-1154 (2003)

- n.b. The celestial ephemeris origin (CEO) was renamed "celestial
    intermediate origin" (CIO) by IAU 2006 Resolution 2.

- McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

"""
s00

"""
    s06(date1, date2, x, y)

The CIO locator s, positioning the Celestial Intermediate Origin on
the equator of the Celestial Intermediate Pole, given the CIP's X,Y
coordinates.  Compatible with IAU 2006/2000A precession-nutation.

### Given ###

- `date1`, `date2`: TT as a 2-part Julian Date (Note 1)
- `x`, `y`: CIP coordinates (Note 3)

### Returned ###

- The CIO locator s in radians (Note 2)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` |                    |
   |:----------|:------------|:-------------------|
   | 2450123.7 |         0.0 | JD method          |
   | 2451545.0 |     -1421.3 | J2000 method       |
   | 2400000.5 |     50123.2 | MJD method         |
   | 2450123.5 |         0.2 | date & time method |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The CIO locator s is the difference between the right ascensions
   of the same point in two systems:  the two systems are the GCRS
   and the CIP,CIO, and the point is the ascending node of the
   CIP equator.  The quantity s remains below 0.1 arcsecond
   throughout 1900-2100.

3. The series used to compute s is in fact for s+XY/2, where X and Y
   are the x and y components of the CIP unit vector;  this series
   is more compact than a direct series for s would be.  This
   function requires X,Y to be supplied by the caller, who is
   responsible for providing values that are consistent with the
   supplied date.

4. The model is consistent with the "P03" precession (Capitaine et
   al. 2003), adopted by IAU 2006 Resolution 1, 2006, and the
   IAU 2000A nutation (with P03 adjustments).

### Called ###

- `eraFal03`: mean anomaly of the Moon
- `eraFalp03`: mean anomaly of the Sun
- `eraFaf03`: mean argument of the latitude of the Moon
- `eraFad03`: mean elongation of the Moon from the Sun
- `eraFaom03`: mean longitude of the Moon's ascending node
- `eraFave03`: mean longitude of Venus
- `eraFae03`: mean longitude of Earth
- `eraFapa03`: general accumulated precession in longitude

### References ###

- Capitaine, N., Wallace, P.T. & Chapront, J., 2003, Astron.
    Astrophys. 432, 355

- McCarthy, D.D., Petit, G. (eds.) 2004, IERS Conventions (2003),
    IERS Technical Note No. 32, BKG

"""
s06

for name in ("s00",
             "s06")
    f = Symbol(name)
    fc = "era" * uppercasefirst(name)
    @eval ($f)(d1, d2, t1, t2) = ccall(($fc, liberfa), Cdouble, (Cdouble, Cdouble, Cdouble, Cdouble), d1, d2, t1, t2)
end


text:
function operator_fn(c :: Char) :: Symbol
    operators = Dict(
        'c' => :change,
        'y' => :yank,
        'd' => :delete
    )
    return operators[c]
end

function change(buf :: IOBuffer, motion :: Motion) #, motion_type :: MotionType)
    delete(buf, motion)
end

# also, when there is whitespace following a word,
# "dw" deletes that whitespace
# and "cw" only removes the inner word
function delete(buf :: IOBuffer, motion :: Motion) #, motion_type :: MotionType)
    yank(buf, motion)
    move(buf, motion) #, motion_type)
    @log motion.stop
    @log edit_splice!(buf, motion.start => motion.stop)
    return nothing
end

# function yank(buf :: IOBuffer, motion :: Motion) #, motion_type :: MotionType)
#     yank(buf, motion)

#     # '0' is the "yank" register
#     vim.registers['0'] = cut(buf, motion)
# end
function move(buf :: IOBuffer, motion :: Motion)#, motion_type :: MotionType)
    seek(buf, motion.stop)
end

function yank(buf :: IOBuffer, motion :: Motion) :: Union{String, Nothing}
    text = cut(buf, motion)
    reg = vim.register
    # '_' is the "black hole" register
    vim.register = '"'
    if reg == '_'
        return nothing
    else
        vim.registers[reg] = text
        return text
    end
end

function insert(buf :: IOBuffer, pos :: Integer, text)
    s = string(text)
    edit_splice!(buf, pos => pos, s)
end

function cut(buf :: IOBuffer, motion :: Motion) :: String
    mark(buf)
    seek(buf, min(motion))
    chars = [ read(buf, Char) for i in 1:length(motion) ]
    reset(buf)
    return String(chars)
end

function paste(buf :: IOBuffer, reg :: Char)
    if !(reg in vim.registers.keys)
        return
    end
    pos = position(buf)
    edit_splice!(buf, pos => pos, vim.registers[reg])
end


text:
julia > function mycompare(a, b)::Cint
      (a < b) ? -1 : ((a > b) ? +1 : 0)
  end
mycompare (generic function with 1 method)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
T = Float64

# options
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,live_plotting=:off)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=false,
    opts_uncon=opts_ilqr,
    cost_tolerance=1.0e-4,
    cost_tolerance_intermediate=1.0e-3,
    constraint_tolerance=max_con_viol,
    penalty_scaling=50.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-3);

opts_ipopt = DIRCOLSolverOptions{T}(verbose=verbose,
    nlp=Ipopt.Optimizer(),
    feasibility_tolerance=max_con_viol)

# AL-iLQR
prob = copy(Problems.cartpole_static)
U0 = deepcopy(controls(prob))
alilqr = AugmentedLagrangianSolver(prob, opts_al)
solve!(alilqr)
max_violation(alilqr)

@btime begin
    initial_controls!($alilqr, $U0)
    solve!($alilqr)
end

# ALTRO
altro = ALTROSolver(prob, opts_altro)
initial_controls!(altro, U0)
solve!(altro)
max_violation(altro)

@btime begin
    initial_controls!($altro, $U0)
    solve!($altro)
end

# Ipopt
prob_ipopt = copy(Problems.cartpole_static)
prob_ipopt = TO.change_integration(prob_ipopt, HermiteSimpson)
rollout!(prob_ipopt)
ipopt = DIRCOLSolver(prob_ipopt, opts_ipopt)
ipopt.opts.verbose = false
solve!(ipopt)
max_violation(ipopt)

@btime solve!($ipopt)


text:
module Pitches

import Base: show, +, -, *, convert, zero, isless, isequal, isapprox

export Interval, IntervalClass, Pitch
export octave, unison, ic, isstep, chromsemi, embed
export intervaltype, intervalclasstype
export topitch, tointerval, pc

# Pitch: basic Types and Interfaces
# =================================

"""
   abstract type Interval end

Any interval type should be a subtype of `Interval`.
Intervals should implement the following operations as far as possible:
- `ic`
- `octave(T)`
- `isstep`
- `chromsemi(T)`
- `intervalclasstype(T)`
- `Base.+`
- `Base.-` (negation and substraction)
- `Base.*` (with integers, both sides)
- `Base.zero(T)`
- `Base.sign`
Where `(T)` marks operations on the type itself.
"""
abstract type Interval end

"""
    abstract type IntervalClass <: Interval end

Any interval class type should be a subtype of `IntervalClass`.
In addition to the methods on intervals, interval classes should implement:
- `embed`
- `intervaltype(T)`
`intervalclasstype(T)` and `ic` should be identities.
"""
abstract type IntervalClass <: Interval end

# interfaces
# ----------

"""
    unison(T)

Returns the interval of a unison for interval type `T`.
Alias for `Base.zero`.
New interval types should implement `Base.zero`,
but user code should call `unison` for better interpretability.
"""
const unison = Base.zero

"""
    octave(T, [n=1])

Returns the interval corresponding to an octave for interval type `T`.
For interval classes, this should return `zero(T)`
(a default method is provided).

If `n` is specified, the octave is multiplied by `n` first.
This is equivalent to `octave(T) * n`.

For convenience, a fallback for `octave(p::T, [n])` is provided.
Only `octave(T)` needs to be implemented.
"""
function octave end
octave(T::Type{PC}) where {PC<:IntervalClass} = zero(T)
octave(T, n::Int) = octave(T) * n
octave(p::Interval) = octave(typeof(p))

"""
    ic(i)

Returns the interval class of an interval, removing the octave
"""
function ic end

"""
    embed(ic, [oct=0])
    embed(pc, [oct=0])

Converts an interval class to an interval in the canonical octave,
adding `oct` octaves, if supplied.
Also works for pitches.
"""
function embed end
embed(ic, oct) = embed(ic) + octave(intervaltype(typeof(ic)), oct)

"""
    intervaltype(IC::Type)

Returns for an interval class type `IC` the corresponding interval type.
For convenience, `intervaltype(ic::IC)` is also provided.
"""
function intervaltype end
intervaltype(::Any) = nothing
intervaltype(::IC) where {IC<:IntervalClass} = intervaltype(IC)

"""
    intervalclasstype(I::Type)

Returns for an interval type `I` the corresponding interval class type.
For convenience, `intervalclasstype(p::P)` is also provided.
"""
function intervalclasstype end
intervalclasstype(::Any) = nothing
intervalclasstype(::I) where {I<:Interval} = intervalclasstype(I)

"""
    isstep(p)

For diatonic intervals, indicates whether `p` is a step.
"""
function isstep end

"""
    chromsemi(I::Type)

Returns a chromatic semitone of type `I`.
"""
function chromsemi end

# pitches
# =======

"""
    Pitch{I}

Represents a pitch for the interval type `I`.
The interval is interpreted as an absolute pitch
by assuming a reference pitch.
The reference pitch is type dependent and known from context.
"""
struct Pitch{I<:Interval}
    pitch :: I
end

topitch(i::I) where {I<:Interval} = Pitch(i)

tointerval(p::Pitch{I}) where {I<:Interval} = p.pitch

+(p::Pitch{I}, i::I) where {I<:Interval} = Pitch(p.pitch + i)
+(i::I, p::Pitch{I}) where {I<:Interval} = Pitch(p.pitch + i)
-(p::Pitch{I}, i::I) where {I<:Interval} = Pitch(p.pitch - i)
-(p1::Pitch{I}, p2::Pitch{I}) where {I<:Interval} = p1.pitch - p2.pitch

"""
    pc(p)

Return the pitch class that corresponds to `p`.
"""
pc(p::Pitch{I}) where {I<:Interval} = Pitch(ic(p.pitch))
embed(p::Pitch{I}, octs::Int=0) where {I<:Interval} = Pitch(embed(p.pitch, octs))

Base.isequal(p1::Pitch{I}, p2::Pitch{I}) where {I<:Interval} = p1.pitch == p2.pitch
Base.isapprox(p1::Pitch, p2::Pitch; kwargs...) =
    Base.isapprox(p1.pitch, p2.pitch; kwargs...)
Base.isless(p1::Pitch{I}, p2::Pitch{I}) where {I<:Interval} = p1.pitch < p2.pitch
Base.hash(p::Pitch, x::UInt) = hash(p.pitch, x)

# specific interval types
# =======================

include("pitches/midi.jl")
include("pitches/spelled.jl")
include("pitches/logfreq.jl")

end # module


text:
macro afCall_Arr_Arr(method, holder, func)
    quote
        function $(esc(method)){B}(arr::AFArray{B})
            verifyAccess(arr)
            af = arr.af
            result = af.results.ptr
            err = ccall(af.$holder.$func,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}),
                result, arr.ptr)
            assertErr(err)
            array(af, result[])
        end
    end
end

macro afCall_Arr_Arr_Arr_Unsigned(method, holder, func)
    quote
        function $(esc(method)){B}(arr1::AFArray{B}, arr2::AFArray{B}, value)
            verifyAccess(arr1)
            verifyAccess(arr2)
            af = arr1.af
            result = af.results.ptr
            err = ccall(af.$holder.$func,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Void}, UInt32),
                result, arr1.ptr, arr2.ptr, UInt32(value))
            assertErr(err)
            array(af, result[])
        end
    end
end

macro afCall_Arr_Arr_Bool(method, holder, func, default)
    quote
        function $(esc(method))(arr::AFArray, value::Bool = $default)
            verifyAccess(arr)
            af = arr.af
            result = af.results.ptr
            err = ccall(af.$holder.$func,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Bool),
                result, arr.ptr, value)
            assertErr(err)
            array(af, result[])
        end
    end
end

macro afCall_Arr_Arr_DimT(method, holder, func, default)
    quote
        function $(esc(method))(arr::AFArray, value::DimT = $default)
            verifyAccess(arr)
            af = arr.af
            result = af.results.ptr
            if value == -1
                value = firstDim(af, arr.ptr)
            end
            err = ccall(af.$holder.$func,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, DimT),
                result, arr.ptr, value)
            assertErr(err)
            array(af, result[])
        end
    end
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# This file was generated, do not modify it. # hide
mach_iterated_pipe = machine(iterated_pipe, X, y)
fit!(mach_iterated_pipe);

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
    function profile!(name, results; kwargs...)

Compute the profiled likelihood over a single parameter of
interest given by `name`.

# Arguments
- `name::String`: parameter of interest
- `results::Results`: fit results
- `stop::Float64=5.0`: Maximum ΔLog(Likelihood)
- `step::Float64=1.0`: Incremental step of the parameter of interest
- `prior=nothing`: Prior probability of the parameter of interest
"""
function profile!(name, results; kwargs...)
  kwargs = Dict(kwargs)
  stop      = get(kwargs, :stop, 5.0)
  step      = get(kwargs, :step, 1.0)
  prior     = get(kwargs, :prior, nothing)
  verbose   = get(kwargs, :verbose, false)
  max_steps = get(kwargs, :maxsteps, 1000)

  fitterOptions = get(kwargs, :fitteroptions, nothing)
  if fitterOptions != nothing
    for (k, v) in fitterOptions
      setproperty!(results.opt, k, v)
    end
  end

  idx = 0
  for (i, p) in enumerate(results.model.params)
    if p.name == name
      idx = i
    end
  end
  if idx == 0
    @warn idx
    return 0, 0
  end
  ## Store these and report if verbose
  verb_count = 0
  verb_steps = Array{Float64}(undef, 0)
  #component = results.model.component_dict[name]
  param = results.model.params[idx]
  ## Go left
  x_left = []
  nll_left = []
  nlow = copy(results.lower_bounds)
  nhigh = copy(results.upper_bounds)
  p0 = copy(results.min_parameters)
  val = copy(results.min_objective)
  xeval = p0[idx]
  count = 0
  while (val-results.min_objective) < stop
    nlow[idx] = xeval
    nhigh[idx] = xeval
    results.opt.lower_bounds = nlow
    results.opt.upper_bounds = nhigh
    p0[idx] = xeval
    init_step = get(kwargs, :init_step, p0./100.0)
    init_step[idx] = xeval * 1e-6
    init_step[(init_step .== 0)] .= 0.1
    results.opt.initial_step=init_step
    val, minx, ret = optimize!(results.opt, p0)
    likelihood_value = exp(-(val-results.min_objective))
    push!(verb_steps, results.opt.numevals)
    if (likelihood_value == Inf) || (likelihood_value < 0)
      println("CYA", likelihood_value)
      break
    end
    push!(nll_left, likelihood_value )
    push!(x_left, xeval)
    xeval -= step
    count += 1
    if count > max_steps
      @warn "Profiling hit limit" max_steps
      break
    end
  end
  verb_count += count
  ## Go right
  x_right = []
  nll_right = []
  nlow = copy(results.lower_bounds)
  nhigh = copy(results.upper_bounds)
  p0 = copy(results.min_parameters)
  val = copy(results.min_objective)
  xeval = p0[idx]
  count = 0
  while (val-results.min_objective) < stop
    nlow[idx] = xeval
    nhigh[idx] = xeval
    results.opt.lower_bounds = nlow
    results.opt.upper_bounds = nhigh
    p0[idx] = xeval
    init_step = get(kwargs, :init_step, p0./100.0)
    init_step[idx] = xeval * 1e-6
    init_step[(init_step .== 0)] .= 0.1
    results.opt.initial_step=init_step
    val, minx, ret = optimize!(results.opt, p0)
    likelihood_value = exp(-(val-results.min_objective))
    push!(verb_steps, results.opt.numevals)
    if (likelihood_value == Inf) || (likelihood_value < 0)
      break
    end
    push!(nll_right, likelihood_value )
    push!(x_right, xeval)
    xeval += step
    count += 1
    if count > max_steps
      @warn "Profiling hit limit" max_steps
      break
    end
  end
  verb_count += count
  x = append!(reverse(x_left), x_right[2:end])
  nll = append!(reverse(nll_left), nll_right[2:end])
  ## update the component with the likelihood results
  if prior != nothing
    nll = nll .* prior.(x)
  end
  param.likelihood_x = x
  param.likelihood_y = nll
  if verbose
    @show verb_count
    @show verb_steps
  end
  return x, nll
end

function compute_profiled_uncertainties!(results; kwargs...)
  kwargs = Dict(kwargs)
  CL = get(kwargs, :CL, nothing)
  σ = get(kwargs, :σ, nothing)
  mode = get(kwargs, :mode, "FC")

  arraystep = get(kwargs, :arraystep, nothing)
  if CL != nothing
    α = CL
  elseif σ != nothing
    α = erf(σ/sqrt(2))
  else
    α = 0.90
  end
  for (idx,param) in enumerate(results.model.params)
    if arraystep != nothing
      profile!(param.name, results; step=arraystep[idx], kwargs...)
    else
      profile!(param.name, results; kwargs...)
    end
    uncertainty!(param.name, results; CL=α, mode=mode)
  end
end

function uncertainty!(name, results;
                      CL=nothing, σ=nothing, mode="FC", prior=nothing)
  if CL != nothing
    α = CL
  elseif σ != nothing
    α = erf(σ/sqrt(2))
  else
    α = erf(1/sqrt(2)) # 1 Sigma default
  end
  idx = 0
  for (i, p) in enumerate(results.model.params)
    if p.name == name
      idx = i
    end
  end
  param = results.model.params[idx]
  x, y = param.likelihood_x, param.likelihood_y
  # Apply the arbitrary prior function f(x)
  if prior != nothing
    y = y .* prior.(x)
  end
  # Mode can be FC (Feldman-Cousins), Mode-centered, mean-centered, left, right, etc.
  # Return cumulative, x_low, x_high
  y = y / sum(y)
  cumulative = cumsum(y)
  β = 1 - α
  if mode == "FC"
    cy = sortperm(y; rev=true)
    cum_y = cumsum(y[cy])
    cum_x = x[cy]
    region = cum_x[ cum_y .<= α ]
    left = minimum(region)
    right = maximum(region)
  end
  param.low, param.high = left, right
  left, right
end


text:
# q-SPC/Fw water model.

# Paesani, F., Zhang, W., Case, D. A., Cheatham III, T. E., & Voth, G. A.
# (2006). An accurate and simple quantum model for liquid water. The Journal of
# Chemical Physics, 125(18), 184507.

"""
    QSPCFW

Implementation of the q-SPC/Fw flexible water model.
"""
struct QSPCFW <: WaterModel
    "Potential terms (returning kJ / mol)."
    potential::TermContainer
    "Force terms (returning kJ / mol nm)."
    force::TermContainer
end

function QSPCFW(N::Int)
    build_model(QSPCFW,
                qspcfw_potential_intra, qspcfw_potential_inter,
                qspcfw_force_intra, qspcfw_force_inter,
                N)
end

const qspcfw_k_b = 1059.162 * kcal_to_kJ / A_to_nm^2 # kJ / mol nm^2
const qspcfw_r_OH_eq = 1.0 * A_to_nm # nm
const qspcfw_k_a = 75.90 * kcal_to_kJ # kJ / mol
const qspcfw_theta_HOH_eq = 112.0 * deg_to_rad

const qspcfw_epsilon_OO = 0.1554252 * kcal_to_kJ # kJ / mol
const qspcfw_sigma_OO = 3.165492 * A_to_nm # nm

const qspcfw_epsilon_0 = (8.8541878128e-12 / e_to_C^2 / N_A / to_kilo
                          / to_nano) # e^2 mol / kJ nm
const qspcfw_q_O = -0.84 # e
const qspcfw_q_H = 0.42 # e
const qspcfw_qs = Vector{Float64}(undef, length(MS))
qspcfw_qs[AO] = qspcfw_q_O
qspcfw_qs[AH1] = qspcfw_q_H
qspcfw_qs[AH2] = qspcfw_q_H

function qspcfw_potential_intra(n::Int)
    function term_f(j::Int, s::State)
        rvec_OH_1 = s.qs[j, :, AO, n] .- s.qs[j, :, AH1, n] # nm
        rvec_OH_2 = s.qs[j, :, AO, n] .- s.qs[j, :, AH2, n] # nm
        r_OH_1 = norm(rvec_OH_1) # nm
        r_OH_2 = norm(rvec_OH_2) # nm
        theta_HOH = acos_(dot(rvec_OH_1 ./ r_OH_1, rvec_OH_2 ./ r_OH_2))

        result = 0.0 # kJ / mol
        result += 0.5 * qspcfw_k_b * (r_OH_1 - qspcfw_r_OH_eq)^2
        result += 0.5 * qspcfw_k_b * (r_OH_2 - qspcfw_r_OH_eq)^2
        result += 0.5 * qspcfw_k_a * (theta_HOH - qspcfw_theta_HOH_eq)^2

        result
    end
end

function qspcfw_elec_pot(n1::Int, a1::Int, n2::Int, a2::Int, j::Int, s::State)
    R = norm(s.qs[j, :, a1, n1] .- s.qs[j, :, a2, n2]) # nm
    qspcfw_qs[a1] * qspcfw_qs[a2] / (4pi * qspcfw_epsilon_0 * R) # kJ / mol
end

function qspcfw_potential_inter(n1::Int, n2::Int)
    function term_f(j::Int, s::State)
        result = 0.0 # kJ / mol

        # O-O (LJ)
        R = norm(s.qs[j, :, AO, n1] .- s.qs[j, :, AO, n2]) # nm
        diff = (qspcfw_sigma_OO / R)^12 - (qspcfw_sigma_OO / R)^6
        result += 4 * qspcfw_epsilon_OO * diff

        # O-O (electrostatic)
        result += qspcfw_elec_pot(n1, AO, n2, AO, j, s)

        # O-H (electrostatic)
        result += qspcfw_elec_pot(n1, AO, n2, AH1, j, s)
        result += qspcfw_elec_pot(n1, AO, n2, AH2, j, s)
        result += qspcfw_elec_pot(n2, AO, n1, AH1, j, s)
        result += qspcfw_elec_pot(n2, AO, n1, AH2, j, s)

        # H-H (electrostatic)
        result += qspcfw_elec_pot(n1, AH1, n2, AH1, j, s)
        result += qspcfw_elec_pot(n1, AH1, n2, AH2, j, s)
        result += qspcfw_elec_pot(n1, AH2, n2, AH1, j, s)
        result += qspcfw_elec_pot(n1, AH2, n2, AH2, j, s)

        result
    end
end

function qspcfw_force_intra(n::Int)
    function term_f!(result::ForceContainer, j::Int, s::State)
        rvec_OH_1 = s.qs[j, :, AO, n] .- s.qs[j, :, AH1, n] # nm
        rvec_OH_2 = s.qs[j, :, AO, n] .- s.qs[j, :, AH2, n] # nm
        r_OH_1 = norm(rvec_OH_1) # nm
        r_OH_2 = norm(rvec_OH_2) # nm
        theta_HOH = acos_(dot(rvec_OH_1 ./ r_OH_1, rvec_OH_2 ./ r_OH_2))

        # kJ / nm mol
        F_OH_1 = -qspcfw_k_b * (1.0 - qspcfw_r_OH_eq / r_OH_1) .* rvec_OH_1
        result.f[j, :, AO, n] .+= F_OH_1
        result.f[j, :, AH1, n] .-= F_OH_1

        # kJ / nm mol
        F_OH_2 = -qspcfw_k_b * (1.0 - qspcfw_r_OH_eq / r_OH_2) .* rvec_OH_2
        result.f[j, :, AO, n] .+= F_OH_2
        result.f[j, :, AH2, n] .-= F_OH_2

        x = (dot(rvec_OH_1, rvec_OH_1) * dot(rvec_OH_2, rvec_OH_2)
             - dot(rvec_OH_1, rvec_OH_2)^2) # nm^4
        # kJ / nm^2 mol
        pre_F_HOH = qspcfw_k_a * (theta_HOH - qspcfw_theta_HOH_eq) / sqrt(x)
        x_HOH_1 = (dot(rvec_OH_1, rvec_OH_2) / dot(rvec_OH_1, rvec_OH_1)
                       .* rvec_OH_1
                   .- rvec_OH_2) # nm
        x_HOH_2 = (dot(rvec_OH_1, rvec_OH_2) / dot(rvec_OH_2, rvec_OH_2)
                       .* rvec_OH_2
                   .- rvec_OH_1) # nm
        @. result.f[j, :, AO, n] -= pre_F_HOH * (x_HOH_1 + x_HOH_2)
        @. result.f[j, :, AH1, n] += pre_F_HOH * x_HOH_1
        @. result.f[j, :, AH2, n] += pre_F_HOH * x_HOH_2

        nothing
    end
end

function qspcfw_elec_force!(result::ForceContainer, n1::Int, a1::Int,
                            n2::Int, a2::Int, j::Int, s::State)
    Rvec = s.qs[j, :, a1, n1] .- s.qs[j, :, a2, n2] # nm
    R = norm(Rvec) # nm
    # kJ / nm mol
    F = qspcfw_qs[a1] * qspcfw_qs[a2] / (4pi * qspcfw_epsilon_0 * R^3) .* Rvec
    result.f[j, :, a1, n1] .+= F
    result.f[j, :, a2, n2] .-= F

    nothing
end

function qspcfw_force_inter(n1, n2)
    function term_f!(result::ForceContainer, j::Int, s::State)
        # O-O (LJ)
        Rvec = s.qs[j, :, AO, n1] .- s.qs[j, :, AO, n2] # nm
        R = norm(Rvec) # nm
        # 1 / nm^2
        diff = 12 * qspcfw_sigma_OO^12/R^14 - 6 * qspcfw_sigma_OO^6/R^8
        F = 4 * qspcfw_epsilon_OO * diff .* Rvec # kJ / nm mol
        result.f[j, :, AO, n1] .+= F
        result.f[j, :, AO, n2] .-= F

        # O-O (electrostatic)
        qspcfw_elec_force!(result, n1, AO, n2, AO, j, s)

        # O-H (electrostatic)
        qspcfw_elec_force!(result, n1, AO, n2, AH1, j, s)
        qspcfw_elec_force!(result, n1, AO, n2, AH2, j, s)
        qspcfw_elec_force!(result, n2, AO, n1, AH1, j, s)
        qspcfw_elec_force!(result, n2, AO, n1, AH2, j, s)

        # H-H (electrostatic)
        qspcfw_elec_force!(result, n1, AH1, n2, AH1, j, s)
        qspcfw_elec_force!(result, n1, AH1, n2, AH2, j, s)
        qspcfw_elec_force!(result, n1, AH2, n2, AH1, j, s)
        qspcfw_elec_force!(result, n1, AH2, n2, AH2, j, s)

        nothing
    end
end

registered_models["q-SPC/Fw"] = QSPCFW


text:
using LowRankIntegrators, Plots
N = 100
Ws = []
for k in 1:2
    W = randn(N,N)
    for i in 1:N
        W[i,i] = 0
        for j in 1:i
            W[i,j] = - W[j,i]
        end
    end
    push!(Ws, W)
end

D = Diagonal([2.0^(-j) for j in 1:N])
Y(t) = exp(t*Ws[1])*exp(t)*D*exp(t*Ws[2])
Y0 = Y(0.0)
Yf = Y(1.0)

# initialization of low rank approximation based on the SVD
r_range = [4,8,16]
dt_range = [0.1,0.05,0.01,0.005]#,0.0001]
tol = 1e-8 # needed for rank adaptation
S_kwargs, K_kwargs, L_kwargs = Dict(), Dict(), Dict()
alg = Euler()
solvers = [ProjectorSplitting(PrimalLieTrotter()), ProjectorSplitting(DualLieTrotter()), ProjectorSplitting(Strang()), UnconventionalAlgorithm(), RankAdaptiveUnconventionalAlgorithm(tol, rmax = 20)]
error = Dict()
for r in r_range
    X0 = truncated_svd(Y0, r)
    prob = MatrixDataProblem(Y, X0, (0.0, 1.0))
    for dt in dt_range
        for solver in solvers
            sol = LowRankIntegrators.solve(prob, solver, dt)
            error[(r,dt,solver)] = norm(Matrix(sol.Y[end]) - Yf)
        end
    end
end

# visualization of data
solver_names = Dict(solvers[1] => "Primal Lie Trotter Projector Splitting",
                    solvers[2] => "Dual Lie Trotter Projector Splitting",
                    solvers[3] => "Strang Projector Splitting",
                    solvers[4] => "Unconventional Algorithm",
                    solvers[5] => "Rank Adaptive Unconventional Algorithm")
error_comparison = plot(yscale = :log, xscale = :log, xlabel = "step size", ylabel = "error", legend = :bottomleft)

colors = Dict(solvers[1] => :blue, solvers[2] => :red, solvers[3] => :green, solvers[4] => :dodgerblue, solvers[5] => :orange)
marker_range = [:v, :hex, :o, :star, :cross, :^]
markers = Dict(r_range[k] => marker_range[k] for k in 1:length(r_range))
for r in r_range
    for solver in solvers
        try
            plot!(error_comparison, dt_range, [error[(r, dt, solver)] for dt in dt_range], label = (r == r_range[1] ? solver_names[solver] : nothing), linestyle=:dash, marker=markers[r], color = colors[solver])
        catch
            continue
        end
    end
end
display(error_comparison)


text:
using JeszenszkiBasis
using Test

@testset "JeszenszkiBasis" begin

@testset "Szbasis" begin
    @testset for (K, N, D) in [(50, 0, 1), (1, 1, 1), (2, 3, 4), (3, 2, 6), (4, 4, 35)]
        sb = Szbasis(K, N)

        @test sb.K == K
        @test sb.N == N
        @test sb.D == D
        @test num_vectors(N, K) == D
        @test size(sb.vectors) == (sb.K, sb.D)
        @test length(sb) == D
        @test site_max(sb) == N

        # Occupations are correct.
        @test all(sum(sb.vectors, dims=1) .== N)

        # Each vector is unique.
        @test sb.vectors == unique(sb.vectors, dims=2)

        # Iteration works, serial numbers match up.
        d = 0
        for (i, v) in enumerate(sb)
            @test v in sb
            @test v == sb.vectors[:, i]
            @test v == sb[i]
            @test i == serial_num(sb, v)
            d += 1
        end
        @test d == D

        # Invalid vectors.
        v = zeros(Int, K)
        if N > 0
            @test !(v in sb)
        else
            @test v in sb
        end
        for n in [-1, N+1]
            v[1] = n
            @test !(v in sb)
        end
    end

    @test_throws DomainError Szbasis(0, 5)
    @test_throws DomainError Szbasis(5, -1)
end

@testset "RestrictedSzbasis" begin
    @testset for (K, N, M, D) in [(50, 0, 0, 1), (1, 1, 1, 1), (2, 3, 2, 2), (3, 2, 1, 3), (4, 4, 2, 19)]
        sb = RestrictedSzbasis(K, N, M)

        @test sb.K == K
        @test sb.N == N
        @test sb.M == M
        @test sb.D == D
        @test num_vectors(N, K, M) == D
        @test size(sb.vectors) == (sb.K, sb.D)
        @test length(sb) == D
        @test site_max(sb) == M

        # Occupations are correct.
        @test all(sum(sb.vectors, dims=1) .== N)

        # Each vector is unique.
        @test sb.vectors == unique(sb.vectors, dims=2)

        # Iteration works, serial numbers match up.
        d = 0
        for (i, v) in enumerate(sb)
            @test v in sb
            @test v == sb.vectors[:, i]
            @test v == sb[i]
            @test i == serial_num(sb, v)
            d += 1
        end
        @test d == D

        # Invalid vectors.
        v = zeros(Int, K)
        if N > 0
            @test !(v in sb)
        else
            @test v in sb
        end
        for n in [-1, N+1, M+1]
            v[1] = n
            @test !(v in sb)
        end
    end

    @test_throws DomainError RestrictedSzbasis(2, 5, 2)
end

@testset "sub_serial_num" begin
    let sb = Szbasis(3, 2)
        # Split each vector into 2 sites on the left and 1 on the right.
        left = []
        right = []

        for v in sb
            # Try range indexing.
            push!(left, sub_serial_num(sb, v[1:2]))
            # Try a view.
            push!(right, sub_serial_num(sb, view(v, 3:3)))
        end

        @test sort(left) == [1, 2, 3, 4, 5, 6]
        @test sort(right) == [1, 1, 1, 2, 2, 3]
    end

    let sb = RestrictedSzbasis(5, 2, 2)
        # Split each vector into 2 sites on the left and 3 on the right.
        left = []
        right = []

        for v in sb
            # Try range indexing.
            push!(left, sub_serial_num(sb, v[1:2]))
            # Try a view.
            push!(right, sub_serial_num(sb, view(v, 3:5)))
        end

        @test sort(left) == [1, 1, 1, 1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6]
        @test sort(right) == [1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 6, 7, 8, 9, 10]
    end
end

@testset "string" begin
    @testset for (s1, v, s2) in [(sz"", Int[], ""), (sz"0", [0], "0"), (sz"123", [123], "123"), (sz"0 0", [0, 0], "0 0"), (sz"1 23 456", [1, 23, 456], "1 23 456"), (sz"0 0 1 0", [0, 0, 1, 0], "0 0 1 0")]
        @test s1 == v
        @test s2 == to_str(v)
        @test s2 == to_str(v[1:end])
        @test s2 == to_str(@view v[1:end])
    end
end

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct AVRTypeI <: AVR
        Ka::Float64
        Ke::Float64
        Kf::Float64
        Ta::Float64
        Te::Float64
        Tf::Float64
        Tr::Float64
        Vr_max::Float64
        Vr_min::Float64
        Ae::Float64
        Be::Float64
        V_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int64
        internal::InfrastructureSystemsInternal
    end

Parameters of an Automatic Voltage Regulator Type I - Resembles IEEE Type DC1

# Arguments
- `Ka::Float64`: Amplifier Gain, validation range: (0, nothing)
- `Ke::Float64`: Field circuit integral deviation, validation range: (0, nothing)
- `Kf::Float64`: Stabilizer Gain in s * pu/pu, validation range: (0, nothing)
- `Ta::Float64`: Amplifier Time Constant in s, validation range: (0, nothing)
- `Te::Float64`: Field Circuit Time Constant in s, validation range: (0, nothing)
- `Tf::Float64`: Stabilizer Time Constant in s, validation range: (0, nothing)
- `Tr::Float64`: Voltage Measurement Time Constant in s, validation range: (0, nothing)
- `Vr_max::Float64`: Maximum regulator voltage in pu, validation range: (0, nothing)
- `Vr_min::Float64`: Minimum regulator voltage in pu, validation range: (0, nothing)
- `Ae::Float64`: 1st ceiling coefficient, validation range: (0, nothing)
- `Be::Float64`: 2nd ceiling coefficient, validation range: (0, nothing)
- `V_ref::Float64`: Reference Voltage Set-point, validation range: (0, nothing)
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int64`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct AVRTypeI <: AVR
    "Amplifier Gain"
    Ka::Float64
    "Field circuit integral deviation"
    Ke::Float64
    "Stabilizer Gain in s * pu/pu"
    Kf::Float64
    "Amplifier Time Constant in s"
    Ta::Float64
    "Field Circuit Time Constant in s"
    Te::Float64
    "Stabilizer Time Constant in s"
    Tf::Float64
    "Voltage Measurement Time Constant in s"
    Tr::Float64
    "Maximum regulator voltage in pu"
    Vr_max::Float64
    "Minimum regulator voltage in pu"
    Vr_min::Float64
    "1st ceiling coefficient"
    Ae::Float64
    "2nd ceiling coefficient"
    Be::Float64
    "Reference Voltage Set-point"
    V_ref::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    n_states::Int64
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function AVRTypeI(Ka, Ke, Kf, Ta, Te, Tf, Tr, Vr_max, Vr_min, Ae, Be, V_ref=1.0, ext=Dict{String, Any}(), )
    AVRTypeI(Ka, Ke, Kf, Ta, Te, Tf, Tr, Vr_max, Vr_min, Ae, Be, V_ref, ext, [:Vf, :Vr1, :Vr2, :Vm], 4, InfrastructureSystemsInternal(), )
end

function AVRTypeI(; Ka, Ke, Kf, Ta, Te, Tf, Tr, Vr_max, Vr_min, Ae, Be, V_ref=1.0, ext=Dict{String, Any}(), )
    AVRTypeI(Ka, Ke, Kf, Ta, Te, Tf, Tr, Vr_max, Vr_min, Ae, Be, V_ref, ext, )
end

# Constructor for demo purposes; non-functional.
function AVRTypeI(::Nothing)
    AVRTypeI(;
        Ka=0,
        Ke=0,
        Kf=0,
        Ta=0,
        Te=0,
        Tf=0,
        Tr=0,
        Vr_max=0,
        Vr_min=0,
        Ae=0,
        Be=0,
        V_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get AVRTypeI Ka."""
get_Ka(value::AVRTypeI) = value.Ka
"""Get AVRTypeI Ke."""
get_Ke(value::AVRTypeI) = value.Ke
"""Get AVRTypeI Kf."""
get_Kf(value::AVRTypeI) = value.Kf
"""Get AVRTypeI Ta."""
get_Ta(value::AVRTypeI) = value.Ta
"""Get AVRTypeI Te."""
get_Te(value::AVRTypeI) = value.Te
"""Get AVRTypeI Tf."""
get_Tf(value::AVRTypeI) = value.Tf
"""Get AVRTypeI Tr."""
get_Tr(value::AVRTypeI) = value.Tr
"""Get AVRTypeI Vr_max."""
get_Vr_max(value::AVRTypeI) = value.Vr_max
"""Get AVRTypeI Vr_min."""
get_Vr_min(value::AVRTypeI) = value.Vr_min
"""Get AVRTypeI Ae."""
get_Ae(value::AVRTypeI) = value.Ae
"""Get AVRTypeI Be."""
get_Be(value::AVRTypeI) = value.Be
"""Get AVRTypeI V_ref."""
get_V_ref(value::AVRTypeI) = value.V_ref
"""Get AVRTypeI ext."""
get_ext(value::AVRTypeI) = value.ext
"""Get AVRTypeI states."""
get_states(value::AVRTypeI) = value.states
"""Get AVRTypeI n_states."""
get_n_states(value::AVRTypeI) = value.n_states
"""Get AVRTypeI internal."""
get_internal(value::AVRTypeI) = value.internal

"""Set AVRTypeI Ka."""
set_Ka!(value::AVRTypeI, val::Float64) = value.Ka = val
"""Set AVRTypeI Ke."""
set_Ke!(value::AVRTypeI, val::Float64) = value.Ke = val
"""Set AVRTypeI Kf."""
set_Kf!(value::AVRTypeI, val::Float64) = value.Kf = val
"""Set AVRTypeI Ta."""
set_Ta!(value::AVRTypeI, val::Float64) = value.Ta = val
"""Set AVRTypeI Te."""
set_Te!(value::AVRTypeI, val::Float64) = value.Te = val
"""Set AVRTypeI Tf."""
set_Tf!(value::AVRTypeI, val::Float64) = value.Tf = val
"""Set AVRTypeI Tr."""
set_Tr!(value::AVRTypeI, val::Float64) = value.Tr = val
"""Set AVRTypeI Vr_max."""
set_Vr_max!(value::AVRTypeI, val::Float64) = value.Vr_max = val
"""Set AVRTypeI Vr_min."""
set_Vr_min!(value::AVRTypeI, val::Float64) = value.Vr_min = val
"""Set AVRTypeI Ae."""
set_Ae!(value::AVRTypeI, val::Float64) = value.Ae = val
"""Set AVRTypeI Be."""
set_Be!(value::AVRTypeI, val::Float64) = value.Be = val
"""Set AVRTypeI V_ref."""
set_V_ref!(value::AVRTypeI, val::Float64) = value.V_ref = val
"""Set AVRTypeI ext."""
set_ext!(value::AVRTypeI, val::Dict{String, Any}) = value.ext = val
"""Set AVRTypeI states."""
set_states!(value::AVRTypeI, val::Vector{Symbol}) = value.states = val
"""Set AVRTypeI n_states."""
set_n_states!(value::AVRTypeI, val::Int64) = value.n_states = val
"""Set AVRTypeI internal."""
set_internal!(value::AVRTypeI, val::InfrastructureSystemsInternal) = value.internal = val


text:
module VisualizeFeatures

using Statistics
using StatsBase
using DataFrames
using Mongoc
using LightGraphs
using Plots
using Interact

include("../utils/misc.jl")
include("../utils/mongo.jl")
include("df.jl")
include("plots.jl")

end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
# module SocketIOs
# export SocketIO, run_taps

import JSON
import Base.read, Base.write, Base.sleep

using Sockets

include("../utils/antenna.jl")
Taps = include("../taps.jl")

using .Taps

mutable struct SocketIO
            host; port; clienthost; clientport; tcp; tcp_server; keep_open; is_tcp_opened; MPI; comm; root; nprc; rank; mFloat64; mUInt8; cache::Dict;
end

function Base.setproperty!(value::SocketIO, name::Symbol, x)
    if name == :host
        x = Sockets.IPv4(host)
    end
    setfield!(value, name, x)
end

SocketIO(;host="127.0.0.1", port=6543, clienthost="127.0.0.1", clientport=6544,
          keep_open=true, is_tcp_opened=false, MPI=nothing, comm=nothing,
          root=0, nprc=nothing, rank=nothing,
          mFloat64=nothing, mUInt8=nothing, cache=Dict(), kwargs...) = begin
    host = Sockets.IPv4(host)
    clienthost = Sockets.IPv4(clienthost)
    tcp_server = rank != root ? nothing : is_tcp_opened ? nothing : listen(host, port)
    is_tcp_opened = true
    # It should be nothing at the beginning
    tcp = nothing

    SocketIO(host, port, clienthost, clientport, tcp, tcp_server, keep_open, is_tcp_opened, MPI, comm, root, nprc, rank, mFloat64, mUInt8, cache)
end

function run_taps(io::SocketIO)
    # Pre define static variables, remove `io.`
    host = io.host; port = io.port; tcp_server = io.tcp_server;
    MPI = io.MPI; comm = io.comm; root = io.root; nprc = io.nprc; rank = io.rank;
    mFloat64 = io.mFloat64; mUInt8 = io.mUInt8
    # MPI set menu
    mpi_args = [MPI, comm, root, nprc, rank]
    mpi_type_args = [mFloat64, mUInt8]
    mpi_full_args = [mpi_args..., mpi_type_args...]
    # socket set menu
    socket_args = [host, port, tcp_server]

    # Start
    while io.keep_open
        io.tcp = rank == root ? accept(tcp_server) : nothing
        len = rank == root ? reinterpret(Int64, read(io.tcp, 8))[1] : nothing
        instruction = rank == root ? read(io.tcp, len) : nothing
        rank == root ? poke(mpi_args...) : sleep(mpi_args...)
        rank == root ? echo(io, Array{UInt8, 1}(reinterpret(UInt8, [len]))) : nothing
        operate(io, instruction)
    end
end

function operate(io::SocketIO, instruction::T) where {T<:Union{Array{UInt8, 1}, Nothing}}
    instruction = instruction == nothing ? b"1" : instruction
    instructiontype = [io.MPI.bcast(instruction[1], io.root, io.comm)]
    if instructiontype == b"*" # WildCard Instruction
        instruction = io.MPI.bcast(instruction, io.root, io.comm)
        eval(Meta.parse(String(instruction[2:end])))
    elseif instructiontype in [b"0", b"1", b"2"] # Semi wild card
        instruction = io.MPI.bcast(instruction, io.root, io.comm)
        howlongisurname = reinterpret(Int64, instruction[2:9])[1]
        method = String(instruction[10:9+howlongisurname])
        argsbytes = instruction[10+howlongisurname:end]
        args, kwargs = unpacking(Array{UInt8, 1}(argsbytes[9:end]))
        if instructiontype == b"0"
            eval(Meta.parse(method))(args...; kwargs...)
        elseif instructiontype == b"1"
            eval(Meta.parse(method))(io, args...; kwargs...)
        elseif instructiontype == b"2"
            eval(Meta.parse(method))(io.cache["paths"], args...; kwargs...)
        end
    end
end

function poke(MPI, comm, root, nprc, rank)
    send_temp = Array{Int64}(undef, 1)
    fill!(send_temp, Int64(0))
    for dst=1:nprc-1 MPI.Isend(send_temp, dst, root+32, comm) end
end

function sleep(MPI, comm, root, nprc, rank)
    while !(MPI.Iprobe(root, root+32, comm)[1]) sleep(0.001) end
    recv_temp = Array{Int64}(undef, 1)
    req = MPI.Irecv!(recv_temp, root, root+32, comm)
end

echo(io::SocketIO, instruction::Array{UInt8, 1}) = write(io.tcp, instruction)
ping(io::SocketIO) = io.rank == io.root ? write(io.tcp, reinterpret(UInt8, [time()])) : nothing
shutdown(io::SocketIO) = (io.keep_open = false)

function set_client(io::SocketIO; clienthost="127.0.0.1", clientport=6544)
    io.clienthost = Sockets.IPv4(clienthost)
    io.clientport = clientport
end

"""
Dict(
"coords"=> class_of_coords
"coords_kwargs" => Dict("coords"=>Array, ...)
"model"=> class_of_model
"model_kwargs"=> Dict("model")

"""
function read(io::SocketIO, args...; kwargs...)
    io.cache["paths"] = Paths(args...; kwargs...)
end

function read_parallel(io::SocketIO, args...; kwargs...)
    host, port = io.clienthost, io.clientport
    args = [args...]
    kwargs = Dict{Symbol, Any}([kwargs...])
    ptcp = connect(host, port)
    write(ptcp, reinterpret(UInt8, [io.rank]))
    len = reinterpret(Int64, read(ptcp, 8))[1]
    bytesarr = read(ptcp, len)
    pargs, pkwargs = unpacking(bytesarr)
    # Merge two lists that element is nothing.
    for i=1:length(args)
        args[i] = args[i] == nothing ? pargs[i] : args[i]
    end
    merge!(kwargs, pkwargs)
    read(io, args...; kwargs...)
end

function update(io::SocketIO, args...; kwargs...)
    for (key, value) in pairs(kwargs)
        setproperty!(io.cache["paths"], key, value)
    end
end

function update_parallel(io::SocketIO, args...; kwargs...)
    host, port = io.clienthost, io.clientport
    args = [args...]
    kwargs = Dict{Symbol, Any}([kwargs...])

    ptcp = connect(host, port)
    write(ptcp, reinterpret(UInt8, [io.rank]))
    len = reinterpret(Int64, read(ptcp, 8))[1]
    bytesarr = read(ptcp, len)
    pargs, pkwargs = unpacking(bytesarr)
    # Merge two lists that element is nothing.
    for i=1:length(args)
        args[i] = args[i] == nothing ? pargs[i] : args[i]
    end
    merge!(kwargs, pkwargs)
    update(io, args...; kwargs...)
end

function write(io::SocketIO, args...; sender=0, kwargs...)
    host, port = io.clienthost, io.clientport
    if io.rank == sender
        paths = io.cache["paths"]
        pargs = Array{Any, 1}([])
        pkwargs = Dict()
        for key in args
            push!(pargs, getproperty(paths, key))
        end
        for (key, value) in pairs(kwargs)
            pkwargs[key] = getproperty(paths, key, value)
        end
        argsbytes = packing(pargs...; pkwargs...)

        tcp = connect(host, port)
        write(tcp, argsbytes)
    end
end

function write_parallel(io::SocketIO, args...; senders="all", kwargs...)
    host, port = io.clienthost, io.clientport
    args = [args...]
    kwargs = Dict{Symbol, Any}([kwargs...])

#    sender = senders == "all" ? io.rank : io.rank in senders ? io.rank : nothing

#    if sender == io.rank
        ptcp = connect(host, port)
        write(ptcp, reinterpret(UInt8, [io.rank]))
        len = reinterpret(Int64, read(ptcp, 8))[1]
        bytesarr = read(ptcp, len)
        pargs, pkwargs = unpacking(bytesarr)

        args = cat(args, pargs, dims=1)
        merge!(kwargs, pkwargs)

        paths = io.cache["paths"]
        pathargs = Array{Any, 1}([])
        pathkwargs = Dict()
        for key in args
            push!(pathargs, getproperty(paths, key))
        end
        for (key, value) in pairs(kwargs)
            pathkwargs[key] = getproperty(paths, key, value)
        end
        argsbytes = packing(pathargs...; pathkwargs...)
        write(ptcp, argsbytes)
end
#
# function get_properties(io::SocketIO, properties, args...; kwargs...)
#     paths = io.cache["paths"]
#     get_properties(paths, properties, args...; kwargs...)
# end

# end

max_stars_repo_name: schinavro/taps

max_issues_repo_name: schinavro/taps

max_forks_repo_name: schinavro/taps


text:
using Documenter, DiffRules

makedocs(modules=[DiffRules],
         doctest = false,
         sitename = "DiffRules",
         pages = ["Documentation" => "index.md"],
         format = Documenter.HTML(
                  prettyurls = get(ENV, "CI", nothing) == "true"
         ),
)

deploydocs(repo = "github.com/JuliaDiff/DiffRules.jl")