text:
struct FINITO_lbfgs_adaptive_DNN_iterable{R<:Real,Tf,Tg, TH, Tf_full,Tp} <: CIAO_iterable
    F::Tf          # smooth term
    F_full::Tf_full
    g::Tg                  # nonsmooth term
    opt_params::Tp                # initial point
    N::Int                    # of data points in the finite sum problem
    L::Maybe{Union{Array{R},R}}  # Lipschitz moduli of nabla f_i
    γ::Maybe{Union{Array{R},R}}  # stepsizes
    η::R                    # ls division parameter for γ
    β::R                    # ls division parameter for τ
    sweeping::Int8          # to only use one stepsize γ
    batch::Int              # batch size
    α::R                    # in (0, 1), e.g.: 0.99
    H::TH
    adaptive::Bool          # adaptive stepsizes
    tol_b::R
    data::Tuple
    DNN_config!
end

mutable struct FINITO_lbfgs_adaptive_DNN_state{R<:Real,Tx, TH}
    γ::Maybe{Tx}                   # stepsize parameter, plays the rolde of hat gamma
    hat_γ::R                # average γ
    av::Tx                  # the running average
    ind::Array{Array{Int}}  # running index set
    d::Int                  # number of batches
    H::TH                   # Hessian approx

    # some extra placeholders
    sum_z::Tx
    sum_fx::R
    sum_fv::R
    sum_nrmx::R             # sum |x_i|^2/γ_i
    sum_nabla::Tx           # sum \nabla f(x_i)
    sum_innprod::R          # sum  < \nabla f(x_i), x_i >
    val_fg::R               # value of g
    z::Tx
    ∇f_temp::Tx             # placeholder for gradients
    temp::Tx
    zbar::Tx                # bar z
    zbar_prev::Maybe{Tx}    # bar z previous
    vbar::Tx                # bar v
    vbar_prev::Maybe{Tx}    # bar v previous
    dir::Tx                 # direction
    xbar::Tx                # linesearch candidate
    inds::Array{Int}        # needed for shuffled only!
    τ::Float64              # number of epochs
end

function FINITO_lbfgs_adaptive_DNN_state(γ::Tx, hat_γ::R, av::Tx, ind, d, H::TH, sum_z, sum_nrmx, sum_nabla, sum_innprod, sum_fx) where {R,Tx,TH}
    return FINITO_lbfgs_adaptive_DNN_state{R,Tx,TH}(
        γ,
        hat_γ,
        av,
        ind,
        d,
        H,
        sum_z,
        sum_fx,
        R(0),
        sum_nrmx,
        sum_nabla,
        sum_innprod,
        R(0),
        copy(av),
        copy(av),
        copy(av),
        copy(av),
        nothing,
        copy(av),
        nothing,
        copy(av),
        copy(av),
        collect(1:d),
        1.0
        )
end

batchmemaybe(x) = tuple(x)
batchmemaybe(x::Tuple) = x

function Base.iterate(iter::FINITO_lbfgs_adaptive_DNN_iterable{R}) where {R}
    N = iter.N
    r = iter.batch # batch size
    # create index sets
    ind = Vector{Vector{Int}}(undef, 0)
    d = Int(floor(N / r))
    for i = 1:d
        push!(ind, collect(r*(i-1)+1:i*r))
    end
    r * d < N && push!(ind, collect(r*d+1:N))

    #initializing the vectors
    # sum_nabla = zero(iter.x0)
    # sum_innprod = R(0)
    # sum_fx = R(0)
    # ∇f = 0
    # fi_x = 0
    # for i = 1:N  # nabla f(x0)
    #     # ∇f, fi_x = gradient(iter.F(i), iter.x0)
    #     ∇f = Flux.gradient(() -> iter.F(i,iter.x0), params(iter.x0))[iter.x0][:,1]
    #     fi_x = iter.F(i,iter.x0)
    #     sum_innprod += real(dot(∇f, iter.x0))
    #     sum_nabla .+= ∇f
    #     sum_fx += fi_x
    # end

    # sum_nabla = Flux.gradient(() -> iter.F_full(iter.opt_params), params(iter.x0))[iter.x0][:,1] * N
    # sum_fx = iter.F_full(iter.x0) * N
    # sum_innprod = real(dot(sum_nabla, iter.x0))

    gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(batchmemaybe(iter.data)...)
    end
    sum_nabla = iter.DNN_config!(gs=gs) * N
    sum_fx = iter.F(batchmemaybe(iter.data)...) * N
    x0 = iter.DNN_config!()
    sum_innprod = real(dot(sum_nabla, x0))

    # println("size testing ∇f: $(size(∇f))")
    # println("size testing fi_x: $(size(fi_x))")
    # println("size testing real.dot: $(size(real(dot(∇f, iter.x0))))")
    # println("size testing sum_nabla: $(size(sum_nabla))")
    # updating the stepsize

    ##-------- γ initialization type 1-------------
    γ = iter.γ
    xeps = x0 .+ one(R)
    #    av_eps = zero(x0)
    #    for i in 1:N
    #         # ∇f, ~ = gradient(iter.F[i], xeps)
    #         ∇f = Flux.gradient(() -> iter.F(i), params(iter.x0))
    #         av_eps .+= ∇f[iter.x0]
    #     end

    iter.DNN_config!(gs=xeps)
    gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(batchmemaybe(iter.data)...)
    end
    av_eps = iter.DNN_config!(gs=gs) * N

    nmg = norm(sum_nabla - av_eps)
    t = 1
    while nmg < eps(R)  # in case xeps has the same gradient
        println("initial upper bound for L is too small")
        xeps = x0 .+ rand(t * [-1, 1], size(x0))
        # av_eps = zero(x0)
        # for i in 1:N
        #     # ∇f, ~ = gradient(iter.F[i], xeps)
        #     ∇f = Flux.gradient(() -> iter.F(i), params(iter.x0))
        #     av_eps .+= ∇f[iter.x0]
        # end
        iter.DNN_config!(gs=xeps)
        gs = Flux.gradient(iter.opt_params) do # sampled gradient
            iter.F(batchmemaybe(iter.data)...)
        end
        av_eps = iter.DNN_config!(gs=gs) * N

        # grad_f_xeps, f_xeps = gradient(iter.F[i], xeps)
        nmg = norm(sum_nabla - av_eps)
        t *= 2
    end
    ###### decide how to initialize!
    L_int = nmg / (t * sqrt(length(x0)))
    L_int /= iter.N # to account for 1/N
    γ = iter.α / (L_int)
    # γ = iter.S(iter.x0)
    println("gamma specified by L_int: ", γ); flush(stdout)

    ##-------- γ initialization type 2-------------
    # γ = 5.0
    # println("gamma specified by hand: ", γ); flush(stdout)
    ##---------------------------------------------

    γ = isa(γ, R) ? N * fill(γ, (N,)) : γ
    hat_γ = 1 / sum(1 ./ γ)
    # update the average
    av = copy(sum_nabla) .* (-hat_γ / N)
    av .+= x0
    println("size testing x0: $(size(x0))"); flush(stdout)

    sum_nrmx = norm(x0)^2 / hat_γ
    sum_z = x0 / hat_γ
    state = FINITO_lbfgs_adaptive_DNN_state(γ, hat_γ, av, ind, cld(N, r), iter.H, sum_z, sum_nrmx, sum_nabla, sum_innprod, sum_fx)

    return state, state
end

function Base.iterate(
    iter::FINITO_lbfgs_adaptive_DNN_iterable{R},
    state::FINITO_lbfgs_adaptive_DNN_state{R},
) where {R}

    # println("outer iter started 1"); flush(stdout)

    if state.zbar_prev === nothing
        state.zbar_prev = zero(state.z)
        state.vbar_prev = zero(state.z)
    end

    γ_init = state.hat_γ

    # full z update
     while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 1, ($(state.hat_γ))"
            return nothing
        end
        prox!(state.zbar, iter.g, state.av, state.hat_γ)

        # sum_fz = 0
        # for i in 1:iter.N
        #     sum_fz += iter.F(i,state.zbar)
        # end
        # sum_fz = iter.F_full(state.zbar) * iter.N
        iter.DNN_config!(gs=state.zbar)
        sum_fz = iter.F(batchmemaybe(iter.data)...) * iter.N

        ls_temp_l = sum_fz
        ls_temp_l -= real(dot(state.sum_nabla, state.zbar))
        ls_temp_l -= state.sum_fx
        ls_temp_l += state.sum_innprod
        temp = (norm(state.zbar)^2)/(2*state.hat_γ) + (state.sum_nrmx)/2
        ls_temp_r =  (1 + 10^-6) * temp - real(dot(state.sum_z, state.zbar))
        ls_temp_r *= iter.N * iter.α

        tol = 10^(-6)  * (1 + abs(ls_temp_r))
        # break
        R(ls_temp_l) <= ls_temp_r + tol && break
        println("ls 1"); flush(stdout)
        # println("x: ", ls_temp_l)
        # println("y: ", ls_temp_r)
        # println("hat_γ: ", state.hat_γ)
        # println("hat_calc_γ: ", 1 / sum(1 ./ state.γ))
        # println("norm: ", state.sum_nrmx)
        # println("p1: ", (norm(state.zbar)^2)/(2*state.hat_γ) + (state.sum_nrmx)/2)
        # println("p2: ", real(dot(state.sum_z, state.zbar)))

        γ_prev = state.hat_γ
        state.hat_γ *= iter.η
        # update av
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla
        reset!(state.H)
    end

    # state.sum_fx = R(0)
    # state.sum_innprod = R(0)
    # state.sum_nabla = zero(state.av)
    # for i = 1:iter.N
    #     state.∇f_temp .= (Flux.gradient(() -> iter.F(i,state.zbar), params(state.zbar)))[state.zbar][:,1]
    #     fi_z = iter.F(i,state.zbar)
    #     # state.∇f_temp, fi_z = gradient(iter.F[i], state.zbar) # update the gradient
    #     state.sum_fx += fi_z
    #     state.sum_nabla .+= state.∇f_temp
    # end
    # state.sum_nabla .= (Flux.gradient(() -> iter.F_full(state.zbar), params(state.zbar)))[state.zbar][:,1] * iter.N
    # state.sum_fx = iter.F_full(state.zbar) * iter.N

    gs = Flux.gradient(iter.opt_params) do # sampled gradient
        iter.F(batchmemaybe(iter.data)...)
    end
    state.sum_nabla .= iter.DNN_config!(gs=gs) * iter.N
    state.sum_fx = iter.F(batchmemaybe(iter.data)...) * iter.N
    state.sum_innprod = real(dot(state.sum_nabla,state.zbar))
    state.av .= state.zbar .- (state.hat_γ / iter.N) .* state.sum_nabla
    state.sum_nrmx = norm(state.zbar)^2 * iter.N

    # println("outer iter started 2"); flush(stdout)

    # full v update
    while true
        if state.hat_γ < iter.tol_b / iter.N
            @warn "parameter `γ` became too small at ls 2, ($(state.hat_γ))"
            return nothing
        end
        state.val_fg = prox!(state.vbar, iter.g, state.av, state.hat_γ) # bar v update
        # state.sum_fv = 0
        # for i in 1:iter.N
        #     state.sum_fv += iter.F(i,state.vbar)
        # end
        # state.sum_fv = iter.F_full(state.vbar) * iter.N

        iter.DNN_config!(gs=state.vbar)
        state.sum_fv = iter.F(batchmemaybe(iter.data)...) * iter.N

        f_model_cnst = state.sum_nrmx / (2 * state.hat_γ) + state.sum_fx - state.sum_innprod

        f_model_z = real(dot(state.vbar, state.av))
        f_model_z -=  norm(state.vbar)^2 / 2
        f_model_z *= iter.N / state.hat_γ

        f_model_z += state.sum_fv
        tol = 10^(-6)  * (1 + abs(f_model_cnst))
        # tol = eps(R)
        # break
        R(f_model_z) <= f_model_cnst + tol && break
        println("ls 2")
        flush(stdout)

        γ_prev = state.hat_γ
        state.hat_γ *= iter.η
        # update av
        state.av .-=  ((state.hat_γ - γ_prev) / iter.N) .* state.sum_nabla
        reset!(state.H)
    end
    state.vbar .-= state.zbar # \bar v- \bar z

    # prepare for linesearch: compute varphi(bar z)
    envVal = state.sum_fx / iter.N
    envVal += state.val_fg
    envVal += real(dot(state.sum_nabla, state.vbar)) / iter.N
    envVal += norm(state.vbar)^2 / (2 *  state.hat_γ)

    # update lbfgs
    update!(state.H, state.zbar - state.zbar_prev, -state.vbar +  state.vbar_prev)
    mul!(state.dir, state.H, state.vbar)
    # store vectors for next update
    copyto!(state.zbar_prev, state.zbar)
    copyto!(state.vbar_prev, state.vbar)
    # println("outer iter started 3"); flush(stdout)

    state.τ = 1
    γ_prev = state.hat_γ
    # while true  ### z_line should be removed using some other placeholder
    for i=1:5
        state.hat_γ = γ_prev
        state.xbar .=  state.zbar .+ (1- state.τ) .* state.vbar + state.τ * state.dir

        # compute varphi(xbar)
        # state.sum_nabla = zero(state.temp)
        # state.sum_fx = 0
        # for i = 1:iter.N
        #     # state.∇f_temp, fi_z = gradient(iter.F[i], state.xbar) # update the gradient
        #     state.∇f_temp .= (Flux.gradient(() -> iter.F(i,state.xbar), params(state.xbar)))[state.xbar][:,1]
        #     fi_z = iter.F(i,state.xbar)
        #     state.sum_fx += fi_z / iter.N
        #     state.sum_nabla .+= state.∇f_temp
        # end
        # state.sum_nabla .= (Flux.gradient(() -> iter.F_full(state.xbar), params(state.xbar)))[state.xbar][:,1] * iter.N
        # state.sum_fx = iter.F_full(state.xbar)

        iter.DNN_config!(gs=state.xbar)
        gs = Flux.gradient(iter.opt_params) do # sampled gradient
            iter.F(batchmemaybe(iter.data)...)
        end
        state.sum_nabla .= iter.DNN_config!(gs=gs) * iter.N
        state.sum_fx = iter.F(batchmemaybe(iter.data)...)

        while true
            if state.hat_γ < iter.tol_b / iter.N
                @warn "parameter `γ` became too small at ls 3, ($(state.hat_γ))"
                return nothing
            end
            state.av .= state.xbar
            state.av .-= (state.hat_γ / iter.N) .* state.sum_nabla
            state.val_fg = prox!(state.z, iter.g, state.av, state.hat_γ)

            # temp = 0
            # for i = 1:iter.N
            #     # ~, fi_z = gradient(iter.F[i], state.z) # update the gradient
            #     fi_z = iter.F(i,state.z)
            #     temp += fi_z / iter.N
            # end
            # temp = iter.F_full(state.z)
            iter.DNN_config!(gs=state.z)
            temp = iter.F(batchmemaybe(iter.data)...)

            envVal_trial = 0
            envVal_trial += state.sum_fx
            state.z .-= state.xbar # r_γ(z^trial)
            envVal_trial += real(dot(state.sum_nabla, state.z)) / iter.N
            envVal_trial += norm(state.z)^2 / (2 *  state.hat_γ)

            tol = 10^(-6)  * (1 + abs(envVal_trial))
            # tol = eps(R)
            # break
            if temp <= envVal_trial + tol
                # println("ls condition met")
                break
            end
            state.hat_γ *= iter.η
            println("ls 3"); flush(stdout)
        end
        envVal_trial += state.val_fg

        tol = 10^(-6) * abs(envVal)
        envVal_trial <= envVal + tol && break
        state.τ *= iter.β
        println("ls on τ"); flush(stdout)
        # println("envVal_trial: ", envVal_trial)
        # println("envVal      : ", envVal)
    end
    # println("outer iter started 4"); flush(stdout)

    state.zbar .= state.xbar # of line search
    ###### at the momoent I am computing the next prox twice (it is already computed in the ls)
    iter.sweeping == 3 && (state.inds = randperm(state.d)) # shuffled
    state.temp .= state.av

    state.sum_fx = R(0)
    state.sum_z = zero(state.temp)
    state.sum_innprod =  R(0)
    state.sum_nrmx =  R(0)
    while_flag = false
    flag = true

    state.γ *= (state.hat_γ / γ_init)
    state.zbar .= state.xbar

    for j in state.inds
        for i in state.ind[j]
            # backtrack γ (warn if γ gets too small)

            # println("γ_i: ", iter.S(state.zbar, i)/iter.N)
            #--------------------------- restarting section ---------------------------
            # hat_γ_prev = state.hat_γ
            # γ_prev = state.γ[i]
            # # gradient!(state.∇f_temp, iter.F[i], state.zbar)
            # # state.γ[i] = iter.S(state.zbar, state.∇f_temp, i)
            # state.γ[i] = iter.S(state.zbar, i)
            # state.hat_γ = 1/(1/state.hat_γ - 1/γ_prev + 1/state.γ[i])
            # # state.hat_γ = 1 / sum(1 ./ state.γ)
            # state.av .+= (hat_γ_prev - state.hat_γ) / iter.N * state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.xbar
            #--------------------------- restarting section ---------------------------

            # print("current gamma: $(state.hat_γ) - current index: $(state.γ[i])")
            temp_x = iter.data[1][:,i]
            temp_y = iter.data[2][:,i]
            # if mod(i,100) == 0
            #     println("inner iter: $(i)")
            # end
            while true
                if state.hat_γ < iter.tol_b / iter.N
                    @warn "parameter `γ` became too small in ls 4 ($(state.hat_γ)) - for index $(state.γ[i])"
                    return nothing
                end
                prox!(state.zbar, iter.g, state.av, state.hat_γ)

                # fi_x = gradient!(state.∇f_temp, iter.F[i], state.xbar) # update the gradient
                # state.∇f_temp .= (Flux.gradient(() -> iter.F[i](state.xbar), params(state.xbar)))[state.xbar][:,1]
                # fi_x = iter.F[i](state.xbar)

                iter.DNN_config!(gs=state.xbar)
                gs = Flux.gradient(iter.opt_params) do # sampled gradient
                    iter.F(temp_x,temp_y)
                end
                state.∇f_temp .= iter.DNN_config!(gs=gs)
                fi_x = iter.F(temp_x,temp_y)

                # state.val_fg = iter.F[i](state.zbar)  ##### change to function value! #####
                iter.DNN_config!(gs=state.zbar)
                state.val_fg = iter.F(temp_x,temp_y)

                fi_model = fi_x + real(dot(state.∇f_temp, state.zbar .- state.xbar)) +
                    (0.5 * iter.α * iter.N/ state.γ[i]) * (norm(state.zbar .- state.xbar)^2)

                tol = 10^(-6)  * (1 + abs(fi_model))
                # break
                R(state.val_fg) <= fi_model + tol && break

                # println("ls 4 γ: $(state.hat_γ) for sample $(i)"); flush(stdout)
                # println("x: ", state.val_fg)
                # println("y: ", fi_model)

                hat_γ_prev = state.hat_γ
                γ_prev = state.γ[i]
                state.γ[i] *= iter.η # detach point!!
                # state.hat_γ = 1 / sum(1 ./ state.γ)
                state.hat_γ = 1/(1/state.hat_γ - 1/γ_prev + 1/state.γ[i]) # bug prone
                state.av .+= (hat_γ_prev - state.hat_γ) / iter.N * state.sum_nabla + (state.hat_γ/state.γ[i] - hat_γ_prev/γ_prev) * state.xbar

                reset!(state.H)
            end

            # iterates
            state.av .+= (state.hat_γ / iter.N) .* state.∇f_temp
            state.sum_nabla .-= state.∇f_temp # update sum nabla f_i for next iter
            # gradient!(state.∇f_temp, iter.F[i], state.zbar) # update the gradient
            # state.∇f_temp .= (Flux.gradient(() -> iter.F[i](state.zbar), params(state.zbar)))[state.zbar][:,1]
            iter.DNN_config!(gs=state.zbar)
            gs = Flux.gradient(iter.opt_params) do # sampled gradient
                iter.F(temp_x,temp_y)
            end
            state.∇f_temp .= iter.DNN_config!(gs=gs)

            state.sum_nabla .+= state.∇f_temp # update sum nabla f_i for next iter
            state.av .-= (state.hat_γ / iter.N) .* state.∇f_temp
            state.av .+= state.hat_γ * (state.zbar .- state.xbar) ./ state.γ[i]

            # updates for the next linesearch
            state.sum_z += state.zbar / state.γ[i]
            state.sum_fx += state.val_fg
            state.sum_nrmx   += norm(state.zbar)^2 / state.γ[i]
            state.sum_innprod += real(dot(state.∇f_temp, state.zbar))
            # println("outer iter started 5"); flush(stdout)

        end
    end
    # println("current stepsize is $(state.hat_γ)")
    return state, state
end
solution(state::FINITO_lbfgs_adaptive_DNN_state) = state.zbar
epoch_count(state::FINITO_lbfgs_adaptive_DNN_state) = state.τ   # number of epochs is 2+ 1/tau , where 1/tau is from ls

#### TODO:
    # careful with the minibatch and adaptive stepsizes!
    # stepsizes if gamma is given...
    # add LS parameter as iter field


text:
import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, PT, finish
import MechanicalSketch: settext, place_image, PALETTE, color_with_lumin
import MechanicalSketch: circle, ∙, arrow_nofill, Point, @layer, sethue
import MechanicalSketch: latexify, arrow
if !@isdefined N
    @import_expand(~m, s, °, N)
end
include("test_functions_39.jl")
empty_figure(filename = joinpath(@__DIR__, "test_39.png");
      backgroundcolor = color_with_lumin(PALETTE[5], 90),
      hue = "black", height = 10m);
pt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
settext("<b>Flow around a cylinder</b>", pt, markup = true)
pt += (0, EM)
pt += (1m, -1m)
r = 1.0m
@layer begin
    sethue(PALETTE[8])
    circle(pt; r, action=:fill)
end
vr =  √2 / 2 ∙(r, r)
arrow_nofill(pt + 0.5vr, pt + vr;  arrowheadlength = 0.5EM)
settext("$r", pt + vr)
arrow(pt, 2EM, 0°, 30°; arrowheadlength = 0.4EM)
pt += (0, 1.5EM)
settext("θ", pt + vr)
#
#expression = :( δ / δr ∙ (δψ / δθ)  + δ / δθ ∙ (-δψ / δr)  = 0)

finish()

ext: jl


text:
using ThreeBodyProblem
using Test

@test AP2a(1e6,1e5) == 5.5e5
@test AP2a(1e6,1e6) == 1e6

p = [1;0;0]
rv = [1;3;4]
@test_broken findrP(rv,p) == 5
rv = [rv,rv]
@test_broken findrP(rv,p) == 5
rv = [1;3;4;5;5;7]
@test_broken findrP(rv,p) == 5
rv = [rv,rv]
@test_broken findrP(rv,p) == 5

# @test_broken M2e

rv = [1,0,0,0,1,0]
μ = 1
a, e, i, Ω, ω, ν, Π, u, l, ℰ = cart2oe(rv, μ, ang_unit=:rad)
@test a == 1.0
@test e == 0.0
@test i == 0.0
@test isnan(Ω)
@test isnan(ω)
@test isnan(ν)
@test isnan(Π)
@test isnan(u)
@test l == 0.0
@test ℰ == -0.5

rv = [1,0,0,0,0,1]
μ = 1
a, e, i, Ω, ω, ν, Π, u, l, ℰ = cart2oe(rv, μ, ang_unit=:rad)
@test a == 1.0
@test e == 0.0
@test i == π/2
@test Ω == 0.0
@test isnan(ω)
@test isnan(ν)
@test isnan(Π)
@test u == 0.0
@test isnan(l)
@test ℰ == -0.5

rv = [1,0,0,0,2,0]
μ = 1
a, e, i, Ω, ω, ν, Π, u, l, ℰ = cart2oe(rv, μ, ang_unit=:rad)
@test a == -0.5
@test e == 3.0
@test i == 0.0
@test isnan(Ω)
@test isnan(ω)
@test ν == 0.0
@test Π == 0.0
@test isnan(u)
@test isnan(l)
@test ℰ == 1.0

@test_throws ErrorException("ang_unit should be :rad or :deg") cart2oe(rv, μ, ang_unit=:foo)

rv = [1,0,0,0,√2,√2]
μ = 1
a, e, i, Ω, ω, ν, Π, u, l, ℰ = cart2oe(rv, μ, ang_unit=:deg)
@test oe2cart(a, e, i, Ω, ω, ν, μ, ang_unit=:deg) == rv
a, e, i, Ω, ω, ν, Π, u, l, ℰ = cart2oe(rv, μ, ang_unit=:rad)
@test oe2cart(a, e, i, Ω, ω, ν, μ, ang_unit=:rad) == rv
@test_throws ErrorException("ang_unit should be :rad or :deg") oe2cart(a, e, i, Ω, ω, ν, μ, ang_unit=:foo)


text:
# Simulations parameters
include("./decision_making.jl")
include("./reaction_time_analysis_functions.jl")
include("./distributions.jl")

param = Dict("T_rest"=>10000,
"mu0"=>30,
"noise_amp"=>0.02,
"I0E1"=>0.3255,
"I0E2"=>0.3255,
"JN11"=>0.2609,
"JN22"=>0.2609,
"JN12"=>0.0497,
"JN21"=>0.0497,
"JAext" => 0.00052,
"Tnmda" => 100,
"Tampa" => 2,
"gamma" => 0.641,
"Tstim"=>5000,
"coh"=>5,
"false_time"=>200,
"Threshold"=>20
)

fixed_param=Dict("a"=>270,"b"=>108,"d"=>0.1540)
#Parameters of the simulation for EUler methods
simu_param=Dict("dt"=>0.5,"time_wind"=>4,"slide_wind"=>4) #4=2/dt

using FileIO
using DataFrames
using CSV

#a=2
Tr=500
println(Tr)
# de 0.02 à 0.025 déjàa sauvegardé (voir 0.02,0.3 50)
df=DataFrame()
df[:Ir]=linspace(0.03,0.1,30)
for j=2:10
#CSV.write("Irsmall.csv",df)
 for i=1:25

for Ir in [0.035,0.047]
println(Ir)
Cohlist=[2,10,20]

 #Need to do the quantile analysis after with r s etc ...

for c in Cohlist[1:end]
c_min=-c
 c_max = c
 number = 2

   distri = make_distributions(c_min,c_max,number,3000)

 param["T_rest"] = Tr
 If=0

  result,r2= sim_mf_decisionmaking_pattern(distri,Ir,If,
    param,fixed_param,simu_param)

 CSV.write("./datatempAll/V$j-Coh$c-Final200Exp$i-T$Tr-Ir$Ir.csv",result)
 CSV.write("./datatempAll/V$j-NondynCoh$c-Final200Exp$i-T$Tr-Ir$Ir.csv",r2)

 end
end

end
end


text:
# This file defines some struct types used in the FrameConnection package.

"""
    ParamStruct(initialdensity::Vector{Float64},
        k_on::Float64, k_off::Float64, k_bleach::Float64, p_miss::Float64,
        nsigmadev::Float64, maxframegap::Int, nnearestclusters::Int)

Structure of parameters needed for frame-connection.

# Fields
- `initialdensity`: Density of emitters at the start of the experiment.
                    (see estimatedensities()) (emitters/pixel^2)
- `nnearestclusters`: Number of nearest preclusters used for local density
                      estimates. (default = 2)(see estimatedensities())
- `k_on`: Rate at which dark emitters convert to the visible state.
          (see estimateparams())(1/frame)
- `k_off`: Rate at which visible emitters are converted to the reversible dark
           state. (see estimateparams())(1/frame)
- `k_bleach`: Rate at which visible emitters are irreversibly photobleached.
              (see estimateparams())(1/frame)
- `p_miss`: Probability of missing a localization of a visible emitter.
- `nsigmadev`: Multiplier of localization errors that defines a pre-clustering
               distance threshold. (default = 5)(see precluster())(pixels)
- `maxframegap`: Maximum frame gap between temporally adjacent localizations in
                 a precluster. (default = 5)(see precluster())(frames)
- `nmaxnn`: Maximum number of nearest-neighbors inspected for precluster
            membership.  Ideally, this would be set to inf, but that's not
            feasible for most data. (default = 2)(see precluster())
"""
mutable struct ParamStruct
    initialdensity::Vector{Float64}
    nnearestclusters::Int
    k_on::Float64
    k_off::Float64
    k_bleach::Float64
    p_miss::Float64
    nsigmadev::Float64
    maxframegap::Int
    nmaxnn::Int
end
ParamStruct() = ParamStruct([], 2, 0.0, 0.0, 0.0, 0.0, 5.0, 5, 2)


text:
Base.size(A::_AbstractRecordArray) = getfield(A, :size)

_buffer(A) = getfield(A, :buffer)
_step(A) = getfield(A, :step)
stepof(A) = valueof(_step(A))

Base.parent(A::_AbstractRecordArray) = _buffer(A)::AbstractArray

@inline Base.pointer(A::_AbstractRecordArray) = getfield(A, :pointer)
@inline Base.pointer(A::_AbstractRecordArray, i::Integer) = pointer(A) + (i - 1) * stepof(A)

@inline Base.getindex(A::_AbstractRecordArray{<:Any,0}) = A[1]
@inline Base.setindex!(A::_AbstractRecordArray{<:Any,0}, x) = A[1] = x

@inline function Base.getindex(A::_AbstractRecordArray, i::Int)
    @boundscheck checkbounds(A, i)
    buf = _buffer(A)
    GC.@preserve buf begin
        return unwrap_union_value(unsafe_load(pointer(A, i)))
    end
end

@inline function Base.setindex!(A::_AbstractRecordArray, x, i::Int)
    @boundscheck checkbounds(A, i)
    buf = _buffer(A)
    v = field_convert(eltype(A), x)
    GC.@preserve buf begin
        unsafe_store!(pointer(A, i), v)
    end
end

@propagate_inbounds Base.getindex(A::_AbstractRecordArray, I::Int...) =
    A[LinearIndices(A)[I...]]

@propagate_inbounds Base.setindex!(A::_AbstractRecordArray, v, I::Int...) =
    A[LinearIndices(A)[I...]] = v

function check_eltype(::Type{T}) where {T}
    if T isa Union
        check_eltype(T.a)
        check_eltype(T.b)
        return
    end
    ismutabletype(T) && noinline() do
        error("mutable types are not supported; got: $T")
    end
    isconcretetype(T) || noinline() do
        error("a concrete type is required; got: $T")
    end
    Base.datatype_pointerfree(T) || noinline() do
        error("containing GC-managed object: $T")
    end
    return
end


text:
# import Particles package
using Particles

# parameters
Δt = 0.1    # time step size
nt = 200    # number of time steps
np = 10000  # number of particles
nx = 16     # number of grid points

# random initial conditions
x₀ = randn(np)
v₀ = randn(np)

# shift x₀ to the interval [0,1]
xmax = ceil(maximum(abs.(x₀)))
x₀ .+= xmax
x₀ ./= 2*xmax

# concatenate x₀ and v₀
z₀ = hcat(x₀, v₀)'

# vector field
function lorentz_force!(ż, z, p::PoissonSolver)
    ż[1] = z[2]
    ż[2] = eval_field(p, z[1])
end

# create Poisson solver
p = PoissonSolver{eltype(z₀)}(nx)
solve!(p, x₀)

# create an Equation instance with callbacks
equ = Equation((ż, z) -> lorentz_force!(ż, z, p), z₀; f_preproc = z -> solve!(p, z[1,:]))

# create Simulation instance
sim = Simulation(equ, ExplicitEuler, Δt, nt)

# run simulation
z = sim()

# load Plots package
using Plots

# compute plot ranges
vmax = ceil(maximum(abs.(v₀)))
xlim = (0, 1)
vlim = (-vmax, +vmax)

# plot initial condition
scatter(mod.(z[1,:,0], 1), z[2,:,0],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
)

# save figure to file
savefig("vlasov_poisson_z₀.png")

# create animation
anim = @animate for n in 0:nt
    scatter(mod.(z[1,:,n], 1), z[2,:,n],
        marker = 3,
        xlim = xlim,
        ylim = vlim,
        title = "Vlasov-Poisson",
        legend = false,
        size = (800, 600)
    )
end

# save animation to file
gif(anim, "vlasov_poisson_anim.gif", fps=10)

max_stars_repo_name: michakraus/JuliaAtHomeAssignments

max_issues_repo_name: michakraus/JuliaAtHomeAssignments

max_forks_repo_name: michakraus/JuliaAtHomeAssignments


text:
function _find_ranges_for_closeback!(ranges, x, y, _fl::F1, _fr::F2, ::Val{T1}, ::Val{T2}; threads = true) where {T1, T2, F1, F2}
    @_threadsfor threads for i in 1:length(x)
        curr_start = ranges[i].start
        if length(ranges[i]) == 0
            ranges[i] = 0:0
        else
            ranges[i] = 1:searchsortedlast_join(_fr, y, _fl(x[i])::T1, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            if ranges[i].stop < curr_start
                ranges[i] = 0:curr_start
            end
        end
    end
end

function _find_ranges_for_closefor!(ranges, x, y, _fl::F1, _fr::F2, ::Val{T1}, ::Val{T2}; threads = true) where {T1, T2, F1, F2}
    @_threadsfor threads for i in 1:length(x)
        cur_stop = ranges[i].stop
        if length(ranges[i]) == 0
            ranges[i] = 0:0
        else
            ranges[i] = 1:searchsortedfirst_join(_fr, y, _fl(x[i])::T1, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            if ranges[i].stop > cur_stop
                ranges[i] = 0:cur_stop
            end
        end
    end
end

function _find_ranges_for_closenearest!(ranges, x, y, _fl::F1, _fr::F2, ::Val{T1}, ::Val{T2}; threads = true) where {T1, T2, F1, F2}
    @_threadsfor threads for i in 1:length(x)
        cur_stop = ranges[i].stop
        curr_start = ranges[i].start
        if length(ranges[i]) == 0
            ranges[i] = 0:0
        else
            #TODO we don't need to search whole range for the second search
            fval = searchsortedfirst_join(_fr, y, _fl(x[i])::T1, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))
            lval = searchsortedlast_join(_fr, y, _fl(x[i])::T1, ranges[i].start, ranges[i].stop, Base.Order.Forward, Val(T2))

            if fval > cur_stop
                ranges[i] = 0:lval
            elseif lval < curr_start
                ranges[i] = 0:fval
            else
                ranges[i] = min(fval, lval):max(fval, lval)
            end
        end
    end
end

_IST_(x,y)=true

function _close_notol_nn!(_res, ranges, x, rnn, lnn, bordervalue; threads = true)
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            r1 = rnn[ranges[i].stop]
            if ranges[i].start == 0
                if ismissing(bordervalue)
                    _res[i] = missing
                else
                    _res[i] = x[ranges[i].stop]
                end
            else
                r2 = rnn[ranges[i].start]
                if isless(abs(r2-lnn[i]), abs(r1-lnn[i]))
                    _res[i] = x[ranges[i].start]
                else
                    _res[i] = x[ranges[i].stop]
                end
            end
        end
    end
end
# instead of finding the nearest, set res = op(x[left], x[right]), in border it sets res = op(missing, x[border])
function _close_notol_nn!(_res, ranges, x, rnn, lnn, bordervalue, op::F; threads = true) where F
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            if ranges[i].start == 0
                _res[i] = op(missing, x[ranges[i].stop])
            else
                _res[i] = op(x[ranges[i].start],x[ranges[i].stop])
            end
        end
    end
end

function _close_notol_noborder!(_res, ranges, x, direction; threads = true)
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            if ranges[i].start == 0
                _res[i] = missing
            else

                if direction == :backward
                    _res[i] = x[ranges[i].start]
                else
                    _res[i] = x[ranges[i].stop]
                end
            end
        end
    end
end
function _close_notol_border!(_res, ranges, x, bordervalue; threads = true)
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            _res[i] = x[ranges[i].stop]
            if !bordervalue && ranges[i].start == 0
                _res[i] = missing
            end
        end
    end
end

function _close_tol_nn!(_res, ranges, x, rnn, lnn, tol, aem::F, bordervalue; threads = true) where F
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            r1 = rnn[ranges[i].stop]
            if ranges[i].start == 0
                if ismissing(bordervalue)
                    _res[i] = missing
                else
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = x[ranges[i].stop]
                    else
                        _res[i] = missing
                    end
                end
            else
                r2 = rnn[ranges[i].start]
                if isless(abs(r2-lnn[i]), abs(r1-lnn[i]))
                    if isless(abs(r2-lnn[i]), tol) && aem(r2, lnn[i])
                        _res[i] = x[ranges[i].start]
                    else
                        _res[i] = missing
                    end
                else
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = x[ranges[i].stop]
                    else
                        _res[i] = missing
                    end
                end
            end
        end
    end
end

# instead of finding the nearest, set res = op(x[left], x[right]), in border it sets res = op(missing, x[border]), op(missing, missing) must be defined
function _close_tol_nn!(_res, ranges, x, rnn, lnn, tol, aem::F1, bordervalue, op::F2; threads = true) where {F1, F2}
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            r1 = rnn[ranges[i].stop]
            if ranges[i].start == 0
                if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                    _res[i] = op(missing, x[ranges[i].stop])
                else
                    _res[i] = op(missing, missing)
                end
            else
                r2 = rnn[ranges[i].start]
                if isless(abs(r2-lnn[i]), tol) && aem(r2, lnn[i])
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = op(x[ranges[i].start], x[ranges[i].stop])
                    else
                        _res[i] = op(x[ranges[i].start], missing)
                    end
                else
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = op(missing, x[ranges[i].stop])
                    else
                        _res[i] = op(missing, missing)
                    end
                end
            end
        end
    end
end

function _close_tol_noborder!(_res, ranges, x, rnn, lnn, tol, aem::F, direction; threads = true) where F
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            if ranges[i].start == 0
                _res[i] = missing
            else
                if direction == :backward
                    r1 = rnn[ranges[i].start]
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = x[ranges[i].start]
                    else
                        _res[i] = missing
                    end
                else
                    r1 = rnn[ranges[i].stop]
                    if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                        _res[i] = x[ranges[i].stop]
                    else
                        _res[i] = missing
                    end
                end
            end
        end
    end
end
function _close_tol_border!(_res, ranges, x, rnn, lnn, tol, aem::F, bordervalue; threads = true) where F
    @_threadsfor threads for i in 1:length(ranges)
        if ranges[i] == 0:0
            _res[i] = missing
        else
            r1 = rnn[ranges[i].stop]
            if isless(abs(r1-lnn[i]), tol) && aem(r1, lnn[i])
                _res[i] = x[ranges[i].stop]
            else
                _res[i] = missing
            end
            if !bordervalue && ranges[i].start == 0
                _res[i] = missing
            end
        end
    end
end

function  _fill_right_cols_table_close!(_res, x, ranges, total, borderval, fill_val, direction; op = nothing, nn = false, rnn = nothing, lnn = nothing, tol = nothing, aem = _IST_, threads = true)
    if borderval == :nearest
        bordervalue = true
    elseif borderval == :none
        bordervalue = missing
    elseif borderval == :missing
        bordervalue = false
    end
    if tol === nothing
        if nn
            if op === nothing
                _close_notol_nn!(_res, ranges, x, rnn, lnn, bordervalue; threads = threads)
            else
                _close_notol_nn!(_res, ranges, x, rnn, lnn, bordervalue, op; threads = threads)
            end
        else
            if ismissing(bordervalue)
                _close_notol_noborder!(_res, ranges, x, direction; threads = threads)
            else
                _close_notol_border!(_res, ranges, x, bordervalue; threads = threads)
            end
        end
    else
        if nn
            if op === nothing
                _close_tol_nn!(_res, ranges, x, rnn, lnn, tol, aem, bordervalue; threads = threads)
            else
                _close_tol_nn!(_res, ranges, x, rnn, lnn, tol, aem, bordervalue, op; threads = threads)
            end
        else
            if ismissing(bordervalue)
                _close_tol_noborder!(_res, ranges, x, rnn, lnn, tol, aem, direction; threads = threads)
            else
                _close_tol_border!(_res, ranges, x, rnn, lnn, tol, aem, bordervalue; threads = threads)
            end
        end
    end

end

function _change_refpool_find_range_for_close!(ranges, dsl, dsr, r_perms, oncols_left, oncols_right, direction, lmf, rmf, j; nsfpaj = true, threads = true)
    var_l = _columns(dsl)[oncols_left[j]]
    var_r = _columns(dsr)[oncols_right[j]]
    l_idx = oncols_left[j]
    r_idx = oncols_right[j]
    if lmf
        _fl = getformat(dsl, l_idx)
    else
        _fl = identity
    end
    if rmf
        _fr = getformat(dsr, r_idx)
    else
        _fr = identity
    end

    T1 = Core.Compiler.return_type(_fl, Tuple{eltype(var_l)})

    if DataAPI.refpool(var_r) !== nothing && nsfpaj
        true && throw(ErrorException("we shouldn't end up here"))
    else
        T2 = Core.Compiler.return_type(_fr, Tuple{eltype(var_r)})
        if direction == :backward
            _find_ranges_for_closeback!(ranges, var_l, view(var_r, r_perms), _fl, _fr, Val(T1), Val(T2); threads = threads)
        elseif direction == :forward
            _find_ranges_for_closefor!(ranges, var_l, view(var_r, r_perms), _fl, _fr, Val(T1), Val(T2); threads = threads)
        elseif direction == :nearest
            _find_ranges_for_closenearest!(ranges, var_l, view(var_r, r_perms), _fl, _fr, Val(T1), Val(T2); threads = threads)
        end
    end
end

# border = :nearest | :missing | :none
function _join_closejoin(dsl, dsr::AbstractDataset, ::Val{T}; onleft, onright, makeunique = false, border = :nearest, mapformats = [true, true], stable = false, alg = HeapSort, accelerate = false, direction = :backward, inplace = false, tol = nothing,  allow_exact_match = true, op = nothing, method = :sort, threads = true) where T
    isempty(dsl) && return copy(dsl)
    if !allow_exact_match
        #aem is the function to check allow_exact_match
        aem = !isequal
        # missing is greater than anything
        if tol == nothing
            tol = missing
        end
    else
        aem = _IST_
    end
    oncols_left = onleft
    oncols_right = onright
    right_cols = setdiff(1:length(index(dsr)), oncols_right)
    if !makeunique && !isempty(intersect(_names(dsl), _names(dsr)[right_cols]))
        throw(ArgumentError("duplicate column names, pass `makeunique = true` to make them unique using a suffix automatically." ))
    end

    nsfpaj = true
    # if the column for close join is a PA we cannot use the fast path
    if DataAPI.refpool(_columns(dsr)[oncols_right[end]]) !== nothing
        nsfpaj = false
    end
    if length(oncols_left) > 1 && method == :hash
        ranges, a, idx, minval, reps, sz, right_cols_2= _find_ranges_for_join_using_hash(dsl, dsr, onleft[1:end-1], onright[1:end-1], mapformats, true, Val(T), threads = threads)
        filter!(!=(0), reps)
        pushfirst!(reps, 1)
        cumsum!(reps, reps)
        pop!(reps)
        grng = GIVENRANGE(idx, reps, Int[], length(reps))
        starts, idx, last_valid_range = _sort_for_join_after_hash(dsr, oncols_right[end], stable, alg, mapformats, nsfpaj, grng, threads = threads)
    else
        ranges = Vector{UnitRange{T}}(undef, nrow(dsl))
        idx, uniquemode = _find_permute_and_fill_range_for_join!(ranges, dsr, dsl, oncols_right, oncols_left, stable, alg, mapformats, accelerate && length(oncols_right) > 1; nsfpaj = nsfpaj, threads = threads)

        for j in 1:(length(oncols_left) - 1)
            _change_refpool_find_range_for_join!(ranges, dsl, dsr, idx, oncols_left, oncols_right, mapformats[1], mapformats[2], j; nsfpaj = nsfpaj, threads = threads)
        end
    end

    # if border = :none , we should use :nearest direction
    _change_refpool_find_range_for_close!(ranges, dsl, dsr, idx, oncols_left, oncols_right, border == :none ? :nearest : direction, mapformats[1], mapformats[2], length(oncols_left); nsfpaj = nsfpaj, threads = threads)
    total_length = nrow(dsl)

    if inplace
        newds = dsl
    else
        res = []
        for j in 1:length(index(dsl))
            push!(res,  _columns(dsl)[j])
        end
        if dsl isa SubDataset
            newds = Dataset(res, copy(index(dsl)))
        else
            newds = Dataset(res, Index(copy(index(dsl).lookup), copy(index(dsl).names), copy(index(dsl).format)))
        end
    end

    for j in 1:length(right_cols)
        _res = allocatecol(_columns(dsr)[right_cols[j]], total_length)
        if DataAPI.refpool(_res) !== nothing
            fill_val = DataAPI.invrefpool(_res)[missing]
            _fill_right_cols_table_close!(_res.refs, view(DataAPI.refarray(_columns(dsr)[right_cols[j]]), idx), ranges, total_length, border, fill_val, direction; nn = direction == :nearest, rnn = view(_columns(dsr)[oncols_right[end]], idx), lnn = _columns(dsl)[oncols_left[end]], tol = tol, aem = aem, op = op, threads = threads)
        else
            _fill_right_cols_table_close!(_res, view(_columns(dsr)[right_cols[j]], idx), ranges, total_length, border, missing, direction; nn = direction == :nearest, rnn = view(_columns(dsr)[oncols_right[end]], idx), lnn = _columns(dsl)[oncols_left[end]], tol = tol, aem = aem, op = op, threads = threads)
        end
        push!(_columns(newds), _res)
        new_var_name = make_unique([_names(dsl); _names(dsr)[right_cols[j]]], makeunique = makeunique)[end]
        push!(index(newds), new_var_name)
        setformat!(newds, index(newds)[new_var_name], getformat(dsr, _names(dsr)[right_cols[j]]))
    end

    if inplace
        _modified(_attributes(newds))
    end
    newds

end


text:
#######################
### Filtered output ###
#######################

type OutputFiltered{N} <: Output
    handles::NTuple{N,FASTQ.Writer}
    filenames::NTuple{N,String}
    outputdir::String
end

function OutputFiltered{N}(protocol::Protocol{N};
                           outputdir::String = ".",
                           kwargs...)

    mkpath(outputdir)

    filenames = map(protocol.readnames) do name
        joinpath(outputdir, name*".fastq")
    end

    handles = map(filenames) do filename
        FASTQ.Writer(open(filename, "w+"))
    end

    OutputFiltered{N}(handles,filenames,outputdir)
end

Base.close(o::OutputFiltered) = map(close,o.handles)

function Base.write{N}(oh::OutputFiltered{N},ir::InterpretedRecord{N})
    if ! ir.unmatched
        for i in 1:N
            write(oh.handles[i],ir.records[i])
        end
    end
end

function mergeoutput{N}(outputs::Vector{OutputFiltered{N}};
                        outputdir::String = ".",
                        kwargs...)
    pmap(1:N) do i
        name = joinpath(outputdir,"$i.fastq")
        filenames = [oh.filenames[i] for oh in outputs]
        catfiles(name,filenames)
        map(rm,filenames)
    end

    pmap(outputs) do oh
        rm(oh.outputdir,force=true,recursive=true)
    end
end

max_stars_repo_name: pwl/FASTQDemultiplexer

max_issues_repo_name: pwl/FASTQDemultiplexer

max_forks_repo_name: pwl/FASTQDemultiplexer


text:
### A Pluto.jl notebook ###
# v0.18.4

using Markdown
using InteractiveUtils

# ╔═╡ eca417a4-a001-11ec-1a38-dbe781d83ace
md"""
# HW 8

Eric Nguyen
"""

# ╔═╡ 67290e16-8033-4c8c-8a14-add2caa56ed8
md"""
### Problem 1

Consider a vector field ``\vec{A} = f \vec{C}`` where ``\vec{C}`` is some arbitrary constant vector and ``f(\vec{r})`` is a scalar field over some volume ``V``.
By Gauss' Theorem,

$\begin{align*}
∮_S \vec{A} ⋅ d\vec{S} &= ∫_V \vec{∇} ⋅ \vec{A} \,dV \\
∮_S (f \vec{C}) ⋅ d\vec{S} &= ∫_V \vec{∇} ⋅ (f \vec{C}) \,dV \\
\vec{C} ⋅ ∮_S f \,d\vec{S} &= \vec{C} ⋅ ∫_V \vec{∇} f \,dV \\
\implies ∮_S f \,d\vec{S} &= ∫_V \vec{∇} f \,dV \\
\end{align*}$

Thus, we have proven that

$∫_V \vec{∇} f \,dV = ∮_S f \,d\vec{S}$
"""

# ╔═╡ 6f759031-498a-48ca-91db-379484f59658
md"""
### Problem 2

The corresponding ``x``- and ``y``-components of the vector field and their respective partial derivatives are

$\begin{align*}
A_x &= \frac{-y}{x^2 + y^2} \implies \frac{∂A_x}{∂y} = \frac{2y^2}{(x^2 + y^2)^2} -
 \frac{1}{x^2 + y^2} = \frac{-x^2 + y^2}{(x^2 + y^2)^2} \\
A_y &= \frac{x}{x^2 + y^2} \implies \frac{∂A_y}{∂x} = -\frac{2x^2}{(x^2 + y^2)^2} + \frac{1}{x^2 + y^2} = \frac{-x^2 + y^2}{(x^2 + y^2)^2}
\end{align*}$

The curl of the vector field in both Cartesian and cylindrical polar coordinates are, respectively,

$\begin{align*}
\text{curl } \vec{A} = \vec{∇} × \vec{A} &= \hat{k} \left(\frac{∂A_y}{∂x} - \frac{∂A_x}{∂y}\right) = 0  &&\text{Cartesian} \\
\text{curl } \vec{A} = \vec{∇} × \vec{A} &= \hat{k} \frac{1}{r} \left(\frac{∂}{∂r} r A_{\phi} - \frac{∂A_r}{∂\phi}\right) = 0 &&\text{Cylindrical polar}
\end{align*}$

The line integral ``\displaystyle ∮ \vec{A} ⋅ d\vec{ℓ}`` around a closed circle of radius ``R`` in the ``xy`` plane, centered at the origin, is

$∮ \vec{A} ⋅ d\vec{ℓ} = ∫_0^{2\pi} \left(\frac{\hat{\phi}}{R}\right) ⋅ \left(\hat{\phi}\,R\right) dt = ∫_0^{2\pi} dt = 2\pi$

The results are inconsistent with the surface theorems because the line integral avoids the singular point at the origin whereas the surface theorems include the origin.
"""

# ╔═╡ 5ec3597d-0aff-4f67-9d6f-90dfd5fca82a
md"""
### Problem 3

The corresponding ``x``-, ``y``-, and ``z``-components of the vector field and their respective partial derivatives are

$\begin{align*}
A_x &= \frac{x}{(x^2 + y^2 + z^2)^{3/2}} \implies \frac{∂A_x}{∂x} = \frac{-2x^2 + y^2 + z^2}{(x^2 + y^2 + z^2)^{5/2}} \\
A_y &= \frac{x}{(x^2 + y^2 + z^2)^{3/2}} \implies \frac{∂A_y}{∂y} = \frac{x^2 - 2y^2 + z^2}{(x^2 + y^2 + z^2)^{5/2}} \\
A_z &= \frac{x}{(x^2 + y^2 + z^2)^{3/2}} \implies \frac{∂A_z}{∂z} = \frac{x^2 + y^2 - 2z^2}{(x^2 + y^2 + z^2)^{5/2}}
\end{align*}$

The divergence of the vector in both Cartesian and spherical polar coordinates are, respectively,

$\begin{align*}
\text{div } \vec{A} = \vec{∇} ⋅ \vec{A} &= \frac{∂A_x}{∂x} + \frac{∂A_y}{∂y} + \frac{∂A_z}{∂z} = 0 &&\text{Cartesian} \\
\text{div } \vec{A} = \vec{∇} ⋅ \vec{A} &= \frac{1}{r} \frac{∂}{∂r} (r A_r) = 0 &&\text{Spherical polar} \\
\end{align*}$

The surface integral ``\displaystyle ∮ \vec{A} ⋅ d\vec{S}`` around a sphere of radius ``R``, centered at the origin, is

$∮ \vec{A} ⋅ d\vec{S} = ∫_0^{2\pi} ∫_0^{\pi} ∫_0^R \left(\frac{1}{r^2}\right) r^2 \sin{\theta} \,dr \,d\theta \,d\phi = 2\pi R$

The results are inconsistent with the surface theorems because the surface integral avoids the singular point at the origin whereas the surface theorems include the origin.
"""

# ╔═╡ 301277aa-7465-46e0-9198-10669d567204
md"""
### Problem 4

First, we write

$u(x,t) = X(x) T(t)$

So the partial differential equation becomes

$4T \frac{d^2 X}{dx^2} = X \frac{dT}{dt} \implies \frac{1}{X} \frac{d^2 X}{dx^2} = \frac{1}{4T} \frac{dT}{dt}$

Set each side to a constant ``-k^2`` to get

$X''(x) = -k^2 X \quad\text{and}\quad T'(t) = -4k^2 T$

The general solutions are

$X(x) = A_x \cos{kx} + B_x \sin{kx} \qquad\text{and}\qquad T(t) = A_y e^{-4k^2t}$

So the solution to the PDE has the form

$u(x,t) = X(x) T(t) = (C \cos{kx} + D \sin{kx}) (e^{-4k^2 t})$

where ``C = A_x A_y`` and ``D = B_x A_y``.
The boundary condition ``u(0,t) = 0`` gives

$(C \cos{0} + D \sin{0}) (e^{-4k^2t}) = 0 \implies C = 0$

The other boundary condition ``u(2,t) = 0`` gives

$(D \sin{2k})(e^{-8k^2 t}) = 0 \implies \sin{2k} = 0 \implies 2k ≡ n\pi \quad\text{where}\quad  n ∈ ℤ \implies k ≡ \frac{n\pi}{2}$

The solution to the PDE has the form

$u(x,t) = D \sin{\left(\frac{n\pi x}{2}\right)} e^{\frac{-4n^2 \pi^2 t}{4}}$

The initial condition ``u(x,0) = \sin{(\pi x / 2)}`` gives

$\sin{\left(\frac{\pi x}{2}\right)} = D \sin{\left(\frac{n\pi x}{2}\right)} \implies D = 1, \,n = 1$

The final solution is

$u(x,t) = \sin{\left(\frac{\pi x}{2}\right)} e^{-\pi^2 t}$
"""

# ╔═╡ 7dd42ca7-1752-49cc-a292-d923c34e697a
md"""
### Problem 5

The Helmholtz Equation is

$\vec{∇}^2 f(r,\phi) + k^2 f(r,\phi) = 0$

We want to show that solutions for ``R(r)`` are the Bessel Functions ``J_m(kr)`` of integer order ``m`` by arriving at the Bessel Equation,

$r^2 R''(r) + r R'(r) + (k^2 r^2 - m^2) R(r) = 0$

Using separation of variables, we write

$f(r,\phi) = R(r) \Phi(\phi)$

So the PDE becomes

$\begin{align*}
\vec{∇}^2 (R(r) \Phi(\phi)) + k^2 R(r) \Phi(\phi) &= 0 \\
\frac{1}{r} R'(r) \Phi(\phi) + R''(r) \Phi(\phi) + \frac{1}{r^2} R(r) \Phi''(\phi) + k^2 R(r) \Phi(\phi) &= 0 \\
r R'(r) \Phi(\phi) + r^2 R''(r) \Phi(\phi) + R(r) \Phi''(\phi) + k^2 r^2 R(r) \Phi(\phi) &= 0 \\
\frac{r R'(r) + r^2 R''(r) + k^2 r^2 R(r)}{R(r)} &= -\frac{\Phi''(\phi)}{\Phi(\phi)}
\end{align*}$

Suppose both sides are equal to a constant ``m^2``.
However in this case, we are only interested in the equation for ``R(r)`` so we can ignore the equation for ``\Phi(\phi)``.
Notice the equation for ``R(r)`` is the Bessel Equation

$r^2 R''(r) + r R'(r) + (k^2 r^2 - m^2) R(r) = 0$

where the solution is the Bessel Function ``J_m(kr)`` of integer order ``m``.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.8.0-beta1"
manifest_format = "2.0"
project_hash = "da39a3ee5e6b4b0d3255bfef95601890afd80709"

[deps]
"""

# ╔═╡ Cell order:
# ╟─eca417a4-a001-11ec-1a38-dbe781d83ace
# ╟─67290e16-8033-4c8c-8a14-add2caa56ed8
# ╟─6f759031-498a-48ca-91db-379484f59658
# ╟─5ec3597d-0aff-4f67-9d6f-90dfd5fca82a
# ╟─301277aa-7465-46e0-9198-10669d567204
# ╟─7dd42ca7-1752-49cc-a292-d923c34e697a
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Define the default algorithms for various pieces of X-ray physics data.
# If you want to change the default data that NeXL uses, make the change in here.

import Unitful: @u_str, ustrip
import PhysicalConstants.CODATA2018: PlanckConstant, SpeedOfLightInVacuum

const plancksConstant = ustrip(PlanckConstant |> u"eV*s")
const hc = ustrip((PlanckConstant * SpeedOfLightInVacuum) |> u"eV*cm") # (plancks⋅speed-of-light)
const speedOfLight = ustrip(SpeedOfLightInVacuum |> u"cm/s")

"""
    edgeenergy(z::Int, ss::Int, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64
    edgeenergy(cxr::CharXRay, ::Type{<:NeXLAlgorithm}=FFASTDB)

Return the minimum energy (in eV) necessary to ionize the specified sub-shell in the specified atom
or the ionized shell for the specified characteristic X-ray.
"""
edgeenergy(z::Int, ss::Int)::Float64 =
    edgeenergy(z, ss, FFASTDB)
edgeenergy(cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB) =
    edgeenergy(cxr.z, cxr.transition.innershell.index, alg)

hasedge(z::Int, ss::Int) = hasedge(z, ss, FFASTDB)

"""
     energy(ass::AtomicSubShell, ty::Type{<:NeXLAlgorithm}=FFASTDB)
     energy(elm::Element, ss::SubShell, ty::Type{<:NeXLAlgorithm}=FFASTDB)

 The edge energy in eV for the specified AtomicSubShell

Example:

    julia> energy(n"Fe L3")
    708.0999999999999
    julia> energy(n"Fe", n"L3")
    708.0999999999999

    energy(elm::Element, tr::Transition)::Float64
    energy(cxr::CharXRay)

The characteristic X-ray energy for the specified element / transition or characteristic X-ray.
"""
energy(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    edgeenergy(ass.z, ass.subshell.index, alg)
energy(elm::Element, ss::SubShell, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    edgeenergy(z(elm), ss.index, alg)
energy(elm::Element, tr::Transition, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64 =
    energy(z(elm), tr.innershell.index, tr.outershell.index, alg)
energy(cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB) = energy(
        cxr.z,
        cxr.transition.innershell.index,
        cxr.transition.outershell.index,
        alg
    )

"""
    eachelement(alg::Type{<:NeXLAlgorithm} = FFASTDB)

Return the range of atomic numbers for which there is a complete set of energy, weight, MAC, ... data
"""
eachelement() = eachelement(FFASTDB)

"""
    strength(cxr::CharXRay)::Float64

The fraction of ionizations of `inner(cxr)` that relax via a characteristic X-ray resulting
from an electronic transition from `outer(cxr)` to `inner(cxr)`.

See also `weight(cxr)`.
"""
strength(cxr::CharXRay)::Float64 = strength(element(cxr), cxr.transition, CullenEADL)

"""
    subshellindices(z::Int, alg::Type{<:NeXLAlgorithm} = FFASTDB)

Return the shells occupied in a neutral, ground state atom of the specified atomic number.
"""
subshellindices(z::Int) = subshellindices(z, FFASTDB)

"""
    jumpratio(z::Int, ss::Int, ::Type{FFASTDB}) =

Compute the jump ratio.
"""
jumpratio(z::Int, ss::Int) = jumpratio(z, ss, FFAST)

"""
    mac(elm::Element, energy::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64
    mac(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::Float64

The mass absorption coefficient for an X-ray of the specified energy (eV) or
characteristic X-ray line in the specified element.  In cm²/g.
"""
mac(elm::Element, energy::Float64)::Float64 =
    mac(elm, energy, FFASTDB)

const userMacs = Dict{Tuple{Element, CharXRay}, Float64}()
mac(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm}=FFASTDB)::Float64 =
    get(userMacs, (elm, cxr), mac(elm, energy(cxr), alg))
set_user_mac!(elm::Element, cxr::CharXRay, mac::Float64) =
    userMacs[(elm, cxr)]=mac
delete_user_mac!(elm::Element, cxr::CharXRay) =
    delete!(userMacs, (elm, cxr))
clear_user_macs!() =
    empty!(userMacs)

"""
    macU(elm::Element, energy::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)
    macU(elm::Element, cxr::Float64, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue
    macU(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue

The mass absorption coefficient (with uncertainty estimate) for an X-ray of the specified energy (eV)
or characteristix X-ray line in the specified element.
"""
macU(elm::Element, energy::Float64)::UncertainValue =
    macU(elm, energy, FFASTDB)
macU(elm::Element, cxr::CharXRay, alg::Type{<:NeXLAlgorithm} = FFASTDB)::UncertainValue =
    macU(elm, energy(cxr), alg)

"""
    fluorescenceyield(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64

The fraction of `inner` sub-shell ionizations that relax via a characteristic X-ray resulting from an
electronic transition from `outer` to `inner`.
"""
fluorescenceyield(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64 =
    totalWeight(z, inner, inner, outer, alg)

"""
    characteristicyield(z::Int, ionized::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64

The fraction of `ionized` sub-shell ionizations that relax via a characteristic X-ray resulting from an
electronic transition from `outer` to `inner`.  This includes both direct transitions (where `outer`==`ionized`)
and cascade (where `outer` != `ionized` due to Coster-Kronig and previous decays.)
"""
characteristicyield(z::Int, ionized::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64 =
    totalWeight(z, ionized, inner, outer, alg)

"""
    characteristicXRayAvailable(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Float64

Is the weight associated with this transition greater than zero?
"""
charactericXRayAvailable(z::Int, inner::Int, outer::Int, alg::Type{<:NeXLAlgorithm} = CullenEADL)::Bool =
    isAvailable(z, inner, outer, alg)

"""
    ionizationcrosssection(ass::AtomicSubShell, energy::AbstractFloat, ty::Type{<:NeXLAlgorithm}=Bote2009)
    ionizationcrosssection(z::Int, ss::Int, energy::AbstractFloat)

Computes the absolute ionization crosssection (in cm²) for the specified AtomicSubShell and
electon energy (in eV) using the default algorithm.

Example:

    julia> (/)(map(e->NeXLCore.ionizationcrosssection(n"Fe K",e),[10.0e3,20.0e3])...)
    0.5672910174711278
"""
ionizationcrosssection(
    ass::AtomicSubShell,
    energy::AbstractFloat,
    ty::Type{<:NeXLAlgorithm} = Bote2009,
) = ionizationcrosssection(ass.z, ass.subshell.index, energy, ty)
ionizationcrosssection(z::Int, ss::Int, energy::AbstractFloat) =
    ionizationcrosssection(z, ss, energy, Bote2009)

    """
    strength(elm::Element, tr::Transition)::Float64

Return the nominal line strenth for the specified transition in the specified element.
The strength differs from the weight by the fluorescence yield.  Assumes an overvoltage of 4.0
"""
strength(elm::Element, tr::Transition, ty::Type{<:NeXLAlgorithm} = CullenEADL)::Float64 =
    ionizationfraction(z(elm), tr.innershell.index, 4.0) *
    fluorescenceyield(z(elm), tr.innershell.index, tr.outershell.index, ty)

"""
    fluorescenceyield(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL)::Float64

The fraction of relaxations from the specified shell that decay via radiative transition
rather than electronic (Auger) transition.  Does not include Coster-Kronig
"""
fluorescenceyield(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL)::Float64 = sum(map(
    s -> fluorescenceyield(ass.z, ass.subshell.index, s, alg),
    ass.subshell.index+1:length(allsubshells),
))

"""
    fluorescenceyieldcc(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL)::Float64

The fraction of relaxations from the specified shell that decay via radiative transition
rather than electronic (Auger) transition.  Includes Coster-Kronig
"""
function fluorescenceyieldcc(ass::AtomicSubShell, alg::Type{<:NeXLAlgorithm}=CullenEADL)::Float64
    f(ss) = sum(map(
        s -> fluorescenceyield(ass.z, ass.subshell.index, s, alg),
        ss.index+1:length(allsubshells),
    ))
    return sum(map(ss -> f(ss), ass.subshell.index+1:lastsubshell(shell(ass)).index))
end


text:
@reexport module _EllipticCurves

# http://jeremykun.com/2014/02/08/introducing-elliptic-curves/
# http://jeremykun.com/2014/02/24/elliptic-curves-as-python-objects/
# http://jeremykun.com/2014/03/19/connecting-elliptic-curves-with-finite-fields-a-reprise/

import Base: +, -, *, show, rand, in
import Nemo: divexact, contains, order, gen

export curve, point, samecurve, ideal, ring, field, isideal

using AutoHashEquals: @auto_hash_equals
using Brainstorm._Math: factor, factors, isprime
using Nemo: FieldElem, FinFieldElem, RingElem, root
using Reexport: @reexport

divexact(x::Number, y::Number) = x / y

include("curve.jl")
include("point.jl")
include("group.jl")

end


text:
@testset "Iterative Refinement for predictor (FourBar)" begin
    @var x a y b x_hat a_hat y_hat b_hat
    gamma = [Variable(Symbol(:gamma, i)) for i = 1:8]
    delta = [Variable(Symbol(:delta, i)) for i = 1:8]
    gamma_hat = [Variable(Symbol(:gamma_hat, i)) for i = 1:8]
    delta_hat = [Variable(Symbol(:delta_hat, i)) for i = 1:8]
    #system of polynomials
    D1 = map(1:8) do i
        (a_hat * x - delta_hat[i] * x) * gamma[i] +
        (a * x_hat - delta[i] * x_hat) * gamma_hat[i] +
        (a_hat - x_hat) * delta[i] +
        (a - x) * delta_hat[i] - delta[i] * delta_hat[i]
    end
    D2 = map(1:8) do i
        (b_hat * y - delta_hat[i] * y) * gamma[i] +
        (b * y_hat - delta[i] * y_hat) * gamma_hat[i] +
        (b_hat - y_hat) * delta[i] +
        (b - y) * delta_hat[i] - delta[i] * delta_hat[i]
    end
    D3 = [gamma[i] * gamma_hat[i] + gamma[i] + gamma_hat[i] for i = 1:8]
    F = System(
        [D1; D2; D3],
        [x; a; y; b; x_hat; a_hat; y_hat; b_hat; gamma; gamma_hat],
        [delta; delta_hat],
    )

    s = Complex{Float64}[
        -2.6370887+0.35329306im,
        -3.7223306+2.5267087im,
        -0.70272787-0.74162847im,
        -0.52112188-1.5249im,
        -0.29444549+0.2698171im,
        0.23877446+1.1593521im,
        -0.2226372+0.24523739im,
        -0.23666974+0.87507397im,
        2.5277031-6.3353057im,
        -0.88010494+1.3133313im,
        -0.73606124+0.15544107im,
        0.80744564-3.0613903im,
        -0.82122921+0.80545901im,
        0.68534751-2.5000172im,
        0.65872999-2.3719828im,
        -0.87767826-0.35435132im,
        -0.93290889+0.12048708im,
        -0.93106365-0.75512927im,
        1.8130785-1.6567022im,
        -0.85699423+0.24221833im,
        -0.73738109-1.1832401im,
        -0.81460307+0.2750148im,
        -0.80200623+0.28313097im,
        -0.12955283+2.5215805im,
    ]
    p = Complex{Float64}[
        7.297148065259444-3.7969299471707454im,
        0.06125202773860465+0.9251492821389162im,
        1.4586512320807385+0.641527879044937im,
        3.3162292880980537-2.979942348702623im,
        -0.6613217458488018+0.8321586519581641im,
        2.6212360706406703-2.6766427578876666im,
        2.4522388848243555-2.5938408784246976im,
        -1.7265529790073868-0.29337325907998707im,
        0.24011427062306998+1.375412535460542im,
        -0.19336068638852033+0.5669754468190339im,
        0.147943349839693-0.2784406507316879im,
        -0.01641513182688631+1.5963735195205522im,
        -0.33880136827427765+0.2772103200265199im,
        -0.19788506286822416+1.6632361622942196im,
        -0.26097545384785903+1.6764055030939693im,
        0.44012684131549085+1.0212717350758722im,
    ]
    q = Complex{Float64}[
        -0.35263257737537096+1.2598825875492277im,
        -0.5353479512983561-0.5963278468670689im,
        1.4096839876427854-1.2227632666172732im,
        -0.2550909822371234-0.641270846350799im,
        0.24180730168363096-0.46611650665985527im,
        0.29953332004980354-0.8469219548644632im,
        0.35386303195558433-0.3683866291866925im,
        1.2021910424897007+1.1148794002014533im,
        -0.35263257737537096-1.2598825875492277im,
        -0.5353479512983561+0.5963278468670689im,
        1.4096839876427854+1.2227632666172732im,
        -0.2550909822371234+0.641270846350799im,
        0.24180730168363096+0.46611650665985527im,
        0.29953332004980354+0.8469219548644632im,
        0.35386303195558433+0.3683866291866925im,
        1.2021910424897007-1.1148794002014533im,
    ]

    @testset "Fourbar - AD: $AD" for AD = 0:3
        tracker = Tracker(
            ParameterHomotopy(F, p, q),
            options = TrackerOptions(automatic_differentiation = AD),
        )
        res = track(tracker, s, 1, 0)
        @test is_success(res)
    end
end


text:
#! /usr/bin/julia

# Rosetta Code, Largest int from concatenated ints

function maxconcat{T<:Integer}(a::Array{T,1})
    b = map(string, a)
    b = sort(b, lt=(x,y)->x*y < y*x, rev=true)
    b = join(b, "")
    try
        b = parseint(b)
    catch
        b = parseint(BigInt, b)
    end
end

tests = {[1, 34, 3, 98, 9, 76, 45, 4],
         [54, 546, 548, 60],
         [1, 34, 3, 98, 9, 76, 45, 4, 54, 546, 548, 60]}

for t in tests
    println("Maxconcating ", t)
    println("    ", maxconcat(t))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module TestPairs

using Test
using LazyGroupBy

@testset begin
    @test pairs.(grouped(isodd, [0, 7, 3, 1, 5, 9, 4, 3, 0, 5])) == Dict(
        false => Dict(7 => 4, 9 => 0, 1 => 0),
        true => Dict(4 => 1, 10 => 5, 2 => 7, 3 => 3, 5 => 5, 8 => 3, 6 => 9),
    )
end

end  # module


text:
using ParticleFilters
using POMDPs
using POMDPModels
using POMDPLinter: @implemented
using Test
using POMDPPolicies
using POMDPSimulators
using POMDPModelTools
using Random
using Distributions
using NBInclude

struct P <: POMDP{Nothing, Nothing, Nothing} end
ParticleFilters.obs_weight(::P, ::Nothing, ::Nothing, ::Nothing, ::Nothing) = 1.0

@testset "implemented" begin
    @test @implemented obs_weight(::P, ::Nothing, ::Nothing, ::Nothing, ::Nothing)
    @test obs_weight(P(), nothing, nothing, nothing, nothing) == 1.0
end

include("example.jl")
include("domain_specific_resampler.jl")
include("beliefs.jl")

struct ContinuousPOMDP <: POMDP{Float64, Float64, Float64} end
@testset "infer" begin
    p = TigerPOMDP()
    filter = BootstrapFilter(p, 10000)
    Random.seed!(filter, 47)
    b = @inferred initialize_belief(filter, initialstate(p))
    @testset "sir" begin
        m = @inferred mode(b)
        m = @inferred mean(b)
        it = @inferred support(b)
        @inferred weighted_particles(b)
    end
    @testset "lowvar" begin
        rs = LowVarianceResampler(1000)
        @inferred resample(rs, b, MersenneTwister(3))
        ps = particles(b)
        ws = ones(length(ps))
        @inferred resample(rs, WeightedParticleBelief(ps, ws, sum(ws)), MersenneTwister(3))
        @inferred resample(rs, WeightedParticleBelief{Bool}(ps, ws, sum(ws), nothing), MersenneTwister(3))
    end
    # test that the special method for ParticleCollections works
    @testset "collection" begin
        rs = LowVarianceResampler(1000)
        b = ParticleCollection(1:1000)
        rb1 = @inferred resample(rs, b, MersenneTwister(3))
        rb2 = @inferred resample(rs, WeightedParticleBelief(particles(b), ones(n_particles(b))), MersenneTwister(3))
        @test all(particles(rb1).==particles(rb2))
    end

    @testset "unweighted" begin
        rng = MersenneTwister(47)
        uf = UnweightedParticleFilter(p, 1000, rng)
        ps = @inferred initialize_belief(uf, initialstate(p))
        a = @inferred rand(rng, actions(p))
        sp, o = @inferred @gen(:sp, :o)(p, rand(rng, ps), a, rng)
        bp = @inferred update(uf, ps, a, o)

        wp1 = @inferred collect(weighted_particles(ParticleCollection([1,2])))
        wp2 = @inferred collect(weighted_particles(WeightedParticleBelief([1,2], [0.5, 0.5])))
        @test wp1 == wp2
    end

    @testset "normal" begin
        pf = BootstrapFilter(ContinuousPOMDP(), 100)
        ps = @inferred initialize_belief(pf, Normal())
    end
end

struct TerminalPOMDP <: POMDP{Int, Int, Float64} end
POMDPs.isterminal(::TerminalPOMDP, s) = s == 1
POMDPs.observation(::TerminalPOMDP, a, sp) = Normal(sp)
POMDPs.transition(::TerminalPOMDP, s, a) = Deterministic(s+a)
@testset "pomdp terminal" begin
    pomdp =  TerminalPOMDP()
    pf = BootstrapFilter(pomdp, 100)
    bp = update(pf, initialize_belief(pf, Categorical([0.5, 0.5])), -1, 1.0)
    @test all(particles(bp) .== 1)
end

@testset "alpha" begin
    # test specific method for alpha vector policies and particle beliefs
    pomdp = BabyPOMDP()
    # the first two alphas were gotten from FIB.jl; the third is from always feeding (see #46)
    alphas = [[-16.0629, -36.5093], [-19.4557, -29.4557], [-50.0, -60.0]]
    amap = collect(ordered_actions(pomdp))
    push!(amap, true)
    policy = AlphaVectorPolicy(pomdp, alphas, amap)

    # initial belief is 100% confidence in baby being hungry
    b = ParticleCollection([true for i=1:100])

    # because baby is hungry, policy should feed (return true)
    @test action(policy, b) == true
    @test isapprox(value(policy, b), -29.4557)
end

is_ci = get(ENV, "CI", "false") == "true"
is_travis = get(ENV, "TRAVIS", "false") == "true"

@show is_ci
@show is_travis

@warn("Notebook smoke testing is disabled on JuliaCI. We should re-enable it asap")

if !is_ci || is_travis
    @testset "data series" begin
        cd("../notebooks") do
            @nbinclude("../notebooks/Filtering-a-Trajectory-or-Data-Series.ipynb")
        end
    end

    @testset "feedback" begin
        @nbinclude("../notebooks/Using-a-Particle-Filter-for-Feedback-Control.ipynb"; softscope=true)
    end

    @testset "pomdps" begin
        @nbinclude("../notebooks/Using-a-Particle-Filter-with-POMDPs-jl.ipynb")
    end
end


text:
#println(PyCall.libpython)
using Comrade
using FFTW
using Statistics
using Test
using Plots

include(joinpath(@__DIR__, "../test_util.jl"))

include(joinpath(@__DIR__, "observation.jl"))
include(joinpath(@__DIR__, "distributions.jl"))
include(joinpath(@__DIR__, "models.jl"))
include(joinpath(@__DIR__, "gains.jl"))
include(joinpath(@__DIR__, "io.jl"))
include(joinpath(@__DIR__, "bayes.jl"))


text:
"""
    TrajectoryOptimization
Primary module for setting up and solving trajectory optimization problems.
"""
module TrajectoryOptimization

using RigidBodyDynamics
using ForwardDiff
using DocStringExtensions
using Interpolations
using RecipesBase
using LinearAlgebra
using Statistics
using Random
using SparseArrays
using SuiteSparse
using StaticArrays
using Logging
using Formatting
using Plots
using BenchmarkTools
using PartedArrays
using Parameters
using Rotations
using BlockArrays
using MathOptInterface
using TimerOutputs

export
    Dynamics,
    Problems

# Primary types
export
    Model,
    QuadraticCost,
    LQRCost,
    LQRObjective,
    GenericCost,
    Trajectory

export
    Problem,
    iLQRSolver,
    iLQRSolverOptions,
    AugmentedLagrangianSolver,
    AugmentedLagrangianSolverOptions,
    AugmentedLagrangianProblem,
    ALTROSolverOptions,
    DIRCOLSolver,
    DIRCOLSolverOptions,
    ProjectedNewtonSolver,
    ProjectedNewtonSolverOptions,
    SequentialNewtonSolver,
    DIRTRELSolver,
    Discrete,
    Continuous,
    Constraint,
    BoundConstraint,
    Equality,
    Inequality,
    ConstraintSet,
    StageConstraintSet,
    TerminalConstraintSet,
    ConstraintSet,
    Objective,
    Constraints

export
    rk3,
    rk4,
    midpoint,
    add_constraints!,
    goal_constraint,
    initial_controls!,
    initial_state!,
    circle_constraint,
    sphere_constraint

# Primary methods
export
    solve,
    solve!,
    rollout!,
    rollout,
    forwardpass!,
    backwardpass!,
    cost,
    max_violation,
    max_violation_direct,
    infeasible_control,
    line_trajectory,
    evaluate!,
    jacobian!

export
    get_sizes,
    num_constraints,
    get_num_constraints,
    get_num_controls,
    init_results,
    to_array,
    get_N,
    to_dvecs,
    quat2rot,
    sphere_constraint,
    circle_constraint,
    plot_trajectory!,
    plot_vertical_lines!,
    convergence_rate,
    plot_obstacles,
    evals,
    reset,
    reset_evals,
    final_time,
    total_time,
    count_constraints,
    inequalities,
    equalities,
    bounds,
    labels,
    terminal,
    stage,
    interp_rows

# Trajectory Types
Trajectory{T} = Vector{T} where T <: AbstractArray
VectorTrajectory{T} = Vector{Vector{T}} where T <: Real
MatrixTrajectory{T} = Vector{Matrix{T}} where T <: Real
AbstractVectorTrajectory{T} = Vector{V} where {V <: AbstractVector{T}, T <: Real}
DiagonalTrajectory{T} = Vector{Diagonal{T,Vector{T}}} where T <: Real
PartedVecTrajectory{T} = Vector{PartedVector{T,Vector{T}}}
PartedMatTrajectory{T} = Vector{PartedMatrix{T,Matrix{T}}}

include("constraints.jl")
include("constraint_sets.jl")
include("cost.jl")
include("model.jl")
include("integration.jl")
include("utils.jl")
include("objective.jl")
include("problem.jl")
include("solvers.jl")
include("rollout.jl")
include("dynamics.jl")
include("problems.jl")
include("logger.jl")

write_ipopt_options()
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
@system Irradiance begin
    PFD ~ hold

    #HACK: should be PPFD from Radiation
    PPFD(PFD): photosynthetic_photon_flux_density ~ track(u"μmol/m^2/s")

    #FIXME: duplicate? already considered in Radiation?
    #FIXME: how is (1 - δ) related to α in EnergyBalance?
    # leaf reflectance + transmittance
    δ: leaf_scattering => 0.15 ~ preserve(parameter)
    f: leaf_spectral_correction => 0.15 ~ preserve(parameter)

    Ia(PPFD, δ): absorbed_irradiance => begin
        PPFD * (1 - δ)
    end ~ track(u"μmol/m^2/s" #= Quanta =#)

    I2(Ia, f): effective_irradiance => begin
        Ia * (1 - f) / 2 # useful light absorbed by PSII
    end ~ track(u"μmol/m^2/s" #= Quanta =#)
end


text:
module Wrapper
    # Disabled precompilation due to segmentation fault when calling functions
    # in Wrapper module -- with __init__() calling @initcxx
    __precompile__(false)

    import Base: iterate, IteratorSize, length

    using CxxWrap
    @wrapmodule(joinpath(@__DIR__, "..", "deps", "usr", "lib", "libpylon_julia_wrapper"), :define_pylon_wrapper)

    notify_async_cond(handle::Ptr{Nothing}) = ccall(:uv_async_send, Cint, (Ptr{Nothing}, ), handle)
    notify_async_cond_safe_c = @safe_cfunction(notify_async_cond, Int32, (Ptr{Nothing}, ))

    retrieve_result(camera::InstantCamera, timeoutMs::UInt32) = retrieve_result(camera, timeoutMs, TimeoutHandling_ThrowException)

    iterate(list::DeviceInfoList) = length(list) > 0 ? (list[1], 2) : nothing
    iterate(list::DeviceInfoList, i) = i <= length(list) ? (list[i], i+1) : nothing
    IteratorSize(list::DeviceInfoList) = length(list)
end


text:
using Reexport

@reexport module IntegrationUncertainData
    import ..CausalityTests
    import CausalityToolsBase: causality
    import UncertainData: resample, AbstractUncertainValue, AbstractUncertainValueDataset, ConstrainedResampling
    import CausalityToolsBase: CausalityTest#, causality

    include("integration_with_uncertaindata.jl")
    include("integration_with_uncertaindata_constrained_resampling.jl")
    include("integration_with_uncertaindata_using_schemes.jl")

    # We may specialize behaviour if needed, but defining `causality`
    # for the generic tests works well for now, so we don't need to
    # import these files. If specific behaviour is needed, just
    # activate the files below and write custom methods for each test.

    #include("integration_uncertaindata_ConvergentCrossMappingTest.jl")
    #include("integration_uncertaindata_CrossMappingTest.jl")
    #include("integration_uncertaindata_VisitationFrequencyTest.jl")
    #include("integration_uncertaindata_TransferOperatorGridTest.jl")
    #include("integration_uncertaindata_JointDistanceDistributionTest.jl")
    #include("integration_uncertaindata_JointDistanceDistributionTTest.jl")

    export causality

end


text:
##### Numeric simplification

"""
simplify(x; rewriter=default_simplifier(),
            threaded=false,
            polynorm=true,
            thread_subtree_cutoff=100)

Simplify an expression (`x`) by applying `rewriter` until there are no changes.
`polynorm=true` applies `polynormalize` in the beginning of each fixpoint iteration.
"""
function simplify(x;
                  polynorm=false,
                  threaded=false,
                  thread_subtree_cutoff=100,
                  rewriter=nothing)
    f = if rewriter === nothing
        if threaded
            threaded_simplifier(thread_subtree_cutoff)
        elseif polynorm
            serial_polynormal_simplifier
        else
            serial_simplifier
        end
    else
        Fixpoint(rewriter)
    end

    PassThrough(f)(x)
end

Base.@deprecate simplify(x, ctx; kwargs...)  simplify(x; rewriter=ctx, kwargs...)

"""
    substitute(expr, dict)

substitute any subexpression that matches a key in `dict` with
the corresponding value.
"""
function substitute(expr, dict; fold=true)
    haskey(dict, expr) && return dict[expr]

    if istree(expr)
        if fold
            canfold=true
            args = map(arguments(expr)) do x
                x′ = substitute(x, dict; fold=fold)
                canfold = canfold && !(x′ isa Symbolic)
                x′
            end
            canfold && return operation(expr)(args...)
            args
        else
            args = map(x->substitute(x, dict), arguments(expr))
        end
        similarterm(expr, operation(expr), args)
    else
        expr
    end
end


text:
using Thrift
using Base.Test

import Thrift.process, Thrift.meta

global srvr

# include the generated module, which in-turn includes our implementation code in `proto_tests_impl.jl`
# prepend with @everywhere when running a TProcessPoolServer on multiple Julia processors
#@everywhere include("gen-jl/proto_tests/proto_tests.jl")
isdefined(:srvcctrl) || include("gen-jl/srvcctrl/srvcctrl.jl");
isdefined(:proto_tests) || include("gen-jl/proto_tests/proto_tests.jl");
import proto_tests: ProtoTestsProcessor, ProtoTestsClient, InvalidOperation, AllTypes, AllTypesDefault, TestEnum
import proto_tests: test_hello, test_exception, test_oneway, ping, test_enum, test_types, test_types_default
import srvcctrl: start_service, stop_service

transport_factory(x) = x
protocol_factory(x) = TBinaryProtocol(x)
#protocol_factory(x) = TCompactProtocol(x)

function make_server()
    # create a server instance with our choice of protocol and transport
    srvr_processor = ProtoTestsProcessor()
    srvr_transport = TServerSocket(19999)

    #srvr = TProcessPoolServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)
    #srvr = TTaskServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)
    srvr = TSimpleServer(srvr_transport, srvr_processor, transport_factory, protocol_factory, transport_factory, protocol_factory)

    println("Transport: $(typeof(srvr_transport))")
    println("Protocol : $(typeof(protocol_factory(srvr_transport)))")
    println("Processor: $(typeof(srvr_processor))")
    println("Server   : $(typeof(srvr))")
    srvr
end

srvr = make_server()
# start serving client requests
println("Starting to serve requests...")

if haskey(ENV, "TEST_SRVR_ASYNC")
    @async try
        start_service(srvr)
    catch ex
        println("Server stopped with $ex")
    finally
        println("Stopped serving requests.")
    end
else
    start_service(srvr)
end


text:
using Test     #src
# # Discrete-time Controlled Invariant Set
#
#md # [![Binder](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/generated/discrete_controlled.ipynb)
#md # [![nbviewer](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/discrete_controlled.ipynb)
#
# An similar example is available in a [codeocean capsule](https://doi.org/10.24433/CO.6396918.v1).
#
# ## Introduction
#
# This example considers the linear control system already introduced in the Example 2 of [LTJ20]; see also Example 2 of [LRJ20]:
# ```math
# \begin{aligned}
# x_{k+1} & = x_k + u_k / 2\\
# u_{k+1} & = u_k'
# \end{aligned}
# ```
# with $(x_k, u_k) \in [-1, 1]^2$.
#
# The system is $x_{k+1} = Ax_k + Bu_k$ where $B = (0, 1)$ and
# ```math
# A = \begin{bmatrix}
# 1 & 1/2\\
# 0 & 0
# \end{bmatrix}.
# ```
# As shown in Example 4 of [LTJ18], a set $S \subseteq [-1, 1]^2$ is controlled invariant for this system if
# ```math
# \begin{bmatrix}
# 1 & 1/2
# \end{bmatrix}
# S \subseteq
# \begin{bmatrix}
# 1 & 0
# \end{bmatrix}
# S
# ```
#
# [LRJ20] Legat, Benoît, Saša V. Raković, and Raphaël M. Jungers.
# *Piecewise semi-ellipsoidal control invariant sets.*
# IEEE Control Systems Letters 5.3 (2020): 755-760.
#
# [LTJ18] B. Legat, P. Tabuada and R. M. Jungers.
# *Computing controlled invariant sets for hybrid systems with applications to model-predictive control*.
# 6th IFAC Conference on Analysis and Design of Hybrid Systems ADHS (2018).

# We need to pick an LP and an SDP solver, see [here](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) for a list of available ones.

using SetProg
import GLPK
lp_solver = optimizer_with_attributes(GLPK.Optimizer, MOI.Silent() => true)
import CSDP
sdp_solver = optimizer_with_attributes(CSDP.Optimizer, MOI.Silent() => true)
A = [1.0 0.5]
E = [1.0 0.0]
using Polyhedra
lib = DefaultLibrary{Float64}(lp_solver)
h = HalfSpace([1, 0], 1.0) ∩ HalfSpace([-1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([0, -1], 1)
□ = polyhedron(h, lib) # [0, 1]^2
v = convexhull([1.0, 0], [0, 1], [-1, 0], [0, -1])
◇ = polyhedron(v, lib); # polar of [0, 1]^2

# ## Polyhedral template
#
# ### Fixed point approach
#
# This section shows that the maximal control invariant set of this simple control system is polyhedral. Moreover, this polyhedron is obtained after only one fixed point iteration of the standard viability
# kernel algorithm. We implement the fixed point iteration with the following function:

function fixed_point_iteration(set::Polyhedron)
    new_set = set ∩ (A \ (E * set))
    removehredundancy!(new_set)
    return new_set
end

# We start with $[-1, 1]^2$ and obtain the polytope `mci` after one iteration.

mci = fixed_point_iteration(□)

# One additional iteration gives the same polytope, showing that this polytope is control invariant.

fixed_point_iteration(mci)

# We plot this polytope `mci` in yellow below along with the safe set $[-1, 1]^2$ in green.

using Plots
function hexcolor(rgb::UInt32)
    r = ((0xff0000 & rgb) >> 16) / 255
    g = ((0x00ff00 & rgb) >>  8) / 255
    b = ((0x0000ff & rgb)      ) / 255
    Plots.RGBA(r, g, b)
end # Values taken from http://www.toutes-les-couleurs.com/code-couleur-rvb.php
lichen = hexcolor(0x85c17e)
canard = hexcolor(0x048b9a)
aurore = hexcolor(0xffcb60)

plot(ratio=:equal, tickfont=Plots.font(12))
plot!(□, color=lichen)
plot!(mci, color=aurore)

# We also plot their respective polars. Note that the inclusion is reversed in the polar space.

polar_mci = polar(mci)
plot(ratio=:equal, tickfont=Plots.font(12))
plot!(polar_mci, color=aurore)
plot!(◇, color=lichen)

function maximal_invariant(template, heuristic::Function)
    solver = if template isa Polytope
        lp_solver
    else
        sdp_solver
    end
    model = Model(solver)
    @variable(model, S, template)
    @constraint(model, S ⊆ □)
    @constraint(model, A * S ⊆ E * S)
    @objective(model, Max, heuristic(volume(S)))
    optimize!(model)
    @show solve_time(model)
    @show termination_status(model)
    @show objective_value(model)
    return value(S)
end
function primal_plot(set; npoints=256, args...)
    plot(ratio=:equal, tickfont=Plots.font(12); args...)
    plot!(□, color=lichen)
    plot!(mci, color=aurore)
    plot!(set, color=canard, npoints=npoints)
end
function polar_plot(set; npoints=256, args...)
    plot(ratio=:equal, tickfont=Plots.font(12); args...)
    plot!(SetProg.Sets.polar(set), color=canard, npoints=npoints)
    plot!(polar_mci, color=aurore)
    plot!(◇, color=lichen)
end

# ### Linear programming approach
#
# As introduced in [R21], given a fixed polyhedral conic partition of the state-space,
# we can search over all polyhedra for which the gauge or support function
# is linear over each piece of the partition.
# The partition can be defined by a polytope as described in [Eq. (4.6), R21].
# In this example, because the sets are represented with the support function,
# the pieces will be used for the support function/polar set and
# the pieces of the gauge function/primal will depend on the value of
# the support function on each piece.
#
# [R21] Raković, S. V.
# *Control Minkowski–Lyapunov functions*
# Automatica, Elsevier BV, 2021, 128, 109598

# Let's start with setting the partition using the cross polytope `◇`.

sol_polytope_◇ = maximal_invariant(Polytope(symmetric=true, piecewise=◇), L1_heuristic)
SetProg.Sets.print_support_function(sol_polytope_◇)

# We can see that the polar set is unbounded so the primal set has zero volume.
# More precisely, the polar is $[-1, 1] \times [-\infty, \infty]$ and the primal
# set is $[-1, 1] \times \{0\}$.

# Let's now try setting the partition using the square `□`.

sol_polytope_□ = maximal_invariant(Polytope(symmetric=true, piecewise=□), L1_heuristic)
SetProg.Sets.print_support_function(sol_polytope_□)

# The primal plot is below:

primal_plot(sol_polytope_□, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# And the polar plot is below:

polar_plot(sol_polytope_□, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# Let's now use a partition with 8 pieces.

pieces8 = convexhull(□, 1.25 * ◇)
sol_polytope_8 = maximal_invariant(Polytope(symmetric=true, piecewise=pieces8), L1_heuristic)
SetProg.Sets.print_support_function(sol_polytope_8)

# The primal plot is below:

primal_plot(sol_polytope_8, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# And the polar plot is below:

polar_plot(sol_polytope_8, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# We can use the polar of the polytope resulting from the first fixed point iteration to generate a refined conic partition.
# With this partition, the maximal piecewise semi-ellipsoidal control invariant set matches the maximal control invariant set.

sol_polytope_mci = maximal_invariant(Polytope(symmetric=true, piecewise=polar_mci), L1_heuristic)
SetProg.Sets.print_support_function(sol_polytope_mci)

# The primal plot is below:

primal_plot(sol_polytope_mci)

# And the polar plot is below:

polar_plot(sol_polytope_mci)

# ## Ellipsoidal template
#
# We now compute the maximal ellipsoidal control invariant set. Here we can either maximize its volume (which corresponds to $\log(\det(Q))$ or $\sqrt[n]{\det(Q)}$ in the objective function) or minimize the sum of the squares of the length of its semi-axes of the polar (which corresponds to the trace of $Q$ in the objective function).

# We can see below that the control invariant set of maximal volume has support function
# $$h^2(S, x) = x_1^2 - \frac{1}{2} x_1x_2 + x_2^2.$$

sol_ell_vol = maximal_invariant(Ellipsoid(symmetric=true, dimension=2), nth_root)
SetProg.Sets.print_support_function(sol_ell_vol)

# We can see below the ellipsoid in blue along with the maximal control invariant set in yellow and the safe set in green.

primal_plot(sol_ell_vol)

# Below is are the corresponding sets in the polar space.

polar_plot(sol_ell_vol)

# We can see below that the control invariant set of minimal sum of squares of the length of the semi-axes of the polar has support function
# $$h^2(S, x) = x_1^2 - \alpha x_1x_2 + x_2^2.$$
# with $\alpha \approx 1.377456$.
# Note that the sum of the squares of the length of the semi-axes of the polar is equal to the integral of $h^2(S, x)$ over the hypercube $[-1, 1]^2$ hence the use of the `L1_heuristic(vol, ones(2))` as objective which means taking the integral over the hyperrectangle with vertex `ones(2) = [1, 1]`.

sol_ell_L1 = maximal_invariant(Ellipsoid(symmetric=true, dimension=2), vol -> L1_heuristic(vol, ones(2)))
SetProg.Sets.print_support_function(sol_ell_L1)

# The primal plot is below:

primal_plot(sol_ell_L1)

# And the polar plot is below:

polar_plot(sol_ell_L1)

# ## Piecewise semi-ellipsoidal template
#
# We now study the maximal piecewise semi-ellipsoidal control invariant sets of a given conic partition.
# The volume is not directly maximized. Instead, for each cone, we compute the sum $s$ of the normalized rays and consider the polytope obtained by intersecting the cone with the halfspace $s^\top x \le \|s\|_2^2$. We integrate the quadratic form corresponding to this cone, i.e. $h^2(S, x)$ over the polytope. The sum of the integrals over each polytope is the objective function we use. This can be seen as the generalization of the sum of the squares of the semi-axes of the polar of the ellipsoid.
#
# Note that the constraint (29) of Program 1 of [LRJ20] is implemented with Proposition 2 of [LRJ20] for all results of this example.

# The conic partition are obtained by considering the conic hull of each facets of a given polytope.
# We first consider the conic partition corresponding to the polar of the safe set $[-1, 1]^2$. This gives the four quadrants as cones of the conic partition.
# The maximal piecewise semi-ellipsoidal control invariant set with this partition has the following support function:
# $$
# h^2(S, x) = \begin{cases}
# (x_1 - x_2)^2 & \text{ if }x_1x_2 \le 0,\\
# x_1^2 - x_1x_2/2 + x_2^2 & \text{ if }x_1x_2 \ge 0.
# \end{cases}
# $$
# For the cones $x_1x_2 \le 0$, the semi-ellipsoid matches the safe set and the maximal control invariant set.
# For the cones $x_1x_2 \ge 0$, the semi-ellipsoid matches the maximal volume control invariant ellipsoid.
# This illustrates one key feature of piecewise semi-ellipsoidal sets, they can combines advantages of both polyhedra and ellipsoids. It can be polyhedral on the directions where the maximal control invariant set is polyhedral and be ellipsoidal on the directions where the maximal control invariant set is smooth or requires many halfspaces in its representation.

sol_piece_◇ = maximal_invariant(Ellipsoid(symmetric=true, piecewise=◇), L1_heuristic)
SetProg.Sets.print_support_function(sol_piece_◇)

# The primal plot is below:

primal_plot(sol_piece_◇, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# And the polar plot is below:

polar_plot(sol_piece_◇, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# Let's now try setting the partition with the square `□`.

sol_piece_□ = maximal_invariant(Ellipsoid(symmetric=true, piecewise=□), L1_heuristic)
SetProg.Sets.print_support_function(sol_piece_□)

# The primal plot is below:

primal_plot(sol_piece_□, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# And the polar plot is below:

polar_plot(sol_piece_□, xlim=(-1.6, 1.6), ylim=(-1.6, 1.6))

# Let's now try with the partition of 8 pieces.

sol_piece_8 = maximal_invariant(Ellipsoid(symmetric=true, piecewise=pieces8), L1_heuristic)
SetProg.Sets.print_support_function(sol_piece_8)

# The primal plot is below:

primal_plot(sol_piece_8, xlim=(-1.05, 1.05), ylim=(-1.05, 1.05))

# And the polar plot is below:

polar_plot(sol_piece_8, xlim=(-1.1, 1.1), ylim=(-1.05, 1.05))

# We can use the polar of the polytope resulting from the first fixed point iteration to generate a refined conic partition.
# With this partition, the maximal piecewise semi-ellipsoidal control invariant set matches the maximal control invariant set. The support function is
# $$
# h^2(S, x) = \begin{cases}
# (x_1 - x_2)^2 & \text{ if }x_1x_2 \le 0,\\
# x_1^2 & \text{ if }x_2(x_1-2x_2) \ge 0,\\
# (x_1/2 + x_2)^2 & \text{ if }x_1(2x_2-x_1) \ge 0.\\
# \end{cases}
# $$

sol_piece_mci = maximal_invariant(Ellipsoid(symmetric=true, piecewise=polar_mci), L1_heuristic)
SetProg.Sets.print_support_function(sol_piece_mci)

# The primal plot is below:

primal_plot(sol_piece_mci)

# And the polar plot is below:

polar_plot(sol_piece_mci)

# ## Polyset template
#
# We now use the polyset templates. As details in [LRJ20], for inclusion
# constraints of the form `A * S ⊆ E * S`, the polar representation is used
# so we need the set to be convex. For this, we set the argument `convex=true`.
# We directly try `degree=4` because `degree=2` would simply give the same as
# `sol_ell_L1`.

sol4 = maximal_invariant(PolySet(symmetric=true, degree=4, convex=true), vol -> L1_heuristic(vol, ones(2)))
SetProg.Sets.print_support_function(sol4)

# The primal plot is below:

primal_plot(sol4)

# And the polar plot is below:

polar_plot(sol4)

# Let's now try degree 6.

sol6 = maximal_invariant(PolySet(symmetric=true, degree=6, convex=true), vol -> L1_heuristic(vol, ones(2)))
SetProg.Sets.print_support_function(sol6)

# The primal plot is below:

primal_plot(sol6)

# And the polar plot is below:

polar_plot(sol6)