text:
export Shamanskii_MA57_BFGS

function Shamanskii_MA57_BFGS(nlp_stop, PData::PDataFact, dₙ::Vector, gt::Vector)
    # printstyled("dans Shamanskii_MA57_BFGS ↓ \n", color = :blue)
    nlp_at_x = nlp_stop.current_state
    ∇fₖ = nlp_at_x.gx
    x = nlp_at_x.x
    xdemi = x .+ dₙ
    ∇fₖ₊₁ = grad(nlp_stop.pb, xdemi)
    gtemp = copy(∇fₖ₊₁)
    Bₖ = nlp_at_x.Hx
    # @show Bₖ
    # @show PData.L
    # @show PData.D
    # @show PData.s
    # @show PData.pp
    # L = convert(SparseMatrixCSC{Float64, Int64}, PData.L)
    # D = convert(SparseMatrixCSC{Float64, Int64}, PData.D)
    # S = sparse(Diagonal(PData.s))
    # P = SparseMatrixCSC(I, length(gtemp), length(gtemp))[PData.pp, :]
    # Bₖ  = Matrix(inv(Matrix(S)) * P * L * D * L' * P' * inv(Matrix(S)))
    ϵ2 = sqrt(eps(eltype(gt)))
    Γ = max.(abs.(PData.Δ), ϵ2)
    # αₖ = 1.0
    sₖ = dₙ

    yₖ = ∇fₖ₊₁ .- ∇fₖ
    Yₖ = yₖ * yₖ'
    yts = yₖ' * sₖ
    Yₖ = Yₖ ./ yts

    Bs = Bₖ * sₖ
    BssBₖ = Bs * Bs'
    sBs = Bs' * sₖ
    BssBₖ = (1.0 ./ sBs) .* BssBₖ

    Bₖ₊₁ = Bₖ .+ Yₖ .- BssBₖ
    B57 = convert(SparseMatrixCSC{Float64,Int64}, Bₖ₊₁)
    M = Ma57(B57)
    ma57_factorize(M)
    (L2, D2, s2, pp2) = ma57_get_factors(M)
    # @show L2
    # @show D2
    # @show s2
    # @show pp2
    L2 = SparseMatrixCSC{Float64,Int64}(L2)
    vD1 = diag(D2)       # create internal representation for block diagonal D
    vD2 = diag(D2, 1)     #
    vQ1 = ones(length(vD1))       # vector representation of orthogonal matrix Q
    vQ2 = zeros(length(vD2))      #
    vQ2m = zeros(length(vD2))     #
    veig = copy(vD1)      # vector of eigenvalues of D, initialized to diagonal of D
    # if D diagonal, nothing more will be computed

    i = 1
    while i < length(vD1)
        if vD2[i] == 0.0
            i += 1
        else
            mA = [vD1[i] vD2[i]; vD2[i] vD1[i+1]] #  2X2 submatrix
            # DiagmA, Qma = eig(mA)                   #  spectral decomposition of mA
            X = eigen(mA)
            DiagmA = X.values
            Qma = X.vectors
            veig[i] = DiagmA[1]
            vQ1[i] = Qma[1, 1]
            vQ2[i] = Qma[1, 2]
            vQ2m[i] = Qma[2, 1]
            vQ1[i+1] = Qma[2, 2]
            veig[i+1] = DiagmA[2]
            i += 2
        end
    end

    Q2 = sparse(SparseArrays.spdiagm(0 => vQ1, 1 => vQ2m, -1 => vQ2))
    dt = L2 \ gtemp[pp2]
    dc = L2' \ (Q2 * (Q2' * dt ./ Γ))
    d = -dc[invperm(pp2)]

    dₕₒ = dₙ .+ d

    return dₕₒ
end


text:
abstract type SolutionState end

struct AntEnding <: SolutionState end

struct NCEnding <: SolutionState end

struct InProgress <: SolutionState end

struct NoChangeEnding <: SolutionState end

mutable struct Solution

    path::Vector{City}
    length::Real

    usedTime::Float64
    usedSpace::Float64

    state:: SolutionState

end # struct

function Solution()
    return Solution(Vector{City}(), Inf, 0., 0., InProgress())
end

#Taille des chemins proportionnels aux Phéromones

ext: jl


text:
using CombUtils

function scale_numerator(numerator, denominator)
    while numerator < denominator
        numerator *= 10
    end
    return numerator
end

function reciprocal_cycle_len{T<:Integer}(n::T)
    num = scale_numerator(10, n)
    while num < n
        num *= 10
    end

    while num > n
        if num < n
            num *= 10
        else
            num %= n
        end
    end
end

function main()
    seq_length = 0
    num = 0

    for i in 1000:-1:1
        if seq_length >= i
            break
        end

        found_remainders = zeros(i)
        value = 1
        position = 0

        while value != 0 && found_remainders[value] == 0
            found_remainders[value] = position;
            value *= 10
            value %= i
            position += 1
        end

        rem_at_value = value == 0 ? 0 : found_remainders[value]

        if position - rem_at_value > seq_length
            num = i
            seq_length = position - rem_at_value
        end
    end
    @printf("num = %d, seq_length = %d\n", num, seq_length)
end

main()


text:
@testset "Differentiating LP; checking gradients for non-active contraints" begin
    # Issue #40 from Gurobi.jl
    # min  x
    # s.t. x >= 0
    #      x >= 3
    qp_test_with_solutions(
        GLPK.Optimizer,
        q = [1.0],
        G = -ones(2, 1),
        h = [0.0, -3.0],
        dzb = -ones(1),
        dhf = [0.0, 1.0],
        # Expected solutions
        z = [3.0],
        λ = [0.0, 1.0],
        dzf = -ones(1),
        dλb = zeros(2),
        dhb = [0.0, 1.0],
        ∇zb = zeros(1),
        ∇λb = [0.0, -1.0],
        dλf = zeros(2),
    )
end


text:
# Load packages
using Dojo
using Plots
using Random
using MeshCat

# Open visualizer
vis= Visualizer()
open(vis)

# Include new files
include(joinpath(module_dir(), "env", "sphere", "deps", "texture.jl"))
include( "../utils.jl")

mech = get_mechanism(:sphere, timestep=0.05, g=-9.81, radius=0.5, friction_coefficient=0.1);
initialize!(mech, :sphere, x=[0,0,0.3], v=[0,0.5,0.], ω=[10,0,0.])
storage = simulate!(mech, 0.5, record=true, verbose=true,
    opts=SolverOptions(btol=1e-6, rtol=1e-6))
visualize(mech, storage, vis=vis)
sphere_texture!(vis, mech)

################################################################################
# Generate & Save Dataset
################################################################################
init_kwargs = Dict(:xlims => [[0,0,0], [1,1,0.2]],
				   :vlims => [-1ones(3), 1ones(3)],
				   :ωlims => [-5ones(3), 5ones(3)])
mech_kwargs = Dict(:friction_coefficient => 0.1, :radius => 0.5)
generate_dataset(:sphere, H=0.75, N=15,
	opts=SolverOptions(btol=3e-4, rtol=3e-4),
	init_kwargs=init_kwargs,
	mech_kwargs=mech_kwargs)

################################################################################
# Load Dataset
################################################################################
params0, trajs0, pairs0 = open_dataset(:sphere; N=15, mech_kwargs...)

data0 = params0[:data]

################################################################################
# CLEAN Optimization Objective: Evaluation & Gradient
################################################################################

clean_loss(:sphere, pairs0, data0, opts=SolverOptions(btol=3e-4, rtol=3e-4))

[clean_loss(:sphere, pairs0, [0.1, 0,0, 0.5+i, 0,0,0], opts=SolverOptions(btol=3e-4, rtol=3e-4))
	for i in Vector(-0.10:0.01:0.1)]
[clean_loss(:sphere, pairs0, [0.1+i, 0,0, 0.5, 0,0,0], opts=SolverOptions(btol=3e-4, rtol=3e-4))
	for i in Vector(-0.10:0.01:0.1)]

plot(hcat([p[1][1:3] for p in pairs0]...)')
plot(hcat([p[1][4:6] for p in pairs0]...)')
plot(hcat([p[1][7:10] for p in pairs0]...)')
plot(hcat([p[1][11:13] for p in pairs0]...)')

################################################################################
# CLEAN Optimization Algorithm: L-BFGS
################################################################################

# Solution for bad cost landscape
	# use L-BFGS
	# use longer horizons to compute the loss (currently the horizon is 1 step)
		# maybe we just need to sum the gradients along the horizon
		# maybe we need to chain them together using th chain rule (not sure how stable this is)

include("../quasi_newton.jl")
function d2data(d)
	data = [d[1]; 0; 0; d[2]; 0; 0; 0]
	return data
end
∇d2data = ForwardDiff.jacobian(d -> d2data(d), zeros(2))

d0 = [0.40, 1.0]
lower = [0.0, 0.0]
upper = [0.80, 2.0]

function f0(d; rot=0)
	return clean_loss(:sphere, pairs0, d2data(d), n_sample=15, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))[1]
end

function fgH0(d; rot=0)
	f, g, H = clean_loss(:sphere, pairs0, d2data(d), n_sample=15, rot=rot, opts=SolverOptions(btol=3e-4, rtol=3e-4))
	return f, ∇d2data' * g, ∇d2data' * H * ∇d2data
end

dsol = quasi_newton_solve(f0, fgH0, d0, iter=200, gtol=1e-8, ftol=3e-5,
	H0=1e1*Matrix(Diagonal(ones(2))), lower=lower, upper=upper)

d0
f0(d0)
f0(dsol)
f0([0.1, 0.50])

using Plots; pyplot()
x=range(0,stop=0.80,length=10)
y=range(0,stop=0.75,length=10)
f(x,y) = log(10, clean_loss(:sphere, pairs0, d2data([x,y]),
	opts=SolverOptions(btol=3e-4, rtol=3e-4))[1])
plot(x,y,f,st=:surface,camera=(20,40))

# We can learn the coefficient of friction and the radius form 15*0.75 seconds of
# recording. We use the simulator to evaluate the loss and its gradients by differentiating
# through the simulator. With gradient information we can use L-BFGS

################################################################################
# Visualization
################################################################################


text:
using Base.Test
using RPOMDPModels, RPOMDPToolbox
using RobustValueIteration
const RPBVI = RobustValueIteration
TOL = 1e-6

@testset "Robust point-based value iteration" begin
    # value function max difference
    rp = BabyRPOMDP()
    alphavecs = [RPBVI.AlphaVec([2.0, 4.0], true),
                    RPBVI.AlphaVec([1.0, 5.0], :nothing),
                    RPBVI.AlphaVec([3.0, 3.0], :nothing),
                    RPBVI.AlphaVec([4.0, 2.0], :nothing),
                    RPBVI.AlphaVec([4.0, 3.0], :feed),
                    RPBVI.AlphaVec([5.0, 1.0], :nothing)]
    av2 = [RPBVI.AlphaVec([2.0, 4.0], true),
                    RPBVI.AlphaVec([1.1, 5.0], :nothing),
                    RPBVI.AlphaVec([3.0, 3.0], :nothing),
                    RPBVI.AlphaVec([4.0, 2.0], :nothing),
                    RPBVI.AlphaVec([4.0, 3.0], :feed),
                    RPBVI.AlphaVec([5.0, 0.9], :nothing)]
    @test RPBVI.diffvalue(alphavecs, av2, rp) ≈ 0.1 atol = TOL

    # minimum probability distribution (minutil)
    srand(429)
    bset = [[0, 1.0], [0.2, 0.8], [0.4, 0.6], [0.6,0.4], [0.8,0.2], [1.0,0.0]]
    alphas = [avec.alpha for avec in alphavecs]
    rip = BabyRIPOMDP()
    b = bset[2]
    a = :feed
    u, pstar = RPBVI.minutil(rip, b, a, alphas)
    @test sum(u) ≈ 4.8925 atol = TOL
    @test all(isapprox.(sum(pstar, (1,2)), 1.0, atol = TOL)) # p contains distributions

    # find αz that optimizes Eq. (5)
    nz = n_observations(rip)
    αz = Array{Array{Float64}}(nz)
    for zind = 1:nz
        αz[zind] = RPBVI.findαz(zind, u, b, pstar, alphas)
    end
    @test αz[1] == [1.0, 5.0]

    # find α*
    s = :hungry
    a = :nothing
    αstar = RPBVI.findαstar(rip, b, s, a, pstar, αz)
    @test αstar ≈ 3.40325 atol = TOL

    # robust point based dp backup
    srand(257349)
    αset = RPBVI.robustdpupdate(alphavecs, bset, rip)
    @test length(αset) == 6
    @test αset[1].alpha[1] ≈ 3.40325 atol = TOL

    # solver
    srand(429)
    bset = [[0, 1.0], [0.2, 0.8], [0.4, 0.6], [0.6,0.4], [0.8,0.2], [1.0,0.0]]
    p = BabyPOMDP()
    ip = BabyIPOMDP()
    rp = BabyRPOMDP()
    rip = BabyRIPOMDP()
    solver = RPBVISolver()
    solp = RPBVI.solve(solver, p)
    solip = RPBVI.solve(solver, ip)
    solrp = RPBVI.solve(solver, rp)
    solrip = RPBVI.solve(solver, rip)
    @test solp.action_map == [:nothing, :nothing, :feed, :feed, :feed, :feed]
    @test solip.action_map == [:feed, :feed, :feed, :feed, :feed, :nothing]
    @test solrp.action_map == [:nothing, :nothing, :feed, :feed, :feed, :feed]
    @test solrip.action_map == [:feed, :feed, :feed, :feed, :feed, :nothing]

    # beleif updates
    b = SparseCat([:hungry, :full], [0.4, 0.6])
    @test update(updater(solp),b,:nothing,:crying).b[1] ≈ 0.8363636363636364 atol = TOL
    @test update(updater(solrp),b,:nothing,:crying).b[1] ≈ 0.8497474114178433 atol = TOL

    # simulate
    sim = RolloutSimulator(max_steps = 1000)
    pbu = updater(solp)
    ipbu = updater(solip)
    rpbu = updater(solrp)
    ripbu = updater(solrip)
    @test simulate(sim, p, solp, pbu)[1] ≈ -20.29439154412 atol = 1
    @test simulate(sim, ip, solip, ipbu)[1] ≈ 10.0 atol = 1
    @test simulate(sim, rp, solrp, rpbu)[1] ≈ -48.60001793571622 atol = 1
    @test simulate(sim, rip, solrip, ripbu)[1] ≈ 4.908130152129587 atol = 1
end # testset

# Oracle values for ambiguity = 0 are from SARSOP
# Oracle values for ambiguity = 0.001, 0.1 from RobustInfoPOMDP/correctness_test.jl
@testset "SARSOP Comparison: RPOMDP" begin
    vtol = 0.1
    srand(8473272)
    p1 = TigerPOMDP(0.95)
    p1_001 = TigerRPOMDP(0.95, 0.001)
    p1_1 = TigerRPOMDP(0.95, 0.1)
    p2 = Baby2POMDP(-5.0, -10.0, 0.9)
    p2_001 = BabyRPOMDP(-5.0, -10.0, 0.9, 0.001)
    p2_1 = BabyRPOMDP(-5.0, -10.0, 0.9, 0.1)
    bs = [[b, 1-b] for b in 0.0:0.05:1.0]
    maxiter = 100
    solver = RPBVISolver(beliefpoints = bs, max_iterations = maxiter)
    pol1 = RPBVI.solve(solver, p1)
    pol1_001 = RPBVI.solve(solver, p1_001)
    pol1_1 = RPBVI.solve(solver, p1_1)
    pol2 = RPBVI.solve(solver, p2)
    pol2_001 = RPBVI.solve(solver, p2_001)
    pol2_1 = RPBVI.solve(solver, p2_1)
    val1 = value(pol1, [0.5, 0.5])
    val1_001 = value(pol1_001, [0.5, 0.5])
    val1_1 = value(pol1_1, [0.5, 0.5])
    val2 = value(pol2, [0.0, 1.0])
    val2_001 = value(pol2_001, [0.0, 1.0])
    val2_1 = value(pol2_1, [0.0, 1.0])
    @test val1 ≈ 19.25 atol = vtol
    @test val1_001 ≈ 18.78 atol = vtol
    @test val1_1 ≈ -15.55 atol = vtol
    @test val2 ≈ -16.30 atol = vtol
    @test val2_001 ≈ -18.49 atol = vtol
    @test val2_1 ≈ -40.86 atol = vtol
end

max_stars_repo_name: ajkeith/RobustValueIteration

max_issues_repo_name: ajkeith/RobustValueIteration

max_forks_repo_name: ajkeith/RobustValueIteration


text:
export EpsilonGreedyExplorer, GreedyExplorer

using Random
using Distributions: Categorical
using Flux

"""
    EpsilonGreedyExplorer{T}(;kwargs...)
    EpsilonGreedyExplorer(ϵ) -> EpsilonGreedyExplorer{:linear}(; ϵ_stable = ϵ)

> Epsilon-greedy strategy: The best lever is selected for a proportion `1 - epsilon` of the trials, and a lever is selected at random (with uniform probability) for a proportion epsilon . [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

Two kinds of epsilon-decreasing strategy are implmented here (`linear` and `exp`).

> Epsilon-decreasing strategy: Similar to the epsilon-greedy strategy, except that the value of epsilon decreases as the experiment progresses, resulting in highly explorative behaviour at the start and highly exploitative behaviour at the finish.  - [Multi-armed_bandit](https://en.wikipedia.org/wiki/Multi-armed_bandit)

# Keywords

- `T::Symbol`: defines how to calculate the epsilon in the warmup steps. Supported values are `linear` and `exp`.
- `step::Int = 1`: record the current step.
- `ϵ_init::Float64 = 1.0`: initial epsilon.
- `warmup_steps::Int=0`: the number of steps to use `ϵ_init`.
- `decay_steps::Int=0`: the number of steps for epsilon to decay from `ϵ_init` to `ϵ_stable`.
- `ϵ_stable::Float64`: the epsilon after `warmup_steps + decay_steps`.
- `is_break_tie=false`: randomly select an action of the same maximum values if set to `true`.
- `rng=Random.GLOBAL_RNG`: set the internal RNG.
- `is_training=true`, in training mode, `step` will not be updated. And the `ϵ` will be set to 0.

# Example

```julia
s = EpsilonGreedyExplorer{:linear}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="linear epsilon")
```
![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/linear_epsilon_greedy_selector.png)

```julia
s = EpsilonGreedyExplorer{:exp}(ϵ_init=0.9, ϵ_stable=0.1, warmup_steps=100, decay_steps=100)
plot([RL.get_ϵ(s, i) for i in 1:500], label="exp epsilon")
```
![](https://github.com/JuliaReinforcementLearning/ReinforcementLearning.jl/raw/master/docs/src/assets/img/exp_epsilon_greedy_selector.png)
"""
mutable struct EpsilonGreedyExplorer{Kind,IsBreakTie,R} <: AbstractExplorer
    ϵ_stable::Float64
    ϵ_init::Float64
    warmup_steps::Int
    decay_steps::Int
    step::Int
    rng::R
    is_training::Bool
end

function EpsilonGreedyExplorer(;
    ϵ_stable,
    kind = :linear,
    ϵ_init = 1.0,
    warmup_steps = 0,
    decay_steps = 0,
    step = 1,
    is_break_tie = false,
    is_training = true,
    rng = Random.GLOBAL_RNG,
)
    EpsilonGreedyExplorer{kind,is_break_tie,typeof(rng)}(
        ϵ_stable,
        ϵ_init,
        warmup_steps,
        decay_steps,
        step,
        rng,
        is_training,
    )
end

Flux.testmode!(p::EpsilonGreedyExplorer, mode = true) = p.is_training = !mode

EpsilonGreedyExplorer(ϵ; kwargs...) = EpsilonGreedyExplorer(; ϵ_stable = ϵ, kwargs...)

function get_ϵ(s::EpsilonGreedyExplorer{:linear}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    elseif step >= (s.warmup_steps + s.decay_steps)
        s.ϵ_stable
    else
        steps_left = s.warmup_steps + s.decay_steps - step
        s.ϵ_stable + steps_left / s.decay_steps * (s.ϵ_init - s.ϵ_stable)
    end
end

function get_ϵ(s::EpsilonGreedyExplorer{:exp}, step)
    if step <= s.warmup_steps
        s.ϵ_init
    else
        n = step - s.warmup_steps
        scale = s.ϵ_init - s.ϵ_stable
        s.ϵ_stable + scale * exp(-1.0 * n / s.decay_steps)
    end
end

get_ϵ(s::EpsilonGreedyExplorer) = s.is_training ? get_ϵ(s, s.step) : 0.0

"""
    (s::EpsilonGreedyExplorer)(values; step) where T

!!! note
    If multiple values with the same maximum value are found.
    Then a random one will be returned!

    `NaN` will be filtered unless all the values are `NaN`.
    In that case, a random one will be returned.
"""
function (s::EpsilonGreedyExplorer{<:Any,true})(values)
    ϵ = get_ϵ(s)
    s.is_training && (s.step += 1)
    rand(s.rng) >= ϵ ? rand(s.rng, find_all_max(values)[2]) : rand(s.rng, 1:length(values))
end

function (s::EpsilonGreedyExplorer{<:Any,false})(values)
    ϵ = get_ϵ(s)
    s.is_training && (s.step += 1)
    rand(s.rng) >= ϵ ? findmax(values)[2] : rand(s.rng, 1:length(values))
end

function (s::EpsilonGreedyExplorer{<:Any,true})(values, mask)
    ϵ = get_ϵ(s)
    s.is_training && (s.step += 1)
    rand(s.rng) >= ϵ ? rand(s.rng, find_all_max(values, mask)[2]) :
    rand(s.rng, findall(mask))
end

function (s::EpsilonGreedyExplorer{<:Any,false})(values, mask)
    ϵ = get_ϵ(s)
    s.is_training && (s.step += 1)
    rand(s.rng) >= ϵ ? findmax(values, mask)[2] : rand(s.rng, findall(mask))
end

Random.seed!(s::EpsilonGreedyExplorer, seed) = Random.seed!(s.rng, seed)

"""
    get_prob(s::EpsilonGreedyExplorer, values) ->Categorical
    get_prob(s::EpsilonGreedyExplorer, values, mask) ->Categorical

Return the probability of selecting each action given the estimated `values` of each action.
"""
function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,true}, values)
    ϵ, n = get_ϵ(s), length(values)
    probs = fill(ϵ / n, n)
    max_val_inds = find_all_max(values)[2]
    for ind in max_val_inds
        probs[ind] += (1 - ϵ) / length(max_val_inds)
    end
    Categorical(probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,true}, values, action::Integer)
    ϵ, n = get_ϵ(s), length(values)
    max_val_inds = find_all_max(values)[2]
    if action in max_val_inds
        ϵ / n + (1 - ϵ) / length(max_val_inds)
    else
        ϵ / n
    end
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,false}, values)
    ϵ, n = get_ϵ(s), length(values)
    probs = fill(ϵ / n, n)
    probs[findmax(values)[2]] += 1 - ϵ
    Categorical(probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,false}, values, action::Integer)
    ϵ, n = get_ϵ(s), length(values)
    if action == findmax(values)[2]
        ϵ / n + 1 - ϵ
    else
        ϵ / n
    end
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,true}, values, mask)
    ϵ, n = get_ϵ(s), length(values)
    probs = zeros(n)
    probs[mask] .= ϵ / sum(mask)
    max_val_inds = find_all_max(values, mask)[2]
    for ind in max_val_inds
        probs[ind] += (1 - ϵ) / length(max_val_inds)
    end
    Categorical(probs)
end

function RLBase.get_prob(s::EpsilonGreedyExplorer{<:Any,false}, values, mask)
    ϵ, n = get_ϵ(s), length(values)
    probs = zeros(n)
    probs[mask] .= ϵ / sum(mask)
    probs[findmax(values, mask)[2]] += 1 - ϵ
    Categorical(probs)
end

# Though we can achieve the same goal by setting the ϵ of [`EpsilonGreedyExplorer`](@ref) to 0,
# the GreedyExplorer is much faster.
struct GreedyExplorer <: AbstractExplorer end

(s::GreedyExplorer)(values) = findmax(values)[2]
(s::GreedyExplorer)(values, mask) = findmax(values, mask)[2]

function RLBase.get_prob(s::GreedyExplorer, values)
    prob = zeros(length(values))
    prob[findmax(values)[2]] = 1.0
    Categorical(prob)
end

RLBase.get_prob(s::GreedyExplorer, values, action::Integer) =
    findmax(values)[2] == action ? 1.0 : 0.0

function RLBase.get_prob(s::GreedyExplorer, values, mask)
    prob = zeros(length(values))
    prob[findmax(values, mask)[2]] = 1.0
    Categorical(prob)
end


text:
function testvertexcache(sd::DelaunayMeshes.SubDivision)
  for (vi, v) in enumerate(sd.vertexCache)
    @assert v >= 1
    @assert sd.edges[v].origin == vi
  end
end

function setuptwopointtesselation(v1::DelaunayMeshes.Point, v2::DelaunayMeshes.Point)
    local tess = DelaunayMeshes.DelaunayTesselation()

    # push two more points and connect them
    push!(tess.vertices, v1)
    push!(tess.vertices, v2)
    push!(tess.vertexCache, -1)
    push!(tess.vertexCache, -1)

    local ea = 1
    local eb = 5
    local ec = 9
    local ed = DelaunayMeshes.makeedge!(tess)
    DelaunayMeshes.endpoints!(tess, ed, 1, 4, 1, 1)
    DelaunayMeshes.splice!(tess, ed, ea)
    local ee = DelaunayMeshes.connect!(tess, ea, DelaunayMeshes.sym(ed), [1.0,1.0])
    local ef = DelaunayMeshes.connect!(tess, eb, DelaunayMeshes.sym(ee), [2.0,2.0])

    local eg = DelaunayMeshes.makeedge!(tess)
    DelaunayMeshes.endpoints!(tess, eg, 4, 5, 3, 3)
    DelaunayMeshes.splice!(tess, eg, DelaunayMeshes.sym(ee))
    local eh = DelaunayMeshes.connect!(tess, DelaunayMeshes.sym(ee), DelaunayMeshes.sym(eg), [3.0,3.0])
    local ei = DelaunayMeshes.connect!(tess, eb, DelaunayMeshes.sym(eh), [4.0,4.0])

    #=
    for (i,e) in enumerate(tess.edges)
      println("$i, $e")
    end
    =#
    return tess
end

function locateandassertintriangle(
  sd::DelaunayMeshes.DelaunayTesselation, x::DelaunayMeshes.Point, ei::Int)
  local e1 = DelaunayMeshes.locatevertex(sd, x, ei)
  local e2 = DelaunayMeshes.lnext(sd, e1)
  local e3 = DelaunayMeshes.lnext(sd, e2)
  #println("$e1")

  return  !DelaunayMeshes.rightof(sd, x, e1) && !DelaunayMeshes.rightof(sd, x, e2) && !DelaunayMeshes.rightof(sd, x, e3)
end

function testdelaunayness(sd::DelaunayMeshes.DelaunayTesselation)
  local triangles = DelaunayMeshes.getalltriangles(sd)
  local faultyTriangles = Vector{Tuple{Int, Int, Int, Int}}()

  for curTriangle in triangles
    local v1 = DelaunayMeshes.org(sd, curTriangle[1])
    local v2 = DelaunayMeshes.org(sd, curTriangle[2])
    local v3 = DelaunayMeshes.org(sd, curTriangle[3])

    if (v1 > 3 && v2 > 3 && v3 > 3)
      for (v, vv) in enumerate(sd.vertices)
        if (v != v1 && v != v2 && v != v3)
          if (DelaunayMeshes.incircle(sd, v1, v2, v3, v))
            push!(faultyTriangles, (v1, v2, v3, v))
          end
        end
      end
    end
  end

  if (length(faultyTriangles) > 0)
    for f in faultyTriangles
      println("$f: $(sd.vertices[f[1]]), $(sd.vertices[f[2]]), $(sd.vertices[f[3]]), $(sd.vertices[f[4]])")
    end
  end

  return faultyTriangles
end

function insertsingleedge!(sd::DelaunayMeshes.SubDivision, v0::Int, v1::Int, lfi::Int, rfi::Int)
  local newEdge = DelaunayMeshes.makeedge!(sd)
  DelaunayMeshes.endpoints!(sd, newEdge, v0, v1, lfi, rfi)
  return newEdge
end

function edgeconnectsvertices(sd::DelaunayMeshes.SubDivision, e0::Int, v0::Int, v1::Int)
  local vo = DelaunayMeshes.org(sd, e0)
  local vd = DelaunayMeshes.dest(sd, e0)
  return (vo == v0 && vd == v1) || (vo == v1 && vd == v0)
end

function findedgeconnectingvertices(sd::DelaunayMeshes.SubDivision, v0::Int, v1::Int)
  local estart = sd.vertexCache[v0]
  local ecur = DelaunayMeshes.onext(sd, estart)

  while (DelaunayMeshes.dest(sd, ecur) != v1 && ecur != estart)
    ecur = DelaunayMeshes.onext(sd, ecur)
  end

  @assert DelaunayMeshes.dest(sd, ecur) == v1
  return ecur
end

function validateconstraints(sd::DelaunayMeshes.SubDivision, constraints::Array{Int})

  # collect all ordered constraints from tesselation
  local foundConstraints = Array{Int, 2}()
  local baseEdges = Array{Int}()
  for i=1:4:length(sd.edges)
    if (DelaunayMeshes.isconstraint(sd, i))
      local baseInd = sd.edges[i].constraintType == DelaunayMeshes.OnRightSide ? i : i + 2
      local ov = DelaunayMeshes.org(sd, baseInd)
      local dv = DelaunayMeshes.dest(sd, baseInd)
      if (size(foundConstraints, 1) == 0)
        foundConstraints = [ov dv]
        baseEdges = [baseInd]
      else
        foundConstraints = [foundConstraints; ov dv]
        baseEdges = [baseEdges; baseInd]
      end
    end
  end

  local sortedFoundConstraints = sortrows(foundConstraints)
  @test sortedFoundConstraints == sortrows(constraints)

  # check that constraints are not intersected by any other edges
  for bi in baseEdges
    local bov = DelaunayMeshes.org(sd, bi)
    local bdv = DelaunayMeshes.dest(sd, bi)
    for i=1:4:length(sd.edges)
      local tov = DelaunayMeshes.org(sd, i)
      local tdv = DelaunayMeshes.dest(sd, i)
      # check for intersection if they do not have common vertices
      if (bov != tov && bov != tdv && bdv != tov && bdv != tdv)
         @test !(DelaunayMeshes.edgeintersects(sd, bi, i))
      end
    end
  end
end

function getrandomnumbersinsideboundingbox(bb::Vector{Float64}, n::Int)
  local sx = (bb[2] - bb[1])
  local sy = (bb[4] - bb[3])
  local pointsX = rand(Float64, n, 1) * sx + bb[1]
  local pointsY = rand(Float64, n, 1) * sy + bb[3]
  return [pointsX pointsY]
end

function pointsinconstrainedregion(
  mesh::DelaunayMeshes.Mesh, points::Array{Float64, 2})

  local sb = DelaunayMeshes.scalepoints(mesh,
    [
      -1.0 -1.0;
      -0.5 -1.0;
      -0.5 0.0;
      -1.0 0.0;
      -1.0 0.0;
       1.0 0.0;
       1.0 1.0;
      -1.0 1.0;
       0.5 -1.0;
       1.0 -1.0;
       1.0 1.0;
       0.5 1.0;
    ]
  )

  local npoints = size(points, 1)
  pointsInRegion = fill(false, npoints)

  for i in 1:npoints
     local cp = points[i,:]
     local inSquare1 =
     [
       cp[1] >= sb[1,1] cp[2] >= sb[1,2];
       cp[1] <= sb[2,1] cp[2] >= sb[2,2];
       cp[1] <= sb[3,1] cp[2] <= sb[3,2];
       cp[1] >= sb[4,1] cp[2] <= sb[4,2];
     ]
    if (inSquare1 == trues(4,2))
      pointsInRegion[i] = true
    end
    local inSquare2 =
     [
       cp[1] >= sb[5,1] cp[2] >= sb[5,2];
       cp[1] <= sb[6,1] cp[2] >= sb[6,2];
       cp[1] <= sb[7,1] cp[2] <= sb[7,2];
       cp[1] >= sb[8,1] cp[2] <= sb[8,2];
     ]
    if (inSquare2 == trues(4,2))
      pointsInRegion[i] = true
    end
    local inSquare3 =
     [
       cp[1] >= sb[9,1] cp[2] >= sb[9,2];
       cp[1] <= sb[10,1] cp[2] >= sb[10,2];
       cp[1] <= sb[11,1] cp[2] <= sb[11,2];
       cp[1] >= sb[12,1] cp[2] <= sb[12,2];
     ]
    if (inSquare3 == trues(4,2))
      pointsInRegion[i] = true
    end
  end

  return pointsInRegion
end

function getplotedges(sd::DelaunayMeshes.DelaunayTesselation)
  local x = Vector{Float64}()
  local y = Vector{Float64}()
  for i in 1:4:length(sd.edges)
    local vo = DelaunayMeshes.org(sd, i)
    local vd = DelaunayMeshes.dest(sd, i)
    if (vo > -1)
      push!(x, sd.vertices[vo][1])
      push!(x, sd.vertices[vd][1])
      push!(x, NaN)
      push!(y, sd.vertices[vo][2])
      push!(y, sd.vertices[vd][2])
      push!(y, NaN)
    end
  end

  return x, y
end

function create_mesh_with_triangle(cx::Float64, cy::Float64, r::Float64, t1::Float64, t2::Float64, t3::Float64)
    local v1 = [cx + r * cos(t1), cy + r * sin(t1)]
    local v2 = [cx + r * cos(t2), cy + r * sin(t2)]
    local v3 = [cx + r * cos(t3), cy + r * sin(t3)]
    local mesh = DelaunayMeshes.Mesh()
    mesh.tesselation.vertices = [v1, v2, v3]

    return mesh
end

function check_conformity(sd::DelaunayMeshes.DelaunayTesselation)
  for ei in 1:2:(size(sd.edges)[1])

    # check that edge is part of a triangle (not more than three vertices)
    local e1 = DelaunayMeshes.lnext(sd, ei)
    local e2 = DelaunayMeshes.lnext(sd, e1)
    local e3 = DelaunayMeshes.lnext(sd, e2)
    @test e3 == ei

    # check that any vertex that is on this segment is either org or dest
    local vo = DelaunayMeshes.org(sd, ei)
    local vd = DelaunayMeshes.dest(sd, ei)
    local x1 = sd.vertices[vo]
    local x2 = sd.vertices[vd]

    for xi in sd.vertices
      if !((xi ≈ x1) || (xi ≈ x2))
        @test !DelaunayMeshes.onlinesegment(xi, x1, x2)
      end
    end
  end
end


text:
function rvrp_data_unit_tests()

end


text:
using BrainWave
#using PlotlyJS

sf = 512
dur = 1.2
nSamp = round(Int, sf*dur)
tArr = collect(0:nSamp-1)./sf
dat = zeros(2, nSamp)
dat[1,:] = sin.(2*pi*75 .* tArr)
dat[2,:] = sin.(2*pi*120 .* tArr)

p1, f1 = getSpectrum(dat[1,:], sf, powerOfTwo=false)
p1b, f1b = getSpectrum(dat[1,:], sf, powerOfTwo=true)

p2, f2 = getSpectrum(dat[2,:], sf, powerOfTwo=false)
p2b, f2b = getSpectrum(dat[2,:], sf, powerOfTwo=true)

# s1 = scatter(;x=f1, y=10*log10.(p1))
# s1b = scatter(;x=f1b, y=10*log10.(p1b))

# s2 = scatter(;x=f2, y=10*log10.(p2))
# s2b = scatter(;x=f2b, y=10*log10.(p2b))


text:
struct bayes_network_model <: distribution_hyper_params
    κ::Float64
    m::AbstractArray{Float64}
    ν::Float64
    ψ::AbstractArray{Float64}
    count::Int64
    ms::Vector{AbstractArray{Float64}}
    λ::Float64
end

mutable struct bayes_network_sufficient_statistics <: sufficient_statistics
    N::Float64
    Ngroups::Vector{Float64}
    points_sum::AbstractArray{Float64,1}
    mu_vector::Vector{AbstractArray{Float64,1}}
    S::AbstractArray{Float64,2}
end

function calc_posterior(prior::bayes_network_model, suff_statistics::bayes_network_sufficient_statistics)
    if suff_statistics.N == 0
        return prior
    end
    κ = prior.κ + suff_statistics.N
    ν = prior.ν + suff_statistics.N
    m = (prior.m.*prior.κ + suff_statistics.points_sum) / κ
    ψ = (prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m'+ suff_statistics.S) / ν
    ψ = Matrix(Hermitian(ψ))
    #ψ = Hermitian((prior.ν * prior.ψ + prior.κ*prior.m*prior.m' -κ*m*m' + suff_statistics.S) ./ ν)
    if  isposdef(ψ) == false
        println(ψ)
        println(m)
        println(κ*m*m')
        println(suff_statistics)
    end
    return bayes_network_model(κ,m,ν,ψ,prior.count,suff_statistics.mu_vector,prior.λ)
end

function sample_distribution(hyperparams::bayes_network_model)
    Σ = rand(Distributions.InverseWishart(hyperparams.ν, hyperparams.ν* hyperparams.ψ))
    mu_vector = Vector{AbstractArray{Float64,1}}()
    μ = rand(Distributions.MvNormal(hyperparams.m, Σ/hyperparams.κ))
    for i=1:hyperparams.count
        push!(mu_vector, rand(Distributions.MvNormal(hyperparams.ms[i], Σ/hyperparams.κ)))
    end
    return mv_group_gaussian(mu_vector,Σ,inv(Σ),logdet(Σ))
end

function create_sufficient_statistics(hyper::bayes_network_model,
        posterior::bayes_network_model,
        points::AbstractArray{Float64,2},
        point_to_group = 0)
    if size(points,2) == 0
        return bayes_network_sufficient_statistics(0,
            [0 for i=1:hyper.count],
            zeros(length(hyper.m)),
            [hyper.m for i=1:hyper.count],
            zeros(length(hyper.m),length(hyper.m)))
    end
    pts_dict = Dict()
    for i=1:hyper.count
        pts_dict[i] = @view points[:,point_to_group .== i]
    end
    mu_vector = create_mu_vector(pts_dict, posterior.ψ,hyper.λ)
    dim = size(posterior.ψ,1)
    S = zeros(dim,dim)
    points_sum = zeros(dim)
    Ngroups = Vector{Float64}()
    N = 0
    for i=1:hyper.count
        pts = Array(pts_dict[i])
        if size(pts,2) == 0
            push!(Ngroups,0)
        else
            movedPts = pts .- mu_vector[i]
            S += movedPts * movedPts'
            points_sum += sum(movedPts, dims = 2)
            push!(Ngroups,size(pts,2))
            N += size(pts,2)
        end
    end
    return bayes_network_sufficient_statistics(N,Ngroups,points_sum[:],mu_vector,S)
end

function log_marginal_likelihood(hyper::bayes_network_model, posterior_hyper::bayes_network_model, suff_stats::bayes_network_sufficient_statistics)
    D = size(suff_stats.S,1)
    logpi = log(pi)
    return -suff_stats.N*D/2*logpi +
        log_multivariate_gamma(posterior_hyper.ν/2, D)-
        log_multivariate_gamma(hyper.ν/2, D) +
         (hyper.ν/2)*logdet(hyper.ψ*hyper.ν)-
         (posterior_hyper.ν/2)*logdet(posterior_hyper.ψ*posterior_hyper.ν) +
         (D/2)*(log(hyper.κ)-(D/2)*log(posterior_hyper.κ))
end

function aggregate_suff_stats(suff_l::bayes_network_sufficient_statistics, suff_r::bayes_network_sufficient_statistics)
    new_suff = deepcopy(suff_l)
    N = 0
    for i=1:length(new_suff.Ngroups)
        new_suff.mu_vector[i] = (new_suff.mu_vector[i] .* (new_suff.Ngroups[i] / (new_suff.Ngroups[i] + suff_r.Ngroups[i]))) +(suff_r.mu_vector[i] .* (suff_r.Ngroups[i] / (new_suff.Ngroups[i] + suff_r.Ngroups[i])))
        new_suff.Ngroups[i] += suff_r.Ngroups[i]
        new_suff.N += suff_r.Ngroups[i]
    end
    new_suff.S += suff_r.S
    new_suff.points_sum += suff_r.points_sum
    return new_suff
end

function create_mu_vector(points::Dict, Σ::AbstractArray{Float64}, λ::Float64)
    group_count = length(keys(points))
    A, b = create_matrix_for_least_squares(points, Σ, λ, group_count)
    xhat = inv(A'*A)*(A'*b)
    dim = size(Σ,1)
    mu_vector = reshape(xhat, dim, group_count)
    return [mu_vector[:,i] for i=1:group_count]
end

function create_matrix_for_least_squares(points::Dict, Σ::AbstractArray{Float64}, λ::Float64, group_count::Int64)
    points_count = zeros(group_count)
    points_sum::Int64 = 0
    dim = size(Σ,1)
    for i=1:length(points)
        points_count[i] = (points[i] == -1 ? 0 : size(points[i],2))
        points_sum += points_count[i]
    end
    points = [x for x in values(points) if x != -1]
    points_arr = vcat(points, [zeros(dim) for x in 1:group_count])

    points_arr = reduce(hcat, points_arr)

    points_arr = points_arr[:]

    A = zeros(length(points_arr),groups_count * dim)
    L = cholesky(Σ).U
    Linv = Matrix(inv(L))

    for i=1:points_sum
        b =  points_arr[(i-1)*dim+1:i*dim]
        tmp = Linv * b

        points_arr[(i-1)*dim+1:i*dim] = Linv * points_arr[(i-1)*dim+1:i*dim]
    end
    count = 1

    for i=1:group_count
        for j=1:points_count[i]
            A[count:count+dim-1,dim*(i-1)+1:dim*i] = Linv
            count+= dim
        end
    end

    λsqrt = λ^2
    count += dim #Skip the first μ

    for i=1:(size(A,2) - dim)
        A[count,i] = -λsqrt
        A[count,i+dim] = λsqrt
        count+=1
    end
    return A, points_arr
end

function create_bayes_model_hyper_from_niw(niw::niw_hyperparams, λ, count)
    ms = [niw.m for i=1:count]
    return bayes_network_model(niw.κ,niw.m,niw.ν,niw.ψ,count,ms,λ)
end


text:
#= tilde, hat, etc, definitions

We use these macros as shorthand notation. For instance

  @tilde_outer("B1")
  @hat_outer("B2")

should expand to

  B1t = B1_x - (Fx + xi_x) * B1p
  B2h = B2_y - (Fy + xi_y) * B2p

etc.

=#
macro tilde_outer(fname::String)
    ft    = Symbol(fname, "t")
    f_x   = Symbol(fname, "_x")
    fp    = Symbol(fname, "p")
    return esc( :($ft = $f_x - (Fx + xi_x) * $fp) )
end
macro hat_outer(fname::String)
    fh    = Symbol(fname, "h")
    f_y   = Symbol(fname, "_y")
    fp    = Symbol(fname, "p")
    return esc( :($fh = $f_y - (Fy + xi_y) * $fp) )
end
macro bar_outer(fname::String)
    fb    = Symbol(fname, "b")
    f_xx  = Symbol(fname, "_xx")
    fpp   = Symbol(fname, "pp")
    fp_x  = Symbol(fname, "p_x")
    return esc( :($fb = $f_xx + (Fx + xi_x) * ( -2*($fp_x) + (Fx + xi_x) * ($fpp) )) )
end
macro star_outer(fname::String)
    fs    = Symbol(fname, "s")
    f_yy  = Symbol(fname, "_yy")
    fpp   = Symbol(fname, "pp")
    fp_y  = Symbol(fname, "p_y")
    return esc( :($fs = $f_yy + (Fy + xi_y) * ( -2*($fp_y) + (Fy + xi_y) * ($fpp) )) )
end
macro cross_outer(fname::String)
    fc    = Symbol(fname, "c")
    f_xy  = Symbol(fname, "_xy")
    fpp   = Symbol(fname, "pp")
    fp_x  = Symbol(fname, "p_x")
    fp_y  = Symbol(fname, "p_y")
    return esc( :($fc = $f_xy  - (Fx + xi_x) * ($fp_y) -
                  (Fy + xi_y) * ( $fp_x - (Fx + xi_x) * ($fpp) ) ) )
end

# assuming
# (A d_uu + B d_u + C Id) f = -S

function S_eq_coeff!(ABCS::Vector, vars::Tuple, ::Outer)
    (phi0, u, xi, B1, B1p, B2, B2p, G, Gp, phi, phip) = vars

    ABCS[1] = *(6, u ^ 4)
    ABCS[2] = *(12, u ^ 3)
    ABCS[3] = Gp ^ 2 + *(3, B2p ^ 2) + *(4, phip ^ 2) + *(B1p ^ 2, cosh(G) ^ 2)
    ABCS[4] = 0

    nothing
end

# this is a coupled equation for Fx and Fy. the notation used is
#
# ( A11 d_uu Fx + A12 d_uu Fy + B11 d_u Fx + B12 d_u Fy + C11 Fx + C12 Fy ) = -S1
# ( A21 d_uu Fx + A22 d_uu Fy + B21 d_u Fx + B22 d_u Fy + C21 Fx + C22 Fy ) = -S2

function Fxy_eq_coeff!(AA::Matrix, BB::Matrix, CC::Matrix, SS::Vector, vars::Tuple, ::Outer)
    (
        phi0, u, xi, xi_x, xi_y,
        B1     ,    B2     ,    G      ,    phi    ,    S      ,
        B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,
        B1pp   ,    B2pp   ,    Gpp    ,                Spp    ,
        B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,
        B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,
        B1p_x  ,    B2p_x  ,    Gp_x   ,                Sp_x   ,
        B1p_y  ,    B2p_y  ,    Gp_y   ,                Sp_y
    ) = vars

    expB1   = exp(B1)
    sinhG   = sinh(G)
    sinh2G  = sinh(*(2, G))
    coshG   = cosh(G)
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2

    AA[1,1] = *(2, S ^ 2, u ^ 4, expB1)
    AA[1,2] = 0
    AA[2,1] = 0
    AA[2,2] = *(2, S ^ 2, u ^ 4)

    BB[1,1] = *(2, S, u ^ 2, *(-1, Sp) + *(-1, B2p, S) + *(2, S, u) + *(-1, B1p, S, coshGsq), expB1)
    BB[1,2] = *(2, S ^ 2, u ^ 2, Gp + *(B1p, coshG, sinhG))
    BB[2,1] = *(2, S ^ 2, u ^ 2, Gp + *(-1, B1p, coshG, sinhG), expB1)
    BB[2,2] = *(S, u ^ 2, *(-2, Sp) + *(S, B1p + *(-2, B2p) + *(4, u)) + *(B1p, S, cosh2G))

    CC[1,1] = *(-2, *(4, Sp ^ 2) + *(-1, B2pp, S ^ 2) + *(-1, Gp ^ 2, S ^ 2) + *(-4, S, Spp) + *(-4, phip ^ 2, S ^ 2) + *(-3, B2p ^ 2, S ^ 2) + *(-1, S, coshGsq, *(B1pp, S) + *(S, B1p ^ 2) + *(3, B1p, Sp)) + *(-3, B2p, S, Sp) + *(-1, B1p, Gp, S ^ 2, sinh2G), expB1)

    CC[1,2] = *(-1, S, *(*(B1pp, S) + *(-1, S, B1p ^ 2) + *(3, B1p, Sp), sinh2G) + *(2, Gpp, S) + *(6, Gp, Sp) + *(-2, B1p, Gp, S) + *(2, B1p, Gp, S, cosh2G))

    CC[2,1] = *(S, *(*(B1pp, S) + *(S, B1p ^ 2) + *(3, B1p, Sp), sinh2G) + *(-6, Gp, Sp) + *(-2, Gpp, S) + *(-2, B1p, Gp, S) + *(2, B1p, Gp, S, cosh2G), expB1)

    CC[2,2] = *(-8, Sp ^ 2) + *(2, B2pp, S ^ 2) + *(2, Gp ^ 2, S ^ 2) + *(6, B2p ^ 2, S ^ 2) + *(8, S, Spp) + *(8, phip ^ 2, S ^ 2) + *(2, S, coshGsq, *(S, B1p ^ 2) + *(-1, B1pp, S) + *(-3, B1p, Sp)) + *(6, B2p, S, Sp) + *(-2, B1p, Gp, S ^ 2, sinh2G)

    SS[1] = *(S ^ 2, *(2, Gp_y) + *(B1p_y + *(xi_y, B1p ^ 2) + *(-1, B1_y, B1p) + *(-1, B1pp, xi_y), sinh2G) + *(-2, B1_y, Gp) + *(-2, Gpp, xi_y) + *(2, B1p, Gp, xi_y) + *(2, B1p, G_y + *(-1, Gp, xi_y), cosh2G)) + *(*(S, *(-8, Sp_x) + *(-1, B2p + *(B1p, coshGsq), *(6, S_x) + *(-6, Sp, xi_x)) + *(8, Spp, xi_x)) + *(S ^ 2, *(-2, B2p_x) + *(-8, phip, phi_x + *(-1, phip, xi_x)) + *(-6, B2p, B2_x + *(-1, B2p, xi_x)) + *(-2, coshGsq, B1p_x + *(-1, B1pp, xi_x)) + *(-2, G_x + *(-1, Gp, xi_x), Gp + *(B1p, sinh2G)) + *(2, B2pp, xi_x) + *(2, B1p, coshGsq, *(-1, B1_x) + *(B1p, xi_x))) + *(8, Sp, S_x + *(-1, Sp, xi_x)), expB1) + *(S, Gp + *(B1p, coshG, sinhG), *(6, S_y) + *(-6, Sp, xi_y))

    SS[2] = *(S, *(-8, Sp_y) + *(-1, B2p + *(-1, B1p, coshGsq), *(6, S_y) + *(-6, Sp, xi_y)) + *(8, Spp, xi_y)) + *(2, S ^ 2, *(-1, B2p_y) + *(B2pp, xi_y) + *(xi_y, Gp ^ 2) + *(coshGsq, B1p_y + *(xi_y, B1p ^ 2) + *(-1, B1_y, B1p) + *(-1, B1pp, xi_y)) + *(-1, G_y, Gp) + *(-4, phi_y, phip) + *(-3, B2_y, B2p) + *(3, xi_y, B2p ^ 2) + *(4, xi_y, phip ^ 2) + *(B1p, G_y + *(-1, Gp, xi_y), sinh2G)) + *(8, Sp, S_y + *(-1, Sp, xi_y)) + *(S, *(S, *(2, Gp_x) + *(-1, B1p_x + *(B1_x, B1p) + *(-1, xi_x, B1pp + B1p ^ 2), sinh2G) + *(-2, xi_x, Gpp + *(B1p, Gp)) + *(2, B1_x, Gp) + *(-2, B1p, G_x + *(-1, Gp, xi_x), cosh2G)) + *(Gp + *(-1, B1p, coshG, sinhG), *(6, S_x) + *(-6, Sp, xi_x)), expB1)

    nothing
end

function Sd_eq_coeff!(ABCS::Vector, vars::Tuple, ::Outer)
    (
        potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,
        B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
        B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
        B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
        B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
        B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
        B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
        B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
        B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                    B2_xy  ,    G_xy   ,                S_xy
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")

    @bar_outer("B1")
    @bar_outer("B2")
    @bar_outer("G")
    @bar_outer("phi")
    @bar_outer("S")

    @star_outer("B1")
    @star_outer("B2")
    @star_outer("G")
    @star_outer("phi")
    @star_outer("S")

    @tilde_outer("Fxp")
    @tilde_outer("Fyp")

    @hat_outer("Fxp")
    @hat_outer("Fyp")

    @cross_outer("B2")
    @cross_outer("G")
    @cross_outer("S")

    VV = -3 -3/2 * phi^2 + phi^4*UU(phi,potential)

   # VVp = -3 * phi + phi^4 * UUp(phi,potential) + 4*phi^3 * UU(phi,potential)

    expB1   = exp(B1)
    expB2   = exp(B2)
    sinh2G  = sinh(*(2, G))
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2
    coshG   = cosh(G)
    sinhG   = sinh(G)

    ABCS[1] = 0

    ABCS[2] = *(-12, S ^ 3, u ^ 2, expB1)

    ABCS[3] = *(24, Sp, S ^ 2, expB1)

    ABCS[4] = *(*(S, *(*(-8, Gh, St) + *(-8, Gt, Sh), coshG) + *(*(-16, Sc) + *(2, Sh, Fxp + *(-4, B2t)) + *(2, Sp, *(4, Fxh) + *(4, Fyt) + *(8, xi_xy)) + *(2, St, Fyp + *(-4, B2h)), sinhG)) + *(S ^ 2, *(*(-4, Gc) + *(-2, Gh, B1t + B2t + *(-1, Fxp)) + *(2, Gp, Fxh + Fyt + *(2, xi_xy)) + *(2, Gt, B1h + Fyp + *(-1, B2h)), coshG) + *(*(-4, B2c) + *(2, Fxph) + *(2, Fypt) + *(-8, phih, phit) + *(-4, Gh, Gt) + *(2, B2h, Fxp + *(-4, B2t)) + *(2, B2p, Fxh + Fyt + *(2, xi_xy)) + *(2, Fyp, B2t + *(-1, Fxp)), sinhG)) + *(8, Sh, St, sinhG), exp(B1 + B2)) + *(*(S, *(*(8, Ss) + *(-2, Sh, Fyp + *(-4, B2h) + *(4, B1h)) + *(-2, Sp, *(4, Fyh) + *(4, xi_yy)), coshG) + *(8, Gh, Sh, sinhG)) + *(S ^ 2, *(*(2, Gs) + *(-2, Gp, Fyh + xi_yy) + *(2, Gh, B2h + *(-1, Fyp) + *(-2, B1h)), sinhG) + *(Fyp ^ 2 + *(-2, B1s) + *(-2, Fyph) + *(2, B2s) + *(2, B1h ^ 2) + *(2, Gh ^ 2) + *(4, B2h ^ 2) + *(4, phih ^ 2) + *(Fyp, *(-2, B2h) + *(2, B1h)) + *(-2, B1h, B2h) + *(2, B1p + *(-1, B2p), Fyh + xi_yy), coshG)) + *(-4, Sh ^ 2, coshG), expB2) + *(*(S, *(*(8, Sb) + *(-8, Fxt, Sp) + *(-8, Sp, xi_xx) + *(-2, Fxp, St) + *(8, B1t, St) + *(8, B2t, St), coshG) + *(8, Gt, St, sinhG)) + *(S ^ 2, *(*(2, Gb) + *(-2, Gp, Fxt + xi_xx) + *(2, Gt, B2t + *(-1, Fxp) + *(2, B1t)), sinhG) + *(Fxp ^ 2 + *(-2, Fxpt) + *(2, B1b) + *(2, B2b) + *(2, B1t ^ 2) + *(2, Gt ^ 2) + *(4, B2t ^ 2) + *(4, phit ^ 2) + *(-1, Fxp, *(2, B1t) + *(2, B2t)) + *(-2, B1p + B2p, Fxt + xi_xx) + *(2, B1t, B2t), coshG)) + *(-4, St ^ 2, coshG), exp(B2 + *(2, B1))) + *(8, VV, S ^ 4, expB1)

    nothing
end

function B2d_eq_coeff!(ABCS::Vector, vars::Tuple, ::Outer)
    (
        phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
        B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
        B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
        B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
        B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
        B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
        B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
        B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
        B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                    B2_xy  ,    G_xy   ,                S_xy
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")

    @bar_outer("B1")
    @bar_outer("B2")
    @bar_outer("G")
    @bar_outer("phi")
    @bar_outer("S")

    @star_outer("B1")
    @star_outer("B2")
    @star_outer("G")
    @star_outer("phi")
    @star_outer("S")

    @tilde_outer("Fxp")
    @tilde_outer("Fyp")

    @hat_outer("Fxp")
    @hat_outer("Fyp")

    @cross_outer("B2")
    @cross_outer("G")
    @cross_outer("S")

    expB1   = exp(B1)
    expB2   = exp(B2)
    sinh2G  = sinh(*(2, G))
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2
    coshG   = cosh(G)
    sinhG   = sinh(G)

    ABCS[1] = 0

    ABCS[2] = *(-12, S ^ 4, u ^ 2, expB1)

    ABCS[3] = *(18, Sp, S ^ 3, expB1)

    ABCS[4] = *(*(S, *(*(2, Gh, St) + *(2, Gt, Sh), coshG) + *(-1, *(-4, Sc) + *(-4, Fxp, Sh) + *(-4, Fyp, St) + *(2, Fxh, Sp) + *(2, Fyt, Sp) + *(4, B2h, St) + *(4, B2t, Sh) + *(4, Sp, xi_xy), sinhG)) + *(S ^ 2, *(*(4, Gc) + *(-2, Gp, Fxh + Fyt + *(2, xi_xy)) + *(-2, Gt, B1h + Fyp + *(2, B2h)) + *(2, Gh, B1t + *(-1, Fxp) + *(-2, B2t)), coshG) + *(*(-8, B2c) + *(-2, Fxph) + *(-2, Fypt) + *(2, Fyp, Fxp + *(2, B2t)) + *(4, B2h, Fxp + *(-1, B2t)) + *(4, B2p, Fxh + Fyt + *(2, xi_xy)) + *(4, Gh, Gt) + *(8, phih, phit), sinhG)) + *(-8, Sh, St, sinhG), exp(B1 + B2)) + *(*(S, *(*(-2, Ss) + *(2, Sh, B1h + *(-2, Fyp) + *(2, B2h)) + *(2, Sp, Fyh + xi_yy), coshG) + *(-2, Gh, Sh, sinhG)) + *(S ^ 2, *(*(-2, Gs) + *(2, Gh, Fyp + *(2, B1h) + *(2, B2h)) + *(2, Gp, Fyh + xi_yy), sinhG) + *(*(-1, Fyp ^ 2) + *(-4, phih ^ 2) + *(-2, B1h ^ 2) + *(-2, Gh ^ 2) + *(2, B1s) + *(2, Fyph) + *(2, B2h ^ 2) + *(4, B2s) + *(-1, Fyp, *(2, B1h) + *(4, B2h)) + *(-1, Fyh + xi_yy, *(2, B1p) + *(4, B2p)) + *(-4, B1h, B2h), coshG)) + *(4, Sh ^ 2, coshG), expB2) + *(*(S, *(*(-2, Sb) + *(-4, Fxp, St) + *(-2, B1t, St) + *(2, Fxt, Sp) + *(2, Sp, xi_xx) + *(4, B2t, St), coshG) + *(-2, Gt, St, sinhG)) + *(S ^ 2, *(*(-2, Gb) + *(2, Gp, Fxt + xi_xx) + *(2, Gt, Fxp + *(-2, B1t) + *(2, B2t)), sinhG) + *(*(-1, Fxp ^ 2) + *(-4, phit ^ 2) + *(-2, B1b) + *(-2, B1t ^ 2) + *(-2, Gt ^ 2) + *(2, Fxpt) + *(2, B2t ^ 2) + *(4, B2b) + *(Fxp, *(-4, B2t) + *(2, B1t)) + *(2, B1p + *(-2, B2p), Fxt + xi_xx) + *(4, B1t, B2t), coshG)) + *(4, St ^ 2, coshG), exp(B2 + *(2, B1))) + *(18, B2p, Sd, S ^ 3, expB1)

    nothing
end

# this is another coupled equation, for B1d and Gd. the notation used is
#
# ( A11 d_uu B1d + A12 d_uu Gd + B11 d_u B1d + B12 d_u Gd + C11 B1d + C12 Gd ) = -S1
# ( A21 d_uu B1d + A22 d_uu Gd + B21 d_u B1d + B22 d_u Gd + C21 B1d + C22 Gd ) = -S2

function B1dGd_eq_coeff!(AA::Matrix, BB::Matrix, CC::Matrix, SS::Vector, vars::Tuple, ::Outer)
    (
        phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
        B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
        B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
        B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
        B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
        B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
        B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
        B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
        B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                    B2_xy  ,    G_xy   ,                S_xy
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")

    @bar_outer("B1")
    @bar_outer("B2")
    @bar_outer("G")
    @bar_outer("phi")
    @bar_outer("S")

    @star_outer("B1")
    @star_outer("B2")
    @star_outer("G")
    @star_outer("phi")
    @star_outer("S")

    @tilde_outer("Fxp")
    @tilde_outer("Fyp")

    @hat_outer("Fxp")
    @hat_outer("Fyp")

    @cross_outer("B2")
    @cross_outer("G")
    @cross_outer("S")

    expB1   = exp(B1)
    expB2   = exp(B2)
    sinh2G  = sinh(*(2, G))
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2
    coshG   = cosh(G)
    tanhG   = tanh(G)
    sinhG   = sinh(G)
    sechG   = sech(G)

    AA[1,1] = 0
    AA[1,2] = 0
    AA[2,1] = 0
    AA[2,2] = 0

    BB[1,1] = *(-12, S ^ 4, u ^ 2, expB1)

    BB[1,2] = 0

    BB[2,1] = 0

    BB[2,2] = *(-12, S ^ 4, u ^ 2, expB1)

    CC[1,1] = *(6, S ^ 3, *(3, Sp) + *(2, Gp, S, tanhG), expB1)

    CC[1,2] = *(12, B1p, S ^ 4, expB1, tanhG)

    CC[2,1] = *(-6, B1p, S ^ 4, expB1, sinh2G)

    CC[2,2] = *(18, Sp, S ^ 3, expB1)

    SS[1] = *(3, *(-4, Sh ^ 2) + *(S, *(2, Ss) + *(-2, Sp, Fyh + xi_yy) + *(2, Sh, B2h + *(2, Fyp))) + *(S ^ 2, Fyp ^ 2 + *(-2, Fyph) + *(2, B2s) + *(4, B2h ^ 2) + *(4, phih ^ 2) + *(-2, B2h, Fyp) + *(-2, B2p, Fyh + xi_yy)), expB2, sechG) + *(3, *(4, St ^ 2) + *(-1, S ^ 2, Fxp ^ 2 + *(-2, Fxpt) + *(2, B2b) + *(4, B2t ^ 2) + *(4, phit ^ 2) + *(-2, B2p, Fxt + xi_xx) + *(-2, B2t, Fxp)) + *(2, S, *(-1, Sb) + *(Fxt, Sp) + *(Sp, xi_xx) + *(-1, B2t, St) + *(-2, Fxp, St)), exp(B2 + *(2, B1)), sechG) + *(6, S, *(Gt, Sh + *(S, B2h + *(-1, Fyp))) + *(-1, Gh, St) + *(Fxp, Gh, S) + *(Fyt, Gp, S) + *(-1, B2t, Gh, S) + *(-1, Fxh, Gp, S), exp(B1 + B2), sechG) + *(18, B1p, Sd, S ^ 3, expB1)

    SS[2] = *(3, *(4, Sh ^ 2) + *(-1, S, *(2, Ss) + *(-2, Sp, Fyh + xi_yy) + *(2, Sh, B2h + *(2, Fyp))) + *(-1, S ^ 2, Fyp ^ 2 + *(-2, Fyph) + *(2, B2s) + *(4, B2h ^ 2) + *(4, phih ^ 2) + *(-2, B2h, Fyp) + *(-2, B2p, Fyh + xi_yy)), expB2, sinhG) + *(3, *(4, St ^ 2) + *(-1, S ^ 2, Fxp ^ 2 + *(-2, Fxpt) + *(2, B2b) + *(4, B2t ^ 2) + *(4, phit ^ 2) + *(-2, B2p, Fxt + xi_xx) + *(-2, B2t, Fxp)) + *(2, S, *(-1, Sb) + *(Fxt, Sp) + *(Sp, xi_xx) + *(-1, B2t, St) + *(-2, Fxp, St)), exp(B2 + *(2, B1)), sinhG) + *(6, *(S, *(2, Sc) + *(Sh, B2t + *(-1, B1t) + *(2, Fxp)) + *(St, B1h + B2h + *(2, Fyp)) + *(-1, Sp, Fxh + Fyt + *(2, xi_xy))) + *(S ^ 2, *(-1, Fxph) + *(-1, Fypt) + *(2, B2c) + *(B1p, Fxh + *(-1, Fyt)) + *(B1t, Fyp) + *(B1h + *(-1, Fyp), B2t + *(-1, Fxp)) + *(-1, B2h, B1t + Fxp + *(-4, B2t)) + *(-1, B2p, Fxh + Fyt + *(2, xi_xy)) + *(4, phih, phit)) + *(-4, Sh, St), coshG, exp(B1 + B2)) + *(18, Gp, Sd, S ^ 3, expB1)

    nothing
end

function phid_eq_coeff!(ABCS::Vector, vars::Tuple, ::Outer)
    (
        potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1     ,    B2     ,    G      ,    phi    ,    S      ,    Fx     ,    Fy     ,  Sd,
        B1p    ,    B2p    ,    Gp     ,    phip   ,    Sp     ,    Fxp    ,    Fyp    ,
        B1pp   ,    B2pp   ,    Gpp    ,    phipp  ,    Spp    ,    Fxpp   ,    Fypp   ,
        B1_x   ,    B2_x   ,    G_x    ,    phi_x  ,    S_x    ,    Fx_x   ,    Fy_x   ,
        B1_y   ,    B2_y   ,    G_y    ,    phi_y  ,    S_y    ,    Fx_y   ,    Fy_y   ,
        B1p_x  ,    B2p_x  ,    Gp_x   ,    phip_x ,    Sp_x   ,    Fxp_x  ,    Fyp_x  ,
        B1p_y  ,    B2p_y  ,    Gp_y   ,    phip_y ,    Sp_y   ,    Fxp_y  ,    Fyp_y  ,
        B1_xx  ,    B2_xx  ,    G_xx   ,    phi_xx ,    S_xx   ,
        B1_yy  ,    B2_yy  ,    G_yy   ,    phi_yy ,    S_yy   ,
                    B2_xy  ,    G_xy   ,    phi_xy,     S_xy
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")

    # @bar_outer("B1")
    # @bar_outer("B2")
    # @bar_outer("G")
    @bar_outer("phi")
    # @bar_outer("S")

    # @star_outer("B1")
    # @star_outer("B2")
    # @star_outer("G")
    @star_outer("phi")
    # @star_outer("S")

    # @tilde_outer("Fxp")
    # @tilde_outer("Fyp")

    # @hat_outer("Fxp")
    # @hat_outer("Fyp")

    # @cross_outer("B2")
    # @cross_outer("G")
    # @cross_outer("S")
    @cross_outer("phi")

   # VV = -3 -3/2 * phi^2 + phi^4*UU(phi,potential)

    VVp = -3 * phi + phi^4 * UUp(phi,potential) + 4*phi^3 * UU(phi,potential)

    expB1   = exp(B1)
    expB2   = exp(B2)
    sinh2G  = sinh(*(2, G))
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2
    coshG   = cosh(G)
    sinhG   = sinh(G)

    ABCS[1] = 0

    ABCS[2] = *(-8, S ^ 3, u ^ 2, expB1)

    ABCS[3] = *(12, Sp, S ^ 2, expB1)

    ABCS[4] = *(*(*(4, phih, Sh + *(-1, S, B1h + Fyp + *(-1, B2h))) + *(4, S, phis + *(-1, phip, Fyh + xi_yy)), coshG) + *(4, Gh, phih, S, sinhG), expB2) + *(*(4, S, *(phib + *(phit, B1t + B2t + *(-1, Fxp)) + *(-1, phip, Fxt + xi_xx), coshG) + *(Gt, phit, sinhG)) + *(4, phit, St, coshG), exp(B2 + *(2, B1))) + *(*(-1, *(4, phih, St) + *(4, phit, Sh), sinhG) + *(S, *(-4, Gh, phit) + *(-4, Gt, phih), coshG) + *(S, *(-8, phic) + *(4, phih, Fxp + *(-1, B2t)) + *(4, phip, Fxh + Fyt + *(2, xi_xy)) + *(4, phit, Fyp + *(-1, B2h)), sinhG), exp(B1 + B2)) + *(-4, S ^ 2, *(S, VVp) + *(-3, phip, Sd), expB1)

    nothing
end

function A_eq_coeff!(ABCS::Vector, vars::Tuple, ::Outer)
    (
        potential, phi0, u, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1   , B2   , G   , phi   , S    , Fx    , Fy    , Sd, B1d, B2d, Gd, phid,
        B1p  , B2p  , Gp  , phip  , Sp   , Fxp   , Fyp   ,
        B1pp , B2pp , Gpp , phipp , Spp  , Fxpp  , Fypp  ,
        B1_x , B2_x , G_x , phi_x , S_x  , Fx_x  , Fy_x  ,
        B1_y , B2_y , G_y , phi_y , S_y  , Fx_y  , Fy_y  ,
        B1p_x, B2p_x, Gp_x, phip_x, Sp_x , Fxp_x , Fyp_x ,
        B1p_y, B2p_y, Gp_y, phip_y, Sp_y , Fxp_y , Fyp_y ,
        B1_xx, B2_xx, G_xx, phi_xx, S_xx ,
        B1_yy, B2_yy, G_yy, phi_yy, S_yy ,
               B2_xy, G_xy, phi_xy, S_xy
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")

    @bar_outer("B1")
    @bar_outer("B2")
    @bar_outer("G")
    @bar_outer("phi")
    @bar_outer("S")

    @star_outer("B1")
    @star_outer("B2")
    @star_outer("G")
    @star_outer("phi")
    @star_outer("S")

    @tilde_outer("Fxp")
    @tilde_outer("Fyp")

    @hat_outer("Fxp")
    @hat_outer("Fyp")

    @cross_outer("B2")
    @cross_outer("G")
    @cross_outer("S")
    @cross_outer("phi")

    VV = -3 -3/2 * phi^2 + phi^4*UU(phi,potential)

   # VVp = -3 * phi + phi^4 * UUp(phi,potential) + 4*phi^3 * UU(phi,potential)

    expB1   = exp(B1)
    expB2   = exp(B2)
    sinh2G  = sinh(*(2, G))
    cosh2G  = cosh(*(2, G))
    coshGsq = cosh(G)^2
    coshG   = cosh(G)
    sinhG   = sinh(G)

    ABCS[1] = *(6, S ^ 4, u ^ 4, expB1)

    ABCS[2] = *(12, S ^ 4, u ^ 3, expB1)

    ABCS[3] = 0

    ABCS[4] = *(*(2, S ^ 4, *(-4, VV) + *(3, Gd, Gp) + *(9, B2d, B2p) + *(12, phid, phip) + *(3, B1d, B1p, coshGsq)) + *(-72, Sd, Sp, S ^ 2), expB1) + *(*(S, *(*(-24, Ss) + *(24, Sh, B1h + *(-1, B2h)) + *(24, Sp, Fyh + xi_yy), coshG) + *(-24, Gh, Sh, sinhG)) + *(S ^ 2, *(*(-6, Gs) + *(-6, B2h, Gh) + *(6, Fyh, Gp) + *(6, Gp, xi_yy) + *(12, B1h, Gh), sinhG) + *(*(-12, B2h ^ 2) + *(-12, phih ^ 2) + *(-6, B2s) + *(-6, B1h ^ 2) + *(-6, Gh ^ 2) + *(3, Fyp ^ 2) + *(6, B1s) + *(-6, B1p + *(-1, B2p), Fyh + xi_yy) + *(6, B1h, B2h), coshG)) + *(12, Sh ^ 2, coshG), expB2) + *(*(S, *(*(-24, Sb) + *(-24, B1t, St) + *(-24, B2t, St) + *(24, Fxt, Sp) + *(24, Sp, xi_xx), coshG) + *(-24, Gt, St, sinhG)) + *(S ^ 2, *(*(-12, B2t ^ 2) + *(-12, phit ^ 2) + *(-6, B1b) + *(-6, B2b) + *(-6, B1t ^ 2) + *(-6, Gt ^ 2) + *(3, Fxp ^ 2) + *(-6, B1t, B2t) + *(6, B1p + B2p, Fxt + xi_xx), coshG) + *(-1, *(6, Gb) + *(-6, Gp, Fxt + xi_xx) + *(6, Gt, B2t + *(2, B1t)), sinhG)) + *(12, St ^ 2, coshG), exp(B2 + *(2, B1))) + *(*(S ^ 2, *(*(12, Gc) + *(-6, Gp, Fxh + Fyt + *(2, xi_xy)) + *(6, Gh, B1t + B2t) + *(6, Gt, B2h + *(-1, B1h)), coshG) + *(*(12, B2c) + *(-6, B2p, Fxh + Fyt + *(2, xi_xy)) + *(-6, Fxp, Fyp) + *(12, Gh, Gt) + *(24, B2h, B2t) + *(24, phih, phit), sinhG)) + *(24, S, *(*(Gh, St) + *(Gt, Sh), coshG) + *(*(2, Sc) + *(B2h, St) + *(B2t, Sh) + *(-1, Sp, Fxh + Fyt + *(2, xi_xy)), sinhG)) + *(-24, Sh, St, sinhG), exp(B1 + B2))

    nothing
end

function xi_t_eq_coeff(vars::Tuple, ::Outer)
    (
        kappa, xi, xi_x, xi_y, xi_xx, xi_yy, xi_xy,
        B1   , B2   , G   , phi  , S    , Fx    , Fy    , Sd ,  B1d  , B2d  , Gd,  phid, A   ,
        B1p  , B2p  , Gp  , phip , Sp   , Fxp   , Fyp   , Sdp,  B1dp , B2dp , Gdp,       Ap  ,
        B1pp , B2pp , Gpp ,        Spp  , Fxpp  , Fypp  ,                                App ,
        B1_x , B2_x , G_x , phi_x, S_x  , Fx_x  , Fy_x  , Sd_x, B1d_x, B2d_x, Gd_x,      A_x ,
	B1_y , B2_y , G_y , phi_y, S_y  , Fx_y  , Fy_y  , Sd_y, B1d_y, B2d_y, Gd_y,      A_y ,
        B1p_x, B2p_x, Gp_x,        Sp_x , Fxp_x , Fyp_x ,                                Ap_x,
        B1p_y, B2p_y, Gp_y,        Sp_y , Fxp_y , Fyp_y ,                                Ap_y,
                                                  Fy_xx ,                                A_xx,
                                          Fx_yy ,                                        A_yy,
                                          Fx_xy , Fy_xy ,                                A_xy,
    ) = vars

    @tilde_outer("B1")
    @tilde_outer("B2")
    @tilde_outer("G")
    @tilde_outer("phi")
    @tilde_outer("S")
    @tilde_outer("Fx")
    @tilde_outer("Fy")
    @tilde_outer("Sd")
    @tilde_outer("B1d")
    @tilde_outer("B2d")
    @tilde_outer("Gd")
    @tilde_outer("A")

    @hat_outer("B1")
    @hat_outer("B2")
    @hat_outer("G")
    @hat_outer("phi")
    @hat_outer("S")
    @hat_outer("Fx")
    @hat_outer("Fy")
    @hat_outer("Sd")
    @hat_outer("B1d")
    @hat_outer("B2d")
    @hat_outer("Gd")
    @hat_outer("A")

    @bar_outer("A")
    @bar_outer("Fy")

    @star_outer("A")
    @star_outer("Fx")

    @tilde_outer("B1p")
    @tilde_outer("B2p")
    @tilde_outer("Gp")
    @tilde_outer("Sp")
    @tilde_outer("Fxp")
    @tilde_outer("Fyp")
    @tilde_outer("Ap")

    @hat_outer("B1p")
    @hat_outer("B2p")
    @hat_outer("Gp")
    @hat_outer("Sp")
    @hat_outer("Fxp")
    @hat_outer("Fyp")
    @hat_outer("Ap")

    @cross_outer("A")
    @cross_outer("Fx")
    @cross_outer("Fy")

    x0 = exp(B1 + B2)
    x1 = cosh(G)
    x2 = S*x1
    x3 = sinh(G)
    x4 = exp(B2)
    x5 = x3*x4
    x6 = 2*x5
    x7 = -B1
    x8 = exp(B2 + x7)
    x9 = B2h*S
    x10 = Fyp*S
    x11 = Sp*xi_y
    x12 = 2*B2p
    x13 = S*x12
    x14 = x13*xi_y
    x15 = exp(B1)
    x16 = S*x15
    x17 = Gt*x16
    x18 = 2*Fx
    x19 = Gp*x16
    x20 = x18*x19
    x21 = 2*xi_x
    x22 = x19*x21
    x23 = Fy + xi_y
    x24 = Gp*x23
    x25 = S*(Gh + 2*x24)
    x26 = B1t*S
    x27 = 2*B1p
    x28 = x12 + x27
    x29 = Fx*S
    x30 = S*x27
    x31 = B2t*S
    x32 = Fxp*S
    x33 = -Fx*Sp - Sp*xi_x + St + x13*xi_x + x31 + x32
    x34 = -Sh
    x35 = B1h*S
    x36 = 6*Sp
    x37 = Fyp*x23
    x38 = x1*x4
    x39 = x37*x38
    x40 = x23^2
    x41 = x38*x40
    x42 = 3*Spp
    x43 = 3*B2p
    x44 = Sp*x41
    x45 = 3*Sp
    x46 = x40*x45*x5
    x47 = Fx + xi_x
    x48 = Fxp*x47
    x49 = exp(2*B1 + B2)
    x50 = x1*x49
    x51 = x48*x50
    x52 = x47^2
    x53 = x50*x52
    x54 = Gp*x3
    x55 = x45*x49*x52
    x56 = Sp*x53
    x57 = B1p + B2p
    x58 = x15*(x1*(B1t + B2t + x47*x57) + x3*(Fx*Gp + Gp*xi_x + Gt))
    x59 = B2p*x23
    x60 = -B2h - x59
    x61 = x1*(-Gh - x24) + x3*x60
    x62 = x58 + x61
    x63 = Sp*x23
    x64 = Sh + x63
    x65 = x3*x64
    x66 = Sp*x47
    x67 = St + x66
    x68 = x1*x15
    x69 = x67*x68
    x70 = -x65 + x69
    x71 = S*x62 + x70
    x72 = x0*x71
    x73 = 2*xi_xy
    x74 = Fxp*x23
    x75 = x0*x3
    x76 = 6*S
    x77 = Fyp*x47
    x78 = x0*x1
    x79 = x15*(Sd*x76 + x5*(2*Fyt + x73 + 2*x77) - x78*(2*Fxt + 2*x48 + 2*xi_xx))
    x80 = -x38*(2*Fyh + 2*x37 + 2*xi_yy) + x75*(2*Fxh + x73 + 2*x74) + x79
    x81 = St - 2*x66
    x82 = x15*x3
    x83 = x81*x82
    x84 = B1h + B1p*x23
    x85 = B2h + x59
    x86 = x1*x85
    x87 = Gh + x24
    x88 = x3*x87
    x89 = B2p*x47 + B2t
    x90 = x82*x89
    x91 = Gp*x47 + Gt
    x92 = x68*x91
    x93 = x1*x84 - x86 - x88 + x90 + x92
    x94 = S*x93 + x1*(x34 - x63) + x83
    x95 = x4*x94
    x96 = Fyp*x95
    x97 = x23*x95
    x98 = Fyh + xi_yy
    x99 = x37 + x98
    x100 = x38*x99
    x101 = 2*Gp
    x102 = x5*x99
    x103 = Fypp*x23
    x104 = Fxpp*x23
    x105 = Fxh + xi_xy
    x106 = x105 + x74
    x107 = x106*x78
    x108 = x106*x75
    x109 = Fypp*x47 + Fypt
    x110 = Fxpp*x47 + Fxpt
    x111 = 2*x78
    x112 = Fyt + xi_xy
    x113 = x112 + x77
    x114 = x113*x5
    x115 = x113*x38
    x116 = Fxt + xi_xx
    x117 = x116 + x48
    x118 = x117*x75
    x119 = Sph + Spp*x23
    x120 = Gp*x1
    x121 = Spp*x47
    x122 = Spt + x121
    x123 = x67*x82
    x124 = B2ph + B2pp*x23
    x125 = Gph + Gpp*x23
    x126 = B1pp*x47 + B1pt
    x127 = B2pp*x47 + B2pt
    x128 = Gpp*x47 + Gpt
    x129 = B1p*x47 + B1t
    x130 = x0*(x18 + x21)
    x131 = x68*x81
    x132 = Fxp*Sp
    x133 = B1ph + B1pp*x23
    x134 = x1*x89
    x135 = Gp*x15
    x136 = x3*x91
    x137 = x4*(2*Fy + 2*xi_y)
    x138 = S*x80 + x137*x94 + x41*x45 + x45*x53
    x139 = -x130*x71 + x138
    x140 = exp(x7)
    x141 = x140/2
    x142 = A*B1p
    x143 = B1d - x142/2
    x144 = x15*(x1*x129 + x134 + x136)
    x145 = x144 + x61
    x146 = S*x145 + x70
    x147 = A*Sp
    x148 = 3*x147
    x149 = A*Gp
    x150 = Gd - x149/2
    x151 = 3*B2d
    x152 = A*B2p
    x153 = x151 - 3*x152/2
    x154 = A/2
    x155 = Ap*Sp
    x156 = 3*Sdp - 3*Spp*x154 - 3*x155/2
    x157 = x150*x3
    x158 = Sd - x147/2
    x159 = x0*x146
    x160 = 2*B2d - x152
    x161 = 2*B1d
    x162 = -x142 + x160 + x161
    x163 = x1*x150
    x164 = Ah + Ap*x23
    x165 = x164/2
    x166 = Sdh + Sdp*x23 - Sp*x165 - x119*x154
    x167 = Ap*x47 + At
    x168 = x167/2
    x169 = Sdp*x47 + Sdt - Sp*x168 - x122*x154
    x170 = B2dh + B2dp*x23 - B2p*x165 - x124*x154
    x171 = Gdh + Gdp*x23 - Gp*x165 - x125*x154
    x172 = B2dp*x47 + B2dt - B2p*x168 - x127*x154
    x173 = Gdp*x47 + Gdt - Gp*x168 - x128*x154
    x174 = 3*A/2
    x175 = x15*x150
    x176 = 2*Gd
    x177 = -x149 + x176
    x178 = 16*S*Sd
    x179 = 3*Sd
    x180 = x1^2
    x181 = 3*B1d*x180
    x182 = S^2
    x183 = 4*x182
    x184 = 4*S
    x185 = 4*Sh
    x186 = 8*x9
    x187 = 4*x32
    x188 = 4*x10
    x189 = Fxs*x184 - Fyc*x184 - x105*x185 + x105*x186 + x112*x185 - x112*x186 + x112*x188 - x187*x98
    x190 = 4*St
    x191 = 8*x31
    x192 = Fxc*x184 - Fyb*x184 - x105*x187 - x105*x190 + x105*x191 + x112*x190 - x112*x191 + x116*x188
    x193 = 6*B2d
    x194 = 8*phid
    x195 = 2*B1dh
    x196 = Ah*B1p
    x197 = x161*x180
    x198 = 2*G
    x199 = sinh(x198)
    x200 = B1d*x199
    x201 = 2*x200
    x202 = S^3
    x203 = 2*x202
    x204 = 4*Gd
    x205 = 2*B1dt
    x206 = S*x199
    x207 = At*B1p
    x208 = B1d*cosh(x198)
    x209 = Sh*x178/2 - x15*x182*(-At*S*x101 - 12*Gd*St - Gdt*x184 + Gt*x184*x208 + 6*St*x200 + x201*x26 - x204*x26 + x205*x206 + x206*x207)/2 + x183*(Ah*Sp - Fyp*x179 - 4*Sdh - Sh*x151 + Sh*x181)/2 - x189*x5/2 + x192*x78/2 - x203*(Ah*B2p - Ap*Fyp + 2*Aph + B1h*x197 + 2*B2dh + B2h*x193 + Gh*x176 - Gh*x201 + phih*x194 - x180*x195 - x180*x196)/2
    x210 = 1/x202
    x211 = x210*x23
    x212 = x141*(6*Sh*x182*(x176 + x200) - x15*(-St*x178 + x183*(-At*Sp + Fxp*x179 + 4*Sdt + St*x151 + St*x181) + x203*(-Ap*Fxp + 2*Apt + At*B2p + B1t*x197 + 2*B2dt + B2t*x193 + Gt*x176 + Gt*x201 + phit*x194 + x180*x205 + x180*x207)) - x189*x38 + x192*x75 - x202*(-Ah*x101 + B1h*x201 + B1h*x204 - 4*Gdh - 4*Gh*x208 - x195*x199 - x196*x199))
    x213 = 3*S
    x214 = Ah*x1
    x215 = Ah*x3
    x216 = x1*x182
    x217 = At*x1
    x218 = x182*x217
    x219 = At*x3
    x220 = Ap*x3
    x221 = x210*x47

    axx = -x0*x2

    axy = S*x6

    ayy = -x2*x8

    bx = x4*(x1*(-x15*(x26 + x28*x29 + x30*xi_x + x33) + x25) + x3*(-Fy*Sp + Fy*x13 + Sh + x10 - x11 + x14 - x17 - x20 - x22 + x9))

    by = x8*(x1*(Fy*(Sp - x13 + x30) - x10 + x11 - x14 + x17 + x20 + x22 + x30*xi_y + x34 + x35 - x9) + x3*(x15*(x12*x29 + x33) - x25))

    cc = x141*(-B1p*x139 - 2*Fxp*x72 + Gp*x46 + S*(B1p*x79 - x100*x12 - x101*x102 + x101*x107 + x108*x28 + x15*(Sd*x36 + Sdp*x76 + x101*x115 - x101*x118 + x109*x6 - x110*x111 - x111*x117*x57 + x114*x12) - x38*(2*Fyph + 2*x103) + x75*(2*Fxph + 2*x104)) + Sp*x80 + x12*x97 - x130*(B1p*x69 + Gp*x123 + S*(B1p*x58 - Gp*x86 - Gp*x88 - x1*x125 - x124*x3 + x15*(x1*x126 + x1*x127 + x120*x91 + x128*x3 + x129*x54 + x54*x89)) + Sp*x62 - x119*x3 - x120*x64 + x122*x68) + x137*(B1p*x83 + Gp*x131 - Gp*x65 + S*(B1p*x90 + B1p*x92 - x1*x124 + x1*x133 - x120*x87 - x125*x3 + x127*x82 + x128*x68 + x134*x135 + x135*x136 + x54*x84 - x54*x85) + Sp*x93 - x1*x119 + x82*(Spt - 2*x121 - 3*x132)) - x28*x47*x72 + x36*x39 + x36*x51 + x41*x42 + x42*x53 + x43*x44 + x54*x55 + x56*(6*B1p + x43) + 2*x96)

    SS = kappa*x139*x141 - x141*x143*(-x130*x146 + x138) + x141*(A*Fxp*x159 - A*x96 + S*(-x100*x160 - x102*x177 + x107*x177 + x108*x162 + x143*x79 + x15*(2*Sd*(-3*x147/2 + x179) + x114*x160 + x115*x177 - x117*x162*x78 - x118*x177 + x213*(-A*Sdp + x140*x210*(-S*x4*(Ap*x2*x98 - As*x2 - Gh*S*x215 - Sh*x214 + x214*x35 - x214*x9) - x0*(S*x3*(Ah*St + At*Sh) + x182*(2*Ac*x3 + B2h*x219 + B2t*x215 + Gh*x217 + Gt*x214 - x105*x220 - x112*x220)) - x15*x202*(-Ap*x179 + B1d^2*S*x180 + B2d^2*x213 + Gd^2*S + phid^2*x184) - x49*(-Ab*x216 + Ap*x116*x216 - At*St*x2 - B1t*x218 - B2t*x218 - Gt*x182*x219) + (Fxh - Fyt)^2*exp(B1 + 2*B2))/3) + x5*(-A*x109 - Fyp*x167 + x209*x221) - x78*(-A*x110 - Fxp*x167 + x212*x221)) - x38*(-A*(Fyph + x103) - Fyp*x164 + x209*x211) + x75*(-A*(Fxph + x104) - Fxp*x164 + x211*x212)) - x130*(S*(-x1*x171 + x143*x144 + x15*(x1*x172 + x1*(B1dp*x47 + B1dt - B1p*x168 - x126*x154) + x129*x157 + x157*x89 + x163*x91 + x173*x3) - x150*x88 + x163*x60 - x170*x3) + x123*x150 + x143*x69 + x145*x158 - x163*x64 - x166*x3 + x169*x68) + x137*(S*(-x1*x170 + x1*(B1dh + B1dp*x23 - B1p*x165 - x133*x154) + x134*x175 + x136*x175 + x143*x90 + x143*x92 + x157*x84 - x157*x85 - x163*x87 - x171*x3 + x172*x82 + x173*x68) - x1*x166 + x131*x150 + x143*x83 - x150*x65 + x158*x93 + x82*(x132*x174 + x169 - x47*(3*Sdp - Spp*x174 - 3*x155/2))) - x148*x39 - x148*x51 + x150*x46 + x153*x44 + x156*x41 + x156*x53 + x157*x55 + x158*x80 - x159*x162*x47 + x160*x97 + x56*(6*B1d - 3*x142 + x153))

    return axx, ayy, axy, bx, by, cc, SS
end


text:
abstract type State end

#### Required Interfaces for State #####
""" queued_count(state)::Int

Returns the total number of jobs in all queues of the system """
function queued_count end

""" transit_count(state)::Int

Returns the total number of jobs in transit between queues in the system """
function transit_count end

""" num_in_queue(q::Int, state)::Int

Returns the number of jobs in the queue with given index """
function num_in_queue end

""" new_job(q, time, state)

Adds a new job to the system in teh given queue
"""
function new_job end

""" new_transit(time, transit_time, state)

Adds a new job to the system in transit
"""
function new_transit end

""" pop_queue(q, time, state)

Removes the first job in the given queue from the system
"""
function pop_queue end

""" queue_to_transit(q, time, transit_time, state)

Moves the first job in the given queue into transit
"""
function queue_to_transit end

""" pop_transit(time, state)

Removes the first job in transit from the system
"""
function pop_transit end

""" transit_to_queue(q, time, state)

Adds the first job in transit to the given queue
"""
function transit_to_queue end

""" update_transit(time, transit_time, state)

Changes the destination time of the first job in transit
"""
function update_transit end

""" failed_arrival(time, state)

Changes the destination time of the first job in transit
"""
function failed_arrival end

""" clear_left(state)

Clears the left_system record
"""
function clear_left end

@with_kw struct NetworkParameters
    L::Int #Number of queues
    gamma_scv::Float64 #This is constant for all scenarios at 3.0
    λ::Float64 #arrival rate
    η::Float64 #This is assumed constant for all scenarios at 4.0
    μ_vector::Vector{Float64} #service rates
    P::Matrix{Float64} #routing matrix
    Q::Matrix{Float64} #overflow matrix
    p_e::Vector{Float64} #external arrival distribution
    K::Vector{Int} #Queue capacity, -1 means infinity
end

""" Constructor to easily change arrival rate"""
NetworkParameters(params::NetworkParameters, λ::Real) = NetworkParameters(
    L = params.L,
    gamma_scv = params.gamma_scv,
    λ = convert(Float64, λ),
    η = params.η,
    μ_vector = copy(params.μ_vector),
    P = copy(params.P),
    Q = copy(params.Q),
    p_e = copy(params.p_e),
    K = copy(params.K))

"""
Returns a Gamma distribution with desired rate (inverse of shape) and SCV.
"""
rate_scv_gamma(desired_rate::Float64, desired_scv::Float64) = Gamma(1/desired_scv, desired_scv/desired_rate)

""" Generates next arrival time from state """
next_arrival_time(state::State)::Float64 = rand(state.arrival_distr)

""" Generates next service time for the given queue from state """
next_service_time(state::State, q::Int)::Float64 = rand(state.service_distrs[q])

""" Generates travel time between queues from state """
travel_time(state::State)::Float64 = rand(state.travel_distr)

""" Randomly selects an entry queue index using the state weights """
function get_entry_queue(state::State)::Int
	prob = rand()
	for i in 1:length(state.params.p_e)
		prob -= state.params.p_e[i]
		prob <= 0 && return i
	end
    error("Entry probabilities do not sum to 1")
end

"""
Randomly selects the next queue for the job from the given routing or overflow matrix weights
Returns -1 corresponding to exiting the system.
"""
function get_next_queue(q::Int, state::State, mode::Symbol)::Int
    (mode ∉ [:routing, :overflow]) && error("$mode is not a valid mode for queue selection")
    # Gets the probability row vector corresponding to the current queue
    next_probs = (mode == :routing) ? state.params.P[q,:] :  state.params.Q[q,:]

    prob = rand()
	for i in 1:length(next_probs)
		prob -= next_probs[i]
		prob <= 0 && return i
	end
    return -1
end

""" Returns total number of jobs in the system """
total_count(state::State)::Int = queued_count(state) + transit_count(state)

""" Returns whether the number of jobs in specified queue is equal to queue capacity """
is_full(q::Int, state::State)::Bool = (state.params.K[q] != -1) && num_in_queue(q, state) == state.params.K[q]


text:
@testset "Statistics" begin

    @testset "Statistics: Ordinary Least Squares" begin
        a = [10, 14,17,21,20, 18, 42, 51, 77, 11, 91]
        b = [0.11,0.7,0.2,0.19,0.09,0.8, 0.71,0.1, 0.6,0.3,0.81]
        OLSbeta(a, b)

        df = DataFrame(a=a, b=b)
        ols = GLM.lm(@formula(a ~ b), df)
        coef(ols) # print coefficients

        # Predict
        predict(ols) # predict the model
    end

    @testset "Statistics: Pearson Correlation" begin
        a = 1:10
        b = 1:10
        @test pearson_correlation(a, b) == 1

        x = [12,15,17,21,18]
        y = [6,9,-1,3,1]

        @test pearson_correlation(x, y) ==  -0.5201361507569889
    end

    @testset "Statistics: Variance" begin
        a = 1:10

        @test variance(a) == 9.166666666666666

        a = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

        @test variance(a) == 9.166666666666666

        a = [12, 15, 17, 20]

        @test variance(a) == 11.333333333333334

    end

end

max_stars_repo_name: ErikQQY/Julia

max_issues_repo_name: ErikQQY/Julia

max_forks_repo_name: ErikQQY/Julia


text:
"""
    ror(x, k)

rotate x to the right (msbs toward lsb) by k bits (lsbs wrap to msbs)
""" ror

"""
    rol(x, k)

rotate x to the left (lsbs toward msb) by k bits (msbs wrap to lsbs)
""" rol

ror(x::UInt32, k::Int) = (x >>> (0x1f & k)) | (x <<  (0x1f & -k))
rol(x::UInt32, k::Int) = (x <<  (0x1f & k)) | (x >>> (0x1f & -k))
ror(x::UInt64, k::Int) = (x >>> (0x3f & k)) | (x <<  (0x3f & -k))
rol(x::UInt64, k::Int) = (x <<  (0x3f & k)) | (x >>> (0x3f & -k))

ror(x::Int32, k::Int) = (x >>> (0x1f & k)) | (x <<  (0x1f & -k))
rol(x::Int32, k::Int) = (x <<  (0x1f & k)) | (x >>> (0x1f & -k))
ror(x::Int64, k::Int) = (x >>> (0x3f & k)) | (x <<  (0x3f & -k))
rol(x::Int64, k::Int) = (x <<  (0x3f & k)) | (x >>> (0x3f & -k))


text:
# GGGenome
# ======
#
# APIs for GGGenome.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioServices.jl/blob/master/LICENSE.md

"""
GGGenome module.

The `GGGenome` module provides a programming interface to GGGenome,
an ultrafast DNA sequence search service, hosted by DBCLS.
Two functions are exported from this module:
1. `gggenome`: DNA sequence search
2. `gggenome_dblist`: List of available databases
See "GGGenome Help" (https://gggenome.dbcls.jp/en/help.html) for more details.
"""
module GGGenome

export
    gggsearch,
    gggdbs

import HTTP

const baseURL = "https://gggenome.dbcls.jp/"
const dblistURL = "https://raw.githubusercontent.com/meso-cacase/GGGenome/master/DBlist.pm"

# APIs of GGGenome
# -------------------

"""
    gggsearch(query::AbstractString;
              db="hg19", k=0, strand=nothing,
              format="html", timeout=5,
              output=nothing, show_url=false)

Retrieve results of gggenome search of a query sequence.

# Arguments
## Required
- `query::String`: Nucleotide sequence, case insensitive.

## Optional
- `db::String`: Target database name. hg19 if not specified. Full list of databases: https://gggenome.dbcls.jp/en/help.html#db_list
- `k::Integer`: Maximum number of mismatches/gaps. 0 if not specified.
- `strand::String`: '+' ('plus') or '-' ('minus') to search specified strand only.
- `format::String`: [html|txt|csv|bed|gff|json]. html if not specified.
- `timeout::Real`: Maximum time allowed for a query.
- `output::String`: If "toString", a `String` object is returned. If
  "extractTopHit", a `String` object containing only top hit is returned
  (Currently, only works with format="txt"). Otherwise, a
  `HTTP.Messages.Response` object is returned.
- `show_url::Bool`: If true, print URL of REST API.
"""
function gggsearch(query::AbstractString;
                    db::AbstractString="hg19", k::Integer=0, strand=nothing, format::AbstractString="html",
                    timeout::Real=5, output=nothing, show_url::Bool=false)
    # Check parameters
    if !checkQuery(query)
        throw(ArgumentError("`query` must be consisted of [A|C|G|T|U|N|R|Y|K|M|S|W|B|D|H|V]"))
    end
    if format ∉ ("html", "txt", "csv", "bed", "gff", "json")
        throw(ArgumentError("`format` must be [html|txt|csv|bed|gff|json]"))
    end
    if strand != nothing && strand ∉ ("+", "-")
        throw(ArgumentError("`output` must be \"+\" or \"-\""))
    end
    if output != nothing && output ∉ ("toString", "extractTopHit")
        throw(ArgumentError("`output` must be \"toString\" or \"extractTopHit\""))
    end

    # Generate URL
    url = baseURL
    url *= db * "/"
    url *= string(k) * "/"
    if strand isa AbstractString
        url *= strand * "/"
    end
    url *= query
    url *= "." * format

    # Show URL
    if show_url
        println(url)
    end

    # Request
    res = HTTP.request("GET", url, timeout=timeout)

    # Output
    if output isa AbstractString
        if output == "toString"
            return gggenomeToString(res)
        elseif output == "extractTopHit" && format == "txt"
            return extractTopHit(gggenomeToString(res))
        else
            # unreachable
            @assert false
        end
    else
        return res
    end
end

function checkQuery(query::AbstractString)
    return all(c -> c ∈ "NRYKMSWBDHVACGTU", query)
end

function gggenomeToString(res::HTTP.Messages.Response)
    return String(res.body)
end

function extractTopHit(res_str::String)
    topResult = "No hit"
    header = ""
    for ln in split(chomp(res_str), "\n")
        if startswith(ln, '#')
            header *= String(ln) * "\n"
        else
            topResult = header * String(ln)
            break
        end
    end
    return topResult
end

"""
    gggdbs()

Retrieve full list of available databases.
Full list of databases: https://gggenome.dbcls.jp/en/help.html#db_list.
"""
function gggdbs()
    res = HTTP.request("GET", dblistURL)
    arr = split(String(res.body), "\n")
    index_l = 0
    index_r = 0
    for i in 1:length(arr)
        if startswith("<<'--EOS--' ;", arr[i])
            index_l = i + 1
        elseif startswith("--EOS--", arr[i])
            index_r = i - 1
            break
        end
    end
    return arr[index_l:index_r]
end

end

max_stars_repo_name: BioJulia/BioServices

max_issues_repo_name: BioJulia/BioServices

max_forks_repo_name: BioJulia/BioServices


text:
import LinearAlgebra.BLAS.@blasfunc
import LinearAlgebra: chkstride1, BlasInt
import LinearAlgebra.LAPACK.chklapackerror

export QuotientSpace

struct QuotientSpace{S,O<:Operator,DD,T,RT} <: Space{DD,T}
    space::S
    bcs::O
    A::Matrix{T}
    x::Vector{T}
    b::Vector{T}
    U::Matrix{T}
    Σ::Diagonal{T}
    VT::Matrix{T}
    work::Vector{T}
    iwork::Vector{BlasInt}
    rwork::Vector{RT}
    info::Ref{BlasInt}
    function QuotientSpace{S,O,DD,T,RT}(space::S, bcs::O) where {S,O,DD,T,RT}
        n = size(bcs, 1)
        @assert isfinite(n)
        A = zeros(T, n, n)
        x = zeros(T, n)
        b = zeros(T, n)
        U = zeros(T, n, n)
        Σ = Diagonal(zeros(T, n))
        VT = zeros(T, n, n)
        work = Vector{T}(undef,max((3n+7)*n,68))
        iwork = Vector{BlasInt}(undef,8*n)
        rwork = Vector{RT}(undef,(5*n+7)*n)
        info = Ref{BlasInt}()
        new{S,O,DD,T,RT}(space, bcs, A, x, b, U, Σ, VT, work, iwork, rwork, info)
    end
end

QuotientSpace(sp::Space, bcs::Operator) = QuotientSpace{typeof(sp), typeof(bcs), domaintype(sp), rangetype(sp), real(rangetype(sp))}(sp, bcs)

QuotientSpace(bcs::Operator) = QuotientSpace(domainspace(bcs), bcs)

domain(QS::QuotientSpace) = domain(QS.space)
dimension(QS::QuotientSpace) = dimension(QS.space)

Conversion(Q::QuotientSpace{SP}, S::SP) where {SP<:Space} = ConcreteConversion(Q, S)

bandwidths(C::ConcreteConversion{QuotientSpace{SP,O,DD,T,RT},SP}) where {SP,O,DD,T,RT} = 0, size(C.domainspace.A, 1)

function getindex(C::ConcreteConversion{QuotientSpace{SP,O,DD,T,RT},SP}, i::Integer, j::Integer) where {SP,O,DD,T,RT}
    sp = domainspace(C)
    n = size(sp.A, 1)
    A = sp.A
    x = sp.x
    b = sp.b
    B = sp.bcs
    U = sp.U
    Σ = sp.Σ
    VT = sp.VT
    work = sp.work
    iwork = sp.iwork
    rwork = sp.rwork
    info = sp.info

    if i == j
        one(T)
    elseif i ≤ n && j ≤ n
        zero(T)
    elseif j-n ≤ i < j
        for jj = 1:n, ii = 1:n
            A[ii,jj] = B[ii,j+jj-n-1]
        end
        for ii = 1:n
            b[ii] = -B[ii,j]
        end
        in_place_gesdd!(A, U, Σ.diag, VT, work, iwork, rwork, info)
        mul!(x, transpose(U), b)
        mul!(b, pinv!(Σ), x)
        mul!(x, transpose(VT), b)
        x[i+n-j+1]
    else
        zero(T)
    end
end

function pinv!(D::Diagonal{T}) where T
    d = D.diag
    @inbounds @simd for i = 1:length(d)
        isfinite(inv(d[i])) ? d[i]=inv(d[i]) : d[i]=zero(T)
    end
    D
end

for (gesdd, elty, relty) in ((:dgesdd_,:Float64,:Float64),
                      (:sgesdd_,:Float32,:Float32),
                      (:zgesdd_,:ComplexF64,:Float64),
                      (:cgesdd_,:ComplexF32,:Float32))
    @eval begin
        #    SUBROUTINE DGESDD( JOBZ, M, N, A, LDA, S, U, LDU, VT, LDVT, WORK,
        #                   LWORK, IWORK, INFO )
        #*     .. Scalar Arguments ..
        #      CHARACTER          JOBZ
        #      INTEGER            INFO, LDA, LDU, LDVT, LWORK, M, N
        #*     ..
        #*     .. Array Arguments ..
        #      INTEGER            IWORK( * )
        #      DOUBLE PRECISION   A( LDA, * ), S( * ), U( LDU, * ),
        #                        VT( LDVT, * ), WORK( * )
        function in_place_gesdd!(A::Matrix{$elty}, U::Matrix{$elty}, S::Vector{$elty}, VT::Matrix{$elty}, work::Vector{$elty}, iwork::Vector{BlasInt}, rwork::Vector{$relty}, info::Ref{BlasInt})
            chkstride1(A)
            m, n   = size(A)
            minmn  = min(m, n)
            job = 'A'
            lwork  = BlasInt(-1)
            cmplx  = $elty != $relty
            for i = 1:2
                if cmplx
                    ccall((@blasfunc($gesdd), liblapack), Nothing,
                          (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                           Ref{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}),
                          job, m, n, A,
                          max(1,stride(A,2)), S, U, max(1,stride(U,2)),
                          VT, max(1,stride(VT,2)), work, lwork,
                          rwork, iwork, info)
                else
                    ccall((@blasfunc($gesdd), liblapack), Nothing,
                          (Ref{UInt8}, Ref{BlasInt}, Ref{BlasInt}, Ptr{$elty},
                           Ref{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{$elty}, Ref{BlasInt}, Ptr{$elty}, Ref{BlasInt},
                           Ptr{BlasInt}, Ptr{BlasInt}),
                          job, m, n, A,
                          max(1,stride(A,2)), S, U, max(1,stride(U,2)),
                          VT, max(1,stride(VT,2)), work, lwork,
                          iwork, info)
                end
                chklapackerror(info[])
                if i == 1
                    # Work around issue with truncated Float32 representation of lwork in
                    # sgesdd by using nextfloat. See
                    # http://icl.cs.utk.edu/lapack-forum/viewtopic.php?f=13&t=4587&p=11036&hilit=sgesdd#p11036
                    # and
                    # https://github.com/scipy/scipy/issues/5401
                    lwork = round(BlasInt, nextfloat(real(work[1])))
                    #work = Vector{$elty}(lwork)
                end
            end
        end
    end
end


text:
function Base.round(::Type{Bit{T}}, x::T) where {T <: Integer}
    x == zero(T) ? Bit(zero(T)) : Bit(one(T))
end

function Base.round(::Type{Bit{T}}, x::T; threshold::T=T(1e-3)) where {T <: AbstractFloat}
    x < threshold ? Bit(zero(T)) : Bit(one(T))
end

function Base.round(::Type{Spin{T}}, x::T) where T
    x > 0 ? Spin(one(T)) : Spin(-one(T))
end

function Base.round(::Type{Half{T}}, x::T) where {T <: AbstractFloat}
    x > 0 ? Half(T(0.5)) : Half(-T(0.5))
end


text:
"""
    random_micrm_params(N,M,leakage)

    Generates a parameter dictionary for micrm system with `N` consumers and `M` resources with random parameters. All values are drawn from uniform random distributions apart from uptake (u) and leakage (l) which are drawn from dirchlet distributions. This means overall uptake and leakage sum to 1 and the value set by the leakage parameter respectively.
"""
function random_micrm_params(N,M,leakage)
    #uptake
    du = Distributions.Dirichlet(M,1.0)
    u = rand(du, N)'

    #cost term
    m = rand(N)

    #inflow + outflow
    ρ,ω = rand(M),rand(M)

    #leakage
    l = rand(du,M)' .* leakage

    return Dict(:N => N, :M => M, :u => u, :m => m, :ρ => ρ, :ω => ω, :l => l)
end


text:
mutable struct GFGenerateState
    trace::GFTrace
    constraints::Any
    score::Float64
    weight::Float64
    visitor::AddressVisitor
    params::Dict{Symbol,Any}
end

function GFGenerateState(constraints, params::Dict{Symbol,Any})
    GFGenerateState(GFTrace(), constraints, 0., 0., AddressVisitor(), params)
end

get_args_change(state::GFGenerateState) = nothing
get_addr_change(state::GFGenerateState, addr) = nothing
set_ret_change!(state::GFGenerateState, value) = begin end

function addr(state::GFGenerateState, dist::Distribution{T}, args, addr) where {T}
    visit!(state.visitor, addr)
    constrained = has_leaf_node(state.constraints, addr)
    if has_internal_node(state.constraints, addr)
        lightweight_got_internal_node_err(addr)
    end
    local retval::T
    if constrained
        retval = get_leaf_node(state.constraints, addr)
    else
        retval = random(dist, args...)
    end
    score = logpdf(dist, retval, args...)
    call = CallRecord(score, retval, args)
    state.trace = assoc_primitive_call(state.trace, addr, call)
    state.trace.has_choices = true
    state.score += score
    if constrained
        state.weight += score
    end
    retval
end

function addr(state::GFGenerateState, gen::Generator{T,U}, args, addr, args_change) where {T,U}
    visit!(state.visitor, addr)
    if has_internal_node(state.constraints, addr)
        constraints = get_internal_node(state.constraints, addr)
    else
        constraints = EmptyAssignment()
    end
    (trace::U, weight) = generate(gen, args, constraints)
    call::CallRecord = get_call_record(trace)
    state.trace = assoc_subtrace(state.trace, addr, trace)
    state.trace.has_choices |= has_choices(trace)
    state.score += call.score
    state.weight += weight
    call.retval::T
end

splice(state::GFGenerateState, gf::GenFunction, args::Tuple) = exec(gf, state, args)

function generate(gen::GenFunction, args, constraints)
    state = GFGenerateState(constraints, gen.params)
    retval = exec(gen, state, args)
    # TODO add return type annotation for gen
    call = CallRecord{Any}(state.score, retval, args)
    state.trace.call = call
    (state.trace, state.weight)
end

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
const should_not_warn = Set{Function}()

"""
    @should_not_warn function foo(x)
      ...
    end

Add `foo` to the list of functions in which no warnings may occur (checkd by `@check`).
"""
macro should_not_warn(expr)
  quote
    fun = $(esc(expr))
    push!(should_not_warn, fun)
    fun
  end
end

"""
    check(f::Function; nowarn=[], except=[], kwargs...)

Run Traceur on `f`, and throw an error if any warnings occur inside functions
tagged with `@should_not_warn` or specified in `nowarn`.

To throw an error if any warnings occur inside any functions, set
`nowarn=:all`.

To throw an error if any warnings occur inside any functions EXCEPT for a
certain set of functions, list the exceptions in the `except` variable,
for example `except=[g,h]`
"""
function check(f; nowarn=Any[], except=Any[], kwargs...)
  if !isempty(except) # if `except` is provided, we ignore the value of `nowarn`
    _nowarn = Any[]
    _nowarn_all = false
    _nowarn_allexcept = true
  elseif nowarn isa Symbol
    _nowarn = Any[]
    _nowarn_all = nowarn == :all
    _nowarn_allexcept = false
  else
    _nowarn = nowarn
    _nowarn_all = false
    _nowarn_allexcept = false
  end
  failed = false
  wp = warning_printer()
  result = trace(f; kwargs...) do warning
    ix = findfirst(warning.stack) do call
      _nowarn_all || call.f in should_not_warn || call.f in _nowarn || (_nowarn_allexcept && !(call.f in except))
    end
    if ix != nothing
      tagged_function = warning.stack[ix].f
      message = "$(warning.message) (called from $(tagged_function))"
      warning = Warning(warning.call, warning.line, message, warning.stack)
      wp(warning)
      failed = true
    end
  end
  @assert !failed "One or more warnings occured inside functions tagged with `@should_not_warn` or specified with `nowarn`"
  result
end

"""
    @check fun(args...) nowarn=[] except=[] maxdepth=typemax(Int)

Run Traceur on `fun`, and throw an error if any warnings occur inside functions
tagged with `@should_not_warn` or specified in `nowarn`.

To throw an error if any warnings occur inside any functions, set
`nowarn=:all`.

To throw an error if any warnings occur inside any functions EXCEPT for a
certain set of functions, list the exceptions in the `except` variable,
for example `except=[g,h]`
"""
macro check(expr, args...)
  quote
      check(() -> $(esc(expr)); $(map(esc, args)...))
    end
end


text:
include("../../Computation/Computation.jl")
import Main.Computation
using StaticArrays

"Mutlidimensional kernel operation"
function compute_convND(kernel, inputs, input_steps, kernel_shape, idx, ::Val{N}, ::Type{T}) where {T, N}
    input_index = idx
    kernel_index = 1
    index_vector = MArray{Tuple{N}, T}(undef)
    sum = 0.0
    dim = 1

    for i in 1:length(index_vector)
        index_vector[i] = 1
    end

    "Multidimensional Iteration. Adapted from HyperSphere.Hyperdimensional.Iteration"
    @label iter
    #Check to see if the index is less then the kernel bound
    if index_vector[dim] <= kernel_shape[end - dim + 1]
        if dim != N
            dim += 1
            @goto iter
        else
            for i in 1:kernel_shape[1]
                sum += kernel[kernel_index] * inputs[input_index]
                kernel_index += 1
                index_vector[end] += 1
                input_index += input_steps[end]
            end
        end

    end
    if dim != 1
        #Reset the input_index back to original position
        input_index -= (index_vector[dim] - 1) * input_steps[dim]
        index_vector[dim] = 1
        dim -= 1
        index_vector[dim] += 1
        input_index += input_steps[dim]
        @goto iter
    end

    return sum
end

"Compile Time Converting the linear indexes from the output_mat index to the input_mat index for mapping in kernel"
function convert_index_matrix(input_shape, output_shape)
    output_len = prod(output_shape)
    output_indexes = reshape(collect(1:output_len), output_shape...)
    input_indexes = reshape(collect(1:prod(input_shape)), input_shape...)

    for idx in CartesianIndices(output_indexes)
        output_indexes[idx] = input_indexes[idx]
    end

    return output_indexes
end

"Compile time experimental determination of the summation deltas between dimensions"
function calculate_input_steps(device, output_shape, input_index_mat, input_type, stride)
    steps = Array{input_type}(undef, length(output_shape))
    first_input_entry_index = fill(1, length(output_shape))

    for dim in 1:length(output_shape)
        init = input_index_mat[first_input_entry_index...]
        first_input_entry_index[dim] = stride + 1
        after = input_index_mat[first_input_entry_index...]
        first_input_entry_index[dim] = 1
        steps[dim] = after - init                              #Calculate the step change
    end

    return Computation.alloc(device, reverse(steps))
end

function test_device_conv_logic(device, input_matrix, kernel_matrix)
#START PRECOMPILE CODE
    stride = 1
    input_shape = size(input_matrix)
    kernel_shape = size(kernel_matrix)
    output_shape = input_shape .- kernel_shape .+ 1
    output_matrix = Computation.alloc(device, Int, output_shape...)

    input_type = Computation.get_min_int(prod(input_shape))
    kernel_type = Computation.get_min_int(prod(kernel_shape))

    #Flatten matrixes for far faster operation
    flat_outputs = reshape(output_matrix, :)
    flat_inputs = reshape(input_matrix, :)
    flat_kernel = reshape(kernel_matrix, :)

    #Determine the mathmatics to perform cartesian operations upon the linear indexes at compile time
    converted_indexes = convert_index_matrix(input_shape, output_shape)
    input_steps = calculate_input_steps(device, output_shape, converted_indexes, input_type, stride)
    converted_indexes = Computation.alloc(device, reshape(input_type.(converted_indexes), length(flat_outputs)))     #Flatten the indexes

    #Create Singleton arguments
    flat_kernel = Ref(flat_kernel)
    flat_inputs = Ref(flat_inputs)
    input_steps = Ref(input_steps)
    dim_size = Val(length(kernel_shape))
    kernel_shape = Ref(Computation.alloc(device, kernel_type.(kernel_shape)))
#END PRECOMPILE CODE

    @time broadcast!(compute_convND, flat_outputs, flat_kernel, flat_inputs, input_steps, kernel_shape, converted_indexes, dim_size, input_type)
    @time broadcast!(compute_convND, flat_outputs, flat_kernel, flat_inputs, input_steps, kernel_shape, converted_indexes, dim_size, input_type)

    return Array(output_matrix)
end

function test()
    println("Starting Device Conv Layer....")

    input_matrix = [[3,2,4,2,7,6] [8,0,2,1,7,8] [1,5,4,6,5,0] [5,4,1,7,5,6] [5,0,2,7,6,8] [1,1,1,1,1,1]]
    kernel_matrix = [[1,1,1] [0,0,0] [-1,-1,-1]]
    true_array = stride_conv_2D(input_matrix, kernel_matrix)
    #input_matrix = reshape(collect(1:125), 5, 5, 5)
    #kernel_matrix = reshape(collect(1:27), 3, 3, 3)
    #true_array = stride_conv_3D(input_matrix, kernel_matrix)

    for device in Computation.get_devices()
        println("\tChecking $device")
        test = test_device_conv_logic(device, Computation.alloc(device, input_matrix), Computation.alloc(device, kernel_matrix))
        @assert test == true_array "Failed $device Dense Layer Test"
        println("\tSuccesfully Completed $device")
    end
    println("Successfully Compeleted Device Conv Layer Test!")
end

function stride_conv_2D(input, filter)
    input_r, input_c = size(input)
    filter_r, filter_c = size(filter)

    result = zeros(Int, (input_r-filter_r) + 1, (input_c-filter_c) + 1)
    result_r, result_c = size(result)

    ir = 1
    ic = 1
    for i in 1:result_r
        for j in 1:result_c
            for k in 1:filter_r
                for l in 1:filter_c
                 #   (i == j && j == 1) && println(input[ic+l-1, ir+k-1])
                    result[i,j] += input[ir+k-1,ic+l-1]*filter[k,l]
                end
            end
            ic += 1
        end
        ir += 1
        ic = 1 # Return back to 1 after finish looping over column
    end

    return result
end

function stride_conv_3D(input, filter)
    input_d, input_r, input_c = size(input)
    filter_d, filter_r, filter_c = size(filter)

    result = zeros(Int, (input_d-filter_d) + 1, (input_r-filter_r) + 1, (input_c-filter_c) + 1)
    result_r, result_c, result_d = size(result)

    ir = 1
    ic = 1
    id = 1
    for d in 1:result_d
        for i in 1:result_r
            for j in 1:result_c
                for k in 1:filter_r
                    for l in 1:filter_c
                        for q in 1:filter_d
                            result[d,i,j] += input[id+q-1, ir+k-1,ic+l-1]*filter[q,k,l]
                        end
                    end
                end
                ic += 1
            end
            ir += 1
            ic = 1 # Return back to 1 after finish looping over column
        end
        id += 1
        ir = 1
    end
    return result
end

test()


text:
export permanent, veroneselift

using Combinatorics

function permanent(A, perms)
    n = LinearAlgebra.checksquare(A)
    sum(σ -> prod(i -> A[i,σ[i]], 1:n), perms)
end

function μ(combi)
    prod = 1
    last = 0
    cur = 1
    for i in combi
        if i == last
            cur += 1
        else
            prod *= factorial(cur)
            last = i
            cur = 1
        end
    end
    prod *= factorial(cur)
    return prod
end

"""
   veroneselift_explicit(A::AbstractMatrix, d::Integer)

Computes the Veronese lift of degree `d` of `A` using the explicit formula
```math
(A^{[d]})_{\\alpha\\beta} = \\frac{\\text{per}\\,A(\\alpha, \\beta)}{\\sqrt{\\mu(\\alpha)\\mu(\\beta)}}
```
given in equation (8) of [PJ08]. Note that this is a lot slower than
[`veroneselift`](@ref) and is implemented mainly for educational purpose.

* [PJ08] P. Parrilo and A. Jadbabaie.
*Approximation of the joint spectral radius using sum of squares*.
Linear Algebra and its Applications, Elsevier, **2008**, 428, 2385-2402
"""
function veroneselift_explicit(A::AbstractMatrix, d::Integer)
    d > 0 || throw(ArgumentError("`d` should be a positive integer, got $d"))

    n = LinearAlgebra.checksquare(A)
    N = binomial(n+d-1, d)

    Ad = Matrix{Float64}(undef, N, N)

    # Computing permutations is the most expensive part so we precompute it
    perms = collect(permutations(1:d))
    for (i,α) in enumerate(with_replacement_combinations(1:n, d))
        for (j,β) in enumerate(with_replacement_combinations(1:n, d))
            Ad[i, j] = permanent(A[α, β], perms) / sqrt(μ(α) * μ(β))
        end
    end

    return Ad
end

function transform_combination(α, m)
    β = zeros(Int, m)
    for i in α
        β[i] += 1
    end
    return β
end

"""
   veroneselift(A::AbstractMatrix, d::Integer)

Computes the Veronese lift of degree `d` of `A` using `DynamicPolynomials`.
"""
function veroneselift(A::AbstractMatrix, d::Integer)
    d > 0 || throw(ArgumentError("`d` should be a positive integer, got $d"))
    m, n = size(A)
    @polyvar x[1:n]
    Ax = A * x
    M = binomial(m + d - 1, d)
    N = binomial(n + d - 1, d)
    Ad = Matrix{Float64}(undef, M, N)
    df = factorial(d)
    scaling(m) = sqrt(df / prod(factorial, exponents(m)))
    X = monomials(x, d)
    # Scaling for the scaled monomial basis
    col_scaling = Float64[scaling(m) for m in X]
    if m == n
        Y = X
        row_scaling = col_scaling
    else
        # Only exponents matter
        @polyvar y[1:m]
        Y = monomials(y, d)
        row_scaling = Float64[scaling(m) for m in Y]
    end
    for (i, mono) in enumerate(Y)
        α = exponents(mono)
        Axα = prod(i -> Ax[i]^α[i], 1:m)
        Ad[i, :] = MP.coefficients(Axα, X) .* row_scaling[i] ./ col_scaling
    end
    return Ad
end