text:
####
#### Lazy wrappers for LHS matrix and RHS vector for augmented Tikhonov-regularized NNLS problems
####

struct PaddedVector{T, Vd <: AbstractVector{T}} <: AbstractVector{T}
    d::Vd
    pad::Int
end
Base.size(x::PaddedVector) = (length(x.d) + x.pad,)

function Base.copyto!(y::AbstractVector{T}, x::PaddedVector{T}) where {T}
    @assert size(x) == size(y)
    @unpack d, pad = x
    m = length(x.d)
    @inbounds @simd ivdep for i in 1:m
        y[i] = d[i]
    end
    @inbounds @simd ivdep for i in m+1:m+pad
        y[i] = zero(T)
    end
    return y
end

struct TikhonovPaddedMatrix{T, MC <: AbstractMatrix{T}} <: AbstractMatrix{T}
    C::MC
    μ::Base.RefValue{T}
end
TikhonovPaddedMatrix(C::AbstractMatrix, μ) = TikhonovPaddedMatrix(C, Ref(μ))
Base.size(A::TikhonovPaddedMatrix) = ((m,n) = size(A.C); return (m+n,n))

function Base.copyto!(B::AbstractMatrix{T}, A::TikhonovPaddedMatrix{T}) where {T}
    @assert size(A) == size(B)
    @unpack C, μ = A
    m, n = size(C)
    @inbounds @simd ivdep for j in 1:n
        @simd ivdep for i in 1:m
            B[i,j] = C[i,j]
        end
        @simd ivdep for i in m+1:m+n
            B[i,j] = ifelse(i == m+j, μ[], zero(T))
        end
    end
    return B
end

cache!(cache::AbstractArray, x, f) = cache !== nothing && push!(cache, (; x, f))
cache!(cache::Dict, x, f) = cache !== nothing && (cache[x] = f)

lin_interp(x, x₁, x₂, y₁, y₂) = y₁ + (y₂ - y₁) * (x - x₁)  / (x₂ - x₁)
exp_interp(x, x₁, x₂, y₁, y₂) = y₁ + log1p(expm1(y₂ - y₁) * (x - x₁)  / (x₂ - x₁))

####
#### Tikhonov regularized NNLS problem
####

struct NNLSTikhonovRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W
end
function NNLSTikhonovRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(TikhonovPaddedMatrix(C, T(NaN)), PaddedVector(d, n))
    NNLSTikhonovRegProblem(C, d, m, n, nnls_work)
end

mu(work::NNLSTikhonovRegProblem) = work.nnls_work.C.μ[]
mu!(work::NNLSTikhonovRegProblem, μ) = work.nnls_work.C.μ[] = μ

function solve!(work::NNLSTikhonovRegProblem, μ)
    # Set regularization parameter and solve NNLS problem
    mu!(work, μ)
    solve!(work.nnls_work)
    return solution(work)
end

solution(work::NNLSTikhonovRegProblem) = solution(work.nnls_work)

resnorm_sq(work::NNLSTikhonovRegProblem) = chi2(work)
resnorm(work::NNLSTikhonovRegProblem) = sqrt(resnorm_sq(work))

seminorm_sq(work::NNLSTikhonovRegProblem) = solution(work) ⋅ solution(work)
seminorm(work::NNLSTikhonovRegProblem) = sqrt(seminorm_sq(work))

loss(work::NNLSTikhonovRegProblem) = chi2(work.nnls_work)

reg(work::NNLSTikhonovRegProblem) = mu(work)^2 * seminorm_sq(work)

chi2(work::NNLSTikhonovRegProblem) = max(loss(work) - reg(work), 0)

# Extract columns spanning solution space of regularized NNLS problem
function subproblem(work::NNLSTikhonovRegProblem{T}) where {T}
    @unpack C, d, m, n, nnls_work = work
    μ = mu(work)
    x = solution(work)
    Jnz = findall(!≈(0), x)
    nnz = length(Jnz)
    return TikhonovPaddedMatrix(C[:,Jnz], μ), PaddedVector(d, nnz), x[Jnz]
end

function ∇reg(work::NNLSTikhonovRegProblem)
    @unpack C, d = work
    μ = mu(work)
    x = solution(work)
    C′, d′, x′ = subproblem(work)
    b′ = PaddedVector(C′.C' \ x′, d′.pad)
    ∇nnls_work = NNLSProblem(C′, b′)
    solve!(∇nnls_work)
    ∇x′ = -solution(∇nnls_work)
    ∇μ² = 2 * ((C′.C * x′ - d′.d) ⋅ (C′.C * ∇x′)) #+ (x′ ⋅ x′) + 2 * μ^2 * (x′ ⋅ ∇x′)
    ∇μ = 2 * μ * ∇μ²
    return ∇μ
end

function chi2factor_relerr!(work::NNLSTikhonovRegProblem, logμ, ∇logμ = nothing; χ²target)
    μ = exp(logμ)
    solve!(work, μ)
    χ² = chi2(work)
    relerr = log(χ² / χ²target) # better behaved than χ² / χ²target - 1 for large χ²
    if ∇logμ !== nothing && length(∇logμ) > 0
        ∂χ²_∂logμ = ∇reg(work)
        ∂relerr_∂logμ = μ * ∂χ²_∂logμ / χ²
        @inbounds ∇logμ[1] = ∂relerr_∂logμ
    end
    return relerr
end
chi2factor_relerr⁻¹(relerr; χ²target) = χ²target * exp(relerr)

function chi2factor_loss!(work::NNLSTikhonovRegProblem, logμ, ∇logμ = nothing; χ²target)
    relerr = chi2factor_relerr!(work, logμ, ∇logμ; χ²target)
    loss = abs(relerr)
    if ∇logμ !== nothing && length(∇logμ) > 0
        @inbounds ∂relerr_∂logμ = ∇logμ[1]
        ∂loss_∂logμ = sign(relerr) * ∂relerr_∂logμ
        @inbounds ∇logμ[1] = ∂loss_∂logμ
    end
    return loss
end

"""
Helper struct which wraps `N` caches of `NNLSTikhonovRegProblem` workspaces.
Useful for optimization problems where the last function call may not be
the optimium, but perhaps it was one or two calls previous and is still in the
`NNLSTikhonovRegProblemCache` and a recomputation can be avoided.
"""
struct NNLSTikhonovRegProblemCache{N,W}
    cache::NTuple{N,W}
    idx::Base.RefValue{Int}
end
function NNLSTikhonovRegProblemCache(C::AbstractMatrix{T}, d::AbstractVector{T}, ::Val{N} = Val(5)) where {T,N}
    cache = ntuple(_ -> NNLSTikhonovRegProblem(C, d), N)
    idx = Ref(1)
    return NNLSTikhonovRegProblemCache(cache, idx)
end
increment_cache_index!(work::NNLSTikhonovRegProblemCache{N}) where {N} = (work.idx[] = mod1(work.idx[] + 1, N))
set_cache_index!(work::NNLSTikhonovRegProblemCache{N}, i) where {N} = (work.idx[] = mod1(i, N))
reset_cache!(work::NNLSTikhonovRegProblemCache{N}) where {N} = foreach(w -> mu!(w, NaN), work.cache)
get_cache(work::NNLSTikhonovRegProblemCache) = work.cache[work.idx[]]

function solve!(work::NNLSTikhonovRegProblemCache, μ)
    i = findfirst(w -> μ == mu(w), work.cache)
    if i === nothing
        increment_cache_index!(work)
        solve!(get_cache(work), μ)
    else
        set_cache_index!(work, i)
    end
    return solution(get_cache(work))
end

####
#### Chi2 method for choosing Tikhonov regularization parameter
####

struct NNLSChi2RegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSChi2RegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSChi2RegProblem(C, d, m, n, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSChi2RegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_chi2(C::AbstractMatrix, d::AbstractVector, Chi2Factor::Real)

Returns the regularized NNLS solution, X, that incurrs an increase in ``\\chi^2`` approximately by a factor of `Chi2Factor`.
The regularized NNLS problem solved internally is:

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||x||_2^2
```

where ``\\mu`` is determined by approximating a solution to the nonlinear equation

```math
\\frac{\\chi^2(\\mu)}{\\chi^2_{min}} = \\mathrm{Chi2Factor}
\\quad
\\text{where}
\\quad
\\chi^2_{min} = \\chi^2(\\mu = 0)
```

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data
- `Chi2Factor::Real`: Desired ``\\chi^2`` increase due to regularization

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Actual increase ``\\chi^2(\\mu)/\\chi^2_{min}``, which will be approximately equal to the input `Chi2Factor`
"""
function lsqnonneg_chi2(C, d, Chi2Factor; kwargs...)
    work = lsqnonneg_chi2_work(C, d)
    lsqnonneg_chi2!(work, Chi2Factor; kwargs...)
end
lsqnonneg_chi2_work(C, d) = NNLSChi2RegProblem(C, d)

function lsqnonneg_chi2!(work::NNLSChi2RegProblem{T}, Chi2Factor::T; bisection = true, legacy = false) where {T}
    # Non-regularized solution
    solve!(work.nnls_work)
    chi2_min = chi2(work.nnls_work)

    # Prepare to solve
    χ²target = Chi2Factor * chi2_min
    reset_cache!(work.nnls_work_smooth_cache)

    if legacy
        # Use the legacy algorithm: double μ starting from an initial guess, then interpolate the root using a cubic spline fit
        mu_final, chi2_final = chi2factor_search_from_minimum(chi2_min, Chi2Factor; legacy = legacy) do μ
            μ == 0 && return chi2_min
            solve!(work.nnls_work_smooth_cache, μ)
            return chi2(get_cache(work.nnls_work_smooth_cache))
        end
        if mu_final == 0
            x_final = solution(work.nnls_work)
        else
            x_final = solve!(work.nnls_work_smooth_cache, mu_final)
        end

    elseif bisection
        # Find bracketing interval containing root, then perform bisection search
        f = function (logμ)
            increment_cache_index!(work.nnls_work_smooth_cache)
            return chi2factor_relerr!(get_cache(work.nnls_work_smooth_cache), logμ; χ²target)
        end
        cache = NamedTuple{(:x, :f), NTuple{2,T}}[]
        a, b, fa, fb = bracketing_interval(f, T(-4.0), T(1.0), T(1.5); maxiters = 6, cache)
        bisect(f, a, b, fa, fb; xtol = T(0.05), ftol = (Chi2Factor-1)/100, cache)

        # Spline rootfinding on evaluated points to improve accuracy
        sort!(cache; by = d -> d.x)
        logmu_root = spline_root([d.x for d in cache], [d.f for d in cache]; deg_spline = 1)
        logmu_root !== nothing && !any(d -> d.x ≈ logmu_root, cache) && cache!(cache, logmu_root, f(logmu_root))

        # Return regularization which minimizes relerr
        _, i = findmin([abs(d.f) for d in cache])
        logmu_final, relerr_final = cache[i]
        mu_final, chi2_final = exp(logmu_final), chi2factor_relerr⁻¹(relerr_final; χ²target)
        x_final = solve!(work.nnls_work_smooth_cache, mu_final)

    else
        # Instead of rootfinding, reformulate as a minimization problem. Solve using NLopt.
        alg = :LN_COBYLA # local, gradient-free, linear approximation of objective
        # alg = :LN_BOBYQA # local, gradient-free, quadratic approximation of objective
        # alg = :GN_AGS # global, gradient-free, hilbert curve based dimension reduction
        # alg = :LN_NELDERMEAD # local, gradient-free, simplex method
        # alg = :LN_SBPLX # local, gradient-free, subspace searching simplex method
        # alg = :LD_CCSAQ # local, first-order (rough ranking: [:LD_MMA, :LD_SLSQP, :LD_LBFGS, :LD_CCSAQ, :LD_AUGLAG])

        opt = NLopt.Opt(alg, 1)
        opt.lower_bounds  = -8.0
        opt.upper_bounds  = 2.0
        opt.xtol_rel      = 0.05
        opt.min_objective = function (logμ, ∇logμ)
            @inbounds _logμ = logμ[1]
            increment_cache_index!(work.nnls_work_smooth_cache)
            loss = chi2factor_loss!(get_cache(work.nnls_work_smooth_cache), _logμ, ∇logμ; χ²target)
            return Float64(loss)
        end
        minf, minx, ret   = NLopt.optimize(opt, [-4.0])

        mu_final = exp(T(minx[1]))
        x_final = solve!(work.nnls_work_smooth_cache, mu_final)
        chi2_final = chi2(get_cache(work.nnls_work_smooth_cache))
    end

    return (x = x_final, mu = mu_final, chi2factor = chi2_final/chi2_min)
end

function chi2factor_search_from_minimum(f, χ²min::T, χ²fact::T, μmin::T = T(1e-3), μfact = T(2.0); legacy = false) where {T}
    # Minimize energy of spectrum; loop to find largest μ that keeps chi-squared in desired range
    μ_cache = T[zero(T)]
    χ²_cache = T[χ²min]
    μnew = μmin
    while true
        # Cache function value at μ = μnew
        χ²new = f(μnew)
        push!(μ_cache, μnew)
        push!(χ²_cache, χ²new)

        # Break when χ²fact reached, else increase regularization
        (χ²new >= χ²fact * χ²min) && break
        μnew *= μfact
    end

    # Solve χ²(μ) = χ²fact * χ²min using a spline fitting root finding method
    if legacy
        # Legacy algorithm fits spline to all (μ, χ²) values observed, including for μ=0.
        # This poses several problems:
        #   1) while unlikely, it is possible for the spline to return a negative regularization parameter
        #   2) the μ values are exponentially spaced, leading to poorly conditioned splines
        μ = spline_root_legacy_slow(μ_cache, χ²_cache, χ²fact * χ²min)
        μ = μ === nothing ? μmin : μ
    else
        if length(μ_cache) == 2
            # Solution is contained in [0,μmin]; `spline_root` with two points performs root finding via simple linear interpolation
            μ = spline_root(μ_cache, χ²_cache, χ²fact * χ²min; deg_spline = 1)
            μ = μ === nothing ? μmin : μ
        else
            # Perform spline fit on log-log scale on data with μ > 0. This solves the above problems with the legacy algorithm:
            #   1) Root is found in terms of logμ, guaranteeing μ > 0
            #   2) logμ is linearly spaced, leading to well-conditioned splines
            logμ = spline_root(log.(μ_cache[2:end]), log.(χ²_cache[2:end]), log(χ²fact * χ²min); deg_spline = 1)
            μ = logμ === nothing ? μmin : exp(logμ)
        end
    end

    # Compute the final regularized solution
    χ² = f(μ)

    return μ, χ²
end

function chi2factor_search_from_guess(f, χ²min::T, χ²fact::T, μ₀::T = T(1e-2), μfact = T(1.5), μmin::T = T(1e-4)) where {T}
    # Find interval containing χ²target = χ²fact * χ²min
    χ²target = χ²fact * χ²min
    μnew = μ₀
    χ²new = f(μnew)
    logμ_cache = T[log(μnew)]
    logχ²_cache = T[log(χ²new)]

    if χ²new > χ²target
        while μnew > μmin && !(μnew ≈ μmin)
            μnew = max(μnew/μfact, μmin)
            χ²new = f(μnew)
            pushfirst!(logμ_cache, log(μnew))
            pushfirst!(logχ²_cache, log(χ²new))
            (χ²new < χ²target && length(logμ_cache) >= 3) && break
        end
    else
        while true
            μnew *= μfact
            χ²new = f(μnew)
            push!(logμ_cache, log(μnew))
            push!(logχ²_cache, log(χ²new))
            (χ²new > χ²target && length(logμ_cache) >= 3) && break
        end
    end

    if logμ_cache[1] ≈ log(μmin) && logχ²_cache[1] > log(χ²target)
        # μ decreased to μmin but χ²target was not reached; linearly interpolate between μ=0 and μ=μmin points
        μ = spline_root([zero(T), μmin], [log(χ²min), logχ²_cache[1]], log(χ²target))
        μ = μ === nothing ? μmin : μ
        χ² = f(μ)
    else
        # Find optimal μ and evaluate at the interpolated solution
        logμ = spline_root(logμ_cache, logχ²_cache, log(χ²target))
        μ = logμ === nothing ? μmin : exp(logμ)
        χ² = f(μ)
    end

    return μ, χ²
end

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

function secant_method(f, xs; xtol = zero(float(real(first(xs)))), xrtol = 8*eps(one(float(real(first(xs))))), maxiters = 1000, cache = nothing)
    if length(xs) == 1 # secant needs a, b; only a given
        a  = float(xs[1])
        h  = eps(one(real(a)))^(1/3)
        da = h * oneunit(a) + abs(a) * h^2 # adjust for if eps(a) > h
        b  = a + da
    else
        a, b = promote(float(xs[1]), float(xs[2]))
    end
    fa, fb = f(a), f(b)
    cache!(cache, a, fa)
    cache!(cache, b, fb)
    secant(f, a, b, fa, fb; xtol, xrtol, maxiters, cache)
end

function secant(f, a::T, b::T, fa::T, fb::T; xtol = zero(T), xrtol = 8eps(T), maxiters = 1000, cache = nothing) where {T}
    # No function change; return arbitrary endpoint
    if fb == fa
        return (a, fa)
    end

    cnt = 0
    mbest = abs(fa) < abs(fb) ? a : b
    fbest = min(abs(fa), abs(fb))
    uatol = xtol / oneunit(xtol) * oneunit(real(a))
    adjustunit = oneunit(real(fb)) / oneunit(real(b))

    while cnt < maxiters
        m = b - (b - a) * fb / (fb - fa)
        fm = f(m)
        cache!(cache, m, fm)

        abs(fm) < abs(fbest) && ((mbest, fbest) = (m, fm))
        iszero(fm) && return (m, fm)
        isnan(fm) || isinf(fm) && return (mbest, fbest) # function failed; bail out
        abs(fm) <= adjustunit * max(uatol, abs(m) * xrtol) && return (m, fm)
        fm == fb && return (m, fm)

        a, b, fa, fb = b, m, fb, fm
        cnt += 1
    end

    return (mbest, fbest) # maxiters reached
end

function bisection_method(f, a::Number, b::Number; xtol = nothing, xrtol = nothing, ftol = nothing, cache = nothing, maxiters = 1000)
    x₁, x₂ = float.((a,b))
    y₁, y₂ = f(x₁), f(x₂)
    cache!(cache, x₁, y₁)
    cache!(cache, x₂, y₂)

    T = eltype(x₁)
    xtol = xtol === nothing ? zero(T) : abs(xtol)
    xrtol = xrtol === nothing ? zero(one(T)) : abs(xrtol)
    ftol = ftol === nothing ? zero(T) : abs(ftol)

    bisect(f, x₁, x₂, y₁, y₂; xtol, xrtol, ftol, cache, maxiters)
end

function bisect(f, x₁::T, x₂::T, y₁::T, y₂::T; xtol = zero(T), xrtol = zero(one(T)), ftol = zero(T), cache = nothing, maxiters = 1000) where {T}
    # No sign change; return arbitrary endpoint
    if y₁ * y₂ >= 0
        return (x₁, y₁)
    end

    if y₂ < 0
        x₁, x₂, y₁, y₂ = x₂, x₁, y₂, y₁
    end

    xₘ = (x₁ + x₂)/2
    yₘ = f(xₘ)
    cache!(cache, xₘ, yₘ)

    cnt = 1
    while cnt < maxiters
        if iszero(yₘ) || isnan(yₘ) || abs(x₁ - x₂) <= xtol + max(abs(x₁), abs(x₂)) * xrtol || abs(yₘ) <= ftol
            return (xₘ, yₘ)
        end

        if yₘ < 0
            x₁, y₁ = xₘ, yₘ
        else
            x₂, y₂ = xₘ, yₘ
        end

        xₘ = (x₁ + x₂)/2
        yₘ = f(xₘ)
        cache!(cache, xₘ, yₘ)

        cnt += 1
    end

    return (xₘ, yₘ)
end

function bracketing_interval(f, a, δ, dilate = 1; maxiters = 1000, cache = nothing)
    # Initialize cache
    fa = f(a)
    cache!(cache, a, fa)
    cnt = 0
    if fa > 0
        b = a - δ
        fb = f(b)
        cache!(cache, b, fb)
        δ *= dilate
        while fb > 0 && cnt < maxiters
            a, fa = b, fb
            b = a - δ
            fb = f(b)
            cache!(cache, b, fb)
            δ *= dilate
            cnt += 1
        end
    else
        b = a + δ
        fb = f(b)
        cache!(cache, b, fb)
        δ *= dilate
        while fb < 0 && cnt < maxiters
            a, fa = b, fb
            b = a + δ
            fb = f(b)
            cache!(cache, b, fb)
            δ *= dilate
            cnt += 1
        end
    end
    return a, b, fa, fb
end

####
#### L-curve method for choosing Tikhonov regularization parameter
####

struct NNLSLCurveRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSLCurveRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSLCurveRegProblem(C, d, m, n, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSLCurveRegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_lcurve(C::AbstractMatrix, d::AbstractVector)

Returns the regularized NNLS solution, X, of the equation

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||L x||_2^2
```

where ``L`` is the identity matrix and ``\\mu`` is chosen by locating the corner of the "L-curve".

Details of L-curve theory and the Generalized Cross-Validation (GCV) method can be found in:
[Hansen, P.C., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580](https://doi.org/10.1137/1034115)

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Resulting increase in ``\\chi^2`` relative to unregularized (``\\mu = 0``) solution
"""
function lsqnonneg_lcurve(C, d)
    work = lsqnonneg_lcurve_work(C, d)
    lsqnonneg_lcurve!(work)
end
lsqnonneg_lcurve_work(C, d) = NNLSLCurveRegProblem(C, d)

function lsqnonneg_lcurve!(work::NNLSLCurveRegProblem{T,N}) where {T,N}
    # Compute the regularization using the L-curve method
    reset_cache!(work.nnls_work_smooth_cache)
    logmu_bounds = (T(-8), T(2))
    logmu_final = lcurve_corner(logmu_bounds...) do μ
        solve!(work.nnls_work_smooth_cache, μ)
        ξ = log(resnorm(get_cache(work.nnls_work_smooth_cache)))
        η = log(seminorm(get_cache(work.nnls_work_smooth_cache)))
        return SA{T}[ξ, η]
    end

    # Return the final regularized solution
    mu_final = exp(logmu_final)
    x_final = solve!(work.nnls_work_smooth_cache, mu_final)
    x_unreg = solve!(work.nnls_work)
    chi2factor_final = chi2(get_cache(work.nnls_work_smooth_cache)) / chi2(work.nnls_work)

    return (x = x_final, mu = mu_final, chi2factor = chi2factor_final)
end

"""
    lcurve_corner(f, xlow, xhigh)

Find the corner of the L-curve via curvature maximization using Algorithm 1 from:

A. Cultrera and L. Callegaro, “A simple algorithm to find the L-curve corner in the regularization of ill-posed inverse problems”.
IOPSciNotes, vol. 1, no. 2, p. 025004, Aug. 2020, doi: 10.1088/2633-1357/abad0d
"""
function lcurve_corner(f, xlow::T = -8.0, xhigh::T = 2.0; xtol::T = 0.05, Ptol::T = 0.05, Ctol::T = 0.01, cache = nothing, refine = false, backtracking = true, verbose = false, kwargs...) where {T}
    # Initialize state
    msg(s, state) = verbose && (@info "$s: [x⃗, P⃗, C⃗] = "; display(hcat(state.x⃗, state.P⃗, [cache[x].C for x in state.x⃗])))
    cache === nothing && (cache = Dict{T, NamedTuple{(:P, :C), Tuple{SVector{2,T}, T}}}())
    state = LCurveCornerState(f, T(xlow), T(xhigh), cache; kwargs...)
    state_cache = backtracking ? [state] : nothing

    # Tolerances are relative to initial curve size
    Ptopleft, Pbottomright = state.P⃗[1], state.P⃗[4]
    Ptol = norm(Ptopleft - Pbottomright) * Ptol
    Ctol = norm(Ptopleft - Pbottomright) * Ctol

    # For very small regularization, points on the L-curve may be extremely close,
    # and curvature calculations will be highly ill-posed. Similarly, while it should
    # not typically occur, additionally check for points on the L-curve becoming
    # extremely close for large regularization. Points which don't satisfy `Pfilter`
    # are assigned curvature -Inf.
    Pfilter = P -> min(norm(P - Ptopleft), norm(P - Pbottomright)) > Ctol
    update_curvature!(state, cache; Pfilter)
    msg("Starting", state)

    while true
        if backtracking
            # Find state with minimum diameter which contains the current best estimate maximum curvature point
            (x, (P, C)), _, _ = mapfindmax(((x, (P, C)),) -> C, collect(cache))
            for s in state_cache
                if (s.x⃗[2] == x || s.x⃗[3] == x) && abs(s.x⃗[4] - s.x⃗[1]) <= abs(state.x⃗[4] - state.x⃗[1])
                    state = s
                end
            end
        end

        # Move state toward region of lower curvature
        if cache[state.x⃗[2]].C > cache[state.x⃗[3]].C
            state = move_left(f, state, cache)
            update_curvature!(state, cache; Pfilter)
            msg("C₂ > C₃; moved left", state)
        else
            state = move_right(f, state, cache)
            update_curvature!(state, cache; Pfilter)
            msg("C₃ ≥ C₂; moved right", state)
        end
        backtracking && push!(state_cache, state)
        is_converged(state; xtol, Ptol) && break
    end

    if refine
        x = refine!(f, state, cache; Pfilter)
        msg("Converged; refined solution", state)
    else
        x = cache[state.x⃗[2]].C > cache[state.x⃗[3]].C ? state.x⃗[2] : state.x⃗[3]
        msg("Converged", state)
    end

    return x
end

struct LCurveCornerState{T, V <: SVector{2,T}}
    x⃗::SVector{4,T} # grid of regularization parameters
    P⃗::SVector{4,V} # points (residual norm, solution seminorm) evaluated at x⃗
end
function LCurveCornerState(f, x₁::T, x₄::T, cache = nothing) where {T}
    x₂   = (T(φ) * x₁ + x₄) / (T(φ) + 1)
    x₃   = x₁ + (x₄ - x₂)
    x⃗    = SA{T}[x₁, x₂, x₃, x₄]
    P⃗    = f.(exp.(x⃗))
    foreach(i -> cache!(cache, x⃗[i], (P = P⃗[i], C = T(-Inf))), 1:4)
    return LCurveCornerState(x⃗, P⃗)
end

is_converged(state::LCurveCornerState; xtol, Ptol) = abs(state.x⃗[4] - state.x⃗[1]) < xtol || norm(state.P⃗[1] - state.P⃗[4]) < Ptol

function maybecall!(f, x::T, state::LCurveCornerState{T}, cache) where {T}
    for (_x, (_P, _C)) in cache
        x ≈ _x && return _P
    end
    P = f(exp(x)) # only recalculate if not cached
    cache!(cache, x, (P = P, C = T(-Inf)))
    return P
end

function move_left(f, state::LCurveCornerState{T}, cache) where {T}
    @unpack x⃗, P⃗ = state
    x⃗ = SA[x⃗[1], (T(φ) * x⃗[1] + x⃗[3]) / (T(φ) + 1), x⃗[2], x⃗[3]]
    P⃗ = SA[P⃗[1], maybecall!(f, x⃗[2], state, cache), P⃗[2], P⃗[3]] # only P⃗[2] is recalculated
    return LCurveCornerState(x⃗, P⃗)
end

function move_right(f, state::LCurveCornerState{T}, cache) where {T}
    @unpack x⃗, P⃗ = state
    x⃗ = SA[x⃗[2], x⃗[3], x⃗[2] + (x⃗[4] - x⃗[3]), x⃗[4]]
    P⃗ = SA[P⃗[2], P⃗[3], maybecall!(f, x⃗[3], state, cache), P⃗[4]] # only P⃗[3] is recalculated
    return LCurveCornerState(x⃗, P⃗)
end

function update_curvature!(state::LCurveCornerState{T}, cache; Pfilter = nothing, menger_curvature = true, global_search = false) where {T}
    @unpack x⃗, P⃗ = state
    Cfun(P₋, P, P₊) = menger_curvature ?
        menger(P₋, P, P₊) :
        360 - rad2deg(kahan_angle(P₋, P, P₊))
    for i in 1:4
        x, P, C = x⃗[i], P⃗[i], T(-Inf)
        if Pfilter === nothing || Pfilter(P)
            if global_search
                # Search for maximum curvature over all neighbours
                for (x₋, (P₋, _)) in cache, (x₊, (P₊, _)) in cache
                    (x₋ < x⃗[i] < x₊) && (C = max(C, Cfun(P₋, P, P₊)))
                end
            else
                # Compute curvature from nearest neighbours
                x₋, x₊ = T(-Inf), T(+Inf)
                P₋, P₊ = P, P
                for (_x, (_P, _)) in cache
                    (x₋ < _x < x ) && ((x₋, P₋) = (_x, _P))
                    (x  < _x < x₊) && ((x₊, P₊) = (_x, _P))
                end
                C = Cfun(P₋, P, P₊)
            end
        end
        cache[x] = (; P, C)
    end
    return state
end

function refine!(f, state::LCurveCornerState{T}, cache; Pfilter = nothing, analytical = false) where {T}
    # Fit spline to (negative) curvature estimates on grid and minimize over the spline
    if analytical
        x_opt, _, _ = maximize_curvature(state)
    else
        C_spl = make_spline(state.x⃗, [-cache[x].C for x in state.x⃗])
        x_opt, _ = spline_opt(C_spl)
    end
    maybecall!(f, x_opt, state, cache)
    update_curvature!(state, cache; Pfilter)
    (x_opt, (_, _)), _, _ = mapfindmax(((x, (P, C)),) -> C, collect(cache))
    return x_opt
end

function menger(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}}
    Δⱼₖ, Δₖₗ, Δₗⱼ = Pⱼ - Pₖ, Pₖ - Pₗ, Pₗ - Pⱼ
    P̄ⱼP̄ₖ, P̄ₖP̄ₗ, P̄ₗP̄ⱼ = Δⱼₖ ⋅ Δⱼₖ, Δₖₗ ⋅ Δₖₗ, Δₗⱼ ⋅ Δₗⱼ
    Cₖ = 2 * (Δⱼₖ × Δₖₗ) / √(P̄ⱼP̄ₖ * P̄ₖP̄ₗ * P̄ₗP̄ⱼ)
end

function menger(xⱼ::T, xₖ::T, xₗ::T, Pⱼ::V, Pₖ::V, Pₗ::V; interp_uniform = true, linear_deriv = true) where {T, V <: SVector{2,T}}
    if interp_uniform
        h  = min(abs(xₖ - xⱼ), abs(xₗ - xₖ)) / T(φ)
        h₋ = h₊ = h
        x₋, x₀, x₊ = xₖ - h, xₖ, xₖ + h
        P₀ = Pₖ
        P₋ = exp_interp.(x₋, xⱼ, xₖ, Pⱼ, Pₖ)
        P₊ = exp_interp.(x₊, xₖ, xₗ, Pₖ, Pₗ)
    else
        P₋, P₀, P₊ = Pⱼ, Pₖ, Pₗ
        x₋, x₀, x₊ = xⱼ, xₖ, xₗ
        h₋, h₊ = x₀ - x₋, x₊ - x₀
    end
    ξ₋, ξ₀, ξ₊ = P₋[1], P₀[1], P₊[1]
    η₋, η₀, η₊ = P₋[2], P₀[2], P₊[2]

    if linear_deriv
        ξ′  = (ξ₊ - ξ₋) / (h₊ + h₋)
        η′  = (η₊ - η₋) / (h₊ + h₋)
    else
        ξ′ = (h₋^2 * ξ₊ + (h₊ + h₋) * (h₊ - h₋) * ξ₀ - h₊^2 * ξ₋) / (h₊ * h₋ * (h₊ + h₋))
        η′ = (h₋^2 * η₊ + (h₊ + h₋) * (h₊ - h₋) * η₀ - h₊^2 * η₋) / (h₊ * h₋ * (h₊ + h₋))
    end

    ξ′′ = 2 * (h₋ * ξ₊ - (h₊ + h₋) * ξ₀ + h₊ * ξ₋) / (h₊ * h₋ * (h₊ + h₋))
    η′′ = 2 * (h₋ * η₊ - (h₊ + h₋) * η₀ + h₊ * η₋) / (h₊ * h₋ * (h₊ + h₋))

    return (ξ′ * η′′ - η′ * ξ′′) / √((ξ′^2 + η′^2)^3)
end

function menger(f; h = 1e-3)
    function menger_curvature_inner(x)
        fⱼ, fₖ, fₗ = f(x-h), f(x), f(x+h)
        Pⱼ, Pₖ, Pₗ = SA[x-h,fⱼ], SA[x,fₖ], SA[x+h,fₗ]
        menger(Pⱼ, Pₖ, Pₗ)
    end
end

function menger(x, y; h = 1e-3)
    function menger_curvature_inner(t)
        x₋, x₀, x₊ = x(t-h), x(t), x(t+h)
        y₋, y₀, y₊ = y(t-h), y(t), y(t+h)
        x′, x′′ = (x₊ - x₋) / 2h, (x₊ - 2x₀ + x₋) / h^2
        y′, y′′ = (y₊ - y₋) / 2h, (y₊ - 2y₀ + y₋) / h^2
        return (x′ * y′′ - y′ * x′′) / √((x′^2 + y′^2)^3)
    end
end

function menger(x::Dierckx.Spline1D, y::Dierckx.Spline1D)
    function menger_curvature_inner(t)
        x′  = Dierckx.derivative(x, t; nu = 1)
        x′′ = Dierckx.derivative(x, t; nu = 2)
        y′  = Dierckx.derivative(y, t; nu = 1)
        y′′ = Dierckx.derivative(y, t; nu = 2)
        return (x′ * y′′ - y′ * x′′) / √((x′^2 + y′^2)^3)
    end
end

function menger(y::Dierckx.Spline1D)
    function menger_curvature_inner(t)
        y′  = Dierckx.derivative(y, t; nu = 1)
        y′′ = Dierckx.derivative(y, t; nu = 2)
        return y′′ / √((1 + y′^2)^3)
    end
end

function maximize_curvature(state::LCurveCornerState{T}) where {T}
    # Maximize curvature and transform back from t-space to x-space
    @unpack x⃗, P⃗ = state
    t₁, t₂, t₃, t₄ = T(0), 1/T(3), 2/T(3), T(1)
    t_opt, P_opt, C_opt = maximize_curvature(P⃗...)
    x_opt =
        t₁ <= t_opt < t₂ ? lin_interp(t_opt, t₁, t₂, x⃗[1], x⃗[2]) :
        t₂ <= t_opt < t₃ ? lin_interp(t_opt, t₂, t₃, x⃗[2], x⃗[3]) :
        lin_interp(t_opt, t₃, t₄, x⃗[3], x⃗[4])
    return x_opt
end

function maximize_curvature(P₁::V, P₂::V, P₃::V, P₄::V; bezier = true) where {T, V <: SVector{2,T}}
    # Analytically maximize curvature of parametric cubic spline fit to data.
    #   see: https://cs.stackexchange.com/a/131032
    a, e = P₁[1], P₁[2]
    b, f = P₂[1], P₂[2]
    c, g = P₃[1], P₃[2]
    d, h = P₄[1], P₄[2]

    if bezier
        ξ = t -> a*(1-t)^3+3*b*(1-t)^2*t+3*c*(1-t)*t^2+d*t^3
        η = t -> e*(1-t)^3+3*f*(1-t)^2*t+3*g*(1-t)*t^2+h*t^3
        P = t -> SA{T}[ξ(t), η(t)]

        # In order to keep the equations in a more compact form, introduce the following substitutions:
        m = d-3*c+3*b-a
        n = c-2*b+a
        o = b-a
        p = h-3*g+3*f-e
        q = g-2*f+e
        r = f-e
    else
        ξ = t -> a*(t-1/3)*(t-2/3)*(t-1/1)/T((0/1-1/3)*(0/1-2/3)*(0/1-1/1)) + # (-9  * a * t^3)/2  + (+18 * a * t^2)/2 + (-11 * a * t)/2 + a
                 b*(t-0/1)*(t-2/3)*(t-1/1)/T((1/3-0/1)*(1/3-2/3)*(1/3-1/1)) + # (+27 * b * t^3)/2  + (-45 * b * t^2)/2 + (+18 * b * t)/2 +
                 c*(t-0/1)*(t-1/3)*(t-1/1)/T((2/3-0/1)*(2/3-1/3)*(2/3-1/1)) + # (-27 * c * t^3)/2  + (+36 * c * t^2)/2 + ( -9 * c * t)/2 +
                 d*(t-0/1)*(t-1/3)*(t-2/3)/T((1/1-0/1)*(1/1-1/3)*(1/1-2/3))   # (+9  * d * t^3)/2  + ( -9 * d * t^2)/2 + ( +2 * d * t)/2
        η = t -> e*(t-1/3)*(t-2/3)*(t-1/1)/T((0/1-1/3)*(0/1-2/3)*(0/1-1/1)) + # (-9  * e * t^3)/2  + (+18 * e * t^2)/2 + (-11 * e * t)/2 + e
                 f*(t-0/1)*(t-2/3)*(t-1/1)/T((1/3-0/1)*(1/3-2/3)*(1/3-1/1)) + # (+27 * f * t^3)/2  + (-45 * f * t^2)/2 + (+18 * f * t)/2 +
                 g*(t-0/1)*(t-1/3)*(t-1/1)/T((2/3-0/1)*(2/3-1/3)*(2/3-1/1)) + # (-27 * g * t^3)/2  + (+36 * g * t^2)/2 + ( -9 * g * t)/2 +
                 h*(t-0/1)*(t-1/3)*(t-2/3)/T((1/1-0/1)*(1/1-1/3)*(1/1-2/3))   # (+9  * h * t^3)/2  + ( -9 * h * t^2)/2 + ( +2 * h * t)/2
        P = t -> SA{T}[ξ(t), η(t)]

        # In order to keep the equations in a more compact form, introduce the following substitutions:
        m = ( -9 * a + 27 * b - 27 * c + 9 * d)/2
        n = (+18 * a - 45 * b + 36 * c - 9 * d)/6
        o = (-11 * a + 18 * b -  9 * c + 2 * d)/6
        p = ( -9 * e + 27 * f - 27 * g + 9 * h)/2
        q = (+18 * e - 45 * f + 36 * g - 9 * h)/6
        r = (-11 * e + 18 * f -  9 * g + 2 * h)/6
    end

    # This leads to the following simplified derivatives:
    ξ′  = t -> 3*(m*t^2+2*n*t+o)
    ξ′′ = t -> 6*(m*t+n)
    η′  = t -> 3*(p*t^2+2*q*t+r)
    η′′ = t -> 6*(p*t+q)

    # Curvature and its derivative:
    k  = t -> ((3*(m*t^2+2*n*t+o)) * (6*(p*t+q)) - (3*(p*t^2+2*q*t+r)) * (6*(m*t+n))) / ((3*(m*t^2+2*n*t+o))^2 + (3*(p*t^2+2*q*t+r))^2)^(3/2)
    k′ = t -> (-18*m*(p*t^2+2*q*t+r)+18*p*(m*t^2+2*n*t+o)-18*(m*t+n)*(2*p*t+2*q)+18*(2*m*t+2*n)*(p*t+q))/(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(3/2)-(3*(18*(p*t+q)*(m*t^2+2*n*t+o)-18*(m*t+n)*(p*t^2+2*q*t+r))*(18*(2*p*t+2*q)*(p*t^2+2*q*t+r)+18*(2*m*t+2*n)*(m*t^2+2*n*t+o)))/(2*(9*(p*t^2+2*q*t+r)^2+9*(m*t^2+2*n*t+o)^2)^(5/2))

    # Solve analytically
    coeffs = MVector{6,Complex{T}}(
        (1296*m*p^2+1296*m^3)*q-1296*n*p^3-1296*m^2*n*p,
        (1620*m*p^2+1620*m^3)*r+3240*m*p*q^2+(3240*m^2*n-3240*n*p^2)*q-1620*o*p^3+((-1620*m^2*o)-3240*m*n^2)*p,
        (5184*m*p*q+1296*n*p^2+6480*m^2*n)*r+1296*m*q^3-1296*n*p*q^2+((-6480*o*p^2)-1296*m^2*o+1296*m*n^2)*q+((-5184*m*n*o)-1296*n^3)*p,
        1296*m*p*r^2+(1944*m*q^2+6480*n*p*q-1296*o*p^2+1296*m^2*o+8424*m*n^2)*r-8424*o*p*q^2-6480*m*n*o*q+((-1296*m*o^2)-1944*n^2*o)*p,
        2592*n*p*r^2+(3888*n*q^2-2592*o*p*q+2592*m*n*o+3888*n^3)*r-3888*o*q^3+((-2592*m*o^2)-3888*n^2*o)*q,
        -324*m*r^3+(1944*n*q+324*o*p)*r^2+((-1944*o*q^2)-324*m*o^2+1944*n^2*o)*r-1944*n*o^2*q+324*o^3*p,
    )
    roots = MVector{6,Complex{T}}(undef)
    PolynomialRoots.roots5!(roots, coeffs, eps(T), true)

    # Check roots and return maximum
    t₁, t₄ = zero(T), one(T)
    k₁, k₄ = k(t₁), k(t₄)
    tmax, Pmax, kmax = k₁ > k₄ ? (t₁, P₁, k₁) : (t₄, P₄, k₄)
    for rᵢ in roots
        _t, _s = real(rᵢ), imag(rᵢ)
        !(_s ≈ 0) && continue # real roots only
        !(t₁ <= _t <= t₄) && continue # filter roots within range
        ((_k = k(_t)) > kmax) && ((tmax, Pmax, kmax) = (_t, P(_t), _k))
    end

    return tmax, Pmax, kmax
end

function directed_angle(v₁::V, v₂::V) where {T, V <: SVector{2,T}}
    α = atan(v₁[2], v₁[1]) - atan(v₂[2], v₂[1])
    α < 0 ? 2*T(π) + α : α
end
directed_angle(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}} = directed_angle(Pⱼ - Pₖ, Pₗ - Pₖ)

function kahan_angle(v₁::V, v₂::V) where {T, V <: SVector{2,T}}
    # Kahan's method for computing the angle between v₁ and v₂.
    #   see: https://scicomp.stackexchange.com/a/27694
    a, b, c = norm(v₁), norm(v₂), norm(v₁ - v₂)
    a, b = max(a,b), min(a,b)
    μ = b ≥ c ? c - (a - b) : (b - (a - c))
    num = ((a - b) + c) * max(μ, zero(μ))
    den = (a + (b + c)) * ((a - c) + b)
    α = 2 * atan(√(num / den))
    v₁ × v₂ > 0 ? 2*T(π) - α : α
end
kahan_angle(Pⱼ::V, Pₖ::V, Pₗ::V) where {V <: SVector{2}} = kahan_angle(Pⱼ - Pₖ, Pₗ - Pₖ)

####
#### GCV method for choosing Tikhonov regularization parameter
####

struct NNLSGCVRegProblem{T, MC <: AbstractMatrix{T}, Vd <: AbstractVector{T}, W1, W2}
    C::MC
    d::Vd
    m::Int
    n::Int
    Aμ::Matrix{T}
    C_buf::Matrix{T}
    Ct_buf::Matrix{T}
    CtC_buf::Matrix{T}
    nnls_work::W1
    nnls_work_smooth_cache::W2
end
function NNLSGCVRegProblem(C::AbstractMatrix{T}, d::AbstractVector{T}) where {T}
    m, n = size(C)
    Aμ = zeros(T, m, m)
    C_buf = zeros(T, m, n)
    Ct_buf = zeros(T, n, m)
    CtC_buf = zeros(T, n, n)
    nnls_work = NNLSProblem(C, d)
    nnls_work_smooth_cache = NNLSTikhonovRegProblemCache(C, d)
    NNLSGCVRegProblem(C, d, m, n, Aμ, C_buf, Ct_buf, CtC_buf, nnls_work, nnls_work_smooth_cache)
end

solution(work::NNLSGCVRegProblem) = solution(get_cache(work.nnls_work_smooth_cache))

"""
    lsqnonneg_gcv(C::AbstractMatrix, d::AbstractVector)

Returns the regularized NNLS solution, X, of the equation

```math
X = \\mathrm{argmin}_{x \\ge 0} ||Cx - d||_2^2 + \\mu^2 ||L x||_2^2
```

where ``L`` is the identity matrix and ``\\mu`` is chosen by the Generalized Cross-Validation (GCV) method.

Details of the GCV method and L-curve theory can be can be found in:
[Hansen, P.C., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580](https://doi.org/10.1137/1034115)

# Arguments
- `C::AbstractMatrix`: Decay basis matrix
- `d::AbstractVector`: Decay curve data

# Outputs
- `X::AbstractVector`: Regularized NNLS solution
- `mu::Real`: Resulting regularization parameter ``\\mu``
- `Chi2Factor::Real`: Resulting increase in ``\\chi^2`` relative to unregularized (``\\mu = 0``) solution
"""
function lsqnonneg_gcv(C, d)
    work = lsqnonneg_gcv_work(C, d)
    lsqnonneg_gcv!(work)
end
lsqnonneg_gcv_work(C, d) = NNLSGCVRegProblem(C, d)

function lsqnonneg_gcv!(work::NNLSGCVRegProblem{T,N}) where {T,N}
    # Find μ by minimizing the function G(μ) (GCV method)
    reset_cache!(work.nnls_work_smooth_cache)

    # opt = NLopt.Opt(:LN_COBYLA, 1) # local, gradient-free, linear approximation of objective
    opt = NLopt.Opt(:LN_BOBYQA, 1) # local, gradient-free, quadratic approximation of objective
    opt.lower_bounds  = -8.0
    opt.upper_bounds  = 2.0
    opt.xtol_rel      = 0.05
    opt.min_objective = (logμ, ∇logμ) -> Float64(gcv!(work, logμ[1]))
    minf, minx, ret   = NLopt.optimize(opt, [-4.0])

    # Return the final regularized solution
    mu_final = exp(T(minx[1]))
    x_final  = solve!(work.nnls_work_smooth_cache, mu_final)
    x_unreg  = solve!(work.nnls_work)
    chi2factor_final = chi2(get_cache(work.nnls_work_smooth_cache)) / chi2(work.nnls_work)

    return (x = x_final, mu = mu_final, chi2factor = chi2factor_final)
end

# Implements equation (32) from:
#
#   Analysis of Discrete Ill-Posed Problems by Means of the L-Curve
#   Hansen et al. 1992 (https://epubs.siam.org/doi/10.1137/1034115)
#
# where here A = C, b = d, λ = μ, and L = identity.
function gcv!(work::NNLSGCVRegProblem, logμ; extract_subproblem = false)
    # Unpack buffers
    @unpack C, d, m, n, Aμ, C_buf, Ct_buf, CtC_buf = work

    # Solve regularized NNLS problem and record chi2 = ||C*X_reg - d||^2 which is returned
    μ = exp(logμ)
    solve!(work.nnls_work_smooth_cache, μ)
    χ² = chi2(get_cache(work.nnls_work_smooth_cache))
    x = solution(get_cache(work.nnls_work_smooth_cache))

    if extract_subproblem
        # Extract equivalent unconstrained least squares subproblem from NNLS problem
        # by extracting columns of C which correspond to nonzero components of x
        n′ = 0
        for (j,xⱼ) in enumerate(x)
            xⱼ ≈ 0 && continue
            n′ += 1
            @inbounds @simd ivdep for i in 1:m
                C_buf[i,n′] = C[i,j]
            end
        end
        C′ = reshape(uview(C_buf, 1:m*n′), m, n′)
        Ct′ = reshape(uview(Ct_buf, 1:n′*m), n′, m)
        CtC′ = reshape(uview(CtC_buf, 1:n′*n′), n′, n′)
    else
        # Use full matrix
        C′ = C
        Ct′ = Ct_buf
        CtC′ = CtC_buf
    end

    # Efficient compution of
    #   Aμ = C * (C'C + μ^2*I)^-1 * C'
    # where the matrices have sizes
    #   C: (m,n), Aμ: (m,m), Ct: (n,m), CtC: (n,n)
    mul!(CtC′, C′', C′) # C'C
    @inbounds @simd ivdep for i in 1:n
        CtC′[i,i] += μ^2 # C'C + μ^2*I
    end
    ldiv!(Ct′, cholesky!(Symmetric(CtC′)), C′') # (C'C + μ^2*I)^-1 * C'
    mul!(Aμ, C′, Ct′) # C * (C'C + μ^2*I)^-1 * C'

    # Return Generalized cross-validation. See equations 27 and 32 in
    #   Hansen, P.C., 1992. Analysis of Discrete Ill-Posed Problems by Means of the L-Curve. SIAM Review, 34(4), 561-580
    #   https://doi.org/10.1137/1034115
    trace = m - tr(Aμ) # tr(I - Aμ) = m - tr(Aμ) for m x m matrix Aμ
    gcv = χ² / trace^2 # ||C*X_reg - d||^2 / tr(I - Aμ)^2

    return gcv
end

# Equation (27) from Hansen et al. 1992 (https://epubs.siam.org/doi/10.1137/1034115),
# specialized for L = identity:
#
#   tr(I_m - A * (A'A + λ^2 * L'L)^-1 * A') = m - n + sum_i λ^2 / (γ_i^2 + λ^2)
#
# where γ_i are the generalized singular values, which are equivalent to ordinary
# singular values when L = identity, and size(A) = (m,n).
function gcv_tr!(A, λ)
    m, n = size(A)
    γ  = svdvals!(A)
    γ .= λ.^2 ./ (γ.^2 .+ λ.^2)
    max(m - n, 0) + sum(γ)
end
gcv_tr(A, λ) = gcv_tr!(copy(A), λ)
gcv_tr_brute(A, λ) = tr(I - A * ((A'A + λ^2 * I) \ A'))


text:
mutable struct Scaling
    dscale::Float64
    tscale::Float64
    uscale::Float64
    function Scaling()
        new(1e6,1e3,1e3)
    end
end

@inline function scale_rv(p::Scaling,r::SVector{3,T},v::SVector{3,T}) where T
    r = r/p.dscale
    v = v/(p.dscale/p.tscale)
    return r,v
end

@inline function unscale_rv(p::Scaling,r::SVector{3,T},v::SVector{3,T}) where T
    r = r*p.dscale
    v = v*(p.dscale/p.tscale)
    return r,v
end

@inline function scale_va(p::Scaling,v::SVector{3,T},a::SVector{3,T}) where T
    v = v/(p.dscale/p.tscale)
    a = a/(p.dscale/p.tscale^2)
    return v,a
end

@inline function unscale_va(p::Scaling,v::SVector{3,T},a::SVector{3,T}) where T
    v = v*(p.dscale/p.tscale)
    a = a*(p.dscale/p.tscale^2)
    return v,a
end

function unscale_X(p::Scaling,X::Vector{SVector{7,Float64}})
    N = length(X)
    Xu = [@SVector zeros(7) for i = 1:N]
    for i = 1:N
        r,v = unscale_rv(p,X[i][SA[1,2,3]],X[i][SA[4,5,6]])
        σ = X[i][7]
        Xu[i] = SVector{7}([r;v;σ])
    end
    return Xu
end

ext: jl


text:
using LOLTools.SpectatorV4

api_key = get(ENV, "RIOT_TOKEN", "")
platform = "kr"
featured = SpectatorV4.featured_games(api_key, platform)
# @info featured
@info sort(featured.gameList, by = x -> x.gameId)[1].participants


text:
# # Calibrate zonally-averaged channel from 3D eddying channel output

# ## Install dependencies

# ```julia
# using Pkg
# pkg"add ParameterEstimocean, Oceananigans, Distributions, CairoMakie"
# ```

using ParameterEstimocean
using ParameterEstimocean.Observations: FieldTimeSeriesCollector

using Oceananigans
using Oceananigans.Units
using Oceananigans.Models.HydrostaticFreeSurfaceModels: SliceEnsembleSize
using Oceananigans.TurbulenceClosures: FluxTapering

using LinearAlgebra, CairoMakie, DataDeps, Distributions, JLD2

using ElectronDisplay

architecture = CPU()

# filepath = "/Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/Ny200Nx200_Lx2000_Ly2000/eddying_channel_catke_zonal_average.jld2"
# filepath = "/Users/navid/Research/mesoscale-parametrization-OSM2022/eddying_channel/eddying_channel_convadj_zonal_time_average_1year.jld2"

filedir = @__DIR__
filename = "eddying_channel_convadj_zonal_time_average_1year.jld2"
filepath = joinpath(filedir, filename)
Base.download("https://www.dropbox.com/s/8wa05l2iqnbck1z/$filename", filepath)

file = jldopen(filepath)

# number of grid points
Nx, Ny, Nz = file["grid/Nx"], file["grid/Ny"], file["grid/Nz"]

# Domain
const Lx, Ly, Lz = file["grid/Lx"], file["grid/Ly"], file["grid/Lz"]

close(file)

grid = RectilinearGrid(architecture;
                       topology = (Periodic, Bounded, Bounded),
                       size = (Nx, Ny, Nz),
                       halo = (3, 3, 3),
                       x = (0, Lx),
                       y = (0, Ly),
                       z = (-Lz, 0)) # z_faces)

function get_field_timeseries(filepath, name, times)
    file = jldopen(filepath)

    iterations = parse.(Int, keys(file["timeseries/t"]))
    final_iteration = iterations[end]

    field = file["timeseries/$name/$final_iteration"]

    LX, LY, LZ = file["timeseries/$name/serialized/location"]

    close(file)

    field_timeseries = FieldTimeSeries{LX, LY, LZ}(grid, times)

    for n in 1:length(times)
        field_timeseries[n] .= field
    end

    return field_timeseries
end

end_time = 60day

times = [0, end_time]

u_timeseries = get_field_timeseries(filepath, "u", times)
v_timeseries = get_field_timeseries(filepath, "v", times)
η_timeseries = get_field_timeseries(filepath, "η", times)
b_timeseries = get_field_timeseries(filepath, "b", times)
c_timeseries = get_field_timeseries(filepath, "c", times)

#=
field_names = (:b, :c, :u, :v, :η)

transformation = (b = ZScore(),
                  c = ZScore(),
                  u = ZScore(),
                  v = RescaledZScore(1e-2),
                  η = RescaledZScore(1e-2))

field_time_serieses = (b = b_timeseries, c = c_timeseries, u = u_timeseries, v = v_timeseries, η = η_timeseries)
=#

# let's try fewer fields
field_names = (:b, :u)

transformation = (b = ZScore(),
                  u = ZScore())

field_time_serieses = (b = b_timeseries, u = u_timeseries)

observations = SyntheticObservations(; transformation, times, field_names, field_time_serieses)

#####
##### Simulation
#####

file = jldopen(filepath)
coriolis = file["serialized/coriolis"]
close(file)

Nensemble = 5
slice_ensemble_size = SliceEnsembleSize(size=(Ny, Nz), ensemble=Nensemble)

ensemble_grid = RectilinearGrid(architecture,
                                size = slice_ensemble_size,
                                topology = (Flat, Bounded, Bounded),
                                y = (0, Ly),
                                z = (-Lz, 0),
                                halo=(3, 3))

gent_mcwilliams_diffusivity = IsopycnalSkewSymmetricDiffusivity(slope_limiter = FluxTapering(1e-2))

file = jldopen(filepath)
closures = file["serialized/closure"]
close(file)

closure_ensemble = ([deepcopy(gent_mcwilliams_diffusivity) for _ = 1:Nensemble], closures[1], closures[2])

# closure_ensemble = [deepcopy(gent_mcwilliams_diffusivity) for k = 1:Nensemble]

ensemble_model = HydrostaticFreeSurfaceModel(grid = ensemble_grid,
                                             tracers = (:b, :c),
                                             buoyancy = BuoyancyTracer(),
                                             coriolis = coriolis,
                                             closure = closure_ensemble,
                                             free_surface = ImplicitFreeSurface())

Δt = 5minutes
simulation = Simulation(ensemble_model; Δt, stop_time=observations.times[end])

priors = (
    κ_skew = ScaledLogitNormal(bounds = (300, 1200)),
    κ_symmetric = ScaledLogitNormal(bounds = (300, 1200))
)

free_parameters = FreeParameters(priors)

calibration = InverseProblem(observations, simulation, free_parameters)

eki = EnsembleKalmanInversion(calibration; noise_covariance = 1e-2)

#=
collected_fields = (b = simulation.model.tracers.b,
                    w = simulation.model.velocities.w)
time_series_collector = FieldTimeSeriesCollector(collected_fields, observations.times)
calibration = InverseProblem(observations, simulation, free_parameters; time_series_collector)
eki = EnsembleKalmanInversion(calibration; noise_covariance = 1e-2)
=#

iterate!(eki; iterations = 5)

# Last, we visualize few metrics regarding how the EKI calibration went about.

θ̅(iteration) = [eki.iteration_summaries[iteration].ensemble_mean...]
varθ(iteration) = eki.iteration_summaries[iteration].ensemble_var

weight_distances = [norm(θ̅(iter)) for iter in 1:eki.iteration]
output_distances = [norm(forward_map(calibration, θ̅(iter))[:, 1] - y) for iter in 1:eki.iteration]
ensemble_variances = [varθ(iter) for iter in 1:eki.iteration]

f = Figure()
lines(f[1, 1], 1:eki.iteration, weight_distances, color = :red, linewidth = 2,
      axis = (title = "Parameter norm",
              xlabel = "Iteration",
              ylabel="|θ̅ₙ|",
              yscale = log10))
lines(f[1, 2], 1:eki.iteration, output_distances, color = :blue, linewidth = 2,
      axis = (title = "Output distance",
              xlabel = "Iteration",
              ylabel="|G(θ̅ₙ) - y|",
              yscale = log10))
ax3 = Axis(f[2, 1:2], title = "Parameter convergence",
           xlabel = "Iteration",
           ylabel = "Ensemble variance",
           yscale = log10)

for (i, pname) in enumerate(free_parameters.names)
    ev = getindex.(ensemble_variances, i)
    lines!(ax3, 1:eki.iteration, ev / ev[1], label = String(pname), linewidth = 2)
end

axislegend(ax3, valign = :top, halign = :right)
save("summary_channel.svg", f); nothing #hide

# ![](summary_channel.svg)

# And also we plot the the distributions of the various model ensembles for few EKI iterations to see
# if and how well they converge to the true diffusivity values.

f = Figure()

axtop = Axis(f[1, 1])

axmain = Axis(f[2, 1],
              xlabel = "κ_skew [m² s⁻¹]",
              ylabel = "κ_symmetric [m² s⁻¹]")

axright = Axis(f[2, 2])
scatters = []
labels = String[]

for iteration in [0, 1, 2, 5]
    ## Make parameter matrix
    parameters = eki.iteration_summaries[iteration].parameters
    Nensemble = length(parameters)
    Nparameters = length(first(parameters))
    parameter_ensemble_matrix = [parameters[i][j] for i=1:Nensemble, j=1:Nparameters]

    label = iteration == 0 ? "Initial ensemble" : "Iteration $iteration"
    push!(labels, label)
    push!(scatters, scatter!(axmain, parameter_ensemble_matrix))
    density!(axtop, parameter_ensemble_matrix[:, 1])
    density!(axright, parameter_ensemble_matrix[:, 2], direction = :y)
end

vlines!(axmain, [κ_skew], color = :red)
vlines!(axtop, [κ_skew], color = :red)

hlines!(axmain, [κ_symmetric], color = :red)
hlines!(axright, [κ_symmetric], color = :red)

colsize!(f.layout, 1, Fixed(300))
colsize!(f.layout, 2, Fixed(200))

rowsize!(f.layout, 1, Fixed(200))
rowsize!(f.layout, 2, Fixed(300))

Legend(f[1, 2], scatters, labels, valign = :bottom, halign = :left)

hidedecorations!(axtop, grid = false)
hidedecorations!(axright, grid = false)

xlims!(axmain, 350, 1350)
xlims!(axtop, 350, 1350)
ylims!(axmain, 650, 1750)
ylims!(axright, 650, 1750)
xlims!(axright, 0, 0.025)
ylims!(axtop, 0, 0.025)

save("distributions_channel.svg", f); nothing #hide

# ![](distributions_channel.svg)


text:
module TestDeprecated

using Test, DataFrames

const ≅ = isequal

@testset "by and aggregate" begin
    @test_throws ArgumentError by()
    @test_throws ArgumentError aggregate()
end

@testset "All indexing" begin
    df = DataFrame(a=1, b=2, c=3)

    @test select(df, All(1, 2)) == df[:, 1:2]
    @test select(df, All(1, :b)) == df[:, 1:2]
    @test select(df, All(:a, 2)) == df[:, 1:2]
    @test select(df, All(:a, :b)) == df[:, 1:2]
    @test select(df, All(2, 1)) == df[:, [2, 1]]
    @test select(df, All(:b, 1)) == df[:, [2, 1]]
    @test select(df, All(2, :a)) == df[:, [2, 1]]
    @test select(df, All(:b, :a)) == df[:, [2, 1]]

    @test df[:, All(1, 2)] == df[:, 1:2]
    @test df[:, All(1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2)] == df[:, 1:2]
    @test df[:, All(:a, :b)] == df[:, 1:2]
    @test df[:, All(2, 1)] == df[:, [2, 1]]
    @test df[:, All(:b, 1)] == df[:, [2, 1]]
    @test df[:, All(2, :a)] == df[:, [2, 1]]
    @test df[:, All(:b, :a)] == df[:, [2, 1]]

    @test df[:, All(1, 1, 2)] == df[:, 1:2]
    @test df[:, All(:a, 1, :b)] == df[:, 1:2]
    @test df[:, All(:a, 2, :b)] == df[:, 1:2]
    @test df[:, All(:a, :b, 2)] == df[:, 1:2]
    @test df[:, All(2, 1, :a)] == df[:, [2, 1]]

    @test select(df, All(1, "b")) == df[:, 1:2]
    @test select(df, All("a", 2)) == df[:, 1:2]
    @test select(df, All("a", "b")) == df[:, 1:2]
    @test select(df, All("b", 1)) == df[:, [2, 1]]
    @test select(df, All(2, "a")) == df[:, [2, 1]]
    @test select(df, All("b", "a")) == df[:, [2, 1]]

    @test df[:, All(1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2)] == df[:, 1:2]
    @test df[:, All("a", "b")] == df[:, 1:2]
    @test df[:, All("b", 1)] == df[:, [2, 1]]
    @test df[:, All(2, "a")] == df[:, [2, 1]]
    @test df[:, All("b", "a")] == df[:, [2, 1]]

    @test df[:, All("a", 1, "b")] == df[:, 1:2]
    @test df[:, All("a", 2, "b")] == df[:, 1:2]
    @test df[:, All("a", "b", 2)] == df[:, 1:2]
    @test df[:, All(2, 1, "a")] == df[:, [2, 1]]

    df = DataFrame(a1=1, a2=2, b1=3, b2=4)
    @test df[:, All(r"a", Not(r"1"))] == df[:, [1, 2, 4]]
    @test df[:, All(Not(r"1"), r"a")] == df[:, [2, 4, 1]]
end

@testset "indicator in joins" begin
    name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
    job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

    @test outerjoin(name, job, on = :ID, indicator=:source) ≅
          outerjoin(name, job, on = :ID, source=:source)
    @test leftjoin(name, job, on = :ID, indicator=:source) ≅
          leftjoin(name, job, on = :ID, source=:source)
    @test rightjoin(name, job, on = :ID, indicator=:source) ≅
          rightjoin(name, job, on = :ID, source=:source)

    @test_throws ArgumentError outerjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
    @test_throws ArgumentError leftjoin(name, job, on = :ID,
                                       indicator=:source, source=:source)
    @test_throws ArgumentError rightjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
end

end # module


text:
include("./imports.jl")

@testset ExtendedTestSet "Normalize" begin
    @testset ExtendedTestSet "normalize" begin
        image = rand(Float32, 10, 10, 3)
        means = SVector{3}(rand(Float32, 1, 1, 3))
        stds = SVector{3}(rand(Float32, 1, 1, 3))
        @test denormalize(normalize(copy(image), means, stds), means, stds) ≈ image
    end

    @testset ExtendedTestSet "Normalize" begin
        item = ArrayItem(zeros(10, 10, 3))
        means =  SVector{3}(ones(1, 1, 3))
        stds =  SVector{3}(ones(1, 1, 3))
        tfm = Normalize(means, stds)
        @test itemdata(apply(tfm, item)) ≈ -1 .* ones(10, 10, 3)
    end
end

@testset ExtendedTestSet "NormalizeIntensity" begin
    @testset ExtendedTestSet "NormalizeIntensity" begin
        item = ArrayItem(Float32.(collect(1:5)))
        ground_truth = [-1.2649, -0.6325, 0, 0.6325, 1.2649]
        tfm = NormalizeIntensity()
        @test itemdata(apply(tfm, item)) ≈ ground_truth
    end

    @testset ExtendedTestSet "NormalizeIntensity" begin
    item = Image(Float32.(collect(1:5)))
    ground_truth = [-1.2649, -0.6325, 0, 0.6325, 1.2649]
    tfm = NormalizeIntensity()
    @test itemdata(apply(tfm, item)) ≈ ground_truth
    end
end

@testset ExtendedTestSet "AddChannel" begin
    @testset ExtendedTestSet "AddChannel" begin
    item = Image(rand(Float32, 20, 20, 20, 2))
    tfm = AddChannel()
    @test size(itemdata(apply(tfm, item))) == (20, 20, 20, 2, 1)
    end

    @testset ExtendedTestSet "AddChannel" begin
    item = Image(rand(Float32, 20, 2))
    tfm = AddChannel()
    @test size(itemdata(apply(tfm, item))) == (20, 2, 1)
    end

    @testset ExtendedTestSet "AddChannel" begin
    item = MaskBinary(rand(Bool, 20, 2))
    tfm = AddChannel()
    @test size(itemdata(apply(tfm, item))) == (20, 2, 1)
    end
end

@testset ExtendedTestSet "ToEltype" begin
    tfm = ToEltype(RGB)
    data = rand(RGB, 10, 10)
    imagergb = Image(data)
    imagegray = Image((Gray.(data)))
    @test itemdata(apply(ToEltype(Gray), imagegray)) === itemdata(imagegray)
    @test itemdata(apply(ToEltype(Gray), imagergb)) ≈ itemdata(imagegray)
end

@testset ExtendedTestSet "ImageToTensor" begin
    @testset ExtendedTestSet "imagetotensor,tensortoimage" begin
        data = rand(RGB, 10, 10)
        @test size(imagetotensor(data)) == (10, 10, 3)
        @test data |> imagetotensor |> tensortoimage ≈ data

    end

    @testset ExtendedTestSet "ImageToTensor" begin
        image = Image(rand(RGB, 10, 10))
        tfm = ImageToTensor()
        @test_nowarn apply(tfm, image)
        a = itemdata(apply(tfm, image))
        @test size(a) == (10, 10, 3)
    end
end

@testset ExtendedTestSet "OneHot" begin
    tfm = OneHot()
    mask = rand(1:4, 10, 10)
    item = MaskMulti(mask, 1:4)
    @test_nowarn apply(tfm, item)
    aitem = apply(tfm, item)
    @test size(itemdata(aitem)) == (10, 10, 4)

    item2 = MaskMulti(rand(1:3, 10, 10), 1:4)
    buf = itemdata(aitem)
    bufcopy = copy(buf)
    apply!(aitem, tfm, item2)
    @test itemdata(item) == itemdata(item2) || itemdata(aitem) != bufcopy

end

@testset ExtendedTestSet "Image pipeline" begin
    image = Image(rand(RGB, 150, 150))

    tfm = compose(
        ToEltype(RGB),
        RandomResizeCrop((128, 64)),
        ImageToTensor(),
        Normalize([0, 0, 0], [1, 1, 1])
    )

    @test_nowarn apply(tfm, image)

    res = apply(tfm, image)
    a = itemdata(res)
    @test size(a) == (128, 64, 3)
    @test eltype(a) == Float32
end


text:
import Base: isvalid, iterate, show

struct ISBN <: AbstractString
    string::AbstractString

    function ISBN(string::AbstractString)
        if !isvalid(ISBN, string)
            throw(DomainError("Invalid ISBN"))
        end
        string = replace(string, r"-" => "")
        new(string)
    end
end

macro isbn_str(s) ISBN(s) end

iterate(isbn::ISBN) = iterate(isbn.string)

iterate(isbn::ISBN, i::Integer) = iterate(isbn.string, i)

show(io::IO, isbn::ISBN) =
    show(io, replace(isbn.string, r"(.)(.{3})(.{5})(.)" => s"\1-\2-\3-\4"))

function isvalid(::Type{ISBN}, string::String)
    string = replace(string, r"-" => "")

    if length(string) != 10
        return false
    end

    tot = 0
    for (i, ch) in enumerate(string)
        if i == 10 && ch == 'X'
            x = 10
        elseif ch in '0':'9'
            x = Int(ch) - 48
        else
            return false
        end

        tot += x * i
    end

    return tot % 11 == 0
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
#This defines the tree structure.
mutable struct Tree{String <: A <: String, Vector{Int64}<:B <: Vector{Int64}, Vector{Vector{Int64}} <:C<: Vector{Vector{Int64}}, D<: Vector{Float64}}
    name::String                 # name of the tree
    parent::B               # parents of nodes in the tree
    children::C             # successor nodes of each parent
    state::D                # states of nodes in the tree
    probability::D          # probability to go from one node to another

    """
    	child(parent::Vector{Int64})

    Returns a vector of successors of each of the parent nodes.
    """
    function child(parent::B) where {B}
        allchildren = Vector{B}([])
        for node in unique(parent)
            push!(allchildren, [i for i = 1 : length(parent) if parent[i] == node])
        end
        return allchildren
    end

    Tree(name::A,parent::B, children::C,
    state::D, probability::D) where {A,B,C,D} =new{A, B,C,D}(name, parent, child(parent), state, probability)
    Tree() = Tree("", Int64[], child(Int64[]), zeros(Float64, 0), ones(Float64, 0))

    """
    	Tree(identifier::Int64)

    Returns some examples of predefined trees.
    - These are (0, 302, 303, 304, 305, 306, 307, 402, 404, 405).
    - You can call any of the above tree and plot to see the properties of the tree.
    """
    function Tree(identifier::Int64)
        if identifier==0
            name = "Empty Tree"
            parent = zeros(Int64 , 0)
            children = child(parent)
            state = zeros(Float64, 0)
            probability = ones(Float64, 0)
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 302
            name = "Tree 1x2x2"
            parent =[0,1,1,2,2,3,3]
            children = child(parent)
            state = vec([2.0 2.1 1.9 4.0 1.0 3.0])
            probability = vec([1.0 0.5 0.5 0.5 0.5 0.5 0.5])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 303
            name = "Tree 1x1x4"
            parent = [0,1,2,2,2,2]
            children = child(parent)
            state = vec([3.0 3.0 6.0 4.0 2.0 0.0])
            probability = vec([1.0 1.0 0.25 0.25 0.25 0.25])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 304
            name = "Tree 1x4x1x1"
            parent = [0,1,2,0,4,5,0,7,8,0,10,11]
            children = child(parent)
            state = vec([0.1 2.1 3.0 0.1 1.9 1.0 0.0 -2.9 -1.0 -0.1 -3.1 -4.0])
            probability = vec([0.14 1.0 1.0 0.06 1.0 1.0 0.48 1.0 1.0 0.32 1.0 1.0])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 305
            name = "Tree 1x1x4"
            parent = [0,1,2,2,2,2]
            children = child(parent)
            state = vec([0.0 10.0 28.0 22.0 21.0 20.0])
            probability = vec([1.0 1.0 0.25 0.25 0.25 0.25])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 306
            name = "Tree 1x2x2"
            parent = [0,1,1,2,2,3,3]
            children = child(parent)
            state = vec([0.0 10.0 10.0 28.0 22.0 21.0 20.0])
            probability = vec([1.0 0.5 0.5 0.5 0.5 0.5 0.5])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 307
            name = "Tree 1x4x1"
            parent = [0,1,1,1,1,2,3,4,5]
            children = child(parent)
            state = vec([0.0 10.0 10.0 10.0 10.0 28.0 22.0 21.0 20.0])
            probability = vec([1.0 0.25 0.25 0.25 0.25 1.0 1.0 1.0 1.0])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 401
            name = "Tree 1x1x2x2"
            parent = [0,1,2,2,3,3,4,4]
            children = child(parent)
            state = vec([10.0 10.0 8.0 12.0 9.0 6.0 10.0 13.0])
            probability = vec([1.0 1.0 0.66 0.34 0.24 0.76 0.46 0.54])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 402
            name = "Tree 1x2x2x2"
            parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            children = child(parent)
            state = vec([10.0 12.0 8.0 15.0 11.0 9.0 5.0 18.0 16.0 13.0 11.0 10.0 7.0 6.0 3.0])
            probability = vec([1.0 0.8 0.2 0.3 0.7 0.8 0.2 0.6 0.4 0.5 0.5 0.4 0.6 0.7 0.3])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier ==404
            name = "Tree 1x2x2x2"
            parent = [0,1,1,2,2,3,3,4,4,5,5,6,6,7,7]
            children = child(parent)
            state = (0.2+0.6744).*vec([0.0 1.0 -1.0 2.0 0.1 0.0 -2.0 3.0 1.1 0.9 -1.1 1.2 -1.2 -0.8 -3.2])
            probability = vec([1.0 0.3 0.7 0.2 0.8 0.1 0.9 0.5 0.5 0.6 0.4 0.4 0.6 0.3 0.7])
            self = Tree(name, parent, child(parent), state, probability)
        elseif identifier == 405
            name = "Tree 5"
            parent = [0,1,2,3,3,2,6,6,6,1,10,10,11,12,12,12]
            children = child(parent)
            state = vec([7.0 12.0 14.0 15.0 13.0 9.0 11.0 10.0 8.0 4.0 6.0 2.0 5.0 0.0 1.0 3.0])
            probability = vec([1.0 0.7 0.4 0.7 0.3 0.6 0.2 0.3 0.5 0.3 0.2 0.8 1.0 0.3 0.5 0.2])
            self = Tree(name, parent, child(parent), state, probability)
        end
        return self
    end

    """
    	Tree(bs::Vector{Int64})

    Returns a tree according to the specified branching vector and the dimension.

    Args:
    - bs - the branching structure of the scenario tree.

    The branching vector must start with 1 for the root node.
    """
    function Tree(bs::Vector{B}) where {B}
        #self = Tree("Tree $(bstructure)",  )
        self = Tree()
        self.name = "Tree $(bs)"

        leaves = cumprod(bs[1:end-1]) #probably not needed, as immediately overwritten in loop
        self.parent = Int32.([0])
        self.state = [1.0]
        self.children = [[1]]
        idx = 2
        for stage = 2 : length(bs)
            newleaves = vec(ones(Int64, bs[stage]) .* transpose(length(self.parent) .+ (1 .- leaves[stage-1] : 0)))
            self.parent =  vcat(self.parent, newleaves)
            for j=1:cumprod(bs)[stage-1]
                push!(self.children,collect(idx:idx+bs[stage]-1) ); idx += bs[stage];
            end
        end

        #self.children = child(self.parent)
        self.state = ones(length(self.parent));
        self.probability = similar(self.state);
        return self
    end
end

function children(parent)
    allchildren = Vector([])
    for node in unique(parent)
        push!(allchildren, [i for i = 1 : length(parent) if parent[i] == node])
    end
    return allchildren
end

"""
	stage(trr::Tree, node=Int64[])

Returns the stage of each node in the tree.

# Arguments:
- trr - an instance of a Tree.
- node - the number of node in the scenario tree you want to know its stage.
"""
function stage(trr::Tree{A,B,C,D}  , node = Int64[]) where {A,B,C,D}
    if isempty(node)
        node = 1 : length(trr.parent)
    elseif isa(node, Int64)
        node = Int64[node]
    end
    #stage = zeros(length(collect(node))) #you can also use this
    stage = zero(node)
    for i = 1 : length(node)
        pred = node[i]
        while pred > 0 && trr.parent[pred] > 0
            pred = trr.parent[pred]
            stage[i] += 1
        end
    end
    return stage
end

"""
	height(trr::Tree)

Returns the height of the tree `trr` which is just the maximum number of the stages of each node.
"""
function height(trr::Tree{A,B,C,D}) where {A,B,C,D}
    return maximum(stage(trr))
end

"""
	leaves(trr::Tree,node=Int64[])

Returns the leaf nodes, their indexes and the conditional probabilities based on parent `node`.

# Arguments:
- trr - an instance of a Tree.
- node - a node for which all leaf nodes are found.

# Values
- leaves - a vector of leaf nodes.
- omegas - a vector of indizes indicating the order of the elements in `leaves`.
- prob - the conditional probability of each leaf node conditional on `node`.

If `node` is not given, then `node` is treated as the root node.

"""
function leaves(trr::Tree{A,B,C,D}, node = Int64[]) where {A,B,C,D}
    nodes = 1 : length(trr.parent)
    leaves = setdiff(nodes, trr.parent)
    #leaves are all nodes which are not parents
    #setdiff(A,B) finds all the elements that belongs to A and are not in B
    omegas = 1 : length(leaves)
    if !isempty(node) && isa(node,Int64)
        node = Int64[node]
    end
    if !isempty(node) && (0 ∉ node)
        omegas = Set{Int64}()
        nodes = leaves
        while any(j !=0 for j ∈ nodes)
            omegas = union(omegas, (ind for (ind, j) ∈ enumerate(nodes) if j ∈ node))
            nodes = Int64[trr.parent[max(0 , j)] for j ∈ nodes]
        end
        omegas = collect(omegas)
    end
    leaves = Int64[leaves[j] for j ∈ omegas]
    prob = ones(Float64 , length(leaves))
    nodes = leaves
    while any(j !=0 for j ∈ nodes)
        prob = prob .* trr.probability[[j for j in nodes]]
        nodes = Int64[trr.parent[max(0 , j)] for j ∈ nodes]
    end
    return leaves, omegas, prob
end

"""
	nodes(trr::Tree,t=Int64[])

Returns the nodes in the tree `trr` at stage `t`.

# Example:

```julia-repl
julia> trr = Tree(404);
julia> nodes(trr,1)
2-element Vector{Int64}:
 2
 3
```
"""
function nodes(trr::Tree{A,B,C,D}, t = Int64[]) where {A,B,C,D}
    nodes = 1 : length(trr.parent)
    if isempty(t) #if stage t is not given, return all nodes of the tree
        return nodes
    else # else return nodes at the given stage t
        stg = stage(trr)
        return Int64[i for i in nodes if stg[i] == t]
    end
end

"""
	root(trr::Tree,nodes=Int64[])
Returns the root of the tree if the node is not specified.

#Arguments:
- trr - an instance of Tree.
- nodes - target node from root.

If `nodes` is not specified, it returns the root of the tree.

If `nodes` is specified, it returns a sequence of nodes from the root to the specified node.
"""
function root(trr::Tree{A,B,C,D}, nodes = Int64[]) where {A,B,C,D}
    if isempty(nodes)
        nodes = trr.children[1]
    elseif isa(nodes , Int64)
        nodes = Int64[nodes]
    end
    root = 1 : length(trr.parent)
    for i in nodes
        iroot = Vector{Int64}([])       # parameterize with B
        tmp = i
        while tmp > 0
            push!(iroot , tmp)
            tmp = trr.parent[tmp]
        end
        root = Int64[i for i in reverse(iroot) if i in root]
    end
    return root
end

"""
	part_tree(trr::Tree,node)

Returns a reduced tree starting at `node`. If `node` is the root node then original tree is returned.

# Arguments:
- trr - an instance of Tree
- node - a node of trr which becomes the new root node.
"""
function part_tree(trr, node::Int64=0)
    if node == 0 || node == 1
        return(trr)
    end
    trr2 = deepcopy(trr)
    trr2.state[1] = trr2.state[node]
    tmp = reverse(unique(trr2.parent)); pop!(tmp)
    for i in tmp
        if !(node ∈ root(trr2,[i]))
            # delete i from parents, children and states
            deleteat!(trr2.state, i .∈ trr2.parent)
            deleteat!(trr2.probability, i .∈ trr2.parent)
            deleteat!(trr2.parent, i .∈ trr2.parent)
        end
    end

    tmp = unique(trr2.parent[2:end])
    for j=1:length(tmp)
        trr2.parent[trr2.parent .== tmp[j]] .= j
    end

    trr2.children = children(trr2.parent)
    return(trr2)

end

"""
	combine_tree(target:Tree,trr::Tree, trrs)

combines a base tree `trr` with a list of trees in `trrs` such that at every
leaf node of `trr` a tree in `trrs` is added. In a sense this reverses `part_tree`.

# Arguments:
- target - an (empty) tree providing the overall tree structure for the combined tree.
- trr - an instance of Tree
- trrs - a vector of trees. The length of trrs must be equal to the number of leaf nodes in trr.

The target tree is used here to make the algorithm easier but strictly speaking not necessary and
will be deleted in future versions.

#Example
trr = Tree(404)
trrs = [];
for i=1:length(nodes(trr,2))
    push!(trrs, part_tree(trr,nodes(trr,2)[i]))
end
target2 = deepcopy(trr); target2.state .= 0.0;
trr.state=trr.state[1:7];
trr.probability=trr.probability[1:7];
trr.parent=trr.parent[1:7];
trr.children=trr.children[1:7];

target3 = combine_tree(target2,trr,trrs)

"""

##NOT WORKING!!! trr must be what is deleted by part_tree....
function combine_tree(target, trr, trrs)
    @assert length(trrs) == length(leaves(trr)[1])

    newStates = [];
    newProb  = [];
    # first stage of new trrs not needed as equal to leaf nodes of trr
    for i=1:height(trrs[1])  # for every additional stage
        stage = i
        for j=1:length(trrs)        # for every tree to be added
            tmpNodes = nodes(trrs[j],stage)
            tmpStates = trrs[j].state[tmpNodes]; # trrs[j].states for states in correct stage
            tmpProb = trrs[j].probability[tmpNodes]
            append!(newStates, tmpStates);
            append!(newProb, tmpProb)
        end
    end

    println(length(newStates))
    println(length(trr.state))
    target.state = append!(trr.state,newStates)
    target.probability = append!(trr.probability,newProb)

    return(target)
end

"""
	build_probabilities!(trr::Tree,probabilities::Array{Float64,2})

Returns the probabilities of the nodes without probabilities
if the array of probabilities is less than the length of parents
in the stochastic approximation procedure.
"""
function build_probabilities!(trr::Tree{A,B,C,D}, probabilities::Array{Float64,2}) where {A,B,C,D}
    leaf, omegas, probaLeaf = leaves(trr)
    if length(probabilities) == length(nodes(trr))
        trr.probability .= probabilities
    else
        probabilities .= max.(0.0,probabilities)
        j = leaf; i = trr.parent[j]; j .= j[i .> 0]; i .= i[i .> 0]
        trr.probability .= zeros(length(trr.state), 1)[:,1]
        trr.probability[j] = probabilities
        while !isempty(i)
            for k = 1:length(i)
                trr.probability[i[k]] = trr.probability[i[k]] + trr.probability[j[k]]
            end
            trr.probability[j] .= trr.probability[j] ./ trr.probability[i]
            trr.probability[isnan.(trr.probability)] .= 0.0
            j = unique(i); i = trr.parent[j]; j = j[i .> 0] ; i = i[i .> 0]
        end
    end
    return trr.probability
end

# Checks if probabilities are strictly positive and if not adjusts slightly.
# What adjustments can be justified have to be determined in each instance.
"""
    checkTree(tree0::Tree{A,B,C,D}, eps) where {A,B,C,D}

Transitions with zero probability are given positive probability `eps` such that conditional probabilities sum up to 1.0.

The correct adjustment can only be justified for each individual tree. Be careful in using this function as it can drastically alter
the tree process described by `tree0`.

"""
function checkTree(trr::Tree{A,B,C,D}, eps) where {A,B,C,D}

    # For every parent node going backwards
    tmp = reverse(unique(trr.parent)); pop!(tmp)
    for i in tmp
        #Note that the probs are the conditional ones
        prob = trr.probability[trr.children[i+1]]
        if any(x->x>0,prob)
          #do nothing
        else
          prob = prob .+ eps
          prob = prob./ sum(prob)
          trr.probability[trr.children[i+1]] = prob
        end

     end

end

"""
    cumulProb!(trr::Tree{A,B,C,D}) where {A,B,C,D}

modifies the probability vector of `trr` to show cumulative probabilities instead of conditional ones.
"""
function cumulProb!(trr::Tree{A,B,C,D}) where {A,B,C,D}
    for (index, node) in enumerate(trr.parent)
        if node==0
            #do nothing
        else
            trr.probability[index] = trr.probability[index] * trr.probability[node]
        end
    end
    return(trr)
end


text:
using WorkshopWizard
using Test, Pkg

function with_temp_env(f, env_name::AbstractString = "Dummy"; rm = true)
    prev_active = Base.ACTIVE_PROJECT[]
    env_path = joinpath(mktempdir(), env_name)
    Pkg.generate(env_path)
    Pkg.activate(env_path)
    if length(LOAD_PATH) == 3
        deleteat!(LOAD_PATH, 2)
        deleteat!(LOAD_PATH, 2)
    end
    try
        applicable(f, env_path) ? f(env_path) : f()
    finally
        Base.ACTIVE_PROJECT[] = prev_active
        try
            rm && Base.rm(env_path; force = true, recursive = true)
        catch err
            # Avoid raising an exception here as it will mask the original exception
            println(
                Base.stderr,
                "Exception in finally: $(sprint(showerror, err))",
            )
        end
        push!(LOAD_PATH, "@v#.#")
        push!(LOAD_PATH, "@stdlib")
    end
end

function test_load_pkg(pkg::Symbol)
    try
        @eval import $pkg
        return true
    catch
        return false
    end
end

function rm_global_IJulia()
    WorkshopWizard.with_pkg_env("", globalenv = true) do
        pkg"rm IJulia"
    end
end

@testset "WorkshopWizard.jl" begin

    @testset "Defaults" begin

        _, latest = findmax(map(
            w -> parse(Int, w[end-1:end]),
            WorkshopWizard.WORKSHOPS,
        ))
        latest_workshop = WorkshopWizard.WORKSHOPS[latest]
        @test WorkshopWizard.default_workshop() == latest_workshop
        @test WorkshopWizard.default_repo() == "https://github.com/crstnbr/$(latest_workshop)"
        if Sys.iswindows()
            @test WorkshopWizard.default_path() == joinpath(
                homedir(),
                "Desktop",
            )
        elseif Sys.islinux()
            @test WorkshopWizard.default_path() == homedir()
        else
            @test WorkshopWizard.default_path() == homedir()
        end
    end

    @testset "IJulia" begin
        with_temp_env() do
            @test WorkshopWizard._check_IJulia() == false
            WorkshopWizard.install_IJulia()
            @test WorkshopWizard._check_IJulia() == true
            with_temp_env() do
                @test WorkshopWizard._check_IJulia() == true
            end
            WorkshopWizard.uninstall_IJulia()
            @test WorkshopWizard._check_IJulia() == false
        end
    end

    @testset "Download" begin
        cd(mktempdir()) do
            WorkshopWizard.download(
                repo = "https://github.com/crstnbr/JuliaTestWorkshop",
                path = pwd(),
            )
            @test isdir("JuliaTestWorkshop")
            @test isfile("JuliaTestWorkshop/README.md")
        end
    end

    @testset "Install" begin
        cd(mktempdir()) do
            WorkshopWizard.install(
                repo = "https://github.com/crstnbr/JuliaTestWorkshop",
                path = pwd(),
            )
            @test isdir("JuliaTestWorkshop")
            @test isfile("JuliaTestWorkshop/README.md")
            @test test_load_pkg(:IJulia)
            WorkshopWizard.with_pkg_env("JuliaTestWorkshop") do
                @test test_load_pkg(:IJulia)
                @test test_load_pkg(:Colors)
                @test test_load_pkg(:BenchmarkTools)
                @test test_load_pkg(:GenericLinearAlgebra)
            end
            WorkshopWizard.uninstall_IJulia()
        end

        # without global IJulia
        cd(mktempdir()) do
            WorkshopWizard.install(
                repo = "https://github.com/crstnbr/JuliaTestWorkshop",
                path = pwd(),
                global_IJulia = false,
            )
            @test isdir("JuliaTestWorkshop")
            @test isfile("JuliaTestWorkshop/README.md")
            @test test_load_pkg(:IJulia) == false
            WorkshopWizard.with_pkg_env("JuliaTestWorkshop") do
                @test test_load_pkg(:IJulia) # assuming that JuliaTestWorkshop depends on IJulia
                @test test_load_pkg(:Colors)
                @test test_load_pkg(:BenchmarkTools)
                @test test_load_pkg(:GenericLinearAlgebra)
            end
            WorkshopWizard.uninstall_IJulia()
        end
    end

end # main


text:
### HMC holds the fields of the HMC sampler
### These fields represent the initial user-defined state of the sampler
### These sampler fields are copied to the corresponding heap type fields, where the latter can be tuned

immutable HMC <: HMCSampler
  nleaps::Int # Number of leapfrog steps
  leapstep::Float64 # Leapfrog stepsize

  function HMC(nl::Int, l::Float64)
    @assert nl > 0 "Number of leapfrog steps is not positive."
    @assert l > 0 "Leapfrog stepsize is not positive."
    new(nl, l)
  end
end

HMC(nl::Int) = HMC(nl, 0.1)
HMC(l::Float64) = HMC(10, l)

HMC(; nleaps::Int=10, leapstep::Float64=0.1) = HMC(nleaps, leapstep)

### Sample type for HMC

type HMCSample <: MCSample{FirstOrder}
  sample::Vector{Float64} # sample position, i.e. model parameters
  logtarget::Float64
  gradlogtarget::Vector{Float64}
  momentum::Vector{Float64}
  hamiltonian::Float64
end

HMCSample(s::Vector{Float64}) = HMCSample(s, NaN, Float64[], Float64[], NaN)
HMCSample() = HMCSample(Float64[], NaN, Float64[], Float64[], NaN)

HMCSample(l::Int) = HMCSample(fill(NaN, l), NaN, fill(NaN, l), fill(NaN, l), NaN)

hamiltonian!(s::HMCSample) = (s.hamiltonian = -s.logtarget+0.5*Base.dot(s.momentum, s.momentum))

function leapfrog(s::HMCSample, f::Function, leapstep::Float64)
  lsample = deepcopy(s)
  lsample.momentum += 0.5*lsample.gradlogtarget*leapstep
  lsample.sample += leapstep*lsample.momentum
  gradlogtargetall!(lsample, f)
  lsample.momentum += 0.5*lsample.gradlogtarget*leapstep
  hamiltonian!(lsample)
  lsample
end

### HMCHeap type holds the internal state ("local variables") of the HMC sampler

type HMCHeap <: MCHeap{HMCSample}
  instate::MCState{HMCSample} # Monte Carlo state used internally by the sampler
  outstate::MCState{HMCSample} # Monte Carlo state outputted by the sampler
  tune::MCTune
  count::Int # Current number of iterations
  nleaps::Int # Number of leapfrog steps for the current Monte Carlo iteration (possibly tuned)
  leapstep::Float64 # Leapfrog stepsize for the current Monte Carlo iteration (possibly tuned)
end

HMCHeap() = HMCHeap(MCState(HMCSample(), HMCSample()), MCState(HMCSample(), HMCSample()), VanillaMCTune(), 0, 0, NaN)

HMCHeap(l::Int, t::MCTune=VanillaMCTune()) =
  HMCHeap(MCState(HMCSample(l), HMCSample(l)), MCState(HMCSample(l), HMCSample(l)), t, 0, 0, NaN)

### Initialize HMC sampler

function initialize_heap(m::MCModel, s::HMC, r::MCRunner, t::MCTuner)
  @assert hasgradient(m) "HMC sampler requires model with gradient function."
  heap::HMCHeap = HMCHeap(m.size)

  heap.instate.current = HMCSample(copy(m.init))
  gradlogtargetall!(heap.instate.current, m.evalallg)
  @assert isfinite(heap.instate.current.logtarget) "Initial values out of model support."

  if isa(t, VanillaMCTuner)
    heap.tune = VanillaMCTune()
  elseif isa(t, EmpiricalMCTuner)
    heap.tune = EmpiricalMCTune(s.leapstep, s.nleaps)
  end

  heap.count = 1

  heap
end

function reset!(heap::HMCHeap, x::Vector{Float64}, m::MCModel)
  heap.instate.current = HMCSample(copy(x))
  gradlogtargetall!(heap.instate.current, m.evalallg)
end

function initialize_task!(heap::HMCHeap, m::MCModel, s::HMC, r::MCRunner, t::MCTuner)
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, (x::Vector{Float64})->reset!(heap, x, m))

  while true
    iterate!(heap, m, s, r, t, produce)
  end
end

### Perform iteration for HMC sampler

function iterate!(heap::HMCHeap, m::MCModel, s::HMC, r::MCRunner, t::MCTuner, send::Function)
  if isa(t, VanillaMCTuner)
    heap.nleaps, heap.leapstep = s.nleaps, s.leapstep
    if t.verbose
      heap.tune.proposed += 1
    end
  elseif isa(t, EmpiricalMCTuner)
    heap.nleaps, heap.leapstep = heap.tune.nsteps, heap.tune.step
    heap.tune.proposed += 1
  end

  heap.instate.current.momentum = randn(m.size)
  hamiltonian!(heap.instate.current)
  heap.instate.successive = deepcopy(heap.instate.current)

  for j = 1:heap.nleaps
    heap.instate.successive = leapfrog(heap.instate.successive, m.evalallg, heap.leapstep)
  end

  if rand() < exp(heap.instate.current.hamiltonian-heap.instate.successive.hamiltonian)
    heap.outstate = MCState(heap.instate.successive, heap.instate.current, Dict{Any, Any}("accept" => true))
    heap.instate.current = deepcopy(heap.instate.successive)

    if isa(t, VanillaMCTuner) && t.verbose
      heap.tune.accepted += 1
    elseif isa(t, EmpiricalMCTuner)
      heap.tune.accepted += 1
    end
  else
    heap.outstate = MCState(heap.instate.current, heap.instate.current, Dict{Any, Any}("accept" => false))
  end

  if isa(t, VanillaMCTuner) && t.verbose && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    rate!(heap.tune)
    println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
  elseif isa(t, EmpiricalMCTuner) && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    adapt!(heap.tune, t)
    reset!(heap.tune)
    if t.verbose
      println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
    end
  end

  heap.count += 1

  send(heap.outstate)
end


text:
# Lies!  Lies and Slander!  Supposedly ccache doesn't make things fail, but
# WHAT DO WE HAVE?! Literal proof that this is not the case.  This whole thing
# fails to do the proper symbol renaming with `ccache`, which is TERRIBLE.
ENV["BINARYBUILDER_USE_CCACHE"] = "false"
using BinaryBuilder

# Collection of sources required to build Arpack
name = "Arpack"
version = v"3.7.0"
sources = [
    "https://github.com/opencollab/arpack-ng.git" =>
    "a580f6fc103c4bba1630ddd957b0b5f58dec284c",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
mkdir arpack-build
cd arpack-build

# arpack tests require finding libgfortran when linking with C linkers,
# and gcc doesn't automatically add that search path.  So we do it for it.
EXE_LINK_FLAGS="${LDFLAGS}"
if [[ ${target} != *darwin* ]]; then
    EXE_LINK_FLAGS="${EXE_LINK_FLAGS} -Wl,-rpath-link,/opt/${target}/${target}/lib -Wl,-rpath-link,/opt/${target}/${target}/lib64"
fi

# Symbols that have float32, float64, complexf32, and complexf64 support
SDCZ_SYMBOLS="axpy copy gemv geqr2 lacpy lahqr lanhs larnv lartg \
              lascl laset scal trevc trmm trsen gbmv gbtrf gbtrs \
              gttrf gttrs pttrf pttrs"

# All symbols that have float32/float64 support (including the SDCZ_SYMBOLS above)
SD_SYMBOLS="${SDCZ_SYMBOLS} dot ger labad laev2 lamch lanst lanv2 \
            lapy2 larf larfg lasr nrm2 orm2r rot steqr swap"

# All symbols that have complexf32/complexf64 support (including the SDCZ_SYMBOLS above)
CZ_SYMBOLS="${SDCZ_SYMBOLS} dotc geru unm2r"

# Add in (s|d)*_64 symbol remappings:
for sym in ${SD_SYMBOLS}; do
    SYMBOL_DEFS="${SYMBOL_DEFS} -Ds${sym}=s${sym}_64 -Dd${sym}=d${sym}_64"
done

# Add in (c|z)*_64 symbol remappings:
for sym in ${CZ_SYMBOLS}; do
    SYMBOL_DEFS="${SYMBOL_DEFS} -Dc${sym}=c${sym}_64 -Dz${sym}=z${sym}_64"
done

# Add one-off symbol mappings; things that don't fit into any other bucket:
for sym in scnrm2 dznrm2 csscal zdscal dgetrf dgetrs; do
    SYMBOL_DEFS="${SYMBOL_DEFS} -D${sym}=${sym}_64"
done

# Set up not only lowercase symbol remappings, but uppercase as well:
SYMBOL_DEFS="${SYMBOL_DEFS} ${SYMBOL_DEFS^^}"

FFLAGS="${FFLAGS} -O2 -fPIC -ffixed-line-length-none -cpp"
LIBOPENBLAS=openblas
if [[ ${nbits} == 64 ]]; then
    LIBOPENBLAS=openblas64_
    FFLAGS="${FFLAGS} -fdefault-integer-8 ${SYMBOL_DEFS}"
fi

cmake ../arpack-ng -DCMAKE_INSTALL_PREFIX="$prefix" -DCMAKE_TOOLCHAIN_FILE="/opt/$target/$target.toolchain" -DBUILD_SHARED_LIBS=ON -DBLAS_LIBRARIES="-L$prefix/lib -l${LIBOPENBLAS}" -DLAPACK_LIBRARIES="-L$prefix/lib -l${LIBOPENBLAS}" -DCMAKE_Fortran_FLAGS="${FFLAGS}" -DCMAKE_EXE_LINKER_FLAGS="${EXE_LINK_FLAGS}"

make -j${nproc} VERBOSE=1
make install VERBOSE=1

# For now, we'll have to adjust the name of the OpenBLAS library on macOS.
# Eventually, this should be fixed upstream
if [[ ${target} == "x86_64-apple-darwin14" ]]; then
    echo "-- Modifying library name for OpenBLAS"
    install_name_tool -change libopenblas64_.0.3.0.dev.dylib libopenblas64_.dylib ${prefix}/lib/libarpack.2.1.0.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line.  We enable the full
# combinatorial explosion of GCC versions because this package most
# definitely links against libgfortran.
platforms = expand_gcc_versions(supported_platforms())

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libarpack", :libarpack)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    "https://github.com/JuliaLinearAlgebra/OpenBLASBuilder/releases/download/v0.3.0-3/build_OpenBLAS.v0.3.0.jl"
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: daviehh/ArpackBuilder

max_issues_repo_name: daviehh/ArpackBuilder

max_forks_repo_name: daviehh/ArpackBuilder