text:
# Translation from our variables to Nocedal and Wright's
# JMW's dx <=> NW's s
# JMW's dg <=> NW' y

struct BFGS{IL, L, H, T, TM} <: FirstOrderOptimizer
    alphaguess!::IL
    linesearch!::L
    initial_invH::H
    initial_stepnorm::T
    manifold::TM
end

Base.summary(::BFGS) = "BFGS"

"""
# BFGS
## Constructor
```julia
BFGS(; alphaguess = LineSearches.InitialStatic(),
       linesearch = LineSearches.HagerZhang(),
       initial_invH = x -> Matrix{eltype(x)}(I, length(x), length(x)),
       manifold = Flat())
```

## Description
The `BFGS` method implements the Broyden-Fletcher-Goldfarb-Shanno algorithm as
described in Nocedal and Wright (sec. 8.1, 1999) and the four individual papers
Broyden (1970), Fletcher (1970), Goldfarb (1970), and Shanno (1970). It is a
quasi-Newton method that updates an approximation to the Hessian using past
approximations as well as the gradient. See also the limited memory variant
`LBFGS` for an algorithm that is more suitable for high dimensional problems.

## References
 - Wright, S. J. and J. Nocedal (1999), Numerical optimization. Springer Science 35.67-68: 7.
 - Broyden, C. G. (1970), The convergence of a class of double-rank minimization algorithms, Journal of the Institute of Mathematics and Its Applications, 6: 76–90.
 - Fletcher, R. (1970), A New Approach to Variable Metric Algorithms, Computer Journal, 13 (3): 317–322,
 - Goldfarb, D. (1970), A Family of Variable Metric Updates Derived by Variational Means, Mathematics of Computation, 24 (109): 23–26,
 - Shanno, D. F. (1970), Conditioning of quasi-Newton methods for function minimization, Mathematics of Computation, 24 (111): 647–656.
"""
function BFGS(; alphaguess = LineSearches.InitialStatic(), # TODO: benchmark defaults
                linesearch = LineSearches.HagerZhang(),  # TODO: benchmark defaults
                initial_invH = nothing,
                initial_stepnorm = nothing,
                manifold::Manifold=Flat())
    BFGS(alphaguess, linesearch, initial_invH, initial_stepnorm, manifold)
end

mutable struct BFGSState{Tx, Tm, T,G} <: AbstractOptimizerState
    x::Tx
    x_previous::Tx
    g_previous::G
    f_x_previous::T
    dx::Tx
    dg::Tx
    u::Tx
    invH::Tm
    s::Tx
    @add_linesearch_fields()
end

function initial_state(method::BFGS, options, d, initial_x::AbstractArray{T}) where T
    n = length(initial_x)
    initial_x = copy(initial_x)
    retract!(method.manifold, initial_x)

    value_gradient!!(d, initial_x)

    project_tangent!(method.manifold, gradient(d), initial_x)

    if method.initial_invH == nothing
        if method.initial_stepnorm == nothing
            invH0 = Matrix{T}(I, n, n)
        else
            initial_scale = method.initial_stepnorm * inv(norm(gradient(d), Inf))
            invH0 = Matrix{T}(initial_scale*I, n, n)
        end
    else
        invH0 = method.initial_invH(initial_x)
    end
    # Maintain a cache for line search results
    # Trace the history of states visited
    BFGSState(initial_x, # Maintain current state in state.x
              copy(initial_x), # Maintain previous state in state.x_previous
              copy(gradient(d)), # Store previous gradient in state.g_previous
              real(T)(NaN), # Store previous f in state.f_x_previous
              similar(initial_x), # Store changes in position in state.dx
              similar(initial_x), # Store changes in gradient in state.dg
              similar(initial_x), # Buffer stored in state.u
              invH0, # Store current invH in state.invH
              similar(initial_x), # Store current search direction in state.s
              @initial_linesearch()...)
end

function update_state!(d, state::BFGSState, method::BFGS)
    n = length(state.x)
    T = eltype(state.s)
    # Set the search direction
    # Search direction is the negative gradient divided by the approximate Hessian
    mul!(vec(state.s), state.invH, vec(gradient(d)))
    rmul!(state.s, T(-1))
    project_tangent!(method.manifold, state.s, state.x)

    # Maintain a record of the previous gradient
    copyto!(state.g_previous, gradient(d))

    # Determine the distance of movement along the search line
    # This call resets invH to initial_invH is the former in not positive
    # semi-definite
    lssuccess = perform_linesearch!(state, method, ManifoldObjective(method.manifold, d))

    # Update current position
    state.dx .= state.alpha.*state.s
    state.x .= state.x .+ state.dx
    retract!(method.manifold, state.x)

    lssuccess == false # break on linesearch error
end

function update_h!(d, state, method::BFGS)
    n = length(state.x)
    # Measure the change in the gradient
    state.dg .= gradient(d) .- state.g_previous

    # Update the inverse Hessian approximation using Sherman-Morrison
    dx_dg = real(dot(state.dx, state.dg))
    if dx_dg == 0.0
        return true # force stop
    end
    mul!(vec(state.u), state.invH, vec(state.dg))

    c1 = (dx_dg + real(dot(state.dg, state.u))) / (dx_dg' * dx_dg)
    c2 = 1 / dx_dg

    # TODO BLASify this
    # invH = invH + c1 * (s * s') - c2 * (u * s' + s * u')
    for i in 1:n
        @simd for j in 1:n
            @inbounds state.invH[i, j] += c1 * state.dx[i] * state.dx[j]' - c2 * (state.u[i] * state.dx[j]' + state.u[j]' * state.dx[i])
        end
    end
end

function trace!(tr, d, state, iteration, method::BFGS, options, curr_time=time())
    dt = Dict()
    dt["time"] = curr_time
    if options.extended_trace
        dt["x"] = copy(state.x)
        dt["g(x)"] = copy(gradient(d))
        dt["~inv(H)"] = copy(state.invH)
        dt["Current step size"] = state.alpha
    end
    g_norm = norm(gradient(d), Inf)
    update!(tr,
    iteration,
    value(d),
    g_norm,
    dt,
    options.store_trace,
    options.show_trace,
    options.show_every,
    options.callback)
end


text:
#=
Copyright (c) 2018-2022 Chris Coey, Lea Kapelevich, and contributors

This Julia package Hypatia.jl is released under the MIT license; see LICENSE
file in the root directory or at https://github.com/chriscoey/Hypatia.jl

common code for native examples
=#

abstract type ExampleInstanceNative{T <: Real} <: ExampleInstance{T} end

# fallback: just check optimal status
function test_extra(
    inst::ExampleInstanceNative,
    solve_stats::NamedTuple,
    ::NamedTuple,
    )
    @test solve_stats.status == Solvers.Optimal
end

function run_instance(
    ex_type::Type{<:ExampleInstanceNative{T}}, # an instance of a native example
    inst_data::Tuple,
    inst_options::NamedTuple = NamedTuple(),
    solver_type::Type{<:Solvers.Solver} = Solvers.Solver{T};
    default_options::NamedTuple = NamedTuple(),
    test::Bool = true,
    rseed::Int = 1,
    verbose::Bool = true,
    ) where {T <: Real}
    new_options = merge(default_options, inst_options)

    verbose && println("setup model")
    setup_time = @elapsed begin
        Random.seed!(rseed)
        inst = ex_type(inst_data...)
        model = build(inst)
        model_stats = get_model_stats(model)
        solver = Solvers.Solver{T}(; default_options..., inst_options...)
    end
    flush(stdout); flush(stderr)

    verbose && println("solve and check")
    check_time = @elapsed begin
        Solvers.load(solver, model)
        Solvers.solve(solver)
        flush(stdout); flush(stderr)

        (solve_stats, solution) = process_result(model, solver)

        test && test_extra(inst, solve_stats, solution)
    end
    flush(stdout); flush(stderr)

    return (; model_stats..., solve_stats..., setup_time,
        check_time, :script_status => "Success")
end

ext: jl


text:
import Base: show

function _show_indented(io::IO, dep::LibraryDependency, indent, lib)
    print_indented(io,"- Library \"$(dep.name)\"",indent+1)
    if !applicable(dep)
        println(io," (not applicable to this system)")
    else
        println(io)
        if !isempty(lib)
            print_indented(io,"- Satisfied by:\n",indent+4)
            for (k,v) in lib
                print_indented(io,"- $(k[1]) at $v\n",indent+6)
            end
        end
        if length(dep.providers) > 0
            print_indented(io,"- Providers:\n",indent+4)
            for (p,opts) in dep.providers
                show_indented(io,p,indent+6)
                if !can_provide(p,opts,dep)
                    print(io," (can't provide)")
                end
                println(io)
            end
        end
    end
end
show_indented(io::IO, dep::LibraryDependency, indent) = _show_indented(io,dep,indent, applicable(dep) ? _find_library(dep) : nothing)
show(io::IO, dep::LibraryDependency) = show_indented(io, dep, 0)

function show(io::IO, deps::LibraryGroup)
    print(io," - Library Group \"$(deps.name)\"")
    all = allf(deps)
    providers = satisfied_providers(deps,all)
    if !(isempty(providers))
        print(io," (satisfied by ",join(providers,", "),")")
    end
    if !applicable(deps)
        println(io," (not applicable to this system)")
    else
        println(io)
        for dep in deps.deps
            _show_indented(io,dep,4,haskey(all,dep)? all[dep] : nothing)
        end
    end
end

function debug_context(pkg::String)
    info("Reading build script...")
    dir = Pkg.dir(pkg)
    file = joinpath(dir,"deps/build.jl")
    context = BinDeps.PackageContext(false,dir,pkg,{})
    m = Module(:__anon__)
    body = Expr(:toplevel,:(ARGS=[$context]),:(include($file)))
    eval(m,body)
    context
end

function debug(io,pkg::String)
    context = debug_context(pkg)
    println(io,"The package declares $(length(context.deps)) dependencies.")
    for dep in context.deps
        show(io,dep)
    end
end
debug(pkg::String) = debug(STDOUT,pkg)


text:
"""
    sinkhorn_unbalanced(μ, ν, C, λ1::Real, λ2::Real, ε; kwargs...)

Compute the optimal transport plan for the unbalanced entropically regularized optimal
transport problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, entropic regularization parameter `ε`, and marginal relaxation
terms `λ1` and `λ2`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma)
+ \\lambda_1 \\operatorname{KL}(\\gamma 1 | \\mu)
+ \\lambda_2 \\operatorname{KL}(\\gamma^{\\mathsf{T}} 1 | \\nu),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term and ``\\operatorname{KL}`` is the Kullback-Leibler divergence.

The keyword arguments supported here are the same as those in the `sinkhorn_unbalanced`
for unbalanced optimal transport problems with general soft marginal constraints.
"""
function sinkhorn_unbalanced(
    μ, ν, C, λ1::Real, λ2::Real, ε; proxdiv_F1=nothing, proxdiv_F2=nothing, kwargs...
)
    if proxdiv_F1 !== nothing && proxdiv_F2 !== nothing
        Base.depwarn(
            "keyword arguments `proxdiv_F1` and `proxdiv_F2` are deprecated",
            :sinkhorn_unbalanced,
        )

        # have to wrap the "proxdiv" functions since the signature changed
        # ε was fixed in the function, so we ignore it
        proxdiv_F1_wrapper(s, p, _) = copyto!(s, proxdiv_F1(s, p))
        proxdiv_F2_wrapper(s, p, _) = copyto!(s, proxdiv_F2(s, p))

        return sinkhorn_unbalanced(
            μ, ν, C, proxdiv_F1_wrapper, proxdiv_F2_wrapper, ε; kwargs...
        )
    end

    # define "proxdiv" functions for the unbalanced OT problem
    proxdivF!(s, p, ε, λ) = (s .= (p ./ s) .^ (λ / (ε + λ)))
    proxdivF1!(s, p, ε) = proxdivF!(s, p, ε, λ1)
    proxdivF2!(s, p, ε) = proxdivF!(s, p, ε, λ2)

    return sinkhorn_unbalanced(μ, ν, C, proxdivF1!, proxdivF2!, ε; kwargs...)
end

"""
    sinkhorn_unbalanced(
        μ, ν, C, proxdivF1!, proxdivF2!, ε;
        atol=0, rtol=atol > 0 ? 0 : √eps, check_convergence=10, maxiter=1_000,
    )

Compute the optimal transport plan for the unbalanced entropically regularized optimal
transport problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, entropic regularization parameter `ε`, and soft marginal
constraints ``F_1`` and ``F_2`` with "proxdiv" functions `proxdivF!` and `proxdivG!`.

The optimal transport plan `γ` is of the same size as `C` and solves
```math
\\inf_{\\gamma} \\langle \\gamma, C \\rangle
+ \\varepsilon \\Omega(\\gamma)
+ F_1(\\gamma 1, \\mu)
+ F_2(\\gamma^{\\mathsf{T}} 1, \\nu),
```
where ``\\Omega(\\gamma) = \\sum_{i,j} \\gamma_{i,j} \\log \\gamma_{i,j}`` is the entropic
regularization term and ``F_1(\\cdot, \\mu)`` and ``F_2(\\cdot, \\nu)`` are soft marginal
constraints for the source and target marginals.

The functions `proxdivF1!(s, p, ε)` and `proxdivF2!(s, p, ε)` evaluate the "proxdiv"
functions of ``F_1(\\cdot, p)`` and ``F_2(\\cdot, p)`` at ``s`` for the entropic
regularization parameter ``\\varepsilon``. They have to be mutating and overwrite the first
argument `s` with the result of their computations.

Mathematically, the "proxdiv" functions are defined as
```math
\\operatorname{proxdiv}_{F_i}(s, p, \\varepsilon)
= \\operatorname{prox}^{\\operatorname{KL}}_{F_i(\\cdot, p)/\\varepsilon}(s) \\oslash s
```
where ``\\oslash`` denotes element-wise division and
``\\operatorname{prox}_{F_i(\\cdot, p)/\\varepsilon}^{\\operatorname{KL}}`` is the proximal
operator of ``F_i(\\cdot, p)/\\varepsilon`` for the Kullback-Leibler
(``\\operatorname{KL}``) divergence.  It is defined as
```math
\\operatorname{prox}_{F}^{\\operatorname{KL}}(x)
= \\operatorname{argmin}_{y} F(y) + \\operatorname{KL}(y|x)
```
and can be computed in closed-form for specific choices of ``F``. For instance, if
``F(\\cdot, p) = \\lambda \\operatorname{KL}(\\cdot | p)`` (``\\lambda > 0``), then
```math
\\operatorname{prox}_{F(\\cdot, p)/\\varepsilon}^{\\operatorname{KL}}(x)
= x^{\\frac{\\varepsilon}{\\varepsilon + \\lambda}} p^{\\frac{\\lambda}{\\varepsilon + \\lambda}},
```
where all operators are acting pointwise.[^CPSV18]

Every `check_convergence` steps it is assessed if the algorithm is converged by checking if
the iterates of the scaling factor in the current and previous iteration satisfy
`isapprox(vcat(a, b), vcat(aprev, bprev); atol=atol, rtol=rtol)` where `a` and `b` are the
current iterates and `aprev` and `bprev` the previous ones. The default `rtol` depends on
the types of `μ`, `ν`, and `C`. After `maxiter` iterations, the computation is stopped.

[^CPSV18]: Chizat, L., Peyré, G., Schmitzer, B., & Vialard, F.-X. (2018). [Scaling algorithms for unbalanced optimal transport problems](https://doi.org/10.1090/mcom/3303). Mathematics of Computation, 87(314), 2563–2609.

See also: [`sinkhorn_unbalanced2`](@ref)
"""
function sinkhorn_unbalanced(
    μ,
    ν,
    C,
    proxdivF1!,
    proxdivF2!,
    ε;
    tol=nothing,
    atol=tol,
    rtol=nothing,
    max_iter=nothing,
    maxiter=max_iter,
    check_convergence::Int=10,
)
    # deprecations
    if tol !== nothing
        Base.depwarn(
            "keyword argument `tol` is deprecated, please use `atol` and `rtol`",
            :sinkhorn_unbalanced,
        )
    end
    if max_iter !== nothing
        Base.depwarn(
            "keyword argument `max_iter` is deprecated, please use `maxiter`",
            :sinkhorn_unbalanced,
        )
    end

    # compute Gibbs kernel
    K = @. exp(-C / ε)

    # set default values of squared tolerances
    T = float(Base.promote_eltype(μ, ν, K))
    sqatol = atol === nothing ? 0 : atol^2
    sqrtol = rtol === nothing ? (sqatol > zero(sqatol) ? zero(T) : eps(T)) : rtol^2

    # initialize iterates
    a = similar(μ, T)
    sum!(a, K)
    proxdivF1!(a, μ, ε)
    b = similar(ν, T)
    mul!(b, K', a)
    proxdivF2!(b, ν, ε)

    # caches for convergence checks
    a_old = similar(a)
    b_old = similar(b)

    isconverged = false
    _maxiter = maxiter === nothing ? 1_000 : maxiter
    for iter in 1:_maxiter
        # update cache if necessary
        ischeck = iter % check_convergence == 0
        if ischeck
            copyto!(a_old, a)
            copyto!(b_old, b)
        end

        # compute next iterates
        mul!(a, K, b)
        proxdivF1!(a, μ, ε)
        mul!(b, K', a)
        proxdivF2!(b, ν, ε)

        # check convergence of the scaling factors
        if ischeck
            # compute norm of current and previous scaling factors and their difference
            sqnorm_a_b = sum(abs2, a) + sum(abs2, b)
            sqnorm_a_b_old = sum(abs2, a_old) + sum(abs2, b_old)
            a_old .-= a
            b_old .-= b
            sqeuclidean_a_b = sum(abs2, a_old) + sum(abs2, b_old)
            @debug "Sinkhorn algorithm (" *
                   string(iter) *
                   "/" *
                   string(_maxiter) *
                   ": squared Euclidean distance of iterates = " *
                   string(sqeuclidean_a_b)

            # check convergence of `a`
            if sqeuclidean_a_b < max(sqatol, sqrtol * max(sqnorm_a_b, sqnorm_a_b_old))
                @debug "Sinkhorn algorithm ($iter/$_maxiter): converged"
                isconverged = true
                break
            end
        end
    end

    if !isconverged
        @warn "Sinkhorn algorithm ($_maxiter/$_maxiter): not converged"
    end

    return K .* a .* b'
end

"""
    sinkhorn_unbalanced2(μ, ν, C, λ1, λ2, ε; plan=nothing, kwargs...)
    sinkhorn_unbalanced2(μ, ν, C, proxdivF1!, proxdivF2!, ε; plan=nothing, kwargs...)

Compute the optimal transport plan for the unbalanced entropically regularized optimal
transport problem with source and target marginals `μ` and `ν`, cost matrix `C` of size
`(length(μ), length(ν))`, entropic regularization parameter `ε`, and marginal relaxation
terms `λ1` and `λ2` or soft marginal constraints with "proxdiv" functions `proxdivF1!` and
`proxdivF2!`.

A pre-computed optimal transport `plan` may be provided. The other keyword arguments
supported here are the same as those in the [`sinkhorn_unbalanced`](@ref) for unbalanced
optimal transport problems with general soft marginal constraints.

See also: [`sinkhorn_unbalanced`](@ref)
"""
function sinkhorn_unbalanced2(
    μ, ν, C, λ1_or_proxdivF1, λ2_or_proxdivF2, ε; plan=nothing, kwargs...
)
    γ = if plan === nothing
        sinkhorn_unbalanced(μ, ν, C, λ1_or_proxdivF1, λ2_or_proxdivF2, ε; kwargs...)
    else
        # check dimensions
        size(C) == (length(μ), length(ν)) ||
            error("cost matrix `C` must be of size `(length(μ), length(ν))`")
        size(plan) == size(C) || error(
            "optimal transport plan `plan` and cost matrix `C` must be of the same size",
        )
        plan
    end
    return dot(γ, C)
end


text:
"""
    ExplicitRungeKuttaCache <: AbstractRungeKuttaCache

A composite type for the [`AbstractRungeKuttaCache`](@ref) of an [`ExplicitRungeKuttaSolver`](@ref).

# Constructors
```julia
ExplicitRungeKuttaCache(n, m, v, k)
ExplicitRungeKuttaCache(problem, solver)
```

# Arguments
- `n       :: Integer`                                               : step counter.
- `m       :: Integer`                                               : adaptive correction counter.
- `v       :: AbstractVector{Union{Number, AbstractVector{Number}}}` : temp for `solution.u[n]`.
- `k       :: AbstractVector{Union{Number, AbstractVector{Number}}}` : stages.
- `problem :: AbstractInitialValueProblem`                           : initial value problem..
- `solver  :: ExplicitRungeKuttaSolver`                              : explicit Runge-Kutta solver.
"""
mutable struct ExplicitRungeKuttaCache{n_T, m_T, v_T, k_T} <: AbstractRungeKuttaCache
    n::n_T
    m::m_T
    v::v_T
    k::k_T
end

function ExplicitRungeKuttaCache(problem::AbstractInitialValueProblem, solver::ExplicitRungeKuttaSolver)
    @↓ u0 = problem
    @↓ s = solver.tableau
    n = m = 1
    v = similar(u0)
    k = Vector{eltype(u0)}(undef, s, length(u0))
    return ExplicitRungeKuttaCache(n, m, v, k)
end

"""
    ImplicitRungeKuttaCache <: AbstractRungeKuttaCache

A composite type for the [`AbstractRungeKuttaCache`](@ref) of an [`ImplicitRungeKuttaSolver`](@ref).

# Constructors
```julia
ImplicitRungeKuttaCache(n, m, v, k, Δk, J)
ImplicitRungeKuttaCache(problem, solver)
```

# Arguments
- `n       :: Integer`                                               : step counter.
- `m       :: Integer`                                               : adaptive correction counter.
- `v       :: AbstractVector{Union{Number, AbstractVector{Number}}}` : temp for `solution.u[n]`.
- `k       :: AbstractVector{Union{Number, AbstractVector{Number}}}` : stages.
- `Δk      :: AbstractVector{Union{Number, AbstractVector{Number}}}` : stages' correction.
- `J       :: AbstractMatrix{Union{Number, AbstractMatrix{Number}}}` : Jacobian of RHS function.
- `problem :: AbstractInitialValueProblem`                           : initial value problem..
- `solver  :: ImplicitRungeKuttaSolver`                              : implicit Runge-Kutta solver.
"""
mutable struct ImplicitRungeKuttaCache{n_T, m_T, v_T, k_T, Δk_T, J_T} <: AbstractRungeKuttaCache
    n::n_T
    m::m_T
    v::v_T
    k::k_T
    Δk::Δk_T
    J::J_T
end

function ImplicitRungeKuttaCache(problem::AbstractInitialValueProblem, solver::ImplicitRungeKuttaSolver)
    @↓ u0 = problem
    @↓ s = solver.tableau
    n = m = 1
    v = similar(u0)
    u0_T = eltype(u0)
    L = length(u0)
    k = Vector{u0_T}(undef, s, L)
    Δk = Vector{u0_T}(undef, s, L)
    J = Matrix{u0_T}(undef, L, L)
    return ImplicitRungeKuttaCache(n, m, v, k, Δk, J)
end

mutable struct ExplicitExponentialRungeKuttaCache{n_T, Q_T, α_T, β_T, γ_T, E_T, E2_T} <: AbstractRungeKuttaCache
    n::n_T
    Q::Q_T
    α::α_T
    β::β_T
    γ::γ_T
    E::E_T
    E2::E2_T
end

function ExplicitExponentialRungeKuttaCache(problem::AbstractInitialValueProblem, solver::ExplicitExponentialRungeKuttaSolver)
    @↓ L = problem.rhs
    @↓ h = solver.stepsize
    #
    n = 1
    #
    M = 16
    r = @. exp(π * 1im * ((1:M) - 0.5) / M)
    #
    Z = h * L .+ r'
    Q = vec(h/M * real.(sum( @. (exp(Z/2) - 1.0) / Z                         ; dims = 2)))
    α = vec(h/M * real.(sum( @. (-4.0 - Z + exp(Z) * (4.0 - 3Z + Z^2)) / Z^3 ; dims = 2)))
    β = vec(h/M * real.(sum( @. (2.0 + Z + exp(Z) * (-2.0 + Z)) / Z^3        ; dims = 2)))
    γ = vec(h/M * real.(sum( @. (-4.0 - 3Z - Z^2 + exp(Z) * (4.0 - Z)) / Z^3 ; dims = 2)))
    # L ≡ diagonal matrix
    E = @. exp(h * L)
    E2 = @. exp(h/2 * L)
    return ExplicitExponentialRungeKuttaCache(n, Q, α, β, γ, E, E2)
end

RungeKuttaCache(problem, solver::ExplicitRungeKuttaSolver) = ExplicitRungeKuttaCache(problem, solver)
RungeKuttaCache(problem, solver::ImplicitRungeKuttaSolver) = ImplicitRungeKuttaCache(problem, solver)
RungeKuttaCache(problem, solver::ExplicitExponentialRungeKuttaSolver) = ExplicitExponentialRungeKuttaCache(problem, solver)


text:
# ================================================================
#     AbstractGraphKernel
# ================================================================

"""
    abstract type AbstractGraphKernel <: KernelFunctions.Kernel

A kernel function between two graphs.

Subtypes of `AbstractGraphKernel` should implement `preprocessed_form` and `apply_preprocessed`.
When `(k::AbstractGraphKernel)(g1, g2)` is invoked on two graphs, then
```
apply_preprocessed(k, preprocessed_form(k, g1), preprocessed_form(k, g2))
```
is called to calculate the kernel function. Therefore one should implement `preprocessed_form`
that transforms a single graph into a suitable representation and `apply_preprocessed` that
takes the representations for both graphs and calculates the kernel function.

### See also
[`preprocessed_form`](@ref), [`apply_preprocessed`](@ref), [`KernelFunctions.Kernel`](@ref), [`kernel_matrix`](@ref), [`kernel_matrix_diag`](@ref)

"""
abstract type AbstractGraphKernel <: KernelFunctions.Kernel end

"""
    preprocessed_form(k::AbstractGraphKernel, g::AbstractGraph) = g

Transform a graph `g` into a suitable form for a graph kernel `k`

When calculating a pairwise kernel matrix for multiple graphs, this preprocessed form
allows us to calculate the transformation only once for each graph, so that we can cache
the result. By default this simply returns `g` without any transformation.

When implementing a custom graph kernel, it might be a good idea to implement this
method.

### See also
[`AbstractGraphKernel`](@ref), [`apply_preprocessed`](@ref)
"""
preprocessed_form(::AbstractGraphKernel, g::AbstractGraph) = g

function (kernel::AbstractGraphKernel)(g1::AbstractGraph, g2::AbstractGraph)

    return apply_preprocessed(kernel, preprocessed_form(kernel, g1), preprocessed_form(kernel, g2))
end

## ---------------------------------------------------------------
##       kernelmatrix & kernelmatrix_diag
## ---------------------------------------------------------------

function _map_preprocessed_form(kernel::AbstractGraphKernel, graphs)

    # TODO we should be able to avoid collecting the graphs
    # but currently ThreadX cannot split them otherwise,
    # maybe we can create some wrapper type that is splittable around graphs
    return ThreadsX.map(g -> preprocessed_form(kernel, g), collect(graphs))
end

"""
    kernelmatrix(kernel, graphs)
Return a matrix of running the kernel on all pairs of graphs.

### See also
[`kernelmatrix_diag`](@ref)
"""
function KernelFunctions.kernelmatrix(kernel::AbstractGraphKernel, graphs::AbstractVector)

    pre = _map_preprocessed_form(kernel, graphs)

    # this simply a guard to make the code more type save, maybe we can get
    # rid of it at some point
    return _kernelmatrix_from_preprocessed(kernel, pre)
end

function _kernelmatrix_from_preprocessed(kernel, pre)

    n = length(pre)

    # TODO maybe we should make the matrix only symmetric afterwards
    # so that we avoid false sharing when using multiple threads
    # TODO create some triangle generator instead of allocating a vector
    # TODO apparently ThreadsX can do load balancing so we should consider that here
    G = Matrix{Float64}(undef, n, n)
    indices = [(i, j) for i in 1:n for j in i:n]
    Threads.@threads for idx in indices
        i, j = idx
        @inbounds v = apply_preprocessed(kernel, pre[i], pre[j])
        @inbounds G[i, j] = v
        @inbounds G[j, i] = v
    end

    return G
end

function KernelFunctions.kernelmatrix_diag(kernel::AbstractGraphKernel, graphs::AbstractVector)

    n = length(graphs)
    pre = _map_preprocessed_form(kernel, graphs)

    D = Vector{Float64}(undef, n)
    Threads.@threads for i in 1:n
        @inbounds D[i] = apply_preprocessed(kernel, pre[i], pre[i])
    end
    return D
end

function KernelFunctions.kernelmatrix(kernel::AbstractGraphKernel, graphs1::AbstractVector, graphs2::AbstractVector)

    n_rows = length(graphs1)
    n_cols = length(graphs2)

    M = Matrix{Float64}(undef, n_rows, n_cols)

    pre1 = _map_preprocessed_form(kernel, graphs1)
    pre2 = _map_preprocessed_form(kernel, graphs2)

    Threads.@threads for i in 1:n_rows
        for j in 1:n_cols
            @inbounds M[i, j] = apply_preprocessed(kernel, pre1[i], pre2[j])
        end
    end

    return M
end


text:
# This module is needed until AbstractTrees.jl#37 is fixed.
# (PR: https://github.com/Keno/AbstractTrees.jl/pull/38)
# because currently `print_tree` does not respect `maxdepth`.
# This just implements the changes in the above PR.
# Code in this file is modified from AbstractTrees.jl
# See LICENSE for a copy of its MIT license.
module TreePrint

using AbstractTrees: printnode, treekind, IndexedTree, children

# Printing
struct TreeCharSet
    mid
    terminator
    skip
    dash
    ellipsis
end

# Default charset
TreeCharSet() = TreeCharSet('├','└','│','─','…')

function print_prefix(io, depth, charset, active_levels)
    for current_depth in 0:(depth-1)
        if current_depth in active_levels
            print(io,charset.skip," "^(textwidth(charset.dash)+1))
        else
            print(io," "^(textwidth(charset.skip)+textwidth(charset.dash)+1))
        end
    end
end

@doc raw"""
# Usage
Prints an ASCII formatted representation of the `tree` to the given `io` object.
By default all children will be printed up to a maximum level of 5, though this
valud can be overriden by the `maxdepth` parameter. The charset to use in
printing can be customized using the `charset` keyword argument.

# Examples
```julia
julia> print_tree(STDOUT,Dict("a"=>"b","b"=>['c','d']))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
├─ b
│  ├─ c
│  └─ d
└─ a
   └─ b

julia> print_tree(STDOUT,Dict("a"=>"b","b"=>['c','d']);
        charset = TreeCharSet('+','\\','|',"--"))
Dict{String,Any}("b"=>['c','d'],"a"=>"b")
+-- b
|   +-- c
|   \-- d
\-- a
   \-- b
```

"""
print_tree

function _print_tree(printnode::Function, io::IO, tree, maxdepth = 5; depth = 0, active_levels = Int[],
    charset = TreeCharSet(), withinds = false, inds = [], from = nothing, to = nothing, roottree = tree)
    nodebuf = IOBuffer()
    isa(io, IOContext) && (nodebuf = IOContext(nodebuf, io))
    if withinds
        printnode(nodebuf, tree, inds)
    else
        tree != roottree && isa(treekind(roottree), IndexedTree) ?
            printnode(nodebuf, roottree[tree]) :
            printnode(nodebuf, tree)
    end
    str = String(take!(isa(nodebuf, IOContext) ? nodebuf.io : nodebuf))
    for (i,line) in enumerate(split(str, '\n'))
        i != 1 && print_prefix(io, depth, charset, active_levels)
        println(io, line)
    end
    depth > maxdepth && return
    c = isa(treekind(roottree), IndexedTree) ?
        childindices(roottree, tree) : children(roottree, tree)
    if c !== ()
        s = Iterators.Stateful(from === nothing ? pairs(c) : Iterators.Rest(pairs(c), from))
        while !isempty(s)
            ind, child = popfirst!(s)
            ind === to && break
            active = false
            child_active_levels = active_levels
            print_prefix(io, depth, charset, active_levels)
            if isempty(s)
                print(io, charset.terminator)
            else
                print(io, charset.mid)
                child_active_levels = push!(copy(active_levels), depth)
            end
            print(io, charset.dash, ' ')
            print_tree(depth == maxdepth ? (io, val) -> print(io, charset.ellipsis) : printnode,
                io, child, maxdepth; depth = depth + 1,
                active_levels = child_active_levels, charset = charset, withinds=withinds,
                inds = withinds ? [inds; ind] : [], roottree = roottree)
        end
    end
end
print_tree(f::Function, io::IO, tree, args...; kwargs...) = _print_tree(f, io, tree, args...; kwargs...)
print_tree(io::IO, tree, args...; kwargs...) = print_tree(printnode, io, tree, args...; kwargs...)
print_tree(tree, args...; kwargs...) = print_tree(stdout::IO, tree, args...; kwargs...)
end


text:
#takes a dataframe and returns a dataframe with upper and lower CI
function confidence_interval(data::DataFrame, num_particles::Int64,
    n_possible_objects::Int64, n_boot::Int64, cil::Float64)

    col_names = []
    for obj = 1:n_possible_objects
        push!(col_names, "fa_" * string(obj) * "_lower")
        push!(col_names, "fa_" * string(obj) * "_upper")
        push!(col_names, "m_" * string(obj) * "_lower")
        push!(col_names, "m_" * string(obj) * "_upper")
    end
    data_ci = DataFrame()
    for i=1:length(col_names)
        data_ci[:, col_names[i]] = Float64[]
    end

    for v = 1:num_videos
        weights = Array{Float64}(undef, num_particles)
        for j = 1:num_particles
            str = "weight_" * string(j)
            weights[j] = data[v, str]
        end
        normalized_log_weights = weights ./ (sum(weights))

        new_row = Float64[] #get bootstrapped values
        for obj = 1:n_possible_objects
            for e = 1:2 #false alarms and misses
                values = Float64[]
                for j = 1:num_particles
                    str = e==1 ? "fa_" * string(obj) * "_" * string(j) : "m_" * string(obj) * "_" * string(j)
                    push!(values, data[v, str])
                end
                df = DataFrame(weights = exp.(normalized_log_weights), values = values)
                bs = bootstrap(weighted_average, df, BasicSampling(n_boot))
                bci = confint(bs, PercentileConfInt(cil))
                lower_ci = bci[1][2]
                upper_ci = bci[1][3]
                push!(new_row, lower_ci)
                push!(new_row, upper_ci)
            end
        end
        push!(data_ci, new_row)
    end
    return data_ci
end

#dataframe has a column named "weights" and a column named "values" and returns a weighted average of the values
function weighted_average(data::DataFrame)
    return sum(data[!, "values"] .* (data[!, "weights"]./sum(data[!, "weights"])))
end

################################################################################
#takes a dataframe. returns the MSE of the weighted average V
#each row of the dataframe is a particle. columns are the weight of the particle,
#the fa and m, and the gt_fa, gt_m
function make_MSE(data::DataFrame)

    n_possible_objects = convert(Int64,(size(df)[2]-1)/4)
    num_particles = size(df)[1]

    weights = data[!, "weights"]
    v_matrixes = Array{Matrix{Float64}}(undef, num_particles)
    for j = 1:num_particles
        v_matrix = zeros(params.n_possible_objects, 2)
        for obj = 1:n_possible_objects
            for e = 1:2 #false alarms and misses
                str = e==1 ? "fa_" * string(obj) : "m_" * string(obj)
                v_matrix[obj, e] = data[j, str]
            end
        end
        v_matrixes[j] = v_matrix
    end
    normalized_log_weights = weights ./ (sum(weights))
    avg_v = sum(exp.(normalized_log_weights)./sum(exp.(normalized_log_weights)) .* v_matrixes)#change out of log

    #get the gt_v_matrix
    gt_v_matrix = zeros(params.n_possible_objects, 2)
    for obj = 1:n_possible_objects
        for e = 1:2
            gt_str = e==1 ? "gt_fa_" * string(obj)  : "gt_m_" * string(obj)
            gt_v_matrix[obj, e] = data[1, gt_str] #every row should have the same thing for gt, so just take the first one
        end
    end

    #now calculate the MSE for this avg_v
    return sum((avg_v .- gt_v_matrix).^2) / length(gt_v_matrix)
end

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

#takes two dataframes: the dataframe with particles, and the ground truth
#returns a dataframe with three columns: mean, upper, and lower CI for overall MSE of V
function MSE_and_confidence_interval(data::DataFrame, gt_data::DataFrame, num_particles::Int64,
    n_possible_objects::Int64, n_boot::Int64, cil::Float64)

    data_ci = DataFrame(MSE = Float64[], upper_MSE = Float64[], lower_MSE = Float64[])

    for v = 1:num_videos

        #set up the dataframe where each row is a particle
        col_names = []
        push!(col_names, "weights")
        for obj = 1:n_possible_objects
            push!(col_names, "fa_" * string(obj))
            push!(col_names, "m_" * string(obj))
            push!(col_names, "gt_fa_" * string(obj))
            push!(col_names, "gt_m_" * string(obj))
        end
        df = DataFrame()
        for i=1:length(col_names)
            df[:, col_names[i]] = Float64[]
        end

        for j = 1:num_particles
            new_row = []
            str = "weight_" * string(j)
            push!(new_row, data[v, str])
            for obj = 1:n_possible_objects
                str = "fa_" * string(obj) * "_" * string(j)
                push!(new_row, data[v, str])
                str = "m_" * string(obj) * "_" * string(j)
                push!(new_row, data[v, str])
                str = "gt_fa_" * string(obj)
                push!(new_row, gt_data[v, str])
                str = "gt_m_" * string(obj)
                push!(new_row, gt_data[v, str])
            end
            push!(df, new_row)
        end

        new_row_ci = []
        MSE = make_MSE(df)
        push!(new_row_ci, MSE)
        bs = bootstrap(make_MSE, df, BasicSampling(n_boot))
        bci = confint(bs, PercentileConfInt(cil))
        lower_ci = bci[1][2]
        upper_ci = bci[1][3]
        push!(new_row_ci, lower_ci)
        push!(new_row_ci, upper_ci)

        push!(data_ci, new_row_ci)
    end
    return data_ci
end

max_stars_repo_name: zhangir-azerbayev/MetaGen

max_issues_repo_name: zhangir-azerbayev/MetaGen

max_forks_repo_name: zhangir-azerbayev/MetaGen


text:
using Documenter
using CCDReduction

DocMeta.setdocmeta!(CCDReduction, :DocTestSetup, :(using CCDReduction); recursive = true)

makedocs(modules = [CCDReduction],
    sitename = "CCDReduction.jl",
    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
    authors = "Siddharth Lal",
    strict = true,
    pages = [
        "Home" => "index.md",
        "API/Reference" => "api.md"
    ],
)

deploydocs(;
    repo = "github.com/JuliaAstro/CCDReduction.jl.git",
    push_preview = true,
    devbranch = "main"
)


text:
export iso, iso2, iso_matrix, is_iso, info_map, uhash, degdeg, fast_iso_test

const iso_err_msg = "The graphs are not isomorphic"

"""
`iso_matrix(G,H)` returns a permutation matrix `P` such that
`A*P==P*B` where `A` is the adjacency matrix of `G` and `B` is the
adjacency matrix of `H`.  If the graphs are not isomorphic, an
error is raised.
"""
function iso_matrix(G::SimpleGraph, H::SimpleGraph)
    n = NV(G)

    if !fast_iso_test_basic(G, H)
        error(iso_err_msg)
    end

    m = Model(get_solver())

    @variable(m, P[1:n, 1:n], Bin)
    A = adjacency(G)
    B = adjacency(H)

    for i = 1:n
        for k = 1:n
            @constraint(
                m,
                sum(A[i, j] * P[j, k] for j = 1:n) == sum(P[i, j] * B[j, k] for j = 1:n)
            )
        end
        @constraint(m, sum(P[i, j] for j = 1:n) == 1)
        @constraint(m, sum(P[j, i] for j = 1:n) == 1)
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error(iso_err_msg)
    end

    return Int.(value.(P))
end

"""
`is_iso(G,H,d)` checks if `d` is an isomorphism from `G` to `H`.
"""
function is_iso(G::SimpleGraph, H::SimpleGraph, d::Dict)::Bool
    n = NV(G)

    # standard quick check
    if n != NV(H) || NE(G) != NE(H) || deg(G) != deg(H) || length(d) != n
        return false
    end

    # check keys and values in d match VG and VH, respectively
    for k in keys(d)
        if !has(G, k)
            return false
        end
    end

    for v in values(d)
        if !has(H, v)
            return false
        end
    end

    # Check edges match. Only need to go one way since we determined
    # that the two graphs have the same number of edges.
    EE = elist(G)
    for e in EE
        v = e[1]
        w = e[2]
        if !has(H, d[v], d[w])
            return false
        end
    end

    return true
end

"""
`is_iso(G,H)` returns `true` if the two graphs are isomorphic and `false`
otherwise.
"""
function is_iso(G::SimpleGraph, H::SimpleGraph)::Bool
    try
        f = iso(G, H)
    catch
        return false
    end
    return true
end

"""
`fast_iso_test_basic(G,H)` is a very quick test to rule out graphs being
isomorphic by considering the number of vertices, number of edges, and
degree sequences. Returns `false` if the graphs fail this very basic
test of isomorphism. A `true` result does *not* imply the graphs are
isomorphic.
"""
function fast_iso_test_basic(G::SimpleGraph, H::SimpleGraph)::Bool
    if NV(G) != NV(H) || NE(G) != NE(H) || deg(G) != deg(H)
        return false
    end
    return true
end

"""
`fast_iso_test(G,H)` applies various heuristics to see if
the graphs *might* be isomorphic. A `false` return certifies the
graphs are **not** isomorphic; a `true` result indicates the
graphs might be (indeed, likely are) isomorphic.
"""
function fast_iso_test(G::SimpleGraph, H::SimpleGraph)::Bool
    return fast_iso_test_basic(G, H) && uhash(G) == uhash(H)
end

"""
`degdeg(G)` creates an n-by-n matrix where the nonzero entries in a
row are the degrees of the neighbors of that vertex. The rows are
lexicographically sorted.
"""
function degdeg(G::SimpleGraph)
    if cache_check(G, :degdeg)
        return cache_recall(G, :degdeg)
    end

    n = NV(G)
    result = zeros(Int, n, n)
    VV = vlist(G)

    for k = 1:n
        v = VV[k]
        dv = [deg(G, w) for w in G[v]]
        dv = sort(dv)
        dv = [dv; zeros(Int, n - deg(G, v))]
        result[k, :] = dv
    end

    result = sortslices(result, dims = 1)
    cache_save(G, :degdeg, result)
    return result
end

"""
`info_map(G)`:
We create a dictionary mapping the vertices of `G` to 128-bit integer
values in such a way that twin vertices *will* have the same value
but, we hope, nontwin vertices will have different values. (By *twin*
we mean a pair of vertices such that there is an automorphism of the
graph mapping one to the other.)
"""
function info_map(G::SimpleGraph)
    if cache_check(G, :info_map)
        return cache_recall(G, :info_map)
    end
    n = NV(G)
    VV = vlist(G)
    H = G'

    # Section 1: Neighborhood degrees
    M1 = zeros(Int, n, n - 1)
    for k = 1:n
        v = VV[k]
        rv = [deg(G, w) for w in G[v]]
        rv = [sort(rv); zeros(Int, n - 1 - deg(G, v))]
        M1[k, :] = rv
    end

    # Section 2: Complement degrees
    M2 = zeros(Int, n, n - 1)
    for k = 1:n
        v = VV[k]
        rv = [deg(H, w) for w in H[v]]
        rv = [sort(rv); zeros(Int, n - 1 - deg(H, v))]
        M2[k, :] = rv
    end

    # Section 3: Distances
    M3 = zeros(Int, n, n - 1)
    for k = 1:n
        v = VV[k]
        dv = dist(G, v)
        M3[k, :] = sort(collect(values(dv)))[1:n-1]
    end

    # Section 4: Distances in complement
    M4 = zeros(Int, n, n - 1)
    for k = 1:n
        v = VV[k]
        dv = dist(H, v)
        M4[k, :] = sort(collect(values(dv)))[1:n-1]
    end

    M = [M1 M2 M3 M4]

    # println(sortrows(M))  # debug

    T = eltype(G)
    d = Dict{T,Int128}()
    for k = 1:n
        v = VV[k]
        d[v] = Int128(hash(M[k, :]))
    end
    cache_save(G, :info_map, d)
    return d
end

"""
`iso(G,H)` finds an isomorphism from `G` to `H` (or throws an error if
the graphs are not isomorphic). Returns a dictionary mapping the
vertices of `G` to `H`.

See also: `iso2`.
"""
function iso(G::SimpleGraph, H::SimpleGraph)

    # quickly rule out some easily detected nonisomorphic graphs
    if !fast_iso_test(G, H) || uhash(G) != uhash(H)
        error(iso_err_msg)
    end

    VG = vlist(G)
    VH = vlist(H)
    n = NV(G)

    TG = eltype(G)
    TH = eltype(H)

    dG = info_map(G)
    dH = info_map(H)

    # Create mappings from vertex key values back to the vertices themselves

    xG = Dict{Int128,Set{TG}}()
    xH = Dict{Int128,Set{TH}}()

    valsG = sort(collect(values(dG)))
    valsH = sort(collect(values(dH)))

    for x in valsG
        xG[x] = Set{TG}()
        xH[x] = Set{TH}()
    end

    for v in VG
        x = dG[v]
        push!(xG[x], v)
    end

    for v in VH
        x = dH[v]
        push!(xH[x], v)
    end

    MOD = Model(get_solver())

    @variable(MOD, x[VG, VH], Bin)

    for v in VG
        @constraint(MOD, sum(x[v, VH[k]] for k = 1:n) == 1)
    end

    for w in VH
        @constraint(MOD, sum(x[VG[k], w] for k = 1:n) == 1)
    end

    for v in VG
        for w in VH
            @constraint(
                MOD,
                sum(has(G, v, VG[k]) * x[VG[k], w] for k = 1:n) ==
                sum(x[v, VH[k]] * has(H, VH[k], w) for k = 1:n)
            )
        end
    end

    # Add contraints based on vertex numbers

    for val in unique(valsG)
        SG = collect(xG[val])
        SH = collect(xH[val])
        s = length(SG)
        @constraint(MOD, sum(x[u, v] for u in SG for v in SH) == s)
    end

    optimize!(MOD)
    status = Int(termination_status(MOD))

    if status != 1
        error(iso_err_msg)
    end

    X = value.(x)

    result = Dict{eltype(G),eltype(H)}()

    for v in VG
        for w in VH
            if X[v, w] > 0
                result[v] = w
            end
        end
    end
    return result
end

"""
`iso2(G,H)` is another version of `iso(G,H)` that does much less preprocessing.
It may be faster for highly symmetric graphs (e.g., vertex transitive graphs).
"""
function iso2(G::SimpleGraph{S}, H::SimpleGraph{T}) where {S,T}
    if NV(G) != NV(H)
        error(iso_err_msg)
    end

    m = Model(get_solver())

    @variable(m, P[G.V, H.V], Bin)

    for v in G.V
        for x in H.V
            @constraint(
                m,
                sum(G[v, w] * P[w, x] for w in G.V) == sum(P[v, y] * H[y, x] for y in H.V)
            )
        end
    end

    for v in G.V
        @constraint(m, sum(P[v, x] for x in H.V) == 1)
    end

    for x in H.V
        @constraint(m, sum(P[v, x] for v in G.V) == 1)
    end

    optimize!(m)
    status = Int(termination_status(m))

    if status != 1
        error(iso_err_msg)
    end

    PP = Int.(value.(P))

    result = Dict{S,T}()

    for v in G.V
        for x in H.V
            if PP[v, x] > 0
                result[v] = x
            end
        end
    end

    return result

end

using LinearAlgebra

function matrix_moments(A, max_exp::Int = 10)
    result = zeros(Int, max_exp)
    B = copy(A)

    for t = 1:max_exp
        B = A * B
        result[t] = tr(B)
    end
    return result
end

"""
`uhash(G)` creates an `UInt64` hash of the graph such that isomorphic
graphs will produce the same value. We hope that nonisomorphic graphs
will create different values, but, alas, that need not be the case.
"""
function uhash(G::SimpleGraph)
    if cache_check(G, :uhash)
        return cache_recall(G, :uhash)
    end

    v1 = sort(collect(values(info_map(G))))

    A = adjacency(G)
    B = deepcopy(A)
    depth = min(10, NV(G))

    v2 = matrix_moments(adjacency(G))
    v3 = matrix_moments(laplace(G))

    vv = [v1; v2; v3]

    x = hash(vv)
    cache_save(G, :uhash, x)
    return x
end

max_stars_repo_name: scheinerman/SimpleGraphAlgorithms

max_issues_repo_name: scheinerman/SimpleGraphAlgorithms

max_forks_repo_name: scheinerman/SimpleGraphAlgorithms


text:
"""
updates B.ss_params.κs_β for the spin system portion of compound B.
"""
# function updateβ!(B::NMRSpectraSimulator.FIDModelType{T,SST},
#     ME::Vector{Vector{Vector{Int}}},
#     p::Vector{T},
#     j_st::Int)::Int where {T <: Real, SST}
#
#     # must have the same number of spin systems.
#     @assert length(ME) == length(B.ss_params.κs_β)
#
#     j = j_st
#     for i = 1:length(B.ss_params.κs_β)
#         #for l = 1:length(B.ss_params.κs_β[i])
#
#         if isempty(ME[i])
#             for l = 1:length(B.ss_params.κs_β[i])
#                 j += 1
#
#                 B.ss_params.κs_β[i][l] = p[j]
#             end
#         else
#
#             for dest_inds in ME[i]
#                 j += 1
#
#                 B.ss_params.κs_β[i][l] = p[j]
#             end
#         end
#     end
#
#     return j
# end

function calibrateregions( y::Vector{Complex{T}},
    U_y, P_y, cost_inds_set, Δ_shifts, As, fs, SW::T, w;
    max_iters = 5000,
    xtol_rel = 1e-7,
    ftol_rel = 1e-12,
    κ_lb_default = 0.2,
    κ_ub_default = 50.0,
    λ_each_lb = 0.9,
    λ_each_ub = 1.1 ) where T <: Real

    #
    N_compounds = length(As)
    N_regions = length(cost_inds_set)

    Es = collect( NMRSpectraSimulator.κCompoundFIDType(As[i]) for i = 1:N_compounds )

    p_star_set = Vector{Vector{T}}(undef, N_regions)
    κ_BLS_set = Vector{Vector{T}}(undef, N_regions)

    d_star_set = Vector{Vector{T}}(undef, N_regions)
    β_star_set = Vector{Vector{T}}(undef, N_regions)
    λ_star_set = Vector{Vector{T}}(undef, N_regions)
    proxies_set = Vector{Function}(undef, N_regions)

    for r = 1:N_regions

        y_cost = y[cost_inds_set[r]]
        U_cost = U_y[cost_inds_set[r]]
        P_cost = P_y[cost_inds_set[r]]

        ###
        LS_inds = 1:length(U_cost)

        println("Timing: runalignment(), r = $(r)")
        @time p_star, q, κ_BLS, getshiftfunc, getβfunc, getλfunc,
        obj_func, N_vars_set = NMRCalibrate.runalignment(Δ_shifts,
        U_cost, y_cost, LS_inds, Es, As, fs, SW;
        max_iters = max_iters,
        xtol_rel = xtol_rel,
        ftol_rel = ftol_rel,
        w = w,
        κ_lb_default = κ_lb_default,
        κ_ub_default = κ_ub_default,
        λ_each_lb = λ_each_lb,
        λ_each_ub = λ_each_ub)

        p_star_set[r] = p_star
        κ_BLS_set[r] = κ_BLS
        proxies_set[r] = q

        d_star_set[r] = getshiftfunc(p_star)
        β_star_set[r] = getβfunc(p_star)
        λ_star_set[r] = getλfunc(p_star)
    end

    return cost_inds_set, p_star_set, κ_BLS_set, d_star_set, β_star_set, λ_star_set, proxies_set
end

function savefigfitresult(save_path::String,
    title_string::String,
    q_U::Vector{T},
    P,
    P_cost,
    y_cost::Vector{T};
    initial_fit::Vector{T} = zeros(T, 0),
    P_y = P_cost,
    y = y_cost) where T

    plot_obj = Plots.plot( P_y,
        y,
        title = title_string,
        label = "Data",
        seriestype = :line,
        ticks = :native,
        xlims = (P[1],P[end]),
        hover = P,
        linewidth = 4,
        xflip = true,
        size = (1600, 900))

    Plots.plot!(plot_obj, P, q_U, label = "Fitted model",
        seriestype = :line,
        linestyle = :dot,
        xflip = true,
        linewidth = 4)

    if !isempty(initial_fit)
        Plots.plot!(plot_obj, P, initial_fit, label = "Initial model",
            seriestype = :line,
            linestyle = :dash,
            xflip = true,
            linewidth = 4)
    end

    if !isempty(y_cost)
        Plots.plot!(plot_obj, P_cost, y_cost, label = "Fit positions",
            markershape = :circle,
            seriestype = :scatter)
    end

    Plots.savefig(plot_obj, save_path)
end


text:
function save_single_features_instance!(db::SQLite.DB, feature_names, use_network, name, scenario, source_type, source_path, pfn_path)
    println("Saving instance features: ($name, $scenario)")
    network = nothing
    if use_network
        network = ismissing(pfn_path) ? PowerFlowNetwork(source_path, source_type) : deserialize(pfn_path)
    end
    feature_values = Dict()
    query = "UPDATE instances SET "
    for feature_name in feature_names
        feature_info = _feature_info_dict[feature_name]
        if feature_info[1] == :graph
            g = SimpleGraph(network)
            feature_value = get_single_feature_graph(feature_name, g)
        elseif feature_info[1] == :network
            feature_value = get_single_feature_network(feature_name, network)
        elseif feature_info[1] == :source_path
            feature_value = get_single_feature_source_path(feature_name, source_path)
        end
        if feature_info[2]
            query *= """
            $(feature_name)_max = $(feature_value[1]),
            $(feature_name)_min = $(feature_value[2]),
            $(feature_name)_mean = $(feature_value[3]),
            """
        else
            query *= "$(feature_name) = $(feature_value),"
        end
    end
    query = query[begin:end - 1] * " WHERE name = '$name' AND scenario = $scenario"
    execute_query(db, query)
end

function save_single_features_instance_dfrow!(db::SQLite.DB, feature_names, use_network, row)
    save_single_features_instance!(db, feature_names, use_network, row[:name], row[:scenario], row[:source_type], row[:source_path], row[:pfn_path])
end

function save_single_features_instances!(db::SQLite.DB;
                                         feature_names, min_nv=typemin(Int), max_nv=typemax(Int),
                                         recompute=false, subset=nothing)
    query = "SELECT name, scenario, source_type, source_path, pfn_path FROM instances WHERE nbus >= $min_nv AND nbus <= $max_nv"
    if !recompute
        for feature_name in feature_names
            query *= " AND $feature_name IS NULL"
        end
    end
    if !isnothing(subset)
        query *= " AND id IN ($(join(subset, ',')))"
    end
    use_network = any([_feature_info_dict[f][1] in (:graph, :network) for f in feature_names])
    results = DBInterface.execute(db, query) |> DataFrame
    println("Saving instance features: feature_names=$(feature_names), min_nv=$min_nv, max_nv=$max_nv, recompute=$recompute")
    println("subset\n$subset\nend subset")
    save_function!(row) = save_single_features_instance_dfrow!(db, feature_names, use_network, row)
    save_function!.(eachrow(results[!, [:name, :scenario, :source_type, :source_path, :pfn_path]]))
end

max_stars_repo_name: charlyalizadeh/PowerFlowNetworks

max_issues_repo_name: charlyalizadeh/PowerFlowNetworks

max_forks_repo_name: charlyalizadeh/PowerFlowNetworks


text:
using NumericFuns
using Base.Test
using Compat

println("    macros")
test_unary(x) = -x
test_binary(x, y) = x + y
@functor1(TestUnary, test_unary, Number)
@functor2(TestBinary, test_binary, Number)

@test evaluate(TestUnary(), 3) == -3
@test evaluate(TestBinary(), 3, 4) == 7

println("    arithmetic operators")

@test evaluate(Negate(), 2) == -2

for (F, sf) in [(Add, +),
                (Subtract, -),
                (Multiply, *),
                (Divide, /),
                (RDivide, \),
                (Pow, ^)]

    @test evaluate(F(), 6, 2) == sf(6, 2)
    @test evaluate(F(), 5.0, 2.0) == sf(5.0, 2.0)
end

println("    comparison operators")

for (F, sf) in [(LT, <),
                (GT, >),
                (LE, <=),
                (GE, >=),
                (EQ, ==),
                (NE, !=)]

    @test evaluate(F(), 1, 1) == sf(1, 1)
    @test evaluate(F(), 1, 2) == sf(1, 2)
    @test evaluate(F(), 2, 1) == sf(2, 1)
end

println("    logical operators")

@test evaluate(Not(), true) == false
@test evaluate(Not(), false) == true

@test evaluate(And(), false, false) == false
@test evaluate(And(), false, true)  == false
@test evaluate(And(), true,  false) == false
@test evaluate(And(), true,  true)  == true

@test evaluate(Or(), false, false) == false
@test evaluate(Or(), false, true)  == true
@test evaluate(Or(), true,  false) == true
@test evaluate(Or(), true,  true)  == true

println("    bitwise operators")

@test evaluate(BitwiseNot(), 5) == ~5

for (F, sf) in [(BitwiseAnd, &),
                (BitwiseOr,  |),
                (BitwiseXor, $)]

    @test evaluate(F(), 5, 9) == sf(5, 9)
end

println("    arithmetic functions")

for (F, sf) in [(DivFun, div),
                (FldFun, fld),
                (RemFun, rem),
                (ModFun, mod),
                (MaxFun, max),
                (MinFun, min)]

    @test evaluate(F(),  14, 3) == sf( 14, 3)
    @test evaluate(F(), -14, 3) == sf(-14, 3)
end

for (F, sf) in [(AbsFun, abs),
                (Abs2Fun, abs2),
                (RealFun, real),
                (ImagFun, imag)]

    @test evaluate(F(), 5.0) == sf(5.0)
    @test evaluate(F(), -5.0) == sf(-5.0)
    @test evaluate(F(), 3.0 + 4.0im) == sf(3.0 + 4.0im)
end

for (F, sf) in [(SignFun, sign),
                (SignbitFun, signbit)]

    @test evaluate(F(),  0) == sf(0)
    @test evaluate(F(),  1) == sf(1)
    @test evaluate(F(), -1) == sf(-1)
end

println("    rounding functions")

for (F, sf) in [(FloorFun, floor),
                (CeilFun, ceil),
                (TruncFun, trunc),
                (RoundFun, round),
                (IfloorFun, x->floor(Integer, x)),
                (IceilFun, x->ceil(Integer, x)),
                (ItruncFun, x->trunc(Integer, x)),
                (IroundFun, x->round(Integer, x))]

    @test evaluate(F(),  3.4) === sf(3.4)
    @test evaluate(F(),  3.8) === sf(3.8)
    @test evaluate(F(), -3.4) === sf(-3.4)
    @test evaluate(F(), -3.8) === sf(-3.8)
end

println("    algebraic functions")

for (F, sf) in [(SqrtFun, sqrt), (CbrtFun, cbrt)]
    @test evaluate(F(), 64.0) == sf(64.0)
end

@test evaluate(HypotFun(), 3.0, 4.0) == hypot(3.0, 4.0)

println("    trancendental functions")

for (F, sf) in [(ExpFun, exp),
                (Exp2Fun, exp2),
                (Exp10Fun, exp10),
                (Expm1Fun, expm1),
                (LogFun, log),
                (Log2Fun, log2),
                (Log10Fun, log10),
                (Log1pFun, log1p),
                (SinFun, sin),
                (SinhFun, sinh),
                (SindFun, sind),
                (SincFun, sinc),
                (SinpiFun, sinpi),
                (AtanFun, atan),
                (AtandFun, atand)]

    @test evaluate(F(), 3.0) == sf(3.0)
end

@test evaluate(Atan2Fun(), 3.0, 4.0) == atan2(3.0, 4.0)

println("    special functions")

for (Fun, sf) in [(ErfFun(), erf),
                  (GammaFun(), gamma),
                  (AiryFun(), airy),
                  (Besselj0Fun(), besselj0),
                  (BesseljFun(2), x -> besselj(2, x))]

    @test evaluate(Fun, 2.0) == sf(2.0)
    @test evaluate(Fun, 1.0 + 2.0im) == sf(1.0 + 2.0im)
end


text:
#####
##### Five nearest neighbor stuff
#####

#=

#########################
#       #       #       #
#       #       #       #
# (0,0) # (0,1) # (0,2) #
#       #       #       #
#       #       #       #
#############################
    #       #       #       #
    #       #       #       #
    # (1,0) # (1,1) # (1,2) #
    #       #       #       #
    #       #       #       #
    #########################
    #       #       #       #
    #       #       #       #
    # (2,0) # (2,1) # (2,2) #
    #       #       #       #
    #       #       #       #
#############################
#       #       #       #
#       #       #       #
# (3,0) # (3,1) # (3,2) #
#       #       #       #
#       #       #       #
#########################

=#

struct Rect5 <: AbstractRectangular
    links::Int
    io_links::Int
end

# Encode hexagonal directions using angles.
directions(::Rect5) = ("north_left", "north_right", "north", "east", "west", "south_left", "south_right", "south")
function procrules(style::Rect5)
    # Set things up by row
    mod0_skeleton = (
        ((1, 0), "north_right_out", "south_left_in"),
        ((1,-1), "north_left_out", "south_right_in"),
        ((-1,0), "south_out", "north_in"),
        ((0, 1), "east_out", "west_in"),
        ((0,-1), "west_out", "east_in"),
    )
    mod0_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in mod0_skeleton, i in 0:style.links-1
    ])
    mod0_rule = ConnectionRule(
        mod0_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> mod(x[1], 4) == 0
    )

    # Mod 1 row
    mod1_skeleton = (
        ((-1,0), "south_left_out", "north_right_in"),
        ((-1,1), "south_right_out", "north_left_in"),
        ((1,0), "north_out", "south_in"),
        ((0, 1), "east_out", "west_in"),
        ((0,-1), "west_out", "east_in"),
    )
    mod1_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in mod1_skeleton, i in 0:style.links-1
    ])
    mod1_rule = ConnectionRule(
        mod1_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> mod(x[1], 4) == 1
    )

    # mod 2 row
    mod2_skeleton = (
        ((1, 1), "north_right_out", "south_left_in"),
        ((1, 0), "north_left_out", "south_right_in"),
        ((-1,0), "south_out", "north_in"),
        ((0, 1), "east_out", "west_in"),
        ((0,-1), "west_out", "east_in"),
    )
    mod2_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in mod2_skeleton, i in 0:style.links-1
    ])
    mod2_rule = ConnectionRule(
        mod2_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> mod(x[1], 4) == 2
    )

    # Mod 3 row
    mod3_skeleton = (
        ((-1,-1), "south_left_out", "north_right_in"),
        ((-1,0), "south_right_out", "north_left_in"),
        ((1,0), "north_out", "south_in"),
        ((0, 1), "east_out", "west_in"),
        ((0,-1), "west_out", "east_in"),
    )
    mod3_offsets = squash([
        Offset(a, "$b[$i]", "$c[$i]")
        for (a,b,c) in mod3_skeleton, i in 0:style.links-1
    ])
    mod3_rule = ConnectionRule(
        mod3_offsets,
        source_filter = filter_proc,
        dest_filter = filter_proc,
        address_filter = x -> mod(x[1], 4) == 3
    )

    return [mod0_rule, mod1_rule, mod2_rule, mod3_rule]
end

function port_boundaries(::Rect5, orientation)
    coords = Dict(
        "east" => ((1,0), (1,1)),

        "north_right"    => ((1, 1), (0.67, 1)),
        "north"         => ((0.67 ,1), (0.33,1)),
        "north_left"   => ((0.33 ,1), (0,1)),

        "west"          => ((0,0), (0,1)),

        "south_right"    => ((1, 0), (0.67, 0)),
        "south"         => ((0.67,0), (0.33,0)),
        "south_left"   => ((0.33,0), (0, 0))
    )

    return coords[orientation]
end

function initial_offset(::Rect5, orientation, direction)
    # Select where to start between the start and stop coordinates depending
    # on the direction selected
    starts = Dict(
        ( "north_left", Output ) => 0.55,
        ( "north", Output ) => 0.55,
        ( "north_right", Output ) => 0.55,
        ( "east" , Output ) => 0.55,
        ( "south_left", Input ) => 0.55,
        ( "south", Input ) => 0.55,
        ( "south_right", Input ) => 0.55,
        ( "west" , Input ) => 0.55,

        ( "north_left", Input ) => 0.55,
        ( "north", Input ) => 0.05,
        ( "north_right", Input ) => 0.55,
        ( "east" , Input ) => 0.05,
        ( "south_left", Output ) => 0.55,
        ( "south", Output ) => 0.05,
        ( "south_right", Output ) => 0.55,
        ( "west" , Output ) => 0.05,
    )

    return starts[(orientation, direction)]
end

# Extend `cartesian` to get the shifting of boxes correct
function cartesian(::Rect5, x, y)
    if mod(x, 4) == 1 || mod(x, 4) == 2
        return (1.5 * x, 1.5 * y + 0.75)
    end
    return (1.5*x, 1.5 * y)
end


text:
# # Threaded Assembly
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`threaded_assembly.ipynb`](@__NBVIEWER_ROOT_URL__/examples/threaded_assembly.ipynb)

# ## Example of a colored grid
#
# Creates a simple 2D grid and colors it.
# Save the example grid to a VTK file to show the coloring.
# No cells with the same color has any shared nodes (dofs).
# This means that it is safe to assemble in parallel as long as we only assemble
# one color at a time.

using Ferrite, SparseArrays

function create_example_2d_grid()
    grid = generate_grid(Quadrilateral, (10, 10), Vec{2}((0.0, 0.0)), Vec{2}((10.0, 10.0)))
    cell_colors, colors = Ferrite.create_coloring(grid)
    vtk_grid("colored", grid) do vtk
        Ferrite.vtk_cell_data_colors(vtk, grid, colors)
    end
end;

create_example_2d_grid()

# ![](coloring.png)

# ## Cantilever beam in 3D with threaded assembly
# We will now look at an example where we assemble the stiffness matrix using multiple
# threads. We set up a simple grid and create a coloring, then create a DofHandler,
# and define the material stiffness

# #### Grid for the beam
function create_colored_cantilever_grid(celltype, n)
    grid = generate_grid(celltype, (10*n, n, n), Vec{3}((0.0, 0.0, 0.0)), Vec{3}((10.0, 1.0, 1.0)))
    cell_colors, final_colors = Ferrite.create_coloring(grid)
    return grid, final_colors
end;

# #### DofHandler
function create_dofhandler(grid::Grid{dim}) where {dim}
    dh = DofHandler(grid)
    push!(dh, :u, dim) # Add a displacement field
    close!(dh)
end;

# ### Stiffness tensor for linear elasticity
function create_stiffness(::Val{dim}) where {dim}
    E = 200e9
    ν = 0.3
    λ = E*ν / ((1+ν) * (1 - 2ν))
    μ = E / (2(1+ν))
    δ(i,j) = i == j ? 1.0 : 0.0
    g(i,j,k,l) = λ*δ(i,j)*δ(k,l) + μ*(δ(i,k)*δ(j,l) + δ(i,l)*δ(j,k))
    C = SymmetricTensor{4, dim}(g);
    return C
end;

# ## Threaded data structures
#
# ScratchValues is a thread-local collection of data that each thread needs to own,
# since we need to be able to mutate the data in the threads independently
struct ScratchValues{T, CV <: CellValues, FV <: FaceValues, TT <: AbstractTensor, dim, Ti}
    Ke::Matrix{T}
    fe::Vector{T}
    cellvalues::CV
    facevalues::FV
    global_dofs::Vector{Int}
    ɛ::Vector{TT}
    coordinates::Vector{Vec{dim, T}}
    assembler::Ferrite.AssemblerSparsityPattern{T, Ti}
end;

# Each thread need its own CellValues and FaceValues (although, for this example we don't use
# the FaceValues)
function create_values(refshape, dim, order::Int)
    ## Interpolations and values
    interpolation_space = Lagrange{dim, refshape, 1}()
    quadrature_rule = QuadratureRule{dim, refshape}(order)
    face_quadrature_rule = QuadratureRule{dim-1, refshape}(order)
    cellvalues = [CellVectorValues(quadrature_rule, interpolation_space) for i in 1:Threads.nthreads()];
    facevalues = [FaceVectorValues(face_quadrature_rule, interpolation_space) for i in 1:Threads.nthreads()];
    return cellvalues, facevalues
end;

# Create a `ScratchValues` for each thread with the thread local data
function create_scratchvalues(K, f, dh::DofHandler{dim}) where {dim}
    nthreads = Threads.nthreads()
    assemblers = [start_assemble(K, f) for i in 1:nthreads]
    cellvalues, facevalues = create_values(RefCube, dim, 2)

    n_basefuncs = getnbasefunctions(cellvalues[1])
    global_dofs = [zeros(Int, ndofs_per_cell(dh)) for i in 1:nthreads]

    fes = [zeros(n_basefuncs) for i in 1:nthreads] # Local force vector
    Kes = [zeros(n_basefuncs, n_basefuncs) for i in 1:nthreads]

    ɛs = [[zero(SymmetricTensor{2, dim}) for i in 1:n_basefuncs] for i in 1:nthreads]

    coordinates = [[zero(Vec{dim}) for i in 1:length(dh.grid.cells[1].nodes)] for i in 1:nthreads]

    return [ScratchValues(Kes[i], fes[i], cellvalues[i], facevalues[i], global_dofs[i],
                         ɛs[i], coordinates[i], assemblers[i]) for i in 1:nthreads]
end;

# ## Threaded assemble

# The assembly function loops over each color and does a threaded assembly for that color
function doassemble(K::SparseMatrixCSC, colors, grid::Grid, dh::DofHandler, C::SymmetricTensor{4, dim}) where {dim}

    f = zeros(ndofs(dh))
    scratches = create_scratchvalues(K, f, dh)
    b = Vec{3}((0.0, 0.0, 0.0)) # Body force

    for color in colors
        ## Each color is safe to assemble threaded
        Threads.@threads for i in 1:length(color)
            assemble_cell!(scratches[Threads.threadid()], color[i], K, grid, dh, C, b)
        end
    end

    return K, f
end

# The cell assembly function is written the same way as if it was a single threaded example.
# The only difference is that we unpack the variables from our `scratch`.
function assemble_cell!(scratch::ScratchValues, cell::Int, K::SparseMatrixCSC,
                        grid::Grid, dh::DofHandler, C::SymmetricTensor{4, dim}, b::Vec{dim}) where {dim}

    ## Unpack our stuff from the scratch
    Ke, fe, cellvalues, facevalues, global_dofs, ɛ, coordinates, assembler =
         scratch.Ke, scratch.fe, scratch.cellvalues, scratch.facevalues,
         scratch.global_dofs, scratch.ɛ, scratch.coordinates, scratch.assembler

    fill!(Ke, 0)
    fill!(fe, 0)

    n_basefuncs = getnbasefunctions(cellvalues)

    ## Fill up the coordinates
    nodeids = grid.cells[cell].nodes
    for j in 1:length(coordinates)
        coordinates[j] = grid.nodes[nodeids[j]].x
    end

    reinit!(cellvalues, coordinates)

    for q_point in 1:getnquadpoints(cellvalues)
        for i in 1:n_basefuncs
            ɛ[i] = symmetric(shape_gradient(cellvalues, q_point, i))
        end
        dΩ = getdetJdV(cellvalues, q_point)
        for i in 1:n_basefuncs
            δu = shape_value(cellvalues, q_point, i)
            fe[i] += (δu ⋅ b) * dΩ
            ɛC = ɛ[i] ⊡ C
            for j in 1:n_basefuncs
                Ke[i, j] += (ɛC ⊡ ɛ[j]) * dΩ
            end
        end
    end

    celldofs!(global_dofs, dh, cell)
    assemble!(assembler, global_dofs, fe, Ke)
end;

function run_assemble()
    refshape = RefCube
    quadrature_order = 2
    dim = 3
    n = 20
    grid, colors = create_colored_cantilever_grid(Hexahedron, n);
    dh = create_dofhandler(grid);

    K = create_sparsity_pattern(dh);
    C = create_stiffness(Val{3}());
    ## compilation
    doassemble(K, colors, grid, dh, C);
    b = @elapsed @time K, f = doassemble(K, colors, grid, dh, C);
    return b
end

run_assemble()

# Running the code with different number of threads give the following runtimes:
# * 1 thread  2.46 seconds
# * 2 threads 1.19 seconds
# * 3 threads 0.83 seconds
# * 4 threads 0.75 seconds

#md # ## [Plain Program](@id threaded_assembly-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [threaded_assembly.jl](threaded_assembly.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using DataStructures
using Parameters

@with_kw struct Cell
    reactants::Array{Int64}
    rates::Array{Float64}
    A::Array{Float64}=zeros(1,34)
    Tk::Array{Float64}=zeros(1,34)
    Pk::Array{Float64}=zeros(1,34)
    pq::PriorityQueue{Int,Float64}=PriorityQueue{Int,Float64}()
    pq_delayed::Array{Array{Float64}}=[[Inf],[Inf],[Inf],[Inf],[Inf],[Inf],[Inf],[Inf]]
    pos::Array{Int}=[1,1,1,1,1,1,1,1]
    levels::Vector{Int}=Vector{Int}()
    Time::Vector{Float64}=Vector{Float64}()
    t::Vector{Float64}=[0.0]
end

function store!(c::Cell)
    push!(c.levels, c.reactants[1])
    push!(c.Time, c.t[1])
end

function update_ndelayed!(c::Cell, n::Int)
    c.pos[n] += 1
end

struct Model
    cells::Vector{Cell}
    T::Vector{Float64}
    propRules::Array{Function}
    stoichiometry::Array{Function}
end

max_stars_repo_name: syedtaz/project-cait

max_issues_repo_name: syedtaz/project-cait

max_forks_repo_name: syedtaz/project-cait


text:
module TestGraphvizWiringDiagrams

using Test
import JSON

using Catlab.Theories, Catlab.WiringDiagrams, Catlab.Graphics
using Catlab.Programs: @relation
import Catlab.Graphics: Graphviz
using Catlab.Graphics.WiringDiagramLayouts: position, normal

function stmts(graph::Graphviz.Graph, type::Type)
  [ stmt for stmt in graph.stmts if stmt isa type ]
end
function stmts(graph::Graphviz.Graph, type::Type, attr::Symbol)
  [ stmt.attrs[attr] for stmt in graph.stmts
    if stmt isa type && haskey(stmt.attrs, attr) ]
end

# Directed drawing
##################

f = singleton_diagram(Box(:f, [:A], [:B]))
g = singleton_diagram(Box(:g, [:B], [:A]))

# Simple wiring diagrams, with default settings.
graph = to_graphviz(f)
@test graph isa Graphviz.Graph && graph.directed
@test stmts(graph, Graphviz.Node, :comment) == ["f"]
@test stmts(graph, Graphviz.Edge, :comment) == ["A","B"]
@test stmts(graph, Graphviz.Edge, :label) == []
@test stmts(graph, Graphviz.Edge, :xlabel) == []
@test length(stmts(graph, Graphviz.Subgraph)) == 2

graph = to_graphviz(singleton_diagram(Box(:h, Symbol[], [:A])))
@test stmts(graph, Graphviz.Node, :comment) == ["h"]
@test length(stmts(graph, Graphviz.Subgraph)) == 1

graph = to_graphviz(singleton_diagram(Box(:h, Symbol[], Symbol[])))
@test stmts(graph, Graphviz.Node, :comment) == ["h"]
@test length(stmts(graph, Graphviz.Subgraph)) == 0

graph = to_graphviz(compose(f,g))
@test stmts(graph, Graphviz.Node, :comment) == ["f","g"]

# Layout orientation.
for orientation in instances(LayoutOrientation)
  graph = to_graphviz(compose(f,g); orientation=orientation)
  @test stmts(graph, Graphviz.Node, :comment) == ["f","g"]
end

# Junction nodes.
graph = to_graphviz(add_junctions!(compose(f, mcopy(Ports([:B])))))
@test stmts(graph, Graphviz.Node, :comment) == ["f","junction"]

# Edge labels.
graph = to_graphviz(f; labels=true)
@test stmts(graph, Graphviz.Edge, :label) == ["A","B"]
@test stmts(graph, Graphviz.Edge, :xlabel) == []

graph = to_graphviz(f; labels=true, label_attr=:xlabel)
@test stmts(graph, Graphviz.Edge, :label) == []
@test stmts(graph, Graphviz.Edge, :xlabel) == ["A","B"]

# Anchoring of outer ports.
graph = to_graphviz(otimes(f,g); anchor_outer_ports=true)
@test stmts(graph, Graphviz.Node, :comment) == ["f","g"]
graph = to_graphviz(otimes(f,g); anchor_outer_ports=false)
@test stmts(graph, Graphviz.Node, :comment) == ["f","g"]

# Undirected drawing
####################

d = singleton_diagram(UndirectedWiringDiagram, 2)
graph = to_graphviz(d)
@test stmts(graph, Graphviz.Node, :id) ==
  ["box1", "outer1", "outer2", "junction1", "junction2"]
@test length(stmts(graph, Graphviz.Edge)) == 4

graph = to_graphviz(d, implicit_junctions=true)
@test stmts(graph, Graphviz.Node, :id) == ["box1", "outer1", "outer2"]
@test length(stmts(graph, Graphviz.Edge)) == 2

d = @relation ((x,z) where (x,y,z)) -> (R(x,y); S(y,z))
graph = to_graphviz(d)
@test stmts(graph, Graphviz.Node, :id) ==
  ["box1", "box2", "outer1", "outer2", "junction1", "junction2", "junction3"]
@test length(stmts(graph, Graphviz.Edge)) == 6

graph = to_graphviz(d, implicit_junctions=true)
@test stmts(graph, Graphviz.Node, :id) == ["box1", "box2", "outer1", "outer2"]
@test length(stmts(graph, Graphviz.Edge)) == 3

# Box, port, and junction labels.
graph = to_graphviz(d, box_labels=:name, junction_labels=:variable,
                    port_labels=false)
@test filter(!isempty, stmts(graph, Graphviz.Node, :label)) == ["R", "S"]
@test stmts(graph, Graphviz.Node, :xlabel) == ["x", "y", "z"]
@test stmts(graph, Graphviz.Edge, :taillabel) == []

graph = to_graphviz(d, box_labels=true, junction_labels=true, port_labels=true)
@test filter(!isempty, stmts(graph, Graphviz.Node, :label)) == ["1", "2"]
@test stmts(graph, Graphviz.Node, :xlabel) == ["1", "2", "3"]
@test stmts(graph, Graphviz.Edge, :taillabel) == ["1","1","2","1","2","2"]

# Directed layout
#################

diagram = include(joinpath("data", "graphviz_wiring_diagram.jl"))
doc = open(JSON.parse,
  joinpath(@__DIR__, "data", "graphviz_wiring_diagram.json"), "r")
graph = parse_graphviz(doc)
layout = graphviz_layout(diagram, graph)

# Is original data preserved?
values(xs) = map(x -> x.value, xs)
@test values(input_ports(layout)) == input_ports(diagram)
@test values(output_ports(layout)) == output_ports(diagram)
@test values(values(boxes(layout))) == values(boxes(diagram))

# Basic geometry.
positions = map(position, boxes(layout))
@test positions[1][1] < positions[2][1]
@test positions[1][2] < positions[3][2]

end


text:
using Revise
using Nexus
using Test
using Dictionaries
using StatsBase

@testset "Nexus.jl" begin
    @testset "graph.jl" begin
        include("graph.jl")
    end
    @testset "arc.jl" begin
        include("arc.jl")
    end
    @testset "path.jl" begin
        include("path.jl")
    end
    @testset "node.jl" begin
        include("node.jl")
    end
end


text:
using SenseHat, Colors, Dates, DataFrames, UnicodePlots
########################
const LED = led_matrix()

function julia_logo()
	LED[:] = SenseHat.JULIA_LOGO
	sleep(3)
	led_clear()
end
########################
function usgs_logo()
	 USGS_LOGO = begin
		W = RGB565(.4,.4,.4)
		G = RGB565(.1,.296, .09)

	permutedims([G G G G G G W G;
        	 G W W G G W G W;
                 W G G W W G W G;
                 G W W G G W G W;
                 W G G W W G W G;
                 G W W G G W G G;
                 W G G W W G G G;
                 G G G G G G G G], (2,1)) # transpose to correct orientation
	end

	LED[:] = USGS_LOGO
	sleep(3)
	led_clear()
end
##############################
df = DataFrame(DateTime = DateTime[],
			   Temperature = Float64[],
			   Humidity = Float64[],
			   Pressure = Float64[])

function get_data()

	d = Dates.now()
	t = round((temperature() * 9/5) + 32, digits=2)
	h = round(humidity(), digits=2)
	p = round(pressure()*0.750062, digits=2)

	println("Temp: $t C", "   ", "Humidity: $h %", "   ", "Pressure: $p mm/hg")
	push!(df::DataFrame, [d t p h])
	#sleep(1)
end

function print_data()

	t = round((temperature() * 9/5) + 32, digits=2)
	show_message("Temp: $t", 0.1)

	h = round(humidity(), digits=2)
	show_message("Humid: $h", 0.1)

	p = round(pressure(), digits=2)
	show_message("Pres: $p",0.1)
end

i = 0
while true
	#julia_logo()
	#usgs_logo()
	#print_data()
	get_data()

	if i % 10 === 0
		display(lineplot(df[:DateTime],df[:Temperature], color=:red))
		println("")
	end

	sleep(1)

	if i >= 100000
		break
	end
	global i += 1
end

max_stars_repo_name: atott/Sense-Julia

max_issues_repo_name: atott/Sense-Julia

max_forks_repo_name: atott/Sense-Julia


text:
using AssignOnce
using Test

struct InvocationCounter{F}
    f::F
    counter::Ref{Int}
end

InvocationCounter(f) = InvocationCounter(f, Ref(0))

function (c::InvocationCounter)(args...; kwargs...)
    c.counter[] += 1
    c.f(args...; kwargs...)
end

@testset "check test purity 1" begin
    y = 3
end

@testset "check test purity 2" begin
    @test !isdefined(@__MODULE__, :y)
end

@testset "InvocationCounter" begin
    f = InvocationCounter(()->3)
    @test f.counter[] == 0
    @test f() == 3
    @test f.counter[] == 1
    @test f() == 3
    @test f.counter[] == 2
end

# TODO maybe generate a new module, so there is a new global scope, for each test set
@testset "does assignment" begin
    @assert !isdefined(@__MODULE__, :x)
    @assignonce x = 3
    @test isdefined(@__MODULE__, :x)
    @test x == 3
end

# `x` is used up.
# TODO, same thing. would be nice to have a new module scope so we can define `ff` in there

ff = InvocationCounter(()->3)
@testset "does assignment only once" begin
    @assert !isdefined(@__MODULE__, :y)
    @assert ff.counter[] == 0
    @assignonce y = ff()
    @test ff.counter[] == 1
    @assignonce y = ff()
    @test ff.counter[] == 1
    @test isdefined(@__MODULE__, :y)
    @test y == 3
end


text:
################################################################################
## Defines KKT system solvers
## KKT system solvers modify the orginal linear system to make it possible to use a linear system solver.
## For example,
## - schur.jl does a primal schur complement which can then be solved with a cholseky linear system solver.
## - symmetric.jl forms a symmetric system which can be solved using LDL
################################################################################

@compat abstract type abstract_KKT_system_solver
    # list variables
end
# The following functions should be defined for any KKT_solver <: abstract_KKT_system_solver:
# - form_system!(kkt_solver::abstract_schur_solver, iter::Class_iterate, timer::class_advanced_timer)
# - factor_implementation!(kkt_solver::KKT_solver, timer::class_advanced_timer)
# - compute_direction_implementation!(kkt_solver::KKT_solver, timer::class_advanced_timer)
# - update_delta_vecs!(kkt_solver::KKT_solver, delta_x_vec::Array{Float64,1}, delta_s_vec::Array{Float64,1}, timer::class_advanced_timer)

@compat abstract type abstract_schur_solver <: abstract_KKT_system_solver end

function initialize!(kkt_solver::abstract_KKT_system_solver, intial_it::Class_iterate)
    # call this before running using the kkt_solver
    initialize!(kkt_solver.ls_solver)
    kkt_solver.dir = zero_point(dim(intial_it),ncon(intial_it))
end

function predicted_lag_change(kkt_solver::abstract_KKT_system_solver)
    #################################################################################
    # INPUT:
    # kkt_solver with a direction and delta_x_vec representing how we perturbed the system to make the inertia correct.
    # OUTPUT:
    # vector denoting how do we expect gradient of the Lagragian to change
    #################################################################################

    fi = kkt_solver.factor_it

    dir_x = kkt_solver.dir.x
    dir_y = kkt_solver.dir.y

    tmp1 = eval_jac_prod(fi, dir_x)
    tmp2 = kkt_solver.delta_s_vec .* tmp1
    tmp3 = eval_jac_T_prod(fi, tmp2)
    delta_err = kkt_solver.delta_x_vec .* dir_x + tmp3
    J_err = eval_jac_T_prod(fi, dir_y)
    H_err = hess_product(kkt_solver.factor_it, dir_x)
    return delta_err + H_err - J_err
end

type Class_kkt_error
    # when we compute the directions how much error is there in the linear system.
    error_D::Float64 # norm of error in dual feasibility
    error_P::Float64 # norm error of error in primal feasibility
    error_mu::Float64 # norm error of error in complementarity
    overall::Float64 # total error
    rhs_norm::Float64 # norm of rhs of linear system
    ratio::Float64 # ratio of total error to norm of rhs. If this is less than zero then the direction is improving things.

    function Class_kkt_error(error_D::Float64,error_P::Float64,error_mu::Float64,overall::Float64, rhs_norm::Float64, ratio::Float64)
          return new(error_D,error_P,error_mu,overall,rhs_norm, ratio)
    end

    function Class_kkt_error()
        return new(0.0,0.0,0.0,0.0)
    end
end

function update_kkt_error!(ss::abstract_KKT_system_solver, p::Float64, timer::class_advanced_timer)
    factor_iter = ss.factor_it;
    rhs = ss.rhs;
    dir = ss.dir;

    dual_scaling_can = 1.0 #dual_scale(candidate)

    start_advanced_timer(timer, "update_kkt_error/compute/D")
    error_D = (predicted_lag_change(ss) - rhs.dual_r) * dual_scaling_can;
    pause_advanced_timer(timer, "update_kkt_error/compute/D")

    start_advanced_timer(timer, "update_kkt_error/compute/P")
    error_P = eval_jac_prod(factor_iter,dir.x) - dir.s - rhs.primal_r;
    pause_advanced_timer(timer, "update_kkt_error/compute/P")

    start_advanced_timer(timer, "update_kkt_error/compute/mu")
    error_mu = get_s(factor_iter) .* dir.y + get_y(factor_iter) .* dir.s  - rhs.comp_r;
    pause_advanced_timer(timer, "update_kkt_error/compute/mu")

    #@show norm(error_D), norm(error_P), norm(error_mu)
    overall = norm([error_D; error_P; error_mu], p)

    dual_scaling_org = 1.0 #dual_scale(iter)
    rhs_norm = norm([rhs.dual_r * dual_scaling_org; rhs.primal_r; rhs.comp_r], p)

    ratio = overall / rhs_norm;
    #@show ratio

    ss.kkt_err_norm = Class_kkt_error(norm(error_D,p), norm(error_P,p), norm(error_mu,p), overall, rhs_norm, ratio)
end

function factor!(kkt_solver::abstract_KKT_system_solver, delta_x::Float64, delta_s::Float64, timer::class_advanced_timer)
    # I don't think this is used anymore
    update_delta!(kkt_solver, delta_x, delta_s, timer)
    factor!(kkt_solver, timer)
end

function factor!(kkt_solver::abstract_KKT_system_solver, delta_x::Float64, timer::class_advanced_timer)
    # factorize the linear system for a particular perturbation to the Hessian of I * delta_x.
    factor!(kkt_solver, delta_x, 0.0, timer)
end

function update_delta!(kkt_solver::abstract_KKT_system_solver, delta_x::Float64, delta_s::Float64, timer::class_advanced_timer)
  delta_x_vec = delta_x * ones(dim(kkt_solver.factor_it))
  delta_s_vec = delta_s * get_s(kkt_solver.factor_it).^(-2)
  update_delta_vecs!(kkt_solver, delta_x_vec, delta_s_vec, timer)
end

function factor_at_approx_min_eigenvalue!(kkt_solver::abstract_KKT_system_solver, iter::Class_iterate)
    # Delete: I don't think this is used

    start_advanced_timer("KKT/factor");
    form_system!(kkt_solver, iter)

    max_it = 100

    delta_min = 1e-8

    inertia = factor!(kkt_solver, 0.0)

    j = 1;
    if inertia == 1
      set_delta(iter, 0.0)
    else
      set_delta(iter, max(delta_min, get_delta(iter)))
      for i = 1:max_it
        inertia = factor!(kkt_solver, get_delta(iter))

        if inertia == 1
            set_delta(iter, get_delta(iter) / 20.0 )
        else
          break
        end
      end

      for j = 1:max_it
        inertia = factor!(kkt_solver, get_delta(iter))

        if inertia == 1
            break
        else
           set_delta(iter, get_delta(iter) * 3.0 )
        end
      end
    end

    if j == max_it
      @show get_delta(iter)
      error("delta too large!")
    end

    pause_advanced_timer("KKT/factor");
end

function kkt_associate_rhs!(kkt_solver::abstract_KKT_system_solver, iter::Class_iterate, eta::Class_reduction_factors, timer::class_advanced_timer)
    # create right hand side for linear system and attach it to iter.
    start_advanced_timer(timer, "KKT/rhs");

    kkt_solver.rhs = System_rhs(iter, eta)
    kkt_solver.dir.mu = -(1.0 - eta.mu) * get_mu(iter)
    kkt_solver.dir.primal_scale = -(1.0 - eta.P) * iter.point.primal_scale

    pause_advanced_timer(timer, "KKT/rhs");
end

function compute_direction!(kkt_solver::abstract_KKT_system_solver, timer::class_advanced_timer)
    # compute a direction and store it in the kkt_solver object
    start_advanced_timer(timer, "KKT/compute_direction");
    if kkt_solver.ready != :factored
        error("kkt solver not ready to compute direction!")
    end

    compute_direction_implementation!(kkt_solver, timer)
    check_for_nan(kkt_solver.dir)
    pause_advanced_timer(timer, "KKT/compute_direction");
end

function factor!(kkt_solver::abstract_KKT_system_solver, timer::class_advanced_timer)
    # factorize the linear system
    # stores the factorize inside the kkt_solver
    start_advanced_timer(timer, "KKT/factor");

    if kkt_solver.ready != :delta_updated
        error("kkt solver not ready to factor kkt_solver.ready = $(kkt_solver.ready) != :delta_updated")
    else
        kkt_solver.ready = :factored
    end
    inertia = factor_implementation!(kkt_solver, timer)
    pause_advanced_timer(timer, "KKT/factor");

    return inertia # is the inertia correct?
end

function compute_eigenvector!(kkt_solver::abstract_KKT_system_solver, iter::Class_iterate, timer::class_advanced_timer)
    # Delete: I think this is old code not really used.
    start_advanced_timer(timer, "KKT/direction/eig");

    approx_eigvec = randn(dim(iter));

    for i = 1:20
      rhs = System_rhs(iter)
      rhs.dual_r = approx_eigvec
      kkt_solver.rhs = rhs
      kkt_solver.dir.mu = 0.0

      compute_direction!(kkt_solver, timer)
      approx_eigvec = kkt_solver.dir.x / norm(kkt_solver.dir.x,2)
    end

    lambda = norm(kkt_solver.dir.x,2)
    eig_est = (1 - kkt_solver.delta_x_vec[1] * lambda) / lambda

    kkt_solver.dir = scale_direction(kkt_solver.dir, 1.0 / lambda)

    pause_advanced_timer(timer, "KKT/direction/eig");

    return eig_est
end

function pick_KKT_solver(pars::Class_parameters)
    # selects the KKT_solver based on the parameters
  kkt_solver_type = pars.kkt.kkt_solver_type
  linear_solver_type = pars.kkt.linear_solver_type
  safe = pars.kkt.linear_solver_safe_mode
  recycle = pars.kkt.linear_solver_recycle

  if kkt_solver_type == :symmetric
        my_kkt_solver = Symmetric_KKT_solver()
        if linear_solver_type == :julia
          my_kkt_solver.ls_solver = linear_solver_JULIA(:symmetric, safe, recycle)
        elseif linear_solver_type == :mumps
          my_kkt_solver.ls_solver = linear_solver_MUMPS(:symmetric, safe, recycle)
        elseif linear_solver_type == :HSL
            my_kkt_solver.ls_solver = linear_solver_HSL(:symmetric, safe, recycle)
        else
            error("pick a valid solver!")
        end
  elseif kkt_solver_type == :clever_symmetric
      my_kkt_solver = Clever_Symmetric_KKT_solver()
      if linear_solver_type == :julia
        my_kkt_solver.ls_solver = linear_solver_JULIA(:symmetric, safe, recycle)
      elseif linear_solver_type == :mumps
        my_kkt_solver.ls_solver = linear_solver_MUMPS(:symmetric, safe, recycle)
      elseif linear_solver_type == :HSL
          my_kkt_solver.ls_solver = linear_solver_HSL(:symmetric, safe, recycle)
      else
          error("pick a valid solver!")
      end
  elseif kkt_solver_type == :schur
    my_kkt_solver = Schur_KKT_solver()
    if linear_solver_type == :julia
      my_kkt_solver.ls_solver = linear_solver_JULIA(:definite, safe, recycle)
    elseif linear_solver_type == :mumps
      my_kkt_solver.ls_solver = linear_solver_MUMPS(:definite, safe, recycle)
    else
      error("pick a valid solver!")
    end
  elseif kkt_solver_type == :schur_direct
    my_kkt_solver = Schur_KKT_solver_direct()
    if linear_solver_type == :julia
      my_kkt_solver.ls_solver = linear_solver_JULIA(:definite, safe, recycle)
    elseif linear_solver_type == :mumps
      my_kkt_solver.ls_solver = linear_solver_MUMPS(:definite, safe, recycle)
  else
    error("pick a valid solver!")
  end
  else
      error("pick a solver!")
  end

  my_kkt_solver.kkt_err_norm = Class_kkt_error()
  my_kkt_solver.pars = pars

  return my_kkt_solver
end

## helper functions

function diag_min(kkt_solver::abstract_KKT_system_solver)
    # find the smallest element of kkt_solver.schur_diag
    return minimum(kkt_solver.schur_diag)
end

function compute_schur_diag(iter::Class_iterate)
    # compute the diagonal of the primal schur complement matrix
    schur_scaling = get_y(iter) ./ get_s(iter);
    return diag(get_lag_hess(iter)) + eval_diag_J_T_J(iter,schur_scaling)
end


text:
"""
    abstract type AbstractReferenceShape{N}

A reference domain/shape in `ℜᴺ`.

Used mostly for defining more complex shapes as transformations mapping an
`AbstractReferenceShape` into some region of `ℜᴹ`.

See e.g. [`ReferenceLine`](@ref) or [`ReferenceTriangle`](@ref) for some
examples of concrete subtypes.
"""
abstract type AbstractReferenceShape{N} end

ambient_dimension(::SType{<:AbstractReferenceShape{N}}) where {N}    = N
geometric_dimension(::SType{<:AbstractReferenceShape{N}}) where {N}  = N
dimension(::SType{<:AbstractReferenceShape{N}}) where {N}  = N

"""
    struct ReferencePoint

Singleton type representing a reference zero-dimensional entitty (i.e. a point).
"""
struct ReferencePoint <: AbstractReferenceShape{0}
end
number_of_nodes(::SType{ReferencePoint}) = 1

"""
    struct ReferenceLine

Singleton type representing the `[0,1]` segment.
"""
struct ReferenceLine <: AbstractReferenceShape{1}
end
Base.in(x,::ReferenceLine)    = 0 ≤ x[1] ≤ 1
center(::Type{ReferenceLine}) = 0.5
center(::ReferenceLine)       = 0.5
number_of_nodes(::Type{ReferenceLine}) = 2
number_of_nodes(::ReferenceLine)       = 2

vertices(ln::ReferenceLine) = SVector(0), SVector(1)

"""
    struct ReferenceTriangle

Singleton type representing the triangle with vertices `(0,0),(0,1),(1,0)`
"""
struct ReferenceTriangle <: AbstractReferenceShape{2}
end
Base.in(x,::ReferenceTriangle) = 0 ≤ x[1] ≤ 1 && 0 ≤ x[2] ≤ 1 - x[1]
number_of_nodes(::Type{ReferenceTriangle}) = 3
number_of_nodes(::ReferenceTriangle)       = 3

vertices(::ReferenceTriangle) = SVector(0,0), SVector(1,0), SVector(0,1)

"""
    struct ReferenceHyperCube{N}

Singleton type representing the axis-aligned hypercube in `N` dimensions with
the lower corner at the origin and the upper-corner at `(1,1,…,1)`.
"""
struct ReferenceHyperCube{N} <: AbstractReferenceShape{N}
end

"""
    const ReferenceSquare = ReferenceHyperCube{2}

Singleton type representing the square with vertices `(0,0),(0,1),(1,1),(1,0)`
"""
const ReferenceSquare = ReferenceHyperCube{2}
Base.in(x,::ReferenceSquare)    = 0 ≤ x[1] ≤ 1 && 0 ≤ x[2] ≤ 1
center(::Type{ReferenceSquare}) = SVector(0.5,0.5)
center(::ReferenceSquare)       = SVector(0.5,0.5)
number_of_nodes(::Type{ReferenceSquare}) = 4
number_of_nodes(::ReferenceSquare)       = 4

vertices(sq::ReferenceSquare) = SVector(0,0), SVector(1,0), SVector(1,1), SVector(0,1)

"""
    struct ReferenceTetrahedron

Singleton type representing the tetrahedron with vertices `(0,0,0),(0,0,1),(0,1,0),(1,0,0)`
"""
struct ReferenceTetrahedron <: AbstractReferenceShape{3}
end
Base.in(x,::ReferenceTetrahedron) = 0 ≤ x[1] ≤ 1 &&
                                    0 ≤ x[2] ≤ 1 - x[1] &&
                                    0 ≤ x[3] ≤ 1 - x[1] - x[2]
number_of_nodes(::Type{ReferenceTetrahedron}) = 4
number_of_nodes(::ReferenceTetrahedron)       = 4

# TODO: generalize structs above to `ReferenceSimplex{N}` and
# `ReferenceCuboid{N}`

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
using Test, ColorTypes, NES

@testset "Load Cartridge Test" begin
  c = loadgame("test/Super_Tilt_Bro.nes")
  actual = screen(c)
  expected = open("test/blank_screen.bin") do f
    reinterpret(RGB24, reinterpret(UInt32, reshape(read(f), 4 * 240, 256)))
  end
  @test actual == expected
end

@testset "Single CPU Instruction Test" begin
  c = loadgame("test/Super_Tilt_Bro.nes")
  step!(c)
  @test c.cpu.cycles == 2
  @test c.cpu.PC == 0x8036
  @test c.cpu.SP == 0xFD
  @test c.cpu.A == 0x00
  @test c.cpu.X == 0x00
  @test c.cpu.Y == 0x00
  @test c.cpu.C == 0x00
  @test c.cpu.Z == 0x00
  @test c.cpu.I == 0x01
  @test c.cpu.D == 0x00
  @test c.cpu.B == 0x00
  @test c.cpu.U == 0x01
  @test c.cpu.V == 0x00
  @test c.cpu.N == 0x00
  @test Int(c.cpu.interrupt) == 0
  @test c.cpu.stall == 0
end

@testset "Single CPU Frame Tests" begin
  c = loadgame("test/Super_Tilt_Bro.nes")
  stepframe!(c)
  @test c.cpu.cycles == 0x0000000000000955
  @test c.cpu.PC == 0x8050
  @test c.cpu.SP == 0xFF
  @test c.cpu.A == 0xFE
  @test c.cpu.X == 0x31
  @test c.cpu.Y == 0x00
  @test c.cpu.C == 0x00
  @test c.cpu.Z == 0x00
  @test c.cpu.I == 0x01
  @test c.cpu.D == 0x00
  @test c.cpu.B == 0x00
  @test c.cpu.U == 0x01
  @test c.cpu.V == 0x00
  @test c.cpu.N == 0x00
  @test Int(c.cpu.interrupt) == 0
  @test c.cpu.stall == 0
end

@testset "Screen Rendering Tests" begin
  c = loadgame("test/Super_Tilt_Bro.nes")
  stepframes!(c, 120)
  actual = screen(c)
  expected = open("test/screen.bin") do f
    reinterpret(RGB24, reinterpret(UInt32, reshape(read(f), 4 * 240, 256)))
  end
  @test actual == expected
end

@testset "Controller Input Tests" begin
  c = loadgame("test/Super_Tilt_Bro.nes")
  inputs = [(0x00, 0x60), (0x08, 0x05), (0x00, 0x60), (0x08, 0x05),
            (0x00, 0x60), (0x08, 0x05), (0x00, 0x60), (0x08, 0x05),
            (0x00, 0x50), (0x80, 0x15), (0x81, 0x05)]
  stepframes!(c, inputs)
  actual = screen(c)
  expected = open("test/inputs.bin") do f
    reinterpret(RGB24, reinterpret(UInt32, reshape(read(f), 4 * 240, 256)))
  end
  @test actual == expected
end