text:
struct NullStep end
struct CartesianVIndex{N,T<:Tuple{Vararg{Union{Integer,NullStep},N}}} <: Base.AbstractCartesianIndex{N}
    I::T
    @inline CartesianVIndex(I::T) where {N, T <: Tuple{Vararg{Union{Integer,NullStep},N}}} = new{N,T}(I)
end
Base.length(::CartesianVIndex{N}) where {N} = N
ArrayInterface.known_length(::Type{<:CartesianVIndex{N}}) where {N} = N
Base.Tuple(i::CartesianVIndex) = getfield(i, :I)
function Base.:(:)(I::CartesianVIndex{N}, J::CartesianVIndex{N}) where {N}
   CartesianIndices(map((i,j) -> i:j, getfield(I, :I), getfield(J,:I)))
end
Base.@propagate_inbounds Base.getindex(I::CartesianVIndex, i) = getfield(I,:I)[i]
_ndim(::Type{<:Base.AbstractCartesianIndex{N}}) where {N} = N
@inline gesp(p::AbstractStridedPointer{T,N}, i::Tuple{CartesianVIndex{N}}) where {T,N} = gesp(p, getfield(getfield(i,1,false),:I))
# _ndim(::Type{<:AbstractArray{N}}) where {N} = N
@generated function CartesianVIndex(I::T) where {T <: Tuple{Vararg{Union{Integer,CartesianIndex,CartesianVIndex,NullStep}}}}
    iexpr = Expr(:tuple)
    Tp = T.parameters
    q = Expr(:block)
    for i in eachindex(Tp)
        I_i = Symbol(:I_, i)
        push!(q.args, Expr(:(=), I_i, Expr(:ref, :I, i)))
        Tpᵢ = Tp[i]
        if Tpᵢ <: Base.AbstractCartesianIndex
            for n in 1:_ndim(Tpᵢ)
                push!(iexpr.args, Expr(:ref, I_i, n))
            end
        # elseif Tpᵢ === NullStep
        #     push!(iexpr.args, :(Zero()))
        else
            push!(iexpr.args, I_i)
        end
    end
    push!(q.args, Expr(:call, :CartesianVIndex, iexpr))
    Expr(:block, Expr(:meta, :inline), Expr(:macrocall, Symbol("@inbounds"), LineNumberNode(@__LINE__, Symbol(@__FILE__)), q))
end

# @inline Base.CartesianIndex(I::Tuple{Vararg{Union{Integer,CartesianIndex,CartesianVIndex,StaticInt}}}) = CartesianVIndex(I)

@generated function _maybestaticfirst(a::Tuple{Vararg{Any,N}}) where {N}
    quote
        $(Expr(:meta,:inline))
        Base.Cartesian.@ntuple $N n -> maybestaticfirst(a[n])
    end
end
@generated function _maybestaticlast(a::Tuple{Vararg{Any,N}}) where {N}
    quote
        $(Expr(:meta,:inline))
        Base.Cartesian.@ntuple $N n -> maybestaticlast(a[n])
    end
end
@inline maybestaticfirst(A::CartesianIndices) = CartesianVIndex(_maybestaticfirst(A.indices))
@inline maybestaticlast(A::CartesianIndices) = CartesianVIndex(_maybestaticlast(A.indices))

for (op,f) ∈ [(:(+),:vadd_nsw), (:(-),:vsub_nsw), (:(*),:vmul_nsw)]
  @eval begin
    @inline Base.$op(a::CartesianVIndex, b) = CartesianVIndex(fmap($f, getfield(a,:I), b))
    @inline Base.$op(a, b::CartesianVIndex) = CartesianVIndex(fmap($f, a, getfield(b,:I)))
    @inline Base.$op(a::CartesianVIndex, b::CartesianVIndex) = CartesianVIndex(fmap($f, getfield(a,:I), getfield(b,:I)))
  end
end


text:
import Base: *, ∈, isempty

export LinearMap,
       an_element,
       constraints_list

"""
    LinearMap{N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}} <: LazySet{N}

Type that represents a linear transformation ``M⋅S`` of a convex set ``S``.

### Fields

- `M` -- matrix/linear map
- `X` -- convex set

### Notes

This type is parametric in the elements of the linear map, `NM`, which is
independent of the numeric type of the wrapped set (`N`).
Typically `NM = N`, but there may be exceptions, e.g., if `NM` is an interval
that holds numbers of type `N`, where `N` is a floating point number type such
as `Float64`.
"""
struct LinearMap{N<:Real, S<:LazySet{N},
                 NM, MAT<:AbstractMatrix{NM}} <: LazySet{N}
    M::MAT
    X::S

    # default constructor with dimension match check
    function LinearMap{N, S, NM, MAT}(M::MAT, X::S) where {N<:Real,
            S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}}
        @assert dim(X) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(X))"
        return new{N, S, NM, MAT}(M, X)
    end
end

# convenience constructor without type parameter
LinearMap(M::MAT,
          X::S) where {N<:Real, S<:LazySet{N}, NM, MAT<:AbstractMatrix{NM}} =
    LinearMap{N, S, NM, MAT}(M, X)

# constructor from a linear map: perform the matrix multiplication immediately
LinearMap(M::MAT, lm::LinearMap{N, S, NM, MAT}) where {N<:Real, S<:LazySet{N},
        NM, MAT<:AbstractMatrix{NM}} =
    LinearMap{N, S, NM, MAT}(M * lm.M, lm.X)

"""
```
    *(M::AbstractMatrix{N}, X::LazySet{N}) where {N<:Real}
```

Return the linear map of a convex set.

### Input

- `M` -- matrix/linear map
- `X` -- convex set

### Output

A lazy linear map, i.e. a `LinearMap` instance.
"""
function *(M::AbstractMatrix{N}, X::LazySet{N}) where {N<:Real}
    return LinearMap(M, X)
end

function *(M::AbstractVector{N}, X::LazySet{N}) where {N<:Real}
    return LinearMap(reshape(M, length(M), 1), X)
end

"""
```
    *(α::N, X::LazySet{N}) where {N<:Real}
```

Return the linear map of a convex set scaled by a given number.

### Input

- `α` -- scalar
- `X` -- convex set

### Output

The lazy linear map of the convex set, `X ↦ M * X`, such that `M` is a sparse
matrix proportional to the identity and whose non-zero entries equal `α`.
"""
function *(α::N, X::LazySet{N}) where {N<:Real}
    n = dim(X)
    return LinearMap(sparse(α * I, n, n), X)
end

"""
```
    *(α::N, lm::LinearMap{N}) where {N<:Real, LM<:LinearMap{N}}
```

Return the linear map scaled by a given value.

### Input

- `α`  -- scalar
- `lm` -- linear map

### Output

The scaled linear map.
"""
function *(α::N, lm::LinearMap{N}) where {N<:Real}
    return LinearMap(α * lm.M, lm.X)
end

"""
```
    *(M::AbstractMatrix{N}, Z::ZeroSet{N})::ZeroSet{N} where {N<:Real}
```

A linear map of a zero set, which is simplified to a zero set (the absorbing
element).

### Input

- `M` -- abstract matrix
- `Z` -- zero set

### Output

The zero set with the output dimension of the linear map.
"""
function *(M::AbstractMatrix{N}, Z::ZeroSet{N})::ZeroSet{N} where {N<:Real}
    @assert dim(Z) == size(M, 2) "a linear map of size $(size(M)) cannot " *
            "be applied to a set of dimension $(dim(Z))"
    return ZeroSet{N}(size(M, 1))
end

"""
    dim(lm::LinearMap)::Int

Return the dimension of a linear map.

### Input

- `lm` -- linear map

### Output

The ambient dimension of the linear map.
"""
function dim(lm::LinearMap)::Int
    return size(lm.M, 1)
end

"""
    σ(d::AbstractVector{N}, lm::LinearMap{N}) where {N<:Real}

Return the support vector of the linear map.

### Input

- `d`  -- direction
- `lm` -- linear map

### Output

The support vector in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it follows
that ``σ(d, L) = M⋅σ(M^T d, S)`` for any direction ``d``.
"""
function σ(d::AbstractVector{N}, lm::LinearMap{N}) where {N<:Real}
    return lm.M * σ(_At_mul_B(lm.M, d), lm.X)
end

"""
    ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real}

Return the support function of the linear map.

### Input

- `d`      -- direction
- `lm`     -- linear map
- `kwargs` -- additional arguments that are passed to the support function
              algorithm

### Output

The support function in the given direction.
If the direction has norm zero, the result depends on the wrapped set.

### Notes

If ``L = M⋅S``, where ``M`` is a matrix and ``S`` is a convex set, it follows
that ``ρ(d, L) = ρ(M^T d, S)`` for any direction ``d``.
"""
function ρ(d::AbstractVector{N}, lm::LinearMap{N}; kwargs...) where {N<:Real}
    return ρ(_At_mul_B(lm.M, d), lm.X; kwargs...)
end

"""
    isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool

Determine whether a linear map is bounded.

### Input

- `lm`       -- linear map
- `cond_tol` -- (optional) tolerance of matrix condition (used to check whether
                the matrix is invertible)

### Output

`true` iff the linear map is bounded.

### Algorithm

We first check if the matrix is zero or the wrapped set is bounded.
If not, we perform a sufficient check whether the matrix is invertible.
If the matrix is invertible, then the map being bounded is equivalent to the
wrapped set being bounded, and hence the map is unbounded.
Otherwise, we check boundedness via [`isbounded_unit_dimensions`](@ref).
"""
function isbounded(lm::LinearMap; cond_tol::Number=DEFAULT_COND_TOL)::Bool
    if iszero(lm.M) || isbounded(lm.X)
        return true
    end
    if isinvertible(lm.M; cond_tol=cond_tol)
        return false
    end
    return isbounded_unit_dimensions(lm)
end

"""
    ∈(x::AbstractVector{N}, lm::LinearMap{N})::Bool where {N<:Real}

Check whether a given point is contained in a linear map of a convex set.

### Input

- `x`  -- point/vector
- `lm` -- linear map of a convex set

### Output

`true` iff ``x ∈ lm``.

### Algorithm

Note that ``x ∈ M⋅S`` iff ``M^{-1}⋅x ∈ S``.
This implementation does not explicitly invert the matrix, which is why it also
works for non-square matrices.

### Examples

```jldoctest
julia> lm = LinearMap([2.0 0.0; 0.0 1.0], BallInf([1., 1.], 1.));

julia> ∈([5.0, 1.0], lm)
false
julia> ∈([3.0, 1.0], lm)
true
```

An example with non-square matrix:
```jldoctest
julia> B = BallInf(zeros(4), 1.);

julia> M = [1. 0 0 0; 0 1 0 0]/2;

julia> ∈([0.5, 0.5], M*B)
true
```
"""
function ∈(x::AbstractVector{N}, lm::LinearMap{N})::Bool where {N<:Real}
    return ∈(lm.M \ x, lm.X)
end

"""
    an_element(lm::LinearMap{N})::Vector{N} where {N<:Real}

Return some element of a linear map.

### Input

- `lm` -- linear map

### Output

An element in the linear map.
It relies on the `an_element` function of the wrapped set.
"""
function an_element(lm::LinearMap{N})::Vector{N} where {N<:Real}
    return lm.M * an_element(lm.X)
end

"""
    isempty(lm::LinearMap)::Bool

Return if a linear map is empty or not.

### Input

- `lm` -- linear map

### Output

`true` iff the wrapped set is empty.
"""
function isempty(lm::LinearMap)::Bool
    return isempty(lm.X)
end

"""
    vertices_list(lm::LinearMap{N})::Vector{Vector{N}} where {N<:Real}

Return the list of vertices of a (polyhedral) linear map.

### Input

- `lm` -- linear map

### Output

A list of vertices.

### Algorithm

We assume that the underlying set `X` is polyhedral.
Then the result is just the linear map applied to the vertices of `X`.
"""
function vertices_list(lm::LinearMap{N})::Vector{Vector{N}} where {N<:Real}
    # for a zero map, the result is just the list containing the origin
    if iszero(lm.M)
        return [zeros(N, dim(lm))]
    end

    # collect low-dimensional vertices lists
    vlist_X = vertices_list(lm.X)

    # create resulting vertices list
    vlist = Vector{Vector{N}}()
    sizehint!(vlist, length(vlist_X))
    for v in vlist_X
        push!(vlist, lm.M * v)
    end

    return vlist
end

"""
    constraints_list(lm::LinearMap{N}) where {N<:Real}

Return the list of constraints of a (polyhedral) linear map.

### Input

- `lm` -- linear map

### Output

The list of constraints of the linear map.

### Notes

We assume that the underlying set `X` is polyhedral, i.e., offers a method
`constraints_list(X)`.

### Algorithm

We fall back to a concrete set representation and apply `linear_map`.
"""
function constraints_list(lm::LinearMap{N}) where {N<:Real}
    return constraints_list(linear_map(lm.M, lm.X))
end

"""
    linear_map(M::AbstractMatrix{N}, lm::LinearMap{N}) where {N<:Real}

Return the linear map of a lazy linear map.

### Input

- `M`  -- matrix
- `lm` -- linear map

### Output

The polytope representing the linear map of the lazy linear map of a set.
"""
function linear_map(M::AbstractMatrix{N}, lm::LinearMap{N}) where {N<:Real}
     return linear_map(M * lm.M, lm.X)
end

ext: jl


text:
struct ExploitabilityHistory
    x::Vector{Int}
    y::Vector{Float64}
    ExploitabilityHistory() = new(Int[], Float64[])
end

function Base.push!(h::ExploitabilityHistory, x, y)
    push!(h.x, x)
    push!(h.y, y)
end

mutable struct ExploitabilityCallback{SOL<:AbstractCFRSolver, ESOL<:ExploitabilitySolver}
    sol::SOL
    e_sol::ESOL
    n::Int
    state::Int
    hist::ExploitabilityHistory
end

function ExploitabilityCallback(sol::AbstractCFRSolver, n::Int=1; p::Int=1)
    e_sol = ExploitabilitySolver(sol, p)
    return ExploitabilityCallback(sol, e_sol, n, 0, ExploitabilityHistory())
end

function (cb::ExploitabilityCallback)()
    if iszero(rem(cb.state, cb.n))
        e = exploitability(cb.e_sol, cb.sol)
        push!(cb.hist, cb.state, e)
    end
    cb.state += 1
end

@recipe function f(hist::ExploitabilityHistory)
    xlabel := "Training Steps"
    @series begin
        subplot := 1
        ylabel := "Exploitability"
        label := ""
        hist.x, hist.y
    end
end

@recipe f(cb::ExploitabilityCallback) = cb.hist

mutable struct Throttle{F}
    f::F
    n::Int
    state::Int
end

function Throttle(f::Function, n::Int)
    return Throttle(f, n, 0)
end

function (t::Throttle)()
    iszero(rem(t.state, t.n)) && t.f()
    t.state += 1
end

struct CallbackChain{T<:Tuple}
	t::T
	CallbackChain(args...) = new{typeof(args)}(args)
end

Base.iterate(chain::CallbackChain, s=1) = iterate(chain.t, s)

function (chain::CallbackChain)()
	for cb in chain
		cb()
	end
end

max_stars_repo_name: WhiffleFish/HelloCFR

max_issues_repo_name: WhiffleFish/HelloCFR

max_forks_repo_name: WhiffleFish/HelloCFR


text:
## Contributed by @doublestrong

## Need 6 npz file and timing

using Caesar
# using RoME, Distributions
# using GraphPlot
using RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

# using Rebugger
# using Gadfly
using NPZ
using MAT

# # add multiproc support
using Distributed
addprocs(3)
using Caesar
@everywhere using Caesar, RoME

# # 90% precompile
for i in 1:5
  warmUpSolverJIT()
end

##

# # [OPTIONAL] add more julia processes to speed up inference
# using Distributed
# nprocs() < 3 ? addprocs(4-nprocs()) : nothing

# @everywhere using RoME
N = 1000
# start with an empty factor graph object
fg = initfg()

# fieldnames(typeof(getSolverParams(fg)))
# getSolverParams(fg).graphinit = true  # init on graph first then solve on tree (default)
# getSolverParams(fg).graphinit = false # init and solve on tree
getSolverParams(fg).useMsgLikelihoods = true
getSolverParams(fg).N = 150

# getSolverParams(fg).N = N
# # Add the first pose :x0
# addVariable!(fg, :x0, Pose2)
# # Add a few more poses
# for i in 1:5
#   addVariable!(fg, Symbol("x$(i)"), Pose2)
# end
# # Add landmarks
# for i in 1:2
#   addVariable!(fg, Symbol("l$(i)"), Point2)
# end

timing = zeros(0)

## Step 0

# Add at a fixed location Prior to pin :x0 to a starting location (0,0,0)
s_time = time_ns()
addVariable!(fg, :x0, Pose2)
addVariable!(fg, :l1, Point2)
addFactor!(fg, [:x0], PriorPose2( MvNormal([0; 0; 0], diagm([0.3;0.3;0.3].^2)) ))
ppr = Pose2Point2Range(MvNormal([7.323699045815659], diagm([0.3].^2)))
addFactor!(fg, [:x0; :l1], ppr )
# Gadfly.draw(PDF("fg.pdf", 10cm, 10cm),pl)  # or PNG(...)
# ensureAllInitialized!(fg)
# dfgplot(fg)

## Perform inference
tree, smt, hist = solveTree!(fg)
e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

# pl = plotKDE(fg, [:x0,:l1],dims=[1;2], levels=4)
pl = plotKDE(fg, [:x0,:l1],dims=[1;2], levels=4)

##

Gadfly.draw(PDF("batch_1.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0) ,N)
L1 = rand(getBelief(fg, :l1) ,N)
batch1 = vcat(X0,L1)
npzwrite("batch1.npz",batch1)

# plX1 |> PNG("/tmp/testX1.png")
# X2 = rand(getBelief(fg, :x2) ,N)
# X3 = rand(getBelief(fg, :x3) ,N)
# X4 = rand(getBelief(fg, :x4) ,N)
# X5 = rand(getBelief(fg, :x5) ,N)
# res = vcat(X0,X1,L1)
# using NPZ
# npzwrite("res.npz",res)
# Gadfly.draw(PDF("batch.pdf", 10cm, 10cm),pl)  # or PNG(...)
# error("here's a breakpoint")

## Step 1
i = 1
s_time = time_ns()
addVariable!(fg, :x1, Pose2)
pp = Pose2Pose2(MvNormal([9.82039106655709;-0.040893643507463134;0.7851856362659602], diagm([0.3;0.3;0.05].^2)))
ppr = Pose2Point2Range(MvNormal([6.783355788017825], diagm([0.3].^2)))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); :l1], ppr )

##

tree, smt, hist = solveTree!(fg, tree)

e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

pl = plotKDE(fg, [:x0,:x1,:l1],dims=[1;2], levels=4)

##

Gadfly.draw(PDF("batch_2.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0) ,N)
L1 = rand(getBelief(fg, :l1) ,N)
X1 = rand(getBelief(fg, :x1) ,N)
batch2 = vcat(X0,X1,L1)
npzwrite("batch2.npz",batch2)

# error("here's a breakpoint")
## Step 2

i = i + 1
s_time = time_ns()
addVariable!(fg, :x2, Pose2)
addVariable!(fg, :l2, Point2)
pp = Pose2Pose2(MvNormal([9.824301116341609;-0.1827443713759503;0.7586281983933238], diagm([0.3;0.3;0.05].^2)))
ppr = Pose2Point2Range(MvNormal([6.768864478545091], diagm([0.3].^2)))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); :l2], ppr )

##

tree, smt, hist = solveTree!(fg, tree)
# tree, smt, hist = solveTree!(fg; smtasks, recordcliqs=ls(fg))

e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

pl = plotKDE(fg, [:x0,:x1,:x2,:l1,:l2],dims=[1;2], levels=4)

##
Gadfly.draw(PDF("batch_3.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0) ,N)
L1 = rand(getBelief(fg, :l1) ,N)
L2 = rand(getBelief(fg, :l2) ,N)
X1 = rand(getBelief(fg, :x1) ,N)
X2 = rand(getBelief(fg, :x2) ,N)
batch3 = vcat(X0,X1,X2,L1,L2)
npzwrite("batch3.npz",batch3)

#
## Step 3

i = i + 1
s_time = time_ns()
addVariable!(fg, :x3, Pose2)
pp = Pose2Pose2(MvNormal([9.776502885351334; -0.010587078502017132; 1.5591793408311467], diagm([0.3;0.3;0.05].^2)))
ppr = Pose2Point2Range(MvNormal([7.401417053438512], diagm([0.3].^2)))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); :l2], ppr )

##
tree, smt, hist = solveTree!(fg, tree)
e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

pl = plotKDE(fg, [:x0,:x1,:x2,:x3,:l1,:l2],dims=[1;2], levels=4)

##
Gadfly.draw(PDF("batch_4.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0) ,N)
L1 = rand(getBelief(fg, :l1) ,N)
L2 = rand(getBelief(fg, :l2) ,N)
X1 = rand(getBelief(fg, :x1) ,N)
X2 = rand(getBelief(fg, :x2) ,N)
X3 = rand(getBelief(fg, :x3) ,N)
batch4 = vcat(X0,X1,X2,X3,L1,L2)
npzwrite("batch4.npz",batch4)

## Step 4 (:l2 has triple range measurements)
i = i + 1
s_time = time_ns()
addVariable!(fg, Symbol("x4"), Pose2)
pp = Pose2Pose2(MvNormal([9.644657137571507; 0.5847494762836476; 0.7422440549101994], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([7.331883435735532], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l2")], ppr )

##
tree, smt, hist = solveTree!(fg, tree)
e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

pl = plotKDE(fg, [:x0,:x1,:x2,:x3,:x4,:l1,:l2],dims=[1;2], levels=4)

##
Gadfly.draw(PDF("batch_5.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0) ,N)
L1 = rand(getBelief(fg, :l1) ,N)
L2 = rand(getBelief(fg, :l2) ,N)
X1 = rand(getBelief(fg, :x1),N)
X2 = rand(getBelief(fg, :x2) ,N)
X3 = rand(getBelief(fg, :x3) ,N)
X4 = rand(getBelief(fg, :x4),N)
batch5 = vcat(X0,X1,X2,X3,X4,L1,L2)
npzwrite("batch5.npz",batch5)

## Step 5 (triple range to :l1)
i = i + 1
s_time = time_ns()
addVariable!(fg, Symbol("x5"), Pose2)
pp = Pose2Pose2(MvNormal([9.725096752125593; 0.6094800276622434; 0.7750400402527422], Matrix(Diagonal([0.3;0.3;0.05].^2))))
ppr = Pose2Point2Range(MvNormal([6.476782344345081], Matrix(Diagonal([0.3].^2))))
addFactor!(fg, [Symbol("x$(i-1)"); Symbol("x$(i)")], pp )
addFactor!(fg, [Symbol("x$(i)"); Symbol("l1")], ppr )

##
tree, smt, hist = solveTree!(fg, tree)
e_time = time_ns()
elapsed = (e_time - s_time)/10^9
append!(timing, elapsed)
println("elapsed time: ", elapsed)

pl = plotKDE(fg, [:x0,:x1,:x2,:x3,:x4,:x5,:l1,:l2],dims=[1;2], levels=4)

##
Gadfly.draw(PDF("batch_6.pdf", 20cm, 10cm),pl)
X0 = rand(getBelief(fg, :x0),N)
L1 = rand(getBelief(fg, :l1),N)
L2 = rand(getBelief(fg, :l2),N)
X1 = rand(getBelief(fg, :x1),N)
X2 = rand(getBelief(fg, :x2),N)
X3 = rand(getBelief(fg, :x3),N)
X4 = rand(getBelief(fg, :x4),N)
X5 = rand(getBelief(fg, :x5),N)
batch6 = vcat(X0,X1,X2,X3,X4,X5,L1,L2)
npzwrite("batch6.npz",batch6)
npzwrite("timing.npz",timing)

#
#
# pl = plotSLAM2D(fg)
#
# # For scripting use-cases you can export the image
# Gadfly.draw(PDF("/tmp/test.pdf", 20cm, 10cm),pl)  # or PNG(...)
#
# L1 = rand(getBelief(fg, :l1) ,N)
# # L2 = rand(getBelief(fg, :l2) ,N)
# X0 = rand(getBelief(fg, :x0) ,N)
# X1 = rand(getBelief(fg, :x1) ,N)
# # X2 = rand(getBelief(fg, :x2) ,N)
# # X3 = rand(getBelief(fg, :x3) ,N)
# # X4 = rand(getBelief(fg, :x4) ,N)
# # X5 = rand(getBelief(fg, :x5) ,N)
# res = vcat(X0,X1,L1)
# using NPZ
# npzwrite("res.npz",res)
# # pl2 = Gadfly.plot(x=X0[1,:],y=X0[2,:], Geom.hexbin);
# # p1H = hstack(pl2)
# #
# # # convert to file
# # p1H |> PNG("/home/chad/Pictures/julia_x0_step0.png")
#
# # using RoMEPlotting
#
# #drawPoses(fg)
# # If you have landmarks, you can instead call
# # drawPosesLandms(fg)
# #
# # # # Draw the KDE for x0
# # # plotKDE(fg, :x0)
# # # # Draw the KDE's for x0 and x1
# # plotKDE(fg, [:x1])
# # plotKDE(fg, [:x2])
# # plotKDE(fg, [:x3])
# # plotKDE(fg, [:x4])
# # plotKDE(fg, [:x5])
# # plotKDE(fg, [:l1])
# # plotKDE(fg, [:l2])


text:
function launch()

    @static if Sys.isapple()
        # OpenGL 3.2 + GLSL 150
        glsl_version = 150
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 2)
        GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # required on Mac
    else
        # OpenGL 3.0 + GLSL 130
        glsl_version = 130
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MAJOR, 3)
        GLFW.WindowHint(GLFW.CONTEXT_VERSION_MINOR, 0)
        # GLFW.WindowHint(GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE) # 3.2+ only
        # GLFW.WindowHint(GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE) # 3.0+ only
    end

    # setup GLFW error callback
    error_callback(err::GLFW.GLFWError) = @error "GLFW ERROR: code $(err.code) msg: $(err.description)"
    GLFW.SetErrorCallback(error_callback)

    # create window
    window = GLFW.CreateWindow(1280, 720, "ElectroDermal Activity Analysis")
    @assert window != C_NULL
    GLFW.MakeContextCurrent(window)
    GLFW.SwapInterval(1)  # enable vsync

    # setup Dear ImGui context
    ctx = CImGui.CreateContext()

    # setup Dear ImGui style
    # CImGui.StyleColorsDark()
    # CImGui.StyleColorsClassic()
    CImGui.StyleColorsLight()

    # load Fonts
    # - If no fonts are loaded, dear imgui will use the default font. You can also load multiple fonts and use `CImGui.PushFont/PopFont` to select them.
    # - `CImGui.AddFontFromFileTTF` will return the `Ptr{ImFont}` so you can store it if you need to select the font among multiple.
    # - If the file cannot be loaded, the function will return C_NULL. Please handle those errors in your application (e.g. use an assertion, or display an error and quit).
    # - The fonts will be rasterized at a given size (w/ oversampling) and stored into a texture when calling `CImGui.Build()`/`GetTexDataAsXXXX()``, which `ImGui_ImplXXXX_NewFrame` below will call.
    # - Read 'fonts/README.txt' for more instructions and details.
    # fonts_dir = joinpath(@__DIR__, "..", "fonts")
    # fonts = CImGui.GetIO().Fonts
    # default_font = CImGui.AddFontDefault(fonts)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Cousine-Regular.ttf"), 15)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "DroidSans.ttf"), 16)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Karla-Regular.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "ProggyTiny.ttf"), 10)
    # CImGui.AddFontFromFileTTF(fonts, joinpath(fonts_dir, "Roboto-Medium.ttf"), 16)
    # @assert default_font != C_NULL

    # setup Platform/Renderer bindings
    ImGui_ImplGlfw_InitForOpenGL(window, true)
    ImGui_ImplOpenGL3_Init(glsl_version)
    clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]

    open_file_dialog = OpenFileDialog(pwd(),"", pwd(),"", false, false)
    open_file_dialog2 = OpenFileDialog(pwd(),"",pwd(),"",false,false)
    data = DataFrame()
    dt = Float32[]
    data2 = DataFrame()
    dt2 = Float32[]
    Open_files = false
    Open_files2 = false
    data_length = 0
    #eda_data = nothing
    #isnothing(eda_data)
    while !GLFW.WindowShouldClose(window)
        GLFW.PollEvents()
        # start the Dear ImGui frame
        ImGui_ImplOpenGL3_NewFrame()
        ImGui_ImplGlfw_NewFrame()
        CImGui.NewFrame()

        if CImGui.BeginMainMenuBar()
            if CImGui.BeginMenu("File")
                populate_file_menu!(open_file_dialog,open_file_dialog2) #visiblity = true
                #populate_file_menu!(open_file_dialog2)
                CImGui.EndMenu()
            end
            CImGui.EndMainMenuBar()
        end
        if isvisible(open_file_dialog)
            display_dialog!(open_file_dialog)
            if has_pending_action(open_file_dialog)
                # Check filename and perform different action.
                data = perform_dialog_action(open_file_dialog)
                # eda_data = ...
                # hr_data = ...
                consume_action!(open_file_dialog)
                #dt = Cfloat.(data[3:end,1])
                Base.display(data)
                Open_files = true
            end
        end

        if isvisible(open_file_dialog2)
            display_dialog!(open_file_dialog2)
            if has_pending_action(open_file_dialog2)
                data2 = perform_dialog_action(open_file_dialog2)
                consume_action!(open_file_dialog2)
                #dt = Cfloat.(data[3:end,1])
                Base.display(data2)
                Open_files2 = true
            end
        end
            #@c CImGui.Checkbox("Open files", &Open_files)
        if Open_files
            if !Open_files2
                # data_length = ...
                plotData(data)
            elseif Open_files2
                plotData2(data,data2)
            end
            CImGui.End()
        end

        # rendering
        CImGui.Render()
        GLFW.MakeContextCurrent(window)
        display_w, display_h = GLFW.GetFramebufferSize(window)
        glViewport(0, 0, display_w, display_h)
        glClearColor(clear_color...)
        glClear(GL_COLOR_BUFFER_BIT)
        ImGui_ImplOpenGL3_RenderDrawData(CImGui.GetDrawData())

        GLFW.MakeContextCurrent(window)
        GLFW.SwapBuffers(window)
    end

    # cleanup
    ImGui_ImplOpenGL3_Shutdown()
    ImGui_ImplGlfw_Shutdown()
    CImGui.DestroyContext(ctx)

    GLFW.DestroyWindow(window)
end

function perform_dialog_action(dialog::OpenFileDialog)
    directory = get_directory(dialog, ConfirmedStatus())
    file_name = get_file(dialog, ConfirmedStatus())
    path = joinpath(directory, file_name)
    data = CSV.File(path ; header= ["EDA"]) |> DataFrame
end

function populate_file_menu!(dialog::AbstractDialog,dialog2::AbstractDialog)
    if CImGui.MenuItem("Open", "Ctrl+O")
        set_visibility!(dialog, true)
    end
    if CImGui.MenuItem("Open2", "Ctrl+O")
        set_visibility!(dialog2, true)
    end
    if CImGui.MenuItem("Save", "Ctrl+s")

    end
    if CImGui.MenuItem("Save as")

    end
    if CImGui.MenuItem("Quit", "Alt+F4")

    end
end

function plotData(data::DataFrame)
    begin
        #df = CSV.read("F:\\julia\\CSV\\EDA.csv", header = ["EDA"])
        #df = CSV.read("C:\\Users\\msi-\\Desktop\\EDA\\EDA.jl\\src\\HR.csv", header = ["HR"])
        CImGui.Begin("Plot")
        CImGui.Text("Plot")

        df_data =sort(Cfloat.(data[3:end,1]))
        max_value = df_data[end,1]
        #max = Cfloat(max_value)
        #length = size(df,1)-2
        len = size(data,1)-2
        #@show typeof(len)

        st = data[1,1]
        freq = data[2,1]
        start_time, end_time = @cstatic start_time=Cint(1) end_time=Cint(2) begin
            @c CImGui.SliderInt("Start Time", &start_time, 1,len)
            CImGui.SameLine(0.0,CImGui.GetStyle().ItemInnerSpacing.x)
            time1=Dates.unix2datetime(st+(start_time-1)/freq)
            CImGui.Text(string(Time(time1)))
            @c CImGui.SliderInt("End Time", &end_time, 2,len)
            CImGui.SameLine(0.0,CImGui.GetStyle().ItemInnerSpacing.x)
            time2=Dates.unix2datetime(st+(end_time-1)/freq)
            CImGui.Text(string(Time(time2)))
            if start_time > end_time -1
                start_time = end_time - Cint(1)
            elseif start_time < end_time - 2000*freq
                start_time = end_time - Cint(2000*freq)
            end
        end

        plot_data = Cfloat.(data[(start_time+2):(end_time+2),1])
        # st = data[1,1]
        # freq = data[2,1]

        CImGui.Text("Primitives")
            sz, thickness, col = @cstatic sz=Cfloat(36.0) thickness=Cfloat(4.0) col=Cfloat[1.0,0.0,0.4,0.2] begin
                CImGui.ColorEdit4("Color", col)
            end
            p = CImGui.GetCursorScreenPos()
            col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))

            begin
                width = 1200
                height = 200

                CImGui.PlotLines("EDA Measurements", plot_data, length(plot_data), 0 , "EDA", 0, Cfloat(max_value*1.2), (width,height))
                draw_list = CImGui.GetWindowDrawList()
                x::Cfloat = p.x
                y::Cfloat = p.y
                spacing = 8.0
                    # Draws (almost transparent) horizontal bars
                for yₙ in range(y, step = 40, stop = y + height - 40)
                    CImGui.AddRectFilled(draw_list, ImVec2(x, yₙ), ImVec2(x+width, yₙ+20), col32);
                end
            end

            #draw the x axis
            p = CImGui.GetCursorScreenPos()
            begin
                width = 1200
                col = Cfloat[0.0,0.0,0.0,1.0]
                col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
                draw_list = CImGui.GetWindowDrawList()
                x = p.x
                y = p.y

                time = Dates.unix2datetime(st+ (start_time-1)/freq)
                f = Cfloat(1 / freq)
                CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+width, y), col32, Cfloat(1));
                for xₙ in range(x, step = 1200/ (end_time-start_time+1), stop = x + width)
                    hou = hour(time);
                    min = minute(time);
                    sec = second(time);
                    mil = millisecond(time);
                    if sec ==0 && mil ==0
                        CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y-5), col32, Cfloat(1));
                    end
                        #if mil ==0
                    if sec ==0 && mil ==0
                            #CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(min), ":", string(sec)));
                        CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(hou), ":", string(min)));
                    end
                    next_time = time + Dates.Millisecond(f*1000);
                    time = next_time;
                end
            end
    end
end

function plotData2(data::DataFrame,data2::DataFrame)
    begin
        CImGui.Begin("Plot")
        CImGui.Text("Plot")

        df_data =sort(Cfloat.(data[3:end,1]))
        df_data2 =sort(Cfloat.(data2[3:end,1]))
        max_value = df_data[end,1]
        max_value2 = df_data2[end,1]
        len = size(data,1)-2
        len2 = size(data2,1) -2
        st = data[1,1]
        freq = data[2,1]
        st2 = data2[1,1]
        freq2 = data2[2,1]

        start_time, end_time = @cstatic start_time=Cint(1) end_time=Cint(2) begin
            @c CImGui.SliderInt("Start Time", &start_time, 1,len)
            CImGui.SameLine(0.0,CImGui.GetStyle().ItemInnerSpacing.x)
            time1=Dates.unix2datetime(st+(start_time-1)/freq)
            CImGui.Text(string(Time(time1)))
            @c CImGui.SliderInt("End Time", &end_time, 2,len)
            CImGui.SameLine(0.0,CImGui.GetStyle().ItemInnerSpacing.x)
            time2=Dates.unix2datetime(st+(end_time-1)/freq)
            CImGui.Text(string(Time(time2)))
            if start_time > end_time -1
                start_time = end_time - Cint(1)
            elseif start_time < end_time - 2000*freq
                start_time = end_time - Cint(2000*freq)
            end
        end

        start_time2, end_time2 = @cstatic start_time2=Cint(1) end_time2=Cint(2) begin
            @c CImGui.SliderInt("Start Time2", &start_time2, 1,len2)
            CImGui.SameLine()
            time3=Dates.unix2datetime(st2+(start_time2-1)/freq2)
            CImGui.Text(string(Time(time3)))
            @c CImGui.SliderInt("End Time2", &end_time2, 2,len2)
            CImGui.SameLine()
            time4=Dates.unix2datetime(st2+(end_time2-1)/freq2)
            CImGui.Text(string(Time(time4)))
            if start_time2 > end_time2 -1
                start_time2 = end_time2 - Cint(1)
            elseif start_time2 < end_time2 - 2000*freq2
                start_time2 = end_time2 - Cint(2000*freq2)
            end
        end

        plot_data = Cfloat.(data[(start_time+2):(end_time+2),1])
        # st = data[1,1]
        # freq = data[2,1]

        plot_data2 = Cfloat.(data2[(start_time2+2):(end_time2+2),1])
        # st2 = data2[1,1]
        # freq2 = data2[2,1]

        CImGui.Text("Primitives")
            sz, thickness, col = @cstatic sz=Cfloat(36.0) thickness=Cfloat(4.0) col=Cfloat[1.0,0.0,0.4,0.2] begin
                CImGui.ColorEdit4("Color", col)
            end
            p = CImGui.GetCursorScreenPos()
            col32 = CImGui.ColorConvertFloat4ToU32(ImVec4(col...))
        begin
            width = 1200
            height = 200

            CImGui.PlotLines("EDA Measurements", plot_data, length(plot_data), 0 , "EDA", 0, Cfloat(max_value*1.2), (width,height))
            draw_list = CImGui.GetWindowDrawList()
            x::Cfloat = p.x
            y::Cfloat = p.y
            spacing = 8.0
            # Draws (almost transparent) horizontal bars
            for yₙ in range(y, step = 40, stop = y + height - 40)
                CImGui.AddRectFilled(draw_list, ImVec2(x, yₙ), ImVec2(x+width, yₙ+20), col32);
            end
        end

        #draw the x axis
        p = CImGui.GetCursorScreenPos()
        begin
            width = 1200
            col_ = Cfloat[0.0,0.0,0.0,1.0]
            col32_ = CImGui.ColorConvertFloat4ToU32(ImVec4(col_...))
            draw_list = CImGui.GetWindowDrawList()
            x = p.x
            y = p.y

            time = Dates.unix2datetime(st+ (start_time-1)/freq)
            f = Cfloat(1 / freq)
            CImGui.AddLine(draw_list, ImVec2(x, y), ImVec2(x+width, y), col32_, Cfloat(1));
            for xₙ in range(x, step = 1200/ (end_time-start_time+1), stop = x + width)
                hou = hour(time);
                min = minute(time);
                sec = second(time);
                mil = millisecond(time);
                if sec ==0 && mil ==0
                    CImGui.AddLine(draw_list, ImVec2(xₙ, y), ImVec2(xₙ, y-5), col32_, Cfloat(1));
                end
                    #if mil ==0
                if sec ==0 && mil ==0
                        #CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(min), ":", string(sec)));
                    CImGui.AddText(draw_list, ImVec2(xₙ, y), col32_, string(string(hou), ":", string(min)));
                end
                next_time = time + Dates.Millisecond(f*1000);
                time = next_time;
            end
        end

        CImGui.NewLine()
        p = CImGui.GetCursorScreenPos()
        begin
            width = 1200
            height = 200
            CImGui.PlotLines("HR Measurements", plot_data2, length(plot_data2), 0 , "HR", 0, Cfloat(max_value2*1.2), (width,height))
            draw_list2 = CImGui.GetWindowDrawList()
            x = p.x
            y = p.y
            spacing = 8.0
            # Draws (almost transparent) horizontal bars
            for yₙ in range(y, step = 40, stop = y + height - 40)
                CImGui.AddRectFilled(draw_list2, ImVec2(x, yₙ), ImVec2(x+width, yₙ+20), col32);
            end
        end

        #draw the x axis
        p = CImGui.GetCursorScreenPos()
        begin
            width = 1200
            col_ = Cfloat[0.0,0.0,0.0,1.0]
            col32_ = CImGui.ColorConvertFloat4ToU32(ImVec4(col_...))
            draw_list2 = CImGui.GetWindowDrawList()
            x = p.x
            y = p.y

            time = Dates.unix2datetime(st2+ (start_time2-1)/freq2)
            f = Cfloat(1 / freq2)
            CImGui.AddLine(draw_list2, ImVec2(x, y), ImVec2(x+width, y), col32_, Cfloat(1));
            for xₙ in range(x, step = 1200/ (end_time2-start_time2+1), stop = x + width)
                hou = hour(time);
                min = minute(time);
                sec = second(time);
                mil = millisecond(time);
                if sec ==0 && mil ==0
                    CImGui.AddLine(draw_list2, ImVec2(xₙ, y), ImVec2(xₙ, y-5), col32_, Cfloat(1));
                end
                    #if mil ==0
                if sec ==0 && mil ==0
                        #CImGui.AddText(draw_list, ImVec2(xₙ, y), col32, string(string(min), ":", string(sec)));
                    CImGui.AddText(draw_list2, ImVec2(xₙ, y), col32_, string(string(hou), ":", string(min)));
                end
                next_time = time + Dates.Millisecond(f*1000);
                time = next_time;
            end
        end

    end
end


text:
# This file is a part of MJDSigGen, licensed under the MIT License (MIT).

struct CartPoint{T<:Real}
    x::T
    y::T
    z::T
end

CartPoint{T}() where {T} = CartPoint{T}(0, 0, 0)

struct CylPoint{T<:Real}
    r::T
    phi::T
    z::T
end

CylPoint{T}() where {T} = CylPoint{T}(0, 0, 0)

mutable struct VelocityLookup
    e::Cfloat
    e100::Cfloat
    e110::Cfloat
    e111::Cfloat
    h100::Cfloat
    h110::Cfloat
    h111::Cfloat
    ea::Cfloat
    eb::Cfloat
    ec::Cfloat
    ebp::Cfloat
    ecp::Cfloat
    ha::Cfloat
    hb::Cfloat
    hc::Cfloat
    hbp::Cfloat
    hcp::Cfloat
    hcorr::Cfloat
    ecorr::Cfloat

    VelocityLookup() = new()
end

mutable struct SigGenSetup
    # general
    verbosity::Cint

    # geometry
    xtal_length::Cfloat
    xtal_radius::Cfloat
    pc_length::Cfloat
    pc_radius::Cfloat
    wrap_around_radius::Cfloat
    ditch_depth::Cfloat
    ditch_thickness::Cfloat
    bottom_taper_length::Cfloat
    hole_length::Cfloat
    hole_radius::Cfloat
    outer_taper_length::Cfloat
    taper_angle::Cfloat
    inner_taper_length::Cfloat
    outer_taper_width::Cfloat
    inner_taper_width::Cfloat
    top_bullet_radius::Cfloat
    bottom_bullet_radius::Cfloat
    Li_thickness::Cfloat

    # electric fields & weighing potentials
    xtal_grid::Cfloat
    impurity_z0::Cfloat
    impurity_gradient::Cfloat
    impurity_quadratic::Cfloat
    impurity_surface::Cfloat
    impurity_radial_add::Cfloat
    impurity_radial_mult::Cfloat
    impurity_rpower::Cfloat
    xtal_HV::Cfloat
    max_iterations::Cint
    write_field::Cint
    write_WP::Cint
    bulletize_PC::Cint

    # file names
    config_name::NTuple{256,Cchar}
    drift_name::NTuple{256,Cchar}
    field_name::NTuple{256,Cchar}
    wp_name::NTuple{256,Cchar}

    # signal calculation
    xtal_temp::Cfloat
    preamp_tau::Cfloat
    time_steps_calc::Cint
    step_time_calc::Cfloat
    step_time_out::Cfloat
    charge_cloud_size::Cfloat
    use_acceleration::Cint
    use_repulsion::Cint
    use_diffusion::Cint
    energy::Cfloat

    coord_type::Cint
    ntsteps_out::Cint

    # data for fields.c
    rmin::Cfloat
    rmax::Cfloat
    rstep::Cfloat
    zmin::Cfloat
    zmax::Cfloat
    zstep::Cfloat
    rlen::Cint
    zlen::Cint
    v_lookup_len::Cint
    v_lookup::Ptr{VelocityLookup}
    efld::Ptr{Ptr{CylPoint{Cfloat}}}
    wpot::Ptr{Ptr{Cfloat}}

    # data for calc_signal.c
    dpath_e::Ptr{CartPoint{Cfloat}}
    dpath_h::Ptr{CartPoint{Cfloat}}
    instant_vel_e::Ptr{CartPoint{Cfloat}}
    instant_vel_h::Ptr{CartPoint{Cfloat}}
    instant_charge_size_e::Ptr{Cfloat}
    instant_charge_size_h::Ptr{Cfloat}
    initial_vel::Cfloat
    final_vel::Cfloat
    dv_dE::Cfloat
    v_over_E::Cfloat
    final_charge_size::Cdouble

    function SigGenSetup()
        s = new()
        fill!(unsafe_wrap(Array, Ptr{UInt8}(pointer_from_objref(s)), sizeof(s)), 0)
        return s
    end
end

SigGenSetup(config_filename::AbstractString) = signal_calc_init(config_filename)

function tup2str(tup::NTuple{N,C}) where {N, C<:Union{AbstractChar, Cchar}}
    out = Vector{UInt8}(undef, N)

    for i in 1:N
        c = @inbounds tup[i]
        c == 0 && (resize!(out, i-1); break)
        @inbounds out[i] = tup[i]
    end

    return String(out)
end

function Base.getproperty(setup::SigGenSetup, sym::Symbol)
    if sym in (:config_name, :drift_name, :field_name, :wp_name)
        return tup2str(getfield(setup, sym))
    elseif sym in (:dpath_e, :dpath_h, :instant_vel_e, :instant_vel_h)
        ptr = getfield(setup, sym)
        N = setup.time_steps_calc
        v = unsafe_wrap(Vector{CartPoint{Cfloat}}, ptr, N)
        return transpose(reshape(reinterpret(Cfloat, v), (3, N)))
    elseif sym in (:instant_charge_size_e, :instant_charge_size_h)
        ptr = getfield(setup, sym)
        return unsafe_wrap(Vector{Cfloat}, ptr, setup.time_steps_calc)
    else
        return getfield(setup, sym)
    end
end


text:
using DCELs
using Documenter

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

makedocs(;
    modules=[DCELs],
    authors="Rodolfo Figueroa Soriano <4rodolfofigueroa2@gmail.com>",
    repo="https://github.com/RodolfoFigueroa/DCELs.jl/blob/{commit}{path}#{line}",
    sitename="DCELs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://RodolfoFigueroa.github.io/DCELs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/RodolfoFigueroa/DCELs.jl",
    devbranch="main",
)


text:
module Qubism

using LinearAlgebra

include("register.jl")
include("hamiltonian.jl")
include("gates.jl")
include("helpers.jl")

# Unicode kron operator $\otimes$
const ⊗ = kron

# Export to user
export ⊗

end


text:
# Copyright (c) 2017: Miles Lubin and contributors
# Copyright (c) 2017: Google Inc.
#
# Use of this source code is governed by an MIT-style license that can be found
# in the LICENSE.md file or at https://opensource.org/licenses/MIT.

_lower_set(set::MOI.Interval) = MOI.GreaterThan(set.lower)
function _lower_set(set::MOI.Interval{T}) where {T<:AbstractFloat}
    if set.lower == typemin(T)
        return nothing
    else
        return MOI.GreaterThan(set.lower)
    end
end
_upper_set(set::MOI.Interval) = MOI.LessThan(set.upper)
function _upper_set(set::MOI.Interval{T}) where {T<:AbstractFloat}
    if set.upper == typemax(T)
        return nothing
    else
        return MOI.LessThan(set.upper)
    end
end
_lower_set(set::MOI.EqualTo) = MOI.GreaterThan(set.value)
_upper_set(set::MOI.EqualTo) = MOI.LessThan(set.value)
_lower_set(set::MOI.Zeros) = MOI.Nonnegatives(set.dimension)
_upper_set(set::MOI.Zeros) = MOI.Nonpositives(set.dimension)

"""
    SplitIntervalBridge{T, F, S, LS, US}

The `SplitIntervalBridge` splits a `F`-in-`S` constraint into a `F`-in-`LS` and
a `F`-in-`US` constraint where we have either:
* `S = MOI.Interval{T}`, `LS = MOI.GreaterThan{T}` and `US = MOI.LessThan{T}`,
* `S = MOI.EqualTo{T}`, `LS = MOI.GreaterThan{T}` and `US = MOI.LessThan{T}`, or
* `S = MOI.Zeros`, `LS = MOI.Nonnegatives` and `US = MOI.Nonpositives`.

For instance, if `F` is `MOI.ScalarAffineFunction` and `S` is `MOI.Interval`,
it transforms the constraint ``l ≤ ⟨a, x⟩ + α ≤ u`` into the constraints
``⟨a, x⟩ + α ≥ l`` and ``⟨a, x⟩ + α ≤ u``.

!!! note
    If `T<:AbstractFloat` and `S` is `MOI.Interval{T}` then no lower (resp.
    upper) bound constraint is created if the lower (resp. upper) bound is
    `typemin(T)` (resp. `typemax(T)`). Similarly, when
    [`MathOptInterface.ConstraintSet`](@ref) is set, a lower or upper bound
    constraint may be deleted or created accordingly.
"""
mutable struct SplitIntervalBridge{
    T,
    F<:MOI.AbstractFunction,
    S<:MOI.AbstractSet,
    LS<:MOI.AbstractSet,
    US<:MOI.AbstractSet,
} <: AbstractBridge
    lower::Union{Nothing,MOI.ConstraintIndex{F,LS}}
    upper::Union{Nothing,MOI.ConstraintIndex{F,US}}
    # To allow the user to do
    # ```jl
    # x = MOI.add_variable(model)
    # c = MOI.add_constraint(model, x, MOI.Interval(-Inf, Inf))
    # MOI.set(model, MOI.ConstraintSet(), c, MOI.Interval(0.0, Inf))
    # ```
    # we need to store the function to create the lower bound constraint.
    func::Union{Nothing,F}
end

const SplitInterval{T,OT<:MOI.ModelLike} =
    SingleBridgeOptimizer{SplitIntervalBridge{T},OT}

function bridge_constraint(
    ::Type{SplitIntervalBridge{T,F,S,LS,US}},
    model::MOI.ModelLike,
    f::F,
    set::S,
) where {T,F,S,LS,US}
    lower_set = _lower_set(set)
    if lower_set === nothing
        lower = nothing
    else
        lower = MOI.add_constraint(model, f, _lower_set(set))
    end
    upper_set = _upper_set(set)
    if upper_set === nothing
        upper = nothing
    else
        upper = MOI.add_constraint(model, f, _upper_set(set))
    end
    if lower === nothing && upper === nothing
        func = f
    else
        func = nothing
    end
    return SplitIntervalBridge{T,F,S,LS,US}(lower, upper, func)
end

function MOI.supports_constraint(
    ::Type{SplitIntervalBridge{T}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{<:Union{MOI.Interval{T},MOI.EqualTo{T}}},
) where {T}
    return true
end

function MOI.supports_constraint(
    ::Type{SplitIntervalBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.Zeros},
) where {T}
    return MOI.Utilities.is_coefficient_type(F, T)
end

function MOI.Bridges.added_constrained_variable_types(
    ::Type{<:SplitIntervalBridge},
)
    return Tuple{Type}[]
end

function MOI.Bridges.added_constraint_types(
    ::Type{SplitIntervalBridge{T,F,S,LS,US}},
) where {T,F,S,LS,US}
    return Tuple{Type,Type}[(F, LS), (F, US)]
end

function concrete_bridge_type(
    ::Type{<:SplitIntervalBridge},
    F::Type{<:MOI.AbstractScalarFunction},
    S::Type{<:Union{MOI.Interval{T},MOI.EqualTo{T}}},
) where {T}
    return SplitIntervalBridge{T,F,S,MOI.GreaterThan{T},MOI.LessThan{T}}
end

function concrete_bridge_type(
    ::Type{<:SplitIntervalBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{MOI.Zeros},
) where {T}
    return SplitIntervalBridge{T,F,MOI.Zeros,MOI.Nonnegatives,MOI.Nonpositives}
end

function MOI.get(
    bridge::SplitIntervalBridge{T,F,S,LS},
    ::MOI.NumberOfConstraints{F,LS},
)::Int64 where {T,F,S,LS}
    if bridge.lower === nothing
        return 0
    else
        return 1
    end
end

function MOI.get(
    bridge::SplitIntervalBridge{T,F,S,LS,US},
    ::MOI.NumberOfConstraints{F,US},
)::Int64 where {T,F,S,LS,US}
    if bridge.upper === nothing
        return 0
    else
        return 1
    end
end

function MOI.get(
    bridge::SplitIntervalBridge{T,F,S,LS},
    ::MOI.ListOfConstraintIndices{F,LS},
) where {T,F,S,LS}
    if bridge.lower === nothing
        return MOI.ConstraintIndex{F,LS}[]
    else
        return [bridge.lower]
    end
end

function MOI.get(
    bridge::SplitIntervalBridge{T,F,S,LS,US},
    ::MOI.ListOfConstraintIndices{F,US},
) where {T,F,S,LS,US}
    if bridge.upper === nothing
        return MOI.ConstraintIndex{F,US}[]
    else
        return [bridge.upper]
    end
end

function MOI.delete(model::MOI.ModelLike, bridge::SplitIntervalBridge)
    if bridge.lower !== nothing
        MOI.delete(model, bridge.lower)
    end
    if bridge.upper !== nothing
        MOI.delete(model, bridge.upper)
    end
    return
end

function MOI.supports(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimalStart,MOI.ConstraintDualStart},
    ::Type{SplitIntervalBridge{T,F,S,LS,US}},
) where {T,F,S,LS,US}
    ci_1 = MOI.ConstraintIndex{F,LS}
    ci_2 = MOI.ConstraintIndex{F,US}
    return MOI.supports(model, attr, ci_1) && MOI.supports(model, attr, ci_2)
end

function _error_double_inf(attr)
    return error(
        "Cannot get `$attr` for a constraint in the interval `[-Inf, Inf]`.",
    )
end

function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintPrimalStart},
    bridge::SplitIntervalBridge,
)
    # lower and upper should give the same value
    if bridge.lower !== nothing
        return MOI.get(model, attr, bridge.lower)
    end
    if bridge.upper !== nothing
        return MOI.get(model, attr, bridge.upper)
    end
    return _error_double_inf(attr)
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintPrimalStart,
    bridge::SplitIntervalBridge,
    value,
)
    if bridge.lower !== nothing
        MOI.set(model, attr, bridge.lower, value)
    end
    if bridge.upper !== nothing
        MOI.set(model, attr, bridge.upper, value)
    end
    return
end

# The map is:
# x ∈ S <=> [1 1]' * x ∈ LS × US
# So the adjoint map is
# [1 1] * y ∈ S* <=> y ∈ (LS × US)*
# where [1 1] * y = y[1] + y[2]
# so we can just sum the dual values.
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintDual,MOI.ConstraintDualStart},
    bridge::SplitIntervalBridge{T},
) where {T}
    if bridge.lower === nothing
        if bridge.upper === nothing
            return zero(T)
        else
            return MOI.get(model, attr, bridge.upper)
        end
    else
        if bridge.upper === nothing
            return MOI.get(model, attr, bridge.lower)
        else
            return MOI.get(model, attr, bridge.lower) +
                   MOI.get(model, attr, bridge.upper)
        end
    end
end

function _split_dual_start(value)
    if value < 0
        return zero(value), value
    else
        return value, zero(value)
    end
end

function _split_dual_start(value::Vector)
    lower = similar(value)
    upper = similar(value)
    for i in eachindex(value)
        lower[i], upper[i] = _split_dual_start(value[i])
    end
end

function MOI.set(
    model::MOI.ModelLike,
    attr::MOI.ConstraintDualStart,
    bridge::SplitIntervalBridge{T},
    value,
) where {T}
    if bridge.lower === nothing
        if bridge.upper !== nothing
            MOI.set(model, attr, bridge.upper, value)
        end
    else
        if bridge.upper === nothing
            MOI.set(model, attr, bridge.lower, value)
        else
            lower, upper = _split_dual_start(value)
            MOI.set(model, attr, bridge.lower, lower)
            MOI.set(model, attr, bridge.upper, upper)
        end
    end
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintBasisStatus,
    bridge::SplitIntervalBridge,
)
    if bridge.upper !== nothing
        upper_stat = MOI.get(model, attr, bridge.upper)
        if upper_stat == MOI.NONBASIC
            return MOI.NONBASIC_AT_UPPER
        end
    end
    if bridge.lower === nothing
        if bridge.upper === nothing
            # The only case where the interval `[-∞, ∞]` is allowed is for
            # `VariableIndex` constraints but `ConstraintBasisStatus` is not
            # defined for `VariableIndex` constraints.
            _error_double_inf(attr)
        else
            return upper_stat
        end
    else
        lower_stat = MOI.get(model, attr, bridge.lower)
        if lower_stat == MOI.NONBASIC
            return MOI.NONBASIC_AT_LOWER
        end
        # Both statuses must be BASIC or SUPER_BASIC, so just return the lower.
        return lower_stat
    end
end

function MOI.modify(
    model::MOI.ModelLike,
    bridge::SplitIntervalBridge,
    change::MOI.AbstractFunctionModification,
)
    if bridge.lower !== nothing
        MOI.modify(model, bridge.lower, change)
    end
    if bridge.upper !== nothing
        MOI.modify(model, bridge.upper, change)
    end
    return
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintFunction,
    bridge::SplitIntervalBridge{T,F},
    func::F,
) where {T,F}
    if bridge.lower !== nothing
        MOI.set(model, MOI.ConstraintFunction(), bridge.lower, func)
    end
    if bridge.upper !== nothing
        MOI.set(model, MOI.ConstraintFunction(), bridge.upper, func)
    end
    return
end

function MOI.set(
    model::MOI.ModelLike,
    ::MOI.ConstraintSet,
    bridge::SplitIntervalBridge{T,F,S},
    change::S,
) where {T,F,S}
    lower_set = _lower_set(change)
    upper_set = _upper_set(change)
    if lower_set === nothing && upper_set === nothing
        # The constraints are going to be deleted, we store the function before
        # it is lost.
        bridge.func = MOI.get(model, MOI.ConstraintFunction(), bridge)
    end
    if bridge.lower === nothing
        if lower_set !== nothing
            func = MOI.get(model, MOI.ConstraintFunction(), bridge)
            bridge.lower = MOI.add_constraint(model, func, lower_set)
        end
    else
        if lower_set === nothing
            MOI.delete(model, bridge.lower)
            bridge.lower = nothing
        else
            MOI.set(model, MOI.ConstraintSet(), bridge.lower, lower_set)
        end
    end
    if bridge.upper === nothing
        if upper_set !== nothing
            func = MOI.get(model, MOI.ConstraintFunction(), bridge)
            bridge.upper = MOI.add_constraint(model, func, upper_set)
        end
    else
        if upper_set === nothing
            MOI.delete(model, bridge.upper)
            bridge.upper = nothing
        else
            MOI.set(model, MOI.ConstraintSet(), bridge.upper, upper_set)
        end
    end
    return
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintFunction,
    bridge::SplitIntervalBridge,
)
    if bridge.lower !== nothing
        return MOI.get(model, attr, bridge.lower)
    end
    if bridge.upper !== nothing
        return MOI.get(model, attr, bridge.upper)
    end
    return bridge.func
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SplitIntervalBridge{T,F,MOI.Interval{T}},
) where {T,F}
    if bridge.lower === nothing
        lower = typemin(T)
    else
        lower = MOI.get(model, attr, bridge.lower).lower
    end
    if bridge.upper === nothing
        upper = typemax(T)
    else
        upper = MOI.get(model, attr, bridge.upper).upper
    end
    return MOI.Interval(lower, upper)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SplitIntervalBridge{T,F,MOI.EqualTo{T}},
) where {T,F}
    return MOI.EqualTo(MOI.get(model, attr, bridge.lower).lower)
end

function MOI.get(
    model::MOI.ModelLike,
    attr::MOI.ConstraintSet,
    bridge::SplitIntervalBridge{T,F,MOI.Zeros},
) where {T,F}
    return MOI.Zeros(MOI.get(model, attr, bridge.lower).dimension)
end

max_stars_repo_name: JuliaOpt/MathOptInterface

max_issues_repo_name: JuliaOpt/MathOptInterface

max_forks_repo_name: JuliaOpt/MathOptInterface


text:
module BetterBanner
function banner(io = stdout)
    GIT_VERSION_INFO = Base.GIT_VERSION_INFO
    if GIT_VERSION_INFO.tagged_commit
        commit_string = Base.TAGGED_RELEASE_BANNER
    elseif isempty(GIT_VERSION_INFO.commit)
        commit_string = ""
    else
        days = Int(floor((ccall(:jl_clock_now, Float64, ()) - GIT_VERSION_INFO.fork_master_timestamp) / (60 * 60 * 24)))
        days = max(0, days)
        unit = days == 1 ? "day" : "days"
        distance = GIT_VERSION_INFO.fork_master_distance
        commit = GIT_VERSION_INFO.commit_short

        if distance == 0
            commit_string = "Commit $(commit) ($(days) $(unit) old master)"
        else
            branch = GIT_VERSION_INFO.branch
            commit_string = "$(branch)/$(commit) (fork: $(distance) commits, $(days) $(unit))"
        end
    end
    empty_str = ""
    logoplusspace = 26
    commit_date = isempty(Base.GIT_VERSION_INFO.date_string) ? "" : "($(split(Base.GIT_VERSION_INFO.date_string)[1]))"
    doc = "Documentation:"
    doc_link = "https://docs.julialang.org"
    general_help = """Type "?" for help"""
    pkg_help = """"]?" for Pkg help."""
    _,width = displaysize(io)
    if width >= 67
        t1 = empty_str
        t2 = doc * " " * doc_link
        t3 = empty_str
        t4 = general_help * ", " * pkg_help
        t5 = empty_str
        t6 = "Version $(Base.VERSION) $(commit_date)"
        t7 = commit_string
        t8 = empty_str
        extra = """

        """
    elseif 44 <= width <= 67
        t1 = doc
        if width >= 52
            t2 = doc_link
        else
            t2 = "docs.julialang.org"
        end
        t3 = empty_str
        t4 = general_help * ","
        t5 = pkg_help
        t6 = empty_str
        version_str = "Ver. $(Base.VERSION) $commit_date"
        if length(version_str) + logoplusspace < width
            t7 = version_str
            t8 = empty_str
        else
            t7 = "Ver. $(Base.VERSION)"
            t8 = commit_date
        end
        extra = """

        $commit_string

        """
    else
        t1 = empty_str
        t2 = empty_str
        t3 = empty_str
        t4 = empty_str
        t5 = empty_str
        t6 = empty_str
        t7 = empty_str
        t8 = empty_str
        #expect overflowing, but the logo is unaffected
        extra = """

        $doc
        $doc_link

        $(general_help), $pkg_help

        Ver. $(Base.VERSION) $commit_date

        $commit_string

        """
    end

    separator = ifelse(44 <= width, " | ",empty_str)
    tt1 = separator * t1
    tt2 = separator * t2
    tt3 = separator * t3
    tt4 = separator * t4
    tt5 = separator * t5
    tt6 = separator * t6
    tt7 = separator * t7
    tt8 = separator * t8
    if get(io, :color, false)
        c = Base.text_colors
        tx = c[:normal] # text
        jl = c[:normal] # julia
        d1 = c[:bold] * c[:blue]    # first dot
        d2 = c[:bold] * c[:red]     # second dot
        d3 = c[:bold] * c[:green]   # third dot
        d4 = c[:bold] * c[:magenta] # fourth dot

    print("""

                    $(d3)⢰⣿⣿⡆$(tx)   $(tt1)
                    $(d3)⠈⠛⠛⠁$(tx)   $(tt2)
     $(d1)⢰⣿⣿⡆$(jl)       ⣴⣾$(d2)⢰⣿⣿⡆$(d4)⢰⣿⣿⡆$(tx) $(tt3)
     $(d1)⠈⠛⠛⠁$(jl)       ⣿⣿$(d2)⠈⠛⠛⠁$(d4)⠈⠛⠛⠁$(tx) $(tt4)
      ⣿⣿ ⣿⣿  ⣿⣿ ⣿⣿ ⣴⣾ ⣾⠟⠛⣿⣷$(tt5)
      ⣿⣿ ⣿⣿  ⣿⣿ ⣿⣿ ⣿⣿⢀⣴⡾⠋⣿⣿$(tt6)
      ⣿⣿ ⠙⢿⣶⠞⣿⣿ ⣿⣿ ⣿⣿⠈⠻⣷⣶⡿⣿$(tt7)
    ⢶⣦⣿⠟                   $(tt8)
    $(extra)""")
    else
        print("""
                    ⢰⠉⠉⡆   $(tt1)
                    ⠈⠒⠒⠁   $(tt2)
     ⢰⠉⠉⡆       ⣴⣾⢰⠉⠉⡆⢰⠉⠉⡆ $(tt3)
     ⠈⠒⠒⠁       ⣿⣿⠈⠒⠒⠁⠈⠒⠒⠁ $(tt4)
      ⣿⣿ ⣿⣿  ⣿⣿ ⣿⣿ ⣴⣾ ⣾⠟⠛⣿⣷$(tt5)
      ⣿⣿ ⣿⣿  ⣿⣿ ⣿⣿ ⣿⣿⢀⣴⡾⠋⣿⣿$(tt6)
      ⣿⣿ ⠙⢿⣶⠞⣿⣿ ⣿⣿ ⣿⣿⠈⠻⣷⣶⡿⣿$(tt7)
    ⢶⣦⣿⠟                   $(tt8)
    $(extra)""")
    end

end
__precompile__()
function __init__()
    @eval begin
        Base.banner(io::IO=stdout) = BetterBanner.banner(io)
    end
end
end #module


text:
module DPSGD

export DifferentialPrivacy, privacy_spent
export solve_niterations, solve_noise_multiplier

using
    Flux,
    LogarithmicNumbers,
    Parameters,
    Random,
    Statistics,
    Zygote,
    LinearAlgebra

using SpecialFunctions: logerfc, gamma

Base.binomial(n, k) = gamma(n + 1) / (gamma(k + 1) * gamma(n - k + 1))

"""
    DifferentialPrivacy(nexamples; ...)

Clip gradients when their absolute value exceeds `clip_threshold`.
"""
@with_kw mutable struct DifferentialPrivacy
    clip_threshold::Float32 = 1.0
    noise_multiplier::Float32 = 1.0
    lotsize::Int = 0
    nexamples::Int
    niterations::Int = 0
    first_parameters = nothing
    norm::Vector{<:AbstractFloat} = zeros(1)
end

function DifferentialPrivacy(nexamples; lotsize = 0, kw...)
    DifferentialPrivacy(; nexamples = nexamples, lotsize = lotsize, kw...)
end

function calculate_norms!(opt, xs::Params, gs)
    per_parameter_norms = nothing
    num_samples = 0
    grad_float_type = AbstractFloat # Specific type will be inferred from gs below
    for x in xs
        isnothing(gs[x]) && continue
        if num_samples == 0
            num_samples = size(gs[x])[1]
            grad_float_type = eltype(gs[x])
        else
            @assert num_samples == size(gs[x])[1]
        end
        parameter_norm_per_sample = [norm(gs[x][i,:]) for i in 1:num_samples]
        if isnothing(per_parameter_norms)
            per_parameter_norms = parameter_norm_per_sample
        else
            per_parameter_norms = hcat(per_parameter_norms, parameter_norm_per_sample)
        end

    end

    # Norm Across Params
    per_sample_norm = zeros(grad_float_type, num_samples)
    for i in 1:num_samples
        per_sample_norm[i] = norm(per_parameter_norms[i,:])
    end

    # Store in DifferentialPrivacy Struct
    for sub_opt in opt.os
        if isa(sub_opt, DifferentialPrivacy)
            sub_opt.norm = per_sample_norm
        end
    end
end

"""
    train!(loss, ps, data, opt::DifferentialPrivacy; ...)

Calculate the per example gradient when the loss
is a vector with an element per example.
"""
function train!(loss, ps, data, opt; cb = () -> ())
    ps = Zygote.Params(ps)
    cb = Flux.Optimise.runall(cb)

    for d in data
        try
            gs = Zygote.jacobian(ps) do
                loss(Flux.Optimise.batchmemaybe(d)...)
            end
            calculate_norms!(opt, ps, gs)
            Flux.update!(opt, ps, gs)
            cb()
        catch ex
            if ex isa Flux.Optimise.StopException
                break
            elseif ex isa Flux.Optimise.SkipException
                continue
            else
                rethrow(ex)
            end
        end
    end
end

Flux.train!(l, p, d, opt::DifferentialPrivacy; cb = () -> ()) = train!(l, p, d, opt; cb)

# This is definitely type piracy...but not in *spirit*, right? It's like, Samuel
# Bellamy-style type piracy. Flux.Optimiser is the thieving rich, and we're just
# stealing DifferentialPrivacy back.
function Flux.train!(l, p, d, opt::Flux.Optimiser; cb = () -> ())
    f = any(x -> isa(x, DifferentialPrivacy), opt.os) ? train! : Flux.train!
    f(l, p, d, opt; cb)
end

function Flux.Optimise.apply!(opt::DifferentialPrivacy, x, Δ)
    if opt.lotsize ≠ size(Δ, 1)
        if opt.lotsize == 0
            opt.lotsize = size(Δ, 1)
            @debug "Setting lot size to $(opt.lotsize)"
        else
            @warn "Lot size appears to be dynamic, make sure to specify it in the constructor or results will be invalid." maxlog = 1
        end
    end

    # Gradient clipping
    for i in 1:length(opt.norm)
        # Δ_norm = Flux.norm(Δ)
        Δ_norm = opt.norm[i]
        if Δ_norm > opt.clip_threshold
            Δ[i,:] = Flux.rmul!(Δ[i,:], opt.clip_threshold / Δ_norm)
        end
    end
    # Accumulate gradients across examples
    Δ = reshape(Statistics.mean(Δ; dims = 1), size(x))

    # Add noise
    a = similar(Δ)
    Random.randn!(a)
    a *= (opt.noise_multiplier * opt.clip_threshold)
    Δ += a

    # Only increment iterations once per model iteration, not for every
    # parameter set
    if isnothing(opt.first_parameters)
        opt.first_parameters = x
    end

    if opt.first_parameters === x
        opt.niterations += 1
    end

    return Δ
end

function privacy_spent(dp::DifferentialPrivacy, target_delta; kwargs...)
    privacy_spent(target_delta;
        lotsize = dp.lotsize,
        nexamples = dp.nexamples,
        noise_multiplier = dp.noise_multiplier,
        niterations = dp.niterations,
        kwargs...)
end

function privacy_spent(target_delta;
    nexamples, noise_multiplier, niterations,
    lotsize = 1,
    orders = 1:64)
    q = lotsize / nexamples
    rdp = renyi_differential_privacy.((q,), (noise_multiplier,), orders; niterations)

    if length(orders) ≠ length(rdp)
        throw(ArgumentError("Input lists must have the same length."))
    end

    eps = (rdp
           .-
           ((log.(target_delta) .+ log.(orders)) ./ (orders .- 1))
           .+
           (log.((orders .- 1) ./ orders)))

    # Special case when there is no privacy.
    all(isnan, eps) && return Inf

    replace!(eps, NaN => Inf)
    return minimum(eps)
end

function renyi_differential_privacy(q, σ, α; niterations = 1)
    iszero(q) && return 0.0
    iszero(σ) && return Inf # No privacy.
    α == Inf && return Inf
    isone(q) && return α / (2 * σ^2)

    return log(compute_A(q, σ, α)) / (α - 1) * niterations
end

# Refers to $A_α$ in RDP.
function compute_A(q, σ, α)
    if isinteger(α)
        compute_A(Integer, q, σ, Integer(α))
    else
        compute_A(Real, q, σ, α)
    end
end

# Notes on `Logarithmic`:
# - `Logarithmic(x)` is the number `x`, but represented inside the structure
#    by `log(x)` for the sake of accurate arithmetic when `x` is large.
# - Calling `log` on a `Logarithmic{T}` number simply unwraps the value already
#   being used to represent the number (it returns a value of type `T`).
# - If you know `a = log(x)`, then `Logarithmic(x)` can be constructed
#   with no calculation using `exp(Logarithmic, a)`.
# - Negative values *can* be represented with `Logarithmic` without doing
#   anything special--it handles managing the sign for you.

# https://arxiv.org/pdf/1908.10530.pdf
# Section 3.3, Case I: Integer α
function compute_A(::Type{Integer}, q, σ, α::Integer)
    q = Logarithmic(q)
    sum(0:α; init = Logarithmic(0)) do i
        c_i = binomial(α, i) * q^i * (1 - q)^(α - i)
        c_i * exp(Logarithmic, (i^2 - i) / (2σ^2))
    end
end

# Section 3.3, Case II: Fractional
function compute_A(::Type{Real}, q, σ, α; tol = 1e-14)
    q = Logarithmic(q)
    # The two parts of A_alpha are integrals over (-∞, z] and [z, ∞).
    z = 1 / 2 + ((1 - q) / q)σ^2

    # For integers, `binomial(n, k)` for `k>n` is zero, but that is not true for
    # non-integer values of `k`. So this is an infinite sum, carried out to the
    # desired precision.
    i = 0
    total = Logarithmic(0.0)
    last_term = tol + 1.0
    while last_term > tol
        b = binomial(α, i)

        c1 = b * q^i * (1 - q)^(α - i)
        c2 = b * q^(α - i) * (1 - q)^i

        e1 = 1 / 2 * exp(Logarithmic, logerfc((i - z) / (√2σ)))
        e2 = 1 / 2 * exp(Logarithmic, logerfc((z - (α - i)) / (√2σ)))

        s1 = c1 * e1 * exp(Logarithmic, (i^2 - i) / (2σ^2))
        s2 = c2 * e2 * exp(Logarithmic, ((α - i)^2 - (α - i)) / (2σ^2))

        last_term = s1 + s2
        total += last_term
        i += 1
    end
    return total
end

function solve_noise_multiplier(; epsilon, delta,
    nexamples, niterations,
    lotsize = 1, orders = 1:64, lo = 1e-15, hi = 1e16, tol = 1e-6,
    kwargs...)
    # Adding less noise to get the same privacy is better, so find the minimum σ
    # over all `orders`.
    minimum(orders) do order
        # The necessary noise for this `order` is somewhere in between these
        # bounds, so we're forced to keep σ_hi.
        σ_lo, σ_hi = binary_search(epsilon; lo = lo, hi = hi,
                                            tol = tol, lt = (>), kwargs...) do σ
            privacy_spent(delta;
                lotsize = lotsize,
                nexamples = nexamples,
                noise_multiplier = σ,
                orders = order:order,
                niterations = niterations)
        end
        σ_hi
    end
end

function solve_niterations(; epsilon, delta,
    nexamples, noise_multiplier,
    lotsize = 1, orders = 1:64, lo = 0, hi = 2^63 - 1, tol = 1,
    kwargs...)
    # See comments for `solve_noise_multiplier`. The logic is the same, except
    # reversed because being able to do more iterations is better.
    maximum(orders) do order
        n_lo, n_hi = binary_search(epsilon; lo = lo, hi = hi, tol = tol, kwargs...) do n
            privacy_spent(delta;
                lotsize = lotsize,
                nexamples = nexamples,
                noise_multiplier = noise_multiplier,
                orders = order:order,
                niterations = round(Int, n))
        end
        isfinite(n_lo) ? floor(Int, n_lo) : -Inf
    end
end

function binary_search(f, y;
    lo = -1e16,
    hi = 1e16,
    maxtries = 10000,
    tol = 1e-12,
    atol = 1e-12,
    lt = isless)
    for i in 1:maxtries
        (hi - lo ≤ tol) && break

        x = (lo + hi) / 2

        yhat = f(x)
        if isapprox(yhat, y; atol = atol)
            break
        end

        # Since the caller can override `lt`, assume a monotonically increasing
        # function.
        if lt(yhat, y)
            # f(x) is too small, so x is the new lower bound
            lo = x
        else
            hi = x
        end
    end

    return if lt(f(lo), y) && lt(y, f(hi))
        (lo, hi)
    elseif lt(f(hi), y)
        (hi, Inf)
    else
        (-Inf, lo)
    end
end

end # module


text:
if Sys.iswindows()
    url = "https://github.com/simulkade/PhreeqcRM/releases/download/3.7/IPhreeqcd.dll"
    lib_path = joinpath(@__DIR__, "..", "lib", "libiphreeqc.dll")
    if isfile(lib_path)
        println("libiphreeqc.dll exists. Skiping download ...")
    else
        try
            println("Downloading the library files ...")
            download(url, joinpath(@__DIR__, "..", "lib", "libiphreeqc.dll"))
            println("PhreeqcRM.dll downloaded successfully.")
        catch
            @warn "Could not downlaod the library file PhreeqcRM.dll"
            @warn "Download it from: $url"
            @warn "and copy it into the IPhreeqc/lib folder"
        end
    end
elseif Sys.islinux()
    lib_path = joinpath(@__DIR__, "..", "lib", "libphreeqcrm.so")
    url = "https://github.com/simulkade/PhreeqcRM/releases/download/3.7/libiphreeqc-3.7.0.so"
    if isfile(lib_path)
        println("PhreeqcRM.dll exists. Skiping download ...")
    else
        try
            println("Downloading the library files ...")
            download(url, joinpath(@__DIR__, "..", "lib", "libiphreeqc.so"))
            println("libphreeqcrm.so downloaded successfully.")
        catch
            @warn "Could not downlaod the library file libiphreeqc.so"
            @warn "Download it from: $url"
            @warn "and copy it into the IPhreeqc/lib folder"
        end
    end
else
    @error "No binary files for your OS. Only Windows and Linux at the moment!"
end


text:
#####
##### dynamics model
#####
struct AircraftDynamics{TF, TAF}
    aircraft::AS.Aircraft{TF,TF,TF,TF,TAF}
    parameters::AS.ParamsSolveVLMBEM{TF}
    environment::AS.Environment{TF}
    stepsymbol::Union{String, LS.LaTeXString}
    cl_stall::TF
    nx::Int64
    nu::Int64
end

"""
    dynamics(model::Aircraft, X, U)

Describes dynamic response of an aircraft built in AircraftSystems.

# Inputs

* `model::Aircraft` inherits from `RobotDynamics.AbstractModel`
* `X::Vector{Float64}` vector of states in the following order:

    * `X[1]` x; position in the global frame
    * `X[2]` y; position in the global frame
    * `X[3]` vx; x component of velocity in the body frame
    * `X[4]` vy; y component of velocity in the body frame
    * `X[5]` theta; angle from the horizon to the aircraft's attitude
    * `X[6]` theta_dot; time rate of change of theta
    * `X[7]` energy expended [J]

* `U::Vector{Float64}` vector of control inputs in the following order

    * `U[1]` radians/second of the lower rotor
    * `U[2]` radians/second of the upper rotor

# Returns

* `X_dot::Vector{Float64}` vector of time derivatives of states as follows:

    * `X[1]` x_dot; x component of velocity in the global frame
    * `X[2]` y_dot; y component of velocity in the global frame
    * `X[3]` vx_dot; x component of acceleration in the body frame
    * `X[4]` vy_dot; y component of acceleration in the body frame
    * `X[5]` theta_dot; time rate of change of theta
    * `X[6]` theta_dot_dot; time acceleration of theta
    * `X[7]` instantaneous power [W]

"""
function dynamics(model, X, U; rotors_on = true)
    # unpack parameters
    aircraft = model.aircraft
    environment = model.environment # assumed constant
    g = environment.g
    mass = aircraft.inertia_system.mass
    inertia = aircraft.inertia_system.inertia_y

    # calculate resultants in the body frame
    force, M_y, power = aerodynamics!(model, X, U; rotors_on, frame = Body())

    # # calculate power
    # Qs = parameters.Qs[:,1]
    # Ps = Qs .* U
    # P = 2 * sum(Ps)

    return dynamics(X, mass, inertia, g, force, M_y, power)
end

@inline function dynamics(X, mass, inertia, g, force, moment, power)
    theta = X[5]

    s, c = sincos(theta)
    R = ST.@SArray [c -s; s c] # rotation from body to global frame
    v = ST.@SVector [X[3], X[4]]
    theta_dot = X[6]

    p_dot = R * v
    v_dot = ( -g * R' * (ST.@SVector [0,1])  # apply gravity in inertial -z direction
            + force / mass                  # aerodynamic forces
            - theta_dot .* (ST.@SArray [0. -1.; 1. 0.]) * v) # coriolis forces
    theta_dot_dot = moment / inertia

    return ST.@SVector [p_dot[1], p_dot[2], v_dot[1], v_dot[2], theta_dot, theta_dot_dot, power]
end

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
# Test cases of semidefinite programming
using NonconvexSemidefinite, NonconvexIpopt, LinearAlgebra, Test
using Distributions, ChainRulesTestUtils, Random
Random.seed!(1)

# Test setting
mat_dim = 2
mat_length = mat_dim^2
n_sample = 300

@testset "Test matrix to vectors function" begin
    test_rrule(NonconvexSemidefinite.rearrange_x, rand((mat_dim^2-mat_dim)÷2), rand(mat_dim))
end

function random_psd_mat(mat_dim)
    _mat = randn(mat_dim, mat_dim)
    return _mat' * _mat
end

function sd_constraint((x_L, x_D))
    decompress_symmetric(x_L, x_D)
end

@testset "Single semi-definite constraint" begin
    # Randomize groundtruth
    μ, Σ = randn(mat_dim), random_psd_mat(mat_dim)

    # Generate samples
    samples = rand(MvNormal(μ, Σ), n_sample)

    # Objective function
    function f((x_L, x_D))
        -loglikelihood(MvNormal(μ, decompress_symmetric(x_L, x_D)), samples)
    end

    model = Model(f)

    mat_x0 = random_psd_mat(mat_dim)
    x0 = [mat_x0[NonconvexSemidefinite.lowertriangind(mat_x0)], diag(mat_x0)]
    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
    addvar!(model, lbs, ubs)

    add_sd_constraint!(model, sd_constraint)

    alg = SDPBarrierAlg(sub_alg=IpoptAlg())
    options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=200))
    result = optimize(model, alg, x0, options = options)

    minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
    _Σ = sd_constraint(minimizer)

    println("result: \n $result")

    println("minimum: \n $minimum")
    println("minimizer: \n $minimizer")
    println("_Σ: \n $(_Σ)")

    println("Σ: \n $Σ")
    println("abs(_Σ - Σ): \n $(abs.(_Σ - Σ))")
    println("mean(abs(_Σ - Σ)): \n $(mean(abs.(_Σ - Σ)))")

    @test Σ ≈ _Σ rtol = 0.1
end

@testset "Two semi-definite constraints" begin
    # Randomize groundtruth
    μ1, Σ1 = randn(mat_dim), random_psd_mat(mat_dim)
    μ2, Σ2 = randn(mat_dim), random_psd_mat(mat_dim)

    # Generate samples
    samples1 = rand(MvNormal(μ1, Σ1), n_sample)
    samples2 = rand(MvNormal(μ2, Σ2), n_sample)

    # Objective function
    function f((x_L1, x_D1, x_L2, x_D2))
        -loglikelihood(MvNormal(μ1, decompress_symmetric(x_L1, x_D1)), samples1) -
            loglikelihood(MvNormal(μ2, decompress_symmetric(x_L2, x_D2)), samples2)
    end
    model = Model(f)

    mat_x0 = random_psd_mat(mat_dim)
    x0 = [mat_x0[NonconvexSemidefinite.lowertriangind(mat_x0)], diag(mat_x0)]
    x0 = [copy.(x0)..., copy.(x0)...]

    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
    lbs = [copy.(lbs)..., copy.(lbs)...]

    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
    ubs = [copy.(ubs)..., copy.(ubs)...]

    addvar!(model, lbs, ubs)

    add_sd_constraint!(model, x -> sd_constraint(x[1:2]))
    add_sd_constraint!(model, x -> sd_constraint(x[3:4]))

    alg = SDPBarrierAlg(sub_alg=IpoptAlg())
    options = SDPBarrierOptions(sub_options=IpoptOptions(max_iter=200),
                                keep_all=false,
                                c_init=[1.5, 2.0],
                                c_decr=[0.2, 0.1])
    result = optimize(model, alg, x0, options = options)

    minimum, minimizer, optimal_ind = result.minimum, result.minimizer, result.optimal_ind
    _Σ1 = sd_constraint(minimizer[1:2])
    _Σ2 = sd_constraint(minimizer[3:4])

    println("result: \n $result")

    println("minimum: \n $minimum")
    println("minimizer: \n $minimizer")
    println("_Σ1: \n $(_Σ1)")
    println("_Σ2: \n $(_Σ2)")

    println("Σ1: \n $Σ1")
    println("Σ2: \n $Σ2")
    println("abs(_Σ1 - Σ1): \n $(abs.(_Σ1 - Σ1))")
    println("abs(_Σ2 - Σ2): \n $(abs.(_Σ2 - Σ2))")
    println("mean(abs(_Σ1 - Σ1)): \n $(mean(abs.(_Σ1 - Σ1)))")
    println("mean(abs(_Σ2 - Σ2)): \n $(mean(abs.(_Σ2 - Σ2)))")

    @test Σ1 ≈ _Σ1 rtol = 0.1
    @test Σ2 ≈ _Σ2 rtol = 0.1
end

@testset "Different algorithm" begin
    # Randomize groundtruth
    μ, Σ = randn(mat_dim), random_psd_mat(mat_dim)

    # Generate samples
    samples = rand(MvNormal(μ, Σ), n_sample)

    # Passing sd_constraint but not using SDPBarrierOptions
    # Objective function
    function f((x_L, x_D))
        -loglikelihood(MvNormal(μ, decompress_symmetric(x_L, x_D)), samples)
    end
    model = Model(f)

    mat_x0 = random_psd_mat(mat_dim)
    x0 = [mat_x0[NonconvexSemidefinite.lowertriangind(mat_x0)], diag(mat_x0)]
    lbs = [fill(-Inf, length(x0[1])), zeros(length(x0[2]))]
    ubs = [fill(Inf, length(x0[1])), fill(Inf, length(x0[2]))]
    addvar!(model, lbs, ubs)
    add_sd_constraint!(model, sd_constraint)

    @test_logs (:warn, "The solver used does not support semidefinite constraints so they will be ignored.") optimize(model, IpoptAlg(), x0, options=IpoptOptions(max_iter=1, first_order=true))
end


text:
import Mimi.compinstance

"""
    compute_scc(m::Model=get_model(); year::Int = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights_normalization_region::Int=0)

Deprecated function for calculating the social cost of carbon for a MimiFUND model. Use `compute_sc` or gas-specific functions `compute_scco2`, `compute_scch4`, `compute_scn2o`, or `compute_scsf6` instead.
"""
function compute_scc(m::Model=get_model(); year::Union{Int, Nothing} = nothing, gas::Symbol = :CO2, last_year::Int = 3000, equity_weights::Bool = false, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights_normalization_region::Int=0)
    Base.depwarn("`compute_scc` is deprecated. Use `compute_sc` or other gas-specific functions instead.", :MimiFUND)
    year === nothing ? error("Must specify an emission year. Try `compute_scc(year=2020)`.") : nothing
    return compute_sc(m, year = year, gas = gas, last_year = last_year, equity_weights = equity_weights, eta = eta, prtp = prtp, equity_weights_normalization_region = equity_weights_normalization_region)
end

"""
    compute_scco2(m::Model=get_model();
        year::Union{Int, Nothing} = nothing,
        eta::Float64 = 1.45,
        prtp::Float64 = 0.015,
        equity_weights::Bool = false,
        equity_weights_normalization_region::Int = 0,
        last_year::Int = 3000,
        pulse_size::Float64 = 1e7,
        return_mm::Bool = false,
        n::Union{Int, Nothing} = nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Returns the Social Cost of CO2 for the specified `year` for the provided MimiFUND model `m`.
If no model is provided, the default model from MimiFUND.get_model() is used.
Units of the returned value are 1995\$ per metric tonne of CO2.

This is a wrapper function that calls the generic social cost function `compute_sc(m, gas = :CO2, args...)`. See docstring for
`compute_sc` for a full description of the available keyword arguments.
"""
function compute_scco2(m::Model=get_model(); year::Union{Int, Nothing} = nothing, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights::Bool = false, equity_weights_normalization_region::Int=0, last_year::Int = 3000, pulse_size::Float64 = 1e7,
    return_mm::Bool = false, n::Union{Int, Nothing} = nothing, trials_output_filename::Union{String, Nothing} = nothing, seed::Union{Int, Nothing} = nothing)

    year === nothing ? error("Must specify an emission year. Try `compute_scco2(year=2020)`.") : nothing
    return compute_sc(m, gas = :CO2, year = year, eta = eta, prtp = prtp, equity_weights = equity_weights, last_year = last_year, pulse_size = pulse_size,
                        return_mm = return_mm, n = n, trials_output_filename = trials_output_filename, seed = seed, equity_weights_normalization_region=equity_weights_normalization_region)
end

"""
    compute_scch4(m::Model=get_model();
        year::Union{Int, Nothing} = nothing,
        eta::Float64 = 1.45,
        prtp::Float64 = 0.015,
        equity_weights::Bool = false,
        equity_weights_normalization_region::Int=0,
        last_year::Int = 3000,
        pulse_size::Float64 = 1e7,
        return_mm::Bool = false,
        n::Union{Int, Nothing} = nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Returns the Social Cost of CH4 for the specified `year` for the provided MimiFUND model `m`.
If no model is provided, the default model from MimiFUND.get_model() is used.
Units of the returned value are 1995\$ per metric tonne of CH4.

This is a wrapper function that calls the generic social cost function `compute_sc(m, gas = :CH4, args...)`. See docstring for
`compute_sc` for a full description of the available keyword arguments.
"""
function compute_scch4(m::Model=get_model(); year::Union{Int, Nothing} = nothing, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights::Bool = false, last_year::Int = 3000, pulse_size::Float64 = 1e7,
    return_mm::Bool = false, n::Union{Int, Nothing} = nothing, trials_output_filename::Union{String, Nothing} = nothing, seed::Union{Int, Nothing} = nothing, equity_weights_normalization_region::Int=0)

    year === nothing ? error("Must specify an emission year. Try `compute_scch4(year=2020)`.") : nothing
    return compute_sc(m, gas = :CH4, year = year, eta = eta, prtp = prtp, equity_weights = equity_weights, last_year = last_year, pulse_size = pulse_size,
                        return_mm = return_mm, n = n, trials_output_filename = trials_output_filename, seed = seed, equity_weights_normalization_region=equity_weights_normalization_region)
end

"""
    compute_scn2o(m::Model=get_model();
        year::Union{Int, Nothing} = nothing,
        eta::Float64 = 1.45,
        prtp::Float64 = 0.015,
        equity_weights::Bool = false,
        equity_weights_normalization_region::Int = 0,
        last_year::Int = 3000,
        pulse_size::Float64 = 1e7,
        return_mm::Bool = false,
        n::Union{Int, Nothing} = nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Returns the Social Cost of N2O for the specified `year` for the provided MimiFUND model `m`.
If no model is provided, the default model from MimiFUND.get_model() is used.
Units of the returned value are 1995\$ per metric tonne of N2O.

This is a wrapper function that calls the generic social cost function `compute_sc(m, gas = :N2O, args...)`. See docstring for
`compute_sc` for a full description of the available keyword arguments.
"""
function compute_scn2o(m::Model=get_model(); year::Union{Int, Nothing} = nothing, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights::Bool = false, last_year::Int = 3000, pulse_size::Float64 = 1e7,
    return_mm::Bool = false, n::Union{Int, Nothing} = nothing, trials_output_filename::Union{String, Nothing} = nothing, seed::Union{Int, Nothing} = nothing, equity_weights_normalization_region::Int=0)

    year === nothing ? error("Must specify an emission year. Try `compute_scn2o(year=2020)`.") : nothing
    return compute_sc(m, gas = :N2O, year = year, eta = eta, prtp = prtp, equity_weights = equity_weights, last_year = last_year, pulse_size = pulse_size,
                        return_mm = return_mm, n = n, trials_output_filename = trials_output_filename, seed = seed, equity_weights_normalization_region=equity_weights_normalization_region)
end

"""
    compute_scsf6(m::Model=get_model();
        year::Union{Int, Nothing} = nothing,
        eta::Float64 = 1.45,
        prtp::Float64 = 0.015,
        equity_weights::Bool = false,
        equity_weights_normalization_region::Int = 0,
        last_year::Int = 3000,
        pulse_size::Float64 = 1e7,
        return_mm::Bool = false,
        n::Union{Int, Nothing} = nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Returns the Social Cost of SF6 for the specified `year` for the provided MimiFUND model `m`.
If no model is provided, the default model from MimiFUND.get_model() is used.
Units of the returned value are 1995\$ per metric tonne of SF6.

This is a wrapper function that calls the generic social cost function `compute_sc(m, gas = :SF6, args...)`. See docstring for
`compute_sc` for a full description of the available keyword arguments.
"""
function compute_scsf6(m::Model=get_model(); year::Union{Int, Nothing} = nothing, eta::Float64 = 1.45, prtp::Float64 = 0.015, equity_weights::Bool = false, last_year::Int = 3000, pulse_size::Float64 = 1e7,
    return_mm::Bool = false, n::Union{Int, Nothing} = nothing, trials_output_filename::Union{String, Nothing} = nothing, seed::Union{Int, Nothing} = nothing, equity_weights_normalization_region::Int=0)

    year === nothing ? error("Must specify an emission year. Try `compute_scsf6(year=2020)`.") : nothing
    return compute_sc(m, gas = :SF6, year = year, eta = eta, prtp = prtp, equity_weights = equity_weights, last_year = last_year, pulse_size = pulse_size,
                        return_mm = return_mm, n = n, trials_output_filename = trials_output_filename, seed = seed, equity_weights_normalization_region=equity_weights_normalization_region)
end

"""
    compute_sc(m::Model=get_model();
        gas::Symbol = :CO2,
        year::Union{Int, Nothing} = nothing,
        eta::Float64 = 1.45,
        prtp::Float64 = 0.015,
        equity_weights::Bool = false,
        equity_weights_normalization_region::Int=0,
        last_year::Int = 3000,
        pulse_size::Float64 = 1e7,
        return_mm::Bool = false,
        n::Union{Int, Nothing} = nothing,
        trials_output_filename::Union{String, Nothing} = nothing,
        seed::Union{Int, Nothing} = nothing)

Returns the Social Cost of CO2 (or other gas if specified) for the specified `year`
for the provided MimiFUND model `m`. If no model is provided, the default model from MimiFUND.get_model() is used.
Units of the returned value are 1995\$ per metric tonne of the specified gas.

The discount factor is computed from the specified `eta` and pure rate of time preference `prtp`.
Optional regional equity weighting can be used by specifying `equity_weights = true`. If equity weights are used, one can specify
the normalization region used for equity weighting with the `equity_weights_normalization_region` parameter. A value of `0` uses world
average per capita values for the normalization, any other value uses per capita values of that region as the normalization.
By default, the social cost includes damages through the year 3000, but this time horizon can be modified
by setting the keyword `last_year` to some other year within the model's time index.
The size of the marginal emission pulse can be modified with the `pulse_size` keyword argument, in metric
tonnes of the specified gas (this does not change the units of the returned value, which is always normalized
by the `pulse_size` used).

If `return_mm` is set to `true`, then a NamedTuple (sc = sc, mm = mm) of the social cost value and the MarginalModel
used to compute it is returned.

By default, `n = nothing`, and a single value for the "best guess" social cost is returned. If a positive
value for keyword `n` is specified, then a Monte Carlo simulation with sample size `n` will run, sampling from
all of FUND's random variables, and a vector of `n` social cost values will be returned.
Optionally providing a CSV file path to `trials_output_filename` will save all of the sampled trial data as a CSV file.
Optionally providing a `seed` value will set the random seed before running the simulation, allowing the
results to be replicated.
"""
function compute_sc(m::Model=get_model();
    gas::Symbol = :CO2,
    year::Union{Int, Nothing} = nothing,
    eta::Float64 = 1.45,
    prtp::Float64 = 0.015,
    equity_weights::Bool = false,
    equity_weights_normalization_region::Int = 0,
    last_year::Int = 3000,
    pulse_size::Float64 = 1e7,
    return_mm::Bool = false,
    n::Union{Int, Nothing} = nothing,
    trials_output_filename::Union{String, Nothing} = nothing,
    seed::Union{Int, Nothing} = nothing
    )

    year === nothing ? error("Must specify an emission year. Try `compute_sc(year=2020)`.") : nothing
    !(last_year in 1950:3000) ? error("Invlaid value for `last_year`: $last_year. `last_year` must be within the model's time index 1950:3000.") : nothing
    !(year in 1950:last_year) ? error("Invalid value for `year`: $year. `year` must be within the model's time index 1950:$last_year.") : nothing

    mm = get_marginal_model(m; year = year, gas = gas, pulse_size = pulse_size)

    ntimesteps = getindexfromyear(last_year)

    if n === nothing
        # Run the "best guess" social cost calculation
        run(mm; ntimesteps = ntimesteps)
        sc = _compute_sc_from_mm(mm, year = year, gas = gas, ntimesteps = ntimesteps, equity_weights = equity_weights, eta = eta, prtp = prtp, equity_weights_normalization_region=equity_weights_normalization_region)
    elseif n < 1
        error("Invalid n = $n. Number of trials must be a positive integer.")
    else
        # Run a Monte Carlo simulation
        simdef = getmcs()
        payload = (Array{Float64, 1}(undef, n), year, gas, ntimesteps, equity_weights, equity_weights_normalization_region, eta, prtp) # first item is an array to hold SC values calculated in each trial
        Mimi.set_payload!(simdef, payload)
        seed !== nothing ? Random.seed!(seed) : nothing
        si = run(simdef, mm, n, ntimesteps = ntimesteps, post_trial_func = _fund_sc_post_trial, trials_output_filename = trials_output_filename)
        sc = Mimi.payload(si)[1]
    end

    if return_mm
        return (sc = sc, mm = mm)
    else
        return sc
    end
end

# helper function for computing SC from a MarginalModel that's already been run, not to be exported
function _compute_sc_from_mm(mm::MarginalModel; year::Int, gas::Symbol, ntimesteps::Int, equity_weights::Bool, equity_weights_normalization_region::Int, eta::Float64, prtp::Float64)

    # Calculate the marginal damage between run 1 and 2 for each year/region
    marginaldamage = mm[:impactaggregation, :loss]

    ypc = mm.base[:socioeconomic, :ypc]

    # Compute discount factor with or without equityweights
    df = zeros(ntimesteps, 16)
    if !equity_weights
        for r = 1:16
            x = 1.
            for t = getindexfromyear(year):ntimesteps
                df[t, r] = x
                gr = (ypc[t, r] - ypc[t - 1, r]) / ypc[t - 1,r]
                x = x / (1. + prtp + eta * gr)
            end
        end
    else
        normalization_ypc = equity_weights_normalization_region==0 ? mm.base[:socioeconomic, :globalypc][getindexfromyear(year)] : ypc[getindexfromyear(year), equity_weights_normalization_region]
        df = Float64[t >= getindexfromyear(year) ? (normalization_ypc / ypc[t, r]) ^ eta / (1.0 + prtp) ^ (t - getindexfromyear(year)) : 0.0 for t = 1:ntimesteps, r = 1:16]
    end

    # Compute global social cost
    sc = sum(marginaldamage[2:ntimesteps, :] .* df[2:ntimesteps, :])   # need to start from second value because first value is missing
    return sc
end

# Post trial function used for computing monte carlo vector of social cost values
function _fund_sc_post_trial(sim::SimulationInstance, trialnum::Int, ntimesteps::Int, tup::Union{Tuple, Nothing})
    mm = sim.models[1]  # get the already-run MarginalModel
    (sc_results, year, gas, ntimesteps, equity_weights, equity_weights_normalization_region, eta, prtp) = Mimi.payload(sim)  # unpack the payload information
    sc = _compute_sc_from_mm(mm, year = year, gas = gas, ntimesteps = ntimesteps, equity_weights = equity_weights, eta = eta, prtp = prtp, equity_weights_normalization_region=equity_weights_normalization_region)
    sc_results[trialnum] = sc
end

"""
    get_marginal_model(m::Model = get_model(); gas::Symbol = :CO2, year::Int = nothing, pulse_size::Float64 = 1e7)

Creates a Mimi MarginalModel where the provided m is the base model, and the marginal model has additional emissions of gas `gas` in year `year`.
If no year is provided, the marginal emissions component will be added without an additional pulse.
If no Model m is provided, the default model from MimiFUND.get_model() is used as the base model.
The size of the marginal emission pulse can be modified with the `pulse_size` keyword argument, in metric tonnes.
"""
function get_marginal_model(m::Model = get_model(); gas::Symbol = :CO2, year::Union{Int, Nothing} = nothing, pulse_size::Float64 = 1e7)
    year !== nothing && !(year in 1950:3000) ? error("Cannot add marginal emissions in $year, year must be within the model's time index 1950:3000.") : nothing

    mm = create_marginal_model(m, pulse_size)
    add_marginal_emissions!(mm.modified, year; gas = gas, pulse_size = pulse_size)

    return mm
end

# A component for an emissions pulse to be used in social cost calculations. Computes the `output` vector by adding
#   `add` to `input`. This is similar to the Mimi.adder component, except that it allows missing values to be passed through.
@defcomp emissionspulse begin
    add    = Parameter(index=[time])
    input  = Parameter(index=[time])
    output = Variable(index=[time])

    function run_timestep(p, v, d, t)
        v.output[t] = Mimi.@allow_missing(p.input[t]) + p.add[t]
    end
end

"""
Adds an emissionspulse component to `m`, and sets the additional emissions if a year is specified.
The size of the marginal emission pulse can be modified with the `pulse_size` keyword argument, in metric tonnes.
"""
function add_marginal_emissions!(m, year::Union{Int, Nothing} = nothing; gas::Symbol = :CO2, pulse_size::Float64 = 1e7)

    # Add additional emissions to m
    add_comp!(m, emissionspulse, before = :climateco2cycle)
    nyears = length(Mimi.time_labels(m))
    addem = zeros(nyears)
    if year != nothing
        # pulse is spread over ten years, and emissions components is in Mt so divide by 1e7, and convert from CO2 to C if gas==:CO2 because emissions component is in MtC
        addem[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / 1e7 * (gas == :CO2 ? 12/44 : 1)
    end
    set_param!(m, :emissionspulse, :add, addem)

    # Reconnect the appropriate emissions in m
    if gas == :CO2
        connect_param!(m, :emissionspulse, :input, :emissions, :mco2)
        connect_param!(m, :climateco2cycle, :mco2, :emissionspulse, :output)
    elseif gas == :CH4
        connect_param!(m, :emissionspulse, :input, :emissions, :globch4)
        connect_param!(m, :climatech4cycle, :globch4, :emissionspulse, :output)
    elseif gas == :N2O
        connect_param!(m, :emissionspulse, :input, :emissions, :globn2o)
        connect_param!(m, :climaten2ocycle, :globn2o, :emissionspulse, :output)
    elseif gas == :SF6
        connect_param!(m, :emissionspulse, :input, :emissions, :globsf6)
        connect_param!(m, :climatesf6cycle, :globsf6, :emissionspulse, :output)
    else
        error("Unknown gas: $gas")
    end

end

"""
Helper function to set the marginal emissions in the specified year.
"""
function perturb_marginal_emissions!(m::Model, year; comp_name::Symbol = :emissionspulse, pulse_size::Float64 = 1e7, gas::Symbol = :CO2)

    ci = compinstance(m, comp_name)
    emissions = Mimi.get_param_value(ci, :add)

    nyears = length(Mimi.dimension(m, :time))
    new_em = zeros(nyears)
    new_em[getindexfromyear(year):getindexfromyear(year) + 9] .= pulse_size / 1e7 * (gas == :CO2 ? 12/44 : 1)
    emissions[:] = new_em

end

"""
Returns a matrix of marginal damages per one metric tonne of additional emissions of the specified gas in the specified year.
"""
function getmarginaldamages(; year=2020, parameters = nothing, gas = :CO2, pulse_size::Float64 = 1e7)

    # Get marginal model
    m = get_model(params = parameters)
    mm = get_marginal_model(m, year = year, gas = gas, pulse_size = pulse_size)
    run(mm)

    # Return marginal damages
    return mm[:impactaggregation, :loss]
end


text:
struct combo_item
    p::Clonglong
    w::Clonglong
    x::Cint
    i::Cint
end

function solve_monoBinary(prob::Problem, λ::Vector{Float64}, assignment::Assignment)

	weightedProfits = [sum(λ .* prob.profits[1:end, iter]) for iter=1:prob.nbVar]

	if prob.nbVar - assignment.assignEndIndex == 0
		return Solution(Float64.(assignment.assign), assignment.profit, assignment.weight, true), true
	else

		test = true
		for iter = (assignment.assignEndIndex+1):prob.nbVar
			test = test && prob.weights[iter] > (prob.maxWeight - assignment.weight)
		end

		if test
			return Solution(append!(Float64.(assignment.assign[1:assignment.assignEndIndex]), zeros(Float64, prob.nbVar-assignment.assignEndIndex)),
							assignment.profit,
							assignment.weight,
							true), true
		elseif sum(prob.weights[(assignment.assignEndIndex+1):end]) <= (prob.maxWeight - assignment.weight)

			return Solution(append!(Float64.(assignment.assign[1:assignment.assignEndIndex]), ones(Float64, prob.nbVar-assignment.assignEndIndex)),
							assignment.profit + [sum(prob.profits[iter, (assignment.assignEndIndex+1):end]) for iter=1:prob.nbObj],
							assignment.weight + sum(prob.weights[(assignment.assignEndIndex+1):end]),
							true), true
		end
	end

	items = [combo_item(weightedProfits[iter], prob.weights[iter], 0, iter) for iter = (assignment.assignEndIndex+1):prob.nbVar]

	z = ccall((:solve, comboPath),
    Clonglong,
    (Ref{combo_item}, Cint, Clonglong, Clonglong, Clonglong),
    items, (prob.nbVar - assignment.assignEndIndex), (prob.maxWeight - assignment.weight), 0, 0)

    if z == 0
        z = dot(weightedProfits[(assignment.assignEndIndex+1):end], broadcast(it->it.x, items))
    end

    x = Float64.(append!(assignment.assign[1:assignment.assignEndIndex], falses(prob.nbVar-assignment.assignEndIndex)))
    for it in items
    	x[it.i] = Float64(it.x)
    end

    return Solution(x,
    				assignment.profit + [sum(x[(assignment.assignEndIndex+1):end] .* prob.profits[iter, (assignment.assignEndIndex+1):end]) for iter = 1:prob.nbObj],
    				assignment.weight + sum(prob.weights[(assignment.assignEndIndex+1):end] .* x[(assignment.assignEndIndex+1):end]),
    				true), true
end

max_stars_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_issues_repo_name: LucasBaussay/BranchAndBoundBiObjKp

max_forks_repo_name: LucasBaussay/BranchAndBoundBiObjKp


text:
abstract AbstractTimedata
abstract AbstractTimenum <: AbstractTimedata
abstract AbstractTimematr <: AbstractTimenum