text:
type D2Function
    value::Float64
    gradient::Array{Float64,1}
    hessian::Array{Float64,2}
end

###{{{ Newton-Raphson

function NR(f, x; iter=200,tol=1e-12,verbose=false,method=0,gamma=1,trace=0)
    val = f(x)::D2Function
    # S = val.grad["grad"]';  H = val["hess"]; loglik=val["val"]
    S = val.gradient::Vector{Float64};  H = val.hessian::Matrix{Float64}; loglik=val.value::Float64
    p = size(x,1)
    counter = 1
    for j=1:iter
        if verbose print("S=",S') end
        s = s0 = mean(S.*S)
        if isnan(s) error("Lack of convergence (bad starting values?)") end
        if s<tol break end
        if method==1
            W = pinv(-H + gamma*(S'S)[1]*eye(p))
        elseif method==2
            W = pinv(-H + gamma*S*S')
        else
            W = pinv(-H)
        end
        if trace>0
            print(S')
        end
        delta = W*S
        x0 = x; loglik0 = loglik
        lambda=1;
        while s>=s0 || loglik<loglik0
            x = x0+lambda*delta
            val = f(x)
            counter += 1
            ##S = val["grad"]';  H = val["hess"]; loglik=val["val"]
            S = val.gradient::Vector{Float64};  H = val.hessian::Matrix{Float64}; loglik=val.value::Float64
            s = mean(S.*S)
            if isnan(s) s=s0 end
            lambda = lambda/2.0
            if (lambda<1e-12)
                error("Lack of convergence (bad starting values?)")
            end
        end
    end
    return x,val,counter
end

###}}} Newton-Raphson


text:
using LazySets, MathematicalPredicates, SparseArrays

function helicopter_specification()
    # initial set:
    # xᵢ ∈ [-0.1, 0.1] for i ∈ [1, 8]
    # xᵢ = 0 otherwise
    X0 = Hyperrectangle(spzeros(28), sparsevec(1:8, fill(0.1, 8), 28))

    # safety property: x1 ≤ 0.12
    property = is_contained_in(HalfSpace(sparsevec([1], [1.0], 28), 0.12))

    # time horizon: 20 time units
    time_horizon = 20.0

    # specification
    O = Dict(:T => time_horizon, :property => property)

    return X0, O
end

ext: jl

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
struct NewtonOptimizerCache{T, AT <: AbstractArray{T}}
    x̄::AT
    x::AT
    δ::AT
    g::AT

    function NewtonOptimizerCache(x::AT) where {T, AT <: AbstractArray{T}}
        new{T,AT}(zero(x), zero(x), zero(x), zero(x))
    end
end

function update!(cache::NewtonOptimizerCache, x::AbstractVector)
    cache.x̄ .= x
    cache.x .= x
    cache.δ .= 0
    return cache
end

function update!(cache::NewtonOptimizerCache, x::AbstractVector, g::AbstractVector)
    update!(cache, x)
    cache.g .= g
    return cache
end

function initialize!(cache::NewtonOptimizerCache, x::AbstractVector)
    cache.x̄ .= eltype(x)(NaN)
    cache.x .= x
    cache.δ .= eltype(x)(NaN)
    cache.g .= eltype(x)(NaN)
    return cache
end

"create univariate objective for linesearch algorithm"
function linesearch_objective(objective::MultivariateObjective, cache::NewtonOptimizerCache{T}) where {T}
    function f(α)
        cache.x .= cache.x̄ .+ α .* cache.δ
        value(objective, cache.x)
    end

    function d(α)
        cache.x .= cache.x̄ .+ α .* cache.δ
        dot(gradient!(objective, cache.x), cache.δ)
    end

    UnivariateObjective(f, d, one(T))
end

struct NewtonOptimizerState{OBJ <: MultivariateObjective, HES <: Hessian, LS <: LinesearchState, NOC <: NewtonOptimizerCache} <: OptimizationAlgorithm
    objective::OBJ
    hessian::HES
    linesearch::LS
    cache::NOC

    function NewtonOptimizerState(objective::OBJ, hessian::HES, linesearch::LS, cache::NOC) where {OBJ, HES, LS, NOC}
        new{OBJ,HES,LS,NOC}(objective, hessian, linesearch, cache)
    end
end

function NewtonOptimizerState(x::VT, objective::MultivariateObjective; hessian = HessianAD, linesearch = Bisection) where {XT, VT <: AbstractVector{XT}}
    cache = NewtonOptimizerCache(x)
    hess = hessian(objective, x)
    ls = LinesearchState(linesearch, linesearch_objective(objective, cache))

    NewtonOptimizerState(objective, hess, ls, cache)
end

NewtonOptimizer(args...; kwargs...) = NewtonOptimizerState(args...; kwargs...)
BFGSOptimizer(args...; kwargs...) = NewtonOptimizerState(args...; hessian = HessianBFGS, kwargs...)
DFPOptimizer(args...; kwargs...) = NewtonOptimizerState(args...; hessian = HessianDFP, kwargs...)

OptimizerState(algorithm::Newton, objective, x, y; kwargs...) = NewtonOptimizerState(x, objective; kwargs...)
OptimizerState(algorithm::BFGS, objective, x, y; kwargs...) = NewtonOptimizerState(x, objective; hessian = HessianBFGS, kwargs...)
OptimizerState(algorithm::DFP, objective, x, y; kwargs...) = NewtonOptimizerState(x, objective; hessian = HessianDFP, kwargs...)

objective(newton::NewtonOptimizerState) = newton.objective
hessian(newton::NewtonOptimizerState) = newton.hessian
cache(newton::NewtonOptimizerState) = newton.cache
linesearch(newton::NewtonOptimizerState) = newton.linesearch

function initialize!(newton::NewtonOptimizerState, x::AbstractVector)
    initialize!(hessian(newton), x)
    initialize!(cache(newton), x)
end

# function solver_step!(opt::Optimizer{<:Newton})
function solver_step!(x, newton::NewtonOptimizerState)
    # shortcuts
    # x̄ = cache(newton).x̄
    δ = cache(newton).δ
    g = cache(newton).g

    # update cache
    update!(cache(newton), x, gradient!(objective(newton), x))

    # update Hessian
    update!(hessian(newton), x)

    # solve H δx = - ∇f
    ldiv!(δ, hessian(newton), g)
    rmul!(δ, -1)

    # apply line search
    α = newton.linesearch()

    x .+= α .* δ
end


text:
using Documenter, MathML

makedocs(sitename="My Documentation")

deploydocs(
    repo = "github.com/anandijain/MathML.jl.git",
    devbranch = "main",
)


text:
# When working with rotation we have three distinct sets of indexes.
#
# 1. Original piece indexes: they come from the instance and are
#    rotation-unaware (the instance may be used to solve the fixed or
#    rotation-allowed variants).
# 2. Rotation-aware indexes: the indexes of the dummy pieces used to
#    simulate rotation. One original piece may: (i) create just one dummy
#    piece because it is a square or the rotated piece does not fit the
#    large object; (ii) have two dummies but the two dummies are also
#    associated and share demand with another original piece that is
#    the exact rotation of the first original piece; (iii) have two
#    dummies but they are not shared with any other original piece.
#    Consequently the dummy pieces are at least as many as the original
#    pieces and at most the double.
# 3. Shared demand indexes: the indexes of the demand vector in a
#    problem with rotation are always fewer or the same as the number of
#    original (or dummy) pieces. They are smaller if original pieces have
#    their rotation among other original pieces. They may be the same size
#    if all pieces are squares (for example). The smallest this set can go
#    is half the amount of original pieces (in the case every original
#    piece is a rotation of other).
#
#	In the structure below, `rai2sdi` maps (2) to (3), and `sdi2rai` does the
#	inverse (this is, maps (3) to (2)). `shared_demand` is the `d` of a
#	problem with rotation allowed. `sdi2opi` maps (3) to (1), therefore
#	allowing to map (2) to (1) by using `rai2sdi` and then `sdi2opi`.
#	There is no `opi2sdi` just because the code until now has no necessity of
#	this mapping (and it can be deduced from the other if really necessary).
struct RotationAwareData{D, S}
	"Maps rotation-aware indexes to the respective `shared_demand` index."
	rai2sdi :: Vector{D} # rotation-aware piece index to shared demand index
	"Maps `shared_demand` indexes to the respective rotation-aware indexes."
	sdi2rai :: Vector{Union{D, Tuple{D, D}}}
	"Maps `shared_demand` indexes to the respective original piece indexes."
	sdi2opi :: Vector{Union{D, Tuple{D, D}}}
	"The demand of pairs/single rotation-aware pieces."
	shared_demand :: Vector{D}
	# Not used until now, but saved to avoid unnecessary rework.
	"The lengths and widths of the rotation-aware pieces (dummies)."
	ra_sllw :: SortedLinkedLW{D, S}
	# In `get_cut_pattern`, a CutPattern using the dummy pieces is converted to a
	# CutPattern using the original pieces. In the case of perfect rotations
	# within the original pieces, we need to know how much each original piece
	# contributed to the shared demand otherwise the procedure cannot do the
	# conversion in a way that guarantee a valid demand for the original pieces.
	"The demand of the original pieces. Needed for `get_cut_pattern`."
	original_demand :: Vector{D}
end

function _find_perfect_rotation(
	pil :: S, piw :: S, pip :: P, sllw :: SortedLinkedLW{D, S}, p
) :: Int where {D, S, P}
	@assert pil != piw # This code assumes pil x piw is different from piw x pil
	# Search for pieces with an width equal to the piece length and length
	# equal to the piece width (in O(log n) as we are using sorted vectors).
	# (sorted) width/length index range
	swir = searchsorted(sllw.sw, pil)
	slir = searchsorted(sllw.sl, piw)
	# If no piece as pil as width or piw as length then there is no rotation.
	(isempty(swir) || isempty(slir)) && return 0
	# Get the found pieces to the the same numbering (i.e., the instance ordering)
	wir = sllw.swi2pii[swir]
	lir = sllw.sli2pii[slir]
	# Get only the piece indexes that have BOTH pil as width AND piw as length.
	r = intersect(wir, lir)
	# If this set is empty then there is not rotation.
	isempty(r) && return 0
	# Get the profit of the pieces that are rotations.
	pr = p[r]
	# Get the indexes of the right profits in the rotation range.
	cpri = findall(isequal(pip), pr) # correct profit range indexes
	# If no rotation has the same profit, then there are no "perfect rotation"s
	isempty(cpri) && return 0
	# If more than a piece is a perfect rotation then the instance is invalid:
	# there are two pieces with the same dimensions and profits (they should
	# be the same piece with a higher demand).
	if length(cpri) > 1
		idxs = r[cpri]
		@assert allsame(sllw.l[idxs])
		@assert allsame(sllw.w[idxs])
		@assert first(sllw.l[idxs]) == piw
		@assert first(sllw.w[idxs]) == pil
		error(
			"The rotation procedure detected two pieces in the instance that have the same dimensions ($(piw)x$(pil)) and the same profit ($pip). This is considered invalid, both pieces should be aggregated as a single piece with the summed demand. This was only detected because there is a third piece that is a rotation of the other two and this broke some assumptions of the rotation procedure."
		)
	end

	return r[only(cpri)]
end

# If we generalize this to multiple large object (i.e., large objects)
# with distinct dimensions we need to deal with the "rotated small object does
# not fit into the large object" in a better way or just ignore it (i.e.,
# create the dummy rotated object, even if it will not be used).
function build_RAD(
	L :: S, W :: S, sllw :: SortedLinkedLW{D, S}, d :: AbstractVector{D},
	p :: AbstractVector{P}
) :: RotationAwareData{D, S} where {D, S, P}
	@assert length(d) == length(p)
	@assert length(sllw.l) == length(d)

	ral, raw = S[], S[]
	shared_demand = D[]

	rai2sdi = D[]
	sdi2rai = Union{D, Tuple{D, D}}[]
	sdi2opi = Union{D, Tuple{D, D}}[]

	processed = falses(length(d))

	for opi in 1:length(d)
		processed[opi] && continue

		# First, let us add the unrotated dummy piece.
		curr_l = sllw.l[opi]
		curr_w = sllw.w[opi]
		push!(ral, curr_l)
		push!(raw, curr_w)
		# We start assuming the demand of the original piece.
		push!(shared_demand, d[opi])
		push!(rai2sdi, lastindex(shared_demand)) # this cannot change later

		# If the piece is a square or the rotation do not fit, then there
		# is no rotated dummy piece.
		if curr_l == curr_w || curr_l > W || curr_w > L
			push!(sdi2rai, lastindex(rai2sdi))
			push!(sdi2opi, opi)
			processed[opi] = true
			continue
		end

		pr_opi = _find_perfect_rotation(last(ral), last(raw), p[opi], sllw, p)
		if !iszero(pr_opi) # a perfect rotation was found
			# Create the rotated dummy piece. Share its demand.
			push!(ral, curr_w)
			push!(raw, curr_l)
			shared_demand[end] += d[pr_opi]
			push!(rai2sdi, lastindex(shared_demand)) # this one is for the rotation

			# map the demand index to the two dummies, an to the two original
			curr_rai = lastindex(rai2sdi)
			push!(sdi2rai, (curr_rai - 1, curr_rai))
			push!(sdi2opi, (opi, pr_opi))

			# mark both original counterparts as processed
			processed[opi] = true
			processed[pr_opi] = true
			continue
		end

		# If the piece has nothing that hinders rotation, nor another original
		# that is a perfect rotation, then it just spawns a second rotated dummy.
		push!(ral, curr_w)
		push!(raw, curr_l)
		# there is no change to the shared demand here
		push!(rai2sdi, lastindex(shared_demand)) # this one is for the rotation

		# map the demand index to the two dummies, and to the single original
		curr_rai = lastindex(rai2sdi)
		push!(sdi2rai, (curr_rai - 1, curr_rai))
		push!(sdi2opi, opi)

		processed[opi] = true
	end

	return RotationAwareData{D, S}(
		rai2sdi, sdi2rai, sdi2opi, shared_demand, SortedLinkedLW(D, ral, raw),
		deepcopy(d)
	)
end


text:
@inline area(face::AbstractFace2Node1D,nodes::AbstractArray{<:Vec1D}) = 1
@inline area(face::AbstractFace2Node2D,nodes::AbstractArray{<:Vec2D}) =  normal_to_2Dline(face.ind,nodes)

@inline function center(face::AbstractFace2Node1D,nodes::AbstractArray{<:Vec1D}); @inbounds r = nodes[face[1]]; return r; end
@inline center(face::AbstractFace2Node2D,nodes::AbstractArray{<:Vec2D}) =  line_center(face.ind,nodes)

struct Face2Node{N} <: AbstractFace2Node{N}
    ind::NTuple{N,UInt}
end

function Base.show(io::IO,f::Face2Node{N}) where N
    compact = get(io, :compact, false)

    if !compact
        print(io,"$(N)DFace")
    end
    print(io,"(",join(f.ind,","),")")
end

abstract type AbstractFace2Cell <: NTupleWrap{2,UInt} end

abstract type AbstractBFace2Cell <: NTupleWrap{1,UInt} end

struct Face2Cell <: AbstractFace2Cell
    ind::Tuple{UInt,UInt}
end

function Base.show(io::IO,f::Face2Cell)
    compact = get(io, :compact, false)

    if !compact
        print(io,"Interface of cells ",f[1]," and ",f[2])
    else
        print(io,"(",join(f.ind,","),")")
    end
end

struct BFace2Cell <: AbstractBFace2Cell
    ind::Tuple{UInt}
end

function Base.show(io::IO,f::BFace2Cell)
    compact = get(io, :compact, false)

    if !compact
        print(io,"Boundary face of cell ",f[1])
    else
        print(io,f[1])
    end
end

struct Face2CellLoop{NBF,NF,VecType,ArrayVecType}
    bf2c::Vector{BFace2Cell}
    bft::Vector{UInt} # List of of boundary faces types
    bfn::ArrayVecType # List of normals of boundary faces
    bfc::ArrayVecType # List of centers of boundary faces
    bcv::Vector{Float64} # List of 1/volumes of boundary faces owner cells
    bccenter::ArrayVecType # List of centers of boundary faces owner cells
    f2c::Vector{Face2Cell}
    fn::ArrayVecType # List of normals of faces
    fc::ArrayVecType # List of centers of faces
    cv::Vector{NTuple{2,Float64}} # List of 1/volumes of faces owner cells
    ccenter::Vector{NTuple{2,VecType}} # List of centers of faces owner cells
    AoL::Vector{Float64}
end

nbfaces(::Type{Face2CellLoop{NBF,NF,Vectype,AV}}) where {NBF,NF,Vectype,AV} = NBF
@inline nbfaces(a::A) where A<:Face2CellLoop = nbfaces(A)
nfaces(::Type{Face2CellLoop{NBF,NF,Vectype,AV}}) where {NBF,NF,Vectype,AV} = NF
@inline nfaces(a::A) where A<:Face2CellLoop = nfaces(A)

struct FaceSimpleInterpolation{T,VecT,NF} <: AbstractVector{T}
    Tc::VecT
    f2c::Vector{Face2Cell}
end

FaceSimpleInterpolation(Tc,mesh) = FaceSimpleInterpolation{eltype(Tc),typeof(Tc),length(mesh.f2cloops.f2c)}(Tc,mesh.f2cloops.f2c)

Base.size(a::Type{FaceSimpleInterpolation{T,VT,NF}}) where {T,VT,NF} = (NF,)
@inline Base.size(a::A) where {A<:FaceSimpleInterpolation} = size(A)
Base.length(a::Type{FaceSimpleInterpolation{T,VT,NF}}) where {T,VT,NF} = NF
@inline Base.length(a::A) where {A<:FaceSimpleInterpolation} = length(A)
Base.linearindices(a::FaceSimpleInterpolation{T,VT,NF}) where {T,VT,NF} = Base.OneTo(length(a))
Base.IndexStyle(::Type{FaceSimpleInterpolation}) = IndexLinear()
@inline Base.@propagate_inbounds function Base.getindex(a::FaceSimpleInterpolation,I)
    j1,j2 = a.f2c[I]
    Tc = a.Tc
    return 0.5*(Tc[j1]+Tc[j2])
end

struct FaceVector{T,NF,NBF,L} <: AbstractVector{T}
    data::Vector{T}
end

Base.length(::Type{<:FaceVector{T,NF,NBF,L}}) where {T,NF,NBF,L} = L
Base.size(::Type{<:FaceVector{T,NF,NBF,L}}) where {T,NF,NBF,L} = (L,)
Base.linearindices(::Type{<:FaceVector{T,NF,NBF,L}}) where {T,NF,NBF,L} = Base.OneTo(L)
Base.IndexStyle(::Type{<:FaceVector{T,NF,NBF,L}}) where {T,NF,NBF,L} = IndexLinear()
Base.similar(a::FaceVector{T,NF,NBF,L}) where {T,NF,NBF,L} = FaceVector{T,NF,NBF,L}(zeros(T,L))

for op in (:(Base.length),:(Base.size),:(Base.linearindices))
    @eval begin
        @inline function ($op)(a::FaceVector)
            return ($op)(typeof(a))
        end
    end
end

@inline function Base.getindex(a::FaceVector,i::Integer)
    d = a.data
    @boundscheck checkbounds(a,i)
    @inbounds r = d[i]
    return r
end

@inline function Base.setindex!(a::FaceVector,x,i::Integer)
    d = a.data
    @boundscheck checkbounds(a,i)
    @inbounds d[i] = x
    return d
end

@inline @generated function Base.getindex(a::FaceVector{T,NF,NBF,L},s::Symbol,i1::Integer) where {T,NF,NBF,L}
    n = NF
    quote
        d = a.data
        i = i1 + $n
        @boundscheck checkbounds(a,i)
        @inbounds r = d[i]
        return r
    end
end

@inline @generated function Base.setindex!(a::FaceVector{T,NF,NBF,L},x,s::Symbol,i1::Integer) where {T,NF,NBF,L}
    n = NF
    quote
        d = a.data
        i = i1 + $n
        @boundscheck checkbounds(a,i)
        @inbounds d[i] = x
        return d
    end
end


text:
module optim

using LinearAlgebra, Statistics
export SGD, ADAM, Newton, optim_step, done

mutable struct SGD
    α
    x
    max_iter
    grad_tol
    log_file
    iter
    done

    function SGD(α :: Float64, x0 :: Array{Float64}, max_iter :: Int, grad_tol :: Float64, log_name :: String = nothing)
        if log_name !== nothing
            log_file = open(log_name, "w+")
        else
            log_file = nothing
        end

        new(α, x0, max_iter, grad_tol, log_file, 0, false)
    end
end

function log(o :: SGD, s :: String)
    if o.log_file !== nothing
        write(o.log_file, s)
        flush(o.log_file)
    end
end

function optim_step(o :: SGD, ss :: Matrix{Float64}, target_Es :: Vector{Float64})
    g = (mean(ss, dims=2) - target_Es)[:, 1]
    optim_step(o, g)
end

function optim_step(o :: SGD, g :: Vector{Float64})
    gn = norm(g)

    if o.iter > o.max_iter
        log(o, "Maximum iterations ($(o.max_iter)) reached. Terminating.\n")
        finish(o)
    elseif gn ≤ o.grad_tol
        log(o, "Gradient tolerance (|g| = $gn ≤ $(o.grad_tol)) reached. Terminating.\n")
        finish(o)
    else
        log(o, "Iteration $(o.iter): |g| = $gn.\n")
        o.x -= o.α * g
        o.iter += 1
    end

    o.x
end

function finish(o :: SGD)
    o.done = true
    close(o.log_file)
end

function done(o :: SGD)
    o.done
end

mutable struct ADAM
    α
    β1
    β2
    ε
    x
    m
    v
    max_iter
    grad_tol
    log_file
    iter
    done

    function ADAM(
        α :: Float64, β1 :: Float64, β2 :: Float64, ε :: Float64,
        x0 :: Array{Float64}, max_iter :: Int, grad_tol :: Float64,
        log_name :: String = nothing
    )
        if log_name !== nothing
            log_file = open(log_name, "w+")
        else
            log_file = nothing
        end

        m0 = zeros(size(x0))
        v0 = zeros(size(x0))

        new(
            α, β1, β2, ε,
            x0, m0, v0,
            max_iter, grad_tol, log_file,
            0, false
        )
    end
end

function log(o :: ADAM, s :: String)
    if o.log_file !== nothing
        write(o.log_file, s)
        flush(o.log_file)
    end
end

function optim_step(o :: ADAM, ss :: Matrix{Float64}, target_Es :: Vector{Float64})
    G = ss .- target_Es
    g = mean(G, dims=2)
    gn = norm(g)

    if o.iter > o.max_iter
        log(o, "Maximum iterations ($(o.max_iter)) reached. Terminating.\n")
        finish(o)
    elseif gn ≤ o.grad_tol
        log(o, "Gradient tolerance (|g| = $gn ≤ $(o.grad_tol)) reached. Terminating.\n")
        finish(o)
    else
        log(o, "Iteration $(o.iter): |g| = $gn.\n")
        o.iter += 1
        o.m = o.β1 * o.m + (1 - o.β1) * g
        o.v = o.β2 * o.v + (1 - o.β2) * mean(G .^ 2, dims=2)
        mc = o.m / (1 - o.β1 ^ o.iter)
        vc = o.v / (1 - o.β2 ^ o.iter)
        o.x -= o.α * mc ./ (sqrt.(vc) .+ o.ε)
    end

    o.x
end

function finish(o :: ADAM)
    o.done = true
    close(o.log_file)
end

function done(o :: ADAM)
    o.done
end

mutable struct Newton
    x
    max_iter
    grad_tol
    log_file
    iter
    done

    function Newton(x0 :: Vector{Float64}, max_iter :: Int, grad_tol :: Float64, log_name :: String = nothing)
        if log_name !== nothing
            log_file = open(log_name, "w+")
        else
            log_file = nothing
        end

        new(x0, max_iter, grad_tol, log_file, 0, false)
    end
end

function log(o :: Newton, s :: String)
    if o.log_file !== nothing
        write(o.log_file, s)
        flush(o.log_file)
    end
end

function optim_step(o :: Newton, ss :: Matrix{Float64}, target_Es :: Vector{Float64})
    p, n = size(ss)
    g = mean(ss, dims=2) - target_Es
    gn = norm(g)

    if o.iter > o.max_iter
        log(o, "Maximum iterations ($(o.max_iter)) reached. Terminating.\n")
        finish(o)
    elseif gn ≤ o.grad_tol
        log(o, "Gradient tolerance (|g| = $gn ≤ $(o.grad_tol)) reached. Terminating.\n")
        finish(o)
    else
        log(o, "Iteration $(o.iter): |g| = $gn.\n")
        H = 1 / (n - 1) * ss * ss'
        o.x -= H^-1 * g
        o.iter += 1
    end

    o.x
end

function finish(o :: Newton)
    o.done = true
    close(o.log_file)
end

function done(o :: Newton)
    o.done
end

end


text:
#!/usr/bin/env julia=".."
using RobotOS
@rosimport std_msgs.msg.String
rostypegen()
using .std_msgs.msg: StringMsg
function callback(msg)
    loginfo("$(RobotOS.get_caller_id()) I heard $(msg.data)")
end
function listener()
    init_node("listener")  # node name
    # subscribe to the topic "chatter"; the keyword argument queue_size is optional
    sub = Subscriber{StringMsg}("chatter", callback; queue_size=10)
    # spin() simply keeps Julia from exiting until this node is stopped
    spin()
end
if !isinteractive()
    listener()
end

max_stars_repo_name: mattborghi/ros-julia-examples

max_issues_repo_name: mattborghi/ros-julia-examples

max_forks_repo_name: mattborghi/ros-julia-examples


text:
"""
    x = solve_pcg(A,b; maxiter=5000,eps=1e-15,reset=50,verbose=false)

Solve the symmetric positive definite system of equations:
`A` as sparse left-hand side and `b` as a vector of right-hand side.
The Jacobi preconditioner (`diag(A)`) is used.
The maximum number of iterations is specified with `maxiter`.
Convergence criterion `eps` is `norm(b-A*x)/norm(b)`.
The residual vector will be recalculated every `reset` rounds.
The details will be shown with `verbose=true`.

```juliadoctest
julia> x = solve_pcg(A,b)
```
"""
function solve_pcg(A, b::Vector{Tv};
                   maxiter=5000,eps=1e-15,reset=50,verbose=false) where {Tv}
   n = size(A,1)
   x = zeros(Tv,n)
   r = zeros(Tv,n)
   p = zeros(Tv,n)
   z = zeros(Tv,n)
   w = zeros(Tv,n)
   Mi = Diagonal(1 ./ diag(A))
   r .= b
   oldtau = 1.0
   for i=1:maxiter
      z .= Mi*r
      tau = z'*r
      if i==1
         beta = 0
         p .= z
      else
         beta = tau/oldtau
         p .= z + beta*p
      end
      w .= A*p
      alpha = tau/(p'*w)
      x .= x + alpha*p
      if i % reset == 0
         r .= b - A*x
      else
         r .= r - alpha*w
      end
      conv = norm(r)^2/norm(b)^2
      if verbose
         print(@sprintf("round %5d: conv %8.3e\n",i,conv))
      end
      if conv<eps; break; end
      oldtau = tau
   end
   x
end

"""
   tr = approximate_trace_of_inverse(A)
"""
function approximate_trace_of_inverse(A; rounds=1000, msglev=0)
   if rounds<1
      throw( ArgumentError("negative rounds") )
   end
   tr = Threads.Atomic{Float64}(0.0)
   done = Threads.Atomic{Int}(0)
   n = size(A,1)
   @threads for iter=1:rounds
      atomic_add!(done, 1)
      z = [ifelse(rand()>0.5,1.0,-1.0) for i in 1:n]
      #tr = tr + z'*y
      if msglev>0
         @time y = solve_pcg(A,z)
         println("iter $(done.value) thread $(Threads.threadid())")
      else
         y = solve_pcg(A,z)
      end
      atomic_add!(tr, z'*y)
   end
   return tr.value/rounds
end


text:
using DiffUtils
using Documenter

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

makedocs(;
    modules=[DiffUtils],
    authors="Simeon Schaub <simeondavidschaub99@gmail.com> and contributors",
    repo="https://github.com/simeonschaub/DiffUtils.jl/blob/{commit}{path}#L{line}",
    sitename="DiffUtils.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://simeonschaub.github.io/DiffUtils.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/simeonschaub/DiffUtils.jl",
    devbranch="main",
)


text:
# Erode and dilate support 3x3 regions only (and higher-dimensional generalizations).

"""
```
imgd = dilate(img, [region])
```

perform a max-filter over nearest-neighbors. The
default is 8-connectivity in 2d, 27-connectivity in 3d, etc. You can specify the
list of dimensions that you want to include in the connectivity, e.g., `region =
[1,2]` would exclude the third dimension from filtering.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(5, 5); img[3, 3] = 1.; img
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> dilate(img)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> dilate(img, 1)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
dilate(img::AbstractArray, region=coords_spatial(img)) = dilate!(copy(img), region)

"""
```
imge = erode(img, [region])
```

perform a min-filter over nearest-neighbors. The
default is 8-connectivity in 2d, 27-connectivity in 3d, etc. You can specify the
list of dimensions that you want to include in the connectivity, e.g., `region =
[1,2]` would exclude the third dimension from filtering.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(5, 5); img[2:4, 2:4] .= 1.; img
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> erode(img)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> erode(img, 1)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
erode(img::AbstractArray, region=coords_spatial(img)) = erode!(copy(img), region)

dilate!(maxfilt, region=coords_spatial(maxfilt)) = extremefilt!(maxfilt, max, region)
erode!(minfilt, region=coords_spatial(minfilt)) = extremefilt!(minfilt, min, region)
function extremefilt!(A::AbstractArray, select::Function, region=coords_spatial(A))
    inds = axes(A)
    for d = 1:ndims(A)
        if size(A, d) == 1 || !in(d, region)
            continue
        end
        Rpre = CartesianIndices(inds[1:d-1])
        Rpost = CartesianIndices(inds[d+1:end])
        _extremefilt!(A, select, Rpre, inds[d], Rpost)
    end
    A
end

@noinline function _extremefilt!(A, select, Rpre, inds, Rpost)
    # TODO: improve the cache efficiency
    for Ipost in Rpost, Ipre in Rpre
        # first element along dim
        i1 = first(inds)
        a2, a3 = A[Ipre, i1, Ipost], A[Ipre, i1+1, Ipost]
        A[Ipre, i1, Ipost] = select(a2, a3)
        # interior along dim
        for i = i1+2:last(inds)
            a1, a2 = a2, a3
            a3 = A[Ipre, i, Ipost]
            A[Ipre, i-1, Ipost] = select(select(a1, a2), a3)
        end
        # last element along dim
        A[Ipre, last(inds), Ipost] = select(a2, a3)
    end
    A
end

"""
`imgo = opening(img, [region])` performs the `opening` morphology operation, equivalent to `dilate(erode(img))`.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(5, 5); img[1, 1] = 1.; img[3:5, 3:5] .= 1.; img
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0

julia> opening(img)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
```
"""
opening(img::AbstractArray, region=coords_spatial(img)) = opening!(copy(img), region)
opening!(img::AbstractArray, region=coords_spatial(img)) = dilate!(erode!(img, region),region)

"""
`imgc = closing(img, [region])` performs the `closing` morphology operation, equivalent to `erode(dilate(img))`.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(7, 7); img[3:5, 3:5] .= 1.; img[4, 4] = 0.; img
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> closing(img)
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
closing(img::AbstractArray, region=coords_spatial(img)) = closing!(copy(img), region)
closing!(img::AbstractArray, region=coords_spatial(img)) = erode!(dilate!(img, region),region)

"""
`imgth = tophat(img, [region])` performs `top hat` of an image,
which is defined as the image minus its morphological opening.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(5, 5); img[1, 1] = 1.; img[3:5, 3:5] .= 1.; img
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0

julia> tophat(img)
5×5 Array{Float64,2}:
 1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
```
"""
tophat(img::AbstractArray, region=coords_spatial(img)) = img - opening(img, region)

"""
`imgbh = bothat(img, [region])` performs `bottom hat` of an image,
which is defined as its morphological closing minus the original image.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(7, 7); img[3:5, 3:5] .= 1.; img[4, 4] = 0.; img
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> bothat(img)
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
bothat(img::AbstractArray, region=coords_spatial(img)) = closing(img, region) - img

"""
`imgmg = morphogradient(img, [region])` returns morphological gradient of the image,
which is the difference between the dilation and the erosion of a given image.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(7, 7); img[3:5, 3:5] .= 1.; img
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> morphogradient(img)
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  0.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
morphogradient(img::AbstractArray, region=coords_spatial(img)) = dilate(img, region) - erode(img, region)

"""
`imgml = morpholaplace(img, [region])` performs `Morphological Laplacian` of an image,
which is defined as the arithmetic difference between the internal and the external gradient.
`region` allows you to control the dimensions over which this operation is performed.

# Examples

```jldoctest; setup = :(using ImageMorphology), filter = r"Array{Float64,2}|Matrix{Float64}"
julia> img = zeros(7, 7); img[3:5, 3:5] .= 1.; img[4, 4] = 0.; img
7×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0  1.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0

julia> morpholaplace(img)
7×7 Array{Float64,2}:
 0.0  0.0   0.0   0.0   0.0  0.0  0.0
 0.0  1.0   1.0   1.0   1.0  1.0  0.0
 0.0  1.0  -1.0  -1.0  -1.0  1.0  0.0
 0.0  1.0  -1.0   1.0  -1.0  1.0  0.0
 0.0  1.0  -1.0  -1.0  -1.0  1.0  0.0
 0.0  1.0   1.0   1.0   1.0  1.0  0.0
 0.0  0.0   0.0   0.0   0.0  0.0  0.0
```
"""
morpholaplace(img::AbstractArray, region=coords_spatial(img)) = dilate(img, region) + erode(img, region) - 2img


text:
module Target

function baz(y::Float64)
    return y + 40.0
end

function foo(x::Float64)
    return baz(x + 20.0)
end

end

struct StateMix <: CompilationContext end

allow(ctx::StateMix, m::Module, fn, args...) = m == Target

show_after_inference(ctx::StateMix) = false
show_after_optimization(ctx::StateMix) = false
debug(ctx::StateMix) = false

mutable struct Recorder
    d::Dict
    ret::Any
    Recorder() = new(Dict(), nothing)
    Recorder(d, ret) = new(d, ret)
end

swap(r, e) = e
function swap(r, e::Expr)
    e.head == :call || return e
    return Expr(:call, r, e.args[1:end]...)
end

function transform(::StateMix, b)
    pushfirst!(b, Expr(:call, Recorder))
    for (v, st) in b
        e = swap(Core.SSAValue(1), st)
        v == 1 || replace!(b, v, e)
    end
    return b
end

@testset "Insert state" begin
    Mixtape.@load_call_interface()

    function (r::Recorder)(f::Function, args...)
        args = map(a -> a isa Recorder ? a.ret : a, args)
        rec = call(StateMix(), f, args...)
        if rec isa Recorder
            merge!(r.d, rec.d)
            r.d[(f, args...)] = rec.ret
            r.ret = rec.ret
        else
            r.d[(f, args...)] = rec
            r.ret = rec
        end
        return r
    end

    rec = call(StateMix(), Target.foo, 5.0)
    @test rec.d[(Target.baz, 25.0)] == 65.0
    @test rec.d[(Base.:(+), 5.0, 20.0)] == 25.0
    @test rec.d[(Base.:(+), 25, 40.0)] == 65.0
    @test rec.ret == 65.0
end


text:
import Pkg
Pkg.add("Documenter")
using Documenter
using FundamentosDMC
push!(LOAD_PATH,"../src/")
makedocs(
    modules=[FundamentosDMC],
    sitename="FundamentosDMC.jl",
    pages = [
        "Início" => "index.md",
        "Sistema simulado" => "sistema.md",
        "Simulación microcanónica" => "simple.md",
        "Termostato isocinético" => "isokinetic.md",
        "Termostato de Berendsen" => "berendsen.md",
        "Termostato de Langevin" => "langevin.md",
        "Monte-Carlo" => "montecarlo.md",
        "Comparando MD y MC" => "md_vs_mc.md",
        "Ejecución en 3D" => "running3D.md",
        "Distribución de velocidades" => "velocities.md",
        "Help entries" => "help.md",
    ]
)
deploydocs(
    repo = "github.com/m3g/FundamentosDMC.jl.git",
    target = "build",
    branch = "gh-pages",
    versions = ["stable" => "v^", "v#.#" ],
)


text:
#--------------------------------------
## LOAD PACKAGES
using RegNeuralDE, OrdinaryDiffEq, Flux, DiffEqFlux, Tracker
using YAML, Dates, BSON, Random, Statistics, Printf
using CUDA
using RegNeuralDE: accuracy
using Flux: @functor
using Tracker: TrackedReal
import Base.show

CUDA.allowscalar(false)
#--------------------------------------

#--------------------------------------
## CONFIGURATION
config_file = joinpath(pwd(), "experiments", "configs", "mnist_node.yml")
config = YAML.load_file(config_file)

Random.seed!(config["seed"])

hparams = config["hyperparameters"]
BATCH_SIZE = hparams["batch_size"]
EPOCHS = hparams["epochs"]
REGULARIZE = hparams["regularize"]
REG_TYPE = hparams["type"]
STEER = hparams["steer"]
identifier =
    REGULARIZE ? "$(string(now()))_$(REGULARIZE)_$(REG_TYPE)" : "$(string(now()))_vanilla"
EXPERIMENT_LOGDIR = joinpath(pwd(), "results", "mnist_node", identifier)
MODEL_WEIGHTS = joinpath(EXPERIMENT_LOGDIR, "weights.bson")
FILENAME = joinpath(EXPERIMENT_LOGDIR, "results.yml")

# Create a directory to store the results
isdir(EXPERIMENT_LOGDIR) || mkpath(EXPERIMENT_LOGDIR)
cp(config_file, joinpath(EXPERIMENT_LOGDIR, "config.yml"))
#--------------------------------------

#--------------------------------------
## NEURAL NETWORK
# The dynamics of the Neural ODE are time dependent
struct MLPDynamics{D1,D2}
    dense_1::D1
    dense_2::D2
end

@functor MLPDynamics

MLPDynamics(in::Integer, hidden::Integer) =
    MLPDynamics(Dense(in + 1, hidden, CUDA.tanh), Dense(hidden + 1, in, CUDA.tanh))

function (mlp::MLPDynamics)(x::AbstractMatrix, t::TrackedReal)
    _t = CUDA.ones(Float32, 1, size(x, 2)) .* t
    return mlp.dense_2(vcat(mlp.dense_1(vcat(x, _t)), _t))
end
#--------------------------------------

#--------------------------------------
## SETUP THE MODELS + DATASET + TRAINING UTILS
# Get the dataset
train_dataloader, test_dataloader = load_mnist(BATCH_SIZE, x -> cpu(x))

if REG_TYPE == "error_est"
    # Anneal the regularization so that it doesn't overpower the
    # the main objective
    λ₀ = 1.0f2
    λ₁ = 1.0f1
    save_func(u, t, integrator) = integrator.EEst * integrator.dt
    solver = Tsit5()
    global agg = mean
elseif REG_TYPE == "stiff_est"
    # No annealing is generally needed for stiff_est
    λ₀ = 0.1f0
    λ₁ = 0.1f0
    const stability_size =
        Tracker.TrackedReal(1 / Float32(OrdinaryDiffEq.alg_stability_size(Tsit5())))
    function save_func(u, t, integrator)
        stiff_est = abs(integrator.eigen_est)
        return stability_size * ((iszero(stiff_est) || isnan(stiff_est)) ? 0 : stiff_est)
    end
    global agg = maximum
    solver = AutoTsit5(Tsit5())
elseif REG_TYPE == "error_stiff_est"
    λ₀ = 1.0f1
    λ₁ = 1.0f1
    const mul_val = Tracker.TrackedReal(1.0f0)
    const stability_size =
        Tracker.TrackedReal(1 / Float32(OrdinaryDiffEq.alg_stability_size(Tsit5())))
    function save_func(u, t, integrator)
        err_est = integrator.EEst * integrator.dt
        eest = Tracker.data(err_est)
        stiff_est = integrator.eigen_est
        sest = Tracker.data(stiff_est)
        return (
            ((iszero(eest) || isnan(eest)) ? 0 : err_est) +
            0.1f0 * stability_size * ((iszero(sest) || isnan(sest)) ? 0 : stiff_est)
        ) * mul_val
    end
    global agg = mean
    solver = AutoTsit5(Tsit5())
else
    global agg = mean
    solver = Tsit5()
end
b = 0.5f0
sample_tspan_ubound() = [0.0f0, 1.0f0 - (2 * rand(Float32) - 1.0f0) * b]
k = log(λ₀ / λ₁) / EPOCHS
# Exponential Decay
λ_func(t) = λ₀ * exp(-k * t)

# Define the models
mlp_dynamics = MLPDynamics(784, 100)

node = ClassifierNODE(
    Chain(x -> reshape(x, 784, :)) |> track |> gpu,
    TrackedNeuralODE(
        mlp_dynamics |> track |> gpu,
        [0.0f0, 1.0f0],
        true,
        REGULARIZE,
        solver,
        save_everystep = false,
        reltol = 1.4f-8,
        abstol = 1.4f-8,
        save_start = false,
    ),
    Dense(784, 10) |> track |> gpu,
)
ps = Flux.trainable(node)

opt = Flux.Optimise.Optimiser(InvDecay(1.0e-5), Momentum(0.1, 0.9))

function loss_function(x, y, model, p1, p2, p3; λ = 1.0f2, notrack = false)
    tspan = STEER ? sample_tspan_ubound() : nothing
    pred, _, sv = model(x, p1, p2, p3; func = save_func, tspan = tspan)
    cross_entropy = Flux.Losses.logitcrossentropy(pred, y)
    reg = REGULARIZE ? λ * agg(sv.saveval) : zero(eltype(pred))
    total_loss = cross_entropy + reg
    if !notrack
        ce_un = cross_entropy |> untrack
        reg_un = reg |> untrack
        total_loss_un = total_loss |> untrack
        logger(
            false,
            Dict(
                "Total Loss" => total_loss_un,
                "Cross Entropy Loss" => ce_un,
                "Regularization" => reg_un,
            ),
        )
    end
    return total_loss
end
#--------------------------------------

#--------------------------------------
## LOGGING UTILITIES
nfe_counts = Vector{Float64}(undef, EPOCHS + 1)
train_accuracies = Vector{Float64}(undef, EPOCHS + 1)
test_accuracies = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes = Vector{Float64}(undef, EPOCHS + 1)
inference_runtimes = Vector{Float64}(undef, EPOCHS + 1)
train_runtimes[1] = 0

logger = table_logger(
    [
        "Epoch Number",
        "NFE Count",
        "Train Accuracy",
        "Test Accuracy",
        "Train Runtime",
        "Inference Runtime",
    ],
    ["Total Loss", "Cross Entropy Loss", "Regularization"],
)
#--------------------------------------

#--------------------------------------
## RECORD DETAILS BEFORE TRAINING STARTS
dummy_data = train_dataloader.data[1][:, :, :, 1:BATCH_SIZE] |> gpu |> track
stime = time()
_, _nfe, _ = node(dummy_data; func = save_func)
inference_runtimes[1] = time() - stime
train_runtimes[1] = 0.0
nfe_counts[1] = _nfe
train_accuracies[1] = accuracy(node, train_dataloader)
test_accuracies[1] = accuracy(node, test_dataloader)

logger(
    false,
    Dict(),
    0.0,
    nfe_counts[1],
    train_accuracies[1],
    test_accuracies[1],
    train_runtimes[1],
    inference_runtimes[1],
)
#--------------------------------------

#--------------------------------------
## WARMSTART THE GRADIENT COMPUTATION
y_ = zeros(Float32, 10, 1)
y_[1, 1] = 1.0
_ = Tracker.gradient(
    (p1, p2, p3) -> loss_function(
        rand(Float32, 28, 28, 1, 1) |> gpu |> track,
        y_ |> gpu |> track,
        node,
        p1,
        p2,
        p3;
        notrack = true,
    ),
    ps...,
)
#--------------------------------------

#--------------------------------------
## TRAINING
for epoch = 1:EPOCHS
    λ = λ_func(epoch - 1)
    timing = 0

    for (i, (x_, y_)) in enumerate(train_dataloader)
        x = x_ |> gpu |> track
        y = Float32.(y_) |> gpu |> track

        start_time = time()
        gs = Tracker.gradient(
            (p1, p2, p3) -> loss_function(x, y, node, p1, p2, p3; λ = λ),
            ps...,
        )
        update_parameters!(ps, gs, opt)
        timing += time() - start_time

        x = y = nothing
        GC.gc(true)
    end

    # Record the time per epoch
    train_runtimes[epoch+1] = timing

    # Record the NFE count
    start_time = time()
    _, nfe, _ = node(dummy_data)
    inference_runtimes[epoch+1] = time() - start_time
    nfe_counts[epoch+1] = nfe

    # Test and Train Accuracy
    train_accuracies[epoch+1] = accuracy(node, train_dataloader)
    test_accuracies[epoch+1] = accuracy(node, test_dataloader)

    logger(
        false,
        Dict(),
        epoch,
        nfe_counts[epoch+1],
        train_accuracies[epoch+1],
        test_accuracies[epoch+1],
        train_runtimes[epoch+1],
        inference_runtimes[epoch+1],
    )
end
logger(true, Dict())
#--------------------------------------

#--------------------------------------
## STORE THE RESULTS
results = Dict(
    :nfe_counts => nfe_counts,
    :train_accuracies => train_accuracies,
    :test_accuracies => test_accuracies,
    :train_runtimes => train_runtimes,
    :inference_runtimes => inference_runtimes,
)

weights = Flux.params(node) .|> cpu .|> untrack
BSON.@save MODEL_WEIGHTS weights

YAML.write_file(FILENAME, results)
#--------------------------------------


text:
"""
    Line

This is a line in the plane, constructed by two complex numbers or four real numbers:
* `Line(w,z)`
* `Line(x1,y1,x2,y2)`
"""
struct Line <: Cline
    a::Complex{Float64}
    b::Complex{Float64}

    function Line(aa::Number, bb::Number)
        if aa == bb || isinf(aa) || isinf(bb)
            error("Line must be specified by two distinct finite points")
        end
        new(aa, bb)
    end
end

Line(x::Real, y::Real, xx::Real, yy::Real) = Line(Complex(x, y), Complex(xx, yy))

export Line

import Base: (==)

function (==)(L::Line, M::Line)
    return collinear(L.a, L.b, M.a) && collinear(L.a, L.b, M.b)
end

import Base: in
"""
    in(z,L::Line)
Check if the point `z` lies on the line `L`.
"""
function in(z::Number, L::Line)
    return collinear(Complex(z), L.a, L.b)
end

export slope
function slope(L::Line)::Float64
    x, y = reim(L.a)
    xx, yy = reim(L.b)
    return (y - yy) / (x - xx)
end

"""
    sort_points(a,b,c)

Given three collinear points (expressed as complex numbers) return them as a
3-tuple ensuring the 2nd entry in the tuple lies between the other two points
on the line that contains them.
"""
function sort_points(a::Number, b::Number, c::Number)
    if !collinear(a, b, c)
        @warn "The three points are not collinear. Nonsense may ensue."
    end

    rat = real((c - a) / (b - a))

    if rat >= 1
        return a, b, c
    end

    if rat <= 0
        return b, a, c
    end

    return a, c, b
end

"""
    three_points(L::Line)
Return a list of three points on the line `L`.
"""
function three_points(L::Line)
    x = L.a
    y = L.b
    z = (x + y) / 2
    return [x, z, y]
end

"""
    dilate(L::Line, factor=2)
Create a new line in which the distance between the defining points, `X.a` and `X.b`,
is magnified by `factor` from those of `L`. The resulting line is equal to `L`.
"""
function dilate(L::Line, factor = 2)
    a = L.a
    b = L.b
    m = (a + b) / 2

    aa = m + factor * (a - m)
    bb = m + factor * (b - m)

    return Line(aa, bb)
end
export dilate

# THIS IS A BAD DRAWING METHOD FOR LINES, BUT WILL HAVE TO DO FOR NOW
# It simply draws an arrowed line segment between L.a and L.b
# See: dilate

function draw(L::Line; args...)
    a = L.a
    c = L.b
    b = (a + c) / 2
    draw_segment(b, a, arrow = true, linecolor = :black, args...)
    draw_segment(b, c, arrow = true, linecolor = :black, args...)
end


text:
module BaseTraits
using SimpleTraits

export IsLeafType, IsConcrete, IsBits, IsImmutable, IsContiguous, IsIndexLinear,
       IsAnything, IsNothing, IsCallable, IsIterator

"Trait which contains all types"
@traitdef IsAnything{X}
@traitimpl IsAnything{X} <- (x->true)(X)

"Trait which contains no types"
@traitdef IsNothing{X}
@traitimpl IsNothing{X} <- (x->false)(X)

"Trait of all isbits-types"
@traitdef IsBits{X}
Base.@pure _isbits(X) = isbitstype(X)
@traitimpl IsBits{X} <- _isbits(X)

"Trait of all immutable types"
@traitdef IsImmutable{X}
if VERSION >= v"1.7.0-DEV.1279"
    Base.@pure _isimmutable(X) = !(Base.ismutabletype(X))
else
    Base.@pure _isimmutable(X) = !X.mutable
end
@traitimpl IsImmutable{X}  <- _isimmutable(X)

"Trait of all callable objects"
@traitdef IsCallable{X}
@traitimpl IsCallable{X} <- (X->(X<:Function ||  length(methods(X))>0))(X)

"Trait of all concrete types types"
@traitdef IsConcrete{X}
@traitimpl IsConcrete{X} <- isconcretetype(X)

Base.@deprecate_binding IsLeafType IsConcrete

"Types which have contiguous memory layout"
@traitdef IsContiguous{X} # https://github.com/JuliaLang/julia/issues/10889
@traitimpl IsContiguous{X} <- Base.iscontiguous(X)

"Array indexing trait."
@traitdef IsIndexLinear{X} # https://github.com/JuliaLang/julia/pull/8432
function isindexlinear(X)
    if IndexStyle(X)==IndexLinear()
        return true
    elseif  IndexStyle(X)==IndexCartesian()
        return false
    else
        error("Not recognized")
    end
end
@traitimpl IsIndexLinear{X} <- isindexlinear(X)

Base.@deprecate_binding IsFastLinearIndex IsIndexLinear

# TODO
## @traitdef IsArray{X} # use for any array like type in the sense of container
##                   # types<:AbstractArray are automatically part

## @traitdef IsMartix{X} # use for any LinearOperator
##                    # types<:AbstractArray are automatically part

"""
Trait of all iterator types.

NOTE: using this will lead to dynamic dispatch, see
https://github.com/mauro3/SimpleTraits.jl/issues/40 for context.
"""
@traitdef IsIterator{X}
function SimpleTraits.trait(::Type{IsIterator{X}}) where {X}
    hasmethod(iterate, Tuple{X}) ? IsIterator{X} : Not{IsIterator{X}}
end

end # module


text:
const FacElemQ = Union{FacElem{fmpq, FlintRationalField}, FacElem{fmpz, FlintIntegerRing}}

@doc Markdown.doc"""
    evaluate(x::FacElem{fmpq}) -> fmpq
    evaluate(x::FacElem{fmpz}) -> fmpz

Expands or evaluates the factored element, i.e. actually computes the
the element.
Works by first obtaining a simplified version of the power product
into coprime base elements.
"""
function evaluate(x::FacElem{fmpq})
  return evaluate_naive(simplify(x))
end

function evaluate(x::FacElem{fmpz})
  return evaluate_naive(simplify(x))
end
@doc Markdown.doc"""
    simplify(x::FacElem{fmpq}) -> FacElem{fmpq}
    simplify(x::FacElem{fmpz}) -> FacElem{fmpz}

Simplfies the factored element, i.e. arranges for the base to be coprime.
"""
function simplify(x::FacElem{fmpq})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify(x::FacElem{fmpz})
  y = deepcopy(x)
  simplify!(y)
  return y
end

function simplify!(x::FacElem{fmpq})
  if length(x.fac) <= 1
    return nothing
  end
  cp = vcat([denominator(y) for (y, v) in x if !iszero(v)], [numerator(y) for (y, v) in x if !iszero(v)])
  ev = Dict{fmpq, fmpz}()
  if isempty(cp)
    ev[fmpq(1)] = 0
    x.fac = ev
    return nothing
  end
  cp = coprime_base(cp)
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for (b, vb) in x
      if !iszero(vb)
        v += valuation(b, abs(p))*vb
      end
    end
    if v != 0
      ev[fmpq(abs(p))] = v
    end
  end
  f = b -> b < 0 && isodd(x.fac[b]) ? -1 : 1
  s = prod((f(v) for v in base(x)))
  if s == -1
    ev[fmpq(-1)] = 1
  else
    if length(ev) == 0
      ev[fmpq(1)] = 0
    end
  end
  x.fac = ev
  return nothing
end

function simplify!(x::FacElem{fmpz})
  if length(x.fac) == 0
    x.fac[fmpz(1)] = 0
    return
  end
  if length(x.fac) <= 1
    k,v = first(x.fac)
    if isone(k)
      x.fac[k] = 0
    elseif k == -1
      if isodd(v)
        x.fac[k] = 1
      else
        delete!(x.fac, k)
        x.fac[fmpz(1)] = 0
      end
    end
    return
  end
  cp = coprime_base(collect(base(x)))
  ev = Dict{fmpz, fmpz}()
  for p = cp
    if p == 1 || p == -1
      continue
    end
    v = fmpz(0)
    for (b, vb) in x
      v += valuation(b, abs(p))*vb
    end
    if v < 0
      throw(DomainError(v, "Negative valuation in simplify!"))
    end
    if v != 0
      ev[abs(p)] = v
    end
  end
  f = b -> b < 0 && isodd(x.fac[b]) ? -1 : 1
  s = prod(f(v) for v in base(x))
  if s == -1
    ev[-1] = 1
  else
    if length(ev) == 0
      ev[fmpz(1)] = 0
    end
  end
  x.fac = ev
  nothing
end

@doc Markdown.doc"""
    isone(x::FacElem{fmpq}) -> Bool
    isone(x::FacElem{fmpz}) -> Bool
Tests if $x$ represents $1$ without an evaluation.
"""
function isone(x::FacElem{fmpq})
  y = simplify(x)
  return all(iszero, values(y.fac)) || all(isone, keys(y.fac))
end

function isone(x::FacElem{fmpz})
  y = simplify(x)
  return all(iszero, values(y.fac)) || all(isone, keys(y.fac))
end

@doc Markdown.doc"""
    factor_coprime(x::FacElem{fmpz}) -> Fac{fmpz}
Computed a partial factorisation of $x$, ie. writes $x$ as a product
of pariwise coprime integers.
"""
function factor_coprime(x::FacElem{fmpz})
  x = deepcopy(x)
  simplify!(x)
  d = Dict(abs(p) => Int(v) for (p,v) = x.fac)
  if haskey(d, fmpz(-1))
    delete!(d, fmpz(-1))
    return Fac(fmpz(-1), d)
  else
    return Fac(fmpz(1), d)
  end
end

function abs(A::FacElemQ)
  B = empty(A.fac)
  for (k,v) = A.fac
    ak = abs(k)
    add_to_key!(B, ak, v)
  end
  if length(B) == 0
    return FacElem(Dict(one(base_ring(A)) => fmpz(1)))
  end
  return FacElem(B)
end

function ==(A::T, B::T) where {T <: FacElemQ}
  C = A*B^-1
  simplify!(C)
  return isone(C)
end

function isone(A::FacElemQ)
  C = simplify(A)
  return all(iszero, values(C.fac)) || all(isone, keys(C.fac))
end


text:
export Area, SOC, Linearize, Formulation, Unary, Binary,
	   BinaryGray, BinaryBlack, Partition4Bit,
	   Metric, Lowerbounds, Objective, CutCount,
	   base_model

abstract Metric
immutable Lowerbounds <: Metric end
immutable Objective   <: Metric end

immutable CutCount{R}
	num_rounds::Int
	metric::Metric

	CutCount() = new(1, Objective())
	CutCount(n::Int,metric::Metric) = new(n,metric)
end

abstract Area
immutable SOC <: Area end
immutable Linearize <: Area
	ε::Float64
end
Linearize() = Linearize(0.005)

abstract Formulation{A<:Area}
immutable Unary{A<:Area} <: Formulation{A}
	area::A
end
abstract Binary{A<:Area} <: Formulation{A}
immutable BinaryGray{A<:Area} <: Binary{A}
	area::A
end
immutable BinaryBlack{A<:Area} <: Binary{A}
	area::A
end
immutable Partition4Bit{A<:Area} <: Formulation{A}
	area::A
end

(::Type{F}){F<:Formulation}() = F(SOC())

get_metrics(::Objective,   prob::Problem, i, j) = prob.c[i,j]
get_metrics(::Lowerbounds, prob::Problem, i, j) = prob.wlb[i] + prob.wlb[j] + prob.hlb[i] + prob.hlb[j]

function binary_variables(model::Model, N, ::Unary; redundant=false)
	@variable(model, v[i=1:N,j=1:N,1:2;i!=j], Bin)
	# @variable(model, v[1:N,1:N,1:2], Bin)
	@constraint(model, c1[i=1:N,j=(i+1):N], v[i,j,1] + v[j,i,1] + v[i,j,2] + v[j,i,2] == 1)

	arrx, arry = Array(Variable, N, N), Array(Variable, N, N)
	for i in 1:N, j in (i+1):N
		arrx[i,j] = v[i,j,1]
		arrx[j,i] = v[j,i,1]
		arry[i,j] = v[i,j,2]
		arry[j,i] = v[j,i,2]
	end
	return v, arrx, arry
	# return v, v[:,:,1], v[:,:,2]
end

function binary_variables(model::Model, N, ::Partition4Bit; redundant=false)
	@variable(model, v[1:N,1:N,1:2], Bin)
	@constraint(model, c1[i=1:N,j=(i+1):N], v[i,j,1] + v[i,j,2] + v[j,i,1] + v[j,i,2] ≥ 1)
	@constraint(model, c1[i=1:N,j=(i+1):N,s=1:2], v[i,j,s] + v[j,i,s] ≤ 1)

	return v, v[:,:,1], v[:,:,2]
end

##############
# Implements the binary (logarithmic) gray encoding; that is, for each pair of boxes,
# the binary code assigned to the two polyhedra in direction x (resp. y) differ by
# a single bit. This assignment is unique up to symmetry, and we choose
#		i j <= v = (0,0)
#		j i <= v = (1,0)
# 		 j
#		 i  <= v = (0,1)
# 		 i
#		 j  <= v = (1,1)
# In other words, the second bit encodes which direction the separation occurs, while
# the first encodes which box proceeds the other in that direction.
#
# We also return "helper" collections arrˣ and arrʸ such that arrˢ[i,j] = 1 iff
# i precedes j in direction s, and arrˢ[i,j] ≤ 0 otherwise. This can be encoded as an
# affine map ℒᵍ that takes {0,1}² ⟶ R⁴ in the way described above; in particular, this
# means that any valid inequality for the unary encoding of the form aᵀx + bᵀv ≤ c
# for b ≥ 0 can be transformed to a valid inequality for the binary gray encoding.

ℒ{B<:BinaryGray}(v, ::Type{B}) = begin
	[-v[1] - v[2] + 1
	  v[1] - v[2]
	 -v[1] + v[2]
	  v[1] + v[2] - 1]
end

function binary_variables{B<:BinaryGray}(model::Model, N, ::B; redundant=false)
	if redundant
		@variable(model, z⁺[1:N,1:N], Bin)
		@variable(model, z⁻[1:N,1:N], Bin)
	else
		@variable(model, z⁺[i=1:N,(i+1):N], Bin)
		@variable(model, z⁻[i=1:N,(i+1):N], Bin)
	end
	arrˣ = Array(AffExpr, N, N)
	arrʸ = Array(AffExpr, N, N)
	for i in 1:N
		for j in (i+1):N
			if redundant
				@constraint(model, z⁺[i,j] + z⁺[j,i] == 1)
				@constraint(model, z⁻[i,j] == z⁻[j,i])
			end
			val = ℒ([z⁺[i,j],z⁻[i,j]], B)
			arrˣ[i,j] = val[1]
			arrˣ[j,i] = val[2]
			arrʸ[i,j] = val[3]
			arrʸ[j,i] = val[4]
		end
	end
	return (z⁺,z⁻), arrˣ, arrʸ
end

########
#		i j <= v = (0,0)
#		j i <= v = (1,1)
# 		 j
#		 i  <= v = (1,0)
# 		 i
#		 j  <= v = (0,1)

ℒ{B<:BinaryBlack}(v, ::Type{B}) = begin
	[-v[1] - v[2] + 1
	  v[1] + v[2] - 1
	  v[1] - v[2]
	 -v[1] + v[2]]
end

function binary_variables{B<:BinaryBlack}(model::Model, N, ::B; redundant=true)
	if redundant
		@variable(model, z⁺[1:N,1:N], Bin)
		@variable(model, z⁻[1:N,1:N], Bin)
		for i in 1:N, j in (i+1):N
			@constraint(model, z⁺[i,j] + z⁺[j,i] == 1)
			@constraint(model, z⁻[i,j] + z⁻[j,i] == 1)
		end
	else
		@variable(model, tmp⁺[i=1:N,(i+1):N], Bin)
		@variable(model, tmp⁻[i=1:N,(i+1):N], Bin)
		z⁺ = Array(AffExpr, N, N)
		z⁻ = Array(AffExpr, N, N)
		for i in 1:N, j in (i+1):N
			z⁺[i,j] =     tmp⁺[i,j]
			z⁺[j,i] = 1 - tmp⁺[i,j]
			z⁻[i,j] =     tmp⁻[i,j]
			z⁻[j,i] = 1 - tmp⁻[i,j]
		end
	end
	arrˣ = Array(AffExpr, N, N)
	arrʸ = Array(AffExpr, N, N)
	for i in 1:N, j in (i+1):N
		val = ℒ([z⁺[i,j],z⁻[i,j]], B)
		arrˣ[i,j] = val[1]
		arrˣ[j,i] = val[2]
		arrʸ[i,j] = val[3]
		arrʸ[j,i] = val[4]
	end
	return (z⁺,z⁻), arrˣ, arrʸ
end

function base_model{F<:Formulation}(model::Model, prob::Problem, form::F; redundant=false, symbreak=true)
	N   = prob.N
	lbˣ = prob.wlb
	lbʸ = prob.hlb
	ubˣ = prob.wub
	ubʸ = prob.hub
	p   = prob.c
	α   = prob.area
	Lˣ  = prob.W
	Lʸ  = prob.H

	@variable(model, 0 ≤ cˣ[i=1:N] ≤ Lˣ)
	@variable(model, 0 ≤ cʸ[i=1:N] ≤ Lʸ)
	@variable(model, lbˣ[i] ≤ ℓˣ[i=1:N] ≤ ubˣ[i])
	@variable(model, lbʸ[i] ≤ ℓʸ[i=1:N] ≤ ubʸ[i])

	@variable(model, tˣ[i=1:N,(i+1):N] ≥ 0)
	@variable(model, tʸ[i=1:N,(i+1):N] ≥ 0)
	dˣ = Array(Variable, N, N)
	dʸ = Array(Variable, N, N)
	for i in 1:N
		for j in (i+1):N
			dˣ[i,j] = tˣ[i,j]
			dˣ[j,i] = tˣ[i,j]
			dʸ[i,j] = tʸ[i,j]
			dʸ[j,i] = tʸ[i,j]
		end
	end

	@constraints(model, begin
		c3[i=1:N,j=(i+1):N], dˣ[i,j] ≥ cˣ[i] - cˣ[j]
		c4[i=1:N,j=(i+1):N], dˣ[i,j] ≥ cˣ[j] - cˣ[i]
		c5[i=1:N,j=(i+1):N], dʸ[i,j] ≥ cʸ[i] - cʸ[j]
		c6[i=1:N,j=(i+1):N], dʸ[i,j] ≥ cʸ[j] - cʸ[i]
	end)
	@objective(model, Min, sum(p[i,j]*(dˣ[i,j]+dʸ[i,j]) for i=1:N, j=(i+1):N))

	v, zˣ, zʸ = binary_variables(model, N, form; redundant=redundant)

	vars = Variables{F}(cˣ, cʸ, ℓˣ, ℓʸ, dˣ, dʸ, v, zˣ, zʸ)

	add_stay_in_the_box(model, prob, form, vars)
	add_area(model, prob, form, vars)

	add_nonoverlap(model, prob, form, vars)

	if symbreak
		# symmetry breaking: box 1 sits above and to the right of box 2
		maxval,linidx = findmax(p)
		ii, jj = ind2sub(size(p),linidx)
		@assert 1 ≤ ii ≤ N
		@assert 1 ≤ jj ≤ N
		@assert p[ii,jj] == maxval
		@assert ii != jj
		@constraints(model, begin
			cˣ[ii] ≤ cˣ[jj]
			cʸ[ii] ≤ cʸ[jj]
			zˣ[jj,ii] ≤ 0
			zʸ[jj,ii] ≤ 0
			(cˣ[jj]-cˣ[ii]) + (cʸ[jj]-cʸ[ii]) ≥ 0.5min(lbˣ[ii]+lbˣ[jj], lbʸ[ii]+lbʸ[jj])
		end)
	end

	return vars
end

function add_area(model::Model, prob::Problem, ::Formulation{SOC}, vars::Variables)
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	@assert (N = length(ℓˣ)) == length(ℓʸ)
	α = prob.area
	@constraint(model, c[i=1:N], ℓˣ[i]*ℓʸ[i] ≥ α[i])
end

function add_area(model::Model, prob::Problem, form::Formulation{Linearize}, vars::Variables)
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	@assert (N = length(ℓˣ)) == length(ℓʸ)
	lbˣ, lbʸ = prob.wlb, prob.hlb
	ubˣ, ubʸ = prob.wub, prob.hub
	α = prob.area
	ε = form.area.ε
    for i in 1:N
        Cᵢ = iceil(log(ubˣ[i]/lbˣ[i]) / log((1+√(ε))/(1-√(ε))))
        fac = (ubˣ[i]/lbˣ[i])^(1/Cᵢ)
        bnd = lbˣ[i]
        for k = 0:Cᵢ
            bnd *= fac
            @constraint(model, -ℓʸ[i] - α[i]/bnd^2*ℓˣ[i] <= -2*α[i]/bnd)
        end
    end

end

function add_nonoverlap(model::Model, prob::Problem, form::Formulation, vars::Variables)
	N = prob.N
	Lˣ, Lʸ = prob.W, prob.H
	cˣ, cʸ = vars.cˣ, vars.cʸ
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	for i in 1:N, j in 1:N
		i == j && continue
		@constraints(model, begin
			cˣ[i] - cˣ[j] + 0.5*(ℓˣ[i]+ℓˣ[j]) ≤ Lˣ*(1 - zˣ[i,j])
			cʸ[i] - cʸ[j] + 0.5*(ℓʸ[i]+ℓʸ[j]) ≤ Lʸ*(1 - zʸ[i,j])
		end)
	end
end

function add_nonoverlap{B<:Partition4Bit}(model::Model, prob::Problem, form::B, vars::Variables)
	N = prob.N
	Lˣ, Lʸ = prob.W, prob.H
	lbˣ, lbʸ = prob.wlb, prob.hlb
	cˣ, cʸ = vars.cˣ, vars.cʸ
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	zˣ, zʸ = vars.zˣ, vars.zʸ
	for i in 1:N, j in (i+1):N
		@constraints(model, begin
			cˣ[i] - cˣ[j] + 0.5*(ℓˣ[i]+ℓˣ[j]) ≤ Lˣ*(1-zˣ[i,j])
			cˣ[j] - cˣ[i] + 0.5*(ℓˣ[i]+ℓˣ[j]) ≤ Lˣ*(1-zˣ[j,i])
			cʸ[i] - cʸ[j] + 0.5*(ℓʸ[i]+ℓʸ[j]) ≤ Lʸ*(1-zʸ[i,j])
			cʸ[j] - cʸ[i] + 0.5*(ℓʸ[i]+ℓʸ[j]) ≤ Lʸ*(1-zʸ[j,i])

			cˣ[i] + 0.5ℓˣ[i] + Lˣ*zˣ[i,j] ≥ cˣ[j] - 0.5ℓˣ[j] + (lbˣ[i]+lbˣ[j])*(zˣ[i,j]+zˣ[j,i])
			cˣ[j] + 0.5ℓˣ[j] + Lˣ*zˣ[j,i] ≥ cˣ[i] - 0.5ℓˣ[i] + (lbˣ[i]+lbˣ[j])*(zˣ[i,j]+zˣ[j,i])
			cʸ[i] + 0.5ℓʸ[i] + Lʸ*zʸ[i,j] ≥ cʸ[j] - 0.5ℓʸ[j] + (lbʸ[i]+lbʸ[j])*(zʸ[i,j]+zʸ[j,i])
			cʸ[j] + 0.5ℓʸ[j] + Lʸ*zʸ[j,i] ≥ cʸ[i] - 0.5ℓʸ[i] + (lbʸ[i]+lbʸ[j])*(zʸ[i,j]+zʸ[j,i])
		end)
	end
end

function add_stay_in_the_box(model::Model, prob::Problem, f::Formulation, vars::Variables)
	cˣ, cʸ = vars.cˣ, vars.cʸ
	ℓˣ, ℓʸ = vars.ℓˣ, vars.ℓʸ
	@assert (N = prob.N) == length(cˣ) == length(cʸ) == length(ℓˣ) == length(ℓʸ)
	Lˣ, Lʸ = prob.W, prob.H
	for i in 1:N
		@constraints(model, begin
			cˣ[i] ≥      0.5ℓˣ[i]
			cˣ[i] ≤ Lˣ - 0.5ℓˣ[i]
			cʸ[i] ≥      0.5ℓʸ[i]
			cʸ[i] ≤ Lʸ - 0.5ℓʸ[i]
		end)
	end
end


text:
solve(i::BudgetedBipartiteMatchingInstance{T}, ::LagrangianAlgorithm, ε; kwargs...) where T = matching_hungarian_budgeted_lagrangian_search(i, ε; kwargs...)
solve(i::BudgetedBipartiteMatchingInstance{T}, ::LagrangianRefinementAlgorithm; kwargs...) where T = matching_hungarian_budgeted_lagrangian_refinement(i; kwargs...)
solve(i::BudgetedBipartiteMatchingInstance{T}, ::IteratedLagrangianRefinementAlgorithm; kwargs...) where T = matching_hungarian_budgeted_lagrangian_approx_half(i; kwargs...)

approximation_term(::BudgetedBipartiteMatchingInstance{T}, ::LagrangianAlgorithm) where T = NaN
approximation_ratio(::BudgetedBipartiteMatchingInstance{T}, ::LagrangianAlgorithm) where T = NaN

approximation_term(i::BudgetedBipartiteMatchingInstance{T}, ::LagrangianRefinementAlgorithm) where T = maximum(values(i.matching.rewards))
approximation_ratio(::BudgetedBipartiteMatchingInstance{T}, ::LagrangianRefinementAlgorithm) where T = NaN

approximation_term(::BudgetedBipartiteMatchingInstance{T}, ::IteratedLagrangianRefinementAlgorithm) where T = NaN
approximation_ratio(::BudgetedBipartiteMatchingInstance{T}, ::IteratedLagrangianRefinementAlgorithm) where T = 0.5

struct BudgetedBipartiteMatchingLagrangianSolution{T, U} <: BudgetedBipartiteMatchingSolution{T, U}
  # Used to store important temporary results from solving the Lagrangian dual.
  instance::BudgetedBipartiteMatchingInstance{T, U}
  solution::Vector{Edge{T}}
  λ::Float64 # Optimum dual multiplier.
  value::Float64 # Optimum value of the dual problem (i.e. with penalised constraint).
  λmax::Float64 # No dual value higher than this is useful (i.e. they all yield the same solution).
end

struct SimpleBudgetedBipartiteMatchingSolution{T, U} <: BudgetedBipartiteMatchingSolution{T, U}
  instance::BudgetedBipartiteMatchingInstance{T, U}
  solution::Vector{Edge{T}}
  value::Float64

  function SimpleBudgetedBipartiteMatchingSolution(instance::BudgetedBipartiteMatchingInstance{T, U}, solution::Vector{Edge{T}}) where {T, U}
    if length(solution) > 0
      return new{T, U}(instance, solution, _budgeted_bipartite_matching_compute_value(instance, solution))
    else
      return new{T, U}(instance, solution, -Inf)
    end
  end

  function SimpleBudgetedBipartiteMatchingSolution(instance::BudgetedBipartiteMatchingInstance{T, U}) where {T, U}
    # No feasible solution.
    return new{T, U}(instance, edgetype(instance.matching.graph)[], -Inf)
  end
end

function matching_hungarian_budgeted_lagrangian(i::BudgetedBipartiteMatchingInstance{T, U}, λ::Float64) where {T, U}
  # Solve the subproblem for one value of the dual multiplier λ:
  #     l(λ) = \max_{x matching} (rewards + λ weights) x - λ budget.
  bmi_rewards = Dict{Edge{T}, Float64}(e => i.matching.reward[e] + λ * i.weight[e] for e in keys(i.matching.reward))
  bmi = BipartiteMatchingInstance(i.matching, bmi_rewards)
  bmi_sol = matching_hungarian(bmi)
  return _budgeted_bipartite_matching_compute_value(bmi, bmi_sol.solution) - λ * i.budget, bmi_sol.solution
end

function matching_hungarian_budgeted_lagrangian_search(i::BudgetedBipartiteMatchingInstance{T, U}, ε::Float64) where {T, U}
  # Approximately solve the problem \min_{l ≥ 0} l(λ), where
  #     l(λ) = \max_{x smatching} (rewards + λ weights) x - λ budget.
  # This problem is the Lagrangian dual of the budgeted maximum spanning-tree problem:
  #     \max_{x matching} rewards x  s.t.  weights x >= budget.
  # This algorithm provides no guarantee on the optimality of the solution.

  # TODO: generalise this function? Highly similar to the one for spanning trees.

  # Initial set of values for λ. The optimum is guaranteed to be contained in this interval.
  weights_norm_inf = maximum(values(i.weight)) # Maximum weight.
  m = min(i.matching.n_left, i.matching.n_right) # Maximum number of items in a solution. Easy to compute for a matching!
  λmax = Float64(weights_norm_inf * m + 1)

  λlow = 0.0
  λhigh = λmax

  # Perform a golden-ratio search.
  while (λhigh - λlow) > ε
    λmidlow = λhigh - (λhigh - λlow) / MathConstants.φ
    λmidhigh = λlow + (λhigh - λlow) / MathConstants.φ
    vmidlow, _ = matching_hungarian_budgeted_lagrangian(i, λmidlow)
    vmidhigh, _ = matching_hungarian_budgeted_lagrangian(i, λmidhigh)

    if vmidlow < vmidhigh
      λhigh = λmidhigh
      vhigh = vmidhigh
    else
      λlow = λmidlow
      vlow = vmidlow
    end
  end

  vlow, bmlow = matching_hungarian_budgeted_lagrangian(i, λlow)
  vhigh, bmhigh = matching_hungarian_budgeted_lagrangian(i, λhigh)
  if vlow < vhigh
    return BudgetedBipartiteMatchingLagrangianSolution(i, bmlow, λlow, vlow, λmax)
  else
    return BudgetedBipartiteMatchingLagrangianSolution(i, bmhigh, λhigh, vhigh, λmax)
  end
end

function matching_hungarian_budgeted_lagrangian_refinement(i::BudgetedBipartiteMatchingInstance{T, U};
                                                           ε::Float64=1.0e-3, ζ⁻::Float64=0.2, ζ⁺::Float64=5.0,
                                                           stalling⁻::Float64=1.0e-5) where {T, U}
  # Approximately solve the following problem:
  #     \max_{x matching} rewards x  s.t.  weights x >= budget
  # This algorithm provides an additive approximation to this problem. If x* is the optimum solution and x~ the one
  # returned by this algorithm,
  #     weights x* >= budget   and   weights x~ >= budget                 (the returned solution is feasible)
  #     rewards x~ >= rewards x* - 2 \max{e edge} reward[e]               (additive approximation)
  # This is based on http://people.idsia.ch/~grandoni/Pubblicazioni/BBGS08ipco.pdf, Lemma 4.

  # Check assumptions.
  if ζ⁻ >= 1.0
    error("ζ⁻ must be strictly less than 1.0: the dual multiplier λ is multiplied by ζ⁻ to reach an infeasible solution by less penalising the budget constraint.")
  end
  if ζ⁺ <= 1.0
    error("ζ⁺ must be strictly greater than 1.0: the dual multiplier λ is multiplied by ζ⁺ to reach a feasible solution by penalising more the budget constraint.")
  end

  # Ensure the problem is feasible by only considering the budget constraint.
  feasible_rewards = Dict{Edge{T}, Float64}(e => i.weight[e] for e in keys(i.weight))
  feasible_instance = BipartiteMatchingInstance(i.matching, feasible_rewards)
  feasible_solution = matching_hungarian(feasible_instance)
  if _budgeted_bipartite_matching_compute_value(feasible_instance, feasible_solution.solution) < i.budget
    # By maximising the left-hand side of the budget constraint, impossible to reach the target budget. No solution!
    return SimpleBudgetedBipartiteMatchingSolution(i)
  end

  # Solve the Lagrangian relaxation to optimality.
  lagrangian = matching_hungarian_budgeted_lagrangian_search(i, ε)
  λ0, v0, bm0 = lagrangian.λ, lagrangian.value, lagrangian.solution
  λmax = lagrangian.λmax
  b0 = _budgeted_bipartite_matching_compute_weight(i, bm0) # Budget consumption of this first solution.

  # If already respecting the budget constraint exactly, done!
  if b0 == i.budget
    return SimpleBudgetedBipartiteMatchingSolution(i, bm0)
  end

  # Find two solutions: one above the budget x⁺ (i.e. respecting the constraint), the other not x⁻.
  x⁺, x⁻ = nothing, nothing

  λi = λ0
  if b0 > i.budget
    x⁺ = bm0
    @assert _budgeted_bipartite_matching_compute_weight(i, x⁺) > i.budget

    stalling = false
    while true
      # Penalise less the constraint: it should no more be satisfied.
      λi *= ζ⁻
      _, sti = matching_hungarian_budgeted_lagrangian(i, λi)
      if _budgeted_bipartite_matching_compute_weight(i, sti) < i.budget
        x⁻ = sti
        break
      end

      # Is the process stalling?
      if λi <= stalling⁻
        stalling = true
        break
      end
    end

    # Specific handling of stallings.
    if stalling # First test: don't penalise the constraint at all.
      _, sti = matching_hungarian_budgeted_lagrangian(i, 0.0)
      new_budget = _budgeted_bipartite_matching_compute_weight(i, sti)
      if new_budget < i.budget
        x⁻ = sti
        stalling = false
      end
    end

    if stalling # Second test: minimise the left-hand side of the budget constraint, in hope of finding a feasible solution.
      # This process is highly similar to the computation of feasible_solution, but with a reverse objective function.
      infeasible_rewards = Dict{Edge{T}, Float64}(e => - i.weight[e] for e in keys(i.weight))
      infeasible_solution = matching_hungarian(BipartiteMatchingInstance(i.matching.graph, infeasible_rewards)).solution

      if _budgeted_bipartite_matching_compute_weight(i, infeasible_solution) < i.budget
        x⁻ = infeasible_solution
        stalling = false
      end
    end

    if stalling # Third: decide there is no solution strictly below the budget. No refinement is possible.
      # As x⁺ is feasible, return it.
      return SimpleBudgetedBipartiteMatchingSolution(i, x⁺)
    end
  else
    x⁻ = bm0
    @assert _budgeted_bipartite_matching_compute_weight(i, x⁻) < i.budget

    while true
      # Penalise more the constraint: it should become satisfied at some point.
      λi *= ζ⁺
      _, sti = matching_hungarian_budgeted_lagrangian(i, λi)
      if _budgeted_bipartite_matching_compute_weight(i, sti) >= i.budget
        x⁺ = sti
        break
      end

      # Is the process stalling? If so, reuse feasible_solution, which is guaranteed to be feasible.
      if λi >= λmax
        x⁺ = feasible_solution.solution
        break
      end
    end
  end

  # Normalise the solutions: the edges point from the left part to the right part of the bipartite graph.
  sort_edge(e::Edge{T}) where T = (i.matching.partition[src(e)] == 1) ? e : reverse(e)
  x⁺ = [sort_edge(e) for e in x⁺]
  x⁻ = [sort_edge(e) for e in x⁻]

  already_seen_solutions = Set{Set{Edge{T}}}()
  push!(already_seen_solutions, Set(x⁺))
  push!(already_seen_solutions, Set(x⁻))

  # Iterative refinement. Stop as soon as there is a difference of at most one edge between the two solutions.
  while _solution_symmetric_difference_size(x⁺, x⁻) > 4
    # Enforce the loop invariant.
    @assert x⁺ !== nothing
    @assert x⁻ !== nothing
    @assert _budgeted_bipartite_matching_compute_weight(i, x⁺) >= i.budget # Feasible.
    @assert _budgeted_bipartite_matching_compute_weight(i, x⁻) < i.budget # Infeasible.

    # Switch elements from one solution to another: use the XOR between the two solutions as an augmenting path/cycle.
    only_in_x⁺, only_in_x⁻ = _solution_symmetric_difference(x⁺, x⁻)
    xor = vcat(only_in_x⁺, only_in_x⁻)

    current_edge = first(xor)
    current_vertex = src(current_edge)
    all_vertices = Set{T}(current_vertex)
    xor = xor[2:end]
    new_x = copy(x⁺)

    while true
      # Exit conditions:
      # - If the XOR is a path: consume all edges; at some point, will no more
      #   be possible to find the next one
      # - Otherwise, once a node is met twice, break.

      # If the current edge is in the solution, remove it. Otherwise, add it.
      if current_edge in new_x
        filter!(e -> e != current_edge, new_x)
      else
        push!(new_x, current_edge)
      end

      # Go to the next edge, following a cycle or a path.
      if current_vertex == src(current_edge)
        current_vertex = dst(current_edge)
      else
        current_vertex = src(current_edge)
      end
      current_edge_idx = findfirst(e -> src(e) == current_vertex || dst(e) == current_vertex, xor)
      if isnothing(current_edge_idx)
        break
      end
      current_edge = xor[current_edge_idx]
      deleteat!(xor, current_edge_idx)

      # Check for a cycle.
      if current_vertex in all_vertices
        break
      end
      push!(all_vertices, current_vertex)
    end

    # If this procedure cannot propose a new solution, it has stalled and will not recover.
    if Set(new_x) in already_seen_solutions
      break
    else
      push!(already_seen_solutions, Set(new_x))
    end

    # Replace one of the two solutions, depending on whether this solution is feasible (x⁺) or not (x⁻).
    if _budgeted_bipartite_matching_compute_weight(i, new_x) >= i.budget
      x⁺ = new_x
    else
      x⁻ = new_x
    end
  end

  # Done!
  return SimpleBudgetedBipartiteMatchingSolution(i, x⁺)
end

function _edge_any_end_match(e::Edge{T}, t::Edge{T}) where T
  return src(e) == src(t) || dst(e) == dst(t) || src(e) == dst(t) || dst(e) == src(t)
end

function matching_hungarian_budgeted_lagrangian_approx_half(i::BudgetedBipartiteMatchingInstance{T, U}; kwargs...) where {T, U}
  # Approximately solve the following problem:
  #     \max_{x matching} rewards x  s.t.  weights x >= budget
  # This algorithm provides a multiplicative approximation to this problem. If x* is the optimum solution and x~ the one
  # returned by this algorithm,
  #     weights x* >= budget   and   weights x~ >= budget                 (the returned solution is feasible)
  #     rewards x~ >= (1 - ε) rewards x*                                  (multiplicative approximation)
  # The parameter ε is not tuneable, but rather fixed to 1/2.

  # This is based on http://people.idsia.ch/~grandoni/Pubblicazioni/BBGS08ipco.pdf, Theorem 1.

  # If there are too few vertices, not much to do. The smallest side of the bipartite graph must have at least 4 vertices,
  # so that 4 of them can be fixed.
  if i.matching.n_left <= 4 || i.matching.n_right <= 4
    # TODO: Should rather perform an exhaustive exploration.
    return matching_hungarian_budgeted_lagrangian_refinement(i; kwargs...)
  end

  # For each combination of four distinct edges, force these four edges to be part of the solution and discard all edges with a higher value.
  # If selecting twice the same edge: no need to go further.
  # If any end of two selected edges are the same: this cannot lead to a feasible solution.
  # Both tests are implemented in _edge_any_end_match.
  best_sol = nothing
  for e1 in edges(i.matching.graph)
    for e2 in edges(i.matching.graph)
      if _edge_any_end_match(e1, e2)
        continue
      end

      for e3 in edges(i.matching.graph)
        if _edge_any_end_match(e1, e3) || _edge_any_end_match(e2, e3)
          continue
        end

        for e4 in edges(i.matching.graph)
          if _edge_any_end_match(e1, e4) || _edge_any_end_match(e2, e4) || _edge_any_end_match(e3, e4)
            continue
          end

          # Filter out the edges that have a higher value than any of these two edges. Give a very large reward to them both.
          cutoff = min(i.matching.reward[e1], i.matching.reward[e2], i.matching.reward[e3], i.matching.reward[e4])
          reward = copy(i.matching.reward)
          filter!(kv -> kv[2] <= cutoff, reward)
          filter!(kv -> ! _edge_any_end_match(kv[1], e1), reward)
          filter!(kv -> ! _edge_any_end_match(kv[1], e2), reward)
          filter!(kv -> ! _edge_any_end_match(kv[1], e3), reward)
          filter!(kv -> ! _edge_any_end_match(kv[1], e4), reward)
          reward[e1] = reward[e2] = reward[e3] = reward[e4] = prevfloat(Inf)

          if length(keys(reward)) == 4 # Nothing left? Skip.
            continue
          end

          graph = SimpleGraph(nv(i.matching.graph))
          for e in keys(reward)
            add_edge!(graph, e)
          end

          weight = Dict(e => i.weight[e] for e in keys(reward))

          # Solve this subproblem.
          bbmi = BudgetedBipartiteMatchingInstance(graph, reward, weight, i.budget)
          sol = matching_hungarian_budgeted_lagrangian_refinement(bbmi; kwargs...)

          # This subproblem is infeasible. Maybe it's because the overall problem is infeasible or just because too many
          # edges were removed.
          if length(sol.solution) == 0
            continue
          end

          # Impossible to have a feasible solution with these two edges, probably because of the budget constraint.
          # It's very unlikely that this happens, due to the checks already performed when looping.
          if ! (e1 in sol.solution) || ! (e2 in sol.solution) || ! (e3 in sol.solution) || ! (e4 in sol.solution)
            continue
          end

          # Only keep the best solution.
          if best_sol == nothing || _budgeted_bipartite_matching_compute_value(i, sol.solution) > _budgeted_bipartite_matching_compute_value(i, best_sol.solution)
            # sol's instance is the one used internally for the subproblems.
            best_sol = SimpleBudgetedBipartiteMatchingSolution(i, sol.solution)
          end
        end
      end
    end
  end

  if best_sol !== nothing
    return best_sol
  else
    return SimpleBudgetedBipartiteMatchingSolution(i)
  end
end


text:
using AsapSim

# Stdlib imports
using Test
using Random
using Profile

# External imports
import DataStructures
using MacroTools
using BenchmarkTools
using LightDES

# Include the TestPrograms package, which contains a collection of programs
# testing portions of the simulator.
include("TestPrograms/src/TestPrograms.jl")
using .TestPrograms

# Assembler tests
include("assembler/macro.jl")

# Fifo Tests
# include("cores/fifos/fifo.jl")

# IO Handlers
include("cores/io/io.jl")

# Core Component Tests
# include("cores/core/stall.jl")
include("cores/core/alu.jl")
include("cores/core/programs.jl")
include("cores/core/sort.jl")

include("sim/snakesort.jl")


text:
@testset "type system" begin
    using DSJulia.Solutions
    @test Solutions.bunchofnumbers_parser("4, 5.5, 8\n") ≈ 4 + 5.5 + 8

    @testset "Wizarding currency" begin

        @test Solutions.WizCur(1, 1, 1) == Solutions.WizCur(0, 17, 30)
        @test Solutions.galleons(Solutions.WizCur(0, 17, 0)) == 1
        @test Solutions.knuts(Solutions.WizCur(0, 17, 78)) == 78 % 29
        @test Solutions.moneyinknuts(Solutions.WizCur(1, 1, 1)) == 1 + 29 + 17 * 29
        @test Solutions.WizCur(1, 1, 1) + Solutions.WizCur(0, 17, 30) == Solutions.WizCur(2, 2, 2)
    end

    @testset "Vandermonde" begin
        using LinearAlgebra
        α = [1, 2, 3, 4]
        V = Solutions.Vandermonde(α, 4)
        @test V == [ones(Int, 4) α α.^2 α.^3 ]
        @test det(V) == det([ones(Int, 4) α α.^2 α.^3 ])
    end
end


text:
# replace all indices by a function of that index
function replaceindices(f, ex::Expr)
    if istensor(ex)
        if ex.head == :ref || ex.head == :typed_hcat
            if length(ex.args) == 1
                return ex
            elseif isa(ex.args[2], Expr) && ex.args[2].head == :parameters
                arg2 = ex.args[2]
                return Expr(ex.head, ex.args[1],
                            Expr(arg2.head, map(f, arg2.args)...),
                            (f(ex.args[i]) for i = 3:length(ex.args))...)
            else
                return Expr(ex.head, ex.args[1],
                            (f(ex.args[i]) for i = 2:length(ex.args))...)
            end
            return ex
        else #if ex.head == :typed_vcat
            arg2, arg3 = map((ex.args[2], ex.args[3])) do arg
                if isa(arg, Expr) && (arg.head == :row || arg.head == :tuple)
                    return Expr(arg.head, map(f, arg.args)...)
                else
                    return f(arg)
                end
            end
            return Expr(ex.head, ex.args[1], arg2, arg3)
        end
    else
        return Expr(ex.head, (replaceindices(f, e) for e in ex.args)...)
    end
end
replaceindices(f, ex) = ex

function normalizeindex(ex)
    if isa(ex, Symbol) || isa(ex, Int)
        return ex
    elseif isa(ex, Expr) && ex.head == prime && length(ex.args) == 1
        return Symbol(normalizeindex(ex.args[1]), "′")
    else
        error("not a valid index: $ex")
    end
end

# replace all tensor objects by a function of that object
function replacetensorobjects(f, ex::Expr)
    # first try to replace ex completely
    ex2 = f(ex, nothing, nothing)
    ex2 !== ex && return ex2
    if istensor(ex)
        obj, leftind, rightind = decomposetensor(ex)
        return Expr(ex.head, f(obj, leftind, rightind), ex.args[2:end]...)
    else
        return Expr(ex.head, (replacetensorobjects(f, e) for e in ex.args)...)
    end
end
replacetensorobjects(f, ex) = f(ex, nothing, nothing)

# expandconj: conjugate individual terms or factors instead of a whole expression
function expandconj(ex::Expr)
    if isgeneraltensor(ex) || isscalarexpr(ex)
        return ex
    elseif ex.head == :call && ex.args[1] == :conj
        @assert length(ex.args) == 2
        return conjexpr(expandconj(ex.args[2]))
    else
        return Expr(ex.head, map(expandconj, ex.args)...)
    end
end
expandconj(ex) = ex

function conjexpr(ex::Expr)
    if ex.head == :call && ex.args[1] == :conj
        return ex.args[2]
    elseif isgeneraltensor(ex) || isscalarexpr(ex)
        return Expr(:call, :conj, ex)
    elseif ex.head == :call && (ex.args[1] == :* || ex.args[1] == :+ || ex.args[1] == :-)
        return Expr(ex.head, ex.args[1], map(conjexpr, ex.args[2:end])...)
    elseif ex.head == :call && (ex.args[1] == :/ || ex.args[1] == :\)
        return Expr(ex.head, ex.args[1], map(conjexpr, ex.args[2:end])...)
    else
        error("cannot conjugate expression: $ex")
    end
end
conjexpr(ex::Number) = conj(ex)
conjexpr(ex::Symbol) = Expr(:call, :conj, ex)
conjexpr(ex) = ex

# explicitscalar: wrap all tensor expressions with zero output indices in scalar call
function explicitscalar(ex::Expr)
    ex = Expr(ex.head, map(explicitscalar, ex.args)...)
    if istensorexpr(ex) && isempty(getindices(ex))
        return Expr(:call, :scalar, ex)
    else
        return ex
    end
end
explicitscalar(ex) = ex

# extracttensorobjects: replace all tensor objects with newly generated symbols, and assign
# them before the expression and after the expression as necessary
function extracttensorobjects(ex)
    inputtensors = getinputtensorobjects(ex)
    outputtensors = getoutputtensorobjects(ex)
    newtensors = getnewtensorobjects(ex)
    existingtensors = unique!(vcat(inputtensors, outputtensors))
    alltensors = unique!(vcat(existingtensors, newtensors))
    tensordict = Dict{Any,Any}(a => gensym() for a in alltensors)
    pre = Expr(:block, [Expr(:(=), tensordict[a], a) for a in existingtensors]...)
    ex = replacetensorobjects((obj,leftind,rightind)->get(tensordict, obj, obj), ex)
    post = Expr(:block, [Expr(:(=), a, tensordict[a]) for a in newtensors]...)
    pre2 = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), pre)
    post2 = Expr(:macrocall, Symbol("@notensor"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), post)
    return Expr(:block, pre2, ex, post2)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "XzBuilder"
version = v"5.2.4"

# Collection of sources required to build XzBuilder
sources = [
    "https://tukaani.org/xz/xz-$(version).tar.xz" =>
    "9717ae363760dedf573dad241420c5fea86256b65bc21d2cf71b2b12f0544f4b",
    "https://tukaani.org/xz/xz-$(version)-windows.zip" =>
    "9a5163623f435b6fa0844b6b884babd6bf4f8d876ae2d8134deeb296afd49c61",
]

# Bash recipe for building across all platforms
script = raw"""
if [ ${target} = "x86_64-w64-mingw32" ]; then
    mkdir -p ${WORKSPACE}/destdir/bin/
    cp bin_x86-64/liblzma.dll ${WORKSPACE}/destdir/bin/
elif [ ${target} = "i686-w64-mingw32" ]; then
    mkdir -p ${WORKSPACE}/destdir/bin/
    cp bin_i686/liblzma.dll ${WORKSPACE}/destdir/bin/
else
    cd $WORKSPACE/srcdir/xz-*
    ./configure --prefix=${prefix} --host=${target}
    make
    make install
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "liblzma", :liblzma)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
using Base.Test
using PointClouds
using StaticArrays

include("cloud_test.jl")


text:
## Coordinate descent

x,y = [1;1.5]
for i=1:10
  global x,y
  x = (1-y)/2
  y = (1.5-x)/2
  @show [x,y]
end

## Animate the algorithm
using Plots
using LaTeXStrings
using Printf
pyplot() # nicer output
theme(:dark)

A = [2 1; 1 2]
b = [1,1.5]
quadratic_ls(A,b) = x-> 0.5*x'*A*x - x'*b
ezsurf(x, y, f) = begin
    X = repeat(x', length(y), 1)
    Y = repeat(y, 1, length(x))
    # Evaluate each f(x, y)
    Z = map((x,y) -> f([x,y]), X, Y)
    plot(x, y, Z, st=:surface)
end
ezsurf(-1.6:0.05:1.6,-1.6:0.05:1.6,
  quadratic_ls(A,b))
plot!(colorbar=false,size=(500,450))
plot!(camera=(170,60))
p = plot3d!(1, label="", marker = 2)
x,y = [1;1.5]
anim = @animate for i=1:20
  global x, y
  if i > 1
    if mod(i,2)==0
      x = (1-y)/2
    else
      y = (1.5-x)/2
    end
  end
  push!(p[1][2], x, y, quadratic_ls(A,b)([x,y]))
  title!(@sprintf("Iteration %3i", i))
end
gif(anim, "coordinate-descent-example.gif")
##
ezsurf(-1.5:0.05:1.5,-1.5:0.05:1.5,
  quadratic_ls(A,b))
x,y = [1;1.5]
plot3d!([x, (1-y)/2],[y, (1.5-x)/2],[quadratic_ls(A,b)([x,y]),
  quadratic_ls(A,b)([(1-y)/2,(1.5-x)/2])],marker=2)
plot!(camera=(170,30))


text:
using PSIS
using Test

@testset "plotting backend" begin
    PSIS._paretoshapeplot(::Val{:TestBackend}, args...; kwargs...) = 10
    PSIS._paretoshapeplot!(::Val{:TestBackend}, args...; kwargs...) = 11
    old_backend = PSIS.PLOTTING_BACKEND[]
    PSIS.plotting_backend!(:TestBackend)
    @test PSIS.PLOTTING_BACKEND[] == :TestBackend
    values = randn(100)
    @test paretoshapeplot(values) === 10
    @test paretoshapeplot!(values) === 11
    PSIS.plotting_backend!(old_backend)
end


text:
# algorithm

struct SinkhornStabilized{T<:Real} <: Sinkhorn
    absorb_tol::T
end

"""
    SinkhornStabilized(; absorb_tol::Real=1_000)

Construct a log-domain stabilized Sinkhorn algorithm with absorption tolerance `absorb_tol`
for solving an entropically regularized optimal transport problem.

# References

Schmitzer, B. (2019).
[Stabilized Sparse Scaling Algorithms for Entropy Regularized Transport Problems](https://doi.org/10.1137/16m1106018).
SIAM Journal on Scientific Computing, 41(3), A1443–A1481.
"""
function SinkhornStabilized(; absorb_tol::Real=1_000)
    absorb_tol > zero(absorb_tol) ||
        throw(ArgumentError("absorption tolerance `absorb_tol` must be positive"))
    return SinkhornStabilized(absorb_tol)
end

function Base.show(io::IO, alg::SinkhornStabilized)
    return print(
        io,
        "Log-domain stabilized Sinkhorn algorithm (absorption tolerance = ",
        alg.absorb_tol,
        ")",
    )
end

# caches

struct SinkhornStabilizedCache{U,V,KT}
    u::U
    v::V
    alpha::U
    beta::V
    K::KT
    Kv::U
end

function build_cache(
    ::Type{T},
    ::SinkhornStabilized,
    size2::Tuple,
    μ::AbstractVecOrMat,
    ν::AbstractVecOrMat,
    C::AbstractMatrix,
    ε::Real,
) where {T}
    # compute Gibbs kernel (has to be mutable for ε-scaling algorithm)
    K = similar(C, T, size(C)..., size2...)
    @. K = exp(-C / ε)

    # create and initialize dual potentials
    u = similar(μ, T, size(μ, 1), size2...)
    v = similar(ν, T, size(ν, 1), size2...)
    α = similar(u)
    β = similar(v)
    fill!(u, one(T))
    fill!(v, one(T))
    fill!(α, zero(T))
    fill!(β, zero(T))

    # cache for next iterate of `u`
    Kv = similar(u)

    return SinkhornStabilizedCache(u, v, α, β, K, Kv)
end

# interface

function prestep!(solver::SinkhornSolver{<:SinkhornStabilized}, iter::Int)
    # unpack
    absorb_tol = solver.alg.absorb_tol
    cache = solver.cache
    u = cache.u
    v = cache.v

    # absorption step
    if maximum(abs, u) > absorb_tol || maximum(abs, v) > absorb_tol
        @debug string(solver.alg) *
               " (" *
               string(iter) *
               "/" *
               string(maxiter) *
               ": absorbing `u` and `v` into `α` and `β`"

        # absorb `u` and `v` into `α` and `β`
        absorb!(solver)

        # re-compute cache for next iterate of `u`
        A_batched_mul_B!(cache.Kv, cache.K, v)
    end

    return nothing
end

# absorption step
function absorb!(solver::SinkhornSolver{<:SinkhornStabilized})
    # unpack
    ε = solver.eps
    cache = solver.cache
    u = cache.u
    v = cache.v
    α = cache.alpha
    β = cache.beta

    # update dual potentials
    α .+= ε .* log.(u)
    β .+= ε .* log.(v)

    # reset `α` and `β`
    fill!(u, one(eltype(u)))
    fill!(v, one(eltype(v)))

    # update kernel
    update_K!(solver)

    return nothing
end

# update kernel
function update_K!(solver::SinkhornSolver{<:SinkhornStabilized})
    cache = solver.cache
    α = add_singleton(cache.alpha, Val(2))
    β = add_singleton(cache.beta, Val(1))
    @. cache.K = exp(-(solver.C - α - β) / solver.eps)
    return nothing
end

# obtain plan
function plan(solver::SinkhornSolver{<:SinkhornStabilized})
    absorb!(solver)
    return copy(solver.cache.K)
end

# deprecations

"""
    sinkhorn_stabilized(μ, ν, C, ε; absorb_tol=1_000, kwargs...)

This method is deprecated, please use
```julia
sinkhorn(
    μ, ν, C, ε, SinkhornStabilized(; absorb_tol=absorb_tol); kwargs...
)
```
instead.

See also: [`sinkhorn`](@ref), [`SinkhornStabilized`](@ref)
"""
function sinkhorn_stabilized(
    μ, ν, C, ε; tol=nothing, atol=tol, return_duals=nothing, absorb_tol=1_000, kwargs...
)
    # deprecation warnings
    if tol !== nothing
        Base.depwarn(
            "keyword argument `tol` is deprecated, please use `atol` and `rtol`",
            :sinkhorn_stabilized,
        )
    end
    if return_duals !== nothing
        Base.depwarn(
            "keyword argument `return_duals` is deprecated and will be removed",
            :sinkhorn_stabilized,
        )
    end

    Base.depwarn(
        "`sinkhorn_stabilized` is deprecated, " *
        "please use `sinkhorn` with `SinkhornStabilized` instead",
        :sinkhorn_stabilized_epsscaling,
    )

    # construct algorithm
    algorithm = SinkhornStabilized(; absorb_tol=absorb_tol)

    # return "dual potentials" (in log space and scaled by 1/ε) if requested
    if return_duals !== nothing && return_duals
        solver = build_solver(μ, ν, C, ε, algorithm; atol=atol, kwargs...)
        solve!(solver)
        absorb!(solver)
        cache = solver.cache
        return cache.alpha, cache.beta
    end

    return sinkhorn(μ, ν, C, ε, algorithm; atol=atol, kwargs...)
end