text:
# Tests for conditional HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random

# Random seed
Random.seed!(11)

#######################################################################################################################
# Test invertibility

# Input
nx = 16
ny = 16
n_channel = 8
n_hidden = 32
batchsize = 2

function inv_test(nx, ny, n_channel, batchsize, rev, logdet, permute)
    print("\nConditional HINT invertibility test with permute=$(permute), logdet=$(logdet), reverse=$(rev)\n")
    CH = ConditionalLayerHINT(nx, ny, n_channel, n_hidden, batchsize; logdet=logdet, permute=permute)
    rev && (CH = reverse(CH))

    # Input image and data
    X = randn(Float32, nx, ny, n_channel, batchsize)
    Y = randn(Float32, nx, ny, n_channel, batchsize)

    # Test inverse
    Zx, Zy = CH.forward(X, Y)[1:2]
    X_, Y_ = CH.inverse(Zx, Zy)

    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-5)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-5)

    # Test backward
    ΔZx = randn(Float32, size(Zx))  # random derivative
    ΔZy = randn(Float32, size(Zx))
    ΔX_, ΔY_, X_, Y_ = CH.backward(ΔZx, ΔZy, Zx, Zy)

    @test isapprox(norm(X - X_)/norm(X), 0f0; atol=1f-5)
    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-5)

    # Test inverse Y only
    Zy = CH.forward_Y(Y)
    Y_ = CH.inverse_Y(Zy)

    @test isapprox(norm(Y - Y_)/norm(Y), 0f0; atol=1f-5)
end

function loss(CH, X, Y)
    Zx, Zy = CH.forward(X, Y)
    Z = tensor_cat(Zx, Zy)
    f = .5*norm(Z)^2
    ΔZx, ΔZy = tensor_split(Z)
    ΔX, ΔY = CH.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    isnothing(CH.C_X) ? v1grad = nothing : v1grad=CH.C_X.v1.grad
    return f, ΔX, ΔY, CH.CL_X.CL[1].RB.W1.grad, v1grad
end

function loss_logdet(CH, X, Y)
    Zx, Zy, logdet = CH.forward(X, Y)
    Z = tensor_cat(Zx, Zy)
    f = -log_likelihood(Z) - logdet
    ΔZ = -∇log_likelihood(Z)
    ΔZx, ΔZy = tensor_split(ΔZ)
    ΔX, ΔY = CH.backward(ΔZx, ΔZy, Zx, Zy)[1:2]
    isnothing(CH.C_X) ? v1grad = nothing : v1grad=CH.C_X.v1.grad
    return f, ΔX, ΔY, CH.CL_X.CL[1].RB.W1.grad, v1grad
end

function grad_test_X(nx, ny, n_channel, batchsize, rev, logdet, permute)
    print("\nConditional HINT gradient test for input with permute=$(permute), logdet=$(logdet), reverse=$(rev)\n")
    logdet ? lossf = loss_logdet : lossf = loss
    CH = ConditionalLayerHINT(nx, ny, n_channel, n_hidden, batchsize; logdet=logdet, permute=permute)
    rev && (CH = reverse(CH))

    # Input image
    X0 = randn(Float32, nx, ny, n_channel, batchsize)
    dX = randn(Float32, nx, ny, n_channel, batchsize)

    # Input data
    Y0 = randn(Float32, nx, ny, n_channel, batchsize)
    dY = randn(Float32, nx, ny, n_channel, batchsize)

    f0, gX, gY = lossf(CH, X0, Y0)[1:3]

    maxiter = 5
    h = 0.1f0
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)

    for j=1:maxiter
        f = lossf(CH, X0 + h*dX, Y0 + h*dY)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX, gX) - h*dot(dY, gY))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

function grad_test_par(nx, ny, n_channel, batchsize, rev, logdet, permute)
    print("\nConditional HINT gradient test for weights with permute=$(permute), logdet=$(logdet), reverse=$(rev)\n")
    logdet ? lossf = loss_logdet : lossf = loss
    CH0 = ConditionalLayerHINT(nx, ny, n_channel, n_hidden, batchsize; logdet=logdet, permute=permute)
    rev && (CH0 = reverse(CH0))
    CHini = deepcopy(CH0)

    # Perturbation
    X = randn(Float32, nx, ny, n_channel, batchsize)
    Y = randn(Float32, nx, ny, n_channel, batchsize)
    dW = randn(Float32, size(CH0.CL_X.CL[1].RB.W1))

    f0, gW, gv = lossf(CH0, X, Y)[[1,4,5]]

    maxiter = 5
    h = 0.1f0

    err3 = zeros(Float32, maxiter)
    err4 = zeros(Float32, maxiter)

    for j=1:maxiter
        CH0.CL_X.CL[1].RB.W1.data = CHini.CL_X.CL[1].RB.W1.data + h*dW
        f = lossf(CH0, X, Y)[1]
        err3[j] = abs(f - f0)
        err4[j] = abs(f - f0 - h*dot(gW, dW))
        print(err3[j], "; ", err4[j], "\n")
        h = h/2f0
    end

    @test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

for logdet in [true, false]
    for permute in [true, false]
        for rev in [true, false]
            inv_test(nx, ny, n_channel, batchsize, rev, logdet, permute)
            grad_test_X(nx, ny, n_channel, batchsize, rev, logdet, permute)
            grad_test_par(nx, ny, n_channel, batchsize, rev, logdet, permute)
        end
    end
end


text:
using Documenter, PowerGraphics

makedocs(
    modules = [PowerGraphics],
    format = Documenter.HTML(),
    sitename = "PowerGraphics.jl",
    authors = "Clayton Barrows",
    pages = [
        "Home" => "index.md",
        "Examples" => "examples.md", # TODO: make examples page
        "Function Index" => "api.md",
    ],
)

Documenter.deploydocs(
    repo = "github.com/NREL-SIIP/PowerGraphics.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    push_preview = true,
    versions = ["stable" => "v^", "v#.#"],
)

ext: jl

max_stars_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_issues_repo_name: NREL-SIIP/InfrastructureSystemsPlots

max_forks_repo_name: NREL-SIIP/InfrastructureSystemsPlots


text:
import Bridge: kernelr3!, R3!, target, auxiliary, constdiff, llikelihood, _b!, B!, σ!, b!

"""
    (t, obsinfo, mT, P, ITER, subsamples, xinit, pars, priorθ, priormom, updatescheme, outdir, printskip)

This is the main function call for doing either MCMC or SGD for landmark models
Backward ODEs used are in terms of the LMμ-parametrisation

## Arguments
- `t`:      time grid
- `obsinfo`: struct of type `ObsInfo`
- `mT`: vector of momenta at time T used for constructing guiding term
- `P`: target process
- `ITER`: number of iterations
- `subsamples`: vector of indices of iterations that are to be saved
- `xinit`: initial value on state at time 0 (`x0`)
- `pars`: tuning pars for updates (see `?Pars_ms()` and `?Pars_ahs()` for documentation)
- `priorθ`: prior on θ = (a,c,γ)
- `priormom`: prior on initial momenta
- `updatescheme`: vector specifying mcmc-updates
- `outdir` output directory for animation
- `printskip` skip every `printskip` iterations in writing output to console

## Returns:
- `Xsave`: saved iterations of all states at all times in tt_
- `parsave`: saved iterations of all parameter updates
- `initendstates_save`: saved iterations of initial and endstate
- `accinfo`: acceptance indicators of mcmc-updates
- `δ`: value of `(δpos, δmom)` af the final iteration of the algorithm (these are stepsize parameters for updating initial positions and momenta respectively)
- `ρ`: value of `ρinit`  at the final iteration of the algorithm
- `δa`: value of `δa` at the final iteration of the algorithm (step size for updating Hamiltonian kernel parameter `a`)
"""
function lm_mcmc(t, obsinfo, mT, P, ITER, subsamples, xinit, pars, priorθ, priormom, updatescheme, outdir, printskip)
    lt = length(t)
    StateW = PointF
    dwiener = dimwiener(P)
    nshapes = obsinfo.nshapes

    if P.n < 25
        chunksize = d*P.n
    elseif 25 <= P.n < 100
        chunksize = P.n
    else
        chunksize = P.n÷2
    end

    guidrec = [GuidRecursions(t,obsinfo)  for _ ∈ 1:nshapes]  # initialise guiding terms
    gramT_container = [gram_matrix(obsinfo.xobsT[k], P) for k ∈ 1:nshapes]  # gram matrices at observations at time T

    Paux = [auxiliary(P, State(obsinfo.xobsT[k],mT), gramT_container[k]) for k ∈ 1:nshapes] # auxiliary process for each shape
    Q = GuidedProposal(P,Paux,t,obsinfo.xobs0,obsinfo.xobsT,guidrec,nshapes,[mT for _ ∈ 1:nshapes])
    Q = update_guidrec!(Q, obsinfo)   # compute backwards recursion

    X = [initSamplePath(t, xinit) for _ ∈ 1:nshapes]
    W = [initSamplePath(t,  zeros(StateW, dwiener)) for _ ∈ 1:nshapes]
    for k ∈ 1:nshapes   sample!(W[k], Wiener{Vector{StateW}}())  end

    x = deepvec(xinit)
    q, p = split_state(xinit)
    qᵒ = deepcopy(q); pᵒ = deepcopy(q); ∇ = deepcopy(q); ∇ᵒ = deepcopy(q)

    # memory allocations, actual state at each iteration is (X, W, Q, x, ∇x) (x, ∇x are initial state and its gradient)
    Xᵒ = deepcopy(X)
    Wᵒ = initSamplePath(t,  zeros(StateW, dwiener))
    Wnew = initSamplePath(t,  zeros(StateW, dwiener))
    # sample guided proposal and compute loglikelihood (write into X)
    At = zeros(UncF, 2P.n, 2P.n)  # container for to compute a=σσ' for auxiliary process.
    ll, X = gp!(LeftRule(), X, xinit, W, Q, At; skip=sk)

    # setup containers for saving objects
    Xsave = typeof(zeros(length(t) * P.n * 2 * d * nshapes))[]
    push!(Xsave, convert_samplepath(X))
    parsave = Vector{Float64}[]
    push!(parsave, getpars(Q))
    initendstates_save = [extract_initial_and_endstate(0,X[1]) ]
    accinfo = DataFrame(fill([], length(updatescheme)+1), [updatescheme...,:iteration]) # columns are all update types + 1 column for iteration number
    acc = zeros(ncol(accinfo))

    δ = [pars.δpos, pars.δmom, pars.δsgd_mom]
    δa, δγ = pars.δa, pars.δγ
    ρ = pars.ρlowerbound
    adapt(i) = (i > 1.5*pars.adaptskip) & (mod(i,pars.adaptskip)==0)

    x0 = X[1].yy[1]
    dK = gramkernel(x0.q, P)
    if :rmrw_mom ∈ updatescheme ||  :rmmala_mom  ∈ updatescheme
        if :parameter in updatescheme
            @error "Cannot use :rmrw_mom or :rmmala_mom in case parameterupdating is also done (because this requries inv_dK te be recomputed at each iteration)"
        else
        inv_dK = inv(dK)
        end
    else
        inv_dK = 0.0 # not used then
    end

    for i ∈ 1:ITER
        k = 1
        for update ∈ updatescheme
            if update == :innov
                accinfo_, X, W, ll = update_path!(X, W, ll, Xᵒ, Wᵒ, Wnew, Q, ρ, At)
                if adapt(i)
                    ρ = adaptpcnstep(ρ, i, accinfo[!,update], pars.η, pars.adaptskip)
                end
            elseif update ∈ [:mala_mom, :tmala_mom, :rmmala_mom, :rmrw_mom]
                accinfo_, X, x, ∇, ll = update_initialstate!(X,Xᵒ,W,ll, x, qᵒ, pᵒ,∇, ∇ᵒ, Q, δ, update, priormom, (dK, inv_dK), At, chunksize)
                if adapt(i)
                    δ[2] = adapt_pospar_step(δ[2], i, accinfo[!,update], pars.η, pars.adaptskip)
                end
            elseif update ∈ [:mala_pos, :tmala_pos, :rmmala_pos]
                accinfo_ , X, x, ∇, ll= update_initialstate!(X,Xᵒ,W,ll, x, qᵒ, pᵒ,∇, ∇ᵒ, Q, δ, update, priormom, (dK, inv_dK), At, chunksize)
                if adapt(i)
                    δ[1] = adapt_pospar_step(δ[1], i, accinfo[!,update], pars.η, pars.adaptskip)
                end
            elseif update == :parameter
                Q, accinfo_, X, ll = update_pars!(X, ll, Xᵒ,W, Q, priorθ, obsinfo, δa, δγ, gramT_container, At)
                if adapt(i)
                    δa = adapt_pospar_step(δa, i, accinfo[!,update], pars.η, pars.adaptskip)
                end
            elseif update == :matching
                obsinfo, accinfo_, Q, X, ll = update_cyclicmatching(X, ll, obsinfo, Xᵒ, W, Q, At)
            elseif update == :sgd_mom
                accinfo_, X, x, ∇, ll = update_initialstate!(X,Xᵒ,W,ll, x, qᵒ, pᵒ,∇, ∇ᵒ, Q, δ, update, priormom, (dK, inv_dK), At)
            else @error "This update scheme has not been implemented"
            end
            acc[k] = accinfo_
            k += 1
        end
        acc[end] = i
        push!(accinfo, acc)

        # if mod(i,pars.adaptskip)==0 && i < 100  # adjust mT (the momenta at time T used in the construction of the guided proposal)
        #     mTvec = [X[k][lt][2].p  for k in 1:nshapes]     # extract momenta at time T for each shape
        #    Q = update_mT!(Q, mTvec, obsinfo)
        # end

        # save some of the results
        if i ∈ subsamples
            push!(Xsave, convert_samplepath(X))
        end
        push!(parsave, getpars(Q))
        push!(initendstates_save, extract_initial_and_endstate(i,X[1]))

        if mod(i,printskip) == 0
            #println();  println("iteration $i, ρ = $ρ,  δ = $δ,  δa = $δa")
            #print("Average acceptance rates over last ",  printskip, " iterations: ")
            if i >= pars.adaptskip
                ac = accinfo[end-pars.adaptskip+1:end, 1:end-1]
                println(round.([mean(x) for x in eachcol(ac)];digits=2))
            end
            #println("parameter a and γ equal: ", getpars(Q))
        end
    end
    Xsave, parsave, initendstates_save, accinfo, δ, ρ, δa
end

"""
    update_path!(X,Xᵒ,W,Wᵒ,Wnew,ll,x, Q, ρ)

Update bridges for all shapes using Crank-Nicholsen scheme with parameter `ρ` (only in case the method is mcmc).
Newly accepted bridges are written into `(X,W)`, loglikelihood on each segment is written into vector `ll`
All Wiener increments are always updated.

## Write into
- `X`: diffusion paths
- `W`: innovations
- `ll`: vector of loglikehoods (k-th element is for the k-th shape)
"""
function update_path!(X, W, ll, Xᵒ,Wᵒ, Wnew, Q, ρlowerbound, At)
    acc = Int64[]
    x0 = X[1].yy[1]
    ρ = rand(Uniform(ρlowerbound, 1.0))
    ρ_ = sqrt(1.0-ρ^2)
    # From current state (x,W) with loglikelihood ll, update to (x, Wᵒ)
    for k ∈ 1:Q.nshapes
        sample!(Wnew, Wiener{Vector{PointF}}())
        for i ∈ eachindex(W[1].yy)
            for j ∈ eachindex(W[1].yy[1])
                Wᵒ.yy[i][j] = ρ * W[k].yy[i][j] + ρ_ * Wnew.yy[i][j]
            end
        end
        llᵒ, Xᵒ[k] = gp!(LeftRule(), Xᵒ[k], x0, Wᵒ, Q, k, At; skip=sk)
        diff_ll = llᵒ - ll[k]
        if log(rand()) <= diff_ll
            for i ∈ eachindex(W[1].yy)
                X[k].yy[i] .= Xᵒ[k].yy[i]
                W[k].yy[i] .= Wᵒ.yy[i]
            end
            ll[k] = llᵒ
            push!(acc,1)
        else
            push!(acc,0)
        end
    end
    mean(acc), X, W, ll
end

"""
    slogρ!(q, p, Q, W, X, priormom,llout)

Main use of this function is to get gradient information of the loglikelihood with respect ot the initial state.

## Arguments
- `q` and `p` make up the initial state (which could be obtained from `x0 = merge_state(q,p)`)
- `Q`: GuidePropsoal!
- `W`: vector of innovations
- `X`: sample path
- `priormom`: prior on the initial momenta
- `llout`: vector where loglikelihoods for each shape are written into

## Writes into
- The guided proposal is written into `X`.
- The computed loglikelihood for each shape is written into `llout`

## Returns
- loglikelihood (summed over all shapes) + the logpriordensity of the initial momenta
"""
function slogρ!(q,p, Q, W, X, priormom, llout, At)
    lltemp, X = gp!(LeftRule(), X, q, p, W, Q, At; skip=sk)   # writes into X
    llout .= ForwardDiff.value.(lltemp)
    sum(lltemp) + logpdf(priormom, p)
end

"""
    slogρ_pos!(p, Q, W, X,priormom, llout) = (q) -> slogρ!(q, p , Q, W, X, priormom, llout)
"""
slogρ_pos!(p, Q, W, X,priormom, llout, At) = (q) -> slogρ!(q, p , Q, W, X, priormom, llout, At)

"""
    slogρ_mom!(q, Q, W, X,priormom, llout) = (p) -> slogρ!(q, p , Q, W, X, priormom, llout)
"""
slogρ_mom!(q, Q, W, X,priormom, llout, At) = (p) -> slogρ!(q, p , Q, W, X, priormom, llout, At)

"""
    update_initialstate!(X,Xᵒ,W,ll,x,xᵒ,∇, ∇ᵒ, Q::GuidedProposal, δ, update,priormom, (dK, inv_dK), At)

## Arguments
- `X`:  current iterate of vector of sample paths
- `Xᵒ`: vector of sample paths to write proposal into
- `W`:  current vector of Wiener increments
- `ll`: current value of loglikelihood
- `x`, `xᵒ`, `∇x`, `∇xᵒ`: allocated vectors for initial state and its gradient (both actual and proposed values)
- `sampler`: either sgd (not checked yet) or mcmc
- `Q`: GuidedProposal
- `δ`: vector with MALA stepsize for initial state positions (δ[1]) and initial state momenta (δ[2])
- `update`:  can be `:mala_pos`, `:mala_mom`, :rmmala_pos`, `:rmmala_mom`, ``:rmrw_mom`
- `priormom`: prior on initial momenta
- (dK, inv_dK): gramkernel at x0.q and its inverse

## Writes into / modifies
- `X`: landmarks paths
- `x`: initial state of the landmark paths
- ``∇`: gradient of loglikelihood with respect to initial state `x`
- `ll`: a vector with as k-th element the loglikelihood for the k-th shape

## Returns
-  0/1 indicator (reject/accept)
"""
function update_initialstate!(X,Xᵒ,W,ll, x, qᵒ, pᵒ,∇, ∇ᵒ,
                 Q::GuidedProposal, δ, update, priormom, (dK, inv_dK), At, chunksize)

    P = Q.target;   n = P.n
    x0 = deepvec2state(x)
    q, p = split_state(x0)
    dn = d * n

    if update ==:sgd_mom
        # @warn "Option :sgd_mom has not been checked/tested so far; don't use as yet."
        for k ∈ 1:Q.nshapes
            sample!(W[k], Wiener{Vector{PointF}}())
        end
        u = slogρ_mom!(q, Q, W, X, priormom,ll, At)
        cfg = ForwardDiff.GradientConfig(u, p, ForwardDiff.Chunk{chunksize}()) # d*P.n is maximal
        ForwardDiff.gradient!(∇, u, p, cfg) # X gets overwritten but does not change
        pᵒ = p .- δ[3] * ∇
        x0ᵒ = merge_state(q, pᵒ)
        ll, X = gp!(LeftRule(), X, x0ᵒ, W, Q, At)
        x .= deepvec(x0ᵒ)
        accepted = 1.0
    else
        accinit = 0.0
        llᵒ = copy(ll)
        if update in [:mala_pos, :tmala_pos, :rmmala_pos]
            u = slogρ_pos!(p, Q, W, X, priormom,ll, At)
            uᵒ = slogρ_pos!(p, Q, W, Xᵒ, priormom,llᵒ, At)
            cfg = ForwardDiff.GradientConfig(u, q, ForwardDiff.Chunk{chunksize}()) # d*P.n is maximal
            ForwardDiff.gradient!(∇, u, q, cfg) # X and ll get overwritten but do not change
            stepsize = sample(δ[1])
        end
        if update in [:mala_mom, :tmala_mom, :rmmala_mom]
            u = slogρ_mom!(q, Q, W, X, priormom,ll, At)
            uᵒ = slogρ_mom!(q, Q, W, Xᵒ, priormom,llᵒ, At)
            cfg = ForwardDiff.GradientConfig(u, p, ForwardDiff.Chunk{chunksize}()) # d*P.n is maximal
            ForwardDiff.gradient!(∇, u, p, cfg) # X and ll get overwritten but do not change
            #stepsize = δ[2]
            stepsize = sample(δ[2])
        end
        if update == :rmrw_mom
            u = slogρ_mom!(q, Q, W, X, priormom,ll, At)
            uᵒ = slogρ_mom!(q, Q, W, Xᵒ, priormom,llᵒ, At)
            #stepsize = δ[2]
            stepsize = sample(δ[2])
        end

        ## compute proposal and log(acceptance ratio)
        if update == :mala_pos
            ndistr = MvNormal(dn,sqrt(stepsize))
            qᵒ .= q .+ .5 * stepsize *  ∇ .+ rand(ndistr)
            cfgᵒ = ForwardDiff.GradientConfig(uᵒ, qᵒ, ForwardDiff.Chunk{chunksize}())
            ForwardDiff.gradient!(∇ᵒ, uᵒ, qᵒ, cfgᵒ)
            accinit = sum(llᵒ) - sum(ll) -
                      logpdf(ndistr, qᵒ - q - .5*stepsize * ∇) +
                      logpdf(ndistr, q - qᵒ - .5*stepsize * ∇ᵒ)
        elseif update == :tmala_pos
              ndistr = MvNormal(dn,sqrt(stepsize))
              qᵒ .= q .+ .5 * tame(stepsize *  ∇) .+ rand(ndistr)
              cfgᵒ = ForwardDiff.GradientConfig(uᵒ, qᵒ, ForwardDiff.Chunk{chunksize}())
              ForwardDiff.gradient!(∇ᵒ, uᵒ, qᵒ, cfgᵒ)
              accinit = sum(llᵒ) - sum(ll) -
                        logpdf(ndistr, qᵒ - q - .5 * tame(stepsize * ∇)) +
                        logpdf(ndistr, q - qᵒ - .5 * tame(stepsize * ∇ᵒ))
        elseif update == :mala_mom
            ndistr = MvNormal(dn, sqrt(stepsize))
            pᵒ .= p .+ .5 * stepsize * ∇ .+ rand(ndistr)
            cfgᵒ = ForwardDiff.GradientConfig(uᵒ, pᵒ, ForwardDiff.Chunk{chunksize}()) # d*P.n is maximal
            ForwardDiff.gradient!(∇ᵒ, uᵒ, pᵒ, cfgᵒ)
            accinit = sum(llᵒ) - sum(ll) -
                      logpdf(ndistr, pᵒ - p - .5 * stepsize * ∇) +
                      logpdf(ndistr, p - pᵒ - .5 * stepsize * ∇ᵒ)
      elseif update == :tmala_mom
          ndistr = MvNormal(dn, sqrt(stepsize))
          pᵒ .= p .+ .5 * tame(stepsize * ∇) .+ rand(ndistr)
          cfgᵒ = ForwardDiff.GradientConfig(uᵒ, pᵒ, ForwardDiff.Chunk{chunksize}())
          ForwardDiff.gradient!(∇ᵒ, uᵒ, pᵒ, cfgᵒ)
          accinit = sum(llᵒ) - sum(ll) -
                    logpdf(ndistr, pᵒ - p - .5 * tame(stepsize * ∇)) +
                    logpdf(ndistr, p - pᵒ - .5 * tame(stepsize * ∇ᵒ))
        elseif update == :rmmala_pos
            #dK = gramkernel(x0.q, P)
            dK = gramkernel(reinterpret(PointF,q), P)
            ndistr = MvNormal(zeros(d*n),stepsize*dK)
            qᵒ .= q .+ .5 * dK * tame(stepsize * ∇) .+ rand(ndistr)
            cfgᵒ = ForwardDiff.GradientConfig(uᵒ, qᵒ, ForwardDiff.Chunk{chunksize}())
            ForwardDiff.gradient!(∇ᵒ, uᵒ, qᵒ, cfgᵒ)
            dKᵒ = gramkernel(reinterpret(PointF,qᵒ), P)
            ndistrᵒ = MvNormal(zeros(d*n),stepsize*dKᵒ)
            accinit = sum(llᵒ) - sum(ll) -
                     logpdf(ndistr, qᵒ - q - .5 * dK * tame(stepsize * ∇)) +
                     logpdf(ndistrᵒ, q - qᵒ - .5 * dKᵒ * tame(stepsize * ∇ᵒ))
        elseif update == :rmmala_mom
             ndistr = MvNormal(stepsize*inv_dK)
             pᵒ .= p .+ .5 * stepsize * inv_dK * ∇ .+  rand(ndistr)
             cfgᵒ = ForwardDiff.GradientConfig(uᵒ, pᵒ, ForwardDiff.Chunk{chunksize}())
             ForwardDiff.gradient!(∇ᵒ, uᵒ, pᵒ, cfgᵒ) # Xᵒ gets overwritten but does not change
             accinit = sum(llᵒ) - sum(ll) -
                       logpdf(ndistr, pᵒ - p - .5*stepsize * inv_dK * ∇) +
                       logpdf(ndistr, p - pᵒ - .5*stepsize * inv_dK * ∇ᵒ)
       elseif update == :rmrw_mom
            ndistr = MvNormal(zeros(d*n),stepsize*inv_dK)
            pᵒ .= p  .+ rand(ndistr)
            uᵒ(pᵒ)  # writes into llᵒ, don't need gradient info here
            accinit = sum(llᵒ) - sum(ll)   # proposal is symmetric
        end

        # MH acceptance decision
        if log(rand()) <= accinit
            X = copypaths!(X,Xᵒ)
            ll .= llᵒ
            if update in [:mala_pos, :rmmala_pos]
                x .= deepvec(merge_state(qᵒ, p))
            end
            if update in [:mala_mom, :rmmala_mom, :rmrw_mom]
                x .= deepvec(merge_state(q, pᵒ))
            end
            accepted = 1.0
        else
            accepted = 0.0
        end
    end
    accepted, X, x, ∇, ll
end

"""
    trun(x,h) = x/max(1.0,0.5*h*norm(x))

useful function for truncated version of mala
"""
trun(x,h) = x/max(1.0,0.5*h*norm(x))

"""
    tame(x) = x/(I + norm(x))

Useful for tamed version of Langevin algorithm
"""
tame(x) = x/(I + norm(x))

"""
    update_pars!(X, Xᵒ,W, Q, x, ll, priorθ, covθprop, obsinfo, At)

For fixed Wiener increments and initial state, update parameters by random-walk-MH.

## Writes into / modifies
`X`: landmark paths under landmarks process with parameter θ
`ll`: vector of loglikelihoods (one for each shape)

## Returns
`Q, accept`: where `accept` is 1/0 according to accept/reject in the MH-step.
"""
function update_pars!(X, ll, Xᵒ, W, Q , priorθ,  obsinfo, δa, δγ, gramT_container, At)
    θ = getpars(Q)
    a = θ[1];  γ = θ[2]
    aᵒ = a * exp(δa*randn())
    γᵒ = γ * exp(δγ*randn())
    θᵒ = [aᵒ, γᵒ]
    Qᵒ = adjust_to_newpars(Q, θᵒ, obsinfo, gramT_container)
    llᵒ, Xᵒ = gp!(LeftRule(), Xᵒ, X[1].yy[1], W, Qᵒ, At; skip=sk)
    A = sum(llᵒ) - sum(ll) +
        logpdf(priorθ, θᵒ) - logpdf(priorθ, θ) +
        log(γᵒ) - log(γ) + log(aᵒ) - log(a)
    if log(rand()) <= A
        ll .= llᵒ
        X = copypaths!(X,Xᵒ)
        Q = Qᵒ
        #Q = deepcopy(Qᵒ)
        accept = 1.0
    else
        accept = 0.0
    end
    Q, accept, X, ll
end

"""
    function copypaths!(X,Xᵒ)

Write Xᵒ into X
"""
function copypaths!(X,Xᵒ)
    for k ∈ eachindex(X), i ∈ eachindex(X[1].yy)
        X[k].yy[i] .= Xᵒ[k].yy[i]
    end
    X
end

######### adaptation of mcmc tuning parameters #################

"""
    adapt_pospar_step!(δ, n, accinfo, η, update; adaptskip = 15, targetaccept=0.5)

Adapt the step size for updates, where tuning par has to be positive.
The adaptation is multiplication/divsion by exp(η(n)).
"""
function adapt_pospar_step(δ, n, accinfo, η, adaptskip; targetaccept=0.5)
    recent_mean = mean(accinfo[end-adaptskip+1:end])
    if recent_mean > targetaccept
        return δ * exp(η(n))
    else
        return δ * exp(-η(n))
    end
end

"""
    adaptpcnstep(n, accpcn, ρ, nshapes, η; adaptskip = 15, targetaccept=0.5)

Adjust pcN-parameter `ρ` adaptive every `adaptskip` steps.
For that `ρ` is first tranformed to (-∞, ∞), then updated, then mapped back to (0,1)
"""
function adaptpcnstep(ρ, n, accinfo, η, adaptskip; targetaccept=0.5)
    recentmean = mean(accinfo[end-adaptskip+1:end])
    if recentmean > targetaccept
        out =  logistic(logit(ρ) - η(n))
    else
        out =  logistic(logit(ρ) + η(n))
    end
    out
end

"""
    show_updates()

Convenience function to show the available implemted updates.
"""
function show_updates()
    println("Available implemented update steps are:")
    println(":innov")
    println(":sgd_mom")
    println(":mala_mom, :tmala_mom, :rmmala_mom, :rmrw_mom")
    println(":mala_pos, :tmala_pos, :rmmala_pos")
    println(":parameter, :matching")
end


text:
"""
    normalize_units(val)
On normal numbers, it is the identity, but on numbers or vectors of `Unitful.Quantity`,it converts the unit to an equivalent SI unit and strips the unit information.

# Examples
```julia-repl
julia> normalize_units(0.0u"°C")
273.15

```
```julia-repl
julia> normalize_units(273.15)
273.15
```
"""
normalize_units(x) =  Unitful.ustrip.(Unitful.upreferred.(x))

#upreferred, but the standard unit with just numbers is transformed to kg/mol
function mw_mul(x,mw::Unitful.Quantity)
    return upreferred(x*mw)
end

function mw_mul(x,mw)
    return 0.001*x*mw
end

function mw_div(x,mw::Unitful.Quantity)
    return upreferred(x/mw)
end

function mw_div(x,mw)
    return 1000.0*x/mw
end

"""
    convert_unit(from,to,val)
Converts an unit from the unit stored in `from` to the unit stored in `to`.

When both units are equal, it justs returns `val`.

If `val` itself is an `unit`, then it converts the from the unit in `val` to the unit in `to`.

# Examples

```julia-repl
julia> convert_unit(u"Pa",u"kPa",1000.0)
1.0
```

```julia-repl
julia> convert_unit(u"Pa",u"kPa",1u"atm")
4053//40
```
"""
function convert_unit end
function convert_unit(from::T,to::T,val::N) where {T,N}
    return val
end

function convert_unit(from::T1,to::T2,val::N) where {T1,T2,N<:Number}
    return Unitful.ustrip(Unitful.uconvert(to,val*from))
end

function convert_unit(from::T1,to::T2,val::N) where {T1,T2,N<:AbstractVector{<: Number}}
    return Unitful.ustrip.(Unitful.uconvert.(to,val*from))
end

function convert_unit(from::T1,to::T2,val::N) where {T1,T2,N<:AbstractVector{<:Unitful.Quantity}}
    return Unitful.ustrip.(Unitful.uconvert.(to,val))
end

is_real(x) = false
is_real(x::Real) = true
is_real(x::Bool) = false
is_real(x::AbstractVector{T} where T<:Real) = true
is_unitful(x) = false
is_unitful(x::Unitful.Quantity) = true
function is_unitful(x::T) where T <: AbstractArray{T2} where T2 <: Unitful.Quantity
    return true
end


text:
function Base.:+(a::GainObj, b::GainObj)
    new_minus = a.minus + b.minus
    new_plus = a.plus + b.plus
    new_minus_counter = a.minus_counter + b.minus_counter
    new_plus_counter = a.plus_counter + b.plus_counter
    new_inf_counter = a.inf_counter + b.inf_counter
    return GainObj(new_minus, new_plus, new_minus_counter, new_plus_counter, new_inf_counter)
end

function Base.copy(a::GainObj)
    return GainObj(a.minus, a.plus, a.minus_counter, a.plus_counter, a.inf_counter)
end

function init_gains(num_disc_var)
    gains_m = zeros(num_disc_var)
    gains_p = zeros(num_disc_var)
    gains_mc = zeros(Int64,num_disc_var)
    gains_pc = zeros(Int64,num_disc_var)
    gains_inf = zeros(Int64,num_disc_var)
    return GainObj(gains_m, gains_p, gains_mc, gains_pc, gains_inf)
end

"""
    update_gains!(tree::BnBTreeObj, parent::BnBNode, l_nd, r_nd)

Update the objective gains for the branch variable used for node
"""
function update_gains!(tree::BnBTreeObj, parent::BnBNode, l_nd, r_nd)
    gain_l = sigma_minus(tree.options, parent, l_nd, parent.solution[parent.var_idx])
    gain_r = sigma_plus( tree.options, parent, r_nd, parent.solution[parent.var_idx])
    idx = tree.var2disc_idx[parent.var_idx]

    gain = 0.0
    gain_c = 0
    if !isinf(gain_l)
        tree.obj_gain.minus[idx] = gain_l
        tree.obj_gain.minus_counter[idx] += 1
        gain += gain_l
        gain_c += 1
    else
        tree.obj_gain.inf_counter[idx] += 1
    end

    if !isinf(gain_r)
        tree.obj_gain.plus[idx] = gain_r
        tree.obj_gain.plus_counter[idx] += 1
        gain += gain_r
        gain_c += 1
    else
        tree.obj_gain.inf_counter[idx] += 1
    end

    if !isinf(gain_l) && !isinf(gain_r) && tree.obj_gain.inf_counter[idx] > 0
        tree.obj_gain.inf_counter[idx] -= 1
    end
    if gain_c == 0
        return 0
    else
        return gain/gain_c
    end
end

"""
    upd_gains_step!(tree, step_obj)

Update the gains using the step_obj if using StrongPseudoCost or PseudoCost
"""
function upd_gains_step!(tree, step_obj)
    branch_strat = tree.options.branch_strategy
    opts = tree.options
    if step_obj.upd_gains == :GainsToTree || (branch_strat == :StrongPseudoCost && step_obj.counter <= opts.strong_branching_nsteps)
        tree.obj_gain += step_obj.obj_gain
        tree.obj_gain.inf_counter[tree.obj_gain.inf_counter .< 0] .= 0
        if step_obj.counter == 1
            cum_counter = tree.obj_gain.minus_counter .+ tree.obj_gain.plus_counter
            strong_disc_vars = findall(cum_counter .> 0)
            step_obj.strong_disc_vars = strong_disc_vars
            # all other variables that haven't been checked get the median value of the others
            med_gain_m = median(tree.obj_gain.minus[strong_disc_vars])
            med_gain_p = median(tree.obj_gain.plus[strong_disc_vars])
            rest = filter(i->!(i in strong_disc_vars),1:tree.m.num_disc_var)
            tree.obj_gain.minus[rest] .= med_gain_m
            tree.obj_gain.plus[rest] .= med_gain_p
            tree.obj_gain.minus_counter[rest] .= 1
            tree.obj_gain.plus_counter[rest] .= 1
        end
    elseif step_obj.upd_gains == :GuessAndUpdate || branch_strat == :PseudoCost || (branch_strat == :StrongPseudoCost && step_obj.counter > opts.strong_branching_nsteps)
        upd_start = time()
        guess_gain_val = guess_gain(tree, step_obj)
        gain = update_gains!(tree, step_obj.node, step_obj.l_nd, step_obj.r_nd)
        step_obj.gain_gap = abs(guess_gain_val-gain)/abs(gain)
        if isnan(step_obj.gain_gap) # 0/0
            step_obj.gain_gap = 0.0
        end
        step_obj.upd_gains_time = time()-upd_start
    end
end

function guess_gain(tree, step_obj)
    i = step_obj.var_idx
    disc_i = tree.var2disc_idx[i]
    x = step_obj.node.solution
    g_minus, g_minus_c = tree.obj_gain.minus,tree.obj_gain.minus_counter
    g_plus, g_plus_c = tree.obj_gain.plus,tree.obj_gain.plus_counter
    g_minus_c += map(i -> (i == 0) && (i = 1), g_minus_c)
    g_plus_c += map(i -> (i == 0) && (i = 1), g_plus_c)
    mu = tree.options.gain_mu
    return score(f_minus(x[i])*g_minus[disc_i]/g_minus_c[disc_i],f_plus(x[i])*g_plus[disc_i]/g_plus_c[disc_i],mu)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

using JuliaFEM
using FEMBase
using Logging
Logging.configure(level=DEBUG)

type Coupling <: BoundaryProblem
    reference_node :: Element{Poi1}
end

function Coupling()
    reference_node = Element(Poi1, Int64[])
    return Coupling(reference_node)
end

d = 3.0
X = Dict(1 => [0.0, 0.0],
         2 => [d, 0.0],
         3 => [d, d],
         4 => [0.0, d],
         5 => [3*d, 1/2*d])
element1 = Element(Quad4, [1, 2, 3, 4])
update!(element1, "geometry", X)
update!(element1, "youngs modulus", 288.0)
update!(element1, "poissons ratio", 1/3)

problem = Problem(Elasticity, "test problem", 2)
problem.properties.formulation = :plane_strain
add_elements!(problem, [element1])

bc = Problem(Dirichlet, "fixed", 2, "displacement")
element2 = Element(Seg2, [4, 1])
update!(element2, "geometry", X)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
add_elements!(bc, [element2])

"""
    add_coupling_nodes!(problem, elements)

Add new coupling nodes into the problem. Nodes must be defined with Poi1
elements.
"""
function add_coupling_nodes!(problem, elements)
    for element in elements
        push!(problem.elements, element)
    end
    return nothing
end

"""
    add_reference_node!(problem, element)

Add new reference node into the problem. Node must be defined with a Poi1
element. If reference node is already defined, it will be replaced with new one.
"""
function add_reference_node!(problem, element)
    problem.properties.reference_node = element
    return nothing
end

function FEMBase.assemble_elements!(problem::Problem{Coupling},
                                    assembly::Assembly,
                                    elements::Vector{Element{Poi1}}, time)

    ref_node = problem.properties.reference_node
    if length(get_connectivity(ref_node)) == 0
        error("Reference node node defined. Define reference node using add_reference_node!")
    end
    ref_node_id = first(get_connectivity(ref_node))
    info("Reference node id = $ref_node_id")
    X_r = first(ref_node("geometry", time))
    info("Reference node geometry = $X_r")
    fe = zeros(2)
    for coupling_node in elements
        fill!(fe, 0.0)
        coupling_node_id = first(get_connectivity(coupling_node))
        X_n = first(coupling_node("geometry", time))
        info("Coupling node id = $coupling_node_id, geometry = $X_n")
        gdofs = get_gdofs(problem, coupling_node)
        info("gdofs = $gdofs")
        if haskey(ref_node, "point moment 3")
            M = ref_node("point moment 3")
            info("calculate point moment")
            # fe += ...
        end
        if haskey(ref_node, "point load 1")
            M = ref_node("point load 1")
            info("calculate point moment")
            # fe += ...
        end
        if haskey(ref_node, "point load 2")
            M = ref_node("point load 2")
            info("calculate point moment")
            # fe += ...
        end
        add!(assembly.f, gdofs, fe)
    end
end

coupling = Problem(Coupling, "test", 2, "displacement")
# "slave nodes"
element3 = Element(Poi1, [2])
element4 = Element(Poi1, [3])
update!([element3, element4], "geometry", X)
element5 = Element(Poi1, [5])
update!(element5, "geometry", X)
update!(element5, "point moment 3", 10.0/d)
add_coupling_nodes!(coupling, [element3, element4])
add_reference_node!(coupling, element5)

step = Analysis(Nonlinear)
#xdmf = Xdmf("example1_results"; overwrite=true)
#add_results_writer!(step, xdmf)
add_problems!(step, [problem, bc, coupling])
run!(step)
#close(xdmf.hdf)
time = 0.0
u = element1("displacement", time)
println("displacement: $u")


text:
module adapters_postgresql_bigserial_pk_test

using Test # @test
using Octo.Adapters.PostgreSQL # Repo Raw Schema.model DROP TABLE IF EXISTS

Repo.debug_sql()

include("options.jl")

Repo.connect(;
    adapter = Octo.Adapters.PostgreSQL,
    Options.for_postgresql...
)

struct Post
end
Schema.model(Post, table_name="posts")

Repo.execute([DROP TABLE IF EXISTS Post])

Repo.execute(Raw("""
CREATE TABLE posts (
  id BIGSERIAL PRIMARY KEY,
  body varchar,
  created_at timestamp DEFAULT current_timestamp
);
"""))

posts = from(Post, :posts)

Repo.insert!(Post, (body="hello",))

df = Repo.query(Post)
@test df[1].body == "hello"
@test typeof(df[1].id) === Int64 # 8 bytes large autoincrementing integer

end # module adapters_postgresql_bigserial_pk_test


text:
# Autogenerated wrapper script for nPth_jll for x86_64-apple-darwin
export libnpth

JLLWrappers.@generate_wrapper_header("nPth")
JLLWrappers.@declare_library_product(libnpth, "@rpath/libnpth.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnpth,
        "lib/libnpth.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using NeuralQuantum, Test
num_types = [Float32, Float64]

#=
@testset "NDMSymm" begin
    N = 4
    permutations = [[1,2,3,4],[2,3,4,1],[3,4,1,2],[4,1,2,3]]
    for T=num_types
        net  = NDMSymm(T,N,2,3, permutations)
        cnet = cached(net)
        v    = DoubleState(NAryState(T, 2, N))

        cder = grad_cache(cnet)
        vals = []; cvals = [];
        grads = [];
        for i=1:spacedimension(v)
            set_index!(v, i)
            push!(vals, net(v))
            push!(cvals, cnet(v))

            der  = ∇logψ(net,  v)
            ∇logψ!(cder, cnet, v)
            for f=fieldnames(typeof(der))
                ∇  = getfield(der,  f)
                c∇ = getfield(cder, f)
                push!(grads, ∇ ≈ c∇)
            end
        end
        @test vals ≈ cvals
        @test all(grads)
        @test NeuralQuantum.input_type(net) == real(T)
        @test NeuralQuantum.out_type(net) == Complex{real(T)}
        @test NeuralQuantum.input_shape(net) == (N, N)
        @test NeuralQuantum.is_analytic(net)
    end

end
=#


text:
using DataDrivenDiffEq
using DataDrivenDiffEq.Optimize
using ModelingToolkit
using LinearAlgebra
using SafeTestsets

@info "Loading OrdinaryDiffEq"
using OrdinaryDiffEq
using Test
@info "Finished loading packages"

const GROUP = get(ENV, "GROUP", "All")

@time begin
    if GROUP == "All" || GROUP == "DataDrivenDiffEq" || GROUP == "Standard"
        include("./basis.jl")
        include("./koopman.jl")
        include("./sindy.jl")
        include("./isindy.jl")
        include("./utils.jl")
        include("./optimize.jl")
    end

    # These are excluded right now, until the deps are figured out
    if GROUP == "Integration"
        @safetestset "Partial Lotka Volterra Discovery " begin include("./applications/partial_lotka_volterra.jl") end
    end
end


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(
    model::InfiniteModel,
    index::FiniteVariableIndex
    )::FiniteVariableRef
    return FiniteVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(
    model::InfiniteModel,
    object::VariableData{<:JuMP.ScalarVariable}
    )::FiniteVariableIndex
    return MOIUC.add_item(model.finite_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(
    model::InfiniteModel,
    ::Type{FiniteVariable}
    )::MOIUC.CleverDict{FiniteVariableIndex, VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}}
    return model.finite_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(
    vref::FiniteVariableRef
    )::MOIUC.CleverDict{FiniteVariableIndex, VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}}
    return JuMP.owner_model(vref).finite_vars
end

# Extend _data_object
function _data_object(
    vref::FiniteVariableRef
    )::VariableData{JuMP.ScalarVariable{Float64, Float64, Float64, Float64}}
  object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
  object === nothing && error("Invalid finite variable reference, cannot find " *
                        "corresponding variable in the model. This is likely " *
                        "caused by using the reference of a deleted variable.")
  return object
end

# Extend _core_variable_object
function _core_variable_object(
    vref::FiniteVariableRef
    )::JuMP.ScalarVariable{Float64, Float64, Float64, Float64}
    return _data_object(vref).variable
end

################################################################################
#                             DEFINTION METHODS
################################################################################
# Define _check_and_make_variable_ref (used by JuMP.add_variable)
function _check_and_make_variable_ref(
    model::InfiniteModel,
    v::JuMP.ScalarVariable{Float64, Float64, Float64, Float64},
    name::String
    )::FiniteVariableRef
    data_object = VariableData(v, name)
    vindex = _add_data_object(model, data_object)
    vref = FiniteVariableRef(model, vindex)
    return vref
end

# Define for scalar variables without strictly Float64 info
function _check_and_make_variable_ref(
    model::InfiniteModel,
    v::JuMP.ScalarVariable,
    name::String
    )::FiniteVariableRef
    new_v = JuMP.ScalarVariable(_make_float_info(v.info))
    return _check_and_make_variable_ref(model, new_v, name)
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Set info for finite variables
function _update_variable_info(
    vref::FiniteVariableRef,
    info::JuMP.VariableInfo
    )::Nothing
    _set_core_variable_object(vref, JuMP.ScalarVariable(info))
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::FiniteVariableRef)::Nothing
    # remove variable info constraints associated with vref
    _delete_info_constraints(vref)
    return
end


text:
module ExceptionRequest

export StatusError

import ..Layer, ..request
import ..HTTP
using ..Messages: iserror

"""
    request(ExceptionLayer, ::URI, ::Request, body) -> HTTP.Response

Throw a `StatusError` if the request returns an error response status.
"""
abstract type ExceptionLayer{Next <: Layer} <: Layer end
export ExceptionLayer

function request(::Type{ExceptionLayer{Next}}, a...; kw...) where Next

    res = request(Next, a...; kw...)

    if iserror(res)
        throw(StatusError(res.status, res.request.method, res.request.target, res))
    end

    return res
end

"""
    StatusError <: Exception

The `Response` has a `4xx`, `5xx` or unrecognised status code.

Fields:
 - `status::Int16`, the response status code.
 - `response` the [`HTTP.Response`](@ref)
"""
struct StatusError <: Exception
    status::Int16
    method::String
    target::String
    response::HTTP.Response
end

# for backwards compatibility
StatusError(status, response::HTTP.Response) = StatusError(status, response.request.method, response.request.target, response)

end # module ExceptionRequest


text:
using ImageMagick, ImageMetadata, ImageTransformations

using Random: bitrand
using Base.CoreLogging: SimpleLogger, with_logger

include("constructed_images.jl")
include("readremote.jl")
include("badimages.jl")

workdir = joinpath(tempdir(), "Images")
try
    rm(workdir, recursive=true)
catch
end


text:
# A module template for Koala supervised models

## INTRODUCTION

# This file briefly describes the low-level methods that must be
# implemented for each new Koala supervised learning algorithm. If
# this algorithm is called "SomeAlgorithm", then this implementation
# is conventionally contained in a module called `KoalaSomeAlgorithm`,
# eg, `KoalaTrees`.  The reader implementing a new algorithm in Koala
# may find the present file useful as a template.

# Normally regressor and classifier models with similar underlying
# algorithms should be placed in the same module.

# For each supervised learning algorithm one defines a supervised
# model type, namely a mutatable structure (a subtype of
# `Koala.SupervisedModel`) whose instances, hereafter called "models",
# store parameters controlling the outcome of training. A parameter
# controlling *how*, rather than *what*, the algorithm does is not a
# model parameter for present purposes. When a model parameter is
# changed during training, a machine wrapped around the model (which
# includes transformed versions of the learning data - see Koala.jl)
# does not need to be reconstructed. Parameters describing the
# transformations to be applied to input or target data (eg, whether
# to standardize target values) are *not* viewed as model parameters,
# as changing them generally necessitates wrapping the model in a new
# machine.

# An example of a model parameter is `min_patterns_split`; this is a
# field of `TreeRegressor` in `KoalaTrees` which controls the degree
# of pruning of decision trees. The number of trees in a random forest
# model is another model parameter, but a flag saying whether or not
# the trees should be trained in parallel is not a model
# parameter.

# Fitting a model to training data (see below) produces a "predictor"
# which is not part of the model. However, the *type* of the predictor
# must be declared implicitly when defining the model type, by an
# appropriately subtyping as shown below. For example, the predictor
# type for the `ConstantRegressor` model type (for simply predicting
# the target mean for any input pattern) is `Float64`.

# Some general guidelines:

# 1. Avoid giving models inner constructors, keeping in mind that
# these are *mutable* types. (Any desired invariants can and should be
# be enforced in the `clean!` method.)

# 2. With the exception of the `fit` method for transformers (see
# `TransformerTemplate.jl.`) avoid type and dimension checking in
# low-level methods. They should not be needed.

################################
# MODULE TEMPLATE STARTS BELOW #
################################

module KoalaSomeAlgorithm # eg KoalaRidge

# new:
export SomeSupervisedModelType # eg, RidgeRegressor

# needed in this module:
import Koala: Regressor, Classifier, softwarn
import DataFrames: AbstractDataFrame # the form of all untransformed input data

# to be extended (but not explicitly rexported):
import Koala: setup, fit, predict
import Koala: default_transformer_X, default_transformer_y, clean!
import Koala: transform, inverse_transform

## MODEL TYPE DEFINITIONS

mutable struct SomeSupervisedModelType <: Regressor{CorrespondingPredictorType}
    # eg ConstantRegressor <: Regressor{Float64}
    param1::Float64
    param2::Bool
    ...and so on...
end
# or
mutable struct SomeSupervisedModelType <: Classifier{CorrespondingPredictorType}
    param1::Float64
    param2::Bool
    ...and so on...
end

# lazy keywork constructor:
function SomeSupervisedModelType(; param1=default1, parmam2=defalut2, etc)
    model = SomeSupervisedModelType(param1, param2, etc)
    softwarn(clean!(model)) # only issues warning if `clean!` changes `model`
    return model
end

# The following optional method may change the fields of `model` to
# attempt to enforce natural invariants and avoid errors in later
# calls to `fit`. Any changes should be reported to the returned
# string (empty if the model parameters are fine). As a last resort, it
# may throw an exception. The high-level `fit!` method always calls
# this method on the model before doing anything else.
function clean!(model::SomeSupervisedModelType) -> message

## TRANSFORMERS

# For each model default `Transformer` objects must be specified with
# the following methods. The transformers are used to transform the
# input and output patterns into a form acceptable to the model's
# corresponding `fit` method. If the fall-backs
# `transformer_X=FeatureTruncater()` and
# `transformer_y=IdentityTransformer()` given in `Koala` suffice, then
# omit the two lines below; otherwise ammend appropriately:
default_transformer_X(model::SomeSupervisedModelType) -> transformer_X::Transformer
default_transformer_y(model::SomeSupervisedModelType) -> transformer_y::Transformer

# Notes:

# 1. Of course if `transformer_X` and `transformer_y` are of new
#   transformer types, then corresponding `fit` and `transform`
#   methods will need to be provided (see TransformerTemplate.jl). An
#   `inverse_transform` method is required only in the case of
#   `transformer_y`. At the very least, the scheme that encodes input
#   transformations should record the labels of the features of the
#   dataframe `X` given as input, so that only these are retained in
#   calls to `transform` on new data.

# 2. Any other "metadata" needed for reporting, in calls to the
#   supervised model's `fit` method defined below, should also be
#   stored in the scheme, `scheme_X`, output by `fit(transformer_X,
#   etc...)`.  The metadata should be extracted from `scheme_X` in
#   `setup` below and bundled in the `cache` part of that method's
#   return value. In this way `fit` will have access to the metadata
#   through its `cache` input.

# 3. To avoid data leakage, it is understood that transformers are
#   only fit on training data.

## TRAINING AND PREDICTION METHODS

# Any preliminary part of training that does *not* depend on the value
# of `model` (ie it's field values) can be placed in
# `setup`. Typically, `setup` places in to its `cache` return value
# the training input and target data needed for `fit`, together with
# some metadata. (It is assumed that the data received by `setup` is
# already transformed.)  Note that every time new training rows are
# given to the high-level method `fit!`, `setup` will be called anew.
# In iterative training algorithms, any such calculation which needn't
# be repeated when iterations are added should go in `setup`. All
# other training should go in `fit`. All results of `setup` needed for
# the rest of training must be returned as `cache` for passing to
# `fit`.
setup(model::SomeSupervisedModelType, Xt, yt,
      scheme_X, parallel, verbosity) -> cache

# The value of `predictor` returned below need only include the
# minimum data needed to furnish predictions on new patterns. In
# iterative training algorithms, data that must be saved for further
# iterations (eg, `Xt`, `yt`) should be stored and returned in
# `cache`. Store in `report` post-fit statistics, such as feature
# rankings or training (in sample) error.
#
# Calling `fit` with `add=true` for an iterative training algorithm
# should add iterations.
#
# IMPORTANT: Care should be taken that successive calls to `fit` with
# `add=false` return *distinct* `predictor` objects, in the sense that
# any change to one has no effect on the other. If this is not
# possible, then uncomment the following line:
# Koala.supports_ensembling(model::SomeSupervisedModelType) = false

# We allow `fit` to change the value of `cache` (but do not use the
# name "fit!" to distiguish the method from the high-level method with
# that name). You may design `fit` to take keyword arguments, `args`
# (to implement special training instructions) and if the high-level
# method `fit!` is called with these arguments on a corresponding
# `SupervisedMachine` object, then these will be passed along to the
# low-level `fit` method below.
function fit(model::SomeSupervisedModelType, cache, add, parallel, verbosity; args...)
    return predictor, report, cache
end

# The following uses `predictor` to predict *transformed* target
# values on new *transformed* input patterns.
predict(model::SomeSupervisedModelType,
        predictor::CorrespondingPredictorType,
        Xt, parallel, verbosity) -> transformed_predictions

# Note that `CorrespondingPredictorType` must be the type appearing in
# the type declaration of `SomeSupervisedModelType` above.

end # of module


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

kernelabstractions = Pkg.PackageSpec(path = root_directory)

BACKEND = get(ENV, "KERNELABSTRACTIONS_TEST_BACKEND", "all")
BUILDKITE = parse(Bool, get(ENV, "BUILDKITE", "false"))

@info "Develop..." BUILDKITE

Pkg.develop(kernelabstractions)
if !(VERSION < v"1.6-")
    if !BUILDKITE  || BACKEND == "ROCM"
        rockernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "ROCKernels"))
        Pkg.develop(rockernels)
    end

    if !BUILDKITE  || BACKEND == "CUDA"
        cudakernels = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "CUDAKernels"))
        Pkg.develop(cudakernels)
    end

    kernelgradients = Pkg.PackageSpec(path = joinpath(root_directory, "lib", "KernelGradients"))
    Pkg.develop(kernelgradients)
end
Pkg.build()
Pkg.precompile()


text:
function complementary(v::AbstractArray{Float64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Vector{Float64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Vector{Int}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::Int64, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end

function complementary(v::AbstractArray{Int64}, n::Int64)
  ns = collect(1:n)
  setdiff(1:n, v)
end


text:
@testset "Testing Transforms..." begin

    for inplace in [true, false]
        @testset "ScalePointCloud inplace=$(inplace)" begin
            p = rand(Float32, 3, 8, 2)
            t1 = ScalePointCloud(2.0; inplace = inplace)
            t2 = ScalePointCloud(0.5; inplace = inplace)
            pc1 = PointCloud(p)
            pc2 = t2(t1(pc1))
            if inplace
                @test pc1 == pc2
            else
                @test pc1 != pc2
            end
            @test all(isapprox.(p, pc2.points, rtol = 1e-5, atol = 1e-5))
        end
    end

    for inplace in [true, false]
        @testset "RotatePointCloud inplace=$(inplace)" begin
            p = rand(Float32, 3, 8, 2)
            rotmat = 2 .* one(rand(Float32, 3, 3))
            rotmat_inv = inv(rotmat)

            t1 = RotatePointCloud(rotmat; inplace = inplace)
            t2 = RotatePointCloud(rotmat_inv; inplace = inplace)
            pc1 = PointCloud(p)
            pc2 = t2(t1(pc1))
            if inplace
                @test pc1 == pc2
            else
                @test pc1 != pc2
            end
            @test all(isapprox.(p, pc2.points, rtol = 1e-5, atol = 1e-5))
        end
    end

    for inplace in [true, false]
        @testset "ReAlignPointCloud inplace=$(inplace)" begin
            p1 = rand(Float32, 3, 8, 2)
            p2 = rand(Float32, 3, 8, 1)
            src = PointCloud(p1)
            tgt = PointCloud(p2)
            t = ReAlignPointCloud(tgt; inplace = inplace)
            src_1 = t(src)
            if inplace
                @test src == src_1
            else
                @test src != src_1
            end
            # Transformed PointCloud should be inside the bounding box defined by `tgt` PointCloud
            @test all(
                reshape(maximum(tgt.points, dims = 2), (:, 1)) .>=
                src_1.points .>=
                reshape(minimum(tgt.points, dims = 2), (:, 1)),
            )
        end
    end

    for inplace in [true, false]
        @testset "NormalizePointCloud inplace=$(inplace)" begin
            p = rand(Float32, 3, 8, 2)
            t = NormalizePointCloud(; inplace = inplace)
            pc1 = PointCloud(p)
            pc2 = t(pc1)
            if inplace
                @test pc1 == pc2
            else
                @test pc1 != pc2
            end
            @test all(isapprox.(mean(pc2.points; dims = 2), 0, rtol = 1e-5, atol = 1e-5))
            @test all(isapprox.(std(pc2.points; dims = 2), 1, rtol = 1e-5, atol = 1e-5))
        end
    end

    @testset "Chain" begin
        p = rand(Float32, 3, 8, 2)
        t = Chain(ScalePointCloud(0.5), RotatePointCloud(rand(3, 3)), NormalizePointCloud())
        pc1 = PointCloud(p)
        pc1 = t(pc1)
        @test pc1.points isa Array
        @test all(isapprox.(mean(pc1.points; dims = 2), 0, rtol = 1e-5, atol = 1e-5))
        @test all(isapprox.(std(pc1.points; dims = 2), 1, rtol = 1e-5, atol = 1e-4))
    end

    _mesh = load_trimesh(
        joinpath(@__DIR__, "../assets/teapot.obj"),
        joinpath(@__DIR__, "../assets/sphere.obj"),
    )

    for inplace in [true, false]
        @testset "ScaleTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            t1 = ScaleTriMesh(2.0; inplace = inplace)
            t2 = ScaleTriMesh(0.5; inplace = inplace)
            m2 = t2(t1(m))

            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for inplace in [true, false]
        @testset "RotateTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            rotmat = 2 .* one(rand(Float32, 3, 3))
            rotmat_inv = inv(rotmat)
            t1 = RotateTriMesh(rotmat; inplace = inplace)
            t2 = RotateTriMesh(rotmat_inv; inplace = inplace)
            m2 = t2(t1(m))
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for inplace in [true, false]
        @testset "ReAlignTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            src = Flux3D.scale(m, 2.0)
            tgt = Flux3D.scale(m, 1.0)
            t = ReAlignTriMesh(tgt, 1; inplace = inplace)
            src_1 = t(src)
            if inplace
                @test src === src_1
            else
                @test src !== src_1
            end
            src_1 = t(src)
            @test all(
                maximum(get_verts_list(_mesh)[1]; dims = 2) .>=
                get_verts_list(src_1)[1] .>=
                minimum(get_verts_list(_mesh)[1]; dims = 2),
            )
            @test all(
                maximum(get_verts_list(_mesh)[1]; dims = 2) .>=
                get_verts_list(src_1)[2] .>=
                minimum(get_verts_list(_mesh)[1]; dims = 2),
            )
        end
    end

    for inplace in [true, false]
        @testset "NormalizeTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            t = NormalizeTriMesh(; inplace = inplace)
            m2 = t(m)
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    mean(get_verts_list(m2)[1]; dims = 2),
                    0.0,
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    std(get_verts_list(m2)[1]; dims = 2),
                    1.0,
                    rtol = 1e-4,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    mean(get_verts_list(m2)[2]; dims = 2),
                    0.0,
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
            @test all(
                isapprox.(
                    std(get_verts_list(m2)[2]; dims = 2),
                    1.0,
                    rtol = 1e-4,
                    atol = 1e-5,
                ),
            )
        end
    end

    for inplace in [true, false]
        @testset "TranslateTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            t1 = TranslateTriMesh(1.0; inplace = inplace)
            t2 = TranslateTriMesh(-1.0; inplace = inplace)
            m2 = t2(t1(m))
            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    for inplace in [true, false]
        @testset "OffsetTriMesh inplace=$(inplace)" begin
            m = deepcopy(_mesh)
            _offset = ones(size(get_verts_packed(m)))
            t1 = OffsetTriMesh(_offset; inplace = inplace)
            t2 = OffsetTriMesh(_offset .* -1; inplace = inplace)
            m2 = t2(t1(m))

            if inplace
                @test m2 === m
            else
                @test m2 !== m
            end
            @test all(
                isapprox.(
                    get_verts_packed(_mesh),
                    get_verts_packed(m2),
                    rtol = 1e-5,
                    atol = 1e-5,
                ),
            )
        end
    end

    _voxels = zeros(Float32, 32, 32, 32, 2)
    _voxels[1:15, 2:10, 18:32, :] .= 1

    _v = VoxelGrid(_voxels)
    _m = load_trimesh(
        joinpath(@__DIR__, "../assets/teapot.obj"),
        joinpath(@__DIR__, "../assets/sphere.obj"),
    )
    _p = PointCloud(sample_points(_m, 1024))

    res = 28
    points = 512
    thresh = 0.9

    @testset "TriMeshToVoxelGrid" begin
        m = deepcopy(_m)
        t = TriMeshToVoxelGrid(res)
        v = t(m)
        @test v isa VoxelGrid{Float32}
        @test size(v.voxels) == (res, res, res, 2)
    end

    @testset "PointCloudToVoxelGrid" begin
        p = deepcopy(_p)
        t = PointCloudToVoxelGrid(res)
        v = t(p)
        @test v isa VoxelGrid{Float32}
        @test size(v.voxels) == (res, res, res, 2)
    end

    for algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets]
        @testset "VoxelGridToTriMesh algo=$(algo)" begin
            v = deepcopy(_v)
            t = VoxelGridToTriMesh(thresh = thresh, algo = algo)
            m = t(v)
            @test m isa TriMesh{Float32,UInt32,Array}
        end
    end

    for algo in [:Exact, :MarchingCubes, :MarchingTetrahedra, :NaiveSurfaceNets]
        @testset "PointCloudToTriMesh algo=$(algo)" begin
            p = deepcopy(_p)
            t = PointCloudToTriMesh(res, algo = algo)
            m = t(p)
            @test m isa TriMesh{Float32,UInt32,Array}
        end
    end

    @testset "TriMeshToPointCloud" begin
        m = deepcopy(_m)
        t = TriMeshToPointCloud(points)
        p = t(m)
        @test p isa PointCloud{Float32}
        @test size(p.points) == (3, points, 2)
    end

    for algo in [:Exact, :MarchingCubes, :MarchingTetrahedra]#, :NaiveSurfaceNets]
        @testset "VoxelGridToPointCloud algo=$(algo)" begin
            v = deepcopy(_v)
            t = VoxelGridToPointCloud(points, thresh = thresh, algo = algo)
            p = t(v)
            @test p isa PointCloud{Float32}
            @test size(p.points) == (3, points, 2)
        end
    end
end # transforms


text:
#
# This file is part of module
#   GenerateMediumDict (ModiaMedia/dict/GenerateMediumDict.jl)
#

const Modelica_Media_IdealGases_Common_DataRecord = ModiaMedia.SingleGasNasaData
const singleGasesData = Dict{AbstractString, ModiaMedia.SingleGasNasaData}()
const fluidData       = Dict{AbstractString, ModiaMedia.IdealGasFluidConstants}()

# Include data from Modelica.Media.IdealGases.Common.SingleGasesData.jl into dict
include(joinpath(dirname(@__FILE__), "Modelica.Media.IdealGases.Common.SingleGasesData.jl"))

# Include data from Modelica.Media.IdealGases.Common.FluidData.jl into dict
include(joinpath(dirname(@__FILE__), "Modelica.Media.IdealGases.Common.FluidData.jl"))

# Construct Medium
"""
    medium = SimpleGasNasa(name, fluidConstants, data)

Return an ideal gas medium of a single substance, or an ideal gas medium
consisting of several substances where the mass fractions are fixed.
Density is a function of T and p. All other quantities are solely a function of T.

# Sources for model and literature

Original Data: *Computer program for calculation of complex chemical equilibrium compositions and applications*.
Part 1: Analysis Document ID: 19950013764 N (95N20180)
File Series: NASA Technical Reports Report
Number: NASA-RP-1311 E-8017 NAS 1.61:1311
Authors: Gordon, Sanford (NASA Lewis Research Center) Mcbride, Bonnie J. (NASA Lewis Research Center)
Published: Oct 01, 1994.

The data to computpe the specific enthalpy is from
McBride B.J., Zehe M.J., and Gordon S. (2002):
*NASA Glenn Coefficients for Calculating Thermodynamic Properties of Individual Species*.
NASA report TP-2002-211556.
https://www.grc.nasa.gov/WWW/CEAWeb/TP-2002-211556.pdf

Known limits of validity: The data is valid for temperatures between 200K and 6000K.
A few of the data sets for monatomic gases have a discontinuous 1st derivative at 1000K,
but this never caused problems so far.

This model has been adapted from the Modelica.Media.IdealGases.Common.SingleGasNasa
Modelica package, which in turn was copied from the ThermoFluid library and
adapted to Modelica.
"""
function SingleGasNasa(name::AbstractString, fluidConstants::ModiaMedia.IdealGasFluidConstants, data::ModiaMedia.SingleGasNasaData)
    ModiaMedia.SingleGasNasa(mediumName     = name,
                             fluidConstants = fluidConstants,
                             fluidLimits    = ModiaMedia.FluidLimits(TMIN=200.0, TMAX=6000.0),
                             data           = data)
end

# Construct SingleGasNasa objects and store them in the medium dict
function storeSingleGasNasaMedium!(mediumDict)
    global fluidData
    global singleGasesData
    for (name,fluidConstants) in fluidData
        dict[name] = SingleGasNasa(name, fluidConstants, singleGasesData[name])
    end
end

# Manually insert air into the dictionary according to model from Modelica.Media.Air.MoistAir; should probably fix somehow
dict["Air"] = SingleGasNasa("Air", fluidData["N2"], singleGasesData["Air"])

storeSingleGasNasaMedium!(dict)


text:
using NonlinearTestProblems
using Test

const SQUARE_PROBLEMS = [F_NWp281(), Rosenbrock(), PowellSingular(), PowellBadlyScaled(),
                         HelicalValley()]

const RECTANGULAR_PROBLEMS = [Beale()]

@testset "basic consistency checks for test functions." begin
    for f in vcat(SQUARE_PROBLEMS, RECTANGULAR_PROBLEMS)
        x0 = root(f)
        @test length(x0) == domain_dimension(f)
        @test all(lower_bounds(f) .≤ x0 .≤ upper_bounds(f))
        fx0 = @inferred f(x0)
        @test length(fx0) == range_dimension(f)
        @test sum(abs2, fx0) ≤ eps()

        x1 = starting_point(f)
        @test length(x1) == domain_dimension(f)
        @test sum(abs2, f(x1)) > eps()

        @test starting_point(f, 0) ≈ x1
        x2 = starting_point(f, 1)
        @test sum(abs2, x2 .- x0) ≥ sum(abs2, x1 .- x0)
    end
end


text:
module Hydrozoa

using CUDAnative

include("device.jl")

using MPI
include("libmp.jl")

end # module


text:
@deprecate einsumexp(args...) loop_einsum(args...)
@deprecate einsumexp!(args...) loop_einsum!(args...)

@deprecate einsumopt einsum

# deprecated manually
function Iterators.flatten(c::Union{NestedEinsum,EinCode})
    @warn "`Iterators.flatten(eincode)` has been deprecated, use `OMEinsum.flatten` instead."
    flatten(c)
end

@deprecate dynamic_einsum(ixs, xs, iy; size_info=nothing) einsum(DynamicEinCode(ixs, iy), xs; size_info=size_info)
@deprecate dynamic_einsum(code::EinCode, xs; size_info=nothing) code(xs...; size_info=size_info)
@deprecate dynamic_einsum(code::NestedEinsum, xs; size_info=nothing) code(xs...; size_info=size_info)


text:
function show_trees(bbl::BlackBoxLearner)
    for lnr in bbl.learners
        IAI.show_in_browser(lnr)
    end
    return
end

show_trees(gm::GlobalModel) = show_trees.(gm.bbls)


text:
include("corey.jl")
include("grid.jl")
include("theis.jl")


text:
cd(@__DIR__)

using Pkg

using Test, TestSetExtensions, SafeTestsets
using SearchLight

@testset ExtendedTestSet "SearchLight tests" begin
  @includetests ARGS
end


text:
export vrr,hrr,chrr,cvrr,all_twoe_ints_chrr, vrr_autogen
# HGP - A (hopefully) fast, simple implementation of Head-Gordon, Pople's [^1]
# ERI recurrence relations.
#
# References
# [^1]: A method for two-electron Gaussian integral and integral derivative
#       evaluation using recurrence relations. Martin Head-Gordon and John
#       A. Pople. JCP, 89 (9), 5777, 1988.
# [^2]: Molecular Integrals Over Gaussian Basis Functions. Peter M. W. Gill. Adv.
#       Q. Chem., 25, 141 (1994).
# [^3]: The Prism Algorithm for Two-Electron Integrals. Peter M. W. Gill and John
#       A. Pople. IJQC, 40, 753 (1991).

"""
    all_twoe_ints_chrr(bfs)

Make multiple calls to `chrr` to form all required two-electron
integrals for a molecular basis set `bfs`.

Returns a 1d array of these integrals [ijkl] with i<j, k<l, and ij<kl.
"""
function all_twoe_ints_chrr(bfs)
    nbf = length(bfs)
    nints = length(collect(MolecularIntegrals.iiterator(nbf)))
    fetcher = eri_fetcher(bfs)
    ints = zeros(Float64,nints)
    Threads.@threads for (ishell,jshell,kshell,lshell) in collect(keys(fetcher))
        hrrs = chrr(bfs.shells[ishell],bfs.shells[jshell],bfs.shells[kshell],bfs.shells[lshell])
        for (ijkl,hi,hj,hk,hl) in fetcher[ishell,jshell,kshell,lshell]
            ints[ijkl] = hrrs[hi,hj,hk,hl]
        end
    end
    return ints
end

"""
    cvrr (ash,bsh,csh,dsh)

Compute and contract the vertical recurrence relations
between shells ash,bsh,csh,dsh.

This routine takes advantage of the fact that the optimal time for
contracting integrals is after the VRR step.

This function returns a nao[ash.L+bsh.L] x nao[csh.L+dsh.L]
with the contracted coefficients (a0|b0). It calls the
primitive function vrr multiple times using different
exponents and contraction coefficients.
"""
function cvrr(ash::Shell,bsh::Shell,csh::Shell,dsh::Shell)
    amax,cmax = ash.L+bsh.L,csh.L+dsh.L
    cvrrs = zeros(Float64,nao[amax],nao[cmax])
    mmax=amax+cmax+1
    vrrs = zeros(Float64,mmax,nao[amax],nao[cmax])

    A,B,C,D = ash.xyz,bsh.xyz,csh.xyz,dsh.xyz
    for (dexpn,dcoef) in zip(dsh.expns,dsh.coefs), (cexpn,ccoef) in zip(csh.expns,csh.coefs)
        for (bexpn,bcoef) in zip(bsh.expns,bsh.coefs), (aexpn,acoef) in zip(ash.expns,ash.coefs)
            #fill!(vrrs,0) # I don't know if this is really necessary
            vrr!(vrrs, amax,cmax, aexpn,bexpn,cexpn,dexpn,A,B,C,D)
            #vrr_turbo!(vrrs, amax,cmax, aexpn,bexpn,cexpn,dexpn,A,B,C,D)
            for j in 1:nao[cmax], i in 1:nao[amax]
                cvrrs[i,j] += (acoef*bcoef*ccoef*dcoef)*vrrs[1,i,j]
            end
        end
    end
    return cvrrs
end

"""
    chrr (ash,bsh,csh,dsh)

Compute the contracted horizontal recurrence relations
between shells ash,bsh,csh,dsh.

This routine makes calls `cvrr` to form the input VRRs,
and then applies the HRRs to form the full set of two
electron integrals for the four shells ash,bsh,csh,dsh.

The routine returns a 4d array where each dimensional
ranges over the number of orbitals in the shell.
"""
function chrr(ash::Shell,bsh::Shell,csh::Shell,dsh::Shell)
    ashell,bshell,cshell,dshell = ash.L,bsh.L,csh.L,dsh.L
    A,B,C,D = ash.xyz,bsh.xyz,csh.xyz,dsh.xyz
    vrrs = cvrr(ash,bsh,csh,dsh)
    hrrs = zeros(Float64,nao[ashell+bshell],nao[bshell],nao[cshell+dshell],nao[dshell])
    for c in 1:nao[cshell+dshell], a in 1:nao[ashell+bshell]
        hrrs[a,1,c,1] = vrrs[a,c]
    end

    # First build (ab,c0) from (a0,c0)
    for c in 1:nao[cshell+dshell], bplus in 2:nao[bshell]
        j = shift_direction[bplus]
        b = shift_index[bplus,j]
        bs = shell_number[b]
        for a in 1:nao[ashell+bshell-bs-1]
            aplus = shift_index_plus[a,j]
            hrrs[a,bplus,c,1] = hrrs[aplus,b,c,1] + (A[j]-B[j])*hrrs[a,b,c,1]
        end
    end

     # now build (ab,cd) from (ab,c0)
    for dplus in 2:nao[dshell]
       j = shift_direction[dplus]
       d = shift_index[dplus,j]
       ds = shell_number[d]
       for c in 1:nao[cshell+dshell-ds-1]
        cplus = shift_index_plus[c,j]
           for b in 1:nao[bshell],a in 1:nao[ashell]
                hrrs[a,b,c,dplus] = hrrs[a,b,cplus,d] +(C[j]-D[j])*hrrs[a,b,c,d]
            end
        end
    end
    return hrrs
end

"vrr!(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns an `n`x`m` array, where `n` is the number
of aos in the `a+b` shell, and `m` is the number of aos in the
`c+d` shell.

Equations refer to
[A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations](https://doi.org/10.1063/1.455553).
Martin Head-Gordon and John A. Pople. JCP, 89 (9), 5777, 1988."
function vrr!(vrrs,amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1

	# Prefactors arising from multiplying Gaussians together
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze # "ooz = one over zeta, etc."
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*oortze*exp(
		-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12 using an array of calls to
	# the incomplete gamma function. This is ~40% of the time here,
	# but can be replaced with table lookup or interpolation. The
	# simple version is included here.
    # First generate (1,1,m) using eq 12
    Tcrit=20.0 # Most code uses a much higher Tcrit (117)
    #boys_array = boys_array_Fgamma(mmax,T)
    for m in 1:mmax
        vrrs[m,1,1] = KabKcd_rtze*Fgamma(m-1,T)
    end

    if (mmax == 1) return nothing end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        if aminus > 0
            a_i = 0.5*ooz*ao2m[a][i]
            for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1] +
					a_i*(vrrs[m,aminus,1]-eta*ooze*vrrs[m+1,aminus,1])
            end
        else
            for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1]
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        cminus = shift_index[c,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            aminus = shift_index[a,i]
            if cminus > 0
                c_i = 0.5*ooe*ao2m[c][i]
                if aminus > 0
                    a_i = 0.5*ooze*ao2m[a][i]
                    for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            a_i*vrrs[m+1,aminus,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                else
                    for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                end
            elseif aminus > 0
                a_i = 0.5*ooze*ao2m[a][i]
                for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                        a_i*vrrs[m+1,aminus,c]
                end
            else
                for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c]
                end
            end
        end
    end
    return nothing
end

"vrr_simd!(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`amax` and `cmax` are related to the sum of the shell angular
momenta for the `a+b`, and `c+d` shells, respectively.

The function returns an `n`x`m` array, where `n` is the number
of aos in the `a+b` shell, and `m` is the number of aos in the
`c+d` shell.

Equations refer to
[A method for two-electron Gaussian integral and integral derivative evaluation using recurrence relations](https://doi.org/10.1063/1.455553).
Martin Head-Gordon and John A. Pople. JCP, 89 (9), 5777, 1988."
function vrr_simd!(vrrs,amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1

	# Prefactors arising from multiplying Gaussians together
    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze # "ooz = one over zeta, etc."
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*oortze*exp(
		-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)
    PA = P-A
    WP = W-P
    QC = Q-C
    WQ = W-Q

    # HGP equation 6, with b=d=0:
    #   [a+1,c]m = (Pi-Ai)[a,c]m + (Wi-Pi)[a,c]m+1
    #        + a_i/2zeta ([a-1,c]m - eta/zeta+eta[a-1,c]m+1)        # eq 6a
    #        + ci/2(zeta+eta)[a,c-1]m+1

    # First generate (1,1,m) using eq 12 using an array of calls to
	# the incomplete gamma function. This is ~40% of the time here,
	# but can be replaced with table lookup or interpolation. The
	# simple version is included here.
    # First generate (1,1,m) using eq 12
    Tcrit=20.0 # Most code uses a much higher Tcrit (117)
    #boys_array = boys_array_Fgamma(mmax,T)
    for m in 1:mmax
        vrrs[m,1,1] = KabKcd_rtze*Fgamma(m-1,T)
    end

    if (mmax == 1) return nothing end

    # Generate (A,1,m)
    # Eq 6a, with c=0 also:
    #   [a+1,0]m = (Pi-Ai)[a,1]m + (Wi-Pi)[a,1]m+1
    #        + a_i/2zeta ([a-1,0]m - eta/zeta+eta[a-1,0]m+1)        # eq 6b

    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        if aminus > 0
            a_i = 0.5*ooz*ao2m[a][i]
            @inbounds @simd for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1] +
					a_i*(vrrs[m,aminus,1]-eta*ooze*vrrs[m+1,aminus,1])
            end
        else
            @inbounds @simd for m in 1:lim
                vrrs[m,aplus,1] = PA[i]*vrrs[m,a,1] + WP[i]*vrrs[m+1,a,1]
            end
        end
    end

    # Now build (A,C,m)
    # The c-based version of 6a is:
    #   [a,c+1]m = (Qj-Bi)[a,c]m + (Wj-Qj)[a,c]m+1
    #       + c_j/2eta ([a,c-1]m - zeta/zeta+eta[a,c-1]m+1)         # eq 6d
    #       + a_j/2(zeta+eta)[a-1,c]m+1
    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        cminus = shift_index[c,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            aminus = shift_index[a,i]
            if cminus > 0
                c_i = 0.5*ooe*ao2m[c][i]
                if aminus > 0
                    a_i = 0.5*ooze*ao2m[a][i]
                    @inbounds @simd for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            a_i*vrrs[m+1,aminus,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                else
                    @inbounds @simd for m in 1:lim
                        vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                            c_i*(vrrs[m,a,cminus]-zeta*ooze*vrrs[m+1,a,cminus])
                    end
                end
            elseif aminus > 0
                a_i = 0.5*ooze*ao2m[a][i]
                @inbounds @simd for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c] +
                        a_i*vrrs[m+1,aminus,c]
                end
            else
                @inbounds @simd for m in 1:lim
                    vrrs[m,a,cplus] = QC[i]*vrrs[m,a,c]+WQ[i]*vrrs[m+1,a,c]
                end
            end
        end
    end
    return nothing
end

function vrr(amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    mmax=amax+cmax+1
    vrrs = zeros(Float64,mmax,nao[amax],nao[cmax])
    vrr!(vrrs, amax,cmax, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    return vrrs
end

"""
hrr(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's horizontal recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.
`ashell`, `bshell`, `cshell` and `dshell` are the shell angular
momenta for the `a`, `b`, `c`, and `d` shells, respectively.

The function returns a (k,l,m,n)-dimensional array, where the
dimensions correspond to the number of aos in the a,b,c,d shells.
"""
function hrr(ashell,bshell,cshell,dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    # Get the relevant vrr terms.
    vrrs = vrr(ashell+bshell,cshell+dshell, aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    hrrs = zeros(Float64,nao[ashell+bshell],nao[bshell],nao[cshell+dshell],nao[dshell])
    for j in 1:nao[cshell+dshell], i in 1:nao[ashell+bshell]
        hrrs[i,1,j,1] = vrrs[1,i,j]
    end

     # First build (ab,c0) from (a0,c0)
    for c in 1:nao[cshell+dshell], bplus in 2:nao[bshell]
        j = shift_direction[bplus]
        b = shift_index[bplus,j]
        bs = shell_number[b]
        for a in 1:nao[ashell+bshell-bs-1]
            aplus = shift_index_plus[a,j]
            hrrs[a,bplus,c,1] = hrrs[aplus,b,c,1] + (A[j]-B[j])*hrrs[a,b,c,1]
        end
     end

     # now build (ab,cd) from (ab,c0)
    for dplus in 2:nao[dshell]
        j = shift_direction[dplus]
        d = shift_index[dplus,j]
        ds = shell_number[d]
        for c in 1:nao[cshell+dshell-ds-1]
            cplus = shift_index_plus[c,j]
            for b in 1:nao[bshell], a in 1:nao[ashell]
                hrrs[a,b,c,dplus] = hrrs[a,b,cplus,d] +(C[j]-D[j])*hrrs[a,b,c,d]
            end
        end
    end
    return hrrs
end

"""
vrr_autogen(amax,cmax)

Autogenerate the unbranched, unrolled vrr routine between shells
`amax` and `cmax`.
"""
function vrr_autogen(amax,cmax)
    mmax = amax+cmax+1
    sa = llabel[amax]
    sc = llabel[cmax]
    indent = "    "
    asize = nao[amax]
    csize = nao[cmax]
    lines = []
    for m in 1:mmax
        push!(lines,"$(indent)vrrs[1,1,$m] = KabKcd_rtze*Fgamma($(m-1),T)")
    end
    slines = join(lines,"\n")

    palines = amax>0 ? "$(indent)PA = P-A\n$(indent)WP = W-P" : ""
    qclines = cmax>0 ? "$(indent)QC = Q-C\n$(indent)WQ = W-Q" : ""

    lines = []
    for aplus in 2:nao[amax]
        ashell = shell_number[aplus]
        i = shift_direction[aplus]
        a = shift_index[aplus,i]
        lim = mmax-ashell
        aminus = shift_index[a,i]
        for m in 1:lim
            line = ["$(indent)vrrs[$aplus,1,$m] = PA[$i]*vrrs[$a,1,$m] + WP[$i]*vrrs[$a,1,$(m+1)]"]
            if aminus > 0
                a_i = ao2m[a][i]
                #push!(line," +\n$indent$indent") # newline + indent
                push!(line," + ") # append to old line
                push!(line,"$a_i*0.5*ooz*(vrrs[$aminus,1,$m]-eta*ooze*vrrs[$aminus,1,$(m+1)])")
            end
            push!(lines,join(line,""))
        end
    end

    for cplus in 2:nao[cmax]
        cshell = shell_number[cplus]
        i = shift_direction[cplus]
        c = shift_index[cplus,i]
        for a in 1:nao[amax]
            ashell = shell_number[a]
            lim = mmax-cshell-ashell
            cminus = shift_index[c,i]
            aminus = shift_index[a,i]
            for m in 1:lim
                line = ["$(indent)vrrs[$a,$cplus,$m] = QC[$i]*vrrs[$a,$c,$m] + WQ[$i]*vrrs[$a,$c,$(m+1)]"]
                if cminus > 0
                    c_i = ao2m[c][i]
                    #push!(line," +\n$indent$indent") # newline + indent
                    push!(line," + ") # append to old line
                    push!(line,"$c_i*0.5*ooe*(vrrs[$a,$cminus,$m]-zeta*ooze*vrrs[$a,$cminus,$(m+1)])")
                end
                if aminus > 0
                    a_i = ao2m[a][i]
                    #push!(line," +\n$indent$indent") # newline + indent
                    push!(line," + ") # append to old line
                    push!(line,"$a_i*0.5*ooze*vrrs[$aminus,$c,$(m+1)]")
                end
                push!(lines,join(line,""))
            end
        end
    end
    glines = join(lines,"\n")

    function_template = """
"
vrr_$sa$sc(aexpn,bexpn,cexpn,dexpn, A,B,C,D)

Use Head-Gordon/Pople's vertical recurrence relations to compute
an array of two-electron integrals.

`A`, `B`, `C`, `D` are the centers of four Gaussian functions.
`aexpn`, `bexpn`, `cexpn`, `dexpn` are their exponents.

This is an auto generated routine specific to when the A shell
has $sa-type angular momentum, and the C shell has $sc-type
angular momentum.
"
function vrr_$sa$sc(aexpn,bexpn,cexpn,dexpn, A,B,C,D)
    vrrs = zeros(Float64,$asize,$csize,$mmax)

    zeta,eta = aexpn+bexpn,cexpn+dexpn
    ze = zeta+eta
    ooz,ooe,ooze = 1/zeta,1/eta,1/ze
    oortze = sqrt(ooze)
    P = (aexpn*A + bexpn*B)*ooz
    Q = (cexpn*C + dexpn*D)*ooe
    W = (zeta*P + eta*Q)*ooze
    rab2 = dist2(A-B)
    rcd2 = dist2(C-D)
    rpq2 = dist2(P-Q)
    T = zeta*eta*rpq2*ooze
    KabKcd_rtze = 2pi*pi*sqrt(pi)*ooze*exp(-aexpn*bexpn*rab2*ooz-cexpn*dexpn*rcd2*ooe)*oortze
$palines
$qclines
$slines
$glines
    return vrrs[:,:,1]
end
"""

    return function_template
end


text:
function LSlider(parent::Scene; bbox = nothing, kwargs...)

    attrs = merge!(Attributes(kwargs), default_attributes(LSlider))

    decorations = Dict{Symbol, Any}()

    @extract attrs (
        halign, valign, linewidth, buttonradius, horizontal,
        startvalue, value, color_active, color_active_dimmed, color_inactive,
        buttonstrokewidth, buttoncolor_inactive
    )

    sliderrange = attrs.range

    protrusions = Node(GridLayoutBase.RectSides{Float32}(0, 0, 0, 0))
    layoutobservables = LayoutObservables(LSlider, attrs.width, attrs.height,
        halign, valign, attrs.alignmode; suggestedbbox = bbox, protrusions = protrusions)

    onany(buttonradius, horizontal, buttonstrokewidth) do br, hori, bstrw
        layoutobservables.autosize[] = if hori
            (nothing, 2 * (br + bstrw))
        else
            (2 * (br + bstrw), nothing)
        end
    end

    subarea = lift(layoutobservables.computedbbox) do bbox
        IRect2D_rounded(bbox)
    end
    subscene = Scene(parent, subarea, camera=campixel!)

    sliderbox = lift(bb -> Rect{2, Float32}(zeros(eltype(bb.origin), 2), bb.widths), layoutobservables.computedbbox)

    endpoints = lift(sliderbox, horizontal) do bb, horizontal

        if horizontal
            y = bottom(bb) + height(bb) / 2
            [Point2f0(left(bb), y),
            Point2f0(right(bb), y)]
        else
            x = left(bb) + width(bb) / 2
            [Point2f0(x, bottom(bb)),
            Point2f0(x, top(bb))]
        end
    end

    # this is the index of the selected value in the slider's range
    # selected_index = Node(1)
    # add the selected index to the attributes so it can be manipulated later
    attrs.selected_index = 1
    selected_index = attrs.selected_index

    # the fraction on the slider corresponding to the selected_index
    # this is only used after dragging
    sliderfraction = lift(selected_index, sliderrange) do i, r
        (i - 1) / (length(r) - 1)
    end

    dragging = Node(false)

    # what the slider actually displays currently (also during dragging when
    # the slider position is in an "invalid" position given the slider's range)
    displayed_sliderfraction = Node(0.0)

    on(sliderfraction) do frac
        # only update displayed fraction through sliderfraction if not dragging
        # dragging overrides the value so there is clear mouse interaction
        if !dragging[]
            displayed_sliderfraction[] = frac
        end
    end

    on(selected_index) do i
        value[] = sliderrange[][i]
    end

    # initialize slider value with closest from range
    selected_index[] = closest_index(sliderrange[], startvalue[])

    buttonpoint = lift(sliderbox, horizontal, displayed_sliderfraction, buttonradius,
            buttonstrokewidth) do bb, horizontal, sf, brad, bstw

        pad = brad + bstw

        if horizontal
            [Point2f0(left(bb) + pad + (width(bb) - 2pad) * sf, bottom(bb) + height(bb) / 2)]
        else
            [Point2f0(left(bb) + 0.5f0 * width(bb), bottom(bb) + pad + (height(bb) - 2pad) * sf)]
        end
    end

    linepoints = lift(endpoints, buttonpoint) do eps, bp
        [eps[1], bp[1], bp[1], eps[2]]
    end

    linecolors = lift(color_active_dimmed, color_inactive) do ca, ci
        [ca, ci]
    end

    linesegs = linesegments!(subscene, linepoints, color = linecolors, linewidth = linewidth, raw = true)[end]

    linestate = addmousestate!(subscene, linesegs)

    bsize = Node{Float32}(buttonradius[] * 2f0)

    bcolor = Node{Any}(buttoncolor_inactive[])

    button = scatter!(subscene, buttonpoint, markersize = bsize, color = bcolor,
        strokewidth = buttonstrokewidth, strokecolor = color_active_dimmed, raw = true)[end]

    buttonstate = addmousestate!(subscene, button)

    # on(buttonstate) do state
    #     typ = typeof(state.typ)
    #     if typ in (MouseDown, MouseDrag, MouseDragStart, MouseDragStop)
    #         bcolor[] = color_active[]
    #     end
    # end

    onmouseleftdown(buttonstate) do state
        bcolor[] = color_active[]
    end

    onmouseleftup(buttonstate) do state
        bcolor[] = buttoncolor_inactive[]
    end

    onmouseleftdrag(buttonstate) do state

        pad = buttonradius[] + buttonstrokewidth[]

        dragging[] = true
        dif = state.pos - state.prev
        fraction = if horizontal[]
            dif[1] / (width(sliderbox[]) - 2pad)
        else
            dif[2] / (height(sliderbox[]) - 2pad)
        end
        if fraction != 0.0f0
            @async begin
                newfraction = min(max(displayed_sliderfraction[] + fraction, 0f0), 1f0)
                displayed_sliderfraction[] = newfraction

                newindex = closest_fractionindex(sliderrange[], newfraction)
                if selected_index[] != newindex
                    selected_index[] = newindex
                end
            end
        end
    end

    onmouseleftdragstop(buttonstate) do state
        dragging[] = false
        # adjust slider to closest legal value
        sliderfraction[] = sliderfraction[]
    end

    scenestate = addmousestate!(subscene)

    onmouseleftclick(scenestate) do state

        pad = buttonradius[] + buttonstrokewidth[]

        pos = state.pos
        dim = horizontal[] ? 1 : 2
        frac = (pos[dim] - endpoints[][1][dim] - pad) / (endpoints[][2][dim] - endpoints[][1][dim] - 2pad)
        selected_index[] = closest_fractionindex(sliderrange[], frac)
    end

    onmouseleftdoubleclick(scenestate) do state
        selected_index[] = closest_index(sliderrange[], startvalue[])
    end

    onmouseenter(scenestate) do state
        # bcolor[] = color_active[]
        linecolors[] = [color_active[], color_inactive[]]
        button.strokecolor = color_active[]
    end

    onmouseout(scenestate) do state
        # bcolor[] = color_inactive[]
        linecolors[] = [color_active_dimmed[], color_inactive[]]
        button.strokecolor = color_active_dimmed[]
    end

    onany(buttonradius, horizontal) do br, horizontal
        protrusions[] = if horizontal
            GridLayoutBase.RectSides{Float32}(br, br, 0, 0)
        else
            GridLayoutBase.RectSides{Float32}(0, 0, br, br)
        end
    end

    # trigger protrusions using one observable
    buttonradius[] = buttonradius[]

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    LSlider(parent, layoutobservables, attrs, decorations)
end

function valueindex(sliderrange, value)
    for (i, val) in enumerate(sliderrange)
        if val == value
            return i
        end
    end
    nothing
end

function closest_fractionindex(sliderrange, fraction)
    n = length(sliderrange)
    onestepfrac = 1 / (n - 1)
    i = round(Int, fraction / onestepfrac) + 1
    min(max(i, 1), n)
end

function closest_index(sliderrange, value)
    for (i, val) in enumerate(sliderrange)
        if val == value
            return i
        end
    end
    # if the value wasn't found this way try inexact
    closest_index_inexact(sliderrange, value)
end

function closest_index_inexact(sliderrange, value)
    distance = Inf
    selected_i = 0
    for (i, val) in enumerate(sliderrange)
        newdist = abs(val - value)
        if newdist < distance
            distance = newdist
            selected_i = i
        end
    end
    selected_i
end

"""
Set the `slider` to the value in the slider's range that is closest to `value`.
"""
function set_close_to!(slider, value)
    closest = closest_index(slider.range[], value)
    slider.selected_index = closest
end


text:
using Documenter, SlimOptim

makedocs(sitename="Slim Optimization toolbox",
         doctest=false, clean=true,
         authors="Mathias Louboutin",
         pages = Any[
             "Home" => "index.md",
             "Tutorials" => Any[
                 "tutorials/01_denoising.md",
                 "tutorials/02_simple_constrained.md",
                 "tutorials/03_constr_fwi_judi.md"
             ]
         ])

deploydocs(repo="github.com/slimgroup/SlimOptim.jl")


text:
using Test
using Scruff
using Scruff.SFuncs
using Scruff.Operators

@testset "net" begin
    @testset "Build NetworkSFunc" begin
        x = NetworkInput{Int}()

        struct AddC <: SFunc{Tuple{Int},Int}
            c::Int
        end

        addc = AddC(2)

        struct MulC <: SFunc{Tuple{Int},Int}
            c::Int
        end

        struct Add <: SFunc{Tuple{Int,Int},Int}
        end

        mulc = MulC(-1)

        add = Add()

        sfuncs = (addc, mulc, add)
        parents = Dict{SFunc,Any}(addc=>[x], mulc=>[x], add=>[addc, mulc])
        outputs = (add,)

        net = NetworkSFunc((x,),
                           sfuncs,
                           parents,
                           outputs)

        @test 1 == 1
    end

    @testset "Netwalk Ops" begin
        x = NetworkInput{Int}()

        x_y = DiscreteCPT([1,2], Dict((1,) => [0.9, 0.1], (2,) => [0.1, 0.9]))
        y_z = DiscreteCPT([1,2], Dict((1,) => [0.1, 0.9], (2,) => [0.9, 0.1]))

        sfuncs = (x_y, y_z)
        parents = Dict{SFunc,Any}(x_y=>[x],y_z=>[x_y])
        outputs = (y_z,)

        net = NetworkSFunc((x,),
                           sfuncs,
                           parents,
                           outputs)

        y1 = sample(net, (1,))
        y2 = sample(net, (1,))
        cpdf2 = logcpdf(net, (1,), y2)
        # cpdf2, y2 = sample_logcpdf(net, 1)

        N = 100000
        correct = [0.18 0.82
                   0.82 0.18]
        for x in 1:2
            for z in 1:2
                cpdf3 = 0.0
                for i in 1:N
                    cpdf3 += exp(logcpdf(net, (x,), z))
                end
                cpdf3 = cpdf3/N
                # println("P(z=$z|x=$x) = $cpdf3")
                @test isapprox(correct[x,z], cpdf3, rtol=5e-2)
            end
        end
    end
end


text:
# Putting the override in a expr because it needs to be evaluated in this process and the
# distributed process.
const show_richest_override = :(
    # Override for `show_richest`.
    # Override the full method because allmimes was replaced by the compiler.
    function PlutoRunner.show_richest(io::IO, @nospecialize(x))::Tuple{<:Any,MIME}
        nonplutomimes = filter(m -> !occursin("pluto.tree", string(m)), PlutoRunner.allmimes)
        # ugly code to fix an ugly performance problem
        local mime = nothing
        for m in nonplutomimes
            if PlutoRunner.pluto_showable(m, x)
                mime = m
                break
            end
        end

        if mime in PlutoRunner.imagemimes
            show(io, mime, x)
            nothing, mime
        elseif mime isa MIME"application/vnd.pluto.table+object"
            PlutoRunner.table_data(x, IOContext(io, :compact => true)), mime
        elseif mime isa MIME"text/latex"
            # Some reprs include $ at the start and end.
            # We strip those, since Markdown.LaTeX should contain the math content.
            # (It will be rendered by MathJax, which is math-first, not text-first.)
            texed = repr(mime, x)
            Markdown.html(io, Markdown.LaTeX(strip(texed, ('$', '\n', ' '))))
            nothing, MIME"text/html"()
        else
            # the classic:
            show(io, mime, x)
            nothing, mime
        end
    end
)

# Giving `val` back so that the stacktrace can be printed easily.
const format_output_exception_override = :(
    function PlutoRunner.format_output(val::CapturedException; context=default_iocontext)
        msg = sprint(PlutoRunner.try_showerror, val.ex)
        mime = MIME"application/vnd.pluto.stacktrace+object"()
        return (Dict{Symbol,Any}(:msg => msg, :stacktrace => val), mime)
    end
)

const OLD_PROCESS_PREAMBLE = WorkspaceManager.process_preamble()

# Override the preamble to disable Pluto's pretty printing.
WorkspaceManager.process_preamble() = quote
    $OLD_PROCESS_PREAMBLE
    $(show_richest_override)
    $(format_output_exception_override)
end

# Yes. These overrides are used when `use_distributed=false`.
eval(show_richest_override)
eval(format_output_exception_override)