text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/tensorrt/secure/8.0.3/tars/TensorRT-8.0.3.4.Linux.x86_64-gnu.cuda-11.3.cudnn8.2.tar.gz",
                      "3177435024ff4aa5a6dba8c1ed06ab11cc0e1bf3bb712dfa63a43422f41313f3")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/tensorrt/secure/8.0.3/zip/TensorRT-8.0.3.4.Windows10.x86_64.cuda-11.3.cudnn8.2.zip",
                      "a347d6e7981d0497ba60c5de78716101d73105946e1ff745f0f426f51ea691b0")],
)

max_stars_repo_name: cmhamel/Yggdrasil

max_issues_repo_name: cmhamel/Yggdrasil

max_forks_repo_name: cmhamel/Yggdrasil


text:
################################################################################
# Algebraic identities
#
# FillArrays.jl has `Zero` and `Eye`, which are good, but:
#  - vec + Zero(2) creates a new `vec` instead of returning `vec`, which is wasteful
#  - They include type info and length, which are annoying to specify.
#    We just need a clean algebraic object.

struct Identity end
Base.:*(::Identity, x) = x
Base.:*(x, ::Identity) = x
Base.:*(g::Gaussian, ::Identity) = g  # disambiguation
Base.:*(::Identity, g::Gaussian) = g  # disambiguation
Base.transpose(::Identity) = Identity()
Base.adjoint(::Identity) = Identity()

struct Zero end
Base.:+(x, ::Zero) = x
Base.:+(::Zero, x) = x
Base.:*(x, ::Zero) = zero(x)
Base.:*(::Zero, x) = zero(x)
Base.transpose(z::Zero) = z
Base.adjoint(z::Zero) = z

################################################################################

""" `make_full` is a helper """
make_full(v) = v
make_full(g::Gaussian) = Gaussian(make_full(mean(g)), make_full(cov(g)))
make_full(d::Diagonal) = convert(Matrix, d)
# See StaticArrays#468. We should probably use Diagonal() in 0.7
make_full(d::SDiagonal{1, Float64}) = @SMatrix [d[1,1]]
make_full(d::SDiagonal{2, Float64}) = @SMatrix [d[1,1] 0.0; 0.0 d[2,2]]
make_full(d::SDiagonal{N}) where N =  # this version allocates on 0.6!!!
    convert(SMatrix{N,N}, Diagonal(diag(d)))

################################################################################

no_noise() = Gaussian(Zero(), Zero())
white_noise2(sigma2s...) =
    Gaussian(zero(SVector{length(sigma2s), typeof(zero(sqrt(sigma2s[1])))}),
             SDiagonal(sigma2s))
# fast special-cases
white_noise2(a) =
    # sqrt(zero(a)) is non-differentiable, but zero(sqrt(a)) is.
    # Must take the square root to have consistent units
    Gaussian(SVector(zero(sqrt(a))), SMatrix{1,1}(a))
white_noise2(a, b) = Gaussian(SVector(zero(sqrt(a)), zero(sqrt(b))), SDiagonal(a, b))
white_noise1(args...) = white_noise2((args.^2)...)

################################################################################

Random.rand(RNG::AbstractRNG, g::Gaussian{<:SVector{1}}) =  # type piracy!
    # otherwise calls Cholesky, which fails with unitful Gaussians
    SVector(rand(RNG, Gaussian(mean(g)[1], cov(g)[1])))

parameters(g::Gaussian) = (mean(g), cov(g))   # convenience

# function Base.lufact(m::SMatrix)
#     # Necessary for kalman_smoother until StaticArrays#73 (... I guess?)
#     return lufact(convert(Matrix, m))
#     #return Base.LinAlg.LU(convert(typeof(m), lu.factors), lu.ipiv, lu.info)
# end

# Type piracy! This (mathematically correct) definition improves filtering speed by 3X!
# I believe that it could also easily support a diagonal A, too.
# See definitions in Base.
Base.:\(A::StaticArrays.SArray{Tuple{1,1},<:Any,2,1},
        B::StaticArrays.SArray) = B ./ A[1]

################################################################################
# Marginal variance and standard-deviation

# TODO: use GaussianDistributions.conditional? Would be slower.
marginal_var(g::Gaussian) = diag(cov(g))
marginal_var(g::Gaussian, i::Int) = cov(g) isa Number ? cov(g)[i] : diag(cov(g))[i]
marginal_std(args...) = sqrt(marginal_var(args...))
marginal(g::Gaussian, i::Int) = Gaussian(mean(g)[i], marginal_var(g, i))
is_marginal(g::Gaussian) = dim(g) == 1  # a bit wonky

ext: jl


text:
"""
    download_word_lists()

Downloads and parses Wordle game code for word lists. Returns the list of known
Wordles and the list of valid words.

**Important!** This function will never return future Wordles to prevent
*spoilers (even though they are present in the game code).

# Examples

```julia-repl
wordles, valid_words = download_word_lists()
```
"""
function download_word_lists()
    js_url = let
        html = HTTP.request(:GET, join((WORDLE_URL, "index.html"), '/')).body |> String

        filename = match(r"src=\"(main.\w+.js)\"", html)[1]
        join((WORDLE_URL, filename), '/')
    end

    js = replace(HTTP.request(:GET, js_url).body |> String, "\n" => "")

    # find array of words
    matches = eachmatch(r"\[([\"\w,]+)\]", js) |> collect .|> only .|> String
    matches = replace.(matches, "\"" => "")
    matches = split.(matches, ',')

    # There are several arrays in the game source, but there are two criteria we
    # can use to identify the list answers. First, the number of words in the
    # list must be longer than the number of wordles we've had. Second, the list
    # of words is not in alphabetical order. I happen to know that the first
    # five words are enough to tell if the word list is in alphabetical order.
    target_lists = begin
        lists = filter(matches) do candidate_list
            length(candidate_list) > LATEST_WORDLE_NUMBER
        end

        map(lists) do list
            string.(list)
        end
    end

    is_sorted = map(target_lists) do list
        first_five = list[1:5]
        all(sort(first_five) .== first_five)
    end

    # the original game code started with the second word in the list
    wordles = target_lists[findall(.!(is_sorted)) |> only][2:LATEST_WORDLE_NUMBER + 1]

    # the original valid word list doesn't include wordles so if we do not add
    # ALL of the wordles into the valid word list then future wordles will not
    # be considered valid words
    valid_words = sort(vcat(target_lists...))

    wordles, valid_words
end


text:
doc"""
    Rayleigh(σ)

The *Rayleigh distribution* with scale `σ` has probability density function

$f(x; \sigma) = \frac{x}{\sigma^2} e^{-\frac{x^2}{2 \sigma^2}}, \quad x > 0$

It is related to the [`Normal`](:func:`Normal`) distribution via the property that if $X, Y \sim \operatorname{Normal}(0,\sigma)$, independently, then $\sqrt{X^2 + Y^2} \sim \operatorname{Rayleigh}(\sigma)$.

```julia
Rayleigh()       # Rayleigh distribution with unit scale, i.e. Rayleigh(1.0)
Rayleigh(s)      # Rayleigh distribution with scale s

params(d)        # Get the parameters, i.e. (s,)
scale(d)         # Get the scale parameter, i.e. s
```

External links

* [Rayleigh distribution on Wikipedia](http://en.wikipedia.org/wiki/Rayleigh_distribution)

"""
immutable Rayleigh <: ContinuousUnivariateDistribution
    σ::Float64

    Rayleigh(σ::Real) = (@check_args(Rayleigh, σ > zero(σ)); new(σ))
    Rayleigh() = new(1.0)
end

@distr_support Rayleigh 0.0 Inf

#### Parameters

scale(d::Rayleigh) = d.σ
params(d::Rayleigh) = (d.σ,)

#### Statistics

mean(d::Rayleigh) = sqrthalfπ * d.σ
median(d::Rayleigh) = 1.177410022515474691 * d.σ   # sqrt(log(4.0)) = 1.177410022515474691
mode(d::Rayleigh) = d.σ

var(d::Rayleigh) = 0.429203673205103381 * d.σ^2   # (2.0 - π / 2) = 0.429203673205103381
std(d::Rayleigh) = 0.655136377562033553 * d.σ

skewness(d::Rayleigh) = 0.631110657818937138
kurtosis(d::Rayleigh) = 0.245089300687638063

entropy(d::Rayleigh) = 0.942034242170793776 + log(d.σ)

#### Evaluation

function pdf(d::Rayleigh, x::Float64)
	σ2 = d.σ^2
	x > 0.0 ? (x / σ2) * exp(- (x^2) / (2.0 * σ2)) : 0.0
end

function logpdf(d::Rayleigh, x::Float64)
	σ2 = d.σ^2
	x > 0.0 ? log(x / σ2) - (x^2) / (2.0 * σ2) : -Inf
end

logccdf(d::Rayleigh, x::Float64) = - (x^2) / (2.0 * d.σ^2)
ccdf(d::Rayleigh, x::Float64) = exp(logccdf(d, x))

cdf(d::Rayleigh, x::Float64) = 1.0 - ccdf(d, x)
logcdf(d::Rayleigh, x::Float64) = log1mexp(logccdf(d, x))

quantile(d::Rayleigh, p::Float64) = sqrt(-2.0 * d.σ^2 * log1p(-p))

#### Sampling

rand(d::Rayleigh) = d.σ * sqrt(2.0 * randexp())


text:
module BiobakeryUtils

export
    metaphlan,
    metaphlan_merge,
    metaphlan_profile,
    metaphlan_profiles,
    humann_profile,
    humann_profiles,
    taxfilter,
    taxfilter!,
    parsetaxa,
    parsetaxon,
    rm_strat!,
    permanova,
    humann,
    humann_regroup,
    humann_rename,
    humann_renorm,
    humann_join,
    read_pcl,
    write_pcl,
    humann_barplot,
    humann_barplots,
    kneaddata,
    kneaddata_database,
    kneaddata_read_count_table

using Reexport
@reexport using Microbiome
using CSV
using Tables
using SparseArrays
using Conda
using ReTest

include("utils.jl")
include("metaphlan.jl")
include("humann.jl")
include("kneaddata.jl")

end

max_stars_repo_name: BioJulia/BiobakeryUtils

max_issues_repo_name: BioJulia/BiobakeryUtils

max_forks_repo_name: BioJulia/BiobakeryUtils


text:
using Compat.Random
using Compat.Statistics
using Compat.Test
using PyCall
using Keras
using StatsBase

import Keras: metric
import Keras.Layers: Dense, Activation, LSTM

mse(actual, pred) = Statistics.mean(square(actual - pred))
mae(actual, pred) = Statistics.mean(abs(actual - pred))
rmse(actual, pred) = sqrt(mse(actual, pred))

@testset "Keras" begin
    @testset "Basic Usage" begin
        model = Sequential()
        add!(model, Dense(80, input_dim=30))
        add!(model, Activation(:relu))
        add!(model, Dense(10))
        add!(model, Activation(:softmax))
        compile!(
            model;
            loss=:categorical_crossentropy,
            optimizer=:sgd,
            metrics=[:accuracy]
        )

        h = fit!(model, rand(100, 30), rand(100, 10); epochs=20, batch_size=10, verbose=0)

        @test haskey(h[:history], "acc")
        @test haskey(h[:history], "loss")

        evaluate(model, rand(10, 30), rand(10, 10); batch_size=5, verbose=0)
        predict(model, rand(10, 30); batch_size=5, verbose=0)
    end

    if VERSION >= v"0.7.0"
        @testset "Custom Objectives & Metrics" begin
            model = Sequential()
            add!(model, Dense(20, input_dim=30))
            add!(model, Activation(:relu))
            add!(model, Dense(10))
            add!(model, Activation(:softmax))

            compile!(
                model;
                loss=metric(mse),
                optimizer=:sgd,
                metrics=[:accuracy, metric(mae), metric(rmse)]
            )

            h = fit!(model, rand(100, 30), rand(100, 10); epochs=10, batch_size=10, verbose=0)

            @test haskey(h[:history], "acc")
            @test haskey(h[:history], "loss")
            @test haskey(h[:history], "mae")
            @test haskey(h[:history], "rmse")

            evaluate(model, rand(10, 30), rand(10, 10); batch_size=5, verbose=0)
            predict(model, rand(10, 30); batch_size=5, verbose=0)
        end
    end

    @testset "Tensor Operations" begin
        @testset "Single Tensor Operations" begin
            x = rand(Float32, 4, 3)
            x_t = variable(x)

            @testset "Testing $op" for op in [-, transpose, sqrt, exp, log, round, sin, cos]
                expected, result = if op in [sqrt, exp, log, round, sin, cos]
                    op_result = Keras.eval(op(x_t))
                    bc_result = Keras.eval(broadcast(op, x_t))
                    @test op_result == bc_result

                    broadcast(op, x), bc_result
                else
                    op(x), Keras.eval(op(x_t))
                end

                @test size(expected) == size(result)
                @test op === transpose || typeof(expected) == typeof(result)
                match = all(map(isapprox, expected, result))
                if !match
                    println(expected)
                    println(result)
                end
                @test match
            end

            @testset "Testing $op" for op in [maximum, minimum, sum, prod, Statistics.var, Statistics.std, Statistics.mean]
                expected = Float32(op(x))
                result = Float32(Keras.eval(op(x_t))[1])

                @test isapprox(expected, result; atol=0.1)
            end

            x = rand(Bool, 4)
            x_t = variable(x)
            @testset "Testing $op" for op in [all, any]
                expected = op(x)
                result = Bool(Keras.eval(op(x_t))[1])
                @test expected == result
            end

            expected = broadcast(~, x)
            result = map(i -> i[1], (~Tensor(x)).o)
            @test expected == result
        end

        @testset "Two Tensor Operations" begin
            a = rand(Float32, 4, 3)
            a_t = variable(a)

            b = rand(Float32, 4, 3)
            b_t = variable(b)

            @testset "$op" for op in [==, !=, >, <, >=, <=]
                expected = broadcast(op, a, b)
                result = BitArray(Keras.eval(broadcast(op, a_t, b_t)))
                @test expected == result
            end

            @testset "Testing $op" for op in [mod, ^, -, +, *, /]
                expected = broadcast(op, a, b)
                result = Keras.eval(broadcast(op, a_t, b_t))
                @test size(expected) == size(result)
                @test typeof(expected) == typeof(result)
                @test all(map(isapprox, expected, result))
            end

            b = rand(Float32, 3, 4)
            b_t = variable(b)

            @testset "Testing $op" for op in [*]
                expected = op(a, b)
                result = Keras.eval(op(a_t, b_t))

                @test size(expected) == size(result)
                @test typeof(expected) == typeof(result)
                @test all(map(isapprox, expected, result))
            end

            a = bitrand(4)
            b = bitrand(4)
            expected = broadcast(&, a, b)
            result = map(i -> i[1], (Tensor(a) & Tensor(b)).o)
            @test expected == result

            expected = broadcast(|, a, b)
            result = map(i -> i[1], (Tensor(a) | Tensor(b)).o)
            @test expected == result
        end

        @testset "Custom Tensor Operations" begin
            x = Float32[-1.1, 0.8, -0.5, 1.2]
            x_t = variable(x)

            @test Keras.eval(clip(x_t, -1.0, 1.0)) == Float32[-1.0, 0.8, -0.5, 1.0]
            square(x_t)
            Keras.cast(x_t, :float64)
        end
    end

    @testset "Layers" begin
        @testset "Basic Usage" begin
            a = Keras._layers[:Input](shape=(140, 256))

            lstm = LSTM(32)
            encoded_a = lstm(a)

            W = weights(lstm)
            weights!(lstm, W)
            config(lstm)
            input(lstm)
            output(lstm)
            input_shape(lstm)
            output_shape(lstm)
        end
    end

    @testset "Regularizers" begin
        @testset "Regularizer classes: $rc" for rc in Keras.Regularizers.keras_regularizer_classes
            q = getfield(Keras.Regularizers, Symbol(rc))()
            @test isa(q.obj, PyCall.PyObject)
        end

        @testset "Creating $reg regularizer" for reg in Keras.Regularizers.keras_regularizer_aliases
            q = getfield(Keras.Regularizers, Symbol(reg))()
            @test isa(q, PyCall.PyObject)
        end
    end

    @testset "PyDoc" begin
        output = sprint(show, Keras.PyDoc(Keras._models, :Sequential))
        @test startswith(output, "Linear stack of layers.")
        text = Base.Docs.catdoc(Keras.PyDoc(Keras._models, :Sequential), Keras.PyDoc(Keras._layers, :Dense))
        output = string(text)
        @test startswith(output, "Linear stack of layers.")
        @test endswith(strip(output), "the output would have shape `(batch_size, units)`.")
    end
end


text:
using PlantGeom
using Documenter

DocMeta.setdocmeta!(PlantGeom, :DocTestSetup, :(using PlantGeom; using MultiScaleTreeGraph; using JSServe); recursive = true)

makedocs(;
    modules = [PlantGeom],
    authors = "remi.vezy <VEZY@users.noreply.github.com> and contributors",
    repo = "https://github.com/VEZY/PlantGeom.jl/blob/{commit}{path}#{line}",
    sitename = "PlantGeom.jl",
    format = Documenter.HTML(;
        prettyurls = get(ENV, "CI", "false") == "true",
        canonical = "https://VEZY.github.io/PlantGeom.jl",
        assets = String[]
    ),
    pages = [
        "Home" => "index.md",
        "Plot recipes" => [
            "MTG diagrams" => [
                "Makie.jl" => "plot_diagram/makie_diagram.md",
                "Plots.jl" => "plot_diagram/plots_diagram.md",
            ],
            "3D Makie recipes" => "makie_3d.md"
        ],
        "API" => "API.md"
    ]
)

deploydocs(;
    repo = "github.com/VEZY/PlantGeom.jl",
    devbranch = "main"
)


text:
"""
Returns all rules of a specific type used in a RuleNode
"""
function rulesoftype(node::RuleNode, ruleset::Set{Int})
	retval = Set()

	if node.ind in ruleset
		union!(retval, [node.ind])
	end

	if isempty(node.children)
		return retval
	else
		for child in node.children
			union!(retval, rulesoftype(child, ruleset))
		end

		return retval
	end
end

rulesoftype(node::RuleNode, grammar::GrammarType, ruletype::Symbol) = rulesoftype(node, Set(grammar[ruletype]))

"""
Returns all rules of a specific type used in a RuleNode but not in the ignoreNode
"""
function rulesoftype(node::RuleNode, ruleset::Set{Int}, ignoreNode::RuleNode)
	retval = Set()

	if node == ignoreNode
		return retval
	end

	if node.ind in ruleset
		union!(retval, [node.ind])
	end

	if isempty(node.children)
		return retval
	else
		for child in node.children
			union!(retval, rulesoftype(child, ruleset))
		end

		return retval
	end
end

rulesoftype(node::RuleNode, grammar::GrammarType, ruletype::Symbol, ignoreNode::RuleNode) = rulesoftype(node, Set(grammar[ruletype]), ignoreNode)

"""
Replace a node in expr, specified by path, with new_expr.
Path is a sequence of child indices, starting from the node
"""
function swap_node(expr::RuleNode, new_expr::RuleNode, path::Vector{Int})
	if length(path) == 1
		expr.children[path[begin]] = new_expr
	else
		swap_node(expr.children[path[begin]], new_expr, path[2:end])
	end
end

"""
	Replace child i of a node,  a part of larger expr, with new_expr
"""
function swap_node(expr::RuleNode, node::RuleNode, child_index::Int, new_expr::RuleNode)
	if expr == node
		node.children[child_index] = new_expr
	else
		for child in expr.children
			swap_node(child, node, child_index, new_expr)
		end
	end
end

"""
	Extract derivation sequence from path  (sequence of child indices)

	if the path is deeper than the deepest node, it returns what it has
"""
function get_rulesequence(node::RuleNode, path::Vector{Int})
	if node.ind == 0 # sign for empty node
		return Vector{Int}()
	elseif isempty(node.children) # no childnen, nowehere to follow the path; still return the index
		return [node.ind]
	elseif isempty(path)
		return [node.ind]
	elseif isassigned(path, 2)
		# at least two items are left in the path
		# need to access the child with get because it can happen that the child is not yet built
		return append!([node.ind], get_rulesequence(get(node.children, path[begin], RuleNode(0)), path[2:end]))
	else
		# if only one item left in the path
		# need to access the child with get because it can happen that the child is not yet built
		return append!([node.ind], get_rulesequence(get(node.children, path[begin], RuleNode(0)), Vector{Int}()))
	end
end

"""
	Extracts rules in the left subtree defined by the path
"""
function rulesonleft(expr::RuleNode, path::Vector{Int})
	if isempty(expr.children)
		# if the encoutered node is terminal or non-expanded non-terminal, return node id
		Set{Int}(expr.ind)
	elseif isempty(path)
		# if path is empty, collect the entire subtree
		ruleset = Set{Int}(expr.ind)
		for ch in expr.children
			union!(ruleset, rulesonleft(ch, Vector{Int}()))
		end
		return ruleset
	elseif length(path) == 1
		# if there is only one element left in the path, collect all children except the one indicated in the path
		ruleset = Set{Int}(expr.ind)
		for i in 1:path[begin]-1
			union!(ruleset, rulesonleft(expr.children[i], Vector{Int}()))
		end
		return ruleset
	else
		# collect all subtrees up to the child indexed in the path
		ruleset = Set{Int}(expr.ind)
		for i in 1:path[begin]-1
			union!(ruleset, rulesonleft(expr.children[i], Vector{Int}()))
		end
		union!(ruleset, rulesonleft(expr.children[path[begin]], path[2:end]))
		return ruleset
	end
end


text:
module MultiCellMapsTests

using Test
using Gridap
import Gridap: ∇

u1fun(x) = x[1] + x[2]
u2fun(x) = x[1] - x[2]

u1fun_grad(x) = VectorValue(1.0,1.0)
u2fun_grad(x) = VectorValue(1.0,-1.0)

∇(::typeof(u1fun)) = u1fun_grad
∇(::typeof(u2fun)) = u2fun_grad

b1fun(x) = u2fun(x)
b2fun(x) = 0.0

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

order = 1
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

V1 = TestFESpace(fespace)
V2 = V1

V = [V1,V2]
U = V

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

b1field = CellField(trian,b1fun)
b2field = CellField(trian,b2fun)

function a(v,u)
  a11 = varinner(∇(v[1]), ∇(u[1]))
  a12 = varinner(v[1],u[2])
  a22 = varinner(∇(v[2]), ∇(u[2]))
  return CellMap[a11,a12,a22], [(1,1),(1,2),(2,2)]
end

function b(v)
  b1 = varinner(v[1],b1field)
  b2 = varinner(v[2],b2field)
  return CellMap[b1,b2], [(1,),(2,)]
end

v = [ CellBasis(Vi) for Vi in V ]
u = [ CellBasis(Ui) for Ui in U ]

blocks, fieldids = a(v,u)

mcm = MultiCellMap(blocks,fieldids)

mca = integrate(mcm,trian,quad)

@test isa(mca,MultiCellArray{Float64,2})

u1 = CellField(trian,u1fun)
u2 = CellField(trian,u2fun)
u = [u1,u2]

blocks, fieldids = a(v,u)

fieldids = [ (fi[1],) for fi in fieldids ]

mcm = MultiCellMap(blocks,fieldids)
mcm1 = mcm

mca = integrate(mcm,trian,quad)

@test isa(mca,MultiCellArray{Float64,1})

blocks, fieldids = b(v)

mcm = MultiCellMap(blocks,fieldids)
mcm2 = mcm

mca = integrate(mcm,trian,quad)

@test isa(mca,MultiCellArray{Float64,1})

mcm = mcm1 + mcm2
mca = integrate(mcm,trian,quad)
@test isa(mca,MultiCellArray{Float64,1})

mcm = mcm1 - mcm2
mca = integrate(mcm,trian,quad)
@test isa(mca,MultiCellArray{Float64,1})

end


text:
function generictest_dixon_Fp(G, p = Characters.dixon_prime(G))
    F = Characters.FiniteFields.GF{p}
    ccG = Characters.conjugacy_classes(G)
    Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)]
    @test isdiag(Characters.common_esd(Ns, F))
    esd = Characters.common_esd(Ns, F)

    let W = esd.basis
        @test inv(W) isa Matrix{F}
        @test isone(inv(W) * W)
        iW = inv(W)
        for N in Ns
            # basis actually diagonalizes all Ns
            m = W * N * iW
            @test isdiag(m)
        end
    end

    tbl = Characters.CharacterTable(F, G, ccG)
    @test Characters.irreducible_characters(tbl) isa
          Vector{<:Characters.Character}

    chars = Characters.irreducible_characters(tbl)

    # checking the degrees
    @test sum(Int.(SymbolicWedderburn.degree.(chars)) .^ 2) == order(G)

    # orthogonality of characters over F
    @test [dot(χ, ψ) for χ in chars, ψ in chars] == I
end

function generictest_dixon_C(G, p = Characters.dixon_prime(G))
    F = Characters.FiniteFields.GF{p}
    ccG = Characters.conjugacy_classes(G)
    tbl = Characters.CharacterTable(F, G, ccG)
    chars_Fp = Characters.irreducible_characters(tbl)

    degrees = SymbolicWedderburn.degree.(chars_Fp)

    m = Characters._multiplicities(chars_Fp)
    for i = 1:size(m, 1)
        @test all(Int.(m[i, :, :]) .<= degrees[i])
    end

    @test Characters.complex_character_table(Rational{Int}, tbl) isa
          Characters.CharacterTable{<:Group, <:Cyclotomic{Rational{Int}}}

    tblCC = Characters.complex_character_table(Rational{Int}, tbl)
    chars_CC = Characters.irreducible_characters(tblCC)

    id = one(G)
    @test [ψ(id) for ψ in chars_CC] == degrees

    # orthogonality of characters over ℂ
    @test [dot(χ, ψ) for χ in chars_CC, ψ in chars_CC] == I

    @test sum(χ -> SymbolicWedderburn.degree(χ)^2, chars_CC) == order(G)
end

@testset "Dixon Algorithm" begin
    @testset "DixonPrimes" begin
        @test Characters.dixon_prime(20, 20) == 41

        @testset "random" begin
            for n in rand(2:1000, 10)
                F = Characters.Primes.factor(n)
                e = lcm(collect(keys(F)))
                p = Characters.dixon_prime(n, e)
                @test mod(p, e) == 1
                @test p > 2 * floor(sqrt(n))
            end
        end

        @testset "DixonPrimeGroups" begin
            G = PermutationGroups.SymmetricGroup(4)
            ccG = Characters.conjugacy_classes(G)
            @test exponent(G) == 12
            @test exponent(ccG) == 12
            @test Characters.dixon_prime(G) ==
                  Characters.dixon_prime(ccG)
        end
    end

    @testset "Dixon over GF{p}" begin
        @testset "Example: Alt(4)" begin
            G = PermGroup([perm"(1,2,3)(4)", perm"(2,3,4)"])
            ccG = let S = gens(G)
                [
                    Orbit(S, one(G)),
                    Orbit(S, G(perm"(2,3,4)")),
                    Orbit(S, G(perm"(2,4,3)")),
                    Orbit(S, G(perm"(1,2)(3,4)")),
                ] # the order is taken from GAP
            end

            let ccG = ccG, p = 29
                F = Characters.FiniteFields.GF{p}
                Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)]
                @test_throws AssertionError Characters.common_esd(Ns, F)
            end

            let ccG = ccG, p = 31
                F = Characters.FiniteFields.GF{p}
                Ns = [Characters.CMMatrix(ccG, i) for i = 1:length(ccG)]
                esd = Characters.EigenSpaceDecomposition(F.(Ns[1]))
                esd = Characters.refine(esd, F.(Ns[2]))
                esd = Characters.refine(esd, F.(Ns[3]))
                @test isdiag(esd)

                @test isdiag(Characters.common_esd(Ns, F))
            end

            # generic tests
            generictest_dixon_Fp(G)

            # tests specific to G
            p = Characters.dixon_prime(ccG)
            F = Characters.FiniteFields.GF{p}

            tbl = Characters.CharacterTable(F, G)
            chars = Characters.irreducible_characters(tbl)

            @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 1, 3]

            @test Set(Int.(values(χ)) for χ in chars) ==
                  Set([[3, 0, 0, 6], [1, 4, 2, 1], [1, 2, 4, 1], [1, 1, 1, 1]])

        end
    end

    @testset "Lifting characters to ℂ" begin
        @testset "Example: Alt(4)" begin
            G = PermGroup([perm"(1,2,3)(4)", perm"(2,3,4)"])
            ccG = let S = gens(G)
                [
                    Orbit(S, one(G)),
                    Orbit(S, G(perm"(2,3,4)")),
                    Orbit(S, G(perm"(2,4,3)")),
                    Orbit(S, G(perm"(1,2)(3,4)")),
                ] # the order of cclasses is taken from GAP
            end

            @testset "PowerMap" begin
                pmG = Characters.PowerMap(ccG)
                @test size(pmG) == (4, 6)
                @test pmG[:, 0] == ones(Int, 4)
                @test pmG[:, 1] == 1:4
                @test pmG[:, 2] == [1, 3, 2, 1]
                @test pmG[:, 3] == [1, 1, 1, 4]
                @test pmG[:, 4] == [1, 2, 3, 1]
                @test pmG[:, 5] == [1, 3, 2, 4]
                @test all(isone.(first.(ccG) .^ 6))
            end

            # generic tests
            generictest_dixon_C(G)

            chars_C = Characters.irreducible_characters(G, ccG)
            E = Characters.Cyclotomics.E

            @test [collect(values(χ)) for χ in chars_C] == [
                E(3, 0) .* [3, 0, 0, -1],
                E(3, 0) .* [1, E(3, 2), E(3, 1), 1],
                E(3, 0) .* [1, E(3, 1), E(3, 2), 1],
                E(3, 0) .* [1, 1, 1, 1],
            ]
        end

        @testset "Example: Sym(4)" begin
            G = PermGroup([perm"(1,2)", perm"(1,2,3,4)"])
            ccG = let S = gens(G)
                ccG = [
                    Orbit(S, one(G)),
                    Orbit(S, G(perm"(1,2)(4)")),
                    Orbit(S, G(perm"(1,2)(3,4)")),
                    Orbit(S, G(perm"(1,2,3)(4)")),
                    Orbit(S, G(perm"(1,2,3,4)")),
                ] # the order of cclasses is taken from GAP
            end
            generictest_dixon_Fp(G)
            generictest_dixon_C(G)

            chars = Characters.irreducible_characters(G, ccG)

            @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 2, 3, 3]
            @test [collect(values(χ)) for χ in chars] == [
                [2, 0, 2, -1, 0],
                [3, 1, -1, 0, -1],
                [1, 1, 1, 1, 1],
                [1, -1, 1, 1, -1],
                [3, -1, -1, 0, 1],
            ]
        end

        @testset "example: C₅ ⋊ C₄" begin
            G = PermGroup(perm"(1,2,4,5,3)", perm"(2,5,3,4)")
            S = gens(G)

            ccG = [
                Orbit(S, one(G)),
                Orbit(S, G(perm"(2,3)(4,5)")),
                Orbit(S, G(perm"(2,4,3,5)")),
                Orbit(S, G(perm"(2,5,3,4)")),
                Orbit(S, G(perm"(1,2,4,5,3)")),
            ]
            # the order of cclasses is taken from GAP
            generictest_dixon_Fp(G)
            generictest_dixon_C(G)

            chars = Characters.irreducible_characters(G, ccG)

            @test sort(SymbolicWedderburn.degree.(chars)) == [1, 1, 1, 1, 4]
            @test [collect(values(χ)) for χ in chars] == [
                E(4, 0) .* [4, 0, 0, 0, -1],
                E(4, 0) .* [1, 1, 1, 1, 1],
                E(4, 0) .* [1, 1, -1, -1, 1],
                E(4, 0) .* [1, -1, E(4), -E(4), 1],
                E(4, 0) .* [1, -1, -E(4), E(4), 1],
            ]
        end
    end

    @testset "Different base rings for characters_dixon" begin
        G = PermGroup(perm"(1,2,3,4)")
        @test eltype(Characters.irreducible_characters(Rational{Int}, G)) <:
            Characters.Character{<:Cyclotomic{Rational{Int}}}
        @test eltype(Characters.irreducible_characters(Rational{BigInt}, G)) <:
            Characters.Character{<:Cyclotomic{Rational{BigInt}}}
    end

    @time @testset "SmallPermGroups" begin
        for (ord, groups) in SmallPermGroups
            @testset "SmallGroup($ord, $n)" for (n, G) in enumerate(groups)
                @test Characters.irreducible_characters(G) isa
                      Vector{<:Characters.Character}
                generictest_dixon_Fp(G)
                generictest_dixon_C(G)
            end
        end
    end
end


text:
"""
    Quat{T} <: Rotation{T}
    Quat(w, x, y, z)

The `Quat` type is a 3×3 matrix representation of a normalized quaternion.
They allow you to transparently use (fast) quaternion algebra to store, compose
and invert 3D rotations, while at the same time letting you apply rotations
through matrix-vector multiplication.

Note: the constructor will always renormalize the input so that the quaternion
has length 1 (w² + x² + y² + z² = 1), and the rotation matrix is orthogonal.
"""
immutable Quat{T} <: Rotation{T}
    w::T
    x::T
    y::T
    z::T

    # For the moment we ensure that it is normalized upon construction.
    function Quat(w, x, y, z)
        norm = copysign(sqrt(w*w + x*x + y*y + z*z), w)
        # Should this be an error or warning, if it isn't approximately normalized? E.g.:
        #if norm !≈ 1
        #    error("Expected a normalized quaternion") # or warn() ?
        #end
        new(w/norm, x/norm, y/norm, z/norm)
    end
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 4 inputs (not 9).
@inline (::Type{Quat}){W,X,Y,Z}(w::W, x::X, y::Y, z::Z) = Quat{promote_type(promote_type(promote_type(W, X), Y), Z)}(w, x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
function (::Type{Q}){Q<:Quat}(t::NTuple{9})
    q = Q(sqrt(abs(1  + t[1] + t[5] + t[9])) / 2,
               copysign(sqrt(abs(1 + t[1] - t[5] - t[9]))/2, t[6] - t[8]),
               copysign(sqrt(abs(1 - t[1] + t[5] - t[9]))/2, t[7] - t[3]),
               copysign(sqrt(abs(1 - t[1] - t[5] + t[9]))/2, t[2] - t[4]))
end

function Base.getindex(q::Quat, i::Integer)
    if i == 1
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww + xx - yy - zz
    elseif i == 2
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy + zw)
    elseif i == 3
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz - yw)
    elseif i == 4
        xy = (q.x * q.y)
        zw = (q.w * q.z)

        2 * (xy - zw)
    elseif i == 5
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx + yy - zz
    elseif i == 6
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz + xw)
    elseif i == 7
        xz = (q.x * q.z)
        yw = (q.y * q.w)

        2 * (xz + yw)
    elseif i == 8
        yz = (q.y * q.z)
        xw = (q.w * q.x)

        2 * (yz - xw)
    elseif i == 9
        ww = (q.w * q.w)
        xx = (q.x * q.x)
        yy = (q.y * q.y)
        zz = (q.z * q.z)

        ww - xx - yy + zz
    else
        throw(BoundsError(r,i))
    end
end

# This speeds up some StaticArray methods (such as conversion to RotMatrix)
function Base.convert(::Type{Tuple}, q::Quat)
    ww = (q.w * q.w)
    xx = (q.x * q.x)
    yy = (q.y * q.y)
    zz = (q.z * q.z)
    xy = (q.x * q.y)
    zw = (q.w * q.z)
    xz = (q.x * q.z)
    yw = (q.y * q.w)
    yz = (q.y * q.z)
    xw = (q.w * q.x)

    # initialize rotation part
    return (ww + xx - yy - zz,
            2 * (xy + zw),
            2 * (xz - yw),
            2 * (xy - zw),
            ww - xx + yy - zz,
            2 * (yz + xw),
            2 * (xz + yw),
            2 * (yz - xw),
            ww - xx - yy + zz)
end

function Base.:*(q::Quat, v::StaticVector)
    if length(v) != 3
        throw("Dimension mismatch: cannot rotate a vector of length $(length(v))")
    end

    qo = (-q.x * v[1] - q.y * v[2] - q.z * v[3],
           q.w * v[1] + q.y * v[3] - q.z * v[2],
           q.w * v[2] - q.x * v[3] + q.z * v[1],
           q.w * v[3] + q.x * v[2] - q.y * v[1])

    T = promote_type(eltype(q), eltype(v))

    return similar_type(v, T)(-qo[1] * q.x + qo[2] * q.w - qo[3] * q.z + qo[4] * q.y,
                              -qo[1] * q.y + qo[2] * q.z + qo[3] * q.w - qo[4] * q.x,
                              -qo[1] * q.z - qo[2] * q.y + qo[3] * q.x + qo[4] * q.w)
end

# TODO consider Ac_mul_B, etc...
function Base.:*(q1::Quat, q2::Quat)
    Quat(q1.w*q2.w - q1.x*q2.x - q1.y*q2.y - q1.z*q2.z,
            q1.w*q2.x + q1.x*q2.w + q1.y*q2.z - q1.z*q2.y,
            q1.w*q2.y - q1.x*q2.z + q1.y*q2.w + q1.z*q2.x,
            q1.w*q2.z + q1.x*q2.y - q1.y*q2.x + q1.z*q2.w)
end

function Base.inv(q::Quat)
    Quat(q.w, -q.x, -q.y, -q.z)
end

@inline Base.eye(::Type{Quat}) = Quat(1.0, 0.0, 0.0, 0.0)
@inline Base.eye{T}(::Type{Quat{T}}) = Quat{T}(one(T), zero(T), zero(T), zero(T))

################################################################################
################################################################################
"""
    immutable SPQuat{T} <: Rotation{T}
    SPQuat(x, y, z)

An `SPQuat` is a 3D rotation matrix represented by the "stereographic projection" of a normalized quaternion (shortened to "SPQuat"), which is
a 3-element parametrization of a unit quaternion Q formed by the intersection of a line from [-1,0,0,0] to Q, with a plane containing the origin and with normal direction [1,0,0,0]. This
is a compact representation of rotations where the derivitives of the rotation matrix's elements w.r.t. the SPQuat parameters are rational functions (making them useful for optimization).

See:

    [1] "A Recipe on the Parameterization of Rotation Matrices for Non-Linear Optimization using Quaternions",
    Terzakis, Culverhouse, Bugmann, Sharma, Sutton,
    MIDAS technical report, 2011
    http://www.tech.plymouth.ac.uk/sme/springerusv/2011/publications_files/Terzakis%20et%20al%202012,%20A%20Recipe%20on%20the%20Parameterization%20of%20Rotation%20Matrices...MIDAS.SME.2012.TR.004.pdf

    Note 1: the singularity (origin) has been moved from [0,0,0,-1] in Ref[1] to [-1,0,0,0], so the 0 rotation quaternion [1,0,0,0] maps to [0,0,0] as opposed of to [1,0,0].
    Note 2: the quaternion rotation ambiguity q = -q means there will be a rotation with ||SPQuat|| <= 1 and an equivilent rotation with ||SPQuat|| > 1.  This package
            uses the solution with ||SPQuat|| <= 1
    Note 3: it is safe to assume that the corresponding matrix is orthogonal/unitary for any input x, y, z.

"""
immutable SPQuat{T} <: Rotation{T}
    x::T
    y::T
    z::T

    # TODO should we enforce norm <= 1?
end

# StaticArrays will take over *all* the constructors and put everything in a tuple...
# but this isn't quite what we mean when we have 3 inputs (not 9).
@inline (::Type{SPQuat}){X,Y,Z}(x::X, y::Y, z::Z) = SPQuat{promote_type(promote_type(X, Y), Z)}(x, y, z)

# These 2 functions are enough to satisfy the entire StaticArrays interface:
@inline (::Type{SPQ}){SPQ <: SPQuat}(t::NTuple{9}) = SPQ(Quat(t))
@inline Base.getindex(spq::SPQuat, i::Integer) = Quat(spq)[i]

@inline function Base.convert{Q <: Quat}(::Type{Q}, spq::SPQuat)
    # Both the sign and norm of the Quat is automatically dealt with in its inner constructor
    return Q(1 - (spq.x*spq.x + spq.y*spq.y + spq.z*spq.z), 2*spq.x, 2*spq.y, 2*spq.z)
end

@inline function Base.convert{SPQ <: SPQuat}(::Type{SPQ}, q::Quat)
    alpha2 = (1 - q.w) / (1 + q.w) # <= 1 since q.w >= 0
    spq = SPQ(q.x * (alpha2 + 1)*0.5,  q.y * (alpha2 + 1)*0.5, q.z * (alpha2 + 1)*0.5)
end

@inline Base.convert(::Type{Tuple}, spq::SPQuat) = Tuple(Quat(spq))

@inline Base.:*(spq::SPQuat, x::StaticVector) = Quat(spq) * x

@inline Base.:*(spq::SPQuat, r::Rotation) = Quat(spq) * r
@inline Base.:*(spq::SPQuat, r::RotMatrix) = Quat(spq) * r
@inline Base.:*(r::Rotation, spq::SPQuat) = r * Quat(spq)
@inline Base.:*(r::RotMatrix, spq::SPQuat) = r * Quat(spq)
@inline Base.:*(spq1::SPQuat, spq2::SPQuat) = Quat(spq1) * Quat(spq2)

@inline Base.inv(spq::SPQuat) = SPQuat(-spq.x, -spq.y, -spq.z)

@inline Base.eye(::Type{SPQuat}) = SPQuat(0.0, 0.0, 0.0)
@inline Base.eye{T}(::Type{SPQuat{T}}) = SPQuat{T}(zero(T), zero(T), zero(T))

# rotation properties
@inline rotation_angle(spq::SPQuat) = rotation_angle(Quat(spq))
@inline rotation_axis(spq::SPQuat) = rotation_axis(Quat(spq))


text:
module BenchCircularBuffer

using DataStructures
using BenchmarkTools
using Random

const SUITE = BenchmarkGroup()
const CAPS = [1024, 1024^2]

# Empirically, 1 second is sufficient for consistent timings here.
BenchmarkTools.DEFAULT_PARAMETERS.seconds = 1

function init_cb(a)
    len = length(a)
    cb = CircularBuffer{Int}(len)
    cb.first = 1 + div(len, 2)
    return append!(cb, a)
end

function perf_chained_getindex(cb)
    i = 1
    for k = 1:capacity(cb)
        i = cb[i]
    end
    return i
end

function perf_first(cb)
    total = 0
    for i = 1:capacity(cb)
        cb.first = i
        total += first(cb)
    end
    return total
end

function perf_last(cb)
    total = 0
    for i = 1:capacity(cb)
        cb.first = i
        total += last(cb)
    end
    return total
end

g = addgroup!(SUITE, "access")

# Creating a new copy for each sample seems to eliminate an unknown but substantial source
# of timing variation on some machines. This is why each benchmark has `setup=(cb = deepcopy($cb))`.
for cap in CAPS
    Random.seed!(0)
    cb = init_cb(randcycle(cap))

    g["chained_getindex", cap] = @benchmarkable perf_chained_getindex(cb) setup=(cb = deepcopy($cb))

    g["sum", cap] = @benchmarkable sum(cb) setup=(cb = deepcopy($cb))
    g["foldl", cap] = @benchmarkable foldl(+, cb) setup=(cb = deepcopy($cb))
    g["foldr", cap] = @benchmarkable foldr(+, cb) setup=(cb = deepcopy($cb))

    g["first", cap] = @benchmarkable perf_first(cb) setup=(cb = deepcopy($cb))
    g["last", cap] = @benchmarkable perf_last(cb) setup=(cb = deepcopy($cb); cb.length = div($cap, 2))

    g["convert", cap] = @benchmarkable convert(Array, cb) setup=(cb = deepcopy($cb))
end

function perf_setindex!(cb, indices)
    for i in indices
        cb[i] = i
    end
    return cb
end

function perf_push!(cb)
    for i = 1:capacity(cb)
        push!(cb, i)
    end
    return cb
end

function perf_pushfirst!(cb)
    for i = 1:capacity(cb)
        pushfirst!(cb, i)
    end
    return cb
end

g = addgroup!(SUITE, "mutate")

for cap in CAPS
    cb = init_cb(-1:-1:-cap)

    Random.seed!(0)
    g["random_setindex!", cap] = @benchmarkable perf_setindex!(cb, $(randperm(cap))) setup=(cb = deepcopy($cb))
    g["forward_setindex!", cap] = @benchmarkable perf_setindex!(cb, 1:$cap) setup=(cb = deepcopy($cb))
    g["reverse_setindex!", cap] = @benchmarkable perf_setindex!(cb, $cap:-1:1) setup=(cb = deepcopy($cb))

    g["full_push!", cap] = @benchmarkable perf_push!(cb) setup=(cb = deepcopy($cb))
    g["full_pushfirst!", cap] = @benchmarkable perf_pushfirst!(cb) setup=(cb = deepcopy($cb))
end

perf_empty_push!(cb) = (empty!(cb); perf_push!(cb))
perf_empty_pushfirst!(cb) = (empty!(cb); perf_pushfirst!(cb))

function perf_pop!(cb)
    cap = capacity(cb)
    cb.length = cap
    total = 0
    for i = 1:cap
        total += pop!(cb)
    end
    return total
end

function perf_popfirst!(cb)
    cap = capacity(cb)
    cb.length = cap
    total = 0
    for i = 1:cap
        total += popfirst!(cb)
    end
    return total
end

perf_fill!(cb) = (empty!(cb); fill!(cb, 1))

for cap in CAPS
    cb = init_cb(-1:-1:-cap)
    g["empty_push!", cap] = @benchmarkable perf_empty_push!(cb) setup=(cb = deepcopy($cb))
    g["empty_pushfirst!", cap] = @benchmarkable perf_empty_pushfirst!(cb) setup=(cb = deepcopy($cb))
    g["pop!", cap] = @benchmarkable perf_pop!(cb) setup=(cb = deepcopy($cb))
    g["popfirst!", cap] = @benchmarkable perf_popfirst!(cb) setup=(cb = deepcopy($cb))
    g["fill!", cap] = @benchmarkable perf_fill!(cb) setup=(cb = deepcopy($cb))
end

end  # module BenchCircularBuffer

if @isdefined(SUITE)
    BenchCircularBuffer.SUITE
else
    # This `else` branch allows this file to be called directly by PkgBenchmark via the
    # `script` keyword in `PkgBenchmark.benchmarkpkg`.
    using BenchmarkTools
    const SUITE = BenchmarkGroup()
    SUITE["CircularBuffer"] = BenchCircularBuffer.SUITE
end


text:
# Autogenerated wrapper script for Severo_jll for i686-linux-musl-cxx03
export libcell

using OpenBLAS_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("Severo")
JLLWrappers.@declare_library_product(libcell, "libcell.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libcell,
        "lib/libcell.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using RNGPool
using Plots
Plots.gr()

include("demo.jl")
include("lgModel.jl")

verbose = false

setRNGs(12345)
model, lM, _ = setupLGModel(1000, 1.0)
@time b11 = CCSMCBoundaries(model, lM, 2, :BS, 10000); plot(b11)
@time b12 = CCSMCBoundaries(model, lM, 4, :BS, 10000); plot!(b12)
@time b13 = CCSMCBoundaries(model, lM, 8, :BS, 10000); plot!(b13)
@time b14 = CCSMCBoundaries(model, lM, 16, :BS, 10000); plot!(b14)
@time b15 = CCSMCBoundaries(model, lM, 32, :BS, 10000); plot!(b15)
plot(b11, label="N=2")
plot!(b12, label="N=4")
plot!(b13, label="N=8")
plot!(b14, label="N=16")
plot!(b15, label="N=32")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for obsVariance = 1.0")
savefig("output/lgDemoBoundaries1.pdf")

setRNGs(12345)
model, lM, _ = setupLGModel(1000, 100.0)
@time b21 = CCSMCBoundaries(model, lM, 2, :BS, 10000); plot(b21)
@time b22 = CCSMCBoundaries(model, lM, 4, :BS, 10000); plot!(b22)
@time b23 = CCSMCBoundaries(model, lM, 8, :BS, 10000); plot!(b23)
@time b24 = CCSMCBoundaries(model, lM, 16, :BS, 10000); plot!(b24)
@time b25 = CCSMCBoundaries(model, lM, 32, :BS, 10000); plot!(b25)
plot(b21, label="N=2")
plot!(b22, label="N=4")
plot!(b23, label="N=8")
plot!(b24, label="N=16")
plot!(b25, label="N=32")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for obsVariance = 100.0")
savefig("output/lgDemoBoundaries2.pdf")

setRNGs(12345)
model, lM, _ = setupLGModel(1000, 1000.0)
@time b35 = CCSMCBoundaries(model, lM, 32, :BS, 500); plot(b35)
@time b36 = CCSMCBoundaries(model, lM, 64, :BS, 500); plot!(b36)
@time b37 = CCSMCBoundaries(model, lM, 128, :BS, 500); plot!(b37)
@time b38 = CCSMCBoundaries(model, lM, 256, :BS, 500); plot!(b38)
plot(b35, label="N=32")
plot!(b36, label="N=64")
plot!(b37, label="N=128")
plot!(b38, label="N=256")
xlabel!("iteration")
ylabel!("boundary")
title!("Boundary against iteration for obsVariance = 1000.0")
savefig("output/lgDemoBoundaries3.pdf")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct TagFilter <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    op::Any # spec type: Union{ Nothing, String } # spec name: op
    value::Any # spec type: Union{ Nothing, String } # spec name: value

    function TagFilter(;name=nothing, op=nothing, value=nothing)
        o = new()
        validate_property(TagFilter, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(TagFilter, Symbol("op"), op)
        setfield!(o, Symbol("op"), op)
        validate_property(TagFilter, Symbol("value"), value)
        setfield!(o, Symbol("value"), value)
        o
    end
end # type TagFilter

const _property_map_TagFilter = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("op")=>Symbol("op"), Symbol("value")=>Symbol("value"))
const _property_types_TagFilter = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("op")=>"String", Symbol("value")=>"String")
Base.propertynames(::Type{ TagFilter }) = collect(keys(_property_map_TagFilter))
Swagger.property_type(::Type{ TagFilter }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TagFilter[name]))}
Swagger.field_name(::Type{ TagFilter }, property_name::Symbol) =  _property_map_TagFilter[property_name]

function check_required(o::TagFilter)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("op")) === nothing) && (return false)
    (getproperty(o, Symbol("value")) === nothing) && (return false)
    true
end

function validate_property(::Type{ TagFilter }, name::Symbol, val)
    if name === Symbol("name")
        Swagger.validate_param(name, "TagFilter", :maxLength, val, 128)
        Swagger.validate_param(name, "TagFilter", :minLength, val, 1)
    end
    if name === Symbol("value")
        Swagger.validate_param(name, "TagFilter", :maxLength, val, 256)
        Swagger.validate_param(name, "TagFilter", :minLength, val, 0)
    end
end


text:
using Test: occursin
using ComputationManager
using Test, IOCapture

@testset "ComputationManager.jl" begin
    ComputationManager.main([])

    # Store a job. The function must return a callable object
    savejob("job1") do name
        () -> @info "Running job: $(name)"
    end

    # By default, we print out all the job names
    c = IOCapture.capture() do
        ComputationManager.main([])
    end
    @test occursin("job1", c.output)

    # Let's add another job
    savejob("job2") do name
        () -> @info "Running job: $(name)"
    end
    c = IOCapture.capture() do
        ComputationManager.main([])
    end
    @test occursin("job1", c.output)
    @test occursin("job2", c.output)

    # Let's run jobs 1 and 2 separately
    c = IOCapture.capture() do; ComputationManager.main(["job1"]); end
    @test occursin("Running job: job1", c.output)
    @test !occursin("Running job: job2", c.output)
    c = IOCapture.capture() do; ComputationManager.main(["job2"]); end
    @test !occursin("Running job: job1", c.output)
    @test occursin("Running job: job2", c.output)
    c = IOCapture.capture() do; ComputationManager.main(["job1", "job2"]); end
    @test occursin("Running job: job1", c.output)
    @test occursin("Running job: job2", c.output)
    c = IOCapture.capture() do; ComputationManager.main(["--runall"]); end
    @test occursin("Running job: job1", c.output)
    @test occursin("Running job: job2", c.output)

    # Test the savejob(::NamedTuple) signature
    empty!(ComputationManager._JOBS) # reset global state
    savejob( (name = "job3", f = () -> @info "Running job: job3") )
    c = IOCapture.capture() do; ComputationManager.main(["--runall"]); end
    @test occursin("Running job: job3", c.output)
end


text:
__precompile__()

module RowTables

import DataFrames
import DataFrames.rename!
import DataStructures.OrderedDict
import Random
import JSON
using Printf

export AbstractRowTable, RowTable, rows, columns, rowdict

include("index.jl")
include("abstractrowtable.jl")
include("rowtable.jl")
include("display.jl")
include("sort.jl")

end  # module RowTables


text:
#####
#
# Example script for TintiNet.jl
# Model IGBT (InceptGOR-Bert): 3-state SS Classifier
#
# To run training, revise your code and uncomment the last line.
# On slower GPUs, an epoch may take several minutes with the default configuration.
# Adjust model size using the `hyperpar` and `config` dictionaries.
#
#####

#####
# 0. Loading Packages and defining Hyperparameters
#####

using CUDA, Flux, Transformers, JSON, BSON, Random, Dates, Printf
using IterTools: ncycle;
using Transformers.Basic, Transformers.Pretrain, Transformers.Datasets, Transformers.Stacks, Transformers.BidirectionalEncoder
using BSON: @load, @save
using Flux: @functor, onehot, onecold, mse, shuffle, sigmoid, gradient, unsqueeze
using Flux.Optimise: update!
using Statistics: mean, std

using TintiNet
dformat = Dates.DateFormat("yyyy-m-d-HH-MM-SS")
Random.seed!(102528)

##
# Algorithm Hyperparameters

hyperpar = Dict(
    "max_seq_len"       => 128,     # Maximum sequence length on the model (default 192)
    "batch_size"        => 128,     # Batch size for training (default 128)
    "number_of_epochs"  => 1000,      # Number of epochs for training (10 for the example file)
    "learning_rate"     => 1e-3,    # Learning rate for optimizer (default 0.001)
    "smooth_epsilon"    => 1e-6     # Small value to smooth onehot result before logcrossentropy (default 1e-6)
)

##
# Model Configuration

const config = Dict(
  "input_vocab_size"                => 22,      # Size of the SEQUENCE/FASTA vocabulary, usually 20 aminoacids + "X" + "-"
  "output_vocab_size"               => 4,       # Size of the STRUCTURE/SS3 vocabulary, usually 3 classes + "-"
  "embed_dimension"                 => 128,      # Size of the individual Embed vectors
  "bert_num_hidden_layers"          => 2,       # Number of consecutive Transformer layers on the BERT encoder
  "bert_transformer_layer_heads"    => 8,       # Number of heads in each Transformer layer
  "bert_transformer_layer_hidden"   => 128,     # Hidden size of the positionwise linear layer in each Transformer layer
  "bert_transformer_dropout_prob"   => 0.1,     # Dropout probability for the Transformer layers
  "bert_transformer_activation"     => gelu,    # Activation function for each Transformer layer (default: gelu)
  "detection_num_filters"           => 32       # Number of filters in each convolutional decoder layer
)

#####
# 2. Obtaining and shaping the data
#####

## 2.1. Training sequences, structures and lengths

training_set_datapath = "../data/fold01training_dataset.json"

open(training_set_datapath, "r") do file

    entries = JSON.parse(file)
    num_entries = length(entries)

    parsed_sequences = Vector{Vector{String}}(undef, num_entries)
    parsed_structures = Vector{Vector{String}}(undef, num_entries)
    parsed_lengths = Vector{Int}(undef, num_entries)

    for (i,j) in enumerate(entries)
        parsed_sequences[i] = string.(j["fasta_seq"])
        parsed_structures[i] = string.(j["dssp_ss3"])
        parsed_lengths[i] = length(parsed_sequences[i])
    end

    order_idx = sortperm(parsed_lengths)

    global train_sequences_ordered = parsed_sequences[order_idx]
    global train_structures_ordered = parsed_structures[order_idx]
    global train_lengths_ordered = parsed_lengths[order_idx]

end

## 2.2. Testing sequences, structures and lengths

testing_set_datapath = "../data/fold01testing_dataset.json"

open(testing_set_datapath, "r") do file

    entries = JSON.parse(file)
    num_entries = length(entries)

    parsed_sequences = Vector{Vector{String}}(undef, num_entries)
    parsed_structures = Vector{Vector{String}}(undef, num_entries)
    parsed_lengths = Vector{Int}(undef, num_entries)

    for (i,j) in enumerate(entries)
        parsed_sequences[i] = string.(j["fasta_seq"])
        parsed_structures[i] = string.(j["dssp_ss3"])
        parsed_lengths[i] = length(parsed_sequences[i])
    end

    order_idx = sortperm(parsed_lengths)

    global test_sequences_ordered = parsed_sequences[order_idx]
    global test_structures_ordered = parsed_structures[order_idx]
    global test_lengths_ordered = parsed_lengths[order_idx]

end

## 2.3. Lazy DataLoaders

training_dataset = Flux.Data.DataLoader((train_sequences_ordered, train_structures_ordered, train_lengths_ordered); batchsize=hyperpar["batch_size"], shuffle=true)
eval_train_dataset = Flux.Data.DataLoader((train_sequences_ordered, train_structures_ordered, train_lengths_ordered); batchsize=64,shuffle=false)
eval_test_dataset = Flux.Data.DataLoader((test_sequences_ordered, test_structures_ordered, test_lengths_ordered); batchsize=64,shuffle=false)

const train_acc_denominator = floor(sum(train_lengths_ordered)*0.96169) # Number of non-gap residues in the training samples. Adjust when using other dataset.
const test_acc_denominator = floor(sum(test_lengths_ordered)*0.96265) # Number of non-gap residues in the testing samples. Adjust when using other dataset.

#####
# 3. Network Structure
#####

model = Stack(
        @nntopo(x → enc → enc → enc → inc → inc → inc → inc → inc → inc → inc → inc → inc → inc → pe:(inc, pe) → bert_input → bert_output → bert_output → bert_output → y → y → y → y → y),
         # Embedding Layers
        Embed(config["embed_dimension"], config["input_vocab_size"]),
        unsqueeze(3),
        x -> permutedims(x, [3,2,1,4]),
        # Processing Layers
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        Inception8(config["embed_dimension"]),
        Dropout(config["bert_transformer_dropout_prob"]),
        x -> permutedims(x, [3,2,1,4]),
        x -> reshape(x, size(x,1), size(x,2), size(x,4)),
        PositionEmbedding(config["embed_dimension"], hyperpar["max_seq_len"]; trainable = false),
        (e, pe) -> e .+ pe,
        # Processing Layers
        Bert(
            config["embed_dimension"],
            config["bert_transformer_layer_heads"],
            config["bert_transformer_layer_hidden"],
            config["bert_num_hidden_layers"];
            act = config["bert_transformer_activation"],
            pdrop = config["bert_transformer_dropout_prob"],
            attn_pdrop = config["bert_transformer_dropout_prob"]
        ),
        # Detection Layers
        unsqueeze(3),
        bert_output -> permutedims(bert_output, [3,2,1,4]),
        Conv((1,5), config["embed_dimension"] => config["detection_num_filters"], sigmoid; pad = (0,2)),
        Conv((1,5), config["detection_num_filters"] => config["output_vocab_size"], identity; pad = (0,2)),
        x -> permutedims(x, [3,2,1,4]),
        x -> reshape(x, size(x,1), size(x,2), size(x,4)),
        logsoftmax
    ) |> gpu

#####
# 4. Model training
#####

smooth(et) = generic_smooth(et, hyperpar["smooth_epsilon"], config["input_vocab_size"])     # Creates an alias for smooth() obeying our hyperparameters
opt = Flux.Optimise.ADAM(hyperpar["learning_rate"])     # Initializes the ADAM optimizer
ps = params(model)      # Collects model parameters

#define training loop
function train!(hyperpar, training_dataset,
                model, ps, opt,
                smooth=smooth, train_loss=train_loss, eval_loss=eval_loss, accuracy_hits=accuracy_hits,
                eval_train_dataset=eval_train_dataset, eval_test_dataset=eval_test_dataset,
                fastaVocab=fastaVocab, ss3Vocab=ss3Vocab)
    @info "start training"
    for epoch in 1:hyperpar["number_of_epochs"]
    i = 0
        for batch in training_dataset
            i += 1
            x, y, l = batch

            train_x_processed = fastaVocab(preprocess(x, hyperpar["max_seq_len"])) |> gpu
            train_y_processed = ss3Vocab(preprocess(y, hyperpar["max_seq_len"])) |> gpu
            train_mask_processed = prepare_mask(y, hyperpar["max_seq_len"]) |> gpu

            grad = gradient(()->train_loss(model, train_x_processed, train_y_processed, train_mask_processed, smooth), ps)
            update!(opt, ps, grad)

            if i%64 == 0

                Flux.testmode!(model)

                let
                batch_loss_ideal_train_set = 0.0
                batch_loss_ideal_test_set = 0.0
                batch_hits_ideal_train_set = 0
                batch_hits_ideal_test_set = 0

                for eval_train_batch in eval_train_dataset

                    eval_x, eval_y, eval_length = eval_train_batch

                    eval_x_processed = fastaVocab(preprocess(eval_x, hyperpar["max_seq_len"])) |> gpu
                    eval_y_processed = ss3Vocab(preprocess(eval_y, hyperpar["max_seq_len"])) |> gpu
                    eval_mask_processed = prepare_mask(eval_y, hyperpar["max_seq_len"]) |> gpu

                    network_outputs = model(eval_x_processed)
                    network_predictions = Flux.onecold(network_outputs)

                    batch_loss_ideal_train_set += eval_loss(network_outputs, eval_y_processed, eval_mask_processed, smooth)
                    batch_hits_ideal_train_set += accuracy_hits(cpu(network_predictions), cpu(eval_y_processed), eval_length)

                end #end IDEAL_TRAIN_EVAL loop

                for eval_test_batch in eval_test_dataset

                    eval_x, eval_y, eval_length = eval_test_batch

                    eval_x_processed = fastaVocab(preprocess(eval_x, hyperpar["max_seq_len"])) |> gpu
                    eval_y_processed = ss3Vocab(preprocess(eval_y, hyperpar["max_seq_len"])) |> gpu
                    eval_mask_processed = prepare_mask(eval_y, hyperpar["max_seq_len"]) |> gpu

                    network_outputs = model(eval_x_processed)
                    network_predictions = Flux.onecold(network_outputs)

                    batch_loss_ideal_test_set += eval_loss(network_outputs, eval_y_processed, eval_mask_processed, smooth)
                    batch_hits_ideal_test_set += accuracy_hits(cpu(network_predictions), cpu(eval_y_processed), eval_length)

                end #end IDEAL_TEST_EVAL loop

            ideal_train_accuracy = batch_hits_ideal_train_set/train_acc_denominator
            ideal_test_accuracy = batch_hits_ideal_test_set/test_acc_denominator
            ideal_train_loss = batch_loss_ideal_train_set
            ideal_test_loss = batch_loss_ideal_test_set

            Flux.trainmode!(model)

            checkpoint_model = cpu(model)
            @save "TintiNet-IGBT-Classifier-Checkpoint-$(Dates.format(Dates.now(), dformat)).bson" checkpoint_model opt
            @printf("Epoch: %d\t Batch: %d\t IdTrLos: %2.6f\t IdTsLoss: %2.6f\t IdTrAcc: %2.6f\t IdTsAcc: %2.6f\t Time: %s\n",
                    epoch, i, ideal_train_loss, ideal_test_loss, ideal_train_accuracy, ideal_test_accuracy, Dates.format(now(), "yyyy-mm-dd HH:MM:SS"))

            end # end let

            end # endif

        end # end for batch

    end # end for epoch

end # end function train!()

#####
# Uncomment next line to allow the script to undergo training.
# train!(hyperpar, training_dataset, model, ps, opt)


text:
module CompressedArraysTests

using Test
using Gridap.Arrays
using FillArrays

values = [10,20,31]
ptrs = [1,2,3,3,2,2]
a = CompressedArray(values,ptrs)
r = values[ptrs]
test_array(a,r)

b = lazy_map(-,a)
test_array(b,-r)
@test isa(b,CompressedArray)

c = lazy_map(*,a,b)
test_array(c,r.*(-r))
@test isa(c,CompressedArray)

c = lazy_map(*,Int,a,b)
test_array(c,r.*(-r))
@test isa(c,CompressedArray)

b = Fill(4,length(ptrs))
c = lazy_map(*,a,b)
test_array(c,a.*b)
@test isa(c,CompressedArray)
@test c.ptrs === a.ptrs

c = lazy_map(*,b,a)
test_array(c,a.*b)
@test isa(c,CompressedArray)
@test c.ptrs === a.ptrs

k = CompressedArray([zero,+,-],copy(ptrs))
r = collect(CompressedArray([0,20,-31],ptrs))
c = lazy_map(evaluate,k,a)
test_array(c,r)
@test isa(c,CompressedArray)

k = CompressedArray([zero,+,-],ptrs)
r = collect(CompressedArray([0,20,-31],ptrs))
c = lazy_map(evaluate,k,a)
test_array(c,r)
@test isa(c,CompressedArray)

values = [10,20]
ptrs = [1,2,1,1,2,2]
b = CompressedArray(values,ptrs)

c = lazy_map(*,a,b)
test_array(c,a.*b)
@test ! isa(c,CompressedArray)

c = lazy_map(*,Int,a,b)
test_array(c,a.*b)
@test ! isa(c,CompressedArray)

c = lazy_map(*,Int,b,a)
test_array(c,a.*b)
@test ! isa(c,CompressedArray)

end # module


text:
module EvolvingGraphs

import Base: ==, show, issorted, deepcopy, length, eltype, push!, append!, reverse
using Requires

export

   # base
   Node, Edge, TimeNode, AttributeNode, TimeEdge, WeightedTimeEdge,
   make_node, node_index, node_key, edge_reverse, node_timestamp, edge_timestamp, find_node,
   node_attributes, edge_weight,
   AbstractPath, TemporalPath,

   # graph types
   AbstractGraph, AbstractEvolvingGraph, AbstractStaticGraph,
   EvolvingGraph, DiGraph, IntAdjacencyList, MatrixList,

   evolving_graph_from_arrays, add_bunch_of_edges!,
   adjacency_matrix, sparse_adjacency_matrix, block_adjacency_matrix,
   evolving_graph_to_adj, evolving_graph_to_matrices,

   # graph functions
   ## modify graph
   add_node!, add_edge!, add_graph!, rm_edge!, add_graph!, undirected!,

   has_edge, has_node, in_edges, in_degree, out_edges, out_degree, has_active_node,
   nodes, num_nodes, edges, num_edges,
   source, target, matrices, num_matrices, unique_timestamps,
   timestamps, num_timestamps, active_nodes, num_active_nodes,
   forward_neighbors, backward_neighbors, is_directed, undirected, attributes_values, aggregate_graph,
   random_graph, random_evolving_graph,

   # io
   egread, egwrite, load_graphml,

   # algorithms
   random_time_graph, random_evolving_graph, breadth_first_visit,

   # metric
   breadth_first_impl, depth_first_impl,
   shortest_temporal_path, shortest_temporal_distance,
   shortest_path, shortest_distance, temporal_efficiency, global_temporal_efficiency,

   # components
   temporal_connected, weak_connected, weak_connected_components,

   # sort slice
   issorted, sort_timestamps!, sort_timestamps, slice_timestamps!, slice_timestamps

include("base.jl")

# graph IO
include("read_write/evolving_graph_io.jl")
# need to import EzXML
@require EzXML include("read_write/graphml.jl")

# graph types
## static graphs
include("graph_types/static_dgraph.jl")

## evolving graphs
include("graph_types/evolving_graph.jl")

# Other types
include("graph_types/matrix_list.jl")
include("graph_types/adjacency_list.jl")
include("graph_types/incidence_list.jl")

include("show.jl")

# LightGraphs interface
@require LightGraphs include("read_write/lg_interface.jl")

# algorithms
include("algorithms/bfs.jl")
include("algorithms/dfs.jl")
include("algorithms/sort_slice.jl")
include("algorithms/components.jl")
include("algorithms/random.jl")

# metric
include("algorithms/shortest_distance.jl")
include("algorithms/shortest_temporal_distance.jl")

# Centrality module
include("centrality/centrality.jl")

end # module


text:
#---------------------------------------------------------
# # [PERK overview](@id 01-overview)
#---------------------------------------------------------

#=
This page illustrates the Julia package
[`PERK`](https://github.com/StevenWhitaker/PERK.jl).

This page was generated from a single Julia file:
[01-overview.jl](@__REPO_ROOT_URL__/01-overview.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`01-overview.ipynb`](@__NBVIEWER_ROOT_URL__/01-overview.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`01-overview.ipynb`](@__BINDER_ROOT_URL__/01-overview.ipynb).

# ### Setup

# Packages needed here.

using PERK
using MIRTjim: jim, prompt
using Plots; default(markerstrokecolor = :auto, label="")
using InteractiveUtils: versioninfo

#src The following line is helpful when running this file as a script;
#src this way it will prompt user to hit a key after each figure is displayed.

#src isinteractive() ? jim(:prompt, true) : prompt(:draw);

# ### Overview

#=
Although neural networks are very popular,
low-dimensional nonlinear regression problems
can be handled quite efficiently
by kernel ridge regression (KRR).
Training KRR does not require iterative algorithms
and is more interpretable than a deep network.
It is simply a nonlinear lifting
followed by ridge regression.
=#

#=
### Example

Here is an example of using KRR
to learn the function ``y = x^3``
from noisy training data.
=#

fun(x) = x^3
xtrain = LinRange(-1, 1, 101) * 3
ytrain = fun.(xtrain) + 1 * randn(size(xtrain))
p0 = scatter(xtrain, ytrain, xlabel="x", ylabel="y", label="training data")
plot!(p0, fun, label="y = x^3", legend=:top, color=:black)

# Here is the key training step
ρ = 1e-5
λ = 0.5
kernel = GaussianKernel(λ)
train = PERK.krr_train(ytrain, xtrain, kernel, ρ);
#src jim(train.K, "PERK K")

# Now examine the fit using (exhaustive) test data.
# The fit is very good within the range of the training data,
# and regresses to the mean outside of that range.
xtest = LinRange(-1, 1, 200) * 4
yhat = PERK.krr(xtest, train, kernel) # todo: remove kernel eventually
plot!(p0, xtest, yhat, label="KRR prediction", color=:magenta)

# The (only!) two parameters ρ and λ can be selected automatically
# using cross validation.
# todo: need example here

# ### Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer()
versioninfo(io)
split(String(take!(io)), '\n')

# And with the following package versions

import Pkg; Pkg.status()


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later
include("compile.jl")

# env var to force compilation from source, for testing purposes
const forcecompile = get(ENV, "FORCE_COMPILE_ZMQ", "no") == "yes"

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS || forcecompile
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libzmq"], :libzmq),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/JuliaInterop/ZMQBuilder/releases/download/v4.2.5+6"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/ZMQ.aarch64-linux-gnu.tar.gz", "ce637adc7f8bbd0a628d1d468e790236b64dc0f41e4c3ead991842a564b413b9"),
    Linux(:aarch64, :musl) => ("$bin_prefix/ZMQ.aarch64-linux-musl.tar.gz", "13a1695fa47fb4d7ffdad55dfcd0fa4a66aad16235331aa7af4995be83f4e556"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/ZMQ.arm-linux-gnueabihf.tar.gz", "c7050cc3eb8e1ad81b9eedbb41a4b2a89383d698dd4931cc1b90239c3640ffe5"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/ZMQ.arm-linux-musleabihf.tar.gz", "e80977a6e6695846de6553448aa9d9d6482bfa52e5fcf536e227162e16322fb8"),
    Linux(:i686, :glibc) => ("$bin_prefix/ZMQ.i686-linux-gnu.tar.gz", "18f043e1179a2db4dc2867c9da73cbccf1b2fc6851280b07891898951fd8a950"),
    Linux(:i686, :musl) => ("$bin_prefix/ZMQ.i686-linux-musl.tar.gz", "9f2425ec361ecadc0e836ea8cb7f378839ae5dcef8205a85d42f513368f2b6a7"),
    Windows(:i686) => ("$bin_prefix/ZMQ.i686-w64-mingw32.tar.gz", "960cd3ee0443389e35614dbd5954f68fe969abc42b9b5d26c0339bf449854473"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/ZMQ.powerpc64le-linux-gnu.tar.gz", "51eb2313e6552a7d2eaf4bb2167b3498fbc16dfbab070bfb9a0a3788a2e64f6d"),
    MacOS(:x86_64) => ("$bin_prefix/ZMQ.x86_64-apple-darwin14.tar.gz", "ee0aa673e4fdb1574ac2d6c7c3a22c3c03f8ca75928968f051f78a2c32941b20"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/ZMQ.x86_64-linux-gnu.tar.gz", "59266a7b62b74da29389e058fe014bef31ff47ed83c021a6d0fd59053f428cd3"),
    Linux(:x86_64, :musl) => ("$bin_prefix/ZMQ.x86_64-linux-musl.tar.gz", "7a88e5c5f933dd1b67b23ae70a4c78e27324170140f0f42a355b5ae53a8794c4"),
    FreeBSD(:x86_64) => ("$bin_prefix/ZMQ.x86_64-unknown-freebsd11.1.tar.gz", "34c1c5077c114417b893ae8048b71ed3b7812875b9c789b33407e56fed5e3646"),
    Windows(:x86_64) => ("$bin_prefix/ZMQ.x86_64-w64-mingw32.tar.gz", "ba1e0394919d6e6eeda13e201278f754069d33967ee96043cb52505e05a9128b"),
)

# source code tarball and hash for fallback compilation
source_url = "https://github.com/zeromq/libzmq/releases/download/v4.2.5/zeromq-4.2.5.tar.gz"
source_hash = "cc9090ba35713d59bb2f7d7965f877036c49c5558ea0c290b0dcc6f2a17e489f"

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key()) && !forcecompile
    url, tarball_hash = download_info[platform_key()]
    if !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)

        # check again whether the dependency is satisfied, which
        # may not be true if dlopen fails due to a libc++ incompatibility (#176)
        unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
    end
end

if unsatisfied || forcecompile
    # Fall back to building from source, giving the library a different name
    # so that it is not overwritten by BinaryBuilder downloads or vice-versa.
    libname = "libzmq_from_source"
    products = [ LibraryProduct(prefix, [libname], :libzmq) ]
    source_path = joinpath(prefix, "downloads", "src.tar.gz")
    if !isfile(source_path) || !verify(source_path, source_hash; verbose=verbose) || !satisfied(products[1]; verbose=verbose)
        compile(libname, source_url, source_hash, prefix=prefix, verbose=verbose)
    end
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
"""
	makecpt(cmd0::String="", arg1=[]; kwargs...)

Make static color palette tables (CPTs).

Full option list at [`makecpt`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html)

- **A** : **alpha** : **transparency** : -- Str --

    Sets a constant level of transparency (0-100) for all color slices.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#a)
- $(GMT.opt_C)
- **D** : -- Str or [] --			Flags = [i|o]

    Select the back- and foreground colors to match the colors for lowest and highest
    z-values in the output CPT.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#d)
- **E** : **data_levels** : -- Int or [] --		Flags = [nlevels]

    Implies reading data table(s) from file or arrays. We use the last data column to
    determine the data range
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#e)
- **F** : **force_rgb** : -- Str or [] --		Flags = [R|r|h|c][+c]]

    Force output CPT to written with r/g/b codes, gray-scale values or color name.
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#f)
- **G** : **truncate** : -- Str --              Flags = zlo/zhi

    Truncate the incoming CPT so that the lowest and highest z-levels are to zlo and zhi.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#g)
- **I** : **inverse** : **reverse** : -- Str --	Flags = [c][z]

    Reverse the sense of color progression in the master CPT.
    [`-I`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#i)
- **M** : **overrule_bg** -- Bool or [] --

    Overrule background, foreground, and NaN colors specified in the master CPT with the values of
    the parameters COLOR_BACKGROUND, COLOR_FOREGROUND, and COLOR_NAN.
    [`-M`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#m)
- **N** : **no_bg** -- Bool or [] --

    Do not write out the background, foreground, and NaN-color fields.
- **Q** : **log** : -- Bool or [] or Str --			Flags = [i|o]

    Selects a logarithmic interpolation scheme [Default is linear].
    [`-Q`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#q)
- **S** : **auto** : -- Bool or [] or Str --			Flags = [mode]

    Determine a suitable range for the -T option from the input table(s) (or stdin).
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#s)
- **T** : **range** : -- Str --			Flags = [min/max/inc[+b|l|n]|file|list]

    Defines the range of the new CPT by giving the lowest and highest z-value and interval.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#t)
- **W** : **wrap** : **categorical** : -- Bool or Str or [] --      Flags = [w]

    Do not interpolate the input color table but pick the output colors starting at the
    beginning of the color table, until colors for all intervals are assigned.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#w)
- **Z** : **continuous** : -- Bool or [] --

    Creates a continuous CPT [Default is discontinuous, i.e., constant colors for each interval].
    [`-Z`](http://gmt.soest.hawaii.edu/doc/latest/makecpt.html#z)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_i)
"""
function makecpt(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("makecpt", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd = parse_V_params("", d)
	cmd, opt_bi = parse_bi(cmd, d)
	cmd, opt_di = parse_bi(cmd, d)
	cmd, opt_i = parse_i(cmd, d)

	# If file name sent in, read it and compute a tight -R if this was not provided
	cmd, arg1, opt_R, = read_data(d, cmd0, cmd, arg1, " ", opt_i, opt_bi, opt_di)
	cmd, arg1, arg2, = add_opt_cpt(d, cmd, [:C :color :cmap], 'C', 0, arg1, [])

	for sym in [:E :data_levels]
		if (haskey(d, sym))
			if (isempty_(arg1) && isempty_(data))
				error("E option requires that a data table is provided as well")
			else
				cmd = cmd * " -E" * arg2str(d[sym])
			end
			break
		end
	end

	cmd = add_opt(cmd, 'A', d, [:A :alpha :transparency])
	cmd = add_opt(cmd, 'D', d, [:D :bg :background])
	cmd = add_opt(cmd, 'F', d, [:F :force_rgb])
	cmd = add_opt(cmd, 'G', d, [:G :truncate])
	cmd = add_opt(cmd, 'I', d, [:I :inverse :reverse])
	cmd = add_opt(cmd, 'M', d, [:M :overrule_bg])
	cmd = add_opt(cmd, 'N', d, [:N :no_bg])
	cmd = add_opt(cmd, 'Q', d, [:Q :log])
	cmd = add_opt(cmd, 'S', d, [:S :auto])
	cmd = add_opt(cmd, 'T', d, [:T :range])
	cmd = add_opt(cmd, 'W', d, [:W :wrap :categorical])
	cmd = add_opt(cmd, 'Z', d, [:Z :continuous])

	if (haskey(d, :cptname))
		cmd = cmd * " > " * d[:cptname]
		C = gmt("makecpt " * cmd)
		(haskey(d, :Vd)) && println(@sprintf("\tmakecpt %s", cmd))
	else
		(haskey(d, :Vd)) && println(@sprintf("\tmakecpt %s", cmd))
		if (isempty_(arg1))
			C = gmt("makecpt " * cmd)
		else
			C = gmt("makecpt " * cmd, arg1)
		end
	end
end

# ---------------------------------------------------------------------------------------------------
# Version to use with the -E option
#makecpt(arg1=[]; kw...) = makecpt("", arg1; kw...)


text:
using MocosSim
using Test
using Random

import MocosSim: time, numdead, numdetected

tests = [
  "age_coupling",
  "matrix_alias_sampler",
  "population_grouping",
  "infection_modulations",
  "household_grouping",
]

if length(ARGS) > 0
    tests = ARGS
end

@testset "MocosSim" begin

for t in tests
    fp = joinpath(dirname(@__FILE__), "test_$t.jl")
    println("$fp ...")
    include(fp)
end

end # @testset


text:
using StaticArrays

# computes the MV energy
# From MV: Omega = sum_n <r2>n - |<r>n|^2
# <r>n = -1/N sum_k,b wb b Im ln(Mnn,kb)
# <r2>n = 1/N sum_k,b wb [(1 - |Mnn,kb|^2) + (Im ln(Mnn,kb))^2]
struct Omega_res
    Ωtot :: Float64
    ΩI :: Float64
    ΩOD :: Float64
    ΩD :: Float64 #not implemented yet
    Ωtilde :: Float64
    frozen_weight :: Float64
    spreads :: Array{Float64,1} #nband
    centers :: Array{Float64,2} #3 x nband
    gradient :: Array{ComplexF64, 5}#nband x nband x n1 x n2 x n3
end

imaglog(z) = atan(imag(z), real(z))

@views function omega(p,A,compute_grad=false,only_r2=false)
    mu=0.
    nfrozen=0 #keep in case we want to do this later on
    if compute_grad
        if !only_r2
            centers = omega(p,A,false).centers
        end
    end
    grad = zeros(ComplexF64,p.nband,p.nwannier,p.N1,p.N2,p.N3)
    r = zeros(3,p.nwannier)
    r2 = zeros(p.nwannier)
    ΩI = 0.
    ΩOD = 0.
    ΩD = 0.
    frozen_weight = 0.
    R = zeros(ComplexF64,p.nband,p.nwannier)
    T = zeros(ComplexF64,p.nband,p.nwannier)
    b = zeros(3)
    Mkb = zeros(ComplexF64,p.nwannier,p.nwannier)
    Okb = zeros(ComplexF64,p.nband,p.nwannier)
    for i=1:p.N1,j=1:p.N2,k=1:p.N3
        K = p.ijk_to_K[i,j,k]

        frozen_weight -= mu*sum(abs2, A[1:nfrozen,:,i,j,k])
        if compute_grad
            grad[1:nfrozen,:,i,j,k] = -2*mu*A[1:nfrozen,:,i,j,k]
        end

        for ib = 1:p.nntot
            neighbor = p.neighbors[K,ib]
            i_n,j_n,k_n = p.K_to_ijk[neighbor,:]
            # Mkb = overlap_A([i,j,k],[i_n,j_n,k_n],p)
            Okb .= overlap(K,neighbor,p)*A[:,:,i_n,j_n,k_n]
            # @assert overlap(K,neighbor,p)' ≈ overlap(neighbor,K,p) #compute-intensive, but should be true
            Mkb .= A[:,:,i,j,k]'*Okb
            b .= p.recp_unit_cell*([p.t1[i_n];p.t2[j_n];p.t3[k_n]] .+ p.displ_vecs[:,K,ib] .- [p.t1[i];p.t2[j];p.t3[k]])

            if compute_grad
                # #MV way
                # A(B) = (B-B')/2
                # S(B) = (B+B')/(2*im)
                # q = imaglog.(diag(Mkb)) + centers'*b
                # for m=1:p.nwannier,n=1:p.nwannier
                #     R[m,n] = Mkb[m,n]*conj(Mkb[n,n])
                #     T[m,n] = Mkb[m,n]/Mkb[n,n]*q[n]
                # end
                # grad[i,j,k,:,:] += 4*p.wb*(A(R) .- S(T))

                q = imaglog.(diag(Mkb))
                if !only_r2
                    q += centers'*b
                end
                for n=1:p.nwannier
                    if abs(Mkb[n,n]) < 1e-10 # error if division by zero. Should not happen if the initial gauge is not too bad
                        println("Mkbnn too large! $j $k -> $(j_n) $(k_n)")
                        display(Mkb)
                        error()
                    end
                    @assert abs(Mkb[n,n]) > 1e-10
                    Tfac = -im*q[n]/Mkb[n,n]
                    for m=1:p.nband
                        R[m,n] = -Okb[m,n]*conj(Mkb[n,n])
                        # T[m,n] = -im*Okb[m,n]/(Mkb[n,n])*q[n]
                        T[m,n] = Tfac*Okb[m,n]
                    end
                end
                grad[:,:,i,j,k] .+= 4*p.wb.*(R.+T)
            end

            ΩI += p.wb*(p.nwannier - sum(abs2,Mkb))
            ΩOD += p.wb*sum(abs2,Mkb .- diagm(0=>diag(Mkb)))
            for n=1:p.nwannier
                if !only_r2
                    r[:,n] -= p.wb*imaglog(Mkb[n,n])*b
                end
                r2[n] += p.wb*(1-abs(Mkb[n,n])^2 + imaglog(Mkb[n,n])^2)
                # r2[n] += p.wb*2*(1 - real(Mkb[n,n]))
            end
        end
    end
    Ntot = (p.N1*p.N2*p.N3)
    r /= Ntot
    r2 /= Ntot
    ΩI /= Ntot
    ΩOD /= Ntot
    ΩD /= Ntot
    frozen_weight /= Ntot
    (grad /= Ntot)

    spreads = (r2 - dropdims(sum(abs.(r).^2,dims=1),dims=1))
    Ωtot = sum(spreads) + frozen_weight
    Ωtilde = Ωtot - ΩI
    return Omega_res(Ωtot,ΩI,ΩOD,ΩD,Ωtilde,frozen_weight,spreads,r,grad)
end

# local part of the contribution to r^2
function omega_loc(p,A)
    r = zeros(3,p.nwannier)
    r2 = zeros(p.nwannier)
    loc = zeros(Float64,p.N1,p.N2,p.N3)
    for i=1:p.N1,j=1:p.N2,k=1:p.N3
        K = p.ijk_to_K[i,j,k]
        for ib = 1:p.nntot
            neighbor = p.neighbors[K,ib]
            i_n,j_n,k_n = p.K_to_ijk[neighbor,:]
            # Mkb = overlap_A([i,j,k],[i_n,j_n,k_n],p)
            Okb = overlap(K,neighbor,p)*A[:,:,i_n,j_n,k_n]
            Mkb = A[:,:,i,j,k]'*Okb
            for n=1:p.nwannier
                loc[i,j,k] += p.wb*(1-abs(Mkb[n,n])^2 + imaglog(Mkb[n,n])^2)
            end
        end
    end
    return loc
end

max_stars_repo_name: antoine-levitt/wannier

max_issues_repo_name: antoine-levitt/wannier

max_forks_repo_name: antoine-levitt/wannier


text:
export
    Separable,
    SepCPTs

"SepCPTs = Array{Dict{I, Array{Float64, 1}} where I}"
SepCPTs = Array{Dict{I, Array{Float64, 1}} where I}

# TODO: Generalize this to separable additive decompositions in general, not just CPTs
"""
    function Separable(range::Vector{O}, probabilities :: Vector{Float64}, compparams :: SepCPTs) where O

Constructs an sfunc representing separable models, defined by additive decompositon of a conditional probability distribution into
separate distributions depending on each of the parents.

`Separable`s are implemented as a `Mixture` of `Extend` sfuncs that extend `DiscreteCPT`s.

To construct a `Separable`, this method is passed the `range` of output values, the `probabilities`
of each of the underlying `DiscreteCPT` (which are the internal sfuncs of the `Mixture`), and
the parameters for each of the `DiscreteCPT`s.  For example,

```
alphas = [0.2, 0.3, 0.5]
cpd1 = Dict((1,) => [0.1, 0.9], (2,) => [0.2, 0.8])
cpd2 = Dict((1,) => [0.3, 0.7], (2,) => [0.4, 0.6], (3,) => [0.5, 0.5])
cpd3 = Dict((1,) => [0.6, 0.4], (2,) => [0.7, 0.3])
cpds :: Array{Dict{I,Array{Float64,1}} where I,1} = [cpd1, cpd2, cpd3]
s = Separable([1, 2], alphas, cpds)
```

See also:  [`Mixture`](@ref), [`Extend`](@ref), [`DiscreteCPT`](@ref), [`Table`](@ref)
"""
function Separable(range::Vector{O}, probabilities :: Vector{Float64}, compparams :: SepCPTs) where O
    N = length(compparams)
    @assert length(probabilities) == N

    IS = []
    for i = 1:N
        ks = collect(keys(compparams[i]))
        push!(IS, typeof(ks[1][1]))
    end
    J = Tuple{IS...}

    # Explicit typing is necessary to ensure that the cpts passed to Mixture all have the same type.
    function make_cpt(i,I)::SFunc{J,O}
        cpt::Table{1, I, Tuple{}, I, O, <:Dist{O}} =
           DiscreteCPT(range, compparams[i])
        extended::SFunc{J,O} = Extend(J, cpt, i)
        return extended
    end

    cpts::Vector{SFunc{J,O}} =
        [make_cpt(i,Tuple{IS[i]}) for i in 1:N]
    return Mixture(cpts, probabilities)
end


text:
import FractionalTimeDG: Store, ODEdG!, reconstruction_pts,
                         evaluate_pcwise_poly!, jumps, dG_error_estimator

using PyPlot
#using Polynomials: polyval

T = Float64
#T = BigFloat

λ = 1 / parse(T, "2")
u0 = one(T)
ampl = one(T)
f(t) = ampl * cos(π*t)
c =  λ^2 + big(π)^2
u(t) = ( ( u0 - ampl * λ / c ) * exp(-λ*t)
        + ampl * ( λ * cos(π*t) + π * sin(π*t) ) / c )
tmax = parse(T, "2")

N = 4
r = 3
M = 4
pts_per_interval = 50
store = Store(one(T), r, M, pts_per_interval)
t, U = ODEdG!(λ, tmax, f, u0, N, r, M)
pcwise_t, pcwise_U = evaluate_pcwise_poly!(U, t, pts_per_interval, store)
pcwise_u = T[ u(t) for t in pcwise_t ]
figure(1)
line1 = plot(pcwise_t, pcwise_u, "k:")
line2 = plot(pcwise_t, pcwise_U, "c")
legend((line1[1], line2[1]), (L"$u$", L"$U$"))
xlabel(L"$t$", size=14)
ylabel(L"$U$", size=14)
xticks((0.0, 1/2, 1, 3/2, 2), ("0.0", "0.5", "1.0", "1.5", "2.0"))
grid(true)

figure(2)
JU = jumps(U, t, u0)
pcwise_tx, pcwise_approx_err = dG_error_estimator(JU, t, r, pts_per_interval)
line1 = plot(pcwise_t, pcwise_U-pcwise_u, "C0")
line2 = plot(pcwise_tx, pcwise_approx_err, "C1:")
xticks((0.0, 1/2, 1, 3/2, 2), ("0.0", "0.5", "1.0", "1.5", "2.0"))
legend((line1[1], line2[1]), (L"$U-u$", "error estimator"))
grid(true)


text:
typealias TCDict Dict{Symbol, TestCase}

# put test cases into buckets
tc_all = TCDict()
tc_nonstiff = TCDict()
tc_stiff = TCDict()
tc_dae = TCDict()
tc_dae1 = TCDict() # index 1
tc_dae2 = TCDict() # index 2
tc_pde = TCDict() # stemming from discretized PDEs
tc_ide = TCDict() # implicit differential equations

# # problems with implicit equation
# tc_impl_stiff = TCDict()
# tc_impl_dae = TCDict()

# # non-stiff ODEs
include("threebody.jl")
include("plei.jl")

# # stiff odes
include("hires.jl")
include("rober.jl")
include("vdpol.jl")
include("bruss1d.jl")

# # DAEs index 1
include("chemakzo.jl")


text:
using NeuralProcesses
using NeuralProcesses.Experiment
using Distributions
using FiniteDifferences
using Flux
using LinearAlgebra
using NNlib
using StatsFuns
using Stheno
using Test
using Tracker

@testset "NeuralProcesses.jl" begin
    include("util.jl")
    include("discretisation.jl")
    include("setconv.jl")
    include("attention.jl")
    include("data.jl")
    include("model.jl")
end


text:
@testset "ModelFrame (legacy API)" begin

    @testset "#133" begin
        df = (x = rand(12),
              y = categorical(repeat(1:3, inner=4)),
              z = categorical(repeat(1:2, outer=6)));
        f = @formula(x ~ y * z);
        mf = ModelFrame(f, df)
        mm = ModelMatrix(mf)
        @test mm.assign == [1, 2, 2, 3, 4, 4]

        # `q` isn't in df  - should throw ArgumentError
        f2 = @formula(x ~ y * q)
        @test_throws ArgumentError ModelFrame(f2, df)

    end

end


text:
"""
    pie(fractions; kwargs...)

Creates a pie chart with the given `fractions`.

## Attributes
$(ATTRIBUTES)
"""
@recipe(Pie, values) do scene
    Theme(
        normalize = true,
        color = :gray,
        strokecolor = :black,
        strokewidth = 1,
        vertex_per_deg = 1,
        radius = 1,
        inner_radius = 0,
        offset = 0,
        inspectable = theme(scene, :inspectable)
    )
end

function plot!(plot::Pie)

    values = plot[1]

    polys = lift(values, plot.vertex_per_deg, plot.radius, plot.inner_radius, plot.offset, plot.normalize) do vals, vertex_per_deg, radius, inner_radius, offset, normalize

        T = eltype(vals)

        # find start and end angles of all pie pieces
        summed = cumsum([zero(T); vals])
        boundaries = if normalize
            summed ./ summed[end] .* 2pi
        else
            summed
        end

        # create vector of a vector of points for each piece
        vertex_arrays = map(boundaries[1:end-1], boundaries[2:end]) do sta, en
            distance = en - sta
            # how many vertices are needed for the curve?
            nvertices = max(2, ceil(Int, rad2deg(distance) * vertex_per_deg))

            # curve points
            points = map(LinRange(sta, en, nvertices)) do rad
                Point2f0(cos(rad + offset), sin(rad + offset)) .* radius
            end

            # add inner points (either curve or one point)
            if inner_radius != 0
                inner_points = map(LinRange(en, sta, nvertices)) do rad
                    Point2f0(cos(rad + offset), sin(rad + offset)) .* inner_radius
                end

                append!(points, inner_points)
            else
                push!(points, Point2f0(0, 0))
            end

            points
        end
    end

    # plot pieces as polys
    poly!(plot, polys, color = plot.color, strokewidth = plot.strokewidth,
        strokecolor = plot.strokecolor, inspectable = plot.inspectable)

    plot
end


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el)
        el
    end
end

# ╔═╡ 40ce97fc-aed8-11ec-381d-39c142fc0c6b
begin
	using Pkg
	using Revise
	Pkg.activate(Base.current_project(@__DIR__))
	using Portinari, PlutoUI, HypertextLiteral, PlutoDevMacros
end

# ╔═╡ a8a51537-c5c9-4f88-8c0b-3eb8d5a3c85d
@bind sz Slider(1:100; default=10)

# ╔═╡ 27d8844f-48fd-4daa-a9fd-97fb1ee088f7
x = rand(sz)

# ╔═╡ 0c9f5641-b917-422d-abd6-ae3b3d1ed0cf
y = rand(sz)

# ╔═╡ d2a0eb54-5962-4df2-92dd-b47f3cf4b0c1
y2 = [v^2 for v ∈ y]

# ╔═╡ 6abb5eef-6400-48b2-b526-f5e5a3ae44af
@bind lineevents line = Line(x, y, "lineid"; d3Attributes=D3Attributes(
	style=Dict(
		"fill" => "none",
		"stroke" => "green",
		"stroke-width" => "6px"
	),
	events=["click"]
),curveType=Portinari.BasisClosed)

# ╔═╡ a19dd476-b70e-4847-9f40-15f9f5715262
@bind events circles = Shape(x, y, y2 .* 1000, "shapes"; d3Attributes=D3Attributes(attributes=Dict("fill" => "rgba(0, 100, 255, 0.3)"), events=["click", "mouseover"]))

# ╔═╡ 3984e4f8-a172-48dc-a135-5ad8d8c6244d
md"""
# Report
###### X: $(events["click"]["data"]["x"])
###### Y: $(events["click"]["data"]["y"])
###### Size: $(events["click"]["data"]["size"])
###### Mouseover size! $(events["mouseover"]["data"]["size"])
"""

# ╔═╡ a4ac8702-c21e-41b0-8d60-ebff0b680214
axis = Axis2D(Axis(x, 400, 50, Portinari.Bottom; show=true), Axis(y, 500, 50, Portinari.Right; show=true), [line, circles], D3Attributes(style=Dict("width" => "600px", "height" => "600px")), [Portinari.zoom], "id")

# ╔═╡ a5648729-28cf-4296-9490-2f07d50720e5
D3Canvas([axis, line], "canv"; d3Attributes=D3Attributes(style=Dict("width" => "600px", "height" => "600px"), attributes=Dict("viewBox" => "0 0 600 600")))

# ╔═╡ Cell order:
# ╠═40ce97fc-aed8-11ec-381d-39c142fc0c6b
# ╠═27d8844f-48fd-4daa-a9fd-97fb1ee088f7
# ╠═0c9f5641-b917-422d-abd6-ae3b3d1ed0cf
# ╟─d2a0eb54-5962-4df2-92dd-b47f3cf4b0c1
# ╟─6abb5eef-6400-48b2-b526-f5e5a3ae44af
# ╠═a19dd476-b70e-4847-9f40-15f9f5715262
# ╟─3984e4f8-a172-48dc-a135-5ad8d8c6244d
# ╟─a8a51537-c5c9-4f88-8c0b-3eb8d5a3c85d
# ╠═a5648729-28cf-4296-9490-2f07d50720e5
# ╠═a4ac8702-c21e-41b0-8d60-ebff0b680214


text:
# JuMP instance tests

module TestJuMP

using Test
using LinearAlgebra
import MathOptInterface
const MOI = MathOptInterface
import JuMP
import MOIPajarito
import PajaritoExtras
import PajaritoExtras: Prim, Dual, svec_idx
import Hypatia
import Hypatia.Cones: vec_length, svec_length, vec_copyto!
import Hypatia.PolyUtils

include("model_utils.jl")

# all instances
inst_all = String[
    "possemideftri1",
    "possemideftri2",
    "possemideftrisparse1",
    "possemideftrisparse2",
    "epinormeucl1",
    "epinormeucl2",
    "epipersquare1",
    "epipersquare2",
    "epinorminf1",
    "epinorminf2",
    "epinormspectraltri1",
    "epinormspectraltri2",
    "epinormspectral1",
    "epinormspectral2",
    "hypogeomean1",
    "hypogeomean2",
    "hyporootdettri1",
    "hyporootdettri2",
    "vector_epipersepspectral1",
    "vector_epipersepspectral2",
    "matrix_epipersepspectral1",
    "matrix_epipersepspectral2",
    "wsosinterpnonnegative1",
    "wsosinterpnonnegative2",
    "specialorderedset",
]

# instances to test with extended formulations option off
inst_noextend = String[
    "epinormeucl1",
    "epinormeucl2",
    "epipersquare1",
    "epipersquare2",
    "epinorminf1",
    "epinorminf2",
    "hypogeomean1",
    "hypogeomean2",
    "vector_epipersepspectral1",
    "vector_epipersepspectral2",
]

function runtests(oa_solver)
    @testset "solve conic" begin
        @testset "iterative" begin
            @testset "all" begin
                run_jump_tests(inst_all, true, true, true, oa_solver)
            end
            @testset "no ext" begin
                run_jump_tests(inst_noextend, false, true, true, oa_solver)
            end
        end
        @testset "one tree" begin
            @testset "all" begin
                run_jump_tests(inst_all, true, false, true, oa_solver)
            end
            @testset "no ext" begin
                run_jump_tests(inst_noextend, false, false, true, oa_solver)
            end
        end
    end
    @testset "sep only" begin
        @testset "iterative" begin
            @testset "all" begin
                run_jump_tests(inst_all, true, true, false, oa_solver)
            end
            @testset "no ext" begin
                run_jump_tests(inst_noextend, false, true, false, oa_solver)
            end
        end
        @testset "one tree" begin
            @testset "all" begin
                run_jump_tests(inst_all, true, false, false, oa_solver)
            end
            @testset "no ext" begin
                run_jump_tests(inst_noextend, false, false, false, oa_solver)
            end
        end
    end
    return
end

function run_jump_tests(
    insts::Vector{String},
    use_extended_form::Bool,
    use_iterative_method::Bool,
    solve_relax_subp::Bool,
    oa_solver,
)
    hypatia = MOI.OptimizerWithAttributes(
        Hypatia.Optimizer,
        MOI.Silent() => true,
        "near_factor" => 1000,
        "tol_feas" => 1e-10,
        "tol_infeas" => 1e-12,
        "tol_rel_opt" => 1e-9,
        "tol_abs_opt" => 1e-8,
        "tol_illposed" => 1e-9,
        "tol_slow" => 2e-2,
        "tol_inconsistent" => 1e-7,
        "tol_tau_infeas" => 1e-7, # TODO tune
    )

    opt = JuMP.optimizer_with_attributes(
        MOIPajarito.Optimizer,
        "verbose" => true,
        "oa_solver" => oa_solver,
        "conic_solver" => hypatia,
        "sep_solver" => hypatia,
        "use_extended_form" => use_extended_form,
        "use_iterative_method" => use_iterative_method,
        "solve_relaxation" => solve_relax_subp,
        "solve_subproblems" => solve_relax_subp,
        "iteration_limit" => 200,
        # "time_limit" => 120.0,
    )

    @testset "$inst" for inst in insts
        println()
        @info inst
        eval(Symbol(inst))(opt)
    end
    return
end

include("JuMP_instances.jl")

end


text:
abstract type RecordingLayer{Next <: Layer} <: Layer{Next} end

before(::Type{<:RecordingLayer}, storage, path) = nothing

function HTTP.request(::Type{RecordingLayer{Next}}, resp) where Next
    state = get_state()
    push!(state.responses, deepcopy(resp))
    return request(Next, resp)
end

function after(::Type{<:RecordingLayer}, storage, path)
    state = get_state()
    for resp in state.responses
        filter!(drop_keys(state.ignore_headers), resp.request.headers)
        filter!(drop_keys(state.ignore_headers), resp.headers)
        resp.request.target = filter_query(resp.request, state.ignore_query)
    end
    mkpath(dirname(path))
    store(storage, path; responses=state.responses, format=FORMAT)
end

function filter_query(request, ignore)
    uri = URI(request.target)
    query = filter(drop_keys(ignore), parse_query(uri))
    return if isempty(query)
        uri.path
    else
        uri.path * '?' * join((join(p, '=') for p in query), '&')
    end
end


text:
using Test, SimpleGraphs, SimpleGraphRepresentations

@test 1==1


text:
export share
export share_replay
export share_recent

"""
    share(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

`share()` is a shortcut for `publish() + ref_count()`

See also: [`AbstractOperator`](@ref), [`multicast`](@ref), [`publish`](@ref), [`ref_count`](@ref)
"""
share(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler } = publish(scheduler = scheduler) + ref_count()

"""
    share_replay(size::Int; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

`share_replay(size)` is a shortcut for `publish_replay(size) + ref_count()`

See also: [`AbstractOperator`](@ref), [`multicast`](@ref), [`publish`](@ref), [`publish_replay`](@ref), [`ref_count`](@ref)
"""
share_replay(size::Int; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler } = publish_replay(size, scheduler = scheduler) + ref_count()

"""
    share_recent(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler }

`share_recent()` is a shortcut for `publish_recent() + ref_count()`

See also: [`AbstractOperator`](@ref), [`multicast`](@ref), [`publish`](@ref), [`publish_recent`](@ref), [`ref_count`](@ref)
"""
share_recent(; scheduler::H = AsapScheduler()) where { H <: AbstractScheduler } = publish_recent(scheduler = scheduler) + ref_count()


text:
const sqrt3 = sqrt(3.)

# upper bound on distance between x and y in the uncertainty region centered at R, T
function upperdist(x, y, R, T)
    return norm(R*x .+ T .- y)
end

# greatest distance that x can be moved within the uncertainty region
# this can be used as a less expensive, but looser, calculation of the lowerbound; ldist = uppderdist - uncertaintydist
function uncertaintydist(x, σᵣ, σₜ)
    γᵣ = 2*sin(min(sqrt3*σᵣ, π/2))*norm(x)
    γₜ = sqrt3 * σₜ
    return γᵣ + γₜ
end

function lowerdist(x, y, R::AngleAxis, T, σᵣ, σₜ)
    # return Inf for bounds if the rotation lies outside the π-sphere
    if R.theta > π
        inf = typemax(promote_type(eltype(x), eltype(y)))
        return inf, inf
    end

    # magnitudes of x and (y-T)
    xnorm, ynorm = norm(x), norm(y-T)

    # α is the angle between R*x (the initial rotation), the origin, and y-T
    if xnorm*ynorm == 0
        cosα = one(promote_type(eltype(x), eltype(y)))
    else
        cosα = dot(R*x, y-T)/(xnorm*ynorm)
    end

    # β is the largest allowed angle of rotation
    cosβ = cos(min(sqrt3*σᵣ, π))

    # if the α > β, then the two points can be rotated within the uncertainty region to be colinear with the origin
    if cosα >= cosβ
        lbdist = max(abs(xnorm-ynorm) - sqrt3*σₜ, 0)
    # otherwise, use the law of cosines to calculate the minimum distance
    else
        lbdist = try max(√(xnorm^2 + ynorm^2 - 2*xnorm*ynorm*(cosα*cosβ+√((1-cosα^2)*(1-cosβ^2)))) - sqrt3*σₜ, 0)  # law of cosines
        catch e     # when the argument for the square root is negative (within machine precision of 0, usually)
            zero(promote_type(eltype(x), eltype(y)))
        end
    end

    # return the lowerbound on distance between x and y within the uncertainty region
    return lbdist
end

# these let you modify the objective function, which takes x, y, and dist as arguments
upperbound(x, y, R, T, objective=p->p[3]) = objective(x,y,upperdist(x,y,R,T))
lowerbound(x, y, R, T, σᵣ, σₜ, objective=p->p[3]) = objective(x,y,lowerdist(x,y,R,T,σᵣ,σₜ))