text:
test_date = Date("2000-03-31", "yyyy-mm-dd")
test_date_plus2qtr = Date("2000-09-30", "yyyy-mm-dd")
df_date = DataFrame()
df_date[!,:date] = ["2000-03-31", "2000-04-30"]
@testset "Test quarter and date utility functions" begin
    @test DSGE.prev_quarter(test_date) == Date("1999-12-31", "yyyy-mm-dd")
    @test DSGE.next_quarter(test_date) == Date("2000-06-30", "yyyy-mm-dd")
    qtr_end_ans = Vector{Date}(undef,3)
    qtr_end_ans[1] = test_date
    qtr_end_ans[2] = Date("2000-06-30", "yyyy-mm-dd")
    qtr_end_ans[3] = test_date_plus2qtr
    @test DSGE.get_quarter_ends(test_date, test_date_plus2qtr) == qtr_end_ans
    @test DSGE.quartertodate("2000-Q1") == test_date
    @test DSGE.datetoquarter(test_date) == 1
    @test DSGE.vinttodate("010331") == Date("2001-03-31", "yyyy-mm-dd")
    @test DSGE.subtract_quarters(test_date_plus2qtr, test_date) == 2
    @test DSGE.format_dates!(:date, df_date) == qtr_end_ans[1:2]
    @test DSGE.iterate_quarters(test_date, 2) == qtr_end_ans[3]
end

custom_settings = Dict{Symbol, Setting}(
    :data_vintage             => Setting(:data_vintage, "160812"),
    :cond_vintage             => Setting(:cond_vintage, "160812"),
    :cond_id                  => Setting(:cond_id, 0),
    :use_population_forecast  => Setting(:use_population_forecast, true),
    :date_forecast_start      => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :date_conditional_end     => Setting(:date_forecast_start, DSGE.quartertodate("2016-Q3")),
    :n_anticipated_shocks     => Setting(:n_anticipated_shocks, 6))

m = Model990(custom_settings = custom_settings, testing = true)
m <= Setting(:rate_expectations_source, :ois)
@testset "Test NaN and missing data checks" begin
    @test sum(isnan.(DSGE.missing2nan(missings(3)))) == 3
    df_nan2miss = DataFrame(:a => vcat(1., NaN))
    DSGE.nan2missing!(df_nan2miss)
    @test sum(ismissing.(df_nan2miss[!,:a])) == 1

    # Can we actually test? Require that FRED API key exists
    if haskey(ENV, "FRED_API_KEY") || isfile(joinpath(ENV["HOME"],".freddatarc"))
        @info "The following warnings are expected."
        df_full = load_data(m; cond_type = :full, check_empty_columns = false,
                            verbose = :none, summary_statistics = :none)
        df_semi = load_data(m; cond_type = :semi, check_empty_columns = false,
                            verbose = :none, summary_statistics = :none)

        @test_throws ErrorException load_data(m; cond_type = :full, summary_statistics = :none,
                                              verbose = :none)
        @test_throws ErrorException load_data(m; cond_type = :semi, summary_statistics = :none,
                                              verbose = :none)
    else
        @warn "Skipping portion of util test on checks for empty columns when calling load_data because FRED_API_KEY not present"
    end

    df1, df2 = DSGE.reconcile_column_names(DataFrame(a = [1, 2]), DataFrame(b = [2, 3]))
    @test sum(ismissing.(df1[!,:b])) == 2 && sum(ismissing.(df2[!,:a])) == 2

    m <= Setting(:dataroot, "")
    @test get_data_filename(m, :semi) == "data/data_cdid=00_cdvt=160812_dsid=02_vint=160812.csv"
    @test get_data_filename(m, :full) == "data/data_cdid=00_cdvt=160812_dsid=02_vint=160812.csv"
    @test get_data_filename(m, :none) == "data/data_dsid=02_vint=160812.csv"
end

nothing


text:
using DelimitedFiles
using Plots, Plots.PlotMeasures
gr()

blackscholes_f16_small  = readdlm("blackscholes_f16_small.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f32_small  = readdlm("blackscholes_f32_small.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f64_small  = readdlm("blackscholes_f64_small.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f16a_small = readdlm("blackscholes_f16a_small.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f8_small   = readdlm("blackscholes_f8_small.txt", ' ', Float64; skipstart=1)[:]

blackscholes_float_small  = readdlm("blackscholes_float_small.txt", ' ', Float64; skipstart=1)[:]
blackscholes_double_small = readdlm("blackscholes_double_small.txt", ' ', Float64; skipstart=1)[:]

blackscholes_f16_large  = readdlm("blackscholes_f16_large.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f32_large  = readdlm("blackscholes_f32_large.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f64_large  = readdlm("blackscholes_f64_large.txt", ' ', Float64; skipstart=1)[:]
blackscholes_f16a_large = readdlm("blackscholes_f16a_large.txt", ' ', Float64; skipstart=1)[:]

blackscholes_float_large  = readdlm("blackscholes_float_large.txt", ' ', Float64; skipstart=1)[:]
blackscholes_double_large = readdlm("blackscholes_double_large.txt", ' ', Float64; skipstart=1)[:]

min_small = round(minimum(blackscholes_double_small))
max_small = round(maximum(blackscholes_double_small)/5)*5

min_large = round(minimum(blackscholes_double_large))
max_large = round(maximum(blackscholes_double_large)/5)*5

function blackscholes_small(w::Int, vec::AbstractVector, sub::AbstractChar; name::String="$w-bit", fid::String="$w")
    idx = vec .> max_small
    a = plot([0, max_small], [0, max_small]; title="($sub) simsmall: $name FPnew vs 64-bit native", color=:red, legend=:none, linestyle=:dash, left_margin=80px)
    a = scatter!(a, blackscholes_double_small[.~idx], vec[.~idx]; color=:steelblue, limit=(min_small, max_small), ywiden=true, xwiden=false)
    a = scatter!(a, blackscholes_double_small[idx], repeat([max_large]; inner=sum(idx)); color=:red, m=:x)
    b = sticks(blackscholes_double_small, blackscholes_double_small-vec; linecolor=:black, fillcolor=nothing, legend=:none, title="Residual", xlimit=(min_small, max_small))

    c = plot(a, b; layout=(2,1), link=:x, size=(1200, 300))

    savefig(c, "blackscholes_small_$fid.png")
end

function blackscholes_large(w::Int, vec::AbstractVector, sub::AbstractChar; name::String="$w-bit", fid::String="$w")
    idx = vec .> max_large
    a = plot([0, max_large], [0, max_large]; title="($sub) simlarge: $name FPnew vs 64-bit native", color=:red, legend=:none, linestyle=:dash, left_margin=80px)
    a = scatter!(a, blackscholes_double_large[.~idx], vec[.~idx]; color=:steelblue, limit=(min_large, max_large), ywiden=true, xwiden=false)
    a = scatter!(a, blackscholes_double_large[idx], repeat([max_large]; inner=sum(idx)); color=:red, m=:x)
    b = sticks(blackscholes_double_large, blackscholes_double_large-vec; linecolor=:black, fillcolor=nothing, legend=:none, title="Residual", xlimit=(min_large, max_large))

    c = plot(a, b; layout=(2,1), link=:x, size=(1200, 300))

    savefig(c, "blackscholes_large_$fid.png")
end

blackscholes_small(16, blackscholes_f16_small, 'a')
blackscholes_small(32, blackscholes_f32_small, 'b')
blackscholes_small(64, blackscholes_f64_small, 'c')
blackscholes_small(16, blackscholes_f16a_small, 'd'; name = "16-bit alternate", fid = "16a")
blackscholes_small(8, blackscholes_f8_small, 'e')

blackscholes_large(16, blackscholes_f16_large, 'a')
blackscholes_large(32, blackscholes_f32_large, 'b')
blackscholes_large(64, blackscholes_f64_large, 'c')
blackscholes_large(16, blackscholes_f16a_large, 'd'; name = "16-bit alternate", fid = "16a")

f16s_idx = .~((blackscholes_f16_small .> 30) .| ((blackscholes_f16_small .== 0.) .& (blackscholes_double_small .!= 0.)) .| (blackscholes_double_small .== 0.))
err16s = (blackscholes_double_small[f16s_idx].-blackscholes_f16_small[f16s_idx])./blackscholes_double_small[f16s_idx]
bs_mape16s = 100sum(abs.(err16s))/length(err16s)
f32s_idx = .~((blackscholes_f32_small .> 30) .| ((blackscholes_f32_small .== 0.) .& (blackscholes_double_small .!= 0.)) .| (blackscholes_double_small .== 0.))
err32s = (blackscholes_double_small[f32s_idx].-blackscholes_f32_small[f32s_idx])./blackscholes_double_small[f32s_idx]
bs_mape32s = 100sum(abs.(err32s))/length(err32s)
f64s_idx = .~((blackscholes_f64_small .> 30) .| ((blackscholes_f64_small .== 0.) .& (blackscholes_double_small .!= 0.)) .| (blackscholes_double_small .== 0.))
err64s = (blackscholes_double_small[f64s_idx].-blackscholes_f64_small[f64s_idx])./blackscholes_double_small[f64s_idx]
bs_mape64s = 100sum(abs.(err64s))/length(err64s)
f16as_idx = .~((blackscholes_f16a_small .> 30) .| ((blackscholes_f16a_small .== 0.) .& (blackscholes_double_small .!= 0.)) .| (blackscholes_double_small .== 0.))
err16as = (blackscholes_double_small[f16as_idx].-blackscholes_f16a_small[f16as_idx])./blackscholes_double_small[f16as_idx]
bs_mape16as = 100sum(abs.(err16as))/length(err16as)
f8s_idx = .~((blackscholes_f8_small .> 30) .| ((blackscholes_f8_small .== 0.) .& (blackscholes_double_small .!= 0.)) .| (blackscholes_double_small .== 0.))
err8s = (blackscholes_double_small[f8s_idx].-blackscholes_f8_small[f8s_idx])./blackscholes_double_small[f8s_idx]
bs_mape8s = 100sum(abs.(err8s))/length(err8s)

f16l_idx = .~((blackscholes_f16_large .> 30) .| ((blackscholes_f16_large .== 0.) .& (blackscholes_double_large .!= 0.)) .| (blackscholes_double_large .== 0.))
err16l = (blackscholes_double_large[f16l_idx].-blackscholes_f16_large[f16l_idx])./blackscholes_double_large[f16l_idx]
bs_mape16l = 100sum(abs.(err16l))/length(err16l)
f32l_idx = .~((blackscholes_f32_large .> 30) .| ((blackscholes_f32_large .== 0.) .& (blackscholes_double_large .!= 0.)) .| (blackscholes_double_large .== 0.))
err32l = (blackscholes_double_large[f32l_idx].-blackscholes_f32_large[f32l_idx])./blackscholes_double_large[f32l_idx]
bs_mape32l = 100sum(abs.(err32l))/length(err32l)
f64l_idx = .~((blackscholes_f64_large .> 30) .| ((blackscholes_f64_large .== 0.) .& (blackscholes_double_large .!= 0.)) .| (blackscholes_double_large .== 0.))
err64l = (blackscholes_double_large[f64l_idx].-blackscholes_f64_large[f64l_idx])./blackscholes_double_large[f64l_idx]
bs_mape64l = 100sum(abs.(err64l))/length(err64l)
f16al_idx = .~((blackscholes_f16a_large .> 30) .| ((blackscholes_f16a_large .== 0.) .& (blackscholes_double_large .!= 0.)) .| (blackscholes_double_large .== 0.))
err16al = (blackscholes_double_large[f16al_idx].-blackscholes_f16a_large[f16al_idx])./blackscholes_double_large[f16al_idx]
bs_mape16al = 100sum(abs.(err16al))/length(err16al)

bs_rmsd16s = sqrt(sum((blackscholes_double_small.-blackscholes_f16_small).^2)/length(blackscholes_double_small))
bs_rmsd32s = sqrt(sum((blackscholes_double_small.-blackscholes_f32_small).^2)/length(blackscholes_double_small))
bs_rmsd64s = sqrt(sum((blackscholes_double_small.-blackscholes_f64_small).^2)/length(blackscholes_double_small))
bs_rmsd16as = sqrt(sum((blackscholes_double_small.-blackscholes_f16a_small).^2)/length(blackscholes_double_small))
bs_rmsd8s = sqrt(sum((blackscholes_double_small.-blackscholes_f8_small).^2)/length(blackscholes_double_small))

bs_rmsd16l = sqrt(sum((blackscholes_double_large.-blackscholes_f16_large).^2)/length(blackscholes_double_large))
bs_rmsd32l = sqrt(sum((blackscholes_double_large.-blackscholes_f32_large).^2)/length(blackscholes_double_large))
bs_rmsd64l = sqrt(sum((blackscholes_double_large.-blackscholes_f64_large).^2)/length(blackscholes_double_large))
bs_rmsd16al = sqrt(sum((blackscholes_double_large.-blackscholes_f16a_large).^2)/length(blackscholes_double_large))

print("16 & ", bs_mape16s, " & ", bs_rmsd16s, " & ", bs_mape16l, " & ", bs_rmsd16l, " \\\\\n")
print("32 & ", bs_mape32s, " & ", bs_rmsd32s, " & ", bs_mape32l, " & ", bs_rmsd32l, " \\\\\n")
print("64 & ", bs_mape64s, " & ", bs_rmsd64s, " & ", bs_mape64l, " & ", bs_rmsd64l, " \\\\\n")
print("16a & ", bs_mape16as, " & ", bs_rmsd16as, " & ", bs_mape16al, " & ", bs_rmsd16al, " \\\\\n")
print("8 & ", bs_mape8s, " & ", bs_rmsd8s, " & {---} & {---} \\\\\n")

swaptions_f16_tiny = readdlm("swaptions_f16_tiny.txt", ' ', Float64; skipstart=3)[:,1]
swaptions_f32_tiny = readdlm("swaptions_f32_tiny.txt", ' ', Float64; skipstart=3)[:,1]
swaptions_f64_tiny = readdlm("swaptions_f64_tiny.txt", ' ', Float64; skipstart=3)[:,1]

for (i, (f16, f32, f64)) in enumerate(zip(swaptions_f16_tiny, swaptions_f32_tiny, swaptions_f64_tiny))
    sf16 = (isnan(f16) || isinf(f16)) ? "& {\\texttt{$f16}} " : "& $f16 "
    sf32 = (isnan(f32) || isinf(f32)) ? "& {\\texttt{$f32}} " : "& $f32 "
    sf64 = (isnan(f64) || isinf(f64)) ? "& {\\texttt{$f64}} " : "& $f64 "
    print(i, ' ', sf16, sf32, sf64, "\\\\\n")
end

diff16 = filter(x -> !isnan(x)&&!isinf(x), swaptions_f64_tiny.-swaptions_f16_tiny)
mape16 = 100sum(abs.(diff16)./filter(x -> !isnan(x)&&!isinf(x), swaptions_f16_tiny))/length(diff16)
rmsd16 = sqrt(sum(diff16.^2)/length(diff16))

diff32 = filter(x -> !isnan(x)&&!isinf(x), swaptions_f64_tiny.-swaptions_f32_tiny)
mape32 = 100sum(abs.(diff32)./filter(x -> !isnan(x)&&!isinf(x), swaptions_f32_tiny))/length(diff32)
rmsd32 = sqrt(sum(diff32.^2)/length(diff32))

print("\\midrule\n")
print("MAPE & {\\SI[round-mode=places,round-precision=2]{", mape16, "}{\\percent}\\tnote{\\textdagger}} & {\\SI[round-mode=places,round-precision=2]{", mape32, "}{\\percent}} & {---} \\\\\n")
print("RMSD & {\\num[round-mode=places,round-precision=2]{", rmsd16, "}\\tnote{\\textdagger}} & {\\num[round-mode=places,round-precision=2]{", rmsd32, "}} & {---} \\\\\n")

max_swaptions = maximum(swaptions_f64_tiny)
min_swaptions = 0
a = plot([0, max_swaptions], [0, max_swaptions]; title="FPnew vs 64-bit native", color=:red, legend=:none, linestyle=:dash)
a = scatter!(a, swaptions_f64_tiny, swaptions_f16_tiny; color=:steelblue, xlimit=(min_swaptions, max_swaptions), widen=true)
b = sticks(swaptions_f64_tiny, swaptions_f64_tiny.-swaptions_f16_tiny; linecolor=:black, fillcolor=nothing, legend=:none, title="Residual", xlimit=(min_swaptions, max_swaptions), widen=true)
sw_errors = swaptions_f64_tiny[isnan.(swaptions_f16_tiny) .| isinf.(swaptions_f16_tiny)]
b = scatter!(b, sw_errors, repeat([0], inner=length(sw_errors)); m = :x, color=:red, legend=:none)

c = plot(a, b; layout=(2,1), link=:x, size=(1200, 500))

savefig(c, "swaptions_tiny.png")

ext: jl


text:
using Test
using TypeSortedCollections

using TypeSortedCollections: indices

module M
f(x::Int64) = 3 * x
f(x::Float64) = round(Int64, x / 2)

g(x::Int64, y1::Float64, y2::Int64) = x * y1 * y2
g(x::Float64, y1::Float64, y2::Int64) = x + y1 + y2
g(x::Int64, y1::Float64, y2::Float64) = x * y1 - y2
g(x::Float64, y1::Float64, y2::Float64) = x + y1 - y2
g(x::Int64, y1::Int64, y2::Float64) = x - y1 * y2

h(x::Int64) = x > 4
h(x::Float64) = x <= 4

struct Foo end
end

macro test_noalloc(expr)
    quote
        $expr
        allocs = @allocated $expr
        @test allocs == 0
    end |> esc
end

@testset "ambiguities" begin
    @test isempty(detect_ambiguities(Base, Core, TypeSortedCollections))
end

@testset "general collection interface" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    @test length(sortedx) == length(x)
    @test !isempty(sortedx)
    @test @allocated(length(sortedx)) == 0

    empty!(sortedx)
    @test length(sortedx) == 0
    @test isempty(sortedx)
end

@testset "empty input" begin
    @test isempty(TypeSortedCollection(Number[]))
    @test isempty(TypeSortedCollection(Float64[], ()))
end

@testset "map! no args" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    @test length(sortedx.data) == 2

    results = similar(x, Int64)
    map!(M.f, results, sortedx)
    allocations = @allocated map!(M.f, results, sortedx)
    @test allocations == 0
    for (index, element) in enumerate(x)
        @test results[index] == M.f(element)
    end
end

@testset "map! with args" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand(length(x))
    y2 = rand(Int, length(x))
    results = similar(x, Float64)
    map!(M.g, results, sortedx, y1, y2)
    for (index, element) in enumerate(x)
        @test results[index] == M.g(element, y1[index], y2[index])
    end
    @test_noalloc map!(M.g, results, sortedx, y1, y2)

    y2 = Number[7.; 8; 9]
    sortedy2 = TypeSortedCollection(y2)
    map!(M.g, results, sortedx, y1, sortedy2)
    for (index, element) in enumerate(x)
        @test results[index] == M.g(element, y1[index], y2[index])
    end
    @test_noalloc map!(M.g, results, sortedx, y1, sortedy2)
end

@testset "map! indices mismatch" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand(length(x))
    y2 = Number[8; 9; Float32(7)]
    sortedy2 = TypeSortedCollection(y2)
    results = similar(x, Float64)
    @test_throws ArgumentError map!(M.g, results, sortedx, y1, sortedy2)
end

@testset "map! length mismatch" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand(length(x) + 1)
    y2 = rand(length(x))
    results = similar(x, Float64)
    @test_throws DimensionMismatch map!(M.g, results, sortedx, y1, y2)
end

@testset "foreach" begin
    x = Number[4.; 5; 3.]
    sortedx = TypeSortedCollection(x)
    @test length(sortedx.data) == 2
    results = []
    foreach(sortedx) do x
        push!(results, x * 4.)
    end
    for (index, element) in enumerate(x)
        @test element * 4. in results
    end

    y1 = rand(length(x))
    y2 = Number[7.; 8; 9.]
    sortedy2 = TypeSortedCollection(y2)
    foreach(M.g, sortedx, y1, sortedy2)
    @test_noalloc foreach(M.g, sortedx, y1, sortedy2)
end

@testset "append!" begin
    x = Number[4.; 5; 3.]
    sortedx = TypeSortedCollection(x)
    @test_throws ArgumentError append!(sortedx, [Float32(6)])
    append!(sortedx, x)
    @test length(sortedx) == 2 * length(x)
end

@testset "mapreduce" begin
    x = Number[4.; 5; 3.]
    let sortedx = TypeSortedCollection(x), v0 = 2. # required to achieve zero allocations.
        result = mapreduce(M.f, +, x, init=v0)
        @test isapprox(result, mapreduce(M.f, +, sortedx, init=v0); atol = 1e-18)
        @test_noalloc mapreduce(M.f, +, sortedx, init=v0)
    end
end

@testset "any/all" begin
    x = Number[4.; 5; 3.]
    y = [missing, true, false]
    z = []
    let sortedx = TypeSortedCollection(x), sortedy = TypeSortedCollection(y), sortedz = TypeSortedCollection(z)
        @test any(M.h, x) == any(M.h, sortedx) && all(M.h, x) == all(M.h, sortedx)
        @test any(y) == any(sortedy) && all(y) == all(sortedy)
        @test any(z) == any(sortedz) && all(z) == all(sortedz)
        @test_noalloc any(M.h, sortedx)
        @test_noalloc all(M.h, sortedx)
        @test_noalloc any(sortedy)
        @test_noalloc all(sortedy)
        @test_noalloc any(sortedz)
        @test_noalloc all(sortedz)
    end
end

@testset "matching indices" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = [7.; 8.; 9.]
    sortedy1 = TypeSortedCollection(y1, indices(sortedx))
    @test length(sortedy1.data) == length(sortedx.data)
    y2 = rand(Int, length(x))
    foreach(M.g, sortedx, sortedy1, y2)
end

@testset "matching indices, empty index vector" begin
    x = [3., 4.]
    sortedx = TypeSortedCollection(x, ([1, 2], Int[]))
    results = similar(x)
    map!(identity, results, sortedx)
    @test results == x
end

@testset "preserve order" begin
    x = Number[3.; 4; 5; 6.]
    sortedx1 = TypeSortedCollection(x)
    sortedx2 = TypeSortedCollection(x, true)
    @test num_types(sortedx1) == 2
    @test num_types(sortedx2) == 3
    results = Number[]
    foreach(x -> push!(results, x), sortedx2)
    @test all(x .== results)
end

@testset "broadcast! consecutive TypeSortedCollections" begin
    # strangely, having this test set appear later in the code results in different behavior
    # see https://github.com/JuliaLang/julia/pull/23800
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y = [7.; 8.; 9.]
    sortedy = TypeSortedCollection(y, indices(sortedx))
    z = 3
    results = similar(y, Float64)
    broadcast!(M.g, results, sortedx, sortedy, z)
    @test all(results .== M.g.(x, y, z))
    @test_noalloc broadcast!(M.g, results, sortedx, sortedy, z)
end

@testset "broadcast! TSC Vec Number" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand(length(x))
    y2 = rand(Int)
    results = similar(x, Float64)
    broadcast!(M.g, results, sortedx, y1, y2)
    @test all(results .== M.g.(x, y1, y2))

    results = similar(x, Float64)
    results .= M.g.(sortedx, y1, y2)
    @test all(results .== M.g.(x, y1, y2))
    let results = results # needed on 0.7 for some reason; TODO: investigate more
        @test_noalloc broadcast!(M.g, results, sortedx, y1, y2)
    end
end

@testset "broadcast! with scalars and TSC as second arg" begin
    x = 3
    y = Number[3.; 4; 5]
    z = 5.
    sortedy = TypeSortedCollection(y)
    results = similar(y, Float64)
    results .= M.g.(x, sortedy, z)
    @test all(results .== M.g.(x, y, z))
    @test_noalloc broadcast!(M.g, results, x, sortedy, z)
end

@testset "broadcast! consecutive scalars" begin
    x = 3
    y = 4.
    z = Number[3.; 4; 5.]
    sortedz = TypeSortedCollection(z)
    results = similar(z, Float64)
    results .= M.g.(x, y, sortedz)
    @test all(results .== M.g.(x, y, z))
    @test_noalloc broadcast!(M.g, results, x, y, sortedz)
end

@testset "broadcast! Array first" begin
    x = rand(Int, 3)
    y = Number[3.; 4; 5.]
    z = rand()
    sortedy = TypeSortedCollection(y)
    results = similar(y, Float64)
    results .= M.g.(x, sortedy, z)
    @test all(results .== M.g.(x, y, z))
    @test_noalloc broadcast!(M.g, results, x, sortedy, z)
end

@testset "broadcast! indices mismatch" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand()
    y2 = Number[8; 9; Float32(7)]
    sortedy2 = TypeSortedCollection(y2)
    results = similar(x, Float64)
    @test_throws ArgumentError broadcast!(M.g, results, sortedx, y1, sortedy2)
end

@testset "broadcast! length mismatch" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand(length(x) + 1)
    y2 = rand(length(x))
    results = similar(x, Float64)
    @test_throws DimensionMismatch results .= M.g.(sortedx, y1, y2)

    y1 = rand()
    y2 = rand(length(x) + 1)
    @test_throws DimensionMismatch results .= M.g.(sortedx, y1, y2)

    results = rand(length(x) + 1)
    y1 = rand()
    y2 = rand(Int)
    @test_throws DimensionMismatch results .= M.g.(sortedx, y1, y2)
end

@testset "broadcast! matching indices" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    y1 = rand()
    y2 = [7.; 8.; 9.]
    sortedy2 = TypeSortedCollection(y2, indices(sortedx))
    results = similar(x, Float64)
    broadcast!(M.g, results, sortedx, y1, sortedy2)
    @test all(results .== M.g.(x, y1, y2))
    @test_noalloc broadcast!(M.g, results, sortedx, y1, sortedy2)
end

@testset "broadcast! TSC destination" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    results = typeof(sortedx)(indices(sortedx))
    results .= 3 .* sortedx
    @test results.data[1][1] === 3 * 3.
    @test results.data[2][1] === 3 * 4
    @test results.data[2][2] === 3 * 5
end

@testset "eltype" begin
    x = [4.; 5; 3.; Int32(2); Int16(1); "foo"]
    let sortedx = TypeSortedCollection(x)
        @test eltype(sortedx) == Union{Float64, Int64, Int32, Int16, String}
        @test_noalloc eltype(sortedx)
    end
end

@testset "push!" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    @test length(sortedx) == 3
    @test_throws ArgumentError push!(sortedx, "foo")
    push!(sortedx, 8)
    @test length(sortedx) == 4
    @test mapreduce(x -> x == 8, (a, b) -> a || b, sortedx, init=false)
    results = similar(x, length(sortedx))
    map!(identity, results, sortedx)
    @test last(results) == 8
    @test @inferred(push!(sortedx, 1)) isa typeof(sortedx)
end

@testset "eltypes/vectortypes" begin
    x = Number[3.; 4; 5]
    sortedx = TypeSortedCollection(x)
    let T = typeof(sortedx)
        @test @inferred eltypes(T) == Tuple{Float64, Int}
        elt = eltypes(T)
        @test @inferred vectortypes(elt) == Tuple{Vector{Float64}, Vector{Int}}
        @test typeof(sortedx.data) == vectortypes(elt)
    end
end

@testset "broadcast with user types" begin
    x = fill(M.Foo(), 3)
    sortedx = TypeSortedCollection(x)
    sortedx .= Ref(M.Foo())
end


text:
struct blob
    k::Int64
    n_below::Int64
    cost::Float64
    neighbor::Set{Int64}
    bin_BB_Tree::bin_BB_Tree
    function blob(k::Int64, n_below::Int64, neighbor_in, bb_tree::bin_BB_Tree, scale::Float64)
        bbb = bb_tree.box
        bbb = OBB(bbb, bbb)
        cost = blobCost(bbb, n_below, scale)
        return new(k, n_below, cost, Set{Int64}(neighbor_in), bb_tree)
    end
end

function createBlobDictionary(eM::eMesh, vec_obb_leaf::Vector{OBB}, vec_tri_tet, scale::Float64)
    vec_neighbor, is_abort = extractTriTetNeighborInformation(vec_tri_tet)
    dict_blob = Dict{Int64,blob}()
    if is_abort
        return dict_blob, true
    end
    for (k, ind_k) = enumerate(vec_tri_tet)
        aabb_k = vec_obb_leaf[k]
        tree_k = bin_BB_Tree(k, aabb_k)
        dict_blob[k] = blob(k, 1, vec_neighbor[k], tree_k, scale)
    end
    return dict_blob, false
end

function extractTriTetNeighborInformation(vec_tri_tet::Vector{SVector{N,Int64}}) where {N}
    dict_vert_pair, is_abort = createSharedEdgeFaceDict(vec_tri_tet)
    vec_neighbor = [Vector{Int64}() for _ = vec_tri_tet]
    if is_abort
        return vec_neighbor, true
    end
    for key_k = keys(dict_vert_pair)
        i_pair = dict_vert_pair[key_k]
        if (N == 3) || all(i_pair .!= -9999)  # tets will not have a partner if on the boundary
            push!(vec_neighbor[i_pair[1]], i_pair[2])
            push!(vec_neighbor[i_pair[2]], i_pair[1])
        end
    end
    return vec_neighbor, false
end

function createSharedEdgeFaceDict(vec_tri_tet::Vector{SVector{N,Int64}}) where {N}
    dict_vert_pair = Dict{SVector{N-1, Int64}, MVector{2,Int64}}()
    for (k_tet_tri, i_vertices) = enumerate(vec_tri_tet)
        for j = 1:N
            sv_new = sortEdgeFace(i_vertices, j)
            if haskey(dict_vert_pair, sv_new)
                mv_existing = dict_vert_pair[sv_new]
                if mv_existing[2] != -9999
                    (N == 3) && error("three triangles share the same edge something is wrong")
                    (N == 4) && error("three tetrahedrons share the same face something is wrong")
                end
                mv_existing[2] = k_tet_tri
            else
                dict_vert_pair[sv_new] = MVector{2,Int64}(k_tet_tri, -9999)
            end
        end
    end
    if N == 3  # confirm that every edge has a partner
        for (key_k, val_k) = dict_vert_pair
            if any(val_k .== -9999)
                #  && error("edge $key_k is connected to a single triangle $(val_k[1])")
                return dict_vert_pair, true
            end
        end
    end
    return dict_vert_pair, false
end

function blobCost(aabb::OBB, n_below::Int64, scale::Float64)
    cA = 1.0
    cV = 1.0
    V = 0.0
    V += n_below * log2(2 * n_below)
    V += cA * area(aabb) / (scale^2)
    V += cV * volume(aabb) / (scale^3)
    return V  # PriorityQueue returns lowest first
end

function doCombineBlob(dict_blob, a::blob, b::blob, k_next::Int64, scale::Float64)
    delete!(a.neighbor, b.k)
    delete!(b.neighbor, a.k)
    c_neighbor = union(a.neighbor, b.neighbor)  # add to neighbor_c
    tree_c = bin_BB_Tree(a.bin_BB_Tree, b.bin_BB_Tree)
    n_below_c = a.n_below + b.n_below
    blob_c = dict_blob[k_next] = blob(k_next, n_below_c, c_neighbor, tree_c, scale)
    return k_next + 1, blob_c
end

function calcMarginalCost(a::blob, b::blob, scale::Float64)
    c_cost = blobCost(OBB(a.bin_BB_Tree.box, b.bin_BB_Tree.box), a.n_below + b.n_below, scale)
    return c_cost - a.cost - b.cost
end

function refreshCostQueue!(dict_blob, pq_delta_cost, blob_a::blob, blob_c::blob, scale::Float64)
    a_k = blob_a.k
    for b_k = blob_a.neighbor  # for each neighbor in blob_a
        delete!(pq_delta_cost, minmax(b_k, a_k))  # delete old cost
        pq_delta_cost[minmax(b_k, blob_c.k)] = calcMarginalCost(dict_blob[b_k], blob_c, scale)  # add new cost
        blob_k = dict_blob[b_k]
        (a_k in blob_k.neighbor) || error("something is wrong")
        delete!(blob_k.neighbor, a_k)
        push!(blob_k.neighbor, blob_c.k)
    end
    delete!(dict_blob, a_k)
end

function createBlobPriorityQueue(dict_blob, scale)
    pq_delta_cost = PriorityQueue{Tuple{Int64,Int64},Float64}()
    for k_a = keys(dict_blob)
        blob_a = dict_blob[k_a]
        for k_b = blob_a.neighbor
            (k_a < k_b) && (pq_delta_cost[(k_a, k_b)] = calcMarginalCost(blob_a, dict_blob[k_b], scale))
        end
    end
    return pq_delta_cost
end

function bottomUp!(dict_blob, pq_delta_cost, scale)
    k_next = length(dict_blob) + 1
    while !isempty(pq_delta_cost)
        key_lowest, val_lowest = dequeue_pair!(pq_delta_cost)  # automatically removes entry
        a_k, b_k = key_lowest
        blob_a = dict_blob[a_k]
        blob_b = dict_blob[b_k]
        k_next, blob_c = doCombineBlob(dict_blob, blob_a, blob_b, k_next, scale)
        refreshCostQueue!(dict_blob, pq_delta_cost, blob_a, blob_c, scale)  # deal with blob_a
        refreshCostQueue!(dict_blob, pq_delta_cost, blob_b, blob_c, scale)  # deal with blob_b
    end
    return nothing
end

function eMesh_to_tree(eM::eMesh{T1,T2}) where {T1,T2}
    isa(eM, eMesh{Tri,Tet}) && error("Cannot create tree for eMesh{Tri,Tet} use as_tri_eMesh or as_tet_eMesh on input first.")
    n_leaf = ifelse(T1==Tri, n_tri, n_tet)(eM)
    if (n_leaf == 1)
        if T1 == Tri
            obb = calc_obb(eM.point[eM.tri[1]])
        else
            obb = calc_obb(eM.point[eM.tet[1]])
        end
        return bin_BB_Tree(1, obb)
    end
    point = eM.point
    vec_tri_tet = ifelse(isa(eM, eMesh{Tri,Nothing}), eM.tri, eM.tet)
    scale = sum(calc_obb(point).e) / 3
    if T1 == Tri
        vec_obb_leaf = [calc_obb(eM.point[k]) for k = eM.tri]
    else
        vec_obb_leaf = [calc_obb(eM.point[k]) for k = eM.tet]
    end
    dict_blob, is_abort = createBlobDictionary(eM, vec_obb_leaf, vec_tri_tet, scale)
    is_abort && error("not implemented error: disconnected mesh")

    # if is_abort  # non-connected meshes cannot be created "bottom up"
        # return top_down(point, vec_tri_tet)
    # end

    pq_delta_cost = createBlobPriorityQueue(dict_blob, scale)
    bottomUp!(dict_blob, pq_delta_cost, scale)
    all_tree = collect(values(dict_blob))
    if (length(all_tree) != 1)
        error("not implemented error: disconnected mesh")
    else

        tree = all_tree[1].bin_BB_Tree

        tight_fit_leaves!(eM, tree)

        return tree

        # return all_tree[1].bin_BB_Tree
    end
end

function tight_fit_leaf!(eM::eMesh{Tri,Nothing}, eT::bin_BB_Tree{OBB})
    eT.box = fit_tri_obb(vertex_pos_for_tri_ind(eM, eT.id))
end

function tight_fit_leaf!(eM::eMesh{Nothing,Tet}, eT::bin_BB_Tree{OBB})
    eT.box = fit_tet_obb(vertex_pos_for_tet_ind(eM, eT.id), eM.ϵ[eM.tet[eT.id]])
end

function tight_fit_leaves!(eM::eMesh, eT::bin_BB_Tree{OBB})
    if is_leaf(eT)
        tight_fit_leaf!(eM, eT)
    else
        tight_fit_leaves!(eM, eT.node_1)
        tight_fit_leaves!(eM, eT.node_2)
    end
end


text:
using FastGaussQuadrature

abstract type AbstractQuadrature <: JalerkinType end;

#--- Quadrature

struct Quadrature{N} <: AbstractQuadrature
    """
    Quadrature rule with N nodes and weights of type T
    """
    nodes :: SVector{N, fp_precision()}
    weights :: SVector{N, fp_precision()}
end

@inline function get_nodes(q::Quadrature)
    q.nodes
end

@inline function get_weights(q::Quadrature)
    q.weights
end

@inline function quad(q::AbstractQuadrature, values::AbstractArray)
    """
    Compute integral of a function using a concrete quadrature rule,
    where 'values' are the values of the function on q.nodes
    """
    get_weights(q) ⋅ values
end

@inline function quad(q::AbstractQuadrature, f)
    """
    Compute integral of a function using a concrete quadrature formula
    """
    quad(q, f.(get_nodes(q)))
end

#--- Constructors for concrete quadrature rules

"""Gaussian quadrature rule on [-1,1] with n points"""
function GaussianQuadrature(n)
    nodes, weights = gausslegendre(n)
    Quadrature{n}(nodes, weights)
end

"""
Trapezoid quadrature rule on interval [-1, 1]
"""
function TrapezoidQuadrature()
    nodes = [-1, 1]
    weights = [1, 1]
    Quadrature{2}(nodes, weights)
end

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
function freezing_time(x, θ)
    if θ*x >= 0 # sic!
        return Inf
    else
        return -x/θ
    end
end

function queue_time!(Q, t, x, θ, i, b, f, κ)
    trefl = poisson_time(b[i], rand())
    tfreeze = freezing_time(x[i], θ[i])
    if tfreeze <= trefl
        f[i] = true
        Q[i] = t[i] + tfreeze
    else
        f[i] = false
        Q[i] = t[i] + trefl
    end
    return
end

function sspdmp_inner!(Ξ, G, G2, ∇ϕ, t, x, θ, Q, c, b, t_old, f, θf, κ, (acc, num),
     F::Union{ZigZag,FactBoomerang}, args...; factor=1.5, adapt=false)
    n = length(x)
    while true
        ii, t′ = peek(Q)
        refresh = ii > n
        i = ii - refresh*n
        refresh && error("refreshment not implemented")
        t, x, θ = smove_forward!(G, i, t, x, θ, t′, F)
        if f[i] # to be frozen
            if abs(x[i]) > 1e-8
                error("x[i] = $(x[i]) !≈ 0")
            end
            x[i] = -0*θ[i]
            θf[i], θ[i] = θ[i], 0.0 # stop and save speed
            t_old[i] = t[i]
            f[i] = false
            Q[i] = t[i] - log(rand())/κ
            t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
            for j in neighbours(G, i)
                if θ[j] != 0 # only non-frozen, especially not i
                    b[j] = ab(G, j, x, θ, c, F)
                    t_old[j] = t[j]
                    queue_time!(Q, t, x, θ, j, b, f, κ)
                end
            end
        elseif x[i] == 0 && θ[i] == 0 # was frozen
            θ[i], θf[i] = θf[i], 0.0 # unfreeze, restore speed
            t_old[i] = t[i]
            t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
            for j in neighbours(G, i)
                if θ[j] != 0 # only non-frozen
                    b[j] = ab(G, j, x, θ, c, F)
                    t_old[j] = t[j]
                    queue_time!(Q, t, x, θ, j, b, f, κ)
                end
            end
        else
            l, lb = sλ(∇ϕ, i, t, x, θ, t′, F, args...), sλ̄(b[i], t[i] - t_old[i])
            num += 1
            if rand()*lb < l
                acc += 1
                if l >= lb
                    !adapt && error("Tuning parameter `c` too small.")
                    acc = num = 0
                    adapt!(c, i, factor)
                end
                t, x, θ = smove_forward!(G2, i, t, x, θ, t′, F)
                θ = reflect!(i, x, θ, F)
                for j in neighbours(G, i)
                    if θ[j] != 0
                        b[j] = ab(G, j, x, θ, c, F)
                        t_old[j] = t[j]
                        queue_time!(Q, t, x, θ, j, b, f, κ)
                    end
                end
            else
                b[i] = ab(G, i, x, θ, c, F)
                t_old[i] = t[i]
                queue_time!(Q, t, x, θ, i, b, f, κ)
                continue
            end
        end
        push!(Ξ, event(i, t, x, θ, F))
        return t, x, θ, t′, (acc, num), c, b, t_old
    end
end

function sspdmp(∇ϕ, t0, x0, θ0, T, c, κ, F::Union{ZigZag,FactBoomerang}, args...;
        factor=1.8, adapt=false)
    n = length(x0)
    t′ = t0
    t = fill(t′, size(θ0)...)
    t_old = copy(t)
    f = zeros(Bool, n)
    G = [i => rowvals(F.Γ)[nzrange(F.Γ, i)] for i in eachindex(θ0)]
    G2 = [i => setdiff(union((G[j].second for j in G[i].second)...), G[i].second) for i in eachindex(G)]
    x, θ = copy(x0), copy(θ0)
    θf = zero(θ)
    num = acc = 0
    Q = SPriorityQueue{Int,Float64}()
    b = [ab(G, i, x, θ, c, F) for i in eachindex(θ)]
    for i in eachindex(θ)
        trefl = poisson_time(b[i], rand())
        tfreez = freezing_time(x[i], θ[i])
        if trefl > tfreez
            f[i] = true
            enqueue!(Q, i => tfreez)
        else
            f[i] = false
            enqueue!(Q, i => trefl)
        end
    end
    if hasrefresh(F)
        for i in eachindex(θ)
            enqueue!(Q, (n + i) => waiting_time_ref(F))
        end
    end
    Ξ = Trace(t0, x0, θ0, F)
    while t′ < T
        t, x, θ, t′, (acc, num), c,  b, t_old = sspdmp_inner!(Ξ, G, G2, ∇ϕ, t, x, θ, Q,
                    c, b, t_old, f, θf, κ, (acc, num), F, args...; factor=factor, adapt=adapt)
    end
    #t, x, θ = ssmove_forward!(t, x, θ, T, F)
    Ξ, (t, x, θ), (acc, num), c
end


text:
using NNlib: batched_mul

"""
GeometricFlux isn't written to sustain many parallel graph computations, and for each graph
it relies on a sequential implementation that isn't supported by GPUs.
Thus we propose our own version of the GeometricFlux API, using only matrices as input, eventually in 3D, to
make it possible to compute GNNs on many graphs simultaneously.
"""

function (g::GeometricFlux.GraphConv)(fgs::BatchedFeaturedGraph{Float32})
    A, X = fgs.graph, fgs.nf
    return BatchedFeaturedGraph{Float32}(
        fgs.graph;
        nf = g.σ.(g.weight1 ⊠ X .+ g.weight2 ⊠ X ⊠ A .+ g.bias),
        ef = fgs.ef,
        gf = fgs.gf
    )
end

struct GraphNorm <: GeometricFlux.GraphNet
    batchNormLayer::Flux.BatchNorm

    GraphNorm(channels::Integer, σ = identity; initβ = zeros, initγ = ones, ϵ = Float32(1e-8), momentum = Float32(.1)) = new(Flux.BatchNorm(channels, σ; initβ = initβ, initγ = initγ, ϵ = ϵ, momentum = momentum))
    GraphNorm(bn::Flux.BatchNorm) = new(bn)
end

Flux.@functor GraphNorm
function (gn::GraphNorm)(fgs::BatchedFeaturedGraph{Float32})
    X = fgs.nf
    Y = reshape(X, size(X)[1:end - 2]..., :)
    return BatchedFeaturedGraph{Float32}(
        fgs.graph;
        nf = reshape(gn.batchNormLayer(Y), size(X)),
        ef = fgs.ef,
        gf = fgs.gf,
    )
end


text:
## Discrete Time
function DiffEqBase.solve(prob::AbstractDiscreteProb{dType, true}, alg::AbstractKoopmanAlgorithm;
    B::AbstractArray = [], digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X[:,1:end-1]
    DX = prob.X[:,2:end]

    # Create a basis
    @variables x[1:size(X, 1)] t
    b = Basis(x, x,  iv = t)

    C = diagm(ones(dType, size(X,1)))

    k = alg(X, DX)
    B = Matrix{dType}(undef, 0, 0)
    # Updateable for all measurements
    Q = DX*X'
    P = DX*X'

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, BitVector((true for i in 1:size(X,1))), b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstractDiscreteProb{dType, false}, alg::AbstractKoopmanAlgorithm;
    B::AbstractArray = [], digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X[:,1:end-1]
    DX = prob.X[:,2:end]
    U = prob.U[:, 1:end-1]

    # Create a basis
    @variables x[1:size(X, 1)] u[1:size(U, 1)] t
    b = Basis([x; u], x, controls = u, iv = t)

    inds = BitVector([ones(Bool, size(X,1)); zeros(Bool, size(U,1))])
    C = diagm(ones(dType, size(X,1)))

    # The input maps

    if isempty(B)
        k, B = alg(X, DX, U)
    else
        k, B = alg(X, DX, U, B)
    end
    # Updateable for all measurements
    Q = DX*[X;U]'
    P = DX*[X;U]'

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstractDiscreteProb{dType, true}, b::Basis, alg::AbstractKoopmanAlgorithm;
    digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X[:,1:end-1]
    DX = prob.X[:,2:end]
    p = prob.p
    t = prob.t

    Ψ₀ = b(X, p, t[1:end-1])
    Ψ₁ = b(DX, p, t[2:end])

    k = alg(Ψ₀, Ψ₁)

    Q = Ψ₁*Ψ₀'
    P = Ψ₀*Ψ₀'
    B = zeros(dType, 0, 0)

    # Outpumap -> just the state dependent
    C = DX / Ψ₁

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, BitVector((true for i in 1:size(Ψ₀,1))), b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstractDiscreteProb{dType, false}, b::Basis, alg::AbstractKoopmanAlgorithm;
    digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X[:,1:end-1]
    DX = prob.X[:,2:end]
    U = prob.U[:, 1:end-1]
    p = prob.p
    t = prob.t

    Ψ₀ = b(X, p, t[1:end-1], U)
    Ψ₁ = b(DX, p,t[2:end], U)

    # Find the indexes of the control states
    inds = .! is_dependent(map(eq->Num(eq.rhs),equations(b)), Num.(controls(b)))[1,:]

    k, B = alg(Ψ₀[inds, :], Ψ₁[inds, :], Ψ₀[.!inds, :])

    Q = Ψ₁[inds, :]*Ψ₀'
    P = Ψ₀*Ψ₀'

    # Outpumap -> just the state dependent
    C = DX / Ψ₁[inds,:]

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
end

## Continouos Time
function DiffEqBase.solve(prob::AbstracContProb{dType, true}, alg::AbstractKoopmanAlgorithm;
    B::AbstractArray = [], digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X
    DX = prob.DX

    # Create a basis
    @variables x[1:size(X, 1)] t
    b = Basis(x, x,  iv = t)

    C = diagm(ones(dType, size(X,1)))

    k = alg(X, DX)
    B = Matrix{dType}(undef, 0, 0)
    # Updateable for all measurements
    Q = DX*X'
    P = DX*X'

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, BitVector((true for i in 1:size(X,1))), b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstracContProb{dType, false}, alg::AbstractKoopmanAlgorithm;
    B::AbstractArray = [], digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X
    DX = prob.DX
    U = prob.U

    # Create a basis
    @variables x[1:size(X, 1)] u[1:size(U, 1)] t
    b = Basis([x; u], x, controls = u, iv = t)

    inds = BitVector([ones(Bool, size(X,1)); zeros(Bool, size(U,1))])
    C = diagm(ones(dType, size(X,1)))

    # The input maps

    if isempty(B)
        k, B = alg(X, DX, U)
    else
        k, B = alg(X, DX, U, B)
    end
    # Updateable for all measurements
    Q = DX*[X;U]'
    P = DX*[X;U]'

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstracContProb{dType, true}, b::Basis, alg::AbstractKoopmanAlgorithm;
    digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X
    DX = prob.DX
    p = prob.p
    t = prob.t

    Ψ₀ = b(X, p, t)
    Ψ₁ = similar(Ψ₀)

    J = jacobian(b)

    for i in 1:size(DX, 2)
        Ψ₁[:, i] .= J(X[:, i], p, t[i])*DX[:, i]
    end

    k = alg(Ψ₀, Ψ₁)

    Q = Ψ₁*Ψ₀'
    P = Ψ₀*Ψ₀'
    B = zeros(dType, 0, 0)

    # Outpumap -> just the state dependent
    C = prob.DX / Ψ₁

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, BitVector((true for i in 1:size(Ψ₀,1))), b, alg, digits = digits, eval_expression = eval_expression)
end

function DiffEqBase.solve(prob::AbstracContProb{dType, false}, b::Basis, alg::AbstractKoopmanAlgorithm;
    digits::Int = 10, operator_only::Bool = false,
    eval_expression = false,
    kwargs...) where {dType <: Number}
    # Check the validity
    @assert is_valid(prob) "The problem seems to be ill-defined. Please check the problem definition."

    X = prob.X
    DX = prob.DX
    U = prob.U
    p = prob.p
    t = prob.t

    Ψ₀ = b(X, p, t, U)
    Ψ₁ = similar(Ψ₀)

    J = jacobian(b)

    for i in 1:size(DX, 2)
        Ψ₁[:, i] .= J(X[:, i], p, t[i], U[:,i])*DX[:, i]
    end

    # Find the indexes of the control states
    inds = .! is_dependent(map(eq->Num(eq.rhs),equations(b)), Num.(controls(b)))[1,:]

    k, B = alg(Ψ₀[inds, :], Ψ₁[inds, :], Ψ₀[.!inds, :])

    Q = Ψ₁[inds, :]*Ψ₀'
    P = Ψ₀*Ψ₀'

    # Outpumap -> just the state dependent
    C = prob.DX / Ψ₁[inds,:]

    operator_only && return (K = k, C = C, B = B, Q = Q, P = P)

    return DataDrivenSolution(prob, k, C, B, Q, P, inds, b, alg, digits = digits, eval_expression = eval_expression)
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains tests for SubGraph operations.

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

@testset "Subgraph" begin
   gtype = completegraph(10)

   g = Graph(completeindxs(10))
   eit = edges(g)

   # Set Labels
   ls = map(string, 1:10)
   setlabel!(g, ls)

   # Set vertex properties
   setvprop!(g, :, 1:10, :p1)
   setvprop!(g, :, 1:10, :p2)
   setvprop!(g, :, 1:10, :p3)
   setvprop!(g, :, 1:10, :p4)

   # Set edge properties
   seteprop!(g, :, 1:90, :p1)
   seteprop!(g, :, 1:90, :p2)
   seteprop!(g, :, 1:90, :p3)
   seteprop!(g, :, 1:90, :p4)

   # Copy
   g1 = subgraph(g)
   @test isequal(g1, g)

   # VS
   g2 = subgraph(g, 1:5)
   @test nv(g2) == 5
   @test ne(g2) == 20
   @test size(vdata(g2)) == (5, 4)
   @test size(edata(g2)) == (20, 4)
   @test encode(g2) == ls[1:5]

   # VPROPS
   g3 = subgraph(g, :, [:p1, :p3])
   @test nv(g3) == 10
   @test names(vdata(g3)) == [:p1, :p3]
   @test size(vdata(g3)) == (10, 2)

   # VS & VPROPS
   g4 = subgraph(g, 1:7, [:p2, :p4])
   @test nv(g4) == 7
   @test ne(g4) == 42
   @test names(vdata(g4)) == [:p2, :p4]
   @test size(vdata(g4)) == (7,2)
   @test size(edata(g4)) == (42,4)
   @test encode(g4) == ls[1:7]

   # ES
   es = eit[1:45]
   g5 = subgraph(g, es)
   @test nv(g5) == 10
   @test ne(g5) == 45
   @test size(edata(g5)) == (45,4)

   # EPROPS
   g6 = subgraph(g, :, :, [:p1,:p4])
   @test ne(g6) == 90
   @test names(edata(g6)) == [:p1,:p4]
   @test size(edata(g6)) == (90,2)

   # ES & EPROPS
   es = eit[35:70]
   g7 = subgraph(g, es, [:p1])
   @test nv(g7) == 10
   @test ne(g7) == 36
   @test names(edata(g7)) == [:p1]
   @test size(edata(g7)) == (36,1)

   # VS & ES
   es = [1=>2, 2=>3, 3=>1]
   g8 = subgraph(g, 1:3, es)
   @test nv(g8) == 3
   @test ne(g8) == 3
   @test encode(g8) == ls[1:3]

   es = [1=>2, 2=>3, 3=>4, 4=>5, 5=>1]
   g9 = subgraph(g, 1:5, es, [:p1], [:p2])
   @test nv(g9) == 5
   @test ne(g9) == 5
   @test names(vdata(g9)) == [:p1]
   @test size(vdata(g9)) == (5,1)
   @test names(edata(g9)) == [:p2]
   @test size(edata(g9)) == (5,1)
   @test encode(g9) == ls[1:5]
end


text:
# test parsing
@test Traits.get_fname(:(sin(x::T1) = sin(x.t1)))==:sin
@test Traits.get_fname(:(sin{X<:Int}(x::X) = sin(x.t1)))==:(sin{X<:Int})
@test Traits.get_fname_only(:(sin{X<:Int}(x::X) = sin(x.t1)))==:sin

# Testing @traitimpl
@traitdef Tr100{X,Y} begin
    fun1(X,Y)
    fun2(X)
    fun3(Y)
end
type A100
end
type B100
end
# manual implementation:
fun1(a::A100, b::B100) = 5
fun2(a::A100) = 6
fun3(b::B100) = 7

head = :(Tr100{A100,B100})

body = quote
    fun1(a::A100, b::B100) = 5
    fun2(a::A100) = 6
    fun3(b::B100) = 7
end
trait = Tr100{A100,B100}

sigs = Any[
        Any[:(a::A100), :(b::B100)],
        Any[:(a::A100)],
        Any[:(b::B100)]
        ]
i = 1
for ln in enumerate(body.args)
    fn = ln[2]
    if isa(fn, Expr) && !(fn.head==:line)
        @test Traits.get_fsig(fn)==sigs[i]
        i+=1
    end
end

implfs = Traits.parse_body(body)
#@test Traits.check_macro_body(body.args, implfs, trait)

# check @traitimpl
type A110
end
type B110
end
#  implementation
@traitimpl Tr100{A110,B110} begin
    fun1(a::A110, b::B110) = 5
    fun2(a::A110) = 6
    fun3(b::B110) = 7
end

### issue 2
@traitdef Tr101{X} begin
    Y = getTr101(X)
    getTr101(Type{X}) -> DataType
    fun11(X, Vector{Y})
end

@traitimpl Tr101{Int} begin
    getTr101(::Type{Int}) = Integer
    fun11{Y<:Integer}(x::Int, y::Vector{Y}) = y[x]
end


text:
######################################################################
# includes stuff I don't want to need distributions for...
######################################################################
export traj2mat

function my_pdf(x::NTuple{2,Float64},m::Vector{Float64},S::Matrix{Float64})
	dx = x[1] - m[1]
	dy = x[2] - m[2]

	Si = inv(S)
	e_stuff = dx*Si[1,1]*dx + dx*Si[1,2]*dy + dy*Si[2,1]*dx + dy*Si[2,2]*dy
	return det(sqrtm(Si)) * exp(-0.5 * e_stuff) / (2.0 * pi)
end

function my_pdf(x::VF, m::VF, S::MF)
	k = length(x)
	xm = x - m
	return exp(-0.5 * dot(xm, inv(S)*xm)) / sqrt((2pi)^k * det(S))
end

# In order to be faster, precompute inverse and denominator.
# That way, I don't have to recompute them for every cell in domain.
function fast_pdf(x::VF, m::VF, inv_S::MF, den::Float64)
	xm = x - m
	return exp(-0.5 * dot(xm, inv_S*xm)) / den
end

function centroid(d::Matrix{Float64}, L::Float64)
	x_val = 0.0; y_val = 0.0
	n = size(d,1)
	cell_size = L / n
	d_sum = 0.0
	for x = 1:n
		for y = 1:n
			x_val += (x-.5) * d[x,y]
			y_val += (y-.5) * d[x,y]
			d_sum += d[x,y]
		end
	end
	return x_val*cell_size/d_sum, y_val*cell_size/d_sum
end

function covariance(d::Matrix{Float64}, L::Float64)
	mu_x, mu_y = centroid(d, L)
	c_xx = c_xy = c_yy = 0.0
	n = size(d,1)
	cell_size = L / n
	d_sum = 0.0
	for xi = 1:n
		for yi = 1:n
			x = (xi-0.5)*cell_size
			y = (yi-0.5)*cell_size

			c_xx += d[xi,yi] * x * x
			c_yy += d[xi,yi] * y * y
			c_xy += d[xi,yi] * (x - mu_x) * (y - mu_y)
			d_sum += d[xi,yi]
		end
	end
	c_xx = c_xx/d_sum - (mu_x * mu_x)
	c_yy = c_yy/d_sum - (mu_y * mu_y)
	c_xy = c_xy/d_sum

	# add 1e-3 to ensure we are positive definite
	return [c_xx+1e-3 c_xy; c_xy c_yy+1e-3]
end

# unscaled
# g_f1 = grad_1 f
# g_f2 = grad_2 f
# u1, u2 are components of direction
# basically,
# [g_f1, g_f2]' * [u1, u2]'
function directional_derivative(g_f1::VVF, g_f2::VVF, u1::VVF, u2::VVF)
	N = length(g_f2)
	dd = 0.0
	for i = 1:N
		dd += dot(g_f1[i], u1[i]) + dot(g_f2[i], u2[i])
	end
	dd += dot(g_f1[N+1], u1[n+1])
	return dd
end
function directional_derivative(g_f1::Matrix{Float64}, g_f2::Matrix{Float64}, u1::VVF, u2::VVF)
	#N = length(g_f1)
	N = size(g_f2, 2)
	#n = size(u1[1])
	#m = size(u2[1])
	#println("N = ", N)
	dd = 0.0
	for i = 1:N
		# old way
		#dd += g_f1[1,i]*u1[i][1] + g_f1[2,i]*u1[i][2]
		#dd += g_f2[1,i]*u2[i][1] + g_f2[2,i]*u2[i][2]

		dd += dot(g_f1[:,i], u1[i]) + dot(g_f2[:,i], u2[i])
	end
	dd += dot(g_f1[:,N+1], u1[N+1])
	return dd
end
# TODO: turn this one into somepin correct for matrix{Float64} version
function directional_derivative(g_f1::Matrix{Float64}, g_f2::Matrix{Float64}, u1::Matrix{Float64}, u2::Matrix{Float64})
	#N = length(g_f1)
	N = size(g_f1, 2)
	dd = 0.0
	for i = 1:N
		#dd += dot(g_f1[i], u1[i]) + dot(g_f2[i], u2[i])
		dd += g_f1[1,i]*u1[1,i] + g_f1[2,i]*u1[2,i]
		dd += g_f2[1,i]*u2[1,i] + g_f2[2,i]*u2[2,i]
	end
	return dd
end
# like above, but scales the direction first
function scaled_dd(g_f1::VVF, g_f2::VVF, u1::VVF, u2::VVF)
	dd = directional_derivative(g_f1, g_f2, u1, u2)
	norm_factor = sqrt( dot(u1,u1) + dot(u2,u2) )
	return dd / norm_factor
end

function scaled_dd(g_f1::Matrix{Float64}, g_f2::Matrix{Float64}, u1::VVF, u2::VVF)
	dd = directional_derivative(g_f1, g_f2, u1, u2)
	norm_factor = sqrt( dot(u1,u1) + dot(u2,u2) )
	return dd / norm_factor
end

function scaled_dd(g_f1::Matrix{Float64}, g_f2::Matrix{Float64}, u1::Matrix{Float64}, u2::Matrix{Float64})
	dd = directional_derivative(g_f1, g_f2, u1, u2)
	norm_factor = sqrt( sum(u1.*u1) + sum(u2.*u2) )
	return dd / norm_factor
end

# TODO: why exactly do I do this?
# "normalizes" a matrix so (bins * bins) / sum(mat) = 1.0
function normalize!(mat::Matrix{Float64}, dA::Float64)
	num_x, num_y = size(mat)
	#c = (bins * bins) / sum(mat)
	c = 1.0 / (sum(mat) * dA)	# LD 3/02/2017
	for xi = 1:num_x
		for yi = 1:num_y
			mat[xi,yi] *= c
		end
	end
end
function normalize!(mat::Array{Float64,3}, dV::Float64)
	num_x, num_y, num_z = size(mat)
	#c = (bins * bins * bins) / sum(mat)
	c = 1.0 / (sum(mat) * dV)		# LD change 3/02/2017
	for xi = 1:num_x
		for yi = 1:num_y
			for zi = 1:num_z
				mat[xi,yi,zi] *= c
			end
		end
	end
end
export normalize!

"""
`mat2traj(mat::Matrix{Float64})`

Converts a matrix to a vector of Float64 vectors (VVF's).
Assumes the matrix is `N x n`, where `n` is the dimensionality of the state space and `N` is the number of points in the trajectory.
"""
function mat2traj(mat::Matrix{Float64})
	N,n = size(mat)
	traj = Array(Vector{Float64}, N)
	for i = 1:N
		traj[i] = vec(mat[i,:])
	end
	return traj
end

"""
`traj2mat(traj::VVF)`

Converts vector of vector of floats into `N x n` matrix, where `N` is the number of points in trajectory and `n` is the dimensionality of the state.
"""
function traj2mat(traj::VVF)
	N = length(traj)
	n = length(traj[1])
	mat = zeros(N, n)
	for i = 1:N
		mat[i,:] = traj[i]
	end
	return mat
end

export vvvf2vvf
function vvvf2vvf(xd::VVVF, ud::VVVF)
	N = length(ud[1])
	xdcat = hcat(xd...)
	udcat = hcat(ud...)
	x = VVF(N+1)
	u = VVF(N)
	for n = 1:N
		x[n] = vcat(xdcat[n,:]...)
		u[n] = vcat(udcat[n,:]...)
	end
	x[N+1] = vcat(xdcat[N+1,:]...)
	return x,u
end


text:
function load!(ctl::Control, file)
    isfile(file) || error("$file does not exist.")
    LibClingo.clingo_control_load(ctl, file)
end

function add!(ctl::Control, program::String; name::String = "base", parameters::Vector{String} = String[])
    LibClingo.clingo_control_add(ctl, name, parameters, length(parameters), program)
end

function ground!(ctl::Control; parts::Vector{LibClingo.clingo_part} = [LibClingo.clingo_part("base",C_NULL, 0)], callback = C_NULL, callback_data = C_NULL)
    LibClingo.clingo_control_ground(ctl, parts, length(parts), callback, callback_data)
end

function solve!(ctl::Control)
    return Handle(LibClingo.clingo_control_solve(ctl, LibClingo.clingo_solve_mode_yield, C_NULL, 0, C_NULL, C_NULL))
end

function solve!(f::Function, ctl)
    handle = solve!(ctl::Control)
    try
        f(handle)
    finally
        LibClingo.clingo_solve_handle_close(handle)
    end
end


text:
#using Test

#=
Test FPN operations:
=#

@testset "FPN" begin
end


text:
export MatPolynomial

import MultivariateMoments: trimat, SymMatrix, getmat
export getmat

#struct SymMatrix{T} <: AbstractMatrix{T}
#    Q::Vector{T}
#    n
#end
#
## i < j
#function trimap(i, j, n)
#    div(n*(n+1), 2) - div((n-i+1)*(n-i+2), 2) + j-i+1
#end
#
#function trimat(::Type{T}, f, n, σ) where {T}
#    Q = Vector{T}(trimap(n, n, n))
#    for i in 1:n
#        for j in i:n
#            Q[trimap(i, j, n)] = f(σ[i], σ[j])
#        end
#    end
#    SymMatrix{T}(Q, n)
#end
#
#Base.size(Q::SymMatrix) = (Q.n, Q.n)
#
#function Base.getindex(Q::SymMatrix, i, j)
#    Q.Q[trimap(min(i, j), max(i, j), Q.n)]
#end
#function Base.getindex(Q::SymMatrix, k)
#    i, j = divrem(k-1, Q.n)
#    Q[i+1, j+1]
#end
#Base.getindex(Q::SymMatrix, I::Tuple) = Q[I...]
#Base.getindex(Q::SymMatrix, I::CartesianIndex) = Q[I.I]

struct MatPolynomial{T, MT <: AbstractMonomial, MVT <: AbstractVector{MT}} <: AbstractPolynomialLike{T} # should be AbstractPolynomialLike{eltype(T)} but it doesn't work
    Q::SymMatrix{T}
    x::MVT
end
# When taking the promotion of a MatPolynomial of JuMP.Variable with a Polynomial JuMP.Variable, it should be a Polynomial of AffExpr
MP.coefficienttype(::Type{<:MatPolynomial{T}}) where {T} = Base.promote_op(+, T, T)
MP.polynomialtype(::Type{MatPolynomial{T, MT, MVT}}) where {T, MT, MVT} = polynomialtype(MT, coefficienttype(MatPolynomial{T, MT, MVT}))
MP.polynomialtype(::Type{MatPolynomial{T, MT, MVT}}, ::Type{S}) where {S, T, MT, MVT} = polynomialtype(MT, S)

Base.zero(::Type{MatPolynomial{T, MT, MVT}}) where {T, MT, MVT} = MatPolynomial{T, MT, monovectype(MT)}(SymMatrix{T}(T[], 0), emptymonovec(MT))
Base.iszero(p::MatPolynomial) = iszero(polynomial(p))

Base.getindex(p::MatPolynomial, I...) = getindex(p.Q, I...)
Base.copy(p::MatPolynomial) = MatPolynomial(copy(p.Q), copy(p.x))

MultivariateMoments.getmat(p::MatPolynomial{T}) where {T} = p.Q

function _matplus(p::MatPolynomial, q::MatPolynomial)
    @assert p.x == q.x
    MatPolynomial(p.Q+q.Q, p.x)
end

function MatPolynomial{T}(f::Function, x::AbstractVector{MT}, σ) where {T, MT<:AbstractMonomial}
    MatPolynomial{T, MT, monovectype(x)}(trimat(T, f, length(x), σ), x)
end
MatPolynomial{T}(f::Function, x::AbstractVector, σ) where T = MatPolynomial{T}(f, monomial.(x), σ)
function MatPolynomial{T}(f::Function, x::AbstractVector) where T
    σ, X = sortmonovec(x)
    MatPolynomial{T}(f, X, σ)
end
MatPolynomial(f::Function, x) = MatPolynomial{Base.promote_op(f, Int, Int)}(f, x)

function MatPolynomial(Q::AbstractMatrix{T}, x, σ) where T
    MatPolynomial{T}((i,j) -> Q[σ[i], σ[j]], x)
end
function MatPolynomial(Q::AbstractMatrix, x)
    σ, X = sortmonovec(x)
    MatPolynomial(Q, X, σ)
end

#function Base.convert{T, PT <: AbstractPolynomial{T}}(::Type{PT}, p::MatPolynomial)
#    # coefficienttype(p) may be different than T and polynomial(p) may be different than PT (different module)
#    convert(PT, polynomial(p))
#end
function MP.polynomial(p::MatPolynomial)
    polynomial(getmat(p), p.x)
end
function MP.polynomial(p::MatPolynomial, ::Type{S}) where {S}
    polynomial(getmat(p), p.x, S)
end

const APL = AbstractPolynomialLike

Base.:(+)(x::APL, y::MatPolynomial) = x + polynomial(y)
Base.:(+)(x::MatPolynomial, y::APL) = polynomial(x) + y
Base.:(+)(x::MatPolynomial, y::MatPolynomial) = polynomial(x) + polynomial(y)
Base.:(-)(x::APL, y::MatPolynomial) = x - polynomial(y)
Base.:(-)(x::MatPolynomial, y::APL) = polynomial(x) - y
Base.:(-)(x::MatPolynomial, y::MatPolynomial) = polynomial(x) - polynomial(y)
Base.:(==)(p::APL, q::MatPolynomial) = p == polynomial(q)
Base.:(==)(p::MatPolynomial, q::MatPolynomial) = iszero(p - q)
function Base.isapprox(p::MatPolynomial, q::MatPolynomial; kwargs...)
    p.x == q.x && isapprox(p.Q, q.Q; kwargs...)
end


text:
abstract type Particles end

export HardSpheres

struct HardSpheres <: Particles

    n
    q
    v
    ϵ

    function HardSpheres(q :: Vector{Vector{Float64}}, v, ϵ :: AbstractFloat)

        n = length(q)

        @assert n == length(v)

        new(n, q, v, ϵ)

    end

    function HardSpheres(rng, n, ϵ)

        q = Vector{Float64}[]

        J4 = [1 0; 0 1]

        push!(q, [0.5000, 0.5000])
        push!(q, [0.1622, 0.4505])
        push!(q, [0.3112, 0.2290])
        push!(q, [0.5285, 0.9133])
        push!(q, [0.1656, 0.1524])
        push!(q, [0.6020, 0.8258])
        push!(q, [0.2630, 0.5383])
        push!(q, [0.3000, 0.8000])
        push!(q, [0.6892, 0.0782])
        push!(q, [0.7482, 0.4427])

        p = zeros(2)

        for klm in 11:n
            frein = 0
            overlap = 1
            while (overlap == 1) && frein < 10000

                p0 = [2ϵ + (1 - 4ϵ) * rand(rng), 2ϵ + (1 - 4ϵ) * rand(rng)]
                p = J4 * p0
                overlap2 = 0
                for subh in 1:(klm-1)
                    if norm(p - q[subh], 2) < 2ϵ
                        overlap2 = 1
                    end
                end

                overlap = overlap2
                frein += 1
            end

            frein == 10000 && (@error "Echec tirage initial")

            push!(q, p)

        end

        v = [ randn(rng, 2) for j in 1:n ]

        new( n, q, v, ϵ )

    end

end

const offset = [[0, 0], [1,  0], [-1,  0],
                [0, 1], [0, -1], [-1,  1],
                [1, 1], [1, -1], [-1, -1]]

function compute_dt( p :: HardSpheres, i, j, k = 1 )

    ϵ = p.ϵ

    δz = p.q[i] .+ offset[k] .- p.q[j]
    δv = p.v[i] .- p.v[j]

    c = (δv'δz).^2 - (δv'δv) * (δz'δz - (2ϵ)^2)

    if δv'δz >= 0
        δt = Inf
    elseif c < 0
        δt = Inf
    else
        δt = -(δv'δz + sqrt(c)) / (δv'δv)
    end

    return δt

end

function compute_dt(hs :: HardSpheres, i)

    ϵ = hs.ϵ
    x = hs.q[i]
    v = hs.v[i]

    if v[1] > 0
        s1 = (1 - ϵ - x[1]) / v[1]
    elseif v[1] < 0
        s1 = (ϵ - x[1]) / v[1]
    elseif v[1] == 0
        s1 = Inf
    end

    if v[2] > 0
        s2 = (1 - ϵ - x[2]) / v[2]
    elseif v[2] < 0
        s2 = (ϵ - x[2]) / v[2]
    elseif v[2] == 0
        s2 = Inf
    end

    r = min(s1, s2)

    if r == s1
        return s1, 1
    else
        return s2, 2
    end

end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ 81cdc1aa-04fd-11eb-3235-bff0e3f52535
using Pkg, DrWatson

# ╔═╡ 81cdee8c-04fd-11eb-29f5-85927396b6dc
begin
	@quickactivate "StatisticalRethinkingTuring"
	using Turing
	using StatisticalRethinking
end

# ╔═╡ 5e06a11c-04f6-11eb-0116-1fd51d614f99
md"## clip-04-72-79t"

# ╔═╡ 82a6656e-04fd-11eb-08c9-3bd9d2d8af20
md"### snippet 4.72"

# ╔═╡ 8ad73da4-0529-11eb-3e51-910e74043bfe
df = CSV.read(sr_datadir("cherry_blossoms.csv"), DataFrame; missingstring = "NA");

# ╔═╡ af9e141e-0529-11eb-3c28-b986f5939349
scatter(df.year, df.doy, leg=false)

# ╔═╡ 82ba1b72-04fd-11eb-33ed-0db4307c66fc
md"### snippet 4.73"

# ╔═╡ d9b4a9fa-0529-11eb-21b5-2da845b95fe9
begin
	dropmissing!(df, :doy)
	num_knots = 15
	knot_list = quantile(df.year, range(0, 1, length = num_knots))
	basis = BSplineBasis(4, knot_list)
	B = basismatrix(basis, df.year)
end;

# ╔═╡ b0b1e776-0575-11eb-093b-f1cb2d9c4b2e
begin
	plot(legend = false, xlabel = "year", ylabel = "basis value")
	for y in eachcol(B)
		plot!(df.year, y)
	end
	plot!()
end

# ╔═╡ 832a285e-04fd-11eb-061b-49a40155abe1
md"## snippet 4.76"

# ╔═╡ cdb02de2-0575-11eb-121e-df6ecd57ba7c
@model function spline_1(D, B = B)
    α ~ Normal(100, 10)
    w ~ filldist(Normal(0, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = α .+ B * w
    D ~ MvNormal(μ, σ)
    return μ
end

# ╔═╡ 8343fe0a-04fd-11eb-2a3e-79daed2e6b98
md"### snippet 4.77"

# ╔═╡ 1233198e-0576-11eb-3c17-672034589062
begin
	quap4_7t = quap(spline_1(df.doy))
	quap4_7t.coef
end

# ╔═╡ 833793e0-04fd-11eb-0706-7567eada0aed
begin
	w_str = ["w[$i]" for i in 1:length(basis)]
	post_3 = DataFrame(rand(quap4_7t.distr, 1000)', ["α"; w_str; "σ"])
	w_3 = mean.(eachcol(post_3[:, w_str]))              # either
	w_3 = [mean(post_3[:, col]) for col in w_str]       # or
	plot(legend = false, xlabel = "year", ylabel = "basis * weight")
	for y in eachcol(B .* w_3')
		plot!(df3.year, y)
	end
	plot!()
end

# ╔═╡ 83b54bb4-04fd-11eb-3208-4162cb3eabc7
md"### snippet 4.78"

# ╔═╡ d9fbed6e-0572-11eb-29bd-1d30ae286045
begin
	mu_3 = post_3.α' .+ B * Array(post_3[!, w_str])'
	mu_3 = meanlowerupper(mu_3)
	scatter(df.year, df.doy, alpha = 0.3)
	plot!(df.year, mu_3.mean, ribbon = (mu_3.mean .- mu_3.lower, mu_3.upper .- mu_3.mean))
end

# ╔═╡ 83f2d7f4-04fd-11eb-3114-d3faa79abb1a
md"### snippet 4.79"

# ╔═╡ 21d2a6a0-0573-11eb-1043-adeebe1e8fde
@model function spline_2(D, B = B)
    α ~ Normal(100, 10)
    w ~ filldist(Normal(0, 10), size(B, 2))
    σ ~ Exponential(1)
    μ = [α + sum(Brow .* w) for Brow in eachrow(B)]
    D ~ MvNormal(μ, σ)
end

# ╔═╡ 840196ae-04fd-11eb-1c8a-3511b18e2835
begin
	quap_alt_4_7t = quap(spline_2(df.doy))
	quap_alt_4_7t.coef
end

# ╔═╡ 4f386d5c-0573-11eb-084c-1758217d06f9
md"## End of clip-04-72-79t.jl"

# ╔═╡ Cell order:
# ╟─5e06a11c-04f6-11eb-0116-1fd51d614f99
# ╠═81cdc1aa-04fd-11eb-3235-bff0e3f52535
# ╠═81cdee8c-04fd-11eb-29f5-85927396b6dc
# ╟─82a6656e-04fd-11eb-08c9-3bd9d2d8af20
# ╠═8ad73da4-0529-11eb-3e51-910e74043bfe
# ╠═af9e141e-0529-11eb-3c28-b986f5939349
# ╟─82ba1b72-04fd-11eb-33ed-0db4307c66fc
# ╠═d9b4a9fa-0529-11eb-21b5-2da845b95fe9
# ╠═b0b1e776-0575-11eb-093b-f1cb2d9c4b2e
# ╟─832a285e-04fd-11eb-061b-49a40155abe1
# ╠═cdb02de2-0575-11eb-121e-df6ecd57ba7c
# ╟─8343fe0a-04fd-11eb-2a3e-79daed2e6b98
# ╠═1233198e-0576-11eb-3c17-672034589062
# ╠═833793e0-04fd-11eb-0706-7567eada0aed
# ╟─83b54bb4-04fd-11eb-3208-4162cb3eabc7
# ╠═d9fbed6e-0572-11eb-29bd-1d30ae286045
# ╟─83f2d7f4-04fd-11eb-3114-d3faa79abb1a
# ╠═21d2a6a0-0573-11eb-1043-adeebe1e8fde
# ╠═840196ae-04fd-11eb-1c8a-3511b18e2835
# ╟─4f386d5c-0573-11eb-084c-1758217d06f9


text:
# distribution(f::F, args) where F = djltype(F)(args...)
# djltype(::typeof(normal)) = Distributions.Normal
# djltype(::typeof(betarv)) = Distributions.Beta
# djltype(::typeof(uniform)) = Distributions.Uniform


text:
@testset "304.range-sum-query-2d-immutable.jl" begin
    matrix = NumMatrix([
        [3, 0, 1, 4, 2],
        [5, 6, 3, 2, 1],
        [1, 2, 0, 1, 5],
        [4, 1, 0, 1, 7],
        [1, 0, 3, 0, 5],
    ])
    @test sum_region(matrix, 3, 2, 5, 4) == 8
    @test sum_region(matrix, 2, 2, 3, 3) == 11
    @test sum_region(matrix, 2, 3, 3, 5) == 12
end


text:
using Base: RefValue
using Base.Broadcast: Style, DefaultArrayStyle, Broadcasted, flatten
using LinearAlgebra
using Zygote: pullback

abstract type Model end

# Basic operations

_push_allparams!(::Nothing, a::AbstractArray{T}) where T <: Number = AbstractArray{T}[a]
_push_allparams!(c::Vector, a::AbstractArray{T}) where T <: Number = push!(c, a)
_push_allparams!(c, ::Function) = c
_push_allparams!(c, t::Tuple) = begin
    for el in t
        c = _push_allparams!(c, el)
    end
    return c
end
_push_allparams!(c, m::T) where T <: Model = begin
    for k in fieldnames(T)
        c = _push_allparams!(c, getfield(m, k))
    end
    return c
end
allparams(m::Model) = _push_allparams!(nothing, m)

# NOTE: looks like the above definition is faster than what follows
# _yield_allparams(::Function) = ()
# _yield_allparams(t::Tuple) = Iterators.flatten((_yield_allparams(el) for el in t))
# _yield_allparams(a::AbstractArray{T}) where T <: Number = (a,)
# _yield_allparams(m::T) where T <: Model = Iterators.flatten((_yield_allparams(getfield(m, k)) for k in fieldnames(T)))
# allparams(m::Model) = _yield_allparams(m)

Base.vec(m::Model) = vcat((vec(p) for p in allparams(m))...)

overlap(m1::Model, m2::Model) = [
    p1 for (p1, p2) in zip(ProtoGrad.allparams(m1), ProtoGrad.allparams(m2))
    if objectid(p1) == objectid(p2)
]

map_recursively(fun, f::Function) = f
map_recursively(fun, a::AbstractArray) = fun(a)
map_recursively(fun, t::Tuple) = fun.(t)
map_recursively(fun, m::T) where T <: Model = T((map_recursively(fun, getfield(m, k)) for k in fieldnames(T))...)

for fun in (:similar, :copy, :zero)
    _fun = Symbol("_", fun)
    @eval begin
        $_fun(x) = $fun(x)
        $_fun(f::Function) = f
        Base.$fun(m::Model) = map_recursively($_fun, m)
    end
end

# Iteration

Base.IteratorSize(::Type{<:Model}) = Base.HasLength()
Base.length(m::Model) = sum(length(p) for p in allparams(m))

function Base.iterate(m::Model, itr=Iterators.flatten(allparams(m)))
    y = iterate(itr)
    if y === nothing
        return nothing
    end
    val, s = y
    return val, Iterators.rest(itr, s)
end

Base.show(io::IO, m::T) where T <: Model = print(io, T)

# Indexing

Base.firstindex(m::Model) = 1
Base.lastindex(m::Model) = length(m)

Base.getindex(m::Model, i) = begin
    # TODO optimize (no need to call allparams)
    ps = allparams(m)
    j = 1
    while i > length(ps[j])
        i -= length(ps[j])
        j += 1
    end
    return ps[j][i]
end

Base.setindex!(m::Model, v, i) = begin
    # TODO optimize (no need to call allparams)
    ps = allparams(m)
    j = 1
    while i > length(ps[j])
        i -= length(ps[j])
        j += 1
    end
    ps[j][i] = v
end

# Broadcast

Base.axes(::Model) = nothing
Base.Broadcast.broadcastable(m::Model) = m

Base.Broadcast.BroadcastStyle(::Type{T}) where T <: Model = Style{T}()
Base.Broadcast.BroadcastStyle(s::Style{T}, ::DefaultArrayStyle{0}) where T <: Model = s

Base.Broadcast.instantiate(bc::Broadcasted{Style{T}}) where T <: Model = bc

# The idea here is:
#   - have a recursive procedure that goes down to the leaves of the struct
#   - where to descend is dictated by the structure of the destination container
#   - scalars get passed down as they are
#   - once the leaves (numerical arrays or functions) are reached, just use default broadcasting

_getindex(a, idx) = getindex(a, idx)
_getindex(rv::RefValue, idx) = rv
_getindex(n::Number, idx) = n
_getindex(f::Function, idx) = f

_getfield(s, k) = getfield(s, k)
_getfield(rv::RefValue, k) = rv
_getfield(n::Number, k) = n
_getfield(f::Function, k) = f

function recursive_copyto!(a::AbstractArray, f, args)
    a .= f.(args...)
end

recursive_copyto!(f::Function, _, args) = return

function recursive_copyto!(t::Tuple, f, args)
    for (i, el) in Iterators.enumerate(t)
        args_ith = [_getindex(arg, i) for arg in args]
        recursive_copyto!(el, f, args_ith)
    end
end

function recursive_copyto!(m::T, f, args) where T <: Model
    for k in fieldnames(T)
        dest_k = getfield(m, k)
        args_k = [_getfield(arg, k) for arg in args]
        recursive_copyto!(dest_k, f, args_k)
    end
end

function Base.copyto!(dest::T, bc::Broadcasted{Style{T}}) where T <: Model
    flat_bc = flatten(bc)
    recursive_copyto!(dest, flat_bc.f, flat_bc.args)
    dest
end

find_model(bc::Base.Broadcast.Broadcasted) = find_model(bc.args)
find_model(args::Tuple) = find_model(find_model(args[1]), Base.tail(args))
find_model(x) = x
find_model(::Tuple{}) = nothing
find_model(m::Model, rest) = m
find_model(::Any, rest) = find_model(rest)

function Base.copy(bc::Broadcasted{Style{T}}) where T <: Model
    m = find_model(bc)
    dest = similar(m)
    copyto!(dest, bc)
end

# Vector space operations

for op in (:+, :-)
    @eval begin
        function Base.$op(a::T, b::T) where T <: Model
            Base.broadcast($op, a, b)
        end

        function Base.$op(a::Model, b::Number)
            Base.broadcast($op, a, b)
        end

        function Base.$op(a::Number, b::Model)
            Base.broadcast($op, a, b)
        end
    end
end

for op in (:*, :/)
    @eval function Base.$op(a::Model, b::Number)
        Base.broadcast($op, a, b)
    end
end

function Base.:*(a::Number, b::Model)
    Base.broadcast(*, a, b)
end

function Base.:\(a::Number, b::Model)
    Base.broadcast(/, a, b)
end

# Linear algebra

LinearAlgebra.dot(m1::T, m2::T) where T <: Model = sum((dot(p1, p2) for (p1, p2) in zip(allparams(m1), allparams(m2))))


text:
module CurveProximityQueries

using LinearAlgebra
using StaticArrays
using IntervalArithmetic
using DataStructures
import ConvexBodyProximityQueries: closest_points, minimum_distance, collision_detection, tolerance_verification

export closest_points, minimum_distance, collision_detection, tolerance_verification

include("parametric.jl")
include("bounds.jl")
include("bnb.jl")
include("plot.jl")

end # module


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
include("cpu_stencils.jl")

nx = 9
ny = 9
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

original_data = CreateData(dx,dy,dz)
## Stencil Definition
t_cases = [32 16 8 4]
i = 16
radius = i
t_steps = Int(max(t_cases...)/i)
coefs = [Float64(round(1/i, digits=2)) for i = (radius+1):-1:1]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 16, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def)

## Input Data size Definition
data = PadData(radius, original_data)
@time dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)

# t_cases = [32 16 8 4]
# i = t_cases[2]
# radius = i
# t_steps = Int(max(t_cases...)/i)
# coefs = [Float64(round(1/i, digits=2)) for i = (radius+1):-1:1]
# star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 16, c[i]*(
#             D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
#             D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
# st_def = CreateStencilDefinition(star_stencil, coefs)
# st_inst = NewStencilInstance(st_def)
#
# ## Input Data size Definition
# data = PadData(radius, original_data)
# @time dout = ApplyStencil(st_inst, data, t_steps)
# gpu_out = Array(dout)
#
# t_cases = [32 16 8 4]
i = t_cases[3]
radius = i
t_steps = Int(max(t_cases...)/i)
coefs = [Float64(round(1/i, digits=2)) for i = (radius+1):-1:1]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 8, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def)

## Input Data size Definition
data = PadData(radius, original_data)
@time dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)

t_cases = [32 16 8 4]
i = t_cases[4]
radius = i
t_steps = Int(max(t_cases...)/i)
coefs = [Float64(round(1/i, digits=2)) for i = (radius+1):-1:1]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 4, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_inst = NewStencilInstance(st_def)

## Input Data size Definition
data = PadData(radius, original_data)
@time dout = ApplyStencil(st_inst, data, t_steps)
gpu_out = Array(dout)


text:
using DependencyWalker
using Test
using ObjectFile, Pango_jll

@testset "DependencyWalker.jl" begin
    pango = Library(Pango_jll.libpango_path)
    @show pango
    if Sys.islinux() && Sys.WORD_SIZE == 32
        @test_broken pango isa Library{<:ObjectHandle}
    else
        @test pango isa Library{<:ObjectHandle}
    end
    @test Library("this does not exist.foo") isa Library{Missing}
end


text:
"""
$(TYPEDEF)

Pose2 is a SE(2) mechanization of two Euclidean translations and one Circular rotation, used for general 2D SLAM.
"""
struct Pose2 <: IncrementalInference.InferenceVariable end
getDimension(::Pose2) = 3
getManifolds(::Pose2) = (:Euclid,:Euclid,:Circular)


text:
using Term
using Test

nlines(x) = length(split(x, "\n"))
lw(x) = max(length.(split(x, "\n"))...)

using TimerOutputs: TimerOutputs, @timeit
const TIMEROUTPUT = TimerOutputs.TimerOutput()

macro timeit_include(path::AbstractString)
    :(@timeit TIMEROUTPUT $path include($path))
end

@timeit_include("00_empty.jl")

print("\n\n")
@timeit_include("01_test_text_utils.jl")

print("\n\n")
@timeit_include("02_test_ansi.jl")  # ansi & color

print("\n\n")
@timeit_include("03_test_measure.jl")

print("\n\n")
@timeit_include("04_test_markup.jl")

print("\n\n")
@timeit_include("05_test_macros.jl")

print("\n\n")
@timeit_include("06_test_renderables.jl")

print("\n\n")
@timeit_include("07_test_layout.jl")

print("\n\n")
@timeit_include("08_test_panel.jl")

print("\n\n")
@timeit_include("09_test_inspect.jl")

show(TIMEROUTPUT; compact = true, sortby = :firstexec)
println("\n")


text:
module MiniFB
using Colors
using MiniFB_jll

using CEnum

include("ctypes.jl")
export Ctm, Ctime_t, Cclock_t

include("libminifb_common.jl")
include("libminifb_api.jl")
include("libminifb_docs.jl")
include("helpers.jl")

end


text:
#=
Based on the file:
http://ls11-www.cs.uni-dortmund.de/_media/rudolph/hypervolume/hv_python.zip
originally written in python by Simon Wessing (2010).

This implementation is for the variant 3 detailed in:
C. M. Fonseca, L. Paquete, and M. Lopez-Ibanez. An improved dimension-sweep
algorithm for the hypervolume indicator. In IEEE Congress on Evolutionary
Computation, pages 1157-1163, Vancouver, Canada, July 2006.
=#

mutable struct Node
    cargo::Vector
    next::Vector{Union{Nothing,Node}}
    prev::Vector{Union{Nothing,Node}}
    ignore::Int
    area::Vector{Float64}
    volume::Vector{Float64}
end

#=
function Base.show(io::IO, node::Node)
    print(io, "Node( ")
    print(io, "cargo = ", node.cargo)
    print(io, ", ")
    for n in node.next
        if isnothing(n)
            print(io, n, " -> ")
            continue
        end

        print(io, n.cargo, " -> ")
    end
    print(io, ")")

end
=#

function Node(numberLists, cargo = zeros(0);
        next   = fill(nothing, numberLists),
        prev   = fill(nothing, numberLists),
        ignore = 0,
        area   = zeros(numberLists),
        volume = zeros(numberLists),
    )

    Node(cargo, next, prev, ignore, area, volume)

end

function extend!(sentinel::Node, nodes, index)
    for node in nodes
        lastButOne = sentinel.prev[index]
        node.next[index] = sentinel
        node.prev[index] = lastButOne
        # set the last element as the new one
        sentinel.prev[index] = node
        lastButOne.next[index] = node
    end
end

"""
    Removes and returns 'node' from all lists in 1,2,..., 'index'.
"""
function remove!(node::Node, index, bounds)
    for i in 1:index-1
        predecessor = node.prev[i]
        successor = node.next[i]
        predecessor.next[i] = successor
        successor.prev[i] = predecessor

        if bounds[i] > node.cargo[i]
            bounds[i] = node.cargo[i]
        end
    end
    return node
end

"""
Inserts 'node' at the position it had in all lists in [1, 'index']
before it was removed. This method assumes that the next and previous
nodes of the node that is reinserted are in the list.

"""
function reinsert!(node::Node, index, bounds)
    for i in 1:index
        node.prev[i].next[i] = node
        node.next[i].prev[i] = node
        if bounds[i] > node.cargo[i]
            bounds[i] = node.cargo[i]
        end
    end

    node
end


text:
function plot_rwd_dists(net_num,xQ,xP,c_data,t_dist,o;fldr::String="figs",fname::String="",outcome::Symbol=:fail)
    cm = mean(c_data)
    cs = std(c_data)

    tm = t_dist[:X3_1][1]
    ts = t_dist[:X3_2][1]

    pygui(false)
    f,ax = subplots(1,2)
    f[:set_size_inches](8,4)

    ax[1][:hist](c_data)
    ax[1][:set_title](@sprintf("Candidate Dist. Drew: %0.2f",o))
    mc = @sprintf("Mean: %0.2f, Var: %0.2f",cm,cs)
    ax[1][:axvline](cm,color=:red,label=mc)
    ax[1][:legend]()

    ax[2][:hist](rand(Normal(tm,ts),1000))
    mt = @sprintf("Mean: %0.2f, Var: %0.2f",tm,ts)
    ax[2][:axvline](tm,color=:red,label=mt)
    ax[2][:set_title]("Trusted Distribution")
    ax[2][:legend]()

    f[:suptitle](@sprintf("Network: %d---Outcome %s---xQ: %0.2f, xP: %0.2f",net_num,string(outcome),xQ,xP))
    #  f[:tight_layout]()
    if fname == ""
        fname = "net_num$net_num"
    end

    if ~isdir(fldr)
        println("making folder $fldr")
        Base.Filesystem.mkdir(fldr)
    end

    file_name = joinpath(fldr,fname)
    println("saving to: $file_name")
    PyPlot.savefig("$file_name.pdf",dpi=300,transparent=true)
end

max_stars_repo_name: COHRINT/FaMSeC

max_issues_repo_name: COHRINT/FaMSeC

max_forks_repo_name: COHRINT/FaMSeC


text:
using PNML, EzXML
using AbstractTrees, Test, SafeTestsets
using PrettyPrinting
using IfElse
# Run the tests embedded in docstrings.
using Documenter, LabelledArrays
using JET

using PNML: PNML,
    parse_pnml, parse_str, parse_file, parse_node, @xml_str,
    parse_net, parse_page, parse_place, parse_transition, parse_arc,
    parse_refPlace, parse_refTransition,
    parse_toolspecific, parse_graphics, parse_structure,
    parse_initialMarking, parse_hlinitialMarking,
    parse_inscription, parse_hlinscription,
    parse_condition,
    parse_declaration, parse_sort, parse_term, parse_label,
    parse_tokengraphics, parse_tokenposition, parse_name

using PNML: tag, xmlnode, Maybe, name,
    Place, Transition, Arc, Page, PnmlNet,
    pid, place_ids, transition_ids, arc_ids,
    reftransition_ids, refplace_ids,
    nets, pages, places, transitions, arcs,
    place, transition, arc,
    has_place, has_transition, has_arc, has_name,
    first_net, firstpage,
    tools, has_tools, get_toolinfo,
    marking, condition, conditions, inscription,
    labels, has_label, has_labels,
    nettype,
    default_marking, default_condition, default_inscription,
    AnyElement, ToolInfo, version,
    SimpleNet

const GROUP        = get(ENV, "GROUP", "All")
const PRINT_PNML   = parse(Bool, get(ENV, "PRINT_PNML", "true"))
const VERBOSE_PNML = parse(Bool, get(ENV, "VERBOSE_PNML", "true"))

# Use default display width for printing.
if !haskey(ENV, "COLUMNS")
    ENV["COLUMNS"] = 180
end

const testdir = dirname(@__FILE__)
const pnml_dir = joinpath(@__DIR__, "data")

"Turn string into XML node."
to_node(s::AbstractString) = root(EzXML.parsexml(s))

"Print `node` prepended by optional label string."
function printnode(io::IO, node; label=nothing, kw...)
    if PRINT_PNML
        !isnothing(label) && print(io, label, " ")
        show(io, MIME"text/plain"(), node)
        println(io, "\n")
    end
end
function printnode(n; kw...)
    printnode(stdout, n; kw...)
end

function printnodeln(io::IO, n; kw...)
    printnode(io, n; kw...)
    PRINT_PNML && println(io)
end

header(s) = if VERBOSE_PNML
    println("##### ", s)
end

const SHOW_SUMMARYSIZE = get(ENV, "SHOW_SUMMARYSIZE", "false") == "true"
function showsize(ob,k)
    if SHOW_SUMMARYSIZE && PRINT_PNML
        summarysz = Base.summarysize(ob[k])
        @show k,summarysz
    end
end

"Return true if one of the GROUP environment variable's values if found in 'v'."
select(v...) = any(any(==(g), v) for g in split(GROUP))

if select("None")
    return
end

@testset verbose=false "PNML.jl" begin
    if select("All", "Doc")
        header("Doctests")
        @testset "doctest" begin doctest(PNML, manual = true) end
    end
    if select("All", "Base")
        header("Base")
        @testset "maps"         begin include("maps.jl") end
        @testset "utils"        begin include("utils.jl") end
    end
    if select("All", "IR")
        header("IR")
        @testset "nodes"        begin include("nodes.jl") end
        @testset "labels"       begin include("labels.jl") end
        @testset "graphics"     begin include("graphics.jl") end
        @testset "parse_labels" begin include("parse_labels.jl") end
        @testset "toolspecific" begin include("toolspecific.jl") end
        @testset "exceptions"   begin include("exceptions.jl") end
        @testset "parse_tree"   begin include("parse_tree.jl") end
        @testset "pages"        begin include("pages.jl") end
        @testset "flatten"      begin include("flatten.jl") end
        @testset "example file" begin include("parse_examples.jl") end
    end
    if select("All", "Net")
        header("Net")
        @testset "document"     begin include("document.jl") end
        @testset "simplenet"    begin include("simplenet.jl") end
    end
end


text:
const ASCII = raw"""
  _
 | |
 | |    _   _  ___ ___ _   _ _ __ ___
 | |   | | | |/ __/ _ \ | | | '_ ` _ \
 | |___| |_| | (_|  __/ |_| | | | | | |
 |______\__, |\___\___|\__,_|_| |_| |_|
         __/ |
        |___/
"""

function str2vec(s::String, len::Int)
    str = zeros(UInt8, len)
    str[1:length(s)] = codeunits(s)
    return str
end

function init_figsensor!(figsensor::Ref{mjvFigure})
    figsensor[].flg_extend = 1
    figsensor[].flg_barplot = 1
    figsensor[].title = str2vec("Sensor Data", length(figsensor[].title))
    figsensor[].yformat = str2vec("%.0f", length(figsensor[].yformat))
    figsensor[].gridsize = [2, 3]
    figsensor[].range = [[0 1], [-1 1]]
    return figsensor
end

@inline inc(x::Integer, min::Integer, max::Integer) = ifelse(x == max, min, x + 1)
@inline dec(x::Integer, min::Integer, max::Integer) = ifelse(x == min, max, x - 1)

function startffmpeg(w::Integer, h::Integer, rout::Integer; squashoutput::Bool = true)
    w > 0 && h > 0 || error("w and h must be > 0")

    dst = tempname() * ".mp4"
    arg = `-y -f rawvideo -pixel_format rgb24 -video_size $(w)x$(h) -use_wallclock_as_timestamps true -i pipe:0 -c:v libx264 -preset ultrafast -tune animation -vf "vflip" -r $rout $dst`

    FFMPEG.@ffmpeg_env begin
        in = Base.PipeEndpoint()
        if squashoutput
            p = Base._spawn(`$(FFMPEG.ffmpeg) $arg`, Any[in, devnull, devnull])
        else
            p = Base._spawn(`$(FFMPEG.ffmpeg) $arg`, Any[in, stdout, stderr])
        end
        p.in = in
        return p, dst
    end
end

function safe_unlock(lck::ReentrantLock)
    if islocked(lck) && current_task() === lck.locked_by
        unlock(lck)
    end
    return
end

"""
    spinwait(delay)

Spin in a tight loop for at least `delay` seconds.

Note this function is only accurate on the order of approximately `@elapsed time()`
seconds.
"""
@inline function spinwait(dt::Real)
    dt > 0 || error("dt must be > 0")
    t0 = time()
    while time() - t0 < dt end
    return nothing
end

@inline function str2unicode(s::AbstractString)
    length(s) == 1 || error("s must be a single length string")
    return Int(first(s))
end

# From PrettyTables.jl: https://github.com/ronisbr/PrettyTables.jl/blob/d8a948a2dd14df66852c9702387c13f4dfce30d8/src/backends/text/print.jl#L50
function get_terminalsize()
    return withenv("LINES" => -1, "COLUMNS" => -1) do
        displaysize(stdout)
    end
end


text:
function EGMFormatReader(arg0::JString, arg1::jboolean)
    return EGMFormatReader((JString, jboolean), arg0, arg1)
end

function EGMFormatReader(arg0::JString, arg1::jboolean, arg2::jboolean)
    return EGMFormatReader((JString, jboolean, jboolean), arg0, arg1, arg2)
end

function get_max_available_degree(obj::PotentialCoefficientsReader)
    return jcall(obj, "getMaxAvailableDegree", jint, ())
end

function get_max_available_order(obj::PotentialCoefficientsReader)
    return jcall(obj, "getMaxAvailableOrder", jint, ())
end

function get_max_parse_degree(obj::PotentialCoefficientsReader)
    return jcall(obj, "getMaxParseDegree", jint, ())
end

function get_max_parse_order(obj::PotentialCoefficientsReader)
    return jcall(obj, "getMaxParseOrder", jint, ())
end

function get_provider(obj::EGMFormatReader, arg0::jboolean, arg1::jint, arg2::jint)
    return jcall(obj, "getProvider", RawSphericalHarmonicsProvider, (jboolean, jint, jint), arg0, arg1, arg2)
end

function get_supported_names(obj::PotentialCoefficientsReader)
    return jcall(obj, "getSupportedNames", JString, ())
end

function load_data(obj::EGMFormatReader, arg0::InputStream, arg1::JString)
    return jcall(obj, "loadData", void, (InputStream, JString), arg0, arg1)
end

function missing_coefficients_allowed(obj::PotentialCoefficientsReader)
    return jcall(obj, "missingCoefficientsAllowed", jboolean, ())
end

function set_max_parse_degree(obj::PotentialCoefficientsReader, arg0::jint)
    return jcall(obj, "setMaxParseDegree", void, (jint,), arg0)
end

function set_max_parse_order(obj::PotentialCoefficientsReader, arg0::jint)
    return jcall(obj, "setMaxParseOrder", void, (jint,), arg0)
end

function still_accepts_data(obj::PotentialCoefficientsReader)
    return jcall(obj, "stillAcceptsData", jboolean, ())
end


text:
"""this is used to download the data from biomd and the sbml-test-suite"""
module SBMLBioModelsRepository

const datadir = joinpath(@__DIR__, "../data")

# function __init__()
# using Pkg
# Pkg.add(url="https://github.com/paulflang/SBML.jl/", rev="pl/mk-species-units")
# end

using CSV, DataFrames, JSON3, JSONTables, Glob
using Base.Threads, Base.Iterators, Downloads

# using SBML
using ModelingToolkit, OrdinaryDiffEq, CSV, DataFrames

# include("lower.jl")
include("biomd.jl")
include("suite.jl")

# "do it for me" functions
export sbml_test_suite, biomodels

export curl_biomd_xmls
export curl_biomd_metadata, biomd_metadata, curl_biomd_zips, biomd_zip_urls, unzip_biomd
export get_sbml_suite_fns, jsonfn_to_df
export datadir
# export goodbad, test_suite, lower_one, lower_fns

end


text:
# Autogenerated wrapper script for CRlibm_jll for powerpc64le-linux-gnu
export libcrlibm

JLLWrappers.@generate_wrapper_header("CRlibm")
JLLWrappers.@declare_library_product(libcrlibm, "libcrlibm.so")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcrlibm,
        "lib/libcrlibm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
module MPI_Tools

using FFTW, Optim, HDF5, PyPlot, Unitful
const μ₀ = (4*π*1e-7)u"T*m/A"
const ϵ₀ = (8.85e-12)u"F/m"
mu0 = μ₀
include("MyFFT.jl")
include("THD.jl")
include("MyBisect.jl")
include("meshgrid.jl")
UL(X) = X |> NoUnits

export meshgrid,
       FFT_1s,
       FFT_1s_to2s,
       μ₀,mu0,ϵ₀,
       mybisect,
       THD_TimeSeries,
       UL

# Write your package code here.

end


text:
@inline column_name(query_info::OraQueryInfo) = unsafe_string(query_info.name, query_info.name_length)

function ncol(stmt::QueryStmt)
    if stmt.columns_info == nothing
        init_columns_info!(stmt)
    end
    return length(stmt.columns_info)
end

function columns_info(stmt::QueryStmt) :: Vector{OraQueryInfo}
    if stmt.columns_info == nothing
        init_columns_info!(stmt)
    end

    return stmt.columns_info
end

function column_info(stmt::QueryStmt, column_index::Integer) :: OraQueryInfo
    if stmt.columns_info == nothing
        init_columns_info!(stmt)
    end

    return stmt.columns_info[column_index]
end

@inline function oracle_type(stmt::QueryStmt, column_index::Integer) :: OraOracleTypeNum
    return column_info(stmt, column_index).type_info.oracle_type_num
end

function init_columns_info!(stmt::QueryStmt)

    function get_stmt_num_columns(ctx::Context, stmt_handle::Ptr{Cvoid})
        num_columns_ref = Ref{UInt32}(0)
        result = dpiStmt_getNumQueryColumns(stmt_handle, num_columns_ref)
        error_check(ctx, result)
        return num_columns_ref[]
    end

    function get_query_column_info(ctx::Context, stmt_handle::Ptr{Cvoid}, column_index::Integer)
        query_info_ref = Ref{OraQueryInfo}()
        result = dpiStmt_getQueryInfo(stmt_handle, UInt32(column_index), query_info_ref)
        error_check(context(stmt), result)
        return query_info_ref[]
    end

    columns_info = Vector{OraQueryInfo}()
    for col in 1:get_stmt_num_columns(context(stmt), stmt.handle)
        push!(columns_info, get_query_column_info(context(stmt), stmt.handle, col))
    end

    stmt.columns_info = columns_info
    nothing
end

function Stmt(connection::Connection, sql::AbstractString; scrollable::Bool=false, tag::AbstractString="", fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE)
    stmt_handle_ref = Ref{Ptr{Cvoid}}()
    result = dpiConn_prepareStmt(connection.handle, scrollable, String(sql), String(tag), stmt_handle_ref)
    error_check(connection.context, result)
    return Stmt(connection, stmt_handle_ref[], scrollable, fetch_array_size=fetch_array_size)
end

function stmt(f::Function, connection::Connection, sql::AbstractString; scrollable::Bool=false, tag::AbstractString="", fetch_array_size::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE)
    stmt = Stmt(connection, sql, scrollable=scrollable, tag=tag, fetch_array_size=fetch_array_size)

    try
        f(stmt)
    finally
        close(stmt)
        destroy!(stmt)
    end

    nothing
end

"Number of affected rows in a DML statement."
function row_count(stmt::Stmt)
    row_count_ref = Ref{UInt64}()
    result = dpiStmt_getRowCount(stmt.handle, row_count_ref)
    error_check(context(stmt), result)
    return row_count_ref[]
end

"""
    execute(stmt::Stmt; exec_mode::dpiExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32

Returns the number of columns which are being queried.
If the statement does not refer to a query, the value is set to 0.
"""
function execute(stmt::Stmt; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    num_columns_ref = Ref{UInt32}(0)
    result = dpiStmt_execute(stmt.handle, exec_mode, num_columns_ref)
    error_check(context(stmt), result)
    return num_columns_ref[]
end

"""
    execute(connection::Connection, sql::AbstractString;
        scrollable::Bool=false,
        tag::AbstractString="",
        exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT
    ) :: UInt32

Execute a single `sql` statement.

Returns the number of columns which are being queried.
If the statement does not refer to a query, the value is set to 0.
"""
function execute(connection::Connection, sql::AbstractString; scrollable::Bool=false, tag::AbstractString="", exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    local result::UInt32

    stmt(connection, sql, scrollable=scrollable, tag=tag) do stmt
        result = execute(stmt, exec_mode=exec_mode)
    end

    return result
end

function execute_script(connection::Connection, filepath::AbstractString; scrollable::Bool=false, tag::AbstractString="", exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    @assert isfile(filepath) "Couldn't find script file $filepath."

    local sql::String
    open(filepath, "r") do io
        sql = read(io, String)
    end

    return execute(connection, sql, scrollable=scrollable, tag=tag, exec_mode=exec_mode)
end

# execute many
function execute_many(connection::Connection, sql::AbstractString, columns::Vector; scrollable::Bool=false, tag::AbstractString="", exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    local result::UInt32

    stmt(connection, sql, scrollable=scrollable, tag=tag) do stmt
        result = execute(stmt, columns, exec_mode=exec_mode)
    end

    return result
end

# execute many
function execute_many(stmt::Stmt, input_columns::Vector; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    @assert !isempty(input_columns) "Cannot bind empty columns to statement."
    @assert eltype(input_columns) <: Vector "`input_columns` argument is expected to be a vector of vectors."

    function check_columns_length(columns::Vector)
        columns_count = length(columns)

        if columns_count <= 1
            return
        end

        first_column_length = length(columns[1])

        for i in 2:length(columns)
            @inbounds @assert length(columns[i]) == first_column_length
        end

        nothing
    end

    check_columns_length(input_columns)

    for (c, column) in enumerate(input_columns)
        stmt[c] = Variable(stmt.connection, column)
    end

    # execute
    num_iters = length(input_columns[1])
    execute_many(stmt, num_iters, exec_mode | ORA_MODE_EXEC_ARRAY_DML_ROWCOUNTS)

    # get row counts
    function row_counts(stmt::Stmt) :: Vector{UInt64}
        row_counts_length_ref = Ref{UInt32}()
        row_counts_array_ref = Ref{Ptr{UInt64}}()
        result = dpiStmt_getRowCounts(stmt.handle, row_counts_length_ref, row_counts_array_ref)
        error_check(context(stmt), result)

        row_counts_length = row_counts_length_ref[]
        row_counts_array_as_ptr = row_counts_array_ref[]

        row_counts = Vector{UInt64}(undef, row_counts_length)
        for i in 1:row_counts_length
            @inbounds row_counts[i] = unsafe_load(row_counts_array_as_ptr, i)
        end

        return row_counts
    end

    return sum(row_counts(stmt))
end

# DEPRECATE
function execute(connection::Connection, sql::AbstractString, columns::Vector; scrollable::Bool=false, tag::AbstractString="", exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    execute_many(connection, sql, columns, scrollable=scrollable, tag=tag, exec_mode=exec_mode)
end

# DEPRECATE
function execute(stmt::Stmt, input_columns::Vector; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) :: UInt32
    execute_many(stmt, input_columns, exec_mode=exec_mode)
end

@inline function execute_many(stmt::Stmt, num_iters::Integer, exec_mode::OraExecMode)
    result = dpiStmt_executeMany(stmt.handle, exec_mode, UInt32(num_iters))
    error_check(context(stmt), result)
    return
end

function execute_many(stmt::Stmt,
        num_iters::Integer,
        variables::Dict{K,V}
        ; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) where {K<:Union{Integer, Symbol}, V<:Variable}

    for (k, v) in variables
        stmt[k] = v
    end

    execute_many(stmt, num_iters, exec_mode)
end

function execute_many(stmt::Stmt,
        num_iters::Integer,
        variables::Vector{V}
        ; exec_mode::OraExecMode=ORA_MODE_EXEC_DEFAULT) where {V<:Variable}

    for (k, v) in enumerate(variables)
        stmt[k] = v
    end

    execute_many(stmt, num_iters, exec_mode)
end

function close(stmt::Stmt; tag::AbstractString="")
    if stmt.is_open
        result = dpiStmt_close(stmt.handle, tag=String(tag))
        error_check(context(stmt), result)
        stmt.is_open = false
    end
    nothing
end

function fetch_array_size(stmt::Stmt) :: UInt32
    array_size_ref = Ref{UInt32}()
    result = dpiStmt_getFetchArraySize(stmt.handle, array_size_ref)
    error_check(context(stmt), result)
    return array_size_ref[]
end

"""
    fetch_array_size!(stmt::Stmt, new_size::Integer)

Sets the array size used for performing fetches.
All variables defined for fetching must have this many (or more) elements allocated for them.
The higher this value is the less network round trips are required to fetch rows from the database
but more memory is also required.

A value of zero will reset the array size to the default value of DPI_DEFAULT_FETCH_ARRAY_SIZE.
"""
function fetch_array_size!(stmt::Stmt, new_size::Integer)
    result = dpiStmt_setFetchArraySize(stmt.handle, UInt32(new_size))
    error_check(context(stmt), result)
    nothing
end

reset_fetch_array_size!(stmt::Stmt) = fetch_array_size!(stmt, UInt32(0))

"""
    fetch(stmt::Stmt) :: FetchResult

Fetches a single row from the statement.
"""
function fetch(stmt::Stmt) :: FetchResult
    found_ref = Ref{Int32}(0)
    buffer_row_index_ref = Ref{UInt32}(0) # This index is used as the array position for getting values from the variables that have been defined for the statement.
    result = dpiStmt_fetch(stmt.handle, found_ref, buffer_row_index_ref)
    error_check(context(stmt), result)

    local found::Bool = false
    if found_ref[] != 0
        found = true
    end
    return FetchResult(found, buffer_row_index_ref[])
end

function fetchrows(stmt::Stmt, max_rows::Integer=ORA_DEFAULT_FETCH_ARRAY_SIZE) :: FetchRowsResult
    buffer_row_index_ref = Ref{UInt32}()
    num_rows_fetched_ref = Ref{UInt32}()
    more_rows_ref = Ref{Int32}()

    result = dpiStmt_fetchRows(stmt.handle, UInt32(max_rows), buffer_row_index_ref, num_rows_fetched_ref, more_rows_ref)
    error_check(context(stmt), result)
    return FetchRowsResult(buffer_row_index_ref[], num_rows_fetched_ref[], more_rows_ref[])
end

function query_oracle_value(stmt::Stmt, column_index::Integer) :: ExternOracleValue
    native_type_ref = Ref{OraNativeTypeNum}()
    data_handle_ref = Ref{Ptr{OraData}}()
    result = dpiStmt_getQueryValue(stmt.handle, UInt32(column_index), native_type_ref, data_handle_ref)
    error_check(context(stmt), result)

    return ExternOracleValue(stmt, oracle_type(stmt, column_index), native_type_ref[], data_handle_ref[]; use_add_ref=true)
end