text:
export GbmModel, GbmRegressor, GbmClassifier, GbmRanker

@with_kw mutable struct GbmModel <: BaseEstimator
    pyo::PyObject = PyNULL()
    name::String = "lightgbm"
    max_bin::Int = 16
    objective::String = "binary:logistic"
    max_depth::Int = 5
    min_child_weight::Float32 = 1
    gamma::Float32 = 0
    num_leaves::Int = 31
    min_child_samples::Int = 1000
    max_position::Int = 200
    label_gain::String = join(0:30, ',')
    n_estimators::Int = 50
    reg_lambda::Float32 = 0.1
    subsample::Float32 = 1
    colsample_bytree::Float32 = 1
    colsample_bylevel::Float32 = 1
    bagging_temperature::Float32 = 1
end

is_classifier(m::GbmModel) = occursin(r"binary|multi", m.objective)
is_ranker(m::GbmModel) = occursin("rank", m.objective)

function gridparams(m::GbmModel)
    @unpack objective, name = m
    rank = is_ranker(m)
    grid = [
        "n_estimators" => [20, 50, 100],
        "reg_lambda" => [0.1, 1, 10],
        [
            "name" => "lightgbm",
            "objective" => rank ? ["rank:ndcg"] : [objective],
            "num_leaves" => [15, 31],
            "max_position" => rank ? [100, 200, 500, 1000] : [1],
            "label_gain" => rank ? [join(0:30, ','), ""] : [""],
            "min_child_samples" => [1000, 100, 5000],
            "subsample" => [1, 0.8],
            "colsample_bytree" => [1, 0.8]
        ],
        [
            "name" => ["xgboost"],
            "objective" => rank ? ["rank:pairwise", "rank:map", "rank:ndcg"] : [objective],
            "max_depth" => [5, 3, 7],
            "min_child_weight" => [1, 10, 100],
            "gamma" => [0, 0.01, 0.1, 1],
            "subsample" => [1, 0.8],
            "colsample_bytree" => [1, 0.8],
        ],
        [
            "name" => ["catboost"],
            "objective" => !rank ? [objective] :
                ["rank:queryrmse", "rank:yetirank",
                "rank:yetirankpairwise", "rank:pairlogit",
                "rank:pairlogitpairwise"],
            "max_depth" => [5, 3, 7],
            "min_child_samples" => !usegpu() ? [1000] : [1000, 100, 5000],
            "bagging_temperature" => [0, 1, 2, 5, 10],
            "colsample_bylevel" => usegpu() ? [1] : [1, 0.8],
        ]
    ]
    return griparams(grid)
end

function fit!(m::GbmModel, x, y, w = nothing; group = nothing, columns = nothing)
    columns = something(columns, string.(1:size(x, 1)))
    if is_ranker(m) && isnothing(group) && ndims(x) == 3
        group = repeat([size(x, 2)], size(x, 3))
    end
    x, y, w, group = pymat(x), vec(y), vec(w), vec(group)
    x = DataFrame(x, columns = columns)
    @unpack name, max_bin, objective, max_depth, min_child_weight, gamma = m
    @unpack num_leaves, min_child_samples, max_position, label_gain = m
    @unpack max_depth, n_estimators, reg_lambda, subsample, colsample_bytree = m
    @unpack colsample_bylevel, bagging_temperature = m
    num_class = occursin("reg", objective) ? 1 : length(unique(y))
    objective = occursin("binary", objective) && num_class > 2 ? "multi:softmax" : objective
    num_class = occursin("multi", objective) ? num_class : 1
    y = occursin("binary", objective) && name == "xgboost" ? (y .+ 1) ./ 2 : y
    objective = denormalize_objective(objective, name)
    min_child_samples = min(length(y) ÷ 10, min_child_samples)
    if name == "xgboost"
        @from xgboost imports XGBModel, XGBRanker
        XGBModel = !is_ranker(m) ? XGBModel : XGBRanker
        tree_method = usegpu() && is_ranker(m) ? "gpu_hist" : "auto"
        pyo = XGBModel(n_jobs = 20, max_bin = max_bin, tree_method = tree_method, single_precision_histogram = true, n_gpus = countgpus(), num_class = num_class)
        pyo.set_params(;@NT(objective, max_depth, min_child_weight, n_estimators, reg_lambda, subsample, gamma)...)
        tree_method == "auto" && pyo.set_params(colsample_bytree = colsample_bytree)
        if isnothing(group)
            @time pyo.fit(x, y, sample_weight = w, eval_set = [(x, y)], sample_weight_eval_set = bra(w), verbose = true)
        else
            if length(group) != length(w)
                @imports numpy as np
                w = map(sum, np.split(w, cumsum(group[1:end - 1])))
                @show length(w), length(group)
            end
            @time pyo.fit(x, y, group, sample_weight = w, eval_set = [(x, y)], sample_weight_eval_set = bra(w), eval_group = bra(group), verbose = true)
        end
    elseif name == "lightgbm"
        @from lightgbm imports LGBMModel
        @from unidecode imports unidecode
        device = usegpu() ? "gpu" : "cpu"
        pyo = LGBMModel(n_jobs = 20, max_bin = max_bin - 1, device = device, label_gain = label_gain, num_class = num_class)
        pyo.set_params(;@NT(objective, num_leaves, min_child_samples, max_position, n_estimators, reg_lambda, subsample, colsample_bytree)...)
        x.columns = columns′ = [replace(filter(!isspace, unidecode(c)), ":" => "-") for c in columns]
        @time pyo.fit(x, y, w, group = group, eval_set = [(x, y)], eval_sample_weight = bra(w), eval_group = bra(group), feature_name = columns′)
    elseif name == "catboost"
        @from catboost imports CatBoost
        @imports numpy as np
        min_child_samples = usegpu() ? min_child_samples : 1
        colsample_bylevel = usegpu() ? 1 : colsample_bylevel
        task_type = usegpu() ? "GPU" : "CPU"
        devices = string("0:", countgpus() - 1)
        loss_function = objective
        pyo = CatBoost(pairs((thread_count = 20, max_bin = max_bin - 1, task_type = task_type, devices = devices, learning_rate = 0.1)))
        pyo.set_params(;@NT(loss_function, max_depth, n_estimators, reg_lambda, colsample_bylevel, min_child_samples, bagging_temperature)...)
        group = isnothing(group) ? nothing : np.repeat(axes(group, 1), group)
        @time pyo.fit(x, y, sample_weight = w, group_id = group)
    end
    @pack! m = pyo
    BSON.bson("model.bson", model = m)
    visualize(m, columns)
    return m
end

function predict_gbm(m::GbmModel, x)
    @unpack name, pyo = m
    if name == "xgboost"
        pyo.predict(pymat(x), output_margin = true, validate_features = false)
    elseif name == "lightgbm"
        pyo.predict(pymat(x), raw_score = true)
    else
        pyo.predict(pymat(x))
    end
end

function predict_proba(m::GbmModel, x)
    ŷ = predict_gbm(m, x)
    if ndims(ŷ) == 1
        ŷ = sigmoid.(ŷ)
        hcat(1 .- ŷ,  ŷ)
    else
        softmax(ŷ, dims = 2)
    end
end

function predict(m::GbmModel, x)
    ŷ = predict_gbm(m, x)
    !is_classifier(m) && return ŷ
    if ndims(ŷ) == 1
        signone.(ŷ)
    else
        [I[2] - 1 for I in argmax(ŷ, dims = 2)]
    end
end

istree(m::GbmModel) = true

function visualize(m::GbmModel, columns)
    @unpack pyo = m
    mname, tname = pymodulename(pyo), pytypename(pyo)
    plot_tree = parseenv("PLOT_TREE", false)
    graphs, w = [], []
    if occursin("lightgbm", mname)
        bst = tname == "Booster" ? pyo : pyo.booster_
        w = bst.feature_importance()
        @from lightgbm imports create_tree_digraph
        plot_tree && for n in 1:min(5, bst.num_trees())
            push!(graphs, create_tree_digraph(bst, tree_index = n - 1))
        end
    elseif occursin("xgboost", mname)
        bst = tname == "Booster" ? pyo : pyo.get_booster()
        score = bst.get_score()
        columns = collect(keys(score))
        w = collect(values(score))
        @from xgboost imports to_graphviz
        plot_tree && for n in 1:min(5, bst.best_iteration)
            push!(graphs, to_graphviz(bst, num_trees = n - 1))
        end
    else
        return
    end
    write_feaimpt(w, columns)
    plot_tree && for (n, g) in enumerate(graphs)
        try g.render("gbm/$n.pdf", cleanup = "true") catch e; println(e) end
    end
    plot_tree && try bash("pdfcat -o gbm.pdf gbm/*.pdf && rm -rf gbm") catch end
end

function save(m::GbmModel, dst = "model.gbm")
    @unpack pyo = m
    mname, tname = pymodulename(pyo), pytypename(pyo)
    if occursin("lightgbm", mname)
        bst = tname == "Booster" ? pyo : pyo.booster_
    elseif occursin("xgboost", mname)
        bst = tname == "Booster" ? pyo : pyo.get_booster()
    elseif occursin("catboost", mname)
        bst = pyo
    else
        return
    end
    bst.save_model(dst)
end

function denormalize_objective(objective, name)
    mapping = Dict(v => k for (k, v) in objective_mapping[name])
    get(mapping, objective, objective)
end

function normalize_objective(objective)
    mapping = merge(values(objective_mapping)...)
    get(mapping, objective, objective)
end

const objective_mapping = Dict(
    "xgboost" => Dict(
        "reg:squarederror" => "reg:mse",
        "reg:squaredlogerror" => "reg:msle",
        "reg:gamma" => "reg:gamma",
        "reg:tweedie" => "reg:tweedie",
        "reg:logistic" => "reg:logistic",
        "binary:logistic" => "binary:logistic",
        "binary:logitraw" => "binary:logitraw",
        "binary:hinge" => "binary:hinge",
        "count:poisson" => "reg:poisson",
        "survival:cox" => "reg:cox",
        "multi:softmax" => "multi:softmax",
        "multi:softprob" => "multi:softprob",
        "rank:pairwise" => "rank:pairwise",
        "rank:ndcg" => "rank:ndcg",
        "rank:map" => "rank:map"
    ),
    "lightgbm" => Dict(
        "regression" => "reg:mse",
        "regression_l1" => "reg:mae",
        "huber" => "reg:huber",
        "fair" => "reg:fair",
        "poisson" => "reg:poisson",
        "quantile" => "reg:quantile",
        "mape" => "reg:mape",
        "gamma" => "reg:gamma",
        "tweedie" => "reg:tweedie",
        "binary" => "binary:logistic",
        "cross_entropy" => "binary:cross_entropy",
        "cross_entropy_lambda" => "binary:cross_entropy_lambda",
        "multiclass" => "multi:softmax",
        "multiclassova" => "multi:ova",
        "lambdarank" => "rank:ndcg"
    ),
    "catboost" => Dict(
        "MAE" => "reg:mae",
        "MAPE" => "reg:mape",
        "Poisson" => "reg:poisson",
        "Quantile" => "reg:quantile",
        "RMSE" => "reg:mse",
        "LogLinQuantile" => "reg:loglinquantile",
        "Huber" => "reg:huber",
        "Expectile" => "reg:expectile",
        "Logloss" => "binary:logistic",
        "CrossEntropy" => "binary:crossentropy",
        "MultiClass" => "multi:softmax",
        "MultiClassOneVsAll" => "multi:ova",
        "PairLogit" => "rank:pairlogit",
        "PairLogitPairwise" => "rank:pairlogitpairwise",
        "YetiRank" => "rank:yetirank",
        "YetiRankPairwise" => "rank:yetirankpairwise",
        "QueryCrossEntropy" => "rank:querycrossentropy",
        "QueryRMSE" => "rank:queryrmse",
        "QuerySoftMax" => "rank:querysoftmax"
    )
)

GbmRegressor(;ka...) = GbmModel(;objective = "reg:mse", ka...)

GbmClassifier(;ka...) = GbmModel(;objective = "binary:logistic", ka...)

GbmRanker(;ka...) = GbmModel(;objective = "rank:ndcg", ka...)


text:
# Autogenerated wrapper script for Libffi_jll for aarch64-apple-darwin
export libffi

JLLWrappers.@generate_wrapper_header("Libffi")
JLLWrappers.@declare_library_product(libffi, "@rpath/libffi.6.dylib")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libffi,
        "lib/libffi.6.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# simple benchmars

using Luxor, BenchmarkTools

using BenchmarkTools

function bm1()
    Drawing(1000, 1000, "luxor-test1.pdf")
    origin()
    t = Tiler(1000, 1000, 10, 10, margin=10)
    for (pos, n) in t
        randomhue()
        circle(pos, t.tilewidth/2, :fill)
    end
    finish()
end

function bm2()
    Drawing(1000, 1000, "luxor-test2.pdf")
    origin()
    t = Tiler(1000, 1000, 10, 10, margin=10)
    for (pos, n) in t
        setcolor"red"
        setopacity(rand())
        box(pos, t.tilewidth/2, t.tilewidth/2, :fill)
    end
    finish()
end

@benchmark bm1()

@benchmark bm2()

#=

1: Julia v1.0.0 Luxor v1.0.0 iMac 2018-10-05 12:37:11

julia-1.0> @benchmark bm1()
BenchmarkTools.Trial:
  memory estimate:  75.00 KiB
  allocs estimate:  2886
  --------------
  minimum time:     1.975 ms (0.00% GC)
  median time:      2.151 ms (0.00% GC)
  mean time:        2.271 ms (1.59% GC)
  maximum time:     56.407 ms (93.63% GC)
  --------------
  samples:          2196
  evals/sample:     1

julia-1.0> @benchmark bm2()
BenchmarkTools.Trial:
  memory estimate:  131.19 KiB
  allocs estimate:  4689
  --------------
  minimum time:     830.662 μs (0.00% GC)
  median time:      1.049 ms (0.00% GC)
  mean time:        1.169 ms (2.81% GC)
  maximum time:     61.274 ms (97.08% GC)
  --------------
  samples:          4247
  evals/sample:     1

2:  Julia v1.0.0 Luxor v1.1.0 iMac 2018-10-05 12:40:11

after some type work, there's a slight improvement in allocations

  julia-1.0> @benchmark bm1()
  BenchmarkTools.Trial:
    memory estimate:  38.94 KiB
    allocs estimate:  782
    --------------
    minimum time:     1.950 ms (0.00% GC)
    median time:      2.136 ms (0.00% GC)
    mean time:        2.302 ms (1.57% GC)
    maximum time:     68.908 ms (92.70% GC)
    --------------
    samples:          2165
    evals/sample:     1

  julia-1.0> @benchmark bm2()
  BenchmarkTools.Trial:
    memory estimate:  92.00 KiB
    allocs estimate:  2285
    --------------
    minimum time:     663.034 μs (0.00% GC)
    median time:      736.955 μs (0.00% GC)
    mean time:        806.080 μs (2.46% GC)
    maximum time:     49.676 ms (96.68% GC)
    --------------
    samples:          6164
    evals/sample:     1

3:  Julia v1.1 Luxor v1.3.0-master iMac 2019-03-03 17:30:20

  julia-1.1> @benchmark bm1()
    BenchmarkTools.Trial:
      memory estimate:  39.50 KiB
      allocs estimate:  787
      --------------
      minimum time:     1.917 ms (0.00% GC)
      median time:      2.145 ms (0.00% GC)
      mean time:        2.247 ms (1.18% GC)
      maximum time:     47.751 ms (91.55% GC)
      --------------
      samples:          2217
      evals/sample:     1

    julia-1.1> @benchmark bm2()
    BenchmarkTools.Trial:
      memory estimate:  50.27 KiB
      allocs estimate:  1190
      --------------
      minimum time:     623.612 μs (0.00% GC)
      median time:      781.187 μs (0.00% GC)
      mean time:        856.280 μs (1.62% GC)
      maximum time:     48.690 ms (95.28% GC)
      --------------
      samples:          5793
      evals/sample:     1

 =#


text:
__precompile__()
"""
Module for defining and checking properties.
"""
module Properties

using LazySets

import LazySets.dim

# ==============================
# Property struct and evaluation
# ==============================
include("Property.jl")
export Property,
       dim,
       check

include("Conjunction.jl")
include("Disjunction.jl")
export Conjunction,
       Disjunction

include("BadStatesProperty.jl")
export BadStatesProperty

include("SafeStatesProperty.jl")
export SafeStatesProperty

end # module


text:
# -*- coding: utf-8 -*-
# + {}
module ExpStateModels
import Base: show, getindex, setindex!, put!, replace!

using SemanticModels.Parsers
using SemanticModels.ModelTools
import SemanticModels.ModelTools: model

export ExpStateModel, ExpStateTransition, model

"""    ExpStateModel

represents an agent based model symbolically, with a collection of stats, agents, and transitions.
The agents are a collection of states and the transitions map from state to state. This structure allows you to represent an agent based model at the semantic level and apply transformations to that model.

The common transformations for an agent based model are adding, removing, or replacing, states, agents, or transitions. See [`model`](@ref), [`put!`](@ref), [`setindex`](@ref), [`getindex`](@ref), [`put!`](@ref), [`replace!`](@ref).
"""
struct ExpStateModel <: AbstractModel
    expr::Expr
    states
    agents
    transitions
end

"""    ExpStateTransiton

represents a state-transition function for an agent based model. `ExpStateTransition(s,x)` represents the transition from state `s` to any other states. The expression `x` should define a function that takes any number of arguments and returns a value representing the new state.
"""
struct ExpStateTransition
  state::Symbol
  expr::Expr
end

function model(::Type{ExpStateModel}, expr::Expr)
    constructor = callsites(expr, :StateModel)[1]
    states = constructor.args[2]
    #TODO: find out why this identifies extra assignments possibly from the (i,a) = enumerate(sm.agents) part.
    agents = findassign(expr, constructor.args[3])
    transitions = findassign(expr, constructor.args[4])
    return ExpStateModel(expr, states, agents, transitions)
end

function show(io::IO, m::ExpStateModel)
    write(io, "ExpStateModel(\n  states=$(repr(m.states)),\n  agents=$(repr(m.agents)),\n  transitions=$(repr(m.transitions))\n)")
end

""" put!(m::ExpStateModel, transition::ExpStateTransition)

Store an [`ExpStateTransition`](@ref) `transition` into an [`ExpStateModel`](@ref) `m`. A `put!` on an already set state throws an `Exception`.

Returns the transition function inserted into `m`.
"""
function put!(m::ExpStateModel, transition::ExpStateTransition)
  q = QuoteNode(transition.state)
  any(x->x.value==transition.state, m.states.args) && error("Symbol $(transition.state) already exists")
  any(x->typeof(x) == Expr && x.args[2] == q, m.transitions[1].args[2].args) && error("Symbol $sym has a transition, but is not in the states")
  push!(m.states.args, QuoteNode(transition.state))
  push!(m.transitions[1].args[2].args, :($q=>$(transition.expr)))
  return transition.expr
end

""" replace!(m::ExpStateModel, transition::ExpStateTransition)

Store an [`ExpStateTransition`](@ref) `transition` into an [`ExpStateModel`(@ref) `m`. A `replace!` on an already set state will replace the current transition function.

Returns the transition function inserted into `m`.
"""
function replace!(m::ExpStateModel, transition::ExpStateTransition)
  !any(x->x.value==transition.state, m.states.args) && error("Symbol $sym doesn't exist")
  found=filter(x->typeof(x) == Expr && x.args[2] == QuoteNode(transition.state), m.transitions[1].args[2].args)
  found[1].args[3] = transition.expr
  return transition.expr
end

""" setindex!(m::ExpStateModel, expr::Expr, sym::Symbol)

Store an [`Expr`](@ref) `expr` as the transition function for state `sym` in [`ExpStateModel`](@ref) `m`. A `setindex!` on an already set state will replace the current transition function.

Returns the transition function inserted into `m`.
"""
function setindex!(m::ExpStateModel, expr::Expr, sym::Symbol)
  any(x->x.value==sym, m.states.args) && return replace!(m, ExpStateTransition(sym, expr))
  return put!(m, ExpStateTransition(sym, expr))
end

""" getindex(m::ExpStateModel, sym::Symbol)
Returns the transition function for state `sym` in [`ExpStateModel`](@ref) `m`. A `getindex` on a state that doesn't exist in `m` throws an `Exception`.
"""
function getindex(m::ExpStateModel, sym::Symbol)
  found=filter(x->typeof(x) == Expr && x.args[2] == QuoteNode(sym), m.transitions[1].args[2].args)
  length(found) == 0 && error("Symbol $sym doesn't exist")
  return found[1].args[3]
end

end


text:
function vehicle_paths(x_coefs::Vector{N}, y_coefs::Vector{N},
        s_max::AbstractVector; T=Float64) where {N <: AbstractVector}
    """
    Take in path coefficients as they are passed from python and convert to a list of VehiclePaths.
    """

    @assert length(x_coefs) == length(y_coefs) == length(s_max) "nplayers path mismatch"
    @assert length(x_coefs[1]) == length(y_coefs[1]) "ncoefs path mismatch"

    P = length(x_coefs)
    D = length(x_coefs[1])
    paths = [VehiclePath(SVector{D,T}(x_coefs[i]), SVector{D}(y_coefs[i]), s_max[i]) for i = 1:P]
    return paths
end

function vehicle_paths(x_coefs::AbstractMatrix, y_coefs::AbstractMatrix,
        s_max::AbstractVector; T=Float64)
    """
    Take in path coefficients as they are passed from python and convert to a list of VehiclePaths.
    """

    @assert size(x_coefs)[1] == size(y_coefs)[1] == length(s_max) "nplayers path mismatch"
    @assert size(x_coefs)[2] == size(y_coefs)[2] "ncoefs path mismatch"

    P,D = size(x_coefs)
    paths = [VehiclePath(SVector{D,T}(x_coefs[i,:]), SVector{D}(y_coefs[i,:]), s_max[i]) for i = 1:P]
    return paths
end

function initial_state(paths::Vector{VehiclePath{M}}, s_sdot::Vector{N}) where {N <: AbstractVector, M}
    """
    Take in initial s and sdot as they are passed from python and convert to an initial state
    """
    @assert length(paths) == length(s_sdot) "nplayers initial state mismatch"
    P = length(paths)
    T = eltype(paths[1].s_max)

    s = SVector{P,T}([s_sdot[i][1] for i=1:P])
    v = SVector{P,T}([s_sdot[i][2] for i=1:P])
    xs = SVector{P,T}([powerseries(paths[i].x_coefs, s[i]) for i=1:P])
    ys = SVector{P,T}([powerseries(paths[i].y_coefs, s[i]) for i=1:P])
    return SVector{4P,T}(xs..., ys..., s..., v...)
end

function initial_state(paths::Vector{VehiclePath{M}}, s_sdot::AbstractMatrix) where {M}
    """
    Take in initial s and sdot as they are passed from python and convert to an initial state
    """
    @assert length(paths) == size(s_sdot)[1] "nplayers initial state mismatch"
    P = length(paths)
    T = eltype(paths[1].s_max)

    s = SVector{P,T}([s_sdot[i,1] for i=1:P])
    v = SVector{P,T}([s_sdot[i,2] for i=1:P])
    xs = SVector{P,T}([powerseries(paths[i].x_coefs, s[i]) for i=1:P])
    ys = SVector{P,T}([powerseries(paths[i].y_coefs, s[i]) for i=1:P])
    return SVector{4P,T}(xs..., ys..., s..., v...)
end

function final_state(paths::Vector{VehiclePath{M}}) where {M}
    """
    Generate final state as s=paths.smax and v=0
    """
    P = length(paths)
    T = eltype(paths[1].s_max)
    s = SVector{P}([path.s_max for path in paths])
    v = zeros(SVector{P,T})

    xf = [SVector{4,T}(powerseries(paths[i].x_coefs, paths[i].s_max),
            powerseries(paths[i].y_coefs, paths[i].s_max),
            paths[i].s_max, 0.) for i = 1:P]
    return xf
end

function setup_problem(ssdot, xpoly, ypoly, smax, uc::Int; N::Int=20,
        dt::AbstractFloat=0.1, q::AbstractFloat=10., r::AbstractFloat=0.1,
        a_lower::AbstractFloat=-4.5, a_upper::AbstractFloat=1.5, radius::AbstractFloat = 4.,
        vel_control::Bool=false, vel_des::AbstractFloat=Inf,
        solver_opts::Options=Options(inner_print=false,outer_print=false))
    """
    Generate roundabout problem
    """

    # Set Up RoundaboutGame
    T = Float64
    paths = vehicle_paths(xpoly, ypoly, 1.0*smax; T=T)
    x0 = initial_state(paths, ssdot)

    model = RoundaboutGame(paths, uc)
    p, n, m = model.p, model.n, model.m

    # Define Problem Size
    probsize = ProblemSize(N, model)

    # Define the objective of each player
    # We use a LQR cost

    # State Cost
    if vel_control # velocity target
        @assert 0. <= vel_des < Inf "desired velocity invalid"
        xf = [SVector{4,T}(0., 0., 0., vel_des) for i=1:p]
        Q = [Diagonal(q*SVector{model.ni[i],T}(0,0,0,1)) for i=1:p] # Quadratic state cost
    else # final state target
        xf = [SVector{4,T}(0., 0., paths[i].s_max, 0.) for i=1:p]
        Q = [Diagonal(q*SVector{model.ni[i],T}(0,0,1,0)) for i=1:p] # Quadratic state cost
    end

    # Control Cost
    uf = [zeros(SVector{model.mi[i],T}) for i=1:p]
    R = [Diagonal(r*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost

    # Objectives of the game
    game_obj = GameObjective(Q,R,xf,uf,N,model)

    # Define the constraints that each player must respect
    game_con = GameConstraintValues(probsize)

    # Add collision avoidance
    add_collision_avoidance!(game_con, radius)

    # Add control bounds
    u_max = [a_upper*ones(SVector{m-model.uc,T}); zeros(SVector{model.uc,T})]
    u_min = [a_lower*ones(SVector{m-model.uc,T}); zeros(SVector{model.uc,T})]

    add_control_bound!(game_con, u_max, u_min)

    # Add velocity bounds ## TODO

    # Add control constraint on final $uc players ## TODO (currently doing with u_max=u_min=0)

    # Define the game problem
    prob = GameProblem(N,dt,x0,model,solver_opts,game_obj,game_con);
    return prob
end

max_stars_repo_name: lassepe/DecNashPlanning

max_issues_repo_name: lassepe/DecNashPlanning

max_forks_repo_name: lassepe/DecNashPlanning


text:
@testset "Shared catchability" begin
    run_sims(SharedQSpec, 1;
             prep_file = "prep.h5",
             csv = true,
             feather = true)

    # Construct all of the expected output file names
    output_files = joinpath.("sharedq",
                             "repl_01",
                             "sharedq_" .*
                             string.(1:length(sim_values(SharedQSpec)),
                                     pad = 2) .*
                             ["_popstate.h5" "_popstate.csv" "_catch.csv" "_popstate.feather" "_catch.feather"])

    @test isdir("sharedq")
    @test isdir("sharedq", "repl_01")
    @test all(isfile.(output_files))
end


text:
################################################################################
# Dynamics and simulation generation
################################################################################
using Pkg
pkg_dir = joinpath(dirname(pathof(ContactImplicitMPC)), "..")
Pkg.activate(pkg_dir)
include(joinpath(pkg_dir, "src/code_gen_loader.jl"))

# Generate dynamics
dynamicsdir = joinpath(dirname(pathof(ContactImplicitMPC)), "..", "src", "dynamics")
include(joinpath(dynamicsdir, "generate_dynamics.jl"))

# Generate simulation
simulationdir = joinpath(dirname(pathof(ContactImplicitMPC)), "..", "src", "simulation")
include(joinpath(simulationdir, "generate_simulation.jl"))


text:
using Plots
include("model.jl")
include("second_order_cone.jl")
include("step.jl")
include("simulate.jl")
include("visualize.jl")

# model
model = Particle(2 * 3, 3, 0, 1.0, 9.81, 0.5, 3, 3)
joinpath(@__DIR__,"model.jl")

# horizon
T = 1001

# time step
h = 0.01
t = range(0, stop = h * (T - 1), length = T)

# initial conditions
v1 = [10.0; 1.0; 10.0]
q1 = [0.0; 0.0; 1.0]

v2 = v1 - gravity(model, q1) * h
q2 = q1 + 0.5 * (v1 + v2) * h

# # simulate
q_sol, y_sol, b_sol, Δq1, Δq2, Δu1 = simulate(q1, q2, T, h)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

visualize!(vis, model,
    q_sol,
    Δt = h)

settransform!(vis["/Cameras/default"],
    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(pi))))

settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, 0.0, 3.0),LinearMap(RotY(-pi/2.5))))

## trajectory optimization
T = 26
h = 0.1
ul, uu = control_bounds(model, T, zeros(model.m), zeros(model.m))
ul[1] = [-Inf; -Inf; 0.0]
uu[1] = [Inf; Inf; 0.0]
x1 = [q1; q2]
xl, xu = state_bounds(model, T, x1 = x1)
qT = [1.0; 2.0; 0.0]

# Objective
obj = quadratic_tracking_objective(
        [t < T ? Diagonal(1.0e-2 * ones(model.n)) : 100.0 * Diagonal(ones(model.n)) for t = 1:T],
        [Diagonal(1.0e-3 * ones(model.m)) for t = 1:T-1],
        [[qT; qT] for t = 1:T],
        [zeros(model.m) for t = 1:T])

# Constraints
include("simulator_dynamics.jl")
con_sim = simulator_constraints(model, T)

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con_sim,
			   dynamics = false)

# Trajectory initialization
x0 = configuration_to_state([q1, linear_interpolation(q1, qT, T)...]) #linear_interpolation(x1, x1, T) # linear interpolation on state
u0 = [1.0e-1 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem
include_snopt()
@time z̄, info = solve(prob, copy(z0),
	nlp = :SNOPT7,
	tol = 1.0e-3, c_tol = 1.0e-3)

x̄, ū = unpack(z̄, prob)
q̄ = state_to_configuration(x̄)
t = range(0, stop = h * (T - 1), length = T)

plot(t, hcat(q̄[2:end]...)', xlabel = "time (s)", ylabel = "")
plot(t, hcat(ū...,ū[end])',
	xlabel = "time (s)" , ylabel = "", linetype = :steppost)

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

visualize!(vis, model,
    [[q̄[1] for i = 1:10] ..., q̄...],
    Δt = h)

@show q̄[end]

settransform!(vis["/Cameras/default"],
    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(-pi))))

settransform!(vis["/Cameras/default"],
	compose(Translation(0.0, 0.0, 1.0),LinearMap(RotY(-pi/2.5))))

setobject!(vis["particle_goal"],
	Rect(Vec(0, 0, 0),Vec(0.2, 0.2, 0.2)),
	MeshPhongMaterial(color = RGBA(0.0, 1.0, 1.0, 0.5)))

settransform!(vis["particle_goal"], Translation(qT...))

open(vis)


text:
using WECOrules
using Test
using DataFrames, CSV, Dates

@testset "WECOrules.jl" begin
    test = CSV.read(joinpath(@__DIR__,"tests.csv"),DataFrame)

    res = WECOrules.weco(test[:,[:Date, :Value]], 10.0, 1.0)

    @test all(test[:,:Rule1].==res[:,:Rule1])

    @test all(ismissing.(res[1:8,:Rule2]))
    @test all(test[9:end,:Rule2].==res[9:end,:Rule2])

    @test all(ismissing.(res[1:6,:Rule3]))
    @test all(test[7:end,:Rule3] .== res[7:end, :Rule3])

    @test all(ismissing.(res[1:15,:Rule4]))
    @test all(test[16:end,:Rule4] .== res[16:end, :Rule4])

    @test all(ismissing.(res[1:2,:Rule5]))
    @test all(test[3:end,:Rule5].==res[3:end,:Rule5])

    @test all(ismissing.(res[1:4,:Rule6]))
    @test all(test[5:end,:Rule6].==res[5:end,:Rule6])

    @test all(ismissing.(res[1:14,:Rule7]))
    @test all(test[15:end,:Rule7].==res[15:end,:Rule7])

    @test all(ismissing.(res[1:7,:Rule8]))
    @test all(test[8:end,:Rule8].==res[8:end,:Rule8])
end


text:
# Word Problems Logical Reasoning Tests
# http://www.graduatewings.co.uk/word-problems-logical-tests

# Questions 11-17

"""
A cruise line is scheduling seven week-long voyages for the ship Freedom. Each
voyage will occur in exactly one of the first seven weeks of the season:
weeks 1 through 7. Each voyage will be to exactly one of four destinations:
Guadeloupe, Jamaica, Martinique, or Trinidad. Each destination will be scheduled
for at least one of the weeks. The following conditions apply to Freedom’s schedule:
"""
julia> using AbstractLogic
Start the repl in command prompt by typing `=`.

# exactly one of four destinations: Guadeloupe, Jamaica, Martinique, or Trinidad
w1, w2, w3, w4, w5, w6, w7  ∈ Guadeloupe, Jamaica, Martinique, Trinidad [clear]

# Each destination will be scheduled for at least one of the weeks.
{{i}} = 'Guadeloupe' {{1,}}
{{i}} = 'Jamaica' {{1,}}
{{i}} = 'Martinique' {{1,}}
{{i}} = 'Trinidad' {{1,}}

# Jamaica will not be its destination in week 4.
w4 != 'Jamaica'

# Trinidad will be its destination in week 7.
w7 = 'Trinidad'

# Freedom will make exactly two voyages to Martinique,
{{i}} = 'Martinique' {{2}}

# and at least one voyage to Guadeloupe will occur in some week between those two voyages.
# Unfortunately this can be a bit tedious.
{{<i}} = 'Martinique' &&& {{>i}} = 'Martinique' &&& {{i}} = 'Guadeloupe' {{1,}}

# Guadeloupe will be its destination in the week preceding any voyage it makes to Jamaica.
{{i!}} = 'Jamaica' ==> {{i-1}} = 'Guadeloupe'

# No destination will be scheduled for consecutive weeks.
{{i}} != {{i+1}}

export lsat4

######
# 11. Which one of the following is an acceptable schedule of destinations for
# Freedom, in order from week 1 through week 7?

# (A) Guadeloupe, Jamaica, Martinique, Trinidad, Guadeloupe, Martinique, Trinidad
# I'll just check the first 4 weeks at first
check w1 = 'Guadeloupe'; w2 = 'Jamaica'; w3 = 'Martinique'; w4 = 'Trinidad'; w5 = 'Guadeloupe'; w6 = 'Martinique'; w7 = 'Trinidad'
# true, 1 out of 50 possible combinations 'true'.

# (B) Guadeloupe, Martinique, Trinidad, Martinique, Guadeloupe, Jamaica, Trinidad
check w1 = 'Guadeloupe'; w2 = 'Martinique'; w3 = 'Trinidad'; w4 = 'Martinique'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'
# false, 0 out of 50 possible combinations 'true'.

# (C) Jamaica, Martinique, Guadeloupe, Martinique, Guadeloupe, Jamaica, Trinidad
check w1 = 'Jamaica'; w2 = 'Martinique'; w3 = 'Guadeloupe'; w4 = 'Martinique'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Trinidad'
# false, 0 out of 50 possible combinations 'true'.

# (D) Martinique, Trinidad, Guadeloupe, Jamaica, Martinique, Guadeloupe, Trinidad
check w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Jamaica'; w5 = 'Martinique'; w6 = 'Guadeloupe'; w7 = 'Trinidad'
# false, 0 out of 50 possible combinations 'true'.

# (E) Martinique, Trinidad, Guadeloupe, Trinidad, Guadeloupe, Jamaica, Martinique
check w1 = 'Martinique'; w2 = 'Trinidad'; w3 = 'Guadeloupe'; w4 = 'Trinidad'; w5 = 'Guadeloupe'; w6 = 'Jamaica'; w7 = 'Martinique'
# false, 0 out of 50 possible combinations 'true'.

######
# 12. Which one of the following CANNOT be true about Freedom’s schedule of voyages?
# (A) Freedom makes a voyage to Trinidad in week 6
check w6 = 'Trinidad'
# false, 0 out of 50 possible combinations 'true'.  - Answer

# (B) Freedom makes a voyage to Martinique in week 5.
check w5 = 'Martinique'
# possible,  11 out of 25 possible combinations 'true'.

# (C) Freedom makes a voyage to Jamaica in week 6.
check w6 = 'Jamaica'
# possible,  3 out of 25 possible combinations 'true'.

# (D) Freedom makes a voyage to Jamaica in week 3.
check w3 = 'Jamaica'
# possible,  8 out of 25 possible combinations 'true'.

# (E) Freedom makes a voyage to Guadeloupe in week 3
check w3 = 'Guadeloupe'
# possible,  4 out of 25 possible combinations 'true'.

######
# 13. If Freedom makes a voyage to Trinidad in week 5, which one of the following
# could be true?
w5 = 'Trinidad'

#  (A) Freedom makes a voyage to Trinidad in week 1.
check w1 = 'Trinidad'
#false, 0 out of 3 possible combinations 'true'.

# (B) Freedom makes a voyage to Martinique in week 2.
check w2 = 'Martinique'
# false, 0 out of 3 possible combinations 'true'.

# (C) Freedom makes a voyage to Guadeloupe in week 3.
check w3 = 'Guadeloupe'
# false, 0 out of 3 possible combinations 'true'.

# (D) Freedom makes a voyage to Martinique in week 4.
check w4 = 'Martinique'
# possible,  1 out of 3 possible combinations 'true'.

# (E) Freedom makes a voyage to Jamaica in week 6.
check w6 = 'Jamaica'
# false, 0 out of 3 possible combinations 'true'.

import lsat4

######
# 14. If Freedom makes a voyage to Guadeloupe in week 1 and a voyage to Jamaica in
# week 5, which one of the following must be true?
w1 = 'Guadeloupe'; w5 = 'Jamaica'

# (A) Freedom makes a voyage to Jamaica in week 2.
check w2 = 'Jamaica'
# possible,  1 out of 3 possible combinations 'true'.

# (B) Freedom makes a voyage to Trinidad in week 2.
check w2 = 'Trinidad'
# possible,  1 out of 3 possible combinations 'true'.

# (C) Freedom makes a voyage to Martinique in week 3.
check w3 = 'Martinique'
# possible,  2 out of 3 possible combinations 'true'.

# (D) Freedom makes a voyage to Guadeloupe in week 6.
check w6 = 'Guadeloupe'
# false, 0 out of 3 possible combinations 'true'.

# (E) Freedom makes a voyage to Martinique in week 6.
check w6 = 'Martinique'
#true, 3 out of 3 possible combinations 'true'.
import lsat4

######
# 15. If Freedom makes a voyage to Guadeloupe in week 1 and to Trinidad in
# week 2, which one of the following must be true?
w1 = 'Guadeloupe'; w2 = 'Trinidad'

# (A) Freedom makes a voyage to Martinique in week 3.
prove w3 = 'Martinique'
# true, 1 out of 1 possible combinations 'true'.

# (B) Freedom makes a voyage to Martinique in week 4.
prove w4 = 'Martinique'
# false, 0 out of 1 possible combinations 'true'.

# (C) Freedom makes a voyage to Martinique in week 5.
prove w5 = 'Martinique'
# false, 0 out of 1 possible combinations 'true'.

# (D) Freedom makes a voyage to Guadeloupe in week 3.
prove w3 = 'Guadeloupe'
# false, 0 out of 1 possible combinations 'true'.

# (E) Freedom makes a voyage to Guadeloupe in week 5
prove w5 = 'Guadeloupe'
# false, 0 out of 1 possible combinations 'true'.
import lsat4

######
# 16. If Freedom makes a voyage to Martinique in week 3, which one of the
# following could be an accurate list of Freedom’s destinations in week 4
# and week 5,respectively?
w3 = 'Martinique'

# (A) Guadeloupe, Trinidad
any w4 = 'Guadeloupe'; w5 = 'Trinidad'
# true,  1 out of 8 possible combinations 'true'.

# (B) Jamaica, Guadeloupe
any w4 = 'Jamaica'; w5 = 'Guadeloupe'
# false, 0 out of 8 possible combinations 'true'.

# (C) Martinique, Trinidad
any w4 = 'Martinique'; w5 = 'Trinidad'
# false, 0 out of 8 possible combinations 'true'.

# (D) Trinidad, Jamaica
any w4 = 'Trinidad'; w5 = 'Jamaica'
# false, 0 out of 8 possible combinations 'true'.

# (E) Trinidad, Martinique
any w4 = 'Trinidad'; w5 = 'Martinique'
# false, 0 out of 8 possible combinations 'true'.
import lsat4

######
# 17. Which one of the following must be true about Freedom’s schedule of voyages?

# (A) Freedom makes a voyage to Guadeloupe either in week 1 or else in week 2.
check w1|w2 = 'Guadeloupe'
# false, 22 out of 25 possible combinations 'true'.
w1,w2 != 'Guadeloupe'

# (B) Freedom makes a voyage to Martinique either in week 2 or else in week 3.
check w2|w3 = 'Martinique'
# false, 12 out of 25 possible combinations 'true'.

# (C) Freedom makes at most two voyages to Guadeloupe.
check {{i}} = 'Guadeloupe' {{,2}}
# false, 20 out of 25 possible combinations 'true'.

# (D) Freedom makes at most two voyages to Jamaica.
check {{i}} = 'Jamaica' {{,2}}
# true, 25 out of 25 possible combinations 'true'.

# (E) Freedom makes at most two voyages to Trinidad.
check {{i}} = 'Trinidad' {{,2}}
# false, 24 out of 25 possible combinations 'true'.


text:
using Test
using LogicCircuits

module TestNodes

    using Test
    using ..LogicCircuits

    mutable struct TestINode <: Dag
        id::Int
        children::Vector{Dag}
        data
        counter::UInt32
        TestINode(i,c) = new(i,c,nothing,false)
    end

    mutable struct TestLNode <: Dag
        id::Int
        data
        counter::UInt32
        TestLNode(i) = new(i,nothing,false)
    end

    LogicCircuits.NodeType(::Type{<:TestINode}) = Inner()
    LogicCircuits.NodeType(::Type{<:TestLNode}) = Leaf()
    LogicCircuits.children(n::TestINode) = n.children

    @testset "Graphs utils for TestNodes" begin

        l1 = TestLNode(1)
        l2 = TestLNode(2)

        @test !has_children(l1)
        @test num_children(l1) == 0
        @test isleaf(l1)
        @test !isinner(l1)

        i1 = TestINode(1,[l1])
        i2 = TestINode(2,[l2])
        i12 = TestINode(3,[l1,l2])

        @test has_children(i1)
        @test has_children(i12)

        @test num_children(i1) == 1
        @test num_children(i12) == 2

        j1 = TestINode(1,[i1,i12])
        j2 = TestINode(2,[i2])
        j12 = TestINode(3,[i1,i2])

        r = TestINode(5,[j1,j2,j12])

        @test has_children(r)
        @test num_children(r) == 3

        reset_counter(r,5)
        @test r.counter == 5
        @test l1.counter == 5
        @test i12.counter == 5

        reset_counter(r)
        @test r.counter == 0
        @test l1.counter == 0
        @test i12.counter == 0

        foreach(r) do n
            n.id += 1
        end
        @test l1.id == 2
        @test l2.id == 3
        @test i12.id == 4
        @test j2.id == 3
        @test r.id == 6
        @test r.counter == 0
        @test l1.counter == 0
        @test i12.counter == 0

        foreach(r, l -> l.id += 1, i -> i.id -= 1)
        @test l1.id == 2+1
        @test l2.id == 3+1
        @test i12.id == 4-1
        @test j2.id == 3-1
        @test r.id == 6-1
        @test r.counter == 0
        @test l1.counter == 0
        @test i12.counter == 0

        @test filter(n -> iseven(n.id), r) == [l2,i2,j2]

        lastvisited = nothing
        foreach(n -> lastvisited=n,r)
        @test lastvisited === r

        lastvisited = nothing
        foreach_down(n -> lastvisited=n,r)
        @test isleaf(lastvisited)

        @test num_nodes(r) == 9
        @test num_edges(r) == 12

        @test isempty(inodes(l1))
        @test leafnodes(l1) == [l1]

        @test issetequal(inodes(r), [i1,i2,i12,j1,j2,j12,r])
        @test issetequal(innernodes(r), [i1,i2,i12,j1,j2,j12,r])
        @test issetequal(leafnodes(r), [l1,l2])

        @test tree_num_edges(r) == 14 # unverified

        @test linearize(r)[end] == r
        @test linearize(r)[1] == l1 || linearize(r)[1] == l2
        @test linearize(l2) == [l2]
        @test length(linearize(i12)) == 3

        @test eltype(linearize(r)) == Dag
        @test eltype(linearize(l1)) == TestLNode
        @test eltype(linearize(r, Any)) == Any

        @test left_most_descendent(r) == l1
        @test right_most_descendent(r) == l2

    end

end

@testset "Node stats tests" begin

    lc = load_logic_circuit(zoo_lc_file("little_4var.circuit"));

    lstats = leaf_stats(lc);
    istats = inode_stats(lc);
    nstats = node_stats(lc);

    @test !(PlainTrueNode in keys(lstats));
    @test !(PlainFalseNode in keys(lstats));
    @test lstats[PlainLiteralNode] == 8;

    @test istats[(Plain⋀Node, 2)] == 9;
    @test istats[(Plain⋁Node, 1)] == 10;
    @test istats[(Plain⋁Node, 4)] == 2;

    for t in keys(lstats)
        @test lstats[t] == nstats[t]
    end
    for t in keys(istats)
        @test istats[t] == nstats[t]
    end

    @test num_nodes(lc) == 29
    @test num_variables(lc) == 4;
    @test num_edges(lc) == 36;

end


text:
# A NamedArray is like a list in R or a DataFrame in Julia without the
# requirement that columns be of equal length. The main reason for its
# existence is to allow creation of a DataFrame from unequal column
# lengths like the following:
#   DataFrame(quote
#       a = 1
#       b = [1:5]
#       c = [1:10]
#   end)
type NamedArray <: Associative{Any,Any}
    data::Vector{Any}
    idx::AbstractIndex
    function NamedArray(data::Vector, idx::AbstractIndex)
        if length(idx) != length(data)
            error("index/names must be the same length as the data")
        end
        new(data, idx)
    end
end
NamedArray() = NamedArray({}, Index())

Base.length(x::NamedArray) = length(x.idx)
Base.names(x::NamedArray) = names(x.idx)

Base.getindex(x::NamedArray, c) = x[x.idx[c]]
Base.getindex(x::NamedArray, c::Integer) = x.data[c]
Base.getindex(x::NamedArray, c::Vector{Int}) = NamedArray(x.data[c], names(x)[c])

function Base.setindex!(x::NamedArray, newdata, ipos::Integer)
    if ipos > 0 && ipos <= length(x)
        x.data[ipos] = newdata
    else
        throw(ArgumentError("Can't replace a non-existent array position"))
    end
    x
end
function Base.setindex!(x::NamedArray, newdata, name)
    ipos = get(x.idx.lookup, name, 0)
    if ipos > 0
        # existing
        setindex!(x, newdata, ipos)
    else
        # new
        push!(x.idx, name)
        push!(x.data, newdata)
    end
    x
end

# Associative methods:
Base.has(x::NamedArray, key) = has(x.idx, key)
Base.get(x::NamedArray, key, default) = has(x, key) ? x[key] : default
Base.keys(x::NamedArray) = keys(x.idx)
Base.values(x::NamedArray) = x.data
# Collection methods:
Base.start(x::NamedArray) = 1
Base.done(x::NamedArray, i) = i > length(x.data)
Base.next(x::NamedArray, i) = ((x.idx.names[i], x[i]), i + 1)
Base.length(x::NamedArray) = length(x.data)
Base.isempty(x::NamedArray) = length(x.data) == 0


text:
#------------------------------------------------------------------------------

#=
    Lucky primes are prime numbers used for modular computation.
    The sequence of lucky primes starts with FIRST_LUCKY_PRIME

    Good primes are prime numbers used for checking correctness.
    The sequence of good primes starts with FIRST_GOOD_PRIME

    There are 50697537 primes between
        FIRST_LUCKY_PRIME and FIRST_GOOD_PRIME
=#

const FIRST_LUCKY_PRIME = 2^31 - 1
const FIRST_GOOD_PRIME  = 2^30 + 3

# const FIRST_LUCKY_PRIME =  5
# const FIRST_GOOD_PRIME  = 16411

mutable struct PrimeTracker
    # integer generator coefficients from the input
    coeffs::Vector{Vector{CoeffZZ}}
    # buffer for BigInts
    buf::BigInt

    # current lucky prime
    luckyprime::UInt64
    # current good prime
    goodprime::UInt64

    # all used lucky primes
    primes::Vector{UInt64}
    # product of all used lucky primes
    modulo::BigInt

    function PrimeTracker(coeffs_zz::Vector{Vector{CoeffZZ}})
        new(coeffs_zz, BigInt(), FIRST_LUCKY_PRIME, FIRST_GOOD_PRIME,
            UInt64[], BigInt(1))
    end
end

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

# check if `prime` is lucky w.r.t. coefficients from tracker
function isluckyprime(tracker::PrimeTracker, prime::UInt64)
    buf = tracker.buf
    p   = BigInt(prime)
    for poly in tracker.coeffs
        for c in poly
            if Base.GMP.MPZ.cmp_ui(Base.GMP.MPZ.tdiv_r!(buf, c, p), 0) == 0
                 return false
            end
        end
    end
    return true
end

function nextluckyprime!(tracker::PrimeTracker)
    prime = tracker.luckyprime

    while !isluckyprime(tracker, prime)
        prime = nextprime(prime + 1)
    end

    tracker.luckyprime = nextprime(prime + 1)
    push!(tracker.primes, prime)

    return prime
end

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

function updatemodulo!(tracker::PrimeTracker)
    Base.GMP.MPZ.mul_ui!(tracker.modulo, last(tracker.primes))
end

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

function nextgoodprime!(tracker::PrimeTracker)
    prime = tracker.goodprime

    while !isluckyprime(tracker, prime)
        prime = nextprime(prime + 1)
    end

    if prime >= FIRST_LUCKY_PRIME
        # @warn "Too large coefficient encountered in the basis. This should not happen normally, please sumbit an issue"
    end

    tracker.goodprime = nextprime(prime + 1)

    return prime
end

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


text:
using PyPlot

function hashf(args)
    return mod1(sum(args), length(args))
end

rands = [[rand(1:10) for i in 1:10] for k in 1:10]
hashed_rands = [hashf(args) for args in rands]

println(hash(rands[0]...))


text:
function eval_ast(interpstate::InterpState, mod, expr)
    interpstate.options.debug && @show :eval_ast mod expr
    try
        # ans = Core.eval(mod, expr)
        ans = @eval mod $expr
        interpstate.options.debug && @show :eval_ast ans
        return ans
    catch exception
        interpstate.options.debug && @show :eval_ast mod expr exception
        rethrow()
    end
end

function eval_lower_ast(interpstate::InterpState, mod, expr)
    interpstate.options.debug && @show :eval_lower_ast mod expr
    try
        lwr = Meta.lower(mod, expr)::Expr
        if lwr.head == :thunk
            codestate = code_state_from_thunk(interpstate, mod, lwr.args[1]::Core.CodeInfo)
            ans = interpret_lower(codestate)
        else
            @show lwr mod expr
            # ans = Core.eval(mod, expr)
            ans = @eval mod $expr
        end
        interpstate.options.debug && @show :eval_lower_ast ans
        return ans
    catch exception
        interpstate.options.debug && @show :eval_lower_ast mod expr exception
        rethrow()
    end
end

function interpret_ast_comparison(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :comparison args
    expr = Expr(:comparison, args...)
    eval_lower_ast(interpstate, mod, expr)
end

function interpret_ast_or(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :(||) args
    expr = Expr(:(||), args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_and(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :(&&) args
    expr = Expr(:(&&), args...)
    eval_lower_ast(interpstate, mod, expr)
end

function interpret_ast_primitive(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :primitive args
    expr = Expr(:primitive, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_call(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :call args
    expr = Expr(:call, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_macrocall(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :macrocall args
    expr = Expr(:macrocall, args...)
    # Macros be better handled by the compiler?
    eval_ast(interpstate, mod, expr)
end
function interpret_ast_broadcast(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :(.=) args
    expr = Expr(:(.=), args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_assign(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :(=) args
    expr = Expr(:(=), args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_const(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :const args
    expr = Expr(:const, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_local(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :local args
    expr = Expr(:local, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_global(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :global args
    expr = Expr(:global, args...)
    eval_ast(interpstate, mod, expr)
end
function interpret_ast_let(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :let args
    expr = Expr(:let, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_block(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :block args
    local linenumbernode
    try
        local ans
        for arg in args
            ans = interpret_ast_node(interpstate, mod, arg)
            if ans isa LineNumberNode
                linenumbernode = ans
            end
        end
        ans
    catch
        @show linenumbernode
        rethrow()
    end
end
function interpret_ast_try(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :try args
    expr = Expr(:try, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_do(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :do args
    expr = Expr(:do, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_if(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :if args
    expr = Expr(:if, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_for(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :for args
    expr = Expr(:for, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_while(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :while args
    expr = Expr(:while, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_function(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :function args
    expr = Expr(:function, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_where(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :function args
    expr = Expr(:where, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_macro(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :macro args
    expr = Expr(:macro, args...)
    # Macros be better handled by the compiler?
    eval_ast(interpstate, mod, expr)
end
function interpret_ast_struct(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :struct args
    expr = Expr(:struct, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_curly(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :curly args
    expr = Expr(:curly, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_abstract(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :abstract args
    expr = Expr(:abstract, args...)
    eval_lower_ast(interpstate, mod, expr)
end
function interpret_ast_using(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :using args
    expr = Expr(:using, args...)
    eval_ast(interpstate, mod, expr)
end
function interpret_ast_import(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :import args
    expr = Expr(:import, args...)
    eval_ast(interpstate, mod, expr)
end
function interpret_ast_export(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :export args
    expr = Expr(:export, args...)
    # Can we really lower exports?
    eval_lower_ast(interpstate, mod, expr)
end

function interpret_ast_error(interpstate::InterpState, mod, args)
    @show mod :error args
end
function interpret_ast_incomplete(interpstate, mod, args)
    @show mod :incomplete args
end
function interpret_ast_toplevel(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :toplevel args
    local linenumbernode
    try
        local ans
        for arg in args
            ans = interpret_ast_node(interpstate, mod, arg)
            if ans isa LineNumberNode
                linenumbernode = ans
            end
        end
        ans
    catch
        @show linenumbernode
        rethrow()
    end
end
function interpret_ast_module(interpstate::InterpState, mod, args)
    interpstate.options.debug && @show mod :module args
    expr = Expr(:module, args...)
    eval_ast(interpstate, mod, expr)
end

function interpret_ast_expr(interpstate::InterpState, mod, expr::Expr)
    if expr.head == :comparison
        interpret_ast_comparison(interpstate, mod, expr.args)
    elseif expr.head == :(||)
        interpret_ast_or(interpstate, mod, expr.args)
    elseif expr.head == :(&&)
        interpret_ast_and(interpstate, mod, expr.args)
    elseif expr.head == :primitive
        interpret_ast_primitive(interpstate, mod, expr.args)
    elseif expr.head == :call
        interpret_ast_call(interpstate, mod, expr.args)
    elseif expr.head == :macrocall
        interpret_ast_macrocall(interpstate, mod, expr.args)
    elseif expr.head == :(.=)
        interpret_ast_broadcast(interpstate, mod, expr.args)
    elseif expr.head == :(=)
        interpret_ast_assign(interpstate, mod, expr.args)
    elseif expr.head == :const
        interpret_ast_const(interpstate, mod, expr.args)
    elseif expr.head == :local
        interpret_ast_local(interpstate, mod, expr.args)
    elseif expr.head == :global
        interpret_ast_global(interpstate, mod, expr.args)
    elseif expr.head == :let
        interpret_ast_let(interpstate, mod, expr.args)
    elseif expr.head == :block
        interpret_ast_block(interpstate, mod, expr.args)
    elseif expr.head == :try
        interpret_ast_try(interpstate, mod, expr.args)
    elseif expr.head == :do
        interpret_ast_do(interpstate, mod, expr.args)
    elseif expr.head == :if
        interpret_ast_if(interpstate, mod, expr.args)
    elseif expr.head == :for
        interpret_ast_for(interpstate, mod, expr.args)
    elseif expr.head == :while
        interpret_ast_while(interpstate, mod, expr.args)
    elseif expr.head == :function
        interpret_ast_function(interpstate, mod, expr.args)
    elseif expr.head == :where
        interpret_ast_where(interpstate, mod, expr.args)
    elseif expr.head == :macro
        interpret_ast_macro(interpstate, mod, expr.args)
    elseif expr.head == :struct
        interpret_ast_struct(interpstate, mod, expr.args)
    elseif expr.head == :curly
        interpret_ast_curly(interpstate, mod, expr.args)
    elseif expr.head == :abstract
        interpret_ast_abstract(interpstate, mod, expr.args)
    elseif expr.head == :using
        interpret_ast_using(interpstate, mod, expr.args)
    elseif expr.head == :import
        interpret_ast_import(interpstate, mod, expr.args)
    elseif expr.head == :export
        interpret_ast_export(interpstate, mod, expr.args)
    elseif expr.head == :error
        interpret_ast_error(interpstate, mod, expr.args)
    elseif expr.head == :incomplete
        interpret_ast_incomplete(interpstate, mod, expr.args)
    elseif expr.head == :toplevel
        interpret_ast_toplevel(interpstate, mod, expr.args)
    elseif expr.head == :module
        interpret_ast_module(interpstate, mod, expr.args)
    else
        @show expr
        @assert false
    end
end

function interpret_ast_node(interpstate::InterpState, mod, node)
    @nospecialize node
    if node isa Expr
        expr = node
        interpret_ast_expr(interpstate, mod, expr)
    elseif node isa Symbol
        symbol = node
        eval_ast(interpstate, mod, symbol)
    elseif node isa LineNumberNode
        linenumbernode = node
        linenumbernode
    else
        @show typeof(node) node
        @assert false
    end
end

function interpret_ast(mod, expr, options)
    interpstate = interp_state(mod, options)
    try
        interpret_ast_node(interpstate, mod, expr)
    finally
        #= for (wt, (meth, sparam_vals, src)) in interpstate.meths
            @show wt meth sparam_vals
        end =#
        if interpstate.options.stats
            for (mod_or_meth, stats) in sort(collect(interpstate.stats), by=kv -> kv[2].elapsed)
                @show mod_or_meth Int(stats.calls) Int(stats.elapsed)
            end
        end
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 8be48d92-2916-40a9-bd99-c3967465e135
using Statistics, Plots, StatsPlots, DataFrames, XLSX, Gadfly, Query

# ╔═╡ 92541b63-eea6-46bc-9961-6c10edf9faea
Gadfly.set_default_plot_size(25cm, 15cm)

# ╔═╡ 50544e4d-5291-46bb-956b-153cb65d527a
"""create stacked histograms showing time results of the different methods"""
function plot_df(df)
	Gadfly.plot(df,
		x = :total_elapsed_time,
		ygroup=:method,
		Geom.subplot_grid(Geom.histogram()),
		Guide.title("Method over time histograms"))
end

# ╔═╡ 68b8dc02-b2aa-42a2-8f47-99a1353f511f
"""turn an array into a head=>tail pair, used to add titles to the XLSX sheet when
making it a DataFrame"""
function col_to_pair(col)
	h, t = Iterators.peel(col)
	h => collect(t)
end

# ╔═╡ 943a8e8f-56dd-4a12-855c-3c18449509a7
begin
	f = XLSX.readxlsx("../df.xlsx")
	sh = f["Sheet1"]
	dfm = DataFrame(map(col_to_pair, eachcol(DataFrame(sh[:]))))
end

# ╔═╡ c7b835d6-abff-42b3-b2d6-6b17ab42a9c0
begin
	g1 = groupby(dfm, [:method])
	md"# Method Group Summary", combine(g1, nrow,
		:lowest_gap=>mean,
		:total_elapsed_time=>mean=>:time_mean)
end

# ╔═╡ ae2ef880-85d6-4362-acd4-812a8722ff32
begin
	groups = groupby(dfm, [:dataset, :case])
	dc_sum = combine(groups, nrow,
		:lowest_gap=>mean,
		:total_elapsed_time=>mean=>:time_mean, :method)
	(md"# Dataset/Case Group Summary", dc_sum)
end

# ╔═╡ c9708e1a-2e9d-43cf-b0f6-768bdf5364a0
plot_df(dfm)

# ╔═╡ Cell order:
# ╠═c7b835d6-abff-42b3-b2d6-6b17ab42a9c0
# ╠═ae2ef880-85d6-4362-acd4-812a8722ff32
# ╠═c9708e1a-2e9d-43cf-b0f6-768bdf5364a0
# ╠═92541b63-eea6-46bc-9961-6c10edf9faea
# ╠═50544e4d-5291-46bb-956b-153cb65d527a
# ╠═943a8e8f-56dd-4a12-855c-3c18449509a7
# ╠═68b8dc02-b2aa-42a2-8f47-99a1353f511f
# ╠═8be48d92-2916-40a9-bd99-c3967465e135

max_stars_repo_name: Dibillilia/JuliaOptHeuristics

max_issues_repo_name: Dibillilia/JuliaOptHeuristics

max_forks_repo_name: Dibillilia/JuliaOptHeuristics


text:
function _columns(table)
    try
        return Tables.columns(table)
    catch
        @warn "encountered error during internal invocation of `Tables.columns`; provided table might not be Tables.jl-compliant"
        rethrow()
    end
end

#####
##### validate tables
#####

"""
    Legolas.extract_schema(table)

Attempt to extract Arrow metadata from `table` via `Arrow.getmetadata(table)`.

If Arrow metadata is present and contains `\"$LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY\" => s`, return [`Legolas.Schema(s)`](@ref).

Otherwise, return `nothing`.
"""
function extract_schema(table)
    metadata = Arrow.getmetadata(table)
    if !isnothing(metadata)
        for (k, v) in metadata
            k == LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY && return Schema(v)
        end
    end
    return nothing
end

"""
    Legolas.validate(table, legolas_schema::Legolas.Schema)

Attempt to determine `s::Tables.Schema` from `table` and return `Legolas.validate(s, legolas_schema)`.

If a `Tables.Schema` cannot be determined, a warning message is logged and `nothing` is returned.
"""
function validate(table, legolas_schema::Schema)
    columns = _columns(table)
    Tables.rowcount(columns) > 0 || return nothing
    tables_schema = Tables.schema(columns)
    if tables_schema isa Tables.Schema
        try
            validate(tables_schema, legolas_schema)
        catch
            @warn "provided table's `Tables.Schema` does not appear to match provided `Legolas.Schema`. Run `[Legolas.Row($legolas_schema, r) for r in Tables.rows(t)]` to try converting the table `t` to a compatible representation."
            rethrow()
        end
    else
        @warn "could not determine `Tables.Schema` from provided table; skipping schema validation"
    end
    return nothing
end

"""
    Legolas.validate(table)

If [`Legolas.extract_schema(table)`](@ref) returns a valid `Legolas.Schema`, return `Legolas.validate(table, Legolas.extract_schema(table))`.

Otherwise, if a `Legolas.Schema` isn't found or is invalid, an `ArgumentError` is thrown.
"""
function validate(table)
    schema = Legolas.extract_schema(table)
    isnothing(schema) && throw(ArgumentError("""
                                             could not extract valid `Legolas.Schema` from provided Arrow table;
                                             is it missing the expected custom metadata and/or the expected
                                             \"$LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY\" field?
                                             """))
    return validate(table, schema)
end

#####
##### read/write tables
#####

"""
    Legolas.read(io_or_path; validate::Bool=true)

Read and return an `Arrow.Table` from `io_or_path`.

If `validate` is `true`, `Legolas.validate` will be called on the table before it is returned.

Note that `io_or_path` may be any type that supports `Base.read(io_or_path)::Vector{UInt8}`.
"""
function read(io_or_path; validate::Bool=true)
    table = read_arrow(io_or_path)
    validate && Legolas.validate(table)
    return table
end

"""
    Legolas.write(io_or_path, table, schema::Schema; validate::Bool=true, kwargs...)

Write `table` to `io_or_path`, inserting the appropriate `$LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY`
field in the written out Arrow metadata.

If `validate` is `true`, `Legolas.validate` will be called on the table before it written out.

Any other provided `kwargs` are forwarded to an internal invocation of `Arrow.write`.

Note that `io_or_path` may be any type that supports `Base.write(io_or_path, bytes::Vector{UInt8})`.
"""
function write(io_or_path, table, schema::Schema; validate::Bool=true,
               metadata=Arrow.getmetadata(table), kwargs...)
    validate && Legolas.validate(table, schema)
    schema_metadata = LEGOLAS_SCHEMA_QUALIFIED_METADATA_KEY => schema_qualified_string(schema)
    if isnothing(metadata)
        metadata = (schema_metadata,)
    else
        metadata = Set(metadata)
        push!(metadata, schema_metadata)
    end
    write_arrow(io_or_path, table; metadata=metadata, kwargs...)
    return table
end

"""
    Legolas.tobuffer(args...; kwargs...)

A convenience function that constructs a fresh `io::IOBuffer`, calls
`Legolas.write(io, args...; kwargs...)`, and returns `seekstart(io)`.

Analogous to the `Arrow.tobuffer` function.
"""
function tobuffer(args...; kwargs...)
    io = IOBuffer()
    Legolas.write(io, args...; kwargs...)
    seekstart(io)
    return io
end

#####
##### read/write Arrow content to generic path types
#####
# It would be better if Arrow.jl supported a generic API for nonstandard path-like types so that
# we can avoid potential intermediate copies here, but its documentation is explicit that it only
# supports `Union{IO,String}`.
#
# TODO: upstream improvements to Arrow.jl to obviate these?

write_full_path(path::AbstractString, bytes) = (mkpath(dirname(path)); Base.write(path, bytes))
write_full_path(path, bytes) = Base.write(path, bytes)

read_arrow(io_or_path::Union{IO,String,Vector{UInt8}}) = Arrow.Table(io_or_path)
read_arrow(path) = read_arrow(Base.read(path))

write_arrow(path::String, table; kwargs...) = Arrow.write(path, table; kwargs...)
write_arrow(io::IO, table; kwargs...) = Arrow.write(io, table; file=get(kwargs, :file, true), kwargs...)
write_arrow(path, table; kwargs...) = (io = IOBuffer(); write_arrow(io, table; kwargs...); write_full_path(path, take!(io)))

#####
##### Tables.jl operations
#####

"""
    locations(collections::Tuple)

Return a `Dict` whose keys are the set of all elements across all provided collections,
and whose values are the indices that locate each corresponding element across all
provided collecitons.

Specifically, `locations(collections)[k][i]` will return a `Vector{Int}` whose elements
are the index locations of `k` in `collections[i]`. If `!(k in collections[i])`, this
`Vector{Int}` will be empty.

For example:

```
julia> Legolas.locations((['a', 'b', 'c', 'f', 'b'],
                          ['d', 'c', 'e', 'b'],
                          ['f', 'a', 'f']))
Dict{Char, Tuple{Vector{Int64}, Vector{Int64}, Vector{Int64}}} with 6 entries:
  'f' => ([4], [], [1, 3])
  'a' => ([1], [], [2])
  'c' => ([3], [2], [])
  'd' => ([], [1], [])
  'e' => ([], [3], [])
  'b' => ([2, 5], [4], [])
```

This function is useful as a building block for higher-level tabular operations
that require indexing/grouping along specific sets of elements.
"""
function locations(collections::T) where {T<:Tuple}
    N = fieldcount(T)
    K = promote_type(eltype.(collections)...)
    results = Dict{K,NTuple{N,Vector{Int}}}()
    for (c, collection) in enumerate(collections)
        for (i, item) in enumerate(collection)
            push!(get!(() -> ntuple(_ -> Int[], N), results, item)[c], i)
        end
    end
    return results
end

function _iterator_for_column(table, c)
    Tables.columnaccess(table) && return Tables.getcolumn(_columns(table), c)
    # there's not really a need to actually materialize this iterable
    # for the caller, but doing so allows the caller to more usefully
    # employ `eltype` on this function's output (since e.g. a generator
    # would just return `Any` for the eltype)
    return [Tables.getcolumn(r, c) for r in Tables.rows(table)]
end

"""
    gather(column_name, tables...; extract=((table, idxs) -> view(table, idxs, :)))

Gather rows from `tables` into a unified cross-table index along `column_name`. Returns
a `Dict` whose keys are the unique values of `column_name` across `tables`, and whose
values are tuples of the form:

    (rows_matching_key_in_table_1, rows_matching_key_in_table_2, ...)

The provided `extract` function is used to extract rows from each table; it takes
as input a table and a `Vector{Int}` of row indices, and returns the corresponding
subtable. The default definition is sufficient for `DataFrames` tables.

Note that this function may internally call `Tables.columns` on each input table, so
it may be slower and/or require more memory if `any(!Tables.columnaccess, tables)`.

Note that we intend to eventually migrate this function from Legolas.jl to a more appropriate package.
"""
function gather(column_name, tables::Vararg{Any,N};
                extract=((cols, idxs) -> view(cols, idxs, :))) where {N}
    iters = ntuple(i -> _iterator_for_column(tables[i], column_name), N)
    return Dict(id => ntuple(i -> extract(tables[i], locs[i]), N) for (id, locs) in locations(iters))
end

"""
    materialize(table)

Return a fully deserialized copy of `table`.

This function is useful when `table` has built-in deserialize-on-access or
conversion-on-access behavior (like `Arrow.Table`) and you'd like to pay
such access costs upfront before repeatedly accessing the table.

Note that we intend to eventually migrate this function from Legolas.jl to a more appropriate package.
"""
materialize(table) = map(collect, Tables.columntable(table))


text:
# This file provides most of the implementation of multi-robot target *coverage*
# problems.
#
# These coverage problems involve controlling robots to cover a number of
# targets e.g. to keep them within range of onboard sensors.
#
# This problem is applicable to cinematography whereas robots may have good
# estimations of position but otherwise seek to "cover" targets with sensors or
# cameras.

export MultiRobotTargetCoverageConfigs, MultiRobotTargetCoverageProblem

# Stores configuration variables for multi-robot target tracking
struct MultiRobotTargetCoverageConfigs
  grid::Grid
  sensor::CoverageSensor
  horizon::Int64

  # Configuration for the MCTS solver
  solver_iterations::Int64

  robot_target_range_limit::Float64

  function MultiRobotTargetCoverageConfigs(;
                                           grid,
                                           horizon=default_horizon,
                                           sensor=CoverageSensor(),
                                           solver_iterations =
                                           default_num_iterations[horizon],
                                           robot_target_range_limit = Inf
                                          )

    new(grid, sensor, horizon, solver_iterations, robot_target_range_limit)
  end
end

# Constructor that can infer parameters based on the number of robots
#
# Note that inferred parameters can still be overriden, but we probably don't
# want to do so
function MultiRobotTargetCoverageConfigs(num_robots;
                                         kwargs...)
  grid=Grid(num_robots=num_robots)
  MultiRobotTargetCoverageConfigs(;grid=grid, kwargs...)
end

# Copy constructor with overrides
function
  MultiRobotTargetCoverageConfigs(configs::MultiRobotTargetCoverageConfigs;
                                  kwargs...)

  fields = fieldnames(MultiRobotTargetCoverageConfigs)
  old_settings = Dict(field => getfield(configs, field) for field in fields)
  settings = merge(old_settings, kwargs)

  MultiRobotTargetCoverageConfigs(; settings...)
end

# Solution elements consist of the robot index and the associated trajectory
# We construct solutions as such because the output may not have the same
# ordering as the robots
#
# Warning: this "problem" object may become invalid if any of the underlying
# objects change and may copy some but not of the inputs.
struct MultiRobotTargetCoverageProblem <: AbstractMultiRobotProblem
  # Target tracking problems are defined by vectors of robot states
  partition_matroid::Vector{State}

  target_states::Vector{State}

  configs::MultiRobotTargetCoverageConfigs
end

# Construct a target coverage problem with configs.
function MultiRobotTargetCoverageProblem(robot_states::Vector{State},
                                         target_states::Vector{State};
                                         kwargs...)
  configs = MultiRobotTargetCoverageConfigs(;kwargs...)

  MultiRobotTargetCoverageProblem(robot_states, target_states, configs)
end

function objective(p::MultiRobotTargetCoverageProblem, X)
  configs = p.configs

  trajectories = map(last, X)

  sum(p.target_states) do target_state
    finite_horizon_coverage(configs.grid, target_state, configs.sensor,
                            trajectories;
                           ).reward
  end
end

function solve_block(p::MultiRobotTargetCoverageProblem, block::Integer,
                     selections::Vector)
  configs = p.configs

  trajectories = map(last, selections)

  problem = SingleRobotTargetCoverageProblem(configs.grid,
                                             configs.sensor,
                                             configs.horizon,
                                             p.target_states,
                                             prior_trajectories=trajectories)

  state = get_state(p, block)

  solution = solve_single_robot(problem, state,
                                n_iterations=configs.solver_iterations)

  (block, solution.trajectory)
end

max_stars_repo_name: mcorah/MultiAgentSensing

max_issues_repo_name: mcorah/MultiAgentSensing

max_forks_repo_name: mcorah/MultiAgentSensing


text:
"""
    energy_balance!_(leaf::LeafModels{I,<:Missing,A,Gs,S},meteo::AbstractAtmosphere,constants = Constants())

Method for when energy balance is missing (do nothing).

# Arguments

- `leaf::LeafModels{I,<:Missing,A,Gs,S}`: a [`LeafModels`](@ref) struct with a missing energy model.
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

"""
function energy_balance!_(leaf::LeafModels{I,<:Missing,A,Gs,S}, meteo::AbstractAtmosphere, constants = Constants()) where {I,A,Gs,S}
    nothing
end

"""
    energy_balance!_(object::Component{I,<:Missing,S},meteo::AbstractAtmosphere,constants = Constants())

Method for when energy balance is missing (do nothing).

# Arguments

- `object::Component{I,<:Missing,S}`: a [`Component`](@ref) struct with a missing energy model.
- `meteo`: meteorology structure, see [`Atmosphere`](@ref)
- `constants = Constants()`: physical constants. See [`Constants`](@ref) for more details

"""
function energy_balance!_(object::Component{I,<:Missing,S}, meteo::AbstractAtmosphere, constants = Constants()) where {I,A,Gs,S}
    nothing
end


text:
################################################################################
# MIT License
#
# Copyright (c) 2021 Hajime Nakagami
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.
################################################################################
using Test, Firebird, TimeZones, DBInterface, Tables, Dates, Decimals

const DEBUG_PRIVATE_KEY = big"0x60975527035CF2AD1989806F0407210BC81EDC04E2762A56AFD529DDDA2D4393"
const DEBUG_SALT = hex2bytes("02E268803000000079A478A700000002D1A6979000000026E1601C000000054F")

@testset "connection" begin
    user = if haskey(ENV, "ISC_USER")
        ENV["ISC_USER"]
    else
        "sysdba"
    end
    password = if haskey(ENV, "ISC_PASSWORD")
        ENV["ISC_PASSWORD"]
    else
        "masterkey"
    end

    conn = DBInterface.connect(Firebird.Connection, "localhost", user, password, "/tmp/julia_test.fdb"; create_new=true, timezone="Asia/Tokyo")
    DBInterface.execute(
        conn, raw"""
            CREATE TABLE foo (
                a INTEGER NOT NULL,
                b VARCHAR(30) NOT NULL UNIQUE,
                c VARCHAR(1024),
                d DECIMAL(16,3) DEFAULT -0.123,
                e DATE DEFAULT '1967-08-11',
                f TIMESTAMP DEFAULT '1967-08-11 23:45:01.1234',
                g TIME DEFAULT '23:45:01.1234',
                h BLOB SUB_TYPE 0,
                i DOUBLE PRECISION DEFAULT 0.0,
                j FLOAT DEFAULT 0.0,
                PRIMARY KEY (a),
                CONSTRAINT CHECK_A CHECK (a <> 0)
            )"""
    )
    DBInterface.close!(conn)

    conn = DBInterface.connect(Firebird.Connection, "localhost", user, password, "/tmp/julia_test.fdb")
    DBInterface.execute(
        conn, raw"insert into foo(a, b, c, h) values (1, 'a', 'b', 'This is a pen')")
    DBInterface.execute(
        conn, raw"insert into foo(a, b, c, h) values (2, 'A', 'B', NULL)")

    expected = (
        A = Int64[1, 2],
        B = String["a", "A"],
        C = Union{Missing, String}["b", "B"],
        D = Union{Missing, Decimal}[Decimal(-0.123), Decimal(-0.123)],
        E = Union{Missing, Dates.Date}[Dates.Date("1967-08-11"), Dates.Date("1967-08-11")],
        F = Union{Missing, Dates.DateTime}[Dates.DateTime("1967-08-11T23:45:01.123"), Dates.DateTime("1967-08-11T23:45:01.123")],
        G = Union{Missing, Dates.Time}[Dates.Time("23:45:01.123"), Dates.Time("23:45:01.123")],
        H = Union{Missing, Vector{UInt8}}[Vector{UInt8}("This is a pen"), missing],
        I = Union{Missing, Float64}[Float64(0.0), Float64(0.0)],
        J = Union{Missing, Float32}[Float32(0.0), Float32(0.0)],
    )

    cursor = DBInterface.execute(conn, "SELECT * FROM foo")
    @test eltype(cursor) == Firebird.Row
    @test Tables.istable(cursor)
    @test Tables.rowaccess(cursor)
    @test Tables.rows(cursor) === cursor
    @test Tables.schema(cursor) == Tables.Schema(propertynames(expected), eltype.(collect(expected)))
    @test Base.IteratorSize(typeof(cursor)) == Base.HasLength()
    @test length(cursor) == 2

    row = first(cursor)
    @test Base.IndexStyle(typeof(row)) == Base.IndexLinear()
    @test length(row) == length(expected)
    @test propertynames(row) == collect(propertynames(expected))
    for (i, prop) in enumerate(propertynames(row))
        @test getproperty(row, prop) == row[prop] == row[i] == expected[prop][1]
    end

    res = DBInterface.execute(conn, raw"SELECT * FROM foo where a=?", (1, )) |> columntable
    @test length(res[1]) == 1

    # as a prepared statement
    stmt = DBInterface.prepare(conn, raw"SELECT * FROM foo")
    cursor = DBInterface.execute(stmt)
    @test eltype(cursor) == Firebird.Row
    @test Tables.istable(cursor)
    @test Tables.rowaccess(cursor)
    @test Tables.rows(cursor) === cursor
    @test Tables.schema(cursor) == Tables.Schema(propertynames(expected), eltype.(collect(expected)))
    @test Base.IteratorSize(typeof(cursor)) == Base.HasLength()
    @test length(cursor) == 2

    row = first(cursor)
    @test Base.IndexStyle(typeof(row)) == Base.IndexLinear()
    @test length(row) == length(expected)
    @test propertynames(row) == collect(propertynames(expected))
    for (i, prop) in enumerate(propertynames(row))
        @test getproperty(row, prop) == row[prop] == row[i] == expected[prop][1]
    end

    res = DBInterface.execute(stmt) |> columntable
    @test length(res) == 10
    @test length(res[1]) == 2
    @test isequal(res, expected)

    stmt = DBInterface.prepare(conn, raw"SELECT * FROM foo where a=?")
    cursor = DBInterface.execute(stmt, (1,))
    row = first(cursor)
    @test length(row) == length(expected)
    res = DBInterface.execute(stmt, (1, )) |> columntable
    @test length(res[1]) == 1

    DBInterface.close!(stmt)
    DBInterface.close!(conn)
    @test !isopen(conn)
end

@testset "timezone" begin
    user = if haskey(ENV, "ISC_USER")
        ENV["ISC_USER"]
    else
        "sysdba"
    end
    password = if haskey(ENV, "ISC_PASSWORD")
        ENV["ISC_PASSWORD"]
    else
        "masterkey"
    end

    conn = DBInterface.connect(Firebird.Connection, "localhost", user, password, "/tmp/julia_test.fdb"; create_new=true, timezone="Asia/Tokyo")

    cursor = DBInterface.execute(conn, raw"SELECT rdb$get_context('SYSTEM', 'ENGINE_VERSION') version from rdb$database")
    row = first(cursor)
    major_version = parse(Int, split(row.VERSION, ".")[1])
    if major_version < 4
        DBInterface.close!(conn)
        return
    end

    DBInterface.execute(
        conn, raw"""
            CREATE TABLE tz_test (
                id INTEGER NOT NULL,
                t TIME WITH TIME ZONE DEFAULT '12:34:56' NOT NULL,
                ts TIMESTAMP WITH TIME ZONE DEFAULT '1967-08-11 23:45:01' NOT NULL,
                PRIMARY KEY (id)
            )"""
    )

    DBInterface.execute(
        conn, raw"insert into tz_test (id) values (1)")
    DBInterface.execute(
        conn, raw"insert into tz_test (id, t, ts) values (2, '12:34:56 Asia/Seoul', '1967-08-11 23:45:01.0000 Asia/Seoul')")
    DBInterface.execute(
        conn, raw"insert into tz_test (id, t, ts) values (3, '03:34:56 UTC', '1967-08-11 14:45:01.0000 UTC')")

    expected = (
        ID = Int64[1, 2, 3],
        T = TimeZones.ZonedDateTime[TimeZones.ZonedDateTime(0, 1, 1, 12, 53, 55, tz"Asia/Tokyo"), TimeZones.ZonedDateTime(0, 1, 1, 12, 2, 48, tz"Asia/Seoul"), TimeZones.ZonedDateTime(0, 1, 1, 3, 34, 56, tz"UTC")],
        TS = TimeZones.ZonedDateTime[TimeZones.ZonedDateTime(1967, 8, 11, 23, 45, 1, tz"Asia/Tokyo"), TimeZones.ZonedDateTime(1967, 8, 11, 23, 45, 1, tz"Asia/Seoul"), TimeZones.ZonedDateTime(1967, 8, 11, 14, 45, 1, tz"UTC")]
    )

    cursor = DBInterface.execute(conn, "select * from tz_test")
    @test eltype(cursor) == Firebird.Row
    @test Tables.istable(cursor)
    @test Tables.rowaccess(cursor)
    @test Tables.rows(cursor) === cursor
    @test Tables.schema(cursor) == Tables.Schema(propertynames(expected), eltype.(collect(expected)))
    @test Base.IteratorSize(typeof(cursor)) == Base.HasLength()
    @test length(cursor) == 3

    DBInterface.close!(conn)
end

@testset "srp" begin
    user = "SYSDBA"
    password = "masterkey"

    keyA, keya = Firebird.get_client_seed()
    salt = Firebird.get_salt()
    v = Firebird.get_verifier(user, password, salt)
    keyB, keyb = Firebird.get_server_seed(v)
    server_key = Firebird.get_server_session(user, password, salt, keyA, keyB, keyb)
    _, client_key = Firebird.get_client_proof(user, password, salt, keyA, keyB, keya, "Srp")
    @test server_key == client_key

    _, client_key = Firebird.get_client_proof(user, password, salt, keyA, keyB, keya, "Srp256")
    @test server_key == client_key

    keyA, keya = Firebird.get_client_seed(DEBUG_PRIVATE_KEY)
    salt = DEBUG_SALT
    v = Firebird.get_verifier(user, password, salt)
    keyB, keyb = Firebird.get_server_seed(v, DEBUG_PRIVATE_KEY)
    server_key = Firebird.get_server_session(user, password, salt, keyA, keyB, keyb)
    keyM, client_key = Firebird.get_client_proof(user, password, salt, keyA, keyB, keya, "Srp")
    @test server_key == client_key
    @test keyM == hex2bytes("8c12324bb6e9e683a3ee62e13905b95d69f028a9")

    keyM, client_key = Firebird.get_client_proof(user, password, salt, keyA, keyB, keya, "Srp256")
    @test server_key == client_key
    @test keyM == hex2bytes("4675c18056c04b00cc2b991662324c22c6f08bb90beb3677416b03469a770308")
end

@testset "arc4" begin
    session_key = Vector{UInt8}("a key")
    src = Vector{UInt8}("plain text")

    a1 = Firebird.Arc4(session_key)
    enc = Firebird.translate(a1, src)
    @test enc == UInt8[0x4b, 0x4b, 0xdc, 0x65, 0x02, 0xb3, 0x08, 0x17, 0x48, 0x82]
    a2 = Firebird.Arc4(session_key)
    plain = Firebird.translate(a2, enc)
    @test plain == src
end

@testset "chacha20" begin
    session_key = hex2bytes("23AD52B15FA7EBDC4672D72289253D95DC9A4324FC369F593FDCC7733AD77617")
    nonce = hex2bytes("000000005A5F6C13C1F12653")
    enc = hex2bytes("6bd00ba222523f58de196fb471eea08d9fff95b5bbe6123dd3a8b9026ac0fa84")
    chacha = Firebird.ChaCha20(session_key, nonce, Int32(0))
    @test Firebird.translate(chacha, enc) == b"TMCTF{Whose_garden_is_internet?}"

    src = Vector{UInt8}("plain text")
    chacha1 = Firebird.ChaCha20(session_key, nonce, Int32(123))
    enc = Firebird.translate(chacha1, src)
    @test enc == UInt8[0x39, 0xdf, 0x7f, 0xdf, 0xcd, 0xd6, 0x6c, 0x56, 0xe7, 0x62]
    chacha2 = Firebird.ChaCha20(session_key, nonce, Int32(123))
    plain = Firebird.translate(chacha2, enc)
    @test plain == src
end


text:
using PyCall

function __init__()
    py"""
    import pyvista as pv
    import numpy as np

    def show_surf(surf_mat):
        mesh = pv.UniformGrid()
        if surf_mat.shape[1] == 3:
            # 数据为点的三个坐标分量
            mesh.origin = (surf_mat[0, 0], surf_mat[0, 1], surf_mat[0, 2])
            n_x = len(np.unique(surf_mat[:, 0]))
            n_y = len(np.unique(surf_mat[:, 1]))
            Z = surf_mat[:, 2].reshape(n_y, n_x)
            if surf_mat[0, 0] == surf_mat[0, 1]:
                # X = surf_mat[:, 0].reshape(n_y, n_x)
                delta_x = surf_mat[0+n_y, 0]-surf_mat[0, 0]
                # Y = surf_mat[:, 1].reshape(n_y, n_x)
                delta_y = surf_mat[1, 1]-surf_mat[0, 1]
            mesh.spacing = (delta_x, delta_y, 0)
        else:
            # 仅包含Z方向高程值, 则默认间距为1, 初始点为坐标原点(0, 0, 0)
            mesh.origin = (0, 0, 0)
            mesh.spacing = (1, 1, 0)
            n_x, n_y = surf_mat.shape
            Z = surf_mat
        mesh.dimensions = np.array([n_y, n_x, 1])
        mesh.point_arrays["values"] = Z.flatten()
        mesh = mesh.warp_by_scalar()
        plotter = pv.Plotter()
        plotter.add_mesh(mesh, show_edges=True, opacity=0.7)
        plotter.show()
    """
end # function


text:
using Test
using Escape

using DataFrames
using Distributions
using FASTX
using HypothesisTests
using Graphs
using Loo
using MetaGraphs
using Plots
using Serialization
using StanInterface
using StaticArrays
using Statistics
using Suppressor
using LinearAlgebra

include("alleles.jl")
include("abstract_hla_model.jl")
include("hla_data.jl")
include("split_hla_data.jl")
include("hla_dataset.jl")
include("dataset_rousseau.jl")
include("replacement.jl")
include("phylogenetic_tree.jl")
include("epitope_prediction.jl")
include("hla_model_result.jl")
include("loo.jl")
include("epitope_map.jl")
include("hla_model.jl")
include("fisher_test.jl")
include("posterior_predictive_checks.jl")


text:
# Parts of this code were taken / derived from Graphs.jl. See LICENSE for
# licensing details.

#################################################
#
#  Minimum Cut
#
#################################################

"""
    mincut(g, distmx=weights(g))

Return a tuple `(parity, bestcut)`, where `parity` is a vector of integer
values that determines the partition in `g` (1 or 2) and `bestcut` is the
weight of the cut that makes this partition. An optional `distmx` matrix may
be specified; if omitted, edge distances are assumed to be 1.
"""
function mincut(g::AbstractGraph,
    distmx::AbstractMatrix{T}=weights(g)) where T <: Real

    U = eltype(g)
    colormap = zeros(UInt8, nv(g))   ## 0 if unseen, 1 if processing and 2 if seen and closed
    parities = falses(nv(g))
    bestweight = typemax(T)
    cutweight = zero(T)
    visited = zero(U)               ## number of vertices visited
    pq = PriorityQueue{U,T}(Base.Order.Reverse)

    # Set number of visited neighbors for all vertices to 0
    for v in vertices(g)
        pq[v] = zero(T)
    end

    # make sure we have at least two vertices, otherwise, there's nothing to cut,
    # in which case we'll return immediately.
    (haskey(pq, one(U)) && nv(g) > one(U)) || return (Vector{Int8}([1]), cutweight)

    #Give the starting vertex high priority
    pq[one(U)] = one(T)

    while !isempty(pq)
        u = dequeue!(pq)
        colormap[u] = 1

        for v in outneighbors(g, u)
            # if the target of e is already marked then decrease cutweight
            # otherwise, increase it
            ew = distmx[u, v]
            if colormap[v] != 0
                cutweight -= ew
            else
                cutweight += ew
            end
            if haskey(pq, v)
                pq[v] += distmx[u, v]
            end
        end

        colormap[u] = 2
        visited += one(U)
        if cutweight < bestweight && visited < nv(g)
            bestweight = cutweight
            for u in vertices(g)
                parities[u] = (colormap[u] == 2)
            end
        end
    end
    return(convert(Vector{Int8}, parities) .+ one(Int8), bestweight)
end

"""
    maximum_adjacency_visit(g[, distmx][, log][, io])

Return the vertices in `g` traversed by maximum adjacency search. An optional
`distmx` matrix may be specified; if omitted, edge distances are assumed to
be 1. If `log` (default `false`) is `true`, visitor events will be printed to
`io`, which defaults to `STDOUT`; otherwise, no event information will be
displayed.
"""
function maximum_adjacency_visit(g::AbstractGraph,
    distmx::AbstractMatrix{T},
    log::Bool=false,
    io::IO=stdout) where T <: Real

    U = eltype(g)
    pq = PriorityQueue{U,T}(Base.Order.Reverse)
    vertices_order = Vector{U}()
    has_key = ones(Bool, nv(g))
    sizehint!(vertices_order, nv(g))
    # if the graph only has one vertex, we return the vertex by itself.
    nv(g) > one(U) || return collect(vertices(g))

    # Setting intial count to 0
    for v in vertices(g)
        pq[v] = zero(T)
    end

    #Give vertex `1` maximum priority
    pq[one(U)] = one(T)

    #start traversing the graph
    while !isempty(pq)
        u = dequeue!(pq)
        has_key[u] = false
        push!(vertices_order, u)
        log && println(io, "discover vertex: $u")
        for v in outneighbors(g, u)
            log && println(io, " -- examine neighbor from $u to $v")
            if has_key[v]
                ed = distmx[u, v]
                pq[v] += ed
            end
        end
        log && println(io, "close vertex: $u")
    end
    return vertices_order
end

maximum_adjacency_visit(g::AbstractGraph) = maximum_adjacency_visit(g,
    weights(g),
    false,
    stdout)


text:
import .Plots: @recipe, @series, RecipesBase

@recipe function plot(x::AOI{T}) where T
  tl = (x.topleft.lon, x.topleft.lat)
  br = (x.bottomright.lon, x.bottomright.lat)
  mapimg = x.mapimg
  size --> (1086,610)
  yflip --> false
  xlims --> (tl[1], br[1])
  ylims --> (br[2], tl[2])
  framestyle --> :none
  ticks --> nothing
  @series begin
    seriestype := :image
    range(tl[1], br[1]; length=size(mapimg,1)),
    range(tl[2], br[2]; length=size(mapimg,2)),
    mapimg[end:-1:1,:]
  end
end

@recipe function plot(x::AbstractVector{<:LLA})
  @series begin
    [v.lon for v ∈ x], [v.lat for v ∈ x]
  end
end

@recipe function plot(x::LLA)
  @series begin
    [x.lon], [x.lat]
  end
end


text:
using KhepriLibrary
using Test

@testset "KhepriLibrary.jl" begin
    # Write your tests here.
end


text:
using StaticArrays
using CLIMA.VariableTemplates
import CLIMA.DGmethods: BalanceLaw,
                        vars_aux, vars_state, vars_gradient, vars_diffusive,
                        flux_nondiffusive!, flux_diffusive!, source!,
                        gradvariables!, diffusive!,
                        init_aux!, init_state!,
                        boundary_state!, wavespeed, LocalGeometry,
                        vars_gradient_laplacian, vars_hyperdiffusive,
                        hyperdiffusive!
using CLIMA.DGmethods.NumericalFluxes: NumericalFluxNonDiffusive,
                                       NumericalFluxDiffusive

abstract type HyperDiffusionProblem end
struct HyperDiffusion{dim, P} <: BalanceLaw
  problem::P
  function HyperDiffusion{dim}(problem::P) where {dim, P <: HyperDiffusionProblem}
    new{dim, P}(problem)
  end
end

vars_aux(::HyperDiffusion, FT) = @vars(D::SMatrix{3, 3, FT, 9})
#
# Density is only state
vars_state(::HyperDiffusion, FT) = @vars(ρ::FT)

# Take the gradient of density
vars_gradient(::HyperDiffusion, FT) = @vars(ρ::FT)
# Take the gradient of laplacian of density
vars_gradient_laplacian(::HyperDiffusion, FT) = @vars(ρ::FT)

vars_diffusive(::HyperDiffusion, FT) = @vars()
# The hyperdiffusion DG auxiliary variable: D ∇ Δρ
vars_hyperdiffusive(::HyperDiffusion, FT) = @vars(σ::SVector{3,FT})

function flux_nondiffusive!(m::HyperDiffusion, flux::Grad, state::Vars,
                            aux::Vars, t::Real)
end

"""
    flux_diffusive!(m::HyperDiffusion, flux::Grad, state::Vars,
                     auxDG::Vars, auxHDG::Vars, aux::Vars, t::Real)

Computes diffusive flux `F` in:

```
∂ρ
-- = - ∇ • (σ) = - ∇ • F
∂t
```
Where

 - `σ` is hyperdiffusion DG auxiliary variable (`σ = D ∇ Δρ` with D being the hyperdiffusion tensor)
"""
function flux_diffusive!(m::HyperDiffusion, flux::Grad, state::Vars,
                         auxDG::Vars, auxHDG::Vars, aux::Vars, t::Real)
  σ = auxHDG.σ
  flux.ρ += σ
end

"""
    gradvariables!(m::HyperDiffusion, transform::Vars, state::Vars,
                   aux::Vars, t::Real)

Set the variable to take the gradient of (`ρ` in this case)
"""
function gradvariables!(m::HyperDiffusion, transform::Vars, state::Vars,
                        aux::Vars, t::Real)
  transform.ρ = state.ρ
end

diffusive!(m::HyperDiffusion, _...) = nothing
function hyperdiffusive!(m::HyperDiffusion, auxHDG::Vars, gradvars::Grad,
                         state::Vars, aux::Vars, t::Real)
  ∇Δρ = gradvars.ρ
  D = aux.D
  auxHDG.σ = D * ∇Δρ
end

"""
    source!(m::HyperDiffusion, _...)

There is no source in the hyperdiffusion model
"""
source!(m::HyperDiffusion, _...) = nothing

"""
    init_aux!(m::HyperDiffusion, aux::Vars, geom::LocalGeometry)

initialize the auxiliary state
"""
function init_aux!(m::HyperDiffusion, aux::Vars, geom::LocalGeometry)
  init_hyperdiffusion_tensor!(m.problem, aux, geom)
end

function init_state!(m::HyperDiffusion, state::Vars, aux::Vars,
                     coords, t::Real)
  initial_condition!(m.problem, state, aux, coords, t)
end

boundary_state!(nf, ::HyperDiffusion, _...) = nothing

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
using Test

using MathOptInterface
const MOI = MathOptInterface

include("dummy.jl")

"""
    MutLessThan{T<:Real} <: MOI.AbstractScalarSet

A mutable `LessThan`-like set to test `copy` of indicator set
"""
mutable struct MutLessThan{T<:Real} <: MOI.AbstractScalarSet
    upper::T
    MutLessThan(v::T) where {T<:Real} = new{T}(v)
end

Base.copy(mlt::MutLessThan) = MutLessThan(Base.copy(mlt.upper))

@testset "Sets" begin
    @testset "==" begin
        # By default, `==` redirects to `===`, it works for bits type
        # but not for `BigInt`s. We define functions creating different
        # instances so that `a() !== a()`.
        a() = big(1)
        b() = big(2)
        @test a() !== a()
        @test b() !== b()
        for S in [MOI.LessThan, MOI.GreaterThan, MOI.EqualTo, MOI.PowerCone, MOI.DualPowerCone]
            @test S(a()) == S(a())
            @test S(a()) != S(b())
            @test S(b()) == S(b())
            @test S(b()) != S(a())
        end
        for S in [MOI.Interval, MOI.Semicontinuous, MOI.Semiinteger]
            @test S(a(), b()) == S(a(), b())
            @test S(a(), b()) != S(b(), a())
            @test S(a(), b()) != S(b(), b())
            @test S(a(), b()) != S(a(), a())
            @test S(a(), a()) != S(b(), b())
            @test S(a(), a()) == S(a(), a())
        end
        S = MOI.IndicatorSet
        A() = MOI.LessThan(a())
        B() = MOI.LessThan(b())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) == S{MOI.ACTIVATE_ON_ZERO}(A())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ONE}(A())
        @test S{MOI.ACTIVATE_ON_ZERO}(A()) != S{MOI.ACTIVATE_ON_ZERO}(B())
        @test S{MOI.ACTIVATE_ON_ONE}(A()) != S{MOI.ACTIVATE_ON_ONE}(B())
    end
    @testset "Copy" begin
        @testset "for $S" for S in [MOI.SOS1, MOI.SOS2]
            s = S([1.0])
            s_copy = copy(s)
            s_copy.weights[1] = 2.0
            @test s.weights[1] == 1.0
        end
        @testset "IndicatorSet" begin
            s1 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(4.0))
            s2 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MOI.GreaterThan(4.0))
            s1_copy = copy(s1)
            s2_copy = copy(s2)
            @test s1_copy isa MOI.IndicatorSet{MOI.ACTIVATE_ON_ONE}
            @test s1 == s1_copy
            @test s2_copy isa MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}
            @test s2 == s2_copy
            s3 = MOI.IndicatorSet{MOI.ACTIVATE_ON_ZERO}(MutLessThan(4.0))
            s3_copy = copy(s3)
            @test s3.set.upper ≈ 4.0
            s3_copy.set.upper = 5.0
            @test s3.set.upper ≈ 4.0
            @test s3_copy.set.upper ≈ 5.0
        end
    end
    @testset "Broadcast" begin
        model = DummyModelWithAdd()
        x = MOI.add_variables(model, 3)
        cis = MOI.add_constraint.(model, x, MOI.EqualTo(0.0))
        @test cis isa Vector{MOI.ConstraintIndex{MOI.SingleVariable,
                                                 MOI.EqualTo{Float64}}}
        @test length(cis) == 3
    end
end


text:
printstyled(stdout, "  mathematical properties\n", color=:light_green)
H = randSeisHdr()
V = randSeisEvent()
S = randSeisData()
T = convert(EventTraceData, S)
C = randSeisChannel()
D = convert(EventChannel, C)

printstyled(stdout, "    reflexivity\n", color=:light_green)
@test C==C
@test D==D
@test H==H
@test S==S
@test T==T
@test V==V

printstyled(stdout, "    commutativity\n", color=:light_green)

printstyled(stdout, "      S1 + S2 == S2 + S1\n", color=:light_green)

# SeisData + SeisData
S1 = randSeisData()
S2 = randSeisData()
@test S1 + S2 == S2 + S1

# EventTraceData + EventTraceData
S1 = convert(EventTraceData, randSeisData())
S2 = convert(EventTraceData, randSeisData())
@test S1 + S2 == S2 + S1

printstyled(stdout, "      S + C == C + S\n", color=:light_green)

# SeisData + SeisChannel
S = randSeisData()
C = randSeisChannel()
U = deepcopy(S)
push!(U, C)
U = sort(U)
@test C + S == S + C == U
@test findid(C, U) == findid(C, S + C) == findid(C, C + S)

# EventTraceData + EventChannel
S = convert(EventTraceData, randSeisData())
C = convert(EventChannel, randSeisChannel())
U = deepcopy(S)
push!(U, C)
U = sort(U)
@test C + S == S + C == U
@test findid(C, U) == findid(C, S + C) == findid(C, C + S)

printstyled(stdout, "      C1 + C2 == C2 + C1\n", color=:light_green)

# SeisChannel + SeisChannel
C1 = randSeisChannel()
C2 = randSeisChannel()
@test C1 + C2 == C2 + C1

# EventChannel + EventChannel
C1 = convert(EventChannel,randSeisChannel())
C2 = convert(EventChannel,randSeisChannel())
@test C1 + C2 == C2 + C1

printstyled(stdout, "      S + U - U == S (for sorted S)\n", color=:light_green)
# SeisData + SeisData - SeisData
(S, T) = mktestseis()
S = sort(S)
@test (S + T - T) == S

# EventTraceData + EventTraceData - EventTraceData
(S, T) = mktestseis()
S = convert(EventTraceData, S)
T = convert(EventTraceData, T)
@test (S + T - T) == sort(S)

printstyled(stdout, "    associativity\n", color=:light_green)

printstyled(stdout, "      (S1 + S2) + S3 == S1 + (S2 + S3)\n", color=:light_green)
# SeisData + SeisData + SeisData
S1 = randSeisData()
S2 = randSeisData()
S3 = randSeisData()
@test S1 + (S2 + S3) == (S1 + S2) + S3

printstyled(stdout, "      (S1 + S2) + C == S1 + (S2 + C)\n", color=:light_green)
# SeisData + SeisChannel
C1 = randSeisChannel()
@test S1 + (S2 + C1) == (S1 + S2) + C1

# EventTraceData + EventTraceData + EventTraceData
S1 = convert(EventTraceData, randSeisData())
S2 = convert(EventTraceData, randSeisData())
S3 = convert(EventTraceData, randSeisData())
@test S1 + (S2 + S3) == (S1 + S2) + S3

# EventTraceData + EventChannel
C1 = convert(EventChannel, randSeisChannel())
@test S1 + (S2 + C1) == (S1 + S2) + C1