text:
export iauAtco13
"""
ICRS RA,Dec to observed place.  The caller supplies UTC, site
coordinates, ambient air conditions and observing wavelength.

SOFA models are used for the Earth ephemeris, bias-precession-
nutation, Earth orientation and refraction.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rc,dc  double   ICRS right ascension at J2000.0 (radians, Note 1)
   pr     double   RA proper motion (radians/year; Note 2)
   pd     double   Dec proper motion (radians/year)
   px     double   parallax (arcsec)
   rv     double   radial velocity (km/s, +ve if receding)
   utc1   double   UTC as a 2-part...
   utc2   double   ...quasi Julian Date (Notes 3-4)
   dut1   double   UT1-UTC (seconds, Note 5)
   elong  double   longitude (radians, east +ve, Note 6)
   phi    double   latitude (geodetic, radians, Note 6)
   hm     double   height above ellipsoid (m, geodetic, Notes 6,8)
   xp,yp  double   polar motion coordinates (radians, Note 7)
   phpa   double   pressure at the observer (hPa = mB, Note 8)
   tc     double   ambient temperature at the observer (deg C)
   rh     double   relative humidity at the observer (range 0-1)
   wl     double   wavelength (micrometers, Note 9)

Returned:
   aob    double*  observed azimuth (radians: N=0,E=90)
   zob    double*  observed zenith distance (radians)
   hob    double*  observed hour angle (radians)
   dob    double*  observed declination (radians)
   rob    double*  observed right ascension (CIO-based, radians)
   eo     double*  equation of the origins (ERA-GST)

Returned (function value):
         int      status: +1 = dubious year (Note 4)
                           0 = OK
                           -1 = unacceptable date

Notes:

   1.  Star data for an epoch other than J2000.0 (for example from the
   Hipparcos catalog, which has an epoch of J1991.25) will require
   a preliminary call to iauPmsafe before use.

   2.  The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

   3.  utc1+utc2 is quasi Julian Date (see Note 2), apportioned in any
   convenient way between the two arguments, for example where utc1
   is the Julian Day Number and utc2 is the fraction of a day.

   However, JD cannot unambiguously represent UTC during a leap
   second unless special measures are taken.  The convention in the
   present function is that the JD day represents UTC days whether
   the length is 86399, 86400 or 86401 SI seconds.

   Applications should use the function iauDtf2d to convert from
   calendar date and time of day into 2-part quasi Julian Date, as
   it implements the leap-second-ambiguity convention just
   described.

   4.  The warning status "dubious year" flags UTCs that predate the
   introduction of the time scale or that are too far in the
   future to be trusted.  See iauDat for further details.

   5.  UT1-UTC is tabulated in IERS bulletins.  It increases by exactly
   one second at the end of each positive UTC leap second,
   introduced in order to keep UT1-UTC within +/- 0.9s.  n.b. This
   practice is under review, and in the future UT1-UTC may grow
   essentially without limit.

   6.  The geographical coordinates are with respect to the WGS84
   reference ellipsoid.  TAKE CARE WITH THE LONGITUDE SIGN:  the
   longitude required by the present function is east-positive
   (i.e. right-handed), in accordance with geographical convention.

   7.  The polar motion xp,yp can be obtained from IERS bulletins.  The
   values are the coordinates (in radians) of the Celestial
   Intermediate Pole with respect to the International Terrestrial
   Reference System (see IERS Conventions 2003), measured along the
   meridians 0 and 90 deg west respectively.  For many
   applications, xp and yp can be set to zero.

   8.  If hm, the height above the ellipsoid of the observing station
   in meters, is not known but phpa, the pressure in hPa (=mB),
   is available, an adequate estimate of hm can be obtained from
   the expression

         hm = -29.3 * tsl * log ( phpa / 1013.25 );

   where tsl is the approximate sea-level air temperature in K
   (See Astrophysical Quantities, C.W.Allen, 3rd edition, section
   52).  Similarly, if the pressure phpa is not known, it can be
   estimated from the height of the observing station, hm, as
   follows:

         phpa = 1013.25 * exp ( -hm / ( 29.3 * tsl ) );

   Note, however, that the refraction is nearly proportional to
   the pressure and that an accurate phpa value is important for
   precise work.

   9.  The argument wl specifies the observing wavelength in
   micrometers.  The transition from optical to radio is assumed to
   occur at 100 micrometers (about 3000 GHz).

   10. The accuracy of the result is limited by the corrections for
   refraction, which use a simple A*tan(z) + B*tan^3(z) model.
   Providing the meteorological parameters are known accurately and
   there are no gross local effects, the predicted observed
   coordinates should be within 0.05 arcsec (optical) or 1 arcsec
   (radio) for a zenith distance of less than 70 degrees, better
   than 30 arcsec (optical or radio) at 85 degrees and better
   than 20 arcmin (optical) or 30 arcmin (radio) at the horizon.

   Without refraction, the complementary functions iauAtco13 and
   iauAtoc13 are self-consistent to better than 1 microarcsecond
   all over the celestial sphere.  With refraction included,
   consistency falls off at high zenith distances, but is still
   better than 0.05 arcsec at 85 degrees.

   11. "Observed" Az,ZD means the position that would be seen by a
   perfect geodetically aligned theodolite.  (Zenith distance is
   used rather than altitude in order to reflect the fact that no
   allowance is made for depression of the horizon.)  This is
   related to the observed HA,Dec via the standard rotation, using
   the geodetic latitude (corrected for polar motion), while the
   observed HA and RA are related simply through the Earth rotation
   angle and the site longitude.  "Observed" RA,Dec or HA,Dec thus
   means the position that would be seen by a perfect equatorial
   with its polar axis aligned to the Earth's axis of rotation.

   12. It is advisable to take great care with units, as even unlikely
   values of the input parameters are accepted and processed in
   accordance with the models used.

Called:
   iauApco13    astrometry parameters, ICRS-observed, 2013
   iauAtciq     quick ICRS to CIRS
   iauAtioq     quick CIRS to observed

This revision:   2016 February 2

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# int iauAtco13(double rc, double dc,
#               double pr, double pd, double px, double rv,
#               double utc1, double utc2, double dut1,
#               double elong, double phi, double hm, double xp, double yp,
#               double phpa, double tc, double rh, double wl,
#               double *aob, double *zob, double *hob,
#               double *dob, double *rob, double *eo)

function iauAtco13(rc::Real, dc::Real, pr::Real, pd::Real,
                   px::Real, rv::Real, utc1::Real, utc2::Real,
                   dut1::Real, elong::Real, phi::Real, hm::Real,
                   xp::Real, yp::Real, phpa::Real, tc::Real,
                   rh::Real, wl::Real)

   # Allocate return values
   ref_aob = Ref{Float64}(0.0)
   ref_zob = Ref{Float64}(0.0)
   ref_hob = Ref{Float64}(0.0)
   ref_dob = Ref{Float64}(0.0)
   ref_rob = Ref{Float64}(0.0)
   ref_eo  = Ref{Float64}(0.0)

   status = ccall((:iauAtco13, libsofa_c), Cint,
            (Cdouble, Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble, Cdouble, Cdouble,
            Cdouble, Cdouble,
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble},
            Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, rc), convert(Float64, dc),
            convert(Float64, pr), convert(Float64, pd),
            convert(Float64, px), convert(Float64, rv),
            convert(Float64, utc1), convert(Float64, utc2),
            convert(Float64, dut1), convert(Float64, elong),
            convert(Float64, phi), convert(Float64, hm),
            convert(Float64, xp), convert(Float64, yp),
            convert(Float64, phpa), convert(Float64, tc),
            convert(Float64, rh), convert(Float64, wl),
            ref_aob, ref_zob, ref_hob, ref_dob, ref_rob, ref_eo)

   return status, ref_aob[], ref_zob[], ref_hob[], ref_dob[], ref_rob[], ref_eo[]
end


text:
include("molzahn.jl")

function minimize_sliwak(cliques::AbstractVector, cliquetree::AbstractVector, alpha, beta, m, nb_lc)
    costs = [sliwak_heuristic(cliques, edge, alpha, beta, m, nb_lc) for edge in cliquetree]
    index = argmin(costs)
    return cliquetree[index]
end

function sliwak_heuristic(cliques::AbstractVector, edge, alpha, beta, m, nb_lc)
    clique1 = cliques[edge[1]]
    clique2 = cliques[edge[2]]
    size_intersect = length(intersect(clique1, clique2))
    size_clique1 = length(clique1)
    size_clique2 = length(clique2)
    size_fusion_clique = size_clique1 + size_clique2 - size_intersect
    nb_lc_new = nb_lc - size_intersect * (2 * size_intersect + 1)
    q2 = (m + nb_lc_new)^3
    cost = alpha * (- size_clique1^3 - size_clique2^3 + size_fusion_clique^3)
           + beta * q2
    return cost
end

function merge_sliwak(cliques::AbstractVector, cliquetree::AbstractVector;
                      alpha, beta, m, nb_lc,
                      treshold::Union{Int,Function}=treshold_percent(cliques, 0.1),
                      kmax=3, kwargs...)
    edges = Set()
    if typeof(treshold) <: Int
        treshold_value = treshold
        treshold = x -> length(x) >= treshold_value
    end
    k = 0
    while treshold(cliques)
        edge = k <= kmax ?
                    minimize_sliwak(cliques, cliquetree, alpha, beta, nb_lc, m) :
                    minimize_molzahn(cliques, cliquetree)
        size_intersect = length(intersect(cliques[edge[1]], cliques[edge[2]]))
        nb_lc -= size_intersect * (2 * size_intersect + 1)
        edges = vcat(edges..., get_added_edges(cliques[edge[1]], cliques[edge[2]])...)
        cliques, cliquetree = merge_cliques(cliques, cliquetree, edge)
        k += 1
    end
    return cliques, cliquetree, edges
end

ext: jl

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
# =================
# Dependencies
# =================

using ReachabilityAnalysis
using CarlemanLinearization
using Plots
using LaTeXStrings
using Plots.PlotMeasures
using LinearAlgebra
using SparseArrays
using DifferentialEquations # to plot trajectories on top

using CarlemanLinearization: _error_bound_specabs_R

include("../utils.jl")

# =================
# Model definition
# =================

# model params
r_vac = 0.19
r_tra = 0.13
T_lat = 5.2
T_inf = 2.3

# Model params

Λ = 0.0
P = 1e7 # total

@show Ptot = P
@show PS0 = Ptot * 0.6
@show PI0 = Ptot * 0.03
@show PE0 = Ptot - PS0 - PI0

function seir_model_carlin()

    F1 = zeros(3, 3)
    F1[1, 1] = -Λ/P -r_vac
    F1[2, 2] = -Λ/P -1/T_lat
    F1[3, 2] = 1/T_lat
    F1[3, 3] = -Λ/P -1/T_inf

    F2 = zeros(3, 9) # [x, x⊗x]
    F2[1, 3] = -r_tra/P
    F2[2, 3] = r_tra/P;

    return F1, F2
end

# =================
# Solution method
# =================

## Solution with CARLIN

function _solve_seir_carlin(; N=4, T=30.0, δ=0.1, radius0=0, bloat=false, resets=nothing)

    x0c = [PS0, PE0, PI0] # initial condition

    F1, F2 = seir_model_carlin()
    R, Re_λ1 = _error_bound_specabs_R(x0c, F1, F2; check=true)

    n = 3
    dirs = _template(n=n, N=N)
    alg = LGG09(δ=δ, template=dirs)

    if radius0 == 0
        X0 = convert(Hyperrectangle, Singleton(x0c))
    else
        X0 = Hyperrectangle(x0c, radius0)
    end

    if isnothing(resets)
        @time sol = _solve_CARLIN(X0, F1, F2; alg=alg, N=N, T=T, bloat=bloat);
    else
        @time sol = _solve_CARLIN_resets(X0, F1, F2; resets=resets, alg=alg, N=N, T=T, bloat=bloat);
    end
    return sol
end

## Solution with TMJETS

@taylorize function seir!(dx, x, p, t)
    P_S, P_E, P_I = x

    aux1 = Λ / P

    dx[1] = -aux1 * P_S - r_vac * P_S + Λ - r_tra * P_S * P_I / P
    dx[2] = -aux1 * P_E - P_E / T_lat + r_tra * P_S * P_I / P
    dx[3] = -aux1 * P_I + P_E / T_lat - P_I / T_inf
end

function _solve_seir_carlin_TM(; T=30.0, radius0=0, trajectories=-1)
    x0c = [PS0, PE0, PI0] # initial condition
    if radius0 == 0
        X0 = convert(Hyperrectangle, Singleton(x0c))
    else
        X0 = Hyperrectangle(x0c, radius0)
    end
    prob = @ivp(x' = seir!(x), x(0) ∈ X0, dim=3)

    if trajectories == -1
        @time sol_tm = solve(prob, T=T, alg=TMJets())
    else
        sol_tm = solve(prob, T=T, alg=TMJets(), ensemble=true, trajectories=trajectories)
    end
    return sol_tm
end

# ===============
# Results
# ===============

# parameters
Tmax = 30.0
rr0  = fill(1.0 * 1e5, 3)

# taylor models solution
_solve_seir_carlin_TM(T=Tmax, radius0=rr0) # warm-up
time_TM = @elapsed _solve_seir_carlin_TM(T=Tmax, radius0=rr0)

# no error bounds, N = 2
_solve_seir_carlin(N=2, T=Tmax, δ=0.1, radius0=rr0, bloat=false) # warm-up
time_NoError_N2 = @elapsed _solve_seir_carlin(N=2, T=Tmax, δ=0.1, radius0=rr0, bloat=false)

# including error bounds, N = 5
_solve_seir_carlin(N=5, T=Tmax, δ=0.1, radius0=rr0, bloat=true, resets=[4.0]) # warm-up
time_Error_N5 = @elapsed _solve_seir_carlin(N=5, T=Tmax, δ=0.1, radius0=rr0, bloat=true, resets=[4.0])

# save times
print(io, "SEIR, Carleman, no error bound, N=2, $(time_NoError_N2)\n")
print(io, "SEIR, Carleman, including error bound, N=5, $(time_Error_N5)\n")
print(io, "SEIR, TMJets, -, -, $(time_TM)\n")

# ===============
# Figures
# ===============

# figure with NO error bounds
function figure_seir()

    Tmax = 30.0
    rr0  = fill(1.0 * 1e5, 3)
    solN2 = _solve_seir_carlin(N=2, T=Tmax, δ=0.1, radius0=rr0, bloat=false);

    # taylor models solution
    sol_tm = _solve_seir_carlin_TM(T=Tmax, radius0=rr0)

    # trajectories
    sol_tm_traj = _solve_seir_carlin_TM(T=Tmax, radius0=rr0, trajectories=5)

    fig = plot(legend=:topright, xlab=L"\textrm{Time (days)}", ylab=L"\textrm{Individuals } (\times 10^6)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick = ([0, 10, 20, 30], [L"0", L"10", L"20", L"30"]),
               ytick = ([0, 1e6, 2e6, 3e6, 4e6, 5e6, 6e6],
                        [L"0", L"1", L"2", L"3", L"4", L"5", L"6"]),
               xlims=(0, 30),
               ylims=(0, 6e6),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    # carleman linearization solution
    plot!(fig, solN2, vars=(0, 1), lab=L"S~(N = 2)", alpha=1., c=:blue, lc=:blue)
    plot!(fig, sol_tm, vars=(0, 1), lab=L"S~(TM)", alpha=1., c=:aquamarine, lc=:aquamarine)

    plot!(fig, solN2, vars=(0, 2), lab=L"E~(N=2)", alpha=1., c=:red, lc=:red)
    plot!(fig, sol_tm, vars=(0, 2), lab=L"E~(TM)", alpha=1., c=:darksalmon, lc=:darksalmon)

    plot!(fig, solN2, vars=(0, 3), lab=L"I~(N=2)", alpha=1., c=:orange, lc=:orange)
    plot!(fig, sol_tm, vars=(0, 3), lab=L"I~(TM)", alpha=1., c=:darkseagreen, lc=:darkseagreen)

    # lens
    lens!(fig, [18, 22], [7.5e4, 1.75e5], inset = (1, bbox(0.32, 0.34, 0.28, 0.24)),
           tickfont=font(14, "Times"),
           subplot=2,
           xticks=xticklatex([18., 20., 22.], 2),
           yticks=([7.5e4, 1.25e5, 1.75e5], [L"0.075", L"0.125", L"0.175"]))

    return fig
end

fig = figure_seir()
savefig(fig, joinpath(TARGET_FOLDER, "figure_seir_noerror.pdf"))

# figure with error bounds
function figure_seir_errors()

    Tmax = 30.0
    rr0  = fill(1.0 * 1e5, 3)

    N = 5
    solN5 = _solve_seir_carlin(N=N, T=Tmax, δ=0.1, radius0=rr0, bloat=false);

    solN5bloat = _solve_seir_carlin(N=N, T=Tmax, δ=0.1, radius0=rr0, bloat=true, resets=[4.0]);

    fig = plot(legend=:topright, xlab=L"\textrm{Time (days)}", ylab=L"\textrm{Individuals } (\times 10^6)",
               legendfontsize=25,
               tickfont=font(25, "Times"),
               guidefontsize=25,
               xguidefont=font(35, "Times"),
               yguidefont=font(35, "Times"),
               xtick = ([0, 10, 20, 30], [L"0", L"10", L"20", L"30"]),
               ytick = ([0, 1e6, 2e6, 3e6, 4e6, 5e6, 6e6],
                        [L"0", L"1", L"2", L"3", L"4", L"5", L"6"]),
               xlims=(0, 30),
               ylims=(0, 6e6),
               bottom_margin=10mm, left_margin=6mm, right_margin=6mm, top_margin=3mm, size=(900, 600))

    # carleman linearization solution

    plot!(fig, solN5bloat, vars=(0, 1), lab=L"err(S)", alpha=1., c=:blue, lc=:blue)

    plot!(fig, solN5bloat, vars=(0, 2), lab=L"err(E)", alpha=1., c=:red, lc=:red)

    plot!(fig, solN5bloat, vars=(0, 3), lab=L"err(I)", alpha=1., c=:orange, lc=:orange)

    plot!(fig, solN5, vars=(0, 1), lab=L"S~(N = %$N)", alpha=1., c=:aquamarine, lc=:aquamarine)
    plot!(fig, solN5, vars=(0, 2), lab=L"E~(N = %$N)", alpha=1., c=:darksalmon, lc=:darksalmon)
    plot!(fig, solN5, vars=(0, 3), lab=L"I~(N = %$N)", alpha=1., c=:darkseagreen, lc=:darkseagreen)

    ylims!(0, 6e6)

    a = [4.0, 0.0]; b = [4.0, 1.3e6]
    plot!(LineSegment(a, b), markershape=:none, seriestype=:shape, alpha=1., c=:black, ls=:dash, lw=2.0)

    # lens
    lens!(fig, [3.5, 4.5], [0.6e6, 1.3e6], inset = (1, bbox(0.32, 0.34, 0.28, 0.24)),
           tickfont=font(14, "Times"),
           subplot=2,
           xticks=xticklatex([3.5, 4.0, 4.5], 2),
           yticks=([0.7e6, 1.0e6, 1.3e6], [L"0.7", L"1.0", L"1.3"]))

    return fig
end

fig = figure_seir_errors()
savefig(fig, joinpath(TARGET_FOLDER, "figure_seir_error.pdf"))


text:
# using Flux
# using CuArrays

"""
    StateBuffer(size::Int, state_size)
A cicular buffer for states. Typically used for images, can be used for state shapes up to 4d.
"""
mutable struct StateBuffer{T, N, S}
    buffer::Array{T, N}
    cur_idx::Int
    capacity::Int
    full::Bool
    state_size::S
end
StateBuffer{T}(size, state_size) where {T} = begin
    s_bfr = zeros(T, state_size..., size)
    StateBuffer(s_bfr,
                1,
                size,
                false,
                state_size)
end

ImageBuffer(size, img_shape) = StateBuffer{UInt8}(size, img_shape)

Base.isempty(sb::StateBuffer) = (sb.cur_idx == 1 && !sb.full)
isfull(sb::StateBuffer) = sb.full
capacity(sb::StateBuffer) = sb.capacity
Base.length(sb::StateBuffer) = isfull(sb) ? capacity(sb) : (sb.cur_idx - 1)
Base.eltype(sb::StateBuffer) = Base.eltype(sb.buffer)
hist_length(sb::StateBuffer) = 1

@inline Base.@propagate_inbounds function laststate(sb::StateBuffer)
    @boundscheck if length(sb) == 0
        throw(BoundsError(sb, 1))
    end

    if sb.cur_idx == 1
        capacity(sb)
    else
        sb.cur_idx - 1
    end
end

# TODO: Generate w/ macros.
@inline Base.@propagate_inbounds function Base.getindex(sb::StateBuffer, idx) where {A}
    if sb.state_size isa Int || sb.state_size isa Tuple{Int}
        sb.buffer[:, idx]
    elseif sb.state_size isa Tuple{Int, Int}
        sb.buffer[ :, :, idx]
    elseif sb.state_size isa Tuple{Int, Int, Int}
        sb.buffer[:, :, :, idx]
    elseif sb.state_size isa Tuple{Int, Int, Int, Int}
        sb.buffer[:, :, :, :, idx]
    end
end

@inline Base.@propagate_inbounds function Base.view(sb::StateBuffer, idx) where {A}
    if sb.state_size isa Int || sb.state_size isa Tuple{Int}
        view(sb.buffer, :, idx)
    elseif sb.state_size isa Tuple{Int, Int}
        view(sb.buffer, :, :, idx)
    elseif sb.state_size isa Tuple{Int, Int, Int}
        view(sb.buffer, :, :, :, idx)
    elseif sb.state_size isa Tuple{Int, Int, Int, Int}
        view(sb.buffer, :, :, :, :, idx)
    end
end

function Base.push!(sb::StateBuffer, state; new_episode=false)
    view(sb, sb.cur_idx) .= state
    sb.cur_idx += 1
    if sb.cur_idx > sb.capacity
        sb.cur_idx = 1
        sb.full = true
    end
    return sb
end

"""
    HistStateBuffer
A buffer which manages the state buffer, giving lastindex as a history buffer. This is useful when the user wants to deal w/ partial observability by using a short domain.
"""

mutable struct HistStateBuffer{T, N, S, SQ<:Val}
    buffer::StateBuffer{T, N, S}
    hist_length::Int
    hist::Array{Int, 1}
    squeeze::SQ
end

HistStateBuffer{T}(size, state_size, hist_length, squeeze=false) where {T} = if squeeze isa Val
    HistStateBuffer(StateBuffer{T}(size + hist_length, state_size), hist_length, ones(Int, hist_length), squeeze)
else
    HistStateBuffer(StateBuffer{T}(size + hist_length, state_size), hist_length, ones(Int, hist_length), Val(squeeze))
end

@forward HistStateBuffer.buffer Base.length, isfull, capacity, Base.isempty

hist_length(sb::HistStateBuffer) = sb.hist_length

@inline Base.@propagate_inbounds function laststate(sb::HistStateBuffer)
    @boundscheck if length(sb) == 0
        throw(BoundsError(sb, 1))
    end
    copy(sb.hist)
end

@inline Base.@propagate_inbounds function Base.getindex(sb::HistStateBuffer, idx::Array{<:Integer, 1})
    if sb.squeeze isa Val{true}
        # sb.buffer[idx]
        data = sb.buffer[idx]
        reshape(data, length(data))
    elseif sb.squeeze isa Val{false}
        sb.buffer[idx]
    end
end

@inline Base.@propagate_inbounds function Base.getindex(sb::HistStateBuffer, idx::Array{<:Integer, 2})
    if sb.squeeze isa Val{true}
        data = sb.buffer[idx]
        reshape(data, sum(Base.size(data)[1:end-1]), Base.size(data)[end])
    elseif sb.squeeze isa Val{false}
        sb.buffer[idx]
    end
end

function Base.push!(sb::HistStateBuffer, state; new_episode=false)
    push!(sb.buffer, state; new_episode=new_episode)
    hist_length = sb.hist_length
    if new_episode
        sb.hist .= laststate(sb.buffer)
    else
        sb.hist[1:(end - 1)] = sb.hist[2:end]
        sb.hist[end] = laststate(sb.buffer)
    end
    return sb
end


text:
# Note time to run is 30mins to 1.5 hours depending on computer

using OneDimensionalCorrelation

##
b4 = BasisLobatto(-12, 12, 4, 24)
b6 = BasisLobatto(-12, 12, 6, 16)
b8 = BasisLobatto(-12, 12, 8, 12)

# Nuclear potential
V(x) = -2( exp(-0.5(x-2)^2) + exp(-0.5(x+2)^2) )

##

E4, Eorb4, orbitals4 = solve_hartree_fock(b4, V)
E6, Eorb6, orbitals6 = solve_hartree_fock(b6, V)
E8, Eorb8, orbitals8 = solve_hartree_fock(b8, V)

##

# Takes some time ~ 3*5mins on good computer
ci4 = full_ci(b4, orbitals4, V)
ci6 = full_ci(b6, orbitals6, V)
ci8 = full_ci(b8, orbitals8, V)

## Remove boundary states

# A little faster than previous one
rorbs4 = reduced_ci_orbitals(b4, orbitals4; pointlike=true)
cir4 = full_ci(b4, rorbs4, V)

rorbs6 = reduced_ci_orbitals(b6, orbitals8; pointlike=true)
cir6 = full_ci(b6, rorbs6, V)

rorbs8 = reduced_ci_orbitals(b8, orbitals8; pointlike=true)
cir8 = full_ci(b8, rorbs8, V)

##

println("E_hf            = ", round(E6; digits=7))
println("E_fci           = ", round(ci6["energies"][1]; digits=7))

println("E_rci (4 elems) = ", round(cir4["energies"][1]; digits=7))
println("E_rci (6 elems) = ", round(cir6["energies"][1]; digits=7))
println("E_rci (8 elems) = ", round(cir8["energies"][1]; digits=7))

##

# Visualizations

using Plots

# HF ground state
plot(b6, orbitals6[:,1])

# Nuclear potential
plot(b6, V)

# Some virtual orbitals used in RCI6
# Note that orthogonalization makes them not to be
# completely point like
plot(b6, rorbs6[:, [2,20,40,60]])

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
using JuMP, Gurobi

function tcm_estimateForest(N, assortments, transaction_counts, time_limit_overall, perTol, depth_limit, initial_left, initial_right, initial_product, initial_isLeaf)
    M = size(assortments, 1)
    @show M

    time_limit_EM = 60.0

    start_time = time()

    forest_left = deepcopy(initial_left)
    forest_right = deepcopy(initial_right)
    forest_product = deepcopy(initial_product)
    forest_isLeaf = deepcopy(initial_isLeaf)

    K = length(forest_left)

    #
    initial_lambda = 1/K * ones(K)
    lambda, v, loglik, elapsed_time = tcm_forest_EM(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, initial_lambda)

    # @show loglik

	u = transaction_counts ./ v
	u[ isinf.(u)] = 0.0
	u[ isnan.(u)] = 0.0

    # NB - re. isinf.(...), this should not be necessary, if the initial model is such that every item assortment pair has a positive choice probability;
    # this is guaranteed for the ranking model because it starts from a set of rankings corresponding to an independent demand model. If the forest
    # is initialized with the ranking model, then the same condition will hold.
    # for isnan.(...), this is definitely necessary, because transaction_counts and v will both be zero at items that are not part of the assortment.

	# perTol = 1e-7
    iter = 0;

    println("CG iter ", iter, " -- time elapsed = ", 0.0, " log likelihood = ",loglik)

	while (time() - start_time < time_limit_overall)
        iter += 1

		u = transaction_counts ./ v
		u[ isinf.(u)] = 0.0
		u[ isnan.(u)] = 0.0
		single_left, single_right, single_product, single_isLeaf, objval = tcm_solveTreeSubproblem(N, u, assortments, depth_limit)
        # error("Stop here")

		#orderings = [orderings; single_ordering']
        push!(forest_left, single_left)
        push!(forest_right, single_right)
        push!(forest_product, single_product)
        push!(forest_isLeaf, single_isLeaf)

		K += 1

        initial_fudge = 0.999
        initial_lambda = [ initial_fudge*lambda; (1-initial_fudge)]
		lambda, v, new_loglik, elapsed_time = tcm_forest_EM(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, initial_lambda)

        println("CG iter ", iter, " -- time elapsed = ", time() - start_time, " log likelihood = ",new_loglik)

        if ( (new_loglik - loglik) / abs(loglik) < perTol )
        	temp = (new_loglik - loglik) / abs(loglik)
        	@show temp
        	@show new_loglik
        	@show loglik

        	loglik = new_loglik
			break
		end

        loglik = new_loglik

	end

	end_time = time()

	elapsed_time = end_time - start_time

	return lambda, forest_left, forest_right, forest_product, forest_isLeaf, loglik, elapsed_time
end

function tcm_forest_EM(N, K, forest_left, forest_right, forest_product, forest_isLeaf, assortments, transaction_counts, time_limit_EM, initial_lambda)
	lambda = copy(initial_lambda); #1/K * ones(K)
    v, v_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
    logv = log.(v)
    logv[ isinf.(logv)] = 0.0;
    loglik = dot( convert(Array{Float64,}, transaction_counts), logv)

    # @show v
    # @show loglik

	v = zeros(size(transaction_counts))

	perTol_EM = 1e-5
    iter = 0;

	start_time = time()

	while (time() - start_time < time_limit_EM)
		iter += 1;
		h, v, v_by_class = tcm_forest_E(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
		lambda = tcm_forest_M_update_lambda(K, h, transaction_counts)

		v, v_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
		logv = log.(v)
		logv[ isinf.(logv)] = 0.0;
		new_loglik = dot( convert(Array{Float64,}, transaction_counts), logv)

		println("\t Ranking EM: iter ", iter, " -- time elapsed = ", time() - start_time, " log likelihood = ", new_loglik)

	    if ( (new_loglik - loglik) / abs(loglik) < perTol_EM )
	    	temp = (new_loglik - loglik) / abs(loglik)
	    	# @show temp
	    	# @show new_loglik
	    	# @show loglik

	    	loglik = new_loglik
			break
		end

        # if (iter > 100)
        #     error("Stop here")
        # end

		loglik = new_loglik
	end

	elapsed_time = time() - start_time

	return lambda, v, loglik, elapsed_time

end

function tcm_forest_M_update_lambda(K, h, transaction_counts)
	sum_of_h = zeros(K)

	for k in 1:K
		sum_of_h[k] = dot( h[k], convert(Array{Float64}, transaction_counts))
	end

	lambda = sum_of_h ./ sum(sum_of_h)

	return lambda
end

function tcm_forest_E(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)
	v, v_by_class = tcm_predictForest(N, K, lambda, forest_left, forest_right, forest_product, forest_isLeaf, assortments)

	h = Array{ Array{Float64,2}, 1}(K)
	for k in 1:K
		single_h = lambda[k] * v_by_class[k] ./ v
		single_h[ isnan.(single_h)] = 0.0;

		h[k] = single_h
	end

	return h, v, v_by_class
end

function tcm_solveTreeSubproblem(N, u, assortments, depth_limit)
    # @objective(m_sub, Max, sum(u[m,p] * a[m,p] for m in 1:M, p in 1:N))
    # solve(m_sub)

    # x_val = getvalue(x)
    # objval = getobjectivevalue(m_sub)

    # single_ranking = sum(x_val, 2)
 #    single_ranking = single_ranking[:]
    # single_ordering = sortperm(single_ranking, rev = true)

    # return single_ordering, objval

    tol_RC = 1e-4

    M = size(assortments,1)

    single_left = Int64[0]
    single_right = Int64[0]
    single_product = Int64[N]
    single_isLeaf = Bool[true]
    single_depth = Int64[1]
    single_include_set = Array{Int64}[ [N] ];
    single_assortment_binary = Array{Int64}[ ones(M) ]
    single_reduced_cost = Float64[0.0];
    # single_reduced_cost = Float64[ sum(u[:,N]) - 2*tol_RC ] #
    single_isClosed = Bool[false]

    iter = 0 ;

    while (true)
        iter += 1
        current_reduced_cost = sum( single_reduced_cost[single_isLeaf])
        println("\t\t Subproblem Iter : ", iter, " -- reduced cost: ", current_reduced_cost, " -- num leaves: ", sum(single_isLeaf))
        candidate_leaves = find(single_isLeaf .& (single_depth .< depth_limit) .&  (.!single_isClosed) )

        if (isempty(candidate_leaves))
            println("\t\t\t No more eligible leaves to split; exiting procedure...")
            break;
        end

        num_candL = length(candidate_leaves)
        candidate_reduced_costs = -Inf* ones(num_candL)
        best_split_p_by_leaf = zeros(Int64, num_candL)
        best_left_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_right_leaf_p_by_leaf  = zeros(Int64, num_candL)
        best_left_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_right_reduced_cost_by_leaf  = zeros(Float64, num_candL)
        best_left_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)
        best_right_assortment_binary_by_leaf = Array{Array{Int64}}(num_candL)

        for ell_ind in 1:num_candL

            leaf = candidate_leaves[ell_ind]

            assortment_binary = copy(single_assortment_binary[leaf])

            candidate_split_products = setdiff(collect(1:(N-1)), single_include_set[leaf])

            for split_p in candidate_split_products

                candidate_left_leaf_products = union(single_include_set[leaf], split_p)
                best_left_reduced_cost = -Inf
                best_left_p = 0;
                best_left_assortment_binary = zeros(M)
                for left_p in candidate_left_leaf_products
                    left_assortment_binary = assortment_binary .* assortments[:, split_p][:]
                    temp = dot( left_assortment_binary, u[:, left_p])
                    if (temp > best_left_reduced_cost)
                        best_left_reduced_cost = temp
                        best_left_p = left_p
                        best_left_assortment_binary = copy(left_assortment_binary)
                    end
                end

                candidate_right_leaf_products = copy(single_include_set[leaf])
                best_right_reduced_cost = -Inf
                best_right_p = 0;
                best_right_assortment_binary = zeros(M)
                for right_p in candidate_right_leaf_products
                    right_assortment_binary = assortment_binary .* (1 - assortments[:,split_p][:])
                    temp = dot(right_assortment_binary, u[:, right_p])
                    if (temp > best_right_reduced_cost)
                        best_right_reduced_cost = temp
                        best_right_p = right_p
                        best_right_assortment_binary = copy(right_assortment_binary)
                    end
                end

                new_reduced_cost = best_left_reduced_cost + best_right_reduced_cost

                if (new_reduced_cost > candidate_reduced_costs[ell_ind] )
                    candidate_reduced_costs[ell_ind] = new_reduced_cost
                    best_split_p_by_leaf[ell_ind] = split_p
                    best_left_leaf_p_by_leaf[ell_ind] = best_left_p
                    best_right_leaf_p_by_leaf[ell_ind] = best_right_p
                    best_left_reduced_cost_by_leaf[ell_ind] = best_left_reduced_cost
                    best_right_reduced_cost_by_leaf[ell_ind] = best_right_reduced_cost
                    best_left_assortment_binary_by_leaf[ell_ind] = best_left_assortment_binary # VVM: should there be a copy here?
                    best_right_assortment_binary_by_leaf[ell_ind] = best_right_assortment_binary # VVM: should there be a copy here?
                end
            end
        end

        incumbent_reduced_costs = single_reduced_cost[candidate_leaves]

        tol_RC = 1e-4

        for ell_ind in 1:num_candL
            # @show incumbent_reduced_costs[ell_ind]
            # @show candidate_reduced_costs[ell_ind]

            if (incumbent_reduced_costs[ell_ind] + tol_RC < candidate_reduced_costs[ell_ind])
                # Splitting this leaf does improve things, so go ahead with

                leaf = candidate_leaves[ell_ind]

                # Grow tree.
                numNodes = length(single_left)
                push!(single_left, 0, 0)
                push!(single_right, 0, 0)
                push!(single_isLeaf, true, true)
                push!(single_product, best_left_leaf_p_by_leaf[ell_ind], best_right_leaf_p_by_leaf[ell_ind])
                push!(single_depth, single_depth[leaf]+1, single_depth[leaf]+1)
                push!(single_include_set, union(single_include_set[leaf],best_split_p_by_leaf[ell_ind]), copy(single_include_set[leaf]) )
                push!(single_assortment_binary, copy(best_left_assortment_binary_by_leaf[ell_ind]), copy(best_right_assortment_binary_by_leaf[ell_ind]))
                push!(single_reduced_cost, best_left_reduced_cost_by_leaf[ell_ind], best_right_reduced_cost_by_leaf[ell_ind])
                push!(single_isClosed, false, false)

                # Update the leaf
                single_left[leaf] = numNodes+1
                single_right[leaf] = numNodes+2
                single_product[leaf] = best_split_p_by_leaf[ell_ind]
                single_isClosed[leaf] = true
                single_isLeaf[leaf] = false
            else
                # Leaf does not improve
                leaf = candidate_leaves[ell_ind]
                single_isClosed[leaf] = true
            end
        end
    end

    optimal_reduced_cost = sum( single_reduced_cost[single_isLeaf] )

    return single_left, single_right, single_product, single_isLeaf, optimal_reduced_cost

end

function tcm_convertRankingToForest(N, orderings)
    K = size(orderings,1)

    forest_left = Array{ Array{Int64,1} }(K)
    forest_right = Array{ Array{Int64,1} }(K)
    forest_product = Array{ Array{Int64,1} }(K)
    forest_isLeaf = Array{ Array{Bool,1} }(K)

    for k in 1:K
        cn = 1;
        single_left = Int64[]
        single_right = Int64[]
        single_product = Int64[]
        single_isLeaf = Bool[]

        for i in 1:N
            if (orderings[k,i] == N)
                push!(single_left, 0)
                push!(single_right, 0)
                push!(single_product, N)
                push!(single_isLeaf, true)
                break;
            else
                # Grow tree by one leaf.
                push!(single_left, cn+1, 0 )
                push!(single_right, cn+2, 0 )
                push!(single_product, orderings[k,i], orderings[k,i])
                push!(single_isLeaf, false, true)
                cn += 2;
            end
        end

        forest_left[k] = single_left
        forest_right[k] = single_right
        forest_product[k] = single_product
        forest_isLeaf[k] = single_isLeaf
    end

    return forest_left, forest_right, forest_product, forest_isLeaf
end

function tcm_createBasicForest(N)
    K = N
    orderings = zeros(Int64, K, N)
    for k in 1:(N-1)
        single_ordering = [k; N; collect(1:(k-1)); collect( (k+1):(N-1) )]
        orderings[k,:] = single_ordering
    end
    orderings[N,:] = [N; collect(1:(N-1))]

    forest_left, forest_right, forest_product, forest_isLeaf = tcm_convertRankingToForest(N, orderings)
    return forest_left, forest_right, forest_product, forest_isLeaf
end

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_licenses: MIT Unlicense


text:
### wiener-time-series.jl
#
# Copyright (C) 2016 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: deconvolution, signal processing
#
# This file is a part of Deconvolution.jl.
#
# License is MIT "Expat".
#
### Commentary:
#
# This file produces an example image used in Deconvolution.jl documentation.
#
### Code:

srand(7)

using LombScargle, Deconvolution, Plots
t = linspace(0, 10, 1000) # observation times
x = sinpi(t) .* cos(5t) - 1.5cospi(t) .* sin(2t) # the original signal
n = rand(length(x)) # noise to be added
y = x + 3(n - mean(n)) # observed noisy signal
# Lomb-Scargle periodogram
p = lombscargle(t, y, maximum_frequency=2, samples_per_peak=10)
plot(freqpower(p)...)

m1 = LombScargle.model(t, y, findmaxfreq(p, [0, 0.5])[1]) # first model
m2 = LombScargle.model(t, y, findmaxfreq(p, [0.5, 1])[1]) # second model
m3 = LombScargle.model(t, y, findmaxfreq(p, [1, 1.5])[1]) # third model

signal = m1 + m2 + m3
noise = rand(length(y)) # noise for `wiener`
polished = wiener(y, signal, noise)
plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
plot!(t, y, label="Observed signal")
savefig("wiener-time-series-observed.png")
plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
plot!(t, polished, label="Recovered with Wiener")
plot!(t, signal, label="Lomb–Scargle model")
savefig("wiener-time-series-recovered.png")


text:
using Distributions
using Test
using DataFrames
using CSVFiles
using DataStructures
# using VegaLite

import GlobalSensitivityAnalysis: borehole
import StatsBase: ordinalrank

################################################################################
## SETUP
################################################################################

# define the (uncertain) parameters of the problem and their distributions
data = DeltaData(
    params = OrderedDict(
        :rw => Normal(0.10, 0.0161812),
        :r => LogNormal(7.71, 1.0056),
        :Tu => Uniform(63070, 115600),
        :Hu => Uniform(990, 1110),
        :Tl => Uniform(63.1, 116),
        :Hl => Uniform(700, 820),
        :L => Uniform(1120, 1680),
        :Kw => Uniform(9855, 12045)
    ),
    N = 1000
)

N = data.N
D = length(data.params)

@testset "Borehole Function - Sampling" begin

    py_samples = load("data/delta/py_borehole/py_samples.csv", header_exists=false, colnames = ["rw", "r", "Tu", "Hu", "Tl", "Hl", "L", "Kw"]) |> DataFrame
    julia_samples = sample(data)

    quants = [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]
    sig_level = 0.05
    output_dir = joinpath(@__DIR__, "../output/LHS_QuantileTesting/borehole")
    mkpath(output_dir)

    for i = 1:D
        sampleA = py_samples[:,i]
        sampleB = julia_samples[:,i]

        # Make dataframe
        # df = DataFrame(A = sampleA, B = sampleB)
        # save(joinpath(output_dir, "samples$(N)_p$i.csv"), df)

        # # Quick plot comparison
        # p1 = df |> @vlplot(:bar, x = {:A, bin = {step=0.1}}, y="count()")
        # save(joinpath(output_dir, "sampleA_$(N)_p$i.png"), p1)
        # p2 = df |> @vlplot(:bar, x = {:B, bin = {step=0.1}}, y="count()", background=:white)
        # save(joinpath(output_dir, "sampleB_$(N)_p$i.png"), p2)

        # Run WRS quantile matching
        results = WRS.pb2gen(sampleA, sampleB, quantiles = quants) |> DataFrame
        # save(joinpath(output_dir, "LHS Sampling Quantile Comparison N$(N)_p$i.csv"), results)

        @test sum(results[!,:signif]) == 0
    end

end

@testset "Borehole Function - Analysis" begin

    # Get the samples
    samples = load("data/delta/py_borehole/py_samples.csv", header_exists=false, colnames = ["rw", "r", "Tu", "Hu", "Tl", "Hl", "L", "Kw"]) |> DataFrame
    samples = Matrix(samples)

    # check borehole
    python_Y = load("data/delta/py_borehole/py_borehole.csv", header_exists=false) |> DataFrame
    julia_Y = borehole(samples)
    @test python_Y[:,1] ≈ julia_Y atol = ATOL_sample

    # julia results
    julia_results = analyze(data, samples, julia_Y; num_resamples = 1_000)

    # python results
    py_firstorder = load("data/delta/py_borehole/py_firstorder.csv", header_exists=false) |> DataFrame
    py_delta = load("data/delta/py_borehole/py_delta.csv", header_exists=false) |> DataFrame
    py_firstorder_conf = load("data/delta/py_borehole/py_firstorder_conf.csv", header_exists=false) |> DataFrame
    py_delta_conf = load("data/delta/py_borehole/py_delta_conf.csv", header_exists=false) |> DataFrame

    # test indices
    @test julia_results[:firstorder] ≈ Matrix(py_firstorder) atol = ATOL_delta
    @test ordinalrank(julia_results[:firstorder]) == ordinalrank(py_firstorder[!,:Column1])

    @test julia_results[:delta] ≈ Matrix(py_delta) atol = 0.05 # TODO - this seems too high?
    @test ordinalrank(julia_results[:delta]) == ordinalrank(py_delta[!,:Column1])

    # test confidence intervals
    @test julia_results[:firstorder_conf] ≈ Matrix(py_firstorder_conf) atol = ATOL_CI
    @test julia_results[:delta_conf] ≈ Matrix(py_delta_conf) atol = ATOL_CI

end

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
module FeatureExtraction

# package code goes here

importall MLKernels

export lda_components!,
       lda_matrices!,
       lda!,
       transform_lda,
       klda!,
       transform_klda,
       pca,
       transform_pca,
       kpca,
       transform_kpca

include("common.jl")
include("kpca.jl")
include("klda.jl")

end # module


text:
function _test_matrix_equal(A::SparseMatrixCSC, B::SparseMatrixCSC)
    @test A.m == B.m
    @test A.n == B.n
    @test A.nzval ≈ B.nzval atol=ATOL rtol=RTOL
    @test A.rowval == B.rowval
    @test A.colptr == B.colptr
end
function _test_matrix_equal(A::MatOI.SparseMatrixCSRtoCSC{Tv, Ti, I}, B::SparseMatrixCSC) where {Tv, Ti, I}
    @test A.m == B.m
    @test A.n == B.n
    @test A.nzval ≈ B.nzval atol=ATOL rtol=RTOL
    if I <: MatOI.OneBasedIndexing
        @test A.rowval == B.rowval
        @test A.colptr == B.colptr
    else
        @test A.rowval == B.rowval .- 1
        @test A.colptr == B.colptr .- 1
    end
    sA = convert(typeof(B), A)
    @test typeof(sA) == typeof(B)
    _test_matrix_equal(sA, B)
end

# _psd1test: https://github.com/jump-dev/MathOptInterface.jl/blob/master/src/Test/contconic.jl#L2417
function psd1(::Type{T}, ::Type{I}) where {T, I}
    # We use `MockOptimizer` to have indices xor'ed so that it tests that we don't assumes they are `1:n`.
    model = MOIU.MockOptimizer(MOIU.Model{T}())

    X = MOI.add_variables(model, 6)
    x = MOI.add_variables(model, 3)

    vov = MOI.VectorOfVariables(X)

    cX = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{T}(vov), MOI.PositiveSemidefiniteConeTriangle(3)
    )

    cx = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction{T}(MOI.VectorOfVariables(x)), MOI.SecondOrderCone(3)
    )

    c1 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.(ones(T, 4), [X[1], X[3], X[end], x[1]])),
            [-one(T)]
        ),
        MOI.Zeros(1)
    )

    c2 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1:1, MOI.ScalarAffineTerm.(T[1, 2, 1, 2, 2, 1, 1, 1], [X; x[2]; x[3]])),
            [-inv(T(2))]
        ),
        MOI.Zeros(1)
    )

    objXidx = [1:3; 5:6]
    objXcoefs = 2ones(T, 5)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(),
        MOI.ScalarAffineFunction(
            MOI.ScalarAffineTerm.([objXcoefs; one(T)], [X[objXidx]; x[1]]),
            zero(T)
        )
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MIN_SENSE)

    conic_form = MatOI.GeometricConicForm{T, MatOI.SparseMatrixCSRtoCSC{T, Int, I}, Vector{T}}([MOI.PositiveSemidefiniteConeTriangle, MOI.SecondOrderCone, MOI.Zeros])
    index_map = MOI.copy_to(conic_form, model)

    @test conic_form.c' ≈ T[2 2 2 0 2 2 1 0 0]
    @test conic_form.b' ≈ T[0 0 0 0 0 0 0 0 0 -1 -inv(T(2))]
    _test_matrix_equal(
        conic_form.A,
        SparseMatrixCSC(
            11, 9,
            [1, 4, 6, 9, 11, 13, 16, 18, 20, 22],
            [1, 10, 11, 2, 11, 3, 10, 11, 4, 11, 5, 11, 6, 10, 11, 7, 10, 8, 11, 9, 11],
            T[-1, -1, -1, -1, -2, -1, -1, -1, -1, -2, -1, -2, -1, -1, -1, -1, -1, -1, -1, -1, -1],
        )
    )
end

# Taken from `MOI.Test.psdt2test`.
# find equivalent diffcp program here - https://github.com/AKS1996/jump-gsoc-2020/blob/master/diffcp_sdp_3_py.ipynb
function psd2(::Type{T}, ::Type{I}, η::T = T(10), α::T = T(4)/T(5), δ::T = T(9)/T(10)) where {T, I}
    # We use `MockOptimizer` to have indices xor'ed so that it tests that we don't assumes they are `1:n`.
    model = MOIU.MockOptimizer(MOIU.Model{T}())

    x = MOI.add_variables(model, 7)

    c1  = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(-one(T), x[1:6])),
            [η]
        ),
        MOI.Nonnegatives(1)
    )
    c2 = MOI.add_constraint(model, MOI.VectorAffineFunction(MOI.VectorAffineTerm.(1:6, MOI.ScalarAffineTerm.(one(T), x[1:6])), zeros(T, 6)), MOI.Nonnegatives(6))

    c3 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(
                [fill(1, 7); fill(2, 5);     fill(3, 6)],
                MOI.ScalarAffineTerm.(
                [ δ/2,       α,   δ, δ/4, δ/8,      0, -1,
                 -δ/(2*√2), -δ/4, 0,     -δ/(8*√2), 0,
                  δ/2,     δ-α,   0,      δ/8,      δ/4, -1],
                [x[1:7];     x[1:3]; x[5:6]; x[1:3]; x[5:7]])),
                zeros(T, 3)
        ),
        MOI.PositiveSemidefiniteConeTriangle(2)
    )
    c4 = MOI.add_constraint(
        model,
        MOI.VectorAffineFunction(
            MOI.VectorAffineTerm.(1, MOI.ScalarAffineTerm.(zero(T), [x[1:3]; x[5:6]])),
            [zero(T)]
        ),
        MOI.Zeros(1)
    )

    MOI.set(model, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{T}}(), MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(one(T), x[7])], zero(T)))
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)

    conic_form = MatOI.GeometricConicForm{T, MatOI.SparseMatrixCSRtoCSC{T, Int, I}, Vector{T}}([MOI.Nonnegatives, MOI.Zeros, MOI.PositiveSemidefiniteConeTriangle])
    index_map = MOI.copy_to(conic_form, model)

    @test conic_form.c ≈ [zeros(T, 6); one(T)]
    @test conic_form.b ≈ [T(10); zeros(T, 10)]
    _test_matrix_equal(
        conic_form.A,
        SparseMatrixCSC(
            11, 7,
            [1, 6, 11, 14, 17, 22, 25, 27],
            [1, 2, 9, 10, 11, 1, 3, 9, 10, 11, 1, 4, 9, 1, 5, 9, 1, 6, 9, 10, 11, 1, 7, 11, 9, 11],
            T[1.0, -1.0, -0.45, 0.318198, -0.45, 1.0, -1.0, -0.8, 0.225, -0.1, 1.0, -1.0, -0.9, 1.0, -1.0, -0.225, 1.0, -1.0, -0.1125, 0.0795495, -0.1125, 1.0, -1.0, -0.225, 1.0, 1.0],
        )
    )
end

@testset "PSD $T, $I" for T in [Float64, BigFloat], I in [MatOI.ZeroBasedIndexing, MatOI.OneBasedIndexing]
    psd1(T, I)
    psd2(T, I)
end


text:
#
# Subgradient descent types and methods.
#
"""
    SubgradientMethod

The parameters for solving a Lagrangian dual using subgradients.

initialbound   Initial upper or lower bound for the Lagrangian dual
tol            Tolerance for terminating
wait           A parameter for reducing step sizes
maxit          A cap on iterations
"""
mutable struct SubgradientMethod{T<:Tolerance} <: AbstractLagrangianMethod
    initialbound::Float64
    tol::T
    wait::Int
    maxit::Int
end
function SubgradientMethod(; initialbound=0.0, tol=Unit(1e-6), wait=30, maxit=10_000)
    SubgradientMethod(initialbound, tol, wait, maxit)
end

function correctsign!(π::Array{Float64}, senses::Array{Symbol})
    for (i, sense) in enumerate(senses)
        if sense == :le
            π[i] = abs(π[i])
        elseif sense == :ge
            π[i] = -abs(π[i])
        end
    end
end

"""
    lagrangian_method!(lp::LinearProgramData{SubgradientMethod}, m::JuMP.Model, π::Vector{Float64})

Solve the Lagrangian dual of a linear program using subgradient descent.

# Arguments
* lp        Information about the primal problem
* m         The primal problem
* π         Initial iterate

# Returns
* status, objective, and modifies π
"""
function lagrangian_method!(lp::LinearProgramData{SubgradientMethod{T}}, m::JuMP.Model, π::Vector{Float64}) where T

    subgradient = lp.method
    tol = subgradient.tol

    # Stepping parameter
    μ = 2.0
    incumbent = subgradient.initialbound
    # A bound, the strongest bound so far, and a bound we can check against every lp.wait iterations
    bound = best_bound = cached_bound = -Inf

    iteration = 0
    timesunchanged = 0

    reversed = false
    if getobjectivesense(m) == :Max
        incumbent *= -1
        lp.obj *= -1
        setobjectivesense(m, :Min)
        π .= -π
        reversed = true
    end

    correctsign!(π, lp.senses) # if user started with the wrong signs

    while iteration < subgradient.maxit

        iteration += 1

        # Wait for some number of iterations before checking if we improved
        if mod(iteration, subgradient.wait) == 0
            # If we found a stronger lower bound
            if cached_bound < best_bound
                cached_bound = best_bound
            else
                # If we didn't improve, take smaller steps
                μ /= 2
                timesunchanged += 1
            end
        end

        # Let L = f + πᵀ(Ax-b), set this as the objective in the original model
        m.internalModelLoaded = false
        bound, direction = solve_primal(m, lp, π)
        # Update if we found a stronger bound
        (bound > best_bound) && (best_bound = bound)

        # Check if we have a new incumbent
        feasible = true
        for (slack, sense) in zip(direction, lp.senses)
            if !issatisfied(slack, sense)
                feasible = false
                break
            end
        end

        if feasible
            # The value of cᵀx, the primal objective, is an upper bound
            candidate = getvalue(lp.obj)
            # Check πᵀ(Ax-b) = candidate - bound
            if isclose(candidate, bound, tol) || timesunchanged > 30 # need to think about how to terminate with duality gap
                # Complementary slackness => we are done
                if reversed
                    # Undo transformation
                    π .*= -1
                    lp.obj *= -1
                    setobjectivesense(m, :Max)
                    return :Optimal, -bound
                else
                    return :Optimal, bound
                end
            end
            # Update incumbent if we strictly improved
            (candidate < incumbent) && (incumbent = candidate)
        end

        # Calculate the stepsize
        step = μ * (incumbent - bound) / sum(direction.^2)
        # Update multipliers
        for i = eachindex(π)
            π[i] += step * direction[i]
        end
        # Set correct sign for all multipliers
        for (i, s) in enumerate(lp.senses)
            if s == :le
                π[i] = max(0, π[i])
            elseif s == :ge
                π[i] = min(0, π[i])
            end
        end
    end
    warn("Lagrangian relaxation did not solve properly.")
    return :IterationLimit, bound
end


text:
function test_recall(truth::Array{Int,1}, pred::Array{Int,1}, k::Int)
    actual = measure(Recall(), truth, pred, k)
    expected = 0.333
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_precision(truth::Array{Int,1}, pred::Array{Int,1}, k::Int)
    actual = measure(Precision(), truth, pred, k)
    expected = 0.5
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_map(truth::Array{Int,1}, pred::Array{Int,1})
    actual = measure(MAP(), truth, pred)
    expected = 0.756
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_auc(truth::Array{Int,1}, pred::Array{Int,1})
    actual = measure(AUC(), truth, pred)
    expected = 0.667
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_mrr(truth::Array{Int,1}, pred::Array{Int,1})
    actual = measure(ReciprocalRank(), truth, pred)
    expected = 1.0
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_mpr(truth::Array{Int,1}, pred::Array{Int,1})
    actual = measure(MPR(), truth, pred)
    expected = 33.333
    eps = 0.001

    @test abs(actual - expected) < eps
end

function test_ndcg(truth::Array{Int,1}, pred::Array{Int,1}, k::Int)
    actual = measure(NDCG(), truth, pred, k)
    expected = 0.613
    eps = 0.001

    @test abs(actual - expected) < eps
end

println("-- Testing ranking metrics")

truth = [1, 2, 4]
pred = [1, 3, 2, 6, 4, 5]
k = 2

test_recall(truth, pred, k)
test_precision(truth, pred, k)
test_map(truth, pred)
test_auc(truth, pred)
test_mrr(truth, pred)
test_mpr(truth, pred)
test_ndcg(truth, pred, k)


max_stars_repo_name: AnAverageHuman/competitive

max_issues_repo_name: AnAverageHuman/competitive

max_forks_repo_name: AnAverageHuman/competitive


text:
##### Multi-physics types

using CLIMAParameters.Planet: T_freeze, cv_l

using ..Microphysics_0M
using ..Microphysics

export Subsidence
struct Subsidence{PV, FT} <: TendencyDef{Source, PV}
    D::FT
end

# Subsidence includes tendencies in Mass, Energy and TotalMoisture equations:
Subsidence(D::FT) where {FT} = (
    Subsidence{Mass, FT}(D),
    Subsidence{Energy, FT}(D),
    Subsidence{TotalMoisture, FT}(D),
)

subsidence_velocity(subsidence::Subsidence{PV, FT}, z::FT) where {PV, FT} =
    -subsidence.D * z

struct PressureGradient{PV <: Momentum} <: TendencyDef{Flux{FirstOrder}, PV} end
struct Pressure{PV <: Energy} <: TendencyDef{Flux{FirstOrder}, PV} end

struct Advect{PV} <: TendencyDef{Flux{FirstOrder}, PV} end
struct Diffusion{PV} <: TendencyDef{Flux{SecondOrder}, PV} end

export RemovePrecipitation
"""
    RemovePrecipitation{PV} <: TendencyDef{Source, PV}
A sink to `q_tot` when cloud condensate is exceeding a threshold.
The threshold is defined either in terms of condensate or supersaturation.
The removal rate is implemented as a relaxation term
in the Microphysics_0M module.
The default thresholds and timescale are defined in CLIMAParameters.jl.
"""
struct RemovePrecipitation{PV <: Union{Mass, Energy, TotalMoisture}} <:
       TendencyDef{Source, PV}
    " Set to true if using qc based threshold"
    use_qc_thr::Bool
end

RemovePrecipitation(b::Bool) = (
    RemovePrecipitation{Mass}(b),
    RemovePrecipitation{Energy}(b),
    RemovePrecipitation{TotalMoisture}(b),
)

function remove_precipitation_sources(
    s::RemovePrecipitation{PV},
    atmos,
    args,
) where {PV <: Union{Mass, Energy, TotalMoisture}}
    @unpack state, aux = args
    @unpack ts = args.precomputed

    FT = eltype(state)

    q = PhasePartition(ts)
    λ::FT = liquid_fraction(ts)
    I_l::FT = internal_energy_liquid(ts)
    I_i::FT = internal_energy_ice(ts)
    Φ::FT = gravitational_potential(atmos.orientation, aux)

    S_qt::FT = 0
    if s.use_qc_thr
        S_qt = remove_precipitation(atmos.param_set, q)
    else
        q_vap_sat = q_vap_saturation(ts)
        S_qt = remove_precipitation(atmos.param_set, q, q_vap_sat)
    end

    S_e::FT = (λ * I_l + (1 - λ) * I_i + Φ) * S_qt

    return (S_ρ_qt = state.ρ * S_qt, S_ρ_e = state.ρ * S_e)
end

export WarmRain_1M
"""
    WarmRain_1M{FT} <: AbstractSource
A collection of source/sink terms related to 1-moment warm rain microphysics.
"""
struct WarmRain_1M{
    PV <: Union{Mass, Energy, TotalMoisture, LiquidMoisture, Rain},
} <: TendencyDef{Source, PV} end

WarmRain_1M() = (
    WarmRain_1M{Mass}(),
    WarmRain_1M{Energy}(),
    WarmRain_1M{TotalMoisture}(),
    WarmRain_1M{LiquidMoisture}(),
    WarmRain_1M{Rain}(),
)

function warm_rain_sources(atmos, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed

    FT = eltype(state)

    q = PhasePartition(ts)
    T::FT = air_temperature(ts)
    I_l::FT = internal_energy_liquid(ts)
    q_rai::FT = state.precipitation.ρq_rai / state.ρ
    Φ::FT = gravitational_potential(atmos.orientation, aux)

    # autoconversion
    src_q_rai_acnv = conv_q_liq_to_q_rai(atmos.param_set.microphys.rai, q.liq)
    # accretion
    src_q_rai_accr = accretion(
        atmos.param_set,
        atmos.param_set.microphys.liq,
        atmos.param_set.microphys.rai,
        q.liq,
        q_rai,
        state.ρ,
    )
    # rain evaporation
    src_q_rai_evap = evaporation_sublimation(
        atmos.param_set,
        atmos.param_set.microphys.rai,
        q,
        q_rai,
        state.ρ,
        T,
    )

    S_qr::FT = src_q_rai_acnv + src_q_rai_accr + src_q_rai_evap
    S_ql::FT = -src_q_rai_acnv - src_q_rai_accr
    S_qt::FT = -S_qr
    S_e::FT = S_qt * (I_l + Φ)

    return (
        S_ρ_qt = state.ρ * S_qt,
        S_ρ_ql = state.ρ * S_ql,
        S_ρ_qr = state.ρ * S_qr,
        S_ρ_e = state.ρ * S_e,
    )
end

export RainSnow_1M
"""
    RainSnow_1M{FT} <: AbstractSource
A collection of source/sink terms related to 1-moment rain and snow microphysics
"""
struct RainSnow_1M{PV <: Union{Mass, Energy, Moisture, Rain, Snow}} <:
       TendencyDef{Source, PV} end

RainSnow_1M() = (
    RainSnow_1M{Mass}(),
    RainSnow_1M{Energy}(),
    RainSnow_1M{TotalMoisture}(),
    RainSnow_1M{LiquidMoisture}(),
    RainSnow_1M{IceMoisture}(),
    RainSnow_1M{Rain}(),
    RainSnow_1M{Snow}(),
)

function rain_snow_sources(atmos, args)
    @unpack state, aux = args
    @unpack ts = args.precomputed

    FT = eltype(state)

    q_rai::FT = state.precipitation.ρq_rai / state.ρ
    q_sno::FT = state.precipitation.ρq_sno / state.ρ
    q = PhasePartition(ts)
    T::FT = air_temperature(ts)
    I_d::FT = internal_energy_dry(ts)
    I_v::FT = internal_energy_vapor(ts)
    I_l::FT = internal_energy_liquid(ts)
    I_i::FT = internal_energy_ice(ts)
    _T_freeze::FT = T_freeze(atmos.param_set)
    _L_f::FT = latent_heat_fusion(ts)
    _cv_l::FT = cv_l(atmos.param_set)
    Φ::FT = gravitational_potential(atmos.orientation, aux)

    # temporary vars for summming different source terms
    S_qr::FT = FT(0)
    S_qs::FT = FT(0)
    S_ql::FT = FT(0)
    S_qi::FT = FT(0)
    S_qt::FT = FT(0)
    S_e::FT = FT(0)

    # source of rain via autoconversion
    tmp = conv_q_liq_to_q_rai(atmos.param_set.microphys.rai, q.liq)
    S_qr += tmp
    S_ql -= tmp
    S_e -= tmp * (I_l + Φ)

    # source of snow via autoconversion
    tmp = conv_q_ice_to_q_sno(
        atmos.param_set,
        atmos.param_set.microphys.ice,
        q,
        state.ρ,
        T,
    )
    S_qs += tmp
    S_qi -= tmp
    S_e -= tmp * (I_i + Φ)

    # source of rain water via accretion cloud water - rain
    tmp = accretion(
        atmos.param_set,
        atmos.param_set.microphys.liq,
        atmos.param_set.microphys.rai,
        q.liq,
        q_rai,
        state.ρ,
    )
    S_qr += tmp
    S_ql -= tmp
    S_e -= tmp * (I_l + Φ)

    # source of snow via accretion cloud ice - snow
    tmp = accretion(
        atmos.param_set,
        atmos.param_set.microphys.ice,
        atmos.param_set.microphys.sno,
        q.ice,
        q_sno,
        state.ρ,
    )
    S_qs += tmp
    S_qi -= tmp
    S_e -= tmp * (I_i + Φ)

    # sink of cloud water via accretion cloud water - snow
    tmp = accretion(
        atmos.param_set,
        atmos.param_set.microphys.liq,
        atmos.param_set.microphys.sno,
        q.liq,
        q_sno,
        state.ρ,
    )
    if T < _T_freeze # cloud droplets freeze to become snow)
        S_qs += tmp
        S_ql -= tmp
        S_e -= tmp * (I_i + Φ)
    else # snow melts, both cloud water and snow become rain
        α::FT = _cv_l / _L_f * (T - _T_freeze)
        S_ql -= tmp
        S_qs -= tmp * α
        S_qr += tmp * (1 + α)
        S_e -= tmp * ((1 + α) * I_l - α * I_i + Φ)
    end

    # sink of cloud ice via accretion cloud ice - rain
    tmp1 = accretion(
        atmos.param_set,
        atmos.param_set.microphys.ice,
        atmos.param_set.microphys.rai,
        q.ice,
        q_rai,
        state.ρ,
    )
    # sink of rain via accretion cloud ice - rain
    tmp2 = accretion_rain_sink(
        atmos.param_set,
        atmos.param_set.microphys.ice,
        atmos.param_set.microphys.rai,
        q.ice,
        q_rai,
        state.ρ,
    )
    S_qi -= tmp1
    S_e -= tmp1 * (I_i + Φ)
    S_qr -= tmp2
    S_e += tmp2 * _L_f
    S_qs += tmp1 + tmp2

    # accretion rain - snow
    if T < _T_freeze
        tmp = accretion_snow_rain(
            atmos.param_set,
            atmos.param_set.microphys.sno,
            atmos.param_set.microphys.rai,
            q_sno,
            q_rai,
            state.ρ,
        )
        S_qs += tmp
        S_qr -= tmp
        S_e += tmp * _L_f
    else
        tmp = accretion_snow_rain(
            atmos.param_set,
            atmos.param_set.microphys.rai,
            atmos.param_set.microphys.sno,
            q_rai,
            q_sno,
            state.ρ,
        )
        S_qs -= tmp
        S_qr += tmp
        S_e -= tmp * _L_f
    end

    # rain evaporation sink (it already has negative sign for evaporation)
    tmp = evaporation_sublimation(
        atmos.param_set,
        atmos.param_set.microphys.rai,
        q,
        q_rai,
        state.ρ,
        T,
    )
    S_qr += tmp
    S_e -= tmp * (I_l + Φ)

    # snow sublimation/deposition source/sink
    tmp = evaporation_sublimation(
        atmos.param_set,
        atmos.param_set.microphys.sno,
        q,
        q_sno,
        state.ρ,
        T,
    )
    S_qs += tmp
    S_e -= tmp * (I_i + Φ)

    # snow melt
    tmp = snow_melt(
        atmos.param_set,
        atmos.param_set.microphys.sno,
        q_sno,
        state.ρ,
        T,
    )
    S_qs -= tmp
    S_qr += tmp
    S_e -= tmp * _L_f

    # total qt sink is the sum of precip sources
    S_qt = -S_qr - S_qs

    return (
        S_ρ_qt = state.ρ * S_qt,
        S_ρ_ql = state.ρ * S_ql,
        S_ρ_qi = state.ρ * S_qi,
        S_ρ_qr = state.ρ * S_qr,
        S_ρ_qs = state.ρ * S_qs,
        S_ρ_e = state.ρ * S_e,
    )
end


text:
using Documenter, StateMachines

makedocs(
    modules = [StateMachines],
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    sitename = "StateMachines.jl",
    authors  = "Paul Bayer",
    pages = [
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Usage" => "usage.md",
        "Examples" => "examples.md",
        "Internals" => "internals.md"
    ]
)

deploydocs(
    repo   = "github.com/pbayer/StateMachines.jl.git",
    target = "build",
    deps   = nothing,
    make   = nothing
)


text:
using Aqua
using Configurations
Aqua.test_all(Configurations)

include("test_utils.jl")
include("test_basic.jl")
include("test_error.jl")
include("test_from_dict.jl")
include("test_parametric.jl")
include("test_reflection.jl")
include("test_issue50.jl")
include("test_reflected.jl")
include("test_auto_maybe_default.jl")
include("test_partial.jl")


text:
#=
see description in examples/polymin/native.jl
=#

include(joinpath(@__DIR__, "data_real.jl"))

struct PolyMinJuMP{T <: Real} <: ExampleInstanceJuMP{T}
    interp_vals::Vector{T}
    Ps::Vector{Matrix{T}}
    true_min::Real
    use_primal::Bool # solve primal, else solve dual
    use_wsos::Bool # use wsosinterpnonnegative cone, else PSD formulation
end

function PolyMinJuMP{Float64}(
    poly_name::Symbol,
    halfdeg::Int,
    args...)
    return PolyMinJuMP{Float64}(get_interp_data(Float64, poly_name, halfdeg)..., args...)
end

function PolyMinJuMP{Float64}(
    n::Int,
    halfdeg::Int,
    args...)
    return PolyMinJuMP{Float64}(random_interp_data(Float64, n, halfdeg)..., args...)
end

function build(inst::PolyMinJuMP{T}) where {T <: Float64}
    (interp_vals, Ps, use_primal) = (inst.interp_vals, inst.Ps, inst.use_primal)
    U = length(interp_vals)

    model = JuMP.Model()
    if use_primal
        JuMP.@variable(model, a)
        JuMP.@objective(model, Max, a)
    else
        JuMP.@variable(model, μ[1:U])
        JuMP.@objective(model, Min, dot(μ, interp_vals))
        JuMP.@constraint(model, sum(μ) == 1)
    end

    if inst.use_wsos
        cone = Hypatia.WSOSInterpNonnegativeCone{Float64, Float64}(U, Ps, !use_primal)
        aff_expr = (use_primal ? interp_vals .- a : μ)
        JuMP.@constraint(model, aff_expr in cone)
    else
        if use_primal
            psd_vars = []
            for Pr in Ps
                Lr = size(Pr, 2)
                psd_r = JuMP.@variable(model, [1:Lr, 1:Lr], Symmetric)
                if Lr == 1
                    # Mosek cannot handle 1x1 PSD constraints
                    JuMP.@constraint(model, psd_r[1, 1] >= 0)
                else
                    JuMP.@constraint(model, psd_r in JuMP.PSDCone())
                end
                push!(psd_vars, psd_r)
            end
            coeffs_lhs = JuMP.@expression(model, [u in 1:U], sum(sum(Pr[u, k] * Pr[u, l] * psd_r[k, l] * (k == l ? 1 : 2) for k in 1:size(Pr, 2) for l in 1:k) for (Pr, psd_r) in zip(Ps, psd_vars)))
            JuMP.@constraint(model, coeffs_lhs .== interp_vals .- a)
        else
            for Pr in Ps
                Lr = size(Pr, 2)
                psd_r = [JuMP.@expression(model, sum(Pr[u, i] * Pr[u, j] * μ[u] for u in 1:U)) for i in 1:Lr for j in 1:i]
                if Lr == 1
                    # Mosek cannot handle 1x1 PSD constraints
                    JuMP.@constraint(model, psd_r[1] >= 0)
                else
                    JuMP.@constraint(model, psd_r in MOI.PositiveSemidefiniteConeTriangle(Lr))
                end
            end
        end
    end

    return model
end

function test_extra(inst::PolyMinJuMP{T}, model::JuMP.Model) where T
    stat = JuMP.termination_status(model)
    @test stat == MOI.OPTIMAL
    if (stat == MOI.OPTIMAL) && !isnan(inst.true_min)
        # check objective value is correct
        tol = eps(T)^0.1
        @test JuMP.objective_value(model) ≈ inst.true_min atol=tol rtol=tol
    end
end


text:
export @cuconvertarrays, @createcuarrays

macro cuconvertarrays(vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = CuArray($(esc(var))); ) for var in vars])
  expr
end

macro createcuarrays(T, dims, vars...)
  expr = Expr(:block)
  append!(expr.args, [:($(esc(var)) = CuArray(zeros($(esc(T)), $(esc(dims)))); ) for var in vars])
  expr
end


text:
"""
    Token(form, head=-1, label=nothing)

A token in a dependency tree.
"""
struct Token{F,H<:Union{Int,Set{Int}},L}
    form::F
    head::H
    label::L
    data # Union{Nothing,Dict} ?
end

Token() = Token(nothing, -1, nothing, nothing)
function Token(form, head::Int=-1, label=nothing; kwargs...)
    return Token(form, head, label, isempty(kwargs) ? nothing : Dict(kwargs))
end
Token(form, head, label=nothing; kwargs...) =
    Token(form, Set(head), label, isempty(kwargs) ? nothing : Dict(kwargs))
Token(token::Token; head=token.head, label=token.label, kwargs...) =
    Token(token.form, head, label, isempty(kwargs) ? nothing : Dict(kwargs))

const Token1H{F,L} = Token{F,Int,L}
const TokenNH{F,L}  = Token{F,Set{Int},L}

const ROOT = Token("ROOT", 0, id=0)

"""
    deptoken(form, head=-1, label=nothing)

Create a token in a dependency tree.
"""
deptoken(form=nothing, head=-1, label=nothing; kwargs...) =
    Token(form, head, label; kwargs...)

function from_indices(x; form::Int=1, head=2, label=3, kw...)
    ks, is = isempty(kw) ? ((), ()) : zip(pairs(kw)...)
    new_kw = NamedTuple{ks}(Tuple(x[i] for i in is))
    return Token(x[form], x[head], x[label], new_kw...)
end

isroot(t::Token) = t === ROOT

has_head(token::Token1H) = token.head >= 0
has_head(token::TokenNH) = !isempty(token.head)

has_head(token::Token1H, h::Int) = token.head == h
has_head(token::TokenNH, h::Int) = h in token.head

function Base.getproperty(t::Token, f::Symbol)
    try
        getfield(t, f)
    catch
        try
            t.data[f]
        catch
            error("$t has no field $f")
        end
    end
end

==(a::Token, b::Token) =
    a.form == b.form && a.head == b.head && a.label == b.label && a.data == b.data

headisroot(t::Token1H) = t.head == 0
headisroot(t::TokenNH) = 0 in t.head


text:
type DiscretePolicy{VT, QT, PT} <: Policy

    V::DenseArray{VT} # utility function
    Q::DenseArray{QT} # Q-matrix
    policy::DenseArray{PT} # policy computed from the utility
end

function DiscretePolicy(;V::DenseArray=[], Q::DenseArray=[], P::DenseArray=[])

    vt = None
    qt = None
    pt = None

    !isempty(V) ? (vt = typeof(V[1])) : (nothing)
    !isempty(Q) ? (qt = typeof(Q[1])) : (nothing)
    !isempty(P) ? (pt = typeof(P[1])) : (nothing)

    return DiscretePolicy{vt, qt, pt}(V, Q, P)
end

function action(p::DiscretePolicy, s::Int64)

    policy = p.policy
    if isempty(policy)
        error("$(typeof(p)) has an empty policy array")
    end

    return policy[s]
end

function value(p::DiscretePolicy, s::Int64, a::Int64)

    Q = p.Q
    if isempty(Q)
        error("$(typeof(p)) has an empty Q-matrix")
    end

    return Q[a,s]
end

function value(p::DiscretePolicy, s::Int64)

    V = p.V
    if isempty(V)
        error("$(typeof(p)) has an empty optimal value array")
    end

    return V[s]
end

# computes policy given value function (can be done in parallel)
function computePolicy(mdp::DiscreteMDP, V::DenseArray)
    nStates  = numStates(mdp)
    nActions = numActions(mdp)

    @assert nStates == length(V)

    policy = zeros(Int64, nStates)

    for si = 1:nStates
        qHi = 0.0
        for ai = 1:nActions
            states, probs = nextStates(mdp, si, ai)
            qNow = reward(mdp, si, ai)
            for sp = 1:length(states)
                spi = states[sp]
                qNow += probs[sp] * V[spi]
            end
            if ai == 1 || qNow > qHi
                qHi = qNow
                policy[si] = ai
            end
        end
    end
    return policy
end


text:
function load()
end

max_stars_repo_name: ucatbas/seqnet-julia

max_issues_repo_name: ucatbas/seqnet-julia

max_forks_repo_name: ucatbas/seqnet-julia


text:
using NeXLParticle
using Test
using DataDeps

@testset "ZepTest" begin
    zeptest = datadep"ZepTestArtifact"
    zep = Zeppelin(joinpath(zeptest,"test.hdz"))
    @test eachparticle(zep)==1:250
    hdz = header(zep)
    @test hdz["MAG0"]=="347 78 250 98.24 5.112"
    @test hdz["INSTRUMENT"]=="NIST's TESCAN MIRA3 in 217 F101"
    es = elms(zep)
    @test length(es)==32
    @test n"Sb" in es
    @test beamenergy(zep)==20.0e3
    @test probecurrent(zep)==0.961832
    @test magdata(zep,0)["Fields"]==78.0
    @test magdata(zep,0)["Area"]==5.112
    @test length(classes(zep))==27
    @test maxparticle(zep,1:10)[123] == 42.0
    refs = references( [
        reference(n"Ag", joinpath(zeptest,"Standards","Ag std.msa")),
        reference(n"Al", joinpath(zeptest,"Standards","Al std.msa")),
        reference(n"Ba", joinpath(zeptest,"Standards","Ba std.msa")),
        reference(n"Bi", joinpath(zeptest,"Standards","Bi std.msa")),
        reference(n"Br", joinpath(zeptest,"Standards","Br std.msa")),
        reference(n"C", joinpath(zeptest,"Standards","C std.msa")),
        reference(n"Ca", joinpath(zeptest,"Standards","Ca std.msa")),
        reference(n"Ce", joinpath(zeptest,"Standards","Ce std.msa")),
        reference(n"Cl", joinpath(zeptest,"Standards","Cl std.msa")),
        reference(n"Co", joinpath(zeptest,"Standards","Co std.msa")),
        reference(n"Cr", joinpath(zeptest,"Standards","Cr std.msa")),
        reference(n"Cu", joinpath(zeptest,"Standards","Cu std.msa")),
        reference(n"F", joinpath(zeptest,"Standards","F std.msa")),
        reference(n"Fe", joinpath(zeptest,"Standards","Fe std.msa")),
        reference(n"K", joinpath(zeptest,"Standards","K std.msa")),
        reference(n"Mg", joinpath(zeptest,"Standards","Mg std.msa")),
        reference(n"Mn", joinpath(zeptest,"Standards","Mn std.msa")),
        reference(n"Na", joinpath(zeptest,"Standards","Na std.msa")),
        reference(n"Nd", joinpath(zeptest,"Standards","Nd std.msa")),
        reference(n"Ni", joinpath(zeptest,"Standards","Ni std.msa")),
        reference(n"O", joinpath(zeptest,"Standards","O std.msa")),
        reference(n"P", joinpath(zeptest,"Standards","P std.msa")),
        reference(n"Pb", joinpath(zeptest,"Standards","Pb std.msa")),
        reference(n"S", joinpath(zeptest,"Standards","S std.msa")),
        reference(n"Sb", joinpath(zeptest,"Standards","Sb std.msa")),
        reference(n"Si", joinpath(zeptest,"Standards","Si std.msa")),
        reference(n"Sn", joinpath(zeptest,"Standards","Sn std.msa")),
        reference(n"Sr", joinpath(zeptest,"Standards","Sr std.msa")),
        reference(n"Ti", joinpath(zeptest,"Standards","Ti std.msa")),
        reference(n"V", joinpath(zeptest,"Standards","V std.msa")),
        reference(n"W", joinpath(zeptest,"Standards","W std.msa")),
        reference(n"Zn", joinpath(zeptest,"Standards","Zn std.msa")),
        reference(n"Zr", joinpath(zeptest,"Standards","Zr std.msa")) ],
        132.0
    )
    @time quantify(zep, refs, withUncertainty=true)
    res = @time quantify(zep, refs, withUncertainty=false)
    @test isapprox(res[6,:FE], 94.87, atol=0.01)
    @test isapprox(res[26,:CA], 23.48, atol=0.01)
    @test isapprox(res[169, :BA], 39.38, atol=0.01)
    @test res[160,:FIRSTELM]==n"Si"
    @test res[172,:SECONDELM]==n"Ti"
    @test res[179,:THIRDELM]==n"Zn"
    @test ismissing(res[177,:FOURTHELM])
    @test res[195,:FOURTHELM]==n"K"
    @test classes(zep)==classes(res)
    outfile = tempname()*".hdz"
    NeXLParticle.writeZep(res, outfile)
    res2 = Zeppelin(outfile)
    @test classes(zep) == classes(res2)
    @test isapprox(res2[6,:FE], 94.87, atol=0.01)
    @test isapprox(res2[26,:CA], 23.48, atol=0.01)
    @test isapprox(res2[169, :BA], 39.38, atol=0.01)
    @test res2[160,:FIRSTELM]==n"Si"
    @test res2[172,:SECONDELM]==n"Ti"
    @test res2[179,:THIRDELM]==n"Zn"
    @test ismissing(res2[177,:FOURTHELM])
    @test res2[195,:FOURTHELM]==n"K"
    @test all(res[:,"CLASS"] .== zep[:,"CLASS"])
    @test all(res[:,"CLASS"] .== res2[:,"CLASS"])
    res3 = @time classify(res, NeXLParticle.BaseRules)
    writeZep(res3, joinpath(tempdir(),"tmp.hdz"))
    @test repr(res3[2, "CLASS"]) == "Pyrite (FeS2)"
    @test (r->(r.FE > 20) && (r.S > 20) && (r.FE + r.S > 80))(res3[2,:])
    @test repr(res3[20, "CLASS"]) == "Dolomite"
    @test (r -> (r.CA > 40) && (r.MG > 10) && (r.CA + r.MG > 70))(res3[20,:])
    rm(joinpath(tempdir(),"tmp.hdz"))
end

max_stars_repo_name: NicholasWMRitchie/NeXLParticle

max_issues_repo_name: NicholasWMRitchie/NeXLParticle

max_forks_repo_name: NicholasWMRitchie/NeXLParticle


text:
mutable struct Leadfield{T <: AbstractFloat, S <: AbstractString}
    L::Array{T, 3}
    x::Vector{T}
    y::Vector{T}
    z::Vector{T}
    sensors::Vector{S}
end

function Base.show(io::IO, l::Leadfield)
    @printf "Leadfield\n"
    @printf "  Number of sources: %d\n" size(l.L, 1)
    @printf "  Number of dimensions + nulls: %d\n" size(l.L, 2)
    @printf "  Number of sensors: %d\n" size(l.L, 3)
end

function match_leadfield(l::Leadfield, s::SSR)

    @info("Matching leadfield to SSR")

    idx = [findfirst(l.sensors, name) for name = channelnames(s)]

    if length(unique(idx)) < length(idx)
        error("Not all SSR channels mapped to sensor #SSR=$(length(channelnames(s))), #L=$(length(l.sensors))")
    end

    l.L = l.L[:,:,idx]
    l.sensors = l.sensors[idx]

    @debug("matched $(length(idx)) sensors")

    return l
end


text:
using MultipleScattering
using Plots
pyplot()

# A very sharp sinc function (dicrete delta)
w_arr = collect(LinRange(0.0,1.0,1000))
x0 = 10.;
t_arr = 4.:0.01:22;
fs = reshape([exp(im*w*x0) for w in w_arr],length(w_arr),1)
fts = frequency_to_time(fs,w_arr, t_arr);

# A gaussian covering same frequency range which leads to a time pulse with height 1
a = 3./maximum(w_arr)^2;
gauss_impulse = get_gaussian_freq_impulse(maximum(w_arr))
ft2s = frequency_to_time(fs,w_arr, t_arr; impulse = gauss_impulse);

# true impulse in time
true_fts = [exp(-(t-x0)^2/(4a)) for t in t_arr]
plot(t_arr, [fts,true_fts,ft2s])


text:
"""
    @abstractmethod

Macro used in generic functions that must be overloaded by derived types.
"""
macro abstractmethod()
  quote
    error("This function belongs to an interface definition and cannot be used.")
  end
end

"""
    @notimplemented
    @notimplemented "Error message"

Macro used to raise an error, when something is not implemented.
"""
macro notimplemented(message="This function is not yet implemented")
  quote
    error($(esc(message)))
  end
end

"""
    @notimplementedif condition
    @notimplementedif condition "Error message"

Macro used to raise an error if the `condition` is true
"""
macro notimplementedif(condition,message="This function is not yet implemented")
  quote
    if $(esc(condition))
      @notimplemented $(esc(message))
    end
  end
end

"""
    @unreachable
    @unreachable "Error message"

Macro used to make sure that a line of code is never reached.
"""
macro unreachable(message="This line of code cannot be reached")
  quote
    error($(esc(message)))
  end
end


text:
# # [Bezier curves and their acceleration](@id BezierCurvesTutorial)
#
# This tutorial illustrates how Bézier curves are generalized to manifolds and how to
# minimize their acceleration, i.e. how to get a curve that is as straight or as geodesic
# while fulfilling constraints
#
# This example also illustrates how to apply the minimization on the corresponding [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html) manifold using a [`gradient_descent`](@ref) with [`ArmijoLinesearch`](@ref).
#
# ## Table of contents
# * [Setup](@ref SetupTB)
# * [de Casteljau algorithm on manifolds](@ref Casteljau)
# * [Composite Bézire curves](@ref CompositeBezier)
# * [Minimizing the acceleration of a Bézier curve](@ref MinAccBezier)
# * [Literature](@ref LiteratureBT)
#
# ## [Setup](@id SetupTB)
#
# We first initialize the necessary packages
export_folder = joinpath( #src
    @__DIR__, #src
    "..", #src
    "..", #src
    "docs", #src
    "src", #src
    "assets", #src
    "images", #src
    "tutorials", #src
) #src
using Manopt, Manifolds
# and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)
using Colors
black = RGBA{Float64}(colorant"#000000")
TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
TolVibrantTeal = RGBA{Float64}(colorant"#009988")
geo_pts = collect(range(0.0, 1.0; length=101)) #hide
bezier_pts = collect(range(0.0, 3.0; length=201)) #hide
camera_position = (-1.0, -0.7, 0.3) #hide
nothing #hide
# Then we load our data, see [`artificial_S2_composite_bezier_curve`](@ref), a composite
# Bezier curve consisting of 3 segments on the Sphere
B = artificial_S2_composite_bezier_curve();
#
# ## [De Casteljau algorithm on manifolds](@id Casteljau)
# This curve can be evaluated using de [Casteljau's algorithm](https://en.wikipedia.org/wiki/De_Casteljau%27s_algorithm)[^Casteljau1959][^Casteljau1963] named after [Paul de Casteljau](https://en.wikipedia.org/wiki/Paul_de_Casteljau)(*1930).
# To simplify the idea and understand this algorithm, we will first only look at the points of the first segment
M = Sphere(2)
b = B[2].pts
# On Euclidean spaces Bézier curves of these $n=4$ so called control points like this segment yield polynomials of degree $3$.
# The resulting curve $\gamma: [0,1] \to \mathbb R^m$ is called [Bezier curve](https://en.wikipedia.org/wiki/Bézier_curve) or Bézier spline and is named after [Piérre Bezier](https://en.wikipedia.org/wiki/Pierre_Bézier) (1910–1999).
# They can be evaluated by the de Casteljau algorithm by evaluating line segments between points.
# While it is not easy to evaluate polynomials on a manifold, evaluating line segments generalizes to the evaluation of [`shortest_geodesic`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any})s
#
# We will illustrate this using these points $b=(b_1,b_2,b_3,b_4)$ on the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/sphere.html) $\mathbb S^2$.
# Let's evaliuate this at the point $t=\frac{1}{4}\in [0,1]$. We first compute
t = 0.66
pts1 = shortest_geodesic.(Ref(M), b[1:3], b[2:4], Ref(t))
# We obtain 3 points on the geodesics connecting the control points. Repeating this again twice
pts2 = shortest_geodesic.(Ref(M), pts1[1:2], pts1[2:3], Ref(t))
p = shortest_geodesic(M, pts2[1], pts2[2], t)
# we obtain the point on the Bézier curve $c(t)$.
# This procedure is illustrated in the following image:
#
curves1 = [shortest_geodesic(M, b[i], b[i + 1], geo_pts) for i in 1:3] #src
curves2 = [shortest_geodesic(M, pts1[i], pts1[i + 1], geo_pts) for i in 1:2] #src
curves3 = [shortest_geodesic(M, pts2[1], pts2[2], geo_pts)] #src
asymptote_export_S2_signals( #src
    export_folder * "/Casteljau-illustr.asy"; #src
    points=[b, pts1, pts2, [p]], #src
    curves=[curves1..., curves2..., curves3..., de_casteljau(M, B[2], geo_pts)], #src
    colors=Dict( #src
        :points => [TolVibrantBlue, TolVibrantCyan, TolVibrantTeal, TolVibrantOrange], #src
        :curves => [ #src
            TolVibrantBlue, #src
            TolVibrantBlue, #src
            TolVibrantBlue, #src
            TolVibrantCyan, #src
            TolVibrantCyan, #src
            TolVibrantTeal, #src
            black, #src
        ], #src
    ), #src
    dot_size=3.5, #src
    line_widths=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.5], #src
    camera_position=camera_position, #src
) #src
render_asymptote(export_folder * "/Casteljau-illustr.asy"; render=2) #src
#md # ![Illustration of de Casteljau's algorithm on the Sphere.](../assets/images/tutorials/Casteljau-illustr.png)
#
# From the control points (blue) and their geodesics, ont evaluation per geodesic yields three interims points (cyan), their two successive geodeics another two points (teal) and at its geodesic at $t=0.66$ we obtain the point on the curve.
#
# In Manopt.jl, to evaluate a Bézier curve knowing its [`BezierSegment`](@ref), use [`de_casteljau`](@ref).
#
# There are a few nice observations to make, that hold also for these Bézier curves on manifolds:
# * The curve starts in the first controlpoint $b_0$ and ends in the last controlpoint $b_3$
# * The tangent vector to the curve at the start $\dot c(0)$ is equal to $\log_{b_0}b_1 = \dot\gamma_{b_0,b_0}(0)$, where $\gamma_{a,b}$ denotes the shortest geodesic.
# * The tangent vector to the curve at the end $\dot c(1)$ is equal to $-\log_{b_3}b_2 = -\dot\gamma_{b_3,b_2}(0) = \dot\gamma_{b_2,b_3}(1)$.
# * the curve is differentiable.
#
# For more details on these properties, see for example [^PopielNoakes2007].
# ## [Composite Bézier curves](@id CompositeBezier)
# With the properties of a single Bézier curve, also called Bézier segment, we can “stitch” curves together. Let $a_0,\ldots,a_n$ and $b_0,\ldots,b_m$ be two sets of controlpoints for the Bézier segments $c(t)$ and $d(t)$, respectively.
# We define the composite Bézier curve by $B(t) = \begin{cases} c(t) & \text{ if } 0\leq t < 1, \\ d(t-1) & \text{ if } 1\leq t \leq 2,\end{cases}$ where $t\in[0,2]$.
# This can of course be generalised straight forward to more than two cases.
# With the properties from the previous section we can now state that
#
# * the curve $B(t)$ is continuous if $c(1)=d(0)$ or in other words $a_n=b_0$
# * the curve $B(t)$ is differentiable if additionally $\dot c(1)=\dot d(0)$ or in other words $-\log_{a_n}a_{n-1} = \log_{b_0}b_1$. This is equivalent to $a_n=b_0 = \gamma_{a_{n-1}b_1}(\tfrac{1}{2})$.
#
# One nice interpretation of the last characterization is, that the tangents $\log_{a_n}a_{n-1}$ and $\log_{b_0}b_1$ point into opposite directions.
# For a continuous curve, the first point of every segment (except for the first segment) can be ommitted, for a differentiable curve the first two points (except for the first segment) can be ommitted.
# You can reduce storage by calling [`get_bezier_points`](@ref), though for econstruciton with [`get_bezier_segments`](@ref) you also need [`get_bezier_degrees`](@ref).
# The reduced storage is represented as an array of points, i.e. an element of the corresponding [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html).
#
# For the three segment example from the beginning this looks as follows[^1]
#
asymptote_export_S2_signals( #src
    export_folder * "/Bezier-composite-curve.asy"; #src
    curves=[de_casteljau(M, B, bezier_pts)], #src
    points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)], #src
    tangent_vectors=[[ #src
        Tuple(a) #src
        for #src
        a in #src
        zip(get_bezier_junctions(M, B, true), get_bezier_junction_tangent_vectors(M, B)) #src
    ]], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantBlue, TolVibrantTeal], #src
        :tvectors => [TolVibrantCyan], #src
    ), #src
    camera_position=camera_position, #src
    arrow_head_size=10.0, #src
    line_widths=[1.5, 1.5], #src
    dot_size=4.0, #src
) #src
render_asymptote(export_folder * "/Bezier-composite-curve.asy"; render=2) #src
#md # ![Illustration of a differentiable composite Bézier curve with 3 segments.](../assets/images/tutorials/Bezier-composite-curve.png)
#
# ## [Minimizing the acceleration of a composite Bézier curve](@id MinAccBezier)
# The motivation to minimize the acceleration of the composite Bézier curve is, that the curve should get “straighter” or more geodesic like.
# If we discretize the curve $B(t)$ with its control points denoted by $b_{i,j}$ for the $j$th note in the $i$th segment, the discretized model for equispaced $t_i$, $i=0,\ldots,N$ in the domain of $B$ reads[^BergmannGousenbourger2018]
#
# ````math
# A(b) \coloneqq\sum_{i=1}^{N-1}\frac{\mathrm{d}^2_2 \bigl[ B(t_{i-1}), B(t_{i}), B(t_{i+1}) \bigr]}{\Delta_t^3},
# ````
#
# where $\mathrm{d}_2$ denotes the second order finite difference using the mid point approach, see [`costTV2`](@ref)[^BacakBergmannSteidlWeinmann2016],
#
# ````math
# d_2(x,y,z) := \min_{c ∈ \mathcal C_{x,z}} d_{\mathcal M}(c,y),\qquad x,y,z∈\mathcal M.
# ````
#
# Another model is based on logarithmic maps, see [^BoumalAbsil2011], but that is not considered here.
# An advantage of the model considered here is, that it only consist of the evaluation of geodesics.
# This yields a gradient of $A(b)$ with respect to $b$ [`adjoint_Jacobi_field`](@ref)s. The following image shows the negative gradient (scaled)
#
gradFullB = Manopt._∇acceleration_bezier( #src
    M, #src
    get_bezier_points(M, B, :differentiable), #src
    [3, 3, 3], #src
    collect(range(0.0, 3.0; length=151)), #src
) #src
asymptote_export_S2_signals( #src
    export_folder * "/Bezier-composite-curve-gradient.asy"; #src
    curves=[de_casteljau(M, B, bezier_pts)], #src
    points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)], #src
    tangent_vectors=[[ #src
        Tuple(a) #src
        for #src
        a in #src
        zip( #src
            get_bezier_points(M, B, :continuous), #src
            -0.05 .* get_bezier_points(M, gradFullB, :continuous), #src
        ) #src
    ]], #src
    colors=Dict( #src
        :curves => [black], #src
        :points => [TolVibrantBlue, TolVibrantTeal], #src
        :tvectors => [TolVibrantOrange], #src
    ), #src
    camera_position=camera_position, #src
    arrow_head_size=10.0, #src
    line_widths=[1.5, 1.5], #src
    dot_size=4.0, #src
) #src
render_asymptote(export_folder * "/Bezier-composite-curve-gradient.asy"; render=2) #src
#md # ![Illustration of the gradient of the acceleration with respect to the control points.](../assets/images/tutorials/Bezier-composite-curve-gradient.png)
#
# In the following we consider two cases: Interpolation, which fixes the junction and end points of $B(t)$
# and approximation, where a weight and a dataterm are additionally introduced.
#
# ### Interpolation
# For interpolation, the junction points are fixed and their gradient entries are hence set to zero.
# After transferring to the already mentioned [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html), we can then perform a [`gradient_descent`](@ref) as follows
curve_samples = collect(range(0.0, 3.0; length=151)) #exactness of approximating d^2
pB = get_bezier_points(M, B, :differentiable)
N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
F(pB) = cost_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples)
∇F(pB) = ∇acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples)
x0 = pB
pB_opt_ip = gradient_descent(
    N,
    F,
    ∇F,
    x0;
    stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.0001),
    stopping_criterion=StopWhenChangeLess(5 * 10.0^(-7)),
    debug=[
        :Iteration,
        " | ",
        :Cost,
        " | ",
        DebugGradientNorm(),
        " | ",
        DebugStepsize(),
        " | ",
        :Change,
        "\n",
        :Stop,
        10,
    ],
)
#
# and the result looks like
#
B_opt_ip = get_bezier_segments(M, pB_opt_ip, get_bezier_degrees(M, B), :differentiable) #src
asymptote_export_S2_signals( #src
    export_folder * "/Bezier-IP-Min.asy"; #src
    curves=[de_casteljau(M, B_opt_ip, bezier_pts), de_casteljau(M, B, bezier_pts)], #src
    points=[get_bezier_junctions(M, B_opt_ip), get_bezier_inner_points(M, B_opt_ip)], #src
    tangent_vectors=[[ #src
        Tuple(a) #src
        for #src
        a in zip( #src
            get_bezier_junctions(M, B_opt_ip, true), #src
            get_bezier_junction_tangent_vectors(M, B_opt_ip), #src
        ) #src
    ]], #src
    colors=Dict( #src
        :curves => [TolVibrantBlue, black], #src
        :points => [TolVibrantBlue, TolVibrantTeal], #src
        :tvectors => [TolVibrantCyan], #src
    ), #src
    camera_position=camera_position, #src
    arrow_head_size=10.0, #src
    line_widths=[1.5, 0.75, 1.5], #src
    dot_size=4.0, #src
) #src
render_asymptote(export_folder * "/Bezier-IP-Min.asy"; render=2) #src
#md # ![Interpolation Min Acc](../assets/images/tutorials/Bezier-IP-Min.png)
#
# ### Approximation
# Similarly if we introduce the junction points as data fixed given $d_i$ and set (for simplicity) $p_i=b_{i,0}$ and $p_{n+1}=b_{n,4}$
# and set $λ=3$ in
#
# ````math
# \frac{\lambda}{2}\sum_{k=0}^3 d_{\mathcal M}(d_i,p_i)^2 + A(b),
# ````
#
# then $λ$ models how important closeness to the data $d_i$ is.
#
λ = 3.0
d = get_bezier_junctions(M, B)
F(pB) = cost_L2_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples, λ, d)
∇F(pB) = ∇L2_acceleration_bezier(M, pB, get_bezier_degrees(M, B), curve_samples, λ, d)
x0 = pB
pB_opt_appr = gradient_descent(
    N,
    F,
    ∇F,
    x0;
    stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001),
    stopping_criterion=StopWhenChangeLess(10.0^(-5)),
    debug=[
        :Iteration,
        " | ",
        :Cost,
        " | ",
        DebugGradientNorm(),
        " | ",
        DebugStepsize(),
        " | ",
        :Change,
        "\n",
        :Stop,
        50,
    ],
)
#
# and the result looks like
#
B_opt_appr = get_bezier_segments(M, pB_opt_appr, get_bezier_degrees(M, B), :differentiable) #src
asymptote_export_S2_signals( #src
    export_folder * "/Bezier-Appr-Min.asy"; #src
    curves=[de_casteljau(M, B_opt_appr, bezier_pts), de_casteljau(M, B, bezier_pts)], #src
    points=[ #src
        d, #src
        get_bezier_junctions(M, B_opt_appr), #src
        get_bezier_inner_points(M, B_opt_appr), #src
    ], #src
    tangent_vectors=[[ #src
        Tuple(a) #src
        for #src
        a in zip( #src
            get_bezier_junctions(M, B_opt_appr, true), #src
            get_bezier_junction_tangent_vectors(M, B_opt_appr), #src
        ) #src
    ]], #src
    colors=Dict( #src
        :curves => [TolVibrantBlue, black], #src
        :points => [TolVibrantOrange, TolVibrantBlue, TolVibrantTeal], #src
        :tvectors => [TolVibrantCyan], #src
    ), #src
    camera_position=camera_position, #src
    arrow_head_size=10.0, #src
    line_widths=[1.5, 0.75, 1.5], #src
    dot_size=4.0, #src
) #src
render_asymptote(export_folder * "/Bezier-Appr-Min.asy"; render=2) #src
#
#md # ![Approximation min Acc](../assets/images/tutorials/Bezier-Appr-Min.png)
#
# The role of $\lambda$ can be interpreted as follows: for large values of $\lambda$, the
# minimizer, i.e. the resulting curve, is closer to the original Bézier junction points.
# For small $\lambda$ the resting curve is closer to a geodesic and the control points are closer to the curve.
# For $\lambda=0$ _any_ (not necessarily shortest) geodesic is a solution and the problem is ill-posed.
# To illustrate the effect of $\lambda$, the following image contains 1000 runs for $\lambda=10$ in dark currant to $\lambda=0.01$ in bright yellow.
#
#md # ![Approximation min Acc](../assets/images/tutorials/Bezier_Approximation_video-Summary-result.png)
#
# The effect of the data term can also be seen in the following video, which starts a little slow and takes about 40 seconds.
#
#md # ![Video of the effect of lambda, the weight of the dataterm](../assets/videos/tutorials/Bezier_Approximation_video-movie.mp4)
#
#
# ## [Literature](@id LiteratureBT)
#
# [^BacakBergmannSteidlWeinmann2016]:
#     > Bačák, M., Bergmann, R., Steidl, G. and Weinmann, A.: _A second order nonsmooth
#     > variational model for restoring manifold-valued images_,
#     > SIAM Journal on Scientific Computations, Volume 38, Number 1, pp. A567–597,
#     > doi: [10.1137/15M101988X](https://doi.org/10.1137/15M101988X),
#     > arXiv: [1506.02409](https://arxiv.org/abs/1506.02409)
# [^BergmannGousenbourger2018]:
#     > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
#     > by minimizing the acceleration of a Bézier curve_.
#     > Frontiers in Applied Mathematics and Statistics, 2018.
#     > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
#     > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
# [^BoumalAbsil2011]:
#     > Boumal, N. and Absil, P.-A.: _A discrete regression method on manifolds and its application to data on SO(n)._
#     > In: IFAC Proceedings Volumes (IFAC-PapersOnline). Vol. 18. Milano (2011). p. 2284–89.
#     > doi: [10.3182/20110828-6-IT-1002.00542](https://dx.doi.org/10.3389/10.3182/20110828-6-IT-1002.00542),
#     > web: [www](https://web.math.princeton.edu/~nboumal/papers/Boumal_Absil_A_discrete_regression_method_on_manifolds_and_its_application_to_data_on_SOn.pdf)
# [^Casteljau1959]:
#     > de Casteljau, P.: _Outillage methodes calcul_, Enveloppe Soleau 40.040 (1959),
#     > Institute National de la Propriété Industrielle, Paris.
# [^Casteljau1963]:
#     > de Casteljau, P.: _Courbes et surfaces à pôles_, Microfiche P 4147-1,
#     > André Citroën Automobile SA, Paris, (1963).
# [^PopielNoakes2007]:
#     > Popiel, T. and Noakes, L.: _Bézier curves and $C^2$ interpolation in Riemannian
#     > manifolds_. Journal of Approximation Theory (2007), 148(2), pp. 111–127.-
#     > doi: [10.1016/j.jat.2007.03.002](https://doi.org/10.1016/j.jat.2007.03.002).
# [^1]:
#     > The images are rendered using [`asymptote_export_S2_signals`](@ref).
#     > For code examples, see [Get started: Optimize!](@ref Optimize).


text:
# Autogenerated wrapper script for Tar_jll for x86_64-unknown-freebsd
export tar

using acl_jll
using Attr_jll
JLLWrappers.@generate_wrapper_header("Tar")
JLLWrappers.@declare_executable_product(tar)
function __init__()
    JLLWrappers.@generate_init_header(acl_jll, Attr_jll)
    JLLWrappers.@init_executable_product(
        tar,
        "bin/tar",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()