text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.331138821460632 <= q <= 11.287379993959185)

                     add_NL_constraint(m, :(gelu(0.16893740747119512 + 0.361142465903193*gelu(0.3759145594542668 + 0.23452173305395663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.21914525995597334*gelu(-0.017859443401365915 + -0.7776638835685663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.7287835676612944*gelu(0.6615910141517634 + 0.18731902906517428*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + gelu(-0.9123133183860701 + 0.273513749523993*gelu(0.3759145594542668 + 0.23452173305395663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.052599029592083735*gelu(-0.017859443401365915 + -0.7776638835685663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.2642412570560455*gelu(0.6615910141517634 + 0.18731902906517428*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) + gelu(0.5514504616591211 + 0.27732859795053555*gelu(0.3759145594542668 + 0.23452173305395663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + 0.192500882174786*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.6554806276835787*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + -0.08702067567597149*gelu(-0.017859443401365915 + -0.7776638835685663*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.3390165730930925*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.16349638009140932*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3])))) + 0.40713382739201665*gelu(0.6615910141517634 + 0.18731902906517428*gelu(-0.8645963938709578 + -0.1467408759040536*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.7957413591360925*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + 0.6670298352855966*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.853933007586007*gelu(-0.9609056399980664 + 0.06060702618594416*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + -0.49474851412125664*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.5046120388520428*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))) + -0.9912964912419922*gelu(-0.9973666264873402 + 0.1126936245130743*gelu(0.5833709147085329 + -0.6819730710434229*$(x[1]) + 0.5963453169666697*$(x[2]) + -0.1905949830790421*$(x[3])) + 0.2338493567409592*gelu(0.2875838686662693 + -0.4272251449908575*$(x[1]) + 0.9704596097299736*$(x[2]) + -0.7281111839596055*$(x[3])) + -0.8532702233572196*gelu(0.428536919496763 + 0.9920670296914205*$(x[1]) + 0.8818625631698058*$(x[2]) + -0.43136616052165655*$(x[3]))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#######################################################################################
# Original correlation sum
#######################################################################################
using Distances, Roots
export kernelprob, correlationsum, grassberger

"""
    kernelprob(X, ε, norm = Euclidean()) → p
Associate each point in `X` (`Dataset` or timesries) with a probability `p` using the
"kernel estimation" (also called "nearest neighbor kernel estimation" and other names):
```math
p_j = \\frac{1}{N}\\sum_{i=1}^N I(||X_i - X_j|| < \\epsilon)
```
where ``N`` is its length and ``I`` gives 1 if the argument is `true`.
Because ``p`` is further normalized, it can be used as
an alternative for the [`genentropy`](@ref) function (using the second method).
"""
function kernelprob(X, ε, norm = Euclidean())
    N = length(X)
    p = zeros(eltype(X), N)
    @inbounds for i in 1:N
        p[i] = count(evaluate(norm, X[i], X[j]) < ε for j in 1:N)
    end
    p ./= sum(p)
    return p
end

"""
    correlationsum(X, ε::Real; w = 1, norm = Euclidean()) → C(ε)
Calculate the correlation sum of `X` (`Dataset` or timeseries)
for a given radius `ε` and `norm`, using the formula:
```math
C(\\epsilon) = \\frac{2}{(N-w)(N-w-1)}\\sum_{i=1}^{N}\\sum_{j=1+w+i}^{N} I(||X_i - X_j|| < \\epsilon)
```
where ``N`` is its length and ``I`` gives 1 if the argument is `true`.
`w` is the Theiler window, a correction to the correlation sum that skips points
that are temporally close with each other, with the aim of removing spurious correlations.

See the book "Nonlinear Time Series Analysis", Ch. 6, for a discussion
around `w` and choosing best values.

See [`grassberger`](@ref) for more.
See also [`takens_best_estimate`](@ref).
"""
function correlationsum(X, ε::Real; norm = Euclidean(), w = 1)
    N, C = length(X), 0
    @inbounds for i in 1:N
        C += count(evaluate(norm, X[i], X[j]) < ε for j in i+1+w:N)
    end
    return 2C/((N-w)*(N-1-w))
end

"""
    correlationsum(X, εs::AbstractVector; kwargs...) → Cs
Calculate the correlation sum for every `ε ∈ εs` using an optimized version.
"""
function correlationsum(X, εs::AbstractVector; norm = Euclidean(), w = 1)
    @assert issorted(εs) "Sorted εs required for optimized version."
    d = distancematrix(X, norm)
    Cs = zeros(length(εs))
    N = length(X)
    factor = 2/((N-w)*(N-1-w))
    for k in length(εs)÷2:-1:1
        ε = εs[k]
        for i in 1:N
            @inbounds Cs[k] += count(d[j, i] < ε for j in i+1+w:N)
        end
        Cs[k] == 0 && break
    end
    for k in (length(εs)÷2 + 1):length(εs)
        ε = εs[k]
        for i in 1:N
            @inbounds Cs[k] += count(d[j, i] < ε for j in i+1+w:N)
        end
        if Cs[k] ≈ 1/factor
            Cs[k:end] .= 1/factor
            break
        end
    end
    return Cs .* factor
end

function distancematrix(X, norm = Euclidean())
    N = length(X)
    d = zeros(eltype(X), N, N)
    @inbounds for i in 1:N
        for j in i+1:N
            d[j, i] = evaluate(norm, X[i], X[j])
        end
    end
    return d
end

"""
    grassberger(data, εs = estimate_boxsizes(data); kwargs...) → D_C
Use the method of Grassberger and Proccacia[^Grassberger1983], and the correction by
Theiler[^Theiler1986], to estimate the correlation dimension `D_C` of the given `data`.

This function does something extrely simple:
```julia
cm = correlationsum(data, εs; kwargs...)
return linear_region(log.(sizes), log(cm))[2]
```
i.e. it calculates [`correlationsum`](@ref) for various radii and then tries to find
a linear region in the plot of the log of the correlation sum versus log(ε).
See [`generalized_dim`](@ref) for a more thorough explanation.

See also [`takens_best_estimate`](@ref).

[^Grassberger1983]: Grassberger and Proccacia, [Characterization of strange attractors, PRL 50 (1983)](https://journals-aps-org.e-bis.mpimet.mpg.de/prl/abstract/10.1103/PhysRevLett.50.346)

[^Theiler1986]: Theiler, [Spurious dimension from correlation algorithms applied to limited time-series data. Physical Review A, 34](https://doi.org/10.1103/PhysRevA.34.2427)
"""
function grassberger(data::AbstractDataset, εs = estimate_boxsizes(data); kwargs...)
    cm = correlationsum(data, εs; kwargs...)
    return linear_region(log.(εs), log.(cm))[2]
end

#######################################################################################
# Takens' best estimate
#######################################################################################
export takens_best_estimate

"""
    takens_best_estimate(X, εmax, metric = Chebyshev(),εmin = 0) → D_C, D_C_95u, D_C_95l
Use the so-called "Takens' best estimate" [^Takens1985][^Theiler1988]
method for estimating the correlation dimension
`D_C` and the upper (`D_C_95u`) and lower (`D_C_95l`) confidence limit for the given dataset `X`.

The original formula is
```math
D_C \\approx \\frac{C(\\epsilon_\\text{max})}{\\int_0^{\\epsilon_\\text{max}}(C(\\epsilon) / \\epsilon) \\, d\\epsilon}
```
where ``C`` is the [`correlationsum`](@ref) and ``\\epsilon_\\text{max}`` is an upper cutoff.
Here we use the later expression
```math
D_C \\approx - \\frac{1}{\\eta},\\quad \\eta = \\frac{1}{(N-1)^*}\\sum_{[i, j]^*}\\log(||X_i - X_j|| / \\epsilon_\\text{max})
```
where the sum happens for all ``i, j`` so that ``i < j`` and ``||X_i - X_j|| < \\epsilon_\\text{max}``. In the above expression, the bias in the original paper has already been corrected, as suggested in [^Borovkova1999].

The confidence limits are estimated from the log-likelihood function by finding
the values of `D_C` where the function has fallen by 2 from its maximum, see e.g.
[^Barlow] chapter 5.3
Because the CLT does not apply (no independent measurements), the limits are not
neccesarily symmetric.

According to [^Borovkova1999], introducing a lower cutoff `εmin` can make the
algorithm more stable (no divergence), this option is given but defaults to zero.

If `X` comes from a delay coordinates embedding of a timseries `x`, a recommended value
for ``\\epsilon_\\text{max}`` is `std(x)/4`.

[^Takens1985]: Takens, On the numerical determination of the dimension of an attractor, in: B.H.W. Braaksma, B.L.J.F. Takens (Eds.), Dynamical Systems and Bifurcations, in: Lecture Notes in Mathematics, Springer, Berlin, 1985, pp. 99–106.
[^Theiler1988]: Theiler, [Lacunarity in a best estimator of fractal dimension. Physics Letters A, 133(4–5)](https://doi.org/10.1016/0375-9601(88)91016-X)
[^Borovkova1999]: Borovkova et al., [Consistency of the Takens estimator for the correlation dimension. The Annals of Applied Probability, 9, 05 1999.](https://doi.org/10.1214/aoap/1029962747)
[^Barlow]: Barlow, R., Statistics - A Guide to the Use of Statistical Methods in the Physical Sciences. Vol 29. John Wiley & Sons, 1993
"""
function takens_best_estimate(X, εmax, metric = Chebyshev(); εmin=0)
    n, η, N = 0, zero(eltype(X)), length(X)
    @inbounds for i in 1:N
        for j in i+1:N
            d = evaluate(metric, X[i], X[j])
            if εmin < d < εmax
                n += 1
                η += log(d/εmax)
            end
        end
    end
    # bias-corrected version (log-likelihood function shifted on x-axis)
    α = -(n-1)/η
    # biased version (maximum of original log-likelihood function)
    α_b = -n/η
    # value of maximum of original log-likelihood function
    mxl = n*log(α_b) + α_b * η
    # at the 95%-confidence interval, the log-l function has dropped by 2
    # -> log_l(x) - mxl + 2 = 0
    # this is a result of the invariance of the MLE, a really nice property
    # these limits are not going to be perfectly symmetric (CLT does not apply)
    mn, mx = fzeros(x-> n * log(x) + η * x - mxl +2 , 0,2*α)

    # Since the bias-correction is just a shift of the log-l function on the
    # x-axis, we can easily shift the confidence limits by the bias α-α_b
    α95u = α - α_b + mn
    α95l = α - α_b + mx

    return α, α95u, α95l
end

ext: jl


text:
export AbstractNode, NullNode, ListNode, LinkedList

abstract type AbstractNode end

struct NullNode <: AbstractNode end

value(node::NullNode) = nothing

next(node::NullNode) = NullNode()

prev(node::NullNode) = NullNode()

Base.show(io::IO, node::NullNode) = nothing

mutable struct ListNode{T} <: AbstractNode
    data::T
    prev::AbstractNode
    next::AbstractNode

    ListNode{T}(data::T) where T = new(data, NullNode(), NullNode())
end

value(node::ListNode) = node.data

next(node::ListNode) = node.next

prev(node::ListNode) = node.prev

Base.show(io::IO, node::ListNode) = print(io, value(node))

mutable struct LinkedList{T}
    head_node::AbstractNode
    length::Int64
end

LinkedList{T}() where T = LinkedList{T}(NullNode(), 0)

Base.length(ll::LinkedList) = ll.length

is_empty(ll::LinkedList) = (ll.length < 1)

head(ll::LinkedList) = ll.head_node

function Base.show(io::IO, ll::LinkedList{T}) where T
    print(io, "LinkedList{$T}([")

    node = head(ll)
    while !(next(node) isa NullNode)
        print(io, "$node, ")
        node = next(node)
    end

    print(io, "$(node)])")
end

function _getindex(ll::LinkedList, n::Int64)
    if is_empty(ll) throw(BoundsError("List is empty.")) end
    if !(0 < n <= ll.length) throw(BoundsError("Out of boundary.")) end

    node = head(ll)
    for index=2:n
        node = next(node)
    end

    return node
end

Base.getindex(ll::LinkedList, i::Int64) = value(_getindex(ll, i))

function Base.setindex!(ll::LinkedList{T}, data::T, i::Int64) where T
    _getindex(ll, i).data = data
end

function Base.pushfirst!(ll::LinkedList, data::T) where T
    node = ListNode{T}(data)
    head_node = head(ll)
    ll.head_node = node

    if !is_empty(ll)
        head_node.prev = node
        node.next = head_node
    end

    ll.length += 1
end

function Base.insert!(ll::LinkedList, data::T, i::Int64) where T
    if i < 1
        throw(BoundsError("Out of boundary."))
    elseif i < 2
        Base.pushfirst!(ll, data)
    else
        node = ListNode{T}(data)
        current_node = _getindex(ll, i-1)
        next_node = next(current_node)

        node.prev = current_node
        current_node.next = node
        node.next = next_node
        !(next_node isa NullNode) && (next_node.prev = node)

        ll.length += 1
    end
end

function Base.push!(ll::LinkedList, data::T) where T
    is_empty(ll) ? Base.pushfirst!(ll, data) : Base.insert!(ll, data, ll.length+1)
end

function Base.popfirst!(ll::LinkedList)
    is_empty(ll) && throw(BoundsError("List is empty."))

    head_node = head(ll)
    if !(head_node.next isa NullNode) head_node.next.prev = NullNode() end

    ll.head_node = next(head_node)
    ll.length -= 1

    return value(head_node)
end

function Base.delete!(ll::LinkedList, i::Int64)
    if i < 1
        throw(BoundsError("Out of boundary."))
    elseif i < 2
        return Base.popfirst!(ll)
    else
        current_node = _getindex(ll, i)
        prev_node = prev(current_node)
        next_node = next(current_node)

        !(next_node isa NullNode) && (next_node.prev = prev_node)
        !(prev_node isa NullNode) && (prev_node.next = next_node)

        ll.length -= 1

        return value(current_node)
    end
end

Base.pop!(ll::LinkedList) = Base.delete!(ll, ll.length)


text:
module TrackedTests

using ReverseDiff, Test
using ReverseDiff: TrackedReal, TrackedArray

include(joinpath(dirname(@__FILE__), "utils.jl"))

samefields(a, b) = a === b

function samefields(a::TrackedReal, b::TrackedReal)
    same_origin = false
    if isdefined(a, :origin) === isdefined(b, :origin)
        if isdefined(a, :origin)
            same_origin = a.origin === b.origin
        else
            same_origin = true
        end
    end
    return (ReverseDiff.value(a) == ReverseDiff.value(b) &&
            ReverseDiff.deriv(a) == ReverseDiff.deriv(b) &&
            ReverseDiff.tape(a) === ReverseDiff.tape(b) &&
            a.index === b.index && same_origin)
end

function samefields(a::TrackedArray, b::TrackedArray)
    return (ReverseDiff.value(a) == ReverseDiff.value(b) &&
            ReverseDiff.deriv(a) == ReverseDiff.deriv(b) &&
            ReverseDiff.tape(a) === ReverseDiff.tape(b))
end

################
# Constructors #
################

# TrackedReal #
#-------------#

v = rand()
d = rand(Int)
o = [v]
tp = InstructionTape()
t = TrackedReal(v, d, tp, 1, o)

@test t.value === v
@test t.deriv === d
@test t.tape === tp
@test t.index === 1
@test t.origin === o
@test isa(t, TrackedReal{Float64,Int,Vector{Float64}})

t = TrackedReal(v, d, tp)

@test t.value === v
@test t.deriv === d
@test t.tape === tp
@test t.index === ReverseDiff.NULL_INDEX
@test t.origin === nothing
@test isa(t, TrackedReal{Float64,Int,Nothing})

t = TrackedReal(v, d)

@test t.value === v
@test t.deriv === d
@test t.tape === ReverseDiff.NULL_TAPE
@test t.index === ReverseDiff.NULL_INDEX
@test t.origin === nothing
@test isa(t, TrackedReal{Float64,Int,Nothing})

# TrackedArray #
#--------------#

v = rand(3, 2, 1)
d = rand(Int, 3, 2, 1)
tp = InstructionTape()

t = TrackedArray(v, d, tp)

@test t.value === v
@test t.deriv === d
@test t.tape === tp
@test isa(t, TrackedArray{Float64,Int,3,Array{Float64,3},Array{Int,3}})

###########
# getters #
###########

# istracked #
#-----------#

@test !(ReverseDiff.istracked(nothing))
@test !(ReverseDiff.istracked(rand()))
@test !(ReverseDiff.istracked([1]))
@test ReverseDiff.istracked(TrackedArray(rand(1), rand(1), InstructionTape()))
@test ReverseDiff.istracked(TrackedReal(1, 1))
@test ReverseDiff.istracked(Any[1])
@test ReverseDiff.istracked([TrackedReal(1, 1)])

# value #
#-------#

v, d = rand(), rand()
varr, darr = [v], [d]
any_varr = Any[v]

@test ReverseDiff.value(nothing) === nothing
@test ReverseDiff.value(v) === v
@test ReverseDiff.value(TrackedArray(varr, darr, InstructionTape())) === varr
@test ReverseDiff.value(TrackedReal(v, d)) === v
@test ReverseDiff.value(varr) === varr
@test ReverseDiff.value(any_varr) !== any_varr
@test ReverseDiff.value(any_varr) == any_varr
@test ReverseDiff.value([TrackedReal(v, d)]) == varr

# deriv #
#-------#

v, d = rand(), rand()
varr, darr = [v], [d]
any_varr = Any[v]

@test ReverseDiff.deriv(TrackedArray(varr, darr, InstructionTape())) === darr
@test ReverseDiff.deriv(TrackedReal(v, d)) === d

# valtype #
#---------#

@test ReverseDiff.valtype(TrackedReal(1,0)) === Int
@test ReverseDiff.valtype(TrackedReal(1.0, 0)) === Float64
@test ReverseDiff.valtype(typeof(TrackedReal(1,0))) === Int
@test ReverseDiff.valtype(typeof(TrackedReal(1.0, 0))) === Float64
@test ReverseDiff.valtype(TrackedArray([1], [0], InstructionTape())) === Int
@test ReverseDiff.valtype(TrackedArray([1.0], [0], InstructionTape())) === Float64
@test ReverseDiff.valtype(typeof(TrackedArray([1], [0], InstructionTape()))) === Int
@test ReverseDiff.valtype(typeof(TrackedArray([1.0], [0], InstructionTape()))) === Float64

# derivtype #
#-----------#

@test ReverseDiff.derivtype(TrackedReal(1,0)) === Int
@test ReverseDiff.derivtype(TrackedReal(1, 0.0)) === Float64
@test ReverseDiff.derivtype(typeof(TrackedReal(1,0))) === Int
@test ReverseDiff.derivtype(typeof(TrackedReal(1, 0.0))) === Float64
@test ReverseDiff.derivtype(TrackedArray([1], [0], InstructionTape())) === Int
@test ReverseDiff.derivtype(TrackedArray([1], [0.0], InstructionTape())) === Float64
@test ReverseDiff.derivtype(typeof(TrackedArray([1], [0], InstructionTape()))) === Int
@test ReverseDiff.derivtype(typeof(TrackedArray([1], [0.0], InstructionTape()))) === Float64

# hasorigin #
#-----------#

@test !(ReverseDiff.hasorigin(rand()))
@test !(ReverseDiff.hasorigin(TrackedReal(rand(), rand())))
@test ReverseDiff.hasorigin(TrackedReal(rand(), rand(), InstructionTape(), 1, [rand()]))

# tape/hastape #
#--------------#

tp = InstructionTape()

null_tape_items = (nothing, rand(), [1], Any[1], [TrackedReal(1, 1)], TrackedReal(1, 1),
                   TrackedArray(rand(1), rand(1), ReverseDiff.NULL_TAPE))
tape_items = ([3, TrackedReal(1, 1, tp)], TrackedReal(1, 1, tp),
              TrackedArray(rand(1), rand(1), tp))

for i in null_tape_items
    @test ReverseDiff.tape(i) === ReverseDiff.NULL_TAPE
    @test !(ReverseDiff.hastape(i))
end

for i in tape_items
    @test ReverseDiff.tape(i) === tp
    @test ReverseDiff.hastape(i)
    for j in null_tape_items
        @test ReverseDiff.tape(i, j) === tp
        @test ReverseDiff.tape(j, i) === tp
        @test ReverseDiff.tape(i, j, i) === tp
        @test ReverseDiff.tape(j, i, i) === tp
        @test ReverseDiff.tape(i, i, j) === tp
    end
    for j in tape_items
        @test ReverseDiff.tape(i, j) === tp
        @test ReverseDiff.tape(i, j, i) === tp
    end
end

###########
# setters #
###########

v, d = rand(), rand()
v2, d2 = rand(), rand()
varr, darr = rand(3), rand(3)
varr2, darr2 = rand(3), rand(3)

tr = TrackedReal(v, d)
ta = TrackedArray(varr, darr, InstructionTape())

# value! #
#--------#

ReverseDiff.value!(tr, v2)
@test ReverseDiff.value(tr) === v2

ReverseDiff.value!(ta, varr2)
@test ReverseDiff.value(ta) == varr2
@test ReverseDiff.value(ta) === varr

ReverseDiff.value!(tr, v2)
@test ReverseDiff.value(tr) === v2

varr_copy = copy(varr)
ReverseDiff.value!((ta, tr), (varr, v))
@test ReverseDiff.value(tr) === v
@test ReverseDiff.value(ta) == varr_copy
@test ReverseDiff.value(ta) === varr

# deriv! #
#--------#

ReverseDiff.deriv!(tr, d2)
@test ReverseDiff.deriv(tr) === d2

ReverseDiff.deriv!(ta, darr2)
@test ReverseDiff.deriv(ta) == darr2
@test ReverseDiff.deriv(ta) === darr

ReverseDiff.deriv!(tr, d2)
@test ReverseDiff.deriv(tr) === d2

darr_copy = copy(darr)
ReverseDiff.deriv!((ta, tr), (darr, d))
@test ReverseDiff.deriv(tr) === d
@test ReverseDiff.deriv(ta) == darr_copy
@test ReverseDiff.deriv(ta) === darr

# pulling values from origin #
#----------------------------#

ta = TrackedArray(rand(3), rand(3), InstructionTape())
ta1 = ta[1]
v_old = ReverseDiff.value(ta)[1]
v_new = rand()
tr = TrackedReal(v_old, v_old)
varr_copy = copy(ReverseDiff.value(ta))
trs = Any[tr, ta[1], ta[2]]

@test ReverseDiff.pull_value!(nothing) === nothing
@test ReverseDiff.pull_value!([1,2,3]) === nothing
@test ReverseDiff.pull_value!(1.0) === nothing

ReverseDiff.value(ta)[1] = v_new
@test ReverseDiff.value(ta1) === v_old
ReverseDiff.pull_value!(ta1)
@test ReverseDiff.value(ta1) === v_new
ReverseDiff.value(ta)[1] = v_old
ReverseDiff.pull_value!(ta1)

@test ReverseDiff.value(tr) === v_old
ReverseDiff.pull_value!(tr)
@test ReverseDiff.value(tr) === v_old

ReverseDiff.pull_value!(ta)
@test ReverseDiff.value(ta) == varr_copy

ReverseDiff.value(ta)[1] = v_new + 1
ReverseDiff.value(ta)[2] = v_new + 2

@test ReverseDiff.value(trs[1]) === ReverseDiff.value(tr)
@test ReverseDiff.value(trs[2]) === varr_copy[1]
@test ReverseDiff.value(trs[3]) === varr_copy[2]

ReverseDiff.pull_value!(trs)

@test ReverseDiff.value(trs[1]) === ReverseDiff.value(tr)
@test ReverseDiff.value(trs[2]) === ReverseDiff.value(ta)[1]
@test ReverseDiff.value(trs[3]) === ReverseDiff.value(ta)[2]

# pulling derivs from origin #
#----------------------------#

ta = TrackedArray(rand(3), rand(3), InstructionTape())
ta1 = ta[1]
d_old = ReverseDiff.deriv(ta)[1]
d_new = rand()
tr = TrackedReal(d_old, d_old)
darr_copy = copy(ReverseDiff.deriv(ta))
trs = Any[tr, ta[1], ta[2]]

@test ReverseDiff.pull_deriv!(nothing) === nothing
@test ReverseDiff.pull_deriv!([1,2,3]) === nothing
@test ReverseDiff.pull_deriv!(1.0) === nothing

ReverseDiff.deriv(ta)[1] = d_new
@test ReverseDiff.deriv(ta1) === d_old
ReverseDiff.pull_deriv!(ta1)
@test ReverseDiff.deriv(ta1) === d_new
ReverseDiff.deriv(ta)[1] = d_old
ReverseDiff.pull_deriv!(ta1)

@test ReverseDiff.deriv(tr) === d_old
ReverseDiff.pull_deriv!(tr)
@test ReverseDiff.deriv(tr) === d_old

ReverseDiff.pull_deriv!(ta)
@test ReverseDiff.deriv(ta) == darr_copy

ReverseDiff.deriv(ta)[1] = d_new + 1
ReverseDiff.deriv(ta)[2] = d_new + 2

@test ReverseDiff.deriv(trs[1]) === ReverseDiff.deriv(tr)
@test ReverseDiff.deriv(trs[2]) === darr_copy[1]
@test ReverseDiff.deriv(trs[3]) === darr_copy[2]

ReverseDiff.pull_deriv!(trs)

@test ReverseDiff.deriv(trs[1]) === ReverseDiff.deriv(tr)
@test ReverseDiff.deriv(trs[2]) === ReverseDiff.deriv(ta)[1]
@test ReverseDiff.deriv(trs[3]) === ReverseDiff.deriv(ta)[2]

# push derivs from origin #
#-------------------------#

ta = TrackedArray(rand(3), rand(3), InstructionTape())
ta1 = ta[1]
d_old = ReverseDiff.deriv(ta)[1]
d_new = rand()
tr = TrackedReal(d_old, d_old)
darr_copy = copy(ReverseDiff.deriv(ta))
trs = Any[tr, ta[1], ta[2]]

@test ReverseDiff.push_deriv!(nothing) === nothing
@test ReverseDiff.push_deriv!([1,2,3]) === nothing
@test ReverseDiff.push_deriv!(1.0) === nothing

ReverseDiff.deriv!(ta1, d_new)
@test ReverseDiff.deriv(ta)[1] === d_old
ReverseDiff.push_deriv!(ta1)
@test ReverseDiff.deriv(ta)[1] === d_new
ReverseDiff.deriv(ta)[1] = d_old
ReverseDiff.pull_deriv!(ta1)

@test ReverseDiff.deriv(tr) === d_old
ReverseDiff.push_deriv!(tr)
@test ReverseDiff.deriv(tr) === d_old

ReverseDiff.push_deriv!(ta)
@test ReverseDiff.deriv(ta) == darr_copy
@test ReverseDiff.deriv(ta1) === d_old

ReverseDiff.deriv!(trs[2], d_new + 1)
ReverseDiff.deriv!(trs[3], d_new + 2)

@test ReverseDiff.deriv(ta)[1] === darr_copy[1]
@test ReverseDiff.deriv(ta)[2] === darr_copy[2]

ReverseDiff.push_deriv!(trs)

@test ReverseDiff.deriv(ta)[1] === ReverseDiff.deriv(trs[2])
@test ReverseDiff.deriv(ta)[2] === ReverseDiff.deriv(trs[3])

# seed!/unseed! #
#---------------#

v, d = rand(), rand()
varr, darr = rand(3), rand(3)
varr_copy, darr_copy = copy(varr), copy(darr)
ta = TrackedArray(varr, darr, InstructionTape())
tr = TrackedReal(v, d)
trs = Any[tr, ta[1], ta[2]]

@test ReverseDiff.seed!(nothing) === nothing
@test ReverseDiff.seed!([1,2,3]) === nothing
@test ReverseDiff.seed!(1.0) === nothing
@test ReverseDiff.unseed!(nothing) === nothing
@test ReverseDiff.unseed!([1,2,3]) === nothing
@test ReverseDiff.unseed!(1.0) === nothing

ReverseDiff.seed!(ta)
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy
ReverseDiff.unseed!(ta)
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == fill!(similar(darr), 0)
ReverseDiff.deriv!(ta, darr_copy)
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy

ReverseDiff.seed!(tr)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === 1.0
ReverseDiff.unseed!(tr)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === 0.0
ReverseDiff.deriv!(tr, d)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d

ReverseDiff.seed!(trs)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy
@test ReverseDiff.value(trs[2]) === varr[1]
@test ReverseDiff.deriv(trs[2]) === darr[1]
@test ReverseDiff.value(trs[3]) === varr[2]
@test ReverseDiff.deriv(trs[3]) === darr[2]

ReverseDiff.unseed!(trs)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === 0.0
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == [0.0, 0.0, darr_copy[3]]
@test ReverseDiff.value(trs[2]) === varr[1]
@test ReverseDiff.deriv(trs[2]) === 0.0
@test ReverseDiff.value(trs[3]) === varr[2]
@test ReverseDiff.deriv(trs[3]) === 0.0

ReverseDiff.deriv!(ta, darr_copy)
ReverseDiff.deriv!(tr, d)
ReverseDiff.pull_deriv!(trs)
@test trs[1] === tr
@test ReverseDiff.deriv(tr) === d
@test ReverseDiff.deriv(ta) === darr == darr_copy
@test ReverseDiff.deriv(trs[2]) === darr[1]
@test ReverseDiff.deriv(trs[3]) === darr[2]

ReverseDiff.seed!((tr, ta, trs))
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy
@test ReverseDiff.value(trs[2]) === varr[1]
@test ReverseDiff.deriv(trs[2]) === darr[1]
@test ReverseDiff.value(trs[3]) === varr[2]
@test ReverseDiff.deriv(trs[3]) === darr[2]

ReverseDiff.unseed!((tr, ta, trs))
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === 0.0
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == fill!(similar(darr), 0)
@test ReverseDiff.value(trs[2]) === varr[1]
@test ReverseDiff.deriv(trs[2]) === 0.0
@test ReverseDiff.value(trs[3]) === varr[2]
@test ReverseDiff.deriv(trs[3]) === 0.0

# increment deriv #
#-----------------#

x = rand(3)
v, d = rand(), rand()
varr, darr = rand(3), rand(3)
varr_copy, darr_copy = copy(varr), copy(darr)
ta = TrackedArray(varr, darr, InstructionTape())
tr = TrackedReal(v, d)
trs = Any[tr, ta[1], ta[2]]

ReverseDiff.increment_deriv!(tr, x[1])
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d + x[1]
ReverseDiff.deriv!(tr, d)

ta1 = ta[1]
ReverseDiff.deriv!(ta, x)
ReverseDiff.increment_deriv!(ta1, d)
@test ReverseDiff.value(ta1) === ReverseDiff.value(ta)[1] === varr_copy[1]
@test ReverseDiff.deriv(ta1) === ReverseDiff.deriv(ta)[1] === x[1] + d
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.increment_deriv!(ta, x)
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy + x
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.increment_deriv!(ta, x[1])
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy .+ x[1]
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.increment_deriv!(trs, x)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d + x[1]
@test ReverseDiff.value(trs[2]) === varr[1] === varr_copy[1]
@test ReverseDiff.deriv(trs[2]) === darr[1] === darr_copy[1] + x[2]
@test ReverseDiff.value(trs[3]) === varr[2] === varr_copy[2]
@test ReverseDiff.deriv(trs[3]) === darr[2] === darr_copy[2] + x[3]
ReverseDiff.deriv!(tr, d)
ReverseDiff.deriv!(ta, darr_copy)
ReverseDiff.pull_deriv!(trs)

ReverseDiff.increment_deriv!(trs, x[1])
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d + x[1]
@test ReverseDiff.value(trs[2]) === varr[1] === varr_copy[1]
@test ReverseDiff.deriv(trs[2]) === darr[1] === darr_copy[1] + x[1]
@test ReverseDiff.value(trs[3]) === varr[2] === varr_copy[2]
@test ReverseDiff.deriv(trs[3]) === darr[2] === darr_copy[2] + x[1]
ReverseDiff.deriv!(tr, d)
ReverseDiff.deriv!(ta, darr_copy)
ReverseDiff.pull_deriv!(trs)

# decrement deriv #
#-----------------#

x = rand(3)
v, d = rand(), rand()
varr, darr = rand(3), rand(3)
varr_copy, darr_copy = copy(varr), copy(darr)
ta = TrackedArray(varr, darr, InstructionTape())
tr = TrackedReal(v, d)
trs = Any[tr, ta[1], ta[2]]

ReverseDiff.decrement_deriv!(tr, x[1])
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d - x[1]
ReverseDiff.deriv!(tr, d)

ta1 = ta[1]
ReverseDiff.deriv!(ta, x)
ReverseDiff.decrement_deriv!(ta1, d)
@test ReverseDiff.value(ta1) === ReverseDiff.value(ta)[1] === varr_copy[1]
@test ReverseDiff.deriv(ta1) === ReverseDiff.deriv(ta)[1] === x[1] - d
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.decrement_deriv!(ta, x)
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy - x
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.decrement_deriv!(ta, x[1])
@test ReverseDiff.value(ta) === varr == varr_copy
@test ReverseDiff.deriv(ta) === darr == darr_copy .- x[1]
ReverseDiff.deriv!(ta, darr_copy)

ReverseDiff.decrement_deriv!(trs, x)
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d - x[1]
@test ReverseDiff.value(trs[2]) === varr[1] === varr_copy[1]
@test ReverseDiff.deriv(trs[2]) === darr[1] === darr_copy[1] - x[2]
@test ReverseDiff.value(trs[3]) === varr[2] === varr_copy[2]
@test ReverseDiff.deriv(trs[3]) === darr[2] === darr_copy[2] - x[3]
ReverseDiff.deriv!(tr, d)
ReverseDiff.deriv!(ta, darr_copy)
ReverseDiff.pull_deriv!(trs)

ReverseDiff.decrement_deriv!(trs, x[1])
@test ReverseDiff.value(tr) === v
@test ReverseDiff.deriv(tr) === d - x[1]
@test ReverseDiff.value(trs[2]) === varr[1] === varr_copy[1]
@test ReverseDiff.deriv(trs[2]) === darr[1] === darr_copy[1] - x[1]
@test ReverseDiff.value(trs[3]) === varr[2] === varr_copy[2]
@test ReverseDiff.deriv(trs[3]) === darr[2] === darr_copy[2] - x[1]
ReverseDiff.deriv!(tr, d)
ReverseDiff.deriv!(ta, darr_copy)
ReverseDiff.pull_deriv!(trs)

###########
# capture #
###########

v, d = rand(), rand()
varr, darr = rand(3), rand(3)
ta = TrackedArray(varr, darr, InstructionTape())
tr = TrackedReal(v, d)
trs = Any[tr, ta[1], ta[2]]

@test ReverseDiff.capture(ta) === ta
@test ReverseDiff.capture(tr) === v
@test samefields(ReverseDiff.capture(ta[1]), ta[1])
@test all(map(samefields, ReverseDiff.capture(trs), Any[v, ta[1], ta[2]]))

########################
# Conversion/Promotion #
########################

v, d = rand(), rand()
varr, darr = rand(3), rand(3)
tp = InstructionTape()
ta = TrackedArray(varr, darr, tp)
tr = TrackedReal(v, d)
ta1 = ta[1]
A = TrackedArray{BigInt,Float64,1,Array{BigInt,1},Array{Float64,1}}
T = TrackedReal{BigInt,Float64,A}

t2 = convert(TrackedReal{BigFloat,BigFloat,Nothing}, ta1)
@test length(tp) == 1
instr = tp[1]
@test instr.func === convert
@test instr.input === ta1
@test samefields(instr.output, TrackedReal(big(varr[1]), big(darr[1]), tp))
@test instr.cache === nothing
empty!(tp)

@test convert(Float64, tr) === v
@test convert(BigFloat, tr) == v

@test samefields(convert(T, 1), T(big(1), 0.0))
@test samefields(convert(TrackedReal{BigInt,Float64,Nothing}, 1), TrackedReal(big(1), 0.0))

@test convert(typeof(tr), tr) === tr
@test convert(typeof(ta), ta) === ta
@test convert(typeof(ta1), ta1) === ta1

@test promote_type(T, BigFloat) === TrackedReal{BigFloat,Float64,A}
@test promote_type(T, Float64) === TrackedReal{BigFloat,Float64,A}
@test promote_type(T, TrackedReal{BigFloat,BigFloat,Nothing}) === TrackedReal{BigFloat,BigFloat,Nothing}
@test promote_type(T, T) === T

###########################
# AbstractArray Interface #
###########################

varr, darr = rand(3, 3), rand(3, 3)
tp = InstructionTape()
ta = TrackedArray(varr, darr, tp)

@test isa(similar(ta), Matrix{eltype(ta)})

@test samefields(ta[2], TrackedReal(varr[2], darr[2], tp, 2, ta))

ta_sub = ta[:,:]
idx = ReverseDiff.index_iterable(axes(ta), (:, :))
@test collect(idx) == [(i, j) for i in 1:3, j in 1:3]
@test samefields(ta_sub, ta)
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, ta)
@test instr.cache === nothing
empty!(tp)

ta_sub = ta[:,1:2]
idx = ReverseDiff.index_iterable(axes(ta), (:, 1:2))
@test collect(idx) == [(i, j) for i in 1:3, j in 1:2]
@test samefields(ta_sub, TrackedArray(varr[:,1:2], darr[:,1:2], tp))
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, TrackedArray(varr[:,1:2], darr[:,1:2], tp))
@test instr.cache === nothing
empty!(tp)

ta_sub = ta[2:3,:]
idx = ReverseDiff.index_iterable(axes(ta), (2:3, :))
@test collect(idx) == [(i, j) for i in 2:3, j in 1:3]
@test samefields(ta_sub, TrackedArray(varr[2:3,:], darr[2:3,:], tp))
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, TrackedArray(varr[2:3,:], darr[2:3,:], tp))
@test instr.cache === nothing
empty!(tp)

ta_sub = ta[1:2,2:3]
idx = ReverseDiff.index_iterable(axes(ta), (1:2, 2:3))
@test collect(idx) == [(i, j) for i in 1:2, j in 2:3]
@test samefields(ta_sub, TrackedArray(varr[1:2,2:3], darr[1:2,2:3], tp))
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, TrackedArray(varr[1:2,2:3], darr[1:2,2:3], tp))
@test instr.cache === nothing
empty!(tp)

ta_sub = ta[2:6]
idx = ReverseDiff.index_iterable(axes(ta), (2:6,))
@test collect(idx) == [(i,) for i in 2:6]
@test samefields(ta_sub, TrackedArray(varr[2:6], darr[2:6], tp))
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, TrackedArray(varr[2:6], darr[2:6], tp))
@test instr.cache === nothing
empty!(tp)

ta_sub = ta[:]
idx = ReverseDiff.index_iterable(axes(ta), (:,))
@test collect(idx) == [(i, j) for i in 1:3, j in 1:3]
@test samefields(ta_sub, TrackedArray(varr[:], darr[:], tp))
@test length(tp) == 1
instr = tp[1]
@test instr.func === getindex
@test instr.input === (ta, idx)
@test samefields(instr.output, TrackedArray(varr[:], darr[:], tp))
@test instr.cache === nothing
empty!(tp)

@test IndexStyle(ta) === IndexLinear()

@test Base.size(ta) === size(varr)

@test Base.copy(ta) === ta

@test all(samefields.(ta, copyto!(similar(ta), ta)))

####################
# `Real` Interface #
####################

v_int, v_float, d = rand(Int), rand(), rand()
tp = InstructionTape()
tr_int = TrackedReal(v_int, d, tp)
tr_float = TrackedReal(v_float, d, tp)

@test hash(tr_float) === hash(v_float)
@test hash(tr_float, hash(1)) === hash(v_float, hash(1))

@test deepcopy(tr_float) === tr_float
@test copy(tr_float) === tr_float

@test samefields(float(tr_int), TrackedReal{Float64,Float64,Nothing}(float(v_int)))
@test float(tr_float) === tr_float

@test samefields(one(tr_float), typeof(tr_float)(one(v_float)))

@test samefields(zero(tr_float), typeof(tr_float)(zero(v_float)))

tr_rand = rand(TrackedReal{Int,Float64,Nothing})
@test samefields(tr_rand, TrackedReal{Int,Float64,Nothing}(ReverseDiff.value(tr_rand)))

tr_rand = rand(MersenneTwister(1), TrackedReal{Int,Float64,Nothing})
@test samefields(tr_rand, TrackedReal{Int,Float64,Nothing}(ReverseDiff.value(tr_rand)))

@test eps(tr_float) === eps(v_float)
@test eps(typeof(tr_float)) === eps(Float64)

@test floor(tr_float) === floor(v_float)
@test floor(Int, tr_float) === floor(Int, v_float)

@test ceil(tr_float) === ceil(v_float)
@test ceil(Int, tr_float) === ceil(Int, v_float)

@test trunc(tr_float) === trunc(v_float)
@test trunc(Int, tr_float) === trunc(Int, v_float)

@test round(tr_float) === round(v_float)
@test round(Int, tr_float) === round(Int, v_float)

################
# track/track! #
################

v, d = rand(), rand()
varr, darr = rand(3), rand(3)
tp = InstructionTape()

@test samefields(ReverseDiff.track(v, tp), TrackedReal(v, zero(v), tp))
@test samefields(ReverseDiff.track(v, Int, tp), TrackedReal(v, zero(Int), tp))

@test samefields(ReverseDiff.track(varr, tp), TrackedArray(varr, fill!(similar(varr), 0), tp))
@test samefields(ReverseDiff.track(varr, Int, tp), TrackedArray(varr, fill(0, size(varr)), tp))

tr = TrackedReal(v, d, tp)
x = rand()
ReverseDiff.track!(tr, x)

@test samefields(tr, TrackedReal(x, zero(x), tp))

ta = TrackedArray(varr, darr, tp)
x = rand(3)
ReverseDiff.track!(ta, x)

@test samefields(ta, TrackedArray(x, fill!(similar(x), 0), tp))

ta = TrackedArray(varr, darr, tp)
trs = similar(ta)
tp2 = InstructionTape()
@test isa(trs, Vector{eltype(ta)})
trs = similar(ta, TrackedReal{Float64,Float64,Nothing})
track!(trs, ta, tp2)
for i in eachindex(trs)
    @test samefields(trs[i], track(varr[i], tp2))
end

end


text:
module TestValue
using Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, refcode,  reftype, leveltype

@testset "leveltype on non CategoricalValue types" begin
    @test leveltype("abc") === String
    @test leveltype(String) === String
    @test leveltype(1.0) === Float64
    @test leveltype(Float64) === Float64
end

@testset "CategoricalValue on DefaultRefType pool in sorted order" begin
    pool = CategoricalPool(
        Dict(
            "a" => DefaultRefType(1),
            "b" => DefaultRefType(2),
            "c" => DefaultRefType(3),
        )
    )

    for i in 1:3
        x = CategoricalValue(pool, i)

        @test leveltype(x) === String
        @test leveltype(typeof(x)) === String
        @test reftype(x) === DefaultRefType
        @test reftype(typeof(x)) === DefaultRefType
        @test x isa CategoricalValue{String, DefaultRefType}

        @test refcode(x) === DefaultRefType(i)
        @test CategoricalArrays.pool(x) === pool

        @test typeof(x)(x) === x
        @test CategoricalValue(pool, UInt8(i)) == x
    end
end

@testset "CategoricalValue on UInt8 pool in custom order" begin
    pool = CategoricalPool(
        Dict(
            "a" => UInt8(3),
            "b" => UInt8(2),
            "c" => UInt8(1),
        )
    )

    for i in 1:3
        x = CategoricalValue(pool, i)

        @test leveltype(x) === String
        @test leveltype(typeof(x)) === String
        @test reftype(x) === UInt8
        @test reftype(typeof(x)) === UInt8
        @test x isa CategoricalValue{String, UInt8}

        @test refcode(x) === UInt8(i)
        @test CategoricalArrays.pool(x) === pool

        @test typeof(x)(x) === x
        @test CategoricalValue(pool, UInt32(i)) == x
    end
end

@testset "constructor from other value" begin
    pool = CategoricalPool([2, 3, 1])
    arr = CategoricalVector{Int}(DefaultRefType[2, 1, 3], pool)
    for x in (CategoricalValue(pool, 1), arr, view(arr, 2:3))
        for (i, v) in enumerate(levels(pool))
            @test CategoricalValue(v, x) ===
                CategoricalValue(float(v), x) ===
                CategoricalValue(CategoricalValue(pool, i), x) ===
                CategoricalValue(pool, i)
        end

        @test_throws ArgumentError CategoricalValue(4, x)
        @test_throws ArgumentError CategoricalValue(missing, x)
    end
end

end


text:
export SGD

immutable SGD <: SolverMethod
end

const defaultDict = @compat Dict(:lr_policy => LRPolicy.Fixed(0.01), :mom_policy => MomPolicy.Fixed(0.))

validate_parameters(method::SGD, params::SolverParameters) = validate_parameters(params, :lr_policy, :mom_policy)
make_solver_parameters(method::SGD; kwargs...) = merge(make_solver_parameters(),
                                                       defaultDict,
                                                       SolverParameters(kwargs))

type SGDSolverState <: InternalSolverState
    learning_rate :: Float64
    momentum      :: Float64
    param_states  :: Vector{LayerState}
    param_history :: Vector{Vector{Blob}}
    last_momentum :: Float64
end

type SGDSolverSnapshot <: SolverStateSnapshot
    iteration     :: Int
    obj_val       :: Float64
    learning_rate :: Float64
    momentum      :: Float64
end

############################################################
#  API functions to be implemented by each solver instance
############################################################

function snapshot(state::SolverState{SGDSolverState})
    SGDSolverSnapshot(state.iter, state.obj_val,
                      state.internal.learning_rate, state.internal.momentum)
end

solver_state(net::Net, snapshot::SGDSolverSnapshot) = begin
    SolverState{SGDSolverState}(snapshot.iteration,
                                snapshot.obj_val,
                                Dict(),
                                SGDSolverState(net, snapshot.learning_rate, snapshot.momentum))
end

solver_state(method::SGD, net::Net, params::SolverParameters) = begin
    SolverState(SGDSolverState(net,
                               get_learning_rate(params[:lr_policy]),
                               get_momentum(params[:mom_policy])))
end

SGDSolverState(net::Net, learning_rate::Float64, momentum::Float64) = begin
  param_states = updatable_layer_states(net)

  param_history = Array(Vector{Blob}, length(param_states))

  for i = 1:length(param_states)
    state = param_states[i]
    param_history[i] = [make_zero_blob(net.backend, eltype(x.blob),size(x.blob)...) for x in state.parameters]
  end
a
  return SGDSolverState(learning_rate, momentum, param_states, param_history, momentum)
end

function shutdown(state::SolverState{SGDSolverState})
  map(x -> map(destroy, x), state.internal.param_history)
end

function update(solver::Solver{SGD}, net::Net, state::SolverState{SGDSolverState})
  solver_state.internal.learning_rate = get_learning_rate(solver.params[:lr_policy], solver_state)
  solver_state.internal.momentum = get_momentum(solver.params[:mom_policy], solver_state)

  for i = 1:length(state.internal.param_states)
    layer_state   = state.internal.param_states[i]
    history = state.internal.param_history[i]
    for j = 1:length(layer_state.parameters)
      hist_blob = history[j]
      gradient  = layer_state.parameters[j].gradient
      data_type = eltype(hist_blob)

      # to keep iteration <-> momentum correspondence consistent with the Nesterov solver,
      # we use the momentum from the last iteration
      update_parameters!(net, solver.method, layer_state.parameters[j].learning_rate * state.internal.learning_rate,
          state.internal.last_momentum, layer_state.parameters[j].blob, hist_blob, gradient, data_type)
    end
  end
  state.internal.last_momentum = state.internal.momentum
end

function update_parameters!(net::Net{CPUBackend}, method::SGD, learning_rate, momentum, param_blob, hist_blob, gradient, data_type)
  # hist_blob = momentum * hist_blob
  BLAS.scal!(length(hist_blob), convert(data_type, momentum), pointer(hist_blob.data), 1)
  # hist_blob = -learning_rate * gradient + hist_blob
  BLAS.axpy!(length(hist_blob), convert(data_type, -learning_rate), pointer(gradient.data), 1, pointer(hist_blob.data), 1)

  # update parameter
  # param_blob += hist_blob
  BLAS.axpy!(length(hist_blob), convert(data_type, 1), pointer(hist_blob.data), 1, pointer(param_blob.data), 1)
end


text:
using GeoMakie, Makie, NASAEarthObservations, Glob
using ImageMagick

source = LonLat()
dest = WinkelTripel()

imgdir = observations("rgb/MYDAL2_M_AER_RA")

imgpaths = sort(Glob.glob(joinpath(relpath(imgdir), "*.PNG"))) # change this to your requirements

img = ImageMagick.load(imgpaths[1])

re = GeoMakie.date_regex("MYDAL2_M_AER_RA", "PNG")
titletext = Observable(join(match(re, basename(imgpaths[1])).captures, '-'))

lons = LinRange(-179.5, 179.5, size(img)[2])
lats = LinRange(89.5, -89.5, size(img)[1])

xs = [lon for lat in lats, lon in lons]
ys = [lat for lat in lats, lon in lons]

Proj4.transform!(source, dest, vec(xs), vec(ys))

scene = surface(xs, ys, zeros(size(xs)); color = img, shading = false, show_axis = false)

geoaxis!(scene, -180, 180, -90, 90; crs = (src = src, dest = dest,));

visible = node(:visible, false)

poprect = lift(scene.events.mouseposition) do mp
    FRect((mp .+ 5), 250, 40)
end
textpos = lift(scene.events.mouseposition) do mp
    Vec3f0((mp .+ 5 .+ (250/2, 40 / 2))..., 120)
end
popup = poly!(campixel(scene), poprect, raw = true, color = :white, strokewidth = 2,
        strokecolor = :black, visible = visible)

rect = popup[end]

translate!(rect, Vec3f0(0, 0, 100))

textnode = Observable(" ")

text!(popup, textnode, textsize = 30, position = textpos, color = :darkred, align = (:center, :center), raw = true, visible = visible)

on(scene.events.mouseposition) do pos
    if ispressed(scene, Mouse.left)
        visible[] = true
        pos = to_world(scene, Point2f0(scene.events.mouseposition[]))
        tpos = transform(dest, source, pos)
        textnode[] = sprint(io-> print(io, round.(Float64.(tpos), digits = 3)))
    else
        visible[] = false
    end
end

fullsc = title(scene, titletext)

save("a.png", scene); lines(rand(10))

record(fullsc, "test.mp4", imgpaths) do img

    scene.plots[1].color = GeoMakie.imflip(ImageMagick.load(img))

    # titletext[] = join(match(re, basename(img)).captures, '-')

end
lines(rand(10))


text:
module StatFiles

using ReadStat, IteratorInterfaceExtensions, TableTraits, TableTraitsUtils
using DataValues, FileIO, TableShowUtils
import IterableTables

export load, File, @format_str

struct StatFile
    filename::String
end

function Base.show(io::IO, source::StatFile)
    file_ext = lowercase(splitext(source.filename)[2])
    filetype = if file_ext==".dta"
        "STATA file"
    elseif file_ext==".sav"
        "SPSS file"
    elseif file_ext==".sas7bdat"
        "SAS file"
    else
        "file"
    end
    TableShowUtils.printtable(io, getiterator(source), filetype)
end

function Base.show(io::IO, ::MIME"text/html", source::StatFile)
    TableShowUtils.printHTMLtable(io, getiterator(source))
end

Base.Multimedia.showable(::MIME"text/html", source::StatFile) = true

function Base.show(io::IO, ::MIME"application/vnd.dataresource+json", source::StatFile)
    TableShowUtils.printdataresource(io, getiterator(source))
end

Base.Multimedia.showable(::MIME"application/vnd.dataresource+json", source::StatFile) = true

function fileio_load(f::FileIO.File{FileIO.format"Stata"})
    return StatFile(f.filename)
end

function fileio_load(f::FileIO.File{FileIO.format"SPSS"})
    return StatFile(f.filename)
end

function fileio_load(f::FileIO.File{FileIO.format"SAS"})
    return StatFile(f.filename)
end

IteratorInterfaceExtensions.isiterable(x::StatFile) = true
TableTraits.isiterabletable(x::StatFile) = true

function IteratorInterfaceExtensions.getiterator(file::StatFile)
    filename, extension = splitext(file.filename)
    if extension==".dta"
        df = read_dta(file.filename)
    elseif extension==".por"
        df = read_por(file.filename)
    elseif extension==".sav"
        df = read_sav(file.filename)
    elseif extension==".sas7bdat"
        df = read_sas7bdat(file.filename)
    else
        error("Unknown file type.")
    end

    it = TableTraitsUtils.create_tableiterator(df.data, df.headers)

    return it
end

end # module


text:
"""
    FieldVector

Simple vector wrapper around a padded multi-dimensional array. Used to
cast field/image data as a vector for linear algebra.
"""
struct FieldVector{T,N,Dt<:AbstractArray{T,N}} <: AbstractVector{T}
    data :: Dt # data array
    R :: CartesianIndices{N,NTuple{N,UnitRange{Int}}}
    FieldVector(data::AbstractArray{T,N},R=inside(data)) where {T,N} =
        new{T,N,typeof(data)}(data,R)
end
inside(dims::NTuple{N}) where {N} = CartesianIndices(ntuple(i-> 2:dims[i]-1,N))
inside(a::AbstractArray) = inside(size(a))

Base.size(x::FieldVector) = (length(x.R),)
Base.IndexStyle(::Type{<:FieldVector}) = IndexLinear()
Base.@propagate_inbounds Base.getindex(x::FieldVector, i::Int) = x.data[x.R[i]]
Base.@propagate_inbounds Base.getindex(x::FieldVector, I::CartesianIndex) = x.data[I]
Base.@propagate_inbounds Base.setindex!(x::FieldVector, y , i::Int) = x.data[x.R[i]]=y
Base.@propagate_inbounds Base.setindex!(x::FieldVector, y , I::CartesianIndex) = x.data[I]=y
Base.similar(x::FieldVector) = FieldVector(similar(x.data),x.R)
Base.zero(x::FieldVector,T=eltype(x)) = FieldVector(zeros(T,size(x.data)),x.R)

import LinearAlgebra: dot,norm
@fastmath function dot(x::FieldVector,y::FieldVector)
    s = zero(eltype(x))
    @inbounds @simd for I ∈ x.R
        s+= y[I]*x[I]
    end
    s
end
norm(x::FieldVector) = √dot(x,x)

macro loop(ex)
    a,I = ex.args[1].args
    return quote
        @inbounds @simd for $I ∈ $a.R
            $ex
        end
    end |> esc
end


text:
# All 3D Pumi mesh tests go here

# MPI Declarations
opts = Dict{ASCIIString, Any}()
opts["order"] = 1
opts["dimensions"] = 3
opts["use_DG"] = true
opts["Tsbp"] = Float64
opts["Tmsh"] = Complex128
opts["operator_type"] = "SBPOmega"
opts["smb_name"] = "../src/mesh_files/tet8cube.smb"
opts["dmg_name"] = ".null"

opts["numBC"] = 3
opts["BC1"] = [ 0, 1, 2, 3]
opts["BC1_name"] = "ExpBC"
opts["BC2"] = [4]
opts["BC2_name"] = "noPenetrationBC"
opts["BC3"] = [10]
opts["BC3_name"] = "noPenetrationBC"

opts["coloring_distance"] = 2 # 0 For CG Mesh 2 for DG Mesh
opts["jac_type"] = 2
opts["jac_method"] = 2
opts["run_type"] = 5

Tsbp = opts["Tsbp"]
Tmsh = opts["Tmsh"]
shape_type = 2
sbp = getTetSBPOmega(degree=opts["order"], Tsbp=Tsbp)
ref_verts = sbp.vtx
face_verts = SummationByParts.SymCubatures.getfacevertexindices(sbp.cub)
topo = ElementTopology{3}(face_verts)
sbpface = TetFace{Tsbp}(opts["order"], sbp.cub, ref_verts)
mesh = PumiMeshDG3(Tmsh, sbp, opts, sbpface, topo, dofpernode=dofpernode,
                     shape_type=shape_type)

#mesh = PumiMeshDG3{Tmsh}(opts["dmg_name"], opts["smb_name"], opts["order"], sbp,
#                         opts, sbpface, topo; dofpernode=1,
#                         coloring_distance=opts["coloring_distance"],
#                         shape_type=shape_type)

orig_vert_coords = deepcopy(mesh.vert_coords)

@testset "--- Checking FFD on 3D serial DG Pumi meshes ---" begin

  @testset "Check control point manipulation with nonlinear mapping on full mesh" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0., 0., 0.] # No offset in the X & Y direction

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, true)
    writeControlPointsVTS(map)

    # Rigid body rotation
    theta = -20*pi/180  # Rotate wall coordinates by 20 degrees
    rotMat = [cos(theta) -sin(theta) 0
              sin(theta) cos(theta)  0
              0          0           1] # Rotation matrix
    # Rotate the control points
    for k = 1:map.nctl[3]
      for j = 1:map.nctl[2]
        for i = 1:map.nctl[1]
          map.cp_xyz[:,i,j,k] = rotMat*map.cp_xyz[:,i,j,k]
        end
      end
    end

    # Rigid body translation
    map.cp_xyz[1,:,:,:] += 0.2
    map.cp_xyz[2,:,:,:] += 0.3
    map.cp_xyz[3,:,:,:] += 0.5

    vertices = evalVolume(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)

    writeVisFiles(mesh, "FFD_full_body_3D")

    fname = "./testvalues/translation_plus_rotation_DG_3D_tet8cube.dat"
    # f = open(fname, "w")
    # for i = 1:length(mesh.vert_coords)
    #   println(f, mesh.vert_coords[i])
    # end
    # close(f)

    test_values = readdlm(fname)
    for i = 1:length(test_values)
      err = abs.(test_values[i] - mesh.vert_coords[i])
      @test  err  < 1e-14
    end

  end # End context("Check control point manipulation with nonlinear mapping on full mesh")

  # Reset the coordinates and mesh to the original value
  for i = 1:mesh.numEl
    update_coords(mesh, i, real(orig_vert_coords[:,:,i]))
  end
  commit_coords(mesh, sbp, opts)

  @testset "Check control point manipulation on a geometry face" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0.1, 0.1, 0.1] # No offset in the X & Y direction
#    geom_faces = opts["BC2"]
    bc_nums = [2]

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, false, bc_nums)
    writeControlPointsVTS(map)

    # Rigid body rotation
    theta = -2*pi/180  # Rotate wall coordinates by 2 degrees
    rotMat = [cos(theta) -sin(theta) 0
              sin(theta) cos(theta)  0
              0          0           1] # Rotation matrix
    # Rotate the control points
    for k = 1:map.nctl[3]
      for j = 1:map.nctl[2]
        for i = 1:map.nctl[1]
          map.cp_xyz[:,i,j,k] = rotMat*map.cp_xyz[:,i,j,k]
        end
      end
    end

    # Rigid body translation
    map.cp_xyz[1,:,:,:] += 0.02
    map.cp_xyz[2,:,:,:] += 0.03
    map.cp_xyz[3,:,:,:] += 0.05

    vertices = evalSurface(map, mesh)
    commitToPumi(map, mesh, sbp, vertices, opts)

    writeVisFiles(mesh, "FFD_perturbed_face_BC2")

    fname = "./testvalues/translation_plus_rotation_DG_3D_tet8cube_face4.dat"
    # f = open(fname, "w")
    # for i = 1:length(mesh.vert_coords)
    #  println(f, mesh.vert_coords[i])
    # end
    # close(f)

    test_values = readdlm(fname)
    for i = 1:length(test_values)
      err = abs.(test_values[i] - mesh.vert_coords[i])
      @test  err  < 1e-14
    end

  end # End context("Check control point manipulation on a geometry face")

  # Reset the coordinates and mesh to the original value
  for i = 1:mesh.numEl
    update_coords(mesh, i, real(orig_vert_coords[:,:,i]))
  end
  commit_coords(mesh, sbp, opts)

  @testset "--- Checking evaldXdControlPointProduct for 3D DG Mesh ---" begin

    ndim = mesh.dim
    order = [4,4,2]  # Order of B-splines in the 3 directions
    nControlPts = [4,4,2]
    offset = [0.5, 0.5, 0.5] # No offset in the X & Y direction
#    geom_faces = opts["BC2"]
    bc_nums = [2]

    # Create a mapping object using nonlinear mapping
    map, box = initializeFFD(mesh, sbp, order, nControlPts, offset, false, bc_nums)

    fill!(map.work, 0.0)
    pert = 1e-6

    # Create seed vector
    # - Get original wall coordinates
    orig_wallCoords = FFD.getUniqueWallCoordsArray(mesh, bc_nums)
    nwall_faces = FFD.getnWallFaces(mesh, bc_nums)
    # Xs_bar = randn(3, size(orig_wallCoords,2))
    Xs_bar = ones(3, size(orig_wallCoords,2))
    cp_xyz_bar = zeros(map.cp_xyz)
    evaldXdControlPointProduct(map, mesh, vec(Xs_bar))
    for i = 1:size(map.work, 4)
      for j = 1:size(map.work, 3)
        for k = 1:size(map.work, 2)
            cp_xyz_bar[1:3,k,j,i] = map.work[1:3, k,j,i]
        end
      end
    end

    # fname = "./testvalues/evaldXdControlPointProduct_tet8cube.dat"
    # f = open(fname, "w")
    # for i = 1:length(map.work)
    #   println(f, map.work[i])
    # end
    # close(f)

    # Check against finite difference
    cp_jacobian = zeros(length(orig_wallCoords), length(map.cp_xyz))
    for i = 1:length(map.cp_xyz)
      map.cp_xyz[i] += pert
      vertices = evalSurface(map, mesh)
      commitToPumi(map, mesh, sbp, vertices, opts)
      new_wallCoords = FFD.getUniqueWallCoordsArray(mesh, bc_nums)
      cp_jacobian[:,i] = (vec(new_wallCoords) - vec(orig_wallCoords))/pert
      map.cp_xyz[i] -= pert
    end # End for i = 1:length(map.cp_xyz)
    prod_val = transpose(cp_jacobian)*vec(Xs_bar)

    error = vec(cp_xyz_bar) - prod_val
    for i = 1:length(error)
      @test isapprox( error[i], 0.0) atol=1e-8
    end

  end # End context("--- Checking evaldXdControlPointProduct for 2D DG Mesh ---")

end # End facts("--- Checking FFD on 3D serial DG Pumi meshes ---")


text:
include("TaylorModels.jl")
using TaylorModels, IntervalArithmetic

t1 = make_Taylor_model(exp, 3, 0, -0.5..0.5)
t2 = make_Taylor_model(sin, 3, 0, -0.5..0.5)

t1 + t2
t1 * t2

t = taylor_var(3, 0, -0.5..0.5)
compute_bound(t)

t = t1 + t2
compute_bound(t)


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.

struct SPResult
    x̄::Vector{Float64}
    ȳ::Dict{Int, Vector{Float64}}
    VRP::Float64
    EWS::Float64
    EVPI::Float64
    VSS::Float64
    EV::Float64
    EEV::Float64
end

problems = Vector{Tuple{StochasticModel,Vector{<:AbstractScenario},SPResult,String}}()
@info "Loading test problems..."
@info "Loading simple..."
include("simple.jl")
@info "Loading instant simple..."
include("instant_simple.jl")
@info "Loading vectorized simple..."
include("vectorized_simple.jl")
@info "Loading infeasible..."
include("infeasible.jl")
@info "Loading vectorized infeasible..."
include("vectorized_infeasible.jl")
@info "Loading farmer..."
include("farmer.jl")
@info "Loading sampler..."
include("sampler.jl")
@info "Test problems loaded. Starting test sequence."


text:
using LJH
using Test
const WT = false # run @code_warntype

# clean out the artifacts directory
isdir("artifacts") && rm("artifacts",recursive=true)
mkdir("artifacts")

testfiles = [
    "ljh",
    "ljh3",
    "ljhutil",
    ]

# In later versions of Julia (≥1.3 maybe?), you can use
# Pkg.test("Pope"; test_args=["ljh", "ljh3"]) to select a subset of
# tests, but this won't affect ARGS in earlier versions.
if length(ARGS) > 0
    tests = ARGS
end

@testset "LJH" begin

for t in testfiles
    fp = joinpath(dirname(@__FILE__), "$t.jl")
    println("Running $fp ...")
    include(fp)
end

end # @testset


text:
# an implementation of Uniform Manifold Approximation and Projection
# for Dimension Reduction, L. McInnes, J. Healy, J. Melville, 2018.

export UMAP, optimize_embedding!

"""
    struct UMAP

The UMAP model struct

# Properties

- `graph`: The fuzzy simplicial set that represents the all knn graph
- `embedding`: The embedding projection
- `k` the number of neighbors used to create the model
- `a` and `b`: parameters to ensure an well distributed and smooth projection (from `min_dist` and `spread` arguments in `fit`)
- `index`: the search index, it can be nothing if the model is handled directly with precomputed `knns` and `dists` matrices
"""
struct UMAP
    graph::SparseMatrixCSC{Float32,Int32}
    embedding::Matrix{Float32}
    k::Int
    a::Float32
    b::Float32
    index::Union{AbstractSearchContext,Nothing}   # this object is used with function barriers
end

const SMOOTH_K_TOLERANCE = 1e-5

"""
    fit(Type{UMAP}, knns, dists; <kwargs>) -> UMAP object

Create a model representing the embedding of data `(X, dist)` into `maxoutdim`-dimensional space.
Note that `knns` and `dists` jointly specify the all `k` nearest neighbors of ``(X, dist)``,
these results must not include self-references. See the `allknn` method in `SimilaritySearch`.

# Arguments

- `knns`: A ``(k, n)`` matrix of integers (identifiers).
- `dists`: A ``(k, n)`` matrix of floating points (distances).

It uses all available threads for the projection.

# Keyword Arguments
- `maxoutdim::Integer=2`: The number of components in the embedding
- `n_epochs::Integer = 300`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `learning_rate_decay::Real = 0.9`: how much `learning_rate` is updated on each epoch `(learning_rate *= learning_rate_decay)` (a minimum value is also considered as 1e-6)
- `layout::AbstractLayout = SpectralLayout()`: how to initialize the output embedding
- `min_dist::Real = 0.1`: the minimum spacing of points in the output embedding
- `spread::Real = 1`: the effective scale of embedded points. Determines how clustered embedded points are in combination with `min_dist`.
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
- `a = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.
- `b = nothing`: this controls the embedding. By default, this is determined automatically by `min_dist` and `spread`.

"""
function fit(::Type{UMAP},
    knns::Matrix{<:Integer},
    dists::Matrix{<:AbstractFloat};
    maxoutdim::Integer = 2,
    n_epochs::Integer = 50,
    learning_rate::Real = 1f0,
    learning_rate_decay::Real = 0.9f0,
    layout::AbstractLayout = SpectralLayout(),
    min_dist::Real = 0.1f0,
    spread::Real = 1f0,
    set_operation_ratio::Real = 1f0,
    local_connectivity::Integer = 1,
    repulsion_strength::Float32 = 1f0,
    neg_sample_rate::Integer = 5
)
    min_dist = convert(Float32, min_dist)
    spread = convert(Float32, spread)
    set_operation_ratio = convert(Float32, set_operation_ratio)
    learning_rate = convert(Float32, learning_rate)
    learning_rate_decay = convert(Float32, learning_rate_decay)
    repulsion_strength = convert(Float32, repulsion_strength)
    size(knns) == size(dists) || throw(ArgumentError("knns and dists must have the same size"))
    maxoutdim > 0 || throw(ArgumentError("maxoutdim must be greater than 0"))

    # argument checking
    n_epochs > 0 || throw(ArgumentError("n_epochs must be greater than 0"))
    learning_rate > 0 || throw(ArgumentError("learning_rate must be greater than 0"))
    min_dist > 0 || throw(ArgumentError("min_dist must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1 || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0 || throw(ArgumentError("local_connectivity must be greater than 0"))

    n_neighbors, n = size(knns)
    println(stderr, "*** computing graph")
    timegraph = @elapsed graph = fuzzy_simplicial_set(knns, dists, n, local_connectivity, set_operation_ratio)
    println(stderr, "*** layout embedding $(typeof(layout))")
    timeinit = @elapsed embedding = initialize_embedding(layout, graph, knns, dists, maxoutdim)
    println(stderr, "*** fit ab / embedding")
    a, b = fit_ab(min_dist, spread, nothing, nothing)
    println(stderr, "*** opt embedding")
    timeopt = @elapsed embedding = optimize_embedding(graph, embedding, embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, a, b; learning_rate_decay)
    # TODO: if target variable y is passed, then construct target graph
    #       in the same manner and do a fuzzy simpl set intersection
    println(stderr,
    """
UMAP construction time cost report:
- fuzzy graph: $timegraph
- embedding init: $timeinit
- embedding opt: $timeopt
""")
    UMAP(graph, embedding, n_neighbors, a, b, nothing)
end

"""
    fit(::Type{<:UMAP}, index_or_data;
        k=15,
        dist::SemiMetric=L2Distance,
        kwargs...)

Wrapper for `fit` that computes `n_nearests` nearest neighbors on `index_or_data` and passes these and `kwargs` to regular `fit`.

# Arguments

- `index_or_data`: an already constructed index (see `SimilaritySearch`), a matrix, or an abstact database (`SimilaritySearch`)
- `k=15`: number of neighbors to compute
- `dist=L2Distance()`: A distance function (see `Distances.jl`)
"""
function fit(
        t::Type{<:UMAP},
        index_or_data::Union{<:AbstractSearchContext,<:AbstractDatabase,<:AbstractMatrix};
        k::Integer=15,
        dist::SemiMetric=L2Distance(),
        kwargs...
        )
    index = if index_or_data isa AbstractMatrix
        db = convert(AbstractDatabase, index_or_data)
        ExhaustiveSearch(; dist, db)
    elseif index_or_data isa AbstractDatabase
        ExhaustiveSearch(; dist, db=index_or_data)
    else
        index_or_data
    end

    0 < k < length(index) || throw(ArgumentError("number of neighbors must be in 0 < k < number of points"))

    @time knns, dists = allknn(index, k; parallel=true)
    m = fit(t, knns, dists; kwargs...)
    UMAP(m.graph, m.embedding, m.k, m.a, m.b, index)
end

"""
    optimize_embedding!(model::UMAP; <kwargs>)

Improves the internal embedding of the model refining with more epochs

# Keyword arguments
- `n_epochs=50`
- `learning_rate::Real = 0.1f0`
- `learning_rate_decay::Real = 0.9f0`
- `repulsion_strength::Float32 = 1f0`
- `neg_sample_rate::Integer = 5`

"""
function optimize_embedding!(U::UMAP;
    n_epochs=50,
    learning_rate::Real = 0.1f0,
    learning_rate_decay::Real = 0.9f0,
    repulsion_strength::Float32 = 1f0,
    neg_sample_rate::Integer = 5,
)
    optimize_embedding(U.graph, U.embedding, U.embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, U.a, U.b; learning_rate_decay)
    U
end

"""
    fit(UMAP::UMAP, maxoutdim; <kwargs>)

Reuses a previously computed model with a different number of components

# Keyword arguments

- `n_epochs=50`: number of epochs to run
- `learning_rate::Real = 1f0`: initial learning rate
- `learning_rate_decay::Real = 0.9f0`: how learning rate is adjusted per epoch `learning_rate *= learning_rate_decay`
- `repulsion_strength::Float32 = 1f0`: repulsion force (for negative sampling)
- `neg_sample_rate::Integer = 5`: how many negative examples per object are used.

"""
function fit(
        U::UMAP, maxoutdim::Integer;
        n_epochs=50,
        learning_rate::Real = 1f0,
        learning_rate_decay::Real = 0.9f0,
        repulsion_strength::Float32 = 1f0,
        neg_sample_rate::Integer = 5,
        graph = U.graph,
        a = U.a,
        b = U.b
    )

    k, n = size(U.embedding)
    embedding = if k >= maxoutdim
        U.embedding[1:maxoutdim, :]
    else
        vcat(U.embedding, rand(-10f0:eps(Float32):10f0, maxoutdim-k, n))
    end

    learning_rate = convert(Float32, learning_rate)
    learning_rate_decay = convert(Float32, learning_rate_decay)
    repulsion_strength = convert(Float32, repulsion_strength)
    embedding = optimize_embedding(graph, embedding, embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, a, b; learning_rate_decay)
    UMAP(graph, embedding, U.k, a, b, U.index)
end

"""
    predict(model::UMAP)

Returns the internal embedding (the entire dataset projection)
"""
predict(model::UMAP) = model.embedding

"""
    predict(model::UMAP, Q::AbstractDatabase; k::Integer=15, kwargs...)
    predict(model::UMAP, knns, dists; <kwargs>) -> embedding

Use the given model to embed new points ``Q`` into an existing embedding produced by ``(X, dist)``.
The second function represent `Q` using its `k` nearest neighbors in `X` under some distance function (`knns` and `dists`)
See `searchbatch` in `SimilaritySearch` to compute both (also for `AbstractDatabase` objects).

# Arguments
- `model`: The fitted model
- `knns`: matrix of identifiers (integers) of size ``(k, |Q|)``
- `dists`: matrix of distances (floating point values) of size ``(k, |Q|)``

Note: the number of neighbors `k` (embedded into knn matrices) control the embedding. Larger values capture more global structure in the data, while small values capture more local structure.

# Keyword Arguments

- `n_epochs::Integer = 30`: the number of training epochs for embedding optimization
- `learning_rate::Real = 1`: the initial learning rate during optimization
- `learning_rate_decay::Real = 0.8`: A decay factor for the `learning_rate` param (on each epoch)
- `set_operation_ratio::Real = 1`: interpolates between fuzzy set union and fuzzy set intersection when constructing the UMAP graph (global fuzzy simplicial set). The value of this parameter should be between 1.0 and 0.0: 1.0 indicates pure fuzzy union, while 0.0 indicates pure fuzzy intersection.
- `local_connectivity::Integer = 1`: the number of nearest neighbors that should be assumed to be locally connected. The higher this value, the more connected the manifold becomes. This should not be set higher than the intrinsic dimension of the manifold.
- `repulsion_strength::Real = 1`: the weighting of negative samples during the optimization process.
- `neg_sample_rate::Integer = 5`: the number of negative samples to select for each positive sample. Higher values will increase computational cost but result in slightly more accuracy.
"""
function predict(model::UMAP,
                knns::AbstractMatrix{<:Integer},
                dists::AbstractMatrix{<:AbstractFloat};
                n_epochs::Integer = 30,
                learning_rate::Real = 1.0,
                learning_rate_decay::Real = 0.8,
                set_operation_ratio::Real = 1.0,
                local_connectivity::Integer = 1,
                repulsion_strength::Real = 1.0,
                neg_sample_rate::Integer = 5
    )

    set_operation_ratio = convert(Float32, set_operation_ratio)
    learning_rate = convert(Float32, learning_rate)
    repulsion_strength = convert(Float32, repulsion_strength)
    learning_rate_decay = convert(Float32, learning_rate_decay)

    # argument checking
    learning_rate > 0                                || throw(ArgumentError("learning_rate must be greater than 0"))
    0 ≤ set_operation_ratio ≤ 1                      || throw(ArgumentError("set_operation_ratio must lie in [0, 1]"))
    local_connectivity > 0                           || throw(ArgumentError("local_connectivity must be greater than 0"))

    n_epochs = max(1, n_epochs)
    # main algorithm
    n = size(model.embedding, 2)
    graph = fuzzy_simplicial_set(knns, dists, n, local_connectivity, set_operation_ratio, false)
    E = initialize_embedding(graph, model.embedding)
    optimize_embedding(graph, E, model.embedding, n_epochs, learning_rate, repulsion_strength, neg_sample_rate, model.a, model.b; learning_rate_decay)
end

function predict(model::UMAP, Q;
        k::Integer=15,
        kwargs...
    )
    model.index === nothing && throw(ArgumentError("this UMAP model doesn't support solving knn queries since `model.index == nothing` please use the alternative function that accepts `knns` and `dists` matrices"))
    Q = convert(AbstractDatabase, Q)
    0 < k <= length(Q) || throw(ArgumentError("number of neighbors must be in 0 < k <= number of points"))
    knns, dists = searchbatch(model.index, Q, k; parallel=true)
    predict(model, knns, dists; kwargs...)
end

"""
    fuzzy_simplicial_set(knns, dists, n_points, local_connectivity, set_op_ratio, apply_fuzzy_combine=true) -> membership_graph::SparseMatrixCSC,

Construct the local fuzzy simplicial sets of each point represented by its distances
to its `k` nearest neighbors, stored in `knns` and `dists`, normalizing the distances
on the manifolds, and converting the metric space to a simplicial set.
`n_points` indicates the total number of points of the original data, while `knns` contains
indices of some subset of those points (ie some subset of 1:`n_points`). If `knns` represents
neighbors of the elements of some set with itself, then `knns` should have `n_points` number of
columns. Otherwise, these two values may be inequivalent.
If `apply_fuzzy_combine` is true, use intersections and unions to combine
fuzzy sets of neighbors (default true).

The returned graph will have size (`n_points`, size(knns, 2)).
"""
function fuzzy_simplicial_set(knns::AbstractMatrix,
                              dists::AbstractMatrix,
                              n_points::Integer,
                              local_connectivity,
                              set_operation_ratio,
                              apply_fuzzy_combine=true)
    # @time σs, ρs = smooth_knn_dists(dists, k, local_connectivity)
    # @time rows, cols, vals = compute_membership_strengths(knns, dists, σs, ρs)
    rows, cols, vals = compute_membership_strengths(knns, dists, local_connectivity)
    # transform uses n_points != size(knns, 2)
    fs_set = sparse(rows, cols, vals, n_points, size(knns, 2))

    if apply_fuzzy_combine
        dropzeros!(combine_fuzzy_sets(fs_set, convert(Float32, set_operation_ratio)))
    else
        dropzeros!(fs_set)
    end
end

"""
    smooth_knn_dists(dists, k, local_connectivity; <kwargs>) -> knn_dists, nn_dists

Compute the distances to the nearest neighbors for a continuous value `k`. Returns
the approximated distances to the kth nearest neighbor (`knn_dists`)
and the nearest neighbor (nn_dists) from each point.
"""
function smooth_knn_dists_vector(col::AbstractVector, k::Integer, local_connectivity::Integer; niter::Integer=64, bandwidth::Float32=1f0)
    local_connectivity = max(1, min(k, local_connectivity))
    ρ = _find_first_non_zero(col, local_connectivity) #col[local_connectivity]
    σ = smooth_knn_dist_opt_binsearch(col, ρ, k, bandwidth, niter)
    ρ, σ
end

# check if this is necessary with SimilaritySearch, maybe just adding eps for near duplicates, d(u,v) = 0, is enough.
# also check tests
function _find_first_non_zero(v, sp)
    @inbounds for i in sp:length(v)
        v[i] != 0 && return v[i]
    end

    v[1]
end

# calculate sigma for an individual point
function smooth_knn_dist_kernel(dists, ρ, mid)
    D::Float32 = 0.0
    invmid = -1f0/mid
    @fastmath @inbounds @simd for d in dists
        d = d - ρ
        D += d > 0 ? exp(d * invmid) : 1f0
    end

    D
end

@fastmath function smooth_knn_dist_opt_binsearch(dists::AbstractVector, ρ, k, bandwidth, niter)
    target = bandwidth * log2(k)
    lo::Float32 = 0
    mid::Float32 = 1
    hi::Float32 = Inf32

    for _ in 1:niter
        psum = smooth_knn_dist_kernel(dists, ρ, mid)
        abs(psum - target) < SMOOTH_K_TOLERANCE && break
        if psum > target
            hi = mid
            mid = (lo + hi) * 0.5f0
        else
            lo = mid
            if hi === Inf32
                mid += mid
            else
                mid = (lo + hi) * 0.5f0
            end
        end
    end

    # TODO: set according to min k dist scale
    mid
end

"""
    compute_membership_strengths(knns, dists, local_connectivity) -> rows, cols, vals

Compute the membership strengths for the 1-skeleton of each fuzzy simplicial set.
"""
function compute_membership_strengths(knns::AbstractMatrix, dists::AbstractMatrix, local_connectivity)
    n = length(knns)
    rows = Vector{Int32}(undef, n)
    cols = Vector{Int32}(undef, n)
    vals = Vector{Float32}(undef, n)
    n_neighbors, n = size(knns) # WARN n is now different

    Threads.@threads for i in 1:n
        D = @view dists[:, i]
        ρ, σ = smooth_knn_dists_vector(D, n_neighbors, local_connectivity)
        invσ = -1f0 / σ
        ii = (i-1) * n_neighbors
        @inbounds for k in 1:n_neighbors
            # if i == knns[j, i] # THIS CONDITION NEVER HAPPENS WITH SimilaritySearch
            d = D[k] - ρ
            d = d > 0 ? exp(d * invσ) : 1f0
            iii = ii + k
            cols[iii] = i
            rows[iii] = knns[k, i]
            vals[iii] = d
        end
    end

    rows, cols, vals
end


text:
###############################################################################
#
# Method to fit soil albedo
#
###############################################################################
"""
    abstract type AbstractAlbedoFitting

Hierachy of AbstractAlbedoFitting
- [`FourBandsFittingCurve`](@ref)
- [`FourBandsFittingHybrid`](@ref)
- [`FourBandsFittingPoint`](@ref)
- [`TwoBandsFittingCurve`](@ref)
- [`TwoBandsFittingHybrid`](@ref)
- [`TwoBandsFittingPoint`](@ref)
"""
abstract type AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use all four GSV PCA bands to fit two flat lines
"""
struct FourBandsFittingCurve <: AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use all four GSV PCA bands to fit a mean and a flat curve
"""
struct FourBandsFittingHybrid <: AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use all four GSV PCA bands to fit 2 mean values
"""
struct FourBandsFittingPoint <: AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use 2 GSV PCA bands to fit two flat lines
"""
struct TwoBandsFittingCurve <: AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use 2 GSV PCA bands to fit a mean and a flat curve
"""
struct TwoBandsFittingHybrid <: AbstractAlbedoFitting end

"""
$(TYPEDEF)

Method to use 2 GSV PCA bands to fit 2 mean values
"""
struct TwoBandsFittingPoint <: AbstractAlbedoFitting end

###############################################################################
#
# Calculate soil albedo based on color class and soil moisture
#
###############################################################################
"""
    soil_albedos(color::Int, swc::FT) where {FT<:AbstractFloat}

Calculate soil broad band albedos for PAR and NIR, given
- `color` Soil color class, 1~20
- `swc` Soil volumetric water content
"""
function soil_albedos(color::Int, swc::FT) where {FT<:AbstractFloat}
    # soil class must be from 1 to 20
    @assert 1 <= color <=20;

    # calculate albedos for PAR and NIR
    _delta   = max(0, FT(0.11) - FT(0.4) * swc);
    _alb_par::FT = max(SOIL_BNDS[color,1], SOIL_BNDS[color,3] + _delta);
    _alb_nir::FT = max(SOIL_BNDS[color,2], SOIL_BNDS[color,4] + _delta);

    return _alb_par,_alb_nir
end

function soil_albedos(
            color::Int,
            rswc::FT,
            rel::Bool
) where {FT<:AbstractFloat}
    # soil class must be from 1 to 20
    @assert 1 <= color <=20;

    # calculate albedos for PAR and NIR
    _alb_par::FT = SOIL_BNDS[color,1] * (1 - rswc) + rswc * SOIL_BNDS[color,3];
    _alb_nir::FT = SOIL_BNDS[color,2] * (1 - rswc) + rswc * SOIL_BNDS[color,4];

    return _alb_par,_alb_nir
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                swc::FT,
                method::AbstractAlbedoFitting;
                clm::Bool = false
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `swc` Soil volumetric water content
- `method` [`AbstractAlbedoFitting`](@ref) type fitting method
- `clm` If true, use CLM method, else use new method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            swc::FT,
            method::AbstractAlbedoFitting;
            clm::Bool = false
) where {FT<:AbstractFloat}
    @unpack iWLF = wls;

    # fit the curve only if the values mismatch
    if clm
        _ref_par,_ref_nir = soil_albedos(soil.color, swc);
    else
        _ref_par,_ref_nir = soil_albedos(soil.color, swc, true);
    end;
    if soil.ρ_PAR !== _ref_par || soil.ρ_NIR !== _ref_nir
        fit_soil_mat!(soil, wls, _ref_par, _ref_nir, method);
    end

    # update the soil information
    soil.ρ_SW_SIF .= view(soil.ρ_SW, iWLF);
    soil.ε_SW .= 1 .- soil.ρ_SW;

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::FourBandsFittingCurve
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`FourBandsFittingCurve`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::FourBandsFittingCurve
) where {FT<:AbstractFloat}
    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights using pinv
    soil.ρ_SW_raw[32:end] .= ref_NIR;
    soil.ρ_SW_raw[1:31] .= ref_PAR;
    soil.SW_vec_4 .= pinv(soil.SW_mat_raw_4) * soil.ρ_SW_raw;

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_4, soil.SW_vec_4);

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::FourBandsFittingHybrid
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`FourBandsFittingHybrid`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::FourBandsFittingHybrid
) where {FT<:AbstractFloat}
    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights using pinv
    soil.ρ_SW_raw[32:end] .= ref_NIR;
    soil.ρ_SW_raw[1:31] .= ref_PAR;
    soil.SW_vec_4 .= pinv(soil.SW_mat_raw_4) * soil.ρ_SW_raw;

    # solve for weights
    @inline _fit(x::Vector{FT}) where {FT<:AbstractFloat} = (
        mul!(soil.ρ_SW_raw, soil.SW_mat_raw_4, x);
        _diff = ( mean(soil.ρ_SW_raw[1:31]) - ref_PAR ) ^ 2 +
                mean( abs.(soil.ρ_SW_raw[32:end] .- ref_NIR) ) ^ 2;
        return -_diff
    );

    _ms = ReduceStepMethodND{FT}(x_mins = FT[-2,-2,-2,-2],
                                 x_maxs = FT[2,2,2,2],
                                 x_inis = soil.SW_vec_4,
                                 Δ_inis = FT[0.1,0.1,0.1,0.1]);
    _tol = SolutionToleranceND{FT}(FT[0.001,0.001,0.001,0.001], 50);
    _sol = find_peak(_fit, _ms, _tol);
    soil.SW_vec_4 .= _sol;

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_4, soil.SW_vec_4);

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::FourBandsFittingPoint
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`FourBandsFittingPoint`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::FourBandsFittingPoint
) where {FT<:AbstractFloat}
    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights using pinv
    soil.ρ_SW_raw[32:end] .= ref_NIR;
    soil.ρ_SW_raw[1:31] .= ref_PAR;
    soil.SW_vec_4 .= pinv(soil.SW_mat_raw_4) * soil.ρ_SW_raw;

    # solve for weights
    @inline _fit(x::Vector{FT}) where {FT<:AbstractFloat} = (
        mul!(soil.ρ_SW_raw, soil.SW_mat_raw_4, x);
        _diff = ( mean(soil.ρ_SW_raw[1:31]) - ref_PAR ) ^ 2 +
                ( mean(soil.ρ_SW_raw[32:end]) - ref_NIR ) ^ 2;
        return -_diff
    );

    _ms = ReduceStepMethodND{FT}(x_mins = FT[-2,-2,-2,-2],
                                 x_maxs = FT[2,2,2,2],
                                 x_inis = soil.SW_vec_4,
                                 Δ_inis = FT[0.1,0.1,0.1,0.1]);
    _tol = SolutionToleranceND{FT}(FT[0.001,0.001,0.001,0.001], 50);
    _sol = find_peak(_fit, _ms, _tol);
    soil.SW_vec_4 .= _sol;

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_4, soil.SW_vec_4);

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::TwoBandsFittingCurve
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`TwoBandsFittingCurve`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::TwoBandsFittingCurve
) where {FT<:AbstractFloat}
    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights using pinv
    soil.ρ_SW_raw[32:end] .= ref_NIR;
    soil.ρ_SW_raw[1:31] .= ref_PAR;
    soil.SW_vec_2 .= pinv(soil.SW_mat_raw_2) * soil.ρ_SW_raw;

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_2, soil.SW_vec_2);

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::TwoBandsFittingHybrid
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`TwoBandsFittingHybrid`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::TwoBandsFittingHybrid
) where {FT<:AbstractFloat}
    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights using pinv
    soil.ρ_SW_raw[32:end] .= ref_NIR;
    soil.ρ_SW_raw[1:31] .= ref_PAR;
    soil.SW_vec_2 .= pinv(soil.SW_mat_raw_2) * soil.ρ_SW_raw;

    # solve for weights
    @inline _fit(x::Vector{FT}) where {FT<:AbstractFloat} = (
        mul!(soil.ρ_SW_raw, soil.SW_mat_raw_2, x);
        _diff = ( mean(soil.ρ_SW_raw[1:31]) - ref_PAR ) ^ 2 +
                mean( abs.(soil.ρ_SW_raw[32:end] .- ref_NIR) ) ^ 2;
        return -_diff
    );

    _ms = ReduceStepMethodND{FT}(x_mins = FT[-2,-2],
                                 x_maxs = FT[2,2],
                                 x_inis = soil.SW_vec_2,
                                 Δ_inis = FT[0.1,0.1]);
    _tol = SolutionToleranceND{FT}(FT[0.001,0.001], 50);
    _sol = find_peak(_fit, _ms, _tol);
    soil.SW_vec_2 .= _sol;

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_2, soil.SW_vec_2);

    return nothing
end

"""
    fit_soil_mat!(
                soil::SoilOpticals{FT},
                wls::WaveLengths{FT},
                ref_PAR::FT,
                ref_NIR::FT,
                method::TwoBandsFittingPoint
    ) where {FT<:AbstractFloat}

Fit soil albedo bands parameters, given
- `soil` [`SoilOpticals`](@ref) type structure
- `wls` [`Wavelengths`](@ref) type structure
- `ref_PAR` Target albedo for PAR
- `ref_NIR` Target albedo for NIR
- `method` [`TwoBandsFittingPoint`](@ref) type fitting method
"""
function fit_soil_mat!(
            soil::SoilOpticals{FT},
            wls::WaveLengths{FT},
            ref_PAR::FT,
            ref_NIR::FT,
            method::TwoBandsFittingPoint
) where {FT<:AbstractFloat}
    @unpack dry_NIR, dry_PAR, wet_NIR, wet_PAR = soil;

    # update soil PAR and NIR albedo
    soil.ρ_PAR = ref_PAR;
    soil.ρ_NIR = ref_NIR;

    # solve for weights
    soil.SW_vec_2[1] = (ref_PAR * wet_NIR - ref_NIR * wet_PAR) /
                       (wet_NIR * dry_PAR - wet_PAR * dry_NIR);
    soil.SW_vec_2[2] = (ref_PAR * dry_NIR - ref_NIR * dry_PAR) /
                       (dry_NIR * wet_PAR - dry_PAR * wet_NIR);

    # update vectors in soil
    mul!(soil.ρ_SW, soil.SW_mat_2, soil.SW_vec_2);

    return nothing
end