text:
using StatsBase, Printf, Distances

export EditDistance, EditDist
export FastEditDistance, FastEditDist
export FixedPenaltyEditDistance, FixPenEditDist
export AvgSizeFpEditDistance
export print_info, get_info, get_info_deep

# EditDistance
# ---
struct EditDistance{T<:Metric} <: Metric
    ground_dist::T
end

const EditDist = EditDistance

function (d::EditDistance)(S1::Vector{T}, S2::Vector{T}) where {T}
    if length(S1) < length(S2)  # This ensures first seq is longest
        d(S2, S1)
    else
        d₀ = d.ground_dist
        prev_row = pushfirst!(cumsum([d₀(nothing, p) for p in S2]), 0.0);
        curr_row = zeros(Float64, length(S2) + 1);

        for i = 1:length(S1)
            curr_row[1] = prev_row[1] + length(S1[i])
            for j = 1:(length(S2))
                # @show i, j, prev_row[j], d.ground_dist(S1[i], S2[j])
                curr_row[j+1] = min(prev_row[j] + d₀(S1[i], S2[j]),
                                        prev_row[j+1] + d₀(nothing, S1[i]),
                                        curr_row[j] + d₀(nothing, S2[j]))
            end
            # @show curr_row
            copy!(prev_row, curr_row)
        end
        return curr_row[end]
    end
end

# EditDistance with Memory
# ---------------

struct FastEditDistance{T<:Metric} <: Metric
    ground_dist::T
    curr_row::Vector{Float64}
    prev_row::Vector{Float64}
    function FastEditDistance(ground_dist::S, K::Int) where {S<:Metric}
        new{S}(ground_dist, zeros(Float64, K), zeros(Float64, K))
    end
end

const FastEditDist = FastEditDistance

function Base.show(io::IO, d::FastEditDistance)
    print(io, "EditDistance (max num. interactions $(length(d.curr_row))) with $(d.ground_dist) ground distance.")
end

function (d::FastEditDistance)(S1::Vector{T}, S2::Vector{T}) where {T}
    if length(S1) < length(S2)  # This ensures first seq is longest
        d(S2, S1)
    else
        d₀ = d.ground_dist
        prev_row = view(d.prev_row, 1:(length(S2)+1))
        curr_row = view(d.curr_row, 1:(length(S2)+1))
        # prev_row = d.prev_row
        # curr_row = d.curr_row

        prev_row[1] = 0.0
        for i in 1:length(S2)
            prev_row[i+1] = prev_row[i] + d₀(nothing, S2[i])
        end
        curr_row .= 0.0

        @views for i = 1:length(S1)
            curr_row[1] = prev_row[1] + length(S1[i])
            for j = 1:(length(S2))
                # @show i, j, prev_row[j], d.ground_dist(S1[i], S2[j])
                curr_row[j+1] = min(
                    prev_row[j] + d₀(S1[i], S2[j]),
                    prev_row[j+1] + d₀(nothing, S1[i]),
                    curr_row[j] + d₀(nothing, S2[j])
                    )
            end
            # @show curr_row
            copy!(prev_row, curr_row)
        end
        return curr_row[length(S2) + 1]
    end
end

function (d::Union{EditDist,FastEditDist})(X::Nothing, Y::Vector{T})::Float64 where {T}
    return sum(p->d.ground_dist(nothing,p), Y)
end
function (d::Union{EditDist,FastEditDist})(X::Vector{T}, Y::Nothing)::Float64 where {T}
    return d(Y, X)
end
function (d::Union{EditDist,FastEditDist})(X::Nothing, Y::Nothing)::Float64 where {T}
    return 0.0
end

function print_info(
    d::Union{EditDistance,FastEditDistance},
    S1::Vector{T}, S2::Vector{T}
    ) where {T}

    d₀ = d.ground_dist
    # First find the substitution matrix
    C = zeros(Float64, length(S1)+1, length(S2)+1)
    C[:,1] = pushfirst!(cumsum([d₀(Λ, p) for p in S1]), 0.0);
    C[1,:] = pushfirst!(cumsum([d₀(Λ, p) for p in S2]), 0.0);

    for j in 1:length(S2)
        for i in 1:length(S1)
            C[i+1,j+1] = minimum([
                C[i,j] + d₀(S1[i], S2[j]),
                C[i,j+1] + d₀(nothing, S2[j]),
                C[i+1,j] + d₀(nothing, S1[i])
            ])
        end
    end
    # Now retrace steps to determine an optimal matching
    i, j = size(C)
    pairs = Tuple{Int,Int}[]
    while (i ≠ 1) & (j ≠ 1)
        if C[i,j] == (C[i,j-1] + d₀(nothing, S2[j-1]))
            pushfirst!(pairs, (0,j-1))
            j = j-1
        elseif C[i,j] == (C[i-1,j] + d₀(nothing, S1[i-1]))
            pushfirst!(pairs, (i-1,0))
            i = i-1
        else
            pushfirst!(pairs,(i-1,j-1))
            i = i-1; j = j-1
        end
    end
    for k in Iterators.reverse(1:(i-1))
        pushfirst!(pairs, (k,0))
    end
    for k in Iterators.reverse(1:(j-1))
        pushfirst!(pairs, (0,k))
    end
    max_len = maximum(map(x->length(@sprintf("%s",x)), S1))
    # @show outputs
    title = "\nOptimal Matching"
    println(title)
    println("-"^length(title), "\n")
    for (k,l) in pairs
        if k == 0
            tmp_S1 = "Λ"
            tmp_S2 = S2[l]
        elseif l == 0
            tmp_S2 = "Λ"
            tmp_S1 = S1[k]
        else
            tmp_S1, tmp_S2 = (S1[k], S2[l])
        end
        pad = max_len - length(@sprintf("%s", tmp_S1))
        println("$tmp_S1" * " "^pad, " → $tmp_S2")
    end

end

function get_info(
    d::Union{EditDistance,FastEditDistance},
    x::Vector{T}, y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    # First find the substitution matrix
    C = zeros(Float64, length(x)+1, length(y)+1)
    C[:,1] = pushfirst!(cumsum([d_g(nothing, p) for p in x]), 0.0);
    C[1,:] = pushfirst!(cumsum([d_g(nothing, p) for p in y]), 0.0);

    for j in 1:length(y)
        for i in 1:length(x)
            C[i+1,j+1] = minimum([
                C[i,j] + d_g(x[i], y[j]),
                C[i,j+1] + d_g(nothing, x[i]),
                C[i+1,j] + d_g(nothing, y[j])
            ])
        end
    end
    # Now retrace steps to find indices of matched paths
    i, j = size(C)
    indx, indy = (zeros(Bool, length(x)), zeros(Bool, length(y)))
    while (i ≠ 1) & (j ≠ 1)
        if C[i,j] == (C[i,j-1] + d_g(nothing, y[j-1]))
            j = j-1
        elseif C[i,j] == (C[i-1,j] + d_g(nothing, x[i-1]))
            i = i-1
        else
            i = i-1; j = j-1
            indx[i], indy[j] = (true,true)
        end
    end
    return indx, indy
end

struct FixedPenaltyEditDistance{T<:Metric} <: Metric
    ground_dist::T
    ρ::Real
end

const FixPenEditDist = FixedPenaltyEditDistance

function (d::FixPenEditDist)(S1::Vector{T}, S2::Vector{T}) where {T}
    if length(S1) < length(S2)
        d(S2, S1)
    else
        prev_row = d.ρ * collect(0:length(S2));
        curr_row = zeros(Float64, length(S2) + 1);

        for i = 1:length(S1)
            curr_row[1] = i*d.ρ
            for j = 1:(length(S2))
                curr_row[j+1] = minimum([prev_row[j] + d.ground_dist(S1[i], S2[j]),
                                        prev_row[j+1] + d.ρ,
                                        curr_row[j] + d.ρ])
            end
            prev_row = copy(curr_row)
        end
        return curr_row[end]
    end
end

function (d::FixPenEditDist)(X::Nothing, Y::Vector{T})::Float64 where {T}
    return d.ρ * length(Y)
end
function (d::FixPenEditDist)(X::Vector{T}, Y::Nothing)::Float64 where {T}
    d(Y, X)
end
function (d::FixPenEditDist)(X::Nothing, Y::Nothing)::Float64 where {T}
    return 0.0
end

function get_info(
    d::FixPenEditDist,
    x::Vector{T}, y::Vector{T}
    ) where {T}
    C = zeros(Float64, length(x)+1, length(y)+1)

    C[:,1] = [d.ρ * i for i = 0:length(x)]
    C[1,:] = [d.ρ * i for i = 0:length(y)]

    for j = 1:length(y)
        for i = 1:length(x)
            C[i+1,j+1] = minimum([
            C[i,j] + d.ground_dist(x[i], y[j]),
            C[i,j+1] + d.ρ,
            C[i+1,j] + d.ρ
            ])
        end
    end
    # Now retrace steps to determine an optimal matching
    i, j = size(C)
    indx, indy = (zeros(Bool, length(x)), zeros(Bool, length(y)))
    while (i ≠ 1) | (j ≠ 1)
        if C[i,j] == (C[i-1,j] + d.ρ)
            i = i-1
        elseif C[i,j] == (C[i,j-1] + d.ρ)
            j = j-1
        else
            i = i-1; j = j-1
            indx[i], indy[j] = (true,true)
        end
    end
    return indx, indy
end

function print_info(
    d::Union{EditDistance,FastEditDistance,FixPenEditDist},
    x::Vector{T}, y::Vector{T}
    ) where {T}

    indx, indy = get_info(d, x, y)
    max_len = maximum(map(z->length(@sprintf("%s",z)), x))
    ix , iy = (0,0)
    title = "\nOptimal Matching"
    println(title)
    println("-"^length(title), "\n")
    while true
        ixerr = isnothing(findnext(indx, ix+1)) ? length(x) : findnext(indx, ix+1)-1
        iyerr = isnothing(findnext(indy, iy+1)) ? length(y) : findnext(indy, iy+1)-1
        for j in (ix+1):ixerr
            tmp1 = x[j]
            tmp2 = "Null"
            pad = max_len - length(@sprintf("%s", tmp1))
            println(" "^pad * "$tmp1" , " → $tmp2")
        end
        for j in (iy+1):iyerr
            tmp1 = "Null"
            tmp2 = y[j]
            pad = max_len - length(@sprintf("%s", tmp1))
            println(" "^pad * "$tmp1", " → $tmp2")
        end
        if isnothing(findnext(indx, ix+1))
            break
        end
        ix = findnext(indx, ix+1)
        iy = findnext(indy, iy+1)
        tmp1,tmp2 = (x[ix],y[iy])
        pad = max_len - length(@sprintf("%s", tmp1))
        println(" "^pad * "$tmp1", " → $tmp2")
    end
end

function get_info_deep(
    d::Union{EditDistance,FastEditDistance,FixPenEditDist},
    x::Vector{T}, y::Vector{T}
    ) where {T}

    d_g = d.ground_dist
    indx, indy = get_info(d, x, y)
    outx, outy = (
        [zeros(Bool,i) for i in length.(x)],
        [zeros(Bool,i) for i in length.(y)]
    )
    ix , iy = (0,0)
    while true
        if isnothing(findnext(indx, ix+1))
            break
        end
        ix = findnext(indx, ix+1)
        iy = findnext(indy, iy+1)
        outx[ix], outy[iy] = get_info(d_g, x[ix], y[iy])
    end
    return outx, outy
end

struct AvgSizeFpEditDistance{T<:Metric} <: Metric
    ground_dist::T
    ρ::Real
end

function (d::AvgSizeFpEditDistance)(S1::Vector{T}, S2::Vector{T}) where {T}

    d_ed = FixPenEditDist(d.ground_dist, d.ρ)(S1, S2)

    return d_ed + (mean(length.(S1)) - mean(length.(S2)))^2

end

ext: jl


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
VERSION

# %%
using Plots
default(size=(400, 250))

using Random
Random.seed!(4649373)

# %%
x = range(0, 1; length=10)

f(x) = sinpi(2x)
noise = 0.2randn(length(x))
t = f.(x) + noise

xs = range(extrema(x)...; length=1000)
plot(xs, f.(xs); label="", xlabel="x", ylabel="t")
scatter!(x, t; label="", color=3)

# %%
PP = []
for d in (0, 1, 3, 9)
    X = x.^(0:d)'
    w = X\t
    g(x) = evalpoly(x, w)
    P = plot(xs, f.(xs); label="")
    plot!(xs, g.(xs); label="degree $d")
    scatter!(x, t; label="")
    push!(PP, P)
end
plot(PP...; size=(800, 500))

# %%
PP = []
for d in (0, 1, 3, 9)
    X = x.^(0:d)'
    w = X'X\X't
    g(x) = evalpoly(x, w)
    P = plot(xs, f.(xs); label="")
    plot!(xs, g.(xs); label="degree $d")
    scatter!(x, t; label="")
    push!(PP, P)
end
plot(PP...; size=(800, 500))

# %%
using SymPy
n = 5
x = collect(symbols("x1:$(n+1)", real=true))

# %%
t = collect(symbols("t1:$(n+1)", real=true))

# %%
d = 3
X = x.^(0:d)'

# %%
A = X'X

# %%
T = X't

# %%

max_stars_repo_path: 0005/PRML Chap.1 Ex.1.1.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0005/PRML Chap.1 Ex.1.1.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0005/PRML Chap.1 Ex.1.1.jl

max_forks_repo_name: genkuroki/public


text: mapvals(f, dictlike) = [k => f(v) for (k, v) in dictlike]


text:
# Lik.jl

export llik, llik_oos!, Qresid!, Qresid_cond!;

function ldensy_marg(y::Vector{T}, μ::T, σ2::T, Cov::PDMat) where T <: Real
    ## assumes mean zero
    n = length(y)
    Σ = Cov + σ2*I
    ldetΣ = logdet(Σ)
    d = (y - μ)
    Σinvd = Σ \ d

    -0.5 * ( n*log(2π) + ldetΣ + d'Σinvd )
end

## naive is faster than:
# logpdf( MvNormal( Cov + σ2*I ), y)

function ldensy_marg(y::Vector{T}, μ::T, σ2::T, κ::T, Cor::PDMat) where T <: Real
    ## assumes mean zero
    ## κ is the signal-to-noise ratio so that Σ = κ*σ2*Cor
    n = length(y)
    Σ0 = (κ .* Cor + I)
    ldetΣ = n*log(σ2) + logdet(Σ)
    d = (y - μ)
    Σinvd = Σ \ d

    -0.5 * ( n*log(2π) + ldetΣ + d'Σinvd )
end

function llik(model::Model_GPMTD)

    llik = 0.0
    lp = Vector{Float64}(undef, model.L+1)

    for i = 1:model.n

        lp[1] = model.state.lλ[1] - 0.5*log(2π*model.state.intercept.σ2) -
                0.5 * (model.y[i] - model.state.intercept.μ)^2 / model.state.intercept.σ2

        for j in 1:model.L
            lp[j+1] = model.state.lλ[j+1] -
                      0.5*log(2π*model.state.mixcomps[j].σ2) -
                      0.5 * (model.y[i] - model.state.mixcomps[j].μ - model.state.mixcomps[j].fx[i])^2 / model.state.mixcomps[j].σ2
        end

        # m = maximum(lp)
        # llik += m + log(sum(exp.(lp .- m))) # logsumexp
        llik += SparseProbVec.logsumexp(lp)
    end

    return llik
end

function llik_oos!(model::Model_GPMTD, postsamp::Dict{Symbol, Any},
                  y_oos::Vector{T}, D_oos::Vector{Array{T, 2}},
                  D_oos_is::Vector{Array{T, 2}}; tol=1.0e-9) where T <: Real

    ## Modifes the state of the rng in the model object

    # D_oos[j] = pairDistMat( X_oos[:,[j]] )
    # D_oos_is[j] = pairDistMat( X_oos[:,[j]],  model.X[:,[j]] )

    llik = 0.0
    lp = Vector{Float64}(undef, model.L+1)

    Fstar = [ postsamp[:mixcomps][j][:μ] .+ rfullcond_fstar!(postsamp[:mixcomps][j][:fx],
        covMat(model.D[j], postsamp[:mixcomps][j][:κ]*postsamp[:mixcomps][j][:σ2], postsamp[:mixcomps][j][:corParams]),
        covMat(D_oos[j], postsamp[:mixcomps][j][:κ]*postsamp[:mixcomps][j][:σ2], postsamp[:mixcomps][j][:corParams]),
        covMat(D_oos_is[j], postsamp[:mixcomps][j][:κ]*postsamp[:mixcomps][j][:σ2], postsamp[:mixcomps][j][:corParams]),
        model.state.mixcomps[j].rng, tol) for j = 1:model.L ]

    for i = 1:length(y_oos)
        lp[1] = postsamp[:lλ][1] - 0.5*log(2π*postsamp[:intercept].σ2) -
                0.5 * (y_oos[i] - postsamp[:intercept].μ)^2 / postsamp[:intercept].σ2

        for j in 1:model.L
            lp[j+1] = postsamp[:lλ][j+1] -
                    0.5*log(2π*postsamp[:mixcomps][j][:σ2]) -
                    0.5 * (y_oos[i] - Fstar[j][i])^2 / postsamp[:mixcomps][j][:σ2]
        end

        llik += SparseProbVec.logsumexp(lp)
    end

    return llik
end

function Qresid!(model::Model_GPMTD, sims::Vector{Any})
    ### Add quantile residuals (conditional model CDF evaluated at each observation) to posterior simulations

    nsim = length(sims)
    lCDF = Vector{Float64}(undef, model.n)
    lP = Vector{Float64}(undef, model.L+1)

    for ii = 1:nsim
        for i = 1:model.n

            lP[1] = sims[ii][:lλ][1] + logcdf( Normal(sims[ii][:intercept].μ, sqrt(sims[ii][:intercept].σ2)), model.y[i] )

            for j in 1:model.L
                lP[j+1] = sims[ii][:lλ][j+1] + logcdf( Normal(sims[ii][:mixcomps][j][:μ] + sims[ii][:mixcomps][j][:fx][i],
                                                               sqrt(sims[ii][:mixcomps][j][:σ2])),
                                                        model.y[i] )
            end

            lCDF[i] = SparseProbVec.logsumexp(lP)
        end

        sims[ii][:CDF] = exp.(lCDF)
    end

    return nothing
end

function Qresid_cond!(model::Model_GPMTD, sims::Vector{Any})
    ### Add quantile residuals (conditional component CDF evaluated at each observation) to posterior simulations

    nsim = length(sims)
    lCDF = Vector{Float64}(undef, model.n)
    lP = 0.0

    for ii = 1:nsim
        for i = 1:model.n

            j = deepcopy(sims[ii][:ζ][i])
            if sims[ii][:ζ][i] == 0
                lP = logcdf( Normal(sims[ii][:intercept].μ, sqrt(sims[ii][:intercept].σ2)), model.y[i] )
            else
                lP = logcdf( Normal(sims[ii][:mixcomps][j][:μ] + sims[ii][:mixcomps][j][:fx][i],
                                                               sqrt(sims[ii][:mixcomps][j][:σ2])),
                                                        model.y[i] )
            end

            lCDF[i] = deepcopy(lP)
        end

        sims[ii][:condCDF] = exp.(lCDF)
    end

    return nothing
end


text:
# Particularly useful or beautiful grammar of graphics invocations.

## Return a DataFrame with x, y column suitable for plotting a function.
#
# Args:
#  f: Function/Expression to be evaluated.
#  a: Lower bound.
#  b: Upper bound.
#  n: Number of points to evaluate the function at.
#
# Returns:
#  A data frame with "x" and "f(x)" columns.
#
function evalfunc(f::Function, a, b, n)
    xs = [x for x in a:(b - a)/n:b]
    df = DataFrame(xs, map(f, xs))
    # NOTE: 'colnames!' is the older deprecated name. 'names!' was also defined
    # but threw an error.
    try
        names!(df, [:x, symbol("f(x)")])
    catch
        colnames!(df, ["x", "f(x)"])
    end
    df
end

evalfunc(f::Expr, a, b, n) = evalfunc(eval(:(x -> $f)), a, b, n)

# Create a dataframe and mapping from a list of functions or expressions.
function datafy(fs::Array, a, b)
    df = DataFrame()
    name_levels = ASCIIString[]
    for (i, f) in enumerate(fs)
        df_i = evalfunc(f, a, b, 250)
        name = typeof(f) == Expr ? string(f) : @sprintf("f<sub>%d</sub>", i)
        try
            df_i[:f] = fill(name, size(df_i, 1))
        catch
            df_i["f"] = fill(name, size(df_i, 1))
        end
        push!(name_levels, name)
        df = vcat(df, df_i)
    end
    df[:f] = PooledDataArray(df[:f], name_levels)

    mapping = {:x => "x", :y => "f(x)"}
    if length(fs) > 1
        mapping[:color] = "f"
    end

    df, mapping
end

# A convenience plot function for quickly plotting functions or expressions.
#
# Args:
#   fs: An array in which each object is either a single argument function or an
#       expression computing some value on x.
#   a: Lower bound on x.
#   b: Upper bound on x.
#   elements: One ore more grammar elements.
#
# Returns:
#   A plot objects.
#
function plot(fs::Array, a, b, elements::ElementOrFunction...; mapping...)
    df, mappingdict = datafy(fs, a, b)
    for (k, v) in mapping
        mappingdict[k] = v
    end
    plot(df, mappingdict, Geom.line, elements...)
end

# Plot a single function.
function plot(f::Function, a, b, elements::ElementOrFunction...; mapping...)
    plot([f], a, b, elements...; mapping...)
end

# Plot a single expression.
function plot(f::Expr, a, b, elements::ElementOrFunction...; mapping...)
    plot([f], a, b, elements...; mapping...)
end

# Plot an expression from a to b.
macro plot(expr, a, b)
    quote
        plot(x -> $(expr), $(a), $(b))
    end
end

# Create a layer from a list of functions or expressions.
function layer(fs::Array, a, b)
    df, mapping = datafy(fs, a, b)
    layer(df, Stat.nil(), Geom.line; mapping...)
end

# Create a layer from a single function.
function layer(f::Function, a, b)
    layer([f], a, b)
end

# Create a layer from a single expression.
function layer(f::Expr, a, b)
    layer([f], a, b)
end

# Create a layer from an expression from a to b.
macro layer(expr, a, b)
    quote
        layer(x -> $(expr), $(a), $(b))
    end
end

# Simple heatmap plots of matrices.
#
# Args:
#   M: A matrix.
#
# Returns:
#   A plot object.
#
function spy(M::AbstractMatrix, elements::ElementOrFunction...; mapping...)
    is, js, values = findnz(M)
    df = DataFrame(i=is, j=js, value=values)
    plot(df, x="j", y="i", color="value",
         Coord.cartesian(yflip=true),
         Scale.continuous_color,
         Scale.x_discrete,
         Scale.y_discrete,
         Geom.rectbin,
         Stat.identity,
         elements...;
         mapping...)
end


text:
struct LabelArgCheck end
const LABEL_ARGCHECK = LabelArgCheck()

const LABEL_BEGIN_CHECK = Expr(:meta, :begin_optional, LABEL_ARGCHECK)
const LABEL_END_CHECK = Expr(:meta, :end_optional, LABEL_ARGCHECK)

function mark_check(code)
    # Mark a code block as check, for usage with OptionalArgChecks.jl
    Expr(:block,
        LABEL_BEGIN_CHECK,
        code,
        LABEL_END_CHECK,
     )
end

struct CheckError <: Exception
    msg::String
end
Base.showerror(io::IO, err::CheckError) = print(io, "CheckError: $(err.msg)")

abstract type AbstractCheckFlavor end
struct ArgCheckFlavor <: AbstractCheckFlavor end
struct CheckFlavor    <: AbstractCheckFlavor end

abstract type AbstractCodeFlavor end
struct CallFlavor       <: AbstractCodeFlavor end
struct ComparisonFlavor <: AbstractCodeFlavor end
struct FallbackFlavor   <: AbstractCodeFlavor end

struct Checker
    code
    checkflavor::AbstractCheckFlavor
    codeflavor::AbstractCodeFlavor
    options
end

abstract type AbstractErrorInfo end
struct CallErrorInfo <: AbstractErrorInfo
    code
    checkflavor::AbstractCheckFlavor
    argument_expressions::Vector
    argument_values::Vector
    options::Tuple
end
struct ComparisonErrorInfo <: AbstractErrorInfo
    code
    checkflavor::AbstractCheckFlavor
    argument_expressions::Vector
    argument_values::Vector
    options::Tuple
end
struct FallbackErrorInfo <: AbstractErrorInfo
    code
    checkflavor::AbstractCheckFlavor
    options::Tuple
end

"""
    @argcheck

Check invariants on function arguments and produce
a nice exception message if they are violated.
Usage is as follows:
```Julia
function myfunction(k,n,A,B)
    @argcheck k > n
    @argcheck size(A) == size(B) DimensionMismatch
    @argcheck det(A) < 0 DomainError()
    # doit
end
```
See also [`@check`](@ref).
"""
macro argcheck(ex, options...)
    check(ex, ArgCheckFlavor(), options...)
end

"""
    @check

Check that a condition holds
and produce a nice exception message, if it does not.
Usage is as follows:
```Julia
@check k > n
@check size(A) == size(B) DimensionMismatch
@check det(A) < 0 DomainError()
```
See also [`@argcheck`](@ref).
"""
macro check(ex, options...)
    check(ex, CheckFlavor(), options...)
end

function check(ex, checkflavor, options...)
    codeflavor = if isexpr(ex, :comparison)
        ComparisonFlavor()
    elseif isexpr(ex, :call)
        CallFlavor()
    else
        FallbackFlavor()
    end
    checker = Checker(ex, checkflavor, codeflavor, options)
    inner = check(checker, codeflavor)
    mark_check(inner)
end

function check(c, ::FallbackFlavor)
    info = Expr(:call, :FallbackErrorInfo,
                QuoteNode(c.code),
                c.checkflavor,
                Expr(:tuple, esc.(c.options)...))

    condition = esc(c.code)
    expr_error_block(info, condition)
end

const SPLAT = Symbol("...")

function analyze_call_arg(expr, isparameter)
    if isexpr(expr, :kw)
        (kind=:kw, expr=expr.args[2], key=expr.args[1], symbol=gensym("kw"), isparameter=isparameter)
    elseif isexpr(expr, SPLAT)
        (kind=:splat, expr=expr.args[1], symbol=gensym("splat"), isparameter=isparameter)
    elseif isparameter
        if expr isa Symbol
            # e.g. ;foo
            (kind=:kw, expr=expr, key=expr, symbol=gensym("arg"), isparameter=isparameter)
        else
            # e.g. ;foo => bar
            (kind=:ghost, expr=expr, isparameter=isparameter)
        end
    else
        (kind=:arg, expr=expr, symbol=gensym("arg"), isparameter=isparameter)
    end
end

function analyze_call(expr)
    @assert Meta.isexpr(expr, :call)
    args = []
    kw_args = []
    for item in expr.args[2:end]
        if Meta.isexpr(item, :parameters)
            append!(kw_args, analyze_call_arg.(item.args, true))
        else
            push!(args, analyze_call_arg(item, false))
        end
    end
    pushfirst!(args, (kind=:calle, expr=expr.args[1], symbol=gensym("calle")))
    append!(args, kw_args)
    return (args=args, expr=expr)
end

function build_call(ana)
    @assert :args in propertynames(ana)
    @assert :expr in propertynames(ana)
    arglist = []
    parameterlist = []
    for item in ana.args[2:end]
        list =  item.isparameter ? parameterlist : arglist
        if item.kind == :arg
            push!(list, item.symbol)
        elseif item.kind == :kw
            push!(list, Expr(:kw, item.key, item.symbol))
        elseif item.kind == :splat
            push!(list, Expr(SPLAT, item.symbol))
        elseif item.kind == :ghost
            push!(list, esc(item.expr))
        else
            error("Unreachable $item $ana")
        end
    end
    callargs = copy(arglist)
    if !isempty(parameterlist)
        pushfirst!(callargs, Expr(:parameters, parameterlist...))
    end
    item = ana.args[1]
    @assert item.kind == :calle
    f = item.symbol
    return Expr(:call, f, callargs...)
end

function check(c, ::CallFlavor)
    ana = analyze_call(c.code)
    variables = []
    argument_expressions = []
    assignments = []
    for item in ana.args
        if item.kind == :ghost
            nothing
        else
            push!(variables, item.symbol)
            push!(argument_expressions, item.expr)
            push!(assignments, :($(item.symbol) = $(esc(item.expr))))
        end
    end
    condition = build_call(ana)
    info = Expr(:call, :CallErrorInfo,
                QuoteNode(c.code),
                c.checkflavor,
                QuoteNode(argument_expressions),
                Expr(:vect, variables...),
                Expr(:tuple, esc.(c.options)...))

    expr_error_block(info, condition, assignments...)
end

function check(c::Checker, ::ComparisonFlavor)
    exprs = c.code.args[1:2:end]
    ops = c.code.args[2:2:end]
    variables = [gensym() for _ in 1:length(exprs)]
    ret = []
    rhs = exprs[1]
    vrhs = variables[1]
    assignment = Expr(:(=), vrhs, esc(rhs))
    push!(ret, assignment)
    for i in eachindex(ops)
        op = ops[i]
        lhs = rhs
        vlhs = vrhs
        rhs = exprs[i+1]
        vrhs = variables[i+1]
        assignment = Expr(:(=), vrhs, esc(rhs))
        condition = Expr(:call, esc(op), vlhs, vrhs)
        code = Expr(:call, op, lhs, rhs)
        info = Expr(:call, :ComparisonErrorInfo,
                    QuoteNode(code),
                    c.checkflavor,
                    QuoteNode([lhs, rhs]),
                    Expr(:vect, vlhs, vrhs),
                    Expr(:tuple, esc.(c.options)...))

        reti = expr_error_block(info, condition, assignment)
        append!(ret, reti.args)
    end
    Expr(:block, ret...)
end

function expr_error_block(info, condition, preamble...)
    quote
        $(preamble...)
        if $condition
            nothing
        else
            throw_check_error($info)
        end
    end |> Base.remove_linenums!
end

@noinline function throw_check_error(info...)::Union{}
    # the compiler sometimes deoptimizes functions with inline errors
    @nospecialize
    err = build_error(info...)
    throw(err)
end

default_exception_type(::ArgCheckFlavor) = ArgumentError
default_exception_type(::CheckFlavor) = CheckError

function build_error(info)
    build_error(info, info.checkflavor, info.options...)
end
function build_error(info, checkflavor, msg::AbstractString)
    E = default_exception_type(checkflavor)
    E(msg)
end

function build_error(info, checkflavor,
                     T::Type{<:Exception}=default_exception_type(checkflavor))
    msg = error_message(info)
    T(msg)
end
function build_error(info, checkflavor, err::Exception)
    err
end

error_message(info::FallbackErrorInfo) = "$(info.code) must hold."
error_message(info::CallErrorInfo) = fancy_error_message(info)
error_message(info::ComparisonErrorInfo) = fancy_error_message(info)

function pretty_string(data)
    io = IOBuffer()
    ioc = IOContext(io, :limit=>true, :compact=>true)
    show(ioc, data)
    seekstart(io)
    String(take!(io))
end
pretty_string(x::Number) = string(x)
pretty_string(s::Symbol) = string(s)
pretty_string(ex::Expr) = string(ex)

function fancy_error_message(info)
    code = info.code
    exprs = info.argument_expressions
    values = info.argument_values
    lines = String[]
    foreach(exprs, values) do ex, val
        pex = pretty_string(ex)
        pval = pretty_string(val)
        if pex != pval
            push!(lines, "$pex => $pval")
        end
    end
    firstline = if isempty(lines)
        "$code must hold."
    else
        "$code must hold. Got"
    end
    pushfirst!(lines, firstline)
    join(lines, '\n')
end


text:
#@ Model: Knapsack Model
#@ Description: Sample knapsack model in JuMP $$\max p^T x$$$$w^T x \leq k$$$$x\in\{0,1\}^5$$

using JuMP, Cbc, JSON

# Labels
items = Dict()
for i in 1:5
    items[i] = "item_$(i)"
end

#@ Input File: profit_file
#@ Description: This is a sample input file
profit_file = open("./examples/jump/knapsack/profit.json")

#@ Input Object: weight
#@ Description: Weight of items, w
weight = [2, 8, 4, 2, 5]

#@ Input Object: capacity
#@ Description: Capacity of knapsack, k
capacity = 10

#@ Helper object: profit
#@ Description: Profit of each item, p
profit = convert(Array{Int64,1}, JSON.parse(profit_file))

#@ Solver: solver
solver = Cbc.Optimizer

#@ Problem: model
model = Model(optimizer_with_attributes(solver, "print_level" => 1))

#@ Variable: x
#@ Description: 1 if an item is included in knapsack, 0 otherwise$$x\in\{0,1\}^5$$
#@ Labels: items
@variable(model, x[1:5], Bin)

#@ Function: objectivefn
#@ Description: objective function $$p^T x$$
@expression(model, objectivefn, profit' * x)

#@ Constraint: cap_constr
#@ Description: Items placed in knapsack may not exceed capacity $$w^T x \leq k$$
@constraint(model,
            cap_constr,
            convert(Array{Int64,1}, weight)' * x <= capacity)

# Execution:
@objective(model, Max, objectivefn)
JuMP.optimize!(model)

#@ Output object: sum_x
#@ Description: This is a sample output objet
sum_x = sum(JuMP.value.(x))

#@ Output file: output.txt
#@ Description: This is a sample output file
open("output.txt", "w") do f
    xval = JuMP.value.(x)
    for i=1:5
        write(f, "$(xval[i])\n")
    end
end

max_stars_repo_name: Fuinn/mos-examples

max_issues_repo_name: Fuinn/mos-examples

max_forks_repo_name: Fuinn/mos-examples


text:
export transition, custom_transition,
    AbstractLine, GenericLine, BroadLine, NarrowLine, BroadBaseLine, AsymmTailLine, MultiCompLine, LineComponent

const transitions_db = DataFrame()

function load_transitions()
    global transitions_db
    (nrow(transitions_db) > 0)  &&  (return nothing)

    d, c = csvread(joinpath(@__DIR__, "atomic_line_list.vac"), '|', commentchar='#', header_exists=true)
    out = DataFrame(collect(d), Symbol.(strip.(c)))
    for i in 1:ncol(out)
        if isa(out[1, i], String)
            out[:, i] .= strip.(out[:,i])
        end
    end
    delete!(out, findall((out.tid .== "")  .|
                         (out.Label .== "")))

    out[!, :tid] .= Symbol.(out.tid)

    append!(transitions_db, out)
    return nothing
end

function transition(tid::Symbol)
    load_transitions()
    i = findall(transitions_db.tid .== tid)
    if length(i) == 0
        error("No known transition with ID: $tid")
    end
    return transitions_db[i, :]
end

function custom_transition_default_tid(λ::Float64)
    tid = "T" * string(λ)
    tid = join(split(tid, "."), "p")
    return Symbol(tid)
end

function custom_transition(λ_vac_ang::Float64; tid::Union{Symbol, Nothing}=nothing)
    load_transitions()
    isnothing(tid)  &&  (tid = custom_transition_default_tid(λ_vac_ang))
    i = findall(transitions_db.tid .== tid)
    if length(i) == 0
        push!(transitions_db, [tid, string(tid), λ_vac_ang, fill("", ncol(transitions_db)-3)...])
    end
    return tid
end

transition(λ_vac_ang::Float64) = transition(custom_transition_default_tid(λ_vac_ang))

# Line type descriptors associated to a specific transition
abstract type AbstractLine end

macro define_line(name, suffix, group)
    return esc(:(
        struct $name <: AbstractLine; tid::Symbol; end;
        suffix(::$(name), multicomp::Bool) = multicomp  ?  $(QuoteNode(suffix))  :  Symbol("");
        group( ::$(name)) = $(QuoteNode(group));
    ))
end

@define_line GenericLine   _gen GenericLines
@define_line BroadLine     _br  BroadLines
@define_line NarrowLine    _na  NarrowLines
@define_line BroadBaseLine _bb  BroadBaseLines

struct AsymmTailLine <: AbstractLine
    tid::Symbol
    side::Symbol

    function AsymmTailLine(tid::Symbol, side::Symbol)
        @assert side in (:red, :blue)
        new(tid, side)
    end
end
suffix(line::AsymmTailLine, multicomp::Bool) = line.side == :blue ?  :_bw  :  :_rw
group(::AsymmTailLine) = :AsymmTailLines

struct MultiCompLine <: AbstractLine
    tid::Symbol
    types::Vector{Type}
end

struct LineComponent
    line::AbstractLine
    comp::AbstractComponent
    multicomp::Bool
end

function collect_LineComponent(source::QSO)
    out = OrderedDict{Symbol, LineComponent}()
    for line in known_spectral_lines(source)
        if isa(line, MultiCompLine)
            for t in line.types
                nn = Symbol(line.tid, suffix(t(line.tid), true))
                (nn in source.options[:skip_lines])  &&  continue
                @assert !haskey(out, nn)
                lc = LineComponent(source, t(line.tid), true)
                out[nn] = lc
            end
        else
            nn = Symbol(line.tid, suffix(line, false))
            (nn in source.options[:skip_lines])  &&  continue
            @assert !haskey(out, nn)
            out[nn] = LineComponent(source, line, false)
        end
    end
    return out
end

LineComponent(source::QSO, line::AbstractLine, multicomp::Bool) =
    error("No constructor method has been defined for LineComponent($(typeof(source)), $(typeof(line)), ::Bool)")


text:
module _XFuzzEn
export XFuzzEn
using Statistics: mean, std
    """
        XFuzz, Ps1, Ps2 = XFuzzEn(Sig)

    Returns the cross-fuzzy entropy estimates (`XFuzz`) and the average fuzzy
    distances (m:Ps1, m+1:Ps2) for m = [1,2] estimated for the data sequences
    contained in 'Sig', using the default parameters: embedding dimension = 2,
    time delay = 1, fuzzy function (Fx) = 'default',
    fuzzy function parameters (r) = [0.2, 2], logarithm = natural

        XFuzz, Ps1, Ps2 = XFuzzEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1, r::Union{Real,Tuple{Real,Real}}=(.2,2), Fx::String="default", Logx::Real=exp(1))

    Returns the cross-fuzzy entropy estimates (`XFuzz`) for dimensions = [1,...,m]
    estimated for the data sequences in `Sig` using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Embedding Dimension, a positive integer   [default: 2]    \n
    `tau`   - Time Delay, a positive integer            [default: 1]    \n
    `Fx`    - Fuzzy funtion name, one of the following:
              {`default`,`sigmoid`,`modsampen`,`gudermannian`,`linear`}\n
    `r`     - Fuzzy function parameters, a scalar or a 2 element tuple
              of positive values. The `r` parameters for each fuzzy
              function are defined as follows:\n
              sigmoid:       r(1) = divisor of the exponential argument
                             r(2) = value subtracted from argument (pre-division)
              modsampen:     r(1) = divisor of the exponential argument
                             r(2) = value subtracted from argument (pre-division)
              default:       r(1) = divisor of the exponential argument
                             r(2) = argument exponent (pre-division)
              gudermannian:  r    = a scalar whose value is the numerator of
                                    argument to gudermannian function:
                                    GD(x) = atan(tanh(`r`/x)).
              linear:        r    = an integer value. When `r` = 0, the
                                    argument of the exponential function is
                                    normalised between [0 1]. When `r` = 1,
                                    the minimuum value of the exponential
                                    argument is set to 0.     \n
    `Logx`  - Logarithm base, a positive scalar  \n

    For further information on the 'keyword' arguments, see the EntropyHub guide.

    # See also `FuzzEn`, `XSampEn`, `XApEn`, `FuzzEn2D`, `XMSEn`, `MSEn`

    # References:
        [1] Hong-Bo Xie, et al.,
            "Cross-fuzzy entropy: A new method to test pattern synchrony of
            bivariate time series."
            Information Sciences
            180.9 (2010): 1715-1724.

    """
    function XFuzzEn(Sig::AbstractArray{T,2} where T<:Real; m::Int=2, tau::Int=1,
         r::Union{Real,Tuple{Real,Real}}=(.2,2), Fx::String="default", Logx::Real=exp(1))

    (size(Sig,2) > size(Sig,1)) ? Sig = transpose(Sig) : nothing
    N = size(Sig,1)
    (N>10 && size(Sig,2)==2) ? nothing :  error("Sig:   must be a 2-columns matrix")
    (m > 0) ? nothing : error("m:     must be an integer > 0")
    (tau>0) ? nothing : error("tau:   must be an integer > 0")
    (minimum(r)>=0 && length(r)<=2) ? nothing :
        error("r:  must be a scalar or 2 element vector of positive values")
    (lowercase(Fx) in ["default","sigmoid","modsampen","gudermannian","linear"]) ?
        nothing : error("Fx:    must be one of the following strings -
        'default', 'sigmoid', 'modsampen', 'gudermannian', 'linear'")
    (Logx>0) ? nothing : error("Logx:   must be a positive number > 0")

    if length(r) == 2 && lowercase(Fx)=="linear"
        r = 0;
        print("Multiple values for r entered. Default value (0) used.\n")
    elseif length(r) == 2 && lowercase(Fx)=="gudermannian"
        r = r[1]
        print("Multiple values for r entered. First value used.\n")
    end

    S1 = Sig[:,1]; S2 = Sig[:,2]
    m += 1
    Fun = getfield(_XFuzzEn,Symbol(lowercase(Fx)))
    Sx1 = zeros(N,m)
    Sx2 = zeros(N,m)
    for k = 1:m
        Sx1[1:N-(k-1)*tau,k] = S1[1 + (k-1)*tau:N]
        Sx2[1:N-(k-1)*tau,k] = S2[1 + (k-1)*tau:N]
    end

    Ps1 = zeros(m)
    Ps2 = zeros(m-1)
    Ps1[1] = 1
    for k = 2:m
        N1 = N - k*tau
        N2 = N - (k-1)*tau
        A = Sx1[1:N2,1:k] .- mean(Sx1[1:N2,1:k],dims=2)
        B = Sx2[1:N2,1:k] .- mean(Sx2[1:N2,1:k],dims=2)
        d2 = zeros(N2,N2)
        for p = 1:N2
            Mu2 = maximum(abs.(transpose(A[p,:]) .- B),dims=2)
            d2[p,:] = Fun(Mu2,r)
        end
        Ps1[k] = mean(d2[1:N1,1:N1])
        Ps2[k-1] = mean(d2)
    end

    XFuzz = log.(Logx, Ps1[1:end-1]) .- log.(Logx, Ps2)
    return XFuzz, Ps1, Ps2
    end

    function sigmoid(x,r)
        if length(r) == 1
            error("When Fx = 'Sigmoid', r must be a two-element vector.")
        end
        y = inv.(1 .+ exp.((x.-r[2])/r[1]))
        return y
    end

    function modsampen(x,r)
        if length(r) == 1
            error("When Fx = 'Modsampen', r must be a two-element vector.")
        end
        y = inv.(1 .+ exp.((x.-r[2])/r[1]))
        return y
    end

    function default(x,r)
        if length(r) == 1
            error("When Fx = 'Default', r must be a two-element vector.")
        end
        y = exp.(-(x.^r[2])/r[1])
        return y
    end

    function gudermannian(x,r)
        if r <= 0
            error("When Fx = 'Gudermannian', r must be a scalar > 0.")
        end
        y = atan.(tanh.(r[1]./x))
        y ./= maximum(y)
        return y
    end

    function linear(x,r)
        if r == 0 && length(x)>1
            y = exp.(-(x .- minimum(x))/(maximum(x)-minimum(x)))
        elseif r == 1
            y = exp.(-(x .- minimum(x)))
        elseif r == 0 && length(x)==1
            y = [0]
        else
            error("When Fx = 'Linear', r must be 0 or 1.")
        end
        return y
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

     http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
addprocs(3)

using SharedDicts
using Base.Test

@everywhere begin
    function check(DS::SharedDicts.SharedDict, idx::Int)
        Base.Test.@test DS[idx] == idx * 2
        nothing
    end

    function checkiter(DS::SharedDicts.SharedDict)
        Base.Test.@test !isempty(DS)
        for idx in 1:length(DS)
            Base.Test.@test DS[idx] == idx * 2
        end
        nothing
    end
end

function dotests()
    N = 1000
    D = Dict{Int, Int}()

    for idx in 1:N
        D[idx] = idx * 2
    end

    DS = share(D)

    @sync @parallel for idx in 1:N
        check(DS, idx)
    end
    for w in workers()
        remotecall_fetch(checkiter, w, DS)
    end
end

dotests()


text:
#GracePlot functions to control Grace/xmgrace
#-------------------------------------------------------------------------------

#==Type definitions
===============================================================================#

#Map an individual attribute to a setter function:
const AttributeFunctionMap = Dict{Symbol, Function}

#Map an "AttributeList" element to a setter function:
#NOTE: Unlike individual attributes, typed obects do not need to be "set"
#      using keyword arguments
#TODO: Find way to restrict Dict to DataTypes inherited from AttributeList
#const AttributeListFunctionMap = Dict{DataType, Function}
#Use IdDict: Compatibility with precompile.
const AttributeListFunctionMap = IdDict #{DataType, Function}

#Map attribute fields to grace commands:
const AttributeCmdMap = Dict{Symbol, String}

#==Helper functions
===============================================================================#

#Copy in only attributes that are not "nothing" (thus "new"):
function copynew!(dest::T, newlist::T) where T<:AttributeList
	for attrib in fieldnames(newlist)
		v = getfield(newlist, attrib)

		if v != nothing
			setfield!(dest, attrib, v)
		end
	end
end

function setattrib(p::Plot, cmd::String, value::Any) #Catchall
	sendcmd(p, "$cmd $value")
end
function setattrib(p::Plot, cmd::String, value::Nothing)
	#Nothing to do
end
function setattrib(p::Plot, cmd::String, value::String)
	sendcmd(p, "$cmd \"$value\"") #Add quotes around string
end
function setattrib(p::Plot, cmd::String, value::GraceConstLitteral)
	sendcmd(p, "$cmd $(value.s)") #Send associated string, unquoted
end
function setattrib(p::Plot, cmd::String, value::Tuple{Number, Number})
	v1 = value[1]; v2 = value[2]
	sendcmd(p, "$cmd $v1, $v2") #Send out tuple
end
setattrib(p::Plot, cmd::String, value::Symbol) =
	setattrib(p, cmd, graceconstmap[value])

#Set plot attributes for a given element:
#-------------------------------------------------------------------------------
function setattrib(p::Plot, fmap::AttributeCmdMap, prefix::String, data::Any)
	for attrib in fieldnames(typeof(data))
		v = getfield(data, attrib)

		if v != nothing
			subcmd = get(fmap, attrib, nothing)

			if subcmd != nothing
				setattrib(p, "$prefix$subcmd", v)
			else
				dtype = typeof(data)
				@warn("Attribute \"$attrib\" of $dtype not currently supported.")
			end
		end
	end
end

#Set graph attributes for a given element:
#-------------------------------------------------------------------------------
function setattrib(g::GraphRef, fmap::AttributeCmdMap, prefix::String, data::Any)
	setactive(g)
	return setattrib(g.plot, fmap, prefix, data)
end
function setattrib(g::GraphRef, prefix::String, value::String)
	setactive(g)
	return setattrib(g.plot, prefix, value)
end

#Set dataset attribute:
#-------------------------------------------------------------------------------
function setattrib(ds::DatasetRef, fmap::AttributeCmdMap, data::Any)
	dsid = ds.id
	setattrib(ds.graph, fmap, "S$dsid ", data)
end
function setattrib(ds::DatasetRef, cmd::String, value::String)
	dsid = ds.id
	setattrib(ds.graph, "S$dsid $cmd", value)
end

#Core algorithm for "set" interface:
#-------------------------------------------------------------------------------
function _set(obj::Any, listfnmap::AttributeListFunctionMap, fnmap::AttributeFunctionMap, args...; kwargs...)
	for value in args
		setfn = get(listfnmap, typeof(value), nothing)

		if setfn != nothing
			setfn(obj, value)
		else
			argstr = string(typeof(value))
			objtype = typeof(obj)
			@warn("Argument \"$argstr\" not recognized by \"set(::$objtype, ...)\"")
		end
	end

	for (arg, value) in kwargs
		setfn = get(fnmap, arg, nothing)

		if setfn != nothing
			setfn(obj, value)
		else
			argstr = string(arg)
			objtype = typeof(obj)
			@warn("Argument \"$argstr\" not recognized by \"set(::$objtype, ...)\"")
		end
	end
	return
end

#Core algorithm for "get" interface:
#-------------------------------------------------------------------------------
function getattrib(obj::Any, fnmap::AttributeFunctionMap, attrib::Symbol)
	getfn = get(fnmap, attrib, nothing)

	if nothing == getfn
		argstr = string(attrib)
		objtype = typeof(obj)
		throw("Argument \"$argstr\" not recognized by \"set(::$objtype, ...)\"")
	end

	return getfn(obj)
end

#Core algorithm for "addannotation" interface:
#(Just use set in the background... it does what you would want...)
#-------------------------------------------------------------------------------
function _addannotation(obj::Any, listfnmap::AttributeListFunctionMap, fnmap::AttributeFunctionMap, args...; kwargs...)
	return _set(obj::Any, listfnmap, fnmap, args...; kwargs...)
end

#==Plot-level functionality
===============================================================================#

#-------------------------------------------------------------------------------
updateall(p::Plot) = sendcmd(p, "UPDATEALL")
function redraw(p::Plot; update=true)
	if !p.guimode; return; end #Causes issues when !guimode
	result = sendcmd(p, "REDRAW")
	if update; updateall(p); end
	return result
end

#-------------------------------------------------------------------------------
function clearall(p::Plot; update::Bool=true, killdata::Bool=true)
	#Delete graphs in reverse order
	#(Grace wants to keep around lower numbered graphs)
	for gidx in (length(p.graphs):-1:1) .- 1
		clearall(graph(p, gidx), update=false, killdata=killdata)
	end

	#Trash old graph info:
	p.graphs = Graph[]

	#Sync up UI, if desired:
	if update
		updateall(p)
	end
end

#-------------------------------------------------------------------------------
function setpagesize(p::Plot, a::CanvasAttributes)
	p.canvas = a
	width = round(Int, val(TPoint(a.width)))
	height = round(Int, val(TPoint(a.height)))
	sendcmd(p, "PAGE SIZE $width, $height")
	#Will also stretch entire plot if RESIZE is used instead:
	#sendcmd(p, "PAGE RESIZE $width, $height")
	return a
end

#Return cached canvas width/height
#-------------------------------------------------------------------------------
getwcanvas(p::Plot) = p.canvas.width
gethcanvas(p::Plot) = p.canvas.height

#Obtain "VIEW" width/height (given cached canvas width/height)
#-------------------------------------------------------------------------------
function getaspectratio(p::Plot)
	w = val(TPoint(p.canvas.width))
	h = val(TPoint(p.canvas.height))
	return w/h
end

#Return canvas width in normalized "view" coordinates:
function getwview(p::Plot)
	ar=getaspectratio(p)
	return (ar>1) ? ar : 1
end

#Return canvas width in normalized "view" coordinates:
function gethview(p::Plot)
	ar=getaspectratio(p)
	return (ar>1) ? 1 : (1/ar)
end

#-------------------------------------------------------------------------------
setnumcols(p::Plot, cols::Int) = (p.ncols = cols)

#-------------------------------------------------------------------------------
function arrange(p::Plot, gdim::GraphCoord; offset=0.15, hgap=0.15, vgap=0.15)
	(rows, cols) = gdim
	sendcmd(p, "ARRANGE($rows, $cols, $offset, $hgap, $vgap)")
	setnumcols(p, cols)
	newsize = rows*cols
	delta = newsize - length(p.graphs)

	for i in 1:delta
		push!(p.graphs, Graph())
	end

	if delta < 0
		resize!(p.graphs, newsize)
		if p.activegraph >= newsize
			setactive(graph(p, newsize-1))
		end
	end
end

#Add new graph to a plot
#NOTE: update = true sends UPDATEALL command to avoid having the GUI out of
#      sync with the plot itself. User expected to call updateall manually when
#      using update=false
#-------------------------------------------------------------------------------
function add(p::Plot, args...; update=true, kwargs...)
	gidx = length(p.graphs)
	push!(p.graphs, Graph())
	g = graph(p, gidx)
	setenable(g, true)
#	setactive(g)
	if update; updateall(p); end
	set(g, args...; kwargs...)

	return g
end

#Add plot annotation:
#-------------------------------------------------------------------------------
function addtext(p::Plot, a::TextAttributes)
	prefix = "    STRING"
	sendcmd(p, "WITH STRING")
	sendcmd(p, "$prefix ON")
	setattrib(p::Plot, "$prefix LOCTYPE", a.loctype)
	setattrib(p::Plot, "$prefix", a.loc)
	setattrib(p::Plot, "$prefix FONT", a.font)
	setattrib(p::Plot, "$prefix CHAR SIZE", a.size)
	setattrib(p::Plot, "$prefix COLOR", a.color)
	setattrib(p::Plot, "$prefix JUST", a.just)
	setattrib(p::Plot, "$prefix DEF", a.value) #This has to go last, for some reason
	#...so cannot use attribcmdmap method....
end

#-------------------------------------------------------------------------------
const defaults_attribcmdmap = AttributeCmdMap(
	:linewidth  => "LINEWIDTH",
	:linestyle  => "LINESTYLE",
	:color      => "COLOR",
	:pattern    => "PATTERN",
	:font       => "FONT",
	:charsize   => "CHAR SIZE",
	:symbolsize => "SYMBOL SIZE",
	:sformat    => "SFORMAT",
)

setdefaults(p::Plot, a::DefaultAttributes) = setattrib(p, defaults_attribcmdmap, "DEFAULT ", a)

#==Graph-level functionality
===============================================================================#

#-------------------------------------------------------------------------------
function Base.kill(g::GraphRef)
	gidx = graphindex(g)
	sendcmd(g.plot, "KILL G$gidx")
end

#-------------------------------------------------------------------------------
function setactive(g::GraphRef)
	gidx = graphindex(g)
	if g.plot.activegraph != gidx
		g.plot.activegraph = gidx
		sendcmd(g.plot, "WITH G$gidx")
	end
	return gidx
end

#-------------------------------------------------------------------------------
function setfocus(g::GraphRef)
	gidx = graphindex(g)
	sendcmd(g.plot, "FOCUS G$gidx")
end

#Convenience functions: Enables the use of set(::Plot, ...) interface:
#(Plot argument is otherwise redundant)
#-------------------------------------------------------------------------------
function setactive(p::Plot, g::GraphRef)
	_ensure(p==g.plot,
		ArgumentError("setactive: GraphRef does not match Plot to control."))
	setactive(g)
end

function setfocus(p::Plot, g::GraphRef)
	_ensure(p==g.plot,
		ArgumentError("setfocus: GraphRef does not match Plot to control."))
	setfocus(g)
end
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
function setenable(g::GraphRef, value::Bool)
	v = graceconstmap[value ? (:on) : (:off)]
	gidx = graphindex(g)
	setattrib(g.plot, "G$gidx", v)
end

#-------------------------------------------------------------------------------
function clearall(g::GraphRef; update::Bool=true, killdata::Bool=true)
	setenable(g, false)
	gdata = graphdata(g)

	if killdata
		#Delete datasets in reverse order
		#(Grace wants to keep around lower numbered datasets)
		for idx in (gdata.datasetcount-1:-1:0)
			killdataset(g, idx)
		end
	end

	#Sync up UI, if desired:
	if update
		updateall(p)
	end
end

#NOTE: AUTOSCALE X/Y does not seem to work...
#-------------------------------------------------------------------------------
function _autofit(g::GraphRef, x::Bool=false, y::Bool=false)
	cmd = "AUTOSCALE"
	if x && y
		; #Send command by itself
	elseif x
		cmd *= " XAXES"
	elseif y
		cmd *= " YAXES"
	else
		return 0
	end

	setactive(g)
	sendcmd(g.plot, cmd)
end
_autofit(g::GraphRef, x::Nothing, y::Nothing) = _autofit(g, true, true)
_autofit(g::GraphRef, x::Nothing, y::Bool) = _autofit(g, false, y)
_autofit(g::GraphRef, x::Bool, y::Nothing) = _autofit(g, x, false)
autofit(g::GraphRef; x=nothing, y=nothing) = _autofit(g, x, y)

#-------------------------------------------------------------------------------
const limits_attribcmdmap = AttributeCmdMap(
	:xmin => "XMIN",
	:xmax => "XMAX",
	:ymin => "YMIN",
	:ymax => "YMAX",
)

setview(g::GraphRef, a::CartesianLimAttributes) = setattrib(g, limits_attribcmdmap, "VIEW ", a)

#-------------------------------------------------------------------------------
const title_attribcmdmap = AttributeCmdMap(
	:value => "",
	:font  => "FONT",
	:size  => "SIZE",
	:color => "COLOR",
)

settitle(g::GraphRef, a::TextAttributes) = setattrib(g, title_attribcmdmap, "TITLE ", a)
setsubtitle(g::GraphRef, a::TextAttributes) = setattrib(g, title_attribcmdmap, "SUBTITLE ", a)
settitle(g::GraphRef, title::String) = settitle(g, text(title))
setsubtitle(g::GraphRef, title::String) = setsubtitle(g, text(title))

#-------------------------------------------------------------------------------
const label_attribcmdmap = AttributeCmdMap(
	:value => "",
	:font  => "FONT",
	:size  => "CHAR SIZE",
	:color => "COLOR",
)

setxlabel(g::GraphRef, a::TextAttributes) = setattrib(g, title_attribcmdmap, "XAXIS LABEL ", a)
setylabel(g::GraphRef, a::TextAttributes) = setattrib(g, title_attribcmdmap, "YAXIS LABEL ", a)
setxlabel(g::GraphRef, label::String) = setxlabel(g, text(label))
setylabel(g::GraphRef, label::String) = setylabel(g, text(label))

#-------------------------------------------------------------------------------
const frameline_attribcmdmap = AttributeCmdMap(
	:style => "LINESTYLE",
	:width => "LINEWIDTH",
)
setframeline(g::GraphRef, a::LineAttributes) = setattrib(g, frameline_attribcmdmap, "FRAME ", a)

const axes_attribcmdmap = AttributeCmdMap(
	:xmin => "WORLD XMIN", :xmax => "WORLD XMAX",
	:ymin => "WORLD YMIN", :ymax => "WORLD YMAX",
	:xscale  => "XAXES SCALE",
	:yscale  => "YAXES SCALE",
	:invertx => "XAXES INVERT",
	:inverty => "YAXES INVERT",
)
setaxes(g::GraphRef, a::AxesAttributes) = setattrib(g, axes_attribcmdmap, "", a)

const legend_attribcmdmap = AttributeCmdMap(
	:loctype   => "LOCTYPE",
	:loc       => "",
	:display   => "",
	:font      => "FONT",
	:charsize  => "CHAR SIZE",
	:color     => "COLOR",
	:length    => "LENGTH",
	:invert    => "INVERT",
	:hgap      => "HGAP",
	:vgap      => "VGAP",
	:boxcolor       => "BOX COLOR",
	:boxpattern     => "BOX PATTERN",
	:boxlinewidth   => "BOX LINEWIDTH",
	:boxlinestyle   => "BOX LINESTYLE",
	:boxfillcolor   => "BOX FILL COLOR",
	:boxfillpattern => "BOX FILL PATTERN",
)
setlegend(g::GraphRef, a::LegendAttributes) = setattrib(g, legend_attribcmdmap, "LEGEND ", a)

#==Dataset-level functionality
===============================================================================#

#Sadly, appears to be no way to kill data without linestyle, etc...
#-------------------------------------------------------------------------------
function killdataset(g::GraphRef, index::Int)
	gidx = graphindex(g)
	sendcmd(g.plot, "KILL G$gidx.S$index")
end

#-------------------------------------------------------------------------------
function add(g::GraphRef, x::DataVec, y::DataVec, args...; kwargs...)
	_ensure(length(x) == length(y),
		ArgumentError("GracePlot.add(): x & y vlengths must match."))
	p = g.plot
	gidx = graphindex(g)
	gdata = graphdata(g)
	dsid = gdata.datasetcount
	gdata.datasetcount += 1
	prefix = "G$(gidx).S$dsid"

	sendcmd(p, "$prefix TYPE XY")
	for (_x, _y) in zip(x, y)
		sendcmd(p, "$prefix POINT $_x, $_y")
	end

	ds = DatasetRef(g, dsid)
	set(ds, args...; kwargs...)
	return ds
end

#-------------------------------------------------------------------------------
const dsline_attribcmdmap = AttributeCmdMap(
	:_type => "LINE TYPE",
	:style => "LINE LINESTYLE",
	:width => "LINE LINEWIDTH",
	:color => "LINE COLOR",
)
setline(ds::DatasetRef, a::LineAttributes) = setattrib(ds, dsline_attribcmdmap, a)

#-------------------------------------------------------------------------------
const glyph_attribcmdmap = AttributeCmdMap(
	:shape       => "SYMBOL",
	:size        => "SYMBOL SIZE",
	:color       => "SYMBOL COLOR",
	:pattern     => "SYMBOL PATTERN",
	:fillcolor   => "SYMBOL FILL COLOR",
	:fillpattern => "SYMBOL FILL PATTERN",
	:linewidth   => "SYMBOL LINEWIDTH",
	:linestyle   => "SYMBOL LINESTYLE",
	:char        => "SYMBOL CHAR",
	:charfont    => "SYMBOL CHAR FONT",
	:skip        => "SYMBOL SKIP",
)
setglyph(ds::DatasetRef, a::GlyphAttributes) = setattrib(ds, glyph_attribcmdmap, a)

setdatasetid(ds::DatasetRef, id::String) = setattrib(ds, "LEGEND", id)

#==Define cleaner "set" interface (minimize # of "export"-ed functions)
===============================================================================#

#-------------------------------------------------------------------------------
const empty_listfnmap = AttributeListFunctionMap()
const empty_fnmap = AttributeFunctionMap()

#-------------------------------------------------------------------------------
const setplot_listfnmap = AttributeListFunctionMap(
	CanvasAttributes  => setpagesize,
	DefaultAttributes => setdefaults,
)
const setplot_fnmap = AttributeFunctionMap(
	:ncols  => setnumcols,
	:active => setactive,
	:focus  => setfocus,
)
set(p::Plot, args...; kwargs...) = _set(p, setplot_listfnmap, setplot_fnmap, args...; kwargs...)

#-------------------------------------------------------------------------------
const setgraph_listfnmap = AttributeListFunctionMap(
	AxesAttributes   => setaxes,
	LegendAttributes => setlegend,
)
const setgraph_fnmap = AttributeFunctionMap(
	:enable    => setenable,
	:view      => setview,
	:title     => settitle,
	:subtitle  => setsubtitle,
	:xlabel    => setxlabel,
	:ylabel    => setylabel,
	:frameline => setframeline,
)
set(g::GraphRef, args...; kwargs...) = _set(g, setgraph_listfnmap, setgraph_fnmap, args...; kwargs...)

#-------------------------------------------------------------------------------
const setds_listfnmap = AttributeListFunctionMap(
	LineAttributes  => setline,
	GlyphAttributes => setglyph,
)
const setds_fnmap = AttributeFunctionMap(
	:id => setdatasetid,
)
set(ds::DatasetRef, args...; kwargs...) = _set(ds, setds_listfnmap, setds_fnmap, args...; kwargs...)

#==Define cleaner "get" interface (minimize # of "export"-ed functions)
===============================================================================#
const getplot_fnmap = AttributeFunctionMap(
	:wcanvas  => getwcanvas,
	:hcanvas  => gethcanvas,
	:wview    => getwview,
	:hview    => gethview,
)
Base.get(p::Plot, attrib::Symbol) = getattrib(p, getplot_fnmap, attrib)

#==Define cleaner "addannotation" interface (minimize # of "export"-ed functions)
===============================================================================#

#-------------------------------------------------------------------------------
const addannot_listfnmap = AttributeListFunctionMap(
	TextAttributes  => addtext,
)
addannotation(p::Plot, args...; kwargs...) = _addannotation(p, addannot_listfnmap, empty_fnmap, args...; kwargs...)

#Last line


text:
@testset "3.longest-substring-without-repeating-characters.jl" begin
    @test length_of_longest_substring("abcabcbb") == 3
    @test length_of_longest_substring("b") == 1
    @test length_of_longest_substring("bbbb") == 1
    @test length_of_longest_substring("pwwkew") == 3
    @test length_of_longest_substring("") == 0
end


text:
@testset "DataValueArray: Map" begin

    using Test
    using DataValues

# create m random arrays each with N dimensions and length dims[i] along
# dimension i for i=1:N
    N = rand(2:5)
    dims = Int[ rand(3:8) for i in 1:N ]
    m = rand(4:6)
    n = rand(1:3)
    As = Array[ [rand(dims...) for i in 1:n] ; [rand(Int, dims...) for i in n + 1:m] ]

    Ms = [ rand(Bool, dims...) for i in 1:m ]
    Xs = Array{DataValueArray,1}()
    for i in 1:m
        push!(Xs, DataValueArray(As[i]))
    end
    Ys = Array{DataValueArray,1}()
    for i in 1:m
        push!(Ys, DataValueArray(As[i], Ms[i]))
    end

    C = Array{Float64}(dims...)
    Z = DataValueArray{Float64}(dims...)

    R = map(|, Ms...)

    f(x...) = sum(x)
    dests = (C, Z)

# 1 arg
    for (args, masks) in (((As, Xs), fill((), m)), ((As, Ys), [ (Ms[i],) for i in 1:m ]))
        for i in 1:m
        # map!
            map!(f, args[1][i], args[1][i]) # map!(f, As[i], As[i])
            map!(f, args[2][i], args[2][i]) # map!(f, Xs[i], Xs[i])
            @test isequal(args[2][i], DataValueArray(args[1][i], masks[i]...))
        # map
            A = map(f, args[1][i])
            X = map(f, args[2][i])
            @test isequal(X, DataValueArray(A, masks[i]...))
        end
    end
# 2 arg
    i, j = rand(1:m), rand(1:m)
    S = map(|, Ms[i], Ms[j])
    for (args, mask) in (((As, Xs), ()), ((As, Ys), (S,)))
    # map!
        map!(f, dests[1], args[1][i], args[1][j])
        map!(f, dests[2], args[2][i], args[2][j])
        @test isequal(dests[2], DataValueArray(dests[1], mask...))
    # map
        map(f, args[1][i], args[1][j])
        map(f, args[2][i], args[2][j])
        @test isequal(dests[2], DataValueArray(dests[1], mask...))
    end
# n arg
    for (args, mask) in (((As, Xs), ()), ((As, Ys), (R,)))
    # map!
        map!(f, dests[1], args[1]...)
        map!(f, dests[2], args[2]...)
        @test isequal(dests[2], DataValueArray(dests[1], mask...))
    # map
        map(f, args[1]...)
        map(f, args[2]...)
        @test isequal(dests[2], DataValueArray(dests[1], mask...))
    end

# test map over empty DataValueArrays
    X = DataValueArray(Int[])
    h1(x) = 5.0 * x
    h2(x) = x
    h2(x...) = prod(x)

    Z1 = map(h1, X)
    Z2 = map(h2, X)
    @test isempty(Z1)
    @test isa(Z1, DataValueArray{Float64})
    @test isempty(Z2)
    @test isa(Z2, DataValueArray{Int})

# if a function has no method for inner eltype of empty DataValueArray,
# result should be empty DataValueArray{Any}() for consistency with generic map()
    h3(x::Float64...) = prod(x)
    Z3 = map(h3, X)
    @test isempty(Z3)
    @test isa(Z3, DataValueArray{Any})
    Z3 = map(h3, X, X)
    @test isempty(Z3)
    @test isa(Z3, DataValueVector{Any})
    Z3 = map(h3, X, X, X)
    @test isempty(Z3)
    @test isa(Z3, DataValueArray{Any})

# test map over all null DataValueArray
    n = rand(10:100)
    Ys = [ DataValueArray(rand(Int, n), fill(true, n)) for i in 1:rand(3:5) ]

    Z2 = map(h2, Ys[1])
    @test isequal(Z2, DataValueArray{Int}(n))
    @test isa(Z2, DataValueArray{Int})
    Z2 = map(h2, Ys[1], Ys[2])
    @test isequal(Z2, DataValueArray{Int}(n))
    @test isa(Z2, DataValueArray{Int})
    Z2 = map(h2, Ys...)
    @test isequal(Z2, DataValueArray{Int}(n))
    @test isa(Z2, DataValueArray{Int})

end


text:
"""
    resnextblock(inplanes, outplanes, cardinality, width, downsample = false)

Create a basic residual block as defined in the paper for ResNeXt
([reference](https://arxiv.org/abs/1611.05431)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: the number of output feature maps
- `cardinality`: the number of groups to use for the convolution
- `width`: the number of feature maps in each group in the bottleneck
- `downsample`: set to `true` to downsample the input
"""
function resnextblock(inplanes, outplanes, cardinality, width, downsample = false)
  stride = downsample ? 2 : 1
  hidden_channels = cardinality * width
  return Chain(conv_bn((1, 1), inplanes, hidden_channels; stride = 1, bias = false),
               conv_bn((3, 3), hidden_channels, hidden_channels;
                        stride = stride, pad = 1, bias = false, groups = cardinality),
               conv_bn((1, 1), hidden_channels, outplanes; stride = 1, bias = false))
end

"""
    resnext(cardinality, width, widen_factor = 2, connection = (x, y) -> @. relu(x) + relu(y);
          block_config, nclasses = 1000)

Create a ResNeXt model
([reference](https://arxiv.org/abs/1611.05431)).

# Arguments
- `cardinality`: the number of groups to use for the convolution
- `width`: the number of feature maps in each group in the bottleneck
- `widen_factor`: the factor by which the width of the bottleneck is increased after each stage
- `connection`: the binary function applied to the output of residual and skip paths in a block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnext(cardinality, width, widen_factor = 2, connection = (x, y) -> @. relu(x) + relu(y);
                 block_config, nclasses = 1000)
  inplanes = 64
  baseplanes = 128
  layers = []
  push!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = (3, 3)))
  push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
  for (i, nrepeats) in enumerate(block_config)
    # output planes within a block
    outplanes = baseplanes * widen_factor
    # push first skip connection on using first residual
    # downsample the residual path if this is the first repetition of a block
    push!(layers, Parallel(connection, resnextblock(inplanes, outplanes, cardinality, width, i != 1),
                                       skip_projection(inplanes, outplanes, i != 1)))
    # push remaining skip connections on using second residual
    inplanes = outplanes
    for _ in 2:nrepeats
        push!(layers, Parallel(connection, resnextblock(inplanes, outplanes, cardinality, width, false),
                                           skip_identity(inplanes, outplanes, false)))
    end
    baseplanes = outplanes
    # double width after every cluster of blocks
    width *= widen_factor
  end

  return Chain(Chain(layers),
               Chain(AdaptiveMeanPool((1, 1)), MLUtils.flatten, Dense(inplanes, nclasses)))
end

"""
    ResNeXt(cardinality, width; block_config, nclasses = 1000)

Create a ResNeXt model
([reference](https://arxiv.org/abs/1611.05431)).

# Arguments
- `cardinality`: the number of groups to use for the convolution
- `width`: the number of feature maps in each group in the bottleneck
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
struct ResNeXt
  layers
end

function ResNeXt(cardinality, width; block_config, nclasses = 1000)
  layers = resnext(cardinality, width; block_config, nclasses)
  ResNeXt(layers)
end

@functor ResNeXt

(m::ResNeXt)(x) = m.layers(x)

backbone(m::ResNeXt) = m.layers[1]
classifier(m::ResNeXt) = m.layers[2]

const resnext_config = Dict(
  50 => (3, 4, 6, 3),
  101 => (3, 4, 23, 3),
  152 => (3, 8, 36, 3)
)

"""
    ResNeXt(config::Int = 50; cardinality = 32, width = 4, pretrain = false, nclasses = 1000)

Create a ResNeXt model with specified configuration. Currently supported values for `config` are (50, 101).
([reference](https://arxiv.org/abs/1611.05431)).
Set `pretrain = true` to load the model with pre-trained weights for ImageNet.

!!! warning
  `ResNeXt` does not currently support pretrained weights.

See also [`Metalhead.resnext`](#).
"""
function ResNeXt(config::Int = 50; cardinality = 32, width = 4, pretrain = false, nclasses = 1000)
  @assert config in keys(resnext_config) "`config` must be one of $(sort(collect(keys(resnext_config))))"

  model = ResNeXt(cardinality, width; block_config = resnext_config[config], nclasses)
  pretrain && loadpretrain!(model, string("ResNeXt", config))
  model
end


text:
struct SpectralParameters{T}
    # plan::Base.DFT.FFTW.cFFTWPlan{Complex{T},-1,true,1}
    # iplan::Base.DFT.ScaledPlan{Complex{T},Base.DFT.FFTW.cFFTWPlan{Complex{T},1,true,1},T}
    plan
    iplan
    A::T
    B::Vector{T}
    force::T
    vv::Vector{Complex{T}}
    T0::T
end

"""
    conv(p,ip,u,v)

Convolution of two vectors. Uses a planned FFT algorithm.

This function directly borrows from the conv in DSP.jl:
https://github.com/JuliaDSP/DSP.jl/blob/master/src/dspbase.jl
"""
function DSP.conv(p,ip,u,v)
    T = eltype(u)
    nu = length(u)
    nv = length(v)
    n = nu+nv-1
    np2 = n>1024 ? nextprod([2,3,5],n) : nextpow(2,n)
    upad = [u;zeros(T,np2-nu)]
    vpad = [v;zeros(T,np2-nv)]

    y = ip*((p*upad).*(p*vpad))
    return y[1:n]
end

"""
    conv_mat(u)
Calculate a matrix `A` such that `A*v = conv(u,v)`.
"""
function conv_mat(u)
    nn = length(u)
    nd2 = round(Int,nn/2)
    inds = (-nd2:nd2-1).+nn.+1

    δ = zeros(eltype(u),nn)
    δ[1] = one(eltype(u))
    v1 = conv(δ,u)[inds][1:end]

    δ *= zero(eltype(u))
    δ[end] = one(eltype(u))
    v2 = conv(δ,u)[inds][end:-1:1]

    v1[1] = v2[1]

    ToeplitzMatrices.Toeplitz(v1,v2)
end

"""
    conv_mat(p,ip,u)
Same as `conv_mat(u)` except it reqires a forward plan `p` and an inverse plan `ip`.
Both plans need to be appropriately sized so that `p*conv(u,u)` will work.
"""
function conv_mat(p,ip,u)
    nn = length(u)
    nd2 = round(Int,nn/2)
    inds = (-nd2:nd2-1).+nn.+1

    δ = zeros(eltype(u),nn)
    δ[1] = one(eltype(u))
    v1 = conv(p,ip,δ,u)[inds][1:end]

    δ *= zero(eltype(u))
    δ[end] = one(eltype(u))
    v2 = conv(p,ip,δ,u)[inds][end:-1:1]

    v1[1] = v2[1]

    ToeplitzMatrices.Toeplitz(v1,v2)
end

function spectral_rhs!(::Type{Val{:jac}},jac,ϕ,params,t)
    nn = length(params.vv)
    p = view(ϕ,1:nn)
    tt = view(ϕ,nn+1:2nn)

    density_pjac!(view(jac,1:nn,1:nn),p,tt,params)
    density_ttjac!(view(jac,1:nn,nn+1:2nn),p,tt,params)
    temperature_pjac!(view(jac,nn+1:2nn,1:nn),tt,p,params)
    temperature_ttjac!(view(jac,nn+1:2nn,nn+1:2nn),tt,p,params)
    jac
end

function density!(dp,p,tt,params)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = -nd2:nd2-1
    dp[:] .= (-im*2π*force*kk.*p
              .-4π^2*( conv(plan,iplan,p.*kk.^2,tt)[kk.+nn.+1]
                      .+conv(plan,iplan,p.*kk,tt.*kk)[kk.+nn.+1])
              .-4π^2*( conv(plan,iplan,p,vv.*kk.^2)[kk.+nn.+1]
                      .+conv(plan,iplan,p.*kk,vv.*kk)[kk.+nn.+1]))
end

function density_pjac!(pjac,p,tt,params)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = (-nd2:nd2-1)*one(eltype(p)) |> collect

    K = Diagonal(kk)
    pjac[:,:] .= -im*2π*force*K.-4π^2*(  conv_mat(plan,iplan,tt)*K.^2 .+conv_mat(plan,iplan,tt.*kk)*K
                                       .+conv_mat(plan,iplan,vv.*kk.^2).+conv_mat(plan,iplan,vv.*kk)*K)
end

function density_ttjac!(ttjac,p,tt,params)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = (-nd2:nd2-1)*one(eltype(p)) |> collect

    K = Diagonal(kk)
    ttjac[:,:] .= -4π^2*(conv_mat(plan,iplan,p.*kk.^2).+conv_mat(plan,iplan,p.*kk)*K)
end

function temperature!(dtt,tt,p,params::SpectralParameters{T}) where T <: Number
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = -nd2:nd2-1
    inds = (-nd2:nd2-1).+nn.+1
    delta_k_0 = zeros(nn)
    delta_k_0[nd2+1] = T0
    dtt[:] .= ( A*force^2*p-2π*im*A*force*conv(plan,iplan,kk.*p,tt)[kk.+nn.+1]
               -4π*im*A*force*conv(plan,iplan,p,kk.*vv)[kk.+nn.+1]
               -4π^2*A*conv(plan,iplan,tt,conv(plan,iplan,kk.*p,kk.*vv)[inds])[inds]
               -4π^2*A*conv(plan,iplan,p,conv(plan,iplan,vv.*kk,vv.*kk)[inds])[inds]
               -4π^2*B[1]*tt.*(kk.^2)
               +B[2]*(delta_k_0-tt))
end

function temperature_pjac!(pjac,tt,p,params)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = (-nd2:nd2-1)*one(eltype(p)) |> collect
    inds = (-nd2:nd2-1).+nn.+1

    K = Diagonal(kk)
    pjac[:,:] .= (A*force^2*Diagonal(ones(nn)).-2π*im*A*force*conv_mat(plan,iplan,tt)*K
                  .-4π*im*A*force*conv_mat(plan,iplan,vv.*kk)
                  .-4π^2*A*(  (conv_mat(plan,iplan,conv(plan,iplan,tt,kk.*vv)[inds]))*K
                            .+conv_mat(plan,iplan,conv(plan,iplan,vv.*kk,vv.*kk)[inds]))
    )
end

function temperature_ttjac!(ttjac,tt,p,params)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = length(axes(p)[1])
    nd2 = round(Int,nn/2)
    kk = (-nd2:nd2-1)*one(eltype(p)) |> collect
    inds = (-nd2:nd2-1).+nn.+1

    K = Diagonal(kk)
    ttjac[:,:] .= (-2π*im*A*force*conv_mat(kk.*p)
                   .-4π^2*A*conv_mat(plan,iplan,conv(plan,iplan,kk.*p,kk.*vv)[inds])
                   .-4π^2*B[1]*K.^2
                   .-B[2]*Diagonal(ones(nn)))
end

function spectral_rhs!(dϕ,ϕ,params,t)
    @unpack plan,iplan,A,B,force,vv,T0 = params
    nn = round(Int,length(ϕ)/2)
    dϕ[1:nn] .= density!(dϕ[1:nn],ϕ[1:nn],ϕ[nn+1:2nn],params)
    dϕ[nn+1:2nn] .= temperature!(dϕ[nn+1:2nn],ϕ[nn+1:2nn],ϕ[1:nn],params)
    dϕ
end

function real_current(ϕ,params)
    nn = round(Int,length(ϕ)/2)
    nd2 = round(Int,nn/2)
    kk = -nd2:nd2-1
    p = ϕ[1:nn]
    tt = ϕ[nn+1:2nn]

    ∇p = -im*2π*kk.*p
    ∇tt = -im*2π*kk.*tt
    ∇vv = -im*2π*kk.*params.vv

    P = fft(fftshift(p))
    T = fft(fftshift(tt))
    ∇P = fft(fftshift(∇p))
    ∇T = fft(fftshift(∇tt))
    ∇V = fft(fftshift(∇vv))

    J = -(P.*(∇V.+params.force).+T.*∇P)
end


text:
function bisection(f::Function, a::Number, b::Number;
                   xtol::AbstractFloat = 1e-1, maxiter::Int = 1000)
    fa = f(a)
    sign(fa) == sign(f(b)) && throw("No real root in [a, b]")
    i = 0
    c = 0.0
    while b-a > xtol
        i += 1
        i == maxiter && throw("Max number of iterations exceeded")
        c  = (a+b)/2
        fc = f(c)
        if fc ≈ 0
            break
        elseif sign(fa) == sign(fc)
            # Root is in the right half of [a, b]
            a = c
            fa = fc
        else
            # Root is in the left half of [a, b]
            b = c
        end
    end
    return c
end

function fast_mvnormal_pdf(x::Vector{Float64})
    coeff_term = (2*pi)^(-length(x)/2)
    exp_term   = exp(-(1/2) * dot(x, x))
    return coeff_term*exp_term
end

function fast_mvnormal_pdf(x::Vector{Float64}, det_Σ::Float64, inv_Σ::Matrix{Float64})
    coeff_term = (2*pi)^(-length(x)/2) * det_Σ^(-1/2)
    exp_term   = exp(-(1/2) * dot(x, inv_Σ * x))
    return coeff_term*exp_term
end


text:
using Test
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil, @is_initialized, @get_package, @get_numbertype, @get_ndims, SUPPORTED_PACKAGES, PKG_CUDA, PKG_NONE, NUMBERTYPE_NONE, NDIMS_NONE
import ParallelStencil: @require, @symbols, longnameof
@static if (PKG_CUDA in SUPPORTED_PACKAGES) import CUDA end

@static for package in SUPPORTED_PACKAGES  eval(:(
    @testset "$(basename(@__FILE__)) (package: $(nameof($package)))" begin
        @testset "1. Reset of ParallelStencil" begin
            @testset "Reset if not initialized" begin
                @require !@is_initialized()
                @reset_parallel_stencil()
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
                @test @get_ndims() == $NDIMS_NONE
            end;
            @testset "Reset if initialized" begin
                @require !@is_initialized()
                @init_parallel_stencil($package, Float64, 3)
                @require @is_initialized() && @get_package() == $package
                @reset_parallel_stencil()
                @test length(@symbols($(@__MODULE__), Data)) == 1
                @test !@is_initialized()
                @test @get_package() == $PKG_NONE
                @test @get_numbertype() == $NUMBERTYPE_NONE
                @test @get_ndims() == $NDIMS_NONE
            end;
        end;
    end;
)) end == nothing || true;


text:
using HorizonSideRobots

function goToSide!(r::Robot, side::HorizonSide; save::Bool=false, mark::Bool = false)::Union{Tuple{Int, HorizonSide}, Nothing}
    if save
        shags = 0
        while !isborder(r, side)
            move!(r, side)
            if mark putmarker!(r) end
            shags+=1
        end
        return (shags, reverseSide(side))
    else
        while !isborder(r, side)
            move!(r, side)
            if mark putmarker!(r) end
        end
    end
end

function goToSideWithShag!(r::Robot, side::HorizonSide, shags::Int; save::Bool = false, checkBorder::Bool = false, mark::Bool = false)::Union{Tuple{Int, HorizonSide}, Nothing}
    if save
        shag = 0
        for _ = 1:shags
            if checkBorder && isborder(r, side) break end
            move!(r, side)
            if mark putmarker!(r) end
            shag+=1
        end
        return (shag, reverseSide(side))
    else
        for _ = 1:shags
            if checkBorder && isborder(r, side) break end
            move!(r, side)
            if mark putmarker!(r) end
        end
    end
end

function aroundBorder!(r::Robot, side::HorizonSide; save::Bool=false)
    listOfShags = []
    k = 0
    while isborder(r, side)
        push!(listOfShags, goToSideWithShag!(r, perSide(side), 1, save = true))
        k+=1
    end
    push!(listOfShags, goToSideWithShag!(r, side, 1, save = true))
    while isborder(r, reverseSide(perSide(side)))
        push!(listOfShags, goToSideWithShag!(r, side, 1, save = true))
    end
    for _ = 1:k
        push!(listOfShags, goToSideWithShag!(r, reverseSide(perSide(side)), 1, save = true))
    end
    if save
        return listOfShags
    end
end

function isRect(r::Robot, side::HorizonSide)::Bool
    if (!isborder(r, side))
        return false
    end
    result = false
    k = 0
    while !isborder(r, perSide(side))
        move!(r, perSide(side))
        k+=1
        if !isborder(r, side)
            result = true
            break
        end
    end
    for _ = 1:k
        move!(r, reverseSide(perSide(side)))
    end
    return result
end

function reverseSide(side::HorizonSide)::HorizonSide
    return HorizonSide(
        (Integer(side)+2) % 4
    )
end

function perSide(side::HorizonSide)::HorizonSide
    return HorizonSide(
        (Integer(side)+1) % 4
    )
end

function goToStart!(r::Robot, listOfShags)::Nothing
    reverse!(listOfShags)
    for tup = listOfShags
        for _ = 1:tup[1]
            move!(r, tup[2])
        end
    end
end

function goToWestSud!(r::Robot)::Nothing
    goToSide!(r, Sud)
    goToSide!(r, West)
end

max_stars_repo_name: Trilonka/MireaAl

max_issues_repo_name: Trilonka/MireaAl

max_forks_repo_name: Trilonka/MireaAl


text:
# import packages
using RoME, Distributions

# [OPTIONAL] add more julia processes to speed up inference
using Distributed
nprocs() < 3 ? addprocs(4-nprocs()) : nothing

@everywhere using RoME

# start with an empty factor graph object
fg = initfg()

# Add the first pose :x0
addVariable!(fg, :x0, Pose2)

# Add at a fixed location PriorPose2 to pin :x0 to a starting location
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3), diagm( [0.1;0.1;0.1].^2 ))))

# Drive around in a hexagon a number of times
for i in 0:12
    psym = Symbol("x$i")
    nsym = Symbol("x$(i+1)")
    addVariable!(fg, nsym, Pose2)
    pp = Pose2Pose2(MvNormal([10.0;0;pi/3], diagm( [0.1;0.1;0.1].^2 ) ))
    addFactor!(fg, [psym;nsym], pp )
end

# perform inference, and remember first runs are slower owing to Julia's just-in-time compiling
solveTree!(fg);

## Inter-operating visualization packages for Caesar/RoME/IncrementalInference exist
using Cairo, RoMEPlotting
Gadfly.set_default_plot_size(35cm,20cm)

# For Juno/Jupyter style use
pl = plotSLAM2D(fg, contour=false)
# For scripting use-cases you can export the image
pl |> Gadfly.PDF("/tmp/test1.pdf", 20cm, 10cm)  # or PNG(...)

# Add a landmark l1
addVariable!(fg, :l1, Point2, tags=[:LANDMARK;])

bear1 = atan(10,20)
bear2 = atan(10,10) - pi/3

# Add landmarks with Bearing range measurements at x0, x6, x12, x18, x24, x30...
# Get variables
vars = ls(fg, r"x") |> sortDFG
for xIndex in 1:6:length(vars)
    @info("Creating factor between $(vars[xIndex]) and l1...")
    p2br = Pose2Point2Bearing(Normal(bear1,0.05))
    addFactor!(fg, [vars[xIndex]; :l1], p2br)
    p2br = Pose2Point2Bearing(Normal(bear2,0.05))
    addFactor!(fg, [vars[xIndex+1]; :l1], p2br)
end

# plotKDE(fg, [:x0;:x6;:x12], dims=[1;2])
# plotKDE(fg, [:x0;:x6;:x12], dims=[3])
# plotKDE(fg, [:x1;:x7;:x13], dims=[1;2])

# drawGraph(fg, show=true)

# Initialize :l1 numerical values but do not rerun solver
ensureAllInitialized!(fg)
pl = plotSLAM2D(fg, contour=false)
pl |> Gadfly.PDF("/tmp/test2.pdf", 20cm, 10cm)  # or PNG(...)

# solve, at time of writing compute time should less than 45s
solveTree!(fg);

# redraw to see that additional information aids in refining the results
pl = plotSLAM2D(fg, contour=false)
pl |> Gadfly.PDF("/tmp/test3.pdf", 20cm, 10cm)  # or PNG(...)


text:
# to compable to most X64 Computers ,you need to RUN this JL with julia -C x86-64
#julia %USERPROFILE%\.julia\packages\PackageCompiler\oT98U\juliac.jl -C haswell -vRe 31-OnlineApp.jl
#julia ~/.julia/packages/PackageCompiler/oT98U/juliac.jl -C haswell -vRe 31-OnlineApp.jl
#julia ~/.julia/packages/PackageCompiler/oT98U/juliac.jl -C ivybridge -vRe 31-OnlineApp.jl
#julia %USERPROFILE%\.julia\packages\PackageCompiler\oT98U\juliac.jl -C ivybridge -vRe 31-OnlineApp.jl

using PackageCompiler

function genbatFileetc( AppDir::String,srcName::String,ExecFolder::String,DLLFolder::String )
    startUpScript = replace( srcName ,Pair(".jl",".bat" ))
    scriptFile = "$AppDir/$startUpScript"

    ExecFolder= replace(ExecFolder,Pair(AppDir,"%cd%" ))
    ExecFolder = replace(ExecFolder,Pair("/","\\") );
    DLLFolder= replace(DLLFolder,Pair(AppDir,"%cd%" ));
    DLLFolder = replace(DLLFolder,Pair("/","\\") );

    f1 =open(scriptFile,"w" )
    write(f1, "set path=%cd%;$ExecFolder;$DLLFolder;%path%\r\n" )
    write(f1,"mkdir $ExecFolder\\Log \r\n")
    write(f1,"pushd $ExecFolder\r\n")

    write(f1,"start cmd /K\r\n")
    close(f1)
end
function CpDlls(tgrDLLFolder::String; DLLS=[] ::String[])
    if Base.Sys.islinux()
        DLLext=".so"
    else
        DLLext=".dll"
    end
    for dllFile in DLLS
        srcDllName = joinpath( ".",dllFile*DLLext)
        tgrDllName = joinpath( tgrDLLFolder ,dllFile*DLLext )
        Base.Filesystem.cp(srcDllName,tgrDllName,force=true)
    end
end

function BuildExecute(;AppDir::String,appname="exec" ,srcName="" )
    glbCfgFolder="glbCfg"
    buildFolder   = "$AppDir/$appname/"
    DLLFolder  = "$AppDir/DLL/"
    ExecFolder = buildFolder
    tgrcfgFolder = "$ExecFolder/$glbCfgFolder/"
    scfCfgFolder = "./$glbCfgFolder/"

    mkpath( tgrcfgFolder)
    mkpath( ExecFolder)
    mkpath( DLLFolder)

    Base.Filesystem.cp(scfCfgFolder,tgrcfgFolder,force=true)
    genbatFileetc( AppDir,srcName,ExecFolder,DLLFolder )

    build_executable(
        srcName,
        # snoopfile = "call_functions.jl", # Julia script which calls functions that you want to make sure to have precompiled [optional]
        builddir = buildFolder
        )
end
##################################################
AppDir="../Tem/"

tgrDLLFolder  = "$AppDir/DLL/"

mkpath( tgrDLLFolder)
CpDlls(tgrDLLFolder,DLLS=["RawProc","RawDataSave","partSnap"])

#BuildExecute(;AppDir=AppDir,appname="OfflineApp" ,srcName="03-OfflineApp.jl" )
#exit(-1)
BuildExecute(;AppDir=AppDir,appname="OnlineApp" ,srcName="31-OnlineApp.jl" )
exit(-1)

# BuildExecute(;AppDir=AppDir,appname="hw2PTFileTransfer" ,srcName="01-Apphw2PTFileTransfer.jl" )
# exit(-1)

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
######
##FEM Heat Method Tests
######
using FiniteElementDiffEq, DiffEqDevTools

#Define a parabolic problem
prob = prob_femheat_moving7

#Solve it with a bunch of different algorithms, plot solution
println("Euler")
sol = solve(prob,FEMDiffEqHeatEuler())
T=1
dt = 1//2^(4) #Make faster for tests
fem_mesh = parabolic_squaremesh([0 1 0 1],dx,dt,T,:dirichlet)
println("Direct")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),solver=:Direct)

println("LU")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),solver=:LU)

println("QR")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),solver=:QR)

println("SVD")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),solver=:SVD)

println("Direct")
sol = solve(prob,FEMDiffEqHeatCrankNicholson(),solver=:Direct)

println("Cholesky")
sol = solve(prob,FEMDiffEqHeatCrankNicholson(),solver=:Cholesky)

println("CG")
sol = solve(prob,FEMDiffEqHeatCrankNicholson(),solver=:CG)

println("CG")
sol = solve(prob,FEMDiffEqHeatImplicitEuler(),solver=:CG)

println("GMRES")
sol = solve(prob,FEMDiffEqHeatCrankNicholson(),solver=:GMRES)

#Define another parabolic problem
prob = prob_femheat_diffuse #also try heatProblemExample_pure() or heatProblemExample_diffuse()

#Solve it with a bunch of different algorithms, plot solution
println("Euler")
sol = solve(prob,FEMDiffEqHeatEuler())

#Define a different parabolic problem
prob = prob_femheat_pure

#Solve with Euler Method
println("Euler")
sol = solve(prob,FEMDiffEqHeatEuler())

#Choose a finer mesh, solve with Euler, and add this result to the previous as
#an approximately true solution.
prob = prob_femheat_pure11
sol2 = solve(prob,FEMDiffEqHeatEuler())
sol3 = appxtrue(sol,sol2)
#plot(sol,plot_analytic=true,cbar=false)

sol3.errors[:l2]<.005 #Returns true if res solution is near the apprxTrue res2 solution


text:
argmaxrow(A::AbstractArray{T}) where {T<:Real} = [argmax(A[i, :]) for i = 1:size(A, 1)]

function lagou(a::AbstractArray; dims = 1, kwargs...)
    T = size(a, dims)
    mapslices(x -> lagou(x, T; kwargs...), a, dims = dims)
end

function lagou(hmm::PeriodicHMM; kwargs...)
    T = size(hmm.A, 3)
    lagou(hmm, T; kwargs...)
end

#! ERROR: LoadError: LoadError: TypeError: in MixtureModel, in CT, expected CT<:(Distributions.Categorical{P, Ps} where {P<:Real, Ps<:AbstractVector{P}}), got Type{Float64}
# function lagou(mix::Array{MixtureModel{Univariate, Continuous, Distribution{Univariate, Continuous}, Float64}, 3}; kwargs...)
#     T = size(mix, 2)
#     lagou(mix, T; kwargs...)
# end

function lagou(hmm::HierarchicalPeriodicHMM; kwargs...)
    T = size(hmm.A, 3)
    lagou(hmm, T; kwargs...)
end

function lagou(a::AbstractVector, T::Int; window = -15:15, kernel = :step)
    cycle = CyclicArray(1:T, "1D")
    if collect(window) == [0]
        return a
    elseif kernel == :step
        return [mean(a[cycle[t+tt]] for tt in window) for t = 1:T]
    elseif kernel == :step_fast # faster but produces for small positive values small negative numbers
        m = similar(a)
        N = length(window)
        m[1] = mean(a[tt] for tt in cycle[(window).+1])
        for t in eachindex(a)[2:end]
            m[t] = (N * m[t-1] + a[cycle[t+window.stop]] - a[cycle[t-1+window.start]]) / N
        end
        return m
    elseif kernel == :Epanechnikov
        ker(u) = 3 / 4 * (1 - abs(u)^2)
        s = [sum(ker(tt / window.stop) for tt in window) for t = 1:T]
        return [sum(a[cycle[t+tt]] * ker(tt / window.stop) for tt in window) for t = 1:T] ./ s
    end
end

function lagou(B::Vector{F}, T::Int; kwargs...) where {F<:Bernoulli}
    K = size(B, 1)
    return Bernoulli.(lagou(succprob.(B), T; kwargs...))
end

function lagou(B::Vector{Product{Discrete,F,Vector{F}}}, T::Int; kwargs...) where {F<:Categorical}
    K = size(B, 1)
    D = length(B[1])
    smooth_param = lagou([probs.(B[t].v) for t = 1:T], T; kwargs...)
    d = [product_distribution([Categorical(abs.(smooth_param[t][j])) for j = 1:D]) for t = 1:T]
    return d
end

function lagou(B::Vector{Product{Discrete,F,Vector{F}}}, T::Int; kwargs...) where {F<:Bernoulli}
    K = size(B, 1)
    D = length(B[1])
    smooth_param = lagou([succprob.(B[t].v) for t = 1:T], T; kwargs...)
    d = [product_distribution([Bernoulli(abs.(smooth_param[t][j])) for j = 1:D]) for t = 1:T]
    return d
end

function lagou(hmm::PeriodicHMM{Multivariate}, T::Int; kwargs...)
    K = size(hmm.A, 1)
    # abs to ensure positivness (sometime means of super small nuber produce negative number small in abs value)
    Q = lagou(hmm.A, dims = 3; kwargs...)

    d = lagou(hmm.B, dims = 2; kwargs...)

    return PeriodicHMM(Q, d)
end

function lagou(B::AbstractVector{F} where {F<:Bernoulli}, T; kwargs...)
    param = succprob.(B)
    lparam = lagou(param, T; kwargs...)
    Bernoulli.(abs.(lparam[t]))
end

function lagou(hmm::HierarchicalPeriodicHMM, T::Int; kwargs...)
    cycle = CyclicArray(1:T, "1D")
    K = size(hmm.A, 1)

    Q = lagou(hmm.A, dims = 3; kwargs...)

    d = lagou(hmm.B, dims = 2; kwargs...)

    return HierarchicalPeriodicHMM(Q, d)
end

function lagou(mix::Vector{F}, T::Int; kwargs...) where {F<:MixtureModel{Univariate}}
    data_type = typeof.(Distributions.components(mix[1]))
    C = length(data_type)
    mix_s = similar(mix)
    p = params.(mix)
    smooth_α = lagou(getindex.(p, 2); kwargs...)
    [p[t][2][:] = smooth_α[t] for t = 1:T]
    for c = 1:C
        len = length(p[1][1][c])
        smooth_c = zeros(T, len)
        for l = 1:len
            smooth_c[:, l] = lagou([p[t][1][c][l] for t = 1:T]; kwargs...)
        end
        [p[t][1][c] = tuple([smooth_c[t, l] for l = 1:len]...) for t = 1:T]
    end
    for t = 1:T
        mix_s[t] = MixtureModel([data_type[c](p[t][1][c]...) for c = 1:C], p[t][2])
    end
    return mix_s
end

# function lagou(mix::Array{MixtureModel{Univariate, Continuous, Distribution{Univariate, Continuous}, Float64}, 3}, T::Int; kwargs...)
#     data_type = typeof.(Distributions.components(mix[1]))
#     C = length(data_type)
#     K, D = size(mix, 1), size(mix, 3)
#     mix_s = similar(mix)
#     for k in 1:K, j in 1:D
#         p = params.(mix[k,:,j])
#         ps = lagou(p; kwargs...)
#         for t in 1:T
#             mix_s[k,t,j] = MixtureModel([data_type[c](ps[t][1][c]...) for c in 1:C], ps[t][2])
#         end
#     end
#     return mix_s
# end
#
# index_cum(c, s) = c == 1 ? 1 : s[c-1]+1
#
# function lagou(p::Vector{Tuple{Vector{Tuple{Float64, Vararg{Float64, N} where N}}, Vector{Float64}}}, T::Int; kwargs...)
#     component = length(p[1][1])
#     size_component = [length(p[1][1][c]) for c in 1:component]
#     ca = unfold.(p)
#
#     cas = lagou(ca, T, kwargs...)
#     [tuple([Tuple(cas[t][index_cum(c, size_component):cumsum(size_component)[c]]) for c in 1:component], cas[t][(sum(size_component)+1):end]) for t in 1:T]
# end
##

# x[first] = previous day -> n-1
# x[last] = farest day -> n-> c
bin2digit(x) = sum(x[length(x)-i+1] * 2^(i - 1) for i = 1:length(x)) + 1
bin2digit(x::Tuple) = bin2digit([x...])

function dayx(lag_obs)
    memory = length(lag_obs)
    t = tuple.([lag_obs[m] for m = 1:memory]...)
    bin2digit.(t)
end

function conditional_to(observations::AbstractArray, memory::Int;
    past = [0 1 0 1 1 0 1 0 0 0
        1 1 0 1 1 1 1 1 1 1
        1 1 0 1 1 1 0 1 1 1
        1 1 0 1 1 0 0 0 1 0
        1 1 0 1 1 0 0 1 0 1]
)
    if memory == 0
        return ones(Int, size(observations))
    elseif memory < 5
        lag_obs = [copy(lag(observations, m)) for m = 1:memory]  # transform dry=0 into 1 and wet=1 into 2 for array indexing
        [lag_obs[m][1:m, :] .= reverse(past[1:m, :], dims = 1) for m = 1:memory] # avoid the missing first row

        return dayx(lag_obs)
    else
        println("Not implemented need to add initial conditions (or fake ones)")
    end
end

function conditional_to(observations::AbstractArray, memory::AbstractVector;
    past = [0 1 0 1 1 0 1 0 0 0
        1 1 0 1 1 1 1 1 1 1
        1 1 0 1 1 1 0 1 1 1
        1 1 0 1 1 0 0 0 1 0
        1 1 0 1 1 0 0 1 0 1]
)
    D = size(observations, 2)
    lag_cat = zeros(Int, size(observations))
    for j = 1:D
        lag_cat[:, j] = conditional_to(observations[:, j], memory[j], past = past[:, j])
    end
    return lag_cat
end

# function conditional_to(observations::AbstractArray, memory::Float64; past =  [0  1  0  1  1  0  1  0  0  0;
#                                                                 1  1  0  1  1  1  1  1  1  1;
#                                                                 1  1  0  1  1  1  0  1  1  1;
#                                                                 1  1  0  1  1  0  0  0  1  0;
#                                                                 1  1  0  1  1  0  0  1  0  1]
#                         )
#     if memory == 0
#         return ones(Int, size(observations))
#     elseif memory<5
#         lag_obs = [copy(lag(observations, m)) for m in 1:memory]  # transform dry=0 into 1 and wet=1 into 2 for array indexing
#         [lag_obs[m][1:m,:] .= reverse(past[1:m,:], dims=1) for m in 1:memory] # avoid the missing first row

#         return dayx(lag_obs)
#     else
#         println("Not implemented need to add initial conditions (or fake ones)")
#     end
# end

function idx_observation_of_past_cat(lag_cat, n2t, T, K, size_memory)
    # Matrix(T,D) of vector that give the index of data of same past.
    # ie. size_memory = 1 (no memory) -> every data is in category 1
    # ie size_memory = 2 (memory on previous day) -> idx_tj[t,j][1] = vector of index of data where previous day was dry, idx_tj[t,j][2] = index of data where previous day was wet
    D = size(lag_cat, 2)
    idx_tj = Matrix{Vector{Vector{Int}}}(undef, T, D)
    n_in_t = [findall(n2t .== t) for t = 1:T] # could probably be speeded up e.g. recusivly suppressing already classified label with like sortperm
    @inbounds for t in OneTo(T)
        n_t = n_in_t[t]
        for i in OneTo(K)
            for j = 1:D
                # apparently the two following takes quite long ~29ms for all the loops
                # TODO improve time!
                n_tm = [findall(lag_cat[n_t, j] .== m) for m = 1:size_memory]
                idx_tj[t, j] = [n_t[n_tm[m]] for m = 1:size_memory]
                ##
            end
        end
    end
    return idx_tj
end

## Trigo part

function polynomial_trigo(t::Number, β; T = 366)
    d = (length(β) - 1) ÷ 2
    if d == 0
        return β[1]
    else
        f = 2π / T
        # everything is shifted from 1 from usual notation due to array starting at 1
        return β[1] + sum(β[2*l] * cos(f * l * t) + β[2*l+1] * sin(f * l * t) for l = 1:d)
    end
end

function polynomial_trigo(t::AbstractArray, β; T = 366)
    d = (length(β) - 1) ÷ 2
    if d == 0
        return β[1]
    else
        f = 2π / T
        # everything is shifted from 1 from usual notation due to array starting at 1
        return β[1] .+ sum(β[2*l] * cos.(f * l * t) + β[2*l+1] * sin.(f * l * t) for l = 1:d)
    end
end

interleave2(args...) = collect(Iterators.flatten(zip(args...))) # merge two vector with alternate elements

function fit_Q!(p::AbstractArray, A::AbstractArray{N,2} where {N}; silence = true)
    T, K = size(A, 2), size(A, 1)
    @assert K - 1 == size(p, 1)
    d = (size(p, 2) - 1) ÷ 2
    model = Model(Ipopt.Optimizer)
    silence && set_silent(model)
    f = 2π / T
    cos_nj = [cos(f * j * t) for t = 1:T, j = 1:d]
    sin_nj = [sin(f * j * t) for t = 1:T, j = 1:d]

    trig = [[1; interleave2(cos_nj[t, :], sin_nj[t, :])] for t = 1:T]

    @variable(model, p_jump[k = 1:(K-1), j = 1:(2d+1)])
    set_start_value.(p_jump, p)
    # Polynomial P_kl

    @NLexpression(model, Pol[t = 1:T, k = 1:K-1], sum(trig[t][j] * p_jump[k, j] for j = 1:length(trig[t])))

    @NLobjective(
        model,
        Min,
        sum((A[k, t] - exp(Pol[t, k]) / (1 + sum(exp(Pol[t, l]) for l = 1:K-1)))^2 for k = 1:K-1, t = 1:T)
        +
        sum((A[K, t] - 1 / (1 + sum(exp(Pol[t, l]) for l = 1:K-1)))^2 for t = 1:T)
    )
    optimize!(model)
    p[:, :] = value.(p_jump)
end

m_Bernoulli(t, p; T = 366) = 1 ./ (1 .+ exp.(polynomial_trigo(t, p; T = T)))
function fit_Y!(p::AbstractVector, B::AbstractVector)
    T = size(B, 1)
    p[:] = curve_fit((t, p) -> m_Bernoulli(t, p, T = T), collect(1:T), B, convert(Vector, p)).param
end

function complete_loglikelihood(hmm::HierarchicalPeriodicHMM, y::AbstractArray, z::AbstractVector, n2t::AbstractArray)
    N, size_memory, D = size(y, 1), size(hmm, 4), size(y, 2)
    memory = Int(log(size_memory) / log(2))
    lag_cat = conditional_to(y, memory)

    return sum(log(hmm.A[z[n], z[n+1], n2t[n]]) for n = 1:N-1) + sum(logpdf(product_distribution(hmm.B[CartesianIndex.(z[n], n2t[n], 1:D, lag_cat[n, :])]), y[n, :]) for n = 1:N)
end

function complete_loglikelihood(hmm::PeriodicHMM, y::AbstractArray, z::AbstractVector, n2t::AbstractArray)
    N, D = size(y, 1), size(y, 2)

    return sum(log(hmm.A[z[n], z[n+1], n2t[n]]) for n = 1:N-1) + sum(logpdf(product_distribution(hmm.B[CartesianIndex.(z[n], n2t[n], 1:D)]), y[n, :]) for n = 1:N)
end

function complete_loglikelihood(hmm::HMM, y::AbstractArray, z::AbstractVector, n2t::AbstractArray)
    N, D = size(y, 1), size(y, 2)

    return sum(log(hmm.A[z[n], z[n+1]]) for n = 1:N-1) + sum(logpdf(product_distribution(hmm.B[CartesianIndex.(z[n], 1:D)]), y[n, :]) for n = 1:N)
end

# function pml_loglikelihood(hmm::HierarchicalPeriodicHMM, y::AbstractArray, n2t::AbstractArray)
#     N, size_memory, D = size(y, 1), size(hmm, 4), size(y, 2)
#     memory = Int(log(size_memory)/log(2))
#     lag_cat = conditional_to(y, memory)

#     # return sum( sum(logpdf(product_distribution(hmm.B[CartesianIndex.(i, n2t[n], 1:D, lag_cat[n,:])]), y[n,:]) for n in 1:N)
# end

# slower
# function fit_Y2!(p::AbstractVector, B::AbstractVector; silence = true)
#     T = size(B,1)
#     d = (size(p,1)-1)÷2
#     model = Model(Ipopt.Optimizer)
#     silence && set_silent(model)
#     f = 2π/T
#     cos_nj = [cos(f*j*t) for t in 1:T, j in 1:d]
#     sin_nj = [sin(f*j*t) for t in 1:T, j in 1:d]

#     trig = [[1; interleave2(cos_nj[t,:], sin_nj[t,:])] for t in 1:T]

#     @variable(model, p_jump[j = 1:(2d+1)])
#     set_start_value.(p_jump, p)
#     # Polynomial P_kl

#     @NLexpression(model, Pol[t = 1:T], sum(trig[t][j] * p_jump[j] for j in 1:length(trig[t])))

#     @NLobjective(
#     model,
#     Min,
#         sum( (B[t] - 1/(1+exp(Pol[t])))^2 for t in 1:T)
#     )
#     optimize!(model)
#     p[:] = value.(p_jump)
# end

# function fit_Y2(p::AbstractVector, B::AbstractVector; silence = true)
#     T = size(B,1)
#     d = (size(p,1)-1)÷2
#     model = Model(Ipopt.Optimizer)
#     silence && set_silent(model)
#     f = 2π/T
#     cos_nj = [cos(f*j*t) for t in 1:T, j in 1:d]
#     sin_nj = [sin(f*j*t) for t in 1:T, j in 1:d]

#     trig = [[1; interleave2(cos_nj[t,:], sin_nj[t,:])] for t in 1:T]

#     @variable(model, p_jump[j = 1:(2d+1)])
#     set_start_value.(p_jump, p)
#     # Polynomial P_kl

#     @NLexpression(model, Pol[t = 1:T], sum(trig[t][j] * p_jump[j] for j in 1:length(trig[t])))

#     @NLobjective(
#     model,
#     Min,
#         sum( (B[t] - 1/(1+exp(Pol[t])))^2 for t in 1:T)
#     )
#     return model
# end


text:
using AtomsBase
using Unitful
using Test

@testset "Printing atomic systems" begin
    at = Atom(:Si, zeros(3) * u"m", extradata=42)
    @test repr(at) == "Atom(Si, [0.0, 0.0, 0.0]u\"m\")"

    atoms = [:Si => [0.0, -0.125, 0.0],
             :C  => [0.125, 0.0, 0.0]]
    box = [[10, 0.0, 0.0], [0.0, 5, 0.0], [0.0, 0.0, 7]]u"Å"

    flexible_system = periodic_system(atoms, box; fractional=true)
    @test repr(flexible_system) == """
    FlexibleSystem(CSi, Periodic, box=[[10.0, 0.0, 0.0], [0.0, 5.0, 0.0], [0.0, 0.0, 7.0]]u"Å")"""

    fast_system = FastSystem(flexible_system)
    @test repr(fast_system) == """
    FastSystem(CSi, Periodic, box=[[10.0, 0.0, 0.0], [0.0, 5.0, 0.0], [0.0, 0.0, 7.0]]u"Å")"""
end