text:
"""
    RecursiveApply

This module contains operators to recurse over nested `Tuple`s or `NamedTuple`s.

To extend to another type `T`, define `RecursiveApply.rmap(fn, args::T...)`
"""
module RecursiveApply

export ⊞, ⊠, ⊟, tuplemap

"""
    tuplemap(fn::Function, tup)

A `map` impl for mapping function `fn` a tuple argument `tup`

Currently just calls `Base.map` behind the scenes but is left
stubbed out for potential specialization in the future.
"""
@inline function tuplemap(fn::F, tup::Tuple) where {F}
    map(fn, tup)
end

"""
    tuplemap(fn::Function, tup1, tup2)

A `map` impl for mapping function `fn` over `tup1`, `tup2` tuple arguments.

Currently just calls `Base.map` behind the scenes but is left
stubbed out for potential specialization in the future.
"""
@inline function tuplemap(fn::F, tup1::Tuple, tup2::Tuple) where {F}
    map(fn, tup1, tup2)
end

"""
    rmap(fn, X...)

Recursively apply `fn` to each element of `X`
"""
rmap(fn::F, X) where {F} = fn(X)
rmap(fn::F, X, Y) where {F} = fn(X, Y)
rmap(fn::F, X::Tuple) where {F} = tuplemap(x -> rmap(fn, x), X)
rmap(fn::F, X::Tuple, Y::Tuple) where {F} =
    tuplemap((x, y) -> rmap(fn, x, y), X, Y)
rmap(fn::F, X::NamedTuple{names}) where {F, names} =
    NamedTuple{names}(rmap(fn, Tuple(X)))
rmap(fn::F, X::NamedTuple{names}, Y::NamedTuple{names}) where {F, names} =
    NamedTuple{names}(rmap(fn, Tuple(X), Tuple(Y)))

"""
    rmaptype(fn, T)

The return type of `rmap(fn, X::T)`.
"""
rmaptype(fn::F, ::Type{T}) where {F, T} = fn(T)
rmaptype(fn::F, ::Type{T}) where {F, T <: Tuple} =
    Tuple{tuplemap(fn, tuple(T.parameters...))...}
rmaptype(
    fn::F,
    ::Type{T},
) where {F, T <: NamedTuple{names, tup}} where {names, tup} =
    NamedTuple{names, rmaptype(fn, tup)}

"""
    rmul(X, Y)
    X ⊠ Y

Recursively scale each element of `X` by `Y`.
"""
rmul(X, Y) = rmap(*, X, Y)
rmul(w::Number, X) = rmap(x -> w * x, X)
rmul(X, w::Number) = rmap(x -> x * w, X)
rmul(w1::Number, w2::Number) = w1 * w2
const ⊠ = rmul

"""
    radd(X, Y)
    X ⊞ Y

Recursively add elements of `X` and `Y`.
"""
radd(X, Y) = rmap(+, X, Y)
radd(w::Number, X) = rmap(x -> w + x, X)
radd(X, w::Number) = rmap(x -> x + w, X)
radd(w1::Number, w2::Number) = w1 + w2
const ⊞ = radd

# Adapted from Base/operators.jl for general nary operator fallbacks
for op in (:rmul, :radd)
    @eval begin
        ($op)(a, b, c, xs...) = Base.afoldl($op, ($op)(($op)(a, b), c), xs...)
    end
end

"""
    rsub(X, Y)
    X ⊟ Y

Recursively subtract elements of `Y` from `X`.
"""
rsub(X) = rmap(-, X)
rsub(X, Y) = rmap(-, X, Y)
rsub(X, w::Number) = rmap(x -> x - w, X)
rsub(w::Number, X) = rmap(x -> w - x, X)
rsub(w1::Number, w2::Number) = w1 - w2
const ⊟ = rsub

"""
    rdiv(X, Y)

Recursively divide each element of `X` by `Y`
"""
rdiv(X, Y) = rmap(/, X, Y)
rdiv(X, w::Number) = rmap(x -> x / w, X)
rdiv(w::Number, X) = rmap(x -> w / x, X)
rdiv(w1::Number, w2::Number) = w1 / w2

"""
    rmuladd(w, X, Y)

Recursively add elements of `w * X + Y`.
"""
rmuladd(w::Number, X, Y) = rmap((x, y) -> muladd(w, x, y), X, Y)
rmuladd(X, w::Number, Y) = rmap((x, y) -> muladd(x, w, y), X, Y)
rmuladd(w::Number, x::Number, y::Number) = muladd(w, x, y)

"""
    rmatmul1(W, S, i, j)

Recursive matrix product along the 1st dimension of `S`. Equivalent to:

    mapreduce(⊠, ⊞, W[i,:], S[:,j])

"""
function rmatmul1(W, S, i, j)
    Nq = size(W, 2)
    r = W[i, 1] ⊠ S[1, j]
    for ii in 2:Nq
        r = rmuladd(W[i, ii], S[ii, j], r)
    end
    return r
end

"""
    rmatmul2(W, S, i, j)

Recursive matrix product along the 2nd dimension `S`. Equivalent to:

    mapreduce(⊠, ⊞, W[j,:], S[i, :])

"""
function rmatmul2(W, S, i, j)
    Nq = size(W, 2)
    r = W[j, 1] ⊠ S[i, 1]
    for jj in 2:Nq
        r = rmuladd(W[j, jj], S[i, jj], r)
    end
    return r
end

end # module


text:
using rendahl
using Test
using Parameters

function rbc_eq(Xl, X, Xf, ϵ, ϵ_sd, params)
    Cf, Rf, Kf, Yf, Zf = Xf
    C, R, K, Y, Z = X
    Cl, Rl, Kl, Yl, Zl = Xl

    @unpack β, α, γ, δ, ρ = params

    ϵ = ϵ[1] * ϵ_sd[1]

    # RBC Equations
    residual = [1.0 - β * Rf * Cf^(-γ) * C^(γ);
                R - α*Z*Kl^(α-1) - 1 + δ;
                K - (1-δ)*Kl - Y + C;
                Y - Z*Kl^α;
                log(Z) - ρ*log(Zl) - ϵ]

    return residual

end

rbc = @with_kw (
    β = 0.96,
    α = 0.33,
    γ = 2.0,
    δ = 0.1,
    ρ = 0.9
)

params = rbc()

function rbc_ss(params)

    #===========
    Get the steady state values
    for X
    ============#

    @unpack β, α, γ, δ, ρ = params

    Z = 1
    R = 1/β
    K = ((R - 1 + δ)/α)^(1/(α - 1))
    Y = K^α
    C = Y - δ*K

    ss = [C, R, K, Y, Z]

    return ss
end

func_ss = rendahl.get_ss(rbc_eq, ones(5), 0.8, [params])
closedform_rbc_ss = rbc_ss(params)

# Simple ODE
function ode_eqn(Xl, X, Xf, ϵ, ϵ_sd ,params)

    @unpack α_1, α_2, α_3 = params

    residual = [ α_1 * Xl[1] + α_2 * X[1] + α_3 * Xf[1] + ϵ[1] * ϵ_sd[1]]

    return residual

end

ode_params = @with_kw (
    α_1 = 0.75, α_2 = -2.0, α_3 = 1.0
)

simple_ode = ode_params()
ode_shocks = [0.0]

ode_sol = solve(ode_eqn, [simple_ode], ode_shocks,
            xinit = [1.0])

@testset "rendahl.jl" begin

    @test isapprox(func_ss, closedform_rbc_ss)
    @test isapprox(ode_sol.F[1], 0.5, atol = 1e-5)

end

ext: jl


text:
module IniPG

    using DataFrames
    using Random
    using StatsModels
    using ..DataMod, ..Utils

    export inipg

    mutable struct Particle
        position::Vector{Float64}
        velocity::Vector{Float64}
        pbest::Vector{Float64}
        pbest_fit::Float64
        pbest_nvars::Integer
    end

    mutable struct Swarm
        particles::Vector{Particle}
        gbest::Vector{Float64}
        gbest_fit::Float64
        gbest_nvars::Integer
        gbest_index::Integer
        Swarm() = new(Vector{Particle}(), Vector{Float64}(), Inf, typemax(Int), -1)
    end
    Base.getindex(S::Swarm, i::Int)::Particle = getindex(S.particles, i)
    Base.iterate(S::Swarm) = iterate(S.particles)
    Base.iterate(S::Swarm, state) = iterate(S.particles, state)
    Base.push!(S::Swarm, p::Particle) = push!(S.particles, p)
    Base.length(S::Swarm) = length(S.particles)

    count_selected_vars(position::Vector{Float64}, θ::Float64)::Integer =
        sum(position .> θ)
    count_selected_vars(part::Particle, θ::Float64)::Integer =
        count_selected_vars(part.position, θ)

    get_selected_vars(position::Vector{Float64}, xs::Vector{Symbol}, θ::Float64)::Vector{Symbol} =
        xs[position .> θ]
    get_selected_vars(part::Particle, xs::Vector{Symbol}, θ::Float64)::Vector{Symbol} =
        get_selected_vars(part.position, xs, θ)

    function fitness(position::Vector{Float64}, data::Data, θ::Float64)::Float64
        selected = get_selected_vars(position, data.xs, θ)
        models = get_models(data, selected)
        return prediction_error(models, data)
    end
    fitness(part::Particle, data::Data, θ::Float64)::Float64 =
        fitness(part.position, data, θ)

    function set_gbest!(swarm::Swarm, verbose::Bool)::Swarm
        minind = argmin([p.pbest_fit for p in swarm])
        bestpart = @inbounds swarm[minind]
        if bestpart.pbest_fit < swarm.gbest_fit || (bestpart.pbest_fit == swarm.gbest_fit && bestpart.pbest_nvars < swarm.gbest_nvars)
            verbose && println("gbest updated -- gbest is now particle #$minind")
            swarm.gbest = bestpart.pbest
            swarm.gbest_fit = bestpart.pbest_fit
            swarm.gbest_nvars = bestpart.pbest_nvars
            swarm.gbest_index = minind
        end
        return swarm
    end

    function initialize_swarm(n::Integer, data::Data, rng::AbstractRNG, θ::Float64, low_frac::Float64,
                              low_var_frac::Float64, high_var_frac_range::Tuple{Float64, Float64}, verbose::Bool)::Swarm
        swarm = Swarm()
        nlow = n * low_frac
        nvars = length(data.xs)
        nvarlow = round(Int, low_var_frac * nvars)
        @inbounds for i in 1:n
            nsel = i <= nlow ? nvarlow : rand(rng, round(Int, high_var_frac_range[1] * nvars):round(Int, high_var_frac_range[2] * nvars))
            vals = shuffle(rng, [(1 - θ) * rand(rng, nsel) .+ θ; θ * rand(rng, nvars - nsel)])
            pbest_fit = fitness(vals, data, θ)
            pbest_nvars = count_selected_vars(vals, θ)
            part = Particle(vals, zeros(nvars), vals, pbest_fit, pbest_nvars)
            push!(swarm, part)
        end
        set_gbest!(swarm, verbose)
        return swarm
    end

    function pso_loop!(swarm::Swarm, data::Data, rng::AbstractRNG, θ::Float64,
                       max_iter::Integer, inertia::Float64, c1::Float64, c2::Float64, verbose::Bool)::Swarm
        for i in 1:max_iter
            verbose && println("PSO iteration $i")
            r1 = rand(rng)
            r2 = rand(rng)
            p = 1
            for part in swarm
                part.velocity = inertia * part.velocity + c1 * r1 * (part.pbest - part.position) + c2 * r2 * (swarm.gbest - part.position)
                part.position += part.velocity
                score = fitness(part, data, θ)
                nvars = count_selected_vars(part, θ)
                if score < part.pbest_fit || (score == part.pbest_fit && nvars < part.pbest_nvars)
                    verbose && println("pbest updated #$p")
                    part.pbest = part.position
                    part.pbest_fit = score
                    part.pbest_nvars = nvars
                end
                p += 1
            end
            set_gbest!(swarm, verbose)
        end
        return swarm
    end

    function inipg(data::Data, rng::Union{AbstractRNG, Nothing} = nothing;
                   n::Integer = 30, θ::Float64 = 0.6, low_frac::Float64 = 2/3, low_var_frac::Float64 = 0.1,
                   high_var_frac_range::Tuple{Float64, Float64} = (0.5, 1.0), max_iter::Integer = 100,
                   inertia::Float64 = 0.7298, c1::Float64 = 1.49618, c2::Float64 = 1.49618, verbose::Bool = true)::Vector{RegressionModel}
        isnothing(rng) && (rng = MersenneTwister())
        swarm = initialize_swarm(n, data, rng, θ, low_frac, low_var_frac, high_var_frac_range, verbose)
        pso_loop!(swarm, data, rng, θ, max_iter, inertia, c1, c2, verbose)
        verbose && println("Best solution was found by particle #$(swarm.gbest_index)")
        return get_models(data, get_selected_vars(swarm.gbest, data.xs, θ))
    end

end

using .IniPG
export inipg

max_stars_repo_name: KasperNooteboom/thesis-rvfl-fs

max_issues_repo_name: KasperNooteboom/thesis-rvfl-fs

max_forks_repo_name: KasperNooteboom/thesis-rvfl-fs


text:
abstract type AbstractCommand end

struct Path
    entries::Vector{String}
end

Base.convert(::Type{Path}, x::AbstractVector{<:AbstractString}) = Path(x)
Base.vcat(p::Path, s...) = Path(vcat(p.entries, s...))
Base.show(io::IO, p::Path) = print(io, string('/', join(p.entries, '/')))

struct SetObject{O <: AbstractObject} <: AbstractCommand
    object::O
    path::Path
end

struct SetTransform{T <: Transformation} <: AbstractCommand
    tform::T
    path::Path
end

struct Delete <: AbstractCommand
    path::Path
end


text:
using Test
using ExecutableSpecifications.Gherkin: hastag, parsefeature, issuccessful

@testset "Tags                 " begin
    @testset "Feature tags" begin
        @testset "@tag1 is applied to a feature; The parsed feature has @tag1" begin
            text = """
            @tag1
            Feature: Some description
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test hastag(feature, "@tag1")
        end

        @testset "Feature without tags; The parsed feature does not have @tag1" begin
            text = """
            Feature: Some description
            """

            result = parsefeature(text)

            @test issuccessful(result)
            @test hastag(result.value, "@tag1") == false
        end

        @testset "Feature with multiple tags; The parsed feature has all tags" begin
            text = """
            @tag1 @tag2 @tag3
            Feature: Some description
            """

            result = parsefeature(text)

            @test issuccessful(result)
            @test hastag(result.value, "@tag1")
            @test hastag(result.value, "@tag2")
            @test hastag(result.value, "@tag3")
        end
    end

    @testset "Scenario tags" begin
        @testset "Scenario has one tag; The parsed scenario has tag1" begin
            text = """
            Feature: Some description

                @tag1
                Scenario: Some description
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test hastag(feature.scenarios[1], "@tag1")
        end

        @testset "Scenario has no tags; The parsed scenario does not have tag1" begin
            text = """
            Feature: Some description

                Scenario: Some description
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            @test hastag(result.value.scenarios[1], "@tag1") == false
        end

        @testset "Feature has tag1, but no the scenario; The parsed scenario does not have tag1" begin
            text = """
            @tag1
            Feature: Some description

                Scenario: Some description
                    Given a precondition
            """

            result = parsefeature(text)

            @test issuccessful(result)
            @test hastag(result.value.scenarios[1], "@tag1") == false
        end
    end

    @testset "Robustness" begin
        @testset "Tag @tag1-2 contains a hyphen; Tag is read as @tag1-2" begin
            text = """
            @tag1-2
            Feature: Some description
            """

            result = parsefeature(text)

            @test issuccessful(result)
            feature = result.value
            @test hastag(feature, "@tag1-2")
        end
    end
end


text:
module ModelPredictiveControl

import JuMP
import Ipopt
import OSQP
import ControlSystems
import LinearAlgebra
import ForwardDiff

export linear_mpc
export mpc_problem
export solve_mpc
export terminal_weight

include("problem_design.jl")
include("optim_design.jl")
include("solve.jl")

end


text:
using BinaryBuilder

name = "Ogg"
version = v"1.3.4"

sources = [
    ArchiveSource("https://downloads.xiph.org/releases/ogg/libogg-$(version).tar.xz",
                  "c163bc12bc300c401b6aa35907ac682671ea376f13ae0969a220f7ddf71893fe"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libogg-*/

# Don't let `configure` use `-ffast-math` ANYWHERE
sed -i.bak -e 's/-ffast-math//g' ./configure
./configure --prefix=$prefix --host=${target} --build=${MACHTYPE}
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
]

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

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
mutable struct LDA <: TopicModel
	"Latent Dirichlet allocation."

	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	alpha::Vector{Float64}
	beta::Matrix{Float64}
	beta_old::Matrix{Float64}
	beta_temp::Matrix{Float64}
	Elogtheta::VectorList{Float64}
	Elogtheta_old::VectorList{Float64}
	gamma::VectorList{Float64}
	phi::MatrixList{Float64}
	elbo::Float64

	function LDA(corp::Corpus, K::Integer)
		check_corp(corp)
		K > 0 || throw(ArgumentError("Number of topics must be a positive integer."))

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		alpha = ones(K)
		beta = rand(Dirichlet(V, 1.0), K)'
		beta_old = copy(beta)
		beta_temp = zeros(K, V)
		Elogtheta = [-eulergamma * ones(K) .- digamma(K) for _ in 1:M]
		Elogtheta_old = deepcopy(Elogtheta)
		gamma = [ones(K) for _ in 1:M]
		phi = [ones(K, N[d]) / K for d in 1:min(M, 1)]
		elbo = 0

		model = new(K, M, V, N, C, copy(corp), topics, alpha, beta, beta_old, beta_temp, Elogtheta, Elogtheta_old, gamma, phi, elbo)
		return model
	end
end

function Elogptheta(model::LDA, d::Int)
	"Compute E_q[log(P(theta))]."

	x = finite(loggamma(sum(model.alpha))) - finite(sum(loggamma.(model.alpha))) + dot(model.alpha .- 1, model.Elogtheta[d])
	return x
end

function Elogpz(model::LDA, d::Int)
	"Compute E_q[log(P(z))]."

	counts = model.corp[d].counts
	x = dot(model.phi[1] * counts, model.Elogtheta[d])
	return x
end

function Elogpw(model::LDA, d::Int)
	"Compute E_q[log(P(w))]."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi[1] .* log.(@boink model.beta[:,terms]) * counts)
	return x
end

function Elogqtheta(model::LDA, d::Int)
	"Compute E_q[log(q(theta))]."

	x = -entropy(Dirichlet(model.gamma[d]))
	return x
end

function Elogqz(model::LDA, d::Int)
	"Compute E_q[log(q(z))]."

	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[1][:,n])) for (n, c) in enumerate(counts)])
	return x
end

function update_elbo!(model::LDA)
	"Update the evidence lower bound."

	model.elbo = 0
	for d in 1:model.M
		terms = model.corp[d].terms
		@positive model.phi[1] = model.beta_old[:,terms] .* exp.(model.Elogtheta_old[d])
		model.phi[1] ./= sum(model.phi[1], dims=1)
		model.elbo += Elogptheta(model, d) + Elogpz(model, d) + Elogpw(model, d) - Elogqtheta(model, d) - Elogqz(model, d)
	end

	return model.elbo
end

function update_alpha!(model::LDA, niter::Integer, ntol::Real)
	"Update alpha."
	"Interior-point Newton's method with log-barrier and back-tracking line search."

	Elogtheta_sum = sum([model.Elogtheta[d] for d in 1:model.M])

	nu = model.K
	for _ in 1:niter
		rho = 1.0
		alpha_grad = nu ./ model.alpha + model.M * (digamma(sum(model.alpha)) .- digamma.(model.alpha)) + Elogtheta_sum
		alpha_invhess_diag = -1 ./ (model.M * trigamma.(model.alpha) + nu ./ model.alpha.^2)
		p = (alpha_grad .- dot(alpha_grad, alpha_invhess_diag) * (model.M * trigamma(sum(model.alpha)) + sum(alpha_invhess_diag))) .* alpha_invhess_diag

		while minimum(model.alpha - rho * p) < 0
			rho *= 0.5
		end
		@finite model.alpha -= rho * p

		if (rho * norm(alpha_grad) < ntol) & (nu / model.K < ntol)
			break
		end
		nu *= 0.5
	end
	@positive model.alpha
end

function update_beta!(model::LDA)
	"Reset beta variables."

	model.beta_old = model.beta
	model.beta = model.beta_temp ./ sum(model.beta_temp, dims=2)
	model.beta_temp = zeros(model.K, model.V)
end

function update_beta!(model::LDA, d::Int)
	"Update beta."
	"Analytic."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.beta_temp[:,terms] += model.phi[1] .* counts'
end

function update_Elogtheta!(model::LDA, d::Int)
	"Update E[log(theta)]."
	"Analytic."

	model.Elogtheta_old[d] = model.Elogtheta[d]
	model.Elogtheta[d] = digamma.(model.gamma[d]) .- digamma(sum(model.gamma[d]))
end

function update_gamma!(model::LDA, d::Int)
	"Update gamma."
	"Analytic."

	counts = model.corp[d].counts
	@positive model.gamma[d] = model.alpha + model.phi[1] * counts
end

function update_phi!(model::LDA, d::Int)
	"Update phi."
	"Analytic."

	terms = model.corp[d].terms
	@positive model.phi[1] = model.beta[:,terms] .* exp.(model.Elogtheta[d])
	model.phi[1] ./= sum(model.phi[1], dims=1)
end

function train!(model::LDA; iter::Integer=150, tol::Real=1.0, niter::Integer=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, checkelbo::Real=1, printelbo::Bool=true)
	"Coordinate ascent optimization procedure for latent Dirichlet allocation variational Bayes algorithm."

	check_model(model)
	all([tol, ntol, vtol] .>= 0)										|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, niter, viter] .>= 0)										|| throw(ArgumentError("Iteration parameters must be nonnegative."))
	(isa(checkelbo, Integer) & (checkelbo > 0)) | (checkelbo == Inf)	|| throw(ArgumentError("checkelbo parameter must be a positive integer or Inf."))
	all([isempty(doc) for doc in model.corp]) && (iter = 0)
	(checkelbo <= iter) && update_elbo!(model)

	for k in 1:iter
		for d in 1:model.M
			for v in 1:viter
				update_phi!(model, d)
				update_gamma!(model, d)
				update_Elogtheta!(model, d)
				if norm(model.Elogtheta[d] - model.Elogtheta_old[d]) < vtol
					break
				end
			end
			update_beta!(model, d)
		end
		update_beta!(model)
		update_alpha!(model, niter, ntol)

		if check_elbo!(model, checkelbo, printelbo, k, tol)
			break
		end
	end

	model.topics = [reverse(sortperm(vec(model.beta[i,:]))) for i in 1:model.K]
	nothing
end


text:
using StaticArrays

abstract type AbstrUTriMatrix end

"""
    UTriMatrix(dim, x)

A upper-triangular matrix of dimension `dim X dim`, whose entries have
the same type as `x`. The parameter `x` is used only for its type. Its
value is not used. A matrix is constructed with a call such as

    a = UTriMatrix(4, BigInt(1))

which generates a 4X4 upper triangular matrix whose entries are of
type BigInt. Then, a[4,3] returns the entry at row 4, column 3, and
a[4,3] throws a BoundsError. The entries of the matrix are initialized
with zeroes.

For values such that `isbits(x)` is true, use BitsUTriMatrix, which is
faster.

Entries are stored in column-major order.

    x x x x  offset of col 1 = 0 = 1*0/2
      x x x  offset of col 2 = 1 = 2*1/2
        x x  offset of col 3 = 3 = 3*2/2
          x  offset of col 4 = 6 = 4*3/2
             offset of col i = (i*(i-1))/2

    A[i, j] = data[i + offset[j]]
    A[4, 4] = data[4 + 6] = data[10]
"""
mutable struct UTriMatrix{N,M,T} <: AbstrUTriMatrix
    data :: Vector{T}        # elements in matrix
    offset :: SVector{N,Int} # num elements before column j

    function UTriMatrix(dim, x)
        t = typeof(x) # element type
        n = Int(dim)
        m = (n*(n+1)) ÷ 2
        mat = new{n, m, t}()

        mat.data = zeros(t, m)
        offset = MVector{n,Int}(undef)

        for i in 1:n
            offset[i] = (i*(i-1)) ÷ 2
        end
        mat.offset = SVector(offset)
        mat
    end
end

"""
    BitsUTriMatrix(dim, x)

Like UTriMatrix, but stores data in a static array for speed. The
elements of the matrix must be a bits type.  The entries of the matrix
are not initialized.

An uppeer-triangular matrix of dimension `dim X dim`, whose entries have
the same type as `x`. The parameter `x` is used only for its type. Its
value is not used. The type of x must be a bits type, so
that `isbits(x)` returns `true`.

A matrix is constructed with a call such as

    a = UTriMatrix(4, 1.0)

which generates a 4X4 upper triangular matrix whose entries are of
type Float64. Then, a[3,4] returns the entry at row 3, column 4, and
a[4,3] throws a BoundsError.
"""
mutable struct BitsUTriMatrix{N,M,T} <: AbstrUTriMatrix
    data :: MVector{M,T}     # elements in matrix
    offset :: SVector{N,Int} # num elements before column j

    function BitsUTriMatrix(dim, x)
        t = typeof(x) # element type
        !isbitstype(t) && throw(DomainError(t,
            "BitsUTriMatrix can only hold bitstype values."))
        n = Int(dim)
        m = (n*(n+1)) ÷ 2
        mat = new{n, m, t}()

        mat.data = MVector{m, t}(undef)
        offset = MVector{n,Int}(undef)

        for i in 1:n
            offset[i] = (i*(i-1)) ÷ 2
        end
        mat.offset = SVector(offset)
        mat
    end
end

@inline function Base.size(A::AbstrUTriMatrix) where {N,M,T}
    dim = length(A.offset)
    (dim, dim)
end

@inline Base.length(A::AbstrUTriMatrix) where {N,M,T} = length(A.data)

@inline Base.IndexStyle(::AbstrUTriMatrix) where {N,M,T} = IndexCartesian()

@inline function Base.checkbounds(A::AbstrUTriMatrix, i::Int,
                          j::Int) where {N,M,T}
    i > j && throw(BoundsError(A, (i,j)))
end

@inline function Base.getindex(A::AbstrUTriMatrix, i::Int,
                               j::Int) where {N,M,T}
    @boundscheck checkbounds(A, i, j)
    A.data[i + A.offset[j]]
end

# Fails if value cannot be converted to type T.
@inline function Base.setindex!(A::AbstrUTriMatrix, value, i::Int,
                                j::Int) where {N,M,T}
    @boundscheck checkbounds(A, i, j)
    A.data[i + A.offset[j]] = value
end

"""
Write a representation of an UTriMatrix to io.
"""
function Base.show(io::IO, A::AbstrUTriMatrix) where {N,M,T}
    n = length(A.offset)
    for i in 1:n
        for j in 1:(i-1)
            print(io, j>1 ? " " : "", zero(typeof(A[1,1])))
        end
        for j in i:n
            print(io, j>1 ? " " : "", A[i,j])
        end
        println(io)
    end
end


text:
module Transaction

using ..Lib

export transaction, TransactionAborted, TransactionEnded

struct TransactionAborted <: Exception end
struct TransactionEnded <: Exception end

# This is taken from the C++ implementation:
#
# https://github.com/pmem/libpmemobj-cpp/blob/master/include/libpmemobj%2B%2B/transaction.hpp
#
# Pass a closure to the transaction for a given pool pointer. Should automatically deal with
# cleaning up and rolling back
function transaction(f, pool)
    # Start transaction
    Lib.tx_begin(pool)

    local obj
    try
        obj = f()
    catch err
        # Error was thrown while executing.
        # Need to call `_tx_abort` and then throw the exception up the stack
        if Lib.tx_stage() == Lib.TX_STAGE_WORK
            Lib.tx_abort(-1)
        end
        Lib.tx_end()
        throw(err)
    end

    # Handle cleanup
    stage = Lib.tx_stage()

    # Optimistic case - everything went well, just need to commit and end
    if stage == Lib.TX_STAGE_WORK
        Lib.tx_commit()

    # Abort was called in "f"
    elseif stage == Lib.TX_STAGE_ONABORT
        Lib.tx_end()
        throw(TransactionAborted())

    # "tx_end" was called in "f"
    elseif stage == Lib.TX_STAGE_NONE
        throw(TransactionEnded())
    end

    Lib.tx_end()
    return obj
end

end # module


text:
### Boundary layer routines

export u, v, vorticity, dudy, d2udy2, eta, η, d99, dstar, theta, Cf, falknerskan

struct BLSolution{ETAT,UT}
  β :: Float64
  Vw :: Float64
  ηmax :: Float64
  η :: ETAT
  f :: UT
  fp :: UT
  fpp :: UT
  h0 :: Float64
  d99coeff :: Float64
  dstarcoeff :: Float64
  θcoeff :: Float64
  cfcoeff :: Float64
end

"""
    u(sol::BLSolution) -> Vector{Float64}

Return the streamwise velocity array, scaled by the external velocity, ``u/U_e``
"""
u(sol::MAE250B.BLSolution) = sol.fp;

"""
    η(sol::BLSolution) -> Vector{Float64}
    eta(sol::BLSolution) -> Vector{Float64}

Return the array of wall-normal coordinates, scaled by the
boundary layer thickness factor ``\\eta = y/\\delta(x) = y/(\\nu x/U_e)^{1/2}``
"""
η(sol::MAE250B.BLSolution) = sol.η
eta(sol::MAE250B.BLSolution) = η(sol)

"""
    v(sol::BLSolution) -> Vector{Float64}

Return the array of wall-normal velocities, ``v``, scaled by the
boundary layer factor: ``v Re_x^{1/2}/U_e``
"""
v(sol::MAE250B.BLSolution) = 0.5*(sol.η.*sol.fp - sol.f)

"""
    vorticity(sol::BLSolution) -> Vector{Float64}

Return the array of vorticity, ``\\omega``, scaled by the
boundary layer factor: ``\\omega \\delta(x)/U_e``
"""
vorticity(sol::MAE250B.BLSolution) = -sol.fpp

"""
    dudy(sol::BLSolution) -> Vector{Float64}

Return the array of ``du/dy``, scaled by the
boundary layer factor: ``du/dy (\\delta(x)/U_e)``
"""
dudy(sol::MAE250B.BLSolution) = sol.fpp

"""
    d2udy2(sol::BLSolution) -> Vector{Float64}

Return the array of ``d^2u/dy^2``, scaled by the
boundary layer factor: ``d^2u/dy^2 (\\delta^2(x)/U_e)``
"""
function d2udy2(sol::MAE250B.BLSolution)
  uyy = zero(sol.fpp)
  uyy[1] = (sol.fpp[2] - sol.fpp[1])/(sol.η[2] - sol.η[1])
  for i in 2:lastindex(sol.fpp)-1
    uyy[i] = (sol.fpp[i+1] - sol.fpp[i-1])/(sol.η[i+1] - sol.η[i-1])
  end
  uyy[end] = (sol.fpp[end] - sol.fpp[end-1])/(sol.η[end] - sol.η[end-1])
  return uyy
end

"""
    d99(sol::BLSolution) -> Vector{Float64}

Return the boundary layer 99 percent thickness divided by
the boundary layer thickness factor ``\\delta_{99}(x)/\\delta(x)``
"""
d99(sol::MAE250B.BLSolution) = sol.d99coeff

"""
    dstar(sol::BLSolution) -> Vector{Float64}

Return the boundary layer displacement thickness divided by
the boundary layer thickness factor ``\\delta^{*}(x)/\\delta(x)``
"""
dstar(sol::MAE250B.BLSolution) = sol.dstarcoeff

"""
    theta(sol::BLSolution) -> Vector{Float64}

Return the boundary layer displacement thickness divided by
the boundary layer thickness factor ``\\theta(x)/\\delta(x)``
"""
theta(sol::MAE250B.BLSolution) = sol.θcoeff

"""
    Cf(sol::BLSolution) -> Vector{Float64}

Return the boundary layer skin friction coefficient scaled by
the boundary layer factor ``C_f(x)Re_x^{1/2}``
"""
Cf(sol::MAE250B.BLSolution) = sol.cfcoeff

### Solver routines ###

function fsrhs!(du,u,p,t)
  _, m, _ = p

  du[1] = u[2]
  du[2] = u[3]
  du[3] = -0.5*(m+1)*u[1]*u[3] - m*(1-u[2]^2)
end

function fs_integrate(h0,p)
  Vw, m, ηmax = p

  f0 = -2Vw
  g0 = 0.0
  gL = 1.0

  u0 = [f0;g0;h0]
  ηspan = (0.0,ηmax)

  prob = ODEProblem(fsrhs!,u0,ηspan,p)
  sol = OrdinaryDiffEq.solve(prob,Tsit5(),reltol=1e-12,abstol=1e-15)

  resid = sol[2,end] - gL

  return resid, sol

end

"""
    falknerskan(β[,Vw=0][,ηmax=10][,h0init=1.3])

Compute the Falkner-Skan boundary layer solution for parameter `β`,
where `β` can be a value larger than -0.1999 (the separation case).
It returns the data in a solution structure which can be used
to return individual results, e.g.,

`sol = falknerskan(0)`

for ``\\beta=0`` and then `u(sol)` to return the streamwise velocity profile.

You can use the optional `Vw=` argument to set a wall-normal velocity
at the wall. (This value represents ``V_w Re_x^{1/2}/U_e``.)
"""
function falknerskan(β;Vw = 0.0,ηmax = 10.0,h0init=1.3,kwargs...)

  m = β/(2-β)
  p = [Vw,m,ηmax]

  h0 = find_zero(x -> fs_integrate(x,p)[1],h0init;kwargs...)
  resid,sol = fs_integrate(h0,p)

  sol = BLSolution(convert(Float64,β),convert(Float64,Vw),convert(Float64,ηmax),_fs_eta(sol),_fs_streamfunction(sol),_fs_velocity(sol),
                    _fs_vorticity(sol),
                    convert(Float64,h0),blthickness_99(sol), blthickness_displacement(sol),
                    blthickness_momentum(sol), blskinfriction(sol))

  return sol

end

_fs_streamfunction(sol::OrdinaryDiffEq.ODESolution) = sol[1,:]
_fs_velocity(sol::OrdinaryDiffEq.ODESolution) = sol[2,:]
_fs_vorticity(sol::OrdinaryDiffEq.ODESolution) = sol[3,:]
_fs_eta(sol::OrdinaryDiffEq.ODESolution) = sol.t

function blthickness_99(sol::OrdinaryDiffEq.ODESolution)

    u, η = _fs_velocity(sol), _fs_eta(sol)

    sign_diff = sign.(u .- 0.99)  # +1 where u > 0.99, -1 otherwise
    i0 = findfirst(diff(sign_diff) .== 2.0)
    return sol.t[i0]+ (0.99-u[i0])/(u[i0+1]-u[i0])*(η[i0+1]-η[i0])
end

function blthickness_displacement(sol::OrdinaryDiffEq.ODESolution)
    u, η = _fs_velocity(sol), _fs_eta(sol)

    uhalf = 0.5*(u[2:end].+u[1:end-1])
    return sum((1 .- uhalf).*diff(η))
end

function blthickness_momentum(sol::OrdinaryDiffEq.ODESolution)
    u, η = _fs_velocity(sol), _fs_eta(sol)

    uhalf = 0.5*(u[2:end].+u[1:end-1])
    return sum(uhalf.*(1 .- uhalf).*diff(η))
end

blskinfriction(sol::OrdinaryDiffEq.ODESolution) = 2*sol[3,1]


text:
# BraidingTensor:
# special (2,2) tensor that implements a standard braiding operation
#====================================================================#
"""
    struct BraidingTensor{S<:IndexSpace} <: AbstractTensorMap{S, 2, 2}

Specific subtype of [`AbstractTensorMap`](@ref) for representing the braiding tensor that
braids the first input over the second input; its inverse can be obtained as the adjoint.
"""
struct BraidingTensor{S<:IndexSpace, A} <: AbstractTensorMap{S, 2, 2}
    V1::S
    V2::S
    adjoint::Bool
    function BraidingTensor(V1::S, V2::S, adjoint::Bool = false, ::Type{A} = Matrix{ComplexF64}) where
                {S<:IndexSpace, A<:DenseMatrix}
        for a in sectors(V1)
            for b in sectors(V2)
                for c in (a ⊗ b)
                    Nsymbol(a, b, c) == Nsymbol(b, a, c) ||
                        throw(ArgumentError("Cannot define a braiding between $a and $b"))
                end
            end
        end
        return new{S,A}(V1, V2, adjoint)
        # partial construction: only construct rowr and colr when needed
    end
end

Base.adjoint(b::BraidingTensor{S,A}) where {S<:IndexSpace, A<:DenseMatrix} =
    BraidingTensor(b.V1, b.V2, !b.adjoint, A)

domain(b::BraidingTensor) = b.adjoint ? b.V2 ⊗ b.V1 : b.V1 ⊗ b.V2
codomain(b::BraidingTensor) = b.adjoint ? b.V1 ⊗ b.V2 : b.V2 ⊗ b.V1

storagetype(::Type{BraidingTensor{S,A}}) where {S<:IndexSpace, A<:DenseMatrix} = A

blocksectors(b::BraidingTensor) = blocksectors(b.V1 ⊗ b.V2)
hasblock(b::BraidingTensor, s::Sector) = s ∈ blocksectors(b)

function fusiontrees(b::BraidingTensor)
    codom = codomain(b)
    dom = domain(b)
    I = sectortype(b)
    F = fusiontreetype(I, 2)
    rowr = SectorDict{I, FusionTreeDict{F, UnitRange{Int}}}()
    colr = SectorDict{I, FusionTreeDict{F, UnitRange{Int}}}()
    for c in blocksectors(codom)
        rowrc = FusionTreeDict{F, UnitRange{Int}}()
        colrc = FusionTreeDict{F, UnitRange{Int}}()
        offset1 = 0
        for s1 in sectors(codom)
            for f1 in fusiontrees(s1, c, map(isdual, codom.spaces))
                r = (offset1 + 1):(offset1 + dim(codom, s1))
                push!(rowrc, f1 => r)
                offset1 = last(r)
            end
        end
        dim1 = offset1
        offset2 = 0
        for s2 in sectors(dom)
            for f2 in fusiontrees(s2, c, map(isdual, dom.spaces))
                r = (offset2 + 1):(offset2 + dim(dom, s2))
                push!(colrc, f2 => r)
                offset2 = last(r)
            end
        end
        dim2 = offset2
        push!(rowr, c=>rowrc)
        push!(colr, c=>colrc)
    end
    return TensorKeyIterator(rowr, colr)
end

function Base.getindex(b::BraidingTensor{S}) where S
    sectortype(S) == Trivial || throw(SectorMismatch())
    (V1, V2) = domain(b)
    d = (dim(V2), dim(V1), dim(V1), dim(V2))
    return sreshape(StridedView(block(b, Trivial())), d)
end

@inline function Base.getindex(b::BraidingTensor, f1::FusionTree{I,2}, f2::FusionTree{I,2}) where {I<:Sector}
    I == sectortype(b) || throw(SectorMismatch())
    c = f1.coupled
    V1, V2 = domain(b)
    @boundscheck begin
        c == f2.coupled || throw(SectorMismatch())
        ((f1.uncoupled[1] ∈ sectors(V2)) && (f2.uncoupled[1] ∈ sectors(V1))) || throw(SectorMismatch())
        ((f1.uncoupled[2] ∈ sectors(V1)) && (f2.uncoupled[2] ∈ sectors(V2))) || throw(SectorMismatch())
    end
    @inbounds begin
        d = (dims(V2 ⊗ V1, f1.uncoupled)..., dims(V1 ⊗ V2, f2.uncoupled)...)
        n1 = d[1]*d[2]
        n2 = d[3]*d[4]
        data = fill!(storagetype(b)(undef, (n1, n2)), zero(eltype(b)))
        a1, a2 = f2.uncoupled
        if f1.uncoupled == (a2, a1)
            braiddict = artin_braid(f2, 1; inv = b.adjoint)
            r = get(braiddict, f1, zero(valtype(braiddict)))
            si = 1 + d[1]*d[2]*d[3]
            sj = d[1] + d[1]*d[2]
            @inbounds for i = 1:d[1], j = 1:d[2]
                data[(i-1)*si + (j-1)*sj + 1] = r
            end
        end
        return sreshape(StridedView(data), d)
    end
end

Base.copy(b::BraidingTensor) = copy!(similar(b), b)
function Base.copy!(t::TensorMap, b::BraidingTensor)
    space(t) == space(b) || throw(SectorMismatch())
    fill!(t, zero(eltype(t)))
    for (f1, f2) in fusiontrees(t)
        data = t[f1, f2]
        if sectortype(t) == Trivial
            r = one(eltype(t))
        else
            a1, a2 = f2.uncoupled
            c = f2.coupled
            f1.uncoupled == (a2, a1) || continue
            braiddict = artin_braid(f2, 1; inv = b.adjoint)
            r = convert(eltype(t), get(braiddict, f1, zero(valtype(braiddict))))
        end
        for i = 1:size(data, 1), j = 1:size(data, 2)
            data[i, j, j, i] = r
        end
    end
    return t
end
TensorMap(b::BraidingTensor) = copy(b)
Base.convert(::Type{TensorMap}, b::BraidingTensor) = copy(b)

function block(b::BraidingTensor, s::Sector)
    sectortype(b) == typeof(s) || throw(SectorMismatch())
    (V1, V2) = domain(b)
    if sectortype(b) == Trivial
        d1, d2 = dim(V1), dim(V2)
        n = d1*d2
        data = fill!(storagetype(b)(undef, (n, n)), zero(eltype(b)))
        si = 1 + d2*d1*d1
        sj = d2 + d2*d1
        @inbounds for i = 1:d2, j = 1:d1
            data[(i-1)*si + (j-1)*sj + 1] = one(eltype(b))
        end
        return data
    end
    n = blockdim(domain(b), s)
    data = fill!(storagetype(b)(undef, (n, n)), zero(eltype(b)))
    iter = fusiontrees(b) # actually contains information about ranges as well
    for (f2, r2) in iter.colr[s]
        for (f1, r1) in iter.rowr[s]
            a1, a2 = f2.uncoupled
            d1 = dim(V1, a1)
            d2 = dim(V2, a2)
            f1.uncoupled == (a2, a1) || continue
            braiddict = artin_braid(f2, 1; inv = b.adjoint)
            r = convert(eltype(b), get(braiddict, f1, zero(valtype(braiddict))))
            si = 1 + n*d1
            sj = d2 + n
            start = first(r1) + (first(r2)-1) * n
            @inbounds for i = 1:d2, j=1:d1
                data[(i-1)*si + (j-1)*sj + start] = r
            end
        end
    end
    return data
end

blocks(b::BraidingTensor) = blocks(TensorMap(b))

function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S},
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple{2}, cindA::IndexTuple{2},
                            oindB::IndexTuple, cindB::IndexTuple{2},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])'

    if BraidingStyle(sectortype(B)) isa Bosonic
        return add!(α, B, β, C, reverse(cindB), oindB)
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindA[1]] > braidingtensor_levels[cindA[2]]
    for (f1, f2) in fusiontrees(B)
        local newtrees
        for ((f1′, f2′), coeff′) in transpose(f1, f2, cindB, oindB)
            for (f1′′, coeff′′) in artin_braid(f1′, 1, inv = inv_braid)

                f12 = (f1′′, f2′)
                coeff = coeff′*coeff′′
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        for ((f1′,f2′), coeff) in newtrees
            TO._add!(coeff*α, B[f1, f2], true, C[f1′,f2′], (reverse(cindB)...,oindB...))
        end
    end
    return C
end
function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor,
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple, cindA::IndexTuple{2},
                            oindB::IndexTuple{2}, cindB::IndexTuple{2},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}
    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])'

    if BraidingStyle(sectortype(A)) isa Bosonic
        return add!(α, A, β, C, oindA, reverse(cindA))
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindB[1]] > braidingtensor_levels[cindB[2]]
    for (f1, f2) in fusiontrees(A)
        local newtrees
        for ((f1′,f2′), coeff′) in transpose(f1, f2, oindA, cindA)
            for (f2′′, coeff′′) in artin_braid(f2′, 1, inv = inv_braid)
                f12 = (f1′,f2′′)
                coeff = coeff′*conj(coeff′′)
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        for ((f1′,f2′), coeff) in newtrees
            TO._add!(coeff*α, A[f1, f2], true, C[f1′,f2′], (oindA...,reverse(cindA)...))
        end
    end
    C
end

function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S},
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple{0}, cindA::IndexTuple{4},
                            oindB::IndexTuple, cindB::IndexTuple{4},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])' &&
                space(B, cindB[3]) == space(A, cindA[3])' &&
                space(B, cindB[4]) == space(A, cindA[4])'

    if BraidingStyle(sectortype(B)) isa Bosonic
        return trace!(α, B, β, C, (), oindB, (cindB[1], cindB[2]), (cindB[3], cindB[4]))
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    I = sectortype(B)
    u = one(I)
    f₀ = FusionTree{I}((), u, (), (), ())
    braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]
    for (f1, f2) in fusiontrees(B)
        local newtrees
        for ((f1′,f2′), coeff′) in transpose(f1, f2, cindB, oindB)
            f1′.coupled == u || continue
            a = f1′.uncoupled[1]
            b = f1′.uncoupled[2]
            f1′.uncoupled[3] == dual(a) || continue
            f1′.uncoupled[4] == dual(b) || continue
            # should be automatic by matching spaces:
            # f1′.isdual[1] != f1′.isdual[3] || continue
            # f1′.isdual[2] != f1′.isdual[4] || continue
            for (f1′′, coeff′′) in artin_braid(f1′, 2, inv = inv_braid)
                f1′′.innerlines[1] == u || continue
                coeff = coeff′ * coeff′′ * sqrtdim(a) * sqrtdim(b)
                if f1′′.isdual[1]
                    coeff *= frobeniusschur(a)
                end
                if f1′′.isdual[3]
                    coeff *= frobeniusschur(b)
                end
                f12 = (f₀, f2′)
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        @isdefined(newtrees) || continue
        for ((f1′, f2′), coeff) in newtrees
            TO._trace!(coeff*α, B[f1, f2], true, C[f1′,f2′], oindB,
                                            (cindB[1], cindB[2]), (cindB[3], cindB[4]))
        end
    end
    return C
end

function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor,
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple, cindA::IndexTuple{4},
                            oindB::IndexTuple{0}, cindB::IndexTuple{4},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])' &&
                space(B, cindB[3]) == space(A, cindA[3])' &&
                space(B, cindB[4]) == space(A, cindA[4])'

    if BraidingStyle(sectortype(B)) isa Bosonic
        return trace!(α, A, β, C, oindA, (), (cindA[1], cindA[2]), (cindA[3], cindA[4]))
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    I = sectortype(B)
    u = one(I)
    f₀ = FusionTree{I}((), u, (), (), ())
    braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindB[2]] > braidingtensor_levels[cindB[3]]
    for (f1, f2) in fusiontrees(A)
        local newtrees
        for ((f1′, f2′), coeff′) in transpose(f1, f2, oindA, cindA)
            f2′.coupled == u || continue
            a = f2′.uncoupled[1]
            b = f2′.uncoupled[2]
            f2′.uncoupled[3] == dual(a) || continue
            f2′.uncoupled[4] == dual(b) || continue
            # should be automatic by matching spaces:
            # f2′.isdual[1] != f2′.isdual[3] || continue
            # f2′.isdual[3] != f2′.isdual[4] || continue
            for (f2′′, coeff′′) in artin_braid(f2′, 2, inv = inv_braid)
                f2′′.innerlines[1] == u || continue
                coeff = coeff′ * conj(coeff′′ * sqrtdim(a) * sqrtdim(b))
                if f2′′.isdual[1]
                    coeff *= conj(frobeniusschur(a))
                end
                if f2′′.isdual[3]
                    coeff *= conj(frobeniusschur(b))
                end
                f12 = (f1′, f₀)
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        @isdefined(newtrees) || continue
        for ((f1′, f2′), coeff) in newtrees
            TO._trace!(coeff*α, A[f1, f2], true, C[f1′,f2′], oindA,
                                            (cindA[1], cindA[2]), (cindA[3], cindA[4]))
        end
    end
    return C
end

function planar_contract!(α, A::BraidingTensor, B::AbstractTensorMap{S},
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple{1}, cindA::IndexTuple{3},
                            oindB::IndexTuple, cindB::IndexTuple{3},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])' &&
                space(B, cindB[3]) == space(A, cindA[3])'

    if BraidingStyle(sectortype(B)) isa Bosonic
        return trace!(α, B, β, C, (cindB[2],), oindB, (cindB[1],), (cindB[3],))
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    I = sectortype(B)
    u = one(I)
    braidingtensor_levels = A.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindA[2]] > braidingtensor_levels[cindA[3]]
    for (f1, f2) in fusiontrees(B)
        local newtrees
        for ((f1′,f2′), coeff′) in transpose(f1, f2, cindB, oindB)
            a = f1′.uncoupled[1]
            b = f1′.uncoupled[2]
            b == f1′.coupled || continue
            a == dual(f1′.uncoupled[3]) || continue
            # should be automatic by matching spaces:
            # f1′.isdual[1] != f1.isdual[3] || continue
            for (f1′′, coeff′′) in artin_braid(f1′, 2, inv = inv_braid)
                f1′′.innerlines[1] == u || continue
                coeff = coeff′ * coeff′′ * sqrtdim(a)
                if f1′′.isdual[1]
                    coeff *= frobeniusschur(a)
                end
                f1′′′ = FusionTree{I}((b,), b, (f1′′.isdual[3],), (), ())
                f12 = (f1′′′, f2′)
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        @isdefined(newtrees) || continue
        for ((f1′,f2′), coeff) in newtrees
            TO._trace!(coeff*α, B[f1, f2], true, C[f1′,f2′],
                                        (cindB[2], oindB...) , (cindB[1],), (cindB[3],))
        end
    end
    return C
end

function planar_contract!(α, A::AbstractTensorMap{S}, B::BraidingTensor,
                            β, C::AbstractTensorMap{S},
                            oindA::IndexTuple, cindA::IndexTuple{3},
                            oindB::IndexTuple{1}, cindB::IndexTuple{3},
                            p1::IndexTuple, p2::IndexTuple,
                            syms::Union{Nothing, NTuple{3, Symbol}}) where {S}

    codA, domA = codomainind(A), domainind(A)
    codB, domB = codomainind(B), domainind(B)
    oindA, cindA, oindB, cindB =
        reorder_indices(codA, domA, codB, domB, oindA, cindA, oindB, cindB, p1, p2)

    @assert space(B, cindB[1]) == space(A, cindA[1])' &&
                space(B, cindB[2]) == space(A, cindA[2])' &&
                space(B, cindB[3]) == space(A, cindA[3])'

    if BraidingStyle(sectortype(A)) isa Bosonic
        return trace!(α, A, β, C, oindA, (cindA[2],), (cindA[1],), (cindA[3],))
    end

    if iszero(β)
        fill!(C, β)
    elseif β != 1
        rmul!(C, β)
    end
    I = sectortype(B)
    u = one(I)
    braidingtensor_levels = B.adjoint ? (1,2,2,1) : (2,1,1,2)
    inv_braid = braidingtensor_levels[cindB[2]] > braidingtensor_levels[cindB[3]]
    for (f1, f2) in fusiontrees(A)
        local newtrees
        for ((f1′,f2′), coeff′) in transpose(f1, f2, oindA, cindA)
            a = f2′.uncoupled[1]
            b = f2′.uncoupled[2]
            b == f2′.coupled || continue
            a == dual(f2′.uncoupled[3]) || continue
            # should be automatic by matching spaces:
            # f2′.isdual[1] != f2.isdual[3] || continue
            for (f2′′, coeff′′) in artin_braid(f2′, 2, inv = inv_braid)
                f2′′.innerlines[1] == u || continue
                coeff = coeff′ * conj(coeff′′ * sqrtdim(a))
                if f2′′.isdual[1]
                    coeff *= conj(frobeniusschur(a))
                end
                f2′′′ = FusionTree{I}((b,), b, (f2′′.isdual[3],), (), ())
                f12 = (f1′, f2′′′)
                if @isdefined newtrees
                    newtrees[f12] = get(newtrees, f12, zero(coeff)) + coeff
                else
                    newtrees = Dict(f12 => coeff)
                end
            end
        end
        @isdefined(newtrees) || continue
        for ((f1′,f2′), coeff) in newtrees
            TO._trace!(coeff*α, A[f1, f2], true, C[f1′,f2′],
                                            (oindA...,cindA[2]) , (cindA[1],), (cindA[3],))
        end
    end
    return C
end

has_shared_permute(t::BraidingTensor, args...) = false
function cached_permute(sym::Symbol, t::BraidingTensor, p1, p2; copy=false)
    tp = TO.cached_similar_from_indices(sym, eltype(t), p1, p2, t, :N)
    return add!(true, t, false, tp, p1, p2)
end


text:
using Documenter
using AutoPages
using AutoPages: gather_pages, replace_reverse, transform_page

tutorials_dir = joinpath(@__DIR__, "..", "tutorials") # tutorials directory
mkpath(tutorials_dir)
tutorials, tutorials_list = gather_pages(;
  directory=tutorials_dir,
  extension_filter=x->endswith(x, ".jl"),
  transform_extension=x->replace_reverse(x, ".jl" => ".md"; count=1),
  remove_first_level=true)

# #### Let's suppose we've used Literate to generate markdown:
# for tutorial in tutorials_list
#     mkpath(joinpath(@__DIR__,"src",dirname(tutorial)))
# end
# for tutorial in tutorials_list
#     title = transform_page(basename(tutorial))
#     open(joinpath(@__DIR__,"src",tutorial), "w") do io
#     print(io, "# $(title)")
#     end
# end
# ####

format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true",
        collapselevel = 1,
        )

makedocs(
    sitename = "AutoPages",
    format = format,
    clean = true,
    modules = [Documenter, AutoPages],
    pages = Any["Home" => "index.md",
                "Tutorials" => tutorials],
)

deploydocs(
    repo = "github.com/charleskawczynski/AutoPages.jl.git",
    target = "build",
)


text:
using SEAL

include("1_bfv_basics.jl")
include("2_encoders.jl")
include("3_levels.jl")
include("4_ckks_basics.jl")
include("5_rotation.jl")
include("6_serialization.jl")
include("7_performance.jl")

function seal_examples()
  println("Microsoft SEAL version: ", version())
  while true
    println("+---------------------------------------------------------+")
    println("| The following examples should be executed while reading |")
    println("| comments in associated files in examples/.              |")
    println("+---------------------------------------------------------+")
    println("| Examples                   | Source Files               |")
    println("+----------------------------+----------------------------+")
    println("| 1. BFV Basics              | 1_bfv_basics.jl            |")
    println("| 2. Encoders                | 2_encoders.jl              |")
    println("| 3. Levels                  | 3_levels.jl                |")
    println("| 4. CKKS Basics             | 4_ckks_basics.jl           |")
    println("| 5. Rotation                | 5_rotation.jl              |")
    println("| 6. Serialization           | 6_serialization.jl         |")
    println("| 7. Performance Test        | 7_performance.jl           |")
    println("+----------------------------+----------------------------+")

    megabytes = alloc_byte_count(memory_manager_get_pool()) >> 20
    @printf("[% 7d MB] Total allocation from the memory pool\n", megabytes)

    selection = 0
    invalid = true
    while true
      println()
      print("> Run example (1 ~ 7) or exit (0): ")
      input = readline()
      if !isopen(stdin)
        return
      end
      try
        selection = parse(Int, input)
        break
      catch
        println("  [Beep~~] Invalid option: type 0 ~ 7")
        continue
      end
    end

    if selection == 1
      example_bfv_basics()
    elseif selection == 2
      example_encoders()
    elseif selection == 3
      example_levels()
    elseif selection == 4
      example_ckks_basics()
    elseif selection == 5
      example_rotation()
    elseif selection == 6
      example_serialization()
    elseif selection == 7
      example_performance_test()
    elseif selection == 0
      return
    end
  end

  return
end


text:
### Combination methods for p-values ###

"""
    combine(PValues, PValueCombination)

Combine p-values

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, FisherCombination())
0.007616871850449092
julia> combine(pvals, StoufferCombination())
0.007098326181265917
```

```jldoctest
julia> subtypes(PValueCombination)
7-element Array{Union{DataType, UnionAll},1}:
 MultipleTesting.FisherCombination
 MultipleTesting.LogitCombination
 MultipleTesting.MinimumCombination
 MultipleTesting.SimesCombination
 MultipleTesting.StoufferCombination
 MultipleTesting.TippettCombination
 MultipleTesting.WilkinsonCombination
```

# See also

`PValueCombination`s:

[`FisherCombination`](@ref) [`LogitCombination`](@ref) [`StoufferCombination`](@ref) [`TippettCombination`](@ref) [`SimesCombination`](@ref) [`WilkinsonCombination`](@ref) [`MinimumCombination`](@ref)

"""
function combine end

function combine(pValues::AbstractVector{T}, method::M)::T where {T<:AbstractFloat, M<:PValueCombination}
    combine(PValues(pValues), method)
end

## Fisher combination ##

"""
Fisher's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, FisherCombination())
0.007616871850449092
```

# References

Fisher, R.A. (1925).
Statistical methods for research workers
(Genesis Publishing Pvt Ltd).

"""
struct FisherCombination <: PValueCombination
end

function combine(pValues::PValues{T}, method::FisherCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0
        return NaN
    end
    x = -2 * sum(log.(pValues))
    p = ccdf(Chisq(2n), x)
    return p
end

## Logit combination ##

"""
Logit p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, LogitCombination())
0.006434494635148462
```

# References

Mudholkar, G.S., and George, E.O. (1977).
The Logit Statistic for Combining Probabilities - An Overview
(Rochester University NY, Dept of Statistics).

"""
struct LogitCombination <: PValueCombination
end

function combine(pValues::PValues{T}, method::LogitCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    c = sqrt( (5n+2)*n*pi^2 / ((5n+4)*3) )
    x = -sum(log.(pValues./(1 .- pValues))) / c
    p = ccdf(TDist(5n+4), x)
    return p
end

## Stouffer combination ##

"""
Stouffer's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, StoufferCombination())
0.007098326181265917

julia> weights = [1.0, 2.0, 0.4, 1.5];

julia> combine(pvals, weights, StoufferCombination())
0.007331653763696742
```

# References

Stouffer, S.A. (1949).
The American soldier. Vol. 1: Adjustment during army life
(Princeton University Press).

Liptak, T. (1958).
On the combination of independent tests.
Magyar Tud Akad Mat Kutato Int Kozl 3, 171–197.

"""
struct StoufferCombination <: PValueCombination
end

function combine(pValues::PValues{T}, method::StoufferCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    z = cquantile.(Ref(Normal()), pValues)
    z = sum(z) ./ sqrt(n)
    p = ccdf(Normal(), z)
    return p
end

function combine(pValues::PValues{T}, weights::AbstractVector{R}, method::StoufferCombination)::T where {T<:AbstractFloat, R<:Real}
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    if minimum(pValues) == 0 || maximum(pValues) == 1
        return NaN
    end
    z = cquantile.(Ref(Normal()), pValues) .* weights
    z = sum(z) ./ sqrt(sum(abs2, weights))
    p = ccdf(Normal(), z)
    return p
end

function combine(pValues::AbstractVector{T}, weights::AbstractVector{R}, method::StoufferCombination)::T where {T<:AbstractFloat, R<:Real}
    combine(PValues(pValues), weights, method)
end

## Tippett combination ##

"""
Tippett's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, TippettCombination())
0.039403990000000055
```

# References

Tippett, L.H.C. (1931). The Methods of Statistics. An introduction mainly for
workers in the biological sciences.

"""
struct TippettCombination <: PValueCombination
end

function combine(pValues::PValues{T}, method::TippettCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    p = 1 - (1 - minimum(pValues))^n
    return p
end

## Simes combination ##

"""
Simes's p-value combination

# Examples

```jldoctest
julia> pvals = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pvals, SimesCombination())
0.04
```

# References

Simes, R.J. (1986). An improved Bonferroni procedure for multiple tests of
significance. Biometrika 73, 751–754.

"""
struct SimesCombination <: PValueCombination
end

function combine(pValues::PValues{T}, method::SimesCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    pValues = sort(pValues)
    p = n * minimum(pValues./(1:n))
    return p
end

## Wilkinson combination ##

"""
Wilkinson's p-value combination

# Examples

```jldoctest
julia> pv = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pv, WilkinsonCombination(1))  # combination with rank 1
0.03940399000000003

julia> combine(pv, WilkinsonCombination(4))  # combination with rank 4
0.0625
```

# References

Wilkinson, B. (1951). A statistical consideration in psychological research.
Psychological Bulletin 48, 156.

"""
struct WilkinsonCombination <: PValueCombination
    rank::Integer

    function WilkinsonCombination(rank)
        if rank < 1
            throw(ArgumentError("Rank must be positive."))
        end
        return new(rank)
    end
end

function combine(pValues::PValues{T}, method::WilkinsonCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    rank = method.rank
    if rank < 1 || rank > n
        throw(ArgumentError("Rank must be in 1,..,$(n)"))
    end
    p_rank = partialsort(pValues, rank)
    p = cdf(Beta(rank, n-rank+1), p_rank)
    return p
end

## Generalised minimum combination ##

"""
Minimum of adjusted p-value combination

# Examples

```jldoctest
julia> pv = PValues([0.01, 0.02, 0.3, 0.5]);

julia> combine(pv, MinimumCombination(BenjaminiHochberg()))
0.04

julia> combine(pv, MinimumCombination(ForwardStop()))
0.01005033585350145
```

"""
struct MinimumCombination <: PValueCombination
    adjustment::PValueAdjustment
end

function combine(pValues::PValues{T}, method::MinimumCombination)::T where T<:AbstractFloat
    n = length(pValues)
    if n == 1
        return pValues[1]
    end
    padj = adjust(pValues, method.adjustment)
    p = minimum(padj)
    return p
end


text:
# # Loading, converting, and saving models

#md # [![](https://mybinder.org/badge_logo.svg)](@__BINDER_ROOT_URL__/notebooks/@__NAME__.ipynb)
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/notebooks/@__NAME__.ipynb)

# `COBREXA` can load models stored in `.mat`, `.json`, and `.xml` formats (with
# the latter denoting SBML formatted models).

# We will primarily use the *E. Coli* "core" model to demonstrate the utilities
# found in `COBREXA`. First, let's download the model in several formats.

## Downloads the model files if they don't already exist
!isfile("e_coli_core.mat") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.mat", "e_coli_core.mat");
!isfile("e_coli_core.json") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.json", "e_coli_core.json");
!isfile("e_coli_core.xml") &&
    download("http://bigg.ucsd.edu/static/models/e_coli_core.xml", "e_coli_core.xml");

# Now, load the package:

using COBREXA

#md # !!! tip "Save bandwidth!"
#md #     The published models usually do not change very often. It is
#md #     therefore pretty useful to save them to a central location and load
#md #     them from there. That saves your time, and does not unnecessarily
#md #     consume the connectivity resources of the model repository.

# ## Loading models

# Load the models using the [`load_model`](@ref) function. Each model is able to
# "pretty-print" itself, hiding the inner complexity.

mat_model = load_model("e_coli_core.mat")
#

json_model = load_model("e_coli_core.json")
#

sbml_model = load_model("e_coli_core.xml")
#

#md # !!! note "Note: `load_model` infers the output type from the file extension"
#md #       Notice how each model was read into memory as a model type corresponding
#md #       to its file type, i.e. the file ending with `.json` loaded as a
#md #       [`JSONModel`](@ref), the file ending with `.mat` loaded as [`MATModel`](@ref), and the
#md #       file ending with `.xml` loaded as an [`SBMLModel`](@ref).

# You can directly inspect the model objects, although only with a specific way
# for each specific type.

# JSON models contain their corresponding JSON:

json_model.json

# SBML models contain a complicated structure from [`SBML.jl`
# package](https://github.com/LCSB-BioCore/SBML.jl):

typeof(sbml_model.sbml)

# MAT models contain MATLAB data:

mat_model.mat

# ## Using the generic interface to access model details

# To prevent the complexities of object representation, `COBREXA.jl` uses a set
# of generic interface functions that extract various important information
# from all supported model types. This approach ensures that the analysis
# functions can work on any data.

# For example, you can check the reactions and metabolites contained in SBML
# and JSON models using the same accessor:

reactions(json_model)

#

reactions(sbml_model)

#

issetequal(reactions(json_model), reactions(mat_model)) # do models contain the same reactions?

# All accessors are defined in a single file in COBREXA source code; you may
# therefore get a list of all accessors as follows:
using InteractiveUtils

for method in filter(
    x -> endswith(string(x.file), "MetabolicModel.jl"),
    InteractiveUtils.methodswith(MetabolicModel, COBREXA),
)
    println(method.name)
end

# ## Converting between model types

# It is possible to convert model types to-and-fro. To do this, use the
# `convert` function, which is overloaded from Julia's `Base`.

#md # !!! danger "Data loss may occur when converting between models"
#md #     The conversion of models only uses the data accessible through the
#md #     generic accessors. Other data may get lost.

m = convert(MATModel, json_model)

# `m` will now contain the MATLAB-style matrix representation of the model:

Matrix(m.mat["S"])

# The loading and conversion can be combined using a shortcut:

m = load_model(MATModel, "e_coli_core.json")

# ## Saving and exporting models

# `COBREXA.jl` supports exporting the models in JSON and MAT format, using [`save_model`](@ref).

save_model(m, "converted_model.json")
save_model(m, "converted_model.mat")

# If you need a non-standard suffix, use the type-specific saving functions:

save_json_model(m, "file.without.a.good.suffix")
save_mat_model(m, "another.file.matlab")

# If you are saving the models only for future processing in Julia environment,
# it is often wasteful to encode the models to external formats and decode them
# back. Instead, you can use the "native" Julia data format, accessible with
# package `Serialization`.
#
# This way, you can use `serialize` to save even the [`StandardModel`](@ref)
# that has no file format associated:

using Serialization

sm = convert(StandardModel, m)

open(f -> serialize(f, sm), "myModel.stdmodel", "w")

# The models can then be loaded back using `deserialize`:

sm2 = deserialize("myModel.stdmodel")
issetequal(metabolites(sm), metabolites(sm2))

# This form of loading operation is usually pretty quick:
t = @elapsed deserialize("myModel.stdmodel")
@info "Deserialization took $t seconds"
# Notably, large and complicated models with thousands of reactions and
# annotations can take seconds to decode properly. Serialization allows you to
# almost completely remove this overhead, and scales well to tens of millions
# of reactions.


text:
include("BinDicts.jl"); using .BinDicts;
include("GSDicts.jl"); using .GSDicts;
include("S3Dicts.jl"); using .S3Dicts;


text:
abstract type Proposal{P} end

struct StaticProposal{issymmetric,P} <: Proposal{P}
    proposal::P
end
const SymmetricStaticProposal{P} = StaticProposal{true,P}

StaticProposal(proposal) = StaticProposal{false}(proposal)
function StaticProposal{issymmetric}(proposal) where {issymmetric}
    return StaticProposal{issymmetric,typeof(proposal)}(proposal)
end

struct RandomWalkProposal{issymmetric,P} <: Proposal{P}
    proposal::P
end
const SymmetricRandomWalkProposal{P} = RandomWalkProposal{true,P}

RandomWalkProposal(proposal) = RandomWalkProposal{false}(proposal)
function RandomWalkProposal{issymmetric}(proposal) where {issymmetric}
    return RandomWalkProposal{issymmetric,typeof(proposal)}(proposal)
end

# Random draws
Base.rand(p::Proposal, args...) = rand(Random.GLOBAL_RNG, p, args...)
Base.rand(rng::Random.AbstractRNG, p::Proposal{<:Distribution}) = rand(rng, p.proposal)
function Base.rand(rng::Random.AbstractRNG, p::Proposal{<:AbstractArray})
    return map(x -> rand(rng, x), p.proposal)
end

# Densities
Distributions.logpdf(p::Proposal{<:Distribution}, v) = logpdf(p.proposal, v)
function Distributions.logpdf(p::Proposal{<:AbstractArray}, v)
    # `mapreduce` with multiple iterators requires Julia 1.2 or later
    return mapreduce(((pi, vi),) -> logpdf(pi, vi), +, zip(p.proposal, v))
end

###############
# Random Walk #
###############

function propose(
    rng::Random.AbstractRNG,
    proposal::RandomWalkProposal{issymmetric,<:Union{Distribution,AbstractArray}},
    ::DensityModel
) where {issymmetric}
    return rand(rng, proposal)
end

function propose(
    rng::Random.AbstractRNG,
    proposal::RandomWalkProposal{issymmetric,<:Union{Distribution,AbstractArray}},
    model::DensityModel,
    t
) where {issymmetric}
    return t + rand(rng, proposal)
end

function q(
    proposal::RandomWalkProposal{issymmetric,<:Union{Distribution,AbstractArray}},
    t,
    t_cond
) where {issymmetric}
    return logpdf(proposal, t - t_cond)
end

##########
# Static #
##########

function propose(
    rng::Random.AbstractRNG,
    proposal::StaticProposal{issymmetric,<:Union{Distribution,AbstractArray}},
    model::DensityModel,
    t=nothing
) where {issymmetric}
    return rand(rng, proposal)
end

function q(
    proposal::StaticProposal{issymmetric,<:Union{Distribution,AbstractArray}},
    t,
    t_cond
) where {issymmetric}
    return logpdf(proposal, t)
end

############
# Function #
############

# function definition with abstract types requires Julia 1.3 or later
for T in (:StaticProposal, :RandomWalkProposal)
    @eval begin
        function (p::$T{issymmetric,<:Function})() where {issymmetric}
            return $T{issymmetric}(p.proposal())
        end
        function (p::$T{issymmetric,<:Function})(t) where {issymmetric}
            return $T{issymmetric}(p.proposal(t))
        end
    end
end

function propose(
    rng::Random.AbstractRNG,
    proposal::Proposal{<:Function},
    model::DensityModel
)
    return propose(rng, proposal(), model)
end

function propose(
    rng::Random.AbstractRNG,
    proposal::Proposal{<:Function},
    model::DensityModel,
    t
)
    return propose(rng, proposal(t), model)
end

function q(
    proposal::Proposal{<:Function},
    t,
    t_cond
)
    return q(proposal(t_cond), t, t_cond)
end

####################
# Multiple proposals
####################

function propose(
    rng::Random.AbstractRNG,
    proposals::AbstractArray{<:Proposal},
    model::DensityModel,
)
    return map(proposals) do proposal
        return propose(rng, proposal, model)
    end
end
function propose(
    rng::Random.AbstractRNG,
    proposals::AbstractArray{<:Proposal},
    model::DensityModel,
    ts,
)
    return map(proposals, ts) do proposal, t
        return propose(rng, proposal, model, t)
    end
end

@generated function propose(
    rng::Random.AbstractRNG,
    proposals::NamedTuple{names},
    model::DensityModel,
) where {names}
    isempty(names) && return :(NamedTuple())
    expr = Expr(:tuple)
    expr.args = Any[:($name = propose(rng, proposals.$name, model)) for name in names]
    return expr
end

@generated function propose(
    rng::Random.AbstractRNG,
    proposals::NamedTuple{names},
    model::DensityModel,
    ts,
) where {names}
    isempty(names) && return :(NamedTuple())
    expr = Expr(:tuple)
    expr.args = Any[
        :($name = propose(rng, proposals.$name, model, ts.$name)) for name in names
    ]
    return expr
end

"""
    logratio_proposal_density(proposal, state, candidate)

Compute the log-ratio of the proposal densities in the Metropolis-Hastings algorithm.

The log-ratio of the proposal densities is defined as
```math
\\log \\frac{g(x | x')}{g(x' | x)},
```
where ``x`` is the current state, ``x'`` is the proposed candidate for the next state,
and ``g(y' | y)`` is the conditional probability of proposing state ``y'`` given state
``y`` (proposal density).
"""
function logratio_proposal_density(proposal::Proposal, state, candidate)
    return q(proposal, state, candidate) - q(proposal, candidate, state)
end

# ratio is always 0 for symmetric proposals
logratio_proposal_density(::RandomWalkProposal{true}, state, candidate) = 0
logratio_proposal_density(::StaticProposal{true}, state, candidate) = 0

# type stable implementation for `NamedTuple`s
function logratio_proposal_density(
    proposals::NamedTuple{names}, states::NamedTuple, candidates::NamedTuple
) where {names}
    if @generated
        args = map(names) do name
            :(logratio_proposal_density(
                proposals[$(QuoteNode(name))],
                states[$(QuoteNode(name))],
                candidates[$(QuoteNode(name))],
            ))
        end
        return :(+($(args...)))
    else
        return sum(names) do name
            return logratio_proposal_density(
                proposals[name], states[name], candidates[name]
            )
        end
    end
end

# use recursion for `Tuple`s to ensure type stability
logratio_proposal_density(proposals::Tuple{}, states::Tuple, candidates::Tuple) = 0
function logratio_proposal_density(
    proposals::Tuple{<:Proposal}, states::Tuple, candidates::Tuple
)
    return logratio_proposal_density(first(proposals), first(states), first(candidates))
end
function logratio_proposal_density(proposals::Tuple, states::Tuple, candidates::Tuple)
    valfirst = logratio_proposal_density(first(proposals), first(states), first(candidates))
    valtail = logratio_proposal_density(
        Base.tail(proposals), Base.tail(states), Base.tail(candidates)
    )
    return valfirst + valtail
end

# fallback for general iterators (arrays etc.) - possibly not type stable!
function logratio_proposal_density(proposals, states, candidates)
    return sum(zip(proposals, states, candidates)) do (proposal, state, candidate)
        return logratio_proposal_density(proposal, state, candidate)
    end
end


text:
using JSON
using Pkg
using Test

# Setup GHA logger in CI
using Logging: global_logger
using GitHubActions: GitHubActionsLogger
get(ENV, "GITHUB_ACTIONS", "false") == "true" && global_logger(GitHubActionsLogger())

for exercise_type in ("concept", "practice")
    for exercise in readdir(joinpath("exercises", exercise_type))
        # Allow only testing specified exercises
        if !isempty(ARGS) && !(exercise in ARGS)
            continue
        end

        exercise_path = joinpath("exercises", exercise_type, exercise)
        isdir(exercise_path) || continue

        # Skip exercise tests if the Julia version doesn't meet the required version as specified in .meta/config.json
        cfg = JSON.parsefile(joinpath(exercise_path, ".meta", "config.json"))
        required_version_spec = Pkg.Types.semver_spec(get(cfg, "language_versions", "1.0"))
        if VERSION ∉ required_version_spec
            @info "$exercise requires Julia $required_version_spec, skipping tests."
            println()
            continue
        end

        # Concept exercises have an exemplar, practice exercises merely an example for CI purposes
        example_or_exemplar = exercise_type == "concept" ? "exemplar" : "example"

        # Create an anonymous module so that exercises are tested in separate scopes
        m = Module()

        # When testing the example solution, all tests must pass, even those that are marked as skipped or broken.
        # The student will not be affected by this.
        # Overwrite @test_skip and @test_broken with @test
        @eval m using Test
        @eval m $(Symbol("@test_skip")) = $(Symbol("@test"))
        @eval m $(Symbol("@test_broken")) = $(Symbol("@test"))

        # runtests.jl includes the solution by calling `include("slug.jl")`
        # Our anonymous module doesn't have `include(s::String)` defined,
        # so we define our own.
        @eval m include(s) = Base.include($m, joinpath(".meta", $example_or_exemplar * ".jl"))
        @info "[$(uppercase(exercise_type))] Testing $exercise"
        Base.include(m, joinpath(exercise_path, "runtests.jl"))

        println() # to make the output more readable
    end
end


text:
using Test, WeakRefStrings

# convenience constructor for testing
function str(x::Union{Missing, String}, e=nothing)
    data = Vector{UInt8}(x)
    poslen = PosLen(1, sizeof(coalesce(x, "")), x === missing, e !== nothing && UInt8(e) in data)
    return PosLenString(data, poslen, UInt8(something(e, 0x00)))
end

# convenience constructor for testing
function strs(x::Vector, e=nothing)
    len = sum(x -> sizeof(coalesce(x, "")), x)
    data = Vector{UInt8}(undef, len)
    poslens = Vector{PosLen}(undef, length(x))
    pos = 1
    anymissing = false
    for (i, s) in enumerate(x)
        if s !== missing
            slen = sizeof(s)
            bytes = codeunits(s)
            copyto!(data, pos, bytes, 1, slen)
            poslens[i] = PosLen(pos, slen, false, e !== nothing && e in bytes)
            pos += slen
        else
            anymissing = true
            poslens[i] = PosLen(pos, 0, true, false)
        end
    end
    return PosLenStringVector{anymissing ? Union{Missing, PosLenString} : PosLenString}(data, poslens, something(e, 0x00))
end

@testset "PosLenStrings" begin

# PosLen
poslen = PosLen(0, 0, false, false)
@test Base.bitcast(UInt64, poslen) == UInt64(0)

poslen = PosLen(0, 0, true, false)
@test poslen.missingvalue

poslen = PosLen(0, 0, false, true)
@test poslen.escapedvalue

poslen = PosLen(1, 3, false, false)
@test poslen.pos == 1
@test poslen.len == 3

@test_throws ArgumentError PosLen(1, 1048576, false, false)

# PosLenString
for y in ("", "hey", "🍕", "\\\\hey\\\\hey")
    x = str(y, '\\')
    y = unescape_string(y)
    @test codeunits(x) == codeunits(y)
    @test sizeof(x) == sizeof(y)
    @test ncodeunits(x) == ncodeunits(y)
    @test length(x) == length(y)
    @test codeunit(x) == UInt8
    @test lastindex(x) == lastindex(y)
    @test isempty(x) == isempty(y)
    @test String(x) == y
    @test Symbol(x) == Symbol(y)
    @test Vector{UInt8}(x) == Vector{UInt8}(y)
    @test Array{UInt8}(x) == Array{UInt8}(y)
    @test isascii(x) == isascii(y)
    @test x * x == y * y
    @test x^5 == y^5
    @test string(x) == string(y)
    @test join([x, x]) == join([y, y])
    @test reverse(x) == reverse(y)
    y != "" && @test startswith(x, x[1]) == startswith(y, x[1])
    y != "" && @test endswith(x, x[1]) == endswith(y, x[1])
    y != "" && @test findfirst(==(x[1]), x) === findfirst(==(x[1]), y)
    y != "" && @test findlast(==(x[1]), x) === findlast(==(x[1]), y)
    y != "" && @test findnext(==(x[1]), x, 1) === findnext(==(x[1]), y, 1)
    y != "" && @test findprev(==(x[1]), x, length(x)) === findprev(==(x[1]), y, length(x))
    @test lpad(x, 12) == lpad(y, 12)
    @test rpad(x, 12) == rpad(y, 12)
    y != "" && @test replace(x, x[1] => 'a') == replace(y, x[1] => 'a')
    r1 = match(Regex(x), x)
    r2 = match(Regex(y), y)
    @test r1 === r2 || r1.match == r2.match
    for i = 1:ncodeunits(x)
        @test codeunit(x, i) == codeunit(y, i)
        @test isvalid(x, i) == isvalid(y, i)
        @test thisind(x, i) == thisind(y, i)
        @test nextind(x, i) == nextind(y, i)
        @test prevind(x, i) == prevind(y, i)
        @test iterate(x, i) == iterate(y, i)
    end
    for i = 1:length(x)
        if isvalid(x, i)
            @test x[i] == y[i]
        end
    end
    @test x == x
    @test x == y
    @test y == x
    @test hash(x) == hash(y)
    @test cmp(x, x) == 0
    @test cmp(x, y) == 0
    @test cmp(y, x) == 0
end

# PosLenStringVector
x = strs(["hey", "there", "sailor", "esc\"aped"], UInt8('\\'))
@test x == ["hey", "there", "sailor", "esc\"aped"]
@test x[1] == "hey"

x = strs(["hey", "there", "sailor", "esc\"aped", missing], UInt8('\\'))
@test isequal(x, ["hey", "there", "sailor", "esc\"aped", missing])

x = strs(["hey", "there", "sailor", "esc\"aped"], UInt8('"'))
@test x == ["hey", "there", "sailor", "escaped"]
@test isassigned(x, 1)
@test x[1:3] == ["hey", "there", "sailor"]

end


text:
"""
    op = JopImag(dom)

Extract the imaginary part of a complex input array where `dom::JetSpace{<:Complex}`.
"""
JopImag(dom::JetSpace{Complex{T}}) where {T} = JopLn(dom = dom, rng = JetSpace(T,size(dom)), df! = JopImag_df!, df′! = JopImag_df′!)
export JopImag

JopImag_df!(d::AbstractArray{T}, m::AbstractArray{Complex{T}}; kwargs...) where {T} = d .= imag.(m)
JopImag_df′!(m::AbstractArray{Complex{T}}, d::AbstractArray{T}; kwargs...) where {T} = m .= im .* d


text:
"""
Google Cloud Container Engine API
"""
module _container

export container

using ..api
using ...root

"""
Google Cloud Container Engine API root.
"""
container = APIRoot(
    "https://container.googleapis.com/v1/projects/{project}/zones/{zone}",
    Dict(
        "cloud-platform" => "Full access to all resources and services in the specified Cloud Platform project.",
    );
    ServerConfig=APIResource("";
        getServerconfig=APIMethod(:GET, "serverconfig", "Returns configuration info about the Container Engine service."),
    ),
    Cluster=APIResource("clusters";
        create=APIMethod(:POST, "", "Creates a cluster, consisting of the specified number and type of Google Compute Engine instances."),
        delete=APIMethod(:DELETE, "{cluster}", "Deletes the cluster, including the Kubernetes endpoint and all worker nodes."),
        get=APIMethod(:GET, "{cluster}", "Gets the details of a specific cluster."),
        list=APIMethod(:GET, "", "Lists all clusters owned by a project in either the specified zone or all zones."),
        update=APIMethod(:PUT, "{cluster}", "Updates the settings of a specific cluster."),
    ),
    NodePool=APIResource("clusters/{cluster}/nodePools";
        create=APIMethod(:POST, "", "Creates a node pool for a cluster."),
        delete=APIMethod(:DELETE, "{nodePool}", "Deletes a node pool from a cluster."),
        get=APIMethod(:GET, "{nodePool}", "Retrieves the node pool requested."),
        list=APIMethod(:GET, "", "Lists the node pools for a cluster."),
    ),
    Operation=APIResource("operations";
        get=APIMethod(:GET, "{operation}", "Gets the specified operation."),
        list=APIMethod(:GET, "", "Lists all operations in a project in a specific zone or all zones."),
    ),
)

end


text:
# This file was generated, do not modify it. # hide
fprs, tprs, thresholds = roc(ŷ, y[test])

figure(figsize=(8,6))
plot(fprs, tprs)

xlabel("False Positive Rate", fontsize=14)
ylabel("True Positive Rate", fontsize=14)
xticks(fontsize=12)
yticks(fontsize=12)

savefig(joinpath(@OUTPUT, "ISL-lab-4-roc.svg")) # hide


text:
"""
```
solvect(m::AbstractCTModel; reduction_settings_check::Bool = false)
```

Driver to compute the model solution and augment transition matrices.

### Inputs

- `m`: the model object

### Outputs
 - TTT, RRR, and CCC matrices of the state transition equation:
    ```
    S_t = TTT*S_{t-1} + RRR*ϵ_t + CCC
    ```
"""
function solve(m::AbstractCTModel; sparse_mat::Bool = true, check_Γ0::Bool = true,
               complex_decomp::Bool = false, alt_policy::Bool = false)

    # Get equilibrium condition matrices
    Γ0, Γ1, Ψ, Π, C = eqcond(m)

    # Check that Γ0 is approximately the identity, otherwise clean it by solving static conditions
    if check_Γ0
        try
            n_total = get_setting(m, :n_jump_vars) + get_setting(m, :n_state_vars)
            @assert Matrix{Float64}(I, n_total, n_total) ≈ Γ0[1:n_total, 1:n_total]
        catch
            Γ0, Γ1, Ψ, Π, C, basis_redundant, inv_basis_redundant = solve_static_conditions(Γ0, Γ1, Ψ, Π, C)
        end
    end
    if !(@isdefined inv_basis_redundant)
        inv_basis_redundant = SparseMatrixCSC{Float64}(I, size(Γ0,1), size(Γ0,1))
        basis_redundant = SparseMatrixCSC{Float64}(I, size(Γ0,1), size(Γ0,1))
    end

    if sparse_mat
        Γ0 = sparse(Γ0); Γ1 = sparse(Γ1); Ψ = sparse(Ψ); Π = sparse(Π)
        C = sparse(reshape(C, length(C), 1))
    end

    # krylov reduction
    if get_setting(m, :reduce_state_vars)
        Γ0, Γ1, Ψ, Π, C, basis_kry, inv_basis_kry = krylov_reduction(m, Γ0, Γ1, Ψ, Π, C)
    end

    # value function reduction via spline projection
    if get_setting(m, :reduce_v)
        Γ0, Γ1, Ψ, Π, C, basis_spl, inv_basis_spl = valuef_reduction(m, Γ0, Γ1, Ψ, Π, C)
    end

    # Compute inverse basis for Z matrix and IRFs transformation
    if get_setting(m, :reduce_v) && get_setting(m, :reduce_state_vars)
        inverse_basis = inv_basis_redundant * inv_basis_kry * inv_basis_spl
        basis = basis_spl * basis_kry * basis_redundant
    elseif get_setting(m, :reduce_state_vars)
        inverse_basis = inv_basis_redundant * inv_basis_kry
        basis = basis_kry * basis_redundant
    else
        inverse_basis = inv_basis_redundant
        basis = basis_redundant
    end

    # Solve LRE model
    TTT_gensys, CCC_gensys, RRR_gensys, ~, ~, eu = gensysct!(Matrix{Float64}(Γ1), Matrix{Float64}(C),
                                                             Matrix{Float64}(Ψ), Matrix{Float64}(Π),
                                                             complex_decomposition = complex_decomp)

    # Check for LAPACK exception, existence and uniqueness
    if eu[1] != 1 || eu[2] != 1
        throw(GensysError())
    end

    TTT_gensys = real(TTT_gensys)
    RRR_gensys = real(RRR_gensys)
    CCC_gensys = vec(real(CCC_gensys))

    # If you wanted to augment states, you would do it here

    return TTT_gensys, RRR_gensys, CCC_gensys, inverse_basis, basis
end


text:
#! /usr/bin/env julia

using Arcus
using Base.Test

import Base.isapprox
isapprox(a::Float64, b::Float64) = isapprox(a, b, atol = 2eps(typeof(a))) || isapprox(1/a, 1/b, atol = 2eps())

@testset begin

@testset "creation" begin

  @test Arc(0.0) != nothing
  @test Arc(pi/2) != nothing
  @test Arcd(44) != nothing
  @test Arc(1.0, 1.0) != nothing
  @test_throws MethodError Arc(1.0, 2.0, 3.0)

end

@testset "dedicated values" begin
  @test Arc(0.0) == zero(Arc)
  @test Float64(Arcd(90)) == pi / 2
  @test Float64(Arcd(180)) == 1pi
  @test Float64(Arcd(270)) == -pi / 2
  @test Arcd(360) == Arcd(0)
end

# @testset "arithmetic" begin
  @test -Arc(pi*3/8) == Arc(-pi*3/8)
  @test -Arc(pi*3/4) == Arc(-pi*3/4)
  @test -Arc(pi*5/4) == Arc(-pi*5/4)
  @test -Arc(pi*7/4) == Arc(-pi*7/4)
  @test -Arc(pi/3) == Arc(5pi/3)
  @test Arc(1) + Arc(1) == 2Arc(1)
  @test Arc(1) + Arc(2) ≈ Arc(3)
  @test Arc(3) - Arc(2) ≈ Arc(1)
  @test -Arc(1) + Arc(3) ≈ Arc(2)
  @test Arcd(190) / 170 ≈ Arcd(-1.0)
  @test Arc(0.3) * 30 ≈ Arc(mod(9.0, 2pi))
  @test 20Arcd(260) == Arcd(260) * 20
# end

@testset "representation as sine 0" begin
  @test reinterpret(Float64, Arc(0)) == 0
end

@testset "representation as sine $(w)" for w in linspace(-pi/4, pi/4, 11)
  @test reinterpret(Float64, Arc(w)) ≈ sin(w)
end

const logsp = logspace(-16, 306, 5) * (nextfloat(0.0)*1e16)

@testset "representation as sine $(w)" for w in logsp
  @test reinterpret(Float64, Arc(w)) == sin(w)
end

const ep = eps(2.0)
@testset "trigonometric functions $w" for w in 0:45:360
  @test sin(Arcd(w)) ≈ sind(w)
  @test cos(Arcd(w)) ≈ cosd(w)
  @test tan(Arcd(w)) ≈ tand(w)
  @test csc(Arcd(w)) ≈ cscd(w)
  @test sec(Arcd(w)) ≈ secd(w)
  @test cot(Arcd(w)) ≈ cotd(w)
end

@testset "trigonometric functions $w" for w in 2/7:2/3:2pi
  @test sin(Arc(w)) ≈ sin(w)
  @test cos(Arc(w)) ≈ cos(w)
  @test tan(Arc(w)) ≈ tan(w)
  @test csc(Arc(w)) ≈ csc(w)
  @test sec(Arc(w)) ≈ sec(w)
  @test cot(Arc(w)) ≈ cot(w)
end

end


text:
using LinearAlgebra: cross, dot, norm
using Unitful
using Random
include("readData.jl") # Read field in from file
include("renders.jl") # Functions to generate outputs
include("interpolation.jl")

# -----------------------------------------------------------------------
mutable struct Particle
    """Properties of a single particle in the system."""
    r::Vector # position vector [x, y, z]
    v::Vector # velocity vector
    q # charge
    m # mass
end

# Create a way to generate random particles
Random.rand(rng::AbstractRNG, ::Random.SamplerType{Particle}) = Particle(
    (rand(rng, typeof(1.0u"m"), 3) - rand(rng, typeof(1.0u"m"), 3)) * 1e-7,
    rand(rng, typeof(1.0u"m/s"), 3) - rand(rng, typeof(1.0u"m/s"), 3),
    rand(rng, typeof(1.0u"C")) - rand(rng, typeof(1.0u"C")),
    rand(rng, typeof(1.0u"kg"))
)

# More limited random values, makes cooler plots
# Random.rand(rng::AbstractRNG, ::Random.SamplerType{Particle}) = Particle(
#     [0u"m", 0u"m", 0u"m"],
#     rand(rng, typeof(1.0u"m/s"), 3),
#     rand(rng, typeof(1.0u"C")),
#     rand(rng, typeof(1.0u"kg"))
# )

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

# Define the electric and magnetic fields.
# They are always defined as functions of position - even for
# constant fields - so that you don't have to go through changing
# the Boris method code whenever you switch to a position dependent field.

# The field functions can be anything of the form E(r::Vector)::Vector{T}
# The output is a vector of length 3 and in units apprioate to the field type
# There must be an E function for the electric field and a B function for the magnetic field

# Units are denoted using Unitful's notation.
# e.g. 1u"V/m" is 1 Volt per metre
# Just uncomment whichever block you want to use.

# Uniform fields
E(r::Vector)::Vector = [0u"V/m", 0u"V/m", 0.01u"V/m"]
B(r::Vector)::Vector = [0u"T", 0u"T", 1000000000u"T"]

# Load magnetic field in from HDF5 data
# Interpolation is slow so don't use it with lots of particles
# E(r::Vector)::Vector = [0u"V/m", 0u"V/m", 0.0u"V/m"]
# Bx, By, Bz = loadHDF5()
# # We'll make up the spatial dimensions, in metres
# x = range(-128, 128, 256) * 1e-8
# y = range(-128, 128, 256) * 1e-8
# z = range(-128, 128, 256) * 1e-8
# function B(r::Vector)::Vector
#     r = ustrip.(u"m", r)
#     fx = interpolate3D(Bx, x, y, z, r...)
#     fy = interpolate3D(By, x, y, z, r...)
#     fz = interpolate3D(Bz, x, y, z, r...)
#     return [fx, fy, fz] * 100000u"T" # Scale factor was chosen to get clear behaviour for default particle generation
# end

# Interpolate over a uniform field
# E(r::Vector)::Vector = [0u"V/m", 0u"V/m", 0.0u"V/m"]
# Bx = fill(100, (256,256,256))
# By = fill(100, (256,256,256))
# Bz = fill(100, (256,256,256))
# # We'll make up the spatial dimensions, in metres
# x = range(-128, 128, 256) * 1e-7
# y = range(-128, 128, 256) * 1e-7
# z = range(-128, 128, 256) * 1e-7
# function B(r::Vector)::Vector
#     r = ustrip.(u"m", r)
#     fx = interpolate3D(Bx, x, y, z, r...)
#     fy = interpolate3D(By, x, y, z, r...)
#     fz = interpolate3D(Bz, x, y, z, r...)
#     return [fx, fy, fz] * 100000u"T" # Scale factor was chosen to get clear behaviour for default particle generation
# end

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

function stepVelocity!(part::Particle, dt::Quantity)
    """Update velocity with change over timestep dt using Boris method."""
    halfEV = part.q / part.m * E(part.r) * dt / 2 # half electric acceleration

    t = part.q / part.m * B(part.r) * dt / 2
    s = 2 * t / (1 + dot(t, t)) # s= 2t/(1+|t|^2)

    vMinus = part.v + halfEV
    vPrime = vMinus + cross(vMinus, t)
    vPlus = vMinus + cross(vPrime, s)
    part.v = vPlus + halfEV
end

function stepPosition!(part::Particle, dt::Quantity)
    """Update position over timestep dt with current velocity."""
    part.r += part.v * dt
end

function lamorRadius(part::Particle)::Quantity
    """Calculate lamorRadius of particle."""
    # First want the tangential velocity
    n = B(part.r) / norm(B(part.r)) # unit vector of magnetic field
    vTang = norm(part.v - dot(part.v, n) * n) # v⟂ = v - v∥ = v - (v⋅n)n
    return part.m * vTang / part.q / norm(B(part.r))
end

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

# Init system
print("Initialising system.\n")
dt = 3e-11u"s" # timestep
iterations = 10000
# Define particles, in inital state
# The second argument to rand sets the number of particles.
particles = rand(Particle, 100)
print("There are $(length(particles)) particles.\n")
# Offset initial velocity back 1/2 step so it leapfrogs with position
stepVelocity!.(particles, -dt / 2)
# We'll store history in the following 3D array with indices [time, particle, coordinate]
positions = Array{Quantity,3}(undef, iterations, length(particles), 3)

# Main loop
print("Init complete. Beginning loop.\n")
for i in 1:iterations
    stepVelocity!.(particles, dt)
    stepPosition!.(particles, dt)
    positions[i, :, :] = [p.r[j] for p in particles, j in 1:3]
end
print("Loop complete.\n")

print("Making graphics.\n")
plotTrajectories(positions)
print("Close plot and press enter to load next plot.")
readline()
plotAtTime(positions, iterations, dt)
print("Close plot and press enter to load next plot.")
readline()
plotParticle(positions, 1)

# This is buggy but the output is the most interesting of the various options
# makeRender(iterations, particles, positions, size=100)

print("Graphics complete. Press enter to exit.")
readline() # Doesn't interact well with running in REPL.


text:
using LinearAlgebra
using Test

using SetProg, SetProg.Sets
using Polyhedra
using MultivariatePolynomials

using JuMP
const MOIT = MOI.Test

@testset "Controlled invariant" begin
    config = MOIT.Config()
    @testset "Ellipsoid" begin
        # Q = [1 0
        #      0 1]
        # t = √det(Q) = 1
        Q = [1.0, -1/4, 1.0]
        t = √15/4
        @testset "Homogeneous" begin
            Tests.ci_ell_homogeneous_test(bridged_mock(mock -> MOI.Utilities.mock_optimize!(mock, [Q; t])),
                                          config)
        end
        @testset "Non-homogeneous" begin
            @testset "Ellipsoid" begin
                β = -1.0
                b = [0.0, 0.0]
                Tests.ci_ell_nonhomogeneous_test(bridged_mock(mock -> MOI.Utilities.mock_optimize!(mock, [Q; β; b; t])),
                                                 config)
            end
            @testset "PolySet" begin
                β = -1.0
                b = [0.0, 0.0]
                Q = [1.0, -0.5, 1.0]

                Tests.ci_quad_nonhomogeneous_test(bridged_mock(mock -> begin
                         # β+1 b[1] b[2]
                         #  .  Q[1] Q[2]
                         #  .   .   Q[3]
                         MOI.Utilities.mock_optimize!(mock, [β+1; b[1]; Q[1]; b[2]; Q[2]; Q[3]; 2Q])
                     end),
                config)
            end
        end
        @testset "Piecewise" begin
            @testset "Diamond/2D cross-polytope" begin
                Q1 = [1.0, -0.25, 1.0]
                Q2 = [1.0, -1.0, 1.0]
                Tests.ci_piecewise_semiell_homogeneous_test(
                    bridged_mock(mock -> MOI.Utilities.mock_optimize!(
                        mock,
                        [Q1; Q2; Q2; Q1; collect(1:MOI.get(mock, MOI.NumberOfVariables()) - 12)])),
                    config)
            end
            @testset "mci" begin
                Q1 = [1.0, -1.0, 1.0]
                Q2 = [1.0, 0.0, 0.0]
                Q3 = [0.25, 0.5, 1.0]
                # FIXME See https://github.com/blegat/SetProg.jl/runs/4384813153?check_suite_focus=true
                if Sys.WORD_SIZE != 32
                    Tests.ci_piecewise_semiell_mci_homogeneous_test(
                        bridged_mock(mock -> MOI.Utilities.mock_optimize!(
                            mock,
                            [Q1; Q1; Q2; Q3; Q2; Q3; collect(1:MOI.get(mock, MOI.NumberOfVariables()) - 18)])),
                        config)
                end
            end
        end
    end
    @testset "Quartic" begin
        # The PSD matrix for the variable  is 3 x 3 so 3 * (3+1) / 2 = 6
        # The PSD matrix for the convexity is 6 x 6 so 6 * (6+1) / 2 = 21
        # entries
        # 1 variable for t
        # hence 28 variables
        ci_quartic_α = -1/8
        ci_quartic_β =  1/4
        ci_quartic_hess = 6 * [2.0, ci_quartic_α, 2.0, ci_quartic_α, 2.0,
                               2.0, 2.0, ci_quartic_α, ci_quartic_α, 2.0]
        sol = [1.0; ci_quartic_α; 6 - 2ci_quartic_β; ci_quartic_β; ci_quartic_α; 1.0;
               ci_quartic_hess]
        Tests.ci_quartic_homogeneous_test(bridged_mock(mock -> MOI.Utilities.mock_optimize!(mock, sol)),
                                          config)
    end
end