text:
module Scapin

using Base.Iterators
using LinearAlgebra

"""
    dimensionality(op)

Return the number of dimensions of the physical space that `op` operates on.
"""
function dimensionality end

"""
   eltype_real(type)

Return the type of the (scalar) elements the operator of given `type` operates
on in the *real* space (as opposed to the *Fourier* space). Note that this type
may be complex!

The definition `eltype_real(x) = eltype_real(typeof(x))` is provided for
convenience so that instances can be passed instead of types. However the form
that accepts a type argument should be defined for new types.
"""
function eltype_real end
eltype_real(op) = eltype_real(typeof(op))

"""
   eltype_fourier(type)

Return the type of the (scalar) elements the operator of given `type` operates
on in the *Fourier* space (as opposed to the *real* space). Note that this type
may be real!

The definition `eltype_fourier(x) = eltype_fourier(typeof(x))` is provided for
convenience so that instances can be passed instead of types. However the form
that accepts a type argument should be defined for new types.
"""
function eltype_fourier end
eltype_fourier(op) = eltype_fourier(typeof(op))

"""
    apply_fourier!(ŷ, ℱ, k, x̂) -> ŷ

Apply in-place operator `ℱ` to `x̂` in Fourier space and return the modified array `ŷ`.
More specifically,

    ŷ[:] = ℱ̂(k)⋅x̂,

where `k`, `x̂` and `ŷ` are column vectors: `k` is the vector of spatial
frequencies and `x̂` (resp. `ŷ`) is the Fourier mode of the input (resp. output),
for the spatial frequency `k`. The following must hold

    size(k, 1) == dimensionality(ℱ),
    size(x̂, 1) == size(ℱ, 2),
    size(ŷ, 1) == size(ℱ, 1).

The input and output vectors `x̂` and `ŷ`  are arrays of type `T` or `Complex{T}`,
where `T = eltype(ℱ)`.

!!! note "Continuous and discrete operators"

    For continuous operators, `k` would typically be a real-valued vector of wave
    numbers, while for discrete operators, `k` would be a integer-valued vector
    of indices.

!!! danger "Overriding the input vector"

    The present method *must* be implemented in such a way that aliasing `ŷ` with
    `x̂` is permitted. In other words, calling `apply_fourier!(x̂, ℱ, x̂, k)` *must*
    always deliver the correct answer.

"""
function apply_fourier! end

"""
    apply_fourier(ℱ, x̂, k)

Apply in-place operator `ℱ` to `x̂` in Fourier space and return the result.

See [apply_fourier!](@ref)

"""
function apply_fourier(ℱ, x̂, k)
    T = promote_type(eltype(ℱ), eltype(x̂), eltype(k))
    ŷ = Array{T}(undef, size(ℱ, 1))
    return apply_fourier!(ŷ, ℱ, x̂, k)
end

"""
    fourier_matrix(ℱ, k)

Return the `k`-th mode of `ℱ`, `ℱ̂(k)`, as a matrix.

!!! note "Performance of the default implementation"

    The default implementation relies on [apply_fourier!](@ref) to build the
    matrix column-by-column. It might be inefficient.
"""
function fourier_matrix(ℱ, k)
    nrows, ncols = size(ℱ)
    T = eltype_fourier(ℱ)
    mat = zeros(T, nrows, ncols)
    x̂ = zeros(T, ncols)
    for i = 1:ncols
        x̂[i] = one(T)
        apply_fourier!(view(mat, :, i), ℱ, k, x̂)
        x̂[i] = zero(T)
    end
    return mat

end

export dimensionality, eltype_real, eltype_fourier, apply_fourier!, apply_fourier, fourier_matrix

include("Elasticity.jl")
include("Grid.jl")
include("Brick.jl")
include("Bri17.jl")

end


text:
module PhysicsTutorials
# inspiration taken from DiffEqTutorials.jl

using Weave, Literate, Pkg

repo_directory = joinpath(@__DIR__,"..")
cssfile = joinpath(@__DIR__, "..", "templates", "skeleton_css.css")
latexfile = joinpath(@__DIR__, "..", "templates", "julia_tex.tpl")
htmlfile = joinpath(@__DIR__, "..", "templates", "julia_html.tpl")

function open_notebooks()
    Base.eval(Main, Meta.parse("import IJulia"))
    path = joinpath(repo_directory,"tutorials")
    IJulia.notebook(;dir=path)
end

function with_localenv(f, folder)
    activate_env(folder)
    f()
    activate_env(repo_directory)
end

function activate_env(folder)
    if isabspath(folder)
        Pkg.activate(folder)
    else
        Pkg.activate(joinpath(repo_directory, folder))
    end
    nothing
end

function format_header!(weavefile)
    lines = split(read(weavefile, String), "\n")
    idx = 1
    while strip(lines[idx]) == ""
        idx += 1
    end
    if startswith(lines[idx], "#") && startswith(lines[idx+1], "###")
        title = strip(lines[2][2:end])
        author = strip(lines[3][4:end])
        header = strip.(split("""---
        title: $title
        author: $author
        ---""", "\n"))
        write(weavefile, join(vcat(lines[1], header, lines[4:end]), "\n"))
    end
    nothing
end

abstract type Source end
struct NotebookSource <: Source end
struct LiterateSource <: Source end
struct WeaveSource <: Source end

convert_tutorial(cat, tut; kwargs...) = convert_tutorial(cat, tut, NotebookSource(); kwargs...)

function convert_tutorial(cat, tut, source::NotebookSource;
                            overwrite=false,
                            quick=false,
                            markdown=false,
                            keep_weave=false,
                            use_weave=false)
    tut_folder = joinpath(repo_directory, "tutorials", cat, tut)
    nbfile = joinpath(tut_folder, string(tut, ".ipynb"))
    weavefile = joinpath(tut_folder, string(tut, ".jmd"))

    # weave
    @info "Preparing"
    if use_weave
        !isfile(weavefile) && error("use_weave=true but no Weave file present!")
    else
        convert_doc(nbfile, weavefile)
    end

    # markdown
    if markdown && (!isfile(joinpath(tut_folder, string(tut, ".md"))) || overwrite)
        @info "Converting to markdown"
        if quick
            read(`jupyter nbconvert --to markdown $nbfile`)
        else
            with_localenv(tut_folder) do
                weave(weavefile, doctype="github")
            end
        end
    end

    format_header!(weavefile)
    # html
    if !isfile(joinpath(tut_folder, string(tut, ".html"))) || overwrite
        @info "Converting to html"
        if quick
            read(`jupyter nbconvert --to HTML $nbfile`)
        else
            with_localenv(tut_folder) do
                weave(weavefile, doctype="md2html"; css=cssfile, template=htmlfile)
            end
        end
    end
    # pdf
    if !isfile(joinpath(tut_folder, string(tut, ".pdf"))) || overwrite
        @info "Converting to pdf"
        if quick
            read(`jupyter nbconvert --to PDF $nbfile`)
        else
            with_localenv(tut_folder) do
                weave(weavefile, doctype="md2pdf"; template=latexfile)
            end
        end
    end
    keep_weave || (isfile(weavefile) && rm(weavefile))
    @info "Done"
    nothing
end

function convert_tutorial(cat, tut, source::LiterateSource; kwargs...)
    tut_folder = joinpath(repo_directory, "tutorials", cat, tut)
    literatefile = joinpath(tut_folder, string(tut, ".jl"))
    nbfile = joinpath(tut_folder, string(tut, ".ipynb"))

    # notebook
    @info "Converting to notebook"
    with_localenv(tut_folder) do
        Literate.notebook(literatefile, tut_folder; documenter=false)
    end
    convert_tutorial(cat, tut, NotebookSource(); kwargs...)
end

function convert_tutorial(cat, tut, source::WeaveSource; kwargs...)
    tut_folder = joinpath(repo_directory, "tutorials", cat, tut)
    weavefile = joinpath(tut_folder, string(tut, ".jmd"))
    nbfile = joinpath(tut_folder, string(tut, ".ipynb"))

    # notebook
    @info "Converting to notebook"
    with_localenv(tut_folder) do
        Weave.notebook(weavefile; out_path=tut_folder)
    end
    convert_tutorial(cat, tut, NotebookSource(); keep_weave=true, use_weave=true, kwargs...)
end

function _md2literate(s)
    # prepocessing
    s = replace(s, "\$\$ "=>"\$\$")
    s = replace(s, " \$\$"=>"\$\$")
    lines = split(s, "\n")
    # markdown to julia comments and julia blocks to plain julia
    juliablock = false
    for i in 1:length(lines)
        l = strip(lines[i])
        l == "" && continue
        if startswith(l, "```julia")
            juliablock = true
            lines[i] = ""
            continue
        elseif startswith(l, "```")
            juliablock = false
            lines[i] = ""
            continue
        end

        if !juliablock
            lines[i] = string("# ", l)
        else
            if startswith(l, "#") # julia code comment
                lines[i] = string("#", l)
            elseif startswith(l, "?")
                lines[i] = string("## ", l)
            end
        end
    end
    # remove subsequent empty lines
    out = String[]
    i = 1
    rmempty = false
    for i in 1:length(lines)
        l = lines[i]
        if l == "" && !rmempty
            rmempty = true
        elseif l =="" && rmempty
            continue # do not push line
        else
            rmempty = false
        end
        push!(out, l)
    end
    return join(out, "\n")
end

function nb2literate(cat, tut)
    tut_folder = joinpath(repo_directory, "tutorials", cat, tut)
    nbfile = joinpath(tut_folder, string(tut, ".ipynb"))
    mdfile = joinpath(tut_folder, string(tut, ".md"))
    jlfile = joinpath(tut_folder, string(tut, ".jl"))
    @info "Converting to markdown"
    read(`jupyter nbconvert --to markdown $nbfile`)
    mv(mdfile, jlfile, force=true)
    @info "Converting markdown to literate"
    s = read(jlfile, String)
    write(jlfile, _md2literate(s))
    asset_folder = joinpath(tut_folder, string(tut, "_files"))
    isdir(asset_folder) && rm(asset_folder, recursive=true)
    nothing
end

end # module

ext: jl

max_issues_repo_name: JuliaPhysics/PhysicsTutorials

max_forks_repo_name: JuliaPhysics/PhysicsTutorials


text:
using Documenter, GalacticOptim

makedocs(
    sitename="GalacticOptim.jl",
    authors="Chris Rackauckas, Vaibhav Kumar Dixit et al.",
    clean=true,
    doctest=false,
    modules=[GalacticOptim],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://galacticoptim.sciml.ai/stable/"),

    pages=[
        "GalacticOptim.jl: Unified Global Optimization Package" => "index.md",

        "Tutorials" => [
            "Basic usage" => "tutorials/intro.md",
            "Rosenbrock function" => "tutorials/rosenbrock.md",
            "Minibatch" => "tutorials/minibatch.md",
            "Symbolic Modeling" => "tutorials/symbolic.md"
        ],

        "API" => [
            "OptimizationProblem" => "API/optimization_problem.md",
            "OptimizationFunction" => "API/optimization_function.md",
            "solve" => "API/solve.md",
            "ModelingToolkit Integration" => "API/modelingtoolkit.md"
        ],

        "Local Optimizers" => [
            "local_optimizers/local_gradient.md",
            "local_optimizers/local_derivative_free.md",
            "local_optimizers/local_hessian.md",
            "local_optimizers/local_hessian_free.md"
        ],

        "Global Optimizers" => [
            "global_optimizers/global.md",
            "global_optimizers/global_constrained.md"
        ]
    ]
)

deploydocs(
    repo="github.com/SciML/GalacticOptim.jl";
    push_preview=true
)


text:
c3 = complement(g3)
c4 = complement(g4)

@test nv(c3) == 5
@test ne(c3) == 6
@test nv(c4) == 5
@test ne(c4) == 16

g = reverse(g4)
@test re1 in edges(g)
reverse!(g)
@test g == g4

g = blkdiag(g3, g3)
@test nv(g) == 10
@test ne(g) == 8

h = PathGraph(2)
@test intersect(g3, h) == h

h = PathGraph(4)
z = difference(g3, h)
@test nv(z) == 5
@test ne(z) == 1
z = difference(h, g3)
@test nv(z) == 4
@test ne(z) == 0

z = symmetric_difference(h,g3)
@test z == symmetric_difference(g3,h)
@test nv(z) == 5
@test ne(z) == 1

h = Graph(6)
add_edge!(h, 5, 6)
e = Edge(5, 6)
z = union(g3, h)
@test has_edge(z, e)
@test z == PathGraph(6)

h = DiGraph(6)
add_edge!(h, 5, 6)
e = Edge(5, 6)
z = union(g4, h)
@test has_edge(z, e)
@test z == PathDiGraph(6)

g10 = CompleteGraph(2)
h10 = CompleteGraph(2)
z = blkdiag(g10, h10)
@test nv(z) == nv(g10) + nv(h10)
@test ne(z) == ne(g10) + ne(h10)
@test has_edge(z, 1, 2)
@test has_edge(z, 3, 4)
@test !has_edge(z, 1, 3)
@test !has_edge(z, 1, 4)
@test !has_edge(z, 2, 3)
@test !has_edge(z, 2, 4)

g10 = Graph(2)
h10 = Graph(2)
z = join(g10, h10)
@test nv(z) == nv(g10) + nv(h10)
@test ne(z) == 4
@test !has_edge(z, 1, 2)
@test !has_edge(z, 3, 4)
@test has_edge(z, 1, 3)
@test has_edge(z, 1, 4)
@test has_edge(z, 2, 3)
@test has_edge(z, 2, 4)

p = PathGraph(10)
x = p*ones(10)
@test  x[1]==1.0 && all(x[2:end-1].==2.0) && x[end]==1.0

@test size(p) == (10,10)
@test size(p, 1) == size(p, 2) == 10
@test size(p, 3) == 1

@test g5 * ones(nv(g5)) == [2.0, 1.0, 1.0, 0.0]
@test sum(g5, 1) ==  [0, 1, 2, 1]
@test sum(g5, 2) ==  [2, 1, 1, 0]
@test sum(g5) == 4
@test sum(p,1) == sum(p,2)
@test_throws ErrorException sum(p,3)

@test sparse(p) == adjacency_matrix(p)
@test eltype(p) == Float64
@test length(p) == 100
@test ndims(p) == 2
@test issymmetric(p)
@test !issymmetric(g5)

g22 = CompleteGraph(2)
h = cartesian_product(g22, g22)
@test nv(h) == 4
@test ne(h)== 4

g22 = CompleteGraph(2)
h = tensor_product(g22, g22)
@test nv(h) == 4
@test ne(h) == 1

nx = 20; ny = 21
G = PathGraph(ny); H = PathGraph(nx)
c = cartesian_product(G, H)
g = crosspath(ny, PathGraph(nx));
@test g == c

function crosspath_slow(len, h)
    g = h
    m = nv(h)
    for i in 1:len-1
        k = nv(g)
        g = blkdiag(g,h)
        for v in 1:m
            add_edge!(g, v+(k-m), v+k)
        end
    end
    return g
end
@test crosspath_slow(2, g22) == crosspath(2,g22)


text:
@everywhere module KernelTrick
import Combinatorics: combinations

export bandwidth, GaussianKernelInRKHS_M, κ_KAMH

function bandwidth(particles)
    """
    median value of the norms.
    In the SMC-ABC context, unique function is instructive to reduce the computation cost.
    """
    values_list = [particle.x_present for particle in particles]
    params = hcat(values_list...)'
    params = unique(params, 1)
    N = size(params)[1]
    bw = 0.0
    if N > 3
        combs = hcat(collect(combinations(collect(1:N),2))...)'
        num_comb = size(combs)[1]
        norms = zeros(num_comb)
        for i=1:num_comb
            norms[i] = norm(params[combs[i,1],:]-params[combs[i,2],:])
        end
        bw = mean(norms)
    else
        bw = -1
    end

    return bw
end

function GaussianKernelInRKHS_M(particles, y, σ)
    """
    y:: current state particle
    particles:: particle subset of history (all particles in SMC context)
    """
    κ(x1, x2) = exp(-norm(x1-x2)^2 / (2σ^2))
    ∇κ(x1, x2) = κ(x1, x2)*(x2-x1) / σ^2
    N = length(particles)
    M = zeros(N, length(particles[1].x_present))
    for i=1:N
        M[i,:] .= 2*∇κ(y.x_present, particles[i].x_present)
    end
    return M
end

function κ_KAMH(particles, y, C, M; γ=1., nu=1.)
    """
    proposed function
    """
    _Σ = γ^2 * eye(size(M)[2]) + nu^2 * M'*C*M
    Σ = Hermitian(_Σ)
    return Σ
end

end


text:
@benchmark random_phase_check!(ψ, A, 𝓀ˣ, 𝓀ʸ, x, y, φ)

@benchmark begin
    P⁻¹ * ψ
    P⁻¹ * θ
    P⁻¹ * u
    P⁻¹ * v
    P⁻¹ * ∂ˣθ
    P⁻¹ * ∂ʸθ
    P⁻¹ * κΔθ
end

@benchmark begin
    @. u = filter * -1.0 * (∂y * ψ)
    @. v = filter * (∂x * ψ)
    # ∇θ
    @. ∂ˣθ = filter * ∂x * θ
    @. ∂ʸθ = filter * ∂y * θ
    @. κΔθ = κ * Δ * θ
end

@benchmark begin
    φ̇ .= arraytype((rand(size(A)...) .- 0.5))
    @. θ̇ = -u * ∂ˣθ - v * ∂ʸθ + κΔθ + s
end

@benchmark begin
    @. φ += sqrt(Δt) * φ̇
    @. θ += Δt * θ̇
    t[1] += Δt
end

@benchmark begin
    # spectral space representation
    P * ψ # in place fft
    P * θ # in place fft
    # ∇ᵖψ
    @. u = filter * -1.0 * (∂y * ψ)
    @. v = filter * (∂x * ψ)
    # ∇θ
    @. ∂ˣθ = filter * ∂x * θ
    @. ∂ʸθ = filter * ∂y * θ
    @. κΔθ = κ * Δ * θ
    # go back to real space
    P⁻¹ * ψ
    P⁻¹ * θ
    P⁻¹ * u
    P⁻¹ * v
    P⁻¹ * ∂ˣθ
    P⁻¹ * ∂ʸθ
    P⁻¹ * κΔθ
    # Assemble RHS
    φ̇ .= arraytype((rand(size(A)...) .- 0.5))
    @. θ̇ = -u * ∂ˣθ - v * ∂ʸθ + κΔθ + s
    # Euler step
    @. φ += sqrt(Δt) * φ̇
    @. θ += Δt * θ̇
    t[1] += Δt
end

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
using BenchmarkTools, StatsBase

function summarize(results, benchmark_set)

  for m in benchmark_set
    file = open("./stochkit/$(m).txt", "a")

    result = results[m]
    adjustment = adjust[m]

    for i in eachindex(result.times)
      # convert to ms and correct for overhead
      result.times[i] = result.times[i] * 1e-6 - adjustment
    end

    println(file, "model: $(m)")
    println(file)
    println(file, "n = $(length(result.times))")
    println(file, "mean:               $(mean(result.times)) ms")
    println(file, "std:                $(std(result.times)) ms")
    println(file, "median:             $(median(result.times)) ms")
    println(file, "mad (normal=true):  $(mad(result.times, normalize=true)) ms")
    println(file, "mad (normal=false): $(mad(result.times, normalize=false)) ms")
    println(file)

    close(file)
  end

  return nothing
end

# initialize test suite
suite = BenchmarkGroup()

# load the following variables:
# - parameters
# - MODELS
# - SEED
include("benchmark_parameters.jl")

# build handy collections
const benchmark_set = intersect(ARGS, MODELS)
const adjust = Dict(
  "kendall"     => 80.0,
  "mmek"        => 100.0,
  "autoreg"     => 115.0,
  "dimer-decay" => 100.0,
  "yeast"       => 130.0
)

# build up the benchmark suite
for m in benchmark_set
  # unpack parameters
  t_final, n_saves, n_trial, n_sample, t_limit = parameters[m]

  # build filename
  fname = "./$(m)/stochkit2.xml"

  # build command
  cmd = `ssa -m $(fname) -t $(t_final) -i $(n_saves) -r $(n_trial) -f --seed $(SEED) '&>' /dev/null`

  suite[m] = @benchmarkable run($cmd) samples=n_sample seconds=t_limit
end

results = run(suite, verbose = true)

# perform adjustments
for b in leaves(results)
  b_info, trial = b

  m = b_info[1]

  adjustment = adjust[m]

  for i in eachindex(trial.times)
    trial.times[i] -= adjustment
  end
end

# summarize(results, benchmark_set)
BenchmarkTools.save("stochkit.json", results)


text:
"""
    grdmath(cmd::String, args...)

Call grdmath with all commands in a single string 'cmd'.
This is not useful in itself as compared to call gmt("grdmath ....") but it's very useful
in 'movie' because it can generate shell scripts from the julai command

Full option list at [`grdmath`]($(GMTdoc)grdmath.html)
"""
function grdmath(cmd::String, args...)

	d = KW()
	cmd = "grdmath " * cmd
	(dbg_print_cmd(d, cmd) !== nothing) && return cmd		# The real useful bit od this fun
	gmt(cmd, args...)
end

"""
    gmtmath(cmd::String, args...)

Call gmtmath with all commands in a single string 'cmd'.
This is not useful in itself as compared to call gmt("gmtmath ....") but it's very useful
in 'movie' because it can generate shell scripts from the julai command

Full option list at [`grdmath`]($(GMTdoc)grdmath.html)
"""
function gmtmath(cmd::String, args...)

	d = KW()
	cmd = "gmtmath " * cmd
	(dbg_print_cmd(d, cmd) !== nothing) && return cmd
	gmt(cmd, args...)
end


text:
using jInv.InverseSolve
export getHessian

"""
function H = getHessian(...)

builds and returns Hessian matrix of misfit functions.

Input:

	sig   - model
	pMis  - misfit param (supported: single, array of misfit params, remote references)
	d2F   - Hessian of distance (supported: single, array, remote references)

Output:

	H     - Hessian

Examples:

	H  = getHessian(sigma, pMis, d2F)
	Ha = getHessian(sigma, [pMis1;pMis2], [d2F1; d2F2])

"""
function getHessian(sig::Vector,  # conductivity on inv mesh (active cells only)
                    pMis::MisfitParam,
                    d2F)
	sigma,dsigma = pMis.modelfun(sig)

	P = pMis.gloc.PForInv
	J = getSensMat(sigma,pMis.pFor)

	dr = J*P*dsigma

	return dr'*sdiag(d2F)*dr
end

function getHessian(sig,pMis::Array,d2F::Array,workerList=workers())

	Hs = Array{Any}(undef, length(pMis))
	i=1; nextidx() = (idx = i; i+=1; idx)

	workerList = intersect(workers(),workerList)
	if isempty(workerList)
		error("getHessian: specified workers do not exist!")
	end

	sigRef = Array{Future}(undef, maximum(workerList))

	@sync begin
		for p=workerList
			@async begin
				sigRef[p] = remotecall(identity,p,sig)
				while true
					idx = nextidx()
					if idx > length(pMis)
						break
					end
					Hs[idx]    = remotecall_fetch(getHessian,p,sigRef[p],pMis[idx],d2F[idx])
				end
			end
		end
	end
	H = Hs[1]
	for k=2:length(Hs)
		H += Hs[k]
	end

	return H
end

getHessian(sig::Union{RemoteChannel,Future}, pMis::Union{RemoteChannel,Future}, d2F::Union{RemoteChannel,Future}) = getHessian(fetch(fetch(sig)),fetch(pMis),fetch(d2F))
getHessian(sig::Array, pMis::Union{RemoteChannel,Future}, d2F::Union{RemoteChannel,Future}) = getHessian(sig,fetch(pMis),fetch(d2F))
getHessian(sig::Union{RemoteChannel,Future}, pMis::MisfitParam, d2F::Array) = getHessian(fetch(sig), pMis, d2F)


text:
module CellNumberApply

using Gridap
using Gridap.Helpers

export apply
import Base: iterate
import Base: length
import Base: size
import Base: getindex
import Base: IndexStyle

function apply(k::NumberKernel,v::Vararg{<:CellValue})
  CellNumberFromKernel(k,v...)
end

function apply(k::NumberKernel,v::Vararg{<:IndexCellValue})
  IndexCellNumberFromKernel(k,v...)
end

function apply(f::Function,v::Vararg;broadcast=false)
  _apply(f,v,Val(broadcast))
end

function _apply(f,v,::Val{false})
  t = _eltypes(v)
  T = Base._return_type(f,t)
  @assert T <: NumberLike
  k = NumberKernelFromFunction(f)
  apply(k,v...)
end

function _eltypes(v)
  tuple([ eltype(vi) for vi in v ]...)
end

struct CellNumberFromKernel{T,K,V} <: IterCellNumber{T}
  kernel::K
  cellvalues::V
end

function CellNumberFromKernel(k::NumberKernel,v::Vararg{<:CellValue})
  _checks(v)
  T = _compute_type(k,v)
  K = typeof(k)
  V = typeof(v)
  CellNumberFromKernel{T,K,V}(k,v)
end

function _checks(v)
  @assert length(v) > 0
  v1, = v
  l = length(v1)
  @assert all( ( length(vi) == l for vi in v ) )
end

function _compute_type(k,v)
  t = _eltypes(v)
  T = compute_type(k,t...)
  @assert T <: NumberLike
  T
end

function length(self::CellNumberFromKernel)
  vi, = self.cellvalues
  length(vi)
end

@inline function iterate(self::CellNumberFromKernel)
  zipped = zip(self.cellvalues...)
  znext = iterate(zipped)
  _iterate(self,znext,zipped)
end

@inline function iterate(self::CellNumberFromKernel,state)
  zipped, zstate = state
  znext = iterate(zipped,zstate)
  _iterate(self,znext,zipped)
end

@inline function _iterate(self::CellNumberFromKernel,znext,zipped)
  if znext === nothing; return nothing end
  a, zstate = znext
  r = compute_value(self.kernel,a...)
  state = (zipped,zstate)
  (r, state)
end

struct IndexCellNumberFromKernel{T,K,V} <: IndexCellNumber{T,1}
  kernel::K
  cellvalues::V
end

function IndexCellNumberFromKernel(k::NumberKernel,v::Vararg{<:IndexCellValue})
  _checks(v)
  T = _compute_type(k,v)
  K = typeof(k)
  V = typeof(v)
  IndexCellNumberFromKernel{T,K,V}(k,v)
end

function length(self::IndexCellNumberFromKernel)
  vi, = self.cellvalues
  length(vi)
end

size(self::IndexCellNumberFromKernel) = (length(self),)

function getindex(self::IndexCellNumberFromKernel,i::Integer)
  vals = _getvalues(i,self.cellvalues...)
  compute_value(self.kernel,vals...)
end

# TODO use a generated function?
_getvalues(i,v...) = @notimplemented
_getvalues(i,v1) = (v1[i],)
_getvalues(i,v1,v2) = (v1[i],v2[i])
_getvalues(i,v1,v2,v3) = (v1[i],v2[i],v3[i])
_getvalues(i,v1,v2,v3,v4) = (v1[i],v2[i],v3[i],v4[i])
_getvalues(i,v1,v2,v3,v4,v5) = (v1[i],v2[i],v3[i],v4[i],v5[i])
_getvalues(i,v1,v2,v3,v4,v5,v6) = (v1[i],v2[i],v3[i],v4[i],v5[i],v6[i])

end # module CellNumberApply


text:
# Will update interface utilities with upgrade to Generalized parametric bisection

#=
function kdel(i,j)
  i == j ? 1.0 : 0.0
end

function Expr_SubSymbols(expr,symdict)
      stack = []
      exprlist = []
      push!(stack,expr)
      cont_flag = false
      while length(stack)>0
            a = pop!(stack)
            if (isempty(exprlist))
                  push!(exprlist,a)
            end
            if ~(typeof(a)<:Number) && ~(typeof(a)==Int64)
              for i=2:length(a.args)
                  if (typeof(a.args[i])==Symbol)
                        a.args[i] = symdict[a.args[i]]
                  elseif typeof(a.args[i])==Expr
                        push!(stack,a.args[i])
                  end
              end
            end
      end
      return exprlist[1]
end

function GenerateJacobianX(ExprArr,SymXArr,SymPArr)

  # creates symbol dictionary
  symdict = Dict{Any,Any}()
  nc = length(SymXArr)
  ns = length(SymPArr)
  for i=1:nc
      symdict[SymXArr[i]] = :(z[$i])
  end
  for i=1:ns
      symdict[SymPArr[i]] = :(p[$i])
  end

  # creates calculates Jacobian symbolically and subs in symbols
  StringArr = string.(ExprArr)
  dexpr_arr = Expr(:vcat)
  for i = 1:length(StringArr)
    out = Calculus.differentiate(StringArr[i],SymXArr)
    expr_out = copy(out)
    for j = 1:length(out)
      expr_out[j] = Expr_SubSymbols(out[j],symdict)
    end
    temp_expr = Expr(:row)
    for j = 1:length(out)
      push!(temp_expr.args,out[j])
    end
    push!(dexpr_arr.args,temp_expr) # populations Jacobian storage
  end

  # puts expression object into function
  @eval Jacob(z,p) = $dexpr_arr
  return Jacob
end

function GenerateH(ExprArr, SymXArr,SymPArr)
    # creates list of symbols
    symdict = Dict{Any,Any}()
    for i=1:length(SymXArr)
        symdict[SymXArr[i]] = :(x[$i])
    end
    for i=1:length(SymPArr)
        symdict[SymPArr[i]] = :(p[$i])
    end
    # creates array of expressions
    htemp1 = deepcopy(ExprArr)
    htemp2 = []
    for i = 1:length(ExprArr)
        push!(htemp2,Expr_SubSymbols(htemp1[i],symdict))
    end
    #
    hexpr = Expr(:vect)
    for i = 1:length(htemp2)
        temp = htemp2[i]
        push!(hexpr.args,Expr(:call,:eval,temp))
    end
    @eval hnew(x,p) = $hexpr
    return hnew
end

# Unpack function
function Unpack(stack)
  cnode = pop!(stack)
  return cnode[1],cnode[2]
end
=#

"""
    Strict_XinY(X::Vector{Interval{T}},Y::Vector{Interval{T}})

Returns true if X is strictly in Y (X.lo>Y.lo && X.hi<Y.hi).
"""
function Strict_XinY(X::Vector{Interval{T}},Y::Vector{Interval{T}}) where {T<:AbstractFloat}
  k::Bool = true
  for i=1:length(X)
    if ((X[i].lo<=Y[i].lo)||
        (X[i].hi>=Y[i].hi))
      k = false
    end
  end
  return k
end
function Strict_XinY(X::Vector{MCInterval{T}},Y::Vector{MCInterval{T}}) where {T<:AbstractFloat}
  k::Bool = true
  for i=1:length(X)
    if ((X[i].lo<=Y[i].lo)||
        (X[i].hi>=Y[i].hi))
      k = false
    end
  end
  return k
end

"""
    Strict_XinY(X::Interval{T},Y::Interval{T})

Returns true if X is strictly in Y (X.lo>Y.lo && X.hi<Y.hi).
"""
function Strict_XinY(X::Interval{T},Y::Interval{T}) where {T<:AbstractFloat}
  (X.lo<=Y.lo) && return false
  (X.hi>=Y.hi) && return false
  return true
end
function Strict_XinY(X::MCInterval{T},Y::MCInterval{T}) where {T<:AbstractFloat}
  (X.lo<=Y.lo) && return false
  (X.hi>=Y.hi) && return false
  return true
end

"""
    isEqual(X1::Vector{Interval{T}},X2::Vector{Interval{T}},atol::Float64)

Returns true if X1 and X2 are equal to within tolerance atol in all dimensions.
"""
function isEqual(X1::Vector{Interval{T}},X2::Vector{Interval{T}},
                 atol::Float64) where {T<:AbstractFloat}
  out::Bool = true
  for i=1:length(X1)
    if (abs(X1[i].lo-X2[i].lo)>=atol ||
        abs(X1[i].hi-X2[i].hi)>=atol )
        out = false
        break
    end
  end
  return out
end
function isEqual(X1::Vector{MCInterval{T}},X2::Vector{MCInterval{T}},
                 atol::Float64) where {T<:AbstractFloat}
  out::Bool = true
  for i=1:length(X1)
    if (abs(X1[i].lo-X2[i].lo)>=atol ||
        abs(X1[i].hi-X2[i].hi)>=atol )
        out = false
        break
    end
  end
  return out
end

"""
    extDivide(A::Interval{T},B::Interval{T},C::Interval{T})

Subfunction to generate output for extended division.
"""
function extDivide(A::Interval{T}) where {T<:AbstractFloat}
  #println("start extdivide")
  if ((A.lo == -0.0) && (A.hi == 0.0))
    B::Interval{T} = Interval(-Inf,Inf)
    C::Interval{T} = B
    return 0,B,C
  end
  if (A.lo == 0.0)
    B = Interval(1.0/A.hi,Inf)
    C = Interval(Inf,Inf)
    return 1,B,C
  elseif (A.hi == 0.0)
    B = Interval(-Inf,1.0/A.lo)
    C = Interval(-Inf,-Inf)
    return 2,B,C
  else
    B = Interval(-Inf,1.0/A.lo)
    C = Interval(1.0/A.hi,Inf)
    return 3,B,C
  end
end
function extDivide(A::MCInterval{T}) where {T<:AbstractFloat}
  #println("start extdivide")
  if ((A.lo == -0.0) && (A.hi == 0.0))
    B::MCInterval{T} = MCInterval{T}(-Inf,Inf)
    C::MCInterval{T} = B
    return 0,B,C
  end
  if (A.lo == 0.0)
    B = MCInterval{T}(1.0/A.hi,Inf)
    C = MCInterval{T}(Inf,Inf)
    return 1,B,C
  elseif (A.hi == 0.0)
    B = MCInterval{T}(-Inf,1.0/A.lo)
    C = MCInterval{T}(-Inf,-Inf)
    return 2,B,C
  else
    B = MCInterval{T}(-Inf,1.0/A.lo)
    C = MCInterval{T}(1.0/A.hi,Inf)
    return 3,B,C
  end
end

"""
    extProcess(N::Interval{T},X::Interval{T},Mii::Interval{T},
                    S1::Interval{T},S2::Interval{T},B::Interval{T},rtol::Float64)

Generates output boxes for extended division and flag.
"""
# Error & Case Handling for Extended IA in NewtonGS (CHECKED C++)
function extProcess(N::Interval{T},X::Interval{T},Mii::Interval{T},
                    S1::Interval{T},S2::Interval{T},B::Interval{T},rtol::Float64) where {T<:AbstractFloat}
  v = 1
  Ntemp::Interval{T} = copy(N)
  M::Interval{T} = (B+S1+S2)+Interval(-rtol,rtol)
  if (M.lo<=0 && M.hi>=0)
    return 0, Interval(-Inf,Inf), Ntemp
  end
  if (v == 1)
    k,IML::Interval{T},IMR::Interval{T} = extDivide(Mii)
    if (k == 1)
      return 0, (mid(X)-M*IML), Ntemp
    elseif (k == 2)
      return 0, (mid(X)-M*IMR), Ntemp
    elseif (k == 3)
      NR = mid(X)-M*IMR
      NL = mid(X)-M*IML
      if (~isdisjoint(NL,X) && isdisjoint(NR,X))
        return 0, NL, Ntemp
      elseif (~isdisjoint(NR,X) && isdisjoint(NL,X))
        return 0, NR, Ntemp
      elseif (~isdisjoint(NL,X) && ~isdisjoint(NR,X))
        N = NL
        Ntemp = NR
        return 1, NL, NR
      else
        return -1, N, Ntemp
      end
    end
  end
  return 0, N, Ntemp
end

function extProcess(N::MCInterval{T},X::MCInterval{T},Mii::MCInterval{T},
                    S1::MCInterval{T},S2::MCInterval{T},B::MCInterval{T},rtol::Float64) where {T<:AbstractFloat}
  v = 1
  Ntemp::MCInterval{T} = copy(N)
  M::MCInterval{T} = (B+S1+S2)+MCInterval(-rtol,rtol)
  if (M.lo<=0 && M.hi>=0)
    return 0, MCInterval(-Inf,Inf), Ntemp
  end
  if (v == 1)
    k,IML::MCInterval{T},IMR::MCInterval{T} = extDivide(Mii)
    if (k == 1)
      return 0, (mid(X)-M*IML), Ntemp
    elseif (k == 2)
      return 0, (mid(X)-M*IMR), Ntemp
    elseif (k == 3)
      NR = mid(X)-M*IMR
      NL = mid(X)-M*IML
      if (~isdisjoint(NL,X) && isdisjoint(NR,X))
        return 0, NL, Ntemp
      elseif (~isdisjoint(NR,X) && isdisjoint(NL,X))
        return 0, NR, Ntemp
      elseif (~isdisjoint(NL,X) && ~isdisjoint(NR,X))
        N = NL
        Ntemp = NR
        return 1, NL, NR
      else
        return -1, N, Ntemp
      end
    end
  end
  return 0, N, Ntemp
end


text:
using QuantumStateDistributions
using QuantumStateBase
using LinearAlgebra
using Test

const QSB = QuantumStateBase
const QSD = QuantumStateDistributions

@testset "QuantumStateDistributions.jl" begin
    include("pdf.jl")
    include("mle.jl")
    include("gaussian_state_bhd.jl")
end


text:
using FactorPostProcessing, Test, Random
using BeastUtils.Logs, LinearAlgebra

Random.seed!(666)

const TOL = 1e-10
const POS_COLS = [9, 5, 13]

function absmax(x::AbstractArray{T}, y::AbstractArray{T}) where T <: Real
    return maximum(abs.(x - y))
end

function svd_process(L_data::Matrix{Float64}, F_data::Matrix{Float64}, row::Int, k::Int, p::Int, n::Int)
    L_vec = @view L_data[row, :]
    L = Matrix(reshape(L_vec, p, k))
    Lsvd = svd(L)
    F_vec = @view F_data[row, :]
    F = reshape(F_vec, k, n)
    F_rotated = Lsvd.Vt * F
    S = Lsvd.S
    Vt = Lsvd.Vt
    U = Lsvd.U
    F_rotated = Vt * F
    return (sv = Lsvd.S, V = Lsvd.U, F_rotated)
end

states = 1000
K = 3
P = 20
N = 100

L_data = randn(states, K * P)
F_data = randn(states, K * N)

L_labels = vec(["L$j$i" for i = 1:K, j=1:P])
F_labels = vec(["factors.taxon_$i.$j" for j=1:K, i = 1:N])

cols = [["posterior", "prior", "likelihood"];
        ["factorPrecision$i" for i = 1:P]; L_labels; F_labels]
data = [randn(states, 3 + P) L_data F_data]

log_path = joinpath(@__DIR__, "testFactorPostProcessing.log")
svd_path = joinpath(@__DIR__, "testFactorPostProcessingSVD.log")

@assert !isfile(log_path)
@assert !isfile(svd_path)

Logs.make_log(log_path, data, cols)
svd_logs(log_path, svd_path, rotate_factors=true)

svd_cols, svd_data = get_log(svd_path, burnin=0.0)

rm(log_path)
rm(svd_path)

L_cols = findall(startswith("L"), svd_cols)
V_cols = findall(startswith("V"), svd_cols)
sv_cols = findall(startswith("sv"), svd_cols)
F_cols = findall(startswith("factors"), svd_cols)

for state = 1:states
    L_actual = reshape(svd_data[state, L_cols], P, K)
    V_actual = reshape(svd_data[state, V_cols], P, K)
    sv_actual = svd_data[state, sv_cols]
    F_actual = reshape(svd_data[state, F_cols], K, N)
    S_expected, V_expected, F_expected = svd_process(L_data, F_data, state, K, P, N)
    L_expected = V_expected * Diagonal(S_expected)
    @test absmax(sv_actual, S_expected) < TOL

    signs = ones(K)
    for k = 1:K
        if L_expected[POS_COLS[k], k] < 0.0
            signs[k] = -1.0
        end
    end

    D = Diagonal(signs)
    L_expected = L_expected * D
    V_expected = V_expected * D
    F_expected = D * F_expected

    @test absmax(V_actual, V_expected) < TOL
    @test absmax(L_actual, L_expected) < TOL
    @test absmax(F_actual, F_expected) < TOL
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "OpenLDAPBuilder"
version = v"2.4.46"

# Collection of sources required to build OpenLDAPBuilder
sources = [
    "ftp://ftp.openldap.org/pub/OpenLDAP/openldap-release/openldap-2.4.46.tgz" =>
    "9a90dcb86b99ae790ccab93b7585a31fbcbeec8c94bf0f7ab0ca0a87ea0c4b2d",

    "patches",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd openldap-2.4.46/
#wget http://www.linuxfromscratch.org/patches/blfs/svn/openldap-2.4.46-consolidated-1.patch
patch -Np1 -i ../openldap_patches/openldap-2.4.46-consolidated-1.patch
autoconf
./configure --prefix=$prefix --host=$target --with-yielding-select=yes --disable-static --enable-dynamic --disable-debug --disable-slapd
apk update
apk add groff binutils
make depend
make -j${nproc}
make install || true # This is a hack

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:x86_64, :glibc),
    MacOS(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "liblber", :liblber),
    LibraryProduct(prefix, "libldap", :libldap),
    LibraryProduct(prefix, "libldap_r", :libldap_r)
]

# 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: jpsamaroo/OpenLDAPBuilder

max_issues_repo_name: jpsamaroo/OpenLDAPBuilder

max_forks_repo_name: jpsamaroo/OpenLDAPBuilder


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

@doc raw"""Represents a source location of a volume to mount, managed by an external CSI driver

    IoK8sApiCoreV1CSIVolumeSource(;
        driver=nothing,
        fsType=nothing,
        nodePublishSecretRef=nothing,
        readOnly=nothing,
        volumeAttributes=nothing,
    )

    - driver::String : Driver is the name of the CSI driver that handles this volume. Consult with your admin for the correct name as registered in the cluster.
    - fsType::String : Filesystem type to mount. Ex. \"ext4\", \"xfs\", \"ntfs\". If not provided, the empty value is passed to the associated CSI driver which will determine the default filesystem to apply.
    - nodePublishSecretRef::IoK8sApiCoreV1LocalObjectReference : NodePublishSecretRef is a reference to the secret object containing sensitive information to pass to the CSI driver to complete the CSI NodePublishVolume and NodeUnpublishVolume calls. This field is optional, and  may be empty if no secret is required. If the secret object contains more than one secret, all secret references are passed.
    - readOnly::Bool : Specifies a read-only configuration for the volume. Defaults to false (read/write).
    - volumeAttributes::Dict{String, String} : VolumeAttributes stores driver-specific properties that are passed to the CSI driver. Consult your driver's documentation for supported values.
"""
mutable struct IoK8sApiCoreV1CSIVolumeSource <: SwaggerModel
    driver::Any # spec type: Union{ Nothing, String } # spec name: driver
    fsType::Any # spec type: Union{ Nothing, String } # spec name: fsType
    nodePublishSecretRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1LocalObjectReference } # spec name: nodePublishSecretRef
    readOnly::Any # spec type: Union{ Nothing, Bool } # spec name: readOnly
    volumeAttributes::Any # spec type: Union{ Nothing, Dict{String, String} } # spec name: volumeAttributes

    function IoK8sApiCoreV1CSIVolumeSource(;driver=nothing, fsType=nothing, nodePublishSecretRef=nothing, readOnly=nothing, volumeAttributes=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1CSIVolumeSource, Symbol("driver"), driver)
        setfield!(o, Symbol("driver"), driver)
        validate_property(IoK8sApiCoreV1CSIVolumeSource, Symbol("fsType"), fsType)
        setfield!(o, Symbol("fsType"), fsType)
        validate_property(IoK8sApiCoreV1CSIVolumeSource, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        setfield!(o, Symbol("nodePublishSecretRef"), nodePublishSecretRef)
        validate_property(IoK8sApiCoreV1CSIVolumeSource, Symbol("readOnly"), readOnly)
        setfield!(o, Symbol("readOnly"), readOnly)
        validate_property(IoK8sApiCoreV1CSIVolumeSource, Symbol("volumeAttributes"), volumeAttributes)
        setfield!(o, Symbol("volumeAttributes"), volumeAttributes)
        o
    end
end # type IoK8sApiCoreV1CSIVolumeSource

const _property_map_IoK8sApiCoreV1CSIVolumeSource = Dict{Symbol,Symbol}(Symbol("driver")=>Symbol("driver"), Symbol("fsType")=>Symbol("fsType"), Symbol("nodePublishSecretRef")=>Symbol("nodePublishSecretRef"), Symbol("readOnly")=>Symbol("readOnly"), Symbol("volumeAttributes")=>Symbol("volumeAttributes"))
const _property_types_IoK8sApiCoreV1CSIVolumeSource = Dict{Symbol,String}(Symbol("driver")=>"String", Symbol("fsType")=>"String", Symbol("nodePublishSecretRef")=>"IoK8sApiCoreV1LocalObjectReference", Symbol("readOnly")=>"Bool", Symbol("volumeAttributes")=>"Dict{String, String}")
Base.propertynames(::Type{ IoK8sApiCoreV1CSIVolumeSource }) = collect(keys(_property_map_IoK8sApiCoreV1CSIVolumeSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1CSIVolumeSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1CSIVolumeSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1CSIVolumeSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1CSIVolumeSource[property_name]

function check_required(o::IoK8sApiCoreV1CSIVolumeSource)
    (getproperty(o, Symbol("driver")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1CSIVolumeSource }, name::Symbol, val)
end


text:
@doc """

""" ->
function lfS_test()

  return lfS_stat, lfS_pvalue
end


text:
using Memoize, LRUCache
using SpinGlassPEPS, TensorOperations, TensorCast
using LinearAlgebra

@memoize function left_env2(ϕ::AbstractMPS, idx::NTuple)
    l = length(idx)
    if l == 0
        L = [1.]
    else
        m = idx[l]
        new_idx = idx[1:l-1]
        L_old = left_env2(ϕ, new_idx)
        M = ϕ[l]
        @reduce L[x] := sum(α) L_old[α] * M[α, $m, x]
    end
    return L
end

@memoize function right_env2(ϕ::AbstractMPS, W::AbstractMPO, idx::NTuple)
    l = length(idx)
    L = length(ϕ)
    if l == 0
        R = fill(1., 1, 1)
    else
        m = idx[1]
        new_idx = idx[2:l]
        R_old = right_env2(ϕ, W, new_idx)
        M = ϕ[L-l+1]
        M̃ = W[L-l+1]
        @reduce R[x, y] := sum(α, β, γ) M̃[y, $m, β, γ] * M[x, γ, α] * R_old[α, β]
    end
    return R
end

@memoize function left_env3(ϕ::AbstractMPS, ψ::AbstractMPS)
    l = length(ψ)
    T = promote_type(eltype(ψ), eltype(ϕ))

    L = Vector{Matrix{T}}(undef, l+1)
    L[1] = ones(eltype(ψ), 1, 1)

    for i ∈ 1:l
        M = ψ[i]
        M̃ = conj.(ϕ[i])

        C = L[i]
        @tensor C[x, y] := M̃[β, σ, x] * C[β, α] * M[α, σ, y] order = (α, β, σ)
        L[i+1] = C
    end
    L
end

ϕ = randn(MPS{Float64}, 10, 10, 10);
ψ = randn(MPS{Float64}, 10, 10, 10);
W = randn(MPO{Float64}, 10, 10, 10);


text:
mutable struct RunningStats
  num_infected::Int
  num_detected::Int
  num_dead::Int
  RunningStats() = new(0, 0, 0)
end

numinfected(stats::RunningStats) = stats.num_infected
numdetected(stats::RunningStats) = stats.num_detected
numdead(stats::RunningStats) = stats.num_dead

function  update!(stats::RunningStats, event::Event)
  ek = kind(event)
  if isdetection(ek)
    stats.num_detected += 1
  elseif istransmission(ek)
    stats.num_infected += 1
  elseif isdeath(ek)
    stats.num_dead += 1
  end
  nothing
end


text:
matrixregression_insts = [
    vcat((4, 3), [(n, m) for n in ns]) # includes compile run
    for (m, ns) in ((15, vcat(50:50:150, 250, 500:500:3000)),
        (30, vcat(50:50:150, 250:250:1500)))
    ]

insts = OrderedDict()
insts["nat"] = (nothing, matrixregression_insts)
insts["ext"] = (:SOCExpPSD, matrixregression_insts)
return (MatrixRegressionJuMP, insts)


text:
using Parameters
using ReverseDiff
using ForwardDiff

abstract type AbstractRatevec end
export AbstractRatevec

@with_kw struct Arrheniusvec{N<:Array,K<:Array,Q<:Array} <: AbstractRatevec
        A::N
        n::K
        Ea::Q
end
function Arrheniusvec(arrs::T) where {T<:AbstractArray}
    A = zeros(length(arrs))
    n = zeros(length(arrs))
    Ea = zeros(length(arrs))
    for (i,aval) in enumerate(arrs)
        A[i] = aval.A
        n[i] = aval.n
        Ea[i] = aval.Ea
    end
    return Arrheniusvec(A=A,n=n,Ea=Ea)
end
@inline (arr::Arrheniusvec)(;T::Q,P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
@inline (arr::Arrheniusvec)(T::Q;P::N=0.0,C::S=0.0,phi=0.0) where {Q<:Real,N<:Real,S<:Real} = @fastmath @inbounds arr.A.*exp.(arr.n.*log(T).-arr.Ea.*(1.0/(R*T)))
export Arrheniusvec

@with_kw struct Chebyshevvec{T<:AbstractArray,Q<:Real,S<:Real,V<:Real,B<:Real} <: AbstractRate
    coefs::T
    Tmin::Q
    Tmax::S
    Pmin::V
    Pmax::B
end

function Chebyshevvec(chevs::T) where {T<:AbstractArray}
    coefs = zeros(length(chevs),size(chevs[1].coefs)[1],size(chevs[1].coefs)[2])
    for k in 1:size(chevs[1].coefs)[2]
        for j in 1:size(chevs[1].coefs)[1]
            for i in 1:length(chevs)
                coefs[i,j,k] = chevs[i].coefs[j,k]
            end
        end
    end
    return Chebyshevvec(coefs=coefs,Tmin=chevs[1].Tmin,Tmax=chevs[1].Tmax,Pmin=chevs[1].Pmin,Pmax=chevs[1].Pmax)
end
export Chebyshevvec

@inline function evalChebyshevPolynomial(ch::Chebyshevvec,n::N,x::Q) where {N<:Integer,Q<:Real}
    """
    evaluate the nth order Chebyshev Polynomial at x
    """
    if n==0
        return 1.0
    elseif n == 1
        return x
    else
        T = 0.0
        T0 = 1.0
        T1 = x
        for i in 1:(n-1)
            @fastmath T = 2*x*T1-T0
            T0 = T1
            T1 = T
        end
        return T
    end
end
export evalChebyshevPolynomial

@inline function getredtemp(ch::Chebyshevvec,T::N) where {N<:Real}
    """
    return a reduced temperature corresponding to the given tempeprature
    for the Chebyshev polynomial, maps the inverse of temperautre onto [-1,1]
    """
    return @fastmath (2.0/T-1.0/ch.Tmin-1.0/ch.Tmax)/(1.0/ch.Tmax-1.0/ch.Tmin)
end
export getredtemp

@inline function getredpress(ch::Chebyshevvec,P::N) where {N<:Real}
    """
    return a reduced pressure corresponding to the given temperature
    for the Chebyshev polynomial maps the logarithm of pressure onto [-1,1]
    """
    return @fastmath (2.0*log10(P)-log10(ch.Pmin)-log10(ch.Pmax))/(log10(ch.Pmax)-log10(ch.Pmin))
end
export getredpress

@inline function (ch::Chebyshevvec)(;T::N,P::Q=0.0,C::B=0.0,phi=0.0) where {N<:Real,B<:Real,Q<:Real}
    k = zeros(N,size(ch.coefs)[1])
    Tred = getredtemp(ch,T)
    Pred = getredpress(ch,P)
    klen,Tlen,Plen = size(ch.coefs)
    @simd for j = 1:Plen
        @simd for i = 1:Tlen
            @fastmath @inbounds k .+= ch.coefs[:,i,j].*(evalChebyshevPolynomial(ch,i-1,Tred)*evalChebyshevPolynomial(ch,j-1,Pred))
        end
    end
    return @fastmath 10.0.^k
end

@with_kw struct Troevec{P<:AbstractArray,Q<:AbstractArray,F<:AbstractArray,L<:AbstractArray,N<:Integer,K<:AbstractFloat,R<:AbstractRateUncertainty} <: AbstractFalloffRate
    arrhigh::Arrheniusvec
    arrlow::Arrheniusvec
    a::P
    T3::Q
    T1::F
    T2::L
    efficiencies::Array{Dict{N,K},1}
    unc::R = EmptyRateUncertainty()
end

function Troevec(troes::T) where {T<:AbstractArray}

    Ahigh = zeros(length(troes))
    nhigh = zeros(length(troes))
    Ehigh = zeros(length(troes))
    Alow = zeros(length(troes))
    nlow = zeros(length(troes))
    Elow = zeros(length(troes))
    a = zeros(length(troes))
    T3 = zeros(length(troes))
    T1 = zeros(length(troes))
    T2 = zeros(length(troes))
    efficiencies = [Dict{Int64,Float64}() for x in troes]
    for (i,falloff) in enumerate(troes)
        if isa(falloff, ThirdBody)
            Alow[i] = falloff.arr.A
            nlow[i] = falloff.arr.n
            Elow[i] = falloff.arr.Ea
            Ahigh[i] = 1e20  #very high limited by energy transfer process
            T3[i] = Inf
            T1[i] = Inf
            if length(efficiencies) > 0
                efficiencies[i] = falloff.efficiencies
            end
        elseif isa(falloff,Lindemann)
            Alow[i] = falloff.arrlow.A
            nlow[i] = falloff.arrlow.n
            Elow[i] = falloff.arrlow.Ea
            Ahigh[i] = falloff.arrhigh.A
            nhigh[i] = falloff.arrhigh.n
            Ehigh[i] = falloff.arrhigh.Ea
            T3[i] = Inf
            T1[i] = Inf
            if length(efficiencies) > 0
                efficiencies[i] = falloff.efficiencies
            end
        elseif isa(falloff,Troe)
            Alow[i] = falloff.arrlow.A
            nlow[i] = falloff.arrlow.n
            Elow[i] = falloff.arrlow.Ea
            Ahigh[i] = falloff.arrhigh.A
            nhigh[i] = falloff.arrhigh.n
            Ehigh[i] = falloff.arrhigh.Ea
            a[i] = falloff.a
            T3[i] = falloff.T3
            T1[i] = falloff.T1
            T2[i] = falloff.T2
            if length(efficiencies) > 0
                efficiencies[i] = falloff.efficiencies
            end
        else
            val = typeof(falloff)
            error("could not process type $val in Troevec creation")
        end
    end
    return Troevec(arrhigh=Arrheniusvec(Ahigh,nhigh,Ehigh),arrlow=Arrheniusvec(Alow,nlow,Elow),
        a=a,T3=T3,T1=T1,T2=T2,efficiencies=efficiencies)
end
export Troevec

@inline function (tr::Troevec)(;T::Q=nothing,P::R=0.0,C::S=nothing,phi=0.0) where {Q<:Real,R<:Real,S<:Real}
    k0 = tr.arrlow(T=T)
    kinf = tr.arrhigh(T=T)
    @fastmath Pr = k0.*C./kinf
    @fastmath log10Pr = log10.(Pr)
    @fastmath log10Fcent = log10.((1.0.-tr.a).*exp.(-T./tr.T3).+tr.a.*exp.(-T./tr.T1).+exp.(-tr.T2./T))
    d = 0.14
    @fastmath n = 0.75.-1.27.*log10Fcent
    @fastmath c = -0.4.-0.67.*log10Fcent
    F = 10.0.^((.!isinf.(tr.T1)) .* @fastmath (log10Fcent./(1.0.+((log10Pr.+c)./(n.-d.*(log10Pr))).^2)))
    return @fastmath ((k0.*C)./(1.0.+Pr)).*F
end

@with_kw struct PdepArrheniusvec{T<:Real,Q<:AbstractRateUncertainty,Z<:Arrheniusvec} <: AbstractRate
    Ps::Array{T,1}
    arrvecs::Array{Z,1}
    unc::Q = EmptyRateUncertainty()
end

function PdepArrheniusvec(pdeparrs::T) where {T<:AbstractArray}
    Ps = pdeparrs[1].Ps
    arrs = [Array{Arrhenius,1}() for i = 1:length(Ps)]
    for ind in 1:length(pdeparrs)
        for i =1:length(Ps)
            push!(arrs[i],pdeparrs[ind].arrs[i])
        end
    end
    return PdepArrheniusvec(;Ps=Ps,arrvecs=Arrheniusvec.(arrs))
end
export PdepArrheniusvec

@inline function (parr::PdepArrheniusvec)(;T::Q=nothing,P::V=nothing,C::S=0.0,phi=0.0) where {Q<:Real,V<:Real,S<:Real}
    inds = getBoundingIndsSorted(P,parr.Ps)::Tuple{Int64,Int64}
    if inds[2] == -1
        return @inbounds parr.arrvecs[inds[1]](T=T)
    else
        @inbounds highk = parr.arrvecs[inds[2]](T=T)
        @inbounds lowk = parr.arrvecs[inds[1]](T=T)
        @inbounds Plow = parr.Ps[inds[1]]
        @inbounds Phigh = parr.Ps[inds[2]]
        return @inbounds @fastmath lowk.*10.0.^(log10.(P./Plow)/log10.(Phigh./Plow)*log10.(highk./lowk))
    end
end
export PdepArrheniusvec


text:
export Extent, extent

"""
A pair of `Point`s representing a rectangle.
"""
struct Extent
    upperleft::Point
    lowerright::Point
end

"""
Compute the center of an `Extent` as a `Point`.
"""
function center(ext::Extent)
    (ext.upperleft + ext.lowerright) / 2
end

"""
Compute an extent which encompasses the given two `Extent`s.
"""
function combine_extents(ext1::Extent, ext2::Extent)
    Extent(min.(ext1.upperleft, ext2.upperleft), max.(ext1.lowerright, ext2.lowerright))
end

"""
Compute an extent which encompasses the given `Path`.
"""
function extent(path::Path)
    upperleft = lowerright = path[1]
    for point in path
        upperleft = min.(upperleft, point)
        lowerright = max.(lowerright, point)
    end
    Extent(upperleft, lowerright)
end

"""
Compute an extent which encompasses the given `MultiPath`.
"""
function extent(paths::MultiPath)
    reduce(combine_extents, [extent(path) for path in paths])
end

"""
Compute an extent which encompasses the given `Layer`.
"""
function extent(layer::Layer)
    extent(layer.paths)
end

"""
Compute an extent which encompasses the given `PenPlot`.
"""
function extent(plot::PenPlot)
    reduce(combine_extents, [extent(layer) for layer in plot.layers])
end


text:
import Boltzmann: split_evenly

if :unit in TEST_GROUPS
    @test split_evenly(10, 2) == [(1,2), (3,4), (5,6), (7,8), (9,10)]
    @test split_evenly(11, 2) == [(1,2), (3,4), (5,6), (7,8), (9,10), (11,11)]
    @test split_evenly(3, 5) == [(1,3)]
end

if :integration in TEST_GROUPS
    # Ensure that the test methods all work correctly
    df1 = Boltzmann.compare(GRBM(10, 5); n_obs=10)
    @test size(df1, 1) > 0

    # Try a comparison with 2 different contexts
    # in this case just use 2 different reporters
    context1 = Boltzmann.DefaultContext()
    context1[:reporter].log = true
    context2 = Boltzmann.DefaultContext()
    context2[:reporter] = Boltzmann.TextReporter()
    df2 = Boltzmann.compare(GRBM(10, 5), context1, context2; input_size=10, n_obs=10)
    @test size(df2, 1) > 0

    # Test the benchmark function (since we won't be running
    # it under coverage normally.
    df3 = Boltzmann.benchmark(GRBM(10, 5), input_size=10, n_obs=10, debug=true)
    @test size(df3, 1) > 0
end


text:
module GalerkinSparseGrids

using HCubature
using JLD2
if haskey(ENV, "MKLROOT")
    using MKLSparse
end
using ODE
using Printf
using SparseArrays
using StaticArrays
using Plots

# The following script files are used:

include("schemes.jl")                    # CartesianIndex Manipulation Methods
include("additional_tools.jl")            # Helper functions not specific to any part of the package
include("dg_basis.jl")                    # Gram-Schmidt procedure for DG basis functions in 1-D
include("1d_dg_functions.jl")             # Explicitly building the 1-D Galerkin Basis
include("dg_methods.jl")                # Multidimensional hierarchical & sparse coefficients
include("dg_vmethods.jl")                # Going between a dictionary & a vector of coeffs
include("derivative_matrix_elements.jl")# 1-D symbolic piecewise derivative
include("1d_derivative.jl")                # Constructing ideal 1D derivative matrix using boundary terms
include("multidim_derivative.jl")        # Multidimensional DG Derivatives in full & sparse bases
include("pdes.jl")                        # Solving the n-D wave equation with periodic boundary
include("1d_nodal_basis.jl")            # Construct 1D nodal basis for multiplication
include("multidim_nodal_basis.jl")        # Construct Multidimensional nodal basis for multiplication

include("basic_function_exact_coeffs.jl") # Exact coefficients for certain easy functions
include("error_measure.jl")                # Monte Carlo Methods to measure error
include("tensor_construct.jl")            # Quickly calculates coeffs of simple tensors of functions
include("traveling_wave.jl")              # Constructs initial data for traveling wave

export

coeffs_DG,
reconstruct_DG,

get_size,
V2D,
D2V,
V2Dref,
D2Vref,
vcoeffs_DG,

transform,
make_modal2point_matrices,
make_point2modal_matrices,

D_matrix,
grad_matrix,
laplacian_matrix,

nodal2points_1D,
points2nodal_1D,
nodal2heir_1D,
hier2points_1D,
points2hier_1D,

get_one_modal,
get_xi_modal,
get_r2_modal,
get_xi_point,
get_r2_point,

mcerr,
mcerr2,

tensor_construct,

wave_evolve_1D,
wave_evolve,
vlasov_evolve,
traveling_wave,

periodic_pos_DLF_matrix,
periodic_hier_DLF_matrix,
periodic_DLF_matrix,
hier2pos,
spzeros,
get_size,
cutoff

end # module


text:
################################################################################
# Types of Mixed Potentials
################################################################################

include("../mixed/linear_combination_potential.jl")


text:
struct GateKernel
	f::Function
	paramCount::Int
	name::String
	GateKernel(f::Function, name::String="Anonymous") = new(f, 0, name)
	GateKernel(f::Function, paramCount::Int,name::String="Anonymous") = new(f, paramCount, name)
end

func(kernel::GateKernel) = kernel.f
paramCount(kernel::GateKernel) = kernel.paramCount
name(kernel::GateKernel) = kernel.name

function control(kernel::GateKernel)
	function controlFunc(p::Real...)
		controlled = zeros(ComplexF64,2,2,2,2)
		controlled[1,:,1,:] = diagm(0=>ones(ComplexF64, 2))
		controlled[2,:,2,:] = func(kernel)(p...)
		return controlled
	end
	return GateKernel(controlFunc, paramCount(kernel), "Control-"*name(kernel))
end

function inverse(kernel::GateKernel)
	# TODO: need optimization
	function inverse(p::Real...)
		tensor = conj(kernel.f(p...))
		shape = size(tensor)
		n = Int(√length(tensor))
		tensor = transpose(reshape(tensor, n,n))
		return Array(reshape(tensor, shape...))
	end
	return GateKernel(inverse, paramCount(kernel), name(kernel)*"†")
end


text:
"""
Terra Object

- `class`      :
- `dimensions` :
- `resolution` :
- `extent`     :
- `coord.ref.` :
- `source`     :
- `names`      : `[bandnames]`
"""

# https://stackoverflow.com/questions/32609570/how-to-set-the-band-description-option-tag-of-a-geotiff-file-using-gdal-gdalw

gdal_close(ds::Ptr{Nothing}) = GDAL.gdalclose(ds)

# gdal_open(file::AbstractString) = ArchGDAL.read(file)
function gdal_open(f::AbstractString, mode = GDAL.GA_ReadOnly, args...)
  GDAL.gdalopen(f, mode, args...)
end

function gdal_open(func::Function, args...; kwargs...)
  ds = gdal_open(args...; kwargs...)
  try
    func(ds)
  finally
    gdal_close(ds)
  end
end

function nband(f)
  gdal_open(f) do ds
    GDAL.gdalgetrastercount(ds)
  end
end
# function nband(file::AbstractString)
#     # ArchGDAL.unsafe_read(file) do ds
#     ArchGDAL.read(file) do ds
#         ArchGDAL.nraster(ds)
#     end
# end
nratser = nband
nlyr = nband

function bandnames(f)
  n = nband(f)
  gdal_open(f) do ds

    map(iband -> begin
      band = GDAL.gdalgetrasterband(ds, iband)
      GDAL.gdalgetdescription(band)
    end, 1:n)
  end
end

# works
function set_bandnames(f, bandnames)
  n = nband(f)
  gdal_open(f, GDAL.GA_Update) do ds
    map(iband -> begin
        band = GDAL.gdalgetrasterband(ds, iband)
        GDAL.gdalsetdescription(band, bandnames[iband])
      end, 1:n)
  end
  nothing
end

function gdalinfo(file::AbstractString)
    ds = ArchGDAL.read(file)
    gt = ArchGDAL.getgeotransform(ds)
    # band = ArchGDAL.getband(ds, 1)
    w, h = ArchGDAL.width(ds), ArchGDAL.height(ds)
    dx, dy = gt[2], -gt[end]
    x0 = gt[1] #+ dx/2
    x1 = x0 + w* dx
    y1 = gt[4] #- dy/2
    y0 = y1 - h*dy
    b = bbox(x0, y0, x1, y1)

    lon = x0 + dx/2 : dx: x1
    lat = reverse(y0 + dy/2 : dy: y1)
    nband = ArchGDAL.nraster(ds)

    Dict(
        # "file"     => basename(file),
        "bbox"     => b,
        "cellsize" => [dx, dy],
        "coords"   => [lon, lat],
        "dim"      => (Int64.([w, h, nband])..., ) # convert to tuple
    )
end

"""
    gdalinfo(ga::GeoArray)
    gdalinfo(file::AbstractString)

get detailed GDAL information

## return
- `range`    : [lon_min, lon_max, lat_min, lat_max]
- `cellsize` : [dx, dy]
- `coords`   : [lon, lat]
- `dim`      : (width, height, ntime)
"""
function gdalinfo(ga::GeoArray)
    # mid = [1, 1];
    # if length(mid) == 1; mid = [mid, mid]; end
    dx = ga.f.linear[1]
    dy = abs(ga.f.linear[4])
    dy2 = ga.f.linear[4]

    mid = [1, 1] # default mid is true
    delta = [dx, dy]/2 .* mid

    b = st_bbox(ga)
    lon = b.xmin + delta[1]:dx:b.xmax
    lat = b.ymin + delta[2]:dy:b.ymax
    if dy2 < 0; lat = reverse(lat); end

    Dict(
        "bbox"     => b,
        "cellsize" => [dx, dy],
        "coords"   => [lon, lat],
        "dim"      => size(ga.A))
end

export gdal_open, gdal_close, nband, bandnames, set_bandnames
export gdalinfo


text:
"Holds the tableau of an Specialised Partitioned Additive Runge-Kutta method for Variational systems."
struct TableauSLRK{DT <: Number} <: AbstractTableau{DT}
    name::Symbol
    o::Int
    s::Int

    q::CoefficientsRK{DT}
    p::CoefficientsRK{DT}

    q̃::CoefficientsRK{DT}
    p̃::CoefficientsRK{DT}

    ω::Matrix{DT}
    d::Vector{DT}

    function TableauSLRK(name::Symbol, o::Int, s::Int,
                        q::CoefficientsRK{DT}, p::CoefficientsRK{DT},
                        q̃::CoefficientsRK{DT}, p̃::CoefficientsRK{DT},
                        ω::Matrix{DT}, d=DT[]) where {DT}

        @assert s > 0 "Number of stages s must be > 0"

        @assert s==q.s==p.s==q̃.s==p̃.s
        @assert size(ω,1)==s
        @assert size(ω,2)==s+1

        @assert length(d)==0 || length(d)==s

        new{DT}(name, o, s, q, p, q̃, p̃, ω, d)
    end
end

"Parameters for right-hand side function of Specialised Partitioned Additive Runge-Kutta methods for Variational systems."
mutable struct ParametersSLRK{DT,TT,D,S,θT,FT,GT,G̅T,ϕT,ψT,tabType} <: Parameters{DT,TT}
    f_ϑ::θT
    f_f::FT
    f_g::GT
    f_g̅::G̅T
    f_ϕ::ϕT
    f_ψ::ψT

    Δt::TT

    tab::tabType

    t::TT
    q::Vector{DT}
    p::Vector{DT}
    λ::Vector{DT}

    function ParametersSLRK{DT,D,S}(f_ϑ::θT, f_f::FT, f_g::GT, f_g̅::G̅T, f_ϕ::ϕT, f_ψ::ψT, Δt::TT, tab::tabType) where {DT,TT,D,S,θT,FT,GT,G̅T,ϕT,ψT,tabType}
        # create solution vectors
        q = zeros(DT,D)
        p = zeros(DT,D)
        λ = zeros(DT,D)

        new{DT,TT,D,S,θT,FT,GT,G̅T,ϕT,ψT,tabType}(f_ϑ, f_f, f_g, f_g̅, f_ϕ, f_ψ, Δt, tab, zero(TT), q, p, λ)
    end
end

function update_params!(params::ParametersSLRK, sol::AtomicSolutionPDAE)
    # set time for nonlinear solver and copy previous solution
    params.t  = sol.t
    params.q .= sol.q
    params.p .= sol.p
    params.λ .= sol.λ
end

@doc raw"""
Specialised Partitioned Additive Runge-Kutta integrator for degenerate
variational systems with projection on secondary constraint.

This integrator solves the following system of equations for the internal stages,
```math
\begin{align}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a^1_{ij} V_{n,j} + h \sum \limits_{j=1}^{\sigma} a^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., s , \\
P_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{a}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{a}^2_{ij} F^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{a}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., s , \\
0 &= \Phi_{n,i} , & i &= 1, ..., s ,
\end{align}
```
the projective stages
```math
\begin{align}
\tilde{Q}_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} \alpha^1_{ij} \tilde{V}_{n,j} + h \sum \limits_{j=1}^{\sigma} \alpha^2_{ij} \tilde{\Lambda}_{n,j} , & i &= 1, ..., \sigma , \\
\tilde{P}_{n,i} &= p_{n} + h \sum \limits_{j=1}^{s} \bar{\alpha}^1_{ij} F^1_{n,j} + h \sum \limits_{j=1}^{s} \bar{\alpha}^2_{ij} \tilde{F}^2_{n,j} + h \sum \limits_{j=1}^{\sigma} \bar{\alpha}^3_{ij} \tilde{F}^3_{n,j} , & i &= 1, ..., \sigma , \\
0 &= \tilde{\Phi}_{n,i} , & i &= 1, ..., \sigma , \\
0 &= \sum \limits_{j=1}^{\sigma} \omega_{ij} \tilde{\Psi}_{n,i} , & i &= 1, ..., \sigma-1 ,
\end{align}
```
and update rule
```math
\begin{align}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b^1_{i} V_{n,i} + h \sum \limits_{i=1}^{\sigma} b^2_{i} \tilde{\Lambda}_{n,i} , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} b^1_{i} F^1_{n,i} + h \sum \limits_{i=1}^{s} b^2_{i} F^2_{n,i} + h \sum \limits_{i=1}^{\sigma} b^3_{i} \tilde{F}^3_{n,i} , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{align}
```
with definitions
```math
\begin{align}
F^1_{n,i} &= \nabla H (Q_{n,i}) , & i &= 1, ..., s , \\
F^2_{n,i} &= \nabla \vartheta (Q_{n,i}) \cdot V_{n,i} , & i &= 1, ..., s , \\
\tilde{F}^3_{n,i} &= \nabla \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) \cdot \tilde{\Lambda}_{n,i} , & i &= 1, ..., \sigma , \\
\Phi_{n,i} &= \phi (Q_{n,i}, P_{n,i}) = P_{n,i} - \vartheta (Q_{n,i}) , & i &= 1, ..., s , \\
\tilde{\Phi}_{n,i} &= \phi (\tilde{Q}_{n,i}, \tilde{P}_{n,i}) = \tilde{P}_{n,i} - \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma , \\
\tilde{\Psi}_{n,i} &= \psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) = \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \vartheta (\tilde{Q}_{n,i}) , & i &= 1, ..., \sigma ,
\end{align}
```
so that
```math
\begin{align}
F^1_{n,i} + F^2_{n,i} &= \frac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\phi (Q_{n,i}, P_{n,i}) &= P_{n,i} - \frac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , & i &= 1, ..., s , \\
\psi (\tilde{Q}_{n,i}, \tilde{V}_{n,i}, \tilde{P}_{n,i}, \tilde{F}_{n,i}) &= \tilde{F}_{n,i} - \tilde{V}_{n,i} \cdot \nabla \frac{\partial L}{\partial v} (\tilde{Q}_{n,i}, \tilde{V}_{n,i}) , & i &= 1, ..., \sigma .
\end{align}
```
"""
struct IntegratorSLRK{DT, TT, ET <: VDAE{DT,TT},
                              PT <: ParametersSLRK{DT,TT},
                              ST <: NonlinearSolver{DT},
                              IT <: InitialGuessPODE{DT,TT}, D, S} <: AbstractIntegratorVSPARK{DT, TT}
    equation::ET
    tableau::TableauSLRK{TT}

    params::PT
    solver::ST
    iguess::IT
    cache::IntegratorCacheSPARK{DT,TT,D,S,S}
end

function IntegratorSLRK(equation::VDAE{DT,TT},
                        tableau::TableauSLRK{TT}, Δt::TT) where {DT,TT}
    D = equation.d
    S = tableau.s

    # println(tableau)

    N = 4*D*S

    if isdefined(tableau, :d) && length(tableau.d) > 0
        N += D
    end

    # create params
    params = ParametersSLRK{DT,D,S}(equation.ϑ, equation.f, equation.g, equation.g̅, equation.ϕ, equation.ψ, Δt, tableau)

    # create solver
    solver = create_nonlinear_solver(DT, N, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache
    cache = IntegratorCacheSPARK{DT,TT,D,S,S}()

    # create integrator
    IntegratorSLRK{DT, TT, typeof(equation), typeof(params), typeof(solver), typeof(iguess), D, S}(
                              equation, tableau, params, solver, iguess, cache)
end

pstages(int::IntegratorSLRK) = int.tableau.s

function compute_stages!(x::Vector{ST}, cache::IntegratorCacheSPARK{ST,TT,D,S},
                                        params::ParametersSLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    local t::TT

    for i in 1:S
        for k in 1:D
            # copy y to Y, Z and Λ
            cache.Yp[i][k] = x[4*(D*(i-1)+k-1)+1]
            cache.Zp[i][k] = x[4*(D*(i-1)+k-1)+2]
            cache.Vp[i][k] = x[4*(D*(i-1)+k-1)+3]
            cache.Λp[i][k] = x[4*(D*(i-1)+k-1)+4]

            # compute Q and P
            cache.Qp[i][k] = params.q[k] + params.Δt * cache.Yp[i][k]
            cache.Pp[i][k] = params.p[k] + params.Δt * cache.Zp[i][k]
        end

        # compute f(X)
        t = params.t + params.Δt * params.tab.p.c[i]
        params.f_f(t, cache.Qp[i], cache.Vp[i], cache.Fp[i])
        params.f_g(t, cache.Qp[i], cache.Λp[i], cache.Gp[i])
        params.f_ϕ(t, cache.Qp[i], cache.Pp[i], cache.Φp[i])
        params.f_ψ(t, cache.Qp[i], cache.Pp[i], cache.Vp[i], cache.Fp[i], cache.Ψp[i])
    end

    if length(params.tab.d) > 0
        for k in 1:D
            cache.μ[k] = x[4*D*S+k]
        end
    end

    # compute q and p
    cache.q̃ .= params.q
    cache.p̃ .= params.p
    for i in 1:S
        cache.q̃ .+= params.Δt .* params.tab.q.b[i] .* cache.Vp[i]
        cache.q̃ .+= params.Δt .* params.tab.q̃.b[i] .* cache.Λp[i]
        cache.p̃ .+= params.Δt .* params.tab.p.b[i] .* cache.Fp[i]
        cache.p̃ .+= params.Δt .* params.tab.p̃.b[i] .* cache.Gp[i]
    end

    # compute ϕ(q,p)
    t = params.t + params.Δt
    params.f_ϕ(t, cache.q̃, cache.p̃, cache.ϕ̃)
end

"Compute stages of specialised partitioned additive Runge-Kutta methods for variational systems."
@generated function function_stages!(y::Vector{ST}, b::Vector{ST}, params::ParametersSLRK{DT,TT,D,S}) where {ST,DT,TT,D,S}
    cache = IntegratorCacheSPARK{ST,TT,D,S,S}()

    quote
        compute_stages!(y, $cache, params)

        # compute b = - [(Y-AV-AU), (Z-AF-AG), Φ, ωΨ]
        for i in 1:S
            for k in 1:D
                b[4*(D*(i-1)+k-1)+1] = - $cache.Yp[i][k]
                b[4*(D*(i-1)+k-1)+2] = - $cache.Zp[i][k]
                b[4*(D*(i-1)+k-1)+3] = - $cache.Φp[i][k]
                b[4*(D*(i-1)+k-1)+4] = params.tab.ω[i,S+1] * $cache.ϕ̃[k]
                for j in 1:S
                    b[4*(D*(i-1)+k-1)+1] += params.tab.q.a[i,j] * $cache.Vp[j][k]
                    b[4*(D*(i-1)+k-1)+1] += params.tab.q̃.a[i,j] * $cache.Λp[j][k]
                    b[4*(D*(i-1)+k-1)+2] += params.tab.p.a[i,j] * $cache.Fp[j][k]
                    b[4*(D*(i-1)+k-1)+2] += params.tab.p̃.a[i,j] * $cache.Gp[j][k]
                    b[4*(D*(i-1)+k-1)+4] += params.tab.ω[i,j]   * $cache.Ψp[j][k]
                end
            end
        end

        if length(params.tab.d) > 0
            for i in 1:S
                for k in 1:D
                    b[4*(D*(i-1)+k-1)+2] += $cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
                    b[4*(D*(i-1)+k-1)+3] += $cache.μ[k] * params.tab.d[i] / params.tab.p.b[i]
                end
            end

            for k in 1:D
                b[4*D*S+k] = 0
                for i in 1:S
                    b[4*D*S+k] -= $cache.Vp[i][k] * params.tab.d[i]
                end
            end
        end
    end
end

function initial_guess!(int::IntegratorSLRK, sol::AtomicSolutionPDAE)
    for i in 1:nstages(int)
        evaluate!(int.iguess, sol.q, sol.p, sol.v, sol.f,
                              sol.q̅, sol.p̅, sol.v̅, sol.f̅,
                              int.cache.q̃, int.cache.p̃, int.cache.ṽ, int.cache.f̃,
                              tableau(int).q̃.c[i], tableau(int).p̃.c[i])

        for k in 1:ndims(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+1] = (int.cache.q̃[k] - sol.q[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+2] = (int.cache.p̃[k] - sol.p[k])/timestep(int)
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+3] = int.cache.ṽ[k]
            int.solver.x[4*(ndims(int)*(i-1)+k-1)+4] = 0
        end
    end

    if isdefined(tableau(int), :d) && length(tableau(int).d) > 0
        for k in 1:ndims(int)
            int.solver.x[4*ndims(int)*pstages(int)+k] = 0
        end
    end
end

function update_solution!(int::IntegratorSLRK{DT,TT}, sol::AtomicSolutionPDAE{DT,TT}) where {DT,TT}
    # compute final update
    update_solution!(sol.q, int.cache.Vp, int.params.tab.q.b, timestep(int))
    update_solution!(sol.q, int.cache.Λp, int.params.tab.q̃.b, timestep(int))
    update_solution!(sol.p, int.cache.Fp, int.params.tab.p.b, timestep(int))
    update_solution!(sol.p, int.cache.Gp, int.params.tab.p̃.b, timestep(int))
end


text:
using Replay

lemma = joinpath(@__DIR__, "lemma.jl")

instructions = [
    "using UnicodePlots",
    raw"""
    function clearline(; move_up::Bool = false)
        buf = IOBuffer()
        print(buf, "\x1b[2K") # clear line
        print(buf, "\x1b[999D") # rollback the cursor
        move_up && print(buf, "\x1b[1A") # move up
        print(buf |> take! |> String)
    end
    """,
    raw"""
    function clearlines(H::Integer)
        for i = 1:H
            clearline(move_up = true)
        end
    end
    """,
    raw"""
    function animate(frames::Vector{<:UnicodePlots.Plot}; duration::Float64 = 0.5)
        print("\x1b[?25l") # hidecursor
        for (n, f) in enumerate(frames)
            print(f)
            str = string(f)
            if n != length(frames)
                sleep(duration)
                clearlines(length(collect(eachmatch(r"\n", str))))
            end
        end
        print("\u001B[?25h") # unhide cursor
    end
    """,
    "frames = [lineplot(x -> sin(x - t), width = 40) for t = -2π:0.25π:2π];",
    "animate(frames)",
]

replay(instructions; julia_project=@__DIR__, use_ghostwriter=false)


text:
import DrWatson: quickactivate
quickactivate(@__DIR__, "Chemostat_Heerden2013")

## -------------------------------------------------------------------
# Run add https://github.com/josePereiro/Chemostat_Heerden2013.jl in the Julia Pkg REPL to install the
# package, then you must activate the package enviroment (see README)
import Chemostat_Heerden2013
const Hd  = Chemostat_Heerden2013.HeerdenData;
const ECC = Chemostat_Heerden2013.EColiCore

## -------------------------------------------------------------------
# run add "https://github.com/josePereiro/Chemostat" in the
# julia Pkg REPL for installing the package
import Chemostat
const ChU = Chemostat.Utils
const ChSS = Chemostat.SteadyState
const ChLP = Chemostat.LP
const ChEP = Chemostat.MaxEntEP
const ChSU = Chemostat.SimulationUtils

## -------------------------------------------------------------------
# Simulation
DATA = Dict()
let
    model = ChU.MetNet(ChU.load_data(ECC.BASE_MODELS_FILE))
    objidx = ChU.rxnindex(model, ECC.BIOMASS_IDER)
    M, N = size(model)

    for (exp, cGLC) in enumerate(Hd.val("cGLC"))

        @info "Doing" exp cGLC
        D = get!(DATA, exp, Dict())

        # prepare model
        exp_growth = Hd.val("D", exp)
        expξ = Hd.val("xi", exp)
        intake_info = deepcopy(ECC.base_intake_info)
        intake_info["EX_glc__D_e"]["c"] = cGLC
        ChSS.apply_bound!(model, expξ, intake_info)
        # ChU.bounds!(model, ECC.BIOMASS_IDER, 0.0, exp_growth * 1.1)

        # maxent ep params
        alpha = Inf
        beta_vec = zeros(size(model, 2))
        epsconv = 1e-7
        alpha = Inf
        verbose = false

        # gradien descent
        epouts = Dict()
        x0 = [0.0]
        x1 = [10.0]
        C = [5e2]
        th = 1e-4
        epout_seed = nothing
        target = [Hd.val("D", exp)]
        expβ = ChSU.grad_desc(;x0, x1, th, C, target, maxiters = 1500) do betas
            beta = first(betas)
            if haskey(epouts, beta)
                epout = epouts[beta]
            else
                beta_vec[objidx] = beta
                solution = epout_seed
                epout = ChEP.maxent_ep(model; alpha, beta_vec, epsconv, verbose, solution)
                epouts[beta] = epout
            end
            epout_seed = epout
            r = [ChU.av(epout)[objidx]]
        end

        # collect results
        D["model"] = model
        D["exp_beta"] = first(expβ)
        D["exp_xi"] = expξ
        D["epouts"] = epouts
        D["fbaout"] = ChLP.fba(model, ECC.BIOMASS_IDER)

        println()
    end
end

## -------------------------------------------------------------------
# save results
ChU.save_data(ECC.MAXENT_RES_FILE, DATA)