text:
# https://adrhill.github.io/ExplainableAI.jl/stable/generated/advanced_lrp/#How-it-works-internally
abstract type AbstractLRPRule end

# Generic LRP rule. Since it uses autodiff, it is used as a fallback for layer types without custom implementation.
function lrp!(rule::R, layer::L, Rₖ, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule,L}
    lrp_autodiff!(rule, layer, Rₖ, aₖ, Rₖ₊₁)
    return nothing
end

function lrp_autodiff!(
    rule::R, layer::L, Rₖ::T1, aₖ::T1, Rₖ₊₁::T2
) where {R<:AbstractLRPRule,L,T1,T2}
    layerᵨ = modify_layer(rule, layer)
    c::T1 = only(
        gradient(aₖ) do a
            z::T2 = layerᵨ(a)
            s = Zygote.@ignore Rₖ₊₁ ./ modify_denominator(rule, z)
            z ⋅ s
        end,
    )
    Rₖ .= aₖ .* c
    return nothing
end

# For linear layer types such as Dense layers, using autodiff is overkill.
function lrp!(rule::R, layer::Dense, Rₖ, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule}
    lrp_dense!(rule, layer, Rₖ, aₖ, Rₖ₊₁)
    return nothing
end

function lrp_dense!(rule::R, l, Rₖ, aₖ, Rₖ₊₁) where {R<:AbstractLRPRule}
    ρW, ρb = modify_params(rule, get_params(l)...)
    ãₖ₊₁ = modify_denominator(rule, ρW * aₖ .+ ρb)
    @tullio Rₖ[j, b] = aₖ[j, b] * ρW[k, j] / ãₖ₊₁[k, b] * Rₖ₊₁[k, b]
    return nothing
end

# Other special cases that are dispatched on layer type:
lrp!(::AbstractLRPRule, ::DropoutLayer, Rₖ, aₖ, Rₖ₊₁) = (Rₖ .= Rₖ₊₁)
lrp!(::AbstractLRPRule, ::ReshapingLayer, Rₖ, aₖ, Rₖ₊₁) = (Rₖ .= reshape(Rₖ₊₁, size(aₖ)))

# To implement new rules, we can define two custom functions `modify_params` and `modify_denominator`.
# If this isn't done, the following fallbacks are used by default:
"""
    modify_params(rule, W, b)

Function that modifies weights and biases before applying relevance propagation.
Returns modified weights and biases as a tuple `(ρW, ρb)`.
"""
modify_params(::AbstractLRPRule, W, b) = (W, b) # general fallback

"""
    modify_denominator(rule, d)

Function that modifies zₖ on the forward pass, e.g. for numerical stability.
"""
modify_denominator(::AbstractLRPRule, d) = stabilize_denom(d, 1.0f-9) # general fallback

"""
    modify_layer(rule, layer)

Function that modifies a layer before applying relevance propagation.
Returns a new, modified layer.
"""
modify_layer(::AbstractLRPRule, layer) = layer # skip layers without modify_params
function modify_layer(rule::R, layer::L) where {R<:AbstractLRPRule,L<:Union{Dense,Conv}}
    return set_params(layer, modify_params(rule, get_params(layer)...)...)
end

"""
    ZeroRule()

Constructor for LRP-0 rule. Commonly used on upper layers.
"""
struct ZeroRule <: AbstractLRPRule end

"""
    GammaRule(; γ=0.25)

Constructor for LRP-``γ`` rule. Commonly used on lower layers.

Arguments:
- `γ`: Optional multiplier for added positive weights, defaults to 0.25.
"""
struct GammaRule{T} <: AbstractLRPRule
    γ::T
    GammaRule(; γ=0.25) = new{Float32}(γ)
end
function modify_params(r::GammaRule, W, b)
    ρW = W + r.γ * relu.(W)
    ρb = b + r.γ * relu.(b)
    return ρW, ρb
end

"""
    EpsilonRule(; ϵ=1f-6)

Constructor for LRP-``ϵ`` rule. Commonly used on middle layers.

Arguments:
- `ϵ`: Optional stabilization parameter, defaults to `1f-6`.
"""
struct EpsilonRule{T} <: AbstractLRPRule
    ϵ::T
    EpsilonRule(; ϵ=1.0f-6) = new{Float32}(ϵ)
end
modify_denominator(r::EpsilonRule, d) = stabilize_denom(d, r.ϵ)

"""
    ZBoxRule()

Constructor for LRP-``z^{\\mathcal{B}}``-rule.
Commonly used on the first layer for pixel input.
"""
struct ZBoxRule <: AbstractLRPRule end

# The ZBoxRule requires its own implementation of relevance propagation.
lrp!(::ZBoxRule, layer::Dense, Rₖ, aₖ, Rₖ₊₁) = lrp_zbox!(layer, Rₖ, aₖ, Rₖ₊₁)
lrp!(::ZBoxRule, layer::Conv, Rₖ, aₖ, Rₖ₊₁) = lrp_zbox!(layer, Rₖ, aₖ, Rₖ₊₁)

function lrp_zbox!(layer::L, Rₖ::T1, aₖ::T1, Rₖ₊₁::T2) where {L,T1,T2}
    W, b = get_params(layer)
    l, h = fill.(extrema(aₖ), (size(aₖ),))

    layer⁺ = set_params(layer, max.(0, W), max.(0, b)) # W⁺, b⁺
    layer⁻ = set_params(layer, min.(0, W), min.(0, b)) # W⁻, b⁻

    c::T1, cₗ::T1, cₕ::T1 = gradient(aₖ, l, h) do a, l, h
        f::T2 = layer(a)
        f⁺::T2 = layer⁺(l)
        f⁻::T2 = layer⁻(h)

        z = f - f⁺ - f⁻
        s = Zygote.@ignore safedivide(Rₖ₊₁, z, 1e-9)
        z ⋅ s
    end
    Rₖ .= aₖ .* c + l .* cₗ + h .* cₕ
    return nothing
end

ext: jl


text:
function MMI.fit(detector::UnsupervisedDetector, verbosity::Int, X)
    model = fit(detector, X)
    model, nothing, (scores = model.scores,)
end

function MMI.fit(detector::SupervisedDetector, verbosity::Int, X, y)
    model = fit(detector, X, y)
    model, nothing, (scores = model.scores,)
end

function MMI.transform(detector::Detector, fitresult::Fit, X)
    score(detector, fitresult, X)
end

function MMI.predict(detector::Detector, fitresult::Fit, X)
    score(detector, fitresult, X)
end

function MMI.transform(ev::Class, _, scores::Result...) # _ because there is no fitresult
    MMI.categorical(detect(ev, scores...))
end

function MMI.transform(ev::Score, _, scores::Result...) # _ because there is no fitresult
    scores = detect(ev, scores...) # if not normalized, return the raw scores, otherwise return a distribution
    isnothing(ev.normalize) ? scores : MMI.UnivariateFinite([CLASS_NORMAL, CLASS_OUTLIER],
                                                            hcat(1 .- scores, scores); pool=missing, ordered=false)
end

# specify scitypes
MMI.input_scitype(::Type{<:Detector}) = Union{MMI.Table(MMI.Continuous),AbstractMatrix{MMI.Continuous}}
MMI.output_scitype(::Type{<:Detector}) = Tuple{AbstractVector{<:MMI.Continuous},AbstractVector{<:MMI.Continuous}}

# data front-end for fit (supervised):
MMI.reformat(::SupervisedDetector, X, y) = (MMI.matrix(X, transpose=true), y)
MMI.reformat(::SupervisedDetector, X, y, w) = (MMI.matrix(X, transpose=true), y, w)
MMI.selectrows(::SupervisedDetector, I, Xmatrix, y) = (view(Xmatrix, :, I), view(y, I))
MMI.selectrows(::SupervisedDetector, I, Xmatrix, y, w) = (view(Xmatrix, :, I), view(y, I), view(w, I))

# data front-end for fit (unsupervised)/predict/transform
MMI.reformat(::Detector, X) = (MMI.matrix(X, transpose=true),)
MMI.selectrows(::Detector, I, Xmatrix) = (view(Xmatrix, :, I),)

MODELS = (ABOD, COF, DNN, KNN, LOF, AE, DeepSAD, ESAD, Class, Score)
MMI.metadata_pkg.(MODELS,
    package_name="OutlierDetection.jl",
    package_uuid="262411bb-c475-4342-ba9e-03b8c0183ca6",
    package_url="https://github.com/davnn/OutlierDetection.jl",
    is_pure_julia=true,
    package_license="MIT",
    is_wrapper=false)


text:
module Int

include("Computation.jl")
using .Computation

include("Format.jl")
using .Format

export int, vint

# Determines upper bound alpha with Interval Method, optionally outputs proof
# Authors: Antara Hebbar and Jason Liu
function int(m::Int64, s::Int64; output::Int64 = 2)
    output > 0 && printHeader(center("INTERVAL METHOD"))

    if m <= 0 || s <= 0
        if output > 0
            printf("Interval Method does not apply", line = true)
            printEnd()
        end
        return 1 // 1
    elseif m < s
        alpha = m // s * int(s, m, output = 0)
        if output > 1
            printfT(
                "Duality Theorem",
                "muffins($m,$s) = $m/$s · muffins($s,$m)",
                "Bounding muffins($s,$m) instead",
            )
            int(s, m, output = output)
        elseif output > 0
            printfT(
                "Interval Method",
                "Upper bound α of muffins($m,$s) is $(formatFrac(alpha))",
            )
            printEnd()
        end
        return alpha
    elseif m % s == 0
        if output > 0
            printfT("Interval Method", "Since $m % $s = 0, muffins($m,$s) = 1")
            printEnd()
        end
        return 1 // 1
    end

    (V, W, sV, sW) = sv(m, s)
    numV = (V)sV
    numW = (W)sW

    # Derive alpha using Int Algorithm
    if numW > numV
        f = Int64(floor(numV / sW))                                     # Upper bound of minimum # of large W-shs a student can have
        g = Int64(floor((numW - numV) / sW))                            # Lower bound of minimum # of small W-shs a student can have
        alpha = min(
            ((W - f)W - (W - f + 1)m // s + f) // ((W - f - 1)W + 2f),  # Value for alpha derived by solving f(1-α) + (W-f)(1-y) = m/s
            ((g - 1)W + (W - 2g + 1)m // s) // (W^2 - g),               # Value for alpha derived by solving gy + (W-g)(1-x) = m/s
        )
        alpha = alpha < 1 / 3 ? 1 // 3 : alpha
        if vint(m, s, alpha, output = output)
            return alpha
        end
    elseif numW < numV
        f = Int64(floor((numV - numW) / sV))                            # Upper bound of minimum # of large V-shs a student can have
        g = Int64(floor(numW / sV))                                     # Lower bound of minimum # of small V-shs a student can have
        alpha = min(
            ((V - f)W + (2f - V - 1)m // s) // ((V - f)W + (f - 1)V),   # Value for alpha derived by solving f(1-x) + (V-f)(1-y) = m/s
            ((V - g + 1)m // s + g - V) // ((V - g - 1)V + 2g),         # Value for alpha derived by solving gα + (V-g)(1-x) = m/s
        )
        alpha = alpha < 1 / 3 ? 1 // 3 : alpha
        if vint(m, s, alpha, output = output)
            return alpha
        end
    end

    if output > 0
        printf("Interval Method inconclusive", line = true)
        printEnd()
    end
    1 // 1
end

# Helper function for int -- verifies whether int(m, s, alpha) is conclusive
function vint(m::Int64, s::Int64, alpha::Rational{Int64}; output::Int64 = 2)
    if !(m > s > 0) || m % s == 0
        output > 0 && printf("VInt does not apply", line = true)
        false
    elseif alpha < 1 / 3
        output > 0 && printfT("No piece size < 1/3", "For m ≥ s, α must be ≥ 1/3")
        false
    elseif alpha > 1
        output > 0 && printfT("No piece size > 1", "α must be ≤ 1")
        false
    else
        (V, W, sV, sW) = sv(m, s)
        numV = (V)sV
        numW = (W)sW

        # Initialize Interval Method proof
        if output > 1
            # Define and format variables for proof
            alphaF = formatFrac(alpha)
            size = formatFrac(m // s)
            a = formatFrac(m // (V + 1)s)
            b = max(1 - m // (W - 1)s, 0)
            b1 = formatFrac(1 - b)
            b = formatFrac(b)

            # Establish assumptions and premises
            printHeader("CLAIM")
            printfT(
                "Claim",
                "For all ($m, $s) procedures where $s students each receive $size muffins, the smallest piece size is bounded above by α = $alphaF",
                "Proven by contradicting the assumption",
            )

            printHeader("ASSUMPTION")
            printfT("Assumption", "Assume muffins($m,$s) > α ≥ 1/3")
            printfT(
                "Property of Buddies",
                "Since muffins($m,$s) > 1/3, each muffin must be cut into 2 shs, totaling $(2m) shs",
                "",
                "Since each muffin is cut into 2 shs that are buddies, each sh of size x must imply the existence of a sh of size 1-x",
            )

            # Begin casework
            printHeader("CASEWORK")

            printfT(
                "Case 1",
                "Alice has ≥ $(V+1) shs",
                "One of them is ≤ ($size × 1/$(V+1)) = $a",
                "",
                "Contradicts assumption if α ≥ $a",
            )

            if W - 1 <= 1 && m / s > 1
                printfT(
                    "Case 2",
                    "Bob cannot have ≤ $(W-1) shs since $size > 1, so this case is impossible",
                )
            else
                printfT(
                    "Case 2",
                    "Bob has ≤ $(W-1) shs",
                    "One of them is ≥ ($size × 1/$(W-1)) = $b1",
                    "Its buddy is ≤ (1 - $b1) = $b",
                    "",
                    "Contradicts assumption if α ≥ $b",
                )
            end
        end

        # Check if V-Conjecture applies
        if m // s * 1 // (V + 1) > alpha || 1 - m // s * 1 // (W - 1) > alpha
            output > 0 && printf("V-Conjecture does not apply, VInt failed", line = true)
            return false
        end

        ((_, x), (y, _)) = findend(m, s, alpha, V)
        xOriginal = m // s - alpha * (V - 1)
        yOriginal = m // s - (1 - alpha) * (V - 2)

        # Check if FindEnd works
        if x != xOriginal && y != yOriginal
            output > 0 && printf("FindEnd inconclusive, VInt failed", line = true)
            return false
        end

        # Define and format variables for proof
        cd = lcm(s, denominator(x), denominator(y))
        alphaF = formatFrac(alpha, cd)
        alpha1 = formatFrac(1 - alpha, cd)
        size = formatFrac(m // s)
        a = formatFrac(m // s - x, cd)
        b = formatFrac(m // s - y, cd)
        xF = formatFrac(x, cd)
        x1 = formatFrac(1 - x, cd)
        yF = formatFrac(y, cd)
        y1 = formatFrac(1 - y, cd)
        diff = abs(numV - numW)

        # Continue proof
        if output > 1
            # Continue casework
            printfT(
                "Note",
                "The remaining cases deal with everyone having either $W or $V shs, so:",
                "",
                "$(W)·s_$W + $(V)·s_$V = $(2m)  (total shs)",
                "s_$W + s_$V = $s  (total students)",
                "where s_N = # of students with N shs",
                "",
                "The solution to the system is:",
                "s_$W = $sW, s_$V = $sV",
                "So there are $numW $W-shs and $numV $V-shs",
            )

            if x == xOriginal
                printfT(
                    "Case 3",
                    "Alice has a $V-sh ≥ $xF",
                    "Her other $(V-1) $V-shs sum to ≤ ($size - $xF) = $a",
                    "One of them is ≤ ($a × 1/$(V-1)) = $alphaF",
                    "",
                    "Contradicts assumption if α ≥ $alphaF",
                )
            else
                printfT("Case 3", "FindEnd did not produce a conclusive bound for $V-shs")
            end

            if y == yOriginal
                printfT(
                    "Case 4",
                    "Bob has a $W-sh ≤ $yF",
                    "His other $(W-1) $W-shs sum to ≥ ($size - $yF) = $b",
                    "One of them is ≥ ($b × 1/$(W-1)) = $alpha1",
                    "Its buddy is ≤ (1 - $alpha1) = $alphaF",
                    "",
                    "Contradicts assumption if α ≥ $alphaF",
                )
            else
                printfT("Case 4", "FindEnd did not produce a conclusive bound for $W-shs")
            end

            printHeader("CASE 5: INTERVAL ANALYSIS")
            println()
            printf("The following intervals capture the negation of the previous cases:")
            if alpha < x <= y < 1 - alpha
                println(
                    "\n",
                    interval(
                        ["(", alphaF],
                        [")[", xF],
                        ["](", yF],
                        [")", alpha1],
                        labels = ["$numV $V-shs", "0", "$numW $W-shs"],
                    ),
                )
                printLine()
            else
                if xF != alphaF
                    println(
                        "\n",
                        interval(
                            ["(", alphaF],
                            [")[", xF],
                            [")", alpha1],
                            labels = ["$numV $V-shs", "0 $V-shs"],
                        ),
                    )
                else
                    println("\n", "No negation interval for $V-shs")
                end

                if yF != alpha1
                    println(
                        "\n",
                        interval(
                            ["(", alphaF],
                            ["](", yF],
                            [")", alpha1],
                            labels = ["0 $W-shs", "$numW $W-shs"],
                        ),
                    )
                else
                    println("\n", "No negation interval for $W-shs")
                end

                printfT(
                    "Case 5",
                    "The Interval Method is inconclusive on these intervals, VInt failed",
                )
            end
        end

        # Fail if interval inconclusive
        if x > y || x == alpha || y == 1 - alpha
            output == 1 &&
                printf("Could not generate disjoint intervals, VInt failed", line = true)
            return false
        end

        # Define and format variables for proof
        (vMin, vMax, sMin, sMax) = (V, W, sV, sW)
        (f, g) = (Int64(floor(numV / sW)), Int64(floor(diff / sW)))
        (lim1, lim2) = ((vMin)sMin, diff)
        (i, j, k, l) = (xF, yF, y1, x1)
        (rngMin, rngMax, rngS, rngL) = ((alphaF, xF), (x1, alpha1), (yF, y1), (x1, alpha1))
        if numV > numW
            (vMin, vMax, sMin, sMax) = (W, V, sW, sV)
            (f, g) = (Int64(floor(diff / sV)), Int64(floor(numW / sV)))
            (lim1, lim2) = (diff, (vMin)sMin)
            (i, j, k, l) = (y1, x1, xF, yF)
            (rngMin, rngMax, rngS, rngL) =
                ((yF, alpha1), (alphaF, y1), (alphaF, y1), (x1, xF))
        end
        numMin = (vMin)sMin

        # Conclude interval case
        if output > 1
            printfT(
                "Property of Buddies",
                "Because 0 shs lie in [$xF,$yF], 0 shs must also lie in [$y1,$x1]",
                "",
                "Similary, the $numMin $vMin-shs in ($(rngMin[1]),$(rngMin[2])) implies $numMin $vMax-shs in ($(rngMax[1]),$(rngMax[2]))",
            )
            println()
            printf("The following intervals capture the previous statements:")
            println(
                "\n",
                interval(
                    ["(", alphaF],
                    [")[", i],
                    ["](", j],
                    [")[", k],
                    ["](", l],
                    [")", alpha1],
                    labels = [
                        "$numMin $V-shs",
                        "0",
                        "$diff $vMax-shs",
                        "0",
                        "$numMin $W-shs",
                    ],
                ),
            )
            printLine()

            printfT(
                "Note",
                "Let the $vMax-shs that lie in ($(rngS[1]),$(rngS[2])) be small $vMax-shs",
                "Let the $vMax-shs that lie in ($(rngL[1]),$(rngL[2])) be large $vMax-shs",
            )
        end

        upperB = (vMax - f) * toFrac(rngS[2]) + (f) * toFrac(rngL[2])
        lowerB = (vMax - g) * toFrac(rngL[1]) + (g) * toFrac(rngS[1])
        u = formatFrac(upperB, cd)
        l = formatFrac(lowerB, cd)

        if upperB <= m // s
            if output > 1
                printfT(
                    "Case 5.1",
                    "Alice has ≤ $f large $vMax-shs",
                    "She must also have ≥ $(vMax-f) small $vMax-shs, so her maximum amount of muffins is:",
                    "< ($(vMax-f) × $(rngS[2])) + ($f × $(rngL[2])) = $u ≤ $size",
                    "",
                    "Since Alice's muffin amount can not reach $size, this case is impossible",
                )
                printfT(
                    "Case 5.2",
                    "All $sMax $vMax-sh students have at least $(f+1) large $vMax-shs",
                    "",
                    "This is impossible because $(f+1) × $sMax = $((f+1)sMax) ≥ $lim1",
                )
            end
        elseif lowerB >= m // s
            if output > 1
                printfT(
                    "Case 5.1",
                    "Bob has ≤ $g small $vMax-shs",
                    "He must also have ≥ $(vMax-g) large $vMax-shs, so his minimum amount of muffins is:",
                    "> ($g × $(rngS[1])) + ($(vMax-g) × $(rngL[1])) = $l ≥ $size",
                    "",
                    "Since Bob's muffin amount can not reach $size, this case is impossible",
                )
                printfT(
                    "Case 5.2",
                    "All $sMax $vMax-sh students have at least $(g+1) small $vMax-shs",
                    "",
                    "This is impossible because $(g+1) × $sMax = $((g+1)sMax) ≥ $lim2",
                )
            end
        else
            if output > 1
                printfT(
                    "Bound # of large $vMax-shs",
                    "All $sMax $vMax-sh students have at least $(f+1) large $vMax-shs",
                    "",
                    "This is impossible because $(f+1) × $sMax = $((f+1)sMax) ≥ $lim1",
                )
                printfT(
                    "Negation of prev. bound",
                    "∃ student Alice with ≤ $f large $vMax-shs",
                    "She must also have ≥ $(vMax-f) small $vMax-shs, so her maximum amount of muffins is:",
                    "< ($(vMax-f) × $(rngS[2])) + ($f × $(rngL[2])) = $u",
                    "",
                    "Since $size lies below the interval, this bound is inconclusive",
                )

                printfT(
                    "Bound # of small $vMax-shs",
                    "All $sMax $vMax-sh students have at least $(g+1) small $vMax-shs",
                    "",
                    "This is impossible because $(g+1) × $sMax = $((g+1)sMax) ≥ $lim2",
                )
                printfT(
                    "Negation of prev. bound",
                    "∃ student Bob with ≤ $g small $vMax-shs",
                    "He must also have ≥ $(vMax-g) large $vMax-shs, so his minimum amount of muffins is:",
                    "> ($g × $(rngS[1])) + ($(vMax-g) × $(rngL[1])) = $l",
                    "",
                    "Since $size lies above the interval, this bound is inconclusive",
                )

                printf(
                    "Bounding # of large and small $vMax-shs inconclusive, restarting with individual case analysis",
                    line = true,
                )
            end

            for k = 0:vMax
                upperB = (vMax - k) * toFrac(rngL[2]) + (k) * toFrac(rngS[2])
                lowerB = (vMax - k) * toFrac(rngL[1]) + (k) * toFrac(rngS[1])
                u = formatFrac(upperB, cd)
                l = formatFrac(lowerB, cd)
                if upperB <= m // s || lowerB >= m // s
                    if output > 1
                        printfT(
                            "Case 5.$(k+1)",
                            "Alice has $k small $vMax-shs and $(vMax-k) large $vMax-shs",
                            "Her possible muffin amount lies in:",
                            "",
                            "($k × $(rngS[1]) + $(vMax-k) × $(rngL[1]),",
                            "$k × $(rngS[2]) + $(vMax-k) × $(rngL[2]))",
                            "= ($l, $u)",
                            "",
                            "Since $size lies outside this interval, this case is impossible",
                        )
                    end
                else
                    # Fail if interval analysis inconclusive
                    if output > 1
                        printfT(
                            "Case 5.$(k+1)",
                            "Bob has $k small $vMax-shs and $(vMax-k) large $vMax-shs",
                            "His possible muffin amount lies in:",
                            "",
                            "($k × $(rngS[1]) + $(vMax-k) × $(rngL[1]),",
                            "$k × $(rngS[2]) + $(vMax-k) × $(rngL[2]))",
                            "= ($l, $u)",
                            "",
                            "Since $size lies inside this interval, this case is inconclusive",
                        )
                    elseif output > 0
                        printf("Interval analysis inconclusive, VInt failed", line = true)
                    end
                    return false
                end
            end
        end

        # Conclude proof
        alpha = formatFrac(alpha)
        if output > 1
            # Conclude with alpha's value
            printHeader("CONCLUSION")
            printfT(
                "Conclusion",
                "All possible cases contradict the assumption iff. α ≥ $alpha",
                "muffins($m,$s) ≤ α, ∀ α ≥ $alpha",
                "",
                "muffins($m,$s) ≤ $alpha",
            )
            printEnd()
        elseif output > 0
            printfT("Interval Method", "Upper bound α of muffins($m,$s) is $alpha")
            printEnd()
        end

        true
    end
end

end


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

sm = @stochastic_model begin
    @stage 1 begin
        @decision(model, x¹ <= 2)
        @decision(model, w¹)
        @decision(model, y¹)
        @objective(model, Min, 100*x¹ + 3*w¹ + 0.5*y¹)
        @constraint(model, x¹ + w¹ - y¹ == 1)
    end
    @stage 2 begin
        @uncertain d
        @variable(model, x² <= 2)
        @variable(model, w²)
        @decision(model, y²)
        @objective(model, Min, x² + 3*w² + 0.5*y²)
        @constraint(model, y¹ + x² + w² - y² == d)
    end
    @stage 3 begin
        @uncertain d
        @variable(model, x³ <= 2)
        @variable(model, w³)
        @recourse(model, y³)
        @objective(model, Min, x³ + 3*w³)
        @constraint(model, y² + x³ + w³ - y³ == d)
    end
end
s₂₁ = @scenario d = 1.0 probability = 0.5
s₂₂ = @scenario d = 3.0 probability = 0.5

s₃₁ = @scenario d = 1.0 probability = 0.25
s₃₂ = @scenario d = 3.0 probability = 0.25
s₃₃ = @scenario d = 1.0 probability = 0.25
s₃₄ = @scenario d = 3.0 probability = 0.25

sp = instantiate(sm, ([s₂₁,s₂₂], [s₃₁,s₃₂,s₃₃,s₃₄]))


text:
module RelativeBenchmarks

using Random
using LinearAlgebra

import MacroTools: postwalk
using MacroTools: MacroTools, prewalk, postwalk, @capture
using BenchmarkTools

#=
    You set this directory with the default const or through the environment.
    Ensure `NewImplementation.jl` and `OldImplementation.jl` are kept there.
=#
const DefaultBenchmarkingDir = pwd()
const BenchmarkingDir = haskey(ENV, "BenchmarkingDir") ? ENV["BenchmarkingDir"] : DefaultBenchmarkingDir

cd(BenchmarkingDir)
isfile("NewImplmentation.jl") && isfile("OldImplmentation.jl") ||
    error("Both `NewImplmentation.jl` and `OldImplementation.jl` must be found in the BenchmarkingDir ($BenchmarkingDir)")

include("setup.jl")
include("relspeeds.jl")

end  # RelativeBenchmarks


text:
# test stacks and queues

# Stack

s = Stack(Int, 5)
n = 100

@test isa(s, Stack{Int})
@test length(s) == 0
@test isempty(s)
@test_throws ArgumentError top(s)
@test_throws ArgumentError pop!(s)

for i = 1 : n
    push!(s, i)
    @test top(s) == i
    @test !isempty(s)
    @test length(s) == i
end

for i = 1 : n
    x = pop!(s)
    @test x == n - i + 1
    if i < n
        @test top(s) == n - i
    else
        @test_throws ArgumentError top(s)
    end
    @test isempty(s) == (i == n)
    @test length(s) == n - i
end

#test that iter returns a LIFO collection

stk = Stack(Int, 10)
#an array to check iteration sequence against
arr = Int64[]

for i = 1:n
    push!(stk,i)
    push!(arr,i)
end

#test iterator

index = length(arr)
for i in stk
    @test(arr[index] == i)
    index -= 1
end

index = 1
for i in reverse_iter(stk)
    @test(arr[index] == i)
    index += 1
end

@test arr == [i for i in reverse_iter(stk)]
@test reverse(arr) == [i for i in stk]

# Queue

s = Queue(Int, 5)
n = 100

@test isa(s, Queue{Int})
@test length(s) == 0
@test isempty(s)
@test_throws ArgumentError front(s)
@test_throws ArgumentError back(s)
@test_throws ArgumentError dequeue!(s)

for i = 1 : n
    enqueue!(s, i)
    @test front(s) == 1
    @test back(s) == i
    @test !isempty(s)
    @test length(s) == i
end

for i = 1 : n
    x = dequeue!(s)
    @test x == i
    if i < n
        @test front(s) == i + 1
        @test back(s) == n
    else
        @test_throws ArgumentError front(s)
        @test_throws ArgumentError back(s)
    end
    @test isempty(s) == (i == n)
    @test length(s) == n - i
end

#test that iter returns a FIFO collection

q = Queue(Int, 10)
#an array to check iteration sequence against
arr = Int64[]

for i = 1:n
    enqueue!(q,i)
    push!(arr,i)
end

#test iterator

index = 1
for i in q
    @test(arr[index] == i)
    index += 1
end

index = length(arr)
for i in reverse_iter(q)
    @test(arr[index] == i)
    index -= 1
end

@test arr == [i for i in q]
@test reverse(arr) == [i for i in reverse_iter(q)]


text:
#!/usr/bin/env julia
using RobotOS
@rosimport geometry_msgs.msg: Point, Pose, Pose2D, PoseStamped, Vector3, Twist
@rosimport std_srvs.srv: Empty, SetBool
@rosimport nav_msgs.srv.GetPlan
@rosimport gazebo_msgs.msg: ModelState, ModelStates
@rosimport gazebo_msgs.srv: SetModelState, GetModelState, GetWorldProperties

rostypegen()
using geometry_msgs.msg
using std_srvs.srv
using nav_msgs.srv
using nav_msgs.srv.GetPlan
using gazebo_msgs.srv
using gazebo_msgs.msg

function loop(get_state, pub)
    loop_rate = Rate(100.0)
    while ! is_shutdown()
        modelName = "hmmwv"  # TODO make this a parameter

        # Get the current position of the Gazebo model
        gs = GetModelStateRequest()
        gs.model_name = modelName
        gs_r = get_state(gs)
        if !gs_r.success
            error(string(" calling /gazebo/get_model_state service: ", gs_r.status_message))
        end

        npt = gs_r.pose
        publish(pub, npt)
        rossleep(loop_rate)
    end
end

function main()
    init_node("rosjl_position")

    # Set up service to get Gazebo model state
    const get_state = ServiceProxy("/gazebo/get_model_state",GetModelState)
    println("Waiting for 'gazebo/get_model_state' service...")
    wait_for_service("gazebo/get_model_state")
    # TODO use getparam to get robot name
    pub = Publisher{Pose}("/robot/pose", queue_size=10)

    loop(get_state, pub)
end

if ! isinteractive()
    main()
end


text:
#=
Returns an empty Associative Array
=#
function emptyAssoc()
    Assoc()
end

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
using Test

m = test_page_model()
include("../src/components/CO2forcing.jl")

add_comp!(m, co2forcing)

set_param!(m, :co2forcing, :c_CO2concentration, readpagedata(m,"test/validationdata/c_co2concentration.csv"))

##running Model
run(m)

forcing=m[:co2forcing,:f_CO2forcing]
forcing_compare=readpagedata(m,"test/validationdata/f_co2forcing.csv")

@test forcing ≈ forcing_compare rtol=1e-3


text:
################################################################################
#
#	ABSTRACT TYPE
#
#   BondHoppingHamiltonianSimpleNN <: AbstractBondHamiltonian{L,1}
#   --> L is the label type of bonds
#   --> N is the dimension the bond term matrix (NxN matrix)
#
#   FILE CONTAINS
#       - abstract type definition
#       - interface definition
#       - generator function
#
################################################################################

################################################################################
#
#   ABSTRACT TYPE DEFINITION
#
################################################################################
mutable struct BondHoppingHamiltonianSimpleNN{L} <: AbstractBondHamiltonian{L,1}

    # coupling
    coupling :: Float64

    # bond label that identifies NN
    label :: L

end

# export the type
export BondHoppingHamiltonianSimpleNN

################################################################################
#
#	INTERFACING / ACCESSING BOND HAMILTONIANS
#	(functions have to be overwritten by concrete types)
#
################################################################################

# get bond term of Hamiltonian
function bondterm(
            h :: BondHoppingHamiltonianSimpleNN{L},
            b :: AbstractBond{L,NB}
        ) :: Float64 where {N,L,NB,NS}

    # get the bond label
    l = label(b)

    # check if it is the NN label
    if l == h.label
        return h.coupling
    else
        return 0.0
    end
end

# export the function
export bondterm

function getHoppingHamiltonianSimpleNN(
            unitcell :: U
        ) :: BondHoppingHamiltonianSimpleNN{L} where {L,NB,S,B<:AbstractBond{L,NB},U<:AbstractUnitcell{S,B}}

    # obtain all couplings
    couplings = unique!(label.(bonds(unitcell)))

    # check if the list contains a default label
    if getDefaultLabel(L) in couplings
        return BondHoppingHamiltonianSimpleNN{L}(1.0, getDefaultLabel(L))
    else
        return BondHoppingHamiltonianSimpleNN{L}(1.0, couplings[1])
    end
end

export getHoppingHamiltonianSimpleNN


text:
#####
##### PrescribedVelocityFields
#####

using Oceananigans.Grids: Center, Face
using Oceananigans.Fields: AbstractField, FunctionField

import Oceananigans.BoundaryConditions: fill_halo_regions!
import Oceananigans.Models.IncompressibleModels: extract_boundary_conditions

using Adapt

struct PrescribedVelocityFields{U, V, W, P}
    u :: U
    v :: V
    w :: W
    parameters :: P
end

@inline Base.getindex(U::PrescribedVelocityFields, i) = getindex((u=U.u, v=U.v, w=U.w), i)

zerofunc(args...) = 0

"""
    PrescribedVelocityFields(; u=zerofunc, v=zerofunc, w=zerofunc, parameters=nothing)

Builds `PrescribedVelocityFields` with prescribed functions `u`, `v`, and `w`.

If `isnothing(parameters)`, then `u, v, w` are called with the signature

    `u(x, y, z, t) = # something interesting`

If `!isnothing(parameters)`, then `u, v, w` are called with the signature

    `u(x, y, z, t, parameters) = # something parameterized and interesting`

In the constructor for `HydrostaticFreeSurfaceModel`, the functions `u, v, w` are wrapped
in `FunctionField` and associated with the model's `grid` and `clock`.
"""
PrescribedVelocityFields(; u=zerofunc, v=zerofunc, w=zerofunc, parameters=nothing) =
    PrescribedVelocityFields(u, v, w, parameters)

PrescribedField(X, Y, Z, f::Function,      grid; kwargs...) = FunctionField{X, Y, Z}(f, grid; kwargs...)
PrescribedField(X, Y, Z, f::AbstractField, grid; kwargs...) = f

function HydrostaticFreeSurfaceVelocityFields(velocities::PrescribedVelocityFields, arch, grid, clock, bcs)

    u = PrescribedField(Face, Center, Center, velocities.u, grid; clock=clock, parameters=velocities.parameters)
    v = PrescribedField(Center, Face, Center, velocities.v, grid; clock=clock, parameters=velocities.parameters)
    w = PrescribedField(Center, Center, Face, velocities.w, grid; clock=clock, parameters=velocities.parameters)

    return PrescribedVelocityFields(u, v, w, velocities.parameters)
end

@inline fill_halo_regions!(::PrescribedVelocityFields, args...) = nothing
@inline fill_halo_regions!(::FunctionField, args...) = nothing

ab2_step_velocities!(::PrescribedVelocityFields, args...) = [NoneEvent()]
ab2_step_free_surface!(::Nothing, args...) = NoneEvent()
compute_w_from_continuity!(::PrescribedVelocityFields, args...) = nothing

validate_velocity_boundary_conditions(::PrescribedVelocityFields) = nothing
extract_boundary_conditions(::PrescribedVelocityFields) = NamedTuple()

FreeSurfaceDisplacementField(::PrescribedVelocityFields, arch, grid) = nothing
HorizontalVelocityFields(::PrescribedVelocityFields, arch, grid) = nothing, nothing
FreeSurface(free_surface::ExplicitFreeSurface{Nothing}, ::PrescribedVelocityFields, arch, grid) = nothing

hydrostatic_prognostic_fields(::PrescribedVelocityFields, free_surface, tracers) = tracers
calculate_hydrostatic_momentum_tendencies!(tendencies, ::PrescribedVelocityFields, args...) = []

apply_flux_bcs!(::Nothing, c, arch, events, barrier, clock, model_fields) = nothing

Adapt.adapt_structure(to, velocities::PrescribedVelocityFields) =
    PrescribedVelocityFields(Adapt.adapt(to, velocities.u),
                             Adapt.adapt(to, velocities.v),
                             Adapt.adapt(to, velocities.w),
                             nothing)


text:
module chapters_interactive_graphing

using Dash, DashHtmlComponents, DashCoreComponents

export examples

include("../../utils.jl")

include("../../reusable_components/Header.jl")
include("../../reusable_components/Example.jl")
include("../../reusable_components/Syntax.jl")

crossfilter_recipe = LoadExampleCode("./dash_docs/chapters/graph_crossfiltering/examples/crossfilter-recipe.jl")

interactive_graph_1 = LoadExampleCode("./dash_docs/chapters/graph_crossfiltering/examples/interactive-graph-1.jl")

interactive_graph_2 = LoadExampleCode("./dash_docs/chapters/graph_crossfiltering/examples/interactive-graph-2.jl")

examples = [crossfilter_recipe, interactive_graph_1, interactive_graph_2]

app =  dash()

crossfilter_recipe.callback!(app)
interactive_graph_1.callback!(app)
interactive_graph_2.callback!(app)

app.layout = html_div() do
    Header("Interactive Visualizations"),
    html_blockquote(dcc_markdown("This is the 4th chapter of the [Dash Tutorial](/).
    The [previous chapter](/basic-callbacks) covered basic callback usage and the [next chapter](/state)
    describes how to share data between callbacks. Just getting started? Make sure to [install the necessary
    dependencies](/installation)")),
    dcc_markdown("
    The `DashCoreComponents` package includes a component called `Graph`.exec

    `Graph` renders interactive visualizations using the open source [plotly.js](https://github.com/plotly/plotly.js) JavaScript
    graphing library. Plotly.js supports over 35 chart types and renders charts in both vector-quality SVG and high-performance
    WebGL.

    The `figure` argument in the `dcc_graph` component is the same `figure` argument that is used by
    `plotlyjs.jl`, an open source Julia graphing library. Check out
    the [plotlyjs.jl documentation](http://juliaplots.org/PlotlyJS.jl/stable/) and gallery to learn more.

    Dash components are described declaratively by a set of attributes. All of these attributes can be updated by callback
    functions, but only a subset of these attributes are updated through user interaction, such as when you click on an
    option in a `Dropdown` component and the `value` property of that component changes.

    The `Graph` component has four attributes that can change through user interaction: `hoverData`, `clickData`, `selectedData`,
    and `relayoutData`. THese properties update when you hover over points, click on points, or select regions in a graph.

    Here's a simple example that prints these attributes to the screen.
    "),
    interactive_graph_1.source_code,
    Example(interactive_graph_1.layout),
    Header("Update Graphs on Hover"),
    dcc_markdown("""
    Let's update our world indicators example from the previous chapter by updating time series
    when we hover over points in our scatter plot.
    """),

    interactive_graph_2.source_code,
    Example(interactive_graph_2.layout),
    dcc_markdown("""

    Try mousing over the points in the scatter plot on the left.
    Notice how the line graphs on the right update based off of the point that you are hovering over."""),

    Header("Generic Crossfilter Recipe"),

    dcc_markdown("""
    Here's a slightly more generic example for crossfiltering across a six-column dataset.
    Each scatter plot's selection filters the underlying dataset.

    """),
    crossfilter_recipe.source_code,
    Example(crossfilter_recipe.layout),
    dcc_markdown("""

    Try clicking and dragging in any of the plots to filter different regions.
    On every selection, the three graph callbacks are fired with the latest selected regions of each plot.
    A dataframe is filtered based off of the selected points and the graphs are replotted
    with the selected points highlighted and the selected region drawn as a dashed rectangle.

    """),

    Header("Current Limitations"),

    dcc_markdown("""

    There are a few limitations in graph interactions right now.

    * It is not currently possible to customize the style of the hover interactions
    or the select box. This issue is being worked on in https://github.com/plotly/plotly.js/issues/1847.

    There's a lot that you can do  with these interactive plotting features. If you need
    help exploring your use case, open up a thread in the [Dash Community Forum](https://community.plotly.com/c/dash/julia/20)

    The next chapter of the Dash User Guide explains how to share data between
    callbacks.

    [Dash Tutorial Part 5: Sharing Data Between Callbacks](/sharing-data-between-callbacks)
    """),

    html_hr(),
    dcc_markdown("[Back to the table of contents](/)")

end

end

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

compute(n::Int)::Int = (((4n + 3)n + 8)n - 9) ÷ 6

compute(5)

compute(1_001)

@benchmark compute(1_001)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
@testset "infill" begin
    using Statistics

    #KD-tree for the plan as well as the samples from the plan
    search_range = [(-10.0,10.0), (0.0,2.0), (25.0,40.0)]
    n_samples = 6
    plan = SurrogateModelOptim.scaled_LHC_sampling_plan(search_range,n_samples,1000;trace=false)
    samples = mapslices(sum,plan,dims=1)

    combined_points_samples = [plan;samples]
    kdtree = SurrogateModelOptim.KDTree(combined_points_samples)

    #Infill objectives
    @test SurrogateModelOptim._nearest_point(kdtree,x->sum(x),plan[:,1]) == 0
    @test SurrogateModelOptim.distance_infill(plan,samples,x->sum(x))(plan[:,1]) == 0
    @test SurrogateModelOptim.minimum_infill(x->(x, 2x, 4x))(1) == 1
    @test SurrogateModelOptim.median_infill(x->(x, 2x, 4x))(1) == 2
    @test SurrogateModelOptim.mean_infill(x->(x, 2x, 4x))(1) == mean([1,2,4])
    @test SurrogateModelOptim.std_infill(x->(x, 2x, 4x))(1) == -1*std([1 2 4])

    #Test that try-catch block correctly shows exception when enough iterations have passed
    @test_throws ArgumentError SurrogateModelOptim.infill_opt(search_range,2,2,sum,[:std],plan,sum,SurrogateModelOptim.Options())

    #Make sure correct input format is used
    @test_throws ErrorException model_infill(search_range,permutedims(plan),samples,sum;
            options = SurrogateModelOptim.Options(trace=false))

    #Make sure that input format is used
    @test_throws ErrorException model_infill(search_range,permutedims(plan),samples,
            sum; options = SurrogateModelOptim.Options( infill_funcs = [:std,:median],
                                                        num_interpolants = 1,
                                                        trace = false))
end


text:
"""
    Saturation function for quadratic saturation models for machines
        Se(x) = B * (x - A)^2 / x
"""
function saturation_function(
    machine::Union{PSY.RoundRotorQuadratic, PSY.SalientPoleQuadratic},
    x::Real,
)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(machine)
    return Sat_B * (x - Sat_A)^2 / x
end

"""
    Saturation function for exponential saturation models for machines
        Se(x) = B * x^A
"""
function saturation_function(
    machine::Union{PSY.RoundRotorExponential, PSY.SalientPoleExponential},
    x::Real,
)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(machine)
    return Sat_B * x^Sat_A
end

function rectifier_function(I::Float64)
    if I <= 0.0
        return 1.0
    elseif I <= 0.433
        return 1.0 - 0.577 * I
    elseif I < 0.75
        return sqrt(0.75 - I^2)
    elseif I <= 1.0
        return 1.732 * (1.0 - I)
    else
        return 0.0
    end
end

function saturation_function(avr::PSY.ESAC1A, x::Real)
    Sat_A, Sat_B = PSY.get_saturation_coeffs(avr)
    return Sat_B * (x - Sat_A)^2 / x
end

function rectifier_function(I::T) where {T <: Real}
    if I <= 0.0
        return one(T)
    elseif I <= 0.433
        return 1.0 - 0.577 * I
    elseif I < 0.75
        return sqrt(0.75 - I^2)
    elseif I <= 1.0
        return 1.732 * (1.0 - I)
    else
        return zero(T)
    end
end

function output_pss_limiter(
    V_ss::X,
    V_ct::Y,
    V_cl::Float64,
    V_cu::Float64,
) where {X <: Real, Y <: Real}
    # Bypass limiter block if one limiter parameter is set to zero.
    if V_cl == 0.0 || V_cu == 0.0
        return V_ss
    end
    if V_cl <= V_ct <= V_cl
        return V_ss
    elseif V_ct < V_cl
        return zero(X)
    elseif V_ct > V_cu
        return zero(X)
    else
        error("Logic error in PSS")
    end
end

function deadband_function(x::T, db_low::Float64, db_high::Float64) where {T <: Real}
    if x > db_high
        return x - db_high
    elseif x < db_low
        return x - db_low
    else
        return zero(T)
    end
end

function current_limit_logic(
    inner_control::PSY.RECurrentControlB,
    ::Type{Base.RefValue{0}}, #PQ_Flag = 0: Q Priority
    Vt_filt::X,
    Ip_cmd::Y,
    Iq_cmd::Z,
) where {X, Y, Z <: Real}
    # X,Y,Z should be the same. Future work: Correct bug to work always with Float64 or Dual.
    I_max = PSY.get_I_max(inner_control)
    Iq_max = I_max
    Iq_min = -Iq_max
    Ip_min = 0.0
    local_I = I_max^2 - Iq_cmd^2
    if local_I < 0
        local_I = 0
    else
        local_I = sqrt(local_I)
    end
    if local_I < I_max
        Ip_max = local_I
    else
        Ip_max = I_max
    end
    return Ip_min, Ip_max, Iq_min, Iq_max
end

function current_limit_logic(
    inner_control::PSY.RECurrentControlB,
    ::Type{Base.RefValue{1}}, #PQ_Flag = 1: P Priority
    Vt_filt::X,
    Ip_cmd::Y,
    Iq_cmd::Z,
) where {X, Y, Z <: Real}
    # X,Y,Z should be the same. Future work: Correct bug to work always with Float64 or Dual.
    I_max = PSY.get_I_max(inner_control)
    Ip_max = I_max
    Ip_min = 0.0
    local_I = I_max^2 - Ip_cmd^2
    if local_I < 0
        local_I = 0
    else
        local_I = sqrt(local_I)
    end
    if local_I < Iq_max
        Iq_max = local_I
    else
        Iq_max = I_max
    end
    Iq_min = -Iq_max
    return Ip_min, Ip_max, Iq_min, Iq_max
end

function get_LVPL_gain(
    Vmeas::T,
    Zerox::Float64,
    Brkpt::Float64,
    Lvpl1::Float64,
) where {T <: Real}
    if Vmeas < Zerox
        return zero(T)
    elseif Vmeas > Brkpt
        return one(T) * Lvpl1
    else
        return Lvpl1 / (Brkpt - Zerox) * (Vmeas - Zerox)
    end
end

function get_LV_current_gain(V_t::T, Lv_pnt0::Float64, Lv_pnt1::Float64) where {T <: Real}
    if V_t < Lv_pnt0
        return zero(T)
    elseif V_t > Lv_pnt1
        return one(T)
    else
        return 1.0 / (Lv_pnt0 - Lv_pnt1) * (V_t - Lv_pnt0)
    end
end


text:
import GLFW
using ModernGL, GeometryTypes
# We're going to use only one feature of GLAbstraction---a slightly
# more "julian" version of glShaderSource. There are many other places
# where it offers simplifications, but in this tutorial we're being
# deliberately low-level.
using GLAbstraction

# A more comprehensive setting of GLFW window hints. Setting all
# window hints reduces platform variance.
# In future files, we'll make use of GLWindow to handle this automatically.
window_hint = [
    (GLFW.SAMPLES,      4),
    (GLFW.DEPTH_BITS,   0),

    (GLFW.ALPHA_BITS,   8),
    (GLFW.RED_BITS,     8),
    (GLFW.GREEN_BITS,   8),
    (GLFW.BLUE_BITS,    8),
    (GLFW.STENCIL_BITS, 0),
    (GLFW.AUX_BUFFERS,  0),
    (GLFW.CONTEXT_VERSION_MAJOR, 3),# minimum OpenGL v. 3
    (GLFW.CONTEXT_VERSION_MINOR, 0),# minimum OpenGL v. 3.0
    (GLFW.OPENGL_PROFILE, GLFW.OPENGL_ANY_PROFILE),
    (GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE),
]

for (key, value) in window_hint
    GLFW.WindowHint(key, value)
end

# Create the window
window = GLFW.CreateWindow(800, 600, "Drawing polygons 1")
GLFW.MakeContextCurrent(window)
# Retain keypress events
GLFW.SetInputMode(window, GLFW.STICKY_KEYS, GL_TRUE)

# Create the Vertex Array Object (VAO) and make it current
# Note that while the tutorial describes this after the attributes (below),
# we need to make vao current before calling glVertexAttribPointer.
# You should also do this before creating any element arrays (see
# drawing_polygons4.jl)
vao = Ref(GLuint(0))
glGenVertexArrays(1, vao)
glBindVertexArray(vao[])

# The vertices of our triangle
vertices = Point2f0[(0, 0.5), (0.5, -0.5), (-0.5, -0.5)] # note Float32

# Create the Vertex Buffer Object (VBO)
vbo = Ref(GLuint(0))   # initial value is irrelevant, just allocate space
glGenBuffers(1, vbo)
glBindBuffer(GL_ARRAY_BUFFER, vbo[])
glBufferData(GL_ARRAY_BUFFER, sizeof(vertices), vertices, GL_STATIC_DRAW)

# The shaders. Here we do everything manually, but life will get
# easier with GLAbstraction. See drawing_polygons5.jl for such an
# implementation.

# The vertex shader
vertex_source = """
#version 150

in vec2 position;

void main()
{
    gl_Position = vec4(position, 0.0, 1.0);
}
"""

# The fragment shader
fragment_source = """
# version 150

out vec4 outColor;

void main()
{
    outColor = vec4(1.0, 1.0, 1.0, 1.0);
}
"""

# Compile the vertex shader
vertex_shader = glCreateShader(GL_VERTEX_SHADER)
glShaderSource(vertex_shader, vertex_source)  # nicer thanks to GLAbstraction
glCompileShader(vertex_shader)
# Check that it compiled correctly
status = Ref(GLint(0))
glGetShaderiv(vertex_shader, GL_COMPILE_STATUS, status)
if status[] != GL_TRUE
    buffer = Array(UInt8, 512)
    glGetShaderInfoLog(vertex_shader, 512, C_NULL, buffer)
    error(bytestring(buffer))
end

# Compile the fragment shader
fragment_shader = glCreateShader(GL_FRAGMENT_SHADER)
glShaderSource(fragment_shader, fragment_source)
glCompileShader(fragment_shader)
# Check that it compiled correctly
status = Ref(GLint(0))
glGetShaderiv(fragment_shader, GL_COMPILE_STATUS, status)
if status[] != GL_TRUE
    buffer = Array(UInt8, 512)
    glGetShaderInfoLog(fragment_shader, 512, C_NULL, buffer)
    error(bytestring(buffer))
end

# Connect the shaders by combining them into a program
shader_program = glCreateProgram()
glAttachShader(shader_program, vertex_shader)
glAttachShader(shader_program, fragment_shader)
glBindFragDataLocation(shader_program, 0, "outColor") # optional

glLinkProgram(shader_program)
glUseProgram(shader_program)

# Link vertex data to attributes
pos_attribute = glGetAttribLocation(shader_program, "position")
glVertexAttribPointer(pos_attribute, length(eltype(vertices)),
                      GL_FLOAT, GL_FALSE, 0, C_NULL)
glEnableVertexAttribArray(pos_attribute)

# Draw while waiting for a close event
glClearColor(0,0,0,0)
while !GLFW.WindowShouldClose(window)
    glClear(GL_COLOR_BUFFER_BIT)
    glDrawArrays(GL_TRIANGLES, 0, length(vertices))
    GLFW.SwapBuffers(window)
    GLFW.PollEvents()
    if GLFW.GetKey(window, GLFW.KEY_ESCAPE) == GLFW.PRESS
        GLFW.SetWindowShouldClose(window, true)
    end
end

GLFW.DestroyWindow(window)  # needed if you're running this from the REPL


text:
using SimpleSMC
using ProgressMeter

function _ccXpf!(model::SMCModel, ccsmcio::CCSMCIO{Particle}, lM::F,
  rngCouple::Bool) where {Particle, F<:Function}
  zetas1 = ccsmcio.smcio1.zetas
  zetaAncs1 = ccsmcio.smcio1.internal.zetaAncs
  lws1 = ccsmcio.smcio1.internal.lws
  ws1 = ccsmcio.smcio1.ws
  as1 = ccsmcio.smcio1.internal.as
  pScratch = ccsmcio.smcio1.internal.particleScratch
  logZhats1 = ccsmcio.smcio1.logZhats
  lZ1 = 0.0

  zetas2 = ccsmcio.smcio2.zetas
  zetaAncs2 = ccsmcio.smcio2.internal.zetaAncs
  lws2 = ccsmcio.smcio2.internal.lws
  ws2 = ccsmcio.smcio2.ws
  as2 = ccsmcio.smcio2.internal.as
  logZhats2 = ccsmcio.smcio2.logZhats
  lZ2 = 0.0

  for p = 1:ccsmcio.n
    if p > 1
      _copyParticles!(zetaAncs1, zetas1, as1)
      _copyParticles!(zetaAncs2, zetas2, as2)
    end

    if rngCouple
      _rngCoupledMutateParticles!(zetas1, zetas2, model.M!, p, zetaAncs1,
        zetaAncs2, pScratch, ccsmcio.ref1[p], ccsmcio.ref2[p], ccsmcio.extraRNG)
    else
      _indexCoupledMutateParticles!(zetas1, zetas2, model.M!, p, zetaAncs1,
        zetaAncs2, pScratch, ccsmcio.ref1[p], ccsmcio.ref2[p])
    end

    _logWeightParticles!(lws1, p, model.lG, zetas1, pScratch)
    _logWeightParticles!(lws2, p, model.lG, zetas2, pScratch)

    maxlw::Float64 = maximum(lws1)
    ws1 .= exp.(lws1 .- maxlw)
    mws::Float64 = mean(ws1)
    lZ1 += maxlw + log(mws)
    ws1 ./= mws
    @inbounds logZhats1[p] = lZ1

    maxlw = maximum(lws2)
    ws2 .= exp.(lws2 .- maxlw)
    mws = mean(ws2)
    lZ2 += maxlw + log(mws)
    ws2 ./= mws
    @inbounds logZhats2[p] = lZ2

    if p < ccsmcio.n
      _coupledResample!(ccsmcio)
      lM != error && _ancestorSample!(ccsmcio, p, lM)
    end

    _intermediateOutput!(ccsmcio.smcio1, p)
    _intermediateOutput!(ccsmcio.smcio2, p)
  end
end

function ccXpf!(model::SMCModel, lM::F, ccsmcio::CCSMCIO{Particle},
  algorithm::Symbol = :BS, rngCouple::Bool = false) where {F<:Function, Particle}
  if algorithm == :AS
    _ccXpf!(model, ccsmcio, lM, rngCouple)
    _pickParticles!(ccsmcio)
  else
    _ccXpf!(model, ccsmcio, error, rngCouple)
    if algorithm == :BS
      _pickParticlesBS!(ccsmcio, lM)
    else
      _pickParticles!(ccsmcio)
    end
  end
end

function ccXpf!(model::SMCModel, ccsmcio::CCSMCIO{Particle},
  rngCouple::Bool = false) where Particle
  ccXpf!(model, error, ccsmcio, :AT, rngCouple)
end

function initializeCCSMC(model::SMCModel, lM::F, ccsmcio::CCSMCIO{Particle},
  algorithm::Symbol = :BS, independent::Bool = false) where
  {F<:Function, Particle}
  ref1 = ccsmcio.ref1
  ref2 = ccsmcio.ref2

  smcio1 = ccsmcio.smcio1
  smcio2 = ccsmcio.smcio2

  SimpleSMC.smc!(model, smcio1)
  SimpleSMC.pickParticle!(ref1, smcio1)

  if independent
    SimpleSMC.smc!(model, smcio2)
    SimpleSMC.pickParticle!(ref2, smcio2)
  else
    SimpleSMC._copyParticles!(ref2, ref1)
  end
  if algorithm == :BS || algorithm == :AS
    SimpleSMC.csmc!(model, lM, smcio1, ref1, ref1)
  else
    SimpleSMC.csmc!(model, smcio1, ref1, ref1)
  end
end

function initializeCCSMC(model::SMCModel, ccsmcio::CCSMCIO{Particle},
  independent::Bool = false) where Particle
  initializeCCSMC(model, error, ccsmcio, :AT, independent)
end


text:
"""
    RegularRectilinearCollocatedGrid{FT, TX, TY, TZ, R} <: AbstractRectilinearGrid{FT, TX, TY, TZ}

A rectilinear grid with collocated points and constant grid spacings `Δx`, `Δy`, and `Δz` between cell centers
and cell faces, elements of type `FT`, topology `{TX, TY, TZ}`, and coordinate ranges
of type `R`.
"""
# NOTE: we can later add info here about the local & global extent of the grid, the halo, the neigboring processors and the type of calculation (CPU/MPI-CPU/GPU)
mutable struct RegularRectilinearCollocatedGrid{FT, D, B} <: AbstractRectilinearCollocatedGrid{FT, D, B}
    # number of dimensions
    dim :: Int

    # Number of global grid points in (x,y,z) [is specified by user]
    N  :: NTuple{D,Int}

    # Actual number of global grid points in (x,y,z) (includes halo if using IGG)
    Ng :: NTuple{D,Int}

    # Local number of local grid points in (x,y,z) on everty processor (w/oout halo)
    Nl :: NTuple{D,Int}

    # Topology
    topology :: NTuple{D,DataType}

    # Stencilwidth
    stencilwidth :: Int

    # Domain size
    L  :: NTuple{D,FT}

    # Grid spacing
    Δ  :: NTuple{D,FT}

    # Range of (local) coordinates at the centers of the cells.
    Center  :: NTuple{D,StepRangeLen}

    # Range of (local) coordinates  (+halo in case of IGG) at the faces of the cells.
    Face  :: NTuple{D,StepRangeLen}

    # Range of (global) coordinates at the centers of the cells.
    Center_global  :: NTuple{D,StepRangeLen}

    # Range of (global) coordinates at the faces of the cells.
    Face_global  :: NTuple{D,StepRangeLen}

    # Indices of local portion of grid in x/y/z direction (w/out halo)
    ind_local :: NTuple{D,UnitRange{Int}}

    # Function which computes local indices from CartesianIndex
    idx

    # Corners not including ghost points
    corners :: NamedTuple

    # Corners with ghost points
    ghostcorners :: NamedTuple

    # Backend employed
    backend :: B

    # store PETSc specific info (if used)
    PETSc   :: petsc_data

    # store data in case we employ ImplicitGlobalGrid
    IGG

    # fields
    fields

end

"""
        RegularRectilinearCollocatedGrid(;
                            size=nothing,
                            local_size=nothing,
                            x = nothing, y = nothing, z = nothing,
                            topology= (Bounded, Bounded, Bounded),
                             extent = nothing,
                             fields = NamedTuple(),
                        stencilwidth= 1,
                         stenciltype= :Star,
                                opts= (),
                        )

Creates a `RegularRectilinearCollocatedGrid` with `size = (Nx, Ny, Nz)` grid cells.
All variables are collocated, that is, placed at the same points

Keyword arguments
=================

- `size` : A tuple prescribing the number of grid points all directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and a 1-tuple for 1D.
- `local_size` : A tuple prescribing the number of grid points all directions.
                     `size` is a 3-tuple for 3D models, a 2-tuple for 2D models, and a 1-tuple for 1D.

- `extent`: A tuple prescribing the physical extent of the grid.
            The origin for three-dimensional domains is the lower left corner `(-Lx/2, -Ly/2, -Lz)`.

- `x`, `y`, and `z`: Each of `x, y, z` are 2-tuples that specify the end points of the domain
                     in their respect directions.

- `topology`: A tuple prescribing the topology of the grid in the different directions.
              This is also used to set the type of the `DMDA` boundary conditions in case of using the `PETSc` backend.
              Examples are `Bounded`, `Ghost`, `Periodic`

- `fields`: A NamedTuple with the names of fields that are initialized using the correct backend usinng a specifiable, constant, value (e.g., `fields=(T=0, P=11.1)`)

*Note*:
- _Either_ `extent`, or all of `x`, `y`, and `z` must be specified.
- _Either_ `size`, or `local_size` (the dimensions per processor; preferred by IGG) must be specified.

Optional PETSc arguments
========================

- `stencilwidth`: the width of the stencil

- `stenciltype`: the type of the stencil (`:Star` or `:Box`)

- `opts`: additional PETSc options

Optional ImplicitGlobalGrid arguments
=====================================

- `stencilwidth`=1: the width of the stencil (`overlapx`=`overlapy`=`overlapz`= `stencilwidth+1`)

- `opts::Dict`: a dictionary that contains additional options that can be passed to IGG, such as `opts=Dict(:dimy=>2,:dimx=>1)`

The physical extent of the domain can be specified via `x`, `y`, and `z` keyword arguments
indicating the left and right endpoints of each dimensions, e.g. `x=(0, 10)` or via
the `extent` argument, e.g. `extent=(Lx, Ly, Lz)` which specifies the extent of each dimension
in which case -Lx/2 ≤ x ≤ Lx/2, -Ly/2 ≤ y ≤ Ly/2, and -Lz ≤ z ≤ 0.

Constants are stored using floating point values of type `FT`, which is the type specified in the backend structure

Grid properties
===============

- `(Nx, Ny, Nz)::Int`: Number of physical points in the (x, y, z)-direction

- `(Lx, Ly, Lz)::FT`: Physical extent of the grid in the (x, y, z)-direction

- `(Δx, Δy, Δz)::FT`: Center width in the (x, y, z)-direction

- `(xC, yC, zC)`: (x, y, z) coordinates of cell centers.

- `(xF, yF, zF)`: (x, y, z) coordinates of cell faces.

Examples
========

* A default grid with Float64 type:

```jldoctest
julia> using CompGrids, ParallelStencil
julia> @init_backend(ParallelStencil, Threads, false, Float64);
julia> grid = RegularRectilinearCollocatedGrid(size=(32, 32, 32), extent=(1, 2, 3))
RegularRectilinearCollocatedGrid{Float64, 3, Backend{BackendParallelStencil}}
        Backend: ParallelStencil ( Threads )
       gridtype: Collocated
         domain: x ∈ [-0.5, 0.5], y ∈ [-1.0, 1.0], z ∈ [-3.0, 0.0]
       topology: (Bounded, Bounded, Bounded)
     resolution: (32, 32, 32)
 grid spacing Δ: (0.03125, 0.0625, 0.09375)
```

"""
function RegularRectilinearCollocatedGrid(;
                                    size=(),
                                    local_size=(),
                                     x = nothing, y = nothing, z = nothing,
                               topology= (Bounded, Bounded, Bounded),
                                extent = nothing,
                                fields = NamedTuple(),
                           stencilwidth= 1,
                            stenciltype= :Star,
                                   opts= ()
                              )
    FT = backend.Scalar;        # Scalar type

    # Checks
    if !isempty(local_size) & !isempty(size)
        error("You need to specify either size or local_size, but you cannot specify both at the same time!")
    end
    if isempty(local_size) & isempty(size)
        error("You need to specify either size or local_size")
    end

    # Create Tuples if needed
    if typeof(size)==Int64
        size = (size,)
    end
    if typeof(extent)==FT
        extent = (extent,)
    end

    size = globalfromlocalsize(size, local_size, opts, stencilwidth, topology, backend)   # Compute global size from local size
    dim =   length(size)                                  # dimensions of the grid [1-3]
    L, X₁ = validate_regular_grid_domain(FT, extent, x, y, z)

    # Unpacking
    N = size
    Δ = L ./ (N .- FT(1))

    # Face-node limits in x, y, z
    XF₋ = @. X₁
    XF₊ = @. XF₋ + L

    # Center-node limits in x, y, z
    XC₋ = @. XF₋ + Δ / 2
    XC₊ = @. XC₋ + L - Δ

    # Generate 1D coordinate arrays of faces & centers in all directionbs
    Face=Center=()
    for idim=1:dim
        Face   = (Face...,   range(XF₋[idim], XF₊[idim]; length = N[idim]  ))
        Center = (Center..., range(XC₋[idim], XC₊[idim]; length = N[idim]-1))
    end

    # Local coordinate index arrays (updated for parallel simulations)
    ind_local = ()
    for idim=1:dim
        ind_local = (ind_local..., 1:N[idim])
    end

    # Initialize corners (to iterate over local portion of grid; note that these values are overwritten in parallel cases)
    corners = (lower = CartesianIndex((ones(Int64,dim)...,)), upper=CartesianIndex(N), size=N );
    ghostcorners = (lower = CartesianIndex((ones(Int64,dim)...,)), upper=CartesianIndex(N), size=N );

    # Initialise backend and create grid structure
    grid = RegularRectilinearCollocatedGrid{FT, dim, typeof(backend)}(
          dim, N, N, N,                                 # dimensions, and sizes of grid
          topology[1:dim],                              # boundary conditions
          stencilwidth,                                 # size of halo (in parallel)
          L, Δ,                                         # domain size and (regular) spacing
          Center, Face, Center, Face,                   # data related to the 1D local/global coordinate grids
          ind_local, nothing,                           # indices of local grid
          corners, ghostcorners,                        # corners of local grid
          backend, petsc_data(), nothing,               # data related to backend
          NamedTuple())                                 # fields

    # Initialize grid - this will reset some of the default parameters above
    initialize_grid!(grid,
                    stencilwidth= stencilwidth,
                    stenciltype = stenciltype,
                            opts= opts,
                          fields= fields)

    return grid
end

function validate_regular_grid_domain(FT, extent, x, y, z)

    # Find domain endpoints or domain extent, depending on user input:
    if !isnothing(extent) # the user has specified an extent
        dim   = length(extent)
        x,y,z = nothing, nothing, nothing

        L = extent;

        # Default domain:
        X₁ = -1 .* collect(L) ./ 2
        if dim>2
            X₁[3] = -L[3]
        end
        X₁ = (X₁...,)

    else # isnothing(extent) === true implies that user has not specified a length

        if isnothing(y)
            L = (x[2] - x[1],)
            X₁= (x[1], )
        elseif isnothing(z)
            L = (x[2] - x[1], y[2] - y[1])
            X₁= (x[1], y[1])
        else
            L = (x[2] - x[1], y[2] - y[1], z[2] - z[1])
            X₁= (x[1], y[1], z[1])
        end

    end

    return FT.(L), FT.(X₁)
end

function domain_string(grid)

    xₗ, xᵣ = grid.Face_global[1][1], grid.Face_global[1][end]
    if grid.dim>1
        yₗ, yᵣ = grid.Face_global[2][1], grid.Face_global[2][end]
    end
    if grid.dim>2
        zₗ, zᵣ = grid.Face_global[3][1], grid.Face_global[3][end]
    end
    if grid.dim==1
        return "x ∈ [$xₗ, $xᵣ]"
    elseif grid.dim==2
        return "x ∈ [$xₗ, $xᵣ], y ∈ [$yₗ, $yᵣ]"
    elseif grid.dim==3
        return "x ∈ [$xₗ, $xᵣ], y ∈ [$yₗ, $yᵣ], z ∈ [$zₗ, $zᵣ]"
    end
end

function resolution_string(g)
    if g.backend.mpi
        if g.backend.type==:ParallelStencil
            str = "$(g.N) (global, no halo)
                 $(g.Nl) (local  + halo)
                 $(g.Ng) (global + halo)"
        else
            str = "$(g.N) (global, no halo)
                 $(g.Nl) (local,  no halo)"
        end

    else
        str="$(g.N)"
    end
    return str
end

# view grid object
function show(io::IO, g::RegularRectilinearCollocatedGrid{FT, DIM, B}) where {FT, DIM, B}

    if g.backend.mpi
        mpi_type = "| MPI "
    else
        mpi_type = ""
    end

    print(io, "RegularRectilinearCollocatedGrid{$FT, $DIM, $B}\n",
              "        Backend: $(g.backend.type) ( $(g.backend.arch) $mpi_type) \n",
              "       gridtype: Collocated \n",
              "         domain: $(domain_string(g)) \n",
              "       topology: ", g.topology, '\n',
              "     resolution: $(resolution_string(g))\n",
              " grid spacing Δ: ", g.Δ, "\n",
              "         fields: $(keys(g.fields)) \n" )
end

"""
    initialize_grid(grid::RegularRectilinearCollocatedGrid{FT, 2, backend{BackendPETSc}})

Initializes a `DMDA` object using the `PETSc` backend. Note that this `DMDA` has one degree of freedom,
and is used to store information. If you want to base a solver on this,
"""
function initialize_grid!(grid::RegularRectilinearCollocatedGrid{FT, D, Backend{BackendPETSc,FT}};
        stencilwidth=1, stenciltype=:Star, opts=(), fields::NamedTuple=NamedTuple()) where {FT, D}

    # initialize backend
    petsclib = check_backend(grid.backend, Scalar=backend.Scalar);

    # Transform local boundaries to PETSc boundary conditions
    bcs = Vector{PETSc.LibPETSc.DMBoundaryType}(undef, D)
    for idim=1:D
        bcs[idim] = bcs_translate(grid.topology[idim]);
    end

    # Stencil type
    if stenciltype==:Star
        stenciltype=PETSc.DMDA_STENCIL_STAR;
    else
        stenciltype=PETSc.DMDA_STENCIL_BOX;
    end

    # Create the DMDA
    da = PETSc.DMDA(
        petsclib,
        comm,
        (bcs...,),               # boundary conditions
        grid.N,                  # Global grid size
        1,                       # Number of DOF per node [set to 1 as this DMDA is just to store vectors with data]
        stencilwidth,            # Stencil width
        stenciltype;             # Stencil type
        opts...,
    )

    # Set regular coordinates
    c_start = -1*ones(D)
    c_end = ones(D)
    for idim=1:D
        c_start[idim] = grid.Face[idim][1];
        c_end[idim] = grid.Face[idim][end];
    end
    PETSc.setuniformcoordinates!(da, (c_start...,), (c_end...,))

    # Determine number of local grid points
    corners = PETSc.getcorners(da)
    grid.Nl = corners.size[1:D]         # local with no halo
    grid.Ng = PETSc.getinfo(da).global_size[1:D]     # global size

    # Add local coordinates of Faces (can likely be done more elegantly)
    coord = PETSc.getlocalcoordinatearray(da)   # local coordinate array
    Face_local=Center_local=()
    for iDim=1:D
        Face = coord[iDim,corners.lower]:grid.Δ[iDim]:coord[iDim,corners.upper]
        Center = (Face[2:end] .+ Face[1:end-1])./2.0

        Face_local = (Face_local..., Face)
        Center_local = (Center_local..., Center)
    end
    grid.Face = Face_local
    grid.Center = Center_local

    # initialize fields - they are initialized as global vectors with 1 dof
    names = keys(fields)
    val   = values(fields)
    for ifield=1:length(fields)
        # Create global vector
        v = PETSc.DMGlobalVec(da)
        v .= val[ifield]
        PETSc.assemble!(v)
        grid.fields = add_field(grid.fields, names[ifield], v )

        # Set constant value to global vector
        #fill!(grid.fields[ifield], val[ifield])
    end

    # Define function which computes the local index from a CartesianIndex
    idx_petsc(id::CartesianIndex) = LinearIndices(grid.Nl)[id-grid.corners.lower+oneunit(grid.corners.lower)] + PETSc.ownershiprange(grid.fields[1], false)[1]
    grid.idx = idx_petsc

    # Corners of grid
    grid.corners = PETSc.getcorners(da)
    grid.ghostcorners = PETSc.getghostcorners(da)

    # Store info in a PETSc data object
    grid.PETSc = petsc_data(petsclib, da)

    return nothing
end

"""
    initialize_grid(grid::RegularRectilinearCollocatedGrid{FT, 2, backend{BackendParallelStencil}})

Initializes a grid when we use `ParallelStencil` (and potentially `ImplicitGlobalGrid`) as a Backend
"""
function initialize_grid!(grid::RegularRectilinearCollocatedGrid{FT, D, Backend{BackendParallelStencil,FT}};
        stencilwidth=1, stenciltype=:Star, opts=nothing, fields::NamedTuple=NamedTuple()) where {FT, D}

    # initialize backend
    check_backend(grid.backend);

    # retrieve global & local grid dimensions
    if grid.backend.mpi
        initialize_backend!(grid, grid.backend, stencilwidth, opts, fields)
    end

    # Initialize fields
    initialize_fields!(grid, grid.backend, fields)

    # Define function which computes the local index from a CartesianIndex
    if backend.mpi
        width = 1+grid.stencilwidth
    else
        width = 1
    end
    idx_ps(id::CartesianIndex) = LinearIndices(grid.Nl)[id - grid.corners.lower + oneunit(grid.corners.lower)*width]
    grid.idx = idx_ps

    return nothing
end

"""
    local_coordinates(b::Backend{BackendParallelStencil}, grid)

Computes local coordinate vectors, including the halo points
"""
function local_coordinates!(grid::RegularRectilinearCollocatedGrid{FT, D, Backend{BackendParallelStencil,FT}}, b::Backend{BackendParallelStencil}) where {FT, D}

    Face_local=Center_local=()
    for idim=1:grid.dim
        ind_nohalo = grid.ind_local[idim];

        # As we have a constant grid-spacing this is easy
        ind_s, ind_e = ind_nohalo[1], ind_nohalo[end]
        stencilwidth = grid.stencilwidth

        range_local = -stencilwidth:(ind_e-ind_s)+stencilwidth;

        Face = range_local*grid.Δ[idim] .+ grid.Face[idim][ind_s]
        Center = (Face[2:end] .+ Face[1:end-1])./2.0

        Face_local = (Face_local..., Face)
        Center_local = (Center_local..., Center)
    end

    grid.Face   = Face_local;
    grid.Center = Center_local;

    return nothing

end


text:
using Glob
using Statistics
using LinearAlgebra
include("../src/StochasticBlockModel.jl")

println("Start execution")

println(pwd())
file_num_start = 0

### -------------------------------------------------------------------
all_datasets = glob("instances/S1/in/n8-*.in")[1:10]

out_folder = "experiments/"

ampl_path   = "/opt/ampl.linux64/ampl"

symm_break_const_list = [true]
symm_w_const_list = [true]

### -------------------------------------------------------------------

method = "minlp"

mod_file = "src/minlp.mod"

for (num_file, filepath) in enumerate(all_datasets)
    if num_file < file_num_start
        continue
    end

    dataset_path = split(filepath, "/")[end]
    println("File number: $num_file")
    println(filepath)

    # Load Dataset
    dataset = StochasticBlockModel.Dataset(filepath)
    println(string("Size: ", dataset.n))

    dataset_name = split(basename(filepath),".")[1]

    StochasticBlockModel.MINLP_AMPL(dataset, dataset_name, mod_file, out_folder, ampl_path)

end

println("End of execution.")

max_stars_repo_name: vidalt/Exact-SBM

max_issues_repo_name: vidalt/Exact-SBM

max_forks_repo_name: vidalt/Exact-SBM


text:
# BEGIN PARTICLE TREE
"""
    Coords{B}

Leaf structure in galaxy tree consists of position and momentum coordinates of a specific particle
"""
struct Coords{B}
    values::ArrayPartition{B}

    function Coords(values::ArrayPartition{B}) where B <: Number
        @assert length(values.x[1]) == 3 "Position must have length == 3"
        @assert length(values.x[2]) == 3 "Momentum must have length == 3"
        new{B}(values)
    end

end # Coords

"""
    Species{T<:AbstractMultiScaleArray,B<:Number}

Internally consists of a particle ID and an array of particle coordinates
"""
struct Species{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
    id::ParticleID
end # Species

"""
    Snapshot{T<:AbstractMultiScaleArray,B<:Number}

Internally consists of a collection of all species at one particular time
"""
struct Snapshot{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

"""
    Timescale{T<:AbstractMultiScaleArray,B<:Number}

Internally a collection of all Snapshots with the same time stepping
"""
struct Timescale{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArray{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
    times::Vector{Float64}
    dt::Float64

    function Timescale(nodes::Vector{T},values::Vector{B},end_idxs::Vector{Int},times::Vector{Float64},dt::Float64 = -1.0) where {T,B}
        @assert length(nodes) == length(times) "Number of nodes must equal the number of times"
        new{T,B}(nodes,values,end_idxs,times,dt)
    end
end

"""
    ParticleTree{T<:AbstractMultiScaleArray,B<:Number}

Collection of all timescales
"""
struct ParticleTree{T<:AbstractMultiScaleArray,B<:Number} <: AbstractMultiScaleArrayHead{B}
    nodes::Vector{T}
    values::Vector{B}
    end_idxs::Vector{Int}
end

export Coords, Species, Snapshot, Timescale, ParticleTree
# END PARTICLE TREE

###############################################################################
"""
    ParticleTree()

Make an empty ParticleTree
"""
function ParticleTree()
    # Get ParticleTree leaf type
    coord = Coords(ArrayPartition([1.,2.,3.],[4.,5.,6.]))
    emptycoords = typeof(coord)[]

    # stick Coords in a Species
    # sp = construct(Species,emptycoords,Float64[],ParticleID(0,0,0))
    # emptysp = typeof(sp)[]

    sp = construct(Species,Vector([coord]),Float64[],ParticleID(0,0,0))

    # stick Species in a Snapshot
    # snap = construct(Snapshot,emptysp)
    # emptysnap = typeof(snap)[]
    snap = construct(Snapshot,Vector([sp]))

    # stick Snapshot in a Timescale
    time_init = Vector([0.0])
    dt = -1.0
    scale = construct(Timescale,Vector([snap]),Float64[],time_init,dt)
    emptyscale = typeof(scale)[]
    # stick Timescale in a ParticleTree
    particletree = construct(ParticleTree,Vector([scale]))
    remove_node!(particletree,1)
    return particletree
end

# FUNCTIONS TO HELP GENERATE PARTICLE TREES
"""
    ParticleTree(particletype::ParticleID,energy::Quantity,xyz::Vector{T},aim::Vector,dt::Quantity) where {T<:Quantity}

Return a ParticleTree with a single particle in it. Applies unit conversion prior to storing tree.

Input:
- particletype: ParticleID
- energy: initial particle energies (must have units ∈ Unitful.Energy)
- xyz: starting Cartesian coordinates with Earth along the -x axis (must have units ∈ Unitful.Distance)
- aim: Vector with length=3 with initial direction of particle trajectory.
- dt: initial time step for simulation (must have units ∈ Unitful.Time)

Returns:
        ParticleTree
"""
function ParticleTree(particletype::ParticleID,energy::Quantity,xyz::Vector{T},aim::Vector,dt::Quantity) where {T<:Quantity}
    @assert typeof(energy) <: Unitful.Energy "Need to be units of energy"
    @assert typeof(xyz[1]) <: Unitful.Length "Need to be units of length"
    @assert typeof(dt)     <: Unitful.Time   "Need to be units of time"
    @assert length(xyz) == 3                 "Coordinates require 3 spatial dimensions"
    @assert length(aim) == 3                 "Direction requires 3 spatial dimensions"

    energy_in = usistrip(Float64(energy))
    xyz_in = Float64.(xyz)
    xyz_in = usistrip.(xyz_in)
    aim_in = Float64.(aim)
    dt_in = usistrip(Float64(dt))

    ParticleTree(particletype,energy_in,xyz_in,aim_in,dt_in)
end

"""
    ParticleTree(particletype::ParticleID, energy::Float64, xyz::Vector{Float64}, aim::Vector{Float64}, dt::Float64)

Return a ParticleTree with a single particle in it. Units are assumed to be SI.

Input:
- particletype: ParticleID
- energy: initial particle energies
- xyz: starting Cartesian coordinates with Earth along the -x axis
- aim: Vector with length=3 with initial direction of particle trajectory.
- dt: initial time step for simulation

Returns:
        ParticleTree
"""
function ParticleTree(particletype::ParticleID, energy::Float64, xyz::Vector{Float64}, aim::Vector{Float64}, dt::Float64)

    # Get ParticleTree leaf type
    leafType = typeof(Coords(ArrayPartition([1.,2.,3.],[4.,5.,6.])))

    # Ensure that the initial direction is normalized for subsequent
    # momentum scaling.
    particleAim = normalize(aim)

    # starts at time 0.0
    timeInit = Vector([0.0])

    particleTree = nothing

    # first make Coords
    dir = normalize(aim)
    pmag = momentum_si(energy,particletype)
    mom = pmag*dir
    coords = Coords(ArrayPartition(xyz,mom))

    # stick Coords in a Species
    sp = construct(Species,Vector([coords]),Float64[],particletype)

    # stick Species in a Snapshot
    snap = construct(Snapshot,Vector([sp]))

    # stick Snapshot in a Timescale
    time_init = Vector([0.0])
    scale = construct(Timescale,Vector([snap]),Float64[],time_init,dt)

    # stick Timescale in a ParticleTree
    particletree = construct(ParticleTree,Vector([scale]))

    return particletree

end # ParticleTree

function __push!(pt::ParticleTree,c::Union{Timescale,Snapshot,Species,Coords},I::Vector{T}) where {T<:Int}
    @assert length(I) < 5 "Something went wrong"
    test_type = length(I) == 1 ? typeof(pt.nodes[end]) : typeof(pt[I...])
    if typeof(c) == test_type
        if length(I) == 1
            add_node!(pt,c)
        else
            add_node!(pt,c,I[1:end-1]...)
        end
    else
        if length(I) == 1
            push!(I,length(pt.nodes[end].nodes))
        else
            push!(I,length(pt[I...].nodes))
        end
        __push!(pt,c,I)
    end

end

"""
    push!(pt::ParticleTree,val::Union{ParticleTree,Timescale,Snapshot,Species,Coords})

Push contents of val to last container of same type.
"""
function Base.push!(pt::ParticleTree,val::Union{ParticleTree,Timescale,Snapshot,Species,Coords})
    if isa(val,ParticleTree)
        if length(val.nodes) == 0
            return nothing
        elseif length(pt.nodes) == 0
            append!(pt.nodes,val.nodes)
            append!(pt.values,val.values)
            append!(pt.end_idxs,val.end_idxs)
            return nothing
        else
            for ts in val.nodes
                add_node!(pt,ts)
            end
            return nothing
        end
    else
        I = Vector([length(pt.nodes)])
        __push!(pt,val,I)
    end
    return nothing
end

"""
    push(pt::ParticleTree,val::Union{ParticleTree,Timescale,Snapshot,Species,Coords})

No-side-effect version of push!.
"""
function push(pt::ParticleTree,val::Union{ParticleTree,Timescale,Snapshot,Species,Coords})
    if isa(val,ParticleTree)
        if length(val.nodes) == 0
            return deepcopy(pt)
        elseif length(pt.nodes) == 0
            return deepcopy(val)
        else
            out = deepcopy(pt)
            for ts in val.nodes
                add_node!(out,ts)
            end
            return out
        end
    else
        out = deepcopy(pt)
        I = Vector([length(out.nodes)])
        __push!(pt,val,I)
    end
    return nothing
end
export push

"""
    initParticles(particleTypes, energies, rStart, particleAim, gpFrac, bfield)

Return a ParticleTree consisting of all combinations of particles, energies, radii, and initial direction.
Timestep is calculated as a fraction of the gyroperiod at the particle point of insertion.

Input:
- particleTypes: Vector of ParticleID(s)
- energies: Vector of initial particle energies (must have units ∈ Unitful.Energy)
- rStart: starting radius from galactic center (must have units ∈ Unitful.Distance)
- particleAim: Vector with length=3 with initial direction of particle trajectory.
- gpFrac: Float giving the fraction of the gyroperiod to use as a time step
- bfield: Function that returns the value of the galactic magnetic field when given spatial coorinates

Returns:
        ParticleTree
"""
function initParticles(particleTypes, energies, rStart, particleAim, gpFrac, bfield)

    # Get ParticleTree leaf type
    leafType = typeof(Coords(ArrayPartition([1.,2.,3.],[4.,5.,6.])))

    # Ensure that the initial direction is normalized for subsequent
    # momentum scaling.
    particleAim = map(x->normalize(x),particleAim)

    # starts at time 0.0
    timeInit = Vector([0.0])

    particleTree = nothing
    # loop over all arrays to get combination of all parameters
    for r in rStart
        rVec = Vector([usistrip(r),0.0,0.0])
        bMag = norm(bfield(rVec))
        for en in energies
            for p in particleTypes
                pMag = momentum(en,p)
                vMag = vFromP(pMag,mass(p))
                coordsVec = nothing
                for d in particleAim
                    println("Test 1")
                    pVec = usistrip(pMag)*d
                    println(d)
                    println(pVec)
                    println(rVec)
                    coords = Coords(ArrayPartition(rVec,pVec))
                    println(coords)
                    if coordsVec === nothing
                        coordsVec = Vector([coords])
                    else
                        push!(coordsVec,coords)
                    end
                end
                # each species will have a different gyroradius, so each
                # one also gets its own time step -> Timescale
                sp = construct(Species,coordsVec,Float64[],p)
                spVec = Vector([sp])
                snap = construct(Snapshot,spVec)
                snapVec = Vector([snap])
                gCirc = 2*π*gyroradius(
                                       ustrip(u"GeV/c",pMag),
                                       abs(p.p-p.e),
                                       usistrip(bMag)*1.0e4
                                       )
                gPrd = gCirc/usistrip(vMag)
                dt = gPrd*gpFrac
                scale = construct(Timescale,snapVec,Float64[],timeInit,dt)
                if particleTree == nothing
                    particleTree = construct(ParticleTree,Vector([scale]))
                else
                    add_node!(particleTree,scale)
                end
            end
        end
    end
    return particleTree

end # initParticles
export initParticles

# TODO: generate readable lookup tree if inputs are given as dictionaries

# TODO: generate from a (.json?) file

###############################################################################
# FUNCTIONS TO HELP READ PARTICLE TREES

"""
        getCoordArray(pt::ParticleTree)

Returns all (position,momentum) coordinates in an array of size (6,nCoords)
"""
function getCoordArray(pt::ParticleTree)
    out = Vector(pt)
    out = transpose(reshape(pt,(6,length(pt)÷6)))
    return out
end

export getCoordArray

"""
        getTimescaleIndices(pt::ParticleTree)

Returns internal MultiScaleArray indices denoting the last index of each
Timescale.
"""
function getTimescaleIndices(pt::ParticleTree)
    return pt.end_idxs
end

export getTimescaleIndices

"""
        getSnapshotIndices(pt::ParticleTree)

Returns internal MultiScaleArray indices denoting the last index of each
Snapshot.
"""
function getSnapshotIndices(pt::ParticleTree)
    out = nothing
    for snapshot in LevelIter(2,pt)
        if out == nothing
            out = Vector([snapshot[1].end_idxs])
        else
            push!(out,snapshot[1].end_idxs)
        end
    end
    return out
end

export getSnapshotIndices

"""
        getSpeciesIndices(pt::ParticleTree)

Returns internal MultiScaleArray indices denoting the last index of each
Species.
"""
function getSpeciesIndices(pt::ParticleTree)
    out = nothing
    for species in LevelIter(3,pt)
        if out === nothing
            out = Vector([species[1].end_idxs])
        else
            push!(out,species[1].end_idxs)
        end
    end
    return out
end

export getSpeciesIndices

"""
        getParticleTrajectory(pt::ParticleTree,iTs::Int64,iSs::UnitRange,iSp::Int64,iCd::Int64)

Returns particle trajectory over given range iSs
"""
function getParticleTrajectory(pt::ParticleTree,iTs::Int64,iSs::UnitRange,iSp::Int64,iCd::Int64)
    return pt.nodes[iTs].times, Vector([pt[iTs,i,iSp,iCd].values for i = iSs])
end

"""
        getParticleTrajectory(pt::ParticleTree,iTs::Int64,iSs::UnitRange,iSp::Int64,iCd::Int64)

Returns entire particle trajectory
"""
function getParticleTrajectory(pt::ParticleTree,iTs::Int64,iSp::Int64,iCd::Int64)
    nTime = length(pt.nodes[iTs].nodes)
    getParticleTrajectory(pt,iTs,1:nTime,iSp,iCd)
end

export getParticleTrajectory
###############################################################################
# SAVE PARTICLE TREES TO FILE


text:
"""
    eom_logistic3(u, p, t)

Equations of motion for a system consisting of three coupled logistic map
representing the response of two independent dynamical variables to the
forcing from a common driver. The dynamical influence goes in the directions
Z -> X and Z -> Y.

Dynamical noise may be added to each of the dynamical variables by tuning the
parameters `σz`, `σx` and `σz`. Default values for the parameter s
`r₁` and `r₂` are set such that the system exhibits chaotic behaviour.

# References
1. Runge, Jakob. Causal network reconstruction from time series: From
theoretical assumptions to practical estimation, Chaos 28, 075310 (2018);
doi: 10.1063/1.5025050
"""
function eom_logistic3(u, p, t)
    r₁, r₂, r₃, σx, σy, σz = (p...,)
    x, y, z = (u...,)

    # Independent dynamical noise for each variable.
    ηx = rand()
    ηy = rand()
    ηz = rand()

    dx = (x*(r₁ - r₁*x - z + σx*ηx)) % 1
    dy = (y*(r₂ - r₂*y - z + σy*ηy)) % 1
    dz = (z*(r₃ - r₃*z + σz*ηz)) % 1
    return SVector{3}(dx, dy, dz)
end

function logistic3(u₀, r₁, r₂, r₃, σx, σy, σz)
    p = [r₁, r₂, r₃, σx, σy, σz]
    DiscreteDynamicalSystem(eom_logistic3, u₀, p)
end

"""
    logistic3(;u₀ = rand(3), r = 4,
        σx = 0.05, σy = 0.05, σz = 0.05) -> DiscreteDynamicalSystem

Initialise a dynamical system consisting of three coupled logistic map
representing the response of two independent dynamical variables to the
forcing from a common driver. The dynamical influence goes in the directions
Z -> X and Z -> Y.

The equations of motion are

```math
dx = (x(r - r_1 x - z + σ_x η_x)) \\mod 1
dy = (y(r - r_2 y - z + σ_y η_y)) \\mod 1
dz = (z(r - r_3 z + σ_z η_z)) \\mod 1
```

Dynamical noise may be added to each of the dynamical variables by tuning the
parameters `σz`, `σx` and `σz`. Default values for the parameters
`r₁`, `r₂` and `r₃` are set such that the system exhibits chaotic behaviour,
with `r₁ = r₂ = r₃ = 4`.

# References
1. Runge, Jakob. Causal network reconstruction from time series: From theoretical assumptions to practical estimation, Chaos 28, 075310 (2018); doi: 10.1063/1.5025050
"""
logistic3(;u₀ = rand(3), r₁ = 4, r₂ = 4, r₃ = 4,
    σx = 0.05, σy = 0.05, σz = 0.05) = logistic3(u₀, r₁, r₂, r₃, σx, σy, σz)


text:
# Nautilus.jl 🐙 Demo
# Miguel Raz-Guzmán Macedo 2022
# TO RUN THIS DEMO:
# 1. Download this repo and then
# $ cd Nautilus
# $ julia --project=.
# $ julia> include("test/demo.jl")
# 2. Click the toggles.
# 3. Go for coffee because JITs hurt. Seriously.
#
# ###########################################
# VERY ALPHA - feedback welcome, see bottom.#
# (Yes, it's a script with tons of globals, #
#   schade Marmelade 🍓)                    #
# ###########################################
#
# Nautilus' mission:
# "If there exists instrumentation to measure compiler internals,
#       we should have a
#           ♦ visual, interactive ♦
#                   interface for it."
#
# The demo works by having:
# * A vector of LLVM optimization passes.
# * A GPUCompiler workflow around hooks that emit `String`s of compiled `kernel()` closures.
# * A set of Makie `Toggle`s that when on/off change the compilation passes.
# # A Makie plotting window to display the results and offer interactivity.
#
# Here are the basic imports.
using Nautilus
using GLMakie
using GPUCompiler
using LLVM
using LLVM.Interop

# These are the optimization passes in LLVM/src/interop.jl".
# We will turn these on and off with Makie `Toggle`s later on.
global PASSES = [
    alloc_opt!,
    speculative_execution_if_has_branch_divergence!,
    loop_unroll!,
    instruction_combining!,
    licm!,
    early_csemem_ssa!,
    dead_store_elimination!,
    cfgsimplification!,
    global_dce!
]

# We're going to pirate this function now for the purposes of this demo. Be warned 💀!
function GPUCompiler.optimize_module!(@nospecialize(job::CompilerJob{PTXCompilerTarget}),
    mod::LLVM.Module)
    tm = llvm_machine(job.target)
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)

        # This is the new logic:
        # Only activate the passes if the toggles are set to `true`
        # You have to check the `toggle.active.val` field for that.
        for (pass, toggle) in zip(PASSES, TOGGLES)
            if toggle.active.val
                pass(pm)
            end
        end
        run!(pm, mod)
    end
end

# The next few lines of code are GPUCompiler boilerplate.
# They help setup the compilation job with the proper hooks, more on them at another demo.
module TestRuntime
signal_exception() = return
malloc(sz) = C_NULL
report_oom(sz) = return
report_exception(ex) = return
report_exception_name(ex) = return
report_exception_frame(idx, func, file, line) = return
end

struct TestCompilerParams <: AbstractCompilerParams end
GPUCompiler.runtime_module(::CompilerJob{<:Any,TestCompilerParams}) = TestRuntime

# This is the function that we will actually be compiling and recompiling
# upon toggles being flipped on and off.
kernel() = begin
    r = Ref{Int}(42)
    r[] = 43
    nothing
end

function main()
    source = FunctionSpec(kernel)
    target = NativeCompilerTarget()
    params = TestCompilerParams()
    job = CompilerJob(target, source, params)

    # This `menu.selection[]` gets updated whenever the `menu` Observable
    # is updated.
    jobtype = menu.selection[]
    jobtype == :llvm && return GPUCompiler.compile(:llvm, job)[1] |> string
    jobtype == :asm && return GPUCompiler.compile(:asm, job)[1] |> string
end

############################################################################################
# The Actual Nautilus.jl GUI demo begins:                                                  #
# 1. Setup a Figure, a vector of `Toggle`s, a vector of `Label`s to know the `Toggle names.#
# 2. Setup a Menu to choose between `llvmir`/`asm`.                                        #
# 3. Setup a TextBox to paint the compiled code text in.                                   #
# 4. Lift all the `Toggle`s and `Menu` to rerun `main()` upon clicks                       #
# 5. Paint the text in the textbox                                                         #
# ##########################################################################################
#
# 1, 2, 3 Setups
# include("addsomecolor.jl") # this kinda works, we just need to synchronize the Observables!
fig = Figure(resolution = (1800, 2600), backgroundcolor = "#748AA6")
txtsize = Observable(0.08)
toggles = [Toggle(fig) for _ in PASSES]
labels = [Label(fig, lift(x -> x ? "$(string(pass))" : "$(string(pass))", t.active); halign = :left)
          for (t, pass) in zip(toggles, PASSES)]
menu = Menu(fig, options = [:llvm, :asm])
fig[1, 2] = vgrid!(
    Label(fig, "Passes", width = nothing),
    menu; tellheight = false, width = 150)
fig[2, 2] = grid!(hcat(toggles, labels), tellheight = false)

textbox = fig[1:2, 1] = Axis(fig[1:2, 1], title = "Compiled code")
sl = Slider(fig[1, 2, Left()], range = 0.05:0.001:0.1, startvalue = 0.075, horizontal = false)
connect!(txtsize, sl.value)
# 4. Lift
#strtxt = Observable("n")
#colorsObs = Observable(RGBA{Float32}[RGBA{Float32}(0,0,0,1)])
str = @lift begin
    $(toggles[1].active)
    $(toggles[2].active)
    $(toggles[3].active)
    $(toggles[4].active)
    $(toggles[5].active)
    $(toggles[6].active)
    $(toggles[7].active)
    $(toggles[8].active)
    $(toggles[9].active)
    $(menu.selection)
    strtmp = string(main(menu))
    #highlight_text(strtmp, Themes.DefaultTheme)
    strtmp
end
#strtxt.val = @lift($str[1])[]
#colorsObs[] = @lift($str[2])[]

s = text!(textbox, str, color = :black, #colorsObs,
    font = "JuliaMono", align = (:left, :top),
    justification = :left, textsize = txtsize, space = :data)
#delete!(textbox.scene, s.val)
xlims!(textbox, -0.5, 10)
ylims!(textbox, -10, 0.5)
fig

# KNOWN ISSUES:
# 1. The text should be centered at the top left. I tried passing in `position` and using `campixel(fig.scene)`
#       but the sizes returned were smaller than those of the total window and I didn't know how to change it.
# 2. The `@lift` is surely suboptimal - I'm open to hearing about a better pattern.
# 3. The bounding box for the `menu` looks *very* wonky. I wish it was a smaller floating menu just above the `Toggle`s.
# 4. 👷 Me - I'm sure my coding can be improved/be redesigned for easier use and avoid some pain. Please advise.
#
# WANTED FEATURES:
# 1. Drag and Drop-ing the passes into a different order.
# 2. Add/delete! passes in a certain order and have recompilation be triggered.
#       🚧 If you have any tips about how to make these work, please reach out! 🚧
#
#                                                                       Till next time,
#                                                                         Cap'n Nemo 🐙