text:
using PrettyTables,Plots,Statistics,Parameters

"""
        describe(gl::ModelGL,gl_2::ModelGL)

Print solution status as well as values and description of parameters and aggregate values:
Comparison of initial and terminal steady-state.
"""
function describe(gl::ModelGL,gl_2::ModelGL)
    print_status(gl,gl_2 )
    print_params(gl,gl_2 )
end

"""
        describe(gl::ModelGL)

Print solution status as well as values and description of parameters and aggregate values.
"""
function describe(gl::ModelGL)
    print_status(gl )
    print_params(gl )
end

"""
        print_params(gl::ModelGL)

Print values and description of parameters and aggregate values.
"""
function print_params(gl::ModelGL)

        params = ["β","ν", "ϕ","γ","η","fin","sep"]
        desc   = ["discount","UI benefits", "borr. constraint",
                        "risk aversion", "eta", "job-finding rate",
                        "job-separation rate", ]
        vals   = [ getfield(gl,Symbol(i) )   for i in params ]

        data = hcat((params,vals,desc)... )

        tgt_nam  = ["NE","νY", "B_4Y","D1_4Y","D2_4Y"]
        tgt_dsc  = ["avg. hours of employed","UI benefit / GDP",
                        "liquid wealth / GDP", "HH debt-to-GDP (initial ss)",
                        "HH debt-to-GDP (terminal ss)" ]
        tgt_vals = [ getfield(gl,Symbol(i) )   for i in tgt_nam ]

        if gl.steadystate_type == "initial"
                actual_nam  = ["NE_actual","νY_actual", "B_4Y_actual","D_4Y_actual" ]
                actual_vals = [[ getfield(gl,Symbol(i) )   for i in actual_nam ]...; NaN ]
        elseif  gl.steadystate_type == "terminal"
                actual_nam  = ["NE_actual","νY_actual", "B_4Y_actual", "D_4Y_actual"]
                actual_vals = [[ getfield(gl,Symbol(i) )   for i in actual_nam[1:3] ]...; NaN; getfield(gl,Symbol(actual_nam[end])) ]
        end

        params = ["β","ν", "ϕ","γ","η","fin","sep","ρ","σϵ"]
        desc   = ["Discount factor","Unemployment benefits", "Borrowing constraint",
                        "Risk aversion", "Leisure curvature", "Job-finding rate",
                        "Job-separation rate","Persistance income","SD income shocks"  ]
        vals   = [ getfield(gl,Symbol(i) )   for i in params ]
        data = hcat((params,vals ,fill("",length(params)),desc)... )

        tgt_nam  = ["NE","νY", "B_4Y","D1_4Y","D2_4Y" ]
        tgt_dsc  = ["Average hours of employed","Unemployment benefits per GDP",
                "Liquid wealth per GDP", "Household debt-to-GDP (initial ss)",
                "Household debt-to-GDP (terminal ss)"  ]
        tgt_vals =  [ getfield(gl,Symbol(i) )   for i in tgt_nam  ]
        actual_nam      = ["NE_actual","νY_actual", "B_4Y_actual","D_4Y_actual" ]

        if (gl.steadystate_type == "initial")
                actual_vals  = [[ getfield(gl,Symbol(i) )   for i in actual_nam ]...; NaN;   ]
        elseif (gl.steadystate_type == "terminal")
                actual_vals  =  [[ getfield(gl ,Symbol(i) )   for i in actual_nam[1:3]]...; NaN; getfield(gl ,Symbol(actual_nam[end]))  ]
        end

        data2 = hcat((tgt_nam,actual_vals ,tgt_vals,tgt_dsc)... )

        filler1 = ["Parameters" ""  "" "" ]
        filler2 = [" " ""   "" ""; "Aggregates" ""   "" ""]
        pretty_table(vcat((filler1,data,filler2,data2)...), tf = tf_simple,
        formatters = ft_printf("%5.3f") ,alignment=[:l,:r,:r, :l],
        header = (["Field", "Value",  "Target","Description" ],
                  [" ", gl.steadystate_type , " "," " ] ),
        highlighters = ( hl_cell( [(1,1);(length(params)+3,1)], crayon"bold") ),
        border_crayon = crayon"bold yellow",  header_crayon = crayon"bold green",  subheader_crayon = crayon"green",
        crop = :none)

end

"""
        print_params(gl::ModelGL,gl_2::ModelGL)

Print values and description of parameters and aggregate values:
Comparison of initial and terminal steady-state.
"""
function print_params(gl::ModelGL,gl_2::ModelGL)

        params = ["β","ν", "ϕ","γ","η","fin","sep","ρ","σϵ"]
        desc   = ["Discount factor","Unemployment benefits", "Borrowing constraint",
                "Risk aversion", "Leisure curvature", "Job-finding rate",
                "Job-separation rate","Persistance income","SD income shocks" ]
        vals_1   = [ getfield(gl,Symbol(i) )   for i in params ]
        vals_2   = [ getfield(gl_2,Symbol(i) )   for i in params ]

        data = hcat((params,vals_1,vals_2,fill("",length(params)),desc)... )

        tgt_nam  = ["NE","νY", "B_4Y","D1_4Y","D2_4Y" ]
        tgt_dsc  = ["Average hours of employed","Unemployment benefits per GDP",
                        "Liquid wealth per GDP", "Household debt-to-GDP (initial ss)",
                        "Household debt-to-GDP (terminal ss)"  ]
        tgt_vals =  [ getfield(gl,Symbol(i) )   for i in tgt_nam  ]

        actual_nam      = ["NE_actual","νY_actual", "B_4Y_actual","D_4Y_actual" ]
        actual_vals_1 = [[ getfield(gl,Symbol(i) )   for i in actual_nam ]...; NaN;   ]
        actual_vals_2 = [[ getfield(gl_2,Symbol(i) )   for i in actual_nam[1:3]]...; NaN; getfield(gl_2,Symbol(actual_nam[end]))  ]

        data2 = hcat((tgt_nam,actual_vals_1,actual_vals_2,tgt_vals,tgt_dsc)... )

        filler1 = ["Parameters" "" ""  "" ""]
        filler2 = [" " "" ""  "" ""; "Aggregates" "" ""  "" ""]
        pretty_table(vcat((filler1,data,filler2,data2)...), tf = tf_simple,
        formatters = ft_printf("%5.3f"),alignment=[:l,:r,:r,:r,:l],
        header = (["Field", "Value 1)", "Value 2)","Target","Description" ],
                  [" ", gl.steadystate_type, gl_2.steadystate_type, "(1)", "" ] ),
        highlighters = ( hl_cell( [(1,1);(length(params)+3,1)], crayon"bold") ),

        border_crayon = crayon"yellow",  header_crayon = crayon"bold green",  subheader_crayon = crayon"green",
        crop = :none)

end

"""
        print_status(gl::ModelGL)

Print solution status.
"""
function print_status(gl::ModelGL)

        data = ["Model" ""   "";
                "Steady-state solved    ?" "."^10 "$(gl.steadystate_solved)";
                "Calibrated             ?" "."^10 "$(gl.calibrated)";
                "Steady-state type      ?" "."^10 "$(gl.steadystate_type)";
                "" ""   "" ]

        pretty_table(data, tf = tf_borderless,
                        noheader = true,
                        cell_alignment = Dict( (1,1) => :l, (2,1) => :l, (3,1) => :l, (4,1) => :l ),
                        formatters = ft_printf("%10.1f", 2),
                        highlighters = (hl_cell( [(1,1);(6,1)], crayon"bold"),
                                        hl_col(2, crayon"dark_gray")),
                        body_hlines = [1,6],
                        body_hlines_format = Tuple('─' for _ = 1:4) )

end

"""
        print_status(gl::ModelGL,gl_2::ModelGL)

Print solution status of both models.
"""
function print_status(gl::ModelGL,gl_2::ModelGL)

        data = ["1) Model" ""   "";
                "Steady-state solved    ?" "."^10 "$(gl.steadystate_solved)";
                "Calibrated             ?" "."^10 "$(gl.calibrated)";
                "Steady-state type      ?" "."^10 "$(gl.steadystate_type)";
                "" ""   "" ]
        data2 = ["2) Model" ""   "";
                "Steady-state solved    ?" "."^10 "$(gl_2.steadystate_solved)";
                "Calibrated             ?" "."^10 "$(gl_2.calibrated)";
                "Steady-state type      ?" "."^10 "$(gl_2.steadystate_type)";
                "" ""   "" ]
        data = hcat(data,data2)
        pretty_table(data, tf = tf_borderless,
                        noheader = true,
                        cell_alignment = Dict( (1,1) => :l, (2,1) => :l, (3,1) => :l, (4,1) => :l ,
                                               (1,4) => :l, (2,4) => :l, (3,4) => :l, (4,4) => :l   ),
                        formatters = ft_printf("%10.1f", 2),
                        highlighters = (hl_cell( [(1,1);(1,4)], crayon"bold"),
                                        hl_col(2, crayon"dark_gray")),
                        body_hlines = [1,6],
                        body_hlines_format = Tuple('─' for _ = 1:4) )

end

"""
        (F1,p1,p2) = plots_figure_1(gl::ModelGL)

Plots figure I graphs.
"""
function plots_figure_1(gl::ModelGL)
        plotly(linewidth = 2.)

        @unpack b_grid,pr = gl                  # objects that are identical across both steady states

        # Figure I (from paper)

        Y1 = gl.Y_actual
        p1 = plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.c_pol[2,b_grid .>= -gl.ϕ],label = "θ = 2")
        plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.c_pol[8,b_grid .>= -gl.ϕ],
                linestyle = :dash,label = "θ = 8")
        title!("consumption")
        xaxis!([-gl.ϕ,12.5])

        p2 = plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.n_pol[2,gl.b_grid .>= -gl.ϕ],label =  "θ = 2")
        plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), gl.n_pol[8,gl.b_grid .>= -gl.ϕ],
                linestyle = :dash,label = "θ = 8")
        title!("labor supply")
        xaxis!([-gl.ϕ,12.5])

        F1 = plot(p1, p2, layout = (1, 2), legend = false)
        display(F1)
        return (F1,p1,p2)
end

"""
        (F4,p1,p2) = plots_figure_4(gl::ModelGL,gl_tss::ModelGL)

Plots figure IV graphs.
"""
function plots_figure_4(gl::ModelGL,gl_tss::ModelGL)
        plotly(linewidth = 2.)

        @unpack b_grid,pr = gl                  # objects that are identical across both steady states

        # Figure IV (from paper)
        Y1 = gl.Y_actual
        b_acc1 =  ((pr'*gl.b_pol[:,b_grid     .>= -gl.ϕ])'     .- b_grid[b_grid .>= -gl.ϕ])/(4*Y1)
        b_acc2 =  ((pr'*gl_tss.b_pol[:,b_grid .>= -gl.ϕ])' .- b_grid[b_grid .>= -gl.ϕ])/(4*Y1)
        p1 =plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1), b_acc1)
        plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1), b_acc2,
                linestyle = :dash)
        title!("bond accumulation policy")
        xaxis!([-2,14],-2:2:14)
        yaxis!([-0.4,0.6],-0.4:0.2:0.6)

        bond_distribution1 = dropdims( pr'    *gl.JD[:,b_grid .>= -gl.ϕ]     ,dims=1)
        bond_distribution2 = dropdims( pr'*gl_tss.JD[:,b_grid .>= -gl.ϕ] ,dims=1)
        p2 =plot( b_grid[b_grid .>= -gl.ϕ]/(4*Y1),  bond_distribution1,label = "initial")
        plot!(b_grid[b_grid .>= -gl.ϕ]/(4*Y1),  bond_distribution2,label = "terminal",
                linestyle = :dash)
        title!("bond distribution")
        xaxis!([-2,14],-2:2:14)
        yaxis!([0,0.004])
        xaxis!([-gl.ϕ,12.5])

        F4 = plot(p1, p2, layout = (2, 1), legend = false)
        display(F4)

        return (F4,p1,p2)
end

"""
        (F3,p1,p2,p3,p4)  = plots_figure_4(gl::ModelGL,gl_tss::ModelGL)

Plots figure III graphs.
"""
function plots_figure_3(gl::ModelGL,Tgl::TransGL)
        plotly(linewidth = 2.)
        Y1 = gl.Y_actual

        # Figure III (from paper)
        Tp = 24    # number of periods plotted

        p1 = plot(0:Tp, Tgl.ϕ_t[1:Tp+1]./(4*Tgl.Y_t[1:Tp+1]),ylims = (0.5,1), yticks = 0.5:0.1:10, legend = false, title = "borrowing limit")         # borrowing limit
        p2 = plot(0:Tp, Tgl.D_4Y_t[1:Tp+1], ylims = (0.08,0.2), yticks = 0:0.02:0.2, legend = false, title = "household debt-to-GDP ratio")             # debt2gdp ratio
        p3 = plot(0:Tp, [gl.r;Tgl.r_t[1:Tp]].*400, ylims = (-2,2.5), yticks = -2:0.5:2, legend = false, title = "interest rate")       # annualized interest rate
        p4 = plot(0:Tp, [0, 100*(Tgl.Y_t[1:Tp+1]./Y1.-1)], ylims = (-1.2,0), yticks = -1:0.2:0, legend = false, title = "output")   # output deviation from steady state
        F3 = plot(p1, p2, p3, p4, layout = (2, 2), legend = false)

        display(F3)

        return (F3,p1,p2,p3,p4)
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
abstract type TechsExportCurtail end
struct CanExport <: TechsExportCurtail end
struct CanCurtail <: TechsExportCurtail end
struct CanExportAndCurtail <: TechsExportCurtail end
struct CanNotExportCurtail <: TechsExportCurtail end

struct Techs{T <: TechsExportCurtail}
    type::T  # have to use the T parameter
    # names::Array{String, 1}
    # classes::AbstractArray{String, 1}
end

struct Parameter
    t::Techs{T} where {T}
end

# the following goes in utils.jl
ExportTiers = []
TechsCanCurtail = []

if !isempty(ExportTiers)
    if !isempty(TechsCanCurtail)
        techs = Techs(CanExportAndCurtail())
    else
        techs = Techs(CanExport())
    end
else
    if !isempty(TechsCanCurtail)
        techs = Techs(CanCurtail())
    else
        techs = Techs(CanNotExportCurtail())
    end
end

p = Parameter(techs)

# the following will be in reopt.jl
add_constraint(p::Parameter) = add_constraint(p.t)

function add_constraint(::Techs{CanExportAndCurtail})
    println("CanExportAndCurtail")
end

function add_constraint(::Techs{CanExport})
    println("CanExport")
end

function add_constraint(::Techs{CanCurtail})
    println("CanCurtail")
end

function add_constraint(::Techs{CanNotExportCurtail})
    println("CanNotExportCurtail")
end

add_constraint(p)

ext: jl


text:
function hessian_matrix(poly, x)
    h = hessian(poly, x)
    return [
        h[1] h[2]
        h[2] h[3]
    ]
end

function saye_newton_iterate(
    xguess,
    xq,
    func,
    grad,
    hess,
    cellmap,
    tol,
    r;
    maxiter = 20,
    condtol = 1e5eps(),
)
    dim = length(xguess)
    jac = jacobian(cellmap)

    x0 = copy(xguess)
    gp = grad(x0)
    l0 = gp' * ((xq - cellmap(x0)) .* jac) / (gp' * gp)

    x1 = copy(x0)
    l1 = l0

    for counter = 1:maxiter
        vp = func(x0)
        gp = grad(x0)
        hp = hess(x0)

        gf = vcat(((cellmap(x0) - xq) .* jac) + l0 * gp, vp)
        hf = [
            diagm(jac .^ 2)+l0*hp gp
            gp' 0.0
        ]

        if inv(cond(hf)) > condtol
            δ = hf \ gf
            normδx = norm(δ[1:dim])
            if normδx > 0.5r
                δ *= 0.5r / normδx
            end
            x1 = x0 - δ[1:dim]
            l1 = l0 - δ[end]
        else
            error("Chopp method not implemented")
        end

        if norm(x1 - xguess) > r
            error("Did not converge in ball of radius $r")
        elseif norm(x1 - x0) < tol
            return x1
        else
            x0 = x1
            l0 = l1
        end
    end
    error("Did not converge in $maxiter iterations")
end

function project_on_zero_levelset(
    xguess,
    func,
    grad,
    tol,
    r;
    maxiter = 20,
    normtol = 1e-8,
    perturbation = 0.1,
)
    x0 = copy(xguess)
    x1 = copy(x0)
    dim = length(xguess)

    for counter = 1:maxiter
        vf = func(x0)
        gf = grad(x0)

        if norm(gf) < normtol
            x0 = x0 + perturbation * (rand(dim) .- 0.5)
        else
            δ = vf / (gf' * gf) * gf
            normδ = norm(δ)
            if normδ > 0.5r
                δ *= 0.5r / normδ
            end

            x1 = x0 - δ

            if abs(func(x1)) < tol
                flag = norm(x1 - xguess) < r
                return x1, flag
            else
                x0 = x1
            end
        end
    end
    error("Did not converge after $maxiter iterations")
end

function reference_seed_points(n)
    @assert n > 0
    xrange = range(-1.0, stop = 1.0, length = n + 2)
    points = ImplicitDomainQuadrature.tensor_product_points(xrange[2:n+1]', xrange[2:n+1]')
end

function seed_cell_zero_levelset(xguess, func, grad; tol = 1e-12, r = 2.5)
    dim, nump = size(xguess)
    pf = [project_on_zero_levelset(xguess[:, i], func, grad, tol, r) for i = 1:nump]
    flags = [p[2] for p in pf]
    valididx = findall(flags)
    validpoints = [p[1] for p in pf[valididx]]
    return hcat(validpoints...)
end

function seed_zero_levelset(nump, levelset, levelsetcoeffs, cutmesh)
    refpoints = reference_seed_points(nump)
    refseedpoints = []
    spatialseedpoints = []
    seedcellids = Int[]
    cellsign = cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)
    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        nodeids = nodal_connectivity(cutmesh.mesh, cellid)
        update!(levelset, levelsetcoeffs[nodeids])

        xk = seed_cell_zero_levelset(refpoints, levelset, x -> vec(gradient(levelset, x)))

        numseedpoints = size(xk)[2]

        append!(seedcellids, repeat([cellid], numseedpoints))
        push!(refseedpoints, xk)
        push!(spatialseedpoints, cellmap(xk))
    end

    refseedpoints = hcat(refseedpoints...)
    spatialseedpoints = hcat(spatialseedpoints...)
    return refseedpoints, spatialseedpoints, seedcellids
end

function seed_zero_levelset_with_interfacequads(interfacequads, cutmesh)
    cellsign = cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)

    totalnumqps = sum([length(interfacequads[1, cellid]) for cellid in cellids])

    refseedpoints = zeros(2, totalnumqps)
    spatialseedpoints = zeros(2, totalnumqps)
    seedcellids = zeros(Int, totalnumqps)

    start = 1
    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        refpoints = interfacequads[1, cellid].points
        spatialpoints = cellmap(refpoints)
        numqps = size(refpoints)[2]

        stop = start + numqps - 1

        refseedpoints[:, start:stop] = refpoints
        spatialseedpoints[:, start:stop] = spatialpoints
        seedcellids[start:stop] = repeat([cellid], numqps)

        start = stop + 1
    end
    return refseedpoints, spatialseedpoints, seedcellids
end

function collect_interface_normals(interfacequads,cutmesh)
    cellsign = cell_sign(cutmesh)
    cellids = findall(cellsign .== 0)

    totalnumqps = sum([length(interfacequads[1,cellid]) for cellid in cellids])

    interfacenormals = zeros(2,totalnumqps)
    start = 1

    for cellid in cellids
        normals = interface_normals(interfacequads,cellid)
        numqps = size(normals)[2]

        stop = start + numqps - 1

        interfacenormals[:,start:stop] = normals
        start = stop + 1
    end
    return interfacenormals
end

function closest_reference_points_on_zero_levelset(
    querypoints,
    refseedpoints,
    spatialseedpoints,
    seedcellids,
    levelset,
    levelsetcoeffs,
    mesh,
    tol,
    boundingradius,
)

    dim, numquerypoints = size(querypoints)
    refclosestpoints = zeros(dim, numquerypoints)
    refclosestcellids = zeros(Int, numquerypoints)
    refgradients = zeros(dim, numquerypoints)

    tree = KDTree(spatialseedpoints)
    seedidx, seeddists = nn(tree, querypoints)

    for (idx, sidx) in enumerate(seedidx)
        xguess = refseedpoints[:, sidx]
        xquery = querypoints[:, idx]
        guesscellid = seedcellids[sidx]
        cellmap = cell_map(mesh, guesscellid)
        update!(levelset, levelsetcoeffs[nodal_connectivity(mesh, guesscellid)])

        refcp = saye_newton_iterate(
            xguess,
            xquery,
            levelset,
            x -> vec(gradient(levelset, x)),
            x -> hessian_matrix(levelset, x),
            cellmap,
            tol,
            boundingradius,
        )

        refclosestpoints[:, idx] = refcp
        refclosestcellids[idx] = guesscellid
        refgradients[:, idx] = gradient(levelset, refcp)
    end
    return refclosestpoints, refclosestcellids, refgradients
end

function distance_to_zero_levelset(
    querypoints,
    refseedpoints,
    spatialseedpoints,
    seedcellids,
    levelset,
    levelsetcoeffs,
    mesh,
    tol,
    boundingradius,
)

    dim, numquerypoints = size(querypoints)
    signeddistance = zeros(numquerypoints)

    refclosestpoints, refclosestcellids, refgradients =
        closest_reference_points_on_zero_levelset(
            querypoints,
            refseedpoints,
            spatialseedpoints,
            seedcellids,
            levelset,
            levelsetcoeffs,
            mesh,
            tol,
            boundingradius,
        )

    for i = 1:numquerypoints
        refcp = refclosestpoints[:, i]
        cellmap = cell_map(mesh, refclosestcellids[i])

        spatialcp = cellmap(refcp)
        xquery = querypoints[:, i]

        g = vec(transform_gradient(refgradients[:, i]', jacobian(cellmap)))
        s = sign(g' * (xquery - spatialcp))

        signeddistance[i] = s * norm(spatialcp - xquery)
    end

    return signeddistance
end


text:
function one_point_crossover(circuit1,circuit2)
    karva1 = circuit1.karva
    karva2 = circuit2.karva
    crossover_point = rand(1:(length(karva1)-1))
    new_karva = karva1[1:crossover_point]*karva2[crossover_point+1:end]
    new_parameters =  vcat(circuit1.parameters[1:crossover_point],circuit2.parameters[crossover_point+1:end])
    return Circuit(new_karva,new_parameters,nothing)
end

function two_point_crossover(circuit1,circuit2)
    karva1 = circuit1.karva
    karva2 = circuit2.karva
    crossover_points = sort(rand(1:length(karva1)-1,2))
    new_karva = karva1[1:crossover_points[1]]*karva2[crossover_points[1]+1:crossover_points[2]]*karva1[crossover_points[2]+1:end]
    new_parameters = vcat(circuit1.parameters[1:crossover_points[1]],circuit2.parameters[crossover_points[1]+1:crossover_points[2]],circuit1.parameters[crossover_points[2]+1:end])
    return Circuit(new_karva,new_parameters,nothing)
end

function mutate(circuit,terminals = "RCLP")
    karva = circuit.karva ; parameters = copy(circuit.parameters)
    total_length = length(karva)
    head_length = (total_length-1)/2 ; tail_length = head_length+1
    operations = "+-"; new_element = ""
    position = rand(1:total_length)
    mutable_karva = collect(karva)
    if position == 1
        new_element = rand(operations)
        mutable_karva[position] = new_element
    elseif position <= head_length
        new_element = rand(operations*terminals)
        mutable_karva[position] = new_element
        parameters[position] = replacementparameter(new_element)
    elseif position > head_length
        new_element = rand(terminals)
        mutable_karva[position] = new_element
        parameters[position] = replacementparameter(new_element)
    end
    return Circuit(String(mutable_karva),parameters,nothing)
end

function replacementparameter(element)
    ranges= Dict('R'=>1000,'C'=>0.001,'L'=>1,'P'=>[100,0.99999],'+'=>0,'-'=>0)
    return rand().*ranges[element]
end

function multipoint_mutation(circuit,N,terminals="RCLP")
    for i in 1:N
        circuit = mutate(circuit,terminals)
    end
    return circuit
end

function transposition(circuit)
    parameters = copy(circuit.parameters)
    coding_length = length(circuit.karva)
    startsite = rand(2:coding_length-1)
    max_transposonlength = coding_length-startsite
    transposonlength = rand(1:max_transposonlength)
    targetsite = rand(2:startsite)
    mutable_karva = collect(circuit.karva)
    mutable_karva[targetsite:targetsite+transposonlength] = mutable_karva[startsite:startsite+transposonlength]
    parameters[targetsite:targetsite+transposonlength] = parameters[startsite:startsite+transposonlength]
    return  Circuit(String(mutable_karva),parameters,nothing)
end

function tournamentselection(population,mating_pool_size,tournament_size)
    selected = []
    for i in 1:mating_pool_size
        tournament = rand(population,tournament_size)
        push!(selected,minimum(tournament))
    end
    return selected
end

function truncationselection(population,truncation_size)
    sort!(population)
    return population[1:truncation_size]
end


text:
using LinearAlgebra, InvertibleNetworks, ForwardDiff, Test, PyPlot

# Input dimension
n = 10:10:300
t_noAD = Array{Float32, 1}(undef, length(n))
t_AD = Array{Float32, 1}(undef, length(n))
n_in = 3
batchsize = 5

# # of evaluations
neval = 10

# Actnorm parameters
s0 = randn(Float32, 1, 1, 1, n_in, 1)
b0 = randn(Float32, 1, 1, 1, n_in, 1)
ds = randn(Float32, size(s0))
db = randn(Float32, size(b0))

# Actnorm jacobian (w/out AD)
function jacobian(ΔX::Array{Float32, 5}, Δs::Array{Float32, 5}, Δb::Array{Float32, 5}, X::Array{Float32, 5}, s::Array{Float32, 5}, b::Array{Float32, 5})
    Y = X .* s .+ b
    ΔY = X .* Δs .+ ΔX .* s .+ Δb
    return ΔY, Y
end

# Loop over input size
for i = 1:length(n)

    # Input
    nx = n[i]
    ny = n[i]
    nz = n[i]
    X0 = randn(Float32, nx, ny, nz, n_in, batchsize)
    dX = randn(Float32, size(X0))

    # Actnorm jacobian (w/ AD)
    function fun(X, s, b)
        return X .* s .+ b
    end
    jacobianAD = t -> ForwardDiff.derivative(t -> fun(X0+t*dX, s0+t*ds, b0+t*db), t)

    # Eval
    dfun = jacobianAD(0f0)
    dfun_ = jacobian(dX, ds, db, X0, s0, b0)[1]

    # Discrepancy
    @test isapprox(dfun, dfun_; rtol=1f-3)

    # Timings
    t_AD[i] = (@elapsed for i = 1:neval jacobianAD(0f0); end;)/neval
    t_noAD[i] = (@elapsed for i = 1:neval jacobian(dX, ds, db, X0, s0, b0)[1]; end;)/neval

    # Print msg
    println("[", i, "/", length(n), "] --- size= ", n[i], ", time= AD:", t_AD[i], ", no AD: ", t_noAD[i])

end

# Plotting
figure()
title(string("Timings for ActNorm jacobian evaluations (AD vs no AD); n_ch=", n_in, ", batchsize=", batchsize))
loglog(n, t_AD)
loglog(n, t_noAD)
legend(["AD", "no AD"])
xlabel("n (size(input)=n^3)")
grid("on")
savefig("AD_vs_noAD.png")


text:
function voxel_structure(	particle_type::String,
						R::Array{Float64, 2},
						Lx::Float64,
						Ly::Float64,
						Lz::Float64,
						X::Array{Float64, 1},
						Y::Array{Float64, 1},
						Z::Array{Float64, 1},
						Q0::Array{Float64, 1},
						Q1::Array{Float64, 1},
						Q2::Array{Float64, 1},
						Q3::Array{Float64, 1},
						voxel_size::Float64)

	Nx::Int64 = convert(Int64, round(Lx/voxel_size))
	Ny::Int64 = convert(Int64, round(Ly/voxel_size))
	Nz::Int64 = convert(Int64, round(Lz/voxel_size))

	M::BitArray{3} = falses(Nx, Ny, Nz)

    X /= voxel_size
	Y /= voxel_size
	Z /= voxel_size
    R /= voxel_size
	Lx = convert(Float64, Nx)
	Ly = convert(Float64, Ny)
	Lz = convert(Float64, Nz)

	number_of_particles::Int64 = length(X)

	# Pre-computed max radii (i.e. radius of bounding sphere).
	RMAX::Array{Float64, 1} = zeros(number_of_particles)
	if particle_type == "sphere"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = R[current_particle, 1]
		end
	elseif particle_type == "ellipsoid"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = maximum( R[current_particle, 1:3] )
		end
	elseif particle_type == "cuboid"
		for current_particle = 1:number_of_particles
			RMAX[current_particle] = sqrt(R[current_particle, 1]^2 + R[current_particle, 2]^2 + R[current_particle, 3]^2)
		end
	end

	x::Float64 = 0.0
	y::Float64 = 0.0
	z::Float64 = 0.0
	x_min::Int64 = 0
	x_max::Int64 = 0
	y_min::Int64 = 0
	y_max::Int64 = 0
	z_min::Int64 = 0
	z_max::Int64 = 0
	a11::Float64 = 0.0
	a12::Float64 = 0.0
	a13::Float64 = 0.0
	a21::Float64 = 0.0
	a22::Float64 = 0.0
	a23::Float64 = 0.0
	a31::Float64 = 0.0
	a32::Float64 = 0.0
	a33::Float64 = 0.0
	f::Float64 = 0.0

	x_prim::Float64 = 0.0
	y_prim::Float64 = 0.0
	z_prim::Float64 = 0.0
	X_prim::Array{Float64, 1} = zeros(3)

	for current_particle = 1:number_of_particles
		for i = -1:1
			for j = -1:1
				for k = -1:1
					x = X[current_particle] + convert(Float64, i) * Lx
					y = Y[current_particle] + convert(Float64, j) * Ly
					z = Z[current_particle] + convert(Float64, k) * Lz
					if all([x, y, z] .< [Lx, Ly, Lz] + RMAX[current_particle]) && all([x, y, z] .> - RMAX[current_particle])
						x_min = max(1, convert(Int64, floor(x - RMAX[current_particle])))
						x_max = min(Nx, convert(Int64, ceil(x + RMAX[current_particle])))
						y_min = max(1, convert(Int64, floor(y - RMAX[current_particle])))
						y_max = min(Ny, convert(Int64, ceil(y + RMAX[current_particle])))
						z_min = max(1, convert(Int64, floor(z - RMAX[current_particle])))
						z_max = min(Nz, convert(Int64, ceil(z + RMAX[current_particle])))

						for xx = x_min:x_max
							for yy =  y_min:y_max
								for zz =  z_min:z_max
									if particle_type == "sphere"
										if (x - convert(Float64, xx))^2 + (y - convert(Float64, yy))^2 + (z - convert(Float64, zz))^2 <= R[current_particle]^2
											M[xx, yy, zz] = true
										end
									elseif particle_type == "ellipsoid"
										(a11, a12, a13, a21, a22, a23, a31, a32, a33) = characteristic_matrix_ellipsoid(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle], R[current_particle, 1], R[current_particle, 2], R[current_particle, 3])
										f = [x - convert(Float64, xx), y - convert(Float64, yy), z - convert(Float64, zz)]' * ([a11 a12 a13; a21 a22 a23; a31 a32 a33] \ [x - convert(Float64, xx), y - convert(Float64, yy), z - convert(Float64, zz)])
										if f <= 1.0
											M[xx, yy, zz] = true
										end
									elseif particle_type == "cuboid"
										(a11, a12, a13, a21, a22, a23, a31, a32, a33) = rotation_matrix(Q0[current_particle], Q1[current_particle], Q2[current_particle], Q3[current_particle])
										X_prim = [a11 a12 a13; a21 a22 a23; a31 a32 a33] \ [x - convert(Float64, xx), y - convert(Float64, yy), z - convert(Float64, zz)]

										x_prim = X_prim[1]#a11 * (x - convert(Float64, xx)) + a12 * (y - convert(Float64, yy)) + a13 * (z - convert(Float64, zz))
										y_prim = X_prim[2]#a21 * (x - convert(Float64, xx)) + a22 * (y - convert(Float64, yy)) + a23 * (z - convert(Float64, zz))
										z_prim = X_prim[3]#a31 * (x - convert(Float64, xx)) + a32 * (y - convert(Float64, yy)) + a33 * (z - convert(Float64, zz))
										if abs(x_prim) <= R[current_particle, 1] && abs(y_prim) <= R[current_particle, 2] && abs(z_prim) <= R[current_particle, 3]
											M[xx, yy, zz] = true
										end
									end
								end
							end
						end
					end
				end
			end
		end
	end

	return M
end

max_stars_repo_name: roding/whitefish

max_issues_repo_name: roding/whitefish

max_forks_repo_name: roding/whitefish


text:
# Aliases (types)
const IntVector    = Array{Int64,1};
const IntMatrix    = Array{Int64,2};
const IntArray     = Array{Int64};
const FloatVector  = Array{Float64,1};
const FloatMatrix  = Array{Float64,2};
const FloatArray   = Array{Float64};
const JVector{T}   = Array{Union{Missing, T}, 1};
const JMatrix{T}   = Array{Union{Missing, T}, 2};
const JArray{T, N} = Array{Union{Missing, T}, N};
const DiagMatrix   = Diagonal{Float64,Array{Float64,1}};
const SymMatrix    = Symmetric{Float64,Array{Float64,2}};

#=
--------------------------------------------------------------------------------------------------------------------------------
Kalman structures
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    MSeries(...)

Convenient mutable struct used in KalmanSettings for describing the data in a dynamic fashion.

# Arguments
- `data`: Observed measurements (`nxT`)
- `n`: Number of series
- `T`: Number of observations
"""
mutable struct MSeries
    data::Union{FloatMatrix, JMatrix{Float64}}
    n::Int64
    T::Int64
end

"""
    KalmanSettings(...)

Define an immutable structure that includes all the Kalman filter and smoother inputs.

# Model
The state space model used below is,

``Y_{t} = B*X_{t} + e_{t}``

``X_{t} = C*X_{t-1} + D*U_{t}``

where ``e_{t} ~ N(0_{nx1}, R)`` and ``U_{t} ~ N(0_{mx1}, Q)``.

# Arguments
- `Y`: Observed measurements (`nxT`)
- `B`: Measurement equations' coefficients
- `R`: Covariance matrix of the measurement equations' error terms
- `C`: Transition equations' coefficients
- `D`: Transition equations' coefficients associated to the error terms
- `Q`: Covariance matrix of the transition equations' error terms
- `X0`: Mean vector for the states at time t=0
- `P0`: Covariance matrix for the states at time t=0
- `DQD`: Covariance matrix of D*U_{t} (i.e., D*Q*D')
- `m`: Number of latent states
- `compute_loglik`: Boolean (true for computing the loglikelihood in the Kalman filter)
- `store_history`: Boolean (true to store the history of the filter and smoother)
"""
struct KalmanSettings
    Y::MSeries
    B::FloatMatrix
    R::SymMatrix
    C::FloatMatrix
    D::FloatMatrix
    Q::SymMatrix
    X0::FloatVector
    P0::SymMatrix
    DQD::SymMatrix
    m::Int64
    compute_loglik::Bool
    store_history::Bool
end

#=
KalmanSettings constructors
=#

function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B, 2);
    D = Matrix(1.0I, m, m);
    X0 = zeros(m);
    P0 = solve_discrete_lyapunov(C, Q);

    # Return KalmanSettings
    return KalmanSettings(MSeries(Y, n, T), B, R, C, D, Q, X0, P0, Q, m, compute_loglik, store_history);
end

function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B, 2);
    X0 = zeros(m);
    DQD = Symmetric(D*Q*D');
    P0 = solve_discrete_lyapunov(C, DQD);

    # Return KalmanSettings
    return KalmanSettings(MSeries(Y, n, T), B, R, C, D, Q, X0, P0, DQD, m, compute_loglik, store_history);
end

function KalmanSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, B::FloatMatrix, R::SymMatrix, C::FloatMatrix, D::FloatMatrix, Q::SymMatrix, X0::FloatVector, P0::SymMatrix; compute_loglik::Bool=true, store_history::Bool=true)

    # Compute default value for missing parameters
    n, T = size(Y);
    m = size(B, 2);
    DQD = Symmetric(D*Q*D');

    # Return KalmanSettings
    return KalmanSettings(MSeries(Y, n, T), B, R, C, D, Q, X0, P0, DQD, m, compute_loglik, store_history);
end

abstract type KalmanStatus end

"""
    OnlineKalmanStatus(...)

Define a mutable structure to handle minimal Kalman filter output.

The Kalman filter history is not stored. This makes it ideal for online filtering problems. However, it is incompatible with the kalman smoother implementation.

# Arguments
- `t`: Current point in time
- `loglik`: Loglikelihood
- `X_prior`: Latest a-priori X
- `X_post`: Latest a-posteriori X
- `P_prior`: Latest a-priori P
- `P_post`: Latest a-posteriori P
- `e`: Forecast error
- `inv_F`: Inverse of the forecast error covariance
- `L`: Convenient shortcut for the filter and smoother
- `buffer_J1`, `buffer_J2`, `buffer_m_m`, `buffer_m_n_obs`: Buffer arrays used for low-level matrix operations
"""
mutable struct OnlineKalmanStatus <: KalmanStatus
    t::Int64
    loglik::Union{Float64, Nothing}
    X_prior::Union{FloatVector, Nothing}
    X_post::Union{FloatVector, Nothing}
    P_prior::Union{SymMatrix, Nothing}
    P_post::Union{SymMatrix, Nothing}
    e::Union{FloatVector, Nothing}
    inv_F::Union{SymMatrix, Nothing}
    L::Union{FloatMatrix, Nothing}
    buffer_J1::Union{FloatVector, Nothing}
    buffer_J2::Union{FloatMatrix, Nothing}
    buffer_m_m::Union{FloatMatrix, Nothing}
    buffer_m_n_obs::Union{FloatMatrix, Nothing}
end

# OnlineKalmanStatus constructor
OnlineKalmanStatus() = OnlineKalmanStatus(0, [nothing for i=1:12]...);

"""
    DynamicKalmanStatus(...)

Define a mutable structure to handle any type of Kalman filter output.

The Kalman filter history is stored when `store_history` is set to true in the filter settings.

# Arguments
- `t`: Current point in time
- `loglik`: Loglikelihood
- `X_prior`: Latest a-priori X
- `X_post`: Latest a-posteriori X
- `P_prior`: Latest a-priori P
- `P_post`: Latest a-posteriori P
- `e`: Forecast error
- `inv_F`: Inverse of the forecast error covariance
- `L`: Convenient shortcut for the filter and smoother
- `buffer_J1`, `buffer_J2`, `buffer_m_m`, `buffer_m_n_obs`: Buffer arrays used for low-level matrix operations
- `history_X_prior`: History of a-priori X
- `history_X_post`: History of a-posteriori X
- `history_P_prior`: History of a-priori P
- `history_P_post`: History of a-posteriori P
- `history_e`: History of the forecast error
- `history_inv_F`: History of the inverse of the forecast error covariance
- `history_L`: History of the shortcut L
"""
mutable struct DynamicKalmanStatus <: KalmanStatus
    t::Int64
    loglik::Union{Float64, Nothing}
    X_prior::Union{FloatVector, Nothing}
    X_post::Union{FloatVector, Nothing}
    P_prior::Union{SymMatrix, Nothing}
    P_post::Union{SymMatrix, Nothing}
    e::Union{FloatVector, Nothing}
    inv_F::Union{SymMatrix, Nothing}
    L::Union{FloatMatrix, Nothing}
    buffer_J1::Union{FloatVector, Nothing}
    buffer_J2::Union{FloatMatrix, Nothing}
    buffer_m_m::Union{FloatMatrix, Nothing}
    buffer_m_n_obs::Union{FloatMatrix, Nothing}
    history_X_prior::Union{Array{FloatVector,1}, Nothing}
    history_X_post::Union{Array{FloatVector,1}, Nothing}
    history_P_prior::Union{Array{SymMatrix,1}, Nothing}
    history_P_post::Union{Array{SymMatrix,1}, Nothing}
    history_e::Union{Array{FloatVector,1}, Nothing}
    history_inv_F::Union{Array{SymMatrix,1}, Nothing}
    history_L::Union{Array{FloatMatrix,1}, Nothing}
end

# DynamicKalmanStatus constructor
DynamicKalmanStatus() = DynamicKalmanStatus(0, [nothing for i=1:19]...);

"""
    SizedKalmanStatus(...)

Define an immutable structure that always store the filter history up to time T.

# Arguments
- `online_status`: `OnlineKalmanStatus` struct
- `history_length`: History length
- `history_X_prior`: History of a-priori X
- `history_X_post`: History of a-posteriori X
- `history_P_prior`: History of a-priori P
- `history_P_post`: History of a-posteriori P
- `history_e`: History of the forecast error
- `history_inv_F`: History of the inverse of the forecast error covariance
- `history_L`: History of the shortcut L
"""
struct SizedKalmanStatus <: KalmanStatus
    online_status::OnlineKalmanStatus
    history_length::Int64
    history_X_prior::Array{FloatVector,1}
    history_X_post::Array{FloatVector,1}
    history_P_prior::Array{SymMatrix,1}
    history_P_post::Array{SymMatrix,1}
    history_e::Array{FloatVector,1}
    history_inv_F::Array{SymMatrix,1}
    history_L::Array{FloatMatrix,1}
end

# SizedKalmanStatus constructor
function SizedKalmanStatus(T::Int64)
    history_X_prior = Array{FloatVector,1}(undef, T);
    history_X_post = Array{FloatVector,1}(undef, T);
    history_P_prior = Array{SymMatrix,1}(undef, T);
    history_P_post = Array{SymMatrix,1}(undef, T);
    history_e = Array{FloatVector,1}(undef, T);
    history_inv_F = Array{SymMatrix,1}(undef, T);
    history_L = Array{FloatMatrix,1}(undef, T);
    return SizedKalmanStatus(OnlineKalmanStatus(), T, history_X_prior, history_X_post, history_P_prior, history_P_post, history_e, history_inv_F, history_L);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
UC models: structures
--------------------------------------------------------------------------------------------------------------------------------
=#

abstract type UCSettings end

"""
    VARIMASettings(...)

Define an immutable structure to manage VARIMA specifications.

# Arguments
- `Y_levels`: Observed measurements (`nxT`) - in levels
- `Y`: Observed measurements (`nxT`) - differenced and demeaned
- `μ`: Sample average (per series)
- `n`: Number of series
- `d`: Degree of differencing
- `p`: Order of the autoregressive model
- `q`: Order of the moving average model
- `nr`: n*max(p, q+1)
- `np`: n*p
- `nq`: n*q
- `nnp`: (n^2)*p
- `nnq`: (n^2)*q
"""
struct VARIMASettings <: UCSettings
    Y_levels::Union{FloatMatrix, JMatrix{Float64}}
    Y::Union{FloatMatrix, JMatrix{Float64}}
    μ::FloatVector
    n::Int64
    d::Int64
    p::Int64
    q::Int64
    nr::Int64
    np::Int64
    nq::Int64
    nnp::Int64
    nnq::Int64
end

# VARIMASettings constructor
function VARIMASettings(Y_levels::Union{FloatMatrix, JMatrix{Float64}}, d::Int64, p::Int64, q::Int64)

    # Initialise
    n = size(Y_levels, 1);
    r = max(p, q+1);

    # Differenciate data
    Y = copy(Y_levels);

    if d > 0
        for i=1:d
            Y = diff(Y, dims=2);
        end
    end

    # Mean
    μ = mean_skipmissing(Y)[:,1];

    # Demean data
    Y = demean(Y);

    # VARIMASettings
    return VARIMASettings(Y_levels, Y, μ, n, d, p, q, n*r, n*p, n*q, (n^2)*p, (n^2)*q);
end

"""
    ARIMASettings(...)

Define an alias of VARIMASettings for arima models.
"""
ARIMASettings(Y_levels::Union{FloatMatrix, JMatrix{Float64}}, d::Int64, p::Int64, q::Int64) = VARIMASettings(Y_levels, d, p, q);


text:
include("WebsocketFrame.jl")
const ioTypes = Union{Nothing, WebsocketFrame, HTTP.ConnectionPool.Transaction, Timer, Closereason, Array{UInt8, 1}}

struct WebsocketConnection
    id::String
    config::NamedTuple
    io::Dict{Symbol, ioTypes}
    callbacks::Dict{Symbol, Union{Bool, Function}}
    buffers::NamedTuple
    closed::Condition
    keepalive::Dict{Symbol, Union{String, Bool}}
    clients::Union{Array{WebsocketConnection,1}, Nothing}
    messageChannel::Channel{Union{String, Array{UInt8,1}}}
    pongChannel::Channel{Union{String, Array{UInt8,1}}}
    validate::Dict{String, Any}

    function WebsocketConnection(
        stream::HTTP.ConnectionPool.Transaction,
        config::NamedTuple,
        clients::Union{Array{WebsocketConnection,1}, Nothing} = nothing
    )
        @debug "WebsocketConnection"
        buffers = (
            maskBytes = IOBuffer(; maxsize = 4),
            frameHeader = IOBuffer(; maxsize = 10),
            outBuffer = config.fragmentOutgoingMessages ? IOBuffer(; maxsize = Int(config.fragmentationThreshold) + 10) : IOBuffer(),
            inBuffer = IOBuffer(; maxsize = getmaxwithheaders(config)),
            fragmentBuffer = IOBuffer(; maxsize = Int(config.maxReceivedMessageSize))
        )

        initmessage = "keepalive-"*string(rand(UInt32))
        keepalive = Dict{Symbol, Union{String, Bool}}(
            :pingmessage => initmessage,
            :pongmessage => initmessage,
            :isopen => true,
            :isalive => true,
        )
        interval = config.keepaliveTimeout
        interval !== false && (keepaliveTimer = Timer(timer -> (
            try
                if !keepalive[:isopen]
                    close(timer)
                    return
                end
                if keepalive[:pingmessage] !== keepalive[:pongmessage]
                    self.io[:closeReason] = Closereason(CLOSE_REASON_ABNORMAL, "could not ping the $(config.type === "client" ? "server" : "client").")
                    close(self.io[:stream])
                    gracefulEnd(self)
                    close(timer)
                elseif !keepalive[:isalive]
                    keepalive[:pingmessage] = "keepalive-"*string(rand(UInt32))
                    ping(self, keepalive[:pingmessage])
                end
                keepalive[:isalive] = false
            catch
            end
        ), interval; interval = interval))

        @async begin
            reason = wait(self.closed)
            clients !== nothing && filter!(client -> (client.id !== self.id), clients)
            for buffer in collect(self.buffers)
                close(buffer)
            end
            closecallback = self.callbacks[:close]
            if closecallback isa Function
                closecallback((; code = reason.code, description = reason.description))
            elseif config.type === "client"
                @warn "$(config.type) websocket connection closed." code = reason.code description = reason.description
            end
        end

        atexit(() -> (
            if isopen(self.io[:stream])
                self.io[:closeReason] = Closereason(CLOSE_REASON_ABNORMAL, "julia process exited.")
                close(self.io[:stream])
                wait(self.closed)
                sleep(0.001)
            end
        ))

        self = new(
            requestHash(),
            config,                                                     #config
            Dict{Symbol, ioTypes}(                                      #io
                :stream => stream,
                :currentFrame => WebsocketFrame(config, buffers),
                :closeTimeout => nothing,
                :closeReason => nothing,
                :stash => Array{UInt8, 1}()
            ),
            Dict{Symbol, Union{Bool, Function}}(                        #callbacks
                :message => false,
                :error => false,
                :close => false,
                :pong => false
            ),
            buffers,                                                    #buffers
            Condition(),                                                #closed
            keepalive,                                                  #keepalive
            clients,                                                    #clients
            Channel{Union{String, Array{UInt8,1}}}(Inf),                #messageChannel
            Channel{Union{String, Array{UInt8,1}}}(Inf),                #pongChannel
            Dict{String, Any}("valid" => true)                          #validate
        )
    end
end
"""
    listen(callback::Function, ws::SimpleWebsockets.WebsocketConnection, event::Symbol)
Register event callbacks onto a `WebsocketConnection`. The callback must be a function with exactly one argument.

Valid events are:
- :message
- :pong
- :error
- :close

# Example
```julia
listen(ws::WebsocketConnection, :message) do message
    #...
end
```
"""
function listen(
    cb::Function,
    self::WebsocketConnection,
    key::Symbol
)
    if !haskey(self.callbacks, key)
        return @warn "WebsocketConnection has no listener for :$key."
    end

    if haskey(self.callbacks, key)
        function errorHandler(err, trace)
            err = CallbackError(err, trace)
            errcallback = self.callbacks[:error]
            if errcallback !== false && key !== :error
                errcallback(err)
            else
                err.log()
                exit()
            end
        end
        if !(self.callbacks[key] isa Function)
            if key === :message
                self.callbacks[key] = data -> ()
                @async while isopen(self.messageChannel)
                    message = take!(self.messageChannel)
                    try
                        cb(message)
                    catch err
                        errorHandler(err, catch_backtrace())
                    end
                end
            elseif key === :pong
                self.callbacks[key] = data -> ()
                @async while isopen(self.pongChannel)
                    message = take!(self.pongChannel)
                    try
                        cb(message)
                    catch err
                        errorHandler(err, catch_backtrace())
                    end
                end
            else
                self.callbacks[key] = data -> (
                    @async try
                        cb(data)
                    catch err
                        errorHandler(err, catch_backtrace())
                    end
                )
            end
        end
    end
end

function gethandles(connection::WebsocketConnection, io::HTTP.Streams.Stream)
    if connection.config.type === "client"
        (; file = io, available = io,)
    else
        (; file = io.stream, available = io.stream.c.io)
    end
end
function startConnection(self::WebsocketConnection, io::HTTP.Streams.Stream)
    handle = gethandles(self, io)
    while !eof(handle.file)
        data = readavailable(handle.available)
        isopen(self.io[:stream]) && try
            if length(data) > self.buffers.inBuffer.maxsize
                throw(error("[$CLOSE_REASON_MESSAGE_TOO_BIG]Maximum message size of $(self.config.maxReceivedMessageSize) Bytes exceeded"))
            end
            seekend(self.buffers.inBuffer)
            unsafe_write(self.buffers.inBuffer, pointer(data), length(data))
            processReceivedData(self)
        catch err
            err = FrameError(err, catch_backtrace())
            if self.callbacks[:error] isa Function
                self.callbacks[:error](err)
            elseif self.config.type === "client"
                err.log()
            else
                @debug "CLOSE_REASON_INVALID_DATA" error = err.msg
            end
            closeConnection(self, CLOSE_REASON_INVALID_DATA, err.msg)
        end
    end
    close(self.messageChannel)
    close(self.pongChannel)
    self.keepalive[:isopen] && gracefulEnd(self)
end

function gracefulEnd(self::WebsocketConnection)
    self.keepalive[:isopen] = false
    @async begin
        (   self.io[:closeTimeout] !== nothing &&
            isopen(self.io[:closeTimeout])
        ) && sleep(0.001)

        if self.io[:closeReason] === nothing
            self.io[:closeReason] = Closereason(CLOSE_REASON_ABNORMAL)
        end
        notify(self.closed, self.io[:closeReason]; all = true)
    end
end

function closeConnection(self::WebsocketConnection, reasonCode::Int, reason::String)
    closereason = Closereason(reasonCode, reason)
    self.io[:closeReason] = closereason
    !closereason.valid && throw(error("invalid connection close code."))

    if isopen(self.io[:stream])
        sendCloseFrame(self, closereason.code, closereason.description)

        self.io[:closeTimeout] = Timer(timer -> (
            try
                isopen(self.io[:stream]) && close(self.io[:stream])
            catch
            end
        ), self.config.closeTimeout)

    elseif self.io[:closeTimeout] isa Timer && isopen(self.io[:closeTimeout])
        close(self.io[:closeTimeout])
    end
end

# Send data
"""
    send(ws::WebsocketConnection, data::Union{String, Number, Array{UInt8,1}})
Send the given `data` as a message over the wire.
"""
send(self::WebsocketConnection, data::Number) = send(self, string(data))
send(self::WebsocketConnection, data::String) = send(self, textbuffer(data))
function send(self::WebsocketConnection, data::Array{UInt8,1})
    @debug "WebsocketConnection.send"
    try

        #BUGFIX
        #Julia version < 1.2 gets confused with allocation under fringe circumstances
            VERSION < v"1.2.0" && (data = copy(data))
        #END

        frame = WebsocketFrame(self.config, self.buffers, data)
        frame.inf[:opcode] = self.config.binary ? BINARY_FRAME : TEXT_FRAME
        fragmentAndSend(self, frame)
    catch err
        err = FrameError(err, catch_backtrace())
        if self.callbacks[:error] isa Function
            self.callbacks[:error](err)
        else
            err.log()
        end
    end
end

function sendCloseFrame(self::WebsocketConnection, reasonCode::Int, description::String)
    frame = WebsocketFrame(self.config, self.buffers, textbuffer(description));
    frame.inf[:opcode] = CONNECTION_CLOSE_FRAME
    frame.inf[:fin] = true
    frame.inf[:closeStatus] = reasonCode
    sendFrame(self, frame)
end
"""
    ping(ws::WebsocketConnection, data::Union{String, Number})
Send `data` as ping message to the `ws` peer.

`data` is limited to 125Bytes, and will automatically be truncated if over this limit.
"""
ping(self::WebsocketConnection, data::Number) = ping(self, String(data))
function ping(self::WebsocketConnection, data::String)
    try
        data = textbuffer(data)
        size(data, 1) > 125 && (data = data[1:125, 1])
        frame = WebsocketFrame(self.config, self.buffers, data)
        frame.inf[:fin] = true
        frame.inf[:opcode] = PING_FRAME
        sendFrame(self, frame)
    catch err
        err = FrameError(err, catch_backtrace())
        if self.callbacks[:error] isa Function
            self.callbacks[:error](err)
        else
            err.log()
        end
    end
end
pong(self::WebsocketConnection, data::Union{String, Number}) = pong(self, textbuffer(data))
function pong(self::WebsocketConnection, data::Array{UInt8,1})
    size(data, 1) > 125 && (data = data[1:125, 1])
    frame = WebsocketFrame(self.config, self.buffers, data)
    frame.inf[:fin] = true
    frame.inf[:opcode] = PONG_FRAME
    sendFrame(self, frame)
end
function fragmentAndSend(
    self::WebsocketConnection,
    frame::WebsocketFrame,
    numFragments::Int,
    fragment::Int
)
    fragment +=1
    binaryPayload = frame.inf[:binaryPayload]
    endIndex = Int(self.config.fragmentationThreshold)

    if length(binaryPayload) > endIndex
        partialPayload = splice!(binaryPayload, 1:endIndex)
    else
        partialPayload = binaryPayload
    end
    partframe = WebsocketFrame(self.config, self.buffers, partialPayload)
    if fragment === 1
        partframe.inf[:opcode] = frame.inf[:opcode]
    else
        partframe.inf[:opcode] = CONTINUATION_FRAME
    end
    partframe.inf[:fin] = fragment === numFragments
    sendFrame(self, partframe)
    fragment < numFragments && fragmentAndSend(self, frame, numFragments, fragment)
end
function fragmentAndSend(self::WebsocketConnection, frame::WebsocketFrame)
    @debug "WebsocketConnection.fragmentAndSend"
    threshold = self.config.fragmentationThreshold
    len = length(frame.inf[:binaryPayload])
    if !self.config.fragmentOutgoingMessages || len <= threshold
        frame.inf[:fin] = true
        sendFrame(self, frame)
        return
    end
    numFragments = Int(ceil(len / threshold))
    fragmentAndSend(self, frame, numFragments, 0)
end

function sendFrame(self::WebsocketConnection, frame::WebsocketFrame)
    frame.inf[:mask] = self.config.maskOutgoingPackets

    proceed = toBuffer(frame)

    outBuffer = self.buffers.outBuffer
    @debug "WebsocketConnection.sendFrame" (;
        type = self.config.type,
        message = String(copy(frame.inf[:binaryPayload])),
        id = self.id,
        proceed = proceed,
        ptr = outBuffer.ptr,
        buffersize = outBuffer.size,
        parseState = frame.inf[:parseState],
        opcode = frame.inf[:opcode],
        length = Int(frame.inf[:length]),
        fin = frame.inf[:fin],
        mask = frame.inf[:mask],
    )...

    !proceed && return
    isopen(self.io[:stream]) && write(self.io[:stream], outBuffer)
end
# End send data

# Receive data
function processReceivedData(self::WebsocketConnection)
    frame = self.io[:currentFrame]
    inBuffer = self.buffers.inBuffer

    continued = addData(frame)

    @debug "WebsocketConnection.processReceivedData" (;
        type = self.config.type,
        message = String(copy(frame.inf[:binaryPayload])),
        id = self.id,
        continued = continued,
        ptr = inBuffer.ptr,
        buffersize = inBuffer.size,
        parseState = frame.inf[:parseState],
        opcode = frame.inf[:opcode],
        length = Int(frame.inf[:length]),
        fin = frame.inf[:fin],
        mask = frame.inf[:mask],
    )...

    !continued && return

    processFrame(self, frame)
    self.io[:currentFrame] = WebsocketFrame(
        self.config,
        self.buffers;
            ptr = frame.inf[:ptr]
    )

    if (inBuffer.ptr - 1) < inBuffer.size
        processReceivedData(self)
    else
        self.io[:currentFrame].inf[:ptr] = 1
        isopen(inBuffer) && truncate(inBuffer, 0)
    end
end

function processFrame(self::WebsocketConnection, frame::WebsocketFrame)
    inf = frame.inf
    opcode = inf[:opcode]
    fragmentBuffer = self.buffers.fragmentBuffer
    binary = self.config.binary
    data = inf[:binaryPayload]
    if fragmentBuffer.size > 0 && (opcode > 0x00 && opcode < 0x08)
        throw(error("illegal frame opcode $opcode received in middle of fragmented message."))
    end
    self.keepalive[:isalive] = true
    if opcode === TEXT_FRAME || opcode === BINARY_FRAME
        if length(data) > self.config.maxReceivedMessageSize
            throw(error("[$CLOSE_REASON_MESSAGE_TOO_BIG]Maximum message size of $(self.config.maxReceivedMessageSize) Bytes exceeded"))
        end
        if frame.inf[:fin]
            put!(self.messageChannel, binary ? data : String(data))
        else
            unsafe_write(fragmentBuffer, pointer(data), length(data))
        end
    elseif opcode === CONTINUATION_FRAME

        if fragmentBuffer.size + length(data) > self.config.maxReceivedMessageSize
            throw(error("[$CLOSE_REASON_MESSAGE_TOO_BIG]Maximum message size of $(self.config.maxReceivedMessageSize) Bytes exceeded"))
        end
        unsafe_write(fragmentBuffer, pointer(data), length(data))
        if inf[:fin]
            seekstart(fragmentBuffer)
            data = binary ? read(fragmentBuffer) : read(fragmentBuffer, String)
            isopen(fragmentBuffer) && truncate(fragmentBuffer, 0)

            put!(self.messageChannel, binary ? data : String(data))
        end

    elseif opcode === PING_FRAME
        pong(self, inf[:binaryPayload])
    elseif opcode === PONG_FRAME
        callback = self.callbacks[:pong]
        message = String(frame.inf[:binaryPayload])

        if message === self.keepalive[:pingmessage]
            self.keepalive[:pongmessage] = message
        elseif callback isa Function
            put!(self.pongChannel, binary ? textbuffer(message) : message)
        end
    elseif opcode === CONNECTION_CLOSE_FRAME
        description = String(inf[:binaryPayload])
        self.io[:closeReason] = Closereason(inf[:closeStatus], description)
        if !self.io[:closeReason].valid
            @warn "received invalid close reason" code = inf[:closeStatus] description = description
            self.io[:closeReason] = Closereason(CLOSE_REASON_NOT_PROVIDED, description)
        end
        if self.io[:closeTimeout] isa Timer && isopen(self.io[:closeTimeout])
            close(self.io[:closeTimeout])
        end
        isopen(self.io[:stream]) && close(self.io[:stream])
    end
end
# End receive data

"""
    close(ws::WebsocketConnection [, reasonCode::Int, description::String])
Closes a websocket connection.

Sends the close frame to the peer, waits for the response close frame,
or times out on `closeTimeout` before closing the underlying TCP stream.

Optional `reasonCode` must be a valid [rfc6455](https://tools.ietf.org/html/rfc6455#section-7.4) code,
suitable for sending over the wire.

Defaults: `1000` : "Normal connection closure"
"""
function Base.close(ws::WebsocketConnection, reasonCode::Int = CLOSE_REASON_NORMAL, description::String = "")
    closeConnection(ws, reasonCode, description)
end
"""
    broadcast(client::WebsocketConnection, data::Union{Array{UInt8,1}, String, Number})
[`send`](@ref) the given `data` as a message to all connected clients, except the given `client`.

In a SERVER context, to communicate with all clients on the SERVER, use [`emit`](@ref)

Only used if the `client` is in a SERVER context, otherwise NOOP.
"""
function Base.broadcast(client::WebsocketConnection, data::Union{Array{UInt8,1}, String, Number})
    client.clients === nothing && return
    for otherclient in client.clients
        otherclient.validate["valid"] !== true && continue
        client.id !== otherclient.id && send(otherclient, data)
    end
end

max_stars_repo_name: citkane/SimpleWebsockets

max_issues_repo_name: citkane/SimpleWebsockets

max_forks_repo_name: citkane/SimpleWebsockets


text:
@testset "map_partitions" begin

# test of map function

sc = SparkContext(master="local")

rdd = parallelize(sc, 1:10, n_split=2)
partitions = map_partitions(rdd, partition -> sum(partition))
values = collect(partitions)

@test values == [15, 40]

close(sc)

end


text:
function local_update!(model::XY, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    nbonds = numbonds(model)
    mbeta = -1.0/T

    @inbounds for site in 1:nsites
        center = model.spins[site]
        new_center = rand(rng)
        de = 0.0
        for (n,b) in neighbors(model, site)
            de += cospi(2(center - model.spins[n])) * Js[bondtype(model,b)]
            de -= cospi(2(new_center - model.spins[n])) * Js[bondtype(model,b)]
        end
        if rand(rng) < exp(mbeta*de)
            model.spins[site] = new_center
        end
    end

    return nothing
end

function SW_update!(model::XY, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    nbonds = numbonds(model)
    m2bJ = (-2/T).*Js
    m = 0.5*rand(rng)
    pspins = zeros(nsites)
    @inbounds for s in 1:nsites
        pspins[s] = cospi(2(model.spins[s]-m))
    end
    uf = UnionFind(nsites)
    @inbounds for bond in bonds(model)
        s1,s2 = source(bond), target(bond)
        bt = bondtype(bond)
        if rand(rng) < -expm1(m2bJ[bt]*pspins[s1]*pspins[s2])
            unify!(uf, s1,s2)
        end
    end
    nc = clusterize!(uf)
    toflips = rand(rng, Bool, nc)
    clustersize = zeros(Int, nc)
    @inbounds for site in 1:nsites
        id = clusterid(uf, site)
        s = model.spins[site]
        s = ifelse(toflips[id], s, 2m+0.5-s)
        s = mod(s, 1.0)
        clustersize[id] += 1
        model.spins[site] = s
    end

    return nothing
end

function Wolff_update!(model::XY, T::Real, Js::AbstractArray)
    rng = model.rng
    nsites = numsites(model)
    b2J = (2.0/T).*Js

    clustersize = 0
    st = Stack(Deque{Int}())

    m = 0.5*rand(rng)
    center = rand(rng,1:nsites)
    push!(st, center)

    model.spins[center] = mod(2m+0.5-model.spins[center], 1.0)

    @inbounds while !isempty(st)
        clustersize += 1
        center = pop!(st)
        cp = cospi(2(model.spins[center]-m))

        for (n,b) in neighbors(model, center)
            np = cospi(2(model.spins[n]-m))
            bt = bondtype(model,b)
            ## Note: center already flipped
            if rand(rng) < -expm1(b2J[bt]*cp*np)
                push!(st, n)
                model.spins[n] = mod(2m+0.5-model.spins[n], 1.0)
            end
        end
    end

    return nothing
end

@gen_convert_parameter(XY, ("T", 1, 1.0), ("J", numbondtypes, 1.0))


text:
TypeTuple{N} = NTuple{N, Type}

astupleoftypes(x::TypeTuple) = x
astupleoftypes(::Type{T}) where {T <: Tuple} = Tuple(T.parameters)

@inline foldlargs(op, x) = x
@inline foldlargs(op, x1, x2, xs...) =
    foldlargs(op,
              @return_if_reduced(op(x1, x2)),
              xs...)

@inline foldltupletype(op, T, ::Type{<:Tuple{}}) = T
@inline foldltupletype(op, T, ::Type{S}) where {S <: Tuple} =
    foldltupletype(op,
                   @return_if_reduced(op(T, Base.tuple_type_head(S))),
		   Base.tuple_type_tail(S))

@inline foldltupletype(op, ::Type{T}, ::Type{<:Tuple{}}) where T = T
@inline foldltupletype(op, ::Type{T}, ::Type{S}) where {T, S <: Tuple} =
    foldltupletype(op,
                   @return_if_reduced(op(T, Base.tuple_type_head(S))),
		   Base.tuple_type_tail(S))

asunion(T::Type{<:Tuple}) = foldltupletype((T, s) -> Union{T, s}, Union{}, T)

struct Padded{T, N}
    value::T
    pad::NTuple{N, UInt8}
end

Padded{T, N}(value::T) where {T, N} = Padded(value, zeropad(N))
zeropad(N) = ntuple(_ -> UInt8(0), N)
addpadding(N::Integer, ::Type{T}) where T = Padded{T, N}
addpadding(N::Integer, value::T) where T = Padded{T, N}(value)

Base.convert(::Type{P}, x::T) where {T, P <: Padded{T}} = P(x)

unpad(x) = x
unpad(x::Padded) = x.value
unpad(::Type{<:Padded{T}}) where T = T

paddedtype(::T) where T = paddedtype(T)
paddedtype(::Type{<:Padded{T}}) where T = T
paddedtype(::Type{T}) where T = T

# not sure relying on `sizeof` is safe; so:
sizeoftype(::Type{T}) where T = sizeof(T)
sizeoftype(::T) where T = sizeof(T)

ofsamesize(bigger::Type, smaller) =
    if sizeof(bigger) < sizeoftype(smaller)
        error("Target type is not big enough")
    elseif sizeof(bigger) == sizeoftype(smaller)
        smaller
    else
        addpadding(sizeof(bigger) - sizeoftype(smaller), smaller)
    end