text:
function QRDecomposer(arg0::jdouble)
    return QRDecomposer((jdouble,), arg0)
end

function decompose(obj::QRDecomposer, arg0::RealMatrix)
    return jcall(obj, "decompose", DecompositionSolver, (RealMatrix,), arg0)
end


text:
#= ===================
   LOGISTIC CLASSIFIER
   =================== =#

"""
Logistic Classifier (typically called "Logistic Regression"). This model is
a standard classifier for both binary and multiclass classification.
In the binary case it corresponds to the LogisticLoss, in the multiclass to the
Multinomial (softmax) loss. An elastic net penalty can be applied with
overall objective function

``L(y, Xθ) + n⋅λ|θ|₂²/2 + n⋅γ|θ|₁``

where ``L`` is either the logistic or multinomial loss and ``λ`` and ``γ`` indicate
the strength of the L2 (resp. L1) regularisation components and
``n`` is the number of samples `size(X, 1)`.
With `scale_penalty_with_samples = false` the objective function is
``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁``

## Parameters

$TYPEDFIELDS

$(example_docstring("LogisticClassifier", nclasses = 2))
"""
@with_kw_noshow mutable struct LogisticClassifier <: MMI.Probabilistic
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::LogisticClassifier, nclasses::Integer) =
    LogisticRegression(m.lambda, m.gamma;
                       penalty=Symbol(m.penalty),
                       multi_class=(nclasses > 2),
                       fit_intercept=m.fit_intercept,
                       penalize_intercept=m.penalize_intercept,
                       scale_penalty_with_samples=m.scale_penalty_with_samples,
                       nclasses=nclasses)

descr(::Type{LogisticClassifier}) = "Classifier corresponding to the loss function ``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the logistic loss."

#= ======================
   MULTINOMIAL CLASSIFIER
   ====================== =#

"""
See `LogisticClassifier`, it's the same except that multiple classes are assumed
by default. The other parameters are the same.

## Parameters

$TYPEDFIELDS

$(example_docstring("LogisticClassifier", nclasses = 3))
"""
@with_kw_noshow mutable struct MultinomialClassifier <: MMI.Probabilistic
    "strength of the regulariser if `penalty` is `:l2` or `:l1`.
    Strength of the L2 regulariser if `penalty` is `:en`."
    lambda::Real             = 1.0
    "strength of the L1 regulariser if `penalty` is `:en`."
    gamma::Real              = 0.0
    "the penalty to use, either `:l2`, `:l1`, `:en` (elastic net) or `:none`."
    penalty::SymStr          = :l2
    "whether to fit the intercept or not."
    fit_intercept::Bool      = true
    "whether to penalize the intercept."
    penalize_intercept::Bool = false
    "whether to scale the penalty with the number of samples."
    scale_penalty_with_samples::Bool = true
    "type of solver to use, default if `nothing`."
    solver::Option{Solver}   = nothing
end

glr(m::MultinomialClassifier, nclasses::Integer) =
    MultinomialRegression(m.lambda, m.gamma;
                          penalty=Symbol(m.penalty),
                          fit_intercept=m.fit_intercept,
                          penalize_intercept=m.penalize_intercept,
                          scale_penalty_with_samples=m.scale_penalty_with_samples,
                          nclasses=nclasses)

descr(::Type{MultinomialClassifier}) =
    "Classifier corresponding to the loss function " *
    "``L(y, Xθ) + λ|θ|₂²/2 + γ|θ|₁`` where `L` is the multinomial loss."

ext: jl


text:
try
    using FoldsCUDATests
    true
catch
    false
end || begin
    let path = joinpath(@__DIR__, "FoldsCUDATests")
        path in LOAD_PATH || push!(LOAD_PATH, path)
    end
    let path = joinpath(@__DIR__, "../benchmark/FoldsCUDABenchmarks")
        path in LOAD_PATH || push!(LOAD_PATH, path)
    end
    using FoldsCUDATests
end


text:
# This file is part of IntegerTriangles.
# Copyright Peter Luschny. License is MIT.

# Version of: UTC 2021-05-22 22:01:34
# 8228b27e-bb38-11eb-3e8e-73d8fe660356

# Do not edit this file, it is generated from the modules and will be overwritten!
# Edit the modules in the src directory and build this file with BuildTriangles.jl!

tstdir = realpath(joinpath(dirname(@__FILE__)))
srcdir = joinpath(dirname(tstdir), "src")
tstdir ∉ LOAD_PATH && push!(LOAD_PATH, tstdir)
srcdir ∉ LOAD_PATH && push!(LOAD_PATH, srcdir)
module perftests
using IntegerTriangles, Dates, InteractiveUtils
InteractiveUtils.versioninfo()
start = Dates.now()
# +++ TrianglesBase.jl +++
# +++ TrianglesExamples.jl +++
T = LaguerreTriangle(8)
Println.(PolyTriangle(T))
Println.(PolyArray(T))
Println.(Inverse(PolyTriangle(T)))
Println.(FubiniTriangle(8))
Println.(DArcaisTriangle(8))
# +++ TrianglesExplorer.jl +++
# +++ TrianglesPlot.jl +++
# +++ TrianglesTables.jl +++
# +++ TrianglesTraitCard.jl +++
# +++ TrianglesUtils.jl +++
GetSeqnum(ℤInt[1, 1, -2, 3, -3, 3, -5, 7, -6, 6, -10,
12, -11, 13, -17, 20, -21, 21, -27, 34, -33, 36, -46, 51,
-53, 58, -68, 78, -82, 89, -104]) |> println
GetSeqnum(ℤInt[0, 1, 1, 1, 2, 1, 2, 1, 5, 2, 2, 1, 5, 1,
2, 1, 14, 1, 5, 1, 5, 2, 2, 1, 15, 2, 2, 5, 4, 1, 4, 1, 51,
1, 2, 1, 14, 1, 2, 2, 14, 1, 6, 1, 4, 2, 2, 1, 52]) |> println
GetSeqnum(ℤInt[1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
29668681, 410815351, 6238931821, 101560835377, 1765092183037,
32838929702671, 644215775792401]) |> println
GetSeqnum(ℤInt[1, 1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
29668681, 410815351, 6238931821, 101560835377, 1765092183037,
32838929702671, 644215775792401]) |> println
GetSeqnum(ℤInt[0, 1, 2, 7, 44, 361, 3654, 44207, 622552, 10005041,
180713290, 3624270839, 79914671748, 1921576392793, 50040900884366,
1403066801155039, 42142044935535536]) |> println
GetSeqnum(ℤInt[0, 70, 3783, 338475, 40565585, 6061961733,
1083852977811, 225615988054171, 53595807366038234, 14308700593468127485,
4241390625289880226714]) |> println
GetSeqnumUri(ℤInt[1, 1, 7, 37, 241, 2101, 18271, 201097, 2270017,
29668681, 410815351, 6238931821, 101560835377, 1765092183037,
32838929702671, 644215775792401]) |> println
stop = Dates.now()
tdiff = stop - start
println("\nJulia version: " * string(VERSION) )
println(start)
println("Total test time: ", tdiff)
end # module


text:
# load packages
using CSV;
using DataFrames;
using StatsBase;

# setup filepaths
path_source = string(@__DIR__,"\\..\\source");
path_dev = string(@__DIR__,"\\..\\dev");

df_symp_appended = CSV.read(joinpath(path_dev,"19-21VAERSSYMPTOMS.csv"), DataFrame);

# Each row in VAERSSYMPTOMS.csv is limited to 5 MedDRA symptom terms
# This means VAERS reports with > 5 symptoms will have multiple rows
# Collapse the 5 MedDRA symptom columns so that we have a single column containing an array of symptoms for each row

test = copy(df_symp_appended);

# rows with 1 symptom
test1 = filter(row -> any(ismissing, row[names(test)[4:5]]), test);
test1[!,:SYMPTOMS] = map(collect, zip(test1[!,:SYMPTOM1]));
select!(test1, [:VAERS_ID, :SYMPTOMS])

# rows with 2 symptoms
test2 = filter(row -> !any(ismissing, row[names(test)[2:5]]) && any(ismissing, row[names(test)[6:7]]), test)
test2[!,:SYMPTOMS] = map(collect, zip(test2[!,:SYMPTOM1],
                                      test2[!,:SYMPTOM2]));
select!(test2, [:VAERS_ID, :SYMPTOMS])

# rows with 3 symptoms
test3 = filter(row -> !any(ismissing, row[names(test)[2:7]]) && any(ismissing, row[names(test)[8:9]]), test)
test3[!,:SYMPTOMS] = map(collect, zip(test3[!,:SYMPTOM1],
                                      test3[!,:SYMPTOM2],
                                      test3[!,:SYMPTOM3]));
select!(test3, [:VAERS_ID, :SYMPTOMS])

# rows with 4 symptoms
test4 = filter(row -> !any(ismissing, row[names(test)[2:9]]) && any(ismissing, row[names(test)[10:11]]), test)
test4[!,:SYMPTOMS] = map(collect, zip(test4[!,:SYMPTOM1],
                                      test4[!,:SYMPTOM2],
                                      test4[!,:SYMPTOM3],
                                      test4[!,:SYMPTOM4]));
select!(test4, [:VAERS_ID, :SYMPTOMS])

# rows with 5 symptoms (no missing symptoms)
test5 = filter(row -> !any(ismissing, row[names(test)[2:11]]), test);
test5[!,:SYMPTOMS] = map(collect, zip(test5[!,:SYMPTOM1],
                                      test5[!,:SYMPTOM2],
                                      test5[!,:SYMPTOM3],
                                      test5[!,:SYMPTOM4],
                                      test5[!,:SYMPTOM5]));
select!(test5, [:VAERS_ID, :SYMPTOMS]);

# append and sort
test = reduce(vcat, [test1, test2, test3, test4, test5]);
sort!(test, :VAERS_ID);

# Collapse the rows so that we get a dict of symptoms for each VAERS_ID
# Step 1: Create the dict
vaers_id_to_symptoms_dict = Dict{Int, Set{String}}()
# Step 2: Populate the keys (VAERS_ID) of the dict
for rownumber in 1:size(test, 1)
    vaers_id = test[rownumber, :VAERS_ID]
    if !haskey(vaers_id_to_symptoms_dict, vaers_id)
        # this is the set where we will store all of the symptoms for this VAERS ID
        vaers_id_to_symptoms_dict[vaers_id] = Set{String}()
    end
end
# Step 3: Populate the values (SYMPTOMS) of the dict
for rownumber in 1:size(test, 1)
    vaers_id = test[rownumber, :VAERS_ID]
    symptoms = test[rownumber, :SYMPTOMS]
    for symptom in symptoms
        push!(vaers_id_to_symptoms_dict[vaers_id], symptom)
    end
end

# write to dev
CSV.write(joinpath(path_dev,"19-21VAERSSYMPTOMS_dev.csv"), vaers_id_to_symptoms_dict);


text:
function walk_dom(f, session::Session, x::JSCode, visited = IdDict())
    walk_dom(f, session, x.source, visited)
end

walk_dom(f, session::Session, x, visited = IdDict()) = f(x)

function walk_dom(f, session::Session, x::Union{Tuple, AbstractVector, Pair}, visited = IdDict()) where T
    get!(visited, x, nothing) !== nothing && return
    for elem in x
        walk_dom(f, session, elem, visited)
    end
end

function walk_dom(f, session::Session, x::Node, visited = IdDict())
    get!(visited, x, nothing) !== nothing && return
    for elem in children(x)
        walk_dom(f, session, elem, visited)
    end
    for (name, elem) in Hyperscript.attrs(x)
        walk_dom(f, session, elem, visited)
    end
end

const mime_order = MIME.((
    "text/html", "text/latex", "image/svg+xml", "image/png",
    "image/jpeg", "text/markdown", "application/javascript", "text/plain"
))

function richest_mime(val)
    for mimetype in mime_order
        showable(mimetype, val) && return mimetype
    end
    error("value not writable for any mimetypes")
end

repr_richest(x) = repr(richest_mime(x), x)


text:
export
ruleSPAdditionOutNGG,
ruleSPAdditionOutNGP,
ruleSPAdditionOutNPG,
ruleSPAdditionOutNPP,
ruleSPAdditionIn1GNG,
ruleSPAdditionIn1PNG,
ruleSPAdditionIn1GNP,
ruleSPAdditionIn1PNP,
ruleSPAdditionIn2GGN,
ruleSPAdditionIn2PGN,
ruleSPAdditionIn2GPN,
ruleSPAdditionIn2PPN,
ruleSPAdditionOutNSP,
ruleSPAdditionOutNPS,
ruleSPAdditionIn1PNS,
ruleSPAdditionIn1SNP,
ruleSPAdditionIn2PSN,
ruleSPAdditionIn2SPN,
ruleMAdditionNGG

function ruleSPAdditionOutNGG(
    msg_out::Nothing,
    msg_in1::Message{F1, V},
    msg_in2::Message{F2, V}) where {F1<:Gaussian, F2<:Gaussian, V<:Union{Univariate, Multivariate}}

    d_in1 = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_in1.dist)
    d_in2 = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_in2.dist)

    Message(V, GaussianMeanVariance, m=d_in1.params[:m] + d_in2.params[:m], v=d_in1.params[:v] + d_in2.params[:v])
end

function ruleSPAdditionIn2GGN(
    msg_out::Message{F1, V},
    msg_in1::Message{F2, V},
    msg_in2::Nothing) where {F1<:Gaussian, F2<:Gaussian, V<:Union{Univariate, Multivariate}}

    d_in1 = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_in1.dist)
    d_out = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_out.dist)

    Message(V, GaussianMeanVariance, m=d_out.params[:m] - d_in1.params[:m], v=d_out.params[:v] + d_in1.params[:v])
end

function ruleSPAdditionIn1GNG(
    msg_out::Message{F1, V},
    ::Nothing,
    msg_in2::Message{F2, V}) where {F1<:Gaussian, F2<:Gaussian, V<:Union{Univariate, Multivariate}}

    ruleSPAdditionIn2GGN(msg_out, msg_in2, nothing)
end

function ruleSPAdditionOutNGP(
    msg_out::Nothing,
    msg_in1::Message{F, V},
    msg_in2::Message{PointMass, V}) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    d_in1 = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_in1.dist)

    Message(V, GaussianMeanVariance, m=d_in1.params[:m] + msg_in2.dist.params[:m], v=d_in1.params[:v])
end

function ruleSPAdditionOutNPG(
    ::Nothing,
    msg_in1::Message{PointMass, V},
    msg_in2::Message{F, V}) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    ruleSPAdditionOutNGP(nothing, msg_in2, msg_in1)
end

function ruleSPAdditionIn1PNG(
    msg_out::Message{PointMass, V},
    msg_in1::Nothing,
    msg_in2::Message{F, V}) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    d_in2 = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_in2.dist)

    Message(V, GaussianMeanVariance, m=msg_out.dist.params[:m] - d_in2.params[:m], v=d_in2.params[:v])
end

function ruleSPAdditionIn2PGN(
    msg_out::Message{PointMass, V},
    msg_in1::Message{F, V},
    msg_in2::Nothing) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    ruleSPAdditionIn1PNG(msg_out, nothing, msg_in1)
end

function ruleSPAdditionIn1GNP(
    msg_out::Message{F, V},
    msg_in1::Nothing,
    msg_in2::Message{PointMass, V}) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    d_out = convert(ProbabilityDistribution{V, GaussianMeanVariance}, msg_out.dist)

    Message(V, GaussianMeanVariance, m=d_out.params[:m] - msg_in2.dist.params[:m], v=d_out.params[:v])
end

function ruleSPAdditionIn2GPN(
    msg_out::Message{F, V},
    msg_in1::Message{PointMass, V},
    msg_in2::Nothing) where {F<:Gaussian, V<:Union{Univariate, Multivariate}}

    ruleSPAdditionIn1GNP(msg_out, nothing, msg_in1)
end

function ruleSPAdditionOutNPP(
    msg_out::Nothing,
    msg_in1::Message{PointMass, V},
    msg_in2::Message{PointMass, V}) where V<:Union{Univariate, Multivariate}

    Message(V, PointMass, m=msg_in1.dist.params[:m] + msg_in2.dist.params[:m])
end

function ruleSPAdditionIn2PPN(
    msg_out::Message{PointMass, V},
    msg_in1::Message{PointMass, V},
    msg_in2::Nothing) where V<:Union{Univariate, Multivariate}

    Message(V, PointMass, m=msg_out.dist.params[:m] - msg_in1.dist.params[:m])
end

function ruleSPAdditionIn1PNP(
    msg_out::Message{PointMass, V},
    msg_in1::Nothing,
    msg_in2::Message{PointMass, V}) where V<:Union{Univariate, Multivariate}

    Message(V, PointMass, m=msg_out.dist.params[:m] - msg_in2.dist.params[:m])
end

function ruleSPAdditionOutNSP(
    msg_out::Nothing,
    msg_in1::Message{SampleList, Univariate},
    msg_in2::Message{PointMass, Univariate})

    new_s = msg_in2.dist.params[:m] .+ msg_in1.dist.params[:s]

    Message(Univariate, SampleList, s=new_s)
end

function ruleSPAdditionOutNSP(
    msg_out::Nothing,
    msg_in1::Message{SampleList, Multivariate},
    msg_in2::Message{PointMass, Multivariate})

    new_s = [msg_in2.dist.params[:m] .+ sample for sample in msg_in1.dist.params[:s]]

    Message(Multivariate, SampleList, s=new_s)
end

function ruleSPAdditionOutNPS(
    ::Nothing,
    msg_in1::Message{PointMass, V},
    msg_in2::Message{SampleList, V}) where {V<:Union{Univariate, Multivariate}}

    ruleSPAdditionOutNSP(nothing, msg_in2, msg_in1)
end

function ruleSPAdditionIn1PNS(
    msg_out::Message{PointMass, Univariate},
    msg_in1::Nothing,
    msg_in2::Message{SampleList, Univariate})

    new_s = msg_in2.dist.params[:m] .- msg_in1.dist.params[:s]

    Message(Univariate, SampleList, s=new_s)
end

function ruleSPAdditionIn1PNS(
    msg_out::Message{PointMass, Multivariate},
    msg_in1::Nothing,
    msg_in2::Message{SampleList, Multivariate})

    new_s = [msg_in2.dist.params[:m] .- sample for sample in msg_in1.dist.params[:s]]

    Message(Multivariate, SampleList, s=new_s)
end

function ruleSPAdditionIn2PSN(
    msg_out::Message{PointMass, V},
    msg_in1::Message{SampleList, V},
    msg_in2::Nothing) where {V<:Union{Univariate, Multivariate}}

    ruleSPAdditionIn1PNS(msg_out, nothing, msg_in1)
end

function ruleSPAdditionIn1SNP(
    msg_out::Message{SampleList, Univariate},
    msg_in1::Nothing,
    msg_in2::Message{PointMass, Univariate})

    new_s = -(msg_in2.dist.params[:m] .- msg_in1.dist.params[:s])

    Message(Univariate, SampleList, s=new_s)
end

function ruleSPAdditionIn1SNP(
    msg_out::Message{SampleList, Multivariate},
    msg_in1::Nothing,
    msg_in2::Message{PointMass, Multivariate})

    new_s = [-(msg_in2.dist.params[:m] .- sample) for sample in msg_in1.dist.params[:s]]

    Message(Multivariate, SampleList, s=new_s)
end

function ruleSPAdditionIn2SPN(
    msg_out::Message{SampleList, V},
    msg_in1::Message{PointMass, V},
    msg_in2::Nothing) where {V<:Union{Univariate, Multivariate}}

    ruleSPAdditionIn1SNP(msg_out, nothing, msg_in1)
end

function ruleMAdditionNGG(
    msg_out::Message{<:Gaussian, V},
    msg_in1::Message{<:Gaussian, V},
    msg_in2::Message{<:Gaussian, V}) where V<:Union{Univariate, Multivariate}

    d_out = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_out.dist)
    d_in1 = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_in1.dist)
    d_in2 = convert(ProbabilityDistribution{V, GaussianWeightedMeanPrecision}, msg_in2.dist)

    xi_out = d_out.params[:xi]
    W_out = d_out.params[:w]
    xi_in1 = d_in1.params[:xi]
    W_in1 = d_in1.params[:w]
    xi_in2 = d_in2.params[:xi]
    W_in2 = d_in2.params[:w]

    return ProbabilityDistribution(Multivariate, GaussianWeightedMeanPrecision, xi=[xi_in1+xi_out; xi_in2+xi_out], w=[W_in1+W_out W_out; W_out W_in2+W_out])
end


text:
## ----
# Find isolating intervals for real roots of a polynomial
# following
# Computing real roots of real polynomials
# Michael Sagraloff, Kurt Mehlhorn
# https://doi.org/10.1016/j.jsc.2015.03.004

## ---- Admissible point

# psuedo admissible point uses randomization to cut down number of operations
#
# This simple function allocates too much:
# * the evalpoly call allocates
# * the allocations in evalpoly are due to allocations in the interval
#   arithmetic implementation of `rmul!`.
# * these allocations can be removed by passing in 3 temporary variables
# * these can't be recycled here, as we don't have a mechanism to change
#   the working precision of a BigFloat number. (`set_prec!` will cause
#   segmentation faults)
# This is an issue then for: evalpoly and the various tests
# zerotest, onetest, newtontest, boundarytest
function admissiblepoint(p, m, ϵ, n, L′ = 53)
    λ = 2 + ceil(Int, log2(n))
    L = L′
    L = max(L′, maximum(precision∘float, p), precision(m))
    while true
        n′ = λ*n
        i = rand(-n′:n′)
        δ = epsᴸ(L-2) # 4*2⁻ᴸ
        ξ = setprecision(() -> m + i * ϵ/n′, L)
        𝐯ᵢ = evalpolyᴸ(ξ, p, L)
        (𝐯ᵢ.lo > δ || 𝐯ᵢ.hi < -δ) && return ξ
        L = 2L
    end

end

## ---- Tests (onetest, zerotest, newtontest, boundarytest, signtest)

# log2(a)
function tₐ(pa,L=precision(pa))
    apa = abs(pa)
    ϵ = epsᴸ(L)
    apa > ϵ ? ceil(Int, log2(apa)) : 2L
end

τ(p::Container{T}) where {T} = max(one(T), log2(norm(p, Inf)))
Mlog(x::T) where {T} = log2(max(one(T), abs(x)))

## Test if the interval 𝑰=[a,b] has a zero in it
function zerotest(p, I)
    a,b = I
    a >= b && return true
    L′ = max(maximum(precision∘float, p), maximum(precision, (a,b)))

    ta, tb = tₐ(a,L′), tₐ(b,L′)

    n = length(p) - 1
    n′ = sum(divrem(n,2))

    L′′ = max(24, L′, max(1, -min(ta-1, tb-1) + 2(n+1) + 1))
    L = ceil(Int, n + τ(p) + n * Mlog(a) + n * Mlog(b) + L′′)

    m = setprecision(() -> a + (b-a)/2, L)
    ϵ = epsᴸ(L)
    𝐚, 𝐦, 𝐛 = ball.((a,m,b), ϵ)

    𝐪 = ball.(deepcopy.(p), ϵ)
    bnda = descartesbound!(𝐪, 𝐚, 𝐦, L)
    !iszero(bnda) && return false

    ball!.(𝐪, p, ϵ)
    bndb = descartesbound!(𝐪, 𝐦, 𝐛, L)
    iszero(bndb) && return true

    return false

end

## test if the interval 𝑰 = [a,b] has exactly 1 root in it
function onetest(p, I)
    a,b = I
    a >= b && return false

    L′ = max(maximum(precision∘float, p), maximum(precision, (a,b)))

    n = length(p) - 1
    n′ = sum(divrem(n,2))

    ϵ = (b-a) * 2.0^(-ceil(Int, log2(n)+2)) #   (b-a) / (4n)
    m = setprecision(() ->a + (b-a)/2, L′)

    mstar = admissiblepoint(p, m, ϵ, n′)

    ta, tb,t = tₐ(a,L′), tₐ(b,L′),tₐ(mstar,L′)
    L′′ = max(24, L′, max(1, -min(ta-1,tb-1, t-1)) + 4(n+2))
    L = ceil(Int, n + τ(p) + n * Mlog(a) + n * Mlog(b) + L′′)

    ϵ = epsᴸ(L)
    𝐚,𝐦,𝐛 = ball.((a ,mstar, b), ϵ)
    𝐪 = ball.(deepcopy.(p), ϵ)

    bnda = descartesbound!(𝐪, 𝐚, 𝐦, L)

    ball!.(𝐪, p, ϵ)
    bndb = descartesbound!(𝐪, 𝐦, 𝐛, L)

    (bnda == -1 || bndb == -1) && throw(ArgumentError("L not big enough for descartesbound"))
    iszero(bnda) && isone(bndb) && return true
    isone(bnda) && iszero(bndb) && return true

    return false

end

# The newton test over I=[a,b] succeeds (informally) if:
# * a cluster of k roots (k ≈ λ below) is in a smaller interval J=[c,d]
# * there are no other roots nearby J
# Success means a new interval with size ≈≤ w(I)/N is returned
# Lemma 20 has this condition on success
# w(J) ≤ 2⁻¹³ w(I)/N
# k roots are in the one-circe region of Δ(J)
# the disk centered at m(I) with radius 2^(log(n)+10) * N * w(I) contains no other roots
function newtontest(p, p′, I, nᴵ)

    a, b = I
    nᴵ = min(5, nᴵ)
    N = 2^(2^nᴵ)

    n = length(p) - 1
    n′ = sum(divrem(n, 2))
    ϵ = 2.0^(-ceil(Int, (5 + log(n))))
    L::Int = Lᵐ = maximum(precision, (a,b))

    for (i,j) ∈ ((1,2),(1,3),(2,3))
        discard_ij = false
        double_cnt = 0
        L = Lᵐ
        while !discard_ij

            # double L until either (25) or (26) is satisfied for the pair i,j
            m, w = setprecision(() -> (a + (b-a)/2, b-a), L)
            uᵢ,uⱼ = setprecision(() -> (a + i*w/4,a + j*w/4), L)
            ξᵢ = admissiblepoint(p, uᵢ, ϵ * w, n′, L)
            ξⱼ = admissiblepoint(p, uⱼ, ϵ * w, n′, L)

            𝐩ᵢ, 𝐩ⱼ  = evalpolyᴸ(ξᵢ, p, L),  evalpolyᴸ(ξⱼ, p, L)
            𝐩′ᵢ,𝐩′ⱼ = evalpolyᴸ(ξᵢ, p′, L), evalpolyᴸ(ξⱼ, p′, L)

            L′ = max(L, maximum(precision, (ξᵢ, ξⱼ)))

            rdiv!(𝐩ᵢ,𝐩′ᵢ); rdiv!(𝐩ⱼ, 𝐩′ⱼ)
            𝚫 = 𝐩ᵢ - 𝐩ⱼ
            𝐯ᵢ, 𝐯ⱼ =  𝐩ᵢ, 𝐩ⱼ

            pᵢ, pⱼ = maxabs(𝐩ᵢ), maxabs(𝐩ⱼ) # store, used to compute new L

            if (min(minabs(𝐯ᵢ), minabs(𝐯ⱼ)) > w) || (4n*maxabs(𝚫) < w)  # (25)
                # discard_too_big_difference_too_small
                discard_ij = true # stop checking this pair

            elseif (max(maxabs(𝐯ᵢ), maxabs(𝐯ⱼ)) < 2w) && (8n*minabs(𝚫) > w) # (26)

                discard_ij = true # stop checking this pair

                # now set L to be big enough and estimate λ
                # L = ....²¹ (leaves diff between L an dactual < 1/(32N)

                L = max(24, ceil(Int, log2(n) + log2(N) +
                          Mlog(1/pᵢ) + Mlog(1/pⱼ) + max(0, -log2(w))))

                𝐩ᵢ,  𝐩ⱼ  = evalpolyᴸ(ξᵢ, p, L),  evalpolyᴸ(ξⱼ, p, L)
                𝐩′ᵢ, 𝐩′ⱼ = evalpolyᴸ(ξᵢ, p′, L), evalpolyᴸ(ξⱼ, p′, L)

                rdiv!(𝐩ᵢ,𝐩′ᵢ); rdiv!(𝐩ⱼ,𝐩′ⱼ)
                𝐯ᵢ,𝐯ⱼ = 𝐩ᵢ,𝐩ⱼ

                λ = 𝐯ᵢ
                rdiv!(λ, 𝐯ᵢ - 𝐯ⱼ)
                rmul!(λ, ξⱼ - ξᵢ)
                radd!(λ, ξᵢ)

                if λ.hi < a || λ.lo > b
                    continue
                elseif isinf(λ) || isnan(λ)
                    continue
                else

                    lᵢⱼ = setprecision(() ->floor(Int, (midpoint(λ) - a)/w * 4N), L)
                    aᵢⱼ = setprecision(() -> a + (lᵢⱼ - 1)  * w/(4N), L)
                    bᵢⱼ = setprecision(() -> a + (lᵢⱼ + 2) * w/(4N), L)

                    aᵢⱼᵅ = admissiblepoint(p, aᵢⱼ, ϵ *w/N, n′, L)
                    bᵢⱼᵅ = admissiblepoint(p, bᵢⱼ, ϵ *w/N, n′, L)

                    lᵢⱼ - 1 <= 0  && (aᵢⱼᵅ = Big(a, precision=L))
                    lᵢⱼ + 2 >= 4N && (bᵢⱼᵅ = Big(b, precision=L))

                    aᵢⱼᵅ, bᵢⱼᵅ = aᵢⱼᵅ < bᵢⱼᵅ  ? (aᵢⱼᵅ, bᵢⱼᵅ) : (bᵢⱼᵅ, aᵢⱼᵅ)

                    b-a < bᵢⱼᵅ - aᵢⱼᵅ && continue

                    bnda = zerotest(p, (a,  aᵢⱼᵅ ))
                    bndb = zerotest(p, (bᵢⱼᵅ, b))

                    if  bnda && bndb
                        return true, 𝐈(aᵢⱼᵅ,  bᵢⱼᵅ)
                    end

                end
            elseif double_cnt > 4
                # abandon and use  bisection
                discard_ij = true
            end

            double_cnt += 1
            L = 2L
        end
    end
    return false, 𝐈(a,b)

end

# The boundary test tests for a cluster of roots near an endpoint of I
# that is not well separated from other roots in Δ(I)
#
# return (true, J) can reduce interval with J ≈≤ w(I)/N
# return (false, I) can't trim, use linear
# no root in I should throw an error
function boundarytest(p, I, nᴵ)
    a, b = I
    nᴵ = min(5, nᴵ)
    N = 2^(2^nᴵ)

    n = length(p) - 1
    n′ = sum(divrem(n,2))
    ϵ = 1/2.0^(2 + ceil(Int, log(n)))

    L = maximum(precision, (a,b))

    w, m, mₗ, mᵣ = setprecision(L) do
        w = b-a
        m = a + w/2
        mₗ, mᵣ = a + w/(2N), b - w/(2N)
        w, m, mₗ, m
    end

    ϵ′ = ϵ*w/N
    mₗ⁺, mᵣ⁺ = admissiblepoint(p, mₗ, ϵ′, n′), admissiblepoint(p, mᵣ, ϵ′, n′)

    zₗ, zᵣ = zerotest(p, (mₗ⁺, b)), zerotest(p, (a, mᵣ⁺))
    zₗ && zᵣ && return (false, 𝐈(a,b)) #error("No root in I $a, $mᵣ⁺, $mₗ⁺  $b")
    zₗ && return (true, 𝐈(a, mₗ⁺))
    zᵣ && return (true, 𝐈(mᵣ⁺, b))
    return (false, 𝐈(a,b))

end

# true if `p` crosses in interval `I`
# false if not
# nothing if ambiguous
function signtest(p, I)
    a,b = I
    L = maximum(precision, (a,b))

    Σ = ball(0, L)
    𝐚 = evalpolyᴸ!(Σ, a, p, L)
    sa = 𝐚.hi < 0 ? -1 : 𝐚.lo > 0 ? 1 : return nothing
    𝐛 = evalpolyᴸ!(Σ, b, p, L)
    sb = 𝐛.hi < 0 ? -1 : 𝐛.lo > 0 ? 1 : return nothing

    sa * sb < 0

end

# split interval at an admissible point
function  bisect_interval(p, I, n′, L)
    setprecision(L) do
        a,b = I
        c = setprecision(L) do
            w = b-a
            m = 4 + ceil(Int, log2(n′))
            admissiblepoint(p, a + w/2, w/2^m, n′,L)
        end

        𝐈(a,c), 𝐈(c,b)
    end
end

## ----

## Upper bound on size of real roots that is tighter than cauchy
## titan.princeton.edu/papers/claire/hertz-etal-99.ps
## δ is buffer ensure interval endpoint do not include roots
function upperbound(p, δ = 1/2)

    T = BigFloat
    n = length(p) - 1
    n > 0 || return δ
    L = 2 + n + ceil(Int, τ(p)) #+ 53

    setprecision(L) do
        descartescount(p) == 0 && return zero(T) + δ
        q, d = p./p[end], n

        d == 0 && return δ
        d == 1 && return max(0, -q[1]) + δ

        a1 = abs(q[d])
        B = maximum([abs(q[i]) for i in 1:(d-1)])

        a,b,c = 1, -(1+a1), a1-B
        out = (-b + sqrt(b^2 - 4a*c))/2
        out + δ
    end

end

function lowerbound(p)
    q = [(isodd(i) ? 1 : -1)*p[i] for i ∈ eachindex(p)]
    return -upperbound(q)
end

root_bound(p) = (lowerbound(p), upperbound(p))

# Descartes bound on real roots in [0,oo) based on the sign changes
# treat 0 (|aᵢ| <= eps(T)) as no sign change
function descartescount(q::Container{T}, tol::T=zero(T)) where {T}
    n = length(q)

    i₀ = 1 # leading 0s
    while abs(q[i₀]) <= tol
        i₀ += 1
    end

    q₀ = q[i₀]
    flag = q₀ < 0 ? true : false

    cnt = 0

    for i ∈ (i₀+1):n
        qᵢ = q[i]
        abs(qᵢ) <= tol && continue
        flag′ = qᵢ < 0 ? true : false
        if flag′ != flag
            cnt += 1
            flag = flag′
        end
    end

    cnt

end

## ----

# hold the state
struct State{N,T,S}
    Isol::Vector{𝐈{S}}                        # DesBound == 1
    Unresolved::Vector{𝐈{S}}
    p::NTuple{N,T}
end

(st::State)(x) = evalpoly(x, st.p)

## iterate over Isol
Base.length(st::State) = length(st.Isol)
Base.iterate(st::State, state=nothing) = state==nothing ? iterate(st.Isol) : iterate(st.Isol, state)

function Base.show(io::IO, st::State)

    if !isempty(st.Unresolved)
        println(io, "There are unresolved intervals:")
        println.(Ref(io), st.Unresolved)
        println(io, "")
    end

    n = length(st.Isol)

    msg = n == 0 ? "No isolating intervals found." :
        n == 1 ?  "There was 1 isolating interval found:" :
        "There were $n isolating intervals found:"

    println(io, msg)
    n > 0 && println.(Ref(io), st.Isol)

    return nothing
end

## ----

"""
    ANewDsc(p; root_bound=(lowerbound(p), upperbound(p)), max_depth=96)
    refine_interval(p, a, b, L)
    refine_roots(st::State)

A method to find isolating intervals for the real roots of a
square-free polynomial specified by the cofficients stored in `p`.

* `p`: the polynomial coefficients, `[a₀, a₁, …, aₙ]`, of a **square-free** polynomial.
* `root_bound`: a lower bound and upper bound for the real roots.

Returns a `State` instance which has components:

* `Isol` holding the isolating intervals. Iteration over a `State` object will iterate over `Isol`.
* `Unresolved` holding any unresolved intervals. The show method alerts the presence of any such intervals.

The algorithm has a random step included, which leads to small variations in the output.

Examples:

```jldoctest anewdsc
julia> using RealPolynomialRoots

julia> ps = [-1, 254, -16129, 0, 0, 0, 0, 1] # mignotte polynomial with two nearby roots

julia> st = ANewDsc(ps)
There were 3 isolating intervals found:
[4.0…, 7.0…]₅₃
[0.00787401571051…, 0.00787401669368…]₆₄
[0.00787401494608…, 0.00787401571051…]₆₄

julia> ps = [3628800, -10628640, 12753576, -8409500, 3416930, -902055, 157773, -18150, 1320, -55, 1]; # πᵢ₌₁¹⁰ (x-i)

julia> st = ANewDsc(ps)
There were 10 isolating intervals found:
[9.0…, 10.5…]₅₃
[8.25…, 9.25…]₅₃
[7.62…, 8.25…]₅₃
[6.25…, 7.75…]₅₃
[5.5…, 6.38…]₅₃
[4.25…, 5.75…]₅₃
[3.44…, 4.38…]₅₃
[2.38…, 3.5…]₅₃
[1.5…, 2.38…]₅₃
[-0.0…, 1.5…]₅₃

julia> ps =[ # from https://discourse.julialang.org/t/root-isolation-of-real-rooted-integer-polynomials/51421/1
                      942438915208811912419937422298363203125
                   164182217245953398816894035758761902846875
                  4584900574568933770264468813466870772155175
                 48995332393110515074735075708247882042540865
                266674183150777010544241114017741621823207005
                852443280934837985352088128423887894438557515
               1738546146302892245736990844587487000484756535
               2381558158813900978436742173305983349418813145
               2262003889258248241081177038309445610985409335
               1516025051068561122302699855213604175083575145
                720810987764796866354279087485114863858085005
                241341213302325116160821160849326697595681275
                 55563538585119205063994483187179167406616375
                  8363912237256094118085070946083688310200625
                   740493466743082745510080711751444519503125
                    29215606371473169285018060091249259296875];

julia> ANewDsc(ps)
There were 15 isolating intervals found:
[-0.062195…, 0.0…]₅₁₂
[-0.1421…, -0.06226…]₅₁₂
[-0.2529…, -0.1426…]₅₁₂
[-0.5078…, -0.2539…]₅₁₂
[...]
```

## Refinement

The `refine_interval` method can be used to refine an interval to have
width smaller than ``2^{-L}`` where `L` may be specified, but
otherwise comes from the intervals precision.

Alternatively, a package like `Roots` could be used; e.g:
`[find_zero(st, I) for I ∈ st]` (where `st` is a `State` object
returned by `ANewDsc`). If refinement over `Float64` values is desired
and appropriate given the root separation, then that call can be
modified, as with `[find_zero(st, Float64.(I)) for I ∈ st]`. (This
should produce roots with a sign change between `nextfloat` and
`prevfloat`.)

## Comparisons

Don't judge the algorithm by its implementation here. This
implementation is not as performant as it could be.

Comparing to some alternatives, we have that the functionality from
Hecke.jl (`Hecke._roots`) is **much** better. (The last example is 33
≈ 0.14s/0.0042 times faster)

However, compared to other alternatives this implementation could be seen as useful:

```
julia> x = variable(Polynomial);

julia> p = -1 + 254*x - 16129*x^2 + x^15;

julia> ANewDsc(coeffs(p))  # ≈ 0.05 seconds;
There were 3 isolating intervals found:
[0.75…, 4.25…]₅₃
[0.00787401574803149653139…, 0.00787401574803149972047…]₂₁₂
[0.0078740157480314937167…, 0.00787401574803149653139…]₂₁₂

julia> filter(isreal, roots(p)) # much faster, but misses two roots with imaginary part ~ 1e-10
1-element Array{Complex{Float64},1}:
 2.1057742291764914 + 0.0im

julia> filter(isreal, AMRVW.roots(Float64.(coeffs(p)))) # using AMRVW. Similarly misses two roots
1-element Array{Complex{Float64},1}:
 2.1057742291764407 + 0.0im

julia> filter(isreal, AMRVW.roots(BigFloat.(coeffs(p)))) # this works here
3-element Vector{Complex{BigFloat}}:
 0.007874015748031494751793842937491860399146218747427882112208862187156603046408902 + 0.0im
 0.007874015748031497374190409031015351762713667398750832139185145345899098824178322 + 0.0im
    2.105774229176482954331883383107195997983004314462374928263620342390986189650864 + 0.0im

julia> filter(isreal, PolynomialRoots.roots(coeffs(p))) # using PolynomialRoots. Misses 2.105?
2-element Array{Complex{Float64},1}:
   0.0078740158234482 + 0.0im
 0.007874015672614794 + 0.0im

julia> IntervalRootFinding.roots(x->p(x), IntervalArithmetic.Interval(0.0, 5.0)) # using IntervalRootFinding, IntervalArithmetic
8-element Array{Root{IntervalArithmetic.Interval{Float64}},1}:
 Root([2.10577, 2.10578], :unique)
 Root([0.00787418, 0.00787422], :unknown)
 Root([0.00787403, 0.00787409], :unknown)
[...]

julia> @syms x::real # using SymPy

julia> @time  rts = sympy.real_roots(p(x)); # correctly identifies 3.
  0.003896 seconds (518 allocations: 13.359 KiB)

julia> sympy.N(rts[2]) # takes a long time! (162s)
0.00787401574803150
```

The algorithm used is a partial implementation of one presented in:

Computing Real Roots of Real Polynomials ... and now For Real!
by Alexander Kobel, Fabrice Rouillier, Michael Sagraloff
arXiv:1605.00410; DOI:	10.1145/2930889.2930937

and

Computing real roots of real polynomials
Michael Sagraloff, Kurt Mehlhorn
https://doi.org/10.1016/j.jsc.2015.03.004

The algorithm relies on Descartes' rule of signs, which gives a bound
on the number of positive real roots of a polynomial, `p(x)`. By
considering the polynomial `p((ax+b)/(x+1))` (a mapping of `[a,b]` to
`[0,∞)` a bound on the number of roots in `[a,b]` can be found. A
simple bisection algorithm for a square-free polynomial would be to
start with an interval large enough to contain all the real roots,
then subdivide at the midpoint throwing out sub intervals which are
found to have no root; saving intervals known to have 1 root, and
repeating the subdivision otherwise. Issues with this are the need for
many subdivisions when clusters of roots are present and the numeric
issues that arise in computing the mapping.

The work of Sagaraloff, Melhorn, Kobel, and Rouillier improves this by
introducing a Newton test and boundary test for rapidly decreasing the
size of an interval when possible, and the ability to use finite
precision arithmetic, instead of exact arithmetic, to compute the
Descartes' bound, in addition to other algorithmic improvements (not
all implemented here).

!!! Note

    A square free polynomial can be found through `p/gcd(p, p')`,
    though in practice this calculation is numerically unstable.

!!! Note

    This implementation is **much** slower than the `Hecke.roots`
    function provided through `arblib` in `Hecke.jl`, which itself
    says is not competitive with more specialized implementations,
    such as provided by the paper authors
    (http://anewdsc.mpi-inf.mpg.de/). There are several reasons: The
    `mobius_transform!` function is 𝑶(n²), and could be 𝑶(n⋅log(n))
    with more effort; the polynomial evaluation in `admissiblepoint`
    could, similarly, be made more efficient; despite using tricks
    learned from the `MutableArithmetics.jl` package to reduce
    allocations with the `BigFloat` type, there are still *far* too
    many allocations as each time the precision is changed new
    (allocating) values must be created, as the old ones can't be
    modified in place (using `set_prec!` causes segfaults); the
    significant engineering speedups suggested by Kobel, Rouillier,
    and Sagraloff are not implemented; etc.

"""
ANewDsc, refine_interval, refine_roots

function ANewDsc(p::Container{<:Real}; root_bound=root_bound(p), max_depth=32)

    T = BigFloat

    n = length(p) - 1
    iszero(n) && return State(𝐈{T}[],  𝐈{T}[], NTuple{n+1,eltype(p)}(p))

    max_depth *= ceil(Int, log2(n+2))

    n′ = sum(divrem(n,2))
    p′ = [i*p[i+1] for i ∈ 1:n]

    m′,M′ = T.(root_bound)
    I = 𝐈(m′, M′)

    Internal = [(I,1,1)]
    Isol = 𝐈{T}[]
    Unresolved = 𝐈{T}[]
    ctr = [0,0,0]
    while !isempty(Internal)

        I, nᴵ, depth = pop!(Internal)
        depth += 1
        L = precision(I)

        # 01-test is decisive
        zerotest(p, I) && continue
        onetest(p, I) && (push!(Isol, I); continue)

        ## pump the brakes, if needed, otherwise can go on and on in a newtonstep
        ## with non-square input
        depth >= 10max_depth && (push!(Unresolved, I); continue)

        # shrink or divide;
        val, J = newtontest(p, p′, I, nᴵ)
        if val
            push!(Internal, (J,nᴵ+1, depth))
            continue
        end

        val, J = boundarytest(p, I, nᴵ)
        if val
            push!(Internal, (J,nᴵ,depth))
            continue
        end

        # bisection
        Iₗ,Iᵣ = bisect_interval(p, I, n′, L)
        nᴵ = max(0, nᴵ-1)
        append!(Internal, ((Iₗ,nᴵ,depth), (Iᵣ,nᴵ,depth)))

    end

    return State(Isol, Unresolved, NTuple{n+1,eltype(p)}(p))

end


text:
using Documenter, BasicBSpline

# Setup for doctests in docstrings
DocMeta.setdocmeta!(BasicBSpline, :DocTestSetup, :(using LinearAlgebra, BasicBSpline))

function generate_indexmd_from_readmemd()
    path_readme = "README.md"
    path_index = "docs/src/index.md"

    # open README.md
    f = open(path_readme)
    text_readme = read(f,String)
    close(f)

    # generate text for index.md
    text_index = replace(text_readme,"![](docs/src/img" => "![](img")

    # save index.md
    open(path_index, "w") do f
        write(f,text_index)
    end
end

generate_indexmd_from_readmemd()

makedocs(;
    modules = [BasicBSpline],
    format = Documenter.HTML(
        canonical = "https://hyrodium.github.io/BasicBSpline.jl/stable/",
        assets = ["assets/favicon.ico"],
    ),
    pages = [
        "Home" => "index.md",
        "Mathematical properties of B-spline" => "math.md",
        "Benchmark" => "benchmark.md",
        "Example 1 - Poisson's equation" => "poisson.md",
        "Example 2 - Linear elasticity" => "elasticity.md",
        "Docstrings" => "detail.md",
        "Future work" => "future.md",
        "Contributing" => "contributing.md",
    ],
    repo = "https://github.com/hyrodium/BasicBSpline.jl/blob/{commit}{path}#L{line}",
    sitename = "BasicBSpline.jl",
    authors = "hyrodium <hyrodium@gmail.com>",
    assets = String[],
)

deploydocs(; repo = "github.com/hyrodium/BasicBSpline.jl")


text:
struct RegistryInstance
    path::String
    name::String
    uuid::UUID
    url::Union{String, Nothing}
    repo::Union{String, Nothing}
    description::Union{String, Nothing}
    pkgs::Dict{AbstractString, PkgEntry}
end

function RegistryInstance(path::AbstractString)
    d = parsefile(joinpath(path, "Registry.toml"))
    pkgs = Dict{AbstractString, PkgEntry}()

    for (uuid, info) in d["packages"]
        uuid = UUID(uuid)
        info
        name = info["name"]
        pkgpath = info["path"]
        pkg = PkgEntry(pkgpath, path, name, uuid)
        pkgs[name] = pkg
    end

    return RegistryInstance(
        path,
        d["name"],
        UUID(d["uuid"]),
        get(d, "url", nothing),
        get(d, "repo", nothing),
        get(d, "description", nothing),
        pkgs
    )
end


text:
using ResumableFunctions
abstract type AbstractTestClass end
@resumable function generator_func()
    num = 1
    @yield num
    num = 5
    @yield num
    num = 10
    @yield num
end

@resumable function generator_func_loop()
    num = 0
    for n = 0:2
        @yield num + n
    end
end

@resumable function generator_func_loop_using_var()
    num = 0
    end_ = 2
    end_ = 3
    for n = 0:end_-1
        @yield num + n
    end
end

@resumable function generator_func_nested_loop()
    for n = 0:1
        for i = 0:1
            @yield (n, i)
        end
    end
end

@resumable function file_reader(file_name::String)
    for file_row in readline(file_name)
        @yield file_row
    end
end

@resumable function testgen()
    println("first")
    @yield 1
    println("second")
    @yield 2
end

@resumable function fib()
    a = 0
    b = 1
    while true
        @yield a
        a, b = (b, a + b)
    end
end

mutable struct TestClass <: AbstractTestClass
end
@resumable function generator_func(self::TestClass)
    num = 123
    @yield num
    num = 5
    @yield num
    num = 10
    @yield num
end

if abspath(PROGRAM_FILE) == @__FILE__
    arr1 = []
    for i in generator_func()
        push!(arr1, i)
    end
    @assert(arr1 == [1, 5, 10])
    arr2 = []
    for i in generator_func_loop()
        push!(arr2, i)
    end
    @assert(arr2 == [0, 1, 2])
    arr3 = []
    for i in generator_func_loop_using_var()
        push!(arr3, i)
    end
    @assert(arr3 == [0, 1, 2])
    arr4 = []
    testClass1::TestClass = TestClass()
    for i in generator_func(testClass1)
        push!(arr4, i)
    end
    @assert(arr4 == [123, 5, 10])
    arr5 = []
    for i in generator_func_nested_loop()
        push!(arr5, i)
    end
    @assert(arr5 == [(0, 0), (0, 1), (1, 0), (1, 1)])
    arr6 = []
    for res in file_reader("C:/Users/Miguel Marcelino/Desktop/test.txt")
        push!(arr6, res)
    end
    @assert(arr6 == ["test\n", "test\n", "test"])
    arr7 = []
    res = fib()
    for i = 0:5
        push!(arr7, res)
    end
    @assert(arr7 == [0, 1, 1, 2, 3, 5])
    for i in testgen()
        println(i)
    end
end


text:
module Log

export pssfss_logger, ScreenLevel, FileLevel, @logfile, @logscreen

using Logging: LogLevel, global_logger, @logmsg
using LoggingExtras: FormatLogger, TeeLogger, EarlyFilteredLogger

const ScreenLevel = LogLevel(-500)
const FileLevel = LogLevel(-400)

macro logscreen(x)
    :(@logmsg(ScreenLevel, $(esc(x))))
end

macro logfile(x)
    :(@logmsg(FileLevel, $(esc(x))))
end

function newlevels_filter(log_args)
    log_args.level ∈ (FileLevel, ScreenLevel)
end

function notnew_filter(log_args)
    log_args.level ∉ (ScreenLevel,FileLevel)
end

function plain_logger(logfile)
    isfile(logfile) && rm(logfile)
    FormatLogger() do io, args
        open(logfile,"a") do f
            println(f, args.message)
        end
        args.level == ScreenLevel && println(args.message)
    end
end

pssfss_logger(logfile) =
    TeeLogger(
        EarlyFilteredLogger(notnew_filter, global_logger()),
        EarlyFilteredLogger(newlevels_filter, plain_logger(logfile))
             )

end # module


text:
####################################################################
#  JSON validation
####################################################################

macro doassert(asserts::Symbol, key::String, block::Expr)
  quote
    if haskey($(esc(asserts)), $key)
      $(esc(:( keyval = ($asserts)[$key]) ))
      $(esc(block))
    end
  end
end

macro report(what, wher, msg)
  quote
    return SingleIssue($(esc(what)), $(esc(wher)), $(esc(msg)))
  end
end

abstract type AbstractIssue end

struct SingleIssue <: AbstractIssue  # validation issue reporting structure
  x                     # subJSON with validation issue
  path::Vector{String}  # JSPointer to x
  msg::String           # error message
end

struct OneOfIssue <: AbstractIssue # validation issue reporting structure
  path::Vector{String}  # JSPointer to x
  issues::Vector{<:AbstractIssue} # potential issues, collectively making the present issue
end

### type assertion
function type_asserts(x, typ::String, path)
  if typ == "string"
    x isa String || @report x path "is not a string"
  elseif typ == "number"
    (x isa Int) || (x isa Float64) || @report x path "is not a number"
  elseif typ == "integer"
    x isa Int || @report x path "is not an integer"
  elseif typ == "array"
    x isa Array || @report x path "is not an array"
  elseif typ == "boolean"
    x isa Bool || @report x path "is not a boolean"
  elseif typ == "object"
    x isa Dict || @report x path "is not an object"
  elseif typ == "null"
    (x == nothing) || @report x path "is not null"
  end
  nothing
end

### check for assertions applicable to arrays
function array_asserts(x, asserts, path)
  @doassert asserts "items" begin
    if (keyval isa Bool) && ! keyval
      (length(x) > 0) && @report x path "schema (=false) does not allow non-empty arrays"

    elseif keyval isa Dict
      for i in 1:length(x)
        ret = validate(x[i], keyval, [path; "$(i-1)"])
        (ret==nothing) || return ret
      end

    elseif keyval isa Array
      for (i, iti) in enumerate(keyval)
        i > length(x) && break
        ret = validate(x[i], iti, [path; "$(i-1)"])
        (ret==nothing) || return ret
      end

      if haskey(asserts, "additionalItems") && (length(keyval) < length(x))
        addit = asserts["additionalItems"]
        if addit isa Bool
          (addit == false) &&  @report x path "additional items not allowed"
        else
          for i in length(keyval)+1:length(x)
            ret = validate(x[i], addit, [path; "$(i-1)"])
            (ret==nothing) || return ret
          end
        end
      end
    end
  end

  @doassert asserts "maxItems" begin
    (length(x) <= keyval) ||
      @report x path "array longer than $keyval"
  end

  @doassert asserts "minItems" begin
    (length(x) >= keyval) ||
      @report x path "array shorter than $keyval"
  end

  @doassert asserts "uniqueItems" begin
    if keyval
      xt = [(xx, typeof(xx)) for xx in x]
      (length(unique(hash.(xt))) == length(x)) ||
        @report x path "non unique elements"
    end
  end

  @doassert asserts "contains" begin
    is = AbstractIssue[]
    oneOK = false
    for (i,el) in enumerate(x)
      ret = validate(el, keyval, [path; "$(i-1)"])
      (ret==nothing) && (oneOK = true; break)
      push!(is, ret)
    end
    oneOK || return OneOfIssue(path, is)
  end

  nothing
end

### check for assertions applicable to numbers
function number_asserts(x, asserts, path)
  @doassert asserts "multipleOf" begin
    δ = abs(x - keyval*round(x / keyval))
    (δ < eps(Float64)) || @report x path "not a multipleOf of $keyval"
  end

  @doassert asserts "maximum" begin
    (x <= keyval) || @report x path "not <= $keyval"
  end

  @doassert asserts "exclusiveMaximum" begin
    if (keyval isa Bool) && haskey(asserts, "maximum")  # draft 4
      val2 = asserts["maximum"]
      if keyval
        (x < val2) || @report x path "not < $val2"
      else
        (x <= val2) || @report x path "not <= $val2"
      end
    elseif keyval isa Number  # draft 6
      (x < keyval) || @report x path "not < $keyval"
    end
  end

  @doassert asserts "minimum" begin
    (x >= keyval) || @report x path "not >= $keyval"
  end

  @doassert asserts "exclusiveMinimum" begin
    if (keyval isa Bool) && haskey(asserts, "minimum")  # draft 4
      val2 = asserts["minimum"]
      if keyval
        (x > val2) || @report x path "not > $val2"
      else
        (x >= val2) || @report x path "not >= $val2"
      end
    elseif keyval isa Number
      (x > keyval) || @report x path "not > $keyval"
    end
  end

  nothing
end

### check for assertions applicable to strings
function string_asserts(x, asserts, path)
  @doassert asserts "maxLength" begin
    (length(x) > keyval) && @report x path "string longer than $keyval characters"
  end

  @doassert asserts "minLength" begin
    (length(x) < keyval) && @report x path "string shorter than $keyval characters"
  end

  @doassert asserts "pattern" begin
    occursin(Regex(keyval), x) || @report x path "string does not match pattern $keyval"
  end

  nothing
end

### check for assertions applicable to objects
function object_asserts(x, asserts, path)
  if haskey(asserts, "dependencies")
    dep = asserts["dependencies"]
    for (k,v) in dep
      if haskey(x, k)
        if (v isa Dict) || (v isa Bool)
          ret = validate(x, v, path)
          (ret == nothing) || return ret
        else
          for rk in v
            haskey(x, rk) || @report x path "required property '$rk' missing (dependency $k)"
          end
        end
      end
    end
  end

  @doassert asserts "propertyNames" begin
    if isa(keyval, Bool)
      if ! keyval
        (length(x) > 0) && @report x path "schema (=false) allows only empty objects here"
      end
    else
      for propname in keys(x)
        ret = string_asserts(propname, keyval , path)
        (ret==nothing) || return ret
      end
    end
  end

  @doassert asserts "maxProperties" begin
    (length(x) <= keyval) || @report x path "nb of properties > $keyval"
  end

  @doassert asserts "minProperties" begin
    (length(x) >= keyval) || @report x path "nb of properties < $keyval"
  end

  remainingprops = collect(keys(x))
  @doassert asserts "required" begin
    for k in keyval
      (k in remainingprops) || @report x path "required property '$k' missing"
    end
  end

  matchedprops = String[]
  if haskey(asserts, "patternProperties") && length(remainingprops) > 0
    patprop = asserts["patternProperties"]
    for k in remainingprops
      hasmatch = false
      for (rk, rtest) in patprop
        if occursin(Regex(rk), k)
          hasmatch = true
          ret = validate(x[k], rtest, [path; k])
          (ret == nothing) || return ret
          # "property $k matches pattern $rk but is not a $rtest"
        end
      end
      hasmatch && push!(matchedprops, k)
    end
  end

  if haskey(asserts, "properties")
    prop = asserts["properties"]
    for k in intersect(remainingprops, keys(prop))
      ret = validate(x[k], prop[k], [path; k])
      (ret == nothing) || return ret
      push!(matchedprops, k)
    end
  end

  remainingprops = setdiff(remainingprops, matchedprops)
  if haskey(asserts, "additionalProperties") && length(remainingprops) > 0
    addprop = asserts["additionalProperties"]
    if addprop isa Bool
      rpstr = join([ "`$p`" for p in remainingprops], ", ", " & ")
      (addprop == false) && @report x path "additional property(ies) $rpstr not allowed"
    else
      for k in remainingprops
        ret = validate(x[k], addprop, [path; k])
        (ret == nothing) || return ret
      end
    end
  end

  nothing
end

# s = spec
# for cases where there is no schema but a true/false (in allOf, etc..)
function validate(x, s::Bool, path=String[])
  s || @report x path "schema (=false) does not allow any value here"
  nothing
end

validate(x, s::Schema, path=String[]) = validate(x, s.data, path)
function validate(x, asserts::Dict, path=String[])
  # if a ref is present, it should supersede all sibling properties
  refhistory = Any[asserts];
  while isa(asserts, Dict) && haskey(asserts, "\$ref") # resolve nested refs until an end is found
    asserts = asserts["\$ref"]
    (asserts in refhistory) && error("circular references in schema")
    push!(refhistory, asserts)
  end
  isa(asserts, Bool) && return validate(x, asserts, path) # in case ref turns out to be a boolean

  @doassert asserts "type" begin
    if keyval isa Array
      any( type_asserts(x, typ2, path)==nothing for typ2 in keyval ) ||
        @report x path "is not any of the allowed types $keyval"
    else
      ret = type_asserts(x, keyval, path)
      (ret==nothing) || return ret
    end
  end

  @doassert asserts "enum" begin
    any(x == e for e in keyval) ||
      @report x path "expected to be one of $keyval"
  end

  @doassert asserts "const" begin
    expected = keyval==nothing ? "'nothing'" : keyval
    (x == keyval) || @report x path "expected to be equal to $expected"
  end

  if isa(x, Array)
    ret = array_asserts(x, asserts, path)
    (ret==nothing) || return ret
  end

  if isa(x, Dict)
    ret = object_asserts(x, asserts, path)
    (ret==nothing) || return ret
  end

  if (x isa Int) || (x isa Float64)  # a 'number' for JSON
    ret = number_asserts(x, asserts, path)
    (ret==nothing) || return ret
  end

  if x isa String
    ret = string_asserts(x, asserts, path)
    (ret==nothing) || return ret
  end

  @doassert asserts "allOf" begin
    for subsch in keyval
      ret = validate(x, subsch, path)
      (ret==nothing) || return ret
    end
  end

  @doassert asserts "anyOf" begin
    is = AbstractIssue[]
    oneOK = false
    for subsch in keyval
      ret = validate(x, subsch, path)
      (ret==nothing) && (oneOK = true; break)
      push!(is, ret)
    end
    oneOK || return OneOfIssue(path, is)
  end

  @doassert asserts "oneOf" begin
    is = AbstractIssue[]
    oneOK, moreOK = false, false
    for subsch in keyval
      ret = validate(x, subsch, path)
      if ret==nothing
        oneOK && (moreOK = true)
        oneOK = true
      else
        push!(is, ret)
      end
    end
    oneOK || return OneOfIssue(path, is)
    # TODO : improve reporting if more than one matches
    moreOK && @report x path "more than one match in a 'oneOf' clause"
  end

  @doassert asserts "not" begin
    ret = validate(x, keyval, path)
    (ret==nothing) && @report x path "does not satisfy 'not' assertion $keyval"
  end

  nothing
end

##############   user facing functions #########################

"""
`isvalid(x, s::Schema)`

Check that a document `x` is valid against the Schema `s`.

## Examples

```julia
julia> myschema = Schema("
  {
    \"properties\": {
      \"foo\": {},
      \"bar\": {}
    },
    \"required\": [\"foo\"]
  }
  ")

julia> isvalid( JSON.parse("{ \"foo\": true }"), myschema)
true
julia> isvalid( JSON.parse("{ \"bar\": 12.5 }"), myschema)
false
```
"""
function isvalid(x, s::Schema)
  validate(x, s.data) == nothing
end

flatten(ofi::JSONSchema.OneOfIssue) = vcat([flatten(i) for i in ofi.issues]...)
flatten(si::JSONSchema.SingleIssue) = [si;]

singleissuerecap(si::JSONSchema.SingleIssue) =
    "in [$(join(si.path, '.'))] : $(si.msg)"

"""
`diagnose(x, s::Schema)`

Check that a document `x` is valid against the Schema `s`. If
valid return `nothing`, and if not, return a diagnostic String containing a
selection of one or more likely causes of failure.

## Examples

```julia
julia> diagnose( JSON.parse("{ \"foo\": true }"), myschema)
nothing
julia> diagnose( JSON.parse("{ \"bar\": 12.5 }"), myschema)
"in [] : required property 'foo' missing"
```
"""
function diagnose(x, s::Schema)
    hyps = JSONSchema.validate(x, s)
    (hyps == nothing) && return nothing

    hyps2 = flatten(hyps)

    # The selection heuristic is to keep only the issues appearing deeper in
    # the tree. This will trim out the 'oneOf' assertions that were not
    # intended in the first place in 'x' (hopefully).
    lmax = maximum(e -> length(e.path), hyps2)
    filter!(e -> length(e.path) == lmax, hyps2)

    if length(hyps2) == 1
        return singleissuerecap(hyps2[1])
    else
        msg = ["One of :";
               map(x -> "  - " * singleissuerecap(x), hyps2) ]
        return join(msg, "\n")
    end
    nothing
end


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking, CmdStan
#gr(size=(600,600));

ProjDir = rel_path("..", "scripts", "05")
cd(ProjDir)

# ### snippet 5.1

wd = CSV.read(rel_path("..", "data", "WaffleDivorce.csv"), delim=';');
df = convert(DataFrame, wd);
df[!, :A] = scale(df[!, :MedianAgeMarriage]);
df[!, :D] = scale(df[!, :Divorce]);
first(df, 5)

# ### snippet 5.1

std(df[!, :MedianAgeMarriage])

ad = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] D; // Outcome
 vector[N] A; // Predictor
}

parameters {
 real a; // Intercept
 real bA; // Slope (regression coefficients)
 real < lower = 0 > sigma; // Error SD
}

model {
  vector[N] mu;
  a ~ normal(0, 0.2);
  bA ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bA * A;
  D ~ normal(mu , sigma);
}
";

# Define the Stanmodel and set the output format to :mcmcchains.

m5_1s = Stanmodel(name="MedianAgeDivorce", model=ad);

# Input data for cmdstan

data = Dict("N" => length(df[!, :D]), "D" => df[!, :Divorce],
    "A" => df[!, :A]);

# Sample using cmdstan

rc, chn, cnames = stan(m5_1s, data, ProjDir, diagnostics=false,
  summary=true, CmdStanDir=CMDSTAN_HOME);

# Describe the draws

MCMCChains.describe(chn)

# Plot the density of posterior draws

plot(chn)

# Result rethinking

rethinking = "
       mean   sd  5.5% 94.5% n_eff Rhat
a      9.69 0.22  9.34 10.03  2023    1
bA    -1.04 0.21 -1.37 -0.71  1882    1
sigma  1.51 0.16  1.29  1.79  1695    1
"

# Plot regression line using means and observations

xi = -3.0:0.01:3.0
rws, vars, chns = size(chn)
alpha_vals = convert(Vector{Float64}, reshape(chn.value[:, 1, :], (rws*chns)))
beta_vals = convert(Vector{Float64}, reshape(chn.value[:, 2, :], (rws*chns)))
yi = mean(alpha_vals) .+ mean(beta_vals)*xi

scatter(df[!, :A], df[!, :D], color=:darkblue,
  xlab="Standardized median age of marriage",
  ylab="Standardize divorce rate")
plot!(xi, yi, lab="Regression line")

# shade(), abline() and link()

mu = link(xi, chn, [1, 2], mean(xi));
yl = [minimum(mu[i]) for i in 1:length(xi)];
yh =  [maximum(mu[i]) for i in 1:length(xi)];
ym =  [mean(mu[i]) for i in 1:length(xi)];
pi = hcat(xi, yl, ym, yh);
pi[1:5,:]

plot!((xi, yl), color=:lightgrey, leg=false)
plot!((xi, yh), color=:lightgrey, leg=false)
for i in 1:length(xi)
  plot!([xi[i], xi[i]], [yl[i], yh[i]], color=:lightgrey, leg=false)
end
scatter!(df[!, :A], df[!, :D], color=:darkblue)
plot!(xi, yi, lab="Regression line")

# End of `05/m5.1s.jl`


text:
struct Dual{T}
    x::T # primal ?
    der::T # derivative
end

# """ a and b are functions
# $ c(x) = a(x) + b(x) $
# $ dc/dx = da/dx + db/dx $
# """

Base.:+(a::Dual, b::Dual) = Dual(a.x + b.x, a.der + b.der)

Base.:+(a::Dual, b::Number) = Dual(a.x + b, a.der)
Base.:+(a::Number, b::Dual) = b + a

Base.:-(a::Dual, b::Dual) = Dual(a.x - b.x, a.der - b.der)

Base.:-(a::Dual, b::Number) = Dual(a.x - b, a.der)
Base.:-(a::Number, b::Dual) = Dual(a - b.x, -b.der)

Base.:-(a::Dual) = Dual(-a.x, -a.der)
# Base.:inv(a::Dual) = Dual(1/a.x, -a.der)

Base.:*(a::Dual, b::Dual) = Dual(a.x * b.x, a.der * b.x + a.x * b.der) # product rule

Base.:*(a::Dual, b::Number) = Dual(a.x * b, a.der * b)
Base.:*(a::Number, b::Dual) = b * a

Base.:/(a::Dual, b::Dual) = Dual(a.x / b.x, (a.der * b.x - a.x * b.der) / (b.x ^ 2))
Base.:/(a::Dual, b::Number) = a * inv(b)
Base.:/(a::Number, b::Dual) = Dual(a / b.x, ( -a * b.der ) / ( b.x ^ 2 ) ) #

Base.:^(a::Dual, b::Integer) = Base.power_by_squaring(a, b)
Base.:^(a::Number, b::Dual) = Dual(Float64(a^b.x), a^b.x * log(a))

Base.:≈(a::Dual, b::Dual) = a.x ≈ b.x && a.der ≈ b.der

derivative(f, x) = f(Dual(x, one(x))).der

# primatives
Base.sin(a::Dual) = Dual(sin(a.x), a.der * cos(a.x))
Base.cos(a::Dual) = Dual(cos(a.x), -a.der * sin(a.x))
Base.tan(a::Dual) = sin(a)/cos(a)

Base.csc(a::Dual) = Dual(-csc(a.x)*cot(a.x), -a.der * sin(a.x))
Base.sec(a::Dual) = Dual(tan(a.x), -a.der * sin(a.x))
Base.cot(a::Dual) = 1/tan(a)

Base.exp(a::Dual) = Dual(exp(a.x), exp(a.x) * a.der)
Base.log(a::Dual) = Dual(log(a.x), a.der/a.x)

# used with ∇ and N-dimensional scalar valued functions
# struct MultiDual{N, T}
#     x::T
#     ders::SVector{N, T}
# end

# Base.:+(a::MultiDual, b::MultiDual) = MultiDual(a.x + b.x, a.ders .+ b.ders)
# Base.:*(a::MultiDual, b::MultiDual) = MultiDual(a.x * b.x, a.ders .* b.x .+ a.x .* b.ders)

# ∇(f, x) = f(MultiDual(x, ones(length(x)))).ders

# struct MultiDual2{M, N, T}
#     x::SVector{M, T}
#     J::SMatrix{M, N, T}
# end


text:
using CompilerPluginTools: rm_code_coverage_effect
using IBMQClient.Schema
using YaoCompiler
using YaoCompiler.Intrinsics
using YaoTargetQobj
using Test

@device function basic()
    1 => X
    2 => Y
    3 => Z
    4 => H
    4 => T
    4 => S
    4 => Rx(1.2)
    4 => Ry(1.2)
    4 => Rz(1.2)
    4 => shift(1.2)
end

@testset "basic" begin
    ci, type = @yao_code_typed(basic())[1]
    ci = rm_code_coverage_effect(ci)
    @test measure_ssa_uses(ci) == Set{Int}()
    mi = MemoryInfo(ci)
    display(mi)
    @test mi.qubits == Dict(1=>0, 2=>1, 3=>2, 4=>3)
    @test mi.n_qubits == 4

    @test code_qobj(typeof(basic())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(basic), Tuple{}}",
            "n_qubits" => 4,
        ),
        config = ExpConfig(;
            memory_slots = 0,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "x",
                qubits = [0],
            ),
            Gate(;
                name = "y",
                qubits = [1],
            ),
            Gate(;
                name = "z",
                qubits = [2],
            ),
            Gate(;
                name = "h",
                qubits = [3],
            ),
            Gate(;
                name = "t",
                qubits = [3],
            ),
            Gate(;
                name = "s",
                qubits = [3],
            ),
            Gate(;
                name = "rx",
                qubits = [3],
                params = [1.2],
            ),
            Gate(;
                name = "ry",
                qubits = [3],
                params = [1.2],
            ),
            Gate(;
                name = "rz",
                qubits = [3],
                params = [1.2],
            ),
            Gate(;
                name = "u1",
                qubits = [3],
                params = [1.2],
            ),
        ],
    )
end

@device function just_measure()
    1 => H
    c = @measure 1
    return c
end

@testset "just_measure" begin
    ci, type = @yao_code_typed(just_measure())[1]
    ci = rm_code_coverage_effect(ci)
    @test measure_ssa_uses(ci) == Set{Int}()
    mi = MemoryInfo(ci)
    display(mi)
    @test isempty(mi.registers)

    @test code_qobj(typeof(just_measure())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(just_measure), Tuple{}}",
            "n_qubits" => 1,
        ),
        config = ExpConfig(;
            memory_slots = 1,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Measure(;
                qubits = [0],
                memory = [0],
            ),
        ],
    )
end

@device function reset_qubit()
    1 => H
    c = @measure 1
    if c == 1
        1 => X
    end
    return
end

@testset "reset_qubit" begin
    ci, type = @yao_code_typed(reset_qubit())[1]
    ci = rm_code_coverage_effect(ci)
    @test measure_ssa_uses(ci) == Set(2)
    @test measure_ssa_returns(ci) == Set()

    mi = MemoryInfo(ci)
    display(mi)
    @test mi.registers == Dict(2=>[0])
    @test mi.n_qubits == 1
    @test mi.qubits == Dict(1=>0)
    @test code_qobj(typeof(reset_qubit())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(reset_qubit), Tuple{}}",
            "n_qubits" => 1,
        ),
        config = ExpConfig(;
            memory_slots = 0,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Reset(;
                qubits = [0],
            ),
        ],
    )
end

@device function reset_qubit_return()
    1 => H
    c = @measure 1
    if c == 1
        1 => X
    end
    return c
end

@testset "reset_qubit_return" begin
    ci, type = @yao_code_typed(reset_qubit_return())[1]
    ci = rm_code_coverage_effect(ci)
    @test measure_ssa_uses(ci) == Set(2)
    # @test measure_ssa_returns(ci) == Set(2)
    @test code_qobj(typeof(reset_qubit_return())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(reset_qubit_return), Tuple{}}",
            "n_qubits" => 1,
        ),
        config = ExpConfig(;
            memory_slots = 1,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Measure(;
                qubits = [0],
                memory = [0],
                register = [0],
            ),
            Gate(;
                name = "x",
                qubits = [0],
                conditional = 0,
            ),
        ],
    )
end

@device function reset_with_return()
    1 => H
    m1 = @measure 2
    m2 = @measure 1
    if m1 == 1
        1 => X
    end
    return m1, m2
end

@testset "reset_with_return" begin
    ci, type = @yao_code_typed(reset_with_return())[1]
    ci = rm_code_coverage_effect(ci)
    @test measure_ssa_uses(ci) == Set(2)
    # @test measure_ssa_returns(ci) == Set([2, 3])
    mi = MemoryInfo(ci)
    @test mi.n_qubits == 2
    @test mi.registers == Dict(2=>[0])
    @test mi.qubits == Dict(1=>0, 2=>1)
    @test code_qobj(typeof(reset_with_return())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(reset_with_return), Tuple{}}",
            "n_qubits" => 2,
        ),
        config = ExpConfig(;
            memory_slots = 2,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Measure(;
                qubits = [1],
                memory = [0],
                register = [0],
            ),
            Measure(;
                qubits = [0],
                memory = [1],
            ),
            Gate(;
                name = "x",
                qubits = [0],
                conditional = 0,
            ),
        ],
    )
end

@device function multi_qubit_measure()
    1 => H
    3 => X
    c = @measure (1, 2)
    return (c=c)
end

@testset "multi_qubit_measure" begin
    ci, type = @yao_code_typed(multi_qubit_measure())[1]
    ci = rm_code_coverage_effect(ci)
    mi = MemoryInfo(ci)
    @test mi.n_qubits == 3
    @test isempty(mi.registers)
    @test mi.qubits == Dict(1=>0, 2=>2, 3=>1)
    @test code_qobj(typeof(multi_qubit_measure())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(multi_qubit_measure), Tuple{}}",
            "n_qubits" => 3,
        ),
        config = ExpConfig(;
            memory_slots = 2,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Gate(;
                name = "x",
                qubits = [1],
            ),
            Measure(;
                qubits = [0, 2],
                memory = [0, 1],
            ),
        ],
    )
end

@device function multi_if_measure()
    1 => H
    3 => X
    c = @measure 1
    if c == 1
        1 => Y
        2 => Z
        3 => H
    end
    return
end

@testset "multi if stmt" begin
    @test code_qobj(typeof(multi_if_measure())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(multi_if_measure), Tuple{}}",
            "n_qubits" => 3,
        ),
        config = ExpConfig(;
            memory_slots = 1,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Gate(;
                name = "x",
                qubits = [1],
            ),
            Measure(;
                qubits = [0],
                memory = [0],
                register = [0],
            ),
            Gate(;
                name = "y",
                qubits = [0],
                conditional = 0,
            ),
            Gate(;
                name = "z",
                qubits = [2],
                conditional = 0,
            ),
            Gate(;
                name = "h",
                qubits = [1],
                conditional = 0,
            ),
        ],
    )
end

@device function conditional_ctrl()
    1 => H
    3 => X
    c = @measure 1
    if c == 1
        @ctrl 2 1 => Y
    end
    return
end

@testset "conditional ctrl" begin
    @test code_qobj(typeof(conditional_ctrl())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(conditional_ctrl), Tuple{}}",
            "n_qubits" => 3,
        ),
        config = ExpConfig(;
            memory_slots = 1,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Gate(;
                name = "x",
                qubits = [1],
            ),
            Measure(;
                qubits = [0],
                memory = [0],
                register = [0],
            ),
            Gate(;
                name = "cy",
                qubits = [0, 2],
                conditional = 0,
            ),
        ],
    )
end

@device function test_barrier()
    1 => H
    3 => X
    @barrier 1,3
    c = @measure 1
    if c == 1
        @ctrl 2 1 => Y
    end
    return
end

@testset "test_barrier" begin
    @test code_qobj(typeof(test_barrier())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(test_barrier), Tuple{}}",
            "n_qubits" => 3,
        ),
        config = ExpConfig(;
            memory_slots = 1,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "h",
                qubits = [0],
            ),
            Gate(;
                name = "x",
                qubits = [1],
            ),
            Barrier(;
                qubits = [0, 1],
            ),
            Measure(;
                qubits = [0],
                memory = [0],
                register = [0],
            ),
            Gate(;
                name = "cy",
                qubits = [0, 2],
                conditional = 0,
            ),
        ],
    )
end

@device function test_ctrl()
    @ctrl 4 1 => H
    @ctrl 4 3 => X
    @ctrl 4 3 => Y
    @ctrl 4 3 => Z
    @ctrl (4,5) 3 => X
    @ctrl 4 3 => Rx(2.1)
    @ctrl 4 3 => Ry(2.1)
    @ctrl 4 3 => Rz(2.1)
    @ctrl 4 3 => shift(2.1)
    @ctrl 4 3 => Intrinsics.UGate(1.2, 2.2, 3.3)
    return
end

@testset "test_ctrl" begin
    @test code_qobj(typeof(test_ctrl())) == Experiment(;
        header = Dict{String, Any}(
            "description" => "Operation{typeof(test_ctrl), Tuple{}}",
            "n_qubits" => 4,
        ),
        config = ExpConfig(;
            memory_slots = 0,
            seed = 1,
            max_credits = 3,
        ),
        instructions = [
            Gate(;
                name = "ch",
                qubits = [0, 1],
            ),
            Gate(;
                name = "cx",
                qubits = [2, 1],
            ),
            Gate(;
                name = "cy",
                qubits = [2, 1],
            ),
            Gate(;
                name = "cz",
                qubits = [2, 1],
            ),
            Gate(;
                name = "ccx",
                qubits = [2, 1, 3],
            ),
            Gate(;
                name = "crx",
                qubits = [2, 1],
                params = [2.1],
            ),
            Gate(;
                name = "cry",
                qubits = [2, 1],
                params = [2.1],
            ),
            Gate(;
                name = "crz",
                qubits = [2, 1],
                params = [2.1],
            ),
            Gate(;
                name = "cu1",
                qubits = [2, 1],
                params = [2.1],
            ),
            Gate(;
                name = "cu3",
                qubits = [2, 1],
                params = [1.2, 2.2, 3.3],
            ),
        ],
    )
end


text:
module PetriNet

using Hecke, Francy

mutable struct Data
  name::Symbol
  data::Dict{Symbol, Any}
  id::Symbol
  function Data(n::Symbol)
    r = new()
    r.name = n
    r.data = Dict{Symbol, Any}()
    return r
  end
  function Data(n::Symbol, d::Pair{Symbol, <:Any}...)
    r = new()
    r.name = n
    r.data = Dict(d...)
    return r
  end
end

mutable struct Transition
  name::Symbol
  action::Vector{Function}
  possible::Vector{Function}
  id::Symbol
  In::Vector{Symbol}
  Out::Vector{Symbol}
  Read::Vector{Symbol}
  function Transition(n::Symbol)
    r = new()
    r.name = n
    return r
  end
end

mutable struct Petri
  d::Dict{Symbol, Data}
  t::Dict{Symbol, Transition}
  change::Set{Symbol}
  function Petri()
    r = new()
    r.d = Dict{Symbol, Data}()
    r.t = Dict{Symbol, Transition}()
    return r
  end
end

function Base.push!(p::Petri, d::Data)
  p.d[d.name] = d
end
function Base.push!(p::Petri, t::Transition)
  p.t[t.name] = t
end

function consume(d::Data, s::Symbol)
  a = d.data[s]
  if isa(a, Array)
    return pop!(a)
  else
    delete!(d.data, s)
    return a
  end
end

function consume(P::Petri, s::Symbol, t::Symbol = s)
  push!(P.change, s)
  return consume(P.d[s], t)
end

function Base.length(d::Data, s::Symbol)
  return length(d.data[s])
end

function read(d::Data, s::Symbol)
  return d.data[s]
end

function put(d::Data, s::Symbol, a::Any)
  if !haskey(d.data, s)
    d.data[s] = a
  else
    push!(d.data[s], a)
  end
end

function put(P::Petri, s::Symbol, t::Symbol, a::Any)
  push!(P.change, s)
  return put(P.d[s], t, a)
end

function put(P::Petri, s::Symbol, a::Any)
  return put(P, s, s, a)
end

function _try()  #parallel(?) sum(1:10)
  P = Petri()
  push!(P, Data(:p, :p =>collect(1:10)))
  push!(P, Data(:r, :r =>[]))
  push!(P, Data(:s, :s =>[]))
  push!(P, Data(:c1, :c => 1))
  push!(P, Data(:c2))

  T = Transition(:T)
  T.action = [() -> put(P, :s, consume(P, :p)),
              () -> put(P, :c2, :c, consume(P, :c1, :c))]
  T.possible = [()-> length(P.d[:p], :p)>0,
                () -> haskey(P.d[:c1].data, :c)]
  T.In = [:p, :c1]
  T.Out = [:s, :c2]
  push!(P, T)

  T = Transition(:S)
  T.action = [() -> put(P, :r, consume(P, :p)),
              () -> put(P, :c1, :c, consume(P, :c2, :c))]
  T.possible = [()-> length(P.d[:p].data[:p])>0,
                () -> haskey(P.d[:c2].data, :c)]
  T.In = [:p, :c2]
  T.Out = [:r, :c1]
  push!(P, T)

  T = Transition(:+)
  T.action = [ () -> put(P, :p, consume(P, :s) + consume(P, :r))]
  T.possible = [ () -> length(P.d[:r].data[:r]) > 0 && length(P.d[:s].data[:s]) > 0]
  T.In = [:s, :r]
  T.Out = [:p]
  push!(P, T)

  return P
end

function fire(P::Petri)
  P.change = Set{Symbol}()
  p = [k for k = keys(P.t) if all(i->i(), P.t[k].possible)]
  for k = p
    [x() for x = P.t[k].action]
  end
#  println(stderr, p)
#  println(stderr, P.d)
end

function fire()
  global pe, M, c, pe
  fire(pe)
  for (k,v) = M
    delete!(c.c[:messages], v.m[:id])
  end

  for d = values(pe.d)
    if false && haskey(d.data, d.name)
      m = message("$(d.name): $(d.data[d.name])")
      push!(c, m)
      M[d.name] = m
    else
      m = message("$(d.name): $(d.data)")
      push!(c, m)
      M[d.name] = m
    end
  end
  display(Francy.FrancyMimeString, c)
end

function fire(t::Symbol)
  global pe, M, c, pe, T, D
  pe.change = Set{Symbol}()
  if all(i->i(), pe.t[t].possible)
    for i = pe.t[t].action
      i()
    end
  end
  for (k,v) = M
    delete!(c.c[:messages], v.m[:id])
  end

  for d = values(pe.d)
    if false && haskey(d.data, d.name)
      m = message("$(d.name): $(d.data[d.name])")
      push!(c, m)
      M[d.name] = m
    else
      m = message("$(d.name): $(d.data)")
      push!(c, m)
      M[d.name] = m
    end
  end
  for (k,v) in T
    if all(i -> i(), pe.t[k].possible)
      v.s[:color] = :green
    else
      v.s[:color] = :blue
    end
  end
  for (d, n) in D
    if haskey(pe.d[d].data, d)
      da = pe.d[d].data[d]
      if length(da) == 0
        n.s[:title] = "$(string(d)):[]"
      else
        n.s[:title] = "$(string(d)):[..$(da[end])]"
      end
    else
      k = keys(pe.d[d].data)
      if length(k) == 0
        n.s[:title] = "$(string(d)):"
      else
        n.s[:title] = "$(string(d)):$(first(k))"
      end
    end
  end
  display(Francy.FrancyMimeString, c)
end

function petri(P::Petri)
  global c,A, n, g, pe, M, T, D
  pe = P
  c = canvas("Petri")
  g = graph()
  D = Dict{Symbol, Francy.Node}()
  T = Dict{Symbol, Francy.Node}()
  M = Dict{Symbol, Francy.Message}()
  for d = keys(P.d)
    n = D[d] = node(string(d), :circle)
    if haskey(P.d[d].data, d)
      da = P.d[d].data[d]
      if length(da) == 0
        n.s[:title] = "$(string(d)):[]"
      else
        n.s[:title] = "$(string(d)):..$(da[end])]"
      end
    else
      k = keys(P.d[d].data)
      if length(k) == 0
        n.s[:title] = "$(string(d)):"
      else
        n.s[:title] = "$(string(d)):$(first(k))"
      end
    end
    push!(g, n)
  end
  for d = keys(P.t)
    n = T[d] = node(string(d), :square)
    push!(n, callback("PetriNet.fire", ":($d)"))
    if all(i->i(), P.t[d].possible)
      n.s[:color] = :green
    else
      n.s[:color] = :blue
    end
    push!(g, n)
  end

  for t = values(P.t)
    for i = t.In
      push!(g, link(D[i], T[t.name], length = 5))
    end
    for i = t.Out
      push!(g, link(T[t.name], D[i], length = 5))
    end
  end
  for d = values(P.d)
    if false && haskey(d.data, d.name)
      m = message("$(d.name): $(d.data[d.name])")
      push!(c, m)
      M[d.name] = m
    else
      m = message("$(d.name): $(d.data)")
      push!(c, m)
      M[d.name] = m
    end
  end
#  push!(n[2], callback("PetriNet.trigger", "(:"*n[2].s[:id]*", 1)"))
  push!(c, g)
  return c
end

node_index(s::Symbol) = findfirst(x -> x.s[:id] == String(s), n)

function in_nodes(s::Symbol)
  global g
  res = []
  for i = values(g.g[:links])
    if i[:target] == String(s)
      push!(res, Symbol(i[:source]))
    end
  end
  return res
end

function out_nodes(s::Symbol)
  global g
  res = []
  for i = values(g.g[:links])
    if i[:source] == String(s)
      push!(res, Symbol(i[:target]))
    end
  end
  return res
end

function trigger(t::Tuple)
  global c, n, g

  for i in in_nodes(t[1])
    j = node_index(i)
    n[j].s[:title] = string(Char(Int(n[j].s[:title][1])-1))
  end
  for i in out_nodes(t[1])
    j = node_index(i)
    n[j].s[:title] = string(Char(Int(n[j].s[:title][1])+1))
  end

  display(Francy.FrancyMimeString, c)
end

end


text:
abstract type MessagePassing <: GraphNet end

"""
    message(mp::MessagePassing, x_i, x_j, e_ij)

Message function for the message-passing scheme,
returning the message from node `j` to node `i` .
In the message-passing scheme. the incoming messages
from the neighborhood of `i` will later be aggregated
in order to [`update`](@ref) the features of node `i`.

By default, the function returns `x_j`.
Layers subtyping [`MessagePassing`](@ref) should
specialize this method with custom behavior.

# Arguments

- `mp`: message-passing layer.
- `x_i`: the features of node `i`.
- `x_j`: the features of the nighbor `j` of node `i`.
- `e_ij`: the features of edge (`i`, `j`).

See also [`update`](@ref).
"""
@inline message(mp::MessagePassing, x_i, x_j, e_ij) = x_j
@inline message(mp::MessagePassing, i::Integer, j::Integer, x_i, x_j, e_ij) = x_j

"""
    update(mp::MessagePassing, m, x)

Update function for the message-passing scheme,
returning a new set of node features `x′` based on old
features `x` and the incoming message from the neighborhood
aggregation `m`.

By default, the function returns `m`.
Layers subtyping [`MessagePassing`](@ref) should
specialize this method with custom behavior.

# Arguments

- `mp`: message-passing layer.
- `m`: the aggregated edge messages from the [`message`](@ref) function.
- `x`: the node features to be updated.

See also [`message`](@ref).
"""
@inline update(mp::MessagePassing, m, x) = m
@inline update(mp::MessagePassing, i::Integer, m, x) = m

@inline apply_batch_message(mp::MessagePassing, i, js, edge_idx, E::AbstractMatrix, X::AbstractMatrix, u) =
    mapreduce(j -> GeometricFlux.message(mp, _view(X, i), _view(X, j), _view(E, edge_idx[(i,j)])), hcat, js)

@inline update_batch_vertex(mp::MessagePassing, M::AbstractMatrix, X::AbstractMatrix, u) =
    mapreduce(i -> GeometricFlux.update(mp, _view(M, i), _view(X, i)), hcat, 1:size(X,2))

function propagate(mp::MessagePassing, fg::FeaturedGraph, aggr=+)
    E, X = propagate(mp, adjacency_list(fg), fg.ef, fg.nf, aggr)
    FeaturedGraph(fg, nf=X, ef=E, gf=Fill(0.f0, 0))
end

function propagate(mp::MessagePassing, adj::AbstractVector{S}, E::R, X::Q, aggr) where {S<:AbstractVector,R,Q}
    E, X, u = propagate(mp, adj, E, X, Fill(0.f0, 0), aggr, nothing, nothing)
    E, X
end


text:
using Documenter, MCPhyloTree

makedocs(
    modules = MCPhyloTree,
    sitename="MCPhyloTree.jl",
    format = Documenter.HTML(),
    pages = [
        "Index" => "index.md",
        "Basics" => "basics.md",
        "Building" => "building.md",
        "Distance" => "distance.md",
        "Moves" => "moves.md",
        "Plotting" => "plotting.md"
        ]
   )

deploydocs(
    repo = "github.com/erathorn/MCPhyloTree.jl.git",
    devbranch="main",
    versions = ["stable" => "v^", "v#.#.#"]
     )


text:
@testset "969.pancake-sorting.jl" begin
    @test issorted(pancake_sort!([3, 2, 4, 1], pancake_sort([3, 2, 4, 1])))
    @test issorted(pancake_sort!([1, 2, 3], pancake_sort([1, 2, 3])))
end


text:
using LinearAlgebra
import ClimaCore:
    Fields,
    Domains,
    Topologies,
    Meshes,
    DataLayouts,
    Operators,
    Geometry,
    Spaces

using OrdinaryDiffEq: ODEProblem, solve, SSPRK33

import Logging
import TerminalLoggers
Logging.global_logger(TerminalLoggers.TerminalLogger())

const FT = Float64

a = Geometry.ZPoint(FT(-20.0))
b = Geometry.ZPoint(FT(20.0))
n = 64
α = FT(0.1)

function heaviside(pt)
    0.5 * (sign(pt.z) + 1)
end

domain = Domains.IntervalDomain(a, b, boundary_names = (:left, :right))
mesh = Meshes.IntervalMesh(domain, nelems = n)

cs = Spaces.CenterFiniteDifferenceSpace(mesh)
fs = Spaces.FaceFiniteDifferenceSpace(cs)

V = Geometry.WVector.(ones(FT, fs))
θ = heaviside.(Fields.coordinate_field(cs))

# Solve advection Equation: ∂θ/dt = -∂(vθ)

# upwinding
function tendency1!(dθ, θ, _, t)
    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(sin(a.z - t)),
        right = Operators.SetValue(sin(b.z - t)),
    )
    ∂ = Operators.DivergenceF2C()

    return @. dθ = -∂(UB(V, θ))
end
function tendency2!(dθ, θ, _, t)
    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    UB = Operators.UpwindBiasedProductC2F(
        left = Operators.SetValue(sin(a.z - t)),
        right = Operators.SetValue(sin(b.z - t)),
    )
    ∂ = Operators.DivergenceF2C()
    return @. dθ = -∂(UB(V, θ)) + fcc(V, θ)
end
# use the advection operator
function tendency3!(dθ, θ, _, t)

    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(sin(-t)),
        right = Operators.Extrapolate(),
    )
    return @. dθ = -A(V, θ)
end
# use the advection operator
function tendency4!(dθ, θ, _, t)

    fcc = Operators.FluxCorrectionC2C(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    fcf = Operators.FluxCorrectionF2F(
        left = Operators.Extrapolate(),
        right = Operators.Extrapolate(),
    )
    A = Operators.AdvectionC2C(
        left = Operators.SetValue(sin(-t)),
        right = Operators.Extrapolate(),
    )
    return @. dθ = -A(V, θ) + fcc(V, θ)
end

# use the advection operator

@show tendency1!(similar(θ), θ, nothing, 0.0)
# Solve the ODE operator
Δt = 0.001
prob1 = ODEProblem(tendency1!, θ, (0.0, 5.0))
prob2 = ODEProblem(tendency2!, θ, (0.0, 5.0))
prob3 = ODEProblem(tendency3!, θ, (0.0, 5.0))
prob4 = ODEProblem(tendency4!, θ, (0.0, 5.0))
sol1 = solve(
    prob1,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol2 = solve(
    prob2,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol3 = solve(
    prob3,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);
sol4 = solve(
    prob4,
    SSPRK33(),
    dt = Δt,
    saveat = 10 * Δt,
    progress = true,
    progress_message = (dt, u, p, t) -> t,
);

ENV["GKSwstype"] = "nul"
import ClimaCorePlots, Plots
Plots.GRBackend()

dirname = "advect_step_function"
path = joinpath(@__DIR__, "output", dirname)
mkpath(path)

anim = Plots.@animate for u in sol1.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "UBP_advect_step_function.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol1.u[end], xlim = (-1, 1)),
    joinpath(path, "sol1_advect_step_function_end.png"),
)

anim = Plots.@animate for u in sol2.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "UBP_advect_step_function_fc.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol2.u[end], xlim = (-1, 1)),
    joinpath(path, "sol2_advect_step_function_end.png"),
)

anim = Plots.@animate for u in sol3.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "C2C_advect_step_function.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol3.u[end], xlim = (-1, 1)),
    joinpath(path, "sol3_advect_step_function_end.png"),
)

anim = Plots.@animate for u in sol4.u
    Plots.plot(u, xlim = (-1, 1))
end
Plots.mp4(anim, joinpath(path, "C2C_advect_step_function_fc.mp4"), fps = 10)
Plots.png(
    Plots.plot(sol4.u[end], xlim = (-1, 1)),
    joinpath(path, "sol4_advect_step_function_end.png"),
)

p = Plots.plot(sol1.u[end], xlim = (-1, 1), label = "UBP")
p = Plots.plot!(sol2.u[end], xlim = (-1, 1), label = "UBP_FC")
p = Plots.plot!(sol3.u[end], xlim = (-1, 1), label = "C2C")
p = Plots.plot!(sol4.u[end], xlim = (-1, 1), label = "C2C_FC")
Plots.png(p, joinpath(path, "all_advect_step_function_end.png"))

function linkfig(figpath, alt = "")
    # buildkite-agent upload figpath
    # link figure in logs if we are running on CI
    if get(ENV, "BUILDKITE", "") == "true"
        artifact_url = "artifact://$figpath"
        print("\033]1338;url='$(artifact_url)';alt='$(alt)'\a\n")
    end
end

linkfig(
    "examples/column/output/$(dirname)/advect_step_function_end.png",
    "Advect End Simulation",
)

linkfig(
    relpath(
        joinpath(path, "advect_step_function_end.png"),
        joinpath(@__DIR__, "../.."),
    ),
    "Advect End Simulation",
)


text:
module Progress
using Logging

using Atom: msg

const progs = Set()

"""
    JunoProgressLogger()

Logger that only handles messages with `progress` argument set and passes all
others through to the `global_logger`.

If the `progress` argument is set, then a progress bar will be shown in the Juno frontend.

Possible keyword arguments to logging messages consumed by this logger are:
  - `progress`:
    - `0 <= progress < 1`: create or update progress bar
    - `progress == nothing || progress = NaN`: set progress bar to indeterminate progress
    - `progress > 1 || progress == "done"`: destroy progress bar
  - `right_text`: Shown instead of a approximation of remaining time.
  - `message`: Shown in a tooltip.
  - `_id`: The progress bars' ID. Should be set (to a symbol or string) in most cases.
"""
struct JunoProgressLogger <: AbstractLogger end

function Logging.handle_message(j::JunoProgressLogger, level, message, _module,
                                group, id, file, line; kwargs...)
  if haskey(kwargs, :progress)
    if !(id in progs)
      push!(progs, id)
      msg("progress", "add", id)
    end
    msg("progress", "leftText", id, message)

    prog = kwargs[:progress]
    if prog isa Real
      if 0 <= prog < 1
        msg("progress", "progress", id, prog)
      elseif prog >= 1
        msg("progress", "delete", id)
        delete!(progs, id)
      end
    elseif prog isa Nothing
      msg("progress", "progress", id, NaN)
    elseif prog isa AbstractString && prog == "done"
      msg("progress", "delete", id)
      delete!(progs, id)
    else
      msg("progress", "progress", id, NaN)
    end

    if haskey(kwargs, :right_text)
      msg("progress", "rightText", id, kwargs[:right_text])
    end

    if haskey(kwargs, :message)
      msg("progress", "message", id, kwargs[:message])
    end
  else
    previous_logger = Logging.global_logger()
    if (Base.invokelatest(Logging.min_enabled_level, previous_logger) <= Logging.LogLevel(level) ||
        Base.CoreLogging.env_override_minlevel(group, _module)) &&
        Base.invokelatest(Logging.shouldlog, previous_logger, level, _module, group, id)
      Logging.handle_message(previous_logger, level, message, _module,
                             group, id, file, line; kwargs...)
    end
  end
  return nothing
end

Logging.shouldlog(::JunoProgressLogger, level, _module, group, id) = true

Logging.catch_exceptions(::JunoProgressLogger) = true

function Logging.min_enabled_level(j::JunoProgressLogger)
  min(Base.invokelatest(Logging.min_enabled_level, Logging.global_logger()), Logging.LogLevel(-1))
end

end # module


text:
# -------------------------------------------------------------------
# AST trees

function add_ast(adjlist, names, depthdict, depthlists, nodetypes, ex::Expr, parent_idx)
    idx = length(names)+1
    iscall = ex.head == :call
    push!(names, iscall ? string(ex.args[1]) : string(ex.head))
    push!(nodetypes, iscall ? :call : :expr)
    l = Int[]
    push!(adjlist, l)

    depth = parent_idx==0 ? 1 : depthdict[parent_idx] + 1
    depthdict[idx] = depth
    while length(depthlists) < depth
        push!(depthlists, Int[])
    end
    push!(depthlists[depth], idx)

    for arg in (iscall ? ex.args[2:end] : ex.args)
        if isa(arg, LineNumberNode)
            continue
        end
        push!(l, add_ast(adjlist, names, depthdict, depthlists, nodetypes, arg, idx))
    end
    idx
end

function add_ast(adjlist, names, depthdict, depthlists, nodetypes, x, parent_idx)
    push!(names, string(x))
    push!(nodetypes, :leaf)
    push!(adjlist, Int[])
    idx = length(names)

    depth = parent_idx==0 ? 1 : depthdict[parent_idx] + 1
    depthdict[idx] = depth
    while length(depthlists) < depth
        push!(depthlists, Int[])
    end
    push!(depthlists[depth], idx)

    idx
end

@recipe function f(ex::Expr)
    names = String[]
    adjlist = Vector{Int}[]
    depthdict = Dict{Int,Int}()
    depthlists = Vector{Int}[]
    nodetypes = Symbol[]
    add_ast(adjlist, names, depthdict, depthlists, nodetypes, ex, 0)
    names := names
    # method := :tree
    method := :buchheim
    root --> :top

    markercolor --> Symbol[(nt == :call ? :pink : nt == :leaf ? :white : :lightgreen) for nt in nodetypes]

    # # compute the y-values from the depthdict dict
    # n = length(depthlists)-1
    # layers = Float64[(depthdict[i]-1)/n for i=1:length(names)]
    # # add_noise --> false
    #
    # positions = zeros(length(names))
    # for (depth, lst) in enumerate(depthlists)
    #     n = length(lst)
    #     pos = n > 1 ? linspace(0, 1, n) : [0.5]
    #     for (i, idx) in enumerate(lst)
    #         positions[idx] = pos[i]
    #     end
    # end
    #
    # layout_kw := Dict{Symbol,Any}(:layers => layers, :add_noise => false, :positions => positions)

    GraphPlot(get_source_destiny_weight(adjlist))
end

# -------------------------------------------------------------------
# Type trees

function add_subs!(nodes, source, destiny, ::Type{T}, supidx) where T
    for sub in subtypes(T)
        push!(nodes, sub)
        subidx = length(nodes)
        push!(source, supidx)
        push!(destiny, subidx)
        add_subs!(nodes, source, destiny, sub, subidx)
    end
end

# recursively build a graph of subtypes of T
@recipe function f(::Type{T}; namefunc = node-> isa(node, UnionAll) ? split(string(node), '.')[end] : node.name.name) where T
    # get the supertypes
    sups = Any[T]
    sup = T
    while sup != Any
        sup = supertype(sup)
        pushfirst!(sups,sup)
    end

    # add the subtypes
    n = length(sups)
    nodes = copy(sups)
    source, destiny = collect(1:n-1), collect(2:n)
    add_subs!(nodes, source, destiny, T, n)

    # set up the graphplot
    names := map(namefunc, nodes)
    method --> :buchheim
    root --> :top
    GraphPlot((source, destiny))
end


text:
using DataFrames, CSV # For loading data
using Plots # For plotting the total error in training process
using Distributions # For generating weights

argA = 0 # An argument value in binary
argB = 1 # Another arugment value in binary
bias = 1 # The bias term is used to help the neural network fit best for the input data

inputNum = 5  # Include bias term
hiddenNum = 4 # Include bias term
outputNum = 1
momentum = 0.9 # A term used to increase the neural network's tranining speed
learningRate = 0.2

###  Activation Function (Sigmoid)  ###
function sigmoid(x::Float64)
    result = (argB - argA) / (1 + exp(-x)) + argA
    return result
end

###  Activation Function (ReLU)  ###
function ReLU(x::Float64)
    result = max(0.0, x)
    return result
end

###  Training Function  ###
function train(input::Array{Float64,1}, expectedOutput::Float64)
    trainOutput = outputFor(input) # forward propagation
    updateWeight(trainOutput, expectedOutput) # backpropagation
    errorRate = 0.5 * (trainOutput - expectedOutput)^2
    return errorRate
end

###  Forward Propagation Function  ###
function outputFor(input::Array{Float64,1})
    for i = 1:inputNum - 1
        inputNeuron[i] = input[i]
    end

    inputNeuron[inputNum] = bias
    hiddenNeuron[hiddenNum] = bias

    for i = 1:hiddenNum - 1
        for j = 1:inputNum
           hiddenNeuron[i] += inputWeights[j, i] * inputNeuron[j]
        end
        hiddenNeuron[i] = sigmoid(hiddenNeuron[i])
    end

    for k = 1:hiddenNum
        outputNeuron[1] += hiddenWeights[k] * hiddenNeuron[k]
    end
    outputNeuron[1] = sigmoid(outputNeuron[1])
    return outputNeuron[1]
end

###  Backpropagation Function  ###
function updateWeight(trainOutput::Float64, expectedOutput::Float64)
    deltaOutput[1] = trainOutput * (1 - trainOutput) * (expectedOutput - trainOutput)

    for k = 1:hiddenNum
        hiddenWeights[k] +=  (momentum * hiddenWeightsDelta[k]) + (learningRate * deltaOutput[1] * hiddenNeuron[k]);
        hiddenWeightsDelta[k] = (momentum * hiddenWeightsDelta[k]) + (learningRate * deltaOutput[1] * hiddenNeuron[k]);
    end

    for k = 1:hiddenNum - 1
        sum = 0.0
        deltaHidden[k] = hiddenNeuron[k] * (1 - hiddenNeuron[k])
        sum += deltaOutput[1] * hiddenWeights[k]
        deltaHidden[k] *= sum
    end

    for j = 1:hiddenNum - 1
        for i = 1:inputNum
            inputWeights[i, j] += (momentum * inputWeightsDelta[i, j]) + (learningRate * deltaHidden[j] * inputNeuron[i])
            inputWeightsDelta[i, j] = (momentum * inputWeightsDelta[i, j]) + (learningRate * deltaHidden[j] * inputNeuron[i])
        end
    end
end

###  Initilize Neurons and Weights  ###
inputNeuron = zeros(inputNum)
hiddenNeuron = zeros(hiddenNum)
outputNeuron = zeros(outputNum)

inputWeights = rand(Uniform(-0.5, 0.5), (inputNum, hiddenNum - 1))
hiddenWeights = rand(Uniform(-0.5, 0.5), (hiddenNum, outputNum))

### Error Signal for Neurons and Weights (Used in Backpropagation)  ###
inputWeightsDelta = zeros(inputNum, hiddenNum - 1)
hiddenWeightsDelta = zeros(hiddenNum, outputNum)

deltaHidden = zeros(hiddenNum)
deltaOutput = zeros(outputNum)

epoch = 0
errorRate = 0.0
acceptError = 0.05 # The threshold value of error rate

###  Read Data from the CSV File  ###
df = DataFrame(CSV.File("data/germandata.csv"))
head(df)

inputATrain = df[1:600,1]  # Credit history
inputBTrain = df[1:600,2]  # Other installment plans
inputCTrain = df[1:600,3]  # Property
inputDTrain = df[1:600,4]  # Telephone
expectedOutput = df[1:600,5] # Expected training output

epochVector = zeros(0) # Store epoch number for plotting
errorRateVector = zeros(0) # Store error rate in each epoch for plotting

###  Training Process  ###

while(epoch == 0 || errorRate > acceptError)
    errorRate = 0.0
    for index in 1:600
        input = convert(Array{Float64,1}, [inputATrain[index], inputBTrain[index], inputCTrain[index], inputDTrain[index]])
        output = convert(Float64, expectedOutput[index])
        errorRate += train(input, output)
        end
    println("Error at epoch $epoch is $errorRate")
    append!(epochVector, convert(Float64,epoch))
    append!(errorRateVector, errorRate)
    epoch += 1
end

### Plot the Training Process  ###
plot(epochVector, errorRateVector, title = "Total Squared Error in Training Process", lw = 3, legend = false)
xlabel!("Epoch")
ylabel!("Error")

###  Testing Process  ###
inputATest = df[601:1000,1]  # Credit history
inputBTest = df[601:1000,2]  # Other installment plans
inputCTest = df[601:1000,3]  # Property
inputDTest = df[601:1000,4]  # Telephone
expectedOutputTest = df[601:1000,5] # Expected training output

count = 0

for index in 1:400
    testInput = convert(Array{Float64,1}, [inputATest[index], inputBTest[index], inputCTest[index], inputDTest[index]])
    testOutput = outputFor(testInput)
    realOutput = convert(Float64, expectedOutputTest[index])

    # The trail can be considered as successful only if the testing value is
    # within the accepted error range
    if testOutput > realOutput - 0.05 && testOutput < realOutput + 0.05
        count += 1
    end

end

### Calculate the Accuracy  ###
accuracy = (count / 400) * 100
println("The accuracy of this trained neural netowrk is $accuracy%")