text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#    Tests of table lines.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Horizontal lines" begin
    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
├────────┼────────┼────────┼────────┤
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
├────────┼────────┼────────┼────────┤
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(
        String,
        data;
        body_hlines = vcat(findall(x->x == true, data[:,2]))
    )
    @test result == expected

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
├........+........+........+........┤
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
├........+........+........+........┤
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""
    result = pretty_table(
        String,
        data;
        body_hlines = vcat(findall(x -> x == true, data[:, 2])),
        body_hlines_format = ('├','+','┤','.')
    )
    @test result == expected

    expected = """
┌───┬───────┬─────┬───┐
│ 1 │ false │ 1.0 │ 1 │
│ 2 │  true │ 2.0 │ 2 │
├───┼───────┼─────┼───┤
│ 3 │ false │ 3.0 │ 3 │
│ 4 │  true │ 4.0 │ 4 │
├───┼───────┼─────┼───┤
│ 5 │ false │ 5.0 │ 5 │
│ 6 │  true │ 6.0 │ 6 │
└───┴───────┴─────┴───┘
"""
    result = pretty_table(
        String,
        data;
        noheader = true,
        body_hlines = vcat(findall(x -> x == true, data[:, 2]))
    )
    @test result == expected

    # Test the case when `hlines` is a symbol
    # --------------------------------------------------------------------------

    expected = """
┌────────┬────────┬────────┬────────┐
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
├────────┼────────┼────────┼────────┤
│      1 │  false │    1.0 │      1 │
├────────┼────────┼────────┼────────┤
│      2 │   true │    2.0 │      2 │
├────────┼────────┼────────┼────────┤
│      3 │  false │    3.0 │      3 │
├────────┼────────┼────────┼────────┤
│      4 │   true │    4.0 │      4 │
├────────┼────────┼────────┼────────┤
│      5 │  false │    5.0 │      5 │
├────────┼────────┼────────┼────────┤
│      6 │   true │    6.0 │      6 │
└────────┴────────┴────────┴────────┘
"""

    result = pretty_table(String, data; hlines = :all)
    @test result == expected

    expected = """
│ Col. 1 │ Col. 2 │ Col. 3 │ Col. 4 │
│      1 │  false │    1.0 │      1 │
│      2 │   true │    2.0 │      2 │
│      3 │  false │    3.0 │      3 │
│      4 │   true │    4.0 │      4 │
│      5 │  false │    5.0 │      5 │
│      6 │   true │    6.0 │      6 │
"""

    result = pretty_table(String, data; hlines = :none)
    @test result == expected
end

@testset "Vertical lines" begin
    expected = """
────────────────────────
  C1     C2     C3   C4
 Int   Bool  Float  Hex
────────────────────────
   1  false    1.0    1
   2   true    2.0    2
   3  false    3.0    3
   4   true    4.0    4
   5  false    5.0    5
   6   true    6.0    6
────────────────────────
"""

    result1 = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = []
    )
    result2 = pretty_table(
        String, data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = :none
    )
    @test result1 == expected
    @test result2 == expected

    expected = """
┌────────────────────────┐
│  C1     C2     C3   C4 │
│ Int   Bool  Float  Hex │
├────────────────────────┤
│   1  false    1.0    1 │
│   2   true    2.0    2 │
│   3  false    3.0    3 │
│   4   true    4.0    4 │
│   5  false    5.0    5 │
│   6   true    6.0    6 │
└────────────────────────┘
"""

    result1 = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = [:begin, :end]
    )
    result2 = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = [0, 4]
    )
    @test result1 == expected
    @test result2 == expected

    expected = """
┌─────┬────────────────────────┐
│ Row │  C1     C2     C3   C4 │
│     │ Int   Bool  Float  Hex │
├─────┼────────────────────────┤
│   1 │   1  false    1.0    1 │
│   2 │   2   true    2.0    2 │
│   3 │   3  false    3.0    3 │
│   4 │   4   true    4.0    4 │
│   5 │   5  false    5.0    5 │
│   6 │   6   true    6.0    6 │
└─────┴────────────────────────┘
"""

    result = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = [:begin, 1 ,:end],
        show_row_number = true
    )
    @test result == expected

    expected = """
┌─────┬───────────┬────────────────────────┐
│ Row │ Row names │  C1     C2     C3   C4 │
│     │           │ Int   Bool  Float  Hex │
├─────┼───────────┼────────────────────────┤
│   1 │   Row 1   │   1  false    1.0    1 │
│   2 │   Row 2   │   2   true    2.0    2 │
│   3 │   Row 3   │   3  false    3.0    3 │
│   4 │   Row 4   │   4   true    4.0    4 │
│   5 │   Row 5   │   5  false    5.0    5 │
│   6 │   Row 6   │   6   true    6.0    6 │
└─────┴───────────┴────────────────────────┘
"""

    result = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        vlines = [:begin, 1, 2, :end],
        show_row_number = true,
        row_names = ["Row $i" for i in 1:6],
        row_name_column_title = "Row names",
        row_name_alignment = :c
    )
    @test result == expected

    expected = """
┌─────┬───────────┬────────────┬───
│ Row │ Row names │  C1     C2 │  ⋯
│     │           │ Int   Bool │  ⋯
├─────┼───────────┼────────────┼───
│   1 │   Row 1   │   1  false │  ⋯
│   2 │   Row 2   │   2   true │  ⋯
│   3 │   Row 3   │   3  false │  ⋯
│  ⋮  │     ⋮     │  ⋮     ⋮   │  ⋱
└─────┴───────────┴────────────┴───
       2 columns and 3 rows omitted
"""

    result = pretty_table(
        String,
        data;
        header = (
            ["C1",  "C2",   "C3",    "C4"],
            ["Int", "Bool", "Float", "Hex"]
        ),
        crop = :both,
        display_size = (12, 35),
        show_row_number = true,
        row_names = ["Row $i" for i in 1:6],
        row_name_column_title = "Row names",
        row_name_alignment = :c,
        vlines = [:begin, 1, 2, 4, :end]
    )
    @test result == expected
end


text:
module AROW

export fit, predict, composition, composition!

type Classifier
    dim::UInt64
    mean::Array{Float64, 1}
    cov::Array{Float64, 1}
    r::Float64

    function Classifier(dimension, param = 1.0)
        new(dimension, rand(dimension), rand(dimension), param)
    end
end

type MultiClassifier
    n_class::UInt64
    arows::Array{Classifier, 1}

    function MultiClassifier(dimension, num_classes, param = 1.0)
        arows = []
        for i=1:num_classes
            push!(arows, Classifier(dimension, param))
        end
        new(num_classes, arows)
    end
end

function fit{T<:AbstractArray}(classifier::MultiClassifier, x::T, label::Int)
    for i=1:classifier.n_class
        if i == label
            fit(classifier.arows[i], x, 1)
        elseif rand() <= 1.0 / (classifier.n_class - 1)
            fit(classifier.arows[i], x, -1)
        end
    end
end

function predict{T<:AbstractArray}(classifier::MultiClassifier, x::T)
    results = []

    for i=1:classifier.n_class
        score = dot(classifier.arows[i].mean, x)
        push!(results, score)
    end

    indmax(results)
end

function fit{T<:AbstractArray}(arow::Classifier, x::T, label::Int)
    assert(label == 1 || label == -1);

    margin = dot(arow.mean, x)

    if margin * label >= 1
        return 0
    end

    confidence = dot(arow.cov, x .* x)
    beta       = 1.0 / (confidence + arow.r)
    alpha      = (1.0 - label * margin) * beta;

    # update mean
    BLAS.axpy!(alpha * label, arow.cov .* x, arow.mean)

    # update covariance
    @inbounds for i=1:arow.dim
        arow.cov[i] = 1.0 / ((1.0 / arow.cov[i]) + x[i]^2 / arow.r)
    end

    ifelse(margin * label < 0, 1, 0)
end

function predict{T<:AbstractArray}(arow::Classifier, x::T)
    ifelse(dot(arow.mean, x) > 0, 1, -1)
end

function composition!(x::MultiClassifier, y::MultiClassifier)
    for i=1:x.n_class
        composition!(x.arows[i], y.arows[i])
    end
end

function composition!(x::Classifier, y::Classifier)
    for i=1:x.dim
        x.mean[i] = (x.mean[i] + y.mean[i]) / 2
        x.cov[i] = (x.cov[i] + y.cov[i]) / 2
    end

    return
end

end # module

ext: jl


text:
using BenchmarkTools, Plots, Plots.PlotMeasures, LaTeXStrings
using BenchmarkTools: minimum, median

SUITE = BenchmarkGroup()
model = "SPACECRAFT"
cases = ["NA01", "NA01-discrete",
         "A01", "A01-discrete",
         "A02", "A02-discrete",
         "A03", "A03-discrete",
         "A04", # "A04-discrete",
         #"A05", "A05-discrete",
         #"A06", "A06-discrete",
         #"A07", "A07-discrete",
         #"A08", "A08-discrete",
         "U01", # "U01-discrete",
         "U02", "U02-discrete"]

SUITE[model] = BenchmarkGroup()

include("Spacecraft.jl")
validation = []

boxdirs = BoxDirections{Float64, Vector{Float64}}(5)

# ----------------------------------------
#  NA01 (dense time)
# ----------------------------------------

prob_NA01 = spacecraft(abort_time=-1.)

sol_NA01 = solve(prob_NA01, alg=BOX(δ=0.04),
                 clustering_method=LazyClustering(),
                 intersection_method=TemplateHullIntersection(boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))
property = SR02_specification(sol_NA01)
push!(validation, Int(property))
SUITE[model][cases[1]] = @benchmarkable solve($prob_NA01, alg=BOX(δ=0.04),
                 clustering_method=LazyClustering(),
                 intersection_method=TemplateHullIntersection($boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))

# ----------------------------------------
#  NA01 (discrete time)
# ----------------------------------------

sol_NA01 = solve(prob_NA01, alg=BOX(δ=0.1, approx_model=NoBloating()),
                 clustering_method=LazyClustering(),
                 intersection_method=TemplateHullIntersection(boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))
property = SR02_specification(sol_NA01)
push!(validation, Int(property))
SUITE[model][cases[2]] = @benchmarkable solve($prob_NA01, alg=BOX(δ=0.1, approx_model=NoBloating()),
                 clustering_method=LazyClustering(),
                 intersection_method=TemplateHullIntersection($boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))

 # ----------------------------------------
 #  A01 (dense time)
 # ----------------------------------------

 prob_A01 = spacecraft(abort_time=120.)

 sol_A01 = solve(prob_A01, alg=BOX(δ=0.04),
                 clustering_method=LazyClustering(3),
                 intersection_method=TemplateHullIntersection(boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))
 property = SR02_specification(sol_A01)
 push!(validation, Int(property))
 SUITE[model][cases[3]] = @benchmarkable solve($prob_A01, alg=BOX(δ=0.04),
                  clustering_method=LazyClustering(3),
                  intersection_method=TemplateHullIntersection($boxdirs),
                  intersect_source_invariant=false,
                  tspan = (0.0 .. 300.0))

 # ----------------------------------------
 #  A01 (discrete time)
 # ----------------------------------------

 sol_A01 = solve(prob_A01, alg=BOX(δ=0.1, approx_model=NoBloating()),
                 clustering_method=LazyClustering(3),
                 intersection_method=TemplateHullIntersection(boxdirs),
                 intersect_source_invariant=false,
                 tspan = (0.0 .. 300.0))
 property = SR02_specification(sol_A01)
 push!(validation, Int(property))
 SUITE[model][cases[4]] = @benchmarkable solve($prob_A01, alg=BOX(δ=0.1, approx_model=NoBloating()),
                  clustering_method=LazyClustering(3),
                  intersection_method=TemplateHullIntersection($boxdirs),
                  intersect_source_invariant=false,
                  tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A02 (dense time)
# ----------------------------------------

prob_A02 = spacecraft(abort_time=[120., 125.])

sol_A02 = solve(prob_A02, alg=BOX(δ=0.04),
              clustering_method=LazyClustering(3),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A02)
push!(validation, Int(property))
SUITE[model][cases[5]] = @benchmarkable solve($prob_A02, alg=BOX(δ=0.04),
               clustering_method=LazyClustering(3),
               intersection_method=TemplateHullIntersection($boxdirs),
               intersect_source_invariant=false,
               tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A02 (discrete time)
# ----------------------------------------

sol_A02 = solve(prob_A02, alg=BOX(δ=0.1, approx_model=NoBloating()),
              clustering_method=LazyClustering(3),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A02)
push!(validation, Int(property))
SUITE[model][cases[6]] = @benchmarkable solve($prob_A02, alg=BOX(δ=0.1, approx_model=NoBloating()),
               clustering_method=LazyClustering(3),
               intersection_method=TemplateHullIntersection($boxdirs),
               intersect_source_invariant=false,
               tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A03 (dense time)
# ----------------------------------------

prob_A03 = spacecraft(abort_time=[120., 145.])

sol_A03 = solve(prob_A03, alg=BOX(δ=0.04),
                clustering_method=LazyClustering(4),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                intersect_source_invariant_method=TemplateHullIntersection(boxdirs),
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A03)
push!(validation, Int(property))
SUITE[model][cases[7]] = @benchmarkable solve($prob_A03, alg=BOX(δ=0.04),
              clustering_method=LazyClustering(4),
              intersection_method=TemplateHullIntersection($boxdirs),
              intersect_source_invariant=false,
              intersect_source_invariant_method=TemplateHullIntersection($boxdirs),
              tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A03 (discrete time)
# ----------------------------------------

sol_A03 = solve(prob_A03, alg=BOX(δ=0.1, approx_model=NoBloating()),
             clustering_method=LazyClustering(4),
             intersection_method=TemplateHullIntersection(boxdirs),
             intersect_source_invariant=false,
             intersect_source_invariant_method=TemplateHullIntersection(boxdirs),
             tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A03)
push!(validation, Int(property))
SUITE[model][cases[8]] = @benchmarkable solve($prob_A03, alg=BOX(δ=0.1, approx_model=NoBloating()),
              clustering_method=LazyClustering(4),
              intersection_method=TemplateHullIntersection($boxdirs),
              intersect_source_invariant=false,
              intersect_source_invariant_method=TemplateHullIntersection($boxdirs),
              tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A04 (dense time)
# ----------------------------------------

prob_A04 = spacecraft(abort_time=240.)

sol_A04 = solve(prob_A04, alg=BOX(δ=0.01),
              clustering_method=LazyClustering(40),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              intersect_source_invariant_method=TemplateHullIntersection(boxdirs),
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A04)
push!(validation, Int(property))
SUITE[model][cases[9]] = @benchmarkable solve($prob_A04, alg=BOX(δ=0.01),
            clustering_method=LazyClustering(40),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            intersect_source_invariant_method=TemplateHullIntersection($boxdirs),
            tspan = (0.0 .. 300.0))

# ----------------------------------------
#  A04 (discrete time)
# ----------------------------------------

#=
sol_A04 = solve(prob_A04, alg=BOX(δ=0.1, approx_model=NoBloating()),
           clustering_method=LazyClustering(16),
           intersection_method=TemplateHullIntersection(boxdirs),
           intersect_source_invariant=false,
           intersect_source_invariant_method=TemplateHullIntersection(boxdirs),
           tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A04)
push!(validation, Int(property))
SUITE[model][cases[10]] = @benchmarkable solve($prob_A04, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            intersect_source_invariant_method=TemplateHullIntersection($boxdirs),
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A05 (dense time)
# ----------------------------------------
#=
prob_A05 = spacecraft(abort_time=[235., 240.])

sol_A05 = solve(prob_A05, alg=BOX(δ=0.04),
              clustering_method=LazyClustering(13),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A05)
push!(validation, Int(property))
SUITE[model][cases[11]] = @benchmarkable solve($prob_A05, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(13),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A05 (discrete time)
# ----------------------------------------
#=
sol_A05 = solve(prob_A05, alg=BOX(δ=0.1, approx_model=NoBloating()),
           clustering_method=LazyClustering(13),
           intersection_method=TemplateHullIntersection(boxdirs),
           intersect_source_invariant=false,
           tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A05)
push!(validation, Int(property))
SUITE[model][cases[12]] = @benchmarkable solve($prob_A05, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(13),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A06 (dense time)
# ----------------------------------------
#=
prob_A06 = spacecraft(abort_time=[230., 240.])

sol_A06 = solve(prob_A06, alg=BOX(δ=0.04),
              clustering_method=LazyClustering(16),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A06)
push!(validation, Int(property))
SUITE[model][cases[13]] = @benchmarkable solve($prob_A06, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A06 (discrete time)
# ----------------------------------------
#=
sol_A06 = solve(prob_A06, alg=BOX(δ=0.1, approx_model=NoBloating()),
           clustering_method=LazyClustering(16),
           intersection_method=TemplateHullIntersection(boxdirs),
           intersect_source_invariant=false,
           tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A06)
push!(validation, Int(property))
SUITE[model][cases[14]] = @benchmarkable solve($prob_A06, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A07 (dense time)
# ----------------------------------------
#=
prob_A07 = spacecraft(abort_time=[50., 150.])

sol_A07 = solve(prob_A07, alg=BOX(δ=0.04),
              clustering_method=LazyClustering(50),
              intersection_method=TemplateHullIntersection(boxdirs),
              intersect_source_invariant=false,
              tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A07)
push!(validation, Int(property))
SUITE[model][cases[15]] = @benchmarkable solve($prob_A07, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(50),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A07 (discrete time)
# ----------------------------------------
#=
sol_A07 = solve(prob_A07, alg=BOX(δ=0.1, approx_model=NoBloating()),
           clustering_method=LazyClustering(50),
           intersection_method=TemplateHullIntersection(boxdirs),
           intersect_source_invariant=false,
           tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A07)
push!(validation, Int(property))
SUITE[model][cases[16]] = @benchmarkable solve($prob_A07, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(50),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A08 (dense time)
# ----------------------------------------
#=
prob_A08 = spacecraft(abort_time=[0., 240.])

sol_A08 = solve(prob_A08, alg=BOX(δ=0.04),
                clustering_method=LazyClustering(800),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A08)
push!(validation, Int(property))
SUITE[model][cases[17]] = @benchmarkable solve($prob_A08, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(800),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  A08 (discrete time)
# ----------------------------------------
#=
sol_A08 = solve(prob_A08, alg=BOX(δ=0.1, approx_model=NoBloating()),
                clustering_method=LazyClustering(900),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_A08)
push!(validation, Int(property))
SUITE[model][cases[18]] = @benchmarkable solve($prob_A08, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(900),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  U01 (dense time)
# ----------------------------------------

prob_U01 = spacecraft(abort_time=260.)

sol_U01 = solve(prob_U01, alg=BOX(δ=0.04),
                clustering_method=LazyClustering(16),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_U01)
push!(validation, Int(property))
SUITE[model][cases[10]] = @benchmarkable solve($prob_U01, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))

# ----------------------------------------
#  U01 (discrete time)
# ----------------------------------------

#=
sol_U01 = solve(prob_U01, alg=BOX(δ=0.1, approx_model=NoBloating()),
                clustering_method=LazyClustering(16),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_U01)
push!(validation, Int(property))
SUITE[model][cases[12]] = @benchmarkable solve($prob_U01, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))
=#

# ----------------------------------------
#  U02 (dense time)
# ----------------------------------------

prob_U02 = spacecraft(abort_time=[0., 260.])

sol_U02 = solve(prob_U02, alg=BOX(δ=0.04),
                clustering_method=LazyClustering(16),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_U02)
push!(validation, Int(property))
SUITE[model][cases[11]] = @benchmarkable solve($prob_U02, alg=BOX(δ=0.04),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))

# ----------------------------------------
#  U02 (discrete time)
# ----------------------------------------

sol_U02 = solve(prob_U02, alg=BOX(δ=0.1, approx_model=NoBloating()),
                clustering_method=LazyClustering(16),
                intersection_method=TemplateHullIntersection(boxdirs),
                intersect_source_invariant=false,
                intersect_source_invariant_method=TemplateHullIntersection(boxdirs),
                tspan = (0.0 .. 300.0))
property = SR02_specification(sol_U02)
push!(validation, Int(property))
SUITE[model][cases[12]] = @benchmarkable solve($prob_U02, alg=BOX(δ=0.1, approx_model=NoBloating()),
            clustering_method=LazyClustering(16),
            intersection_method=TemplateHullIntersection($boxdirs),
            intersect_source_invariant=false,
            tspan = (0.0 .. 300.0))

sol_NA01 = nothing
sol_A01 = nothing
sol_A02 = nothing
sol_A03 = nothing
sol_A04 = nothing
#sol_A05 = nothing
#sol_A06 = nothing
#sol_A07 = nothing
#sol_A08 = nothing
sol_U01 = nothing
sol_U02 = nothing
GC.gc()

# ==============================================================================
# Execute benchmarks and save benchmark results
# ==============================================================================

# tune parameters
tune!(SUITE)

# run the benchmarks
results = run(SUITE, verbose=true)

# return the sample with the smallest time value in each test
println("minimum time for each benchmark:\n", minimum(results))

# return the median for each test
println("median time for each benchmark:\n", median(results))

# export runtimes
runtimes = Dict()
for (i, c) in enumerate(cases)
    tm = median(results[model][c]).time * 1e-9
    runtimes[c] = round(tm, digits=4)
end

for (i, c) in enumerate(cases)
    print(io, "JuliaReach, $model, $c, $(validation[i]), $(runtimes[c])\n")
end

# ==============================================================================
# Plot
# ==============================================================================

function plot_SRNA01()
    prob_NA01 = spacecraft(abort_time=-1.)
    sol = solve(prob_NA01, alg=BOX(δ=0.04),
                     clustering_method=LazyClustering(),
                     intersection_method=TemplateHullIntersection(boxdirs),
                     intersect_source_invariant=false,
                     tspan = (0.0 .. 300.0))

    idx_approaching = findall(x -> x == 1, location.(sol))
    idx_attempt = findall(x -> x == 2, location.(sol))
    idx_aborting = findall(x -> x == 3, location.(sol))

    fig = plot(legend=:bottomright, tickfont=font(30, "Times"), guidefontsize=45,
               xlab=L"x", ylab=L"y",
               xtick=[-800, -600, -400, -200.0, 0.0], ytick=[-400, -300, -200, -100, 0.],
               xlims=(-1000.0, 0.0), ylims=(-450.0, 0.0),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
               size=(1000, 1000))

    for idx in idx_approaching
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1.)
    end
    for idx in idx_attempt
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:red, alpha=1.)
    end
    for idx in idx_aborting
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:cyan, alpha=1.)
    end
    fig
end

fig = plot_SRNA01()
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Spacecraft-SRNA01"))

# ----

function plot_SRA01()
    prob_A01 = spacecraft(abort_time=120.)
    @time sol = solve(prob_A01, alg=BOX(δ=0.04),
                      clustering_method=LazyClustering(3),
                      intersection_method=TemplateHullIntersection(boxdirs),
                      intersect_source_invariant=false,
                      tspan = (0.0 .. 300.0))

    idx_approaching = findall(x -> x == 1, location.(sol))
    idx_attempt = findall(x -> x == 2, location.(sol))
    idx_aborting = findall(x -> x == 3, location.(sol))

    fig = plot(legend=:bottomright, tickfont=font(30, "Times"), guidefontsize=45,
               xlab=L"x", ylab=L"y", lw=0.0,
               xtick=[-750, -500, -250, 0, 250.], ytick=[-400, -300, -200, -100, 0.],
               xlims=(-1000.0, 400.0), ylims=(-450.0, 0.0),
               bottom_margin=6mm, left_margin=2mm, right_margin=4mm, top_margin=3mm,
               size=(1000, 1000))

    for idx in idx_approaching
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:lightgreen, alpha=1.)
    end
    for idx in idx_attempt
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:red, alpha=1.)
    end
    for idx in idx_aborting
        plot!(fig, sol[idx], vars=(1, 2), lw=0.0, color=:cyan, alpha=1.)
    end
    fig
end

fig = plot_SRA01()
savefig(fig, joinpath(@__DIR__, "ARCH-COMP21-JuliaReach-Spacecraft-SRA01"))

GC.gc()


text:
## OPERATIONS ------------------------------------------------------------------
"""Check whether two `GeoArrays`s `a` and `b` are
geographically equal, although not necessarily in content."""
function equals(a::GeoArray, b::GeoArray)
    size(a) == size(b) && a.f == b.f && a.crs == b.crs
end

function Base.:-(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .- b.A, a.f, a.crs)
end

function Base.:+(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .+ b.A, a.f, a.crs)
end

function Base.:*(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A .* b.A, a.f, a.crs)
end

function Base.:/(a::GeoArray, b::GeoArray)
    equals(a, b) || throw(DimensionMismatch("Can't operate on non-geographic-equal `GeoArray`s"))
    GeoArray(a.A ./ b.A, a.f, a.crs)
end

Base.:+(a::GeoArray, b::Real) = GeoArray(a.A + b, a.f, a.crs)
Base.:-(a::GeoArray, b::Real) = GeoArray(a.A - b, a.f, a.crs)
Base.:*(a::GeoArray, b::Real) = GeoArray(a.A * b, a.f, a.crs)
Base.:/(a::GeoArray, b::Real) = GeoArray(a.A / b, a.f, a.crs)

## Math OPERATIONS -------------------------------------------------------------
function Base.:sum(ga::GeoArray, dims = 3)
    arr = sum(ga.A, dims = dims)
    GeoArray(arr, ga.f, ga.crs)
end

function Statistics.:mean(ga::GeoArray, dims = 3)
    arr = mean(ga.A; dims = dims)
    GeoArray(arr, ga.f, ga.crs)
end

function Base.:maximum(ga::GeoArray, dims = 3)
    arr = maximum(ga.A, dims = dims)
    GeoArray(arr, ga.f, ga.crs)
end

function Base.:minimum(ga::GeoArray, dims = 3)
    arr = minimum(ga.A, dims = dims)
    GeoArray(arr, ga.f, ga.crs)
end

export sum, mean, maximum, minimum


text:
"""
    IStringValued(nexid)

An interface for string valued widgets.
"""
struct IStringValued <: AbstractINextion
    _nid::NexID

    function IStringValued(nid::NexID)
        new(nid)
    end
end

function setproperty!(obj::IStringValued, property::Symbol, new_val)

    if property == :value
        setnexproperty!(NexID(obj), :txt, new_val)
    else
        setfield!(obj, property, new_val)
    end

end

function getproperty(obj::IStringValued, property::Symbol)

    if property == :value
        getnexproperty(NexID(obj), :txt, String)
    else
        getfield(obj, property)
    end

end


text:
using Random
using Statistics
using Test
using XiCor

@testset "Matches reference implementation in R" begin
    @test xicor(1:50, 1:50) ≈ 0.9411765 atol=0.00001
    @test xicor(1:100, 1:100) ≈ 0.970297 atol=0.00001
    @test xicor([1, 2, 3, 4, 5], [1, 1, 2, 2, 3]) ≈ 0.53125 atol=0.00001
    @test xicor([1, 2], [1, 0]) ≈ 0 atol=0.00001
    @test xicor([3.4, 2.3, 4.5, 1.2], [9.8, 7.6, 8.7, 6.5]) ≈ 0.2 atol=0.00001

    # The next test has ties in the x variable, which require randomization
    # The tolerance for matching against the R implementation is going to be
    # much more lenient for this particular test.
    x = div.(1:100, 10)
    y = rand(MersenneTwister(0), 100)
    @test mean([xicor(x, y, true, MersenneTwister(i)) for i in 1:100_000]) ≈ 0.01972526 atol=0.001
end


text:
"""
    PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,Parallel}

The coefficients of a derivative operator on a periodic grid.
"""
struct PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients} <: AbstractDerivativeCoefficients{T}
    # coefficients defining the operator and its action
    lower_coef::SVector{LowerOffset, T}
    central_coef::T
    upper_coef::SVector{UpperOffset, T}
    parallel::Parallel
    # corresponding orders etc.
    derivative_order::Int
    accuracy_order  ::Int
    symmetric       ::Bool
    source_of_coeffcients::SourceOfCoefficients

    function PeriodicDerivativeCoefficients(lower_coef::SVector{LowerOffset, T}, central_coef::T, upper_coef::SVector{UpperOffset, T},
                                            parallel::Parallel, derivative_order::Int, accuracy_order::Int,
                                            source_of_coeffcients::SourceOfCoefficients) where {T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients}
        symmetric = LowerOffset == UpperOffset
        if symmetric
            @inbounds for i in Base.OneTo(LowerOffset)
                symmetric = symmetric && lower_coef[i] == upper_coef[i]
            end
        end
        new{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients}(lower_coef, central_coef, upper_coef, parallel, derivative_order, accuracy_order, symmetric)
    end
end

"""
    mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients, u::AbstractVector, α, β)

Compute `α*D*u + β*dest` and store the result in `dest`.
"""
function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α, β) where {T,LowerOffset,UpperOffset}
    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > LowerOffset+UpperOffset DimensionMismatch
    end

    @unpack lower_coef, central_coef, upper_coef, parallel = coefficients
    convolve_periodic_boundary_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, β, parallel)
end

"""
    mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients, u::AbstractVector, α)

Compute `α*D*u` and store the result in `dest`.
"""
function mul!(dest::AbstractVector, coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset}, u::AbstractVector, α) where {T,LowerOffset,UpperOffset}
    @boundscheck begin
        @argcheck length(dest) == length(u) DimensionMismatch
        @argcheck length(u) > LowerOffset+UpperOffset DimensionMismatch
    end

    @unpack lower_coef, central_coef, upper_coef, parallel = coefficients
    convolve_periodic_boundary_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α)
    convolve_interior_coefficients!(dest, lower_coef, central_coef, upper_coef, u, α, parallel)
end

@generated function convolve_periodic_boundary_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                             u::AbstractVector, α, β) where {LowerOffset, UpperOffset}
    ex_lower = :( nothing )
    for i in 1:LowerOffset
        ex = :( lower_coef[$LowerOffset]*u[end-$(LowerOffset-i)] )
        for j in LowerOffset-1:-1:1
            if i-j < 1
                ex = :( $ex + lower_coef[$j]*u[end-$(j-i)] )
            else
                ex = :( $(ex) + lower_coef[$j]*u[$(i-j)] )
            end
        end
        ex = :( $ex + central_coef*u[$i] )
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[$(i+j)] )
        end
        ex_lower = quote
            $ex_lower
            @inbounds dest[$i] = β*dest[$i] + α*$ex
        end
    end

    ex_upper = :( nothing )
    for i in (UpperOffset-1):-1:0
        ex = :( lower_coef[$LowerOffset]*u[end-$(i+LowerOffset)] )
        for j in LowerOffset-1:-1:1
            ex = :( $ex + lower_coef[$j]*u[end-$(j+i)] )
        end
        ex = :( $ex + central_coef*u[end-$i] )
        for j in 1:UpperOffset
            if i-j < 0
                ex = :( $ex + upper_coef[$j]*u[$(j-i)] )
            else
                ex = :( $ex + upper_coef[$j]*u[end-$(i-j)] )
            end
        end
        ex_upper = quote
            $ex_upper
            @inbounds dest[end-$i] = β*dest[end-$i] + α*$ex
        end
    end

    quote
        $ex_lower
        $ex_upper
    end
end

@generated function convolve_periodic_boundary_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                             u::AbstractVector, α) where {LowerOffset, UpperOffset}
    ex_lower = :( nothing )
    for i in 1:LowerOffset
        ex = :( lower_coef[$LowerOffset]*u[end-$(LowerOffset-i)] )
        for j in LowerOffset-1:-1:1
            if i-j < 1
                ex = :( $ex + lower_coef[$j]*u[end-$(j-i)] )
            else
                ex = :( $(ex) + lower_coef[$j]*u[$(i-j)] )
            end
        end
        ex = :( $ex + central_coef*u[$i] )
        for j in 1:UpperOffset
            ex = :( $ex + upper_coef[$j]*u[$(i+j)] )
        end
        ex_lower = quote
            $ex_lower
            @inbounds dest[$i] = α*$ex
        end
    end

    ex_upper = :( nothing )
    for i in (UpperOffset-1):-1:0
        ex = :( lower_coef[$LowerOffset]*u[end-$(i+LowerOffset)] )
        for j in LowerOffset-1:-1:1
            ex = :( $ex + lower_coef[$j]*u[end-$(j+i)] )
        end
        ex = :( $ex + central_coef*u[end-$i] )
        for j in 1:UpperOffset
            if i-j < 0
                ex = :( $ex + upper_coef[$j]*u[$(j-i)] )
            else
                ex = :( $ex + upper_coef[$j]*u[end-$(i-j)] )
            end
        end
        ex_upper = quote
            $ex_upper
            @inbounds dest[end-$i] = α*$ex
        end
    end

    quote
        $ex_lower
        $ex_upper
    end
end

@generated function convolve_interior_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                    u::AbstractVector, α, β, parallel) where {LowerOffset, UpperOffset}
    ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] )
    for j in LowerOffset-1:-1:1
        ex = :( $ex + lower_coef[$j]*u[i-$j] )
    end
    ex = :( $ex + central_coef*u[i] )
    for j in 1:UpperOffset
        ex = :( $ex + upper_coef[$j]*u[i+$j] )
    end

    quote
        @inbounds for i in $(LowerOffset+1):(length(dest)-$UpperOffset)
            dest[i] = β*dest[i] + α*$ex
        end
    end
end

function convolve_interior_coefficients!(dest::AbstractVector{T}, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                         u::AbstractVector, α, β, parallel::Val{:threads}) where {T, LowerOffset, UpperOffset}
    Threads.@threads for i in (LowerOffset+1):(length(dest)-UpperOffset) @inbounds begin
        tmp = zero(T)
        for j in Base.OneTo(LowerOffset)
            tmp += lower_coef[j]*u[i-j]
        end
        tmp += central_coef*u[i]
        for j in Base.OneTo(UpperOffset)
            tmp += upper_coef[j]*u[i+j]
        end
        dest[i] = β*dest[i] + α*tmp
    end end
end

@generated function convolve_interior_coefficients!(dest::AbstractVector, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                                    u::AbstractVector, α, parallel) where {LowerOffset, UpperOffset}
    ex = :( lower_coef[$LowerOffset]*u[i-$LowerOffset] )
    for j in LowerOffset-1:-1:1
        ex = :( $ex + lower_coef[$j]*u[i-$j] )
    end
    ex = :( $ex + central_coef*u[i] )
    for j in 1:UpperOffset
        ex = :( $ex + upper_coef[$j]*u[i+$j] )
    end

    quote
        @inbounds for i in $(LowerOffset+1):(length(dest)-$UpperOffset)
            dest[i] = α*$ex
        end
    end
end

function convolve_interior_coefficients!(dest::AbstractVector{T}, lower_coef::SVector{LowerOffset}, central_coef, upper_coef::SVector{UpperOffset},
                                         u::AbstractVector, α, parallel::Val{:threads}) where {T, LowerOffset, UpperOffset}
    Threads.@threads for i in (LowerOffset+1):(length(dest)-UpperOffset) @inbounds begin
        tmp = zero(T)
        for j in Base.OneTo(LowerOffset)
            tmp += lower_coef[j]*u[i-j]
        end
        tmp += central_coef*u[i]
        for j in Base.OneTo(UpperOffset)
            tmp += upper_coef[j]*u[i+j]
        end
        dest[i] = α*tmp
    end end
end

"""
    BeljaddLeFlochMishraParés2017

Coefficients of the periodic operators given in
  Beljadid, LeFloch, Mishra, Parés (2017)
  Schemes with Well-Controlled Dissipation. Hyperbolic Systems in
    Nonconservative Form.
  Communications in Computational Physics 21.4, pp. 913-946.
"""
struct BeljaddLeFlochMishraParés2017 <: SourceOfCoefficients end

function Base.show(io::IO, ::BeljaddLeFlochMishraParés2017)
    print(io,
        "  Beljadid, LeFloch, Mishra, Parés (2017) \n",
        "  Schemes with Well-Controlled Dissipation. Hyperbolic Systems in \n",
        "    Nonconservative Form. \n",
        "  Communications in Computational Physics 21.4, pp. 913-946. \n")
end

"""
    periodic_central_derivative_coefficients(derivative_order, accuracy_order, T=Float64, parallel=Val{:serial}())

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_central_derivative_coefficients(derivative_order, accuracy_order, T=Float64, parallel=Val{:serial}())
    @argcheck derivative_order > 0
    @argcheck accuracy_order > 0
    @argcheck typeof(parallel) <: Union{Val{:serial}, Val{:threads}}

    if isodd(accuracy_order)
        @warn("Central derivative operators support only even orders of accuracy.")
    end

    if derivative_order == 1
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if iseven(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            coef[j] =  tmp / j
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = zero(T)
        lower_coef = -upper_coef
        source = BeljaddLeFlochMishraParés2017()
    elseif derivative_order == 2
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if iseven(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            coef[j] =  2*tmp / j^2
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = -2*sum(upper_coef)
        lower_coef = upper_coef
        source = BeljaddLeFlochMishraParés2017()
    elseif derivative_order == 3
        # Exact evaluation of the coefficients, see
        # Beljadid, LeFloch, Mishra, Parés (2017)
        # Schemes with Well-Controlled Dissipation. Hyperbolic Systems in Nonconservative Form.
        # Communications in Computational Physics 21.4, pp. 913-946.
        n = (accuracy_order+1) ÷ 2
        coef = Vector{T}(undef, n)
        for j in 1:n
            tmp = one(Rational{Int128})
            if isodd(j)
                tmp *= -1
            end
            for i in 1:j
                tmp *= (n+1-i) // (n+i)
            end
            fac = zero(Rational{Int128})
            for k in 1:n
                k == j && continue
                fac += 1 // k^2
            end
            coef[j] =  6*tmp*fac/ j
        end
        upper_coef = SVector{n,T}(coef)
        central_coef = zero(T)
        lower_coef = -upper_coef
        source = BeljaddLeFlochMishraParés2017()
    else
        throw(ArgumentError("Derivative order $derivative_order not implemented yet."))
    end

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, parallel, derivative_order, accuracy_order, source)
end

"""
    Fornberg1998

Coefficients of the periodic operators given in
  Fornberg (1998)
  Calculation of Weights in Finite Difference Formulas.
  SIAM Rev. 40.3, pp. 685-691.
"""
struct Fornberg1998 <: SourceOfCoefficients end

function Base.show(io::IO, ::Fornberg1998)
    print(io,
        "  Fornberg (1998) \n",
        "  Calculation of Weights in Finite Difference Formulas. \n",
        "  SIAM Rev. 40.3, pp. 685-691. \n")
end

"""
    fornberg(x::Vector{T}, m::Int) where {T}

Calculate the weights of a finite difference approximation of the `m`th derivative
with maximal order of accuracy at `0` using the nodes `x`, see
Fornberg (1998)
Calculation of Weights in Finite Difference Formulas
SIAM Rev. 40.3, pp. 685-691.
"""
function fornberg(xx::Vector{T}, m::Int) where {T}
    x = sort(xx)
    @argcheck x[1] <= 0
    @argcheck x[end] >= 0
    @argcheck length(unique(x)) == length(x)
    @argcheck zero(T) in x
    @argcheck m >= 1

    z = zero(T)
    n = length(x) - 1
    c = fill(zero(T), length(x), m+1)
    c1 = one(T)
    c4 = x[1] - z
    c[1,1] = one(T)
    for i in 1:n
        mn = min(i,m)
        c2 = one(T)
        c5 = c4
        c4 = x[i+1] - z
        for j in 0:i-1
            c3 = x[i+1] - x[j+1]
            c2 = c2 * c3
            if j == i-1
                for k in mn:-1:1
                    c[i+1,k+1] = c1 * (k*c[i,k]-c5*c[i,k+1]) / c2
                end
                c[i+1,1] = -c1*c5*c[i,1] / c2
            end
            for k in mn:-1:1
                c[j+1,k+1] = (c4*c[j+1,k+1]-k*c[j+1,k]) / c3
            end
            c[j+1,1] = c4*c[j+1,1] / c3
        end
        c1 = c2
    end

    c[:,end]
end

"""
    function periodic_derivative_coefficients(derivative_order, accuracy_order, left_offset=-(accuracy_order+1)÷2, T=Float64, parallel=Val{:serial}())

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T` where
the leftmost grid point used is determined by `left_offset`.
The evaluation of the derivative can be parallised using threads by chosing
parallel=Val{:threads}())`.
"""
function periodic_derivative_coefficients(derivative_order, accuracy_order, left_offset::Int=-(accuracy_order+1)÷2, T=Float64, parallel=Val{:serial}())
    @argcheck -accuracy_order <= left_offset <= 0

    x = Rational{Int128}.(left_offset:left_offset+accuracy_order)
    c = fornberg(x, derivative_order)

    z = zero(eltype(x))
    lower_idx = x .< z
    central_idx = first(findall(xx->xx==z, x))
    upper_idx = x .> z

    LowerOffset = sum(lower_idx)
    UpperOffset = sum(upper_idx)

    lower_coef = SVector{LowerOffset, T}(T.(reverse(c[lower_idx])))
    central_coef = T(c[central_idx])
    upper_coef = SVector{UpperOffset, T}(T.(c[upper_idx]))

    source = Fornberg1998()

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, parallel, derivative_order, accuracy_order, source)
end

"""
    Holoborodko2008

Coefficients of the periodic operators given in
  Holoborodko (2008)
  Smooth Noise Robust Differentiators.
  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/
"""
struct Holoborodko2008 <: SourceOfCoefficients end

function Base.show(io::IO, ::Holoborodko2008)
    print(io,
        "  Holoborodko (2008) \n",
        "  Smooth Noise Robust Differentiators. \n",
        "  http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/ \n")
end

"""
    periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
                                     T=Float64, parallel=Val{:serial}(),
                                     stencil_width=accuracy_order+3)

Create the `PeriodicDerivativeCoefficients` approximating the `derivative_order`-th
derivative with an order of accuracy `accuracy_order` and scalar type `T` given
by Holoborodko2008.
The evaluation of the derivative can be parallised using threads by chosing
parallel=Val{:threads}())`.
"""
function periodic_derivative_coefficients(source::Holoborodko2008, derivative_order, accuracy_order;
                                          T=Float64, parallel=Val{:serial}(),
                                          stencil_width=accuracy_order+3)
    method_exists = true
    if derivative_order == 1
        if accuracy_order == 2
            if stencil_width == 5
                lower_coef = SVector{2, T}([
                    -2//8, -1//8
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 7
                lower_coef = SVector{3, T}([
                    -5//32, -4//32, -1//32
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -14//128, -14//128, -6//128, -1//128
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 11
                lower_coef = SVector{5, T}([
                    -42//512, -48//512, -27//512, -8//512, -1//512
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            else
                method_exists = false
            end
        elseif accuracy_order == 4
            if stencil_width == 7
                lower_coef = SVector{3, T}([
                    -39//96, -12//96, 5//96
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -27//96, -16//96, 1//96, 2//96
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            elseif stencil_width == 11
                lower_coef = SVector{5, T}([
                    -322//1536, -256//1536, -39//1536, 32//1536, 11//1536
                ])
                upper_coef = -lower_coef
                central_coef = T(0)
            else
                method_exists = false
            end
        else
            method_exists = false
        end
    elseif derivative_order == 2
        if accuracy_order == 2
            if stencil_width == 5
                lower_coef = SVector{2, T}([
                    0, 1//4
                ])
                upper_coef = lower_coef
                central_coef = T(-2//4)
            elseif stencil_width == 7
                lower_coef = SVector{3, T}([
                    -1//16, 2//16, 1//16
                ])
                upper_coef = lower_coef
                central_coef = T(-4//16)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -4//64, 4//64, 4//64, 1//64
                ])
                upper_coef = lower_coef
                central_coef = T(-10//64)
            else
                method_exists = false
            end
        elseif accuracy_order == 4
            if stencil_width == 7
                lower_coef = SVector{3, T}([
                    1//12, 5//12, -1//12
                ])
                upper_coef = lower_coef
                central_coef = T(-10//12)
            elseif stencil_width == 9
                lower_coef = SVector{4, T}([
                    -12//192, 52//192, 12//192, -7//192
                ])
                upper_coef = lower_coef
                central_coef = T(-90//192)
            else
                method_exists = false
            end
        else
            method_exists = false
        end
    else
        method_exists = false
    end
    if method_exists == false
        throw(ArgumentError("Method with derivative_order=$derivative_order, accuracy_order=$accuracy_order, stencil_width=$stencil_width not implemented/derived."))
    end

    PeriodicDerivativeCoefficients(lower_coef, central_coef, upper_coef, parallel, derivative_order, accuracy_order, source)
end

"""
    PeriodicDerivativeOperator

A derivative operator on a uniform periodic grid.
"""
struct PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid} <: AbstractPeriodicDerivativeOperator{T}
    coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients}
    grid_compute::Grid   # N-1 nodes, including the left and excluding the right boundary
    grid_evaluate::Grid #  N  nodes, including both boundaries
    Δx::T
    factor::T

    function PeriodicDerivativeOperator(coefficients::PeriodicDerivativeCoefficients{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients},
                                        grid_evaluate::Grid) where {T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}
        @argcheck length(grid_evaluate) > LowerOffset+UpperOffset DimensionMismatch
        grid_compute = range(grid_evaluate[1], stop=grid_evaluate[end-1], length=length(grid_evaluate)-1)

        Δx = step(grid_evaluate)
        factor = inv(Δx^coefficients.derivative_order)
        new{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}(coefficients, grid_compute, grid_evaluate, Δx, factor)
    end
end

@inline source_of_coeffcients(D::PeriodicDerivativeOperator) = source_of_coeffcients(D.coefficients)

function Base.show(io::IO, D::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset}) where {T,LowerOffset,UpperOffset}
    x = D.grid_evaluate
    if derivative_order(D) == 1
        print(io, "Periodic 1st derivative operator of order ")
    elseif  derivative_order(D) == 2
        print(io, "Periodic 2nd derivative operator of order ")
    elseif  derivative_order(D) == 3
        print(io, "Periodic 3rd derivative operator of order ")
    else
        print(io, "Periodic ", derivative_order(D), "th derivative operator of order ")
    end
    print(io, accuracy_order(D), " {T=", T, ", Parallel=", typeof(D.coefficients.parallel), "} \n")
    print(io, "on a grid in [", first(x), ", ", last(x), "] using ", length(x), " nodes, \n")
    print(io, "stencils with ", LowerOffset, " nodes to the left, ", UpperOffset,
                " nodes to the right, and coefficients from \n", source_of_coeffcients(D))
end

"""
    mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α, β)

Compute `α*D*u + β*dest` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α, β)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor, β)
end

"""
    mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α)

Compute `α*D*u` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, D::PeriodicDerivativeOperator, u::AbstractVector, α)
    @boundscheck begin
        @argcheck size(D, 2) == length(u) DimensionMismatch
        @argcheck size(D, 1) == length(dest) DimensionMismatch
    end
    @inbounds mul!(dest, D.coefficients, u, α*D.factor)
end

"""
    integrate(func, u, D::PeriodicDerivativeOperator)

Map the function `func` to the coefficients `u` and integrate with respect to
the quadrature rule associated with the periodic derivative operator `D`.
"""
function integrate(func, u::AbstractVector, D::PeriodicDerivativeOperator)
    @boundscheck begin
        length(u) == length(grid(D))
    end
    @unpack Δx = D

    @inbounds res = sum(func, u)

    Δx * res
end

"""
    periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, parallel=Val{:serial}())

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_central_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, parallel=Val{:serial}())
    grid = range(xmin, stop=xmax, length=N) # N includes two identical boundary nodes
    coefficients = periodic_central_derivative_coefficients(derivative_order, accuracy_order, eltype(grid), parallel)
    PeriodicDerivativeOperator(coefficients, grid)
end

"""
    periodic_central_derivative_operator(derivative_order, accuracy_order, grid, parallel=Val{:serial}())

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on the uniform `grid` up to order of accuracy `accuracy_order`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_central_derivative_operator(derivative_order, accuracy_order, grid::Union{LinRange,StepRangeLen}, parallel=Val{:serial}())
    coefficients = periodic_central_derivative_coefficients(derivative_order, accuracy_order, eltype(grid), parallel)
    PeriodicDerivativeOperator(coefficients, grid)
end

"""
    periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, left_offset=-(accuracy_order+1)÷2, parallel=Val{:serial}())

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order` where the leftmost grid point used is
determined by `left_offset`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_derivative_operator(derivative_order::Int, accuracy_order, xmin, xmax, N, left_offset::Int=-(accuracy_order+1)÷2,
                                      parallel::Union{Val{:serial},Val{:threads}}=Val{:serial}())
    grid = range(xmin, stop=xmax, length=N) # N includes two identical boundary nodes
    coefficients = periodic_derivative_coefficients(derivative_order, accuracy_order, left_offset, eltype(grid), parallel)
    PeriodicDerivativeOperator(coefficients, grid)
end

@inline function periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N,
                                              parallel::Union{Val{:serial},Val{:threads}}, left_offset::Int=-(accuracy_order+1)÷2)
    periodic_derivative_operator(derivative_order, accuracy_order, xmin, xmax, N, left_offset, parallel)
end

"""
    periodic_derivative_operator(source::Holoborodko2008, derivative_order, accuracy_order,
                                 xmin, xmax, N; parallel=Val{:serial}(), kwargs...)

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on a uniform grid between `xmin` and `xmax` with `N` grid points up
to order of accuracy `accuracy_order` where the leftmost grid point used is
determined by `left_offset`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_derivative_operator(source::Holoborodko2008, derivative_order, accuracy_order,
                                      xmin, xmax, N; parallel=Val{:serial}(), kwargs...)
    grid = range(xmin, stop=xmax, length=N) # N includes two identical boundary nodes
    coefficients = periodic_derivative_coefficients(source, derivative_order, accuracy_order;
                                                    kwargs..., T=eltype(grid), parallel=parallel)
    PeriodicDerivativeOperator(coefficients, grid)
end

"""
    periodic_derivative_operator(derivative_order, accuracy_order, grid, left_offset=-(accuracy_order+1)÷2, parallel=Val{:serial}())

Create a `PeriodicDerivativeOperator` approximating the `derivative_order`-th
derivative on thr uniform `grid` up to order of accuracy `accuracy_order` where
the leftmost grid point used is determined by `left_offset`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function periodic_derivative_operator(derivative_order, accuracy_order, grid::Union{LinRange,StepRangeLen},
                                      left_offset::Int=-(accuracy_order+1)÷2, parallel=Val{:serial}())
    coefficients = periodic_derivative_coefficients(derivative_order, accuracy_order, left_offset, eltype(grid), parallel)
    PeriodicDerivativeOperator(coefficients, grid)
end

@inline function periodic_derivative_operator(derivative_order, accuracy_order, grid::Union{LinRange,StepRangeLen},
                                              parallel::Union{Val{:serial},Val{:threads}}, left_offset::Int=-(accuracy_order+1)÷2)
    periodic_derivative_operator(derivative_order, accuracy_order, grid, left_offset, parallel)
end

"""
    PeriodicDissipationOperator

A dissipation operator on a periodic finite difference grid.
"""
struct PeriodicDissipationOperator{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid} <: AbstractPeriodicDerivativeOperator{T}
    factor::T
    Di::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}

    function PeriodicDissipationOperator(factor::T, Di::PeriodicDerivativeOperator{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}) where {T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}
        # check validity
        @argcheck iseven(derivative_order(Di))
        @argcheck (-1)^(1+derivative_order(Di)÷2)*factor >= 0
        new{T,LowerOffset,UpperOffset,Parallel,SourceOfCoefficients,Grid}(factor, Di)
    end
end

grid(Di::PeriodicDissipationOperator) = grid(Di.Di)
derivative_order(Di::PeriodicDissipationOperator) = derivative_order(Di.Di)
accuracy_order(Di::PeriodicDissipationOperator) = accuracy_order(Di.Di)
source_of_coeffcients(Di::PeriodicDissipationOperator) = MattssonSvärdNordström2004()

function Base.show(io::IO, Di::PeriodicDissipationOperator{T}) where {T}
    if  derivative_order(Di) == 2
        print(io, "SBP 2nd derivative dissipation operator of order ")
    else
        print(io, "SBP ", derivative_order(Di), "th derivative dissipation operator of order ")
    end
    x = Di.Di.grid_evaluate
    print(io, accuracy_order(Di), " {T=", T, ", Parallel=", typeof(Di.Di.coefficients.parallel), "} \n")
    print(io, "on a grid in [", first(x), ", ", last(x), "] using ", length(x), " nodes \n")
    print(io, "and coefficients given in \n")
    print(io, source_of_coeffcients(Di))
end

"""
    mul!(dest::AbstractVector, D::PeriodicDissipationOperator, u::AbstractVector, α, β)

Compute `α*D*u + β*dest` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                       u::AbstractVector, α, β)
    mul!(dest, Di.Di, u, Di.factor*α, β)
end

"""
    mul!(dest::AbstractVector, D::PeriodicDissipationOperator, u::AbstractVector, α)

Compute `α*D*u` and store the result in `dest`.
"""
Base.@propagate_inbounds function mul!(dest::AbstractVector, Di::PeriodicDissipationOperator,
                                       u::AbstractVector, α)
    mul!(dest, Di.Di, u, Di.factor*α)
end

"""
    dissipation_operator(D::PeriodicDerivativeOperator;
                         strength=one(eltype(D)),
                         order=accuracy_order(D),
                         parallel=D.coefficients.parallel)

Create a negative semidefinite `DissipationOperator` using undivided differences
approximating a `order`-th derivative with strength `strength` adapted to the
derivative operator `D`.
The evaluation of the derivative can be parallised using threads by chosing
`parallel=Val{:threads}())`.
"""
function dissipation_operator(D::PeriodicDerivativeOperator;
                              strength=eltype(D)(1),
                              order=accuracy_order(D),
                              parallel=D.coefficients.parallel)
    @argcheck iseven(order) ArgumentError("Dissipation operators require even derivatives.")
    # account for the grid spacing
    factor = strength * (-1)^(1 + order÷2) * D.Δx^order
    x = D.grid_evaluate
    Di = periodic_derivative_operator(order, order, first(x), last(x), length(x), parallel)
    PeriodicDissipationOperator(factor, Di)
end


text:
typealias NSE Union{Number, Symbol, Expr}
typealias SE Union{Symbol, Expr}

type OdeModel
  file::AbstractString

  states::Dict{AbstractString, Float64}
  parameters::Dict{AbstractString, Float64}
  odes::Dict{AbstractString, NSE}
end

OdeModel(file::AbstractString) = parse_model(file::AbstractString)


text:
using Documenter, WAHVectors

makedocs(
    format = :html,
    modules = [WAHVectors],
    sitename = "WAHVectors",
    doctest = false,
    strict = false,
    pages = [
        "Home" => "index.md",
        "API" => [
            "WAHElements" => "man/api_elements.md"
        ]
    ],
    authors = "Ben J. Ward"
)
deploydocs(repo = "github.com/BioJulia/WAHVectors.jl.git")


text:
using SafeTestsets

@safetestset "PF00011" begin include("PF00011.jl") end


text:
# struct Diagonal{T, A <: AbstractVector{T}} <: AbstractArray{T}
#     diag::A
#
#     function Diagonal{T, A}(diag) where {T, A <: AbstractArray{T}}
#         Base.require_one_based_indexing(diag)
#         return new{T, A}(diag)
#     end
# end
# Diagonal(a::AbstractArray{T, N}) where {T, N} = Diagonal{T, typeof(a)}(a)
# Diagonal{T, N}(a::AbstractArray) where {T, N} = Diagonal(convert(AbstractArray{T, N}, a)::AbstractArray{T, N})

function LinearAlgebra.diagind(A::AbstractArray)
    return CartesianIndex{ndims(A)}[i for i in CartesianIndices(A) if allequal(i.I)]
end

LinearAlgebra.diag(A::AbstractArray) = A[diagind(A)]

using Einsum

#=
If I were to extend LinearAlgebra.diag for an arbitrary dimension, what would it look like?  Would it be all of the indices that are the same (thus returning a vector)?  I.e., (1, 1, …, 1), (2, 2, …, 2), …, (n, n, …, n)?  Or would it be a diagonal hyperplane that is one dimension less than the given matrix?
=#

"""
```julia
diag_alt(A::AbstractArray)
```

Given an array of arbitrary dimesions, select all of the indices that are the same, thus returning a vector.  E.g.:
```julia
julia> A = rand(Int8, 3, 3)
3×3 Array{Int8,2}:
 93  123  -32
 71  -68  -75
 99   -7  -57

julia> diag_alt(A)
3-element Array{Int8,1}:
  93
 -68
 -57

julia> B = rand(Int8, 3, 3, 3)
3×3×3 Array{Int8,3}:
[:, :, 1] =
  16   -23    53
  13    -7  -115
 -96  -105    29

[:, :, 2] =
 -80  17    53
  88  40    15
  37  51  -109

[:, :, 3] =
  -81    88   37
   38  -115  -29
 -112    73  -41

julia> diag_alt(B)
3-element Array{Int8,1}:
  16
  40
 -41
```
"""
@generated function diag_alt(A::AbstractArray)
      I = Symbol[:i for _ in 1:ndims(A)]
      return :(@einsum B[i] := A[$(I...)]) # e.g., @einsum y[i] := x[i,i,i]
end

"""
```julia
diag_hyperplaneplane(A::AbstractArray)
```

Given an array of arbitrary dimesions, return the diagonal hyperplane.
"""
@generated function diag_hyperplaneplane(A::AbstractArray{T, N}) where {T, N}
      N < 3 && return :(diag_alt(A))
      I = vcat(Symbol[:i for _ in 1:(ndims(A) - 1)], :k)
      K = vcat(Symbol[:i for _ in 1:(ndims(A) - 2)], :k)
      return :(@einsum B[$(K...)] := A[$(I...)]) # e.g., @einsum y[i,k] := x[i,i,k]
end

function diagind_alt(A::AbstractArray)
      ishypercube(A) || throw(DimensionMismatch("Error trying to access the diagonal indices from a non-diagonalisable array (i.e., input array is not a hypercube)"))
      return map(i -> CartesianIndex(ntuple(d -> i, ndims(A))), axes(A, 1))
end

function diag_alt_alt(A::AbstractArray)
      return A[diagind_alt(A)]
end

#=
julia> @pretty @matmul M[i,j] := sum(k) A[i, k] * B[j, k]
begin
    local lobster = permutedims(B, (2, 1))
    local boar = A * lobster
    M = boar
end
julia> @pretty @reduce M[i,j] := sum(k) A[i, k] * B[j, k]
begin
    local lobster = orient(PermuteDims(B), (*, :, :))
    M = dropdims(sum(@__dot__(A * lobster), dims = 2), dims = 2)
end
=#


text:
endpoint_tuple = [
        (:get_oauth, :get_media_upload, "upload.json", nothing),
        (:post_oauth, :post_media_upload, "upload.json", :Media),
        (:post_oauth, :post_media_metadata_create, "metadata/create.json", nothing),
        (:post_oauth, :post_media_subtitles_delete, "subtitles/delete.json", nothing),
        (:post_oauth, :post_media_subtitles_create, "subtitles/create.json", nothing)
]

#dynamically build
for (verb, func, endp, t) in endpoint_tuple
    @eval begin
            function ($func)(;kwargs...)
                options = Dict{String, Any}()
                for arg in kwargs
                    options[string(arg[1])] = string(arg[2])
                end
                r = ($verb)($"https://upload.twitter.com/1.1/media/$endp", options)
                if r.status == 200
                    success = JSON.parse(String(r.body))
                    if ($t) != nothing && ($t) <: TwitterType
                        return ($t)(success)
                    else
                        return success
                    end
                else
                    error("Twitter API returned $(r.status) status")
                end
            end
    end
end


text:
module sparse

using GraphBLASInterface, SuiteSparseGraphBLAS

include("bfs_primitive.jl")
include("GraphBLAS_bfs.jl")
include("GraphBLAS_sparse_square_matrix.jl")
include("sparse_square_matrix.jl")

export bfs_primitive, sparse_square_matrix,
    GraphBLAS_bfs!, GraphBLAS_sparse_square_matrix!

end # module


text:
import Singular.PolynomialRing
import Singular.QQ

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

funcname(x) = PyCall.PyObject(x).func.__name__

fn_map = Dict(
              "Add" => :+,
              "Sub" => :-,
              "Mul" => :*,
              "Div" => :/,
              "Pow" => :^
)

val_map = Dict(
               "Zero"             => (0),
               "One"              => (1),
               "NegativeOne"      => (-1),
               "Half"             => QQ(1)//QQ(2),
            #    "Pi"               => :pi,
            #    "Exp1"             => :e,
            #    "Infinity"         => :Inf,
            #    "NegativeInfinity" => :(-Inf),
            #    "ComplexInfinity"  => :Inf, # error?
            #    "ImaginaryUnit"    => :im,
            #    "BooleanTrue"      => :true,
            #    "BooleanFalse"     => :false
)

# map_fn(fn_map, key)   = get(fn_map, key, key)
# map_var(var_map, key) = get(var_map, key, key)

function replace_sympy(ex, vars::Dict)
    fn = funcname(ex)

    if fn == "Symbol"
        return get(vars, string(ex), ex)
    elseif fn in ["Integer" , "Float"]
        return N(ex)
    elseif fn == "Rational"
        return QQ(convert(Int, numer(ex)))//QQ(convert(Int, denom(ex)))
    elseif haskey(val_map, fn)
        return val_map[fn]
    end
    error("should not be visited")
    # Expr(:call, map_fn(fns_map, fn), [replace_sympy(a, values=values) for a in args(ex)]...)
end

@syms x y
const ADD = (x+y).func
const MUL = (x*y).func
const POW = (x^y).func

function sym2spoly(b::Array{Sym, 1})
    fvars = union(free_symbols.(b)...)
    svars = string.(fvars)
    R, pvars = PolynomialRing(QQ, svars)
    dict = Dict(zip(string.(fvars), pvars))
    basis = [shallow_replace(p, dict) for p in b]
    basis = eval.(basis)
    return R, basis, dict
end

function shallow_replace(expr::Sym, vars::Dict)
    fn = funcname(expr)
    # println("Function: ", fn==ADD)
    a = [shallow_replace(arg, vars) for arg in expr.args]
    # println("Start: ", expr)
    # println("Start: ", a)
    if fn == "Add"
        return sum(a)
    elseif fn == "Mul"
        return prod(a)
    elseif fn == "Pow"
        return a[1]^a[2]
    end
    return replace_sympy(expr, vars)
end

# function sym2spoly(p::Sym; vars = Sym[])
#     sym2spoly([p], vars=vars)[1]
# end

function sym2spoly(poly::Sym, varmap::Dict)
    # fns_map = merge(fn_map, fns)
    # vals_map = merge(val_map, varmap)
    eval(shallow_replace(poly, varmap))
end

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

# function Ideal(basis::Array{Sym, 1}; vars = Sym[])
#     R, sbasis, varmap = sym2spoly(basis, vars)
#     return Singular.Ideal(R, sbasis...), varmap
# end

function Ideal(basis::Array{Sym, 1})
    R, sbasis, varmap = sym2spoly(basis)
    return Singular.Ideal(R, sbasis...), varmap
end

groebner(ideal::Singular.sideal) = std(ideal)

function eliminate(basis::Array{Sym, 1}, elim::Array{Sym, 1})
    ideal, varmap = Ideal(basis)
    vars = [varmap[string(v)] for v in elim]
    return Singular.eliminate(ideal, vars...)
end


text:
# phylo/branch_basics.jl
# ======================
#
# Types and methods for phylogenetic trees.
#
# Part of the Bio.Phylo module.
#
# This file contains methods for querying and setting attributes on branches,
# as well as adding and removing branches.
#
# This file is a part of BioJulia. License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

function branchlength{C}(tree::Phylogeny{C, Float64}, edge::Edge)
    return get(tree.edgedata, edge, -1.0)
end

function branchlength!{C}(tree::Phylogeny{C, Float64}, edge::Edge, value::Float64)
    tree.edgedata[edge] = value
    return tree
end

function rem_branch!(tree::Phylogeny, branch::Edge)
    rem_edge!(tree.graph, src(branch), dst(branch))
    delete!(tree.edgedata, branch)
    return tree
end

function add_branch!(tree::Phylogeny, branch::Edge, branchlength::Float64)
    add_edge!(tree.graph, branch)
    branchlength!(tree, branch, branchlength)
    return tree
end

function parent_branch(tree::Phylogeny, vertex::Int)
    @assert hasparent(tree, vertex) error("Vertex is not connected to a parent.")
    return in_edges(tree.graph, vertex)[1]
end

function child_branches(tree::Phylogeny, vertex::Int)
    return out_edges(tree.graph, vertex)
end


text:
##
## dynamcis
##

function gbm(x0::R, μ::R, σ::R, nsteps::Z=1_000, Δt::R=0.01) where R <: Real where Z <: Integer
    xs = zeros(nsteps)
    ξ  = randn(nsteps)
    x  = x0
    for i in 1:nsteps
        xs[i] = x
        Δx = x * (μ * Δt + σ * ξ[i] * sqrt(Δt))
        x += Δx
    end
    return xs
end

function gbms(x0::R, μ::R, σ::R, nsteps::Z=1_000, Δt::R=0.01, nsims::Z=100) where R <: Real where Z <: Integer
    xs = []
    for i in 1:nsims
        push!(xs, gbm(x0, μ, σ, nsteps, Δt))
    end
    return hcat(xs...)
end

# function g_est(x::Array{R,2}, tgrid::Array{R,1}, t::R, n::Z, x0::R) where R <: Real where Z <: Integer
#     xx = deepcopy(x)
#     xx ./= x0
#     t_idx = findall(t .== tgrid)
#     nsamples = Int(floor(size(x,2)/n, digits=0))
#     g = sum(
#             log( sum(xx[t_idx, (i-1)*n+1:(i*n)])/n ) / t for i in 1:nsamples
#             ) / nsamples
#     return g
# end
function g_est(x::Array{R,2}, t_idx::Z, n::Z, x0::R) where R <: Real where Z <: Integer
    xx = deepcopy(x)
    xx ./= x0
    nsamples = Int(floor(size(x,2)/n, digits=0))
    g = sum( log( sum(xx[t_idx, (i-1)*n+1:(i*n)])/n ) / t_idx for i in 1:nsamples ) / nsamples
    return g
end
function g_est(x::Array{R,2}, t::R, n::Z, x0::R) where R <: Real where Z <: Integer
    xx = deepcopy(x)
    xx ./= x0
    t_idx = findall(t .== tgrid)
    nsamples = Int(floor(size(x,2)/n, digits=0))
    g = sum( log( sum(xx[t_idx, (i-1)*n+1:(i*n)])/n ) / t for i in 1:nsamples ) / nsamples
    return g
end

function cointoss(x0::R, p::Array{R,1}, v::Array{R,1}, nsteps::Z=1_000, bet_pct::R=0.5) where R <: Real where Z <: Integer
    """
    we will see how gbm is an `attractor` for different types of multiplicative growth
    """
    c = cumsum(p)
    cc = [0; c[1:end-1]]
    xs = zeros(nsteps)
    ξ  = rand(nsteps)
    x  = x0
    for i in 1:nsteps
        xs[i] = x
        idx = first(findall((ξ[i] .< c) .* (ξ[i] .>= cc)))
        Δx = (bet_pct * x) * v[idx]
        x += Δx
    end
    return xs
end

function cointosses(x0::R, p::Array{R,1}, v::Array{R,1}, nsteps::Z=1_000, bet_pct::R=0.5, nsims::Z=100) where R <: Real where Z <: Integer
    xs = []
    for i in 1:nsims
        push!(xs, cointoss(x0, p, v, nsteps, bet_pct))
    end
    return hcat(xs...)
end

##
## simulations
##

using Pkg; Pkg.activate(".")
using PyPlot
using Random
Random.seed!(12345)

## simulate
const nsteps = 10_000
const Δt = 0.01
const nsims = 1_000
const μ = 1.0
const σ = 2.0
const x0 = 100.0
tgrid = collect(range(0.0, length=nsteps, step=Δt))
xs = gbms(x0, μ, σ, nsteps, Δt, nsims)
xs[xs .< 0] .= 1e-321

## plot
fig, ax = subplots()
ax.cla()
ax.plot(xs)
ax.set_yscale("log")

const p = [0.2, 0.8]
const v = [5.1, -1.0]
const bet_pct = 0.05
ys = cointosses(x0, p, v, nsteps, bet_pct, nsims)
ys10 = cointosses(x0, p, v, 10nsteps, bet_pct, nsims)

##
## averages
##

const nt = 10
const nn = 20
ntgrid = round.(collect(range(1.0, stop=maximum(tgrid), length=nt)), digits=2)
nngrid = Int.(round.(collect(range(1, stop=nsims, length=nn)), digits=0))
gxs = zeros(nt, nn)
gys = zeros(nt, nn)
for i in enumerate(ntgrid)
    for j in enumerate(nngrid)
        gxs[i[1],j[1]] = g_est(xs, i[2], j[2], x0)
        gys[i[1],j[1]] = g_est(ys, first(findall(i[2] .== tgrid)), j[2], x0)
    end
end

gxbar = μ - σ^2/2
gxens = μ

gybar = sum(p .*v ) - (p[1]*p[2])/bet_pct

max_stars_repo_name: jacob-roth/ergodicity

max_issues_repo_name: jacob-roth/ergodicity

max_forks_repo_name: jacob-roth/ergodicity


text:
module GeometryWrapper

using JavaCall

include("FovWrapper/FovWrapper.jl")

end


text:
using DelimitedFiles

function neighbours(mat, i, j)
	t = vec(mat[max(i-1,1):min(i+1,size(mat,1)), max(j-1,1):min(j+1,size(mat,2))])
	count(c->c=='#', t) - ((mat[i,j] == '#') ? 1 : 0)
end

function visible(mat, i, j)
	n = size(mat,1)
	m = size(mat,2)
	dirs = [[1,0],[1,1],[0,1],[-1,1],[-1,0],[-1,-1],[0,-1],[1,-1]]
	c = 0
	for dir in dirs
		pos = [i,j]
		while true
			pos += dir
			if pos[1] < 1 || pos[2] < 1 || pos[1] > n || pos[2] > m
				break
			elseif mat[pos...] == '#'
				c += 1
				break
			elseif mat[pos...] == 'L'
				break
			end
		end
	end
	c
end

function x(func, num)
	t = vec(readdlm("11.txt", String))
	n = length(t)
	m = length(t[1])
	mat = Array{Char,2}(undef, n,m)
	for i=1:n, j=1:m
		mat[i,j] = t[i][j]
	end

	while true
		mat2 = copy(mat)
		for i=1:n, j=1:m
			if mat[i,j] == 'L' && func(mat, i, j) == 0
				mat2[i,j] = '#'
			elseif mat[i,j] == '#' && func(mat, i, j) >= num
				mat2[i,j] = 'L'
			end
		end
		if mat == mat2
			break
		end
		mat=copy(mat2)
		for i=1:n, j=1:m
			if mat2[i,j] == 'L' && func(mat2, i, j) == 0
				mat[i,j] = '#'
			elseif mat2[i,j] == '#' && func(mat2, i, j) >= num
				mat[i,j] = 'L'
			end
		end
		if mat == mat2
			break
		end
	end
	count(c->c=='#', mat)
end
#part 1
println(x(neighbours, 4))
#part 2
println(x(visible, 5))


text:
module Plotly

function newPlot(graphdiv,data,layout,options)
    Pages.broadcast("script","""Plotly.newPlot("$(id)",$(data),$(layout),$(options));""")
end

end


text:
using EwaldSum
using Test

nr = 3
nk = 3

# 2D square lattice
axes = [1 0; 0 1]
alpha = 1.0
val = EwaldSum.Madelung.madelung_constant(axes, alpha, nr, nk)
@test isapprox(val, -1.9501324732857146)

# 3D cubic lattice
axes = [1 0 0; 0 1 0; 0 0 1]
alpha = 1.8
val = EwaldSum.Madelung.madelung_constant(axes, alpha, nr, nk)
@test isapprox(val, -1.4186487397403096)


text:
using PLY

fn = "../libplyfile/facemodel.ply"

ply,nelems,elem_names,file_type,version = plyopen(fn)

nelems,nprops,datatype = plygetvertexproperties(ply)

println("nelems: ", nelems, ", nprops: ", nprops, ", datatype: ", datatype)

vertices = zeros(Float32, 3, nelems)
plygetvertices(ply, nelems, datatype, vertices)
println("first 10: \n", vertices[:, 1:10])

vertexnormals = zeros(Float32, 3, nelems)
plygetvertexnormals(ply, nelems, datatype, vertexnormals)
println("first 10: \n", vertexnormals[:, 1:10])

plyclose(ply)


text:
module class
    include("./fixed_parameters.jl")
    using .fixed_parameters

    mutable struct Class
        warning::Bool
        num_students::Int64
        num_teachers::Int64

        pick::Any
        in_class::Any
        lessoning::Any

        crew::Array
        class_border::Tuple

        function Class(num_students::Int64, num_teachers::Int64)
            new(false, num_students, num_teachers, pick, in_class, lessoning)
        end
    end

    function pick(self::Class, crew::Array, start::Array)
        my_crew = [Array{Any}(undef, self.num_students), Array{Any}(undef, self.num_teachers)]
        for i in 1:2
            for j in 1:length(my_crew[i])
                my_crew[i][j] = crew[i][j+start[i]]
            end
        end
        self.crew = my_crew
    end

    function in_class(self::Class, position::Tuple)
        if position[1]>=self.class_border[1][1] && position[1]<=self.class_border[2][1] && position[2]>=self.class_border[2][2] && position[2]<=self.class_border[1][2]
            return true
        else
            return false
        end
    end

    function lessoning(self::Class, walls::Array, fakewalls::Array)
        self.class_border = ((walls[1][1][1],walls[1][2][2]), (walls[2][1][1],walls[2][2][1]))
        for i in self.crew
            for j in i
                if j.position==(nothing,nothing)
                    j.position = (rand(d:0.01:lund_size[1]-d), rand(cwh:0.01:lund_size[2]-hwh))
                end
                if !self.in_class(self, j.position) && (j.destination==nothing || !self.in_class(self, j.destination))
                    #println(self.class_border[2][1], " ", self.class_border[1][1])
                    j.destination = (rand(self.class_border[1][1]+d:0.01:self.class_border[2][1]-d), rand(self.class_border[2][2]+d:0.01:self.class_border[1][2]-d))
                    j.step = move_speed
                elseif self.in_class(self, j.position) && j.destination!=nothing && sqrt((j.destination[1]-j.position[1])^2+(j.destination[2]-j.position[2])^2)<j.step
                    j.destination = nothing
                    j.step = indoor_speed
                elseif self.in_class(self, j.position) && j.destination==nothing
                    j.step = indoor_speed
                end
                loc_1 = j.position[2]>lund_size[2]-hwh || j.position[2]<cwh
                j.move(j, fakewalls)
                loc_2 = j.position[2]>lund_size[2]-hwh || j.position[2]<cwh
                if loc_1!=loc_2
                    if !loc_2 && (self.warning || rand()<=1-j.break_rules)
                        j.radius = radius_with_masks
                    else
                        j.radius = radius_without_masks
                    end
                end
            end
        end
    end

    export Class

end


text:
module Controlz

using Polynomials
using DifferentialEquations
using PyPlot
using Interpolations
using Printf

# hipster plot theme
PyPlot.matplotlib.style.use(normpath(joinpath(pathof(Controlz), "..", "hipster.mplstyle")))

include("tf.jl")
include("systems.jl")
include("sim.jl")
include("viz.jl")
include("controls.jl")
include("show.jl")

export
    # tf.jl
    TransferFunction, zeros_poles_gain, zeros_poles_k, zero_frequency_gain, evaluate, s, proper, strictly_proper, pole_zero_cancellation, zpk_form,
    # systems.jl
    characteristic_polynomial,
    # show.jl
    # sim.jl
    simulate, interpolate,
    # viz.jl
    viz_response, viz_poles_and_zeros, nyquist_diagram, root_locus, bode_plot,
    # controlz.jl
    PController, PIController, PIDController
end