text:
#---------------------------------------------------------
# # [Random Disks](@id 09-disk)
#---------------------------------------------------------

#=
This page illustrates the `disk_phantom_params` method in the Julia package
[`ImagePhantoms`](https://github.com/JuliaImageRecon/ImagePhantoms.jl).

This page was generated from a single Julia file:
[09-disk.jl](@__REPO_ROOT_URL__/09-disk.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`09-disk.ipynb`](@__NBVIEWER_ROOT_URL__/09-disk.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`09-disk.ipynb`](@__BINDER_ROOT_URL__/09-disk.ipynb).

# ### Setup

# Packages needed here.

using ImagePhantoms: Ellipse, phantom, disk_phantom_params
using ImageGeoms: ImageGeom
using MIRTjim: jim, prompt
using Plots # @animate, gif

# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

#=
For training machine-learning methods for image reconstruction,
it can be helpful to have a way to generate a family of phantoms
drawn from some common distribution,
especially for debugging
or when real ground-truth data is limited.
The `disk_phantom_params` function is one example of how
one can use the tools in this package to generate such phantoms.
=#

# ### A single disk phantom image

# Make a digital image of it using `phantom` and display it.
function disk_phantom(title::String)
    (dx,dy) = (1, 1)
    (M,N) = (2^8,2^8)
    x = (-M÷2:M÷2-1) * dx
    y = (-N÷2:N÷2-1) * dy
    params = disk_phantom_params( ; rhead = () -> rand(100:105))
    objects = Ellipse(params) # vector of Ellipse objects
    oversample = 3
    img = phantom(x, y, objects, oversample)
    jim(x, y, img; title, clim=(0,1300))
end
disk_phantom("A single disk phantom realization")

# ### Several realizations

anim = @animate for i in 1:8
    disk_phantom("Realization $i")
end
gif(anim, "disk.gif", fps = 6)


text:
struct AlgebraicWrapper
    _source
end

function algebraic_version(m::Model)
    if m._jump_model===nothing
        return AlgebraicWrapper(build(m))
    else
        return AlgebraicWrapper(m._jump_model)
    end
end

function constraint_values(m::Model)
    for j=1:length(m._nlexpressions)
        var_values = []
        println("Constraint $j: ", value(m._nlexpressions[j],i->begin
            val = MPSGE.Complementarity.result_value(i)
            push!(var_values, "  $i = $val")
            return val
        end))

        for s in var_values
            println(s)
        end
    end
end

function Base.show(io::IO, m::AlgebraicWrapper)
    println(io, "Mixed complementarity problem with $(length(m._source.ext[:MCP])) constraints:")
    constraint_strings = [JuMP.nl_expr_string(m._source, JuMP.REPLMode, m._source.nlp_data.nlexpr[c.F.index]) for c in m._source.ext[:MCP]]

    column1_width = maximum(textwidth.(constraint_strings))

    for (constraint_string, c) in zip(constraint_strings, m._source.ext[:MCP])

        print(io, "  ")

        print(io, rpad(constraint_string, column1_width))

        print(io, "  ┴  ")

        if !isinf(c.lb) && c.lb==c.ub
            print(io, c.var_name, " = $(c.ub)")
        else
            print(io, isinf(c.lb) ? "" : "$(c.lb) < ", c.var_name, isinf(c.ub) ? "" : " < $(c.ub)")
        end

        println(io)
    end
end

function Base.show(io::IO, ::MIME"text/latex", m::AlgebraicWrapper)
    println(io, raw"$$ \begin{alignat*}{3}\\")
    for c in m._source.ext[:MCP]

        print(io, "& ")

        println(io, JuMP.nl_expr_string(m._source, JuMP.IJuliaMode, m._source.nlp_data.nlexpr[c.F.index]))

        print(io, raw"\quad && \perp \quad && ")

        if !isinf(c.lb) && c.lb==c.ub
            print(io, c.var_name, " = $(c.ub)")
        else
            print(io, isinf(c.lb) ? "" : "$(c.lb) <", c.var_name, isinf(c.ub) ? "" : " < $(c.ub)")
        end

        print(io, "\\\\")
    end
    println(io, raw"\end{alignat*}")
    println(io, raw" $$")
end

ext: jl


text:
function EOPEntry(arg0::jint, arg1::jdouble, arg10::AbsoluteDate, arg2::jdouble, arg3::jdouble, arg4::jdouble, arg5::jdouble, arg6::jdouble, arg7::jdouble, arg8::jdouble, arg9::ITRFVersion)
    return EOPEntry((jint, jdouble, AbsoluteDate, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, jdouble, ITRFVersion), arg0, arg1, arg10, arg2, arg3, arg4, arg5, arg6, arg7, arg8, arg9)
end

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

function get_dd_eps(obj::EOPEntry)
    return jcall(obj, "getDdEps", jdouble, ())
end

function get_dd_psi(obj::EOPEntry)
    return jcall(obj, "getDdPsi", jdouble, ())
end

function get_dx(obj::EOPEntry)
    return jcall(obj, "getDx", jdouble, ())
end

function get_dy(obj::EOPEntry)
    return jcall(obj, "getDy", jdouble, ())
end

function get_itrf_type(obj::EOPEntry)
    return jcall(obj, "getITRFType", ITRFVersion, ())
end

function get_lod(obj::EOPEntry)
    return jcall(obj, "getLOD", jdouble, ())
end

function get_mjd(obj::EOPEntry)
    return jcall(obj, "getMjd", jint, ())
end

function get_ut1_minus_utc(obj::EOPEntry)
    return jcall(obj, "getUT1MinusUTC", jdouble, ())
end

function get_x(obj::EOPEntry)
    return jcall(obj, "getX", jdouble, ())
end

function get_y(obj::EOPEntry)
    return jcall(obj, "getY", jdouble, ())
end


text:
# Defines the modules referred to in the example docs (under src/) and then builds them.
# It can be called separately to build the examples/, or as part of the test suite.
#
# It defines a set of variables (`examples_*`) that can be used in the tests.
# The `examples_root` should be used to check whether this file has already been included
# or not and should be kept unique.
isdefined(@__MODULE__, :examples_root) && error("examples_root is already defined\n$(@__FILE__) included multiple times?")

# The `Mod` and `AutoDocs` modules are assumed to exist in the Main module.
(@__MODULE__) === Main || error("$(@__FILE__) must be included into Main.")

# DOCUMENTER_TEST_EXAMPLES environment variable can be used to control which
# builds actually run. E.g. to only build the HTML deployment example build, you
# could call the make.jl file as follows:
#
#     DOCUMENTER_TEST_EXAMPLES=html julia --project test/examples/make.jl
#
EXAMPLE_BUILDS = if haskey(ENV, "DOCUMENTER_TEST_EXAMPLES")
    split(ENV["DOCUMENTER_TEST_EXAMPLES"])
else
    ["markdown", "html", "html-mathjax3", "html-local"]
end

# Modules `Mod` and `AutoDocs`
module Mod
    """
        func(x)

    [`T`](@ref)
    """
    func(x) = x

    """
        T

    [`func(x)`](@ref)
    """
    mutable struct T end
end

"`AutoDocs` module."
module AutoDocs
    module Pages
        include("pages/a.jl")
        include("pages/b.jl")
        include("pages/c.jl")
        include("pages/d.jl")
        include("pages/e.jl")
    end

    "Function `f`."
    f(x) = x

    "Constant `K`."
    const K = 1

    "Type `T`."
    mutable struct T end

    "Macro `@m`."
    macro m() end

    "Module `A`."
    module A
        "Function `A.f`."
        f(x) = x

        "Constant `A.K`."
        const K = 1

        "Type `B.T`."
        mutable struct T end

        "Macro `B.@m`."
        macro m() end
    end

    "Module `B`."
    module B
        "Function `B.f`."
        f(x) = x

        "Constant `B.K`."
        const K = 1

        "Type `B.T`."
        mutable struct T end

        "Macro `B.@m`."
        macro m() end
    end

    module Filter
        "abstract super type"
        abstract type Major end

        "abstract sub type 1"
        abstract type Minor1 <: Major end

        "abstract sub type 2"
        abstract type Minor2 <: Major end

        "random constant"
        qq = 3.14

        "random function"
        function qqq end
    end
end

# Helper functions
function withassets(f, assets...)
    src(asset) = joinpath(@__DIR__, asset)
    dst(asset) = joinpath(@__DIR__, "src/assets/$(basename(asset))")
    for asset in assets
        isfile(src(asset)) || error("$(asset) is missing")
    end
    for asset in assets
        cp(src(asset), dst(asset))
    end
    try
        f()
    finally
        @debug "Cleaning up assets" assets
        for asset in assets
            rm(dst(asset))
        end
    end
end

# Build example docs
using Documenter, DocumenterMarkdown
isdefined(@__MODULE__, :TestUtilities) || (include("../TestUtilities.jl"); using .TestUtilities)

examples_root = @__DIR__
builds_directory = joinpath(examples_root, "builds")
ispath(builds_directory) && rm(builds_directory, recursive=true)

expandfirst = ["expandorder/AA.md"]
htmlbuild_pages = Any[
    "**Home**" => "index.md",
    "Manual" => [
        "man/tutorial.md",
        "man/style.md",
    ],
    hide("hidden.md"),
    "Library" => [
        "lib/functions.md",
        "lib/autodocs.md",
        "lib/editurl.md",
    ],
    hide("Hidden Pages" => "hidden/index.md", Any[
        "Page X" => "hidden/x.md",
        "hidden/y.md",
        "hidden/z.md",
    ]),
    "Expandorder" => [
        "expandorder/00.md",
        "expandorder/01.md",
        "expandorder/AA.md",
    ],
    "unicode.md",
    "latex.md",
    "example-output.md",
]

function html_doc(build_directory, mathengine)
    @quietly withassets("images/logo.png", "images/logo.jpg", "images/logo.gif") do
        makedocs(
            debug = true,
            root  = examples_root,
            build = "builds/$(build_directory)",
            doctestfilters = [r"Ptr{0x[0-9]+}"],
            sitename = "Documenter example",
            pages = htmlbuild_pages,
            expandfirst = expandfirst,
            doctest = false,
            format = Documenter.HTML(
                assets = [
                    "assets/favicon.ico",
                    "assets/custom.css",
                    asset("https://example.com/resource.js"),
                    asset("http://example.com/fonts?param=foo", class=:css),
                    asset("https://fonts.googleapis.com/css?family=Nanum+Brush+Script&display=swap", class=:css),
                ],
                prettyurls = true,
                canonical = "https://example.com/stable",
                mathengine = mathengine,
                highlights = ["erlang", "erlang-repl"],
                footer = "This footer has been customized.",
            )
        )
    end
end

# Build with pretty URLs and canonical links and a PNG logo
examples_html_doc = if "html" in EXAMPLE_BUILDS
    @info("Building mock package docs: HTMLWriter / deployment build")
    html_doc("html",
        MathJax2(Dict(
            :TeX => Dict(
                :equationNumbers => Dict(:autoNumber => "AMS"),
                :Macros => Dict(
                    :ket => ["|#1\\rangle", 1],
                    :bra => ["\\langle#1|", 1],
                    :pdv => ["\\frac{\\partial^{#1} #2}{\\partial #3^{#1}}", 3, ""],
                ),
            ),
        )),
    )
else
    @info "Skipping build: HTML/deploy" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

# same as HTML, but with MathJax3
examples_html_mathjax3_doc = if "html-mathjax3" in EXAMPLE_BUILDS
    @info("Building mock package docs: HTMLWriter / deployment build using MathJax v3")
    html_doc("html-mathjax3",
        MathJax3(Dict(
            :loader => Dict("load" => ["[tex]/physics"]),
            :tex => Dict(
                "inlineMath" => [["\$","\$"], ["\\(","\\)"]],
                "tags" => "ams",
                "packages" => ["base", "ams", "autoload", "physics"],
            ),
        )),
    )
else
    @info "Skipping build: HTML/deploy MathJax v3" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

# HTML: local build with pretty URLs off
examples_html_local_doc = if "html-local" in EXAMPLE_BUILDS
    @info("Building mock package docs: HTMLWriter / local build")
    @quietly makedocs(
        debug = true,
        root  = examples_root,
        build = "builds/html-local",
        doctestfilters = [r"Ptr{0x[0-9]+}"],
        sitename = "Documenter example",
        pages = htmlbuild_pages,
        expandfirst = expandfirst,

        linkcheck = true,
        linkcheck_ignore = [r"(x|y).md", "z.md", r":func:.*"],
        format = Documenter.HTML(
            assets = ["assets/custom.css"],
            prettyurls = false,
            edit_branch = nothing,
            footer = nothing,
        ),
    )
else
    @info "Skipping build: HTML/local" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

# Markdown
examples_markdown_doc = if "markdown" in EXAMPLE_BUILDS
    @info("Building mock package docs: MarkdownWriter")
    @quietly makedocs(
        format = Markdown(),
        debug = true,
        root  = examples_root,
        build = "builds/markdown",
        doctest = false,
        expandfirst = expandfirst,
    )
else
    @info "Skipping build: Markdown" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

# PDF/LaTeX
examples_latex_simple_doc = if "latex_simple" in EXAMPLE_BUILDS
    @info("Building mock package docs: LaTeXWriter/simple")
    @quietly makedocs(
        format = Documenter.Writers.LaTeXWriter.LaTeX(platform = "docker"),
        sitename = "Documenter LaTeX Simple",
        root  = examples_root,
        build = "builds/latex_simple",
        source = "src.latex_simple",
        pages = ["Main section" => ["index.md"]],
        doctest = false,
        debug = true,
    )
else
    @info "Skipping build: LaTeXWriter/simple" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

examples_latex_doc = if "latex" in EXAMPLE_BUILDS
    @info("Building mock package docs: LaTeXWriter/latex")
    @quietly makedocs(
        format = Documenter.Writers.LaTeXWriter.LaTeX(platform = "docker"),
        sitename = "Documenter LaTeX",
        root  = examples_root,
        build = "builds/latex",
        pages = htmlbuild_pages = Any[
            "General" => [
                "index.md",
                "latex.md",
                "unicode.md",
                hide("hidden.md"),
                "example-output.md",
            ],
            # SVG images nor code blocks in footnotes are allowed in LaTeX
            # "Manual" => [
            #     "man/tutorial.md",
            #     "man/style.md",
            # ],
            hide("Hidden Pages" => "hidden/index.md", Any[
                "Page X" => "hidden/x.md",
                "hidden/y.md",
                "hidden/z.md",
            ]),
            "Library" => [
                "lib/functions.md",
                "lib/autodocs.md",
                "lib/editurl.md",
            ],
            "Expandorder" => [
                "expandorder/00.md",
                "expandorder/01.md",
                "expandorder/AA.md",
            ],
        ],
        doctest = false,
        debug = true,
    )
else
    @info "Skipping build: LaTeXWriter/latex" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

examples_latex_simple_nondocker_doc = if "latex_simple_nondocker" in EXAMPLE_BUILDS
    @info("Building mock package docs: LaTeXWriter/latex_simple_nondocker")
    @quietly makedocs(
        format = Documenter.Writers.LaTeXWriter.LaTeX(),
        sitename = "Documenter LaTeX Simple Non-Docker",
        root  = examples_root,
        build = "builds/latex_simple_nondocker",
        source = "src.latex_simple",
        pages = ["Main section" => ["index.md"]],
        doctest = false,
        debug = true,
    )
else
    @info "Skipping build: LaTeXWriter/latex_simple_nondocker" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end

examples_latex_texonly_doc = if "latex_texonly" in EXAMPLE_BUILDS
    @info("Building mock package docs: LaTeXWriter/latex_texonly")
    @quietly makedocs(
        format = Documenter.Writers.LaTeXWriter.LaTeX(platform = "none"),
        sitename = "Documenter LaTeX",
        root  = examples_root,
        build = "builds/latex_texonly",
        pages = htmlbuild_pages = Any[
            "General" => [
                "index.md",
                "latex.md",
                "unicode.md",
                hide("hidden.md"),
                "example-output.md",
            ],
            # SVG images nor code blocks in footnotes are allowed in LaTeX
            # "Manual" => [
            #     "man/tutorial.md",
            #     "man/style.md",
            # ],
            hide("Hidden Pages" => "hidden/index.md", Any[
                "Page X" => "hidden/x.md",
                "hidden/y.md",
                "hidden/z.md",
            ]),
            "Library" => [
                "lib/functions.md",
                "lib/autodocs.md",
                "lib/editurl.md",
            ],
            "Expandorder" => [
                "expandorder/00.md",
                "expandorder/01.md",
                "expandorder/AA.md",
            ],
        ],
        doctest = false,
        debug = true,
    )
else
    @info "Skipping build: LaTeXWriter/latex_texonly" EXAMPLE_BUILDS get(ENV, "DOCUMENTER_TEST_EXAMPLES", nothing)
    nothing
end


text:
#Demonstrate the designLCBP function
#=WARNING:
designLCBP is very ill-behaved due to deficiencies in MATLAB's constr() and
minimax(). These functions frequently take a viable set of parameters as a
starting point and turn them into an unstable modulator. I should provide a
more robust implementation...
=#
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import Printf: @sprintf
j=im

#==
===============================================================================#
println("\n*** Continuous-Time LC Modulator Design")

n = 3
OSR = 64
opt = 2
Hinf = 1.7
f0 = 1/16
t = [0.5 1]
form = :FB
dbg = true
(param, H, L0, ABCD) = designLCBP(n, OSR=OSR, opt=opt,
	Hinf=Hinf, f0=f0, t=t, form=form, dbg=dbg
)

:END_OF_DEMO


text:
using GR

const viridis = [0x440154, 0x472c7a, 0x3b518b, 0x2c718e, 0x21908d, 0x27ad81, 0x5cc863, 0xaadc32, 0xfde725]
const inferno = [0x000004, 0x1f0c48, 0x550f6d, 0x88226a, 0xa83655, 0xe35933, 0xf9950a, 0xf8c932, 0xfcffa4]
const plasma = [0x0c0887, 0x4b03a1, 0x7d03a8, 0xa82296, 0xcb4679, 0xe56b5d, 0xf89441, 0xfdc328, 0xf0f921]
const magma = [0x000004, 0x1c1044, 0x4f127b, 0x812581, 0xb5367a, 0xe55964, 0xfb8761, 0xfec287, 0xfbfdbf]

red(a)   = ((a .>> 16) .& 0xff ) ./ 256
green(a) = ((a .>>  8) .& 0xff ) ./ 256
blue(a)  = ( a         .& 0xff ) ./ 256

function show_colormaps()
    setwindow(1, 5, 0, 1)
    for (x, cmap) in enumerate((viridis, inferno, plasma, magma))
        setcolormapfromrgb(red(cmap), green(cmap), blue(cmap))
        cellarray(x, x + 0.25, 0, 1, 1, 255, 1000:1255)
        setcolormap(GR.COLORMAP_VIRIDIS + x - 1)
        cellarray(x + 0.5, x + 0.75, 0, 1, 1, 255, 1000:1255)
    end
    updatews()
end

show_colormaps()


text:
# Autogenerated wrapper script for SoXResampler_jll for aarch64-linux-musl
export libsoxr, libsoxr_lsr

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SoXResampler")
JLLWrappers.@declare_library_product(libsoxr, "libsoxr.so.0")
JLLWrappers.@declare_library_product(libsoxr_lsr, "libsoxr-lsr.so.0")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libsoxr,
        "lib/libsoxr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsoxr_lsr,
        "lib/libsoxr-lsr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function _get_minutes_per_period(container::OptimizationContainer)
    resolution = get_resolution(container)
    if resolution > Dates.Minute(1)
        minutes_per_period = Dates.value(Dates.Minute(resolution))
    else
        @warn("Not all formulations support under 1-minute resolutions. Exercise caution.")
        minutes_per_period = Dates.value(Dates.Second(resolution)) / 60
    end
    return minutes_per_period
end

function _get_ramp_constraint_devices(
    container::OptimizationContainer,
    devices::Union{Vector{U}, IS.FlattenIteratorWrapper{U}},
) where {U <: PSY.Component}
    minutes_per_period = _get_minutes_per_period(container)
    filtered_device = Vector{U}()
    for d in devices
        ramp_limits = PSY.get_ramp_limits(d)
        if ramp_limits !== nothing
            p_lims = PSY.get_active_power_limits(d)
            max_rate = abs(p_lims.min - p_lims.max) / minutes_per_period
            if (ramp_limits.up >= max_rate) & (ramp_limits.down >= max_rate)
                @debug "Generator has a nonbinding ramp limits. Constraints Skipped" PSY.get_name(
                    d,
                )
                continue
            else
                push!(filtered_device, d)
            end
        end
    end
    return filtered_device
end

@doc raw"""
Constructs allowed rate-of-change constraints from variables, initial condtions, and rate data.

If t = 1:

``` variable[name, 1] - initial_conditions[ix].value <= rate_data[1][ix].up ```

``` initial_conditions[ix].value - variable[name, 1] <= rate_data[1][ix].down ```

If t > 1:

``` variable[name, t] - variable[name, t-1] <= rate_data[1][ix].up ```

``` variable[name, t-1] - variable[name, t] <= rate_data[1][ix].down ```

# LaTeX

`` r^{down} \leq x_1 - x_{init} \leq r^{up}, \text{ for } t = 1 ``

`` r^{down} \leq x_t - x_{t-1} \leq r^{up}, \forall t \geq 2 ``

"""
function add_linear_ramp_constraints!(
    container::OptimizationContainer,
    T::Type{<:ConstraintType},
    U::Type{S},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {
    S <: Union{PowerAboveMinimumVariable, ActivePowerVariable},
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    parameters = built_for_recurrent_solves(container)
    time_steps = get_time_steps(container)
    variable = get_variable(container, U(), V)
    ramp_devices = _get_ramp_constraint_devices(container, devices)
    minutes_per_period = _get_minutes_per_period(container)
    IC = _get_initial_condition_type(T, V, W)
    initial_conditions_power = get_initial_condition(container, IC(), V)
    expr_dn = get_expression(container, ActivePowerRangeExpressionLB(), V)
    expr_up = get_expression(container, ActivePowerRangeExpressionUB(), V)

    set_name = [PSY.get_name(r) for r in ramp_devices]
    con_up = add_constraints_container!(container, T(), V, set_name, time_steps, meta="up")
    con_down =
        add_constraints_container!(container, T(), V, set_name, time_steps, meta="dn")

    for ic in initial_conditions_power
        name = get_component_name(ic)
        # This is to filter out devices that dont need a ramping constraint
        name ∉ set_name && continue
        ramp_limits = PSY.get_ramp_limits(get_component(ic))
        ic_power = get_value(ic)
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        con_up[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            expr_up[name, 1] - ic_power <= ramp_limits.up * minutes_per_period
        )
        con_down[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            ic_power - expr_dn[name, 1] >= -1 * ramp_limits.down * minutes_per_period
        )
        for t in time_steps[2:end]
            con_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                expr_up[name, t] - variable[name, t - 1] <=
                ramp_limits.up * minutes_per_period
            )
            con_down[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t - 1] - expr_dn[name, t] >=
                -1 * ramp_limits.down * minutes_per_period
            )
        end
    end
    return
end

function add_linear_ramp_constraints!(
    container::OptimizationContainer,
    T::Type{<:ConstraintType},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.Component, W <: AbstractDeviceFormulation}
    parameters = built_for_recurrent_solves(container)
    time_steps = get_time_steps(container)
    variable = get_variable(container, U(), V)
    ramp_devices = _get_ramp_constraint_devices(container, devices)
    minutes_per_period = _get_minutes_per_period(container)
    IC = _get_initial_condition_type(T, V, W)
    initial_conditions_power = get_initial_condition(container, IC(), V)

    set_name = [PSY.get_name(r) for r in ramp_devices]
    con_up = add_constraints_container!(container, T(), V, set_name, time_steps, meta="up")
    con_down =
        add_constraints_container!(container, T(), V, set_name, time_steps, meta="dn")

    for ic in initial_conditions_power
        name = get_component_name(ic)
        # This is to filter out devices that dont need a ramping constraint
        name ∉ set_name && continue
        ramp_limits = PSY.get_ramp_limits(get_component(ic))
        ic_power = get_value(ic)
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        con_up[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            variable[name, 1] - ic_power <= ramp_limits.up * minutes_per_period
        )
        con_down[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            ic_power - variable[name, 1] <= ramp_limits.down * minutes_per_period
        )
        for t in time_steps[2:end]
            con_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t] - variable[name, t - 1] <=
                ramp_limits.up * minutes_per_period
            )
            con_down[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t - 1] - variable[name, t] <=
                ramp_limits.down * minutes_per_period
            )
        end
    end
    return
end

@doc raw"""
Constructs allowed rate-of-change constraints from variables, initial condtions, start/stop status, and rate data

# Equations
If t = 1:

``` variable[name, 1] - initial_conditions[ix].value <= rate_data[1][ix].up + rate_data[2][ix].max*varstart[name, 1] ```

``` initial_conditions[ix].value - variable[name, 1] <= rate_data[1][ix].down + rate_data[2][ix].min*varstop[name, 1] ```

If t > 1:

``` variable[name, t] - variable[name, t-1] <= rate_data[1][ix].up + rate_data[2][ix].max*varstart[name, t] ```

``` variable[name, t-1] - variable[name, t] <= rate_data[1][ix].down + rate_data[2][ix].min*varstop[name, t] ```

# LaTeX

`` r^{down} + r^{min} x^{stop}_1 \leq x_1 - x_{init} \leq r^{up} + r^{max} x^{start}_1, \text{ for } t = 1 ``

`` r^{down} + r^{min} x^{stop}_t \leq x_t - x_{t-1} \leq r^{up} + r^{max} x^{start}_t, \forall t \geq 2 ``
"""
function add_semicontinuous_ramp_constraints!(
    container::OptimizationContainer,
    T::Type{<:ConstraintType},
    U::Type{S},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {
    S <: Union{PowerAboveMinimumVariable, ActivePowerVariable},
    V <: PSY.Component,
    W <: AbstractDeviceFormulation,
}
    parameters = built_for_recurrent_solves(container)
    time_steps = get_time_steps(container)
    variable = get_variable(container, U(), V)
    varstart = get_variable(container, StartVariable(), V)
    varstop = get_variable(container, StopVariable(), V)

    ramp_devices = _get_ramp_constraint_devices(container, devices)
    minutes_per_period = _get_minutes_per_period(container)
    IC = _get_initial_condition_type(T, V, W)
    initial_conditions_power = get_initial_condition(container, IC(), V)
    expr_dn = get_expression(container, ActivePowerRangeExpressionLB(), V)
    expr_up = get_expression(container, ActivePowerRangeExpressionUB(), V)

    set_name = [PSY.get_name(r) for r in ramp_devices]
    con_up = add_constraints_container!(container, T(), V, set_name, time_steps, meta="up")
    con_down =
        add_constraints_container!(container, T(), V, set_name, time_steps, meta="dn")

    for ic in initial_conditions_power
        name = get_component_name(ic)
        # This is to filter out devices that dont need a ramping constraint
        name ∉ set_name && continue
        device = get_component(ic)
        ramp_limits = PSY.get_ramp_limits(device)
        power_limits = PSY.get_active_power_limits(device)
        ic_power = get_value(ic)
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        con_up[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            expr_up[name, 1] - ic_power <=
            ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, 1]
        )
        con_down[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            ic_power - expr_dn[name, 1] <=
            ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, 1]
        )
        for t in time_steps[2:end]
            con_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                expr_up[name, t] - variable[name, t - 1] <=
                ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, t]
            )
            con_down[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t - 1] - expr_dn[name, t] <=
                ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, t]
            )
        end
    end
    return
end

function add_semicontinuous_ramp_constraints!(
    container::OptimizationContainer,
    T::Type{<:ConstraintType},
    U::Type{<:VariableType},
    devices::IS.FlattenIteratorWrapper{V},
    model::DeviceModel{V, W},
    X::Type{<:PM.AbstractPowerModel},
) where {V <: PSY.Component, W <: AbstractDeviceFormulation}
    parameters = built_for_recurrent_solves(container)
    time_steps = get_time_steps(container)
    variable = get_variable(container, U(), V)
    varstart = get_variable(container, StartVariable(), V)
    varstop = get_variable(container, StopVariable(), V)

    ramp_devices = _get_ramp_constraint_devices(container, devices)
    minutes_per_period = _get_minutes_per_period(container)
    IC = _get_initial_condition_type(T, V, W)
    initial_conditions_power = get_initial_condition(container, IC(), V)

    set_name = [PSY.get_name(r) for r in ramp_devices]
    con_up = add_constraints_container!(container, T(), V, set_name, time_steps, meta="up")
    con_down =
        add_constraints_container!(container, T(), V, set_name, time_steps, meta="dn")

    for ic in initial_conditions_power
        name = get_component_name(ic)
        # This is to filter out devices that dont need a ramping constraint
        name ∉ set_name && continue
        device = get_component(ic)
        ramp_limits = PSY.get_ramp_limits(device)
        power_limits = PSY.get_active_power_limits(device)
        ic_power = get_value(ic)
        @debug "add rate_of_change_constraint" name ic_power
        @assert (parameters && isa(ic_power, PJ.ParameterRef)) || !parameters
        con_up[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            variable[name, 1] - ic_power <=
            ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, 1]
        )
        con_down[name, 1] = JuMP.@constraint(
            container.JuMPmodel,
            ic_power - variable[name, 1] <=
            ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, 1]
        )
        for t in time_steps[2:end]
            con_up[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t] - variable[name, t - 1] <=
                ramp_limits.up * minutes_per_period + power_limits.min * varstart[name, t]
            )
            con_down[name, t] = JuMP.@constraint(
                container.JuMPmodel,
                variable[name, t - 1] - variable[name, t] <=
                ramp_limits.down * minutes_per_period + power_limits.min * varstop[name, t]
            )
        end
    end
    return
end


text:
# unit test the bGrad

include("../src/DMD_Util.jl");
include("../src/DMD_Type.jl");
include("../src/DMD_Func.jl");

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, X, t, lossFunc, lossGrad);
#
# Initialize with true parameter
#------------------------------------------------------------------------------
copyto!(params.a, at);
copyto!(params.B, Bt);
#
# Evaluate bGrad
#------------------------------------------------------------------------------
update_P!(params);
id = rand(1:n);
#
bt = copy(params.br[id]);
Random.randn!(params.br[id]);
#
projb(params, id);
#
err = sum(abs2, bt - params.br[id]);
#
if err < 1e-6
	println("projb: OK");
else
	println("projb: Wrong, err: $err");
end


text:
using MAT
cd(dirname(Base.source_path()))
include("../src/utils.jl")
include("../src/1d/cpab1dIntegration.jl")
file = matopen("cpa_data_nCells_5.mat")
As = read(file, "As")
pts_src = read(file, "pts_src")
pts_src = squeeze(pts_src,2)
tw_args = read(file, "tw_args")
nCols = tw_args["nCols"]
nCells = size(As, 1)
close(file)

nSteps = 100
Nsteps = 10
tess = getTess(As, nCells, nCols)
v = getVelocityField1d(pts_src, As, nCells, nCols)

pts_fwd1 = solver1dParallelIntegration(pts_src, 0.1, tess, As, nCells, nCols, Nsteps, nSteps)
pts_fwd2 = solver1dParallelIntegration(pts_src, 0.5, tess, As, nCells, nCols, Nsteps, nSteps)
@time pts_fwd3 = solver1dParallelIntegration(pts_src, 1.0, tess, As, nCells, nCols, Nsteps, nSteps)
pts_fwd4 = solver1dParallelIntegration(pts_src, 2.0, tess, As, nCells, nCols, Nsteps, nSteps)
pts_fwd5 = solver1dParallelIntegration(pts_src, 8.0, tess, As, nCells, nCols, Nsteps, nSteps)

using PyPlot
PyPlot.figure("1D Integration",figsize=(12,6))
PyPlot.subplot(1,2,1)
PyPlot.title("Velocity Field")
PyPlot.plot(pts_src, v, color="blue")
PyPlot.xlabel("x")
PyPlot.ylabel("v")

PyPlot.subplot(1,2,2)
PyPlot.title("1D Integration for different t")
PyPlot.plot(pts_src, pts_src, label="t=0", color="black")
PyPlot.plot(pts_src, pts_fwd1, label="t=0.1", color="red")
PyPlot.plot(pts_src, pts_fwd2, label="t=0.5", color="orange")
PyPlot.plot(pts_src, pts_fwd3, label="t=1.0", color="yellow")
PyPlot.plot(pts_src, pts_fwd4, label="t=2.0", color="green")
PyPlot.plot(pts_src, pts_fwd5, label="t=8.0", color="blue")
PyPlot.xlabel("x")
PyPlot.ylabel("Transformed x")
PyPlot.legend(loc="lower right",fancybox="true")


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct RLFilter <: Filter
        rf::Float64
        lf::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
    end

Parameters of RL series filter in algebraic representation

# Arguments
- `rf::Float64`: Series resistance in p.u. of converter filter to the grid, validation range: `(0, nothing)`
- `lf::Float64`: Series inductance in p.u. of converter filter to the grid, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: RLFilter has zero states
- `n_states::Int`: RLFilter has zero states
"""
mutable struct RLFilter <: Filter
    "Series resistance in p.u. of converter filter to the grid"
    rf::Float64
    "Series inductance in p.u. of converter filter to the grid"
    lf::Float64
    ext::Dict{String, Any}
    "RLFilter has zero states"
    states::Vector{Symbol}
    "RLFilter has zero states"
    n_states::Int
end

function RLFilter(rf, lf, ext=Dict{String, Any}(), )
    RLFilter(rf, lf, ext, Vector{Symbol}(), 0, )
end

function RLFilter(; rf, lf, ext=Dict{String, Any}(), states=Vector{Symbol}(), n_states=0, )
    RLFilter(rf, lf, ext, states, n_states, )
end

# Constructor for demo purposes; non-functional.
function RLFilter(::Nothing)
    RLFilter(;
        rf=0,
        lf=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`RLFilter`](@ref) `rf`."""
get_rf(value::RLFilter) = value.rf
"""Get [`RLFilter`](@ref) `lf`."""
get_lf(value::RLFilter) = value.lf
"""Get [`RLFilter`](@ref) `ext`."""
get_ext(value::RLFilter) = value.ext
"""Get [`RLFilter`](@ref) `states`."""
get_states(value::RLFilter) = value.states
"""Get [`RLFilter`](@ref) `n_states`."""
get_n_states(value::RLFilter) = value.n_states

"""Set [`RLFilter`](@ref) `rf`."""
set_rf!(value::RLFilter, val) = value.rf = val
"""Set [`RLFilter`](@ref) `lf`."""
set_lf!(value::RLFilter, val) = value.lf = val
"""Set [`RLFilter`](@ref) `ext`."""
set_ext!(value::RLFilter, val) = value.ext = val


text:
# Julia wrapper for header: include/declarations.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0 and then modified as some places

export read_prob

function triu(A::blockmatrix)
    ccall((:triu,CSDP.csdp),Void,(blockmatrix,),A)
end

function store_packed(A::blockmatrix,B::blockmatrix)
    ccall((:store_packed,CSDP.csdp),Void,(blockmatrix,blockmatrix),A,B)
end

function store_unpacked(A::blockmatrix,B::blockmatrix)
    ccall((:store_unpacked,CSDP.csdp),Void,(blockmatrix,blockmatrix),A,B)
end

function alloc_mat_packed(A::blockmatrix,pB::Ptr{blockmatrix})
    ccall((:alloc_mat_packed,CSDP.csdp),Void,(blockmatrix,Ptr{blockmatrix}),A,pB)
end

function free_mat_packed(A::blockmatrix)
    ccall((:free_mat_packed,CSDP.csdp),Void,(blockmatrix,),A)
end

function structnnz(n::Cint,k::Cint,C::blockmatrix,constraints::Ptr{constraintmatrix})
    ccall((:structnnz,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{constraintmatrix}),n,k,C,constraints)
end

function actnnz(n::Cint,lda::Cint,A::Ptr{Cdouble})
    ccall((:actnnz,CSDP.csdp),Cint,(Cint,Cint,Ptr{Cdouble}),n,lda,A)
end

function bandwidth(n::Cint,lda::Cint,A::Ptr{Cdouble})
    ccall((:bandwidth,CSDP.csdp),Cint,(Cint,Cint,Ptr{Cdouble}),n,lda,A)
end

function qreig(n::Cint,maindiag::Ptr{Cdouble},offdiag::Ptr{Cdouble})
    ccall((:qreig,CSDP.csdp),Void,(Cint,Ptr{Cdouble},Ptr{Cdouble}),n,maindiag,offdiag)
end

function sort_entries(k::Cint,C::blockmatrix,constraints::Ptr{constraintmatrix})
    ccall((:sort_entries,CSDP.csdp),Void,(Cint,blockmatrix,Ptr{constraintmatrix}),k,C,constraints)
end

function norm2(n::Cint,x::Ptr{Cdouble})
    ccall((:norm2,CSDP.csdp),Cdouble,(Cint,Ptr{Cdouble}),n,x)
end

function norm1(n::Cint,x::Ptr{Cdouble})
    ccall((:norm1,CSDP.csdp),Cdouble,(Cint,Ptr{Cdouble}),n,x)
end

function norminf(n::Cint,x::Ptr{Cdouble})
    ccall((:norminf,CSDP.csdp),Cdouble,(Cint,Ptr{Cdouble}),n,x)
end

function Fnorm(A::blockmatrix)
    ccall((:Fnorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function Knorm(A::blockmatrix)
    ccall((:Knorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function mat1norm(A::blockmatrix)
    ccall((:mat1norm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function matinfnorm(A::blockmatrix)
    ccall((:matinfnorm,CSDP.csdp),Cdouble,(blockmatrix,),A)
end

function calc_pobj(C::blockmatrix,X::blockmatrix,constant_offset::Cdouble)
    ccall((:calc_pobj,CSDP.csdp),Cdouble,(blockmatrix,blockmatrix,Cdouble),C,X,constant_offset)
end

function calc_dobj(k::Cint,a::Ptr{Cdouble},y::Ptr{Cdouble},constant_offset::Cdouble)
    ccall((:calc_dobj,CSDP.csdp),Cdouble,(Cint,Ptr{Cdouble},Ptr{Cdouble},Cdouble),k,a,y,constant_offset)
end

function trace_prod(A::blockmatrix,B::blockmatrix)
    ccall((:trace_prod,CSDP.csdp),Cdouble,(blockmatrix,blockmatrix),A,B)
end

function linesearch(n::Cint,dX::blockmatrix,work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,cholinv::blockmatrix,q::Ptr{Cdouble},z::Ptr{Cdouble},workvec::Ptr{Cdouble},stepfrac::Cdouble,start::Cdouble,printlevel::Cint)
    ccall((:linesearch,CSDP.csdp),Cdouble,(Cint,blockmatrix,blockmatrix,blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cdouble,Cdouble,Cint),n,dX,work1,work2,work3,cholinv,q,z,workvec,stepfrac,start,printlevel)
end

function pinfeas(k::Cint,constraints::Ptr{constraintmatrix},X::blockmatrix,a::Ptr{Cdouble},workvec::Ptr{Cdouble})
    ccall((:pinfeas,CSDP.csdp),Cdouble,(Cint,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},Ptr{Cdouble}),k,constraints,X,a,workvec)
end

function dinfeas(k::Cint,C::blockmatrix,constraints::Ptr{constraintmatrix},y::Ptr{Cdouble},Z::blockmatrix,work1::blockmatrix)
    ccall((:dinfeas,CSDP.csdp),Cdouble,(Cint,blockmatrix,Ptr{constraintmatrix},Ptr{Cdouble},blockmatrix,blockmatrix),k,C,constraints,y,Z,work1)
end

function dimacserr3(k::Cint,C::blockmatrix,constraints::Ptr{constraintmatrix},y::Ptr{Cdouble},Z::blockmatrix,work1::blockmatrix)
    ccall((:dimacserr3,CSDP.csdp),Cdouble,(Cint,blockmatrix,Ptr{constraintmatrix},Ptr{Cdouble},blockmatrix,blockmatrix),k,C,constraints,y,Z,work1)
end

function op_a(k::Cint,constraints::Ptr{constraintmatrix},X::blockmatrix,result::Ptr{Cdouble})
    ccall((:op_a,CSDP.csdp),Void,(Cint,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble}),k,constraints,X,result)
end

function op_at(k::Cint,y::Ptr{Cdouble},constraints::Ptr{constraintmatrix},result::blockmatrix)
    ccall((:op_at,CSDP.csdp),Void,(Cint,Ptr{Cdouble},Ptr{constraintmatrix},blockmatrix),k,y,constraints,result)
end

function makefill(k::Cint,C::blockmatrix,constraints::Ptr{constraintmatrix},pfill::Ptr{constraintmatrix},work1::blockmatrix,printlevel::Cint)
    ccall((:makefill,CSDP.csdp),Void,(Cint,blockmatrix,Ptr{constraintmatrix},Ptr{constraintmatrix},blockmatrix,Cint),k,C,constraints,pfill,work1,printlevel)
end

function op_o(k::Cint,constraints::Ptr{constraintmatrix},byblocks::Ptr{Ptr{sparseblock}},Zi::blockmatrix,X::blockmatrix,O::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix)
    ccall((:op_o,CSDP.csdp),Void,(Cint,Ptr{constraintmatrix},Ptr{Ptr{sparseblock}},blockmatrix,blockmatrix,Ptr{Cdouble},blockmatrix,blockmatrix),k,constraints,byblocks,Zi,X,O,work1,work2)
end

function addscaledmat(A::blockmatrix,scale::Cdouble,B::blockmatrix,C::blockmatrix)
    ccall((:addscaledmat,CSDP.csdp),Void,(blockmatrix,Cdouble,blockmatrix,blockmatrix),A,scale,B,C)
end

function zero_mat(A::blockmatrix)
    ccall((:zero_mat,CSDP.csdp),Void,(blockmatrix,),A)
end

function add_mat(A::blockmatrix,B::blockmatrix)
    ccall((:add_mat,CSDP.csdp),Void,(blockmatrix,blockmatrix),A,B)
end

function sym_mat(A::blockmatrix)
    ccall((:sym_mat,CSDP.csdp),Void,(blockmatrix,),A)
end

function make_i(A::blockmatrix)
    ccall((:make_i,CSDP.csdp),Void,(blockmatrix,),A)
end

function copy_mat(A::blockmatrix,B::blockmatrix)
    ccall((:copy_mat,CSDP.csdp),Void,(blockmatrix,blockmatrix),A,B)
end

function mat_mult(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix)
    ccall((:mat_mult,CSDP.csdp),Void,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix),scale1,scale2,A,B,C)
end

function mat_multspa(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspa,CSDP.csdp),Void,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_multspb(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspb,CSDP.csdp),Void,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_multspc(scale1::Cdouble,scale2::Cdouble,A::blockmatrix,B::blockmatrix,C::blockmatrix,fill::constraintmatrix)
    ccall((:mat_multspc,CSDP.csdp),Void,(Cdouble,Cdouble,blockmatrix,blockmatrix,blockmatrix,constraintmatrix),scale1,scale2,A,B,C,fill)
end

function mat_mult_raw(n::Cint,scale1::Cdouble,scale2::Cdouble,ap::Ptr{Cdouble},bp::Ptr{Cdouble},cp::Ptr{Cdouble})
    ccall((:mat_mult_raw,CSDP.csdp),Void,(Cint,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),n,scale1,scale2,ap,bp,cp)
end

function mat_mult_rawatlas(n::Cint,scale1::Cdouble,scale2::Cdouble,ap::Ptr{Cdouble},bp::Ptr{Cdouble},cp::Ptr{Cdouble})
    ccall((:mat_mult_rawatlas,CSDP.csdp),Void,(Cint,Cdouble,Cdouble,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble}),n,scale1,scale2,ap,bp,cp)
end

function matvec(A::blockmatrix,x::Ptr{Cdouble},y::Ptr{Cdouble})
    ccall((:matvec,CSDP.csdp),Void,(blockmatrix,Ptr{Cdouble},Ptr{Cdouble}),A,x,y)
end

function alloc_mat(A::blockmatrix,pB::Ptr{blockmatrix})
    ccall((:alloc_mat,CSDP.csdp),Void,(blockmatrix,Ptr{blockmatrix}),A,pB)
end

function free_mat(A::blockmatrix)
    ccall((:free_mat,CSDP.csdp),Void,(blockmatrix,),A)
end

function initparams(params::Ptr{paramstruc},pprintlevel::Ptr{Cint})
    ccall((:initparams,CSDP.csdp),Void,(Ptr{paramstruc},Ptr{Cint}),params,pprintlevel)
end

function initsoln(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix})
    X = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    y = Ref{Ptr{Cdouble}}(C_NULL)
    Z = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    ccall((:initsoln,CSDP.csdp),Void,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{blockmatrix}),n,k,C,a,constraints,X,y,Z)
    X[], y[], Z[]
end

function trans(A::blockmatrix)
    ccall((:trans,CSDP.csdp),Void,(blockmatrix,),A)
end

function chol_inv(A::blockmatrix,B::blockmatrix)
    ccall((:chol_inv,CSDP.csdp),Void,(blockmatrix,blockmatrix),A,B)
end

function chol(A::blockmatrix)
    ccall((:chol,CSDP.csdp),Cint,(blockmatrix,),A)
end

function solvesys(m::Cint,ldam::Cint,A::Ptr{Cdouble},rhs::Ptr{Cdouble})
    ccall((:solvesys,CSDP.csdp),Cint,(Cint,Cint,Ptr{Cdouble},Ptr{Cdouble}),m,ldam,A,rhs)
end

function user_exit(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},dobj::Cdouble,pobj::Cdouble,constant_offset::Cdouble,constraints::Ptr{constraintmatrix},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix,params::paramstruc)
    ccall((:user_exit,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{Cdouble},Cdouble,Cdouble,Cdouble,Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},blockmatrix,paramstruc),n,k,C,a,dobj,pobj,constant_offset,constraints,X,y,Z,params)
end

function read_sol(fname::Ptr{UInt8},n::Cint,k::Cint,C::blockmatrix,pX::Ptr{blockmatrix},py::Ptr{Ptr{Cdouble}},pZ::Ptr{blockmatrix})
    ccall((:read_sol,CSDP.csdp),Cint,(Ptr{UInt8},Cint,Cint,blockmatrix,Ptr{blockmatrix},Ptr{Ptr{Cdouble}},Ptr{blockmatrix}),fname,n,k,C,pX,py,pZ)
end

function read_prob(fname::String,printlevel::Integer=0)
    n = Ref{Cint}(0)
    k = Ref{Cint}(0)
    C = Ref{blockmatrix}(blockmatrix(0, C_NULL))
    a = Ref{Ptr{Cdouble}}(C_NULL)
    constraints = Ref{Ptr{constraintmatrix}}(C_NULL)
    ccall((:read_prob,CSDP.csdp),Cint,(Ptr{UInt8},Ref{Cint},Ref{Cint},Ref{blockmatrix},Ref{Ptr{Cdouble}},Ref{Ptr{constraintmatrix}},Cint),fname,n,k,C,a,constraints,printlevel)
    C = mywrap(C[])
    a = mywrap(a[], k[])
    constraints = mywrap(constraints[], k[])
    # Nothing to free as the constraintmatrix is immutable
    # the array is allocated as an array of bitstype in CSDP
    # the array of blocks, .. inside should be free'd though
    #constraints = map(c->ConstraintMatrix(c, k[]), constraints)
    @assert n[] == size(C, 1)
    C, a, constraints
end

function write_prob(fname::String,
                    n::Cint,
                    k::Cint,
                    C::blockmatrix,
                    a::Ptr{Cdouble},
                    constraints::Ptr{constraintmatrix})
    ccall((:write_prob,CSDP.csdp),Cint,(Ptr{UInt8},Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix}),
          fname,n,k,C,a,constraints)
end

function write_sol(fname::String,
                   n::Cint,
                   k::Cint,
                   X::blockmatrix,
                   y::Ptr{Cdouble},
                   Z::blockmatrix)
    ccall((:write_sol,CSDP.csdp),
          Cint,
          (Ptr{UInt8},
           Cint,
           Cint,
           blockmatrix,
           Ptr{Cdouble},
           blockmatrix),
          fname,n,k,X,y,Z)
end

function free_prob(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix)
    ccall((:free_prob,CSDP.csdp),Void,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},blockmatrix,Ptr{Cdouble},blockmatrix),n,k,C,a,constraints,X,y,Z)
end

function sdp(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constant_offset::Cdouble,constraints::Ptr{constraintmatrix},byblocks::Ptr{Ptr{sparseblock}},fill::constraintmatrix,X::blockmatrix,y::Ptr{Cdouble},Z::blockmatrix,cholxinv::blockmatrix,cholzinv::blockmatrix,pobj::Ptr{Cdouble},dobj::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,workvec1::Ptr{Cdouble},workvec2::Ptr{Cdouble},workvec3::Ptr{Cdouble},workvec4::Ptr{Cdouble},workvec5::Ptr{Cdouble},workvec6::Ptr{Cdouble},workvec7::Ptr{Cdouble},workvec8::Ptr{Cdouble},diagO::Ptr{Cdouble},bestx::blockmatrix,besty::Ptr{Cdouble},bestz::blockmatrix,Zi::blockmatrix,O::Ptr{Cdouble},rhs::Ptr{Cdouble},dZ::blockmatrix,dX::blockmatrix,dy::Ptr{Cdouble},dy1::Ptr{Cdouble},Fp::Ptr{Cdouble},printlevel::Cint,parameters::paramstruc)
    ccall((:sdp,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{Cdouble},Cdouble,Ptr{constraintmatrix},Ptr{Ptr{sparseblock}},constraintmatrix,blockmatrix,Ptr{Cdouble},blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},blockmatrix,Ptr{Cdouble},blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cint,paramstruc),n,k,C,a,constant_offset,constraints,byblocks,fill,X,y,Z,cholxinv,cholzinv,pobj,dobj,work1,work2,work3,workvec1,workvec2,workvec3,workvec4,workvec5,workvec6,workvec7,workvec8,diagO,bestx,besty,bestz,Zi,O,rhs,dZ,dX,dy,dy1,Fp,printlevel,parameters)
end

function easy_sdp(n::Cint,k::Cint,C::blockmatrix,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},constant_offset::Cdouble,pX::Ptr{blockmatrix},py::Ptr{Ptr{Cdouble}},pZ::Ptr{blockmatrix})
    pobj = Ref{Cdouble}(.0)
    dobj = Ref{Cdouble}(.0)
    status = ccall((:easy_sdp,CSDP.csdp),Cint,(Cint,Cint,blockmatrix,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,Ptr{blockmatrix},Ptr{Ptr{Cdouble}},Ptr{blockmatrix},Ref{Cdouble},Ref{Cdouble}),n,k,C,a,constraints,constant_offset,pX,py,pZ,pobj,dobj)
    status, pobj[], dobj[]
end

function tweakgap(n::Cint,k::Cint,a::Ptr{Cdouble},constraints::Ptr{constraintmatrix},gap::Cdouble,Z::blockmatrix,dZ::blockmatrix,y::Ptr{Cdouble},dy::Ptr{Cdouble},work1::blockmatrix,work2::blockmatrix,work3::blockmatrix,work4::blockmatrix,workvec1::Ptr{Cdouble},workvec2::Ptr{Cdouble},workvec3::Ptr{Cdouble},workvec4::Ptr{Cdouble},printlevel::Cint)
    ccall((:tweakgap,CSDP.csdp),Void,(Cint,Cint,Ptr{Cdouble},Ptr{constraintmatrix},Cdouble,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},blockmatrix,blockmatrix,blockmatrix,blockmatrix,Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Cint),n,k,a,constraints,gap,Z,dZ,y,dy,work1,work2,work3,work4,workvec1,workvec2,workvec3,workvec4,printlevel)
end

function bisect_(n::Ptr{Cint},eps1::Ptr{Cdouble},d::Ptr{Cdouble},e::Ptr{Cdouble},e2::Ptr{Cdouble},lb::Ptr{Cdouble},ub::Ptr{Cdouble},mm::Ptr{Cint},m::Ptr{Cint},w::Ptr{Cdouble},ind::Ptr{Cint},ierr::Ptr{Cint},rv4::Ptr{Cdouble},rv5::Ptr{Cdouble})
    ccall((:bisect_,CSDP.csdp),Cint,(Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cint},Ptr{Cint},Ptr{Cdouble},Ptr{Cdouble}),n,eps1,d,e,e2,lb,ub,mm,m,w,ind,ierr,rv4,rv5)
end


text:
__precompile__(true)
module CPEG

# greet() = print("Hello World!")

using LinearAlgebra
using StaticArrays
using ForwardDiff
using SparseArrays
using SuiteSparse
using Printf

include("qp_solver.jl")
include("atmosphere.jl")
include("scaling.jl")
include("gravity.jl")
include("aero_forces.jl")
include("vehicle.jl")
include("dynamics.jl")
include("post_process.jl")

end # module


text:
using DependencyTrees: dependents, leftdeps, rightdeps, leftmostdep, rightmostdep
using DependencyTrees: token, tokens

@testset "Graphs" begin

    line = "1	From	_	ADP	IN	_	3	case	_	_"
    d = CoNLLU(line)
    @test strip(DependencyTrees.toconllu(d)) == line

    sent = [
        ("Economic", "ATT", 2),
        ("news", "SBJ", 3),
        ("had", "PRED", 0),
        ("little", "ATT", 5),
        ("effect", "OBJ", 3),
        ("on", "ATT", 5),
        ("financial", "ATT", 8),
        ("markets", "PC", 6),
        (".", "PU", 3),
    ]

    graph = DependencyTree(TypedDependency, sent, add_id=true)

    @test [t.form for t in tokens(graph)] == [first(t) for t in sent]
    for (i, t) in enumerate(sent)
        @test token(graph, i).form == first(t)
        @test [t.form for t in tokens(graph, [1, i])] == ["Economic", first(t)]
    end

    @test DependencyTrees.toconllu(graph) |> strip == """
1	Economic	_	_	_	_	2	ATT	_	_
2	news	_	_	_	_	3	SBJ	_	_
3	had	_	_	_	_	0	PRED	_	_
4	little	_	_	_	_	5	ATT	_	_
5	effect	_	_	_	_	3	OBJ	_	_
6	on	_	_	_	_	5	ATT	_	_
7	financial	_	_	_	_	8	ATT	_	_
8	markets	_	_	_	_	6	PC	_	_
9	.	_	_	_	_	3	PU	_	_
""" |> strip

    @test showstr(graph) |> strip == """
DependencyTree{TypedDependency}
1	Economic	ATT	2
2	news	SBJ	3
3	had	PRED	0
4	little	ATT	5
5	effect	OBJ	3
6	on	ATT	5
7	financial	ATT	8
8	markets	PC	6
9	.	PU	3
""" |> strip

    @test eltype(graph) == DependencyTrees.deptype(graph) == TypedDependency
    @test isroot(root(graph))

    @test length(graph) == length(sent) == 9
    @test isroot(graph[0])
    for (i, (word, tag, id_)) in enumerate(sent)
        @test i == id(graph[i])
        sent_deps = filter(x -> x[3] == i, sent)
        deps_ = dependents(graph, i)
        @test length(deps_) == length(sent_deps)
        @test Set(form(graph, id) for id in deps_) == Set([d[1] for d in sent_deps])
    end
    @test isprojective(graph)

    @test leftdeps(graph, graph[0]) == []
    @test rightdeps(graph, graph[0]) == [graph[3]]

    @test leftdeps(graph, 0) == []
    @test rightdeps(graph, 0) == [3]
    @test leftdeps(graph, 1) == []
    @test rightdeps(graph, 1) == []
    @test leftdeps(graph, 2) == [1]
    @test rightdeps(graph, 2) == []
    @test leftdeps(graph, 3) == [2]
    @test rightdeps(graph, 3) == [5,9]
    @test leftdeps(graph, 4) == []
    @test rightdeps(graph, 4) == []
    @test leftdeps(graph, 5) == [4]
    @test rightdeps(graph, 5) == [6]
    @test leftdeps(graph, 6) == []
    @test rightdeps(graph, 6) == [8]
    @test leftdeps(graph, 7) == []
    @test rightdeps(graph, 7) == []
    @test leftdeps(graph, 8) == [7]
    @test rightdeps(graph, 8) == []
    @test leftdeps(graph, 9) == []
    @test rightdeps(graph, 9) == []

    @test leftmostdep(graph, 0) == -1
    @test leftmostdep(graph, graph[0]).id == -1
    @test leftmostdep(graph.tokens, graph[0]).id == -1
    @test rightmostdep(graph, 0) == 3
    @test rightmostdep(graph, graph[0]).id == 3
    @test rightmostdep(graph.tokens, graph[0]).id == 3

    @test leftmostdep(graph, 1) == -1
    @test leftmostdep(graph, graph[1]).id == -1
    @test leftmostdep(graph.tokens, graph[1]).id == -1
    @test rightmostdep(graph, 1) == -1
    @test rightmostdep(graph, graph[1]).id == -1
    @test rightmostdep(graph.tokens, graph[1]).id == -1

    @test leftmostdep(graph, 2) == 1
    @test leftmostdep(graph, graph[2]).id == 1
    @test leftmostdep(graph.tokens, graph[2]).id == 1
    @test rightmostdep(graph, 2) == -1
    @test rightmostdep(graph, graph[2]).id == -1
    @test rightmostdep(graph.tokens, graph[2]).id == -1

    @test leftmostdep(graph, 3) == 2
    @test leftmostdep(graph, graph[3]).id == 2
    @test leftmostdep(graph.tokens, graph[3]).id == 2
    @test rightmostdep(graph, 3) == 9
    @test rightmostdep(graph, graph[3]).id == 9
    @test rightmostdep(graph.tokens, graph[3]).id == 9

    @test leftmostdep(graph, 4) == -1
    @test leftmostdep(graph, graph[4]).id == -1
    @test leftmostdep(graph.tokens, graph[4]).id == -1
    @test rightmostdep(graph, 4) == -1
    @test rightmostdep(graph, graph[4]).id == -1
    @test rightmostdep(graph.tokens, graph[4]).id == -1

    @test leftmostdep(graph, 5) == 4
    @test leftmostdep(graph, graph[5]).id == 4
    @test leftmostdep(graph.tokens, graph[5]).id == 4
    @test rightmostdep(graph, 5) == 6
    @test rightmostdep(graph, graph[5]).id == 6
    @test rightmostdep(graph.tokens, graph[5]).id == 6

    @test leftmostdep(graph, 6) == -1
    @test leftmostdep(graph, graph[6]).id == -1
    @test leftmostdep(graph.tokens, graph[6]).id == -1
    @test rightmostdep(graph, 6) == 8
    @test rightmostdep(graph, graph[6]).id == 8
    @test rightmostdep(graph.tokens, graph[6]).id == 8

    @test leftmostdep(graph, 7) == -1
    @test leftmostdep(graph, graph[7]).id == -1
    @test leftmostdep(graph.tokens, graph[7]).id == -1
    @test rightmostdep(graph, 7) == -1
    @test rightmostdep(graph, graph[7]).id == -1
    @test rightmostdep(graph.tokens, graph[7]).id == -1

    @test leftmostdep(graph, 8) == 7
    @test leftmostdep(graph, graph[8]).id == 7
    @test leftmostdep(graph.tokens, graph[8]).id == 7
    @test rightmostdep(graph, 8) == -1
    @test rightmostdep(graph, graph[8]).id == -1
    @test rightmostdep(graph.tokens, graph[8]).id == -1

    @test leftmostdep(graph, 9) == -1
    @test leftmostdep(graph, graph[9]).id == -1
    @test leftmostdep(graph.tokens, graph[9]).id == -1
    @test rightmostdep(graph, 9) == -1
    @test rightmostdep(graph, graph[9]).id == -1
    @test rightmostdep(graph.tokens, graph[9]).id == -1

    sent1 = [
        (1, "Economic", "ATT", 2),
        (2, "news", "SBJ", 3),
        (3, "had", "PRED", 0),
        (4, "little", "ATT", 5),
        (5, "effect", "OBJ", 3),
        (6, "on", "ATT", 5),
        (7, "financial", "ATT", 8),
        (8, "markets", "PC", 6),
        (9, ".", "PU", 3),
    ]
    @test graph == DependencyTree(TypedDependency, sent1, add_id=false)

    sent = [
        ("Pierre", "NNP", 2),
        ("Vinken", "NNP", 8),
        (",", ",", 2),
        ("61", "CD", 5),
        ("years", "NNS", 6),
        ("old", "JJ", 2),
        (",", ",", 2),
        ("will", "MD", 0),
        ("join", "VB", 8),
        ("the", "DT", 11),
        ("board", "NN", 9),
        ("as", "IN", 9),
        ("a", "DT", 15),
        ("nonexecutive", "JJ", 15),
        ("director", "NN", 12),
        ("Nov.", "NNP", 9),
        ("29", "CD", 16),
        (".", ".", 8)
    ]

    graph = DependencyTree(TypedDependency, sent, add_id=true)

    @test length(graph) == length(sent) == 18
    @test isroot(graph[0])
    for i in 1:length(sent)
        @test i == id(graph[i])
        sent_deps = filter(x -> x[3] == i, sent)
        deps_ = dependents(graph, i)
        @test length(deps_) == length(sent_deps)
        @test Set(form(graph, id) for id in deps_) == Set([d[1] for d in sent_deps])
    end
    @test isprojective(graph)

    @testset "Projectivity" begin

        # mcdonald & pereira 2005
        # https://www.seas.upenn.edu/~strctlrn/bib/PDF/nonprojectiveHLT-EMNLP2005.pdf
        sent = [
            ("john", 2),      # 1
            ("saw", 0),       # 2
            ("a", 4),         # 3
            ("dog", 2),       # 4
            ("yesterday", 2), # 5
            ("which", 7),     # 6
            ("was", 4),       # 7
            ("a", 9),         # 8
            ("yorkshire", 10),# 9
            ("terrier", 7)    # 10
        ]
        graph = DependencyTree(UntypedDependency, sent, add_id=true)
        @test !isprojective(graph)

        @test DependencyTrees.toconllu(graph) |> strip == """
1	john	_	_	_	_	2	_	_	_
2	saw	_	_	_	_	0	_	_	_
3	a	_	_	_	_	4	_	_	_
4	dog	_	_	_	_	2	_	_	_
5	yesterday	_	_	_	_	2	_	_	_
6	which	_	_	_	_	7	_	_	_
7	was	_	_	_	_	4	_	_	_
8	a	_	_	_	_	9	_	_	_
9	yorkshire	_	_	_	_	10	_	_	_
10	terrier	_	_	_	_	7	_	_	_
""" |> strip

        @test showstr(graph) |> strip == """
DependencyTree{UntypedDependency}
1	john	2
2	saw	0
3	a	4
4	dog	2
5	yesterday	2
6	which	7
7	was	4
8	a	9
9	yorkshire	10
10	terrier	7
""" |> strip

        # jurafsky & martin, speech & language processing (3ed)
        sent = [
            ("jetblue", "nsubj", 2), # 1
            ("canceled", "root", 0), # 2
            ("our", "det", 4),       # 3
            ("flight", "dobj", 2),   # 4
            ("this", "det", 2),      # 5
            ("morning", "nmod", 7),  # 6
            ("which", "case", 4),    # 7
            ("was", "mod", 9),       # 8
            ("already", "adv", 10),  # 9
            ("late", "mod", 7)       # 10
        ]

        graph = DependencyTree(TypedDependency, sent, add_id=true)
        @test !isprojective(graph)
        @test DependencyTrees.deprel(graph, 1) == "nsubj"
    end

    @testset "Errors" begin

        using DependencyTrees: GraphConnectivityError, RootlessGraphError, MultipleRootsError, NonProjectiveGraphError

        noroot = [("no", 2), ("root", 1)]
        @test_throws RootlessGraphError DependencyTree(UntypedDependency, noroot, add_id=true)

        tworoots = [("two", 0), ("roots", 0)]
        @test_throws MultipleRootsError DependencyTree(UntypedDependency, tworoots, add_id=true)
        @test_throws GraphConnectivityError begin
            DependencyTree(UntypedDependency, [(1, "the", 0), (2, "cat", -1)])
        end

        sent = [
            ("john", 2),      # 1
            ("saw", 0),       # 2
            ("a", 4),         # 3
            ("dog", 2),       # 4
            ("yesterday", 2), # 5
            ("which", 7),     # 6
            ("was", 4),       # 7
            ("a", 9),         # 8
            ("yorkshire", 10),# 9
            ("terrier", 7)    # 10
        ]
        @test_throws NonProjectiveGraphError DependencyTree(UntypedDependency, sent, add_id=true, check_projective=true)
    end

end


text:
# Tests for KLMANB

noOfArms = 5

@testset "KLMANB" begin
    agent = KLMANB( noOfArms, 3.0, 4.0 )

    @test agent.noOfArms == noOfArms
    @test agent.σₒ² == 3.0
    @test agent.σₜ² == 4.0
    @test agent.noOfSteps == 0

    @test info_str( agent ) == "KLMANB (σₒ² = 3.000, σₜ² = 4.000)"
    @test info_str( agent, false ) == "KLMANB (σₒ² = 3.000, σₜ² = 4.000)"
    @test info_str( agent, true ) == "KLMANB \$(\\sigma_{o}^2 = 3.000, \\sigma_{t}^2 = 4.000)\$"

    for t = 1:10
        get_arm_index( agent ); update_reward!( agent, randn() );
    end
    @test agent.noOfSteps == 10

    reset!( agent )
    @test agent.noOfArms == noOfArms &&
            agent.σₒ² == 3.0 &&
            agent.σₜ² == 4.0 &&
            agent.noOfSteps == 0
end


text:
using DemoProject
using Documenter

DocMeta.setdocmeta!(DemoProject, :DocTestSetup, :(using DemoProject); recursive=true)

makedocs(;
    modules=[DemoProject],
    authors="Stephan Scholz <stephan.scholz@rwu.de> and contributors",
    repo="https://github.com/stephans3/DemoProject.jl/blob/{commit}{path}#{line}",
    sitename="DemoProject.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://stephans3.github.io/DemoProject.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/stephans3/DemoProject.jl",
    devbranch="main",
)


text:
struct JLSOFile
    version::VersionNumber
    julia::VersionNumber
    format::Symbol
    compression::Symbol
    image::String
    project::Dict{String, Any}
    manifest::Dict{String, Any}
    objects::Dict{Symbol, Vector{UInt8}}
    lock::ReentrantLock
end

"""
    JLSOFile(data; format=:julia_serialize, compression=:gzip, kwargs...)

Stores the information needed to write a .jlso file.

# Arguments

- `data` - The objects to be stored in the file.

# Keywords

- `image=""` - The docker image URI that was used to generate the file
- `julia=$VERSION` - The julia version used to write the file
- `version=v"4"` - The file schema version
- `format=:julia_serialize` - The format to use for serializing individual objects. While `:bson` is
    recommended for longer term object storage, `:julia_serialize` tends to be the faster choice
    for adhoc serialization.
- `compression=:gzip`, what form of compression to apply to the objects.
    Use `:none`, to not compress. `:gzip_fastest` for the fastest gzip compression,
    `:gzip_smallest` for the most compact (but slowest), or `:gzip` for a generally good compromise.
    Due to the time taken for disk IO, `:none` is not normally as fast as using some compression.
"""
function JLSOFile(
    data::Dict{Symbol, <:Any};
    version=v"4",
    julia=VERSION,
    format=:julia_serialize,
    compression=:gzip,
    image=_image(),
)
    if format === :serialize
        # Deprecation warning
        @warn "The `:serialize` format has been renamed to `:julia_serialize`."
        format = :julia_serialize
    end

    _versioncheck(version, WRITEABLE_VERSIONS)
    project_toml, manifest_toml = _env()
    jlso = JLSOFile(
        version,
        julia,
        format,
        compression,
        image,
        Pkg.TOML.parse(project_toml),
        Pkg.TOML.parse(manifest_toml),
        Dict{Symbol, Vector{UInt8}}(),
        ReentrantLock()
    )

    @sync for (key, val) in data
        @spawn jlso[key] = val
    end

    return jlso
end

function JLSOFile(;
    version=v"4",
    julia=VERSION,
    format=:julia_serialize,
    compression=:gzip,
    image=_image(),
    kwargs...
)
    return JLSOFile(
        Dict(kwargs);
        version=version,
        julia=julia,
        format=format,
        compression=compression,
        image=image,
    )
end

JLSOFile(data::Pair{Symbol}...; kwargs...) = JLSOFile(Dict(data); kwargs...)

function Base.show(io::IO, jlso::JLSOFile)
    variables = join(names(jlso), ", ")
    kwargs = join(
        [
            "version=\"$(jlso.version)\"",
            "julia=\"$(jlso.julia)\"",
            "format=:$(jlso.format)",
            "compression=:$(jlso.compression)",
            "image=\"$(jlso.image)\"",
        ],
        ", "
    )

    print(io, "JLSOFile([$variables]; $kwargs)")
end

function Base.:(==)(a::JLSOFile, b::JLSOFile)
    return (
        a.version == b.version &&
        a.julia == b.julia &&
        a.image == b.image &&
        a.manifest == b.manifest &&
        a.format == b.format &&
        a.compression == b.compression &&
        a.objects == b.objects
    )
end

Base.names(jlso::JLSOFile) = collect(keys(jlso.objects))


text:
using LsqFit
using StatsBase

#####
# Multiple ways to estimate an autocorrelation time (tau).
#####
"""
Integrated autocorrelation time. Defined as `tau_int = sum(autocor(M))+0.5`.
"""
function tau_integrated(X::Vector{T}) where T<:Real
    return sum(autocor(X))+0.5
end

"""
Autocorrelation time obtained from fitting autocorrelation C(t) to `exp(-t/tau_fitting)`.
"""
function tau_fitting(X::Vector{T}) where T<:Real
    logautocor = T[]
    Cx = autocor(X)

    firstneg = findfirst(x->x<=0, Cx)
    logCx = firstneg != 0 ? log.(Cx[1:firstneg-1]) : log.(Cx)
    length(logCx) <= 2 ? warn("Fitting log autocorrelation with less than 3 datapoints.") : nothing

    # Fitting
    model(x,p) = p[1].*x+p[2]
    pstart = [-1., 1.]
    fitres = curve_fit(model, 1.0:1.0:length(logCx), logCx, pstart)

    return -1/fitres.param[1]
end

"""
Autocorrelation time obtained from either static or dynamic binning analysis.

Optional keyword `binsize`. Decides wether we do static (fixed, finite binsize)
or dynamic (varying binsize) binning. (Based on Ch. 3.4)
"""
function tau_binning(X::Vector{T}; binsize=0) where T<:Real

    if binsize == 0
        #dynamic binning with varying binsize
        bss, R = R_function(X)
        length(R) > 0 || error("Binning failed. Maybe min_nbins to large?")
        Rplateau = R[end]
        for i in 1:length(bss)-2
            width = abs(maximum(R[i:i+2]) - minimum(R[i:i+2]))
            if width <= 0.05
                Rplateau = R[i]
                break;
            end
        end
        Rplateau != R[end] || warn("Couldn't find plateau. Maybe not converged?")
    else
        #static binning with fixed, single binsize
        Rplateau = R_value(X, binsize)
    end

    return 0.5*(Rplateau-1)
end

#####
# Statistical error (effective standard deviation) calculated from explicit estimate of autocorrelation time
#####
## CORRECT?

"""
Statistical error for correlated data from tau_binning analysis (find plateau (default) or fixed bin size). (Eq. 3.14)
"""
tau_binning_error(X::Vector{T}; binsize=0) where T<:Real = sqrt(var(X)*(1+2*tau_binning(X, binsize=binsize))) # 1/length(X) factor "left out" to have tau_binning_error = std for uncorrelated data
tau_binning_error(X::Vector{T}; binsize=0) where T<:Complex = sqrt(tau_binning_error(real(X), binsize=binsize)^2 + tau_binning_error(imag(X), binsize=binsize)^2)
tau_binning_error(X::Array{T}; binsize=0) where T<:Number = squeeze(mapslices(ts->tau_binning_error(ts, binsize=binsize), X, ndims(X)), ndims(X))

Neff_binning(X::Vector{T}; binsize=0) where T<:Real = floor(Int, length(X)/(1+2*tau_binning(X, binsize=binsize)))

# function tau_binning_error(X::Array{T}; binsize=0) where T<:Real
#     N = size(X, ndims(X)) # length of time series
#     linX = reshape(X, (:, N))
#     Nel = size(linX, 1) # number of elements/time series
#     el_shape = size(X)[1:end-1]

#     errs = zeros(T, Nel)
#     for i in 1:size(errs, 1)
#         errs[i] = tau_binning_error(linX[i,:], binsize=binsize)
#     end
#     return reshape(errs, el_shape)
# end

"""
Statistical error for correlated data from integrated autocorrelation time. (Eq. 3.14)
"""
tau_integrated_error(X::Vector{T}) where T<:Real = sqrt(var(X) *(2*tau_integrated(X)))
tau_integrated_error(X::Vector{T}) where T<:Complex = sqrt(tau_integrated_error(real(X))^2 + tau_integrated_error(imag(X))^2)
tau_integrated_error(X::Array{T}) where T<:Number = squeeze(mapslices(ts->tau_integrated_error(ts), X, ndims(X)), ndims(X))

Neff_integrated(X::Vector{T}) where T<:Real = floor(Int, length(X)/(2*tau_integrated(X)))

"""
Statistical error for correlated data from fitting `C(t) ~ exp(-t/tau)`. (Eq. 3.14)
"""
tau_fitting_error(X::Vector{T}) where T<:Real = sqrt(var(X)*(1+2*tau_fitting(X)))
tau_fitting_error(X::Vector{T}) where T<:Complex = sqrt(tau_fitting_error(real(X))^2 + tau_fitting_error(imag(X))^2)
tau_fitting_error(X::Array{T}) where T<:Number = squeeze(mapslices(ts->tau_fitting_error(ts), X, ndims(X)), ndims(X))

Neff_fitting(X::Vector{T}) where T<:Real = floor(Int, length(X)/(1+2*tau_fitting(X)))


text:
# Tests for validated_integ

using TaylorModels
using LinearAlgebra: norm
using Test

setformat(:full)

# NOTE: IntervalArithmetic v0.16.0 includes this function; but
# IntervalRootFinding is bounded to use v0.15.x
interval_rand(X::Interval{T}) where {T} = X.lo + rand(T) * (X.hi - X.lo)
interval_rand(X::IntervalBox) = interval_rand.(X)

@testset "Tests for `validated_integ`" begin
    @testset "Forward integration" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) =
            (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 0.0, 10.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end

        # initializaton with a Taylor model
        X0tm = qTM[:, 1]
        @assert X0tm isa Vector{TaylorModel1{TaylorN{Float64}, Float64}}
        tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
        @test all(iszero, (qTM - qTM2))

        @taylorize function x_square(dx, x, p, t)
            dx[1] = x[1]^2
            nothing
        end

        exactsol(t, x0) = x0[1] / (1 - x0[1] * t)

        tini, tend = 0., 0.45
        abstol = 1e-15
        orderQ = 5
        orderT = 20
        q0 = [2.]
        δq0 = IntervalBox(-0.1 .. 0.1, Val(1))
        X0 = IntervalBox(q0 .+ δq0)
        ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)

        tTM, qv, qTM = validated_integ(x_square, X0, tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501
        normalized_box = IntervalBox(-1 .. 1, Val(1))

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1, n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:, n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, q0 .+ q0ξ) .∈ q)
            end
        end

    end

    @testset "Forward integration for validated_integ2" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) =
            (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 0.0, 10.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ2(falling_ball!, X0,
            tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end

        @taylorize function x_square(dx, x, p, t)
            dx[1] = x[1]^2
            nothing
        end

        exactsol(t, x0) = x0[1] / (1 - x0[1] * t)

        tini, tend = 0., 0.45
        abstol = 1e-15
        orderQ = 5
        orderT = 20
        q0 = [2.]
        δq0 = IntervalBox(-0.1 .. 0.1, Val(1))
        X0 = IntervalBox(q0 .+ δq0)
        ξ = set_variables("ξₓ", numvars=1, order=2*orderQ)
        tTM, qv, qTM = validated_integ2(x_square, X0, tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501
        normalized_box = IntervalBox(-1 .. 1, Val(1))

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1, n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:, n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, q0 .+ q0ξ) .∈ q)
            end
        end

    end

    @testset "Backward integration" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) = (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 10.0, 0.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ(falling_ball!, X0, tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end

        # initializaton with a Taylor model
        X0tm = qTM[:, 1]
        @assert X0tm isa Vector{TaylorModel1{TaylorN{Float64}, Float64}}
        tTM2, qv2, qTM2 = validated_integ(falling_ball!, X0tm, tini, tend, orderQ, orderT, abstol)
        @test all(iszero, (qTM - qTM2))
    end

    @testset "Backward integration for validated_integ2" begin
        @taylorize function falling_ball!(dx, x, p, t)
            dx[1] = x[2]
            dx[2] = -one(x[1])
            nothing
        end
        exactsol(t, t0, x0) = (x0[1] + x0[2]*(t-t0) - 0.5*(t-t0)^2, x0[2] - (t-t0))

        # Initial conditions
        tini, tend = 10.0, 0.0
        q0 = [10.0, 0.0]
        δq0 = IntervalBox(-0.25 .. 0.25, 2)
        X0 = IntervalBox(q0 .+ δq0)

        # Parameters
        abstol = 1e-20
        orderQ = 2
        orderT = 4
        ξ = set_variables("ξₓ ξᵥ", order=2*orderQ, numvars=length(q0))
        normalized_box = IntervalBox(-1 .. 1, Val(orderQ))

        tTM, qv, qTM = validated_integ2(falling_ball!, X0,
                                        tini, tend, orderQ, orderT, abstol)

        @test length(qv) == length(qTM[1, :]) == length(tTM)
        @test length(tTM) < 501

        for n = 2:length(tTM)
            for it = 1:10
                δt = interval_rand(domain(qTM[1,n]))
                q0ξ = interval_rand(δq0)
                q = evaluate.(evaluate.(qTM[:,n], δt), (normalized_box,))
                @test all(exactsol(tTM[n-1]+δt, tini, q0 .+ q0ξ) .∈ q)
            end
        end
    end

    @testset "Pendulum with constant torque" begin
        @taylorize function pendulum!(dx, x, p, t)
            si = sin(x[1])
            aux = 2 *  si
            dx[1] = x[2]
            dx[2] = aux + 8*x[3]
            dx[3] = zero(x[1])
            nothing
        end
        # Conseerved quantity
        ene_pendulum(x) = x[2]^2/2 + 2 * cos(x[1]) - 8 * x[3]

        # Initial conditions
        tini, tend = 0.0, 12.0
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, 0..0)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        # Parameters
        abstol = 1e-20
        orderQ = 3
        orderT = 10
        ξ = set_variables("ξ", order=2*orderQ, numvars=length(q0))

        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            maxsteps=1800);
        @test all(ene0 .⊆ ene_pendulum.(qv))

        # tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
        #     maxsteps=1800);
        # @test all(ene0 .⊆ ene_pendulum.(qv))

        # Initial conditions 2
        q0 = [1.1, 0.1, 0.0]
        δq0 = IntervalBox(-0.1 .. 0.1, -0.1 .. 0.1, -0.01 .. 0.01)
        X0 = IntervalBox(q0 .+ δq0)
        ene0 = ene_pendulum(X0)

        tTM, qv, qTM = validated_integ(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
            maxsteps=1800);
        @test all(ene0 .⊆ ene_pendulum.(qv))

        # tTM, qv, qTM = validated_integ2(pendulum!, X0, tini, tend, orderQ, orderT, abstol,
        #     maxsteps=1800);
        # @test all(ene0 .⊆ ene_pendulum.(qv))
    end
end


text:
Base.convert(::Type{Expression}, ex::Expr) = Expression(ex)
function Expression(ex;mod=Main)
    ex.head === :if && (ex = Expr(:call, ifelse, ex.args...))
    ex.head === :call || throw(ArgumentError("internal representation does not support non-call Expr"))

    op = getproperty(mod,Symbol(ex.args[1]))
    args = convert.(Expression, ex.args[2:end])

    return Operation(op, args)
end
Base.convert(::Type{Expression}, x::Expression) = x
Base.convert(::Type{Expression}, x::Number) = Constant(x)
Base.convert(::Type{Expression}, x::Bool) = Constant(x)
Base.convert(::Type{Expression}, x::Variable) = convert(Operation,x)
Base.convert(::Type{Expression}, x::Operation) = x
Base.convert(::Type{Expression}, x::Symbol) = Operation(Variable(x),[])
Expression(x::Bool) = Constant(x)

function build_expr(head::Symbol, args)
    ex = Expr(head)
    append!(ex.args, args)
    ex
end

# used in parsing
isblock(x) = length(x) == 1 && x[1] isa Expr && x[1].head == :block
function flatten_expr!(x)
    isb = isblock(x)
    if isb
        x = MacroTools.striplines(x[1])
        filter!(z->z isa Symbol || z.head != :line, x.args)
        x = (x.args...,)
    end
    x
end

function detime_dvs(op::Operation)
  if op.op isa Variable
    Operation(Variable{vartype(op.op)}(op.op.name),Expression[])
  else
    Operation(op.op,detime_dvs.(op.args))
  end
end
detime_dvs(op::Constant) = op

function retime_dvs(op::Operation,dvs,iv)
  if op.op isa Variable && op.op ∈ dvs
    Operation(Variable{vartype(op.op)}(op.op.name),Expression[iv])
  else
    Operation(op.op,retime_dvs.(op.args,(dvs,),iv))
  end
end
retime_dvs(op::Constant,dvs,iv) = op

is_constant(::Constant) = true
is_constant(::Any) = false

is_operation(::Operation) = true
is_operation(::Any) = false

is_derivative(O::Operation) = isa(O.op, Differential)
is_derivative(::Any) = false

Base.occursin(t::Expression) = Base.Fix1(occursin, t)
Base.occursin(t::Expression, x::Operation ) = isequal(x, t) || any(occursin(t), x.args)
Base.occursin(t::Expression, x::Expression) = isequal(x, t)

vars(exprs) = foldl(vars!, exprs; init = Set{Variable}())
function vars!(vars, O)
    isa(O, Operation) || return vars
    O.op isa Variable && push!(vars, O.op)
    for arg ∈ O.args
        if isa(arg, Operation)
            isa(arg.op, Variable) && push!(vars, arg.op)
            vars!(vars, arg)
        end
    end

    return vars
end

# variable extraction
is_singleton(e) = false
is_singleton(e::Operation) = e.op isa Variable

"""
get_variables(O::Operation)

Returns the variables in the Operation
"""
get_variables!(vars, e::Constant, varlist=nothing) = vars
get_variables(e::Constant, varlist=nothing) = get_variables!(Operation[], e, varlist)

function get_variables!(vars, e::Operation, varlist=nothing)
    if is_singleton(e)
      (isnothing(varlist) ? true : (e.op in varlist)) && push!(vars, e)
    else
        foreach(x -> get_variables!(vars, x, varlist), e.args)
    end
    return unique(vars)
end
get_variables(e::Operation, varlist=nothing) = get_variables!(Operation[], e, varlist)

function get_variables!(vars, e::Equation, varlist=nothing)
  get_variables!(vars, e.rhs, varlist)
end
get_variables(e::Equation, varlist=nothing) = get_variables!(Operation[],e,varlist)

modified_states!(mstates, e::Equation, statelist=nothing) = get_variables!(mstates, e.lhs, statelist)

# variable substitution
"""
substitute(expr::Operation, s::Pair)
substitute(expr::Operation, s::Dict)
substitute(expr::Operation, s::Vector)

Performs the substitution `Operation => val` on the `expr` Operation.
"""
substitute(expr::Constant, s) = expr
substitute(expr::Operation, s::Pair) = _substitute(expr, [s[1]], [s[2]])
substitute(expr::Operation, dict::Dict) = _substitute(expr, keys(dict), values(dict))
substitute(expr::Operation, s::Vector) = _substitute(expr, first.(s), last.(s))

function _substitute(expr, ks, vs)
    _substitute(expr, Dict(map(Pair, map(to_symbolic, ks), map(to_symbolic, vs))))
end

function _substitute(expr, dict::Dict)
    simplify(SymbolicUtils.substitute(expr, dict))
end

@deprecate substitute_expr!(expr,s) substitute(expr,s)


text:
# PackageCompiler sysimg compiling.
using Pkg, InstantiateFromURL
github_project("QuantEcon/quantecon-notebooks-julia") # for versioning
pkg"add GR Plots Plots StatsPlots DataFrames StatsPlots CSV PlotUtils GeometryTypes Tables PackageCompiler#sd-notomls CategoricalArrays IteratorInterfaceExtensions PooledArrays WeakRefStrings"
pkg"add Images DualNumbers Unitful Compat LaTeXStrings UnicodePlots DataValues IterativeSolvers VisualRegressionTests"

# new sysimg
using PackageCompiler
syso, sysold = PackageCompiler.compile_incremental(:Plots, install = true)
cp(syso, sysold, force = true)

max_stars_repo_name: Shunsuke-Hori/lecture-source-jl

max_issues_repo_name: Shunsuke-Hori/lecture-source-jl

max_forks_repo_name: Shunsuke-Hori/lecture-source-jl


text: I(parse(Int, readline(stdin)))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=

Definition
==========

https://en.wikipedia.org/wiki/Cauchy_distribution

Function proportional to the "Cauchy distribution".
The formula for cauchy pdf is 1/(pi(x^2/gamma^2 + 1)),

"Proportional to":
- the constant pi is removed from the usual formula.
- negative log is applied
These alterations speed up the computation with no loss of accuracy for the purpose.

negativeLog(x) == log(1/x)

A parameter is the gamma.
AKA the spread/dispersion/deviation.
so-called "autism" by the original author.

"Standard cauchy pdf":
- parameter for the center, Xsub0, is zero
- parameter for the spread, gamma, is 1

As defined, the max value is ~4.3 at max difference of 1.0

Parameters of the PDF
=====================
A parameter is the gamma.
AKA the spread/dispersion/deviation.
so-called "autism" by the original author.

It specifies the spikiness of the PDF peak.
Or, the "fatness" of the tails of the PDF.

Implementation notes
====================

This is a run-time computation when <value> is a variable.

When <value> is a compile-time constant,
does the compiler optimize this away?

When values (the domain of the function) are discrete
this can be evaluated using a table
instead of using run-time arithmetic.

Resynthesizer specific notes
============================

For a discussion of why the standard cauchy PDF is used, instead of the Gaussian,
and other discussion,
see Paul Harrison's thesis https://www.logarithmic.net/pfh-files/thesis/dissertation.pdf.

In the original Resynthesizer typeof value was byte integer (0-254)
i.e. a pixelel or channel value is unsigned byte.
So discrete, and computable using a table.
And the default (usually the best results) gamma parameter was 0.117.

!!! This is at the deepest point of the Resynthesizer computation iteration,
i.e. critical for performance.

Default gamma for resynthesizer
===============================

Should be relative to the domain of the function.
For example, if the max value is 255 (the difference between two byte values)
a gamma of 0.117 gives a certain very spiky curve
(and good results for the Resynthesizer algorithm, by experimental trials.)
When the max value is 1.0 (the float difference between to float values
in the range 0-1)
Then the corresponding default gamma should be 0.000459
to give a similarly shaped PDF curve.

=#

#=
TODO general for type of value
This is for float values in range 0-1

TODO table driven

FUTURE gamma a parameter
=#

#  As in C, log is natural log

#=
Specialized for Real in range 0,1.
The constant derives from the range.
=#
function negLnCauchy(value::Real)

    return log( (value/0.000459)^2 + 1 )
end

max_stars_repo_name: bootchk/resynthesizerJulia

max_issues_repo_name: bootchk/resynthesizerJulia

max_forks_repo_name: bootchk/resynthesizerJulia


text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

"""SphericalPolygon uses StaticArrays to represent vertices. Expect performance degradation for polygons with large numbers of vertices. Performance appears to be good up to perhaps 100 vertices, perhaps as much as 1000 vertices. By 10,000 vertices performance is terrible."""
struct SphericalPolygon{N,T<:Real}
    ptvectors::SMatrix{3,N,T}
    spherecenter::SVector{3,T}
    radius::T

    function SphericalPolygon(points::SMatrix{3,N,T},spherecenter::SVector{3,T},radius::T) where{T<:Real,N}
        normpoints = MMatrix{3,N,T}(undef)

        for i in 1:N
            normpoints[:,i] = normalize(points[:,i] - spherecenter)
        end
        return new{N,T}(SMatrix{3,N,T}(normpoints),SVector{3,T}(spherecenter),radius)
    end
end
export SphericalPolygon

"""returns the spherical angle formed by the cone with centervector at its center with neighbor1,neighbor2 the edges"""
function sphericalangle(neighbor1::SVector{3,T}, centervector::SVector{3,T}, neighbor2::SVector{3,T}) where{T<:Real}
    vec1 = normalize(neighbor1 - (dot(neighbor1,centervector)*centervector)) #subtract off the component of the neighbor vectors from the center vector. This leaves only the component orthogonal to center vector
    vec2 = normalize(neighbor2 - (dot(neighbor2,centervector)*centervector))

    return acos(dot(vec1,vec2))
end

struct SphericalTriangle{T<:Real}
    ptvectors::SMatrix{3,3,T}
    spherecenter::SVector{3,T}
    radius::T

    """Constructor which uses SMatrix to store vectors. For compatibility with LazySets this is more efficient"""
    function SphericalTriangle(vectors::SMatrix{3,3,T},spherecenter::SVector{3,T},radius::T) where{T}
        temp = MMatrix{3,3,T}(undef)
        for col in 1:3
            temp[:,col] = normalize(vectors[:,col]-spherecenter)
        end

        new{T}(temp,spherecenter,radius)
    end
end
export SphericalTriangle

vector(a::SphericalTriangle,i::Int) = a.ptvectors[:,i]

function area(tri::SphericalTriangle{T}) where{T<:Real}
    vec1,vec2,vec3 = vector(tri,1),vector(tri,2),vector(tri,3)
    sum = T(0)
    sum += sphericalangle(vec1,vec3,vec2)
    sum += sphericalangle(vec2,vec1,vec3)
    sum += sphericalangle(vec3,vec2,vec1)
    return (sum-π) * tri.radius^2
end

"""Conceptually breaks the convex spherical polygon into spherical triangles and computes the sum of the angles of all the triangles. The sum of all the angles around the centroid is 2π. Have to subtract π for each of the N triangles. Rather than compute the angles of triangles formed by taking edges from the centroid to each vertex, can instead just compute the internal angle of neighboring edges. Total polygon area is 2π -Nπ + ∑(interior angles)."""
function area(poly::SphericalPolygon{N,T}) where{T<:Real,N}
    return area(poly.ptvectors,poly.radius)
end

function area(ptvecs::SMatrix{3,N,T},radius::T) where{N,T<:Real}
    accum = T(0)

    for i in 2:N-1
        v1 = SVector{3,T}(ptvecs[1,i-1],ptvecs[2,i-1],ptvecs[3,i-1])
        v2 = SVector{3,T}(ptvecs[1,i],ptvecs[2,i],ptvecs[3,i])
        v3 = SVector{3,T}(ptvecs[1,i+1],ptvecs[2,i+1],ptvecs[3,i+1])
        accum += sphericalangle(v1,v2,v3)
    end

    # finish up first and last interior angles which have different indexing because of wraparound
    v1 = SVector{3,T}(ptvecs[1,N-1],ptvecs[2,N-1],ptvecs[3,N-1])
    v2 = SVector{3,T}(ptvecs[1,N],ptvecs[2,N],ptvecs[3,N])
    v3 = SVector{3,T}(ptvecs[1,1],ptvecs[2,1],ptvecs[3,1])

    accum += sphericalangle(v1,v2,v3)

    v1 = SVector{3,T}(ptvecs[1,2],ptvecs[2,2],ptvecs[3,2])
    v2 = SVector{3,T}(ptvecs[1,1],ptvecs[2,1],ptvecs[3,1])
    v3 = SVector{3,T}(ptvecs[1,N],ptvecs[2,N],ptvecs[3,N])

    accum += sphericalangle(v1,v2,v3)

    accum = (2π -N*π + accum)*radius^2
    return accum
end

function circlepoly(nsides; offset = [0.0,0.0,1.0])
    step = 2π/nsides
    result = MMatrix{3,nsides,Float64}(undef)
    for i in 0:nsides-1
        y,x = sincos(step*i)
        result[:,i+1] = [x,y,0.0] .+ offset
    end
    return SMatrix{3,nsides,Float64}(result)
end
export circlepoly

"""creates a circular polygon that subtends a half angle of θ"""
function sphericalcircle(θ, nsides = 10)
    temp = MMatrix{3,nsides,Float64}(undef)
    for i in 0:1:(nsides-1)
        ϕ = i*2π/nsides
        temp[1,i+1] = sin(θ)*cos(ϕ)
        temp[2,i+1] = cos(θ)
        temp[3,i+1] = sin(θ)*sin(ϕ)
    end
    return SphericalPolygon(SMatrix{3,nsides,Float64}(temp),SVector(0.0,0.0,0.0),1.0)
end
export sphericalcircle

oneeigthsphere() = SphericalTriangle(SMatrix{3,3,Float64}(
    0.0,1.0,0.0,
    1.0,0.0,0.0,
    0.0,0.0,1.0),
    SVector(0.0,0.0,0.0),
    1.0)
export oneeigthsphere

onesixteenthphere() = SphericalTriangle(SMatrix{3,3,Float64}(
    0.0,1.0,0.0,
    1.0,1.0,0.0,
    0.0,0.0,1.0),
    SVector(0.0,0.0,0.0),
    1.0)
export onesixteenthphere

threesidedpoly() = SphericalPolygon(SMatrix{3,3,Float64}(
    0.0,1.0,0.0,
    1.0,0.0,0.0,
    0.0,0.0,1.0),
    SVector(0.0,0.0,0.0),
    1.0)
export threesidedpoly

foursidedpoly() = SphericalPolygon(SMatrix{3,4,Float64}(
    0.0,1.0,0.0,
    1.0,1.0,-.9,
    1.0,0.0,0.0,
    0.0,0.0,1.0),
    SVector(0.0,0.0,0.0),
    1.0)
export foursidedpoly


text:
module MangalTestCache
    using Mangal
    using Test
    include("setup.jl")

    # random node
    N = nodes()[1]

    n = node(N.id)

    @test length(Mangal._MANGAL_CACHES[MangalNode]) != 0

end


text:
function dnf_and(exprs::BoolExprVec)::Or
    dnf_exprs = map(dnf, exprs)
    # Each item in dnf_exprs is in DNF form, that is
    # * A variable
    # * A disjunction (Or) of conjunctions (And) of variables
    childDNFs = map(children, dnf_exprs)
    prod = cartprod(childDNFs)
    Or(map(a -> And(children(a)), prod))
end

function dnf_or(exprs::BoolExprVec)::Or
    dnf_exprs = map(dnf, exprs)
    Or(flatmap(children, dnf_exprs))
end

dnf_not(parent::Not, expr::BoolExpr) = parent
dnf_not(parent::Not, expr::And) = dnf(Or(negate(expr.exprs))) # DeMorgan
dnf_not(parent::Not, expr::Or)  = dnf(And(negate(expr.exprs))) # DeMorgan
dnf_not(parent::Not, expr::Not) = dnf(expr.expr) # Double negation
dnf_not(parent::Not, expr::True) = f
dnf_not(parent::Not, expr::False) = t

dnf(expr::Not) = dnf_not(expr, expr.expr)
dnf(expr::And) = dnf_and(expr.exprs)
dnf(expr::Or) = dnf_or(expr.exprs)
dnf(expr::BoolExpr) = expr