text:
export CircularCompactSARTSATrajectory

const CircularCompactSARTSATrajectory = Trajectory{
    SARTSA,
    T1,
    NamedTuple{RTSA,T2},
} where {T1,T2<:Tuple{Vararg{<:CircularArrayBuffer}}}

"""
    CircularCompactSARTSATrajectory(;kwargs...)

Similar to [`VectorialCompactSARTSATrajectory`](@ref),
instead of using `Vector`s as containers, [`CircularArrayBuffer`](@ref)s are used here.

# Key word arguments

- `capacity`::Int, the maximum length of each trace.
- `state_type` = Int
- `state_size` = ()
- `action_type` = Int
- `action_size` = ()
- `reward_type` = Float32
- `reward_size` = ()
- `terminal_type` = Bool
- `terminal_size` = ()
"""
function CircularCompactSARTSATrajectory(;
    capacity,
    state_type = Int,
    state_size = (),
    action_type = Int,
    action_size = (),
    reward_type = Float32,
    reward_size = (),
    terminal_type = Bool,
    terminal_size = (),
)
    capacity > 0 || throw(ArgumentError("capacity must > 0"))
    reward = CircularArrayBuffer{reward_type}(reward_size..., capacity)
    terminal = CircularArrayBuffer{terminal_type}(terminal_size..., capacity)
    state = CircularArrayBuffer{state_type}(state_size..., capacity + 1)
    action = CircularArrayBuffer{action_type}(action_size..., capacity + 1)
    ts = NamedTuple{RTSA}((reward, terminal, state, action))

    CircularCompactSARTSATrajectory{
        Tuple{frame_type(state),frame_type(action),map(frame_type, ts)...},
        typeof(ts).parameters[2],
    }(
        ts,
    )
end

isfull(t::CircularCompactSARTSATrajectory) = isfull(t[:action])

Base.length(t::CircularCompactSARTSATrajectory) = nframes(t[:terminal])


text:
export abs_upper_bound

function _arb_get_fmpq(x::arb)
  mid = ccall((:arb_mid_ptr, :libarb), Ptr{arf_struct}, (Ref{arb}, ), x)
  e = fmpz()
  m = fmpz()
  b = ccall((:arf_get_fmpz_2exp, :libarb), Cint, (Ref{fmpz}, Ref{fmpz}, Ref{arf_struct}), m, e, mid)
  @assert abs(e) < typemax(Int)
  ee = Int(e)
  return fmpq(m, fmpz(1))*fmpq(2)^(ee)
end

function mul2exp!(z::arb, x::arb, y::Int)
  ccall((:arb_mul_2exp_si, :libarb), Nothing, (Ref{arb}, Ref{arb}, Int), z, x, y)
  return nothing
end

function muleq!(z::arb, x::arb, y::arb)
  q = max(bits(x), bits(y))
  ccall((:arb_mul, :libarb), Nothing, (Ref{arb}, Ref{arb}, Ref{arb}, Int), z, x, y, q)
  return nothing
end

function muleq!(z::arb, x::arb, y::fmpz)
  ccall((:arb_mul_fmpz, :libarb), Nothing, (Ref{arb}, Ref{arb}, Ref{fmpz}, Int), z, x, y, bits(x))
  return nothing
end

function addmul!(z::arb, x::arb, y::fmpz)
  q = max(bits(z), bits(x))
  ccall((:arb_addmul_fmpz, :libarb), Nothing, (Ref{arb}, Ref{arb}, Ref{fmpz}, Int), z, x, y, q)
  return nothing
end

function abs!(z::arb, x::arb)
  ccall((:arb_abs, :libarb), Nothing, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
  return nothing
end

function abs!(z::arb, x::acb)
  ccall((:acb_abs, :libarb), Nothing, (Ref{arb}, Ref{acb}, Int), z, x, parent(x).prec)
  return nothing
end

function log!(z::arb, x::arb)
  ccall((:arb_log, :libarb), Nothing, (Ref{arb}, Ref{arb}, Int), z, x, parent(x).prec)
  return nothing
end

################################################################################
#
#  Get upper bounds for absolute value
#
################################################################################

# 1 = ARF_RND_UP = round away from zero
# 2 = ARF_RND_FLOOR = round towards -infinity
# 3 = ARF_RND_CEIL = round towards +infinity

@doc Markdown.doc"""
    abs_upper_bound(x::arb, ::Type{fmpz}) -> fmpz

Returns a positive integer $b$ of type `fmpz` such that $\lvert x \rvert \leq
b$. It is not guarenteed that $b$ is as tight as possible.
"""
function abs_upper_bound(x::arb, ::Type{fmpz})
  tarf = arf_struct(0, 0, 0, 0)

  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), tarf)
  ccall((:arb_get_abs_ubound_arf, :libarb), Nothing,
        (Ref{arf_struct}, Ref{arb}, Int), tarf, x, 64)

  bound = fmpz()

  # round towards +oo
  ccall((:arf_get_fmpz, :libarb), Nothing,
        (Ref{fmpz}, Ref{arf_struct}, Cint), bound, tarf, 3)

  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), tarf)

  return bound
end

@doc Markdown.doc"""
    abs_upper_bound(x::arb, ::Type{Float64}) -> Float64

Returns a positive double $b$ such that $\lvert x \rvert \leq b$. It is not
guarenteed that $b$ is as tight as possible.
"""
function abs_upper_bound(x::arb, ::Type{Float64})
  tarf = arf_struct(0, 0, 0, 0)

  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), tarf)
  ccall((:arb_get_abs_ubound_arf, :libarb), Nothing,
        (Ref{arf_struct}, Ref{arb}, Int),
          tarf, x, 64)

  # round towards -oo
  bound = ccall((:arf_get_d, :libarb), Cdouble,
                (Ref{arf_struct}, Cint), tarf, 3)

  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), tarf)

  return bound
end

@doc Markdown.doc"""
    upper_bound(x::arb, ::Type{fmpz}) -> fmpz

Returns an integer $b$ of type `fmpz` such that $x \leq
b$. It is not guarenteed that $b$ is as tight as possible.
"""
function upper_bound(x::arb, ::Type{fmpz})
  tarf = arf_struct(0, 0, 0, 0)

  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), tarf)
  ccall((:arb_get_ubound_arf, :libarb), Nothing,
        (Ref{arf_struct}, Ref{arb}, Int), tarf, x, 64)

  bound = fmpz()

  # round towards +oo
  ccall((:arf_get_fmpz, :libarb), Nothing,
        (Ref{fmpz}, Ref{arf_struct}, Cint), bound, tarf, 3)

  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), tarf)

  return bound
end

@doc Markdown.doc"""
    lower_bound(x::arb, ::Type{fmpz}) -> fmpz

Returns an integer $b$ of type `fmpz` such that $x \geq
b$. It is not guarenteed that $b$ is as tight as possible.
"""
function lower_bound(x::arb, ::Type{fmpz})
  tarf = arf_struct(0, 0, 0, 0)

  ccall((:arf_init, :libarb), Nothing, (Ref{arf_struct}, ), tarf)
  ccall((:arb_get_lbound_arf, :libarb), Nothing,
        (Ref{arf_struct}, Ref{arb}, Int), tarf, x, 64)

  bound = fmpz()

  # round towards +oo
  ccall((:arf_get_fmpz, :libarb), Nothing,
        (Ref{fmpz}, Ref{arf_struct}, Cint), bound, tarf, 2)

  ccall((:arf_clear, :libarb), Nothing, (Ref{arf_struct}, ), tarf)

  return bound
end

ext: jl


text:
""" custom add_remote_console_cb callback that apply GtkIDE's style"""
function add_remote_console_cb(id, port)
    @info "GtkIDE: Starting console for port $port with id $id"

    c = try
        w = connect(port)

        lang = main_window.style_and_language_manager.languageDefinitions[".jl"]
        c = Console{GtkSourceView, GtkSourceBuffer}(id, main_window, w, (v, b)->init_console!(v, b, main_window), (lang, ))
        GtkREPL.init!(c)

        c.worker_port = port
        GtkREPL.init!(c)

        #for some reason I need to warm-up things here, otherwise it bugs later on.
        @assert remotecall_fetch(identity, GtkREPL.worker(c), 1) == 1

        showall(main_window)
        c
    catch err
        warn(err)
    end

    RemoteGtkREPL.remotecall_fetch(println, worker(c), "Worker connected")
    "done"
end

#hook into GtkREPL `on_command_done`
function on_command_done(main_window::MainWindow, console)
    on_path_change(main_window)
    update!(workspacepanel)
end

#here the index in the notebook isn't updated yet, so it's important to pass `console`
GtkREPL.on_console_mng_switch_page(cm::ConsoleManager, console::Console) = begin
    on_path_change(cm.main_window, false, console)
end


text:
module CellFieldsMocks

using Gridap
using Gridap.CellFieldsOperations: _merge_val_and_grad

using ..CellValuesMocks
using ..FieldsMocks

export IterCellFieldMock
export IndexCellFieldMock
export IterCellBasisMock
export IndexCellBasisMock
export IterCellGeomapMock
export IndexCellGeomapMock

function IterCellFieldMock(D::Int,X::Type,l::Int)
  f = MockField(D,X)
  g = gradient(f)
  cf = TestIterCellValue(f,l)
  cg = TestIterCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

function IndexCellFieldMock(D::Int,X::Type,l::Int)
  f = MockField(D,X)
  g = gradient(f)
  cf = TestIndexCellValue(f,l)
  cg = TestIndexCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

function IterCellBasisMock(D::Int,X::Type,l::Int)
  f = MockBasis(D,X)
  g = gradient(f)
  cf = TestIterCellValue(f,l)
  cg = TestIterCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

function IndexCellBasisMock(D::Int,X::Type,l::Int)
  f = MockBasis(D,X)
  g = gradient(f)
  cf = TestIndexCellValue(f,l)
  cg = TestIndexCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

function IterCellGeomapMock(D::Int,X::Type,l::Int)
  f = MockGeomap(D,X)
  g = gradient(f)
  cf = TestIterCellValue(f,l)
  cg = TestIterCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

function IndexCellGeomapMock(D::Int,X::Type,l::Int)
  f = MockGeomap(D,X)
  g = gradient(f)
  cf = TestIndexCellValue(f,l)
  cg = TestIndexCellValue(g,l)
  _merge_val_and_grad(cf,cg)
end

end # module


text:
@info "Testing ModelNet10..."
@testset "ModelNet10" begin
    _categories = ["sofa", "table"]
    _root = normpath(@__DIR__, "../assets")
    for (split, train) in [("train", true), ("test", false)]
        dset = ModelNet10(root = _root, train = train, categories = _categories)
        @test dset isa Flux3D.Dataset.AbstractDataset
        @test dset.root == normpath(_root)
        @test dset.path == normpath(_root, "ModelNet10")
        @test dset.train == train
        @test dset.transform isa Nothing

        @test dset.length == 2
        @test size(dset) == (2,)
        @test length(dset) == 2

        dpoint = dset[1]
        @test dpoint isa Flux3D.Dataset.DataPoint
        @test dpoint.data isa TriMesh{Float32,UInt32}
        @test dpoint.category_name in _categories

        points = 32
        t = Chain(TriMeshToPointCloud(points))
        dset =
            ModelNet10(root = _root, train = train, categories = _categories, transform = t)
        @test dset.transform isa Chain
        dpoint2 = dset[2]
        @test dpoint2 isa Flux3D.Dataset.DataPoint
        @test dpoint2.data isa PointCloud{Float32}
        @test dpoint2.category_name in _categories
    end
end

@info "Testing ModelNet40..."
@testset "ModelNet40" begin
    _categories = ["desk", "monitor"]
    _root = normpath(@__DIR__, "../assets")
    for (split, train) in [("train", true), ("test", false)]
        dset = ModelNet40(root = _root, train = train, categories = _categories)
        @test dset isa Flux3D.Dataset.AbstractDataset
        @test dset.root == normpath(_root)
        @test dset.path == normpath(_root, "ModelNet40")
        @test dset.train == train
        @test dset.transform isa Nothing

        @test dset.length == 2
        @test size(dset) == (2,)
        @test length(dset) == 2

        dpoint = dset[1]
        @test dpoint isa Flux3D.Dataset.DataPoint
        @test dpoint.data isa TriMesh{Float32,UInt32}
        @test dpoint.category_name in _categories

        points = 32
        t = Chain(TriMeshToPointCloud(points))
        dset =
            ModelNet40(root = _root, train = train, categories = _categories, transform = t)
        @test dset.transform isa Chain
        dpoint2 = dset[2]
        @test dpoint2 isa Flux3D.Dataset.DataPoint
        @test dpoint2.data isa PointCloud{Float32}
        @test dpoint2.category_name in _categories
    end
end


text:
using DiffusionMCMCTools
using Documenter

makedocs(;
    modules=[DiffusionMCMCTools],
    authors="Sebastiano Grazzi, Frank van der Meulen, Marcin Mider, Moritz Schauer",
    repo="https://github.com/JuliaDiffusionBayes/DiffusionMCMCTools.jl/blob/{commit}{path}#L{line}",
    sitename="DiffusionMCMCTools.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaDiffusionBayes.github.io/DiffusionMCMCTools.jl",
        mathengine = Documenter.MathJax(
            Dict(
                :TeX => Dict(
                    :equationNumbers => Dict(
                        :autoNumber => "AMS"
                    ),
                    :Macros => Dict(
                        :dd => "{\\textrm d}",
                        :RR => "\\mathbb{R}",
                        :wt => ["\\widetilde{#1}", 1]
                    ),
                )
            )
        ),
        collapselevel = 1,
    ),
    pages=[
        "Home" => "index.md",
        "Get started" => joinpath("get_started", "overview.md"),
        "User manual" => Any[
            "SamplingUnit" => joinpath("manual", "sampling_unit.md"),
            "SamplingPair" => joinpath("manual", "sampling_pair.md"),
            "SamplingEnsemble" => joinpath("manual", "sampling_ensemble.md"),
            "Block" => joinpath("manual", "block.md"),
            "BiBlock" => joinpath("manual", "biblock.md"),
            "BlockCollection" => joinpath("manual", "block_collection.md"),
            "BlockEnsemble" => joinpath("manual", "block_ensemble.md"),
            "Containers for parameter names" => joinpath("manual", "param_names_collections.md"),
        ],
        "How to..." => Any[
            "..." => joinpath("how_to_guides", "under_construction.md")
        ],
        "Tutorials" => Any[
            "⚠ Preamble ⚠" => joinpath("tutorials", "preamble.md"),
            "BiBlocks" => Any[
                "Smoothing" => joinpath("tutorials", "biblock", "smoothing.md"),
                "Smoothing with blocking" => joinpath("tutorials", "biblock", "smoothing_with_blocking.md"),
                "Inference" => joinpath("tutorials", "biblock", "inference.md"),
                "Inference with blocking" => joinpath("tutorials", "biblock", "inference_with_blocking.md"),
            ],
            "Containers for parameter names" => Any[
                "Inference with BiBlocks" => joinpath("tutorials", "pnames", "inference_with_biblock.md"),
            ],
            "BlockCollection" => Any[
                "Inference" => joinpath("tutorials", "block_collection", "inference.md"),
                "Inference with blocking" => joinpath("tutorials", "block_collection", "inference_with_blocking.md"),
            ],
            "BlockEnsemble" => Any[
                "Inference" => joinpath("tutorials", "block_ensemble", "inference.md"),
            ],
        ],
        "Index" => "module_index.md",
    ],
)

deploydocs(;
    repo="github.com/JuliaDiffusionBayes/DiffusionMCMCTools.jl",
)


text:
struct CalibrationResult{T}
    status::Symbol
    residual::T
    calibration_params::Dict{<:Joint{T}, Vector{T}}
    configurations::Vector{TreeJointSegmentedVector{T}}
    marker_positions::Dict{RigidBody{T}, Vector{Point3DS{T}}}
end

num_poses(result::CalibrationResult) = length(result.configurations)
num_calibration_params(result::CalibrationResult) = sum(length, values(result.calibration_params))
num_markers(result::CalibrationResult) = sum(length, values(result.marker_positions))
RigidBodyDynamics.num_bodies(result::CalibrationResult) = length(result.marker_positions)

function Base.show(io::IO, result::CalibrationResult{T}) where {T}
    println(io, "CalibrationResult{$T}: $(string(result.status)), residual = $(result.residual). Calibration parameters:")
    num_joints = length(result.calibration_params)
    n = 0
    for (joint, params) in result.calibration_params
        print(io, "$(joint.name): $params")
        (n += 1) < num_joints && println(io)
    end
end


text:
# In this file, methods for the descent step calculation are defined.
# The method `get_criticality` takes all the same arguments as `iterate!`
# • `mop` to have access to constraint information
# • `iter_data` for iteration data (x, fx, …)
# • `data_base` (just in case, e.g. if some simplex gradients have to be done on the fly)
# • `sc` to evaluate the surrogate model
# • `algo_config` to have access to `desc_cfg`, an AbstractDescentConfiguration
#
# For the individual descent step methods a method with signature
# `get_criticality( :: typeof{desc_cfg}, :: AbstractMOP, ::AbstractIterData,
#  ::AbstractDB, sc :: SurrogateContainer )`
# should be defined, because this is then called internally
#
# Return: a criticality measure the new trial point (same eltype as x),
# the new trial point model values and the steplength (inf-norm)

# method to only compute criticality
function get_criticality( mop :: AbstractMOP, iter_data :: AbstractIterData,
    data_base :: AbstractDB, sc :: SurrogateContainer, algo_config :: AbstractConfig )
    return get_criticality( descent_method( algo_config ), mop, iter_data, data_base, sc, algo_config )
end

# methods to compute a local descent step
function compute_descent_step( mop :: AbstractMOP, iter_data :: AbstractIterData,
    data_base :: AbstractDB, sc :: SurrogateContainer, algo_config :: AbstractConfig, args... )
    return compute_descent_step( descent_method( algo_config ), mop, iter_data, data_base, sc, algo_config, args... )
end

# fallback, if `get_criticality` allready does the whole job: simply return `ω` and `data`
function compute_descent_step( cfg :: AbstractDescentConfig, mop :: AbstractMOP, iter_data :: AbstractIterData,
    data_base :: AbstractDB, sc :: SurrogateContainer, algo_config :: AbstractConfig )
    ω, data = get_criticality(cfg, mop, iter_data, data_base, sc, algo_config)
    return (ω, data...)
end

# STEEPEST DESCENT AND BACKTRACKING

@with_kw struct SteepestDescentConfig <: AbstractDescentConfig

    "Require a descent in all model objective components."
    strict_backtracking :: Bool = true

    armijo_const_rhs :: Float64 = 1e-3
    armijo_const_shrink :: Float64 = .5

    max_loops :: Int = 30
    min_stepsize :: Float64 = -1.0
end

"Provided `x` and the (surrogate) jacobian `∇F` at `x`, as well as bounds `lb` and `ub`, return steepest multi descent direction."
function _steepest_descent_direction( x :: AbstractVector{F}, ∇F :: Mat, lb :: Vec, ub :: Vec ) where F <: AbstractFloat

    n = length(x);
    try
        opt_problem = JuMP.Model( OSQP.Optimizer );
        JuMP.set_silent(opt_problem);

        JuMP.set_optimizer_attribute( opt_problem, "eps_rel", 1e-5 );
        JuMP.set_optimizer_attribute( opt_problem, "polish", true );

        JuMP.@variable(opt_problem, α)  # criticality measure with flipped sign
        JuMP.@objective(opt_problem, Min, α)

        JuMP.@variable(opt_problem, d[1:n]) # steepest descent direction

        JuMP.@constraint(opt_problem, descent_constraints, ∇F*d .<= α)
        JuMP.@constraint(opt_problem, norm_constraints, -1 .<= d .<= 1 );
        JuMP.@constraint(opt_problem, global_scaled_var_bounds, lb .<= x .+ d .<= ub);

        JuMP.optimize!(opt_problem)
        return F.(JuMP.value.(d)), -F(JuMP.value(α))
    catch e
        println(e)
        @warn("Could not optimize steepest descent subproblem.\n")
        return zeros(F, n), -F(Inf)
    end
end

function _armijo_condition( strict :: Val{true}, Fx, Fx₊, step_size, ω, const_rhs  )
    return all( (Fx .- Fx₊) .>= step_size * const_rhs * ω )
end

function _armijo_condition( strict :: Val{false}, Fx, Fx₊, step_size, ω, const_rhs  )
    return maximum(Fx) - maximum(Fx₊) >= step_size * const_rhs * ω
end

@doc """
Perform a backtracking loop starting at `x` with an initial step of
`step_size .* dir` and return trial point `x₊`, the surrogate value-vector `m_x₊`
and the final step `s = x₊ .- x`.
"""
function _backtrack( x :: AbstractVector{F}, dir, ω, sc, cfg ) where F<:AbstractFloat

    MIN_STEPSIZE = cfg.min_stepsize >= 0 ? cfg.min_stepsize : eps(F)
    MAX_LOOPS = cfg.max_loops
    strict_backtracking = cfg.strict_backtracking
    α = cfg.armijo_const_shrink
    c = cfg.armijo_const_rhs

    step_size = 1

    # values at iterate (required for _armijo_condition )
    mx = eval_models(sc, x)
    # first trial point and its value vector
    x₊ = x .+ step_size .* dir
    mx₊ = eval_models( sc, x₊ )

    for i = 1 : MAX_LOOPS
        if _armijo_condition( Val(strict_backtracking), mx, mx₊, step_size, ω, c )
            break
        end
        if step_size <= MIN_STEPSIZE
            @warn "Could not find a descent by backtracking."
            break
        end
        step_size *= α
        x₊ .= x .+ step_size .* dir
        mx₊ .= eval_models( sc, x₊ )
    end

    step = step_size .* dir;
    return x₊, mx₊, step
end

function get_criticality( desc_cfg :: SteepestDescentConfig , mop, iter_data, data_base, sc, algo_config )

    @logmsg loglevel3 "Calculating steepest descent direction."

    x = get_x( iter_data )
    ∇m = get_jacobian(sc, x)

    lb, ub = full_bounds_internal( mop )

    # compute steepest descent direction and criticality
    d, ω = _steepest_descent_direction( x, ∇m, lb, ub )
    return ω, d
end

function compute_descent_step( desc_cfg :: SteepestDescentConfig , mop, iter_data, data_base, sc, algo_config, ω, d )
    @logmsg loglevel4 "Calculating steepest stepsize."

    x = get_x( iter_data )
    Δ = get_delta( iter_data )

    # scale direction for backtracking as in paper
    norm_d = norm(d, Inf)
    if norm_d > desc_cfg.min_stepsize
        d_normed = d ./ norm_d
        σ = intersect_bounds( mop, x, Δ, d_normed; return_vals = :pos )
        # Note: For scalar Δ the above should equal
        # `σ = norm(d,Inf) < 1 || Δ <= 1 ? min( d, norm(d,Inf) ) : Δ`

        if σ > desc_cfg.min_stepsize
            x₊, mx₊, step = _backtrack( x, σ * d_normed, ω, sc, desc_cfg )
            return ω, x₊, mx₊, norm(step, Inf)
        end
    end

    return 0, copy(x), eval_models( sc, x ), 0
end

function compute_descent_step( desc_cfg :: SteepestDescentConfig , mop, iter_data, data_base, sc, algo_config )
    ω, d = get_criticality( desc_cfg, mop, iter_data, data_base, sc, algo_config )
    return compute_descent_step(desc_cfg, mop, iter_data, data_base, sc, algo_config, ω, d )
end

# PASCOLETTI-SERAFINI
# TODO would it matter to allow single precision reference_point?
# Does NLopt honour precision?
@with_kw struct PascolettiSerafiniConfig <: AbstractDescentConfig

    "Local utopia point to guide descent search in objective space."
    reference_point :: Vector{Float64} = Float64[]

    "Objective space direction to guide local subproblem solution."
    reference_direction :: Vector{Float64} = Float64[]

    trust_region_factor :: Float64 = 1.0

    max_ps_problem_evals :: Int = -1
    max_ps_polish_evals :: Int = -1

    max_ideal_point_problem_evals :: Int = -1

    "Algorithm used for the Pascoletti Serafini subproblem minimization."
    main_algo :: Symbol = :GN_ISRES

    "Algorithm to compute local reference point if no point or direction is given."
    reference_algo :: Symbol = :GN_ISRES
    reference_trust_region_factor :: Float64 = 1.1

    "Specify local algorithm to polish Pascoletti-Serafini solution. Uses 1/4 of maximum allowed evals."
    ps_polish_algo :: Union{Nothing, Symbol} = nothing
end

function _get_local_dir( cfg :: PascolettiSerafiniConfig, fx, sorting_indices )
    if !isempty( cfg.reference_direction )
        return cfg.reference_direction[ sorting_indices ]
    elseif !isempty( cfg.reference_point )
        return fx .- cfg.reference_point[ sorting_indices ]
    else
        return nothing
    end
end

function _min_component( x, n_vars, lb, ub, algo_name, MAX_EVALS, opt_handle )
    opt = NLopt.Opt( algo_name, n_vars )
    opt.lower_bounds = lb
    opt.upper_bounds = ub
    opt.xtol_rel = 1e-3
    opt.maxeval = MAX_EVALS
    #opt.min_objective = get_optim_handle( sc, mop, l )
    opt.min_objective = opt_handle
    minf, _ = NLopt.optimize( opt, x );
    return minf
end

function _local_ideal_point( x :: AbstractVector{F}, lb, ub, optim_handles, max_evals, algo_name ) where F
    @logmsg loglevel4 "Computing local ideal point. This can take a bit…"

    num_objectives = length( optim_handles )
    n_vars = length(x)

    MAX_EVALS = max_evals < 0 ? 500 * (n_vars+1) : max_evals

    # preallocate local ideal point:
    ȳ = fill( -F(Inf), num_objectives )

    # minimize each individual scalar surrogate output
    for (l, opt_handle) in enumerate( optim_handles )
        ȳ[l] = _min_component(x, n_vars, lb, ub, algo_name, MAX_EVALS, opt_handle)
    end
    return ȳ
end

function get_criticality( desc_cfg :: PascolettiSerafiniConfig, mop, iter_data, data_base, sc, algo_config )

    @logmsg loglevel3 "Calculating Pascoletti-Serafini descent."
    x = get_x( iter_data )
    fx = get_fx( iter_data )
    Δ = get_delta( iter_data )
    n_vars = length( x )
    n_out = num_objectives(mop);

    r = _get_local_dir( desc_cfg, fx, mop )
    if isnothing(r)
        opt_handles = [ get_optim_handle(sc, l) for l = 1 : num_objectives ]
        lb, ub = local_bounds(mop, x, desc_cfg.reference_trust_region_factor * Δ )
        ip = _local_ideal_point(x,lb, ub, opt_handles, desc_cfg.MAX_EVALS, desc_cfg.algo_name)
        r = fx .- ip
    end

    @logmsg loglevel4 "Local image direction is $(_prettify(r))"

    # If any component is not positive, we are critical
    if any( r .<= 0 )
        return 0, copy(x), eval_models(sc, x), 0
    end

    # total number of sub solver evaluations
    MAX_EVALS = desc_cfg.max_ps_problem_evals < 0 ? 500 * (n_vars+1) : desc_cfg.max_ps_problem_evals
    # number of evaluations for the global optimization algo
    if isnothing( desc_cfg.ps_polish_algo )
        MAX_EVALS_global = MAX_EVALS
        MAX_EVALS_local = 0
    else
        if desc_cfg.max_ps_polish_evals < 0
            MAX_EVALS_global = Int( floor( MAX_EVALS*3/4 ) )
            MAX_EVALS_local = MAX_EVALS - MAX_EVALS_global
        else
            MAX_EVALS_global = MAX_EVALS
            MAX_EVALS_local = desc_cfg.max_ps_polish_evals
        end
    end

    mx = eval_models( sc, x );
    lb, ub = local_bounds(mop, x, decs_cfg.trust_region_factor * Δ)

    τ, χ_min, ret = _ps_optimization(sc, desc_cfg.main_algo, lb, ub,
        MAX_EVALS_global, [-0.5;x], mx, r, n_vars, n_out)

    if MAX_EVALS_local > 0
        @logmsg loglevel4 "Local polishing for PS enabled."

        τₗ, χ_minₗ, retₗ = _ps_optimization(sc, desc_cfg.ps_polish_algo, lb, ub,
            MAX_EVALS_local, χ_min, mx, r, n_vars, n_out);
        if !(retₗ == :FAILURE || isinf(τₗ) || isnan(τₗ) || any(isnan.(χ_minₗ)) )
            τ, χ_min, ret = τₗ, χ_minₗ, retₗ
        end
    end

    if (ret == :FAILURE || isinf(τ) || isnan(τ) || any(isnan.(χ_min)) )
        return 0, copy(x), eval_models(sc, x), 0
    else
        ω = abs( τ );
        x₊ = χ_min[2 : end];
        mx₊ = eval_models( sc, x₊ );
        sl = norm( x .- x₊, Inf );

        return ω, x₊, mx₊, sl
    end
end

"Construct and solve Pascoletti Serafini subproblem using surrogates from `sc`."
function _ps_optimization( sc :: SurrogateContainer, algo :: Symbol,
    lb :: Vec, ub :: Vec, MAX_EVALS :: Int, χ0 :: AbstractVector{F},
    mx :: Vec, r :: Vec, n_vars :: Int, n_out :: Int ) where F

    opt = NLopt.Opt( algo, n_vars + 1 );
    opt.lower_bounds = [-1.0 ; lb ];
    opt.upper_bounds = [ 0.0 ; ub ];
    opt.xtol_rel = 1e-3;
    opt.maxeval = MAX_EVALS;
    opt.min_objective = _get_ps_objective_func();

    #NLopt.inequality_constraint!( opt, _get_ps_constraint_func(sc,mx,r), 1e-12 );
    for l = 1 : n_out
        NLopt.inequality_constraint!( opt, _get_ps_constraint_func(sc, mx, r, l), eps(F) )
    end

    @logmsg loglevel4 "Starting PS optimization."
    τ, χ_min, ret = NLopt.optimize(opt, χ0 );
    @logmsg loglevel4 "Finished with $(ret) after $(opt.numevals) model evals."

    return τ, χ_min, ret
end

#=
NOTE Does not work using the whole jacobian somehow
function _get_ps_constraint_func( sc :: SurrogateContainer, mx :: RVec, dir :: RVec ) :: Function
    # return the l-th constraint functions for pascoletti_serafini
    # dir .>= 0 is the image direction
    # χ = [t;x] is the augmented variable vector
    # g is a matrix of gradients (columns are grads, hence g is a jacobian transposed)
    return function(χ, g)
        if !isempty(g)
            g[1, :] .= -dir;
            g[2:end, :] .= transpose(get_jacobian( sc, χ[2:end]));
        end
        return eval_models(sc, χ[2:end]) .- mx .- χ[1] .* dir
    end
end
=#

"""
    _get_ps_constraint_func( sc :: SurrogateContainer, mx, dir, l )

Return the `l`-th (possibly non-linear) constraint function
for Pascoletti-Serafini.
`dir` .>= 0 is the image direction;
`χ = [t;x]` is the augmented variable vector;

"""
function _get_ps_constraint_func( sc :: SurrogateContainer, mx, dir, l )
    return function(χ, g)
        if !isempty(g)
            g[1] = -dir[l]
            g[2:end] .= get_gradient(sc, χ[2:end], l)
        end
        ret_val = eval_models(sc, χ[2:end], l) - mx[l] - χ[1] * dir[l]
        return ret_val
    end
end

"""
Return objective function for Pascoletti-Serafini, modifying gradient in place.
"""
function _get_ps_objective_func()
    function( χ, g )
        if !isempty(g)
            g[1] = 1.0;
            g[2:end] .= 0.0;
        end
        return χ[1]
    end
end

# Directed Search
#=
function compute_descent_step(::Val{:directed_search}, algo_config :: AbstractConfig,
    mop :: AbstractMOP, id :: AbstractIterData, sc :: SurrogateContainer )
    pascoletti_serafini(Val(ds), algo_config, mop, id, sc )
end

function compute_descent_step(::Val{:ds}, algo_config :: AbstractConfig,
    mop :: AbstractMOP, id :: AbstractIterData, sc :: SurrogateContainer )
    @logmsg loglevel3 "Calculating Pascoletti-Serafini descent."
    x = xᵗ(id);
    fx = fxᵗ(id);
    n_vars = num_vars(mop);
    n_out = num_objectives(mop);

    # very similar to PS, but inverted sign(s)
    r = begin
        if !isempty( reference_direction(algo_config) )
            dir = apply_internal_sorting( reference_direction(algo_config), mop )
            if all( dir .>= 0 )
                dir *= -1
            end
        elseif !isempty( reference_point(algo_config))
            dir = apply_internal_sorting(reference_point(algo_config), mop) .- fx
        else
            dir = _local_ideal_point(mop, algo_config, id, sc) .- fx
        end
        dir
    end

    @logmsg loglevel4 "Local image direction is $(_prettify(r))"

    if any( r .>= 0 )
        return 0, copy(x), eval_models(sc, x), 0
    end

    lb, ub = full_bounds_internal( mop );
    mx = eval_models( sc, x );
    ∇m = get_jacobian(sc, x);

    ∇m⁺ = pinv( ∇m )
    CONSTRAINED = !isempty(MOI.get(mop,MOI.ListOfConstraints()))
    if !CONSTRAINED
        dir = ∇m⁺*r;
    else
        dir_prob = JuMP.Model(OSQP.Optimizer);
        JuMP.set_silent(dir_prob)

        # TODO are parameters sensible?
        JuMP.set_optimizer_attribute(dir_prob,"eps_rel",1e-5)
        JuMP.set_optimizer_attribute(dir_prob, "polish", true)

        # slightly adapted version of directed search (allow only non-anscent directions)
        JuMP.@variable(dir_prob, d[1:n_vars] )
        JuMP.@objective(dir_prob, Min, sum(( ∇m * d .- r ).^2) )
        JuMP.@constraint(dir_prob, norm_constraint, -1.0 .<= d .<= 1.0)
        JuMP.@constraint(dir_prob, descent, ∇m*d .<= 0)
        JuMP.@constraint(dir_prob, global_constraint, lb .<= x .+ d .<= ub )

        JuMP.optimize!(dir_prob)
        dir = JuMP.value.(d)
    end
    ω = - maximum(∇m * dir);

    # scale direction for backtracking as in paper
    norm_d = norm(dir,Inf);
    if norm_d > 0
        d_normed = dir ./ norm_d;
        σ = intersect_bounds( mop, x, Δᵗ(id), d_normed; return_vals = :pos )
        # Note: For scalar Δ the above should equal
        # `σ = norm(d,Inf) < 1 || Δ <= 1 ? min( d, norm(d,Inf) ) : Δ`
        @assert σ >= 0

        x₊, mx₊, step = _backtrack( x, d_normed, σ, ω, sc, strict_backtracking(algo_config) );
        return ω, x₊, mx₊, norm(step, Inf)
    else
        return ω, copy(x), eval_models( sc, x ), 0
    end
end

TODO Re-enable Directed Search
=#

####### Symbol values as quick configuration

function _cfg_from_symbol( desc_cfg :: Symbol, F :: Type )
    if desc_cfg == :steepest || desc_cfg == :sd || desc_cfg == :steepest_descent
        return SteepestDescentConfig()
    elseif desc_cfg == :ps || desc_cfg == :pascoletti_serafini
        return PascolettiSerafiniConfig()
    end
end

# DefaultDescent -> Fallback to SteepestDescent
function compute_descent_step( desc_cfg :: Symbol, mop, iter_data, data_base, sc, algo_config, args... )
    F = eltype(get_x(iter_data))
    true_cfg = _cfg_from_symbol( desc_cfg, F )
    return compute_descent_step( true_cfg, mop, iter_data, data_base, sc, algo_config, args... )
end

# DefaultDescent -> Fallback to SteepestDescent
function get_criticality( desc_cfg :: Symbol, mop, iter_data, data_base, sc, algo_config, args... )
    F = eltype(get_x(iter_data))
    true_cfg = _cfg_from_symbol( desc_cfg, F )
    return get_criticality( true_cfg, mop, iter_data, data_base, sc, algo_config, args... )
end


text:
# defines API for discrete MDPs and their solvers
module DiscreteMDPs

# basic types
export DiscreteMDP, Solver, Simulator, Policy, MDPFile, PolicyFile

abstract DiscreteMDP
abstract Solver
abstract Simulator
abstract Policy
abstract MDPFile
abstract PolicyFile

# module function
export solve, simulate

solve(solver::Solver, mdp::DiscreteMDP) = error("$(typeof(solver)) does not implement solve for MDP $(typeof(mdp))")
simulate(simulator::Simulator, mdp::DiscreteMDP, policy::Policy) = error("$(typeof(simulator)) does not implement simulate for MDP $(typeof(mdp)) and policy $(typeof(policy))")

# file io
export load, save

load(file::PolicyFile) = error("$(typeof(file)) does not implement load")
load(file::MDPFile) = error("$(typeof(file)) does not implement load")
save(file::PolicyFile, policy::Policy) = error("$(typeof(file)) does not implement save for policy $(typeof(policy))")
save(file::MDPFile, mdp::DiscreteMDP) = error("$(typeof(file)) does not implement save for DiscreteMDP $(typeof(mdp))")

# model function
export states, actions, reward, transition, nextStates, numStates, numActions

states(mdp::DiscreteMDP) = error("$(typeof(mdp)) does not implement states") # returns something that might be iterable (e.g., 1:100)
actions(mdp::DiscreteMDP) = error("$(typeof(mdp)) does not implement actions") # returns something that might be iterable (e.g., 1:100)
reward(mdp::DiscreteMDP, state, action) = error("$(typeof(mdp)) does not implement reward") # returns reward
transition(mdp::DiscreteMDP, state, action, nextState) = error("$(typeof(mdp)) does not implement transition") # returns probability
nextStates(mdp::DiscreteMDP, state, aciton) = error("$(typeof(mdp)) does not implement transition") # returns neighbor states and prob vectors

numStates(mdp::DiscreteMDP) = length(collect(states(mdp)))
numActions(mdp::DiscreteMDP) = length(collect(actions(mdp)))

# policy functions
export action, value

action(p::Policy, state) = error("$(typeof(p)) does not implement action") # returns the action according to the policy p
value(p::Policy, state, action) = error("$(typeof(p)) does not implement value") # returns the expected value for a (s,a) pair
value(p::Policy, state) = error("$(typeof(p)) does not implement value") # returns the optimal expeted value for state s

end # module


text:
using VortexDistributions, Test

# @testset "Multivortex creation and detection" begin
@test found_near(30)
# end


text:
## Rejuvenation moves for SMC algorithms ##
export pf_rejuvenate!, pf_move_accept!, pf_move_reweight!
export move_reweight

"""
    pf_rejuvenate!(state::ParticleFilterState, kern, kern_args::Tuple=(),
                   n_iters::Int=1; method=:move)

Rejuvenates particles by repeated application of a kernel `kern`. `kern`
should be a callable which takes a trace as its first argument, and returns
a tuple with a trace as the first return value. `method` specifies the
rejuvenation method: `:move` for MCMC moves without a reweighting step,
and `:reweight` for rejuvenation with a reweighting step.
"""
function pf_rejuvenate!(state::ParticleFilterView, kern, kern_args::Tuple=(),
                        n_iters::Int=1; method::Symbol=:move)
    if method == :move
        return pf_move_accept!(state, kern, kern_args, n_iters)
    elseif method == :reweight
        return pf_move_reweight!(state, kern, kern_args, n_iters)
    else
        error("Method $method not recognized.")
    end
end

"""
    pf_move_accept!(state::ParticleFilterState, kern,
                    kern_args::Tuple=(), n_iters::Int=1)

Rejuvenates particles by repeated application of a MCMC kernel `kern`. `kern`
should be a callable which takes a trace as its first argument, and returns
a tuple `(trace, accept)`, where `trace` is the (potentially) new trace, and
`accept` is true if the MCMC move was accepted. Subsequent arguments to `kern`
can be supplied with `kern_args`. The kernel is repeatedly applied to each trace
for `n_iters`.
"""
function pf_move_accept!(state::ParticleFilterView,
                         kern, kern_args::Tuple=(), n_iters::Int=1)
    # Potentially rejuvenate each trace
    for (i, trace) in enumerate(state.traces)
        for k = 1:n_iters
            trace, accept = kern(trace, kern_args...)
            @debug "Accepted: $accept"
        end
        state.new_traces[i] = trace
    end
    update_refs!(state)
    return state
end

"""
    pf_move_reweight!(state::ParticleFilterState, kern,
                      kern_args::Tuple=(), n_iters::Int=1)

Rejuvenates and reweights particles by repeated application of a reweighting
kernel `kern`, as described in [1]. `kern` should be a callable which takes a
trace as its first argument, and returns a tuple `(trace, rel_weight)`,
where `trace` is the new trace, and `rel_weight` is the relative log-importance
weight. Subsequent arguments to `kern` can be supplied with `kern_args`.
The kernel is repeatedly applied to each trace for `n_iters`, and the weights
accumulated accordingly. Both the [`move_reweight`](@ref) function and
[symmetric trace translators](https://www.gen.dev/stable/ref/trace_translators/)
can serve as reweighting kernels.

[1] R. A. G. Marques and G. Storvik, "Particle move-reweighting strategies for
online inference," Preprint series. Statistical Research Report, 2013.
"""
function pf_move_reweight!(state::ParticleFilterView,
                           kern, kern_args::Tuple=(), n_iters::Int=1)
    # Move and reweight each trace
    for (i, trace) in enumerate(state.traces)
        weight = 0
        for k = 1:n_iters
            trace, rel_weight = kern(trace, kern_args...)
            weight += rel_weight
            @debug "Rel. Weight: $rel_weight"
        end
        state.new_traces[i] = trace
        state.log_weights[i] += weight
    end
    update_refs!(state)
    return state
end

"""
    move_reweight(trace, selection)
    move_reweight(trace, proposal, proposal_args)
    move_reweight(trace, proposal, proposal_args, involution)
    move_reweight(trace, proposal_fwd, args_fwd,
                  proposal_bwd, args_bwd, involution)

Move-reweight MCMC kernel, which takes in a `trace` and returns a new trace
along with a relative importance weight. This can be used for rejuvenation
within a particle filter, as described in [1].

Several variants of `move_reweight` exist, differing in the complexity
involved in proposing and re-weighting random choices:

1. The `selection` variant regenerates the addresses specified by `selection`
   using the model's default proposal.
2. A `proposal` generative function and arguments can be provided, to propose
   new random choices using a custom proposal distribution. `proposal` must
   take the original trace as its first input argument.
3. An `involution` can also be provided, to handle more complex proposal
   distributions which add/remove trace addresses.
4. Separate forward and backward proposal distributions, `proposal_fwd` and
   `proposal_bwd`, can be provided, as long they share the same support. This
   adjusts the computation of the relative importance weight by scoring
   the backward choices under the backward proposal.

[1] R. A. G. Marques and G. Storvik, "Particle move-reweighting strategies for
online inference," Preprint series. Statistical Research Report, 2013.
"""
function move_reweight(trace::Trace, selection::Selection)
    args = get_args(trace)
    argdiffs = map((_) -> NoChange(), args)
    new_trace, rel_weight = regenerate(trace, args, argdiffs, selection)
    return new_trace, rel_weight
end

function move_reweight(trace::Trace, proposal::GenerativeFunction,
                       proposal_args::Tuple)
    model_args = Gen.get_args(trace)
    argdiffs = map((_) -> NoChange(), model_args)
    fwd_choices, fwd_score, fwd_ret =
        propose(proposal, (trace, proposal_args...,))
    new_trace, weight, _, discard =
        update(trace, model_args, argdiffs, fwd_choices)
    bwd_score, bwd_ret =
        assess(proposal, (new_trace, proposal_args...), discard)
    rel_weight = weight - fwd_score + bwd_score
    return new_trace, rel_weight
end

function move_reweight(trace::Trace, proposal::GenerativeFunction,
                       proposal_args::Tuple, involution)
    fwd_choices, fwd_score, fwd_ret =
        propose(proposal, (trace, proposal_args...,))
    new_trace, bwd_choices, weight =
        involution(trace, fwd_choices, fwd_ret, proposal_args)
    bwd_score, bwd_ret =
        assess(proposal, (new_trace, proposal_args...), bwd_choices)
    rel_weight = weight - fwd_score + bwd_score
    return new_trace, rel_weight
end

function move_reweight(trace::Trace, proposal_fwd::GenerativeFunction,
                       args_fwd::Tuple, proposal_bwd::GenerativeFunction,
                       args_bwd::Tuple, involution)
    fwd_choices, fwd_score, fwd_ret =
        propose(proposal_fwd, (trace, args_fwd...,))
    new_trace, bwd_choices, weight =
        involution(trace, fwd_choices, fwd_ret, args_fwd)
    bwd_score, bwd_ret =
        assess(proposal_bwd, (new_trace, args_bwd...), bwd_choices)
    rel_weight = weight - fwd_score + bwd_score
    return new_trace, rel_weight
end


text:
using Test

@testset "broadcasting" begin

    @testset "number_to_bdint" begin
        a = BD{Array{Int}}(zeros(Int,3,3))

        y, dy = (1 .+ a).f
        @test dy(1) == (1.0, [1 1 1; 1 1 1; 1 1 1])

        y, dy = (1.0 .+ a).f
        @test dy(1.0) == (1.0, [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0])
    end

    @testset "number_to_bdfloat" begin
        a = BD{Array{Float64}}(zeros(Float64,3,3))

        y, dy = (1 .+ a).f
        @test dy(1) == (1.0, [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0])

        y, dy = (1.0 .+ a).f
        @test dy(1.0) == (1.0, [1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0])
    end

    @testset "bdint_to_number" begin
        a = BD{Array{Int}}(zeros(Int,3,3))

        y, dy = (a .+ 1).f
        @test dy(1) == ([1 1 1; 1 1 1; 1 1 1], 1)

        y, dy = (a .+ 1.0).f
        @test dy(1.0) == ([1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0], 1.0)
    end

    @testset "bdfloat_to_number" begin
        a = BD{Array{Float64}}(zeros(Float64,3,3))

        y, dy = (a .+ 1).f
        @test dy(1) == ([1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0], 1.0)

        y, dy = (a .+ 1.0).f
        @test dy(1.0) == ([1.0 1.0 1.0; 1.0 1.0 1.0; 1.0 1.0 1.0], 1.0)
    end
end


text:
export ELFHeader, ELFHeader32, ELFHeader64

import Base: show

"""
    ELFHeader

ELF Header type, containing metadata about the ELF object itself, such as its
type, offsets for the Program and Section headers, the number of other header
entries, etc...
"""
abstract type ELFHeader{H <: ObjectHandle} end
@io struct ELFHeader32{H <: ObjectHandle} <: ELFHeader{H}
    e_type::UInt16
    e_machine::UInt16
    e_version::UInt32
    e_entry::UInt32
    e_phoff::UInt32
    e_shoff::UInt32
    e_flags::UInt32
    e_ehsize::UInt16
    e_phentsize::UInt16
    e_phnum::UInt16
    e_shentsize::UInt16
    e_shnum::UInt16
    e_shstrndx::UInt16
end
@io struct ELFHeader64{H <: ObjectHandle} <: ELFHeader{H}
    e_type::UInt16
    e_machine::UInt16
    e_version::UInt32
    e_entry::UInt64
    e_phoff::UInt64
    e_shoff::UInt64
    e_flags::UInt32
    e_ehsize::UInt16
    e_phentsize::UInt16
    e_phnum::UInt16
    e_shentsize::UInt16
    e_shnum::UInt16
    e_shstrndx::UInt16
end

function filetype(e_type)
    if haskey(ET_TYPES, e_type)
        return ET_TYPES[e_type]
    end
    return string("Unknown (0x",string(e_type, base=16),")")
end

function machinetype(e_machine)
    if haskey(EM_MACHINES, e_machine)
        return EM_MACHINES[e_machine]
    end
    return string("Unknown (0x",string(e_machine, base=16),")")
end

function show(io::IO, header::ELFHeader)
    println(io, "ELF Header")
    println(io, "  Type ", filetype(header.e_type))
    println(io, "  Machine ", machinetype(header.e_machine))
    # Skip e_version (not particularly useful)
    println(io, "  Entrypoint ", "0x", string(header.e_entry, base=16))
    println(io, "  PH Offset ", "0x", string(header.e_phoff, base=16))
    println(io, "  SH Offset ", "0x", string(header.e_shoff, base=16))
    # Skip flags
    println(io, "  Header Size ", "0x", string(header.e_ehsize, base=16))
    println(io, "  PH Entry Size ", "0x", string(header.e_phentsize, base=16))
    println(io, "  PH Entry Count ", string(header.e_phnum, base=10))
    println(io, "  SH Entry Size ", "0x", string(header.e_shentsize, base=16))
    println(io, "  SH Entry Count ", string(header.e_shnum, base=10))
    println(io, "  Strtab Index ", string(header.e_shstrndx, base=10))
end


text:
#=
The following codes are based on old pmm related functions written by sammy-suyama.
https://github.com/sammy-suyama/BayesBook/blob/master/src/PoissonMixtureModel.jl
=#

const PMM = MixtureModel{Univariate, Discrete, Poisson}

struct BayesianPMM
    α
    prior::Vector{Gamma}
    BayesianPMM(α, prior) = new(α, prior)
end

function rand(model::BayesianPMM)
    K = length(model.prior)
    a = rand(Dirichlet(K, model.α))
    λ = rand.(model.prior)
    MixtureModel(Poisson, λ, a)
end

function one_hot(model::PMM, N::Int)
    K = length(model.components)
    a = probs(model) # mixture weight

    # one hot vector of (K, N)
    onehotbatch(rand(Categorical(a), N), collect(1:K))
end

function generate_toy_data(model::PMM, N::Int)
    D = length(model)

    data = Array{eltype(model.components[1])}(undef, N, D)
    # convert to vector if D == 1
    D == 1 ? data = dropdims(data, dims = 2) : data
    S = one_hot(model, N)
    @inbounds for n = 1:N
        s = argmax(S[:, n])
        data[n, :] .= rand(model.components[s])
    end
    data, S
end

function init_S(data, prior_model::BayesianPMM)
    N = size(data, 1)
    K = length(prior_model.prior)
    onehotbatch(rand(Categorical(ones(K)./K), N), collect(1:K))
end

function add_stats(prior_model::BayesianPMM, data, S)
    @argcheck length(prior_model.prior) == size(S, 1)
    D = length(prior_model.prior[1])
    K = size(S, 1)
    sum_S = sum(S, dims=2)
    α = [prior_model.α + sum_S[k] for k = 1:K]
    XS = data'*S'

    a = Vector{Float64}(undef, K)
    b = Vector{Float64}(undef, K)
    for k = 1:K
        a[k] = prior_model.prior[k].α + XS[d, k] for d = 1:D
        b[k] = prior_model.prior[k].θ + sum_S[k]
    end
    BayesianPMM(α, Gamma.(a, b))
end

max_stars_repo_name: SosUts/SuyamaBayes

max_issues_repo_name: SosUts/SuyamaBayes

max_forks_repo_name: SosUts/SuyamaBayes


text:
"""
    SNLI

A corpus of 570k human-written English sentence pairs for NLI.

SNLI sentence pairs are manually labeled with labels entailment, contradiction,
and neutral.

For details, see the [SNLI home page](https://nlp.stanford.edu/projects/snli/)
or read the [2015 paper](https://nlp.stanford.edu/pubs/snli_paper.pdf) "A large annotated corpus for learning natural language inference"
by Samuel R. Bowman, Gabor Angeli, Christopher Potts, and Christopher Manning.

Included data:
```julia
SNLI.train_tsv()
SNLI.train_jsonl()
SNLI.dev_tsv()
SNLI.dev_jsonl()
SNLI.test_tsv()
SNLI.test_jsonl()
```
"""
module SNLI

import ...register_data
using DataDeps: @datadep_str

snli_file(filename) = joinpath(datadep"SNLI", "snli_1.0", filename)

train_tsv()   = snli_file("snli_1.0_train.txt")
train_jsonl() = snli_file("snli_1.0_train.jsonl")
dev_tsv()     = snli_file("snli_1.0_dev.txt")
dev_jsonl()   = snli_file("snli_1.0_dev.jsonl")
test_tsv()    = snli_file("snli_1.0_test.txt")
test_jsonl()  = snli_file("snli_1.0_test.jsonl")

function __init__()
    register_data(
        "SNLI",
        """
        The SNLI corpus (version 1.0) is a collection of 570k human-written
        English sentence pairs manually labeled for balanced classification
        with the labels entailment, contradiction, and neutral, supporting
        the task of natural language inference (NLI), also known as recognizing
        textual entailment (RTE).

        https://nlp.stanford.edu/projects/snli/

        The dataset is distributed as a 100MB zip file.

        The Stanford Natural Language Inference Corpus by The Stanford NLP Group
        is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
        Based on a work at http://shannon.cs.illinois.edu/DenotationGraph/.
        """,
        "https://nlp.stanford.edu/projects/snli/snli_1.0.zip",
        "afb3d70a5af5d8de0d9d81e2637e0fb8c22d1235c2749d83125ca43dab0dbd3e")
end

end # module


text:
# Types

export Error
export error_none, error_invalid_argument, error_system_error,
       error_runtime_error, error_exception
"""
    @enum Error begin
        error_none
        error_invalid_argument
        error_system_error
        error_runtime_error
        error_exception
    end

`Error` return types for all ADIOS2 functions

Based on the [library C++ standardized
exceptions](https://en.cppreference.com/w/cpp/error/exception). Each
error will issue a more detailed description in the standard error
output, stderr
"""
@enum Error begin
    error_none = 0
    error_invalid_argument = 1
    error_system_error = 2
    error_runtime_error = 3
    error_exception = 4
end

@enum AType begin
    type_unknown = -1
    type_string = 0
    type_float = 1
    type_double = 2
    type_float_complex = 3
    type_double_complex = 4
    type_int8_t = 5
    type_int16_t = 6
    type_int32_t = 7
    type_int64_t = 8
    type_uint8_t = 9
    type_uint16_t = 10
    type_uint32_t = 11
    type_uint64_t = 12
    type_long_double = 13
end

# We omit `type_long_double` that we cannot handle
const adios_types = AType[type_string, type_float, type_double,
                          type_float_complex, type_double_complex, type_int8_t,
                          type_int16_t, type_int32_t, type_int64_t,
                          type_uint8_t, type_uint16_t, type_uint32_t,
                          type_uint64_t]
adios_type(::Type{String}) = type_string
adios_type(::Type{Float32}) = type_float
adios_type(::Type{Float64}) = type_double
adios_type(::Type{Complex{Float32}}) = type_float_complex
adios_type(::Type{Complex{Float64}}) = type_double_complex
adios_type(::Type{Int8}) = type_int8_t
adios_type(::Type{Int16}) = type_int16_t
adios_type(::Type{Int32}) = type_int32_t
adios_type(::Type{Int64}) = type_int64_t
adios_type(::Type{UInt8}) = type_uint8_t
adios_type(::Type{UInt16}) = type_uint16_t
adios_type(::Type{UInt32}) = type_uint32_t
adios_type(::Type{UInt64}) = type_uint64_t

const julia_types = Type[String, Float32, Float64, Complex{Float32},
                         Complex{Float64}, Int8, Int16, Int32, Int64, UInt8,
                         UInt16, UInt32, UInt64]
julia_type(type::AType) = julia_types[Int(type) + 1]

export AdiosType
"""
    const AdiosType = Union{AbstractString,
                            Float32,Float64,
                            Complex{Float32},Complex{Float64},
                            Int8,Int16,Int32,Int64,
                            UInt8,UInt16,UInt32,UInt64}

A Union of all scalar types supported in ADIOS files.
"""
const AdiosType = Union{AbstractString,Float32,Float64,Complex{Float32},
                        Complex{Float64},Int8,Int16,Int32,Int64,UInt8,UInt16,
                        UInt32,UInt64}

export Mode
export mode_undefined, mode_write, mode_read, mode_append, mode_deferred,
       mode_sync
"""
    @enum Mode begin
        mode_undefined
        mode_write
        mode_read
        mode_append
        mode_deferred
        mode_sync
    end

Mode specifies for various functions. `write`, `read`, `append` are
used for file operations, `deferred`, `sync` are used for get and put
operations.
"""
@enum Mode begin
    mode_undefined = 0
    mode_write = 1
    mode_read = 2
    mode_append = 3
    mode_deferred = 4
    mode_sync = 5
end

export StepMode
export step_mode_append, step_mode_update, step_mode_read
"""
    @enum StepMode begin
        step_mode_append
        step_mode_update
        step_mode_read
    end
"""
@enum StepMode begin
    step_mode_append = 0
    step_mode_update = 1
    step_mode_read = 2
end

export StepStatus
export step_status_other_error, step_status_ok, step_status_not_ready,
       step_status_end_of_stream
"""
    @enum StepStatus begin
        step_status_other_error
        step_status_ok
        step_status_not_ready
        step_status_end_of_stream
    end
"""
@enum StepStatus begin
    step_status_other_error = -1
    step_status_ok = 0
    step_status_not_ready = 1
    step_status_end_of_stream = 2
end

export ShapeId
export shapeid_unknown, shapeid_global_value, shapeid_global_array,
       shapeid_joined_array, shapeid_local_value, shapeid_local_array
"""
    @enum ShapeId begin
        shapeid_unknown
        shapeid_global_value
        shapeid_global_array
        shapeid_joined_array
        shapeid_local_value
        shapeid_local_array
    end
"""
@enum ShapeId begin
    shapeid_unknown = -1
    shapeid_global_value = 0
    shapeid_global_array = 1
    shapeid_joined_array = 2
    shapeid_local_value = 3
    shapeid_local_array = 4
end
const shapeid_strings = String["shapeid_unknown", "shapeid_global_value",
                               "shapeid_global_array", "shapeid_joined_array",
                               "shapeid_local_value", "shapeid_local_array"]
shapeid_string(shapeid::ShapeId) = shapeid_strings[Int(shapeid) + 2]
Base.show(io::IO, shapeid::ShapeId) = print(io, shapeid_string(shapeid))

const string_array_element_max_size = 4096

const local_value_dim = typemax(Csize_t) - 2
export LocalValue
struct LocalValue end


text:
# Execute this julia file to generate the notebooks from ../notebook.jl

env = @__DIR__
joinpath(env, "..", "..", "generate.jl") |> include
generate(env, execute=false, pluto=false)


text:
################################
## Helper methods for Sampler ##
################################

is_licensor(f :: Feature) = length(f) > 1 && f[1] == '+'
is_licensee(f :: Feature) = length(f) > 1 && f[1] == '-'

###Finds selector features to branch children
function find_features(lex)
    selectors = Feature[]
    licensors = Feature[]
    licensees = Feature[]
    for f in lex.features
        if is_selector(f)
            push!(selectors,f[3:length(f)])
        elseif is_licensor(f)
            push!(licensors,f[2:length(f)])
        elseif is_licensee(f)
            push!(licensees,f[2:length(f)])
        end
    end
    return (selectors, licensors, licensees)
end

#### Recursively builds a sample tree rooted at category C
function treeRootC_move(L, c, parent, prev_licensors)
    move_feat_mismatch = true
    move_feats = []
    id = ""
    selectors = []
    licensors = []
    licensees = []
    while move_feat_mismatch
        move_feats = copy(prev_licensors)
        id = sample(L.MG.lexitem_cond.dists[c])
        lex = L.lexicon[id]
        selectors, licensors, licensees = find_features(lex)
        if licensees != []
            for licensee in licensees
                if licensee in move_feats
                     deleteat!(move_feats ,findfirst(move_feats,licensee))
                else
                     break
                end
                move_feat_mismatch = false
            end
        end
        if selectors == []
            if move_feats != []
                move_feat_mismatch = true
            else
                move_feat_mismatch = false
            end
        else
            move_feat_mismatch = false
        end
    end
    node = TreeNode(id, parent)
    if length(selectors) > 0
        len = length(subroots)
        for sel in subroots[1:(len-1)]
            insert_child!(node, treeRootC(L,sel,node, []))
        end
        insert_child!(node, treeRootC(L,subroots[len],node, append!(move_feats,licensors)))
        return node
    else
        return node
    end
end

###################
##Forward Sampler##
###################

function samplerMGderivation_move(L::SampleGrammar{String, LogProb, SimpleCond})
    i = rand(1:length(L.selectees))
    c = L.selectees[i]
    root = EmptyTree{String}()
    d = treeRootC_move(L,c,root,[])
    prob = probDerivation(d, L)
    return SampleDerivation(d,prob)
end

####################
## Make a corpus  ##
####################

function getchain_fromderiv_move(d, L)
    chain = []
    if ! isterminal(d)
        lex = L.lexicon[d.data]
        selectors, licensors, licensees = find_features(lex)
        for i in 1:(length(selectors))
            f = selectors[i]
            if contains(f, "L")
                chain = append!(getchain_fromderiv_move(d.children[i], L), [lex])
            else
                chain = append!([lex], getchain_fromderiv_move(d.children[i], L))
            end
        end
    else
        lex = L.lexicon[d.data]
        chain = [lex]
    end
    return chain
end

function apply_move(chain)
    sent = []
    while chain != []
        selectors, licensors, licensees = find_features(chain[1])
        while licensors != []
            for i in (2:length(chain))
                next_selectors, next_licensors, next_licensees = find_features(chain[i])
                if next_licensees != []
                    if licensors[1] == next_licensees[1]
                        append!(sent, [chain[i]])
                        deleteat!(licensors,1)
                        deleteat!(chain, i )
                        break
                    end
                end
            end
        end
        push!(sent, chain[1])
        deleteat!(chain, 1)
    end
    return sent
end

function getsent_fromderiv_move(d, L)
    chain = apply_move(getchain_fromderiv_move(d, L))
    sent = ""
    for lex in chain
        sent = string(sent, lex.phon , " ")
    end
    return sent
end

function sample_corpus_move(corpus_size, L, result_dir)
    corpus = Array{Any}(corpus_size)
    for i in 1:corpus_size
        corpus[i] = samplerMGderivation_move(L)
    end
    name1 = string(result_dir,"/corpus.txt")
    name2 = string(result_dir,"/corpus_probs.txt")
    open(name1, "w") do f
        for deriv in corpus
            sent = getsent_fromderiv_move(deriv.d, L)
            write(f, sent)
            write(f, "\n")
        end
    end
    open(name2, "w") do f
        for deriv in corpus
            sent = getsent_fromderiv_move(deriv.d, L)
            prob = deriv.prob
            p = string(exp(prob))
            write(f, sent)
            write(f, string(";", p))
            write(f, "\n")
        end
    end
end


text:
using Juqst

@testset "Random Channel" begin

testChannel = [1.0         0.0           0.0          0.0
               0.00101911  0.971286     -0.0706597   -0.000750742
               0.0126756   0.0711915     0.964136     0.00751068
               -0.0142784   0.000199791  -0.00735239   0.953949 ]

@test round(fidelity(testChannel),digits=10) == 0.9815618333
@test round(unitarity(testChannel),digits=10) == .9310485031
@test round(unitarityPercent(testChannel),digits=10) ==0.0475365736

testChannel = randomFidelityNoise()
@test size(testChannel) == (4,4)
@test iscp(pauliliou2liou(testChannel))
testChannel = randomPrepNoise()
@test size(testChannel) == (4,4)
@test iscp(pauliliou2liou(testChannel))
testChannel = randomMeasureNoise()
@test size(testChannel) == (4,4)
@test iscp(pauliliou2liou(testChannel))

end
@testset "Open Systems" begin

testChannel = [1.0         0.0           0.0          0.0
               0.00101911  0.971286     -0.0706597   -0.000750742
               0.0126756   0.0711915     0.964136     0.00751068
               -0.0142784   0.000199791  -0.00735239   0.953949 ]

t1 = pauliliou2liou(testChannel)
c_s = [0.969835+0.0im         9.98955e-5+0.0036762im   9.98955e-5-0.0036762im    0.0158863+0.0im
       0.000134184+0.0100931im     0.967711+0.0709256im     0.003575+0.0002659im  0.000884926+0.00258246im
       0.000134184-0.0100931im     0.003575-0.0002659im     0.967711-0.0709256im  0.000884926-0.00258246im
       0.0301647+0.0im        -9.98955e-5-0.0036762im  -9.98955e-5+0.0036762im     0.984114+0.0im  ]
@test isapprox(round.(t1,digits=6),round.(c_s,digits=6))
c_s =  [0.969835+0.0im        0.000134184-0.0100931im   9.98955e-5-0.0036762im      0.967711-0.0709256im
 0.000134184+0.0100931im    0.0301647+0.0im           0.003575+0.0002659im   -9.98955e-5+0.0036762im
  9.98955e-5+0.0036762im     0.003575-0.0002659im    0.0158863+0.0im         0.000884926-0.00258246im
    0.967711+0.0709256im  -9.98955e-5-0.0036762im  0.000884926+0.00258246im     0.984114+0.0im ]
@test isapprox(round.(liou2choi(t1),digits=6),round.(c_s,digits=6))
c_s =   [0.969835-0.0im         9.98955e-5-0.0036762im   0.000134184-0.0100931im     0.967711-0.0709256im
         9.98955e-5+0.0036762im    0.0158863-0.0im            0.003575-0.0002659im  0.000884926-0.00258246im
         0.000134184+0.0100931im     0.003575+0.0002659im     0.0301647-0.0im        -9.98955e-5+0.0036762im
         0.967711+0.0709256im  0.000884926+0.00258246im  -9.98955e-5-0.0036762im     0.984114-0.0im ]
@test isapprox(round.(liou2choiX(t1),digits=6),round.(c_s,digits=6))
@test isapprox(round.(testChannel,digits=6),round.(liou2pauliliou(choi2liou(liou2choi(pauliliou2liou(testChannel)))),digits=6))

c_s =   [3.88937+0.0im            0.00101911-0.0148631im   -0.0126756+0.000950533im    -0.0142784+0.141851im
        0.00101911+0.0148631im        0.053201+0.0im         -0.0005318-0.0142784im    -0.000550951+0.0126756im
        -0.0126756-0.000950533im    -0.0005318+0.0142784im     0.038901+0.0im           -0.00015829+0.00101911im
        -0.0142784-0.141851im     -0.000550951-0.0126756im  -0.00015829-0.00101911im       0.018527+0.0im  ]

@test isapprox(round.(c_s,digits=5),round.(choi2chi(liou2choi(pauliliou2liou(testChannel))),digits=5))

end


text:
# Autogenerated wrapper script for Graphviz_jll for x86_64-linux-musl-cxx11
export acyclic, bcomps, ccomps, circo, cluster, dijkstra, dot, dot2gxl, dot_builtins, edgepaint, fdp, gc, gml2gv, graphml2gv, gv2gml, gv2gxl, gvcolor, gvgen, gvmap, gvpack, gvpr, gxl2dot, gxl2gv, libcdt, libcgraph, libgvc, libgvpr, liblab_gamut, libpathplan, libxdot, mm2gv, neato, nop, osage, patchwork, prune, sccmap, sfdp, tred, twopi, unflatten

using Cairo_jll
using Expat_jll
using Pango_jll
JLLWrappers.@generate_wrapper_header("Graphviz")
JLLWrappers.@declare_library_product(libcdt, "libcdt.so.5")
JLLWrappers.@declare_library_product(libcgraph, "libcgraph.so.6")
JLLWrappers.@declare_library_product(libgvc, "libgvc.so.6")
JLLWrappers.@declare_library_product(libgvpr, "libgvpr.so.2")
JLLWrappers.@declare_library_product(liblab_gamut, "liblab_gamut.so.1")
JLLWrappers.@declare_library_product(libpathplan, "libpathplan.so.4")
JLLWrappers.@declare_library_product(libxdot, "libxdot.so.4")
JLLWrappers.@declare_executable_product(acyclic)
JLLWrappers.@declare_executable_product(bcomps)
JLLWrappers.@declare_executable_product(ccomps)
JLLWrappers.@declare_executable_product(circo)
JLLWrappers.@declare_executable_product(cluster)
JLLWrappers.@declare_executable_product(dijkstra)
JLLWrappers.@declare_executable_product(dot)
JLLWrappers.@declare_executable_product(dot2gxl)
JLLWrappers.@declare_executable_product(dot_builtins)
JLLWrappers.@declare_executable_product(edgepaint)
JLLWrappers.@declare_executable_product(fdp)
JLLWrappers.@declare_executable_product(gc)
JLLWrappers.@declare_executable_product(gml2gv)
JLLWrappers.@declare_executable_product(graphml2gv)
JLLWrappers.@declare_executable_product(gv2gml)
JLLWrappers.@declare_executable_product(gv2gxl)
JLLWrappers.@declare_executable_product(gvcolor)
JLLWrappers.@declare_executable_product(gvgen)
JLLWrappers.@declare_executable_product(gvmap)
JLLWrappers.@declare_executable_product(gvpack)
JLLWrappers.@declare_executable_product(gvpr)
JLLWrappers.@declare_executable_product(gxl2dot)
JLLWrappers.@declare_executable_product(gxl2gv)
JLLWrappers.@declare_executable_product(mm2gv)
JLLWrappers.@declare_executable_product(neato)
JLLWrappers.@declare_executable_product(nop)
JLLWrappers.@declare_executable_product(osage)
JLLWrappers.@declare_executable_product(patchwork)
JLLWrappers.@declare_executable_product(prune)
JLLWrappers.@declare_executable_product(sccmap)
JLLWrappers.@declare_executable_product(sfdp)
JLLWrappers.@declare_executable_product(tred)
JLLWrappers.@declare_executable_product(twopi)
JLLWrappers.@declare_executable_product(unflatten)
function __init__()
    JLLWrappers.@generate_init_header(Cairo_jll, Expat_jll, Pango_jll)
    JLLWrappers.@init_library_product(
        libcdt,
        "lib/libcdt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcgraph,
        "lib/libcgraph.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgvc,
        "lib/libgvc.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libgvpr,
        "lib/libgvpr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        liblab_gamut,
        "lib/liblab_gamut.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libpathplan,
        "lib/libpathplan.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libxdot,
        "lib/libxdot.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        acyclic,
        "bin/acyclic",
    )

    JLLWrappers.@init_executable_product(
        bcomps,
        "bin/bcomps",
    )

    JLLWrappers.@init_executable_product(
        ccomps,
        "bin/ccomps",
    )

    JLLWrappers.@init_executable_product(
        circo,
        "bin/circo",
    )

    JLLWrappers.@init_executable_product(
        cluster,
        "bin/cluster",
    )

    JLLWrappers.@init_executable_product(
        dijkstra,
        "bin/dijkstra",
    )

    JLLWrappers.@init_executable_product(
        dot,
        "bin/dot",
    )

    JLLWrappers.@init_executable_product(
        dot2gxl,
        "bin/dot2gxl",
    )

    JLLWrappers.@init_executable_product(
        dot_builtins,
        "bin/dot_builtins",
    )

    JLLWrappers.@init_executable_product(
        edgepaint,
        "bin/edgepaint",
    )

    JLLWrappers.@init_executable_product(
        fdp,
        "bin/fdp",
    )

    JLLWrappers.@init_executable_product(
        gc,
        "bin/gc",
    )

    JLLWrappers.@init_executable_product(
        gml2gv,
        "bin/gml2gv",
    )

    JLLWrappers.@init_executable_product(
        graphml2gv,
        "bin/graphml2gv",
    )

    JLLWrappers.@init_executable_product(
        gv2gml,
        "bin/gv2gml",
    )

    JLLWrappers.@init_executable_product(
        gv2gxl,
        "bin/gv2gxl",
    )

    JLLWrappers.@init_executable_product(
        gvcolor,
        "bin/gvcolor",
    )

    JLLWrappers.@init_executable_product(
        gvgen,
        "bin/gvgen",
    )

    JLLWrappers.@init_executable_product(
        gvmap,
        "bin/gvmap",
    )

    JLLWrappers.@init_executable_product(
        gvpack,
        "bin/gvpack",
    )

    JLLWrappers.@init_executable_product(
        gvpr,
        "bin/gvpr",
    )

    JLLWrappers.@init_executable_product(
        gxl2dot,
        "bin/gxl2dot",
    )

    JLLWrappers.@init_executable_product(
        gxl2gv,
        "bin/gxl2gv",
    )

    JLLWrappers.@init_executable_product(
        mm2gv,
        "bin/mm2gv",
    )

    JLLWrappers.@init_executable_product(
        neato,
        "bin/neato",
    )

    JLLWrappers.@init_executable_product(
        nop,
        "bin/nop",
    )

    JLLWrappers.@init_executable_product(
        osage,
        "bin/osage",
    )

    JLLWrappers.@init_executable_product(
        patchwork,
        "bin/patchwork",
    )

    JLLWrappers.@init_executable_product(
        prune,
        "bin/prune",
    )

    JLLWrappers.@init_executable_product(
        sccmap,
        "bin/sccmap",
    )

    JLLWrappers.@init_executable_product(
        sfdp,
        "bin/sfdp",
    )

    JLLWrappers.@init_executable_product(
        tred,
        "bin/tred",
    )

    JLLWrappers.@init_executable_product(
        twopi,
        "bin/twopi",
    )

    JLLWrappers.@init_executable_product(
        unflatten,
        "bin/unflatten",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using ITensors
using ITensorGaussianMPS
using LinearAlgebra

# Half filling
N = 20
Nf_up = N ÷ 2
Nf_dn = N ÷ 2
Nf = Nf_up + Nf_dn

@show N, Nf

# Maximum MPS link dimension
_maxlinkdim = 50

@show _maxlinkdim

# DMRG cutoff
_cutoff = 1e-8

# Hopping
t = 1.0

# Electron-electron on-site interaction
U = 1.0

@show t, U

# Make the free fermion Hamiltonian for the up spins
os_up = OpSum()
for n in 1:(N - 1)
  os_up .+= -t, "Cdagup", n, "Cup", n + 1
  os_up .+= -t, "Cdagup", n + 1, "Cup", n
end

# Make the free fermion Hamiltonian for the down spins
os_dn = OpSum()
for n in 1:(N - 1)
  os_dn .+= -t, "Cdagdn", n, "Cdn", n + 1
  os_dn .+= -t, "Cdagdn", n + 1, "Cdn", n
end

# Hopping Hamiltonians for the up and down spins
h_up = hopping_hamiltonian(os_up)
h_dn = hopping_hamiltonian(os_dn)

# Get the Slater determinant
Φ_up = slater_determinant_matrix(h_up, Nf_up)
Φ_dn = slater_determinant_matrix(h_dn, Nf_dn)

# Create an MPS from the slater determinants.
s = siteinds("Electron", N; conserve_qns=true)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(
  s, Φ_up, Φ_dn; eigval_cutoff=1e-4, cutoff=_cutoff, maxdim=_maxlinkdim
)
@show maxlinkdim(ψ0)

# The total interacting Hamiltonian
os = os_up + os_dn
for n in 1:N
  os .+= U, "Nupdn", n
end
H = MPO(os, s)

println("Free fermion starting state energy")
@show flux(ψ0)
@show inner(ψ0, H, ψ0)

println("\nStart from free fermion state")
sweeps = Sweeps(5)
maxdim!(sweeps, _maxlinkdim)
cutoff!(sweeps, _cutoff)
e, ψ = @time dmrg(H, ψ0, sweeps)
@show e
@show flux(ψ)

using ITensorGaussianMPS: correlation_matrix_to_gmps, correlation_matrix_to_mps, entropy

Λ_up = correlation_matrix(ψ, "Cdagup", "Cup")
Λ_dn = correlation_matrix(ψ, "Cdagdn", "Cdn")
ψ̃0 = correlation_matrix_to_mps(
  s, Λ_up, Λ_dn; eigval_cutoff=1e-2, maxblocksize=4
)
@show inner(ψ̃0, ψ)
@show inner(ψ̃0, H, ψ̃0)


text:
# utility/z-test.jl

using Test: @test

@test downsample(:test)
@test interp1(:test)
@test jinc(:test)
@test map_many(:test)
@test mask_test()
@test max_percent_diff(:test)
@test ndgrid(:test)
@test rect(:test)
@test eql_root(:test)
@test rmsd100(:test)


text:
# Autogenerated wrapper script for Bison_jll for x86_64-apple-darwin
export bison

using Libiconv_jll
JLLWrappers.@generate_wrapper_header("Bison")
JLLWrappers.@declare_executable_product(bison)
function __init__()
    JLLWrappers.@generate_init_header(Libiconv_jll)
    JLLWrappers.@init_executable_product(
        bison,
        "bin/bison",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    TrajectoryOptimization
Primary module for setting up and solving trajectory optimization problems with
iterative Linear Quadratic Regulator (iLQR). Module supports unconstrained and
constrained optimization problems. Constrained optimization problems are solved
using Augmented Lagrangian methods. Supports automatic differentiation via the
`ForwardDiff` package.
"""
module TrajectoryOptimization

using RigidBodyDynamics
using ForwardDiff
using JuMP
using DocStringExtensions
using Interpolations
using RecipesBase
using LinearAlgebra
using Statistics
using Random
using SparseArrays
using StaticArrays
using Logging
using Formatting
using Plots
using BenchmarkTools

export
    Dynamics

# Primary types
export
    Model,
    Solver,
    SolverResults,
    ConstrainedObjective,
    UnconstrainedObjective,
    LQRObjective,
    QuadraticCost,
    LQRCost,
    GenericCost,
    ConstrainedVectorResults,
    UnconstrainedVectorResults,
    SolverOptions,
    Trajectory

# Primary methods
export
    solve,
    rollout!,
    rollout,
    forwardpass!,
    backwardpass!,
    cost,
    max_violation,
    update_objective,
    infeasible_control,
    line_trajectory

export
    get_sizes,
    get_num_constraints,
    get_num_controls,
    init_results,
    to_array,
    get_N,
    quat2rot,
    sphere_constraint,
    circle_constraint,
    plot_trajectory!,
    plot_vertical_lines!,
    convergence_rate,
    plot_obstacles,
    generate_controller,
    lqr,
    evals,
    reset,
    reset_evals

include("objective.jl")
include("model.jl")
include("integration.jl")
include("solver.jl")
include("results.jl")
include("results_dircol.jl")
include("backwardpass.jl")
include("forwardpass.jl")
include("constraints.jl")
include("rollout.jl")
#include("newton.jl")
include("infeasible.jl")
include("minimum_time.jl")
include("ilqr_methods.jl")
include("augmented_lagrangian.jl")
include("solve.jl")
include("utils.jl")
# include("dynamics.jl")
include("logger.jl")
include("controller.jl")
include("partitioning.jl")

using Ipopt

# DIRCOL methods
export
solve_dircol,
gen_usrfun,
DircolResults,
DircolVars,
collocation_constraints,
collocation_constraints!,
cost_gradient,
cost_gradient!,
constraint_jacobian,
constraint_jacobian!,
get_weights,
get_initial_state

export
packZ,
unpackZ

include("dircol.jl")
include("dircol_ipopt.jl")
write_ipopt_options()

# if "Snopt" in keys(Pkg.installed())
#     using Snopt # not safe for precompilation
#     include("dircol_snopt.jl")
# end

end


text:
import Pkg
Pkg.activate("RbfInterpolationTools")
using RbfInterpolationTools
using Plots

# Lets generate some observed points
X = hcat(rand(30) .* 6.28)
# We will try to interpolate the sin function
y = vec(map(x -> sin(x), X))

# Lets create our interpolation environment
inter = RbfInterpolationTools.Interpolator(X = X, y = (y))
# Lets regenerate our interpolation domain
X_new = hcat(0:6.28/(999):6.28)

# And we simply get our answer:
y_new = RbfInterpolationTools.interpolateMesh(inter, X_new)

# Lets plot the results
plot()
Plots.scatter!(X, y, label = "Shown points")
Plots.plot!(X_new, y_new, label = "Interpolation")

# The derivative shall be the cosine function. Lets verify it:

∂y = vec(map(x -> RbfInterpolationTools.gradientPoint(inter, [x])[1], X_new))
Plots.plot!(X_new, ∂y, label = "Interpolated Derivative")

# We can actually compare it with the real values
Plots.plot!(X_new, vec(map(x -> cos(x), X_new)), label = "Real Derivative")


text:
using NumericIO
using Test

#==Generate test code
===============================================================================#
function _gentestcode()
	tf = open("testgen.jl", "w")
	_ndigits=5
	println(tf, "SI(v) = formatted(v, :SI, ndigits=$_ndigits)")

	for _exp in -30:30
		v = e*10^Float64(_exp)
		vstr = formatted(v, :SI, ndigits=_ndigits)
		println(tf, "@test SI($v) == \"$vstr\"")
	end

	for _ndigits in 0:5
		println(tf, "SI(v) = formatted(v, :SI, ndigits=$_ndigits)")

		for _exp in 6:9
			v = e*10^Float64(_exp)
			vstr = formatted(v, :SI, ndigits=_ndigits)
			println(tf, "@test SI($v) == \"$vstr\"")
		end
	end

	_exp = 6
	v = e*10^Float64(_exp)
	for _ndigits in 0:5
		println(tf, "dflt(v) = formatted(v, ndigits=$_ndigits)")
		vstr = formatted(v, ndigits=_ndigits)
		println(tf, "@test dflt($v) == \"$vstr\"")
	end
	close(tf)
end
#_gentestcode()

#==Helper functions
===============================================================================#
newscope(fn::Function) = (fn()) #Allow re-definition of SI(v) / dflt(v)

#==Exponent tests: Formatting with SI prefixes
===============================================================================#
newscope() do
SI(v) = formatted(v, :SI, ndigits=5)
@test SI(-2.7182818284590455e-30) == "−2.7183×10⁻³⁰" #Negative
@test SI(-2.7182818284590453e-24) == "−2.7183y"
@test SI(-2.718281828459045e30) == "−2.7183×10³⁰"

@test SI(2.7182818284590455e-30) == "2.7183×10⁻³⁰"
@test SI(2.7182818284590454e-29) == "27.183×10⁻³⁰"
@test SI(2.718281828459045e-28) == "271.83×10⁻³⁰"
@test SI(2.7182818284590452e-27) == "2.7183×10⁻²⁷"
@test SI(2.718281828459045e-26) == "27.183×10⁻²⁷"
@test SI(2.718281828459045e-25) == "271.83×10⁻²⁷"
@test SI(2.7182818284590453e-24) == "2.7183y"
@test SI(2.718281828459045e-23) == "27.183y"
@test SI(2.718281828459045e-22) == "271.83y"
@test SI(2.7182818284590455e-21) == "2.7183z"
@test SI(2.7182818284590455e-20) == "27.183z"
@test SI(2.718281828459045e-19) == "271.83z"
@test SI(2.7182818284590453e-18) == "2.7183a"
@test SI(2.718281828459045e-17) == "27.183a"
@test SI(2.718281828459045e-16) == "271.83a"
@test SI(2.718281828459045e-15) == "2.7183f"
@test SI(2.7182818284590452e-14) == "27.183f"
@test SI(2.718281828459045e-13) == "271.83f"
@test SI(2.718281828459045e-12) == "2.7183p"
@test SI(2.7182818284590454e-11) == "27.183p"
@test SI(2.718281828459045e-10) == "271.83p"
@test SI(2.718281828459045e-9) == "2.7183n"
@test SI(2.718281828459045e-8) == "27.183n"
@test SI(2.718281828459045e-7) == "271.83n"
@test SI(2.718281828459045e-6) == "2.7183μ"
@test SI(2.718281828459045e-5) == "27.183μ"
@test SI(0.00027182818284590454) == "271.83μ"
@test SI(0.002718281828459045) == "2.7183m"
@test SI(0.027182818284590453) == "27.183m"
@test SI(0.27182818284590454) == "271.83m"
@test SI(2.718281828459045) == "2.7183"
@test SI(27.18281828459045) == "27.183"
@test SI(271.8281828459045) == "271.83"
@test SI(2718.2818284590453) == "2.7183k"
@test SI(27182.818284590452) == "27.183k"
@test SI(271828.1828459045) == "271.83k"
@test SI(2.718281828459045e6) == "2.7183M"
@test SI(2.718281828459045e7) == "27.183M"
@test SI(2.718281828459045e8) == "271.83M"
@test SI(2.718281828459045e9) == "2.7183G"
@test SI(2.718281828459045e10) == "27.183G"
@test SI(2.718281828459045e11) == "271.83G"
@test SI(2.718281828459045e12) == "2.7183T"
@test SI(2.718281828459045e13) == "27.183T"
@test SI(2.718281828459045e14) == "271.83T"
@test SI(2.718281828459045e15) == "2.7183P"
@test SI(2.7182818284590452e16) == "27.183P"
@test SI(2.718281828459045e17) == "271.83P"
@test SI(2.718281828459045e18) == "2.7183E"
@test SI(2.7182818284590453e19) == "27.183E"
@test SI(2.7182818284590452e20) == "271.83E"
@test SI(2.718281828459045e21) == "2.7183Z"
@test SI(2.718281828459045e22) == "27.183Z"
@test SI(2.7182818284590448e23) == "271.83Z"
@test SI(2.718281828459045e24) == "2.7183Y"
@test SI(2.7182818284590455e25) == "27.183Y"
@test SI(2.7182818284590452e26) == "271.83Y"
@test SI(2.7182818284590453e27) == "2.7183×10²⁷"
@test SI(2.718281828459045e28) == "27.183×10²⁷"
@test SI(2.7182818284590454e29) == "271.83×10²⁷"
@test SI(2.718281828459045e30) == "2.7183×10³⁰"
end

#==Number of digits: Formatting with SI prefixes
===============================================================================#
newscope() do
SI(v) = formatted(v, :SI, ndigits=0)
@test SI(2.718281828459045e6) == "2.718281828459045M"
@test SI(2.718281828459045e7) == "27.18281828459045M"
@test SI(2.718281828459045e8) == "271.8281828459045M"
@test SI(2.718281828459045e9) == "2.718281828459045G"
end
newscope() do
SI(v) = formatted(v, :SI, ndigits=1)
@test SI(2.718281828459045e6) == "3M"
@test SI(2.718281828459045e7) == "30M"
@test SI(2.718281828459045e8) == "300M"
@test SI(2.718281828459045e9) == "3G"
end
newscope() do
SI(v) = formatted(v, :SI, ndigits=2)
@test SI(2.718281828459045e6) == "2.7M"
@test SI(2.718281828459045e7) == "27M"
@test SI(2.718281828459045e8) == "270M"
@test SI(2.718281828459045e9) == "2.7G"
end
newscope() do
SI(v) = formatted(v, :SI, ndigits=3)
@test SI(2.718281828459045e6) == "2.72M"
@test SI(2.718281828459045e7) == "27.2M"
@test SI(2.718281828459045e8) == "272M"
@test SI(2.718281828459045e9) == "2.72G"
end
newscope() do
SI(v) = formatted(v, :SI, ndigits=4)
@test SI(2.718281828459045e6) == "2.718M"
@test SI(2.718281828459045e7) == "27.18M"
@test SI(2.718281828459045e8) == "271.8M"
@test SI(2.718281828459045e9) == "2.718G"
end
newscope() do
SI(v) = formatted(v, :SI, ndigits=5)
@test SI(2.718281828459045e6) == "2.7183M"
@test SI(2.718281828459045e7) == "27.183M"
@test SI(2.718281828459045e8) == "271.83M"
@test SI(2.718281828459045e9) == "2.7183G"
end

#==Number of digits: Formatting with scientific notation
===============================================================================#
newscope() do
dflt(v) = formatted(v, ndigits=0)
@test dflt(2.718281828459045e6) == "2.718281828459045×10⁶"
end
newscope() do
dflt(v) = formatted(v, ndigits=1)
@test dflt(2.718281828459045e6) == "3×10⁶"
end
newscope() do
dflt(v) = formatted(v, ndigits=2)
@test dflt(2.718281828459045e6) == "2.7×10⁶"
end
newscope() do
dflt(v) = formatted(v, ndigits=3)
@test dflt(2.718281828459045e6) == "2.72×10⁶"
end
newscope() do
dflt(v) = formatted(v, ndigits=4)
@test dflt(2.718281828459045e6) == "2.718×10⁶"
end
newscope() do
dflt(v) = formatted(v, ndigits=5)
@test dflt(2.718281828459045e6) == "2.7183×10⁶"
end

#==Integer values: Formatting
===============================================================================#
@test formatted(0, :SI, ndigits=0) == "0" #Ensure not "0.00000000..."
@test formatted(27_182_818_284, :SI, ndigits=3) == "27.2G"
@test formatted(27_182_818_284, :ENG, ndigits=3) == "27.2×10⁹"
@test formatted(27_182_818_284, :SCI, ndigits=3) == "2.72×10¹⁰"
@test formatted(27_182_818_284, :SCI, ndigits=3, charset=:ASCII) == "2.72E10"

#==Max ndigits, fixed decpos
===============================================================================#
fmt = NumericIO.IOFormattingReal(NumericIO.UEXPONENT_SI,
	ndigits=0, decpos=9, decfloating=false, eng=true,
	minus=NumericIO.UTF8_MINUS_SYMBOL, inf=NumericIO.UTF8_INF_STRING
)
newscope() do
SI(v) = formatted(v, fmt)
@test SI(0) == "0G" #Ensure not "0.00000000...G"
@test SI(2.718281828459045e5) == "0.0002718281828459045G"
@test SI(2.718281828459045e6) == "0.002718281828459045G"
@test SI(2.718281828459045e7) == "0.02718281828459045G"
@test SI(2.718281828459045e8) == "0.2718281828459045G"
@test SI(2.718281828459045e9) == "2.718281828459045G"
@test SI(2.718281828459045e10) == "27.18281828459045G"
@test SI(2.718281828459045e11) == "271.8281828459045G"
@test SI(2.718281828459045e12) == "2718.281828459045G"
@test SI(2.718281828459045e13) == "27182.81828459045G"
end

#==Fixed ndigits, fixed decpos
===============================================================================#
#Use lower-level structures for fixed decimal position:
fmt = NumericIO.IOFormattingReal(NumericIO.UEXPONENT_SI,
	ndigits=4, decpos=9, decfloating=false, eng=true,
	minus=NumericIO.UTF8_MINUS_SYMBOL, inf=NumericIO.UTF8_INF_STRING
)
newscope() do
SI(v) = formatted(v, fmt)
@test SI(0.75e9) == "0.750G"
@test SI(2.718281828459045e5) == "0.000G"
@test SI(2.718281828459045e6) == "0.003G"
@test SI(2.718281828459045e7) == "0.027G"
@test SI(2.718281828459045e8) == "0.272G"
@test SI(2.718281828459045e9) == "2.718G"
@test SI(2.718281828459045e10) == "27.18G"
@test SI(2.718281828459045e11) == "271.8G"
@test SI(2.718281828459045e12) == "2718G"
@test SI(2.718281828459045e13) == "27180G"
end

:Tests_Complete


text:
using Plasmo
using Plots
graph = OptiGraph()

sg1 = OptiGraph()
sg2 = OptiGraph()
add_subgraph!(graph,sg1)
add_subgraph!(graph,sg2)

@optinode(sg1,ng1[1:5])
@optinode(sg2,ng2[1:5])

for node in getnodes(sg1)
    @variable(node,0 <= x <= 2)
    @objective(node,Max,x)
end

for node in getnodes(sg2)
    @variable(node,x>=2)
    @objective(node,Min,x)
end

@linkconstraint(sg1,sum(ng1[i][:x] for i = 1:5) <= 4)
@linkconstraint(sg2,sum(ng2[i][:x] for i = 1:5) >= 4)
@linkconstraint(graph,ng1[3][:x]==ng2[3][:x])

plt_G = plt_G = Plots.plot(graph,node_labels = true, markersize=60, labelsize=30, linewidth=4, layout_options = Dict(:tol => 0.01, :iterations => 2))
savefig("004.png")

max_stars_repo_name: BoyangSheng/Skill-studing

max_issues_repo_name: BoyangSheng/Skill-studying

max_forks_repo_name: BoyangSheng/Skill-studying


text:
# OPTIMO_AUGLAG
#
# PART OF OPTIMO

export AugLagOptiModel, AugLagUpdate!

"""
Let an AbstractOptiModel `prob` represent the problem
   optimize     f(x) + g(x)
   subject to   c(x) in S

Given an AbstractVector `mu` of positive penalty parameters and
an AbstractVector `y` of Lagrange multipliers,

   AugLagOptiModel(prob, mu, y)

returns an AugLagOptiModel ´alprob´ which represents the problem

   optimize     L(X) + g(X)

where
  L(x) = f(x) + 0.5 * sum_i [(c_i(x) - p_i(x))^2 / mu_i] - 0.5 * sum_i [mu_i * y_i^2]
  p(x) = proj_S( c(x) )

Note that both ´mu´ and `y` are vectors of dimension `prob.meta.ncon`.
"""
mutable struct AugLagOptiModel{R <: Real} <: AbstractOptiModel
    meta::OptiModelMeta
    base::AbstractOptiModel             # original problem
    mu::AbstractVector{R}               # penalty parameter
    y::AbstractVector{R}                # dual variable estimate
    # additional storage
    _muy::AbstractVector{R}              # mu .* y
    _musqy::R                            # 0.5 * sum( mu .* y.^2 )
    _cx::AbstractVector{R}               # c(x)
    _wz::AbstractVector{R}               # w - z
    _yu::AbstractVector{R}               # yu = ( w - z )./mu

    function AugLagOptiModel{R}( meta::OptiModelMeta,
                                 prob::TP,
                                 mu::Tx,
                                 y::Tx,
                        ) where {R <: Real, TP <: AbstractOptiModel, Tx <: AbstractVector{R}}
        ny = prob.meta.ncon
        muy = mu .* y
        musqy = 0.5 * sum( muy .* y )
        new( meta, prob, mu, y, muy, musqy, zeros(R, ny), zeros(R, ny), zeros(R, ny) )
    end
end

function AugLagOptiModel{R}( prob::TP,
                             mu::Tx,
                             y::Tx;
                             name=prob.meta.name * "-al",
                        ) where {R <: Real, TP <: AbstractOptiModel, Tx <: AbstractVector{R}}
    @assert prob.meta.ncon > 0
    @lencheck prob.meta.ncon mu y
    @assert all( mu .> 0 )
    meta = OptiModelMeta( prob.meta.nvar, 0, x0=prob.meta.x0, minimize=prob.meta.minimize, name=name )
    al = AugLagOptiModel{R}( meta, prob, mu, y )
    finalizer( prob -> finalize(prob.base), al )
    return al
end

AugLagOptiModel( prob, mu, y, ::Type{R}; kwargs...) where {R} =
    AugLagOptiModel{R}( prob, mu, y; kwargs...)
AugLagOptiModel( prob, mu, y; kwargs...) =
    AugLagOptiModel( prob, mu, y, Float64; kwargs...)

"""
    AugLagUpdate!( al, mu, y )
"""
function AugLagUpdate!( al::TA, mu::Tx, y::Tx,
            ) where {R <: Real, Tx <: AbstractVector{R}, TA <: AugLagOptiModel{R}}
    @lencheck al.base.meta.ncon mu y
    @assert all( mu .> 0 )
    al.mu = mu
    al.y  = y
    al._muy = al.mu .* al.y
    al._musqy = 0.5 * sum( al._muy .* al.y )
    return nothing
end

"""
    lx = obj( al, x )
"""
function obj( al::AugLagOptiModel, x::AbstractVector )
    @lencheck al.meta.nvar x
    lx = objcons!( al.base, x, al._cx )     # fx, cx
    al._wz .= al._cx .+ al._muy             # w = cx + mu .* y
    proj!( al.base, al._wz, al._yu )        # z = proj( cx + mu .* y )
    al._wz .-= al._yu                       # w - z
    al._yu .= al._wz ./ al.mu               # yu = (w - z) ./ mu
    lx += 0.5 * sum( al._wz .* al._yu )     # fx + 1/(2 mu) dist[S]^2( w )
    lx -= al._musqy
    return lx
end

"""
    grad!( al, x, dlx )
"""
function grad!( al::AugLagOptiModel, x::AbstractVector, dlx::AbstractVector )
    @lencheck al.meta.nvar x dlx
    cons!( al.base, x, al._cx )             # cx
    al._wz .= al._cx .+ al._muy             # w = cx + mu .* y
    proj!( al.base, al._wz, al._yu )        # z = proj( cx + mu .* y )
    al._wz .-= al._yu                       # w - z
    al._yu .= al._wz ./ al.mu               # yu = (w - z) ./ mu
    grad!( al.base, x, dlx )                # dfx
    dlx .+= jtprod( al.base, x, al._yu )    # dfx + dcx' * yu
    return nothing
end

"""
    lx = objgrad!( al, x, dlx )
"""
function objgrad!( al::AugLagOptiModel, x::AbstractVector, dlx::AbstractVector )
    @lencheck al.meta.nvar x dlx
    cons!( al.base, x, al._cx )             # cx
    al._wz .= al._cx .+ al._muy             # w = cx + mu .* y
    proj!( al.base, al._wz, al._yu )        # z = proj( cx + mu .* y )
    al._wz .-= al._yu                       # w - z
    al._yu .= al._wz ./ al.mu               # yu = (w - z) ./ mu
    lx = objgrad!( al.base, x, dlx )        # fx, dfx
    lx += 0.5 * sum( al._wz .* al._yu )     # fx + 1/(2 mu) dist[S]^2( w )
    lx -= al._musqy
    dlx .+= jtprod( al.base, x, al._yu )    # dfx + dcx' * yu
    return lx
end

"""
    prox!( al, x, a, z )
"""
function prox!(al::AugLagOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
    @lencheck al.meta.nvar x z
    prox!( al.base, x, a, z )
    return nothing
end

"""
    gz = objprox!( al, x, a, z )
"""
function objprox!(al::AugLagOptiModel, x::AbstractVector, a::Real, z::AbstractVector)
    @lencheck al.meta.nvar x z
    gz = objprox!( al.base, x, a, z )
    return gz
end


text:
function parse_file(cur_file::AbstractString)

  file_shards = Array{Expr}(0)

  opened_file = open(cur_file)

  cur_expression = parse(opened_file)

  while cur_expression != nothing
    isa(cur_expression, Number) &&
      ( cur_expression = :( cur_expression + 0 ) )

    push!(file_shards, cur_expression)
    cur_expression = parse(opened_file)
  end

  close(opened_file)

  file_shards

end


text:
@testset "297.serialize-and-deserialize-binary-tree.jl" begin
    root = TreeNode{Int}([1, 2, 3, nothing, nothing, 4, 5])
    @test deserialize(serialize(root)) == root
    @test deserialize(serialize(nothing)) === nothing
    @test deserialize(serialize(TreeNode(1))) == TreeNode(1)
    @test deserialize(serialize(TreeNode{Int}([1, 2]))) == TreeNode{Int}([1, 2])
    @test deserialize(serialize(TreeNode{Int}([3, 2, 4, 3]))) == TreeNode{Int}([3, 2, 4, 3])
end


text:
"""
    varN(BinningAnalysis[, lvl = 0])

Calculates the variance/N of a given level in the Binning Analysis.
"""
function varN(B::LogBinner, lvl::Int64 = 0)

    n = B.count[lvl+1]
    var(B, lvl) / n
end

"""
    var(BinningAnalysis[, lvl = 0])

Calculates the variance of a given level in the Binning Analysis.
"""
function var(
        B::LogBinner{N, T},
        lvl::Int64 = 0
    ) where {N, T <: Real}

    n = B.count[lvl+1]
    X = B.x_sum[lvl+1]
    X2 = B.x2_sum[lvl+1]

    # lvl = 1 <=> original values
    # correct variance:
    # (∑ xᵢ^2) / (N-1) - (∑ xᵢ)(∑ xᵢ) / (N(N-1))
    X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(
        B::LogBinner{N, T},
        lvl::Int64 = 0
    ) where {N, T <: Complex}

    n = B.count[lvl+1]
    X = B.x_sum[lvl+1]
    X2 = B.x2_sum[lvl+1]

    # lvl = 1 <=> original values
    (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

function var(
        B::LogBinner{N, <: AbstractArray{T, D}},
        lvl::Int64 = 0
    ) where {N, D, T <: Real}

    n = B.count[lvl+1]
    X = B.x_sum[lvl+1]
    X2 = B.x2_sum[lvl+1]

    @. X2 / (n - 1) - X^2 / (n*(n - 1))
end

function var(
        B::LogBinner{N, <: AbstractArray{T, D}},
        lvl::Int64 = 0
    ) where {N, D, T <: Complex}

    n = B.count[lvl+1]
    X = B.x_sum[lvl+1]
    X2 = B.x2_sum[lvl+1]

    @. (real(X2) + imag(X2)) / (n - 1) - (real(X)^2 + imag(X)^2) / (n*(n - 1))
end

"""
    all_vars(BinningAnalysis)

Calculates the variance for each level of the Binning Analysis.
"""
function all_vars(B::LogBinner{N}) where {N}
    [var(B, lvl) for lvl in 0:N-1 if B.count[lvl+1] > 1]
end

"""
    all_varNs(BinningAnalysis)

Calculates the variance/N for each level of the Binning Analysis.
"""
function all_varNs(B::LogBinner{N}) where {N}
    [varN(B, lvl) for lvl in 0:N-1 if B.count[lvl+1] > 1]
end

################################################################################

# NOTE works for all
"""
    mean(BinningAnalysis[, lvl=0])

Calculates the mean for a given level in the Binning Analysis.
"""
function mean(B::LogBinner, lvl::Int64 = 0)
    B.x_sum[lvl+1] / B.count[lvl+1]
end

# NOTE works for all
"""
    all_means(BinningAnalysis)

Calculates the mean for each level of the Binning Analysis.
"""
function all_means(B::LogBinner{N}) where {N}
    [mean(B, lvl) for lvl in 0:N-1 if B.count[lvl+1] > 1]
end

################################################################################

"""
    tau(BinningAnalysis[, lvl=0])

Calculates the autocorrelation time tau for a given binning level.
"""
function tau(B::LogBinner, lvl::Int64 = 0)
    var_0 = varN(B, 0)
    var_l = varN(B, lvl)
    0.5 * (var_l / var_0 - 1)
end

"""
    all_taus(BinningAnalysis)

Calculates the autocorrelation time tau for each level of the Binning Analysis.
"""
function all_taus(B::LogBinner{N}) where {N}
    [tau(B, lvl) for lvl in 0:N-1 if B.count[lvl+1] > 1]
end

################################################################################

"""
    std_error(BinningAnalysis[, lvl=0])

Calculates the standard error for a given level.
"""
function std_error(B::LogBinner{N, T}, lvl::Int64=0) where {N, T <: Number}
    sqrt(varN(B, lvl))
end
function std_error(B::LogBinner{N, T}, lvl::Int64=0) where {N, T <: AbstractArray}
    sqrt.(varN(B, lvl))
end

"""
    all_std_errors(BinningAnalysis)

Calculates the standard error for each level of the Binning Analysis.
"""
function all_std_errors(B::LogBinner{N, T}) where {N, T <: Number}
    sqrt.(all_varNs(B))
end
function all_std_errors(B::LogBinner{N, T}) where {N, T <: AbstractArray}
    (x -> sqrt.(x)).all_varNs(B)
end

"""
    convergence(BinningAnalysis, lvl)

Computes the difference between the variance of this lvl and the last,
normalized to the last lvl. If this value tends to 0, the Binning Analysis has
converged.
"""
function convergence(B::LogBinner{N, T}, lvl::Int64) where {N, T <: Number}
    abs((varN(B, lvl) - varN(B, lvl-1)) / varN(B, lvl-1))
end
function convergence(B::LogBinner{N, T}, lvl::Int64) where {N, T <: AbstractArray}
    mean(abs.((varN(B, lvl) .- varN(B, lvl-1)) ./ varN(B, lvl-1)))
end

"""
    has_converged(BinningAnalysis, lvl[, threshhold = 0.05])

Returns true if the Binning Analysis has converged for a given lvl.
"""
function has_converged(B::LogBinner, lvl::Int64, threshhold::Float64 = 0.05)
    convergence(B, lvl) <= threshhold
end


text:
@testset "solve_hadc.jl" begin

end


text:
include("test_database.jl")
include("test_models.jl")
include("test_methods.jl")


text:
module OPFInstances

function getinstancepath(archive::String, subfolder::String, instance::String)
    relpath = ""
    if archive == "Matpower"
        if subfolder == "m"
            relpath = joinpath("instances", "data_"*archive, "matpower", instance*".m")
        elseif subfolder == "QCQP"
            relpath = joinpath("instances", "data_"*archive, "matpower_QCQP", instance*".dat")
        else
            error("getinstancepath(): incorrect matpower instance kind (should be either `m` or `QCQP`).")
        end

    elseif archive == "GOC"
        relpath = joinpath("instances", "data_"*archive, subfolder, instance)
    else
        error("getinstancepath(): incorrect instance archive name.")
    end

    return joinpath(Pkg.dir("OPFInstances"), relpath)
end

export getinstancepath

end