text:
#Using the convention that generator satisfies ϕ(0) = 1 (this is opposite to e.g. https://en.wikipedia.org/wiki/Copula_(probability_theory))

# We should follow closely
# https://www.uni-ulm.de/fileadmin/website_uni_ulm/mawi.inst.zawa/forschung/preprintmariushofert.pdf
# for the implementation of sampling methods.

abstract type ArchimedeanCopula{d} <: Copula{d} end
function Distributions.cdf(C::CT,u) where {CT<:ArchimedeanCopula}
    @assert length(C) == length(u)
    sum_ϕ⁻¹u = 0.0
    for us in u
        sum_ϕ⁻¹u += ϕ⁻¹(C,us)
    end
    return ϕ(C,sum_ϕ⁻¹u)
end
ϕ⁽¹⁾(C::CT, t) where {CT<:ArchimedeanCopula} = ForwardDiff.derivative(x -> ϕ(C,x), t)
function ϕ⁽ᵈ⁾(C::ArchimedeanCopula{d},t) where d
    X = Taylor1(eltype(t),d)
    taylor_expansion = ϕ(C,t+X)
    coef = getcoeff(taylor_expansion,d) # gets the dth coef.
    der = coef * factorial(d) # gets the dth derivative of $\phi$ taken in t.
    return der
end
function Distributions._logpdf(C::CT, u) where {CT<:ArchimedeanCopula}
    d = length(C)
    # @assert d == length(u) "Dimension mismatch"
    sum_ϕ⁻¹u = 0.0
    sum_logϕ⁽¹⁾ϕ⁻¹u = 0.0
    for us in u
        ϕ⁻¹u = ϕ⁻¹(C,us)
        sum_ϕ⁻¹u += ϕ⁻¹u
        sum_logϕ⁽¹⁾ϕ⁻¹u += log(-ϕ⁽¹⁾(C,ϕ⁻¹u)) #log of negative here because ϕ⁽¹⁾ is necessarily negative
    end

    numer = ϕ⁽ᵈ⁾(C, sum_ϕ⁻¹u)
    dimension_sign = iseven(d) ? 1.0 : -1.0 #need this for log since (-1.0)ᵈ ϕ⁽ᵈ⁾ ≥ 0.0
    return log(dimension_sign*numer) - sum_logϕ⁽¹⁾ϕ⁻¹u
end

ϕ(C::ArchimedeanCopula{d},x) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
ϕ⁻¹(C::ArchimedeanCopula{d},x) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
radial_dist(C::ArchimedeanCopula{d}) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
τ(C::ArchimedeanCopula{d}) where d  = @error "Archimedean interface not implemented for $(typeof(C)) yet."
τ⁻¹(::ArchimedeanCopula{d},τ) where d = @error "Archimedean interface not implemented for $(typeof(C)) yet."
# radial_dist(C::ArchimedeanCopula) = laplace_transform(t -> ϕ(C,t))

function Distributions._rand!(rng::Distributions.AbstractRNG, C::CT, x::AbstractVector{T}) where {T<:Real, CT<:ArchimedeanCopula}
    r = rand(rng,radial_dist(C))
    Random.rand!(rng,x)
    for i in 1:length(C)
        x[i] = ϕ(C,-log(x[i])/r)
    end
    return x
end
function Base.rand(rng::Distributions.AbstractRNG,C::CT) where CT<: ArchimedeanCopula
    x = rand(rng,length(C))
    r = rand(rng,radial_dist(C))
    for i in 1:length(C)
        x[i] = ϕ(C,-log(x[i])/r)
    end
    return x
end
function Distributions.fit(::Type{CT},u) where {CT <: ArchimedeanCopula}
    # @info "Archimedean fits are by default through inverse kendall tau."
    d = size(u,1)
    τ = StatsBase.corkendall(u')
    # Then the off-diagonal elements of the matrix should be averaged:
    avgτ = (sum(τ) .- d) / (d^2-d)
    θ = τ⁻¹(CT,avgτ)
    return CT(d,θ)
end


text:
export pca

using mlpack._Internal.io

import mlpack_jll
const pcaLibrary = mlpack_jll.libmlpack_julia_pca

# Call the C binding of the mlpack pca binding.
function pca_mlpackMain()
  success = ccall((:pca, pcaLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module pca_internal
  import ..pcaLibrary

end # module

"""
    pca(input; [decomposition_method, new_dimensionality, scale, var_to_retain, verbose])

This program performs principal components analysis on the given dataset using
the exact, randomized, randomized block Krylov, or QUIC SVD method. It will
transform the data onto its principal components, optionally performing
dimensionality reduction by ignoring the principal components with the smallest
eigenvalues.

Use the `input` parameter to specify the dataset to perform PCA on.  A desired
new dimensionality can be specified with the `new_dimensionality` parameter, or
the desired variance to retain can be specified with the `var_to_retain`
parameter.  If desired, the dataset can be scaled before running PCA with the
`scale` parameter.

Multiple different decomposition techniques can be used.  The method to use can
be specified with the `decomposition_method` parameter, and it may take the
values 'exact', 'randomized', or 'quic'.

For example, to reduce the dimensionality of the matrix `data` to 5 dimensions
using randomized SVD for the decomposition, storing the output matrix to
`data_mod`, the following command can be used:

```julia
julia> using CSV
julia> data = CSV.read("data.csv")
julia> data_mod = pca(data; decomposition_method="randomized",
            new_dimensionality=5)
```

# Arguments

 - `input::Array{Float64, 2}`: Input dataset to perform PCA on.
 - `decomposition_method::String`: Method used for the principal
      components analysis: 'exact', 'randomized', 'randomized-block-krylov',
      'quic'.  Default value `exact`.

 - `new_dimensionality::Int`: Desired dimensionality of output dataset. If
      0, no dimensionality reduction is performed.  Default value `0`.

 - `scale::Bool`: If set, the data will be scaled before running PCA, such
      that the variance of each feature is 1.  Default value `false`.

 - `var_to_retain::Float64`: Amount of variance to retain; should be
      between 0 and 1.  If 1, all variance is retained.  Overrides -d.  Default
      value `0`.

 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Float64, 2}`: Matrix to save modified dataset to.

"""
function pca(input;
             decomposition_method::Union{String, Missing} = missing,
             new_dimensionality::Union{Int, Missing} = missing,
             scale::Union{Bool, Missing} = missing,
             var_to_retain::Union{Float64, Missing} = missing,
             verbose::Union{Bool, Missing} = missing,
             points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, pcaLibrary), Nothing, ());

  IORestoreSettings("Principal Components Analysis")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  if !ismissing(decomposition_method)
    IOSetParam("decomposition_method", convert(String, decomposition_method))
  end
  if !ismissing(new_dimensionality)
    IOSetParam("new_dimensionality", convert(Int, new_dimensionality))
  end
  if !ismissing(scale)
    IOSetParam("scale", convert(Bool, scale))
  end
  if !ismissing(var_to_retain)
    IOSetParam("var_to_retain", convert(Float64, var_to_retain))
  end
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  pca_mlpackMain()

  return IOGetParamMat("output", points_are_rows)
end

ext: jl


text:
using SparseArrays, DiffEqOperators, LinearAlgebra, Random,
      Test, BandedMatrices, FillArrays

function second_derivative_stencil(N)
  A = zeros(N,N+2)
  for i in 1:N, j in 1:N+2
      (j-i==0 || j-i==2) && (A[i,j]=1)
      j-i==1 && (A[i,j]=-2)
  end
  A
end

# Analytic solutions to higher order operators.
# Do not modify unless you are completely certain of the changes.
function fourth_deriv_approx_stencil(N, scale)
    A = zeros(N,N+2)
    A[1,1:8] = [3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3]
    A[2,1:8] = [2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6]
    A[N-1,N-5:end] = reverse([2/3 -11/6 0.0 31/6 -22/3 4.5 -4/3 1/6], dims=2)
    A[N,N-5:end] = reverse([3.5 -56/3 42.5 -54.0 251/6 -20.0 5.5 -2/3], dims=2)
    for i in 3:N-2
        A[i,i-2:i+4] = [-1/6 2.0 -13/2 28/3 -13/2 2.0 -1/6]
    end
    return A*scale
end

function second_deriv_fourth_approx_stencil(N, scale)
    A = zeros(N,N+2)
    A[1,1:6] = [5/6 -1.25 -1/3 7/6 -0.5 5/60]
    A[N,N-3:end] = reverse([5/6 -1.25 -1/3 7/6 -0.5 5/60], dims=2)
    for i in 2:N-1
        A[i,i-1:i+3] = [-1/12 4/3 -5/2 4/3 -1/12]
    end
    return A*scale
end

# Tests the corrrectness of stencils.
# Do not modify the following test-set unless you are completely certain of your changes.
@testset "Correctness of Stencils on Uniform grid" begin
    N = 20
    dx = ones(Float64, N+1)

    L = CenteredDifference(2,2, dx, N)
    correct = convert(Array{Float64, 1}, [1,-2,1])

    for stencil in L.stencil.coefs
        @test stencil ≈ correct
    end

    L = CenteredDifference(4,4, dx, N)
    correct = fourth_deriv_approx_stencil(N, 1.0)

    # Check that stencils (according to convert_by_multiplication) agree with correct
    @test convert_by_multiplication(Array, L, N) ≈ correct

    # Check that concretization agrees correct
    @test Array(L) ≈ correct
    @test sparse(L) ≈ correct
    @test BandedMatrix(L) ≈ correct

    scale = 0.1
    L = CenteredDifference(2,4, dx/scale, N)
    correct = second_deriv_fourth_approx_stencil(N, scale)

    # Check that stencils (according to convert_by_multiplication) agree with correct
    @test convert_by_multiplication(Array, L, N) ≈ correct

    # Check that concretization agrees correct
    @test Array(L) ≈ correct
    @test sparse(L) ≈ correct
    @test BandedMatrix(L) ≈ correct
end

# Cross check from http://web.media.mit.edu/~crtaylor/calculator.html
@testset "Correctness of Stencils on Non-Uniform grid" begin
    N = 4
    dx = [0.1, 0.2, 0.3, 0.4, 0.5]

    L = CenteredDifference(2,2, dx, N)
    correct_interior = [[200/3,  -300/2,  100/3],
                        [60/3,   -100/3,  40/3],
                        [200/21, -350/21, 150/21],
                        [50/9,   -90/9,   40/9]]

    for i in length(L.stencil_coefs)
        @test L.stencil_coefs[i] ≈ correct_interior[i]
    end
end

# tests for full and sparse function
@testset "Full and Sparse functions:" begin
    N = 10
    d_order = 2
    approx_order = 2
    correct = second_derivative_stencil(N)
    A = CenteredDifference(d_order,approx_order,1.0,N)

    @test convert_by_multiplication(Array,A,N) == correct
    @test Array(A) == second_derivative_stencil(N)
    @test sparse(A) == second_derivative_stencil(N)
    @test BandedMatrix(A) == second_derivative_stencil(N)
    @test opnorm(A, Inf) == opnorm(correct, Inf)

    # testing higher derivative and approximation concretization
    N = 20
    d_order = 4
    approx_order = 4
    A = CenteredDifference(d_order,approx_order,1.0,N)
    correct = convert_by_multiplication(Array,A,N)

    @test Array(A) ≈ correct
    @test sparse(A) ≈ correct
    @test BandedMatrix(A) ≈ correct

    N = 26
    d_order = 8
    approx_order = 8
    A = CenteredDifference(d_order,approx_order,1.0,N)
    correct = convert_by_multiplication(Array,A,N)

    @test Array(A) ≈ correct
    @test sparse(A) ≈ correct
    @test BandedMatrix(A) ≈ correct

    # testing correctness of multiplication
    N = 1000
    d_order = 4
    approx_order = 10
    y = collect(1:1.0:N+2).^4 - 2*collect(1:1.0:N+2).^3 + collect(1:1.0:N+2).^2;
    y = convert(Array{BigFloat, 1}, y)

    A = CenteredDifference(d_order,approx_order,one(BigFloat),N)
    correct = convert_by_multiplication(Array,A,N)
    @test Array(A) ≈ correct
    @test sparse(A) ≈ correct
    @test BandedMatrix(A) ≈ correct
    @test A*y ≈ Array(A)*y
end

@testset "Indexing tests" begin
    N = 1000
    d_order = 4
    approx_order = 10

    A = CenteredDifference(d_order,approx_order,1.0,N)
    @test A[1,1] == Array(A)[1,1]
    @test A[10,20] == 0

    correct = Array(A)
    for i in 1:N
        @test A[i,i] == correct[i,i]
    end

    # Indexing Tests
    N = 1000
    d_order = 2
    approx_order = 2

    A = CenteredDifference(d_order,approx_order,1.0,N)
    M = Array(A,1000)
    @test A[1,1] == M[1,1]
    @test A[1:4,1] == M[1:4,1]
    @test A[5,2:10] == M[5,2:10]
    @test A[60:100,500:600] == M[60:100,500:600]
end

@testset begin "Operations on matrices"
    N = 51
    M = 101
    d_order = 2
    approx_order = 2

    xarr = range(0,stop=1,length=N)
    yarr = range(0,stop=1,length=M)
    dx = xarr[2]-xarr[1]
    dy = yarr[2]-yarr[1]
    F = [x^2+y for x = xarr, y = yarr]

    A = CenteredDifference(d_order,approx_order,dx,length(xarr)-2)
    B = CenteredDifference(d_order,approx_order,dy,length(yarr))

    @test A*F ≈ 2*ones(N-2,M) atol=1e-2
    F*B
    A*F*B

    G = [x^2+y^2 for x = xarr, y = yarr]

    @test A*G ≈ 2*ones(N-2,M) atol=1e-2
    G*B
    A*G*B
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

using Test, Profile

function busywait(t, n_tries)
    iter = 0
    while iter < n_tries && Profile.len_data() == 0
        iter += 1
        tend = time() + t
        while time() < tend end
    end
end

Profile.clear()
@profile busywait(1, 20)

let r = Profile.retrieve()
    mktemp() do path, io
        serialize(io, r)
        close(io)
        open(path) do io
            @test isa(deserialize(io), Tuple{Vector{UInt},Dict{UInt64,Vector{StackFrame}}})
        end
    end
end

let iobuf = IOBuffer()
    Profile.print(iobuf, format=:tree, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:tree, maxdepth=2)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, C=true)
    str = String(take!(iobuf))
    @test !isempty(str)
    truncate(iobuf, 0)
    Profile.print(iobuf)
    @test !isempty(String(take!(iobuf)))
    truncate(iobuf, 0)
    Profile.print(iobuf, format=:flat, sortedby=:count)
    @test !isempty(String(take!(iobuf)))
    Profile.clear()
    @test isempty(Profile.fetch())
end

# issue #13229
module I13229
using Test, Profile
global z = 0
@timed @profile for i = 1:5
    function f(x)
        return x + i
    end
    global z = f(i)
end
@test z == 10
end

max_stars_repo_name: YingboMa/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using Test
using NeXLCore

@testset "Compton" begin
    cxr = n"U L3-M5"
    @test isapprox(comptonShift(π/4, cxr), 1.0/(1.0+(energy(cxr)/0.511e6)*(1.0-cos(π/4))), rtol=1.0e-5)
    @test isapprox(energy(cxr)*(1.0 - comptonShift(π, cxr)), 689.0, atol=1.0)
    @test isapprox(comptonEnergy(π/3, cxr), energy(cxr)/(1.0+(energy(cxr)/0.511e6)*(1.0-cos(π/3))), rtol=1.0e-5)
    @test isapprox(energy(cxr) - comptonEnergy(π, cxr), 689.0, atol=1.0)
    @test energy(cxr) - comptonEnergy(π/2, cxr) < 689.0
    # These test values are taken from Wikipedia
    @test isapprox(comptonAngular(π/2, 2.75) / comptonAngular(0.0, 2.75), 0.5, atol=0.0001)
    @test isapprox(comptonAngular(π/2, 60.0e3) / comptonAngular(0.0, 60.0e3), 0.4, atol=0.02)
    @test isapprox(comptonAngular(π/2, 511.0e3) / comptonAngular(0.0, 511.0e3), 0.2, atol=0.02)
end


text:
eom!(dv, v, u, p, t) = total_acceleration!(dv, p)
system_callbacks(s) = CallbackSet(boundary_conditions!(s), thermostats!(s), neighbours(s))

function DiffEqBase.SecondOrderODEProblem{true}(system::AbstractNBodySystem)
    du0, u0, tspan = initial_conditions(system)
    callbacks = system_callbacks(system)

    SecondOrderODEProblem{true}(eom!, du0, u0, tspan, callbacks)
end


text:
using Test

using Astro

print("leap years")
@test is_leap_year(2004) == true
@test is_leap_year(333, false) == false
println(" passed")

print("7.a Convert Gregorian date to Julian day number")
jd = cal_to_jd(1957, 10, 4.81)
@test jd == 2436116.31
println(" passed")

print("7.a Convert Gregorian date to Julian day number — meeus'")
jd = cal_to_jd_m(1957, 10, 4.81)
@test jd == 2436116.31
println(" passed")

print("7.b Convert Julian date to Julian day number")
jd = cal_to_jd(333, 1, 27.5, false)
# currently fails, is off by one
#    @test isapprox(jd, 1842713.0, atol=0.1)
println(" passed")

print("7.b Convert Julian date to Julian day number — meeus")
jd = cal_to_jd_m(333, 1, 27.5, false)
# currently fails, is off by one
#    @test isapprox(jd, 1842713.0, atol=0.1)
println(" passed")

print("7.c Convert Julian day number to Gregorian date")
yr, mo, day = jd_to_cal(2436116.31)
@test yr == 1957
@test mo == 10
@test isapprox(day, 4.81, atol=0.1)
println(" passed")

print("7.c(1) Convert Julian day number to Julian date")
yr, mo, day = jd_to_cal(1842713.0, false)
@test yr == 333
@test mo == 1
@test isapprox(day, 27.5, atol=0.01)
println(" passed")

print("7.c(2) Convert Julian day number to Julian date")
yr, mo, day = jd_to_cal(1507900.13, false)
@test yr == -584
@test mo == 5
@test isapprox(day, 28.63, atol=0.01)
println(" passed")

print("7.d Time interval in days")
jd0 = cal_to_jd(1910, 4, 20.0)
jd1 = cal_to_jd(1986, 2, 9.0)
@test jd1 - jd0 == 27689
println(" passed")

print("7.d(1) Time interval in days")
jd = cal_to_jd(1991, 7, 11)
jd = jd + 10000
yr, mo, day = jd_to_date(jd)
@test yr == 2018
@test mo == 11
@test day == 26
dow = jd_to_day_of_week(jd)
@test dow == 1
println(" passed")

print("7.f Day of the year")
N = cal_to_day_of_year(1978, 11, 14)
@test N == 318
println(" passed")

print("7.g Day of the year")
N = cal_to_day_of_year(1988, 4, 22)
@test N == 113
println(" passed")

print("7(pg 66-1) Day of the year to calendar")
mo, day = day_of_year_to_cal(1978, 318)
@test mo == 11
@test day == 14
println(" passed")

print("7(pg 66-2) Day of the year to calendar")
mo, day = day_of_year_to_cal(1988, 113)
@test mo == 4
@test day == 22
println(" passed")

tbl = Array[
    [1991, 3, 31],
    [1992, 4, 19],
    [1993, 4, 11],
    [1954, 4, 18],
    [2000, 4, 23],
    [1818, 3, 22]
]

print("8(pg 68) Gregorian Easter (6 times)")
for (yr, mo, day) in tbl
    # calculate xmo and xday given yr
    xmo, xday = easter(yr)
    map(print, [" ", yr, " ", mo, " ", day, " -> ", yr, " ", xmo, " ", xday])
    println()
    @test xmo == mo
    @test xday == day
end
println(" passed")

print("8(pg 69) Julian Easter (3 times)")
for yr in [179, 711, 1243]
    mo, day = easter(yr, false)
    @test mo == 4
    @test day == 12
end
println(" passed")

print("9.a Jewish Pesach")
mo, day = pesach(1990)
@test mo == 4
@test day == 10
println(" passed")

print("9.b Julian date of Moslem year 1421 first day")
yr, mo, day = moslem_to_christian(1421, 1, 1)
@test yr == 2000
@test mo == 4
@test day == 6
println(" passed")

print("9.c Moslem date for Aug 13 1991 gregorian")
yr, mo, day = christian_to_moslem(1991, 8, 13)
@test yr == 1412
@test mo == 2
@test day == 2
println(" passed")

print("Julian day to date")
(y, month, d, h, m, s) = jd_to_date(2436116.31, true)
@test y == 1957
@test month == 10
@test d == 4
@test h == 19
@test m == 26
@test isapprox(s, 24, atol=0.1)
println(" passed")

# 2014 September 5 09:41:24  2456905.90375 Friday
print("is September 5 DST?")
@test is_dst(2456905.90375) == true
println(" passed")

print("number of Julian centuries since J2000.0")
@test jd_to_jcent(2456905.90375) == 0.14677354551677085
println(" passed")

print("Julian date to formatted string 1")
@test lt_to_str(2456905.90375, "GMT", "second") == "Fri, 05 Sep 2014 09:41:24"
println(" passed")

print("Julian date to formatted string 2")
@test lt_to_str(2456905.90375, "GMT", "day") == "Fri, 05 Sep 2014"
println(" passed")


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Structured optimizer interface
# ========================== #
"""
    supports_structure(optimizer::StochasticProgramOptimizerType, structure::AbstractStochasticStructure)

Return a `Bool` indicating whether `optimizer` supports the stochastic `structure`. That is, `load_structure!(optimizer, structure)` will not throw `UnsupportedStructure`
"""
function supports_structure(optimizer::StochasticProgramOptimizerType, structure::AbstractStochasticStructure)
    return false
end
"""
    check_loadable(optimizer::AbstractStructuredOptimizer, structure::AbstractStochasticStructure)

Throws an `UnloadableStructure` exception if `structure` is not loadable by `optimizer`.

    See also: [`load_structure!`](@ref)
"""
function check_loadable end
"""
    load_structure!(optimizer::AbstractStructuredOptimizer, structure::AbstractStochasticStructure, x₀::AbstractVector)

Instantiate the `optimizer` with the stochastic program represented in memory by the given `structure` and initial decision `x₀`.

See also: [`optimize!`](@ref)
"""
function load_structure! end
"""
    restore_structure!(optimizer::AbstractStructuredOptimizer)

Restore the stochastic program to the state it was in before a call to `optimize!`

See also: [`load_structure!`](@ref)
"""
function restore_structure! end
"""
    optimize!(optimizer::AbstractStructuredOptimizer)

Start the solution procedure for `optimizer` after a call to [`load_structure!`](@ref).

See also: [`load_structure!`](@ref)
"""
function optimize! end
"""
    optimizer_name(optimizer::AbstractStructuredOptimizer)

Optionally, return a string identifier of `AbstractStructuredOptimizer`.
"""
function optimizer_name(::AbstractStructuredOptimizer)
    return "SolverName() attribute not implemented by the optimizer."
end
"""
    master_optimizer(optimizer::AbstractStructuredOptimizer)

Return a MOI optimizer constructor
"""
function master_optimizer end
"""
   subproblem_optimizer(optimizer::AbstractStructuredOptimizer)

Return a MOI optimizer constructor for solving subproblems
"""
function subproblem_optimizer end
"""
    num_iterations(optimizer::AbstractStructuredOptimizer)

Return the number of iterations ran by the structured optimizer

See also: [`optimize!`](@ref)
"""
function num_iterations end
# Sample-based solver interface
# ========================== #
"""
    load_model!(optimizer::AbstractSampledOptimizer, model::StochasticModel, x₀::AbstractVector)

Instantiate the `optimizer` with the stochastic model and initial decision `x₀`.

See also: [`optimize!`](@ref)
"""
function load_model! end
"""
    optimizer_name(optimizer::AbstractSampledOptimizer)

Optionally, return a string identifier of `AbstractSampledOptimizer`.
"""
function optimizer_name(::AbstractSampledOptimizer)
    return "SolverName() attribute not implemented by the optimizer."
end
"""
    optimal_instance(optimizer::AbstractSampledOptimizer)

Return a stochastic programming instance of the stochastic model after a call to [`optimize!`](@ref).
"""
function optimal_instance end


text:
cd(joinpath(@__DIR__, "../"))
using Pkg
pkg"activate ."

using VectorAutoregressions
using CSV, DataFrames, DataFramesMeta
using Dates

cd(@__DIR__)
datadf = CSV.read("lutkepohl_data.csv")

# Take the subset from 1960q1 to 1978q4, and use log-difference for each series
Y = @linq datadf |>
      where(:date .< Date(1979, 3, 1)) |>
      select(
          investment = diff(log.(:invest)),
          income = diff(log.(:income)),
          consumption = diff(log.(:cons)),
      )

# Create quarterly dummies
q1 = repeat([1, 0, 0, 0], Int(size(datadf, 1) / 4))
q2 = repeat([0, 1, 0, 0], Int(size(datadf, 1) / 4))
q3 = repeat([0, 0, 1, 0], Int(size(datadf, 1) / 4))
q4 = repeat([0, 0, 0, 1], Int(size(datadf, 1) / 4))
qdums = hcat(q2, q3, q4) # Array
qdums_df = DataFrame(q2 = q2, q3 = q3, q4 = q4) # DataFrame
qdums_df = qdums_df[18:end, :]

lags = 2
constant = true
trend = false

v = VAR(Y, lags, trend = trend, exog = qdums_df)


text:
# TODO(tim): support NA
# TODO(tim): support Nullable
# TODO(tim): support default value for missing keys

struct CategoricalDiscretizer{N,D} <: AbstractDiscretizer{N,D}
    n2d :: Dict{N,D}    # maps natural to discrete
    d2n :: Dict{D,N}    # maps discrete to natural
end

function CategoricalDiscretizer(natural_to_discrete::Dict{N,D}) where {N,D}
    d2n = Dict{D,N}()
    for (k,v) in natural_to_discrete
        d2n[v] = k
    end
    CategoricalDiscretizer{N,D}(natural_to_discrete, d2n)
end
function CategoricalDiscretizer(data::AbstractArray{N}, ::Type{D}=Int) where {N, D<:Integer}
    # build a label mapping N -> D <: Integer
    i = zero(D)
    n2d = Dict{N, D}()
    for x in data
        if !haskey(n2d,x)
            n2d[x] = (i += 1)
        end
    end
    CategoricalDiscretizer(n2d)
end

supports_encoding(cd::CategoricalDiscretizer{N,D}, x::N) where {N,D} = haskey(cd.n2d, x)
supports_decoding(cd::CategoricalDiscretizer{N,D}, d::D) where {N,D} = 1 ≤ d ≤ nlabels(cd)

encode(cd::CategoricalDiscretizer{N,D}, x::N) where {N,D} = cd.n2d[x]::D
encode(cd::CategoricalDiscretizer{Symbol,D}, x::AbstractString) where {D} = cd.n2d[Symbol(x)]::D
encode(cd::CategoricalDiscretizer{N,D}, x) where {N,D} = cd.n2d[convert(N,x)]::D
function encode(cd::CategoricalDiscretizer{N,D}, data::AbstractArray) where {N,D}
    arr = Array{D}(undef, length(data))
    for (i,x) in enumerate(data)
        arr[i] = encode(cd, x)
    end
    reshape(arr, size(data))
end

decode(cd::CategoricalDiscretizer{N,D}, x::D) where {N,D} = cd.d2n[x]::N
decode(cd::CategoricalDiscretizer{N,D}, x) where {N,D} = cd.d2n[convert(D,x)]::N
function decode(cd::CategoricalDiscretizer{N,D}, data::AbstractArray{D}) where {N,D}
    arr = Array{N}(undef, length(data))
    for (i,d) in enumerate(data)
        arr[i] = decode(cd, d)
    end
    reshape(arr, size(data))
end

nlabels(cd::CategoricalDiscretizer) = length(cd.n2d)


text:
module CoordinateGA

using Random
using LinearAlgebra: dot
using ..GAFramework
import ..GAFramework: fitness, crossover!, mutation!, selection, randcreature

export CoordinateCreature, FunctionModel, SumModel

# This is a CoordinateCreature, T is some coordinate type like Vector
# and objvalue is some objective value
struct CoordinateCreature{T}
    value :: T
    objvalue :: Float64
end

# Negative, since we are minimizing the objective value
fitness(x::CoordinateCreature{T}) where {T} = -x.objvalue

# We can use multiple models to optimize this creature.
# To demonstrate, below we use the FunctionModel model
# as well as the SumModel

# The following are some crossover functions for the CoordinateCreature
# Each model can choose its own crossover function
struct AverageCrossover end
function crossover!(::AverageCrossover, z::CoordinateCreature{T},
    x::CoordinateCreature{T}, y::CoordinateCreature{T},
    st::GAState, aux, rng::AbstractRNG) where {T}
    z.value .= 0.5 .* (x.value .+ y.value)
    CoordinateCreature(z.value, st.model)
end

struct SinglePointCrossover end
function crossover!(::SinglePointCrossover, z::CoordinateCreature{T},
    x::CoordinateCreature{T}, y::CoordinateCreature{T},
    st::GAState, aux, rng::AbstractRNG) where {T}
    N = length(x.value)
    i = rand(rng, 1:N)
    if rand(rng) < 0.5
        z.value[1:i] = x.value[1:i]
        z.value[i+1:end] = y.value[i+1:end]
    else
        z.value[1:i] = y.value[1:i]
        z.value[i+1:end] = x.value[i+1:end]
    end
    CoordinateCreature(z.value, st.model)
end

struct TwoPointCrossover end
function crossover!(::TwoPointCrossover, z::CoordinateCreature{T},
    x::CoordinateCreature{T}, y::CoordinateCreature{T},
    st::GAState, aux, rng::AbstractRNG) where {T}
    N = length(x.value)
    i,j = rand(rng, 1:N, 2)
    i,j = i > j ? (j,i) : (i,j)
    if rand(rng) < 0.5
        z.value[:] = x.value
        z.value[i+1:j] = y.value[i+1:j]
    else
        z.value[:] = y.value
        z.value[i+1:j] = x.value[i+1:j]
    end
    CoordinateCreature(z.value, st.model)
end

printfitness(curgen::Int, x::CoordinateCreature{T}) where {T} =
    println("curgen: $curgen value: $(x.value) obj. value: $(x.objvalue)")

# The following is the SumModel
# It finds a d-dimensional coordinate x such that dot(x, 1:length(x)) ≈ y
"""
    model = SumModel(5, 42.)
    state = GAState(model, ngen=500, npop=6_000, elite_fraction=0.1,
                       params=Dict(:mutation_rate=>0.9, :print_fitness_iter=>1))
    ga!(state)
"""
struct SumModel{T} <: GAModel
    d::Int
    target::T
end

function CoordinateCreature(value::Vector{T}, m::SumModel{T}) where {T}
    objval = Float64(abs(dot(value, 1:length(value)) - m.target))
    CoordinateCreature(value, objval)
end

function randcreature(m::SumModel{T}, aux, rng::AbstractRNG) where {T}
    value = rand(rng, m.d)
    CoordinateCreature(value, m)
end

function mutation!(x::CoordinateCreature,
    st::GAState{SumModel{T}}, aux, rng::AbstractRNG) where {T}
    gp = st.params
    if rand(rng) < get(gp, :mutation_rate, 0.1)
        x.value[rand(rng, 1:length(x.value))] += st.model.target * randn(rng,T)
        CoordinateCreature(x.value, st.model)
    else
        x
    end
end

function crossover!(z::CoordinateCreature,
    x::CoordinateCreature, y::CoordinateCreature,
    st::GAState{SumModel{T}}, aux, rng::AbstractRNG) where {T}
    crossover!(TwoPointCrossover(), z, x, y, st, aux, rng)
end

# The following is the FunctionModel
# It's a more complicated model than above
# It minimizes the objective value using a given objective function
"""
# E.g.
    model = FunctionModel(x -> abs(x[1] - 20.0), [-200.0], [200.0])
    state = GAState(model, ngen=500, npop=6_000, elite_fraction=0.1,
                       params=Dict(:mutation_rate=>0.9, :print_fitness_iter=>1))
    ga!(state)

    type T has to have properties
        y-x :: T
        0.25 .* (x+y) :: T
        randn(T) :: T
        z + 0.25*(y-x)*randn(T) :: T

"""
struct FunctionModel{F,T} <: GAModel
    f::F
    xmin::T
    xmax::T
    xspan::T # xmax-xmin
    clamp::Bool
end
function FunctionModel(f::F,xmin,xmax,clamp::Bool=true) where {F}
    xmin,xmax = promote(xmin,xmax)
    ET = eltype(xmin)
    N = length(xmin)
    xspan = xmax .- xmin
    # check that f(xmin), f(xmax) can be converted to Float64 without error
    z1 = Float64(f(xmin))
    z2 = Float64(f(xmax))
    # and that f(xspan), f(xmin), and f(xmax) has sane values maybe
    # z1!=Inf && z2!=Inf && !isnan(z1) && !isnan(z2) ||
    #    error("f(xmin) or f(xmax) objective function is either NaN or Inf")
    all(xspan .>= zero(ET)) || error("xmax[i] < xmin[i] for some i")
    FunctionModel{F,typeof(xspan)}(f,xmin,xmax,xspan,clamp)
end

CoordinateCreature(value::T, m::FunctionModel{F,T}) where {F,T} =
    CoordinateCreature(value, m.f(value))

function randcreature(m::FunctionModel{F,T}, aux, rng::AbstractRNG) where {F,T}
    ET = eltype(T)
    if T <: Vector
        xvalue = m.xmin .+ m.xspan .* rand(rng, ET, length(m.xspan))
    else
        xvalue = m.xmin .+ m.xspan .* rand(rng, T)
    end
    CoordinateCreature(xvalue, m)
end

function crossover!(z::CoordinateCreature{T},
    x::CoordinateCreature{T}, y::CoordinateCreature{T},
    st::GAState{FunctionModel{F,T}}, aux, rng::AbstractRNG) where {F,T}
    crossover!(TwoPointCrossover(), z, x, y, st, aux, rng)
end

# Mutate over all dimensions
function mutatenormal!(x::CoordinateCreature{T}, temp::Real,
    model::FunctionModel{F,T}, rng::AbstractRNG) where {F,T}
    x.value .+= temp .* model.xspan .* randn(rng,length(x.value))
    model.clamp && (x.value .= clamp.(x.value, model.xmin, model.xmax))
    CoordinateCreature(x.value, model)
end

# Mutate through a single dimension
function mutatenormaldim!(x::CoordinateCreature{T}, temp::Real, dim::Integer,
    model::FunctionModel{F,T}, rng::AbstractRNG) where {F,T}
    ET = eltype(T)
    x.value[dim] += temp * model.xspan[dim] * randn(rng,ET)
    model.clamp && (x.value[dim] = clamp(x.value[dim], model.xmin[dim], model.xmax[dim]))
    CoordinateCreature(x.value, model)
end

function mutation!(x::CoordinateCreature{T},
    st::GAState{FunctionModel{F,T}}, aux, rng::AbstractRNG) where {F,T}
    gp = st.params
    if rand(rng) < get(gp, :mutation_rate, 0.1)
        if rand(rng) < get(gp, :sa_rate, 0.0)
            sa(x, st.model, gp[:sa_k], gp[:sa_lambda],
                gp[:sa_maxiter], st.curgen, aux, rng)
        else
            N = length(x.value)
            mutatenormaldim!(x, 0.1, rand(rng, 1:N), st.model, rng)
        end
    else
        x
    end
end

# export sa,satemp,saprob,mutatenormal
function sa(x::CoordinateCreature{T}, model::FunctionModel{F,T},
    k::Real, lambda::Real, maxiter::Integer, curgen::Int,
    aux, rng::AbstractRNG) where {F,T}
    N = length(x.value)
    y = x
    numnoups = 0 # number of consecutive
    #curgen = 0
    iter = curgen*maxiter + 0
    while true
        temp = satemp(iter, k, lambda) + numnoups/maxiter
        #temp = satemp(iter, k, lambda) + lambda * log(1+numnoups)
        dim = rand(rng,1:N)
        yvdim_old = y.value[dim]
        yov_old = y.objvalue
        fitness_old = fitness(y)
        y = mutatenormaldim!(y, temp, dim, model, rng)
        diffe = fitness(y) - fitness_old
        if diffe >= 0
        elseif rand(rng) < saprob(diffe, temp)
        else
            y.value[dim] = yvdim_old
            y = CoordinateCreature(y.value, yov_old)
        end
        #println("temp: $temp newy: $(newy.objvalue) diffe: $diffe prob: $(saprob(diffe, temp))")
        #numnoups = ifelse(diffe > 0, 0, ifelse(diffe < 0, numnoups+1, numnoups))
        numnoups = ifelse(diffe > 0, 0, numnoups + 1)
        iter += 1
        iter > curgen*maxiter + maxiter && break
    end
    #println("pre: $(x.objvalue) post: $(y.objvalue)")
    y
end

satemp(iter::Integer, k::Real, lambda::Real) = k * exp(-lambda * iter)
# diff = newscore - oldscore
saprob(diff::Real, iter::Integer, k::Real, lambda::Real) =
    exp(diff / satemp(iter, k, lambda))
saprob(diff::Real, temp::Real) = exp(diff / temp)

function selection(pop::Vector, n::Integer,
    st::GAState{FunctionModel{F,T}}, rng::AbstractRNG) where {F,T}
    selection(TournamentSelection(2), pop, n, st, rng)
end

end # CoordinateGA


text:
module BCBI_base

using Pkg

export install_all,
       check_installed,
       using_all

const plotting_pkgs = [ "VegaLite", "PGFPlotsX", "PlotlyJS", "Gadfly"]

const datasets_pkg = ["RDatasets", "VegaDatasets"]

const external_pkgs = ["RCall", "PyCall", "Pandas", "PyPlot", "Seaborn"]

const problematic_pkgs = ["GLMNet", "EzXML", "Lasso"]

const predictmd_extra_deps = [  "Atom",
                                "BenchmarkTools",
                                "CSVFiles",
                                "CUDAapi",
                                "Combinatorics",
                                "Documenter",
                                "FluxJS",
                                "GLM",
                                "GPUArrays",
                                "Knet",
                                "Literate",
                                "MLBase",
                                "MLDatasets",
                                "Metalhead",
                                "NBInclude",
                                "NumericalIntegration",
                                "PGFPlots",
                                "PGFPlotsX",
                                "ProgressMeter",
                                "RDatasets",
                                "Requires",
                                "ValueHistories",
                                "Weave"]

const base_pkgs = [ "MySQL",
                    "BioServices",
                    "ClassImbalance",
                    "Clustering",
                    "Compat",
                    "CSV",
                    "DataFrames",
                    "DecisionTree",
                    "GLM",
                    "HTTP",
                    "HypothesisTests",
                    "IJulia",
                    "JLD",
                    "JLD2",
                    "JuliaDB",
                    "LIBSVM",
                    "LightXML",
                    "MixedModels",
                    "Revise",
                    "ROCAnalysis",
                    "StatsBase",
                    "Query"
                    ]

"""
    clone_pkgs
Dictionary of package name and URL to clone.
We use it for unregistered and packages that we wish to checkout the master branch.
For registered packages, we perform `Pkg.clone` instead of `Pkg.checkout` as the later may fail
if the lates tagged version cannot be resoved
"""
const clone_pkgs =Dict("ARules"=>"https://github.com/bcbi/ARules.jl",
                       "ScikitLearn"=>"https://github.com/cstjean/ScikitLearn.jl.git",
                       "BioMedQuery"=>"https://github.com/bcbi/BioMedQuery.jl.git",
                       "PredictMD"=>"https://github.com/bcbi/PredictMD.jl")

"""
    checkout_pkgs
Dictionary of package name to brach
"""
const checkout_pkgs=Dict("PredictMD"=>"master")

"""
    add(pkgs = BCBI_base.base_pkgs)

Call `Pkg.add` and `using` on list of desired packages
"""
function add(pkgs = BCBI_base.base_pkgs)

    Pkg.update()

    failed_pkgs = Vector{String}()

    for pkg in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")

        try
            println("* Adding")
            Pkg.add(pkg)
        catch
            warn("Pkg.add failed")
            push!(failed_pkgs, pkg)
            continue
        end

    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    add_and_precompile(pkgs = BCBI_base.base_pkgs)

Call `Pkg.add` and `using` on list of desired packages
"""
function add_and_precompile(pkgs = BCBI_base.base_pkgs)

    Pkg.update()

    failed_pkgs = Vector{String}()

    for pkg in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")

        try
            println("* Adding")
            Pkg.add(pkg)
        catch
            warn("Pkg.add failed")
            push!(failed_pkgs, pkg)
            continue
        end

        try
            println("* Using")
            pkgsym = Symbol(pkg)
            eval(:(using $pkgsym))
        catch
            warn("using pkg failed")
            push!(failed_pkgs, pkg)
            Pkg.rm(pkg)
        end

    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    clone(pkgs = BCBI_base.clone_pkgs)

Call `Pkg.clone` and `using` on `Dict("name"->"url")` of desired packages
"""
function clone(pkgs = BCBI_base.clone_pkgs)

    failed_pkgs = Vector{String}()

    for (pkg, url) in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")

        try
            v = Pkg.installed(pkg)
            if typeof(v) != VersionNumber
                println("* Clonning")
                Pkg.clone(url)
            end
        catch
            println("* Clonning")
            Pkg.clone(url)
        end

    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    clone_and_precompile(pkgs = BCBI_base.clone_pkgs)

Call `Pkg.clone` and `using` on `Dict("name"->"url")` of desired packages
"""
function clone_and_precompile(pkgs = BCBI_base.clone_pkgs)

    failed_pkgs = Vector{String}()

    for (pkg, url) in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")

        try
            v = Pkg.installed(pkg)
            if typeof(v) != VersionNumber
                println("* Clonning")
                Pkg.clone(url)
            end
        catch
            println("* Clonning")
            Pkg.clone(url)
        end

        try
            println("* Using")
            pkgsym = Symbol(pkg)
            eval(:(using $pkgsym))
        catch
            warn("using pkg failed")
            push!(failed_pkgs, pkg)
            Pkg.rm(pkg)
        end

    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    checkout(pkgs = Dict())

Call `Pkg.checkout` and `using` on `Dict("name"->"branch")` of desired packages. It assumes that
package is already added
"""
function checkout(pkgs = BCBI_base.checkout_pkgs)

    failed_pkgs = Vector{String}()

    for (pkg, branch) in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")
        println("* Checkout ", branch)

        try
            Pkg.checkout(pkg, branch)
        catch
            warn("checking out branch failed")
            push!(failed_pkgs, pkg)
        end
    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
checkout_and_precompile(pkgs = Dict())

Call `Pkg.checkout` and `using` on `Dict("name"->"branch")` of desired packages. It assumes that
package is already added
"""
function checkout_and_precompile(pkgs = BCBI_base.checkout_pkgs)
    failed_pkgs = Vector{String}()

    for (pkg, branch) in pkgs
        println("--------------------------------")
        println("Package: ", pkg)
        println("--------------------------------")
        println("* Checkout ", branch)
        try
            Pkg.checkout(pkg, branch)
            println("* Using")
            pkgsym = Symbol(pkg)
            eval(:(using $pkgsym))
        catch
            warn("using pkg failed")
            push!(failed_pkgs, pkg)
            Pkg.rm(pkg)
        end
    end

    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    install_all()
Run checkout/add/clone functions
"""
function install_all()
    add(base_pkgs)
    add(plotting_pkgs)
    add(datasets_pkg)
    add(external_pkgs)
    clone(clone_pkgs)
    checkout(checkout_pkgs)
end

"""
    check_installed()
Print and return list of missing "desired" packages
"""
function check_installed()
    desired = vcat(base_pkgs, plotting_pkgs, external_pkgs, datasets_pkg, collect(keys(clone_pkgs)), collect(keys(checkout_pkgs)))
    installed = collect(keys(Pkg.installed()))
    println("--------------------------------")
    println("Missing desired packages:")
    println("--------------------------------")
    miss_pkgs = setdiff(desired, installed)
    println(miss_pkgs)

    return miss_pkgs
end

"""
    using_all()
Run `using` an all lists of packages
"""
function using_all()
    desired = vcat(base_pkgs, plotting_pkgs, external_pkgs, datasets_pkg, collect(keys(clone_pkgs)), collect(keys(checkout_pkgs)))
    failed_pkgs = Vector{String}()
    for pkg in desired
        try
            pkgsym = Symbol(pkg)
            eval(:(using $pkgsym))
        catch
            warn("using pkg failed")
            push!(failed_pkgs, pkg)
        end
    end
    println("--------------------------------")
    println("Failed packages: ", length(failed_pkgs))
    map(x->println(x), failed_pkgs)
    println("--------------------------------")

    return failed_pkgs
end

"""
    clean_up_failing()
Removes packages that fail to precompile
"""
function clean_up()
    failed_pkgs = using_all()
    for pkg in failed_pkgs
        println("--------------------------------")
        println("Removing package: ", pkg)
        println("--------------------------------")
        Pkg.rm(pkg)
    end
end

end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6f9be298-5374-11ec-3b9d-fd74a32d6b53
using Colors, PlutoUI, Manopt, Manifolds

# ╔═╡ e701b513-4ad9-474c-9d79-60a24f545ebd
md"""
# Bezier curves and their acceleration

This tutorial illustrates how Bézier curves are generalized to manifolds and how to
minimize their acceleration, i.e. how to get a curve that is as straight or as geodesic
while fulfilling constraints

This example also illustrates how to apply the minimization on the corresponding [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html) manifold using a [`gradient_descent`](https://manoptjl.org/stable/solvers/gradient_descent.html) with [`ArmijoLinesearch`](https://manoptjl.org/stable/plans/index.html#Manopt.ArmijoLinesearch).
"""

# ╔═╡ 087fb979-d68d-4ee7-9be5-fad992fdb792
md"and we define some colors from [Paul Tol](https://personal.sron.nl/~pault/)"

# ╔═╡ 85d8985b-f5ee-4cba-bcc5-528049798330
begin
    black = RGBA{Float64}(colorant"#000000")
    TolVibrantBlue = RGBA{Float64}(colorant"#0077BB")
    TolVibrantOrange = RGBA{Float64}(colorant"#EE7733")
    TolVibrantMagenta = RGBA{Float64}(colorant"#EE3377")
    TolVibrantCyan = RGBA{Float64}(colorant"#33BBEE")
    TolVibrantTeal = RGBA{Float64}(colorant"#009988")
end;

# ╔═╡ e62e36e9-19ce-4a7d-8bda-5634813c1ee9
md"Let's also set up a few parameters"

# ╔═╡ 18d378dd-abc6-4fa6-a2c9-1cc09e4cec46
begin
    geo_pts = collect(range(0.0, 1.0; length=101))
    bezier_pts = collect(range(0.0, 3.0; length=201))
    camera_position = (-1.0, -0.7, 0.3)
    #render asy yes/no. If not, images included w/ markdown are assumed to be prerendered
    render_asy = false
    localpath = join(splitpath(@__FILE__)[1:(end - 1)], "/") # remove file to get files folder
    image_prefix = localpath * "/bezier"
    @info image_prefix
end;

# ╔═╡ dee4eb65-145c-4b2b-98a7-ae3a7045f24b
md"""
We finally load our data, see [`artificial_S2_composite_bezier_curve`](https://manoptjl.org/stable/helpers/data.html#Manopt.artificial_S2_composite_bezier_curve-Tuple{}), a composite Bezier curve consisting of 3 segments on the Sphere. The middle segment consists of the control points
"""

# ╔═╡ 2b485956-459d-40e7-9364-79ea77073804
begin
    B = artificial_S2_composite_bezier_curve()
    b = B[2].pts
end

# ╔═╡ 98c925ac-9ac8-42f3-9855-a01e5eaac350
M = Sphere(2)

# ╔═╡ d22e825e-d8ff-44fe-8b06-b14c09765180
md"""
On Euclidean spaces Bézier curves of these ``n=4`` so called control points like this segment yield polynomials of degree``3``.
The resulting curve ``γ: [0,1] → ℝ^m`` is called [Bezier curve](https://en.wikipedia.org/wiki/Bézier_curve) or Bézier spline and is named after [Piérre Bezier](https://en.wikipedia.org/wiki/Pierre_Bézier) (1910–1999).
They can be evaluated by the de Casteljau algorithm by evaluating line segments between points.
While it is not easy to evaluate polynomials on a manifold, evaluating line segments generalizes to the evaluation of [`shortest_geodesic`](https://juliamanifolds.github.io/Manifolds.jl/latest/interface.html#ManifoldsBase.shortest_geodesic-Tuple{Manifold,Any,Any})s.
We will illustrate this using these points ``b=(b_1,b_2,b_3,b_4)`` on the [`Sphere`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/sphere.html) ``\mathbb S^2``.
Let's evaliuate this at the point ``t=\frac{1}{4}∈[0,1]``. We first compute
"""

# ╔═╡ 157d13c8-7311-4832-a9ff-da22125b2bbf
begin
    t = 0.66
    pts1 = shortest_geodesic.(Ref(M), b[1:3], b[2:4], Ref(t))
end

# ╔═╡ 5af58241-50c7-464e-bfa3-53819e605793
md"""
We obtain 3 points on the geodesics connecting the control points. Repeating this again twice
"""

# ╔═╡ b372fcfc-cad0-48bc-b759-54d15bd2e28b
begin
    pts2 = shortest_geodesic.(Ref(M), pts1[1:2], pts1[2:3], Ref(t))
    p = shortest_geodesic(M, pts2[1], pts2[2], t)
end

# ╔═╡ 62115e3f-cdfe-464c-91a1-eb80e674a610
md"""we obtain the point on the Bézier curve ``c(t)``."""

# ╔═╡ 8653ca13-21e9-4049-a81e-02b6047d2c66
begin
    curves1 = [shortest_geodesic(M, b[i], b[i + 1], geo_pts) for i in 1:3]
    curves2 = [shortest_geodesic(M, pts1[i], pts1[i + 1], geo_pts) for i in 1:2]
    bezier_curve = [shortest_geodesic(M, pts2[1], pts2[2], geo_pts)]
end;

# ╔═╡ 05612a12-1080-4b79-b505-198e079a4b2a
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/Casteljau-illustr.asy";
        points=[b, pts1, pts2, [p]],
        curves=[curves1..., curves2..., bezier_curve..., de_casteljau(M, B[2], geo_pts)],
        colors=Dict(
            :points =>
                [TolVibrantBlue, TolVibrantCyan, TolVibrantTeal, TolVibrantOrange],
            :curves => [
                TolVibrantBlue,
                TolVibrantBlue,
                TolVibrantBlue,
                TolVibrantCyan,
                TolVibrantCyan,
                TolVibrantTeal,
                black,
            ],
        ),
        dot_size=3.5,
        line_widths=[1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.5],
        camera_position=camera_position,
    )
    render_asymptote(image_prefix * "/Casteljau-illustr.asy"; render=2)
end;

# ╔═╡ bdb85fb8-0707-47d1-95ee-9bf2e56500f5
md"""This image summarizes Casteljaus algorithm:"""

# ╔═╡ 38b8f8af-ba81-4d3b-83b7-56b0d41ea708
PlutoUI.LocalResource(image_prefix * "/Casteljau-illustr.png")

# ╔═╡ d6759a8e-a705-4d95-a8da-e172ada644d8
md"""
From the control points (blue) and their geodesics, ont evaluation per geodesic yields three interims points (cyan), their two successive geodeics another two points (teal) and at its geodesic at ``t=0.66`` we obtain the point on the curve.

In Manopt.jl, to evaluate a Bézier curve knowing its [`BezierSegment`](https://manoptjl.org/stable/functions/bezier.html#Manopt.BezierSegment), use [`de_casteljau`](https://manoptjl.org/stable/functions/bezier.html#Manopt.de_casteljau-Tuple{AbstractManifold,%20Vararg{Any,%20N}%20where%20N}).

There are a few nice observations to make, that hold also for these Bézier curves on manifolds:
* The curve starts in the first controlpoint ``b_0`` and ends in the last controlpoint ``b_3``
* The tangent vector to the curve at the start ``\dot c(0)`` is equal to ``\frac{1}{3}\log_{b_0}b_1 = \dot γ_{b_0,b_1}(0)``, where ``γ_{a,b}`` denotes the shortest geodesic between ``a`` and ``b``.
* The tangent vector to the curve at the end ``\dot c(1)`` is equal to ``-\frac{1}{3}\log_{b_3}b_2 = -γ_{b_3,b_2}(0) = \dot γ_{b_2,b_3}(1)``.
* the curve is differentiable.

For more details on these properties, see for example [^PopielNoakes2007].
"""

# ╔═╡ 7c6947a0-7d0b-4459-b999-58ea4f2b0169
md"""
## Composite Bézier curves

With the properties of a single Bézier curve, also called Bézier segment, we can “stitch” curves together. Let ``a_0,…,a_n`` and ``b_0,…,b_m`` be two sets of controlpoints for the Bézier segments ``c(t)`` and ``d(t)``, respectively.
We define the composite Bézier curve by ``B(t) = \begin{cases} c(t) & \text{ if } 0\leq t < 1, \\ d(t-1) & \text{ if } 1\leq t \leq 2,\end{cases}`` where ``t∈[0,2]``.
This can of course be generalised straight forward to more than two cases.
With the properties from the previous section we can now state that

* the curve ``B(t)`` is continuous if ``c(1)=d(0)`` or in other words ``a_n=b_0``
* the curve ``B(t)`` is differentiable if additionally ``\dot c(1)=\dot d(0)`` or in
other words ``-\log_{a_n}a_{n-1} = \log_{b_0}b_1``. This is equivalent to ``a_n=b_0 = \gamma_{a_{n-1}b_1}(\tfrac{1}{2})``.

One nice interpretation of the last characterization is, that the tangents ``\log_{a_n}a_{n-1}`` and ``\log_{b_0}b_1`` point into opposite directions.
For a continuous curve, the first point of every segment (except for the first segment) can be ommitted, for a differentiable curve the first two points (except for the first segment) can be ommitted.
You can reduce storage by calling [`get_bezier_points`](https://manoptjl.org/stable/functions/bezier.html#Manopt.get_bezier_points), though for econstruciton with [`get_bezier_segments`](https://manoptjl.org/stable/functions/bezier.html#Manopt.get_bezier_segments-Union{Tuple{P},%20Tuple{AbstractManifold,%20Vector{P},%20Any},%20Tuple{AbstractManifold,%20Vector{P},%20Any,%20Symbol}}%20where%20P) you also need [`get_bezier_degrees`](https://manoptjl.org/stable/functions/bezier.html#Manopt.get_bezier_degrees-Tuple{AbstractManifold,%20AbstractVector{var%22#s56%22}%20where%20var%22#s56%22%3C:BezierSegment}).
The reduced storage is represented as an array of points, i.e. an element of the corresponding [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html).

For the three segment example from the beginning this looks as follows
"""

# ╔═╡ 3f1d0989-2f97-4945-82f8-22e9c7bda6d2
render_asy && begin
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-composite-curve.asy";
        curves=[de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B, true),
                get_bezier_junction_tangent_vectors(M, B),
            )
        ]],
        colors=Dict(
            :curves => [black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-composite-curve.asy"; render=2)
end;

# ╔═╡ 41ba8ac2-5a81-46f5-9269-404cc6a1d287
PlutoUI.LocalResource(image_prefix * "/Bezier-composite-curve.png")

# ╔═╡ e7fc05b3-a198-4a64-b62a-43cdf8f56238
md"""
## Minimizing the acceleration of a composite Bézier curve
The motivation to minimize the acceleration of the composite Bézier curve is, that the curve should get “straighter” or more geodesic like.
If we discretize the curve ``B(t)`` with its control points denoted by ``b_{i,j}`` for the ``j``th note in the ``i``th segment, the discretized model for equispaced ``t_i``, ``i=0,…,N`` in the domain of ``B`` reads[^BergmannGousenbourger2018]

```math
A(b) :eqq\sum_{i=1}^{N-1}\frac{\mathrm{d}^2_2 \bigl[ B(t_{i-1}), B(t_{i}), B(t_{i+1}) \bigr]}{\Delta_t^3},
```

where $\mathrm{d}_2$ denotes the second order finite difference using the mid point approach, see [`costTV2`](@ref)[^BacakBergmannSteidlWeinmann2016],

```math
d_2(x,y,z) := \min_{c ∈ \mathcal C_{x,z}} d_{\mathcal M}(c,y),\qquad x,y,z∈\mathcal M.
```

Another model is based on logarithmic maps, see [^BoumalAbsil2011], but that is not considered here.
An advantage of the model considered here is, that it only consist of the evaluation of geodesics.
This yields a gradient of ``A(b)`` with respect to ``b`` [`adjoint_Jacobi_field`](https://manoptjl.org/stable/functions/Jacobi_fields.html#Manopt.adjoint_Jacobi_field)s. The following image shows the negative gradient (scaled)
"""

# ╔═╡ 7b545a1e-b9c3-4c57-a5ee-3a629407b67a
render_asy && begin
    gradFullB = Manopt._grad_acceleration_bezier(
        M,
        get_bezier_points(M, B, :differentiable),
        [3, 3, 3],
        collect(range(0.0, 3.0; length=151)),
    )
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-composite-curve-gradient.asy";
        curves=[de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B), get_bezier_inner_points(M, B)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_points(M, B, :continuous),
                -0.05 .* get_bezier_points(M, gradFullB, :continuous),
            )
        ]],
        colors=Dict(
            :curves => [black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantOrange],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-composite-curve-gradient.asy"; render=2)
end;

# ╔═╡ 5157ac76-d359-4f18-883d-8afbe78d3b4d
PlutoUI.LocalResource(image_prefix * "/Bezier-composite-curve-gradient.png")

# ╔═╡ 7640628b-e232-4ba3-94aa-d8a00218ff6a
md"""
In the following we consider two cases: Interpolation, which fixes the junction and end points of ``B(t)``
and approximation, where a weight and a dataterm are additionally introduced.
"""

# ╔═╡ da970042-8326-4d17-b440-ab1960f10283
md"""
### Interpolation

For interpolation, the junction points are fixed and their gradient entries are hence set to zero.
After transferring to the already mentioned [`PowerManifold`](https://juliamanifolds.github.io/Manifolds.jl/latest/manifolds/power.html), we can then perform a  [`gradient_descent`](https://manoptjl.org/stable/solvers/gradient_descent.html)  as follows
"""

# ╔═╡ 7db3e138-6ec1-4a0c-8f3b-279f7251b4c3
begin
    #number of sampling points refers to exactness of approximating d^2
    curve_samples = collect(range(0.0, 3.0; length=151))
    pB = get_bezier_points(M, B, :differentiable)
    N = PowerManifold(M, NestedPowerRepresentation(), length(pB))
    function F(M, pB)
        return cost_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples
        )
    end
    function gradF(M, pB)
        return grad_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples
        )
    end
    x0 = get_bezier_points(M, B, :differentiable)
    pB_opt_ip = gradient_descent(
        N,
        F,
        gradF,
        x0;
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.0001),
        stopping_criterion=StopWhenChangeLess(5 * 10.0^(-7)),
    )
end;

# ╔═╡ 5e8e98f6-8896-462e-b610-4b379c727a72
md"and the resut looks like"

# ╔═╡ 50854504-af68-4d7f-999e-f3054c088fb6
render_asy && begin
    B_opt_ip = get_bezier_segments(M, pB_opt_ip, get_bezier_degrees(M, B), :differentiable)
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-IP-Min.asy";
        curves=[de_casteljau(M, B_opt_ip, bezier_pts), de_casteljau(M, B, bezier_pts)],
        points=[get_bezier_junctions(M, B_opt_ip), get_bezier_inner_points(M, B_opt_ip)],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B_opt_ip, true),
                get_bezier_junction_tangent_vectors(M, B_opt_ip),
            )
        ]],
        colors=Dict(
            :curves => [TolVibrantBlue, black],
            :points => [TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 0.75, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-IP-Min.asy"; render=2)
end;

# ╔═╡ 4f35b3e8-90ca-4d21-a377-007f93eb0efd
PlutoUI.LocalResource(image_prefix * "/Bezier-IP-Min.png")

# ╔═╡ 203a0538-d608-4d3e-b57d-edec20756300
md"""
Where the original courve is shown in black and the interpolating curve with minimized (discretized) acceleration is shown in blue including its junction points (also blue), tangent vectors (light blue) and control points (teal).
"""

# ╔═╡ 3e43f3fe-5c9a-4fb5-9033-db4a34bac5a5
md"""
### Approximation
Similarly if we introduce the junction points as data fixed given ``d_i`` and set (for simplicity) ``p_i=b_{i,0}`` and ``p_{n+1}=b_{n,4}`` and ``λ=3`` in

```math
\frac{λ}{2}\sum_{k=0}^3 d_{\mathcal M}(d_i,p_i)^2 + A(b),
```

then ``λ`` models how important closeness to the data ``d_i`` is.

Then we obtain the folowing code to minimize the acceleration while approximating the original curve
"""

# ╔═╡ 2f9ef252-5046-4af7-aa69-449a04bab26e
begin
    λ = 3.0
    d = get_bezier_junctions(M, B)
    function F2(M, pB)
        return cost_L2_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples, λ, d
        )
    end
    function gradF2(M, pB)
        return grad_L2_acceleration_bezier(
            M.manifold, pB, get_bezier_degrees(M.manifold, B), curve_samples, λ, d
        )
    end
    x1 = get_bezier_points(M, B, :differentiable)
    pB_opt_appr = gradient_descent(
        N,
        F2,
        gradF2,
        x1;
        stepsize=ArmijoLinesearch(1.0, ExponentialRetraction(), 0.5, 0.001),
        stopping_criterion=StopWhenChangeLess(10.0^(-5)),
    )
end

# ╔═╡ 6cfb38a8-c29f-42a6-b642-31c047dd6cb8
render_asy && begin
    B_opt_appr = get_bezier_segments(M, pB_opt_appr, get_bezier_degrees(M, B), :differentiable)
    asymptote_export_S2_signals(
        image_prefix * "/Bezier-Appr-Min.asy";
        curves=[de_casteljau(M, B_opt_appr, bezier_pts), de_casteljau(M, B, bezier_pts)],
        points=[
            d,
            get_bezier_junctions(M, B_opt_appr),
            get_bezier_inner_points(M, B_opt_appr),
        ],
        tangent_vectors=[[
            Tuple(a) for a in zip(
                get_bezier_junctions(M, B_opt_appr, true),
                get_bezier_junction_tangent_vectors(M, B_opt_appr),
            )
        ]],
        colors=Dict(
            :curves => [TolVibrantBlue, black],
            :points => [TolVibrantOrange, TolVibrantBlue, TolVibrantTeal],
            :tvectors => [TolVibrantCyan],
        ),
        camera_position=camera_position,
        arrow_head_size=10.0,
        line_widths=[1.5, 0.75, 1.5],
        dot_size=4.0,
    )
    render_asymptote(image_prefix * "/Bezier-Appr-Min.asy"; render=2)
end;

# ╔═╡ db8110e0-d271-424c-b920-817dc1aab0c9
PlutoUI.LocalResource(image_prefix * "/Bezier-Appr-Min.png")

# ╔═╡ 0855f66e-f9af-4c63-b041-424c9f3fd655
md"""Additionally to the last image, the data points ``d_i`` (junction points of the original curve) are shown in orange, the distance between these and the blue junction points is part of the cost function here."""

# ╔═╡ 67e8e794-a16a-4d39-804a-6bfb350f8247
md"""
The role of ``λ`` can be interpreted as follows: for large values of $λ$, the
minimizer, i.e. the resulting curve, is closer to the original Bézier junction points.
For small ``λ`` the resting curve is closer to a geodesic and the control points are closer to the curve.
For ``λ=0`` _any_ (not necessarily shortest) geodesic is a solution and the problem is ill-posed.
To illustrate the effect of ``λ``, the following image contains 1000 runs for ``λ=10`` in dark currant to ``λ=0.01`` in bright yellow.

![Approximation min Acc](https://manoptjl.org/stable/assets/images/tutorials/Bezier_Approximation_video-Summary-result.png)

The effect of the data term can also be seen in the following video, which starts a little slow and takes about 40 seconds.

![Video of the effect of lambda, the weight of the dataterm](https://manoptjl.org/stable/assets/videos/tutorials/Bezier_Approximation_video-movie.mp4)
"""

# ╔═╡ 446f7276-7ed8-4b5c-a3f1-fdde4fab1479
md"""
## Literature

[^BacakBergmannSteidlWeinmann2016]:
    > Bačák, M., Bergmann, R., Steidl, G. and Weinmann, A.: _A second order nonsmooth
    > variational model for restoring manifold-valued images_,
    > SIAM Journal on Scientific Computations, Volume 38, Number 1, pp. A567–597,
    > doi: [10.1137/15M101988X](https://doi.org/10.1137/15M101988X),
    > arXiv: [1506.02409](https://arxiv.org/abs/1506.02409)
[^BergmannGousenbourger2018]:
    > Bergmann, R. and Gousenbourger, P.-Y.: _A variational model for data fitting on manifolds
    > by minimizing the acceleration of a Bézier curve_.
    > Frontiers in Applied Mathematics and Statistics, 2018.
    > doi: [10.3389/fams.2018.00059](https://dx.doi.org/10.3389/fams.2018.00059),
    > arXiv: [1807.10090](https://arxiv.org/abs/1807.10090)
[^BoumalAbsil2011]:
    > Boumal, N. and Absil, P.-A.: _A discrete regression method on manifolds and its application to data on SO(n)._
    > In: IFAC Proceedings Volumes (IFAC-PapersOnline). Vol. 18. Milano (2011). p. 2284–89.
    > doi: [10.3182/20110828-6-IT-1002.00542](https://dx.doi.org/10.3389/10.3182/20110828-6-IT-1002.00542),
    > web: [www](https://web.math.princeton.edu/~nboumal/papers/Boumal_Absil_A_discrete_regression_method_on_manifolds_and_its_application_to_data_on_SOn.pdf)
[^Casteljau1959]:
    > de Casteljau, P.: _Outillage methodes calcul_, Enveloppe Soleau 40.040 (1959),
    > Institute National de la Propriété Industrielle, Paris.
[^Casteljau1963]:
    > de Casteljau, P.: _Courbes et surfaces à pôles_, Microfiche P 4147-1,
    > André Citroën Automobile SA, Paris, (1963).
[^PopielNoakes2007]:
    > Popiel, T. and Noakes, L.: _Bézier curves and ``C^2`` interpolation in Riemannian
    > manifolds_. Journal of Approximation Theory (2007), 148(2), pp. 111–127.-
    > doi: [10.1016/j.jat.2007.03.002](https://doi.org/10.1016/j.jat.2007.03.002).
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
Manifolds = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
Manopt = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"

[compat]
Colors = "~0.12.8"
Manifolds = "~0.7.1"
Manopt = "~0.3.14"
PlutoUI = "~0.7.22"
"""

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

julia_version = "1.7.2"
manifest_format = "2.0"

[[deps.AbstractFFTs]]
deps = ["LinearAlgebra"]
git-tree-sha1 = "485ee0867925449198280d4af84bdb46a2a404d0"
uuid = "621f4979-c628-5d54-868e-fcf4e3e8185c"
version = "1.0.1"

[[deps.AbstractPlutoDingetjes]]
deps = ["Pkg"]
git-tree-sha1 = "abb72771fd8895a7ebd83d5632dc4b989b022b5b"
uuid = "6e696c72-6542-2067-7265-42206c756150"
version = "1.1.2"

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

[[deps.ArnoldiMethod]]
deps = ["LinearAlgebra", "Random", "StaticArrays"]
git-tree-sha1 = "62e51b39331de8911e4a7ff6f5aaf38a5f4cc0ae"
uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
version = "0.2.0"

[[deps.ArrayInterface]]
deps = ["Compat", "IfElse", "LinearAlgebra", "Requires", "SparseArrays", "Static"]
git-tree-sha1 = "265b06e2b1f6a216e0e8f183d28e4d354eab3220"
uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
version = "3.2.1"

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

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

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

[[deps.ChangesOfVariables]]
deps = ["LinearAlgebra", "Test"]
git-tree-sha1 = "9a1d594397670492219635b35a3d830b04730d62"
uuid = "9e997f8a-9a97-42d5-a9f1-ce6bfc15e2c0"
version = "0.1.1"

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

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

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

[[deps.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 = "dce3e3fea680869eaa0b774b2e8343e9ff442313"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.40.0"

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

[[deps.CovarianceEstimation]]
deps = ["LinearAlgebra", "Statistics", "StatsBase"]
git-tree-sha1 = "bc3930158d2be029e90b7c40d1371c4f54fa04db"
uuid = "587fd27a-f159-11e8-2dae-1979310e6154"
version = "0.2.6"

[[deps.DataAPI]]
git-tree-sha1 = "cc70b17275652eb47bc9e5f81635981f13cea5c8"
uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
version = "1.9.0"

[[deps.DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

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

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

[[deps.DensityInterface]]
deps = ["InverseFunctions", "Test"]
git-tree-sha1 = "80c3e8639e3353e5d2912fb3a1916b8455e2494b"
uuid = "b429d917-457f-4dbc-8f4c-0cc954292b1d"
version = "0.4.0"

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

[[deps.Distributions]]
deps = ["ChainRulesCore", "DensityInterface", "FillArrays", "LinearAlgebra", "PDMats", "Printf", "QuadGK", "Random", "SparseArrays", "SpecialFunctions", "Statistics", "StatsBase", "StatsFuns", "Test"]
git-tree-sha1 = "7f3bec11f4bcd01bc1f507ebce5eadf1b0a78f47"
uuid = "31c24e10-a181-5473-b8eb-7969acd0382f"
version = "0.25.34"

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

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

[[deps.Einsum]]
deps = ["Compat"]
git-tree-sha1 = "4a6b3eee0161c89700b6c1949feae8b851da5494"
uuid = "b7d42ee7-0b51-5a75-98ca-779d3107e4c0"
version = "0.4.1"

[[deps.FillArrays]]
deps = ["LinearAlgebra", "Random", "SparseArrays", "Statistics"]
git-tree-sha1 = "8756f9935b7ccc9064c6eef0bff0ad643df733a3"
uuid = "1a297f60-69ca-5386-bcde-b61e274b549b"
version = "0.12.7"

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

[[deps.Graphs]]
deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
git-tree-sha1 = "92243c07e786ea3458532e199eb3feee0e7e08eb"
uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
version = "1.4.1"

[[deps.HybridArrays]]
deps = ["LinearAlgebra", "Requires", "StaticArrays"]
git-tree-sha1 = "c1d5b1dcdf2140644e1c6beb9ca09fbed601c241"
uuid = "1baab800-613f-4b0a-84e4-9cd3431bfbb9"
version = "0.4.9"

[[deps.Hyperscript]]
deps = ["Test"]
git-tree-sha1 = "8d511d5b81240fc8e6802386302675bdf47737b9"
uuid = "47d2ed2b-36de-50cf-bf87-49c2cf4b8b91"
version = "0.0.4"

[[deps.HypertextLiteral]]
git-tree-sha1 = "2b078b5a615c6c0396c77810d92ee8c6f470d238"
uuid = "ac1192a8-f4b3-4bfe-ba22-af5b92cd3ab2"
version = "0.9.3"

[[deps.IOCapture]]
deps = ["Logging", "Random"]
git-tree-sha1 = "f7be53659ab06ddc986428d3a9dcc95f6fa6705a"
uuid = "b5f81e59-6552-4d32-b1f0-c071b021bf89"
version = "0.2.2"

[[deps.IfElse]]
git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
version = "0.1.1"

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

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

[[deps.InverseFunctions]]
deps = ["Test"]
git-tree-sha1 = "a7254c0acd8e62f1ac75ad24d5db43f5f19f3c65"
uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
version = "0.1.2"

[[deps.IrrationalConstants]]
git-tree-sha1 = "7fd44fd4ff43fc60815f8e764c0f352b83c49151"
uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
version = "0.1.1"

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

[[deps.JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[deps.Kronecker]]
deps = ["LinearAlgebra", "NamedDims", "SparseArrays", "StatsBase"]
git-tree-sha1 = "a51f46415c844dee694cb8b20a3fcbe6dba342c2"
uuid = "2c470bb0-bcc8-11e8-3dad-c9649493f05e"
version = "0.5.1"

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

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

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

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

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

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

[[deps.LogExpFunctions]]
deps = ["ChainRulesCore", "ChangesOfVariables", "DocStringExtensions", "InverseFunctions", "IrrationalConstants", "LinearAlgebra"]
git-tree-sha1 = "be9eef9f9d78cecb6f262f3c10da151a6c5ab827"
uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
version = "0.3.5"

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

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

[[deps.Manifolds]]
deps = ["Colors", "Distributions", "Einsum", "Graphs", "HybridArrays", "Kronecker", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "RecipesBase", "RecursiveArrayTools", "Requires", "SimpleWeightedGraphs", "SpecialFunctions", "StaticArrays", "Statistics", "StatsBase"]
git-tree-sha1 = "3811935fd2549b0f5f9b365d6e7173fcdb98ec9a"
uuid = "1cead3c2-87b3-11e9-0ccd-23c62b72b94e"
version = "0.7.1"

[[deps.ManifoldsBase]]
deps = ["LinearAlgebra", "Markdown"]
git-tree-sha1 = "77a5949567437d185ee929c405e3c6c0768118ea"
uuid = "3362f125-f0bb-47a3-aa74-596ffd7ef2fb"
version = "0.12.9"

[[deps.Manopt]]
deps = ["ColorSchemes", "ColorTypes", "Colors", "DataStructures", "Dates", "LinearAlgebra", "ManifoldsBase", "Markdown", "Random", "Requires", "SparseArrays", "StaticArrays", "Statistics", "Test"]
git-tree-sha1 = "e87c818e19f79444fff1c77880cdb43b51974dfe"
uuid = "0fc0a36d-df90-57f3-8f93-d78a9fc72bb5"
version = "0.3.14"

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

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

[[deps.Missings]]
deps = ["DataAPI"]
git-tree-sha1 = "bf210ce90b6c9eed32d25dbcae1ebc565df2687f"
uuid = "e1d29d7a-bbdc-5cf2-9ac0-f12de2c33e28"
version = "1.0.2"

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

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

[[deps.NamedDims]]
deps = ["AbstractFFTs", "ChainRulesCore", "CovarianceEstimation", "LinearAlgebra", "Pkg", "Requires", "Statistics"]
git-tree-sha1 = "4f34e3ff2fa7f2a1c03fb2c4fdd637380d760bbc"
uuid = "356022a1-0364-5f58-8944-0da4b18d706f"
version = "0.2.42"

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

[[deps.OpenBLAS_jll]]
deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"

[[deps.OpenLibm_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "05823500-19ac-5b8b-9628-191a04bc5112"

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

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

[[deps.PDMats]]
deps = ["LinearAlgebra", "SparseArrays", "SuiteSparse"]
git-tree-sha1 = "ee26b350276c51697c9c2d88a072b339f9f03d73"
uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
version = "0.11.5"

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

[[deps.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"

[[deps.PlutoUI]]
deps = ["AbstractPlutoDingetjes", "Base64", "Dates", "Hyperscript", "HypertextLiteral", "IOCapture", "InteractiveUtils", "JSON", "Logging", "Markdown", "Random", "Reexport", "UUIDs"]
git-tree-sha1 = "565564f615ba8c4e4f40f5d29784aa50a8f7bbaf"
uuid = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
version = "0.7.22"

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

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

[[deps.QuadGK]]
deps = ["DataStructures", "LinearAlgebra"]
git-tree-sha1 = "78aadffb3efd2155af139781b8a8df1ef279ea39"
uuid = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
version = "2.4.2"

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

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

[[deps.RecipesBase]]
git-tree-sha1 = "6bf3f380ff52ce0832ddd3a2a7b9538ed1bcca7d"
uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
version = "1.2.1"

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

[[deps.Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

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

[[deps.Rmath]]
deps = ["Random", "Rmath_jll"]
git-tree-sha1 = "bf3188feca147ce108c76ad82c2792c57abe7b1f"
uuid = "79098fc4-a85e-5d69-aa6a-4863f24498fa"
version = "0.7.0"

[[deps.Rmath_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "68db32dff12bb6127bac73c209881191bf0efbb7"
uuid = "f50d1b31-88e8-58de-be2c-1cc44531875f"
version = "0.3.0+0"

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

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

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

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

[[deps.SimpleWeightedGraphs]]
deps = ["Graphs", "LinearAlgebra", "Markdown", "SparseArrays", "Test"]
git-tree-sha1 = "a6f404cc44d3d3b28c793ec0eb59af709d827e4e"
uuid = "47aef6b3-ad0c-573a-a1e2-d07658019622"
version = "1.2.1"

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

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

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

[[deps.SpecialFunctions]]
deps = ["ChainRulesCore", "IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
git-tree-sha1 = "f0bccf98e16759818ffc5d97ac3ebf87eb950150"
uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
version = "1.8.1"

[[deps.Static]]
deps = ["IfElse"]
git-tree-sha1 = "e7bc80dc93f50857a5d1e3c8121495852f407e6a"
uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
version = "0.4.0"

[[deps.StaticArrays]]
deps = ["LinearAlgebra", "Random", "Statistics"]
git-tree-sha1 = "3c76dde64d03699e074ac02eb2e8ba8254d428da"
uuid = "90137ffa-7385-5640-81b9-e52037218182"
version = "1.2.13"

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

[[deps.StatsAPI]]
git-tree-sha1 = "0f2aa8e32d511f758a2ce49208181f7733a0936a"
uuid = "82ae8749-77ed-4fe6-ae5f-f523153014b0"
version = "1.1.0"

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

[[deps.StatsFuns]]
deps = ["ChainRulesCore", "InverseFunctions", "IrrationalConstants", "LogExpFunctions", "Reexport", "Rmath", "SpecialFunctions"]
git-tree-sha1 = "bedb3e17cc1d94ce0e6e66d3afa47157978ba404"
uuid = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
version = "0.9.14"

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

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

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

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

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

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

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

[[deps.ZygoteRules]]
deps = ["MacroTools"]
git-tree-sha1 = "8c1a8e4dfacb1fd631745552c8db35d0deb09ea0"
uuid = "700de1a5-db45-46bc-99cf-38207098b444"
version = "0.2.2"

[[deps.libblastrampoline_jll]]
deps = ["Artifacts", "Libdl", "OpenBLAS_jll"]
uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

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

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

# ╔═╡ Cell order:
# ╟─e701b513-4ad9-474c-9d79-60a24f545ebd
# ╠═6f9be298-5374-11ec-3b9d-fd74a32d6b53
# ╠═087fb979-d68d-4ee7-9be5-fad992fdb792
# ╠═85d8985b-f5ee-4cba-bcc5-528049798330
# ╟─e62e36e9-19ce-4a7d-8bda-5634813c1ee9
# ╠═18d378dd-abc6-4fa6-a2c9-1cc09e4cec46
# ╟─dee4eb65-145c-4b2b-98a7-ae3a7045f24b
# ╠═2b485956-459d-40e7-9364-79ea77073804
# ╟─98c925ac-9ac8-42f3-9855-a01e5eaac350
# ╟─d22e825e-d8ff-44fe-8b06-b14c09765180
# ╠═157d13c8-7311-4832-a9ff-da22125b2bbf
# ╟─5af58241-50c7-464e-bfa3-53819e605793
# ╠═b372fcfc-cad0-48bc-b759-54d15bd2e28b
# ╟─62115e3f-cdfe-464c-91a1-eb80e674a610
# ╠═8653ca13-21e9-4049-a81e-02b6047d2c66
# ╠═05612a12-1080-4b79-b505-198e079a4b2a
# ╟─bdb85fb8-0707-47d1-95ee-9bf2e56500f5
# ╠═38b8f8af-ba81-4d3b-83b7-56b0d41ea708
# ╟─d6759a8e-a705-4d95-a8da-e172ada644d8
# ╟─7c6947a0-7d0b-4459-b999-58ea4f2b0169
# ╠═3f1d0989-2f97-4945-82f8-22e9c7bda6d2
# ╟─41ba8ac2-5a81-46f5-9269-404cc6a1d287
# ╟─e7fc05b3-a198-4a64-b62a-43cdf8f56238
# ╠═7b545a1e-b9c3-4c57-a5ee-3a629407b67a
# ╟─5157ac76-d359-4f18-883d-8afbe78d3b4d
# ╟─7640628b-e232-4ba3-94aa-d8a00218ff6a
# ╟─da970042-8326-4d17-b440-ab1960f10283
# ╠═7db3e138-6ec1-4a0c-8f3b-279f7251b4c3
# ╟─5e8e98f6-8896-462e-b610-4b379c727a72
# ╠═50854504-af68-4d7f-999e-f3054c088fb6
# ╟─4f35b3e8-90ca-4d21-a377-007f93eb0efd
# ╟─203a0538-d608-4d3e-b57d-edec20756300
# ╟─3e43f3fe-5c9a-4fb5-9033-db4a34bac5a5
# ╠═2f9ef252-5046-4af7-aa69-449a04bab26e
# ╠═6cfb38a8-c29f-42a6-b642-31c047dd6cb8
# ╠═db8110e0-d271-424c-b920-817dc1aab0c9
# ╟─0855f66e-f9af-4c63-b041-424c9f3fd655
# ╟─67e8e794-a16a-4d39-804a-6bfb350f8247
# ╟─446f7276-7ed8-4b5c-a3f1-fdde4fab1479
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
"""
    rasterKernelCircle(radius, raf::Array{Float64, 3})::Array{Float64, 3}

Expand single-pixel points in a raster (as obtained e.g. from `rasterize`) to
small circles with the specified `radius`. Alpha and color is preserved.
"""

function rasterKernelCircle(radius, raf::Raster)::Raster
    # premultiply the radius
    rt = radius * radius
    # maximum required kernel offset
    r = Int64(ceil(radius))
    # resulting raster
    res = zeros(Float64, size(raf))

    for x = 1:size(raf, 2)
        for y = 1:size(raf, 3)
            sr = 0.0
            sg = 0.0
            sb = 0.0
            pia = 1.0 # alpha is blended by multiplying (1-alpha) values
            w = 0.0 # total weight

            # run through the environment of x,y by offsets
            for ox = (-r):r
                for oy = (-r):r
                    # check if the offset is in limits
                    tx = x + ox
                    ty = y + oy
                    if tx < 1 || tx > size(raf, 2) || ty < 1 || ty > size(raf, 3)
                        continue
                    end
                    # cut out the circle (TODO: antialias the border??)
                    if ox * ox + oy * oy > rt
                        continue
                    end

                    # blend it
                    sr += raf[1, tx, ty] * raf[4, tx, ty]
                    sg += raf[2, tx, ty] * raf[4, tx, ty]
                    sb += raf[3, tx, ty] * raf[4, tx, ty]
                    w += raf[4, tx, ty]
                    pia *= 1 - raf[4, tx, ty]
                end
            end
            # compute the weighted mean & convert the alpha back
            if w > 0
                res[1, x, y] = sr / w
                res[2, x, y] = sg / w
                res[3, x, y] = sb / w
                res[4, x, y] = 1 - pia
            end
        end
    end
    res
end

"""
    rasterKernelGauss(sigma, raf::Raster; approximation_limit=3.5)::Raster

Similar to [`rasterKernelCircle`](@ref), but expands smooth gaussian
neighborhoods rather than circle-ish ones. `approximation_limit` sets the size
of the kernel square in multiplies of `sigma`.
"""
function rasterKernelGauss(sigma, raf::Raster; approximation_limit = 3.5)::Raster
    # premultiplied sigma
    st = -1.0 / (sigma * sigma)
    # limit
    r = Int64(ceil(approximation_limit * sigma))
    #result
    res = zeros(Float64, size(raf))

    for x = 1:size(raf, 2)
        for y = 1:size(raf, 3)
            sr = 0.0
            sg = 0.0
            sb = 0.0
            pia = 1.0 # alpha is blended by multiplying (1-alpha) values
            w = 0.0 # total weight

            # run through the environment of x,y by offsets
            for ox = (-r):r
                for oy = (-r):r
                    # check if the offset is in limits
                    tx = x + ox
                    ty = y + oy
                    if tx < 1 || tx > size(raf, 2) || ty < 1 || ty > size(raf, 3)
                        continue
                    end

                    # gaussian alpha
                    ga = raf[4, tx, ty] * exp(st * (ox * ox + oy * oy))

                    # blend it
                    sr += raf[1, tx, ty] * ga
                    sg += raf[2, tx, ty] * ga
                    sb += raf[3, tx, ty] * ga
                    w += ga
                    pia *= 1 - ga
                end
            end
            # compute the weighted mean & convert the alpha back
            if w > 0
                res[1, x, y] = sr / w
                res[2, x, y] = sg / w
                res[3, x, y] = sb / w
                res[4, x, y] = 1 - pia
            end
        end
    end
    res
end


text:
function mark(msg::String)
    active[] && ccall((:nvtxMarkA, libnvtx), Cvoid, (Cstring,), msg)
end


text:
macro addctx(t, func)
    n = length(func.args[2].args)
    insert!(func.args[2].args, n, :(pop!(ps.closer.cc)))
    pushfirst!(func.args[2].args, :(push!(ps.closer.cc, $(t))))
    return esc(func)
end

function accept_rparen(ps)
    if kindof(ps.nt) == Tokens.RPAREN
        return mPUNCTUATION(next(ps))
    else
        ps.errored = true
        return mErrorToken(mPUNCTUATION(Tokens.RPAREN, 0, 0), UnexpectedToken)
    end
end
accept_rparen(ps::ParseState, args) = push!(args, accept_rparen(ps))

function accept_rsquare(ps)
    if kindof(ps.nt) == Tokens.RSQUARE
        return mPUNCTUATION(next(ps))
    else
        ps.errored = true
        return mErrorToken(mPUNCTUATION(Tokens.RSQUARE, 0, 0), UnexpectedToken)
    end
end
accept_rsquare(ps::ParseState, args) = push!(args, accept_rsquare(ps))

function accept_rbrace(ps)
    if kindof(ps.nt) == Tokens.RBRACE
        return mPUNCTUATION(next(ps))
    else
        ps.errored = true
        return mErrorToken(mPUNCTUATION(Tokens.RBRACE, 0, 0), UnexpectedToken)
    end
end
accept_rbrace(ps::ParseState, args) = push!(args, accept_rbrace(ps))

function accept_end(ps::ParseState)
    if kindof(ps.nt) == Tokens.END
        return mKEYWORD(next(ps))
    else
        ps.errored = true
        return mErrorToken(mKEYWORD(Tokens.END, 0, 0), UnexpectedToken)
    end
end
accept_end(ps::ParseState, args) = push!(args, accept_end(ps))

function accept_comma(ps)
    if kindof(ps.nt) == Tokens.COMMA
        return mPUNCTUATION(next(ps))
    else
        return mPUNCTUATION(Tokens.RPAREN, 0, 0)
    end
end
accept_comma(ps::ParseState, args) = push!(args, accept_comma(ps))

function recover_endmarker(ps)
    if kindof(ps.nt) == Tokens.ENDMARKER
        if !isempty(ps.closer.cc)
            closert = last(ps.closer.cc)
            if closert == :block
                ps.errored = true
                return mErrorToken(mKEYWORD(Tokens.END, 0, 0), Unknown)
            elseif closert == :paren
                ps.errored = true
                return mErrorToken(mPUNCTUATION(Tokens.RPAREN, 0, 0), Unknown)
            elseif closert == :square
                ps.errored = true
                return mErrorToken(mPUNCTUATION(Tokens.RSQUARE, 0, 0), Unknown)
            elseif closert == :brace
                ps.errored = true
                return mErrorToken(mPUNCTUATION(Tokens.RBRACE, 0, 0), Unknown)
            end
        end
    end
end

function requires_ws(x, ps)
    if x.span == x.fullspan
        ps.errored = true
        return mErrorToken(x, Unknown)
    else
        return x
    end
end

function requires_no_ws(x, ps)
    if x.span != x.fullspan
        ps.errored = true
        return mErrorToken(x, UnexpectedWhiteSpace)
    else
        return x
    end
end


text:
module LRUCache

export LRU, @get!

include("list.jl")

# Default cache size
const __MAXCACHE__ = 100

type LRU{K,V} <: Associative{K,V}
    ht::Dict{K, LRUNode{K, V}}
    q::LRUList{K, V}
    maxsize::Int

    LRU(m::Int=__MAXCACHE__) = new(Dict{K, V}(), LRUList{K, V}(), m)
end
LRU(m::Int=__MAXCACHE__) = LRU{Any, Any}(m)

Base.show{K, V}(io::IO, lru::LRU{K, V}) = print(io,"LRU{$K, $V}($(lru.maxsize))")

Base.start(lru::LRU) = start(lru.ht)
Base.next(lru::LRU, state) = next(lru.ht, state)
Base.done(lru::LRU, state) = done(lru.ht, state)

Base.length(lru::LRU) = length(lru.q)
Base.isempty(lru::LRU) = isempty(lru.q)
Base.sizehint!(lru::LRU, n::Integer) = sizehint!(lru.ht, n)

Base.haskey(lru::LRU, key) = haskey(lru.ht, key)
Base.get(lru::LRU, key, default) = haskey(lru, key) ? lru[key] : default

macro get!(lru, key, default)
    quote
        if haskey($(esc(lru)), $(esc(key)))
            value = $(esc(lru))[$(esc(key))]
        else
            value = $(esc(default))
            $(esc(lru))[$(esc(key))] = value
        end
        value
    end
end

function Base.get!{K,V}(default::Base.Callable, lru::LRU{K, V}, key::K)
    if haskey(lru, key)
        return lru[key]
    else
        value = default()
        lru[key] = value
        return value
    end
end

function Base.get!{K,V}(lru::LRU{K,V}, key::K, default::V)
    if haskey(lru, key)
        return lru[key]
    else
        lru[key] = default
        return default
    end
end

function Base.getindex(lru::LRU, key)
    node = lru.ht[key]
    move_to_front!(lru.q, node)
    return node.v
end

function Base.setindex!{K, V}(lru::LRU{K, V}, v, key)
    if haskey(lru, key)
        item = lru.ht[key]
        item.v = v
        move_to_front!(lru.q, item)
    elseif length(lru) == lru.maxsize
        # At capacity. Roll the list so last el is now first, remove the old
        # data, and update new data in place.
        rotate!(lru.q)
        item = first(lru.q)
        delete!(lru.ht, item.k)
        item.k = key
        item.v = v
        lru.ht[key] = item
    else
        item = LRUNode{K, V}(key, v)
        unshift!(lru.q, item)
        lru.ht[key] = item
    end
    return lru
end

function Base.resize!(lru::LRU, n::Int)
    n < 0 && error("size must be a positive integer")
    lru.maxsize = n
    for i in 1:(length(lru) - lru.maxsize)
        rm = pop!(lru.q)
        delete!(lru.ht, rm.k)
    end
    return lru
end

function Base.delete!(lru::LRU, key)
    item = lru.ht[key]
    delete!(lru.q, item)
    delete!(lru.ht, key)
    return lru
end

function Base.empty!(lru::LRU)
    empty!(lru.ht)
    empty!(lru.q)
end

end # module


text:
module StaticNets

using Distributions
using StatsBase
using Optim
using LineSearches
using StatsFuns
using Roots
using MLBase
using Statistics
using LinearAlgebra
using Random
using Distributed

using ..AReg
using ..Utilities
using ..Scalings
using ..ErgmRcall

## STATIC NETWORK MODEL
abstract type Ergm end
abstract type ErgmDirBin <: Ergm end

#constants
targetErrValStaticNets = 1e-2
export targetErrValStaticNets

targetErrValStaticNetsW = 1e-5
export targetErrValStaticNetsW

bigConstVal = 10^6
export bigConstVal

maxLargeVal =  1e40# 1e-10 *sqrt(prevfloat(Inf))
export maxLargeVal

minSmallVal = 1e2*eps()
export minSmallVal

# region StaticNet interface (work in progress)

name(x::T where T <: Ergm) = error("to be defined")

type_of_obs(model::T where T <: Ergm) = error("to be defined")

obj_fun(model::T where T <: Ergm, obs_t, N, par) = error("to be defined")

grad_obj_fun(model::T where T <: Ergm, obs_t, N, par) = error("to be defined")

hessian_obj_fun(model::T where T <: Ergm, obs_t, N, par) = error("to be defined")

fisher_info_obj_fun(model::T where T <: Ergm, obs_t, N, par) = error("to be defined")

estimate(model::T where T <: Ergm, A) = error("to be defined")

estimate_sequence(model::T where T <: Ergm, obsT) = error("to be defined")

sample_ergm(model::T where T <: Ergm, N, parVec, nSample) = error("to be defined")  # return a 3D array N x N x nSample

"""
return a 4D array N x N x T x nSample
"""
function sample_ergm_sequence(model::T where T <: Ergm, N, parVecSeq_T::Matrix, nSample)

    T = size(parVecSeq_T)[2]
    A_T_dgp = zeros(Int8, N, N, T, nSample)
    for t=1:T
        A_T_dgp[:, :, t, :] = sample_ergm(model, get_N_t(N,t), parVecSeq_T[:, t], nSample)
    end

    return A_T_dgp
end

#endregion

include("./StaticNets_models/StaticNets_DirBin1.jl")
include("./StaticNets_models/StaticNets_DirBin0Rec0.jl")
include("./StaticNets_models/StaticNets_ErgmPML.jl")
include("./StaticNets_models/StaticNets_conf_intervals.jl")

end