text:
const MJSTATE_FIELDS = (:time, :qpos, :qvel, :act, :mocap_pos, :mocap_quat, :userdata, :qacc_warmstart)
const MJACTION_FIELDS = (:ctrl, :qfrc_applied, :xfrc_applied)

const DEFAULT_SKIP = 1
const DEFAULT_CAN_SHARE_MODEL = true

"""
    MJSim

The `MJSim` type couples a `jlModel` and `jlData` from MuJoCo.jl to provide a full simulation.

The following are the official/internal/minimum set of fields from `jlData` for
state, observation, and action in MuJoCo:

- State: `(time, qpos, qvel, act, mocap_pos, mocap_quat, userdata, qacc_warmstart)`
- Observation: `sensordata`
- Action: `(ctrl, qfrc_applied, xfrc_applied)`

MJSim follows this definition except for actions (e.g. `setaction!``), which is
composed of just `ctrl` by default.

For more information, see the "State and control" section of the
[MuJoCo Documentation](www.mujoco.org/book/programming.html)

# Fields

- `m::jlModel`: contains the model description and is expected to remain constant.
- `d::jlData`: contains all dynamic variables and intermediate results.
- `mn::Tuple`: named-access version of `MJSim.m` provided by AxisArrays.jl.
- `dn::Tuple`: named-access version of `MJSim.d` provided by AxisArrays.jl.
- `initstate::Vector{mjtNum}`: The initial state vector at the time when this MJSim was constructed.
- `skip::Int`: the number of times the simulation is integrated, yielding an
  effective simulation timestep of `skip * m.opt.timestep`.
"""
struct MJSim{MN, DN, S, SE, A}
    m::jlModel
    d::jlData
    mn::MN
    dn::DN
    initstate::Vector{mjtNum}
    skip::Int

    # internal
    statespace::S
    observationspace::SE
    actionspace::A

    @doc """
        $(TYPEDSIGNATURES)

    Construct an MJSim from `m` and `d`, with skip `skip`.
    """
    function MJSim(m::jlModel, d::jlData; skip::Integer = DEFAULT_SKIP)
        check_skip(skip)

        m_named, d_named = namify(m, d)

        nameshapes = map(MJSTATE_FIELDS) do name
            field = getproperty(d, name)
            name => Shape(eltype(field), size(field)...)
        end
        statespace = MultiShape(nameshapes...)

        if m.nsensor > 0
            nameshapes = map(1:m.nsensor) do id
              name = jl_id2name(m, MJCore.mjOBJ_SENSOR, id)
              name = isnothing(name) ? Symbol("unnamed_$id") : Symbol(name)
              dof = m.sensor_dim[id]
              shape = dof == 1 ? ScalarShape(mjtNum) : VectorShape(mjtNum, dof)
              name => shape
            end
            observationspace = MultiShape(nameshapes...)
        else
            observationspace = VectorShape(mjtNum, 0)
        end

        if m.nu > 0
            clampctrl = !jl_disabled(m, MJCore.mjDSBL_CLAMPCTRL)
            nameshapes = map(1:m.nu) do id
                name = jl_id2name(m, MJCore.mjOBJ_ACTUATOR, id)
                name = isnothing(name) ? Symbol("unnamed_$id") : Symbol(name)
                name => ScalarShape(mjtNum) #TODO boundedshape
            end
            actionspace = MultiShape(nameshapes...)
        else
            actionspace = VectorShape(mjtNum, 0)
        end

        initstate = mapreduce(vcat, MJSTATE_FIELDS) do name
            field = getproperty(d, name)
            ndims(field) > 1 ? vec(field) : field
        end

    	sim = new{typeof(m_named), typeof(d_named), typeof(statespace), typeof(observationspace), typeof(actionspace)}(m, d, m_named, d_named,
            initstate, skip,
            statespace, observationspace, actionspace)

        forward!(sim)
    end
end

"""
    $(TYPEDSIGNATURES)

Construct an `MJSim` the MJCF or MJB model file located at `modelpath`, with skip `skip`.
"""
function MJSim(modelpath::String; skip::Integer = DEFAULT_SKIP)
    m = jlModel(modelpath)
    d = jlData(m)
    MJSim(m, d, skip=skip)
end

function LyceumBase.tconstruct(::Type{MJSim}, N::Integer, m::jlModel; skip::Integer = DEFAULT_SKIP, can_share_model::Bool=DEFAULT_CAN_SHARE_MODEL)
    N > 0 || throw(ArgumentError("N must be > 0"))
    if can_share_model
        ntuple(_ -> MJSim(m, jlData(m), skip=skip), N)
    else
        ntuple(N) do _
            mcopy = copy(m)
            MJSim(mcopy, jlData(mcopy), skip=skip)
        end
    end
end

function LyceumBase.tconstruct(::Type{MJSim}, N::Integer, modelpath::AbstractString; kwargs...)
    tconstruct(MJSim, N, jlModel(modelpath); kwargs...)
end

"""
    $(TYPEDSIGNATURES)

Return a description of `sim`'s statespace.
"""
@inline LyceumBase.statespace(sim::MJSim) = sim.statespace

"""
    $(SIGNATURES)

Copy the following state fields from `sim.d` into `state`:
```
($(join(map(string, MJSTATE_FIELDS), ", ")))
```
"""
@propagate_inbounds function LyceumBase.getstate!(state::RealVec, sim::MJSim)
    @boundscheck checkaxes(statespace(sim), state)
    shaped = statespace(sim)(state)
    @uviews shaped begin _copyshaped!(shaped, sim.d) end
    state
end

"""
    $(SIGNATURES)

Return a flattened vector of the following state fields from `sim.d`:
```
($(join(map(string, MJSTATE_FIELDS), ", ")))
```
"""
@propagate_inbounds LyceumBase.getstate(sim::MJSim) = getstate!(allocate(statespace(sim)), sim)

# Two things to note about MuJoCo state:
# 1. MJSTATE_FIELDS are theoretically the minimal fields/sufficient statistics, as in
#    all other fields in mjData are a function of MJSTATE_FIELDS, but this is only true for
#    **dynamic** elements. For example, calling rand!(d.xpos) followed by mj_forward(m, d) will
#    overwrite the entries of d.xpos corresponding to dynamically changing elements, but NOT for
#    e.g. the root world body (d.xpos[:, 1]). If the user inadvertantly overwrites these fields,
#    mjData will be left in a corrupt state. Because of this, we call mj_resetData(m, d) followed
#    by mj_forward(m, d) inside the exported setstate!, while providing the ability to
#    skip this using copystate! for users who know what they're doing and want to squeeze some
#    extra performance out.
# 2. qacc_warmstart is part of state in that it will control the final solution (because MuJoCo's
#    optimizer may not find the exact minima). This means that:
#       setstate!(sim, state)
#       @assert getstate(sim) == state
#    will fail if the optimizer terminated early. To avoid this, we first copy over the state
#    variables (MJSTATE_FIELDS, including qacc_warmstart), call mj_forward, and then re-copy
#    over qacc_warmstart so that the above equality holds true. This means that mjData is left
#    in a slightly inconsitent state (qacc_warmstart will be slightly out of sync).

"""
    $(SIGNATURES)

Copy the components of `state` to their respective fields in `sim.d`, namely:
```
($(join(map(string, MJSTATE_FIELDS), ", ")))
```
"""
@propagate_inbounds function LyceumBase.setstate!(sim::MJSim, state::RealVec)
    mj_resetData(sim.m, sim.d)
    copystate!(sim, state)
    forward!(sim)
    shaped = statespace(sim)(state)
    @uviews shaped begin copyto!(sim.d.qacc_warmstart, shaped.qacc_warmstart) end
    sim
end

@propagate_inbounds function copystate!(sim::MJSim, state::RealVec)
    @boundscheck checkaxes(statespace(sim), state)
    shaped = @inbounds statespace(sim)(state)
    @uviews shaped begin _copyshaped!(sim.d, shaped) end
    sim
end

@inline function _copyshaped!(dst, src)
    @inbounds dst.time = src.time
    @inbounds copyto!(dst.qpos, src.qpos)
    @inbounds copyto!(dst.qvel, src.qvel)
    @inbounds copyto!(dst.act, src.act)
    @inbounds copyto!(dst.mocap_pos, src.mocap_pos)
    @inbounds copyto!(dst.mocap_quat, src.mocap_quat)
    @inbounds copyto!(dst.userdata, src.userdata)
    @inbounds copyto!(dst.qacc_warmstart, src.qacc_warmstart)
    dst
end

"""
    $(TYPEDSIGNATURES)

Return a description of `sim`'s observation space.
"""
@inline LyceumBase.observationspace(sim::MJSim) = sim.observationspace

"""
    $(SIGNATURES)

Copy `sim.d.sensordata` into `obs`.
"""
@propagate_inbounds function LyceumBase.getobservation!(obs::RealVec, sim::MJSim)
    @boundscheck checkaxes(observationspace(sim), obs)
    @inbounds copyto!(obs, sim.d.sensordata)
end

"""
    $(TYPEDSIGNATURES)

Return a copy of `sim.d.sensordata`.
"""
@propagate_inbounds LyceumBase.getobservation(sim::MJSim) = getobservation!(allocate(sim.observationspace), sim)

"""
    $(TYPEDSIGNATURES)

Return a description of `sim`'s action space.
"""
@inline LyceumBase.actionspace(sim::MJSim) = sim.actionspace

"""
    $(TYPEDSIGNATURES)

Copy `sim.d.ctrl` into `action`.
"""
@propagate_inbounds function LyceumBase.getaction!(action::RealVec, sim::MJSim)
    @boundscheck checkaxes(actionspace(sim), action)
    @inbounds copyto!(action, sim.d.ctrl)
end

"""
    $(TYPEDSIGNATURES)

Return a copy of `sim.d.ctrl`.
"""
@propagate_inbounds LyceumBase.getaction(sim::MJSim, action::RealVec) = getaction!(allocate(sim.actionspace), sim)

"""
    $(TYPEDSIGNATURES)

Copy `action` into `sim.d.ctrl` into `action` and compute the new forward dynamics.
"""
@propagate_inbounds function LyceumBase.setaction!(sim::MJSim, action::RealVec)
    forwardskip!(setaction_nofwd!(sim, action), MJCore.mjSTAGE_VEL)
end

@propagate_inbounds function setaction_nofwd!(sim::MJSim, action::RealVec)
    @boundscheck checkaxes(actionspace(sim), action)
    @inbounds copyto!(sim.d.ctrl, action)
    sim
end

@inline LyceumBase.reset!(sim::MJSim) = forward!(reset_nofwd!(sim))
@inline reset_nofwd!(sim::MJSim) = (mj_resetData(sim.m, sim.d); sim)

# typically 2-3x faster than reset!
@inline fastreset!(sim::MJSim) = forward!(fastreset_nofwd!(sim))
@inline function fastreset_nofwd!(sim::MJSim)
    zerofullctrl_nofwd!(sim)
    @inbounds copystate!(sim, sim.initstate)
end

"""
    $(TYPEDSIGNATURES)

Step the simulation by `skip` steps, where `skip` defaults to `sim.skip`.

State-dependent controls (e.g. the `ctrl`, `xfrc_applied`, `qfrc_applied`
fields of `sim.d`) should be set before calling `step!`.
"""
function LyceumBase.step!(sim::MJSim, skip::Integer=sim.skip)
    check_skip(skip)
    # According to MuJoCo docs order must be:
    # 1. mj_step1
    # 2. set_{ctrl, xfrc_applied, qfrc_applied}
    # 3. mj_step2
    # This implies user should call setaction! then step!
    for _=1:skip
        mj_step2(sim.m, sim.d)
        mj_step1(sim.m, sim.d)
    end
    sim
end

"""
    $(TYPEDSIGNATURES)

Zero out `sim.d.ctrl` and compute new forward dynamics.
"""
@inline zeroctrl!(sim::MJSim) = forwardskip!(zeroctrl_nofwd!(sim), MJCore.mjSTAGE_VEL)

@inline zeroctrl_nofwd!(sim::MJSim) = (fill!(sim.d.ctrl, zero(mjtNum)); sim)

"""
    $(TYPEDSIGNATURES)

Zero out all of the fields in `sim.d` that contribute to forward dynamics calculations,
namely `ctrl`, `qfrc_applied`, and `xfrc_applied`, and compute the forward dynamics.
"""
@inline zerofullctrl!(sim::MJSim) = forwardskip!(zerofullctrl_nofwd!(sim), MJCore.mjSTAGE_VEL)

@inline function zerofullctrl_nofwd!(sim::MJSim)
    fill!(sim.d.ctrl, zero(mjtNum))
    fill!(sim.d.qfrc_applied, zero(mjtNum))
    fill!(sim.d.xfrc_applied, zero(mjtNum))
    sim
end

"""
    $(TYPEDSIGNATURES)

Calculate the center of mass of all the bodies in simulation `sim` by computing the
mass-weighted sum of the Cartesian coordinates of each body in the world frame.
"""
function masscenter(sim::MJSim)
    mcntr = zeros(SVector{3, Float64})
    mtotal = 0.0
    body_mass, xipos = sim.m.body_mass, sim.d.xipos
    @uviews xipos for i=1:sim.m.nbody
        bmass = body_mass[i]
        bodycom = SVector{3, Float64}(uview(xipos, :, i))
        mcntr = mcntr + bmass * bodycom
        mtotal += bmass
    end
    mcntr = SVector(mcntr)
    mcntr / mtotal
end

"""
    $(TYPEDSIGNATURES)

Compute the forward dynamics of the simulation and store them in `sim.d`.
Equivalent to `mj_forward(sim.m, sim.d)`.
"""
@inline forward!(sim::MJSim) = (mj_forward(sim.m, sim.d); sim)

"""
    $(TYPEDSIGNATURES)

Compute the forward dynamics of the simulation and store them in `sim.d`, optionally
skipping parts of the dynamics calculation. Equivalent to `mj_forwardSkip(sim.m, sim.d,
skipstage, skipsensor)`.
"""
@inline function forwardskip!(sim::MJSim, skipstage::MJCore.mjtStage=MJCore.mjSTAGE_NONE, skipsensor::Bool=false)
    mj_forwardSkip(sim.m, sim.d, skipstage, skipsensor)
    sim
end

"""
    $(TYPEDSIGNATURES)

Return the effective timestep of `sim`. Equivalent to `sim.skip * sim.m.opt.timestep`.
"""
@inline LyceumBase.timestep(sim::MJSim) = sim.m.opt.timestep * sim.skip

"""
    $(TYPEDSIGNATURES)

Return the current simulation time, in seconds, of `sim`. Equivalent to sim.d.time.
"""
@inline Base.time(sim::MJSim) = sim.d.time

@inline getsim(sim::MJSim) = sim

# TODO(cxs):: Make more informative. This is a temp fix for those gawdawful error messages from sim.{mn, dn}.
Base.show(io::IO, ::MIME"text/plain", sim::Union{<:MJSim, Type{<:MJSim}}) = show(io, sim)
Base.show(io::IO, sim::Union{<:MJSim, Type{<:MJSim}}) = print(io, "MJSim")

@inline check_skip(skip) = skip > 0 || throw(ArgumentError("`skip` must be > 0"))


text:
"""
    psconvert(cmd0::String="", arg1=nothing; kwargs...)

Place images or EPS files on maps.

Full option list at [`psconvert`]($(GMTdoc)psconvert.html)

Parameters
----------

- **A** | **adjust** | **crop** :: [Type => Str or Number]

    Adjust the BoundingBox and HiResBoundingBox to the minimum required by the image content.
    ($(GMTdoc)psconvert.html#a)
- **C** | **gs_option** :: [Type => Str or Array os strings]

    Specify a single, or an araay of, custom option that will be passed on to GhostScript as is.
    ($(GMTdoc)psconvert.html#c)
- **D** | **out_dir** | **output_dir** :: [Type => Str]

    Sets an alternative output directory (which must exist) [Default is the same directory
    as the PS files].
    ($(GMTdoc)psconvert.html#d)
- **E** | **dpi** :: [Type => Number]

    Set raster resolution in dpi [default = 720 for PDF, 300 for others].
    ($(GMTdoc)psconvert.html#e)
- **F** | **:out_name** | **output_name** :: [Type => Str]

    Force the output file name.
    ($(GMTdoc)psconvert.html#f)
- **G** | **ghost_path** :: [Type => Bool]

    Full path to your GhostScript executable.
    ($(GMTdoc)psconvert.html#g)
- **I** | **icc_gray** :: [Type => Bool]

    Enforce gray-shades by using ICC profiles.
    ($(GMTdoc)psconvert.html#i)
- **in_memory** :: [Type => Bool]

    Process a in memory PS file. No other input file should be provided.
    Currently works on Windows only.
- **L** | **list_file** :: [Type => Str]

    The listfile is an ASCII file with the names of the PostScript files to be converted.
    ($(GMTdoc)psconvert.html#l)
- **Q** | **anti_aliasing** :: [Type => Str]

    Set the anti-aliasing options for graphics or text. Append the size of the subsample box
    (1, 2, or 4) [4]. This option is set by default.
    ($(GMTdoc)psconvert.html#q)
- **S** | **gs_command** :: [Type => Bool]

    Print to standard error the GhostScript command after it has been executed.
    ($(GMTdoc)psconvert.html#s)
- **T** | **format** :: [Type => Str]

    b|e|E|f|F|j|g|G|m|s|t Sets the output format, where b = BMP, e = EPS, E = EPS with PageSize command,
    f = PDF, F = multi-page PDF, j = JPEG, g = PNG, G = transparent PNG (untouched regions are
    transparent), m = PPM,  and t = TIFF [default is JPEG].
    Alternatively, the format may be set with the *fmt* keyword, e.g. *fmt=:png*.
    ($(GMTdoc)psconvert.html#t)
- **W** | **world_file** :: [Type => Str]

    Write a ESRI type world file suitable to make (e.g) .tif files be recognized as geotiff by
    software that know how to do it.
    ($(GMTdoc)psconvert.html#w)
- **kml** :: [Type => Str | []]

    Create a minimalist KML file that allows loading the image in GoogleEarth.
    ($(GMTdoc)psconvert.html#w)
- **Z** | **del_input_ps** :: [Type => Bool]

    Remove the input PostScript file(s) after the conversion.
    ($(GMTdoc)psconvert.html#z)
- $(GMT.opt_V)
"""
function psconvert(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("psconvert", cmd0, arg1)

	if (!isempty(cmd0)) && (arg1 === nothing)  arg1 = cmd0  end

	d = init_module(false, kwargs...)[1]		# Also checks if the user wants ONLY the HELP mode
	cmd = add_opt(d, "", "A", [:A :adjust :crop])
	(cmd == " -A") && (cmd = cmd * "1p")			# If just -A default to -A1p
	cmd = parse_these_opts(cmd, d, [[:D :out_dir :output_dir], [:E :dpi], [:F :out_name :output_name],
	                                [:G :ghost_path], [:I :icc_gray], [:L :list_file], [:P :portrait], [:Q :anti_aliasing], [:S :gs_command], [:Z :del_input_ps]])
	cmd = parse_V_params(d, cmd)

	if ((val = find_in_dict(d, [:fmt])[1]) !== nothing)
		fmt = isa(val, Symbol) ? string(val) : val
		if     (fmt == "pdf")  cmd *= " -Tf"
		elseif (fmt == "eps")  cmd *= " -Te"
		elseif (fmt == "png")  cmd *= " -Tg"
		elseif (fmt == "jpg")  cmd *= " -Tj"
		elseif (fmt == "tif")  cmd *= " -Tt"
		end
	else
		cmd = add_opt(d, cmd, "T", [:T :format])
	end

	if ((val = find_in_dict(d, [:C :gs_option])[1]) !== nothing)
		if (isa(val, String) || isa(val, Symbol))
			cmd = string(cmd, " -C", val)
		elseif (isa(val, Array{String}))
			for k = 1:length(val)
				cmd *= " -C" * val[k]
			end
		end
	end

	cmd = add_opt(d, cmd, "W", [:W :world_file])
	(haskey(d, :kml)) && (cmd *= " -W+k" * d[:kml])

	if (haskey(d, :in_memory))
		(arg1 === nothing) ? cmd *= " =" : arg1 = nothing	# in_memory and input file are incompat. File wins
		delete!(d, :in_memory)
	end

	# In case DATA holds a file name, copy it into cmd.
	if (cmd0 != "" || arg1 !== nothing)						# Data was passed as file name
		cmd, got_fname, arg1 = find_data(d, cmd0, cmd, arg1)		# Find how data was transmitted
	end

	if (isempty(cmd))          cmd = "-A1p -Tj -Qg4 -Qt4"  end 	# Means no options were used. Allowed case
	if (!occursin("-Q", cmd))  cmd = cmd * " -Qt4 -Qg4"    end	# We promised to have these as default

	if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end
	gmt("psconvert " * cmd, arg1)
end

# ---------------------------------------------------------------------------------------------------
psconvert(arg1::GMTps; kw...) = psconvert("", arg1; kw...)

ext: jl


text:
using SeisIO, SeisIO.Quake
import Printf
# US FDSNget example: 5 stations, 2 networks, all channels, last 600 seconds
println(stdout, "Beginning real \"use case\" examples...\n\n",
  "Note: ALL data from these examples can be written to file; \n",
  "type `wseis(object_name)` to save to native SeisIO format, \n",
  "type `writesac(object_name)` to write to SAC files,\n",
  "or rerun the data request with added keyword \"w=true\" to\n",
  "write mseed.")

println(stdout, "FDSN: see structure seis_fdsn")
CHA = "CC.PALM, UW.HOOD, CC.TIMB, CC.HIYU, UW.TDH"
s = -600
t = u2d(time())
println(stdout, "Command: seis_fdsn = get_data(\"FDSN\", \"", CHA, "\", src=\"IRIS\", s=", s, ", t=", t)
seis_fdsn = get_data("FDSN", CHA, src="IRIS", s=s, t=t)

# IRISWS example: 6 channels, 30 minutes, synchronized, saved to SAC format"
println(stdout, "IRISWS timeseries: see structure seis_iris")
STA = ["CC.TIMB..EHE", "CC.TIMB..EHN", "CC.TIMB..EHZ", "UW.HOOD..HHE", "UW.HOOD..HHN", "UW.HOOD..HHZ"]
st = -3600
en = -1800
println(stdout, "Command: seis_iris = get_data(\"IRIS\", ", STA, ", s=", st, ", t=", en)
seis_iris = get_data("IRIS", STA, s=st, t=en)

# The Tohoku-Oki great earthquake, from IRIS FDSN, recorded by boreholes in WA (USA)
println(stdout, "FDSN event example: see structure seis_evt")
println(stdout, "seis_evt = FDSNevt(\"201103110547\", \"PB.B004..EH?,PB.B004..BS?,PB.B001..BS?,PB.B001..EH?\")")
seis_evt = FDSNevt("201103110547", "PB.B004..EH?,PB.B004..BS?,PB.B001..BS?,PB.B001..EH?")
show(seis_evt)

# IRIS SeedLink session in TIME mode
println(stdout, "SeedLink example (part 1): see structure seis_sl")
sta = "UW.GRUT,UW.H1K,UW.MDW"
s1 = -120
t1 = 120
println(stdout, "seis_sl = seedlink(\"", sta, "\", mode=\"TIME\", s=", s1, ", t=", t1, ")")
seis_sl = seedlink(sta, mode="TIME", s=s1, t=t1)

# IRIS SeedLink session in DATA mode, to same structure
println(stdout, "SeedLink example (part 2): adds to structure seis_sl")
println(stdout, "SeedLink!(seis_sl, \"SampleFiles/SL_long_test.conf\", mode=\"DATA\")")
seedlink!(seis_sl, "SampleFiles/SL_long_test.conf", mode="DATA")
println(stdout, "When finished, close connections with command \"for conn in seis_sl.c; close(conn); end\"")


text:
function BodiesElements(arg0::AbsoluteDate, arg1::jdouble, arg10::jdouble, arg11::jdouble, arg12::jdouble, arg13::jdouble, arg14::jdouble, arg15::jdouble, arg16::jdouble, arg17::jdouble, arg18::jdouble, arg19::jdouble, arg2::jdouble, arg20::jdouble, arg21::jdouble, arg22::jdouble, arg23::jdouble, arg24::jdouble, arg25::jdouble, arg26::jdouble, arg27::jdouble, arg28::jdouble, arg29::jdouble, arg3::jdouble, arg30::jdouble, arg31::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::jdouble)
    return BodiesElements((AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg10, arg11, arg12, arg13, arg14, arg15, arg16, arg17, arg18, arg19, arg2, arg20, arg21, arg22, arg23, arg24, arg25, arg26, arg27, arg28, arg29, arg3, arg30, arg31, arg4, arg5, arg6, arg7, arg8, arg9)
end

function get_d(obj::DelaunayArguments)
    return jcall(obj, "getD", jdouble, ())
end

function get_d_dot(obj::DelaunayArguments)
    return jcall(obj, "getDDot", jdouble, ())
end

function get_date(obj::DelaunayArguments)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_f(obj::DelaunayArguments)
    return jcall(obj, "getF", jdouble, ())
end

function get_f_dot(obj::DelaunayArguments)
    return jcall(obj, "getFDot", jdouble, ())
end

function get_gamma(obj::DelaunayArguments)
    return jcall(obj, "getGamma", jdouble, ())
end

function get_gamma_dot(obj::DelaunayArguments)
    return jcall(obj, "getGammaDot", jdouble, ())
end

function get_l(obj::DelaunayArguments)
    return jcall(obj, "getL", jdouble, ())
end

function get_l_dot(obj::DelaunayArguments)
    return jcall(obj, "getLDot", jdouble, ())
end

function get_l_ju(obj::BodiesElements)
    return jcall(obj, "getLJu", jdouble, ())
end

function get_l_ju_dot(obj::BodiesElements)
    return jcall(obj, "getLJuDot", jdouble, ())
end

function get_l_ma(obj::BodiesElements)
    return jcall(obj, "getLMa", jdouble, ())
end

function get_l_ma_dot(obj::BodiesElements)
    return jcall(obj, "getLMaDot", jdouble, ())
end

function get_l_me(obj::BodiesElements)
    return jcall(obj, "getLMe", jdouble, ())
end

function get_l_me_dot(obj::BodiesElements)
    return jcall(obj, "getLMeDot", jdouble, ())
end

function get_l_ne(obj::BodiesElements)
    return jcall(obj, "getLNe", jdouble, ())
end

function get_l_ne_dot(obj::BodiesElements)
    return jcall(obj, "getLNeDot", jdouble, ())
end

function get_l_prime(obj::DelaunayArguments)
    return jcall(obj, "getLPrime", jdouble, ())
end

function get_l_prime_dot(obj::DelaunayArguments)
    return jcall(obj, "getLPrimeDot", jdouble, ())
end

function get_l_sa(obj::BodiesElements)
    return jcall(obj, "getLSa", jdouble, ())
end

function get_l_sa_dot(obj::BodiesElements)
    return jcall(obj, "getLSaDot", jdouble, ())
end

function get_l_ur(obj::BodiesElements)
    return jcall(obj, "getLUr", jdouble, ())
end

function get_l_ur_dot(obj::BodiesElements)
    return jcall(obj, "getLUrDot", jdouble, ())
end

function get_l_ve(obj::BodiesElements)
    return jcall(obj, "getLVe", jdouble, ())
end

function get_l_ve_dot(obj::BodiesElements)
    return jcall(obj, "getLVeDot", jdouble, ())
end

function get_le(obj::BodiesElements)
    return jcall(obj, "getLE", jdouble, ())
end

function get_le_dot(obj::BodiesElements)
    return jcall(obj, "getLEDot", jdouble, ())
end

function get_omega(obj::DelaunayArguments)
    return jcall(obj, "getOmega", jdouble, ())
end

function get_omega_dot(obj::DelaunayArguments)
    return jcall(obj, "getOmegaDot", jdouble, ())
end

function get_pa(obj::BodiesElements)
    return jcall(obj, "getPa", jdouble, ())
end

function get_pa_dot(obj::BodiesElements)
    return jcall(obj, "getPaDot", jdouble, ())
end

function get_tc(obj::DelaunayArguments)
    return jcall(obj, "getTC", jdouble, ())
end


text:
# FitzHugh_Nagumo_2D_concentric_plane_spiral_waves_v1.jl
# ======================================================   u.p. 4.8.21

# only concetric and spiral waves - no linear wave

using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot, OrdinaryDiffEq
using Random

function FitzHugh_Nagumo_ODE(du,u,p,t)
    a, b, d, epsilon, hsq6 = p
    uu = @view u[1,:,:]
    vv = @view u[2,:,:]

    # This code can be optimized a lot if time allows
    # Create padded u matrix to incorporate Newman boundary conditions - 9-point stencil
    uuu=[ [uu[2,2] uu[2,:]' uu[2,end-1] ] ; [ uu[:,2] uu uu[:,end-1] ] ; [ uu[end-1,2] uu[end-1,:]' uu[end-1,end-1] ] ]
    diff_term = d .* (
        4 .* (uuu[2:end-1,1:end-2] .+ uuu[2:end-1,3:end] .+
            uuu[3:end,2:end-1] .+ uuu[1:end-2,2:end-1] ) .+
            uuu[3:end,3:end] .+ uuu[3:end,1:end-2] .+ uuu[1:end-2,3:end] .+
            uuu[1:end-2,1:end-2] .- 20 .*uu
    ) ./ hsq6

    du[1,:,:] = a .* uu .*(1 .-uu).*(uu.-b) .- vv .+ diff_term
    du[2,:,:] = epsilon .* (uu .- vv )
end

# --------------------------  main program  -----------------------------

@time begin

# parameters
L = 300   # size of domain
N = 150   # no. of grid points
h = L/N # spatial steps size
hsq6 = 6*h*h

# FHN parameters
a = 3.
b = 0.2
d = 1.   #  diffusion
epsilon = 0.01

# open figure
kplt = 0  # no. of subplot

allsols = []
allts = []

for ic = 1:2    # initial conditions ( = rows of the plot )

println("ic=",ic)

uu = zeros(N,N)
vv = zeros(N,N)
Random.seed!(13371)
uu = uu .+ 0.01*randn(N,N)    # random perturbation

# local stimulius = concentric waves
if ic == 1
   uu[40:41,40:41] .= 0.999
   uu[75:76,75:76] .= 0.999
   tspan = (0., 300.)
   tsamp = [0., 90., 150., 300.]
end

# local stimulius = spiral waves
if ic == 2
   uu[1:end,10:11] .= 0.999
end

# initial values
u0 = zeros(2,N,N)
u0[1,:,:] = uu
u0[2,:,:] = vv

p = a, b, d, epsilon, hsq6

if ic == 1
   prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
   sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat = tsamp)
   tvec = sol.t
   uout = sol.u
   push!(allsols, uout)
   push!(allts, tvec)

elseif ic == 2
   tspan = (0.,300.)
   tsamp = [300.]

   prob = ODEProblem(FitzHugh_Nagumo_ODE, u0, tspan, p)
   sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat = tsamp)

   # second perturbation line
   uout = sol.u
   u = uout[end]
   u[1,55:73,1:70] .= 0.999
   tspan = (300.,1000.)
   tsamp = [300., 400., 780., 781.]
   prob = ODEProblem(FitzHugh_Nagumo_ODE, u, tspan, p)
   sol = solve(prob, Tsit5(), reltol=1e-6, abstol=1e-6, maxiters = 1e7, saveat = tsamp)

   tvec = sol.t
   uout = sol.u

   tvec[2:4] = tvec[1:3]
   tvec[1] = 0.

   uout[2:4] = uout[1:3]
   uout[1] = u0  # initial condition
   push!(allsols, uout)
   push!(allts, tvec)

end
end # ic loop
end # cputime

fig = plt.figure(figsize=(figx,1.58*figy))

kplt = 0
for ic in (1,2)
   uout = allsols[ic]
   tvec = allts[ic]
   nts = length(tvec)

   for n = 1:nts
      ures = uout[n]
      kplt = kplt + 1
      subplot(2,4, kplt)
      im = plt.pcolormesh(ures[1,:,:],cmap = "turbo")   # "gnuplot" ) #"inferno")
      im.set_clim(-0.2, 1.)

      ax = gca()
      ax.set_xticks([])
      ax.set_yticks([])
      plt.title(L"t = "*string(floor(Int,tvec[n])), fontsize = 24)
      ax.set_aspect("equal")
   end

end # end of ic-loop

fig.tight_layout(pad=0.3)
add_identifiers!(fig)

wsave(plotsdir("11", "FitzHugh_Nagumo_2D_concentric_plane_spiral_waves"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
module Monoids
using ..SuiteSparseGraphBLAS: isGxB, isGrB, TypedMonoid, AbstractMonoid, splitconstant
using ..libgb
export Monoid

function _monoidnames(name)
    if isGxB(name) || isGrB(name)
        simple = splitconstant(name)[2]
    else
        simple = name
    end
    containername = Symbol(simple, "_MONOID_T")
    exportedname = Symbol(simple * "_MONOID")
    return containername, exportedname
end

function Monoid(name)
    containername, exportedname = _monoidnames(name)
    structquote = quote
        struct $containername <: AbstractMonoid
            typedops::Dict{DataType, TypedMonoid}
            name::String
            $containername() = new(Dict{DataType, TypedMonoid}(), $name)
        end
    end
    @eval($structquote)
    constquote = quote
        const $exportedname = $containername()
        export $exportedname
    end
    @eval($constquote)
    return getproperty(Monoids, exportedname)
end
end
const MonoidUnion = Union{AbstractMonoid, TypedMonoid}

#TODO: Rewrite
function _createmonoids()
    builtins = [
        "GrB_MIN",
        "GrB_MAX",
        "GrB_PLUS",
        "GrB_TIMES",
        "GxB_ANY",
        "GrB_LOR",
        "GrB_LAND",
        "GrB_LXOR",
        "GrB_LXNOR",
        "GxB_EQ",
        "GxB_BOR",
        "GxB_BAND",
        "GxB_BXOR",
        "GxB_BXNOR",
    ]
    for name ∈ builtins
        Monoids.Monoid(name)
    end
end

#This is adapted from the fork by cvdlab
#NOTE: Terminals do not work with built in BinaryOps.
function _addmonoid(op::AbstractMonoid, binop::BinaryUnion, id::T, terminal = nothing) where {T}
    if terminal === nothing
        terminal = C_NULL
    else
        typeof(terminal) == T || throw(ArgumentError("id and terminal must have the same type."))
    end
    if binop isa AbstractBinaryOp
        binop = binop[T]
    end
    monref = Ref{libgb.GrB_Monoid}()
    if T <: valid_union
        if terminal === C_NULL
            libgb.monoididnew[T](monref, binop, id)
        else
            libgb.monoidtermnew[T](monref, binop, id, terminal)
        end
    else
        libgb.monoidtermnew[Any](monref, binop, Ptr{Cvoid}(id), Ptr{Cvoid}(terminal))
    end
    op.typedops[T] = TypedMonoid{xtype(binop), ytype(binop), ztype(binop)}(monref[])
    return nothing
end

#Monoid Constructors
####################

function Monoids.Monoid(name::String, binop::BinaryUnion, id::T, terminal = nothing) where {T}
    if binop isa AbstractBinaryOp #If this is an AbstractBinaryOp we need to narrow down
        binop = binop[T]
    end
    m = Monoids.Monoid(name)
    _addmonoid(m, binop, id, terminal)
    return m
end
function Monoids.Monoid(name::String, binop::AbstractBinaryOp, id::AbstractVector, terminal = nothing)
    m = Monoids.Monoid(name)
    for i ∈ 1:length(id)
        binop2 = binop[typeof(id[i])]
        if terminal === nothing
            _addmonoid(m, binop2, id[i], nothing)
        else
            _addmonoid(m, binop2, id[i], terminal[i])
        end
    end
    return m
end

Monoids.Monoid(::Function) =
    error("You must construct monoids from an existing BinaryOp and an identity value.")

function _load(monoid::AbstractMonoid)
    booleans = ["GxB_ANY", "GrB_LOR", "GrB_LAND", "GrB_LXOR", "GrB_LXNOR", "GxB_EQ"]
    integers = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY"]
    unsignedintegers = [
        "GrB_MIN",
        "GrB_MAX",
        "GrB_PLUS",
        "GrB_TIMES",
        "GxB_ANY",
        "GxB_BOR",
        "GxB_BAND",
        "GxB_BXOR",
        "GxB_BXNOR",
    ]
    floats = ["GrB_MIN", "GrB_MAX", "GrB_PLUS", "GrB_TIMES", "GxB_ANY"]
    complexes = ["GxB_PLUS", "GxB_TIMES", "GxB_ANY"]
    name = monoid.name

    if name ∈ booleans
        constname = name * ((isGxB(name) ? "_BOOL_MONOID" : "_MONOID_BOOL"))
        monoid.typedops[Bool] = TypedMonoid(load_global(constname, libgb.GrB_Monoid))
    end

    if name ∈ integers
        monoid.typedops[Int8] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_INT8_MONOID" : "_MONOID_INT8"), libgb.GrB_Monoid))
        monoid.typedops[Int16] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_INT16_MONOID" : "_MONOID_INT16"), libgb.GrB_Monoid))
        monoid.typedops[Int32] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_INT32_MONOID" : "_MONOID_INT32"), libgb.GrB_Monoid))
        monoid.typedops[Int64] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_INT64_MONOID" : "_MONOID_INT64"), libgb.GrB_Monoid))
    end

    if name ∈ unsignedintegers
        monoid.typedops[UInt8] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_UINT8_MONOID" : "_MONOID_UINT8"), libgb.GrB_Monoid))
        monoid.typedops[UInt16] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_UINT16_MONOID" : "_MONOID_UINT16"), libgb.GrB_Monoid))
        monoid.typedops[UInt32] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_UINT32_MONOID" : "_MONOID_UINT32"), libgb.GrB_Monoid))
        monoid.typedops[UInt64] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_UINT64_MONOID" : "_MONOID_UINT64"), libgb.GrB_Monoid))
    end

    if name ∈ floats
        monoid.typedops[Float32] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_FP32_MONOID" : "_MONOID_FP32"), libgb.GrB_Monoid))
        monoid.typedops[Float64] =
            TypedMonoid(load_global(name * (isGxB(name) ? "_FP64_MONOID" : "_MONOID_FP64"), libgb.GrB_Monoid))
    end
    name = "GxB_" * name[5:end]
    if name ∈ complexes
        #Complex monoids are always GxB, so "_MONOID" is always at the end.
        monoid.typedops[ComplexF32] = TypedMonoid(load_global(name * "_FC32_MONOID", libgb.GrB_Monoid))
        monoid.typedops[ComplexF64] = TypedMonoid(load_global(name * "_FC64_MONOID", libgb.GrB_Monoid))
    end
end

ztype(::TypedMonoid{X, Y, Z}) where {X, Y, Z} = Z
xtype(::TypedMonoid{X, Y, Z}) where {X, Y, Z} = X
ytype(::TypedMonoid{X, Y, Z}) where {X, Y, Z} = Y


text:
"""
    fo7()
``fo7'' is an MINLPLib test instance with binary and continuous variables
(see http://www.minlplib.org/fo7.html). The best known primal feasible
objective value for this problem is 20.72982507. This example contains 14 nonlinear
terms, all of which are of the form 1/x.
"""

function post_fo7_linear_constraints!(m)

    @variable(m, objvar)
    b_Idx = Any[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]
    @variable(m, b[b_Idx], Bin)
    x_Idx = Any[44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115]
    @variable(m, x[x_Idx])
    set_lower_bound(x[86], 2.0)
    set_upper_bound(x[86], 8.0)
    set_lower_bound(x[87], 2.0)
    set_upper_bound(x[87], 8.0)
    set_lower_bound(x[88], 2.0)
    set_upper_bound(x[88], 8.0)
    set_lower_bound(x[89], 3.0)
    set_upper_bound(x[89], 8.5399)
    set_lower_bound(x[90], 1.5)
    set_upper_bound(x[90], 6.0)
    set_lower_bound(x[91], 1.5)
    set_upper_bound(x[91], 6.0)
    set_lower_bound(x[92], 1.5)
    set_upper_bound(x[92], 6.0)
    set_lower_bound(x[93], 8.54)
    set_upper_bound(x[93], 8.54)
    set_lower_bound(x[94], 2.0)
    set_upper_bound(x[94], 8.0)
    set_lower_bound(x[95], 2.0)
    set_upper_bound(x[95], 8.0)
    set_lower_bound(x[96], 2.0)
    set_upper_bound(x[96], 8.0)
    set_lower_bound(x[97], 4.2155)
    set_upper_bound(x[97], 12.0)
    set_lower_bound(x[98], 1.5)
    set_upper_bound(x[98], 6.0)
    set_lower_bound(x[99], 1.5)
    set_upper_bound(x[99], 6.0)
    set_lower_bound(x[100], 1.5)
    set_upper_bound(x[100], 6.0)
    set_lower_bound(x[101], 13.0)
    set_upper_bound(x[101], 13.0)

    # ----- Constraints ----- #
    @constraint(m, e1, objvar-x[44]-x[45]-x[56]-x[57]-x[66]-x[67]-x[74]-x[75]-x[80]-x[81]-x[84]-x[85] == 0.0)
    @constraint(m, e2, x[102]-x[103] <= 0.0)
    @constraint(m, e3, 0.5*x[86]-x[93]+x[102] <= 0.0)
    @constraint(m, e4, 0.5*x[86]-x[102] <= 0.0)
    @constraint(m, e5, 0.5*x[94]-x[101]+x[109] <= 0.0)
    @constraint(m, e6, 0.5*x[94]-x[109] <= 0.0)
    @constraint(m, e7, 0.5*x[87]-x[93]+x[103] <= 0.0)
    @constraint(m, e8, 0.5*x[87]-x[103] <= 0.0)
    @constraint(m, e9, 0.5*x[95]-x[101]+x[110] <= 0.0)
    @constraint(m, e10, 0.5*x[95]-x[110] <= 0.0)
    @constraint(m, e11, 0.5*x[88]-x[93]+x[104] <= 0.0)
    @constraint(m, e12, 0.5*x[88]-x[104] <= 0.0)
    @constraint(m, e13, 0.5*x[96]-x[101]+x[111] <= 0.0)
    @constraint(m, e14, 0.5*x[96]-x[111] <= 0.0)
    @constraint(m, e15, 0.5*x[89]-x[93]+x[105] <= 0.0)
    @constraint(m, e16, 0.5*x[89]-x[105] <= 0.0)
    @constraint(m, e17, 0.5*x[97]-x[101]+x[112] <= 0.0)
    @constraint(m, e18, 0.5*x[97]-x[112] <= 0.0)
    @constraint(m, e19, 0.5*x[90]-x[93]+x[106] <= 0.0)
    @constraint(m, e20, 0.5*x[90]-x[106] <= 0.0)
    @constraint(m, e21, 0.5*x[98]-x[101]+x[113] <= 0.0)
    @constraint(m, e22, 0.5*x[98]-x[113] <= 0.0)
    @constraint(m, e23, 0.5*x[91]-x[93]+x[107] <= 0.0)
    @constraint(m, e24, 0.5*x[91]-x[107] <= 0.0)
    @constraint(m, e25, 0.5*x[99]-x[101]+x[114] <= 0.0)
    @constraint(m, e26, 0.5*x[99]-x[114] <= 0.0)
    @constraint(m, e27, 0.5*x[92]-x[93]+x[108] <= 0.0)
    @constraint(m, e28, 0.5*x[92]-x[108] <= 0.0)
    @constraint(m, e29, 0.5*x[100]-x[101]+x[115] <= 0.0)
    @constraint(m, e30, 0.5*x[100]-x[115] <= 0.0)
    @constraint(m, e31, -x[44]+x[102]-x[103] <= 0.0)
    @constraint(m, e32, -x[44]-x[102]+x[103] <= 0.0)
    @constraint(m, e33, -x[45]+x[109]-x[110] <= 0.0)
    @constraint(m, e34, -x[45]-x[109]+x[110] <= 0.0)
    @constraint(m, e35, -8.54*b[1]-8.54*b[2]+0.5*x[86]+0.5*x[87]-x[102]+x[103] <= 0.0)
    @constraint(m, e36, -8.54*b[1]+8.54*b[2]+0.5*x[86]+0.5*x[87]+x[102]-x[103] <= 8.54)
    @constraint(m, e37, 13*b[1]-13*b[2]+0.5*x[94]+0.5*x[95]-x[109]+x[110] <= 13.0)
    @constraint(m, e38, 13*b[1]+13*b[2]+0.5*x[94]+0.5*x[95]+x[109]-x[110] <= 26.0)
    @constraint(m, e39, -x[46]+x[102]-x[104] <= 0.0)
    @constraint(m, e40, -x[46]-x[102]+x[104] <= 0.0)
    @constraint(m, e41, -x[47]+x[109]-x[111] <= 0.0)
    @constraint(m, e42, -x[47]-x[109]+x[111] <= 0.0)
    @constraint(m, e43, -8.54*b[3]-8.54*b[4]+0.5*x[86]+0.5*x[88]-x[102]+x[104] <= 0.0)
    @constraint(m, e44, -8.54*b[3]+8.54*b[4]+0.5*x[86]+0.5*x[88]+x[102]-x[104] <= 8.54)
    @constraint(m, e45, 13*b[3]-13*b[4]+0.5*x[94]+0.5*x[96]-x[109]+x[111] <= 13.0)
    @constraint(m, e46, 13*b[3]+13*b[4]+0.5*x[94]+0.5*x[96]+x[109]-x[111] <= 26.0)
    @constraint(m, e47, -x[48]+x[102]-x[105] <= 0.0)
    @constraint(m, e48, -x[48]-x[102]+x[105] <= 0.0)
    @constraint(m, e49, -x[49]+x[109]-x[112] <= 0.0)
    @constraint(m, e50, -x[49]-x[109]+x[112] <= 0.0)
    @constraint(m, e51, -8.54*b[5]-8.54*b[6]+0.5*x[86]+0.5*x[89]-x[102]+x[105] <= 0.0)
    @constraint(m, e52, -8.54*b[5]+8.54*b[6]+0.5*x[86]+0.5*x[89]+x[102]-x[105] <= 8.54)
    @constraint(m, e53, 13*b[5]-13*b[6]+0.5*x[94]+0.5*x[97]-x[109]+x[112] <= 13.0)
    @constraint(m, e54, 13*b[5]+13*b[6]+0.5*x[94]+0.5*x[97]+x[109]-x[112] <= 26.0)
    @constraint(m, e55, -x[50]+x[102]-x[106] <= 0.0)
    @constraint(m, e56, -x[50]-x[102]+x[106] <= 0.0)
    @constraint(m, e57, -x[51]+x[109]-x[113] <= 0.0)
    @constraint(m, e58, -x[51]-x[109]+x[113] <= 0.0)
    @constraint(m, e59, -8.54*b[7]-8.54*b[8]+0.5*x[86]+0.5*x[90]-x[102]+x[106] <= 0.0)
    @constraint(m, e60, -8.54*b[7]+8.54*b[8]+0.5*x[86]+0.5*x[90]+x[102]-x[106] <= 8.54)
    @constraint(m, e61, 13*b[7]-13*b[8]+0.5*x[94]+0.5*x[98]-x[109]+x[113] <= 13.0)
    @constraint(m, e62, 13*b[7]+13*b[8]+0.5*x[94]+0.5*x[98]+x[109]-x[113] <= 26.0)
    @constraint(m, e63, -x[52]+x[102]-x[107] <= 0.0)
    @constraint(m, e64, -x[52]-x[102]+x[107] <= 0.0)
    @constraint(m, e65, -x[53]+x[109]-x[114] <= 0.0)
    @constraint(m, e66, -x[53]-x[109]+x[114] <= 0.0)
    @constraint(m, e67, -8.54*b[9]-8.54*b[10]+0.5*x[86]+0.5*x[91]-x[102]+x[107] <= 0.0)
    @constraint(m, e68, -8.54*b[9]+8.54*b[10]+0.5*x[86]+0.5*x[91]+x[102]-x[107] <= 8.54)
    @constraint(m, e69, 13*b[9]-13*b[10]+0.5*x[94]+0.5*x[99]-x[109]+x[114] <= 13.0)
    @constraint(m, e70, 13*b[9]+13*b[10]+0.5*x[94]+0.5*x[99]+x[109]-x[114] <= 26.0)
    @constraint(m, e71, -x[54]+x[102]-x[108] <= 0.0)
    @constraint(m, e72, -x[54]-x[102]+x[108] <= 0.0)
    @constraint(m, e73, -x[55]+x[109]-x[115] <= 0.0)
    @constraint(m, e74, -x[55]-x[109]+x[115] <= 0.0)
    @constraint(m, e75, -8.54*b[11]-8.54*b[12]+0.5*x[86]+0.5*x[92]-x[102]+x[108] <= 0.0)
    @constraint(m, e76, -8.54*b[11]+8.54*b[12]+0.5*x[86]+0.5*x[92]+x[102]-x[108] <= 8.54)
    @constraint(m, e77, 13*b[11]-13*b[12]+0.5*x[94]+0.5*x[100]-x[109]+x[115] <= 13.0)
    @constraint(m, e78, 13*b[11]+13*b[12]+0.5*x[94]+0.5*x[100]+x[109]-x[115] <= 26.0)
    @constraint(m, e79, -x[56]+x[103]-x[104] <= 0.0)
    @constraint(m, e80, -x[56]-x[103]+x[104] <= 0.0)
    @constraint(m, e81, -x[57]+x[110]-x[111] <= 0.0)
    @constraint(m, e82, -x[57]-x[110]+x[111] <= 0.0)
    @constraint(m, e83, -8.54*b[13]-8.54*b[14]+0.5*x[87]+0.5*x[88]-x[103]+x[104] <= 0.0)
    @constraint(m, e84, -8.54*b[13]+8.54*b[14]+0.5*x[87]+0.5*x[88]+x[103]-x[104] <= 8.54)
    @constraint(m, e85, 13*b[13]-13*b[14]+0.5*x[95]+0.5*x[96]-x[110]+x[111] <= 13.0)
    @constraint(m, e86, 13*b[13]+13*b[14]+0.5*x[95]+0.5*x[96]+x[110]-x[111] <= 26.0)
    @constraint(m, e87, -x[58]+x[103]-x[105] <= 0.0)
    @constraint(m, e88, -x[58]-x[103]+x[105] <= 0.0)
    @constraint(m, e89, -x[59]+x[110]-x[112] <= 0.0)
    @constraint(m, e90, -x[59]-x[110]+x[112] <= 0.0)
    @constraint(m, e91, -8.54*b[15]-8.54*b[16]+0.5*x[87]+0.5*x[89]-x[103]+x[105] <= 0.0)
    @constraint(m, e92, -8.54*b[15]+8.54*b[16]+0.5*x[87]+0.5*x[89]+x[103]-x[105] <= 8.54)
    @constraint(m, e93, 13*b[15]-13*b[16]+0.5*x[95]+0.5*x[97]-x[110]+x[112] <= 13.0)
    @constraint(m, e94, 13*b[15]+13*b[16]+0.5*x[95]+0.5*x[97]+x[110]-x[112] <= 26.0)
    @constraint(m, e95, -x[60]+x[103]-x[106] <= 0.0)
    @constraint(m, e96, -x[60]-x[103]+x[106] <= 0.0)
    @constraint(m, e97, -x[61]+x[110]-x[113] <= 0.0)
    @constraint(m, e98, -x[61]-x[110]+x[113] <= 0.0)
    @constraint(m, e99, -8.54*b[17]-8.54*b[18]+0.5*x[87]+0.5*x[90]-x[103]+x[106] <= 0.0)
    @constraint(m, e100, -8.54*b[17]+8.54*b[18]+0.5*x[87]+0.5*x[90]+x[103]-x[106] <= 8.54)
    @constraint(m, e101, 13*b[17]-13*b[18]+0.5*x[95]+0.5*x[98]-x[110]+x[113] <= 13.0)
    @constraint(m, e102, 13*b[17]+13*b[18]+0.5*x[95]+0.5*x[98]+x[110]-x[113] <= 26.0)
    @constraint(m, e103, -x[62]+x[103]-x[107] <= 0.0)
    @constraint(m, e104, -x[62]-x[103]+x[107] <= 0.0)
    @constraint(m, e105, -x[63]+x[110]-x[114] <= 0.0)
    @constraint(m, e106, -x[63]-x[110]+x[114] <= 0.0)
    @constraint(m, e107, -8.54*b[19]-8.54*b[20]+0.5*x[87]+0.5*x[91]-x[103]+x[107] <= 0.0)
    @constraint(m, e108, -8.54*b[19]+8.54*b[20]+0.5*x[87]+0.5*x[91]+x[103]-x[107] <= 8.54)
    @constraint(m, e109, 13*b[19]-13*b[20]+0.5*x[95]+0.5*x[99]-x[110]+x[114] <= 13.0)
    @constraint(m, e110, 13*b[19]+13*b[20]+0.5*x[95]+0.5*x[99]+x[110]-x[114] <= 26.0)
    @constraint(m, e111, -x[64]+x[103]-x[108] <= 0.0)
    @constraint(m, e112, -x[64]-x[103]+x[108] <= 0.0)
    @constraint(m, e113, -x[65]+x[110]-x[115] <= 0.0)
    @constraint(m, e114, -x[65]-x[110]+x[115] <= 0.0)
    @constraint(m, e115, -8.54*b[21]-8.54*b[22]+0.5*x[87]+0.5*x[92]-x[103]+x[108] <= 0.0)
    @constraint(m, e116, -8.54*b[21]+8.54*b[22]+0.5*x[87]+0.5*x[92]+x[103]-x[108] <= 8.54)
    @constraint(m, e117, 13*b[21]-13*b[22]+0.5*x[95]+0.5*x[100]-x[110]+x[115] <= 13.0)
    @constraint(m, e118, 13*b[21]+13*b[22]+0.5*x[95]+0.5*x[100]+x[110]-x[115] <= 26.0)
    @constraint(m, e119, -x[66]+x[104]-x[105] <= 0.0)
    @constraint(m, e120, -x[66]-x[104]+x[105] <= 0.0)
    @constraint(m, e121, -x[67]+x[111]-x[112] <= 0.0)
    @constraint(m, e122, -x[67]-x[111]+x[112] <= 0.0)
    @constraint(m, e123, -8.54*b[23]-8.54*b[24]+0.5*x[88]+0.5*x[89]-x[104]+x[105] <= 0.0)
    @constraint(m, e124, -8.54*b[23]+8.54*b[24]+0.5*x[88]+0.5*x[89]+x[104]-x[105] <= 8.54)
    @constraint(m, e125, 13*b[23]-13*b[24]+0.5*x[96]+0.5*x[97]-x[111]+x[112] <= 13.0)
    @constraint(m, e126, 13*b[23]+13*b[24]+0.5*x[96]+0.5*x[97]+x[111]-x[112] <= 26.0)
    @constraint(m, e127, -x[68]+x[104]-x[106] <= 0.0)
    @constraint(m, e128, -x[68]-x[104]+x[106] <= 0.0)
    @constraint(m, e129, -x[69]+x[111]-x[113] <= 0.0)
    @constraint(m, e130, -x[69]-x[111]+x[113] <= 0.0)
    @constraint(m, e131, -8.54*b[25]-8.54*b[26]+0.5*x[88]+0.5*x[90]-x[104]+x[106] <= 0.0)
    @constraint(m, e132, -8.54*b[25]+8.54*b[26]+0.5*x[88]+0.5*x[90]+x[104]-x[106] <= 8.54)
    @constraint(m, e133, 13*b[25]-13*b[26]+0.5*x[96]+0.5*x[98]-x[111]+x[113] <= 13.0)
    @constraint(m, e134, 13*b[25]+13*b[26]+0.5*x[96]+0.5*x[98]+x[111]-x[113] <= 26.0)
    @constraint(m, e135, -x[70]+x[104]-x[107] <= 0.0)
    @constraint(m, e136, -x[70]-x[104]+x[107] <= 0.0)
    @constraint(m, e137, -x[71]+x[111]-x[114] <= 0.0)
    @constraint(m, e138, -x[71]-x[111]+x[114] <= 0.0)
    @constraint(m, e139, -8.54*b[27]-8.54*b[28]+0.5*x[88]+0.5*x[91]-x[104]+x[107] <= 0.0)
    @constraint(m, e140, -8.54*b[27]+8.54*b[28]+0.5*x[88]+0.5*x[91]+x[104]-x[107] <= 8.54)
    @constraint(m, e141, 13*b[27]-13*b[28]+0.5*x[96]+0.5*x[99]-x[111]+x[114] <= 13.0)
    @constraint(m, e142, 13*b[27]+13*b[28]+0.5*x[96]+0.5*x[99]+x[111]-x[114] <= 26.0)
    @constraint(m, e143, -x[72]+x[104]-x[108] <= 0.0)
    @constraint(m, e144, -x[72]-x[104]+x[108] <= 0.0)
    @constraint(m, e145, -x[73]+x[111]-x[115] <= 0.0)
    @constraint(m, e146, -x[73]-x[111]+x[115] <= 0.0)
    @constraint(m, e147, -8.54*b[29]-8.54*b[30]+0.5*x[88]+0.5*x[92]-x[104]+x[108] <= 0.0)
    @constraint(m, e148, -8.54*b[29]+8.54*b[30]+0.5*x[88]+0.5*x[92]+x[104]-x[108] <= 8.54)
    @constraint(m, e149, 13*b[29]-13*b[30]+0.5*x[96]+0.5*x[100]-x[111]+x[115] <= 13.0)
    @constraint(m, e150, 13*b[29]+13*b[30]+0.5*x[96]+0.5*x[100]+x[111]-x[115] <= 26.0)
    @constraint(m, e151, -x[74]+x[105]-x[106] <= 0.0)
    @constraint(m, e152, -x[74]-x[105]+x[106] <= 0.0)
    @constraint(m, e153, -x[75]+x[112]-x[113] <= 0.0)
    @constraint(m, e154, -x[75]-x[112]+x[113] <= 0.0)
    @constraint(m, e155, -8.54*b[31]-8.54*b[32]+0.5*x[89]+0.5*x[90]-x[105]+x[106] <= 0.0)
    @constraint(m, e156, -8.54*b[31]+8.54*b[32]+0.5*x[89]+0.5*x[90]+x[105]-x[106] <= 8.54)
    @constraint(m, e157, 13*b[31]-13*b[32]+0.5*x[97]+0.5*x[98]-x[112]+x[113] <= 13.0)
    @constraint(m, e158, 13*b[31]+13*b[32]+0.5*x[97]+0.5*x[98]+x[112]-x[113] <= 26.0)
    @constraint(m, e159, -x[76]+x[105]-x[107] <= 0.0)
    @constraint(m, e160, -x[76]-x[105]+x[107] <= 0.0)
    @constraint(m, e161, -x[77]+x[112]-x[114] <= 0.0)
    @constraint(m, e162, -x[77]-x[112]+x[114] <= 0.0)
    @constraint(m, e163, -8.54*b[33]-8.54*b[34]+0.5*x[89]+0.5*x[91]-x[105]+x[107] <= 0.0)
    @constraint(m, e164, -8.54*b[33]+8.54*b[34]+0.5*x[89]+0.5*x[91]+x[105]-x[107] <= 8.54)
    @constraint(m, e165, 13*b[33]-13*b[34]+0.5*x[97]+0.5*x[99]-x[112]+x[114] <= 13.0)
    @constraint(m, e166, 13*b[33]+13*b[34]+0.5*x[97]+0.5*x[99]+x[112]-x[114] <= 26.0)
    @constraint(m, e167, -x[78]+x[105]-x[108] <= 0.0)
    @constraint(m, e168, -x[78]-x[105]+x[108] <= 0.0)
    @constraint(m, e169, -x[79]+x[112]-x[115] <= 0.0)
    @constraint(m, e170, -x[79]-x[112]+x[115] <= 0.0)
    @constraint(m, e171, -8.54*b[35]-8.54*b[36]+0.5*x[89]+0.5*x[92]-x[105]+x[108] <= 0.0)
    @constraint(m, e172, -8.54*b[35]+8.54*b[36]+0.5*x[89]+0.5*x[92]+x[105]-x[108] <= 8.54)
    @constraint(m, e173, 13*b[35]-13*b[36]+0.5*x[97]+0.5*x[100]-x[112]+x[115] <= 13.0)
    @constraint(m, e174, 13*b[35]+13*b[36]+0.5*x[97]+0.5*x[100]+x[112]-x[115] <= 26.0)
    @constraint(m, e175, -x[80]+x[106]-x[107] <= 0.0)
    @constraint(m, e176, -x[80]-x[106]+x[107] <= 0.0)
    @constraint(m, e177, -x[81]+x[113]-x[114] <= 0.0)
    @constraint(m, e178, -x[81]-x[113]+x[114] <= 0.0)
    @constraint(m, e179, -8.54*b[37]-8.54*b[38]+0.5*x[90]+0.5*x[91]-x[106]+x[107] <= 0.0)
    @constraint(m, e180, -8.54*b[37]+8.54*b[38]+0.5*x[90]+0.5*x[91]+x[106]-x[107] <= 8.54)
    @constraint(m, e181, 13*b[37]-13*b[38]+0.5*x[98]+0.5*x[99]-x[113]+x[114] <= 13.0)
    @constraint(m, e182, 13*b[37]+13*b[38]+0.5*x[98]+0.5*x[99]+x[113]-x[114] <= 26.0)
    @constraint(m, e183, -x[82]+x[106]-x[108] <= 0.0)
    @constraint(m, e184, -x[82]-x[106]+x[108] <= 0.0)
    @constraint(m, e185, -x[83]+x[113]-x[115] <= 0.0)
    @constraint(m, e186, -x[83]-x[113]+x[115] <= 0.0)
    @constraint(m, e187, -8.54*b[39]-8.54*b[40]+0.5*x[90]+0.5*x[92]-x[106]+x[108] <= 0.0)
    @constraint(m, e188, -8.54*b[39]+8.54*b[40]+0.5*x[90]+0.5*x[92]+x[106]-x[108] <= 8.54)
    @constraint(m, e189, 13*b[39]-13*b[40]+0.5*x[98]+0.5*x[100]-x[113]+x[115] <= 13.0)
    @constraint(m, e190, 13*b[39]+13*b[40]+0.5*x[98]+0.5*x[100]+x[113]-x[115] <= 26.0)
    @constraint(m, e191, -x[84]+x[107]-x[108] <= 0.0)
    @constraint(m, e192, -x[84]-x[107]+x[108] <= 0.0)
    @constraint(m, e193, -x[85]+x[114]-x[115] <= 0.0)
    @constraint(m, e194, -x[85]-x[114]+x[115] <= 0.0)
    @constraint(m, e195, -8.54*b[41]-8.54*b[42]+0.5*x[91]+0.5*x[92]-x[107]+x[108] <= 0.0)
    @constraint(m, e196, -8.54*b[41]+8.54*b[42]+0.5*x[91]+0.5*x[92]+x[107]-x[108] <= 8.54)
    @constraint(m, e197, 13*b[41]-13*b[42]+0.5*x[99]+0.5*x[100]-x[114]+x[115] <= 13.0)
    @constraint(m, e198, 13*b[41]+13*b[42]+0.5*x[99]+0.5*x[100]+x[114]-x[115] <= 26.0)
    println("creating lifted variables and constraints")
    @variable(m, y[1:14])
    @constraint(m, e199, 16*y[1]-x[94] <= 0.0)
    @constraint(m, e200, 16*y[2]-x[86] <= 0.0)
    @constraint(m, e201, 16*y[3]-x[95] <= 0.0)
    @constraint(m, e202, 16*y[4]-x[87] <= 0.0)
    @constraint(m, e203, 16*y[5]-x[96] <= 0.0)
    @constraint(m, e204, 16*y[6]-x[88] <= 0.0)
    @constraint(m, e205, 36*y[7]-x[97] <= 0.0)
    @constraint(m, e206, 36*y[8]-x[89] <= 0.0)
    @constraint(m, e207, 9*y[9]-x[98] <= 0.0)
    @constraint(m, e208, 9*y[10]-x[90] <= 0.0)
    @constraint(m, e209, 9*y[11]-x[99] <= 0.0)
    @constraint(m, e210, 9*y[12]-x[91] <= 0.0)
    @constraint(m, e211, 9*y[13]-x[100] <= 0.0)
    @constraint(m, e212, 9*y[14]-x[92] <= 0.0)
    y_indexes = collect(1:14)
    x_indexes = [86, 94, 87, 95, 88, 96, 89, 97, 90, 98, 91, 99, 92, 100]
    partitions = []
    for index in x_indexes
        push!(partitions, [lower_bound(x[index]), upper_bound(x[index])])
    end

    # @constraint(m, e199, 16/x[86]-x[94] <= 0.0)
    # @constraint(m, e200, 16/x[94]-x[86] <= 0.0)
    # @constraint(m, e201, 16/x[87]-x[95] <= 0.0)
    # @constraint(m, e202, 16/x[95]-x[87] <= 0.0)
    # @constraint(m, e203, 16/x[88]-x[96] <= 0.0)
    # @constraint(m, e204, 16/x[96]-x[88] <= 0.0)
    # @constraint(m, e205, 36/x[89]-x[97] <= 0.0)
    # @constraint(m, e206, 36/x[97]-x[89] <= 0.0)
    # @constraint(m, e207, 9/x[90]-x[98] <= 0.0)
    # @constraint(m, e208, 9/x[98]-x[90] <= 0.0)
    # @constraint(m, e209, 9/x[91]-x[99] <= 0.0)
    # @constraint(m, e210, 9/x[99]-x[91] <= 0.0)
    # @constraint(m, e211, 9/x[92]-x[100] <= 0.0)
    # @constraint(m, e212, 9/x[100]-x[92] <= 0.0)

    # ----- Objective ----- #
    @objective(m, Min, objvar)
    return y_indexes, x_indexes, partitions
end

function fo7()
    best_known_objective = 20.72982507
    cplex_optimizer =
        JuMP.optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0)
    error_tolerances = [NaN64, 1e-1, 1e-2]
    # base partition holds the inflection points and the end-points for the trigonometric functions
    milp_statistics = zeros(length(error_tolerances), 5)
    lp_statistics = zeros(length(error_tolerances), 5)
    func = x -> 1/x
    # solve MILP for various error tolerances
    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        milp = Model(cplex_optimizer)
        y_indexes, x_indexes, base_partition = post_fo7_linear_constraints!(milp)
        # construct MILP relaxations
        for j = 1:length(y_indexes)
            var, x_index, y_index = add_milp_relaxation(
                func,
                base_partition[j],
                error_tolerance,
                milp,
                name = string(j),
            )
            @constraint(milp, var[x_index] == milp[:x][x_indexes[j]])
            @constraint(milp, var[y_index] == milp[:y][y_indexes[j]])
        end
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(milp)
        relaxation_objective = objective_value(milp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        milp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    for i = 1:length(error_tolerances)
        error_tolerance = error_tolerances[i]
        lp = Model(cplex_optimizer)
        y_indexes, x_indexes, base_partition = post_fo7_linear_constraints!(lp)
        # construct MILP relaxations
        for j = 1:length(y_indexes)
            var, x_index, y_index = add_lp_relaxation(
                func,
                base_partition[j],
                error_tolerance,
                lp,
                name = string(j),
            )
            @constraint(lp, var[x_index] == lp[:x][x_indexes[j]])
            @constraint(lp, var[y_index] == lp[:y][y_indexes[j]])
        end
        _, solve_time, solve_bytes_alloc, sec_in_gc = @timed optimize!(lp)
        relaxation_objective = objective_value(lp)
        @assert relaxation_objective ≤ best_known_objective
        relative_gap =
            abs(best_known_objective - relaxation_objective) / abs(relaxation_objective)
        lp_statistics[i, :] = [
            error_tolerance,
            best_known_objective,
            relaxation_objective,
            relative_gap,
            round(solve_time, digits = 2),
        ]
    end

    open("csv/fo7_milp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, milp_statistics, ",")
    end

    open("csv/fo7_lp.csv", "w") do io
        writedlm(io, ["error_tol" "g_opt" "relaxation" "relative_gap" "time"], ',')
        writedlm(io, lp_statistics, ",")
    end

end

max_stars_repo_name: sujeevraja/univariate-experiments

max_issues_repo_name: sujeevraja/univariate-experiments

max_forks_repo_name: sujeevraja/univariate-experiments


text:
abstract type AbsLinearOperator{T,N} #=<: AbstractArray{T,N} end=# end
const AbsLinearOp{T} = AbsLinearOperator{T,2}
const AbsLinearSuperOperator{T} = AbsLinearOperator{T,4}

abstract type AbsOpConnection end

Base.eltype(c::AbsOpConnection) = eltype(typeof(c))

# Random sampler for an abstract OpConnection (assuming it implements getindex and length)
function Random.rand(rng::AbstractRNG, sampl::Random.SamplerTrivial{T}) where {T<:NeuralQuantum.AbsOpConnection}
    conns = sampl[]
    N = length(conns)
    i = rand(rng, 1:N)
    @inbounds res = conns[i]
    return res
end

"""
    row_valdiff(op::AbsLinearOperator, v) -> OpConnection

Returns all non-zero elements in the row represented by the state `v` of
operator `op`. The result is an `OpConnection` type.
"""
row_valdiff(op::AbsLinearOperator, v) = row_valdiff!(conn_type(op)(), op, v, init=false)

"""
    row_valdiff!(opconn::OpConnection, op::AbsLinearOperator, v; init=true)

Returns all non-zero elements in the row represented by the state `v` of
operator `op` by mutating in-place the `OpConnection` object `opconn`
"""
function row_valdiff!(opconn, op, v; init=true)
    init && resize!(opconn, 0)
    _row_valdiff!(opconn, op, v)
end

function row_valstate(op::AbsLinearOperator, v)
    conns = row_valdiff(op, v)

    vp    = state_similar(v, length(conns))
    state_copy!(vp, v)

    T     = eltype(conns).parameters[1]
    mels  = zeros(T, length(conns))

    for (i,(mel, cngs)) = enumerate(conns)
        mels[i] = mel
        apply!(unsafe_get_batch(vp,i), cngs)
    end
    return mels, vp
end

row_valdiff_index(op::AbsLinearOperator, v) = row_valdiff_index!(OpConnectionIndex(op), op, v)
function row_valdiff_index! end

# standard functions
Base.:*(op::AbsLinearOperator, α::Number) = α*op
Base.:/(op::AbsLinearOperator, α::Number) = inv(α)*op

conn_type(op::AbsLinearOperator) = conn_type(typeof(op))
## AbstractArray Interface


text:
"""
This module provides utilities to build neural networks with Flux,
along with a library of standard architectures.
"""
module GraphNeuralNetworkLib

export SimpleGNN, SimpleGraphNetHP

using ..AlphaZero

using CUDA
using Base: @kwdef

import Flux

CUDA.allowscalar(false)
array_on_gpu(::Array) = false
array_on_gpu(::CuArray) = true
array_on_gpu(arr) = error("Usupported array type: ", typeof(arr))

using Flux: relu, softmax, flatten
using Flux: Chain, Dense, Conv, BatchNorm, SkipConnection
import Zygote

# include("graph_network.jl")

#####
##### Flux Networks
#####

"""
    FluxNetwork <: AbstractNetwork

Abstract type for neural networks implemented using the _Flux_ framework.

The `regularized_params_` function must be overrided for all layers containing
parameters that are subject to regularization.

Provided that the above holds, `FluxNetwork` implements the full
network interface with the following exceptions:
[`Network.HyperParams`](@ref), [`Network.hyperparams`](@ref),
[`Network.forward`](@ref) and [`Network.on_gpu`](@ref).
"""
abstract type FluxGNN <: GraphNetwork end

function Base.copy(nn::Net) where Net <: FluxGNN
  #new = Net(Network.hyperparams(nn))
  #Flux.loadparams!(new, Flux.params(nn))
  #return new
  return Base.deepcopy(nn)
end

Network.to_cpu(nn::FluxGNN) = Flux.cpu(nn)

function Network.to_gpu(nn::FluxGNN)
  CUDA.allowscalar(false)
  return Flux.gpu(nn)
end

function Network.set_test_mode!(nn::FluxGNN, mode)
  Flux.testmode!(nn, mode)
end

Network.convert_input(nn::FluxGNN, x) =
  Network.on_gpu(nn) ? Flux.gpu(x) : x

Network.convert_output(nn::FluxGNN, x) = Flux.cpu(x)

Network.params(nn::FluxGNN) = Flux.params(nn)

# This should be included in Flux
function lossgrads(f, args...)
  val, back = Zygote.pullback(f, args...)
  grad = back(Zygote.sensitivity(val))
  return val, grad
end

function Network.train!(callback, nn::FluxGNN, opt::Adam, loss, data, n)
  optimiser = Flux.ADAM(opt.lr)
  params = Flux.params(nn)
  for (i, d) in enumerate(data)
    l, grads = lossgrads(params) do
      loss(d...)
    end
    Flux.update!(optimiser, params, grads)
    callback(i, l)
  end
end

function Network.train!(
    callback, nn::FluxGNN, opt::CyclicNesterov, loss, data, n)
  lr = CyclicSchedule(
    opt.lr_base,
    opt.lr_high,
    opt.lr_low, n=n)
  momentum = CyclicSchedule(
    opt.momentum_high,
    opt.momentum_low,
    opt.momentum_high, n=n)
  optimiser = Flux.Nesterov(opt.lr_low, opt.momentum_high)
  params = Flux.params(nn)
  for (i, d) in enumerate(data)
    l, grads = lossgrads(params) do
      loss(d...)
    end
    Flux.update!(optimiser, params, grads)
    optimiser.eta = lr[i]
    optimiser.rho = momentum[i]
    callback(i, l)
  end
end

regularized_params_(l) = []
regularized_params_(l::Flux.Dense) = [l.weight]
regularized_params_(l::Flux.Conv) = [l.weight]

function Network.regularized_params(net::GraphNetwork)
  return (w for l in Flux.modules(net) for w in regularized_params_(l))
end

function Network.gc(::GraphNetwork)
  GC.gc(true)
  # CUDA.reclaim()
end

# Flux.@functor does not work with abstract types
function Flux.functor(nn::Net) where Net <: FluxGNN
  children = (nn.common, nn.vhead, nn.phead, )
  constructor = cs -> Net(nn.gspec, nn.hyper, cs...)
  return (children, constructor)
end

Network.hyperparams(nn::FluxGNN) = nn.hyper

Network.game_spec(nn::FluxGNN) = nn.gspec

Network.on_gpu(nn::FluxGNN) = array_on_gpu(nn.common[end].bias)

#####
##### Include networks library
#####

include("architectures/demographnet.jl")

end


text:
# This function relates to the ansiotropic case

function generate_multi_index(d::S,mu::Array{S,1}) where {S<:Integer}

  nt = num_terms(mu,d)
  multi_index = Array{S,2}(undef,nt,d)
  multi_index[1,:] = ones(S,1,d)

  max_mu = maximum(mu)

  w = Tuple(repeat([max_mu+1],inner = d))
  j = 1
  @inbounds for i = 2:(max_mu+1)^d
    candidate_index = Tuple(CartesianIndices(w)[i])
    if sum(candidate_index) <= d+max_mu && sum(candidate_index .<= mu.+1) == d
      #if j < nt
        j += 1
        multi_index[j,:] = [candidate_index...]
      #else # handles case where nt is under-estimated (this branch is not used?)
      #  multi_index = [multi_index; collect(candidate_index)']
      #end
    end
  end
  if j < nt # handles case where nt is over-estimated
    multi_index = multi_index[1:j,:]
  end

  return multi_index

end

# The function below relates to the isotropic case

function generate_multi_index(d::S,mu::S) where {S<:Integer}

  nt = num_terms(mu,d)
  multi_index      = Array{S,2}(undef,nt,d)
  multi_index[1,:] = ones(S,1,d)

  w = Tuple(repeat([mu+1],inner = d))
  j = 1
  @inbounds for i = 2:(mu+1)^d
    candidate_index = Tuple(CartesianIndices(w)[i])
    if sum(candidate_index) <= d+mu
      j += 1
      multi_index[j,:] = [candidate_index...]
    end
  end

  return multi_index

end

# Following function computes the number of terms in the multi-index for the
# isotropic case (it also computes the number of terms in a complete
# polynominal based on the order and the number of dimensions.

function num_terms(order::S,d::S) where {S <: Integer}

  if d == 1
    return order+1
  else
    return div(num_terms(order,d-1)*(order+d),d)
  end

end

# The following function is a poor approximation to the number of terms in
# the multi-index for the ansiotropic case.

function num_terms(order::Array{S,1},d::S) where {S <: Integer}

  max_mu = maximum(order)
  nt = num_terms(max_mu,d) # Deliberate over-estimate of the number of terms

  return nt

end


text:
using GeometricBase
using GeometricBase.Config
using GeometricIntegrators.Integrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Tableaus
using GeometricProblems.HarmonicOscillator
using Test

using GeometricEquations: _get_v, _get_f, _get_v̄, _get_f̄
using GeometricIntegrators.Integrators: _euler_extrapolation_ode!, _hermite_extrapolation!,
                                        _midpoint_extrapolation_ode!, _midpoint_extrapolation_iode!, _midpoint_extrapolation_pode!

const Δt = 0.1
const t₀ = 0.0
const t₁ = t₀ + Δt
const t₂ = t₁ + Δt

ode  = harmonic_oscillator_ode()
pode = harmonic_oscillator_pode()
iode = harmonic_oscillator_iode()

# Create ODE Solution Arrays

x₀ = ode.q₀[1]
x₁ = zero(x₀)
x₂ = zero(x₀)
xᵢ = zero(x₀)

ẋₚ = zero(x₀)
ẋ₀ = zero(x₀)
ẋ₁ = zero(x₀)
ẋ₂ = zero(x₀)
ẋₙ = zero(x₀)
ẋᵢ = zero(x₀)

# Compute Reference Solution for ODEs

ref_prev = integrate(ode, TableauGauss(8), -Δt, 1)
ref_next = integrate(ode, TableauGauss(8), +Δt, 1)

tₚ = ref_prev.t[end]
xₚ = ref_prev.q[end]

tₙ = ref_next.t[end]
xₙ = ref_next.q[end]

tᵢ = tₙ

ode.v(tₚ, xₚ, ẋₚ, ode.parameters)
ode.v(t₀, x₀, ẋ₀, ode.parameters)
ode.v(tₙ, xₙ, ẋₙ, ode.parameters)

# Hermite Extrapolation

evaluate!(HermiteExtrapolation(tₚ, t₀), xₚ, x₀, ẋₚ, ẋ₀, tᵢ, xᵢ, ẋᵢ)

# println(xᵢ, xₙ, xᵢ .- xₙ)
# println(ẋᵢ, ẋₙ, ẋᵢ .- ẋₙ)

@test xᵢ ≈ xₙ atol=1E-5
@test ẋᵢ ≈ ẋₙ atol=1E-4

@test _hermite_extrapolation!(tₚ, t₀, xₚ, x₀, ẋₚ, ẋ₀, t₁, x₁) == xᵢ
@test _hermite_extrapolation!(tₚ, t₀, xₚ, x₀, ẋₚ, ẋ₀, t₁, x₁, ẋ₁) == (xᵢ, ẋᵢ)
@test _hermite_extrapolation!(_get_v(ode), tₚ, t₀, xₚ, x₀, t₁, x₁) == xᵢ
@test _hermite_extrapolation!(_get_v(ode), tₚ, t₀, xₚ, x₀, t₁, x₁, ẋ₁) == (xᵢ, ẋᵢ)

# Euler Extrapolation for ODEs

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 0)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-1

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 1)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-2

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 2)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-4

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 3)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-6

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 4)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-8

_euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 5)
# println(xᵢ, xₙ, xᵢ .- xₙ)
@test xᵢ ≈ xₙ atol=1E-10

for i in 0:5
    extrap1 = EulerExtrapolation(ode, i)
    extrap2 = EulerExtrapolationODE(_get_v(ode), i)
    @test extrap1 == extrap2
    @test evaluate!(extrap1, t₀, t₁, x₀, xᵢ) == _euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, i)
    @test evaluate!(extrap2, t₀, t₁, x₀, xᵢ) == _euler_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, i)
end

# Midpoint Extrapolation for ODEs

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 0)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-4

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 1)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-8

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 2)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-12

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 3)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-15

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 4)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-16

_midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, 5)
# println(xᵢ, xₙ, qᵢ .- qₙ)
@test xᵢ ≈ xₙ atol=1E-15

for i in 0:5
    extrap1 = MidpointExtrapolation(ode, i)
    extrap2 = MidpointExtrapolationODE(_get_v(ode), i)
    @test extrap1 == extrap2
    @test evaluate!(extrap1, t₀, t₁, x₀, xᵢ) == _midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, i)
    @test evaluate!(extrap2, t₀, t₁, x₀, xᵢ) == _midpoint_extrapolation_ode!(_get_v(ode), t₀, t₁, x₀, xᵢ, i)
end

# Create PODE Solution Arrays

q₀ = pode.q₀[1]
p₀ = pode.p₀[1]

qᵢ = zero(q₀)
pᵢ = zero(p₀)

q̇ₚ = zero(q₀)
q̇₀ = zero(q₀)
q̇ₙ = zero(q₀)
q̇ᵢ = zero(q₀)

ṗₚ = zero(p₀)
ṗ₀ = zero(p₀)
ṗₙ = zero(p₀)
ṗᵢ = zero(p₀)

# Compute Reference Solution for PODEs

ref_prev = integrate(pode, PartitionedTableauGauss(8), -Δt, 1)
ref_next = integrate(pode, PartitionedTableauGauss(8), +Δt, 1)

qₚ = ref_prev.q[end]
pₚ = ref_prev.p[end]

qₙ = ref_next.q[end]
pₙ = ref_next.p[end]

pode.v(tₚ, qₚ, pₚ, q̇ₚ, pode.parameters)
pode.v(t₀, q₀, p₀, q̇₀, pode.parameters)
pode.v(tₙ, qₙ, pₙ, q̇ₙ, pode.parameters)

pode.f(tₚ, qₚ, pₚ, ṗₚ, pode.parameters)
pode.f(t₀, q₀, p₀, ṗ₀, pode.parameters)
pode.f(tₙ, qₙ, pₙ, ṗₙ, pode.parameters)

# Midpoint Extrapolation for PODEs

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 0)
# println(0, qᵢ, qₙ, qᵢ .- qₙ)
# println(0, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-6
@test pᵢ ≈ pₙ atol=1E-4

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 1)
# println(1, qᵢ, qₙ, qᵢ .- qₙ)
# println(1, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-10
@test pᵢ ≈ pₙ atol=1E-8

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 2)
# println(2, qᵢ, qₙ, qᵢ .- qₙ)
# println(2, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-14
@test pᵢ ≈ pₙ atol=1E-12

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 3)
# println(3, qᵢ, qₙ, qᵢ .- qₙ)
# println(3, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-15
@test pᵢ ≈ pₙ atol=1E-16

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 4)
# println(4, qᵢ, qₙ, qᵢ .- qₙ)
# println(4, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-16
@test pᵢ ≈ pₙ atol=1E-16

_midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 5)
# println(5, qᵢ, qₙ, qᵢ .- qₙ)
# println(5, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-15
@test pᵢ ≈ pₙ atol=1E-16

for i in 0:5
    extrap1 = MidpointExtrapolation(pode, i)
    extrap2 = MidpointExtrapolationPODE(_get_v(pode), _get_f(pode), i)
    @test extrap1 == extrap2
    @test evaluate!(extrap1, t₀, t₁, q₀, qᵢ, p₀, pᵢ) == _midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, i)
    @test evaluate!(extrap2, t₀, t₁, q₀, qᵢ, p₀, pᵢ) == _midpoint_extrapolation_pode!(_get_v(pode), _get_f(pode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, i)
end

# Create IODE Solution Arrays

q₀ = iode.q₀[1]
p₀ = iode.p₀[1]

qᵢ = zero(q₀)
pᵢ = zero(p₀)

q̇ₚ = zero(q₀)
q̇₀ = zero(q₀)
q̇ₙ = zero(q₀)
q̇ᵢ = zero(q₀)

ṗₚ = zero(p₀)
ṗ₀ = zero(p₀)
ṗₙ = zero(p₀)
ṗᵢ = zero(p₀)

# Compute Reference Solution for IODEs

ref_prev = integrate(iode, IntegratorVPRKpSymmetric(iode, TableauVPGLRK(8), -Δt), 1)
ref_next = integrate(iode, IntegratorVPRKpSymmetric(iode, TableauVPGLRK(8), +Δt), 1)

qₚ = ref_prev.q[end]
pₚ = ref_prev.p[end]

qₙ = ref_next.q[end]
pₙ = ref_next.p[end]

iode.v̄(tₚ, qₚ, q̇ₚ, iode.parameters)
iode.v̄(t₀, q₀, q̇₀, iode.parameters)
iode.v̄(tₙ, qₙ, q̇ₙ, iode.parameters)

iode.f̄(tₚ, qₚ, q̇ₚ, ṗₚ, iode.parameters)
iode.f̄(t₀, q₀, q̇₀, ṗ₀, iode.parameters)
iode.f̄(tₙ, qₙ, q̇ₙ, ṗₙ, iode.parameters)

# Midpoint Extrapolation for IODEs

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 0)
# println(0, qᵢ, qₙ, qᵢ .- qₙ)
# println(0, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-4
@test pᵢ ≈ pₙ atol=1E-4

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 1)
# println(1, qᵢ, qₙ, qᵢ .- qₙ)
# println(1, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-8
@test pᵢ ≈ pₙ atol=1E-8

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 2)
# println(2, qᵢ, qₙ, qᵢ .- qₙ)
# println(2, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-12
@test pᵢ ≈ pₙ atol=1E-12

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 3)
# println(3, qᵢ, qₙ, qᵢ .- qₙ)
# println(3, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-15
@test pᵢ ≈ pₙ atol=1E-16

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 4)
# println(4, qᵢ, qₙ, qᵢ .- qₙ)
# println(4, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-16
@test pᵢ ≈ pₙ atol=1E-16

_midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, 5)
# println(5, qᵢ, qₙ, qᵢ .- qₙ)
# println(5, pᵢ, pₙ, pᵢ .- pₙ)
@test qᵢ ≈ qₙ atol=1E-15
@test pᵢ ≈ pₙ atol=1E-16

for i in 0:5
    extrap1 = MidpointExtrapolation(iode, i)
    extrap2 = MidpointExtrapolationIODE(_get_v̄(iode), _get_f̄(iode), i)
    @test extrap1 == extrap2
    @test evaluate!(extrap1, t₀, t₁, q₀, qᵢ, p₀, pᵢ) == _midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, i)
    @test evaluate!(extrap2, t₀, t₁, q₀, qᵢ, p₀, pᵢ) == _midpoint_extrapolation_iode!(_get_v̄(iode), _get_f̄(iode), t₀, t₁, q₀, qᵢ, p₀, pᵢ, i)
end


text:
include(joinpath(@__DIR__, "config.jl"))
include(joinpath(@__DIR__, "..", "gmm_utils.jl"))
using JuMP, Gurobi, MLKernels

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)

### data settings ###
data_dir = "gmm"
data_info = "gmm_outlier_detection"

data_gen_dimensions = [2, 4, 6, 8, 10]
data_gen_num_gaussians = [3, 5, 7]
data_gen_repetions = 100
data_gen_num_data = 1000
data_gen_threshold = 0.1
data_gen_num_gaussians_train = [1]
data_gen_settings = DataGenGMM(data_gen_dimensions, data_gen_num_gaussians,
                               data_gen_repetions, data_gen_num_data, data_gen_threshold,
                               data_gen_num_gaussians_train)

### learning scenario ###
initial_pool_strategy = [("Pnin", Dict(:n => 25))]
split_strategy = [("Sl", Dict())]

#### models ####
models = [Dict(:type => :SVDDneg, :param => Dict{Symbol, Any}())]
init_strategies = [WangCombinedInitializationStrategy(solver, 2.0.^range(-4, stop=4, step=1.0), FixedCStrategy(1.0))]
classify_precision = SVDD.OPT_PRECISION

#### oracle ####
# the oracle is initialized during the data generation process

#### query strategies ####
basline_query_strategies = [
    Dict(:type => :RandomQss, :param => Dict{Symbol, Any}(:epsilon => 0.1)),
    Dict(:type => :RandomOutlierQss, :param => Dict{Symbol, Any}(:epsilon => 0.1))
]

fancy_query_strategies = [
    Dict(:type => :DecisionBoundaryQss, :param => Dict{Symbol, Any}()),
    Dict(:type => :ExplorativeMarginQss, :param => Dict{Symbol, Any}(:solver => solver, :lambda => 1.0, :use_penalty => false)),
    Dict(:type => :ExplorativeMarginQss, :param => Dict{Symbol, Any}(:solver => solver, :lambda => 1.0, :use_penalty => true)),
]

qs_optimizers = [
    ParticleSwarmOptimization(eps=1.0, swarmsize=1_000, maxiter=1_000),
    BlackBoxOptimization(:dxnes, eps=1.0),
]

query_strategies = basline_query_strategies

for qs in fancy_query_strategies
    for o in qs_optimizers
        qs_parametrized = deepcopy(qs)
        qs_parametrized[:param][:optimizer] = o
        push!(query_strategies, qs_parametrized)
    end
end

num_al_iterations = 100
num_resamples_initial_pool = 1
exp_name = "eval_part_1"

max_stars_repo_name: englhardt/des-evaluation

max_issues_repo_name: englhardt/des-evaluation

max_forks_repo_name: englhardt/des-evaluation


text:
import Dates
import JSON

function compare_http_date_header(header_value::String, timestamp_request_completed::Dates.DateTime) :: Nothing
    header_value_timestamp::Dates.DateTime = Dates.DateTime(split(header_value, " UTC")[1], "e, d u Y H:M:S")
    @test header_value_timestamp <= timestamp_request_completed
    nothing
end

function compare_http_header(headers::Array, key::String, value::String) :: Nothing
    @test header_get_value(headers::Array, key::String)==value
    nothing
end

function header_get_value(headers::Array, key::String) :: String
    for item in headers
        if item[1]==key
            return item[2]
        end
    end
end

function extract_html_body_content(html_body::Array{UInt8,1}) :: String
    return split(
        split(String(html_body), "
")[2],
        "
"
    )[1]
end

function extract_json_content(json_body::Array{UInt8,1}) :: String
    return String(json_body)
end

function parse_and_test_request(r::HTTP.Messages.Response, status::Int64, headers::Dict{String, String}, content_length::Int64, is_json_body::Bool, body::Any)
    timestamp_request_completed::Dates.DateTime = Dates.now(Dates.UTC)

    @test r.status==status

    for (key, value) in headers
        compare_http_header(r.headers, key, value)
    end
    # TODO: fix Windows HTML file longer due to line final char
    if !Sys.iswindows()
        compare_http_header(r.headers, "Content-Length", string(content_length))
    end
    if header_get_value(r.headers, "Content-Type")=="application/json"
        compare_http_header(r.headers, "Access-Control-Allow-Origin", Dance.Configuration.Settings[:api_access_control_allow_origin])
    end
    compare_http_date_header(header_get_value(r.headers, "Date"), timestamp_request_completed)

    if is_json_body
        if header_get_value(r.headers, "Content-Type")=="text/html; charset=UTF-8"
            @test JSON.parse(extract_html_body_content(r.body))==body
        else
            @test JSON.parse(extract_json_content(r.body))==body
        end
    else
        if header_get_value(r.headers, "Content-Type")=="text/html; charset=UTF-8"
            @test extract_html_body_content(r.body)==body
        else
            # Static file case
            @test r.body==body
        end
    end

    nothing
end


text:
# using BlossomV
# using Metis
# using SparseArrays
# using LightGraphs: add_edge!
# using SimpleWeightedGraphs
# using MetaGraphs
# using ..Utils

# const δINT = 999999
# const MIN_INT = 1
# const MAX_INT = δINT + MIN_INT

# #############
# # Metis top down method
# #############

# struct MetisContext
#     info::Matrix{Int64}
# end

# MetisContext(mi::Matrix{Float64}) = MetisContext(to_long_mi(mi, MIN_INT, MAX_INT))

# #Add edge weights to Metis.jl
# using Metis: idx_t, ishermitian

# struct WeightedGraph
#     nvtxs::idx_t
#     xadj::Vector{idx_t}
#     adjncy::Vector{idx_t}
#     adjwgt::Vector{idx_t} # edge weights
#     WeightedGraph(nvtxs, xadj, adjncy, adjwgt) = new(nvtxs, xadj, adjncy, adjwgt)
# end

# function my_graph(G::SparseMatrixCSC; check_hermitian=true)
#     if check_hermitian
#         ishermitian(G) || throw(ArgumentError("matrix must be Hermitian"))
#     end
#     N = size(G, 1)
#     xadj = Vector{idx_t}(undef, N+1)
#     xadj[1] = 1
#     adjncy = Vector{idx_t}(undef, nnz(G))
#     adjncy_i = 0
#     adjwgt = Vector{idx_t}(undef, nnz(G))
#     @inbounds for j in 1:N
#         n_rows = 0
#         for k in G.colptr[j] : (G.colptr[j+1] - 1)
#             i = G.rowval[k]
#             if i != j # don't include diagonal elements
#                 n_rows += 1
#                 adjncy_i += 1
#                 adjncy[adjncy_i] = i
#                 adjwgt[adjncy_i] = G[i, j]
#             end
#         end
#         xadj[j+1] = xadj[j] + n_rows
#     end
#     resize!(adjncy, adjncy_i)
#     resize!(adjwgt, adjncy_i)
#     return WeightedGraph(idx_t(N), xadj, adjncy, adjwgt)
# end

# function my_partition(G::WeightedGraph, nparts::Integer; alg = :KWAY)
#     part = Vector{Metis.idx_t}(undef, G.nvtxs)
#     edgecut = fill(idx_t(0), 1)
#     if alg === :RECURSIVE
#         Metis.METIS_PartGraphRecursive(G.nvtxs, idx_t(1), G.xadj, G.adjncy, C_NULL, C_NULL, G.adjwgt,
#                                  idx_t(nparts), C_NULL, C_NULL, Metis.options, edgecut, part)
#     elseif alg === :KWAY
#         Metis.METIS_PartGraphKway(G.nvtxs, idx_t(1), G.xadj, G.adjncy, C_NULL, C_NULL, G.adjwgt,
#                             idx_t(nparts), C_NULL, C_NULL, Metis.options, edgecut, part)
#     else
#         throw(ArgumentError("unknown algorithm $(repr(alg))"))
#     end
#     return part
# end

# my_partition(G, nparts; alg = :KWAY) = my_partition(my_graph(G), nparts, alg = alg)

# "Metis top down method"
# function metis_top_down(vars::Vector{Var}, context::MetisContext)::Tuple{Vector{Var}, Vector{Var}}

#     vertices = sort(collect(vars))
#     sub_context = context.info[vertices, vertices]
#     len = length(vertices)
#     for i in 1 : len
#         sub_context[i, i] = 0
#     end
#     g = convert(SparseMatrixCSC, sub_context)
#     partition = my_partition(my_graph(g), 2, alg = :RECURSIVE)

#     subsets = (Vector{Var}(), Vector{Var}())
#     for (index, p) in enumerate(partition)
#         push!(subsets[p], vertices[index])
#     end

#     return subsets
# end

# function metis_top_down_curry(context::MetisContext)
#     f(vars) = metis_top_down(vars, context)
#     return f
# end

# #############
# # Blossom bottom up method
# #############

# # TODO change API to DisjointSet
# mutable struct BlossomContext
#     variable_sets::Vector{Vector{Var}}
#     partition_id::Vector{Int64} # map vars to index in variable_sets
#     info::Matrix
# end

# BlossomContext(vars::Vector{Var}, mi::Matrix{Float64}) =
#     BlossomContext( [[v] for v in sort(collect(vars))],
#                     collect(1 : length(vars)),
#                     round.(Int64, 1000001 .+ to_long_mi(mi, -1, -1000000)))
#                     #mi)

# "Blossom bottom up method, vars are not used"
# function blossom_bottom_up!(vars::Vector{Var}, context::BlossomContext)::Vector{Tuple{Var, Var}}

#     "even number of nodes, use blossomv alg"
#     function blossom_bottom_up_even!(vars::Vector{Var}, context::BlossomContext; update = true)::Tuple{Vector{Tuple{Var, Var}}, Int64}
#         "1. calculate pMI"
#         pMI = set_mutual_information(context.info, context.variable_sets)
#         #pMI = 1000001 .+ to_long_mi(pMI, -1, -1000000)
#         pMI = round.(Int64, pMI)

#         "2. solve by blossomv alg"
#         len = length(context.variable_sets)
#         m = Matching(len)
#         for i in 1 : len, j in i + 1 : len
#             add_edge(m, i - 1, j - 1, pMI[i, j]) # blossomv index start from 0
#         end

#         solve(m)
#         all_matches = Set{Tuple{Var, Var}}()
#         for v in 1 : len
#             push!(all_matches, order_asc(v, get_match(m, v - 1) + 1))
#         end

#         "3. calculate scores, map index to var"
#         all_matches = Vector(collect(all_matches))
#         score = 0

#         for i in 1 : length(all_matches)
#             (x, y) = all_matches[i]
#             score += pMI[x, y]
#             all_matches[i] = (context.variable_sets[x][1], context.variable_sets[y][1])
#         end
#         all_matches

#         "4. update context when called by outer layer"
#         if update
#             update_context(all_matches, context)
#         end

#         return (all_matches, score)
#     end

#     "odd number of nodes, try every 2 combinations"
#     function blossom_bottom_up_odd!(vars::Vector{Var}, context::BlossomContext)::Tuple{Vector{Tuple{Var, Var}}, Int64}

#         "1. try all len - 1 conditions, find best score(minimun cost)"
#         (best_matches, best_score) = (Vector{Tuple{Var, Var}}(), typemax(Int64))

#         for index in 1 : length(context.variable_sets)
#             sub_context = deepcopy(context)
#             set = copy(sub_context.variable_sets[index])
#             deleteat!(sub_context.variable_sets, index)

#             (matches, score) = blossom_bottom_up_even!(vars, sub_context; update = false)
#             if score < best_score
#                 (best_matches, best_score) = (matches, score)
#             end

#             insert!(sub_context.variable_sets, index, set)
#         end

#         "2. update information"
#         update_context(best_matches, context)
#         return (best_matches, best_score)

#     end

#     function update_context(matches::Vector{Tuple{Var, Var}}, context::BlossomContext)
#         for (x, y) in matches
#             y_partition = copy(context.variable_sets[context.partition_id[y]])
#             context.variable_sets[context.partition_id[y]] = Vector()
#             foreach(value -> push!(context.variable_sets[context.partition_id[x]], value), y_partition)
#         end

#         context.variable_sets = [x for x in context.variable_sets if x != []]
#         for index in 1 : length(context.variable_sets)
#             for y in context.variable_sets[index]
#                 context.partition_id[y] = index
#             end
#         end
#     end

#     if length(vars) % 2 == 0
#         (matches, score) = blossom_bottom_up_even!(vars, context)
#     else
#         (matches, score) = blossom_bottom_up_odd!(vars, context)
#     end

#     for (left, right) in matches
#         deleteat!(vars, findall(v->v==right, vars))
#     end
#     return matches
# end

# function blossom_bottom_up_curry(context::BlossomContext)
#     f(vars) = blossom_bottom_up!(vars, context)
#     return f
# end

# # refactor later
# function learn_vtree_bottom_up(train_x::PlainXData; α)
#     (_, mi) = mutual_information(feature_matrix(train_x), Data.weights(train_x); α = α)
#     vars = Var.(collect(1:num_features(train_x)))
#     context = BlossomContext(vars, mi)
#     vtree = bottom_up_vtree(PlainVtreeNode, vars, blossom_bottom_up_curry(context))
# end

# #############
# # Test method
# #############

# "Test top down method, split nodes by ascending order, balanced"
# function test_top_down(vars::Vector{Var})::Tuple{Vector{Var}, Vector{Var}}
#     sorted_vars = sort(collect(vars))
#     len = length(sorted_vars)
#     len1 = Int64(len % 2 == 0 ? len // 2 : (len - 1) // 2)
#     return (sorted_vars[1 : len1], sorted_vars[len1 + 1 : end])
# end

# "Test bottom up method, split nodes by ascending order, balanced"
# function test_bottom_up!(vars::Vector{Var})::Vector{Tuple{Var, Var}}
#     sorted_vars = sort(collect(vars))
#     len = length(sorted_vars)
#     len1 = Int64(len % 2 == 0 ? len // 2 : (len - 1) // 2)
#     matches = Vector{Tuple{Var, Var}}()
#     for i in 1 : len1
#         push!(matches, (sorted_vars[2 * i - 1], sorted_vars[2 * i]))
#         pop!(vars, sorted_vars[2 * i])
#     end
#     return matches
# end


text:
"""
    insert!(tree::RTree, pt::Point, [id], val)
    insert!(tree::RTree, br::Rect, [id], val)
    insert!(tree::RTree, elem::SpatialElem)

Inserts `val` value identified by `br` bounding box (or point `pt`) and
`id` (if tree elememnts support `HasID` trait) into the `tree`.
"""
function Base.insert!(tree::RTree, el::Any)
    check_eltype_rtree(el, tree)
    _insert!(tree, el, 1)
    tree.nelems += 1
    tree.nelem_insertions += 1
    return tree
end

Base.insert!(tree::RTree{T,N,SpatialElem{T,N,K,V}}, br::Rect{T,N}, id::Any, val::Any) where {T,N,K,V} =
    insert!(tree, SpatialElem{T,N,K,V}(br, id, val))

Base.insert!(tree::RTree{T,N,SpatialElem{T,N,K,V}},
             pt::Point{T,N}, id::Any, val::Any) where {T,N,K,V} =
    insert!(tree, SpatialElem{T,N,K,V}(Rect(pt), id, val))

Base.insert!(tree::RTree{T,N,SpatialElem{T,N,Nothing,V}},
             br::Union{Rect{T,N}, Point{T,N}}, val::Any) where {T,N,V} =
    insert!(tree, br, nothing, val)

# inserts the child into the node rebalancing R-tree as needed
# returns true if MBR was changed and upstream R-tree readjustments are required
function _insert!(node::Node, child::Any, con::SubtreeContext)
    if length(node) < capacity(node, con.tree) # underfull node
        oldmbr = mbr(node) # this has to happen before _attach!() modifies node's MBR
        _attach!(node, child, con.tree)
        if !in(mbr(child), oldmbr) && hasparent(node)
            #@debug "_insert!(): mbr change"
            _updatembr!(parent(node), pos_in_parent(node), oldmbr, con)
            return true
        else
            return false
        end
    elseif variant(con.tree) == RTreeStar && (node !== con.tree.root) && !isoverflow(con, level(node))
        return _insert!_fullnode_rstar(node, child, con)
    else
        return _insert!_fullnode(node, child, con)
    end
end

function _insert!(con::SubtreeContext, node::Any, lev::Integer)
    #@debug "_insert!(con): len=$(length(con.tree)) node_lev=$(node isa Node ? level(node) : 0) node_len=$(node isa Node ? length(node) : 0)"
    par = choose_subtree(con.tree.root, lev, mbr(node), con.tree)
    @assert level(par) == lev
    _insert!(par, node, con)
end

_insert!(tree::RTree, node::Any, lev::Int) =
    _insert!(SubtreeContext(tree), node, lev)

# index of the child that would have the least MBR enlargement by br
function find_least_enlargement(node::Node, br::Rect)
    min_enl = Inf
    best_ix = 0
    best_area = Inf
    for i in eachindex(children(node))
        child_mbr = mbr(@inbounds node[i])
        child_area = area(child_mbr)
        enl = combined_area(child_mbr, br) - child_area
        if enl < min_enl || (min_enl == enl && child_area < best_area)
            min_enl = enl
            best_area = child_area
            best_ix = i
        end
    end
    return best_ix
end

# index of the child that would have the least MBR overlap with br
function find_least_overlap(node::Node, br::Rect, tree::RTree)
    min_enl = Inf
    best_area = Inf
    best_ix = 0
    enls = Vector{Float64}(undef, length(node)) # FIXME use pool?
    # find combined region and enlargement of every entry and store it.
    for (i, child) in enumerate(children(node))
        child_area = area(mbr(child))
        enls[i] = enl = combined_area(br, mbr(child)) - child_area
        if enl < min_enl || (enl == min_enl && child_area < best_area)
            min_enl = enl
            best_area = child_area
            best_ix = i
        end
    end

    abs(min_enl) <= eps(min_enl) && return best_ix # no enlargement

    if length(node) > tree.nearmin_overlap
        # sort entries in increasing order of enlargement
        enl_order = sortperm(enls) # FIXME reuse from pool
        @assert min_enl == enls[enl_order[1]] <= enls[enl_order[end]]
        niter = tree.nearmin_overlap
    else
        niter = length(node)
        enl_order = collect(1:niter) # FIXME reuse from pool
    end

    # calculate overlap of most important original entries (near minimum overlap cost).
    min_delta_olap = Inf
    best_ix = 0
    best_area = Inf
    best_enl = Inf
    for i in 1:niter
        node_ix = enl_order[i]
        node_mbr = mbr(node[node_ix])
        comb_mbr = combine(node_mbr, br)
        delta_olap = 0.0

        for (j, child) in enumerate(children(node))
            if node_ix != j
                olap = overlap_area(comb_mbr, mbr(child))
                if olap > 0.0
                    delta_olap += olap - overlap_area(node_mbr, mbr(child))
                end
            end
        end

        node_area = area(node_mbr)
        # keep the one with the least, in the order:
        #   delta overlap, enlargement, area
        if (delta_olap < min_delta_olap) || (delta_olap == min_delta_olap &&
           (enls[node_ix] < best_enl || (enls[node_ix] == best_enl &&
            node_area < best_area)))
            min_delta_olap = delta_olap
            best_ix = node_ix
            best_area = node_area
            best_enl = enls[node_ix]
        end
    end
    return best_ix
end

# choose the node (for inserting the given mbr at specified level)
function choose_subtree(node::Node, lev::Int, br::Rect, tree::RTree)
    while level(node) > lev
        if variant(tree) == RTreeLinear || variant(tree) == RTreeQuadratic
            child_ix = find_least_enlargement(node, br)
        elseif variant(tree) == RTreeStar
            child_ix = level(node) == 1 ?
                    find_least_overlap(node, br, tree) : # for the node pointing to leaves
                    find_least_enlargement(node, br)
        else
            throw(SpatialIndexException("RTree variant not supported"))
        end
        node = node[child_ix]
    end
    level(node) == lev ||
        throw(SpatialIndexException("Couldn't find the node of level $lev in RTree of height $(height(tree))"))
    return node
end

# R*-tree-specific insertion into the full node
# splits the node and rebalances the R-tree
function _insert!_fullnode_rstar(node::Node, child::Any, con::SubtreeContext)
    #@debug "_insert!_fullnode_rstar(): lev=$(level(node)) len=$(length(node))"
    @assert isa(con.level_overflow, BitVector)
    @assert hasparent(node)

    setoverflow!(con, level(node)) # the R* insert could only be done once per level
    oldmbr = mbr(node)
    _attach!(node, child, con.tree, force=true) # force-push the child exceeding the capacity

    node_center = center(mbr(node))
    # calculate relative distance of every entry from the node MBR (ignore square root.)
    children_dist = sqrdistance.(Ref(node_center), center.(mbr.(children(node)))) # FIXME reuse array?
    # return the indices sorted by the increasing order of distances.
    # Since children_dist is sorted in ascending order
    # it would match the reinsertion order suggested in the paper.
    children_order = sortperm(children_dist) # FIXME reuse array?
    nkeep = length(node) - floor(Int, length(node) * con.tree.reinsert_factor)

    newnode = acquire(con.tree, typeof(node), level(node))   # the replacement node
    for i in 1:nkeep
        _attach!(newnode, node[children_order[i]], con.tree)
    end

    # Divertion from R*-Tree algorithm proposed in libspatialindex:
    # First update the path to the root, then start reinserts, to avoid complicated handling
    # of changes to the same node from multiple insertions.
    #@debug "_insert!_fullnode_rstar(): replace the node"
    _replace!(node, newnode, oldmbr, con) # replace the node (so (nkeep+1):... children are detached from the tree)
    #@debug "_insert!_fullnode_rstar(): reinsert $(length((nkeep+1):length(children_order))) children of lev=$(level(node)) node"
    # don't use par during/after reinsertion, it might have changed
    for i in (nkeep+1):length(children_order)
        child = node[children_order[i]]
        #@debug "_insert!_fullnode_rstar(): reinsert #$i child (type=$(typeof(child)) lev=$(child isa Node ? level(child) : "data") parent=$(isa(child, Node) && hasparent(child)))"
        _insert!(con, child, level(node))
        con.tree.nnode_reinsertions += 1
    end

    release(con.tree, node) # return node to the pool
    return true
end

# R-tree-specific insertion into the full node
# splits the node and rebalances the R-tree
function _insert!_fullnode(node::Node, child::Any, con::SubtreeContext)
    #@debug "_insert!_fullnode(): lev=$(level(node)) len=$(length(node)) parent=$(hasparent(node))"
    oldmbr = mbr(node)
    _attach!(node, child, con.tree, force=true) # temporary overflow the node capacity
    n1, n2 = _split!(node, con.tree)

    if node === con.tree.root
        #@debug "_insert!_fullnode(): newroot tree_height=$(height(con.tree)) tree_len=$(length(con.tree))"
        # the node must be the root, make the new root and attach n1, n2 to it
        newroot = acquire(con.tree, Branch, level(node) + 1)
        _attach!(newroot, n1, con.tree)
        _attach!(newroot, n2, con.tree)
        # tree grows
        con.tree.root = newroot
        con.tree.nnodes_perlevel[level(node)] = 2
        push!(con.tree.nnodes_perlevel, 1)
    else # non-root, n1 and n2 replace the node in its parent
        _replace!(node, n1, n2, oldmbr, con)
    end
    return true
end


text:
using NQCDynamics
using Test
using RecursiveArrayTools: ArrayPartition
using ComponentArrays: ComponentVector

atoms = Atoms([:H])
model = NQCModels.Harmonic()
sim = Simulation(atoms, model)

positions = [randn(1, length(atoms)) for i=1:10]
velocities = [randn(1, length(atoms)) for i=1:10]
distribution = DynamicalDistribution(positions, velocities, (1, 1))
tspan = (0.0, 10.0)

@testset "run_ensemble" begin
    out = run_ensemble(sim, tspan, distribution;
        output=(:position), dt=1, trajectories=10)
    @test length(out) == 10
    @test out[1].position isa Vector{<:Matrix}
    @test out[1].t isa Vector{<:Real}
end

@testset "run_ensemble" begin
    out = run_ensemble(sim, tspan, distribution;
        output=Ensembles.OutputFinal(), dt=1, trajectories=10, reduction=:append)
    @test out isa Vector{<:ArrayPartition}
end

@testset "run_ensemble, reduction=$reduction" for reduction ∈ (:sum, :mean)
    out = run_ensemble(sim, tspan, distribution;
        output=Ensembles.OutputFinal(), dt=1, trajectories=10, reduction=reduction)
    @test out isa ComponentVector
end


text:
# abstract type for LaguerrePolynomial and LaguerreFunction
abstract type AbstractLaguerre{T <: Number} end

Base.iterate(la::AbstractLaguerre{T}, state = 1) where T = state > la.len ? nothing : ((la.orders[state], la.coeffs[state]), state + 1)

function lashow(la::AbstractLaguerre{T}) where T
    # return a string representation of a linear Combination of Laguerre polynomials
    st = ""
    if la isa LaguerrePolynomial{T}
        name = "L"
    elseif la isa LaguerreFunction{T}
        name = "Lf"
    end
    for i in 1:la.len - 1
        st *= "($(la.coeffs[i]))" * name * "$(la.orders[i]) + "
    end
    # add the last one
    st *= "($(la.coeffs[end]))" * name * "$(la.orders[end])"
    return st
end

Base.show(io::IO, la::AbstractLaguerre{T}) where T = println(lashow(la))

# implement addition of LaguerrePolynomial or LaguerreFunction .
function Base.:+(self::AbstractLaguerre{T}, other::AbstractLaguerre{S}) where {T,S}
    # addition of two Laguerre Polynomials or functions.
    # need to be the same type.
    nameof(typeof(self)) == nameof(typeof(other)) || throw(ErrorException("cannot add LaguerrePolynomial to a LaguerreFunction"))
    l = addla(self, other)
    return l
end

function addla(l1::AbstractLaguerre{T}, l2::AbstractLaguerre{S}) where {T,S}
    # implementation of addition of Laguerre Polynomials/Functions.
    # need to copy to create a new array.
    orders = copy(l1.orders)
    coefs = copy(l1.coeffs)
    # type promotion
    R = promote_type(T, S)
    # convert coeffs to type R.
    coefs = R.(coefs)
    # iterate over l2 add each to l1
    for (order, coef) in l2
        if order in orders
            id = findall(orders .== order)[1]
            coefs[id] += coef
        else
            append!(orders, order)
            append!(coefs, R(coef))
        end
    end
    # get the symbol of the typename.
    typename = nameof(typeof(l1))
    # constructor expression.
    expr = :($typename($orders, $coefs))
    return eval(expr)
end

# multiplication of LaguerrePolynomial/LaguerreFunction with a number.
function Base.:*(l::AbstractLaguerre{T}, num::S) where {S <: Number} where T
    # product of LaguerrePolynomial with a constant.
    R = promote_type(T, S)
    coefs = R.(copy(l.coeffs))
    coefs *= num
    # get the symbol of the typename.
    typename = nameof(typeof(l))
    # constructor expression.
    expr = :($typename($(l.orders), $coefs))
    return eval(expr)
end

Base.:*(num::S, l::AbstractLaguerre{T}) where T where S <: Number = l * num

function Base.:-(l1::AbstractLaguerre{T}, l2::AbstractLaguerre{S}) where {T,S}
    # subtraction of two LaguerrePolynomial or LaguerreFunction
    # need to be the same type.
    # no need to check it here, because addition will check .
    R = promote_type(T, S)
    minus1 = R(-1)
    return l1 + minus1 * l2
end

function Base.:/(la::AbstractLaguerre{T}, num::S) where {S <: Number} where T
    num != 0 || throw(ErrorException("division by zero not allowed."))
    # divide a LaguerrePolynomial by a number can be done.
    R = promote_type(T, S)
    coefs = R.(copy(la.coeffs))
    coefs /= num
    # get the symbol of the typename.
    typename = nameof(typeof(la))
    # constructor expression.
    expr = :($typename($(la.orders), $coefs))
    return eval(expr)
end

# check equality of two LaguerrePolynomials or LaguerreFunctions
function Base.:(==)(l1::AbstractLaguerre{T}, l2::AbstractLaguerre{S}) where {T,S}
    if (nameof(typeof(l1)) != nameof(typeof(l2)))
        return false
    end
    return (l1.orders == l2.orders) && (l1.coeffs == l2.coeffs)
end

function simplify(l::AbstractLaguerre{T}) where T
    # remove zero elements from LaguerrePolynomial or LaguerreFunction
    # find nonzero coeffs
    idx = findall(l.coeffs .!= 0)
    if isempty(idx)
        # all zeros.
        # just return the zero order term with 0 coefficient.
        coeffs = 0
        orders = 0
    elseif length(idx) == l.len
        # all coefficients are nonzero.
        # nothing needs to be done.
        return l
    else
        # some are zero.
        coeffs = l.coeffs[idx]
        orders = l.orders[idx]
    end
    # get the symbol of the typename.
    typename = nameof(typeof(l))
    # constructor expression.
    expr = :($typename($orders, $coeffs))
    return eval(expr)
end


text:
# This file contains functions/methods useful for transforming coordinates in
# SMLD structures.

"""
    addoffset!(smld::SMLD2D, offset::Vector{Float64})

Add an offset elementwise to each localization `smld.y` and `smld.x`.

# Inputs
-`offset`: Two element vector describing the offset. ([offset_y; offset_x])
"""
function addoffset!(smld::SMLD2D, offset::Vector{Float64})
    smld.y .+= offset[1]
    smld.x .+= offset[2]
end

"""
    smld_offset = addoffset(smld::SMLD2D, offset::Vector{Float64})

Add an offset elementwise to each localization `smld.y` and `smld.x`.

# Inputs
-`offset`: Two element vector describing the offset. ([offset_y; offset_x])
"""
function addoffset(smld::SMLD2D, offset::Vector{Float64})
    smld_offset = deepcopy(smld)
    smld_offset.y .+= offset[1]
    smld_offset.x .+= offset[2]

    return smld_offset
end


text:
@discord_object struct EmbedAuthor
    name::String
    url::String
    icon_url::String
    proxy_icon_url::String
end