text:
using POMDPs
using FiniteHorizonPOMDPs
using POMDPModelTools

struct FHExampleState
    position::Int64
    done::Bool # are we in a terminal state?
end

# initial state constructor
FHExampleState(position::Int64)::FHExampleState = FHExampleState(position, false)

struct FHExample <: MDP{FHExampleState, Symbol} # Note that our MDP is parametarized by the state and the action
    no_states::Int64 # number od states
    actions::Vector{Symbol}
    actionCost::Float64
    actionsImpact::Base.ImmutableDict{Symbol, Int64}
    reward_states::Vector{Int64}
    reward::Float64
    discount_factor::Float64 # discount factor
    noise::Float64
end

function FHExample(;no_states=11, actionCost=1., reward=-10., discount_factor=1., noise=.3)
    actions = [:l, :r]
    actionsImpact = Base.ImmutableDict(:l => -1, :r => 1)
    reward_states = [1, no_states]
    return FHExample(no_states, actions, actionCost, actionsImpact, reward_states, reward, discount_factor, noise)
end

POMDPs.isterminal(mdp::FHExample, s::FHExampleState)::Bool = s.position in mdp.reward_states
POMDPs.isterminal(mdp::FHExample, position::Int64)::Bool = position in mdp.reward_states

POMDPs.initialstate(mdp::FHExample) = SparseCat([FHExampleState(max(1, ceil(Int64, mdp.no_states / 2)))], [1.])

function POMDPs.reward(mdp::FHExample, s::FHExampleState, a::Symbol, sp::FHExampleState)::Float64
    isterminal(mdp, sp.position) ? mdp.reward : mdp.actionCost
end

function POMDPs.actions(mdp::FHExample, s::FHExampleState)
    mdp.actions
end

function POMDPs.actions(mdp::FHExample)
    mdp.actions
end

# Creates mdp.no_states states
function POMDPs.states(mdp::FHExample)::Array{FHExampleState}
    mdp_states = FHExampleState[]
    for i=1:mdp.no_states
        push!(mdp_states, FHExampleState(i, isterminal(mdp, i)))
    end

    return mdp_states
end

POMDPs.stateindex(mdp::FHExample, s::FHExampleState)::Int64 = s.position
POMDPs.actions(mdp::FHExample)::Vector{Symbol} = mdp.actions
POMDPs.actionindex(mdp::FHExample, a::Symbol)::Int64 = findfirst(x->x==a, POMDPs.actions(mdp))
POMDPs.discount(mdp::FHExample)::Number = mdp.discount_factor

# returns transition distributions - works only for 1D Gridworld with possible moves to left and to right
function POMDPs.transition(mdp::FHExample, s::FHExampleState, a::Symbol)::SparseCat{Vector{FHExampleState},Vector{Float64}}
    sp = FHExampleState[]
    prob = Float64[]

    # add original transition target and probability
    position = s.position + mdp.actionsImpact[a]
    push!(sp, FHExampleState(position, isterminal(mdp, position)))
    push!(prob, 1. - mdp.noise)

    # add noise transition target and probability
    noise_action = a == :l ? :r : :l
    position = s.position + mdp.actionsImpact[noise_action]
    push!(sp, FHExampleState(position, isterminal(mdp, position)))
    push!(prob, mdp.noise)

    return SparseCat(sp, prob)
end


text:
"Policy solution where values and Q-values are directly stored."
@kwdef struct PolicyValue <: PolicySolution
    V::Dict{UInt64,Float64} = Dict()
    Q::Dict{UInt64,Dict{Term,Float64}} = Dict()
end

get_action(sol::PolicyValue, state::State) =
    best_action(sol, state)
best_action(sol::PolicyValue, state::State) =
    argmax(sol.Q[hash(state)])
rand_action(sol::PolicyValue, state::State) =
    best_action(sol, state)
get_value(sol::PolicyValue, state::State) =
    sol.V[hash(state)]
get_value(sol::PolicyValue, state::State, action::Term) =
    sol.Q[hash(state)][action]
get_action_values(sol::PolicyValue, state::State) =
    pairs(sol.Q[hash(state)])

ext: jl


text:
function load_json_data(file_path::String)
    println("This is not yet implemented")
end

function save_json_data(c::Configuration, file_path::String)
    println("This is not yet implemented")
end


text:
using KeywordConstructor
using Base.Test

immutable Foo
    a :: Float64
    b
end

@keyword_constructor Foo

obj = Foo(b=4, a=2)
@test obj.a == 2
@test obj.b == 4


text:
@doc raw"""
    MultivariateBasicStructural(y::Matrix{Fl}, s::Int) where Fl

An implementation of a non-homogeneous seemingly unrelated time series equations for basic structural
state-space model consists of trend (local linear trend) and seasonal components.
It is defined by:
```math
\begin{gather*}
    \begin{aligned}
    y_{t} &=  \mu_{t} + \gamma_{t} + \varepsilon_{t} \quad &\varepsilon_{t} \sim \mathcal{N}(0, \Sigma_{\varepsilon})\\
    \mu_{t+1} &= \mu_{t} + \nu_{t} + \xi_{t} \quad &\xi_{t} \sim \mathcal{N}(0, \Sigma_{\xi})\\
    \nu_{t+1} &= \nu_{t} + \zeta_{t} \quad &\zeta_{t} \sim \mathcal{N}(0, \Sigma_{\zeta})\\
    \gamma_{t+1} &= -\sum_{j=1}^{s-1} \gamma_{t+1-j} + \omega_{t} \quad & \omega_{t} \sim \mathcal{N}(0, \Sigma_{\omega})\\
    \end{aligned}
\end{gather*}
```

# Example
```jldoctest
julia> model = MultivariateBasicStructural(rand(100, 2), 12)
MultivariateBasicStructural model
```

# References
 * Durbin, James, & Siem Jan Koopman. (2012). "Time Series Analysis by State Space Methods: Second Edition." Oxford University Press.
"""
mutable struct MultivariateBasicStructural <: StateSpaceModel
    hyperparameters::HyperParameters
    system::LinearMultivariateTimeInvariant
    seasonality::Int
    results::Results

    function MultivariateBasicStructural(y::Matrix{Fl}, s::Int) where Fl
        p = size(y, 2)
        Z = kron(Matrix{Fl}(I, p, p), [1 0 1 zeros(Fl, 1, s - 2)])
        T = kron(
            Matrix{Fl}(I, p, p),
            [
                1 1 zeros(Fl, 1, s - 1)
                0 1 zeros(Fl, 1, s - 1)
                0 0 -ones(Fl, 1, s - 1)
                zeros(Fl, s - 2, 2) Matrix{Fl}(I, s - 2, s - 2) zeros(Fl, s - 2)
            ],
        )

        R = kron(
            Matrix{Fl}(I, p, p),
            [
                Matrix{Fl}(I, 3, 3)
                zeros(Fl, s - 2, 3)
            ],
        )

        d = zeros(Fl, p)
        c = zeros(Fl, p * (s + 1))
        H = kron(one(Fl), Matrix{Fl}(I, p, p))
        Q = kron(ones(Fl, p, p), Matrix{Fl}(I, 3, 3))

        system = LinearMultivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q)

        names = handle_multivariate_basicstructural_names(p)
        hyperparameters = HyperParameters{Fl}(names)
        return new(hyperparameters, system, s, Results{Fl}())
    end
end

function handle_multivariate_basicstructural_names(p::Int)
    # generate \varepsilon  names
    greek_letters_for_states = ["ε", "ξ", "ζ", "ω"]
    names = String[]
    # Walk the lower triangle positions
    for letter in greek_letters_for_states, i in 1:p, j in 1:i
        if i == j
            push!(names, "sigma2_$(letter)$(i)")
        else
            push!(names, "sigma_$(letter)$(i)sigma_$(letter)$(j)")
        end
    end
    return names
end

function default_filter(model::MultivariateBasicStructural)
    Fl = typeof_model_elements(model)
    steadystate_tol = Fl(1e-5)
    a1 = zeros(Fl, num_states(model))
    P1 = Fl(1e6) .* Matrix{Fl}(I, num_states(model), num_states(model))
    return MultivariateKalmanFilter(
        size(model.system.y, 2), a1, P1, num_states(model), steadystate_tol
    )
end

function initial_hyperparameters!(model::MultivariateBasicStructural)
    Fl = typeof_model_elements(model)
    initial_hyperparameters = Dict{String,Fl}()
    for variable in get_names(model)
        if occursin("sigma2", variable)
            initial_hyperparameters[variable] = one(Fl)
        else
            initial_hyperparameters[variable] = zero(Fl)
        end
    end
    set_initial_hyperparameters!(model, initial_hyperparameters)
    return nothing
end

function constrain_hyperparameters!(model::MultivariateBasicStructural)
    Fl = typeof_model_elements(model)
    p = size(model.system.y, 2)
    # H
    H_unconstrained = zeros(Fl, p, p)
    current_name = 1
    for i in 1:p, j in 1:i
        H_unconstrained[i, j] = get_unconstrained_value(
            model, get_names(model)[current_name]
        )
        current_name += 1
    end
    H_constrained = H_unconstrained * H_unconstrained'
    # Q
    Q_unconstrained = zeros(Fl, 3p, 3p)
    sigma2_ξ = filter(x -> occursin("sigma2_ξ", x), get_names(model))
    sigma2_ξ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 1
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma2_ξ[sigma2_ξ_counter]
            )
            sigma2_ξ_counter += 1
        end
    end
    sigma2_ζ = filter(x -> occursin("sigma2_ζ", x), get_names(model))
    sigma2_ζ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 2
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma2_ζ[sigma2_ζ_counter]
            )
            sigma2_ζ_counter += 1
        end
    end
    sigma2_ω = filter(x -> occursin("sigma2_ω", x), get_names(model))
    sigma2_ω_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 3
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma2_ω[sigma2_ω_counter]
            )
            sigma2_ω_counter += 1
        end
    end
    sigma_ξxsigma_ξx = filter(x -> occursin("sigma_ξ", x), get_names(model))
    sigma_ξxsigma_ξx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 1
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma_ξxsigma_ξx[sigma_ξxsigma_ξx_counter]
            )
            sigma_ξxsigma_ξx_counter += 1
        end
    end
    sigma_ζxsigma_ζx = filter(x -> occursin("sigma_ζ", x), get_names(model))
    sigma_ζxsigma_ζx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 2
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma_ζxsigma_ζx[sigma_ζxsigma_ζx_counter]
            )
            sigma_ζxsigma_ζx_counter += 1
        end
    end
    sigma_ωxsigma_ωx = filter(x -> occursin("sigma_ω", x), get_names(model))
    sigma_ωxsigma_ωx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 3
            Q_unconstrained[pos] = get_unconstrained_value(
                model, sigma_ωxsigma_ωx[sigma_ωxsigma_ωx_counter]
            )
            sigma_ωxsigma_ωx_counter += 1
        end
    end
    Q_constrained = Q_unconstrained * Q_unconstrained'

    # Updates
    # H
    current_name = 1
    for i in 1:p, j in 1:i
        update_constrained_value!(
            model, get_names(model)[current_name], H_constrained[i, j]
        )
        current_name += 1
    end
    # Q
    sigma2_ξ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 1
            update_constrained_value!(model, sigma2_ξ[sigma2_ξ_counter], Q_constrained[pos])
            sigma2_ξ_counter += 1
        end
    end
    sigma2_ζ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 2
            update_constrained_value!(model, sigma2_ζ[sigma2_ζ_counter], Q_constrained[pos])
            sigma2_ζ_counter += 1
        end
    end
    sigma2_ω_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 3
            update_constrained_value!(model, sigma2_ω[sigma2_ω_counter], Q_constrained[pos])
            sigma2_ω_counter += 1
        end
    end
    sigma_ξxsigma_ξx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 1
            update_constrained_value!(
                model, sigma_ξxsigma_ξx[sigma_ξxsigma_ξx_counter], Q_constrained[pos]
            )
            sigma_ξxsigma_ξx_counter += 1
        end
    end
    sigma_ζxsigma_ζx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 2
            update_constrained_value!(
                model, sigma_ζxsigma_ζx[sigma_ζxsigma_ζx_counter], Q_constrained[pos]
            )
            sigma_ζxsigma_ζx_counter += 1
        end
    end
    sigma_ωxsigma_ωx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 3
            update_constrained_value!(
                model, sigma_ωxsigma_ωx[sigma_ωxsigma_ωx_counter], Q_constrained[pos]
            )
            sigma_ωxsigma_ωx_counter += 1
        end
    end
    return model
end

function unconstrain_hyperparameters!(model::MultivariateBasicStructural)
    Fl = typeof_model_elements(model)
    p = size(model.system.y, 2)
    # H
    H_constrained = zeros(Fl, p, p)
    current_name = 1
    for i in 1:p, j in 1:i
        H_constrained[i, j] = get_constrained_value(model, get_names(model)[current_name])
        current_name += 1
    end
    H_unconstrained = cholesky(Symmetric(H_constrained, :L)).L
    # Q
    Q_constrained = zeros(Fl, 3p, 3p)
    sigma2_ξ = filter(x -> occursin("sigma2_ξ", x), get_names(model))
    sigma2_ξ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 1
            Q_constrained[pos] = get_constrained_value(model, sigma2_ξ[sigma2_ξ_counter])
            sigma2_ξ_counter += 1
        end
    end
    sigma2_ζ = filter(x -> occursin("sigma2_ζ", x), get_names(model))
    sigma2_ζ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 2
            Q_constrained[pos] = get_constrained_value(model, sigma2_ζ[sigma2_ζ_counter])
            sigma2_ζ_counter += 1
        end
    end
    sigma2_ω = filter(x -> occursin("sigma2_ω", x), get_names(model))
    sigma2_ω_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 3
            Q_constrained[pos] = get_constrained_value(model, sigma2_ω[sigma2_ω_counter])
            sigma2_ω_counter += 1
        end
    end
    sigma_ξxsigma_ξx = filter(x -> occursin("sigma_ξ", x), get_names(model))
    sigma_ξxsigma_ξx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 1
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ξxsigma_ξx[sigma_ξxsigma_ξx_counter]
            )
            sigma_ξxsigma_ξx_counter += 1
        end
    end
    sigma_ζxsigma_ζx = filter(x -> occursin("sigma_ζ", x), get_names(model))
    sigma_ζxsigma_ζx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 2
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ζxsigma_ζx[sigma_ζxsigma_ζx_counter]
            )
            sigma_ζxsigma_ζx_counter += 1
        end
    end
    sigma_ωxsigma_ωx = filter(x -> occursin("sigma_ω", x), get_names(model))
    sigma_ωxsigma_ωx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 3
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ωxsigma_ωx[sigma_ωxsigma_ωx_counter]
            )
            sigma_ωxsigma_ωx_counter += 1
        end
    end
    Q_unconstrained = cholesky(Symmetric(Q_constrained, :L)).L

    # Updates
    # H
    current_name = 1
    for i in 1:p, j in 1:i
        update_unconstrained_value!(
            model, get_names(model)[current_name], H_unconstrained[i, j]
        )
        current_name += 1
    end
    # Q
    sigma2_ξ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 1
            update_unconstrained_value!(
                model, sigma2_ξ[sigma2_ξ_counter], Q_unconstrained[pos]
            )
            sigma2_ξ_counter += 1
        end
    end
    sigma2_ζ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 2
            update_unconstrained_value!(
                model, sigma2_ζ[sigma2_ζ_counter], Q_unconstrained[pos]
            )
            sigma2_ζ_counter += 1
        end
    end
    sigma2_ω_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, 0)))
        if mod1(i, 3) == 3
            update_unconstrained_value!(
                model, sigma2_ω[sigma2_ω_counter], Q_unconstrained[pos]
            )
            sigma2_ω_counter += 1
        end
    end
    sigma_ξxsigma_ξx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 1
            update_unconstrained_value!(
                model, sigma_ξxsigma_ξx[sigma_ξxsigma_ξx_counter], Q_unconstrained[pos]
            )
            sigma_ξxsigma_ξx_counter += 1
        end
    end
    sigma_ζxsigma_ζx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 2
            update_unconstrained_value!(
                model, sigma_ζxsigma_ζx[sigma_ζxsigma_ζx_counter], Q_unconstrained[pos]
            )
            sigma_ζxsigma_ζx_counter += 1
        end
    end
    sigma_ωxsigma_ωx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_unconstrained, -3)))
        if mod1(i, 3) == 3
            update_unconstrained_value!(
                model, sigma_ωxsigma_ωx[sigma_ωxsigma_ωx_counter], Q_unconstrained[pos]
            )
            sigma_ωxsigma_ωx_counter += 1
        end
    end
    return model
end

function fill_model_system!(model::MultivariateBasicStructural)
    Fl = typeof_model_elements(model)
    p = size(model.system.y, 2)
    H_constrained = zeros(Fl, p, p)
    current_name = 1
    for i in 1:p, j in 1:i
        H_constrained[i, j] = get_constrained_value(model, get_names(model)[current_name])
        current_name += 1
    end
    model.system.H = Matrix(Symmetric(H_constrained, :L))
    # Q
    Q_constrained = zeros(Fl, 3p, 3p)
    sigma2_ξ = filter(x -> occursin("sigma2_ξ", x), get_names(model))
    sigma2_ξ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 1
            Q_constrained[pos] = get_constrained_value(model, sigma2_ξ[sigma2_ξ_counter])
            sigma2_ξ_counter += 1
        end
    end
    sigma2_ζ = filter(x -> occursin("sigma2_ζ", x), get_names(model))
    sigma2_ζ_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 2
            Q_constrained[pos] = get_constrained_value(model, sigma2_ζ[sigma2_ζ_counter])
            sigma2_ζ_counter += 1
        end
    end
    sigma2_ω = filter(x -> occursin("sigma2_ω", x), get_names(model))
    sigma2_ω_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, 0)))
        if mod1(i, 3) == 3
            Q_constrained[pos] = get_constrained_value(model, sigma2_ω[sigma2_ω_counter])
            sigma2_ω_counter += 1
        end
    end
    sigma_ξxsigma_ξx = filter(x -> occursin("sigma_ξ", x), get_names(model))
    sigma_ξxsigma_ξx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 1
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ξxsigma_ξx[sigma_ξxsigma_ξx_counter]
            )
            sigma_ξxsigma_ξx_counter += 1
        end
    end
    sigma_ζxsigma_ζx = filter(x -> occursin("sigma_ζ", x), get_names(model))
    sigma_ζxsigma_ζx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 2
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ζxsigma_ζx[sigma_ζxsigma_ζx_counter]
            )
            sigma_ζxsigma_ζx_counter += 1
        end
    end
    sigma_ωxsigma_ωx = filter(x -> occursin("sigma_ω", x), get_names(model))
    sigma_ωxsigma_ωx_counter = 1
    for (i, pos) in enumerate(collect(diagind(Q_constrained, -3)))
        if mod1(i, 3) == 3
            Q_constrained[pos] = get_constrained_value(
                model, sigma_ωxsigma_ωx[sigma_ωxsigma_ωx_counter]
            )
            sigma_ωxsigma_ωx_counter += 1
        end
    end
    model.system.Q = Matrix(Symmetric(Q_constrained, :L))
    return model
end

function reinstantiate(model::MultivariateBasicStructural, y::Matrix{Fl}) where Fl
    return MultivariateBasicStructural(y, model.seasonality)
end

has_exogenous(::MultivariateBasicStructural) = false


text:
using JSONLines
using Test, DataFrames, JSONTables, Pipe, Tables

full_web = LineIndex("testfiles/jsonlwebsite.jsonl") |> DataFrame;
nrow_fw = nrow(full_web)

mtcars = jsontable(read("testfiles/mtcars.json")) |> DataFrame
full_mtcars = LineIndex("testfiles/mtcars.jsonl"; nworkers = 4) |> DataFrame;
@show full_mtcars
# Fix R export differences
rename!(full_mtcars, :_row => :model);
rename!(full_mtcars, names(full_mtcars) .=> lowercase.(names(full_mtcars)))
rename!(mtcars, names(mtcars) .=> lowercase.(names(mtcars)))

# Read without promotion
noprom_mtcars = LineIndex("testfiles/mtcars.jsonl") |> DataFrame;
nrow_mt = nrow(mtcars)

oneline = LineIndex("testfiles/oneline.jsonl")|> DataFrame;
oneline_plus = LineIndex("testfiles/oneline_plus.jsonl")|> DataFrame;

escaped = LineIndex("testfiles/escapedeol.jsonl")|> DataFrame;

@testset "Read" begin
    @test full_web.name == ["Gilbert", "Alexa", "May", "Deloise"]
    @test full_web.wins[1] == [["straight", "7♣"], ["one pair", "10♥"]]
    @test full_web.wins[end] == [["three of a kind", "5♣"]]
    @test full_mtcars.mpg == mtcars.mpg
    @test noprom_mtcars.cyl[32] == 4
    @test noprom_mtcars.wt[30] == 2.77
    @test noprom_mtcars.qsec[16] == 17.82
    @test noprom_mtcars[!, :_row][16] == "Lincoln Continental"
    @test noprom_mtcars[!, :drat] == mtcars[!, :drat]
    @test oneline.name == ["Daniel"]
    @test oneline_plus.name == ["Daniel"]
    @test nrow(escaped) == 4
    @test escaped.name[1] == "Daniel\n"
    @test escaped.age[2] == "}"
end

@testset "LineIterator" begin
    @test [x for x in LineIterator("testfiles/jsonlwebsite.jsonl")] |> DataFrame == full_web
    @test [x for x in LineIterator("testfiles/mtcars.jsonl")] |> DataFrame == noprom_mtcars
    @test [x for x in LineIterator("testfiles/oneline.jsonl")] |> DataFrame== oneline
    @test [x for x in LineIterator("testfiles/oneline_plus.jsonl")]  |> DataFrame == oneline_plus
    @test [x for x in LineIterator("testfiles/escapedeol.jsonl")]  |> DataFrame == escaped
end

@testset "Read arrays" begin
    @test LineIndex("testfiles/array.jsonl", skip = 1)[2,:a] == 4
    @test LineIndex("testfiles/array.jsonl", skip = 1)[2,1] == 4
    @test LineIndex("testfiles/jsonlwebsitearray.jsonl")[1:end, :Score] == [24, 29, 14, 19]
end

const JL = JSONLines
@testset "readcols" begin
    webl = @pipe JL.readcols("testfiles/jsonlwebsite.jsonl", :name) |> materialize  |> DataFrame
    @test webl == full_web[:, [:name]]
    mtl = @pipe JL.readcols("testfiles/mtcars.jsonl", :gear, :hp; nworkers = 4) |> materialize |> DataFrame
    @test mtl == noprom_mtcars[:, [:gear, :hp]]
    onel = @pipe JL.readcols("testfiles/oneline.jsonl", :age) |> materialize |> DataFrame
    @test onel == oneline[:, [:age]]
    onepl = @pipe JL.readcols("testfiles/oneline_plus.jsonl", :name) |> materialize |> DataFrame
    @test onepl == oneline_plus[:, [:name]]
    escl = @pipe JL.readcols("testfiles/escapedeol.jsonl", :name) |> materialize |> DataFrame
    @test escl == escaped[:, [:name]]
end

@testset "Read nworkers" begin
# full file equal
    @test LineIndex("testfiles/jsonlwebsite.jsonl", nworkers = 2) |> DataFrame == full_web
    @test LineIndex("testfiles/mtcars.jsonl", nworkers = Base.Threads.nthreads() + 1) |> DataFrame == noprom_mtcars
    @test LineIndex("testfiles/oneline.jsonl", nworkers = 1) |> DataFrame== oneline
end

@testset "skip & nrows" begin
# skip + nrows = nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = nrow_fw-1)  |> DataFrame == full_web[2:end, :]
    @test LineIndex("testfiles/mtcars.jsonl", skip = 2, nrows = nrow_mt-2)  |> DataFrame == noprom_mtcars[3:end, :]

# skip + nrows < nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = 2)  |> DataFrame == full_web[2:3, :]
    @test LineIndex("testfiles/mtcars.jsonl", skip = 2, nrows = 2)  |> DataFrame == noprom_mtcars[3:4, :]
    @test LineIndex("testfiles/escapedeol.jsonl", skip = 2, nrows = 1)  |> DataFrame == escaped[3:3, :]

# skip + nrows > nrow(file) (through nrow)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", skip = 1, nrows = nrow_fw)  |> DataFrame == full_web[2:end, :]
    @test LineIndex("testfiles/mtcars.jsonl", skip = 12, nrows = nrow_mt + 10)  |> DataFrame == noprom_mtcars[13:end, :]
    @test LineIndex("testfiles/oneline.jsonl", skip = 0, nrows = 5)  |> DataFrame == oneline
    @test LineIndex("testfiles/oneline_plus.jsonl", skip = 0, nrows = 2)  |> DataFrame == oneline_plus
    @test LineIndex("testfiles/escapedeol.jsonl", skip = 2, nrows = 10)  |> DataFrame == escaped[3:end, :]

# skip + nrows > nrow(file) (through skip)
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/jsonlwebsite.jsonl", skip = nrow_fw+1, nrows = 1))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/mtcars.jsonl", skip = nrow_mt +12, nrows = 120))  |> DataFrame ==  DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline.jsonl", skip = 2, nrows = 10))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline_plus.jsonl", skip = 2, nrows = 123))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/escapedeol.jsonl", skip = 5, nrows = 1))  |> DataFrame == DataFrame()

# skip = nrow(file) + nrows > 0
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/jsonlwebsite.jsonl", skip = nrow_fw, nrows = 10))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/mtcars.jsonl", skip = nrow_mt, nrows = 1))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline.jsonl", skip = 1, nrows = 12)) |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline_plus.jsonl", skip = 1, nrows = 1)) |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/escapedeol.jsonl", skip = 4, nrows = 1)) |> DataFrame == DataFrame()
end

@testset "skip" begin
# skip = nrow(file)
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/jsonlwebsite.jsonl", skip = nrow_fw))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/mtcars.jsonl", skip = nrow_mt))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline.jsonl", skip = 1))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline_plus.jsonl", skip = 1))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/escapedeol.jsonl", skip = 4))  |> DataFrame == DataFrame()

# skip > nrow(file)
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/jsonlwebsite.jsonl", skip = nrow_fw + 1))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/mtcars.jsonl", skip = nrow_mt + 42)) |> DataFrame  == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/mtcars.jsonl", skip = typemax(Int)))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline.jsonl", skip = 2))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/oneline_plus.jsonl", skip = 2))  |> DataFrame == DataFrame()
    @test (@test_logs (:warn, "Skipped all lines") LineIndex("testfiles/escapedeol.jsonl", skip = 5)) |> DataFrame  == DataFrame()

# skip < nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", skip = nrow_fw - 1)  |> DataFrame == full_web[end:end, :]
    @test LineIndex("testfiles/mtcars.jsonl", skip = nrow_mt - 12)  |> DataFrame == noprom_mtcars[(end-11):end, :]
    @test LineIndex("testfiles/escapedeol.jsonl", skip = 2)  |> DataFrame == escaped[3:end, :]
end

@testset "nrows" begin
# nrows < nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", nrows = 2)  |> DataFrame == full_web[begin:2, :]
    @test LineIndex("testfiles/mtcars.jsonl", nrows = 12)  |> DataFrame == noprom_mtcars[begin:12, :]
    @test LineIndex("testfiles/escapedeol.jsonl", nrows = 3)  |> DataFrame == escaped[begin:3, :]

# nrows = nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", nrows = nrow_fw)  |> DataFrame == full_web
    @test LineIndex("testfiles/mtcars.jsonl", nrows = nrow_mt)  |> DataFrame == noprom_mtcars
    @test LineIndex("testfiles/oneline.jsonl", nrows = 1)  |> DataFrame == oneline
    @test LineIndex("testfiles/oneline_plus.jsonl", nrows = 1)  |> DataFrame == oneline_plus
    @test LineIndex("testfiles/escapedeol.jsonl", nrows = 4)  |> DataFrame == escaped

# nrows > nrow(file)
    @test LineIndex("testfiles/jsonlwebsite.jsonl", nrows = nrow_fw+1)  |> DataFrame == full_web
    @test LineIndex("testfiles/mtcars.jsonl", nrows = nrow_mt+100)  |> DataFrame == noprom_mtcars
    @test LineIndex("testfiles/oneline.jsonl", nrows = 2)  |> DataFrame == oneline
    @test LineIndex("testfiles/oneline_plus.jsonl", nrows = 2)  |> DataFrame == oneline_plus
    @test LineIndex("testfiles/escapedeol.jsonl", nrows = 5)  |> DataFrame == escaped
end

writelines("full_web.jsonl", full_web)
writelines("full_mtcars.jsonl", full_mtcars)
writelines("oneline2.jsonl", oneline)
writelines("oneline_plus2.jsonl", oneline_plus)
writelines("escaped2.jsonl", escaped)
@testset "write" begin
    @test LineIndex("full_web.jsonl") |> DataFrame == full_web
    @test LineIndex("full_mtcars.jsonl") |> DataFrame == full_mtcars
    @test LineIndex("oneline2.jsonl") |> DataFrame == oneline
    @test LineIndex("oneline_plus2.jsonl") |> DataFrame == oneline_plus
    @test LineIndex("escaped2.jsonl") |> DataFrame == escaped
end

@testset "MStructType" begin
    @MStructType MyType hp gear drat
    @test LineIndex("testfiles/mtcars.jsonl", structtype = MyType) |> DataFrame == full_mtcars[:, [:hp, :gear, :drat]]
    @MStructType EscType name
    @test LineIndex("testfiles/escapedeol.jsonl", structtype = EscType) |> DataFrame == escaped[:, [:name]]
end

@testset "Columnwise" begin
    mt_cur = LineIndex("testfiles/mtcars.jsonl")
    @test columnwise(mt_cur).mpg == noprom_mtcars.mpg
end

@testset "Filter" begin
    @test filter(row -> row[:mpg] > 20, LineIndex("testfiles/mtcars.jsonl")) |> DataFrame == filter(row -> row[:mpg] > 20, noprom_mtcars)
end

# Cleanup
rm("full_web.jsonl")
rm("full_mtcars.jsonl")
rm("oneline2.jsonl")
rm("oneline_plus2.jsonl")
rm("escaped2.jsonl")


text:
# Spinnaker.jl: wrapper for FLIR/Point Grey Spinnaker SDK
# Copyright (C) 2019 Samuel Powell

# acquistion.jl: replicate Acquisition_C Spinnaker SDK example with low-level library
using Spinnaker

outpath = joinpath(@__DIR__, "output")
mkpath(outpath)

using Spinnaker
import Spinnaker: bool8_t, MAX_BUFFER_LEN

# Get system interface and library version
hSystem = Ref(spinSystem(C_NULL))
spinSystemGetInstance(hSystem)

hLibraryVersion = Ref(spinLibraryVersion(0,0,0,0))
spinSystemGetLibraryVersion(hSystem[], hLibraryVersion)

libver = VersionNumber(hLibraryVersion[].major, hLibraryVersion[].minor, hLibraryVersion[].type)
@info "Spinnaker library version: $libver, build $(hLibraryVersion[].build)"

# Retrieve camera list (from system object)
hCameraList = Ref(spinCameraList(C_NULL))
nCameras = Ref(Csize_t(0))
spinCameraListCreateEmpty(hCameraList)
spinSystemGetCameras(hSystem[], hCameraList[])
spinCameraListGetSize(hCameraList[], nCameras)

@info "Discovered $(nCameras[]) cameras"

function readable(hNode, nodeName)
   pbAvailable = Ref(bool8_t(false))
   pbReadable = Ref(bool8_t(false))
   spinNodeIsAvailable(hNode, pbAvailable)
   spinNodeIsReadable(hNode, pbReadable)
   return (pbReadable[] == 1) && (pbAvailable[] == 1)
end

function writable(hNode, nodeName)
  pbAvailable = Ref(bool8_t(false))
  pbReadable = Ref(bool8_t(false))
  spinNodeIsAvailable(hNode, pbAvailable)
  spinNodeIsWritable(hNode, pbReadable)
  return (pbReadable[] == 1) && (pbAvailable[] == 1)
end

function acquire_images(hCam, hNodeMap, hNodeMapTLDevice)

  # Set acquisition mode
  hAcquisitionMode = Ref(spinNodeHandle(C_NULL))
  hAcquisitionModeContinuous = Ref(spinNodeHandle(C_NULL))
  acquisitionModeContinuous = Ref(Int64(0))

  spinNodeMapGetNode(hNodeMap[], "AcquisitionMode", hAcquisitionMode)
  if readable(hAcquisitionMode[], "AcquisitionMode")
    spinEnumerationGetEntryByName(hAcquisitionMode[], "Continuous", hAcquisitionModeContinuous)
  else
    @error "Unable to retrieve acquistion mode enumeration"
  end

  if readable(hAcquisitionModeContinuous[], "AcquisitionModeContinuous")
    spinEnumerationEntryGetIntValue(hAcquisitionModeContinuous[], acquisitionModeContinuous);
  else
    @error "Unable to retrieve continuous acquistion enumeration"
  end

  if writable(hAcquisitionMode[], "AcquisitionMode")
      spinEnumerationSetIntValue(hAcquisitionMode[], acquisitionModeContinuous[])
  else
      @info "Unable to set continuous acquistion mode"
  end

  @info "Acquistion mode set to continuous"

  # Begin acquistion
  spinCameraBeginAcquisition(hCam[])
  @info "Started camera acquistion"

  # Get serial number
  hDeviceSerialNumber = Ref(spinNodeHandle(C_NULL))
  deviceSerialNumber = Vector{UInt8}(undef, MAX_BUFFER_LEN)
  lenDeviceSerialNumber = Ref(Csize_t(MAX_BUFFER_LEN))

   spinNodeMapGetNode(hNodeMapTLDevice[], "DeviceSerialNumber", hDeviceSerialNumber)
   if readable(hDeviceSerialNumber[], "DeviceSerialNumber")
     spinStringGetValue(hDeviceSerialNumber[], deviceSerialNumber, lenDeviceSerialNumber)
   else
     @info "Cannot receive camera serial number"
   end

   @info "Camera serial number: $(unsafe_string(pointer(deviceSerialNumber)))"

   # Retrieve, convert, and save images
   for imageCnt in 1:10

    hResultImage = Ref(spinImage(C_NULL))
    spinCameraGetNextImage(hCam[], hResultImage);

    isIncomplete = Ref(bool8_t(false))
    hasFailed = Ref(bool8_t(false))

    spinImageIsIncomplete(hResultImage[], isIncomplete);

    if (isIncomplete == true)
      imageStatus = Ref(spinImageStatus(IMAGE_NO_ERROR))
      spinImageGetStatus(hResultImage[], imageStatus)
      @info "Image incomplete with error $(imageStatus)"
      hasFailed = true
    end

    if (hasFailed == true)
      spinImageRelease(hResultImage[]);
    end

    # Get image dimensions
    width = Ref(Csize_t(0))
    height = Ref(Csize_t(0))
    spinImageGetWidth(hResultImage[], width);
    spinImageGetHeight(hResultImage[], height);
    @info "Grabbed image $imageCnt, width: $(width[]), height: $(height[])"

    # Convert image to MONO8
    hConvertedImage = Ref(spinImage(C_NULL))
    spinImageCreateEmpty(hConvertedImage);
    spinImageConvert(hResultImage[], Spinnaker.PixelFormat_Mono8, hConvertedImage[]);
    @info "Converted image to MONO 8"

    # Save to disc
    filename = joinpath(outpath, "acquistion_" * unsafe_string(pointer(deviceSerialNumber)) * string(imageCnt) * ".jpg")
    spinImageSave(hConvertedImage[], filename, Spinnaker.JPEG);
    @info "Image Saved to $filename"

    # Destroy converted image and release camera image
    spinImageDestroy(hConvertedImage[]);
    spinImageRelease(hResultImage[]);

  end

  spinCameraEndAcquisition(hCam[])
  @info "Stopped camera acquistion"

  return nothing
end

function run_camera(hCam)
  hNodeMapTLDevice = Ref(spinNodeMapHandle(C_NULL))
  hNodeMap = Ref(spinNodeMapHandle(C_NULL))
  spinCameraGetTLDeviceNodeMap(hCam[], hNodeMapTLDevice);
  spinCameraInit(hCam[]);
  spinCameraGetNodeMap(hCam[], hNodeMap);
  acquire_images(hCam, hNodeMap, hNodeMapTLDevice);
  spinCameraDeInit(hCam[]);
  return nothing
end

if nCameras[] > 0
  hCam = Ref(spinCamera(C_NULL))
  for i in 0:nCameras[]-1
    spinCameraListGet(hCameraList[], i, hCam);
    run_camera(hCam)
    spinCameraRelease(hCam[])
  end
else
  # Finish if there are no cameras
  @info "No cameras, finishing"
end

spinCameraListClear(hCameraList[])
spinCameraListDestroy(hCameraList[])
spinSystemReleaseInstance(hSystem[])


text:
polymin_n_ds = [
    [
    (1, 4), # compile run
    (1, 100),
    (1, 200),
    (1, 500),
    (1, 1000),
    (1, 2000),
    (1, 3000),
    (1, 4000),
    (1, 5000),
    ],
    [
    (2, 3), # compile run
    (2, 15),
    (2, 30),
    (2, 45),
    (2, 60),
    (2, 75),
    ],
    [
    (3, 2), # compile run
    (3, 6),
    (3, 9),
    (3, 12),
    (3, 15),
    (3, 18),
    ],
    [
    (3, 2), # compile run
    (4, 4),
    (4, 6),
    (4, 8),
    (4, 10),
    ],
    [
    (3, 2), # compile run
    (8, 2),
    (8, 3),
    (8, 4),
    ],
    [
    (3, 2), # compile run
    (16, 1),
    (16, 2),
    ],
    [
    (3, 2), # compile run
    (32, 1),
    ],
    [
    (3, 2), # compile run
    (64, 1),
    ],
    ]
polymin_insts(use_nat::Bool) = [
    [(n, d, false, use_nat) for (n, d) in nds]
    for nds in polymin_n_ds
    ]

insts = OrderedDict()
insts["nat"] = (nothing, polymin_insts(true))
insts["ext"] = (:SOCExpPSD, polymin_insts(false))
return (PolyMinJuMP, insts)


text:
# Defining regions on the plane, testing for inclusion, tessellations, and interpolation

using VoronoiDelaunay  # for creating Delaunay meshes on unstructured data

# Region types

"A region on the plane"
abstract type Region end

"A boxed region [`xmin`, `xmax`] ⊗ [`ymin`, `ymax`]"
struct BoundingBox <: Region
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
end

""" A polygon whose corners are given in the vectors `x` and `y`

    For example, Polygon([0, 1, 1, 0], [0, 0, 1, 1]) defines the unit square.
    """
struct Polygon <: Region
    x::Vector{Float64}
    y::Vector{Float64}
    boundingbox::BoundingBox

    # we precalculate the bounding box since we'll be using it a lot
    function Polygon(x, y)
        @assert (length(x) == length(y)) "Uneven number of x and y points"
        new(x, y, BoundingBox(extrema(x)..., extrema(y)...))
    end
end
"Get the corners of a `Polygon`"
corners(p::Polygon) = zip(p.x, p.y)
BoundingBox(p::Polygon) = p.boundingbox

# Inclusion testing

Base.in(xy, r::Region) = in(xy..., r)
Base.in(x, y, bb::BoundingBox) = (bb.xmin ≤ x ≤ bb.xmax) && (bb.ymin ≤ y ≤ bb.ymax)
function Base.in(x, y, p::Polygon)
    if !((x, y) in p.boundingbox)
        return false
    end

    # O(n) polygon inclusion test - this could definitely be faster
    inpoly = false
    lastcx = p.x[end]
    lastcy = p.y[end]
    for (cx, cy) in corners(p)
        if cy < y && lastcy ≥ y || lastcy < y && cy ≥ y
            if (cx + (y - cy)/(lastcy - cy)*(lastcx - cx)) < x
                inpoly = !inpoly
            end
        end
        lastcx, lastcy = cx, cy
    end

    return inpoly
end

# Tesselation utils: mapping to the interval (1+eps(), 2-2*eps())

"Map `x` from the interval [`lb`, `ub`] to (1.0, 2.0)"
norm12(x, lb, ub) = (x - lb) ./ (ub - lb) * (max_coord - min_coord) + min_coord

"Map each element of the vector `x` to the interval (1.0, 2.0)"
norm12(x::Vector) = norm12(x, extrema(x)...)

"Map `x` from the interval (1.0, 2.0) to [`lb`, `ub`]"
unnorm(x, lb, ub) = (x - min_coord) / (max_coord - min_coord) .* (ub - lb) + lb

"Map the log of `x` from [`lb`, `ub`] to the interval (1.0, 2.0)"
lognorm12(x, lb, ub) = norm12(log10.(x), log10(lb), log10(ub))

"Map the log of the vector `x` to the interval (1.0, 2.0)"
lognorm12(x::Vector) = norm12(log10.(x))

"Map the log of `x` from (1.0, 2.0) back to [`lb`, `ub`]"
unlognorm(x, lb, ub) = 10.^unnorm(x, log10(lb), log10(ub))

# Tessellations

const Tessellation = DelaunayTessellation2D

BoundingBox(t::Tessellation) = BoundingBox(min_coord, max_coord, min_coord, max_coord)

function Base.in(xn, yn, tess::Tessellation)
    if !in(xn, yn, BoundingBox(tess))
        return false
    else
        tri = findtriangle(tess, xn, yn)
        return !isexternal(tri)
    end
end

"Make a Delaunay tessellation on the plane using points in `xs` and `ys`"
function get_tessellation(xs::Vector, ys::Vector; uselog=false)
    xn = uselog ? lognorm12(xs) : norm12(xs)
    yn = uselog ? lognorm12(ys) : norm12(ys)
    tess = DelaunayTessellation(length(xs))
    for (x, y) in zip(xn, yn)
        push!(tess, Point(x, y))
    end

    tess
end

"Find the triangle in the tessellation `tess` containing points `xn` and `yn`"
function findtriangle(t::Tessellation, xn, yn)
    triangle = locate(t, Point(xn, yn))
end

"Get the vertices of triangle `tri`"
getvertices(tri) = map(f -> f(tri), [geta, getb, getc])

"Get barycentric coordinates on the triangle `tri` from cartesian (`xn`, `yn`)"
function barycoords(tri, xn, yn)
    a, b, c = getvertices(tri)
    x = map(getx, [a, b, c])
    y = map(gety, [a, b, c])

    r3 = [x[3], y[3]]
    r = [xn, yn]
    T = [x[1]-x[3]  x[2]-x[3];
         y[1]-y[3]  y[2]-y[3]]
    λ1, λ2 = T \ (r - r3)
    λ = [λ1, λ2, 1 - λ1 - λ2]
end

"Linear interpolation from barycentric coordinates `λ` between vertex values `x`"
lininterp(λ, x) = λ⋅x


text:
# ---
# title: 1262. Greatest Sum Divisible by Three
# id: problem1262
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/greatest-sum-divisible-by-three/description/>
# hidden: true
# ---
#
# Given an array `nums` of integers, we need to find the maximum possible sum of
# elements of the array such that it is divisible by three.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [3,6,5,1,8]
#     Output: 18
#     Explanation: Pick numbers 3, 6, 1 and 8 their sum is 18 (maximum sum divisible by 3).
#
# **Example 2:**
#
#
#
#     Input: nums = [4]
#     Output: 0
#     Explanation: Since 4 is not divisible by 3, do not pick any number.
#
#
# **Example 3:**
#
#
#
#     Input: nums = [1,2,3,4,4]
#     Output: 12
#     Explanation: Pick numbers 1, 3, 4 and 4 their sum is 12 (maximum sum divisible by 3).
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 4 * 10^4`
#   * `1 <= nums[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Particles: marginal_log_prior, log_prior, candidates, simulate, add

@testset "Priors on states" begin
    @testset "Chinese restaurant process" begin

        crp = ChineseRestaurantProcess(0.5)
        @test crp.N == Float64[]
        @test crp.α == 0.5

        @test candidates(crp) == 1:1

        crp1, _ = add(crp, 1)
        @test crp1.N == [1.]
        @test crp1.α == crp.α
        @test candidates(crp1) == 1:2

        crp11, _ = add(crp1, 1)
        @test crp11.N == [2.]

        crp111, _ = add(crp11, 1, 0.5)
        @test crp111.N == [2.5]

        @test_throws BoundsError add(crp, 2)
        @test_throws BoundsError add(crp1, 3)

        crp1112, _ = add(crp111, 2)
        @test crp1112.N == [2.5, 1.]
        @test candidates(crp1112) == 1:3

    end

    @testset "Sticky CRP" begin

        @testset "StickyCRP with ρ=0 == CRP" begin
            scrp = StickyCRP(0.5, 0.0)
            crp = ChineseRestaurantProcess(0.5)

            Random.seed!(1)
            scrp_sim, scrp_states = simulate(scrp, 100)
            Random.seed!(1)
            crp_sim, crp_states = simulate(crp, 100)

            @test scrp_states == crp_states
            @test scrp_sim.N == crp_sim.N

            @test log_prior.(Ref(crp_sim), candidates(crp_sim)) ==
                log_prior.(Ref(scrp_sim), candidates(scrp_sim))[2:end]

        end

        crp = StickyCRP(0.5, 0.9)
        crp1, s1 = add(crp, 1, 0.5)
        crp0, s0 = add(crp, 0, 0.5)

        # state is same after `add` with sticky or not
        @test s1 == s0
        @test crp1.N == [0.5]
        # ... but counts are not.
        @test crp0.N == [0.]

    end
end


text:
using SIUnits
using SIUnits.ShortUnits
using Base.Test

if VERSION <= v"0.3-"
    macro test_throws_compat(a,b)
        esc(:(@test_throws $b))
    end
else
    macro test_throws_compat(a,b)
        esc(:(@test_throws $a $b))
    end
end

if VERSION <= v"0.4.0-dev+3338"
    const AssertionError = ErrorException
end

# Test types

for x in (1,1.0f0,1.0,complex(1.0,0.0))
    @test typeof(x*Meter) == quantity(typeof(x),Meter)
    @test SIUnits.unit(x*Meter) == Meter
end

# Basic arithmetic things
@test V*5 == 5*V
@test -(1V) == (-1)V
@test 1V + 2V == 3V
@test (1//2)V - 1V == (-1//2)V
@test_throws_compat ErrorException 1//2V - 1V
@test_throws_compat ErrorException 1V + 2s + 2kg
@test 1V + zero(1V) == 1V
@test 1V + zero(typeof(1V)) == 1V
@test 1V * one(1V) == 1V
@test 1V * one(typeof(1V)) == 1V
@test 3V / 3V == one(3)
@test 3.V / 3.V === one(3.)
@test isa(3.V / 3.V, typeof(one(3.)))
@test isa(one(3.V), typeof(one(3.)))

OneNewton = 1*(kg*m/s^2)
@test OneNewton*(1s)^2 == 1kg*m
@test OneNewton*s^2 == 1*kg*m

@test OneNewton/(kg*m) == 1Hertz^2

@test OneNewton^2 == 1kg^2*m^2/s^4

@test s == sqrt(1s^2)
@test sqrt(2)*m == sqrt(2m^2)

@test_throws_compat InexactError sqrt(1s)

@test 1/s == 1Hz

@test 1Hz*1s == 1
@test 1s/(1s) == 1

@test (1.23V)^0 == one(1.23)
@test (1.23V)^(0//1) == one(1.23)
@test (1.23V)^(0s/(1s)) == one(1.23)

# Issue #2

immutable note{T<:Real}
    pitch::quantity(T,Hz)     #has dimensions of inverse time
    duration::quantity(T,s) #has dimensions of time
    sustained::Bool
end

note{Float64}(1.0Hz,1.0s,true)

@test_throws_compat ErrorException immutable foo{T}
    bar::quantity(T,2s)
end

# Non-SI Units
@test_approx_eq_eps 3mtorr 3*SIUnits.Milli*as(1torr,Pascal) eps()*Pa

# Unitful Powers
@test_throws_compat ErrorException 5.0^(3Pa)
@test_throws_compat ErrorException 5.0^(3torr)
@test 1^((5torr)/3Pa) == 1.0

# Interaction of si and non-si units
@test (1mm)^2*1torr == as(1torr,Pascal)*(1mm)^2
@test_throws_compat AssertionError as(1torr,s)

# Ranges (#4)
r1 = 1Hz:5Hz
@test length(r1) == 5

@test collect(1Hz:5Hz) == collect(1:5)Hz # Tests the iteration protocol
@test r1[1] == 1Hz # Test indexing

@test r1[2:4] == 2Hz:4Hz # Indexing ranges by ranges should be SIRanges

# Ranges: Ensure forwarded methods are working appropriately
@test eltype(r1) == typeof(1Hz)
@test first(r1) == step(r1) == 1Hz
@test last(r1) == 5Hz
@test r1 - 1Hz == 0Hz:4Hz
@test 1Hz - r1 == 0Hz:-1Hz:-4Hz
@test r1 + 1Hz == 1Hz + r1 == 2Hz:6Hz
@test r1 * 2s == 2s * r1 == 2:2:10
@test r1 * 1m == 1m * r1 == 1Hz*m:5Hz*m
@test r1 / 2Hz == 0.5:0.5:2.5

# Others

@test isnan(1.23K) == false
@test isnan(NaN*K) == true
@test isreal(1.23K) == true
@test isreal(1.23im*K) == false
@test real(1.0K + 2.0im*K) == 1.0K
@test imag(1.0K + 2.0im*K) == 2.0K
@test conj(1.0K + 2.0im*K) == 1.0K - 2.0im*K

@test mod(2µm,4µm) == 2µm
@test_throws_compat ErrorException mod(3kg,4s)

# Issue #9
a = [1m 2N]
b = [1m 2N 3V]
@test a[1:2] == b[1:2]

# Issue #10
a = 1m
b = 2m
@test sqrt(a*b) == (a*b)^(1/2) == (a*b)^(1//2) == sqrt(2)*m

# Issue #11
au=[1m 2m 3m]
bu=[2N 3N 4N]
@test au*bu' == bu*au'
@test (au*bu')[1] == dot(vec(au),vec(bu))

# Irrationals
@test (3pi)*m == 3*(pi*m)
a = 1m
b = 2mm
@test (pi*a^2)/b == (pi/2)*(m^2/mm)

# Comparisons with numbers and different units
a = SIUnits.UnitQuantity{Float64}(3.0)
@test a < 4
@test 2.8 < a
@test_throws_compat MethodError 1m < 2kg

# Issue #27
@test ([1s]/(1s))[1] == 1
@test ([1s]/(1.0s))[1] == 1
@test ([1.0s]/(1s))[1] == 1
@test ([1.0s]/(1.0s))[1] == 1

@test ([1s]*(1s))[1] == 1s*s
@test ([1s]*(1.0s))[1] == 1s*s
@test ([1.0s]*(1s))[1] == 1s*s
@test ([1.0s]*(1.0s))[1] == 1s*s

# Issue #49
@test m != kg != s != A != K != mol != rad != sr

# Issue #52
@test show(IOBuffer(), 1*Meter) == nothing

# Test angular units
@test 1rad + 2rad == 3rad
@test 1sr  + 2sr  == 3sr
for func in (sin,cos,tan,cot,sec,csc)
    @test func(1.23rad) == func(1.23)
    @test_approx_eq func(1.23deg) func(as(1.23deg,rad))
end

# Test printing and superscript()
@test repr(SIUnits.SIQuantity{Float64,zeros(8)...,0.4}(1.0)) == "1.0 sr⁰･⁴"


text:
module Dierckx

export Spline1D,
       Spline2D,
       evaluate,
       evalgrid,
       get_knots,
       get_coeffs,
       get_residual

const ddierckx = joinpath(dirname(@__FILE__),
                          "../deps/src/ddierckx/libddierckx")

# Ensure library is available.
if (dlopen_e(ddierckx) == C_NULL)
    error("Dierckx not properly installed. Run Pkg.build(\"Dierckx\")")
end

# ----------------------------------------------------------------------------
# 1-d splines

const _fit1d_messages = [
2=>
"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small.
There is an approximation returned but the corresponding weighted sum
of squared residuals does not satisfy the condition abs(fp-s)/s <
tol.""",
3=>
"""The maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small. There is an approximation returned but the corresponding
weighted sum of squared residuals does not satisfy the condition
abs(fp-s)/s < tol.""",
10=>
"""Error on entry, no approximation returned. The following conditions
must hold:
1<=k<=5
x[1] < x[2] < ... < x[end]
w[i] > 0.0 for all i

Additionally, if spline knots are given:
length(xknots) <= length(x) + k + 1
x[1] < xknots[1] < xknots[k+2] < ... < xknots[end] < x[end]
The schoenberg-whitney conditions: there must be a subset of data points
xx[j] such that t[j] < xx[j] < t[j+k+1] for j=1,2,...,n-k-1"""]

const _eval1d_messages = [
1=>
"""Input point out of range""",
10=>
"""Invalid input data. The following conditions must hold:
length(x) != 0 and xb <= x[1] <= x[2] <= ... x[end] <= xe"""]

_translate_bc(bc::String) = (bc == "extrapolate" ? 0 :
                             bc == "zero" ? 1 :
                             bc == "error" ? 2 :
                             bc == "nearest" ? 3 :
                             error("unknown boundary condition: \"$bc\""))

type Spline1D
    t::Vector{Float64}
    c::Vector{Float64}
    k::Int
    bc::Int
    fp::Float64
end

get_knots(spl::Spline1D) = spl.t[spl.k+1:end-spl.k]
get_coeffs(spl::Spline1D) = spl.c[1:end-spl.k+1]
get_residual(spl::Spline1D) = spl.fp

function Spline1D(x::Vector{Float64}, y::Vector{Float64};
                  w::Vector{Float64}=ones(length(x)),
                  k::Int=3, s::Real=0.0, bc::String="nearest")
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("length(x) must be greater than k")

    # outputs
    nest = m+k+1
    n = Array(Int32, 1)
    t = Array(Float64, nest)
    c = Array(Float64, nest)
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # working space
    lwrk = m * (k+1) + nest * (7 + 3*k)
    wrk = Array(Float64, lwrk)
    iwrk = Array(Int32, nest)

    ccall((:curfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # x, y, w
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Float64},  # k, s
           Ptr{Int32}, Ptr{Int32}, # nest, n
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # t, c, fp
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32},  # wrk, lwrk, iwrk
           Ptr{Int32}),  # ier
          &0, &m, x, y, w, &x[1], &x[end], &k, &float64(s), &nest,
          n, t, c, fp, wrk, &lwrk, iwrk, ier)

    ier[1] <= 0 || error(_fit1d_messages[ier[1]])

    # resize output arrays
    resize!(t, n[1])
    resize!(c, n[1] - k - 1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[1])
end

# version with user-supplied knots
function Spline1D(x::Vector{Float64}, y::Vector{Float64},
                  xknots::Vector{Float64};
                  w::Vector{Float64}=ones(length(x)),
                  k::Int=3, bc::String="nearest")
    m = length(x)
    length(y) == m || error("length of x and y must match")
    length(w) == m || error("length of x and w must match")
    m > k || error("k must be less than length(x)")

    # x knots
    # (k+1) knots will be added on either end of interior knots.
    n = length(xknots) + 2 * (k+1)
    t = Array(Float64, n)  # All knots
    for i in 1:length(xknots)
        t[i+k+1] = xknots[i]
    end

    # outputs
    c = Array(Float64, n)
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # working space
    lwrk = m * (k+1) + n * (7 + 3*k)
    wrk = Array(Float64, lwrk)
    iwrk = Array(Int32, n)

    ccall((:curfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # x, y, w
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Float64},  # k, s
           Ptr{Int32}, Ptr{Int32}, # nest, n
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # t, c, fp
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32},  # wrk, lwrk, iwrk
           Ptr{Int32}),  # ier
          &(-1), &m, x, y, w, &x[1], &x[end], &k, &(-1.0), &n,
          &n, t, c, fp, wrk, &lwrk, iwrk, ier)

    ier[1] <= 0 || error(_fit1d_messages[ier[1]])
    resize!(c, n-k-1)

    return Spline1D(t, c, k, _translate_bc(bc), fp[1])
end

function evaluate(spline::Spline1D, x::Vector{Float64})
    m = length(x)
    y = Array(Float64, m)
    ier = Array(Int32, 1)
    ccall((:splev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32},  # x, y, m
           Ptr{Int32}, Ptr{Int32}),  # e, ier
          spline.t, &length(spline.t), spline.c, &spline.k,
          x, y, &m, &spline.bc, ier)
    ier[1] == 0 || error(_eval1d_messages[ier[1]])
    return y
end

function evaluate(spline::Spline1D, x::FloatingPoint)
    y = Array(Float64, 1)
    ier = Array(Int32, 1)
    ccall((:splev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # t, n
           Ptr{Float64}, Ptr{Int32},  # c, k
           Ptr{Float64}, Ptr{Float64}, Ptr{Int32},  # x, y, m
           Ptr{Int32}, Ptr{Int32}),  # e, ier
          spline.t, &length(spline.t), spline.c, &spline.k,
          &x, y, &1, &spline.bc, ier)
    ier[1] == 0 || error(_eval1d_messages[ier[1]])
    return y[1]
end

# ----------------------------------------------------------------------------
# 2-d splines

# NOTE REGARDING ARGUMENT ORDER: In the "grid" version of the Spline2D
# constructor and evaluators, the fortran functions expects z to have
# shape (my, mx), but we'd rather have x be the fast axis in z.  So,
# in the ccall()s in these methods, all the x and y related inputs are
# swapped with regard to what the Fortran documentation says.

const _fit2d_messages = [
-3=>

"""The coefficients of the spline returned have been computed as the
minimal norm least-squares solution of a (numerically) rank deficient
system (deficiency=%i). If deficiency is large, the results may be
inaccurate. Deficiency may strongly depend on the value of eps.""",

1=>

"""The required storage space exceeds the available storage space:
nxest or nyest too small, or s too small.  The weighted least-squares
spline corresponds to the current set of knots.""",

2=>

"""A theoretically impossible result was found during the iteration
process for finding a smoothing spline with fp = s: s too small or
badly chosen eps.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

3=>

"""the maximal number of iterations maxit (set to 20 by the program)
allowed for finding a smoothing spline with fp=s has been reached: s
too small.  Weighted sum of squared residuals does not satisfy
abs(fp-s)/s < tol.""",

4=>

"""No more knots can be added because the number of b-spline
coefficients (nx-kx-1)*(ny-ky-1) already exceeds the number of data
points m: either s or m too small.  The weighted least-squares spline
corresponds to the current set of knots.""",

5=>

"""No more knots can be added because the additional knot would
(quasi) coincide with an old one: s too small or too large a weight to
an inaccurate data point.  The weighted least-squares spline
corresponds to the current set of knots.""",

10=>

"""Error on entry, no approximation returned. The following conditions
must hold:
xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
If iopt==-1, then
  xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
  yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye"""]

const _eval2d_message = (
"""Invalid input data. Restrictions:
length(x) != 0, length(y) != 0
x[i-1] <= x[i] for i=2,...,length(x)
y[j-1] <= y[j] for j=2,...,length(y)
""")

type Spline2D
    tx::Vector{Float64}
    ty::Vector{Float64}
    c::Vector{Float64}
    kx::Int
    ky::Int
    fp::Float64
end

get_knots(spl::Spline2D) = (spl.tx[spl.kx+1:end-spl.kx],
                            spl.ty[spl.ky+1:end-spl.ky])
get_residual(spl::Spline2D) = spl.fp

# Helper function for calculating required size of work array in surfit.
# Note that x and y here are as in the Fortran documentation.
function calc_surfit_lwrk1(m, kx, ky, nxest, nyest)
    u = nxest - kx - 1
    v = nyest - ky - 1
    km = max(kx, ky) + 1
    ne = max(nxest, nyest)
    bx = kx*v + ky + 1
    by = ky*u + kx + 1
    b1 = b2 = 0
    if (bx<=by)
        b1 = bx
        b2 = bx + v - ky
    else
        b1 = by
        b2 = by + u - kx
    end
    return u*v*(2 + b1 + b2) + 2*(u+v+km*(m+ne)+ne-kx-ky) + b2 + 1
end

# Construct spline from unstructured data
function Spline2D(x::Vector{Float64}, y::Vector{Float64}, z::Vector{Float64};
                  w::Vector{Float64}=ones(length(x)), kx::Int=3, ky::Int=3,
                  s::Real=0.0)

    # arary sizes
    m = length(x)
    (length(y) == length(z) == m) || error("lengths of x, y, z must match")
    (length(w) == m) || error("length of w must match other inputs")

    nxest = max(kx+1+iceil(sqrt(m/2)), 2*(kx+1))
    nyest = max(ky+1+iceil(sqrt(m/2)), 2*(ky+1))
    nmax = max(nxest, nyest)

    eps = 1.0e-16

    # bounds
    xb = minimum(x)
    xe = maximum(x)
    yb = minimum(y)
    ye = maximum(y)

    # return values
    nx = Array(Int32, 1)
    tx = Array(Float64, nxest)
    ny = Array(Int32, 1)
    ty = Array(Float64, nyest)
    c = Array(Float64, (nxest-kx-1) * (nyest-ky-1))
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # work arrays
    # Note: in lwrk1, x and y are swapped on purpose.
    lwrk1 = calc_surfit_lwrk1(m, ky, kx, nyest, nxest)
    lwrk2 = 1
    kwrk = m + (nxest - 2*kx - 1) * (nyest - 2*ky - 1)
    wrk1 = Array(Float64, lwrk1)
    wrk2 = Array(Float64, lwrk2)
    iwrk = Array(Int32, kwrk)

    ccall((:surfit_, ddierckx), Void,
          (Ptr{Int32}, Ptr{Int32},  # iopt, m
           Ptr{Float64}, Ptr{Float64},  # y, x
           Ptr{Float64}, Ptr{Float64},  # z, w
           Ptr{Float64}, Ptr{Float64},  # yb, ye
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32},  # nyest, nxest, nmax
           Ptr{Float64},  # eps
           Ptr{Int32}, Ptr{Float64},  # ny, tx
           Ptr{Int32}, Ptr{Float64},  # nx, tx
           Ptr{Float64}, Ptr{Float64},  # c, fp
           Ptr{Float64}, Ptr{Int32},  # wrk1, lwrk1
           Ptr{Float64}, Ptr{Int32},  # wrk2, lwrk2
           Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),   # iwrk, kwrk, ier
          &0, &m, y, x, z, w, &yb, &ye, &xb, &xe, &ky, &kx, &float64(s),
          &nyest, &nxest, &nmax, &eps, ny, ty, nx, tx, c, fp,
          wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, ier)

    if ier[1] > 10
        # lwrk2 is too small, i.e., there is not enough workspace
        # for computing the minimal least-squares solution of a rank
        # deficient system of linear equations. ier gives the
        # requested value for lwrk2. Rerun with that value in "continue"
        # mode, with iopt = 1.
        lwrk2 = ier[1]
        resize!(wrk2, lwrk2)
        ccall((:surfit_, ddierckx), Void,
              (Ptr{Int32}, Ptr{Int32},  # iopt, m
               Ptr{Float64}, Ptr{Float64},  # y, x
               Ptr{Float64}, Ptr{Float64},  # z, w
               Ptr{Float64}, Ptr{Float64},  # yb, ye
               Ptr{Float64}, Ptr{Float64},  # xb, xe
               Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
               Ptr{Int32}, Ptr{Int32}, Ptr{Int32},  # nyest, nxest, nmax
               Ptr{Float64},  # eps
               Ptr{Int32}, Ptr{Float64},  # ny, tx
               Ptr{Int32}, Ptr{Float64},  # nx, tx
               Ptr{Float64}, Ptr{Float64},  # c, fp
               Ptr{Float64}, Ptr{Int32},  # wrk1, lwrk1
               Ptr{Float64}, Ptr{Int32},  # wrk2, lwrk2
               Ptr{Int32}, Ptr{Int32}, Ptr{Int32}),   # iwrk, kwrk, ier
              &1, &m, y, x, z, w, &yb, &ye, &xb, &xe, &ky, &kx, &s,
              &nyest, &nxest, &nmax, &eps, ny, ty, nx, tx, c, fp,
              wrk1, &lwrk1, wrk2, &lwrk2, iwrk, &kwrk, ier)
    end

    if (ier[1] == 0 || ier[1] == -1 || ier[1] == -2)
        # good values, pass.
    elseif ier[1] < -2
        warn("""
        The coefficients of the spline returned have been
        computed as the minimal norm least-squares solution of a
        (numerically) rank deficient system. The rank is $(-ier[1]).
        The rank deficiency is $((nx[1]-kx-1)*(ny[1]-ky-1)+ier[1]).
        Especially if the rank deficiency is large the results may
        be inaccurate.""")
        # "The results could also seriously depend on the value of
        # eps" (not in message because eps is currently not an input)
    else
        error(_fit2d_messages[ier[1]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[1])
    resize!(ty, ny[1])
    resize!(c, (nx[1] - kx - 1) * (ny[1] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[1])
end

# Construct spline from data on a grid.
function Spline2D(x::Vector{Float64}, y::Vector{Float64}, z::Array{Float64,2};
                  kx::Int=3, ky::Int=3, s::Real=0.0)
    mx = length(x)
    my = length(y)
    assert(size(z, 1) == mx && size(z, 2) == my)

    mx > kx || error("length(x) must be greater than kx")
    my > ky || error("length(y) must be greater than ky")

    # Bounds
    xb = x[1]
    xe = x[end]
    yb = y[1]
    ye = y[end]
    nxest = mx+kx+1
    nyest = my+ky+1

    # Return values
    nx = Array(Int32, 1)
    tx = Array(Float64, nxest)
    ny = Array(Int32, 1)
    ty = Array(Float64, nyest)
    c = Array(Float64, (nxest-kx-1) * (nyest-ky-1))
    fp = Array(Float64, 1)
    ier = Array(Int32, 1)

    # Work arrays.
    # Note that in lwrk, x and y are swapped with respect to the Fortran
    # documentation. See "NOTE REGARDING ARGUMENT ORDER" above.
    lwrk = (4 + nyest * (mx+2*ky+5) + nxest * (2*kx+5) +
            my*(ky+1) + mx*(kx+1) + max(mx, nyest))
    wrk = Array(Float64, lwrk)
    kwrk = 3 + mx + my + nxest + nyest
    iwrk = Array(Int32, kwrk)

    ccall((:regrid_, ddierckx), Void,
          (Ptr{Int32},  # iopt
           Ptr{Int32}, Ptr{Float64},  # my, y
           Ptr{Int32}, Ptr{Float64},  # mx, x
           Ptr{Float64},  # z
           Ptr{Float64}, Ptr{Float64},  # yb, ye
           Ptr{Float64}, Ptr{Float64},  # xb, xe
           Ptr{Int32}, Ptr{Int32}, Ptr{Float64},  # ky, kx, s
           Ptr{Int32}, Ptr{Int32},  # nyest, nxest
           Ptr{Int32}, Ptr{Float64},  # ny, ty
           Ptr{Int32}, Ptr{Float64},  # nx, tx
           Ptr{Float64}, Ptr{Float64},  # c, fp
           Ptr{Float64}, Ptr{Int32},  # wrk, lwrk
           Ptr{Int32}, Ptr{Int32},  # iwrk, lwrk
           Ptr{Int32}),  # ier
          &0f0, &my, y, &mx, x, z, &yb, &ye, &xb, &xe, &ky, &kx, &float64(s),
          &nyest, &nxest, ny, ty, nx, tx, c, fp,
          wrk, &lwrk, iwrk, &kwrk, ier)

    if !(ier[1] == 0 || ier[1] == -1 || ier[1] == -2)
        error(_fit2d_messages[ier[1]])
    end

    # Resize output arrays to the size actually used.
    resize!(tx, nx[1])
    resize!(ty, ny[1])
    resize!(c, (nx[1] - kx - 1) * (ny[1] - ky - 1))

    return Spline2D(tx, ty, c, kx, ky, fp[1])
end

# Evaluate spline at individual points
function evaluate(spline::Spline2D, x::Vector{Float64},
                  y::Vector{Float64})
    m = length(x)
    assert(length(y) == m)

    ier = Array(Int32, 1)
    lwrk = spline.kx + spline.ky + 2
    wrk = Array(Float64, lwrk)
    z = Array(Float64, m)

    ccall((:bispeu_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # y, x, z
           Ptr{Int32},  # m
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),  # wrk, lwrk, ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, y, x, z, &m,
          wrk, &lwrk, ier)

    ier[1] == 0 || error(_eval2d_message)

    return z
end

function evaluate(spline::Spline2D, x::Float64, y::Float64)
    ier = Array(Int32, 1)
    lwrk = spline.kx + spline.ky + 2
    wrk = Array(Float64, lwrk)
    z = Array(Float64, 1)
    ccall((:bispeu_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Float64}, Ptr{Float64},  # y, x, z
           Ptr{Int32},  # m
           Ptr{Float64}, Ptr{Int32}, Ptr{Int32}),  # wrk, lwrk, ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, &y, &x, z, &1,
          wrk, &lwrk, ier)

    ier[1] == 0 || error(_eval2d_message)
    return z[1]
end

# Evaluate spline on the grid spanned by the input arrays.
function evalgrid(spline::Spline2D, x::Vector{Float64},
                  y::Vector{Float64})
    mx = length(x)
    my = length(y)

    lwrk = mx*(spline.kx + 1) + my*(spline.ky + 1)
    wrk = Array(Float64, lwrk)
    kwrk = mx + my
    iwrk = Array(Int32, kwrk)
    ier = Array(Int32, 1)
    z = Array(Float64, mx, my)

    ccall((:bispev_, ddierckx), Void,
          (Ptr{Float64}, Ptr{Int32},  # ty, ny
           Ptr{Float64}, Ptr{Int32},  # tx, nx
           Ptr{Float64},  # c
           Ptr{Int32}, Ptr{Int32},  # ky, kx
           Ptr{Float64}, Ptr{Int32},  # y, my
           Ptr{Float64}, Ptr{Int32},  # x, mx
           Ptr{Float64},  # z
           Ptr{Float64}, Ptr{Int32},  # wrk, lwrk
           Ptr{Int32}, Ptr{Int32},  # iwrk, kwrk
           Ptr{Int32}),  # ier
          spline.ty, &length(spline.ty),
          spline.tx, &length(spline.tx),
          spline.c, &spline.ky, &spline.kx, y, &my, x, &mx, z,
          wrk, &lwrk, iwrk, &kwrk, ier)

    ier[1] == 0 || error(_eval2d_message)

    return z
end

end # module


text:
# This file defines the MultiRunSimulation type. This type bundles all the
#   information of a simulation that runs multiple times.

export MultiRunSimulation, MRS
mutable struct MultiRunSimulation

    mpSim::MPsim
    nRuns::Int
    resultsDBname::String
    maxThreads::Int
    showInfo::Bool

    seedRNG::MersenneTwister
    resultsDB::SQLite.DB
    nComplete::Int

    function MultiRunSimulation( mpSim::MPsim )::MRS

        newMRS = new()
        newMRS.mpSim = deepcopy( mpSim )
        newMRS.nRuns = 1
        newMRS.resultsDBname = ""
        newMRS.maxThreads = Threads.nthreads()
        newMRS.showInfo = false
        newMRS.seedRNG = MersenneTwister()
        newMRS.resultsDB = SQLite.DB( "" )
        newMRS.nComplete = 0
        return newMRS

    end  # MultiRunSimulation( mpSim )

    MultiRunSimulation( name::String="sim" )::MRS = MRS( MPsim( name ) )

end  # MultiRunSimulation

const MRS = MultiRunSimulation


text:
type Image
  name::String
  image::Array{Float64}
  answer::Int
end

type Perceptron
  percept_type::Int
  weights::Array{Float64}
end

@enum PERCEPT_TYPES Happy = 1 Sad = 2 Mischevous = 3 Mad = 4
const LEARNING_RATE = 0.1
const THRESHOLD = 1

function parseline(line)
  values = filter(val -> !isempty(val), split(line, " "))
  map(x -> parse(Int, x)/32, values)
end

function getimages(file = "data/training.txt")
  images = []
  open(file) do f
    tmpimg = Image("", [], 0);
    for line in eachline(f)
      if contains(line, "Image")
        tmpimg = Image(normalize_string(line, stripcc=true), [], 0)

      elseif !isempty(tmpimg.name) && (length(line) == 1)
        append!(tmpimg.image, [rand()])
        append!(images, [tmpimg])

      elseif !startswith(line, "#") && (length(line) > 1)
        append!(tmpimg.image, parseline(normalize_string(line, stripcc=true)))

      end
    end
  end
  return images
end

function addfacit!(images, file = "data/training-facit.txt")
  open(file) do f
    lines = split(readstring(f), "\n")
    lines = filter(line -> !startswith(line, "#") && length(line) > 1, lines)
    for i in range(1, (length(lines) - 1))
      if length(lines[i]) > 1 && !startswith(lines[i], "#")
        images[i].answer = parse(Int, lines[i][end])
      end
    end
  end
end

function sigmoid(x)
  1.0 / (1.0 + e^(-x))
end

function process(perceptron, image)
  out = map(*, perceptron.weights, image.image)
  sigmoid(sum(out))
end

function update_weights!(perceptron, dw)
  perceptron.weights = map((x,y) -> (x+y), perceptron.weights, dw)
end

function accurate(errors)
  if length(errors) == 0
    return false
  end
  mse = (sum(map(x -> x^2, errors)))/2
  mse < THRESHOLD
end

function expose!(image, perceptron)
  err = Int(perceptron.percept_type == image.answer) - process(perceptron, image)
  dw = map(x -> (LEARNING_RATE * err * x), image.image)
  update_weights!(perceptron, dw)
  return err
end

function train!(perceptrons, images)
  errors = []
  while(!accurate(errors))
    empty!(errors)
    for perceptron in perceptrons
      shuffled_images = shuffle(images)
      for image in shuffled_images
        append!(errors, expose!(image, perceptron))
      end
    end
  end
end

function examine(perceptrons, images)
  correct = 0

  for image in images
    activation_max = 0
    p_type = 0
    for perceptron in perceptrons
      activation = process(perceptron, image)
      if activation > activation_max
        activation_max = activation
        p_type = perceptron.percept_type
      end
    end
    if image.answer == p_type
      correct += 1
    end
  end

  println("Got $correct out of $(length(images)) correct. $(round(correct/length(images),2)*100)%")
end

function main()
  training_set = getimages("data/FaceTest/training-A.txt")
  addfacit!(training_set, "data/FaceTest/facit-A.txt")
  examine_set = getimages("data/FaceTest/test-B.txt")
  addfacit!(examine_set, "data/FaceTest/facit-B.txt")

  happy = Perceptron(Int(Happy), rand!(Array{Float64}(401)))
  sad = Perceptron(Int(Sad), rand!(Array{Float64}(401)))
  mischevous = Perceptron(Int(Mischevous), rand!(Array{Float64}(401)))
  mad = Perceptron(Int(Mad), rand!(Array{Float64}(401)))
  perceptrons = (happy, sad, mischevous, mad)

  train!(perceptrons, training_set)
  examine(perceptrons, examine_set)
end

function runner()
  println("Starting")
  for i in range(1, 100)
    main()
  end
end

runner()

max_stars_repo_name: linuslagerhjelm/learning-every-language

max_issues_repo_name: linuslagerhjelm/learning-every-language

max_forks_repo_name: linuslagerhjelm/learning-every-language


text:
module KissABC

macro reexport(ex) # taken from Reexport.jl
    modules = Any[e.args[end] for e in ex.args]
    esc(Expr(:toplevel, ex, [:(eval(Expr(:export, names($mod)...))) for mod in modules]...))
end

import AbstractMCMC
import AbstractMCMC: sample, step, MCMCThreads, MCMCDistributed
using Random
import Base.length

using FLoops # added by ML

@reexport using MonteCarloMeasurements
@reexport using Distributions

include("priors.jl")
include("types.jl")
include("transition.jl")
include("smc.jl")
include("abcde_floop.jl") # added by ML

struct AIS <: AbstractMCMC.AbstractSampler
    nparticles::Int
end

struct AISState{S,L}
    "Sample of the Affine Invariant Sampler."
    sample::S
    "Log-likelihood of the sample."
    loglikelihood::L
    "Current particle"
    i::Int
    AISState(s::S, l::L, i = 1) where {S,L} = new{S,L}(s, l, i)
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::AbstractMCMC.AbstractModel,
    spl::AIS;
    retry_sampling::Int = 100,
    kwargs...,
)
    nparticles = spl.nparticles
    nparticles < length(model) + 5 && error(
        "nparticles = ",
        nparticles,
        " is insufficient, set number of particles in AIS(⋅) atleast to ",
        length(model) + 5,
    )

    particles = [op(float, unconditional_sample(rng, model)) for i = 1:nparticles]
    logdensity = [loglike(model, push_p(model, particles[i])) for i = 1:nparticles]
    retrys = retry_sampling * nparticles
    for i = 1:nparticles
        while !is_valid_logdensity(model, logdensity[i])
            particles[i] = op(float, unconditional_sample(rng, model))
            logdensity[i] = loglike(model, push_p(model, particles[i]))
            retrys -= 1
            retrys < 0 &&
                error("Prior leads to ∞ costs too often, tune the prior or increase `retry_sampling`.")
        end
    end

    push_p(model, particles[end]), AISState(particles, logdensity)
end

function AbstractMCMC.step(
    rng::Random.AbstractRNG,
    model::AbstractMCMC.AbstractModel,
    spl::AIS,
    state::AISState;
    ntransitions::Int = 1,
    kwargs...,
)
    i = state.i
    for reps = 1:ntransitions
        transition!(model, state.sample, state.loglikelihood, i, rng)
    end
    push_p(model, state.sample[i]),
    AISState(state.sample, state.loglikelihood, 1 + (i % spl.nparticles))
end

function AbstractMCMC.bundle_samples(
    samples::Vector{T},
    ::AbstractMCMC.AbstractModel,
    ::AIS,
    ::Any,
    ::Type;
    kwargs...,
) where {T<:Particle}
    l = length(samples[1].x)
    P = map(x -> Particles(x), getindex.(getfield.(samples, :x), i) for i = 1:l)
    length(P) == 1 && return P[1]
    return P
end

function AbstractMCMC.chainsstack(c::Vector{Vector{T}}) where {T<:Particles}
    nc = length(c)
    pl = length(c[1])
    return [Particles(reduce(vcat, c[n][i].particles for n = 1:nc)) for i = 1:pl]
end

function AbstractMCMC.chainsstack(c::Vector{T}) where {T<:Particles}
    return Particles(reduce(vcat, c[i].particles for i in eachindex(c)))
end

"""
    sample(model, AIS(N), Ns[; optional args])
    sample(model, AIS(N), MCMCThreads(), Ns, Nc[; optional keyword args])
    sample(model, AIS(N), MCMCDistributed(), Ns, Nc[; optional keyword args])

# Generalities

This function will run an Affine Invariant MCMC sampler, and will return a `Particles` object for each parameter,
the mandatory parameters are:

`model`: a subtype of `AbstractDensity`, look at `ApproxPosterior`, `ApproxKernelizedPosterior`, `CommonLogDensity`.

`N`: number of particles in the ensemble, this particles will be evolved to generate new samples.

`Ns`: total number of samples which must be recorded.

`Nc`: total number of chains to run in parallel if MCMCThreads or MCMCDistributed is enabled.

the optional arguments available are:

`discard_initial`: number of mcmc particles to discard before saving any sample.

`ntransitions`: number of mcmc steps per particle between each sample.

`retry_sampling`: number of maximum attempts to resample an initial particle whose cost (or log-density) is ±∞ or NaN.

`progress`: a boolean to disable verbosity

# Minimal Example for `CommonLogDensity`:

```julia
using KissABC
D = CommonLogDensity(
    2, #number of parameters
    rng -> randn(rng, 2), # initial sampling strategy
    x -> -100 * (x[1] - x[2]^2)^2 - (x[2] - 1)^2, # rosenbrock banana log-density
)
res = sample(D, AIS(50), 1000, ntransitions = 100, discard_initial = 500, progress = false)
println(res)
```

output:
```
Particles{Float64,1000}[1.43 ± 1.4, 0.99 ± 0.67]
```

# Minimal Example for `ApproxKernelizedPosterior` (`ApproxPosterior`)

```julia
using KissABC
prior = Uniform(-10, 10) # prior distribution for parameter
sim(μ) = μ + rand((randn() * 0.1, randn())) # simulator function
cost(x) = abs(sim(x) - 0.0) # cost function to compare simulations to target data, in this case simply '0'
plan = ApproxPosterior(prior, cost, 0.01) # Approximate model of log-posterior density (ABC)
#                                           ApproxKernelizedPosterior can be used in the same fashion
res = sample(plan, AIS(100), 2000, discard_initial = 10000, progress = false)
println(res)
```

output:
```
0.0 ± 0.46
```

"""
sample

export sample, AIS, MCMCThreads, MCMCDistributed
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

@recipe function f(grid::CartesianGrid{Dim,T}) where {Dim,T}
  or = coordinates(minimum(grid))
  sp = spacing(grid)
  sz = size(grid)

  seriescolor --> :black
  legend --> false

  if Dim == 1
    @series begin
      seriestype --> :scatterpath
      marker --> :vline
      x = range(or[1], step=sp[1], length=sz[1]+1)
      x, fill(zero(T), sz[1]+1)
    end
  elseif Dim == 2
    @series begin
      seriestype --> :path
      aspect_ratio --> :equal
      xs = range(or[1], step=sp[1], length=sz[1]+1)
      ys = range(or[2], step=sp[2], length=sz[2]+1)
      coords = []
      for x in xs
        push!(coords, (x, first(ys)))
        push!(coords, (x, last(ys)))
        push!(coords, (NaN, NaN))
      end
      for y in ys
        push!(coords, (first(xs), y))
        push!(coords, (last(xs), y))
        push!(coords, (NaN, NaN))
      end
      x = getindex.(coords, 1)
      y = getindex.(coords, 2)
      x, y
    end
  elseif Dim == 3
    @series begin
      seriestype --> :path
      aspect_ratio --> :equal
      xs = range(or[1], step=sp[1], length=sz[1]+1)
      ys = range(or[2], step=sp[2], length=sz[2]+1)
      zs = range(or[3], step=sp[3], length=sz[3]+1)
      coords = []
      for y in ys, z in zs
        push!(coords, (first(xs), y, z))
        push!(coords, (last(xs), y, z))
        push!(coords, (NaN, NaN, NaN))
      end
      for x in xs, z in zs
        push!(coords, (x, first(ys), z))
        push!(coords, (x, last(ys), z))
        push!(coords, (NaN, NaN, NaN))
      end
      for x in xs, y in ys
        push!(coords, (x, y, first(zs)))
        push!(coords, (x, y, last(zs)))
        push!(coords, (NaN, NaN, NaN))
      end
      x = getindex.(coords, 1)
      y = getindex.(coords, 2)
      z = getindex.(coords, 3)
      x, y, z
    end
  else
    @error "cannot plot in more than 3 dimensions"
  end
end

@recipe function f(grid::CartesianGrid{Dim,T}, data::AbstractVector) where {Dim,T}
  or = coordinates(minimum(grid))
  sp = spacing(grid)
  sz = size(grid)
  Z  = reshape(data, sz)

  if Dim == 1
    seriestype --> :path
    x = range(or[1]+sp[1]/2, step=sp[1], length=sz[1])
    x, Z
  elseif Dim == 2
    seriestype --> :heatmap
    aspect_ratio --> :equal
    colorbar --> true
    x = range(or[1]+sp[1]/2, step=sp[1], length=sz[1])
    y = range(or[2]+sp[2]/2, step=sp[2], length=sz[2])
    x, y, reverse(rotr90(Z), dims=2)
  elseif Dim == 3
    seriestype --> :heatmap
    aspect_ratio --> :equal
    colorbar --> true
    x = range(or[1]+sp[1]/2, step=sp[1], length=sz[1])
    y = range(or[2]+sp[2]/2, step=sp[2], length=sz[2])
    z = range(or[3]+sp[3]/2, step=sp[3], length=sz[3])
    Z
  else
    @error "cannot plot in more than 3 dimensions"
  end
end


text:
# --- Call FFTW3 directly --- #

using FFTW

function c_G_to_R( Ns::Array{Int,1}, fG::Array{ComplexF64,1} )
    out = zeros(ComplexF64,Ns[1]*Ns[2]*Ns[3])
    ccall( (:fftw_inv_fft3d, "../libs/fft3d.so"), Nothing,
           (Ptr{ComplexF64}, Ptr{ComplexF64}, Int64, Int64, Int64),
            fG, out, Ns[3],Ns[2],Ns[1] )
    return out
end

"""
Multicolumn version
"""
function c_G_to_R( Ns::Array{Int,1}, fG::Array{ComplexF64,2} )
    Npoints = size(fG)[1]
    Ncol = size(fG)[2]
    out = zeros(ComplexF64,Npoints,Ncol)
    for ic = 1:Ncol
       ccall( (:fftw_inv_fft3d, "../libs/fft3d.so"), Nothing,
              (Ptr{ComplexF64}, Ptr{ComplexF64}, Int64, Int64, Int64),
               fG[:,ic], out[:,ic], Ns[3],Ns[2],Ns[1] )
    end
    return out
end

# --- Call FFTW3 directly
function c_R_to_G( Ns::Array{Int,1}, fR::Array{ComplexF64,1} )
    out = zeros(ComplexF64,Ns[1]*Ns[2]*Ns[3])
    ccall( (:fftw_fw_fft3d, "../libs/fft3d.so"), Nothing,
           (Ptr{ComplexF64}, Ptr{ComplexF64}, Int64, Int64, Int64),
            fR, out, Ns[3],Ns[2],Ns[1] )
    return out
end

function c_R_to_G( Ns::Array{Int,1}, fR_::Array{Float64,1} )
    out = zeros(ComplexF64,Ns[1]*Ns[2]*Ns[3])
    fR = convert(Array{ComplexF64,1},fR_)
    ccall( (:fftw_fw_fft3d, "../libs/fft3d.so"), Nothing,
           (Ptr{ComplexF64}, Ptr{ComplexF64}, Int64, Int64, Int64),
            fR, out, Ns[3],Ns[2],Ns[1] )
    return out
end

"""
Multicolumn version
"""
function c_R_to_G( Ns::Array{Int,1}, fR::Array{ComplexF64,2} )
    Npoints = size(fR)[1]
    Ncol = size(fR)[2]
    out = zeros(ComplexF64,Npoints,Ncol)
    for ic = 1:Ncol
        ccall( (:fftw_fw_fft3d, "../libs/fft3d.so"), Nothing,
               (Ptr{ComplexF64}, Ptr{ComplexF64}, Int64, Int64, Int64),
                fR[:,ic], out[:,ic], Ns[3],Ns[2],Ns[1] )
    end
    return out
end

function G_to_R( Ns::Array{Int,1}, fG::Array{ComplexF64,1} )
    out = reshape( ifft( reshape(fG,Ns[1],Ns[2],Ns[3]) ),size(fG) )
end

function G_to_R( Ns::Array{Int,1}, fG::Array{ComplexF64,2} )
    Npoints = prod(Ns)
    out = zeros( ComplexF64, size(fG) ) # Is this safe?
    for ic = 1:size(fG,2)
        out[:,ic] = reshape( ifft( reshape(fG[:,ic],Ns[1],Ns[2],Ns[3]) ), Npoints )
    end
    return out
end

function R_to_G( Ns::Array{Int,1}, fR::Array{ComplexF64,1} )
    out = reshape( fft( reshape(fR,Ns[1],Ns[2],Ns[3]) ), size(fR) )
end

# In case we forget to convert the input, we convert it in this version
function R_to_G( Ns::Array{Int,1}, fR_::Array{Float64,1} )
    fR = convert(Array{ComplexF64,1},fR_)
    out = reshape( fft( reshape(fR,Ns[1],Ns[2],Ns[3]) ), size(fR) )
end

function R_to_G( Ns::Array{Int,1}, fR::Array{ComplexF64,2} )
    Npoints = prod(Ns)
    Ncol = size(fR,2)
    out = zeros( ComplexF64, size(fR) )
    for ic = 1:Ncol
        out[:,ic] = reshape( fft( reshape(fR[:,ic],Ns[1],Ns[2],Ns[3]) ), Npoints )
    end
    return out
end


text:
module PrettyTables

using Formatting
using Reexport
using Tables
using Markdown

@reexport using Crayons

import Base: Dict, ismalformed, isoverlong, @kwdef

# The performance of PrettyTables.jl does not increase by a lot of optimizations
# that is performed by the compiler. Hence, we disable then to improve compile
# time.
if isdefined(Base, :Experimental) && isdefined(Base.Experimental, Symbol("@optlevel"))
       @eval Base.Experimental.@optlevel 1
end

################################################################################
#                                    Types
################################################################################

include("types.jl")
include("backends/text/types.jl")
include("backends/html/types.jl")
include("backends/latex/types.jl")

################################################################################
#                                  Includes
################################################################################

include("configuration.jl")
include("deprecations.jl")
include("files.jl")
include("helpers.jl")
include("misc.jl")
include("predefined_formatters.jl")
include("print.jl")
include("private.jl")
include("tables.jl")

# Backends
# ========

# Text backend
include("backends/text/alignment.jl")
include("backends/text/cell_parse.jl")
include("backends/text/display.jl")
include("backends/text/fill.jl")
include("backends/text/misc.jl")
include("backends/text/predefined_formats.jl")
include("backends/text/predefined_highlighters.jl")
include("backends/text/print.jl")
include("backends/text/print_aux.jl")
include("backends/text/recipes.jl")
include("backends/text/string.jl")
include("backends/text/title.jl")

include("backends/text/custom_text_cells/custom_text_cell.jl")
include("backends/text/custom_text_cells/url_text_cell.jl")

# HTML backend
include("backends/html/cell_parse.jl")
include("backends/html/predefined_formats.jl")
include("backends/html/predefined_highlighters.jl")
include("backends/html/print.jl")
include("backends/html/private.jl")

# LaTeX backend
include("backends/latex/cell_parse.jl")
include("backends/latex/predefined_formats.jl")
include("backends/latex/private.jl")
include("backends/latex/print.jl")
include("backends/latex/string.jl")

if Base.VERSION >= v"1.4.2"
    # This try/catch is necessary in case the precompilation statements do not
    # exists. In this case, PrettyTables.jl will work correctly but without the
    # optimizations.
    try
        include("precompile/precompile_PrettyTables.jl")
        _precompile_()
    catch
    end
end

end # module


text:
S = SeisData(3)
S.fs = [1.0, 2.0, 3.0]
S.x[1] = randn(1024); S.t[1] = [1 0; 2 1092831; 1000 12; 1023 12312421; 1024 0]
S.x[2] = randn(768); S.t[2] = [1 0; 768 0]
S.x[3] = randn(1900); S.t[3] = [1 0; 1900 0]
@test nx_max(S) == 1900


text:
import DataFrames
import PlotlyJS
import Statistics

"""
```julia
chill(::Val{:CumulativeChillPortion}, weatherData::DataFrame)::DataFrame
```

Implements the [Cumulative Chill Portion](chill_CumulativeChillPortion.md) model
(also known as the Dynamic model).

# Arguments

- `weatherData::DataFrames.DataFrame`: The weather data used for calculating
    the accumulated chill portions. `weatherData` must be formatted to match the
    output of [`BiTWeather.read`](@ref) and contain at least the`:temperature`
    field.

# Return

- A `DataFrames.DataFrame` containing the chill portions indexed by
    year/month/day/hour. The returned columns are:  year, month, day, hour,
    temperature, ``x(n)``, ``\\Delta y(n)``, and ``y(n)``.

# Example

```julia
import DataFrames
import Dates
import BiTWeather

configuration = BiTWeather.Configuration(
    dsn = "meteobridge",
    table = "backinthirty",
    fieldMappings = Dict(
        :dateTime =>    BiTWeather.FieldMapping(name = "DateTime",   units = nothing),
        :temperature => BiTWeather.FieldMapping(name = "TempOutNow", units = "F")
    )
)

year::Int = 2021
rangeStart::Dates.Date = Dates.Date(year - 1, 9, 1)
rangeEnd::Dates.Date = Dates.Date(year, 2, 28)
range::Tuple{Float64, Float64} = (rangeStart, rangeEnd)
fields::Vector{Symbol} = [:temperature]
weatherData::DataFrames.DataFrame = BiTWeather.read(configuration, range, fields)
chillData::DataFrames.DataFrame = BiTWeather.chill(Val(:CumulativeChillPortion), weatherData)
println(chillData)
```
"""
function chill(::Val{:CumulativeChillPortion}, weatherData::DataFrames.DataFrame)::DataFrames.DataFrame

    # The chill calculation requires the :temperature field.
    if !(:temperature in DataFrames.propertynames(weatherData))
        error("weatherData must contain the :temperature field")
    end

    #= -------------------------------------------------------------------------
    This function performs the chill calculation for one interval. It is passed
    the temperature as well as the state needed to perform the calculation. It
    returns the result of the calculation and the state needed to perform the
    next interval's calculation.

    Arguments:
        temperature: The temperature.
        x_0_n::Float64:
        y_nMinus1::Float64: The accumulated chill portions after the previous
            interval's calculation
    Returns:
        Result:
            delta_y_n::Float64: The number of chill portions from this interval.
            y_n::Float64: The accumulated chill portions at the completion of
                this interval.
        State:
            x_0_n
            y_nMinus1
    ------------------------------------------------------------------------- =#
    function calculate(temperature::Float64, x_0_n::Float64, y_nMinus1::Float64)

        e_0::Float64     = 4153.5
        e_1::Float64     = 12888.8
        a_0::Float64     = 1.395 * 10.0^5.0
        a_1::Float64     = 2.567 * 10.0^18.0
        f::Float64       = 0.4
        theta_f::Float64 = 4.0 + 273.0

        theta_n::Float64 = temperature + 273

        k_0_n::Float64       = a_0 * exp(-e_0 / theta_n)
        k_1_n::Float64       = a_1 * exp(-e_1 / theta_n)
        x_s_n::Float64       = k_0_n / k_1_n
        x_n::Float64         = x_s_n - (x_s_n - x_0_n) * exp(-k_1_n)
        p_n::Float64         = exp(4 * f * theta_f * (1 - (theta_f / theta_n))) /
                               (1 + exp(4 * f * theta_f * (1 - (theta_f / theta_n))))
        delta_y_n::Float64   = x_n < 1.0 ? 0.0 : p_n * x_n
        y_n::Float64         = y_nMinus1 + delta_y_n
        x_0_nPlus1::Float64  = x_n - delta_y_n

        return x_n, delta_y_n, y_n, x_0_nPlus1, y_n
    end

    ccp::DataFrames.DataFrame = weatherData

    # Determine the hourly temperature.
    gdf::DataFrames.GroupedDataFrame = DataFrames.groupby(ccp, [:year, :month, :day, :hour])
    ccp = DataFrames.combine(gdf,
        [:temperature] .=> Statistics.mean .=> [:temperature]
    )
    DataFrames.select!(ccp, [:year, :month, :day, :hour, :temperature])
    DataFrames.sort!(ccp, [:year, :month, :day, :hour])

    x::Vector{Float64} = []
    delta_y::Vector{Float64} = []
    y::Vector{Float64} = []

    x_0_n::Float64 , y_nMinus1::Float64 = 0.0, 0.0
    for row in DataFrames.eachrow(ccp)
        x_n::Float64, delta_y_n::Float64, y_n::Float64, x_0_n, y_nMinus1 = calculate(row.temperature, x_0_n, y_nMinus1)
        push!(x, x_n)
        push!(delta_y, delta_y_n)
        push!(y, y_n)
    end

    ccp[!, :x_n] = x
    ccp[!, :delta_y_n] = delta_y
    ccp[!, :y_n] = y

    for name in DataFrames.propertynames(ccp)
        if eltype(ccp[!, name]) == Float64
            ccp[!, name] = round.(ccp[!, name], digits = 2)
        end
    end

    return ccp
end

"""
```julia
chillPlot(::Val{:CumulativeChillPortion}, chillData::DataFrame)::PlotlyJS.Plot
```

Generates a plot of the [Cumulative Chill Portion](chill_CumulativeChillPortion.md)
data.

# Arguments

- `chillData::DataFrames.DataFrame`: The chill data used in generating the plot.
    `chillData` must be formated to match the output of
    [`BiTWeather.chill`](@ref) called with `Val(:CumulativeChillPortion)`.

# Return

- A `PlotlyJS.Plot` containing a plot of the temperature and chill portion data.

# Example

```julia
import DataFrames
import Dates
import PlotlyJS
import BiTWeather

configuration = BiTWeather.Configuration(
    dsn = "meteobridge",
    table = "backinthirty",
    fieldMappings = Dict(
        :dateTime =>    BiTWeather.FieldMapping(name = "DateTime",   units = nothing),
        :temperature => BiTWeather.FieldMapping(name = "TempOutNow", units = "F")
    )
)

year::Int = 2021
rangeStart::Dates.Date = Dates.Date(year - 1, 9, 1)
rangeEnd::Dates.Date = Dates.Date(year, 2, 28)
range::Tuple{Float64, Float64} = (rangeStart, rangeEnd)
fields::Vector{Symbol} = [:temperature]
weatherData::DataFrames.DataFrame = BiTWeather.read(configuration, range, fields)
chillData::DataFrames.DataFrame = BiTWeather.chill(Val(:CumulativeChillPortion), weatherData)
plot::PlotlyJS.Plot = BiTWeather.chillPlot(Val(:CumulativeChillPortion), chillData)
PlotlyJS.display(plot)
```
"""
function chillPlot(::Val{:CumulativeChillPortion}, chillData::DataFrames.DataFrame)::PlotlyJS.Plot
    return _chillPlot(chillData, :temperature, :y_n, "Cumulative Chill Portion (Dynamic)", "Portions")
end

max_stars_repo_name: pebender/BiTWeather

max_issues_repo_name: pebender/BiTWeather

max_forks_repo_name: pebender/BiTWeather


text:
# This file includes modified code from https://github.com/odow/SDDP.jl as at
# commit a6662eecff6d5de5c499b760773a8569a0979b9b.

#  Copyright 2017, Eyob Zewdie, Oscar Dowson
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################

#=
    This example concerns the acceptance of booking requests for rooms in a
    hotel in the leadup to a large event.

    Each stage, we receive a booking request and can choose to accept or decline
    it. Once accepted, bookings cannot be terminated.
=#
using SDDP, JuMP, GLPKMathProgInterface, Base.Test, SDDiP

function bookingmanagementmodel(NUMBER_OF_DAYS, NUMBER_OF_ROOMS, NUMBER_OF_REQUESTS)

    # maximum revenue that could be accrued.
    MAX_REVENUE = (NUMBER_OF_ROOMS + NUMBER_OF_REQUESTS) * NUMBER_OF_DAYS * NUMBER_OF_ROOMS

    #=
        BOOKING_REQUESTS is a vector of {0,1} arrays of size
        (NUMBER_OF_DAYS x NUMBER_OF_ROOMS) if the room is requested.
    =#
    BOOKING_REQUESTS = Array{Int, 2}[]
    for room in 1:NUMBER_OF_ROOMS
        for day in 1:NUMBER_OF_DAYS
            # note: length_of_stay is 0 indexed to avoid unncecessary +/- 1
            # on the indexing
            for length_of_stay in 0:(NUMBER_OF_DAYS - day)
                booking_request = zeros(Int, (NUMBER_OF_ROOMS, NUMBER_OF_DAYS))
                booking_request[room:room, day + (0:length_of_stay)] = 1
                push!(BOOKING_REQUESTS, booking_request)
            end
        end
    end

    m = SDDPModel(
                 stages = NUMBER_OF_REQUESTS,
        objective_bound = MAX_REVENUE,
                  sense = :Max,
                 solver = GLPKSolverMIP()
                            ) do sp, stage

        # Seat remaining?
        @binarystate(sp, 0 <= vacancy_after_decision[room=1:NUMBER_OF_ROOMS, day=1:NUMBER_OF_DAYS] <= 1, vacancy_before_decision==1)
        @variables(sp, begin
            # accept request for booking of room for length of time
            0 <= accept_request <= 1
            # accept a booking for an individual room on an individual day
            0 <= room_request_accepted[1:NUMBER_OF_ROOMS, 1:NUMBER_OF_DAYS] <= 1
        end)
        for room in 1:NUMBER_OF_ROOMS, day in 1:NUMBER_OF_DAYS
            @constraints(sp, begin
                # Update vacancy if we accept a room request
                vacancy_after_decision[room, day] == vacancy_before_decision[room, day] - room_request_accepted[room, day]
                # can't accept a request of a filled room
                room_request_accepted[room, day] <= vacancy_before_decision[room, day]
                # can't accept invididual room request if entire request is declined
                room_request_accepted[room, day] <= accept_request

            end)
            @rhsnoises(sp, booking_request = BOOKING_REQUESTS, begin
                # can't accept request if room not requested
                room_request_accepted[room, day] <= booking_request[room, day]
                # accept all individual rooms is entire request is accepted
                room_request_accepted[room, day] + (1-accept_request) >= booking_request[room, day]
            end)
        end
        @stageobjective(sp, sum(
            (room + stage - 1) * room_request_accepted[room, day]
            for room in 1:NUMBER_OF_ROOMS for day in 1:NUMBER_OF_DAYS
            )
        )
        setSDDiPsolver!(sp, method=KelleyMethod())
    end
end

srand(1234)
m_1_2_5 = bookingmanagementmodel(1, 2, 5)
@test solve(m_1_2_5, max_iterations = 10) == :iteration_limit
@test isapprox(getbound(m_1_2_5), 7.25, atol=0.001)

srand(1234)
m_2_2_3 = bookingmanagementmodel(2, 2, 3)
@test solve(m_2_2_3, max_iterations = 90) == :iteration_limit
@test isapprox(getbound(m_2_2_3), 6.13, atol=0.001)


text:
using MoodleQuiz
using Documenter

makedocs(;
    modules=[MoodleQuiz],
    authors="Michal Sojka <michal.sojka@cvut.cz> and contributors",
    repo="https://github.com/wentasah/MoodleQuiz.jl/blob/{commit}{path}#L{line}",
    sitename="MoodleQuiz.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wentasah.github.io/MoodleQuiz.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/wentasah/MoodleQuiz.jl",
)


text:
using Test
import CryptoGroups.Specs: gn_basis_exist, gn_basis_representation_rule

@test gn_basis_exist(4, 3) == true

# C.1 Table of GNB in X9.62
@test gn_basis_representation_rule(161) == 6
@test gn_basis_representation_rule(185) == 8
@test gn_basis_representation_rule(186) == 2
@test gn_basis_representation_rule(190) == 10
@test gn_basis_representation_rule(300) == 19
@test gn_basis_representation_rule(487) == 4
@test gn_basis_representation_rule(628) == 7
@test gn_basis_representation_rule(1380) == 1
@test gn_basis_representation_rule(1703) == 2