text:
using Distributed
# addprocs(8);
# addprocs(4);
# addprocs(16);
addprocs((Sys.CPU_THREADS >> 1) + (1 - nprocs()) # master, and 1 slave for every 2 threads

# @everywhere begin
#
# # using Pkg
# # include(joinpath(dir("AsymptoticPosteriors"), "examples/NGSCov.jl"));
# end
@everywhere begin
# include("/home/chris/.julia/dev/AsymptoticPosteriors/examples/NGSCov.jl")
includet("/home/chriselrod/.julia/dev/AsymptoticPosteriors/examples/NGSCov.jl")
# include("/home/celrod/.julia/dev/AsymptoticPosteriors/examples/NGSCov.jl")
# using DifferentiableObjects

const truth = CorrErrors((0.15, 0.4), (0.9, 0.95), (0.85,0.97), (0.6, 0.2));
const data = NoGoldDataCorr(truth,160,320,52);
const datacorr = NoGoldData(truth,160,320,52);
x = randnsimd(6);
end
@everywhere begin
@time const ap = AsymptoticPosterior(data, x);
end
@everywhere begin
@time quantile(ap, 0.025)

# @time apc = AsymptoticPosterior(data, randn(6), Val(6));
summary(ap)
end

@everywhere begin
function set_up_sim(truth, n_common=100, n1=1000, n2=30)
    S₁, S₂, C₁, C₂, π₁, π₂ = truth.S[1], truth.S[2], truth.C[1], truth.C[2], truth.π[1], truth.π[2]

    init = SizedSIMDVector{6,Float64}(undef)
    for i ∈ 1:4
        init[i] = truth.Θ[i]
    end
    for i ∈ 7:length(truth.Θ)
        init[i-2] = truth.Θ[i]
    end
    data = NoGoldData(truth, n_common, n1, n2)
    # resample!(data, truth, n_common, n1, n2)
    ap = AsymptoticPosterior(undef, data, init)
    data, ap, init
end
function simulation(truth, k, iter, n_common, n1, n2)
    data, ap, init = set_up_sim(truth, n_common, n1, n2)
    simulation_core!(data, ap, init, truth, k, iter, n_common, n1, n2)
end

function simulation_core!(data, ap, init, truth, k, iter, n_common = 100, n1 = 1000, n2 = 30)

    sim_results = Array{Float64}(undef,3,6,k)
    S₁, S₂, C₁, C₂, π₁, π₂ = truth.S[1], truth.S[2], truth.C[1], truth.C[2], truth.π[1], truth.π[2]
    AsymptoticPosteriors.fit!(ap.map, init)
    set_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, 1)
    set_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, 1)
    set_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, 1)
    set_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, 1)
    set_buffer!(sim_results, inv_logit,       π₁, ap, 5, 1)
    set_buffer!(sim_results, inv_logit,       π₂, ap, 6, 1)
    @inbounds for j ∈ 2:iter
        resample!(data, truth, n_common, n1, n2)
        AsymptoticPosteriors.fit!(ap.map, init)
        update_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, 1)
        update_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, 1)
        update_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, 1)
        update_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, 1)
        update_buffer!(sim_results, inv_logit,       π₁, ap, 5, 1)
        update_buffer!(sim_results, inv_logit,       π₂, ap, 6, 1)
    end
    for i ∈ 2:k
        n_common *= 2
        n1 *= 2
        n2 *= 2
        resample!(data, truth, n_common, n1, n2)
        AsymptoticPosteriors.fit!(ap.map, init)
        set_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, i)
        set_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, i)
        set_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, i)
        set_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, i)
        set_buffer!(sim_results, inv_logit,       π₁, ap, 5, i)
        set_buffer!(sim_results, inv_logit,       π₂, ap, 6, i)
        @inbounds for j ∈ 2:iter
            resample!(data, truth, n_common, n1, n2)
            AsymptoticPosteriors.fit!(ap.map, init)
            update_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, i)
            update_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, i)
            update_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, i)
            update_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, i)
            update_buffer!(sim_results, inv_logit,       π₁, ap, 5, i)
            update_buffer!(sim_results, inv_logit,       π₂, ap, 6, i)
        end
    end
    sim_results
end
function simulation(truth, K1,K2,K3, iter, n_common, n1, n2)
    data, ap, init, sim_results = set_up_sim(truth, n_common, n1, n2)
    simulation_core!(data, ap, init, sim_results, truth, K1,K2,K3, iter, n_common, n1, n2)
end

function simulation_core!(data, ap, init, truth, K1,K2,K3, iter, n_common_base, n1_base, n2_base)

    sim_results = Array{Float64}(undef,3,6,K1,K2,K3)
    S₁, S₂, C₁, C₂, π₁, π₂ = truth.S[1], truth.S[2], truth.C[1], truth.C[2], truth.π[1], truth.π[2]
    for k3 ∈ 1:K3, k2 ∈ 1:K2, k1 ∈ 1:K1
        n_common = round(Int, n_common_base * 2^((k1-1)//8))
        n1 = round(Int, n1_base * 2^((k2-1)//8))
        n2 = round(Int, n2_base * 2^((k3-1)//8))
        resample!(data, truth, n_common, n1, n2)
        # try
            AsymptoticPosteriors.fit!(ap.map, init)
            set_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, k1,k2,k3)
            set_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, k1,k2,k3)
            set_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, k1,k2,k3)
            set_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, k1,k2,k3)
            set_buffer!(sim_results, inv_logit,       π₁, ap, 5, k1,k2,k3)
            set_buffer!(sim_results, inv_logit,       π₂, ap, 6, k1,k2,k3)
        # catch err
        #     @show data
        #     @show n_common, n1, n2
        #     rethrow(err)
        # end
        @inbounds for j ∈ 2:iter
            resample!(data, truth, n_common, n1, n2)
            # try
                AsymptoticPosteriors.fit!(ap.map, init)
                update_buffer!(sim_results, inv_uhalf_logit, S₁, ap, 1, k1,k2,k3)
                update_buffer!(sim_results, inv_uhalf_logit, S₂, ap, 2, k1,k2,k3)
                update_buffer!(sim_results, inv_uhalf_logit, C₁, ap, 3, k1,k2,k3)
                update_buffer!(sim_results, inv_uhalf_logit, C₂, ap, 4, k1,k2,k3)
                update_buffer!(sim_results, inv_logit,       π₁, ap, 5, k1,k2,k3)
                update_buffer!(sim_results, inv_logit,       π₂, ap, 6, k1,k2,k3)
            # catch err
            #     @show data
            #     @show n_common, n1, n2
            #     rethrow(err)
            # end
        end
    end
    sim_results
end

function update_buffer!(buffer::AbstractArray{T}, f, true_val, ap, j, k...) where T
    l = f(quantile(ap, 0.025, j))
    isfinite(l) || throw("$l not finite, j = $j")
    u = f(quantile(ap, 0.975, j))
    isfinite(u) || throw("$u not finite, j = $j")
    if (l < true_val) && (u > true_val)
        buffer[1,j,k...] += one(T)
    end
    buffer[2,j,k...] += u - l
    buffer[3,j,k...] += abs(true_val- f(ap.map.θhat[j]))
end
function set_buffer!(buffer::AbstractArray{T}, f, true_val, ap, j, k...) where T
    l = f(quantile(ap, 0.025, j))
    isfinite(l) || throw("$l not finite, j = $j")
    u = f(quantile(ap, 0.975, j))
    isfinite(u) || throw("$u not finite, j = $j")
    buffer[1,j,k...] = ifelse((l < true_val) && (u > true_val), one(T), zero(T)) # need to set it to zero.
    buffer[2,j,k...] = u - l
    buffer[3,j,k...] = abs(true_val- f(ap.map.θhat[j]))
end

@inline inv_uhalf_logit(x) = 0.5+0.5/(1+exp(-x ))

const increments = (
    ((4,0),  (-4,0), (0, 0) ),
    ((4,0),  (0,0),  (-4, 0)),
    ((-4,0), (8,0),  (0, 0) ),
    ((-4,0), (0,0),  (8, 0) ),
    ((0,0),  (4,0),  (0, 0) ),
    ((0,0),  (0,0),  (10, 0)),
    ((2,0),  (0,0),  (0, 0) ),
    ((0,4),  (0,-4), (0, 0) ),
    ((0,4),  (0,0),  (0,-4) ),
    ((0,-4), (0,8),  (0, 0) ),
    ((0,-4), (0,0),  (0,8)  ),
    ((0,0),  (0,4),  (0, 0) ),
    ((0,0),  (0,0),  (0, 4) ),
    ((0,2),  (0,0),  (0, 0) )
)

function increment_n(n_common, n1, n2, i)
    increment = increments[i]
    nctemp = n_common .+ increment[1]
    n1temp = n1 .+ increment[2]
    n2temp = n2 .+ increment[3]
    nctemp, n1temp, n2temp
end
function cost((n1cost,n2cost), n_common, n1, n2)
    n1cost * sum(n_common .+ n1) + n2cost * sum(n_common .+ n2)
end
function increment_cost((n1cost,n2cost), i)
    n_common, n1, n2 = increments[i]
    n1cost * sum(n_common .+ n1) + n2cost * sum(n_common .+ n2)
end

function average_loss(costs, iter, n_common, n1, n2)
    sim_results = Array{Float64}(undef,length(increments))
    S₁, S₂, C₁, C₂, π₁, π₂ = truth.S[1], truth.S[2], truth.C[1], truth.C[2], truth.π[1], truth.π[2]

    init = SizedSIMDVector{6,Float64}(undef)
    for i ∈ 1:4
        init[i] = truth.Θ[i]
    end
    for i ∈ 7:length(truth.Θ)
        init[i-2] = truth.Θ[i]
    end
    # uses global data and ap
    for i ∈ 1:length(increments)
        nctemp, n1temp, n2temp = increment_n(n_common, n1, n2, i)
        if any(nctemp .< 0) || any(n1temp .< 0) || any(n2temp .< 0)
            sim_results[i] = Inf
            continue
        end
        sim_results[i] = increment_cost(costs, i)*iter
        @inbounds for k ∈ 1:iter
            resample!(data, truth, nctemp, n1temp, n2temp)
            AsymptoticPosteriors.fit!(ap.map, init)
            sim_results[i] += loss(ap, init)
        end
    end
    sim_results# ./ iter
end

# @inline lossv2(v, truth, p) = v > truth ? p*(v-truth) : (1-p)*(truth-v)
@inline loss(v, truth, p) = (p - (v < truth))*(v-truth)
loss((l, u)::Tuple{<:Number,<:Number}, truth, pl = 0.025, pu = 0.975) = loss(l, truth, pl) + loss(u, truth, pu)

function loss(ap::AsymptoticPosteriors.AsymptoticPosterior, truth)
    l = 0.0
    @inbounds for i ∈ eachindex(truth)
        l += loss((quantile(ap, 0.025, i),quantile(ap, 0.975, i)), truth[i])
    end
    l
end

end #everywhere

function run_simulation(truth, k = 5, iter = 25, n_common = 200, n1 = 100, n2 = 30)
    @distributed (+) for i ∈ 2:nprocs()
        simulation(truth, k, iter, n_common, n1, n2)
    end
end
function run_simulation2(truth, k = 10, iter = 50, n_common = 100, n1 = 500, n2 = 30)
    @distributed (+) for i ∈ 2:nprocs()
        simulation(truth, k,k,k, iter, n_common, n1, n2)
    end
end

function run_simulation(truth, N=10, k = 5, iter = 25, n_common = 200, n1 = 100, n2 = 30)
    out = simulation_core(truth, k, iter, n_common, n1, n2)
    for i ∈ 2:N
        out += simulation(truth, k, iter, n_common, n1, n2)
    end
    out
end

# @time @everywhere begin
# const data = NoGoldData(truth)
# const ap = AsymptoticPosterior(undef, data, x)
# end #everywhere

function simulation_search( costs, num_iter, n_common = (100,100), n1 = (200,200), n2 = (30,30) )
    nproc = nprocs()
    current_loss = Inf
    losses = @distributed (+) for i ∈ 2:nproc
        average_loss(costs, num_iter, n_common, n1, n2)
    end
    last_loss, minind = findmin(losses)
    iteration = 1
    while last_loss < current_loss
        @show iteration, losses
        @show n_common, n1, n2
        iteration += 1
        current_loss = last_loss
        n_common, n1, n2 = increment_n(n_common, n1, n2, minind)
        losses = @distributed (+) for i ∈ 2:nproc
            average_loss(costs, num_iter, n_common, n1, n2)
        end
        last_loss, minind = findmin(losses)
        if last_loss > current_loss # We'll try again.
            println("Failed to find best argument. Searching broadly.")
            for n ∈ 1:length(increments)
                n_common_temp, n1_temp, n2_temp = increment_n(n_common, n1, n2, n)
                (any(n_common_temp .< 0) || any(n1_temp .< 0) || any(n2_temp .< 0)) && continue
                losses_temp = @distributed (+) for i ∈ 2:nproc
                    average_loss(costs, num_iter, n_common_temp, n1_temp, n2_temp)
                end
                @show n_common_temp, n1_temp, n2_temp
                @show losses_temp
                last_loss_temp, minind_temp = findmin(losses_temp)
                last_loss_temp += increment_cost(costs, n) * num_iter * (nproc-1)
                if last_loss_temp < current_loss
                    last_loss = last_loss_temp
                    n_common, n1, n2 = n_common_temp, n1_temp, n2_temp
                    minind = minind_temp
                    break
                end
            end

        end
    end
    @show iteration, losses
    n_common, n1, n2, current_loss
end

# @time simulation_search((0.0000005,0.005), 1000, (150,150), (5,5), (5,5))
# 4869.788621 seconds (380.42 k allocations: 35.323 MiB, 0.00% gc time)
# ((350, 320), (75, 15), (5, 5), 441712.07542121154)

# 138 137
# 274 235
# 898 275
# 942 283
# 964* 900
# 968* 944
# 971* 966*
# 973* 970*
# 1012 973*
# 1014* 975*
# 1159 1014
# 1295 1016*
# 1335 1161
#      1297
#      1337
#


text:
using Documenter, PowerSystems
using InfrastructureSystems
const PSYPATH = dirname(pathof(PowerSystems))

# This is commented out because the output is not user-friendly. Deliberation on how to best
# communicate this information to users is ongoing.
#include(joinpath(@__DIR__, "src", "generate_validation_table.jl"))

makedocs(
    modules = [PowerSystems],
    format = Documenter.HTML(mathengine = Documenter.MathJax()),
    sitename = "PowerSystems.jl",
    authors = "Jose Daniel Lara, Daniel Thom and Clayton Barrows",
    pages = Any[ # Compat: `Any` for 0.4 compat
        "Introduction" => "index.md",
        "User Guide" => Any[
            "man/parsing.md",
            "man/data.md",
        ],
        "Developer" => Any[
            "Tests" => "developer/tests.md",
            "Logging" => "developer/logging.md",
            "Style Guide" => "developer/style.md",
            "Extending Parsing" => "developer/extending_parsing.md",
        ],
        "API" => Any[
            "PowerSystems" => "api/PowerSystems.md"
        ]
    ]
)

deploydocs(
    repo = "github.com/NREL/PowerSystems.jl.git",
    branch = "gh-pages",
    target = "build",
    make = nothing,
)

ext: jl


text:
using QueryOperators
using DataValues
using Test

@testset "QueryOperators" begin

@testset "Core" begin

source_1 = [1,2,2,3,4]
enum = QueryOperators.query(source_1)

@test collect(QueryOperators.@filter(QueryOperators.query(source_1), i->i>2)) == [3,4]

@test collect(QueryOperators.@map(QueryOperators.query(source_1), i->i^2)) == [1,4,4,9,16]

@test collect(QueryOperators.@take(enum, 2)) == [1,2]
@test collect(QueryOperators.@drop(enum, 2)) == [2,3,4]
@test QueryOperators.@count(enum) == 5
@test QueryOperators.@count(enum, x->x%2==0) == 3

dropped_str = ""
for i in QueryOperators.drop(enum, 2)
    dropped_str *= string(i)
end
@test dropped_str == "234"

dropped_str = ""
for i in QueryOperators.drop(enum, 80)
    dropped_str *= string(i)
end
@test dropped_str == ""

taken_str = ""
for i in QueryOperators.take(enum, 2)
    taken_str *= string(i)
end
@test taken_str == "12"

filtered_str = ""
for i in QueryOperators.@filter(enum, x->x%2==0)
    filtered_str *= string(i)
end
@test filtered_str == "224"

filtered_str = ""
for i in QueryOperators.@filter(enum, x->x>100)
    filtered_str *= string(i)
end
@test filtered_str == ""

@test collect(QueryOperators.@filter(enum, x->x<3)) == [1,2,2]

grouped = []
for i in QueryOperators.@groupby(QueryOperators.query(source_1), i->i, i->i^2)
    push!(grouped, i)
end

@test grouped == [[1],[4,4],[9],[16]]

mapped = []
for i in collect(QueryOperators.@map(enum, i->i*3))
    push!(mapped, i)
end
@test mapped == [3,6,6,9,12]

# ensure that the default value must be of the same type
errored = false
try
    QueryOperators.@default_if_empty(source_1, "string")
catch
    errored = true
end

@test errored == true

# default_if_empty for regular array
d = []
for i in QueryOperators.@default_if_empty(source_1, 0)
    push!(d, i)
end
@test d == [1, 2, 2, 3, 4]

@test collect(QueryOperators.default_if_empty(DataValue{Int}[]))[1] == DataValue{Int}()
@test collect(QueryOperators.default_if_empty(DataValue{Int}[], DataValue{Int}()))[1] == DataValue{Int}()

# passing in a NamedTuple of DataValues
nt = (a=DataValue(2), b=DataValue("test"), c=DataValue(3))
def = QueryOperators.default_if_empty(typeof(nt)[])
@test typeof(collect(def)[1]) == typeof(nt)

ordered = QueryOperators.@orderby(enum, x -> -x)
@test collect(ordered) == [4, 3, 2, 2, 1]

filtered = QueryOperators.@orderby(QueryOperators.@filter(enum, x->x%2 == 0), x->x)
@test collect(filtered) == [2, 2, 4]

ordered = QueryOperators.@orderby_descending(enum, x -> -x)
@test collect(ordered) == [1, 2, 2, 3, 4]

desired = [[1], [2, 2, 3], [4]]
grouped = QueryOperators.@groupby(enum, x -> floor(x/2), x->x)
@test collect(grouped) == desired

group_no_macro = QueryOperators.groupby(enum, x -> floor(x/2), quote x->floor(x/2) end)
@test collect(group_no_macro) == desired

outer = QueryOperators.query([1,2,3,4,5,6])
inner = QueryOperators.query([2,3,4,5])

join_desired = [[3,2], [4,3], [5,4], [6,5]]
@test collect(QueryOperators.@join(outer, inner, x->x, x->x+1, (i,j)->[i,j])) == join_desired

group_desired = [[1, Int64[]], [2, Int64[]], [3, [2]], [4, [3]], [5, [4]], [6, [5]]]
@test collect(QueryOperators.@groupjoin(outer, inner, x->x, x->x+1, (i,j)->[i,j])) == group_desired

many_map_desired =  [[1, 2], [2, 4], [2, 4], [3, 6], [4, 8]]
success = collect(QueryOperators.@mapmany(enum, x->[x*2], (x,y)->[x,y])) == many_map_desired
@test success       # for some reason, this is required to avoid a BoundsError

first = QueryOperators.query([1, 2])
second = [3, 4]
many_map_desired = [(1,3), (1,4), (2,3), (2,4)]
success = collect(QueryOperators.@mapmany(first, i->second, (x,y)->(x,y))) == many_map_desired
@test success

ntups = QueryOperators.query([(a=1, b=2, c=3), (a=4, b=5, c=6)])

@test sprint(show, ntups) == """
2x3 query result
a │ b │ c
──┼───┼──
1 │ 2 │ 3
4 │ 5 │ 6"""

@test sprint(show, enum) == """
5-element query result
 1
 2
 2
 3
 4"""

@test sprint((stream,data)->show(stream, "text/html", data), ntups) ==
    "

a	b	c

1	2	3

4	5	6

"

gather_result1 = QueryOperators.gather(QueryOperators.query([(US=1, EU=1, CN=1), (US=2, EU=2, CN=2), (US=3, EU=3, CN=3)]))
@test sprint(show, gather_result1) == """9x2 query result\nkey │ value\n────┼──────\n:US │ 1    \n:EU │ 1    \n:CN │ 1    \n:US │ 2    \n:EU │ 2    \n:CN │ 2    \n:US │ 3    \n:EU │ 3    \n:CN │ 3    """
gather_result2 = QueryOperators.gather(QueryOperators.query([(Year=2017, US=1, EU=1, CN=1), (Year=2018, US=2, EU=2, CN=2), (Year=2019, US=3, EU=3, CN=3)]), :US, :EU, :CN)
@test sprint(show, gather_result2) == "9x3 query result\nYear │ key │ value\n─────┼─────┼──────\n2017 │ :US │ 1    \n2017 │ :EU │ 1    \n2017 │ :CN │ 1    \n2018 │ :US │ 2    \n2018 │ :EU │ 2    \n2018 │ :CN │ 2    \n2019 │ :US │ 3    \n2019 │ :EU │ 3    \n2019 │ :CN │ 3    "

@test sprint((stream,data)->show(stream, "application/vnd.dataresource+json", data), ntups) ==
    "{\"schema\":{\"fields\":[{\"name\":\"a\",\"type\":\"integer\"},{\"name\":\"b\",\"type\":\"integer\"},{\"name\":\"c\",\"type\":\"integer\"}]},\"data\":[{\"a\":1,\"b\":2,\"c\":3},{\"a\":4,\"b\":5,\"c\":6}]}"

end

include("test_enumerable_unique.jl")

include("test_namedtupleutilities.jl")

end


text:
"""
# Refs
[1] F. L. Lewis, D. Vrabie, and K. Vamvoudakis,	"Reinforcement Learning and Feedback Control: Using Natural Decision Mehods to Design Optimal Adaptive Controllers," IEEE Control Systems, vol. 32, no. 6, pp.76-105, 2012.
# Variables
V̂ ∈ R: the estimate of (state) value function
    - V̂.basis: Φ
    - V̂.param: w
d: polynomial degree
"""
mutable struct CTLinearValueIterationIRL
    Q
    R
    V̂::LinearApproximator
    T::Real
    N::Int
    function CTLinearValueIterationIRL(Q, R,
                                       T=0.04, N=3, d_value::Int=2, d_controller::Int=4;
                                       V̂=LinearApproximator(2, 2; with_bias=false),
                                      )
        @assert T > 0
        new(Q, R, V̂, T, N)
    end
end

function RunningCost(irl::CTLinearValueIterationIRL)
    @unpack Q, R = irl
    return QuadraticCost(Q, R)
end

"""
Infer he approximate optimal input.
"""
function ApproximateOptimalInput(irl::CTLinearValueIterationIRL, B)
    @unpack R = irl
    return function (x, p, t)
        P = ARE_solution(irl)
        K = inv(R) * B' * P
        û = -K * x
    end
end

function ARE_solution(irl::CTLinearValueIterationIRL)
    ŵ = irl.V̂.param
    P = [  ŵ[1] ŵ[2]/2;
         ŵ[2]/2   ŵ[3]]
    P
end

function update_params_callback(irl::CTLinearValueIterationIRL, w_tol)
    stop_conds = function(w, w_prev)
        stop_conds_dict = Dict(
                               :w_tol => norm(w - w_prev) < w_tol,
                              )
    end
    i = 0
    w_prev = nothing
    ∫rs = []
    V̂_nexts = []
    Φs = []
    V̂_nexts = []
    is_terminated = false
    affect! = function (integrator)
        @unpack t = integrator
        X = integrator.u
        @unpack x, ∫r = X
        ∫r = length(∫r) == 1 ? ∫r[1] : error("Invalid ∫r")  # for both Number and Array
        push!(∫rs, ∫r)
        push!(Φs, irl.V̂.basis(x))
        if length(∫rs) > 1
            push!(V̂_nexts, diff(∫rs[end-1:end])[1] + irl.V̂(x)[1])  # initial_affect=true
        end
        if is_terminated == false
            if length(V̂_nexts) >= irl.N
                i += 1
                irl.V̂.param = pinv(hcat(Φs[end-irl.N:end-1]...)') * hcat(V̂_nexts...)'
                if w_prev != nothing
                    stop_conds_dict = stop_conds(irl.V̂.param, w_prev)
                    is_terminated = stop_conds_dict |> values |> any
                end
                w_prev = deepcopy(irl.V̂.param)
                V̂_nexts = []
                @show t, i, irl.V̂.param
            end
        end
    end
    cb_train = PeriodicCallback(affect!, irl.T; initial_affect=true)
end


text:
module Types
    abstract type AbstractSpec end
    abstract type SpecModifier <: AbstractSpec  end

    struct SingleComponent <: AbstractSpec end
    struct OneMol <: AbstractSpec end

    abstract type MassBasisModifier <: SpecModifier end
    struct MOLAR <: MassBasisModifier end
    struct MASS <: MassBasisModifier end
    struct TOTAL <: MassBasisModifier end

    abstract type CompoundModifier <: SpecModifier end
    struct FRACTION <: CompoundModifier end
    struct TOTAL_AMOUNT <: CompoundModifier end

    abstract type VolumeModifier <: SpecModifier end
    struct DENSITY <: VolumeModifier end
    struct VOLUME <: VolumeModifier end

    abstract type AbstractIntensiveSpec{T1<:MassBasisModifier} <: AbstractSpec  end
    abstract type AbstractTotalSpec <: AbstractSpec  end
    abstract type AbstractFractionSpec <: AbstractSpec  end
    abstract type CategoricalSpec <: AbstractSpec  end

    abstract type AbstractEnergySpec{T1} <: AbstractIntensiveSpec{T1} end

    struct Enthalpy{T} <: AbstractEnergySpec{T} end
    struct InternalEnergy{T} <: AbstractEnergySpec{T} end
    struct Gibbs{T} <: AbstractEnergySpec{T} end
    struct Helmholtz{T} <: AbstractEnergySpec{T} end

    struct Entropy{T} <: AbstractIntensiveSpec{T} end

    struct VolumeAmount{T1,T2<:VolumeModifier} <: AbstractIntensiveSpec{T1} end

    struct Pressure <: AbstractSpec end
    struct Temperature <: AbstractSpec end

    # those are vectors,

    struct MaterialCompounds{T1<:MassBasisModifier,T2<:CompoundModifier} <: AbstractTotalSpec end

    struct MaterialAmount{T1<:MassBasisModifier} <: AbstractTotalSpec end

    #humidity, to support wet gas
    abstract type HumidityModifier <: AbstractSpec end
    struct WetBulbTemperature <:HumidityModifier end
    struct RelativeHumidity <:HumidityModifier end
    struct HumidityRatio <:HumidityModifier end
    struct MassHumidity  <:HumidityModifier end
    struct MolarHumidity <:HumidityModifier end
    struct HumidityDewPoint <:HumidityModifier end

    struct HumiditySpec{T<:HumidityModifier} <: AbstractSpec end

    struct PhaseFractions <: AbstractFractionSpec end
    struct VolumeFraction <: AbstractFractionSpec end
    struct VaporFraction{T<:MassBasisModifier} <: AbstractFractionSpec end

    struct PhaseTag <: CategoricalSpec end
    struct TwoPhaseEquilibrium <: CategoricalSpec end

    struct Options <: CategoricalSpec end

    #not defined for now
    struct MolecularWeight <: AbstractSpec end

    #variable spec, signal to make a variable specs
    struct VariableSpec  end

    #base type to define models
    abstract type ThermoModel end

    export AbstractSpec
    export SpecModifier
    export AbstractIntensiveSpec
    export AbstractTotalSpec
    export AbstractFractionSpec
    export CategoricalSpec
    export SingleComponent
    export OneMol
    export MOLAR
    export MASS
    export TOTAL
    export FRACTION
    export TOTAL_AMOUNT
    export DENSITY
    export VOLUME
    export AbstractEnergySpec
    export Enthalpy
    export InternalEnergy
    export Gibbs
    export Helmholtz
    export Entropy
    export VolumeAmount
    export Pressure
    export Temperature
    export MaterialCompounds
    export MaterialAmount
    export PhaseFractions
    export VaporFraction
    export PhaseTag
    export TwoPhaseEquilibrium
    export Options
    export MolecularWeight
    export VariableSpec
    export ThermoModel

    export HumidityModifier
    export WetBulbTemperature
    export RelativeHumidity
    export HumidityRatio
    export MassHumidity
    export MolarHumidity
    export HumidityDewPoint

    export HumiditySpec

end

using .Types

const KW_TO_SPEC = Dict{Symbol,Any}(
:h =>  Enthalpy{MOLAR}()
,:g =>  Gibbs{MOLAR}()
,:a =>  Helmholtz{MOLAR}()
,:u =>  InternalEnergy{MOLAR}()

,:mol_h =>  Enthalpy{MOLAR}()
,:mol_g =>  Gibbs{MOLAR}()
,:mol_a =>  Helmholtz{MOLAR}()
,:mol_u =>  InternalEnergy{MOLAR}()

,:mass_h =>  Enthalpy{MASS}()
,:mass_g =>  Gibbs{MASS}()
,:mass_a =>  Helmholtz{MASS}()
,:mass_u =>  InternalEnergy{MASS}()

,:total_h =>  Enthalpy{TOTAL}()
,:total_g =>  Gibbs{TOTAL}()
,:total_a =>  Helmholtz{TOTAL}()
,:total_u =>  InternalEnergy{TOTAL}()

,:s =>  Entropy{MOLAR}()
,:mol_s =>  Entropy{MOLAR}()
,:mass_s =>  Entropy{MASS}()
,:total_s =>  Entropy{TOTAL}()

,:p =>  Pressure()
,:P =>  Pressure()
,:t =>  Temperature()
,:T => Temperature()

,:v =>  VolumeAmount{MOLAR,VOLUME}()
,:mol_v =>  VolumeAmount{MOLAR,VOLUME}()
,:mass_v =>  VolumeAmount{MASS,VOLUME}()
,:total_v =>  VolumeAmount{TOTAL,VOLUME}()
,:V =>  VolumeAmount{TOTAL,VOLUME}()

,:rho =>  VolumeAmount{MOLAR,DENSITY}()
,:mol_rho =>  VolumeAmount{MOLAR,DENSITY}()
,:mass_rho =>  VolumeAmount{MASS,DENSITY}()

,:ρ =>  VolumeAmount{MOLAR,DENSITY}()
,:mol_ρ =>  VolumeAmount{MOLAR,DENSITY}()
,:mass_ρ =>  VolumeAmount{MASS,DENSITY}()
,:mass =>  MaterialAmount{MASS}()
,:moles =>  MaterialAmount{MOLAR}()
,:xn =>  MaterialCompounds{MOLAR,FRACTION}()
,:xm =>   MaterialCompounds{MASS,FRACTION}()
,:n =>  MaterialCompounds{MOLAR,TOTAL_AMOUNT}()
,:m =>   MaterialCompounds{MASS,TOTAL_AMOUNT}()

,:mw =>  MolecularWeight()

,:quality =>  VaporFraction{MASS}() #looking for better name
,:vfrac =>  VaporFraction{MOLAR}() #looking for better name

,:phase_fracs =>  PhaseFractions() #looking for better name

,:phase =>PhaseTag()

,:sat => TwoPhaseEquilibrium()

,:single_component => SingleComponent()
,:one_mol => OneMol()
,:options => Options()

,:hum_wetbulb => HumiditySpec{WetBulbTemperature}()
,:hum_ratio => HumiditySpec{HumidityRatio}()
,:hum_molfrac => HumiditySpec{MolarHumidity}()
,:hum_massfrac => HumiditySpec{MassHumidity}()
,:rel_hum => HumiditySpec{RelativeHumidity}()
,:hum_dewpoint => HumiditySpec{HumidityDewPoint}()

)

const SPEC_TO_KW = Dict{Any,Symbol}(
Enthalpy{MOLAR}() => :mol_h
,Gibbs{MOLAR}() => :mol_g
,Helmholtz{MOLAR}() => :mol_a
,InternalEnergy{MOLAR}() => :mol_u

,Enthalpy{MASS}() => :mass_h
,Gibbs{MASS}()  => :mass_g
,Helmholtz{MASS}() => :mass_a
,InternalEnergy{MASS}() => :mass_u

,Enthalpy{TOTAL}() => :total_h
,Gibbs{TOTAL}() => :total_g
,Helmholtz{TOTAL}() => :total_a
,InternalEnergy{TOTAL}() => :total_u

,Entropy{MOLAR}() => :mol_s
,Entropy{MASS}() => :mass_s
,Entropy{TOTAL}() => :total_s

,Pressure() => :p
,Temperature() => :t

,VolumeAmount{MOLAR,VOLUME}() => :mol_v
,VolumeAmount{MASS,VOLUME}() => :mass_v
,VolumeAmount{TOTAL,VOLUME}() => :total_v

,VolumeAmount{MOLAR,DENSITY}()  => :rho
,VolumeAmount{MOLAR,DENSITY}()  => :mol_rho
,VolumeAmount{MASS,DENSITY}()  => :mass_rho

,MaterialAmount{MASS}() => :mass
,MaterialAmount{MOLAR}() => :moles
,MaterialCompounds{MOLAR,FRACTION}() => :xn
,MaterialCompounds{MASS,FRACTION}() => :xm
,MaterialCompounds{MOLAR,TOTAL_AMOUNT}() => :n
,MaterialCompounds{MASS,TOTAL_AMOUNT}() =>  :n

,MolecularWeight() => :mw

,VaporFraction{MASS}() => :quality#looking for better name
,VaporFraction{MOLAR}() => :vfrac#looking for better name

,PhaseFractions() => :phase_fracs #looking for better name

,PhaseTag() => :phase

,TwoPhaseEquilibrium() => :sat

,SingleComponent() => :single_component
,OneMol() => :one_mol
,Options() => :options

,HumiditySpec{WetBulbTemperature}()  => :hum_wetbulb
,HumiditySpec{HumidityRatio}() => :hum_ratio
,HumiditySpec{MolarHumidity}() => :hum_molfrac
,HumiditySpec{MassHumidity}() => :hum_massfrac
,HumiditySpec{RelativeHumidity}() => :rel_hum
,HumiditySpec{HumidityDewPoint}() => :hum_dewpoint

)

module QuickStates
using ..Types
const SinglePT = Tuple{Pressure,Temperature,SingleComponent}
const MultiPT = Tuple{Pressure,Temperature,MaterialCompounds}

pt() = (Pressure(),Temperature(),SingleComponent())
ptx() = (Pressure(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
ptn() = (Pressure(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SinglePV = Tuple{Pressure,VolumeAmount,SingleComponent}
const MultiPV = Tuple{Pressure,VolumeAmount,MaterialCompounds}

pv() = (Pressure(),VolumeAmount{MOLAR,VOLUME}(),SingleComponent())
pvx() = (Pressure(),VolumeAmount{MOLAR,VOLUME}(),MaterialCompounds{MOLAR,FRACTION}())
pvn() = (Pressure(),VolumeAmount{MOLAR,VOLUME}(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SingleVT = Tuple{VolumeAmount,Temperature,SingleComponent}
const MultiVT = Tuple{VolumeAmount,Temperature,MaterialCompounds}

vt() = (VolumeAmount{MOLAR,VOLUME}(),Temperature(),SingleComponent())
vtx() = (VolumeAmount{MOLAR,VOLUME}(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
vtn() = (VolumeAmount{MOLAR,VOLUME}(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

ρt() = (VolumeAmount{MOLAR,DENSITY}(),Temperature(),SingleComponent())
ρtx() = (VolumeAmount{MOLAR,DENSITY}(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
ρtn() = (VolumeAmount{MOLAR,VOLUME}(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SinglePS = Tuple{Pressure,Entropy,SingleComponent}
const MultiPS = Tuple{Pressure,Entropy,MaterialCompounds}

st() = (Pressure(),Entropy{MOLAR}(),SingleComponent())
stx() = (Pressure(),Entropy{MOLAR}(),MaterialCompounds{MOLAR,FRACTION}())
stn() = (Pressure(),Entropy{MOLAR}(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SinglePH = Tuple{Pressure,Enthalpy,SingleComponent}
const MultiPH = Tuple{Pressure,Enthalpy,MaterialCompounds}

ph() = (Pressure(),Enthalpy{MOLAR}(),SingleComponent())
phx() = (Pressure(),Enthalpy{MOLAR}(),MaterialCompounds{MOLAR,FRACTION}())
phn() = (Pressure(),Enthalpy{MOLAR}(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SingleSatT = Tuple{TwoPhaseEquilibrium,Temperature,SingleComponent}
const MultiSatT = Tuple{TwoPhaseEquilibrium,Temperature,MaterialCompounds}

sat_t() = (TwoPhaseEquilibrium(),Temperature(),SingleComponent())
sat_tx() = (TwoPhaseEquilibrium(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
sat_tn() = (TwoPhaseEquilibrium(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SingleSatP = Tuple{TwoPhaseEquilibrium,Pressure,SingleComponent}
const MultiSatP = Tuple{TwoPhaseEquilibrium,Pressure,MaterialCompounds}

sat_p() = (TwoPhaseEquilibrium(),Pressure(),SingleComponent())
sat_px() = (TwoPhaseEquilibrium(),Pressure(),MaterialCompounds{MOLAR,FRACTION}())
sat_pn() = (TwoPhaseEquilibrium(),Pressure(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

molar_ϕt() = (VaporFraction{MOLAR}(),Temperature(),SingleComponent())
molar_ϕntx() = (VaporFraction{MOLAR}(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
molar_ϕtn() = (VaporFraction{MOLAR}(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

mass_ϕt() = (VaporFraction{MASS}(),Temperature(),SingleComponent())
mass_ϕntx() = (VaporFraction{MASS}(),Temperature(),MaterialCompounds{MOLAR,FRACTION}())
mass_ϕtn() = (VaporFraction{MASS}(),Temperature(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const SingleΦP = Tuple{VaporFraction,Pressure,SingleComponent}
const MultiΦP = Tuple{VaporFraction,Pressure,MaterialCompounds}
const SingleΦT = Tuple{VaporFraction,Temperature,SingleComponent}
const MultiΦT = Tuple{VaporFraction,Temperature,MaterialCompounds}

const SingleΦnP = Tuple{VaporFraction{MOLAR},Pressure,SingleComponent}
const MultiΦnP = Tuple{VaporFraction{MOLAR},Pressure,MaterialCompounds}
const SingleΦnT = Tuple{VaporFraction{MOLAR},Temperature,SingleComponent}
const MultiΦnT = Tuple{VaporFraction{MOLAR},Temperature,MaterialCompounds}

const SingleΦmP = Tuple{VaporFraction{MASS},Pressure,SingleComponent}
const MultiΦmP = Tuple{VaporFraction{MASS},Pressure,MaterialCompounds}
const SingleΦmT = Tuple{VaporFraction{MASS},Temperature,SingleComponent}
const MultiΦmT = Tuple{VaporFraction{MASS},Temperature,MaterialCompounds}

ϕp() = (VaporQuality(),Pressure(),SingleComponent())
ϕpx() = (VaporQuality(),Pressure(),MaterialCompounds{MOLAR,FRACTION}())
ϕpn() = (VaporQuality(),Pressure(),MaterialCompounds{MOLAR,TOTAL_AMOUNT}())

const MultiT = Union{Tuple{Temperature,AbstractSpec,MaterialCompounds},Tuple{AbstractSpec,Temperature,MaterialCompounds}}
const SingleT = Union{Tuple{Temperature,AbstractSpec,SingleComponent},Tuple{AbstractSpec,Temperature,SingleComponent}}

export SinglePT,MultiPT
export SingleVT,MultiVT
export SinglePS,MultiPS
export SinglePH,MultiPH
export SingleSatT,MultiSatT
export SingleSatP,MultiSatP
export SingleΦT,MultiΦT
export SingleΦP,MultiΦP
export SingleΦmT,MultiΦmT
export SingleΦmP,MultiΦmP
export SingleΦnT,MultiΦnT
export SingleΦnP,MultiΦnP
export SinglePV,MultiPV
export SingleT,MultiT

end


text:
"""
    Strings = Union{Str,String,StringDocument, Vector{String}}

Create a Union of three types of strings that the functions of the package accepts; Str, String, TextAnalysis::StringDocument and  Vector{String}.
"""
const Strings = Union{Str,String,StringDocument, Vector{String}, Vector{Vector{String}}}

"Create a new type that carries all relevant information needed for calculating the lexical gravity index: gdict is a feature that includes a dictionary
of the words and their lexical richness index."
struct RichnessBundle
    gdict::Dict{String, Float64}
    uniquecontext::Vector{String}
    bigrams::Vector{Vector{String}}
    word::String
end


text:
@testset "Processes" begin
  @testset "Basic" begin
    for p in [BinomialProcess(100), PoissonProcess(100.)]
      r = Box((0.,1.), (1.,2.))
      pp = rand(p, r)
      xs = [coordinates(pp[i]) for i in 1:nelements(pp)]
      @test all(0 .≤ first.(xs) .≤ 1)
      @test all(1 .≤ last.(xs) .≤ 2)
    end
  end

  @testset "Binomial" begin
    # TODO
  end

  @testset "Poisson" begin
    # TODO
  end

  @testset "Union" begin
    b  = Box((0.,0.), (100.,100.))
    p₁ = BinomialProcess(50)
    p₂ = BinomialProcess(50)
    p  = p₁ ∪ p₂ # 100 points

    s = rand(p, b, 2)
    @test length(s) == 2
    @test s[1] isa PointSet
    @test s[2] isa PointSet
    @test nelements.(s) == [100,100]
  end
end


text:
using Mimi
using CSV

## Notes:
## Applied PAGE09 weighting is non-stochastic

include("../src/main_model.jl")

for scenario in ["RCP4.5 & SSP2", "RCP8.5 & SSP5", "RCP1.9 & SSP1", "RCP2.6 & SSP1"]

    m0 = getpage(scenario, false, false, true; use_page09weights=true, page09_discontinuity=true, page09_sealevelrise=true) # just updating climate
    set_param!(m0, :NonMarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m0, :SLRDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m0, :MarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    m1 = getpage(scenario, false, false, true; use_page09weights=true, page09_discontinuity=true) # long-tailed SLR
    set_param!(m1, :NonMarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m1, :SLRDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m1, :MarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    m2 = getpage(scenario, true, false, true; use_page09weights=true, page09_discontinuity=true) # no sea-ice, yes perma
    set_param!(m2, :NonMarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m2, :SLRDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m2, :MarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    m3 = getpage(scenario, true, true, true; use_page09weights=true, page09_discontinuity=true) # all arctic
    set_param!(m3, :NonMarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m3, :SLRDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m3, :MarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    m4 = getpage(scenario, true, true, true; use_page09weights=true) # reduced discontinuity
    set_param!(m4, :NonMarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m4, :SLRDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m4, :MarketDamages, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    set_param!(m4, :Discontinuity, :GDP_per_cap_focus_0_FocusRegionEU, 27934.244777382406)
    m5 = getpage(scenario, true, true, true; use_page09weights=true)  # update GDPpc_0
    m6 = getpage(scenario, true, true, true) # update weighting
    m7 = getpage(scenario, true, true, false) # full PAGE-ICE

    samplesize = 50000
    year = 2020
    pulse_size = 75000.

    # Setup of location of final results
    scc_results0 = zeros(samplesize)
    scc_results1 = zeros(samplesize)
    scc_results2 = zeros(samplesize)
    scc_results3 = zeros(samplesize)
    scc_results4 = zeros(samplesize)
    scc_results5 = zeros(samplesize)
    scc_results6 = zeros(samplesize)
    scc_results7 = zeros(samplesize)
    scc_results8 = zeros(samplesize)
    dmg_results0 = zeros(samplesize, 8)
    dmg_results1 = zeros(samplesize, 8)
    dmg_results2 = zeros(samplesize, 8)
    dmg_results3 = zeros(samplesize, 8)
    dmg_results4 = zeros(samplesize, 8)
    dmg_results5 = zeros(samplesize, 8)
    dmg_results6 = zeros(samplesize, 8)
    dmg_results7 = zeros(samplesize, 8)
    dmg_results8 = zeros(samplesize, 8)

    function mc_scc_calculation(sim_inst::SimulationInstance, trialnum::Int, ntimesteps::Int, ignore::Nothing)
        scc_results0[trialnum] = sim_inst.models[1][:EquityWeighting, :tdac_totalimpactandadaptation_2200]
        scc_results1[trialnum] = sim_inst.models[2][:EquityWeighting, :tdac_totalimpactandadaptation_2200]
        scc_results2[trialnum] = sim_inst.models[3][:EquityWeighting, :tdac_totalimpactandadaptation_2200]
        scc_results3[trialnum] = sim_inst.models[4][:EquityWeighting, :tdac_totalimpactandadaptation_2200]
        scc_results4[trialnum] = sim_inst.models[5][:EquityWeighting, :tdac_totalimpactandadaptation_2200]
        scc_results5[trialnum] = sim_inst.models[5][:EquityWeighting, :tdac_totalimpactandadaptation]
        scc_results6[trialnum] = sim_inst.models[6][:EquityWeighting, :tdac_totalimpactandadaptation]
        scc_results7[trialnum] = sim_inst.models[7][:EquityWeighting, :tdac_totalimpactandadaptation]
        scc_results8[trialnum] = sim_inst.models[8][:EquityWeighting, :tdac_totalimpactandadaptation]

        dmg_results0[trialnum, :] = sum(sim_inst.models[1][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated][1:8, :] + sim_inst.models[1][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated][1:8, :], dims=1)
        dmg_results1[trialnum, :] = sum(sim_inst.models[2][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated][1:8, :] + sim_inst.models[2][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated][1:8, :], dims=1)
        dmg_results2[trialnum, :] = sum(sim_inst.models[3][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated][1:8, :] + sim_inst.models[3][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated][1:8, :], dims=1)
        dmg_results3[trialnum, :] = sum(sim_inst.models[4][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated][1:8, :] + sim_inst.models[4][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated][1:8, :], dims=1)
        dmg_results4[trialnum, :] = sum(sim_inst.models[5][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated][1:8, :] + sim_inst.models[5][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated][1:8, :], dims=1)
        dmg_results5[trialnum, :] = sum(sim_inst.models[5][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] + sim_inst.models[5][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated], dims=1)
        dmg_results6[trialnum, :] = sum(sim_inst.models[6][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] + sim_inst.models[6][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated], dims=1)
        dmg_results7[trialnum, :] = sum(sim_inst.models[7][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] + sim_inst.models[7][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated], dims=1)
        dmg_results8[trialnum, :] = sum(sim_inst.models[8][:EquityWeighting, :addt_equityweightedimpact_discountedaggregated] + sim_inst.models[8][:EquityWeighting, :aact_equityweightedadaptation_discountedaggregated], dims=1)
    end

    # get simulation
    mcs = getsim()

    mm0 = get_marginal_model(m0, year=year, pulse_size=pulse_size)
    mm1 = get_marginal_model(m1, year=year, pulse_size=pulse_size)
    mm2 = get_marginal_model(m2, year=year, pulse_size=pulse_size)
    mm3 = get_marginal_model(m3, year=year, pulse_size=pulse_size)
    mm4 = get_marginal_model(m4, year=year, pulse_size=pulse_size)
    mm5 = get_marginal_model(m5, year=year, pulse_size=pulse_size)
    mm6 = get_marginal_model(m6, year=year, pulse_size=pulse_size)
    mm7 = get_marginal_model(m7, year=year, pulse_size=pulse_size)

    # Run
    res = run(mcs, [mm0, mm1, mm2, mm3, mm4, mm5, mm6, mm7], samplesize; post_trial_func=mc_scc_calculation)

    results = DataFrame(page09=scc_results0, bigslr=scc_results1, permaf=scc_results2, seaice=scc_results3, redisc=scc_results4,
                        to2300=scc_results5, gdppc0=scc_results6, newwts=scc_results7, burked=scc_results8)

    if scenario == "1.5 degC Target"
        suffix = "deg1p5"
    elseif scenario == "RCP1.9 & SSP1"
        suffix = "ssp119"
    elseif scenario == "RCP2.6 & SSP1"
        suffix = "ssp126"
    elseif scenario == "RCP4.5 & SSP2"
        suffix = "ssp245"
    elseif scenario == "RCP8.5 & SSP5"
        suffix = "ssp585"
    elseif scenario == "RCP8.5 & SSP2"
        suffix = "ssp285"
    end
    CSV.write("xscc-damages-$suffix.csv", results)

    using Tables
    CSV.write("page09-xreg-$suffix.csv", Tables.table(dmg_results0; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("bigslr-xreg-$suffix.csv", Tables.table(dmg_results1; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("permaf-xreg-$suffix.csv", Tables.table(dmg_results2; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("seaice-xreg-$suffix.csv", Tables.table(dmg_results3; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("redisc-xreg-$suffix.csv", Tables.table(dmg_results4; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("to2300-xreg-$suffix.csv", Tables.table(dmg_results5; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("gdppc0-xreg-$suffix.csv", Tables.table(dmg_results6; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("newwts-xreg-$suffix.csv", Tables.table(dmg_results7; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
    CSV.write("burked-xreg-$suffix.csv", Tables.table(dmg_results8; header=[:EU, :US, :OT, :EE, :CA, :IA, :AF, :LA]))
end


text:
export _run_qflex, _run_multi_qflex

"""
run model for Q flexibility objective with AC Power Flow equations
"""

function _run_qflex(file, model_type::_PM.Type, optimizer; kwargs...)
    return _PM.solve_model(file, model_type, optimizer, _build_qflex; kwargs...)
end

"""
given a JuMP model and a PowerModels network data structure,
builds an "reactive power flexibility" formulation of the given data and returns the JuMP model
"""

function _build_qflex(pm::_PM.AbstractPowerModel)

    _PM.variable_bus_voltage(pm)
    _PM.variable_gen_power(pm)
    _PM.variable_branch_power(pm)
    _PM.variable_dcline_power(pm, bounded = false) # TODO: why false?

    objective_qflex(pm)

    _PM.constraint_model_voltage(pm)

    for i in _PM.ids(pm, :ref_buses)
        _PM.constraint_theta_ref(pm, i)
    end

    for i in _PM.ids(pm, :bus)
        _PM.constraint_power_balance(pm, i)
    end

    for (i, branch) in _PM.ref(pm, :branch)
        _PM.constraint_ohms_yt_from(pm, i)
        _PM.constraint_ohms_yt_to(pm, i)
        _PM.constraint_thermal_limit_from(pm, i)
        _PM.constraint_thermal_limit_to(pm, i)
    end

    for i in _PM.ids(pm, :dcline)
        _PM.constraint_dcline_power_losses(pm, i)
    end
end

function objective_qflex(pm::_PM.AbstractPowerModel)

    return JuMP.@objective(pm.model, Min,
    sum((var(pm, :q, (content["element_index"], content["f_bus"], content["t_bus"])) - content["value"])^2
        for (i, content) in pm.ext[:setpoint_q]))
end

function _run_multi_qflex(file, model_type::_PM.Type, optimizer; kwargs...)
    return _PM.solve_model(file, model_type, optimizer, _build_multi_qflex; multinetwork=true, kwargs...)
end

"""
given a JuMP model and a PowerModels network data structure,
builds an "reactive power flexibility" formulation of the given data and returns the JuMP model
"""

function _build_multi_qflex(pm::_PM.AbstractPowerModel)
    for (n, network) in _PM.nws(pm)
            _PM.variable_bus_voltage(pm, nw=n)
            _PM.variable_gen_power(pm, nw=n)
            _PM.variable_branch_power(pm, nw=n)
            _PM.variable_dcline_power(pm, nw=n)

            _PM.constraint_model_voltage(pm, nw=n)

            for i in ids(pm, :ref_buses, nw=n)
                _PM.constraint_theta_ref(pm, i, nw=n)
            end

            for i in ids(pm, :bus, nw=n)
                _PM.constraint_power_balance(pm, i, nw=n)
            end

            for i in ids(pm, :branch, nw=n)
                _PM.constraint_ohms_yt_from(pm, i, nw=n)
                _PM.constraint_ohms_yt_to(pm, i, nw=n)

                _PM.constraint_voltage_angle_difference(pm, i, nw=n)

                _PM.constraint_thermal_limit_from(pm, i, nw=n)
                _PM.constraint_thermal_limit_to(pm, i, nw=n)
            end

            for i in ids(pm, :dcline, nw=n)
                _PM.constraint_dcline_power_losses(pm, i, nw=n)
            end
        end
        objective_multi_qflex(pm)
end

function objective_multi_qflex(pm::_PM.AbstractPowerModel)
    timestep_ids = [id for id in _PM.nw_ids(pm) if id != 0]
    return JuMP.@objective(pm.model, Min,
        sum(
        sum((var(pm, nw, :q, (content["element_index"], content["f_bus"], content["t_bus"])) - content["value"])^2 for (i, content) in pm.ext[:setpoint_q])
        for nw in timestep_ids)
            )
end


text:
strings = ["this", "is", "an", "array"]
ints = [1, 2, 3, 4]
floats = [1, 2, 3.0, 4]
empty = []
println(length(strings))
println(length(empty))

# Arrays can contain mixed elements
mixed = ["string", 1, 3.0]
println(typeof(mixed)) # Array containing elements of Any type

# Nested Arrays
nest = [[1,2], [3,4,5], [6]]
nest2 = [[1,2], "string", [6]]
nest3 = [[1, [3,4]], "string", [6]]
println(typeof(nest), '\t', length(nest))
println(typeof(nest2), '\t', length(nest2))
println(typeof(nest3), '\t', length(nest3))

# Access and assignment
println("ints[3]: $(ints[3])")
println("nest3[1]: $(nest3[1])")
println("nest3[1][2]: $(nest3[1][2])")
ints[3] = 7
println("ints[3]: $(ints[3])")
println("ints: $(ints)")

println("3 ∈ ints : $(3 ∈ ints)")
println("\"an\" ∈ strings : $("an" ∈ strings)")

# Looping
for num in ints
    print(num, '\t')
end
println()

for i in eachindex(strings)
    strings[i] *= 'z'
end
println(strings)

# Slicing
orig = [1, 2, 3, 4, 5, 6]
copy = orig[:]
println("orig == copy: $(orig == copy)")
println("orig === copy: $(orig === copy)")
println("orig[3:end]: $(orig[3:end])")
println("orig[2:4]: $(orig[2:4])")
orig[3:5] = [10, 11, 12]
println("orig: $(orig)")
println("copy: $(copy)")


text:
@testset "153.find-minimum-in-rotated-sorted-array.jl" begin
    @test find_min_154([3 3 3 4 5 6 6 6 1 1 2 2 2 2 3 3][:]) == 1
    @test find_min_154([1,3,5]) == 1
    @test find_min_154([2,2,2,0,1]) == 0
    @test find_min_154([3,3,1,3]) == 1
end


text:
mutable struct position_3d_t <: LCMType
    translation::vector_3d_t
    rotation::quaternion_t
end

@lcmtypesetup(position_3d_t)


text:
using Compat
using Patchwork
import Patchwork: jsonfmt, refdiff

# Compare diffs produced by
# Patchwork and virtual-dom
# This prints stuff in the JS console as well

macro compare(a, b)
    hint = string("<> ", a, " -> ", b)
    quote
        refdiff($(esc(a)), $(esc(b)), label=$hint)
        println($hint)
        println(jsonfmt(diff($(esc(a)), $(esc(b)))))
    end
end

a = em("x")
b = a & @compat Dict(:className => "cls")

p,q,r = li("p"), li("q"), li("r")

@compare a a
@compare span(a) span(a, "a")
@compare span(a, "a") span(a)
@compare span(a, "a") span(a, "b")
@compare span(a, "a", b) span(a, "a")
@compare a b
@compare a & @compat Dict(:style => [:color => "red"]) a & Dict(:style => [:color => "blue"])
@compare b a
@compare em(a) em(a)
@compare em(a) em(b)
@compare em(b) em(a)
@compare span([a, b]) span([a, b])
@compare span([a, b]) span([b, a])
@compare ul(p, q, r) ul(p, q, r)
@compare ul(p, q, r) ul(p, r)
@compare ul(p, ul(q, q), r) ul(p, ul(q, q), r)
@compare ul(p, ul(q, q), r) ul(p, ul(q, q & attrs(className="clas")), r & attrs(id="id"))
@compare ul(p, ul(q, q, key=symbol("12")), r) ul(p, ul(q, q, key=symbol("12"))  & attrs(className="clas")+ r & attrs(id="id"))
@compare ul(p, q, r) ul(p, r)
@compare ul(p, r) ul(p, q, r)
@compare ul(p, r) ul(p, q)
@compare ul(p, r) ul(q, p)
@compare ul(p, r) ul(q, p & attrs(className="cls"))


text:
using Test

using SwapLiterals
using SwapLiterals: floats_use_rationalize!

using BitIntegers, SaferIntegers

literalswapper(sw::Pair...) = SwapLiterals.literals_swapper(sw)

literalswapper(F, I, I128, B=nothing) =
    literalswapper(Float64=>F, Int=>I, Int128=>I128, BigInt=>B)

makeset(ex) = Expr(:call, :Set, Expr(:vect, ex.args...))

# just a random function, which transform `a := x` into `A := x`
makecoloneq(ex) = Expr(:(=),
                       Symbol(uppercase(String(ex.args[1]))),
                       ex.args[2:end]...)

@testset "swapliterals" begin
    swapbig = literalswapper(:BigFloat, :big, "@big_str")
    @test swapbig(1) == :(big(1))
    @test swapbig(1.2) == :(BigFloat(1.2))

    @swapliterals :BigFloat :big "@big_str" begin
        @test 1 == Int(1)
        @test 1 isa BigInt
        @test 10000000000 isa BigInt # Int64 literals are transformed like Int literals
        @test 1.2 isa BigFloat
        @test 1.2 == big"1.1999999999999999555910790149937383830547332763671875"
        @test 1.0 == Float64(1.0)
        @test $1 isa Int
        @test $10000000000 isa Int64 # even on 32-bits systems
        @test $1.2 isa Float64
    end

    # next three blocs should be equivalent
    @swapliterals begin
        @test 1.2 isa BigFloat
        @test 1.2 == big"1.2"
        @test 1 isa BigInt
        @test 10000000000 isa BigInt # Int64 literals are transformed like Int literals
        @test 11111111111111111111 isa BigInt
    end
    @swapliterals "@big_str" :big :big begin
        @test 1.2 isa BigFloat
        @test 1.2 == big"1.2"
        @test 1 isa BigInt
        @test 11111111111111111111 isa BigInt
    end
    @swapliterals Float64=>"@big_str" Int=>:big Int128=>:big begin
        @test 1.2 isa BigFloat
        @test 1.2 == big"1.2"
        @test 1 isa BigInt
        if Int === Int64
            @test 10000000000 isa BigInt
        else
            @test 10000000000 isa Int64 # Int64 literals are *not* transformed like Int literals
        end
        @test 11111111111111111111 isa BigInt
    end

    @swapliterals Int64 => :big begin
        if Int === Int64
            @test 1 isa BigInt
        else
            @test 1 isa Int
        end
        @test 10000000000 isa BigInt
    end

    # TODO: these tests in loop are dubious
    for T in Base.BitUnsigned_types
        @test typeof(swapbig(T(1))) == T
    end
    for T in [Float32, Float16]
        @test typeof(swapbig(T(1))) == T
    end

    x = eval(swapbig(1.0))
    @test x isa BigFloat && x == 1.0
    x = eval(swapbig(1))
    @test x == 1 && x isa BigInt
    x = eval(swapbig(:11111111111111111111))
    @test x == 11111111111111111111 && x isa BigInt
    x = eval(swapbig(:1111111111111111111111111111111111111111))
    @test x isa BigInt

    @swapliterals :BigFloat :big "@big_str" begin
        x = 1.0
        @test x isa BigFloat && x == Float64(1.0)
        x = 1
        @test x == Int(1) && x isa BigInt
        x = 11111111111111111111
        @test x == big"11111111111111111111" && x isa BigInt
        x = 1111111111111111111111111111111111111111
        @test x isa BigInt

        @test $1.0 isa Float64
        @test $11111111111111111111 isa Int128
        # throws as un-handled quote
        # @test $1111111111111111111111111111111111111111 isa BigInt
    end

    swap128 = literalswapper(:Float64, :Int128, "@int128_str")
    x = eval(swap128(1))
    @test x == 1 && x isa Int128
    x = eval(swap128(:11111111111111111111))
    @test x == 11111111111111111111 && x isa Int128
    x = eval(swap128(:1111111111111111111111111111111111111111))
    @test x isa BigInt

    @swapliterals :Float64 :Int128 "@int128_str" begin
        x = 1
        @test x == Int(1) && x isa Int128
        x = 10000000000
        @test x == big"10000000000"
        @test x isa Int128 # Int64 literals are transformed like Int literals

        x = 11111111111111111111
        @test x == Int128(11111111111111111111) && x isa Int128
        x = 1111111111111111111111111111111111111111
        @test x isa BigInt
    end

    swapnothing = literalswapper(nothing, nothing, nothing)
    x = eval(swapnothing(1.0))
    @test x isa Float64
    x = eval(swapnothing(:11111111111111111111))
    @test x isa Int128
    x = eval(swapnothing(:1111111111111111111111111111111111111111))
    @test x isa BigInt

    @swapliterals nothing nothing nothing begin
        x = 1.0
        @test x isa Float64
        @test 1 isa Int
        @test 10000000000 isa Int64
        x = 11111111111111111111
        @test x isa Int128
        x = 1111111111111111111111111111111111111111
        @test x isa BigInt
    end

    # pass :big instead of a string macro
    swaponly128 = literalswapper(nothing, nothing, :big)
    x = eval(swaponly128(:11111111111111111111))
    @test x isa BigInt

    @swapliterals nothing nothing :big begin
        x = 11111111111111111111
        @test x isa BigInt
    end

    # pass symbol for Int128
    swapBitIntegers = literalswapper(nothing, :Int256, :Int256)
    x = eval(swapBitIntegers(123))
    @test x isa Int256
    x = eval(swapBitIntegers(:11111111111111111111))
    @test x isa Int256

    swapSaferIntegers = literalswapper(nothing, :SafeInt, :SafeInt128)
    x = eval(swapSaferIntegers(123))
    @test x isa SafeInt
    x = eval(swapSaferIntegers(:11111111111111111111))
    @test x isa SafeInt128

    @swapliterals nothing :Int256 :Int256 begin
        x = 123
        @test x isa Int256
        x = 11111111111111111111
        @test x isa Int256
    end

    @swapliterals nothing :SafeInt :SafeInt128 begin
        x = 123
        @test x isa SafeInt
        x = 11111111111111111111
        @test x isa SafeInt128
    end

    # pass symbol for BigInt
    swapbig = literalswapper(nothing, nothing, :Int1024, :Int1024)
    x = eval(swapbig(:11111111111111111111))
    @test x isa Int1024
    x = eval(swapbig(:1111111111111111111111111111111111111111))
    @test x isa Int1024

    @swapliterals nothing nothing :Int1024 :Int1024 begin
        @test 11111111111111111111 isa Int1024
        @test 1111111111111111111111111111111111111111 isa Int1024
        @test $11111111111111111111 isa Int128
        @test $1111111111111111111111111111111111111111 isa BigInt
    end

    swapbig = literalswapper(nothing, nothing, :big, :big)
    x = eval(swapbig(:11111111111111111111))
    @test x isa BigInt
    x = eval(swapbig(:1111111111111111111111111111111111111111))
    @test x isa BigInt

    @swapliterals nothing nothing :big :big begin
        x = 11111111111111111111
        @test x isa BigInt
        x = 1111111111111111111111111111111111111111
        @test x isa BigInt
    end

    # kwargs
    kwswapper = literalswapper(Int=>:big)
    @test eval(kwswapper(1.2)) isa Float64
    @test eval(kwswapper(1)) isa BigInt
    @test eval(kwswapper(:11111111111111111111)) isa Int128
    @test eval(kwswapper(:1111111111111111111111111111111111111111)) isa BigInt

    # Float32
    @swapliterals Float32 => :big begin
        @test 1.2f0 == big"1.2000000476837158203125"
    end
    @swapliterals Float32 => "@big_str" begin
        @test 1.2f0 == big"1.2"
    end

    # string swappers
    @swapliterals Int => "@big_str" UInt8 => "@raw_str" begin
        @test 1 isa BigInt
        @test 0x01 === "1"
    end

    # Int & UInt
    @swapliterals Int => :UInt8 UInt => :Int8 begin
        @test 1 isa UInt8
        if UInt === UInt64
            @test 0x0000000000000001 isa Int8
            @test 0x00000001 isa UInt32
        else
            @test 0x00000001 isa Int8
            @test 0x0000000000000001 isa UInt64
        end
    end
    @swapliterals Int32=>:UInt8 Int64=>:UInt8 UInt64=>:Int8 begin
        @test 1 isa UInt8
        @test 0x0000000000000001 isa Int8
    end

    # unsigned
    @swapliterals UInt8=>:Int UInt16=>:Int UInt32=>:Int UInt64=>:Int UInt128=>:Int128 begin
        @test 0x1 isa Int
        @test 0x0001 isa Int
        @test $0x0001 isa UInt16
        @test 0x00000001 isa Int
        @test $0x00000001 isa UInt32
        @test :0x00000001 isa UInt32
        @test 0x0000000000000001 isa Int
        @test :0x0000000000000001 isa UInt64
        @test 0x00000000000000000000000000000001 isa Int128
        @test $0x00000000000000000000000000000001 isa UInt128
    end
    @swapliterals UInt128=>"@int128_str" begin
        @test 0x00000000000000000000000000000001 isa Int128
        @test $0x00000000000000000000000000000001 isa UInt128
    end

    # strings & chars
    @swapliterals Char=>:string String => "@r_str" begin
        @test "123" isa Regex
        @test 'a' isa String
    end

    @swapliterals Char => :UInt64 String => :Symbol begin
        @test "123" isa Symbol
        @test 'a' === 0x0000000000000061
    end

    @test_throws ArgumentError literalswapper(Array=>:Int)

    # quoted expressions which are not handled by a swapper remain quoted
    swapper = literalswapper(Char=>UInt)
    @test 'a' !== swapper('a') == 0x61
    @test swapper(Expr(:$, 'a')) === 'a'
    @test_throws ErrorException eval(swapper(Expr(:$, 1)))
    @test_throws ErrorException eval(swapper(Expr(:$, :11111111111111111111)))
    @test_throws ErrorException eval(swapper(Expr(:$, :1111111111111111111111111111111111111111)))

    # function swappers
    @swapliterals UInt8 => (x -> x+1) Int => UInt8 Int128 => (ex -> ex.args[3]) begin
        @test 0x01 == 2.0
        @test 1 isa UInt8
        @test 11111111111111111111 == "11111111111111111111"
    end

    # :braces, :tuple, :vect
    @swapliterals [:braces => makeset,
                   :bracescat => makeset,
                   :tuple => makeset,
                   :vect => makeset
                   ] begin
        r = push!(Set{Int}(), 1, 2, 3)
        s = {1, 2, 3}
        @test s isa Set{Int}
        @test s == r
        s = [1, 2, 3]
        @test s isa Set{Int}
        @test s == r
        s = (1, 2, 3)
        @test s isa Set{Int}
        @test s == r
        s = { 1
              2
              3 }
        @test s isa Set{Int}
        @test s == r
    end
    @swapliterals :tuple => :collect begin
        v = (1, 2, 3)
        @test v isa Vector{Int}
        @test v == [1, 2, 3]
    end

    # :(:=)
    @swapliterals :(:=) => makecoloneq begin
        a := 2
        @test A == 2
        @test !@isdefined(a)
        bcd := 1, 2, 3
        @test BCD == (1, 2, 3)
        @test !@isdefined(bcd)
    end

    # test that swapper is applied recursively
    @swapliterals :tuple => makeset Int => :big begin
        @test (1, 2) isa Set{BigInt}
    end
    @swapliterals Int => :big begin
        @test (1, 2) isa Tuple{BigInt,BigInt}
    end
end

# test name resolution for functions
module TestModule

using SwapLiterals, Test

uint8(x) = UInt8(x)

@swapliterals Int => uint8  Char => (x -> uint8(x)+1) begin
    @test 1 isa UInt8
    @test 'a' == 0x62
end
end

## playing with floats_use_rationalize!()

# can't be in a @testset apparently, probably because the parsing
# in @testset is done before floats_use_rationalize!() takes effect

@swapliterals Float32="@big_str" Float64="@big_str" begin
    @test 1.2 == big"1.2"
    @test 1.2f0 == big"1.2"
end

floats_use_rationalize!()
@swapliterals Float32="@big_str" Float64="@big_str" begin
    @test 1.2 == big"1.2"
    @test 1.2f0 == big"1.2"
end

# try again, with explicit `true` arg, and with :BigFloat instead of :big
floats_use_rationalize!(true)
@swapliterals Float32=:BigFloat Float64=:BigFloat begin
    @test 1.2 == big"1.2"
    @test 1.2f0 == big"1.2"
end

floats_use_rationalize!(false)
@swapliterals Float32=:BigFloat Float64=:BigFloat begin
    @test 1.2 == big"1.1999999999999999555910790149937383830547332763671875"
    @test 1.2f0 == big"1.2000000476837158203125"
end


text:
"""
    cleave(seq::AminoAcidSequence, peptidase::Function = trypsin)

Predict cleavage sites. Returns a vector of cleavage site indices (excluding gaps) and a vector containing the resulting peptides. `protease` is a function that takes two arguments (`seq` and an index `i`) and returns `true` if the sequence would be cleaved between position `i` and `i+1` by the given peptidase.
Available (but unexported) options for `protease` are
- `trypsin`
- `proteinaseK`
- `lysC`
- `proline_endopeptidase`
For details on cleavage sites, see https://web.expasy.org/peptide_cutter/peptidecutter_enzymes.html
"""
function cleave(seq::LongAminoAcidSeq, protease::Function = trypsin)
    seq = ungap(seq)
    cleavagesites = Int[]
    last = 0
    for i in 1:length(seq)
        if protease(seq, i)
            push!(cleavagesites, i)
            last = i
        end
    end
    peptides = getpeptides(seq, cleavagesites)
    return (positions = cleavagesites, peptides = peptides)
end

function getpeptides(seq, sites)
    peptides = Array{LongAminoAcidSeq}(undef, length(sites)+1)
    last = 1
    for (i, site) in enumerate(sites)
        peptides[i] = seq[last:site]
        last = site + 1
    end
    peptides[end] = seq[last:end]
    return peptides
end

function trypsin(seq, i)
    if i-1 <= 0 || i+1 > lastindex(seq)
        return false
    end
    window = seq[i-1:i+1]
    if window == aa"MRP" || window == aa"WKP"
        return true
    elseif window == aa"CKD" || window == aa"DKD" ||
           window == aa"CKH" || window == aa"CKY" ||
           window == aa"CRK" || window == aa"RRR" ||
           window == aa"RRH"
        return false
    elseif (window[2] == AA_K || window[2] == AA_R) && window[3] != AA_P
        return true
    end
    return false
end

function proteinaseK(seq, i)
    i >= lastindex(seq) && return false
    if seq[i] in [AA_A, AA_E, AA_F, AA_I, AA_L, AA_J, AA_T, AA_V, AA_W, AA_Y]
        return true
    end
    return false
end

function lysC(seq, i)
    i >= lastindex(seq) && return false
    seq[i] == AA_K && return true
    return false
end

function proline_endopeptidase(seq, i)
    i >= lastindex(seq) && return false
    if seq[i] == AA_P && seq[i-1] in [AA_H, AA_K, AA_R] && seq[i+1] != AA_P
        return true
    end
    return false
end


text:
export SimulationBuilder, set!, build

mutable struct SimulationBuilder{N}
    unit_system::UnitSystem
    kernel::Kernel{N}
    eos::EquationOfState
    mass_distribution::InitialMassDistribution{N}
    velocity_distribution::InitialVelocityDistribution{N}
    energy_distribution::InitialEnergyDistribution{N}
    mass_component_type::Tuple{Type,Any}
    energy_component_type::Tuple{Type,Any}
    diffusion_component_type::Tuple{Type,Any}
    boundary_component_type::Tuple{Type,Any}
    stepper_type::Tuple{Type,Any}

    function SimulationBuilder(
        mass_distribution::InitialMassDistribution{N},
        energy_distribution::InitialEnergyDistribution{N};
        unit_system::UnitSystem = SI,
        kernel::Kernel{N} = TabulatedKernel(CubicSplineKernel{N}()),
        eos::EquationOfState = AdiabaticIdealGasEOS(1.4, NORMALIZED),
        velocity_distribution::InitialVelocityDistribution{N} = StaticVelocityDistribution{
            N,
        }(),
    ) where {N}
        new{N}(
            unit_system,
            kernel,
            eos,
            mass_distribution,
            velocity_distribution,
            energy_distribution,
            (AdaptiveMass, []),
            (StandardEnergy, []),
            (StandardDiffusion, []),
            (NoBoundaries{N}, []),
            (EulerStepper, []),
        )
    end
end

function Base.show(io::IO, builder::SimulationBuilder{N}) where {N}
    print(
        io,
        """
        SimulationBuilder{
            $(N)D
            $(typeof(builder.mass_distribution))
            $(typeof(builder.velocity_distribution))
            $(typeof(builder.energy_distribution))
            $(typeof(builder.eos))
            $(typeof(builder.kernel))
            $(builder.mass_component_type[1])
            $(builder.energy_component_type[1])
            $(builder.diffusion_component_type[1])
            $(builder.boundary_component_type[1])
            $(builder.stepper_type[1])
        }""",
    )
end

(Base.:+)(builder::SimulationBuilder, property) = set!(builder, property)
(Base.:+)(
    builder::SimulationBuilder,
    (property, kwargs)::Tuple{<:Any,NamedTuple},
) = set!(builder, property, kwargs...)

function set!(builder::SimulationBuilder{N}, unit_system::UnitSystem) where {N}
    builder.unit_system = unit_system
    builder
end

function set!(builder::SimulationBuilder{N}, kernel::Kernel{N}) where {N}
    builder.kernel = kernel
    builder
end

function set!(builder::SimulationBuilder{N}, eos::EquationOfState) where {N}
    builder.eos = eos
    builder
end

function set!(
    builder::SimulationBuilder{N},
    mass_distribution::InitialMassDistribution{N},
) where {N}
    builder.mass_distribution = mass_distribution
    builder
end

function set!(
    builder::SimulationBuilder{N},
    velocity_distribution::InitialVelocityDistribution{N},
) where {N}
    builder.velocity_distribution = velocity_distribution
    builder
end

function set!(
    builder::SimulationBuilder{N},
    energy_distribution::InitialEnergyDistribution{N},
) where {N}
    builder.energy_distribution = energy_distribution
    builder
end

function set!(
    builder::SimulationBuilder{N},
    mass_component_type::Type{<:MassComponent};
    kwargs...,
) where {N}
    builder.mass_component_type = (mass_component_type, kwargs)
    builder
end

function set!(
    builder::SimulationBuilder{N},
    energy_component_type::Type{<:EnergyComponent};
    kwargs...,
) where {N}
    builder.energy_component_type = (energy_component_type, kwargs)
    builder
end

function set!(
    builder::SimulationBuilder{N},
    diffusion_component_type::Type{<:DiffusionComponent};
    kwargs...,
) where {N}
    builder.diffusion_component_type = (diffusion_component_type, kwargs)
    builder
end

function set!(
    builder::SimulationBuilder{N},
    boundary_component_type::Type{<:BoundaryComponent{N}};
    kwargs...,
) where {N}
    builder.boundary_component_type = (boundary_component_type, kwargs)
    builder
end

function set!(
    builder::SimulationBuilder{N},
    stepper_type::Type{<:Stepper{N}};
    kwargs...,
) where {N}
    builder.stepper_type = (stepper_type, kwargs)
    builder
end

function build(builder::SimulationBuilder{N}) where {N}
    mass_component, positions = builder.mass_component_type[1](
        builder.kernel,
        builder.mass_distribution;
        builder.mass_component_type[2]...,
    )
    mass_densities = get_mass_densities(mass_component)
    energy_component = builder.energy_component_type[1](
        positions,
        mass_densities,
        builder.energy_distribution,
        builder.eos;
        builder.energy_component_type[2]...,
    )
    diffusion_component =
        builder.diffusion_component_type[1](builder.diffusion_component_type[2]...)
    velocities =
        Velocities(positions, mass_densities, builder.velocity_distribution)
    boundary_component =
        builder.boundary_component_type[1](builder.boundary_component_type[2]...)
    stepper = builder.stepper_type[1](
        positions,
        velocities,
        mass_component,
        energy_component,
        diffusion_component,
        boundary_component,
        builder.eos;
        builder.stepper_type[2]...,
    )
    stepper
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
# utility functions needed for sandshark

datadir = joinpath(ENV["HOME"],"data","sandshark","full_wombat_2018_07_09","extracted")
matcheddir = joinpath(datadir, "matchedfilter", "particles")
beamdir = joinpath(datadir, "beamformer", "particles")

function loaddircsvs(datadir)
  # https://docs.julialang.org/en/v0.6.1/stdlib/file/#Base.Filesystem.walkdir
  datadict = Dict{Int, Array{Float64}}()
  for (root, dirs, files) in walkdir(datadir)
    # println("Files in $root")
    for file in files
      # println(joinpath(root, file)) # path to files
      data = readdlm(joinpath(root, file),',')
      datadict[parse(Int,split(file,'.')[1])/1] = data
    end
  end
  return datadict
end

rangedata = loaddircsvs(matcheddir)
azidata = loaddircsvs(beamdir)
timestamps = intersect(sort(collect(keys(rangedata))), sort(collect(keys(azidata))))

# NAV data
navdata = Dict{Int, Vector{Float64}}()
navfile = readdlm(joinpath(datadir, "nav_data.csv"))
for row in navfile
    s = split(row, ",")
    id = round(Int, 1e9*parse(Float64, s[1]))
    # round(Int, 1000 * parse(s[1])) = 1531153292381
    navdata[id] = parse.(Float64,s)
end
navkeys = sort(collect(keys(navdata)))
# NAV colums are X,Y = 7,8
# lat,long = 9,10
# time,pitch,roll,heading,speed,[Something], internal_x,internal_y,internal_lat,internal_long, yaw_rad

# LBL data - note the timestamps need to be exported as float in future.
lbldata = Dict{Int,  Vector{Float64}}()
lblfile = readdlm(joinpath(datadir, "lbl.csv"))
for row in lblfile
    s = split(row, ",")
    id = round(Int, 1e9*parse(Float64, s[1]))
    if s[2] != "NaN"
        lbldata[id] = parse.(Float64, s)
    end
end
lblkeys = sort(collect(keys(lbldata)))

# function heading2yaw(heading)
# heading = map( x->deg2rad(navdata[x][4]), navkeys )
# wrapheading = TU.wrapRad.(heading)
# end

# GET Y = north,  X = East,  Heading along +Y clockwise [0,360)]
# east = Float64[]
# north = Float64[]
# heading = Float64[]

# WANT X = North,  Y = West,  Yaw is right and rule from +X (0) towards +Y pi/2, [-pi,pi)
# so the drawPoses picture will look flipped from Nicks picture
# remember theta = atan2(y,x)    # this is right hand rule
X = Float64[]
Y = Float64[]
yaw = Float64[]
for id in navkeys
  # push!(east, getindex(navdata[id],7)) # x-column csv
  # push!(north, getindex(navdata[id],8)) # y-column csv
  # push!(heading, getindex(navdata[id],4))
  # push!(yaw, TU.wrapRad(-deg2rad(getindex(navdata[id],4))))

  push!(X, 0.7*getindex(navdata[id],7) ) # 8
  push!(Y, 0.7*getindex(navdata[id],8) ) # 7
  push!(yaw, TU.wrapRad(pi/2-deg2rad(getindex(navdata[id],4))) )  # rotation about +Z
end

lblX = Float64[]
lblY = Float64[]
for id in lblkeys
    push!(lblX, getindex(lbldata[id],2) )
    push!(lblY, getindex(lbldata[id],3) )
end

# Build interpolators for x, y, yaw
interp_x = LinearInterpolation(navkeys, X)
interp_y = LinearInterpolation(navkeys, Y)
interp_yaw = LinearInterpolation(navkeys, yaw)

## Caching factors
ppbrDict = Dict{Int, Pose2Point2BearingRange}()
odoDict = Dict{Int, Pose2Pose2}()
NAV = Dict{Int, Vector{Float64}}()

## what what

# Gadfly.plot(y=180/pi.*yaw, Geom.path)

#


text:
module ExponentialIntegrators

using ExponentialUtilities

import ExponentialUtilities: expv!

abstract type AbstractIntegrator{T} end

order(::AbstractIntegrator) = 1
order(::Type{AbstractIntegrator}) = 1

abstract type AbstractExponentiator{T} end

exact(::AbstractExponentiator) = false
inplace(::AbstractExponentiator) = false
order(::AbstractExponentiator) = 1

# The last argument to `expv!`, the integer, signifies if the
# exponentiation should be performed in the forward direction (1), or
# backward direction (0). This is useful if `E` split further, and is
# part of a symmetric splitting scheme. If the argument is 0, the
# order does not matter (as is the case for the middle term of a
# certain splitting).
expv!(v, E::AbstractExponentiator, Δt, u) = expv!(v, E, Δt, u, 0)

abstract type AbstractOperator{T} end
# Time-dependent operators could be of the kind A(t), A + f(t)B, etc,
# for which the different temporal integration routines would have to
# be special-cased.

checkorder(o,integrator) =
    o < order(integrator) && @warn("Order of $(typeof(o)) lower than required for $(integrator)")

macro swap!(x,y)
    quote
        local tmp = $(esc(x))
        $(esc(x)) = $(esc(y))
        $(esc(y)) = tmp
    end
end

# Temporal integration of operators
include("integrate.jl")

# Exponentiators for single operators
include("krylov.jl")
include("pade.jl")

# Integrators composed of multiple exponentiators
include("strang_splitting.jl")
include("exact_integrator.jl")

export step!, StrangSplitting, ExactIntegrator

end # module


text:
include("test_operators.jl")
include("test_systems.jl")
include("test_fourier.jl")
include("test_basis_utils.jl")
include("test_time_evolution.jl")
include("test_composite_systems.jl")
include("integration_tests.jl")


text:
@testset "Spine2D" begin
    @polyvar x y
    f = x^2*y + y^2 + 3x^2*y^3 + y^4 + x^4*y^4

    S = spine(f)
    @test length(components_complement(S)) == 5
    @test sprint(show, S) ==
        """Spine2D{DynamicPolynomials.Polynomial{true,Tropical{Float64}}} with:
            $(" ")* 4 vertices
            $(" ")* 4 half-rays
            $(" ")* 1 bounded cells"""
    S = spine(f, minimal_component_size=0.05)
    @test length(components_complement(S)) == 5

    S = spine(f, domain=(-5, 5, -5, 5), minimal_component_size=0.01)
    @test length(components_complement(S)) == 5
end


text:
# =====================================
# Crane
#
# system type: linear continuous system
# state dimension: 6
# =====================================
using MathematicalSystems

function crane_model()
    # system matrix
    A = [ 0.0           1.0               0.0   0.0   0.0           0.0;
         -0.417533     -3.1931759963     39.24  0.0 -14.825331     11.123344;
          0.0           0.0               0.0   1.0   0.0           0.0;
          0.0417533     0.31931759963    -4.905 0.0   1.4825331    -1.1123344;
          0.0638407957 -0.32473339016573  0.0   0.0  -3.7332068901 -0.7007592976;
          0.0853437452 -0.72366802635628  0.0   0.0  -5.9714023436 -2.2736115136]

    # linear continuous system
    S = @system(x' = Ax)

    return S
end

max_stars_repo_name: JuliaReach/ReachabilityBenchmarks

max_issues_repo_name: JuliaReach/ReachabilityBenchmarks

max_forks_repo_name: JuliaReach/ReachabilityBenchmarks


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "Spot"
version = v"2.8.1"

# Collection of sources required to build SpotBuilder
sources = [
    "http://www.lrde.epita.fr/dload/spot/spot-2.8.1.tar.gz" =>
    "dcb7aa684725304afb3d435f26f25b51fbd6e9a6ef610e16163cc0030ad5eab4"
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd spot-2.8.1

if [[ "${target}" == *-freebsd* ]] || [[ "${target}" == *-apple-* ]]; then
    export FC=/opt/${target}/bin/${target}-gfortran
    export LD=/opt/${target}/bin/${target}-ld
    export AR=/opt/${target}/bin/${target}-ar
    export AS=/opt/${target}/bin/${target}-as
    export NM=/opt/${target}/bin/${target}-nm
    export OBJDUMP=/opt/${target}/bin/${target}-objdump
fi

./configure --prefix=$prefix --host=${target} --disable-python
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc7)),
    Windows(:x86_64, compiler_abi=CompilerABI(:gcc8)),
    Linux(:x86_64, compiler_abi=CompilerABI(:gcc7)),
    Linux(:x86_64, compiler_abi=CompilerABI(:gcc8)),
    MacOS()
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libspot", :libspot)
]

# Dependencies that must be installed before this package can be built
dependencies = [

]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)


text:
abstract type SauterSchwab3DStrategy end

#= Tetrahedron-Tetrahedron Interaction =#
struct CommonVolume6D{S,Q}      <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonFace6D{S,Q}        <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge6D{S,Q}        <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex6D{S,Q}      <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance6D{S,Q}  <: SauterSchwab3DStrategy sing::S; qps::Q end

struct CommonVolume6D_S{S,Q}     <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonFace6D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge6D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex6D_S{S,Q}     <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance6D_S{S,Q} <: SauterSchwab3DStrategy sing::S; qps::Q end

#= Tetrahedron-Triangle Interaction =#
struct CommonFace5D{S,Q}         <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge5D{S,Q}         <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex5D{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance5D{S,Q}   <: SauterSchwab3DStrategy sing::S; qps::Q end

struct CommonFace5D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge5D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex5D_S{S,Q}     <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance5D_S{S,Q} <: SauterSchwab3DStrategy sing::S; qps::Q end

#= Triangle-Triangle Interaction =#
struct CommonFace4D{S,Q}         <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge4D{S,Q}         <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex4D{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance4D{S,Q}   <: SauterSchwab3DStrategy sing::S; qps::Q end

struct CommonFace4D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonEdge4D_S{S,Q}       <: SauterSchwab3DStrategy sing::S; qps::Q end
struct CommonVertex4D_S{S,Q}     <: SauterSchwab3DStrategy sing::S; qps::Q end
struct PositiveDistance4D_S{S,Q} <: SauterSchwab3DStrategy sing::S; qps::Q end

function _legendre(n,a,b)
    x, w = FastGaussQuadrature.gausslegendre(n)
    w .*= (b-a)/2
    x = (x.+1)/2*(b-a).+a
    collect(zip(x,w))
end

@inbounds function calc_vol(vertices::SMatrix{D,N,T}) where {D,N,T}
    @assert N == D + 1
    X = SMatrix{D,D,T}(vertices[i, j + 1] - vertices[i, 1]
                       for i in 1:D, j in 1:D)
    vol = det(X)
    return vol
end

function get_pts_wts(scheme::TnScheme{N,T},  vertices::SMatrix{D,N,U}) where {N,T,D,U}
    @assert N > 0
    @assert N >= D + 1

    ws = scheme.weights
    ps = scheme.points
    @assert length(ws) == length(ps)

    xs  = Vector{SVector{D,T}}()
    for i in 1:length(ws)

        p = ps[i]

        push!(xs,vertices * p)

    end

    vol = calc_vol(vertices) / factorial(N - 1)

    return xs, ws*vol
end

function get_pts_wts(scheme::TnScheme{N,T}, vertices::SVector{N,SVector{D,U}}) where {N,T,D,U}
    return get_pts_wts( scheme,SMatrix{D,N,U}(vertices[n][d] for d in 1:D, n in 1:N))
end

function get_pts_wts( scheme::TnScheme{N}, vertices::AbstractVector) where {N}
    @assert length(vertices) == N
    @assert N > 0
    D = length(vertices[1])
    @assert N >= D + 1
    vertices′ = SVector{N}(map(SVector{D}, vertices))
    return get_pts_wts(scheme, vertices′)
end

function _grundmannMoeller2D(n)
    scheme = GrundmannMoeller.grundmann_moeller(Float64,Val(2),n)
    vertices=[[0,0 ], [1,0], [1,1]]
    x,w = get_pts_wts(scheme,vertices)

    zip(x,w)
end

function _grundmannMoeller3D(n)
    scheme = GrundmannMoeller.grundmann_moeller(Float64,Val(3),n)
    vertices=[[0 0 0], [1,0,0], [1,1,0], [1,1,1] ]
    x,w = get_pts_wts(scheme,vertices)

    zip(x,w)
end

function _grundmannMoeller4D(n)
    scheme = GrundmannMoeller.grundmann_moeller(Float64,Val(4),n)
    vertices=[[0,0,0,0], [1,0,0,0], [1,1,0,0], [1,1,1,0], [1,1,1,1]]
    x,w = get_pts_wts(scheme,vertices)

    zip(x,w)
end

function _shunnham2D(m)
    vertices =  [[0 0], [1 0], [1 1]]
    scheme = ShunnHamQuadrature.shunnham2D_ref(m)

    x,w = ShunnHamQuadrature.get_pts_wts(scheme,vertices)

    zip(x,w)
end

function _shunnham3D(m)
    vertices =  [[0 0 0], [1 0 0], [1 1 0], [1 1 1]]
    scheme = ShunnHamQuadrature.shunnham3D_ref(m)

    x,w = ShunnHamQuadrature.get_pts_wts(scheme,vertices)

    zip(x,w)
end

function _shunnham4D(m)
    vertices =  [[0 0 0 0], [1 0 0 0], [1 1 0 0], [1 1 1 0], [1 1 1 1]]
    scheme = ShunnHamQuadrature.shunnham4D(m)

    x,w = ShunnHamQuadrature.get_pts_wts(scheme,vertices)

    zip(x,w)
end

#=Tetrahedron-tetrahedron interaction =#
#Common Volume Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonVolume6D)

    qps = method.qps
	sum(w1*w2*w3*w4*w5*w6*k6p_ct(integrand, η1, η2, η3, η4, ξ1, ξ2)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (η4, w4) in qps, (ξ1, w5) in qps, (ξ2, w6) in qps)
end
#Common Volume Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonVolume6D_S)

    qps = method.qps
	sum(w1*w2*w3*k6p_ct_A(integrand, x[1], x[2], x[3], x[4], y, z)
		for (x, w1) in qps[3], (y, w2) in qps[1], (z, w3) in qps[1]) +
    sum(w1*w2*k6p_ct_B(integrand, x[1], x[2], x[3], x[4], y[1], y[2])
		for (x, w1) in qps[3], (y, w2) in qps[2])
end

#Common Face Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonFace6D)

	qps = method.qps
	sum(w1*w2*w3*w4*w5*w6*k6p_cf(integrand, η1, η2, η3, η4, ξ1, ξ2)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (η4, w4) in qps, (ξ1, w5) in qps, (ξ2, w6) in qps)
end

#Common Face Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonFace6D_S)

	qps = method.qps
	sum(w1*w2*w3*k6p_cf_A(integrand, x[1], x[2], x[3], y[1], y[2], z)
		for (x, w1) in qps[3], (y, w2) in qps[2], (z, w3) in qps[1])+
    sum(w1*w2*w3*k6p_cf_B(integrand, x[1], x[2], x[3], y[1], y[2], z)
		for (x, w1) in qps[3], (y, w2) in qps[2], (z, w3) in qps[1])+
    sum(w1*w2*k6p_cf_C(integrand, x[1], x[2], x[3], y[1], y[2], y[3])
		for (x, w1) in qps[3], (y, w2) in qps[3])
end

#Common Edge Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonEdge6D)
	qps = method.qps
	sum(w1*w2*w3*w4*w5*w6*k6p_ce(integrand, η1, η2, η3, η4, ξ1, ξ2)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (η4, w4) in qps, (ξ1, w5) in qps, (ξ2, w6) in qps)
end

#Common Simplex Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonEdge6D_S)
	qps = method.qps
	sum(w1*w2*w3*k6p_ce_A(integrand, x[1], x[2], y[1], y[2], z[1], z[2])
		for (x, w1) in qps[2], (y, w2) in qps[2], (z, w3) in qps[2]) +
    sum(w1*w2*w3*k6p_ce_B(integrand, x[1], x[2], y[1], y[2], y[3], z)
		for (x, w1) in qps[2], (y, w2) in qps[3], (z, w3) in qps[1]) +
    sum(w1*w2*k6p_ce_C(integrand, x[1], x[2], y[1], y[2], y[3], y[4])
		for (x, w1) in qps[2], (y, w2) in qps[4] ) +
    sum(w1*w2*w3*k6p_ce_D(integrand, x[1], x[2], y, z[1], z[2], z[3])
		for (x, w1) in qps[2], (y, w2) in qps[1], (z, w3) in qps[3]) +
    sum(w1*w2*w3*w4*k6p_ce_E(integrand, x[1], x[2], y, z[1], z[2], w)
		for (x, w1) in qps[2], (y, w2) in qps[1], (z, w3) in qps[2], (w, w4) in qps[1])
end

#Common Vertex Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonVertex6D)
	qps = method.qps
    sum(w1*w2*w3*w4*w5*w6*k6p_cv(integrand, η1, η2, η3, η4, ξ1, ξ2)
        for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (η4, w4) in qps, (ξ1, w5) in qps, (ξ2, w6) in qps)
end

#Common Vertex Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonVertex6D_S)
	qps = method.qps
    sum(w1*w2*k6p_cv_A(integrand, x[1], x[2], x[3], y[1], y[2], y[3])
        for (x, w1) in qps, (y, w2) in qps)
end

#Positive Distance Tensor-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance6D)
	qps = method.qps
	sum(w1*w2*w3*w4*w5*w6*k6p_pd(integrand, η1, η2, η3, η4, ξ1, ξ2)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (η4, w4) in qps, (ξ1, w5) in qps, (ξ2, w6) in qps)
end

#Positive Distance Simplex-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance6D_S)

	qps = method.qps
	sum(w1*w2*k6p_pd_A(integrand, x[1], x[2], x[3], y[1], y[2], y[3])
		for (x, w1) in qps, (y, w2) in qps)
end

#=Tetrahedron-triangle interaction =#
#Common Face Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonFace5D)
    qps = method.qps
	sum(w1*w2*w3*w4*w5*k5p_cf(integrand, x1, x2, x3, y1, y2)
		for (x1, w1) in qps, (x2, w2) in qps, (x3, w3) in qps, (y1, w4) in qps, (y2, w5) in qps)
end
#Common Face Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonFace5D_S)
    qps = method.qps
	sum(w1*w2*w3*k5p_cf_A(integrand, x[1], x[2], x[3], y, z)
		for (x, w1) in qps[3], (y, w2) in qps[1], (z, w3) in qps[1]) +
    sum(w1*w2*k5p_cf_B(integrand, x[1], x[2], x[3], y[1], y[2])
        for (x, w1) in qps[3], (y, w2) in qps[2])
end

#Common Edge Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonEdge5D)
    qps = method.qps
	sum(w1*w2*w3*w4*w5*k5p_ce(integrand, x1, x2, x3, y1, y2)
		for (x1, w1) in qps, (x2, w2) in qps, (x3, w3) in qps, (y1, w4) in qps, (y2, w5) in qps)
end

#Common Edge Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonEdge5D_S)
    qps = method.qps
	sum(w1*w2*w3*k5p_ce_A(integrand, x[1], x[2], y[1], y[2], z)
		for (x, w1) in qps[2], (y, w2) in qps[2], (z, w3) in qps[1]) +
    sum(w1*w2*k5p_ce_B(integrand, x[1], x[2], y[1], y[2], y[3])
		for (x, w1) in qps[2], (y, w2) in qps[3]) +
    sum(w1*w2*w3*k5p_ce_C(integrand, x[1], x[2], y[1], y[2], z)
		for (x, w1) in qps[2], (y, w2) in qps[2], (z, w3) in qps[1])
end

#Common Vertex Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonVertex5D)
    qps = method.qps
	sum(w1*w2*w3*w4*w5*k5p_cv(integrand, x1, x2, x3, y1, y2)
		for (x1, w1) in qps, (x2, w2) in qps, (x3, w3) in qps, (y1, w4) in qps, (y2, w5) in qps)
end

#Common Vertex Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonVertex5D_S)
    qps = method.qps
	sum(w1*w2*k5p_cv_A(integrand, x[1], x[2], x[3], y[1], y[2])
		for (x, w1) in qps[1], (y, w2) in qps[2])
end

#Positive Distance Tensor-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance5D)

    qps = method.qps
	sum(w1*w2*w3*w4*w5*k5p_pd(integrand, x1, x2, x3, y1, y2)
		for (x1, w1) in qps, (x2, w2) in qps, (x3, w3) in qps, (y1, w4) in qps, (y2, w5) in qps)
end

#Positive Distance Simplex-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance5D_S)

    qps = method.qps
	sum(w1*w2*k5p_pd_A(integrand, x[1], x[2], x[3], y[1], y[2])
		for (x, w1) in qps[1], (y, w2) in qps[2])
end

#=Triangle-triangle interaction =#
#Common Face Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonFace4D)

	qps = method.qps
	sum(w1*w2*w3*w4*k4p_cf(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end
#Common Face Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonFace4D_S)

	qps = method.qps
	sum(w1*w2*k4p_cf_A(integrand, x[1], x[2], x[3], y)
		for (x, w1) in qps[2], (y, w2) in qps[1])
end

#Common Edge Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonEdge4D)

	qps = method.qps
	sum(w1*w2*w3*w4*k4p_ce(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end

#Common Edge Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonEdge4D_S)

	qps = method.qps
	sum(w1*w2*w3*k4p_ce_A(integrand, x[1], x[2], y, z)
		for (x, w1) in qps[2],  (y, w2) in qps[1], (z, w3) in qps[1]) +
    sum(w1*w2*k4p_ce_B(integrand, x[1], x[2], y[1], y[2])
            for (x, w1) in qps[2], (y, w2) in qps[2])
end

#Common Vertex Tensor-Product
function sauterschwab_parameterized(integrand, method::CommonVertex4D)

	qps = method.qps
	sum(w1*w2*w3*w4*k4p_cv(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end

#Common Vertex Simplex-Product
function sauterschwab_parameterized(integrand, method::CommonVertex4D_S)

	qps = method.qps
	sum(w1*w2*k4p_cv_A(integrand, x[1], x[2], y[1], y[2])
		for (x, w1) in qps, (y, w2) in qps)
end

#Positive Distance Tensor-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance4D)

	qps = method.qps
	sum(w1*w2*w3*w4*k4p_pd(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end

#Positive Distance Simplex-Product
function sauterschwab_parameterized(integrand, method::PositiveDistance4D_S)

	qps = method.qps
	sum(w1*w2*k4p_pd_A(integrand, x[1], x[2], y[1], y[2])
		for (x, w1) in qps, (y, w2) in qps)
end

#=
function sauterschwab_parameterized(integrand, method::CommonEdgeS)

    integral = 0.0
    qps_1 = method.qps[1]
    qps_2 = method.qps[2]
    qps_3 = method.qps[3]
    qps_4 = method.qps[4]

	integral += sum(w1*w2*w3*k3p_ce_1(integrand, x[1], x[2], y[1], y[2], z[1], z[2])
		for (x, w1) in qps_2, (y, w2) in qps_2, (z, w3) in qps_2)

    integral += sum(w1*w2*w3*k3p_ce_2(integrand, x[1], x[2], y[1], y[2], y[3], z[1])
        for (x, w1) in qps_2, (y, w2) in qps_3, (z, w3) in qps_1)

    integral += sum(w1*w2*k3p_ce_3(integrand, x[1], x[2], y[1], y[2], y[3], y[4])
        for (x, w1) in qps_2, (y, w2) in qps_4)

    integral += sum(w1*w2*w3*k3p_ce_4(integrand, x[1], x[2], y[1], z[1], z[2], z[3])
        for (x, w1) in qps_2, (y, w2) in qps_1, (z, w3) in qps_3)

    integral += sum(w1*w2*w3*w4*k3p_ce_5(integrand, x[1], x[2], y[1], z[1], z[2], w[1])
        for (x, w1) in qps_2, (y, w2) in qps_1, (z, w3) in qps_2, (w, w4) in qps_1)

    return integral
end

function sauterschwab_parameterized(integrand, method::CommonVertexS)

    qps = method.qps
	sum(w1*w2*k3p_cv_s( integrand, z[1], z[2], z[3], x[1], x[2], x[3])
		for (x,w1) in qps, (z,w2) in qps)

end

function sauterschwab_parameterized(integrand, method::CommonEdgeS_a)

    integral = 0.0
    qps_1 = method.qps[1]
    qps_2 = method.qps[2]
    qps_3 = method.qps[3]
    qps_4 = method.qps[4]

	integral += sum(w1*w2*w3*k3p_ce_1_slow(integrand, x[1], x[2], y[1], y[2], z[1], z[2])
		for (x, w1) in qps_2, (y, w2) in qps_2, (z, w3) in qps_2)
    #=
    integral += sum(w1*w2*w3*k3p_ce_2(integrand, x[1], x[2], y[1], y[2], y[3], z[1])
        for (x, w1) in qps_2, (y, w2) in qps_3, (z, w3) in qps_1)

    integral += sum(w1*w2*k3p_ce_3(integrand, x[1], x[2], y[1], y[2], y[3], y[4])
        for (x, w1) in qps_2, (y, w2) in qps_4)

    integral += sum(w1*w2*w3*k3p_ce_4(integrand, x[1], x[2], y[1], z[1], z[2], z[3])
        for (x, w1) in qps_2, (y, w2) in qps_1, (z, w3) in qps_3)

    integral += sum(w1*w2*w3*w4*k3p_ce_5(integrand, x[1], x[2], y[1], z[1], z[2], w[1])
        for (x, w1) in qps_2, (y, w2) in qps_1, (z, w3) in qps_2, (w, w4) in qps_1)
     =#
    return integral
end

function sauterschwab_parameterized(integrand, method::CommonVertexGM)

	qps = method.qps
    sum(w1*w2*k3p_cv_gm(integrand, x[1], x[2], x[3], y[1], y[2], y[3])
        for (x, w1) in qps[1], (y, w2) in qps[2])
end

=#

#=
function sauterschwab_parameterized(integrand, method::New)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_cf_new(integrand, η1, η2, ξ1, ξ2)
		for (η1, w1) in qps, (η2, w2) in qps, (ξ1, w3) in qps, (ξ2, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::GM)

	qps = method.qps
	sum(w1*w2*k3p_cf_gm( integrand, y[1], y[2], y[3], x)
		for (x,w1) in qps[1], (y,w2) in qps[2])

end

function sauterschwab_parameterized(integrand, method::Orig2)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_cf_orig2(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::NewEdge)

	qps = method.qps
	sum(w1*w2*w3*w4*k3p_ce_new(integrand, η1, η2, η3, ξ1)
		for (η1, w1) in qps, (η2, w2) in qps, (η3, w3) in qps, (ξ1, w4) in qps)
end

function sauterschwab_parameterized(integrand, method::GMEdge)

	qps = method.qps
	sum(w1*w2*w3*k3p_ce_gm( integrand, z[1], z[2], x, y)
		for (x,w1) in qps[1], (y,w2) in qps[1], (z,w3) in qps[2])

end
=#
#=
function sauterschwab_parameterized(integrand, method::GMVertex)

	qps = method.qps
	sum(w1*w2*k3p_cv_gm( integrand, z[1], z[2], x[1], x[2])
		for (x,w1) in qps[1], (z,w2) in qps[2])

end
=#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct EncryptionSettingsElement <: SwaggerModel
    diskEncryptionKey::Any # spec type: Union{ Nothing, KeyVaultAndSecretReference } # spec name: diskEncryptionKey
    keyEncryptionKey::Any # spec type: Union{ Nothing, KeyVaultAndKeyReference } # spec name: keyEncryptionKey

    function EncryptionSettingsElement(;diskEncryptionKey=nothing, keyEncryptionKey=nothing)
        o = new()
        validate_property(EncryptionSettingsElement, Symbol("diskEncryptionKey"), diskEncryptionKey)
        setfield!(o, Symbol("diskEncryptionKey"), diskEncryptionKey)
        validate_property(EncryptionSettingsElement, Symbol("keyEncryptionKey"), keyEncryptionKey)
        setfield!(o, Symbol("keyEncryptionKey"), keyEncryptionKey)
        o
    end
end # type EncryptionSettingsElement

const _property_map_EncryptionSettingsElement = Dict{Symbol,Symbol}(Symbol("diskEncryptionKey")=>Symbol("diskEncryptionKey"), Symbol("keyEncryptionKey")=>Symbol("keyEncryptionKey"))
const _property_types_EncryptionSettingsElement = Dict{Symbol,String}(Symbol("diskEncryptionKey")=>"KeyVaultAndSecretReference", Symbol("keyEncryptionKey")=>"KeyVaultAndKeyReference")
Base.propertynames(::Type{ EncryptionSettingsElement }) = collect(keys(_property_map_EncryptionSettingsElement))
Swagger.property_type(::Type{ EncryptionSettingsElement }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_EncryptionSettingsElement[name]))}
Swagger.field_name(::Type{ EncryptionSettingsElement }, property_name::Symbol) =  _property_map_EncryptionSettingsElement[property_name]

function check_required(o::EncryptionSettingsElement)
    true
end

function validate_property(::Type{ EncryptionSettingsElement }, name::Symbol, val)
end


text:
# starting with julia
# https://dataframes.juliadata.org/stable/man/comparisons/
using Pkg; Pkg.instantiate(); Pkg.activate(".")
using RCall
using CSV
using DataFrames
using ShiftedArrays
# get the data

dat_raw = DataFrame(CSV.File("data/day02.txt" , header = 0))

rename!(dat_raw, [:direction, :steps])

direction_matrix = DataFrame(direction = ["forward", "up", "down"],
                            axis = [1,2,2],
                            value = [1,-1,1])

combined_data = leftjoin(dat_raw, direction_matrix, on = :direction)

combined_data.stepvalue = combined_data.value .* combined_data.steps

combine(groupby(combined_data, :axis), :stepvalue => sum)


text:
const known_dependencies = [:FunctionPtr, :RefCounter]

function resolve_type(type)
    if type isa Symbol
        type in extension_types ? :(vk.$type) : type
    else
        postwalk(type) do ex
            ex isa Symbol ? resolve_type(ex) : ex
        end
    end
end

resolve_types(ex) = postwalk(resolve_type, ex)

function field_deps(ex)
    @match ex begin
        :($_::$T) => innermost_type(T)
        _ => nothing
    end
end

function raw_dependencies(ex)
    p = deconstruct(ex)
    @match category(ex) begin
        :struct => begin
            # handle structs wrapped with @auto_hash_equals
            p = haskey(p, :macro) ? deconstruct(p[:decl]) : p
            field_deps.(p[:fields])
        end
        :function => vcat(map.(Ref(@λ(begin
            :($arg::$type) => inner_type(type)
            :(::$type) => inner_type(type)
            arg => nothing
        end)), [p[:args], p[:kwargs]])...)
        :const => p[:value] isa Symbol ? [p[:value]] : Symbol[]
        :enum => [p[:type]]
        :doc => raw_dependencies(p[:ex])
    end
end

function dependencies(ex)
    deps = raw_dependencies(ex)
    filter!.(
        [
            !isalias,
            x -> !startswith(string(x), r"(?:Vk|VK_|StdVideo)"),
            !is_vulkan_type,
            !in(known_dependencies),
            !in(extension_types),
            !in(Base.names(Core)),
            !in(Base.names(Base)),
            x -> x isa Symbol,
        ],
        Ref(deps),
    )
    deps
end

function check_is_dag(g, decls)
    if is_cyclic(g) || !is_directed(g)
        cycles = simplecycles_hadwick_james(g)
        inds = unique(vcat(cycles...))
        problematic_decls = getindex.(Ref(decls), inds)
        error(
            """
      Dependency graph is not a directed acyclic graph (is $(is_directed(g) ? "directed" : "undirected") and $(is_cyclic(g) ? "cyclic" : "acyclic"))

      $(is_cyclic(g) ? """
      Cycles:

      $(join(generate.(problematic_decls), "\n"^2))""" : "")
      """,
        )
    end
end

function resolve_dependencies(decls)
    defined_names = names.(decls)
    g = SimpleDiGraph(length(decls))

    for (j, decl) ∈ enumerate(decls)
        for dep ∈ dependencies(decl)
            i = findfirst(Base.Fix1(in, dep), defined_names)
            !isnothing(i) || error("Could not find dependency '$dep' for $decl")
            add_edge!(g, i, j)
        end
    end

    check_is_dag(g, decls)
    topological_sort_by_dfs(g)
end

function check_dependencies(decls)
    encountered_deps = []
    for decl ∈ decls
        deps = dependencies(decl)
        if deps ⊈ encountered_deps
            error("Unsolved dependencies $deps for declaration $decl")
        end
        append!(encountered_deps, names(decl))
    end
    decls
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction