text:
# # Example: Build-a-Yeast
#
# This is an example generated from this source
# file: [`build-a-yeast.jl`](@__REPO_ROOT_URL__examples/build-a-genome/build-a-yeast.jl)
# You are seeing the
#md # online documentation version. The corresponding notebook can be found
#md # here: [`build-a-yeast.ipynb`](@__NBVIEWER_ROOT_URL__examples/build-a-genome/build-a-yeast.ipynb)
#nb # jupyter notebook version. The corresponding online documentation page can
#nb # be found here: [`build-a-yeast.html`](https://bioinfologics.github.io/Pseudoseq.jl/dev/examples/build-a-genome/build-a-yeast),
# and the script can be found here: [`build-a-yeast.jl`](./build-a-yeast.jl)

# We are going to use Pseudoseq to create a simple set of fake genomes, based on
# chromosome 1 of the reference genome for yeast.
#
# 1. Load the sequence of chromosome 1 of the yeast reference from FASTA file.
# 2. Create blank chromosome blueprints for different genome designs.
# 3. Edit the blueprints to achive a:
#    a. Diploid chromosome
#    b. Triploid chromosome
#    c. Tetraploid chromosome

using FASTX, BioSequences, Pseudoseq

# ## Load the seed sequence
#
# First load the sequence from FASTA file. This uses tools from FASTX,
# which is a dependency of Pseudoseq and so the user should have it available,
# as julia's package manager would have installed any dependencies.
# The lines below open a fasta file, with a fasta reader, and load a single
# FASTA record, and get its sequence.

refseq = open(FASTA.Reader, "yeast-chr1.fasta") do rdr
    FASTA.sequence(LongSequence{DNAAlphabet{2}}, read(rdr))
end

reflen = length(refseq)

# This results in a 230,218nt long sequence called `refseq`.
#
# ## Designing a diploid
#
# First lets make a plan for a diploid genome with 1% heterozygosity.
# Pseudoseq has a lot of flexibilty in how you can define the heterozygosity,
# but for ease here I'm going to let it decide where the sites should be, and
# what the alleles at each site should be, by calling [`plan_het`](@ref),
# providing just the proportion of sites I want to be heterozygous, and the
# number of alleles I want at each of these heterozygous sites (2).

dploid = plan_chrom(reflen, 2)
dploidhet = plan_het(dploid, .01, 2)

# Now fabricate the FASTA, with [`fabricate`](@ref), providing a filename, and then a
# series of `blueprint => seed` pairs. In this case we just give 1 blueprint-seed
# pair.

fabricate("build-a-yeast-di.fasta", dploidhet => refseq)

# Thats it for the diploid genome.
#
# ## Designing a triploid
#
# Now let's make a triploid genome. This should also be simple enough. It's much
# like the process for making the diploid genome, but specifying 3 chromosome
# copies instead of 2.

trploid = plan_chrom(reflen, 3)

# When I add some heterozygous sites, just like with the diploid case, I'm going
# to simply ask for a proportion of heterozygous sites (and let Pseudoseq decide
# where they should occur. Like last time also I'm going to make all the sites have 2
# alleles, so for every heterozygous site, two of the three chromosome copies
# willl have the same base, and one of the copies will be different. There are
# other ways to define some heterozygous sites for example all 3 chromosome copies
# could have different bases

trploidhet = plan_het(trploid, .01, 2)

# Now fabricate the FASTA, with `fabricate`, providing a filename, and then a
# series of `blueprint => seed` pairs. In this case we just give 1 blueprint-seed
# pair.

fabricate("build-a-yeast-tri.fasta", trploidhet => refseq)

# That's it for the triploid genome.
#
# ## Designing a tetraploid genome
#
# Lets make a tetraploid genome, that consists of two diploid genomes, an A
# genome and a B genome.
# First I'm going to create a chromosome blueprint with 4 copies:

tetploid = plan_chrom(reflen, 4)

# We want to plan the heterozygosity in this tetraploid such that the
# number of heterozygous sites within the A and B genomes is ~1%. But we want
# ~3% of sites to be heterozygous between the A and B genomes.
#
# First we will say copies 1 + 2 in our blueprint are the A genome, and copies
# 3 + 4 are the B genome. Next consider heterozygous sites with 2 different alleles,
# A and B, 6 different patterns are possible:
#
# | Copy 1 | Copy 2 | Copy 3 | Copy 4 |
# |:------:|:------:|:------:|:------:|
# | A      | A      | B      | B      |
# | B      | B      | A      | A      |
# | A      | B      | A      | B      |
# | B      | A      | B      | A      |
# | A      | B      | B      | A      |
# | B      | A      | A      | B      |
#
# The first 2 patterns in the table constitute fixed mutations between
# genomes A and B,
# but not polymorphism within genomes A and B. This is because Copy 1 and 3 differ,
# as do copies 2 and 4, yet copies 1 and 2 are the same, as are copies 3 and 4.
# To achieve a ~3% divergence between genomes A and B, let's plan some
# divergent mutations by grouping the copy numbers together when planning some
# heterozygous sites:

thet = plan_het(tetploid, .03, [1, 2], [3, 4])

# By using the groups `[1, 2]` and `[3, 4]` we ensure chromosome copy 1 and 2
# will have the same allele, as will chromosomes 3 and 4

# Next we'll just add some of the other heterozygosity patterns to achieve the
# 1% internal heterozygosity desired. We'll use a heterozygosity pattern with
# two states.

thet = plan_het(thet, .01, [1, 3], [2, 4])

fabricate("build-a-yeast-tet.fasta", thet => refseq)


text:
function _error_on_nothing(idx)
    (idx===nothing) && throw(DomainError("x is not in the support of the mesh"))
end

function _get_nodes_coeffs_from_index(cell_idx::Int,i::Int,dq::DiscretisedFluidQueue)
    cellnodes = cell_nodes(dq)[:,cell_idx]
    coeff_idx = (N₋(dq) + (i-1)*n_bases_per_cell(dq) + (cell_idx-1)*n_bases_per_level(dq)) .+ (1:n_bases_per_cell(dq))
    return cellnodes, coeff_idx
end

function _get_coeff_index_pos(x::Float64,i::Int,dq::DiscretisedFluidQueue)
    cell_idx = findlast(x.>=dq.mesh.nodes)
    _error_on_nothing(cell_idx)
    cellnodes, coeff_idx = _get_nodes_coeffs_from_index(cell_idx,i,dq)
    return cell_idx, cellnodes, coeff_idx
end
function _get_coeff_index_neg(x::Float64,i::Int,dq::DiscretisedFluidQueue)
    cell_idx = findfirst(x.<=dq.mesh.nodes) - 1
    _error_on_nothing(cell_idx)
    cellnodes, coeff_idx = _get_nodes_coeffs_from_index(cell_idx,i,dq)
    return cell_idx, cellnodes, coeff_idx
end

function _get_point_mass_data_pos(i::Int,dq)
    cellnodes = dq.mesh.nodes[end]
    coeff_idx = N₋(dq) + total_n_bases(dq) + N₊(dq.model.S[1:i])
    return cellnodes, coeff_idx
end
function _get_point_mass_data_neg(i::Int,dq::DiscretisedFluidQueue)
    cellnodes = dq.mesh.nodes[1]
    coeff_idx = N₋(dq.model.S[1:i])
    return cellnodes, coeff_idx
end
_is_left_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue) =
    (x==dq.mesh.nodes[1])&&_has_left_boundary(dq.model.S,i)
_is_right_point_mass(x::Float64,i::Int,dq::DiscretisedFluidQueue) =
    (x==dq.mesh.nodes[end])&&_has_right_boundary(dq.model.S,i)

function _get_coeffs_index(x::Float64,i::Int,dq::DiscretisedFluidQueue)
    !(i∈phases(dq)) && throw(DomainError("phase i must be in the support of the model"))

    # find which inteval Dₗ,ᵢ x is in
    if _is_left_point_mass(x,i,dq)
        cell_idx = "point mass"
        cellnodes, coeff_idx = _get_point_mass_data_neg(i,dq)
    elseif _is_right_point_mass(x,i,dq)
        cell_idx = "point mass"
        cellnodes, coeff_idx = _get_point_mass_data_pos(i,dq)
    else # not a point mass
        if _has_right_boundary(dq.model.S,i)
            cell_idx, cellnodes, coeff_idx = _get_coeff_index_pos(x,i,dq)
        elseif _has_left_boundary(dq.model.S,i)
            cell_idx, cellnodes, coeff_idx = _get_coeff_index_neg(x,i,dq)
        end
    end

    return cell_idx, cellnodes, coeff_idx
end

"""
    pdf(d::SFMDistribution)

Return a function of two variables (x::Float64,i::Int) which is the
probability distribution function defined by `d`. Approximates the density function of
a fluid queue.
"""
function pdf(d::SFMDistribution)
    throw(DomainError("unknown <:Mesh}"))
end

function pdf(d::SFMDistribution{DGMesh{T}}) where T
    function f(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.dq.mesh
        fxi = 0.0
        if ((x<=mesh.nodes[1])||(x>=mesh.nodes[end]))
            #fxi = 0.0
        else
            cell_idx, cellnodes, coeff_idx = _get_coeffs_index(x,i,d.dq)
            coeffs = d.coeffs[coeff_idx]
            # if not a point mass, then reconstruct solution
            # if !(cell_idx=="point mass")
                if basis(mesh) == "legendre"
                    coeffs = legendre_to_lagrange(coeffs)
                else
                    V = vandermonde(n_bases_per_cell(mesh))
                    coeffs = coeffs
                end
                lp = lagrange_polynomials(cellnodes, x) # interpolating basis
                # transform to lagrange basis where each basis function integrates to 1,
                # i.e. divide the basis functions by their integrals.
                # their integrals are V.w*Δ/2
                basis_values = (lp./V.w)*2/Δ(d.dq,cell_idx)
                fxi = LinearAlgebra.dot(basis_values,coeffs)
            # else
            #     fxi = coeffs
            # end
        end
        return fxi
    end
    return f
end
"""
    pdf(d::SFMDistribution, x, i)

Evaluate pdf(d::SFMDistribution) at (x,i).
"""
pdf(d::SFMDistribution,x,i) =
    throw(DomainError("x must be Real, i must be Int/Array{Int,1}"))
pdf(d::SFMDistribution,x::Real,i::Int) = pdf(d)(x,i)

# abstract type ClosingOperator end

# struct UnnormalisedClosingOperator <: ClosingOperator
#     pdf::Function
#     cdf::Function
# end
# struct NaiveNormalisedClosingOperator <: ClosingOperator
#     cdf::Function
# end
# struct NormalisedClosingOperator <: ClosingOperator end

function unnormalised_closing_operator_pdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    #
    return x->pdf(a,me,x)
end
function unnormalised_closing_operator_cdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    return x->ccdf(a,me,x)
end

function naive_normalised_closing_operator_pdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    #
    return x->(pdf(a,me,x) + pdf(a,me,2.0-x))./cdf(a,me,2.0)
end
function naive_normalised_closing_operator_cdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    #
    return x->(ccdf(a,me,x)-ccdf(a,me,2.0-x))/cdf(a,me,2.0)
end

function normalised_closing_operator_pdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    # a [exp(Sx) + exp(S(2m-x))][I-exp(S2m)]^-1 s
    tmp = LinearAlgebra.I(size(me.S,1)) - exp(me.S*2.0)
    inv_factor = inv(tmp)

    return x->only(a*(exp(me.S*x)+exp(me.S*(2.0-x)))*inv_factor*me.s)
end
function normalised_closing_operator_cdf(a::AbstractArray{Float64,2},
    me::AbstractMatrixExponential)
    #
    tmp = LinearAlgebra.I(size(me.S,1)) - exp(me.S*2.0)
    inv_factor = inv(tmp)

    exp_factor(x) = LinearAlgebra.I(size(me.S,1))-(exp(me.S*x)-exp(me.S*(2.0-x))+exp(me.S*2.0))
    return x->only(1.0.-sum(a*exp_factor(x)*inv_factor,dims=2))
end

function pdf(d::SFMDistribution{FRAPMesh{T}},
    closing_operator::F=normalised_closing_operator_pdf) where {F<:Function, T}

    function f(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.dq.mesh
        fxi = 0.0
        if ((x<=mesh.nodes[1])||(x>=mesh.nodes[end]))
            # fxi = 0.0
        else
            cell_idx, ~, coeff_idx = _get_coeffs_index(x,i,d.dq)
            # if not a point mass, then reconstruct solution
            if _has_right_boundary(d.dq.model.S,i)
                yₖ₊₁ = mesh.nodes[cell_idx+1]
                to_go = (yₖ₊₁-x)./Δ(mesh,cell_idx)
            elseif _has_left_boundary(d.dq.model.S,i)
                yₖ = mesh.nodes[cell_idx]
                to_go = (x-yₖ)./Δ(mesh,cell_idx)
            end
            me = mesh.me
            fxi = closing_operator(Array(transpose(d.coeffs[coeff_idx])),me)(to_go)./Δ(mesh,cell_idx)
        end
        return fxi
    end
    return f
end

function pdf(d::SFMDistribution{FVMesh{T}}) where T
    function f(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.dq.mesh
        fxi = 0.0
        if ((x<=mesh.nodes[1])||(x>=mesh.nodes[end]))
            # fxi = 0.0
        else
            cell_idx, ~, coeff_idx = _get_coeffs_index(x,i,d.dq)
            # if not a point mass, then reconstruct solution
            ptsLHS = (membership(d.dq.model.S,i)<0) ? (Int(ceil(_order(mesh)/2))+1) : Int(ceil(_order(mesh)/2))
            if cell_idx-ptsLHS < 0
                nodesIdx = 1:_order(mesh)
                nodes = cell_nodes(mesh)[nodesIdx]
                poly_vals = lagrange_polynomials(nodes,x)
            elseif cell_idx-ptsLHS+_order(mesh) > n_bases_per_phase(mesh)
                nodesIdx = (n_bases_per_phase(mesh)-_order(mesh)+1):n_bases_per_phase(mesh)
                nodes = cell_nodes(mesh)[nodesIdx]
                poly_vals = lagrange_polynomials(nodes,x)
            else
                nodesIdx =  (cell_idx-ptsLHS) .+ (1:_order(mesh))
                poly_vals = lagrange_polynomials(cell_nodes(mesh)[nodesIdx],x)
            end
            coeff_idx = (N₋(d.dq) + i) .+ ((nodesIdx.-1).*n_phases(d.dq))
            coeffs = d.coeffs[coeff_idx]
            fxi = LinearAlgebra.dot(poly_vals,coeffs)
        end
        return fxi
    end
    return f
end

############
### CDFs ###
############
"""
    cdf(d::SFMDistribution)

Return a function of two variables (x::Real,i::Int) which is the
cumulative distribution distribution function defined by `d`. Approximates the distribution function of
a fluid queue.
"""
function cdf(d::SFMDistribution)
    throw(DomainError("unknown <:Mesh"))
end

"""

    _sum_cells_left(d::SFMDistribution, i::Int, cell_idx::Int, dq::DiscretisedFluidQueue)

Add up all the probability mass in phase `i` in the cells to the left of `cell_idx`.
"""
function _sum_cells_left(d::SFMDistribution, i::Int, cell_idx::Int)
    c = 0.0
    if basis(d.dq.mesh) == "legendre"
        for cell in 1:(cell_idx-1)
            # first legendre basis function =1 & has all the mass
            idx = (N₋(d.dq) + (i-1)*n_bases_per_cell(d.dq) + (cell-1)*n_bases_per_level(d.dq)) .+ 1
            c += d.coeffs[idx]
        end
    else
        for cell in 1:(cell_idx-1)
            idx = (N₋(d.dq) + (i-1)*n_bases_per_cell(d.dq) + (cell-1)*n_bases_per_level(d.dq)) .+ (1:n_bases_per_cell(d.dq))
            c += sum(d.coeffs[idx])
        end
    end
    return c
end

function cdf(d::SFMDistribution{DGMesh{T}}) where T
    # First, get the coeffs and project in to higher dimensional space
    # get coeffs without boundary masses
    coeffs = d.coeffs[N₋(d.dq)+1:end-N₊(d.dq)]
    coeffs = reshape(coeffs,n_bases_per_cell(d.dq),n_phases(d.dq),n_intervals(d.dq))
    # reweight basis to get it into the typicall lagrange interpolating basis
    v = vandermonde(n_bases_per_cell(d.dq))
    for i in 1:n_phases(d.dq)
        coeffs[:,i,:] = (coeffs[:,i,:] ./ Δ(d.dq)') ./ (v.w ./ 2.0)
    end
    # coeffs = lagrange_to_legendre(coeffs) # do this later to save computation power
    # project to higher space (in lagrange this is equiv. to adding a 0 coeff)
    # coeffs = cat(coeffs,zeros(1,n_intervals(d.dq),n_phases(d.dq)),dims=1)# do this later
    V = vandermonde(n_bases_per_cell(d.dq)+1)
    # differentiation operator in legendre basis of degree n+1
    # i.e. ϕ(x) a legendre basis for polynomials of degree n+1
    # p(x) = ϕ(x)a, where a is a vector of coeffs, a polynomial
    # p'(x) = ϕ'(x)a = ϕ(x)Ga
    G = V.inv*V.D
    # We have p'(x) = ϕ(x)a, we want to recover p(x).
    # p(x) is a polynomial of order n+1 so has a representation ϕ(x)b
    # and therefore p'(x)=ϕ'(x)b. To recover p(x) we want to find b.
    # We have, p'(x) = ϕ'(x)b = ϕ(x)a
    # multiplty by ϕᵀ(x) and integrate
    # ⟹ ∫ϕᵀ(x)ϕ'(x)b dx = ∫ϕᵀ(x)ϕ(x)a dx
    # ⟹ Gb=Ma=a (as M=I)
    # but G is not invertible (which makes sense as integrals are only
    # known up to an additive constant). If we have the additional condition that
    # p(-1)=0, then we can invert. i.e. ϕ(-1)b = 0.
    # First, evaluate ϕ(-1)
    ϕ₋₁ = V.V[1,:]
    # add condition to G
    G[end,:] .= ϕ₋₁
    # now G is invertible. We also need to add the constraint to the right-hand side
    # of the system. The system is
    # G⋅b = coeffs, with the additional constraint begin G[end,:]⋅b = 0, i.e.
    # we should set coeffs[:,end,:].=0, but this is already the case (the cat() command above)
    # coeffs[:,end,:] .= 0.0
    #
    # now solve it!
    # Ginv = G\I # = G^-1
    # also, transform back to lagrange basis for easy evaluation,
    # In summary multiply by V(n+1) G⁻¹(n+1) [I(n); zeros(n)] V⁻¹(n)
    # where the order of each operator is in ()
    # precompute transform
    transform = (V.V/G)*[LinearAlgebra.I(n_bases_per_cell(d.dq)); zeros(1,n_bases_per_cell(d.dq))]*v.inv
    integral_coeffs = zeros(n_bases_per_cell(d.dq)+1,n_phases(d.dq),n_intervals(d.dq))
    for i in 1:n_phases(d.dq)
        integral_coeffs[:,i,:] = (transform*coeffs[:,i,:]).*Δ(d.dq)'./2.0 # the factor at the end
        # is actually part of the integral operator G^-1. i.e. we should have Δ(d.dq,cell)G^-1
        # for each cell.
    end
    function F(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq.model)) && throw(DomainError("phase i must be in the support of the model"))
        mesh = d.dq.mesh
        Fxi = 0.0
        if (x<mesh.nodes[1])
            # Fxi = 0.0
        else # x>= mesh.nodes[1]
            # Fxi = 0.0
            # left pm
            if (x>=mesh.nodes[1])&&_has_left_boundary(d.dq.model.S,i)
                ~, left_pm_idx = _get_point_mass_data_neg(i,d.dq)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            if (x>mesh.nodes[1])&&(x<mesh.nodes[end])
                cell_idx, ~, ~ = _get_coeffs_index(x,i,d.dq)
                # add all mass from cells to the left
                Fxi += _sum_cells_left(d, i, cell_idx)

                # integrate up to x in the cell which contains x
                # i.e. evaluate p(x) from above
                cellnodes = gauss_lobatto_points(mesh.nodes[cell_idx],
                    mesh.nodes[cell_idx+1], n_bases_per_cell(d.dq)+1)

                basis_values = lagrange_polynomials(cellnodes, x) # interpolating basis
                Fxi += LinearAlgebra.dot(basis_values,integral_coeffs[:,i,cell_idx])
            elseif x>=mesh.nodes[end] # integrate over the whole space
                Fxi += _sum_cells_left(d, i, n_intervals(d.dq)+1)
                if _has_right_boundary(d.dq.model.S,i)
                    ~, right_pm_idx = _get_point_mass_data_pos(i,d.dq)
                    right_pm = d.coeffs[right_pm_idx]
                    Fxi += right_pm
                end
            end
        end
        return Fxi
    end
    return F
end
"""
    cdf(d::SFMDistribution, x, i)

Evaluate cdf(d::SFMDistribution) at (x,i).
"""
cdf(d::SFMDistribution,x,i) =
    throw(DomainError("x must be Float64/Int/Array{Float64/Int,1}, i must be Int/Array{Int,1}"))
cdf(d::SFMDistribution,x::Real,i::Int) = cdf(d)(x,i)
cdf(d::SFMDistribution,x::Int,i::Int) = cdf(d)(convert(Float64,x),i)

function cdf(d::SFMDistribution{FRAPMesh{T}},
    closing_operator::Function=normalised_closing_operator_cdf) where T

    function F(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.dq.mesh
        Fxi = 0.0
        if (x<mesh.nodes[1])
            # Fxi = 0.0
        else
            # Fxi = 0.0
            # left pm
            if _has_left_boundary(d.dq.model.S,i)
                ~, left_pm_idx = _get_point_mass_data_neg(i,d.dq)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            (x.>=mesh.nodes[end]) ? (xd=mesh.nodes[end]-sqrt(eps())) : xd = x
            cell_idx, ~, coeff_idx = _get_coeffs_index(xd,i,d.dq)
            coeffs = d.coeffs[coeff_idx]

            if !(cell_idx=="point mass")
                # add all mass from cells to the left
                Fxi += _sum_cells_left(d, i, cell_idx)

                # integrate up to x in the cell which contains x
                # me = build_me(cme_params[n_bases_per_cell(mesh)], mean = Δ(mesh)[cell_idx])
                me = mesh.me
                a = Array(coeffs')
                if _has_right_boundary(d.dq.model.S,i)
                    yₖ₊₁ = mesh.nodes[cell_idx+1]
                    mass = sum(a)
                    if mass > 0
                        a = a./mass
                        to_go = (yₖ₊₁-xd)/Δ(mesh,cell_idx)
                        Fxi += mass*closing_operator(a,me)(to_go)
                    end
                elseif _has_left_boundary(d.dq.model.S,i)
                    yₖ = mesh.nodes[cell_idx]
                    mass = sum(a)
                    if mass > 0
                        a = a./mass
                        to_go = (xd-yₖ)/Δ(mesh,cell_idx)
                        Fxi += mass*(1-closing_operator(a,me)(to_go)) # sum(a) - (ccdf(a,me,to_go) - ccdf(a,me,2.0-to_go))/cdf(a,me,2.0)
                    end
                end
            end
            if (x>=mesh.nodes[end])&&_has_right_boundary(d.dq.model.S,i)
                ~, right_pm_idx = _get_point_mass_data_pos(i,d.dq)
                right_pm = d.coeffs[right_pm_idx]
                Fxi += right_pm
            end
        end
        return Fxi
    end
    return F
end

function _sum_cells_left(d::SFMDistribution{FVMesh{T}}, i::Int, cell_idx::Int) where T
    c = 0
    for cell in 1:(cell_idx-1)
        idx = N₋(d.dq) + i + (cell-1)*n_bases_per_level(d.dq)
        c += d.coeffs[idx]*Δ(d.dq,cell)
    end
    return c
end

function cdf(d::SFMDistribution{FVMesh{T}}) where T
    function F(x::Real,i::Int) # the PDF
        # check phase is in support
        !(i∈phases(d.dq)) && throw(DomainError("phase i must be in the support of the model"))
        # if x is not in the support return 0.0
        mesh = d.dq.mesh
        Fxi = 0.0
        if (x<mesh.nodes[1])
            # Fxi = 0.0
        else
            # Fxi = 0.0
            # left pm
            if (x>=mesh.nodes[1])&&_has_left_boundary(d.dq.model.S,i)
                ~, left_pm_idx = _get_point_mass_data_neg(i,d.dq)
                left_pm = d.coeffs[left_pm_idx]
                Fxi += left_pm
            end
            # integral over density
            (x.>=mesh.nodes[end]) ? (xd=mesh.nodes[end]-sqrt(eps())) : xd = x
            cell_idx, ~, ~ = _get_coeffs_index(xd,i,d.dq)
            # if not a point mass, then reconstruct solution
            if !(cell_idx=="point mass")
                # add all mass from cells to the left
                Fxi += _sum_cells_left(d, i, cell_idx)

                # integrate up to x in the cell which contains x
                temp_pdf(y) = pdf(d)(y,i)
                quad = gauss_lobatto_quadrature(temp_pdf,mesh.nodes[cell_idx],xd,_order(mesh))
                Fxi += quad
            end
            if (x>=mesh.nodes[end])&&_has_right_boundary(d.dq.model.S,i)
                ~, right_pm_idx = _get_point_mass_data_pos(i,d.dq)
                right_pm = d.coeffs[right_pm_idx]
                Fxi += right_pm
            end
        end
        return Fxi
    end
    return F
end

function cell_probs(d::SFMDistribution)
    function p(x::Real,i::Int)
        _x_in_bounds = (x>d.dq.mesh.nodes[1])&&(x<d.dq.mesh.nodes[end])
        if _x_in_bounds
            ~, ~, coeff_idx = _get_coeff_index_pos(x,i,d.dq)
            return sum(d.coeffs[coeff_idx])
        else
            return 0.0
        end
    end
    return p
end

function cell_probs(d::SFMDistribution{FVMesh{T}}) where T
    function p(x::Real,i::Int)
        _x_in_bounds = (x>d.dq.mesh.nodes[1])&&(x<d.dq.mesh.nodes[end])
        if _x_in_bounds
            cell_idx, ~, coeff_idx = _get_coeff_index_pos(x,i,d.dq)
            cell_average = d.coeffs[coeff_idx]
            return only(cell_average*Δ(d.dq,cell_idx))
        else
            return 0.0
        end
    end
    return p
end

ext: jl

max_stars_repo_name: angus-lewis/StochasticFluidQueues

max_issues_repo_name: angus-lewis/StochasticFluidQueues

max_forks_repo_name: angus-lewis/StochasticFluidQueues


text:
module Cuddle

using Reexport

for p in Pkg.available()
  eval(:(@reexport using $(symbol(p))))
end

end # module


text:
using Test
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using HDGElasticity

function allapprox(v1,v2)
    return all(v1 .≈ v2)
end

function allapprox(v1,v2,atol)
    @assert length(v1) == length(v2)
    return all([isapprox(v1[i],v2[i],atol=atol) for i = 1:length(v1)])
end

HL = reshape([0.0],1,1)
sfunc(x) = [0.5*x[1]*(x[1]+1)]
vfunc(x) = [(1.0-x[2])*(1.0+x[2])*x[1]*(x[1]-1.0)*0.5]
linemap = HDGElasticity.LineMap([-1.,-1.],[-1.,+1.])
squad = tensor_product_quadrature(1,3)
comp = [1.0]
icomp = ones(length(squad))'
D = reshape([1.0],1,1)

HDGElasticity.HLop!(HL,sfunc,vfunc,squad,linemap,comp,comp,[D],1.)
@test allapprox(HL,[2/15])

fill!(HL,0.0)
HDGElasticity.HLop!(HL,sfunc,vfunc,squad,linemap,comp,comp,[D],0.5)
@test allapprox(HL,0.5*[2/15])

fill!(HL,0.0)
HDGElasticity.HLop!(HL,sfunc,vfunc,squad,linemap,-comp,comp,[D],0.5)
@test allapprox(HL,-0.5*[2/15])

fill!(HL,0.0)
HDGElasticity.HLop!(HL,sfunc,vfunc,squad,linemap,comp,-comp,[D],0.5)
@test allapprox(HL,0.5*[2/15])

imap = InterpolatingPolynomial(2,1,1)
coeffs = [-1.0,-1.,-1.,1.]
update!(imap,coeffs)
fill!(HL,0.0)
scale = ones(length(squad))
HDGElasticity.HLop!(HL,sfunc,vfunc,squad,imap,icomp,icomp,D,scale)
@test allapprox(HL,[2/15])

fill!(HL,0.0)
HDGElasticity.HLop!(HL,sfunc,vfunc,squad,imap,icomp,icomp,D,0.5*scale)
@test allapprox(HL,0.5*[2/15])

sbasis = TensorProductBasis(1,1)
vbasis = TensorProductBasis(2,1)
linemap = HDGElasticity.LineMap([-1.,-1.],[-1.,1.])
components = [0.0,1.0]
normals = [-1.,0.0]
icomponents = repeat(components,inner=(1,length(squad)))
inormals = repeat(normals,inner=(1,length(squad)))

Dhalf = sqrt(HDGElasticity.plane_strain_voigt_hooke_matrix_2d(1.0,2.))
Ek = HDGElasticity.vec_to_symm_mat_converter(2)
ED = [Ek[k]'*Dhalf for k = 1:2]

cellmap = HDGElasticity.CellMap([0.,0.],[2.,1.])
facejacobian = HDGElasticity.face_determinant_jacobian(cellmap)
HL = zeros(4,12)
HDGElasticity.HLop!(HL,sbasis,vbasis,squad,linemap,normals,components,ED,facejacobian[4])

HLl = HDGElasticity.HLop(sbasis,vbasis,squad,linemap,normals,components,Dhalf,facejacobian[4])
@test size(HLl) == (4,12)

iscale = HDGElasticity.scale_area(cellmap,inormals)
HLi = HDGElasticity.HLop(sbasis,vbasis,squad,imap,inormals,icomponents,Dhalf,iscale)
@test size(HLi) == (4,12)
@test allapprox(HLl,HLi)

HLl2 = HDGElasticity.HLop(sbasis,vbasis,squad,linemap,normals,components,Dhalf,0.5facejacobian[4])
@test allapprox(0.5HLl,HLl2)

HU = reshape([0.0],1,1)
HDGElasticity.HUop!(HU,sfunc,vfunc,squad,linemap,comp,1.,1)
@test allapprox(HU,[2/15])

fill!(HU,0.0)
HDGElasticity.HUop!(HU,sfunc,vfunc,squad,linemap,comp,0.5,1)
@test allapprox(HU,0.5*[2/15])

fill!(HU,0.0)
HDGElasticity.HUop!(HU,sfunc,vfunc,squad,linemap,-comp,0.5,1)
@test allapprox(HU,0.5*[2/15])

fill!(HU,0.0)
coeffs = [-1.,-1.,-1.,1.]
update!(imap,coeffs)
scale = ones(length(squad))
HDGElasticity.HUop!(HU,sfunc,vfunc,squad,imap,icomp,scale,1)
@test allapprox(HU,[2/15])

fill!(HU,0.0)
HDGElasticity.HUop!(HU,sfunc,vfunc,squad,imap,icomp,0.5scale,1)
@test allapprox(HU,0.5*[2/15])

HUl = HDGElasticity.HUop(sbasis,vbasis,squad,linemap,components,facejacobian[4])
@test size(HUl) == (4,8)

scale = HDGElasticity.scale_area(cellmap,inormals)
HUi = HDGElasticity.HUop(sbasis,vbasis,squad,imap,icomponents,scale)
@test size(HUi) == (4,8)
@test allapprox(HUl,HUi)

HUl2 = HDGElasticity.HUop(sbasis,vbasis,squad,linemap,components,2*facejacobian[4])
@test size(HUl) == (4,8)
@test allapprox(2HUl,HUl2)

hybloc1 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,vbasis,
    squad,linemap,normals,components,Dhalf,1.,facejacobian[4])
@test size(hybloc1) == (4,20)

hybloc2 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,vbasis,
    squad,linemap,normals,components,Dhalf,0.1,facejacobian[4])
@test allapprox(0.1hybloc1,hybloc2)

iscale = HDGElasticity.scale_area(cellmap,inormals)
hybloc1 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,vbasis,
    squad,imap,inormals,icomponents,Dhalf,1.,iscale)
@test size(hybloc1) == (4,20)

hybloc2 = HDGElasticity.hybrid_local_operator_traction_components(sbasis,vbasis,
    squad,imap,inormals,icomponents,Dhalf,0.2,iscale)
@test allapprox(0.2hybloc1,hybloc2)


text:
using Makie, GeometryTypes, Colors
using AbstractPlotting: slider!, playbutton
using Observables

cd(@__DIR__)
using Pkg
Pkg.pkg"add NRRD"
using FileIO
data = load(joinpath(homedir(), "desktop", "brain.nrrd"))
psps = [rand(50, 50, 50) for i in 1:50]

mini, maxi = mapreduce(extrema, (a, b)-> (min(a[1], b[1]), max(a[2], b[2])), psps)
psps_n = map(psps) do vol
    Float32.((vol .- mini) ./ (maxi - mini))
end

volume100 = last(psps_n)
scene = Scene()
g_area = map(scene.px_area) do a
    IRect(0, 0, widths(a)[1], 50)
end
legend_w = 110
vol_area = map(scene.px_area) do a
    w, h = widths(a)
    w = w - legend_w
    IRect(0, 50, w / 2, h - 50)
end
heat_area = map(scene.px_area) do a
    w, h = widths(a)
    w = w - legend_w
    IRect(w / 2, 50, w / 2, h - 50)
end
gui = Scene(scene, g_area)
campixel!(gui)
theme(gui)[:plot] = NT(raw = true)

scene3d = Scene(scene, vol_area)
hscene = Scene(scene, heat_area)
cam3d!(scene3d)
r = range(-1, stop = 1, length = size(volume100, 1))
c = contour!(
    scene3d, r, r, r, volume100,
    alpha = 0.4, levels = 8, colorrange = (0.0, 1.0)
)[end]

rotate!(hscene, Vec3f0(0, 1, 0), -0.5pi)

volume = c[4]
planes = (:xy, :xz, :yz)
sliders = ntuple(3) do i
    s = slider!(gui, 1:size(volume[], i))[end]
    move!(s, size(volume[], i) ÷ 2)
    idx = s[:value]; plane = planes[i]
    indices = ntuple(3) do j
        planes[j] == plane ? 1 : (:)
    end
    hmap = contour!(
        hscene, r, r, volume[][indices...],
        raw = true, colorrange = (0.0, 1.0), fillrange = true,
        interpolate = true, linewidth = 0.1
    )[end]
    foreach(idx, volume) do _idx, vol
        idx = (i in (1, 2)) ? (size(vol, i) - _idx) + 1 : _idx
        transform!(hmap, (plane, r[_idx]))
        indices = ntuple(3) do j
            planes[j] == plane ? idx : (:)
        end
        if checkbounds(Bool, vol, indices...)
            hmap[3][] = view(vol, indices...)
        end
    end
    s
end
t = nothing
b = playbutton(gui, 1:length(psps_n)) do idx
    println(idx)
    move!(t, idx)
end
t = slider!(gui, 1:length(psps_n))[end];
foreach(t[:value]) do idx
    if checkbounds(Bool, psps_n, idx)
        c[4][] = psps_n[idx]
    end
end

b2 = button(gui, "3d/2d"; dimensions = (60, 40)) do clicks
    if iseven(clicks)
        cam3d!(hscene)
    else
        cam2d!(hscene)
        update_cam!(hscene, FRect(-1.25, -1.25, 2.5, 2.5))
    end
end

AbstractPlotting.vbox(gui.plots)

l = lines!(hscene, FRect(-1.01, -1.01, 2.02, 2.02), color = :black, linewidth = 2, raw = true)
transform!(l, (:yz, -1.0))
pixelarea(scene)[] = IRect(0,0,1500, 1000)
scene


text:
import LinearAlgebra.diagm

## Kalman filter for multivariate linear Gaussian models

struct KalmanMVOut{d}
  predictionMeans::Vector{MVector{d, Float64}}
  predictionVariances::Vector{MMatrix{d, d, Float64}}
  filteringMeans::Vector{MVector{d, Float64}}
  filteringVariances::Vector{MMatrix{d, d, Float64}}
  smoothingMeans::Vector{MVector{d, Float64}}
  smoothingVariances::Vector{MMatrix{d, d, Float64}}
  logZhats::Vector{Float64}
end

function logmvnormpdf(y::StaticVector{d}, μ::StaticVector{d},
  Σ::StaticMatrix{d,d}) where d
  v::SVector{d,Float64} = y-μ
  lnc::Float64 = - 0.5 * d * log(2 * π) - 0.5 * logdet(Σ)
  return lnc - 0.5 * dot(v, Σ \ v)
end

function kalmanMV(theta::MVLGTheta, ys::Vector{SVector{d,Float64}}) where d
  n = length(ys)
  A = theta.A
  C = theta.C
  Q = theta.Q
  R = theta.R
  x0 = theta.x0
  v0 = theta.v0
  predictionMeans = Vector{MVector{d,Float64}}(undef, n)
  predictionVariances = Vector{MMatrix{d,d,Float64}}(undef, n)
  filteringMeans = Vector{MVector{d,Float64}}(undef, n)
  filteringVariances = Vector{MMatrix{d,d,Float64}}(undef, n)
  logZhats = Vector{Float64}(undef, n)
  mutt1 = MVector{d,Float64}(undef)
  mutt = MVector{d,Float64}(undef)
  sigmatt1 = MMatrix{d,d,Float64}(undef)
  sigmatt = MMatrix{d,d,Float64}(undef)
  lZ = 0.0
  for p = 1:n
    if p == 1
      mutt1 = x0
      sigmatt1 = diagm(0 => v0)
    else
      mutt1 = A * mutt
      sigmatt1 = A * sigmatt * A' + Q
    end
    predictionMeans[p] = mutt1
    predictionVariances[p] = sigmatt1
    lZ += logmvnormpdf(ys[p], C*mutt1, C*sigmatt1*C' + R)
    logZhats[p] = lZ
    K = sigmatt1 * C' * inv(C * sigmatt1 * C' + R)
    mutt = mutt1 + K * (ys[p] - C * mutt1)
    sigmatt = sigmatt1 - K * C * sigmatt1
    filteringMeans[p] = mutt
    filteringVariances[p] = sigmatt
  end
  smoothingMeans = Vector{MVector{d,Float64}}(undef, n)
  smoothingVariances = Vector{MMatrix{d,d,Float64}}(undef, n)
  smoothingMeans[n] = filteringMeans[n]
  smoothingVariances[n] = filteringVariances[n]
  for p = n:-1:2
    J = filteringVariances[p-1] * A' * inv(predictionVariances[p])
    smoothingMeans[p-1] = filteringMeans[p-1] +
      J * (smoothingMeans[p] - predictionMeans[p])
    smoothingVariances[p-1] = filteringVariances[p-1] +
      J * (smoothingVariances[p] - predictionVariances[p]) * J'
  end

  return KalmanMVOut(predictionMeans, predictionVariances, filteringMeans,
        filteringVariances, smoothingMeans, smoothingVariances, logZhats)
end


text:
using LyonPassiveEar
using Test

@testset "soscascade" begin
    EPS = 0.01

    input = [-0.50, -0.50, -0.70, -1.10, -1.70]

    coeffs = reshape(0:9, 2, 5) ./ 100

    state = zeros(2, 2)

    output, state = LyonPassiveEar.soscascade(input, coeffs, state)

    exp_state = [-0.07 -0.06; -0.00 -0.00]
    @test all(abs.(state .- exp_state) .< EPS)

    exp_output = [0.00 -0.01 -0.03 -0.03 -0.05;
                  0.00 -0.00 -0.00 -0.00 -0.00]'

    @test all(abs.(output .- exp_output) .< EPS)
end

# this doesn't work yet but doesn't seem to be used anyway in the passive ear function?

# @testset "sosfilters" begin
#     EPS = 0.01
#
#     coeffs = reshape(0:9, 2, 5) ./ 100
#     state = zeros(2, 2)
#     agcOut = reshape(0:9, 2, 5)' .* 0.1
#
#     output = LyonPassiveEar.sosfilters(agcOut, coeffs, state)
#
#     exp_state = [0.04 0.03; 0.06 0.04]'
#     @test all(abs.(state .- exp_state) .< EPS)
#
#     exp_output = [0.00 0.00 0.00 0.02 0.03;
#                   0.00 0.01 0.02 0.04 0.05]'
#     @test all(abs.(output .- exp_output) .< EPS)
# end

@testset "sosfilters 2" begin

    EPS = 1e-5

    decFilt = [0.          0.          0.02356786 -1.69296345  0.71653131]

    inp = [0.         0.        ;
              0.         0.        ;
              0.         0.00116701;
              0.         0.        ;
              0.         0.        ]

    decState = [0.03989792 -0.01688643;
                0.02095108 -0.00886735;
                0.01244632 -0.00526648;
                0.00594191 -0.00251233;
                0.00233343 -0.00097628]

    output2, decState2 = LyonPassiveEar.sosfilters(inp, decFilt, decState)

    outexp = [0.03989792 0.05065929;
              0.02095108 0.02660206;
              0.01244632 0.01580469;
              0.00594191 0.00754711;
              0.00233343 0.00297413]

    decStateexp = [0.05717621 -0.03629896;
                   0.03002421 -0.01906121;
                   0.01783859 -0.01129705;
                   0.00851941 -0.00540774;
                   0.00336311 -0.00213106]

    @test all(abs.(output2 .- outexp) .< EPS)
    @test all(abs.(decState2 .- decStateexp) .< EPS)

end

@testset "agc" begin

    EPS = 0.01
    nSamples = 5
    nChannels = 2
    nStages = 3

    output = (0:nChannels*nSamples-1) .* 0.1
    agcParams = fill(0.5, 2 * nStages)
    agcState = zeros(nChannels * nStages)
    agcOut = zeros(nChannels * nSamples)

    LyonPassiveEar.agc(output, nChannels, nSamples, nStages, agcParams, agcState, agcOut)

    exp_agcState = vec([0.64 0.41 0.30;
                        0.66 0.42 0.30])
    @test all(abs.(agcState .- exp_agcState) .< EPS)

    exp_agcOut = vec([0.00 0.20 0.19 0.15 0.15;
                      0.10 0.22 0.19 0.15 0.15])

    @test all(abs.(agcOut .- exp_agcOut) .< EPS)

end

@testset "design_lyon_filters" begin
    result = LyonPassiveEar.design_lyon_filters(4000.0)

    # from python version
    @test all(abs.(result[2] .- [1897.10751154, 1831.00638006, 1766.69348626, 1704.10601949, 1643.18285418,
        1583.86449022, 1526.0929948, 1469.81194586, 1414.96637697, 1361.50272363, 1309.368771,
        1258.51360287, 1208.88755196, 1160.44215137, 1113.1300873, 1066.90515279, 1021.72220264,
        977.53710928, 934.3067197, 891.98881327, 850.54206056, 809.92598293, 770.10091301,
        731.02795597, 692.66895152, 654.98643664, 617.943609, 581.50429104, 545.63289457,
        510.29438609, 475.45425253, 441.07846755, 407.13345831, 373.5860727, 340.40354693,
        307.55347355, 275.00376981, 242.72264629, 210.67857589, 178.84026303, 147.17661309,
        115.656702, 84.2497461]) .< 1e-6)
end

@testset "lyon_passive_ear" begin
    expected_output = [0.         0.00502543 0.01080264;
                       0.         0.00384371 0.00826242;
                       0.         0.00213897 0.00459793]

    result = lyon_passive_ear([1.0, 0, 0, 0, 0, 0], sample_rate = 400, decimation_factor = 2)

    @test all(abs.(expected_output .- result) .< 1e-5)

    expected_output_2 = [0.         0.00734996 0.01580469;
                         0.         0.00350625 0.00754711;
                         0.         0.00136727 0.00297413]

    result_2 = lyon_passive_ear([1.0, 0, 0, 0, 0, 0], sample_rate = 400, decimation_factor = 2, differ = false)

    @test all(abs.(expected_output_2 .- result_2) .< 1e-5)
end


text:
module PSplines

import LinearAlgebra: diagind, I, Diagonal
import Kronecker: kronecker
import ForwardDiff

# buildmatrices.jl
export D1matrix, D2matrix, D3matrix, Ω1matrix, Ω2matrix, Ω3matrix, Wmatrix
export PD2matrices2dim
include("buildMatrices.jl")

# types.jl
export BSpline, BSpline2dim, knots, vertices, buildKnots
include("types.jl")

#utils.jl
export equiSpacedPosition, buildSystemMatrix, buildSystemIndVect
export evalBSplineBasis, evalBSpline
export evalBSplineBasis2dim, evalBSpline2dim
export evalBSplineBasisD1, evalBSplineD1
export evalBSplineBasis2dimD1x, evalBSpline2dimD1x, evalBSplineBasis2dimD1y, evalBSpline2dimD1y
export expandBSpline, expandBSpline2dim
export nVerticesNeededLeftRight, extendedKnots
include("utils.jl")

#promotion.jl
export promote_rule, convert
include("promotion.jl")
end


text:
using PokeDB
using Base.Test

# mergesorted
it1 = mergesorted(1:10, 5:15)
@assert collect(it1) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

it2 = mergesorted(5:15, 1:10)
@assert collect(it2) == [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]

it3 = mergesorted([:a, :c, :d, :f], [:b, :c, :d, :e])
@assert collect(it3) == [:a, :b, :c, :d, :e, :f]

it4 = mergesorted([], ["aaa", "bbb", "ccc"])
@assert collect(it4) == ["aaa", "bbb", "ccc"]

it5 = mergesorted(["aaa", "bbb", "ccc"], [])
@assert collect(it5) == ["aaa", "bbb", "ccc"]

it6 = mergesorted([], [])
@assert collect(it6) == []


text:
using DiskBackedDicts
const DBD = DiskBackedDicts
using Test

function test_dict_interface(d_candidate, d_test)
    @assert isempty(d_candidate)
    @assert !isempty(d_test)

    @test isempty(d_candidate)
    @test isempty(keys(d_candidate))
    @test isempty(values(d_candidate))
    @test length(d_candidate) == 0

    k, v = first(d_test)
    @test !haskey(d_candidate, k)
    @test v === get(d_candidate, k, v)
    d_candidate[k] = v
    @test !isempty(d_candidate)
    @test !isempty(keys(d_candidate))
    @test !isempty(values(d_candidate))
    @test haskey(d_candidate, k)
    @test d_candidate[k] == v
    @test d_candidate == delete!(d_candidate, k)
    @test_throws KeyError d_candidate[k]
    @test d_candidate == delete!(d_candidate, k)
    @test isempty(d_candidate)
    @test v === get!(d_candidate, k, v)
    get!(error, d_candidate, k)
    get(error, d_candidate, k)
    delete!(d_candidate, k)
    @test v === get(() -> v, d_candidate, k)
    @test !haskey(d_candidate, k)
    @test v === get!(() -> v, d_candidate, k)
    @test d_candidate[k] === v

    merge!(d_candidate, d_test)
    @test length(d_candidate) == length(d_test)
    @test length(d_candidate) == length(keys(d_candidate))
    @test length(d_candidate) == length(values(d_candidate))
    associative_elements_equal(d_candidate, d_test)

    @test !isempty(d_candidate)
    @test isempty(empty!(d_candidate))
end

function associative_elements_equal(d1,d2)
    @test length(d1) == length(d2)
    for (k,v) in d1
        @test haskey(d2, k)
        @test d2[k] == v
    end
end

struct MyString
    a::String
end
Base.:(==)(s1::MyString, s2::MyString) = s1.a == s2.a
struct MyInt
    b::Int
end
struct MyPair
    s::MyString
    t::MyInt
end
struct MyContainer{T}
    inner::T
end

@testset "AbstractDict interface" begin
    test_dicts = []

    d_test = Dict("a"=>1, "b"=>2)
    push!(test_dicts, d_test)

    t = MyInt(1)
    s = MyString("s")
    st = MyPair(s,t)
    d_test = Dict("a"=>1, "b"=>2, "s" => s, st => t)
    push!(test_dicts, d_test)

    d_test = Dict(MyString("a") => MyInt(2), MyString("") => MyInt(0))
    push!(test_dicts, d_test)

    for d_test in test_dicts
        K = eltype(keys(d_test))
        V = eltype(values(d_test))
        candidates = [
            DiskBackedDict{K,V}(tempname()*".jld2"),
            DBD.JLD2BlobDict{K,V}(tempname()*".jld2"),
            DBD.FullyCachedDict(DBD.JLD2BlobDict{K,V}(tempname()*".jld2")),
            DBD.CachedDict{K,V}(Dict{K,V}(), Dict{K,V}()),
            Dict{K,V}(),
            DBD.JLD2FilesDict{K,V}(tempname()*".jld2"),
            DBD.CachedDict{K,V}(Dict{K,V}(), DBD.JLD2FilesDict{K,V}(tempname()*".jld2")),
        ]
        for d_candidate in candidates
            test_dict_interface(d_candidate, d_test)
        end
    end

end

@testset "DiskBackedDict" begin
    d = DiskBackedDict(tempname()*".jld2")
    @test d isa DiskBackedDict{Any, Any}

    d = @inferred DiskBackedDict{MyString,MyInt}(tempname()*".jld2")
    @test d isa AbstractDict{MyString,MyInt}

    path = tempname()*".jld2"
    d = DiskBackedDict{Int, String}(path)
    d[1] = "one"
    @test_throws MethodError DiskBackedDict{MyString,MyInt}(path)

    @assert ispath(path)
    d2 = DiskBackedDict{Int, String}(path)
    @test typeof(d2) == typeof(d)
    @test d2[1] == "one"
    @test length(d2) == 1

    d3 = DiskBackedDict(path)
    @test_broken typeof(d3) == typeof(d2)
    @test d3[1] == "one"
    @test length(d3) == 1
end

function test_long_term_persistence(path, d::Dict{K,V}) where {K,V}
    if !ispath(path)
        d_save = DiskBackedDict{K,V}(path)
        @info("$path does not exist, create it with content $d.")
        merge!(d_save, d)
    else
        @info("Using $path from previous run.")
    end
    d_loaded = DiskBackedDict{K,V}(path)
    associative_elements_equal(d, d_loaded)
end

assetpath(args...) = joinpath(@__DIR__, "assets", args...)
@testset "long term persistence" begin
    for (s, d) ∈ [
                  "empty.jld2" => Dict(),
                  "1.jld2"     => Dict(1 => 1),
                  "mixed.jld2"=> Dict(1 => "1", :two => [1,2]),
                  "custom_types.jld2" => Dict(MyInt(1) => MyString("two")),
                  "parametric_type.jld2" => Dict(MyContainer(1) => MyContainer(MyInt(1))),
                 ]

        path = assetpath(s)
        test_long_term_persistence(path, d)
    end
end

@testset "merge!" begin
    dst_ref = Dict(1=>1, 2=>2)
    src_ref = Dict(2=>2.5, 3=>3)
    res_ref = Dict(1=>1, 2=>2.5, 3=>3)
    function _replace!(dst, src)
        empty!(dst)
        merge!(dst, src)
    end
    srcs = [DBD.JLD2BlobDict{Any,Any}(tempname() * ".jld2"),
            DBD.FullyCachedDict(Dict(),Dict()),
            DBD.DiskBackedDict(tempname() * ".jld2"),
            Dict(),
           ]
    dsts = [DBD.JLD2BlobDict{Any,Any}(tempname() * ".jld2"),
            DBD.FullyCachedDict(Dict(),Dict()),
            DBD.DiskBackedDict(tempname() * ".jld2"),
            Dict(),
           ]
    for src in srcs
        for dst in dsts
            _replace!(src, src_ref)
            _replace!(dst, dst_ref)
            associative_elements_equal(src, src_ref)
            associative_elements_equal(dst, dst_ref)
            res = merge!(dst, src)
            @test res === dst
            @test typeof(res) == typeof(dst)
            associative_elements_equal(res, res_ref)
        end
    end
end


text:
function isValid(input::String)
    checkNine(n) = n > 9 ? n - 9 : n

    input = filter(x -> x != ' ' ,collect(input))

    length(input) <= 1 && return false

    parsedArr = try
        map(x -> parse(Int,x), input)
    catch
        return false
    end

    flag = true
    result = 0
    for i in reverse(parsedArr)
        if flag
            result += i
        else
            result += checkNine(2 * i)
        end
        flag = !flag
    end

    0 == (result % 10)
end

max_stars_repo_name: ccqpein/Arithmetic-Exercises

max_issues_repo_name: ccqpein/Arithmetic-Exercises

max_forks_repo_name: ccqpein/Arithmetic-Exercises


text:
using Plots

function lif!(v,I)
    v = v + 0.2*(I-v-70)
    return v
end

function step_v()
    vhist = []
    v = -50
    I = 0.0
    for i in 1:30
        if i%4==0
            I += 2
        end
        v = lif!(v,I)
        @show v
        push!(vhist, v)
    end
    return vhist
end
# %%
vhist = step_v()
# %%

plot(vhist)

#%%
v= -50
#%%
v = lif!(v,0.0)

# %%

using Plots

include("../../../Modules/Networks/net_arch.jl")
include("../../../Modules/Poisson/PoissonStateSpike.jl")

params = YAML.load(open("Julia/CartPole/LIF_base/cfglif4820.yml"))
n = Network([1,1], params)

# %%
n.v[1] = 0
n.s.e[1][1] = 0.7
n.v[2] = 0

times = [5,7,9,15,17]

vh = []
x = []

for t in 1:30
    spiked = step_LIF!(n, [(1,0.4)], false)
    push!(x,t)
    push!(vh, n.v[1])
end

vh2 = []
n.v[1] = 0

for t in 1:30
    spiked = step_LIF!(n, [(1,0.8)], false)
    push!(vh2, n.v[1])
end

plot(x,vh, line = :blue, label = "Membrane Potential I=0.4")
plot!(x,vh2, line = :red, label = "Membrane Potential I=0.8")
hline!([1.0], line = (3,:dash, :green), label = "Threshold Potential")
xlabel!("Time (ms)")
ylabel!("Potential")

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
# Test Stan.jl examples.

using Test

examples = [

  "Examples/ARM/Ch03/Kid/kidscore.jl",
  "Examples/Bernoulli/bernoulli.jl",
  "Examples/Binomial/binomial.jl",
  "Examples/Binormal/binormal.jl",
  "Examples/Dyes/dyes.jl",
  "Examples/EightSchools/schools8.jl",

  "Examples_Stan_Methods/Diagnose/diagnose.jl",
  "Examples_Stan_Methods/Generate_Quantities/generate_quantities.jl",
  "Examples_Stan_Methods/Optimize/optimize.jl",
  "Examples_Stan_Methods/Parse_and_Interpolate/parse.jl",
  "Examples_Stan_Methods/Variational/variational.jl",

  "Examples_Test_Cases/Diagnostics/diagnostics.jl",
  "Examples_Test_Cases/InitThetaDict/init_dict.jl",
  "Examples_Test_Cases/InitThetaDictArray/init_dict_array.jl",
  "Examples_Test_Cases/InitThetaFile/init_file.jl",
  "Examples_Test_Cases/NamedArray/namedarray.jl",
  "Examples_Test_Cases/ScalarObs/scalar.jl",
  "Examples_Test_Cases/ZeroLengthArray/zerolengtharray.jl",

]

@testset "Stan.jl v6.0" begin

  for example in examples
      #println("\n  * $(example) *")
      @testset "$(example)" begin
        include(example)
      end
      #println("\n * $(example) done! *")
  end

end


text:
using Base.Iterators: partition
using Flux
using Flux.Optimise: update!
using Flux: logitbinarycrossentropy, binarycrossentropy, mae
using Images
using MLDatasets
using Statistics
using Parameters: @with_kw
using Random
using Printf
using CUDAapi
using Zygote
if has_cuda()		# Check if CUDA is available
    @info "CUDA is on"
    import CuArrays		# If CUDA is available, import CuArrays
    CuArrays.allowscalar(false)
end

@with_kw struct HyperParams
    batch_size::Int = 128
    latent_dim::Int = 10
    epochs::Int = 50
    verbose_freq::Int = 1000
    output_x::Int = 6        # No. of sample images to concatenate along x-axis
    output_y::Int = 6        # No. of sample images to concatenate along y-axis
    lr_dscr::Float64 = 0.0002
    lr_gen::Float64 = 0.0002
end

function reparameterization(mu, logvar, args)
    std = exp.(logvar / 2) |> gpu
    sampled_z = Float32.(randn(size(mu)...)) |> gpu
    return sampled_z .* std .+ mu
end

struct encoder
    m
    mu
    logvar
end

function encoder(args)
    m =  Chain(x->reshape(x, :, size(x, 4)), Dense(784, 1024), x-> relu.(x),Dense(1024,1024), BatchNorm(1024), x->relu.(x)) |> gpu
    mu = Dense(1024, args.latent_dim) |> gpu
    logvar = Dense(1024, args.latent_dim) |>gpu
    return encoder(m, mu, logvar)
end

function decoder(args)
    return Chain(Dense(args.latent_dim, 1024), x-> relu.(x), Dense(1024,1024),BatchNorm(1024), x-> relu.(x),
                 Dense(1024, 784, tanh), x-> reshape(x, 28,28,1,size(x,2))) |> gpu
end

function discriminator(args)
    return Chain(Dense(args.latent_dim, 1024), x-> leakyrelu.(x, 0.2f0), Dense(1024, 256),x-> leakyrelu.(x, 0.2f0),Dense(256, 1, sigmoid)) |> gpu
end

function load_data(hparams)
    # Load MNIST dataset
    images, labels = MLDatasets.MNIST.traindata(Float32)
    # Normalize to [-1, 1] and convert it to WHCN
    image_tensor = permutedims(reshape(@.(2f0 * images - 1f0), 28, 28, 1, :), (2, 1, 3, 4))
    # Partition into batches
    data = [image_tensor[:,:,:,r] |> gpu for r in partition(1:60000, hparams.batch_size)]
    return data
end

# Loss functions
function discr_loss(real_output, fake_output)
    real_loss = mean(binarycrossentropy.(real_output, 1f0))
    fake_loss = mean(binarycrossentropy.(fake_output, 0f0))
    return (real_loss + fake_loss)/2
end

generator_loss(fake_output) = mean(binarycrossentropy.(fake_output, 1f0))

function pixelwise_loss(real_output, decoded_output)
    return mean(Flux.mae.(real_output, decoded_output))
end

function train_discr(discr, fake_data, opt_discr, hparams)
    noise = randn!(similar(fake_data, (hparams.latent_dim, hparams.batch_size))) |> gpu
    ps = Flux.params(discr)
    loss, back = Zygote.pullback(ps) do
                      discr_loss(discr(noise), discr(fake_data))
    end
    grads = back(1f0)
    update!(opt_discr, ps, grads)
    return loss
end
Zygote.@nograd train_discr

function train_gan(encoder,decoder, discr, original_data, opt_en_dec, opt_discr, hparams)
    loss = Dict()
    ps = Flux.params(encoder.m, encoder.mu, encoder.logvar, decoder)
    loss["gen"], back = Zygote.pullback(ps) do
                          en_m = encoder.m(original_data)
                          encoded_img = reparameterization(encoder.mu(en_m), encoder.logvar(en_m), hparams)
                          loss["discr"] = train_discr(discr, encoded_img, opt_discr, hparams)
                          decoded_img = decoder(encoded_img)
                          0.001f0*generator_loss(discr(encoded_img)) + 0.999f0*pixelwise_loss(original_data, decoded_img)
    end
    grads = back(1f0)
    update!(opt_en_dec, ps, grads)
    return loss
end

function create_output_image(dec, fixed_noise, hparams)
    @eval Flux.istraining() = false
    fake_images = @. cpu(dec(fixed_noise))
    @eval Flux.istraining() = true
    image_array = dropdims(reduce(vcat, reduce.(hcat, partition(fake_images, hparams.output_y))); dims=(3, 4))
    image_array = @. Gray(image_array + 1f0) / 2f0
    return image_array
end

function train()
    hparams = HyperParams()

    data = load_data(hparams)

    fixed_noise = [randn(hparams.latent_dim, 1) |> gpu for _=1:hparams.output_x*hparams.output_y]

    # Discriminator
    dscr = discriminator(hparams)
    # Generator
    enc =  encoder(hparams)

    dec = decoder(hparams)

    # Optimizers
    opt_dscr = ADAM(hparams.lr_dscr)
    opt_en_dec = ADAM(hparams.lr_gen)

    isdir("output")||mkdir("output")

    # Training
    train_steps = 0
    for ep in 1:hparams.epochs
        @info "Epoch $ep"
        for x in data
            # Update discriminator and generator
            loss = train_gan(enc, dec, dscr, x, opt_en_dec, opt_dscr, hparams)

            if train_steps % hparams.verbose_freq == 0
                @info("Train step $(train_steps), Discriminator loss = $(loss["discr"]), Generator loss = $(loss["gen"])")
                # Save generated fake image
                output_image = create_output_image(dec, fixed_noise , hparams)
                save(@sprintf("output/aae_steps_%06d.png", train_steps), output_image)
            end
            train_steps += 1
        end
    end

    output_image = create_output_image(dec, fixed_noise, hparams)
    save(@sprintf("output/aae_steps_%06d.png", train_steps), output_image)
end

cd(@__DIR__)
train()


text:
"""
TrackedFiles

Convenience type to keep track of files to watch.
"""
const TrackedFiles = Dict{Pair{String, String}, Float64}

"""
Prepare the output directory `PATHS[:pub]`.

## Argument

* `clear=true`: removes the content of the output directory if it exists to
                start from a blank slate
"""
function prepare_output_dir(clear::Bool=true)::Nothing
    if FD_ENV[:STRUCTURE] < v"0.2"
        # if required to start from a blank slate -> remove the output dir
        (clear & isdir(PATHS[:pub])) && rm(PATHS[:pub], recursive=true)
        # create the output dir and the css dir if necessary
        !isdir(path(:pub)) && mkdir(path(:pub))
        !isdir(path(:css)) && mkdir(path(:css))
    else
        (clear & isdir(path(:site))) && rm(path(:site), recursive=true)
        !isdir(path(:site)) && mkdir(path(:site))
    end
    return nothing
end

"""
$(SIGNATURES)

Given a file path split in `(base, file)`, form the output path (where the
output file will be written/copied).
"""
function form_output_path(base::AS, file::AS, case::Symbol)
    # .md -> .html for md pages:
    case == :md && (file = change_ext(file))
    if FD_ENV[:STRUCTURE] < v"0.2"
        outbase = _out_path(base)
    else
        outbase = _out_path2(base)
        if case in (:md, :html)
            # file is index.html or 404.html --> keep the path
            # file is page.html  --> .../page/index.html
            fname = splitext(file)[1]
            if fname ∉ ("index", "404")
                file = joinpath(fname, "index.html")
            end
        end
    end
    outpath = joinpath(outbase, file)
    outdir  = splitdir(outpath)[1]
    isdir(outdir) || mkpath(outdir)
    return outpath
end

# NOTE: LEGACY way of getting the target path
function _out_path(base::String)::String
    if startswith(base, PATHS[:src_css])
        f_out_path = replace(base, PATHS[:src_css] => PATHS[:css])
        !ispath(f_out_path) && mkpath(f_out_path)
        return f_out_path
    end
    len_in = lastindex(joinpath(PATHS[:src], ""))
    length(base) <= len_in && return PATHS[:folder]
    dpath = base[nextind(base, len_in):end]
    # construct the out path
    f_out_path = joinpath(PATHS[:folder], dpath)
    f_out_path = replace(f_out_path, r"([^a-zA-Z\d\s_:])pages" => s"\1pub")
    # if it doesn't exist, make the path
    !ispath(f_out_path) && mkpath(f_out_path)
    return f_out_path
end

function _out_path2(base::String)::String
    if startswith(base, path(:assets)) ||
       startswith(base, path(:css))    ||
       startswith(base, path(:layout)) ||
       startswith(base, path(:libs))   ||
       startswith(base, path(:literate))

       # add a closing separator to folder path
       rbase   = joinpath(path(:folder), "") * "_"
       outpath = base[nextind(base, length(rbase)):end]
       outpath = joinpath(path(:site), outpath)
   else
       # path is not a 'special folder'
       outpath = replace(base, path(:folder) => path(:site))
   end
    # if it doesn't exist, make the path
    !ispath(outpath) && mkpath(outpath)
    return outpath
end

_access(p)   = isa(p, Regex) ? p.pattern : p
_isempty(p)  = isempty(_access(p))
_endswith(p) = endswith(_access(p), '/')

"""
$(SIGNATURES)

Update the dictionaries referring to input files and their time of last change.
The variable `verb` propagates verbosity.
"""
function scan_input_dir!(args...; kw...)
    to_ignore = union(IGNORE_FILES, globvar("ignore"))
    # remove empty patterns or strings
    filter!(!_isempty, to_ignore)
    # differentiate between files and dirs pattern
    dir_indicator = [_endswith(c) for c in to_ignore]
    # ignore "/"
    d2i = filter!(d -> length(d) > 1, to_ignore[dir_indicator])
    f2i = to_ignore[.!dir_indicator]

    if FD_ENV[:STRUCTURE] < v"0.2"
        return _scan_input_dir!(args...; files2ignore=f2i, dirs2ignore=d2i)
    end
    return _scan_input_dir2!(args...; files2ignore=f2i, dirs2ignore=d2i, kw...)
end

function _scan_input_dir!(other_files::TrackedFiles,
                          infra_files::TrackedFiles,
                          md_files::TrackedFiles,
                          html_files::TrackedFiles,
                          literate_files::TrackedFiles,
                          verb::Bool=false;
                          files2ignore=String[],
                          dirs2ignore=String[])::Nothing
    # top level files (src/*)
    for file ∈ readdir(PATHS[:src])
        fpath = joinpath(PATHS[:src], file)
        isfile(fpath) || continue
        # skip if it has to be ignored
        should_ignore(fpath, files2ignore, dirs2ignore) && continue
        fname, fext = splitext(file)
        fpair = (PATHS[:src] => file)
        if file == "config.md"
            add_if_new_file!(infra_files, fpair, verb)
        elseif fext == ".md"
            add_if_new_file!(md_files, fpair, verb)
        else
            add_if_new_file!(html_files, fpair, verb)
        end
    end
    # pages files (src/pages/*)
    for (root, _, files) ∈ walkdir(PATHS[:src_pages])
        for file ∈ files
            fpath = joinpath(root, file)
            isfile(fpath) || continue
            # skip if it has to be ignored
            should_ignore(fpath, files2ignore, dirs2ignore) && continue
            fname, fext = splitext(file)
            fpair = (root => file)
            if fext == ".md"
                add_if_new_file!(md_files, fpair, verb)
            elseif fext == ".html"
                add_if_new_file!(html_files, fpair, verb)
            else
                add_if_new_file!(other_files, fpair, verb)
            end
        end
    end
    # infastructure files (src/_css/* and src/_html_parts/*)
    for d ∈ (:src_css, :src_html), (root, _, files) ∈ walkdir(PATHS[d])
        for file ∈ files
            isfile(joinpath(root, file)) || continue
            fname, fext = splitext(file)
            # skipping files that are not of the type INFRA_FILES
            fext ∉ INFRA_FILES && continue
            add_if_new_file!(infra_files, root=>file, verb)
        end
    end
    # literate script files if any, note that the folder may not exist
    if isdir(PATHS[:literate])
        for (root, _, files) ∈ walkdir(PATHS[:literate])
            for file ∈ files
                isfile(joinpath(root, file)) || continue
                fname, fext = splitext(file)
                # skipping files that are not script file
                fext != ".jl" && continue
                add_if_new_file!(literate_files, root=>file, verb)
            end
        end
    end
    return nothing
end

function _scan_input_dir2!(other_files::TrackedFiles,
                           infra_files::TrackedFiles,
                           md_pages::TrackedFiles,
                           html_pages::TrackedFiles,
                           literate_scripts::TrackedFiles,
                           verb::Bool=false;
                           in_loop::Bool=false,
                           files2ignore::Vector=String[],
                           dirs2ignore::Vector=String[])::Nothing
    # go over all files in the website folder
    for (root, _, files) ∈ walkdir(path(:folder))
        for file in files
            # assemble full path (root is an absolute path)
            fpath = joinpath(root, file)
            fpair = root => file
            fext  = splitext(file)[2]

            opts = (fpair, verb, in_loop)

            # early skips
            !isfile(fpath) && continue
            should_ignore(fpath, files2ignore, dirs2ignore) && continue

            # skip over `__site` folder, `.git` and `.github` folder
            startswith(fpath, path(:site)) && continue
            startswith(fpath, joinpath(path(:folder), ".git")) && continue
            # skip over toml files
            fext == ".toml" && continue

            # assets file --> other
            if startswith(fpath, path(:assets))
                add_if_new_file!(other_files, opts...)
            # infra_files
            elseif startswith(fpath, path(:css))    ||
                   startswith(fpath, path(:layout)) ||
                   startswith(fpath, path(:libs))
                add_if_new_file!(infra_files, opts...)
            # literate_files
            elseif startswith(fpath, path(:literate))
                # ignore files that are not script files
                fext == ".jl" || continue
                add_if_new_file!(literate_scripts, opts...)
            else
                if file == "config.md"
                    add_if_new_file!(infra_files, opts...)
                elseif fext == ".md"
                    add_if_new_file!(md_pages, opts...)
                elseif fext ∈ (".html", ".htm")
                    add_if_new_file!(html_pages, opts...)
                else
                    add_if_new_file!(other_files, opts...)
                end
            end
        end
    end
    return nothing
end

"""
$SIGNATURES

Helper function, if `fpair` is not referenced in the dictionary (new file) add
the entry to the dictionary with the time of last modification as val.
"""
function add_if_new_file!(dict::TrackedFiles, fpair::Pair{String,String},
                          verb::Bool, in_loop::Bool=false)::Nothing
    haskey(dict, fpair) && return nothing
    # it's a new file
    verb && println("tracking new file '$(fpair.second)'.")
    # save it's modification time, set to zero if it's a new file in a loop
    # to force its processing in FS2
    dict[fpair] = ifelse(in_loop, 0, mtime(joinpath(fpair...)))
    return nothing
end

"""
$SIGNATURES

Check if a file path should be ignored. This is a helper function for the
`scan_input_dir` functions. A bit of a similar principle to `match_url`. The
file argument is an absolute path, the list would be a list of relative paths.
Rules:

* `''` or `'/'`  -> ignore
* `'path/fname'` -> ignore exactly that
* `'path/dir/'`  -> ignore everything starting with `path/dir/`
"""
function should_ignore(fpath::AS, files2ignore::Vector,
                       dirs2ignore::Vector)::Bool
    # fpath is necessarily an absolute path so can strip the folder part
    if FD_ENV[:STRUCTURE] < v"0.2"
        fpath = fpath[length(path(:src))+length(PATH_SEP)+1:end]
    else
        fpath = fpath[length(path(:folder))+length(PATH_SEP)+1:end]
    end
    if any(c -> c isa Regex ? match(c, fpath) !== nothing : c == fpath,
               files2ignore)
        return true
    end
    flag = findfirst(c -> startswith(fpath, c), dirs2ignore)
    isnothing(flag) || return true
    return false
end


text:
module IMEXExperiments

# Some explicit RK integrators

function rk2(u::S, f::Function, t::T, dt::T) where {S<:Number, T<:Real}
    t1 = t
    u1 = u
    f1 = f(t1, u1)

    t2 = t + dt
    u2 = u + dt/2 * f1
    f2 = f(t2, u2)

    t + dt, u + dt * f2
end

function rk_ssp2_222(u::S, f::Function, t::T, dt::T) where
        {S<:Number, T<:Real}
    t1 = t
    u1 = u
    f1 = f(t1, u1)

    t2 = t + dt
    u2 = u + dt * f1
    f2 = f(t2, u2)

    t + dt, u + dt/2 * f1 + dt/2 * f2
end

function rk_ssp3_332(u::S, f::Function, t::T, dt::T) where
        {S<:Number, T<:Real}
    t1 = t
    u1 = u
    f1 = f(t1, u1)

    t2 = t + dt
    u2 = u + dt * f1
    f2 = f(t2, u2)

    t3 = t + dt/2
    u3 = u + dt/4 * f1 + dt/4 * f2
    f3 = f(t3, u3)

    t + dt, u + dt/6 * f1 + dt/6 * f2 + dt*2/3 * f3
end

# Some IMEX RK integrators

# Notation: SSP_k(s, σ, p) where
#     k: order of SSP scheme
#     s: number of stages of the implicit scheme
#     σ: number of stages of the explicit scheme
#     p: order of the IMEX scheme

function imex_ssp2_222(u::S, f::Function, r::Function, t::T, dt::T) where
        {S<:Number, T<:Real}
    γ = 1 - 1/√T(2)

    t1 = t
    t1′ = t + γ * dt
    # u1 = u + γ * dt * r(t1′, u1)
    u1, r1 = r(t1′, u, γ * dt)
    f1 = f(t1, u1)

    t2 = t + dt
    t2′ = t + (1 - γ) * dt
    # u2 = u + dt * f1 + (1 - 2γ) * dt * r1 + (1 - γ) * dt * r(t + γ * dt, u2)
    u2, r2 = r(t2′, u + dt * f1 + (1 - 2γ) * dt * r1, γ * dt)
    f2 = f(t2, u2)

    t + dt, u + dt/2 * f1 + dt/2 * f2 + dt/2 * r1 + dt/2 * r2
end

# A coupled harmonic oscillator as example system

struct State{T<:Number} <: Number
    x1::T
    v1::T
    x2::T
    v2::T
end

function map1(f::Function, s::State{T})::State{T} where {T}
    State{T}(f(s.x1), f(s.v1), f(s.x2), f(s.v2))
end

function map2(f::Function, s1::State{T}, s2::State{T})::State{T} where {T}
    State{T}(f(s1.x1, s2.x1), f(s1.v1, s2.v1), f(s1.x2, s2.x2), f(s1.v2, s2.v2))
end

import Base: +, -, *, /
+(s::State{T}) where {T<:Number} = map1(+, s)
-(s::State{T}) where {T<:Number} = map1(-, s)

+(s::State{T}, x::T) where {T<:Number} = map1(y->y+x, s)
-(s::State{T}, x::T) where {T<:Number} = map1(y->y-x, s)
*(s::State{T}, x::T) where {T<:Real} = map1(y->y*x, s)
/(s::State{T}, x::T) where {T<:Real} = map1(y->y/x, s)

+(x::T, s::State{T}) where {T<:Number} = map1(y->x+y, s)
-(x::T, s::State{T}) where {T<:Number} = map1(y->x-y, s)
*(x::T, s::State{T}) where {T<:Real} = map1(y->x*y, s)
/(x::T, s::State{T}) where {T<:Real} = map1(y->x/y, s)

+(s1::State{T}, s2::State{T}) where {T<:Number} = map2(+, s1, s2)
-(s1::State{T}, s2::State{T}) where {T<:Number} = map2(-, s1, s2)
*(s1::State{T}, s2::State{T}) where {T<:Real} = map2(*, s1, s2)
/(s1::State{T}, s2::State{T}) where {T<:Real} = map2(/, s1, s2)

function output(n::Integer, t::T, s::State{T}) where {T}
    global outfile
    e = energy(s)
    println("$n $t $e")
    println(outfile, "$n $t $(s.x1) $(s.v1) $(s.x2) $(s.v2) $e")
end

const omega1 = Float64(1)
const omega2 = Float64(2)
const alpha1 = Float64(0)
const alpha2 = Float64(1000)

function initial(::Type{T}) where {T<:Real}
    t::T = 0
    x1 = -1
    v1 = 0
    x2 = -1
    v2 = 0
    t::T, State{T}(x1, v1, x2, v2)
end

function rhs(t::T, s::State{T}) where {T<:Real}
    x1 = s.v1                   # + alpha1 * (s.x2 - s.x1)
    v1 = - omega1^2 * s.x1
    x2 = s.v2                   # + alpha2 * (s.x1 - s.x2)
    v2 = - omega2^2 * s.x2
    State{T}(x1, v1, x2, v2)
end

function irhs(t::T, s0::State{T}, a::T) where {T<:Real}
    # u1 = u0 + a * r(t1, u1)
    # u1, r1 = r(t1, u0, a)

    # x1 = alpha1 * (s.x2 - s.x1)
    # v1 = 0
    # x2 = alpha2 * (s.x1 - s.x2)
    # v2 = 0
    # r = State{T}(x1, v1, x2, v2)
    # s = s0 + a * r

    offset = a * (alpha2 * s0.x1 + alpha1 * s0.x2)
    denom = 1 + a * (alpha1 + alpha2)

    x1 = (s0.x1 + offset) / denom
    v1 = s0.v1
    x2 = (s0.x2 + offset) / denom
    v2 = s0.v2
    s = State{T}(x1, v1, x2, v2)

    x1rhs = - alpha1 * (s0.x1 - s0.x2) / denom
    v1rhs = 0
    x2rhs = - alpha2 * (s0.x2 - s0.x1) / denom
    v2rhs = 0
    r = State{T}(x1rhs, v1rhs, x2rhs, v2rhs)

    s, r
end

function energy(s::State{T})::T where {T<:Real}
    e1 = T(1)/2 * omega1 * s.x1^2 + T(1)/2 / omega1 * s.v1^2
    e2 = T(1)/2 * omega2 * s.x2^2 + T(1)/2 / omega2 * s.v2^2
    e1 + e2
end

# Driver

function evolve()
    global outfile = open("harmonic.asc", "w")
    println(outfile, "# Iteration time x1 v1 x2 v2 e")

    odeint = rk2
    # odeint = rk_ssp2_222
    # odeint = rk_ssp3_332
    odeimex = imex_ssp2_222

    t::Float64, s::State{Float64} = initial(Float64)
    output(0, t, s)
    niters = 1000
    dt::Float64 = 0.05
    for n in 1:niters
        # t, s = odeint(s, rhs, t, dt)
        t, s = odeimex(s, rhs, irhs, t, dt)
        output(n, t, s)
    end

    close(outfile)
end

end


text:
### A Pluto.jl notebook ###
# v0.17.3

using Markdown
using InteractiveUtils

# ╔═╡ e39a8330-8947-4227-83ea-bcd344601159
# hideall

using Genie, Genie.Router, Genie.Renderer.Html, Genie.Requests

# ╔═╡ e28a7f22-319f-11ec-21e2-45db28a06407
md"""
# Handling file uploads

Genie has built-in support for working with file uploads. The collection of uploaded files (as `POST` variables) can be accessed through the `Requests.filespayload` method. Or, we can retrieve the data corresponding to a given file form input by using `Requests.filespayload(key)` -- where `key` is the name of the file input in the form.

In the following snippet we configure two routes in the root of the app (`/`): the first route, handling `GET` requests, displays an upload form. The second route, handling `POST` requests, processes the uploads, generating a file from the uploaded data, saving it, and displaying the file stats.

### Example
"""

# ╔═╡ 6e9844b4-bd77-4e0b-9908-a5afb2a1dd54
md"""

```julia
using Genie, Genie.Router, Genie.Renderer.Html, Genie.Requests

form = \"""
<form action="/" method="POST" enctype="multipart/form-data">
  <input type="file" name="yourfile" /><br/>
  <input type="submit" value="Submit" />
</form>
\"""

route("/") do
  html(form)
end

route("/", method = POST) do
  if infilespayload(:yourfile)
    write(filespayload(:yourfile))

    stat(filename(filespayload(:yourfile)))
  else
    "No file uploaded"
  end
end

up()
```
"""

# ╔═╡ e9a246a3-c6e7-4c0e-b220-61d21707cf57
# hideall

form = """
<form action="/" method="POST" enctype="multipart/form-data">
  <input type="file" name="yourfile" /><br/>
  <input type="submit" value="Submit" />
</form>
""";

# ╔═╡ 747ac8a7-ba1d-48be-98ec-e66756ca9b02
# hideall

route("/") do
  Genie.Renderer.Html.html(form)
end;

# ╔═╡ 7d3decbf-7276-448b-949a-ceea2ffd0b4b
# hideall

route("/", method = POST) do
  if infilespayload(:yourfile)
    write(filespayload(:yourfile))

    stat(filename(filespayload(:yourfile)))
  else
    "No file uploaded"
  end
end;

# ╔═╡ 968c2c12-b564-432e-b879-70162b070462
# hideall

Genie.Router.routes();

# ╔═╡ 6e511e93-a104-481b-bb8d-35771768a1a2
# hideall

up(async=true);

# ╔═╡ 7a7bcced-c077-4531-9176-23ddc88c6886
# hideall

down(); #shuting down the server

# ╔═╡ 43a94071-cb2b-4815-a25b-614e85bd2a0f
md"""
Upon uploading a file and submitting the form, our app will display the file's stats.
"""

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
Genie = "c43c736e-a2d1-11e8-161f-af95117fbd1e"

[compat]
Genie = "~3.2.0"
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

[[ArgParse]]
deps = ["Logging", "TextWrap"]
git-tree-sha1 = "3102bce13da501c9104df33549f511cd25264d7d"
uuid = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
version = "1.1.4"

[[ArgTools]]
uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"

[[Artifacts]]
uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"

[[Base64]]
uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"

[[CSTParser]]
deps = ["Tokenize"]
git-tree-sha1 = "b2667530e42347b10c10ba6623cfebc09ac5c7b6"
uuid = "00ebfdb7-1f24-5e51-bd34-a7502290713f"
version = "3.2.4"

[[CodeTracking]]
deps = ["InteractiveUtils", "UUIDs"]
git-tree-sha1 = "9aa8a5ebb6b5bf469a7e0e2b5202cf6f8c291104"
uuid = "da1fd8a2-8d9e-5ec2-8556-3022fb5608a2"
version = "1.0.6"

[[CommonMark]]
deps = ["Crayons", "JSON", "URIs"]
git-tree-sha1 = "393ac9df4eb085c2ab12005fc496dae2e1da344e"
uuid = "a80b9123-70ca-4bc0-993e-6e3bcb318db6"
version = "0.8.3"

[[Compat]]
deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "SHA", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
git-tree-sha1 = "31d0151f5716b655421d9d75b7fa74cc4e744df2"
uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
version = "3.39.0"

[[Crayons]]
git-tree-sha1 = "3f71217b538d7aaee0b69ab47d9b7724ca8afa0d"
uuid = "a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"
version = "4.0.4"

[[DataStructures]]
deps = ["Compat", "InteractiveUtils", "OrderedCollections"]
git-tree-sha1 = "7d9d316f04214f7efdbb6398d545446e246eff02"
uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
version = "0.18.10"

[[Dates]]
deps = ["Printf"]
uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"

[[DelimitedFiles]]
deps = ["Mmap"]
uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"

[[Distributed]]
deps = ["Random", "Serialization", "Sockets"]
uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"

[[Downloads]]
deps = ["ArgTools", "LibCURL", "NetworkOptions"]
uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"

[[EzXML]]
deps = ["Printf", "XML2_jll"]
git-tree-sha1 = "0fa3b52a04a4e210aeb1626def9c90df3ae65268"
uuid = "8f5d6c58-4d21-5cfd-889c-e3ad7ee6a615"
version = "1.1.0"

[[FilePathsBase]]
deps = ["Dates", "Mmap", "Printf", "Test", "UUIDs"]
git-tree-sha1 = "7fb0eaac190a7a68a56d2407a6beff1142daf844"
uuid = "48062228-2e41-5def-b9a4-89aafe57970f"
version = "0.9.12"

[[FileWatching]]
uuid = "7b1f6079-737a-58dc-b8bc-7a2ca5c1b5ee"

[[GMP_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"

[[Genie]]
deps = ["ArgParse", "Dates", "Distributed", "EzXML", "FilePathsBase", "HTTP", "HttpCommon", "Inflector", "JSON3", "JuliaFormatter", "Logging", "Markdown", "MbedTLS", "Millboard", "Nettle", "OrderedCollections", "Pkg", "REPL", "Reexport", "Revise", "SHA", "Serialization", "Sockets", "UUIDs", "Unicode", "YAML"]
git-tree-sha1 = "02fde56046c4ebec4771e84c258bc982f702c96c"
uuid = "c43c736e-a2d1-11e8-161f-af95117fbd1e"
version = "3.2.0"

[[HTTP]]
deps = ["Base64", "Dates", "IniFile", "Logging", "MbedTLS", "NetworkOptions", "Sockets", "URIs"]
git-tree-sha1 = "14eece7a3308b4d8be910e265c724a6ba51a9798"
uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
version = "0.9.16"

[[HttpCommon]]
deps = ["Dates", "Nullables", "Test", "URIParser"]
git-tree-sha1 = "46313284237aa6ca67a6bce6d6fbd323d19cff59"
uuid = "77172c1b-203f-54ac-aa54-3f1198fe9f90"
version = "0.5.0"

[[Inflector]]
deps = ["Unicode"]
git-tree-sha1 = "8555b54ddf27806b070ce1d1cf623e1feb13750c"
uuid = "6d011eab-0732-4556-8808-e463c76bf3b6"
version = "1.0.1"

[[IniFile]]
deps = ["Test"]
git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
version = "0.5.0"

[[InteractiveUtils]]
deps = ["Markdown"]
uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"

[[JLLWrappers]]
deps = ["Preferences"]
git-tree-sha1 = "642a199af8b68253517b80bd3bfd17eb4e84df6e"
uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
version = "1.3.0"

[[JSON]]
deps = ["Dates", "Mmap", "Parsers", "Unicode"]
git-tree-sha1 = "8076680b162ada2a031f707ac7b4953e30667a37"
uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
version = "0.21.2"

[[JSON3]]
deps = ["Dates", "Mmap", "Parsers", "StructTypes", "UUIDs"]
git-tree-sha1 = "7d58534ffb62cd947950b3aa9b993e63307a6125"
uuid = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
version = "1.9.2"

[[JuliaFormatter]]
deps = ["CSTParser", "CommonMark", "DataStructures", "Pkg", "Tokenize"]
git-tree-sha1 = "10c95cebcfa37c1f510a726c90886db4745e1238"
uuid = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
version = "0.15.11"

[[JuliaInterpreter]]
deps = ["CodeTracking", "InteractiveUtils", "Random", "UUIDs"]
git-tree-sha1 = "e273807f38074f033d94207a201e6e827d8417db"
uuid = "aa1ae85d-cabe-5617-a682-6adf51b2e16a"
version = "0.8.21"

[[LibCURL]]
deps = ["LibCURL_jll", "MozillaCACerts_jll"]
uuid = "b27032c2-a3e7-50c8-80cd-2d36dbcbfd21"

[[LibCURL_jll]]
deps = ["Artifacts", "LibSSH2_jll", "Libdl", "MbedTLS_jll", "Zlib_jll", "nghttp2_jll"]
uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"

[[LibGit2]]
deps = ["Base64", "NetworkOptions", "Printf", "SHA"]
uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"

[[LibSSH2_jll]]
deps = ["Artifacts", "Libdl", "MbedTLS_jll"]
uuid = "29816b5a-b9ab-546f-933c-edad1886dfa8"

[[Libdl]]
uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"

[[Libiconv_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "42b62845d70a619f063a7da093d995ec8e15e778"
uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
version = "1.16.1+1"

[[LinearAlgebra]]
deps = ["Libdl"]
uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[[Logging]]
uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"

[[LoweredCodeUtils]]
deps = ["JuliaInterpreter"]
git-tree-sha1 = "491a883c4fef1103077a7f648961adbf9c8dd933"
uuid = "6f1432cf-f94c-5a45-995e-cdbf5db27b0b"
version = "2.1.2"

[[Markdown]]
deps = ["Base64"]
uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"

[[MbedTLS]]
deps = ["Dates", "MbedTLS_jll", "Random", "Sockets"]
git-tree-sha1 = "1c38e51c3d08ef2278062ebceade0e46cefc96fe"
uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
version = "1.0.3"

[[MbedTLS_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"

[[Millboard]]
git-tree-sha1 = "ea6a5b7e56e76d8051023faaa11d91d1d881dac3"
uuid = "39ec1447-df44-5f4c-beaa-866f30b4d3b2"
version = "0.2.5"

[[Mmap]]
uuid = "a63ad114-7e13-5084-954f-fe012c677804"

[[MozillaCACerts_jll]]
uuid = "14a3606d-f60d-562e-9121-12d972cd8159"

[[Nettle]]
deps = ["Libdl", "Nettle_jll"]
git-tree-sha1 = "a68340b9edfd98d0ed96aee8137cb716ea3b6dea"
uuid = "49dea1ee-f6fa-5aa6-9a11-8816cee7d4b9"
version = "0.5.1"

[[Nettle_jll]]
deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
git-tree-sha1 = "eca63e3847dad608cfa6a3329b95ef674c7160b4"
uuid = "4c82536e-c426-54e4-b420-14f461c4ed8b"
version = "3.7.2+0"

[[NetworkOptions]]
uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"

[[Nullables]]
git-tree-sha1 = "8f87854cc8f3685a60689d8edecaa29d2251979b"
uuid = "4d1e1d77-625e-5b40-9113-a560ec7a8ecd"
version = "1.0.0"

[[OrderedCollections]]
git-tree-sha1 = "85f8e6578bf1f9ee0d11e7bb1b1456435479d47c"
uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
version = "1.4.1"

[[Parsers]]
deps = ["Dates"]
git-tree-sha1 = "98f59ff3639b3d9485a03a72f3ab35bab9465720"
uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
version = "2.0.6"

[[Pkg]]
deps = ["Artifacts", "Dates", "Downloads", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "Serialization", "TOML", "Tar", "UUIDs", "p7zip_jll"]
uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"

[[Preferences]]
deps = ["TOML"]
git-tree-sha1 = "00cfd92944ca9c760982747e9a1d0d5d86ab1e5a"
uuid = "21216c6a-2e73-6563-6e65-726566657250"
version = "1.2.2"

[[Printf]]
deps = ["Unicode"]
uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"

[[REPL]]
deps = ["InteractiveUtils", "Markdown", "Sockets", "Unicode"]
uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"

[[Random]]
deps = ["Serialization"]
uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"

[[Reexport]]
git-tree-sha1 = "45e428421666073eab6f2da5c9d310d99bb12f9b"
uuid = "189a3867-3050-52da-a836-e630ba90ab69"
version = "1.2.2"

[[Requires]]
deps = ["UUIDs"]
git-tree-sha1 = "4036a3bd08ac7e968e27c203d45f5fff15020621"
uuid = "ae029012-a4dd-5104-9daa-d747884805df"
version = "1.1.3"

[[Revise]]
deps = ["CodeTracking", "Distributed", "FileWatching", "JuliaInterpreter", "LibGit2", "LoweredCodeUtils", "OrderedCollections", "Pkg", "REPL", "Requires", "UUIDs", "Unicode"]
git-tree-sha1 = "41deb3df28ecf75307b6e492a738821b031f8425"
uuid = "295af30f-e4ad-537b-8983-00126c2a3abe"
version = "3.1.20"

[[SHA]]
uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"

[[Serialization]]
uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"

[[SharedArrays]]
deps = ["Distributed", "Mmap", "Random", "Serialization"]
uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"

[[Sockets]]
uuid = "6462fe0b-24de-5631-8697-dd941f90decc"

[[SparseArrays]]
deps = ["LinearAlgebra", "Random"]
uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"

[[Statistics]]
deps = ["LinearAlgebra", "SparseArrays"]
uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"

[[StringEncodings]]
deps = ["Libiconv_jll"]
git-tree-sha1 = "50ccd5ddb00d19392577902f0079267a72c5ab04"
uuid = "69024149-9ee7-55f6-a4c4-859efe599b68"
version = "0.3.5"

[[StructTypes]]
deps = ["Dates", "UUIDs"]
git-tree-sha1 = "d24a825a95a6d98c385001212dc9020d609f2d4f"
uuid = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"
version = "1.8.1"

[[TOML]]
deps = ["Dates"]
uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"

[[Tar]]
deps = ["ArgTools", "SHA"]
uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"

[[Test]]
deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

[[TextWrap]]
git-tree-sha1 = "9250ef9b01b66667380cf3275b3f7488d0e25faf"
uuid = "b718987f-49a8-5099-9789-dcd902bef87d"
version = "1.0.1"

[[Tokenize]]
git-tree-sha1 = "0952c9cee34988092d73a5708780b3917166a0dd"
uuid = "0796e94c-ce3b-5d07-9a54-7f471281c624"
version = "0.5.21"

[[URIParser]]
deps = ["Unicode"]
git-tree-sha1 = "53a9f49546b8d2dd2e688d216421d050c9a31d0d"
uuid = "30578b45-9adc-5946-b283-645ec420af67"
version = "0.4.1"

[[URIs]]
git-tree-sha1 = "97bbe755a53fe859669cd907f2d96aee8d2c1355"
uuid = "5c2747f8-b7ea-4ff2-ba2e-563bfd36b1d4"
version = "1.3.0"

[[UUIDs]]
deps = ["Random", "SHA"]
uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"

[[Unicode]]
uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"

[[XML2_jll]]
deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Pkg", "Zlib_jll"]
git-tree-sha1 = "1acf5bdf07aa0907e0a37d3718bb88d4b687b74a"
uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
version = "2.9.12+0"

[[YAML]]
deps = ["Base64", "Dates", "Printf", "StringEncodings"]
git-tree-sha1 = "3c6e8b9f5cdaaa21340f841653942e1a6b6561e5"
uuid = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6"
version = "0.4.7"

[[Zlib_jll]]
deps = ["Libdl"]
uuid = "83775a58-1f1d-513f-b197-d71354ab007a"

[[nghttp2_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "8e850ede-7688-5339-a07c-302acd2aaf8d"

[[p7zip_jll]]
deps = ["Artifacts", "Libdl"]
uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
"""

# ╔═╡ Cell order:
# ╟─e28a7f22-319f-11ec-21e2-45db28a06407
# ╟─6e9844b4-bd77-4e0b-9908-a5afb2a1dd54
# ╠═e39a8330-8947-4227-83ea-bcd344601159
# ╠═e9a246a3-c6e7-4c0e-b220-61d21707cf57
# ╠═747ac8a7-ba1d-48be-98ec-e66756ca9b02
# ╠═7d3decbf-7276-448b-949a-ceea2ffd0b4b
# ╠═968c2c12-b564-432e-b879-70162b070462
# ╠═6e511e93-a104-481b-bb8d-35771768a1a2
# ╠═7a7bcced-c077-4531-9176-23ddc88c6886
# ╟─43a94071-cb2b-4815-a25b-614e85bd2a0f
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# Helper that returns a random N-dimensional ellipsoid
function random_ellipsoid(N::Integer)
    cov = rand(rng, N, N)
    while abs(det(cov)) < 1e-10
        cov = rand(rng, N, N)
    end
    return Ellipsoid(zeros(N), pinv(cov * cov'))
end

const BOUNDST = [
    Bounds.Ellipsoid,
    Bounds.MultiEllipsoid
]

@testset "interface - $B, $T, D=$D" for B in BOUNDST, T in [Float32, Float64], D in 1:20
    # creation, inspection
    bound = B(T, D)
    @test eltype(bound) == T
    @test ndims(bound) == D

    # sampling
    sample = rand(rng, bound)

    @test eltype(sample) == T
    @test size(sample) == (D,)
    @test sample ∈ bound

    nsamples = 1000
    samples = rand(rng, bound, nsamples)

    @test eltype(samples) == T
    @test size(samples) == (D, nsamples)
    @test all(samples[:, i] ∈ bound for i in axes(samples, 2))

    # fitting
    bound = Bounds.fit(B, samples)
    @test eltype(bound) == T
    @test all(samples[:, i] ∈ bound for i in axes(samples, 2))

    # robust fitting
    pv = Bounds.volume(bound) / size(samples, 2)
    bound2 = Bounds.fit(B, samples; pointvol = pv)
    @test Bounds.volume(bound2) ≈ Bounds.volume(bound) rtol = 1e-3

    # volume and scaling
    volfrac = 0.5
    bound_scaled = Bounds.scale!(deepcopy(bound), volfrac)
    @test Bounds.volume(bound) ≈ Bounds.volume(bound_scaled) / volfrac rtol = 1e-3

    # expected number of points that will fall within inner bound
    npoints = 5000
    expect = volfrac * npoints
    σ = sqrt((1 - volfrac) * expect)
    ninner = count(rand(rng, bound) ∈ bound_scaled for _ in 1:npoints)
    @test ninner ≈ expect atol = 3σ

    # printing
    @test sprint(show, bound) == "$(string(nameof(B))){$T}(ndims=$D)"

    # rand_live
    x = rand(rng, bound, 10)
    point, _bound = Bounds.rand_live(rng, bound, x)
    count(point ∈ x[:, i] for i in axes(x, 2)) == 1
    Btarget = B ∈ [Bounds.MultiEllipsoid] ? Bounds.Ellipsoid : B
    @test _bound isa Btarget
    @test point ∈ _bound && point ∈ bound
end

@testset "interface - NoBounds, $T, D=$D" for T in [Float16, Float32, Float64], D in 1:20
    # creation, inspection
    bound = Bounds.NoBounds(T, D)
    @test bound == Bounds.NoBounds{T}(D)
    @test eltype(bound) == T
    @test ndims(bound) == D

    # sampling
    sample = rand(rng, bound)

    @test eltype(sample) == T
    @test length(sample) == D
    @test sample ∈ bound

    samples = rand(rng, bound, 3)

    @test eltype(samples) == T
    @test size(samples) == (D, 3)
    @test all(samples[:, i] ∈ bound for i in axes(samples, 2))

    # fitting
    samples = randn(rng, T, D, 100)
    @test Bounds.fit(Bounds.NoBounds, samples) == bound
    @test eltype(bound) == T # matches eltype

    # robust fitting
    pv = 1 / size(samples, 2)
    bound_fit = Bounds.fit(Bounds.NoBounds, samples; pointvol = pv)
    @test bound_fit == bound
    @test Bounds.volume(bound_fit) == 1

    # volume and scaling
    bound_scaled = Bounds.scale!(deepcopy(bound), 0.5)
    @test bound_scaled == bound
    @test Bounds.volume(bound_scaled) == 1

    @test Bounds.axes(bound) == Diagonal(ones(D))
end

include("helpers.jl")
include("ellipsoids.jl")


text:
# example 2: add some flags and the help lines for options

using ArgParse

function main(args)

    s = ArgParseSettings("Example 2 for argparse.jl: " *  # description
                         "flags, options help, " *
                         "required arguments.")

    @add_arg_table s begin
        "--opt1"
            help = "an option"     # used by the help screen
        "--opt2", "-o"
            action = :store_true   # this makes it a flag
            help = "a flag"
        "arg1"
            help = "an argument"
            required = true        # makes the argument mandatory
    end

    parsed_args = parse_args(args, s)
    println("Parsed args:")
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end
end

main(ARGS)


text:
function calculateLossReallocation(edges, capacity, demand, flows, T, Tf, k, splittingratios, scenarios, probabilities; progress=false)
    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)
    nscenarios = length(scenarios)

    #CREATE TUNNEL SCENARIO MATRIX
    X  = ones(nscenarios,ntunnels)
    for s in 1:nscenarios
        for t in 1:ntunnels
            if size(T[t],1) == 0
                X[s,t] = 0
            else
                for e in 1:nedges
                    if scenarios[s][e] == 0
                        back_edge = findfirst(x -> x == (edges[e][2],edges[e][1]), edges)
                        if in(e, T[t]) || in(back_edge, T[t])
                            X[s,t] = 0
                        end
                    end
                end
            end
        end
    end

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    as = zeros(nflows,k)
    routed = zeros(nscenarios, nflows, k)
    u = zeros(nscenarios, nflows)
    t = zeros(nscenarios)

    # SCENARIO LOSS PER FLOW
    for s in 1:nscenarios
        for f in 1:nflows
            totalup = 0
            for t in 1:size(Tf[f],1)
                totalup += splittingratios[f,t] * X[s,Tf[f][t]]
            end
            if totalup == 0
                splittingratios[f,:] = splittingratios[f,:] .+ .2
                for t in 1:size(Tf[f],1)
                    totalup += splittingratios[f,t] * X[s,Tf[f][t]]
                end
            end

            for t in 1:size(Tf[f],1)
                # as[f,t] = max(as[f,t], splittingratios[f,t] / totalup * demand[f] * X[s,Tf[f][t]])
                if totalup != 0
                    routed[s,f,t] = splittingratios[f,t] / totalup * demand[f] * X[s,Tf[f][t]]
                end
            end
        end

        t[s] = sum(routed[s,:,:])

        congestion_loss = 0
        for e in 1:nedges
            edge_utilization = 0
            for f in 1:nflows
                edge_utilization += sum(routed[s,f,t] * L[Tf[f][t],e] *  X[s,Tf[f][t]] for t in 1:size(Tf[f],1))
                for t in enumerate(Tf[f])
                    # edge_utilization += sum(as[f,t[1]] * L[t[2],e] * X[s,t[2]])
                end
            end
            # println("Edge: ", edges[e])
            # println(max(0, edge_utilization - capacity[e]))
            congestion_loss += max(0, round((edge_utilization - capacity[e])*1000)/1000)
        end
        t[s] -= congestion_loss
    end

    umax = map(x -> round((1 - x/sum(demand))*100000)/100000, t)
    return umax
end

function PDF(losses, probabilities, sla)
    usorted = []
    psorted = []
    umodified = losses
    pmodified = probabilities
    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    for s in 1:length(usorted)
        loss = usorted[s]
        if loss > sla
            break
        end
        total += psorted[s]

    end
    return total
end

function PDF(losses, sla)
    c = 0
    for s in 1:length(losses)
        if (1-losses[s]) >= sla
            c += 1
        end
    end
    return c / length(losses)
end

function VarUniform(losses, beta)
    usorted = sort(losses)
    probabilities = zeros(length(losses)) .+ 1/length(losses)

    total = 0
    loss = 0
    varindex = 0
    for s in 1:size(usorted, 1)
        total += probabilities[s]
        loss = usorted[s]
        if total >= beta
            break
        end
    end
    return loss
end

function VAR(losses, probabilities, beta)
    usorted = []
    psorted = []
    umodified = losses
    pmodified = probabilities
    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    varindex = 0
    for s in 1:size(usorted, 1)
        total += psorted[s]
        loss = usorted[s]
        if total >= beta
            break
        end
    end
    return loss
end

function CVAR(losses, probabilities, beta)
    usorted = []
    psorted = []
    umodified = losses
    pmodified = probabilities
    while length(umodified) > 0
        s = argmin(umodified)
        push!(usorted, umodified[s])
        push!(psorted, pmodified[s])
        umodified = umodified[1:end .!= s]
        pmodified = pmodified[1:end .!= s]
    end

    total = 0
    loss = 0
    prob_total = 0
    for s in 1:length(usorted)
        total += psorted[s]
        if total >= beta
            prob_total += psorted[s]
            loss += usorted[s]*psorted[s]
        end
    end
    return loss / prob_total
end

function simulateUtilizationNoFailures(edges, capacity, demand, flows, T, Tf, k, a, bandwidth_allowed)
    nedges = length(edges)
    nflows = length(flows)
    ntunnels = length(T)

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    routed = zeros(nflows, k)
    for f in 1:nflows
        totalup = 0
        for t in 1:size(Tf[f],1)
            totalup += a[f,t]
        end
        if totalup == 0
            a[f,:] = a[f,:] .+ .2
            for t in 1:size(Tf[f],1)
                totalup += a[f,t]
            end
        end

        for t in 1:size(Tf[f],1)
            routed[f,t] = a[f,t] / totalup * bandwidth_allowed[f] * demand[f]
        end
    end

    edge_utilization = zeros(nedges)
    edge_utilization_percentage = zeros(nedges)

    for e in 1:nedges
        for f in 1:nflows
            edge_utilization[e] += sum(routed[f,t] * L[Tf[f][t],e] for t in 1:size(Tf[f],1))
        end
        edge_utilization_percentage = edge_utilization[e]/capacity[e]
    end

    return edge_utilization, routed
end


text:
# A complicated way to do argmin x^1/2 - x^1/3
module TestExampleDiagonal

using DiffEqBase
using Setfield
using SteadyStateFit
using Test

struct NullSetter <: Lens end
Setfield.set(obj, ::NullSetter, ::Any) = obj

ode = ODEProblem(
    ODEFunction(
        (u, p, t) -> [u[1]^2, u[2]^3] .- exp.(p);
        jac = (u, p, t) -> [
            2u[1]   0
            0       3u[2]^2
        ],
    ),
    [1.0, 1.0],  # u0
    nothing,
    [1.0],  # p
)

sso = SteadyStateObjective(
    (((x1, x2),), (c,), _) -> x1 - c * x2,  # loss
    ode,
    [1],                           # conditions
    NullSetter(),                  # conditionsetter
    (@lens _),                     # parameterlens
)

result = solve(sso)

@test exp(Optim.minimizer(result)[1]) ≈ 64/729  rtol=0.01
@test occursin("SteadyStateFitResult", sprint(show, "text/plain", result))
@test occursin("SteadyStateFitResult: success",
               sprint(show, result; context=:limit => true))
@test "SteadyStateObjective: 1 conditions 1 parameters" ==
    sprint(show, sso; context=:limit => true)

result2 = solve(sso, Optim.BFGS(), Optim.Options(store_trace=true))
@test exp(Optim.minimizer(result2)[1]) ≈ 64/729  rtol=0.01

end  # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct DeploymentOperationProperties <: SwaggerModel
    provisioningOperation::Any # spec type: Union{ Nothing, String } # spec name: provisioningOperation
    provisioningState::Any # spec type: Union{ Nothing, String } # spec name: provisioningState
    timestamp::Any # spec type: Union{ Nothing, DateTime } # spec name: timestamp
    duration::Any # spec type: Union{ Nothing, String } # spec name: duration
    serviceRequestId::Any # spec type: Union{ Nothing, String } # spec name: serviceRequestId
    statusCode::Any # spec type: Union{ Nothing, String } # spec name: statusCode
    statusMessage::Any # spec type: Union{ Nothing, StatusMessage } # spec name: statusMessage
    targetResource::Any # spec type: Union{ Nothing, TargetResource } # spec name: targetResource
    request::Any # spec type: Union{ Nothing, HttpMessage } # spec name: request
    response::Any # spec type: Union{ Nothing, HttpMessage } # spec name: response

    function DeploymentOperationProperties(;provisioningOperation=nothing, provisioningState=nothing, timestamp=nothing, duration=nothing, serviceRequestId=nothing, statusCode=nothing, statusMessage=nothing, targetResource=nothing, request=nothing, response=nothing)
        o = new()
        validate_property(DeploymentOperationProperties, Symbol("provisioningOperation"), provisioningOperation)
        setfield!(o, Symbol("provisioningOperation"), provisioningOperation)
        validate_property(DeploymentOperationProperties, Symbol("provisioningState"), provisioningState)
        setfield!(o, Symbol("provisioningState"), provisioningState)
        validate_property(DeploymentOperationProperties, Symbol("timestamp"), timestamp)
        setfield!(o, Symbol("timestamp"), timestamp)
        validate_property(DeploymentOperationProperties, Symbol("duration"), duration)
        setfield!(o, Symbol("duration"), duration)
        validate_property(DeploymentOperationProperties, Symbol("serviceRequestId"), serviceRequestId)
        setfield!(o, Symbol("serviceRequestId"), serviceRequestId)
        validate_property(DeploymentOperationProperties, Symbol("statusCode"), statusCode)
        setfield!(o, Symbol("statusCode"), statusCode)
        validate_property(DeploymentOperationProperties, Symbol("statusMessage"), statusMessage)
        setfield!(o, Symbol("statusMessage"), statusMessage)
        validate_property(DeploymentOperationProperties, Symbol("targetResource"), targetResource)
        setfield!(o, Symbol("targetResource"), targetResource)
        validate_property(DeploymentOperationProperties, Symbol("request"), request)
        setfield!(o, Symbol("request"), request)
        validate_property(DeploymentOperationProperties, Symbol("response"), response)
        setfield!(o, Symbol("response"), response)
        o
    end
end # type DeploymentOperationProperties

const _property_map_DeploymentOperationProperties = Dict{Symbol,Symbol}(Symbol("provisioningOperation")=>Symbol("provisioningOperation"), Symbol("provisioningState")=>Symbol("provisioningState"), Symbol("timestamp")=>Symbol("timestamp"), Symbol("duration")=>Symbol("duration"), Symbol("serviceRequestId")=>Symbol("serviceRequestId"), Symbol("statusCode")=>Symbol("statusCode"), Symbol("statusMessage")=>Symbol("statusMessage"), Symbol("targetResource")=>Symbol("targetResource"), Symbol("request")=>Symbol("request"), Symbol("response")=>Symbol("response"))
const _property_types_DeploymentOperationProperties = Dict{Symbol,String}(Symbol("provisioningOperation")=>"String", Symbol("provisioningState")=>"String", Symbol("timestamp")=>"DateTime", Symbol("duration")=>"String", Symbol("serviceRequestId")=>"String", Symbol("statusCode")=>"String", Symbol("statusMessage")=>"StatusMessage", Symbol("targetResource")=>"TargetResource", Symbol("request")=>"HttpMessage", Symbol("response")=>"HttpMessage")
Base.propertynames(::Type{ DeploymentOperationProperties }) = collect(keys(_property_map_DeploymentOperationProperties))
Swagger.property_type(::Type{ DeploymentOperationProperties }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_DeploymentOperationProperties[name]))}
Swagger.field_name(::Type{ DeploymentOperationProperties }, property_name::Symbol) =  _property_map_DeploymentOperationProperties[property_name]

const _allowed_DeploymentOperationProperties_provisioningOperation = ["NotSpecified", "Create", "Delete", "Waiting", "AzureAsyncOperationWaiting", "ResourceCacheWaiting", "Action", "Read", "EvaluateDeploymentOutput", "DeploymentCleanup"]

function check_required(o::DeploymentOperationProperties)
    true
end

function validate_property(::Type{ DeploymentOperationProperties }, name::Symbol, val)
    if name === Symbol("provisioningOperation")
        Swagger.validate_param(name, "DeploymentOperationProperties", :enum, val, _allowed_DeploymentOperationProperties_provisioningOperation)
    end
end