text:
"""
    lagrange_polynomials(nodes::Array{Float64, 1}, evalPt::Float64)

Evaluate the lagrange interpolating polynomials defied by `nodes` at the point `evalPt`.
"""
function lagrange_polynomials(nodes::Array{Float64,1}, evalPt::Float64)
    order = length(nodes)
    poly_coefs = zeros(order)
    if order > 1
        for n in 1:order
            notn = [1:n-1;n+1:order]
            poly_coefs[n] = prod(evalPt.-nodes[notn])./prod(nodes[n].-nodes[notn])
        end
    else
        poly_coefs = [1.0]
    end
    return poly_coefs
end

"""
    gauss_lobatto_points(a::Float64, b::Float64, n_evals::Int)

Return array containing `n_evals` gauss lobatto points in the interval `[a,b]`.
"""
function gauss_lobatto_points(a::Float64,b::Float64,n_evals::Int)
    (b<a)&&throw(DomainError("must have a<b"))
    (n_evals<1)&&throw(DomainError("n_evals must be > 0"))
    if n_evals > 1
        # the GL nodes
        # nodes in [-1,1]
        nodes = Jacobi.zglj(n_evals, 0, 0)
        # shift nodes to [a,b]
        nodes *= 0.5*(b-a)
        nodes .+= 0.5*(a+b)
    else
        nodes = [0.5*(a+b)]
    end
    return nodes
end

"""
    gauss_lobatto_weights(a::Float64, b::Float64, n_evals::Int)

Return array containing `n_evals` weights of the polynomials associated with the
gauss lobatto points in the interval `[a,b]`. i.e. the itegrals over the interval `[a,b]`
of the lagrange polynomials defined by the nodes given by `nodes = gauss_lobatto_points(a,b,n_evals)`
"""
function gauss_lobatto_weights(a::Float64,b::Float64,n_evals::Int)
    (b<a)&&throw(DomainError("must have a<b"))
    (n_evals<1)&&throw(DomainError("n_evals must be > 0"))
    if n_evals > 1
        nodes = Jacobi.zglj(n_evals, 0, 0)
        weights = Jacobi.wglj(nodes,0,0)*(b-a)/2.0
    else
        weights = [(b-a)]
    end
    return weights
end

"""
    lagrange_interpolation(fun::Function, a::Float64, b::Float64, n_evals::Int)

Return a polynomial approximation to `fun` of order `n_evals-1` on the interval `[a,b]` using
the `gauss_lobatto_points` as the nodes of the lagrange polynomials.
"""
function lagrange_interpolation(fun::Function,a::Float64,b::Float64,n_evals::Int)
    nodes = gauss_lobatto_points(a,b,n_evals)

    fun_vals = fun.(nodes)
    interpolant(x) = LinearAlgebra.dot(fun_vals,lagrange_polynomials(nodes,x))
    return interpolant
end

"""
    gauss_lobatto_quadrature(fun::Function, a::Float64, b::Float64, n_evals::Int)

Compute a quadrature approximation of `fun` on the interval `[a,b]` with `n_evals` function
evaluations using the `gauss_lobatto_points`
"""
function gauss_lobatto_quadrature(fun::Function,a::Float64,b::Float64,n_evals::Int)
    nodes = gauss_lobatto_points(a,b,n_evals)
    weights = gauss_lobatto_weights(a,b,n_evals)
    fun_vals = fun.(nodes)
    quad = LinearAlgebra.dot(fun_vals,weights)

    return quad
end
"""
    gauss_lobatto_quadrature(fun::Function, nodes::Array{Float64,1})

Compute a quadrature approximation of `fun` on the interval `nodes[1]` to `nodes[end]`
evaluating the function `fun` at `nodes`
"""
function gauss_lobatto_quadrature(fun::Function,nodes::Array{Float64,1})
    weights = gauss_lobatto_weights(a,b,length(nodes))
    fun_vals = fun.(nodes)
    quad = LinearAlgebra.dot(fun_vals,weights)

    return quad
end

"""
Construct a generalised vandermonde matrix.

    vandermonde( nBases::Int)

Note: requires Jacobi package Pkg.add("Jacobi")

# Arguments
- `nBases::Int`: the degree of the basis

# Output
- a tuple with keys
    - `:V::Array{Float64,2}`: where `:V[:,i]` contains the values of the `i`th
        legendre polynomial evaluate at the GLL nodes.
    - `:inv`: the inverse of :V
    - `:D::Array{Float64,2}`: where `V.D[:,i]` contains the values of the derivative
        of the `i`th legendre polynomial evaluate at the GLL nodes.
"""
function vandermonde(nBases::Int)
    if nBases > 1
        z = Jacobi.zglj(nBases, 0, 0) # the LGL nodes
    elseif nBases == 1
        z = 0.0
    end
    V = zeros(Float64, nBases, nBases)
    DV = zeros(Float64, nBases, nBases)
    if nBases > 1
        for j = 1:nBases
            # compute the polynomials at gauss-labotto quadrature points
            V[:, j] = Jacobi.legendre.(z, j - 1) .* sqrt((2 * (j - 1) + 1) / 2)
            DV[:, j] = Jacobi.dlegendre.(z, j - 1) .* sqrt((2 * (j - 1) + 1) / 2)
        end
        # Compute the Gauss-Lobatto weights for numerical quadrature
        w =
            2.0 ./ (
                nBases *
                (nBases - 1) *
                Jacobi.legendre.(Jacobi.zglj(nBases, 0, 0), nBases - 1) .^ 2
            )
    elseif nBases == 1
        V .= [1.0/sqrt(2)]
        DV .= [0.0]
        w = [2.0]
    end
    return (V = V, inv = inv(V), D = DV, w = w)
end

function legendre_to_lagrange(coeffs::Vector{Float64})
    order = length(coeffs)
    V = vandermonde(order)
    return V.V*coeffs
end

function lagrange_to_legendre(coeffs::Vector{Float64})
    order = length(coeffs)
    V = vandermonde(order)
    return V.inv*coeffs
end

function legendre_to_lagrange(coeffs::Array{Float64,3})
    order = size(coeffs,1)
    V = vandermonde(order)
    for  i in 1:size(coeffs,3)
        coeffs[:,:,i] = V.V*coeffs[:,:,i]
    end
    return coeffs
end

function lagrange_to_legendre(coeffs::Array{Float64,3})
    order = size(coeffs,1)
    V = vandermonde(order)
    for  i in 1:size(coeffs,3)
        coeffs[:,:,i] = V.inv*coeffs[:,:,i]
    end
    return coeffs
end


text:
# shorten float strings
xprint(x) = sprint(show, x, context = :compact => true)
# return plain name of object
wrapper_name(x) = typeof(x).name.wrapper

# general single laser verbose string
function laser_description(s)
    unit_system = fundamental_constants(s, :unit_system)
    rows = String[]
    push!(rows, "$(wrapper_name(s)) with $unit_system units")
    a_0 = xprint(s.a₀)
    push!(rows, "a₀ = "*a_0)
    omega  = xprint(immutable_cache(s, :ω))
    push!(rows, "ω  = "*omega)
    phi_0 = xprint(s.ϕ₀)
    push!(rows, "ϕ₀ = "*phi_0)
    push!(rows, "and temporal $(s.profile)")
    return rows
end

# Abstract Laser
function Base.show(io::IO, ::MIME"text/plain", s::AbstractLaser)
    laser_string = join(laser_description(s),"\n")
    return print(io, laser_string)
end

# Gaussian Laser
function Base.show(io::IO, ::MIME"text/plain", s::GaussLaser)
    rows = laser_description(s)
    lastrow = rows[end]
    w_0 = xprint(s.w₀)
    rows[end] = "w₀ = "*w_0
    push!(rows, "with $(s.polarization)")
    push!(rows, lastrow)
    laser_string = join(rows, "\n")
    return print(io, laser_string)
end

# Laguerre-Gauss Laser
function Base.show(io::IO, ::MIME"text/plain", s::LaguerreGaussLaser)
    rows = laser_description(s)
    lastrow = rows[end]
    w_0 = xprint(s.w₀)
    rows[end] = "w₀ = "*w_0
    push!(rows, "p = $(s.p) and m = $(s.m)")
    push!(rows, "with $(s.polarization)")
    push!(rows, lastrow)
    laser_string = join(rows, "\n")
    return print(io, laser_string)
end

# general single laser brief string
function laser_brief(s)
    a_0 = xprint(s.a₀)
    omega = xprint(immutable_cache(s, :ω))
    phi_0 = xprint(s.ϕ₀)
    temp = wrapper_name(s.profile)
    rows = String[]
    push!(rows, "a₀ = $a_0")
    push!(rows, "ω = $omega")
    push!(rows, "ϕ₀ = "*phi_0)
    push!(rows, "temporal $temp")
    return rows
end

# multiobject/brief print
# AbstractLaser
function Base.show(io::IO, s::AbstractLaser)
    name = wrapper_name(s)
    rows = laser_brief(s)
    laser_string = join(rows, ", ")
    return print(io, "$name: "*laser_string)
end

# GaussLaser
function Base.show(io::IO, s::GaussLaser)
    rows = laser_brief(s)
    lastrow = rows[end]
    rows[end] = "w₀ = $(xprint(s.w₀))"
    push!(rows, lastrow)
    laser_string = join(rows, ", ")
    return print(io, "Gaussian laser: "*laser_string)
end

# LaguerreGaussLaser
function Base.show(io::IO, s::LaguerreGaussLaser)
    rows = laser_brief(s)
    lastrow = rows[end]
    rows[end] = "w₀ = $(xprint(s.w₀))"
    push!(rows, "p = $(s.p)")
    push!(rows, "m = $(s.m)")
    push!(rows, lastrow)
    laser_string = join(rows, ", ")
    return print(io, "Laguerre-Gauss laser: "*laser_string)
end

# polarization show
function Base.show(io::IO, p::LaserPolarization)
    @unpack ξx, ξy = p
    tangent = ξy/ξx
    if abs(ξx) ≈ abs(ξy)
        angle(tangent) ≈  π/2 && return print(io, "right-handed polarization")
        angle(tangent) ≈ -π/2 && return print(io, "left-handed polarization")
    else
        angle(tangent) ≈  0.  && return print(io, "linear polarization")
        angle(tangent) ≈  π   && return print(io, "linear polarization")
        return print(io, "eliptical polarization")
    end
end

# temporal profile show
# AbstractTemporalProfile
function Base.show(io::IO, p::AbstractTemporalProfile)
    print(io, wrapper_name(p))
end

# GaussProfile
function Base.show(io::IO, p::GaussProfile)
    @unpack inv_τ, t₀, z₀ = p
    rows = String[]
    push!(rows, "Gaussian profile")
    z0, t0 = xprint(z₀), xprint(t₀)
    push!(rows, "centered in z₀ = $z0 and t₀ = $t0")
    tau = xprint(inv(inv_τ))
    push!(rows, "with duration of pulse (FWHM) τ = $tau")
    laser_string = join(rows, "\n")
    return print(io, laser_string)
end

# Cos²Profile
function Base.show(io::IO, p::Cos²Profile)
    @unpack τ, t₀, z₀ = p
    rows = String[]
    push!(rows, "Cos² profile")
    z0, t0 = xprint(z₀), xprint(t₀)
    push!(rows, "centered in z₀ = $z0 and t₀ = $t0")
    tau = xprint(τ)
    push!(rows, "with duration of pulse τ = $tau")
    laser_string = join(rows, "\n")
    return print(io, laser_string)
end

# QuasiRectangularProfile
function Base.show(io::IO, p::QuasiRectangularProfile)
    @unpack Δt, inv_τ, t₀, z₀ = p
    rows = String[]
    push!(rows, "Quasi rectangular profile")
    z0, t0 = xprint(z₀), xprint(t₀)
    push!(rows, "centered in z₀ = $z0 and t₀ = $t0")
    tau = xprint(inv(inv_τ))
    push!(rows, "with exponential decay width τ = $tau")
    dt = xprint(Δt)
    push!(rows, "length of the plateau Δt = $dt")
    laser_string = join(rows, "\n")
    return print(io, laser_string)
end

ext: jl


text:
include("cg_sem.jl")
using ForwardDiff: derivative
using LinearAlgebra: eigen

function mms_test(N, K, z0, z1, cfl)
    elem = element_operators(N)
    mesh = create_mesh(K, elem, z0, z1)

    problem = Isothermal{Float64}()

    # mms solution
    f(z, t) = cos(t) * cos(2π * (z - z0) / (z1 - z0))
    ρ′(z, t) = f(1z, 4t)
    w′(z, t) = f(2z, 2t)
    E′(z, t) = f(4z, 1t)

    ∂tρ′(z, t) = derivative(t -> ρ′(z, t), t)
    ∂tw′(z, t) = derivative(t -> w′(z, t), t)
    ∂tE′(z, t) = derivative(t -> E′(z, t), t)

    ∂zρ̄w(z, t) = derivative(z -> ρ̄_(problem, z) * w′(z, t), z)
    ∂zδp(z, t) = derivative(z -> δp_(problem, z, E′(z, t), ρ′(z, t)), z)
    ∂zĒ_p̄w′(z, t) =
        derivative(z -> (Ē_(problem, z) + p̄_(problem, z)) * w′(z, t), z)

    # Forcing for each term
    fρ(z, t) = ∂tρ′(z, t) + ∂zρ̄w(z, t)

    grav_source(z, t) = problem.g * ρ′(z, t)
    fw(z, t) = ∂tw′(z, t) + (grav_source(z, t) + ∂zδp(z, t)) / ρ̄_(problem, z)

    fE(z, t) = ∂tE′(z, t) + ∂zĒ_p̄w′(z, t)

    forcing = (ρ = fρ, w = fw, E = fE)

    # Matrices for upwinding the boundary conditions
    A0 = A_(problem, z0)
    (λ, V) = eigen(A0)
    A0⁺ = V * Diagonal(max.(0, λ)) / V
    A0⁻ = V * Diagonal(min.(0, λ)) / V

    A1 = -A_(problem, z1)
    (λ, V) = eigen(A1)
    A1⁺ = V * Diagonal(max.(0, λ)) / V
    A1⁻ = V * Diagonal(min.(0, λ)) / V

    bc = (
        (A = A0, A⁺ = A0⁺, A⁻ = A0⁻, δw = t -> w′(z0, t)),
        (A = A1, A⁺ = A1⁺, A⁻ = A1⁻, δw = t -> w′(z1, t)),
    )

    # Calculate the wave speed
    c̄ = maximum(c̄_(problem, mesh.zdg))

    # estimate time step
    dt = cfl * mesh.Δz / (N^2 * c̄)

    # Make sure we end at the right time
    tspan = (0, 11π)
    steps = ceil(Int, (tspan[2] - tspan[1]) / dt)
    dt = (tspan[2] - tspan[1]) / steps

    # Initial conditions
    q = (
        δρ = ρ′.(mesh.zdg, tspan[1]),
        δw = w′.(mesh.zdg, tspan[1]),
        δE = E′.(mesh.zdg, tspan[1]),
    )

    # advance in time
    rhs!(∂q, q, t) = tendency!(∂q, q, t, problem, mesh, elem, bc, forcing)
    timestep!(q, rhs!, dt, tspan)

    # Exact solution at the final time
    qe = (
        δρ = ρ′.(mesh.zdg, tspan[2]),
        δw = w′.(mesh.zdg, tspan[2]),
        δE = E′.(mesh.zdg, tspan[2]),
    )

    # Calculate the L2 error
    ϵ = zeros(3)
    for i in 1:3
        ϵ[i] = sqrt(sum(mesh.J .* elem.ω .* (q[i] - qe[i]) .^ 2))
    end

    return ϵ
end

function mms_run(N, Ks; cfl = 1)
    ϵ = zeros(3, length(Ks))
    for (lvl, K) in enumerate(Ks)
        @show lvl, K
        ϵ[:, lvl] .= mms_test(N, K, 0, 100π, cfl)
        @show ϵ[:, lvl]
        if lvl > 1
            rate =
                (log.(ϵ[:, lvl]) - log.(ϵ[:, lvl - 1])) ./
                (log.(Ks[lvl - 1]) - log.(Ks[lvl]))
            @show rate
        end
        println()
    end
end


text:
using DelimitedFiles
using Random

export dump_params, load_params

# chimera
function dump_params(folder, sg::Spinglass{LT}, jt::JT, ht::HT, seed::Int) where {LT,JT<:Jtype, HT<:Jtype}
    jfilename = joinpath(folder, _gen_jfname(sg, jt, seed))
    hfilename = joinpath(folder, _gen_hfname(sg, ht, seed))
    open(jfilename, "w") do f
        write(f, "# $JT\n")
        for ((i,j), J) in zip(sgbonds(sg.lattice), sg.Js)
            write(f, "$i $j $J\n")
        end
    end
    open(hfilename, "w") do f
        write(f, "# $HT\n")
        for (i, h) in zip(sgvertices(sg.lattice), sg.hs)
            write(f, "$i $h\n")
        end
    end
    return sg
end

function _gen_jfname(sg::Spinglass{LT}, ::JT, seed::Int) where {LT,JT, HT}
    sz = size(sg.lattice)
    return "$(LT)_J$(JT)_Lx$(sz[end-1])_Ly$(sz[end])_seed$(seed).dat"
end

function _gen_hfname(sg::Spinglass{LT}, ::HT, seed::Int) where {LT,JT, HT}
    sz = size(sg.lattice)
    return "$(LT)_h$(HT)_Lx$(sz[end-1])_Ly$(sz[end])_seed$(seed).dat"
end

function load_params(folder, sg::Spinglass{LT,T}, jt::JT, ht::HT, seed::Int) where {LT,T,JT<:Jtype, HT<:Jtype}
    jfilename = joinpath(folder, _gen_jfname(sg, jt, seed))
    hfilename = joinpath(folder, _gen_hfname(sg, ht, seed))
    open(jfilename, "r") do f
        readline(f)
        for (k, (i,j)) in enumerate(sgbonds(sg.lattice))
            line = strip(readline(f))
            sg.Js[k] = parse(T, split(line, ' ')[3])
        end
    end
    open(hfilename, "r") do f
        readline(f)
        for i in 1:length(sgvertices(sg.lattice))
            line = strip(readline(f))
            a, b = split(line, ' ')
            sg.hs[i] = parse(T, b)
        end
    end
    return sg
end


text:
##############
using Revise
using FDMIC_Geodynamics
using LoopVectorization, Printf, Base.Threads, Plots, Revise, LinearAlgebra, Statistics, SparseArrays
##############
function SetMarkers!( p::Markers, params::ModelParameters, dom::ModelDomain )
    R = params.user[1]
    # Use this function to set up the model geometry
     for k=1:p.nmark #@tturbo <---------rand does not work in @tturbo!!!!
        in         = ( ( (p.x[k])^2 + p.y[k]^2) < R^2 )
        p.phase[k] = (in==0)* 1.0 + (in==1)*2.0
    end
end
##############
function Rheology!( f, Eps, Tau, params, materials, BC, ncx, ncy )
    # Compute effective viscosity for each phase
    f.etac .= 0.0
    for m=1:materials.nphase
        eta      = materials.eta0[m] * Eps.II.^(1.0/materials.n[m] - 1.0)
        f.etac .+= f.phase_perc[m,:,:] .* eta
    end
    @time CentroidsToVertices!( f.etav, f.etac, ncx, ncy, BC )
    println("min Eii: ", minimum(Eps.II), " --- max Eii: ", maximum(Eps.II))
    println("min eta: ", minimum(f.etac), " --- max eta: ", maximum(f.etac))
end
export Rheology!
##############
@views function main( )
    # Main routine
    println("#####################")
    println("###### M2Di.jl ######")
    println("#####################")
    # Name of the experiement
    experiment       = "PeriodicSimpleShear"
    # Domain
    domain           = ModelDomain()
    domain.xmin      = -1.0
    domain.xmax      =  1.0
    domain.ymin      = -1.0
    domain.ymax      =  1.0
    # parameters
    params           = ModelParameters()
    params.Ebg       = -1.0      # reference strain rate (>0 horizontal compression)
    # Scales
    scale            = CharScales()
    scale.L          = 1.0
    scale.T          = 1.0
    scale.t          = 1.0/abs(params.Ebg)
    # Spatial discretisation
    ncx              = 100; domain.ncx = ncx
    ncy              = 100; domain.ncy = ncy
    # Time discretisation
    params.advection = 1
    params.nt        = 200
    params.dt        = 1e4
    params.Courant   = 0.40  # Courant number
    params.dt_var    = 1     # variable dt
    params.t         = 0.0   # time
    # Boundary conditions
    params.BC_type   = 3     # 1: Pure shear / 2: Simple shear / 3: Simple shear periodic
    params.ALE       = 0     # Deform box
    params.comp      = 0     # Compressible
    # Solver
    params.solver    = 1     # 0: coupled --- 1: decoupled
    params.gamma     = 1e4   # Penalty factor
    params.Dir_scale = 1.0   # Dirichlet scaling factor
    # Non-linear iterations
    params.niter_nl  = 5     # max. number of non-linear iterations
    params.tol_nl    = 1e-7  # non-linear tolerance
    params.JFNK      = 0     # Jacobian-Free Newton_Krylov
    # Visualisation
    params.show_figs = 1     # activates visualisation...
    params.nout      = 10    # ... every nout
    # Initialize coordinates
    domain.L0 = domain.xmax-domain.xmin
    GenerateMesh!( domain )
    dx, dy = domain.dx, domain.dy
    # Material setup
    params.user[1]        = 1.0/3.0  # layer half thickness
    # Material parameters
    materials        = MaterialParameters()
    materials.nphase = 2
    materials.Tref   = zeros(Float64, materials.nphase)
    materials.eta0   = zeros(Float64, materials.nphase)
    materials.n      = zeros(Float64, materials.nphase)
    # Material 1
    materials.Tref[1]   = 2000.0  # reference flow stress
    materials.n[1]      = 1.0
    materials.eta0[1]   = materials.Tref[1] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[1])  # matrix viscosity
    # Material 2
    materials.Tref[2]   = 2 # reference flow stress
    materials.n[2]      = 1.0
    materials.eta0[2]   = materials.Tref[2] * (1.0/2.0) * abs(params.Ebg)^(-1.0/materials.n[2]) # inclusion viscosity
    # Output arrays for monitoring
    Pana = zeros(Float64, params.nt)
    Pnum = zeros(Float64, params.nt)
    # BC definition
    BC = BoundaryConditions()
    if params.BC_type==1
        ## Pure shear
        BC.Vx.type_W = 1;     BC.Vx.Dir_W =-params.Ebg*domain.xmin*ones(ncy+2)
        BC.Vx.type_E = 1;     BC.Vx.Dir_E =-params.Ebg*domain.xmax*ones(ncy+2)
        BC.Vx.type_S =22;     BC.Vx.Dir_S = 0.0*ones(ncx+1)
        BC.Vx.type_N =22;     BC.Vx.Dir_N = 0.0*ones(ncx+1)
        BC.Vy.type_W =22;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E =22;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S = params.Ebg*domain.ymin*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = params.Ebg*domain.ymax*ones(ncx+2)
    elseif params.BC_type==2
        ## Simple shear
        BC.Vx.type_W = 1;     BC.Vx.Dir_W = 0.0*ones(ncy+2)
        BC.Vx.type_E = 1;     BC.Vx.Dir_E =-0.0*ones(ncy+2)
        BC.Vx.type_S =11;     BC.Vx.Dir_S =-1.0*ones(ncx+1)
        BC.Vx.type_N =11;     BC.Vx.Dir_N = 1.0*ones(ncx+1)
        BC.Vy.type_W =22;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E =22;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S =-0.0*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = 0.0*ones(ncx+2)
    elseif params.BC_type==3
        ## Simple shear periodic in x
        BC.periodix  = 1
        BC.Vx.type_W = 0;     BC.Vx.Dir_W = 0.0*ones(ncy+2)
        BC.Vx.type_E = 0;     BC.Vx.Dir_E =-0.0*ones(ncy+2)
        BC.Vx.type_S =11;     BC.Vx.Dir_S =-1.0*ones(ncx+1)
        BC.Vx.type_N =11;     BC.Vx.Dir_N = 1.0*ones(ncx+1)
        BC.Vy.type_W = 0;     BC.Vy.Dir_W = 0.0*ones(ncy+1)
        BC.Vy.type_E = 0;     BC.Vy.Dir_E = 0.0*ones(ncy+1)
        BC.Vy.type_S = 1;     BC.Vy.Dir_S =-0.0*ones(ncx+2)
        BC.Vy.type_N = 1;     BC.Vy.Dir_N = 0.0*ones(ncx+2)
    end
    # Allocate tables
    fields           = Fields2D()
    fields.etac      = zeros(Float64, ncx  , ncy  )
    fields.etav      = zeros(Float64, ncx+1, ncy+1)
    fields.Kc        = zeros(Float64, ncx  , ncy  )
    fields.rhoc      = zeros(Float64, ncx  , ncy  )
    fields.Pc0       = zeros(Float64, ncx+0, ncy+0)
    fields.Pc        = zeros(Float64, ncx+0, ncy+0)
    fields.Vx        = zeros(Float64, ncx+1, ncy+2) # !!! GHOST ROWS
    fields.Vy        = zeros(Float64, ncx+2, ncy+1) # !!! GHOST COLUMNS
    fields.Damc      = ones(Float64, ncx  , ncy  )
    fields.Damv      = ones(Float64, ncx+1, ncy+1)
    Tau              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    Eps              = Tensor2D( zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+1, ncy+1), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0), zeros(Float64, ncx+0, ncy+0) )
    if BC.periodix==0
        fields.Fx    = zeros(Float64, ncx+1, ncy+0)
    else
        fields.Fx    = zeros(Float64, ncx+0, ncy+0)
    end
    fields.Fy        = zeros(Float64, ncx+0, ncy+1)
    fields.Fp        = zeros(Float64, ncx  , ncy  )
    # Initialise markers
    p                = Markers()
    p.nmx            = 4            # 2 marker per cell in x
    p.nmy            = 4            # 2 marker per cell in y
    InitialiseMarkerPositions!( p, domain )
    # Over allocate markers
    p.phase          = zeros(Float64, p.nmark_max)
    p.T              = zeros(Float64, p.nmark_max)
    p.cellx          = zeros(Int64,   p.nmark_max)
    p.celly          = zeros(Int64,   p.nmark_max)
    # Marker 2 cell mapping
    fields.mpc           = zeros(Float64,(ncx  ,ncy  )) # markers per cell
    fields.mpc_th        = zeros(Float64,(nthreads(), ncx  ,ncy   )) # markers per cell
    fields.phase_perc    = zeros(Float64,(materials.nphase, ncx  ,ncy  )) # markers per cell
    fields.phase_perc_th = zeros(Float64,(nthreads(), materials.nphase, ncx  ,ncy   )) # markers per cell
    # Initial configuration
    SetMarkers!( p, params, domain )  # Define phase on markers
    SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
    fields.Vxc    = 0.5.*(fields.Vx[1:end-1,2:end-1] .+ fields.Vx[2:end-0,2:end-1])
    fields.Vyc    = 0.5.*(fields.Vy[2:end-1,1:end-1] .+ fields.Vy[2:end-1,2:end-0])
    # Visualisation
    viz_directory = string( "Figures_", experiment )
    if isdir( viz_directory ) == false
        mkdir( viz_directory )
    end
    path = string( "./", viz_directory, "/" )
    anim = Plots.Animation( path, String[] )
    # TIME LOOP
    for it=1:params.nt
        @printf("-------------------------------------\n")
        @printf("---------- Time step #%04d ----------\n", it)
        @printf("-------------------------------------\n")
        # History
        fields.Pc0 .= fields.Pc
        # Count markers and reseed
        @time CountMarkers3!( p, fields, materials, domain )
        @time ReseedMarkers2!( p, fields.mpc, domain)
        nmark = p.nmark
        @time CountMarkers3!( p, fields, materials, domain )
        # Define dt
        if params.dt_var == 1
            maxVx     = maximum(abs.(fields.Vx))
            maxVy     = maximum(abs.(fields.Vy))
            params.dt = params.Courant*(maximum([domain.dx,domain.dy])/maximum([maxVx,maxVy]))
            @printf("dt = %2.2e\n", params.dt)
        end
        # Update cell info on markers
        @time LocateMarkers2!(p, domain)
        # Number equations
        println("Numbering")
        @time NumberingStokes!( fields, BC, domain )
        # Non-linear iterations
        for iter=1:params.niter_nl

            if iter>10  params.JFNK =1 end
            @printf("-------------------------------------\n")
            @printf("---------- Iteration #%04d ----------\n", iter)
            @printf("-------------------------------------\n")

            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            @time Rheology!( fields, Eps, Tau, params, materials, BC, ncx, ncy )
            @time Stress!( Tau, Eps, fields )
            @time ResidualsComp!( fields, Tau, Eps, BC, domain, params )
            println("|Fx| = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy| = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp| = ", norm(fields.Fp)/length(fields.Fp))
            if ( (norm(fields.Fx)/length(fields.Fx) < params.tol_nl) & (norm(fields.Fy)/length(fields.Fy) < params.tol_nl) & (norm(fields.Fp)/length(fields.Fp) < params.tol_nl) )
                println("Mechanical solver converged")
                break
            end

            # Assemble Stokes matrices
            println("Assembly")
            @time Kuu, Kup, Kpu = StokesAssembly( BC, fields, params.Dir_scale, domain )
            # Call solver
            println("Solver")
            @time StokesSolver!( fields, params, domain, materials, BC, Kuu, Kup, Kpu )

            # Evaluate residuals
            println("Residuals")
            @time SetBCs( fields.Vx, fields.Vy, BC )
            @time StrainRate!( Eps, fields.Vx, fields.Vy, domain )
            @time Stress!( Tau, Eps, fields )
            @time ResidualsComp!( fields, Tau, Eps, BC, domain, params )
            println("|Fx| = ", norm(fields.Fx)/length(fields.Fx))
            println("|Fy| = ", norm(fields.Fy)/length(fields.Fy))
            println("|Fp| = ", norm(fields.Fp)/length(fields.Fp))
        end
        params.t += params.dt
        # Advect particles
        if params.advection == 1
            @time RungeKutta2!( p, fields, params, BC, domain )
            # Deform box
            if params.ALE==1
                dx, dy = PureShearBoxUpdate!( domain, BC, params )
                SetInitialVelocity2!( fields.Vx, fields.Vy, BC, domain )
            end
        end
        # Postprocessing
        # Visualisation
        if params.show_figs==1 && ( mod( it, params.nout )==0 || it==1 )
            p1 = Plots.heatmap(domain.xc,  domain.yc, Array(Tau.II)', aspect_ratio=1, xlims=(domain.xmin, domain.xmax), ylims=(domain.ymin, domain.ymax), c=Plots.cgrad(:roma, rev = true), title="Tii")
            display(Plots.plot( p1, dpi=200 ) ); Plots.frame(anim)
        end
    end
    if params.show_figs==1 Plots.gif(anim, string( path, experiment, ".gif" ), fps = 15) end
    return
end

for it=1:1
   @time main()
end


text:
export Backend, CPUBackend, AbstractGPUBackend
export init, shutdown, registry_reset, registry_put

@compat abstract type Backend end
const ParameterRegistry = Dict{AbstractString, Vector{AbstractParameter}}

import Base.show
export show
show(io::IO, backend::Backend) = show(io, typeof(backend))

function init(backend::Backend)
end
function shutdown(backend::Backend)
  registry_reset(backend)
end
function registry_reset(backend::Backend)
  for (k,params) in backend.param_registry
    map(destroy, params)
  end
  backend.param_registry = ParameterRegistry()
end
function registry_put(backend::Backend, key::AbstractString, params::Vector)
  if haskey(backend.param_registry, key)
    map(destroy, backend.param_registry[key])
  end

  # we keep a reference to the parameters, so that even those the original
  # network is destroyed, we still have valid access to those trained parameters
  backend.param_registry[key] = AbstractParameter[share_parameter(backend, p) for p in params]
end
function registry_get(backend::Backend, key::AbstractString)
  return get(backend.param_registry, key, nothing)
end

type CPUBackend <: Backend
  param_registry :: ParameterRegistry

  CPUBackend() = new(ParameterRegistry())
end

# This is forward declaration to allow some code to compile
# (especially testing codes) even if CUDA module is completely
# disabled. See test/layers/pooling.jl for example.
@compat abstract type AbstractGPUBackend <: Backend end


text:
using ProximityEnrichment
using Test
using GenomicFeatures

intA = Interval("chr1", 10, 100, '.')
intB = Interval("chr1", 80, 200, '.')
intC = Interval("chr1", 110, 500, '.')
intD = Interval("chr2", 100, 200, '.')
intE = Interval("chr2", 400, 600, '.')
intF = Interval("chrX", 10, 600, '-')

ivs = [intA, intB, intC, intD, intE, intF]

@testset "ProximityEnrichment.jl" begin
    # Write your own tests here
    @test ProximityEnrichment.intervaldist(intA, intB) == 0
    @test ProximityEnrichment.intervaldist(intB, intA) == 0
    @test ProximityEnrichment.intervaldist(intA, intC) == 10
    @test ProximityEnrichment.intervaldist(intC, intA) == 10
    @test ProximityEnrichment.intervaldist(intB, intC) == 0

    @test ProximityEnrichment.intervaldist(Interval("chr1", 1:20, '.'), Interval("chr2", 1:20, '.')) == typemax(Int)
    @test ProximityEnrichment.intervaldist(Interval("chr1", 1:20,), Interval("chr1", 100:120)) == 80

    @test ProximityEnrichment.closestinterval(ivs, Interval("chr1", 1, 8, '.')) == 2
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr1", 1, 9, '.')) == 1
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr1", 5, 15, '.')) == 0
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr1", 50, 150, '.')) == 0

    @test ProximityEnrichment.closestinterval(ivs, Interval("chr1", 600, 750, '.')) == 100
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr2", 10, 15, '.')) == 85
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr2", 300, 300, '.')) == 100
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr2", 300, 300, '.')) == 100
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr2", 1000, 2000, '+')) == 400
    @test ProximityEnrichment.closestinterval(ivs, Interval("chr3", 1000, 2000, '+')) == -1
    @test ProximityEnrichment.closestinterval(ivs, Interval("chrY", 1000, 2000, '+')) == -1

end

intA = Interval("chr1", 10, 100, '.')
intB = Interval("chr1", 110, 500, '.')
intC = Interval("chr1", 80, 200, '.')

@show ProximityEnrichment.intervaldist(intA, intB)

max_stars_repo_name: owensnick/ProximityEnrichment

max_issues_repo_name: owensnick/ProximityEnrichment

max_forks_repo_name: owensnick/ProximityEnrichment


text:
#-----------------------------------------------------------------------
# JuMPeR  --  JuMP Extension for Robust Optimization
# http://github.com/IainNZ/JuMPeR.jl
#-----------------------------------------------------------------------
# Copyright (c) 2016: Iain Dunning
# This Source Code Form is subject to the terms of the Mozilla Public
# License, v. 2.0. If a copy of the MPL was not distributed with this
# file, You can obtain one at http://mozilla.org/MPL/2.0/.
#-----------------------------------------------------------------------
# example/portfolio.jl
# Solve a robust portfolio optimization problem. Can choose between a
# polyhedral set (requires a LP solver) or an ellipsoidal set (requires
# a SOCP solver).
#-----------------------------------------------------------------------

using JuMP, JuMPeR  # Modeling
using Distributions, Random  # For generating data
using LinearAlgebra
using Gurobi
using Printf # for pretty printing.

# Number of stocks
const NUM_ASSET = 10
# Number of samples of past returns
const NUM_SAMPLES = 1000
# Uncertainty set type (:Polyhedral or :Ellipsoidal)
#const UNCERTAINY_SET = :Polyhedral
const UNCERTAINY_SET = :Ellipsoidal
# Seed the RNG for reproducibility
Random.seed!(1234)

"""
    generate_data()

Sample returns of the assets from a multivariate normal distribution.
Returns a matrix with the samples in rows & each column is an asset.
"""
function generate_data()
    data = zeros(NUM_SAMPLES, NUM_ASSET)
    # Linking factors to induce correlations
    β = [(i-1.0)/NUM_ASSET for i in 1:NUM_ASSET]
    for i in 1:NUM_SAMPLES
        # Common market factor, μ=3%, σ=5%, truncated at ±3σ
        z = rand(Normal(0.03, 0.05))
        z = max(z, 0.03 - 3*0.05)
        z = min(z, 0.03 + 3*0.05)
        for j in 1:NUM_ASSET
            # Idiosyncratic contribution, μ=0%, σ=5%, truncated at ±3σ
            r = rand(Normal(0.00, 0.05))
            r = max(r, 0.00 - 3*0.05)
            r = min(r, 0.00 + 3*0.05)
            data[i,j] = β[j] * z + r
        end
    end
    return data
end

"""
    solve_portfolio(past_returns, Γ)

Solve the robust portfolio problem given a matrix of past returns.
The particular problem solved is

    max   z
    s.t.  ∑ rᵢ xᵢ ≥ z    ∀ r ∈ R      # Worst-case return
          ∑ xᵢ    = 1                 # Is a portfolio
          xᵢ ≥ 0,

where R is one of the uncertainty sets:

    R = { (r,z) | r = μ + Σ^½ z,
            Polyhedral:  ‖z‖₁ ≤ Γ, ‖z‖∞ ≤ 1
           Ellipsoidal:  ‖z‖₂ ≤ Γ
        }

which is controlled by the global flag UNCERTAINY_SET.
"""
function solve_portfolio(past_returns, Γ)
    # Setup the robust optimization model
    m = RobustModel(solver = GurobiSolver())
    # Each asset is a share of the money to invest...
    @variable(m, 0 <= x[1:NUM_ASSET] <= 1)
    # ... and we must allocate all the money.
    @constraint(m, sum(x) == 1)
    # JuMPeR doesn't support uncertain objective functions directly. Instead,
    # an auxiliary variable should be added and the objective function should
    # be placed as a constraint on the auxiliary (epigraph form).
    @variable(m, obj)
    @objective(m, Max, obj)
    # The uncertain parameters are the returns for each asset
    @uncertain(m, r[1:NUM_ASSET])
    # The uncertainty set requires the "square root" of the covariance
    μ = vec(mean.(past_returns[:,i] for i=1:size(past_returns,2)))  # Want a column vector
    Σ = cov(past_returns)
    L = Matrix(cholesky(Σ))  # Σ^½
    # Define auxiliary uncertain parameters to model the underlying factors
    @uncertain(m, z[1:NUM_ASSET])
    # Link r with z
    @constraint(m, r .== L*z + μ)
    if UNCERTAINY_SET == :Polyhedral
        @constraint(m, norm(z,   1) <= Γ)  # ‖z‖₁ ≤ Γ
        @constraint(m, norm(z, Inf) <= 1)  # ‖z‖∞ ≤ 1
    else
        @constraint(m, norm(z,   2) <= Γ)  # ‖z‖₂ ≤ Γ
    end
    # The objective function: the worst-case return
    @constraint(m, obj <= dot(r, x))
    # Solve it!
    solve(m)
    # Return the allocation and the worst-case return
    return getvalue(x), getvalue(obj)
end

"""
    solve_and_simulate()

Generates past and future returns from the same distribution, builds a
range of portfolios using the past returns, and evaluates them on the future
returns (reporting distribution of returns).
"""
function solve_and_simulate()
    # Generate the simulated data
    past_returns   = generate_data()
    future_returns = generate_data()
    # Generate and store the portfolios
    portfolios = Dict()
    for Γ in 0:NUM_ASSET
        x, _ = solve_portfolio(past_returns, Γ)
        portfolios[Γ] = x
    end
    # Print table header
    println("    Γ |    Min |    10% |    20% |   Mean |    Max | Support")
    # Evaluate and print distributions
    for Γ in 0:NUM_ASSET
        # Evaluate portfolio returns in future
        future_z = future_returns * portfolios[Γ]
        sort!(future_z)
        min_z    = future_z[1]*100
        ten_z    = future_z[div(NUM_SAMPLES,10)]*100
        twenty_z = future_z[div(NUM_SAMPLES, 5)]*100
        mean_z   = mean(future_z)*100
        max_z    = future_z[end]*100
        support  = sum(portfolios[Γ] .> 0.01)  # Number of assets used in portfolio
        @printf(" %4.1f | %6.2f | %6.2f | %6.2f | %6.2f | %6.2f | %7d\n",
                    Γ, min_z, ten_z, twenty_z, mean_z, max_z, support)
    end
end

solve_and_simulate()


text:
import Base.push!  # for DisjointMinSets

"""
    label = label_components(A; bkg = zero(eltype(A)), dims=coords_spatial(A))
    label = label_components(A, connectivity; bkg = zero(eltype(A)))

Find the connected components in an array `A`. Components are defined as
connected voxels that all have the same value distinct from `bkg`, which
corresponds to the "background" component.

Specify connectivity in one of three ways:

- A list indicating which dimensions are used to determine
  connectivity. For example, `dims = (1,3)` would not test neighbors along
  dimension 2 for connectivity. This corresponds to just the nearest neighbors,
  i.e., default 4-connectivity in 2d and 6-connectivity in 3d.

- An iterable `connectivity` object with `CartesianIndex` elements encoding the
  displacement of each checked neighbor.

- A symmetric boolean array of the same dimensionality as `A`, of size 1 or 3
  along each dimension. Each entry in the array determines whether a given
  neighbor is used for connectivity analyses. For example, in two dimensions
  `connectivity = trues(3,3)` would include all pixels that touch the
  current one, even the corners.

The output `label` is an integer array, where `bkg` elements get a value of 0.

# Examples

```jldoctest; setup=:(using ImageMorphology)
julia> A = [true false false true  false;
            true false true  true  true]
2×5 Matrix{Bool}:
 1  0  0  1  0
 1  0  1  1  1

julia> label_components(A)
2×5 Matrix{$Int}:
 1  0  0  2  0
 1  0  2  2  2

julia> label_components(A; dims=2)
2×5 Matrix{$Int}:
 1  0  0  4  0
 2  0  3  3  3
```
With `dims=2`, entries in `A` are connected if they are in the same row, but
not if they are in the same column.
"""
label_components(A::AbstractArray; bkg=zero(eltype(A)), dims=coords_spatial(A)) =
    label_components(A, half_diamond(A, dims); bkg=bkg)
label_components(A::AbstractArray, connectivity::AbstractArray{Bool}; bkg=zero(eltype(A))) =
    label_components(A, half_pattern(A, connectivity); bkg=bkg)
label_components(A::AbstractArray, iter; bkg=zero(eltype(A))) =
    label_components!(similar(A, Int), A, iter; bkg=bkg) # Int is a safe choice

label_components!(out::AbstractArray{<:Integer}, A::AbstractArray; bkg=zero(eltype(A)), dims=coords_spatial(A)) =
    label_components!(out, A, half_diamond(A, dims); bkg=bkg)
label_components!(out::AbstractArray{<:Integer}, A::AbstractArray, connectivity::AbstractArray{Bool}; bkg=zero(eltype(A))) =
    label_components!(out, A, half_pattern(A, connectivity); bkg=bkg)
function label_components!(out::AbstractArray{T}, A::AbstractArray, iter; bkg=zero(eltype(A))) where T<:Integer
    axes(out) == axes(A) || throw_dmm(axes(out), axes(A))
    fill!(out, zero(T))
    sets = DisjointMinSets{T}()
    sizehint!(sets.parents, floor(Int, sqrt(length(A))))
    @inbounds for i in CartesianIndices(A)
        val = A[i]
        val == bkg && continue
        label = typemax(T)    # sentinel value
        for Δi in iter
            ii = i + Δi
            checkbounds(Bool, A, ii) || continue
            if A[ii] == val
                newlabel = out[ii]
                label = ((label == typemax(T)) | (label == newlabel)) ? newlabel : union!(sets, label, newlabel)
            end
        end
        if label == typemax(T)
            label = push!(sets)
        end
        out[i] = label
    end
    # Now parse sets to find the labels
    newlabel = minlabel(sets)
    @inbounds for i in eachindex(A, out)
        if A[i] != bkg
            out[i] = newlabel[find_root!(sets, out[i])]
        end
    end
    return out
end

throw_dmm(ax1, ax2) = throw(DimensionMismatch("axes of input and output must match, got $ax1 and $ax2"))

function half_diamond(A::AbstractArray{T,N}, dims) where {T,N}
    offsets = CartesianIndex{N}[]
    for d = 1:N
        if d ∈ dims
            push!(offsets, CartesianIndex(ntuple(i -> i == d ? -1 : 0, N)))
        end
    end
    return (offsets...,)   # returning as a tuple allows specialization
    # return offsets
end
half_diamond(A::AbstractArray{T,N}, ::Colon) where {T,N} = half_diamond(A, 1:N)

function half_pattern(A::AbstractArray{T,N}, connectivity::AbstractArray{Bool}) where {T,N}
    all(in((1,3)), size(connectivity)) || throw(ArgumentError("connectivity must have size 1 or 3 in each dimension"))
    for d = 1:ndims(connectivity)
        size(connectivity, d) == 1 || reverse(connectivity; dims=d) == connectivity || throw(ArgumentError("connectivity must be symmetric"))
    end
    center = CartesianIndex(map(axes(connectivity)) do ax
        (first(ax) + last(ax)) ÷ 2
    end)
    offsets = CartesianIndex{N}[]
    for i in CartesianIndices(connectivity)
        i == center && break   # we only need the ones that come prior
        if connectivity[i]
            push!(offsets, i - center)
        end
    end
    # return offsets
    return (offsets...,)   # returning as a tuple allows specialization
end

# Copied directly from DataStructures.jl, but specialized
# to always make the parent be the smallest label
struct DisjointMinSets{T<:Integer}
    parents::Vector{T}

    DisjointMinSets{T}(n::Integer) where T = new([T(1):T(n);])
end
DisjointMinSets{T}() where T<:Integer = DisjointMinSets{T}(T(0))
DisjointMinSets() = DisjointMinSets{INt}()

Base.@propagate_inbounds function find_root!(sets::DisjointMinSets, m::Integer)
    p = sets.parents[m]   # don't use @inbounds here, it might not be safe
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    p
end

# an unsafe variant of the above
function find_root_unsafe!(sets::DisjointMinSets, m::Int)
    @inbounds p = sets.parents[m]
    @inbounds if sets.parents[p] != p
        sets.parents[m] = p = find_root_unsafe!(sets, p)
    end
    p
end

Base.@propagate_inbounds function union!(sets::DisjointMinSets, m::Integer, n::Integer)
    mp = find_root!(sets, m)
    np = find_root!(sets, n)
    if mp < np
        sets.parents[np] = mp
        return mp
    elseif np < mp
        sets.parents[mp] = np
        return np
    end
    mp
end

function push!(sets::DisjointMinSets)
    m = length(sets.parents) + 1
    m >= typemax(eltype(sets.parents)) && error("labels exhausted, use a larger integer type")
    push!(sets.parents, m)
    m
end

function minlabel(sets::DisjointMinSets)
    out = Vector{Int}(undef, length(sets.parents))
    k = 0
    @inbounds for i = 1:length(sets.parents)
        if sets.parents[i] == i
            k += 1
        end
        out[i] = k
    end
    out
end

"`component_boxes(labeled_array)` -> an array of bounding boxes for each label, including the background label 0"
function component_boxes(img::AbstractArray{Int})
    nd = ndims(img)
    n = [Vector{Int}[ fill(typemax(Int),nd), fill(typemin(Int),nd) ]
            for i=0:maximum(img)]
    s = CartesianIndices(size(img))
    for i=1:length(img)
        vcur = s[i]
        vmin = n[img[i]+1][1]
        vmax = n[img[i]+1][2]
        for d=1:nd
            vmin[d] = min(vmin[d], vcur[d])
            vmax[d] = max(vmax[d], vcur[d])
        end
    end
    map(x->map(y->tuple(y...),x),n)
end

"`component_lengths(labeled_array)` -> an array of areas (2D), volumes (3D), etc. for each label, including the background label 0"
function component_lengths(img::AbstractArray{Int})
    n = zeros(Int,maximum(img)+1)
    for i=1:length(img)
        n[img[i]+1]+=1
    end
    n
end

"`component_indices(labeled_array)` -> an array of pixels for each label, including the background label 0"
function component_indices(img::AbstractArray{Int})
    n = [Int[] for i=0:maximum(img)]
    for i=1:length(img)
      push!(n[img[i]+1],i)
    end
    n
end

"`component_subscripts(labeled_array)` -> an array of pixels for each label, including the background label 0"
function component_subscripts(img::AbstractArray{Int})
    n = [Tuple[] for i=0:maximum(img)]
    s = CartesianIndices(size(img))
    for i=1:length(img)
      push!(n[img[i]+1],s[i])
    end
    n
end

"`component_centroids(labeled_array)` -> an array of centroids for each label, including the background label 0"
function component_centroids(img::AbstractArray{Int,N}) where N
    len = length(0:maximum(img))
    n = fill(zero(CartesianIndex{N}), len)
    counts = fill(0, len)
    @inbounds for I in CartesianIndices(size(img))
        v = img[I] + 1
        n[v] += I
        counts[v] += 1
    end
    map(v -> n[v].I ./ counts[v], 1:len)
end


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

# %% [markdown]
# # Absorbing boundary condition - 1D case
#
# Gen Kuroki
#
# 2021-06-21
#
# 1D Wave equation on the interval $[a, b]$ with absorbing boundary condition:
#
# $$
# \begin{aligned}
# &
# \frac{\partial^2 u}{\partial t^2} = \frac{\partial^2 u}{\partial x^2} \quad (a\le x \le b),
# \\ &
# \left.\left(\frac{\partial}{\partial t} - \frac{\partial}{\partial x}\right)u\;\right|_{x = a} = 0,
# \\ &
# \left.\left(\frac{\partial}{\partial t} + \frac{\partial}{\partial x}\right)u\;\right|_{x = b} = 0.
# \end{aligned}
# $$
#
# The second equation means backward wave only and the third means forward wave only.
#
# Its discretization: $dx = (b-a)/n$, $x = a + k\,dx$, $k=0,1,\ldots,n$,
#
# $$
# \begin{aligned}
# &
# u_{tt}(t, x) = \frac{u(t, x - dx) + u(t, x + dx) - 2u(t, x)}{dx^2} \quad (a\le x\le b),
# \\ &
# u_t(t, a) - \frac{u(t, a+dx) - u(t, a-dx)}{2dx} = 0,
# \\ &
# u_t(t, b) + \frac{u(t, b+dx) - u(t, b-dx)}{2dx} = 0.
# \end{aligned}
# $$
#
# The second and third equations are equivalent to
#
# $$
# \begin{aligned}
# &
# u(t, a-dx) = u(t, a+dx) - 2u_t(t, a)\,dx,
# \\ &
# u(t, b+dx) = u(t, b-dx) - 2u_t(t, b)\,dx.
# \end{aligned}
# $$
#
# Therefore the discretization above can be rewritten as
#
# $$
# \begin{aligned}
# &
# u_{tt}(t, x) = \frac{u(t, x - dx) + u(t, x + dx) - 2u(t, x)}{dx^2} \quad (a < x < b),
# \\ &
# u_{tt}(t, a) = 2\left(\frac{u(t, a + dx) - u(t, a)}{dx^2} - \frac{u_t(t, a)}{dx}\right),
# \\ &
# u_{tt}(t, b) = 2\left(\frac{u(t, b - dx) - u(t, b)}{dx^2} - \frac{u_t(t, b)}{dx}\right).
# \end{aligned}
# $$

# %%
@show VERSION
@time using DifferentialEquations
@time using ForwardDiff
@time using Plots

# %%
function f2_abc!(dv, v, u, p, t)
    (; dx) = p
    a, b = firstindex(u), lastindex(u)
    @. @views dv[a+1:b-1] = (u[a:b-2] + u[a+2:b] - 2u[a+1:b-1])/dx^2
    dv[a] = 2((u[a+1] - u[a])/dx^2 - v[a]/dx)
    dv[b] = 2((u[b-1] - u[b])/dx^2 - v[b]/dx)
    return
end

# %%
n = 20
x = range(-10, 10; length=n+1)
dx = step(x)
p = (; dx)

U(t, x) = 2/3*exp(-(x - t + 5)^2) + 1/3*exp(-(x + t - 5)^2)
V(t, x) = ForwardDiff.derivative(t -> U(t, x), t)
u0 = U.(0, x)
v0 = V.(0, x)
tspan = (0.0, 40.0)

prob = SecondOrderODEProblem(f2_abc!, v0, u0, tspan, p)
sol = solve(prob)

ts = range(sol.prob.tspan...; length=200)
anim = @animate for t in [fill(ts[begin], 10); ts; fill(ts[end], 10)]
    plot(x, sol(t)[end÷2+1:end]; label="numerical n = $n", ylim=(-0.4, 1.0), size=(600, 300))
    xs = range(extrema(x)...; length=201)
    plot!(xs, U.(t, xs); label="exact solution", ls=:dash)
end
gif(anim, "1d_wave_eq_abc_21.gif")

# %%
n = 40
x = range(-10, 10; length=n+1)
dx = step(x)
p = (; dx)

U(t, x) = 2/3*exp(-(x - t + 5)^2) + 1/3*exp(-(x + t - 5)^2)
V(t, x) = ForwardDiff.derivative(t -> U(t, x), t)
u0 = U.(0, x)
v0 = V.(0, x)
tspan = (0.0, 40.0)

prob = SecondOrderODEProblem(f2_abc!, v0, u0, tspan, p)
sol = solve(prob)

ts = range(sol.prob.tspan...; length=200)
anim = @animate for t in [fill(ts[begin], 10); ts; fill(ts[end], 10)]
    plot(x, sol(t)[end÷2+1:end]; label="numerical n = $n", ylim=(-0.2, 1.0), size=(600, 300))
    xs = range(extrema(x)...; length=201)
    plot!(xs, U.(t, xs); label="exact solution", ls=:dash)
end
gif(anim, "1d_wave_eq_abc_41.gif")

# %%
n = 100
x = range(-10, 10; length=n+1)
dx = step(x)
p = (; dx)

U(t, x) = 2/3*exp(-(x - t + 5)^2) + 1/3*exp(-(x + t - 5)^2)
V(t, x) = ForwardDiff.derivative(t -> U(t, x), t)
u0 = U.(0, x)
v0 = V.(0, x)
tspan = (0.0, 40.0)

prob = SecondOrderODEProblem(f2_abc!, v0, u0, tspan, p)
sol = solve(prob)

ts = range(sol.prob.tspan...; length=200)
anim = @animate for t in [fill(ts[begin], 10); ts; fill(ts[end], 10)]
    plot(x, sol(t)[end÷2+1:end]; label="numerical n = $n", ylim=(-0.1, 1.1), size=(600, 250))
    xs = range(extrema(x)...; length=201)
    plot!(xs, U.(t, xs); label="exact solution", ls=:dash)
end
gif(anim, "1d_wave_eq_abc_101.gif")

# %%
n = 200
x = range(-10, 10; length=n+1)
dx = step(x)
p = (; dx)

U(t, x) = 2/3*exp(-(x - t + 5)^2) + 1/3*exp(-(x + t - 5)^2)
V(t, x) = ForwardDiff.derivative(t -> U(t, x), t)
u0 = U.(0, x)
v0 = V.(0, x)
tspan = (0.0, 20.0)

prob = SecondOrderODEProblem(f2_abc!, v0, u0, tspan, p)
sol = solve(prob)

ts = range(sol.prob.tspan...; length=100)
anim = @animate for t in [fill(ts[begin], 10); ts; fill(ts[end], 10)]
    plot(x, sol(t)[end÷2+1:end]; label="numerical n = $n", ylim=(-0.1, 1.1), size=(600, 250))
    plot!(x, U.(t, x); label="exact solution", ls=:dash)
end
gif(anim, "1d_wave_eq_abc_201.gif")

# %%

max_stars_repo_path: 0006/Absorbing boundary condition - 1D case.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
module Cursor
    # note that  \u1b  is in decimal 033, which is an escape code
    const ERASE_TO_EOL = "\u1b[K"
    const UP1 = "\u1b[A"
    up(n) = "\u1b[$(n)A"
    const MOVE_TO_SOL = "\r" #Carriage Return
end

function move_cursor_up_while_clearing_lines(io, numlinesup)
    for _ in 1:numlinesup
        print(io, Cursor.MOVE_TO_SOL * Cursor.ERASE_TO_EOL * Cursor.UP1)
    end
    print(io, Cursor.MOVE_TO_SOL * Cursor.ERASE_TO_EOL * Cursor.MOVE_TO_SOL)
end


text:
export ilan

using IterativeSolvers
using LinearAlgebra
using Random
using Statistics

# Types specifying which way to B_prod in chebyshev iar
abstract type compute_Bmul_method end;
abstract type compute_Bmul_method_DEP <: compute_Bmul_method end;
abstract type compute_Bmul_method_SPMF_NEP <: compute_Bmul_method end;
abstract type compute_Bmul_method_DerSPMF <: compute_Bmul_method end;
abstract type compute_Bmul_method_Auto <: compute_Bmul_method end;

# Data collected in a precomputation phase.
abstract type IlanAbstractPrecomputeData end
mutable struct IlanPrecomputeDataDEP <: IlanAbstractPrecomputeData
     G; vv; ZZ; QQ; QQ2
end
mutable struct IlanPrecomputeDataSPMF_NEP <: IlanAbstractPrecomputeData
     G; FDH; ZZ; QQ
end
mutable struct IlanPrecomputeDataDerSPMF <: IlanAbstractPrecomputeData
     G; FDH; ZZ; QQ
end

"""
    ilan(nep,[maxit=30,][σ=0,][γ=1,][linsolvecreator=DefaultLinSolverCreator(),][tolerance=eps()*10000,][neigs=6,][errmeasure,][v=rand(size(nep,1),1),][logger=0,][check_error_every=30,][orthmethod=DGKS])

Run the infinite Lanczos method on the symmetric nonlinear eigenvalue problem stored in `nep`. The current implementation supports only `nep`s in `SPMF` format.

The target `σ` is the center around which eiganvalues are computed. The kwarg `errmeasure` is a function handle which can be used to specify how the error is measured to be used in termination (default is absolute residual norm). A Ritz pair `λ` and `v` is flagged a as converged (to an eigenpair) if `errmeasure` is less than `tol`. The vector `v` is the starting vector for constructing the Krylov space. The orthogonalization method, used in contructing the orthogonal basis of the Krylov space, is specified by `orthmethod`, see the package `IterativeSolvers.jl`. The iteration is continued until `neigs` Ritz pairs have converged.
This function throws a `NoConvergenceException` if the wanted eigenpairs are not computed after `maxit` iterations. However, if `neigs` is set to `Inf` the iteration is continued until `maxit` iterations without an error being thrown.

See [`augnewton`](@ref) for other parameters.

# Example
```julia-repl
julia> using NonlinearEigenproblems, LinearAlgebra
julia> nep=nep_gallery("dep_symm_double",10);
julia> v0=ones(size(nep,1));
julia> λ,v=ilan(nep;v=v0,tol=1e-5,neigs=3);
julia> norm(compute_Mlincomb!(nep,λ[1],v[:,1])) # Is it an eigenvalue?
julia> λ    # print the computed eigenvalues
3-element Array{Complex{Float64},1}:
  0.03409997385842267 - 1.425205509434608e-19im
 -0.03100798730589012 + 5.66201058098364e-20im
  -0.0367653644764646 - 1.607494907684445e-19im
```

# References
* Algorithm 2 in Mele, The infinite Lanczos method for symmetric nonlinear eigenvalue problems, https://arxiv.org/abs/1812.07557, 2018
"""
ilan(nep::NEP;params...)=ilan(ComplexF64,nep;params...)
function ilan(
    ::Type{T},
    nep::NEP;
    orthmethod::Type{T_orth}=DGKS,
    maxit=30,
    linsolvercreator=DefaultLinSolverCreator(),
    tol=eps(real(T))*10000,
    neigs=6,
    errmeasure::ErrmeasureType = DefaultErrmeasure(nep),
    σ=zero(T),
    γ=one(T),
    v=randn(real(T),size(nep,1)),
    logger=0,
    check_error_every=30,
    inner_solver_method=DefaultInnerSolver(),
    compute_Bmul_method::Type{T_y0}=compute_Bmul_method_Auto,
    proj_solve=true,
    inner_logger=0
    )where{T<:Number,T_orth<:IterativeSolvers.OrthogonalizationMethod,T_y0<:compute_Bmul_method}

    @parse_logger_param!(logger)
    @parse_logger_param!(inner_logger)

    # Ensure types σ and v are of type T
    σ=T(σ)
    v=Array{T,1}(v)

    if (compute_Bmul_method == compute_Bmul_method_Auto)
        if (isa(nep,DEP))
            compute_Bmul_method=compute_Bmul_method_DEP
        elseif (isa(nep,DerSPMF))
            compute_Bmul_method=compute_Bmul_method_DerSPMF
        elseif (isa(nep,SPMF_NEP))
            compute_Bmul_method=compute_Bmul_method_SPMF_NEP
        else
            println("Error in compute_Bmul_method. You provided and invalid `nep` for which `compute_Bmul_method` is not automatically supported. Provide the `nep` in a compatible format: `DEP`, `SPMF_NEP` or `DerSPMF`. Otherwise you can overload the function `Bmult` for your specif `nep`.")
        end
    end

    n = size(nep,1);
    m = maxit;

    # initialization
    V=zeros(T,n,m+1)
    Q=zeros(T,n,m+1)
    Qp=zeros(T,n,m+1)
    Qn=zeros(T,n,m+1)
    Z=zeros(T,n,m+1)
    H=zeros(T,m+1,m)
    HH=zeros(T,m+1,m)
    ω=zeros(T,m+1)
    a=Vector{T}(γ.^(0:2m+2)); a[1]=zero(T);
    local M0inv::LinSolver = create_linsolver(linsolvercreator,nep,σ)
    err=NaN*ones(m,m);
    W=zeros(T,n,m+1);

    # allocate extra memory for storing the first blocks of the Arnoldi sequence
    # only if the extraction of the eigenpairs is done with the Ritz pairs
    if !proj_solve QQ=zeros(T,n,m+1) end

    # precomputation for exploiting the structure DEP, PEP, GENERAL
    precomp=precompute_data(T,nep,compute_Bmul_method,n,m,σ,γ)

    # setting initial step
    Q[:,1]=v/norm(v)
    ω[1]=Q[:,1]⋅compute_Mlincomb(nep,0,hcat(Q[:,1],Q[:,1]),[0,1]);
    V[:,1]=Q[:,1];

    k=1; conv_eig=0;
    Av=get_Av(nep)

    while (k <= m) && (conv_eig<neigs)
        # store all the (first blocks) vectors of the Arnoldi sequence. This is needed
        # for the extraction of Ritz vectors.
        if !proj_solve
            QQ[:,k]=Q[1:n,1];
        end
        broadcast!(/,view(Qn,:,2:k+1),view(Q,:,1:k),(1:k)')
        Qn[:,1] = compute_Mlincomb!(nep,σ,view(Qn,:,1:k+1),a[1:k+1]);
        Qn[:,1] .= -lin_solve(M0inv,Qn[:,1]);

        # call Bmult
        Bmult!(compute_Bmul_method,k,view(Z,:,1:k+1),Qn,Av,precomp,γ)

        # orthogonalization (three terms recurrence)
        if k>1 β=mat_sum(view(Z,:,1:k),view(Qp,:,1:k)) end
        α=mat_sum(view(Z,:,1:k),view(Q,:,1:k))
        η=mat_sum(view(Z,:,1:k+1),view(Qn,:,1:k+1))

        H[k,k]=α/ω[k]
        if k>1 H[k-1,k]=β/ω[k-1] end
        #Qn[:,1:k] .-= H[k,k]*view(Q,:,1:k);
        mul_and_sub!(view(Qn,:,1:k),view(Q,:,1:k),H[k,k])
        #if k>1 Qn[:,1:k] .-= H[k-1,k]*view(Qp,:,1:k) end
        if k>1 mul_and_sub!(view(Qn,:,1:k),view(Qp,:,1:k),H[k-1,k]) end

        H[k+1,k]=norm(Qn);
        #Qn[:,1:k+1] ./= H[k+1,k]
        #Qn[:] ./= H[k+1,k]
        scal_mul!(view(Qn,:,1:k+1), 1/H[k+1,k])
        #ldiv!(H[k+1,k],view(Qn,:,1:k+1))

        ω[k+1]=η-2*α*H[k,k]+ω[k]*H[k,k]^2;
        if k>1
            ω[k+1]=ω[k+1]-2*β*H[k-1,k]+ω[k-1]*H[k-1,k]^2;
        end
        ω[k+1]=ω[k+1]/H[k+1,k]^2;
        V[:,k+1]=Qn[:,1];

        orthogonalize_and_normalize!(view(V,:,1:k),view(V,:,k+1), view(HH,1:k,k), orthmethod)

        # extract eigenpair approximation
        if ((rem(k,check_error_every)==0)||(k==m))
            if !proj_solve

                # Extract eigenvalues from Hessenberg matrix
                D,W_Ritz = eigen(H[1:k,1:k])

                mul!(view(W,:,1:k),view(QQ,:,1:k),W_Ritz)
                #ZZ = view(QQ,:,1:k)*Z_Ritz
                λ = σ .+ γ ./ D
            else
                VV=view(V,:,1:k+1)
                # Projected solve to extract eigenvalues (otw hessenberg matrix)

                # create the projected NEP
                mm=size(VV,2)
                pnep=create_proj_NEP(nep,mm); # maxsize=mm
                set_projectmatrices!(pnep,VV,VV);

                # solve the projected NEP
                push_info!(logger,2,"Solving the projected problem",continues=true);

                λproj,Wproj=inner_solve(inner_solver_method,T,pnep;
                                        neigs=m,
                                        inner_logger=inner_logger,
                                        maxit=100,
                                        tol=tol);

                #println("λproj=",λproj)

                q=length(λproj)
                λ=λproj[1:q];

                if q>m q=m end
                mul!(view(W,:,1:q),VV,view(Wproj,:,1:q))
                #ZZ=VV*Zproj[:,1:k];

            end
            # compute the errors
            err[k,1:size(λ,1)]=map(s->estimate_error(errmeasure,λ[s],W[:,s]),1:size(λ,1))

            # Log them and compute the converged
            push_iteration_info!(logger,2, k,err=err[k,1:size(λ,1)], λ=λ,
                                 continues=true);
            conv_eig=0
            for s=1:size(λ,1)
                if err[k,s]<tol;
                    conv_eig=conv_eig+1;
                    push_info!(logger,"+", continues=true);
                elseif err[k,s]<tol*10
                    push_info!(logger,"=", continues=true);
                else
                    push_info!(logger,"-", continues=true);
                end
            end
            push_info!(logger,"");

            # Sort the errors
            idx=sortperm(err[k,1:k]); # sort the error
            err[k,1:k]=err[k,idx];

            # extract the converged Ritzpairs
            if (k==m)||(conv_eig>=neigs)
                nrof_eigs = Int(min(conv_eig,neigs))
                λ=λ[idx[1:nrof_eigs]]
                W=W[:,idx[1:length(λ)]]
            end

        end

        k=k+1;
        # shift the vectors
        Qp[:]=Q; Q[:]=Qn; Qn[:] .=0;
    end

    k=k-1

    # NoConvergenceException
    if conv_eig<neigs && neigs != Inf
        # Sort the errors
        idx=sortperm(err[end,1:k]); err=err[end,idx];
        kk=length(λ);
        λ=λ[idx[1:kk]]; W=W[:,idx[1:kk]]; err=err[1:kk];
        msg="Number of iterations exceeded. maxit=$(maxit)."
        if conv_eig<3
            msg=string(msg, "Try to change the inner_solver_method for better performance.")
        end
        throw(NoConvergenceException(λ,Q,err,msg))
    end

    return λ,W,err,V[:,1:k+1], H[1:k,1:k-1], ω[1:k], HH[1:k,1:k]
end

# this function computes V *= h avoiding allocations (overwrites V)
function scal_mul!(V,h)
    n,m=size(V)
    for j=1:m
        for i=1:n
            @inbounds V[i,j] = h*V[i,j]
        end
    end
end

# this function computes V.-= h*W avoiding allocations (overwrites V)
function mul_and_sub!(V,W,h)
    n,m=size(V)
    for j=1:m
        for i=1:n
            @inbounds V[i,j] -= h*W[i,j]
        end
    end
end

# this function set to zero the variable V avoiding allocations
function mat_zero!(V)
    n,m=size(V)
    for j=1:m
        for i=1:n
            @inbounds V[i,j] = 0
        end
    end
end

# this function computes sum(V.*W) avoiding allocations
function mat_sum(V,W)
    TT=eltype(V)
    β::TT=0
    n,m=size(V)
    for j=1:m
        for i=1:n
            @inbounds β += V[i,j]*W[i,j]
        end
    end
    return β
end

# Contructors for the precomputed data
function PrecomputeDataInit(::Type{compute_Bmul_method_DEP})
    return IlanPrecomputeDataDEP(0,0,0,0,0)
end
function PrecomputeDataInit(::Type{compute_Bmul_method_DerSPMF})
    return IlanPrecomputeDataSPMF_NEP(0,0,0,0)
end
function PrecomputeDataInit(::Type{compute_Bmul_method_SPMF_NEP})
    return IlanPrecomputeDataSPMF_NEP(0,0,0,0)
end

# functions that precompute datas
function precompute_data(T,nep::NEPTypes.DEP,::Type{compute_Bmul_method_DEP},n,m,σ,γ)
    τ=nep.tauv
    precomp=PrecomputeDataInit(compute_Bmul_method_DEP)
    precomp.ZZ=zeros(T,n,m+1)       # aux matrix for pre-allocation
    precomp.QQ=zeros(T,n,m+1)       # aux matrix for pre-allocation
    precomp.QQ2=zeros(T,n,m+1)      # aux matrix for pre-allocation
    precomp.vv=zeros(T,m+1,length(τ))
    for j=1:length(τ)
        precomp.vv[:,j]=sqrt(τ[j]*γ)*exp(-σ)*(-τ[j]*γ).^(0:m)
    end
    precomp.G=symmetrizer_coefficients(T,m)

    return precomp
end

# functions that precompute datas
function precompute_data(T,nep::NEPTypes.SPMF_NEP,::Type{compute_Bmul_method_SPMF_NEP},n,m,σ,γ)
    precomp=PrecomputeDataInit(compute_Bmul_method_SPMF_NEP)
    precomp.ZZ=zeros(T,n,m+1)       # aux matrix for pre-allocation
    precomp.QQ=zeros(T,n,m+1)       # aux matrix for pre-allocation

    # getting matrices and functions
    fv=get_fv(nep); p=length(fv); Av=get_Av(nep)

    # precompute derivatives and FDH matrices
    SS=diagm(0 => σ*ones(T,2m+2)) + diagm(-1 => γ*(1:2m+1))
    fD=zeros(T,2*m+2,p)
    for t=1:p fD[:,t]=fv[t](SS)[:,1] end
    precomp.FDH=Vector{Array{T,2}}(undef,p)
    for t=1:p precomp.FDH[t]=zeros(T,m+1,m+1)
        for i=1:m+1 for j=1:m+1
            precomp.FDH[t][i,j]=fD[i+j,t];
        end end
    end
    precomp.G=symmetrizer_coefficients(T,m)
    return precomp
end

# functions that precompute datas
function precompute_data(T,nep::NEPTypes.DerSPMF,::Type{compute_Bmul_method_DerSPMF},n,m,σ,γ)
    precomp=PrecomputeDataInit(compute_Bmul_method_DerSPMF)
    precomp.ZZ=zeros(T,n,m+1)       # aux matrix for pre-allocation
    precomp.QQ=zeros(T,n,m+1)       # aux matrix for pre-allocation

    # getting matrices and functions
    fv=get_fv(nep); p=length(fv); Av=get_Av(nep)

    precomp.FDH=Vector{Array{T,2}}(undef,p)
    for t=1:p precomp.FDH[t]=zeros(T,m+1,m+1)
        for i=1:m+1 for j=1:m+1
            precomp.FDH[t][i,j]=nep.fD[i+j,t];
        end end
    end
    precomp.G=symmetrizer_coefficients(T,m)
    return precomp
end

function Bmult!(::Type{compute_Bmul_method_SPMF_NEP},k,Z,Qn,Av,precomp,γ)
    # B-multiplication
    #Z[:,:]=zero(Z);
    mat_zero!(Z)
    @inbounds for t=1:length(Av)
        mul!(view(precomp.QQ,:,1:k+1),view(Qn,:,1:k+1),view(precomp.G,1:k+1,1:k+1).*view(precomp.FDH[t],1:k+1,1:k+1))
        mul!(view(precomp.ZZ,:,1:k+1),Av[t],view(precomp.QQ,:,1:k+1))
        Z .+= view(precomp.ZZ,:,1:k+1);
    end
end

function Bmult!(::Type{compute_Bmul_method_DerSPMF},k,Z,Qn,Av,precomp,γ)
    Bmult!(compute_Bmul_method_SPMF_NEP,k,Z,Qn,Av,precomp,γ)
end

function Bmult!(::Type{compute_Bmul_method_DEP},k,Z,Qn,Av,precomp,γ)
    # B-multiplication
    T=eltype(Z)
    n=size(Z,1)
#    Z[:,:]=zeros(T,n,k+1)
    mat_zero!(Z)
    # low-rank factorization of G
    tolG=1e-12; U,S,V=svd(precomp.G[1:k+1,1:k+1]);q=sum(S.>tolG*ones(length(S)))
    U=view(U,:,1:q).*sqrt.(S[1:q]')
    V=view(V,:,1:q).*sqrt.(S[1:q]');
    Z[:,1]=-γ*Qn[:,1]

    @inbounds for t=2:length(Av)
        mul!(view(precomp.QQ,:,1:q),view(Qn,:,1:k+1),U.*(view(precomp.vv,1:k+1,t-1)))
        # decide multiplication order based on matrix size
        if q>k+1
            mul!(view(precomp.QQ2,:,1:k+1),view(precomp.QQ,:,1:q),(V.*view(precomp.vv,1:k+1,t-1))')
            mul!(view(precomp.ZZ,:,1:k+1),Av[t],view(precomp.QQ2,:,1:k+1));
        else
            mul!(view(precomp.QQ2,:,1:q),Av[t],view(precomp.QQ,:,1:q));
            mul!(view(precomp.ZZ,:,1:k+1),view(precomp.QQ2,:,1:q),(V.*view(precomp.vv,1:k+1,t-1))')
        end
        Z .-= view(precomp.ZZ,:,1:k+1)
    end
end

function symmetrizer_coefficients(T,m)
    G=zeros(T,m+1,m+1); # symmetrizer coefficients
    for i=1:m+1 G[i,1]=1/i end
    for j=1:m for i=1:m+1
            G[i,j+1]=(G[i,j]*j)/(i+j);
    end end
    return G
end


text:
using DM, Test

@testset "Parameter Set" begin
    d = Dict{String,String}("alpha" => "%.2f", "beta" => "%d")
    v = Dict("alpha" => pi, "beta" => 1)
    activate_param_set(d)
    @test get_param_set() == d
    @test "alpha=3.14_beta=1" == DM.print_params(d, v, ["alpha", "beta"])
    activate_param_set("alpha" => "%.2f", "beta" => "%d")
    @test get_param_set() == d
    activate_param_set("task/params.json")
    exp = Dict{String,String}(
                "alex" => "%.2f",
                "bob" => "%.2f",
                "eve" => "%.2f",
                "sandy" => "%.2f",
                "anna" => "%.2f",
        )
    @test get_param_set() == exp

    activate_param_set("task/params.json")
    d, v = load_config_from_json("task/test_1.json")
    d1 = d["task1"]
    d2 = d["task2"]
    @test DM.get_folder_path(d1, v) ==
              joinpath("./data", "alex=1.20_bob=4.00/eve=25.13_sandy=0.21")
    @test DM.get_group_path(d1, v) == "anna=0.45"
end

@testset "HDF5 file support" begin
    activate_param_set("task/params.json")
    d, v = load_config_from_json("task/test_1.json")
    d1 = d["task1"]
    d2 = d["task2"]
    test_data = collect(range(0, stop=1, length=100))
    @info "Saving test hdf5 file."
    save(d1, v, "data.h5", "x", test_data)
    save(d1, v, "data.h5", "y", [1, 2, 3], "z", "hello")
    testattrs = Dict("test_attr_1"=>"test_val_1", "test_attr_2"=>2)
    testattrs_x = Dict("test_attr_1" => 1.2, "test_attr_2"=> "test_val_2")
    writeattr(d1, v, "data.h5", testattrs)
    writeattr(d1, v, "data.h5", "x", testattrs_x)
    @test isfile("./data/alex=1.20_bob=4.00/eve=25.13_sandy=0.21/data.h5")
    @test check(d1, v, "data.h5", "x") == true
    @test load(d1, v, "data.h5", "x") == test_data
    test_x = load(d1, v, "data.h5")
    @test read(test_x["anna=0.45/x"]) == test_data
    close(test_x)
    loaded_attrs = readattr(d1, v, "data.h5")
    loaded_attrs_x = readattr(d1, v, "data.h5", "x")
    @test loaded_attrs["test_attr_1"] == "test_val_1"
    @test loaded_attrs_x["test_attr_1"] ≈ 1.2
    @test readattr(d1, v, "data.h5", ["test_attr_2"])["test_attr_2"] == 2
    @test readattr(d1, v, "data.h5", "x", ["test_attr_2"])["test_attr_2"] == "test_val_2"
    save(d2, v, "data.h5", "x", test_data)
    @test isfile("./data/alex=1.20_bob=4.00_eve=25.13_sandy=0.21/data.h5")
    delete(d1, v, "data.h5", "y", "z")
    @test !check(d1, v, "data.h5", "y")
    @info "Deleting saved hdf5 file."
    Sys.iswindows() && GC.gc()
    delete(d1, v, "data.h5")
    delete(d2, v, "data.h5")

    save(d1, v, "data-1.h5", "x", test_data)
    save(d1, v, "data-2.h5", "x", test_data)
    load_file_array(d1, v, "data", "x") == [test_data, test_data]
    Sys.iswindows() && GC.gc()
    delete(d1, v, "data-1.h5")
    delete(d1, v, "data-2.h5")

    test_entry = DataEntry("./data", [[]], [[]])
    save(test_entry, Dict(), "test.h5", "x", test_data)
    @test isfile("./data/test.h5")
    @test load(test_entry, Dict(), "test.h5", "x") == test_data
    Sys.iswindows() && GC.gc()
    delete(test_entry, Dict(), "test.h5")
end

@testset "CSV file support" begin
    activate_param_set("task/params.json")
    d, v = load_config_from_json("task/test_1.json")
    df = DataFrame(A=1:4, B=["M", "F", "F", "M"])
    @info "Saving test CSV file."
    save(d["task1"], v, "test.csv", df)
    @test isfile("./data/alex=1.20_bob=4.00/eve=25.13_sandy=0.21/test.csv")
    @test load(d["task1"], v, "test.csv") == df
    Sys.iswindows() && GC.gc()
    @info "Deleting saved CSV file."
    delete(d["task1"], v, "test.csv")
end

@testset "Registry" begin
    activate_param_set("task/params.json")
    d, v = load_config_from_json("task/test_1.json")
    d1 = d["task1"]
    @test_throws ArgumentError load_registry("test_1.csv")
    reg = load_registry("test_1.csv", d1)
    register_entry(reg, d1, v, "data.h5", "res")
    @test reg == load_registry("data/test.csv")
    write_registry(reg, "data/test_1.csv")
    @test isfile("data/test_1.csv")

    v["sandy"] = 0.3
    register_entry(reg, d1, v, "data.h5", "res")
    v1 = delete!(copy(v), "sandy")
    @test query_registry(reg, v) == query_registry(reg, v1, (x)->parse(Float64, x.sandy) > 0.25)

    Sys.iswindows() && GC.gc()
    rm("data/test_1.csv")
end

@testset "Utility function" begin
    test_array = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
    @test convert_recursive_array(test_array) == [[1, 2, 3] [4, 5, 6] [7, 8, 9]]
    test_array = [[[1,2],[3,4]],[[5,6],[7,8]],[[9,10],[10,11]]]
    @test convert_recursive_array(test_array) == cat([hcat(x...) for x in test_array]..., dims=3)
    test_dict = parse_string_to_num(Dict("hello"=>"2*pi"))
    @test test_dict["hello"] ≈ 2π
end


text:
# de Kort, Selvino R. and Correia, Sérgio P. C. and Alexis, Dean M. and
# Dickinson, Anthony and Clayton, Nicola S., The control of food-caching
# behavior by Western scrub-jays (Aphelocoma californica). Journal of
# Experimental Psychology: Animal Behavior Processes 33:361-370, 2007.
# http://dx.doi.org/10.1037/0097-7403.33.4.361
#

function _summarize(::Experiment{:deKort07_exp1}, results)
	combine(groupby(results, [:exp, :group, :action, :trial]),
				 df -> DataFrame(μ = mean(df.cache), sem = sem(df.cache)))
end
function statistical_tests(exp::Experiment{:deKort07_exp1}, data)
    tests = Test[]
    for g in groupby(data, :exp)
        id = "Experiment $(g.exp[1])"
        test = Test(id,
                    splitplotanova(g,
                                   exp.tests[id].locals,
                                   betweenfactors = :group,
                                   withinfactors = :trial,
                                   id = :id, data = :cache),
                    exp.tests)
        push!(tests, test)
    end
    TestCollection(tests)
end
function run!(exp::Experiment{:deKort07_exp1}, models; N = 16, samemodels = false)
    results = DataFrame(exp = String[], group = String[], action = String[],
                        trial = Int[],
						id = Int[], cache = Int[])
	id = 0
	for exp in ("a", "b")		# experiments
		id = samemodels ? 0 : id
		for i in 1:N			# birds
			id += 1
			m = models[id]
			add!(m, MaintenanceDiet)
			for trial in 1:6	# trials
				remove!(m, MaintenanceDiet)
				wait!(m, 3u"hr")
				add!(m, Mealworm, 50)
				trayx = Tray(0, trial)
				add!(m, trayx)
				wait!(m, 15u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				cache = countcache(trayx)
				wait!(m, 28.25u"hr")
				remove!(m, MaintenanceDiet)
				wait!(m, 1u"hr")
				group = "replenish"
				if exp == "a" && id > div(N, 2)
					degrade!(trayx)
					group = "degrade"
				elseif exp == "b" && id % 2 == 1
					pilfer!(trayx)
					group = "pilfer"
				end
				push!(results, [exp, group, "cache", trial, id, cache])
				add!(m, trayx)
				wait!(m, 10u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				wait!(m, 20u"minute" + 15u"hr" + 1u"d")
			end
		end
	end
	results
end

function _summarize(::Experiment{:deKort07_exp2}, results)
	combine(groupby(results, [:exp, :tray, :action, :trial]),
				 df -> DataFrame(μ = mean(df.cache), sem = sem(df.cache)))
end
function statistical_tests(exp::Experiment{:deKort07_exp2}, data)
    groups = groupby(data, [:exp])
    tests = Test[]
    for g in groups
        id = "Experiment $(g.exp[1]) tray × trial"
        test = Test(id,
                    splitplotanova(g,
                                   exp.tests[id].locals,
                                   betweenfactors = [],
                                   withinfactors = [:tray, :trial],
                                   id = :id, data = :cache),
                    exp.tests)
        push!(tests, test)
        subgroups = groupby(g, :tray)
        for (i, sb) in enumerate(subgroups)
            id = "Experiment $(g.exp[1]), tray $(sb.tray[1])"
            test = Test(id,
                        anova(sb, @formula(cache ~ trial + id), exp.tests[id].locals),
                        exp.tests)
            push!(tests, test)
        end
    end
    push!(tests, Test("Experiment a, trial 6",
                      is_significantly_bigger(@where(groups[1], :trial .== 6),
                                         condbigger = "replenish",
                                         condsmaller = "degrade",
                                         conditionkey = :tray,
                                         valuekey = :cache),
                      exp.tests))
    push!(tests, Test("Experiment b, trial 4",
                      is_significantly_bigger(@where(groups[2], :trial .== 4),
                                         condbigger = "replenish",
                                         condsmaller = "pilfer",
                                         conditionkey = :tray,
                                         valuekey = :cache),
                      exp.tests))
    TestCollection(tests)
end
function run!(exp::Experiment{:deKort07_exp2}, models; Na = 10, Nb = 8, N = 0, samemodels = false)
    results = DataFrame(exp = String[], tray = String[], action = String[],
                        trial = Int[],
						id = Int[], cache = Int[])
	id = 0
	for exp in ("a", "b")		# experiments
		id = samemodels ? 0 : id
		for i in 1:(N != 0 ? N : exp == "a" ? Na : Nb)	# birds
			id += 1
			m = models[id]
			add!(m, MaintenanceDiet)
			for trial in 1:(exp == "a" ? 6 : 4)	# trials
				remove!(m, MaintenanceDiet)
				wait!(m, 3u"hr")
				add!(m, Mealworm, 25)
				trayx = Tray(0, trial)
				trayy = Tray(1, trial + 10)
				add!(m, trayx)
				add!(m, trayy)
				wait!(m, 15u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				cachex = countcache(trayx)
				cachey = countcache(trayy)
				wait!(m, 28.25u"hr")
				remove!(m, MaintenanceDiet)
				wait!(m, 1u"hr")
				if exp == "a"
					degrade!(trayx)
					tray = "degrade"
				else
					pilfer!(trayx)
					tray = "pilfer"
				end
				push!(results, [exp, tray, "cache", trial, id, cachex])
				push!(results, [exp, "replenish", "cache", trial, id, cachey])
				add!(m, trayx)
				add!(m, trayy)
				wait!(m, 10u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				wait!(m, 20u"minute" + 15u"hr" + 1u"d")
			end
		end
	end
	results
end

function _summarize(::Experiment{:deKort07_exp3}, results)
	combine(groupby(results, [:tray, :trial, :action]),
				 df -> DataFrame(μ = mean(df.cache), sem = sem(df.cache)))
end
function statistical_tests(exp::Experiment{:deKort07_exp3}, data)
    tests = Test[]
    id = "tray × trial"
    push!(tests, Test(id,
                      splitplotanova(data,
                                     exp.tests[id].locals,
                                     betweenfactors = [:tray],
                                     withinfactors = [:trial],
                                     id = :id, data = :cache),
                      exp.tests))
    groups = groupby(data, :tray)
    for g in groups
        id = "Tray $(g.tray[1])"
        push!(tests, Test(id,
                          anova(g, @formula(cache ~ trial + id), exp.tests[id].locals),
                          exp.tests))
    end
    TestCollection(tests)
end
function run!(exp::Experiment{:deKort07_exp3}, models; N = 7)
    results = DataFrame(tray = String[], trial = Int[], action = String[],
						id = Int[], cache = Int[])
	for id in 1:N			# birds
		m = models[id]
		add!(m, MaintenanceDiet)
		for trial in 1:10	# trials
			remove!(m, MaintenanceDiet)
			wait!(m, 3u"hr")
			add!(m, Waxworm, 30)
			trayx = Tray(0, trial)
			trayy = Tray(1, trial + 10)
			add!(m, trayx)
			add!(m, trayy)
            cover!(trayy)
			wait!(m, 15u"minute")
			add!(m, Waxworm, 10)
			uncover!(trayy)
			cover!(trayx)
			wait!(m, 15u"minute")
			remove!(m, Any)
			add!(m, MaintenanceDiet)
			uncover!(trayx)
			cachex = countcache(trayx)
			cachey = countcache(trayy)
			wait!(m, 28u"hr")
			remove!(m, MaintenanceDiet)
			wait!(m, 1u"hr")
			pilfer!(trayx)
			push!(results, ["period1", trial, "cache", id, cachex])
			push!(results, ["period2", trial, "cache", id, cachey])
			add!(m, trayx)
			add!(m, trayy)
			wait!(m, 10u"minute")
			remove!(m, Any)
			add!(m, MaintenanceDiet)
			wait!(m, 20u"minute" + 15u"hr" + 1u"d")
		end
	end
	results
end

function _summarize(::Experiment{:deKort07_exp4}, results)
	combine(groupby(results, [:exp, :group, :action, :tray]),
				 df -> DataFrame(μ = mean(df.cache), sem = sem(df.cache)))
end
function statistical_tests(exp::Experiment{:deKort07_exp4}, data)
    groups = groupby(data, :exp)
    tests = Test[]
    id = "group × tray"
    test = Test(id,
                splitplotanova(groups[1],
                               exp.tests[id].locals,
                               betweenfactors = :group, withinfactors = :tray,
                               id = :id, data = :cache),
                exp.tests)
    push!(tests, test)
    push!(tests, Test("experimental group caches more in B than control",
                      is_significantly_bigger(@where(groups[1], :tray .== "B"),
                                              conditionkey = :group,
                                              condbigger = "experimental",
                                              condsmaller = "control",
                                              valuekey = :cache, paired = false),
                      exp.tests))
    push!(tests, Test("number of worms cache in A is the same across groups",
                      isnot_significantly_different(@where(groups[1], :tray .== "A"),
                                                    conditionkey = :group,
                                                    condition1 = "experimental",
                                                    condition2 = "control",
                                                    valuekey = :cache, paired = false),
                      exp.tests))
    push!(tests, Test("experimental group cache same numbers in A and B",
                      isnot_significantly_different(@where(groups[1], :group .== "experimental"),
                                                    conditionkey = :tray,
                                                    condition1 = "A",
                                                    condition2 = "B",
                                                    valuekey = :cache,
                                                    paired = false),
                      exp.tests))
    push!(tests, Test("control group caches more in A than in B",
                      is_significantly_bigger(@where(groups[1], :group .== "control"),
                                                     conditionkey = :tray,
                                                     condbigger = "A",
                                                     condsmaller = "B",
                                                     valuekey = :cache,
                                                     paired = false),
                      exp.tests))
    id = "Effect of tray on the worms cache"
    push!(tests, Test(id,
                      anova(groups[2], @formula(cache ~ tray + id), exp.tests[id].locals),
                      exp.tests))
    push!(tests, Test("more in B than A",
                      is_significantly_bigger(groups[2],
                                              conditionkey = :tray,
                                              condbigger = "B",
                                              condsmaller = "A",
                                              valuekey = :cache, paired = true),
                      exp.tests))
    push!(tests, Test("more in C than A",
                      is_significantly_bigger(groups[2],
                                              conditionkey = :tray,
                                              condbigger = "C",
                                              condsmaller = "A",
                                              valuekey = :cache, paired = true),
                      exp.tests))
    push!(tests, Test("same in B and C",
                      isnot_significantly_different(groups[2],
                                                    conditionkey = :tray,
                                                    condition1 = "B",
                                                    condition2 = "C",
                                                    valuekey = :cache, paired = true),
                      exp.tests))
    TestCollection(tests)
end
function run!(exp::Experiment{:deKort07_exp4}, models; Na = 8, Nb = 4, N = 0, samemodels = false)
	results = DataFrame(exp = String[], group = String[], tray = String[],
                        action = String[],
						id = Int[], cache = Int[])
	id = 0
	for exp in ("a", "b")						# experiments
		id = samemodels ? 0 : id
		for i in 1:(N != 0 ? N : exp == "a" ? Na : Nb)		# birds
			id += 1
			m = models[id]
			add!(m, MaintenanceDiet)
            for trial in 3:5	# 2 pretraining (not modeled), 2 training, 1 testing
				remove!(m, MaintenanceDiet)
				wait!(m, 3u"hr")
				add!(m, Waxworm, 30)
				trayx = Tray(0, trial)
				add!(m, trayx)
				if trial > 2
					trayy = Tray(1, trial + 10)
					if trial < 5; cover!(trayy); end
					add!(m, trayy)
                    if exp == "b"
                        trayz = Tray(2, trial + 20)
                        if trial < 5; cover!(trayz); end
                        add!(m, trayz)
                    end
				end
				wait!(m, 15u"minute")
				if trial == 5
					group = id > div(Na, 2) && exp == "a" ? "control" : "experimental"
					push!(results, [exp, group, "A", "cache", id, countcache(trayx)])
					push!(results, [exp, group, "B", "cache", id, countcache(trayy)])
					if exp == "b"
						push!(results, [exp, group, "C", "cache", id, countcache(trayz)])
					end
                    remove!(m, Any)
                    break # end of experiment
				end
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				wait!(m, 28.25u"hr")
				remove!(m, MaintenanceDiet)
				wait!(m, 1u"hr")
				if trial > 2
					uncover!(trayy)
					add!(m, trayy)
					if exp == "b"
						uncover!(trayz)
						add!(m, trayz)
					end
					if id <= div(Na, 2) || exp == "b"
						trayy.eatableitems = deepcopy(trayx.eatableitems)
						empty!(trayx.eatableitems)
					end
				end
				add!(m, trayx)
				wait!(m, 10u"minute")
				remove!(m, Any)
				add!(m, MaintenanceDiet)
				wait!(m, 20u"minute" + 15u"hr" + 1u"d")
			end
			wait!(m, 11u"d")
		end
	end
	results
end


text:
export TagSet,
       addtags,
       hastags,
       Tag

const Tag = SmallString
const maxTagLength = smallLength
const MTagStorage = MSmallStringStorage # A mutable tag storage
const IntTag = IntSmallString  # An integer that can be cast to a Tag
const maxTags = 4
const TagSetStorage = SVector{maxTags,IntTag}
const MTagSetStorage = MVector{maxTags,IntTag}  # A mutable tag storage

struct TagSet
  tags::TagSetStorage
  plev::Int
  length::Int
  function TagSet()
    ts = TagSetStorage(ntuple(_ -> IntTag(0),Val(4)))
    plev = -1
    new(ts,plev,0)
  end
  TagSet(tags::TagSetStorage,plev::Int,len::Int) = new(tags,plev,len)
end

#function TagSet(tags::TagSetStorage,plev::Int=-1)
#  len = 0
#  while tags[len+1] ≠ IntTag(0)
#    len += 1
#  end
#  TagSet(tags,plev,len)
#end

TagSet(ts::TagSet) = ts

function _hastag(ts::MTagSetStorage, ntags::Int, tag::IntTag)
  for n = 1:ntags
    @inbounds ts[n] == tag && return true
  end
  return false
end

function _addtag_ordered!(ts::MTagSetStorage, ntags::Int, tag::IntTag)
  if ntags == 0 || tag > @inbounds ts[ntags]
    @inbounds setindex!(ts,tag,ntags+1)
  else
    # check for repeated tags
    _hastag(ts,ntags,tag) && return ntags
    pos = ntags+1   # position new tag should go
    while pos > 1 && tag < @inbounds ts[pos-1]
      pos -= 1
      @inbounds setindex!(ts,ts[pos],pos+1)
    end
    @inbounds setindex!(ts,tag,pos)
  end
  return ntags+1
end

function _addtag!(ts::MTagSetStorage, plev::Int, ntags::Int, tag::IntTag)
  plnew = plev
  t = Tag(tag)
  if !isnull(t)
    if isint(t)
      error("""Cannot use a bare integer as a tag.
            If you are looking to set the prime level, use the syntax ("x",1)
            (instead of "x,1")""")
      #plev ≥ 0 && error("You can only make a TagSet with one prime level/integer tag.")
      #plnew = parse(Int,t)
    else
      ntags = _addtag_ordered!(ts, ntags,tag)
    end
  end
  return plnew, ntags
end

#isnull(v::MTagStorage) = v[0] == IntChar(0)

function reset!(v::MTagStorage, nchar::Int)
  for i = 1:nchar
    @inbounds v[i] = IntChar(0)
  end
end

function TagSet(str::AbstractString)
  # Mutable fixed-size vector as temporary Tag storage
  current_tag = MTagStorage(ntuple(_ -> IntChar(0),Val(maxTagLength)))
  # Mutable fixed-size vector as temporary TagSet storage
  ts = MTagSetStorage(ntuple(_ -> IntTag(0),Val(maxTags)))
  nchar = 0
  ntags = 0
  plev = -1
  for n = 1:length(str)
    @inbounds current_char = str[n]
    if current_char == ','
      if nchar != 0
        plev, ntags = _addtag!(ts,plev,ntags,cast_to_uint64(current_tag))
        # Reset the current tag
        reset!(current_tag,nchar)
        nchar = 0
      end
    elseif current_char != ' ' # TagSet constructor ignores whitespace
      nchar == maxTagLength && error("Currently, tags can only have up to $maxTagLength characters")
      nchar += 1
      @inbounds current_tag[nchar] = current_char
    end
  end
  # Store the final tag
  if nchar != 0
    plev, ntags = _addtag!(ts,plev,ntags,cast_to_uint64(current_tag))
  end
  return TagSet(TagSetStorage(ts),plev,ntags)
end

function TagSet((str,plev)::Tuple{<:AbstractString,<:Integer})
  ts = TagSet(str)
  return setprime(ts,plev)
end

convert(::Type{TagSet}, str::String) = TagSet(str)

tags(T::TagSet) = T.tags
plev(T::TagSet) = T.plev
Base.length(T::TagSet) = T.length
Base.getindex(T::TagSet,n::Int) = Tag(getindex(tags(T),n))
#Base.setindex!(T::TagSet,val,n::Int) = TagSet(setindex(tags(T),val,n),plev(T),length(T))
Base.copy(ts::TagSet) = TagSet(tags(ts),plev(ts),length(ts))

# A TagSet is considered to have a prime level if
# the prime level is greater than or equal to zero
hasplev(ts::TagSet) = (plev(ts) ≥ 0)

setprime(ts::TagSet,pl::Int) = TagSet(tags(ts),pl,length(ts))
function prime(T::TagSet,plinc::Int=1)
  if !hasplev(T)
    return setprime(T,plinc)
  else
    return setprime(T,plev(T)+plinc)
  end
end

# TODO: define iteration?
#Base.iterate(ts::TagSet,state::Int=1) = iterate(tags(ts),state)

# TODO: define `in` in terms of hastags?
#Base.in(s, ts::TagSet) = hastags(ts, TagSet(s))

# Cast SVector of IntTag of length 4 to SVector of UInt128 of length 2
# This is to make TagSet comparison a little bit faster
function cast_to_uint128(a::TagSetStorage)
  return unsafe_load(convert(Ptr{SVector{2,UInt128}},pointer_from_objref(MTagSetStorage(a))))
end

import Base.==
function ==(ts1::TagSet,ts2::TagSet)
  plev(ts1) ≠ plev(ts2) && return false
  # Block the bits together to make the comparison faster
  return cast_to_uint128(tags(ts1)) == cast_to_uint128(tags(ts2))
end

function hastag(ts::TagSet, tag)
  l = length(ts)
  l < 1 && return false
  for n = 1:l
    @inbounds Tag(tag) == ts[n] && return true
  end
  return false
end

function hastags(ts2::TagSet, tags1)
  ts1 = TagSet(tags1)
  (hasplev(ts1) && plev(ts1) ≠ plev(ts2)) && return false
  l1 = length(ts1)
  l2 = length(ts2)
  l1 > l2 && return false
  for n1 = 1:l1
    @inbounds !hastag(ts2,ts1[n1]) && return false
  end
  return true
end

function addtags(ts::TagSet, tagsadd)
  if length(ts) == maxTags
    throw(ErrorException("Cannot add tag: TagSet already maximum size"))
  end
  tsadd = TagSet(tagsadd)
  (hasplev(ts) && hasplev(tsadd)) && error("In addtags(::TagSet,...), cannot add a prime level")
  res_ts = MVector(tags(ts))
  res_plev = plev(ts)
  ntags = length(ts)
  for n = 1:length(tsadd)
    @inbounds ntags = _addtag_ordered!(res_ts, ntags,IntSmallString(tsadd[n]))
  end
  if !hasplev(ts) && hasplev(tsadd)
    res_plev = plev(tsadd)
  end
  return TagSet(TagSetStorage(res_ts),res_plev,ntags)
end

function _removetag!(ts::MTagSetStorage, ntags::Int, t::Tag)
  for n = 1:ntags
    if @inbounds Tag(ts[n]) == t
      for j = n:ntags-1
        @inbounds ts[j] = ts[j+1]
      end
      @inbounds ts[ntags] = IntTag(0)
      return ntags -= 1
    end
  end
  return ntags
end

#TODO: optimize this function
function removetags(ts::TagSet, tagsremove)
  tsremove = TagSet(tagsremove)
  hasplev(tsremove) && error("In removetags(::TagSet,...), cannot remove a prime level")
  res_ts = MVector(tags(ts))
  res_plev = plev(ts)
  ntags = length(ts)
  for n=1:length(tsremove)
    @inbounds ntags = _removetag!(res_ts, ntags, tsremove[n])
  end
  return TagSet(TagSetStorage(res_ts),res_plev,ntags)
end

#TODO: optimize this function
function replacetags(ts::TagSet, tagsremove, tagsadd)
  tsremove = TagSet(tagsremove)
  tsadd = TagSet(tagsadd)
  if hasplev(tsremove) != hasplev(tsadd)
    error("In replacetags(::TagSet,...), cannot remove or add a prime level")
  end
  res_ts = MVector(tags(ts))
  res_plev = plev(ts)
  ntags = length(ts)
  # The TagSet must have the same prime level as the one being replaced
  (hasplev(tsremove) && (res_plev≠plev(tsremove))) && return ts
  # The TagSet must have the tags to be replaced
  !hastags(ts,tsremove) && return ts
  for n = 1:length(tsremove)
    @inbounds ntags = _removetag!(res_ts, ntags, tsremove[n])
  end
  for n = 1:length(tsadd)
    @inbounds ntags = _addtag_ordered!(res_ts, ntags,IntSmallString(tsadd[n]))
  end
  if hasplev(ts) && (plev(ts)==plev(tsremove))
    res_plev = plev(tsadd)
  end
  return TagSet(TagSetStorage(res_ts),res_plev,ntags)
end

function primestring(ts::TagSet)
  if !hasplev(ts)
    return " (warning: no prime level)"
  end
  pl = plev(ts)
  if pl==0
    return ""
  elseif pl > 3
    return "'$pl"
  else
    return "'"^pl
  end
end

function tagstring(T::TagSet)
  res = ""
  length(T) == 0 && return res
  for n=1:length(T)-1
    res *= "$(Tag(T[n])),"
  end
  res *= "$(Tag(T[length(T)]))"
  return res
end

function show(io::IO, T::TagSet)
  print(io,"(")
  lT = length(T)
  if lT > 0
    print(io,T[1])
    for n=2:lT
      print(io,",$(T[n])")
    end
  end
  print(io,")")
  print(io,primestring(T))
end

function readcpp(io::IO,::Type{TagSet}; kwargs...)
  format = get(kwargs,:format,"v3")
  ts = TagSet()
  if format=="v3"
    mstore = MTagSetStorage(ntuple(_ -> IntTag(0),Val(maxTags)))
    ntags = 0
    for n=1:4
      t = readcpp(io,Tag;kwargs...)
      if t != Tag()
        ntags = _addtag_ordered!(mstore,ntags,IntSmallString(t))
      end
    end
    plev = convert(Int,read(io,Int32))
    ts = TagSet(TagSetStorage(mstore),plev,ntags)
  else
    throw(ArgumentError("read TagSet: format=$format not supported"))
  end
  return ts
end


text:
baremodule SoundStatus
    const STOPPED = 0
    const PAUSED = 1
    const PLAYING = 2
end


text:
using RoME, TransformUtils, KernelDensityEstimate
using IncrementalInference, CloudGraphs, Neo4j

include("BlandAuthDB.jl")

configuration = CloudGraphs.CloudGraphConfiguration(dbaddress, 7474, dbusr, dbpwd, mongoaddress, 27017, false, "", "");
cloudGraph = connect(configuration);
# register types of interest in CloudGraphs
registerGeneralVariableTypes!(cloudGraph)
IncrementalInference.setCloudDataLayerAPI!()

N = 200
fg = initfg()
fg.sessionname="SESSANIM"
fg.cg = cloudGraph

initCov = diagm([0.03;0.03;0.001])
odoCov = diagm([3.0;3.0;0.4])

x1 = initFactorGraph!(fg, ready=0, N=N)
println("Starting $(fg.sessionname) at $(x1), at $(getKDEMax(getVertKDE(fg,x1)))")

r = 3.0/12.0*pi
cDX = [10.0;0;r]
for i in 2:13
	addOdoFG!(fg, "x$(i)", cDX, odoCov, ready=0, N=N)
end

# insert loop closure
v1a = getVert(fg, "x1")
v1b = getVert(fg, "x13")

DX = zeros(3,1)
pp = Pose2Pose2(DX, 0.1*eye(3), [1.0]) #[prev;v],
f = addFactor!(fg, [v1a;v1b], pp, ready=0)

v1c = getVert(fg, "x7")

DX = zeros(3,1)
DX[2] = 38.63
DX[3] = pi
pp = Pose2Pose2(DX, 0.1*eye(3), [1.0]) #[prev;v],
f = addFactor!(fg, [v1a;v1c], pp, ready=0)

v2a = getVert(fg, "x4")
v2b = getVert(fg, "x10")
pp = Pose2Pose2(DX, 0.1*eye(3), [0.1;0.9]) #[prev;v],
f = addFactor!(fg, [v2a;v2b], pp, ready=0)

setDBAllReady!(fg)

#