text:
### File declaring various reaction networks for the tests to be run on ###

#Declares the vectors which contains the various test sets.
reaction_networks_standard = Vector{ReactionSystem}(undef,10)
reaction_networks_hill = Vector{ReactionSystem}(undef,10)
reaction_networks_constraint = Vector{ReactionSystem}(undef,10)
reaction_network_constraints = Vector{Matrix{Int}}(undef, 10)
reaction_networks_real = Vector{ReactionSystem}(undef,4)
reaction_networks_weird = Vector{ReactionSystem}(undef,10)

### Standard reaction networks. ###
reaction_networks_standard[1] = @reaction_network begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
    (k1,k2), X2 ⟷ X1 + 2X3
    (k3,k4), X1 ⟷ X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end p1 p2 p3 k1 k2 k3 k4 d1 d2 d3

reaction_networks_standard[2] = @reaction_network begin
    mmr(X2,v1,K1), ∅ → X1
    mm(X1,v2,K2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 v2 K2 d

reaction_networks_standard[3] = @reaction_network begin
    mm(X2,v1,K1), ∅ → X1
    mm(X3,v2,K2), ∅ → X2
    (k1,k2), X1 ⟷ X3
    (k3,k4), X3 + X2 ⟷ X4 + X1
    d, (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 k1 k2 k3 k4 d

reaction_networks_standard[4] = @reaction_network begin
    mmr(X4,v1,K1), ∅ → X1
    mmr(X1,v2,K2), ∅ → X2
    mmr(X2,v3,K3), ∅ → X3
    mmr(X3,v4,K4), ∅ → X4
    (d1,d2,d3,d4), (X1,X2,X3,X4) → ∅
end v1 K1 v2 K2 v3 K3 v4 K4 d1 d2 d3 d4

reaction_networks_standard[5] = @reaction_network begin
    p, ∅ → X1
    (k1,k2), X1 ⟷ X2
    (k3,k4), X2 ⟷ X3
    (k5,k6), X3 ⟷ X4
    d, X4 → ∅
end p k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[6] = @reaction_network begin
    (p1,p2), ∅ → (X1,X2)
    (k1,k2), 2X1 ⟷ X3
    (k3,k4), X2 + X3 ⟷ 4X4
    (k5,k6), X4 + X1 ⟷ 2X3
    d, (X1,X2,X3,X4) → ∅
end p1 p2 k1 k2 k3 k4 k5 k6 d

reaction_networks_standard[7] = @reaction_network begin
    (p1,p2,p3), ∅ → (X1,X2,X3)
    (k1,k2), X1 + 2X2 ⟷ X4
    (mm(X3,v1,K1),k3), X4 ⟷ X5
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5) → ∅
end p1 p2 p3 k1 k2 k3 v1 K1 d1 d2 d3 d4 d5

reaction_networks_standard[8] = @reaction_network begin
    p, ∅ → 2X1
    k1, X1 → X2
    (k2, k3), X2 → X3
    d, X3 → ∅
end p k1 k2 k3 d

reaction_networks_standard[9] = @reaction_network begin
    (p1,p2,p3), ∅ ⟶ (X1,X2,X3)
    (d1,d2,d3), (X1,X2,X3) ⟶ ∅
    (k1,k2), X1 + X2 ⟷ X3
    (k3,k4), X3 ⟷ X4
    (k5,k6), X4 ⟷ X1 + X2
end p1 p2 p3 d1 d2 d3 k1 k2 k3 k4 k5 k6

reaction_networks_standard[10] = @reaction_network begin
    p, ∅ ⟶ X1
    (k1,k2), X1 → X2
    (k3,k4), X2 → X3
    (k5,k6), X3 → X4
    (k7,k8), X4 → X5
    d, X5  ⟶ ∅
end p k1 k2 k3 k4 k5 k6 k7 k8 d

### Network with hill functions ###.

reaction_networks_hill[1] = @reaction_network begin
    hillr(X2,v1,K1,n1),   ∅ → X1
    hillr(X1,v2,K2,n2),   ∅ → X2
    (d1,d2),          (X1,X2) → ∅
end v1 v2 K1 K2 n1 n2 d1 d2

reaction_networks_hill[2] = @reaction_network begin
    hillr(X3,v1,K1,n1),   ∅ → X1
    hillr(X1,v2,K2,n2),   ∅ → X2
    hillr(X2,v3,K3,n3),   ∅ → X3
    (d1,d2,d3),          (X1,X2,X3) → ∅
end v1 v2 v3 K1 K2 K3 n1 n2 n3 d1 d2 d3

reaction_networks_hill[3] = @reaction_network begin
    hillr(X2,v1,K1,n1), ∅ → X1
    hill(X1,v2,K2,n2), ∅ → X2
    d, X1+X2 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[4] = @reaction_network begin
    hillr(X2,v1,K1,n1)*hillr(X3,v1,K1,n1), ∅ → X1
    hillr(X1,v2,K2,n2)*hillr(X3,v2,K2,n2), ∅ → X2
    hillr(X1,v3,K3,n3)*hillr(X2,v3,K3,n3), ∅ → X3
    (d1,d2,d3), (X1,X2,X3)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 d1 d2 d3

reaction_networks_hill[5] = @reaction_network begin
    hillr(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
    hill(X5,v2,K2,n2), ∅ → X2
    hill(X3,v3,K3,n3), ∅ → X3
    hillr(X1,v4,K4,n4), ∅ → X4
    hill(X2,v5,K5,n5), ∅ → X5
    (k1,k2), X2 ⟷ X1 + 2X4
    (k3,k4), X4 ⟷ X3
    (k5,k6), 3X5 + X1 ⟷ X2
    (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5

reaction_networks_hill[6] = @reaction_network begin
    v/10+hill(X1,v,K,n), ∅ → X1
    d, X1 → ∅
end v K n d

reaction_networks_hill[7] = @reaction_network begin
    v/10 + hill(X1,v,K,n), ∅ → X1 + X2
    (k1,k2), X1 + X2 ↔ X3
    k3, X3 → X1
    d, (X1,X2,X3) → ∅
end v K n k1 k2 k3 d

reaction_networks_hill[8] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillr(X1,v2,K2,n2)*hill(X3,v3,K3,n3), ∅ → X2
    hill(X2,v4,K4,n4), ∅ → X3
    (d1,d2,d3), (X1,X2,X3) → ∅
end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 d1 d2 d3

reaction_networks_hill[9] = @reaction_network begin
    hill(X1,v1,K1,n1)*hillr(X1,v2,K2,n2), ∅ → X1
    d, X1 → ∅
end v1 K1 n1 v2 K2 n2 d

reaction_networks_hill[10] = @reaction_network begin
    hill(X2,v1,K1,n1), ∅ → X1
    hillr(X4,v2,K2,n2), ∅ → X2
    (k1,k2), 2X1 + X2 ⟷ X3
    (k3,k4), 2X2 + X3 ⟷ X4
    (k5,k6), X1 + X2 + X3 + X4 ⟷ X5 + X6
    (d1,d2), (X5,X6) → ∅
end v1 K1 n1 v2 K2 n2 k1 k2 k3 k4 k5 k6 d1 d2

### Reaction networks were some linnear combination concentrations remain fixed (steady state values depends on initial conditions).

reaction_networks_constraint[1] = @reaction_network begin
    (k1,k2), X1 ↔ X2
    (k3,k4), X2 ↔ X3
    (k5,k6), X3 ↔ X1
end k1 k2 k3 k4 k5 k6
reaction_network_constraints[1] = [1 1 1]

reaction_networks_constraint[2] = @reaction_network begin
    (k1,k2), X1 ↔ 2X1
    (k3,k4), X1 + X2 ↔ X3
    (k5,k6), X3 ↔ X2
end k1 k2 k3 k4 k5 k6
reaction_network_constraints[2] = [0 1 1]

reaction_networks_constraint[3] = @reaction_network begin
    (k1,k2*X5), X1 ↔ X2
    (k3*X5,k4), X3 ↔ X4
    (p+k5*X2*X3,d), ∅ ↔ X5
end k1 k2 k3 k4 p k5 d
reaction_network_constraints[3] = [0 0 1 1 0; 1 1 0 0 0]

reaction_networks_constraint[4] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v,K),d), ∅ ↔ X4
end k1 k2 v K d
reaction_network_constraints[4] = [0 1 1 0; -1 1 0 0]

reaction_networks_constraint[5] = @reaction_network begin
    (k1,k2), X1 ↔ 2X2
    (k3,k4), 2X2 ↔ 3X3
    (k5,k6), 3X3 ↔ 4X4
end k1 k2 k3 k4 k5 k6
reaction_network_constraints[5] = [12 6 4 3]

reaction_networks_constraint[6] = @reaction_network begin
    mmr(X1,v1,K1), X1 → X2
    mmr(X2,v2,K2), X2 → X3
    mmr(X3,v3,K3), X3 → X1
end v1 K1 v2 K2 v3 K3
reaction_network_constraints[6] = [1 1 1]

reaction_networks_constraint[7] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v,K),d), ∅ ↔ X2
    (k3,k4), X2 ↔ X4
end k1 k2 k3 k4 v K d
reaction_network_constraints[7] = [1 0 1 0]

reaction_networks_constraint[8] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (mm(X3,v1,K1),mm(X4,v2,K2)), X3 ↔ X4
end k1 k2 v1 K1 v2 K2
reaction_network_constraints[8] = [-1 1 0 0; 0 1 1 1]

reaction_networks_constraint[9] = @reaction_network begin
    (k1,k2), X1 + X2 ↔ X3
    (k3,k4), X3 + X4 ↔ X5
    (k5,k6), X5 + X6 ↔ X7
end k1 k2 k3 k4 k5 k6
reaction_network_constraints[9] = [ 1 0 1 0 1 0 1; -1 1 0 0 0 0 0; 0 0 0 1 1 0 1; 0 0 0 0 0 1 1 ]

reaction_networks_constraint[10] = @reaction_network begin
    kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
    kDeg,       vPp ⟶ phos
    (kBw,kDw),  2w ⟷ w2
    (kB1,kD1),  w2 + v ⟷ w2v
    (kB2,kD2),  w2v + v ⟷ w2v2
    kK1,        w2v ⟶ w2 + vP
    kK2,        w2v2 ⟶ w2v + vP
    (kB3,kD3),  w2 + σB ⟷ w2σB
    (kB4,kD4),  w2σB + v ⟷ w2v + σB
    (kB5,kD5),  vP + phos ⟷ vPp
    kP,         vPp ⟶ v + phos
    v0*((1+F*σB)/(K+σB)),     ∅ ⟶ σB
    λW*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ w
    λV*v0*((1+F*σB)/(K+σB)),  ∅ ⟶ v
end kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV;
reaction_network_constraints[10] = [ 0 0 0 0 0 0 0 0 1 1 ]

### Reaction networks that are actual models that have been used ###

# Brusselator.
reaction_networks_real[1] = @reaction_network begin
    A, ∅ → X
    1, 2X + Y → 3X
    B, X → Y
    1, X → ∅
end A B;

# The B.subtilis σV Lysozyme stress response.
reaction_networks_real[2] = @reaction_network begin
    v0 + hill(σ,v,K,n), ∅ → (σ+A)
    deg, (σ,A,Aσ) → ∅
    (kB,kD), A + σ ↔ Aσ
    S*kC, Aσ → σ
end v0 v K n kD kB kC deg S;

# A cell cycle model
reaction_networks_real[3] = @reaction_network begin
  k1, 0 --> Y
  k2p, Y --> 0
  k2pp*P, Y --> 0
  (k3p+k3pp*A)/(J3+Po), Po-->P
  (k4*m)/(J4+P), Y + P --> Y + Po
end k1 k2p k2pp k3p k3pp A J3 k4 m J4
#@add_constraint cc_network P+Po=1
#reaction_networks_real[3] = cc_network

# A bistable switch
reaction_networks_real[4] = @reaction_network begin
    d,    (X,Y) → ∅
    hillr(Y,v1,K1,n1), ∅ → X
    hillr(X,v2,K2,n2), ∅ → Y
end d v1 K1 n1 v2 K2 n2

### Reaction networks that contain weird functions, stuff, and other oddities ###

reaction_networks_weird[1] = @reaction_network begin
    exp(p), ∅ → X
    d*exp(X), X → ∅
end p d

reaction_networks_weird[2] = @reaction_network begin
    k1, ∅ → X
    k2*log(12+X), X → Y
    k3*log(3+Y), Y → Z
    log(5,6+k4), Z → ∅
end k1 k2 k3 k4

reaction_networks_weird[3] = @reaction_network begin
    d,    (X,Y) → ∅
    hillr(hill(Y,v2,K2,n2),v1,K1,n1), ∅ → X
    hillr(hill(X,v1,K1,n1),v2,K2,n2), ∅ → Y
end d v1 K1 n1 v2 K2 n2

reaction_networks_weird[4] = @reaction_network begin
    p, ∅ → X1
    (k1,X2^X3), X1 ⟷ X2
    (X2/X4,k4), X2 ⟷ X3
    (k5,X1*X2*X4), X3 ⟷ X4
    d*X2, X4 → ∅
end p k1 k2 k3 k4 k5 k6 d

reaction_networks_weird[5] = @reaction_network begin
    (k1*tanh(X1),k2), X1 ↔ 2X1
    (2+sin(X1),3+sin(k3)), X1 + X2 ↔ X3
    (k5,4+cos(X3+X1)+sin(X2)), X3 ↔ X2
end k1 k2 k3 k4 k5 k6

reaction_networks_weird[6] = @reaction_network begin
    (p,d), ∅ ↔ X1
    (p,d), ∅ ↔ X1
    (p,d), ∅ ↔ X1
    (p,d), ∅ ↔ X1
    (p,d), ∅ ↔ X1
    (p,d), ∅ ↔ X1
end p d

reaction_networks_weird[7] = @reaction_network begin
    k1, X1 → X2
    0, X2 → X3
    k2, X3 → X4
    k3, X4 → X5
end k1 k2 k3

reaction_networks_weird[8] = @reaction_network begin
    k1+X3, X1 → X2
    k2*X4, X2 → X3
    k3/(X5+0.01), X3 → X4
    X6/k4, X4 → X5
    X7^k5, X5 → X6
    k6^X8, X6 → X7
    sqrt(abs(k7*X9)), X7 → X8
    cbrt(abs(k8+X1)), X8 → X9
    X2^3+2X2^2+k9, X9 → X1
end k1 k2 k3 k4 k5 k6 k7 k8 k9

reaction_networks_weird[9] = @reaction_network begin
    p, ∅ → X1
    hill(hill(hill(hill(X1,v1,K1,n1),v2,K2,n2),v3,K3,n3),v4,K4,n4), X1 → ∅
end p v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4

reaction_networks_weird[10] = @reaction_network begin
    d, 5X1 → 4X1
end d

### Gathers all netowkrs in a simgle array ###
reaction_networks_all = [reaction_networks_standard...,
                         reaction_networks_hill...,
                         reaction_networks_constraint...,
                         reaction_networks_real...,
                         reaction_networks_weird...]


text:
# an example on mixed integer programming
#
#   maximize x + 2 y + 5 z
#
#   s.t.  x + y + z <= 10
#         x + 2 y + z <= 15
#
#         x is continuous: 0 <= x <= 5
#         y is integer: 0 <= y <= 10
#         z is binary
#

using Gurobi, Test

@testset "GRB Attributes" begin
    env = Gurobi.Env()

    model = Gurobi.Model(env, "mip_01", :maximize)

    add_cvar!(model, 1., 0., 5.)  # x
    add_ivar!(model, 2., 0, 10)   # y
    add_bvar!(model, 5.)          # z
    update_model!(model)

    add_constr!(model, ones(3), '<', 10.)
    add_constr!(model, [1., 2., 1.], '<', 15.)

    update_model!(model)
    @test Gurobi.get_intattr(model, "NumVars") == 3
    @test Gurobi.get_dblattr(model, "ObjCon") == 0.0
    @test Gurobi.get_strattr(model, "ModelName") == "mip_01"

    @test Gurobi.get_dblattrelement(model, "UB", 1) == 5.0
    @test Gurobi.get_charattrelement(model, "Sense", 1) == '<'

    Gurobi.set_dblattrelement!(model, "UB", 1, 4.0)
    Gurobi.set_charattrelement!(model, "Sense", 1, '>')
    update_model!(model)
    @test Gurobi.get_dblattrelement(model, "UB", 1) == 4.0
    @test Gurobi.get_charattrelement(model, "Sense", 1) == '>'
end

ext: jl


text:
module Pool

using ..Utility: random_id

import Base:
    open, isopen, close,
    write, unsafe_write

export FIFO, delete, path
export pool, shutdown

# ------------------------------------------------------------------------
# Named Pipes
# NOTE: assumes a posix environment!
# TODO: use pycall in the case for windows?
struct FIFO
    root::String
    base::String
    mode::Integer

    function FIFO(base, mode=0o666)
        root = tempdir()
        path = joinpath(root,base)
        if ispath(path)
            return nothing
        end

        err  = ccall(:mkfifo, Cint, (Cstring, Cuint), path, mode)
        systemerror("failed to make fifo at $(repr(path))", err != 0)

        return new(root, base, mode)
    end
end

FIFO() = FIFO(random_id())

path(f::FIFO) = joinpath(f.root, f.base)

function delete(f::FIFO)
    err = ccall(:remove, Cint, (Cstring,), path(f))
    systemerror("failed to delete fifo at $(repr(f.base))", err != 0)
end

struct FIFOStream <: IO
    fd :: Int
end

const RDONLY = 0o0
const WRONLY = 0o1
const RDWR   = 0o2

function open(f::FIFO;
    read     :: Union{Bool,Nothing} = nothing,
    write    :: Union{Bool,Nothing} = nothing,
    create   :: Union{Bool,Nothing} = nothing,
    truncate :: Union{Bool,Nothing} = nothing,
    append   :: Union{Bool,Nothing} = nothing,
)
    # TODO: compute correct open flags
    fd = ccall(:open, Cint, (Cstring, Cint), path(f), WRONLY)
    fd < 0 && error("failed to open FIFO '$(path(f))'")

    return FIFOStream(fd)
end

function open(f::FIFO, mode::AbstractString)::FIFOStream
    mode == "r"  ? Base.open(f, read = true)                  :
    mode == "r+" ? Base.open(f, read = true, write = true)    :
    mode == "w"  ? Base.open(f, truncate = true)              :
    mode == "w+" ? Base.open(f, truncate = true, read = true) :
    mode == "a"  ? Base.open(f, append = true)                :
    mode == "a+" ? Base.open(f, append = true, read = true)   :
    throw(ArgumentError("invalid open mode: $mode"))
end

function open(func::Function, f::FIFO, mode::AbstractString)
    io = open(f, mode)
    try
        func(io)
    finally
        close(io)
    end
end

function close(io::FIFOStream)
    err = ccall(:close, Cint, (Cint,), io.fd)
    err != 0 && error("failed to close FIFOStream $(io.fd)")

    return
end

# write a single byte
function write(io::FIFOStream, x::UInt8)
    n = ccall(:write, Cssize_t, (Cint,Ptr{Cvoid},Csize_t), io.fd, Ref{UInt8}(x), 1)
    n != 1 && error("failed to write to FIFOStream $(io.fd), errno = $(Base.Libc.errno())")
    return 1
end

# write many bytes
function unsafe_write(io::FIFOStream, p::Ptr{UInt8}, n::UInt)
    m = ccall(:write, Cssize_t, (Cint,Ptr{Cvoid},Csize_t), io.fd, p, n)
    m != n && error("failed to write to FIFOStream $(io.fd), errno = $(Base.Libc.errno())")
    return Int(n)
end

# --------------------------------
# worker queues

function pool(size)
    chan = Channel{FIFO}(size)
    for i in 1:size
        f = FIFO()
        while isnothing(f)
            f = FIFO()
        end
        put!(chan, f)
    end

    return chan
end

function shutdown(fifos::Channel{FIFO})
    close(fifos)
    for f in fifos
        delete(f)
    end
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
module PhysicalOptics

using ChainRulesCore
using FFTW
using SpecialFunctions
using QuadGK
using Interpolations
using IndexFunArrays
using FourierTools
using Zygote
using CUDA

to_gpu_or_cpu(arr::AbstractArray, x::AbstractArray) = x
to_gpu_or_cpu(::Type{<:AbstractArray}, x::AbstractArray) = x
is_cuda(arr::AbstractArray) = false
if CUDA.functional()
    is_cuda(arr::CuArray) = true
    to_gpu_or_cpu(b::Bool, x::AbstractArray) = b ? CuArray(x) : x
    to_gpu_or_cpu(x::AbstractArray) = CuArray(x)
    to_gpu_or_cpu(::Type{<:CuArray}, x::AbstractArray) = CuArray(x)
    to_gpu_or_cpu(arr::CuArray, x::AbstractArray) = CuArray{eltype(arr)}(x)
else
    to_gpu_or_cpu(b::Bool, x::AbstractArray) = b ? throw("CUDA.functional() = false") : x
    to_gpu_or_cpu(arr::AbstractArray) = arr
end

export Point

 # define point struct
struct Point{T}
    x::T
    y::T
    z::T
end

const λ0 = 550f-9

 # some useful functions
include("utils.jl")
 # some optics related conversions
include("physical_conversions.jl")

 # simple equations to calculate some resolution criterias
include("resolution_equations.jl")

 # fast FFT based convolutions
include("convolutions.jl")

 # apertures
include("apertures.jl")

 # utils to get useful objects to view
include("utils_view.jl")

 # light propagation
include("propagations.jl")

 # # light propagation based on analytical integrals
include("propagation_integrals.jl")

 # some lens functions
include("lenses.jl")

 # point spread functions
include("psf.jl")

# gaussian beams
include("gaussian.jl")

end # module

max_stars_repo_name: JuliaPhysics/PhysicalOptics

max_issues_repo_name: JuliaPhysics/PhysicalOptics

max_forks_repo_name: JuliaPhysics/PhysicalOptics


text:
# Making sure the Taro package is loaded...
if !isdefined( :Taro )
    using Taro.Workbook
    using Taro.CellStyle
    using Taro.createCellStyle
end  # !isdefined( :Taro )

module ExcelWrapper

# Including inside the module!
try
    using Taro
    Taro.init()  # And properly initialized
catch
end  # try

export Taro, Workbook

if !isdefined( :JavaCall )
    using JavaCall
end  # if !isdefined( :JavaCall )

colnum( c::AbstractString ) = Taro.colnum( c ) + 1

const BorderStyle = @jimport org.apache.poi.ss.usermodel.BorderStyle
const CellRangeAddress = @jimport org.apache.poi.ss.util.CellRangeAddress
const ColIndex = Union{AbstractString, Integer}
const HorizontalAlignment = @jimport org.apache.poi.ss.usermodel.HorizontalAlignment
const ExcelMIMEType = "application/vnd.openxmlformats-officedocument.spreadsheetml.sheet"

const allowedVals = Dict{DataType, Vector{String}}()

# Set up the allowed values for all defined enumeration types.
for enumType in [ HorizontalAlignment, BorderStyle ]
    allowedVals[ enumType ] = Vector{String}()

    for value in jcall( enumType, "values", Array{enumType, 1}, (), )
        push!( allowedVals[ enumType ], jcall( value, "toString", JString, (), ) )
    end  # for value in jcall( enumType, "values", Array{JavaObject, 1}, (), )
end  # for enumType in [ HorizontalAlignment, BorderStyle ]

function CellRangeAddress( firstCol::Integer, firstRow::Integer,
    lastCol::Integer, lastRow::Integer )
    return CellRangeAddress( ( jint, jint, jint, jint, ), firstRow - 1,
        lastRow - 1, firstCol - 1, lastCol - 1 )
end  # CellRangeAddress( firstCol, lastCol, firstRow, lastRow )

function CellRangeAddress( firstCol::AbstractString, firstRow::Integer,
    lastCol::AbstractString, lastRow::Integer )
    return CellRangeAddress( colnum( firstCol ), colnum( lastCol ), firstRow,
        lastRow )
end  # CellRangeAddress( firstCol, lastCol, firstRow, lastRow )

function valueOf( valType::DataType, value::AbstractString )
    if !haskey( allowedVals, valType )
        error( valType, " is not a recognised formatting type." )
    end  # !haskey( allwedVals, valType )

    tmpVal = uppercase( value )

    if !in( tmpVal, allowedVals[ valType ] )
        error( valType, " cannot have the value \"", tmpVal, "\"." )
    end  # !in( tmpVal, allowedVals[ valType ] )

    return jcall( valType, "valueOf", valType, ( JString, ), tmpVal )
end  # valueOf( valType, value )

export getSheet, getRow, getCell, getCellValue, getColumnWidth, getRowHeight,
    getCellStyle, createSheet, createRow, cloneCellStyle, createCell, numRows,
    setCellValue, setColumnWidth, setRowHeight, mergeCells, setCellStyle,
    setHorizontalAlignment, setBorder, setBorderBottom, setBorderTop,
    setBorderLeft, setBorderRight, setBorders

# Redefines the getSheet function such that sheet 1 in the code corresponds to
#   the first sheet in the Excel workbook.
getSheet( w::Workbook, s::Integer ) = Taro.getSheet( w, s - 1 )
getSheet( w::Workbook, s::AbstractString ) = Taro.getSheet( w, s )

# Redefines the getRow function such that row 1 in the code corresponds to row 1
#   in the Excel sheet.
getRow( s::Taro.Sheet, r::Integer ) = Taro.getRow( s, r - 1 )

# Redefines the getCell function such that cell 1 in the code corresponds to
#   column A in the Excel sheet.
getCell( r::Taro.Row, c::Integer ) = Taro.getCell( r, c - 1 )
getCell( r::Taro.Row, c::AbstractString ) = getCell( r, colnum( c ) )

# This function return the cell "rc" from the requested worksheet. If the row r
#   doesn't exist, the function return a NULL cell.
function getCell( s::Taro.Sheet, c::Integer, r::Integer )
    rowObj = getRow( s, r )
    return rowObj.ptr === Ptr{Void}( 0 ) ? Taro.Cell( 0 ) : getCell( rowObj, c )
end  # getCell( s, c, r )

getCell( s::Taro.Sheet, c::AbstractString, r::Integer ) =
    getCell( s, colnum( c ), r )

# Redefines the getCellValue function such that it won't throw an error when
#   called on a NULL cell. Instead, it will return "nothing".
getCellValue( cell::Taro.Cell ) = cell.ptr === Ptr{Void}( 0 ) ? nothing :
    Taro.getCellValue( cell )
getCellValue( s::Taro.Sheet, c::ColIndex, r::Integer ) =
    getCellValue( getCell( s, c, r ) )

# This function allows the user to retrieve the contents of a cell with
#   s[ c, r ].
Base.getindex( s::Taro.Sheet, c::ColIndex, r::Integer ) =
    getCellValue( s, c, r )

# This function retrieves the width of a column of the sheet. As before, column
#   1 corresponds to column "A" in the sheet. The argument "inPixels" determines
#   whether the width is returned in char widths (true) or in pixels.
getColumnWidth( s::Taro.Sheet, c::Integer, inChars::Bool = true ) =
    inChars ? jcall( s, "getColumnWidth", jint, ( jint, ), c - 1 ) / 256 :
        jcall( s, "getColumnWidthInPixels", jfloat, ( jint, ), c - 1 )
getColumnWidth( s::Taro.Sheet, c::AbstractString, inChars::Bool = true ) =
    getColumnWidth( s, colnum( c ), inChars )

# This function retrieves the height of a row of the sheet in points. As before,
#   row 1 corresponds to row 1 of the sheet.
getRowHeight( r::Taro.Row ) = jcall( r, "getHeightInPoints", jfloat, () )

function getRowHeight( s::Taro.Sheet, r::Integer )
    tmpRow = getRow( s,r )
    return tmpRow.ptr === Ptr{Void}( 0 ) ?
        jcall( s, "getDefaultRowHeightInPoints", jfloat, () ) :
        getRowHeight( tmpRow )
end  # getRowHeight( s, r )

# This function gets the style of the given cell. As before, the cell in row 1,
#   column 1 corresponds to cell "A1" of the sheet.
getCellStyle( c::Taro.Cell ) = jcall( c, "getCellStyle", Taro.CellStyle, (), )
getCellStyle( s::Taro.Sheet, c::ColIndex, r::Integer ) =
    getCellStyle( createCell( s, c, r ) )

# This function clones the source cell style to the target style. Important:
#   this function has a bug in Java as it doesn't copy border styles.
function cloneCellStyle( targetStyle::Taro.CellStyle,
    sourceStyle::Taro.CellStyle )

    jcall( targetStyle, "cloneStyleFrom", Void, ( Taro.CellStyle, ),
        sourceStyle )
end  # cloneCellStyle( targetStyle, sourceStyle )

cloneCellStyle( targetStyle::Taro.CellStyle, c::Taro.Cell ) =
    cloneCellStyle( targetStyle, getCellStyle( c ) )
cloneCellStyle( targetStyle::Taro.CellStyle, s::Taro.Sheet, c::ColIndex, r::Integer ) =
    cloneCellStyle( targetStyle, getCellStyle( s, c, r ) )

# Redefines the createSheet function such that it only creates a new sheet IF
#   none of that name exists in the workbook. The function then returns the
#   created/existing sheet.
function createSheet( w::Workbook, s::AbstractString )
    tmpSheet = getSheet( w, s )

    if tmpSheet.ptr === Ptr{Void}( 0 )
        tmpSheet = Taro.createSheet( w, s )
    end  # if tmpSheet.ptr === Ptr{Void}( 0 )

    return tmpSheet
end  # cratesheet( w, s )

# Redefines the createRow function such that row 1 in the code corresponds to
#   row 1 in the Excel sheet. The function only creates the row if it doesn't
#   exist yet in the sheet, and returns the created/existing row.
function createRow( s::Taro.Sheet, r::Integer )
    tmpRow = getRow( s, r )

    if tmpRow.ptr === Ptr{Void}( 0 )
        tmpRow = Taro.createRow( s, r - 1 )
    end  # if tmpRow.ptr === Ptr{Void}( 0 )

    return tmpRow
end  # createRow( s, r )

# Redefines the createCell function such that column 1 in the code corresponds
#   to column A in the Excel sheet.The function only creates the cell if it
#   doesn't exist yet in the row, and returns the created/existing cell.
function createCell( r::Taro.Row, c::Integer )
    tmpCell = getCell( r, c )

    if tmpCell.ptr === Ptr{Void}( 0 )
        tmpCell = Taro.createCell( r, c - 1 )
    end  # if tmpCell.ptr === Ptr{Void}( 0 )

    return tmpCell
end  # createCell( r, c )

createCell( r::Taro.Row, c::AbstractString ) = createCell( r, colnum( c ) )

# This function creates a cell in column c and row r. If the row doesn't exist,
#   the function creates the new row as well and returns the new cell.
createCell( s::Taro.Sheet, c::ColIndex, r::Integer ) =
    createCell( createRow( s, r ), c )

# This function counts the number of defined rows in an Excel sheet. If the
#   second argument is not zero, it checks if the first 'nCols' columns have
#   content and returns the last row that has content in one of those columns.
function numRows( s::Taro.Sheet, nCols::T = 0 ) where T <: Integer
    nRows = jcall( s, "getLastRowNum", jint, () ) + 1

    if nCols <= 0
        return nRows
    end  # if nCols <= 0

    # Check the first nCols of each row to see if anything's actually been
    #   defined in the row, starting from the last.
    isRowEmpty = true

    while ( nRows > 0 ) && isRowEmpty
        isRowEmpty = all( ii -> s[ ii, nRows ] === nothing, 1:nCols )
        nRows -= isRowEmpty ? 1 : 0
    end  # ( nRows > 0 ) && !isRowEmpty

    return nRows
end  # numRows( s, nCols )

function numRows( s::Taro.Sheet, colEnd::String )
    return numRows( s, colnum( colEnd ) )
end  # numRows( s, colEnd )

# This function does the same as the previous numRows function, but checks the
#   columns with a number between "colStart" and "colEnd"
function numRows( s::Taro.Sheet, colStart::T1, colEnd::T2 ) where T1 <: Integer where T2 <: Integer
    nRows = jcall( s, "getLastRowNum", jint, () ) + 1

    if colStart > colEnd
        return nRows
    end  # if colStart > colEnd

    # Check all columns between column "colStart" and "colEnd" to see if
    #   anything's defined in those cells of that row, starting from the last.
    isRowEmpty = true

    while ( nRows > 0 ) && isRowEmpty
        isRowEmpty = all( ii -> s[ ii, nRows ] === nothing, colStart:colEnd )
        nRows -= isRowEmpty ? 1 : 0
    end  # ( nRows > 0 ) && !isRowEmpty

    return nRows
end  # numRows( s, colStart, colEnd )

function numRows( s::Taro.Sheet, colStart::String, colEnd::String )
    return numRows( s, colnum( colStart ), colnum( colEnd ) )
end  # numRows( s, colStart, colEnd )

# Redefines the setCellValue function to bundle the setCellValue and
#   setCellFormula functions. If the string starts with '=', it is considered to
#   be a formula. Otherwise, it's an ordinary string.
setCellValue( c::Taro.Cell, val::AbstractString ) = startswith( val, "=" ) ?
    Taro.setCellFormula( c, val[ 2:end ] ) : Taro.setCellValue( c, val )
setCellValue( c::Taro.Cell, val::Union{Real, Date, DateTime} ) =
    Taro.setCellValue( c, val )
setCellValue( s::Taro.Sheet, c::ColIndex, r::Integer,
    val::Union{AbstractString, Real, Date, DateTime} ) =
    setCellValue( createCell( s, c, r ), val )

# This function allows the user to set a cell's value with s[ c, r ] = val.
Base.setindex!( s::Taro.Sheet, val::Union{AbstractString, Real, Date, DateTime},
    c::ColIndex, r::Integer ) =
    setCellValue( s, c, r, val )

# This function sets the width of a column of the sheet to the given value. As
#   before, column 1 corresponds to column "A" in the sheet. The width of the
#   column is given in char widths.
setColumnWidth( s::Taro.Sheet, c::Integer, w::Real ) =
    jcall( s, "setColumnWidth", Void, ( jint, jint, ), c - 1,
        floor( Integer, w * 256 ) )
setColumnWidth( s::Taro.Sheet, c::AbstractString, w::Real ) =
    setColumnWidth( s, colnum( c ), w )

# This function sets the height of a row of the sheet to the given value. As
#   before, row 1 corresponds to row 1 in the sheet. The height of the row is
#   given in points.
setRowHeight( r::Taro.Row, h::Real ) =
    jcall( r, "setHeightInPoints", Void, ( jfloat, ), h )

setRowHeight( s::Taro.Sheet, r::Integer, h::Real ) =
    setRowHeight( createRow( s, r ), h )

# This function merges the cells in the given range. As always, the cell in
#   row 1, column 1 corresponds to cell "A1".
mergeCells( s::Taro.Sheet, firstCol::Integer, firstRow::Integer,
    lastCol::Integer, lastRow::Integer ) =
    jcall( s, "addMergedRegion", jint, ( CellRangeAddress, ),
        CellRangeAddress( firstCol, firstRow, lastCol, lastRow,  ) )
mergeCells( s::Taro.Sheet, firstCol::AbstractString, firstRow::Integer,
    lastCol::AbstractString, lastRow::Integer ) =
    mergeCells( s, colnum( firstCol ), firstRow, colnum( lastCol ), lastRow )

# This function sets the style of the given cell. As before, the cell in row 1,
#   column 1 corresponds to cell "A1" of the sheet.
setCellStyle( c::Taro.Cell, style::Taro.CellStyle ) =
    jcall( c, "setCellStyle", Void, ( Taro.CellStyle, ), style )
setCellStyle( s::Taro.Sheet, c::ColIndex, r::Integer,
    style::Taro.CellStyle ) =
    setCellStyle( createCell( s, c, r ), style )

# Sets the desired border of the given style.
# XXX Directly passing the border as argument doesn't work since the function
#   jcall cannot properly handle an enum type as argument to pass to the Java
#   function.
function setBorder( style::Taro.CellStyle, border::BorderStyle,
    location::AbstractString )

    tmpLoc = titlecase( location )

    if !in( tmpLoc, [ "Top", "Bottom", "Left", "Right" ] )
        return
    end  # if !in( tmpLoc, [ "top", "bottom", "left", "right" ] )

    borderIndex = jcall( border, "ordinal", jint, (), )
    jcall( style, "setBorder" * tmpLoc, Void, ( jshort, ), borderIndex )
end  # setBorder( style, border, location )

function setBorder( style::Taro.CellStyle, border::AbstractString,
    location::AbstractString )

    tmpBorder = valueOf( BorderStyle, border )
    setBorder( style, tmpBorder, location )
end  # setBorder( style, border, location )

# These functions set the desired border of the given style.
setBorderBottom( style::Taro.CellStyle,
    border::Union{BorderStyle, AbstractString} ) =
    setBorder( style, border, "bottom" )
setBorderTop( style::Taro.CellStyle,
    border::Union{BorderStyle, AbstractString} ) =
    setBorder( style, border, "top" )
setBorderLeft( style::Taro.CellStyle,
    border::Union{BorderStyle, AbstractString} ) =
    setBorder( style, border, "left" )
setBorderRight( style::Taro.CellStyle,
    border::Union{BorderStyle, AbstractString} ) =
    setBorder( style, border, "right" )

# This function sets all the borders of the given style.
function setBorders( style::Taro.CellStyle,
    border::Union{BorderStyle, AbstractString} )

    for loc in [ "bottom", "top", "left", "right" ]
        setBorder( style, border, loc )
    end  # for loc in [ "bottom", "top", "left", "right" ]
end  # setBorders( style, border )

end  # module ExcelWrapper


text:
using Ripserer
using Test

using DataStructures
using Ripserer: DisjointSetsWithBirth, find_leaves!

for arr in (1:10, CartesianIndices([1 2 3 4 5; 1 2 3 4 5]))
    @testset "with $(typeof(collect(arr)))" begin
        s = DisjointSetsWithBirth(arr, reshape(1:10, size(arr)))

        @testset "basic tests" begin
            @test length(s) == 10
            for i in 1:10
                @test find_root!(s, arr[i]) == arr[i]
                @test find_leaves!(s, arr[i]) == [arr[i]]
                @test birth(s, arr[i]) == (i, arr[i])
            end
            @test_throws BoundsError find_root!(s, arr[11])
        end

        @testset "union!, find_root!, find_leaves!" begin
            union!(s, arr[2], arr[3])
            @test find_root!(s, arr[3]) == arr[2]
            @test find_leaves!(s, arr[2]) == [arr[2], arr[3]]

            @test union!(s, arr[8], arr[7]) == arr[8]
            @test union!(s, arr[5], arr[6]) == arr[5]
            @test find_leaves!(s, arr[5]) == [arr[5], arr[6]]
            @test union!(s, arr[8], arr[5]) == arr[8]
            @test find_root!(s, arr[6]) == arr[8]
            @test find_leaves!(s, arr[8]) == [arr[5], arr[6], arr[7], arr[8]]
            union!(s, arr[2], arr[6])
            @test find_root!(s, arr[2]) == arr[8]
            root1 = find_root!(s, arr[6])
            root2 = find_root!(s, arr[2])
            @test root_union!(s, root1, root2) == arr[8]
            @test union!(s, arr[5], arr[6]) == arr[8]
        end
    end
end


text:
module Runtests

using Test, MatrixEquations

@testset "Test MatrixEquations.jl" begin
    include("test_clyap.jl")
    include("test_dlyap.jl")
    include("test_meutil.jl")
    include("test_riccati.jl")
    include("test_sylvester.jl")
    include("test_sylvkr.jl")
    include("test_cplyap.jl")
    include("test_dplyap.jl")
    include("test_mecondest.jl")
end

end


text:
"""
    ApproximateInverseRetraction <: AbstractInverseRetractionMethod

An abstract type for representing approximate inverse retraction methods.
"""
abstract type ApproximateInverseRetraction <: AbstractInverseRetractionMethod end

"""
    NLsolveInverseRetraction{T<:AbstractRetractionMethod,TV,TK} <:
        ApproximateInverseRetraction

An inverse retraction method for approximating the inverse of a retraction using `NLsolve`.

# Constructor

    NLsolveInverseRetraction(
        method::AbstractRetractionMethod[, X0];
        project_tangent=false,
        project_point=false,
        nlsolve_kwargs...,
    )

Constructs an approximate inverse retraction for the retraction `method` with initial guess
`X0`, defaulting to the zero vector. If `project_tangent` is `true`, then the tangent
vector is projected before the retraction using `project`. If `project_point` is `true`,
then the resulting point is projected after the retraction. `nlsolve_kwargs` are keyword
arguments passed to `NLsolve.nlsolve`.
"""
struct NLsolveInverseRetraction{TR<:AbstractRetractionMethod,TV,TK} <:
       ApproximateInverseRetraction
    retraction::TR
    X0::TV
    project_tangent::Bool
    project_point::Bool
    nlsolve_kwargs::TK
    function NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
        isdefined(Manifolds, :NLsolve) ||
            @warn "To use NLsolveInverseRetraction, NLsolve must be loaded using `using NLsolve`."
        return new{typeof(m),typeof(X0),typeof(nlsolve_kwargs)}(
            m,
            X0,
            project_point,
            project_tangent,
            nlsolve_kwargs,
        )
    end
end
function NLsolveInverseRetraction(
    m,
    X0=nothing;
    project_tangent::Bool=false,
    project_point::Bool=false,
    nlsolve_kwargs...,
)
    return NLsolveInverseRetraction(m, X0, project_point, project_tangent, nlsolve_kwargs)
end

@decorator_transparent_signature inverse_retract(
    M::AbstractDecoratorManifold,
    p,
    q,
    m::NLsolveInverseRetraction,
)

@decorator_transparent_signature inverse_retract!(
    M::AbstractDecoratorManifold,
    X,
    p,
    q,
    m::NLsolveInverseRetraction,
)

@doc raw"""
    inverse_retract(M, p, q method::NLsolveInverseRetraction; kwargs...)

Approximate the inverse of the retraction specified by `method.retraction` from `p` with
respect to `q` on the [`Manifold`](@ref) `M` using NLsolve. This inverse retraction is
not guaranteed to succeed and probably will not unless `q` is close to `p` and the initial
guess `X0` is close.

If the solver fails to converge, an [`OutOfInjectivityRadiusError`](@ref) is raised.
See [`NLsolveInverseRetraction`](@ref) for configurable parameters.
"""
inverse_retract(::Manifold, p, q, ::NLsolveInverseRetraction; kwargs...)

function inverse_retract!(M::Manifold, X, p, q, method::NLsolveInverseRetraction; kwargs...)
    X0 = method.X0 === nothing ? zero_tangent_vector(M, p) : method.X0
    res = _inverse_retract_nlsolve(
        M,
        p,
        q,
        method.retraction,
        X0,
        method.project_tangent,
        method.project_point,
        method.nlsolve_kwargs;
        kwargs...,
    )
    if !res.f_converged
        @debug res
        throw(OutOfInjectivityRadiusError())
    end
    return copyto!(X, res.zero)
end


text:
# The following structure is adapted from the OpticSim.jl package @ # d492ca0
import Pkg, Libdl, PackageCompiler

function compile(sysimage_path = "JuliaSysimage.$(Libdl.dlext)"; dev=false)
    env_to_precompile = dirname(dirname(@__DIR__))
    precompile_execution_file = joinpath(@__DIR__, "precompile.jl")
    project_filename = joinpath(env_to_precompile, "Project.toml")
    project = Pkg.API.read_project(project_filename)
    used_packages = Symbol.(collect(keys(project.deps)))

    # Remove unneeded packages
    filter!(x -> x ∉ [:Libdl, :PackageCompiler, :Pkg], used_packages)

    if !dev
        push!(used_packages, :QXContexts)
    end

    @info "Creating QXContexts.jl sysimg: $(sysimage_path)"
    PackageCompiler.create_sysimage(
        used_packages,
        sysimage_path = sysimage_path,
        project = env_to_precompile,
        precompile_execution_file = precompile_execution_file
    )
end


text:
export Heun

"""
    Heun()

The [Heun](https://en.wikipedia.org/wiki/List_of_Runge–Kutta_methods#Heun's_method)
predictor of order 2.
"""
struct Heun <: AbstractPredictor end
struct HeunCache{T} <: AbstractPredictorCache
    dt::Vector{T}
    mk₁::Vector{T}
    mk₂::Vector{T}
end

function cache(::Heun, H, x, ẋ, t)
    HeunCache(promote(dt(H, x, t), copy(ẋ), copy(ẋ))...)
end

function predict!(
    xnext,
    cache::HeunCache,
    H::HomotopyWithCache,
    x,
    t,
    Δt,
    ẋ,
    Jac::JacobianMonitor,
)
    dt, mk₁, mk₂ = cache.dt, cache.mk₁, cache.mk₂
    n = length(xnext)
    @inbounds for i = 1:n
        mk₁[i] = -ẋ[i]
        xnext[i] = x[i] + Δt * ẋ[i]
    end

    minus_ẋ!(mk₂, H, xnext, t + Δt, Jac, dt)

    @inbounds for i = 1:n
        xnext[i] = x[i] - 0.5 * Δt * (mk₁[i] + mk₂[i])
    end

    nothing
end

order(::HeunCache) = 3


text:
module TestMOIWrapper

using Juniper
using MathOptInterface
using Test
import Ipopt

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

function test_SolverName()
    @test MOI.get(Juniper.Optimizer(), MOI.SolverName()) == "Juniper"
    return
end

function test_runtests()
    config = MOI.Test.Config(
        atol = 1e-4,
        rtol = 1e-4,
        optimal_status = MOI.LOCALLY_SOLVED,
        exclude = Any[
            MOI.ConstraintDual,
            MOI.DualObjectiveValue,
            MOI.ConstraintBasisStatus,
            MOI.VariableBasisStatus,
            MOI.NLPBlockDual,
        ],
    )
    optimizer = Juniper.Optimizer(
        Pair{String, Any}[
            "log_levels" => [:Table],
            "nl_solver" => MOI.OptimizerWithAttributes(
                Ipopt.Optimizer,
                "print_level" => 0,
                "sb" => "yes",
            ),
        ],
    )
    caching_optimizer = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        optimizer,
    )
    MOI.Test.runtests(
        caching_optimizer,
        config;
        exclude = String[
            # ====================== Unexplained failures ======================
            # LOCALLY_SOLVED instead of INFEASIBLE
            "test_conic_linear_INFEASIBLE_2",
            # NORM_LIMIT instead of LOCALLY_SOLVED
            "test_conic_linear_VectorOfVariables",
            # INVALID_MODEL instead of LOCALLY_SOLVED
            "test_linear_VectorAffineFunction_empty_row",
            # Attributes are not checked properly on copy_to
            "test_model_copy_to_UnsupportedAttribute",
            "test_model_copy_to_UnsupportedConstraint",
            # supports_constraint not checking for element types
            "test_model_supports_constraint_ScalarAffineFunction_EqualTo",
            "test_model_supports_constraint_VariableIndex_EqualTo",
            "test_model_supports_constraint_VectorOfVariables_Nonnegatives",
            # Result index not being checked properly
            "test_solve_result_index",
            # ======================= Explained failures =======================
            # LOCALLY_INFEASIBLE instead of INFEASIBLE
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
            "test_conic_linear_INFEASIBLE",
            "test_constraint_ZeroOne_bounds_3",
            "test_linear_INFEASIBLE",
            "test_linear_INFEASIBLE_2",
            # NORM_LIMIT instead of DUAL_INFEASIBLE
            "test_solve_TerminationStatus_DUAL_INFEASIBLE",
            "test_linear_DUAL_INFEASIBLE",
            "test_linear_DUAL_INFEASIBLE_2",
        ]
    )
    return
end

end

TestMOIWrapper.runtests()


text:
#Ost - "E6" - восток, Nord - "N8" - север, West - "W4" - запад, Sud - "S2" - юг
function size(r)
    countx = 0
    county = 0
    countx = move_to_wall(r,Ost)
    move_to_wall(r,West)
    county = move_to_wall(r, Nord)
    move_to_wall(r,Sud)
    return countx, county
end

function next_side(side)
    if side == Nord
        return Ost
    end
    if side == Sud
        return West
    end
    if side == West
        return Nord
    end
    if side == Ost
        return Sud
    end
end

function move_to_wall(r, side)
    count = 0
    while isborder(r,side)==false
        move!(r,side)
        count += 1
    end
    return count
end

function move_to_wall_nc(r, side, ncount)
    count = 0
    while isborder(r,side)==false
        move!(r,side)
        count += 1
    end
    if count < ncount
        putmarker!(r)
        return 1
    else
        move_to_wall(r,Sud)
        return 0
    end
end

function rev(side)
    if side == Nord
        return Sud
    end
    if side == Sud
        return Nord
    end
    if side == West
        return Ost
    end
    if side == Ost
        return West
    end
end

function corner(r, side1, side2)
    path = []
    while (isborder(r,side1)&&isborder(r,side2))==false
        if isborder(r,side1)==false
            move!(r,side1)
            push!(path,rev(side1))
        end
        if isborder(r,side2)==false
            move!(r,side2)
            push!(path,rev(side2))
        end
    end
    return path
end

function way_back(r,path)
    for i in 0:(length(path)-1)
        move!(r,path[length(path)-i])
    end
end

function mark_arround(r)
    path = corner(r,Sud,West)
    sizex,sizey = size(r)
    while 1 == 1
        f = move_to_wall_nc(r,Nord,sizey)
        if f == 1
            break
        end
        move!(r,Ost)
    end
    while isborder(r,Nord) == true
        putmarker!(r)
        move!(r,Ost)
        putmarker!(r)
    end
    move!(r,Nord)
    while isborder(r,West) == true
        putmarker!(r)
        move!(r,Nord)
        putmarker!(r)
    end
    move!(r,West)
    while isborder(r,Sud) == true
        putmarker!(r)
        move!(r,West)
        putmarker!(r)
    end
    move!(r,Sud)
    while isborder(r,Ost) == true
        putmarker!(r)
        move!(r,Sud)
        putmarker!(r)
    end
    corner(r,West,Sud)
    way_back(r,path)
end


text:
struct nSunsCAP3_RPE <: AbstractProgramme end
typeProgramme = nSunsCAP3_RPE()
nameProgramme = "nSunsCAP3_RPE"
exerProg = Dict()
### Main lifts
BenchT1 = Exercise(
    name = bench["flat"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 92.5,
)
BenchT2 = Exercise(
    name = bench["cg"],
    equipment = equipment["bb"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 75,
)
push!(exerProg, "BenchT1" => (exercise = BenchT1, progression = CAP3_Bench_T1_RPE))
push!(exerProg, "BenchT2" => (exercise = BenchT2, progression = CAP3_Bench_T2_RPE))

DeadliftT1 = Exercise(
    name = deadlift["dl"],
    equipment = equipment["bb"],
    muscles = [
        muscles["hamstrings"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["quads"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 162.5,
)
DeadliftT2 = Exercise(
    name = deadlift["sumo"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["hamstrings"],
        muscles["spinal erectors"],
        muscles["traps"],
        muscles["lats"],
        muscles["finger flexors"],
    ],
    trainingMax = 140,
)
push!(exerProg, "DeadliftT1" => (exercise = DeadliftT1, progression = CAP3_Deadlift_T1_RPE))
push!(exerProg, "DeadliftT2" => (exercise = DeadliftT2, progression = CAP3_Deadlift_T2_RPE))

SquatT1 = Exercise(
    name = squat["low"],
    equipment = equipment["bb"],
    muscles = [
        muscles["quads"],
        muscles["glutes"],
        muscles["spinal erectors"],
        muscles["transverse abdominis"],
    ],
    trainingMax = 115,
)
SquatT2 = Exercise(
    name = squat["front"],
    equipment = equipment["bb"],
    muscles = [muscles["quads"], muscles["glutes"], muscles["core"]],
    trainingMax = 80,
)
push!(exerProg, "SquatT1" => (exercise = SquatT1, progression = CAP3_Squat_T1_RPE))
push!(exerProg, "SquatT2" => (exercise = SquatT2, progression = CAP3_Squat_T2_RPE))

RowT1 = Exercise(
    name = row["pendlay"],
    modality = modality["cheat"],
    equipment = equipment["bb"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 132.5,
)
RowT2 = Exercise(
    name = row["kroc"],
    equipment = equipment["db"],
    muscles = [
        muscles["lats"],
        muscles["biceps"],
        muscles["rear delts"],
        muscles["rhomboids"],
        muscles["finger flexors"],
    ],
    trainingMax = 42.5,
)
push!(exerProg, "RowT1" => (exercise = RowT1, progression = CAP3_Row_T1_RPE))
push!(exerProg, "RowT2" => (exercise = RowT2, progression = CAP3_Row_T2_RPE))

PressT1 = Exercise(
    name = press["ohp"],
    equipment = equipment["bb"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 60,
)
PressT2 = Exercise(
    name = press["ohp"],
    equipment = equipment["bb"],
    muscles = [muscles["front delts"], muscles["triceps"], muscles["traps"]],
    trainingMax = 60,
)
push!(exerProg, "PressT1" => (exercise = PressT1, progression = DCG_CAP3_OHP_T1_RPE))
push!(exerProg, "PressT2" => (exercise = PressT2, progression = DCG_CAP3_OHP_T2_RPE))

## Accessories
# Squat accessories.
lunges = Exercise(;
    name = squat["lunge"],
    modality = [modality["walking"], modality["backstep"], modality["dragon"]],
    equipment = equipment["db"],
    muscles = [muscles["quads"], muscles["glutes"]],
    trainingMax = 47.5,
)
bulgarianSplit = Exercise(;
    name = squat["bulgarian"],
    equipment = equipment["db"],
    muscles = [muscles["quads"], muscles["glutes"]],
    trainingMax = 37.5,
)
calfRaise = Exercise(;
    name = squat["calf raise"],
    modality = modality["standing"],
    equipment = equipment["db"],
    muscles = muscles["calves"],
    trainingMax = 20,
)
ProgLunges = deepcopy(h23_lmh)
ProgBulg = deepcopy(h23_lmh)
ProgCalfRaise = deepcopy(l23_lmh)
push!(exerProg, "lunges" => (exercise = lunges, progression = ProgLunges))
push!(exerProg, "bulgarianSplit" => (exercise = bulgarianSplit, progression = ProgBulg))
push!(exerProg, "calfRaise" => (exercise = calfRaise, progression = ProgCalfRaise))

# Grip
wristExt = Exercise(;
    name = grip["wrist extension"],
    equipment = equipment["db"],
    muscles = muscles["wrist extensors"],
    trainingMax = 10,
    roundBase = 1,
)
wristCurl = Exercise(;
    name = grip["wrist curl"],
    equipment = equipment["db"],
    muscles = muscles["wrist extensors"],
    trainingMax = 10,
    roundBase = 1,
)
platePinch = Exercise(;
    name = grip["plate pinch"],
    equipment = equipment["plate"],
    muscles = muscles["hand"],
    trainingMax = 27.5,
)
fingerCurl = Exercise(;
    name = grip["finger curl"],
    equipment = equipment["db"],
    muscles = muscles["finger flexors"],
    trainingMax = 20,
)
ProgWristExt = deepcopy(AMRAP_Grip)
ProgWristCurl = deepcopy(AMRAP_Grip)
ProgPlatePinch = deepcopy(AMRAP_Grip)
ProgFingerCurl = deepcopy(AMRAP_Grip)
push!(exerProg, "wristExt" => (exercise = wristExt, progression = ProgWristExt))
push!(exerProg, "wristCurl" => (exercise = wristCurl, progression = ProgWristCurl))
push!(exerProg, "platePinch" => (exercise = platePinch, progression = ProgPlatePinch))
push!(exerProg, "fingerCurl" => (exercise = fingerCurl, progression = ProgFingerCurl))

# Pulling
facePull = Exercise(;
    name = pull["face"],
    equipment = equipment["cable"],
    muscles = [muscles["rear delts"], muscles["rhomboids"], muscles["rot cuff"]],
    trainingMax = 30,
)
pullup = Exercise(;
    name = pull["pullup"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["biceps"], muscles["brachialis"]],
    trainingMax = 12.5,
    roundBase = 1.25,
)
hammerCurl = Exercise(;
    name = pull["hammer"],
    equipment = [equipment["db"], equipment["fat grips"]],
    muscles = [muscles["brachialis"], muscles["biceps"]],
    trainingMax = 20,
)
zottmanCurl = Exercise(;
    name = pull["zottman"],
    equipment = [equipment["db"], equipment["fat grips"]],
    muscles = [muscles["biceps"], muscles["brachialis"]],
    trainingMax = 12.5,
)
ProgFacePull = deepcopy(l23_lmh)
ProgPullup = deepcopy(h23_lmh)
ProgHammerCurl = deepcopy(l23_lmh)
ProgZottmanCurl = deepcopy(l23_lmh)
push!(exerProg, "facePull" => (exercise = facePull, progression = ProgFacePull))
push!(exerProg, "pullup" => (exercise = pullup, progression = ProgPullup))
push!(exerProg, "hammerCurl" => (exercise = hammerCurl, progression = ProgHammerCurl))
push!(exerProg, "zottmanCurl" => (exercise = zottmanCurl, progression = ProgZottmanCurl))

# Pushing
inclineDBBench = Exercise(;
    name = bench["db incline"],
    equipment = equipment["db"],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 45,
)
dips = Exercise(;
    name = bench["dips"],
    equipment = [equipment["parallette"], equipment["chain belt"]],
    muscles = [muscles["pecs"], muscles["triceps"], muscles["front delts"]],
    trainingMax = 12.5,
    roundBase = 1.25,
)
sideRaise = Exercise(;
    name = press["side raise"],
    equipment = equipment["db"],
    muscles = muscles["side delts"],
    trainingMax = 12.5,
)
ProgInclineDBBench = deepcopy(m23_lmh)
ProgDips = deepcopy(m23_lmh)
ProgSideRaise = deepcopy(l23_lmh)
push!(
    exerProg,
    "inclineDBBench" => (exercise = inclineDBBench, progression = ProgInclineDBBench),
)
push!(exerProg, "dips" => (exercise = dips, progression = ProgDips))
push!(exerProg, "sideRaise" => (exercise = sideRaise, progression = ProgSideRaise))

# Core
dragonFlag = Exercise(;
    name = core["dragon flag"],
    equipment = equipment["none"],
    muscles = muscles["core"],
    trainingMax = 10,
)
vsit = Exercise(;
    name = core["v sit"],
    equipment = equipment["parallette"],
    muscles = muscles["core"],
    trainingMax = 10,
)
legRaise = Exercise(;
    name = core["leg raise"],
    equipment = equipment["pullup bar"],
    muscles = muscles["core"],
    trainingMax = 12.5,
)
toeToBar = Exercise(;
    name = core["toes to bar"],
    equipment = equipment["pullup bar"],
    muscles = muscles["core"],
    trainingMax = 10,
)
abRollout = Exercise(;
    name = core["ab roll"],
    modality = modality["kneeling"],
    equipment = equipment["ab wheel"],
    muscles = muscles["core"],
    trainingMax = 10,
)
landMine = Exercise(;
    name = core["landmine"],
    equipment = equipment["bb"],
    muscles = muscles["core"],
    trainingMax = 50,
)
paraCrunch = Exercise(;
    name = core["parallette crunch"],
    equipment = [equipment["parallette"], equipment["db"]],
    muscles = muscles["core"],
    trainingMax = 15,
)
ProgDragonFlag = deepcopy(AMRAP_Cali)
ProgVsit = deepcopy(AMRAP_Cali)
ProgLegRaise = deepcopy(m23_lmh)
ProgToeToBar = deepcopy(AMRAP_Cali)
ProgAbRollout = deepcopy(AMRAP_Cali)
ProgLandMine = deepcopy(m23_lmh)
ProgParaCrunch = deepcopy(m23_lmh)
push!(exerProg, "dragonFlag" => (exercise = dragonFlag, progression = ProgDragonFlag))
push!(exerProg, "vsit" => (exercise = vsit, progression = ProgVsit))
push!(exerProg, "legRaise" => (exercise = legRaise, progression = ProgLegRaise))
push!(exerProg, "toeToBar" => (exercise = toeToBar, progression = ProgToeToBar))
push!(exerProg, "abRollout" => (exercise = abRollout, progression = ProgAbRollout))
push!(exerProg, "landMine" => (exercise = landMine, progression = ProgLandMine))
push!(exerProg, "paraCrunch" => (exercise = paraCrunch, progression = ProgParaCrunch))

# Calisthenics Isometrics
frontLever = Exercise(;
    name = isometric["front"],
    modality = modality["one leg out"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["core"]],
    trainingMax = 15,
)
backLever = Exercise(;
    name = isometric["back"],
    modality = modality["straddle"],
    equipment = equipment["pullup bar"],
    muscles = [muscles["lats"], muscles["pecs"], muscles["front delts"], muscles["core"]],
    trainingMax = 15,
)
ProgFrontLever = deepcopy(AMRAP_Cali)
ProgBackLever = deepcopy(AMRAP_Cali)
push!(exerProg, "frontLever" => (exercise = frontLever, progression = ProgFrontLever))
push!(exerProg, "backLever" => (exercise = backLever, progression = ProgBackLever))

## Weeks
import Lifting: makeDays
function makeDays(::nSunsCAP3_RPE, dict::Dict)
    week1 = [[], [], [], [], [], [], []]
    week2 = [[], [], [], [], [], [], []]
    week3 = [[], [], [], [], [], [], []]

    push!(week1[1], dict["BenchT1"][1], dict["BenchT1"][2], 1)
    push!(week1[1], dict["BenchT2"][1], dict["BenchT2"][2], 1)
    push!(week1[1], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week1[1], dict["wristExt"][1], dict["wristExt"][2])
    push!(week1[1], dict["dragonFlag"][1], dict["dragonFlag"][2])
    push!(week1[1], dict["frontLever"][1], dict["frontLever"][2])
    push!(week1[1], dict["lunges"][1], dict["lunges"][2], 2)
    push!(week1[1], dict["facePull"][1], dict["facePull"][2], 2)
    push!(week1[1], dict["legRaise"][1], dict["legRaise"][2], 2)

    push!(week1[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 1)
    push!(week1[2], dict["RowT2"][1], dict["RowT2"][2], 1)
    push!(week1[2], dict["sideRaise"][1], dict["sideRaise"][2], 2)
    push!(week1[2], dict["hammerCurl"][1], dict["hammerCurl"][2], 2)
    push!(week1[2], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 3)
    push!(week1[2], dict["platePinch"][1], dict["platePinch"][2])
    push!(week1[2], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 2)
    push!(week1[2], dict["vsit"][1], dict["vsit"][2])

    push!(week1[3], dict["SquatT1"][1], dict["SquatT1"][2], 1)
    push!(week1[3], dict["PressT2"][1], dict["PressT2"][2], 1)
    push!(week1[3], dict["backLever"][1], dict["backLever"][2])
    push!(week1[3], dict["pullup"][1], dict["pullup"][2], 2)
    push!(week1[3], dict["calfRaise"][1], dict["calfRaise"][2], 1)
    push!(week1[3], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week1[3], dict["dips"][1], dict["dips"][2], 3)
    push!(week1[3], dict["abRollout"][1], dict["abRollout"][2])

    push!(week1[4], dict["BenchT1"][1], dict["BenchT1"][2], 2)
    push!(week1[4], dict["BenchT2"][1], dict["BenchT2"][2], 2)
    push!(week1[4], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week1[4], dict["wristExt"][1], dict["wristExt"][2])
    push!(week1[4], dict["landMine"][1], dict["landMine"][2], 1)
    push!(week1[4], dict["frontLever"][1], dict["frontLever"][2])
    push!(week1[4], dict["bulgarianSplit"][1], dict["bulgarianSplit"][2], 3)
    push!(week1[4], dict["facePull"][1], dict["facePull"][2], 2)
    push!(week1[4], dict["toeToBar"][1], dict["toeToBar"][2])

    push!(week1[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 1)
    push!(week1[5], dict["RowT1"][1], dict["RowT1"][2], 1)
    push!(week1[5], dict["sideRaise"][1], dict["sideRaise"][2], 2)
    push!(week1[5], dict["hammerCurl"][1], dict["hammerCurl"][2], 2)
    push!(week1[5], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 3)
    push!(week1[5], dict["platePinch"][1], dict["platePinch"][2])
    push!(week1[5], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 2)
    push!(week1[5], dict["vsit"][1], dict["vsit"][2])

    push!(week1[6], dict["SquatT2"][1], dict["SquatT2"][2], 1)
    push!(week1[6], dict["PressT1"][1], dict["PressT1"][2], 1)
    push!(week1[6], dict["backLever"][1], dict["backLever"][2])
    push!(week1[6], dict["pullup"][1], dict["pullup"][2], 2)
    push!(week1[6], dict["calfRaise"][1], dict["calfRaise"][2], 1)
    push!(week1[6], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week1[6], dict["dips"][1], dict["dips"][2], 3)
    push!(week1[6], dict["paraCrunch"][1], dict["paraCrunch"][2], 1)

    # Week 2
    push!(week2[1], dict["BenchT1"][1], dict["BenchT1"][2], 3)
    push!(week2[1], dict["BenchT2"][1], dict["BenchT2"][2], 3)
    push!(week2[1], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week2[1], dict["wristExt"][1], dict["wristExt"][2])
    push!(week2[1], dict["dragonFlag"][1], dict["dragonFlag"][2])
    push!(week2[1], dict["frontLever"][1], dict["frontLever"][2])
    push!(week2[1], dict["lunges"][1], dict["lunges"][2], 1)
    push!(week2[1], dict["facePull"][1], dict["facePull"][2], 3)
    push!(week2[1], dict["legRaise"][1], dict["legRaise"][2], 3)

    push!(week2[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 2)
    push!(week2[2], dict["RowT2"][1], dict["RowT2"][2], 2)
    push!(week2[2], dict["sideRaise"][1], dict["sideRaise"][2], 3)
    push!(week2[2], dict["hammerCurl"][1], dict["hammerCurl"][2], 3)
    push!(week2[2], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 1)
    push!(week2[2], dict["platePinch"][1], dict["platePinch"][2])
    push!(week2[2], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 3)
    push!(week2[2], dict["vsit"][1], dict["vsit"][2])

    push!(week2[3], dict["SquatT1"][1], dict["SquatT1"][2], 2)
    push!(week2[3], dict["PressT2"][1], dict["PressT2"][2], 2)
    push!(week2[3], dict["backLever"][1], dict["backLever"][2])
    push!(week2[3], dict["pullup"][1], dict["pullup"][2], 3)
    push!(week2[3], dict["calfRaise"][1], dict["calfRaise"][2], 2)
    push!(week2[3], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week2[3], dict["dips"][1], dict["dips"][2], 1)
    push!(week2[3], dict["abRollout"][1], dict["abRollout"][2])

    push!(week2[4], dict["BenchT1"][1], dict["BenchT1"][2], 4)
    push!(week2[4], dict["BenchT2"][1], dict["BenchT2"][2], 4)
    push!(week2[4], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week2[4], dict["wristExt"][1], dict["wristExt"][2])
    push!(week2[4], dict["landMine"][1], dict["landMine"][2], 3)
    push!(week2[4], dict["frontLever"][1], dict["frontLever"][2])
    push!(week2[4], dict["bulgarianSplit"][1], dict["bulgarianSplit"][2], 1)
    push!(week2[4], dict["facePull"][1], dict["facePull"][2], 3)
    push!(week2[4], dict["toeToBar"][1], dict["toeToBar"][2])

    push!(week2[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 2)
    push!(week2[5], dict["RowT1"][1], dict["RowT1"][2], 2)
    push!(week2[5], dict["sideRaise"][1], dict["sideRaise"][2], 3)
    push!(week2[5], dict["hammerCurl"][1], dict["hammerCurl"][2], 3)
    push!(week2[5], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 1)
    push!(week2[5], dict["platePinch"][1], dict["platePinch"][2])
    push!(week2[5], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 3)
    push!(week2[5], dict["vsit"][1], dict["vsit"][2])

    push!(week2[6], dict["SquatT2"][1], dict["SquatT2"][2], 2)
    push!(week2[6], dict["PressT1"][1], dict["PressT1"][2], 2)
    push!(week2[6], dict["backLever"][1], dict["backLever"][2])
    push!(week2[6], dict["pullup"][1], dict["pullup"][2], 3)
    push!(week2[6], dict["calfRaise"][1], dict["calfRaise"][2], 2)
    push!(week2[6], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week2[6], dict["dips"][1], dict["dips"][2], 1)
    push!(week2[6], dict["paraCrunch"][1], dict["paraCrunch"][2], 3)

    # Week3
    push!(week3[1], dict["BenchT1"][1], dict["BenchT1"][2], 5)
    push!(week3[1], dict["BenchT2"][1], dict["BenchT2"][2], 5)
    push!(week3[1], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week3[1], dict["wristExt"][1], dict["wristExt"][2])
    push!(week3[1], dict["dragonFlag"][1], dict["dragonFlag"][2])
    push!(week3[1], dict["frontLever"][1], dict["frontLever"][2])
    push!(week3[1], dict["lunges"][1], dict["lunges"][2], 3)
    push!(week3[1], dict["facePull"][1], dict["facePull"][2], 1)
    push!(week3[1], dict["legRaise"][1], dict["legRaise"][2], 1)

    push!(week3[2], dict["DeadliftT2"][1], dict["DeadliftT2"][2], 3)
    push!(week3[2], dict["RowT2"][1], dict["RowT2"][2], 3)
    push!(week3[2], dict["sideRaise"][1], dict["sideRaise"][2], 1)
    push!(week3[2], dict["hammerCurl"][1], dict["hammerCurl"][2], 1)
    push!(week3[2], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 2)
    push!(week3[2], dict["platePinch"][1], dict["platePinch"][2])
    push!(week3[2], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 1)
    push!(week3[2], dict["vsit"][1], dict["vsit"][2])

    push!(week3[3], dict["SquatT1"][1], dict["SquatT1"][2], 3)
    push!(week3[3], dict["PressT2"][1], dict["PressT2"][2], 3)
    push!(week3[3], dict["backLever"][1], dict["backLever"][2])
    push!(week3[3], dict["pullup"][1], dict["pullup"][2], 1)
    push!(week3[3], dict["calfRaise"][1], dict["calfRaise"][2], 3)
    push!(week3[3], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week3[3], dict["dips"][1], dict["dips"][2], 2)
    push!(week3[3], dict["abRollout"][1], dict["abRollout"][2])

    push!(week3[4], dict["BenchT1"][1], dict["BenchT1"][2], 6)
    push!(week3[4], dict["BenchT2"][1], dict["BenchT2"][2], 6)
    push!(week3[4], dict["wristCurl"][1], dict["wristCurl"][2])
    push!(week3[4], dict["wristExt"][1], dict["wristExt"][2])
    push!(week3[4], dict["landMine"][1], dict["landMine"][2], 2)
    push!(week3[4], dict["frontLever"][1], dict["frontLever"][2])
    push!(week3[4], dict["bulgarianSplit"][1], dict["bulgarianSplit"][2], 2)
    push!(week3[4], dict["facePull"][1], dict["facePull"][2], 1)
    push!(week3[4], dict["toeToBar"][1], dict["toeToBar"][2])

    push!(week3[5], dict["DeadliftT1"][1], dict["DeadliftT1"][2], 3)
    push!(week3[5], dict["RowT1"][1], dict["RowT1"][2], 3)
    push!(week3[5], dict["sideRaise"][1], dict["sideRaise"][2], 1)
    push!(week3[5], dict["hammerCurl"][1], dict["hammerCurl"][2], 1)
    push!(week3[5], dict["inclineDBBench"][1], dict["inclineDBBench"][2], 2)
    push!(week3[5], dict["platePinch"][1], dict["platePinch"][2])
    push!(week3[5], dict["zottmanCurl"][1], dict["zottmanCurl"][2], 1)
    push!(week3[5], dict["vsit"][1], dict["vsit"][2])

    push!(week3[6], dict["SquatT2"][1], dict["SquatT2"][2], 3)
    push!(week3[6], dict["PressT1"][1], dict["PressT1"][2], 3)
    push!(week3[6], dict["backLever"][1], dict["backLever"][2])
    push!(week3[6], dict["pullup"][1], dict["pullup"][2], 1)
    push!(week3[6], dict["calfRaise"][1], dict["calfRaise"][2], 3)
    push!(week3[6], dict["fingerCurl"][1], dict["fingerCurl"][2])
    push!(week3[6], dict["dips"][1], dict["dips"][2], 2)
    push!(week3[6], dict["paraCrunch"][1], dict["paraCrunch"][2], 2)

    return [
        week1[1],
        week1[2],
        week1[3],
        week1[4],
        "Rest",
        week1[5],
        week1[6],
        week2[1],
        week2[2],
        "Rest",
        week2[3],
        week2[4],
        week2[5],
        week2[6],
        "Rest",
        week3[1],
        week3[2],
        week3[3],
        week3[4],
        "Rest",
        week3[5],
        week3[6],
        "Rest",
    ]
end

days = makeDays(nSunsCAP3_RPE(), exerProg)
push!(
    Lifting_Programmes,
    nameProgramme => Programme(typeProgramme, nameProgramme, exerProg, days),
)

max_stars_repo_name: dcelisgarza/Lifting


text:
push!(LOAD_PATH,"simulator")
using EVQueues, Plots

lambda=120.0;
mu=1.0;
gamma=0.5;

C=60.0;

Tfinal=1000.0;

sim = ev_llf(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
compute_statistics!(sim)

#sim2 = ev_edffixed(lambda,mu,gamma,Tfinal,C,snapshots=[Tfinal])
#compute_statistics!(sim2)


text:
#!/usr/local/bin/julia
# last tested Julia version: 1.6.1
# Hodgkin-Huxley model on a 2D lattice
# FvW 03/2018

using NPZ
using PyCall
using PyPlot
using Statistics
using VideoIO
@pyimport matplotlib.animation as anim

function hh2d(N, T, t0, dt, s, D, C, gNa, gK, gL, VNa, VK, VL, I0, stim, blocks)
    # initialize Hodgkin-Huxley system
	C1 = 1.0/C
    v = -10*ones(Float64,N,N)
    m = zeros(Float64,N,N)
    n = zeros(Float64,N,N)
    h = zeros(Float64,N,N)
    dvdt = zeros(Float64,N,N)
    dmdt = zeros(Float64,N,N)
    dndt = zeros(Float64,N,N)
    dhdt = zeros(Float64,N,N)
    s_sqrt_dt = s*sqrt(dt)
    X = zeros(Float64,T,N,N)
    # stimulation protocol
    I = zeros(Float64,t0+T,N,N)
    for st in stim
        t_on, t_off = st[1]
        x0, x1 = st[2]
        y0, y1 = st[3]
        I[t0+t_on:t0+t_off, x0:x1, y0:y1] .= I0
    end
    # iterate
    for t in range(2, stop=t0+T, step=1)
        (t%100 == 0) && print("    t = ", t, "/", t0+T, "\r")
        # HH equations
        # dV/dt (membrane potential)
        dvdt = C1.*(-gNa.*(m.^3).*h.*(v.-VNa)
                    .- gK.*(n.^4).*(v.-VK)
                    .- gL*(v.-VL) .+ I[t,:,:]) .+ D.*L(v)
        # dm/dt (Na+ channel activation)
        alpha_m = 0.1.*(-v.+25)./(exp.((-v.+25)./10).-1)
        beta_m = 4 * exp.(-v./18)
        dmdt = (alpha_m .* (1 .- m) .- beta_m .* m)
        # dh/dt (Na+ channel inactivation)
        alpha_h = 0.07 * exp.(-v./20)
        beta_h = 1.0 ./ (exp.((-v .+ 30) ./ 10) .+ 1)
        dhdt = (alpha_h .* (1.0 .- h) .- beta_h.*h)
        # dn/dt (K+ channel gating)
        alpha_n = 0.01 .* (-v .+ 10)./(exp.((-v .+ 10)./10) .- 1)
        beta_n = 0.125 .* exp.(-v./80)
        dndt = (alpha_n .* (1.0 .- n) .- beta_n.*n)
        # stochastic integration
        v += (dvdt*dt + s_sqrt_dt*randn(N,N))
        m += (dmdt*dt)
        n += (dndt*dt)
        h += (dhdt*dt)
        # dead block(s):
        for bl in blocks
            v[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
            m[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
            n[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
            h[bl[1][1]:bl[1][2], bl[2][1]:bl[2][2]] .= 0.0
        end
        # store
        (t > t0) && (X[t-t0,:,:] = v)
    end
    println("\n")
    return X
end

function animate_pyplot(fname, data)
    """
    Animate 3D array as .mp4 using PyPlot, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    nt, nx, ny = size(data)
    vmin = minimum(data)
    vmax = maximum(data)
    # setup animation image
    fig = figure(figsize=(6,6))
    axis("off")
    t = imshow(data[1,:,:], origin="lower", cmap=ColorMap("gray"),
			   vmin=vmin, vmax=vmax)
    tight_layout()
    # frame generator
    println("[+] animate")
    function animate(i)
        (i%100 == 0) && print("    t = ", i, "/", nt, "\r")
        t.set_data(data[i+1,:,:])
    end
    # create animation
    ani = anim.FuncAnimation(fig, animate, frames=nt, interval=10)
    println("\n")
    # save animation
    ani[:save](fname, bitrate=-1,
               extra_args=["-vcodec", "libx264", "-pix_fmt", "yuv420p"])
    show()
end

function animate_video(fname, data)
    """
    Animate 3D array as .mp4 using VideoIO, save as `fname`
    array dimensions:
        1: time
        2, 3: space
    """
    # BW
    y = UInt8.(round.(255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    # BW inverted
    #y = UInt8.(round.(255 .- 255*(data .- minimum(data))/(maximum(data)-minimum(data))))
    encoder_options = (color_range=2, crf=0, preset="medium")
    framerate=30
    T = size(data,1)
    open_video_out(fname, y[1,end:-1:1,:], framerate=framerate,
                   encoder_options=encoder_options) do writer
        for i in range(2,stop=T,step=1)
            write(writer, y[i,end:-1:1,:])
        end
    end
end

function L(x)
    # Laplace operator
    # periodic boundary conditions
    xU = circshift(x, [-1 0])
    xD = circshift(x, [1 0])
    xL = circshift(x, [0 -1])
    xR = circshift(x, [0 1])
    Lx = xU + xD + xL + xR - 4x
    # non-periodic boundary conditions
    Lx[1,:] .= 0.0
    Lx[end,:] .= 0.0
    Lx[:,1] .= 0.0
    Lx[:,end] .= 0.0
    return Lx
end

function main()
    println("Hodgkin-Huxley (HH) lattice model\n")
    N = 128
    T = 5000
    t0 = 2000
    dt = 0.01
    s = 1.0
    D = 1.0
    # Hodgkin-Huxley parameters
    C = 1
    gNa = 120
    gK = 36
    gL = 0.3
    VNa = 115
    VK = -12
    VL = 10.6
    I = 45.0
    println("[+] Lattice size N: ", N)
    println("[+] Time steps T: ", T)
    println("[+] Warm-up steps t0: ", t0)
    println("[+] Integration time step dt: ", dt)
    println("[+] Noise intensity s: ", s)
    println("[+] Diffusion coefficient D: ", D)
    println("[+] HH parameter gNa: ", gNa)
    println("[+] HH parameter gK: ", gK)
    println("[+] HH parameter gL: ", gL)
    println("[+] HH parameter VNa: ", VNa)
    println("[+] HH parameter VK: ", VK)
    println("[+] HH parameter VL: ", VL)
    println("[+] HH parameter C: ", C)
    println("[+] Stimulation current I: ", I)

    # stimulation protocol, array of elements [[t0,t1], [x0,x1], [y0,y1]]
    stim = [ [[50,350], [1,5], [1,5]],
             [[1900,2200], [25,30], [1,15]] ]
    #stim = [] # empty, no current

    # dead blocks, array of elementy [[x0,x1], [y0,y1]]
    blocks = [ [[1,15], [5,10]] ]
    #blocks = [] # empty, no dead areas

    # run simulation
    data = hh2d(N, T, t0, dt, s, D, C, gNa, gK, gL, VNa, VK, VL, I, stim, blocks)
    println("[+] Data dimensions: ", size(data))

    # plot mean voltage
    m = mean(reshape(data, (T,N*N)), dims=2)
    plot(m, "-k"); show()

    # save data
    I_str = rpad(I, 4, '0') # stim. current amplitude as 4-char string
    s_str = rpad(s, 4, '0') # noise as 4-char string
    D_str = rpad(D, 4, '0') # diffusion coefficient as 4-char string
    #fname1 = string("ml2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".npy")
    #npzwrite(fname1, data)
    #println("[+] Data saved as: ", data_filename)

    # video
    fname2 = string("hh2d_I_", I_str, "_s_", s_str, "_D_", D_str, ".mp4")
    #animate_pyplot(fname2, data) # slow
    animate_video(fname2, data) # fast
    println("[+] Animation saved as: ", fname2)
end

main()


text:
"""
     struct Table{T,Vd<:AbstractVector{T},Vp<:AbstractVector} <: AbstractVector{Vector{T}}
        data::Vd
        ptrs::Vp
     end

Type representing a list of lists (i.e., a table) in
compressed format.
"""
struct Table{T,Vd<:AbstractVector{T},Vp<:AbstractVector} <: AbstractVector{Vector{T}}
  data::Vd
  ptrs::Vp
  function Table(data::AbstractVector,ptrs::AbstractVector)
    new{eltype(data),typeof(data),typeof(ptrs)}(data,ptrs)
  end
end

"""
    Table(a::AbstractArray{<:AbstractArray})

Build a table from a vector of vectors. If the inputs are
multidimensional arrays instead of vectors, they are flattened.
"""
function Table(a::AbstractArray{<:AbstractArray})
  data, ptrs = generate_data_and_ptrs(a)
  Table(data,ptrs)
end

function Table(a::Table)
  a
end

function Base.convert(::Type{Table{T,Vd,Vp}},table::Table{Ta,Vda,Vpa}) where {T,Vd,Vp,Ta,Vda,Vpa}
  data = convert(Vd,table.data)
  ptrs = convert(Vp,table.ptrs)
  Table(data,ptrs)
end

function Base.convert(::Type{Table{T,Vd,Vp}},table::Table{T,Vd,Vp}) where {T,Vd,Vp}
  table
end

"""
"""
function identity_table(::Type{T},::Type{P},l::Integer) where {T,P}
  data = Vector{T}(1:l)
  ptrs = Vector{P}(1:l+1)
  Table(data,ptrs)
end

"""
    empty_table(::Type{T},::Type{P}, l::Integer) where {T,P}
    empty_table(l::Integer)
"""
function empty_table(::Type{T},::Type{P}, l::Integer) where {T,P}
  data = T[]
  ptrs = ones(P,l+1)
  Table(data,ptrs)
end

empty_table(l::Integer) = empty_table(Int,Int32,l)

size(a::Table) = (length(a.ptrs)-1,)

IndexStyle(::Type{<:Table}) = IndexLinear()

function array_cache(a::Table)
  if length(a.ptrs) > 1
    pini = a.ptrs[1]
    pend = a.ptrs[2]
  else
    pend = a.ptrs[1]
    pini = pend
  end
  T = eltype(a.data)
  l = pend - pini
  r = zeros(T,l)
  CachedArray(r)
end

function getindex!(c,a::Table,i::Integer)
  pini = a.ptrs[i]
  l = a.ptrs[i+1] - pini
  setsize!(c,(l,))
  pini -= 1
  r = c.array
  for j in 1:l
     @inbounds r[j] = a.data[pini+j]
  end
  r
end

function Base.getindex(a::Table,i::Integer)
  cache = array_cache(a)
  getindex!(cache,a,i)
end

function Base.getindex(a::Table,i::UnitRange)
  r = a.ptrs[i.start]:(a.ptrs[i.stop+1]-1)
  data = a.data[r]
  r = i.start:(i.stop+1)
  ptrs = a.ptrs[r]
  o = ptrs[1]-1
  ptrs .-= o
  Table(data,ptrs)
end

function Base.getindex(a::Table,ids::AbstractVector{<:Integer})
  ptrs = similar(a.ptrs,eltype(a.ptrs),length(ids)+1)
  for (i,id) in enumerate(ids)
    ptrs[i+1] = a.ptrs[id+1]-a.ptrs[id]
  end
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = similar(a.data,eltype(a.data),ndata)
  for (i,id) in enumerate(ids)
    n = a.ptrs[id+1]-a.ptrs[id]
    p1 = ptrs[i]-1
    p2 = a.ptrs[id]-1
    for j in 1:n
      data[p1+j] = a.data[p2+j]
    end
  end
  Table(data,ptrs)
end

# Helper functions related with Tables

"""
    rewind_ptrs!(ptrs)

Rewind the given vector of pointers.
"""
function rewind_ptrs!(ptrs::AbstractVector{<:Integer})
  @inbounds for i in (length(ptrs)-1):-1:1
    ptrs[i+1] = ptrs[i]
  end
  ptrs[1] = 1
end

"""
    length_to_ptrs!(ptrs)

Given a vector of integers, mutate it from length state to pointer state.
"""
function length_to_ptrs!(ptrs::AbstractArray{<:Integer})
  ptrs[1] = 1
  @inbounds for i in 1:(length(ptrs)-1)
    ptrs[i+1] += ptrs[i]
  end
end

"""
    data, ptrs = generate_data_and_ptrs(vv)

Given a vector of vectors, compress it and return the corresponding data and and ptrs
"""
function generate_data_and_ptrs(vv::AbstractArray{<:AbstractArray{T}}) where T
  ptrs = Vector{Int32}(undef,length(vv)+1)
  _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  length_to_ptrs!(ptrs)
  ndata = ptrs[end]-1
  data = Vector{T}(undef,ndata)
  _generate_data_and_ptrs_fill_data!(data,vv)
  (data, ptrs)
end

function _generate_data_and_ptrs_fill_ptrs!(ptrs,vv)
  c = array_cache(vv)
  k = 1
  for i in eachindex(vv)
    v = getindex!(c,vv,i)
    ptrs[k+1] = length(v)
    k += 1
  end
end

function _generate_data_and_ptrs_fill_data!(data,vv)
  c = array_cache(vv)
  k = 1
  for i in eachindex(vv)
    v = getindex!(c,vv,i)
    for vi in v
      data[k] = vi
      k += 1
    end
  end
end

"""
    append_ptrs(pa,pb)

Append two vectors of pointers.
"""
function append_ptrs(pa::AbstractVector{T},pb::AbstractVector{T}) where T
  p = copy(pa)
  append_ptrs!(p,pb)
end

"""
"""
function append_ptrs!(pa::AbstractVector{T},pb::AbstractVector{T}) where T
  na = length(pa)-1
  nb = length(pb)-1
  _append_grow!(pa,nb,zero(T))
  _append_count!(pa,pb,na,nb)
  rewind_ptrs!(pa)
  length_to_ptrs!(pa)
  pa
end

function _append_grow!(pa,nb,z)
  for i in 1:nb
    push!(pa,z)
  end
end

function _append_count!(pa,pb,na,nb)
  for ca in 1:na
    pa[ca] = pa[ca+1] - pa[ca]
  end
  for cb in 1:nb
    pa[cb+na] = pb[cb+1] - pb[cb]
  end
end

"""
"""
const UNSET = 0

"""
"""
function find_inverse_index_map(a_to_b, nb=maximum(a_to_b))
  T = eltype(a_to_b)
  b_to_a = fill(T(UNSET),nb)
  find_inverse_index_map!(b_to_a, a_to_b)
  b_to_a
end

"""
"""
function find_inverse_index_map!(b_to_a, a_to_b)
  for (a,b) in enumerate(a_to_b)
    if b != UNSET
      @inbounds b_to_a[b] = a
    end
  end
end

"""
"""
function append_tables_globally(tables::Table{T,Vd,Vp}...) where {T,Vd,Vp}
  first_table, = tables
  data = copy(first_table.data)
  ptrs = copy(first_table.ptrs)
  for (i,table) in enumerate(tables)
    if  i != 1
      append!(data,table.data)
      append_ptrs!(ptrs,table.ptrs)
    end
  end
  Table(data,ptrs)
end

function append_tables_globally()
  @unreachable "At least one table has to be provided"
end

"""
"""
get_ptrs_eltype(::Table{T,Vd,Vp}) where {T,Vd,Vp} = eltype(Vp)
get_ptrs_eltype(::Type{Table{T,Vd,Vp}}) where {T,Vd,Vp} = eltype(Vp)

"""
"""
get_data_eltype(::Table{T,Vd,Vp}) where {T,Vd,Vp} = T
get_data_eltype(::Type{Table{T,Vd,Vp}}) where {T,Vd,Vp} = T

"""
    append_tables_locally(tables::Table...)
"""
function append_tables_locally(tables::Table...)
  n = length(tables)
  offsets = tfill(0,Val{n}())
  append_tables_locally(offsets,tables)
end

"""
"""
function append_tables_locally(offsets::NTuple, tables::NTuple)

  first_table, = tables

  @assert all( map(length,tables) .== length(first_table) ) "All tables must have the same length"
  ndata = sum( (length(table.data) for table in tables) )

  T = get_data_eltype(first_table)
  P = get_ptrs_eltype(first_table)

  ptrs = zeros(P,length(first_table.ptrs))
  data = zeros(T,ndata)

  for table in tables
    _append_tables_locally_count!(ptrs,table)
  end

  length_to_ptrs!(ptrs)

  for (offset,table) in zip(offsets,tables)
    _append_tables_locally_fill!(data,ptrs,offset,table)
  end

  rewind_ptrs!(ptrs)

  Table(data,ptrs)

end

function append_tables_locally(offsets::Tuple{}, tables::Tuple{})
  @unreachable "At least one table has to be provided"
end

function  _append_tables_locally_count!(ptrs,table)
  table_ptrs = table.ptrs
  n = length(table_ptrs)-1
  for i in 1:n
    ptrs[i+1] += table_ptrs[i+1]- table_ptrs[i]
  end
end

function _append_tables_locally_fill!(data,ptrs,offset,table)
  table_data = table.data
  table_ptrs = table.ptrs
  n = length(table_ptrs)-1
  for i in 1:n
    for j in table_ptrs[i]:(table_ptrs[i+1]-1)
      data[ptrs[i]] = table_data[j] + offset
      ptrs[i] += 1
    end
  end
end

"""
    collect1d(a)

Equivalent to

    [a[i] for in 1:length(a)]
"""
collect1d(a) = [a[i] for i in 1:length(a)]

function lazy_map(::typeof(getindex),a::Table,b::AbstractArray{<:Integer})
  LocalItemFromTable(a,b)
end

function get_local_item(a::Table,li::Integer)
  LocalItemFromTable(a,Fill(li,length(a)))
end

struct LocalItemFromTable{T,Vd,Vp,A} <: AbstractVector{T}
  a_to_lb_to_b::Table{T,Vd,Vp}
  a_to_lb::A
end

Base.size(m::LocalItemFromTable) = size(m.a_to_lb_to_b)

Base.IndexStyle(::Type{<:LocalItemFromTable}) = IndexLinear()

@propagate_inbounds function Base.getindex(m::LocalItemFromTable, a::Integer)
  p = m.a_to_lb_to_b.ptrs[a]-1
  m.a_to_lb_to_b.data[p+m.a_to_lb[a]]
end

"""
    find_local_index(a_to_b, b_to_la_to_a)
"""
function find_local_index(a_to_b, b_to_la_to_a)
  @notimplemented "find_local_index only implemented for table"
end

function find_local_index(a_to_b, b_to_la_to_a::Table)
  a_to_la = LocalIndexFromTable(a_to_b, b_to_la_to_a)
  a_to_la
end

struct LocalIndexFromTable{T,Vd,Vp,V<:AbstractVector} <: AbstractVector{T}
  a_to_b::V
  b_to_la_to_a::Table{T,Vd,Vp}
end

Base.size(m::LocalIndexFromTable) = size(m.a_to_b)

Base.IndexStyle(::Type{<:LocalIndexFromTable}) = IndexStyle(Table)

@propagate_inbounds function Base.getindex(m::LocalIndexFromTable{T}, a::Integer) where T
  b = m.a_to_b[a]
  pini = m.b_to_la_to_a.ptrs[b]
  pend = m.b_to_la_to_a.ptrs[b+1]-1
  la = zero(T)
  for (la,p) in enumerate(pini:pend)
    if a == m.b_to_la_to_a.data[p]
      return T(la)
    end
  end
  return T(UNSET)
end

"""
    flatten_partition(a_to_bs::Table,nb::Integer)
    flatten_partition(a_to_bs::Table)
"""
function flatten_partition(a_to_bs::Table,nb::Integer=maximum(a_to_bs.data))
  T = eltype(eltype(a_to_bs))
  b_to_a = zeros(T,nb)
  flatten_partition!(b_to_a,a_to_bs)
  b_to_a
end

function  flatten_partition!(b_to_a,a_to_bs::Table)
  for a in 1:length(a_to_bs)
    pini = a_to_bs.ptrs[a]
    pend = a_to_bs.ptrs[a+1]-1
    for p in pini:pend
      b = a_to_bs.data[p]
      b_to_a[b] = a
    end
  end
end

function to_dict(table::Table)
  dict = Dict{Symbol,Any}()
  dict[:data] = table.data
  dict[:ptrs] = table.ptrs
  dict
end

function from_dict(::Type{Table{T,Vd,Vp}}, dict::Dict{Symbol,Any}) where {T,Vd,Vp}
  data::Vd = dict[:data]
  ptrs::Vp = dict[:ptrs]
  Table(data,ptrs)
end

function Base.copy(a::Table)
  Table(copy(a.data),copy(a.ptrs))
end


text:
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using BenchmarkTools
using CUDA
using DataDeps
using Oceananigans
using Benchmarks
using Statistics
using Oceananigans.Solvers: sparse_approximate_inverse

function benchmark_spai_preconditioner(N, ε, nzrel, inverse)

    grid = RectilinearGrid(CPU(), size=(N, N, 1), extent=(1, 1, 1))

    model = HydrostaticFreeSurfaceModel(
                      grid = grid,
        momentum_advection = VectorInvariant(),
              free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver, preconditioner_method = nothing)
    )

    # to correctly create the matrix
    time_step!(model, 1)

    matrix = model.free_surface.implicit_step_solver.matrix_iterative_solver.matrix

    if !inverse
        trial = @benchmark begin
            CUDA.@sync blocking = true sparse_approximate_inverse($matrix, ε = $ε, nzrel = $nzrel)
        end samples = 5
    else
        trial = @benchmark begin
            CUDA.@sync blocking = true inv(Array($matrix))
        end samples = 5
    end

    return trial
end

N       = [64, 128, 256]
ε       = [0.1, 0.3]
nzrel   = [0.5, 1.0, 2.0]
inverse = (false, )

# Run and summarize benchmarks
print_system_info()
suite = run_benchmarks(benchmark_spai_preconditioner; N, ε, nzrel, inverse)

df = benchmarks_dataframe(suite)
benchmarks_pretty_table(df, title="SPAI preconditioner benchmarks")


text:
"""
    basemap(cmd0::String=""; kwargs...)

Plot base maps and frames.

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

Parameters
----------

- $(GMT.opt_J)
- $(GMT.opt_R)
- **A** : **polygon** : -- Str or [] --

    No plotting is performed. Instead, we determine the geographical coordinates of the polygon
    outline for the (possibly oblique) rectangular map domain.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#a)
- $(GMT.opt_B)
- **D** : **inset** : -- Str --

    Draw a simple map insert box on the map. Requires -F.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#d)
- **F** : **box** : -- Str --

    Without further options, draws a rectangular border around any map insert (D), map scale (L)
    or map rose (T)
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#f)
- $(GMT.opt_Jz)
- **L** : **map_scale** : -- Str --

    Draw a map scale.
    [`-L`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#l)
- $(GMT.opt_P)
- **Td** : **rose`** : -- Str --

    Draws a map directional rose on the map at the location defined by the reference and anchor points.
    [`-Td`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#t)
- **Tm** : **compass** : -- Str --

    Draws a map magnetic rose on the map at the location defined by the reference and anchor points.
    [`-Tm`](http://gmt.soest.hawaii.edu/doc/latest/psbasemap.html#t)
- $(GMT.opt_U)
- $(GMT.opt_V)
- $(GMT.opt_X)
- $(GMT.opt_Y)
- $(GMT.opt_bo)
- $(GMT.opt_f)
- $(GMT.opt_p)
- $(GMT.opt_t)
"""
function basemap(cmd0::String="", arg1=[]; K=false, O=false, first=true, kwargs...)

	length(kwargs) == 0 && return monolitic("psbasemap", cmd0, arg1)	# Speedy mode
	d = KW(kwargs)
	output, opt_T, fname_ext = fname_out(d)		# OUTPUT may have been an extension only

    cmd, opt_B, opt_J, opt_R = parse_BJR(d, "", "", O, " -JX12c/0")
	cmd, = parse_JZ(cmd, d)
	cmd  = parse_UVXY(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_p(cmd, d)
	cmd, = parse_t(cmd, d)
	cmd, = parse_bo(cmd, d)
	cmd  = parse_params(cmd, d)

	cmd, K, O, opt_B = set_KO(cmd, opt_B, first, K, O)		# Set the K O dance

	cmd = add_opt(cmd, 'A', d, [:A :polygon])
	cmd = add_opt(cmd, 'D', d, [:D :inset])
	cmd = add_opt(cmd, 'F', d, [:F :box])
	cmd = add_opt(cmd, 'L', d, [:L :map_scale])
	cmd = add_opt(cmd, "Td", d, [:Td :rose])
	cmd = add_opt(cmd, "Tm", d, [:Td :compass])

	cmd = finish_PS(d, cmd, output, K, O)
    return finish_PS_module(d, cmd, "", output, fname_ext, opt_T, K, "psbasemap", arg1)
end

# ---------------------------------------------------------------------------------------------------
basemap!(cmd0::String="", arg1=[]; K=true, O=true, first=false, kw...) =
	basemap(cmd0, arg1; K=K, O=O, first=first, kw...)
basemap(arg1=[]; K=false, O=false, first=true, kw...) = basemap("", arg1; K=K, O=O, first=first, kw...)
basemap!(arg1=[]; K=true, O=true, first=false, kw...) = basemap("", arg1; K=K, O=O, first=first, kw...)

psbasemap  = basemap 		# Alias
psbasemap! = basemap!		# Alias


text:
#
# Boundaries and Boundary Conditions
#

const coordinates = (:x, :y, :z)
const dims = length(coordinates)
const solution_fields = (:u, :v, :w, :T, :S)
const nsolution = length(solution_fields)

abstract type BCType end
struct Default <: BCType end

abstract type NonDefault <: BCType end
struct Flux <: NonDefault end
struct Gradient <: NonDefault end
struct Value <: NonDefault end

"""
    BoundaryCondition(BCType, condition)

Construct a boundary condition of `BCType` with `condition`,
where `BCType` is `Flux` or `Gradient`. `condition` may be a
number, array, or a function with signature:

    `condition(i, j, grid, t, iteration, u, v, w, T, S) = # function definition`

where `i` and `j` are indices along the boundary.
"""
struct BoundaryCondition{C<:BCType, T}
    condition :: T
end

# Constructors
BoundaryCondition(Tbc, c) = BoundaryCondition{Tbc, typeof(c)}(c)

Adapt.adapt_structure(to, b::BoundaryCondition{C, A}) where {C<:BCType, A<:AbstractArray} =
    BoundaryCondition(C, Adapt.adapt(to, parent(b.condition)))

DefaultBC() = BoundaryCondition{Default, Nothing}(nothing)

"""
    CoordinateBoundaryConditions(c)

Construct `CoordinateBoundaryCondition` to be applied along coordinate `c`, where
`c` is `:x`, `:y`, or `:z`. A CoordinateBoundaryCondition has two fields
`left` and `right` that store boundary conditions on the 'left' (negative side)
and 'right' (positive side) of a given coordinate.
"""
mutable struct CoordinateBoundaryConditions
     left :: BoundaryCondition
    right :: BoundaryCondition
end

CoordinateBoundaryConditions() = CoordinateBoundaryConditions(DefaultBC(), DefaultBC())

const CBC = CoordinateBoundaryConditions

#=
Here we overload setproperty! and getproperty to permit users to call
the 'right' and 'left' bcs in the z-direction 'bottom' and 'top'.

Note that 'right' technically corresponds to face point N+1. Thus
the fact that right == bottom is associated with the reverse z-indexing
convention. With ordinary indexing, right == top.
=#
Base.setproperty!(cbc::CBC, side::Symbol, bc) = setbc!(cbc, Val(side), bc)
setbc!(cbc::CBC, ::Val{S}, bc) where S = setfield!(cbc, S, bc)
setbc!(cbc::CBC, ::Val{:bottom}, bc) = setfield!(cbc, :right, bc)
setbc!(cbc::CBC, ::Val{:top}, bc) = setfield!(cbc, :left, bc)

Base.getproperty(cbc::CBC, side::Symbol) = getbc(cbc, Val(side))
getbc(cbc::CBC, ::Val{S}) where S = getfield(cbc, S)
getbc(cbc::CBC, ::Val{:bottom}) = getfield(cbc, :right)
getbc(cbc::CBC, ::Val{:top}) = getfield(cbc, :left)

"""
    FieldBoundaryConditions()

Construct `FieldBoundaryConditions` for a field.
A FieldBoundaryCondition has `CoordinateBoundaryConditions` in
`x`, `y`, and `z`.
"""
struct FieldBoundaryConditions <: FieldVector{dims, CoordinateBoundaryConditions}
    x :: CoordinateBoundaryConditions
    y :: CoordinateBoundaryConditions
    z :: CoordinateBoundaryConditions
end

"""
    BoundaryConditions()

Construct a boundary condition type full of default
`FieldBoundaryConditions` for u, v, w, T, S.
"""
struct ModelBoundaryConditions <: FieldVector{nsolution, FieldBoundaryConditions}
    u :: FieldBoundaryConditions
    v :: FieldBoundaryConditions
    w :: FieldBoundaryConditions
    T :: FieldBoundaryConditions
    S :: FieldBoundaryConditions
end

FieldBoundaryConditions() = FieldBoundaryConditions(CoordinateBoundaryConditions(),
                                                    CoordinateBoundaryConditions(),
                                                    CoordinateBoundaryConditions())

function ModelBoundaryConditions()
    bcs = (FieldBoundaryConditions() for i = 1:length(solution_fields))
    return ModelBoundaryConditions(bcs...)
end

#
# User API
#

#=
Notes:

- We assume that the boundary tendency has been previously calculated for
  a 'no-flux' boundary condition.

  This means that boudnary conditions take the form of
  an addition/subtraction to the tendency associated with a flux at point (A, A, I) below the bottom cell.
  This paradigm holds as long as consider boundary conditions on (A, A, C) variables only, where A is
  "any" of C or I.

 - We use the physics-based convention that

        flux = -κ * gradient,

    and that

        tendency = ∂ϕ/∂t = Gϕ = - ∇ ⋅ flux

=#

const BC = BoundaryCondition
const FBCs = FieldBoundaryConditions

# Do nothing in default case. These functions are called in cases where one of the
# z-boundaries is set, but not the other.
@inline apply_z_top_bc!(args...) = nothing
@inline apply_z_bottom_bc!(args...) = nothing

# These functions compute vertical fluxes for (A, A, C) quantities. They are not currently used.
@inline ∇κ∇ϕ_t(κ, ϕt, ϕt₋₁, flux, ΔzC, ΔzF) = (      -flux        - κ*(ϕt - ϕt₋₁)/ΔzC ) / ΔzF
@inline ∇κ∇ϕ_b(κ, ϕb, ϕb₊₁, flux, ΔzC, ΔzF) = ( κ*(ϕb₊₁ - ϕb)/ΔzC +       flux        ) / ΔzF

# Multiple dispatch on the type of boundary condition
getbc(bc::BC{C, <:Number}, args...)              where C = bc.condition
getbc(bc::BC{C, <:AbstractArray}, i, j, args...) where C = bc.condition[i, j]
getbc(bc::BC{C, <:Function}, args...)            where C = bc.condition(args...)

Base.getindex(bc::BC{C, <:AbstractArray}, inds...) where C = getindex(bc.condition, inds...)

"""
    apply_z_top_bc!(top_bc, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S)

Add the part of flux divergence associated with a top boundary condition on ϕ.
to Gϕ, where the conservation equation for ϕ is ∂ϕ/∂t = Gϕ.
If `top_bc.condition` is a function, the function must have the signature

    `top_bc.condition(i, j, grid, t, iteration, u, v, w, T, S)`

"""
@inline apply_z_top_bc!(top_flux::BC{<:Flux}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, 1] -= getbc(top_flux, i, j, grid, t, iteration, u, v, w, T, S) / grid.Δz

@inline apply_z_top_bc!(top_gradient::BC{<:Gradient}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, 1] += κ * getbc(top_gradient, i, j, grid, t, iteration, u, v, w, T, S) / grid.Δz

@inline apply_z_top_bc!(top_value::BC{<:Value}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, 1] += 2κ / grid.Δz * (
        getbc(top_value, i, j, grid, t, iteration, u, v, w, T, S) - ϕ[i, j, 1])

"""
    apply_z_bottom_bc!(bottom_bc, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S)

Add the part of flux divergence associated with a bottom boundary condition on ϕ.
to Gϕ, where the conservation equation for ϕ is ∂ϕ/∂t = Gϕ.
If `bottom_bc.condition` is a function, the function must have the signature

    `bottom_bc.condition(i, j, grid, t, iteration, u, v, w, T, S)`

"""
@inline apply_z_bottom_bc!(bottom_flux::BC{<:Flux}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, grid.Nz] += getbc(bottom_flux, i, j, grid, t, iteration, u, v, w, T, S) / grid.Δz

@inline apply_z_bottom_bc!(bottom_gradient::BC{<:Gradient}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, grid.Nz] -= κ * getbc(bottom_gradient, i, j, grid, t, iteration, u, v, w, T, S) / grid.Δz

@inline apply_z_bottom_bc!(bottom_value::BC{<:Value}, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S) =
    Gϕ[i, j, grid.Nz] -= 2κ / grid.Δz * (
        ϕ[i, j, grid.Nz] - getbc(bottom_value, i, j, grid, t, iteration, u, v, w, T, S))