text:
#==========================================================================================#

# TABLE OF ESTIMATES

function etable(args...;
        digits::Int = 4,
        aux::Union{Function, Nothing} = nothing,
        stars::Matrix{Any} = [0.1 "*"; 0.05 "**"; 0.01 "***"],
        titles::Vector{String} = [""]
    )

    N                  = length(args)
    cutpoints, symbols = _parsestars(stars)

    (titles == [""]) && (titles = String["(" * string(i) * ")" for i = 1:N])

    μ = Vector{Vector{String}}(undef, N)
    β = Vector{Vector{String}}(undef, N)
    σ = Vector{Vector{String}}(undef, N)

    for (i, ai) in enumerate(args)

        μ[i] = coefnames(ai)
        β[i] = format.("{:.$(digits)f}", coef(ai))

        if typeof(aux) != Nothing
            σ[i] = "(" .* format.("{:.$(digits)f}", aux(ai)) .* ")"
        end

        if cutpoints != []
            for (j, pj) in enumerate(pval(ai))
                s        = searchsortedlast(cutpoints, pj)
                β[i][j] *= symbols[s]
            end
        end
    end

    inter = union(μ...)

    if typeof(aux) == Nothing
        output = vcat("", inter)
    else
        output          = fill("", 2 * length(inter) + 1)
        output[2:2:end] = inter
    end

    w1 = similar(output)
    w2 = similar(output[2:end])

    _alignatchar!(output)

    for (i, (ni, ti, βi)) in enumerate(zip(μ, titles, β))
        idx = findall((in)(ni), inter)
        (typeof(aux) == Nothing) || (idx .= 2 .* idx .- 1)
        w2      .= ""
        w2[idx] .= βi
        (typeof(aux) == Nothing) || (w2[(idx .+ 1)] .= σ[i])
        _alignatchar!(w2, '.')
        w1[1]      = ti
        w1[2:end] .= w2
        _alignatchar!(w1)
        output .*= "    " .* w1
    end

    for i in output
        println(i)
    end
end

#==========================================================================================#

# PARSING STAR SCHEME

function _parsestars(x::Matrix{Any})
    if size(x, 1) != 0
        s         = sort(x, dims = 1, rev = true)
        cutpoints = vcat(0.0, convert(Vector{Float64}, s[:, 1]), 1.0)
        symbols   = vcat(convert(Vector{String}, s[:, 2]), "")
    else
        cutpoints = Vector{Float64}()
        symbols   = Vector{String}()
    end
    return cutpoints, symbols
end

#==========================================================================================#

# ALIGN STRING VECTOR AT FIRST CHARACTER

function _alignatchar!(x::Vector{String})
    len    = length.(x)
    maxlen = maximum(len)
    len    = maxlen .- len
    add    = " " .^ len
    x     .= x .* add
end

# ALIGN STRING VECTOR AT FIRST OCCURRENCE OF CHARACTER

function _alignatchar!(x::Vector{String}, y::Char)

    pos = fill(0, length(x))

    for (xi, xx) in enumerate(x)
        (xx == "") || (pos[xi] = findlast(z -> z == y, xx))
    end

    maxpos = maximum(pos)
    add    = " " .^ (maxpos .- pos)
    x     .= add .* x
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function trueanom(E::T, e::T) where {T<:AbstractFloat}
    @assert e >= 0 "eccentricity must be in the range [0, 1]"
    return 2 * atan(sqrt((1 + e) / (1 - e)) * tan(E / 2))
end

"""
    trueanom(E, e) -> true anomaly

### Purpose ###

Calculate true anomaly for a particle in elliptic orbit with eccentric anomaly
\$E\$ and eccentricity \$e\$.

### Explanation ###

In the two-body problem, once that the [Kepler's
equation](https://en.wikipedia.org/wiki/Kepler%27s_equation) is solved and
\$E(t)\$ is determined, the polar coordinates \$(r(t), \\theta(t))\$ of the body
at time \$t\$ in the elliptic orbit are given by

`` \\theta(t) = 2\\arctan \\left(\\sqrt{\\frac{1 + e}{1 - e}} \\tan\\frac{E(t)}{2} \\right)``

`` r(t) = \\frac{a(1 - e^{2})}{1 + e\\cos(\\theta(t) - \\theta_{0})}``

in which \$a\$ is the semi-major axis of the orbit, and \$\\theta_0\$ the value
of angular coordinate at time \$t = t_{0}\$.

### Arguments ###

* `E`: eccentric anomaly.
* `e`: eccentricity, in the elliptic motion regime (\$0 \\leq e \\leq 1\$)

### Output ###

The true anomaly.

### Example ###

Plot the true anomaly as a function of mean anomaly for eccentricity \$e = 0\$,
\$0.5\$, \$0.9\$.  Use [PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for
plotting.

```julia
using PyPlot
M = range(0, stop=2pi, length=1001)[1:end-1];
for ecc in (0, 0.5, 0.9)
    plot(M, mod2pi.(trueanom.(kepler_solver.(M, ecc), ecc)))
end
```

### Notes ###

The eccentric anomaly can be calculated with [`kepler_solver`](@ref) function.
"""
trueanom(E::Real, e::Real) = trueanom(promote(float(E), float(e))...)

ext: jl


text:
using OrbitalTrajectories

using Test
using DifferentialEquations

@testset "convert_to_frame" begin
    u0 = [0.8574053516112442, 0., 0., 0., 0.47, 0.]
    system_nbp = EphemerisNBP(:earth, :moon)
    prob = State(system_nbp, SynodicFrame(), u0, (0., 3600.0*24*30))

    # Test that converting to and from inertial frames works
    inert_prob = convert_to_frame(prob, InertialFrame())
    synod_prob = convert_to_frame(inert_prob, SynodicFrame())
    @test synod_prob.prob.u0 ≈ u0

    # Test that converting to/from non-normalised frames also works
    synod_nonorm_prob = convert_to_frame(prob, SynodicFrame(false))
    synod_norm_prob = convert_to_frame(synod_nonorm_prob, SynodicFrame(true))
    @test synod_norm_prob.prob.u0 ≈ prob.prob.u0
    inert_prob2 = convert_to_frame(synod_nonorm_prob, InertialFrame())
    synod_nonorm_prob2 = convert_to_frame(inert_prob2, SynodicFrame(false))
    @test synod_nonorm_prob2.prob.u0 ≈ synod_nonorm_prob.prob.u0

    # Roughly in decreasing order of allocations (30μs~40μs, allocations 8 to 38)
    # @btime convert_to_frame(prob, InertialFrame())  # norm -> inert
    # @btime convert_to_frame(inert_prob, SynodicFrame())  # inert -> norm
    # @btime convert_to_frame(prob, SynodicFrame(false))  # norm -> unorm
    # @btime convert_to_frame(synod_nonorm_prob, SynodicFrame(true))  # unorm -> norm
    # @btime convert_to_frame(synod_nonorm_prob, InertialFrame())  # unorm -> inert
    # @btime convert_to_frame(inert_prob2, SynodicFrame(false))  # inert -> unorm

    # Solve and test that converting the solution to synodic also matches
    sol = solve(prob)
    sol2 = convert_to_frame(sol, SynodicFrame())
    @test sol2.sol.u[begin] ≈ u0
end


text:
#=
Created on Thursday 2 April 2020
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Experiments on the host-parasite datasets.
=#

using DataFrames, CSV
using Plots, Statistics
using NamedArrays

using EcologicalNetworks

all_hp_data = filter(x -> occursin("Hadfield", x.Reference), web_of_life());

ids = getfield.(all_hp_data, :ID);
networks = web_of_life.(ids)

N = reduce(union, convert.(BinaryNetwork, networks))

parasites, hosts = species(N, dims=1), species(N, dims=2);

M = NamedArray(zeros(size(N)...), (parasites, hosts))

networks_completed = NamedArray{Float64,2}[]

for (network, id) in zip(networks, ids)
    Y = similar(M)
    fill!(Y, 0.0)

    for p in species(network, dims=1)
        for h in species(network, dims=2)
            Y[p, h] = network[p, h]
        end
    end
    push!(networks_completed, Y)
end

networks_completed

sum.(networks_completed)

using OTSIN, StatsBase

train_inds = sample(1:51, 26, replace=false)
test_inds = setdiff(1:51, train_inds);

γ = 1e-2

Ps = [i in train_inds ? Y .+ 1e-4 : Y for (i, Y) in enumerate(networks_completed)]
Ps = OTSIN.normalize.(Ps)
@time M = fitM(Matrix.(Ps[train_inds]), fix_a=true, fix_b=true, maxitter=5, γ=γ)
#@time M = fitM(Matrix.(Ps[train_inds]), fix_a=false, fix_b=false, maxitter=10, γ=1.0)

heatmap(M, color=:pu_or)
title!("Global utility matrix host-parasite")
xlabel!("hosts")
ylabel!("parasites")

savefig("figures/hostparM.svg")

kl_ot_ab = Float64[]
kl_ot_a = Float64[]
kl_ot_b = Float64[]
kl_ot_softmax = Float64[]
kl_ot_neutral = Float64[]

cutoff = 1e-6  # remove sp. with no interactions

for P in Ps[test_inds]
    P = Matrix(P)
    a, b = marginals(P)
    # remove empty rows/columns
    P = P[a.>cutoff, b.>cutoff]
    Mss = M[a.>cutoff, b.>cutoff]
    a, b = marginals(P)
    # optimal transport both
    Q = optimaltransport(Mss, a, b)
    push!(kl_ot_ab, KL(P, Q))
    Q .= optimaltransport(Mss, a, nothing)
    push!(kl_ot_a, KL(P, Q))
    Q .= optimaltransport(Mss, nothing, b)
    push!(kl_ot_b, KL(P, Q))
    Q .= optimaltransport(Mss, nothing, nothing)
    push!(kl_ot_softmax, KL(P, Q))
    Q .= optimaltransport(a, b)
    push!(kl_ot_neutral, KL(P, Q))
end

using BSON

Ps = convert.(Array, Ps)
Pstrain = Ps[train_inds]
Pstest = Ps[test_inds]

BSON.@save "results/hp.bson" Pstrain Pstest M kl_ot_ab  kl_ot_a kl_ot_b kl_ot_softmax kl_ot_neutral γ

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
function print_log(sv)
    if lg1()
        msg = string("\niter ", sv.nb_test, "=", sv.nb_reject, "+",
                    sv.nb_move)
        if lg2()
            # affiche coût + ordre des avions
            msg *= string(" => ", to_s(sv.bestsol))
        else
            # affiche seulement le coût
            msg *= string(" => ", sv.bestsol.cost)
        end
        print(msg)
    end
end


text:
doc"""
    OneStepForward!(spt2, spt1, fidMtx)

One time step forward modeling of acoustic wave field, the current snapshot is
kept in spt2. $ p^{n+1} = F \times p^{n} + s^{n+1} $

# Arguments
* `spt2 :: SnapShot`: composite type of SnapShot
* `spt1 :: SnapShot`: composite type of SnapShot
* `fidMtx :: FidMtx`: composite type of sparse partial differential matrix
"""
function OneStepForward!(spt2::SnapShot, spt1::SnapShot, fidMtx::FidMtx)
    spt2.it  = spt1.it + 1
    spt2.vz = fidMtx.MvzBvz * spt1.vz + fidMtx.MvzBp  * (spt1.pz+spt1.px)
    spt2.vx = fidMtx.MvxBvx * spt1.vx + fidMtx.MvxBp  * (spt1.pz+spt1.px)
    spt2.pz = fidMtx.MpzBpz * spt1.pz + fidMtx.MpzBvz *  spt2.vz
    spt2.px = fidMtx.MpxBpx * spt1.px + fidMtx.MpxBvx *  spt2.vx
    return nothing
end

"""
    MultiStepForward(pos_rec, srcs, fidMtx)

Multiple time step forward modeling of acoustic wave field, generate one shot gather, support
inject multiple sources simutaneously

# Arguments
* `pos_rec :: Array{Int64,2}`: receivers' location, first column specify vertical index of receivers
second column specify horizontal index.
* `srcs :: Array{Source,1}`: composite type for point source injection.
* `fidMtx :: FidMtx`: composite type of sparse partial differential matrix

# Output
* `shot :: Shot`: composite type of common shot gather
"""
function MultiStepForward(pos::Array{Int64,2}, src::Array{Source,1}, fidMtx::FidMtx; tmax=2.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    if iflag == 1
       zupper =    ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper =    0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dt = fidMtx.dt
    (tl, tu) = SrcRange(src)
    nt = round(Int64, tmax/dt)+1
    shot = InitShot(0, 0, pos, 0.0, nt, dt)
    spt1 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 2)
    AddSources!(spt1, src)
    spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if tl <= (it-1)*dt <= tu   #Add sources
           AddSources!(spt2, src)
        end
        copy(spt1, spt2)
        spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return shot
end

"""
    MultiStepForward(pos_rec, src, fidMtx)

Multiple time step forward modeling of acoustic wave field, generate one shot gather

# Arguments
* `pos_rec :: Array{Int64,2}`: receivers' location, first column specify vertical index of receivers
second column specify horizontal index.
* `src :: Source`: composite type for point source injection.
* `fidMtx :: FidMtx`: composite type of sparse partial differential matrix

# Output
* `shot :: Shot`: composite type of common shot gather
"""
function MultiStepForward(pos::Array{Int64,2}, src::Source, fidMtx::FidMtx; tmax=2.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dt = fidMtx.dt
    nt = round(Int64, tmax/dt)+1
    isz= src.iz; isx= src.ix;
    shot = InitShot(isz, isx, pos, 0.0, nt, dt)
    tl = src.ot; tu = tl + (src.nt-1)*dt;
    spt1 = InitSnapShot(src.iz, src.ix, nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(src.iz, src.ix, nz, nx, ext, iflag, dt, 2)
    AddSource!(spt1, src)
    spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if tl <= (it-1)*dt <= tu   #Add sources
           AddSource!(spt2, src)
        end
        copy(spt1, spt2)
        spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return shot
end

"""
    MultiStepForward(path, srcs, fidMtx)

Multiple time step forward modeling of acoustic wave field, generate one shot gather, support
inject multiple sources simutaneously

# Arguments
* `path :: ASCIIString`: the directory for writing binary wavefiled
* `srcs :: Array{Source,1}`: array of composite source type
* `fidMtx :: FidMtx`: composite type of sparse partial differential matrix

# Output
* `shot :: Shot`: composite type of common shot gather
"""
function MultiStepForward(path::ASCIIString, src::Array{Source,1}, fidMtx::FidMtx; tmax=2.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    dt = fidMtx.dt ;
    (tl, tu) = SrcRange(src)
    nt = round(Int64, tmax/dt)+1
    spt1 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 2)
    AddSources!(spt1, src)
    fid = writeSnapShot(path, spt1)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if tl <= (it-1)*dt <= tu   #Add sources
           AddSources!(spt2, src)
        end
        copy(spt1, spt2)
        writeSnapShot(fid, spt1)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    close(fid)
    return nothing
end

# inject sources as SrcSpt, return shot
function MultiStepForward(pos::Array{Int64,2}, path::ASCIIString, fidMtx::FidMtx; tmax=1.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dt = fidMtx.dt
    nt = round(Int64, tmax/dt+1)
    (nz1, nx1, ot, dt1, nt1) = InfoSrcSpt(path)
    if nz != nz1 || nx != nx1 || dt != dt1
       error("size dismatch")
    end
    shot = InitShot(0, 0, pos, ot, nt, dt)
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, 1)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, 1)
    AddSrcSpt2Spt(spt1, path)
    spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if it <= nt1
           AddSrcSpt2Spt(spt2, path)
        end
        copy(spt1, spt2)
        spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return shot
end

#  forward modeling for source location
function MultiStepForward(pos::Array{Int64,2}, w::Array{Float64,1}, dis::Array{Float64,1}, fidMtx::FidMtx; tmax=1.0, interval=500, print_flag=false)
    nz = fidMtx.nz ; nx    = fidMtx.nx;
    ext= fidMtx.ext; iflag = fidMtx.iflag;
    if iflag == 1
       zupper = ext
       Nz = nz + 2*ext
    elseif iflag == 2
       zupper = 0
       Nz = nz +   ext
    end
    Nx = nx + 2*ext
    dt = fidMtx.dt
    nt = round(Int64, tmax/dt+1)
    ntsrc = length(w)
    shot = InitShot(0, 0, pos, 0.0, nt, dt)
    spt1 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, 1)
    AddSources!(spt1, w, dis)
    spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
    spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, 1)
    for it = 2: nt
        OneStepForward!(spt2, spt1, fidMtx)
        if it <= ntsrc
           AddSources!(spt2, w, dis)
        end
        copy(spt1, spt2)
        spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
        if print_flag
           if mod(it, interval) == 0
              println("the current step: $it")
           end
        end
    end
    return shot
end

# function MultiStepForward!(path::ASCIIString, src::Array{Source,1}, fidMtx::FidMtx; tmax=2.0, samp=1)
#     nz = fidMtx.nz ; nx = fidMtx.nx;
#     ext= fidMtx.ext; iflag = fidMtx.iflag;
#     dt = fidMtx.dt ;
#     (tl, tu) = SrcRange(src)
#     nt = round(Int64, tmax/dt)+1
#     spt1 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 1)
#     spt2 = InitSnapShot(src[1].iz, src[1].ix, nz, nx, ext, iflag, dt, 2)
#     AddSources!(spt1, src)
#     fid = write(path, spt1, 0.0, dt)
#     for it = 2: nt
#         OneStepForward!(spt2, spt1, fidMtx)
#         if tl <= (it-1)*dt <= tu   #Add sources
#            AddSources!(spt2, src)
#         end
#         copy(spt1, spt2)
#         if mod(it-1, samp) == 0
#            write(fid, spt1)
#         end
#         if mod(it, 500) == 0
#            println("the current step: $it")
#         end
#     end
#     close(fid)
#     return nothing
# end

# function MultiStepForward(pos::Array{Int64,2}, path::ASCIIString, fidMtx::FidMtx; tmax=1.0, interval=500, print_flag=false)
#     nz = fidMtx.nz ; nx    = fidMtx.nx;
#     ext= fidMtx.ext; iflag = fidMtx.iflag;
#     if iflag == 1
#        zupper = ext
#        Nz = nz + 2*ext
#     elseif iflag == 2
#        zupper = 0
#        Nz = nx +   ext
#     end
#     Nx = nx + 2*ext
#     dt = fidMtx.dt
#     nt = round(Int64, tmax/dt+1)
#     ntsrc = numberOfSnapShots(path)
#     shot = InitShot(0, 0, pos, 0.0, nt, dt)
#     spt1 = readSnapShot(path, 1)
#     spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
#     spt2 = InitSnapShot(0, 0, nz, nx, ext, iflag, dt, 1)
#     for it = 2: nt
#         OneStepForward!(spt2, spt1, fidMtx)
#         if it <= ntsrc
#            spt1 = readSnapShot(path, it)
#            spt2 = spt2 + spt1
#         end
#         copy(spt1, spt2)
#         spt2shot!(shot, spt1, Nz, Nx, ext, zupper)
#         if print_flag
#            if mod(it, interval) == 0
#               println("the current step: $it")
#            end
#         end
#     end
#     return shot
# end

max_stars_repo_name: WneleiGao/AcousticWaveLab

max_issues_repo_name: WneleiGao/AcousticWaveLab

max_forks_repo_name: WneleiGao/AcousticWaveLab


text:
module FeedbackParticleFilters

using LinearAlgebra
using Distributions
using Statistics
using Random
using PDMats
using ProgressMeter
using StatsBase

import Statistics.mean
import Statistics.cov
import Statistics.var

####################################
########## CORE FUNCTIONS ##########
####################################

include("Basics/Abstractions.jl")
export
    TimeType,
    DiscreteTime,
    ContinuousTime

include("Basics/AbstractModel.jl")
export
    AbstractModel,
    state_type,
    obs_type,
    time_type

include("Basics/AbstractFilteringProblem.jl")
export
    AbstractFilteringProblem,
    state_model,
    obs_model,
    state_dim,
    obs_dim,
    hidden_time_type,
    obs_time_type

include("Basics/AbstractFilterRepresentation.jl")
export
    AbstractFilterRepresentation,
    represented_type,
    dim

include("Basics/AbstractFilteringAlgorithm.jl")
export
    AbstractFilteringAlgorithm,
    AbstractFilterState,
    update!

include("Basics/HiddenStateModel.jl")
export
    HiddenStateModel,
    state_dim,
    initial_condition,
    state_type,
    time_type

include("Basics/ObservationModel.jl")
export
    ObservationModel,
    state_dim,
    obs_dim,
    state_type,
    obs_type,
    time_type,
    emit!

include("Basics/ParametricRepresentation.jl")
export
    ParametricRepresentation,
    parameter_vector

include("Basics/ParticleRepresentation.jl")
export
    ParticleRepresentation,
    get_pos,
    list_of_pos,
    no_of_particles,
    eff_no_of_particles,
    particle_dim,
    dim,
    propagate!

include("Basics/UnweightedParticleRepresentation.jl")
export
    UnweightedParticleRepresentation

include("Basics/WeightedParticleRepresentation.jl")
export
    WeightedParticleRepresentation,
    get_weight,
    list_of_weights,
    sum_of_weights,
    eff_no_of_particles,
    dim,
    resample!

include("Basics/GainEstimation.jl")
export
    GainEquation,
    GainEstimationMethod,
    solve!,
    update!,
    state_dim,
    obs_dim

include("Basics/FilteringProblem.jl")
export
    FilteringProblem

include("Basics/ParticleFlow.jl")
export
    ParticleFlowMethod,
    flow!

####################################
########### STATE MODELS ###########
####################################

include("StateModels/DiffusionStateModel.jl")
export
    DiffusionStateModel,
    ScalarDiffusionStateModel,
    drift,
    diffusion,
    initial_condition,
    state_dim,
    noise_dim,
    initialize,
    drift_function,
    diffusion_function

include("StateModels/LinearDiffusionStateModel.jl")
export
    LinearDiffusionStateModel,
    find_stationary_variance,
    drift,
    diffusion,
    initial_condition,
    initialize,
    state_dim,
    noise_dim

include("StateModels/BrownianMotionTorus.jl")
export
    BrownianMotionTorus,
    BrownianMotionCircle

####################################
######## OBSERVATION MODELS ########
####################################

include("ObservationModels/DiffusionObservationModel.jl")
export
    DiffusionObservationModel,
    ScalarDiffusionObservationModel,
    initial_condition,
    drift,
    state_dim,
    obs_dim,
    noise_dim,
    observation_function

include("ObservationModels/LinearDiffusionObservationModel.jl")
export
    LinearDiffusionObservationModel,
    find_stationary_variance,
    drift,
    state_dim,
    obs_dim,
    noise_dim,
    observation_function

include("ObservationModels/CountingObservationModel.jl")
export
    CountingObservationModel,
    initial_condition,
    drift,
    state_dim,
    obs_dim,
    noise_dim,
    observation_function

####################################
###### FILTER REPRESENTATIONS ######
####################################

include("FilterRepresentations/MultivariateGaussian.jl")
export
    MultivariateGaussian,
    parameter_vector

include("FilterRepresentations/UnweightedParticleEnsemble.jl")
export
    UnweightedParticleEnsemble

include("FilterRepresentations/WeightedParticleEnsemble.jl")
export
    WeightedParticleEnsemble

####################################
########## GAIN EQUATIONS ##########
####################################

include("GainEquations/PoissonEquation.jl")
export
    PoissonEquation,
    Htilde

####################################
##### GAIN ESTIMATION METHODS ######
####################################

include("GainEstimationMethods/ConstantGainApproximation.jl")
export
    ConstantGainApproximation

include("GainEstimationMethods/ConstantGainEKSPF.jl")
export
    ConstantGainEKSPF

include("GainEstimationMethods/SemigroupMethod.jl")
export
    SemigroupMethod

include("GainEstimationMethods/DifferentialRKHSMethod.jl")
export
    DifferentialRKHSMethod

include("GainEstimationMethods/DifferentialRKHSMethodS1.jl")
export
    DifferentialRKHSMethodS1

####################################
####### FILTERING ALGORITHMS #######
####################################

# Kalman-Bucy filter
include("FilteringAlgorithms/KBF.jl")
export
    KBF,
    KBState,
    initialize

# Bootstrap Particle Filter
include("FilteringAlgorithms/BPF.jl")
export
    BPF,
    BPFState,
    update_weights!

# Feedback Particle Filter (original)
include("FilteringAlgorithms/FPF.jl")
export
    FPF,
    FPFState,
    gain_estimation_method,
    propagate!,
    assimilate!,
    heun!,
    gainxerror,
    add_gainxerror!,
    applygain!

# Ensemble Kushner-Stratonovich-Poisson Filter (Venugopal et al.)
include("FilteringAlgorithms/EKSPF.jl")
export
    EKSPF,
    EKSPFState,
    propagate!,
    assimilate!

# Point-process Feedback Particle Filter
include("FilteringAlgorithms/ppFPF.jl")
export
    ppFPF,
    ppFPFState,
    propagate!,
    assimilate!

####################################
############ SIMULATION ############
####################################

include("Simulation/Simulation.jl")
export
    Simulation,
    ContinuousTimeSimulation,
    run!,
    simulate!

# Feedback Particle Filter (original)
include("Simulation/SimulationState.jl")
export
    SimulationState,
    hidden_state,
    obs,
    filter_state,
    cond_mean,
    cond_cov,
    cond_var,
    propagate!

####################################
###### PARTICLE FLOW METHODS #######
####################################

include("ParticleFlowMethods/DeterministicFlow.jl")
export
    DeterministicFlow

end # module


text:
module Plot2D
using TransD_GP, PyPlot, MCMC_Driver, JLD, StatsBase, Statistics
function get_posterior(opt_in::TransD_GP.Options; fdataname="", decimate=1)

    M = TransD_GP.history(opt_in, stat=:fstar)
    n = TransD_GP.history(opt_in, stat=:nodes)
    x_ft = TransD_GP.history(opt_in, stat=:x_ftrain)
    iter = TransD_GP.history(opt_in, stat=:iter)
    misfit = load("misfit_T0_"*fdataname*".jld", "misfit")
    T0loc = load("misfit_T0_"*fdataname*".jld", "T0loc")

    return M[1:decimate:end], n[1:decimate:end],
    x_ft[1:decimate:end], iter[1:decimate:end], misfit[1:decimate:end], T0loc[1:decimate:end]
end

function plot_posterior(X::Array{Float64,1}, Y::Array{Float64,1}, Xall::Array{Float64, 2},
                        opt_in::TransD_GP.Options;
        burnin=50000,
        nbins=50,
        figsize=(8,8),
        fdataname="",
        fontsize=14)

    M, n, x_ft, iter, misfit, T0loc = get_posterior(opt_in, fdataname=fdataname)

    f2, ax2 = plt[:subplots](2,1, sharex=true, figsize=(6,4))
    ax2[1][:plot](iter,n)
    ax2[1][:grid]()
    ax2[2][:plot](iter, misfit)
    ax2[2][:grid]()
    ax2[2][:set_xlabel]("iterations")
    ax2[2][:set_ylabel]("-log L")
    ax2[1][:set_ylabel]("# training")
    MCMC_Driver.nicenup(gcf(), fsize=13)
    savefig(fdataname*"_conv.png", dpi=300)

    from = findfirst(iter.>burnin)

    xall = opt_in.xall
    sq = zeros(size(M[1]))
    mn = zeros(size(M[1]))
    for mm in M[from:end]
        sq += mm.^2
        mn += mm
    end

    N = length(M[from:end])
    mn = mn/N

    f3, ax3 = plt[:subplots](1,2,figsize=(10,5), sharex=true, sharey=true)
    nmodel = from
    im1 = ax3[1][:imshow](reshape(M[nmodel],length(Y), length(X)), extent=[X[1],X[end],Y[end],Y[1]])
    ax3[1][:scatter](x_ft[nmodel][1:n[nmodel],1], x_ft[nmodel][1:n[nmodel],2], s=20, color="black", alpha=0.5)
    cb1 = colorbar(im1, ax=ax3[1])
    im2 = ax3[2][:imshow](reshape(mn,length(Y), length(X)), extent=[X[1],X[end],Y[end],Y[1]])
    cb2 = colorbar(im2, ax=ax3[2])
    MCMC_Driver.nicenup(gcf(), fsize=14)
    savefig(fdataname*"_final_mean.png", dpi=300)

    vr = sq/N - (mn).^2

    x, y = zeros(sum(n[from:end])), zeros(sum(n[from:end]))
    nlast = 0
    for (i, xt) in enumerate(x_ft[from:end])
        x[nlast+1:nlast+n[from-1+i]]  = xt[1:n[from-1+i],1]
        y[nlast+1:nlast+n[from-1+i]] = xt[1:n[from-1+i],2]
        nlast = nlast+n[from-1+i]
    end

    f2, ax2 = plt[:subplots](2,2, figsize=figsize)
    edgesx = LinRange(opt_in.xbounds[1,1], opt_in.xbounds[1,2], nbins+1)
    edgesy = LinRange(opt_in.xbounds[2,1], opt_in.xbounds[2,2], nbins+1)
    h = fit(Histogram, (y,x), (edgesy, edgesx)).weights
    im3 = ax2[1][:imshow](h, extent=[edgesy[1],edgesy[end],edgesx[end],edgesx[1]], aspect="auto")
    # ax2[1][:set_ylabel]("m")
    # ax2[1][:set_xlabel]("m")
    yhist = sum(h,dims=1)[:]
    yhist = yhist/sum(yhist)/diff(edgesy)[1]
    ax2[2][:bar](0.5*(edgesy[2:end]+edgesy[1:end-1]), yhist, width=diff(edgesy)[1])
    ax2[2][:set_xlabel]("y")
    ax2[2][:set_ylabel]("pdf")
    xhist = sum(h,dims=2)[:]
    xhist = xhist/sum(xhist)/diff(edgesx)[1]
    ax2[3][:barh](0.5*(edgesx[2:end]+edgesx[1:end-1]), xhist, height=diff(edgesx)[1])
    ax2[3][:set_xlim](0, 3*maximum(xhist))
    ax2[3][:set_ylim](ax2[1][:get_ylim]())
    ax2[3][:set_ylabel]("x")
    ax2[3][:set_xlabel]("pdf")
    # ax2[3][:xaxis][:set_label_position]("top")
    ax2[2][:get_shared_y_axes]()[:join](ax2[1], ax2[3])
    ax2[2][:get_shared_x_axes]()[:join](ax2[1], ax2[2])
    ax2[2][:set_xlim](ax2[1][:get_xlim]())
    k = fit(Histogram, n[from:end], (opt_in.nmin-0.5:opt_in.nmax+0.5)).weights
    k = k/sum(k)
    @show mean(k)
    ax2[4][:bar](opt_in.nmin:opt_in.nmax, k, width=1)
    ax2[4][:plot]([opt_in.nmin, opt_in.nmax],1/(opt_in.nmax-opt_in.nmin+1)*[1,1],"--k")
    ax2[4][:set_xlabel]("# training points")
    ax2[4][:set_ylabel]("pdf")
    ax2[4][:set_xlim](opt_in.nmin, opt_in.nmax)
    MCMC_Driver.nicenup(gcf(), fsize=fontsize)
    savefig(fdataname*"_post_k.png", dpi=300)

    f,ax = plt[:subplots](1,2,figsize=(10,5), sharex=true, sharey=true)
    im1 = ax[1][:imshow](reshape(log10.(sqrt.(vr)),length(Y), length(X)), extent=[X[1],X[end],Y[end],Y[1]], cmap="gray_r")
    ax[1][:scatter](Xall[1,:], Xall[2,:], s=5, color="red", alpha=0.25)
    cb1 = colorbar(im1, ax=ax[1])
    im2 = ax[2][:imshow](log10.(h), extent=[edgesy[1],edgesy[end],edgesx[end],edgesx[1]], cmap="gray_r")
    cb2 = colorbar(im2, ax=ax[2])
    @info "lala6"
    # ax[1][:set_ylabel]("m")
    # ax[1][:set_xlabel]("m")
    # ax[2][:set_xlabel]("m")
    MCMC_Driver.nicenup(gcf(), fsize=14)
    savefig(fdataname*"_final_sd", dpi=300)

    return
end

end


text:
using AstroTime: Period
using DataInterpolations: AbstractInterpolation
using SmoothingSplines
using LinearAlgebra

import Base: getindex, lastindex, show
import AstroBase: state

export Trajectory, initial, final, state, events, times,
    LogEntry, count_id, id, epoch, detector

struct TrajectoryComponent{T, IntType<:AbstractInterpolation}
    data::Vector{T}
    interp::IntType
    name::Symbol
end

struct LogEntry
    id::Int
    detector::Symbol
    time::Float64
    ep::Epoch
end

id(l::LogEntry) = l.id
detector(l::LogEntry) = l.detector
epoch(l::LogEntry) = l.epoch
count_id(idx, log) = count(x->id(x) == idx, log)

struct Trajectory
    initial
    final
    times::Vector
    splines::Vector{SmoothingSpline}
    # components::Vector{TrajectoryComponents}
    array::Matrix
    events::Vector{LogEntry}
end

initial(tra::Trajectory) = tra.initial
final(tra::Trajectory) = tra.final
events(tra::Trajectory) = tra.events
times(tra::Trajectory) = tra.times

function Trajectory(initial, final, events=LogEntry[])
    Trajectory(initial, final, Vector(undef, 0), SmoothingSpline[], Matrix(undef, 0, 0), events)
end

function Trajectory(initial, final, t, vectors, events=LogEntry[])
    n = length(vectors[1])
    splines = Array{SmoothingSpline}(undef, n)
    arr = permutedims(hcat(vectors...))
    m = size(arr)[2]
    for i = 1:m
        splines[i] = fit(SmoothingSpline, t, arr[:,i], 0.0)
    end
    Trajectory(initial, final, t, splines, arr, events)
end

function show(io::IO, tra::Trajectory)
    println(io, "Trajectory")
    println(io, " Start date: $(epoch(initial(tra)))")
    println(io, " End date:   $(epoch(final(tra)))")
end

function interpolate(tra::Trajectory, time)
    SmoothingSplines.predict.(tra.splines, float(time))
end

function state(tra::Trajectory, time)
    if isempty(tra.times)
        error("Trajectory does not contain data points.")
    end
    rv = interpolate(tra, time)
    r = rv[1:3]
    v = rv[4:6]
    ep1 = epoch(initial(tra)) + time * seconds
    f = frame(initial(tra))
    b = body(initial(tra))
    return State(ep1, r, v, f, b)
end

function state(tra::Trajectory, ep::Epoch)
    ep0 = epoch(initial(tra))
    time = typeof(ep0)(ep) - ep0
    tra(seconds(time))
end

state(tra::Trajectory, time::Period) = state(tra, get(seconds(time)))

(tra::Trajectory)(time) = state(tra, time)

getindex(tra::Trajectory, idx) = state(tra, tra.times[idx])
lastindex(tra::Trajectory) = length(tra.times)


text:
## Cairo Constants ##

# deprecated exports---do not add more!
export CAIRO_FORMAT_ARGB32,
    CAIRO_FORMAT_RGB24,
    CAIRO_FORMAT_A8,
    CAIRO_FORMAT_A1,
    CAIRO_FORMAT_RGB16_565,
    CAIRO_CONTENT_COLOR,
    CAIRO_CONTENT_ALPHA,
    CAIRO_CONTENT_COLOR_ALPHA,
    CAIRO_FILTER_FAST,
    CAIRO_FILTER_GOOD,
    CAIRO_FILTER_BEST,
    CAIRO_FILTER_NEAREST,
    CAIRO_FILTER_BILINEAR,
    CAIRO_FILTER_GAUSSIAN

const format_t = Int32
const FORMAT_INVALID   = Int32(-1)
const FORMAT_ARGB32    = Int32(0)
const FORMAT_RGB24     = Int32(1)
const FORMAT_A8        = Int32(2)
const FORMAT_A1        = Int32(3)
const FORMAT_RGB16_565 = Int32(4)
const FORMAT_RGB30     = Int32(5)

const status_t = Int32
const STATUS_SUCCESS                   = Int32(0)
const STATUS_NO_MEMORY                 = Int32(1)
const STATUS_INVALID_RESTORE           = Int32(2)
const STATUS_INVALID_POP_GROUP         = Int32(3)
const STATUS_NO_CURRENT_POINT          = Int32(4)
const STATUS_INVALID_MATRIX            = Int32(5)
const STATUS_INVALID_STATUS            = Int32(6)
const STATUS_NULL_POINTER              = Int32(7)
const STATUS_INVALID_STRING            = Int32(8)
const STATUS_INVALID_PATH_DATA         = Int32(9)
const STATUS_READ_ERROR                = Int32(10)
const STATUS_WRITE_ERROR               = Int32(11)
const STATUS_SURFACE_FINISHED          = Int32(12)
const STATUS_SURFACE_TYPE_MISMATCH     = Int32(13)
const STATUS_PATTERN_TYPE_MISMATCH     = Int32(14)
const STATUS_INVALID_CONTENT           = Int32(15)
const STATUS_INVALID_FORMAT            = Int32(16)
const STATUS_INVALID_VISUAL            = Int32(17)
const STATUS_FILE_NOT_FOUND            = Int32(17)
const STATUS_INVALID_DASH              = Int32(18)
const STATUS_INVALID_DSC_COMMENT       = Int32(19)
const STATUS_INVALID_INDEX             = Int32(20)
const STATUS_CLIP_NOT_REPRESENTABLE    = Int32(21)
const STATUS_TEMP_FILE_ERROR           = Int32(22)
const STATUS_INVALID_STRIDE            = Int32(23)
const STATUS_FONT_TYPE_MISMATCH        = Int32(24)
const STATUS_USER_FONT_IMMUTABLE       = Int32(25)
const STATUS_USER_FONT_ERROR           = Int32(26)
const STATUS_NEGATIVE_COUNT            = Int32(27)
const STATUS_INVALID_CLUSTERS          = Int32(28)
const STATUS_INVALID_SLANT             = Int32(29)
const STATUS_INVALID_WEIGHT            = Int32(30)
const STATUS_INVALID_SIZE              = Int32(31)
const STATUS_USER_FONT_NOT_IMPLEMENTED = Int32(32)
const STATUS_DEVICE_TYPE_MISMATCH      = Int32(33)
const STATUS_DEVICE_ERROR              = Int32(34)
const STATUS_INVALID_MESH_CONSTRUCTION = Int32(35)
const STATUS_DEVICE_FINISHED           = Int32(36)
const STATUS_LAST_STATUS               = Int32(37)

const surface_t = Int32
const CAIRO_SURFACE_TYPE_IMAGE            = Int32(0)
const CAIRO_SURFACE_TYPE_PDF,             = Int32(1)
const CAIRO_SURFACE_TYPE_PS,              = Int32(2)
const CAIRO_SURFACE_TYPE_XLIB,            = Int32(3)
const CAIRO_SURFACE_TYPE_XCB,             = Int32(4)
const CAIRO_SURFACE_TYPE_GLITZ,           = Int32(5)
const CAIRO_SURFACE_TYPE_QUARTZ,          = Int32(6)
const CAIRO_SURFACE_TYPE_WIN32,           = Int32(7)
const CAIRO_SURFACE_TYPE_BEOS,            = Int32(8)
const CAIRO_SURFACE_TYPE_DIRECTFB,        = Int32(9)
const CAIRO_SURFACE_TYPE_SVG,             = Int32(10)
const CAIRO_SURFACE_TYPE_OS2,             = Int32(11)
const CAIRO_SURFACE_TYPE_WIN32_PRINTING,  = Int32(12)
const CAIRO_SURFACE_TYPE_QUARTZ_IMAGE,    = Int32(13)
const CAIRO_SURFACE_TYPE_SCRIPT,          = Int32(14)
const CAIRO_SURFACE_TYPE_QT,              = Int32(15)
const CAIRO_SURFACE_TYPE_RECORDING,       = Int32(16)
const CAIRO_SURFACE_TYPE_VG,              = Int32(17)
const CAIRO_SURFACE_TYPE_GL,              = Int32(18)
const CAIRO_SURFACE_TYPE_DRM,             = Int32(19)
const CAIRO_SURFACE_TYPE_TEE,             = Int32(20)
const CAIRO_SURFACE_TYPE_XML,             = Int32(21)
const CAIRO_SURFACE_TYPE_SKIA,            = Int32(22)
const CAIRO_SURFACE_TYPE_SUBSURFACE,      = Int32(23)
const CAIRO_SURFACE_TYPE_COGL             = Int32(24)

const font_slant_t = Int32
const FONT_SLANT_NORMAL  = Int32(0)
const FONT_SLANT_ITALIC  = Int32(1)
const FONT_SLANT_OBLIQUE = Int32(2)

const font_weight_t = Int32
const FONT_WEIGHT_NORMAL = Int32(0)
const FONT_WEIGHT_BOLD = Int32(1)

const CONTENT_COLOR = Int32(0x1000)
const CONTENT_ALPHA = Int32(0x2000)
const CONTENT_COLOR_ALPHA = Int32(0x3000)

const FILTER_FAST = 0
const FILTER_GOOD = 1
const FILTER_BEST = 2
const FILTER_NEAREST = 3
const FILTER_BILINEAR = 4
const FILTER_GAUSSIAN = 5

const EXTEND_NONE = 0
const EXTEND_REPEAT = 1
const EXTEND_REFLECT = 2
const EXTEND_PAD = 3

const ANTIALIAS_DEFAULT = 0
const ANTIALIAS_NONE = 1
const ANTIALIAS_GRAY = 2
const ANTIALIAS_SUBPIXEL = 3
# next 3 are not available on all Cairo installations (only since 1.12)
const ANTIALIAS_FAST = 4
const ANTIALIAS_GOOD = 5
const ANTIALIAS_BEST = 6

const OPERATOR_CLEAR = 0
const OPERATOR_SOURCE = 1
const OPERATOR_OVER = 2
const OPERATOR_IN = 3
const OPERATOR_OUT = 4
const OPERATOR_ATOP = 5
const OPERATOR_DEST = 6
const OPERATOR_DEST_OVER = 7
const OPERATOR_DEST_IN = 8
const OPERATOR_DEST_OUT = 9
const OPERATOR_DEST_ATOP = 10
const OPERATOR_XOR = 11
const OPERATOR_ADD = 12
const OPERATOR_SATURATE = 13
const OPERATOR_MULTIPLY = 14
const OPERATOR_SCREEN = 15
const OPERATOR_OVERLAY = 16
const OPERATOR_DARKEN = 17
const OPERATOR_LIGHTEN = 18

## LaTex Token Dicts ##
const _common_token_dict = Dict(
    "\\{"              => "{",
    "\\}"              => "}",
    "\\_"              => "_",
    "\\^"              => "^",
    "\\-"              => "-",

    ## ignore stray brackets
    "{"                => "",
    "}"                => ""
)

const _text_token_dict = Dict(
    ## non-math symbols (p438)
    "\\S"              => "\ua7",
    "\\P"              => "\ub6",
    "\\dag"            => "\u2020",
    "\\ddag"           => "\u2021"
)

const _math_token_dict = Dict(

    "-"                => "\u2212", # minus sign

    ## spacing
    "\\quad"           => "\u2003", # 1 em
    "\\qquad"          => "\u2003\u2003", # 2 em
    "\\,"              => "\u2006", # 3/18 em
    "\\>"              => "\u2005", # 4/18 em
    "\\;"              => "\u2004", # 5/18 em

    ## lowercase greek
    "\\alpha"          => "\u03b1",
    "\\beta"           => "\u03b2",
    "\\gamma"          => "\u03b3",
    "\\delta"          => "\u03b4",
    "\\epsilon"        => "\u03b5",
    "\\varepsilon"     => "\u03f5",
    "\\zeta"           => "\u03b6",
    "\\eta"            => "\u03b7",
    "\\theta"          => "\u03b8",
    "\\vartheta"       => "\u03d1",
    "\\iota"           => "\u03b9",
    "\\kappa"          => "\u03ba",
    "\\lambda"         => "\u03bb",
    "\\mu"             => "\u03bc",
    "\\nu"             => "\u03bd",
    "\\xi"             => "\u03be",
    "\\omicron"        => "\u03bf",
    "\\pi"             => "\u03c0",
    "\\varpi"          => "\u03d6",
    "\\rho"            => "\u03c1",
    "\\varrho"         => "\u03f1",
    "\\sigma"          => "\u03c3",
    "\\varsigma"       => "\u03c2",
    "\\tau"            => "\u03c4",
    "\\upsilon"        => "\u03c5",
    "\\phi"            => "\u03d5",
    "\\varphi"         => "\u03c6",
    "\\chi"            => "\u03c7",
    "\\psi"            => "\u03c8",
    "\\omega"          => "\u03c9",

    ## uppercase greek
    "\\Alpha"           => "\u0391",
    "\\Beta"            => "\u0392",
    "\\Gamma"           => "\u0393",
    "\\Delta"           => "\u0394",
    "\\Epsilon"         => "\u0395",
    "\\Zeta"            => "\u0396",
    "\\Eta"             => "\u0397",
    "\\Theta"           => "\u0398",
    "\\Iota"            => "\u0399",
    "\\Kappa"           => "\u039a",
    "\\Lambda"          => "\u039b",
    "\\Mu"              => "\u039c",
    "\\Nu"              => "\u039d",
    "\\Xi"              => "\u039e",
    "\\Pi"              => "\u03a0",
    "\\Rho"             => "\u03a1",
    "\\Sigma"           => "\u03a3",
    "\\Tau"             => "\u03a4",
    "\\Upsilon"         => "\u03a5",
    "\\Phi"             => "\u03a6",
    "\\Chi"             => "\u03a7",
    "\\Psi"             => "\u03a8",
    "\\Omega"           => "\u03a9",

    ## miscellaneous
    "\\aleph"           => "\u2135",
    "\\hbar"            => "\u210f",
    "\\ell"             => "\u2113",
    "\\wp"              => "\u2118",
    "\\Re"              => "\u211c",
    "\\Im"              => "\u2111",
    "\\partial"         => "\u2202",
    "\\infty"           => "\u221e",
    "\\prime"           => "\u2032",
    "\\emptyset"        => "\u2205",
    "\\nabla"           => "\u2206",
    "\\surd"            => "\u221a",
    "\\top"             => "\u22a4",
    "\\bot"             => "\u22a5",
    "\\|"               => "\u2225",
    "\\angle"           => "\u2220",
    "\\triangle"        => "\u25b3", # == \bigtriangleup
    "\\backslash"       => "\u2216",
    "\\forall"          => "\u2200",
    "\\exists"          => "\u2203",
    "\\neg"             => "\uac",
    "\\flat"            => "\u266d",
    "\\natural"         => "\u266e",
    "\\sharp"           => "\u266f",
    "\\clubsuit"        => "\u2663",
    "\\diamondsuit"     => "\u2662",
    "\\heartsuit"       => "\u2661",
    "\\spadesuit"       => "\u2660",

    ## large operators
    "\\sum"             => "\u2211",
    "\\prod"            => "\u220f",
    "\\coprod"          => "\u2210",
    "\\int"             => "\u222b",
    "\\oint"            => "\u222e",
    "\\bigcap"          => "\u22c2",
    "\\bigcup"          => "\u22c3",
    "\\bigscup"         => "\u2a06",
    "\\bigvee"          => "\u22c1",
    "\\bigwedge"        => "\u22c0",
    "\\bigodot"         => "\u2a00",
    "\\bigotimes"       => "\u2a02",
    "\\bigoplus"        => "\u2a01",
    "\\biguplus"        => "\u2a04",

    ## binary operations
    "\\pm"              => "\ub1",
    "\\mp"              => "\u2213",
    "\\setminus"        => "\u2216",
    "\\cdot"            => "\u22c5",
    "\\times"           => "\ud7",
    "\\ast"             => "\u2217",
    "\\star"            => "\u22c6",
    "\\diamond"         => "\u22c4",
    "\\circ"            => "\u2218",
    "\\bullet"          => "\u2219",
    "\\div"             => "\uf7",
    "\\cap"             => "\u2229",
    "\\cup"             => "\u222a",
    "\\uplus"           => "\u228c", # 228e?
    "\\sqcap"           => "\u2293",
    "\\sqcup"           => "\u2294",
    "\\triangleleft"    => "\u22b2",
    "\\triangleright"   => "\u22b3",
    "\\wr"              => "\u2240",
    "\\bigcirc"         => "\u25cb",
    "\\bigtriangleup"   => "\u25b3", # == \triangle
    "\\bigtriangledown" => "\u25bd",
    "\\vee"             => "\u2228",
    "\\wedge"           => "\u2227",
    "\\oplus"           => "\u2295",
    "\\ominus"          => "\u2296",
    "\\otimes"          => "\u2297",
    "\\oslash"          => "\u2298",
    "\\odot"            => "\u2299",
    "\\dagger"          => "\u2020",
    "\\ddagger"         => "\u2021",
    "\\amalg"           => "\u2210",

    ## relations
    "\\leq"             => "\u2264",
    "\\prec"            => "\u227a",
    "\\preceq"          => "\u227c",
    "\\ll"              => "\u226a",
    "\\subset"          => "\u2282",
    "\\subseteq"        => "\u2286",
    "\\sqsubseteq"      => "\u2291",
    "\\in"              => "\u2208",
    "\\vdash"           => "\u22a2",
    "\\smile"           => "\u2323",
    "\\frown"           => "\u2322",
    "\\geq"             => "\u2265",
    "\\succ"            => "\u227b",
    "\\succeq"          => "\u227d",
    "\\gg"              => "\u226b",
    "\\supset"          => "\u2283",
    "\\supseteq"        => "\u2287",
    "\\sqsupseteq"      => "\u2292",
    "\\ni"              => "\u220b",
    "\\dashv"           => "\u22a3",
    "\\mid"             => "\u2223",
    "\\parallel"        => "\u2225",
    "\\equiv"           => "\u2261",
    "\\sim"             => "\u223c",
    "\\simeq"           => "\u2243",
    "\\asymp"           => "\u224d",
    "\\approx"          => "\u2248",
    "\\cong"            => "\u2245",
    "\\bowtie"          => "\u22c8",
    "\\propto"          => "\u221d",
    "\\models"          => "\u22a7", # 22a8?
    "\\doteq"           => "\u2250",
    "\\perp"            => "\u27c2",

    ## arrows
    "\\leftarrow"       => "\u2190",
    "\\Leftarrow"       => "\u21d0",
    "\\rightarrow"      => "\u2192",
    "\\Rightarrow"      => "\u21d2",
    "\\leftrightarrow"  => "\u2194",
    "\\Leftrightarrow"  => "\u21d4",
    "\\mapsto"          => "\u21a6",
    "\\hookleftarrow"   => "\u21a9",
    "\\leftharpoonup"   => "\u21bc",
    "\\leftharpoondown" => "\u21bd",
    "\\rightleftharpoons" => "\u21cc",
    "\\longleftarrow"   => "\u27f5",
    "\\Longleftarrow"   => "\u27f8",
    "\\longrightarrow"  => "\u27f6",
    "\\Longrightarrow"  => "\u27f9",
    "\\longleftrightarrow" => "\u27f7",
    "\\Longleftrightarrow" => "\u27fa",
    "\\hookrightarrow"  => "\u21aa",
    "\\rightharpoonup"  => "\u21c0",
    "\\rightharpoondown" => "\u21c1",
    "\\uparrow"         => "\u2191",
    "\\Uparrow"         => "\u21d1",
    "\\downarrow"       => "\u2193",
    "\\Downarrow"       => "\u21d3",
    "\\updownarrow"     => "\u2195",
    "\\Updownarrow"     => "\u21d5",
    "\\nearrow"         => "\u2197",
    "\\searrow"         => "\u2198",
    "\\swarrow"         => "\u2199",
    "\\nwarrow"         => "\u2196",

    ## openings
#    "\\lbrack"          => "[",
#    "\\lbrace"          => "{",
    "\\langle"          => "\u27e8",
    "\\lfloor"          => "\u230a",
    "\\lceil"           => "\u2308",

    ## closings
#    "\\rbrack"          => "]",
#    "\\rbrace"          => "}",
    "\\rangle"          => "\u27e9",
    "\\rfloor"          => "\u230b",
    "\\rceil"           => "\u2309",

    ## alternate names
    "\\ne"              => "\u2260",
    "\\neq"             => "\u2260",
    "\\le"              => "\u2264",
    "\\ge"              => "\u2265",
    "\\to"              => "\u2192",
    "\\gets"            => "\u2192",
    "\\owns"            => "\u220b",
    "\\land"            => "\u2227",
    "\\lor"             => "\u2228",
    "\\lnot"            => "\uac",
    "\\vert"            => "\u2223",
    "\\Vert"            => "\u2225",

    ## extensions
    "\\deg"             => "\ub0",
    "\\degr"            => "\ub0",
    "\\degree"          => "\ub0",
    "\\degrees"         => "\ub0",
    "\\arcdeg"          => "\ub0",
    "\\arcmin"          => "\u2032",
    "\\arcsec"          => "\u2033"
)

const CAIRO_FILL_RULE_WINDING = Int32(0);
const CAIRO_FILL_RULE_EVEN_ODD = Int32(1);

const CAIRO_LINE_CAP_BUTT = Int32(0);
const CAIRO_LINE_CAP_ROUND = Int32(1);
const CAIRO_LINE_CAP_SQUARE = Int32(2);

const CAIRO_LINE_JOIN_MITER = Int32(0);
const CAIRO_LINE_JOIN_ROUND = Int32(1);
const CAIRO_LINE_JOIN_BEVEL = Int32(2);

#typedef enum _cairo_path_data_type {
#    CAIRO_PATH_MOVE_TO,
#    CAIRO_PATH_LINE_TO,
#    CAIRO_PATH_CURVE_TO,
#    CAIRO_PATH_CLOSE_PATH
#} cairo_path_data_type_t;

const CAIRO_PATH_MOVE_TO = Int32(0);
const CAIRO_PATH_LINE_TO = Int32(1);
const CAIRO_PATH_CURVE_TO = Int32(2);
const CAIRO_PATH_CLOSE_PATH = Int32(3);


text:
#---------------------------------------------------------------------#
#This code computes the global matrices using DSS
#Written by F.X. Giraldo on April 24, 2019
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function global_matrices(Me,intma,periodicity,Ne,Np,Npoin,DFloat)

    #Initialize Matrices
    M=zeros(DFloat,Npoin,Npoin)

    #Construct Mass and Differentiation Matrices
    for e=1:Ne
        for l=1:Np, k=1:Np
            J=intma[k,l,e]
            for j=1:Np, i=1:Np
                I=intma[i,j,e]
                M[I,J] += Me[i,j,k,l,e]
            end #j,i
        end #k,l
    end #e

    #Periodicity
    for i=1:Npoin
        j=periodicity[i]
        if (i != j)
            for j=1:Npoin
                M[i,j]=DFloat(0.0)
            end
            M[i,i]=DFloat(1.0)
        end
    end

    return (M)
end


text:
module DataIO

using Match
using GZip
using Printf

using SparseArrays

export root
export read_ply, read_obj
export read_mtx, read_features, read_barcodes, read_matrix, expand_matrix

# ------------------------------------------------------------------------
# globals

const root = "$(homedir())/root/data/seqspace"
panic(msg) = error(msg)

# ------------------------------------------------------------------------
# geometric file formats

"""
    const plytype = Dict{String, Type}

Lookup table for keywords of PLY file format matched to type information.
"""
const plytype = Dict{String, Type}(
    "char"   => Int8,    "int8"    => Int8,
    "short"  => Int16,   "int16"   => Int16,
    "int"    => Int32,   "int32"   => Int32,
    "long"   => Int64,   "int64"   => Int64,

    "uchar"   => UInt8,  "uint8"   => UInt8,
    "ushort"  => UInt16, "uint16"  => UInt16,
    "uint"    => UInt32, "uint32"  => UInt32,
    "ulong"   => UInt64, "uint64"  => UInt64,

    "float"  => Float32, "float32" => Float32,
    "double" => Float64, "float64" => Float64
)

readwords(io::IO) = split(readline(io))

"""
    struct PLYProp
        name :: Symbol
        type :: Type
        len  :: Union{Type,Nothing}
    end

Data structure that encapsulates a property of a PLY file.
If `len` is not nothing, it is assumed to be a collection of properties, i.e. a list.
"""
struct PLYProp
    name :: Symbol
    type :: Type
    len  :: Union{Type,Nothing}
end

PLYProp(word) = if word[1] == "list"
                    PLYProp(Symbol(word[4]), plytype[word[3]], plytype[word[2]])
                else
                    PLYProp(Symbol(word[2]), plytype[word[1]], nothing)
                end

"""
    fill!(prop::Array{PLYProp}, io::IO)

Read and parse a single line from `io` stream.
Fill the interpreted property into `prop`.
"""
function fill!(prop::Array{PLYProp}, io::IO)
    word = readwords(io)
    while word[1] == "property" && length(word) >= 3
        push!(prop, PLYProp(word[2:end]))
        word = readwords(io)
    end

    return word
end

"""
    readprops(io::IO, props::Array{PLYProp,1})

Read a collection of properties from `io` stream into `props` array.
Return a named tuple of parsed properties.
"""
function readprops(io::IO, props::Array{PLYProp,1})
    data = (p) -> begin
        if isnothing(p.len)
            return read(io,p.type)
        else
            len = read(io,p.len)
            return [read(io,p.type) for _ in 1:len]
        end
    end

    return (; (prop.name => data(prop) for prop in props)...)
end

"""
    read_ply(io::IO)

Read an oriented mesh from `io` stream assuming bytes are formatted in PLY file format.
"""
function read_ply(io::IO)
    # magic
    line = readline(io)
    if line != "ply"
        panic("can not interpret byte stream as PLY data")
    end

    # file format
    word = readwords(io)
    if word[1] != "format" || word[3] != "1.0"
        panic("unrecognized format of PLY file")
    end
    # FIXME: we just ignore this for now and assume little endian binary
    fmt = Symbol(word[2])

    # data layout
    n₀, n₂ = 0, 0
    field₀ = PLYProp[] # list of fields of vertices
    field₂ = PLYProp[] # list of fields of faces

    @label getline # ----------------------------
    word = readwords(io)

    @label doline
    if word[1] == "end_header"
        @goto endloop
    end

    if word[1] != "element" || length(word) != 3
        @goto getline
    end

    @match word[2] begin
        "vertex" => begin
            n₀   = parse(Int, word[3])
            word = fill!(field₀, io)
            @goto doline
        end
        "face" => begin
            n₂   = parse(Int, word[3])
            word = fill!(field₂, io)
            @goto doline
        end
    end

    @goto  getline

    @label endloop # ----------------------------

    return (
        verts = [readprops(io,field₀) for _ in 1:n₀],
        faces = [readprops(io,field₂) for _ in 1:n₂]
    )
end

"""
    read_obj(io::IO)

Read an oriented mesh from `io` stream assuming bytes are formatted in OBJ file format.
"""
function read_obj(io::IO)
    bufᵥ = []
    bufₙ = []
    bufₜ = []

    for ln in eachline(io)
        if startswith(ln, "#") || length(ln) == 0
            continue
        end

        w = split(ln)
        @match w[1] begin
            "v"  => push!(bufᵥ, [parse(Float32, x) for x in w[2:end]])
            "vₙ" => push!(bufₙ, [parse(Float32, x) for x in w[2:end]])
            "f"  => push!(bufₜ, [parse(Int32, split(x,"//")[1]) for x in w[2:end]])
             _   => continue
        end
    end

    return (
        verts = hcat(bufᵥ...)',
        norms = hcat(bufₙ...)',
        faces = hcat(bufₜ...)',
    )
end

# ------------------------------------------------------------------------
# matrix market exchange format

"""
    read_mtx(io::IO)

Read a matrix from `io` stream assuming bytes are formatted in matrix exchance file format.
"""
function read_mtx(io::IO)
    header = readwords(io)
    if header[1] != "%%MatrixMarket"
        panic("can not interpret given IO stream as matrix market exchange")
    end

    T =
    if header[4] == "integer"
        Int
    elseif header[4] == "real"
        Float64
    elseif header[4] == "complex"
        Complex
    elseif header[4] == "pattern"
        Bool
    else
        panic("unrecognized element type '$(header[4])' in matrix file")
    end

    D =
    if header[2] == "matrix"
        2
    elseif header[2] == "vector"
        1
    else
        panic("unrecognized data type '$(header[2])' in matrix file")
    end

    data =
    if header[3] == "array"
        line = readwords(io)
        while line[1][1] == '%'
            line = readwords(io)
        end

        m, n  = parse.(Int, line[1:2])
        array = Array{T,D}(undef, m, n)

        for i in 1:m
            for j in 1:n
                line = getline(io)
                array[i,j] = parse(T, line)
            end
        end

        array
    elseif header[3] == "coordinate"
        line = readwords(io)
        while line[1][1] == '%'
            line = readwords(io)
        end

        m, n, numel = parse.(Int, line[1:3])
        row = Array{Int,1}(undef, numel)
        col = Array{Int,1}(undef, numel)
        val = Array{T,1}(undef, numel)

        for (i,line) in enumerate(eachline(io))
            word = split(line)
            row[i] = parse(Int, word[1])
            col[i] = parse(Int, word[2])
            val[i] = parse(T, word[3])
        end

        sparse(row, col, val, m, n)
    else
        panic("unrecognized data format '$(header[3])'")
    end
end

# NOTE: assumes Gene Expression is the only feature type
read_features(io::IO) = [split(line)[2] for line in eachline(io)]
read_barcodes(io::IO) = [line for line in eachline(io)]

# ------------------------------------------------------------------------
# scRNAseq file formats

"""
    read_matrix(io::IO; type=Float64, named_cols=false, named_rows=false, start_cols=1, start_rows=1)

Read a white-space delimited matrix from `io` stream of element type `type`.
If `named_cols` is true, the first line is assumed to be column labels.
If `named_rows` is true, the first column of each row is assumed to be denote the label.
"""
function read_matrix(io::IO; type=Float64, named_cols=false, named_rows=false, start_cols=1, start_rows=1)
    cols = named_cols ? readwords(io)[start_cols:end] : nothing

    x = position(io)
    ncols = if named_cols
        length(cols)
    else
        words = readwords(io)
        n = length(words) - (named_rows ? 1 : 0)
        seek(io, x)

        n
    end
    nrows = sum(1 for row in eachline(io))
    seek(io, x)

    data = Array{type, 2}(undef, nrows, ncols)
    rows = named_rows ? Array{String, 1}(undef, nrows) : nothing

    for (i,row) in enumerate(eachline(io))
        entry = split(row)
        if named_rows
            rows[i] = entry[1]
            entry   = entry[2:end]
        end
        data[i,:] = parse.(type, entry)
    end

    return data, rows, cols
end

"""
    expand_matrix(io::IO, dir::String)

Read a matrix in MTX format from `io` stream and dump as 3 text files into directory `dir`.
"""
function expand_matrix(io::IO, dir::String)
    if isdir(dir)
        error("directory exists")
    end
    mkdir(dir)

    data, genes, barcodes = read_matrix(io; type=Int, named_cols=true, start_cols=2, named_rows=true)

    open("$dir/barcodes.tsv", "w") do fd
        write(fd, barcodes[1])
        for barcode ∈ barcodes[2:end]
            write(fd, "\n", barcode)
        end
    end

    feature(i, gene) = ((@sprintf "Gene%04d" i), "\t", gene, "\t", "Gene Expression")
    open("$dir/features.tsv", "w") do fd
        write(fd, feature(1, genes[1])...)
        for (i, gene) in enumerate(genes[2:end])
            write(fd, "\n", feature(i+1, gene)...)
        end
    end

    nnz(x) = sum(x[:] .> 0)
    open("$dir/matrix.mtx", "w") do fd
        write(fd, "%%MatrixMarket matrix coordinate integer general", "\n")
        write(fd, "%", "\n")
        write(fd, string(size(data,1)), " ", string(size(data,2)), " ", string(nnz(data)))
        for j in 1:size(data,2)
            for i in findall(data[:,j] .> 0)
                write(fd, "\n", string(i), " ", string(j), " ", string(data[i,j]))
            end
        end
    end
end

# ------------------------------------------------------------------------
# point of entry for testing

function test()
    mesh = open("$root/gut/mesh_apical_stab_000153.ply") do io
        read_ply(io)
    end

    bdntp, genes = GZip.open("$root/dvex/bdtnp.txt.gz") do io
        data, _, cols = read_matrix(io; named_cols=true)
        return data, cols
    end

    scrna, genes = GZip.open("$root/dvex/dge_normalized.txt.gz") do io
        data, rows, _ = read_matrix(io; named_rows=true, named_cols=true)
        return data, rows
    end
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
function Maxwell(
            sim::TSIM;
            ky::Number=0,
            kz::Number=0,
            m::Int = 1) where TSIM<:Simulation{1}

    kx, ka, kb, kc = 0, 0, 0, 0
    curlcurl = sim.curlcurl(ky,kz)

    Σ01,Σ02,Σ11,Σ12,Σ2 = _compute_Σs(sim,ky,kz)
    Σs = (Σ01 + Σ11, Σ02 + Σ12, Σ2)
    f = sim.self_energy.f
    D² = Σs[1]*f[1](0,ka,ky,kz) + Σs[2]*f[2](0,ka,ky,kz) + Σs[3] + curlcurl

    αε = _compute_αε(sim)
    αεpFχ = copy(αε)
    Fχ = spdiagm(0=>zeros(ComplexF64,size(αε,1)))

    Fχs = Vector{SparseMatrixCSC{ComplexF64,Int}}(undef,m)
    dFχdψr = zeros(ComplexF64,3length(sim),m,m)
    dFχdψi = copy(dFχdψr)
    dFχdω = copy(dFχdψi)
    dFχdϕ = copy(dFχdω)
    for μ ∈ eachindex(Fχs) Fχs[μ] = spdiagm(0=>zeros(ComplexF64,size(αε,1))) end
    Maxwell{1,typeof(Σs),TSIM}(D²+I,D²,αεpFχ,curlcurl,Σs,αε,Fχ,Fχs,dFχdψr,dFχdψi,dFχdω,dFχdϕ,sim,kx,ky,kz,ka,kb,kc)
end

@inline function (m::Maxwell{1})(ω::Number,args...)
    ka, ky, kz = m.ka, m.ky, m.kz
    f = m.sim.self_energy.f

    f1 = f[1](ω,ka,ky,kz)
    f2 = f[2](ω,ka,ky,kz)

    rows = rowvals(m.D²)
    vals = nonzeros(m.D²)
    _, n = size(m.D²)
    @inbounds for i ∈ 1:n
        col = i
        @fastmath @inbounds @simd for j ∈ nzrange(m.D², i)
            row = rows[j]
            vals[j] = m.Σs[1][row,col]*f1 + m.Σs[2][row,col]*f2 + m.Σs[3][row,col] + m.curlcurl[row,col]
        end
    end

    maxwell_susceptability!(m,ω,args...)
    ω² = ω^2

    rows = rowvals(m.A)
    vals = nonzeros(m.A)
    @inbounds for col ∈ 1:n
        @fastmath @inbounds @simd for j ∈ nzrange(m.A, col)
            row = rows[j]
            vals[j] = m.D²[row,col] - ω²*m.αεpFχ[row,col]
        end
    end
    return m.A
end

################################################################################

# linear
@inline function maxwell_susceptability!(m::Maxwell{1},ω::Number)
    sim = m.sim
    @fastmath @inbounds @simd for i ∈ eachindex(m.Fχ.nzval)
        j = mod1(i,length(sim))
        χ = susceptability(sim.χ[j], ω)
        m.Fχ.nzval[i] = sim.F[j]*χ
    end
    @fastmath @inbounds @simd for i ∈ eachindex(m.αε.nzval) m.αεpFχ.nzval[i] = m.αε.nzval[i] + m.Fχ.nzval[i] end
    return nothing
end

# nonlinear + local jacobian
@inline function maxwell_susceptability!(m::Maxwell{1},ω,ωs::Vector,ψs::ElectricField)
    sim = m.sim
    foreach(d->susceptability(d.χ,ω,ωs,ψs), sim.dispersive_domains)
    @inbounds for μ ∈ eachindex(m.Fχs)
        @inbounds for i ∈ eachindex(m.Fχ.nzval)
            j = mod1(i,length(sim))
            if typeof(sim.χ[j])<:TwoLevelSystem
                χ::TwoLevelSystem = sim.χ[j]
                μ==1 ? m.Fχ.nzval[i] = sim.F[j]*χ.chi[i] : nothing
                m.Fχs[μ].nzval[i] = sim.F[j]*χ.chis[i,μ]
                @inbounds @simd for ν ∈ eachindex(m.Fχs)
                    m.dFχdψr[i,μ,ν] = sim.F[j]*χ.dχdψr[i,μ,ν]
                    m.dFχdψi[i,μ,ν] = sim.F[j]*χ.dχdψi[i,μ,ν]
                    m.dFχdω[i,μ,ν] = sim.F[j]*χ.dχdω[i,μ,ν]
                    m.dFχdϕ[i,μ,ν] = sim.F[j]*χ.dχdϕ[i,μ,ν]
                end
            end
        end
    end
    @fastmath @inbounds @simd for i ∈ eachindex(m.αε.nzval) m.αεpFχ.nzval[i] = m.αε.nzval[i] + m.Fχ.nzval[i] end
    return nothing
end

################################################################################

@inline function _compute_Σs(sim::Simulation{1},ky::Number,kz::Number)

    ky², kz², kykz = ky^2, kz^2, ky*kz

    Σ0 = sim.self_energy.Σ0
    Σ1 = sim.self_energy.Σ1
    Σ2 = sim.self_energy.Σ2

    I1 = sparse([2,3],[2,3],[1,1],3,3)
    I2 = sparse([2,3,1,1],[1,1,2,3],[ky,kz,ky,kz],3,3)

    Σ01 = kron(I1,Σ0[1])
    Σ02 = kron(I1,Σ0[2])
    Σ11 = kron(I2,Σ1[1])
    Σ12 = kron(I2,Σ1[2])
    Σ2 = kron(sparse([1,2,3,2,3],[1,2,2,3,3],[ky²+kz²,kz²,kykz,kykz,ky²],3,3),Σ2)
    return Σ01,Σ02,Σ11,Σ12,Σ2
end

_compute_αε(sim::Simulation{1})::SparseMatrixCSC{ComplexF64,Int} =
    kron(I3,spdiagm(0=>sim.α[1].*sim.ε))


text:
using GLMakie
GLMakie.activate!()

module Perlin
    const PERMS1 = [151, 160, 137, 91, 90, 15, 131, 13, 201, 95, 96, 53, 194, 233, 7,
        225, 140, 36, 103, 30, 69, 142, 8, 99, 37, 240, 21, 10, 23, 190, 6, 148, 247,
        120, 234, 75, 0, 26, 197, 62, 94, 252, 219, 203, 117, 35, 11, 32, 57, 177, 33,
        88, 237, 149, 56, 87, 174, 20, 125, 136, 171, 168, 68, 175, 74, 165, 71, 134,
        139, 48, 27, 166, 77, 146, 158, 231, 83, 111, 229, 122, 60, 211, 133, 230, 220,
        105, 92, 41, 55, 46, 245, 40, 244, 102, 143, 54, 65, 25, 63, 161, 1, 216, 80,
        73, 209, 76, 132, 187, 208, 89, 18, 169, 200, 196, 135, 130, 116, 188, 159, 86,
        164, 100, 109, 198, 173, 186, 3, 64, 52, 217, 226, 250, 124, 123, 5, 202, 38,
        147, 118, 126, 255, 82, 85, 212, 207, 206, 59, 227, 47, 16, 58, 17, 182, 189,
        28, 42, 223, 183, 170, 213, 119, 248, 152, 2, 44, 154, 163, 70, 221, 153, 101,
        155, 167, 43, 172, 9, 129, 22, 39, 253, 19, 98, 108, 110, 79, 113, 224, 232,
        178, 185, 112, 104, 218, 246, 97, 228, 251, 34, 242, 193, 238, 210, 144, 12,
        191, 179, 162, 241, 81, 51, 145, 235, 249, 14, 239, 107, 49, 192, 214, 31, 181,
        199, 106, 157, 184, 84, 204, 176, 115, 121, 50, 45, 127, 4, 150, 254, 138, 236,
        205, 93, 222, 114, 67, 29, 24, 72, 243, 141, 128, 195, 78, 66, 215, 61, 156, 180]

    const PERMS = [PERMS1; PERMS1]

    function fade(x)
        return 6x^5 - 15x^4 + 10x^3
    end

    function smoothstep(a0, a1, x)
        a0 != a1 || throw(ArgumentError("Arguments a0 and a1 cannot be equal"))
        x = clamp((x - a0) / (a1 - a0), 0.0, 1.0)
        return 3x^2 - 2x^3
    end

    function smootherstep(a0, a1, x)
        a0 != a1 || throw(ArgumentError("Arguments a0 and a1 cannot be equal"))
        x = clamp((x - a0) / (a1 - a0), 0.0, 1.0)
        return fade(x)
    end

    # Might refactor this down to a simple lerp
    function interpolate(a0, a1, w, method::Symbol=:linear)
        0 ≤ w ≤ 1 || throw(ArgumentError("Expected 0 ≤ w ≤ 1, got $w"))
        if method ∉ [:linear, :smoothstep, :smootherstep]
            throw(ArgumentError("Unrecognized interpolation method $method"))
        end
        x = method == :linear ? w : eval(method)(a0, a1, w)
        return a0 + x * (a1 - a0)
    end

    function gradient(hash::Int, x, y, z)
        h = hash & 15
        u = h < 8 ? x : y
        v = h < 4 ? y : h == 12 || h == 14 ? x : z
        return ((h & 1) == 0 ? u : -u) + ((h & 2) == 0 ? v : -v)
    end

    function perlin(x, y, z)
        xi = round(Int, x) & 255 + 1
        yi = round(Int, y) & 255 + 1
        zi = round(Int, z) & 255 + 1

        xf = first(modf(x))
        yf = first(modf(y))
        zf = first(modf(z))

        u = fade(xf)
        v = fade(yf)
        w = fade(zf)

        aaa = PERMS[PERMS[PERMS[xi] + yi] + zi]
        aba = PERMS[PERMS[PERMS[xi] + yi + 1] + zi]
        aab = PERMS[PERMS[PERMS[xi] + yi] + zi + 1]
        abb = PERMS[PERMS[PERMS[xi] + yi + 1] + zi + 1]
        baa = PERMS[PERMS[PERMS[xi + 1] + yi] + zi]
        bba = PERMS[PERMS[PERMS[xi + 1] + yi + 1] + zi]
        bab = PERMS[PERMS[PERMS[xi + 1] + yi] + zi + 1]
        bbb = PERMS[PERMS[PERMS[xi + 1] + yi + 1] + zi + 1]

        x1 = interpolate(gradient(aaa, xf, yf, zf), gradient(baa, xf - 1, yf, zf), u)
        x2 = interpolate(gradient(aba, xf, yf - 1, zf), gradient(bba, xf - 1, yf - 1, zf), u)
        y1 = interpolate(x1, x2, v)

        x1 = interpolate(gradient(aab, xf, yf, zf - 1), gradient(bab, xf - 1, yf, zf - 1), u)
        x2 = interpolate(gradient(abb, xf, yf - 1, zf - 1), gradient(bbb, xf - 1, yf - 1, zf - 1), u)
        y2 = interpolate(x1, x2, v)

        return (interpolate(y1, y2, w) + 1) / 2
    end

    function octaveperlin(x, y, z, octaves::Int, persistence)
        total = 0.0
        frequency = 1.0
        amplitude = 1.0
        maxval = 0.0
        for i = 1:octaves
            total += perlin(x * frequency, y * frequency, z * frequency) * amplitude
            maxval += amplitude
            amplitude *= persistence
            frequency *= 2
        end
        return total / maxval
    end
end

begin
    noise = Perlin.octaveperlin.(range(0, 2, length = 300), range(0, 2, length = 300)', 6.0, 10, 0.5)
    noise = noise .+ [x * y for x in range(0, 1, length = 300), y in range(0, 1, length = 300)]
    s = Scene(show_axis = false, camera = cam3d!)
    surf = surface!(s, 0..3, 0..3, noise, colormap = :grays)
    center!(s)
    resize!(s, 1708, 619)

    s.camera.view[] = Float32[0.9745395 -0.22421555 2.8969254f-8 -0.84478927; 0.031377867 0.13638215 0.9901593 -1.1058276; -0.22200912 -0.9649494 0.13994522 -4.068703; 0.0 0.0 0.0 1.0];
    s.camera.projection[] = Float32[4.655257 0.0 0.0 0.0; 0.0 12.8452015 0.0 0.0; 0.0 0.0 -1.0689656 -0.92784494; 0.0 0.0 -1.0 0.0]
    s.camera.projectionview[] = Float32[4.536732 -1.043781 1.3485933f-7 -3.9327114; 0.403055 1.7518562 12.718796 -14.204578; 0.23732011 1.0314977 -0.14959662 3.4214587; 0.22200912 0.9649494 -0.13994522 4.068703];
    s.backgroundcolor = :gray80
    s
end


text:
# ====================================================================
# Component structure
#
mutable struct Gaussian_1D <: AbstractComponent
    norm::Parameter
    center::Parameter
    sigma::Parameter

    function Gaussian_1D(norm::Number, center::Number, sigma::Number)
        @assert norm  > 0
        @assert sigma > 0

        out = new(Parameter(norm), Parameter(center), Parameter(sigma))
        out.norm.low = 0
        out.sigma.low = 0
        return out
    end
end

mutable struct Gaussian_2D <: AbstractComponent
    norm::Parameter
    centerX::Parameter
    centerY::Parameter
    sigmaX::Parameter
    sigmaY::Parameter
    angle::Parameter

    function Gaussian_2D(norm::Number, centerX::Number, centerY::Number, sigmaX::Number, sigmaY::Number, angle::Number)
        @assert norm   > 0
        @assert sigmaX > 0
        @assert sigmaY > 0

        out = new(Parameter(norm), Parameter(centerX), Parameter(centerY), Parameter(sigmaX), Parameter(sigmaY), Parameter(angle))
        out.norm.low = 0
        out.sigmaX.low = 0
        out.sigmaY.low = 0
        return out
    end
end

# ====================================================================
Gaussian(norm, center, sigma) = Gaussian_1D(norm, center, sigma)
Gaussian(norm, centerX, centerY, sigmaX, sigmaY, angle) = Gaussian_2D(norm, centerX, centerY, sigmaX, sigmaY, angle)
function Gaussian(norm, centerX, centerY, sigma)
    out = Gaussian_2D(norm, centerX, centerY, sigma, sigma, 0.)
    out.sigmaX.fixed = true
    out.angle.fixed = true
    return out
end

# ====================================================================
# Evaluate component
function evaluate!(buffer, comp::Gaussian_1D, x::AbstractDomain{1},
                   norm, center, sigma)
    @. (buffer = exp( ((x - center) / sigma)^2. / (-2.)) /
        2.5066282746310002 / sigma * norm) # sqrt(2pi) = 2.5066282746310002
end

function evaluate!(buffer, comp::Gaussian_2D, domain::AbstractDomain{2},
                   norm, centerX, centerY, sigmaX, sigmaY, angle)
    angle *= -pi / 180.
    a =  (cos(angle) / sigmaX)^2 / 2  +  (sin(angle) / sigmaY)^2 / 2
    b = -sin(2angle) / sigmaX^2  / 2  +  sin(2angle) / sigmaY^2  / 2
    c =  (sin(angle) / sigmaX)^2 / 2  +  (cos(angle) / sigmaY)^2 / 2

    x = domain[1]
    y = domain[2]

    @. (buffer = norm *
        exp(
            -(
                a * (x - centerX)^2. +
                b * (x - centerX) * (y - centerY) +
                c *                 (y - centerY)^2.
            )
        ) / 6.283185307179586 / sigmaX / sigmaY) # 2pi = 6.283185307179586
end


text:
export flattenednames

flattenednames(::Type{NamedTuple{(),Tuple{}}}; prefix="") = ()
flattenednames(::Type{T}; prefix="") where {T<:Real} = (prefix, )
flattenednames(::Type{T}; prefix="") where {T<:SVector} = ntuple(i -> "$prefix[$i]", length(T))
function flattenednames(::Type{T}; prefix="") where {T<:SHermitianCompact}
  N = size(T,1)
  [["$prefix[$i,$j]" for i = j:N] for j = 1:N] |> Iterators.flatten |> collect
end
function flattenednames(::Type{T}; prefix="") where {T<:NamedTuple}
  map(1:fieldcount(T)) do i
    Ti = fieldtype(T,i)
    name = fieldname(T,i)
    flattenednames(Ti, prefix= prefix=="" ? string(name) : string(prefix, '.', name))
  end |> Iterators.flatten |> collect
end

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
abstract PartialDN
immutable PartialD{order} <: PartialDN
  n :: Int64
  function PartialD()
    @assert eltype(order) <: Int
    return new(order)
  end
end
export PartialD

abstract PartialDStrategy
abstract PartialDImplicitStrategy <: PartialDStrategy
abstract PartialDExplicitStrategy <: PartialDStrategy

abstract BoundaryStrategy
immutable NaiveBoundary <: BoundaryStrategy end
immutable NaivePeriodicShared <: BoundaryStrategy
  n_shared :: Int64
  function NaivePeriodicShared(n :: Int64 = 3)
    @assert n > 0
    new(n)
  end
end
export NaiveBoundary
export NaivePeriodicShared

immutable PartialD_OP{T1 <: PartialDN, T2<: PartialDStrategy, T3 <: BoundaryStrategy} <: oneDoperator
  pfunctor :: T1
  sfunctor :: T2
  bfunctor :: T3
  function PartialD_OP()
    new(T1(),T2(),T3())
  end
  function PartialD_OP(x1 :: T1, x2 :: T2, x3:: T3)
    new(x1,x2,x3)
  end
  function PartialD_OP(x1 :: Tuple, x2 :: Tuple, x3:: Tuple)
    new(T1(x1...), T2(x2...), T3(x3...))
  end
end
export PartialD_OP

include("implicit/partialD.jl")


text:
# support for device-side exceptions

## exception type

struct KernelException <: Exception
    dev::CuDevice
end

function Base.showerror(io::IO, err::KernelException)
    print(io, "KernelException: exception thrown during kernel execution on device $(name(err.dev))")
end

## exception codegen

# emit a global variable for storing the current exception status
#
# since we don't actually support globals, access to this variable is done by calling the
# julia_exception_flag function (lowered here to actual accesses of the variable)
function emit_exception_flag!(mod::LLVM.Module)
    # add the global variable if it doesn't exist yet (might have come from the runtime)
    # TODO: get!
    if !haskey(globals(mod), "exception_flag")
        T_ptr = convert(LLVMType, Ptr{Cvoid})
        gv = GlobalVariable(mod, T_ptr, "exception_flag")
        initializer!(gv, LLVM.ConstantInt(T_ptr, 0))
        linkage!(gv, LLVM.API.LLVMWeakAnyLinkage)
        extinit!(gv, true)

        set_used!(mod, gv)
    else
        gv = globals(mod)["exception_flag"]
    end

    # lower uses of the getter
    if haskey(functions(mod), "julia_exception_flag")
        buf_getter = functions(mod)["julia_exception_flag"]
        @assert return_type(eltype(llvmtype(buf_getter))) == eltype(llvmtype(gv))

        # find uses
        worklist = Vector{LLVM.CallInst}()
        for use in uses(buf_getter)
            call = user(use)::LLVM.CallInst
            push!(worklist, call)
        end

        # replace uses by a load from the global variable
        for call in worklist
            Builder(JuliaContext()) do builder
                position!(builder, call)
                ptr = load!(builder, gv)
                replace_uses!(call, ptr)
            end
            unsafe_delete!(LLVM.parent(call), call)
        end
    end
end

## exception handling

const exception_flags = Dict{CuContext, Mem.HostBuffer}()

# create a CPU/GPU exception flag for error signalling, and put it in the module
#
# also see compiler/irgen.jl::emit_exception_flag!
function create_exceptions!(mod::CuModule)
    if VERSION >= v"1.5" && !(v"1.6-" <= VERSION < v"1.6.0-DEV.90")
        flag_ptr = CuGlobal{Ptr{Cvoid}}(mod, "exception_flag")
        exception_flag = get!(exception_flags, mod.ctx,
                            Mem.alloc(Mem.Host, sizeof(Int), Mem.HOSTALLOC_DEVICEMAP))
        flag_ptr[] = reinterpret(Ptr{Cvoid}, convert(CuPtr{Cvoid}, exception_flag))
    else
        ctx = mod.ctx
        try
            flag_ptr = CuGlobal{Ptr{Cvoid}}(mod, "exception_flag")
            exception_flag = get!(exception_flags, ctx, Mem.alloc(Mem.Host, sizeof(Int),
                                Mem.HOSTALLOC_DEVICEMAP))
            flag_ptr[] = reinterpret(Ptr{Cvoid}, convert(CuPtr{Cvoid}, exception_flag))
        catch err
            (isa(err, CuError) && err.code == ERROR_NOT_FOUND) || rethrow()
        end
    end

    return
end

# check the exception flags on every API call, similarly to how CUDA handles errors
function check_exceptions()
    for (ctx,buf) in exception_flags
        if isvalid(ctx)
            ptr = convert(Ptr{Int}, buf)
            flag = unsafe_load(ptr)
            if flag != 0
                unsafe_store!(ptr, 0)
                dev = device(ctx)
                throw(KernelException(dev))
            end
        end
    end
    return
end


text:
module EliminateGraphs

using Graphs
import Graphs: vertices, edges, neighbors, ne, nv, degree

include("utils.jl")
include("Core.jl")
include("Graph.jl")
include("iterset.jl")
include("generateset.jl")
include("degrees.jl")
include("eliminate.jl")
include("graphlib.jl")
include("algorithms/algorithms.jl")
end # module


text:
#=
Suppose we represent our file system by a string in the following manner:

The string "dir\n\tsubdir1\n\tsubdir2\n\t\tfile.ext" represents:

dir
    subdir1
    subdir2
        file.ext

The directory dir contains an empty sub-directory subdir1 and a sub-directory subdir2 containing a file file.ext.

The string "dir\n\tsubdir1\n\t\tfile1.ext\n\t\tsubsubdir1\n\tsubdir2\n\t\tsubsubdir2\n\t\t\tfile2.ext" represents:

dir
    subdir1
        file1.ext
        subsubdir1
    subdir2
        subsubdir2
            file2.ext

The directory dir contains two sub-directories subdir1 and subdir2. subdir1 contains a file file1.ext and an empty second-level sub-directory subsubdir1. subdir2 contains a second-level sub-directory subsubdir2 containing a file file2.ext.

We are interested in finding the longest (number of characters) absolute path to a file within our file system. For example, in the second example above, the longest absolute path is "dir/subdir2/subsubdir2/file2.ext", and its length is 32 (not including the double quotes).

Given a string representing the file system in the above format, return the length of the longest absolute path to a file in the abstracted file system. If there is no file in the system, return 0.

Note:

The name of a file contains at least a period and an extension.

The name of a directory or sub-directory will not contain a period.
=#

# Create a data structure for the file system

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
struct MaxSafeAccel
    dt::Float64
end

"""
Calculate the maximum safe acceleration that will allow the car to avoid a collision if the car in front slams on its brakes
"""
function max_safe_acc(safety::MaxSafeAccel, c::Context, idx::Int)
    dt = safety.dt

    fidx = fore_index(c, idx)
    if fidx == 0
        return Inf
    else
        n_brake_acc = nullable_max_safe_acc(headway(c, idx),
                                            vel(state(c, idx)),
                                            vel(state(c, fidx)),
                                            brake_limit(vehicle(c, idx)), dt)
        return get(n_brake_acc, -brake_limit(vehicle(c, idx)))
    end
end

"""
Return max_safe_acc or an empty Nullable if the discriminant is negative.
"""
function nullable_max_safe_acc(gap, v_behind, v_ahead, braking_limit, dt)
    bp = braking_limit
    v = v_behind
    vo = v_ahead
    g = gap
    # VVV see mathematica notebook
    discriminant = 8.*g*bp + bp^2*dt^2 - 4.*bp*dt*v + 4.*vo^2
    if discriminant >= 0.0
        return Nullable{Float64}(- (bp*dt + 2.*v - sqrt(discriminant)) / (2.*dt))
    else
        return Nullable{Float64}()
    end
end

function is_safe(safety::MaxSafeAccel, c::Context, idx::Int, a::AccelSpeed)
    return a.long_accel <= max_safe_acc(safety, c, idx)
end

struct RoadEdges
    dt::Float64
end

function is_safe(safety::RoadEdges, c::Context, idx::Int, a::AccelSpeed)
    y = pos(state(c, idx))[2]
    w = width(vehicle(c, idx))
    rw = c.roadway
    next_y = y + a.lat_speed*safety.dt
    return next_y - w/2 >= -rw.lanewidth/2 && next_y + w/2 <= (rw.n_lanes-0.5)*rw.lanewidth
end


text:
using Documenter
using GadgetUnits

makedocs(
    sitename="GadgetUnits",
    format=Documenter.HTML(),
    modules=[GadgetUnits],
    pages = [
            "Table of Contents" => "index.md",
            "Install"           => "install.md",
            "Basic Unit Conversion" => "conversion_structs.md",
            "Cosmological conversions" => "cosmo.md",
            "API reference"     => "api.md"
            ]
        )

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    repo = "github.com/LudwigBoess/GadgetUnits.jl.git"
)


text:
struct Row{names, T <: Tuple}
    data::T
end

@inline Row(in::Pair{<:Label}...) = Row{getnames(in)}(getvalues(in))
@inline Row{names}(in::Tuple) where {names} = Row{names, typeof(in)}(in)

@inline getindex(row::Row{names}, l::Label) where {names} = row.data[label_to_index(names, l)]

function show(io::IO, row::Row{names}) where names
    print(io, "Row(")
    for i ∈ 1:length(names)
        print(io, "l", '"', names[i], '"', " => ", row.data[i])
        if i != length(names)
            print(io, ", ")
        end
    end
    print(io, ")")
end

@inline names(::Row{_names}) where {_names} = _names
@inline labels(::Row{names}) where {names} = labels(names...)

Base.hash(r::Row, h::UInt) = hash(r.data, h)


text:
module Kernels

using Test
using CellwiseValues

export NumberKernel
export ArrayKernel
export NumberKernelFromFunction
export ArrayKernelFromBroadcastedFunction
export compute_type
export compute_value
export compute_value!
export compute_size
export compute_ndim
export test_number_kernel
export test_array_kernel

# Interfaces

abstract type NumberKernel end

function compute_type(::NumberKernel,::Vararg{<:Type})::Type{<:NumberLike}
  @abstractmethod
end

function compute_value(::NumberKernel,::Vararg)::NumberLike
  @abstractmethod
end

abstract type ArrayKernel end

function compute_type(::ArrayKernel,::Vararg{<:Type})::Type{<:NumberLike}
  @abstractmethod
end

function compute_ndim(::ArrayKernel,::Vararg{Int})::Int
  @abstractmethod
end

function compute_size(::ArrayKernel,::Vararg{<:NTuple})::NTuple
  @abstractmethod
end

function compute_value!(::AbstractArray,::NumberKernel,::Vararg)
  @abstractmethod
end

# Testers

function test_number_kernel(k::NumberKernel,o::T,i::Vararg) where T
  t = [ typeof(ii) for ii in i ]
  S = compute_type(k,t...)
  @test S == T
  r = compute_value(k,i...)
  @test r == o
end

function test_array_kernel(
  k::ArrayKernel,o::AbstractArray{T,N},i::Vararg) where {T,N}
  t = [ eltype(ii) for ii in i ]
  S = compute_type(k,t...)
  @test S == T
  s = [ _size_for_broadcast(ii) for ii in i ]
  si = compute_size(k,s...)
  @test size(o) == si
  d = [length(si) for si in s]
  n = compute_ndim(k,d...)
  @test n == length(si)
  r = Array{T,N}(undef,si)
  compute_value!(r,k,i...)
  @test r == o
end

@inline _size_for_broadcast(a) = size(a)

@inline _size_for_broadcast(a::NumberLike) = ()

# Implementations

struct NumberKernelFromFunction{F<:Function} <: NumberKernel
  fun::F
end

function compute_type(
  self::NumberKernelFromFunction,types::Vararg{<:Type})
  T = Base._return_type(self.fun,types)
  @assert T <: NumberLike
  T
end

function compute_value(self::NumberKernelFromFunction,args::Vararg)
  self.fun(args...)
end

struct ArrayKernelFromBroadcastedFunction{F<:Function} <: ArrayKernel
  fun::F
end

function compute_type(
  self::ArrayKernelFromBroadcastedFunction,etypes::Vararg{<:Type})
  Base._return_type(self.fun,etypes)
end

function compute_size(::ArrayKernelFromBroadcastedFunction,s::Vararg{<:NTuple})
  Base.Broadcast.broadcast_shape(s...)
end

function compute_value!(
  v::AbstractArray, k::ArrayKernelFromBroadcastedFunction, a::Vararg)
  broadcast!(k.fun,v,a...)
end

function compute_ndim(::ArrayKernelFromBroadcastedFunction,d::Vararg{Int})
  n = 0
  for di in d
    n = max(n,di)
  end
  n
end

end # module Kernels


text:
# example matrices courtesy of Pablo San-Jose

psj_1_A = sparse(
[9, 10, 9, 10, 11, 14, 12, 13, 12, 13, 11, 14, 15, 16, 15, 16, 1, 9, 10, 16, 2, 9, 10, 11, 3, 10, 11, 14, 4, 12, 13, 15, 5, 12, 13, 14, 6, 11, 13, 14, 7, 12, 15, 16, 8, 9, 15, 16],
[1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16],
Complex{Float64}[3.7796350217469814 + 0.0im, -3.312563559813306 + 0.0im, -3.3125635598133054 + 0.0im, 3.779635021746982 + 0.0im, 3.7796350217469814 + 0.0im, -3.312563559813306 + 0.0im, 3.779635021746982 + 0.0im, -3.3125635598133054 + 0.0im, -3.312563559813306 + 0.0im, 3.7796350217469814 + 0.0im, -3.3125635598133054 + 0.0im, 3.779635021746982 + 0.0im, 3.7796350217469814 + 0.0im, -3.312563559813306 + 0.0im, -3.3125635598133054 + 0.0im, 3.779635021746982 + 0.0im, 1.0 + 0.0im, 6.418270043493963 + 0.0im, -6.625127119626612 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -6.625127119626611 + 0.0im, 6.418270043493964 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, 6.418270043493963 + 0.0im, -6.625127119626612 + 0.0im, 1.0 + 0.0im, 6.418270043493964 + 0.0im, -6.625127119626611 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -6.625127119626612 + 0.0im, 6.418270043493963 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -6.625127119626611 + 0.0im, -1.0 + 0.0im, 6.418270043493964 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, 6.418270043493963 + 0.0im, -6.625127119626612 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -6.625127119626611 + 0.0im, 6.418270043493964 + 0.0im],16,16)
psj_1_B = sparse(
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 10, 11, 14, 12, 13, 12, 13, 11, 14, 15, 16, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16],
Complex{Float64}[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, -3.7796350217469814 + 0.0im, 3.3125635598133054 + 0.0im, 3.312563559813306 + 0.0im, -3.779635021746982 + 0.0im, -3.7796350217469814 + 0.0im, 3.3125635598133054 + 0.0im, -3.779635021746982 + 0.0im, 3.312563559813306 + 0.0im, 3.3125635598133054 + 0.0im, -3.7796350217469814 + 0.0im, 3.312563559813306 + 0.0im, -3.779635021746982 + 0.0im, -3.7796350217469814 + 0.0im, 3.3125635598133054 + 0.0im, 3.312563559813306 + 0.0im, -3.779635021746982 + 0.0im],16,16)
psj_2_A = sparse(
[20, 22, 24, 13, 15, 17, 1, 13, 19, 20, 2, 14, 20, 21, 3, 15, 21, 22, 4, 16, 22, 23, 5, 17, 23, 24, 6, 18, 19, 24, 7, 13, 18, 19, 8, 13, 14, 20, 9, 14, 15, 21, 10, 15, 16, 22, 11, 16, 17, 23, 12, 17, 18, 24],
[2, 4, 6, 7, 9, 11, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16, 17, 17, 17, 17, 18, 18, 18, 18, 19, 19, 19, 19, 20, 20, 20, 20, 21, 21, 21, 21, 22, 22, 22, 22, 23, 23, 23, 23, 24, 24, 24, 24],
Complex{Float64}[-1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -2.62 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -1.0 + 0.0im, -2.62 + 0.0im],24,24)
psj_2_B = sparse(
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 19, 21, 23, 14, 16, 18],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 15, 17, 20, 22, 24],
Complex{Float64}[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im],24,24)
psj_3_A = sparse(
[9, 10, 9, 10, 11, 14, 12, 13, 12, 13, 11, 14, 15, 16, 15, 16, 1, 9, 10, 16, 2, 9, 10, 11, 3, 10, 11, 14, 4, 12, 13, 15, 5, 12, 13, 14, 6, 11, 13, 14, 7, 12, 15, 16, 8, 9, 15, 16],
[1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10, 11, 11, 11, 11, 12, 12, 12, 12, 13, 13, 13, 13, 14, 14, 14, 14, 15, 15, 15, 15, 16, 16, 16, 16],
Complex{Float64}[0.33748484079831426 + 0.0im, -0.10323794456968927 + 0.0im, -0.10323794456968927 + 0.0im, 0.3374848407983142 + 0.0im, 0.33748484079831426 + 0.0im, -0.10323794456968927 + 0.0im, 0.3374848407983142 + 0.0im, -0.10323794456968927 + 0.0im, -0.10323794456968927 + 0.0im, 0.33748484079831426 + 0.0im, -0.10323794456968927 + 0.0im, 0.3374848407983142 + 0.0im, 0.33748484079831426 + 0.0im, -0.10323794456968927 + 0.0im, -0.10323794456968927 + 0.0im, 0.3374848407983142 + 0.0im, 1.0 + 0.0im, -2.5940303184033713 + 0.0im, -0.20647588913937853 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -0.20647588913937853 + 0.0im, -2.5940303184033713 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -2.5940303184033713 + 0.0im, -0.20647588913937853 + 0.0im, 1.0 + 0.0im, -2.5940303184033713 + 0.0im, -0.20647588913937853 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -0.20647588913937853 + 0.0im, -2.5940303184033713 + 0.0im, -1.0 + 0.0im, 1.0 + 0.0im, -0.20647588913937853 + 0.0im, -1.0 + 0.0im, -2.5940303184033713 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -2.5940303184033713 + 0.0im, -0.20647588913937853 + 0.0im, 1.0 + 0.0im, -1.0 + 0.0im, -0.20647588913937853 + 0.0im, -2.5940303184033713 + 0.0im],16,16)
psj_3_B = sparse(
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 9, 10, 11, 14, 12, 13, 12, 13, 11, 14, 15, 16, 15, 16],
[1, 2, 3, 4, 5, 6, 7, 8, 9, 9, 10, 10, 11, 11, 12, 12, 13, 13, 14, 14, 15, 15, 16, 16],
Complex{Float64}[1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, 1.0 + 0.0im, -0.33748484079831426 + 0.0im, 0.10323794456968927 + 0.0im, 0.10323794456968927 + 0.0im, -0.3374848407983142 + 0.0im, -0.33748484079831426 + 0.0im, 0.10323794456968927 + 0.0im, -0.3374848407983142 + 0.0im, 0.10323794456968927 + 0.0im, 0.10323794456968927 + 0.0im, -0.33748484079831426 + 0.0im, 0.10323794456968927 + 0.0im, -0.3374848407983142 + 0.0im, -0.33748484079831426 + 0.0im, 0.10323794456968927 + 0.0im, 0.10323794456968927 + 0.0im, -0.3374848407983142 + 0.0im],16,16)


text:
function truncate!(P::CuVector{Float64};
                   kwargs...)::Tuple{Float64,Float64,CuVector{Float64}}
  maxdim::Int = min(get(kwargs,:maxdim,length(P)), length(P))
  mindim::Int = min(get(kwargs,:mindim,1), maxdim)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  origm = length(P)
  docut = 0.0
  maxP  = maximum(P)
  if maxP == 0.0
    P = CuArrays.zeros(Float64, 1)
    return 0.,0.,P
  end
  if origm==1
    docut = maxP/2
    return 0., docut, P[1:1]
  end

  @timeit "setup rP" begin
      #Zero out any negative weight
      #neg_z_f = (!signbit(x) ? x : 0.0)
      rP = map(x -> !signbit(x) ? x : 0.0, P)
      n = origm
      truncerr = 0.0
      if n > maxdim
          truncerr = sum(rP[1:n-maxdim])
          n = maxdim
      end
  end
  @timeit "handle cutoff" begin
      if absoluteCutoff
        #Test if individual prob. weights fall below cutoff
        #rather than using *sum* of discarded weights
        sub_arr = rP .- cutoff
        err_rP  = sub_arr ./ abs.(sub_arr)
        flags   = reinterpret(Float64, (signbit.(err_rP) .<< 1 .& 2) .<< 61)
        cut_ind = CuArrays.CUBLAS.iamax(err_rP .* flags) - 1
        n = min(maxdim, length(P) - cut_ind)
        n = max(n, mindim)
        truncerr += sum(rP[cut_ind+1:end])
      else
        scale = 1.0
        @timeit "find scale" begin
            if doRelCutoff
              scale = sum(P)
              scale = scale > 0.0 ? scale : 1.0
            end
        end

        #Continue truncating until *sum* of discarded probability
        #weight reaches cutoff reached (or m==mindim)
        sub_arr = rP .+ truncerr .- cutoff*scale
        err_rP  = sub_arr ./ abs.(sub_arr)
        flags   = reinterpret(Float64, (signbit.(err_rP) .<< 1 .& 2) .<< 61)
        cut_ind = CuArrays.CUBLAS.iamax(err_rP .* flags) - 1
        if cut_ind > 0
            truncerr += sum(rP[cut_ind+1:end])
            n = min(maxdim, length(P) - cut_ind)
            n = max(n, mindim)
            if scale==0.0
              truncerr = 0.0
            else
              truncerr /= scale
            end
        else # all are above cutoff
            truncerr += sum(rP[1:maxdim])
            n = min(maxdim, length(P) - cut_ind)
            n = max(n, mindim)
            if scale==0.0
              truncerr = 0.0
            else
              truncerr /= scale
            end
        end
      end
  end
  if n < 1
    n = 1
  end
  if n < origm
    hP = collect(P)
    docut = (hP[n]+hP[n+1])/2
    if abs(hP[n]-hP[n+1]) < 1E-3*hP[n]
      docut += 1E-3*hP[n]
    end
  end
  @timeit "setup return" begin
      rinds = 1:n
      rrP   = P[rinds]
  end
  return truncerr,docut,rrP
end