text:
using Rectangle
using Test

@test Any[] == detect_ambiguities(Base, Core, Rectangle)

include("bsttest.jl")
include("interval.jl")
include("xytest.jl")

@testset "Utils" begin
    @test pcTol(Int) == 0
    @test pcTol(Rational) == 0
    @test pcTol(Float32) == 1f-3
    @test pcTol(Float64) == 1e-6
end
@testset "Rectangle" begin
    @test Rect(0.0, 0.0, 10, 10) == Rect(0, 0, 10, 10)
    r = Rect(0, 5, 10, 15)
    @test r == Rect(lx(r), ly(r), rx(r), ry(r))
    @test Rect([0 10; 0 10]) == Rect(0, 0, 10, 10)
    @test Rect{Float32}(Rect(1, 2, 3, 4)) == Rect(1f0, 2f0, 3f0, 4f0)
    @test string(Rect(0, 0, 10, 10)) == "Rect:[0 0 10 10]"
    @test union(Rect(0, 0, 10, 5), Rect(11, -1, 12, 13)) == Rect(0, -1, 12, 13)
    @test intersect(Rect(0.0, 0, 10, 5), Rect(11, -1, 12, 13)) == nothing
    @test intersect(Rect(0.0, 0.0, 5.0, 2.0),
                    Rect(4.0, 2.0, 10.0, 10.0)) == nothing
    @test inside((1, 2),   Rect(0,   0, 10, 10))
    @test inside((1.0, 2), Rect(0f0, 0, 10, 10))
    @test inside(Rect(1, 1, 3, 3), Rect(0, 0, 10, 10))
    @test Rect(0, 0, 10, 10) == Rect(0.0000001, 0.0000001, 10.0, 10.0)
    @test Rect(0.0, 0.0, 10.0, 10.0) != Rect(0.000002, 0.0000001, 10.0, 10.0)
    @test perimeter(Rect(0, 0, 1, 10)) == 22
    @test area(Rect(0, 0, 1, 10))      == 10
    @test perimeter(Rect(0, 0, 1//2, 1)) == 3

    @test x(Rect(0, 5, 10, 10)) == (0, 10)
    @test y(Rect(0, 5, 10, 10)) == (5, 10)

    @test  has_x_overlap(Rect(0.0, 0.0, 10.0, 1.0), Rect(3.0, 2.0, 11.0, 4.0))
    @test !has_y_overlap(Rect(0.0, 0.0, 10.0, 1.0), Rect(3.0, 2.0, 11.0, 4.0))

    @test  has_y_overlap(Rect(0.0, 0.0, 1.0, 10.0), Rect(2.0, 3.0, 4.0, 11.0))
    @test !has_x_overlap(Rect(0.0, 0.0, 1.0, 10.0), Rect(2.0, 3.0, 4.0, 11.0))

    @test visibleX(Rect(0.0, 0.0, 10.0, 1.0), Rect(3.0, 2.0, 11.0, 4.0)) ==
        Rect(3.0, 1.0, 10.0, 2.0)
    @test visibleY(Rect(0.0, 0.0, 10.0, 1.0), Rect(3.0, 2.0, 11.0, 4.0)) ==
            nothing
    @test visibleY(Rect(0.0, 0.0, 10.0, 1.0), Rect(10.0, 0.0, 11.0, 4.0)) ==
            Line(10.0, 0.0, 10.0, 1.0)

    @test visibleY(Rect(0.0, 0.0, 1.0, 10.0), Rect(2, 3, 4, 11)) ==
        Rect(1.0, 3.0, 2.0, 10.0)
    @test visibleX(Rect(0.0, 0.0, 1.0, 10.0), Rect(2.0, 3.0, 4.0, 11.0)) ==
        nothing

    @test begin
        r = (8.928571428571429, 4.583333333333334)
        v = avg_min_dist(Rect(0,0,10,10), Rect(15,5,20,20))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end

    @test begin
        r = (4.583333333333334, 8.928571428571429)
        v = avg_min_dist(Rect(0,0,10,10), Rect(5,15,20,20))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end

    @test begin
        r = (0.8333333333333334, 8.928571428571429)
        v = avg_min_dist(Rect(0,0,10,10), Rect(0,15,15,20))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end

    @test begin
        r = (0.0, 9.166666666666668)
        v = avg_min_dist(Rect(0,0,10,10), Rect(0,15,10,20))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end

    @test begin
        v = projectX(Rect(0,0,10,10), Rect(0,15.0,10,20))
        @assert v[1] == (nothing, nothing)
        @assert v[2] == (Rect(0.0,0.0,10.0,10.0), Rect(0.0,15.0,10.0,20.0))
        v[3] == (nothing, nothing)
    end

    @test begin
        r = (7.0, 7.0)
        v = avg_min_dist(Rect(0,0,6,6), Rect(10,10,16,16))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end

    @test begin
        r = (0, 5)
        v = min_dist(Rect(0,0,10,10.0), Rect(0,15,15,20))
        abs(v[1] - r[1]) < 1e-6 && abs(v[2] - r[2]) < 1e-6
    end
    @test min_dist(Rect(0,0,10,10), Rect(11,15,15,20)) == (1, 5)

    @test to_plot_shape(Rect(1, 2, 3, 4)) == ([1, 3, 3, 1], [2, 2, 4, 4])

    rects = [Rect(0, 20*i, 100, 10 + 20*i) for i = 1:10]
    is = [i for i=1:10]
    ormy = create_ordered_map(rects, is, dir=2)
    v = insert_rect!(ormy, Rect(0, 20, 100, 30), 11)
    @test delete_rect!(ormy, Rect(0, 20, 100, 30)) == 11
    @test intersect(ormy, Rect(0, 0, 110, 55)) == ([Rect(0, 40, 100, 50)], [2])

    rects = [Rect(0, 20*i, 100.0, 10 + 20*i) for i = 1:10]
    is = [i for i=1:10]
    @test_throws AssertionError create_ordered_map(rects, is, dir=2, reverseMax=-1)
    ormy = create_ordered_map(rects, is, dir=2, reverseMax=400)
    v = insert_rect!(ormy, Rect(0, 20, 100, 30), 11)
    @test delete_rect!(ormy, Rect(0, 20, 100, 30)) == 11
    @test intersect(ormy, Rect(0, 0, 110, 75)) ==
        ([Rect(0, 60, 100, 70), Rect(0, 40, 100, 50)], [3, 2])

    ormx = create_ordered_map(rects, is, dir=1)
    v = insert_rect!(ormx, Rect(0, 20, 100, 30), 11)
    @test delete_rect!(ormx, Rect(0, 20, 100, 30)) == 11
    @test intersect(ormx, Rect(0, 0, 110, 55)) == ([Rect(0, 40, 100, 50)], [2])

    begin
        rs = [Rect(0, 0, 10, 10), Rect(11, 11, 20, 20),
              Rect(21, 21, 30, 30), Rect(31, 31, 40, 40)]
        ormx = create_ordered_map(rs, fill(1, 4))
        @test intersect(ormx, Rect(0, 0, 10, 10), 1, 1) ==
            ([Rect(0, 0, 10, 10), Rect(11, 11, 20, 20),
              Rect(21, 21, 30, 30), Rect(31, 31, 40, 40)], fill(1, 4))

        @test intersect(ormx, Rect(0, 0, 10, 10), 1, 0) ==
            ([Rect(0, 0, 10, 10)], fill(1, 1))
        @test intersect(ormx, Rect(1, 11, 9, 19), 1, 0) ==
            (Rect{Int}[], Nothing[])
    end

    @test hlines(Rect(0, 0, 10, 10)) == [Line(0, 0, 10, 0), Line(0, 10, 10, 10)]
    @test vlines(Rect(0, 0, 10, 10)) == [Line(0, 0, 0, 10), Line(10, 0, 10, 10)]
    @test lines(Rect(0, 0, 10, 10)) == [Line(0, 0, 10, 0), Line(0, 10, 10, 10),
                                        Line(0, 0, 0, 10), Line(10, 0, 10, 10)]
    @test olines(Rect(0, 0, 10, 10)) == [Line(0, 0, 10, 0), Line(10, 0, 10, 10),
                                         Line(10, 10, 0, 10), Line(0, 10, 0, 0)]

    @test cg(Rect(0, 0, 10, 10)) == (5, 5)
    @test intersects(Rect(0, 0, 10, 10), Rect(5, 5, 15, 15))
    @test intersects(Rect(0, 0, 10, 10.0), Line(-1, -1, 11.0, 11))
    @test intersects(Rect(0, 0, 10, 10), Line(5, -5, 11, 11))
    @test !intersects(Rect(0, 0, 10, 10), Line(11, 0, 11, 11))

    hls = [Line(0, 5, 10, 5), Line(0, 4, 10, 4), Line(0, 3, 10, 3),
           Line(0, 2, 10, 2), Line(0, 1, 10, 1), Line(0, 0, 10, 0f0)]

    @test hline_xsection(Rect(1, 1.5, 9, 4.5), hls) == [2, 3, 4]

    vls = [Line(0, 0, 0, 10), Line(1, 0, 1, 10), Line(2, 0, 2, 10),
           Line(3, 0, 3, 10), Line(4, 0, 4, 10), Line(5, 0, 5, 10f0)]

    @test vline_xsection(Rect(1.5, 1, 4.5, 9), vls) == [3, 4, 5]
end

@testset "Line" begin
    @test begin
        l = Line(0, 0, 10, 20)
        start(l) == (start(l, 1), start(l, 2))
    end
    @test begin
        l = Line(0, 0, 10, 20)
        endof(l) == (endof(l, 1), endof(l, 2))
    end
    @test Line([0 10; 0 10]) == Line(0.0, 0, 10, 10)
    @test string(Line(0, 0, 10, 10)) == "Line:[0 0 10 10]"
    @test isHorizontal(Line(0.0, 0, 10, 0))
    @test !isHorizontal(Line(0.0, 0, 10, 1))
    @test isVertical(Line(10, 0, 10, 10))
    @test !isVertical(Line(10, 0, 11.0, 10.0))
    @test parallelogram_area((0, 0), (10, 0), (0, 20)) == 200
    @test length(Line(0, 0, 3, 4)) == 5.0
    @test Rectangle.ratio(Line(0, 0, 5, 10.0), (1, 2)) == 0.2
    @test Rectangle.ratio(Line(0, 0, 5, 10.0), (2, 3)) == nothing
    @test Rectangle.ratio(Line(0.0, 0, 0, 10), (0, 5)) == 0.5
    @test Rectangle.ratio(Line(0.0, 0, 0, 10), (1, 5)) == nothing
    @test intersects(Line(0, 0, 10, 10), Line(10, 0, 0.0, 10))
    @test intersects(Line(0, 0, 10, 10), Line(5, 5, 0.0, 10))
    @test !intersects(Line(0, 0, 4, 4), Line(5, 5, 0.0, 10))
    @test !intersects(Line(0, 0, 4, 4), Line(10, 0, 0.0, 10))
    @test intersects(Line(0, 0, 10, 10.0), Line(-1, -1, 11.0, 11))
    @test intersects(Line(-1, -1, 11, 11.0), Line(0, 0, 10.0, 10))
    @test reverse(Line(0, 0, 10, 20)) == Line(10, 20, 0, 0)
    @test div(Line(0, 0, 10, 10), 1//20) == (1//2, 1//2)
    @test begin
        l = Line(1, 2, 3, 4)
        l == Line(sx(l), sy(l), ex(l), ey(l))
    end
    @test begin
        l = Line(1, 2, 3, 4)
        start(l) == (sx(l), sy(l)) && endof(l) == (ex(l), ey(l))
    end
    @test begin
        a = [Line(0, 0, 10, 0), Line(10, 0, 20, 0), Line(20, 0, 30, 0), Line(40, 0, 50, 0),
             Line(40, 5, 50, 5),
             Line(0, 10, 10, 10), Line(10, 10, 20, 10), Line(20, 10, 30, 10)]
        merge_axis_aligned(a) == [Line(0, 0, 30, 0), Line(40, 0, 50, 0),
                                  Line(40, 5, 50, 5), Line(0, 10, 30, 10)]
    end
    @test begin
        a = [Line(0, 0, 0, 10), Line(0, 10, 0, 20), Line(0, 20, 0, 30), Line(0, 40, 0, 50),
             Line(5, 40, 5, 50),
             Line(10, 0, 10, 10), Line(10, 10, 10, 20), Line(10, 20, 10, 30)]
        merge_axis_aligned(a, 2) == [Line(0, 0, 0, 30), Line(0, 40, 0, 50),
                                     Line(5, 40, 5, 50), Line(10, 0, 10, 30)]
    end
    @test begin
        a = [Line(0, 40, 0, 50), Line(0, 20, 0, 30), Line(0, 10, 0, 20), Line(0, 0, 0, 10),
             Line(5, 40, 5, 50),
             Line(10, 20, 10, 30), Line(10, 10, 10, 20), Line(10, 0, 10, 10)]
        merge_axis_aligned(a, 2, :decreasing) == [Line(0, 40, 0, 50), Line(0, 0, 0, 30),
                                                  Line(5, 40, 5, 50), Line(10, 0, 10, 30)]
    end
    @test begin
        a = [Line(0, 40, 0, 50), Line(0, 20, 0, 30f0), Line(0, 10, 0, 20), Line(0, 0, 0, 10),
             Line(5, 40, 5, 50),
             Line(10, 20, 10, 30), Line(10, 10, 10, 20), Line(10, 0, 10, 10)]
        merge_axis_aligned(a, 2, :decreasing) == [Line(0, 40, 0, 50f0), Line(0, 0, 0, 30),
                                                  Line(5, 40, 5, 50), Line(10, 0, 10, 30)]
    end
    @test intersect_axis_aligned(Line(0, 0, 10, 0), Line(0, 1, 0, 10), 1) == [0, 0]
    @test intersect_axis_aligned(Line(10, 0, 0, 0), Line(0, 10, 0, 1), 1) == [0, 0]
    @test intersect_axis_aligned(Line(0, 0, 10, 0), Line(0, 1, 0, 10), 1.0) == [0, 0]
    @test intersect_axis_aligned(Line(0, 0, 10, 0), Line(0, 2, 0, 10), 1.0) == []

    @test !horiz_desc(Line(0, 0, 10, 0), Line(0, 0, 10.0, 0.0))
    @test  horiz_desc(Line(0, 0, 10, 0), Line(1, 0, 10, 0))
    @test  horiz_desc(Line(0, 1, 10, 1), Line(1, 0, 10, 0))
    @test  horiz_desc(Line(0, 1f-4, 10, 1f-4), Line(1, 0, 10, 0))
    @test !horiz_desc(Line(2, 1f-4, 10, 1f-4), Line(1, 0, 10, 0))

    @test !vert_asc(Line(0, 0, 0, 10), Line(0, 0, 0.0, 10.0))
    @test  vert_asc(Line(0, 1, 0, 11), Line(0, 0, 0, 10))
    @test  vert_asc(Line(0, 0, 0, 10), Line(1, 0, 1, 8))
    @test  vert_asc(Line(1f-4, 0, 0, 11), Line(0, 0, 0, 10))
    @test !vert_asc(Line(1f-4, 0, 0, 9), Line(0, 0, 0, 10))
end


text:
cd(@__DIR__); include("setups/grid23x22.jl")
gr(dpi = 200)

## build rNGWF
rNGWF, dic_l2x = rngwf_all_vectors(D, 𝚽; σ = 0.2 * maximum(D), thres = 0.15)
rNGWF_dual = (rNGWF * rNGWF') \ rNGWF
Γ = rngwf_lx(dic_l2x)
f = digit_img[:]

important_idx = sortperm((rNGWF' * f).^2; rev = true)
red_box_inds = [41, 43, 66, 69, 71, 73, 75, 84, 85, 86, 89, 99]
plot(layout = Plots.grid(10,10), size = (2300, 2200))
for i = 1:100
    plot!(clims = (-0.002, 0.02))
    grid_vector_plot!(important_idx[i], i, rNGWF)
    if i in red_box_inds
        plot_square!(Nx, Ny; subplot = i)
    end
end
plt = current()
savefig(plt, "../figs/Grid23x22_fdigit3_rngwf_top100.png")

ext: jl


text:
import Literate
using Literate: script, markdown, notebook

# source directory for model files written using Literate.jl format
source_dir = joinpath(@__DIR__, "..", "examples")

# target directory for the tests files (bare .jl files)
target_dir_jl = joinpath(@__DIR__, "..", "test", "models")

# target directory for the markdown files (used by Documenter)
target_dir_md = joinpath(@__DIR__, "src", "models")
mkpath(target_dir_md)

# model files in sub-directories of source_dir
MODELS = ["Brusselator/Brusselator.jl",
          "Building/Building.jl",
          "DuffingOscillator/DuffingOscillator.jl",
          "EMBrake/EMBrake.jl",
          "Heat3D/Heat3D.jl",
          "ISS/ISS.jl",
          "LaubLoomis/LaubLoomis.jl",
          "Lorenz/Lorenz.jl",
          "LotkaVolterra/LotkaVolterra.jl",
          #"LotkaVolterraTangential/LotkaVolterraTangential.jl",
          "OpAmp/OpAmp.jl",
          "SquareWaveOscillator/SquareWaveOscillator.jl",
          "Platoon/Platoon.jl",
          "ProductionDestruction/ProductionDestruction.jl",
          "Quadrotor/Quadrotor.jl",
          "SEIR/SEIR.jl",
          "Spacecraft/Spacecraft.jl",
          "TransmissionLine/TransmissionLine.jl",
          "VanDerPol/VanDerPol.jl"]

for file in MODELS
    source_path = abspath(joinpath(source_dir, file))
    text = script(source_path, target_dir_jl, credit=false)
    code = strip(read(text, String))
    mdpost(str) = replace(str, "@__CODE__" => code)
    markdown(source_path, target_dir_md, postprocess=mdpost, credit=false)
    notebook(source_path, target_dir_md, execute=true, credit=false)
end


text:
module FusionFormulary

include("units.jl")
include("constants.jl")

include("plasma.jl")

end # module


text:
"""
For an initial guess from the index of a Hopf bifurcation point located in ContResult.bifpoint, returns a point which will be refined using `newtonHopf`.
"""
function HopfPoint(br::ContResult, index::Int64)
	@assert br.bifpoint[index].type == :hopf "The provided index does not refer to a Hopf point"
	bifpoint = br.bifpoint[index]							# Hopf point
	eigRes   = br.eig										# eigenvector at the Hopf point
	p = bifpoint.param										# parameter value at the Hopf point
	ω = (imag(eigRes[bifpoint.idx].eigenvals[bifpoint.ind_bif]))	# frequency at the Hopf point
	return BorderedArray(bifpoint.x, [p, ω] )
end

@with_kw struct HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl <: Lens, vectype, S <: AbstractLinearSolver, Sbd <: AbstractBorderedLinearSolver, Sbda <: AbstractBorderedLinearSolver}
	F::TF 								# Function F(x, p) = 0
	J::TJ 								# Jacobian of F wrt x
	Jadjoint::TJa						# Adjoint of the Jacobian of F
	d2F::Td2f = nothing					# Hessian of F
	lens::Tl							# parameter axis for the Hopf point
	a::vectype							# close to null vector of (J - iω I)^*
	b::vectype							# close to null vector of J - iω I
	linsolver::S						# linear solver
	linbdsolver::Sbd					# linear bordered solver
	linbdsolverAdjoint::Sbda			# linear bordered solver for the jacobian adjoint
end

hasHessian(pb::HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl, vectype, S, Sbd, Sbda}) where {TF, TJ, TJa, Td2f, Tl, vectype, S, Sbd, Sbda} = Td2f != Nothing

hasAdjoint(pb::HopfProblemMinimallyAugmented{TF, TJ, TJa, Td2f, Tl, vectype, S, Sbd, Sbda}) where {TF, TJ, TJa, Td2f, Tl, vectype, S, Sbd, Sbda} = TJa != Nothing

HopfProblemMinimallyAugmented(F, J, Ja, d2F, lens::Lens, a, b, linsolve) = HopfProblemMinimallyAugmented(F, J, Ja, d2F, lens, a, b, linsolve, BorderingBLS(linsolve), BorderingBLS(linsolve))

HopfProblemMinimallyAugmented(F, J, Ja, lens::Lens, a, b, linsolve) = HopfProblemMinimallyAugmented(F, J, Ja, nothing, lens, a, b, linsolve)

function (hp::HopfProblemMinimallyAugmented)(x, p::T, ω::T, par) where {T}
	# These are the equations of the minimally augmented (MA) formulation of the Hopf bifurcation point
	# input:
	# - x guess for the point at which the jacobian has a purely imaginary eigenvalue
	# - p guess for the parameter for which the jacobian has a purely imaginary eigenvalue
	# The jacobian of the MA problem is solved with a bordering method
	a = hp.a
	b = hp.b

	# we solve (J+iω)v + a σ1 = 0 with <b, v> = n
	n = T(1)
	v, σ1, flag, it = hp.linbdsolver(hp.J(x, set(par, hp.lens, p)),
							a, b,
							T(0), zero(x), n; shift = Complex{T}(0, ω))

	# we solve (J+iω)'w + b σ2 = 0 with <a, w> = n
	# we find sigma2 = conj(sigma1)
	# w, σ2, _ = fp.linbdsolver(fp.Jadjoint(x, p) - Complex(0, ω) * I, b, a, 0., zeros(N), n)

	# the constraint is σ = <w, Jv> / n
	# σ = -dot(w, apply(fp.J(x, p) + Complex(0, ω) * I, v)) / n
	# we should have σ = σ1

	return hp.F(x, set(par, hp.lens, p)), real(σ1), imag(σ1)
end

function (hopfpb::HopfProblemMinimallyAugmented)(x::BorderedArray, params)
	res = hopfpb(x.u, x.p[1], x.p[2], params)
	return BorderedArray(res[1], [res[2], res[3]])
end

# Struct to invert the jacobian of the Hopf MA problem.
struct HopfLinearSolverMinAug <: AbstractLinearSolver; end

"""
The function solve the linear problem associated with a linearization of the minimally augmented formulation of the Hopf bifurcation point. The keyword `debug_` is used to debug the routine by returning several key quantities.
"""
function hopfMALinearSolver(x, p::T, ω::T, pbMA::HopfProblemMinimallyAugmented, par,
	 						duu, dup, duω;
							debug_ = false) where T
	# N = length(du) - 2
	# The jacobian should be passed as a tuple as Jac_hopf_MA(u0, pb::HopfProblemMinimallyAugmented) = (return (u0, pb, d2F::Bool))
	# The Jacobian J of the vector field is expressed at (x, p)
	# the jacobian expression of the hopf problem Jhopf is
	#					[ J dpF   0
	#					 σx  σp  σω]
	########## Resolution of the bordered linear system ########
	# J * dX	  + dpF * dp		   = du => dX = x1 - dp * x2
	# The second equation
	#	<σx, dX> +  σp * dp + σω * dω = du[end-1:end]
	# thus becomes
	#   (σp - <σx, x2>) * dp + σω * dω = du[end-1:end] - <σx, x1>
	# This 2x2 system is then solved to get (dp, dω)
	############### Extraction of function names #################
	Fhandle = pbMA.F
	J = pbMA.J

	d2F = pbMA.d2F
	a = pbMA.a
	b = pbMA.b

	# parameter axis
	lens = pbMA.lens

	# we define the following jacobian. It is used at least 3 times below. This avoid doing 3 times the possibly costly building of J(x, p)
	J_at_xp = J(x, set(par, lens, p))

	# we do the following to avoid computing J_at_xp twice in case pbMA.Jadjoint is not provided
	if hasAdjoint(pbMA)
		JAd_at_xp = pbMA.Jadjoint(x, set(par, lens, p))
	else
		JAd_at_xp = transpose(J_at_xp)
	end

	δ = T(1e-9)
	ϵ1, ϵ2, ϵ3 = T(δ), T(δ), T(δ)

	# we solve Jv + a σ1 = 0 with <b, v> = n
	n = T(1)
	v, σ1, _, _ = pbMA.linbdsolver(J_at_xp, a, b, T(0), zero(x), n; shift = Complex{T}(0, ω))
	w, σ2, _, _ = pbMA.linbdsolverAdjoint(JAd_at_xp, b, a, T(0), zero(x), n; shift = -Complex{T}(0, ω))

	################### computation of σx σp ####################
	dpF   = (Fhandle(x, set(par, lens, p + ϵ1))	 - Fhandle(x, set(par, lens, p - ϵ1))) / T(2ϵ1)
	dJvdp = (apply(J(x, set(par, lens, p + ϵ3)), v) - apply(J(x, set(par, lens, p - ϵ3)), v)) / T(2ϵ3)
	σp = -dot(w, dJvdp) / n

	# case of sigma_omega
	σω = -dot(w, Complex{T}(0, 1) * v) / n

	x1, x2, _, (it1, it2) = pbMA.linsolver(J_at_xp, duu, dpF)

	# the case of ∂_xσ is a bit more involved
	# we first need to compute the value of ∂_xσ written σx
	# σx = zeros(Complex{T}, length(x))
	σx = similar(x, Complex{T})

	if hasHessian(pbMA) == false
		# We invert the jacobian of the Hopf problem when the Hessian of x -> F(x, p) is not known analytically. We thus rely on finite differences which can be slow for large dimensions
		prod(size(x)) > 1e4 && @warn "You might want to pass the Hessian, finite differences with $(prod(size(x))) unknowns"
		e = zero(x)
		for ii in CartesianIndices(x)
			e[ii] = T(1)
			d2Fve = (apply(J(x + ϵ2 * e, set(par, lens, p)), v) - apply(J(x - ϵ2 * e, set(par, lens, p)), v)) / T(2ϵ2)
			σx[ii] = -dot(w, d2Fve) / n
			e[ii] = T(0)
		end
		σxx1 = dot(σx, x1)
		σxx2 = dot(σx, x2)
	else
		d2Fv = d2F(x, p, v, x1)
		σxx1 = -dot(w, d2Fv) / n
		d2Fv = d2F(x, p, v, x2)
		σxx2 = -dot(w, d2Fv) / n
	end
	# we need to be carefull here because the dot produce conjugates. Hence the + dot(σx, x2) and + imag(dot(σx, x1) and not the opposite
	dp, dω = [real(σp - σxx2) real(σω);
			  imag(σp + σxx2) imag(σω) ] \
			  [dup - real(σxx1), duω + imag(σxx1)]

	if debug_
		return x1 - dp * x2, dp, dω, true, it1 + it2, (σx, σp, σω, dpF)
	else
		return x1 - dp * x2, dp, dω, true, it1 + it2
	end
end

function (hopfl::HopfLinearSolverMinAug)(Jhopf, du::BorderedArray{vectype, T}; debug_ = false)  where {vectype, T}
	out = hopfMALinearSolver((Jhopf.x).u,
				(Jhopf.x).p[1],
				(Jhopf.x).p[2],
				Jhopf.hopfpb,
				Jhopf.param,
				du.u, du.p[1], du.p[2];
				debug_ = debug_)
	if debug_ == false
		return BorderedArray(out[1], [out[2], out[3]]), out[4], out[5]
	else
		return BorderedArray(out[1], [out[2], out[3]]), out[4], out[5], out[6]
	end
end

################################################################################################### Newton / Continuation functions
"""
	newtonHopf(F, J, hopfpointguess::BorderedArray{vectypeR, T}, par, lens::Lens, eigenvec, eigenvec_ad, options::NewtonPar; Jt = nothing, d2F = nothing, normN = norm) where {vectypeR, T}

This function turns an initial guess for a Hopf point into a solution to the Hopf problem based on a Minimally Augmented formulation. The arguments are as follows
- `F   = (x, p) -> F(x, p)` where `p` is a set of parameters.
- `dF  = (x, p) -> d_xF(x, p)` associated jacobian
- `hopfpointguess` initial guess (x_0, p_0) for the Hopf point. It should a `BorderedArray` as returned by the function `HopfPoint`.
- `par` parameters used for the vector field
- `lens` parameter axis used to locate the Hopf point.
- `eigenvec` guess for the  iω eigenvector
- `eigenvec_ad` guess for the -iω eigenvector
- `options::NewtonPar` options for the Newton-Krylov algorithm, see [`NewtonPar`](@ref).

# Optional arguments:
- `Jt = (x, p) -> transpose(d_xF(x, p))` jacobian adjoint, it should be implemented in an efficient manner. For matrix-free methods, `transpose` is not readily available and the user must provide a dedicated method. In the case of sparse based jacobian, `Jt` should not be passed as it is computed internally more efficiently, i.e. it avoid recomputing the jacobian as it would be if you pass `Jt = (x, p) -> transpose(dF(x, p))`
- `d2F = (x, p, v1, v2) ->  d2F(x, p, v1, v2)` a bilinear operator representing the hessian of `F`. It has to provide an expression for `d2F(x,p)[v1,v2]`.
- `normN = norm`
- `kwargs` keywords arguments to be passed to the regular Newton-Krylov solver

# Simplified call:
Simplified call to refine an initial guess for a Hopf point. More precisely, the call is as follows

	newtonHopf(F, J, br::ContResult, ind_hopf::Int64, par, lens::Lens; Jt = nothing, d2F = nothing, normN = norm, options = br.contparams.newtonOptions, kwargs...)

where the optional argument `Jt` is the jacobian transpose and the Hessian is `d2F`. The parameters / options are as usual except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of bifurcation point in `br` you want to refine. You can pass newton parameters different from the ones stored in `br` by using the argument `options`.

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian will be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""
function newtonHopf(F, J, hopfpointguess::BorderedArray{vectypeR, T}, par, lens::Lens, eigenvec, eigenvec_ad, options::NewtonPar; Jt = nothing, d2F = nothing, normN = norm, kwargs...) where {vectypeR, T}
	hopfproblem = HopfProblemMinimallyAugmented(
		F, J, Jt, d2F, lens,
		_copy(eigenvec),
		_copy(eigenvec_ad),
		options.linsolver)

	# Jacobian for the Hopf problem
	Jac_hopf_MA = (x, param) -> (x = x, param = param, hopfpb = hopfproblem)

	# options for the Newton Solver
	opt_hopf = @set options.linsolver = HopfLinearSolverMinAug()

	# solve the hopf equations
	return newton(hopfproblem, Jac_hopf_MA, hopfpointguess, par, opt_hopf, normN = normN, kwargs...)
end

function newtonHopf(F, J, br::ContResult, ind_hopf::Int64, par, lens::Lens; Jt = nothing, d2F = nothing, normN = norm, options = br.contparams.newtonOptions, kwargs...)
	hopfpointguess = HopfPoint(br, ind_hopf)
	bifpt = br.bifpoint[ind_hopf]
	options.verbose && println("--> Newton Hopf, the eigenvalue considered here is ", br.eig[bifpt.idx].eigenvals[bifpt.ind_bif])
	@assert bifpt.idx == bifpt.step + 1 "Error, the bifurcation index does not refer to the correct step"
	eigenvec = geteigenvector(options.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_bif)
	eigenvec_ad = conj.(eigenvec)

	# solve the hopf equations
	return newtonHopf(F, J, hopfpointguess, par, lens, eigenvec_ad, eigenvec, options; Jt = Jt, d2F = d2F, normN = normN, kwargs...)
end

"""
$(SIGNATURES)

codim 2 continuation of Hopf points. This function turns an initial guess for a Hopf point into a curve of Hopf points based on a Minimally Augmented formulation. The arguments are as follows
- `F = (x, p) ->	F(x, p)` where `p` is a set of parameters
- `J = (x, p) -> d_xF(x, p)` associated jacobian
- `hopfpointguess` initial guess (x_0, p1_0) for the Hopf point. It should be a `Vector` or a `BorderedArray`
- `par` set of parameters
- `lens1` parameter axis for parameter 1
- `lens2` parameter axis for parameter 2
- `eigenvec` guess for the iω eigenvector at p1_0
- `eigenvec_ad` guess for the -iω eigenvector at p1_0
- `options_cont` keywords arguments to be passed to the regular [`continuation`](@ref)

# Optional arguments:

- `Jt = (x, p) -> adjoint(d_xF(x, p))` associated jacobian adjoint
- `d2F = p -> ((x, p, v1, v2) -> d2F(x, p, v1, v2))` this is the hessian of `F` computed at `(x, p)` and evaluated at `(v1, v2)`.
- `kwargs` keywords arguments to be passed to the regular [`continuation`](@ref)

# Simplified call:
The call is as follows

	continuationHopf(F, J, br::ContResult, ind_hopf::Int64, par, lens1::Lens, lens2::Lens, options_cont::ContinuationPar ;  Jt = nothing, d2F = nothing, kwargs...)

where the parameters are as above except that you have to pass the branch `br` from the result of a call to `continuation` with detection of bifurcations enabled and `index` is the index of Hopf point in `br` you want to refine.

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6

!!! tip "Jacobian tranpose"
    The adjoint of the jacobian `J` is computed internally when `Jt = nothing` by using `tranpose(J)` which works fine when `J` is an `AbstractArray`. In this case, do not pass the jacobian adjoint like `Jt = (x, p) -> transpose(d_xF(x, p))` otherwise the jacobian would be computed twice!

!!! warning "Hessian"
    The hessian of `F`, when `d2F` is not passed, is computed with Finite differences. This can be slow for many variables, e.g. ~1e6
"""
function continuationHopf(F, J, hopfpointguess::BorderedArray{vectype, Tb}, par, lens1::Lens, lens2::Lens, eigenvec, eigenvec_ad, options_cont::ContinuationPar ; Jt = nothing, d2F = p2 -> nothing, kwargs...) where {T,Tb,vectype}

	# options for the Newton Solver inheritated from the ones the user provided
	options_newton = options_cont.newtonOptions

	hopfPb = HopfProblemMinimallyAugmented(
		F, J, Jt, d2F,
		lens1,
		_copy(eigenvec), #copy(eigenvec),
		_copy(eigenvec_ad),
		options_newton.linsolver)

	# Jacobian for the Hopf problem
	Jac_hopf_MA = (x, param) -> (x = x, param = param, hopfpb = hopfPb)

	opt_hopf_cont = @set options_cont.newtonOptions.linsolver = HopfLinearSolverMinAug()

	# solve the hopf equations
	branch, u, tau = continuation(
		hopfPb, Jac_hopf_MA,
		hopfpointguess, par, lens2,
		opt_hopf_cont;
		printSolution = (u, p) -> u.p[1],
		plotSolution = (x, p; kwargs...) -> (xlabel!("p2", subplot=1); ylabel!("p1", subplot=1)), kwargs...)

	return setproperties(branch; type = :HopfCodim2, functional = hopfPb), u, tau
end

function continuationHopf(F, J, br::ContResult, ind_hopf::Int64, par, lens1::Lens, lens2::Lens, options_cont::ContinuationPar ;  Jt = nothing, d2F = nothing, kwargs...)
	hopfpointguess = HopfPoint(br, ind_hopf)
	bifpt = br.bifpoint[ind_hopf]
	eigenvec = geteigenvector(options_cont.newtonOptions.eigsolver ,br.eig[bifpt.idx].eigenvec, bifpt.ind_bif)
	eigenvec_ad = conj.(eigenvec)
	return continuationHopf(F, J, hopfpointguess, par, lens1, lens2, eigenvec_ad, eigenvec, options_cont ; Jt = Jt, d2F = d2F, kwargs...)
end


text:
using Distributions

#export AbstractConvolutionKernel, GaussianConvolutionKernel, GaussianMixtureConvolutionKernel

"""Convolution kernel using a truncated Gaussian with zero mean and unit width.
Optionally truncate at limit_kernel_width (defaults to 10)."""
struct GaussianConvolutionKernel{T} <: AbstractConvolutionKernel  where {T<:Number}
  dist::Normal{T}
  norm::T
  limit_kernel_width::Float64

  function GaussianConvolutionKernel(width::T; limit_kernel_width::Number=limit_kernel_width_default) where {T<:Number}
      @assert width > 0
      @assert limit_kernel_width > width
      dist = Normal(zero(T), width)
      norm = one(T)/(cdf(dist,limit_kernel_width)-cdf(dist,-limit_kernel_width))
      new{T}(dist, norm, limit_kernel_width)
  end
end

"Evaluate kernel at wavelength(s) x"
function (kernel::GaussianConvolutionKernel)(x::Number)
    x < kernel.limit_kernel_width^2 ? kernel.norm*pdf(kernel.dist,x) : zero(x)
end

function (kernel::GaussianConvolutionKernel)(x::A1) where { T1<:Number, A1<:AbstractArray{T1} }
    kernel.(x)
end

"""Convolution kernel using a mixture of Gaussians with zero mean and various width and weights.
Optionally truncate at limit_kernel_width (defaults to 10)."""
struct GaussianMixtureConvolutionKernel{T1,T2} <: AbstractConvolutionKernel  where {T1<:Number,T2<:Number}
    dists::Array{Normal{T1},1}
    weights::Array{T2,1}
    norm::T2
    limit_kernel_width::Float64

    function GaussianMixtureConvolutionKernel(kernel_widths::AbstractArray{T1,1}, kernel_weights::AbstractArray{T2,1};
            limit_kernel_width::Number=limit_kernel_width_default) where {T1<:Number,T2<:Number}
        @assert size(kernel_widths) == size(kernel_weights)
        @assert all(kernel_widths .> 0)
        @assert all(limit_kernel_width .> kernel_widths)
        @assert sum(kernel_weights) ≈ one(eltype(kernel_weights))
        dists = [Normal(zero(kernel_widths[i]),kernel_widths[i]) for i in 1:length(kernel_widths) ]
        cdf_at_limit(d) = cdf(d,limit_kernel_width)
        norm = sum(kernel_weights./(cdf_at_limit.(dists))) # assumes components symmetric around zero
        new{T1,T2}(dists, kernel_weights, norm, limit_kernel_width)
    end
end

"Evaluate kernel at wavelength(s) x"
function (kernel::GaussianMixtureConvolutionKernel)(x::Number)
    abs(x) < kernel.limit_kernel_width ?  kernel.norm*sum(kernel.weights.*pdf.(kernel.dists,x)) : zero(x)
end

function (kernel::GaussianMixtureConvolutionKernel)(x::AbstractArray{T}) where T<:Number
    kernel.(x)
end


text:
#
#   Tool-Module for calculating the coupling to a beam
#
#   V: 2019-04-05
#
#   Stefan Knirck
#

# This old export doesnt make a whole lot of sense. Can the func
#export matching_ratio, unnormalized_match, antenna
export unnormalized_match

function unnormalized_match(Eout, Eantenna; dx=0.01, dy=0.01)
    return sum(Eout.*conj(Eantenna)).*dx.*dy
end


text:
const Word = UInt32

struct SPIRFormatError <: Exception
    msg
end

defines_extra_operands(arg, category) = is_enum(category) && haskey(extra_operands, typeof(arg)) && haskey(extra_operands[typeof(arg)], arg)
get_extra_operands(arg) = extra_operands[typeof(arg)][arg]

function update_infos!(op_infos, i, arg, category)
    if defines_extra_operands(arg, category)
        foreach(get_extra_operands(arg)) do operand
            insert!(op_infos, i + 1, operand)
        end
    end
end

showerror(io::IO, err::SPIRFormatError) = print(io, "Invalid SPIR-V format: ", err.msg)

invalid_format(msg) = throw(SPIRFormatError(msg))

"""
SPIR-V instruction. Must contain an opcode, and optionally a type id and a result id.
"""
@broadcastref abstract type AbstractInstruction end

info(opcode::OpCode) = classes[opcode][2]
info(opcode::Integer) = info(OpCode(opcode))

function info(inst::AbstractInstruction, skip_ids::Bool = true)
    res = info(inst.opcode)
    skip_ids ? res[start_idx(inst):end] : res
end

start_idx(type_id, result_id) = 1 + !isnothing(type_id) + !isnothing(result_id)
start_idx(inst::AbstractInstruction) = start_idx(inst.type_id, inst.result_id)

function operand_kinds(opcode, skip_ids::Bool)
    op_kinds = getproperty.(info(opcode), :kind)
    if skip_ids
        filter(!in((IdResultType, IdResult)), op_kinds)
    else
        op_kinds
    end
end

operand_kinds(inst::AbstractInstruction, skip_ids::Bool = true) = getproperty.(info(inst, skip_ids), :kind)

"""
SPIR-V instruction in binary format.
"""
struct PhysicalInstruction <: AbstractInstruction
    word_count::UInt16
    opcode::UInt16
    type_id::Optional{Word}
    result_id::Optional{Word}
    operands::Vector{Word}
end

(==)(x::PhysicalInstruction, y::PhysicalInstruction) = x.word_count == y.word_count && x.opcode == y.opcode && x.type_id == y.type_id && x.result_id == y.result_id && x.operands == y.operands

"""
SPIR-V module, as a series of headers followed by a stream of instructions.
The header embeds two magic numbers, one for the module itself and one for the tool that generated it (e.g. [glslang](https://github.com/KhronosGroup/glslang)). It also contains the version of the specification applicable to the module, the maximum ID number and an optional instruction schema.
"""
struct PhysicalModule
    magic_number::Word
    generator_magic_number::Word
    version::Word
    bound::Word
    schema::Word
    instructions::Vector{PhysicalInstruction}
end

(==)(x::PhysicalModule, y::PhysicalModule) = all(prop -> getproperty(x, prop) == getproperty(y, prop), fieldnames(PhysicalModule))

function PhysicalModule(file::AbstractString)
    open(x -> PhysicalModule(x), file)
end

function PhysicalModule(io::IO)
    insts = PhysicalInstruction[]
    _magic_number = read(io, Word)
    next_word = next_word_f(_magic_number)
    version = next_word(io)
    generator_magic_number = next_word(io)
    bound = next_word(io)
    bound > 4_194_303 && invalid_format("ID bound above valid limit ($bound > $(4_194_303))")
    schema = next_word(io)
    while !eof(io)
        push!(insts, next_instruction(io, next_word))
    end

    PhysicalModule(_magic_number, generator_magic_number, version, bound, schema, insts)
end

function next_id(io::IO, next_word, op_kinds, id_type)
    @match id_type begin
        &IdResultType => begin
            if length(op_kinds) ≠ 0 && first(op_kinds) == id_type
                next_word(io)
            else
                nothing
            end
        end
        &IdResult => begin
            if length(op_kinds) == 1 && id_type == first(op_kinds) || length(op_kinds) > 1 && id_type in op_kinds[1:2]
                next_word(io)
            else
                nothing
            end
        end
    end
end

function next_instruction(io::IO, next_word)
    op_data = next_word(io)
    word_count = op_data >> 2^4

    if word_count == 0
        invalid_format("SPIR-V instructions cannot consume 0 words")
    end

    opcode = OpCode(op_data & 0xffff)
    op_kinds = operand_kinds(opcode, false)

    type_id = next_id(io, next_word, op_kinds, IdResultType)
    result_id = next_id(io, next_word, op_kinds, IdResult)

    operands = Word[]
    for i in 1:(word_count - start_idx(type_id, result_id))
        push!(operands, next_word(io))
    end

    PhysicalInstruction(word_count, opcode, type_id, result_id, operands)
end

function next_word_f(_magic_number::Word)
    if _magic_number == 0x30203270
        swap_endianness = ENDIAN_BOM == 0x04030201 ? ntoh : ltoh
        io::IO -> swap_endianness(read(io, Word))
    elseif _magic_number == magic_number
        io::IO -> read(io, Word)
    else
        error("Unknown magic number $_magic_number")
    end
end

"""
Parsed SPIR-V instruction. It represents an instruction of the form `%result_id = %opcode(%arguments...)::%type_id`.
"""
@auto_hash_equals struct Instruction <: AbstractInstruction
    opcode::OpCode
    type_id::Optional{ID}
    result_id::Optional{ID}
    arguments::Vector{Any}
    Instruction(opcode, type_id, result_id, arguments::AbstractVector) = new(convert(OpCode, opcode), convert(Optional{ID}, type_id), convert(Optional{ID}, result_id), convert(Vector{Any}, arguments))
    Instruction(opcode, type_id, result_id, arguments...) = Instruction(opcode, type_id, result_id, collect(arguments))
end

"""
Information regarding the arguments of an `Instruction`, including extra operands.
"""
function info(inst::Instruction, skip_ids::Bool = true)
    op_infos = copy(info(inst.opcode))

    # add extra operands
    foreach(enumerate(inst.arguments)) do (i, arg)
        info = op_infos[i]
        category = kind_to_category[info.kind]
        update_infos!(op_infos, i, arg, category)
    end

    skip_ids ? op_infos[start_idx(inst):end] : op_infos
end

function Instruction(inst::PhysicalInstruction)
    opcode = OpCode(inst.opcode)
    op_infos = copy(info(inst))
    op_kinds = operand_kinds(inst)
    operands = inst.operands

    arguments = []
    i = 1
    while i ≤ length(operands)
        operand = length(arguments) + 1
        info = op_infos[operand]
        category = kind_to_category[info.kind]
        if hasproperty(info, :quantifier)
            quantifier = info.quantifier
            if quantifier == "*"
                push!(arguments, operands[i:end])
                break
            elseif quantifier == "?"
                error("Unhandled '?' quantifier")
            end
        else
            j, arg = next_argument(operands[i:end], info, category)
            push!(arguments, arg)
            update_infos!(op_infos, i, arg, category)
            i += j
        end
    end
    Instruction(opcode, inst.type_id, inst.result_id, arguments)
end

function next_argument(operands, info, category)
    kind = info.kind
    if kind == LiteralString
        bytes = reinterpret(UInt8, operands)
        i, chars = parse_bytes_for_utf8_string(bytes)
        str = GC.@preserve chars unsafe_string(pointer(chars))
        div(i, 4, RoundUp), str
    elseif kind isa Literal
        arg = first(operands)
        sizeof(arg) <= 4 || error("Literals with a size greater than 32 bits are not supported.")
        1, arg
    elseif is_enum(category)
        1, kind(first(operands))
    else
        1, first(operands)
    end
end

function parse_bytes_for_utf8_string(bytes)
    utf8_chars = UInt8[]
    for (i, byte) ∈ enumerate(bytes)
        push!(utf8_chars, byte)
        byte == 0 && return (i, utf8_chars)
    end
    error("String is not NUL-terminated")
end

is_enum(category) = category in ("ValueEnum", "BitEnum")

@auto_hash_equals struct Module
    magic_number::UInt32
    generator_magic_number::UInt32
    version::VersionNumber
    bound::Int
    schema::Int
    instructions::Vector{Instruction}
end

Module(mod::PhysicalModule) = Module(mod.magic_number, mod.generator_magic_number, spirv_version(mod.version), mod.bound, mod.schema, Instruction.(mod.instructions))
Module(source) = Module(PhysicalModule(source))

function spirv_version(word)
    major = (word & 0x00ff0000) >> 16
    minor = (word & 0x0000ff00) >> 8
    VersionNumber(major, minor)
end

function spirv_version(version::VersionNumber)
    version.major << 16 + version.minor << 8
end

show(io::IO, mod::Module) = print(io, "Module(#instructions=$(length(mod.instructions)))")

Base.write(io::IO, mod::Module) = write(io, assemble(mod))

function print_diff(mod1::Module, mod2::Module)
    buff1 = IOBuffer()
    disassemble(buff1, mod1)
    seekstart(buff1)
    buff2 = IOBuffer()
    disassemble(buff2, mod2)
    seekstart(buff2)
    for (l1, l2) in zip(readlines(buff1), readlines(buff2))
        if l1 ≠ l2
            println(l1, " => ", l2)
        end
    end
end


text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "PolyPairBasis" begin

@info("-------- Test PolyPairBasis Implementation ---------")

##

using ACE1
using Printf, Test, LinearAlgebra, JuLIP, JuLIP.Testing
using JuLIP: evaluate, evaluate_d
using JuLIP.Potentials: i2z, numz

randr() = 1.0 + rand()
randcoeffs(B) = rand(length(B)) .* (1:length(B)).^(-2)

##

maxdeg = 8
r0 = 1.0
rcut = 3.0

trans = PolyTransform(1, r0)
Pr = transformed_jacobi(maxdeg, trans, rcut; pcut = 2)
pB = ACE1.PairPotentials.PolyPairBasis(Pr, :W)

##

at = bulk(:W, cubic=true) * 3
rattle!(at, 0.03)
r0 = rnn(:W)
X = copy(positions(at))

E = energy(pB, at)
println(@test length(E) == length(pB))
F = forces(pB, at)
println(@test (length(F) == length(pB)) && all(length.(F) .==  length(at)))

##

@info("test (de-)dictionisation of PairBasis")
println(@test all(JuLIP.Testing.test_fio(pB)))

##

@info("Finite-difference test on PolyPairBasis forces")
for ntest = 1:30
   E = energy(pB, at)
   DE = - forces(pB, at)
   U = [ (rand(JVecF) .- 0.5) for _=1:length(at) ]
   DExU = dot.(DE, Ref(U))
   errs = Float64[]
   for p = 2:10
      h = 0.1^p
      Eh = energy(pB, set_positions!(at, X+h*U))
      DEhxU = (Eh-E) / h
      push!(errs, norm(DExU - DEhxU, Inf))
   end
   success = (/(extrema(errs)...) < 1e-3) || (minimum(errs) < 1e-10)
   print_tf(@test success)
end
println()
##

end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 5f1fe571-3357-49cf-93d6-49ae6b4b9a05
begin
	using OrdinaryDiffEq
	using Plots, LaTeXStrings
	using PlutoUI
	# using CSV
	# using DataFrames, CategoricalArrays
	# using Turing
	# using Dates
end

# ╔═╡ 9997fcae-e8c0-11eb-29d3-19aa2548bb51
md"""
# Within-host modelling of SARS-CoV-2

*C. Contreras, J. Newby, T. Hillen.*
© 2021
"""

# ╔═╡ f02367f9-77ec-4af5-ab95-c216cc61773f

# ╔═╡ a3756580-3ac1-40b5-94ea-9e164c9065e5
md"Importing modules"

# ╔═╡ 4a72dcd8-bb16-4d63-a8a8-e89865b3ff95
md"""
Text
"""

# ╔═╡ 7541a4c5-e8bd-4700-880d-8d7e1b05cdf9

# ╔═╡ 9087b209-8385-4da8-8cfd-6eb7f4a2a0aa
begin
	md"Plot configuration"
	plotly()
	figdir = "figures/"
	# ENV["COLUMNS"]=120
	# ENV["LINES"]=20
end

# ╔═╡ 8e990cf8-3fee-4fde-b044-0e95c90a19a8
md"""
## Standard virus load function

The standard virus load function with three main phases given by

$$V(t)=v_1(t)v_2(t)v_3(t),$$

where $v _1$ describes the initial growth phase between $a _1$ and $a _2$ , $v _2$ the intermediate slow decay
phase between $a _2$ and $b _1$, and $v _3$ the final decay phase between $b _1$ and $b _2$. These functions are given as sigmoid and exponential functions, respectively

$$v_1(t) = 1 + \frac{V_{\max} - 1}{2} \left[ \tanh \left( \frac{6}{a_2 - a_1} \left( t - \frac{a_1 + a_2}{2} \right) \right) - \tanh  \left( -3\frac{a_2 + a_1}{a_2 - a_1} \right) \right]$$

$$v_2(t) = \left\{\begin{aligned}
        & 1 & t < a_2 \\
        & e^{-\alpha (t-a_2)} & t \geq a_2 \\
    \end{aligned} \right.$$

$$v_3(t) = 1 - \frac{1 - V_{\min}}{2} \left[ \tanh \left( \frac{6}{b_2 - b_1} \left( t - \frac{b_1 + b_2}{2} \right) \right) - \tanh  \left( -3\frac{b_2 + b_1}{b_2 - b_1} \right) \right],$$

.
"""

# ╔═╡ 54a73741-2b04-4bb5-9bcb-b594b451d9ec
begin
	H(x) = 0.5*(tanh(1e3*x) + 1.)
	v₁(t, a₁, a₂, Vmax) = 1. + (Vmax - 1.0)*(tanh(6.0*(t - (a₁ + a₂)/2)/(a₂ - a₁)) - tanh(-3.0*(a₂ + a₁)/(a₂ - a₁)))/2
	v₂(t, a₂, α) = 1 - H(t - a₂) + H(t - a₂)*exp(-α*(t - a₂))
	v₃(t, b₁, b₂, Vmin) = 1 - (1.0 - Vmin)*(tanh(6.0*(t - (b₁ + b₂)/2)/(b₂ - b₁)) - tanh(-3.0*(b₂ + b₁)/(b₂ - b₁)))/2
	function VLF(t, a₁, a₂, b₁, b₂, α, Vmax, Vmin)
		v₁(t, a₁, a₂, Vmax)*v₂(t, a₂, α)*v₃(t, b₁, b₂, Vmin)
	end
end

# ╔═╡ 568527c6-0760-405e-a056-98a22f95454c
md"""
## Kinetic model

The within-host model for COVID-19 infections dynamics on the lung epithilium, $E$, is given by

$$\begin{aligned}
    \frac{d E}{dt} & = \mu_E E(1-E) - \delta_E T \frac{V}{V_{max}}.
\end{aligned}$$
"""

# ╔═╡ dea4b784-87c7-4497-b11b-b1ca1cc83c71
function model!(du, u, p, t)
    E, = u
    a₁, a₂, b₁, b₂, α, Vmax, Vmin, μ_E, δ_E = p
    # Lung damage
    du[1] = dE = μ_E*E*(1.0 - E) - δ_E*E*VLF(t, a₁, a₂, b₁, b₂, α, Vmax, Vmin)/Vmax
	# Healthy lung epithelium
end

# ╔═╡ 5a3929ad-5719-4021-993c-acc31db4a674
begin
	a₁_slider = @bind a₁ Slider(0.0:0.1:4.0, default=0.4)
	α_slider = @bind α Slider(0.0:0.1:1.0, default=0.1)
	md"""Parameter

	a1: $(a₁_slider)

	alpha: $(α_slider)"""
end

# ╔═╡ b330f75d-125c-4a49-a830-5672ed18c56b
begin
	# initial
	u0 = [
		1.0   # E: Healthy lung epithilium
	]
	tspan = (0.0, 40.0)
	pars = [
		a₁     # a₁
		4.0     # a₂
		13.0    # b₁
		19.0    # b₂
		α     # α
		1e3  # Vmax
		1e-8    # Vmin
		0.01    # μ_E
		0.2     # δ_E
		]
	parnames = [:a₁, :a₂, :b₁, :b₂, :α, :Vmax, :Vmin, :μ_E, :δ_E]
	thresholds = (E = 0.3,);
	prob = ODEProblem(model!, u0, tspan, pars)
end

# ╔═╡ 2f80e9cb-9bd6-4189-9d84-20d5e1c6b0e3
function plot_sol!(sol; kwargs...)
    tt = range(0.0, tspan[2], length=100)
    pars = sol.prob.p
    vv = VLF.(sol.t, pars[1], pars[2], pars[3], pars[4], pars[5], pars[6], pars[7])
    plot!(sol.t, vv, lw=3, xaxis=("Time (days)"), yaxis=("Viral load", extrema(vv).*1.05), subplot=1; kwargs...)
    plot!(sol, vars=1, lw=3, xaxis=("Time (days)"), yaxis=("Healthy lung epithelium", (0, 1.1)), subplot=2; kwargs...)
    plot!(collect(extrema(sol.t)), [thresholds.E, thresholds.E], color=:black, lw=2, subplot=2)
end

# ╔═╡ 364206c0-954b-45b1-9f20-e5031c34de16
function plot_sol(sol; kwargs...)
    plot(layout=2, legend=:none, grid=:none, size=(600,300))
    plot_sol!(sol; kwargs...)
end

# ╔═╡ 8864fe08-3dd3-4731-ac26-7ae256835f03
begin
	sol = solve(prob, Tsit5(), dtmax=1e-1)
	plot_sol(sol)
end

# ╔═╡ 5df4f25f-ff30-46ea-b7d1-923001d900fa
md"""
$(a₁)
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
LaTeXStrings = "~1.2.1"
OrdinaryDiffEq = "~5.60.1"
Plots = "~1.19.2"
PlutoUI = "~0.7.9"
"""

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

[[Adapt]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "84918055d15b3114ede17ac6a7182f68870c16f7"
uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
version = "3.3.1"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "f87e559f87a45bece9c9ed97458d3afe98b1ebb9"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.1.0"

[[ArrayInterface]]
deps = ["IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "045ff5e1bc8c6fb1ecb28694abba0a0d55b5f4f5"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.1.17"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[Bzip2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c3598e525718abcc440f69cc6d5f60dda0a1b61e"
uuid = "6e34b625-4abd-537c-b88f-471c36dfa7a0"
version = "1.0.6+5"

[[Cairo_jll]]
deps = ["Artifacts", "Bzip2_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "e2f47f6d8337369411569fd45ae5753ca10394c6"
uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
version = "1.16.0+6"

[[ChainRulesCore]]
deps = ["Compat", "LinearAlgebra", "SparseArrays"]
git-tree-sha1 = "f53ca8d41e4753c41cdafa6ec5f7ce914b34be54"
uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
version = "0.10.13"

[[ColorSchemes]]
deps = ["ColorTypes", "Colors", "FixedPointNumbers", "Random", "StaticArrays"]
git-tree-sha1 = "ed268efe58512df8c7e224d2e170afd76dd6a417"
uuid = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
version = "3.13.0"

[[ColorTypes]]
deps = ["FixedPointNumbers", "Random"]
git-tree-sha1 = "024fe24d83e4a5bf5fc80501a314ce0d1aa35597"
uuid = "3da002f7-5984-5a60-b8a6-cbb66c0b333f"
version = "0.11.0"

[[Colors]]
deps = ["ColorTypes", "FixedPointNumbers", "Reexport"]
git-tree-sha1 = "417b0ed7b8b838aa6ca0a87aadf1bb9eb111ce40"
uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
version = "0.12.8"

[[CommonSolve]]
git-tree-sha1 = "68a0743f578349ada8bc911a5cbd5a2ef6ed6d1f"
uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
version = "0.2.0"

[[CommonSubexpressions]]
deps = ["MacroTools", "Test"]
git-tree-sha1 = "7b8a93dba8af7e3b42fecabf646260105ac373f7"
uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
version = "0.3.0"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "dc7dedc2c2aa9faf59a55c622760a25cbefbe941"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.31.0"

[[CompilerSupportLibraries_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"

[[ConstructionBase]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "f74e9d5388b8620b4cee35d4c5a618dd4dc547f4"
uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
version = "1.3.0"

[[Contour]]
deps = ["StaticArrays"]
git-tree-sha1 = "9f02045d934dc030edad45944ea80dbd1f0ebea7"
uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
version = "0.5.7"

[[DataAPI]]
git-tree-sha1 = "ee400abb2298bd13bfc3df1c412ed228061a2385"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.7.0"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "4437b64df1e0adccc3e5d1adbc3ac741095e4677"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.9"

[[DataValueInterfaces]]
git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
uuid = "e2d170a0-9d28-54be-80f0-106bbe20a464"
version = "1.0.0"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[DiffEqBase]]
deps = ["ArrayInterface", "ChainRulesCore", "DataStructures", "DocStringExtensions", "FastBroadcast", "FunctionWrappers", "IterativeSolvers", "LabelledArrays", "LinearAlgebra", "Logging", "MuladdMacro", "NonlinearSolve", "Parameters", "Printf", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "Requires", "SciMLBase", "Setfield", "SparseArrays", "StaticArrays", "Statistics", "SuiteSparse", "ZygoteRules"]
git-tree-sha1 = "932153f62d0508e59733e0fc33361470d293a889"
uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
version = "6.68.1"

[[DiffResults]]
deps = ["StaticArrays"]
git-tree-sha1 = "c18e98cba888c6c25d1c3b048e4b3380ca956805"
uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
version = "1.0.3"

[[DiffRules]]
deps = ["NaNMath", "Random", "SpecialFunctions"]
git-tree-sha1 = "214c3fcac57755cfda163d91c58893a8723f93e9"
uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
version = "1.0.2"

[[Distances]]
deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
git-tree-sha1 = "abe4ad222b26af3337262b8afb28fab8d215e9f8"
uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
version = "0.10.3"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[DocStringExtensions]]
deps = ["LibGit2"]
git-tree-sha1 = "a32185f5428d3986f47c2ab78b1f216d5e6cc96f"
uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
version = "0.8.5"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EarCut_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "92d8f9f208637e8d2d28c664051a00569c01493d"
uuid = "5ae413db-bbd1-5e63-b57d-d24a61df00f5"
version = "2.1.5+1"

[[Expat_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b3bfd02e98aedfa5cf885665493c5598c350cd2f"
uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
version = "2.2.10+0"

[[ExponentialUtilities]]
deps = ["ArrayInterface", "LinearAlgebra", "Printf", "Requires", "SparseArrays"]
git-tree-sha1 = "ad435656c49da7615152b856c0f9abe75b0b5dc9"
uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
version = "1.8.4"

[[FFMPEG]]
deps = ["FFMPEG_jll"]
git-tree-sha1 = "b57e3acbe22f8484b4b5ff66a7499717fe1a9cc8"
uuid = "c87230d0-a227-11e9-1b43-d7ebe4e7570a"
version = "0.4.1"

[[FFMPEG_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "LibVPX_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "Pkg", "Zlib_jll", "libass_jll", "libfdk_aac_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
git-tree-sha1 = "3cc57ad0a213808473eafef4845a74766242e05f"
uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
version = "4.3.1+4"

[[FastBroadcast]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "26be48918640ce002f5833e8fc537b2ba7ed0234"
uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
version = "0.1.8"

[[FastClosures]]
git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
uuid = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
version = "0.3.2"

[[FiniteDiff]]
deps = ["ArrayInterface", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays"]
git-tree-sha1 = "f6f80c8f934efd49a286bb5315360be66956dfc4"
uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
version = "2.8.0"

[[FixedPointNumbers]]
deps = ["Statistics"]
git-tree-sha1 = "335bfdceacc84c5cdf16aadc768aa5ddfc5383cc"
uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
version = "0.8.4"

[[Fontconfig_jll]]
deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "35895cf184ceaab11fd778b4590144034a167a2f"
uuid = "a3f928ae-7b40-5064-980b-68af3947d34b"
version = "2.13.1+14"

[[Formatting]]
deps = ["Printf"]
git-tree-sha1 = "8339d61043228fdd3eb658d86c926cb282ae72a8"
uuid = "59287772-0a20-5a39-b81b-1366585eb4c0"
version = "0.4.2"

[[ForwardDiff]]
deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "NaNMath", "Printf", "Random", "SpecialFunctions", "StaticArrays"]
git-tree-sha1 = "e2af66012e08966366a43251e1fd421522908be6"
uuid = "f6369f11-7733-5829-9624-2563aa707210"
version = "0.10.18"

[[FreeType2_jll]]
deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "cbd58c9deb1d304f5a245a0b7eb841a2560cfec6"
uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
version = "2.10.1+5"

[[FriBidi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "aa31987c2ba8704e23c6c8ba8a4f769d5d7e4f91"
uuid = "559328eb-81f9-559d-9380-de523a88c83c"
version = "1.0.10+0"

[[FunctionWrappers]]
git-tree-sha1 = "241552bc2209f0fa068b6415b1942cc0aa486bcc"
uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
version = "1.1.2"

[[Future]]
deps = ["Random"]
uuid = "9fa8497b-333b-5362-9e8d-4d0656e87820"

[[GLFW_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libglvnd_jll", "Pkg", "Xorg_libXcursor_jll", "Xorg_libXi_jll", "Xorg_libXinerama_jll", "Xorg_libXrandr_jll"]
git-tree-sha1 = "dba1e8614e98949abfa60480b13653813d8f0157"
uuid = "0656b61e-2033-5cc2-a64a-77c0f6c09b89"
version = "3.3.5+0"

[[GR]]
deps = ["Base64", "DelimitedFiles", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Pkg", "Printf", "Random", "Serialization", "Sockets", "Test", "UUIDs"]
git-tree-sha1 = "9f473cdf6e2eb360c576f9822e7c765dd9d26dbc"
uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
version = "0.58.0"

[[GR_jll]]
deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Pkg", "Qt5Base_jll", "Zlib_jll", "libpng_jll"]
git-tree-sha1 = "eaf96e05a880f3db5ded5a5a8a7817ecba3c7392"
uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
version = "0.58.0+0"

[[GeometryBasics]]
deps = ["EarCut_jll", "IterTools", "LinearAlgebra", "StaticArrays", "StructArrays", "Tables"]
git-tree-sha1 = "15ff9a14b9e1218958d3530cc288cf31465d9ae2"
uuid = "5c1252a2-5f33-56bf-86c9-59e7332b4326"
version = "0.3.13"

[[Gettext_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "9b02998aba7bf074d14de89f9d37ca24a1a0b046"
uuid = "78b55507-aeef-58d4-861c-77aaff3498b1"
version = "0.21.0+0"

[[Glib_jll]]
deps = ["Artifacts", "Gettext_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "47ce50b742921377301e15005c96e979574e130b"
uuid = "7746bdde-850d-59dc-9ae8-88ece973131d"
version = "2.68.1+0"

[[Grisu]]
git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
version = "1.0.2"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "c6a1fff2fd4b1da29d3dccaffb1e1001244d844e"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.12"

[[Hwloc]]
deps = ["Hwloc_jll"]
git-tree-sha1 = "92d99146066c5c6888d5a3abc871e6a214388b91"
uuid = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
version = "2.0.0"

[[Hwloc_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "3395d4d4aeb3c9d31f5929d32760d8baeee88aaf"
uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
version = "2.5.0+0"

[[IfElse]]
git-tree-sha1 = "28e837ff3e7a6c3cdb252ce49fb412c8eb3caeef"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.0"

[[Inflate]]
git-tree-sha1 = "f5fc07d4e706b84f72d54eedcc1c13d92fb0871c"
uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
version = "0.1.2"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[IterTools]]
git-tree-sha1 = "05110a2ab1fc5f932622ffea2a003221f4782c18"
uuid = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
version = "1.3.0"

[[IterativeSolvers]]
deps = ["LinearAlgebra", "Printf", "Random", "RecipesBase", "SparseArrays"]
git-tree-sha1 = "1a8c6237e78b714e901e406c096fc8a65528af7d"
uuid = "42fd0dbc-a981-5370-80f2-aaf504508153"
version = "0.9.1"

[[IteratorInterfaceExtensions]]
git-tree-sha1 = "a3f24677c21f5bbe9d2a714f95dcd58337fb2856"
uuid = "82899510-4779-5014-852e-03e436cf321d"
version = "1.0.0"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "81690084b6198a2e1da36fcfda16eeca9f9f24e4"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.1"

[[JpegTurbo_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d735490ac75c5cb9f1b00d8b5509c11984dc6943"
uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
version = "2.1.0+0"

[[LAME_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "f6250b16881adf048549549fba48b1161acdac8c"
uuid = "c1c5ebd0-6772-5130-a774-d5fcae4a789d"
version = "3.100.1+0"

[[LZO_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "e5b909bcf985c5e2605737d2ce278ed791b89be6"
uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
version = "2.10.1+0"

[[LaTeXStrings]]
git-tree-sha1 = "c7f1c695e06c01b95a67f0cd1d34994f3e7db104"
uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
version = "1.2.1"

[[LabelledArrays]]
deps = ["ArrayInterface", "LinearAlgebra", "MacroTools", "StaticArrays"]
git-tree-sha1 = "5e38cfdd771c34821ade5515f782fe00865d60b3"
uuid = "2ee39098-c373-598a-b85f-a56591580800"
version = "1.6.2"

[[Latexify]]
deps = ["Formatting", "InteractiveUtils", "LaTeXStrings", "MacroTools", "Markdown", "Printf", "Requires"]
git-tree-sha1 = "a4b12a1bd2ebade87891ab7e36fdbce582301a92"
uuid = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
version = "0.15.6"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[LibVPX_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "12ee7e23fa4d18361e7c2cde8f8337d4c3101bc7"
uuid = "dd192d2f-8180-539f-9fb4-cc70b1dcf69a"
version = "1.10.0+0"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libffi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "761a393aeccd6aa92ec3515e428c26bf99575b3b"
uuid = "e9f186c6-92d2-5b65-8a66-fee21dc1b490"
version = "3.2.2+0"

[[Libgcrypt_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll", "Pkg"]
git-tree-sha1 = "64613c82a59c120435c067c2b809fc61cf5166ae"
uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
version = "1.8.7+0"

[[Libglvnd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
git-tree-sha1 = "7739f837d6447403596a75d19ed01fd08d6f56bf"
uuid = "7e76a0d4-f3c7-5321-8279-8d96eeed0f29"
version = "1.3.0+3"

[[Libgpg_error_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "c333716e46366857753e273ce6a69ee0945a6db9"
uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
version = "1.42.0+0"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[Libmount_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "9c30530bf0effd46e15e0fdcf2b8636e78cbbd73"
uuid = "4b2f31a3-9ecc-558c-b454-b3730dcb73e9"
version = "2.35.0+0"

[[Libtiff_jll]]
deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "Zlib_jll", "Zstd_jll"]
git-tree-sha1 = "340e257aada13f95f98ee352d316c3bed37c8ab9"
uuid = "89763e89-9b03-5906-acba-b20f662cd828"
version = "4.3.0+0"

[[Libuuid_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7f3efec06033682db852f8b3bc3c1d2b0a0ab066"
uuid = "38a345b3-de98-5d2b-a5d3-14cd9215e700"
version = "2.36.0+0"

[[LightGraphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "432428df5f360964040ed60418dd5601ecd240b6"
uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
version = "1.3.5"

[[LineSearches]]
deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Parameters", "Printf"]
git-tree-sha1 = "f27132e551e959b3667d8c93eae90973225032dd"
uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
version = "7.1.1"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[LogExpFunctions]]
deps = ["DocStringExtensions", "LinearAlgebra"]
git-tree-sha1 = "7bd5f6565d80b6bf753738d2bc40a5dfea072070"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.2.5"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoopVectorization]]
deps = ["ArrayInterface", "DocStringExtensions", "IfElse", "LinearAlgebra", "OffsetArrays", "Polyester", "Requires", "SLEEFPirates", "Static", "StrideArraysCore", "ThreadingUtilities", "UnPack", "VectorizationBase"]
git-tree-sha1 = "e03414d325902b1690d02a528ef46f9b29daaa73"
uuid = "bdcacae8-1622-11e9-2a5c-532679323890"
version = "0.12.51"

[[MacroTools]]
deps = ["Markdown", "Random"]
git-tree-sha1 = "6a8a2a625ab0dea913aba95c11370589e0239ff0"
uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
version = "0.5.6"

[[ManualMemory]]
git-tree-sha1 = "71c64ebe61a12bad0911f8fc4f91df8a448c604c"
uuid = "d125e4d3-2237-4719-b19c-fa641b8a4667"
version = "0.1.4"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Measures]]
git-tree-sha1 = "e498ddeee6f9fdb4551ce855a46f54dbd900245f"
uuid = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
version = "0.3.1"

[[Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "4ea90bd5d3985ae1f9a908bd4500ae88921c5ce7"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.0"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[MuladdMacro]]
git-tree-sha1 = "c6190f9a7fc5d9d5915ab29f2134421b12d24a68"
uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
version = "0.2.2"

[[NLSolversBase]]
deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
git-tree-sha1 = "50608f411a1e178e0129eab4110bd56efd08816f"
uuid = "d41bc354-129a-5804-8e4c-c37616107c6c"
version = "7.8.0"

[[NLsolve]]
deps = ["Distances", "LineSearches", "LinearAlgebra", "NLSolversBase", "Printf", "Reexport"]
git-tree-sha1 = "019f12e9a1a7880459d0173c182e6a99365d7ac1"
uuid = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
version = "4.5.1"

[[NaNMath]]
git-tree-sha1 = "bfe47e760d60b82b66b61d2d44128b62e3a369fb"
uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
version = "0.3.5"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[NonlinearSolve]]
deps = ["ArrayInterface", "FiniteDiff", "ForwardDiff", "IterativeSolvers", "LinearAlgebra", "RecursiveArrayTools", "RecursiveFactorization", "Reexport", "SciMLBase", "Setfield", "StaticArrays", "UnPack"]
git-tree-sha1 = "ef18e47df4f3917af35be5e5d7f5d97e8a83b0ec"
uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
version = "0.3.8"

[[OffsetArrays]]
deps = ["Adapt"]
git-tree-sha1 = "4f825c6da64aebaa22cc058ecfceed1ab9af1c7e"
uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
version = "1.10.3"

[[Ogg_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7937eda4681660b4d6aeeecc2f7e1c81c8ee4e2f"
uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
version = "1.3.5+0"

[[OpenSSL_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "15003dcb7d8db3c6c857fda14891a539a8f2705a"
uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
version = "1.1.10+0"

[[OpenSpecFun_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "13652491f6856acfd2db29360e1bbcd4565d04f1"
uuid = "efe28fd5-8261-553b-a9e1-b2916fc3738e"
version = "0.5.5+0"

[[Opus_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "51a08fb14ec28da2ec7a927c4337e4332c2a4720"
uuid = "91d4177d-7536-5919-b921-800302f37372"
version = "1.3.2+0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[OrdinaryDiffEq]]
deps = ["Adapt", "ArrayInterface", "DataStructures", "DiffEqBase", "DocStringExtensions", "ExponentialUtilities", "FastClosures", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "NLsolve", "Polyester", "RecursiveArrayTools", "Reexport", "SparseArrays", "SparseDiffTools", "StaticArrays", "UnPack"]
git-tree-sha1 = "379ed10814cba5087978bce9df31161a496620d3"
uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
version = "5.60.1"

[[PCRE_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b2a7af664e098055a7529ad1a900ded962bca488"
uuid = "2f80f16e-611a-54ab-bc61-aa92de5b98fc"
version = "8.44.0+0"

[[Parameters]]
deps = ["OrderedCollections", "UnPack"]
git-tree-sha1 = "2276ac65f1e236e0a6ea70baff3f62ad4c625345"
uuid = "d96e819e-fc66-5662-9728-84c9c7592b0a"
version = "0.12.2"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "c8abc88faa3f7a3950832ac5d6e690881590d6dc"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "1.1.0"

[[Pixman_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "b4f5d02549a10e20780a24fce72bea96b6329e29"
uuid = "30392449-352a-5448-841d-b1acce4e97dc"
version = "0.40.1+0"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[PlotThemes]]
deps = ["PlotUtils", "Requires", "Statistics"]
git-tree-sha1 = "a3a964ce9dc7898193536002a6dd892b1b5a6f1d"
uuid = "ccf2f8ad-2431-5c83-bf29-c5338b663b6a"
version = "2.0.1"

[[PlotUtils]]
deps = ["ColorSchemes", "Colors", "Dates", "Printf", "Random", "Reexport", "Statistics"]
git-tree-sha1 = "501c20a63a34ac1d015d5304da0e645f42d91c9f"
uuid = "995b91a9-d308-5afd-9ec6-746e21dbc043"
version = "1.0.11"

[[Plots]]
deps = ["Base64", "Contour", "Dates", "FFMPEG", "FixedPointNumbers", "GR", "GeometryBasics", "JSON", "Latexify", "LinearAlgebra", "Measures", "NaNMath", "PlotThemes", "PlotUtils", "Printf", "REPL", "Random", "RecipesBase", "RecipesPipeline", "Reexport", "Requires", "Scratch", "Showoff", "SparseArrays", "Statistics", "StatsBase", "UUIDs"]
git-tree-sha1 = "f3d4d35b8cb87adc844c05c722f505776ac29988"
uuid = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
version = "1.19.2"

[[PlutoUI]]
deps = ["Base64", "Dates", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "Suppressor"]
git-tree-sha1 = "44e225d5837e2a2345e69a1d1e01ac2443ff9fcb"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.9"

[[Polyester]]
deps = ["ArrayInterface", "IfElse", "ManualMemory", "Requires", "Static", "StrideArraysCore", "ThreadingUtilities", "VectorizationBase"]
git-tree-sha1 = "f5a74523ebc205723baefb95874f708741199d70"
uuid = "f517fe37-dbe3-4b94-8317-1923a5111588"
version = "0.3.4"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[Qt5Base_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Pkg", "Xorg_libXext_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "xkbcommon_jll"]
git-tree-sha1 = "ad368663a5e20dbb8d6dc2fddeefe4dae0781ae8"
uuid = "ea2cea3b-5b76-57ae-a6ef-0a8af62496e1"
version = "5.15.3+0"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[RecipesBase]]
git-tree-sha1 = "b3fb709f3c97bfc6e948be68beeecb55a0b340ae"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.1.1"

[[RecipesPipeline]]
deps = ["Dates", "NaNMath", "PlotUtils", "RecipesBase"]
git-tree-sha1 = "2a7a2469ed5d94a98dea0e85c46fa653d76be0cd"
uuid = "01d81517-befc-4cb6-b9ec-a95719d0359c"
version = "0.3.4"

[[RecursiveArrayTools]]
deps = ["ArrayInterface", "ChainRulesCore", "DocStringExtensions", "LinearAlgebra", "RecipesBase", "Requires", "StaticArrays", "Statistics", "ZygoteRules"]
git-tree-sha1 = "2a76e8f24c67f3ebecaccefa8d4abd27db828407"
uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
version = "2.14.9"

[[RecursiveFactorization]]
deps = ["LinearAlgebra", "LoopVectorization"]
git-tree-sha1 = "2e1a88c083ebe8ba69bc0b0084d4b4ba4aa35ae0"
uuid = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
version = "0.1.13"

[[Reexport]]
git-tree-sha1 = "5f6c21241f0f655da3952fd60aa18477cf96c220"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.1.0"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[SLEEFPirates]]
deps = ["IfElse", "Static", "VectorizationBase"]
git-tree-sha1 = "da6d214ffc85b1292f300649ef86d3c4f9aaf25d"
uuid = "476501e8-09a2-5ece-8869-fb82de89a1fa"
version = "0.6.22"

[[SciMLBase]]
deps = ["ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "RecipesBase", "RecursiveArrayTools", "StaticArrays", "Statistics", "Tables", "TreeViews"]
git-tree-sha1 = "932aaae93e81686e6473f27b5a11c403576a2183"
uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
version = "1.18.0"

[[Scratch]]
deps = ["Dates"]
git-tree-sha1 = "0b4b7f1393cff97c33891da2a0bf69c6ed241fda"
uuid = "6c6a2e73-6563-6170-7368-637461726353"
version = "1.1.0"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[Setfield]]
deps = ["ConstructionBase", "Future", "MacroTools", "Requires"]
git-tree-sha1 = "d5640fc570fb1b6c54512f0bd3853866bd298b3e"
uuid = "efcf1570-3423-57d1-acb7-fd33fddbac46"
version = "0.7.0"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Showoff]]
deps = ["Dates", "Grisu"]
git-tree-sha1 = "91eddf657aca81df9ae6ceb20b959ae5653ad1de"
uuid = "992d4aef-0814-514b-bc4d-f2e9a6c4116f"
version = "1.0.3"

[[SimpleTraits]]
deps = ["InteractiveUtils", "MacroTools"]
git-tree-sha1 = "5d7e3f4e11935503d3ecaf7186eac40602e7d231"
uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
version = "0.9.4"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SortingAlgorithms]]
deps = ["DataStructures"]
git-tree-sha1 = "b3363d7460f7d098ca0912c69b082f75625d7508"
uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
version = "1.0.1"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[SparseDiffTools]]
deps = ["Adapt", "ArrayInterface", "Compat", "DataStructures", "FiniteDiff", "ForwardDiff", "LightGraphs", "LinearAlgebra", "Requires", "SparseArrays", "StaticArrays", "VertexSafeGraphs"]
git-tree-sha1 = "cc8703d09d1e7337fb794dbf6c9634afaa4f49c7"
uuid = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
version = "1.14.1"

[[SpecialFunctions]]
deps = ["ChainRulesCore", "LogExpFunctions", "OpenSpecFun_jll"]
git-tree-sha1 = "a50550fa3164a8c46747e62063b4d774ac1bcf49"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.5.1"

[[Static]]
deps = ["IfElse"]
git-tree-sha1 = "2740ea27b66a41f9d213561a04573da5d3823d4b"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.2.5"

[[StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "1b9a0f17ee0adde9e538227de093467348992397"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.7"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StatsAPI]]
git-tree-sha1 = "1958272568dc176a1d881acb797beb909c785510"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.0.0"

[[StatsBase]]
deps = ["DataAPI", "DataStructures", "LinearAlgebra", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
git-tree-sha1 = "2f6792d523d7448bbe2fec99eca9218f06cc746d"
uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
version = "0.33.8"

[[StrideArraysCore]]
deps = ["ArrayInterface", "ManualMemory", "Requires", "ThreadingUtilities", "VectorizationBase"]
git-tree-sha1 = "2525850b6887e217c3cc149ab29eb2a7a6360d37"
uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
version = "0.1.16"

[[StructArrays]]
deps = ["Adapt", "DataAPI", "StaticArrays", "Tables"]
git-tree-sha1 = "000e168f5cc9aded17b6999a560b7c11dda69095"
uuid = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
version = "0.6.0"

[[SuiteSparse]]
deps = ["Libdl", "LinearAlgebra", "Serialization", "SparseArrays"]
uuid = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"

[[Suppressor]]
git-tree-sha1 = "a819d77f31f83e5792a76081eee1ea6342ab8787"
uuid = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
version = "0.2.0"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[TableTraits]]
deps = ["IteratorInterfaceExtensions"]
git-tree-sha1 = "c06b2f539df1c6efa794486abfb6ed2022561a39"
uuid = "3783bdb8-4a98-5b6b-af9a-565f29a5fe9c"
version = "1.0.1"

[[Tables]]
deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "LinearAlgebra", "TableTraits", "Test"]
git-tree-sha1 = "8ed4a3ea724dac32670b062be3ef1c1de6773ae8"
uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
version = "1.4.4"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[ThreadingUtilities]]
deps = ["ManualMemory"]
git-tree-sha1 = "593009385e6536dc50baedc2c8bc1a3717ef358f"
uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
version = "0.4.5"

[[TreeViews]]
deps = ["Test"]
git-tree-sha1 = "8d0d7a3fe2f30d6a7f833a5f19f7c7a5b396eae6"
uuid = "a2a6695c-b41b-5b7d-aed9-dbfdeacea5d7"
version = "0.3.0"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[UnPack]]
git-tree-sha1 = "387c1f73762231e86e0c9c5443ce3b4a0a9a0c2b"
uuid = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
version = "1.0.2"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[VectorizationBase]]
deps = ["ArrayInterface", "Hwloc", "IfElse", "Libdl", "LinearAlgebra", "Static"]
git-tree-sha1 = "b609848b41158b8cc13a72dc6deb9511f9a96187"
uuid = "3d5dd08c-fd9d-11e8-17fa-ed2836048c2f"
version = "0.20.22"

[[VertexSafeGraphs]]
deps = ["LightGraphs"]
git-tree-sha1 = "b9b450c99a3ca1cc1c6836f560d8d887bcbe356e"
uuid = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
version = "0.1.2"

[[Wayland_jll]]
deps = ["Artifacts", "Expat_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Pkg", "XML2_jll"]
git-tree-sha1 = "3e61f0b86f90dacb0bc0e73a0c5a83f6a8636e23"
uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
version = "1.19.0+0"

[[Wayland_protocols_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll"]
git-tree-sha1 = "2839f1c1296940218e35df0bbb220f2a79686670"
uuid = "2381bf8a-dfd0-557d-9999-79630e7b1b91"
version = "1.18.0+4"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[XSLT_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgcrypt_jll", "Libgpg_error_jll", "Libiconv_jll", "Pkg", "XML2_jll", "Zlib_jll"]
git-tree-sha1 = "91844873c4085240b95e795f692c4cec4d805f8a"
uuid = "aed1982a-8fda-507f-9586-7b0439959a61"
version = "1.1.34+0"

[[Xorg_libX11_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
git-tree-sha1 = "5be649d550f3f4b95308bf0183b82e2582876527"
uuid = "4f6342f7-b3d2-589e-9d20-edeb45f2b2bc"
version = "1.6.9+4"

[[Xorg_libXau_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4e490d5c960c314f33885790ed410ff3a94ce67e"
uuid = "0c0b7dd1-d40b-584c-a123-a41640f87eec"
version = "1.0.9+4"

[[Xorg_libXcursor_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXfixes_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "12e0eb3bc634fa2080c1c37fccf56f7c22989afd"
uuid = "935fb764-8cf2-53bf-bb30-45bb1f8bf724"
version = "1.2.0+4"

[[Xorg_libXdmcp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "4fe47bd2247248125c428978740e18a681372dd4"
uuid = "a3789734-cfe1-5b06-b2d0-1dd0d9d62d05"
version = "1.1.3+4"

[[Xorg_libXext_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "b7c0aa8c376b31e4852b360222848637f481f8c3"
uuid = "1082639a-0dae-5f34-9b06-72781eeb8cb3"
version = "1.3.4+4"

[[Xorg_libXfixes_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "0e0dc7431e7a0587559f9294aeec269471c991a4"
uuid = "d091e8ba-531a-589c-9de9-94069b037ed8"
version = "5.0.3+4"

[[Xorg_libXi_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXfixes_jll"]
git-tree-sha1 = "89b52bc2160aadc84d707093930ef0bffa641246"
uuid = "a51aa0fd-4e3c-5386-b890-e753decda492"
version = "1.7.10+4"

[[Xorg_libXinerama_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll"]
git-tree-sha1 = "26be8b1c342929259317d8b9f7b53bf2bb73b123"
uuid = "d1454406-59df-5ea1-beac-c340f2130bc3"
version = "1.1.4+4"

[[Xorg_libXrandr_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libXext_jll", "Xorg_libXrender_jll"]
git-tree-sha1 = "34cea83cb726fb58f325887bf0612c6b3fb17631"
uuid = "ec84b674-ba8e-5d96-8ba1-2a689ba10484"
version = "1.5.2+4"

[[Xorg_libXrender_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "19560f30fd49f4d4efbe7002a1037f8c43d43b96"
uuid = "ea2f1a96-1ddc-540d-b46f-429655e07cfa"
version = "0.9.10+4"

[[Xorg_libpthread_stubs_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "6783737e45d3c59a4a4c4091f5f88cdcf0908cbb"
uuid = "14d82f49-176c-5ed1-bb49-ad3f5cbd8c74"
version = "0.1.0+3"

[[Xorg_libxcb_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "XSLT_jll", "Xorg_libXau_jll", "Xorg_libXdmcp_jll", "Xorg_libpthread_stubs_jll"]
git-tree-sha1 = "daf17f441228e7a3833846cd048892861cff16d6"
uuid = "c7cfdc94-dc32-55de-ac96-5a1b8d977c5b"
version = "1.13.0+3"

[[Xorg_libxkbfile_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll"]
git-tree-sha1 = "926af861744212db0eb001d9e40b5d16292080b2"
uuid = "cc61e674-0454-545c-8b26-ed2c68acab7a"
version = "1.1.0+4"

[[Xorg_xcb_util_image_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "0fab0a40349ba1cba2c1da699243396ff8e94b97"
uuid = "12413925-8142-5f55-bb0e-6d7ca50bb09b"
version = "0.4.0+1"

[[Xorg_xcb_util_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxcb_jll"]
git-tree-sha1 = "e7fd7b2881fa2eaa72717420894d3938177862d1"
uuid = "2def613f-5ad1-5310-b15b-b15d46f528f5"
version = "0.4.0+1"

[[Xorg_xcb_util_keysyms_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "d1151e2c45a544f32441a567d1690e701ec89b00"
uuid = "975044d2-76e6-5fbe-bf08-97ce7c6574c7"
version = "0.4.0+1"

[[Xorg_xcb_util_renderutil_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "dfd7a8f38d4613b6a575253b3174dd991ca6183e"
uuid = "0d47668e-0667-5a69-a72c-f761630bfb7e"
version = "0.3.9+1"

[[Xorg_xcb_util_wm_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xcb_util_jll"]
git-tree-sha1 = "e78d10aab01a4a154142c5006ed44fd9e8e31b67"
uuid = "c22f9ab0-d5fe-5066-847c-f4bb1cd4e361"
version = "0.4.1+1"

[[Xorg_xkbcomp_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libxkbfile_jll"]
git-tree-sha1 = "4bcbf660f6c2e714f87e960a171b119d06ee163b"
uuid = "35661453-b289-5fab-8a00-3d9160c6a3a4"
version = "1.4.2+4"

[[Xorg_xkeyboard_config_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_xkbcomp_jll"]
git-tree-sha1 = "5c8424f8a67c3f2209646d4425f3d415fee5931d"
uuid = "33bec58e-1273-512f-9401-5d533626f822"
version = "2.27.0+4"

[[Xorg_xtrans_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "79c31e7844f6ecf779705fbc12146eb190b7d845"
uuid = "c5fb5394-a638-5e4d-96e5-b29de1b5cf10"
version = "1.4.0+3"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[Zstd_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "cc4bf3fdde8b7e3e9fa0351bdeedba1cf3b7f6e6"
uuid = "3161d3a3-bdf6-5164-811a-617609db77b4"
version = "1.5.0+0"

[[ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "9e7a1e8ca60b742e508a315c17eef5211e7fbfd7"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.1"

[[libass_jll]]
deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "acc685bcf777b2202a904cdcb49ad34c2fa1880c"
uuid = "0ac62f75-1d6f-5e53-bd7c-93b484bb37c0"
version = "0.14.0+4"

[[libfdk_aac_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "7a5780a0d9c6864184b3a2eeeb833a0c871f00ab"
uuid = "f638f0a6-7fb0-5443-88ba-1cc74229b280"
version = "0.1.6+4"

[[libpng_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Zlib_jll"]
git-tree-sha1 = "94d180a6d2b5e55e447e2d27a29ed04fe79eb30c"
uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
version = "1.6.38+0"

[[libvorbis_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Ogg_jll", "Pkg"]
git-tree-sha1 = "c45f4e40e7aafe9d086379e5578947ec8b95a8fb"
uuid = "f27f6e37-5d2b-51aa-960f-b287f2bc3b7a"
version = "1.3.7+0"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"

[[x264_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "d713c1ce4deac133e3334ee12f4adff07f81778f"
uuid = "1270edf5-f2f9-52d2-97e9-ab00b5d0237a"
version = "2020.7.14+2"

[[x265_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "487da2f8f2f0c8ee0e83f39d13037d6bbf0a45ab"
uuid = "dfaa095f-4041-5dcd-9319-2fabd8486b76"
version = "3.0.0+3"

[[xkbcommon_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Wayland_jll", "Wayland_protocols_jll", "Xorg_libxcb_jll", "Xorg_xkeyboard_config_jll"]
git-tree-sha1 = "ece2350174195bb31de1a63bea3a41ae1aa593b6"
uuid = "d8fb68d0-12a3-5cfd-a85a-d49703b185fd"
version = "0.9.1+5"
"""

# ╔═╡ Cell order:
# ╟─9997fcae-e8c0-11eb-29d3-19aa2548bb51
# ╠═f02367f9-77ec-4af5-ab95-c216cc61773f
# ╟─a3756580-3ac1-40b5-94ea-9e164c9065e5
# ╠═5f1fe571-3357-49cf-93d6-49ae6b4b9a05
# ╠═4a72dcd8-bb16-4d63-a8a8-e89865b3ff95
# ╠═7541a4c5-e8bd-4700-880d-8d7e1b05cdf9
# ╠═9087b209-8385-4da8-8cfd-6eb7f4a2a0aa
# ╠═8e990cf8-3fee-4fde-b044-0e95c90a19a8
# ╠═54a73741-2b04-4bb5-9bcb-b594b451d9ec
# ╠═568527c6-0760-405e-a056-98a22f95454c
# ╠═dea4b784-87c7-4497-b11b-b1ca1cc83c71
# ╠═b330f75d-125c-4a49-a830-5672ed18c56b
# ╠═364206c0-954b-45b1-9f20-e5031c34de16
# ╠═2f80e9cb-9bd6-4189-9d84-20d5e1c6b0e3
# ╠═8864fe08-3dd3-4731-ac26-7ae256835f03
# ╠═5a3929ad-5719-4021-993c-acc31db4a674
# ╠═5df4f25f-ff30-46ea-b7d1-923001d900fa
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct XLSXInput <: File
        name::String
        sheet::String
        range::String
        descriptor::String
    end

Read .xlsx file.

# Arguments
- `name::String`: input file name
- `sheet::String`: input sheet name
- `range::String`: input sheet range
- `descriptor::String`: file descriptor
"""
mutable struct XLSXInput <: File
    "input file name"
    name::String
    "input sheet name"
    sheet::String
    "input sheet range"
    range::String
    "file descriptor"
    descriptor::String
end

function XLSXInput(; name, sheet, range, descriptor, )
    XLSXInput(name, sheet, range, descriptor, )
end

"""Get [`XLSXInput`](@ref) `name`."""
get_name(value::XLSXInput) = value.name
"""Get [`XLSXInput`](@ref) `sheet`."""
get_sheet(value::XLSXInput) = value.sheet
"""Get [`XLSXInput`](@ref) `range`."""
get_range(value::XLSXInput) = value.range
"""Get [`XLSXInput`](@ref) `descriptor`."""
get_descriptor(value::XLSXInput) = value.descriptor

"""Set [`XLSXInput`](@ref) `name`."""
set_name!(value::XLSXInput, val) = value.name = val
"""Set [`XLSXInput`](@ref) `sheet`."""
set_sheet!(value::XLSXInput, val) = value.sheet = val
"""Set [`XLSXInput`](@ref) `range`."""
set_range!(value::XLSXInput, val) = value.range = val
"""Set [`XLSXInput`](@ref) `descriptor`."""
set_descriptor!(value::XLSXInput, val) = value.descriptor = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
using DataFrames, DataFramesMeta, Statistics
using StatsPlots, Plots
using PyPlot
using Plots.PlotMeasures
using CSV

# Functions to process dataframe files (outputs of simulations)
# for heatmap plots (Figure 3 and Figure 6)
# for jitter plots (Figure 4 and Figure 7)

# Top-level call
# Get mean error across seeds and write it in a csv file
# For heatmap plot of Exact Bayes (Figure 3 and Figure 6)
function runmakeheatmap_csv(df, savepath; prefix = "gaussian_")
    df2 = df_mergeseeds_meanerrors(df)
    learners = unique(df2[:learner])
    for i in learners
        dfheat = DataFrame()
        z, rows, cols, colstring = makeheatmapmatrix_csv(df2, i)
        @show i
        dfheat[Symbol(colstring)] = rows
        dfheat[Symbol("changeprobability")] = cols
        dfheat[Symbol(i*"_meanerror")] = z
        CSV.write(joinpath(savepath, prefix*"heatmap_"* i *".csv"), dfheat,
            delim = " ");
    end
end
# Top-level call
# 1. Get mean error of Exact Bayes across seeds
# 2. Get mean error of all other learners across seeds
# 3. Compute their difference and write it in a csv file
# For heatmap plots (Figure 3 and Figure 6)
function runmakeheatmap_diff_csv(df, savepath; learnerbase = "bayesfilter",
                                prefix = "gaussian_")
    df2 = df_mergeseeds_meanerrors(df)
    zbase, rowsbase, colsbase, colstringbase = makeheatmapmatrix_csv(df2, learnerbase)
    learners = unique(df2[:learner])
    for i in learners
        dfheatdiff = DataFrame()
        z, rows, cols, colstring = makeheatmapmatrix_csv(df2, i)
        # Calc difference from benchmarck (Exact Bayes)
        z_diff = z .- zbase
        @show i
        @show z_diff
        dfheatdiff[Symbol(colstring)] = rows
        dfheatdiff[Symbol("changeprobability")] = cols
        dfheatdiff[Symbol(i*"_meanerror")] = z_diff
        CSV.write(joinpath(savepath, prefix*"heatmap_diff_"* i *".csv"), dfheatdiff,
            delim = " ");
    end
end
# Same as above, but calculate also the std of the difference, considering them
# as not independent (because the environments are same).
# For jitter plots (Figure 4 and Figure 7)
function runmakeheatmap_diff_stats_nonindependent_csv(df, savepath; learnerbase = "bayesfilter",
                                prefix = "gaussian_")
    df_meanerrors_acrossseeds = df_meanerrors_acrosseeds(df) # Calculate the mean error within every seed and collect them all across seeds

    zbase_mean_seed, rowsbase_mean_seed, colsbase_mean_seed, colstringbase_mean_seed = makeheatmapmatrix_std_csv(df_meanerrors_acrossseeds,
                                                                                learnerbase, quantity=:allmeanerrors_seed)
    learners = unique(df_meanerrors_acrossseeds[:learner])
    for i in learners
        dfheatdiff = DataFrame()
        z_mean_seed, rows_mean_seed, cols_mean_seed, colstring_mean_seed = makeheatmapmatrix_std_csv(df_meanerrors_acrossseeds, i, quantity=:allmeanerrors_seed)
        # Calc difference from benchmarck (Exact Bayes)
        diff_mean_seeds = z_mean_seed .- zbase_mean_seed
        z_diff_mean = mean.(diff_mean_seeds)
        z_diff_std = std.(diff_mean_seeds)
        z_diff_se = z_diff_std ./ sqrt.(length.(diff_mean_seeds))
        dfheatdiff[Symbol(colstring_mean_seed)] = rows_mean_seed
        dfheatdiff[Symbol("changeprobability")] = cols_mean_seed
        dfheatdiff[Symbol(i*"_meanerror")] = z_diff_mean
        dfheatdiff[Symbol(i*"_diff_std")] = z_diff_std
        dfheatdiff[Symbol(i*"_diff_se")] = z_diff_se
        CSV.write(joinpath(savepath, prefix*"heatmap_diff_std_"* i *".csv"), dfheatdiff,
            delim = " ");
    end
end
# Put information as we want it in order to create the heatmaps in pgfplots
function makeheatmapmatrix_csv(df, learnername)
    df = @where(df, :learner .== learnername)
    z = []; rows = []; cols = [];
    if in(:sigma, names(df))
        colstring = "sigma"
        rows = df[:sigma]; cols = df[:changeprobability]
        z = df[:meanerrors]
    elseif in(:stochasticity, names(df))
        colstring = "stochasticity"
        rows = df[:stochasticity]; cols = df[:changeprobability]
        z = df[:meanerrors]
    end
    z, rows, cols, colstring
end
# Put information as we want it in order to create the heatmaps in pgfplots
function makeheatmapmatrix_std_csv(df, learnername; quantity=:stderrors_seed)
    df = @where(df, :learner .== learnername)
    z = []; rows = []; cols = [];
    if in(:sigma, names(df))
        colstring = "sigma"
        rows = df[:sigma]; cols = df[:changeprobability]
        z = df[Symbol(quantity)]
    elseif in(:stochasticity, names(df))
        colstring = "stochasticity"
        rows = df[:stochasticity]; cols = df[:changeprobability]
        z = df[Symbol(quantity)]
    end
    z, rows, cols, colstring
end
function makeindicesvec(valuesvec)
    indicesvec = zeros(length(valuesvec))
    [indicesvec[i] = j for i in 1:length(valuesvec) for j in 1:length(unique(valuesvec)) if valuesvec[i]==unique(valuesvec)[j]]
end
function df_addcolstringlabel(df)
    df = @transform(df, learnerlabel = :learner * string(formatparamlist(:param)))
end

# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# ------------------------------------------------------------------------------
# --------- Plot heatmaps within julia: not used eventually for the paper ------
# 1. Get mean error of all other learners across seeds and plot it in a matrix
function runmakeheatmap(df, savepath;
                        ylabelstring = "Stochasticity level",
                        clims = [0., 0.2],
                        prefix = "gaussian_")
    df2 = df_mergeseeds_meanerrors(df)
    learners = unique(df2[:learner])
    for i in learners
        makeheatmap(df2, i, savepath, ylabelstring = ylabelstring, clims = clims,
                    prefix = prefix)
    end
end
# 1. Get mean error of Exact Bayes across seeds
# 2. Get mean error of all other learners across seeds
# 3. Compute their difference and plot it in a matrix
function runmakeheatmap_diff(df, savepath;
                        learnerbase = "bayesfilter",
                        ylabelstring = "Stochasticity level",
                        clims = [-0.005, 0.3],
                        prefix = "gaussian_")
    df2 = df_mergeseeds_meanerrors(df)
    zmatbase, rowsbase, rows_intbase, colsbase, cols_intbase = makeheatmapmatrix(df2, learnerbase)
    learners = unique(df2[:learner])
    for i in learners
        zmat, rows, rows_int, cols, cols_int = makeheatmapmatrix(df2, i)
        # Calc diff
        zmat_diff = zmat .- zmatbase
        pyplot()
        heatmap(zmat_diff, y_flip = true,
            xticks = (unique(cols_int), unique(cols)),
            yticks = (unique(rows_int), unique(rows)),
            clims = (clims[1], clims[2]),
            framestyle = :grid,
            xlabel = "Pc",
            ylabel = ylabelstring,
            title = "Diff Mean errors. Learner: " * i,
            titlefont = font(10),
            guidefont = font(10),
            tickfont = font(10),
            bottom_margin = 5mm, top_margin = 5mm, left_margin = 10mm)
        Plots.savefig(joinpath(savepath, prefix*"heatmap_diff_"* i *".png"))
    end
end
function makeheatmap(df, learnername, savepath;
                        ylabelstring = "Stochasticity level",
                        clims = [0., 0.2],
                        prefix = "gaussian_")
    zmat, rows, rows_int, cols, cols_int = makeheatmapmatrix(df, learnername)
    pyplot()
    heatmap(zmat, y_flip = true,
        xticks = (unique(cols_int), unique(cols)),
        yticks = (unique(rows_int), unique(rows)),
        clims = (clims[1], clims[2]),
        framestyle = :grid,
        xlabel = "Pc",
        ylabel = ylabelstring,
        title = "Mean errors. Learner: " * learnername,
        titlefont = font(10),
        guidefont = font(10),
        tickfont = font(10),
        bottom_margin = 5mm, top_margin = 5mm, left_margin = 10mm)
    Plots.savefig(joinpath(savepath, prefix*"heatmap_"* learnername *".png"))
end
function makeheatmapmatrix(df, learnername)
    df = @where(df, :learner .== learnername)
    zmat = []; rows = []; rows_int = []; cols = []; cols_int = []
    if in(:sigma, names(df))
        rows = df[:sigma]; cols = df[:changeprobability]
        rows_int = makeindicesvec(rows); cols_int = makeindicesvec(cols)
        zmat = zeros(length(unique(rows)), length(unique(cols)))
        for k in 1:length(rows)
            zmat[rows_int[k], cols_int[k]] = @where(df, :sigma .== rows[k], :changeprobability .== cols[k])[:meanerrors][1]
        end
    elseif in(:stochasticity, names(df))
        rows = df[:stochasticity]; cols = df[:changeprobability]
        rows_int = makeindicesvec(rows); cols_int = makeindicesvec(cols)
        zmat = zeros(length(unique(rows)), length(unique(cols)))
        for k in 1:length(rows)
            zmat[rows_int[k], cols_int[k]] = @where(df, :stochasticity .== rows[k], :changeprobability .== cols[k])[:meanerrors][1]
        end
    end
    zmat, rows, rows_int, cols, cols_int
end


text:
using MIDI
using Test

include("variablelength.jl")
include("sysexevent.jl")
include("note.jl")
include("midievent.jl")
include("metaevent.jl")
include("miditrack.jl")
include("midiio.jl")
include("negative_delta.jl")


text:
using RICoverAlg

# Queueing problems
m = 2;  n = 40 ; lmd = 1.0
g = ones(n); g /=sum(g)
U = zeros(m, n)
for i in (1:n)
   U[2,i] = 28 - i
end
Z = exp.(U/lmd)

risolve(U, g, lmd)

t1 = @elapsed begin
waitProbs = []  ## will store the unconditional probability of waiting for varying parameter
for parameter in (0:0.01:50)
    for i in (1:n)
       U[2,i] = parameter - i
    end
    Z = exp.(U/lmd)

result = risolve(U, g, lmd)
push!(waitProbs,result[1][2])   ## push!(list,x) adds x to list - here the second entry (lining up probs) of the result
end
end

@test t1 <= 0.43

# using BenchmarkTools
# @btime begin
# waitProbs = []  ## will store the unconditional probability of waiting for varying parameter
# for parameter in (0:0.01:50)
#     for i in (1:n)
#        U[2,i] = parameter - i
#     end
#     Z = exp.(U/lmd)
#
# result = risolve(U, g, lmd)
# push!(waitProbs,result[1][2])   ## push!(list,x) adds x to list - here the second entry (lining up probs) of the result
# end
# end
### output: 96.626 ms (1165614 allocations: 174.65 MiB)


text:
for ff in (exp,)
    f = nameof(ff)
    m = Base.parentmodule(ff)
    for PT in (:Particles, :StaticParticles)
        @eval begin
            @inline function ($m.$f)(p::$PT{Float64,N}) where N
                res = (SLEEFPirates.$f).(p.particles)
                return $PT{Float64,N}(res)
            end
            @inline function ($m.$f)(p::$PT{Float32,N}) where N
                res = (SLEEFPirates.$f).(p.particles)
                return $PT{Float32,N}(res)
            end
        end
    end
end

for ff in (log, sin, cos, asin, acos, atan) # tan is not faster
    f = nameof(ff)
    fs = Symbol(f,"_fast")
    m = Base.parentmodule(ff)
    for PT in (:Particles, :StaticParticles)
        @eval begin
            @inline function ($m.$f)(p::$PT{Float64,N}) where N
                res = (SLEEFPirates.$fs).(p.particles)
                return $PT{Float64,N}(res)
            end
            @inline function ($m.$f)(p::$PT{Float32,N}) where N
                res = (SLEEFPirates.$fs).(p.particles)
                return $PT{Float32,N}(res)
            end
        end
    end
end