text:
function maxpool_nnpack!(y::A, x::A, pdims::PoolDims) where {A<:Array{Float32, 4}}
    check_dims(size(x), size(y), pdims)
    threadpool = select_threadpool(pdims, size(y, 4))
    nnp_max_pooling_output(y, x, kernel_size(pdims), padding = padding(pdims),
                           stride = stride(pdims), threadpool = threadpool)
end

function conv_nnpack!(y::A1, x::A1, w::A1, cdims::ConvDims;
                                       b::A2 = zeros(Float32, size(x, 3)),
                                       algo = UInt32(0)) where {A1<:Array{Float32, 4},
                                                                A2<:Array{Float32, 1}}
    check_dims(size(x), size(w), size(y), cdims)
    threadpool = select_threadpool(cdims, size(y, 4))

    if flipkernel(cdims) == 0
        w = flipweight(w)
    end

    nnp_convolution_output(y, x, w, b, algo = algo, padding = padding(cdims),
                           stride = stride(cdims), threadpool = threadpool)
end

function ∇conv_data_nnpack!(dx::A, dy::A, w::A, cdims::ConvDims;
                                             algo = UInt32(0)) where{A<:Array{Float32, 4}}
    check_dims(size(dx), size(w), size(dy), cdims)
    threadpool = select_threadpool(cdims, size(dy, 4))

    if flipkernel(cdims) == 0
        w = flipweight(w)
    end

    nnp_convolution_input_gradient(dx, dy, w, algo = algo, padding = padding(cdims),
                                   stride = stride(cdims), threadpool = threadpool)
end

function ∇conv_filter_nnpack!(dw::A, x::A, dy::A, cdims::ConvDims;
                                               algo = UInt32(0)) where{A<:Array{Float32, 4}}
    check_dims(size(x), size(dw), size(dy), cdims)
    threadpool = select_threadpool(cdims, size(dy, 4))

    nnp_convolution_kernel_gradient(dw, x, dy, algo = algo, padding = padding(cdims),
                                    stride = stride(cdims), threadpool = threadpool)

    if flipkernel(cdims) == 0
        dw .= flipweight(dw)
    end

    dw
end


text:
############################################################################
### Types
############################################################################

struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V, V2, F4, A4, F5, A5, D1} <: SemImplySymbolic
    Σ_function::F1
    ∇Σ_function::F2
    ∇²Σ_function::F3
    Σ::A1
    ∇Σ::A2
    ∇²Σ::A3
    Σ_symbolic::S1
    ∇Σ_symbolic::S2
    ∇²Σ_symbolic::S3
    n_par::V
    ram_matrices::V2
    μ_function::F4
    μ::A4
    ∇μ_function::F5
    ∇μ::A5
    identifier::D1
end

############################################################################
### Constructors
############################################################################

function RAMSymbolic(;
        specification,
        loss_types = nothing,
        vech = false,
        gradient = true,
        hessian = false,
        kwargs...)

    if specification isa RAMMatrices
        ram_matrices = specification
        identifier = StructuralEquationModels.identifier(ram_matrices)
    elseif specification isa ParameterTable
        ram_matrices = RAMMatrices(specification)
        identifier = StructuralEquationModels.identifier(ram_matrices)
    else
        throw(ErrorException("The RAMSymbolic constructor does not know how to handle your specification object.
        \n Please specify your model as either a ParameterTable or RAMMatrices."))
    end

    n_par = length(ram_matrices.parameters)
    n_var, n_nod = ram_matrices.size_F

    par = (Symbolics.@variables θ[1:n_par])[1]

    A = zeros(Num, n_nod, n_nod)
    S = zeros(Num, n_nod, n_nod)
    !isnothing(ram_matrices.M_ind) ? M = zeros(Num, n_nod) : M = nothing
    F = zeros(ram_matrices.size_F); F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0

    set_RAMConstants!(A, S, M, ram_matrices.constants)
    fill_A_S_M(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par)

    A, S, F = sparse(A), sparse(S), sparse(F)

    if !isnothing(loss_types)
        any(loss_types .<: SemWLS) ? vech = true : nothing
    end

    # Σ
    Σ_symbolic = get_Σ_symbolic_RAM(S, A, F; vech = vech)
    #print(Symbolics.build_function(Σ_symbolic)[2])
    Σ_function = Symbolics.build_function(Σ_symbolic, par, expression=Val{false})[2]
    Σ = zeros(size(Σ_symbolic))
    precompile(Σ_function, (typeof(Σ), Vector{Float64}))

    # ∇Σ
    if gradient
        ∇Σ_symbolic = Symbolics.sparsejacobian(vec(Σ_symbolic), [par...])
        ∇Σ_function = Symbolics.build_function(∇Σ_symbolic, par, expression=Val{false})[2]
        constr = findnz(∇Σ_symbolic)
        ∇Σ = sparse(constr[1], constr[2], fill(1.0, nnz(∇Σ_symbolic)), size(∇Σ_symbolic)...)
        precompile(∇Σ_function, (typeof(∇Σ), Vector{Float64}))
    else
        ∇Σ_symbolic = nothing
        ∇Σ_function = nothing
        ∇Σ = nothing
    end

    if hessian
        n_sig = length(Σ_symbolic)
        n_par = size(par, 1)
        ∇²Σ_symbolic_vec = [Symbolics.sparsehessian(σᵢ, [par...]) for σᵢ in vec(Σ_symbolic)]

        @variables J[1:n_sig]
        ∇²Σ_symbolic = zeros(Num, n_par, n_par)
        for i in 1:n_sig
            ∇²Σ_symbolic += J[i]*∇²Σ_symbolic_vec[i]
        end

        ∇²Σ_function = Symbolics.build_function(∇²Σ_symbolic, J, par, expression=Val{false})[2]
        ∇²Σ = zeros(n_par, n_par)
    else
        ∇²Σ_symbolic = nothing
        ∇²Σ_function = nothing
        ∇²Σ = nothing
    end

    # μ
    if !isnothing(M)
        μ_symbolic = get_μ_symbolic_RAM(M, A, F)
        μ_function = Symbolics.build_function(μ_symbolic, par, expression=Val{false})[2]
        μ = zeros(size(μ_symbolic))
        if gradient
            ∇μ_symbolic = Symbolics.jacobian(μ_symbolic, [par...])
            ∇μ_function = Symbolics.build_function(∇μ_symbolic, par, expression=Val{false})[2]
            ∇μ = zeros(size(F, 1), size(par, 1))
        else
            ∇μ_function = nothing
            ∇μ = nothing
        end
    else
        μ_function = nothing
        μ = nothing
        ∇μ_function = nothing
        ∇μ = nothing
    end

    return RAMSymbolic(
        Σ_function,
        ∇Σ_function,
        ∇²Σ_function,
        Σ,
        ∇Σ,
        ∇²Σ,
        Σ_symbolic,
        ∇Σ_symbolic,
        ∇²Σ_symbolic,
        n_par,
        ram_matrices,
        μ_function,
        μ,
        ∇μ_function,
        ∇μ,
        identifier
    )
end

############################################################################
### functors
############################################################################

function (imply::RAMSymbolic)(par, F, G, H, model)
    imply.Σ_function(imply.Σ, par)
    if G || H
        imply.∇Σ_function(imply.∇Σ, par)
    end
    if !isnothing(imply.μ)
        imply.μ_function(imply.μ, par)
        if G || H
            imply.∇μ_function(imply.∇μ, par)
        end
    end
end

############################################################################
### Recommended methods
############################################################################

identifier(imply::RAMSymbolic) = imply.identifier
n_par(imply::RAMSymbolic) = imply.n_par

function update_observed(imply::RAMSymbolic, observed::SemObs; kwargs...)
    if n_man(observed) == size(imply.Σ)
        return imply
    else
        return RAMSymbolic(;observed = observed, kwargs...)
    end
end

############################################################################
### additional methods
############################################################################

Σ(imply::RAMSymbolic) = imply.Σ
∇Σ(imply::RAMSymbolic) = imply.∇Σ
∇²Σ(imply::RAMSymbolic) = imply.∇²Σ

μ(imply::RAMSymbolic) = imply.μ
∇μ(imply::RAMSymbolic) = imply.∇μ

Σ_function(imply::RAMSymbolic) = imply.Σ_function
∇Σ_function(imply::RAMSymbolic) = imply.∇Σ_function
∇²Σ_function(imply::RAMSymbolic) = imply.∇²Σ_function

############################################################################
### additional functions
############################################################################

function get_Σ_symbolic_RAM(S, A, F; vech = false)
    invia = neumann_series(A)
    Σ_symbolic = F*invia*S*permutedims(invia)*permutedims(F)
    Σ_symbolic = Array(Σ_symbolic)
    # Σ_symbolic = Symbolics.simplify.(Σ_symbolic)
    Threads.@threads for i in eachindex(Σ_symbolic)
        Σ_symbolic[i] = Symbolics.simplify(Σ_symbolic[i])
    end
    if vech Σ_symbolic = Σ_symbolic[tril(trues(size(F, 1), size(F, 1)))] end
    return Σ_symbolic
end

function get_μ_symbolic_RAM(M, A, F)
    invia = neumann_series(A)
    μ_symbolic = F*invia*M
    μ_symbolic = Array(μ_symbolic)
    Threads.@threads for i in eachindex(μ_symbolic)
        μ_symbolic[i] = Symbolics.simplify(μ_symbolic[i])
    end
    return μ_symbolic
end

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::RAMSymbolic)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end

ext: jl


text:
# module Bukdu

immutable Route{AC<:ApplicationController}
    verb::Function # get
    kind::Symbol # :match, :forward
    path::String
    host::String
    controller::Type{AC}
    action::Function
    private::Dict{Symbol,String}
    assigns::Dict{Symbol,String}
end

module RouterRoute

import ..Bukdu: ApplicationController
import ..Bukdu: Route

routes = Vector{Route}()

function build{AC<:ApplicationController}(kind::Symbol, verb::Function, path::String,
             host::String, ::Type{AC}, action::Function, private::Dict{Symbol,String}, assigns::Dict{Symbol,String})
    Route(verb, kind, path, host, AC, action, private, assigns)
end

end # module Bukdu.RouterRoute


text:
# include("../setup.jl")
include("../all_functions.jl")

# process for se_kernel_base
# @vars σ δ λ
# kernel_coder(σ * exp(-δ * δ / (2 * λ * λ)), "se")
@vars δ λ
kernel_coder(se_kernel_base(λ, δ), "se")

#process for matern52_kernel_base
@vars δ λ
kernel_coder(matern52_kernel_base(λ, δ), "m52")

@vars δ, λ
kernel_coder(pp_kernel_base(λ, δ), "pp"; cutoff_var="λ")

@vars δ λ1 λ2 sratio
kernel_coder(matern52_kernel_base(λ1, δ) + sratio * sratio * matern52_kernel_base(λ2, δ), "m52_m52")

@vars δ λ1 λ2 sratio
kernel_coder(se_kernel_base(λ1, δ) + sratio * sratio * se_kernel_base(λ2, δ), "se_se")

# process for periodic_kernel_base
# @vars δ se_P λ
# kernel_coder(periodic_kernel_base([se_P, λ], δ), "periodic")
# @vars δ λ
# kernel_coder(se_kernel_base(λ, δ), "se"; periodic_var="δ")

# process for quasi_periodic_kernel_base
# @vars δ se_λ qp_P p_λ
# kernel_coder(quasi_periodic_kernel_base([se_λ, qp_P, p_λ], δ), "quasi_periodic")
@vars δ δp se_λ p_amp
kernel_coder(se_kernel_base(se_λ, δ) * se_kernel_base(1 / p_amp, δp), "qp"; periodic_var="δp")

#process for rq_kernel_base

# @vars δ k θ
# kernel_coder(rq_kernel_base([k, θ], δ), "rq_kernel")
# @vars δ α β
# kernel_coder(rq_kernel_base([α, β], δ), "rq_kernel")
@vars δ α μ
kernel_coder(rq_kernel_base([α, μ], δ), "rq")

@vars δ α μ
kernel_coder(rm52_kernel_base([α, μ], δ), "rm52")

# @vars δ P λ α
# kernel_coder(periodic_rq_kernel_base([P, λ, α], δ), "periodic_rq_kernel")
@vars δ λ α
kernel_coder(rq_kernel_base([α, λ], δ), "rq"; periodic_var="δ")

@vars δ λ
kernel_coder(cosine_kernel_base(λ, δ), "cos")


text:
using StaticArrays

export billiard_rectangle, billiard_sinai, billiard_polygon, billiard_lorentz,
billiard_raysplitting_showcase, billiard_hexagonal_sinai, billiard_bunimovich,
billiard_stadium, billiard_mushroom, billiard_logo, billiard_vertices,
polygon_vertices

####################################################
## Famous/Standard Billiards
####################################################
"""
    billiard_rectangle(x=1.0, y=1.0; setting = "standard")
Return a vector of obstacles that defines a rectangle billiard of size (`x`, `y`).

### Settings
* "standard" : Specular reflection occurs during collision.
* "periodic" : The walls are `PeriodicWall` type,
  enforcing periodicity at the boundaries
* "random" : The velocity is randomized upon collision.
* "ray-splitting" : All obstacles in the billiard allow for ray-splitting.
"""
function billiard_rectangle(x′ = 1.0, y′ = 1.0;
    x = x′, y = y′, setting::String = "standard")

    x = convert(AbstractFloat, x)
    x, y = promote(x,y)
    o = typeof(x)(0.0)
    if setting == "standard"
        sp = [o,y]; ep = [o, o]; n = [x,o]
        leftw = InfiniteWall(sp, ep, n, "Left wall")
        sp = [x,o]; ep = [x, y]; n = [-x,o]
        rightw = InfiniteWall(sp, ep, n, "Right wall")
        sp = [x,y]; ep = [o, y]; n = [o,-y]
        topw = InfiniteWall(sp, ep, n, "Top wall")
        sp = [o,o]; ep = [x, o]; n = [o,y]
        botw = InfiniteWall(sp, ep, n, "Bottom wall")
    elseif setting == "periodic"
        sp = [o,y]; ep = [o, o]; n = [x,o]
        leftw = PeriodicWall(sp, ep, n, "Left periodic boundary")
        sp = [x,o]; ep = [x, y]; n = [-x,o]
        rightw = PeriodicWall(sp, ep, n, "Right periodic boundary")
        sp = [x,y]; ep = [o, y]; n = [o,-y]
        topw = PeriodicWall(sp, ep, n, "Top periodic boundary")
        sp = [o,o]; ep = [x, o]; n = [o,y]
        botw = PeriodicWall(sp, ep, n, "Bottom periodic boundary")
    elseif setting == "random"
        sp = [o,y]; ep = [o, o]; n = [x,o]
        leftw = RandomWall(sp, ep, n, "Left random wall")
        sp = [x,o]; ep = [x, y]; n = [-x,o]
        rightw = RandomWall(sp, ep, n, "Right random wall")
        sp = [x,y]; ep = [o, y]; n = [o,-y]
        topw = RandomWall(sp, ep, n, "Top random wall")
        sp = [o,o]; ep = [x, o]; n = [o,y]
        botw = RandomWall(sp, ep, n, "Bottom random wall")
    elseif setting == "ray-splitting"
        sp = [o,y]; ep = [o, o]; n = [x,o]
        leftw = SplitterWall(sp, ep, n, "Left ray-splitting wall")
        sp = [x,o]; ep = [x, y]; n = [-x,o]
        rightw = SplitterWall(sp, ep, n, "Right ray-splitting wall")
        sp = [x,y]; ep = [o, y]; n = [o,-y]
        topw = SplitterWall(sp, ep, n, "Top ray-splitting wall")
        sp = [o,o]; ep = [x, o]; n = [o,y]
        botw = SplitterWall(sp, ep, n, "Bottom ray-splitting wall")
    else
        throw(ArgumentError("The given setting=$setting is unknown."))
    end
    return Billiard(botw, rightw, topw, leftw)
end

"""
    billiard_sinai(r=0.25, x=1.0, y=1.0; setting = "standard")
Return a vector of obstacles that defines a Sinai billiard of size (`x`, `y`) with
a disk in its center, of radius `r`.

In the periodic case, the system is also known as "Lorentz Gas".

### Settings
* "standard" : Specular reflection occurs during collision.
* "periodic" : The walls are `PeriodicWall` type,
  enforcing periodicity at the boundaries
* "random" : The velocity is randomized upon collision.
* "ray-splitting" : All obstacles in the billiard allow for ray-splitting.
"""
function billiard_sinai(r′ = 0.25, x′ = 1.0, y′ = 1.0;
    r = r′, x = x′, y = y′, setting = "standard")

    if (setting == "periodic") && (r>=x/2 || r>=y/2)
        es = "Disk radius too big for a periodic Sinai billiard.\n"
        es*= "Obstacles must not overlap with `PeriodicWall`s."
        error(es)
    end
    r, x, y = promote(r,x,y)
    bdr = billiard_rectangle(x, y; setting = setting)

    c = [x/2, y/2]
    if setting == "random"
        centerdisk = RandomDisk(c, r, "Random disk")
    elseif setting == "ray-splitting"
        centerdisk = Antidot(c, r, "Antidot")
    else
        centerdisk = Disk(c, r, "Disk")
    end

    return Billiard(centerdisk, bdr...)
end

"""
    billiard_lorentz(r=0.25, x=1.0, y=1.0)
Alias for `billiard_sinai(r,x,y; setting = "periodic")`.
"""
billiard_lorentz(r=0.25, x=1.0, y=1.0) = billiard_sinai(r,x,y; setting = "periodic")

"""
    billiard_polygon(n::Int, R, center = [0,0]; setting = "standard")
Return a vector of obstacles that defines a regular-polygonal billiard
with `n` sides, radius `r` and given `center`.

Note: `R` denotes the so-called outer radius, not the inner one.

### Settings
* "standard" : Specular reflection occurs during collision.
* "periodic" : The walls are `PeriodicWall` type, enforcing periodicity
  at the boundaries. Only available for `n=4` or `n=6`.
* "random" : The velocity is randomized upon collision.
"""
function billiard_polygon(sides′::Int, r′::Real, center′ = [0,0];
    sides::Int = sides′, r = r′, center = center′, setting = "standard")

    S = typeof(convert(AbstractFloat, r))
    bd = Obstacle{S}[]
    verteces = polygon_vertices(r, sides, center)

    if setting == "standard"
        T = InfiniteWall
        wallname = "wall"
    elseif setting == "periodic"
        if sides != 4 && sides != 6
            error("Polygonal and periodic billiard can exist only for `n=4` or `n=6`")
        end
        T = PeriodicWall
        wallname = "periodic wall"
        inr = S(r*cos(π/sides))
    elseif setting == "random"
        T = RandomWall
        wallname = "random wall"
    end

    for i in eachindex(verteces)
        starting = verteces[i]
        ending = verteces[mod1(i+1, sides)]
        # Normal vector must look at where the particle is coming from
        w = ending - starting
        if setting == "periodic"
            normal = 2inr*normalize([-w[2], w[1]])
            wall = T(starting, ending, normal, wallname*" $i")
        else
            normal = [-w[2], w[1]]
            wall = T(starting, ending, normal, wallname*" $i")
        end
        push!(bd, wall)
    end
    return Billiard(bd)
end

"""
    polygon_vertices(r, sides, center = [0, 0], θ=0) -> v
Return the vertices that define a regular polygon with `sides` sides and
radius `r`, centered at `center`.

Optionally rotate by `θ` degrees.
"""
function polygon_vertices(r, sides, center = [0, 0.0], θ = 0)
    S = typeof(convert(AbstractFloat, r))
    [S[r*cos(θ + 2π*i/sides), r*sin(θ + 2π*i/sides)] .+ center for i in 1:sides]
end

"""
    billiard_hexagonal_sinai(r, R, center = [0,0]; setting = "standard")
Create a sinai-like billiard, which is a hexagon of outer radius `R`, containing
at its center (given by `center`) a disk of radius `r`. The `setting` keyword
is passed to `billiard_polygon`.
"""
function billiard_hexagonal_sinai(r′::Real = 0.5, R′::Real = 1.0, center′ = [0,0];
    r = r′, R = R′, center = center′, setting = "standard")
    r, R = promote(r, R)
    T = typeof(r); center = T[center...]
    bdr = billiard_polygon(6, R, center; setting = setting)
    DT = setting == "random" ? RandomDisk : Disk
    return Billiard(Disk(center, r), bdr...)
end

"""
    billiard_raysplitting_showcase(x=2.0, y=1.0, r1=0.3, r2=0.2) -> bd, rayspl
Showcase example billiard for ray-splitting processes. A rectangle `(x,y)` with a
SplitterWall at `x/2` and two disks at each side, with respective radii `r1`, `r2`.

**Notice**: This function returns a billiard `bd` as well as a `rayspl`
dictionary!
"""
function billiard_raysplitting_showcase(x=2.0, y=1.0, r1=0.3, r2=0.2)

    r1≥x/4 || r2≥x/4 && throw(ArgumentError("Disks overlap with walls!"))

    bdr =  billiard_rectangle(x, y)
    sw = SplitterWall([x/2, 0.0], [x/2,y], [-1,0], true)
    a1 = Antidot([x/4, y/2], r1, "Left Antidot")
    a2 = Antidot([3x/4, y/2], r2, "Right Antidot")

    sa = (φ, pflag, ω) -> pflag ? 0.5φ : 2.0φ
    Tp = (p) -> (φ, pflag, ω) -> begin
        if pflag
            p*exp(-(φ)^2/2(π/8)^2)
        else
            abs(φ) < π/4 ? (1-p)*exp(-(φ)^2/2(π/4)^2) : 0.0
        end
    end
    newoantidot = ((x, bool) -> bool ? -0.5x : -2.0x)
    newowall = ((x, bool) -> bool ? 0.5x : 2.0x)

    # I Need 3 ray-splitters because I want different affect for
    # different antidots
    raywall = RaySplitter([3], Tp(0.5), sa, newowall)
    raya = RaySplitter([1, 2], Tp(0.64), sa, newoantidot)

    return Billiard(a1, a2, sw, bdr...), (raywall, raya)
end

"""
    billiard_mushroom(sl = 1.0, sw = 0.2, cr = 1.0, sloc = 0.0; door = true)
Create a mushroom billiard with stem length `sl`, stem width `sw` and
cap radius `cr`. The center of the cap (which is Semicircle) is always
at `[0, sl]`. The center of the stem is located at `sloc`.

Optionally, the bottom-most `Wall` is a `Door` (see [`escapetime`](@ref)).
"""
function billiard_mushroom(stem_length = 1.0, stem_width=0.2, cap_radious=1.0,
    stem_location = 0.0; sl = stem_length, sw = stem_width, cr = cap_radious,
    sloc = stem_location, door = true)

    abs(sloc) + sw/2 > cr && error("Stem is outside the mushroom cap!")

    sl, sw, cr, sloc = promote(sl, sw, cr, sloc)
    T = eltype(sl)
    leftcorn = SV(-sw/2 + sloc, 0)
    rightcorn = SV(sw/2 + sloc, 0)
    upleftcorn = SV(-sw/2 + sloc, sl)
    uprightcorn = SV(sw/2 + sloc, sl)

    stembot = FiniteWall(leftcorn, rightcorn, SV(0, sw), door, "Stem bottom")
    stemleft = FiniteWall(upleftcorn, leftcorn, SV(sw, 0), false, "Stem left")
    stemright = FiniteWall(rightcorn, uprightcorn, SV(-sw, 0), false, "Stem right")

    farleft = SV(-cr, sl)
    farright = SV(cr, sl)

    capbotleft = FiniteWall(
    farleft, upleftcorn, SV(0, sw), false, "Cap bottom left")
    capbotright = FiniteWall(
    uprightcorn, farright, SV(0, sw), false, "Cap bottom right")

    cap = Semicircle([0, sl], cr, [0, -T(1.0)], "Mushroom cap")
    return Billiard(stembot, stemright, capbotright, cap, capbotleft, stemleft)
end

"""
    billiard_bunimovich(l=1.0, w=1.0)
Return a vector of `Obstacle`s that define a Buminovich billiard, also called a
stadium. The length is considered *without* the attached semicircles, meaning that the
full length of the billiard is `l + w`. The left and right edges of the stadium
are [`Semicircle`](@ref)s.

`billiard_stadium` is an alias of `billiard_bunimovich`.
"""
function billiard_bunimovich(l′=1.0, w′=1.0; l = l′, w = w′)

    l = convert(AbstractFloat, l)
    l, w = promote(l,w)
    o = typeof(l)(0.0)
    bw = InfiniteWall([o,o], [l,o], [o,  w], "Bottom wall")
    tw = InfiniteWall([l,w], [o,w], [o, -w], "Top wall")
    leftc = Semicircle([o, w/2], w/2, [l, o], "Left semicircle")
    rightc = Semicircle([l, w/2], w/2, [-l, o], "Right semicircle")
    return Billiard(bw, rightc, tw, leftc)
end

billiard_stadium = billiard_bunimovich

"""
    billiard_logo(;h=1.0, α=0.8, r=0.18, off=0.25, T = Float64) -> bd, ray
Create the billiard used as logo of `DynamicalBilliards` and return it
along with the tuple of raysplitters.
"""
function billiard_logo(;h=1.0, α=0.8, r=0.18, off=0.25, T = Float64)

    cos6 = cos(π/6)
    β = (h - cos6*α)/cos6
    t = α + 2β
    (α, r, h, off, cos6, β, t) = T.((α, r, h, off, cos6, β, t))
    center_of_mass = T[0.0, √3*t/6]
    startloc = T[-α/2, 0.0]

    # create directions of the hexagonal 6:
    hexvert = [T.((cos(2π*i/6), sin(2π*i/6))) for i in 1:6]
    dirs = [SVector{2}(hexvert[i] .- hexvert[mod1(i+1, 6)]) for i in 1:6]

    frame = Obstacle{T}[]

    sp = startloc
    ep = startloc + α*dirs[1]
    normal = (w = ep .- sp; [-w[2], w[1]])
    push!(frame, InfiniteWall(sp, ep, normal, "frame 1"))

    for i in 2:6
        s = iseven(i) ? β : α
        WT = InfiniteWall #iseven(i) ? RandomWall : InfiniteWall
        sp = frame[i-1].ep
        ep = sp + s*dirs[i]
        normal = (w = ep .- sp; [-w[2], w[1]])
        push!(frame, WT(sp, ep, normal, "frame $(i)"))
    end

    # Radii of circles that compose the Julia logo
    offset = T[0.0, off]

    R = T[cos(2π/3) -sin(2π/3);
         sin(2π/3)  cos(2π/3)]

    green = Disk(center_of_mass .+ offset, r, "green")
    red = Antidot(center_of_mass .+ R*offset, r, "red")
    purple = RandomDisk(center_of_mass .+ R*R*offset, r, "purple")

    bd = Billiard(green, red, purple, frame...)

    # %% Raysplitting functions for the red circle:
    refraction = (φ, pflag, ω) -> pflag ? 0.5φ : 2.0φ
    transmission_p = (p) -> (φ, pflag, ω) -> begin
        if pflag
            p*exp(-(φ)^2/2(π/8)^2)
        else
            abs(φ) < π/4 ? (1-p)*exp(-(φ)^2/2(π/4)^2) : 0.0
        end
    end
    newoantidot = ((x, bool) -> bool ? -2.0x : -0.5x)
    raya = RaySplitter([2], transmission_p(0.8), refraction, newoantidot)
    return bd, raya
end

export billiard_iris

"""
    billiard_iris(a=0.2, b=0.4, w=1.0; setting = "standard")
Return a billiard that is a square of side `w` enclosing at its center an ellipse
with semi axes `a`, `b`.
"""
function billiard_iris(a′ = 0.2, b′ = 0.4, w′ = 1.0;
    setting = "standard", a = a′, b = b′, w = w′)

    rec = billiard_rectangle(w, w; setting = setting)
    e = Ellipse([w/2, w/2], a, b, true, "Ellipse")
    return Billiard(e, rec.obstacles...)
end

"""
    billiard_vertices(v, type = FiniteWall)
Construct a polygon billiard that connects the given vertices `v`
(vector of 2-vectors). The vertices should construct a billiard in a
counter-clockwise orientation (i.e. the normal vector always points to the
left of `v[i+1] - v[i]`.).

`type` decides what kind of walls to use.
Ths function assumes that the first entry of `v` should be connected with the
last.
"""
function billiard_vertices(vertices, type = FiniteWall)
    @assert vertices[end] != vertices[1]
    obs = type[]
    for i in 1:length(vertices)-1
        sp, ep = vertices[i], vertices[i+1]
        push!(obs, type(sp, ep))
    end
    sp, ep = vertices[end], vertices[1]
    push!(obs, type(sp, ep))
    return Billiard(obs)
end


text:
module TestWorkStealingDeque

using ConcurrentCollections
using Test

function test_single_thread_push_pop()
    deque = WorkStealingDeque{Int}()
    xs = 1:50
    foldl(push!, xs; init = deque)
    @test [pop!(deque) for _ in xs] == reverse(xs)
    @test maybepop!(deque) === nothing
    @test maybepopfirst!(deque) === nothing
    push!(deque, 1)
    @test maybepopfirst!(deque) === Some(1)
    @test maybepopfirst!(deque) === nothing
    @test maybepop!(deque) === nothing

    foldl(push!, xs; init = deque)
    n = length(deque.buffer)
    sizehint!(deque, 2 * n)
    @test length(deque.buffer) > n
    @test [pop!(deque) for _ in xs] == reverse(xs)
end

function random_pushpop(xs; ntasks = Threads.nthreads() - 1, sentinel = -1)
    sentinel = convert(eltype(xs), sentinel)  # verify argument
    ntasks = max(1, ntasks)

    deque = WorkStealingDeque{eltype(xs)}()

    local tasks, zs
    try
        tasks = map(1:ntasks) do _
            Threads.@spawn begin
                local ys = eltype(xs)[]
                while true
                    local r = maybepopfirst!(deque)
                    if r === nothing
                        GC.safepoint()
                        continue
                    end
                    local y = something(r)
                    y == $sentinel && break
                    push!(ys, y)
                end
                ys
            end
        end

        zs = eltype(xs)[]
        for (i, x) in enumerate(xs)
            push!(deque, x)
            # continue
            if mod(i, 8) == 0
                r = maybepop!(deque)
                GC.safepoint()
                r === nothing && continue
                push!(zs, something(r))
            end
        end
    finally
        for _ in 1:ntasks
            push!(deque, sentinel)
        end
    end

    return zs, fetch.(tasks)
end

function test_random_push_pop()
    @testset for trial in 1:100
        @testset for T in [Int, Any, Pair{Any,Int}]
            test_random_push_pop(T)
        end
    end
end

function test_random_push_pop(T::Type)
    sentinel = -1
    xs = 1:2^10
    return test_random_push_pop(T::Type, xs, sentinel)
end

function test_random_push_pop(T::Type{Pair{Any,Int}})
    sentinel = T(-1, -1)
    xs = [T(x, x) for x in 1:2^10]
    return test_random_push_pop(T::Type, xs, sentinel)
end

function test_random_push_pop(T::Type, xs, sentinel)
    if T !== eltype(xs)
        xs = collect(T, xs)
    end
    zs, yss = random_pushpop(xs; sentinel)
    @test allunique(zs)
    @test all(allunique, yss)
    @debug "random_pushpop(xs)" length(zs) length.(yss)
    ys = sort!(foldl(append!, yss; init = copy(zs)))
    @test length(ys) == length(xs)
    @test setdiff(ys, xs) == []
    @test setdiff(xs, ys) == []
    @test ys == xs
    return (; zs, yss)
end

end  # module


text:
# This file was generated, do not modify it. # hide
@load PCA pkg=MultivariateStats

pca_mdl = PCA(pratio=1)
pca = machine(pca_mdl, X)
fit!(pca)
PCA
W = transform(pca, X);


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule LLVM_full_assert_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("LLVM_full_assert")
JLLWrappers.@generate_main_file("LLVM_full_assert", UUID("6ec703ca-3f29-566b-9bb1-b5c9e844abaf"))
end  # module LLVM_full_assert_jll


text:
using Distributions

module GKDistribution

  import Distributions:
         ContinuousUnivariateDistribution,
         location,
         scale,
         skewness,
         kurtosis,
         params,
         minimum,
         maximum,
         insupport,
         quantile,
         logpdf,
         fit_mle,
         Normal,
         cdf,
         pdf

  import Optim:
         optimize
  import Roots:
         fzero

  include("gk.jl")
  export GK,
         logpdf,
         quantile,
         location,
         scale,
         skewness,
         kurtosis,
         asymmetry,
         logpdf,
         quantile,
         params,
         minimum,
         maximum,
         insupport,
         pdf

end # module


text:
struct SRT{Force} <: CollisionModel
    τ::Float64
    force::Force
end
SRT(τ) = SRT(τ, nothing)

function CollisionModel(cm::Type{<:SRT}, q::Quadrature, problem::FluidFlowProblem)
    τ = q.speed_of_sound_squared * lattice_viscosity(problem) + 0.5

    if has_external_force(problem)
        force_field = (x_idx, y_idx, t) -> lattice_force(problem, x_idx, y_idx, t)
        return SRT(τ, force_field)
    else
        return SRT(τ, nothing)
    end
end

function collide!(
    collision_model::SRT{Force},
    q::Quadrature,
    f_in::Populations,
    f_out::Populations;
    time = 0.0,
) where {Force, T <: Real, Populations <: AbstractArray{T, 3}}
    τ = collision_model.τ

    feq = Array{T}(undef, size(f_in, 3))
    f = Array{T}(undef, size(f_in, 3))
    u = zeros(dimension(q))
    if !(Force <: Nothing)
        F = zeros(dimension(q))
    end

    nx, ny, nf = size(f_in)
    @inbounds for x in 1:nx, y in 1:ny
        @inbounds for f_idx in 1:nf
            f[f_idx] = f_in[x, y, f_idx]
        end

        ρ = density(q, f)

        # Momentum
        velocity!(q, f, ρ, u)

        # Temperature
        # temperature = temperature(q, f, ρ, u)
        temperature = 1.0

        if Force <: Nothing
            equilibrium!(q, ρ, u, temperature, feq)
        else
            F .= collision_model.force(x, y, time)

            equilibrium!(q, ρ, u + τ * F, temperature, feq)
        end

        @inbounds for f_idx in 1:nf
            f_out[x, y, f_idx] = (1 - 1 / τ) * f[f_idx] + (1 / τ) * feq[f_idx]
        end
    end
    return
end


text:
function excess_gibbs_free_energy(model::ActivityModel,p,T,z)
    lnγ = log.(activity_coefficient(model,p,T,z))
    return sum(z[i]*R̄*T*lnγ[i] for i ∈ @comps)
end

function eos(model::ActivityModel,V,T,z)
    x = z./sum(z)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    g_pure = VT_gibbs_free_energy.(model.puremodel,V,T)
    p      = pressure.(model.puremodel,V,T)
    g_ideal = sum(z[i]*R̄*T*log(x[i]) for i ∈ @comps)
    return g_E+g_ideal+sum(z[i]*g_pure[i] for i ∈ @comps)-sum(x[i]*p[i] for i ∈ @comps)*V
end

function eos_res(model::ActivityModel,V,T,z)
    x = z./sum(z)
    g_E = excess_gibbs_free_energy(model,V,T,z)
    g_pure_res = VT_chemical_potential_res.(model.puremodel,V,T)
    p_res      = pressure.(model.puremodel,V,T).-sum(z)*R̄*T/V
    return g_E+sum(z[i]*g_pure_res[i][1] for i ∈ @comps)-sum(z[i]*p_res[i] for i ∈ @comps)*V
end

function bubble_pressure(model::ActivityModel,T,x)
    sat = saturation_pressure.(model.puremodel,T)
    p_sat = [tup[1] for tup in sat]
    γ     = activity_coefficient(model,1e-4,T,x)
    p     = sum(x.*γ.*p_sat)
    y     = x.*γ.*p_sat ./ p
    return (p,y)
end


text:
function DebugDynamicMatLawLoss(domain::Domain, E_all::Array{Float64}, w∂E∂u_all::Array{Float64},
    F_tot::Array{Float64}, nn::Function)
   # todo, use fint instead of computed F_tot
   # F_tot =  hcat(domain.history["fint"]...)'
   # define variables
   neles = domain.neles
   nGauss = length(domain.elements[1].weights)
   nstrains = size(E_all,3)

   NT = size(E_all,1)-1
   @assert size(E_all)==(NT+1, neles*nGauss, nstrains)
   @assert size(F_tot)==(NT, domain.neqs)
   # @show E_all[2,:,:]
   E_all = constant(E_all)
   F_tot = constant(F_tot)
   w∂E∂u_all = constant(w∂E∂u_all)

   function cond0(i, ta_loss, ta_σ)
       i<=NT+1
       # i<=2
   end

   function body(i, ta_loss, ta_σ)
       E = E_all[i]
       DE = E_all[i-1]
       w∂E∂u = w∂E∂u_all[i]
       σ0 = read(ta_σ, i-1)

       fint, σ = assembleInternalForce(domain,nn,E,DE,w∂E∂u,σ0)

       # op = tf.print(i, fint, summarize=-1)
       # fint = bind(fint, op)

       # op = tf.print(i, F_tot[i-1], summarize=-1)
       # fint = bind(fint, op)

       # op = tf.print("F_tot",F_tot[i-1], summarize=-1)
       # i = bind(i, op)

       ta_σ = write(ta_σ, i, σ)
       ta_loss = write(ta_loss, i, sum((fint-F_tot[i-1])^2))
       i+1, ta_loss, ta_σ
   end

   σ0 = constant(zeros(neles*nGauss, nstrains))
   ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, σ0)
   ta_loss = TensorArray(NT+1); ta_loss = write(ta_loss, 1, constant(0.0))
   i = constant(2, dtype=Int32)
   _, out, _ = while_loop(cond0, body, [i,ta_loss, ta_σ]; parallel_iterations=20)
   stack(out)[2:NT]
end

@doc """
    for debugging
"""->
function DynamicMatLawLoss(domain::Domain, E_all::Array{Float64}, w∂E∂u_all::Array{Float64},
    F_tot::Array{Float64}, σ_all::Array{Float64}, nn::Function)
   F_tot =  hcat(domain.history["fint"]...)'
   # define variables
   neles = domain.neles
   nGauss = length(domain.elements[1].weights)
   nstrains = size(E_all,3)

   NT = size(E_all,1)-1
   @assert size(E_all)==(NT+1, neles*nGauss, nstrains)
   @assert size(F_tot)==(NT, domain.neqs)
#    println( size(σ_all))
   @assert size(σ_all)==(NT+1, neles*nGauss, nstrains)
   # @show E_all[2,:,:]
   E_all = constant(E_all)
   F_tot = constant(F_tot)
   σ_all = constant(σ_all)
   w∂E∂u_all = constant(w∂E∂u_all)

   function cond0(i, ta_loss, ta_σ)
       i<=NT+1
   end

   function body(i, ta_loss, ta_σ)
       E = E_all[i]
       DE = E_all[i-1]
       w∂E∂u = w∂E∂u_all[i]

    #    σ0 = σ_all[i-1] # ! for debugging
       σ0 = read(ta_σ, i-1)

       fint, σ = assembleInternalForce(domain,nn,E,DE,w∂E∂u,σ0)

       current_loss = sum((fint-F_tot[i-1])^2)
       op = tf.print(i, current_loss, summarize=-1)
       ta_loss = write(ta_loss, i, current_loss)
       op3 = tf.print(i-1, "correct σ",   σ_all[i]', summarize=-1)
       op4 = tf.print(i-1, "σ",   σ', summarize=-1)
    #    op5 = tf.print(i-1, "E",   E, summarize=-1)
    #    op6 = tf.print(i-1, "DE",   DE, summarize=-1)
    #    op7 = tf.print(i-1, "w∂E∂u",   w∂E∂u, summarize=-1)
    #    op8 = tf.print(i-1, "σ_",   σ_', summarize=-1)

    #    i = bind(i, op1)
    #    i = bind(i, op2)
    #    i = bind(i, op3)
    #    i = bind(i, op4)
    #    i = bind(i, op5)
    #    i = bind(i, op6)
    #    i = bind(i, op7)
    #    i = bind(i, op8)

       ta_σ = write(ta_σ, i, σ)
       i+1, ta_loss, ta_σ
   end

   σ0 = constant(zeros(neles*nGauss, nstrains))
   ta_loss = TensorArray(NT+1); ta_loss = write(ta_loss, 1, constant(0.0))
   ta_σ = TensorArray(NT+1); ta_σ = write(ta_σ, 1, σ0)
   i = constant(2, dtype=Int32)
   _, out = while_loop(cond0, body, [i,ta_loss,ta_σ]; parallel_iterations=20)
   total_loss = sum(stack(out)[2:NT])
   return total_loss
end

function DebugDynamicMatLawLoss(domain::Domain, globdat::GlobalData, state_history::Array{T}, fext_history::Array{S}, nn::Function, Δt::Float64) where {T, S}
    # todo convert to E_all, Ftot
    domain.history["state"] = state_history
    F_tot, E_all, w∂E∂u_all = preprocessing(domain, globdat, hcat(fext_history...), Δt)
    DebugDynamicMatLawLoss(domain, E_all, w∂E∂u_all, F_tot, nn)
end


text:
immutable Pareto <: ContinuousUnivariateDistribution
    scale::Float64
    shape::Float64
    function Pareto(sc::Real, sh::Real)
        sc > zero(sc) && sh > zero(sh) || error("shape and scale must be positive")
        new(float64(sc), float64(sh))
    end
end

Pareto() = Pareto(1.0, 1.0)
Pareto(scale::Real) = Pareto(scale, 1.0)

islowerbounded(::Union(Pareto, Type{Pareto})) = true
isupperbounded(::Union(Pareto, Type{Pareto})) = false
isbounded(::Union(Pareto, Type{Pareto})) = false

minimum(d::Pareto) = d.scale
maximum(d::Pareto) = Inf
insupport(d::Pareto, x::Number) = isfinite(x) && x >= d.scale

mean(d::Pareto) = d.shape > 1.0 ? (d.scale * d.shape) / (d.shape - 1.0) : Inf

median(d::Pareto) = d.scale * 2.0^(1.0/d.shape)

mode(d::Pareto) = d.scale

function var(d::Pareto)
    α = d.shape
    α > 2.0 ? (d.scale^2 * α) / ((α - 1.0)^2 * (α - 2.0)) : Inf
end

function skewness(d::Pareto)
    α = d.shape
    α > 3.0 ? ((2.0 * (1.0 + α)) / (α - 3.0)) * sqrt((α - 2.0) / α) : NaN
end

function kurtosis(d::Pareto)
    α = d.shape
    α > 4.0 ? (6.0 * (α^3 + α^2 - 6.0 * α - 2.0)) / (α * (α - 3.0) * (α - 4.0)) : NaN
end

entropy(d::Pareto) = log(d.scale / d.shape) + 1. / d.shape + 1.

function pdf(d::Pareto, q::Real)
    q >= d.scale ? d.shape * (d.scale/q)^d.shape / q : 0.0
end

ccdf(d::Pareto, q::Real) = q >= d.scale ? (d.scale / q)^d.shape : 1.0
logccdf(d::Pareto, q::Real) = q >= d.scale ? d.shape*log(d.scale / q) : 0.0
cdf(d::Pareto, q::Real) = 1.0 - ccdf(d,q)
logcdf(d::Pareto, q::Real) = log1p(-ccdf(d,q))

cquantile(d::Pareto, p::Real) = d.scale / p^(1.0 / d.shape)
quantile(d::Pareto, p::Real) = cquantile(d,1.0-p)

rand(d::Pareto) = d.scale*exp(rand(Exponential(1.0/d.shape)))


text:
using PkgBenchmark

results = benchmarkpkg("EAGO")
show(results)

# specify tag and uncommit to benchmark versus prior tagged version
#tag = "9b5ba7d756e3feac2a9ab5d42af464f1139b8900" # v0.4.1
#results = judge("McCormick", tag)
#show(results)

export_markdown("results.md", results)


text:
include("required.jl")

## Conditional arguments
# save_prepdata = true # should prepared data be overwritten (optional)

## EBD data preparation

# Select subset of columns to load from CSV file (overflows memory otherwise)
cols = [
    "COMMON NAME",
    "SCIENTIFIC NAME",
    "LATITUDE",
    "LONGITUDE",
    "OBSERVATION DATE",
    "GROUP IDENTIFIER",
    "APPROVED",
]

# Load data from CSV files (from file cut with terminal)
@time df = CSV.read(joinpath("data", "raw", "ebd_warblers_cut.csv"), DataFrame; select=cols)

# Prepare data (arrange values & columns)
@time prepare_ebd_data!(df)

# Select subset with specific columns
select!(df, [:species, :commonName, :year, :latitude, :longitude, :groupIdentifier])

# Remove 1 Aleutian Islands observation with positive longitude
filter!(:longitude => <(0.0), df)

# Remove duplicates (BUT not ok for counts, so-so for dates, see group-observation.jl in src/test/)
df_nogroups = filter(:groupIdentifier => ismissing, df)
dropmissing!(df, :groupIdentifier)
df_groups_unique = unique(df, [:species, :groupIdentifier])
df = vcat(df_nogroups, df_groups_unique)

# Remove groupIdentifier column
select!(df, Not(:groupIdentifier))

## Export prepared data
# save_prepdata = true # should prepared data be overwritten (optional)
if (@isdefined save_prepdata) && save_prepdata == true
    # Export file
    @info "Data exported to file (data preparation)"
    CSV.write(joinpath("data", "proc", "ebd_warblers_prep.csv"), df; delim="\t")

    # Update placeholder file (as file is too big for version control)
    placeholder_path = joinpath("data", "proc", "ebd_warblers_prep_placeholder.csv")
    open(placeholder_path, "w") do io
        write(io, string(Dates.now()))
    end
else
    @info "Data not exported (data preparation)"
end

max_stars_repo_name: gabrieldansereau/betadiversity-hotspots

max_issues_repo_name: gabrieldansereau/betadiversity-hotspots

max_forks_repo_name: gabrieldansereau/betadiversity-hotspots


text:
# Named Tuple Macro
macro NT(expr...)
    nmexpr = Expr(:curly, esc(NamedTuple), :($expr))
    tpexpr = Expr(:tuple, (esc(e) for e in expr)...)
    Expr(:call, nmexpr, tpexpr)
end

cpu = identity

exp_mean(xs) = logsumexp(xs) - log(length(xs))
ess_exp(xs) = exp(2 * logsumexp(xs) - logsumexp(2 * xs))

pass_threshold(logweights, threshold) =
    ess_exp(logweights) >= threshold * length(logweights)

function adaptive_resample_particles(
    rng,
    logweights,
    ancestors,
    atable,
    old_particles,
    threshold = 1.0,
)
    if pass_threshold(logweights, threshold)
        return old_particles, logweights, false
    else
        resampled_particles =
            resample_particles(rng, cpu(logweights), ancestors, atable, old_particles)
        return resampled_particles, zero(logweights), true
    end
end

dpre_default(Xs_prev, t) = zeros(eltype(Xs_prev), size(Xs_prev, 2))

Zygote.@nograd pass_threshold

function setindex_ng!(xs, i, x)
    xs[i] = x
    return nothing
end

Zygote.@nograd setindex_ng!

function getindex_ng(xs, is)
    return xs[is]
end

Zygote.@nograd getindex_ng
struct SMCModel{record_history,T<:Any}
    θ::T
end

function SMCModel(; kwargs...)
    if haskey(kwargs, :record_history)
        record_history = kwargs[:record_history]
    else
        record_history = false
    end
    SMCModel{record_history,typeof(kwargs.data)}(kwargs.data)
end

function Base.show(io::IO, x::SMCModel{record_history}) where {record_history}
    print(io, "SMCModel(θ=", x.θ, ")\n")
    print(io, "record_history=", record_history)
end

function smc(
    rng,
    end_t,
    rinit,
    rproposal,
    dinit,
    dprior,
    dproposal,
    dtransition,
    dmeasure,
    dpre = dpre_default;
    # record_history = false,
    threshold = 1.0,
)
    m = SMCModel(
        record_history = false,
        rinit = rinit,
        rproposal = rproposal,
        dinit = dinit,
        dprior = dprior,
        dproposal = dproposal,
        dtransition = dtransition,
        dmeasure = dmeasure,
        dpre = dpre,
        threshold = threshold,
    )

    m(rng, end_t)
end

function record_history(::SMCModel{rh,T}) where {rh,T}
    return rh
end

function (model::SMCModel)(rng, end_t)
    θ = model.θ
    particles = θ.rinit(rng)
    D, n_particles = size(particles)
    logK = convert(typeof(particles[1]), log(n_particles))
    # Calculate log-likelihood of first particles
    logweights = θ.dmeasure(particles, 1) + θ.dprior(particles) - θ.dinit(particles)
    ancestors = zeros(Int, n_particles)

    if (record_history(model))
        particle_history = repeat(particles, outer = (1, 1, end_t))
        logweight_history = repeat(logweights, outer = (1, end_t))
        ancestor_history = zeros(Int, n_particles, end_t)
    else
        particle_history = missing
        logweight_history = missing
        ancestor_history = missing
    end

    atable = AliasTable(logweights[1], n_particles)

    loglik = exp_mean(logweights)
    ess = zeros(typeof(logweights[1]), end_t)
    logliks = zeros(typeof(logweights[1]), end_t)
    fq_ratio = zeros(typeof(logweights[1]), end_t)
    resampled = fill(false, end_t)

    setindex_ng!(logliks, 1, loglik)
    setindex_ng!(fq_ratio, 1, zero(fq_ratio[1]))
    setindex_ng!(ess, 1, ess_exp(logweights))

    particles, logweights = let old_particles = particles, logweights = logweights
        for t = 2:end_t
            ## Resampling
            logweights = logweights .- logsumexp(logweights)
            g_pre = θ.dpre(old_particles, t)
            @assert !any(isnan.(g_pre))
            logweights += g_pre
            ll_pre = logsumexp(logweights)

            resampled_particles, logweights_prev, resampled_t =
                adaptive_resample_particles(
                    rng,
                    logweights,
                    ancestors,
                    atable,
                    old_particles,
                    θ.threshold,
                )
            setindex_ng!(resampled, t, resampled_t)
            if resampled_t
                g_pre = getindex_ng(g_pre, ancestors)
            end
            new_particles = θ.rproposal(rng, resampled_particles, t)
            f_trans = θ.dtransition(resampled_particles, new_particles, t)
            @assert !any(isnan.(f_trans))
            g_measure = θ.dmeasure(new_particles, t)
            @assert !any(isnan.(g_measure))
            q_prop = θ.dproposal(resampled_particles, new_particles, t)
            @assert !any(isnan.(q_prop))
            f_over_q = f_trans - q_prop
            logweights = f_trans + g_measure - q_prop - g_pre + logweights_prev

            logliks_t = logsumexp(logweights) - logsumexp(logweights_prev) + ll_pre
            loglik += logliks_t

            if loglik == -Inf
                break
            end

            setindex_ng!(logliks, t, logliks_t)
            setindex_ng!(fq_ratio, t, exp_mean(f_over_q))
            setindex_ng!(ess, t, ess_exp(logweights))

            if record_history(model)
                particle_history[:, :, t] .= new_particles
                logweight_history[:, t] .= logweights
                ancestor_history[:, t] = ancestors
            end

            old_particles = new_particles
        end
        particles, logweights
    end
    return @NT(
        loglik,
        logliks,
        particles,
        logweights,
        ess,
        fq_ratio,
        resampled,
        particle_history,
        logweight_history,
        ancestor_history
    )
end

function resample_particles(rng, logweights, ancestors, atable, old_particles)
    n_particles = length(logweights)
    lw_max = maximum(logweights)
    W = exp.(logweights .- lw_max)
    residual_resample!(rng, 1:n_particles, W, ancestors, atable)

    resampled_particles = (old_particles[:, ancestors])
    return resampled_particles
end

Zygote.@nograd resample_particles

function calc_filtered_mean(smc_out)
    D, K, T = size(smc_out.particle_history)
    mu = fill(zero(eltype(smc_out.particle_history)), D, T)
    w = Vector{eltype(smc_out.particle_history)}(undef, K)
    xw = Vector{eltype(smc_out.particle_history)}(undef, D)
    for t = 1:T
        w .=
            exp.(
                @view(smc_out.logweight_history[:, t]) .-
                maximum(@view(smc_out.logweight_history[:, t])),
            )
        W = sum(w)
        for i = 1:K
            xw .= @view(smc_out.particle_history[:, i, t]) .* w[i] ./ W
            mu[:, t] .+= xw
        end
    end
    return mu
end

function simulate_backward(rng, particle_history, logweight_history, dtransition, n_draws)
    D, n_particles, T = size(particle_history)
    X = particle_history
    LW = logweight_history

    w = exp.(view(LW, :, T) .- maximum(view(LW, :, T)))
    a = residual_resample!(rng, 1:n_particles, w, Vector{Int}(undef, n_draws))

    X_sm = X[:, a, :]

    F = fill(zero(X[1]), n_draws, n_particles)

    for t = (T-1):-1:1
        F = dtransition_outer(view(X, :, :, t), view(X_sm, :, :, t + 1), dtransition, t + 1)
        LW_sm = calculate_ffbs_weights(view(LW, :, t), F)
        for j = 1:n_draws
            lw = view(LW_sm, :, j)
            w = exp.(lw .- maximum(lw))
            aj = wsample(rng, 1:n_particles, w)
            X_sm[:, j, t] = X[:, aj, t]
        end
    end

    return X_sm
end

function dtransition_outer(Xprev, X, dtransition, t)
    n_particles = size(Xprev)[2]
    n_draws = size(X)[2]
    X1 = view(Xprev, :, repeat(1:n_particles, outer = n_draws))
    X2 = view(X, :, repeat(1:n_draws, inner = n_particles))
    F = dtransition(X1, X2, t)
    return reshape(F, n_particles, n_draws)
end

function calculate_ffbs_weights(ws_x::AbstractVector, F::AbstractMatrix)
    ws_x .+ F
end


text:
# Run this in the REPL, not Juno!!!

using PkgBenchmark
bresult = benchmarkpkg("DynamicalBilliards"; retune = true)


text:
__precompile__()
module NetSim

if !isdefined(Main, :FloatXX)
	const FloatXX = Float64
else
	const FloatXX = Main.FloatXX
end
#print_with_color(:blue, "Simulator.FloatXX = $FloatXX\n")

abstract Network
export Network
"""
	abstract Connection

Connections have fields (see `createconnection` in connections.jl)
	* `w`			parameters
	* `pre`			pointer to neurons of pre layer
	* `post`		pointer to neurons of post layer
	* `prename`		name of pre layer
	* `postname`	name of post layer

PlasticSparseConnection <: Connection
StaticSparseConnection <: Connection
PlasticDenseConnection <: Connection
StaticDenseConnection <: Connection
One2OneConnection <: Connection
TransposeDenseConnection <: Connection
"""
abstract Connection
export Connection
"""
	abstract Neuron

Neurons should at least have field `outp` (used by `weightedprerates`).
"""
abstract Neuron
export Neuron
abstract NeuronParameters
export NeuronParameters
include("layers.jl")
include("networks.jl")
include("connections.jl")
include("neurons.jl")
include("monitor.jl")
include("visualize.jl")

end

max_stars_repo_name: jbrea/NetSim

max_issues_repo_name: jbrea/NetSim

max_forks_repo_name: jbrea/NetSim


text:
using Test,
      Zygote,
      DiffImages,
      ImageCore,
      ImageTransformations,
      Interpolations,
      StaticArrays,
      FiniteDifferences,
      ChainRulesCore,
      CoordinateTransformations,
      Rotations

@testset "DiffImages" begin
    @info "Testing Colorspace modules"
    @testset "Colorspace" begin
        include("colors/conversions.jl")
    end

    @info "Testing Geometry modules"
    @testset "Adjoints" begin
        include("geometry/adjoints.jl")
    end
    @testset "Warps" begin
        include("geometry/warp.jl")
    end
end


text:
using Distributions
using Distributions: MvNormal, logpdf
using StatsBase, StatsFuns, ..Random
using LinearAlgebra, Clustering
using Formatting

eye(d) = Matrix(I, d, d)
eff_ss(W) = 1/sum((W./sum(W)).^2)
weight_perp(W) = let w=W/sum(W); -sum(w.* log.(w))/log(length(W)) end

# softmax which accepts kwargs. This is actually slower than `softmax(x')'`
# in many cases, but :shrugs:, I wrote this a while ago, and is not usually the
# bottleneck of the code.
function softmax2(logp; dims=2)
    p = exp.(logp .- maximum(logp, dims=dims))
    p ./= sum(p, dims=dims)
    return p
end

function logsumexprows(X::AbstractArray{T}) where {T<:Real}
    n = size(X,1)
    out = zeros(T, n)
    for i = 1:n
        out[i] = logsumexp(X[i,:])
    end
    return out
end

# Kitagawa resampling (as used e.g. in SMC). This is sometimes useful if one is
# using resampling within a sequential AMIS scheme.
"""
    multicategorical_kitagawa(p::Vector{T}, m::Int64)

Draw `m` samples from the Categorical distribution defined by the probability
vector `p`. (Note that computationally there is no advantage to enforcing that
`p` sums to one, and hence ``\\sum_i p_i`` may be any positive number; it will
be implicitly renormalized.) This may be considered a stratified multinomial
sample, where the sample space is pre-partitioned into `m` intervals, each of
which obtains one random draw. Nevertheless, the number of repetitions of each
particle does not follow a multinomial distribution hence my clunky name,
"multicategorical". This follows the description given in Hol et al. 2006
("On Resampling Algorithms for Particle Filters"), but of course following
Kitagawa, 1996.
"""
function multicategorical_kitagawa(rng::AbstractRNG, p::Vector{T}, m::Int64) where T <: AbstractFloat
    n = length(p)
    x = zeros(Int, m)     # return
    interval = sum(p)/m   # stratified intervals on [0,1) * sum(p) [sum p = 1 normally]
    cs = p[1]             # cum probability (calc online)
    i = 1                 # category index
#     println(" cs     lb    r     smp")
    for mm in 1:m
        lb = (mm-1)*interval
        r = rand(rng)*interval
        while cs < lb + r
#             println(format("{:.3f} {:.3f} {:.3f} {:.3f}", cs, lb, r, lb+r))
            i += 1
            if i == n
                x[mm:end] .= i
                return x
            end
            @inbounds cs += p[i]
        end
#         println(format("{:.3f} {:.3f} {:.3f} {:.3f}", cs, lb, r, lb+r))
#         println("Choose i = ", i)
        @inbounds x[mm] = i
    end
    return x
end
multicategorical_kitagawa(p::Vector{T}, m::Int64) where T <: AbstractFloat = multicategorical_kitagawa(Random.GLOBAL_RNG, p, m)

################################################################################
##                                                                            ##
##                  Gaussian Mixture Model Utils                              ##
##    --------------------------------------------------------------------    ##
##    Custom GMM fitting utils: this uses both conjugate priors over GMM      ##
##    components, and allows for weighted observations: both are reqd for     ##
##    use in AMIS. I wrote this a year or so ago when Julia utils were        ##
##    less mature; probably some of this can be pointed at standard libs.     ##
##                                                                            ##
################################################################################

function gmm_llh(X, weights, pis, mus, sigmas; disp=false)
    n, p = size(X)
    k = length(pis)
    thrsh_comp = 0.005
    inactive_ixs = pis[:] .< thrsh_comp

    P = zeros(n, k)
    for j = 1:k
        P[:,j] = log_gauss_llh(X, mus[j,:], sigmas[:,:,j],
            bypass=inactive_ixs[j]) .+ log(pis[j])
    end
    P .*= weights
    if disp
        display(P)
        display(logsumexprows(P))
    end
    return logsumexprows(P)
end

function gmm_prior_llh(pis, mus, sigmas, pi_prior, mu_prior, cov_prior)
    # Prior over GMM: Using conjugate distributions (Normal-Inverse-Wishart IIRC)
    d = size(cov_prior, 1)
    ν = pi_prior # alias
    k = length(pis)
    out = zeros(k)
    @views for j = 1:k
        out[j] += logpdf(MvNormal(mu_prior[j,:], sigmas[:,:,j]/ν[j]), mus[j,:])
        out[j] += -(ν[j] + d + 1)*logdet(sigmas[:,:,j])/2
        out[j] += -ν[j]*sum(diag(cov_prior[:,:,j]*inv(sigmas[:,:,j])))/2
        out[j] += (ν[j] - 1)*log(pis[j])
    end
    return sum(out)
end

"""
    log_gauss_llh(X, mu, sigma; bypass=false)

loglikelihood of multivariate Gaussian. Probably no need for this anymore; it
was a bit faster than the built-in, but I think these have largely caught up.
The major point of this function is to have a `bypass`, which ignores the
entire calculation when convenient. Probably could be dealt with elsewhere,
but it's easier at the moment not to switch up the code when it works.
"""
function log_gauss_llh(X, mu, sigma; bypass=false)
    if bypass
        return -ones(size(X, 1))*Inf
    else
        retval = try _log_gauss_llh(X, mu, sigma)
            catch e
                return -ones(size(X, 1))*Inf
            end
        return retval
    end
end

function _log_gauss_llh(X, mu, sigma)
    d = size(X,2)
#     invUT = Matrix(cholesky(inv(sigma)).U)
    invUT = inv(cholesky(sigma).L)
    Z = (X .- mu')*invUT'
    exponent = -0.5*sum(Z.^2, dims=2)
    lognormconst = -d*log(2*pi)/2 -0.5*logdet(sigma) #.+ sum(log.(diag(invUT)))
    return exponent .+ lognormconst
end

"""
    gmm_custom(X, weights, pi_prior, mu_prior, cov_prior; max_iter=100, tol=1e-3, verbose=true)

Fitting a Gaussian Mixture Model (GMM) via MAP EM with *weighted* observations, and
priors over each of the parameters ``\\mu_j, \\Sigma_j, \\pi_j``, ``j \\in 1,\\ldots,k``.
The priors can be useful (i) to avoid exploding likelihood issues, and (ii) to provide
faster convergence given previous computation.
"""
function gmm_custom(X, weights, pi_prior, mu_prior, cov_prior; max_iter=100, tol=1e-3, verbose=true)
    n, p = size(X)
    k = length(pi_prior)
    @assert size(weights) == (n,)
    @assert size(mu_prior) == (k, p)
    @assert size(cov_prior) == (p, p, k)
    pis = pi_prior/sum(pi_prior)
    mus = copy(mu_prior)
    sigmas = copy(cov_prior)

    weights = weights / mean(weights)    # diff. to Cappé et al. due to prior

    thrsh_comp = 0.005
    inactive_ixs = pi_prior[:] .< thrsh_comp
    pi_prior = copy(pi_prior)
    Ns = zeros(6)

    for i in range(1, stop=max_iter)
        # E-step
        rs = reduce(hcat, map(j -> log_gauss_llh(X, mus[j,:], sigmas[:,:,j], bypass=inactive_ixs[j]), 1:k))
        try
            rs .+= log.(pis)[:]'
            catch e
            display(rs)
            display(log.(pis))
            rethrow(e)
        end

        rs = softmax2(rs, dims=2)
        # reweight according to importance weights (see Adaptive IS in General Mix. Cappé et al. 2008)
        rs .*= weights

        # M-step
        Ns = sum(rs, dims=1)
        inactive_ixs = Ns[:] .< 1
        active_ixs = .!inactive_ixs
        if any(inactive_ixs)
            pis[inactive_ixs] .= 0.0
            pi_prior[inactive_ixs] .= 0.0
        end
        pis = Ns[:] + pi_prior[:]

        pis /= sum(pis)

        _mus = reduce(vcat, map(j -> sum(X .* rs[:,j], dims=1) .+ pi_prior[j]*mu_prior[j,:]', findall(active_ixs)))
        _mus ./= vec(Ns[active_ixs] + pi_prior[active_ixs])
        mus[active_ixs,:] = _mus

        @views for j in findall(active_ixs)
            Δx = X .- mus[j, :]'
            Δμ = (mus[j,:] - mu_prior[j,:])'
            sigmas[:,:,j] = (Δx.*rs[:,j])'Δx + pi_prior[j]*(Δμ'Δμ + cov_prior[:,:,j])
            sigmas[:,:,j] ./= (Ns[j] + pi_prior[j] + p + 2)     # normalizing terms from Wishart prior
            sigmas[:,:,j] = (sigmas[:,:,j] + sigmas[:,:,j]')/2 + eye(p)*1e-6   # hack: prevent collapse
        end

    end

    return pis, mus, sigmas
end

function sample_from_gmm(n, pis, mus, covs; shuffle=true)
    k, p = size(mus)
    Ns = rand(Multinomial(n, pis[:]))
    active_ixs = findall(Ns[:] .>= 1)

    ixs = hcat(vcat(1, 1 .+ cumsum(Ns[1:end-1], dims=1)), cumsum(Ns, dims=1))
    out = zeros(n, p)
    for j=active_ixs
        out[ixs[j,1]:ixs[j,2],:] = rand(MvNormal(mus[j,:], covs[:,:,j]), Ns[j])'
    end
    if shuffle
        out = out[randperm(n),:]
    end
    return out
end

"""
    GMM_IS(n, pis, mus, covs, log_f)

GMM importance sampling. Sample `n` values from the GMM defined by the vectors
of parameters, `pis`, `mus`, `covs`, and return log importance weights relative
to the target distribution. The (unnormalized) target distribution is defined
by the log pdf `log_f`.
"""
function GMM_IS(n, pis, mus, covs, log_f)
    S = sample_from_gmm(n, pis, mus, covs, shuffle=false)
    logW = log_f(S) - gmm_llh(S, 1., pis, mus, covs, disp=false);
    return S, logW;
end

################################################################################
##                                                                            ##
##                  Adaptive Mixture Importance Sampling                      ##
##    --------------------------------------------------------------------    ##
##    Essentially following Cappé et al. 2008. This is the main (only)        ##
##    inference procedure in this file, which we use extensively for MTDS.    ##
##                                                                            ##
################################################################################

function amis(log_f, p::MvNormal, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)
    S = rand(p, gmm_smps)'
    logW = zeros(eltype(S), gmm_smps)
    amis(log_f, S, logW, k; nepochs=nepochs, gmm_smps=gmm_smps, IS_tilt=IS_tilt, terminate=terminate, debug=debug)
end

function amis(log_f, S, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)
    logW = zeros(size(S, 1))
    amis(log_f, S, logW, k; nepochs=nepochs, gmm_smps=gmm_smps, IS_tilt=IS_tilt, terminate=terminate, debug=debug)
end

function amis(log_f, S, logW, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)
    @assert !debug "debug facilities removed to reduce loading time (plots etc.)"
    n, p = size(S)

    W = softmax2(logW, dims=1)
    km = kmeans(copy(S'), k, weights=W)

    cmus = zeros(k, p)
    ccovs = zeros(p, p,k)
    for i in range(1, stop=k)
        ixs = findall(x -> isequal(x,i), km.assignments)
        cX = S[ixs, :]; cw = ProbabilityWeights(W[ixs])
        cmus[i,:] = cX' * cw/cw.sum
        ccovs[:,:,i] = StatsBase.cov(cX, cw, corrected=true)
    end

    cpis = zeros(k)
    cnts = countmap(km.assignments)

    for i in 1:k
        try
            cpis[i] = cnts[i]
        catch e
            @warn format("Error: {:s}. Commonly happens when cluster ({:d}) has no assigned points.", e.msg, i)
        end
    end
    cpis /= sum(cpis)
    # if debug
    #     f, axs = PyPlot.subplots(5,3, figsize=(8,12))
    #     plot_is_vs_target(S[:,1:2], W, ax=axs[1,1], c_num=7)
    #     for i = 1:k
    #         axs[1,2].plot(columns(AxPlot.utils.gaussian_2D_level_curve_pts(cmus[i,1:2], ccovs[1:2,1:2,i]))...);
    #     end
    # end

    amis(log_f, cpis, cmus, ccovs; nepochs=nepochs, gmm_smps=gmm_smps, IS_tilt=IS_tilt, terminate=terminate, debug=debug)
end

"""
    amis(log_f, pis, mus, covs::AbstractArray; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)

    amis(log_f, S, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)

    amis(log_f, S, logW, k::Int; nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)

Adaptive mixture importance sampling of Cappé et al. 2008. I wrote this a couple
of years before the docstring, so please check the source vs. the paper to
verify correctness; I may have made a few changes, but I don't recall now.
The arguments are:

(1) `log_f`: the target function (distribution). This should take a single
argument: a set of ``n \\times d`` samples, and return a vector of length ``n`` of
the log of the unnormalized target function (distribution) for each sample.

(2) `pis`: ``k``-length probability vector of probabilities of each cluster.

(3) `mus`: matrix of cluster means: each *row* is a mean vector.

(4) `covs`: batched covariance matrices stacked in Tensor d×d×k for `k`
clusters.

Instead of args (2-4), i.e. the parameters of the GMM, one can instead supply a
matrix `S` of ``n`` initial samples (``n \\times d`` matrix), and optionally a
``n``-length vector `logW` of log weights corresponding to each sample. The final
positional argument is then an integer `k` specifying the number of clusters to
use in the approximating GMM.

Optional kwargs are `nepochs`: the number of iterations of iterated IS;
`gmm_smps`, the number of samples to perform on each iteration, `IS_tilt`, the
exponential tilt applied to the proposal distribution at each epoch, `terminate`
specifies the proportion of `gmm_smps` the effective sample size (Owen, §9) must
reach before early stopping. Finally `debug` provides some (2D) plots at each
epoch which I used when developing.

The AMIS procedure is an iterated importance sampling (IIS) technique. At each
epoch, samples are drawn from a proposal distribution, which here is a Gaussian
Mixture Model. Each sample is assigned an importance weight vs. the target, and
a new proposal distribution is fitted to these weighted samples. Since we use
an iterative EM method for this, the outer IIS iterations are called epochs.
This procedure is iterated `nepochs` times.
"""
function amis(log_f, pis, mus, covs::AbstractArray; nepochs=5, gmm_smps=1000, IS_tilt=1.,
    terminate=0.75, debug=false)
    @assert !debug "debug facilities removed to reduce loading time (plots etc.)"

    S = sample_from_gmm(gmm_smps, pis, mus, covs*IS_tilt, shuffle=false)
    W = ones(eltype(mus), gmm_smps) / gmm_smps
    amis(log_f, pis, mus, covs, S, W; nepochs=nepochs, gmm_smps=gmm_smps,
        IS_tilt=IS_tilt, terminate=terminate, debug=debug)
end

function amis(log_f, pis, mus, covs::AbstractArray, S::AbstractMatrix, W::AbstractVector;
    nepochs=5, gmm_smps=1000, IS_tilt=1., terminate=0.75, debug=false)
    @assert !debug "debug facilities removed to reduce loading time (plots etc.)"

    ν_S, ν_W, log_W = S, W, nothing
    for i = 1:nepochs
        ## #TODO: inefficient, since most 'entry points' have that the weighted
        ## sample (S,W) is taken from the current GMM. Hence the following two
        ## lines should be moved to the *end* of the loop body. I'll need to do
        ## some checks to make sure it doesn't break anything.
        pis, mus, covs = gmm_custom(ν_S, ν_W, pis, mus, covs; max_iter=3, tol=1e-3, verbose=false);
        ν_S = sample_from_gmm(gmm_smps, pis, mus, covs*IS_tilt, shuffle=false)

        log_W = log_f(ν_S) - gmm_llh(ν_S, 1, pis, mus, covs*IS_tilt);
        ν_W = softmax(vec(log_W));

        (i == nepochs) || (eff_ss(ν_W) >= terminate * gmm_smps) && break
    end
    return ν_S, log_W, pis, mus, covs
end

# using PyPlot, AxPlot
# function plot_is_vs_target(S, W; ax=Nothing, c_num=1, cm="tab10", kwargs...)
# 	@assert cm == "tab10" || cm == "Blues"
# 	col_matrix = cm == "tab10" ? cols_tab10 : cols_Blues
# 	rgba_colors = repeat(col_matrix[c_num:c_num,:], size(S, 1), 1)
#     rgba_colors[:, 4] = W/maximum(W)   # alpha
# #     print(rgba_colors)
#     ax = ax==Nothing ? gca() : ax
#     # display(S)
# #     plot_level_curves_all(mus, UTs; ax=ax)
#     ax[:scatter](splat(S)..., c=rgba_colors, kwargs...)
# end

max_stars_repo_name: ornithos/SeqAdaptiveIS

max_issues_repo_name: ornithos/SeqAdaptiveIS

max_forks_repo_name: ornithos/SeqAdaptiveIS


text:
using Test
using WavePropBase
using WavePropBase.Geometry
using WavePropBase.Interpolation
using StaticArrays

@testset "Cheb tensor product" begin
    # 1d
    n = 10
    xcheb = [Interpolation.chebnode(i,n) for i in 1:n]
    f = x->cos(x)*exp(x)
    vals = map(f,xcheb)
    domain = HyperRectangle(-1,1)
    p = ChebInterp(vals,domain)
    @test interpolation_nodes(p,1) ≈ xcheb
    xtest = 0.1
    @test f(xtest) ≈ p(xtest)
    # shift
    xcheb =  xcheb .+ (1 + π)
    f = x->cos(x)*exp(x)
    vals = map(f,xcheb)
    domain = HyperRectangle(π,π+2)
    p = ChebInterp(vals,domain)
    @test interpolation_nodes(p,1) ≈ xcheb
    xtest = π+0.1
    @test f(xtest) ≈ p(xtest)
    # 2d
    # domain = HyperRectangle((-1,0),(2,1))
    # xcheb = cheb_nodes(10) .+ 1 .+ π
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
# gemm on various size matrices

function matmultest(n, iter)
    a = rand(n,n)
    b = similar(a)
    for i=1:iter
        A_mul_B!(b, a, a)
    end
    b
end

problemsizes = [(2, 10^6, "tiny"), (2^4, 10^5, "small"), (2^6, 10^4, "medium"), (2^8, 10^2, "large")]
testdata = [(gemvtest, "matmul", "matrix-matrix multiplication", problemsizes)]
include("../perfgeneric.jl")

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
"""
    AbstractNode

At the top of the type hierarchy of `StenoGraphs`.
Anything that might resemble a node.
"""
abstract type AbstractNode end

"""
    Node

Subtype of [`AbstractNode`](@ref).
Any subtype of `Node` must have the field `node` (and no other) that uniquily identifies the node.
"""
abstract type Node <: AbstractNode end

"""
    MetaNode

Subtype of [`AbstractNode`](@ref).
Any subtype of `MetaNode` must have a field `node` of subtype [`Node`](@ref), but may have any number of other fields containing metadata.
"""
abstract type MetaNode <: AbstractNode end

"""
    SimpleNode

Constructs a subtype of [`Node`](@ref) with a symbol as identifier.

# Example

```jldoctest
julia> SimpleNode(:a)
a
```
"""
struct SimpleNode{T <: Symbol} <: Node
    node::T
end

"""
    Node(node)

Alias for [`SimpleNode`](@ref).

# Example

```jldoctest
julia> Node(:a)
a
```
"""
Node(node) = SimpleNode(node)

"""
    AbstractEdge

At the top of the type hierarchy of `StenoGraphs`.
Anything that might resemble an edge.
"""
abstract type AbstractEdge end

"""
    Edge

Subtype of [`AbstractEdge`](@ref).
Any subtype of `Edge` must have the fields `src` and `dst` (and no other), which must be a subtype of [`AbstractNode`](@ref).
Any implementation might be stricter about typing.
"""
abstract type Edge <: AbstractEdge end

"""
    MetaEdge

Subtype of [`AbstractEdge`](@ref).
Any subtype of `MetaEdge` must have a field `edge` of subtype [`Edge`](@ref), but may have any number of other fields containing metadata.
"""
abstract type MetaEdge <: AbstractEdge end

"""
    Edge(src, dst)

Alias for [`DirectedEdge`](@ref).

# Example

```jldoctest
julia> Edge(Node(:a), Node(:b))
a → b
```

"""
Edge(src, dst) = DirectedEdge(src, dst)

"""
    DirectedEdge(src, dst)

Subtype of [`Edge`](@ref).
Directed edge from `src` to `dst`.

# Example

```jldoctest
julia> DirectedEdge(Node(:a), Node(:b))
a → b
```
"""
struct DirectedEdge{T1 <: AbstractNode, T2 <: AbstractNode} <: Edge
    src::T1
    dst::T2
end

"""
    UndirectedEdge(src, dst)

Subtype of [`Edge`](@ref).
Undirected edge from `src` to `dst`. What is what does not matter.

# Example
```jldoctest
julia> e1 = UndirectedEdge(Node(:a), Node(:b))
a ↔ b

julia> e2 = UndirectedEdge(Node(:b), Node(:a))
b ↔ a

julia> isequal(e1, e2)
true

julia> unique([e1, e2])
a ↔ b
```
"""
struct UndirectedEdge{T1 <: AbstractNode, T2 <: AbstractNode} <: Edge
    src::T1
    dst::T2
end

"""
    Modifier

The abstract type that powers EdgeModifier and NodeModifier.
`StenoGraphs` does not implement any concrete modifiers.

"""
abstract type Modifier end

"""
    EdgeModifier

Subtype of [`Modifier`](@ref).
[`ModifiedEdge`](@ref)s require `EdgeModifier`s.
`EdgeModifier`s usually make use of `*` for creating [`ModifiedEdge`](@ref)s/[`ModifyingNode`](@ref)s.
One special application for `EdgeModifier`s is the creation of [`ModifyingNode`](@ref)s.
Since `StenoGraphs` does not implement any `EdgeModifier` users must implement them.
If these are not atomic they must take care to implement comparison methods (see examples in [`NodeModifier`](@ref))

# Example

```jldoctest
# `StenoGraphs` does not implement any `EdgeModifier`s
julia> struct Weight{N <: Number} <: EdgeModifier w::N end

julia> ModifiedEdge(Edge(Node(:a), Node(:b)), Weight(.5)) == # directly create ModifiedEdge
        Edge(Node(:a), Node(:b)) * Weight(.5) == # modify an edge
        Edge(Node(:a), Node(:b) * Weight(.5)) # modify Edge through a ModifyingNode
true

```
"""
abstract type EdgeModifier <: Modifier end

"""
    NodeModifier

Subtype of [`Modifier`](@ref).
[`ModifiedNode`](@ref)s require `NodeModifier`s.
`NodeModifier`s usually make use of `^` for creating [`ModifiedNode`](@ref)s.
Since `StenoGraphs` does not implement any `NodeModifier` users must implement them.
If these may contain any mutable fields (i.e. strings, vectors, arrays, etc.) users must take care to implement comparison methods.

# Example

```jldoctest
# `StenoGraphs` does not implement any `NodeModifier`s
julia> struct Label <: NodeModifier l end

julia> import Base.==

julia> ==(x::Label, y::Label) = x.l == y.l;

julia> ModifiedNode(Node(:a), Label("hi")) == Node(:a)^Label("hi")
true

```
"""
abstract type NodeModifier <: Modifier end
const NodeOrEdgeModifier = Union{EdgeModifier,NodeModifier}

"""
    ModifiedEdge

Subtype of [`MetaEdge`](@ref) that contains two fields (`edge` and `modifiers`).
`modifiers` is a `Dict{Symbol, EdgeModifier}` where the keys are `nameof(typeof(EdgeModifier))`.
Modifiying an edge with several modifiers of the same type will therefore overwrite old modifiers.
A `ModifiedEdge` is created by modifying an edge directly (with `*`) or via [`ModifyingNode`](@ref)s where a node is modified that than modifies the edge.

# Example
```jldoctest
julia> struct Weight{N <: Number} <: EdgeModifier w::N end

julia> struct Start{N <: Number} <: EdgeModifier s::N end

julia> ModifiedEdge(Edge(Node(:a), Node(:b)), Weight(3))
a → b * Weight{Int64}(3)

julia> ModifiedEdge(Edge(Node(:a), Node(:b)), Weight(3)) == Edge(Node(:a), Node(:b)) * Weight(3)
true

julia> ModifiedEdge(Edge(Node(:a), Node(:b)), [Weight(3), Start(2)])
a → b * [Start{Int64}(2), Weight{Int64}(3)]

julia> @StenoGraph a → b * Weight(3) * Start(2)
a → b * [Start{Int64}(2), Weight{Int64}(3)]

julia> @StenoGraph a → b * Weight(3) * Start(2) * Weight(2)
a → b * [Start{Int64}(2), Weight{Int64}(2)]
```
"""
struct ModifiedEdge{E <: AbstractEdge, DM <: AbstractDict{S, M} where {S <: Symbol, M <: EdgeModifier}} <: MetaEdge
    edge::E
    modifiers::DM
end

"""
    ModifyingNode

Subtype of [`MetaNode`](@ref) that contains two fields (`node` and `modifiers`).
`modifiers` is a `Dict{Symbol, EdgeModifier}` where the keys are `nameof(typeof(EdgeModifier))`.
Modifiying an node with several modifiers of the same type will therefore overwrite old modifiers.
A `ModifyingNode` is created by multiplying it (`*`) with an *`EdgeModifier`* since it will modify edges build upon it.
If you want to modify the node use `^` and see [`ModifiedNode`](@ref).

# Example
```jldoctest
julia> struct Weight{N <: Number} <: EdgeModifier w::N end

julia> struct Start{N <: Number} <: EdgeModifier s::N end

julia> @StenoGraph a → b * Weight(3) * Start(2)
a → b * [Start{Int64}(2), Weight{Int64}(3)]

julia> @StenoGraph a → b * Weight(3) * Start(2) * Weight(2)
a → b * [Start{Int64}(2), Weight{Int64}(2)]
```
"""
struct ModifyingNode{N <: AbstractNode, DM <: AbstractDict{S, M} where {S <: Symbol, M <: EdgeModifier}} <: MetaNode
    node::N
    modifiers::DM
end

"""
    ModifiedNode

Subtype of [`MetaNode`](@ref) that contains two fields (`node` and `modifiers`).
`modifiers` is a `Dict{Symbol, NodeModifier}` where the keys are `nameof(typeof(NodeModifier))`.
Modifiying a node with several modifiers of the same type will therefore overwrite old modifiers.

# Example

```jldoctest
julia> struct Observed <: NodeModifier end

julia> struct Label{N <: String} <: NodeModifier s::N end

julia> Node(:b) ^ Label("some label") ^ Observed()
b^[Observed(), Label{String}("some label")]

julia> Node(:b) ^ Label("some label") ^ Observed() ^ Label("some other label")
b^[Observed(), Label{String}("some other label")]
```
"""
struct ModifiedNode{N <: AbstractNode, DM <: AbstractDict{S, M} where {S <: Symbol, M <: NodeModifier}} <: MetaNode
    node::N
    modifiers::DM
end

struct Arrow{VE <: Vector{E} where {E <: AbstractEdge}} <: MetaEdge
    edges::VE
    lhs
    rhs
end


text:
using Pkg.Artifacts
using Pkg.TOML

WEIGHT_DIR = joinpath(@__DIR__, "weights")
BUILD_DIR = joinpath(@__DIR__, "build")
URL = "https://github.com/darsnack/MetalheadWeights/releases/download"
ARTIFACT_TOML = joinpath(@__DIR__, "Artifacts.toml")

function build(version; weights = WEIGHT_DIR,
                        build = BUILD_DIR,
                        artifact_toml = ARTIFACT_TOML,
                        url = URL)
    # Clean up prior builds.
    ispath(build) && rm(build; recursive=true, force=true)
    mkdir(build)

    # Package up weights
    for weight in readdir(weights)
        model = split(basename(weight), ".")[1]
        artifact_filename = "$model-$version.tar.gz"
        product_hash = create_artifact() do artifact_dir
            mkdir(joinpath(artifact_dir, model))
            cp(joinpath(weights, weight), joinpath(artifact_dir, weight); force=true)
        end
        download_hash = archive_artifact(product_hash, joinpath(build, artifact_filename))
        remote_url = "$url/v$version/$artifact_filename"

        @info "Creating:" model version product_hash artifact_filename download_hash remote_url
        bind_artifact!(
            artifact_toml,
            "$model",
            product_hash,
            force=true,
            lazy=true,
            download_info=Tuple[(remote_url, download_hash)])
    end
end

max_stars_repo_name: darsnack/MetalheadWeights

max_issues_repo_name: darsnack/MetalheadWeights

max_forks_repo_name: FluxML/MetalheadWeights


text:
# Make `map!` fast if the inputs and output share tokens

# TODO consider if `map` should respect iteration order or indices??

function Base.map!(f, out::AbstractDictionary, d::AbstractDictionary, d2::AbstractDictionary, ds::AbstractDictionary...)
    if sharetokens(out, d, d2, ds...)
        @inbounds for t in tokens(out)
            settokenvalue!(out, t, f(gettokenvalue(d, t), gettokenvalue(d2, t), map(x -> @inbounds(gettokenvalue(x, t)), ds)...))
        end
    else
        @boundscheck nothing # TODO check that indices match
        @inbounds for i in keys(out)
            out[i] = f(d[i], d2[i], map(x -> @inbounds(x[i]), ds)...)
        end
    end
    return out
end

# (avoid an _apply for the two-input case)
function Base.map!(f, out::AbstractDictionary, d::AbstractDictionary, d2::AbstractDictionary)
    if sharetokens(out, d, d2)
        @inbounds for t in tokens(out)
            settokenvalue!(out, t, f(gettokenvalue(d, t), gettokenvalue(d2, t)))
        end
    else
        @boundscheck nothing # TODO check that indices match
        @inbounds for i in keys(out)
            out[i] = f(d[i], d2[i])
        end
    end
    return out
end

function Base.map!(f, out::AbstractDictionary, d::AbstractDictionary)
    if sharetokens(out, d)
        @inbounds for t in tokens(out)
            settokenvalue!(out, t, f(gettokenvalue(d, t)))
        end
    else
        @boundscheck nothing # TODO check that indices match
        @inbounds for i in keys(out)
            out[i] = f(d[i])
        end
    end
    return out
end

function Base.map!(f, out::AbstractDictionary)
    if istokenizable(out)
        @inbounds for t in tokens(out)
            settokenvalue!(out, t, f())
        end
    else
        @inbounds for i in keys(out)
            out[i] = f()
        end
    end
    return out
end

function Base.map(f, d::AbstractDictionary)
    out = similar(d, Core.Compiler.return_type(f, Tuple{eltype(d)}))
    map!(f, out, d)
    return out
end

function Base.map(f, d::AbstractDictionary, ds::AbstractDictionary...)
    out = similar(d, Core.Compiler.return_type(f, Tuple{eltype(d), map(eltype, ds)...}))
    map!(f, out, d, ds...)
    return out
end