text:
struct GJS_BWG_HNE_DL{K,H} <: UnknownAcvTimeSeriesModel{3,Float64}
    α::Float64
    ωₚ::Float64
    γ::Float64
    r::Float64
    s::Float64
    ϕₘ::Float64
    β::Float64
    ν::Float64
    σₗ::Float64
    σᵣ::Float64

    r_over_s::Float64
    ωₚ²::Float64
    ωₚ³::Float64
    ωₚˢ::Float64
    ωₚˢ⁻¹::Float64
    logγ::Float64
    s⁻¹::Float64
    logωₚminuss⁻¹::Float64
    σᵣ_over3::Float64
    cosϕₘ::Float64
    sinϕₘ::Float64
    cos2ϕₘ::Float64
    sin2ϕₘ::Float64
    function GJS_BWG_HNE_DL{K,H}(α,ωₚ,γ,r,s,ϕₘ,β,ν,σₗ,σᵣ) where {K,H}
        α > 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires α > 0"))
        ωₚ > 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires ωₚ > 0"))
        γ >= 1 || throw(ArgumentError("GJS_BWG_HNE_DL requires γ > 1"))
        r > 1 || throw(ArgumentError("GJS_BWG_HNE_DL requires r > 1"))
        s > 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires s > 0"))
        β >= 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires β > 0"))
        ν >= 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires ν > 0"))
        σₗ >= 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires σₗ > 0"))
        σᵣ >= 0 || throw(ArgumentError("GJS_BWG_HNE_DL requires σᵣ > 0"))
        new{K,H}(α,ωₚ,γ,r,s,ϕₘ,β,ν,σₗ,σᵣ,
        r/s,ωₚ^2,ωₚ^3,ωₚ^s,ωₚ^(s-1),log(γ),inv(s),log(ωₚ)-inv(s),
        σᵣ/3,cos(ϕₘ),sin(ϕₘ),cos(2ϕₘ),sin(2ϕₘ)
        )
    end
    function GJS_BWG_HNE_DL{K,H}(x::AbstractVector{Float64}) where {K,H}
        @boundscheck checkparameterlength(x,GJS_BWG_HNE_DL{K,H})
        @inbounds GJS_BWG_HNE_DL{K,H}(x[1], x[2], x[3], x[4], x[5], x[6], x[7], x[8], x[9], x[10])
    end
end

# functions to throw informative error if type parameter not provided
GJS_BWG_HNE_DL(x::AbstractVector{Float64}) = GJS_BWG_HNE_DL(ones(10)...)
GJS_BWG_HNE_DL(α,ωₚ,γ,r,s,ϕₘ,β,ν,σₗ,σᵣ) = error("GJS_BWG_HNE_DL process requires the ammount of aliasing and water depth to be specified as a type parameter. Use GJS_BWG_HNE_DL{K,H}() where K ∈ N₀ and H ∈ R.")

WhittleLikelihoodInference.npars(::Type{GJS_BWG_HNE_DL{K,H}}) where {K,H} = 10
WhittleLikelihoodInference.nalias(::GJS_BWG_HNE_DL{K,H}) where {K,H} = K

WhittleLikelihoodInference.lowerbounds(::Type{GJS_BWG_HNE_DL{K,H}}) where {K,H} = [0,0,1,1,0,-Inf,0,0,0,0]
WhittleLikelihoodInference.upperbounds(::Type{GJS_BWG_HNE_DL{K,H}}) where {K,H} = [Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]

@propagate_inbounds @fastmath function WhittleLikelihoodInference.add_sdf!(out, model::GJS_BWG_HNE_DL{K,H}, ω::Real) where {K,H}
    @boundscheck checkbounds(out,1:6)
    @inbounds begin
        s_om = sign(ω)
        ω = abs(ω)

        if ω > 1e-10
            tanhkh = tanh(approx_dispersion(ω, H))
            inv_tanhkh = inv(tanhkh)
            inv_tanhkh² = inv_tanhkh^2
            α,ωₚ,γ,r,s,β,ν,σₗ = model.α,model.ωₚ,model.γ,model.r,model.s,model.β,model.ν,model.σₗ

            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            ω_over_ωₚ = ω / ωₚ
            ωₚ_over_ω = 1/ω_over_ωₚ

            δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
            ω⁻ˢ = ω^(-s)
            ωₚˢ_over_ωˢ = model.ωₚˢ * ω⁻ˢ
            ∂S∂α = ω^(-r) * exp(-(model.r_over_s) * ωₚˢ_over_ωˢ) * γ^δ / 2
            sdf = α * ∂S∂α

            PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ_over_ω)^8)
            σ² = σ^2

            cosPS = cos(PeakSep)
            exp2part = (2 * exp(2.0 * σ²))
            cos2part = (model.cos2ϕₘ*cosPS) / exp2part
            imexppart = sdf * 1im * s_om / exp(0.5 * σ²) * cos(PeakSep/2)

            out[1] += sdf
            out[2] += imexppart * model.cosϕₘ * inv_tanhkh
            out[3] += imexppart * model.sinϕₘ * inv_tanhkh
            out[4] += sdf * (0.5 + cos2part) * inv_tanhkh²
            out[5] += sdf * model.sin2ϕₘ * cosPS / exp2part * inv_tanhkh²
            out[6] += sdf * (0.5 - cos2part) * inv_tanhkh²
        end # add zero otherwise
    end
    return nothing
end

@propagate_inbounds @fastmath function WhittleLikelihoodInference.grad_add_sdf!(out, model::GJS_BWG_HNE_DL{K,H}, ω::Real) where {K,H}
    @boundscheck checkbounds(out,1:6,1:10)
    @inbounds begin
        s_om = sign(ω)
        ω = abs(ω)

        if ω > 1e-10
            tanhkh = tanh(approx_dispersion(ω, H))
            inv_tanhkh = inv(tanhkh)
            inv_tanhkh² = inv_tanhkh^2
            α,ωₚ,γ,r,s,β,ν,σₗ,σᵣ = model.α,model.ωₚ,model.γ,model.r,model.s,model.β,model.ν,model.σₗ,model.σᵣ

            σ1² = 0.0049 + 0.0032 * (ω > ωₚ)
            ω_over_ωₚ = ω / ωₚ
            ωₚ_over_ω = 1/ω_over_ωₚ

            δ = exp(-1 / (2σ1²) * (ω / ωₚ - 1)^2)
            ω⁻ˢ = ω^(-s)
            ωₚˢ_over_ωˢ = model.ωₚˢ * ω⁻ˢ
            ∂S∂α = ω^(-r) * exp(-(model.r_over_s) * ωₚˢ_over_ωˢ) * γ^δ / 2
            sdf = α * ∂S∂α

            ∂S∂α = sdf/α
            ∂S∂ωₚ = sdf * (δ*model.logγ * ω / σ1² *(ω-ωₚ) / model.ωₚ³ - r*model.ωₚˢ⁻¹*ω⁻ˢ)
            ∂S∂γ = sdf * δ / γ
            ∂S∂r = sdf * (-log(ω)-ωₚˢ_over_ωˢ*model.s⁻¹)
            ∂S∂s = sdf * model.r_over_s * ωₚˢ_over_ωˢ * (log(ω)-model.logωₚminuss⁻¹)

            PeakSep = β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            σ = σₗ - model.σᵣ_over3*(4*(ωₚ_over_ω)^2 - (ωₚ_over_ω)^8)
            σ² = σ^2

            cosϕₘ = model.cosϕₘ
            sinϕₘ = model.sinϕₘ
            cos2ϕₘ = model.cos2ϕₘ
            sin2ϕₘ = model.sin2ϕₘ
            cosPShalf = cos(PeakSep/2)
            cosPS = cos(PeakSep)
            sinPShalf = sin(PeakSep/2)
            sinPS = sin(PeakSep)

            exp2part = (exp(2.0 * σ²))
            cos2part = (cos2ϕₘ*cosPS) / exp2part
            exphalf = exp(0.5 * σ²)
            imexppart = 1im * s_om / exphalf * cosPShalf

            # hzz = 1
            hxx = (0.5 + cos2part/2)
            hyy = (0.5 - cos2part/2)
            hxz = imexppart * cosϕₘ
            hyz = imexppart * sinϕₘ
            hyx = sin2ϕₘ * cosPS / exp2part / 2

            #### h_derivs
            ## ϕₘ
            hyxpart∂ϕₘ = cosPS / exp2part
            hxx∂ϕₘ = -sin2ϕₘ * hyxpart∂ϕₘ
            hyy∂ϕₘ = -hxx∂ϕₘ
            hzpart∂ϕₘ = -1im * s_om * cosPShalf / exphalf
            hxz∂ϕₘ = sinϕₘ * hzpart∂ϕₘ
            hyz∂ϕₘ = -cosϕₘ * hzpart∂ϕₘ
            hyx∂ϕₘ = cos2ϕₘ * hyxpart∂ϕₘ

            ## β, ν
            PeakSep∂β = PeakSep / β
            PeakSep∂ν = -((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0) * β / exp(ν * ((ω > ωₚ)*((ωₚ_over_ω) - 1.0) + 1.0))
            hx∂βνtemp = sinPS / 2 / exp2part
            hz∂βνtemp = -0.5im * sinPShalf / exphalf * s_om

            hx∂βt = hx∂βνtemp * PeakSep∂β
            hx∂νt = hx∂βνtemp * PeakSep∂ν
            hz∂βt = hz∂βνtemp * PeakSep∂β
            hz∂νt = hz∂βνtemp * PeakSep∂ν

            hxx∂β = -cos2ϕₘ * hx∂βt
            hyy∂β = -hxx∂β
            hxz∂β = cosϕₘ * hz∂βt
            hyz∂β = sinϕₘ * hz∂βt
            hyx∂β = -sin2ϕₘ * hx∂βt

            hxx∂ν = -cos2ϕₘ * hx∂νt
            hyy∂ν = -hxx∂ν
            hxz∂ν = cosϕₘ * hz∂νt
            hyz∂ν = sinϕₘ * hz∂νt
            hyx∂ν = -sin2ϕₘ * hx∂νt

            ## σₗ, σᵣ
            # ∂σ∂σₗ = 1
            ∂σ∂σᵣ = (σ-σₗ)/σᵣ
            hx∂σtemp = cosPS * 2 / exp2part * σ
            hz∂σtemp = -1im * cosPShalf / exphalf * σ * s_om

            hx∂σᵣt = hx∂σtemp * ∂σ∂σᵣ
            hz∂σᵣt = hz∂σtemp * ∂σ∂σᵣ

            hxx∂σₗ = -cos2ϕₘ * hx∂σtemp
            hyy∂σₗ = -hxx∂σₗ
            hxz∂σₗ = cosϕₘ * hz∂σtemp
            hyz∂σₗ = sinϕₘ * hz∂σtemp
            hyx∂σₗ = -sin2ϕₘ * hx∂σtemp

            hxx∂σᵣ = -cos2ϕₘ * hx∂σᵣt
            hyy∂σᵣ = -hxx∂σᵣ
            hxz∂σᵣ = cosϕₘ * hz∂σᵣt
            hyz∂σᵣ = sinϕₘ * hz∂σᵣt
            hyx∂σᵣ = -sin2ϕₘ * hx∂σᵣt

            ## ωₚ
            PeakSep∂ωₚ = -(ω > ωₚ) * ν/ω * β * exp(-ν*ωₚ/ω)
            ∂σ∂ωₚ = (-8σᵣ/3*(ωₚ/(ω^2)-ωₚ^7/(ω^8)))

            hx∂ωₚt = (sinPS / 2 * PeakSep∂ωₚ + 2*cosPS*σ*∂σ∂ωₚ)/ exp2part
            hz∂ωₚt = -1im * (sinPShalf / 2 * PeakSep∂ωₚ + cosPShalf*σ*∂σ∂ωₚ)/ exphalf * s_om

            hxx∂ωₚ = -cos2ϕₘ * hx∂ωₚt
            hyy∂ωₚ = -hxx∂ωₚ
            hxz∂ωₚ = cosϕₘ * hz∂ωₚt
            hyz∂ωₚ = sinϕₘ * hz∂ωₚt
            hyx∂ωₚ = -sin2ϕₘ * hx∂ωₚt

            ## α
            out[1, 1] += ∂S∂α
            out[2, 1] += ∂S∂α * hxz * inv_tanhkh
            out[3, 1] += ∂S∂α * hyz * inv_tanhkh
            out[4, 1] += ∂S∂α * hxx * inv_tanhkh²
            out[5, 1] += ∂S∂α * hyx * inv_tanhkh²
            out[6, 1] += ∂S∂α * hyy * inv_tanhkh²

            ## ωₚ
            out[1, 2] += ∂S∂ωₚ
            out[2, 2] += (∂S∂ωₚ * hxz + sdf * hxz∂ωₚ) * inv_tanhkh
            out[3, 2] += (∂S∂ωₚ * hyz + sdf * hyz∂ωₚ) * inv_tanhkh
            out[4, 2] += (∂S∂ωₚ * hxx + sdf * hxx∂ωₚ) * inv_tanhkh²
            out[5, 2] += (∂S∂ωₚ * hyx + sdf * hyx∂ωₚ) * inv_tanhkh²
            out[6, 2] += (∂S∂ωₚ * hyy + sdf * hyy∂ωₚ) * inv_tanhkh²

            ## γ
            out[1, 3] += ∂S∂γ
            out[2, 3] += ∂S∂γ * hxz * inv_tanhkh
            out[3, 3] += ∂S∂γ * hyz * inv_tanhkh
            out[4, 3] += ∂S∂γ * hxx * inv_tanhkh²
            out[5, 3] += ∂S∂γ * hyx * inv_tanhkh²
            out[6, 3] += ∂S∂γ * hyy * inv_tanhkh²

            ## r
            out[1, 4] += ∂S∂r
            out[2, 4] += ∂S∂r * hxz * inv_tanhkh
            out[3, 4] += ∂S∂r * hyz * inv_tanhkh
            out[4, 4] += ∂S∂r * hxx * inv_tanhkh²
            out[5, 4] += ∂S∂r * hyx * inv_tanhkh²
            out[6, 4] += ∂S∂r * hyy * inv_tanhkh²

            ## s
            out[1, 5] += ∂S∂s
            out[2, 5] += ∂S∂s * hxz * inv_tanhkh
            out[3, 5] += ∂S∂s * hyz * inv_tanhkh
            out[4, 5] += ∂S∂s * hxx * inv_tanhkh²
            out[5, 5] += ∂S∂s * hyx * inv_tanhkh²
            out[6, 5] += ∂S∂s * hyy * inv_tanhkh²

            ## ϕₘ
            # out[1, 5] += 0
            out[2, 6] += sdf * hxz∂ϕₘ * inv_tanhkh
            out[3, 6] += sdf * hyz∂ϕₘ * inv_tanhkh
            out[4, 6] += sdf * hxx∂ϕₘ * inv_tanhkh²
            out[5, 6] += sdf * hyx∂ϕₘ * inv_tanhkh²
            out[6, 6] += sdf * hyy∂ϕₘ * inv_tanhkh²

            ## β
            # out[1, 7] += 0
            out[2, 7] += sdf * hxz∂β * inv_tanhkh
            out[3, 7] += sdf * hyz∂β  * inv_tanhkh
            out[4, 7] += sdf * hxx∂β  * inv_tanhkh²
            out[5, 7] += sdf * hyx∂β * inv_tanhkh²
            out[6, 7] += sdf * hyy∂β * inv_tanhkh²

            ## ν
            # out[1, 8] += 0
            out[2, 8] += sdf * hxz∂ν * inv_tanhkh
            out[3, 8] += sdf * hyz∂ν  * inv_tanhkh
            out[4, 8] += sdf * hxx∂ν  * inv_tanhkh²
            out[5, 8] += sdf * hyx∂ν * inv_tanhkh²
            out[6, 8] += sdf * hyy∂ν * inv_tanhkh²

            ## σₗ
            # out[1, 9] += 0
            out[2, 9] += sdf * hxz∂σₗ * inv_tanhkh
            out[3, 9] += sdf * hyz∂σₗ * inv_tanhkh
            out[4, 9] += sdf * hxx∂σₗ * inv_tanhkh²
            out[5, 9] += sdf * hyx∂σₗ * inv_tanhkh²
            out[6, 9] += sdf * hyy∂σₗ * inv_tanhkh²

            ## σᵣ
            # out[1, 10] += 0
            out[2, 10] += sdf * hxz∂σᵣ * inv_tanhkh
            out[3, 10] += sdf * hyz∂σᵣ * inv_tanhkh
            out[4, 10] += sdf * hxx∂σᵣ * inv_tanhkh²
            out[5, 10] += sdf * hyx∂σᵣ * inv_tanhkh²
            out[6, 10] += sdf * hyy∂σᵣ * inv_tanhkh²

        end # add zero otherwise
    end
    return nothing
end


text:
# Container for log-density values
mutable struct PVals{FT <: AbstractFloat}
    x::FT
    pr_x::FT
end

# Initialise random walk adaptation
function init_rwm_adapt(algorithm, x0)
    if algorithm==:ram
        s = RobustAdaptiveMetropolis(x0)
    elseif algorithm==:am
        s = AdaptiveMetropolis(x0)
    elseif algorithm==:asm
        s = AdaptiveScalingMetropolis(x0)
    elseif algorithm==:aswam
        s = AdaptiveScalingWithinAdaptiveMetropolis(x0)
    else
        error("No such algorithm: $algorithm")
    end
    s
end

# Initialise temperature adaptation based on swap rates
function init_sw_adapt(L::Int, alpha_opt::FT=0.234, sc=FT(1.0),
    gamma=FT(0.66)) where {FT <: AbstractFloat}
    step = PolynomialStepSize(gamma, FT(L-1))
    Rhos = Vector{AdaptiveScalingMetropolis}(undef, L-1)
    for i in 1:L-1
        Rhos[i] = AdaptiveScalingMetropolis(alpha_opt, sc, step)
    end
    (Rhos...,)
end

# Update inverse temperatures based on adapted parameters
@inline function rho2beta!(Betas::Vector{FT}, Rhos) :: Nothing where {
    FT <: AbstractFloat}
    T = one(FT)
    Betas[1] = T
    @inbounds for i = 1:length(Rhos)
        T += exp(Rhos[i].sc[])
        Betas[i+1] = one(FT)/T
    end
    nothing
end

# Single random walk update & possible adaptation
@inline function arwm_step!(r::RWMState, s::AdaptState, p::PVals{FT},
    beta::FT, log_p::Function, log_pr::Function, k::Int, do_adapt::Bool) :: Int where {
    FT <: AbstractFloat }
    # Draw proposal, evaluate log-target, acc.rate
    draw!(r, s)
    p_y = log_p(r.y); p_pr_y = log_pr(r.y)
    dp = beta*(p_y - p.x) + p_pr_y - p.pr_x
    #alpha = min(one(FT), exp(dp))
    rxy = exp(dp); alpha = (rxy<one(FT)) ? rxy : one(FT) # This makes r=NaN -> alpha=1
    if rand(r.rng) <= alpha
        # Accept: interchange current/proposed states & log-target values
        accept!(r)
        p.x = p_y; p.pr_x = p_pr_y
        acc = 1
    else
        acc = 0
    end
    if do_adapt
        adapt!(s, r, alpha, k)
    end
    acc
end

# Attempt swap move between levels lev_ <-> lev & possible adaptation
@inline function swap_step!(R, Rhos, P, Betas::Vector{FT}, lev_::Int, lev::Int,
    k::Real, do_adapt::Bool) :: Int where {FT <: AbstractFloat}
    p_ = P[lev_]; beta_ = Betas[lev_]
    p  = P[lev];  beta = Betas[lev]

    # Acceptance probability:
    dp = (beta_ - beta) * (p.x - p_.x)
    #alpha = min(one(FT), exp(dp))
    rxy = exp(dp); alpha = (rxy<one(FT)) ? rxy : one(FT) # This makes r=NaN -> alpha=1

    # Accept swap:
    if rand(R[lev_].rng) <= alpha
        R[lev], R[lev_] = R[lev_], R[lev]
        P[lev], P[lev_] = p_, p
        acc = 1
    else
        acc = 0
    end

    # Do swap adaptation:
    if do_adapt
        adapt!(Rhos[lev_], R[lev_], alpha, k)
        rho2beta!(Betas, Rhos)
    end
    acc
end

# Initialise the key variables in the adaptive random walk Metropolis algoritm
@inline function init_arwm(x0::T, algorithm, rng, q, L, all_levels, nX,
                           log_p, log_pr) where {FT <: AbstractFloat,
                           T <: AbstractVector{FT}}
    d = length(x0)
    # Initialise output variables
    X = Vector(undef, L)
    D = Vector(undef, L)
    # RWM states, adaptation states...
    R = Vector{RWMState{d,FT,T,typeof(q),typeof(rng)}}(undef, L)
    S = Vector(undef, L)
    P = Vector{PVals{FT}}(undef, L)
    for lev = 1:L
        r = RWMState(x0, rng, q); R[lev] = r
        P[lev] = PVals(log_p(r.x), log_pr(r.x))
        # Initialise adaptation
        if typeof(algorithm) == Symbol
            S[lev] = init_rwm_adapt(algorithm, x0)
        else
            S[lev] = algorithm[lev]
        end
        if lev == 1 || all_levels
            X[lev] = zeros(FT, d, nX)
            D[lev] = zeros(FT, nX)
        else
            X[lev] = missing
            D[lev] = missing
        end
    end
    (X...,), (D...,), R, (S...,), P
end

"""
    out = adaptive_rwm(x0, log_p, n; kwargs)

Generic adaptive random walk Metropolis algorithm from initial state vector
`x0` targetting log probability density `log_p` run for `n` iterations,
including adaptive parallel tempering.

# Arguments
- `x0::Vector{<:AbstractFloat}`: The initial state vector
- `log_p::Function`: Function that returns log probability density values
                     (up to an additive constant) for any state vector.
- `n::Int`: Total number of iterations

# Keyword arguments
- `algorithm::Symbol`: The random walk adaptation algorithm; current choices are
  `:ram` (default), `:am`, `:asm`, `:aswam` and `:rwm`.
  (Alternatively, if algorithm is a vector of AdaptState, then this will be used as an
  initial state for adaptation.)
- `b::Int`: Number of burn-in iterations (before starting to collect samples); default `n/5`
- `thin::Int`: Thinning factor; only every `thin`:th sample is stored; default `1`
- `fulladapt::Bool`: Whether to adapt after burn-in; default `true`
- `rng::AbstractRNG`: Random number generator; default `Random.GLOBAL_RNG`
- `q::Function`: Zero-mean symmetric proposal generator (with arguments `x` and `rng`);
   default `q=randn!(x, rng)`
- `L::Int`: Number of parallel tempering levels
- `acc_sw::AbstractFloat`: Desired acceptance rate between level swaps; default `0.234`
- `all_levels::Bool`: Whether to store output of all levels; default `false`
- `log_pr::Function`: Log-prior density function; default `log_pr(x) = 0.0`.
- `swaps::Symbol`: Swap strategy, one of:
   `:single` (default, single randomly picked swap)
   `:randperm` (swap in random order)
   `:sweep` (up- or downward sweep, picked at random)
   `:nonrev` (alternate even/odd sites as in Syed, Bouchard-Côté, Deligiannidis,
   Doucet, 	arXiv:1905.02939)

Note that if `log_pr` is supplied, then `log_p(x)` is regarded as the
log-likelihood (or, equivalently, log-target is `log_p(x) + log_pr(x)`).
Tempering is only applied to `log_p`, not to `log_pr`.

The output `out.X contains the simulated samples (column vectors).
`out.allX[k]` for `k>=2` contain higher temperature auxiliary chains (if requested)

# Examples
```
log_p(x) = -.5*sum(x.^2)
o = adaptive_rwm(zeros(2), log_p, 10_000; algorithm=:am)
using MCMCChains, StatsPlots # Assuming MCMCChains & StatsPlots are installed...
c = Chains(o.X[1]', start=o.params.b, thin=o.params.thin); plot(c)
```
"""
function adaptive_rwm(x0::T, log_p::Function, n::Int;
    algorithm::Union{Symbol,Vector{<:AdaptState}}=:ram,
    thin::Int=1, b::Int=Int(floor(n/5)), fulladapt::Bool=true,
    q::Function=randn!, L::Int=1, log_pr::Function = (x->zero(FT)),
    all_levels::Bool=false, acc_sw::FT = FT(0.234), swaps::Symbol = :single,
    rng::AbstractRNG=Random.GLOBAL_RNG) where {FT <: AbstractFloat,
    T <: AbstractVector{FT}}

    args = (x0, log_p, n)
    params = (algorithm=algorithm, thin=thin, b=b, fulladapt=fulladapt,
              q=q, L=L, log_pr=log_pr, all_levels=all_levels, acc_sw=acc_sw,
              swaps=swaps, rng=copy(rng))

    # Initialise key variables
    nX = length(b:thin:n) # Number of output
    X, D, R, S, P = init_arwm(x0, algorithm, rng, q, L,
                              all_levels, nX, log_p, log_pr)
    # Acceptance statistics
    accRWM = zeros(L)  # How many are accepted
    accSW = zeros(L-1); nSW = zeros(L-1)

    # Shorthands for all levels, swap levels, stored levels
    allLevels = 1:L
    oddLevels = 1:2:L-1; evenLevels = 2:2:L-1
    upSweep = 1:(L-1); downSweep = L-1:-1:1
    if swaps == :single
        sc_sw = FT(L-1)
    elseif swaps == :randperm
        sc_sw = one(FT)
        sweepLevels = Vector{Int}(undef, L-1)
    elseif swaps == :nonrev
        sc_sw = one(FT)*2
    else # :sweep
        sc_sw = one(FT)
    end
    Rhos = init_sw_adapt(L, acc_sw, sc_sw)

    # Which ones are stored
    storeLevels = all_levels ? (1:L) : [1]

    # Initialise inverse temperatures
    Betas = zeros(FT, L)
    rho2beta!(Betas, Rhos)

    for k = 1:n
        # Random walk moves:
        for lev = allLevels
            accRWM[lev] += arwm_step!(R[lev], S[lev], P[lev], Betas[lev],
                                 log_p, log_pr, k, fulladapt || k <= b)
        end

        # Swap move between random adjacent levels
        if L>1
            if swaps == :single
                lev = rand(rng, upSweep);
                nSW[lev] += 1
                accSW[lev] += swap_step!(R, Rhos, P, Betas, lev, lev+1,
                                      k/L, fulladapt || k <= b)
            else
                if swaps == :nonrev
                    if k % 2 == 1
                        sweepLevels = oddLevels
                    else
                        sweepLevels = evenLevels
                    end
                elseif swaps == :randperm
                    randperm!(rng, sweepLevels)
                else # :sweep
                    if rand(rng) < 0.5
                        sweepLevels = upSweep
                    else
                        sweepLevels = downSweep
                    end
                end
                for lev in sweepLevels
                    nSW[lev] += 1
                    accSW[lev] += swap_step!(R, Rhos, P, Betas, lev, lev+1,
                                          k, fulladapt || k <= b)
                end
            end
        end

        # Store output
        for lev in storeLevels
            if k>=b && rem(k-b, thin) == 0
                i = Int((k-b)/thin)+1
                X[lev][:,i] .= R[lev].x
                D[lev][i] = Betas[lev]*P[lev].x + P[lev].pr_x
            end
        end
    end

    (X=X[1], allX=X, D=D, R=R, S=S, Rhos=Rhos, accRWM=accRWM/n, accSW=accSW./nSW,
    args=args, params=params)
end

ext: jl


text:
import .JuMP: Model, @objective, @variable, @constraint, termination_status, dual, value

"""
    ot_jump(D, P1, P2)

Solve the optimal transport problem using JuMP. This function is only available if `using JuMP, GLPK`.

# Arguments:
- `D`: Distance matrix
- `P1`: Weight vector 1
- `P2`: Weight vector 2
"""
function ot_jump(D, P1, P2; kwargs...)
    n = length(P1)
    @assert all(isfinite, P1) "Got nonfinite P1"
    @assert all(isfinite, P2) "Got nonfinite P2"
    @assert sum(P1) ≈ 1 "Sum(P1) ≠ 1: $(sum(P1))"
    @assert sum(P2) ≈ 1 "Sum(P2) ≠ 1: $(sum(P2))"

    model = Model(GLPK.Optimizer)
    @variable(model, γ[1:n^2])
    @objective(model, Min, γ'vec(D))
    Γ = reshape(γ,n,n)
    con1 = JuMP.@constraint(model, con1,  sum(Γ, dims=1)[:] .== P2)
    con2 = @constraint(model, con2,  sum(Γ, dims=2)[:] .== P1)
    con3 = @constraint(model, con3,  γ .>= 0)
    JuMP.optimize!(model)
    if Int(termination_status(model)) != 1
        @error Int(termination_status(model))
    end
    α, β = -dual.(con2)[:], -dual.(con1)[:]
    α .-= mean(α)
    β .-= mean(β)
    reshape(value.(γ),n,n), α, β #, -dual.(con3)[:]
end


text:
struct PortRange
    "Initial port in the range or the value of the port for single ports"
    start::Int
    "End port in the range"
    stop::Int
    "Is this a range (`true`) or a single port (`false`)"
    r::Bool
end
PortRange(p::Integer) = PortRange(p, -1, false)
PortRange(p::UnitRange) = PortRange(Int(p.start), Int(p.stop), true)

"""
`PortRange(p::AbstractString)`
`PortRange(p::Integer)`
`PortRange(p::UnitRange)`

Creates a `PortRange` object that handles continuous ranges of ports in multiple scanners.

In the DTC Initium, up to 8 scanners can be used simultaneously.
Each scanner can hav 8, 16, 32 or 64 pressure ports.
Since the total number ports can be large, when specifying which ports
should be used ranges are a necessity.

The DTC Initium uses the following convention to name pressure ports:

 * XYZ
 * X is the number of the slot where the scanner is connected (1-8)
 * YZ is the number of the port where 1 ≤ YZ ≤ N where N is the total number of ports in the scanner
 * A range is given by ABC-XYZ
 * If X is different from A, then all ports between ABC and XYZ are included

As an example, consider that 8 scanners with 64 ports are used.
In this case the port range 228-547 consists of the following pressure ports

 * 228-264: ports 28-64 in the scanner located in slot 2
 * 301-364: ports 28-64 in the scanner located in slot 3
 * 401-464: ports 28-64 in the scanner located in slot 4
 * 501-547: ports 28-64 in the scanner located in slot 5

It should be noted that *just a `PortRange` object does not specify the exact ports that
are used*. There is need for context: the exact scanners that are connected to the
DTC Initum.

The function [`strport`](@ref) returns a string that represents the port range and
can be used in DTC Initium commands.

## Examples
```jldoctest
julia> PortRange("101-353") # 101-353
PortRange(101, 353, true)

julia> PortRange(101:353) # 101-353
PortRange(101, 353, true)

julia> PortRange(412) # 412
PortRange(412, -1, false)
```

"""
function PortRange(p::AbstractString)
    p = strip(p)

    r1 = r"^[0-9][0-9][0-9]$"
    r2 = r"^[0-9][0-9][0-9]-[0-9][0-9][0-9]$"

    if occursin(r1, p)
        return PortRange(parse(Int, p), -1, false)
    elseif occursin(r2, p)
        i = findfirst(isequal('-'), p)
        p1 = parse(Int, p[1:(i-1)])
        p2 = parse(Int, p[(i+1):end])
        return PortRange(p1, p2, true)
    else
        throw(ArgumentError(p, "Not a valid port or port range"))
    end

end

"Is the [`PortRange`](@ref) object a range (`true`) or a single port (`false`)"
isrange(p::PortRange) = p.r

"""
`strport(p::PortRange)`

Creates the corresponding string of the [`PortRange`](@ref) object.

## Example
```jldoctest
julia> strport(PortRange(101:353)) # 101-353
"101-353"
```
"""
function strport(p::PortRange)

    if isrange(p)
        return "$(p.start)-$(p.stop)"
    else
        return "$(p.start)"
    end
end

"""
`portlist(ports::AbstractString)`
`portlist(ports...)`

Creates a vector of [`PortRange`](@ref) objects.

The input can be a string with ranges and single ports separated by
spaces (`' '`) or commas (`','`).

A single string or multiple arguments can be passed to this function.
In the second case, the arguments should be compatible with the
constructors available to [`PortRange`](@ref).

## Examples
```jldoctest
julia> portlist("401-464", 505:543)
2-element Vector{PortRange}:
 PortRange(401, 464, true)
 PortRange(505, 543, true)

julia> portlist("401-464", 505:543, 555)
3-element Vector{PortRange}:
 PortRange(401, 464, true)
 PortRange(505, 543, true)
 PortRange(555, -1, false)

julia> portlist("101-234 245-264 ,, 301-325, 401-464")
4-element Vector{PortRange}:
 PortRange(101, 234, true)
 PortRange(245, 264, true)
 PortRange(301, 325, true)
 PortRange(401, 464, true)
```
"""
function portlist(ports::AbstractString)

    ports = replace(strip(ports), "," => " ")

    plst = PortRange[]
    c = ' '
    if occursin(c, ports)
        for p in split(ports, c, keepempty=false)
            push!(plst, PortRange(strip(p)))
        end
    else
        push!(plst, PortRange(strip(ports)))
    end

    return plst
end

function portlist(ports...)
    plst = PortRange[]

    for p in ports
        push!(plst, PortRange(p))
    end
    return plst
end

"""
`strportlist(ports::AbstractVector{PortRange})`

Translates a vector of [`PortRange`](@ref) to a string that can be used
as an argument to the DTC Initium command SD3.

This function is called by function [`SD3`](@ref) to add pressure channels
to the [`Initium`](@ref) device.

## Example
```jldoctest
julia> plst = portlist("101-234 245-264 ,, 301-325, 401-464")
4-element Vector{PortRange}:
 PortRange(101, 234, true)
 PortRange(245, 264, true)
 PortRange(301, 325, true)
 PortRange(401, 464, true)

julia> import DTCInitium: strportlist

julia> strportlist(plst)
" 101-234 245-264 301-325 401-464"
```
"""
function strportlist(ports::AbstractVector{PortRange})
    plst = ""

    for p1 in ports
        if isrange(p1)
            i1 = p1.start
            i2 = p1.stop
            plst *= " $i1-$i2"
        else
            plst *= " $(p1.start)"
        end
    end
    return plst
end

"""
`checkportlist(scanners, ports::AbstractVector{PortRange})`

Check if a set of port ranges (argument `ports`) are compatible with the
scanners configured.

This will depend on which scanners are connected to the DTC Initium and
the number of pressure ports on each scanner.

## Example
```jldoctest
julia> scanners = [(1,64,1), (2,64,1), (3,64,1), (4,64,1)]; # 4 scanners with 64 channels.

julia> checkportlist(scanners, portlist("101-464"))
true

julia> checkportlist(scanners, portlist("101-564"))
false

```
"""
function checkportlist(scanners, ports::AbstractVector{PortRange})

    scn = [s[1] for s in scanners]
    npp = [s[2] for s in scanners]

    ii = sortperm(scn)
    scn = scn[ii]
    npp = npp[ii]

    for p in ports
        p1 = p.start
        if !isrange(p)
            s = floor(Int, (p1-1)/100)
            if s ∉ scn
                return false
            else
                idxs = findfirst(isequal(s), scn)
                nn = npp[idxs]
                if p1 > s*100 + nn
                    return false
                end
            end
        else
            p2 = p.stop
            s1 = floor(Int, (p1-1)/100)
            s2 = floor(Int, (p2-1)/100)
            for s in s1:s2
                if s ∉ scn
                    return false
                end
            end
            idxs = findfirst(isequal(s2), scn)
            nn = npp[idxs]
            if p2 > s2*100+nn
                return false
            end

        end
    end
    return true

end

"""
`defscanlist(scanners, ports::AbstractVector{PortRange})`

Return the list of pressure ports in a vector of [`PortRange`](@ref) objects.

## Example
```jldoctest
julia> scanners = [(1,64,1), (2,64,1), (3,64,1), (4,64,1)];

julia> defscanlist(scanners, portlist("101-104 131-134 315-320"))
14-element Vector{Int64}:
 101
 102
 103
 104
 131
 132
 133
 134
 315
 316
 317
 318
 319
 320
```
"""
function defscanlist(scanners, ports::AbstractVector{PortRange})
    scn = [s[1] for s in scanners]
    npp = [s[2] for s in scanners]

    ii = sortperm(scn)
    scn = scn[ii]
    npp = npp[ii]

    plst = Int[]
    for p in ports
        if !isrange(p)
            pn = p.start
            s = floor(Int, pn/100)  # Scanner number
            if s ∉ scn  # Not valid scanner
                throw(DomainError(pn, "Port $pn not part of any scanner"))
            end
            idx = pn - 100*s
            if idx > npp[findfirst(isequal(s), scn)]
                throw(BoundsError(pn, "Port $pn not valid for scanner $s"))
            end
            push!(plst, pn)
        else
            pstart = p.start
            pstop = p.stop
            s1 = floor(Int, pstart/100)
            s2 = floor(Int, pstop/100)
            for s in s1:s2
                if s ∉ scn
                    throw(DomainError(s, "Ports $p have scanners not configured!"))
                end
            end

            i = findfirst(isequal(s1), scn)
            if !(1 ≤ (pstart-100*s1) ≤ npp[i])
                throw(BoundsError(pstart, "Illegal port number!"))
            end

            i = findfirst(isequal(s2), scn)
            if !(1 ≤ (pstop-100*s2) ≤ npp[i])
                throw(BoundsError(pstart, "Illegal port number!"))
            end

            if s2 != s1
                append!(plst, pstart:(s1*100+npp[i]))
                for s in (s1+1):(s2-1)
                    i = findfirst(isequal(s), scn)
                    append!(plst, (s*100) .+ (1:npp[i]))
                end
                i = findfirst(isequal(s2), scn)
                append!(plst, (s2*100+1):pstop)
            else
                i = findfirst(isequal(s2), scn)
                append!(plst, pstart:pstop)
            end
        end
    end

    return plst

end


text:
module MomentWrapper

using JavaCall

const AbstractStorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.AbstractStorelessUnivariateStatistic
const AbstractUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.AbstractUnivariateStatistic
const AggregatableStatistic = @jimport org.hipparchus.stat.descriptive.AggregatableStatistic
const Class = @jimport java.lang.Class
const DoubleConsumer = @jimport java.util.function.DoubleConsumer
const FirstMoment = @jimport org.hipparchus.stat.descriptive.moment.FirstMoment
const FourthMoment = @jimport org.hipparchus.stat.descriptive.moment.FourthMoment
const GeometricMean = @jimport org.hipparchus.stat.descriptive.moment.GeometricMean
const Iterable = @jimport java.lang.Iterable
const JString = @jimport java.lang.String
const Kurtosis = @jimport org.hipparchus.stat.descriptive.moment.Kurtosis
const Mean = @jimport org.hipparchus.stat.descriptive.moment.Mean
const Object = @jimport java.lang.Object
const SecondMoment = @jimport org.hipparchus.stat.descriptive.moment.SecondMoment
const SemiVariance = @jimport org.hipparchus.stat.descriptive.moment.SemiVariance
const SemiVariance_Direction = @jimport org.hipparchus.stat.descriptive.moment.SemiVariance$Direction
const Skewness = @jimport org.hipparchus.stat.descriptive.moment.Skewness
const StandardDeviation = @jimport org.hipparchus.stat.descriptive.moment.StandardDeviation
const StorelessUnivariateStatistic = @jimport org.hipparchus.stat.descriptive.StorelessUnivariateStatistic
const SumOfLogs = @jimport org.hipparchus.stat.descriptive.summary.SumOfLogs
const ThirdMoment = @jimport org.hipparchus.stat.descriptive.moment.ThirdMoment
const UnivariateStatistic = @jimport org.hipparchus.stat.descriptive.UnivariateStatistic
const Variance = @jimport org.hipparchus.stat.descriptive.moment.Variance
const WeightedEvaluation = @jimport org.hipparchus.stat.descriptive.WeightedEvaluation

export FirstMoment
export FourthMoment
export GeometricMean
export Kurtosis
export Mean
export SecondMoment
export SemiVariance
export Skewness
export StandardDeviation
export ThirdMoment
export Variance
export accept
export aggregate
export and_then
export clear
export copy
export equals
export evaluate
export get_class
export get_data
export get_n
export get_result
export get_variance_direction
export hash_code
export increment
export increment_all
export is_bias_corrected
export notify
export notify_all
export set_data
export to_string
export wait
export with_bias_corrected
export with_bias_correction
export with_variance_direction

include("first_moment.jl")
include("fourth_moment.jl")
include("geometric_mean.jl")
include("kurtosis.jl")
include("mean.jl")
include("second_moment.jl")
include("semi_variance.jl")
include("skewness.jl")
include("standard_deviation.jl")
include("third_moment.jl")
include("variance.jl")

end


text:
using IterativeSolvers

A = Float64[
  16 4 8;
  4 5 -4;
  8 -4 22;
]
a = Float64[16; 18; -22]
x = Float64[1; 1; 1]

maxiters = 100
println("\nJacobi:")
for i in 1:100
  res = jacobi!(x, A, a, tol=1.0e-8, maxiter=1)
  if i < 5 || i in [1; 10:10:maxiters]
    println([i x'])
  end
  if res[2].isconverged
    println([i x'])
    break
  end
end

println("\nGauss-Seidel:")
@show  A * x
println()

x = Float64[1; 1; 1]
for i in 1:100
  res = gauss_seidel!(x, A, a, tol=1.0e-8, maxiter=1)
  if i < 5 || i in [1; 10:10:maxiters]
    println([i x'])
  end
  if res[2].isconverged
    println([i x'])
    break
  end
end

println("\nSequential Succesive Over-Relaxation:")
@show A * x
println()

x = Float64[1; 1; 1]
ω = 1.0
for i in 1:100
  res = ssor!(x, A, a, ω, tol=1.0e-8, maxiter=1)
  if i < 5 || i in [1; 10:10:maxiters]
    println([i x'])
  end
  if res[2].isconverged
    println([i x'])
    break
  end
end

println()
A * x


text:
using MAGMA
using Test, Random, LinearAlgebra

# for simply tests the very beginning codes for gesvd
@testset "dense_linear_algebra" begin
    # Write your own tests here.
    include("dense/runtests.jl")
end


text:
reload("Persa")

using Base.Test
using DecisionTree
using DatasetsCF

# write your own tests here
#@test 1 == 2
###
reload("COFILS")

dataset = DatasetsCF.MovieLens()

holdout = Persa.HoldOut(dataset, 0.9)

(ds_train, ds_test) = Persa.get(holdout)

model = COFILS.Cofils(ds_train, 10)
Persa.train!(model, ds_train)

print(Persa.aval(model, ds_test))


text:
"""
    function is_constraint_violated(
        com::CoM,
        constraint::BoolConstraint,
        fct,
        set::AndSet,
    )

Check if one of the inner constraints is violated
"""
function is_constraint_violated(
    com::CoM,
    constraint::BoolConstraint,
    fct,
    set::AndSet,
)
    return is_lhs_constraint_violated(com, constraint) || is_rhs_constraint_violated(com, constraint)
end

"""
    still_feasible(com::CoM, constraint::AndConstraint, fct, set::AndSet, vidx::Int, value::Int)

Return whether the constraint can be still fulfilled when setting a variable with index `vidx` to `value`.
**Attention:** This assumes that it isn't violated before.
"""
function still_feasible(
    com::CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet,
    vidx::Int,
    value::Int,
)
    lhs_indices = constraint.lhs.indices
    for i in 1:length(lhs_indices)
        if lhs_indices[i] == vidx
            !still_feasible(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set, vidx, value) && return false
        end
    end
    rhs_indices = constraint.rhs.indices
    for i in 1:length(rhs_indices)
        if rhs_indices[i] == vidx
            !still_feasible(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set, vidx, value) && return false
        end
    end
    return true
end

"""
    prune_constraint!(com::CS.CoM, constraint::AndConstraint, fct, set::AndSet; logs = true)

Reduce the number of possibilities given the `AndConstraint` by pruning both parts
Return whether still feasible
"""
function prune_constraint!(
    com::CS.CoM,
    constraint::AndConstraint,
    fct,
    set::AndSet;
    logs = true,
)
    !activate_lhs!(com, constraint) && return false
    !prune_constraint!(com, constraint.lhs, constraint.lhs.fct, constraint.lhs.set; logs=logs) && return false
    !activate_rhs!(com, constraint) && return false
    feasible = prune_constraint!(com, constraint.rhs, constraint.rhs.fct, constraint.rhs.set; logs=logs)
    return feasible
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 05653593-4365-46a2-8e1b-79c3c8e11575
begin
	using Pkg
	cd(joinpath(dirname(@__FILE__),".."))
    Pkg.activate(pwd())
	Pkg.instantiate()
    using NativeSVG
	using Plots
	using Optim
end

# ╔═╡ 282c0f80-015b-11eb-19d8-d9235a2a77d4
md"# Global Search Algorithms"

# ╔═╡ 4d9273a0-015a-11eb-2caa-130d0245a946
md"""The solution of the unconstrained optimization problem

```math
\min f\left(\vec x\right)
```

can always be found by solving

```math
\nabla f\left(\vec x^\star\right) = \vec 0\,.
```

We choose an initial point and the generate a sequence of iterates. Typically, the best we can hope for is that the sequence converges to a local minimizer. For this reason, it is often desirable for the initial point to be close to a global optimizer. Moreover, these methods require first derivatives (and als second derivatives in the case of Newton's method).

In this lecture we discuss a method that is global in nature in the sense that they attempt to search throughout the entire feasible set.

The _Nelder-Mead simplex_ method uses only objective function values and do not require derivatives. Consequently, it is applicable to a much wider class of optimization problems."""

# ╔═╡ bfe6caf0-015a-11eb-23a0-919e0520bd8d
md"""## The Nelder-Mead Simplex Algorithm

The Nelder-Mead algorithm uses the concept of a _simplex_. A simplex is a geometric object determined by an assembly of $n+1$ points, $\vec{p}_{0},\vec{p}_{1},\dots,\vec{p}_{n}$, in the n-dimensional space such that

```math
\det\begin{bmatrix}
\vec{p}_{0} & \vec{p}_{1} & \cdots & \vec{p}_{n}\\
1 & 1 & \cdots & 1
\end{bmatrix}\neq0\,.
```

This condition ensures that two points in $\mathbb R$ do not coincide, three points in $\mathbb R^{2}$ are not collinear, four points in $\mathbb R^{3}$ are not coplanar, and so on. Thus, a simplex in $\mathbb R$ is a line segment, in $\mathbb R^{2}$ it is a triangle, while a simplex in $\mathbb R^{3}$ is a tetrahedron; in each case it encloses a finite $n$-dimensional volume.

Suppose that we wish to minimize $f\left(\vec{x}\right),\,\vec{x}\in\mathbb R^{n}$.
To start the algorithm, we initialize a simplex of $n+1$ points.
A possible way to set up a simplex is to start with an initial point
$\vec{x}^{\left(0\right)}=\vec{p}_{0}$ and generate the remaining
points of the initial simplex as follows:

```math
\vec{p}_{i}=\vec{p}_{0}+\lambda_{i}\vec{e}_{i},\quad i=1,2,\dots,n\:,
```

where the $\vec{e}_{i}$ are unit vectors constituting the natural
basis of $\mathbb R^{n}$. The positive constant coefficients $\lambda_{i}$
are selected in such a way that their magnitudes reflect the length
scale of the optimization problem.

Our objective is to modify the initial simplex stage by stage so that
the resulting simplices converge toward the minimizer. At each iteration
we evaluate the function $f$ at each point of the simplex. In the
function minimization process, the point with the largest function
value is replaced with another point. The process for modifying the
simplex continues until it converges toward the function minimizer.

We now present the rules for modifying the simplex stage by stage.
To aid in our presentation, we use a two-dimensional example to illustrate
the rules. We begin by selecting the initial set of $n+1$ points
that are to form the initial simplex. We next evaluate $f$ at each
point and order the $n+1$ vertices to satisfy

```math
f\left(\vec{p}_{0}\right)\leq f\left(\vec{p}_{1}\right)\leq\cdots\leq f\left(\vec{p}_{n}\right)\,.
```

For the two-dimensional case we let $\vec{p}_{l}$, $\vec{p}_{n}$
and $\vec{p}_{s}$ denote the points of the simplex for which $f$
is largest, next largest and smallest. We next compute $\vec{p}_{g}$,
the _centroid_ (center of gravity) of the best $n$ points:

```math
\vec{p}_{g}=\frac{1}{n}\sum_{i=0}^{n-1}\vec{p}_{i}\,.
```

We then reflect the worst vertex, $\vec{p}_{n}$, in $\vec{p}_{g}$
using a _reflection coefficient_ $\rho>0$ to obtain the reflection
point

```math
\vec{p}_{r}=\vec{p}_{g}+\rho\left(\vec{p}_{g}-\vec{p}_{n}\right)\,.
```"""

# ╔═╡ 0aa85db0-015b-11eb-01ae-694c21b1f7a4
Drawing(width = "12cm", height = "6cm", viewBox = "0 0 600 300", preserveAspectRatio="none") do
    circle(cx = "52", cy = "150", r = "5", fill = "black")
    circle(cx = "250", cy = "50", r = "5", fill = "black")
    circle(cx = "250", cy = "250", r = "5", fill = "black")
    circle(cx = "250", cy = "150", r = "5", fill = "black")
    circle(cx = "450", cy = "150", r = "5", fill = "black")
    polyline(
        points="50, 150 250, 150 250, 50 50, 150 250, 250 250, 150 450 150",
        fill="none", stroke="black", stroke_width="2.5"
    )
    NativeSVG.text(x = "40", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pl") end
    NativeSVG.text(x = "260", y = "40", width = "60", height = "60", font_size = "16pt") do
		str("pn") end
    NativeSVG.text(x = "260", y = "240", width = "60", height = "60", font_size = "16pt") do
		str("ps") end
    NativeSVG.text(x = "210", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pg") end
    NativeSVG.text(x = "440", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pr") end
end

# ╔═╡ c6569c40-016d-11eb-2d6b-4d8a54cff763
md"""A typical value is $\rho=1$.

We proceed to evaluate $f$ at $\vec{p}_{r}$. If $f\left(\vec{p}_{0}\right)\leq f\left(\vec{p}_{r}\right)<f\left(\vec{p}_{n-1}\right)$,
then the point $\vec{p}_{r}$ replaces $\vec{p}_{l}$ to form a new
simplex and we terminate the iteration. We proceed to repeat the process.

If , however, $f\left(\vec{p}_{r}\right)<f\left(\vec{p}_{0}\right)$,
so that the point $\vec{p}_{r}$ yields the smallest function value
among the points of the simplex, we argue that this direction is a
good one. In this case we increase the distance traveled using an
_expansion coefficient_ $\chi>1$, e.g. $\chi=2$, to obtain

```math
\vec{p}_{e}=\vec{p}_{g}+\chi\left(\vec{p}_{r}-\vec{p}_{g}\right)\,.
```
"""

# ╔═╡ 0e6ae94e-016e-11eb-3f19-676f24836ebb
Drawing(width = "16cm", height = "6cm", viewBox = "0 0 800 300", preserveAspectRatio="none") do
    circle(cx = "50", cy = "150", r = "5", fill = "black")
    circle(cx = "250", cy = "50", r = "5", fill = "black")
    circle(cx = "250", cy = "250", r = "5", fill = "black")
    circle(cx = "250", cy = "150", r = "5", fill = "black")
    circle(cx = "450", cy = "150", r = "5", fill = "black")
    circle(cx = "650", cy = "150", r = "5", fill = "black")
    polyline(
        points="50, 150 250, 150 250, 50 50, 150 250, 250 250, 150 650 150",
        fill="none", stroke="black", stroke_width="2.5"
    )
    NativeSVG.text(x = "40", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pl") end
    NativeSVG.text(x = "260", y = "40", width = "60", height = "60", font_size = "16pt") do
		str("pn") end
    NativeSVG.text(x = "260", y = "240", width = "60", height = "60", font_size = "16pt") do
		str("ps") end
    NativeSVG.text(x = "210", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pg") end
    NativeSVG.text(x = "440", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pr") end
    NativeSVG.text(x = "640", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pe") end
end

# ╔═╡ 6bf23d80-016e-11eb-1c64-6d2da5f0f954
md"""This operation yields a new point on the line $\vec{p}_{n}\vec{p}_{g}\vec{p}_{r}$
extend beyond $\vec{p}_{r}$. If $f\left(\vec{p}_{e}\right)<f\left(\vec{p}_{r}\right)$,
the expansion is declared a success and $\vec{p}_{e}$ replaces $\vec{p}_{n}$
in the next simplex. If, on the other hand, $f\left(\vec{p}_{e}\right)\geq f\left(\vec{p}_{r}\right)$,
the expansion is a failure and $\vec{p}_{r}$ replaces $\vec{p}_{n}$.

Finally, if $f\left(\vec{p}_{r}\right)\geq f\left(\vec{p}_{n-1}\right)$,
the reflected point $\vec{p}_{r}$ would constitute the point with
the largest function value in the new simplex. Then in the next step
it would be reflected in $\vec{p}_{g}$, probably an unfruitful operation.
Instead, this case is dealt with by a _contraction_ operation
in one of two ways. First, if $f\left(\vec{p}_{r}\right)<f\left(\vec{p}_{n}\right)$,
then we contract $\left(\vec{p}_{r}-\vec{p}_{g}\right)$ with a _contraction
coefficient_ $0<\gamma<1$, e.g. $\gamma=\frac{1}{2}$, to obtain"""

# ╔═╡ 7cf25a20-016e-11eb-30f1-8394316a43c2
Drawing(width = "16cm", height = "6cm", viewBox = "0 0 800 300", preserveAspectRatio="none") do
    circle(cx = "50", cy = "150", r = "5", fill = "black")
    circle(cx = "250", cy = "50", r = "5", fill = "black")
    circle(cx = "250", cy = "250", r = "5", fill = "black")
    circle(cx = "250", cy = "150", r = "5", fill = "black")
    circle(cx = "450", cy = "150", r = "5", fill = "black")
    circle(cx = "350", cy = "150", r = "5", fill = "black")
    polyline(
        points="50, 150 250, 150 250, 50 50, 150 250, 250 250, 150 450 150",
        fill="none", stroke="black", stroke_width="2.5"
    )
    NativeSVG.text(x = "40", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pl") end
    NativeSVG.text(x = "260", y = "40", width = "60", height = "60", font_size = "16pt") do
		str("pn") end
    NativeSVG.text(x = "260", y = "240", width = "60", height = "60", font_size = "16pt") do
		str("ps") end
    NativeSVG.text(x = "210", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pg") end
    NativeSVG.text(x = "440", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pr") end
    NativeSVG.text(x = "340", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pc") end
end

# ╔═╡ c5157440-016e-11eb-33fd-4b9c155cf832
md"""```math
\vec{p}_{c}=\vec{p}_{g}+\gamma\left(\vec{p}_{r}-\vec{p}_{g}\right)\,.
```

We refer to this operation as the _outside contraction_. If,
on the other hand, $f\left(\vec{p}_{r}\right)\geq f\left(\vec{p}_{n}\right)$,
then $\vec{p}_{l}$ replaces $\vec{p}_{r}$ in the contraction operation
and we get

```math
\vec{p}_{c}=\vec{p}_{g}+\gamma\left(\vec{p}_{l}-\vec{p}_{g}\right)\,.
```

This operation is referred to as the _inside contraction_. """

# ╔═╡ fce70140-016e-11eb-224a-3d863b6817b6
Drawing(width = "16cm", height = "6cm", viewBox = "0 0 800 300", preserveAspectRatio="none") do
    circle(cx = "50", cy = "150", r = "5", fill = "black")
    circle(cx = "250", cy = "50", r = "5", fill = "black")
    circle(cx = "250", cy = "250", r = "5", fill = "black")
    circle(cx = "250", cy = "150", r = "5", fill = "black")
    circle(cx = "150", cy = "150", r = "5", fill = "black")
    polyline(
        points="50, 150 250, 150 250, 50 50, 150 250, 250 250, 150",
        fill="none", stroke="black", stroke_width="2.5"
    )
    NativeSVG.text(x = "40", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pl") end
    NativeSVG.text(x = "260", y = "40", width = "60", height = "60", font_size = "16pt") do
		str("pn") end
    NativeSVG.text(x = "260", y = "240", width = "60", height = "60", font_size = "16pt") do
		str("ps") end
    NativeSVG.text(x = "210", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pg") end
    NativeSVG.text(x = "140", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pc") end
end

# ╔═╡ 18554920-0186-11eb-0cb2-f5d020aff0d9
md"""If , in either case, $f\left(\vec{p}_{c}\right)\leq f\left(\vec{p}_{n}\right)$,
the contraction is considered a success, and we replace $\vec{p}_{n}$
with $\vec{p}_{c}$ in the new simplex. If, however, $f\left(\vec{p}_{c}\right)>f\left(\vec{p}_{n}\right)$,
the contraction is a failure, and in this case a new simplex can be
formed by retaining $\vec{p}_{0}$ only and halving the distance from
$\vec{p}_{0}$ to every point in the simplex. We can refer to this
event as a _shrinkage operation_. In general, the shrink step
produces the $n$ new vertices of the new simplex according to the
formula

```math
\vec{v}_{i}=\vec{p}_{0}+\sigma\left(\vec{p}_{i}-\vec{p}_{0}\right)\,,\quad i=1,2,\dots,n\,,
```"""

# ╔═╡ 392084d2-0186-11eb-11e5-2b5128129e19
Drawing(width = "16cm", height = "6cm", viewBox = "0 0 800 300", preserveAspectRatio="none") do
    circle(cx = "50", cy = "150", r = "5", fill = "black")
    circle(cx = "250", cy = "50", r = "5", fill = "black")
    circle(cx = "250", cy = "250", r = "5", fill = "black")
    circle(cx = "250", cy = "150", r = "5", fill = "black")
    circle(cx = "150", cy = "200", r = "5", fill = "black")
    polyline(
        points="250, 150 250, 50 50, 150 250, 250 250, 150 150, 200",
        fill="none", stroke="black", stroke_width="2.5"
    )
	NativeSVG.text(x = "40", y = "170", width = "60", height = "60", font_size = "16pt") do
		str("pl") end
    NativeSVG.text(x = "260", y = "40", width = "60", height = "60", font_size = "16pt") do
		str("pn") end
    NativeSVG.text(x = "260", y = "240", width = "60", height = "60", font_size = "16pt") do
		str("ps") end
end

# ╔═╡ 6e95360e-0186-11eb-1ed3-770285bc0e74
md"""where $\sigma=\frac{1}{2}$. Hence, the vertices of the new simplex
are $\vec{p}_{0},\vec{v}_{1},\dots,\vec{v}_{n}$.

When implementing the simplex algorithm, we need a tie-breaking rule
to order points in the case of equal function values. Often, the highest
possible index consistent with the relation

```math
f\left(\vec{p}_{0}\right)\leq f\left(\vec{p}_{1}\right)\leq\cdots\leq f\left(\vec{p}_{n}\right)
```

is assigned to the new vertex."""

# ╔═╡ 86bb011e-0186-11eb-2b14-8ddd7f9ca6c2
md"""## Example

To show how the `Optim` package can be used, we minimize the Rosenbrock function, a classical test problem for numerical optimization. We'll assume that you've already installed the Optim package using Julia's package manager. First, we load Optim and define the Rosenbrock function:"""

# ╔═╡ 202572f0-0187-11eb-0172-9bea075d89a9
f(x) = (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2;

# ╔═╡ d332f672-0186-11eb-205a-97aeef39b75d
x, y, z = let
	x = -1:0.007:2
	y = -1:0.007:2
	z = Surface((x,y)->log10.(f([x,y])), x, y)
	x, y, z
end;

# ╔═╡ 316b4b20-0187-11eb-33a2-dda15c1b3ca9
contour(x,y,z)

# ╔═╡ 4680a1e0-0187-11eb-3a9c-5b631ea7c822
md"""Once we've defined this function, we can find the minimizer (the input that minimizes the objective) and the minimum (the value of the objective at the minimizer) using any of our favorite optimization algorithms. With a function defined, we just specify an initial point x and call optimize with a starting point `x0`:"""

# ╔═╡ 5815d290-0187-11eb-391d-ff073eb9dc85
x0 = [0.0, 0.0];

# ╔═╡ 75d97390-0187-11eb-1b3f-69c879780178
optimize(f, x0, NelderMead())

# ╔═╡ 7927b430-0187-11eb-345f-a10499db73fa
md"""Other local solvers are available. Below, we use BFGS, a quasi-Newton method that requires a gradient. If we pass f alone, Optim will construct an approximate gradient for us using central finite differencing:"""

# ╔═╡ 86f5f1d0-0187-11eb-1186-ad0df4a1ef6c
optimize(f, x0, BFGS())

# ╔═╡ 8dd8f1f0-0187-11eb-370a-a1b5053e91e4
md"""For better performance and greater precision, you can pass your own gradient function. If your objective is written in all Julia code with no special calls to external (that is non-Julia) libraries, you can also use automatic differentiation, by using the `autodiff` keyword and setting it to `:forward`:"""

# ╔═╡ 97810580-0187-11eb-38a3-fd83c085e89c
optimize(f, x0, BFGS(); autodiff = :forward)

# ╔═╡ 9fb64f30-0187-11eb-1a26-7f9e577886d9
md"""For the Rosenbrock example, the analytical gradient can be shown to be:"""

# ╔═╡ ad0c2880-0187-11eb-31c1-bdfd311b2650
function g!(G, x)
	G[1] = -2.0 * (1.0 - x[1]) - 400.0 * (x[2] - x[1]^2) * x[1]
	G[2] = 200.0 * (x[2] - x[1]^2)
end

# ╔═╡ ba3db0f0-0187-11eb-2efa-296cff31defa
md"""Note, that the functions we're using to calculate the gradient (and later the Hessian `h!`) of the Rosenbrock function mutate a fixed-sized storage array, which is passed as an additional argument called `G` (or `H` for the Hessian) in these examples. By mutating a single array over many iterations, this style of function definition removes the sometimes considerable costs associated with allocating a new array during each call to the `g!` or `h!` functions. """

# ╔═╡ c818e6e0-0187-11eb-2c4a-b78c74a10231
optimize(f, g!, x0, BFGS())

# ╔═╡ cd8db100-0187-11eb-0ac7-2f4510dc2f08
md"""In addition to providing gradients, you can provide a Hessian function `h!` as well. In our current case this is"""

# ╔═╡ d80597b0-0187-11eb-23ef-7daa3ac30f36
function h!(H, x)
    H[1, 1] = 2.0 - 400.0 * x[2] + 1200.0 * x[1]^2
    H[1, 2] = -400.0 * x[1]
    H[2, 1] = -400.0 * x[1]
    H[2, 2] = 200.0
end

# ╔═╡ e07a5b60-0187-11eb-246d-f18cbf91b1f1
md"""Now we can use Newton's method for optimization by running:"""

# ╔═╡ f2674900-0187-11eb-017d-e5d26651604d
optimize(f, g!, h!, x0)

# ╔═╡ f933b3e0-0187-11eb-113f-876918888691
md"""Which defaults to `Newton()` since a Hessian function was provided. Like gradients, the Hessian function will be ignored if you use a method that does not require it.

Note that `Optim` will not generate approximate Hessians using finite differencing because of the potentially low accuracy of approximations to the Hessians. Other than Newton's method, none of the algorithms provided by the Optim package employ exact Hessians."""

# ╔═╡ 2837d5e2-0188-11eb-1eb3-67c8aeb366d0
md"""## Algorithms

The first choice to be made is that of the order of the method. Zeroth-order methods do not have gradient information, and are very slow to converge, especially in high dimension. First-order methods do not have access to curvature information and can take a large number of iterations to converge for badly conditioned problems. Second-order methods can converge very quickly once in the vicinity of a minimizer. Of course, this enhanced performance comes at a cost: the objective function has to be differentiable, you have to supply gradients and Hessians, and, for second order methods, a linear system has to be solved at each step.

If you can provide analytic gradients and Hessians, and the dimension of the problem is not too large, then second order methods are very efficient. The _Newton method with trust region_ is the method of choice.

When you do not have an explicit Hessian or when the dimension becomes large enough that the linear solve in the Newton method becomes the bottleneck, first order methods should be preferred. _BFGS_ is a very efficient method, but also requires a linear system solve. _L-BFGS_ usually has a performance very close to that of BFGS, and avoids linear system solves (the parameter `m` can be tweaked: increasing it can improve the convergence, at the expense of memory and time spent in linear algebra operations). The _conjugate gradient method_ usually converges less quickly than L-BFGS, but requires less memory. _Gradient descent_ should only be used for testing.

When the objective function is non-differentiable or you do not want to use gradients, use zeroth-order methods. _Nelder-Mead_ is currently the most robust."""

# ╔═╡ 35b66ce0-0188-11eb-2b93-fb6f4e8fe9ec
md"""## Linesearches

Linesearches are used in every first- and second-order method except for the trust-region Newton method. Linesearch routines attempt to locate quickly an approximate minimizer of the univariate function

```math
\alpha \rightarrow f\left(\vec x + \alpha \vec d\right)\,,
```

where $\vec d$ is the descent direction computed by the algorithm. They vary in how accurate this minimization is. Two good linesearches are BackTracking and HagerZhang, the former being less stringent than the latter. For well-conditioned objective functions and methods where the step is usually well-scaled (such as L-BFGS or Newton), a rough linesearch such as BackTracking is usually the most performant. For badly behaved problems or when extreme accuracy is needed (gradients below the square root of the machine epsilon, about `10−8` with Float64), the HagerZhang method proves more robust. An exception is the conjugate gradient method which requires an accurate linesearch to be efficient, and should be used with the HagerZhang linesearch."""

# ╔═╡ 648f6030-0188-11eb-2ca2-7b2e080de91b
md"""## Summary

As a very crude heuristic:

For a low-dimensional problem with analytic gradients and Hessians, use the Newton method with trust region. For larger problems or when there is no analytic Hessian, use L-BFGS, and tweak the parameter `m` if needed. If the function is non-differentiable, use Nelder-Mead.

Use the `HagerZhang` linesearch for robustness and `BackTracking` for speed.

Requires only a function handle:

- `NelderMead()`
- `SimulatedAnnealing()`, see next year
- `ParticleSwarm()`, see next year

Requires a function and gradient (will be approximated if omitted):

- `BFGS()`
- `LBFGS()`
- `ConjugateGradient()`
- `GradientDescent()`

Requires a function, a gradient, and a Hessian (cannot be omitted):

- `Newton()`
- `NewtonTrustRegion()`

Special methods for bounded univariate optimization:

- `Brent()`
- `GoldenSection()`"""

# ╔═╡ Cell order:
# ╟─282c0f80-015b-11eb-19d8-d9235a2a77d4
# ╟─05653593-4365-46a2-8e1b-79c3c8e11575
# ╟─4d9273a0-015a-11eb-2caa-130d0245a946
# ╟─bfe6caf0-015a-11eb-23a0-919e0520bd8d
# ╟─0aa85db0-015b-11eb-01ae-694c21b1f7a4
# ╟─c6569c40-016d-11eb-2d6b-4d8a54cff763
# ╟─0e6ae94e-016e-11eb-3f19-676f24836ebb
# ╟─6bf23d80-016e-11eb-1c64-6d2da5f0f954
# ╟─7cf25a20-016e-11eb-30f1-8394316a43c2
# ╟─c5157440-016e-11eb-33fd-4b9c155cf832
# ╟─fce70140-016e-11eb-224a-3d863b6817b6
# ╟─18554920-0186-11eb-0cb2-f5d020aff0d9
# ╟─392084d2-0186-11eb-11e5-2b5128129e19
# ╟─6e95360e-0186-11eb-1ed3-770285bc0e74
# ╟─86bb011e-0186-11eb-2b14-8ddd7f9ca6c2
# ╠═202572f0-0187-11eb-0172-9bea075d89a9
# ╠═d332f672-0186-11eb-205a-97aeef39b75d
# ╠═316b4b20-0187-11eb-33a2-dda15c1b3ca9
# ╟─4680a1e0-0187-11eb-3a9c-5b631ea7c822
# ╠═5815d290-0187-11eb-391d-ff073eb9dc85
# ╠═75d97390-0187-11eb-1b3f-69c879780178
# ╟─7927b430-0187-11eb-345f-a10499db73fa
# ╠═86f5f1d0-0187-11eb-1186-ad0df4a1ef6c
# ╟─8dd8f1f0-0187-11eb-370a-a1b5053e91e4
# ╠═97810580-0187-11eb-38a3-fd83c085e89c
# ╟─9fb64f30-0187-11eb-1a26-7f9e577886d9
# ╠═ad0c2880-0187-11eb-31c1-bdfd311b2650
# ╟─ba3db0f0-0187-11eb-2efa-296cff31defa
# ╠═c818e6e0-0187-11eb-2c4a-b78c74a10231
# ╟─cd8db100-0187-11eb-0ac7-2f4510dc2f08
# ╠═d80597b0-0187-11eb-23ef-7daa3ac30f36
# ╟─e07a5b60-0187-11eb-246d-f18cbf91b1f1
# ╠═f2674900-0187-11eb-017d-e5d26651604d
# ╟─f933b3e0-0187-11eb-113f-876918888691
# ╟─2837d5e2-0188-11eb-1eb3-67c8aeb366d0
# ╟─35b66ce0-0188-11eb-2b93-fb6f4e8fe9ec
# ╟─648f6030-0188-11eb-2ca2-7b2e080de91b


text:
module Exp

using Base.Test
using Bitmanip, Defs, FloatProps, Rounding

export vexp, vexp10, vexp2

# Note: The Coefficients were determined by a Mathematica script,
# available at <https://bitbucket.org/eschnett/vecmathlib>. They
# minimize the error in a least-square sense over the interval dx ∈
# [-0.5; +0.5]. A better alternative would be to minimize the maximum
# of the error.

# One can use a smaller set of coefficients to increase speed at the
# extent of a slightly larger error.

# TODO: Use @evalpoly for this
@inline function vexp2_poly(dx::Float32)
    # Coefficients for Float32,
    # error=4.55549108005200277750378992345e-9
    r = mkFloat(Float32, 0.000154653240842602623787395880898)
    r = muladd(r, dx, mkFloat(Float32, 0.00133952915439234389712105060319))
    r = muladd(r, dx, mkFloat(Float32, 0.0096180399118156827664944870552))
    r = muladd(r, dx, mkFloat(Float32, 0.055503406540531310853149866446))
    r = muladd(r, dx, mkFloat(Float32, 0.240226511015459465468737123346))
    r = muladd(r, dx, mkFloat(Float32, 0.69314720007380208630542805293))
    r = muladd(r, dx, mkFloat(Float32, 0.99999999997182023878745628977))
    return r
end

@inline function vexp2_poly(dx::Float64)
    # Coefficients for Float64,
    # error=9.32016781355638010975628074746e-18
    r = mkFloat(Float64, 4.45623165388261696886670014471e-10)
    r = muladd(r, dx, mkFloat(Float64, 7.0733589360775271430968224806e-9))
    r = muladd(r, dx, mkFloat(Float64, 1.01780540270960163558119510246e-7))
    r = muladd(r, dx, mkFloat(Float64, 1.3215437348041505269462510712e-6))
    r = muladd(r, dx, mkFloat(Float64, 0.000015252733849766201174247690629))
    r = muladd(r, dx, mkFloat(Float64, 0.000154035304541242555115696403795))
    r = muladd(r, dx, mkFloat(Float64, 0.00133335581463968601407096905671))
    r = muladd(r, dx, mkFloat(Float64, 0.0096181291075949686712855561931))
    r = muladd(r, dx, mkFloat(Float64, 0.055504108664821672870565883052))
    r = muladd(r, dx, mkFloat(Float64, 0.240226506959101382690753994082))
    r = muladd(r, dx, mkFloat(Float64, 0.69314718055994530864272481773))
    r = muladd(r, dx, mkFloat(Float64, 0.9999999999999999978508676375))
    return r
end

@inline function vexp2{T}(x::T)
    # Here it does not matter how ties are broken
    xi = vround(x)
    dx = x - xi
    @vassert mkFloat(T,-0.5) <= dx <= mkFloat(T,0.5)
    r = vexp2_poly(dx)
    # r = vldexp(r, mkInt(T, xi))
    r = vldexp(r, xi)
    r = vifthen(x < mkFloat(T, min_exponent(T)), mkFloat(T,0.0), r)
    r = vifthen(x > mkFloat(T, max_exponent(T)), mkFloat(T,Inf), r)
    return r
end

const M_LOG2E = 1.442695040888963407359924681001892137426645954152985934135449406931109219181187e+00
@inline function vexp{T}(x::T)
    # TODO: Re-calculate the coefficients for this function to
    # increase accuracy and avoid the multiplication
    return vexp2(mkFloat(T,M_LOG2E)*x)
end

const M_LOG10E = 4.342944819032518276511289189166050822943970058036665661144537831658646492088688e-01
@inline function vexp10{T}(x::T)
    # TODO: Re-calculate the coefficients for this function to
    # increase accuracy and avoid the multiplication
    return vexp2(mkFloat(T,M_LOG10E)*x)
end

function selftest_type(T::Type)
    zero = convert(T,0)
    min = realmin(T)
    eps = Base.eps(T)
    max = realmax(T)
    vals = T[zero, nextfloat(zero),
             prevfloat(min), min, nextfloat(min),
             0.1*eps, eps, 10.0*eps,
             0.1*sqrt(eps), sqrt(eps), 10.0*sqrt(eps),
             1.0-eps, 1.0, 1.0+eps,
             2.0, 10.0, 100.0,
             0.1*sqrt(max), sqrt(max), 10.0*sqrt(max),
             0.1*max, prevfloat(max), max, Inf, NaN]
    for uval in vals
        for val in (uval, -uval)
            @test isequal(vexp2(val), exp2(val))
        end
    end
end

function selftest()
    selftest_type(Float32)
    selftest_type(Float64)
end

selftest()

end


text:
using Powseeker
using ForwardDiff
using StaticArrays

n = 1000
rng = MersenneTwister(1)
r = [rand(rng, SVector{2,Float64}) for i in 1:n]

for i in 1:n
    if norm(grad_peaks(r[i]) - ForwardDiff.gradient(peaks, r[i])) > 0.0000001
        println(grad_peaks(r[i]))
        println(ForwardDiff.gradient(peaks, r[i]))
    end
end

@time for i in 1:n
    grad_peaks(r[i])
end

@time for i in 1:n
    ForwardDiff.gradient(peaks, r[i])
end

for i in 1:n
    if norm(grad_scaled_peaks(r[i]) - ForwardDiff.gradient(scaled_peaks, r[i])) > 0.0000001
        println(grad_scaled_peaks(r[i]))
        println(ForwardDiff.gradient(scaled_peaks, r[i]))
    end
end

for i in 1:n
    if norm(grad_peaks2(r[i]) - ForwardDiff.gradient(peaks2, r[i])) > 0.0000001
        @show grad_peaks2(r[i])
        @show ForwardDiff.gradient(peaks2, r[i])
    end
end

@time for i in 1:n
    grad_peaks2(r[i])
end

@time for i in 1:n
    ForwardDiff.gradient(peaks2, r[i])
end

for i in 1:n
    if norm(grad_scaled_peaks2(r[i]) - ForwardDiff.gradient(scaled_peaks2, r[i])) > 0.0000001
        @show grad_scaled_peaks2(r[i])
        @show ForwardDiff.gradient(scaled_peaks2, r[i])
    end
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: workmailmessageflow
using AWS.Compat
using AWS.UUIDs

"""
    get_raw_message_content(message_id)
    get_raw_message_content(message_id, params::Dict{String,<:Any})

Retrieves the raw content of an in-transit email message, in MIME format.

# Arguments
- `message_id`: The identifier of the email message to retrieve.

"""
get_raw_message_content(messageId; aws_config::AbstractAWSConfig=global_aws_config()) = workmailmessageflow("GET", "/messages/$(messageId)"; aws_config=aws_config)
get_raw_message_content(messageId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = workmailmessageflow("GET", "/messages/$(messageId)", params; aws_config=aws_config)

"""
    put_raw_message_content(content, message_id)
    put_raw_message_content(content, message_id, params::Dict{String,<:Any})

Updates the raw content of an in-transit email message, in MIME format. This example
describes how to update in-transit email message. For more information and examples for
using this API, see  Updating message content with AWS Lambda.  Updates to an in-transit
message only appear when you call PutRawMessageContent from an AWS Lambda function
configured with a synchronous  Run Lambda rule. If you call PutRawMessageContent on a
delivered or sent message, the message remains unchanged, even though GetRawMessageContent
returns an updated message.

# Arguments
- `content`: Describes the raw message content of the updated email message.
- `message_id`: The identifier of the email message being updated.

"""
put_raw_message_content(content, messageId; aws_config::AbstractAWSConfig=global_aws_config()) = workmailmessageflow("POST", "/messages/$(messageId)", Dict{String, Any}("content"=>content); aws_config=aws_config)
put_raw_message_content(content, messageId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = workmailmessageflow("POST", "/messages/$(messageId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("content"=>content), params)); aws_config=aws_config)


text:
module FlowNodeCouplingFlowsPlanarFlowTest

using Test
using ReactiveMP
using ReactiveMP: getdim, getu, getb, getw, getall, setu!, setb!, setw!
using ReactiveMP: forward, forward!, jacobian, jacobian!, inv_jacobian, inv_jacobian!
using ReactiveMP: det_jacobian, absdet_jacobian, logdet_jacobian, logdet_jacobian, logabsdet_jacobian

@testset "Planar Flow" begin

    @testset "Constructors" begin

        # check for unspecified input
        f = PlanarFlow()
        @test typeof(f) == ReactiveMP.PlanarFlowPlaceholder

        # check when creating placeholder
        f = ReactiveMP.PlanarFlowPlaceholder()
        @test typeof(f) == ReactiveMP.PlanarFlowPlaceholder

        # check when creating empty planar flow
        f = ReactiveMP.PlanarFlowEmpty(3)
        @test typeof(f) == ReactiveMP.PlanarFlowEmpty{3}
        @test getdim(f) == 3

        # check for specified dimensionality
        f = PlanarFlow(1)
        @test typeof(f.u) == Float64
        @test typeof(f.w) == Float64
        @test typeof(f.b) == Float64

        f = PlanarFlow(5)
        @test typeof(f.u) == Array{Float64,1}
        @test typeof(f.w) == Array{Float64,1}
        @test typeof(f.b) == Float64
        @test length(f.u) == 5
        @test length(f.w) == 5

        # check for prespecified parameters (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        @test typeof(f.u) == Float64
        @test typeof(f.w) == Float64
        @test typeof(f.b) == Float64
        @test f.u == 1.0
        @test f.w == 2.0
        @test f.b == 3.0

        # check for prespecified parameters (multivariates)
        f = PlanarFlow([2.0, 3.0], [5.0, 6.0], 3.0)
        @test typeof(f.u) == Array{Float64,1}
        @test typeof(f.w) == Array{Float64,1}
        @test typeof(f.b) == Float64
        @test f.u == [2.0, 3.0]
        @test f.w == [5.0, 6.0]
        @test f.b == 3.0

        # wrong constructors
        @test_throws MethodError PlanarFlow([2.0, 3.0], 2.0, 1.0)
        @test_throws MethodError PlanarFlow([5.0, 6.0], 1.0, 0.5)
        @test_throws MethodError PlanarFlow(2.0, [2.0, 3.0], 1.0)
        @test_throws MethodError PlanarFlow(1.0, [5.0, 6.0], 0.5)
        @test_throws AssertionError PlanarFlow([2.0, 3.0, 4.0], [2.0, 3.0], 1.0)
        @test_throws AssertionError PlanarFlow([1.0, 6.0, 8.0], [5.0, 9.0], -2.0)
        @test_throws AssertionError PlanarFlow([2.0, 3.0], [2.0, 3.0, 4.0], 1.0)
        @test_throws AssertionError PlanarFlow([5.0, 9.0], [1.0, 6.0, 8.0], -2.0)

    end

    @testset "Prepare-Compile" begin

        f = ReactiveMP.PlanarFlowPlaceholder()
        @test ReactiveMP.prepare(2, f) == ReactiveMP.PlanarFlowEmpty(2)
        @test ReactiveMP.prepare(5, f) == ReactiveMP.PlanarFlowEmpty(5)
        @test ReactiveMP.prepare(9, f) == ReactiveMP.PlanarFlowEmpty(9)

        f = ReactiveMP.PlanarFlowEmpty(1)
        params = [1.0, 2.0, 3.0]
        fc = ReactiveMP.compile(f)
        fcp = ReactiveMP.compile(f, params)
        @test typeof(fc) <: PlanarFlow
        @test typeof(fc.u) == Float64
        @test typeof(fc.w) == Float64
        @test typeof(fc.b) == Float64
        @test typeof(fcp) <: PlanarFlow
        @test typeof(fcp.u) == Float64
        @test typeof(fcp.w) == Float64
        @test typeof(fcp.b) == Float64
        @test fcp.u == 1.0
        @test fcp.w == 2.0
        @test fcp.b == 3.0

        f = ReactiveMP.PlanarFlowEmpty(2)
        params = [1.0, 2.0, 3.0, 4.0, 5.0]
        fc = ReactiveMP.compile(f)
        fcp = ReactiveMP.compile(f, params)
        @test typeof(fc) <: PlanarFlow
        @test typeof(fc.u) == Array{Float64,1}
        @test typeof(fc.w) == Array{Float64,1}
        @test typeof(fc.b) == Float64
        @test typeof(fcp) <: PlanarFlow
        @test typeof(fcp.u) == Array{Float64,1}
        @test typeof(fcp.w) == Array{Float64,1}
        @test typeof(fcp.b) == Float64
        @test fcp.u == [1.0, 2.0]
        @test fcp.w == [3.0, 4.0]
        @test fcp.b == 5.0

    end

    @testset "nr_params" begin

        for k = 1:10
            f = ReactiveMP.PlanarFlowEmpty(k)
            fc = compile(f)
            @test nr_params(f)  == 2*k + 1
            @test nr_params(fc) == 2*k + 1
        end

    end

    @testset "Get-Set" begin

        # check get functions for univariate PlanarFlow
        f = PlanarFlow(1.0, 2.0, 3.0)
        @test getu(f) == f.u
        @test getw(f) == f.w
        @test getb(f) == f.b
        @test getall(f) == (f.u, f.w, f.b)
        @test getu(f) == 1.0
        @test getw(f) == 2.0
        @test getb(f) == 3.0
        @test getall(f) == (1.0, 2.0, 3.0)

        # check get functions for multivariate PlanarFlow
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        @test getu(f) == f.u
        @test getw(f) == f.w
        @test getb(f) == f.b
        @test getall(f) == (f.u, f.w, f.b)
        @test getu(f) == [1.0, 2.0]
        @test getw(f) == [3.0, 4.0]
        @test getb(f) == 5.0
        @test getall(f) == ([1.0, 2.0], [3.0, 4.0], 5.0)

        # check setu! function (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        setu!(f, 4.0)
        @test getu(f) == 4.0
        @test getw(f) == 2.0
        @test getb(f) == 3.0

        # check setu! function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        setu!(f, [10.0, 11.0])
        @test getu(f) == [10.0, 11.0]
        @test getw(f) == [3.0, 4.0]
        @test getb(f) == 5.0

        # check setw! function (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        setw!(f, 4.0)
        @test getu(f) == 1.0
        @test getw(f) == 4.0
        @test getb(f) == 3.0

        # check setw! function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        setw!(f, [10.0, 11.0])
        @test getu(f) == [1.0, 2.0]
        @test getw(f) == [10.0, 11.0]
        @test getb(f) == 5.0

        # check setb! function (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        setb!(f, 4.0)
        @test getu(f) == 1.0
        @test getw(f) == 2.0
        @test getb(f) == 4.0

        # check setb! function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        setb!(f, 6.0)
        @test getu(f) == [1.0, 2.0]
        @test getw(f) == [3.0, 4.0]
        @test getb(f) == 6.0

        # check getdim
        f = ReactiveMP.PlanarFlowEmpty(2)
        @test ReactiveMP.getdim(f) == 2

        # check errors (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        @test_throws MethodError setu!(f, [1.0, 2.0])
        @test_throws MethodError setw!(f, [1.0, 2.0])
        @test_throws MethodError setb!(f, [1.0, 2.0])

        # check errors (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        @test_throws MethodError setu!(f, 3.0)
        @test_throws MethodError setw!(f, 4.0)
        @test_throws MethodError setb!(f, [1.0, 2.0])
        @test_throws AssertionError setu!(f, [1.0, 3.0, 6.0])
        @test_throws AssertionError setw!(f, [1.0, 3.0, 6.0])

    end

    @testset "Base" begin

        # check base functions (univariate)
        f = PlanarFlow(1.0, 2.0, 3.0)
        @test eltype(f) == Float64
        @test eltype(PlanarFlow{Float64,Float64}) == Float64
        @test size(f) == 1
        @test length(f) == 1

        # check base functions (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 5.0)
        @test eltype(f) == Float64
        @test eltype(PlanarFlow{Array{Float64,1},Float64}) == Float64
        @test size(f) == 2
        @test length(f) == 2

        # check base functions empty object
        f = ReactiveMP.PlanarFlowEmpty(2)
        @test size(f)   == 2
        @test length(f) == 2

    end

    @testset "Forward-Backward" begin

        # check forward function (univariate)
        f = PlanarFlow(1.0, 2.0, -3.0)
        @test forward(f, 1.5) == 1.5
        @test forward(f, 2.5) == 3.464027580075817
        @test forward.(f, [1.5, 2.5]) == [1.5, 3.464027580075817]

        # check forward function (univariate function, multivariate input)
        f = PlanarFlow(1.0, 2.0, -3.0)
        @test forward(f, [1.5]) == 1.5
        @test forward(f, [2.5]) == 3.464027580075817
        @test forward.(f, [[1.5], [2.5]]) == [1.5, 3.464027580075817]

        # check forward function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 0.0)
        @test forward(f, [-4.0, 3.0]) == [-4.0, 3.0]
        @test forward(f, [-2.0, 1.5]) == [-2.0, 1.5]
        @test forward.(f, [[-4.0, 3.0], [-2.0, 1.5]]) == [[-4.0, 3.0], [-2.0, 1.5]]

        # check forward! function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 0.0)
        output = zeros(2)
        forward!(output, f, [-4.0, 3.0])
        @test output == [-4.0, 3.0]
        forward!(output, f, [-2.0, 1.5])
        @test output == [-2.0, 1.5]

    end

    @testset "Jacobian" begin

        # check jacobian function (univariate)
        f = PlanarFlow(1.0, 2.0, -3.0)
        @test jacobian(f, 1.5) == 3.0
        @test jacobian(f, 2.5) == 1.1413016497063289
        @test jacobian.(f, [1.5, 2.5]) == [3.0, 1.1413016497063289]

        # check jacobian function (univariate function, multivariate input)
        f = PlanarFlow(1.0, 2.0, -3.0)
        @test jacobian(f, [1.5]) == 3.0
        @test jacobian(f, [2.5]) == 1.1413016497063289
        @test jacobian.(f, [[1.5], [2.5]]) == [3.0, 1.1413016497063289]

        # check jacobian function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 0.0)
        @test jacobian(f, [-4.0, 3.0]) == [4.0 4.0; 6.0 9.0]
        @test jacobian(f, [-2.0, 1.5]) == [4.0 4.0; 6.0 9.0]
        @test jacobian.(f, [[-4.0, 3.0], [-2.0, 1.5]]) == [[4.0 4.0; 6.0 9.0], [4.0 4.0; 6.0 9.0]]

        # check jacobian! function (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 0.0)
        output = zeros(2,2)
        jacobian!(output, f, [-4.0, 3.0])
        @test output == [4.0 4.0; 6.0 9.0]
        jacobian!(output, f, [-2.0, 1.5])
        @test output == [4.0 4.0; 6.0 9.0]

    end

    @testset "Utility Jacobian" begin

        # check utility functions jacobian (univariate)
        f = PlanarFlow(1.0, 2.0, -3.0)
        @test det_jacobian(f, 1.5) == 3.0
        @test inv_jacobian(f, 1.5) == 1/3
        @test absdet_jacobian(f, 1.5) == 3.0
        @test logabsdet_jacobian(f, 1.5) == log(3.0)

        # check utility functions jacobian (multivariate)
        f = PlanarFlow([1.0, 2.0], [3.0, 4.0], 0.0)
        @test det_jacobian(f, [-4.0, 3.0]) == 12.0
        @test inv_jacobian(f, [-4.0, 3.0]) ≈ [0.75 -1/3; -0.5 1/3]
        @test absdet_jacobian(f, [-4.0, 3.0]) == 12.0
        @test logabsdet_jacobian(f, [-4.0, 3.0]) == (log(12.0), 1)

    end

end

end


text:
# This file adds a new syntax to call elementwise operations on Convex expressions
# because overloading broadcast no longer works in Julia v0.6

import LinearAlgebra.dot
export dot

# multiplication
dot(::typeof(*)) = applyDotMultiplyAtom
applyDotMultiplyAtom(x, y) = broadcast(*, x, y)

# division
dot(::typeof(/)) = applyDotDivideAtom
applyDotDivideAtom(x, y) = broadcast(/, x, y)

# exponentiation
dot(::typeof(^)) = applyDotExpAtom
applyDotExpAtom(x, y) = broadcast(^, x, y)

# addition and subtraction are no problem since + and .+ should do the same thing

# alternative syntax

# dot(::typeof(*), x::Variable, y::Variable) = broadcast(*, x, y)
# dot(::typeof(/), x::Variable, y::Variable) = broadcast(/, x, y)
# dot(::typeof(+), x::Variable, y::Variable) = broadcast(+, x, y)
# dot(::typeof(-), x::Variable, y::Variable) = broadcast(-, x, y)
# dot(::typeof(^), x::Variable, y::Variable) = broadcast(^, x, y)


text:
#=
BFSLevelsMultiSource:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#
using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSLevelsMultiSource.jl")
dim = 20000
sourcdim = 3
matrix =  sprand(Bool, dim, dim, 0.005)
matrix = Array(matrix)
matrix = GBMatrix(matrix)
sources =  sprand(Bool, sourcdim, dim, 0.15)
sources = Array(sources)
sources = GBMatrix(sources)
@time ms_bfsl(matrix, sources, dim, sourcdim)


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

## QuasiDiagonal matrices

struct QuasiDiagonal{T,V<:AbstractQuasiVector{T}} <: AbstractQuasiMatrix{T}
    diag::V

    function QuasiDiagonal{T,V}(diag) where {T,V<:AbstractQuasiVector{T}}
        new{T,V}(diag)
    end
end
QuasiDiagonal(v::AbstractQuasiVector{T}) where {T} = QuasiDiagonal{T,typeof(v)}(v)
QuasiDiagonal{T}(v::AbstractQuasiVector) where {T} = QuasiDiagonal(convert(AbstractQuasiVector{T}, v)::AbstractQuasiVector{T})

"""
    QuasiDiagonal(A::AbstractQuasiMatrix)

Construct a matrix from the diagonal of `A`.

# Examples
```jldoctest
julia> A = [1 2 3; 4 5 6; 7 8 9]
3×3 Array{Int64,2}:
 1  2  3
 4  5  6
 7  8  9

julia> QuasiDiagonal(A)
3×3 QuasiDiagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9
```
"""
QuasiDiagonal(A::AbstractQuasiMatrix) = QuasiDiagonal(diag(A))

@deprecate Diagonal(A::AbstractQuasiArray) QuasiDiagonal(A)

"""
    QuasiDiagonal(V::AbstractQuasiVector)

Construct a matrix with `V` as its diagonal.

# Examples
```jldoctest
julia> V = [1, 2]
2-element Array{Int64,1}:
 1
 2

julia> QuasiDiagonal(V)
2×2 QuasiDiagonal{Int64,Array{Int64,1}}:
 1  ⋅
 ⋅  2
```
"""
QuasiDiagonal(V::AbstractQuasiVector)

QuasiDiagonal(D::QuasiDiagonal) = D
QuasiDiagonal{T}(D::QuasiDiagonal{T}) where {T} = D
QuasiDiagonal{T}(D::QuasiDiagonal) where {T} = QuasiDiagonal{T}(D.diag)

AbstractQuasiMatrix{T}(D::QuasiDiagonal) where {T} = QuasiDiagonal{T}(D)
AbstractQuasiArray{T}(D::QuasiDiagonal) where T = AbstractQuasiMatrix{T}(D)

convert(::Type{QuasiDiagonal}, A::QuasiDiagonal) = A
convert(::Type{QuasiDiagonal{T}}, A::QuasiDiagonal{T}) where T = A
convert(::Type{QuasiDiagonal{T}}, A::QuasiDiagonal) where T = QuasiDiagonal(convert(AbstractQuasiVector{T}, A.diag))

# For D<:QuasiDiagonal, similar(D[, neweltype]) should yield a QuasiDiagonal matrix.
# On the other hand, similar(D, [neweltype,] shape...) should yield a sparse matrix.
# The first method below effects the former, and the second the latter.
similar(D::QuasiDiagonal, ::Type{T}) where {T} = QuasiDiagonal(similar(D.diag, T))
# The method below is moved to SparseArrays for now
# similar(D::QuasiDiagonal, ::Type{T}, dims::Union{Dims{1},Dims{2}}) where {T} = spzeros(T, dims...)

copyto!(D1::QuasiDiagonal, D2::QuasiDiagonal) = (copyto!(D1.diag, D2.diag); D1)

size(D::QuasiDiagonal) = (cardinality(D.diag),cardinality(D.diag))
axes(D::QuasiDiagonal) = (axes(D.diag,1), axes(D.diag,1))

function size(D::QuasiDiagonal,d::Integer)
    if d<1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    end
    return d<=2 ? length(D.diag) : 1
end

@inline function getindex(D::QuasiDiagonal, i::Number, j::Number)
    @boundscheck checkbounds(D, i, j)
    if i == j
        @inbounds r = D.diag[i]
    else
        r = diagzero(D, i, j)
    end
    r
end
diagzero(::QuasiDiagonal{T},i,j) where {T} = zero(T)
diagzero(D::QuasiDiagonal{Matrix{T}},i,j) where {T} = zeros(T, size(D.diag[i], 1), size(D.diag[j], 2))

function setindex!(D::QuasiDiagonal, v, i::Int, j::Int)
    @boundscheck checkbounds(D, i, j)
    if i == j
        @inbounds D.diag[i] = v
    elseif !iszero(v)
        throw(ArgumentError("cannot set off-diagonal entry ($i, $j) to a nonzero value ($v)"))
    end
    return D
end

## structured matrix methods ##
function Base.replace_in_print_matrix(A::QuasiDiagonal,i::Integer,j::Integer,s::AbstractString)
    i==j ? s : Base.replace_with_centered_mark(s)
end

parent(D::QuasiDiagonal) = D.diag

ishermitian(D::QuasiDiagonal{<:Real}) = true
ishermitian(D::QuasiDiagonal{<:Number}) = isreal(D.diag)
ishermitian(D::QuasiDiagonal) = all(ishermitian, D.diag)
issymmetric(D::QuasiDiagonal{<:Number}) = true
issymmetric(D::QuasiDiagonal) = all(issymmetric, D.diag)
isposdef(D::QuasiDiagonal) = all(isposdef, D.diag)

factorize(D::QuasiDiagonal) = D

real(D::QuasiDiagonal) = QuasiDiagonal(real(D.diag))
imag(D::QuasiDiagonal) = QuasiDiagonal(imag(D.diag))

iszero(D::QuasiDiagonal) = all(iszero, D.diag)
isone(D::QuasiDiagonal) = all(isone, D.diag)
isdiag(D::QuasiDiagonal) = all(isdiag, D.diag)
isdiag(D::QuasiDiagonal{<:Number}) = true
istriu(D::QuasiDiagonal) = true
istril(D::QuasiDiagonal) = true
function triu!(D::QuasiDiagonal,k::Integer=0)
    n = size(D,1)
    if !(-n + 1 <= k <= n + 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n + 1) and at most $(n + 1) in an $n-by-$n matrix")))
    elseif k > 0
        fill!(D.diag,0)
    end
    return D
end

function tril!(D::QuasiDiagonal,k::Integer=0)
    n = size(D,1)
    if !(-n - 1 <= k <= n - 1)
        throw(ArgumentError(string("the requested diagonal, $k, must be at least ",
            "$(-n - 1) and at most $(n - 1) in an $n-by-$n matrix")))
    elseif k < 0
        fill!(D.diag,0)
    end
    return D
end

(==)(Da::QuasiDiagonal, Db::QuasiDiagonal) = Da.diag == Db.diag
(-)(A::QuasiDiagonal) = QuasiDiagonal(-A.diag)
(+)(Da::QuasiDiagonal, Db::QuasiDiagonal) = QuasiDiagonal(Da.diag + Db.diag)
(-)(Da::QuasiDiagonal, Db::QuasiDiagonal) = QuasiDiagonal(Da.diag - Db.diag)
(/)(D::QuasiDiagonal, x::Number) = QuasiDiagonal(D.diag / x)

conj(D::QuasiDiagonal) = QuasiDiagonal(conj(D.diag))
transpose(D::QuasiDiagonal{<:Number}) = D
transpose(D::QuasiDiagonal) = QuasiDiagonal(transpose.(D.diag))
adjoint(D::QuasiDiagonal{<:Number}) = conj(D)
adjoint(D::QuasiDiagonal) = QuasiDiagonal(adjoint.(D.diag))

function diag(D::QuasiDiagonal, k::Integer=0)
    # every branch call similar(..., ::Int) to make sure the
    # same vector type is returned independent of k
    if k == 0
        return copyto!(similar(D.diag, length(D.diag)), D.diag)
    elseif -size(D,1) <= k <= size(D,1)
        return fill!(similar(D.diag, size(D,1)-abs(k)), 0)
    else
        throw(ArgumentError(string("requested diagonal, $k, must be at least $(-size(D, 1)) ",
            "and at most $(size(D, 2)) for an $(size(D, 1))-by-$(size(D, 2)) matrix")))
    end
end
tr(D::QuasiDiagonal) = sum(D.diag)
det(D::QuasiDiagonal) = prod(D.diag)
logdet(D::QuasiDiagonal{<:Real}) = sum(log, D.diag)
function logdet(D::QuasiDiagonal{<:Complex}) # make sure branch cut is correct
    z = sum(log, D.diag)
    complex(real(z), rem2pi(imag(z), RoundNearest))
end

# Matrix functions
for f in (:exp, :log, :sqrt,
          :cos, :sin, :tan, :csc, :sec, :cot,
          :cosh, :sinh, :tanh, :csch, :sech, :coth,
          :acos, :asin, :atan, :acsc, :asec, :acot,
          :acosh, :asinh, :atanh, :acsch, :asech, :acoth)
    @eval $f(D::QuasiDiagonal) = QuasiDiagonal($f.(D.diag))
end

permutedims(D::QuasiDiagonal) = D
inv(D::QuasiDiagonal) = QuasiDiagonal(inv.(D.diag))
copy(D::QuasiDiagonal) = QuasiDiagonal(copy(D.diag))


text:
import Oceananigans.TimeSteppers: calculate_tendencies!

using Oceananigans: fields, prognostic_fields
using Oceananigans.Utils: work_layout

"""
    calculate_tendencies!(model::IncompressibleModel)

Calculate the interior and boundary contributions to tendency terms without the
contribution from non-hydrostatic pressure.
"""
function calculate_tendencies!(model::HydrostaticFreeSurfaceModel)

    # Calculate contributions to momentum and tracer tendencies from fluxes and volume terms in the
    # interior of the domain
    calculate_hydrostatic_free_surface_interior_tendency_contributions!(model.timestepper.Gⁿ,
                                                                        model.architecture,
                                                                        model.grid,
                                                                        model.advection,
                                                                        model.coriolis,
                                                                        model.buoyancy,
                                                                        model.closure,
                                                                        model.velocities,
                                                                        model.free_surface,
                                                                        model.tracers,
                                                                        model.pressure.pHY′,
                                                                        model.diffusivities,
                                                                        model.auxiliary_fields,
                                                                        model.forcing,
                                                                        model.clock)

    # Calculate contributions to momentum and tracer tendencies from user-prescribed fluxes across the
    # boundaries of the domain
    calculate_hydrostatic_boundary_tendency_contributions!(model.timestepper.Gⁿ,
                                                           model.architecture,
                                                           model.velocities,
                                                           model.free_surface,
                                                           model.tracers,
                                                           model.clock,
                                                           fields(model))

    return nothing
end

function calculate_hydrostatic_momentum_tendencies!(tendencies, velocities, arch, grid, advection, coriolis, closure,
                                                    free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                                    auxiliary_fields, forcings, clock, barrier)

    Gu_event = launch!(arch, grid, :xyz, calculate_hydrostatic_free_surface_Gu!,
                       tendencies.u, grid, advection.momentum, coriolis, closure,
                       velocities, free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                       auxiliary_fields, forcings, clock; dependencies = barrier)

    Gv_event = launch!(arch, grid, :xyz, calculate_hydrostatic_free_surface_Gv!,
                       tendencies.v, grid, advection.momentum, coriolis, closure,
                       velocities, free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                       auxiliary_fields, forcings, clock; dependencies = barrier)

    Gη_event = launch!(arch, grid, :xy, calculate_hydrostatic_free_surface_Gη!,
                       tendencies.η, grid, velocities, free_surface, tracers,
                       auxiliary_fields, forcings, clock; dependencies = barrier)

    events = [Gu_event, Gv_event, Gη_event]

    return events
end

""" Store previous value of the source term and calculate current source term. """
function calculate_hydrostatic_free_surface_interior_tendency_contributions!(tendencies,
                                                                             arch,
                                                                             grid,
                                                                             advection,
                                                                             coriolis,
                                                                             buoyancy,
                                                                             closure,
                                                                             velocities,
                                                                             free_surface,
                                                                             tracers,
                                                                             hydrostatic_pressure_anomaly,
                                                                             diffusivities,
                                                                             auxiliary_fields,
                                                                             forcings,
                                                                             clock)

    barrier = Event(device(arch))

    events = calculate_hydrostatic_momentum_tendencies!(tendencies, velocities, arch, grid, advection, coriolis, closure,
                                                        free_surface, tracers, diffusivities, hydrostatic_pressure_anomaly,
                                                        auxiliary_fields, forcings, clock, barrier)

    for (tracer_index, tracer_name) in enumerate(propertynames(tracers))
        @inbounds c_tendency = tendencies[tracer_name]
        @inbounds c_advection = advection[tracer_name]
        @inbounds forcing = forcings[tracer_name]

        Gc_event = launch!(arch, grid, :xyz, calculate_hydrostatic_free_surface_Gc!,
                           c_tendency, grid, Val(tracer_index),
                           c_advection, closure, buoyancy,
                           velocities, free_surface, tracers, diffusivities,
                           auxiliary_fields, forcing, clock;
                           dependencies=barrier)

        push!(events, Gc_event)
    end

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end

#####
##### Tendency calculators for u, v
#####

""" Calculate the right-hand-side of the u-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gu!(Gu, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gu[i, j, k] = hydrostatic_free_surface_u_velocity_tendency(i, j, k, grid, args...)
end

""" Calculate the right-hand-side of the v-velocity equation. """
@kernel function calculate_hydrostatic_free_surface_Gv!(Gv, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gv[i, j, k] = hydrostatic_free_surface_v_velocity_tendency(i, j, k, grid, args...)
end

#####
##### Tendency calculators for tracers
#####

""" Calculate the right-hand-side of the tracer advection-diffusion equation. """
@kernel function calculate_hydrostatic_free_surface_Gc!(Gc, grid, args...)
    i, j, k = @index(Global, NTuple)
    @inbounds Gc[i, j, k] = hydrostatic_free_surface_tracer_tendency(i, j, k, grid, args...)
end

#####
##### Tendency calculators for an explicit free surface
#####

""" Calculate the right-hand-side of the free surface displacement (η) equation. """
@kernel function calculate_hydrostatic_free_surface_Gη!(Gη, grid, args...)
    i, j = @index(Global, NTuple)
    @inbounds Gη[i, j, 1] = free_surface_tendency(i, j, grid, args...)
end

#####
##### Boundary condributions to hydrostatic free surface model
#####

function apply_flux_bcs!(Gcⁿ, events, c, arch, barrier, clock, model_fields)
    x_bcs_event = apply_x_bcs!(Gcⁿ, c, arch, barrier, clock, model_fields)
    y_bcs_event = apply_y_bcs!(Gcⁿ, c, arch, barrier, clock, model_fields)
    z_bcs_event = apply_z_bcs!(Gcⁿ, c, arch, barrier, clock, model_fields)

    push!(events, x_bcs_event, y_bcs_event, z_bcs_event)

    return nothing
end

""" Apply boundary conditions by adding flux divergences to the right-hand-side. """
function calculate_hydrostatic_boundary_tendency_contributions!(Gⁿ, arch, velocities, free_surface, tracers, clock, model_fields)

    barrier = Event(device(arch))

    events = []

    # Velocity fields
    for i in (:u, :v)
        apply_flux_bcs!(Gⁿ[i], events, velocities[i], arch, barrier, clock, model_fields)
    end

    # Free surface
    apply_flux_bcs!(Gⁿ.η, events, displacement(free_surface), arch, barrier, clock, model_fields)

    # Tracer fields
    for i in propertynames(tracers)
        apply_flux_bcs!(Gⁿ[i], events, tracers[i], arch, barrier, clock, model_fields)
    end

    events = filter(e -> typeof(e) <: Event, events)

    wait(device(arch), MultiEvent(Tuple(events)))

    return nothing
end


text:
using wiki2
using Test

@testset "wiki2.jl" begin
    # Write your tests here.
end