text:
"""
The MSA module of MIToS has utilities for working with Multiple Sequence Alignments of protein Sequences (MSA).

**Features**

- Read and write MSAs in `Stockholm`, `FASTA` or `Raw` format
- Handle MSA annotations
- Edit the MSA, e.g. delete columns or sequences, change sequence order, shuffling...
- Keep track of positions and annotations after modifications on the MSA
- Describe a MSA, e.g. mean percent identity, sequence coverage, gap percentage...

```julia
using MIToS.MSA
```
"""
module MSA

using DataStructures        # OrderedDicts for Annotations
using AutoHashEquals        # Annotations, Clusters
using NamedArrays           # Col and Seq names, basic sequence/MSA object
using FastaIO               # FastaReader (fast)
using Random                # GLOBAL_RNG, shuffle!, rand, Sampler
using Dates                 # Dates.now()
using PairwiseListMatrices  # Percent Identity Matrices
using Clustering            # Used for sequence clustering: ClusteringResult
using StatsBase             # Weights for clustering
using RecipesBase           # Plots for MSAs
using TranscodingStreams    # To solve MethodError seek(::TranscodingStream, ::Int)
using MIToS.Utils

import Clustering: ClusteringResult, nclusters, counts, assignments

export  # Residue
        Residue,
        GAP, XAA,
        @res_str,
        # Alphabet
        ResidueAlphabet,
        GappedAlphabet, UngappedAlphabet, ReducedAlphabet,
        getnamedict,
        # ThreeLetters
        residue2three, three2residue,
        # Annotations
        Annotations,
        # filtersequences!,
        ncolumns,
        filtercolumns!,
        getannotfile,  getannotcolumn,  getannotsequence,  getannotresidue,
        setannotfile!, setannotcolumn!, setannotsequence!, setannotresidue!,
        annotate_modification!, delete_annotated_modifications!, printmodifications,
        # MultipleSequenceAlignment
        AbstractAlignedObject,
        AbstractMultipleSequenceAlignment,
        AbstractAlignedSequence,
        MultipleSequenceAlignment, AnnotatedMultipleSequenceAlignment,
        AlignedSequence, AnnotatedAlignedSequence,
        AnnotatedAlignedObject, UnannotatedAlignedObject,
        NamedResidueMatrix,
        annotations,
        namedmatrix,
        nsequences,
        getresidues, getsequence, getresiduesequences,
        stringsequence,
        getcolumnmapping, getsequencemapping,
        sequencenames, columnnames, # TO DO: sequencenames!(...)
        # MSAStats
        gapfraction,
        residuefraction,
        coverage,
        columngapfraction,
        # MSAEditing
        filtersequences, filtersequences!,
        filtercolumns, filtercolumns!,
        swapsequences!,
        setreference!,
        adjustreference, adjustreference!,
        gapstrip, gapstrip!,
        # GeneralParserMethods
        deletefullgapcolumns, deletefullgapcolumns!,
        # Raw
        Raw,
        # Stockholm
        Stockholm,
        # FASTA
        FASTA,
        # NBRF/PIR
        PIR,
        # PLM
        sequencepairsmatrix, columnpairsmatrix,
        # Identity
        percentidentity, meanpercentidentity, percentsimilarity,
        # Clusters
        ClusteringResult, # from Clustering.jl
        nclusters, counts, assignments, # from Clustering.jl
        NoClustering, Clusters,
        getweight, nelements,
        # Hobohm
        hobohmI,
        # Imported from Base or StdLib (and exported for docs)
        names,
        parse,
        isvalid,
        rand,
        shuffle,
        shuffle!

include("Residues.jl")
include("Alphabet.jl")
include("ThreeLetters.jl")
include("Annotations.jl")
include("MultipleSequenceAlignment.jl")
include("MSAStats.jl")
include("MSAEditing.jl")
include("GeneralParserMethods.jl")
include("Raw.jl")
include("Stockholm.jl")
include("FASTA.jl")
include("PIR.jl")
include("Shuffle.jl")
include("PLM.jl")
include("Identity.jl")
include("Clusters.jl")
include("Hobohm.jl")
include("Plots.jl")

end


text:
# Test error-handling.
#

@testset "error-handling" begin
    @testset "in DualTensor creation" begin
        try
            DualTensor( rand(3, 3), (rand(3, 3), (rand(3, 2))) )
        catch err
            @test err isa DimensionMismatch
        end
    end
end

ext: jl


text:
using Test
using CUDA
using Random
using JSON

include("../brains/continuous_time_rnn.jl")

function kernel_test_initialize(individuals, brains)

    initialize(brains, individuals)

    sync_threads()

    reset(brains)

    sync_threads()

    return

end

function kernel_test_brain_step(input_all, output_all, brains)

    threadID = threadIdx().x
    blockID = blockIdx().x

    offset_shared_memory = 0

    input = @cuDynamicSharedMem(Float32, brains.input_size, offset_shared_memory)
    offset_shared_memory += sizeof(input)

    output = @cuDynamicSharedMem(Float32, brains.output_size, offset_shared_memory)
    offset_shared_memory += sizeof(output)

    sync_threads()

    # Load inputs to shared memory
    if threadID <= brains.input_size
        input[threadID] = input_all[threadID, blockID]
    end

    sync_threads()

    step(brains, input, output, offset_shared_memory)

    sync_threads()

    # Load outputs from shared memory
    if threadID <= brains.output_size
        output_all[threadID, blockID] = output[threadID]
    end

    return
end

@testset "Brain Continuous Time RNN" begin

    for differential_eq in ["LiHoChow2005", "NaturalNet"]
        for set_principle_diagonal_elements_of_W_negative in [false, true]

            config_brain = OrderedDict()
            config_brain["delta_t"] = 0.05
            config_brain["number_neurons"] = 50
            config_brain["differential_equation"] = differential_eq
            config_brain["clipping_range_min"] = -1.0
            config_brain["clipping_range_max"] = 1.0
            config_brain["set_principle_diagonal_elements_of_W_negative"] = set_principle_diagonal_elements_of_W_negative
            config_brain["alpha"] = 0.1

            number_inputs = 30
            number_outputs = 15
            number_individuals = 100

            number_time_steps = 200

            # Initialize brains
            brains = ContinuousTimeRNN(config_brain, number_inputs, number_outputs, number_individuals)

            individual_size = get_individual_size(brains)

            @test individual_size == number_inputs * brains.number_neurons + brains.number_neurons * brains.number_neurons + number_outputs * brains.number_neurons

            v_size = brains.number_neurons * brains.input_size
            w_size = brains.number_neurons * brains.number_neurons
            t_size = brains.output_size * brains.number_neurons

            individuals = rand(number_individuals, individual_size)
            individuals_gpu = CuArray(individuals)

            @cuda threads = brains.number_neurons blocks = number_individuals kernel_test_initialize(individuals_gpu, brains)

            CUDA.synchronize()

            V = zeros(brains.number_neurons, brains.input_size, number_individuals)
            W = zeros(brains.number_neurons, brains.number_neurons, number_individuals)
            T = zeros(brains.output_size, brains.number_neurons, number_individuals)

            for j = 1:number_individuals

                V[:, :, j] = reshape(view(individuals, j, 1:v_size), (brains.number_neurons, brains.input_size))
                W[:, :, j] = reshape(view(individuals, j, v_size+1:v_size+w_size), (brains.number_neurons, brains.number_neurons))
                T[:, :, j] = reshape(view(individuals, j, v_size+w_size+1:v_size+w_size+t_size), (brains.output_size, brains.number_neurons))

                if brains.set_principle_diagonal_elements_of_W_negative == true
                    for i = 1:brains.number_neurons
                        W[i, i, j] = -abs(W[i, i, j])
                    end
                end

            end

            # Test brain initialization
            @test V ≈ Array(brains.V) atol = 0.00001
            @test W ≈ Array(brains.W) atol = 0.00001
            @test T ≈ Array(brains.T) atol = 0.00001

            x = zeros(brains.number_neurons, number_individuals)

            for i = 1:number_time_steps

                input = randn(number_inputs, number_individuals)
                input_gpu = CuArray(input)

                output = zeros(number_outputs, number_individuals)
                output_gpu = CuArray(output)

                shared_memory = get_memory_requirements(brains) + sizeof(Float32) * (brains.input_size + brains.output_size)

                CUDA.@cuda threads = brains.number_neurons blocks = number_individuals shmem = shared_memory kernel_test_brain_step(input_gpu, output_gpu, brains)
                CUDA.synchronize()

                for j = 1:number_individuals

                    # Differential Equation
                    if brains.differential_equation == NaturalNet
                        dx_dt = -brains.alpha * x[:, j] + W[:, :, j] * map(tanh, x[:, j]) + V[:, :, j] * input[:, j]
                    elseif brains.differential_equation == LiHoChow2005
                        dx_dt = -brains.alpha * x[:, j] + W[:, :, j] * map(tanh, (x[:, j] + V[:, :, j] * input[:, j]))
                    end

                    # Euler forward discretization
                    x[:, j] += brains.delta_t * dx_dt

                    # Clip x to state boundaries
                    x[:, j] = clamp.(x[:, j], brains.clipping_range_min, brains.clipping_range_max)

                    # Calculate outputs
                    y = map(tanh, T[:, :, j] * x)

                    @test x[:, j] ≈ Array(brains.x[:, j]) rtol = 0.001
                    @test y[:, j] ≈ Array(output_gpu[:, j]) rtol = 0.001

                end

            end

        end
    end
end


text:
using Test

@time @testset "Mesh" begin include("MeshTests.jl") end
@time @testset "Quadrature" begin include("QuadTests.jl") end

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
#RSDeltaSigma: Text generation
#-------------------------------------------------------------------------------

const ORDINAL_STRING = [
	"Zeroth", "First", "Second", "Third", "Fourth", "Fifth",
	"Sixth", "Seventh", "Eighth", "Nineth", "Tenth",
	"Twelfth", "Thirteenth",
]

function ds_orderString(n::Int, capitalize::Bool=false)
	idx = n+1
	if idx in keys(ORDINAL_STRING)
		s = ORDINAL_STRING[idx]
		capitalize || (s = lowercase(s))
	else
		s = "$(n)th"
	end
	return s
end

function str_modulatortype(dsm::AbstractDSM)
	typestr = (0==dsm.f0) ? "Lowpass" : "Bandpass"
	if isquadrature(dsm)
		typestr = "Quadrature " * typestr
	end
	return ds_orderString(dsm.order, true) * "-Order " * typestr
end

#Last line


text:
video = Video(800, 800)

@testset "Shorthands" begin
    function ground(args...)
        background("white")
        sethue("black")
    end

    Background(1:90, ground)

    line = Object(JLine(Point(100, -250)))
    line1 = Object(JLine(Point(-350, 350), Point(100, -250)))

    circle0 = Object(JCircle(45, action = :fill))
    circle1 = Object(JCircle(O, 100, action = :stroke))
    circle2 = Object(JCircle(Point(200, -200), Point(190, -190), action = :fill))
    circle3 = Object(JCircle(-350, -250, 20, action = :fill))

    box = Object(JBox(Point(-200, 150), Point(10, -150), color = "red", action = :stroke))
    box1 = Object(JBox([Point(20, -150), Point(100, 250)], color = "blue"))
    box2 = Object(JBox(Point(-250, -100), 200, 200, color = "yellow", action = :fill))
    box3 = Object(JBox(250, 300, 40, 40, color = "yellow", action = :fill))
    box4 = Object(JBox(O, 150, 150, 5.0, color = "black", action = :stroke))

    rect = Object(JRect(250, 15, 30, 55, color = "orange", action = :fill))
    rect1 = Object(JRect(Point(250, -150), 30, 55, color = "orange", action = :fill))

    ellipse = Object(JEllipse(-50, 15, 40, 25, color = "blue"))
    ellipse0 = Object(JEllipse(Point(-50, 15), 45, 30, color = "blue"))
    ellipse1 = Object(
        JEllipse(Point(-50, 15), Point(150, -50), 70, color = "red", action = :stroke),
    )
    ellipse2 = Object(
        JEllipse(
            Point(-50, 15),
            Point(50, -50),
            Point(-100, 290),
            color = "red",
            action = :stroke,
        ),
    )

    star = Object(JStar(Point(-50, 15), 30, color = "red", action = :stroke))
    star1 = Object(JStar(50, -105, 45, color = "red", action = :fill))

    poly0 = Object(JPoly(ngon(O, 150, 3, -π / 2, vertices = true)))
    poly1 = Object(JPoly(ngon(O, 250, 4, π / 2, vertices = true)), action = :fill)

    somepath = Object(@JShape begin
        sethue("blue")
        circle.(ngon(O, 150, 3, -π / 2, vertices = true), 20, :fill)
    end)

    action = :stroke
    somepath1 = Object(
        @JShape begin
            sethue(color)
            poly(ngon(O, 400, 11, 5, vertices = true), action, close = true)
        end action = action color = "red" radius = 8
    )

    act!(
        [
            line,
            circle0,
            circle1,
            circle2,
            circle3,
            box,
            box1,
            box2,
            box3,
            rect,
            rect1,
            ellipse,
            ellipse0,
            ellipse1,
            ellipse2,
            star,
            star1,
            poly0,
            poly1,
        ],
        Action(20:30, anim_translate(Point(-150, 150))),
    )

    act!([circle0], Action(1:50, change(:radius, 45 => 10)))
    act!(
        [
            line,
            line1,
            circle0,
            circle1,
            box,
            box1,
            box2,
            box3,
            box4,
            rect,
            rect1,
            ellipse,
            ellipse0,
            ellipse1,
            ellipse2,
            star,
            star1,
            poly0,
            poly1,
        ],
        [
            Action(40:51, change(:color, "red")),
            Action(52:61, change(:color, colorant"red" => colorant"blue")),
        ],
    )
    act!(
        [
            line,
            line1,
            circle0,
            circle1,
            rect,
            rect1,
            ellipse,
            ellipse0,
            ellipse1,
            ellipse2,
            star,
            star1,
            poly0,
            poly1,
        ],
        Action(71:90, change(:linewidth, 1 => 5)),
    )

    render(video; tempdirectory = "images", pathname = "shorthands.gif")

    for i in ["01", 20, 21, 39, 40, 59, 65, 85]
        @test_reference "refs/shorthands$i.png" load("images/00000000$i.png")
        @test isfile("shorthands.gif")
    end

    for i in 1:90
        rm("images/$(lpad(i, 10, "0")).png")
    end

    rm("images/palette.png")
    rm("shorthands.gif")
end


text:
# Autogenerated wrapper script for CPUInfo_jll for i686-w64-mingw32
export libcpuinfo

JLLWrappers.@generate_wrapper_header("CPUInfo")
JLLWrappers.@declare_library_product(libcpuinfo, "libcpuinfo.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libcpuinfo,
        "bin\\libcpuinfo.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
mutable struct RTBasis{T,M,P} <: Space{T}
  geo::M
  fns::Vector{Vector{Shape{T}}}
  pos::Vector{P}
end

RTBasis(geo, fns) = RTBasis(geo, fns, Vector{vertextype(geo)}(undef,length(fns)))

positions(rt) = rt.pos
refspace(space::RTBasis{T}) where {T} = RTRefSpace{T}()
subset(rt::RTBasis,I) = RTBasis(rt.geo, rt.fns[I], rt.pos[I])

mutable struct ValDiv end

"""
    raviartthomas(mesh, cellpairs::Array{Int,2})

Constructs the RT basis on the input `mesh`. The i-th RT basis function will
    represent a current distribution flowing from cell `cellpairs[1,i]` to
    `cellpairs[2,i]` on the mesh.

Returns an object of type `RTBasis`, which comprises both the mesh and pairs of
    Shape objects which corresponds to the cell pairs, containing the necsessary
    coefficients and indices to compute the exact basis functions when required
    by the solver.
"""
function raviartthomas(mesh, cellpairs::Array{Int,2})

    # combine now the pairs of monopolar RWGs in div-conforming RWGs
    numpairs = size(cellpairs,2)
    functions = Vector{Vector{Shape{Float64}}}(undef,numpairs)
    positions = Vector{vertextype(mesh)}(undef,numpairs)
    Cells = cells(mesh)
    for i in 1:numpairs
        if cellpairs[2,i] > 0
            c1 = cellpairs[1,i]; cell1 = Cells[c1] #mesh.faces[c1]
            c2 = cellpairs[2,i]; cell2 = Cells[c2] #mesh.faces[c2]
            e1, e2 = getcommonedge(cell1, cell2)
            functions[i] = [
              Shape(c1, abs(e1), +1.0),
              Shape(c2, abs(e2), -1.0)]
            isct = intersect(cell1, cell2)
            @assert length(isct) == 2
            @assert !(cell1[abs(e1)] in isct)
            @assert !(cell2[abs(e2)] in isct)

            ctr1 = cartesian(center(chart(mesh, cell1)))
            ctr2 = cartesian(center(chart(mesh, cell2)))
            positions[i] = (ctr1 + ctr2) / 2
        else
            c1 = cellpairs[1,i]
            e1 = cellpairs[2,i]
            functions[i] = [
            Shape(c1, abs(e1), +1.0)]
            positions[i] = cartesian(center(chart(mesh, Cells[c1])))
        end
    end

    geo = mesh
    RTBasis(geo, functions, positions)
end

"""
    raviartthomas(mesh)

Conducts pre-processing on the input `mesh` by extracting the cell edges, cell pairs
    and indices required to construct the RT basis on the `mesh`.

Calls raviartthomas(mesh::Mesh, cellpairs::Array{Int,2}), which constructs
    the RT basis on the `mesh`, using the cell pairs identified.

Returns the RT basis object.
"""
function raviartthomas(mesh)
    edges = skeleton(mesh, 1)
    cps = cellpairs(mesh, edges, dropjunctionpair=true)
    ids = findall(x -> x>0, cps[2,:])
    raviartthomas(mesh, cps[:,ids])
end

"""
    raviartthomas(Γ, γ)

Constructs the RT space relative to boundary `γ` of an open surface, only
    selecting cell pairs whose common edge does not lie on `γ` . (This prevents
    the calculation of physically-impossible surface currents, such as those
    flowing 'off the edge' of a surface.)

Calls raviartthomas(Γ, duals) which constructs
    the RT basis on the `mesh`, using the cell pairs identified.
"""
function raviartthomas(Γ, γ)

    in_interior = interior_tpredicate(Γ)
    on_junction = overlap_gpredicate(γ)

    pred = c -> (in_interior(c) || on_junction(chart(Γ,c)))

    edges = skeleton(pred, Γ, 1)
    cps = cellpairs(Γ, edges, dropjunctionpair=true)

    raviartthomas(Γ, cps)
end

raowiltonglisson = raviartthomas

"""
    portcells(Γ::Mesh, γ::Mesh)
returns an array containing cell pairs of mesh Γ
around a boundary edge that overlaps with mesh γ
"""
function portcells(Γ, γ)
  in_interior = interior_tpredicate(Γ)

  overlaps = overlap_gpredicate(γ)
  on_junction = c -> overlaps(simplex(vertices(Γ,c)))

  pred = x -> (!in_interior(x) && on_junction(x)) #check only for exterior overlapping edges
  # - if γ is defined to overlap within the structure Γ, this isn't considered a port -
  edges = skeleton(pred, Γ, 1) #Take only exterior edges overlapping γ segment
  cps = cellpairs(Γ, edges, dropjunctionpair=true)
  return cps
end

"""
    rt_cedge(cps::Array{Int,2}, weight)
Computes single basis function with equally distributed constant current
leaving or entering port defined by cellpairs cps.  weight defines the
total current over the port and its direction (+ve = out, -ve = in)
"""
function rt_cedge(cps::Array{Int,2}, weight)
  numpairs = size(cps,2)
  @assert numpairs > 0
  weight = weight / numpairs #total current leaving and entering equal 1
  functions =  Vector{Shape{Float64}}(undef,numpairs) #note: not a Vector{Vector}
    for i in 1:numpairs
        c1 = cps[1,i]
        e1 = cps[2,i]
        functions[i] =
        Shape(c1, abs(e1), weight)
        # this assumes cellpairs will be adjacent eachother through the loop
      end
  functions
end

"""
    rt_vedge(cps::Array{Int,2}, weight)
Computes n-1 basis function with oscillating current in(leaving and entering)
pairs of half triangles defined over port specified by cellpairs cps.
weight defines the magnitude of individual current in and out the half triangles,
and it's polarity simply defines whether to start with in or out
"""
function rt_vedge(cps::Array{Int,2}, weight)
  numpairs = size(cps,2)
  @assert numpairs > 0
  #adjacent cells are considered a pair, so we have one less numpairs
  functions =  Vector{Vector{Shape{Float64}}}(undef,numpairs - 1)
    for i in 1:numpairs
      if i < numpairs # stop when on last cellpair
        c1 = cps[1,i]; e1 = cps[2,i]
        c2 = cps[1,i+1]; e2 = cps[2,i+1]
        functions[i] = [
          Shape(c1, abs(e1), weight),
          Shape(c2, abs(e2), -weight)]
        # this assumes cellpairs will be adjacent eachother through the loop
      end
    end
 functions
end

"""
    rt_ports(Γ::Mesh, γ::Mesh ...)
Constructs the RT space on `Γ`, relative to boundary pairs in `γ`. `γ` expects
any number of pairs-of-ports as arguments and accepts tuples, arrays,
vectors etc. e.g `rt_ports(Γ, a, b ...);` where a = [γ₁ γ₂], b = (γ₃,γ₄) etc.
`rt_ports` with no pair of ports supplied i.e `rt_ports(Γ)` reduces to the
`raviartthomas(Γ)` function. The RT space ensures current continuity in each
pair of ports. i.e. current leaving mesh Γ through γ₁ is accounted for in γ₂.

Returns the RT basis object.
"""
function rt_ports(Γ, γ...)
    #TODO: can include an extra argument for polarity applied to ports

    T = coordtype(Γ)

    internal = raviartthomas(Γ)

    fns = internal.fns
    pos = internal.pos

    for i = 1:length(γ) #loop through pairs of ports

        port1 = portcells(Γ, γ[i][1])
        port2 = portcells(Γ, γ[i][2])

        ce1 = rt_cedge(port1, +1.0)
        ce2 = rt_cedge(port2, -1.0)

        ffs = Vector{Shape{T}}(undef,length(ce1) + length(ce2))
        ffs = [ce1;ce2]

        ve1 = rt_vedge(port1, +1.0)
        ve2 = rt_vedge(port2, -1.0)

        fns = [fns;[ffs];ve1;ve2]
    end

    RTBasis(Γ, fns)

end

"""
    getindex_rtg(RT::RTBasis)
Returns the indices of the global half RWGs present in `RT`.
`RT` is typically gotten from `rt_ports`
"""
function getindex_rtg(RT::RTBasis)
    #TODO: use more reasonable condition to extract indices
    # or have rt_ports return indices with RTBasis
    idx=[]
    A = RT.fns
    for i in 1:length(A)
        if size(A[i])[1] > 2
            push!(idx,i)
        end
    end
    idx
end

divergence(X::RTBasis, geo, fns) = LagrangeBasis{0,-1,1}(geo, fns, deepcopy(positions(X)))
ntrace(X::RTBasis, geo, fns) = LagrangeBasis{0,-1,1}(geo, fns, deepcopy(positions(X)))


text:
JumpType = Union{VariableRateJump, ConstantRateJump, MassActionJump}

"""
$(TYPEDEF)

A system of jump processes.

# Fields
$(FIELDS)

# Example

```julia
using ModelingToolkit

@parameters β γ t
@variables S I R
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = ConstantRateJump(rate₂,affect₂)
j₃      = MassActionJump(2*β+γ, [R => 1], [S => 1, R => -1])
js      = JumpSystem([j₁,j₂,j₃], t, [S,I,R], [β,γ])
```
"""
struct JumpSystem{U <: ArrayPartition} <: AbstractSystem
    """
    The jumps of the system. Allowable types are `ConstantRateJump`,
    `VariableRateJump`, `MassActionJump`.
    """
    eqs::U
    """The independent variable, usually time."""
    iv::Any
    """The dependent variables, representing the state of the system."""
    states::Vector
    """The parameters of the system."""
    ps::Vector
    pins::Vector
    observed::Vector{Equation}
    """The name of the system."""
    name::Symbol
    """The internal systems."""
    systems::Vector{JumpSystem}
end

function JumpSystem(eqs, iv, states, ps;
                    pins = [],
                    observed = Equation[],
                    systems = JumpSystem[],
                    name = gensym(:JumpSystem))

    ap = ArrayPartition(MassActionJump[], ConstantRateJump[], VariableRateJump[])
    for eq in eqs
        if eq isa MassActionJump
            push!(ap.x[1], eq)
        elseif eq isa ConstantRateJump
            push!(ap.x[2], eq)
        elseif eq isa VariableRateJump
            push!(ap.x[3], eq)
        else
            error("JumpSystem equations must contain MassActionJumps, ConstantRateJumps, or VariableRateJumps.")
        end
    end

    JumpSystem{typeof(ap)}(ap, value(iv), value.(states), value.(ps), pins, observed, name, systems)
end

function generate_rate_function(js, rate)
    build_function(rate, states(js), parameters(js),
                   independent_variable(js),
                   conv = states_to_sym(states(js)),
                   expression=Val{true})
end

function generate_affect_function(js, affect, outputidxs)
    build_function(affect, states(js),
                   parameters(js),
                   conv = states_to_sym(states(js)),
                   independent_variable(js),
                   expression=Val{true},
                   headerfun=add_integrator_header,
                   outputidxs=outputidxs)[2]
end

function assemble_vrj(js, vrj, statetoid)
    rate   = @RuntimeGeneratedFunction(generate_rate_function(js, vrj.rate))
    outputvars = (value(affect.lhs) for affect in vrj.affect!)
    outputidxs = ((statetoid[var] for var in outputvars)...,)
    affect = @RuntimeGeneratedFunction(generate_affect_function(js, vrj.affect!, outputidxs))
    VariableRateJump(rate, affect)
end

function assemble_vrj_expr(js, vrj, statetoid)
    rate   = generate_rate_function(js, vrj.rate)
    outputvars = (value(affect.lhs) for affect in vrj.affect!)
    outputidxs = ((statetoid[var] for var in outputvars)...,)
    affect = generate_affect_function(js, vrj.affect!, outputidxs)
    quote
        rate = $rate
        affect = $affect
        VariableRateJump(rate, affect)
    end
end

function assemble_crj(js, crj, statetoid)
    rate   = @RuntimeGeneratedFunction(generate_rate_function(js, crj.rate))
    outputvars = (value(affect.lhs) for affect in crj.affect!)
    outputidxs = ((statetoid[var] for var in outputvars)...,)
    affect = @RuntimeGeneratedFunction(generate_affect_function(js, crj.affect!, outputidxs))
    ConstantRateJump(rate, affect)
end

function assemble_crj_expr(js, crj, statetoid)
    rate   = generate_rate_function(js, crj.rate)
    outputvars = (value(affect.lhs) for affect in crj.affect!)
    outputidxs = ((statetoid[var] for var in outputvars)...,)
    affect = generate_affect_function(js, crj.affect!, outputidxs)
    quote
        rate = $rate
        affect = $affect
        ConstantRateJump(rate, affect)
    end
end

function numericrstoich(mtrs::Vector{Pair{V,W}}, statetoid) where {V,W}
    rs = Vector{Pair{Int,W}}()
    for (spec,stoich) in mtrs
        if !(spec isa Term) && _iszero(spec)
            push!(rs, 0 => stoich)
        else
            push!(rs, statetoid[value(spec)] => stoich)
        end
    end
    sort!(rs)
    rs
end

function numericnstoich(mtrs::Vector{Pair{V,W}}, statetoid) where {V,W}
    ns = Vector{Pair{Int,W}}()
    for (spec,stoich) in mtrs
        !(spec isa Term) && _iszero(spec) && error("Net stoichiometry can not have a species labelled 0.")
        push!(ns, statetoid[spec] => stoich)
    end
    sort!(ns)
end

# assemble a numeric MassActionJump from a MT MassActionJump representing one rx.
function assemble_maj(maj::MassActionJump, statetoid, subber, invttype)
    rval = subber(maj.scaled_rates)
    rs   = numericrstoich(maj.reactant_stoch, statetoid)
    ns   = numericnstoich(maj.net_stoch, statetoid)
    maj  = MassActionJump(convert(invttype, value(rval)), rs, ns, scale_rates = false)
    maj
end

# For MassActionJumps that contain many reactions
# function assemble_maj(maj::MassActionJump{U,V,W}, statetoid, subber,
#                       invttype) where {U <: AbstractVector,V,W}
#     rval = [convert(invttype,numericrate(sr, subber)) for sr in maj.scaled_rates]
#     rs   = [numericrstoich(rs, statetoid) for rs in maj.reactant_stoch]
#     ns   = [numericnstoich(ns, statetoid) for ns in maj.net_stoch]
#     maj  = MassActionJump(rval, rs, ns, scale_rates = false)
#     maj
# end
"""
```julia
function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan,
                                    parammap=DiffEqBase.NullParameters; kwargs...)
```

Generates a blank DiscreteProblem for a pure jump JumpSystem to utilize as
its `prob.prob`. This is used in the case where there are no ODEs
and no SDEs associated with the system.

Continuing the example from the [`JumpSystem`](@ref) definition:
```julia
using DiffEqBase, DiffEqJump
u₀map = [S => 999, I => 1, R => 0]
parammap = [β => .1/1000, γ => .01]
tspan = (0.0, 250.0)
dprob = DiscreteProblem(js, u₀map, tspan, parammap)
```
"""
function DiffEqBase.DiscreteProblem(sys::JumpSystem, u0map, tspan::Tuple,
                                    parammap=DiffEqBase.NullParameters(); kwargs...)

    (u0map isa AbstractVector) || error("For DiscreteProblems u0map must be an AbstractVector.")
    u0d = Dict( value(u[1]) => u[2] for u in u0map)
    u0 = [u0d[u] for u in states(sys)]
    if parammap != DiffEqBase.NullParameters()
        (parammap isa AbstractVector) || error("For DiscreteProblems parammap must be an AbstractVector.")
        pd  = Dict( value(u[1]) => u[2] for u in parammap)
        p  = [pd[u] for u in parameters(sys)]
    else
        p = parammap
    end
    f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
    sts = states(sys)
    df = DiscreteFunction{true,true}(f, syms=tosymbol.(sts, states=sts, escape=false))
    DiscreteProblem(df, u0, tspan, p; kwargs...)
end

"""
```julia
function DiffEqBase.DiscreteProblemExpr(sys::JumpSystem, u0map, tspan,
                                    parammap=DiffEqBase.NullParameters; kwargs...)
```

Generates a black DiscreteProblem for a JumpSystem to utilize as its
solving `prob.prob`. This is used in the case where there are no ODEs
and no SDEs associated with the system.

Continuing the example from the [`JumpSystem`](@ref) definition:
```julia
using DiffEqBase, DiffEqJump
u₀map = [S => 999, I => 1, R => 0]
parammap = [β => .1/1000, γ => .01]
tspan = (0.0, 250.0)
dprob = DiscreteProblem(js, u₀map, tspan, parammap)
```
"""
function DiscreteProblemExpr(sys::JumpSystem, u0map, tspan::Tuple,
                                    parammap=DiffEqBase.NullParameters(); kwargs...)
    u0 = varmap_to_vars(u0map, states(sys))
    p  = varmap_to_vars(parammap, parameters(sys))
    # identity function to make syms works
    sts = states(sys)
    quote
        f  = DiffEqBase.DISCRETE_INPLACE_DEFAULT
        u0 = $u0
        p = $p
        tspan = $tspan
        df = DiscreteFunction{true,true}(f, syms=$(tosymbol.(sts, states=sts, escape=false)))
        DiscreteProblem(df, u0, tspan, p; kwargs...)
    end
end

"""
```julia
function DiffEqBase.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)
```

Generates a JumpProblem from a JumpSystem.

Continuing the example from the [`DiscreteProblem`](@ref) definition:
```julia
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())
```
"""
function DiffEqJump.JumpProblem(js::JumpSystem, prob, aggregator; kwargs...)

    statetoid = Dict(value(state) => i for (i,state) in enumerate(states(js)))
    eqs       = equations(js)
    invttype  = typeof(1 / prob.tspan[2])

    # handling parameter substition and empty param vecs
    p = (prob.p == DiffEqBase.NullParameters()) ? Num[] : prob.p
    parammap  = map((x,y)->Pair(x,y), parameters(js), p)
    subber    = substituter(parammap)

    majs = MassActionJump[assemble_maj(j, statetoid, subber, invttype) for j in eqs.x[1]]
    crjs = ConstantRateJump[assemble_crj(js, j, statetoid) for j in eqs.x[2]]
    vrjs = VariableRateJump[assemble_vrj(js, j, statetoid) for j in eqs.x[3]]
    ((prob isa DiscreteProblem) && !isempty(vrjs)) && error("Use continuous problems such as an ODEProblem or a SDEProblem with VariableRateJumps")
    jset = JumpSet(Tuple(vrjs), Tuple(crjs), nothing, isempty(majs) ? nothing : majs)

    if needs_vartojumps_map(aggregator) || needs_depgraph(aggregator)
        jdeps = asgraph(js)
        vdeps = variable_dependencies(js)
        vtoj = jdeps.badjlist
        jtov = vdeps.badjlist
        jtoj = needs_depgraph(aggregator) ? eqeq_dependencies(jdeps, vdeps).fadjlist : nothing
    else
        vtoj = nothing; jtov = nothing; jtoj = nothing
    end

    JumpProblem(prob, aggregator, jset; dep_graph=jtoj, vartojumps_map=vtoj, jumptovars_map=jtov, kwargs...)
end

### Functions to determine which states a jump depends on
function get_variables!(dep, jump::Union{ConstantRateJump,VariableRateJump}, variables)
    (jump.rate isa Symbolic) && get_variables!(dep, jump.rate, variables)
    dep
end

function get_variables!(dep, jump::MassActionJump, variables)
    sr = value(jump.scaled_rates)
    (sr isa Symbolic) && get_variables!(dep, sr, variables)
    for varasop in jump.reactant_stoch
        any(isequal(varasop[1]), variables) && push!(dep, varasop[1])
    end
    dep
end

### Functions to determine which states are modified by a given jump
function modified_states!(mstates, jump::Union{ConstantRateJump,VariableRateJump}, sts)
    for eq in jump.affect!
        st = eq.lhs
        any(isequal(st), sts) && push!(mstates, st)
    end
end

function modified_states!(mstates, jump::MassActionJump, sts)
    for (state,stoich) in jump.net_stoch
        any(isequal(state), sts) && push!(mstates, state)
    end
end


text:
using TuringModels

Turing.setadbackend(:reversediff)

d = DataFrame(CSV.read(joinpath(@__DIR__, "..", "..", "data", "UCBadmit.csv"),
delim=';'));

size(d) |> display # Should be 12x5

# Turing model

@model m11_5(admit, applications) = begin
    N=length(applications)
    θ ~ truncated(Exponential(1), 0, Inf)
    α ~ Normal(0,2)

    for i ∈ 1:N
        prob = logistic(α)

        # alpha and beta for the BetaBinomial must be provided.
        # The two parameterizations are related by
        # alpha = prob * theta, and beta = (1-prob) * theta.
        # See https://github.com/rmcelreath/rethinking/blob/master/man/dbetabinom.Rd
        alpha = prob * θ
        beta = (1 - prob) * θ

        admit[i] ~ BetaBinomial(applications[i], alpha, beta)
    end
end

# Sample

chns = sample(m11_5(d[:, :admit],d[:, :applications]), Turing.NUTS(0.65), 1000);

# Result rethinking

m115rethinking = "
         mean   sd  5.5% 94.5% n_eff Rhat
theta  2.74 0.96  1.43  4.37  3583    1
a       -0.37 0.31 -0.87  0.12  3210    1
";

# Show summary

describe(chns)

# End of `11/m811.5t.jl`


text:
using TailRiskScenGen: vec_isless, SurvivorApproximator, below_nonrisk_frontier, above_risk_frontier, decapitate!
using Distributions
using LinearAlgebra: I
using DataStructures

A = [[1.0, 2.0, 3.0] [2.0, 0.0, 6.0] [3.0, 8.0, 9.0]]
@test vec_isless([1.0, 2.0, 3.0], [2.0, 3.0, 4.0])
@test !vec_isless([1.0, 1.0], [1.0, 2.0])
@test !vec_isless(view(A,:,1), view(A,:,2))
@test vec_isless(view(A,:,1), view(A,:,3))
@test !vec_isless(view(A,:,3), [10.0, 9.0, 7.0])

d = 3
dist = MvNormal(Array{Float64}(I, d, d))
sample = rand(dist, 10000)
α = 0.99
surv = SurvivorApproximator(sample, α)
approx_prob, err = surv(zeros(d))
@test 0.5^d ≈ approx_prob atol=err

frontier=list([2.0, 0.0, 1.0], [1.0, 1.0, 1.0], [1.0, 0.0, 2.0])
@test above_risk_frontier([3.0, 3.0, 3.0], frontier)
@test above_risk_frontier([1.5,1.5,1.5], frontier)
@test !above_risk_frontier([3.0, 0.0, 3.0], frontier)
@test !above_risk_frontier([-2.0, -1.0, -3.0], frontier)
@test below_nonrisk_frontier([0.0, 0.0, 0.0], frontier)
@test below_nonrisk_frontier([0.5, -0.1, 1.5], frontier)
@test !below_nonrisk_frontier([1.0, 1.0, 1.0], frontier)
@test !below_nonrisk_frontier([5.0, 5.0, 1.0], frontier)

l = list(1,2,3,4)
l2 = tail(l)
decapitate!(l2)
@test l == list(1,3,4)
@test l2 == list(3,4)


text:
using Test
using Logging
import Dates
import TimeSeries
import UUIDs
import JSON3
using DataStructures: SortedDict
using DataFrames
using Random

import InfrastructureSystems
const IS = InfrastructureSystems
const BASE_DIR =
    abspath(joinpath(dirname(Base.find_package("InfrastructureSystems")), ".."))
const DATA_DIR = joinpath(BASE_DIR, "test", "data")
const FORECASTS_DIR = joinpath(DATA_DIR, "time_series")

const LOG_FILE = "infrastructure-systems.log"

LOG_LEVELS = Dict(
    "Debug" => Logging.Debug,
    "Info" => Logging.Info,
    "Warn" => Logging.Warn,
    "Error" => Logging.Error,
)

include("common.jl")

"""
Copied @includetests from https://github.com/ssfrr/TestSetExtensions.jl.
Ideally, we could import and use TestSetExtensions.  Its functionality was broken by changes
in Julia v0.7.  Refer to https://github.com/ssfrr/TestSetExtensions.jl/pull/7.
"""

"""
Includes the given test files, given as a list without their ".jl" extensions.
If none are given it will scan the directory of the calling file and include all
the julia files.
"""
macro includetests(testarg...)
    if length(testarg) == 0
        tests = []
    elseif length(testarg) == 1
        tests = testarg[1]
    else
        error("@includetests takes zero or one argument")
    end

    quote
        tests = $tests
        rootfile = @__FILE__
        if length(tests) == 0
            tests = readdir(dirname(rootfile))
            tests = filter(
                f ->
                    startswith(f, "test_") && endswith(f, ".jl") && f != basename(rootfile),
                tests,
            )
        else
            tests = map(f -> string(f, ".jl"), tests)
        end
        println()
        for test in tests
            print(splitext(test)[1], ": ")
            include(test)
            println()
        end
    end
end

function get_logging_level(env_name::String, default)
    level = get(ENV, env_name, default)
    log_level = get(LOG_LEVELS, level, nothing)
    if log_level === nothing
        error("Invalid log level $level: Supported levels: $(values(LOG_LEVELS))")
    end

    return log_level
end

function run_tests()
    console_level = get_logging_level("SYS_CONSOLE_LOG_LEVEL", "Error")
    console_logger = ConsoleLogger(stderr, console_level)
    file_level = get_logging_level("SYS_LOG_LEVEL", "Info")

    IS.open_file_logger(LOG_FILE, file_level) do file_logger
        levels = (Logging.Info, Logging.Warn, Logging.Error)
        multi_logger =
            IS.MultiLogger([console_logger, file_logger], IS.LogEventTracker(levels))
        global_logger(multi_logger)

        @time @testset "Begin Systems tests" begin
            @includetests ARGS
        end

        @test length(IS.get_log_events(multi_logger.tracker, Logging.Error)) == 0

        @info IS.report_log_summary(multi_logger)
    end
end

logger = global_logger()

try
    run_tests()
finally
    # Guarantee that the global logger is reset.
    global_logger(logger)
    nothing
end


text:
module TestIO

using Test, DataFrames, CategoricalArrays, Dates
using LaTeXStrings

# Test LaTeX export
@testset "LaTeX export" begin
    df = DataFrame(A = 1:4,
                B = ["\$10.0", "M&F", "A~B", "\\alpha"],
                C = [L"\alpha", L"\beta", L"\gamma", L"\sum_{i=1}^n \delta_i"],
                D = [1.0, 2.0, missing, 3.0],
                E = CategoricalArray(["a", missing, "c", "d"]),
                F = Vector{String}(undef, 4)
                )
    str = """
        \\begin{tabular}{r|cccccc}
        \t& A & B & C & D & E & F\\\\
        \t\\hline
        \t& $(Int) & String & LaTeXStr… & Float64⍰ & Categorical…⍰ & String\\\\
        \t\\hline
        \t1 & 1 & \\\$10.0 & \$\\alpha\$ & 1.0 & a & \\#undef \\\\
        \t2 & 2 & M\\&F & \$\\beta\$ & 2.0 &  & \\#undef \\\\
        \t3 & 3 & A\\textasciitilde{}B & \$\\gamma\$ &  & c & \\#undef \\\\
        \t4 & 4 & \\textbackslash{}\\textbackslash{}alpha & \$\\sum_{i=1}^n \\delta_i\$ & 3.0 & d & \\#undef \\\\
        \\end{tabular}
        """
    @test repr(MIME("text/latex"), df) == str
end

@testset "Huge LaTeX export" begin
    df = DataFrame(a=1:1000)
    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => false)
    show(ioc, "text/latex", df)
    @test length(String(take!(ioc.io))) > 10000

    io = IOBuffer()
    show(io, "text/latex", df)
    @test length(String(take!(io))) < 10000
end

#Test HTML output for IJulia and similar
@testset "HTML output" begin
    df = DataFrame(Fish = ["Suzy", "Amir"], Mass = [1.5, missing])
    io = IOBuffer()
    show(io, "text/html", df)
    str = String(take!(io))
    @test str == "

" *
                "	Fish	Mass
" *
                "
	String	Float64⍰

" *
                "
2 rows × 2 columns
" *
                "
1	Suzy	1.5
" *
                "
2	Amir	missing

"

    df = DataFrame(Fish = Vector{String}(undef, 2), Mass = [1.5, missing])
    io = IOBuffer()
    show(io, "text/html", df)
    str = String(take!(io))
    @test str == "

" *
                "	Fish	Mass
" *
                "
	String	Float64⍰

" *
                "
2 rows × 2 columns
" *
                "
1	#undef	1.5
" *
                "
2	#undef	missing

"

    io = IOBuffer()
    show(io, MIME"text/html"(), df, summary=false)
    str = String(take!(io))
    @test str == "

" *
                "	Fish	Mass
" *
                "
	String	Float64⍰

" *
                "
1	#undef	1.5
" *
                "
2	#undef	missing

"
end

# test limit attribute of IOContext is used
@testset "limit attribute" begin
    df = DataFrame(a=1:1000)
    ioc = IOContext(IOBuffer(), :displaysize => (10, 10), :limit => false)
    show(ioc, "text/html", df)
    @test length(String(take!(ioc.io))) > 10000

    io = IOBuffer()
    show(io, "text/html", df)
    @test length(String(take!(io))) < 10000
end

@testset "printtable" begin
    df = DataFrame(A = 1:3,
                B = 'a':'c',
                C = ["A", "B", "C"],
                D = CategoricalArray(string.('a':'c')),
                E = CategoricalArray(["A", "B", missing]),
                F = Vector{Union{Int, Missing}}(1:3),
                G = missings(3),
                H = fill(missing, 3))

    @test sprint(DataFrames.printtable, df) ==
        """
        "A","B","C","D","E","F","G","H"
        1,"'a'","A","a","A","1",missing,missing
        2,"'b'","B","b","B","2",missing,missing
        3,"'c'","C","c",missing,"3",missing,missing
        """
end

@testset "DataStreams" begin
    using DataStreams
    I = DataFrame(id = Int64[1, 2, 3, 4, 5],
        firstname = Union{String, Missing}["Benjamin", "Wayne", "Sean", "Charles", missing],
        lastname = String["Chavez", "Burke", "Richards", "Long", "Rose"],
        salary = Union{Float64, Missing}[missing, 46134.1, 45046.2, 30555.6, 88894.1],
        rate = Float64[39.44, 33.8, 15.64, 17.67, 34.6],
        hired = Union{Date, Missing}[Date("2011-07-07"), Date("2016-02-19"), missing,
                                     Date("2002-01-05"), Date("2008-05-15")],
        fired = DateTime[DateTime("2016-04-07T14:07:00"), DateTime("2015-03-19T15:01:00"),
                         DateTime("2006-11-18T05:07:00"), DateTime("2002-07-18T06:24:00"),
                         DateTime("2007-09-29T12:09:00")],
        reserved = missings(5)
    )
    sink = DataStreams.Data.close!(DataStreams.Data.stream!(I, deepcopy(I)))
    sch = DataStreams.Data.schema(sink)
    @test size(sch) == (5, 8)
    @test DataStreams.Data.header(sch) == ["id", "firstname", "lastname", "salary",
                                           "rate", "hired", "fired", "reserved"]
    @test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String,
                                          Union{Float64, Missing}, Float64,
                                          Union{Date, Missing}, DateTime, Missing)
    @test sink[:id] == [1:5;]

    transforms = Dict(1=>x->x+1)
    sink = DataStreams.Data.close!(DataStreams.Data.stream!(I, deepcopy(I);
                                                            append=true, transforms=transforms))
    sch = DataStreams.Data.schema(sink)
    @test size(sch) == (10, 8)
    @test DataStreams.Data.header(sch) == ["id", "firstname", "lastname", "salary",
                                           "rate", "hired", "fired", "reserved"]
    @test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String,
                                          Union{Float64, Missing}, Float64,
                                          Union{Date, Missing}, DateTime, Missing)
    @test sink[:id] == [1:5; 2:6]

    sink = DataStreams.Data.close!(Data.stream!(I, DataFrame, deepcopy(I)))
    sch = DataStreams.Data.schema(sink)
    @test size(sch) == (5, 8)
    @test DataStreams.Data.header(sch) == ["id", "firstname", "lastname", "salary",
                                           "rate", "hired", "fired", "reserved"]
    @test DataStreams.Data.types(sch) == (Int64, Union{String, Missing}, String,
                                          Union{Float64, Missing}, Float64,
                                          Union{Date, Missing}, DateTime, Missing)
    @test sink[:id] == [1:5;]

    # test DataFrameStream creation
    dfs = DataFrame(sch)
    DataStreams.Data.close!(dfs)
end

@testset "csv/tsv output" begin
    df = DataFrame(a = [1,2], b = [1.0, 2.0])
    @test sprint(show, "text/csv", df) == """
    "a","b"
    1,1.0
    2,2.0
    """
    @test sprint(show, "text/tab-separated-values", df) == """
    "a"\t"b"
    1\t1.0
    2\t2.0
    """
end

end # module


text:
### Device Functions

export  getDevice,
        setDevice,
        getNumDevices

function getDevice()
    n = Base.Ref{Cint}(0)
    af_get_device(n)
    Int(n[])
end

function setDevice(n::Integer)
    af_set_device(Cint(n))
end

function getNumDevices()
    n = Base.Ref{Cint}(0)
    af_get_device_count(n)
    Int(n[])
end


text:
include("util.jl")
import Cytof5.VB.ADVI

function logprob_lam(s::StateArray, y::Vector{Matrix{Float64}}, c::Cytof5.VB.Constants)
  # Get params
  eps = s.eps
  W = s.W
  sig = sqrt.(s.sig2)
  eta0 = s.eta0
  eta1 = s.eta1

  # reshape mu[z] to have size (1 x 1 x L[z])
  mu0 = reshape(-cumsum(s.delta0), 1, 1, c.L[0])
  mu1 = reshape(cumsum(s.delta1), 1, 1, c.L[1])

  # Compute Z
  H = s.H
  v = s.v
  if c.use_stickbreak
    Z = Int.(reshape(cumprod(v), 1, c.K) .> H)
  else
    Z = Int.(reshape(v, 1, c.K) .> H)
  end

  noisy_sd = sqrt(c.noisy_var)

  # Compute loglike
  function engine(i::Int)
    Ni = size(y[i], 1)

    # Ni x J x Lz
    yi = reshape(y[i], Ni, c.J, 1)

    # Ni x J x 1
    logmix_L0 = ADVI.lpdf_gmm(yi, mu0, sig[i], eta0[i:i, :, :], dims=3, dropdim=false)
    logmix_L1 = ADVI.lpdf_gmm(yi, mu1, sig[i], eta1[i:i, :, :], dims=3, dropdim=false)

    # Ni x J x K -> Ni x K
    Z_rs = reshape(Z, 1, c.J, c.K)
    Z_mix = ADVI.sumdd(Z_rs .* logmix_L1 + (1 .- Z_rs) .* logmix_L0, dims=2)

    # Ni x K
    lpi_quite = (Z_mix .+ log.(W[i:i, :])) .+ log1p(-s.eps[i])

    # Ni x 1
    lpi_noisy = ADVI.sumdd(ADVI.lpdf_normal.(y[i], 0., noisy_sd), dims=2) .+ log(s.eps[i])

    return [lpi_quite lpi_noisy]
  end

  return [engine(i) for i in 1:c.I]
end


text:
dir = @__DIR__
push!(LOAD_PATH, dir * "/../pomdp")
push!(LOAD_PATH, dir * "/../adv")
using adv
using pomdp

mutable struct PofsNode
  x::AbstractArray{Float64, 1}
  u::AbstractArray{Float64, 1}
  b::AbstractArray{Float64, 1}
  adv_x::AbstractArray{Float64, 1}
  c::pomdp.Collision
  r::Float64

  PofsNode() = new(Float64[], Float64[], Float64[], Float64[],
                   pomdp.Collision(), 0.0)
end

const olm_pofs_dt = 0.5
function plan_pofs(x::AbstractArray{Float64, 1}, b::AbstractArray{Float64, 1},
                   agent::Agent, adv_agent::Agent, world::World,
                   reward::Function, ctrl_d::Discretization, depth::Int)
  # backup variables
  adv_agent_controller! = adv_agent.controller!
  adv_agent_custom = adv_agent.custom

  # collision from the perspective of the adversary!!!
  value = PofsNode()
  value.x = x
  # uniform prior belief
  value.b = b
  value.adv_x = adv_agent.x

  adv_agent.controller! = pomdp.adv_controller!
  (c, _) = adv.predict_collision(adv_agent.x, agent.x, world)
  value.c = c

  root = GroupTree(value)

  (rs, us) = select_action_pofs(root, agent, adv_agent, world, reward, ctrl_d,
                                depth)
  ret = root

  # revert variables
  adv_agent.controller! = adv_agent_controller!
  adv_agent.custom = adv_agent_custom

  if us == Float64[]
    error("Action empty, error in POFS")
  end

  return (us, ret)
end

function select_action_pofs(node::GroupTree, agent::Agent, adv_agent::Agent,
                            world::World, reward::Function,
                            ctrl_d::Discretization, depth::Int)
  if depth <= 0
    return (0.0, Float64[])
  end

  node.value.x[1] = mod(node.value.x[1], world.road.path.S[end])
  (adv_c, _) = predict_collision(node.value.adv_x, node.value.x, world)

  # allocate next level of the tree
  la_len = ctrl_d.thr[end] * ctrl_d.pt[end] # number of actions to survey
  node.next = Array{Array{GroupTree, 1}, 1}(undef, la_len)
  for aidx in 1:la_len
    node.next[aidx] = Array{GroupTree, 1}(undef, length(pomdp.ACTIONS))
  end

  # predict adv agent doing three possible actions
  adv_NX = Array{Array{Float64, 1}, 1}(undef, length(pomdp.ACTIONS))
  for o in pomdp.ACTIONS
    oidx = Int(o)

    adv_nx = copy(node.value.adv_x)
    #adv_u = fill(0.0, ctrl_d.dim)
    adv_agent.custom[3] = o
    sim.advance!(sim.default_dynamics!, adv_nx, Pair(adv_agent, world), 0.0,
                 olm_pofs_dt, olm_h)
    adv_nx[1] = mod(adv_nx[1], world.road.path.S[end])

    adv_NX[oidx] = adv_nx
  end

  # for each action, evaluate
  max_r = -Inf
  us = Float64[]
  for la in 0:(la_len - 1)
    u = dis.ls2x(ctrl_d, la)
    agent.custom = u
    nx = copy(node.value.x)
    sim.advance!(sim.default_dynamics!, nx, Pair(agent, world), 0.0,
                 olm_pofs_dt, olm_h)
    nx[1] = mod(nx[1], world.road.path.S[end])

    ra = 0.0
    for o in pomdp.ACTIONS
      oidx = Int(o)
      adv_nx = adv_NX[oidx]
      (nc, _) = adv.predict_collision(nx, adv_nx, world)

      custom = agent.custom
      agent.custom = nc
      r = reward(node.value.x, u, nx, agent, world)
      agent.custom = custom

      bp = pomdp.update_belief(node.value.b, o, adv_c)

      value = PofsNode()
      value.x = nx
      value.u = u
      value.b = bp
      value.adv_x = adv_nx
      value.c = nc

      next_node = GroupTree(value)
      (next_r, _) = select_action_pofs(next_node, agent, adv_agent, world,
                                       reward, ctrl_d, depth - 1)
      r += olm_gamma * next_r
      next_node.value.r = r
      node.next[la + 1][oidx] = next_node

      ra += r * reduce(+, map(i -> node.value.b[i] *
                              pomdp.P_adv(o, adv_c, pomdp.DRIVERS[i]),
                              1:length(pomdp.DRIVERS)))
    end

    if ra > max_r
      max_r = ra
      us = u
    end
  end

  return (max_r, us)
end

function controller_pofs!(u::AbstractArray{Float64},
                          x::AbstractArray{Float64},
                          dx::AbstractArray{Float64},
                          agent_world::Pair{Agent, World}, t::Float64)
  agent = agent_world.first
  u[1] = agent.custom[1]
  u[2] = agent.custom[2]
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Graphics (Modia3D/Graphics/_module.jl)
#
# Content:
#   This file contains mutable visual geometries that are only used for animation
#   (immutable solid geometries that have a volume are defined in module Modia3D.Solids.
#    They can be additionally used for collision handling and other purposes).

mutable struct Sphere <: Modia3D.AbstractGeometry
   Dx::Float64
   material::Material

   function Sphere(Dx::Float64; material=Material())
      @assert(Dx >= 0.0)
      new(Dx, material)
   end
end

mutable struct Ellipsoid <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Dz::Float64
  material::Material

  function Ellipsoid(Dx::Float64,Dy::Float64,Dz::Float64; material=Material())
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Dz >= 0.0)
    new(Dx,Dy,Dz,material)
  end
end

mutable struct Box <: Modia3D.AbstractGeometry
  Lx::Float64
  Ly::Float64
  Lz::Float64
  material::Material

  function Box(Lx::Float64,Ly::Float64,Lz::Float64; material=Material())
    @assert(Lx >= 0.0)
    @assert(Ly >= 0.0)
    @assert(Lz >= 0.0)
    new(Lx,Ly,Lz,material)
  end
end

mutable struct Cylinder <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  material::Material

  function Cylinder(Dx::Float64, Lz::Float64; Dy=Dx, material=Material())
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    new(Dx,Dy,Lz,material)
  end
end

mutable struct Capsule <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  material::Material

  function Capsule(Dx::Float64, Lz::Float64; Dy=Dx, material=Material())
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    new(Dx,Dy,Lz,material)
  end
end

mutable struct Beam <: Modia3D.AbstractGeometry
  Lx::Float64
  Dy::Float64
  Lz::Float64
  material::Material

  function Beam(Lx::Float64,Dy::Float64,Lz::Float64; material=Material())
    @assert(Lx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    new(Lx,Dy,Lz,material)
  end
end

mutable struct Cone <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  material::Material

  function Cone(Dx::Float64, Lz::Float64; Dy=Dx,relativeTipDiameter=0.5,relativeInnerDiameter=1, material=Material())
    # Tip diameter >= 0.0, Tip diameter is relative to base diameter (0 = real cone, 1 = cylinder),
    # Inner diameter any real number, Inner diameter of the cone is proportional to tip size (0..1)
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(relativeTipDiameter>=0.0 && relativeTipDiameter<=1.0)
    @assert(relativeInnerDiameter>=0.0)
    extras=[relativeTipDiameter,relativeInnerDiameter,0.0]
    new(Dx,Dy,Lz,extras,material)
  end
end

mutable struct Pipe <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  material::Material

  function Pipe(Dx::Float64, Lz::Float64; Dy=Dx, relativeInnerDiameter=0.9,relativeTipDiameter=1.0, material=Material())
    #relativeInnerDiameter ... Inner diameter of pipe (0 = rigid, 1 = fully hollow), every real number is allowed
    #relativeTipDiameter ... Tip diameter of pipe relative to base diameter (0 = real cone, 1 = cylinder), every real number >= 0 is possible
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(relativeTipDiameter>=0.0)
    extras=[0.0,relativeInnerDiameter,relativeTipDiameter]
    new(Dx,Dy,Lz,extras,material)
  end
end

mutable struct Spring <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  material::Material

  function Spring(Dx::Float64, Lz::Float64; Dy=Dx,windings=5,springRadius=0.02, material=Material())
    # windings ... Winding number of the spring, integer > 0
    # springRadius ... Radius of spring wire > 0.0
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(isinteger(windings))
    @assert(windings>0.0)
    @assert(springRadius>0.0)
    extras=[windings,springRadius,0.0]
    new(Dx,Dy,Lz,extras,material)
  end
end

mutable struct GearWheel <: Modia3D.AbstractGeometry
  Dx::Float64
  Dy::Float64
  Lz::Float64
  extras::MVector{3,Float64}
  material::Material

  function GearWheel(Dx::Float64, Lz::Float64; Dy=Dx,relativeInnerDiameter=1,teeth=20,angle=0, material=Material())
    # relativeInnerDiameter ... Ratio of inner/outer radius of gearwheel, proportional to radius of gearwheel (0..1)
    # teeth ... Number of gearwheel teeth (negative teeth number for inner gearwheel (not 0)
    # angle ... Angle for beveled gearwheels (0 = cylindrical gearwheel)
    @assert(Dx >= 0.0)
    @assert(Dy >= 0.0)
    @assert(Lz >= 0.0)
    @assert(isinteger(teeth))
    @assert(teeth!=0.0)
    extras=[relativeInnerDiameter,teeth,angle]
    new(Dx,Dy,Lz,extras,material)
  end
end

mutable struct CoordinateSystem <: Modia3D.AbstractVisualElement
  Lx::Float64

  function CoordinateSystem(Lx::Float64)
    @assert(Lx >= 0.0)
    new(Lx)
  end
end

mutable struct Grid <: Modia3D.AbstractVisualElement
  Lx::Float64
  Ly::Float64
  extras::MVector{3,Float64}

  function Grid(Lx::Float64, Ly::Float64; distance=0.05, lineWidth=1.0)
    # distance ... Distance or grid points in x- and in y-direction
    # lineWidth ... Line width of the grid
    @assert(Lx >= 0.0)
    @assert(Ly >= 0.0)
    @assert(distance > 0.0)
    @assert(lineWidth > 0.0)
    extras=[distance,lineWidth,0.0]
    new(Lx,Ly,extras)
  end
end

mutable struct FileMesh <: Modia3D.AbstractGeometry
   filename::AbstractString
   scaleFactor::MVector{3,Float64}
   useMaterialColor::Bool
   smoothNormals::Bool
   material::Material

   function FileMesh(filename::AbstractString; scaleFactor=Basics.onesMVector()  , useMaterialColor::Bool=false, smoothNormals::Bool=false, material=Material())
      @assert(isfile(filename))
      @assert(scaleFactor[1] >= 0.0)
      @assert(scaleFactor[2] >= 0.0)
      @assert(scaleFactor[3] >= 0.0)
      new(filename, scaleFactor, useMaterialColor, smoothNormals, material)
   end
end
FileMesh(filename::AbstractString, scale::Number; useMaterialColor::Bool=false, smoothNormals::Bool=false, material=Material()) =
         FileMesh(filename; scaleFactor=[scale,scale,scale], useMaterialColor=useMaterialColor, smoothNormals=smoothNormals, material=material)


text:
export remove_gapped_sequences, remove_close_seqs, remove_gapped_cols, afa2fasta

##############################################################
"""
    function bool_gaps(seq, max_gaps)

    (fastapath, kwargs..) --> outpath

    USE:
    Given an amino acid sequence and a number of gaps,
    returns "true" if the number of gaps in the sequence is bigger than threshold.

 	INPUT:
    "seq": amino acid sequence in numbers or letters format

    "max_gaps": real number between 0 and the sequence length

    OUTPUT:

    boolean

"""

#--------------------------------------------------------

bool_gaps(seq::String, max_gaps::Real) = sum([1 for l in seq if l == '-']) > max_gaps ? true : false
bool_gaps(seq::Array{<:Integer, 1}, max_gaps::Real) = sum([1 for l in seq if l == 21 ]) > max_gaps ? true : false

##############################################################
"""
    function my_hamming(v1, v2)

    (v1, v2) --> hamming distance

    USE:
    Function used to compute the hamming distance between two vectors.
    Works for strings as well.

    INPUT:

    OUTPUT:

    number between 0 and the length of the two vectors

"""

#--------------------------------------------------------

my_hamming(v1, v2) = sum( [v1[i] .!= el_v2 for (i, el_v2) in enumerate(v2)] )

##############################################################
"""
    function remove_gapped_sequences(fastapath::AbstractString
        ;outpath::AbstractString = "",threshold::Real = 0.1)

    (fastapath, kwargs..) --> outpath

    USE:
    Given the path of an aligned MSA,
    produces a new MSA that does not contain gapped sequence
    (up to the fraction [of the alignment] specified by "threshold").

    INPUT:
    "fastapah": path of the MSA

    "outpath": path of the resulting MSA, if not specified adds "_maxXXgaps"
    to the filename (before .file)

    "threshold": fraction of gaps allowed in the MSA,
    an integer number is also accepted for the total number of gaps

    OUTPUT:

"""

#--------------------------------------------------------

function remove_gapped_sequences(fastapath::AbstractString; outpath::AbstractString = "" ,threshold::Real = 0.1)
    N = 0
    f = FastaReader(fastapath)
    for (desc, seq_string) in f
        N = length(seq_string)
        break
    end

    frac_gaps = 0
    max_gaps = 0
    if threshold <= 1
        frac_gaps = Int64(round(threshold*100, digits = 0))
        max_gaps =  Int64(round(threshold*N, digits = 0))

    else
        frac_gaps = Int64(round((threshold/N)*100, digits = 0))
        max_gaps = threshold

    end

    if outpath == ""
        dir, file = splitdir(fastapath)
        split_file = split(file, ".")
        l_file = length(split_file)
        split_file[end-1] = split_file[end-1]*"_max0$(frac_gaps)rowgaps"
        outpath = joinpath(dir, join(split_file, "."))
    end

    count_removed_seqs = 0
    for (desc, seq_string) in f
        if !bool_gaps(seq_string, max_gaps)
            writefasta(outpath, [(desc, seq_string)], "a")
        else
            count_removed_seqs += 1
        end
    end
    close(f)

    plural_flag = "s"
    verb_flag = "ve"
    if count_removed_seqs == 1
        plural_flag = ""
        verb_flag = "s"
    end

    println("$(count_removed_seqs) sequence$(plural_flag) with more than $(max_gaps) gaps ha$(verb_flag) been removed.")
    return outpath
end

##############################################################
"""
    function remove_close_seqs(fastapath::AbstractString, wtpaths...; outpath::AbstractString = "",
    threshold::Real = 0.8)

    (fastapath, kwargs..) --> outpath

    USE:
    Given the path of an aligned MSA,
    produces a new MSA that does not contain sequences too close
    (in terms of Hamming distance) to the wildtype(s)
    (up to the fraction of proteins lenght specified by "threshold").

    INPUT:
    "fastapah": path of the MSA

    "outpath": path of the resulting MSA, if not specified adds "_noclose"
    to the filename (before .file)

    "threshold": max allowed fractional sequence identity to the wts
    OUTPUT:

"""

#--------------------------------------------------------

function remove_close_seqs(fastapath::AbstractString, wtpaths...; outpath::AbstractString = "",
    threshold::Real = 0.8)
    N = 0
    n_wts = length(wtpaths)
    f = FastaReader(fastapath)
    for (desc, seq_string) in f
        N = length(seq_string)
        break
    end

    if threshold <= 1
        seqID =  Int64(round( threshold*N, digits = 0))
        frac_close = Int64(round(threshold*100, digits = 0))
    end

    if outpath == ""
        dir, file = splitdir(fastapath)
        split_file = split(file, ".")
        l_file = length(split_file)
        split_file[end-1] = split_file[end-1]*"_max0$(frac_close)wtid"
        outpath = joinpath(dir, join(split_file, "."))
    end

    wts = [wt[:, 1] for wt in fasta2matrix.(wtpaths)]
    count_removed_seqs = 0
    for (desc, seq_string) in f
        if prod(  N .- my_hamming.( [string2vec(seq_string) for i in 1:n_wts] , wts) .< [seqID for i in 1:n_wts] )
            writefasta(outpath, [(desc, seq_string)], "a")
        else
            count_removed_seqs += 1
        end
    end

    close(f)

    plural_flag = "s"
    verb_flag = "ve"
    if count_removed_seqs == 1
        plural_flag = ""
        verb_flag = "s"
    end
    plural_flag_wt = ""
    if n_wts > 1
        plural_flag_wt = "s"
    end

    println("$(count_removed_seqs) sequence$(plural_flag) with less than $(N - seqID) mutations wrt to the wildtype$(plural_flag_wt) ha$(verb_flag) been removed.")
    return outpath
end

##############################################################
"""
    function remove_gapped_cols(fastapath::AbstractString; outpath::AbstractString = "",
        threshold::Real = 0.8, only_flanks::Bool = false)

    (fastapath, kwargs..) --> outpath

    USE:
    Given the path of an aligned MSA,
    produces a new MSA that does not contain gapped columns
    (up to the fraction specified by "threshold".)
    "threshold" = xx means that columns are removed if they contain
    a fraction of gaps > than xx.

    INPUT:
    "fastapah": path of the MSA

    "outpath": path of the resulting MSA, if not specified adds "_fracXXcolgaps"
    to the filename (before .file)

    "threshold": fraction of gaps allowed in the MSA columns

    OUTPUT:
    outpupath and columns removed
"""

#--------------------------------------------------------

function remove_gapped_cols(fastapath::AbstractString; outpath::AbstractString = "", threshold::Real = 0.8, only_flanks::Bool = false)
    MSA = fasta2matrix(fastapath)
    M, N = size(MSA)

    col_to_rem = [col for col in 1:N if bool_gaps(MSA[:, col], threshold*M)]
    col_to_keep = [i for i in 1:N if i ∉ col_to_rem]

    if only_flanks == true
        col_to_rem = [col for col in 1:N if bool_gaps(MSA[:, col], threshold*M) && (N - col)%Int64(round(N*9/10, digits = 0)) < N/10    ]
        col_to_keep = [i for i in 1:N if i ∉ col_to_rem]
    end

    max_gaps =  Int64(round(threshold*M, digits = 0))
    frac_gaps = Int64(round(threshold*100, digits = 0))

    if outpath == ""
        dir, file = splitdir(fastapath)
        split_file = split(file, ".")
        l_file = length(split_file)

        if only_flanks == false
            split_file[end-1] = split_file[end-1]*"_max0$(frac_gaps)colgaps"
        else
            split_file[end-1] = split_file[end-1]*"_max0$(frac_gaps)flankcolgaps"
        end

        outpath = joinpath(dir, join(split_file, "."))
    end
    count_removed_cols = length(col_to_rem)

    f = FastaReader(fastapath)
    for (desc, seq_string) in f
        writefasta(outpath, [(desc, seq_string[col_to_keep])], "a")
    end

    close(f)
    plural_flag = "s"
    verb_flag = "ve"
    if count_removed_cols == 1
        plural_flag = ""
        verb_flag = "s"
    end

    println("$(count_removed_cols) column$(plural_flag) with more than $(max_gaps) gaps ha$(verb_flag) been removed.")
    return outpath,col_to_rem
end

##############################################################
"""
    function afa2fasta(path_in::String, path_out::String = "")

    (path_in, path_out) -->

    USE:
    Given the path of a MSA in .afa format, that is with insertions and gaps,
    writes in "path_out" the same MSA by removing insertions.

    INPUT:
    "path_in": path of starting .afa file

    "path_out": path of cleaned (removed gaps) MSA

    OUTPUT:

"""

#--------------------------------------------------------

function afa2fasta(path_in::String, path_out::String = "")
    path_out == "" && (path_out = join( split(path_in, '.')[1:end-1], "."  )*"_removed_dels.fasta")
    f = FastaReader(path_in)
    for (desc, seq_string) in f
        seq_clean = remove_dels(seq_string)
        writefasta(path_out, [(desc, seq_clean)], "a")
    end
    close(f)
end

##############################################################
"""
    function remove_dels(string_wdel::String)

    (string_wdel) --> string without lowecase letters

    USE:
    Given a sequence with lowercase letters (insertions) and '.'
    returns the same string without them.

    INPUT:
    "string_wdel": string with insertions

    OUTPUT:

"""

#--------------------------------------------------------

function remove_dels(string_wdel::String)
    clean_string = ""
    for lettr in string_wdel
        if islowercase(lettr)
            nothing
        elseif (lettr == '.')
            nothing
        else
            clean_string = string(clean_string, lettr)
        end
    end
    return clean_string
end

max_stars_repo_name: matteobisardi/KitMSA

max_issues_repo_name: matteobisardi/KitMSA

max_forks_repo_name: matteobisardi/KitMSA


text:
module UnitfulAtomic

import Unitful
using Unitful: @unit, Dimension, Dimensions, NoDims, NoUnits, Units, dimension, uconvert, ustrip

export auconvert, aunit, austrip

# The following five constants are used as the “base” atomic units
@unit me_au "mₑ"  ElectronRestMass      Unitful.me                     false
@unit e_au  "e"   ElementaryCharge      Unitful.q                      false
@unit ħ_au  "ħ"   ReducedPlanckConstant Unitful.ħ                      false
@unit k_au  "k"   BoltzmannConstant     Unitful.k                      false
@unit a0_au "a₀"  BohrRadius            5.291_772_109_03e-11*Unitful.m false # CODATA 2018 recommended value

# Hartree energy is derived from the base atomic units
@unit Eh_au "Eₕ"  HartreeEnergy         1ħ_au^2/(me_au*a0_au^2)       false

# Units that are not Hartree atomic units, but are commonly used in atomic physics
@unit Ry    "Ry"  RydbergEnergy         Eh_au//2                      false
@unit μ_N   "μ_N" NuclearMagneton       e_au*ħ_au/(2*Unitful.mp)      false

# Aliases for units
const bohr    = a0_au
const hartree = Eh_au

"""
    aunit(x::Unitful.Quantity)
    aunit(x::Unitful.Units)
    aunit(x::Unitful.Dimensions)

Returns the appropriate atomic unit (a `Unitful.Units` object) for the dimension of `x`.

# Examples

```jldoctest
julia> aunit(2.3u"cm")
a₀

julia> aunit(u"T")
a₀^-2 e^-1 ħ
```
"""
aunit(x) = aunit(dimension(x))

# `aunit` for `Dimension` types
aunit(x::Dimension{:Length})      = (a0_au)^x.power
aunit(x::Dimension{:Mass})        = (me_au)^x.power
aunit(x::Dimension{:Time})        = (ħ_au/Eh_au)^x.power
aunit(x::Dimension{:Current})     = (e_au*Eh_au/ħ_au)^x.power
aunit(x::Dimension{:Temperature}) = (Eh_au/k_au)^x.power

# For dimensions not specified above, there is no atomic unit.
aunit(::Dimension{D}) where D = throw(ArgumentError("no atomic unit defined for dimension $D."))

# `aunit` for `Dimensions` types
@generated aunit(::Dimensions{N}) where N = prod(aunit, N)
aunit(::typeof(NoDims)) = NoUnits

# Simplifications for some derived dimensions, so that e.g. `aunit(u"J")` returns `Eₕ`
# instead of `a₀^2 mₑ Eₕ^2 ħ^-2`. The following units/dimensions are considered:
#   * Energy: Eₕ
#   * Momentum: ħ/a₀
#   * Action/angular momentum: ħ
#   * Force: Eₕ/a₀
#   * Pressure: Eₕ/a₀^3
#   * E-field: Eₕ/(e*a₀)
#   * B-field: ħ/(e*a₀^2)
#   * Voltage/electric potential: Eₕ/e
#   * Magnetic dipole moment: e*ħ/mₑ
#   * Entropy: k
#   * 1st hyperpolarizability: e^3*a₀^3/Eₕ^2
#   * 2nd hyperpolarizability: e^4*a₀^4/Eₕ^3
#   * E-field gradient: Eₕ/(e*a₀^2)
#   * Electric polarizability: e^2*a₀^2/Eₕ
#   * Permittivity: e^2/(a₀*Eₕ)
for unit in (:(Eh_au), :(ħ_au/a0_au), :(ħ_au), :(Eh_au/a0_au), :(Eh_au/a0_au^3),
             :(Eh_au/(e_au*a0_au)), :(ħ_au/(e_au*a0_au^2)), :(Eh_au/e_au), :(e_au*ħ_au/me_au),
             :(k_au), :(e_au^3*a0_au^3/Eh_au^2), :(e_au^4*a0_au^4/Eh_au^3),
             :(Eh_au/(e_au*a0_au^2)), :(e_au^2*a0_au^2/Eh_au), :(e_au^2/(a0_au*Eh_au)))
    @eval aunit(::typeof(dimension($unit))) = $unit
end

"""
    auconvert(x::Unitful.Quantity)

Convert a quantity to the appropriate atomic unit.

# Examples

```jldoctest
julia> auconvert(13.6u"eV")
0.4997907858599377 Eₕ

julia> auconvert(20u"nm")
377.94522509156565 a₀
```
"""
auconvert(x) = uconvert(aunit(x), x)

"""
    auconvert(u::Unitful.Units, x::Number)

Interpret `x` as a quantity given in atomic units and convert it to the unit `u`.

# Examples

```jldoctest
julia> auconvert(u"eV", 1)  # convert 1 Eₕ to eV
27.211386013449417 eV

julia> auconvert(u"m", 1)   # convert 1 a₀ to m
5.2917721067e-11 m
```
"""
auconvert(u::Units, x::Number) = uconvert(u, x*aunit(u))

"""
    austrip(x::Unitful.Quantity)

Returns the value of the quantity converted to atomic units as a number type (i.e., with the
units removed). This is equivalent to `Unitful.ustrip(auconvert(x))`.

# Examples

```jldoctest
julia> austrip(13.6u"eV")
0.4997907858599377

julia> austrip(20u"nm")
377.94522509156565
```
"""
austrip(x) = ustrip(auconvert(x))

# In order to enable precompilation, some things need to be set at runtime
const localunits = Unitful.basefactors
function __init__()
    merge!(Unitful.basefactors, localunits)
    Unitful.register(UnitfulAtomic)
end

end # module UnitfulAtomic


text:
using ImageCore
import ImageView

imgdy = normedview(UInt8[y for y=0x00:0xff, x=0:100])
imgdx = normedview(UInt8[x for y=1:100, x=0x00:0xff])
imshow_now(imgdy)
imshow_now(imgdx)

A = zeros(300, 200)
lbl = similar(A, Int)
fill!(lbl, 0)
A[80:120, 50:150] .= 0.8
lbl[80:120, 50:150] .= 1
cx, cy = 200, 100
for j = -30:30, i = -30:30
    if i^2+j^2 <= 900
        A[cx+i,cy+j] = rand()
        lbl[cx+i,cy+j] = 2
    end
end
ImageView.imshowlabeled(A, lbl; name="circle and rectangle")


text:
"""
    OneClassJKNN( Normal, New; J::Int = 1, K::Int = 1, DistanceType = "euclidean" )

Creates a one class JK-NN classifier from `Normal` data and evaluates it on `New` data. This compares the inter sample
distance (`DistanceType`) between a `New` and `Normal` J nearest neighbors to the K nearest neighbors of those
J nearest neighbors in the `Normal` set. No cut off is provided, that should be done by the end-user. A typical
cut off value is 1.0 .
"""
function OneClassJKNN( Normal, New; J::Int = 1, K::Int = 1, DistanceType = "euclidean" )
    #ToDo: Only find distances of K's that matter! This will speed things up and slim memory use.
    Obs = size( Normal )[ 1 ]
    ObsNew = size( New )[ 1 ]
    DistMat = zeros( Obs, ObsNew )
    JNNPool = zeros( ObsNew, J )
    IntraDataMeanDist = zeros(ObsNew)
    #Apply Distance Fn
    if DistanceType == "euclidean"
        DistMat = SquareEuclideanDistance( Normal, New )
    elseif DistanceType == "manhattan"
        DistMat = ManhattanDistance( Normal, New )
    end
    #Find J nearest neighbors between new samples and the old samples
    for obs in 1 : ObsNew
        JNNPool[ obs , : ] = sortperm( DistMat[ : , obs ] )[ 1 : J ]
        IntraDataMeanDist[obs] = sum( DistMat[ JNNPool[ obs , : ] .|> Int, obs] )
    end
    #Find K nearest neighbors in the original/normal data
    InterDataDists = zeros( Obs, Obs )
    if DistanceType == "euclidean"
        InterDataDists = SquareEuclideanDistance( Normal )
    elseif DistanceType == "manhattan"
        InterDataDists = ManhattanDistance( Normal )
    end
    #For each new sample, find it's J nearest neighbors K nearest neighbor distances
    MeanDists = zeros( ObsNew )
    for obs in 1 : ObsNew
        #Find row where J lives ( ObsNew, J )
        Smpls = JNNPool[ obs , : ] .|> Int
        for j in 1 : J
            KNNInds = sortperm( InterDataDists[  : , Smpls[ j ] ] )
            Screened = (KNNInds[1] == Smpls[j]) ? KNNInds[ 2 : ( K + 1 ) ] : KNNInds[ 1 : K ]
            MeanDists[ obs ] += sum( InterDataDists[ Smpls, Screened ] )
        end
    end

    MeanDists ./= J #* K #No need  to divide by K
    return (IntraDataMeanDist ./ MeanDists)
end


text:
# tests for the various kinds of waits
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

MPI.Init()

myrank = MPI.Comm_rank(MPI.COMM_WORLD)
commsize = MPI.Comm_rank(MPI.COMM_WORLD)

nsends = 2
send_arr = [ArrayType{Int}([i]) for i = 1:nsends]
recv_arr = [ArrayType{Int}(undef,1) for i = 1:nsends]

# send to self
send_reqs = [MPI.Isend(send_arr[i], MPI.COMM_WORLD; dest=myrank, tag=i) for i = 1:nsends]
recv_reqs = [MPI.Irecv!(recv_arr[i], MPI.COMM_WORLD; source=myrank, tag=i) for i = 1:nsends]

send_check = zeros(Int, nsends)
recv_check = zeros(Int, nsends)

for i = 1:nsends
    idx = MPI.Waitany(send_reqs)
    @test send_reqs[idx] == MPI.REQUEST_NULL
    send_check[idx] += 1
end
@test send_check == ones(Int, nsends)

for i = 1:nsends
    idx = MPI.Waitany(recv_reqs)
    @test recv_reqs[idx] == MPI.REQUEST_NULL
    recv_check[idx] += 1
end
@test recv_check == ones(Int, nsends)

MPI.Barrier(MPI.COMM_WORLD)
MPI.Finalize()
@test MPI.Finalized()


text:
using Statistics
using Oceananigans.BuoyancyModels: g_Earth
using Oceananigans.Architectures: device_event

using Oceananigans.Models.HydrostaticFreeSurfaceModels:
    ImplicitFreeSurface,
    FreeSurface,
    FFTImplicitFreeSurfaceSolver,
    PCGImplicitFreeSurfaceSolver,
    MatrixImplicitFreeSurfaceSolver,
    implicit_free_surface_step!,
    implicit_free_surface_linear_operation!

function set_simple_divergent_velocity!(model)
    # Create a divergent velocity
    grid = model.grid

    u, v, w = model.velocities
    η = model.free_surface.η

    u .= 0
    v .= 0
    η .= 0

    imid = Int(floor(grid.Nx / 2)) + 1
    jmid = Int(floor(grid.Ny / 2)) + 1
    CUDA.@allowscalar u[imid, jmid, 1] = 1

    update_state!(model)

    return nothing
end

function run_pcg_implicit_free_surface_solver_tests(arch, grid)

    Δt = 900
    Nx = grid.Nx
    Ny = grid.Ny

    # Create a model
    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        momentum_advection = nothing,
                                        free_surface = ImplicitFreeSurface(solver_method=:PreconditionedConjugateGradient,
                                                                           tolerance = 1e-15))

    set_simple_divergent_velocity!(model)
    implicit_free_surface_step!(model.free_surface, model, Δt, 1.5, device_event(arch))

    η = model.free_surface.η
    @info "PCG implicit free surface solver test, norm(η_pcg): $(norm(η)), maximum(abs, η_pcg): $(maximum(abs, η))"

    # Extract right hand side "truth"
    right_hand_side = model.free_surface.implicit_step_solver.right_hand_side

    # Compute left hand side "solution"
    g = g_Earth
    η = model.free_surface.η
    ∫ᶻ_Axᶠᶜᶜ = model.free_surface.implicit_step_solver.vertically_integrated_lateral_areas.xᶠᶜᶜ
    ∫ᶻ_Ayᶜᶠᶜ = model.free_surface.implicit_step_solver.vertically_integrated_lateral_areas.yᶜᶠᶜ

    left_hand_side = Field{Center, Center, Nothing}(grid)
    implicit_free_surface_linear_operation!(left_hand_side, η, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

    # Compare
    extrema_tolerance = 1e-9
    std_tolerance = 1e-9

    CUDA.@allowscalar begin
        @test maximum(abs, interior(left_hand_side) .- interior(right_hand_side)) < extrema_tolerance
        @test std(interior(left_hand_side) .- interior(right_hand_side)) < std_tolerance
    end

    return nothing
end

function run_matrix_implicit_free_surface_solver_tests(arch, grid)

    Δt = 900
    Nx = grid.Nx
    Ny = grid.Ny

    # Create a model
    model = HydrostaticFreeSurfaceModel(grid = grid,
                                        momentum_advection = nothing,
                                        free_surface = ImplicitFreeSurface(solver_method=:HeptadiagonalIterativeSolver,
                                                                           tolerance = 1e-15))

    set_simple_divergent_velocity!(model)
    implicit_free_surface_step!(model.free_surface, model, Δt, 1.5, device_event(arch))

    η = model.free_surface.η
    @info "Matrix implicit free surface solver test, norm(η_mat): $(norm(η)), maximum(abs, η_mat): $(maximum(abs, η))"

    # Extract right hand side "truth"
    right_hand_side = model.free_surface.implicit_step_solver.right_hand_side

    # Compute left hand side "solution"
    g = g_Earth
    η = model.free_surface.η
    ∫ᶻ_Axᶠᶜᶜ = model.free_surface.implicit_step_solver.vertically_integrated_lateral_areas.xᶠᶜᶜ
    ∫ᶻ_Ayᶜᶠᶜ = model.free_surface.implicit_step_solver.vertically_integrated_lateral_areas.yᶜᶠᶜ

    left_hand_side = Field{Center, Center, Nothing}(grid)
    implicit_free_surface_linear_operation!(left_hand_side, η, ∫ᶻ_Axᶠᶜᶜ, ∫ᶻ_Ayᶜᶠᶜ, g, Δt)

    # Compare
    extrema_tolerance = 1e-9
    std_tolerance = 1e-9

    CUDA.@allowscalar begin
        # Note: `right_hand_side` is not a Field but an Array/CuArray
        @test maximum(abs, interior(left_hand_side) .- reshape(right_hand_side, (Nx, Ny, 1))) < extrema_tolerance
        @test std(interior(left_hand_side) .- reshape(right_hand_side, (Nx, Ny, 1))) < std_tolerance
    end

    return nothing
end

@testset "Implicit free surface solver tests" begin
    for arch in archs

        rectilinear_grid = RectilinearGrid(arch, size = (128, 1, 5),
                                                 x = (0, 1000kilometers), y = (0, 1), z = (-400, 0),
                                                 topology = (Bounded, Periodic, Bounded))

        lat_lon_grid = LatitudeLongitudeGrid(arch, size = (90, 90, 5),
                                                   longitude = (-30, 30), latitude = (15, 75), z = (-4000, 0))

        for grid in (rectilinear_grid, lat_lon_grid)
            @info "Testing PreconditionedConjugateGradient implicit free surface solver [$(typeof(arch)), $(typeof(grid).name.wrapper)]..."
            run_pcg_implicit_free_surface_solver_tests(arch, grid)

            @info "Testing Matrix implicit free surface solver [$(typeof(arch)), $(typeof(grid).name.wrapper)]..."
            run_matrix_implicit_free_surface_solver_tests(arch, grid)
        end

        @info "Testing implicit free surface solvers compared to FFT [$(typeof(arch))]..."

        Δt = 900

        mat_free_surface = ImplicitFreeSurface(solver_method = :HeptadiagonalIterativeSolver,    tolerance=1e-15, maximum_iterations=128^3)
        pcg_free_surface = ImplicitFreeSurface(solver_method = :PreconditionedConjugateGradient, tolerance=1e-15, maximum_iterations=128^3)
        fft_free_surface = ImplicitFreeSurface(solver_method = :FastFourierTransform)

        pcg_model = HydrostaticFreeSurfaceModel(grid = rectilinear_grid,
                                                momentum_advection = nothing,
                                                free_surface = pcg_free_surface)

        fft_model = HydrostaticFreeSurfaceModel(grid = rectilinear_grid,
                                                momentum_advection = nothing,
                                                free_surface = fft_free_surface)

        mat_model = HydrostaticFreeSurfaceModel(grid = rectilinear_grid,
                                                momentum_advection = nothing,
                                                free_surface = mat_free_surface)

        @test fft_model.free_surface.implicit_step_solver isa FFTImplicitFreeSurfaceSolver
        @test pcg_model.free_surface.implicit_step_solver isa PCGImplicitFreeSurfaceSolver
        @test mat_model.free_surface.implicit_step_solver isa MatrixImplicitFreeSurfaceSolver

        for m in (mat_model, pcg_model, fft_model)
            set_simple_divergent_velocity!(m)
            implicit_free_surface_step!(m.free_surface, m, Δt, 1.5, device_event(arch))
        end

        mat_η = mat_model.free_surface.η
        pcg_η = pcg_model.free_surface.η
        fft_η = fft_model.free_surface.η

        mat_η_cpu = Array(interior(mat_η))
        pcg_η_cpu = Array(interior(pcg_η))
        fft_η_cpu = Array(interior(fft_η))

        @info "FFT/PCG/MAT implicit free surface solver comparison, " *
            "maximum(abs, η_mat - η_fft): $(maximum(abs, mat_η_cpu .- fft_η_cpu)), " *
            "maximum(abs, η_pcg - η_fft): $(maximum(abs, pcg_η_cpu .- fft_η_cpu)), " *
            "maximum(abs, η_mat): $(maximum(abs, mat_η_cpu)), " *
            "maximum(abs, η_pcg): $(maximum(abs, pcg_η_cpu)), " *
            "maximum(abs, η_fft): $(maximum(abs, fft_η_cpu)), "

        @test all(mat_η_cpu .≈ fft_η_cpu)
        if arch isa CPU
            @test all(pcg_η_cpu .≈ fft_η_cpu)
        else
            # It seems that the PCG algorithm is not always stable on sverdrup's GPU, often leading to failure.
            # This behavior is not observed on tartarus, where this test _would_ pass.
            # Suffice to say that the FFT solver appears to be accurate (as of this writing), and tests pass
            # on the CPU.
            @info "  Skipping comparison between pcg and fft implicit free surface solver"
            @test_skip all(pcg_η_cpu .≈ fft_η_cpu)
        end

        pcg_η = pcg_model.free_surface.η
        fft_η = fft_model.free_surface.η

        Δη_mat = Array(interior(mat_η) .- interior(fft_η))
        Δη_pcg = Array(interior(pcg_η) .- interior(fft_η))

        @test all(isapprox.(Δη_mat, 0, atol=sqrt(eps(eltype(rectilinear_grid)))))
        @test all(isapprox.(Δη_pcg, 0, atol=sqrt(eps(eltype(rectilinear_grid)))))
    end
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Miscellaneous functions related to the Text backend.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _update_column_width(column_width::Int, largest_cell_width::Int, column_width_specification::Int, maximum_column_width::Int, minimum_column_width::Int)

Compute the column width `column_width` considering the largest cell width in
the column `largest_cell_width`, the user specification in
`column_width_specification`, and the maximum and minimum allowed column width
in `maximum_column_width` and `minimum_column_width`, respectively.
"""
function _update_column_width(
    column_width::Int,
    largest_cell_width::Int,
    column_width_specification::Int,
    maximum_column_width::Int,
    minimum_column_width::Int
)
    if column_width_specification ≤ 0
        # The columns width must never be lower than 1.
        column_width = max(column_width, largest_cell_width)

        # Make sure that the maximum column width is respected.
        if (maximum_column_width > 0) && (maximum_column_width < column_width)
            column_width = maximum_column_width
        end

        # Make sure that the minimum column width is respected.
        if (minimum_column_width > 0) && (minimum_column_width > column_width)
            column_width = minimum_column_width
        end
    else
        column_width = column_width_specification
    end

    return column_width
end


text:
# `InputModel` subtypes are types of models for the dendritic input
abstract type InputModel end

# default InputModels are time-homogeneous
update!(x, model::InputModel, dt, t) = update!(x, model, dt)

struct PoissonExpModel{T1, T2, T3} <: InputModel
    ρ::T1    # input firing rate
    τm::T2   # membrane time constant
    ε0::T3   # amplitude of EPSP
    dim::Int # dimensionality
end

# default amplitude is chosen such that the postsynaptic kernel is normalized to one
PoissonExpModel(ρ, τm, dim::Int) = PoissonExpModel(ρ, τm, 1/τm, dim)

function update!(x::AbstractArray, model::PoissonExpModel, dt)
    τm = model.τm
    ε0 = model.ε0
    αm = dt / τm
    λ = dt * model.ρ

    for i in eachindex(x)
        @inbounds x[i] *= exp(-αm)
        @inbounds x[i] += ε0 * rand(Poisson(λ))
    end

    return x
end

# input neurons are active in equal-sized blocks
struct BlockPoissonExpModel{T1, T2, T3, T4} <: InputModel
    ρ::T1           # input firing rate
    τm::T2          # membrane time constant
    ε0::T3          # amplitude of EPSP
    numblocks::Int  # number of blocks
    blocksize::Int  # size of blocks
    τblock::T4      # duration of activation of a single block
end

# default amplitude is chosen such that the postsynaptic kernel is normalized to one
BlockPoissonExpModel(
    ρ, τm, numblocks::Int, blocksize::Int, τblock
) = BlockPoissonExpModel(ρ, τm, 1/τm, numblocks::Int, blocksize::Int, τblock)

function update!(x::AbstractArray, model::BlockPoissonExpModel, dt, t)
    τm = model.τm
    ε0 = model.ε0
    αm = dt / τm
    λ = dt * model.ρ

    # decay for all indices
    x .*= exp(-αm)

    # compute scheduled block id (zero-based)
    block_id = floor(Int, (t / model.τblock) % model.numblocks)

    # sample new spikes only for indices within block:
    istart = block_id * model.blocksize + 1
    istop  = (block_id + 1) * model.blocksize
    for i in istart:istop
        @inbounds x[i] += ε0 * rand(Poisson(λ))
    end

    return x
end

# `SynapseModel` subtypes are types of models for the evolution of synaptic weights
abstract type SynapseModel end

# default SynapseModels are time-homogeneous
update!(x, model::SynapseModel, dt, t) = update!(x, model, dt)

struct OUModel{T1, T2} <: SynapseModel
    τ::T1 # time constant for weight change
    σs::T2 # variance of weights
end

function update!(w::AbstractArray, model::OUModel, dt)
    α = dt / model.τ
    σ = sqrt(2 * dt * model.σs / model.τ)

    for i in eachindex(w)
        # Ornstein-Uhlenbeck process
        @inbounds w[i] += -α * w[i] + σ * randn()
    end

    return w
end

struct State{T1, T2, T3}
    w::T1 # weights
    x::T2 # dendritic input
    η::T3 # adaptation state
end

State(w, x) = State(w, x, [0.])

function State(dim::Integer)
    w = zeros(dim)
    x = zeros(dim)
    return State(w, x)
end

State(imodel::PoissonExpModel) = State(imodel.dim)
State(imodel::BlockPoissonExpModel) = State(imodel.numblocks * imodel.blocksize)

function update!(state::State, model::InputModel, dt)
    update!(state.x, model, dt)
    return state
end

function update!(state::State, model::SynapseModel, dt)
    update!(state.w, model, dt)
    return state
end

function update!(state::State, model::InputModel, dt, t)
    update!(state.x, model, dt, t)
    return state
end

function update!(state::State, model::SynapseModel, dt, t)
    update!(state.w, model, dt, t)
    return state
end

# `OutputModel` subtypes are types of models for the output neuron
abstract type OutputModel end

# default OutputModels are time-homogeneous
update!(x, model::OutputModel, dt, t) = update!(x, model, dt)

struct NeuronObs{T1, T2}
    x::T1 # input spikes / current
    y::T2 # output spikes
end

# must return a NeuronObs
function update!(state::State, model::OutputModel, dt) end

struct ExpGainModel{T1, T2} <: OutputModel
    g0::T1 # baseline firing rate
    β::T2 # gain
end

function update!(state::State, model::ExpGainModel, dt)
    x  = state.x
    w  = state.w
    g0 = model.g0
    β  = model.β

    λ = g0 * dt * exp(β * dot(w, x))
    y = rand(Poisson(λ))

    return NeuronObs(x, y)
end

struct AdaptiveExpGainModel{T1, T2, T3, T4} <: OutputModel
    g0::T1 # baseline firing rate
    β::T2 # gain
    α::T3 # log firing rate change upon spikes
    τr::T4 # adaptation time constant
end

function update!(state::State, model::AdaptiveExpGainModel, dt)
    x  = state.x
    w  = state.w
    η  = state.η # is going to be modified
    g0 = model.g0
    β  = model.β
    α  = model.α
    τr = model.τr

    λ = g0 * dt * exp(β * dot(w, x) + η[1])
    y = rand(Poisson(λ))
    η .= η .* exp(-dt / τr) .+ α * y # exponential adaptation

    return NeuronObs(x, y)
end


text:
using Packing, FreeTypeAbstraction, SignedDistanceFields, FreeType
using GLAbstraction: Texture, render, gpu_data

using GPUArrays: AbstractSampler, JLBackend
using GPUArrays.JLBackend: Sampler

type TextureAtlas{IMG <: AbstractSampler}
    rectangle_packer::RectanglePacker
    mapping         ::Dict{Any, Int} # styled glyph to index in sprite_attributes
    index           ::Int
    images          ::IMG
    attributes      ::Vector{Vec4f0}
    scale           ::Vector{Vec2f0}
    extent          ::Vector{FontExtent{Float64}}
end

function TextureAtlas(initial_size = (2048, 2048))
    images = Sampler(zeros(Float32, initial_size...))
    TextureAtlas(
        RectanglePacker(SimpleRectangle(0, 0, initial_size...)),
        Dict{Any, Int}(),
        1,
        images,
        Vec4f0[],
        Vec2f0[],
        FontExtent{Float64}[]
    )
end

begin #basically a singleton for the textureatlas
    # random list of chars we cache
    # to make runtests fast, until we figure out a better way to cache
    # newly rendered chars.
    const local _tobe_cached = [
        'π','∮','⋅','→','∞','∑','∏','∀','∈','ℝ','⌈','⌉','−','⌊','⌋','α','∧','β','∨','ℕ','⊆','₀',
        '⊂','ℤ','ℚ','ℂ','⊥','≠','≡','≤','≪','⊤','⇒','⇔','₂','⇌','Ω','⌀',
    ]
    const local _atlas_cache = TextureAtlas[]
    const local _cache_path = joinpath(dirname(@__FILE__), "..", ".cache", "texture_atlas.jls")
    const local _default_font = Vector{Ptr{FreeType.FT_FaceRec}}[]
    const local _alternative_fonts = Vector{Ptr{FreeType.FT_FaceRec}}[]

    function defaultfont()
        if isempty(_default_font)
            push!(_default_font, newface(assetpath("fonts", "DejaVuSans.ttf")))
        end
        _default_font[]
    end
    function alternativefonts()
        if isempty(_alternative_fonts)
            alternatives = [
                "DejaVuSans.ttf",
                "NotoSansCJKkr-Regular.otf",
                "NotoSansCuneiform-Regular.ttf",
                "NotoSansSymbols-Regular.ttf",
                "FiraMono-Medium.ttf"
            ]
            for font in alternatives
                push!(_alternative_fonts, newface(assetpath("fonts", font)))
            end
        end
        _alternative_fonts
    end

    function cached_load()
        if isfile(_cache_path)
            return open(_cache_path) do io
                deserialize(io)
            end
        else
            atlas = TextureAtlas()
            info("Caching fonts, this may take a while. Needed only on first run!")
            for c in '\u0000':'\u00ff' #make sure all ascii is mapped linearly
                insert_glyph!(atlas, c, defaultfont())
            end
            for c in _tobe_cached
                insert_glyph!(atlas, c, defaultfont())
            end

            to_cache(atlas) # cache it
            return atlas
        end
    end

    function to_cache(atlas)
        if !ispath(dirname(_cache_path))
            mkdir(dirname(_cache_path))
        end
        open(_cache_path, "w") do io
            serialize(io, atlas)
        end
    end

    function get_texture_atlas()
        if isempty(_atlas_cache)
            atlas = cached_load() # initialize only on demand
            push!(_atlas_cache, atlas)
            atlas
        else
            _atlas_cache[]
        end
    end

end

function glyph_index!(atlas::TextureAtlas, c::Char, font)
    if FT_Get_Char_Index(font[], c) == 0
        for afont in alternativefonts()
            if FT_Get_Char_Index(afont[], c) != 0
                font = afont
            end
        end
    end
    if c < '\u00ff' && font == defaultfont() # characters up to '\u00ff'(255), are directly mapped for default font
        Int(c)+1
    else #others must be looked up, since they're inserted when used first
        return insert_glyph!(atlas, c, font)
    end
end

glyph_scale!(c::Char, scale) = glyph_scale!(get_texture_atlas(), c, defaultfont(), scale)
glyph_uv_width!(c::Char) = glyph_uv_width!(get_texture_atlas(), c, defaultfont())

function glyph_uv_width!(atlas::TextureAtlas, c::Char, font)
    atlas.attributes[glyph_index!(atlas, c, font)]
end
function glyph_scale!(atlas::TextureAtlas, c::Char, font, scale)
    atlas.scale[glyph_index!(atlas, c, font)] .* (scale * 0.02)
end
function glyph_extent!(atlas::TextureAtlas, c::Char, font)
    atlas.extent[glyph_index!(atlas, c, font)]
end

function bearing(extent)
    Point2f0(
        extent.horizontal_bearing[1],
        -(extent.scale[2] - extent.horizontal_bearing[2])
    )
end
function glyph_bearing!(atlas::TextureAtlas, c::Char, font, scale)
    bearing(atlas.extent[glyph_index!(atlas, c, font)]) .* (scale * 0.02)
end
function glyph_advance!(atlas::TextureAtlas, c::Char, font, scale)
    atlas.extent[glyph_index!(atlas, c, font)].advance .* (scale * 0.02)
end

insert_glyph!(atlas::TextureAtlas, glyph::Char, font) = get!(atlas.mapping, (glyph, font)) do
    uv, extent, width_nopadd, pad = render(atlas, glyph, font)
    tex_size       = Vec2f0(size(atlas.images))
    uv_start       = Vec2f0(uv.x, uv.y)
    uv_width       = Vec2f0(uv.w, uv.h)
    real_start     = uv_start + pad - 1 # include padding
    # padd one additional pixel
    relative_start = real_start ./ tex_size # use normalized texture coordinates
    relative_width = (real_start + width_nopadd + 2) ./ tex_size

    uv_offset_width = Vec4f0(relative_start..., relative_width...)
    i = atlas.index
    push!(atlas.attributes, uv_offset_width)
    push!(atlas.scale, Vec2f0(width_nopadd + 2))
    push!(atlas.extent, extent)
    atlas.index = i + 1
    return i
end

function sdistancefield(img, downsample = 8, pad = 8*downsample)
    w, h = size(img)
    wpad = 0; hpad = 0;
    while w % downsample != 0
        w += 1
    end
    while h % downsample != 0
        h += 1
    end
    w, h = w + 2pad, h + 2pad #pad this, to avoid cuttoffs

    in_or_out = Matrix{Bool}(w, h)
    @inbounds for i=1:w, j=1:h
        x, y = i-pad, j-pad
        in_or_out[i,j] = checkbounds(Bool, img, x, y) && img[x,y] > 0.5*255
    end
    yres, xres = div(w, downsample), div(h, downsample)
    sd = sdf(in_or_out, xres, yres)
    map(Float16, sd)
end

function GLAbstraction.render(atlas::TextureAtlas, glyph::Char, font)
    #select_font_face(cc, font)
    if glyph == '\n' # don't render  newline
        glyph = ' '
    end
    downsample = 5
    pad = 8
    bitmap, extent = renderface(font, glyph, (50*downsample, 50*downsample))
    sd = sdistancefield(bitmap, downsample, downsample*pad)
    extent = extent ./ Vec2f0(downsample)
    rect = SimpleRectangle(0, 0, size(sd)...)
    uv = push!(atlas.rectangle_packer, rect) #find out where to place the rectangle
    uv == nothing && error("texture atlas is too small. Resizing not implemented yet. Please file an issue at GLVisualize if you encounter this") #TODO resize surface
    atlas.images[uv.area] = sd
    uv.area, extent, Vec2f0(size(bitmap)) ./ downsample, pad
end


text:
@testset "Array quantizer" begin

m = 2
k = 2
Q = QuantizedArrays.OrthogonalQuantization
for T in [UInt8, UInt16, UInt32, UInt64, Int8, Int16, Int32, Int64,
          Float16, Float32, Float64]
    for mᵢ in 1:2
        X = rand(T[0,1,2,3], m, 10)  # values in 0-3 to keep distances small
        aqm = build_quantizer(X, k=k, m=mᵢ, method=:sample)
        @test aqm isa ArrayQuantizer{Q, UInt8, Distances.SqEuclidean, T, 2}
        @test length(aqm.codebooks) == mᵢ
        qm = QuantizedArrays.quantize_data(aqm, X)
        @test qm isa Matrix{UInt8}
        @test size(qm) == (mᵢ, size(X,2))

        Xᵥ= rand(T[0,1,2,3,4], 10)
        aqv = build_quantizer(Xᵥ, k=k, m=mᵢ, method=:sample)
        @test aqv isa ArrayQuantizer{Q, UInt8, Distances.SqEuclidean, T, 1}
        @test length(aqv.codebooks) == 1
        qv = QuantizedArrays.quantize_data(aqv, Xᵥ)
        @test qv isa Matrix{UInt8}
        @test length(qv) == length(Xᵥ)
    end
end

T = Float32
Twrong = Float16
X = rand(T, m, 100)
m = 2
k = 2
aq = build_quantizer(X, k=k, m=m, method=:sample)
@test QuantizedArrays.codebooks(aq) == aq.codebooks

@test_throws AssertionError QuantizedArrays.quantize_data(aq, rand(T, m+1, 100))
@test_throws MethodError QuantizedArrays.quantize_data(aq, rand(Twrong, m, 100))

# The ArrayQuantizer works with any `m`:
X2 = rand(T, m+1, 100)
aq2 = ArrayQuantizer(X2, k=k, m=m, method=:sample)
@test aq2 isa ArrayQuantizer{Q, UInt8, Distances.SqEuclidean, T, 2}
@test length(aq2.codebooks) == m

aq3 = ArrayQuantizer(X, k=k, m=m+100, method=:sample)
@test aq3 isa ArrayQuantizer{Q, UInt8, Distances.SqEuclidean, T, 2}
@test length(aq3.codebooks) == m

end