text:
# packaged version of align.jl. Run align_prep.jl or prep_full.jl first.

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

import NLopt
import PlotlyJS
using Plots; plotly()

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

# based on one region from pkg_calibrate.jl

r = 1
y_cost = y[cost_inds_set[r]]
U_cost = U_y[cost_inds_set[r]]
P_cost = P_y[cost_inds_set[r]]

local_optim_algorithm = NLopt.LN_BOBYQA
xtol_rel = 1e-9
maxeval = 50
maxtime = Inf
β_optim_algorithm = :GN_DIRECT_L
κ_lb_default = 1e-2
κ_ub_default = 1e2
β_max_iters = 500
β_xtol_rel = 1e-9
β_ftol_rel = 1e-9
β_maxtime = Inf

N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )

U_rad_cost = U_cost .* (2*π)

q, updatedfunc, getshiftfunc, N_vars_set,
run_optim, obj_func_β, E_BLS, κ_BLS, b_BLS, updateβfunc, updateκfunc,
q_β = NMRCalibrate.setupcostnesteddwarp(Es, Bs, As, fs, SW, LS_inds, U_rad_cost,
    y_cost, Δsys_cs, a_setp, b_setp; #κs_β_DOFs, κs_β_orderings;
    w = w,
    optim_algorithm = β_optim_algorithm,
    κ_lb_default = κ_lb_default,
    κ_ub_default = κ_ub_default,
    max_iters = β_max_iters,
    xtol_rel = β_xtol_rel,
    ftol_rel = β_ftol_rel,
    maxtime = β_maxtime)

#
#N_β = sum( NMRCalibrate.getNβ(Bs[n]) for n = 1:length(Bs) )
p_β = zeros(T, N_β) # persistant buffer.

obj_func = pp->NMRCalibrate.costnestedd(U_rad_cost, y_cost, updatedfunc, updateκfunc, pp,
Es, Bs, run_optim, E_BLS, κ_BLS, b_BLS, p_β)

### test.
p_test = rand(N_d)
obj_func(p_test)

#
# optim.
N_starts = 100
prob = MultistartOptimization.MinimizationProblem(obj_func, shift_lb, shift_ub)

local_method = MultistartOptimization.NLoptLocalMethod(; algorithm = local_optim_algorithm,
xtol_rel = xtol_rel,
maxeval = maxeval,
maxtime = maxtime)

multistart_method = MultistartOptimization.TikTak(N_starts)
ret_mo = MultistartOptimization.multistart_minimization(multistart_method,
    local_method, prob)


text:
using Plots

using GCMMesh
using GCMMesh.Mesh
using GCMMesh.Graphs
using GCMMesh.Partition

FT = Float64
# plotting a rectangular mesh
x1min, x1max = FT(0), FT(1)
x2min, x2max = FT(0), FT(1)
n1, n2 = 3, 4
per = (false, false)

msh1 = equispaced_rectangular_mesh(x1min, x2min, x1max, x2max, n1, n2, per)

scatter(msh1.coordinates[:, 1], msh1.coordinates[:, 2], legend = false)

for i = 1:msh1.nfaces
    plot!(
        msh1.coordinates[msh1.face_verts[i, :], 1],
        msh1.coordinates[msh1.face_verts[i, :], 2],
        legend = false,
    )
end

savefig("rectangular_mesh.png")
# plotting a sphere mesh
ne = 4
radius = FT(10)
msh = sphere_mesh(ne, radius, EquiangularSphereWarp())
#msh = sphere_mesh(ne, FT, EquidistantSphereWarp())
#msh = cube_panel_mesh(ne, FT)

scatter(msh.coordinates[:, 1], msh.coordinates[:, 2], msh.coordinates[:, 3], legend = false)

for i = 1:msh.nfaces
    plot!(
        msh.coordinates[msh.face_verts[i, :], 1],
        msh.coordinates[msh.face_verts[i, :], 2],
        msh.coordinates[msh.face_verts[i, :], 3],
        legend = false,
        aspect_ratio = 1,
    )
end

savefig("sphere_mesh.png")

# plotting partitions for a rectangular mesh
FT = Float64

x1min, x1max = FT(0), FT(1)
x2min, x2max = FT(0), FT(1)
limits = (x1min, x2min, x1max, x2max)
n1, n2 = 16, 16 # # of elements along x1 and x2 axis respectively
#n1, n2 = 4, 4 # # of elements along x1 and x2 axis respectively
#npart = 4 # # of partitions
#npart = 5 # # of partitions
#npart = 2 # # of partitions
npart = 3 # # of partitions

mesh = equispaced_rectangular_mesh(limits..., n1, n2, (false, false))

partition = rsb_partition(mesh, npart)
colors = (:red, :blue, :green, :yellow, :brown)

scatter(mesh.coordinates[:, 1], mesh.coordinates[:, 2], legend = false)

for i = 1:mesh.nelems
    verts = mesh.elem_verts[i, :]
    x = mesh.coordinates[verts, 1]
    y = mesh.coordinates[verts, 2]
    color = colors[partition[i]+1]
    shape = Shape(x, y)
    plot!(shape, c = color)
end

#savefig("rectangle_mesh_4_partitions.png")
#savefig("rectangle_mesh_5_partitions.png")
#savefig("rectangle_mesh_2_partitions.png")
savefig("rectangle_mesh_3_partitions.png")

# plotting a sphere mesh
colors = (:red, :blue, :green, :yellow, :brown)
ne = 4
radius = FT(10)
mesh = cube_panel_mesh(ne, FT)

#npart = 4
npart = 2
#npart = 3
partition = rsb_partition(mesh, npart)

elems_per_panel = div(mesh.nelems, 6)

plot()

for i = 1:mesh.nelems
    verts = mesh.elem_verts[i, :]
    x = mesh.coordinates[verts, 1]
    y = mesh.coordinates[verts, 2]
    z = mesh.coordinates[verts, 3]
    color = colors[partition[i]+1]
    panelno = cld(i, elems_per_panel)
    xx, yy = unfold_cube_panel_to_plane(x, y, z, panelno)
    plot!(Shape(xx, yy), c = color, legend = false)
end

savefig("sphere_mesh_2_partitions.png")
#savefig("sphere_mesh_3_partitions.png")
#savefig("sphere_mesh_4_partitions.png")

ext: jl

max_stars_repo_name: CliMA/GCMMesh

max_issues_repo_name: CliMA/GCMMesh

max_forks_repo_name: CliMA/GCMMesh


text:
using ReverseDiffDSGE
using Test

@testset "ReverseDiffDSGE.jl" begin
    # Write your tests here.
end


text:
#
# ERSI Shapefile's index file
#   It contains the same `.shp` 100-bytes header
#   Each index record is 8-bytes long:
#       1. 4-bytes (Int32) offset from the file position 0.
#          the offset is in the unit of 16-bit, i.e. need to multiply 2 to get bytes
#       2. 4-bytes (Int32) content length, without the 8-bytes record header
#          to get total number of bytes: (content length)*2 + 8
#
struct IndexRecord
    offset::Int32
    contentLen::Int32
end

struct IndexHandle
    code::Int32
    length::Int32
    version::Int32
    shapeType::Int32
    MBR::Rect
    zrange::Interval
    mrange::Interval
    indices::Vector{IndexRecord}
end

function Base.read(io::IO,::Type{IndexHandle})
    code = ntoh(read(io,Int32))
    read!(io, Vector{Int32}(undef, 5))
    fileSize = ntoh(read(io,Int32))
    version = read(io,Int32)
    shapeType = read(io,Int32)
    MBR = read(io,Rect)
    zmin = read(io,Float64)
    zmax = read(io,Float64)
    mmin = read(io,Float64)
    mmax = read(io,Float64)
    jltype = SHAPETYPE[shapeType]
    records = Vector{IndexRecord}(undef, 0)
    file = IndexHandle(code,fileSize,version,shapeType,MBR,Interval(zmin,zmax),Interval(mmin,mmax),records)
    while !eof(io)
        push!(records, IndexRecord(ntoh(read(io,Int32)),ntoh(read(io,Int32))))
    end
    file
end


text:
#=
Note that the following styles used values from the matplotlib style library
(https://github.com/matplotlib/matplotlib/tree/master/lib/matplotlib/mpl-data/stylelib):

- ggplot
- fivethirtyeight
- seaborn

=#

struct Cycler
    vals::Vector
end

Base.isempty(c::Cycler) = isempty(c.vals)
Base.length(c::Cycler) = length(c.vals)
Cycler(t::Tuple) = Cycler(collect(t))
Cycler(x::Union{String,Number,Date,Symbol}) = Cycler([x])

function Base.getindex(c::Cycler, ix::Integer)
    n = length(c.vals)
    @inbounds v = c.vals[mod1(ix, n)]
    v
end

function Base.getindex(c::Cycler, ixs::AbstractVector{<:Integer})
    [c[i] for i in ixs]
end

struct Style
    layout::Layout
    global_trace::PlotlyAttribute
    trace::Dict{Symbol,PlotlyAttribute}
end

function Style(;
        color_cycle=[],
        layout=Layout(), global_trace=attr(),
        trace=Dict{Symbol,PlotlyAttribute}(),
    )
    if !isempty(color_cycle)
        msg = """
        `color_cycle` argument deprecated. Set the `marker_color` attrbute
        on `global_trace` instead using a Cycler like this:

        Style(global_trace=attr(marker_color=Cycler($(color_cycle))))
        """
        global_trace[:marker_color] = Cycler(color_cycle)
        warn(msg)
    end
    Style(layout, global_trace, trace)
end

function Style(ps1::Style, ps2::Style)
    la = deepcopy(ps1.layout)
    for (k, v) in ps2.layout.fields
        la[k] = v
    end

    gta = merge(ps1.global_trace, ps2.global_trace)

    ta = deepcopy(ps1.trace)
    for (k, v) in ps2.trace
        ta_k = get(ta, k, attr())
        merge!(ta_k, v)
        ta[k] = ta_k
    end

    Style(layout=la, global_trace=gta, trace=ta)
end

Style(pss::Style...) = foldl(Style, pss[1], pss[2:end])

function Style(base::Style;
        color_cycle=[], layout=Layout(),
        global_trace=attr(), trace=Dict{Symbol,PlotlyAttribute}()
    )
    new_style = Style(
        color_cycle=color_cycle, layout=layout, global_trace=global_trace,
        trace=trace
    )
    Style(base, new_style)
end

function Base.show(io::IO, ::MIME"text/plain", s::Style)
    ctx = IOContext(io, limit=true)
    println(io, "Style with:")

    if !isempty(s.layout)
        print(io, "  - "); show(ctx, MIME"text/plain"(), s.layout)
    end

    if !isempty(s.global_trace)
        print(io, "  - global_trace: ")
        show(ctx, MIME"text/plain"(), s.global_trace)
    end

    if !isempty(s.trace)
        println(io, "  - trace: ")
        for (k, v) in s.trace
            print(io, "    - ", k, ": "); show(ctx, MIME"text/plain"(), v)
        end
    end
end

function ==(s1::Style, s2::Style)
    all(nm -> getfield(s1, nm) == getfield(s2, nm), fieldnames(s1))
end

function ggplot_style()
    axis = attr(showgrid=true, gridcolor="white", linewidth=1.0,
                linecolor="white", titlefont_color="#555555",
                titlefont_size=14, ticks="outside",
                tickcolor="#555555"
                )
    layout = Layout(plot_bgcolor="#E5E5E5",
                    paper_bgcolor="white",
                    font_size=10,
                    xaxis=axis,
                    yaxis=axis,
                    titlefont_size=14)

    gta = attr(marker_line_width=0.5, marker_line_color="#348ABD")

    colors = ["#E24A33", "#348ABD", "#988ED5", "#777777", "#FBC15E",
              "#8EBA42", "#FFB5B8"]
    gta[:marker_color] = Cycler(colors)
    Style(layout=layout, global_trace=gta)
end

function fivethirtyeight_style()
    ta = Dict(:scatter=>attr(line_width=4))
    axis = attr(showgrid=true, gridcolor="#cbcbcb",
                linewidth=1.0, linecolor="#f0f0f0",
                ticklen=0.0, tickcolor="#555555", ticks="outside",
                titlefont_size=12, titlefont_color="#555555")
    layout = Layout(plot_bgcolor="#f0f0f0",
                    paper_bgcolor="#f0f0f0",
                    font_size=14,
                    xaxis=axis,
                    yaxis=axis,
                    legend=attr(borderwidth=1.0,
                                bgcolor="f0f0f0", bordercolor="f0f0f0"),
                    titlefont_size=14)
    colors = ["#008fd5", "#fc4f30", "#e5ae38", "#6d904f",
              "#8b8b8b", "#810f7c"]
    gta = attr(marker_color=colors)
    Style(layout=layout, trace=ta, global_trace=gta)
end

function seaborn_style()
    ta = Dict(:heatmap=>attr(colorscale="Greys"),
              :scatter=>attr(marker=attr(size=7, line_width=0.0),
                             line_width=1.75))
    axis = attr(showgrid=true, gridcolor="white",
                linewidth=1.0, linecolor="white",
                ticklen=0.0, tickcolor="#555555", ticks="outside", tickfont_size=10,
                titlefont_size=12, titlefont_color="#555555")
    # TODO: no concept of major vs minor ticks...
    layout = Layout(plot_bgcolor="EAEAF2",
                    paper_bgcolor="white",
                    width=800,
                    height=550, # TODO: what does font_color=0.15 mean??
                    font=attr(family="Arial", size=14, color=0.15),
                    xaxis=axis,
                    yaxis=axis,
                    legend=attr(font_size=10,
                                bgcolor="white", bordercolor="white"),
                    titlefont_size=14)
    colors = ["#4C72B0", "#55A868", "#C44E52", "#8172B2", "#CCB974", "#64B5CD"]
    gta = attr(marker_color=colors)
    Style(trace=ta, layout=layout, global_trace=gta)
end

# This theme was taken from here:
# https://github.com/dcjones/Gadfly.jl/blob/cb28d6aca6b031d01e44146799e520b8bb0d349b/src/theme.jl#L342-L409
function gadfly_dark_style()
    label_color = "#a1a1a1"
    bgcolor = "#222831"
    grid_color = "#575757"

    color_cycle = ["#FE4365", "#ECA25C", "#3F9778", "#EEDBFF", "#236EAD",
                   "#60F6FF", "#D4EC9F", "#7E674B", "#9E7EC1", "#7CB5FB"]

    axis = attr(showgrid=true, gridcolor=grid_color, gridwidth=0.35,
                linecolor=grid_color, titlefont_color=label_color,
                linewidth=1.2, titlefont_size=14, tickcolor=label_color
                )

    layout = Layout(plot_bgcolor=bgcolor,
                    paper_bgcolor=bgcolor,
                    font_size=10,
                    xaxis=axis,
                    yaxis=axis,
                    font_color=label_color,
                    titlefont_size=14,
                    margin=attr(l=40, r=10, t=10, b=30))

    gta = attr(marker_color=color_cycle)
    Style(layout=layout, global_trace=gta)
end

function tomorrow_night_eighties_style()

    bgcolor = "#2d2d2d"  # Background
    grid_color = "#515151"  # Selection
    label_color = "#cccccc"  # Comment
    color_cycle = [
                    "#cc99cc",
                    "#66cccc",
                    "#f2777a",
                    "#ffcc66",
                    "#99cc99",
                    "#f99157",
                    "#6699cc",
                   ]

    axis = attr(showgrid=true, gridcolor=grid_color, gridwidth=0.35,
                linecolor=grid_color, titlefont_color=label_color,
                linewidth=1.2, titlefont_size=14, tickcolor=label_color
                )

    layout = Layout(plot_bgcolor=bgcolor,
                    paper_bgcolor=bgcolor,
                    font_size=10,
                    xaxis=axis,
                    yaxis=axis,
                    font_color=label_color,
                    titlefont_size=14,
                    margin=attr(l=65, r=65, t=65, b=65))

    gta = attr(marker_color=color_cycle)
    Style(layout=layout, global_trace=gta)
end

function style(sty::Symbol)
    sty == :ggplot ? ggplot_style() :
    sty == :fivethirtyeight ? fivethirtyeight_style() :
    sty == :seaborn ? seaborn_style() :
    sty == :gadfly_dark ? gadfly_dark_style() :
    sty == :tomorrow_night_eighties ? tomorrow_night_eighties_style() :
    sty == :default ? DEFAULT_STYLE[] :
    error("Unknown style $sty")
end

const STYLES = [
    :default, :ggplot, :fivethirtyeight, :seaborn, :gadfly_dark,
    :tomorrow_night_eighties
]

# NOTE: DEFAULT_STYLE is reset in __init__() based on environment vars
const DEFAULT_STYLE = Ref{Style}(Style())
const CURRENT_STYLE = Ref{Style}(DEFAULT_STYLE[])
reset_style!() = CURRENT_STYLE[] = DEFAULT_STYLE[]
use_style!(sty::Symbol) = CURRENT_STYLE[] = style(sty)
use_style!(s::Style) = CURRENT_STYLE[] = s


text:
# Model
include_model("box")

# Horizon
T = 10

# Time step
tf = 1.0
h = tf / (T-1)

# ul <= u <= uu
_uu = Inf * ones(model.m)
_ul = zeros(model.m)
_ul[model.idx_u] .= -Inf
ul, uu = control_bounds(model, T, _ul, _uu)

# Initial and final states
mrp_init = MRP(UnitQuaternion(RotY(0.0) * RotX(0.0)))
mrp_side = MRP(UnitQuaternion(RotY(0.0) * RotX(pi / 4.0)))
mrp_corner = MRP(UnitQuaternion(RotY(-1.0 * atan(1.0 / sqrt(2.0))) * RotX(pi / 4.0)))

q1 = [model.r, model.r, model.r, mrp_init.x, mrp_init.y, mrp_init.z]
qT = [model.r, 0.0, model.r * sqrt(2.0), mrp_side.x, mrp_side.y, mrp_side.z]

# qT = [0.0, 0.0, model.r * sqrt(3.0), mrp_corner.x, mrp_corner.y, mrp_corner.z]

x1 = [q1; q1]
xT = [qT; qT]

xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj_penalty = PenaltyObjective(1000.0, model.m)

Qq = Diagonal(0.0 * ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 0.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-3 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [[qT; qT] for t = 1:T],
    [zeros(model.m) for t = 1:T])

obj = MultiObjective([obj_tracking, obj_penalty])

# Constraints
include_constraints("contact")
con_contact = contact_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_contact)

# Trajectory initialization
q_ref = [q1, linear_interpolation(q1, qT, T)...]
x0 = configuration_to_state(q_ref) # linear interpolation on state
u0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time z̄, info = solve(prob, copy(z0),
    tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

check_slack(z̄, prob)

x̄, ū = unpack(z̄, prob)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model,
    # q_ref,
    state_to_configuration([[x̄[1] for i = 1:10]...,x̄..., [x̄[end] for i = 1:10]...]),
    Δt = h)
# open(vis)

# using Plots
plot(hcat(ū...)[model.idx_u, :]', linetype = :steppost, )

model

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Gridap
using Test

include("CellValuesBench.jl")
include("CellFunctionsBench.jl")
include("GeometryBenchs.jl")


text:
# Utility functions

"""
    all_values(x)

Return a tuple of all the values contained in an instance of a composite
type `x`.

# Example
```
julia> struct A
       x
       y
       end

julia> all_values(A(1, "2"))
(1, "2")
```
"""
@generated function all_values(x)
    quote
        ($([:(x.$f) for f in fieldnames(x)]...),)
    end
end

"""
    transform_name(s::AbstractString) -> s′::Symbol

Transform the name of an attribute or element of an FDSN Station XML
document into a `Symbol` suitable for assignment into a `struct`.

For example:

```julia
julia> name = "EmailType";

julia> transform_name(name)
:email
```
"""
function transform_name(s)
    # Special cases
    if s == "FrequencyDBVariation"
        return :frequency_db_variation
    elseif s == "end"
        return :end_
    end
    s = string(s)
    # CamelCase to Camel_Case
    s = replace(s, r"([a-z])([A-Z])"=>s"\1_\2")
    # lowercase
    s = lowercase(s)
    # Special cases
    s = replace(s, r"^module$"=>"module_name")
    Symbol(s)
end

"""
    xml_element_name(s::Symbol) -> name

Transform a field name `s` of a struct into its equivalent FDSN StationXML
element name.

# Example
```julia
julia> StationXML.xml_element_name(:module_name)
"Module"
```
"""
function xml_element_name(s::Symbol)
    # Special cases
    if s === :module_name
        return "Module"
    elseif s === :module_uri
        return "ModuleURI"
    elseif s === :fir
        return "FIR"
    elseif s === :frequency_db_variation
        return "FrequencyDBVariation"
    elseif s === :uri
        return "URI"
    end
    name = String(s)
    name = replace(name, "_"=>" ")
    name = titlecase(name)
    name = replace(name, " "=>"")
end

"""
    xml_attribute_name(s::Symbol) -> name

Transform a field name `s` of a struct into its equivalent FDSN StationXML
attribute name.

# Example
```julia
julia> StationXML.xml_attribute_name(:schema_version)
"schemaVersion"
```
"""
function xml_attribute_name(s::Symbol)
    if s === :schema_version
        return "schemaVersion"
    elseif s === :unit
        return "unit"
    elseif s === :location_code
        return "locationCode"
    elseif s === :resource_id
        return "resourceId"
    elseif s === :number
        return "number"
    elseif s === :id
        return "id"
    elseif s === :i
        return "i"
    elseif s === :number
        return "number"
    elseif s === :plus_error
        return "plusError"
    elseif s === :minus_error
        return "minusError"
    elseif s === :measurement_method
        return "measurementMethod"
    elseif s === :datum
        return "datum"
    elseif s === :description
        return "description"
    elseif s === :name
        return "name"
    elseif s === :code
        return "code"
    elseif s === :start_date
        return "startDate"
    elseif s === :end_date
        return "endDate"
    elseif s === :source_id
        return "sourceID"
    elseif s === :restricted_status
        return "restrictedStatus"
    elseif s === :alternate_code
        return "alternateCode"
    elseif s === :historical_code
        return "historicalCode"
    elseif s === :subject
        return "subject"
    elseif s === :start
        return "start"
    elseif s === :end_
        return "end"
    elseif s === :maximum_time_tear
        return "maximumTimeTear"
    elseif s === :number_segments
        return "numberSegments"
    elseif s === :type
        return "type"
    else
        throw(ArgumentError("$s is not a known StationXML attribute name"))
    end
end

"""
    xml_unescape(s) -> s′

Replace escaped occurrences of the five XML character entity
references &, <, >, " and ' with their unescaped equivalents.

Reference:
    https://en.wikipedia.org/wiki/Character_encodings_in_HTML#XML_character_references
"""
xml_unescape(s) =
    # Workaround for JuliaLang/julia#28967
                  reduce(replace,
                         ("&"  => "&",
                          "<"   => "<",
                          ">"   => ">",
                          """ => "\"",
                          "'" => "'"),
                         init=s)

"""
    xml_escape(s) -> s′

Replace unescaped occurrences of the five XML characters
&, <, >, " and ' with their escaped equivalents.

Reference:
    https://en.wikipedia.org/wiki/Character_encodings_in_HTML#XML_character_references
"""
xml_escape(s) =
    # Workaround for JuliaLang/julia#28967
                reduce(replace,
                       ("&"  => "&",
                        "<"  => "<",
                        ">"  => ">",
                        "\"" => """,
                        "'"  => "'"),
                       init=s)

"""
    xs_nmtoken_or_throw(s) -> token

Return `token`, which is guaranteed to conform to the `xsd:NMTOKEN` specification
of an XML token.  If `s` cannot be made into a valid token, then an
`ArgumentError` is thrown.

From the XML Schema 1.0 specification:
> The type `xsd:NMTOKEN` represents a single string token. `xsd:NMTOKEN`
> values may consist of letters, digits, periods (`.`), hyphens (`-`),
> underscores (`_`), and colons (`:`). They may start with any of these
> characters. `xsd:NMTOKEN` has a `whiteSpace` facet value of `collapse`,
> so any leading or trailing whitespace will be removed. However, no
> whitespace may appear within the value itself.

!!! note
    This implementation forbids non-ASCII characters, though the specification
    [may tolerate them](http://books.xmlschemata.org/relaxng/ch19-77231.html).
"""
function xs_nmtoken_or_throw(s)
    s = strip(s)
    occursin(r"^[a-zA-Z0-9._\-:]*$", s) ||
        throw(ArgumentError("string \"$s\" is not a valid token"))
    s
end

"""
    abstract type EnumeratedStruct end

Supertype of all structs made with the [`@enumerated_struct`](@ref) macro.
"""
abstract type EnumeratedStruct end

"""
    @enumerated_struct(name, values)

Create a `struct` called `name` with a single `String` field, `:value`, which
must match one of the items in `values`.

Create a keyword constructor with one required keyword argument,
`value`, and add a method for `permitted_values` returning the values
which can be used in `name`.

The struct is a subtype of [`StationXML.EnumeratedStruct`](@ref),
and hence has methods defined for the following functions:

- [`StationXML.attribute_fields`](@ref)
- [`StationXML.element_fields`](@ref)
- [`StationXML.has_text_field`](@ref)
- [`StationXML.text_field`](@ref)
- [`StationXML.parse_node`](@ref)
- [`StationXML.local_parse`](@ref)

See [`StationXML.EnumeratedStruct`](@ref) for more details on other
methods defined for the new type.

# Example

The following code
```
@enumerated_struct(Example, ("A", "B", "D"))
```

is equivalent to
```
struct Example <: StationXML.EnumeratedStruct
    value::String
    function Example(value)
        if value ∉ ("A", "B", "D")
            throw(ArgumentError("value must be one of \"A\", \"B\" or \"D\""))
        end
        new(value)
    end
    Example(; value) = Example(value)
end

# WARNING!
# This method is added to the module from which the macro is called,
# not StationXML!
permitted_values(::Type{Example}) = (\"A\", \"B\", \"D\")
```
"""
macro enumerated_struct(name, values)
    hasfield(typeof(values), :head) && values.head === :tuple ||
        throw(ArgumentError("final argument must be a tuple of values"))
    values_string = string(values)
    # A nice set of Markdown-formatted values.
    values_docstring = join((string("`\"", val, "\"`") for val in values.args), ", ", " or ")
    # For the example in the docstring
    n_values = length(values.args)
    n_values == 0 && throw(ArgumentError("number of permitted values cannot be 0"))
    first_value = string(first(values.args))
    second_value = n_values == 1 ? first_value : string(values.args[2])
    module_name = :($(@__MODULE__).$name)
    name = esc(name)
    quote
        """
            $($name) <: EnumeratedStruct

        # Constructors

            $($name)(value)
            $($name)(; value)

        Enumerated struct containing a single string which must be one
        of the following: $($values_docstring).

        This field is named `value`.

        Note that when a field of another type is a `$($name)`, it
        is not necessary to assign a field of type `$($name)` to the
        field.  Instead, one can simply use a `String`, from which a
        `$($name)` will be automatically constructed.

        For this reason, `$($name)` is not exported even when bringing
        StationXML's types into scope by doing `using StationXML.Types`.

        # Example
        ```
        julia> using StationXML

        julia> mutable struct ExampleStruct
                   field::$($name)
               end

        julia> es = ExampleStruct("$($first_value)")
        ExampleStruct($($name)("$($first_value)"))

        julia> es.field = "$($second_value)"
        "$($second_value)"
        ```
        """
        struct $name <: $(esc(EnumeratedStruct))
            value::String
            function $name(value)
                if value ∉ $values
                    $(esc(throw))($(esc(ArgumentError))("value must be one of $($values_string)"))
                end
                new(value)
            end
            $name(; value) = $name(value)
        end

        # FIXME: Work out correct escaping so permitted_values is always
        #        added to the StationXML module and not the module in which
        #        the macro is called.
        #        In other words, this currently adds a method to the
        #        module where we call @enumerated_struct, though this
        #        is fine if it's only used in this module.
        $(esc(:(permitted_values)))(::$(esc(Type)){$name}) = $values
    end
end

attribute_fields(::Type{<:EnumeratedStruct}) = ()
element_fields(::Type{<:EnumeratedStruct}) = ()
has_text_field(::Type{<:EnumeratedStruct}) = true
text_field(::Type{<:EnumeratedStruct}) = :value
Base.convert(::Type{S}, s::EnumeratedStruct) where {S<:AbstractString} = S(s.value)
Base.convert(::Type{T}, s::AbstractString) where {T<:EnumeratedStruct} = T(s)
parse_node(::Type{T}, node::EzXML.Node, warn::Bool=false) where {T<:EnumeratedStruct} = T(node.content)
local_parse(::Type{T}, s::AbstractString) where {T<:EnumeratedStruct} = T(s)
Base.parse(::Type{T}, s::AbstractString) where {T<:EnumeratedStruct} = T(s)


text:
using Test
using XLSX: readxlsx
using DataFrames
using Mimi
using MimiDICE2016R2
using CSVFiles

using MimiDICE2016R2: getparams

@testset "MimiDICE2016R2" begin

    #------------------------------------------------------------------------------
    #   1. Run tests on the whole model
    #------------------------------------------------------------------------------

    @testset "MimiDICE2016R2-model" begin

        df = load(joinpath(@__DIR__, "../data/DICE2016-R2-GAMSoutputs.csv")) |> DataFrame

        m = MimiDICE2016R2.get_model();
        run(m)

        # TODO: validate this model, noting we will need to update certain parameters
        # and rerun the Mimi version so that we are imitating the optimized pathway
        # in GAMS

        #Test Precision
        Precision = 1.0e-10

        #Time Periods
        T=100

        # #TATM Test (temperature increase)
        # @test maximum(abs, m[:climatedynamics, :TATM] .- True_TATM) ≈ 0. atol = Precision

        # #MAT Test (carbon concentration atmosphere)
        # @test maximum(abs, m[:co2cycle, :MAT] .- True_MAT) ≈ 0. atol = Precision

        # #DAMFRAC Test (damages fraction)
        # @test maximum(abs, m[:damages, :DAMFRAC] .- True_DAMFRAC) ≈ 0. atol = Precision

        # #DAMAGES Test (damages $)
        # @test maximum(abs, m[:damages, :DAMAGES] .- True_DAMAGES) ≈ 0. atol = Precision

        # #E Test (emissions)
        # @test maximum(abs, m[:emissions, :E] .- True_E) ≈ 0. atol = Precision

        # #YGROSS Test (gross output)
        # @test maximum(abs, m[:grosseconomy, :YGROSS] .- True_YGROSS) ≈ 0. atol = Precision

        # #AL test (total factor productivity)
        # @test maximum(abs, m[:totalfactorproductivity, :AL] .- True_AL) ≈ 0. atol = Precision

        # #CPC Test (per capita consumption)
        # @test maximum(abs, m[:neteconomy, :CPC] .- True_CPC) ≈ 0. atol = Precision

        # #FORCOTH Test (exogenous forcing)
        # @test maximum(abs, m[:radiativeforcing, :FORCOTH] .- True_FORCOTH) ≈ 0. atol = Precision

        # #FORC Test (radiative forcing)
        # @test maximum(abs, m[:radiativeforcing, :FORC] .- True_FORC) ≈ 0. atol = Precision

        # #Utility Test
        # @test maximum(abs, m[:welfare, :UTILITY] .- True_UTILITY) ≈ 0. atol = Precision

    end #MimiDICE2016R2-model testset

    #------------------------------------------------------------------------------
    #   2. Run tests on SCC
    #------------------------------------------------------------------------------

    @testset "Standard API" begin

        m = MimiDICE2016R2.get_model()
        run(m)

        # Test the errors
        @test_throws ErrorException MimiDICE2016R2.compute_scc()  # test that it errors if you don't specify a year
        @test_throws ErrorException MimiDICE2016R2.compute_scc(year=2021)  # test that it errors if the year isn't in the time index
        @test_throws ErrorException MimiDICE2016R2.compute_scc(last_year=2299)  # test that it errors if the last_year isn't in the time index
        @test_throws ErrorException MimiDICE2016R2.compute_scc(year=2105, last_year=2100)  # test that it errors if the year is after last_year

        # Test the SCC
        scc1 = MimiDICE2016R2.compute_scc(year=2020)
        @test scc1 isa Float64

        # Test that it's smaller with a shorter horizon
        scc2 = MimiDICE2016R2.compute_scc(year=2020, last_year=2200)
        @test scc2 < scc1

        # Test that it's smaller with a larger prtp
        scc3 = MimiDICE2016R2.compute_scc(year=2020, last_year=2200, prtp=0.02)
        @test scc3 < scc2

        # Test with a modified model
        m = MimiDICE2016R2.get_model()
        update_param!(m, :climatedynamics, :t2xco2, 5)
        scc4 = MimiDICE2016R2.compute_scc(m, year=2020)
        @test scc4 > scc1   # Test that a higher value of climate sensitivty makes the SCC bigger

        # Test compute_scc_mm
        result = MimiDICE2016R2.compute_scc_mm(year=2030)
        @test result.scc isa Float64
        @test result.mm isa Mimi.MarginalModel
        marginal_temp = result.mm[:climatedynamics, :TATM]
        @test all(marginal_temp[1:findfirst(isequal(2030), MimiDICE2016R2.model_years)] .== 0.)
        @test all(marginal_temp[findfirst(isequal(2035), MimiDICE2016R2.model_years):end] .!= 0.)

    end

    # ------------------------------------------------------------------------------
    #   3. Deterministic SCC values
    # ------------------------------------------------------------------------------
    @testset "SCC values" begin

        atol = 1e-6 # TODO what is a reasonable tolerance given we test on a few different machines etc.

        # Test several validation configurations against the pre-saved values from previous version MimiDICE2016R2
        specs = Dict([
            :year => [2020],
            :eta => [0, 1.5],
            :prtp => [0.015, 0.03],
            :last_year => [2200, 2300],
        ])

        results = DataFrame(year = [], eta = [], prtp = [], last_year = [], SC = [])

        for year in specs[:year]
            for eta in specs[:eta]
                for prtp in specs[:prtp]
                    for last_year in specs[:last_year]
                        sc = MimiDICE2016R2.compute_scc(year=Int(year), eta=eta, prtp=prtp, last_year=Int(last_year))
                        push!(results, (year, eta, prtp, last_year, sc))
                    end
                end
            end
        end

        validation_results = load(joinpath(@__DIR__, "..", "data", "SC validation data", "deterministic_sc_values_v0-1-0.csv")) |> DataFrame
        # println("MAXIMUM DIFF IS $(maximum(results[!, :SC] - validation_results[!, :SC]))")
        @test all(isapprox.(results[!, :SC], validation_results[!, :SC], atol = atol))

    end # SCC values testset

end #MimiDICE2016R2 testset

nothing


text:
using Metrics
using Test
using Random

Random.seed!(0)

@testset "Metrics.jl" begin
    include("regression.jl")
    include("classification.jl")
    include("rank.jl")
    include("nlp.jl")
    include("cv.jl")
end


text:
using PyPlot
export noteplotter
maxnotevel(::Note) = 127
maxnotevel(::MoreVelNote) = 157

function plotpianonote!(ax, note, cmap)
    r = PyPlot.matplotlib.patches.Rectangle(
        (note.position, note.pitch - 0.5), note.duration, 1,
        color = cmap(note.velocity/maxnotevel(note))
    )
    ax.add_artist(r)
    return note.pitch
end

"""
    noteplotter(notes::Notes; kwargs...)
Plot the given `notes` on a "piano roll" like setting with
- x axis being the position of the notes
- y axis being the "pitch" or "value" of the notes (see below)
- color being the velocity of the notes

Then return the values of the y axis.

## Keywords
* `st = (notes[1].position ÷ notes.tpq) * notes.tpq`  time to start plotting from
* `fi = st + 16notes.tpq` time to stop plotting at, by default 16 quarter notes, i.e. four bars.
  Give `Inf` if you want to plot until the end of the notes.
* `ax = (PyPlot.figure(); PyPlot.gca())` the axis to plot on.
* `cmap = "viridis"` the colormap to use for the velocity.
* `grid = 0:0.25:1` a grid to plot along with the notes (by default the 16th notes).
  Give nothing if you don't want grid lines to be plotted.
* `names` :
  a dictionary that given the y-axis value returns how this value should be named.
  By default the function [`pitch_to_name`](@ref) is used, along with a
  heuristic to only name around 7 notes.
* `plotnote!` A function with call arguments
  `plotnote!(ax, note, cmap)` (with `cmap` a colormap instance, not a string),
  that actually plots the notes. By default plots a "piano-roll".

The `plotnote!` argument allows for customization. The function is supposed
to plot a note on the given axis **and return** the "value" of the note.
See the official documentation for an example of how this is taken
advantage of to e.g. plot drum notes.
"""
function noteplotter(notes::Notes;
    st = (notes[1].position ÷ notes.tpq) * notes.tpq,
    fi = st + 16notes.tpq,
    ax = (PyPlot.figure(); PyPlot.gca()),
    names = nothing,
    plotnote! = plotpianonote!,
    grid = 0:0.25:1,
    cmap = "viridis"
    )

    cm = PyPlot.cm.get_cmap(cmap)
    fi == Inf && (fi = notes[end].position + notes[end].duration)

    # plot all notes:
    plottedpitches = Int[]
    for note in notes
        st > note.position && continue
        fi < note.position && break
        p = plotnote!(ax, note, cm)
        p ∉ plottedpitches && push!(plottedpitches, p)
    end

    # add gridpoints
    if !isnothing(grid) && isgrid(grid)
        gpoints = notes.tpq*grid[2:end-1]
        grid = UInt[]
        beat = UInt[]
        for i = (st÷notes.tpq):(fi÷notes.tpq + 1)
            append!(grid, i*notes.tpq .+ gpoints)
            push!(beat, i*notes.tpq)
        end
        for γ in grid
            ax.axvline(γ, color=fill(0.7, 3), lw = 1, ls="dashed", zorder = -5)
        end
        for b in beat
            ax.axvline(b, lw=1.5, color=fill(0.5, 3), zorder = -4)
        end
    end

    # limits and labels
    if names == nothing
        yticks, ylabels = heuristic_note_names(plottedpitches)
    else
        yticks = plottedpitches
        ylabels = [names[p] for p in plottedpitches]
    end
    ax.set_xlim(st, fi)
    ax.set_yticks(yticks)
    ax.set_ylim(minimum(plottedpitches)-1, maximum(plottedpitches)+1)
    ax.set_yticklabels(ylabels)
    ax.set_xlabel("time (ticks)")
    return plottedpitches
end

function heuristic_note_names(plottedpitches)
    mi, ma = extrema(plottedpitches)
    notestep = max((ma-mi)÷7, 1)
    pitches = mi:notestep:ma
    ylabels = [pitch_to_name(p) for p in pitches]
    return pitches, ylabels
end


text:
mutable struct SparseSDP{T<:Number}
    obj::SparseSymmetricBlockMatrix{T}
    cons::Dict{Any,SparseSymmetricBlockMatrix{T}}
    rhs::Dict{Any,T}
    maximize::Bool
    normalized::Bool
    freecons::Dict{Any}{Dict{Any}{T}}
    freeobj::Dict{Any}{T}
end

SparseSDP(T::Type; maximize=true, normalized=false) = SparseSDP(SparseSymmetricBlockMatrix(T),
                                                                Dict{Any,SparseSymmetricBlockMatrix{T}}(),
                                                                Dict{Any,T}(),
                                                                maximize,
                                                                normalized,
                                                                Dict{Any}{Dict{Any}{T}}(),
                                                                Dict{Any}{T}())

SparseSDP(; maximize=true, normalized=false) = SparseSDP(Float64, maximize=maximize, normalized=normalized)

ismaximizationproblem(sdp::SparseSDP) = sdp.maximize

setmaximizationproblem!(sdp::SparseSDP, maximize::Bool) = sdp.maximize = maximize

normalized(sdp::SparseSDP) = sdp.normalized

copy(sdp::SparseSDP) = SparseSDP(copy(obj(sdp)), copy(cons(sdp)), copy(rhs(sdp)), sdp.maximize)

# There are (4) ways to set values to the SDP each depending on the input given.
setobj!(sdp::SparseSDP{Float64}, c::SparseSymmetricBlockMatrix{Float64}) = obj(sdp) = obj

function setobj!(sdp::SparseSDP{Float64}, bi, m::SparseSymmetricMatrix{Float64})
    obj(sdp)[bi] = m
end

function setobj!(sdp::SparseSDP{Float64}, bi, i, j, v::Float64)
    obj(sdp)[bi, i, j] = v
end

function setobj!(sdp::SparseSDP{Float64}, bi, m::Matrix{Float64})
    n = size(m)[1]
    for j = 1:n
        for i = 1:j
            if m[i, j] != zero(T)
                setobj!(sdp, bi, i, j, m[i, j])
            end
        end
    end
end

obj(sdp::SparseSDP) = sdp.obj

obj(sdp::SparseSDP, bi) = sdp.obj[bi]

freecons(sdp::SparseSDP) = sdp.freecons

freeobj(sdp::SparseSDP) = sdp.freeobj

function setcon!(sdp::SparseSDP{Float64}, ri, bi, i, j, v::Float64)
    if !isapprox(v, 0.0)
        if haskey(cons(sdp), ri)
            cons(sdp)[ri][bi, i, j] = v
        else
            a = SparseSymmetricBlockMatrix(Float64)
            a[bi, i, j] = v
            cons(sdp)[ri] = a
        end
    end
end

function setcon!(sdp::SparseSDP{Float64}, ri, bi, m::AbstractMatrix{Float64})
    for i = 1:size(m, 1)
        for j = 1:size(m, 2)
            if m[i, j] != zero(Float64)
                setcon!(sdp, ri, bi, i, j, m[i, j])
            end
        end
    end
end

function setfreecon!(sdp::SparseSDP{Float64}, ri, ci, v::Float64)
    if !isapprox(v, 0.0)
        if haskey(D(sdp), ri)
            freecon(sdp)[ri][ci] = Float64
        else
            freecon(sdp)[ri] = Dict{Any,Float64}(ci => v)
        end
    end
end

setfreeobj!(sdp::SparseSDP{Float64}, ri, v::Float64) = freeobj(sdp)[ri] = v

cons(sdp::SparseSDP) = sdp.cons

getcon(sdp::SparseSDP, ri) = sdp.cons[ri]

setrhs!(sdp::SparseSDP{Float64}, ri, v::Float64) = sdp.rhs[ri] = v

rhs(sdp::SparseSDP) = sdp.rhs

rhs(sdp::SparseSDP, ri) = sdp.rhs[ri]

function rhsdense(sdp::SparseSDP{Float64})
    l = zeros(Float64, ncons(sdp))
    for (ri, v) in rhs(sdp)
        l[ri] = v
    end
    l
end

function ncons(sdp::SparseSDP)
    v = 0
    for ri in keys(cons(sdp))
        if ri > v
            v = ri
        end
    end
    for ri in keys(rhs(sdp))
        if ri > v
            v = ri
        end
    end
    v
end

function blocksizes(sdp::SparseSDP)
    blockdict = Dict{Any,Set}()
    for (bi, block) in obj(sdp).blocks
        indices = get(blockdict, bi, Set())
        union!(indices, block.indices)
        blockdict[bi] = indices
    end
    for a in values(cons(sdp))
        for (bi, block) in a
            indices = get(blockdict, bi, Set())
            union!(indices, block.indices)
            blockdict[bi] = indices
        end
    end
    blocksizedict = Dict{Any,Int}()
    for (bi, entries) in blockdict
        blocksizedict[bi] = length(entries)
    end
    blocksizedict
end

max_stars_repo_name: JAndriesJ/SemidefiniteProgramming

max_issues_repo_name: JAndriesJ/SemidefiniteProgramming

max_forks_repo_name: JAndriesJ/SemidefiniteProgramming


text:
using DataFrames, CSV
include("../src/sdistance.jl") # 必ずダブルクオーテーション
include("../test/draw.jl")
include("../src/SignedDistanceFunction.jl")
import .Sdistance: benchmark_floodfill, benchmark_singlecurves_isinside #,signedDistance2D
using .SignedDistanceFunction
import .Draw: parformance_graphs

@enum ExecuteKinds _multicurves _singlecurve _singlecurve_floodfill

function _exe(kinds)
    _execute_times = 2
    _phi = []
    init_N = 100
    increment_N = 200
    runtime = zeros(_execute_times + 1, 2)
    N = [init_N + increment_N * item for item = 0:_execute_times]
    if kinds == 1
        for i = 0:_execute_times
            runtime[i+1, 1] = benchmark_floodfill(init_N + increment_N * i, "./test/mock_csv_data/multiple_curves.csv")
            runtime[i+1, 2] = benchmark_floodfill(init_N + increment_N * i, "./test/mock_csv_data/multiple_curves.csv", false)
        end
        parformance_graphs(N, runtime, "multiple_curves", ["Parallel processing", "Normal processing"])
    elseif kinds == 2
        for i = 0:_execute_times
            runtime[i+1, 1] = benchmark_floodfill(init_N + increment_N * i, "./test/mock_csv_data/interface.csv")
            runtime[i+1, 2] = benchmark_floodfill(init_N + increment_N * i, "./test/mock_csv_data/interface.csv", false)
        end
        parformance_graphs(N, runtime, "interface_floodfill", ["Parallel processing", "Normal processing"])
    elseif kinds == 3
        for i = 0:_execute_times
            runtime[i+1, 1] = benchmark_singlecurves_isinside(init_N + increment_N * i, "./test/mock_csv_data/interface.csv")
            runtime[i+1, 2] = benchmark_singlecurves_isinside(init_N + increment_N * i, "./test/mock_csv_data/interface.csv", false)
        end
        parformance_graphs(N, runtime, "the jordan curve", ["Parallel processing", "Normal processing"])
    end
end

# _exe(1)
# _exe(2)
_exe(3)

#==
Jabelic test it and get it as below.

N2300:  300+200i, ARG = 10
 [300, 500, 700, 900, 1100, 1300, 1500, 1700, 1900, 2100, 2300]
 [2.727736517 1.549037993; 2.430711433 4.540753791;
 4.87217261 8.72733049; 6.970006918 11.296165896;
 10.919117877 16.429877231; 15.524662846 28.576650596;
 20.630861499 38.540759007; 26.182433943 51.252256581;
 33.966169503 64.128344396; 40.489381557 75.419365612;
 50.053971199 91.594363467]

300+200i, ARG = 3, interface  multi-thread, single-thread
 [2.402844242 1.564786665; 2.415855251 4.372297246; 4.636160664 9.164722692; 7.243630318 11.635953353]
signedDistance2D, 300+200i, ARG = 3, interface, multi-thread only. jordan_curve(using isinside), multi_curves(using floodfill)
 [2.126074484 13.096371426; 2.496373027 125.23088408; 5.242287246 596.261015286; 7.987274914 1743.300632586]

1-3-500, interface, inpolygon
[100, 300, 500]
[1.248799153 0.530269039; 2.28104107 4.735589021; 6.293578857 11.423713058]
32

1-3-500, interface, Floodfill
[100, 300, 500][1.291998074 0.123935482; 3.388728651 5.878153884; 69.895109677 84.980423578]
32

1-3-500, multi_curves, Floodfill
[100, 300, 500][2.400877896 2.741995989; 38.038218511 92.339324457; 476.64809553 833.91326708]
32

1-3-500, multi_curves, Floodfill
[100, 300, 500][0.5009443733333333 0.891537955; 12.114947880666668 28.633708910333336; 133.44351538533337 261.82990826500003]

 ==#

max_stars_repo_name: jabelic/Signed-distance-function

max_issues_repo_name: jabelic/Signed-distance-function

max_forks_repo_name: jabelic/Signed-distance-function


text:
@testset "populate_lp_model!" begin
    form = _build_dmip_formulation()
    state = Cerberus.CurrentState()
    let node = Cerberus.Node()
        node_result = Cerberus.NodeResult(node)
        @inferred Cerberus.populate_lp_model!(
            state,
            form,
            node,
            node_result,
            CONFIG,
        )
        model = state.gurobi_model

        @test MOI.get(model, MOI.NumberOfVariables()) == 3
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1
        c1 = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_ET}())[1]
        f1 = MOI.get(model, MOI.ConstraintFunction(), c1)
        _is_equal(f1, 1.0 * _SV(_VI(1)) + 2.1 * _SV(_VI(2)) + 3.0 * _SV(_VI(3)))
        s1 = MOI.get(model, MOI.ConstraintSet(), c1)
        @test s1 == _ET(3.0)

        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 1
        c2 = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_LT}())[1]
        f2 = MOI.get(model, MOI.ConstraintFunction(), c2)
        _is_equal(f2, -3.5 * _SV(_VI(1)) + 1.2 * _SV(_VI(2)))
        s2 = MOI.get(model, MOI.ConstraintSet(), c2)
        @test s2 == _LT(4.0)
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 0

        @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(1), state),
        ) == (0.5, 1.0)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(2), state),
        ) == (-1.3, 2.3)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(3), state),
        ) == (0.0, 1.0)

        @test MOI.get(model, MOI.ObjectiveSense()) == MOI.MIN_SENSE
        obj = MOI.get(model, MOI.ObjectiveFunction{_SAF}())
        _is_equal(obj, 1.0 * _SV(_VI(1)) - 1.0 * _SV(_VI(2)))
    end

    # Safe to do this now that we've populated the model once (undefined reference otherwise).
    model = state.gurobi_model
    # TODO: Test that model does NOT get created if the config/state says not to.
    state.rebuild_model = false
    let node = Cerberus.Node(
            [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
            [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))],
            2,
        )
        node_result = Cerberus.NodeResult(node)
        @inferred Cerberus.populate_lp_model!(
            state,
            form,
            node,
            node_result,
            CONFIG,
        )
        @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(1), state),
        ) == (0.5, 0.0)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(2), state),
        ) == (-1.3, 2.3)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(3), state),
        ) == (1.0, 1.0)
    end

    # Lie to Cerberus so that it doesn't reset the bounds
    state.on_a_dive = true
    let node = Cerberus.Node()
        node_result = Cerberus.NodeResult(node)
        @inferred Cerberus.populate_lp_model!(
            state,
            form,
            node,
            node_result,
            CONFIG,
        )
        @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(1), state),
        ) == (0.5, 0.0)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(2), state),
        ) == (-1.3, 2.3)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(3), state),
        ) == (1.0, 1.0)
    end

    state.on_a_dive = false
    let node = Cerberus.Node()
        node_result = Cerberus.NodeResult(node)
        @inferred Cerberus.populate_lp_model!(
            state,
            form,
            node,
            node_result,
            CONFIG,
        )
        @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(1), state),
        ) == (0.5, 1.0)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(2), state),
        ) == (-1.3, 2.3)
        @test MOI.Utilities.get_bounds(
            model,
            Float64,
            Cerberus.instantiate(_CVI(3), state),
        ) == (0.0, 1.0)
    end
end

@testset "update_node_bounds!" begin
    form = _build_dmip_formulation()
    state = _CurrentState()
    node = Cerberus.Node()
    node_result = Cerberus.NodeResult(node)
    @inferred Cerberus.populate_lp_model!(
        state,
        form,
        node,
        node_result,
        CONFIG,
    )
    model = state.gurobi_model
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 1.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)

    node = Cerberus.Node(
        [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
        [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))],
        Cerberus.AffineConstraint{_LT}[],
        Cerberus.AffineConstraint{_GT}[],
        2,
    )
    @inferred Cerberus.apply_branchings!(state, node)
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 0.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (1.0, 1.0)

    let f = _CSAF([1.2, 3.4], [_CVI(1), _CVI(2)], 5.6), s = _LT(7.8)
        bd = Cerberus.AffineConstraint(f, s)
        Cerberus.apply_branching!(node, bd)
        @inferred Cerberus.apply_branchings!(state, node)
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 2
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 0
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1
        cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_LT}())
        # NOTE: cis[1] should be the constraint from the "base" formulation.
        ci = cis[2]
        f_rt = MOI.get(model, MOI.ConstraintFunction(), ci)
        s_rt = MOI.get(model, MOI.ConstraintSet(), ci)
        # Constraint added with MOIU.normalize_and_add_constraint, which shifts
        # constant over to set.
        @test _is_equal(f_rt, Cerberus.instantiate(f, state) - 5.6)
        @test s_rt == _LT(7.8 - 5.6)
    end

    let f = _CSAF([2.4, 6.4], [_CVI(3), _CVI(1)], 0.0), s = _GT(3.5)
        bd = Cerberus.AffineConstraint(f, s)
        Cerberus.apply_branching!(node, bd)
        @inferred Cerberus.apply_branchings!(state, node)
        # NOTE: We are testing here that the _LT general branching constraint
        # only gets added to the model once.
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 2
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 1
        @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1
        cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_GT}())
        ci = cis[1]
        f_rt = MOI.get(model, MOI.ConstraintFunction(), ci)
        s_rt = MOI.get(model, MOI.ConstraintSet(), ci)
        @test _is_equal(f_rt, Cerberus.instantiate(f, state))
        @test s_rt == s
    end
end

@testset "remove_all_branchings!" begin
    form = _build_dmip_formulation()
    state = Cerberus.CurrentState()
    node = Cerberus.Node()
    Cerberus.create_base_model!(state, form, node, CONFIG)
    f_lt = _CSAF([1.2, 3.4], [_CVI(1), _CVI(2)], 0.0)
    s_lt = _LT(7.8)
    f_gt = _CSAF([2.4, 6.4], [_CVI(3), _CVI(1)], 0.0)
    s_gt = _GT(3.5)
    node_1 = Cerberus.Node(
        [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
        [Cerberus.BoundUpdate(_CVI(3), _GT(1.0))],
        [Cerberus.AffineConstraint{_LT}(f_lt, s_lt)],
        [Cerberus.AffineConstraint{_GT}(f_gt, s_gt)],
        4,
    )
    Cerberus.apply_branchings!(state, node_1)

    model = state.gurobi_model
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 0.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (1.0, 1.0)
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 2
    let lt_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_LT}())
        lt_ci = lt_cis[2]
        f_lt_rt = MOI.get(model, MOI.ConstraintFunction(), lt_ci)
        s_lt_rt = MOI.get(model, MOI.ConstraintSet(), lt_ci)
        @test _is_equal(f_lt_rt, Cerberus.instantiate(f_lt, state))
        @test s_lt_rt == s_lt
    end
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 1
    let gt_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_GT}())
        gt_ci = gt_cis[1]
        f_gt_rt = MOI.get(model, MOI.ConstraintFunction(), gt_ci)
        s_gt_rt = MOI.get(model, MOI.ConstraintSet(), gt_ci)
        @test _is_equal(f_gt_rt, Cerberus.instantiate(f_gt, state))
        @test s_gt_rt == s_gt
    end
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1

    @inferred Cerberus.remove_all_branchings!(state, form)
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 1.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 0
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1

    node_2 = Cerberus.Node(
        [Cerberus.BoundUpdate(_CVI(1), _LT(0.0))],
        Cerberus.BoundUpdate{_GT}[],
        Cerberus.AffineConstraint{_LT}[],
        [Cerberus.AffineConstraint{_GT}(f_gt, s_gt)],
        2,
    )
    Cerberus.apply_branchings!(state, node_2)
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 0.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 1
    let gt_cis = MOI.get(model, MOI.ListOfConstraintIndices{_SAF,_GT}())
        gt_ci = gt_cis[1]
        f_gt_rt = MOI.get(model, MOI.ConstraintFunction(), gt_ci)
        s_gt_rt = MOI.get(model, MOI.ConstraintSet(), gt_ci)
        @test _is_equal(f_gt_rt, Cerberus.instantiate(f_gt, state))
        @test s_gt_rt == s_gt
    end
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1

    @inferred Cerberus.remove_all_branchings!(state, form)
    @test MOI.get(model, MOI.NumberOfConstraints{_SV,_IN}()) == 3
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 1.0)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOI.Utilities.get_bounds(
        model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_LT}()) == 1
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_GT}()) == 0
    @test MOI.get(model, MOI.NumberOfConstraints{_SAF,_ET}()) == 1
end

@testset "apply_branchings!" begin
    form = _build_dmip_formulation()
    state = Cerberus.CurrentState()
    cs = state.constraint_state

    node = Cerberus.Node()
    Cerberus.create_base_model!(state, form, node, CONFIG)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    @inferred Cerberus.apply_branchings!(state, node)
    @test cs.branch_state.num_lt_branches == 0
    @test cs.branch_state.num_gt_branches == 0
    @test isempty(cs.branch_state.lt_general_constrs)
    @test isempty(cs.branch_state.gt_general_constrs)

    node_result = Cerberus.NodeResult(node)
    Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (0.5, 1.0)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)

    bd_1 = Cerberus.BoundUpdate(_CVI(1), _GT(1.0))
    Cerberus.apply_branching!(node, bd_1)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test Cerberus._unattached_bounds(cs, node, _GT) == [bd_1]
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    Cerberus.apply_branchings!(state, node)
    @test cs.branch_state.num_lt_branches == 0
    @test cs.branch_state.num_gt_branches == 1
    @test isempty(cs.branch_state.lt_general_constrs)
    @test isempty(cs.branch_state.gt_general_constrs)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (1.0, 1.0)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 1.0)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    bd_2 = Cerberus.BoundUpdate(_CVI(3), _LT(0.0))
    Cerberus.apply_branching!(node, bd_2)
    @test Cerberus._unattached_bounds(cs, node, _LT) == [bd_2]
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    Cerberus.apply_branchings!(state, node)
    @test cs.branch_state.num_lt_branches == 1
    @test cs.branch_state.num_gt_branches == 1
    @test isempty(cs.branch_state.lt_general_constrs)
    @test isempty(cs.branch_state.gt_general_constrs)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (1.0, 1.0)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 0.0)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    # Now we jump across the tree to a different node at the same depth.
    # Notably, we do not update the constraint_state. In particular, we still
    # believe that there has only been two variable branches. Therefore, the
    # model will still have the same bounds as above, and will NOT correspond
    # to the bounds in node below.
    node = Cerberus.Node()
    Cerberus.apply_branching!(node, Cerberus.BoundUpdate(_CVI(3), _GT(1.0)))
    Cerberus.apply_branching!(node, Cerberus.BoundUpdate(_CVI(1), _LT(0.0)))
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    Cerberus.apply_branchings!(state, node)
    @test cs.branch_state.num_lt_branches == 1
    @test cs.branch_state.num_gt_branches == 1
    @test isempty(cs.branch_state.lt_general_constrs)
    @test isempty(cs.branch_state.gt_general_constrs)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(1), state),
    ) == (1.0, 1.0)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(2), state),
    ) == (-1.3, 2.3)
    @test MOIU.get_bounds(
        state.gurobi_model,
        Float64,
        Cerberus.instantiate(_CVI(3), state),
    ) == (0.0, 0.0)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _LT))
    @test isempty(Cerberus._unattached_constraints(cs, node, _GT))

    lt_gen_1 = Cerberus.AffineConstraint(2.3 * _SV(_VI(1)), _LT(2.3))
    lt_gen_2 = Cerberus.AffineConstraint(2.5 * _SV(_VI(3)), _LT(2.5))
    gt_gen = Cerberus.AffineConstraint(2.4 * _SV(_VI(2)), _GT(2.4))
    Cerberus.apply_branching!(node, lt_gen_1)
    Cerberus.apply_branching!(node, lt_gen_2)
    Cerberus.apply_branching!(node, gt_gen)
    @test isempty(Cerberus._unattached_bounds(cs, node, _LT))
    @test isempty(Cerberus._unattached_bounds(cs, node, _GT))
    @test Cerberus._unattached_constraints(cs, node, _LT) ==
          [lt_gen_1, lt_gen_2]
    @test Cerberus._unattached_constraints(cs, node, _GT) == [gt_gen]

    Cerberus.apply_branchings!(state, node)
    @test cs.branch_state.num_lt_branches == 1
    @test cs.branch_state.num_gt_branches == 1
    @test length(cs.branch_state.lt_general_constrs) == 2
    @test length(cs.branch_state.gt_general_constrs) == 1
    @test MOI.get(state.gurobi_model, MOI.NumberOfConstraints{_SAF,_LT}()) == 3
    @test MOI.get(state.gurobi_model, MOI.NumberOfConstraints{_SAF,_GT}()) == 1
end

@testset "formulate_disjunctions!" begin
    # min  y + 1.2
    # s.t. (-2 <= x <= -1 & y = -x - 1) or
    #          (-1 <= x <= +1 & y = 0) or
    #          (+1 <= x <= +2 & y = x - 1)
    #      x + y >= -0.5
    #      x in [-1.5, 3.0]
    #      y in [0.0, 0.5]
    form = _build_formulation_with_single_disjunction(
        DisjunctiveConstraints.NaiveBigM(
            DisjunctiveConstraints.IntervalArithmetic(),
        ),
    )
    let node = Cerberus.Node()
        # min  y + 1.2
        # s.t. z_1 + z_2 + z_3 = 1
        #      x + y >= -0.5
        #      x >= -1.5 - 0.5 z_1
        #      x >= -1.5 + 0.5 z_2
        #      x >= -1.5 + 2.5 z_3
        #      y >= -x - 1.5 + 0.5 z_1
        #      y >= 0
        # .    y >= x - 3 + 2.0 z_3
        #      x <=  3.0 - 4.0 z_1
        #      x <=  3.0 - 2.0 z_2
        #      x <=  3.0 - 1.0 z_3
        #      y <= -x + 3.5 - 4.5 z_1
        #      y <= 0.5 - 0.5 z_2
        #      y <= x + 2.0 - 3.0 z_3
        #      x in [-1.5, 3.0]
        #      y in [0.0, 0.5]
        #      z in [0,1]^3
        state = Cerberus.CurrentState()
        node_result = Cerberus.NodeResult(node)
        Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
        Cerberus.apply_branchings!(state, node)
        x = [_SV(Cerberus.instantiate(_CVI(i), state)) for i in 1:5]
        @assert [v.variable.value for v in x] == collect(1:5)

        expected_bounds =
            [(-1.5, 3.0), (0.0, 0.5), (0.0, 1.0), (0.0, 1.0), (0.0, 1.0)]
        expected_et_acs = [(1.0 * x[3] + 1.0 * x[4] + 1.0 * x[5], _ET(1.0))]
        expected_gt_acs = [
            (1.0 * x[1] + 1.0 * x[2], _GT(-0.5)),
            (1.0 * x[1] + 0.5 * x[3], _GT(-1.5)),
            (1.0 * x[1] - 0.5 * x[4], _GT(-1.5)),
            (1.0 * x[1] - 2.5 * x[5], _GT(-1.5)),
            (1.0 * x[1] + 1.0 * x[2] - 0.5 * x[3], _GT(-1.5)),
            (1.0 * x[2], _GT(0.0)),
            (-1.0 * x[1] + 1.0 * x[2] - 2.0 * x[5], _GT(-3.0)),
        ]
        expected_lt_acs = [
            (1.0 * x[1] + 4.0 * x[3], _LT(+3.0)),
            (1.0 * x[1] + 2.0 * x[4], _LT(+3.0)),
            (1.0 * x[1] + 1.0 * x[5], _LT(+3.0)),
            (1.0 * x[1] + 1.0 * x[2] + 4.5 * x[3], _LT(+3.5)),
            (1.0 * x[2] + 0.5 * x[4], _LT(0.5)),
            (-1.0 * x[1] + 1.0 * x[2] + 3.0 * x[5], _LT(+2.0)),
        ]

        _test_roundtrip_model(
            state.gurobi_model,
            expected_bounds,
            expected_lt_acs,
            expected_gt_acs,
            expected_et_acs,
        )
    end

    let node = Cerberus.Node(
            [Cerberus.BoundUpdate(_CVI(3), _LT(0.0))],
            Cerberus.BoundUpdate{_GT}[],
            1,
        )
        # min  y + 1.2
        # s.t. z_1 + z_2 + z_3 = 1
        #      x + y >= -0.5
        #      x >= -1.5 + 0.5 z_2
        #      x >= -1.5 + 2.5 z_3
        #      y >= 0
        # .    y >= x - 3 + 2.0 z_3
        #      x <=  3.0 - 2.0 z_2
        #      x <=  3.0 - 1.0 z_3
        #      y <= 0.5 - 0.5 z_2
        #      y <= x + 2.0 - 3.0 z_3
        #      x in [-1.5, 3.0]
        #      y in [0.0, 0.5]
        #      z in [0,1]^3
        #      z_1 = 0
        state = Cerberus.CurrentState()
        node_result = Cerberus.NodeResult(node)
        Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
        Cerberus.apply_branchings!(state, node)
        x = [_SV(Cerberus.instantiate(_CVI(i), state)) for i in 1:5]
        @assert [v.variable.value for v in x] == collect(1:5)

        expected_bounds =
            [(-1.5, 3.0), (0.0, 0.5), (0.0, 0.0), (0.0, 1.0), (0.0, 1.0)]
        expected_et_acs = [(1.0 * x[3] + 1.0 * x[4] + 1.0 * x[5], _ET(1.0))]
        expected_gt_acs = [
            (1.0 * x[1] + 1.0 * x[2], _GT(-0.5)),
            (1.0 * x[1] - 0.5 * x[4], _GT(-1.5)),
            (1.0 * x[1] - 2.5 * x[5], _GT(-1.5)),
            (1.0 * x[2], _GT(0.0)),
            (-1.0 * x[1] + 1.0 * x[2] - 2.0 * x[5], _GT(-3.0)),
        ]
        expected_lt_acs = [
            (1.0 * x[1] + 2.0 * x[4], _LT(+3.0)),
            (1.0 * x[1] + 1.0 * x[5], _LT(+3.0)),
            (1.0 * x[2] + 0.5 * x[4], _LT(0.5)),
            (-1.0 * x[1] + 1.0 * x[2] + 3.0 * x[5], _LT(+2.0)),
        ]
        _test_roundtrip_model(
            state.gurobi_model,
            expected_bounds,
            expected_lt_acs,
            expected_gt_acs,
            expected_et_acs,
        )
    end
end

@testset "MOI.optimize!" begin
    form = _build_dmip_formulation()
    state = _CurrentState()
    node = Cerberus.Node()
    node_result = Cerberus.NodeResult(node)
    Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
    model = state.gurobi_model
    MOI.optimize!(model)
    @test MOI.get(model, MOI.TerminationStatus()) == MOI.OPTIMAL
    @test MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @test MOI.get(model, MOI.VariablePrimal(), _VI(1)) ≈ 0.5
    @test MOI.get(model, MOI.VariablePrimal(), _VI(2)) ≈ 2.5 / 2.1
    @test MOI.get(model, MOI.VariablePrimal(), _VI(3)) ≈ 0.0
end

@testset "_fill_solution!" begin
    form = _build_dmip_formulation()
    state = _CurrentState()
    node = Cerberus.Node()
    node_result = Cerberus.NodeResult(node)
    Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
    model = state.gurobi_model
    MOI.optimize!(model)
    @assert MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    @inferred Cerberus._update_lp_solution!(state, form)
    x = state.current_solution
    @test length(x) == 3
    @test x ≈ [1 / 2, 2.5 / 2.1, 0.0]
end

@testset "get_basis" begin
    form = _build_dmip_formulation()
    state = _CurrentState()
    node = Cerberus.Node()
    node_result = Cerberus.NodeResult(node)
    Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
    model = state.gurobi_model
    MOI.optimize!(model)
    @assert MOI.get(model, MOI.PrimalStatus()) == MOI.FEASIBLE_POINT
    true_basis = @inferred Cerberus.get_basis(state)
    expected_basis = Cerberus.Basis(
        [MOI.NONBASIC_AT_LOWER, MOI.BASIC, MOI.NONBASIC_AT_LOWER],
        [MOI.BASIC],
        MOI.BasisStatusCode[],
        [MOI.NONBASIC],
        MOI.BasisStatusCode[],
        MOI.BasisStatusCode[],
    )
    @test true_basis.base_var_constrs == expected_basis.base_var_constrs
    @test true_basis.base_lt_constrs == expected_basis.base_lt_constrs
    @test true_basis.base_gt_constrs == expected_basis.base_gt_constrs
    @test true_basis.base_et_constrs == expected_basis.base_et_constrs
    @test true_basis.branch_lt_constrs == expected_basis.branch_lt_constrs
    @test true_basis.branch_gt_constrs == expected_basis.branch_gt_constrs
end

function _set_basis_model(basis::Cerberus.Basis)
    form = _build_dmip_formulation()
    state = _CurrentState()
    node = Cerberus.Node(
        Cerberus.BoundUpdate{_LT}[],
        Cerberus.BoundUpdate{_GT}[],
        Cerberus.AffineConstraint{_LT}[],
        Cerberus.AffineConstraint{_GT}[],
        0,
        -Inf,
    )
    node_result = Cerberus.NodeResult(node)
    Cerberus.populate_lp_model!(state, form, node, node_result, CONFIG)
    state.warm_starts[node] = basis
    Cerberus.set_basis_if_available!(state, node)
    return state.gurobi_model
end

@testset "set_basis_if_available!" begin
    # First, seed a suboptimal basis. This will disable presolve. It is only one pivot away from the optimal basis.
    let
        subopt_basis = Cerberus.Basis(
            [MOI.BASIC, MOI.NONBASIC_AT_LOWER, MOI.NONBASIC_AT_LOWER],
            [MOI.BASIC],
            MOI.BasisStatusCode[],
            [MOI.NONBASIC],
            MOI.BasisStatusCode[],
            MOI.BasisStatusCode[],
        )
        model = _set_basis_model(subopt_basis)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 1
    end

    # Now, seed the optimal basis. This will solve the problem without any simplex iterations.
    let
        opt_basis = Cerberus.Basis(
            [MOI.NONBASIC_AT_LOWER, MOI.BASIC, MOI.NONBASIC_AT_LOWER],
            [MOI.BASIC],
            MOI.BasisStatusCode[],
            [MOI.NONBASIC],
            MOI.BasisStatusCode[],
            MOI.BasisStatusCode[],
        )
        model = _set_basis_model(opt_basis)
        MOI.optimize!(model)
        @test MOI.get(model, MOI.SimplexIterations()) == 0
    end
end


text:
using AbstractAlgebra

@testset "Groebner modular" begin
    R, (x, y) = PolynomialRing(QQ, ["x", "y"], ordering=:degrevlex)

    fs = [x, y]
    G = Groebner.groebner(fs)
    @test G ≂ [y, x]

    fs = [5y, 3//4*x]
    G = Groebner.groebner(fs)
    @test G ≂ [y, x]

    fs = [x^2 - (1//6)*y, x*y]
    G = Groebner.groebner(fs)
    @test G ≂ [y^2, x*y, x^2 - (1//6)y]

    fs = [QQ(11, 3)*x^2*y - QQ(2, 4)*x*y - QQ(1, 7)*y,
            QQ(1, 1)*x*y + QQ(7, 13)*x]
    G = Groebner.groebner(fs)
    @test G ≂ [y^2 + 7//13*y,
                x*y + 7//13*x,
                x^2 - 3//22*x + 39//539*y]

    #=
    root = Groebner.change_ordering(Groebner.rootn(3, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)

    root = Groebner.change_ordering(Groebner.rootn(4, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)

    root = Groebner.change_ordering(Groebner.rootn(8, ground=QQ), :degrevlex)
    gb = Groebner.groebner(root)
    @test Groebner.isgroebner(gb)

    noon = Groebner.change_ordering(Groebner.noon3(ground=QQ), :degrevlex)
    gb = Groebner.groebner(noon)
    @test Groebner.isgroebner(gb)

    noon = Groebner.change_ordering(Groebner.noonn(4, ground=QQ), :degrevlex)
    gb = Groebner.groebner(noon)
    @test Groebner.isgroebner(gb)
    =#
end

@testset "Groebner modular corner cases" begin
    R, (x, y, z, w) = PolynomialRing(QQ, ["x", "y", "z", "w"], ordering=:degrevlex)

    fs = [(12345678//12347)x,
          (222222221111123//2131232232097)y + z]
    G = Groebner.groebner(fs)
    @test G == [y + 2131232232097//222222221111123*z, x]

    fs = [
        (12345//12345678)x + y,
        x^2 + z^2,
        y^2 + w^2,
        x + 2*y + 3*z + 4*w
    ]
    G = Groebner.groebner(fs)

    @test G ≂ [w^3,
                z*w + 1692834523553//4063974*w^2,
                z^2 - 16935085031076//16933225*w^2,
                y - 12345//4106996*z - 4115//1026749*w,
                x + 6172839//2053498*z + 4115226//1026749*w]
    #@test Groebner.isgroebner(G)

    # what if we are unlucky to start from initial prime in coefficients
    G = [
        (2^30 + 3)x,
        (1//(2^30 + 3))y
    ]
    G = Groebner.groebner(G)
    @test G == [y, x]

    G = [
        (2^31 - 1)*x + y,
        (1//(2^31 - 1))*y + 1
    ]
    G = Groebner.groebner(G)
    @test G == [y + 2147483647, x - 1]

    R, (x, y, z, w) = PolynomialRing(QQ, ["x", "y", "z", "w"], ordering=:lex)
    fs = [(12345678//12347)x,
          (222222221111123//2131232232097)y + z]
    G = Groebner.groebner(fs)
    @test G == [y + 2131232232097//222222221111123*z, x]
end


text:
export iauLtpecl
"""
Long-term precession of the ecliptic.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   epj     double         Julian epoch (TT)

Returned:
   vec     double[3]      ecliptic pole unit vector

Notes:

   1. The returned vector is with respect to the J2000.0 mean equator
   and equinox.

   2. The Vondrak et al. (2011, 2012) 400 millennia precession model
   agrees with the IAU 2006 precession at J2000.0 and stays within
   100 microarcseconds during the 20th and 21st centuries.  It is
   accurate to a few arcseconds throughout the historical period,
   worsening to a few tenths of a degree at the end of the
   +/- 200,000 year time span.

References:

   Vondrak, J., Capitaine, N. and Wallace, P., 2011, New precession
   expressions, valid for long time intervals, Astron.Astrophys. 534,
   A22

   Vondrak, J., Capitaine, N. and Wallace, P., 2012, New precession
   expressions, valid for long time intervals (Corrigendum),
   Astron.Astrophys. 541, C1

This revision:  2016 February 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauLtpecl(double epj, double vec[3])
function iauLtpecl(epj::Real)
   vec = zeros(Float64, 3)
   ccall((:iauLtpecl, libsofa_c), Cvoid,
         (Cdouble, Ref{Cdouble}),
          convert(Float64, epj),
          vec)

   return SVector{3}(vec)
end


text:
# Autogenerated wrapper script for libfreenect_jll for x86_64-apple-darwin
export libfreenect, libfreenect_sync

using libusb_jll
JLLWrappers.@generate_wrapper_header("libfreenect")
JLLWrappers.@declare_library_product(libfreenect, "@rpath/libfreenect.0.dylib")
JLLWrappers.@declare_library_product(libfreenect_sync, "@rpath/libfreenect_sync.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(libusb_jll)
    JLLWrappers.@init_library_product(
        libfreenect,
        "lib/libfreenect.0.6.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libfreenect_sync,
        "lib/libfreenect_sync.0.6.1.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
gmres_test_pde(n)

PDE test from FR16. Test of kl_gmres with all kinds of preconditioning.
"""
function gmres_test_pde(n; orth = "cgs2", write = false, eta = 9.8 * 1.e-4)
    pdata = pdegminit(n)
    fdata = pdata.fdata
    RHS = pdata.RHS
    ue = pdata.ue
    b = Pfish2d(RHS, fdata)
    u0 = zeros(n * n)
    V = zeros(n * n,20)
    # Solve with left preconditioning hard-wired in
    goutp = kl_gmres(u0, b, pdelpatv, V, eta; pdata = pdata, orth = orth)
    pcres = goutp.reshist
    pcres /= pcres[1]
    sollhw = goutp.sol
    # Solve with right preconditioning hard-wired in
    goutrp = kl_gmres(u0, RHS, pderatv, V, eta; pdata = pdata, orth = orth)
    pcresr = goutrp.reshist
    pcresr /= pcresr[1]
    solrhw = pdeptv(goutp.sol, pdata)
    # Put left preconditioning in the argument list.
    goutpl2 =
        kl_gmres(u0, RHS, pdeatv, V, eta, pdeptv; pdata = pdata, orth = orth, side = "left")
    pcresl2 = goutpl2.reshist
    pcresl2 /= pcresl2[1]
    soll = goutpl2.sol
    # Put right preconditioning in the argument list.
    goutp2 = kl_gmres(
        u0,
        RHS,
        pdeatv,
        V,
        eta,
        pdeptv;
        pdata = pdata,
        orth = orth,
        side = "right",
    )
    pcres2 = goutp2.reshist
    pcres2 /= pcres2[1]
    solr = goutp2.sol
    soldel = norm(solrhw - solr, Inf)
    solrdel = norm(sollhw - soll, Inf)
    solerr = norm(soll - ue, Inf)
    solerr2 = norm(solr - ue, Inf)
    passfull = (
        (soldel == 0) &&
        (solrdel == 0) &&
        (solerr < 1.e-2) &&
        (solerr2 < 1.e-2) &&
        (length(pcresr) == 12) &&
        (length(pcres) == 9)
    )
    if write
        println(soldel, "  ", solrdel, "  ", solerr, "   ", solerr2)
    end
    passfull || println("Linear pde test for GMRES fails.")
# Now for some restarts ...
    V=zeros(n*n,5)
    goutpl = kl_gmres(u0, RHS, pdeatv, V, eta, pdeptv;
           pdata = pdata, orth = orth, side = "left", lmaxit=20)
    goutpr = kl_gmres(u0, RHS, pdeatv, V, eta, pdeptv;
           pdata = pdata, orth = orth, side = "right", lmaxit=20)
    soldelr=norm(goutpl.sol - goutpr.sol,Inf)
    solerrr=norm(goutpr.sol-ue)
    solerrl=norm(goutpl.sol-ue)
    numitsr=length(goutpr.reshist)
    numitsl=length(goutpl.reshist)
    pass_res= (
        (soldelr < 1.e-3) &&
        (solerr < 1.e-2) &&
        (solerr2 < 1.e-2) &&
        (numitsr == 16) &&
        (numitsl == 13)
    )
    pass_res || println("Linear pde test for GMRES(m) fails.")
    pass = passfull && pass_res
    return pass
end

"""
bicgstab_test_pde(n)

PDE test from FR16. Test of kl_bicgstab with all kinds of preconditioning.
"""
function bicgstab_test_pde(n; write = false, eta = 9.8 * 1.e-4)
    pdata = pdegminit(n)
    RHS = pdata.RHS
    ue = pdata.ue
    V = zeros(n*n)
    u0 = zeros(n * n)
#
    # Solve with left preconditioning hard-wired in
#
    fdata = pdata.fdata
    b = Pfish2d(RHS, fdata)
    goutp = kl_bicgstab(u0, b, pdelpatv, V, eta; pdata = pdata, lmaxit=200)
    pcres = goutp.reshist
    pcres /= pcres[1]
    sollhw = goutp.sol
#
    # Solve with right preconditioning hard-wired in
#
    goutrp = kl_bicgstab(u0, RHS, pderatv, V, eta; pdata = pdata, lmaxit=200)
    pcresr = goutrp.reshist
    pcresr /= pcresr[1]
    solrhw = copy(u0)
    solrhw .= pdeptv(goutrp.sol,pdata)
    solldiff=norm(solrhw-sollhw,Inf)
#
    # Solve with right preconditioning
#
    goutrp1 = kl_bicgstab(u0, RHS, pdeatv, V, eta, pdeptv;
             pdata = pdata, lmaxit=200)
    pcresr1 = goutrp1.reshist
    pcresr1 /= pcresr1[1]
    solr = goutrp1.sol
#    solldiff += norm(solr-sollhw,Inf)
    solldiff = max(solldiff,norm(solr-sollhw,Inf))
#
    # Solve with left preconditioning
#
    goutl1 = kl_bicgstab(u0, RHS, pdeatv, V, eta, pdeptv;
             pdata = pdata, lmaxit=200, side="left")
    pcresl1 = goutl1.reshist
    pcresl1 /= pcresl1[1]
    soll = goutl1.sol
#    solldiff += norm(soll-sollhw,Inf)
    solldiff = max(solldiff,norm(soll-sollhw,Inf))
#
# Hardwired and normal give same results?
#
    leftdel= norm(soll-sollhw,Inf) + norm(pcres-pcresl1,Inf)
    leftpass=(leftdel < 1.e-15)
    rightdel = norm(solr-solrhw,Inf) + norm(pcresr1-pcresr,Inf)
    rightpass=(rightdel < 1.e-15)
#
    # Solve with no preconditioning to duplicate fig 3.4 in red book
#
    goutnp = kl_bicgstab(u0, RHS, pdeatv, V, eta; lmaxit=200, pdata = pdata)
    pcresnp = goutnp.reshist
    pcresnp /= pcresnp[1]
    solnone=goutnp.sol
    solldiff = max(solldiff,norm(solnone-sollhw,Inf))
#    solldiff += norm(solnone-sollhw,Inf)
#
# Are the answers close enough?
#
    sollpass=(solldiff < 2.0*eta)
#
# Are the iteration counts correct?
#
    ll=length(pcres)
    lr=length(pcresr1)
    ln=length(pcresnp)
    countok=((ll==7) && (lr==8) && (ln==37))
    pass=sollpass && rightpass && leftpass & countok
    return pass
end

function pdelpatv(u, pdata)
    L = pdata.L
    fdata = pdata.fdata
    au = L * u
    pau = Pfish2d(au, fdata)
    return pau
end

function pderatv(u, pdata)
    L = pdata.L
    fdata = pdata.fdata
    pau = Pfish2d(u, fdata)
    au = L * pau
    return au
end

function pdeptv(u, pdata)
    fdata = pdata.fdata
    ptv = Pfish2d(u, fdata)
end

function pdeatv(u, pdata)
    xc = pdata.xc
    L = pdata.L
    mul!(xc, L, u)
    return xc
end

"""
pdegminit(n)

collects the precomputed data for the linear elliptic pde example.
This is the example on page 54-55 of FR16.

This
includes

- the sparse matrix representation of the operators,
- the right side of the equation,
- the exact solution,
- the data that the fft-based fast Poisson solver (fish2d) needs
"""
function pdegminit(n)
    # Make the grids
    n2 = n * n
    h = 1.0 / (n + 1.0)
    x = collect(h:h:1.0-h)
    o = ones(n)
    Y = o * x'
    y20 = 20.0 * reshape(Y, (n2,))
    DiagY = Diagonal(y20)
    # collect the operators
    D2 = Lap2d(n)
    DX = Dx2d(n)
    DY = Dy2d(n)
    L = D2 + I
    L .+= DX
    LY = copy(DY)
    mul!(LY, DiagY, DY)
    L .+= LY
    # Exact solution and its derivatives
    uexact = solexact(x)
    dxe = dxexact(x)
    dye = dyexact(x)
    d2e = l2dexact(x)
    dxv = reshape(dxe, (n2,))
    dyv = reshape(dye, (n2,))
    d2v = reshape(d2e, (n2,))
    uv = reshape(uexact, (n2,))
    # Preallocate a copy of the unknown for the function
    # and preconditioner evaluation.
    xc = copy(uv)
    fdata = fishinit(n)
    # The right side of the equation
    RHS = d2v + dxv + y20 .* dyv + uv
    # Pack it and ship it.
    pdedata = (L, RHS = RHS, ue = uv, xc = xc, fdata = fdata)
end


text:
struct ImplySparse{
    F1 <: Any,
    F2 <: Any,
    A <: AbstractArray,
    C <: SparseMatrixCSC,
    D <: SparseMatrixCSC,
    E <: SparseMatrixCSC,
    S <: Array{Float64},
    A2 <: Union{Nothing, AbstractArray}} <: Imply
    imp_fun_invia::F1
    imp_fun_S::F2
    F::C
    invia::D
    S::E
    imp_cov::A
    start_val::S
    imp_mean::A2
end

function ImplySparse(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val;
        M = nothing
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC
            }

    invia = neumann_series(A)

    #imp_cov_sym = F*invia*S*permutedims(invia)*permutedims(F)
    #imp_cov_sym = Array(imp_cov_sym)
    invia .= ModelingToolkit.simplify.(invia)

    imp_fun_invia_, imp_fun_invia =
        eval.(ModelingToolkit.build_function(
            invia,
            parameters
            ))

    imp_fun_S_, imp_fun_S =
        eval.(ModelingToolkit.build_function(
            S,
            parameters
            ))

    invia_pre = Base.invokelatest(imp_fun_invia_, start_val)
    S_pre = Base.invokelatest(imp_fun_S_, start_val)
    imp_cov = rand(size(F)[1], size(F)[1])

    if isnothing(M)
        imp_mean = nothing
    else
    end

    return ImplySparse( imp_fun_invia,
                        imp_fun_S,
                        copy(F),
                        invia_pre,
                        S_pre,
                        imp_cov,
                        copy(start_val),
                        imp_mean)
end

function (imply::ImplySparse)(parameters, model)
    imply.imp_fun_invia(imply.invia, parameters)
    imply.imp_fun_S(imply.S, parameters)
    imp_cov = Matrix(imply.invia*imply.S*transpose(imply.invia))
    nvar = size(imply.F, 1)
    imply.imp_cov .= imp_cov[1:nvar, 1:nvar]
end

#################
struct ImplyCommon{
    C <: Real,
    D <: AbstractArray,
    E <: AbstractArray,
    S <: Array{Float64},
    A2 <: Union{Nothing, AbstractArray},
    A <: AbstractArray,
    F1 <: Any,
    F2 <: Any} <: Imply
    n_var::C
    imp_fun_A::F1
    imp_fun_S::F2
    A::D
    S::E
    imp_cov::A
    start_val::S
    imp_mean::A2
end

function ImplyCommon(
        A::Spa1,
        S::Spa2,
        F::Spa3,
        parameters,
        start_val;
        M = nothing
            ) where {
            Spa1 <: SparseMatrixCSC,
            Spa2 <: SparseMatrixCSC,
            Spa3 <: SparseMatrixCSC
            }

    A = Matrix(I-A)
    S = Matrix(S)

    imp_fun_A_, imp_fun_A =
        eval.(ModelingToolkit.build_function(
            A,
            parameters
            ))

    imp_fun_S_, imp_fun_S =
        eval.(ModelingToolkit.build_function(
            S,
            parameters
            ))

    A_pre = Base.invokelatest(imp_fun_A_, start_val)
    S_pre = Base.invokelatest(imp_fun_S_, start_val)
    imp_cov = rand(size(F)[1], size(F)[1])

    n_var = float(size(F, 1))

    if isnothing(M)
        imp_mean = nothing
    else
    end

    return ImplyCommon( n_var,
                        imp_fun_A,
                        imp_fun_S,
                        A_pre,
                        S_pre,
                        imp_cov,
                        copy(start_val),
                        imp_mean)
end

function (imply::ImplyCommon)(parameters, model)
    imply.imp_fun_A(imply.A, parameters)
    imply.imp_fun_S(imply.S, parameters)
    invia = inv(imply.A)
    imp_cov = invia*imply.S*transpose(invia)
    #imply.imp_cov .= imp_cov[1:imply.n_var, 1:imply.n_var]
    copyto!(imply.imp_cov, CartesianIndices(imply.imp_cov), imp_cov, CartesianIndices((imply.imp_cov)))
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 85e3f1ee-9158-11eb-2c00-672502866677
begin
using Pkg
	Pkg.add("JuMP")
	Pkg.add("Juniper")
	Pkg.add("Ipopt") # Juniper + Ipopt = solver for Quadratic optimization problems like this
using JuMP, Ipopt, Juniper
end

# ╔═╡ 2cfafa80-9157-11eb-35c2-51be3c25e814
html"""<style>
main {
    max-width: 1500px;
}"""

# ╔═╡ feb87f00-914f-11eb-2784-8f05d5919eba
md"""# Optimal Weights and Activation Values"""

# ╔═╡ 22f58180-9153-11eb-00c6-6586404966cc
md"""__Data:__  In Sanskrit we observe an alternation between several allomorphic roots for one and the same verb. Let's take _budh_ 'awake, know, $\ldots$' as a concrete example

+ bodh-ati (3sg. pres.)
+ bhot-sya-ti (3sg. fut.)
+ bhut (Noun, Nom.)
+ bhud-bhis (Noun, Instr.)
+ bud-dha (Past participle)

GSC offers a possible solution to this puzzle if we assume that the input for morphological derivations is a blend of different allomorphs with some partial activation values (see GSCNet).

But how do we find suitable values for these activations and how do we find the right constraint weights in a Gradient Harmonic Grammar?
"""

# ╔═╡ 13e8db00-9154-11eb-187a-0d4cd645dcaa
md"""The problem can be phrased as a mathematical optimization problem: We have a function:

$$H_{winner}(x) = c^Tx$$

expressing the Harmony (sum of violations) of the winner, that we want to maximize.
Our $x$s are the Harmonic constraints and the coefficients $c$ are the violations for the winning candidate.
This expression is subject to some constraints: for each other candidate it should hold that

$$H_{winner} > H_{candidate_i} $$ $$ \forall cand_i \in Tableaux$$

We also want the weights to be negative (these are penalties):

$$[x_i, x_2, x_3 ... x_n] \lt 0$$

and we know that the activation values are between 0 and 1: 0 meaning inactive and 1 meaning totally active (discrete).

$$1 \leq [\alpha .... \omega] \geq 0$$

"""

# ╔═╡ 65cb0890-9153-11eb-25b4-b19428afad2c
md""" ### Example: the past participle
__Inputs__:
root allomorphs:

+ (bhud)$^\alpha$
+ (budh)$^\beta$
+ (bhut)$^\gamma$

suffix allomorphs:

+ (-ta)$^\delta$
+ (-dha)$^\epsilon$

__Constraints__:
(all of them $\lt 0$)

+ Max : anti-deletion constraint
+ Dep : anti-epenthesis constraint
+ Voice : Avoid voice contrast ([+ voice][-voice] is bad)
+ License(lar) : Aspirated segments cannot occur before voiceless segments. Also avoid two aspirated next to each other.
+ Lazy : Reduce articulatory effort (avoid aspirated sounds)
+ *h]_W:Aspirated are not allowed in word final position
+ Id(lar) : Aspirated in the input == Aspirated in the output

"""

# ╔═╡ 1f066d00-9153-11eb-00ee-b760c15b6b67
md""" __Toy-example__

| $root_{(\alpha + \beta + \gamma)} + suff_{(\delta + \epsilon)}$ | Max   |   Dep | Voice | Lic | h]W | Lazy | Id(Lar) | Harmony |
|:---:|:---:|:---:|:-----:|:---:|:---:|:---:|:---:|:---:|
|bhud-ta | $\beta + \gamma + \epsilon$ | $(2-\alpha - \delta)$ | 1 | 0 | 0 | 0 | 0 | $M(\beta + \gamma + \epsilon) + D(2-\alpha - \delta) + Voice$ |
|:) bud-dha | $\alpha + \gamma + \delta$ | $(2-\beta - \epsilon)$ | 0 | 0 | 0 | 0 | 1 | $M(\alpha + \gamma + \epsilon) + D(2-\beta - \epsilon) + Id(lar)$ |
|budh-i-ta | $\alpha + \gamma + \epsilon$ | $(3-\beta - \delta)$ | 0 | 0 | 0 | 0 | 1 | $M(\alpha + \gamma + \epsilon) + D(3-\beta - \delta)$ |
| ... |... |... |... |... |... |... |... |... |

Notice that candidate 3 (*budh-i-ta*, built on the model of *kridh-i-ta*) has an extra epenthesis, hence the violation for dep (**3** - $\beta - \delta$)

Let's now solve this problem using the _JuMP_ library for mathematical optimization.
The expression we want to maximize is the Harmony of the winner (candidate 2):

$$max: M(\alpha + \gamma + \epsilon) + D(2-\beta - \epsilon) + Id(lar)$$

subject to:

$M(\alpha + \gamma + \epsilon) + D(2-\beta - \epsilon) + Id(lar) > M(\beta + \gamma + \epsilon) + D(2-\alpha - \delta) + Voice$
$M(\alpha + \gamma + \epsilon) + D(2-\beta - \epsilon) + Id(lar) > M(\alpha + \gamma + \epsilon) + D(3-\beta - \delta)$
$0 \lt [\alpha, \beta, \gamma, \delta, \epsilon] \lt 1$
$M, D, Id(lar), Voice, Lic, h]W < 0$

"""

# ╔═╡ 79976060-915a-11eb-325f-916ca3326950
# Initialize Model
begin
	optimizer = Juniper.Optimizer
	nl_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)
	model = Model(optimizer_with_attributes(optimizer, "nl_solver" => nl_solver))
end

# ╔═╡ 75b962e0-915a-11eb-1b24-5d35a2b179e2
begin
	# Declare variables
	a = 1 # The difference between winner and losers should be greater than a
	b = -1 # Start value for the constraint weights
end

# ╔═╡ 4d9bf512-915b-11eb-2fee-e9a1ff31c9b7
begin
	unregister(model, :M)
	unregister(model, :D)
	unregister(model, :A)
	unregister(model, :L)
	unregister(model, :Li)
	unregister(model, :Coda)
	unregister(model, :Id)
	unregister(model, :α)
	unregister(model, :β)
	unregister(model, :γ)
	unregister(model, :δ)
	unregister(model, :ϵ)
end

# ╔═╡ 1bf6e310-9153-11eb-2fe5-5b4d277cff3a
begin
	@variable(model,  -10 <= Ma <= -1, start = b,integer = true) # MaxIO
	@variable(model,  -10 <= De <= -1, start = b,integer = true) # DepIO
	@variable(model,  -10 <= V <= -1, start = b,integer = true) # Assimilate == Voice
	@variable(model,  -10 <= Lic <= -1, start = b,integer = true) # License(lar)
	@variable(model, -10 <= hW <= -1, start = b,integer = true) # h]W
	@variable(model, -10 <= La <= -1, start = b,integer = true) #Lazy
	@variable(model, -10 <= IdL <= -3, start = 2 * b,integer = true) # Ident(lar)
	@variable(model, 0.1 <= α <= 0.9, start = 0.2) # bhud
	@variable(model, 0.1 <= β <= 0.9, start = 0.9) # budh
	@variable(model, 0.1 <= γ <= 0.9, start = 0.4) # bhut
	@variable(model, 0.1 <= δ <= 0.5, start = 0.5) # dha
	@variable(model, 0.5 <= ϵ <= 0.9, start = 0.7) # ta
end

# ╔═╡ 6d36f150-915f-11eb-252a-17639f320ed2
@objective(model, Max, Ma * (α + γ + δ) + De * (2 - β - ϵ) + IdL)

# ╔═╡ 900302e0-915b-11eb-3b97-855f0f06e4ba
# Objective function : Harmony of the winner
"""@objective(model, Max, Ma * (α + γ + δ) + De * (2 - β - ϵ) + IdL)"""

# ╔═╡ 8856e662-915b-11eb-1903-37e8e9ef1ed0
# Constraints
begin
	# Winner Harmony
	winner = @expression(model, Ma * (α + γ + δ) + De * (2 - β - ϵ) + IdL)

	bhudta = @expression(model, winner - (Ma * (β + γ + δ) + De * (2 - α - ϵ) + V))
	@constraint(model, bhudta >= a)

	budhta = @expression(model, winner - (Ma * (α + γ + δ) + De * (2 - β - ϵ) + V + Lic))
	@constraint(model, budhta >= a)

	budhita = @expression(model, winner - (Ma * (α + γ + δ) + De * (3 - β - ϵ)))
	@constraint(model, budhita >= a)
end

# ╔═╡ 6774d1e0-915c-11eb-1070-e98f93030697
"""$model"""

# ╔═╡ 7d023ed0-915c-11eb-25dc-e538e4c8f878
# Optimize / Solve
optimize!(model)

# ╔═╡ 7ae6d8e0-915c-11eb-39a3-bbfbfb998e5a
# Printing logic
# Objective value (Harmony of the winner after the optimization)
JuMP.objective_value(model)

# ╔═╡ 830bb260-915d-11eb-1afa-ffdeaa98ad5a
# Solver status
JuMP.primal_status(model)

# ╔═╡ 80785d50-915d-11eb-1441-2f013038b3e0
JuMP.termination_status(model)

# ╔═╡ 745246e0-915c-11eb-14f6-593066f40ef8
# MaxIO
value(Ma)

# ╔═╡ 6ecf24e0-915c-11eb-28cc-7d5e29ba29f1
#Dep
value(De)

# ╔═╡ b747204e-915d-11eb-1d2b-b77c4386f7f5
# Voice
value(V)

# ╔═╡ c1afebce-915d-11eb-3ba2-71fa0cc58aa1
# Coda
value(hW)

# ╔═╡ c20dea50-915d-11eb-2046-bb9095f0f080
value(α)

# ╔═╡ c22f06e0-915d-11eb-22da-3d03814961cd
value(β)

# ╔═╡ fcbe5450-915d-11eb-1080-c91bfcde29b5
value(γ)

# ╔═╡ 0187ea52-915e-11eb-2475-11a14ba0bc20
value(δ)

# ╔═╡ 062ce150-915e-11eb-1d7b-7558b1067064
value(ϵ)

# ╔═╡ ff36a752-915d-11eb-0835-33c93eaa8233
md"Of course this was only a toy-example. JuMP evaluated the expressions for just two candidates and some value were arbitrarly fixed but I think it gives the idea, how Mathematical Optimization (Quadratic Programming in this case) can be used to search suitable values for the constraints and the activation values"

# ╔═╡ Cell order:
# ╟─2cfafa80-9157-11eb-35c2-51be3c25e814
# ╟─feb87f00-914f-11eb-2784-8f05d5919eba
# ╟─22f58180-9153-11eb-00c6-6586404966cc
# ╟─13e8db00-9154-11eb-187a-0d4cd645dcaa
# ╟─65cb0890-9153-11eb-25b4-b19428afad2c
# ╟─1f066d00-9153-11eb-00ee-b760c15b6b67
# ╟─85e3f1ee-9158-11eb-2c00-672502866677
# ╠═79976060-915a-11eb-325f-916ca3326950
# ╠═75b962e0-915a-11eb-1b24-5d35a2b179e2
# ╟─4d9bf512-915b-11eb-2fee-e9a1ff31c9b7
# ╠═1bf6e310-9153-11eb-2fe5-5b4d277cff3a
# ╠═6d36f150-915f-11eb-252a-17639f320ed2
# ╠═900302e0-915b-11eb-3b97-855f0f06e4ba
# ╠═8856e662-915b-11eb-1903-37e8e9ef1ed0
# ╠═6774d1e0-915c-11eb-1070-e98f93030697
# ╠═7d023ed0-915c-11eb-25dc-e538e4c8f878
# ╠═7ae6d8e0-915c-11eb-39a3-bbfbfb998e5a
# ╠═830bb260-915d-11eb-1afa-ffdeaa98ad5a
# ╠═80785d50-915d-11eb-1441-2f013038b3e0
# ╠═745246e0-915c-11eb-14f6-593066f40ef8
# ╠═6ecf24e0-915c-11eb-28cc-7d5e29ba29f1
# ╠═b747204e-915d-11eb-1d2b-b77c4386f7f5
# ╠═c1afebce-915d-11eb-3ba2-71fa0cc58aa1
# ╠═c20dea50-915d-11eb-2046-bb9095f0f080
# ╠═c22f06e0-915d-11eb-22da-3d03814961cd
# ╠═fcbe5450-915d-11eb-1080-c91bfcde29b5
# ╠═0187ea52-915e-11eb-2475-11a14ba0bc20
# ╠═062ce150-915e-11eb-1d7b-7558b1067064
# ╟─ff36a752-915d-11eb-0835-33c93eaa8233

max_stars_repo_name: amasotti/grascom

max_issues_repo_name: amasotti/grascom

max_forks_repo_name: amasotti/grascom