text:
using SQLite
using NeXLCore
using IntervalSets
"""
A `DBCampaign` represents a collection of measurements of the same material.  It includes
both the unknown spectra and the reference spectra necessary to fit it as well as some
contextual data.
"""
struct DBCampaign
    database::SQLite.DB
    pkey::Int
    project::DBProject
    detector::DBDetector
    elements::Vector{Element}
    material::Union{Missing, Material}
    disposition::String
    fitspectrum::Vector{DBFitSpectrum}
    refspectrum::Vector{DBReference}
end

function NeXLUncertainties.asa(::Type{DataFrame}, campaign::DBCampaign)
    prop, values = String[], String[]
    push!(prop, "Primary Key"), push!(values, "$(campaign.pkey)")
    push!(prop, "Project"), push!(values, repr(campaign.project))
    push!(prop, "Analyst"), push!(values, campaign.project.createdBy.name)
    push!(prop, "Elements"), push!(values, "$(_elmstostr(campaign.elements))")
    push!(prop, "Material"), push!(values, repr(campaign.material))
    push!(prop, "Disposition"), push!(values, campaign.disposition)
    res = DataFrame(Property=prop, Value=values)
    return vcat(res, asa(DataFrame, campaign.detector))
end

function Base.show(io::IO, campaign::DBCampaign)
    println(io, "      Index = $(campaign.pkey)")
    println(io, "    Project = $(repr(campaign.project))")
    println(io, "   Elements = $(_elmstostr(campaign.elements))")
    println(io, "   Detector = $(repr(campaign.detector))")
    if !ismissing(campaign.material)
        println(io, "   Material = $(repr(campaign.material))")
    end
    println(io, "Disposition = $(campaign.disposition)")
    print(io, "==== Unknowns ====")
    for dbf in campaign.fitspectrum
        print(io, "\n\t$dbf")
    end
    print(io, "\n==== References ====")
    for elm in campaign.elements
        print(io, "\n\t$(symbol(elm)) = $(join(dbreferences(campaign, elm),','))")
    end
end

function NeXLUncertainties.asa(::Type{DataFrame}, camps::AbstractVector{DBCampaign})
    nm(m::Missing) = "Unknown"
    nm(m::Material) = name(m)
    DataFrame(
        Key = [ camp.pkey for camp in camps ],
        Material = [ camp.material for camp in camps ],
        E0 = [ camp.fitspectrum[1].spectrum.beamenergy for camp in camps],
        Unknowns = [ length(camp.fitspectrum) for camp in camps ],
        References = [ join([ nm(ref.spectrum.composition) for ref in camp.refspectrum],", ") for camp in camps ],
        Project = [ repr(camp.project) for camp in camps ],
        Analyst = [ camp.project.createdBy.name for camp in camps ],
        Disposition = [camp.disposition for camp in camps ],
    )
end

"""
    measured(campaign::DBCampaign)::Vector{Spectrum}
Return a vector containing the spectra to be fit.
"""
measured(campaign::DBCampaign)::Vector{Spectrum} = map(fs -> asa(Spectrum, fs), campaign.fitspectrum)

"""
    references(campaign::DBFitSpectra, elm::Element)::Vector{Spectrum})
Return a vector of the spectra which could be used as references for the Element.
"""
NeXLSpectrum.references(campaign::DBCampaign, elm::Element)::Vector{Spectrum} = map(ref->asa(Spectrum,ref), dbreferences(campaign, elm))

"""
    dbreferences(campaign::DBCampaign, elm::Element)::Vector{DBReference}

Return a vector of the spectra which could be used as references for the Element as DBReference objects.
"""
dbreferences(campaign::DBCampaign, elm::Element)::Vector{DBReference} =
    collect(filter(rs -> elm in rs.elements, campaign.refspectrum))

NeXLCore.elms(campaign::DBCampaign) = campaign.elements

_elmstostr(elms::Vector{Element}) = join(symbol.(elms), ',')
_strtoelms(str::String) = parse.(Element, strip.(split(str, ',')))

function Base.write(db::SQLite.DB, ::Type{DBCampaign}, projKey::Int, detKey::Int, elms::Vector{Element}, matkey::Integer, disposition::String="Pending")::Int
    stmt1 = SQLite.Stmt(db, "INSERT INTO CAMPAIGN(PROJECT, DETECTOR, DISPOSITION, ELEMENTS, MATKEY ) VALUES( ?, ?, ?, ?, ?);")
    q = DBInterface.execute(stmt1, (projKey, detKey, disposition, _elmstostr(elms), matkey))
    return DBInterface.lastrowid(q)
end

function Base.read(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int)::DBCampaign
    stmt1 = SQLite.Stmt(db, "SELECT * FROM CAMPAIGN WHERE PKEY=?;")
    q1 = DBInterface.execute(stmt1, (campaign,))
    if SQLite.done(q1)
        error("No fit spectra record with campaign = $campaign.")
    end
    r1 = SQLite.Row(q1)
    @assert r1[:PKEY] == campaign "Mismatching pkey in DBFitSpectrum"
    project = read(db, DBProject, r1[:PROJECT])
    detector = read(db, DBDetector, r1[:DETECTOR])
    elms = _strtoelms(r1[:ELEMENTS])
    material = r1[:MATKEY] != -1 ? read(db, Material, r1[:MATKEY]) : missing
    disposition = r1[:DISPOSITION]
    stmt2 = SQLite.Stmt(db, "SELECT * FROM FITSPECTRUM WHERE CAMPAIGN=?;")
    q2 = DBInterface.execute(stmt2, (campaign,))
    tobefit = DBFitSpectrum[]
    for r2 in q2
        @assert r2[:CAMPAIGN] == campaign
        spec = read(db, DBSpectrum, r2[:SPECTRUM])
        push!(tobefit, DBFitSpectrum(db, campaign, spec))
    end
    stmt3 = SQLite.Stmt(db, "SELECT * FROM REFERENCESPECTRUM WHERE CAMPAIGN=?;")
    q3 = DBInterface.execute(stmt3, (campaign,))
    refs = DBReference[]
    for r3 in q3
        @assert r3[:CAMPAIGN] == campaign
        spec = read(db, DBSpectrum, r3[:SPECTRUM])
        push!(refs, DBReference(db, r3[:PKEY], r3[:CAMPAIGN], spec, _strtoelms(r3[:ELEMENTS])))
    end
    return DBCampaign(db, campaign, project, detector, elms, material, disposition, tobefit, refs)
end

function kratios(campaign::DBCampaign, elm::Union{Element,Nothing}=nothing, mink::Float64=0.0)::Vector{DBKRatio}
    return findall(campaign.database, DBKRatio, campaign=campaign.pkey, elm=elm, mink=mink)
end

"""
    disposition(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int)

Reports the disposition of this set of fitted spectra.

  * `"Pending"` - Entered in the database but not reviewed
  * `"Rejected: "*msg` - Rejected for the reason specified in `msg`
  * `"Accepted: "*msg` - Reviewed and accepted with optional `msg`
"""
function disposition(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int)
    stmt1 = SQLite.Stmt(db, "SELECT DISPOSITION FROM CAMPAIGN WHERE PKEY=?;")
    q1 = DBInterface.execute(stmt1, (campaign,))
    if SQLite.done(q1)
        error("No fit spectra record with campaign = $campaign.")
    end
    return SQLite.Row(q1)[:DISPOSITION]
end
disposition(camp::DBCampaign) = camp.disposition

"""
    disposition!(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int, value::String)

Sets the disposition of this set of fitted spectra.

  * `"Pending"` - Entered in the database but not reviewed
  * `"Rejected: "*msg` - Rejected for the reason specified in `msg`
  * `"Accepted: "*msg` - Reviewed and accepted with optional `msg`
"""
function disposition!(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int, value::AbstractString)
    @assert isequal("Pending", value) || startswith(value, "Rejected") || startswith("Accepted")
        "The disposition must be \"Pending\" or start with \"Rejected\" or \"Accepted\""
    stmt1 = SQLite.Stmt(db, "UPDATE CAMPAIGN SET DISPOSITION = ? WHERE PKEY=?;")
    DBInterface.execute(stmt1, (campaign, value))
    return value
end
function disposition!(campaign::DBCampaign, value::AbstractString)
    return disposition!(campaign.database, DBCampaign, campaign.pkey, value)
end

function Base.delete!(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int)
    stmt1 = SQLite.Stmt(db, "DELETE FROM CAMPAIGN WHERE PKEY=?;")
    stmt2 = SQLite.Stmt(db, "DELETE FROM REFERENCESPECTRUM WHERE CAMPAIGN=?;")
    stmt3 = SQLite.Stmt(db, "DELETE FROM FITSPECTRUM WHERE CAMPAIGN=?;")
    stmt4 = SQLite.Stmt(db, "DELETE FROM KRATIO WHERE FITSPEC=?;")
    DBInterface.execute(stmt4, (campaign,))
    DBInterface.execute(stmt3, (campaign,))
    DBInterface.execute(stmt2, (campaign,))
    DBInterface.execute(stmt1, (campaign,))
end

function Base.delete!(campaign::DBCampaign)
    delete!(campaign.database, DBCampaign, campaign.pkey)
end

function NeXLUncertainties.asa(::Type{DataFrame}, db::SQLite.DB, ::Type{DBCampaign})
    stmt1 = SQLite.Stmt(db, "SELECT PKEY FROM CAMPAIGN")
    q1 = DBInterface.execute(stmt1)
    return SQLite.done(q1) ? DataFrame() : asa(DataFrame, [ read(db, DBCampaign, row[:PKEY]) for row in res ])
end

function Base.findall(db::SQLite.DB, ::Type{DBCampaign}; material::Union{String,Missing}=missing, project::Union{DBProject,Missing} = missing, det::Union{DBDetector,Missing}=missing)::Vector{DBCampaign}
    keys, args = String[], Any[]
    if !ismissing(material)
        matkey = find(db, Material, material)
        @assert matkey!=-1 "Unable to find the material named `$material`."
        push!(keys, "MATKEY=?")
        push!(args, matkey)
    end
    if !ismissing(project)
        push!(keys,"PROJECT=?")
        push!(args, project.pkey)
    end
    if !ismissing(det)
        push!(keys, "DETECTOR=?")
        push!(args, det.pkey)
    end
    if length(keys)==0
        @error "Please specify either a project or a detector."
    end
    stmt=SQLite.Stmt(db, "SELECT PKEY FROM CAMPAIGN WHERE $(join(keys," AND "));")
    q1 = DBInterface.execute(stmt, args)
    return [ read(db, DBCampaign, r[:PKEY] ) for r in q1 ]
end

"""
    NeXLMatrixCorrection.quantify(#
        db::SQLite.DB, #
        ::Type{DBCampaign}, #
        campaign::Int; #
        strip::AbstractVector{Element} = Element[],
        mc::Type{<:MatrixCorrection} = XPP,
        fl::Type{<:FluorescenceCorrection} = ReedFluorescence,
        cc::Type{<:CoatingCorrection} = Coating,
        kro::KRatioOptimizer = SimpleKRatioOptimizer(1.5),
    )::Vector{IterationResult}

Quantify the k-ratios in the database associated with te specified DBCampaign.
"""
function NeXLMatrixCorrection.quantify(#
    db::SQLite.DB, #
    ::Type{DBCampaign}, #
    campaign::Int; #
    args...
)::Vector{IterationResult}
    quantify(findall(db, DBKRatio, campaign=campaign, mink=0.0), args...)
end

function Base.write(
    db::SQLite.DB, #
    ::Type{DBCampaign}, #
    project::DBProject, # The project
    sample::DBSample,   # The sample containing the unknown
    unkComp::Union{Material,Missing},  # The unknown material or missing
    detector::DBDetector, # The detector on which the data was collected
    analyst::DBPerson, # The person who collected the data
    e0::Float64, # The beam energy for the unknown in eV
    measSpectra::Vector{String}, # The files containing the measured spectra
    refSpectra::Vector{Tuple{DBSample, Material{T,U}, String, Float64, Vector{Element}}}, # The reference spectra (sample, material, filename, e0, extra elements)
    extraElms::AbstractVector{Element} = Element[], #
)::Int where { T<:AbstractFloat, U<:AbstractFloat}
    SQLite.transaction(db) do # All or nothing...
        elms, matkey = if ismissing(unkComp)
            extraElms, -1
        else
            unique(append!(collect(keys(unkComp)), extraElms)), write(db, unkComp)
        end
        campaign = write(db, DBCampaign, project.pkey, detector.pkey, elms, matkey)
        dt = now()
        for fn in measSpectra
            format, name = sniffFormat(fn), splitext(splitdir(fn)[2])[1]
            spec = write(db, DBSpectrum, detector, e0, unkComp, analyst, sample, dt, name, fn, format)
            write(db, NeXLDatabase.DBFitSpectrum, campaign, spec)
        end
        for (samp, std, fn, e0s, elms) in refSpectra
            format, name = sniffFormat(fn), splitext(splitdir(fn)[2])[1]
            spec = write(db, DBSpectrum, detector, e0s, std, analyst, sample, dt, name, fn, format)
            refelms = ismissing(std) ? elms : append!(collect(keys(std)), elms)
            ref = write(db, DBReference, campaign, spec, refelms)
        end
        return campaign
    end
end

"""
    NeXLSpectrum.fit_spectrum(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int, unkcomp::Union{Material, Missing}=missing, update=true)

Fit a collection of spectra from the database against the reference spectra associated in the database with it.  If
`unkcomp` is provided then previous k-ratio results are deleted and the new k-ratio results are written to the database.
"""
function NeXLSpectrum.fit_spectrum(db::SQLite.DB, ::Type{DBCampaign}, campaign::Int, unkcomp::Union{Material, Missing}=missing; update::Bool=true)::Vector{FilterFitResult}
    fs = read(db, DBCampaign, campaign)
    unks = measured(fs)
    det = convert(BasicEDS, fs.detector)
    refs =  NeXLSpectrum.ReferencePacket[]
    specs = [ asa(Spectrum, ref.spectrum) for ref in fs.refspectrum ]
    for elm in fs.elements
        # Make sure that the candidate reference material has some of this element
        for spec in filter(sp->haskey(sp,:Composition) && (value(sp[:Composition][elm]) > 0.01) , specs)
			push!(refs, reference(elm, spec))
		end
	end
	ffp = NeXLSpectrum.references(refs, det)
    ress = map(sp->fit_spectrum(sp, ffp), unks)
    update && SQLite.transaction(db) do
        # Remove previous k-ratios for this DBCampaign
        DBInterface.execute(SQLite.Stmt(db, "DELETE FROM KRATIO WHERE CAMPAIGN=?;"), (campaign, ))
        for (unk, res) in zip(unks, ress)
            # @show campaign, res
            write(db, DBKRatio, fs, unk, unkcomp, res)
        end
        true
    end
    return ress
end


text:
# code for generating uniform random unitary matrices
# TODO for now we'll just do SU(2)
# TODO generators should probably be of the Julia Hermitian type

# obviously this does not in itself guarantee unitarity
const UnitaryMatrix = Matrix{Complex{Float64}}

# will eventually need to replace with boilerplate code for random SU(N) generators
const σ₁ = Complex{Float64}[0 1; 1 0]
const σ₂ = Complex{Float64}[0 -im; im 0]
const σ₃ = Complex{Float64}[1 0; 0 -1]

const σ⃗ = [σ₁, σ₂, σ₃]

const one2by2 = eye(Complex{Float64}, 2)

"""
    randunitvec(N::Integer)

Creates are random unit vector in Rᴺ.
"""
function randunitvec(::Type{T}, N::Integer) where T
    x = randn(T, N)
    r = vecnorm(x)
    x/r
end
randunitvec(N::Integer) = randunitvec(Float64, N)

"""
    randazimuth()

Creates a random number on [0, 2π).
"""
randazimuth(::Type{T}) where T = 2π*rand(T)
randazimuth() = randazimuth(Float64)

"""
    randSU2_ineff()

An elegant way of generating uniformly distributed random SU(2)
matrices.  This version appears inefficient, as it probably does
lots of dense matrix multiplication.
"""
randSU2_ineff(ϕ::AbstractFloat, n::AbstractVector) = expm(im*ϕ*n'*σ⃗)
randSU2_ineff() = randSU2_ineff(randazimuth(), randunitvec(3))
export randSU2_ineff

"""
    SU2(ϕ::AbstractFloat, n::AbstractVector)

Creat a vector from SU(2) with phase `ϕ*n`.
"""
function SU2(ϕ::AbstractFloat, n::AbstractVector)
    one2by2*cos(ϕ) + im*(n'*σ⃗)*sin(ϕ)
end
export SU2

"""
    randSU2(ϕmax::AbstractFloat)
    randSU2()

Generate a random matrix uniformly distributed on SU(2).
Only involves scalar multiplication and dense matrix addition.

If Φmax is provided, the magnitude of the phase will lie in the
interval [0, ϕmax) (the randomly generated unit vector can be negative).
"""
randSU2(ϕmax::AbstractFloat) = SU2(ϕmax*rand(), randunitvec(3))
randSU2() = randSU2(randazimuth())
export randSU2

"""
    tr(A)

Returns the trace of a matrix. This is simply a convenient alias
for the Julia function `trace`.
"""
tr(A) = trace(A)
export tr

ext: jl


text:
import UMAP

function plot_noise_estimation_diagnostics(edge_lengths::Vector{Float64}, confidences::Vector{Float64}, d1::T, d2::T; title::String="Noise estimation",
        confidence_nn_id::Union{Int, String}="k", linewidth::Float64=4.0, bins::Int=50) where T <: Distributions.UnivariateDistribution
        x_max = quantile(edge_lengths, 0.99);
        x_vals = range(0, x_max, length=1000)
        n1 = sum(confidences);
        n2 = length(confidences) - n1;

        p_df = estimate_hist(edge_lengths[edge_lengths .< x_max], normalize=true, nbins=bins)
        p_df[!, :intra] = n1 / (n1 + n2) .* pdf.(d1, p_df.s)
        p_df[!, :bg] = n2 / (n1 + n2) .* pdf.(d2, p_df.s)

        return p_df |>
        @vlplot(x={:s, title="Distance to $(confidence_nn_id)'th nearest neighbor"}, title=title, width=400, height=300) +
        @vlplot(:bar, x2=:e, y={:h, title="Density"}, color={datum="Observed", scale={scheme="category10"}, legend={title="Distribution"}}) +
        @vlplot({:line, size=linewidth}, y=:bg, color={datum="Background"}) +
        @vlplot({:line, size=linewidth}, y=:intra, color={datum="Intracellular"})
end

function plot_num_transcript_overview(genes::Vector{Int}, confidences::Vector{Float64}, gene_names::Vector; alpha::Float64=0.3)
    order = sortperm(gene_names)
    return plot_expression_vectors(
        count_array(genes[confidences .>= 0.5], max_value=length(gene_names))[order],
        count_array(genes[confidences .< 0.5], max_value=length(gene_names))[order],
        gene_names=gene_names[order]; labels=["Real", "Noise"], ylabel="Num. molecules",
        min_expr_frac=0.01, alpha=alpha
    )
end

function plot_gene_structure(df_spatial::DataFrame, gene_names::Vector, confidence::Vector{Float64}=df_spatial.confidence)
    adjacent_points, adjacent_weights = build_molecule_graph(df_spatial, filter=false)[1:2];
    cor_mat = pairwise_gene_spatial_cor(df_spatial.gene, df_spatial.confidence, adjacent_points, adjacent_weights);
    cor_vals = vec(cor_mat)
    min_cor = quantile(cor_vals[cor_vals .> 0], 0.001) ./ 5;
    max_cor = quantile(cor_vals, 0.99);
    p_dists = 1 .- max.(min.(cor_mat, max_cor), min_cor) ./ max_cor;
    p_dists[diagind(p_dists)] .= 0.0;

    embedding = UMAP.umap(p_dists, 2; metric=:precomputed, spread=1.0, min_dist=0.1, n_epochs=5000, n_neighbors=max(min(15, length(gene_names) ÷ 2), 2));
    marker_sizes = log.(count_array(df_spatial.gene));

    p_df = DataFrame(Dict(:x => embedding[1,:], :y => embedding[2,:], :gene => Symbol.(gene_names), :size => marker_sizes));

    return p_df |>
    @vlplot(
        x={:x, scale={domain=val_range(p_df.x)}, title="UMAP-1"},
        y={:y, scale={domain=val_range(p_df.y)}, title="UMAP-2"},
        size={:size, scale={range=[5, 10]}, legend=false},
        tooltip={:gene, type="nominal"},
        width=600, height=600, title="Gene local structure",
        config={axis={grid=false, ticks=false, ticklabels=false, labels=false}}
    ) +
    @vlplot(:text, text={:gene, type="nominal"}, selection={view={type=:interval, bind=:scales}}) +
    @vlplot({:point, filled=true})
end

function estimate_panel_plot_size(df_spatial::DataFrame, min_molecules_per_cell::Int, min_pixels_per_cell::Int=7)
    n_cells_per_side = sqrt(size(df_spatial, 1) / min_molecules_per_cell)
    plot_size = min_pixels_per_cell * n_cells_per_side
    x_rng = val_range(df_spatial.x)
    y_rng = val_range(df_spatial.y)
    y_ratio = (y_rng[2] - y_rng[1]) / (x_rng[2] - x_rng[1])

    y_ratio = y_ratio^0.5;
    x_ratio = 1 / y_ratio^0.5;
    plot_size = (x_ratio * plot_size, y_ratio * plot_size)
    return plot_size, n_cells_per_side
end

plot_dataset_colors(df_spatial::DataFrame, colors::Symbol; kwargs...) =
    plot_dataset_colors(df_spatial, df_spatial[!, colors]; kwargs...)

function plot_dataset_colors(df_spatial::DataFrame, color::Union{Vector, Symbol, String}; min_molecules_per_cell::Int, min_pixels_per_cell::Int=7, markersize::Float64=-1., alpha::Union{Float64, Vector{Float64}}=0.5,
        prior_polygons::Array{Matrix{Float64}, 1}=Matrix{Float64}[], polygons::Array{Matrix{Float64}, 1}=Matrix{Float64}[], ticks::Bool=true, axis_kwargs::KWArgT=nothing, kwargs...)

    axis_kwargs = update_args((xticklabelsize=12, yticklabelsize=12), axis_kwargs)
    plot_size = estimate_panel_plot_size(df_spatial, min_molecules_per_cell, min_pixels_per_cell)[1]
    if markersize < 0
        markersize = min(max(50.0 / min_molecules_per_cell, 0.25), 5)
    end

    fig = MK.Figure(resolution=plot_size)
    fig[1, 1] = MK.Axis(fig; xticksvisible=ticks, yticksvisible=ticks, axis_kwargs...);

    if length(prior_polygons) > 0
        MK.poly!([MK.Point2.(eachrow(p)) for p in prior_polygons]; strokecolor="darkred", color=Colors.RGBA(1, 0.65, 0, 0.25), strokewidth=0.5)
    end

    return plot_molecules!(df_spatial, polygons; color=color, markersize=markersize, alpha=alpha, kwargs...)
end

function plot_confidence_distribution(confidence::Vector{Float64}, assignment::Vector{<:Integer}; bins::AbstractVector{Float64}=0.0:0.025:1.025, size=(500, 250))
    v1 = confidence[assignment .!= 0]
    v2 = confidence[assignment .== 0]
    p_df = estimate_hist(v1, bins=bins)
    p_df[!, :h2] = estimate_hist(v2, bins=bins).h;

    return p_df |>
    @vlplot(x={:s, title="Confidence"}, x2=:e, width=size[1], height=size[2], title="Confidence per molecule") +
    @vlplot(:bar, y={:h, title="Num. molecules"}, color={datum="Assigned molecules"}) +
    @vlplot({:bar, opacity=0.5}, y=:h2, color={datum="Noise molecules"})
end

max_stars_repo_name: repushko/Baysor

max_issues_repo_name: repushko/Baysor

max_forks_repo_name: repushko/Baysor


text:
push!(LOAD_PATH,"../src/")

using Test, Documenter, LatticeUtilities

doctest(LatticeUtilities)


text:
# ------------------------------------------------------------------------------------------
# # Callbacks and Events
#
# In working with a differential equation, our system will evolve through many states.
# Particular states of the system may be of interest to us, and we say that an ***"event"***
# is triggered when our system reaches these states. For example, events may include the
# moment when our system reaches a particular temperature or velocity. We ***handle*** these
# events with ***callbacks***, which tell us what to do once an event has been triggered.
#
# These callbacks allow for a lot more than event handling, however. For example, we can use
# callbacks to achieve high-level behavior like exactly preserve conservation laws and save
# the trace of a matrix at pre-defined time points. This extra functionality allows us to
# use the callback system as a modding system for the DiffEq ecosystem's solvers.
#
# This tutorial is an introduction to the callback and event handling system in
# DifferentialEquations.jl, documented in the [Event Handling and Callback
# Functions](http://docs.juliadiffeq.org/latest/features/callback_functions.html) page of
# the documentation. We will also introduce you to some of the most widely used callbacks in
# the [Callback Library](http://docs.juliadiffeq.org/latest/features/callback_library.html),
# which is a library of pre-built mods.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Events and Continuous Callbacks
#
# Event handling is done through continuous callbacks. Callbacks take a function,
# `condition`, which triggers an `affect!` when `condition == 0`. These callbacks are called
# "continuous" because they will utilize rootfinding on the interpolation to find the
# "exact" time point at which the condition takes place and apply the `affect!` at that time
# point.
#
# ***Let's use a bouncing ball as a simple system to explain events and callbacks.*** Let's
# take Newton's model of a ball falling towards the Earth's surface via a gravitational
# constant `g`. In this case, the velocity is changing via `-g`, and position is changing
# via the velocity. Therefore we receive the system of ODEs:
# ------------------------------------------------------------------------------------------

using DifferentialEquations
ball! = @ode_def BallBounce begin
  dy =  v
  dv = -g
end g

# ------------------------------------------------------------------------------------------
# We want the callback to trigger when `y=0` since that's when the ball will hit the Earth's
# surface -- our event. We do this with the condition:
# ------------------------------------------------------------------------------------------

function condition(u,t,integrator) # Event when event_f(u,t) == 0
  u[1]
end

# ------------------------------------------------------------------------------------------
# Recall that the `condition` will trigger when it evaluates to zero, and here it will
# evaluate to zero when `u[1] == 0`, which occurs when `v == 0`. *Now we have to say what we
# want the callback to do.* Callbacks make use of the [Integrator
# Interface](http://docs.juliadiffeq.org/latest/basics/integrator.html). Instead of giving a
# full description, a quick and usable rundown is:
#
# - Values are strored in `integrator.u`
# - Times are stored in `integrator.t`
# - The parameters are stored in `integrator.p`
# - `integrator(t)` performs an interpolation in the current interval between
# `integrator.tprev` and `integrator.t` (and allows extrapolation)
# - User-defined options (tolerances, etc.) are stored in `integrator.opts`
# - `integrator.sol` is the current solution object. Note that `integrator.sol.prob` is the
# current problem
#
# While there's a lot more on the integrator interface page, that's a working knowledge of
# what's there.
#
# What we want to do with our `affect!` is to "make the ball bounce". Mathematically
# speaking, the ball bounces when the sign of the velocity flips. As an added behavior,
# let's also use a small friction constant to dampen the ball's velocity. This way only a
# percentage of the velocity will be retained when the event is triggered and the callback
# is used.  We'll define this behavior in the `affect!` function:
# ------------------------------------------------------------------------------------------

function affect!(integrator)
    integrator.u[2] = -integrator.p[2] * integrator.u[2]
end

# ------------------------------------------------------------------------------------------
# `integrator.u[2]` is the second value of our model, which is `v` or velocity, and
# `integrator.p[2]`, is our friction coefficient.
#
# Therefore `affect!` can be read as follows: `affect!` will take the current value of
# velocity, and multiply it `-1` multiplied by our friction coefficient. Therefore the ball
# will change direction and its velocity will dampen when `affect!` is called.
#
# Now let's build the `ContinuousCallback`:
# ------------------------------------------------------------------------------------------

bounce_cb = ContinuousCallback(condition,affect!)

# ------------------------------------------------------------------------------------------
# Now let's make an `ODEProblem` which has our callback:
# ------------------------------------------------------------------------------------------

u0 = [50.0,0.0]
tspan = (0.0,15.0)
p = (9.8,0.9)
prob = ODEProblem(ball!,u0,tspan,p,callback=bounce_cb)

# ------------------------------------------------------------------------------------------
# Notice that we chose a friction constant of `0.9`. Now we can solve the problem and plot
# the solution as we normally would:
# ------------------------------------------------------------------------------------------

sol = solve(prob,Tsit5())
using Plots; gr()
plot(sol)

# ------------------------------------------------------------------------------------------
# and tada, the ball bounces! Notice that the `ContinuousCallback` is using the
# interpolation to apply the effect "exactly" when `v == 0`. This is crucial for model
# correctness, and thus when this property is needed a `ContinuousCallback` should be used.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### Exercise 1
#
# In our example we used a constant coefficient of friction, but if we are bouncing the ball
# in the same place we may be smoothing the surface (say, squishing the grass), causing
# there to be less friction after each bounce. In this more advanced model, we want the
# friction coefficient at the next bounce to be `sqrt(friction)` from the previous bounce
# (since `friction < 1`, `sqrt(friction) > friction` and `sqrt(friction) < 1`).
#
# Hint: there are many ways to implement this. One way to do it is to make `p` a `Vector`
# and mutate the friction coefficient in the `affect!`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Discrete Callbacks
#
# A discrete callback checks a `condition` after every integration step and, if true, it
# will apply an `affect!`. For example, let's say that at time `t=2` we want to include that
# a kid kicked the ball, adding `20` to the current velocity. This kind of situation, where
# we want to add a specific behavior which does not require rootfinding, is a good candidate
# for a `DiscreteCallback`. In this case, the `condition` is a boolean for whether to apply
# the `affect!`, so:
# ------------------------------------------------------------------------------------------

function condition_kick(u,t,integrator)
    t == 2
end

# ------------------------------------------------------------------------------------------
# We want the kick to occur at `t=2`, so we check for that time point. When we are at this
# time point, we want to do:
# ------------------------------------------------------------------------------------------

function affect_kick!(integrator)
    integrator.u[2] += 50
end

# ------------------------------------------------------------------------------------------
# Now we build the problem as before:
# ------------------------------------------------------------------------------------------

kick_cb = DiscreteCallback(condition_kick,affect_kick!)
u0 = [50.0,0.0]
tspan = (0.0,10.0)
p = (9.8,0.9)
prob = ODEProblem(ball!,u0,tspan,p,callback=kick_cb)

# ------------------------------------------------------------------------------------------
# Note that, since we are requiring our effect at exactly the time `t=2`, we need to tell
# the integration scheme to step at exactly `t=2` to apply this callback. This is done via
# the option `tstops`, which is like `saveat` but means "stop at these values".
# ------------------------------------------------------------------------------------------

sol = solve(prob,Tsit5(),tstops=[2.0])
plot(sol)

# ------------------------------------------------------------------------------------------
# Note that this example could've been done with a `ContinuousCallback` by checking the
# condition `t-2`.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Merging Callbacks with Callback Sets
#
# In some cases you may want to merge callbacks to build up more complex behavior. In our
# previous result, notice that the model is unphysical because the ball goes below zero!
# What we really need to do is add the bounce callback together with the kick. This can be
# achieved through the `CallbackSet`.
# ------------------------------------------------------------------------------------------

cb = CallbackSet(bounce_cb,kick_cb)

# ------------------------------------------------------------------------------------------
# A `CallbackSet` merges their behavior together. The logic is as follows. In a given
# interval, if there are multiple continuous callbacks that would trigger, only the one that
# triggers at the earliest time is used. The time is pulled back to where that continuous
# callback is triggered, and then the `DiscreteCallback`s in the callback set are called in
# order.
# ------------------------------------------------------------------------------------------

u0 = [50.0,0.0]
tspan = (0.0,15.0)
p = (9.8,0.9)
prob = ODEProblem(ball!,u0,tspan,p,callback=cb)
sol = solve(prob,Tsit5(),tstops=[2.0])
plot(sol)

# ------------------------------------------------------------------------------------------
# Notice that we have now merged the behaviors. We can then nest this as deep as we like.
#
# #### Exercise 2
#
# Add to the model a linear wind with resistance that changes the acceleration to `-g + k*v`
# after `t=10`. Do so by adding another parameter and allowing it to be zero until a
# specific time point where a third callback triggers the change.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Integration Termination and Directional Handling
#
# Let's look at another model now -- the model of the [Harmonic
# Oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator). We can write this as:
# ------------------------------------------------------------------------------------------

u0 = [1.,0.]
harmonic! = @ode_def HarmonicOscillator begin
   dv = -x
   dx = v
end
tspan = (0.0,10.0)
prob = ODEProblem(harmonic!,u0,tspan)
sol = solve(prob)
plot(sol)

# ------------------------------------------------------------------------------------------
# Let's instead stop the integration when a condition is met. From the [Integrator Interface
# stepping controls](http://docs.juliadiffeq.org/latest/basics/integrator.html#Stepping-
# Controls-1) we see that `terminate!(integrator)` will cause the integration to end. So our
# new `affect!` is simply:
# ------------------------------------------------------------------------------------------

function terminate_affect!(integrator)
    terminate!(integrator)
end

# ------------------------------------------------------------------------------------------
# Let's first stop the integration when the particle moves back to `x=0`. This means we want
# to use the condition:
# ------------------------------------------------------------------------------------------

function terminate_condition(u,t,integrator)
    u[2]
end
terminate_cb = ContinuousCallback(terminate_condition,terminate_affect!)

# ------------------------------------------------------------------------------------------
# Note that instead of adding callbacks to the problem, we can also add them to the `solve`
# command. This will automatically form a `CallbackSet` with any problem-related callbacks
# and naturally allows you to distinguish between model features and integration controls.
# ------------------------------------------------------------------------------------------

sol = solve(prob,callback=terminate_cb)
plot(sol)

# ------------------------------------------------------------------------------------------
# Notice that the harmonic oscilator's true solution here is `sin` and `cosine`, and thus we
# would expect this return to zero to happen at `t=π`:
# ------------------------------------------------------------------------------------------

sol.t[end]

# ------------------------------------------------------------------------------------------
# This is one way to approximate π! Lower tolerances and arbitrary precision numbers can
# make this more exact, but let's not look at that. Instead, what if we wanted to halt the
# integration after exactly one cycle? To do so we would need to ignore the first zero-
# crossing. Luckily in these types of scenarios there's usually a structure to the problem
# that can be exploited. Here, we only want to trigger the `affect!` when crossing from
# positive to negative, and not when crossing from negative to positive. In other words, we
# want our `affect!` to only occur on upcrossings.
#
# If the `ContinuousCallback` constructor is given a single `affect!`, it will occur on both
# upcrossings and downcrossings. If there are two `affect!`s given, then the first is for
# upcrossings and the second is for downcrossings. An `affect!` can be ignored by using
# `nothing`. Together, the "upcrossing-only" version of the effect means that the first
# `affect!` is what we defined above and the second is `nothing`. Therefore we want:
# ------------------------------------------------------------------------------------------

terminate_upcrossing_cb = ContinuousCallback(terminate_condition,terminate_affect!,nothing)

# ------------------------------------------------------------------------------------------
# Which gives us:
# ------------------------------------------------------------------------------------------

sol = solve(prob,callback=terminate_upcrossing_cb)
plot(sol)

# ------------------------------------------------------------------------------------------
# ## Callback Library
#
# As you can see, callbacks can be very useful and through `CallbackSets` we can merge
# together various behaviors. Because of this utility, there is a library of pre-built
# callbacks known as the [Callback
# Library](http://docs.juliadiffeq.org/latest/features/callback_library.html). We will walk
# through a few examples where these callbacks can come in handy.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### Manifold Projection
#
# One callback is the manifold projection callback. Essentially, you can define any manifold
# `g(sol)=0` which the solution must live on, and cause the integration to project to that
# manifold after every step. As an example, let's see what happens if we naively run the
# harmonic oscillator for a long time:
# ------------------------------------------------------------------------------------------

tspan = (0.0,10000.0)
prob = ODEProblem(harmonic!,u0,tspan)
sol = solve(prob)
gr(fmt=:png) # Make it a PNG instead of an SVG since there's a lot of points!
plot(sol,vars=(1,2))

plot(sol,vars=(0,1),denseplot=false)

# ------------------------------------------------------------------------------------------
# Notice that what's going on is that the numerical solution is drifting from the true
# solution over this long time scale. This is because the integrator is not conserving
# energy.
# ------------------------------------------------------------------------------------------

plot(sol.t,[u[2]^2 + u[1]^2 for u in sol.u]) # Energy ~ x^2 + v^2

# ------------------------------------------------------------------------------------------
# Some integration techniques like [symplectic
# integrators](http://docs.juliadiffeq.org/latest/solvers/dynamical_solve.html#Symplectic-
# Integrators-1) are designed to mitigate this issue, but instead let's tackle the problem
# by enforcing conservation of energy. To do so, we define our manifold as the one where
# energy equals 1 (since that holds in the initial condition), that is:
# ------------------------------------------------------------------------------------------

function g(resid,u,p,t)
  resid[1] = u[2]^2 + u[1]^2 - 1
  resid[2] = 0
end

# ------------------------------------------------------------------------------------------
# Here the residual measures how far from our desired energy we are, and the number of
# conditions matches the size of our system (we ignored the second one by making the
# residual 0). Thus we define a `ManifoldProjection` callback and add that to the solver:
# ------------------------------------------------------------------------------------------

cb = ManifoldProjection(g)
sol = solve(prob,callback=cb)
plot(sol,vars=(1,2))

plot(sol,vars=(0,1),denseplot=false)

# ------------------------------------------------------------------------------------------
# Now we have "perfect" energy conservation, where if it's ever violated too much the
# solution will get projected back to `energy=1`.
# ------------------------------------------------------------------------------------------

u1,u2 = sol[500]
u2^2 + u1^2

# ------------------------------------------------------------------------------------------
# While choosing different integration schemes and using lower tolerances can achieve this
# effect as well, this can be a nice way to enforce physical constraints and is thus used in
# many disciplines like molecular dynamics. Another such domain constraining callback is the
# [`PositiveCallback()`](http://docs.juliadiffeq.org/latest/features/callback_library.html#P
# ositiveDomain-1) which can be used to enforce positivity of the variables.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ### SavingCallback
#
# The `SavingCallback` can be used to allow for special saving behavior. Let's take a linear
# ODE define on a system of 1000x1000 matrices:
# ------------------------------------------------------------------------------------------

prob = ODEProblem((du,u,p,t)->du.=u,rand(1000,1000),(0.0,1.0))

# ------------------------------------------------------------------------------------------
# In fields like quantum mechanics you may only want to know specific properties of the
# solution such as the trace or the norm of the matrix. Saving all of the 1000x1000 matrices
# can be a costly way to get this information! Instead, we can use the `SavingCallback` to
# save the `trace` and `norm` at specified times. To do so, we first define our
# `SavedValues` cache. Our time is in terms of `Float64`, and we want to save tuples of
# `Float64`s (one for the `trace` and one for the `norm`), and thus we generate the cache
# as:
# ------------------------------------------------------------------------------------------

saved_values = SavedValues(Float64, Tuple{Float64,Float64})

# ------------------------------------------------------------------------------------------
# Now we define the `SavingCallback` by giving it a function of `(u,p,t,integrator)` that
# returns the values to save, and the cache:
# ------------------------------------------------------------------------------------------

cb = SavingCallback((u,t,integrator)->(trace(u),norm(u)), saved_values)

# ------------------------------------------------------------------------------------------
# Here we take `u` and save `(trace(u),norm(u))`. When we solve with this callback:
# ------------------------------------------------------------------------------------------

sol = solve(prob, Tsit5(), callback=cb, save_everystep=false, save_start=false, save_end = false) # Turn off normal saving

# ------------------------------------------------------------------------------------------
# Our values are stored in our `saved_values` variable:
# ------------------------------------------------------------------------------------------

saved_values.t

saved_values.saveval

# ------------------------------------------------------------------------------------------
# By default this happened only at the solver's steps. But the `SavingCallback` has similar
# controls as the integrator. For example, if we want to save at every `0.1` seconds, we do
# can so using `saveat`:
# ------------------------------------------------------------------------------------------

saved_values = SavedValues(Float64, Tuple{Float64,Float64}) # New cache
cb = SavingCallback((u,t,integrator)->(trace(u),norm(u)), saved_values, saveat = 0.0:0.1:1.0)
sol = solve(prob, Tsit5(), callback=cb, save_everystep=false, save_start=false, save_end = false) # Turn off normal saving

saved_values.t

saved_values.saveval

# ------------------------------------------------------------------------------------------
# #### Exercise 3
#
# Go back to the Harmonic oscillator. Use the `SavingCallback` to save an array for the
# energy over time, and do this both with and without the `ManifoldProjection`. Plot the
# results to see the difference the projection makes.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
using Documenter, ConfigLH

makedocs(
    modules = [ConfigLH],
    format = :html,
    checkdocs = :exports,
    sitename = "ConfigLH.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/hendri54/ConfigLH.jl.git",
)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to the exponential atmospheric model.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/earth/atmospheric_models/expatmosphere
# ==================================================

# Functions: expatmosphere
# -------------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
#   Example 8-4 [1, p. 567]
#
#   Using the exponential atmospheric model, one gets:
#
#       ρ( 747.2119 km ) = 2.1219854 ⋅ 10⁻⁴ kg / m³
#
# Scenario 02
# ===========
#
#   By definition, one gets `ρ(h₀) = ρ₀(h₀)` for all tabulated values.
#
# Scenario 03
# ===========
#
#   Inside every interval, the atmospheric density must be monotonically
#   decreasing.
#
################################################################################

@testset "Function expatmosphere" begin
    # Scenario 01
    # ===========

    @test expatmosphere(747211.9) ≈ 2.1219854e-14 rtol = 1e-8

    # Scenario 02
    # ===========

    for i = 1:length(SatelliteToolbox._expatmosphere_h₀)
        h = SatelliteToolbox._expatmosphere_h₀[i]*1000

        @test expatmosphere(h) == SatelliteToolbox._expatmosphere_ρ₀[i]
    end

    # Scenario 03
    # ===========

    for i = 2:length(SatelliteToolbox._expatmosphere_h₀)
        h₀i = SatelliteToolbox._expatmosphere_h₀[i-1]*1000
        h₀f = SatelliteToolbox._expatmosphere_h₀[i]*1000
        Δ   = 100
        Δh₀ = (h₀f - h₀i)/Δ

        ρk = expatmosphere(h₀i)

        for k = 2:Δ
            ρk₋₁ = ρk
            h    = h₀i + Δh₀*(k-1)
            ρk   = expatmosphere(h)

            @test ρk < ρk₋₁
        end
    end
end


text:
# Test ParticleContainer

using Turing

import Turing: ParticleContainer, weights, resample!, effectiveSampleSize, Trace, Trace, current_trace, VarName, Sampler, consume, produce

if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :n)
  n[] = 0
else
  const n = Ref(0)
end

alg = PG(5, 1)
spl = Turing.Sampler(alg)
dist = Normal(0, 1)

function fpc()
  t = TArray(Float64, 1);
  t[1] = 0;
  while true
    ct = current_trace()
    vn = VarName(gensym(), :x, "[$n]", 1)
    Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
    produce(0)
    vn = VarName(gensym(), :x, "[$n]", 1)
    Turing.assume(spl, dist, vn, ct.vi); n[] += 1;
    t[1] = 1 + t[1]
  end
end

pc = ParticleContainer{Trace}(fpc)

push!(pc, Trace(pc.model))
push!(pc, Trace(pc.model))
push!(pc, Trace(pc.model))

Base.@assert weights(pc)[1] == [1/3, 1/3, 1/3]
Base.@assert weights(pc)[2] ≈ log(3)
Base.@assert pc.logE ≈ log(1)

Base.@assert consume(pc) == log(1)

resample!(pc)
Base.@assert pc.num_particles == length(pc)
Base.@assert weights(pc)[1] == [1/3, 1/3, 1/3]
Base.@assert weights(pc)[2] ≈ log(3)
Base.@assert pc.logE ≈ log(1)
Base.@assert effectiveSampleSize(pc) == 3

Base.@assert consume(pc) ≈ log(1)
resample!(pc)
Base.@assert consume(pc) ≈ log(1)


text:
#=
    http://cplint.eu/example/inference/dice.swinb

    """
    A six-sided die is repeatedly thrown until the outcome is six. on(T,F)
    means that on the Tth throw the face F came out.
    """

    What is the probability that the die lands on face 1 at time 0?

    Cf ~/blog/dice_6_throws.blog
     ~/psi/dice_6_throws.psi
     ~/webppl/dice_6_throws.wppl

     Also see the related problem rolling_dice4.jl where we investigate
     the probability of the length of this (stopping) sequence.

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function dice_6_throws(ix=1,val=1)

    max_len = 100
    function throws(a)
        # Note: If t is a random variable then we got StackOverflowError
        # t ~  DiscreteUniform(1,6) # Don't work!
        # Hence we draw a "plain" random value
        t = rand(DiscreteUniform(1,6))
        if t == 6
            return a
        else
            return throws(vcat(a,t))
        end
    end

    # len is a pure random variable, though
    len ~ DiscreteUniform(0,max_len)
    a = throws(Int16[])
    len = length(a)

    # true ~ Dirac(length(a) >= 1 && a[1] == 1)
    # true ~ Dirac(length(a) >= 2 && a[2] == 1)

    # return len
    # Probability of
    return length(a) >= ix && a[ix] == val
    # return length(a) >= 1 && a[1] == 1
end

function run_model(ix, val)
    model = dice_6_throws(ix,val)

    num_chains = 4

    # chains = sample(model, Prior(), 10_000)

    # chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
    # chains = sample(model, MH(), MCMCThreads(), 10_000, num_chains)

    chains = sample(model, MH(), 10_000)
    # chains = sample(model, MH(), 1_000)

    # chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)

    # chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
    # chains = sample(model, SMC(1000), 10_000)
    # chains = sample(model, IS(), 10_000)
    #
    #chains = sample(model, Gibbs(HMC(0.1,5,:a,:b),PG(15,:p)), 10_000)
    # chains = sample(model, Gibbs(NUTS(1000,0.65,:a,:b),PG(15,:p)), 10_000)
    # chains = sample(model, Gibbs(HMC(0.1,5,:a,:b),SMC(1000,:p)), 10_000) # Nope

    # display(chains)
    # display(plot(chains))
    # show_var_dist_pct(chains,:len,1000)

    println("\nprob return value:")
    genq = generated_quantities(model, chains)
    show_var_dist_pct(genq,1000)

end

for i in 1:10
    println("\nix:$i val:1")
    run_model(i,1)
end

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
function lif(nstate::Array{Float64}, vstart::Float64, vthresh::Float64)
    outputs = zeros(5)
    if (nstate[1] >= vthresh) && nstate[2] < 0.8
        outputs[1] = vstart - nstate[1]
        outputs[2] = 1.0 - nstate[2]
        outputs[5] = 1.0
    else
        outputs[1] = -0.05 * nstate[1] + nstate[5]
        outputs[2] = -0.05 * nstate[2]
    end

    max.(-1.0, min.(1.0, outputs))
end

function izhikevich(nstate::Array{Float64}, vstart::Float64, vthresh::Float64,
                    vscale::Float64; a::Float64=0.02, b::Float64=0.2,
                    c::Float64=-65.0/vscale, d::Float64=2.0/vscale)
    v = nstate[1]
    u = nstate[3]
    I = nstate[5]

    fired = 0.0
    if nstate[1] >= vthresh
        v = c
        u = u + d
        fired = 1.0
    end

    c1 = 0.04 * vscale
    c2 = 140 / vscale
    dv = c1 * v^2 + 5 * v + c2 - u + I
    du = a * (b * v - u)

    outputs = [dv, 0.0, du, 0.0, fired]

    max.(-1.0, min.(1.0, outputs))
end


text:
# DataValueArray is dense and allows fast linear indexing.
import Base: LinearFast

Base.IndexStyle(::Type{<:DataValueArray}) = IndexLinear()

# resolve ambiguity created by the two definitions that follow.
function Base.getindex{T, N}(X::DataValueArray{T, N})
    return X[1]
end

"""
    getindex{T, N}(X::DataValueArray{T, N}, I::Int...)

Retrieve a single entry from a `DataValueArray`. If the value in the entry
designated by `I` is present, then it will be returned wrapped in a
`DataValue{T}` container. If the value is missing, then this method returns
`DataValue{T}()`.
"""
# Extract a scalar element from a `DataValueArray`.
@inline function Base.getindex{T, N}(X::DataValueArray{T, N}, I::Int...)
    if isbits(T)
        ifelse(X.isnull[I...], DataValue{T}(), DataValue{T}(X.values[I...]))
    else
        if X.isnull[I...]
            DataValue{T}()
        else
            DataValue{T}(X.values[I...])
        end
    end
end

"""
    getindex{T, N}(X::DataValueArray{T, N}, I::DataValue{Int}...)

Just as above, with the additional behavior that this method throws an error if
any component of the index `I` is null.
"""
@inline function Base.getindex{T, N}(X::DataValueArray{T, N},
                                     I::DataValue{Int}...)
    any(isnull, I) && throw(NullException())
    values = [ get(i) for i in I ]
    return getindex(X, values...)
end

"""
    setindex!(X::DataValueArray, v::DataValue, I::Int...)

Set the entry of `X` at position `I` equal to a `DataValue` value `v`. If
`v` is null, then only `X.isnull` is updated to indicate that the entry at
index `I` is null. If `v` is not null, then `X.isnull` is updated to indicate
that the entry at index `I` is present and `X.values` is updated to store the
value wrapped in `v`.
"""
# Insert a scalar element from a `DataValueArray` from a `DataValue` value.
@inline function Base.setindex!(X::DataValueArray, v::DataValue, I::Int...)
    if isnull(v)
        X.isnull[I...] = true
    else
        X.isnull[I...] = false
        X.values[I...] = get(v)
    end
    return v
end

"""
    setindex!(X::DataValueArray, v::Any, I::Int...)

Set the entry of `X` at position `I` equal to `v`. This method always updates
`X.isnull` to indicate that the entry at index `I` is present and `X.values`
to store `v` at `I`.
"""
# Insert a scalar element from a `DataValueArray` from a non-DataValue value.
@inline function Base.setindex!(X::DataValueArray, v::Any, I::Int...)
    X.values[I...] = v
    X.isnull[I...] = false
    return v
end

function unsafe_getindex_notnull(X::DataValueArray, I::Int...)
    return DataValue(getindex(X.values, I...))
end

function unsafe_getvalue_notnull(X::DataValueArray, I::Int...)
    return getindex(X.values, I...)
end

function Base.checkindex(::Type{Bool}, inds::AbstractUnitRange, i::DataValue)
    isnull(i) ? throw(NullException()) : checkindex(Bool, inds, get(i))
end

function Base.checkindex{N}(::Type{Bool}, inds::AbstractUnitRange, I::DataValueArray{Bool, N})
    any(isnull, I) && throw(NullException())
    checkindex(Bool, inds, I.values)
end

function Base.checkindex{T<:Real}(::Type{Bool}, inds::AbstractUnitRange, I::DataValueArray{T})
    any(isnull, I) && throw(NullException())
    b = true
    for i in 1:length(I)
        @inbounds v = unsafe_getvalue_notnull(I, i)
        b &= checkindex(Bool, inds, v)
    end
    return b
end

function Base.to_index(X::DataValueArray)
    any(isnull, X) && throw(NullException())
    Base.to_index(X.values)
end

"""
    nullify!(X::DataValueArray, I...)

This is a convenience method to set the entry of `X` at index `I` to be null
"""
function nullify!(X::DataValueArray, I...)
    setindex!(X, DataValue{eltype(X)}(), I...)
end


text:
abstract type 状态 end
id(s::状态)=String(typeof(s).name.name)
exist(::状态)=nothing
e_exist(::状态)=nothing
c_union(::状态)=true
function Base.union(s1::T,::T,t1::UInt,t2::UInt)where T<:状态
	return Pair(max(t1,t2),s1)
end
upd(::状态,::Entity)=nothing

mutable struct 状态组
	data::Dict{String,Union{Pair{UInt,状态},Vector{Pair{UInt,状态}}}}
end
function Base.push!(ss::状态组,t::UInt,s::状态)
	sid=id(s)
	if haskey(ss.data,sid)
		if c_union(s)
			@inbounds pa=ss.data[sid]
			@inbounds ss.data[sid]=union(pa.second,s,pa.first,t)
		else
			push!(ss.data[sid]::Vector{Pair},Pair(t,s))
		end
	else
		if c_union(s)
			ss.data[sid]=Pair(t,s)
		else
			ss.data[sid]=[Pair(t,s)]
		end
	end
end

struct 自发移动
	x::Float
	y::Float
end
upd(s::自发移动,t::Entity)=move(t,s.x,s.y)
c_union(::自发移动)=false

abstract type 中毒<:状态 end


text:
module ConstrainedDynamics

using LinearAlgebra
using StaticArrays
using ForwardDiff
using StaticArrays: SUnitRange
using Rotations
using Rotations: RotationError, pure_quaternion, params, lmult, rmult, tmat, vmat, hmat, skew
using Colors: RGBA, RGB
using LightXML

using DocStringExtensions

export Origin,
    Body,
    EqualityConstraint,
    InequalityConstraint,
    Mechanism,
    LinearMechanism,
    Controller,
    Storage,
    UnitQuaternion,

    Box,
    Cylinder,
    Sphere,
    Pyramid,
    Mesh,

    Floating,
    Prismatic,
    Spherical,
    Cylindrical,
    Revolute,
    Planar,
    PlanarFree,
    Fixed,
    FixedOrientation,
    CylindricalFree,

    Impact,
    Friction,
    UnitQuaternion,

    setPosition!,
    setVelocity!,
    setForce!,
    addForce!,
    getid,
    getcomponent,
    getbody,
    geteqconstraint,
    getineqconstraint,
    simulate!,
    initializeConstraints!,
    disassemble,
    minimalCoordinates,
    minimalVelocities,
    activate!,
    deactivate!,
    linearsystem,

    RotX,
    RotY,
    RotZ,
    RGBA,

    szeros,
    sones,
    srand

include(joinpath("util", "util.jl"))
include(joinpath("util", "custom_static.jl"))
include(joinpath("util", "customdict.jl"))
include(joinpath("util", "quaternion.jl"))

include(joinpath("optional_components", "shapes.jl"))
include(joinpath("optional_components", "storage.jl"))

include(joinpath("main_components", "component.jl"))
include(joinpath("main_components", "state.jl"))
include(joinpath("main_components", "body.jl"))
include(joinpath("main_components", "equalityconstraint.jl"))
include(joinpath("main_components", "inequalityconstraint.jl"))
include(joinpath("main_components", "graph.jl"))
include(joinpath("main_components", "controller.jl"))
include(joinpath("main_components", "sparseldu.jl"))
include(joinpath("main_components", "mechanism_struct.jl"))
include(joinpath("main_components", "mechanism_functions.jl"))

include(joinpath("joints", "abstract_joint.jl"))

include(joinpath("joints", "joint.jl"))
include(joinpath("joints", "translational.jl"))
include(joinpath("joints", "rotational.jl"))
include(joinpath("joints", "genericjoint.jl"))
include(joinpath("joints", "prototypes.jl"))

include(joinpath("bounds", "bound.jl"))
include(joinpath("bounds", "contact.jl"))
include(joinpath("bounds", "impact.jl"))
include(joinpath("bounds", "friction.jl"))

include(joinpath("solver", "solverfunctions.jl"))
include(joinpath("solver", "system.jl"))
include(joinpath("solver", "initconstraints.jl"))
include(joinpath("solver", "newton.jl"))
include(joinpath("solver", "linesearch.jl"))

include(joinpath("discretization", "SymplecticEuler.jl"))
# include(joinpath("discretization", "ImplicitTrapezoid.jl"))

include(joinpath("ui", "mechanism_ui.jl"))
include(joinpath("ui", "simulate.jl"))
include(joinpath("ui", "initialize.jl"))
include(joinpath("ui", "urdf.jl"))

end


text:
function all_1e_ints(bfs::BasisSet,mol::Molecule)
  n = length(bfs.bfs)
  S = Array(Float64,(n,n))
  T = Array(Float64,(n,n))
  V = Array(Float64,(n,n))
  for (i,j) in pairs(n)
    a,b = bfs.bfs[i],bfs.bfs[j]
    S[i,j] = S[j,i] = overlap(a,b)
    T[i,j] = T[j,i] = kinetic(a,b)
    V[i,j] = V[j,i] = nuclear_attraction(a,b,mol)
  end
  return S,T,V
end

function all_twoe_ints(bflist,ERI=coulomb)
  n = length(bflist.bfs)
  totlen = div(n*(n+1)*(n*n+n+2),8)
  ints2e = Array(Float64,totlen)
  for (i,j,k,l) in iiterator(n)
    ints2e[iindex(i,j,k,l)] = ERI(bflist.bfs[i],bflist.bfs[j],bflist.bfs[k],bflist.bfs[l])
  end
  return ints2e
end

function make2JmK(D::Array{Float64,2},Ints::Array{Float64,1})
  n = size(D,1)
  G = Array(Float64,(n,n))
  D1 = reshape(D,n*n)
  temp = Array(Float64,n*n)
  for (i,j) in pairs(n)
    kl = 1
    for (k,l) in pairs(n,"rect")
      temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]
      kl += 1
    end
    G[i,j] = G[j,i] = dot(D1,temp)
  end
  return G
end

dmat(U::Array{Float64,2},nocc::Int64) = U[:,1:nocc]*U[:,1:nocc]'

function rhf(mol::Molecule,MaxIter::Int64=8,verbose::Bool=false)
  @printf("Starting RHF calculation:")

  @printf("Building basis ...")
  bfs = build_basis(mol)
  @printf("... done\n\n")

  @printf("1e integrals ...")
  S,T,V = all_1e_ints(bfs,mol)
  @printf("... done\n\n")

  @printf("2e integrals ...")
  Ints = all_twoe_ints(bfs)
  @printf("... done\n\n")

  h = T+V
  E,U = eig(h,S)
  Enuke = nuclear_repulsion(mol)
  nclosed,nopen = divrem(nel(mol),2)
  Eold = 0
  Energy = 0
  println("Nel=$(nel(mol)) Nclosed=$nclosed")
  if verbose
    println("S=\n$S")
    println("h=\n$h")
    println("T=\n$T")
    println("V=\n$V")
    println("E: $E")
    println("U: $U")
    println("2e ints:\n$Ints")
  end
  for iter in 1:MaxIter
    D = dmat(U,nclosed)
    if verbose
      println("D=\n$D")
    end
    G = make2JmK(D,Ints)
    H = h+G
    E,U = eig(H,S)
    Eone = trace2(D,h)
    Etwo = trace2(D,H)
    Energy = Enuke + Eone + Etwo
    println("HF: $iter  $Energy : $Enuke  $Eone  $Etwo")
    if isapprox(Energy,Eold)
      break
    end
    Eold  = Energy
  end
  return Energy,E,U
end


text:
using GADM
using Test
using DataDeps
using ArchGDAL
using GeoInterface
using Tables

testfiles = [
    "GADM.jl",
]

@testset "GADM.jl" begin
  for testfile in testfiles
    include(testfile)
  end
end


text:
## First program Hello World with julia

println("*******************")
println("**** Hola mundo ***")
println("*******************")

println(PROGRAM_FILE)
for x in ARGS
    println(x)
end

for T in [Int8,Int16,Int32,Int64,Int128,UInt8,UInt16,UInt32,UInt64,UInt128]
    println("$(lpad(T,7)): [$(typemin(T)),$(typemax(T))]")
end

v(x) = (println(x); x)

max_stars_repo_name: MadaiA/JuliaProjects

max_issues_repo_name: MadaiA/JuliaProjects

max_forks_repo_name: MadaiA/JuliaProjects


text:
@testset "Lattices" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x^2-2
  K, a = number_field(f)
  D = matrix(K, 3, 3, [1//64, 0, 0, 0, 1//64, 0, 0, 0, 1//64])
  gens = [[32, 0, 0], [944*a+704, 0, 0], [16, 16, 0], [72*a+96, 72*a+96, 0], [4*a, 4*a+8, 8], [20*a+32, 52*a+72, 32*a+40]]
  L = @inferred quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  D = matrix(K, 3, 3, [1//64, 0, 0, 0, 1//64, 0, 0, 0, 1//64])
  gens = [[32, 0, 0], [720*a+448, 0, 0], [16, 16, 0], [152*a+208, 152*a+208, 0], [4*a+24, 4*a, 8], [116*a+152, 20*a+32, 32*a+40]]
  M = @inferred quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  p = prime_decomposition(base_ring(L), 2)[1][1]
  @test @inferred islocally_isometric(L, M, p)

  # Smoke test for genus symbol
  Qx, x = PolynomialRing(FlintQQ, "x")
  K, a = NumberField(x^2 - 2, "a")
  L = @inferred quadratic_lattice(K, identity_matrix(K, 10), gram_ambient_space = 35*identity_matrix(K, 10))
  P = prime_decomposition(maximal_order(K), 5)[1][1]
  #GM = @inferred Hecke._genus_symbol_kirschmer(L, P)
  #@test GM._genus_symbol_kirschmer(L, P).data == [(10, 1, 1)]

  #P = prime_decomposition(maximal_order(K), 2)[1][1]
  #GS = @inferred Hecke._genus_symbol_kirschmer(L, P)
  #@test GS.data[1] == [10]
  #@test GS.data[2] == [0]
  #@test GS.data[3] == [2]
  #@test GS.data[4] == [35]
  #@test GS.data[5] == []

  Qx, x = FlintQQ["x"];
  K, a = number_field(x^2 - 2, cached = true);
  Kt, t = K["t"];
  L, b = number_field(t^2 + 11, "b", check = true)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  P = prime_decomposition(maximal_order(L), p)[1][1]
  H = @inferred Hecke.lattice(hermitian_space(K, 2 * identity_matrix(K, 3)), pseudo_matrix(identity_matrix(K, 3), [p, p, p]))
  #@test Hecke._genus_symbol_kirschmer(H, fmpz(2)) == Any[(3, 4, true, 4, -64)]
  @test @inferred islocally_isometric(H, H, fmpz(2))

  H = Hecke.lattice(hermitian_space(L, L(elem_in_nf(uniformizer(p))) * identity_matrix(L, 3)), pseudo_matrix(identity_matrix(L, 3), [P, P, P]))
  #@test Hecke._genus_symbol_kirschmer(H, p) == Any[(3, 3, false)]
  @test @inferred islocally_isometric(H, H, p)
  GH = @inferred genus(H, p)
  GH2 = @inferred genus(HermLat, L, p, [(3, 3, -1)])
  @test GH == GH2

  V = @inferred hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  GM = @inferred genus(M, p)
  GM2 = @inferred genus(HermLat, L, p, [(0, 2, 1), (1, 1, -1)])
  @test GM == GM2
  q = prime_decomposition(maximal_order(K), 11)[1][1]
  @test Hecke.genus(M, q) == genus(HermLat, L, q, [(-1, 2, 1), (0, 1, 1)])

  L, _ = number_field(t^2 - gen(K) * t + 1)
  V = hermitian_space(L, L(a) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, -1, 0)])

  V = hermitian_space(L, L(10) * identity_matrix(L, 3))
  M = @inferred Hecke.maximal_integral_lattice(V)
  @test Hecke.genus(M, p) == genus(HermLat, L, p, [(-2, 2, 1, 0), (0, 1, 1, 0)])

  Qx, x = QQ["x"]
  f = x^3-39*x-65
  K, a = CyclotomicRealSubfield(8, "a")
  Kt, t = K["t"]
  E, b = number_field(t^2 - a * t + 1, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  P = prime_decomposition(maximal_order(E), p)[1][1]
  ME = maximal_order(E)
  pm = pseudo_matrix(matrix(E, 3, 3, [1, 0, 0, b, 1, 0, 0, 0, 1]), [one(E)*ME, inv(P)^2, one(E)*ME])
  V = @inferred hermitian_space(E, identity_matrix(E, 3))
  L = @inferred lattice(V, pm)
  o = @inferred Hecke.automorphism_group_order(L)
  @test o == 1536
end

@testset "Misc" begin
  Qx, x = FlintQQ["x"]
  K, a = NumberField(x - 1, "a")
  Kt, t = K["t"]
  E, b = NumberField(t^2 - 2, "b")
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  u = @inferred Hecke._non_norm_rep(E, K, p)
  @test parent(u) === K
  @test @inferred !islocal_norm(E, u, p)
  @test valuation(u - 1, p) == @inferred normic_defect(E, u, p)
end

@testset "Jordan decomposition" begin
  Qx, x = PolynomialRing(FlintQQ, "x", cached = false)
  f = x - 1;
  K, a = number_field(f)
  D = matrix(K, 3, 3, [3, 2, 1, 2, 3, 1, 1, 1, 1]);
  gens = [[1, -1, 0], [1, -1, 0], [0, 1, -1], [0, 1, -1]]
  L = quadratic_lattice(K, generators = gens, gram_ambient_space = D)
  p = prime_decomposition(maximal_order(K), 2)[1][1]
  B, B, S = jordan_decomposition(L, p)
  @test length(S) == 1
end


text:
immutable MixtureModel{VF,VS,Component<:Distribution} <: Distribution{VF,VS}
    components::Vector{Component}
    probs::Vector{Float64}
    aliastable::AliasTable
    function MixtureModel(c::Vector{Component}, p::Vector{Float64})
        if !(Component <: Distribution{VF,VS})
            throw(TypeError(:MixtureModel,
                            "Mixture components type mismatch",
                            Distribution{VF,VS},
                            Component))
        end
        if length(c) != length(p)
            throw(ArgumentError(string("components and probs must have",
                                       " the same number of elements")))
        end
        sizes = [size(component)::Tuple for component in c]
        if !all(sizes .== sizes[1])
            error("MixtureModel: mixture components have different dimensions")
        end
        sump = 0.0
        for i in 1:length(p)
            if p[i] < 0.0
                error("MixtureModel: probabilities must be non-negative")
            end
            sump += p[i]
        end
        table = AliasTable(p ./ sump)
        new(c, p ./ sump, table)
    end
end

function MixtureModel{C<:Distribution}(cs::Vector{C}, ps::Vector{Float64})
    VF = variate_form(C)
    VS = value_support(C)
    MixtureModel{VF,VS,C}(cs, ps)
end

dim(d::MixtureModel{Multivariate}) = dim(d.components[1])

function mean(d::MixtureModel)
    m = 0.0
    for i in 1:length(d.components)
        m += mean(d.components[i]) * d.probs[i]
    end
    return m
end

function _mixpdf(d::MixtureModel, x)
    p = 0.0
    for i in 1:length(d.components)
        p += pdf(d.components[i], x) * d.probs[i]
    end
    return p
end

function _mixlogpdf(d::MixtureModel, x)
    l = [(logpdf(d.components[i],x)+log(d.probs[i]))::Float64
         for i in find(d.probs .> 0.)]
    m = maximum(l)
    log(sum(exp(l.-m))) + m
end

pdf(d::MixtureModel{Univariate}, x::Real) = _mixpdf(d, x)
logpdf(d::MixtureModel{Univariate}, x::Real) = _mixlogpdf(d, x)

_pdf(d::MixtureModel{Multivariate}, x::AbstractVector) = _mixpdf(d, x)
_logpdf(d::MixtureModel{Multivariate}, x::AbstractVector) = _mixpdf(d, x)

function rand(d::MixtureModel)
    i = rand(d.aliastable)
    return rand(d.components[i])
end

function var(d::MixtureModel)
    m = 0.0
    squared_mean_mixture = mean(d).^2
    for i in 1:length(d.components)
        m += (var(d.components[i]) .- squared_mean_mixture .+ mean(d.components[i]).^2) * d.probs[i]
    end
    return m
end

size(d::MixtureModel) = size(d.components[1])
length(d::MixtureModel) = length(d.components[1])


text:
function springpot_spring_reduce_G()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Ga(t, [1.5, 0.0])
    spring = Spring.Ga(t, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(t))
end
@test springpot_spring_reduce_G()

function springpot_spring_reduce_J()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Ja(t, [1.5, 0.0])
    spring = Spring.Ja(t, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(t))
end
@test springpot_spring_reduce_J()

function springpot_spring_reduce_Gp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Gpa(ω, [1.5, 0.0])
    spring = Spring.Gpa(ω, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(ω))
end
@test springpot_spring_reduce_Gp()

function springpot_spring_reduce_Gpp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Gppa(ω, [1.5, 0.0])
    spring = Spring.Gppa(ω, [1.5])

    all(i -> isapprox(springpot[i], spring[i]), eachindex(ω))
end
@test springpot_spring_reduce_Gpp()

function springpot_dashpot_reduce_G()
    # This gives NaN as defined presently,
    # not test written
    return true
end
@test springpot_dashpot_reduce_G()

function springpot_dashpot_reduce_J()
    t = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Ja(t, [1.5, 1.0])
    dashpot = Dashpot.Ja(t, [1.5])

    all(i -> isapprox(springpot[i], dashpot[i]), eachindex(t))
end
@test springpot_dashpot_reduce_J()

function springpot_dashpot_reduce_Gp(tol)
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Gpa(ω, [1.5, 1.0])
    dashpot = Dashpot.Gpa(ω, [1.5])

    # lower tolerance due to cos(0) innacuracy in Julia
    all(i -> isapprox(springpot[i], dashpot[i], atol=tol), eachindex(ω))
end
@test springpot_dashpot_reduce_Gp(tol)

function springpot_dashpot_reduce_Gpp()
    ω = Vector{RheoFloat}(0.0:0.1:10.0)

    springpot = Springpot.Gppa(ω, [1.5, 1.0])
    dashpot = Dashpot.Gppa(ω, [1.5])

    all(i -> isapprox(springpot[i], dashpot[i]), eachindex(ω))
end
@test springpot_dashpot_reduce_Gpp()


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

module LP

import ..FileFormats
import MathOptInterface
const MOI = MathOptInterface

# Julia 1.6 removes Grisu from Base. Previously, we went
#   _print_shortest(io, x) = Base.Grisu.print_shortest(io, x)
# To avoid adding Grisu as a dependency, use the following printing heuristic.
# TODO(odow): consider printing 1.0 as 1.0 instead of 1, i.e., without the
# rounding branch.
function _print_shortest(io::IO, x::Real)
    x_int = round(Int, x)
    if isapprox(x, x_int)
        print(io, x_int)
    else
        print(io, x)
    end
    return
end

MOI.Utilities.@model(
    Model,
    (MOI.ZeroOne, MOI.Integer),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan, MOI.Interval),
    (),
    (MOI.SOS1, MOI.SOS2),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    ()
)
function MOI.supports(
    ::Model{T},
    ::MOI.ObjectiveFunction{<:MOI.ScalarQuadraticFunction{T}},
) where {T}
    return false
end

struct Options
    maximum_length::Int
    warn::Bool
end

function get_options(m::Model)
    default_options = Options(255, false)
    return get(m.ext, :LP_OPTIONS, default_options)
end

"""
    Model(; kwargs...)

Create an empty instance of FileFormats.LP.Model.

Keyword arguments are:

 - `maximum_length::Int=255`: the maximum length for the name of a variable.
   lp_solve 5.0 allows only 16 characters, while CPLEX 12.5+ allow 255.

 - `warn::Bool=false`: print a warning when variables or constraints are renamed.
"""
function Model(; maximum_length::Int = 255, warn::Bool = false)
    model = Model{Float64}()
    options = Options(maximum_length, warn)
    model.ext[:LP_OPTIONS] = options
    return model
end

function Base.show(io::IO, ::Model)
    print(io, "A .LP-file model")
    return
end

# ==============================================================================
#
#   Base.write
#
# ==============================================================================

const _START_REG = r"^([\.0-9eE])"
const _NAME_REG = r"([^a-zA-Z0-9\!\"\#\$\%\&\(\)\/\,\.\;\?\@\_\`\'\{\}\|\~])"

function _write_function(
    io::IO,
    ::Model,
    func::MOI.VariableIndex,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, variable_names[func])
    return
end

function _write_function(
    io::IO,
    ::Model,
    func::MOI.ScalarAffineFunction{Float64},
    variable_names::Dict{MOI.VariableIndex,String},
)
    is_first_item = true
    if !(func.constant ≈ 0.0)
        _print_shortest(io, func.constant)
        is_first_item = false
    end
    for term in func.terms
        if !(term.coefficient ≈ 0.0)
            if is_first_item
                _print_shortest(io, term.coefficient)
                is_first_item = false
            else
                print(io, term.coefficient < 0 ? " - " : " + ")
                _print_shortest(io, abs(term.coefficient))
            end
            print(io, " ", variable_names[term.variable])
        end
    end
    return
end

function _write_constraint_suffix(io::IO, set::MOI.LessThan)
    print(io, " <= ")
    _print_shortest(io, set.upper)
    println(io)
    return
end

function _write_constraint_suffix(io::IO, set::MOI.GreaterThan)
    print(io, " >= ")
    _print_shortest(io, set.lower)
    println(io)
    return
end

function _write_constraint_suffix(io::IO, set::MOI.EqualTo)
    print(io, " = ")
    _print_shortest(io, set.value)
    println(io)
    return
end

function _write_constraint_suffix(io::IO, set::MOI.Interval)
    print(io, " <= ")
    _print_shortest(io, set.upper)
    println(io)
    return
end

function _write_constraint_prefix(io::IO, set::MOI.Interval)
    _print_shortest(io, set.lower)
    print(io, " <= ")
    return
end

_write_constraint_prefix(::IO, ::Any) = nothing

function _write_constraint(
    io::IO,
    model::Model,
    index::MOI.ConstraintIndex,
    variable_names::Dict{MOI.VariableIndex,String};
    write_name::Bool = true,
)
    func = MOI.get(model, MOI.ConstraintFunction(), index)
    set = MOI.get(model, MOI.ConstraintSet(), index)
    if write_name
        print(io, MOI.get(model, MOI.ConstraintName(), index), ": ")
    end
    _write_constraint_prefix(io, set)
    _write_function(io, model, func, variable_names)
    _write_constraint_suffix(io, set)
    return
end

const _SCALAR_SETS = (
    MOI.LessThan{Float64},
    MOI.GreaterThan{Float64},
    MOI.EqualTo{Float64},
    MOI.Interval{Float64},
)

function _write_sense(io::IO, model::Model)
    if MOI.get(model, MOI.ObjectiveSense()) == MOI.MAX_SENSE
        println(io, "maximize")
    else
        println(io, "minimize")
    end
    return
end

function _write_objective(
    io::IO,
    model::Model,
    variable_names::Dict{MOI.VariableIndex,String},
)
    print(io, "obj: ")
    F = MOI.get(model, MOI.ObjectiveFunctionType())
    f = MOI.get(model, MOI.ObjectiveFunction{F}())
    _write_function(io, model, f, variable_names)
    println(io)
    return
end

function _write_integrality(
    io::IO,
    model::Model,
    key::String,
    ::Type{S},
    variable_names::Dict{MOI.VariableIndex,String},
) where {S}
    indices = MOI.get(model, MOI.ListOfConstraintIndices{MOI.VariableIndex,S}())
    if length(indices) == 0
        return
    end
    println(io, key)
    for index in indices
        f = MOI.get(model, MOI.ConstraintFunction(), index)
        _write_function(io, model, f, variable_names)
        println(io)
    end
    return
end

function _write_constraints(io, model, S, variable_names)
    F = MOI.ScalarAffineFunction{Float64}
    for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        _write_constraint(io, model, index, variable_names; write_name = true)
    end
    return
end

function _write_bounds(io, model, S, variable_names, free_variables)
    F = MOI.VariableIndex
    for index in MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        delete!(free_variables, MOI.VariableIndex(index.value))
        _write_constraint(io, model, index, variable_names; write_name = false)
    end
    return
end

function _write_sos_constraints(io, model, variable_names)
    T, F = Float64, MOI.VectorOfVariables
    sos1_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,MOI.SOS1{T}}())
    sos2_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,MOI.SOS2{T}}())
    if length(sos1_indices) + length(sos2_indices) == 0
        return
    end
    println(io, "SOS")
    for index in sos1_indices
        _write_constraint(io, model, index, variable_names)
    end
    for index in sos2_indices
        _write_constraint(io, model, index, variable_names)
    end
    return
end

_to_string(::Type{MOI.SOS1{Float64}}) = "S1::"
_to_string(::Type{MOI.SOS2{Float64}}) = "S2::"

function _write_constraint(
    io::IO,
    model::Model,
    index::MOI.ConstraintIndex{MOI.VectorOfVariables,S},
    variable_names::Dict{MOI.VariableIndex,String},
) where {S<:Union{MOI.SOS1{Float64},MOI.SOS2{Float64}}}
    f = MOI.get(model, MOI.ConstraintFunction(), index)
    s = MOI.get(model, MOI.ConstraintSet(), index)
    name = MOI.get(model, MOI.ConstraintName(), index)
    if name !== nothing && !isempty(name)
        print(io, name, ": ")
    end
    print(io, _to_string(S))
    for (w, x) in zip(s.weights, f.variables)
        print(io, " ", variable_names[x], ":", w)
    end
    println(io)
    return
end

"""
    Base.write(io::IO, model::FileFormats.LP.Model)

Write `model` to `io` in the LP file format.
"""
function Base.write(io::IO, model::Model)
    options = get_options(model)
    FileFormats.create_unique_names(
        model,
        warn = options.warn,
        replacements = [
            s -> match(_START_REG, s) !== nothing ? "_" * s : s,
            s -> replace(s, _NAME_REG => "_"),
            s -> s[1:min(length(s), options.maximum_length)],
        ],
    )
    variable_names = Dict{MOI.VariableIndex,String}(
        index => MOI.get(model, MOI.VariableName(), index) for
        index in MOI.get(model, MOI.ListOfVariableIndices())
    )
    free_variables = Set(keys(variable_names))
    _write_sense(io, model)
    _write_objective(io, model, variable_names)
    println(io, "subject to")
    for S in _SCALAR_SETS
        _write_constraints(io, model, S, variable_names)
    end
    println(io, "Bounds")
    for S in _SCALAR_SETS
        _write_bounds(io, model, S, variable_names, free_variables)
    end
    # If a variable is binary, it should not be listed as `free` in the bounds
    # section.
    attr = MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}()
    for index in MOI.get(model, attr)
        delete!(free_variables, MOI.VariableIndex(index.value))
    end
    # By default, variables have bounds of [0, ∞), so we need to explicitly
    # declare variables as free.
    for variable in sort(collect(free_variables), by = x -> x.value)
        println(io, variable_names[variable], " free")
    end
    _write_integrality(io, model, "General", MOI.Integer, variable_names)
    _write_integrality(io, model, "Binary", MOI.ZeroOne, variable_names)
    _write_sos_constraints(io, model, variable_names)
    println(io, "End")
    return
end

# ==============================================================================
#
#   Base.read!
#
# ==============================================================================

const _KW_OBJECTIVE = Val{:objective}()
const _KW_CONSTRAINTS = Val{:constraints}()
const _KW_BOUNDS = Val{:bounds}()
const _KW_INTEGER = Val{:integer}()
const _KW_BINARY = Val{:binary}()
const _KW_SOS = Val{:sos}()
const _KW_END = Val{:end}()

const _KEYWORDS = Dict(
    # _KW_OBJECTIVE
    "max" => _KW_OBJECTIVE,
    "maximize" => _KW_OBJECTIVE,
    "maximise" => _KW_OBJECTIVE,
    "maximum" => _KW_OBJECTIVE,
    "min" => _KW_OBJECTIVE,
    "minimize" => _KW_OBJECTIVE,
    "minimise" => _KW_OBJECTIVE,
    "minimum" => _KW_OBJECTIVE,
    # _KW_CONSTRAINTS
    "subject to" => _KW_CONSTRAINTS,
    "such that" => _KW_CONSTRAINTS,
    "st" => _KW_CONSTRAINTS,
    "s.t." => _KW_CONSTRAINTS,
    # _KW_BOUNDS
    "bounds" => _KW_BOUNDS,
    "bound" => _KW_BOUNDS,
    # _KW_INTEGER
    "gen" => _KW_INTEGER,
    "general" => _KW_INTEGER,
    "generals" => _KW_INTEGER,
    # _KW_BINARY
    "bin" => _KW_BINARY,
    "binary" => _KW_BINARY,
    "binaries" => _KW_BINARY,
    # _KW_SOS
    "sos" => _KW_SOS,
    # _KW_END
    "end" => _KW_END,
)

mutable struct _ReadCache
    objective::MOI.ScalarAffineFunction{Float64}
    constraint_function::MOI.ScalarAffineFunction{Float64}
    constraint_name::String
    num_constraints::Int
    name_to_variable::Dict{String,MOI.VariableIndex}
    has_default_bound::Set{MOI.VariableIndex}
    function _ReadCache()
        return new(
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0),
            "",
            0,
            Dict{String,MOI.VariableIndex}(),
            Set{MOI.VariableIndex}(),
        )
    end
end

function _get_variable_from_name(model::Model, cache::_ReadCache, name::String)
    current_variable = get(cache.name_to_variable, name, nothing)
    if current_variable !== nothing
        return current_variable
    end
    options = get_options(model)
    if length(name) > options.maximum_length
        error("Name exceeds maximum length: $name")
    elseif match(r"^([\.0-9])", name) !== nothing
        error("Name starts with invalid character: $name")
    elseif match(_NAME_REG, name) !== nothing
        error("Name contains with invalid character: $name")
    end
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, name)
    # By default, all variables have a lower bound of 0 unless otherwise
    # specified.
    MOI.add_constraint(model, x, MOI.GreaterThan(0.0))
    push!(cache.has_default_bound, x)
    cache.name_to_variable[name] = x
    return x
end

_tokenize(line::AbstractString) = String.(split(line, " "; keepempty = false))

@enum(_TokenType, _TOKEN_VARIABLE, _TOKEN_COEFFICIENT, _TOKEN_SIGN)

function _parse_token(token::String)
    if token == "+"
        return _TOKEN_SIGN, +1.0
    elseif token == "-"
        return _TOKEN_SIGN, -1.0
    end
    coef = tryparse(Float64, token)
    if coef === nothing
        return _TOKEN_VARIABLE, token
    else
        return _TOKEN_COEFFICIENT, coef
    end
end

function _get_term(token_types, token_values, offset)
    coef = 1.0
    if token_types[offset] == _TOKEN_SIGN
        coef = token_values[offset]
        offset += 1
    end
    if token_types[offset] == _TOKEN_COEFFICIENT
        coef *= token_values[offset]
        offset += 1
    elseif token_types[offset] == _TOKEN_SIGN
        error("Invalid line")
    end
    if offset > length(token_types) || token_types[offset] == _TOKEN_SIGN
        return coef, offset  # It's a standalone constant!
    end
    @assert token_types[offset] == _TOKEN_VARIABLE
    x = MOI.VariableIndex(Int64(token_values[offset]))
    return MOI.ScalarAffineTerm(coef, x), offset + 1
end

function _parse_affine_terms(
    f::MOI.ScalarAffineFunction{Float64},
    model::Model,
    cache::_ReadCache,
    tokens::Vector{String},
)
    N = length(tokens)
    token_types = Vector{_TokenType}(undef, N)
    token_values = Vector{Float64}(undef, N)
    for i in 1:length(tokens)
        token_type, token = _parse_token(tokens[i])
        token_types[i] = token_type
        if token_type in (_TOKEN_SIGN, _TOKEN_COEFFICIENT)
            token_values[i] = token::Float64
        else
            @assert token_type == _TOKEN_VARIABLE
            x = _get_variable_from_name(model, cache, token::String)
            # A cheat for type-stability. Store `Float64` of the variable index!
            token_values[i] = Float64(x.value)
        end
    end
    offset = 1
    while offset <= length(tokens)
        term, offset = _get_term(token_types, token_values, offset)
        if term isa MOI.ScalarAffineTerm{Float64}
            push!(f.terms, term::MOI.ScalarAffineTerm{Float64})
        else
            f.constant += term::Float64
        end
    end
    return
end

# _KW_OBJECTIVE

_set_objective_sense(::Any, ::Model, ::String) = nothing

function _set_objective_sense(
    ::typeof(_KW_OBJECTIVE),
    model::Model,
    sense::String,
)
    if sense in ("max", "maximize", "maximise", "maximum")
        MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    else
        @assert sense in ("min", "minimize", "minimise", "minimum")
        MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    end
    return
end

function _parse_section(
    ::typeof(_KW_OBJECTIVE),
    model::Model,
    cache::_ReadCache,
    line::AbstractString,
)
    if occursin(":", line)  # Strip name of the objective
        line = String(match(r"(.*?)\:(.*)", line)[2])
    end
    tokens = _tokenize(line)
    if length(tokens) == 0
        # Can happen if the name of the objective is on one line and the
        # expression is on the next.
        return
    end
    _parse_affine_terms(cache.objective, model, cache, tokens)
    return
end

# _KW_CONSTRAINTS

function _parse_section(
    ::typeof(_KW_CONSTRAINTS),
    model::Model,
    cache::_ReadCache,
    line::AbstractString,
)
    # SOS constraints should be in their own "SOS" section, but we can also
    # recognize them if they're mixed into the constraint section.
    if match(r" S([1-2])\w*:: ", line) !== nothing
        _parse_section(_KW_SOS, model, cache, line)
        return
    end
    if isempty(cache.constraint_name)
        if occursin(":", line)
            m = match(r"(.*?)\:(.*)", line)
            cache.constraint_name = String(m[1])
            line = String(m[2])
        else
            # Give it a temporary name for now
            cache.constraint_name = "R$(cache.num_constraints)"
        end
    end
    tokens = _tokenize(line)
    if length(tokens) == 0
        # Can happen if the name is on one line and the constraint on the next.
        return
    end
    # This checks if the constaint is finishing on this line.
    constraint_set = nothing
    if length(tokens) >= 2 && tokens[end-1] in ("<", "<=", ">", ">=", "=", "==")
        rhs = parse(Float64, pop!(tokens))
        sym = pop!(tokens)
        constraint_set = if sym in ("<", "<=")
            MOI.LessThan(rhs)
        elseif sym in (">", ">=")
            MOI.GreaterThan(rhs)
        else
            @assert sym in ("=", "==")
            MOI.EqualTo(rhs)
        end
    end
    _parse_affine_terms(cache.constraint_function, model, cache, tokens)
    if constraint_set !== nothing
        c = MOI.add_constraint(model, cache.constraint_function, constraint_set)
        MOI.set(model, MOI.ConstraintName(), c, cache.constraint_name)
        cache.num_constraints += 1
        empty!(cache.constraint_function.terms)
        cache.constraint_function.constant = 0.0
        cache.constraint_name = ""
    end
    return
end

# _KW_BOUNDS

function _parse_float(token::String)
    coef = lowercase(token)
    if coef in ("-inf", "-infinity")
        return -Inf
    elseif coef in ("+inf", "+infinity")
        return Inf
    else
        return parse(Float64, coef)
    end
end

# Yes, the last elements here are really accepted by CPLEX...
_is_less_than(token) = token in ("<=", "<", "=<")
_is_greater_than(token) = token in (">=", ">", "=>")
_is_equal_to(token) = token in ("=", "==")

function _parse_section(
    ::typeof(_KW_BOUNDS),
    model::Model,
    cache::_ReadCache,
    line::AbstractString,
)
    tokens = _tokenize(line)
    if length(tokens) == 2 && tokens[2] == "free"
        x = _get_variable_from_name(model, cache, tokens[1])
        _delete_default_lower_bound_if_present(model, cache, x)
        return
    end
    lb, ub, name = -Inf, Inf, ""
    if length(tokens) == 5
        name = tokens[3]
        if _is_less_than(tokens[2]) && _is_less_than(tokens[4])
            lb = _parse_float(tokens[1])
            ub = _parse_float(tokens[5])
        elseif _is_greater_than(tokens[2]) && _is_greater_than(tokens[4])
            lb = _parse_float(tokens[5])
            ub = _parse_float(tokens[1])
        else
            error("Unable to parse bound due to invalid inequalities: $(line)")
        end
    elseif length(tokens) == 3
        name = tokens[1]
        if _is_less_than(tokens[2])
            ub = _parse_float(tokens[3])
            # LP files have default lower bounds of 0, unless the upper bound is
            # less than 0.
            lb = ub > 0.0 ? 0.0 : -Inf
        elseif _is_greater_than(tokens[2])
            lb = _parse_float(tokens[3])
        elseif _is_equal_to(tokens[2])
            lb = ub = _parse_float(tokens[3])
        else
            error("Unable to parse bound due to invalid inequalities: $(line)")
        end
    else
        error("Unable to parse bound: $(line)")
    end
    x = _get_variable_from_name(model, cache, name)
    if lb == ub
        _delete_default_lower_bound_if_present(model, cache, x)
        MOI.add_constraint(model, x, MOI.EqualTo(lb))
    elseif -Inf < lb < ub < Inf
        _delete_default_lower_bound_if_present(model, cache, x)
        MOI.add_constraint(model, x, MOI.Interval(lb, ub))
    elseif -Inf < lb
        _delete_default_lower_bound_if_present(model, cache, x)
        MOI.add_constraint(model, x, MOI.GreaterThan(lb))
    else
        if ub < 0
            # We only need to delete the default lower bound if the upper bound
            # is less than 0.
            _delete_default_lower_bound_if_present(model, cache, x)
        end
        MOI.add_constraint(model, x, MOI.LessThan(ub))
    end
    return
end

function _delete_default_lower_bound_if_present(model, cache, x)
    if !(x in cache.has_default_bound)
        return
    end
    c = MOI.ConstraintIndex{MOI.VariableIndex,MOI.GreaterThan{Float64}}(x.value)
    MOI.delete(model, c)
    delete!(cache.has_default_bound, x)
    return
end

# _KW_INTEGER

function _parse_section(::typeof(_KW_INTEGER), model, cache, line)
    for token in _tokenize(line)
        x = _get_variable_from_name(model, cache, token)
        MOI.add_constraint(model, x, MOI.Integer())
    end
    return
end

# _KW_BINARY

function _parse_section(::typeof(_KW_BINARY), model, cache, line)
    for token in _tokenize(line)
        x = _get_variable_from_name(model, cache, token)
        MOI.add_constraint(model, x, MOI.ZeroOne())
    end
    return
end

# _KW_SOS

function _parse_section(
    ::typeof(_KW_SOS),
    model::Model,
    cache::_ReadCache,
    line::AbstractString,
)
    # SOS constraints can have all manner of whitespace issues with them.
    # Normalize them here before attempting to do anything else.
    line = replace(line, r"\s+:\s+" => ":")
    line = replace(line, r"\s+::" => "::")
    tokens = _tokenize(line)
    if length(tokens) < 3
        error("Malformed SOS constraint: $(line)")
    end
    name = String(split(tokens[1], ":")[1])
    if tokens[2] == "S1::"
        order = 1
    elseif tokens[2] == "S2::"
        order = 2
    else
        error("SOS of type $(tokens[2]) not recognised")
    end
    variables, weights = MOI.VariableIndex[], Float64[]
    for token in tokens[3:end]
        items = String.(split(token, ":"))
        if length(items) != 2
            error("Invalid sequence: $(token)")
        end
        push!(variables, _get_variable_from_name(model, cache, items[1]))
        push!(weights, parse(Float64, items[2]))
    end
    c_ref = if tokens[2] == "S1::"
        MOI.add_constraint(model, variables, MOI.SOS1(weights))
    else
        @assert tokens[2] == "S2::"
        MOI.add_constraint(model, variables, MOI.SOS2(weights))
    end
    MOI.set(model, MOI.ConstraintName(), c_ref, name)
    return
end

# _KW_END

function _parse_section(
    ::typeof(_KW_END),
    ::Model,
    ::_ReadCache,
    line::AbstractString,
)
    return error("Corrupted LP File. You have the lne $(line) after an end.")
end

function _strip_comment(line::String)
    if occursin("\\", line)
        m = match(r"(.*?)\\(.*)", line)
        return strip(String(m[1]))
    else
        return strip(line)
    end
end

"""
    Base.read!(io::IO, model::FileFormats.LP.Model)

Read `io` in the LP file format and store the result in `model`.

This reader attempts to follow the CPLEX LP format, because others like the
lpsolve version are very...flexible...in how they accept input. Read more about
them here: http://lpsolve.sourceforge.net
"""
function Base.read!(io::IO, model::Model)
    if !MOI.is_empty(model)
        error("Cannot read in file because model is not empty.")
    end
    cache = _ReadCache()
    section = Val{:header}()
    while !eof(io)
        line = _strip_comment(string(readline(io)))
        if isempty(line)
            continue
        end
        lower_line = lowercase(line)
        if haskey(_KEYWORDS, lower_line)
            section = _KEYWORDS[lower_line]
            _set_objective_sense(section, model, lower_line)
            continue
        end
        _parse_section(section, model, cache, line)
    end
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        cache.objective,
    )
    return
end

end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface