text:
# This file is part of IntegerSequences.
# Copyright Peter Luschny. License is MIT.

(@__DIR__) ∉ LOAD_PATH && push!(LOAD_PATH, (@__DIR__))

module FigurativeNumbers

export ModuleFigurativeNumbers
export PolygonalNumber, PyramidalNumber
export V014107, V095794, V067998, V080956, V001477, V000217, V000290
export V000326, V000384, V000566, V000567, V001106, V001107
export V005564, V058373, V254749, V000292, V000330, V002411, V002412
export V002413, V002414, V007584, V007585

"""

* PolygonalNumber, PyramidalNumber, V014107, V095794, V067998, V080956, V001477, V000217, V000290, V000326, V000384, V000566, V000567, V001106, V001107, V005564, V058373, V254749, V000292, V000330, V002411, V002412, V002413, V002414, V007584, V007585
"""
const ModuleFigurativeNumbers = ""

"""

Return the polygonal number with shape k.
"""
function PolygonalNumber(n, k)
    s = div(n^2 * (k - 2) - n * (k - 4), 2)
    k < 2 ? -s : s
end

"""

Return the pyramidal number with shape k.
"""
function PyramidalNumber(n, k)
    s = div(3 * n^2 + n^3 * (k - 2) - n * (k - 5), 6)
    k < 2 ? -s : s
end

"""

Return the polygonal numbers of shape -2.
"""
V014107(n) = PolygonalNumber(n, -2)
"""

Return the polygonal numbers of shape -1.
"""
V095794(n) = PolygonalNumber(n, -1)
"""

Return the polygonal numbers of shape 0.
"""
V067998(n) = PolygonalNumber(n, 0)
"""

Return the polygonal numbers of shape 1.
"""
V080956(n) = PolygonalNumber(n, 1)
"""

Return the polygonal numbers of shape 2 (these are the natural numbers).
"""
V001477(n) = PolygonalNumber(n, 2)
"""

Return the polygonal numbers of shape 3 (the triangular numbers).
"""
V000217(n) = PolygonalNumber(n, 3)
"""

Return the polygonal numbers of shape 4 (the squares).
"""
V000290(n) = PolygonalNumber(n, 4)
"""

Return the polygonal numbers of shape 5 (the pentagonal numbers).
"""
V000326(n) = PolygonalNumber(n, 5)
"""

Return the polygonal numbers of shape 6 (the hexagonal numbers).
"""
V000384(n) = PolygonalNumber(n, 6)
"""

Return the polygonal numbers of shape 7 (the heptagonal numbers).
"""
V000566(n) = PolygonalNumber(n, 7)
"""

Return the polygonal numbers of shape 8 (the octagonal numbers).
"""
V000567(n) = PolygonalNumber(n, 8)
"""

Return the polygonal numbers of shape 9 (the nonagonal numbers).
"""
V001106(n) = PolygonalNumber(n, 9)
"""

Return the polygonal numbers of shape 10 (decagonal numbers).
"""
V001107(n) = PolygonalNumber(n, 10)

"""

Return the pyramidal numbers of shape -1.
"""
V005564(n) = PyramidalNumber(n, -1)
"""

Return the pyramidal numbers of shape 0.
"""
V058373(n) = PyramidalNumber(n, 0)
"""

Return the pyramidal numbers of shape 1.
"""
V254749(n) = PyramidalNumber(n, 1)
#"""
#Return the pyramidal numbers of shape 2 (triangular numbers).
#"""
#V000217(n) = PyramidalNumber(n, 2)
"""

Return the pyramidal numbers of shape 3 (tetrahedral numbers).
"""
V000292(n) = PyramidalNumber(n, 3)
"""

Return the pyramidal numbers of shape 4 (square pyramidal numbers).
"""
V000330(n) = PyramidalNumber(n, 4)
"""

Return the pyramidal numbers of shape 5 (pentagonal pyramidal numbers).
"""
V002411(n) = PyramidalNumber(n, 5)
"""

Return the pyramidal numbers of shape 6 (hexagonal pyramidal numbers).
"""
V002412(n) = PyramidalNumber(n, 6)
"""

Return the pyramidal numbers of shape 7 (heptagonal pyramidal numbers).
"""
V002413(n) = PyramidalNumber(n, 7)
"""

Return the pyramidal numbers of shape 8 (octagonal pyramidal numbers).
"""
V002414(n) = PyramidalNumber(n, 8)
"""

Return the pyramidal numbers of shape 9 (enneagonal pyramidal numbers).
"""
V007584(n) = PyramidalNumber(n, 9)
"""

Return the pyramidal numbers of shape 10 (decagonal pyramidal numbers).
"""
V007585(n) = PyramidalNumber(n, 10)

#START-TEST-########################################################

using Test, SeqTests

function test()

    @testset "Figurative" begin
        @test V002411(1000) == 500500000

        if data_installed()

            V = [
                V014107,
                V067998,
                V001477,
                V000217,
                V000290,
                V000326,
                V000384,
                V000566,
                V000567,
                V001106,
                V001107,
                V000292,
                V000330,
                V002411,
                V002412,
                V002413,
                V007584,
                V007585
            ]
                # V080956, V095794, V005564, V058373, V254749, V002414
            for v ∈ V
                SeqTest(v, 'V')
            end
        end
    end
end

function demo()
    for k ∈ -2:10
        V = [PolygonalNumber(n, k) for n ∈ 0:9]
        println(k, " ", V)
    end

    for k ∈ -2:10
        V = [PyramidalNumber(n, k) for n ∈ 0:9]
        println(k, " ", V)
    end
end

"""

"""
function perf()
    @time (for n ∈ 1:1000
        V000326(n)
    end)
end

function main()
    test()
    demo()
    perf()
end

main()

#-2 [0, -1,  2, 9, 20, 35,  54,  77, 104, 135]
#-1 [0, -1,  1, 6, 14, 25,  39,  56,  76,  99]
# 0 [0, -1,  0, 3,  8, 15,  24,  35,  48,  63]
# 1 [0, -1, -1, 0,  2,  5,   9,  14,  20,  27]
# --------------------------------------------
# 2 [0, 1,  2,  3,  4,  5,   6,   7,   8,   9]
# 3 [0, 1,  3,  6, 10, 15,  21,  28,  36,  45]
# 4 [0, 1,  4,  9, 16, 25,  36,  49,  64,  81]
# 5 [0, 1,  5, 12, 22, 35,  51,  70,  92, 117]
# 6 [0, 1,  6, 15, 28, 45,  66,  91, 120, 153]
# 7 [0, 1,  7, 18, 34, 55,  81, 112, 148, 189]
# 8 [0, 1,  8, 21, 40, 65,  96, 133, 176, 225]
# 9 [0, 1,  9, 24, 46, 75, 111, 154, 204, 261]
#10 [0, 1, 10, 27, 52, 85, 126, 175, 232, 297]

# ==============================================
# -2 [0,-1,  1, 10, 30,  65, 119, 196, 300,  435]
# -1 [0,-1,  0,  6, 20,  45,  84, 140, 216,  315] V005564
#  0 [0,-1, -1,  2, 10,  25,  49,  84, 132,  195] V058373
#  1 [0,-1, -2, -2,  0,   5,  14,  28,  48,   75] V254749
# -----------------------------------------------
#  2 [0, 1,  3,  6, 10,  15,  21,  28,  36,   45] V000217
#  3 [0, 1,  4, 10, 20,  35,  56,  84, 120,  165] V000292
#  4 [0, 1,  5, 14, 30,  55,  91, 140, 204,  285] V000330
#  5 [0, 1,  6, 18, 40,  75, 126, 196, 288,  405] V002411
#  6 [0, 1,  7, 22, 50,  95, 161, 252, 372,  525] V002412
#  7 [0, 1,  8, 26, 60, 115, 196, 308, 456,  645] V002413
#  8 [0, 1,  9, 30, 70, 135, 231, 364, 540,  765] V002414
#  9 [0, 1, 10, 34, 80, 155, 266, 420, 624,  885] V007584
# 10 [0, 1, 11, 38, 90, 175, 301, 476, 708, 1005] V007585

end # module

max_stars_repo_name: OpenLibMathSeq/Sequences

max_issues_repo_name: OpenLibMathSeq/Sequences


text:
# stolen from https://github.com/amitmurthy/AWS.jl/blob/master/src/crypto.jl
# "since OpenSSL.jl does not yet support the functions this package
# needs. This will be replaced with calls to OpenSSL.jl when it
# supports the same"

module Crypto

macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, eval(lib))) ), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

@ctypedef EVP_MD Void
@ctypedef EVP_MD_CTX Void
@ctypedef ENGINE Void

typealias size_t Csize_t

@unix_only const libcrypto = "libcrypto"
@windows_only const libcrypto = Pkg.dir("WinRPM","deps","usr","$(Sys.ARCH)-w64-mingw32","sys-root","mingw","bin","libcrypto-10")

@c Ptr{UInt8} HMAC (Ptr{EVP_MD}, Ptr{Void}, Int32, Ptr{UInt8}, size_t, Ptr{UInt8}, Ptr{UInt32}) libcrypto
@c Ptr{UInt8} MD5 (Ptr{UInt8}, size_t, Ptr{UInt8}) libcrypto
@c Ptr{UInt8} SHA256 (Ptr{UInt8}, size_t, Ptr{UInt8}) libcrypto

@c Ptr{EVP_MD} EVP_md5 () libcrypto
@c Ptr{EVP_MD} EVP_sha1 () libcrypto
@c Ptr{EVP_MD} EVP_sha256 () libcrypto

@c Ptr{EVP_MD_CTX} EVP_MD_CTX_create () libcrypto
@c Int32 EVP_DigestInit_ex (Ptr{EVP_MD_CTX}, Ptr{EVP_MD}, Ptr{ENGINE}) libcrypto
@c Int32 EVP_DigestUpdate (Ptr{EVP_MD_CTX}, Ptr{Void}, size_t) libcrypto
@c Int32 EVP_DigestFinal_ex (Ptr{EVP_MD_CTX}, Ptr{UInt8}, Ptr{UInt32}) libcrypto
@c Union{} EVP_MD_CTX_destroy (Ptr{EVP_MD_CTX},) libcrypto

hmacsha256_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}) =  hmacsha_digest(s, k, EVP_sha256(), 32)
export hmacsha256_digest

hmacsha1_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}) = hmacsha_digest(s, k, EVP_sha1(), 20)
export hmacsha1_digest

function hmacsha_digest(s::AbstractString, k::Union{ASCIIString, Vector{UInt8}}, evp, dgst_len)
    if evp == C_NULL error("EVP_sha1() failed!") end

    sig = zeros(UInt8, dgst_len)
    sig_len = zeros(UInt32, 1)

    if isa(k, AbstractString)
        k = convert(Array{UInt8}, k)
    end

    if HMAC(evp, k, length(k), s, length(s), sig, sig_len) == C_NULL error("HMAC() failed!") end
    if (sig_len[1] != dgst_len) error("Wrong length of signature!") end

    return sig
end

function sha256(s::AbstractString)
    sha = zeros(UInt8, 32)
    assert(SHA256(s, length(s), sha) != C_NULL)
    return sha
end

end # Module end

ext: jl


text:
using PatternFormation
using Test

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


text:
type BatchNorm1{T} <: Layer{T}
    # https://arxiv.org/pdf/1502.03167.pdf
    name::Symbol
    γ::Vector{T}
    β::Vector{T}
    ϵ::T
    μ::Vector{T} # train set mean
    ν::Vector{T} # train set variance

    input::Vector{T}
    output::Vector{T}
end
name(a::BatchNorm1) = a.name
output(a::BatchNorm1) = a.output
function forward!(a::BatchNorm1)
    x = (a.input - a.μ) ./ sqrt(a.ν + a.ϵ)
    a.output .= (a.γ .* x) .+ a.β
    a
end
function Base.push!{T}(net::ForwardNet{T}, ::Type{BatchNorm1},
    name::Symbol,
    parent::NameOrIndex,
    ϵ::T=T(1e-5)
    )

    input = output(net[parent])::Vector{T}
    output = Array{T}(length(input))
    γ = Array{T}(length(input))
    β = Array{T}(length(input))
    μ = Array{T}(length(input))
    ν = Array{T}(length(input))

    node = BatchNorm1(name, γ, β, ϵ, μ, ν, input, output)
    push!(net, node, parent)
end
Base.push!{T}(net::ForwardNet{T}, ::Type{BatchNorm1}, parent::NameOrIndex, ϵ::T=T(1e-5)) =
    push!(net, BatchNorm1, create_next_name!(net, "BatchNorm1"), parent, ϵ)
function restore!{T}(a::BatchNorm1{T},
    filename_gamma::String,
    filename_beta::String,
    filename_mean::String,
    filename_variance::String)

    copy!(a.γ, open(io->read_binary_vec(io, T), filename_gamma))
    copy!(a.β, open(io->read_binary_vec(io, T), filename_beta))
    copy!(a.μ, open(io->read_binary_vec(io, T), filename_mean))
    copy!(a.ν, open(io->read_binary_vec(io, T), filename_variance))

    a
end

type BatchNorm3{T} <: Layer{T}
    # https://arxiv.org/pdf/1502.03167.pdf
    name::Symbol
    γ::Vector{T}
    β::Vector{T}
    ϵ::T
    μ::Vector{T} # train set mean
    ν::Vector{T} # train set variance

    input::Array{T, 3}
    output::Array{T, 3}
end
name(a::BatchNorm3) = a.name
output(a::BatchNorm3) = a.output
function forward!(a::BatchNorm3)
    for i in 1 : size(a.output)[3]
        x = (a.input[:, :, i] - a.μ[i]) / sqrt(a.ν[i] + a.ϵ)
        a.output[:, :, i] = a.γ[i] * x + a.β[i]
    end
    a
end
function Base.push!{T}(net::ForwardNet{T}, ::Type{BatchNorm3},
    name::Symbol,
    parent::NameOrIndex,
    ϵ::T=T(1e-5)
    )

    input = output(net[parent_index])::Array{T, 3}
    output = Array{T}(size(input))
    in_c = size(input, 3)
    γ = Array{T}(in_c)
    β = Array{T}(in_c)
    μ = Array{T}(in_c)
    ν = Array{T}(in_c)

    node = BatchNorm3(name, γ, β, ϵ, μ, ν, input, output)
    push!(net, node, parent)
end
Base.push!{T}(net::ForwardNet{T}, ::Type{BatchNorm3}, parent::NameOrIndex, ϵ::T=T(1e-5)) =
    push!(net, BatchNorm3, create_next_name!(net, "BatchNorm3"), parent, ϵ)
function restore!{T}(a::BatchNorm3{T},
    filename_gamma::String,
    filename_beta::String,
    filename_mean::String,
    filename_variance::String)

    copy!(a.γ, open(io->read_binary_vec(io, T), filename_gamma))
    copy!(a.β, open(io->read_binary_vec(io, T), filename_beta))
    copy!(a.μ, open(io->read_binary_vec(io, T), filename_mean))
    copy!(a.ν, open(io->read_binary_vec(io, T), filename_variance))

    a
end


text:
#=
In this file the core computations for skeletonizing a recurrence matrix
are defined.
=#

# Core function, which gets exported

"""
    skeletonize(R) → R_skel

Skeletonize the `RecurrenceMatrix R` by using the algorithm proposed by
Kraemer & Marwan [^Kraemer2019]. This function returns `R_skel`, a recurrence
matrix, which only consists of diagonal lines of "thickness" one.

[^Kraemer2019]: Kraemer, K.H., Marwan, N. (2019). [Border effect corrections for diagonal line based recurrence quantification analysis measures. Physics Letters A 383(34)](https://doi.org/10.1016/j.physleta.2019.125977).
"""
function skeletonize(X::Union{ARM,SparseMatrixCSC})
    if issymmetric(X)
        symm = true
        # convert lower triangle into a close returns map
        X_cl1 = create_close_returns_map(tril(X))
        # get "vertical" line distribution with position indices
        lines1t, lines1l, lines1c = verticalhisto(X_cl1)
        lines_copy1t = deepcopy(lines1t)
        lines_copy1l = deepcopy(lines1l)
        lines_copy1c = deepcopy(lines1c)

        # create a close returns map with horizontal lines represented by
        # numbers, equal to their lengths
        X_vertical1 = create_close_returns_map(lines_copy1t, lines_copy1l, lines_copy1c, size(X_cl1))

    else
        symm = false
        # convert upper and lower triangles into close returns maps
        X_cl1 = create_close_returns_map(tril(X))
        X_cl2 = create_close_returns_map(sparse(triu(X)'))

        # get "vertical" line distribution with position indices
        lines1t, lines1l, lines1c = verticalhisto(X_cl1)
        lines_copy1t = deepcopy(lines1t)
        lines_copy1l = deepcopy(lines1l)
        lines_copy1c = deepcopy(lines1c)
        lines2t, lines2l, lines2c = verticalhisto(X_cl2)
        lines_copy2t = deepcopy(lines2t)
        lines_copy2l = deepcopy(lines2l)
        lines_copy2c = deepcopy(lines2c)

        # create close returns maps with horizontal lines represented by
        # numbers, equal to their lengths
        X_vertical1 = create_close_returns_map(lines_copy1t, lines_copy1l, lines_copy1c, size(X_cl1))
        X_vertical2 = create_close_returns_map(lines_copy2t, lines_copy2l, lines_copy2c, size(X_cl2))
    end

    # scan the lines, start with the longest one and discard all adjacent lines
    lines_final_t, lines_final_l, lines_final_c = get_final_line_matrix(lines1t, lines1l, lines1c, lines_copy1t, lines_copy1l, lines_copy1c, X_vertical1)

    # if not symmetric input RP, than compute for the upper triangle as well
    if !symm
        lines_final2_t, lines_final2_l, lines_final2_c = get_final_line_matrix(lines2t, lines2l, lines2c, lines_copy2t, lines_copy2l, lines_copy2c, X_vertical2)
        # build RP based on the histogramm of the reduced lines
        X_new = build_skeletonized_RP(lines_final_t, lines_final_l, lines_final_c, lines_final2_t, lines_final2_l, lines_final2_c, size(X_vertical1,1), size(X_vertical1,2))
    else
        # build RP based on the histogramm of the reduced lines
        X_new = build_skeletonized_RP(lines_final_t, lines_final_l, lines_final_c, size(X_vertical1,1), size(X_vertical1,2))
    end
    return X_new
end

# Auxiliary functions

# Transforms the standard RP into a close returns map
function create_close_returns_map(R::SparseMatrixCSC; triangle::Bool = true)
    nr = size(R, 1)
    nrplus = nr+1
    rowvalues = rowvals(R)
    diagvalues = nrplus .+ colvals(R) .- rowvalues # LOI at nrplus (nr+1)
    if triangle
        lower = (diagvalues .<= nrplus)  # indices of the lower triangle
        return sparse(rowvalues[lower], diagvalues[lower], trues(count(lower)), nr, nrplus)
    else
        return sparse(rowvalues, diagvalues, trues(nnz(R)), nr, 2nr+1)
    end
end
# from the 3 vectors, which store all horizontal lines in the reverted RP,
# (close returns map) [obtained by calling first `create_close_returns_map` on a
# recurrence Matrix R, and consecutevely `verticalhisto(R)`] construct and return
# a "dummy" close returns map, where all lines are encoded according to their length.
function create_close_returns_map(lines1, lines2, lines3, dims)
    nzv = sum(lines1)
    rowvalues = zeros(Int, nzv)
    colvalues = zeros(Int, nzv)
    lenvalues = zeros(Int, nzv)
    startspan = 0
    for (i, val) in enumerate(lines1)
        span = startspan .+ (1:val)
        colvalues[span] .= lines3[i]
        rowvalues[span] .= lines2[i]-1 .+ (1:val)
        lenvalues[span] .= val
        startspan += val
    end
    return sparse(rowvalues, colvalues, lenvalues, dims...)
end

# Transforms the reverted RP (close returns map) into a normal RP
function revert_close_returns_map(R::SparseMatrixCSC; triangle::Bool = true)
    nr = size(R, 1)
    rowvalues = rowvals(R)
    columnvalues = rowvalues .+ colvals(R) .- (nr + 1)
    if triangle
        lower = (columnvalues .<= rowvalues)
        return sparse(rowvalues[lower], columnvalues[lower], trues(count(lower)), nr, nr)
    else
        return sparse(rowvalues, columnvalues, trues(nnz(R)), nr, nr)
    end
end

# deletes a line, specified in 'l_vec' (line vector, with first line being
# the total line length, the second line the line-index of the starting point
# and the third line the column-index of the starting point) from the close
# returns "RP" 'R'.
function delete_line_from_cl_ret_RP!(R::Union{ARM,SparseMatrixCSC}, line1::Int, line2::Int, line3::Int)
    R[line2 .+ (1:line1).-1 , line3] .= 0
end

# build the skeletonized RP from the line-Matrix
function build_skeletonized_RP(lines1::Vector{Int}, lines2::Vector{Int}, lines3::Vector{Int}, N::Int, M::Int)

    X_cl_new = spzeros(Bool, N, M)
    # fill up close returns map with lines stored in the new line matrix
    @inbounds for i in eachindex(lines1)
        for j in 1:lines1[i]
            X_cl_new[lines2[i]+j-1, lines3[i]] = true
        end
    end
    XX = revert_close_returns_map(X_cl_new) # revert this close returns map into a legal RP
    X_new = XX .+ XX'
    X_new[diagind(X_new)] .= true # LOI

    return X_new
end
function build_skeletonized_RP(lines1::Vector{Int}, lines2::Vector{Int}, lines3::Vector{Int},
    lines11::Vector{Int}, lines22::Vector{Int}, lines33::Vector{Int}, N::Int, M::Int)

    X_cl_new = spzeros(Bool, N, M)
    X_cl2_new = spzeros(Bool, N, M)
    # fill up close returns map with lines stored in the new line matrix
    @inbounds for i in eachindex(lines1)
        for j in 1:lines1[i]
            X_cl_new[lines2[i]+j-1, lines3[i]] = true
        end
    end
    @inbounds for i in eachindex(lines11)
        for j in 1:lines11[i]
            X_cl2_new[lines22[i]+j-1,lines33[i]] = true
        end
    end

    XX = revert_close_returns_map(X_cl_new) # revert this close returns map into a legal RP
    XXX= revert_close_returns_map(X_cl2_new)
    X_new = XX .+ XXX'
    X_new[diagind(X_new)] .= true # LOI

    return X_new
end

function get_final_line_matrix(lines1t::Vector{Int}, lines1l::Vector{Int},
    lines1c::Vector{Int}, lines_copy1t::Vector{Int}, lines_copy1l::Vector{Int},
    lines_copy1c::Vector{Int}, X_vertical1::SparseMatrixCSC)

    N, M = size(X_vertical1)
    # initialize final lines
    lines_final_t = Int[]
    lines_final_l = Int[]
    lines_final_c = Int[]

    Nlines1 = length(lines1t) # number of found lines

    # go through all lines stored in the sorted line matrix
    @inbounds for l_ind in 1:Nlines1
        # check if line is still in the rendered line matrix
        common_ind = intersect(findall(x-> x==true, lines1t[l_ind] .== lines_copy1t),
                        findall(x-> x==true, lines1l[l_ind] .== lines_copy1l),
                        findall(x-> x==true, lines1c[l_ind] .== lines_copy1c))
        isempty(common_ind) ? continue : nothing
        # get index pair for start of the line
        linei, columni = lines1l[l_ind], lines1c[l_ind]
        # copy this line in the final line matrix
        push!(lines_final_t, lines1t[l_ind])
        push!(lines_final_l, lines1l[l_ind])
        push!(lines_final_c, lines1c[l_ind])
        # delete this line from the RP
        delete_line_from_cl_ret_RP!(X_vertical1, lines1t[l_ind], lines1l[l_ind], lines1c[l_ind])
        # go along each point of the line and check for neighbours
        l_max = lines1t[l_ind]

        @inbounds for l in 1:l_max
            # scan each line twice - above and underneth
            for index in -1:2:1
                # make sure not to exceed RP-boundaries
                (columni+index > M || columni+index == 0) ? break : nothing
                # if there is a neighbouring point, call recursive scan-function
                (X_vertical1[linei+l-1, columni+index] != 0) ? scan_lines!(X_vertical1, lines_copy1t,
                                lines_copy1l, lines_copy1c, linei+l-1, columni+index) : nothing
            end
        end
    end
    return lines_final_t, lines_final_l, lines_final_c
end

# in-place extension of diff, adding diffs of initial and final points
function extendeddiff!(d, x::AbstractVector{T}) where T
    n = length(x)
    d[1] = prev = x[1]
    @inbounds for i in 2:n
            next = x[i]
            d[i] = next-prev
            prev = next
    end
    d[n+1] = -prev
    return d
end

# compute and store all vertical lines from a converted recurrence matrix
# (after applying `create_close_returns_map` to an ARM). The
# line lengths of each found line are stored along with the corresponding starting
# indices of the row and column.
function verticalhisto(R::SparseMatrixCSC)
    rows = rowvals(R)
    # early return
    isempty(rows) && return (Int[], Int[], Int[])
    cols = colvals(R)
    n = length(rows)
    lengthvector = [0]
    columnvector = cols[1:1]
    startvector = rows[1:1]
    # Iterate over columns
    cprev = cols[1]
    r1 = rows[1]
    rprev = r1 - 1 # coerce that (a) is not hit in the first iteration
    @inbounds for i in 1:n
        r = rows[i]
        c = cols[i]
        # Search the second and later segments in the column
        if c == cprev
            if r-rprev != 1 # (a): there is a separation between rprev and r
                # update histogram of segments
                extend_skeleton_histogram!(lengthvector, startvector, columnvector, r, c, rprev-r1+1)
                r1 = r # update the start of the next segment
            end
            rprev = r  # update the previous position
        else # hit in the first point of a new column
            # process the last fragment of the previous column
            extend_skeleton_histogram!(lengthvector, startvector, columnvector, r, c, rprev-r1+1)
            # initialize values for searching new fragments
            cprev = c
            r1 = r
            rprev = r
        end
    end
    # Process the latest fragment
    lengthvector[end] = rprev-r1+1

    ordered = sortperm(lengthvector; rev=true)
    return lengthvector[ordered], startvector[ordered], columnvector[ordered]
end

@inline function extend_skeleton_histogram!(lengthvector, startvector, columnvector, r, c, p)
    lengthvector[end] = p
    push!(lengthvector, 0)
    push!(startvector, r)
    push!(columnvector, c)
end

# manipulates XX inplace and returns a view on l_vec
function scan_lines!(XX::SparseMatrixCSC, l_vec1::AbstractVector{<:Integer}, l_vec2::AbstractVector{<:Integer}, l_vec3::AbstractVector{<:Integer}, line::Integer, column::Integer)

    # for the input index tuple look for the start indices
    index = 0
    local del_ind::Int
    while true
        # check whether the input index tuple is a listed index for starting
        # points of line lengths in the line matrix
        loc_line = findall(column .== l_vec3)
        found = findfirst(isequal(line+index), l_vec2[loc_line])
        if found !== nothing
            del_ind = loc_line[found]
            break
        end
        index -= 1
    end
    # delete the line from RP
    delete_line_from_cl_ret_RP!(XX, l_vec1[del_ind], l_vec2[del_ind], l_vec3[del_ind])

    len, li , co =  l_vec1[del_ind], l_vec2[del_ind], l_vec3[del_ind]

    #delete the line from the line matix, i.e. set to zero
    deleteat!(l_vec1, del_ind)
    deleteat!(l_vec2, del_ind)
    deleteat!(l_vec3, del_ind)

    # check for borders of the RP
    for i in 1:len
        newli, newco = neighborindices(XX, i, li, co)
        if newli != 0 && newco != 0
            scan_lines!(XX, l_vec1, l_vec2, l_vec3, newli, newco)
        end
    end
end

function neighborindices(XX, i::T, li::T, co::T) where T<:Integer
    N, M = size(XX)
    for newli in (li+i .+ (-2:0)), newco in (co-1, co+1)
        if (1 ≤ newli ≤ N) && (1 ≤ newco ≤ M)
            (XX[newli, newco] != 0) && return (newli, newco)
        end
    end
    return zero(T), zero(T) # default if the conditions are not met
end


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

# abstract types
mutable struct SteadyStateProcess <: AbstractMarkovProcess end

# properties
const steady_state_process_props = [:renewal, :markovian, :steady_state]

# parameters
function set_rates!(std::AbstractSTD, cls::AbstractMarkovProcess)
    for nt in transitions(std)
        if !has_prop(std, nt, :rate) && has_prop(std, nt, :distr)
            if get_prop(std, nt, :distr) isa Exponential
                set_prop!(std, nt, :rate, rate(get_prop(std, nt, :distr)))
    end end end
end
function set_markov_chain_matrix!(std::AbstractSTD, cls::SteadyStateProcess)
    P = zeros(_MSM.Measurement, ns(std), ns(std))
    for nt in transitions(std)
        rate = get_prop(std, nt, :rate)
        P[_LG.src(nt),_LG.dst(nt)] =
            ifelse(_UF.unit(rate)==_UF.NoUnits, rate, rate * 1.0u"s" |> u"s/s")
    end
    for ns in 1:ns(std)
        P[ns,ns] = 1.0 - sum(P[ns,:])
    end
    set_prop!(std, :P, ifelse(all(_MSM.uncertainty.(P).== 0.0_UF.unit(first(P))),
                              _MSM.value.(P),
                              P))
end
function set_parameters!(std::AbstractSTD, cls::SteadyStateProcess)
    set_rates!(std, cls)
    set_markov_chain_matrix!(std, cls)
end

# stochastic process
solve_steady_state(std) =
    hcat(get_prop(std, :P) .- _LA.I(ns(std)),ones(ns(std)))' \ [zeros(ns(std))..., 1.0]
function solve!(std::AbstractSTD, cls::SteadyStateProcess)
    set_parameters!(std, cls)

    set_prop!(std, :cls, cls)
    set_prop!(std, :time, [Inf])
    set_prop!(std, states(std), :prob, solve_steady_state(std))

    set_info!(std, :solved, true)
end


text:
################################################################################
#
#  UnsafeRationals
#
################################################################################

# The ordinary Rational{Int} checks for overflow in every operations.
# UnsafeRational is like Rational, but with unsafe operations, that is, without
# the overflow checking.

import Base: Int128

################################################################################
#
#  Constructors and conversions
#
################################################################################

function UnsafeRational{T}(x::Int) where {T}
  return UnsafeRational{T}(T(x), one(T))
end

function UnsafeRational{T}(x::T) where {T}
  return UnsafeRational{T}(x, one(T))
end

function UnsafeRational{T}(x::Rational{<: Integer}) where {T}
  return UnsafeRational{T}(T(x.num), T(x.den))
end

function Int(x::UnsafeRational{T}) where {T}
  return Int(x.num)
end

Rational(x::UnsafeRational) = x.num//x.den

function Base.Int128(x::fmpz)
  return Base.Int128(BigInt(x))
end

################################################################################
#
#  Printing
#
################################################################################

function Base.show(io::IO, x::UnsafeRational{T}) where T
  print(io, UnsafeRational{T}, "(", x.num, "//", x.den, ")")
end

################################################################################
#
#  Canonicalization functions
#
################################################################################

@inline function divgcd(x::Integer, y::Integer)
  g = gcd(x, y)
  return div(x, g), div(y, g)
end

@inline function divgcd(x::Int,y::Int)
  #g = gcd(x, y)
  if x < 0
    xx = -x
  else
    xx = x
  end
  if y < 0
    yy = -y
  else
    yy = y
  end
  g = ccall((:n_gcd, libflint), UInt, (UInt, UInt), xx, yy)
  return div(x, g), div(y, g)
end

################################################################################
#
#  Multiplication
#
################################################################################

@inline function Base.:(*)(x::UnsafeRational{T}, y::UnsafeRational{T}) where T
  xn,yd = divgcd(x.num,y.den)
  xd,yn = divgcd(x.den,y.num)
  z = UnsafeRational(xn * yn, xd * yd)
  #@assert Rational(x) * Rational(y) == Rational(z)
  return z
end

################################################################################
#
#  Addition
#
################################################################################

@inline function Base.:(+)(x::UnsafeRational{T}, y::UnsafeRational{T}) where T
  xd, yd = divgcd(x.den, y.den)
  s, t = x.num * yd + y.num * xd, x.den * yd
  sq, tq = divgcd(s, t)
  #@assert Rational(x) + Rational(y) == sq//tq
  return UnsafeRational{T}(sq, tq)
end

################################################################################
#
#  Division
#
################################################################################

@inline function Base.:(//)(x::UnsafeRational{T}, y::UnsafeRational{T}) where {T}
  xn,yn = divgcd(x.num,y.num)
  xd,yd = divgcd(x.den,y.den)
  z = UnsafeRational{T}(xn * yd, xd * yn)
  return z
end

################################################################################
#
#  Rounding
#
################################################################################

@inline function Base.floor(x::UnsafeRational{T}) where {T}
  z = UnsafeRational{T}(fld(x.num, x.den))
  return z
end

@inline function Base.ceil(x::UnsafeRational{T})where {T}
  z = UnsafeRational{T}(cld(x.num, x.den))
  return z
end

################################################################################
#
#  Ad hoc unary and binary operations
#
################################################################################

@inline function Base.:(*)(x::UnsafeRational{T}, y::Int) where {T}
  z = T(y)
  zn, zd = divgcd(x.den, z)
  z = UnsafeRational{T}(x.num * zd, zn)
  return z
end

@inline function Base.:(-)(x::UnsafeRational{T}, y::Int) where {T}
  z = x - UnsafeRational{T}(y)
  return z
end

@inline function Base.:(+)(x::UnsafeRational{T}, y::Int) where {T}
  z = x + UnsafeRational{T}(y)
  return z
end

@inline function Base.:(-)(x::UnsafeRational{T}, y::UnsafeRational{T}) where {T}
  z = x + (-y)
  return z
end

@inline function Base.:(-)(x::UnsafeRational{T}) where {T}
  z = UnsafeRational{T}(-x.num, x.den)
  return z
end

@inline isnegative(x::UnsafeRational) = x.num < 0

################################################################################
#
#  Comparison
#
################################################################################

@inline function Base.:(isless)(x::UnsafeRational{T}, y::Int) where {T}
  z = (x.num < x.den * y)
  return z
end


text:
function toh04()
#
#  generalization of K.C. Toh and L. N. Trefethen
#
    p = poly((10/11)*ones(20))
    z = [10/11 20]

    p, PolyZeros(z)
end


text:
#!/usr/bin/env julia

using Filetimes

# Command line programming for converting DateTimes and Filetimes,
# e.g.
# ./filetime.jl 2018-07-19T18:26:00.766
# ./filetime.jl 131764980100450000
# ./filetime.jl

function main(args)

    if length(args) > 0
        x = args[1]
        y = tryparse(Int64,x)
        if isnull(y)
            println("$x\t$(filetime(x))")
        else
            println("$(datetime(get(y)))\t$x")
        end
    else
        x = now()
        println("$x\t$(filetime(x))")
    end
end

main(ARGS)


text:
#!/usr/bin/env julia

# generate pictures of builtin data
# usage: julia draw-svg-swatches.jl > somewhere/colorschemes.svg

using Colors, ColorSchemes

function draw_swatches()
    imagewidth, imageheight = 800, 3500
    margin = 20
    swatch_height = 15
    # svg information sent to STDIO
    # write svg header
    println("""<svg xmlns="http://www.w3.org/2000/svg" version="1.1" width="$(imagewidth)" height="$(imageheight)"> \n
<rect x="0" y="0" width="$(imagewidth)" height="$(imageheight)" fill="white" /> \n """)
    # title
    println("""<text x="$(margin)" y="$(margin - margin/2)" font-size="14px" font-family="Helvetica-Bold">Colorschemes</text>\n""")
    files = readdir()
    y   = 50
    for f in filter(f -> ismatch(Regex("^[a-z]+\.txt"), f), readdir(Pkg.dir("ColorSchemes", "data")))
        cs = loadcolorscheme(f)
        w = (imagewidth - margin - margin - 50) / length(cs)
        println("""<text x="$(margin)" y="$(y-5)" fill="black" font-size="10px" font-family="Futura">$f</text>\n""")
        x = margin
        #swatch
        for (i, p) in enumerate(cs)
            r, g, b = convert(Int, round(256 * p.r)), convert(Int, round(256 * p.g)), convert(Int, round(256 * p.b))
            println("""<rect x="$(x)" y="$(y)" width="$(w)" height="$(swatch_height)"  stroke="none" fill="rgb($r, $g, $b)" />\n""")
            x = margin + (i * w)
        end
        y += 20
        #blend
        stepping = 0.001
        printwidth = (imagewidth - margin - margin - 50)
        if length(cs) == 0
            continue
        end
        for i in 0:stepping:1
            c = colorscheme(cs, i)
            r, g, b = convert(Int, round(256 * c.r)), convert(Int, round(256 * c.g)), convert(Int, round(256 * c.b))
            x = margin + (i * printwidth)
            println( """<rect x="$(x)" y="$(y)" width="$(1 + printwidth * stepping)" height="$(swatch_height)" stroke="none" fill="rgb($r, $g, $b)" />\n""")
        end
        y += 35
    end
    println("</svg>")
end

draw_swatches()


text:
import DPFEHM
import GaussianRandomFields
import Optim
import PyPlot
import Random
import RegularizationDP
import Zygote

Random.seed!(0)
function getobsnodes(coords, obslocs)
	obsnodes = Array{Int}(undef, length(obslocs))
	for i = 1:length(obslocs)
		obsnodes[i] = findmin(map(j->sum((obslocs[i] .- coords[:, j]) .^ 2), 1:size(coords, 2)))[2]
	end
	return obsnodes
end
mins = [0, 0]#meters
maxs = [100, 100]#meters
ns = [101, 101]
num_eigenvectors = 200
x_true = randn(num_eigenvectors)
x0 = zeros(num_eigenvectors)
sqrtnumobs = 16
obslocs_x = range(mins[1], maxs[1]; length=sqrtnumobs + 2)[2:end - 1]
obslocs_y = range(mins[2], maxs[2]; length=sqrtnumobs + 2)[2:end - 1]
obslocs = collect(Iterators.product(obslocs_x, obslocs_y))[:]
observations = Array{Float64}(undef, length(obslocs))
dist(x, y) = sum((x .- y) .^ 2)

coords, neighbors, areasoverlengths, _ = DPFEHM.regulargrid2d(mins, maxs, ns, 1.0)
Qs = zeros(size(coords, 2))
boundaryhead(x, y) = 5 * (x - maxs[1]) / (mins[1] - maxs[1])
dirichletnodes = Int[]
dirichleths = zeros(size(coords, 2))
for i = 1:size(coords, 2)
	if coords[1, i] == mins[1] || coords[1, i] == maxs[1] || coords[2, i] == mins[2] || coords[2, i] == maxs[2]
		push!(dirichletnodes, i)
		dirichleths[i] = boundaryhead(coords[1:2, i]...)
	end
	if coords[1, i] == 50 && coords[2, i] == 50
		Qs[i] = 1.0
	end
end

lambda = 50.0#meters -- correlation length of log-conductivity
sigma = 1.0#standard deviation of log-conductivity
mu = 0.0#mean of log-permeability can be arbitrary because this is steady-state and there are no fluid sources

cov = GaussianRandomFields.CovarianceFunction(2, GaussianRandomFields.Matern(lambda, 1; σ=sigma))
x_pts = range(mins[1], maxs[1]; length=ns[1])
y_pts = range(mins[2], maxs[2]; length=ns[2])
grf = GaussianRandomFields.GaussianRandomField(cov, GaussianRandomFields.KarhunenLoeve(num_eigenvectors), x_pts, y_pts)
logKs = GaussianRandomFields.sample(grf)
parameterization = copy(grf.data.eigenfunc)
sigmas = copy(grf.data.eigenval)

logKs2Ks_neighbors(Ks) = exp.(0.5 * (Ks[map(p->p[1], neighbors)] .+ Ks[map(p->p[2], neighbors)]))
function solveforh(logKs, dirichleths)
	@assert length(logKs) == length(Qs)
	if maximum(logKs) - minimum(logKs) > 25
		return fill(NaN, length(Qs))#this is needed to prevent the solver from blowing up if the line search takes us somewhere crazy
	else
		Ks_neighbors = logKs2Ks_neighbors(logKs)
		return DPFEHM.groundwater_steadystate(Ks_neighbors, neighbors, areasoverlengths, dirichletnodes, dirichleths, Qs)
	end
end
solveforheigs(x) = solveforh(x2logKs(x), dirichleths)
x2logKs(x) = reshape(parameterization * (sigmas .* x), ns...)
logKs_true = x2logKs(x_true)
h_true = solveforheigs(x_true)

obsnodes = getobsnodes(coords, obslocs)
obssigma = 1e-3
observations .= h_true[obsnodes]#set up the observations
f(logKs) = sum((solveforh(logKs, dirichleths)[obsnodes] - observations) .^ 2 ./ obssigma ^ 2)

#do the optimization
options = Optim.Options(iterations=250, extended_trace=false, store_trace=true, show_trace=true, x_tol=1e-6)
logKs_est, opt = RegularizationDP.optimize(f, x2logKs, x->sum(x .^ 2), zeros(num_eigenvectors), options)
#plot results from the inverse analysis
fig, axs = PyPlot.subplots(2, 3, figsize=(12, 8))
ims = axs[1].imshow(logKs_true)
axs[1].title.set_text("True Conductivity")
fig.colorbar(ims, ax=axs[1])
ims = axs[2].imshow(logKs_est, vmin=minimum(logKs_true), vmax=maximum(logKs_true))
axs[2].title.set_text("Estimated Conductivity")
fig.colorbar(ims, ax=axs[2])
ims = axs[3].imshow(reshape(solveforh(logKs_true, dirichleths), ns...))
axs[3].title.set_text("True Head")
fig.colorbar(ims, ax=axs[3])
ims = axs[4].imshow(reshape(solveforh(logKs_est, dirichleths), ns...))
axs[4].title.set_text("Estimated Head")
fig.colorbar(ims, ax=axs[4])
axs[5].plot([0, 5], [0, 5], "k", alpha=0.5)
axs[5].plot(observations, solveforh(logKs_est, dirichleths)[obsnodes], "k.")
axs[5].set_xlabel("Observed Head")
axs[5].set_ylabel("Predicted Head")
axs[6].semilogy(map(i->opt.trace[i].value, 1:length(opt.trace)))
axs[6].set_xlabel("Iteration")
axs[6].set_ylabel("Loss")
fig.tight_layout()
display(fig)
println()
PyPlot.close(fig)


text:
#! /usr/bin/env julia --project --color=yes

using Base.Filesystem

docsDir = abspath(dirname(@__FILE__))
baseDir = abspath(joinpath(docsDir, "../"))

if !in(LOAD_PATH,baseDir)
  push!(LOAD_PATH,baseDir)
end

using Documenter, Paxos

makedocs(
  sitename="Paxos.jl",
  pages = [
    "index.md",
    "Ballots" => "ballots.md",
    "Clients" => "clients.md",
    "Configurations" => "configurations.md",
    "Leaders" => "leaders.md",
    "Ledgers" => "ledgers.md",
    "Nodes" => "nodes.md",
    "Protocols" => "protocols.md",
    "Transports" => [
      "Common Definitions" => "transports_common.md",
      "Memory" => "transports_memory.md",
      "TCP" => "transports_tcp.md"
    ],
    "Utils" => "utils.md"
  ]
  )


text:
## basic type declarations and utilities
"""
A categorized string, supertype of all token types in this module.

In addition to the AbstractString API, an AbstractToken has a
category which roughly classifies its content.
It is accessed by [`category`](@ref)(t::AbstractToken).
Token categories are typed as [`Nibble`](@ref) and internally stored in 4 bits.
For highest performance and still high readability, use named constants
in your application to indicate semantic meaning in your context.
Even better for readability and debugging is the use of Enum types for
token category. Have a look at the quite generic [`TCategory`](@ref)
as a starting point

For comparisons and hashing, tokens are treated as strings. Two tokens with
equal content but different category are treated equal and have the
same hash value.

# Missing support

Some Missing and Nothing support is already built into all AbstractToken subtypes:
they have an internal encoding for missing and nothing, and all concrete subtypes
must have a constructor with one parameter of type Missing resp. Nothing,
which construct a token with that encoding.
ismissing(t::AbstractToken) and isnothing(t::AbstractToken) will return true for this encoding,
and Tests with == and === will compare missing/nothing against these encodings.

This is sufficient for most applications. For vectors of tokens, you can explicitly specify
a subtype of Union{Token,Missing,Nothing} for element types, enabling explicit
missing and/or nothing values on array read access (internal encodings are
converted to missing or nothing).

# string literals

For the three most important token types  [`DToken`](@ref),
[`BToken`](@ref) and  [`HToken`](@ref), there are special string literal
macros implemented. They use the first letter D,B or H plus a capitalited
hexadecimal character as string prefix to define type and category.
E.g.: D2"hello" is equivalent to DToken(Nibble(2),"hello"), and
HB"hello" will create HToken(Nibble(11),"hello").

# token properties

A token t :: T <: AbstractToken supports the following readonly properties,
implemented by a getproperty method definition

    * t.cat :: Nibble category of token, === category(t)

    * t.ofs :: UInt32 offset into internal buffer, === offset(t)

    * t.len :: UInt64 number of code units, === usize(t)

# interface requirements

Tokens support the AbstractString API. Many methods are defined
with specialized, efficient code. Developing a new AbstractToken subtype
requires overloading of many of them.

Every implementation must use Utf8 encoding: token methods operate
mostly on Utf8 code units and rely on its properties.

There is a Token specific interface besides AbstractString API, in the
source code it is characterized by chapter comments

# implementations in this module

This module supplies a very memory efficient implementation
for very short tokens with [`DToken`](@ref), a
short string plus category in 8 bytes. It can act as a
substitute for short String instances, avoiding any allocation and
indirect (pointer) access. Data locality is improved, gaining further speed
advantages by better CPU cache use.

[`BToken`](@ref) stores larger strings, but needs an additional code
unit buffer. Tokens can share the same buffer, reducing overhead for heap
management and even re-use content, like SubString does.

There is a technical size limit for tokens: at most 2^27-1 code units.
It applies to *all* AbstractToken implementations. It might restrict usage in
some general text processing tasks, eg processing huge log files.
Recommended countermeasure is processing chunks of some KByte to MByte,
using [`IOShared`](@ref). Tokenizing a text, usually leads to small fragments
of a couple of bytes. For tokens in the sense of "atoms of text", a size in
the range of a couple of bytes is expected. Text written by humans, like
source code and published books, has sizes up to some MByte.

[`HToken`](@ref) is a combination of DToken and BToken in
one type, using no buffer for code unit counts below 8.

# related types and APIs in this module

[`IOShared`](@ref) is simular to IOBuffer but designed to interact smoothly
with tokens. Though mutable, its buffer can safely be shared with tokens
(as well as SubString instances).

[`TokenVector`](@ref) acts as vector of tokens, using one shared buffer for
all vector elements. It uses significantly less memory than a Vector{String},
but might require copying token contents into the buffer used by the TokenVector
when a token is put into the vector.

[`TokenTree`](@ref) implements a simple tree structure upon a TokenVector.

[`Lexer`](@ref) implementations have functions for lexical analysis and
tokenizing of text into tokens
"""
abstract type AbstractToken <: AbstractString
end

"""
strings known to have codeunit type UInt8 and Utf8-Encoding.

They support the AbstractString API and additionally
the methods [`usize`](@ref) and [`offset`](@ref)
"""
const Utf8String = Union{String,SubString{String},AbstractToken}

offset(s::String) = 0

offset(s::SubString) = s.offset

#########################################################
################# AbstractToken API #####################
#########################################################

"""
    offset(t::AbstractToken)

number of code units in the buffer associated with an token
before the first code unit belonging to the token.

If code units are stored directly in a FlyToken, offset is 0.

Buffer is not needed in this function, no consistency or bounds
checks are performed.
"""
function offset(t::AbstractToken)
    throw(MethodError(offset, (t,)))
end

"""
    isDirecty(t::AbstractToken) -> Bool

true if t stores its code units directly in its FlyToken.
"""
function isdirect(t::AbstractToken)
    throw(MethodError(isdirect, (t,)))
end

"""
    category(t::AbstractToken) -> Nibble

Current category of given token. A value in 0:15.
Meaning depends on context.

"""
function category(t::AbstractToken)::Nibble
    throw(MethodError(category, (t,)))
end

# Token category default definitions
# use is optional, except in Token constructors
# categories are ordered by groups.
# category constants start with "T_".
# Names in in Tokens which begin with "T_"
# must be a category constant

## Category group 1: used in Token constructors and parsers

"""
# T_END = 0

A token signalling some end of some structure.

Token constructors with argument *nothing* create a token of this category
and empty contents.

An empty T_END token can be returned if an attempt is made to process
beyond end of data.

It may contain a string, e. g. the end of a node
in XML, or an end-of-line sequence.

There are grammars, where a structure end has a lexical
representation and could be identified by a lexer. Common examples are tokens like
"end", ")", "}", "]", ";". But the end of a semantical sequence can also be
defined by context data, like indentation level, or operator precedence rules,
which have no lexical characterization. In such cases, the parser could generate T_END
tokens.

In many grammars, whitespace always terminates a token, and is allowed between tokens.
Usually, lexers will skip whitespace. If you need a lexer which returns whitespace
as a token, consider using T_END tokens for it.

In a  [`TokenTree`](@ref), empty T_END tokens are used internally to terminate a list of children.
For this reason, empty T_END tokens are not allowed in a TokenTree as a node.
"""
const T_END = Nibble(0)

"""
# T_INT = 1

An integer value with decimal encoding. Details depend on the lexer: it could be restricted to
a sequence of digits, or allow a leading sign, like '+' or '-', and even group separators
like ',' in English locales.

Token constructors with an Integer argument create a token of category
T_INT and its decimal string representation as content.
"""
const T_INT = Nibble(1)

"""
## T_REAL = 2

An optionally signed number with a decimal separator and/or decimal exponent.

Token constructors with a Real argument create a token of category
T_REAL and its default string representation as content.
"""
const T_REAL = Nibble(2)

"""
## T_CHAR = 3

Character token. Token content is usually one Utf8-coded character.
Typical lexers recognize a string enclosed in single quotes as category T_CHAR,
and will usually remove the leading and trailing quotes.

Token constructors with an AbstractChar argument create a token of category
T_CHAR and its Utf8 string representation as content.
"""
const T_CHAR = Nibble(3)

"""
## T_TEXT = 4

Some text which was not (fully) tokenized according to the preceding token categories.
This is the default category for tokens used in general string processing.
Example: text entities in XML.

A token constructor with a single AbstractString argument (which is no token)
creates a token of category T_TEXT.

"""
const T_TEXT = Nibble(4)

"""
## T_IDENT = 5

An identifier in a lexer context.

Typical rules are: 1st character is a letter, following characters are
letters or digits. Some special characters like '_' could also
appear in a T_IDENT token.

Lexers working on byte level might treat all Non-ASCII unicode characters
either as letter or as non-letter.

A token constructor with a single Bool argument creates a token of category T_IDENT
and content "false" or "true".
"""
const T_IDENT= Nibble(5)

"""
## T_SPECIAL = 6

A sequence of special characters, which has no syntactical meaning to the
lexer. Special characters which identify another lexical construct are
not reported as T_SPECIAL. Examples are characters , like quotes,
or are recognized as operator or symbol
, and not recognized as symbol
or operator token.

Lexers decide to report each special character as its own token, or
report a contiguous sequence of special characters as one token.

Token constructors with argument *missing* create a T_SPECIAL token
with empty content.

"""
const T_SPECIAL = Nibble(6)

## Category group 2: used in lexers and parsers

"""
## T_QUOTE = 7

A string enclosed in some form of quotes.

Double quotes are a typical case. Parsers will normally remove the leading
and trailing quotes, and resolve escape sequences. Common escape rules are
doubling quotes inside a quoted string, or the use of a dedicated escape character,
like backslash.

This category is simple to detect by a simple lexer (determined by its 1st character),
but may need post processing. See [`QuoteLA`](@ref) for an example.

"""
const T_QUOTE = Nibble(7)

"""
## T_OP = 8

A sequence starting with a special character like '+' and '-', which has several alternative
lexical continuations to a complete token. The name T_OP shall indicate the use as operator,
and the different case of an option (like in command line arguments).
Typical scenarios for '+' and '-' are:

 * number scenario: a sign is immediately followed by a digit sequence.

 * operator scenario: sign and optionally further special characters define an operator Symbol, e.g. '+='

 * option scenario: a sign followed by an identifier, giving a set or unset option.

Recommendation: a lexer should return the char sequence of one of the given scenarios, and
the parser will test the length and check for the other scenarios if length is 1.

The default lexer definition tries to continue with digits, which is probably the
most common case. If the parser context requires to consider the sign as an operator,
and not part of a number, the parser has to break the raw token into two tokens.
If the raw token has length 1, it implies that no number but an operator has to be parsed.
Alternative: define no continuation bytes for T_OP.

Alternatively, a lexer could return the char sequence of one of the given scenarios, and
the parser will test the length and check for the other scenarios if length is 1.
See [`OpLA`](@ref) for an implementation.

T_OP is the first category of tokens which can have children in a TokenTree.
"""
const T_OP = Nibble(8)

"""
## T_COMMENT = 9

comment. Typically contains the pure comment text, without delimiters.
Delimiters may be accessible via lexer context. Many computer languages support
different comment flavours like inline, rest of line or multiline comments.

Comments can be used to implement extensions to some base language.
They are ignored by the base language parser, but parsed in a later stage.
For this reason, a comment token is allowed to have children in a TokenTree.
"""
const T_COMMENT = Nibble(9)

# group (3): higher level categories, intended for parsers and token trees

"""
## T_KEY = 10

Identifier recognized as keyword.

Lexers may support different nonations for a keyword, and return a unique
'canonical' representation,
e. g. converting SQL keywords to uppercase. If all T_KEY strings have a
canonical representation with less than 8 code units, your application could
restrict T_KEY tokens to type DToken. Even if the grammar has longer keywords,
you can define a short form fitting into a DToken as canonical representation,
to enable use of DToken in the following processing steps.

Reserved words in programming languages are usually tokenized as T_KEY.
In a TokenTree, they can have children, e. g. condition and action for
control structures like IF/THEN/ELSE or FUNCTION parameters and code.
"""
const T_KEY = Nibble(10)

"""
## T_SYM = 11

One or more special characters which form a semantically interpreted
symbol, e.g. "*", ">>>" or "+=". A lexer may accept different notations for the
same symbol, e. g. "<>", "!=" and "≠" for inequality,
and may return a normalized notation.

Symbols usually fit into a DToken. Consider to grant this in your
application (and use DToken for symbols).

In TokenTree, symbols can have children, to reflect the use of symbols
as operators in common computer languages.

A token t of category T_SYM and empty content is regarded a missing
value, i.e. ismissing(t) returns true.

"""
const T_SYM = Nibble(11)

"""
## T_CMD = 12

embedded commands and processing instructions.

In a TokenTree, it may have children which represent
the parsed contents of the instruction.

A lexer will typically identify a processing instruction by some unique prefix
and suffix. In the first parser stage, the text between prefix and suffix may be
stored as unparsed text in a T_CMD token. When building a syntax tree,
it might get parsed (probably with a different parser).

Examples are DTD and PI in XML, embedded javascript code in HTML,
processing instructions in text templates like XSLT.

From a lexer-s view, T_CMD is structurally similar to quoted strings and comments:
a sequence of unparsed content with a delimiter sequence at its begin and end,
which is excluded from the returned token.
"""
const T_CMD = Nibble(12)

"""
## T_REC = 13

record structure: a node in TokenTree with children.

T_REC is recommended if children access is based on a key attribute which
identifies them within the children list, like a field name in a julia struct or the
attribute name in an XML item. Often, children have different data types and
a static list of allowed keys exists.

A specialized lexer may "abuse" T_REC for a lexical structure, e.g. a match
against a certain regular expression, like a date format.

"""
const T_REC = Nibble(13)

"""
## T_LIST = 14

list: a node in TokenTree with children.
T_LIST is recommended if children access is based on the sequence order,
like julia vectors or JSON arrays. Often, all children have the same data type.

A specialized lexer may "abuse" T_LIST for a lexical structure, e.g. a match
against a certain regular expression, like a date format.
"""
const T_LIST = Nibble(14)

"""
## T_EXT = 15

Extension category.

No predefined syntax or semantic, define it application-specific. Can have children in
a TokenTree.
"""
const T_EXT = Nibble(15)

"""
    subtoken(offset::UInt32, size::UInt64, t::T<:AbstractToken)

central method to construct a new token of the same type which
holds a sequence of code units of its source t.

It uses an (offset::UInt32,size::UInt64) pair to specify the
code units in the result: skips offset code units and keeps the
following size code units in a newly constructed token.

In most cases, it is a "view" operation like SubString
construction. Exceptions are very short sequences for token
types which allow for direct encoding without buffer.

Specifying concrete Integer subtypes may look uncommon to
Julia, but is for safety and efficieny: in most token
operations, it is exactly the integer type needed, and having
different types for offset and size reduces the risk that these
parameters are mixed-up: it is difficult enough for the user
to distinguish the common string (first,last) index notation
from the performance-oriented lowlevel (offset,size) code unit
segment notation. The need for explicit conversion of offset and
size hopefully helps to to avoid mix-up.

There is another subtoken method which uses a first:last
character index sequence like usual String and SubString
operations. The (offset,size) oriented operation is faster than
the form subtoken(t,first,last) , but less safe: you can construct
tokens which are no valid strings, because there is no check
whether offset+1 is a valid character index of t or
codeunit(t,offset+size) is the last code unit of a character.
"""
subtoken(offset::UInt32, size::UInt64, t::T) where T<:AbstractToken = T(offset,size,t)

"""
    subtoken(t::AbstractToken, first::Int, last::Int)

Like SubString with same parameters, but returns
a new token of same type with the same category
and content SubString(string(t),first,last).

"""
function subtoken(t::T, first::Int, last::Int) where T<:AbstractToken
    first <= last || return T(t.cat)
    @boundscheck begin
        checkbounds(t, first:last)
        @inbounds isvalid(t, first) || string_index_err(t, first)
        @inbounds isvalid(t, last) || string_index_err(t, last)
    end
    nextlast = nextind(t,last)
    subtoken(UInt32(first-1),(nextlast-first)%UInt64,t)
end

#= try with default implementation, using isvalid ...

function Base.thisind(t::AbstractToken, i::Int)
    if i<=0
        @boundscheck i==0 || boundserr(t,i)
        return 0
    else
        lenp1 = ncodeunits(t)+1
        if i>=lenp1
            @boundscheck i==lenp1 || boundserr(t,i)
            return i
        end
        # now we have a normal case: i is inbound index
        # copy/paste from strings.jl
        @inbounds b = codeunit(t, i)
        (b & 0xc0 == 0x80) & (i-1 > 0) || return i
        @inbounds b = codeunit(t, i-1)
        between(b, 0b11000000, 0b11110111) && return i-1
        (b & 0xc0 == 0x80) & (i-2 > 0) || return i
        @inbounds b = codeunit(t, i-2)
        between(b, 0b11100000, 0b11110111) && return i-2
        (b & 0xc0 == 0x80) & (i-3 > 0) || return i
        @inbounds b = codeunit(t, i-3)
        between(b, 0b11110000, 0b11110111) && return i-3
        return i
    end
end
=#

## Token API definition and implementations which need no specialization

# needs overloading!
"""
    function usize(t) :: UInt64

the size in bytes of some text contents. Applicable to token Types,
Strings, IOShared and more.

returned as an UInt64, which is the data type for token
sizes used throughout all token implementations.

Functions ncodeunits and sizeof are derived from usize
"""
function usize end

usize32(s) = usize(s) % UInt32

"""
returns the category of a token as a Nibble
"""
function category end

## pseudo properties

function Base.getproperty(a::AbstractToken, s::Symbol)
    if s === :ofs
        return offset(a)
    elseif s === :len
        return usize(a)
    elseif s === :cat
        return category(a)
    else
        return Core.getproperty(a, s)
    end
end

## AbstractString API implementations which need no specialization

Base.codeunit(t::AbstractToken) = UInt8

Base.ncodeunits(t::AbstractToken) = t.len%Int

# simplified implementation, assuming t is valid UTF8.
# isvalid only checks current byte. It does not check if all code units
# of the unicode character exist and have valid values.
Base.isvalid(t::AbstractToken, i::Int) = codeunit(t, i) & 0xc0 != 0x80

# Not helpful. Julia Base uses optimized c code with pointers.
# Base.cmp(a::Utf8String, b::Utf8String) = cmp_codeunits(a,b)

#= use default?checkbounds
function Base.==(a::AbstractToken, b::AbstractToken)
    cmp_codeunits(a,b) == 0
end
=#


text:
# Problem evaluation #
# ========================== #
function _eval_first_stage(stochasticprogram::StochasticProgram, x::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    return eval_objective(first_stage.obj, x)
end
function _eval_second_stage(stochasticprogram::TwoStageStochasticProgram, x::AbstractVector, scenario::AbstractScenario, solver::MPB.AbstractMathProgSolver)
    outcome = outcome_model(stochasticprogram, x, scenario, solver)
    solve(outcome)
    return probability(scenario)*getobjectivevalue(outcome)
end
function _eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                             x::AbstractVector,
                             solver::MPB.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    return outcome_mean(outcome_generator, scenarios(stochasticprogram))
end
function _eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                             x::AbstractVector,
                             solver::MPB.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    Qs = Vector{Float64}(undef, nworkers())
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    @sync begin
        for (i,w) in enumerate(workers())
            @async Qs[i] = remotecall_fetch((sp,outcome_generator)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return 0.0
                return outcome_mean(outcome_generator, scenarioproblems.scenarios)
            end,
            w,
            stochasticprogram.scenarioproblems[w-1],
            outcome_generator)
        end
    end
    return sum(Qs)
end
function _stat_eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                                  x::AbstractVector,
                                  solver::MPB.AbstractMathProgSolver) where {S, SP <: ScenarioProblems}
    N = nscenarios(stochasticprogram)
    outcome_generator = scenario -> outcome_model(stochasticprogram, x, scenario; solver = solver)
    𝔼Q, σ² = outcome_welford(outcome_generator, scenarios(stochasticprogram))
    return 𝔼Q, sqrt(σ²)
end
function _stat_eval_second_stages(stochasticprogram::TwoStageStochasticProgram{S,SP},
                                  x::AbstractVector,
                                  solver::MPB.AbstractMathProgSolver) where {S, SP <: DScenarioProblems}
    N = nscenarios(stochasticprogram)
    partial_welfords = Vector{Tuple{Float64,Float64,Int}}(undef, nworkers())
    @sync begin
        for (i,w) in enumerate(workers())
            @async partial_welfords[i] = remotecall_fetch((sp,stage_one_generator,stage_two_generator,stage_one_params,stage_two_params,x,solver)->begin
                scenarioproblems = fetch(sp)
                isempty(scenarioproblems.scenarios) && return zero(eltype(x)), zero(eltype(x))
                    outcome_generator = scenario -> begin
                        outcome_model = Model(solver = solver)
                        _outcome_model!(outcome_model,
                                        stage_one_generator,
                                        stage_two_generator,
                                        stage_one_params,
                                        stage_two_params,
                                        x,
                                        scenario)
                        return outcome_model
                   end
                   return (outcome_welford(outcome_generator, scenarioproblems.scenarios)..., length(scenarioproblems.scenarios))
            end,
            w,
            stochasticprogram.scenarioproblems[w-1],
            stochasticprogram.generator[:stage_1_vars],
            stochasticprogram.generator[:stage_2],
            stage_parameters(stochasticprogram, 1),
            stage_parameters(stochasticprogram, 2),
            x,
            solver)
        end
    end
    𝔼Q, σ², _ = reduce(aggregate_welford, partial_welfords)
    return 𝔼Q, sqrt(σ²)
end
function _eval(stochasticprogram::StochasticProgram{2}, x::AbstractVector, solver::MPB.AbstractMathProgSolver)
    xlength = decision_length(stochasticprogram)
    length(x) == xlength || error("Incorrect length of given decision vector, has ", length(x), " should be ", xlength)
    all(.!(isnan.(x))) || error("Given decision vector has NaN elements")
    cᵀx = _eval_first_stage(stochasticprogram, x)
    𝔼Q = _eval_second_stages(stochasticprogram, x, solver)
    return cᵀx+𝔼Q
end
# Mean/variance calculations #
# ========================== #
function outcome_mean(outcome_generator::Function, scenarios::Vector{<:AbstractScenario})
    Qs = zeros(length(scenarios))
    for (i,scenario) in enumerate(scenarios)
        outcome = outcome_generator(scenario)
        status = solve(outcome)
        if status != :Optimal
            error("Outcome model could not be solved, returned status: $status")
        end
        Qs[i] = probability(scenario)*getobjectivevalue(outcome)
    end
    return sum(Qs)
end
function outcome_welford(outcome_generator::Function, scenarios::Vector{<:AbstractScenario})
    Q̄ₖ = 0
    Sₖ = 0
    N = length(scenarios)
    for k = 1:N
        Q̄ₖ₋₁ = Q̄ₖ
        outcome = outcome_generator(scenarios[k])
        status = solve(outcome)
        if status != :Optimal
            error("Outcome model could not be solved, returned status: $status")
        end
        Q = getobjectivevalue(outcome)
        Q̄ₖ = Q̄ₖ + (Q-Q̄ₖ)/k
        Sₖ = Sₖ + (Q-Q̄ₖ)*(Q-Q̄ₖ₋₁)
    end
    return Q̄ₖ, Sₖ/(N-1)
end
function aggregate_welford(left::Tuple, right::Tuple)
    x̄ₗ, σₗ², nₗ = left
    x̄ᵣ, σᵣ², nᵣ = right
    δ = x̄ᵣ-x̄ₗ
    N = nₗ+nᵣ
    x̄ = (nₗ*x̄ₗ+nᵣ*x̄ᵣ)/N
    Sₗ = σₗ²*(nₗ-1)
    Sᵣ = σᵣ²*(nᵣ-1)
    S = Sₗ+Sᵣ+nₗ*nᵣ/N*δ^2
    return (x̄, S/(N-1), N)
end
# Evaluation API #
# ========================== #
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram,
                      decision::AbstractVector;
                      solver = JuMP.UnsetSolver())

Evaluate the first-stage `decision` in `stochasticprogram`.

In other words, evaluate the first-stage objective at `decision` and solve outcome models of `decision` for every available scenario. Optionally, supply a capable `solver` to solve the outcome models. Otherwise, any previously set solver will be used.
"""
function evaluate_decision(stochasticprogram::StochasticProgram{2}, decision::AbstractVector; solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot evaluate decision without a solver.")
    end
    return _eval(stochasticprogram, decision, internal_solver(supplied_solver))
end
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram,
                      decision::AbstractVector,
                      scenario::AbstractScenario;
                      solver = JuMP.UnsetSolver())

Evaluate the result of taking the first-stage `decision` if `scenario` is the actual outcome in `stochasticprogram`.
"""
function evaluate_decision(stochasticprogram::StochasticProgram{2},
                           decision::AbstractVector,
                           scenario::AbstractScenario;
                           solver::SPSolverType = JuMP.UnsetSolver())
    # Use cached solver if available
    supplied_solver = pick_solver(stochasticprogram, solver)
    # Abort if no solver was given
    if isa(supplied_solver, JuMP.UnsetSolver)
        error("Cannot evaluate decision without a solver.")
    end
    outcome = outcome_model(stochasticprogram, decision, scenario; solver = solver)
    status = solve(outcome)
    if status == :Optimal
        return _eval_first_stage(stochasticprogram, decision) + getobjectivevalue(outcome)
    end
    error("Outcome model could not be solved, returned status: $status")
end
"""
    evaluate_decision(stochasticmodel::StochasticModel{2},
                      decision::AbstractVector,
                      sampler::AbstractSampler;
                      solver = JuMP.UnsetSolver(),
                      confidence = 0.9,
                      n = 1000,
                      N = 100,
                      M = 10)

Return a statistical estimate of the objective of the two-stage `stochasticmodel` at `decision` in the form of a confidence interval at level `confidence`, when the underlying scenario distribution is inferred by `sampler`.

In other words, evaluate `decision` on an SAA model of size `n`. Generate an upper bound using the sample variance of the evaluation. The lower bound is calculated from the lower bound around the true optimum generated by `confidence_interval`

See also: [`confidence_interval`](@ref)
"""
function evaluate_decision(stochasticmodel::StochasticModel{2},
                           decision::AbstractVector,
                           sampler::AbstractSampler{S};
                           solver::SPSolverType = JuMP.UnsetSolver(),
                           confidence::AbstractFloat = 0.9,
                           n::Integer = 1000,
                           N::Integer = 100,
                           M::Integer = 10) where S <: AbstractScenario
    eval_model = SAA(stochasticmodel, sampler, n)
    # Condidence level
    α = (1-confidence)/2
    # Compute confidence interval around true optimum for lower bound
    CI = confidence_interval(stochasticmodel, sampler; solver = solver, confidence = 1-α, N = N, M = M)
    L = lower(CI)
    # Upper bound
    cᵀx = _eval_first_stage(eval_model, decision)
    𝔼Q, σ = _stat_eval_second_stages(eval_model, decision, internal_solver(solver))
    U = cᵀx + 𝔼Q + quantile(Normal(0,1), 1-α)*σ
    return ConfidenceInterval(L, U, confidence)
end
"""
    confidence_interval(stochasticmodel::StochasticModel{2},
                        sampler::AbstractSampler;
                        solver = JuMP.UnsetSolver(),
                        confidence = 0.95,
                        N = 100,
                        M = 10)

Generate a confidence interval around the true optimum of the two-stage `stochasticmodel` at level `confidence`, when the underlying scenario distribution is inferred by `sampler`.

`N` is the size of the SAA models used to generate the interval and generally governs how tight it is. `M` is the amount of SAA samples used.
"""
function confidence_interval(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler{S}; solver::SPSolverType = JuMP.UnsetSolver(), confidence::AbstractFloat = 0.95, N::Integer = 100, M::Integer = 10) where {S <: AbstractScenario}
    # Condidence level
    α = 1-confidence
    # Lower bound
    Qs = Vector{Float64}(undef, M)
    for i = 1:M
        saa = SAA(stochasticmodel, sampler, N)
        Qs[i] = VRP(saa, solver = solver)
    end
    Q̂ = mean(Qs)
    σ = std(Qs)
    t = quantile(TDist(M-1), 1-α)
    L = Q̂ - t*σ
    U = Q̂ + t*σ
    return ConfidenceInterval(L, U, 1-α)
end
# ========================== #


text:
module pex

include("utils.jl")

end # module


text:
export PackedUniform

mutable struct PackedUniform <: PackedSamplableBelief
  a::Float64
  b::Float64
  PackedSamplableTypeJSON::String
end

function convert(::Union{Type{<:PackedSamplableBelief},Type{<:PackedUniform}},
                  obj::Distributions.Uniform)
  #
  packed = PackedUniform(obj.a, obj.b,
                        "IncrementalInference.PackedUniform")
  #
  return JSON2.write(packed)
end

convert(::Type{<:SamplableBelief}, obj::PackedUniform) = return Uniform(obj.a, obj.b)

# TODO stop-gap string storage of Distrubtion types, should be upgraded to more efficient storage
function normalfromstring(str::AbstractString)
  meanstr = match(r"μ=[+-]?([0-9]*[.])?[0-9]+", str).match
  mean = split(meanstr, '=')[2]
  sigmastr = match(r"σ=[+-]?([0-9]*[.])?[0-9]+", str).match
  sigma = split(sigmastr, '=')[2]
  Normal{Float64}(parse(Float64,mean), parse(Float64,sigma))
end

function mvnormalfromstring(str::AbstractString)
  means = split(split(split(str, 'μ')[2],']')[1],'[')[end]
  mean = Float64[]
  for ms in split(means, ',')
    push!(mean, parse(Float64, ms))
  end
  sigs = split(split(split(str, 'Σ')[2],']')[1],'[')[end]
  sig = Float64[]
  for ms in split(sigs, ';')
    for m in split(ms, ' ')
      length(m) > 0 ? push!(sig, parse(Float64, m)) : nothing
    end
  end
  len = length(mean)
  sigm = reshape(sig, len,len)
  MvNormal(mean, sigm)
end

function categoricalfromstring(str::AbstractString)
  # pstr = match(r"p=\[", str).match
  psubs = split(str, '=')[end]
  psubs = split(psubs, '[')[end]
  psubsub = split(psubs, ']')[1]
  pw = split(psubsub, ',')
  p = parse.(Float64, pw)
  return Categorical(p ./ sum(p))
end

# NOTE SEE EXAMPLE IN src/Flux/FluxModelsSerialization.jl
function _extractDistributionJson(jsonstr::AbstractString, checkJson::AbstractVector{<:AbstractString})
  # Assume first word after split is the type
  mjs = findfirst(r"([a-zA-Z0-9._]+)\w", checkJson[2])
  maybemodule = split(checkJson[2][mjs], '.')
  # get dedicated Module or revert to Main
  packmodule = 1 < length(maybemodule) ? getfield(Main, Symbol(maybemodule[1])) : Main
  packtype = getfield(packmodule, Symbol(maybemodule[end]))
  packed = JSON2.read(jsonstr, packtype)
  # call the dedicated converter for this packed type using dispatch
  convert(SamplableBelief, packed)
end

# NOTE part of new effort to overhaul the SamplableBelief serialization approach
# maybe it all becomes a JSON struct sort of thing in the long run.
StringThemSamplableBeliefs = Union{Normal, MvNormal, Categorical, DiscreteNonParametric, BallTreeDensity, AliasingScalarSampler}
convert(::Type{<:PackedSamplableBelief}, obj::StringThemSamplableBeliefs) = string(obj)

function convert(::Type{<:SamplableBelief}, str::Union{<:PackedSamplableBelief,<:AbstractString})
  # TODO improve use of multidispatch and packing of Distribution types
  # extractdistribution(str::AS) where {AS <: AbstractString}
  # TODO improve use of multidispatch and packing of Distribution types
  # TODO use startswith
  checkJson = split(str, r"PackedSamplableTypeJSON")
  if str == ""
    return nothing
  elseif length(checkJson) == 2
    # TODO this is the new direction for serializing (pack/unpack) of <:Samplable objects
    # NOTE uses intermediate consolidation keyword search pattern `SamplableTypeJSON`
    return _extractDistributionJson(str, checkJson)
  elseif startswith(str, "DiagNormal")
    # Diags are internally squared, so only option here is to sqrt on input.
    return mvnormalfromstring(str)
  elseif occursin(r"FullNormal", str)
    return mvnormalfromstring(str)
  elseif (occursin(r"Normal", str) )# && !occursin(r"FullNormal", str))
    return normalfromstring(str)
  elseif occursin(r"Categorical", str)
    return categoricalfromstring(str)
  elseif occursin(r"DiscreteNonParametric", str)
    return categoricalfromstring(str)
  elseif occursin(r"KDE:", str)
    return convert(BallTreeDensity, str)
  elseif occursin(r"AliasingScalarSampler", str)
    return convert(AliasingScalarSampler, str)
  else
    error("Don't know how to extract distribution from str=$(str)")
  end
end

function convert(::Type{<:PackedSamplableBelief}, obj::SamplableBelief)
  # FIXME must use string, because unpacking templated e.g. PackedType{T} has problems, see DFG #668
  string(obj)
end

#


text:
"""
    knockout(gene_ids::Vector{String})

A modification that zeroes the bounds of all reactions that would be knocked
out by the specified genes (effectively disables the reactions).
"""
knockout(gene_ids::Vector{String}) =
    (model, optmodel) -> _do_knockout(model, optmodel, gene_ids)

"""
    knockout(gene_id::String)

A helper variant of [`knockout`](@ref) for a single gene.
"""
knockout(gene_id::String) = knockout([gene_id])

"""
    _do_knockout(model::MetabolicModel, opt_model)

Internal helper for knockouts on generic MetabolicModels. This can be
overloaded so that the knockouts may work differently (more efficiently) with
other models.
"""
function _do_knockout(model::MetabolicModel, opt_model, gene_ids::Vector{String})
    for (rxn_num, rxn_id) in enumerate(reactions(model))
        rga = reaction_gene_association(model, rxn_id)
        if !isnothing(rga) &&
           all([any(in.(gene_ids, Ref(conjunction))) for conjunction in rga])
            set_bound(rxn_num, opt_model, ub = 0, lb = 0)
        end
    end
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.Variables (ModiaMath/Variables/_module.jl)
#

# Define SimulationModel struct and functions

"""
    sm = ModiaMath.SimulationModel(model::ModiaMath.AbstractComponentWithVariables;
                                   startTime = 0.0, stopTime = 1.0, tolerance = 1e-4,
                                   interval = (stopTime-startTime)/500.0)

Return a simulationModel `sm` struct that can be simulated with [`ModiaMath.simulate!`](@ref).
The given `model` is typically constructed with the [`@component`](@ref) macro.
As keyword arguments, default values for `startTime`, `stopTime`, `tolerance`, `interval` can be given.
"""
mutable struct SimulationModel <: ModiaMath.AbstractSimulationModel
    modelName::String
    simulationState::ModiaMath.SimulationState
    var::ModelVariables
    model::ModiaMath.AbstractComponentWithVariables

    function SimulationModel(model::ModiaMath.AbstractComponentWithVariables;
                             startTime = 0.0,
                             stopTime  = 1.0,
                             tolerance = 1e-4,
                             interval  = (stopTime-startTime)/500.0,
                             structureOfDAE = ModiaMath.DAE_LinearDerivativesAndConstraints,
							 has_constraintDerivatives = false,
                             hev = 1e-8,
                             scaleConstraintsAtEvents::Bool = true)
        modelName = ModiaMath.componentName(model)
        var = ModelVariables(model)
        x = zeros(var.nx)
        x_fixed   = fill(false,var.nx)
        x_nominal = fill(1.0, var.nx)
        ModiaMath.copy_start_to_x!(var,x,x_fixed,x_nominal)

        # Last nfc equations are the constraint equations
        is_constraint = fill(false, var.nx)
        for i = (var.nx-var.nfc+1):var.nx
            is_constraint[i] = true
        end

        simulationState = ModiaMath.SimulationState(modelName, getModelResidues!, x, getVariableName;
                                                    structureOfDAE = structureOfDAE,
                                                    is_constraint = is_constraint,
                                                    has_constraintDerivatives = has_constraintDerivatives,
                                                    x_fixed = x_fixed, x_nominal=x_nominal,
                                                    getResultNames = getResultNames, storeResult! = storeVariables!,
                                                    getResult        = getResult,
                                                    defaultStartTime = startTime,
                                                    defaultStopTime  = stopTime,
                                                    defaultTolerance = tolerance,
                                                    defaultInterval  = interval,
                                                    hev = hev,
                                                    scaleConstraintsAtEvents = scaleConstraintsAtEvents)

        # model._internal.simulationState = simulationState
        new(String(modelName), simulationState, var, model)
    end
end

print_ModelVariables(simulationModel::SimulationModel) = print_ModelVariables(simulationModel.var)

getVariableName(model,vcat,vindex) = ModiaMath.DAE.getVariableName(model,vcat,vindex;
                                                                   xNames = model.var.x_names)

getResultNames(model::Any) = model.var.result_names

function storeVariables!(res::ModiaMath.RawResult, model, t::Float64, x::Vector{Float64}, derx::Vector{Float64}, w::Vector{Float64})
    res.nt += 1
    i = res.nt
    ntMax = size(res.data,1)
    if i > ntMax
        # Result storage is full, double the current storage
        res.data = [res.data; zeros(ntMax,size(res.data,2))]
    end

    # Copy data in raw result storage
    for v in model.var.result_var
        if ModiaMath.isScalar(v)
            res.data[i, v.iresult] = v.value
        else
            vv = v.value
            for j in 1:length(vv)
                res.data[i, v.iresult+j-1 ] = vv[j]
            end
        end
    end
    return nothing
end

function getResult(model::ModiaMath.AbstractSimulationModel, res::ModiaMath.RawResult)
    nt         = res.nt
    data       = res.data

    seriesDict = Dict{AbstractString,Any}()
    varDict    = Dict{Symbol,ModiaMath.AbstractVariable}()
    for v in model.var.result_var
        name    = instanceName(v)
        len     = length(v.value)
        iresult = v.iresult
        if len == 1
            seriesDict[string(name)] = @view data[1:nt, iresult]
        else
            seriesDict[string(name)] = @view data[1:nt, iresult:iresult+len-1]
        end
        varDict[name] = v
   end
   result = ModiaMath.ResultWithVariables(seriesDict, varDict, model.modelName)
   return result, nt
end

computeVariables!(model::ModiaMath.AbstractComponentWithVariables,
                  simulationState::ModiaMath.SimulationState) =
          error("Function ModiaMath.computeVariables!(model) not defined for model ", typeof(model))

function getModelResidues!(m::SimulationModel, t::Float64, _x::Vector{Float64}, _derx::Vector{Float64}, _r::Vector{Float64}, _w::Vector{Float64})
    # Copy _x and _derx values to variables
    ModiaMath.copy_x_and_derx_to_variables!(t, _x, _derx, m.var)

    # Compute model variables
    m.simulationState.time = t
    computeVariables!(m.model, m.simulationState)

    # Copy variables to residues
    ModiaMath.copy_variables_to_residue!(m.var,_x,_derx,_r)
    return nothing
end

# Return a table of actual variable and residue values from nonlinear solver in case of error
function ModiaMath.getVariableAndResidueValues(m::SimulationModel)
   var = m.var
   v_table = get_variableValueTable(var)
   r_table = get_residueValueTable(var, m.simulationState.derxev)
   return (v_table, r_table)
end


text:
# HTTP request
#
# - method   => valid HTTP method string (e.g. "GET")
# - resource => requested resource (e.g. "/hello/world")
# - headers  => HTTP headers
# - data     => request data
# - state    => used to store various data during request processing
# - sock      => The socket
#

export Request,
       Reader,
       save,
       loadedcl,
       leftcl,
       StringReader,
       event,
       event!,
       CloseRequestException,
       init,
       readinfo,
       readheaders,
       parsedata,
       param,
       body,
       bodes

using JSON

type UploadFile
    name::STR_TYPE
    content_type::STR_TYPE
    size::Integer
    io::IO
    path::STR_TYPE
    headers::Associative{STR_TYPE, STR_TYPE}
end

UploadFile(n::STR_TYPE, c::STR_TYPE, s::Integer, io::IO; headers=Dict{STR_TYPE,STR_TYPE}()) = UploadFile(n, c, s, io, "", headers)

type RequestException <: Exception
    info
end

RequestException(key::Union(STR_TYPE, (STR_TYPE...)), message::STR_TYPE; kwargs...) = begin
    all = append!(Any[(:key, key), (:message, message)], kwargs)
    RequestException(Dict{Any, Any}(all))
end

type Request
    io::AbstractIOSocket
    protocol::Union(N_TYPE,Protocol)
    method::Union(N_TYPE,STR_TYPE)
    resource::Union(N_TYPE,STR_TYPE)
    headers::Headers
    data::Union(N_TYPE,STR_TYPE)
    env::Associative{Union(STR_TYPE,Symbol), Any}
    finished::Bool
    query::Union(N_TYPE,STR_TYPE)
    get::Associative{STR_TYPE, Any}
    post::Associative{STR_TYPE, Any}
    events::Events
    read_bytes::Integer
    header_size::Integer
    l::Associative{Any, Any}
    websocket
    content_length::Integer
    chunk::Bool
    files::Union(N_TYPE,Array{UploadFile})
    content_type::Union(N_TYPE, STR_TYPE)
end

Request(io::AbstractIOSocket) = Request(
    io, N, N, N, headers(), N, Dict{Union(STR_TYPE,Symbol),Any}(), false,
    N, params(), params(), Events(),
    0, 0, Dict{Any,Any}(), N, 0, false, N, N
)

"""
Return number of loaded bytes of content length
"""
loadedcl(r::Request) = r.read_bytes - r.header_size

"""
Return left loaded bytes of content length
"""
leftcl(r::Request) = r.content_length - loadedcl(r)

trigger(r::Request, e::Event) = trigger(r.events, e, r)
listener(r::Request, key::Union(Symbol,STR_TYPE), cb::Function) = listener(r.events, key, cb)

show(io::IO, r::Request) = print(
    io,
    "Request(",
    r.method,
    " ",
    r.resource,
    ", ",
    length(r.headers),
    " Headers, ",
    length(r.env),
    " Enviroment, ",
    isa(r.data, STR_TYPE) ? "$(sizeof(r.data)) Bytes in Body" : r.data,
    ", ",
    r.query,
    " Query STR_TYPE,",
    r.read_bytes,
    " Bytes read,",
    r.header_bytes,
    " Header size",
    ")"
)

soma(a::Integer, b::Integer) = (a + b; true)

function soma(a::Integer, b::Integer)
    a + b
    true
end

header(r::Request, key::STR_TYPE, value::STR_TYPE) = header(r.headers, key, value)

Base.read{T}(r::Request, datatype::Type{T}, size::Integer) = begin
    data = Base.read(r.io.io, datatype, size)
    r.read_bytes += size
    data
end

Base.read(r::Request, size::Integer) = Base.read(r, BYTE_TYPE, size)

Base.readuntil(r::Request, delim) = begin
    data = Base.readuntil(r.io.io, delim)
    r.read_bytes += sizeof(data)
    data
end

Base.readline(r::Request) = readline(r.io)

using StreamReader

function Reader{T,D}(r::Request, bp::PartsIterator; read_type::Type{T} = BYTE_TYPE,
    datatype::Type{D} = N_TYPE, prestart::Union(N_TYPE,Function) = N,
    kwargs...)
    ReaderIterator(bp; kwargs...) do ri
        ri.read = (size) -> Base.read(r, read_type, size)

        if datatype != N_TYPE
            setdatatype(ri, datatype)
        end

        if prestart != N
            prestart(ri)
        end
    end
end

Reader(r::Request, len::Integer, bsize::Integer=0; kwargs...) = Reader(r,
    PartsIterator(len, bsize); kwargs...)

function save(r::Request, out::IO, len::Integer, bsize::Integer=0;
    cb::Function = emptyfn, kwargs...)
    r = Reader(r, len, bsize; kwargs...)
    for data in
        cb(r)
        Base.write(out, data)
    end
end

save(r::Request, out::IO, bsize::Integer=0; kwargs...) = save(r, out, r.content_length, bsize; kwargs...)

>>(i::ReaderIterator, o::IO) = begin
    for data in i
        Base.write(o, data)
    end
end

function StringReader{D<:STR_TYPE}(r::Request, pi::PartsIterator;
    datatype::Type{D} = UTF8String, kwargs...)
    Reader(io, pi; datatype=datatype, kwargs...)
end

StringReader(r::Request, len::Integer, bsize::Integer=0; kwargs...) = StringReader(io,
    StreamReader.PartsIterator(len, bsize); kwargs...)

body(req::Request) = read(req, req.content_length)
bodys(req::Request) =  bytestring(body(req))

function readfirstline(io::AbstractIOSocket)
    line = readlinebarestr(io)
    s = sizeof(line) + 2

    if line == ""
        throw(RequestException("invalid_header", "Empty first line."; line_size=s, line=line))
    end

    s, line
end

function readheaders(io::AbstractIOSocket, lnum::Integer=0; h::Headers=headers())
    s = Integer[0]

    lnum += eachline(io, :s!) do itr, line
        s[1] += sizeof(line) + 2 # + CRLF

        if isempty(line)
            # end headers
            stop()
        end

        parts = split(line, ": ", 2)

        if length(parts) != 2
            throw(RequestException(("invalid_header", "invalid_format"), "Header " *
                repr(line) * ", in line " * string(linum+iter.count) *
                ", not contains a valid format."; ready_bytes=s, line_number=(lnum + itr.count), line=line, line_size=ls))
        end

        key, value = parts
        header(h, key, value)
    end

    s[1], lnum, h
end

function readheaders(r::Request)
    s, lines, r.headers = readheaders(r.io)
    r.read_bytes += s
    r.header_size = s

    cl = get(r.headers, r.protocol.header_msg_size, N)

    if cl != N
        r.content_length = parseint(cl)
    end

    ct = get(r.headers, r.protocol.header_msg_type, N)

    if ct != N
        r.content_type = ct
    end

    r
end

type CloseRequestException <: Exception
    info
end

"""Parse HTTP info on first line:

    METHOD RESOURCE PROTOCOL/PROTOCOL_VERSION

Example:

    GET /index.html
    GET /index.html HTTP/1.1
    GET /directory%20name/file.txt HTTP/1.1
"""
function parseinfo(line::STR_TYPE, default_protocol::Union(N_TYPE, (STR_TYPE, STR_TYPE))=N)
    parts = split(line, " ")
    method = resource = protocol = N

    if length(parts) < 2 || length(parts) > 3
        throw(RequestException(("invalid_header_info", "invalid_format"), "Header " *
            repr(line) * ", not contains a valid format."; line=line))
    elseif length(parts) == 2
        if default_protocol == N
            throw(RequestException(("invalid_header_info", "protocol_not_informed"),
                "Protocol isn't informed."; line=line))
        else
            protocol = default_protocol
            method, resource = parts
        end
    else
        method, resource, protocol = parts
        parts = split(protocol, "/")

        if length(parts) != 2
            throw(RequestException(("invalid_header_info", "invalid_protocol_format"), "Invalid PROTOCOL format " *
                repr(protocol) * "."; protocol=protocol))
        end

        protocol = parts
    end

    method, resource, protocol
end

function readinfo(io::AbstractIOSocket, protocols::Protocols,
    default_protocol::Union(N_TYPE, (STR_TYPE, STR_TYPE))=N)
    line = readlinebarestr(io)
    s = sizeof(line) + 2

    if line == "quit"
        throw(CloseRequestException("close command recieved."))
    end

    method, resource, protoc = parseinfo(line, default_protocol)
    query = ""

    if resource[1] == '/' # is path
        # get query string
        query_pos = search(resource, '?')

        if query_pos > 0
            query = resource[query_pos+1:end]
            resource = resource[1:query_pos-1]
        end
    end

    if protoc != N
        protoc = protocol(protocols, tuple(protoc...))
    end

    (s, method, resource, protoc, query)
end

function readinfo(r::Request, protocols::Protocols,
    default_protocol::Union(N_TYPE, (STR_TYPE, STR_TYPE))=N)
    (r.read_bytes, r.method, r.resource, r.protocol, r.query) = readinfo(r.io,
        protocols, default_protocol)
    r
end

function init(r::Request, protocols, default_protocol::Union(Protocol, N_TYPE)=N)
    readinfo(r, protocols, default_protocol != N? key(default_protocol) : N)
    readheaders(r)

    if r.query != N
        r.get = parseqsr(r.query)
    end

    r.protocol.req_post_init(r)

    res = Response(r.io, r.protocol)
    res.protocol.res_post_init(res)

    trigger(r, Event("post_init"; args=(r, res)))

    r, res
end

"""
Parse content data if content-type is "application/x-www-form-urlencoded",
"application/json" or "multipart/form-data". If parsed, return ``true``,
otherwise, ``false``.
"""
function parsedata(r::Request)
    if r.content_type != N
        ct = r.content_type
        if search(ct, "application/x-www-form-urlencoded").start > 0
            r.post = parseqsr(bodys(r))
        elseif search(ct, "application/json").start > 0
            r.post = JSON.parse(bodys(r))
        elseif search(ct, "multipart/form-data").start > 0
            parsempfd(r)
        end
        return true
    end
    false
end

"""
Return a request param. If ``key`` exists in ``Request.get``, return it, otherwise
``key`` exists in ``Request.post``, return it, otherwise, return ``default``.
For default, value of ``default`` is ``N``.
"""
param(r::Request, key::STR_TYPE, default=N) = (haskey(r.get, key) ?
    r.get[key] : (haskey(r.post, key) ? r.post[key] : default))

"""
Extract multipart boundary of content type.
Returns tuple of ``(boundary, eof_boundary)``.

Example:

.. doctest::

    julia> multipartboundary("multipart/form-data; boundary=AxdFtg1")
    ("AxdFtg1", "AxdFtg1--")

"""
function multipartboundary(content_type::STR_TYPE)
    d = content_type
    b_pos = search(d, "boundary=")
    e_pos = search(d, "; ", b_pos.stop)
    pose = e_pos.start > 0 ? e_pos.start - 1: sizeof(d)
    boundary = d[b_pos.stop+1:pose]
    boundary = string("--", boundary)
    boundary, string(boundary, "--")
end

function parsecontentdisposition(s::STR_TYPE)
    d = Dict{STR_TYPE, STR_TYPE}
    disposition, s = split(s, ";", 2)

    name_p = search(s, " name=\"")
    name_end = search(s, "; ", name_p.stop)

    if name_end.start == 0
        name = s[name_p.stop+1:end-1]
    else
        name = s[name_p.stop+1:name_end.start-2]
    end

    filename = ""

    fname_p = search(s, " filename=\"")

    if fname_p.start > 0
        fname_end = search(s, "; ", fname_p.stop)
        if fname_end.start == 0
            filename = s[fname_p.stop+1:end-1]
        else
            filename = s[fname_p.stop+1:fname_end.start-2]
        end
    end

    disposition, name, filename
end

typealias BoundaryTuple (Array{BYTE_TYPE}, Array{BYTE_TYPE})

function boundarywriter(io::AbstractIOSocket, boundaryd::BoundaryTuple)
    oio = IOBuffer()
    prev_line = BYTE_TYPE[]
    stop = false
    s = 0

    while true
        next_line = readline(io)
        s += sizeof(next_line)

        if next_line == boundaryd[1]
            # write previous line without CRLF
            write(oio, prev_line[1:end-2])
            # stop
            break
        elseif next_line == boundaryd[2]
            # write previous line without CRLF
            write(oio, prev_line[1:end-2])
            stop = true
            # stop
            break
        end

        write(oio, prev_line)
        prev_line = next_line
    end

    seekstart(oio)
    s, stop, oio
end

function parsempfdpart(io::AbstractIOSocket, boundaryd::BoundaryTuple, encoding::STR_TYPE, d, files)
   s, lnum, h = readheaders(io; h=Dict{STR_TYPE, STR_TYPE}())
   stop = true

   if length(h) > 0
        disposition, name, filename = parsecontentdisposition(h["Content-Disposition"])

        if !haskey(d, name)
            values = d[name] = {}
        else
            values = d[name]
        end

        if haskey(h, "Content-Type")
            ls, stop, mio = boundarywriter(io, boundaryd)
            s += ls
            value = UploadFile(filename, h["Content-Type"], sizeof(mio.data), mio; headers=h)
            push!(files, value)
            push!(values, value)
        else
            ls, stop, mio = boundarywriter(io, boundaryd)
            s += ls
            value = takebuf_array(mio)

            if !isempty(encoding)
                value = decode(value, encoding)
            end

            push!(values, value)
        end
    end

    s, stop
end

function parsempfd(io::AbstractIOSocket, boundary::(STR_TYPE, STR_TYPE), encoding::STR_TYPE="")
    # read boundary on first content line
    line = readlinebare(io)
    s = sizeof(line) + 2
    line == boundary[1].data || throw(RequestException("parse_multipart",
        "Expected boundary data.", (s, [:boundary => boundary[1], :line => line])))

    d = params()
    files = UploadFile[]
    boundaryd = tuple(Array{BYTE_TYPE}[(_ * CRLF).data for _ in boundary]...)
    stop = false

    while !stop
        ls, stop = parsempfdpart(io, boundaryd, encoding, d, files)
        s += ls
    end

    s, d, files
end

"""
Parse multipart form data content into ``Request.post`` and ``Request.files``.
"""
function parsempfd!(r::Request)
    # the temporary
    io = IOBuffer()
    # read all content into BUFFER
    write(io, body(r))
    seekstart(io)

    parse = (io) -> begin
        boundary = S.multipartboundary(r.headers["Content-Type"])
        encoding = ""
        s, r.post, r.files = S.parsempfd(IOSocket(io), boundary, encoding)
        # discard read size...
        r.post, r.files # dict of parameters, files list
    end

    parse(io)
end


text:
Base.exit_on_sigint(false)

try
    include("Jus.jl")

    Jus.exec(Jus.serve, ARGS)
catch err
    !(err isa InterruptException) && @warn "Error" exception=(err, catch_backtrace())
    exit(1)
end


text:
using Altro
using TrajectoryOptimization
using Test
using StaticArrays
using ForwardDiff
using RobotDynamics
using LinearAlgebra
const RD = RobotDynamics
const TO = TrajectoryOptimization

T = Float64
n, m, N = 3, 2, 11
tf = 2.0
xf = [1,2,3.]

con_eq = TO.GoalConstraint(xf, SA[1,3])
inds_eq = 1:10
confun_eq(x) = (x - xf)[[1,3]]

xc, yc, rad = [4,2,2,0.], [0,1,2,3.], [0.2, 0.1, 0.3, 0.4]
con_in = TO.CircleConstraint(n, xc, yc, rad)
inds_in = 2:6
confun_in(x) = rad.^2 .- (x[1] .- xc).^2 .- (x[2] .- yc).^2

con_so = TO.NormConstraint(n, m, 3, TO.SecondOrderCone(), :state)
inds_so = N:N
confun_so(x) = [x; 3]

cons = ConstraintList(n, m, N)
add_constraint!(cons, con_eq, inds_eq)
add_constraint!(cons, con_in, inds_in, sig=RD.InPlace())
add_constraint!(cons, con_so, inds_so)

Z = SampledTrajectory{n,m}(randn(n,N), randn(m,N), dt=0.1)
opts = SolverOptions()
alcosts = zeros(N)
conset = Altro.ALConstraintSet{T}()
Altro.initialize!(conset, cons, Z, opts, alcosts)

@test length(conset) == 3
@test length(conset.constraints) == 3
@test length(conset.c_max) == N
@test length(conset.μ_max) == 3
@test conset[begin].con == con_eq
@test conset[2].con == con_in
@test conset[end].con == con_so
@test conset[1].sig == RD.StaticReturn()
@test conset[2].sig == RD.InPlace()

# Generate a random trajectory
X = [randn(T,n) for k = 1:N]
U = [randn(T,m) for k = 1:N-1]
Z = SampledTrajectory{n,m}(X, U, tf=2.0)

# Methods
Altro.evaluate_constraints!(conset, Z)
Altro.constraint_jacobians!(conset, Z)
Altro.alcost(conset)
@test sum(alcosts) != 0
Altro.algrad!(conset)
Altro.alhess!(conset)
Altro.max_violation!(conset[1])
Altro.dualupdate!(conset)
Altro.reset_penalties!(conset)
Altro.penaltyupdate!(conset)
@test Altro.max_penalty(conset) == 10.0

## Reset
opts.penalty_initial = 0.1
Altro.reset!(conset)
@test Altro.max_penalty(conset) ≈ 0.1


text:
using PortAudio

const pstream = PortAudioStream(0, 2)
sendmorsesound(t, f) = write(pstream, SinSource(eltype(stream), samplerate(stream)*0.8, [f]), (t/1000)s)

char2morse = Dict[
      "!" => "---.", "\"" => ".-..-.", "$" => "...-..-", "'" => ".----.",
      "(" => "-.--.", ")" => "-.--.-", "+" => ".-.-.", "," => "--..--",
      "-" => "-....-", "." => ".-.-.-", "/" => "-..-.", "0" => "-----",
      "1" => ".----", "2" => "..---", "3" => "...--", "4" => "....-", "5" => ".....",
      "6" => "-....", "7" => "--...", "8" => "---..", "9" => "----.", ":" => "---...",
      ";" => "-.-.-.", "=" => "-...-", "?" => "..--..", "@" => ".--.-.", "A" => ".-",
      "B" => "-...", "C" => "-.-.", "D" => "-..", "E" => ".", "F" => "..-.",
      "G" => "--.", "H" => "....", "I" => "..", "J" => ".---", "K" => "-.-",
      "L" => ".-..", "M" => "--", "N" => "-.", "O" => "---", "P" => ".--.",
      "Q" => "--.-", "R" => ".-.", "S" => "...", "T" => "-", "U" => "..-",
      "V" => "...-", "W" => ".--", "X" => "-..-", "Y" => "-.--", "Z" => "--..",
      "[" => "-.--.", "]" => "-.--.-", "_" => "..--.-"]

function sendmorsesound(freq, duration)
cpause() = sleep(0.080)
wpause = sleep(0.400)

dit() = sendmorsesound(0.070, 700)
dash() = sensmorsesound(0.210, 700)
sendmorsechar(c) = for d in char2morse(c) d == '.' ? dit(): dash() end end
sendmorseword(w) = for c in w sendmorsechar(c) cpause() end wpause() end
sendmorse(msg) = for word in uppercase(msg) sendmorseword(word) end

sendmorse("sos sos sos")
sendmorse("The case of letters in Morse coding is ignored."

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
@with_kw struct FrequencyDomainMonitor{TQuantity ,GeometryType} <: AbstractMonitor
	geometry::GeometryType
	quantity::TQuantity
	frequency::Float
end