text:
@doc raw"""
Cache of a Specialised Partitioned Additive Runge-Kutta integrator.

### Fields

* `n`: time step number
* `t`: time of current time step
* `t̄`: time of previous time step
* `q`: current solution of q
* `q̄`: previous solution of q
* `p`: current solution of p
* `p̄`: previous solution of p
* `v`: vector field of q
* `v̄`: vector field of q̄
* `f`: vector field of p
* `f̄`: vector field of p̄
* `q̃`: initial guess of q
* `p̃`: initial guess of p
* `ṽ`: initial guess of v
* `f̃`: initial guess of f
* `s̃`: holds shift due to periodicity of solution
* `Q`: internal stages of q
* `P`: internal stages of p
* `V`: internal stages of v
* `F`: internal stages of f
* `Y`: vector field of internal stages of q
* `Z`: vector field of internal stages of p
"""
mutable struct IntegratorCacheSPARK{ST,D,S,R} <: IDAEIntegratorCache{ST,D}
    q::Vector{ST}
    q̄::Vector{ST}
    p::Vector{ST}
    p̄::Vector{ST}
    λ::Vector{ST}
    λ̄::Vector{ST}
    μ::Vector{ST}
    μ̅::Vector{ST}

    v::Vector{ST}
    v̄::Vector{ST}
    f::Vector{ST}
    f̄::Vector{ST}
    u::Vector{ST}
    ū::Vector{ST}
    g::Vector{ST}
    ḡ::Vector{ST}
    ϕ::Vector{ST}
    ϕ̅::Vector{ST}

    q̃::Vector{ST}
    p̃::Vector{ST}
    ṽ::Vector{ST}
    f̃::Vector{ST}
    ϕ̃::Vector{ST}
    s̃::Vector{ST}

    Qi::Vector{Vector{ST}}
    Pi::Vector{Vector{ST}}
    Vi::Vector{Vector{ST}}
    Fi::Vector{Vector{ST}}
    Gi::Vector{Vector{ST}}
    Hi::Vector{Vector{ST}}
    Yi::Vector{Vector{ST}}
    Zi::Vector{Vector{ST}}
    Φi::Vector{Vector{ST}}
    Ψi::Vector{Vector{ST}}

    Qp::Vector{Vector{ST}}
    Pp::Vector{Vector{ST}}
    Vp::Vector{Vector{ST}}
    Λp::Vector{Vector{ST}}
    Up::Vector{Vector{ST}}
    Fp::Vector{Vector{ST}}
    Gp::Vector{Vector{ST}}
    G̅p::Vector{Vector{ST}}
    Hp::Vector{Vector{ST}}
    Yp::Vector{Vector{ST}}
    Zp::Vector{Vector{ST}}
    Φp::Vector{Vector{ST}}
    Ψp::Vector{Vector{ST}}

    function IntegratorCacheSPARK{ST,D,S,R}() where {ST,D,S,R}
        q = zeros(ST,D)
        q̄ = zeros(ST,D)
        p = zeros(ST,D)
        p̄ = zeros(ST,D)
        λ = zeros(ST,D)
        λ̄= zeros(ST,D)
        μ = zeros(ST,D)
        μ̅ = zeros(ST,D)

        # create update vectors
        v = zeros(ST,D)
        v̄ = zeros(ST,D)
        f = zeros(ST,D)
        f̄ = zeros(ST,D)
        u = zeros(ST,D)
        ū = zeros(ST,D)
        g = zeros(ST,D)
        ḡ = zeros(ST,D)
        ϕ = zeros(ST,D)
        ϕ̅ = zeros(ST,D)

        # create temporary vectors
        q̃ = zeros(ST,D)
        p̃ = zeros(ST,D)
        ṽ = zeros(ST,D)
        f̃ = zeros(ST,D)
        ϕ̃ = zeros(ST,D)
        s̃ = zeros(ST,D)

        # create internal stage vectors
        Qi = create_internal_stage_vector(ST, D, S)
        Pi = create_internal_stage_vector(ST, D, S)
        Vi = create_internal_stage_vector(ST, D, S)
        Fi = create_internal_stage_vector(ST, D, S)
        Gi = create_internal_stage_vector(ST, D, S)
        Hi = create_internal_stage_vector(ST, D, S)
        Yi = create_internal_stage_vector(ST, D, S)
        Zi = create_internal_stage_vector(ST, D, S)
        Φi = create_internal_stage_vector(ST, D, S)
        Ψi = create_internal_stage_vector(ST, D, S)

        Qp = create_internal_stage_vector(ST, D, R)
        Pp = create_internal_stage_vector(ST, D, R)
        Vp = create_internal_stage_vector(ST, D, R)
        Λp = create_internal_stage_vector(ST, D, R)
        Up = create_internal_stage_vector(ST, D, R)
        Fp = create_internal_stage_vector(ST, D, R)
        Gp = create_internal_stage_vector(ST, D, R)
        G̅p = create_internal_stage_vector(ST, D, R)
        Hp = create_internal_stage_vector(ST, D, R)
        Yp = create_internal_stage_vector(ST, D, R)
        Zp = create_internal_stage_vector(ST, D, R)
        Φp = create_internal_stage_vector(ST, D, R)
        Ψp = create_internal_stage_vector(ST, D, R)

        new(q, q̄, p, p̄, λ, λ̄, μ, μ̅,
            v, v̄, f, f̄, u, ū, g, ḡ, ϕ, ϕ̅,
            q̃, p̃, ṽ, f̃, ϕ̃, s̃,
            Qi, Pi, Vi, Fi, Gi, Hi, Yi, Zi, Φi, Ψi,
            Qp, Pp, Vp, Λp, Up, Fp, Gp, G̅p, Hp, Yp, Zp, Φp, Ψp)
    end
end

function IntegratorCache(int::AbstractIntegratorSPARK{DT,TT}) where {DT,TT}
    IntegratorCacheSPARK{DT, ndims(int), nstages(int), pstages(int)}()
end

function GeometricBase.reset!(cache::IntegratorCacheSPARK, Δt)
    cache.t̄  = cache.t
    cache.q̄ .= cache.q
    cache.p̄ .= cache.p
    cache.λ̄.= cache.λ
    cache.v̄ .= cache.v
    cache.f̄ .= cache.f
    cache.ū .= cache.u
    cache.ḡ .= cache.g
    cache.ϕ̅ .= cache.ϕ
end


text:
module War

using DataStructures: isempty, collect, popfirst!
using PlayingCards52, DataStructures

import Base: show

export Position, battle, play_war, verbose_war

ace_high()

Pile = Deque{Card}

"""
    make_pile(clist::Vector(Card))

Make a `Pile` from a list of cards.
"""
function make_pile(clist::Vector{Card})::Pile
    result = Pile()
    add_to_pile(result, clist)
    return result
end

"""
    add_to_pile(pile::Pile, clist::Vector{Card})

Add the cards in `clist` to the bottom of `pile`.
"""
function add_to_pile(pile::Pile, clist::Vector{Card})
    for c in clist
        push!(pile, c)
    end
end

function _string(pile::Pile)::String
    result = "[ "
    for c in pile
        result *= "$c "
    end
    result *= "]"
    return result
end

struct Position
    A::Pile
    B::Pile
    function Position(d::Vector{Card})
        new(make_pile(d[1:26]), make_pile(d[27:52]))
    end
    function Position(P1::Pile, P2::Pile)
        new(P1, P2)
    end
end

Position() = Position(deck())

function Position(d1::Vector{Card}, d2::Vector{Card})
    nd = length(d1) + length(d2)
    S = Set{Card}(d1) ∪ Set{Card}(d2)
    ns = length(S)
    if nd != 52 || ns != 52
        @warn "Improper deck: the there are $nd cards in the two piles and $ns cards in their union. Carrying on anyhow."
    end
    A = make_pile(d1)
    B = make_pile(d2)
    Position(A, B)
end

function show(io::IO, P::Position)
    println(io, _string(P.A))
    println(io, _string(P.B))
end

const A_win = -1
const B_win = 1
const Tie = 0

VERBOSE = false

"""
    verbose_war
Use `verbose_war(true)` to have detailed reporting during a game of war and `verbose_war(false)` to supress that.
"""
function verbose_war(tf::Bool = true)
    global VERBOSE = tf
end

"""
    battle_compare(a::Card, b::Card)::Int

Compare two cards by War rules. Returns:
* `A_win` if `a` beats `b`.
* `B_win` if `b` beats `a`.
* `Tie` if `a` and `b` have the same rank.
"""
function battle_compare(a::Card, b::Card)::Int
    if rank(a) == rank(b)
        return Tie
    end
    return a > b ? A_win : B_win
end

function give_spoils(P::Position, used::Set{Card}, winner::Int)
    spoils = shuffle(collect(used))
    if winner == A_win
        add_to_pile(P.A, spoils)
    end
    if winner == B_win
        add_to_pile(P.B, spoils)
    end
    nothing
end

"""
    battle(P::Position, used::Set{Card})

execute one turn in War.
"""
function battle(P::Position, used::Set{Card} = Set{Card}())

    # if a player is out of cards, the spoils go the other player (who then wins)
    if isempty(P.A) || isempty(P.B)
        result = isempty(P.A) ? B_win : A_win
        give_spoils(P, used, result)
        return
    end

    # get the first card from each player and compare
    a = popfirst!(P.A)
    b = popfirst!(P.B)
    push!(used, a)
    push!(used, b)

    result = battle_compare(a, b)
    if result != Tie
        give_spoils(P, used, result)
        return
    end

    if VERBOSE
        println("It's war!\n")
    end

    # need at least 4 cards or the game is lost

    if length(P.A) < 4      # A loses
        while !isempty(P.A)
            c = popfirst!(P.A)
            push!(used, c)
        end
        give_spoils(P, used, B_win)
        return
    end

    if length(P.B) < 4      # B loses
        while !isempty(P.B)
            c = popfirst!(P.B)
            push!(used, c)
        end
        give_spoils(P, used, A_win)
        return
    end

    # both have enough cards
    for k = 1:3
        c = popfirst!(P.A)
        push!(used, c)
        c = popfirst!(P.B)
        push!(used, c)
    end
    battle(P, used)
end

function play_war(P::Position)
    result = [length(P.A)]
    if VERBOSE
        println(P)
    end
    while length(P.A) > 0 && length(P.B) > 0
        battle(P)
        if VERBOSE
            println(P)
        end
        push!(result, length(P.A))
    end
    return result
end

function play_war(d1::Vector{Card}, d2::Vector{Card})
    P = Position(d1, d2)
    play_war(P)
end

function play_war(d::Vector{Card})
    P = Position(d)
    play_war(P)
end

"""
    play_war
Play a game of War. Options:
+ `play_war()`: Play a new game of war with random deck of cards.
+ `play_war(d::Vector{Card})`: Play war with a given deck of cards.
+ `play_war(d1::Vector{Card},d2::Vector{Card})`: Play war in which player 1 gets the cards in `d1` and player 2 gets the cards in `d2`.

The output of `play_war` is vector containing the number of cards held by the first player at each turn. That vector necessarily begins with the number 26 and ends with either 0 or 52. The length of the vector is one more than the number of battles to complete the game. In a normal battle, the difference between consecutive terms is ±1; after an ordinary war, the difference is ±5; after a double war, the difference is ±9, and so forth.

See: `verbose_war`.
"""
play_war() = play_war(deck())

end  #end of module War

ext: jl


text:
include("window.jl")
include("eventloop.jl")


text:
line = readline()
x, y, z = [parse(Int, x) for x in split(line)]
println("$z $x $y")

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using Test
using ManifoldProjections
using LinearAlgebra
using Random
Random.seed!(0)

# test sphere with radius 1
x = randn(4)+im*randn(4)
M = Sphere()
retract!(M,x)
@test norm(x) ≈ 1
v = randn(4)+im*randn(4)
@test project_tangent(M,v,x) == project_tangent!(M,v,x)
@test abs(real(v'*x)) < 1e-8
ε = 1e-8
@test (norm(x + ε*v) - 1) < ε^(3/2) # the tangent vector is such that x+εv is on the manifold up to O(ε^2)

# test sphere with radius r
r = 120.0
x = randn(4)+im*randn(4)
M = Sphere(r)
retract!(M,x)
@test norm(x) ≈ r
v = randn(4)+im*randn(4)
@test project_tangent(M,v,x) == project_tangent!(M,v,x)
@test abs(real(v'*x)) < 1e-8
@test (norm(x + ε*v) - r) < ε^(3/2) # the tangent vector is such that x+εv is on the manifold up to O(ε^2)

for M in (Stiefel_CholQR(), Stiefel_SVD())
    X = randn(4,2) + im*randn(4,2)
    retract!(M,X)
    @test X'X ≈ I
    V = randn(4,2)+im*randn(4,2)
    @test project_tangent(M,V,X) == project_tangent!(M,V,X)
    @test norm((X + ε*V)'*(X+ε*V) - I) < ε^(3/2)
end


text:
mutable struct UpdateTrigger
    execution_wait_count::Int
    current_count::Int
end

function set_execution_wait_count!(trigger::UpdateTrigger, val::Int)
    trigger.execution_wait_count = val
    return
end

function update_count!(trigger::UpdateTrigger)
    trigger.current_count += 1
    return
end

function trigger_update(trigger::UpdateTrigger)
    return trigger.current_count == trigger.execution_wait_count
end

function reset_trigger_count!(trigger::UpdateTrigger)
    trigger.current_count = 0
    return
end

function initialize_trigger_count!(trigger::UpdateTrigger)
    trigger.current_count = trigger.execution_wait_count
    return
end

function get_execution_wait_count(trigger::UpdateTrigger)
    return trigger.execution_wait_count
end

############################Chronologies For FeedForward###################################
@doc raw"""
    Synchronize(periods::Int)
Defines the co-ordination of time between Two stages.

# Arguments
- `periods::Int`: Number of time periods to grab data from
"""
mutable struct Synchronize <: FeedForwardChronology
    periods::Int
    current::Int
    trigger::UpdateTrigger
    function Synchronize(; periods)
        new(periods, 0, UpdateTrigger(-1, -1))
    end
end
# TODO: Add DocString
"""
    RecedingHorizon(period::Int)
"""
mutable struct RecedingHorizon <: FeedForwardChronology
    periods::Int
    trigger::UpdateTrigger
    function RecedingHorizon(; periods::Int = 1)
        new(periods, UpdateTrigger(-1, -1))
    end
end

mutable struct Consecutive <: FeedForwardChronology
    trigger::UpdateTrigger
    function Consecutive()
        new(UpdateTrigger(-1, -1))
    end
end

mutable struct FullHorizon <: FeedForwardChronology
    trigger::UpdateTrigger
    function FullHorizon()
        new(UpdateTrigger(-1, -1))
    end
end

mutable struct Range <: FeedForwardChronology
    range::UnitRange{Int}
    trigger::UpdateTrigger
    function Range(; range::UnitRange{Int})
        new(range, UpdateTrigger(-1, -1))
    end
end

function check_chronology!(sim::Simulation, key::Pair, sync::Synchronize)
    source_stage = get_stage(sim, key.first)
    destination_stage = get_stage(sim, key.second)
    source_stage_horizon = sim.sequence.horizons[key.first]
    destination_stage_horizon = sim.sequence.horizons[key.second]
    source_stage_interval = get_stage_interval(sim, key.first)
    destination_stage_interval = get_stage_interval(sim, key.second)

    source_stage_resolution = get_resolution(source_stage)
    @debug source_stage_resolution, destination_stage_interval
    # How many times the second stages executes per solution retireved from the source_stage.
    # E.g. source_stage_resolution = 1 Hour, destination_stage_interval = 5 minutes => 12 executions per solution
    destination_stage_executions_per_solution =
        Int(source_stage_resolution / destination_stage_interval)
    # Number of periods in the horizon that will be synchronized between the source_stage and the destination_stage
    source_stage_sync = sync.periods

    if source_stage_sync > source_stage_horizon
        throw(IS.ConflictingInputsError("The lookahead length $(source_stage_horizon) in stage is insufficient to syncronize with $(source_stage_sync) feedforward periods"))
    end

    if (source_stage_sync % destination_stage_executions_per_solution) != 0
        throw(IS.ConflictingInputsError("The current configuration implies $(source_stage_sync / destination_stage_executions_per_solution) executions of $(key.second) per execution of $(key.first). The number of Synchronize periods $(sync.periods) in stage $(key.first) needs to be a mutiple of the number of stage $(key.second) execution for every stage $(key.first) interval."))
    end

    return
end

function check_chronology!(sim::Simulation, key::Pair, ::Consecutive)
    source_stage_horizon = sim.sequence.horizons[key.first]
    source_stage_interval = get_stage_interval(sim, key.first)
    if source_stage_horizon != source_stage_interval
        @warn("Consecutive Chronology Requires the same interval and horizon, the parameter horizon = $(source_stage_horizon) in stage $(key.first) will be replaced with $(source_stage_interval). If this is not the desired behviour consider changing your chronology to RecedingHorizon")
    end
    sim.sequence.horizons[key.first] = get_stage_interval(sim, key.first)
    return
end

check_chronology!(sim::Simulation, key::Pair, ::RecedingHorizon) = nothing
check_chronology!(sim::Simulation, key::Pair, ::FullHorizon) = nothing
# TODO: Add missing check
check_chronology!(sim::Simulation, key::Pair, ::Range) = nothing

function check_chronology!(
    sim::Simulation,
    key::Pair,
    ::T,
) where {T <: FeedForwardChronology}
    error("Chronology $(T) not implemented")
    return
end

############################FeedForward Definitions########################################

struct UpperBoundFF <: AbstractAffectFeedForward
    variable_source_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function UpperBoundFF(
        variable_source_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_source_stage), Symbol.(affected_variables), cache)
    end
end

function UpperBoundFF(; variable_source_stage, affected_variables)
    return UpperBoundFF(variable_source_stage, affected_variables, nothing)
end

get_variable_source_stage(p::UpperBoundFF) = p.variable_source_stage

struct RangeFF <: AbstractAffectFeedForward
    variable_source_stage_ub::Symbol
    variable_source_stage_lb::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function RangeFF(
        variable_source_stage_ub::AbstractString,
        variable_source_stage_lb::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(
            Symbol(variable_source_stage_ub),
            Symbol(variable_source_stage_lb),
            Symbol.(affected_variables),
            cache,
        )
    end
end

function RangeFF(; variable_source_stage_ub, variable_source_stage_lb, affected_variables)
    return RangeFF(
        variable_source_stage_ub,
        variable_source_stage_lb,
        affected_variables,
        nothing,
    )
end

get_bounds_source_stage(p::RangeFF) =
    (p.variable_source_stage_lb, p.variable_source_stage_lb)

struct SemiContinuousFF <: AbstractAffectFeedForward
    binary_source_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function SemiContinuousFF(
        binary_source_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(binary_source_stage), Symbol.(affected_variables), cache)
    end
end

function SemiContinuousFF(; binary_source_stage, affected_variables)
    return SemiContinuousFF(binary_source_stage, affected_variables, nothing)
end

get_binary_source_stage(p::SemiContinuousFF) = p.binary_source_stage
get_affected_variables(p::AbstractAffectFeedForward) = p.affected_variables

struct IntegralLimitFF <: AbstractAffectFeedForward
    variable_source_stage::Symbol
    affected_variables::Vector{Symbol}
    cache::Union{Nothing, Type{<:AbstractCache}}
    function IntegralLimitFF(
        variable_source_stage::AbstractString,
        affected_variables::Vector{<:AbstractString},
        cache::Union{Nothing, Type{<:AbstractCache}},
    )
        new(Symbol(variable_source_stage), Symbol.(affected_variables), cache)
    end
end

function IntegralLimitFF(; variable_source_stage, affected_variables)
    return IntegralLimitFF(variable_source_stage, affected_variables, nothing)
end

get_variable_source_stage(p::IntegralLimitFF) = p.variable_source_stage

struct ParameterFF <: AbstractAffectFeedForward
    variable_source_stage::Symbol
    affected_parameters::Any
    function ParameterFF(variable_source_stage::AbstractString, affected_parameters)
        new(Symbol(variable_source_stage), affected_parameters)
    end
end

function ParameterFF(; variable_source_stage, affected_parameters)
    return ParameterFF(variable_source_stage, affected_parameters)
end

####################### Feed Forward Affects ###############################################

@doc raw"""
        ub_ff(psi_container::PSIContainer,
              cons_name::Symbol,
              constraint_infos::Vector{DeviceRangeConstraintInfo},
              param_reference::UpdateRef,
              var_name::Symbol)

Constructs a parametrized upper bound constraint to implement feedforward from other models.
The Parameters are initialized using the uppper boundary values of the provided variables.

# Constraints
``` variable[var_name, t] <= param_reference[var_name] ```

# LaTeX

`` x \leq param^{max}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function ub_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    variable = get_variable(psi_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps
    container = add_param_container!(psi_container, param_reference, set_name)
    param_ub = get_parameter_array(container)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        value = JuMP.upper_bound(variable[name, 1])
        param_ub[name] = PJ.add_parameter(psi_container.JuMPmodel, value)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(expression_ub, variable[name, t])
            end
            con_ub[name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, expression_ub <= param_ub[name])
        end
    end
    return
end

@doc raw"""
        range_ff(psi_container::PSIContainer,
                        cons_name::Symbol,
                        param_reference::NTuple{2, UpdateRef},
                        var_name::Symbol)

Constructs min/max range parametrized constraint from device variable to include feedforward.

# Constraints

``` param_reference[1][var_name] <= variable[var_name, t] ```
``` variable[var_name, t] <= param_reference[2][var_name] ```

where r in range_data.

# LaTeX

`` param^{min} \leq x ``
`` x \leq param^{max}``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* param_reference::NTuple{2, UpdateRef} : Tuple with the lower bound and upper bound parameter reference
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
"""
function range_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::NTuple{2, UpdateRef},
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")

    variable = get_variable(psi_container, var_name)
    # Used to make sure the names are consistent between the variable and the infos
    axes = JuMP.axes(variable)
    set_name = axes[1]
    @assert axes[2] == time_steps

    #Create containers for the constraints
    container_lb = add_param_container!(psi_container, param_reference[1], set_name)
    param_lb = get_parameter_array(container_lb)
    container_ub = add_param_container!(psi_container, param_reference[2], set_name)
    param_ub = get_parameter_array(container_ub)

    #Create containers for the parameters
    con_lb = add_cons_container!(psi_container, lb_name, set_name, time_steps)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        param_lb[name] =
            PJ.add_parameter(psi_container.JuMPmodel, JuMP.lower_bound(variable[name, 1]))
        param_ub[name] =
            PJ.add_parameter(psi_container.JuMPmodel, JuMP.upper_bound(variable[name, 1]))
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(expression_ub, variable[name, t])
            end
            expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_lb
                JuMP.add_to_expression!(expression_lb, variable[name, t], -1.0)
            end
            con_ub[name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, expression_ub <= param_ub[name])
            con_lb[name, t] =
                JuMP.@constraint(psi_container.JuMPmodel, expression_lb >= param_lb[name])
        end
    end

    return
end

@doc raw"""
            semicontinuousrange_ff(psi_container::PSIContainer,
                                    cons_name::Symbol,
                                    var_name::Symbol,
                                    param_reference::UpdateRef)

Constructs min/max range constraint from device variable with parameter setting.

# Constraints
If device min = 0:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

Otherwise:

``` variable[var_name, t] <= r[2].max*param_reference[var_name] ```

``` variable[var_name, t] >= r[2].min*param_reference[var_name] ```

where r in range_data.

# LaTeX

`` 0.0 \leq x^{var} \leq r^{max} x^{param}, \text{ for } r^{min} = 0 ``

`` r^{min} x^{param} \leq x^{var} \leq r^{min} x^{param}, \text{ otherwise } ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* var_name::Symbol : the name of the continuous variable
* param_reference::UpdateRef : UpdateRef of the parameter
"""
function semicontinuousrange_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    constraint_infos::Vector{DeviceRangeConstraintInfo},
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "ub")
    lb_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "lb")
    variable = get_variable(psi_container, var_name)
    # Used to make sure the names are consistent between the variable and the infos
    axes = JuMP.axes(variable)
    set_name = [get_component_name(ci) for ci in constraint_infos]
    @assert axes[2] == time_steps
    container = add_param_container!(psi_container, param_reference, set_name)
    multiplier = get_multiplier_array(container)
    param = get_parameter_array(container)
    con_ub = add_cons_container!(psi_container, ub_name, set_name, time_steps)
    con_lb = add_cons_container!(psi_container, lb_name, set_name, time_steps)

    for constraint_info in constraint_infos
        name = get_component_name(constraint_info)
        ub_value = JuMP.upper_bound(variable[name, 1])
        lb_value = JuMP.lower_bound(variable[name, 1])
        @debug "SemiContinuousFF" name ub_value lb_value
        param[name] = PJ.add_parameter(psi_container.JuMPmodel, 1.0)
        for t in time_steps
            expression_ub = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_ub
                JuMP.add_to_expression!(
                    expression_ub,
                    get_variable(psi_container, val)[name, t],
                )
            end
            expression_lb = JuMP.AffExpr(0.0, variable[name, t] => 1.0)
            for val in constraint_info.additional_terms_lb
                JuMP.add_to_expression!(
                    expression_lb,
                    get_variable(psi_container, val)[name, t],
                    -1.0,
                )
            end
            con_ub[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_ub <= ub_value * param[name]
            )
            con_lb[name, t] = JuMP.@constraint(
                psi_container.JuMPmodel,
                expression_lb >= lb_value * param[name]
            )
        end
    end

    # If the variable was a lower bound != 0, not removing the LB can cause infeasibilities
    for v in variable
        if JuMP.has_lower_bound(v)
            @debug "lb reset" v
            JuMP.set_lower_bound(v, 0.0)
        end
    end

    return
end

@doc raw"""
        integral_limit_ff(psi_container::PSIContainer,
                        cons_name::Symbol,
                        param_reference::UpdateRef,
                        var_name::Symbol)

Constructs a parametrized integral limit constraint to implement feedforward from other models.
The Parameters are initialized using the upper boundary values of the provided variables.

# Constraints
``` sum(variable[var_name, t] for t in time_steps)/length(time_steps) <= param_reference[var_name] ```

# LaTeX

`` \sum_{t} x \leq param^{max}``
TO DO: New formulation when DataStore available
`` \sum_{t} x * DeltaT_lower \leq param^{max} * DeltaT_upper ``
TO DO: New formulation when Commitment is considered: SemiContinuousFF
    `` P_LL - P_max * ON_upper <= 0.0 ``
    `` P_LL - P_min * ON_upper >= 0.0 ``

# Arguments
* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* cons_name::Symbol : name of the constraint
* param_reference : Reference to the PJ.ParameterRef used to determine the upperbound
* var_name::Symbol : the name of the continuous variable
"""
function integral_limit_ff(
    psi_container::PSIContainer,
    cons_name::Symbol,
    param_reference::UpdateRef,
    var_name::Symbol,
)
    time_steps = model_time_steps(psi_container)
    ub_name = middle_rename(cons_name, PSI_NAME_DELIMITER, "integral_limit")
    variable = get_variable(psi_container, var_name)

    axes = JuMP.axes(variable)
    set_name = axes[1]

    @assert axes[2] == time_steps
    container_ub = add_param_container!(psi_container, param_reference, set_name)
    param_ub = get_parameter_array(container_ub)
    con_ub = add_cons_container!(psi_container, ub_name, set_name)

    for name in axes[1]
        value = JuMP.upper_bound(variable[name, 1])

        param_ub[name] = PJ.add_parameter(psi_container.JuMPmodel, value)
        con_ub[name] = JuMP.@constraint(
            psi_container.JuMPmodel,
            sum(variable[name, t] for t in time_steps) / length(time_steps) <=
            param_ub[name]
        )
    end
end

########################## FeedForward Constraints #########################################
function feedforward!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::Nothing,
) where {T <: PSY.Component}
    return
end

function feedforward!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::UpperBoundFF,
) where {T <: PSY.StaticInjection}
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_active_power_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        ub_ff(
            psi_container,
            constraint_name(FEEDFORWARD_UB, T),
            constraint_infos,
            parameter_ref,
            var_name,
        )
    end
end

function feedforward!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    model::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::SemiContinuousFF,
) where {T <: PSY.StaticInjection}
    bin_var = make_variable_name(get_binary_source_stage(ff_model), T)
    parameter_ref = UpdateRef{JuMP.VariableRef}(bin_var)
    constraint_infos = Vector{DeviceRangeConstraintInfo}(undef, length(devices))
    for (ix, d) in enumerate(devices)
        name = PSY.get_name(d)
        limits = PSY.get_active_power_limits(d)
        constraint_info = DeviceRangeConstraintInfo(name, limits)
        add_device_services!(constraint_info, d, model)
        constraint_infos[ix] = constraint_info
    end
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        semicontinuousrange_ff(
            psi_container,
            make_constraint_name(FEEDFORWARD_BIN, T),
            constraint_infos,
            parameter_ref,
            var_name,
        )
    end
end

function feedforward!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, <:AbstractDeviceFormulation},
    ff_model::IntegralLimitFF,
) where {T <: PSY.StaticInjection}
    for prefix in get_affected_variables(ff_model)
        var_name = make_variable_name(prefix, T)
        parameter_ref = UpdateRef{JuMP.VariableRef}(var_name)
        integral_limit_ff(
            psi_container,
            make_constraint_name(FEEDFORWARD_INTEGRAL_LIMIT, T),
            parameter_ref,
            var_name,
        )
    end
end

#########################FeedForward Variables Updating#####################################
# This makes the choice in which variable to get from the results.
function get_stage_variable(
    chron::RecedingHorizon,
    stages::Pair{Stage{T}, Stage{U}},
    device_name::AbstractString,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    step = axes(variable)[2][chron.periods]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_stage_variable(
    ::Consecutive,
    stages::Pair{Stage{T}, Stage{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    step = axes(variable)[2][get_end_of_interval_step(stages.first)]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_stage_variable(
    chron::Synchronize,
    stages::Pair{Stage{T}, Stage{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    e_count = get_execution_count(stages.second)
    wait_count = get_execution_wait_count(get_trigger(chron))
    index = (floor(e_count / wait_count) + 1)
    step = axes(variable)[2][Int(index)]
    var = variable[device_name, step]
    if JuMP.is_binary(var)
        return round(JuMP.value(var))
    else
        return JuMP.value(var)
    end
end

function get_stage_variable(
    ::FullHorizon,
    stages::Pair{Stage{T}, Stage{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    vars = variable[device_name, :]
    if JuMP.is_binary(first(vars))
        return round.(JuMP.value(vars))
    else
        return JuMP.value.(vars)
    end
end

function get_stage_variable(
    chron::Range,
    stages::Pair{Stage{T}, Stage{U}},
    device_name::String,
    var_ref::UpdateRef,
) where {T, U <: AbstractOperationsProblem}
    variable = get_variable(stages.first.internal.psi_container, var_ref.access_ref)
    vars = variable[device_name, chron.range]
    if JuMP.is_binary(first(vars))
        return round.(JuMP.value(vars))
    else
        return JuMP.value.(vars)
    end
end

function feedforward_update!(
    destination_stage::Stage,
    source_stage::Stage,
    chronology::FeedForwardChronology,
    param_reference::UpdateRef{JuMP.VariableRef},
    param_array::JuMPParamArray,
    current_time::Dates.DateTime,
)
    trigger = get_trigger(chronology)
    if trigger_update(trigger)
        for device_name in axes(param_array)[1]
            var_value = get_stage_variable(
                chronology,
                (source_stage => destination_stage),
                device_name,
                param_reference,
            )
            previous_value = PJ.value(param_array[device_name])
            PJ.fix(param_array[device_name], var_value)
            IS.@record :simulation FeedForwardUpdateEvent(
                "FeedForward",
                current_time,
                param_reference,
                device_name,
                var_value,
                previous_value,
                destination_stage,
                source_stage,
            )
        end
        reset_trigger_count!(trigger)
    end
    update_count!(trigger)
    return
end


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

"""
    BlockCrossValidation(sides; loss=Dict())

Cross-validation with blocks of given `sides`. Optionally,
specify `loss` function from `LossFunctions.jl` for some
of the variables. If only one side is provided, then blocks
become cubes.

## References

* Roberts et al. 2017. [Cross-validation strategies for data with
  temporal, spatial, hierarchical, or phylogenetic structure]
  (https://onlinelibrary.wiley.com/doi/10.1111/ecog.02881)
* Pohjankukka et al. 2017. [Estimating the prediction performance
  of spatial models via spatial k-fold cross-validation]
  (https://www.tandfonline.com/doi/full/10.1080/13658816.2017.1346255)
"""
struct BlockCrossValidation{S} <: AbstractErrorEstimator
  sides::S
  loss::Dict{Symbol,SupervisedLoss}
end

BlockCrossValidation(sides; loss=Dict()) =
  BlockCrossValidation{typeof(sides)}(sides, loss)

function error(solver::AbstractLearningSolver,
               problem::LearningProblem,
               eestimator::BlockCrossValidation)
  sdata = sourcedata(problem)
  ovars = outputvars(task(problem))
  loss  = eestimator.loss
  sides = eestimator.sides
  for var in ovars
    if var ∉ keys(loss)
      loss[var] = defaultloss(sdata[var][1])
    end
  end

  bsides    = length(sides) > 1 ? sides : ntuple(i->sides, ncoords(sdata))
  blocks    = partition(sdata, BlockPartitioner(bsides))
  neighbors = metadata(blocks)[:neighbors]
  bsubsets  = subsets(blocks)
  allblocks = 1:length(blocks)

  # error for a block b
  function ε(b)
    # source and target blocks
    sblocks = setdiff(allblocks, [neighbors[b]; b])
    tblocks = [b]

    # source and target indices
    sinds = reduce(vcat, bsubsets[sblocks])
    tinds = reduce(vcat, bsubsets[tblocks])

    # source and target data
    train = view(sdata, sinds)
    hold  = view(sdata, tinds)

    # setup and solve sub-problem
    subproblem = LearningProblem(train, hold, task(problem))
    solution   = solve(subproblem, solver)

    # loss for each variable
    losses = map(ovars) do var
      y = hold[var]
      ŷ = solution[var]
      ℒ = value(loss[var], y, ŷ, AggMode.Mean())
      var => ℒ
    end

    Dict(losses)
  end

  # compute error for each fold in parallel
  εs = foldxt(vcat, Map(ε), allblocks)

  # combine error from different folds
  Dict(var => mean(get.(εs, var, 0)) for var in ovars)
end


text:
module TestTypeCheck
  using TypeCheck, FactCheck

  istype(t) = isa(t,TypeCheck.AType)

  facts("Check Return Types: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.check_return_types(f) => anything # => FunctionSignature([],Symbol)
            [@fact istype(TypeCheck.returntype(e)) => true for e in code_typed(f)]
            [@fact TypeCheck.returntype(e) => istype for e in code_typed(f)]
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  caught(x) = x[2] == true
  notcaught(x) = x[2] == false
  function check_return(f,check)
    @fact length(code_typed(f)) => 1
    @fact TypeCheck.check_return_type(code_typed(f)[1]) => check
  end

  facts("Check Return Types: True Positives") do
    barr(x::Int) = isprime(x) ? x : false
    check_return(barr, caught)
  end

  facts("Check Return Types: False Negatives") do
    foo(x::Any) = isprime(x) ? x : false
    check_return(foo, notcaught)
  end

  facts("Check Loop Types: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.check_loop_types(f) => anything # => LoopResults
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  passed(x) = isempty(x.methods)
  failed(x) = !passed(x)
  function check_loops(f,check)
    @fact length(code_typed(f)) => 1
    @fact check_loop_types(f) => check
  end

  facts("Check Loop Types: True Positives") do
    function f1(x::Int)
      for n in 1:x
        x /= n
      end
      return x
    end
    check_loops(f1,failed)
  end

  facts("Check Loop Types: True Negatives") do
    function g1()
      x::Int = 5
      for i = 1:100
        x *= 2.5
      end
      return x
    end
    check_loops(g1,passed)
    function g2()
      x = 5
      x = 0.2
      for i = 1:10
        x *= 2
      end
      return x
    end
    check_loops(g2,passed)
  end

  facts("Check Method Calls: Make Sure It Runs on Base") do
    for n in names(Base)
      if isdefined(Base,n)
        f = eval(Base,n)
        if isgeneric(f) && typeof(f) == Function
          context(string(n)) do
            @fact TypeCheck.check_method_calls(f) => anything # => FunctionCalls
          end
        end
      else
        @fact n => x->isdefined(Base,x)
      end
    end
  end

  facts("Check Method Calls: True Positives") do
  end

  facts("Check Method Calls: False Negatives") do
  end

  exitstatus()
end


text:
test_data_name = "eth_test.pkl";                                                    # test data set name
test_scene_id = 0;                                                                  # test data id
start_time_idx = 905;                                                               # start time index in test data
ego_pos_init_vec = [5., 0.5] .+ [-5.263534, -5.314636];                             # initial ego position [x, y] [m]
ego_pos_goal_vec = [5., 8.9] .+ [-5.263534, -5.314636];                             # goal ego position [x, y] [m]
target_speed = 0.7;                                                                 # target speed [m/s]
sim_horizon = 10.0;                                                                 # simulation time horizon [s]


text:
# Practical Machine learning
# Decision Trees based learning - Random Forests example
# Chapter 5

module Trees

using StatsBase
using DataFrames

include("example.jl")
include("sort.jl")
include("split.jl")

export Tree, fit, predict

abstract Element

type Node{T} <: Element
    feature::Int
    threshold::T
    impurity::Float64
    n_samples::Int
    left::Int
    right::Int
end

abstract Leaf <: Element

type ClassificationLeaf <: Leaf
    counts::Vector{Int}
    impurity::Float64
    n_samples::Int

    function ClassificationLeaf(example::Example, samples::Vector{Int}, impurity::Float64)
        counts = zeros(Int, example.n_labels)
        for s in samples
            label = example.y[s]
            counts[label] += round(Int, example.sample_weight[s])
        end
        new(counts, impurity, length(samples))
    end
end

majority(leaf::ClassificationLeaf) = indmax(leaf.counts)

type RegressionLeaf <: Leaf
    mean::Float64
    impurity::Float64
    n_samples::Int

    function RegressionLeaf(example::Example, samples::Vector{Int}, impurity::Float64)
        new(mean(example.y[samples]), impurity, length(samples))
    end
end

Base.mean(leaf::RegressionLeaf) = leaf.mean

immutable Undef <: Element; end
const undef = Undef()

# parameters to build a tree
immutable Params
    criterion::Criterion
    splitter  # splitter constructor
    max_features::Int
    max_depth::Int
    min_samples_split::Int
end

# bundled arguments for splitting a node
immutable SplitArgs
    index::Int
    depth::Int
    range::UnitRange{Int}
end

type Tree
    index::Int
    nodes::Vector{Element}

    Tree() = new(0, Element[])
end

getnode(tree::Tree, index::Int) = tree.nodes[index]
getroot(tree::Tree) = getnode(tree, 1)
getleft(tree::Tree, node::Node) = getnode(tree, node.left)
getright(tree::Tree, node::Node) = getnode(tree, node.right)

isnode(tree::Tree, index::Int) = isa(tree.nodes[index], Node)
isleaf(tree::Tree, index::Int) = isa(tree.nodes[index], Leaf)
isnode(node::Element) = isa(node, Node)
isleaf(node::Element) = isa(node, Leaf)

impurity(node::Node) = node.impurity
impurity(leaf::Leaf) = leaf.impurity
n_samples(node::Node) = node.n_samples
n_samples(leaf::Leaf) = leaf.n_samples

function next_index!(tree::Tree)
    push!(tree.nodes, undef)
    tree.index += 1
end

function fit(tree::Tree, example::Example, criterion::Criterion, max_features::Int, max_depth::Int, min_samples_split::Int)
    if isa(criterion, ClassificationCriterion)
        splitter = ClassificationSplitter{typeof(criterion)}
    elseif isa(criterion, RegressionCriterion)
        splitter = RegressionSplitter{typeof(criterion)}
    else
        error("invalid criterion")
    end
    params = Params(criterion, splitter, max_features, max_depth, min_samples_split)
    samples = where(example.sample_weight)
    sample_range = 1:length(samples)
    next_index!(tree)
    args = SplitArgs(tree.index, 1, sample_range)
    build_tree(tree, example, samples, args, params)
    return
end

function where(v::AbstractVector)
    n = countnz(v)
    indices = Array(Int, n)
    i = 1
    j = 0

    while (j = findnext(v, j + 1)) > 0
        indices[i] = j
        i += 1
    end

    indices
end

function leaf(example::Example, samples, criterion::RegressionCriterion)
    RegressionLeaf(example, samples, impurity(samples, example, criterion))
end

function leaf(example::Example, samples, criterion::ClassificationCriterion)
    ClassificationLeaf(example, samples, impurity(samples, example, criterion))
end

function build_tree(tree::Tree, example::Example, samples::Vector{Int}, args::SplitArgs, params::Params)
    n_features = example.n_features
    range = args.range  # shortcut
    n_samples = length(range)

    if args.depth >= params.max_depth || n_samples < params.min_samples_split
        tree.nodes[args.index] = leaf(example, samples[range], params.criterion)
        return
    end

    best_feature = 0
    best_impurity = Inf
    local best_threshold, best_boundary

    for k in sample(1:n_features, params.max_features, replace=false)
        feature = example.x[samples[range], k]
        sort!(samples, feature, range)
        splitter = params.splitter(samples, feature, range, example)

        for s in splitter
            averaged_impurity = (s.left_impurity * s.n_left_samples + s.right_impurity * s.n_right_samples) / (s.n_left_samples + s.n_right_samples)

            if averaged_impurity < best_impurity
                best_impurity = averaged_impurity
                best_feature = k
                best_threshold = s.threshold
                best_boundary = s.boundary
            end
        end
    end

    if best_feature == 0
        tree.nodes[args.index] = leaf(example, samples[range], params.criterion)
    else
        feature = example.x[samples[range], best_feature]
        sort!(samples, feature, range)

        left = next_index!(tree)
        right = next_index!(tree)
        tree.nodes[args.index] = Node(best_feature, best_threshold, best_impurity, n_samples, left, right)

        next_depth = args.depth + 1
        left_node = SplitArgs(left, next_depth, range[1:best_boundary])
        right_node = SplitArgs(right, next_depth, range[best_boundary+1:end])
        build_tree(tree, example, samples, left_node, params)
        build_tree(tree, example, samples, right_node, params)
    end

    return
end

function predict(tree::Tree, x::AbstractVector)
    node = getroot(tree)

    while true
        if isa(node, Node)
            if x[node.feature] <= node.threshold
                # go left
                node = getleft(tree, node)
            else
                # go right
                node = getright(tree, node)
            end
        elseif isa(node, ClassificationLeaf)
            return majority(node)
        elseif isa(node, RegressionLeaf)
            return mean(node)
        else
            error("found invalid type of node (type: $(typeof(node)))")
        end
    end
end

end  # module Trees


text:
function test_multicombination_combinations()
    @test collect(combinations(Int[], Int[], 0)) == Array{Int,1}[[]]
    @test length(combinations(Int[], Int[], 0)) == 1
    @test eltype(combinations(Int[], Int[], 0)) == Array{Int,1}

    @test collect(combinations(Int[], Int[], 1)) == Array{Int,1}[]
    @test length(combinations(Int[], Int[], 1)) == 0
    @test eltype(combinations(Int[], Int[], 1)) == Array{Int,1}

    @test collect(combinations(["a", "b"], [2, 1], 3)) == Array[
        ["a", "a", "b"],
    ]
    @test length(combinations(["a", "b"], [2, 1], 3)) == 1
    @test eltype(combinations(["a", "b"], [2, 1], 3)) == Array{String,1}

    @test collect(combinations(["a", "b"], [3, 2], 3)) == Array[
        ["a", "a", "a"],
        ["a", "a", "b"],
        ["a", "b", "b"],
    ]
    @test length(combinations(["a", "b"], [3, 2], 3)) == 3
    @test eltype(combinations(["a", "b"], [3, 2], 3)) == Array{String,1}

    @test collect(combinations(["P", "O", "E", "S"], [1, 1, 2, 5], 4)) == Array[
        ["P", "O", "E", "E"],
        ["P", "O", "E", "S"],
        ["P", "O", "S", "S"],
        ["P", "E", "E", "S"],
        ["P", "E", "S", "S"],
        ["P", "S", "S", "S"],
        ["O", "E", "E", "S"],
        ["O", "E", "S", "S"],
        ["O", "S", "S", "S"],
        ["E", "E", "S", "S"],
        ["E", "S", "S", "S"],
        ["S", "S", "S", "S"],
    ]
    @test length(combinations(["P", "O", "E", "S"], [1, 1, 2, 5], 4)) == 12
    @test eltype(combinations(["P", "O", "E", "S"], [1, 1, 2, 5], 4)) == Array{String,1}
end

function test_multicombination_all()
    print(rpad("Math.Combinatorics.Combination[Multiset]...", 50, ' '))

    test_multicombination_combinations()

    println("PASS")
end


text:
#=
@enum(Monotonicity, MONO_CONST, MONO_NONINCR, MONO_NONDECR, MONO_NONE, MONO_UNSET)

function *(m::Monotonicity, v::Vexity)
    if m == MONO_NONDECR
        return v
    elseif m == MONO_NONINCR
        return -v
    elseif m == MONO_NONE
        if v == VEX_CONVEX
            return VEX_NONDCP
        elseif v == VEX_CONCAVE
            return VEX_NONDCP
        else
            return v
        end
    end
    return VEX_NONDCP # TODO CHECK FALLBACK RETURN
end

monotonicity(f, x::Interval) = error("monotonicity(f,x) for $f is not implemented.")

function monotonicity(df::Interval{T}, x::Interval{T}) where T
    if zero(T) <= df
        return (MONO_NONDECR,)
    elseif zero(T) >= df
        return (MONO_NONINCR,)
    end
    return (MONO_NONE,)
end

for F in (positive, negative, lower_bnd, upper_bnd, bnd)
    @eval monotonicity(::typeof($F), x::Interval) = MONO_AFFINE
end

function monotonicity(::typeof(+), x::Interval, y::Interval, args...)
    ntuple(i -> MONO_NONDECR, length(args) + 2)
end
monotonicity(::typeof(+), x::Interval) = (MONO_NONDECR)

monotonicity(::typeof(-), x::Interval, y::Interval) = (MONO_NONDECR, MONO_NONINCR)
monotonicity(::typeof(-), x::Interval) = (MONO_NONINCR,)

function monotonicity(::typeof(/), x::Interval{T}, y::Interval{T}) where T
    if v > zero(T)
        return (MONO_AFFINE, MONO_NONINCR)
    elseif v < zero(T)
        return (MONO_AFFINE, MONO_NONDECR)
    end
    return (MONO_AFFINE, MONO_NONE)
end

function monotonicity(::typeof(*), args...)
    return ntuple(i -> prod_exclude(_sign, args, i)*MONO_NONDECR, length(args))
end

for F in (abs, abs2, cosh, sech)
    @eval function monotonicity(::typeof($F), x::Interval{T}) where T
        if zero(T) <= x
            return (MONO_NONDECR,)
        elseif zero(T) >= x
            return (MONO_NONINCR,)
        end
        return (MONO_NONE,)
    end
end

for F in (exp, exp2, exp10, expm1, sinh, tanh,
          log, log2, log10, log1p,
          asin, atan,
          erf, erfinv,
          sqrt, cbrt, xabsx,
          # rectifer & sigmoid activation functions
          relu, leaky_relu, maxtanh, maxsig, pentanh, sigmoid,
          bisigmoid, softsign, softplus, elu)
    @eval monotonicity(::typeof($F), x::Interval) = (MONO_NONDECR,)
end

for F in (acos, acot, erfc, erfcinv)
    @eval monotonicity(::typeof($F), x::Interval) = (MONO_NONINCR,)
end

for F in (coth, csch)
    @eval function monotonicity(::typeof($F), x::Interval{T}) where T
        if zero(T) <= x
            return (MONO_NONINCR,)
        elseif zero(T) >= x
            return (MONO_NONDECR,)
        end
        return (MONO_NONE,)
    end
end

monotonicity(::typeof(sin), x::Interval) = monotonicity(cos(x), x)
monotonicity(::typeof(cos), x::Interval) = monotonicity(sin(x), x)
monotonicity(::typeof(tan), x::Interval) = monotonicity(sec(x)^2, x)

function monotonicity(::typeof(asec), x::Interval{T}) where T
    if one(T) <= x
        return (MONO_NONDECR,)
    elseif -one(T) >= x
        return (MONO_NONINCR,)
    end
    return (MONO_NONE,)
end

function monotonicity(::typeof(acsc), x::Interval{T}) where T
    if one(T) <= x
        return (MONO_NONINCR,)
    elseif -one(T) >= x
        return (MONO_NONDECR,)
    end
    return (MONO_NONE,)
end

function monotonicity(::typeof(xlogx), x::Interval{T}) where T
    if one(T)/MathConstants.e <= x
        return (MONO_NONDECR,)
    elseif one(T)/MathConstants.e >= x
        return (MONO_NONINCR,)
    end
    return (MONO_NONE,)
end

function monotonicity(::typeof(xexpax), x::Interval{T}, a::Interval{T}) where T
    # TODO:
    return (MONO_NONE,)
end

# arity 1 implies EXPRESSION OR SELECT
function monotonicity(::Val{1}, n::NodeType, atype::AtomType, bnds::Interval)
    if n == SELECT
        error("Not implemented yet...")
    end
    if n == EXPRESSION
        # binary_switch on univariate atype
    end
end

# arity >= 2 implies EXPRESSION
function monotonicity(n::NodeType, atype::AtomType, bnds::Interval, i::Int)
    # binary_switch... atype
end

function populate_monotonicity!(n, buffer, indx)
    mono = _mono(buffer, indx)
    if _arity(n) == 1
        @inbounds mono[1] = monotonicity(Val(1),
                                         _node_type(n),
                                         _atom_type(n),
                                         _bnd(buffer, _child_id(buffer, indx)))
    else
        unsafe_map!(monotonicity, mono, n, _arity(n))
    end
    return nothing
end

#=
Gets monotonicity of buffer at i. Recalling from a cached value if
discovered or locked. Otherwise, populating the mono[i] field.
=#
function monotonicity(buffer::ConvexityBuffer, indx)
    if _discovered(buffer, indx)
        return _mono(buffer, indx)
    end
    populate_monotonicity!(_node(dag, indx), buffer, indx)
    return _mono(buffer, i)
end
=#

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
module Device
using CUDA

struct CPU end
struct CUDADevice end

array_type() = CUDA.has_cuda_gpu() ? CuArray : Array
array_device(::Union{Array}) = CPU()
array_device(::CuArray) = CUDADevice()

end

max_stars_repo_name: climate-machine/JRRTMGP

max_forks_repo_name: climate-machine/JRRTMGP


text:
export CircuitComponent
"""
An abstract representation of a circuit component e.g. a capacitor, inductor, or
transmission line.
"""
abstract type CircuitComponent end

# the empty string is used for the ground net
ground = ""

"""
    Blackbox(ω::Vector{<:Real}, comp::CircuitComponent)
    Blackbox(ω::Vector{<:Real}, comp::TransmissionLine)
    Blackbox(ω::Vector{<:Real}, comp::Cascade)

Create a Blackbox model for a CircuitComponent.
"""
function Blackbox(ω::Vector{<:Real}, comp::CircuitComponent)
    return canonical_gauge(Blackbox(ω, PSOModel(comp)))
end

export SeriesComponent
"""
    SeriesComponent(p1::String, p2::String, inv_inductance::Float64, conductance::Float64,
        capacitance::Float64)
    SeriesComponent(p1::String, p2::String, k::Real, g::Real, c::Real)

A parallel inductor, capacitor, and resistor between two ports.
"""
struct SeriesComponent <: CircuitComponent
    p1::String
    p2::String
    inv_inductance::Float64
    conductance::Float64
    capacitance::Float64
    function SeriesComponent(p1::String, p2::String, inv_inductance::Float64,
        conductance::Float64, capacitance::Float64)
        @assert p1 != ground
        @assert p2 != ground
        return new(p1, p2, inv_inductance, conductance, capacitance)
    end
end

SeriesComponent(p1::String, p2::String, k::Real, g::Real, c::Real) =
    SeriesComponent(p1, p2, 1.0 * k, 1.0 * g, 1.0 * c)

"""
    Circuit(comp::SeriesComponent)
    Circuit(comp::ParallelComponent)
    Circuit(comp::TransmissionLine, vertex_prefix::AbstractString,
        stages_per_segment::Union{Nothing, AbstractVector{Int}}=nothing)

Create a Circuit for a CircuitComponent.
"""
function Circuit(comp::SeriesComponent)
    c = Circuit([ground, comp.p1, comp.p2])
    set_inv_inductance!(c, comp.p1, comp.p2, comp.inv_inductance)
    set_conductance!(c, comp.p1, comp.p2, comp.conductance)
    set_capacitance!(c, comp.p1, comp.p2, comp.capacitance)
    return c
end

"""
    PSOModel(comp::SeriesComponent)
    PSOModel(comp::ParallelComponent)
    PSOModel(comp::TransmissionLine)
    PSOModel(comp::Cascade)

Create a PSOModel for a CircuitComponent.
"""
function PSOModel(comp::SeriesComponent)
    return PSOModel(Circuit(comp), [(comp.p1, ground), (comp.p2, ground)], [comp.p1, comp.p2])
end

export ParallelComponent
"""
    ParallelComponent(p::String, inv_inductance::Float64,
        conductance::Float64, capacitance::Float64)
    ParallelComponent(p::String, k::Real, g::Real, c::Real)

A parallel inductor, capacitor, and resistor with one port.
"""
struct ParallelComponent <: CircuitComponent
    p::String
    inv_inductance::Float64
    conductance::Float64
    capacitance::Float64
    function ParallelComponent(p::String, inv_inductance::Float64,
        conductance::Float64, capacitance::Float64)
        @assert p != ground
        return new(p, inv_inductance, conductance, capacitance)
    end
end

ParallelComponent(p::String, k::Real, g::Real, c::Real) =
    ParallelComponent(p, 1.0 * k, 1.0 * g, 1.0 * c)

function Circuit(comp::ParallelComponent)
    c = Circuit([ground, comp.p])
    set_inv_inductance!(c, comp.p, ground, comp.inv_inductance)
    set_conductance!(c, comp.p, ground, comp.conductance)
    set_capacitance!(c, comp.p, ground, comp.capacitance)
    return c
end

PSOModel(comp::ParallelComponent) = PSOModel(Circuit(comp), [(comp.p, ground)], [comp.p])

export TransmissionLine
"""
    TransmissionLine(ports::Vector{String}, propagation_speed::Float64,
        characteristic_impedance::Float64, len::Float64,
        locations::AbstractVector{Float64}, δ::Float64)
    TransmissionLine(ports::Vector{String}, propagation_speed::Real,
        characteristic_impedance::Real, len::Real; locations::Vector{<:Real}=Float64[],
        δ::Real=len/100)

A transmission line with given propagation speed and characteristic impedance. Ports will
be present at the two ends of the transmission line in addition to the `locations`. `δ` is
the maximum length of an LC stage in the circuit approximation to a transmission line. It
is ignored when constructing a Blackbox from a TransmissionLine.
"""
struct TransmissionLine <: CircuitComponent
    ports::Vector{String}
    propagation_speed::Float64
    characteristic_impedance::Float64
    len::Float64 # length, but that's the name of a built in function
    port_locations::UniqueVector{Float64}
    δ::Float64 # ignored when converted directly to Blackbox
    function TransmissionLine(ports::Vector{String}, propagation_speed::Float64,
        characteristic_impedance::Float64, len::Float64, locations::AbstractVector{Float64},
        δ::Float64)
        @assert !(ground in ports)
        @assert propagation_speed > 0
        @assert characteristic_impedance > 0
        @assert len > 0
        @assert all(locations .> 0)
        @assert all(locations .< len)
        port_locations = UniqueVector([0; locations; len])
        @assert length(ports) == length(port_locations)
        @assert δ >= 0
        return new(ports, propagation_speed, characteristic_impedance,
            len, port_locations, δ)
    end
end

TransmissionLine(ports::Vector{String}, propagation_speed::Real,
    characteristic_impedance::Real, len::Real; locations::Vector{<:Real}=Float64[],
    δ::Real=len/100) = TransmissionLine(ports, propagation_speed,
    characteristic_impedance, len, locations, δ)

function Circuit_and_ports(comp::TransmissionLine,
    stages_per_segment::Union{Nothing, AbstractVector{Int}}=nothing)
    sort_inds = sortperm(comp.port_locations)
    port_locations, ports = comp.port_locations[sort_inds], comp.ports[sort_inds]
    ν, Z0, δ = comp.propagation_speed, comp.characteristic_impedance, comp.δ
    segment_lengths = port_locations[2:end] - port_locations[1:end-1]
    if stages_per_segment == nothing
        stages_per_segment = [ceil(Int, len / δ) for len in segment_lengths]
    end
    num_vertices = sum(stages_per_segment) + 2 # ground is 0
    c = Circuit(0:(num_vertices-1))
    v0 = 0
    port_edges = Tuple{Int, Int}[]
    for (segment_length, num_stages) in zip(segment_lengths, stages_per_segment)
        stage_length = 1.0 * segment_length / num_stages
        inductance = stage_length * Z0 / ν
        capacitance = stage_length / (Z0 * ν)
        for v in (v0 + 1):(v0 + num_stages)
            set_inductance!(c, v, v + 1, inductance)
        end
        for v in [v0 + 1, v0 + num_stages + 1]
            set_capacitance!(c, v, 0, 0.5 * capacitance + get_capacitance(c, v, 0))
        end
        for v in (v0 + 2):(v0 + num_stages)
            set_capacitance!(c, v, 0, capacitance)
        end
        push!(port_edges, (v0 + 1, 0))
        v0 += num_stages
    end
    @assert v0 + 1 == num_vertices-1
    push!(port_edges, (num_vertices-1, 0))
    return c, port_edges, ports
end

function Circuit(comp::TransmissionLine, vertex_prefix::AbstractString,
    stages_per_segment::Union{Nothing, AbstractVector{Int}}=nothing)
    c, port_edges, ports = Circuit_and_ports(comp, stages_per_segment)
    vertices = [ground; ["$(vertex_prefix)$v" for v in c.vertices[2:end]]]
    for (e, p) in zip(port_edges, ports)
        vertices[e[1]+1] = p
    end
    return Circuit(matrices(c)..., vertices)
end

PSOModel(comp::TransmissionLine) = PSOModel(Circuit_and_ports(comp)...)

function Blackbox(ω::Vector{<:Real}, comp::TransmissionLine)
    ν = comp.propagation_speed
    Z0 = comp.characteristic_impedance
    function Blackbox_from_params(ports::Vector, len::Real)
        u = exp.(-1im * ω * len/ν)
        # S matrix is [[0, u] [u, 0]]
        # turn this into Y matrix using closed form expression
        denom = (1 .- u.^2) * Z0
        Y11_Y22 = (1 .+ u.^2) ./ denom
        Y12_Y21 = -2 * u ./ denom
        Y = [[[a, b] [b, a]] for (a, b) in zip(Y11_Y22, Y12_Y21)]
        U = eltype(eltype(Y))
        P = Matrix{U}(I, 2, 2)
        return Blackbox(ω, Y, P, P, ports)
    end
    sort_inds = sortperm(comp.port_locations)
    sort_port_locations = comp.port_locations[sort_inds]
    sort_ports = comp.ports[sort_inds]
    params = [(ports = [sort_ports[i], sort_ports[i + 1]],
              len = sort_port_locations[i+1] - sort_port_locations[i])
              for i in 1:(length(comp.ports)-1)]
    return cascade_and_unite([Blackbox_from_params(p...) for p in params])
end

export PortOperation
"""
The functions allowed in a Cascade.
"""
PortOperation = Union{typeof(short_ports),
                      typeof(short_ports_except),
                      typeof(open_ports),
                      typeof(open_ports_except),
                      typeof(unite_ports),
                      typeof(canonical_gauge)}

export Cascade
"""
    Cascade(components::Vector{CircuitComponent},
        operations::Vector{Pair{PortOperation, Vector{String}}})

A system created by cascading and uniting several CircuitComponents and then applying
PortOperations to the ports.
"""
struct Cascade <: CircuitComponent
    components::Vector{CircuitComponent}
    operations::Vector{Pair{PortOperation, Vector{String}}}
end

function PSOModel(comp::Cascade)
    model = cascade_and_unite([PSOModel(m) for m in comp.components])
    for (op, args) in comp.operations
        model = op(model, args...)
    end
    return model
end

function Blackbox(ω::Vector{<:Real}, comp::Cascade)
    model = cascade_and_unite([Blackbox(ω, m) for m in comp.components])
    for (op, args) in comp.operations
        model = op(model, args...)
    end
    return model
end


text:
    """
    Cartesian 3×1 vector

        Vector3()
        Vector3(x,y,z)
        Vector3(vector)

    """
    const Vector3 = SVector{3, Float64}

    Vector3() = Vector3(fill(0.0, 3))

    export ô, î, ĵ, k̂
    """
    Origin vector ``\\hat{o}=\\pmatrix{0\\\\0\\\\0}``

        const ô = Vector3(0,0,0)

    """
    const ô = Vector3(0,0,0)
    """
    X-axis basis vector ``\\hat{\\imath}=\\pmatrix{1\\\\0\\\\0}``

        const î = Vector3(1,0,0)

    """
    const î = Vector3(1,0,0)
    """
    Y-axis basis vector ``\\hat{\\jmath}=\\pmatrix{0\\\\1\\\\0}``

        const ĵ = Vector3(0,1,0)

    """
    const ĵ = Vector3(0,1,0)
    """
    Z-axis basis vector ``\\hat{k}=\\pmatrix{0\\\\0\\\\1}``

        const k̂ = Vector3(0,0,1)

    """
    const k̂ = Vector3(0,0,1)

    function  LinearAlgebra.norm(a::Vector3, p::Real=2)
        return LinearAlgebra.norm(a.data, p)
    end

    """
    Returns a normalized a vector whose magnitude is 1.0 or 0.0

        n = normalize(v)

    """
    function LinearAlgebra.normalize(v::Vector3)::Vector3
        m2 = dot(v,v)
        if m2>0 && m2!=1.0
            return v/sqrt(m2)
        else
            return v
        end
    end

    """
    Dot product between vectors ``\\pmatrix{x \\\\ y \\\\ z} \\cdot \\pmatrix{u \\\\ v \\\\ w} = x*u + y*v + z*w``

        dot(a::Vector3, b::Vector3)
        a ⋅ b => dot(a,b)

    """
    @inline dot(a::Vector3, b::Vector3)::Float64 = a[1]*b[1]+a[2]*b[2]+a[3]*b[3]

    """
    Cross product of vectors

        cross(a::Vector3, b::Vector3)::Vector3
        a × b => cross(a,b)

    """
    cross(a::Vector3, b::Vector3)::Vector3 = Vector3(
        a[2]*b[3]-a[3]*b[2],
        a[3]*b[1]-a[1]*b[3],
        a[1]*b[2]-a[2]*b[1])

    """
    Cross product matrix operator ``\\pmatrix{x \\\\ y \\\\ z}\\times =  \\pmatrix{0 & -z & y \\\\ z & 0 & -x \\\\ -y & x & 0}``

        cross(a::Vector3)
        ×(a::Vector3)

    """
    cross(a::Vector3)::Matrix3 = Matrix3(0,-a[3],a[2],a[3],0,-a[1],-a[2],a[1],0)

    """
    Double cross product matrix operator ``\\pmatrix{x \\\\ y \\\\ z}\\times \\pmatrix{x \\\\ y \\\\ z}\\times =  \\pmatrix{-y^2-z^2 & x y & x z \\\\ x y & -x^2-z^2 & y z \\\\ x z & y z & -x^2-z^2}``

        cross2(a::Vector3)

    """
    cross2(a::Vector3)::Matrix3 = Matrix3(
        -a[2]^2-a[3]^2, a[1]*a[2], a[1]*a[3],
        a[1]*a[2], -a[1]^2-a[3]^2, a[2]*a[3],
        a[1]*a[3], a[2]*a[3], -a[1]^2-a[2]^2)

    outer(a::Vector3, b::Vector3)::Matrix3 = Matrix3(
            a[1]*b[1], a[1]*b[2], a[1]*b[3],
            a[2]*b[1], a[2]*b[2], a[2]*b[3],
            a[3]*b[1], a[3]*b[2], a[3]*b[3])

    function Base.angle(a::Vector3, b::Vector3)::Float64
        x = dot(a,b)
        y = cross(a,b).Magnitude
        return atan(y, x)
    end

    """
    Rotational interpolation between two vectors

        slerp(a,b,t)    # where t=0..1

    """
    function slerp(a::Vector3, b::Vector3, t::Float64) :: Vector3
        θ = angle(a, b)
        y = sin(θ)
        return sin((1-t)*θ)/y*a + sin(t*θ)/y*b
    end

    function Base.getproperty(v::Vector3, name::Symbol)
        if name===:X
            return v[1]
        end
        if name===:Y
            return v[2]
        end
        if name===:Z
            return v[3]
        end
        if name===:SumSquares
            return dot(v,v)
        end
        if name===:Magnitude
            return sqrt(dot(v,v))
        end
        return getfield(v, name)
    end

    function Base.propertynames(v::Vector3)
        return (:X, :Y, :Z, :Magnitude, :SumSquares)
    end


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

typealias PairedSequences{S1,S2} Pair{GappedSequence{S1},GappedSequence{S2}}

type AlignmentResult{T,S1,S2}
    score::T
    seqpair::Nullable{PairedSequences{S1,S2}}
end

AlignmentResult{S1,S2}(::Type{S1}, ::Type{S2}, score) = AlignmentResult(score, Nullable{PairedSequences{S1,S2}}())
function AlignmentResult(score, a′::GappedSequence, b′::GappedSequence)
    if length(a′) != length(b′)
        error("gapped sequences have different length")
    end
    AlignmentResult(score, Nullable(a′ => b′))
end

function Base.show(io::IO, aln::AlignmentResult)
    print(io, "Score: ", aln.score)
    if !isnull(aln.seqpair)
        pair = get(aln.seqpair)
        a′ = pair.first
        b′ = pair.second
        print(io, '\n')
        print(io, "  a: ", a′, '\n')
        print(io, "     ", matching_string(a′, b′), '\n')
        print(io, "  b: ", b′)
    end
end

function Base.getindex(aln::AlignmentResult, i::Integer)
    if isnull(aln.seqpair)
        error("alignment is not available")
    end
    pair = get(aln.seqpair)
    if i == 1
        return pair.first
    elseif i == 2
        return pair.second
    end
    throw(BoundsError(i))
end

function matching_string(a, b)
    @assert length(a) == length(b)
    sa = string(a)
    sb = string(b)
    gmap_a = gapmap(a)
    gmap_b = gapmap(b)
    str = Char[]
    for i in 1:length(a)
        if sa[i] == sb[i] && !gmap_a[i] && !gmap_b[i]
            push!(str, '|')
        else
            push!(str, ' ')
        end
    end
    return ASCIIString(str)
end


text:
export coilcraft1010vsgroup

"""
    coilcraft1010vsgroup(partnumber, <keyword arguments>)

Returns a `Group` for a Coilcraft
[1010VS](http://www.coilcraft.com/1010vs.cfm) series inductor.

**Part Numbers**

- 1010VS-23NME
- 1010VS-46NME
- 1010VS-79NME
- 1010VS-111ME
- 1010VS-141ME

**Keyword Arguments**

- `nhinc`, `nwinc`  -- integer number of filaments in height and width
- `rh`, `rw`        -- ratio in the height and width

Terminals are :a and :b.
"""
function coilcraft1010vsgroup(partnumber::AbstractString; nhinc=0, nwinc=0, rh=NaN, rw=NaN)
  # Assumes units = mm !!!
  pndict = Dict("1010VS-23NME" => (3.55, 0.65, 1.5*2π),
                "1010VS-46NME" => (3.55, 0.65, 2.5*2π),
                "1010VS-79NME" => (3.55, 0.65, 3.5*2π),
                "1010VS-111ME" => (3.55, 0.65, 4.5*2π),
                "1010VS-141ME" => (3.55, 0.65, 5.5*2π))
  hn = helixnodes(pndict[partnumber]...)
  transform!(hn,txyz(0,0,2.025))
  b1 = Node(xyz(hn[1])[1], -xyz(hn[1])[1], xyz(hn[1])[3])
  b2 = Node(xyz(hn[1])[1], -xyz(hn[1])[1], 0.325)
  b3 = Node(xyz(hn[1])[1], xyz(hn[1])[1], 0.325)
  t1 = Node(xyz(hn[end])[1], xyz(hn[end])[1], xyz(hn[end])[3])
  t2 = Node(xyz(hn[end])[1], xyz(hn[end])[1], 0.325)
  t3 = Node(xyz(hn[end])[1], -xyz(hn[end])[1], 0.325)
  nodes = [b3;b2;b1;hn;t1;t2;t3]
  sigma = 1e-3*58.5e6 # conductivity of copper
  segments = connectnodes(nodes, SegmentParameters(w=1.7, h=0.6, sigma=sigma,
                          nhinc=nhinc, nwinc=nwinc, rh=rh, rw=rw))
  Group([nodes;segments],Dict(:a=>b3,:b=>t3),Units("mm"))
end


text:
push!(LOAD_PATH, "$(homedir())/code/VinDsl.jl/src")  # needed if VB not a full module

using PDMats
import Base: ==
using FactCheck
using Distributions
using VinDsl

###################################################
# Utility functions for testing
###################################################
# use these equality definitions for testing purposes
==(x::PDMat, y::PDMat) = x.mat == y.mat
function =={D <: Distribution}(x::D, y::D)
    all(f -> x.(f) == y.(f), fieldnames(x))
end
function =={D <: Distribution}(x::RandomNode{D}, y::RandomNode{D})
    all(f -> x.(f) == y.(f), fieldnames(x))
end

###################################################
# Test files to run
###################################################
tests = [
    "typetests",
    "constrainttests",
    "utiltests",
    "structuretests",
    "expressiontests",
    "hmmtests",
    "integrationtests",
    "distributiontests",
    "advitests",
]

print_with_color(:blue, "Running tests:\n")

srand(12345)

for t in tests
    test_fn = "$t.jl"
    print_with_color(:green, "* $test_fn\n")
    include(test_fn)
end
FactCheck.exitstatus()


text:
# Autogenerated wrapper script for libportaudio_jll for aarch64-apple-darwin
export libportaudio

using alsa_jll
JLLWrappers.@generate_wrapper_header("libportaudio")
JLLWrappers.@declare_library_product(libportaudio, "@rpath/libportaudio.dylib")
function __init__()
    JLLWrappers.@generate_init_header(alsa_jll)
    JLLWrappers.@init_library_product(
        libportaudio,
        "lib/libportaudio.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# ---
# title: Pre-grouped data
# cover: assets/pregrouped_data.png
# description: Working with arrays of arrays.
# author: "[Pietro Vertechi](https://github.com/piever)"
# ---

using AlgebraOfGraphics, CairoMakie
set_aog_theme!() #src

x = [rand(10) for i in 1:3]
y = [rand(10) for i in 1:3]
z = [rand(10) for i in 1:3]
c = ["a", "b", "c"]

m = mapping(x, y, color=c => (t -> "Type " * t ) => "Category")
draw(m)

#

m = mapping(x, (y, z) => (+) => "sum", color=c => (t -> "Type " * t ) => "Category")
draw(m)

#

m = mapping(x, [y z], color=dims(1) => renamer(["a", "b", "c"]))
draw(m)

#

m = mapping(x, [y z], color=["1" "2"])
layers = visual(Scatter) + linear()
fg = draw(m * layers)

# save cover image #src
mkpath("assets") #src
save("assets/pregrouped_data.png", fg) #src


text:
using EngEconomics

# Given
ex1Quant = 1275
ex2Quant = 92
fillQuant = 260
pil2Quant = 683

ex1Lab = 1.76
ex2Lab = 27.87
fillLab = 10.18
pil1Lab = 5688
pil2Lab = 10.27

ex1Eq = 1.87
ex2Eq = 6.54
fillEq = 2.25
pil1Eq = 6420
pil2Eq = 9.61

ex1Mat = 0
ex2Mat = 0
fillMat = 0
pil1Mat = 300
pil2Mat = 54.36

# Find
# a.) Total Cost of Structure Excavation
ex2TotalCost = ex2Quant * (ex2Lab + ex2Eq + ex2Mat)

# b.) Total Tower Cost
pil1TotalCost = pil1Eq + pil1Lab + pil1Mat

# c.) Total Cost of Labor
labTotalCost = ex1Lab * ex1Quant + ex2Quant * ex2Lab + fillQuant * fillLab + pil1Lab + pil2Quant * pil2Lab

# d.) Find Costs for Installation of Pillars in 2010
costIndex2007 = 120
costIndex2010 = 138
i = costIndex2010 / costIndex2007
pil2TotalCost = pil2Quant * (pil2Lab + pil2Eq + pil2Mat)
pil1TotalCost2010 = pil2TotalCost * i

# e.) Incremental Cost of ex1 if following changes
marginalLab = 0.03
marginalEq = -0.01
ex1Incremental = ex1Quant * (marginalLab * ex1Lab + marginalEq * ex1Eq)


text:
struct Disturbances
	D # covariance matrix
	w # deterministically sample disturbance input
end

function disturbances(D)
	T = length(D)
	d = size(D[1], 1)
	w = []
	for t = 1:T
		_w = Array(sqrt(D[t]))
		tmp = []
		for i = 1:d
			push!(tmp, _w[:, i])
			push!(tmp, -1.0 * _w[:, i])
		end

		push!(w, tmp)
	end

	return Disturbances(D, w)
end


text:
@testset "Misc/NumberField" begin
  @testset "is_subfield" begin
    Qx, x = FlintQQ["x"]
    K, a = NumberField(x^2 + 1, "a")
    L, b = NumberField(x^4 + 1, "b")

    c, KtoL = is_subfield(K, L)
    @test c == true
    @test parent(KtoL(a)) == L

    c, KtoL = Hecke.is_subfield_normal(K, L)
    @test c == true
    @test parent(KtoL(a)) == L

    OK = maximal_order(K)
    OL = maximal_order(L)
    c, KtoL = is_subfield(K, L)
    @test c == true
    @test parent(KtoL(a)) == L
  end

  @testset "is_isomorphic" begin
    Qx, x = FlintQQ["x"]
    f = x^5 + 12x - 92
    K, a = NumberField(f, "a")

    g = x^5 - 172x^4 + 7024x^3 + 8656448x^2 + 55735552x + 45796197888
    K2, a2 = NumberField(g, "a2")

    c, KtoK2 = is_isomorphic_with_map(K, K2)
    @test c == true
    @test parent(KtoK2(a)) == K2

    OK = maximal_order(K)
    OK2 = maximal_order(K2)
    c, KtoK2 = is_isomorphic_with_map(K, K2)
    @test c == true
    @test parent(KtoK2(a)) == K2

    h = f - 1
    K3, a3 = NumberField(h, "a3")
    @test !is_isomorphic(K, K3)
  end
end


text:
# https://juliagraphics.github.io/Gtk.jl/latest/manual/textwidgets/#Entry-1

using Gtk

# The entry widget allows the user to enter text. The entered text can be read
# and write using
#
# ent = GtkEntry()
# set_gtk_property!(ent,:text,"My String")
# str = get_gtk_property(ent,:text,String)
# The maximum number of characters can be limited using
# set_gtk_property!(ent,:max_length,10).

win = GtkWindow("Text Widgets", 400, 200)
ent = GtkEntry()
push!(win, ent)
set_gtk_property!(ent, :text, "My String")
set_gtk_property!(ent, :max_length, 10)

showall(win)

#
# Sometimes you might want to make the widget non-editable. This can be done by
# a call
#
# # using the accessor methods
# GAccessor.editable(GtkEditable(ent),false)
# # using the property system
# set_gtk_property!(ent,:editable,false)

#set_gtk_property!(ent, :editable, false)

# If you want to use the entry to retrieve passwords you can hide the
# visibility of the entered text. This can be achieve by calling
#
# set_gtk_property!(ent,:visibility,false)

#set_gtk_property!(ent, :visibility, false)

# To get notified by changes to the entry one can listen the "changed" event.

function entryChanged(ent)
    text = get_gtk_property(ent, :text, String)
    println("Text: $text")
end
sid = signal_connect(entryChanged, ent, "changed")

# TODO: setting progress and setting icons in entry


text:
using Dates
using FFTW
using Glob
using SeisIO
using SeisNoise
using Statistics
using Test
import StatsBase: sample

include("test_NoiseData.jl")
include("test_show.jl")
include("test_ArrayFuncs.jl")
include("test_distance.jl")
# include("test_io.jl")
include("test_tools.jl")
include("test_slicing.jl")
include("test_filter.jl")
include("test_phase_shift.jl")
include("test_stacking.jl")
#
# # import  routines for doin' stuff
include("test_compute_fft.jl")
include("test_correlation.jl")
include("test_rotation.jl")
include("test_plotting.jl")
#
# # out with the old
include("test_deprecated.jl")


text:
##### convert the coordinates into fractional indices at grid cell center and face for a individual at (x,y,z)
##### 1-based indexing for grid faces (excluding halo regions)
@inline get_xf_index(x) = x + 1
@inline get_yf_index(y) = y + 1
@inline get_zf_index(z) = z + 1

##### 1-based indexing for grid centers (excluding halo regions)
@inline get_xc_index(x) = x + 0.5
@inline get_yc_index(y) = y + 0.5
@inline get_zc_index(z) = z + 0.5

##### trilinear interpolation
@inline ψ₀₀₀(xd, yd, zd) = (1 - xd) * (1 - yd) * (1 - zd)
@inline ψ₀₀₁(xd, yd, zd) = (1 - xd) * (1 - yd) *      zd
@inline ψ₀₁₀(xd, yd, zd) = (1 - xd) *      yd  * (1 - zd)
@inline ψ₀₁₁(xd, yd, zd) = (1 - xd) *      yd  *      zd
@inline ψ₁₀₀(xd, yd, zd) =      xd  * (1 - yd) * (1 - zd)
@inline ψ₁₀₁(xd, yd, zd) =      xd  * (1 - yd) *      zd
@inline ψ₁₁₀(xd, yd, zd) =      xd  *      yd  * (1 - zd)
@inline ψ₁₁₁(xd, yd, zd) =      xd  *      yd  *      zd

@inline normalized_trilinear_interpolation_u(u, xd, yd, zd, xi, yi, zi, g::AbstractGrid) =
        @inbounds (ψ₀₀₀(xd, yd, zd) * u[xi,   yi,   zi  ] / ΔxC(xi,   yi,   zi,   g) +
                   ψ₀₀₁(xd, yd, zd) * u[xi,   yi,   zi+1] / ΔxC(xi,   yi,   zi+1, g) +
                   ψ₀₁₀(xd, yd, zd) * u[xi,   yi+1, zi  ] / ΔxC(xi,   yi+1, zi,   g) +
                   ψ₀₁₁(xd, yd, zd) * u[xi,   yi+1, zi+1] / ΔxC(xi,   yi+1, zi+1, g) +
                   ψ₁₀₀(xd, yd, zd) * u[xi+1, yi,   zi  ] / ΔxC(xi+1, yi,   zi,   g) +
                   ψ₁₀₁(xd, yd, zd) * u[xi+1, yi,   zi+1] / ΔxC(xi+1, yi,   zi+1, g) +
                   ψ₁₁₀(xd, yd, zd) * u[xi+1, yi+1, zi  ] / ΔxC(xi+1, yi+1, zi,   g) +
                   ψ₁₁₁(xd, yd, zd) * u[xi+1, yi+1, zi+1] / ΔxC(xi+1, yi+1, zi+1, g))

@inline normalized_trilinear_interpolation_v(v, xd, yd, zd, xi, yi, zi, g::AbstractGrid) =
        @inbounds (ψ₀₀₀(xd, yd, zd) * v[xi,   yi,   zi  ] / ΔyC(xi,   yi,   zi,   g) +
                   ψ₀₀₁(xd, yd, zd) * v[xi,   yi,   zi+1] / ΔyC(xi,   yi,   zi+1, g) +
                   ψ₀₁₀(xd, yd, zd) * v[xi,   yi+1, zi  ] / ΔyC(xi,   yi+1, zi,   g) +
                   ψ₀₁₁(xd, yd, zd) * v[xi,   yi+1, zi+1] / ΔyC(xi,   yi+1, zi+1, g) +
                   ψ₁₀₀(xd, yd, zd) * v[xi+1, yi,   zi  ] / ΔyC(xi+1, yi,   zi,   g) +
                   ψ₁₀₁(xd, yd, zd) * v[xi+1, yi,   zi+1] / ΔyC(xi+1, yi,   zi+1, g) +
                   ψ₁₁₀(xd, yd, zd) * v[xi+1, yi+1, zi  ] / ΔyC(xi+1, yi+1, zi,   g) +
                   ψ₁₁₁(xd, yd, zd) * v[xi+1, yi+1, zi+1] / ΔyC(xi+1, yi+1, zi+1, g))

@inline normalized_trilinear_interpolation_w(w, xd, yd, zd, xi, yi, zi, g::AbstractGrid) =
        @inbounds (ψ₀₀₀(xd, yd, zd) * w[xi,   yi,   zi  ] / ΔzC(xi,   yi,   zi,   g) +
                   ψ₀₀₁(xd, yd, zd) * w[xi,   yi,   zi+1] / ΔzC(xi,   yi,   zi+1, g) +
                   ψ₀₁₀(xd, yd, zd) * w[xi,   yi+1, zi  ] / ΔzC(xi,   yi+1, zi,   g) +
                   ψ₀₁₁(xd, yd, zd) * w[xi,   yi+1, zi+1] / ΔzC(xi,   yi+1, zi+1, g) +
                   ψ₁₀₀(xd, yd, zd) * w[xi+1, yi,   zi  ] / ΔzC(xi+1, yi,   zi,   g) +
                   ψ₁₀₁(xd, yd, zd) * w[xi+1, yi,   zi+1] / ΔzC(xi+1, yi,   zi+1, g) +
                   ψ₁₁₀(xd, yd, zd) * w[xi+1, yi+1, zi  ] / ΔzC(xi+1, yi+1, zi,   g) +
                   ψ₁₁₁(xd, yd, zd) * w[xi+1, yi+1, zi+1] / ΔzC(xi+1, yi+1, zi+1, g))

@inline function u_itpl(u, x, y, z, ac, g::AbstractGrid)
    xi = get_xf_index(x) * ac
    yi = get_yc_index(y) * ac
    zi = get_zc_index(z) * ac
    xd, xi = mod(xi, 1), unsafe_trunc(Int, xi)
    yd, yi = mod(yi, 1), unsafe_trunc(Int, yi)
    zd, zi = mod(zi, 1), unsafe_trunc(Int, zi)
    ##### return fraction of grid spacing
    return normalized_trilinear_interpolation_u(u, xd, yd, zd, xi+g.Hx, yi+g.Hy, zi+g.Hz, g)
end

@inline function v_itpl(v, x, y, z, ac, g::AbstractGrid)
    xi = get_xc_index(x) * ac
    yi = get_yf_index(y) * ac
    zi = get_zc_index(z) * ac
    xd, xi = mod(xi, 1), unsafe_trunc(Int, xi)
    yd, yi = mod(yi, 1), unsafe_trunc(Int, yi)
    zd, zi = mod(zi, 1), unsafe_trunc(Int, zi)
    ##### return fraction of grid spacing
    return normalized_trilinear_interpolation_v(v, xd, yd, zd, xi+g.Hx, yi+g.Hy, zi+g.Hz, g)
end

@inline function w_itpl(w, x, y, z, ac, g::AbstractGrid)
    xi = get_xc_index(x) * ac
    yi = get_yc_index(y) * ac
    zi = get_zf_index(z) * ac
    xd, xi = mod(xi, 1), unsafe_trunc(Int, xi)
    yd, yi = mod(yi, 1), unsafe_trunc(Int, yi)
    zd, zi = mod(zi, 1), unsafe_trunc(Int, zi)
    ##### return fraction of grid spacing
    return normalized_trilinear_interpolation_w(w, xd, yd, zd, xi+g.Hx, yi+g.Hy, zi+g.Hz, g)
end


text:
# Adaptation of some examples in MosekTools
# Refer to:
# https://github.com/JuliaOpt/MosekTools.jl/blob/master/test/jump_soc.jl
#
# MIT License
# Copyright (c) 2017, Ulf Worsøe, Mosek Aps

using KNITRO
using JuMP
using Test
using MathOptInterface
const MOI = MathOptInterface

model = Model(with_optimizer(KNITRO.Optimizer, outlev=0))
@variable(model, x)
@variable(model, y)
@variable(model, t >= 0.0)
@objective(model, Min, t)
@constraint(model, x + y >= 1.0)
@constraint(model, [t, x, y] in MOI.SecondOrderCone(3))

JuMP.optimize!(model)

@test JuMP.has_values(model)
@test JuMP.termination_status(model) == MOI.LOCALLY_SOLVED
@test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
@test JuMP.objective_value(model) ≈ sqrt(1/2) atol=1e-6
@test JuMP.value.([x,y,t]) ≈ [0.5, 0.5, sqrt(1/2)] atol=1e-3


text:
using FFTW
using Distributions
using Random

function fftFilter2D(maxlevel, h)
    #Random.seed!(2342)
    L = 2^maxlevel
    L2 = L ÷ 2
    β = (h+1)/2
    f = zeros(Complex{Float64}, L,L)

    # upper left  and lower right
    rad = [(i^2 + j^2)^(-β) for i=1:L2, j=1:L2]
    rad = rad .* rand(Normal(), L2, L2)
    ϕ = 2π .* rand(L2, L2)
    exp_iϕ = cos.(ϕ) .+ sin.(ϕ) .* im
    Z = rad.*exp_iϕ
    Z[L2,L2] = real(Z[L2, L2])
    Zt = conj(Z[L2:-1:1, L2:-1:1])

    f[2:L2+1,2:L2+1] = Z     # line 1 is dealt with separetely
    f[(L2+1):L,(L2+1):L] = Zt

    # lower left and  upper right
    rad = [(i^2 + j^2)^(-β) for i=1:(L2-1), j=1:(L2-1)]
    rad = rad .* rand(Normal(), (L2-1), (L2-1))
    ϕ = 2π .* rand((L2-1), (L2-1))
    exp_iϕ = cos.(ϕ) .+ sin.(ϕ) .* im
    W = rad.*exp_iϕ
    Wt = conj(W[(L2-1):-1:1, (L2-1):-1:1])

    f[L:-1:(L2+2),2:L2] = W
    f[L2:-1:2,(L2+2):L] = Wt

    # first line
    rad = [j^(-2*β) for j=1:L2]
    rad = rad .* rand(Normal(), L2)
    ϕ = (2*π) .* rand(L2)
    exp_iϕ = cos.(ϕ) .+ sin.(ϕ) .* im
    line = rad.*exp_iϕ
    line[L2] = real(line[L2])
    linet = conj(line[L2:-1:1])

    f[1,2:(L2+1)] = line
    f[1,(L2+1):L] = linet

    # first column
    rad = [i^(-2*β) for i=1:L2]
    rad = rad .* rand(Normal(), L2)
    ϕ = (2*π) .* rand(L2)
    exp_iϕ = cos.(ϕ) .+ sin.(ϕ) .* im
    column = rad.*exp_iϕ
    column[L2] = real(column[L2])
    columnt = conj(column[L2:-1:1])

    f[2:(L2+1),1] = column
    f[(L2+1):L,1] = columnt

    f[1,1] = 1

    g = ifft(f)
    # println(sum(abs.(imag.(g))))
    # h = fftshift(f)
    # display( heatmap(log.(abs.(h.*h)).+ 1.0,size=(600,500), color=ColorGradient([:black,:yellow,:red])))
    z = real(g)
    z = (z .- mean(z)).*L^2
    return z
end

include("dfa.jl")

function fbmSurf(maxlevel, h; atol=0.01)
    """
    - uses fftFilter2D to generate surfaces
    - checks Hurst using DFA
    - returns surface heights z with Hurst within atol
    """
    hurst = 50.0
    n=0
    z = []
    while !isapprox(hurst, h; atol=atol) && n<1e3
        z = fftFilter2D(maxlevel, h)
        hurst = dfa2d(z; scales=3:0.5:maxlevel-2)
        n+=1
    end
    println(n)
    return z
end

function make2DSurfs(;maxlevel=9, hursts=0.3:0.1:0.8, samples=1:1, fbm=nothing, atol=0.01)
    """
    returns a dictionary with the set of selfAffine surfaces
    """

    fbm = fbm == nothing ? Dict() : fbm
    fbm["hursts"]=collect(hursts)
    fbm["samples"]=collect(samples)
    fbm["maxlevel"] = maxlevel
    for h=hursts
        sample = Dict()
        for s=samples
            z = fbmSurf(maxlevel,h, atol=atol)
            sample["$s"]=z
        end
        fbm["$h"]=sample
        println("finished $h")
    end
    return fbm
end

function make2DSurfs(filename::String)
    fbm = load(filename)["fbm"]
end


text:
using Base.Test
ismatch(1:10, isodd, iseven, isodd, matchall=true) == [
	([1], [2], [3]),
	([5], [6], [7]),
]
ismatch("abcffedefgx", x-> x in ('d', 'e', 'f')) == (Any['f','f','e','d','e','f'],) ||
error("second failed")

_replace([1,2,3,4,5,6],
	(isodd,) => x -> 77
)

function test(a, b)
	if a == 10
		return b
	else
		return 77
	end
end

ast = Base.uncompressed_ast(code_lowered(test, (Int, Int))[])
isunless(x::Expr) = x.head == :gotoifnot
isunless(x) = false
islabelnode(x::LabelNode) = true
islabelnode(x) = false
function remove_goto(ast)
	_replace(ast, isunless, (x->true)*(*), islabelnode) do unless, body, label
		condition = unless[].args[1]
		Expr(:if, condition, remove_goto(body)...)
	end
end


text:
module TestMaxdiff
    using Calculus2
    using Base.Test

    @test isequal(
        Calculus2.maxdiff(5, 7),
        2
    )

    @test isequal(
        Calculus2.maxdiff([0, 0, 3, 3], [+1, -1, +1, -1]),
        4
    )
end


text:
using Documenter
using Dagitty

makedocs(
    sitename="Dagitty",
    modules=[Dagitty],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
)

deploydocs(
    repo = "github.com/Shmuma/Dagitty.jl.git",
    devbranch = "main",
)


text:
function stencil1d(m::Int)

    ic = Int(ceil(m/2.0))
    M = zeros((m, m))
    for i=1:m
        for j=1:m
            M[i, j] = (i - ic)^(j-1) / factorial(j-1)
        end
    end
    return inv(M)
end

function periodic_projection(N, m)

    js = [N-m+1:N...]
    append!(js, [1:N...], [1:m...])

    P = sparse([1:N...], [m+1:m+N...], ones(N), N, N+2*m)
    PT = sparse([1:N+2*m...], js, ones(N+2*m), N+2*m, N)

    return P, PT
end

function periodic_derivative_1d(n, m, N)

    Nb = N + 2*m
    D = sparse(zeros(Float64, (Nb, Nb)))
    stencil = stencil1d(m)[n+1,:]
    ic = Int(ceil(m/2.0))

    for i=1:m
        k = i - ic
        vec = stencil[i] * ones(Nb - abs(k))
        D += spdiagm(k => vec)
    end

    P, PT = periodic_projection(N, m)
    return P * D * PT
end

function periodic_derivative(nx, ny, mx, my, Nx, Ny, dx, dy)

    Dx = periodic_derivative_1d(nx, mx, Nx) / dx^nx
    Dy = periodic_derivative_1d(ny, my, Ny) / dy^ny
    return kron(Dy, Dx)
end