text:
## Type definitions for intermediate representations of crystals up to a net

include("specialsolver.jl")

import Base: ==
using Serialization
using Tokenize

"""
    EquivalentPosition

Representation of a symmetry operation in 3D, defined by an affine function.
"""
struct EquivalentPosition
    mat::SMatrix{3,3,Rational{Int}, 9}
    ofs::SVector{3,Rational{Int}}
end

"""
Find the reference identifiers for the three dimensions for the CIF group called
symmetry_equiv_pos_as_xyz or space_group_symop_operation_xyz.
Usually this is simply ("x", "y", "z").
"""
function find_refid(eqs)
    isempty(eqs) && return ("x", "y", "z")
    for eq in eqs # Normally it should be the first one but I'm not sure of the specs here
        ts = collect(tokenize(eq))
        any(Tokenize.Tokens.isoperator∘Tokenize.Tokens.kind, ts) && continue
        refid = [""]
        not_at_the_end = true
        for t in ts
            @assert not_at_the_end
            k = Tokenize.Tokens.kind(t)
            if k === Tokenize.Tokens.IDENTIFIER
                @assert refid[end] == ""
                refid[end] = Tokenize.Tokens.untokenize(t)
            elseif k === Tokenize.Tokens.COMMA || k === Tokenize.Tokens.SEMICOLON
                @assert refid[end] != ""
                push!(refid, "")
            elseif k === Tokenize.Tokens.ENDMARKER
                not_at_the_end = false
            else
                error("Input string {$eq} is not a valid symmetry equivalent")
            end
        end
        if not_at_the_end
            error("Unknown end of line marker for symmetry equivalent {$eq}")
        end
        if length(refid) != 3 || refid[end] == ""
            error("Input string {$eq} is not a valid symmetry equivalent")
        end
        return tuple(refid[1], refid[2], refid[3])
    end
    return ("x", "y", "z")
end

function Base.parse(::Type{EquivalentPosition}, s::AbstractString, refid=("x", "y", "z"))
    const_dict = Dict{String, Int}(refid[1]=>1, refid[2]=>2, refid[3]=>3)
    mat = zeros(Rational{Int}, 3, 3)
    ofs = zeros(Rational{Int}, 3)
    curr_num::Union{Int, Nothing} = nothing
    curr_val::Union{Rational{Int}, Nothing} = nothing
    curr_sign::Union{Bool, Nothing} = nothing
    encountered_div::Bool = false
    i = 1
    something_written = false
    for x in tokenize(s)
        k = Tokenize.Tokens.kind(x)
        k === Tokenize.Tokens.WHITESPACE && continue
        if k === Tokenize.Tokens.INTEGER
            @assert isnothing(curr_val)
            if encountered_div
                curr_val = Rational{Int}(Int(curr_num), parse(Int, x.val))
                curr_num = nothing
                encountered_div = false
            else
                @assert isnothing(curr_num)
                curr_num = parse(Int, x.val)
            end
        else
            @assert !encountered_div
            if k === Tokenize.Tokens.IDENTIFIER
                if !isnothing(curr_num)
                    @assert isnothing(curr_val)
                    curr_val = curr_num
                    curr_num = nothing
                end
                sign = isnothing(curr_sign) ? 1 : 2*curr_sign - 1
                val = isnothing(curr_val)  ? Rational{Int}(1) : curr_val
                j = const_dict[Tokenize.Tokens.untokenize(x)]
                mat[i,j] += sign * val
                curr_val = nothing
                curr_sign = nothing
                something_written = true
            else
                if x.kind === Tokenize.Tokens.FWD_SLASH
                    @assert isnothing(curr_val)
                    @assert !isnothing(curr_num)
                    encountered_div = true
                    continue
                end
                if !isnothing(curr_num)
                    @assert isnothing(curr_val)
                    curr_val = curr_num
                    curr_num = nothing
                end
                if !isnothing(curr_val)
                    sign = isnothing(curr_sign) ? 1 : 2*curr_sign - 1
                    @ifwarn begin if !iszero(ofs[i])
                            @warn "Existing offset already existing for position $i in {$s}"
                        end
                    end
                    ofs[i] += sign * Rational{Int}(curr_val)
                    curr_val = nothing
                    curr_sign = nothing
                else
                    @assert isnothing(curr_sign)
                end
                if x.kind === Tokenize.Tokens.PLUS
                    curr_sign = true
                elseif x.kind === Tokenize.Tokens.MINUS
                    curr_sign = false
                elseif k === Tokenize.Tokens.COMMA || k === Tokenize.Tokens.SEMICOLON
                    i > 2 && error("Too many dimensions specified for symmetry equivalent {$s}")
                    something_written || error("{$s} is not a symmetry equivalent (no dependency expressed in position $i)")
                    something_written = false
                    i += 1
                else
                    k !== Tokenize.Tokens.ENDMARKER && error("Unknown end of line marker for symmetry equivalent {$s}")
                    i!= 3 && error("Input string {$s} is not a valid symmetry equivalent")
                end
            end
        end
    end
    EquivalentPosition(SMatrix{3, 3, Rational{Int}, 9}(mat), SVector{3, Rational{Int}}(ofs))
end

function Base.show(io::IO, eq::EquivalentPosition)
    function rationaltostring(x::Rational{<:Integer}, notofs::Bool, first::Bool)
        if notofs && (x == 1 || x == -1)
            return x < 0 ? '-' : first ? "" : "+"
        end
        sign = x < 0 || first ? "" : "+"
        sign * (x.den == 1 ? string(x.num) : string(x.num)*'/'*string(x.den))
    end
    xyz = ('x', 'y', 'z')
    for i in 1:3
        first = true
        for j in 1:3
            if eq.mat[i,j] != 0
                coeff = rationaltostring(eq.mat[i,j], true, first)
                first = false
                print(io, coeff)
                print(io, xyz[j])
            end
        end
        if eq.ofs[i] != 0
            print(io, rationaltostring(eq.ofs[i], false, first))
        end
        i < 3 && print(io, ',')
    end
    nothing
end

"""
    Cell

Representation of a periodic cell in 3D. Contains information about the cell
(axes lengths and angles) and its symmetry group.
"""
struct Cell
    latticesystem::Symbol
    spacegroup::String
    tablenumber::Int
    mat::SMatrix{3,3,BigFloat,9} # Cartesian coordinates of a, b and c
    equivalents::Vector{EquivalentPosition}

    function Cell(lattice, space, table, a, b, c, α, β, γ, eq)
        cosα = cosd(α); cosβ = cosd(β); cosγ = cosd(γ); sinγ = sind(γ)
        ω = sqrt(1 - cosα^2 - cosβ^2 - cosγ^2 + 2*cosα*cosβ*cosγ)
        mat = SMatrix{3,3,BigFloat,9}([a  b*cosγ  c*cosβ ;
                                      0   b*sinγ  c*(cosα - cosβ*cosγ)/sinγ ;
                                      0   0       c*ω/sinγ ])
        return new(lattice, space, table, mat, eq)
    end

    function Cell(c::Cell, eqs::Vector{EquivalentPosition})
        return new(c.latticesystem, c.spacegroup, c.tablenumber, c.mat, eqs)
    end
end

function ==(c1::Cell, c2::Cell)
    c1.latticesystem == c2.latticesystem && c1.spacegroup == c2.spacegroup &&
    c1.tablenumber == c2.tablenumber && c1.mat == c2.mat && c1.equivalents == c2.equivalents
end

Cell() = Cell(Symbol(""), "P 1", 0, 10, 10, 10, 90, 90, 90, EquivalentPosition[])
function cell_parameters(mat::StaticArray{Tuple{3,3},BigFloat})
    _a, _b, _c = eachcol(mat)
    a = norm(_a)
    b = norm(_b)
    c = norm(_c)
    α = acosd(_b'_c/(b*c))
    β = acosd(_c'_a/(c*a))
    γ = acosd(_a'_b/(a*b))
    return (a, b, c, α, β, γ)
end
cell_parameters(cell::Cell) = cell_parameters(cell.mat)
function Cell(cell::Cell, mat::StaticArray{Tuple{3,3},BigFloat})
    a, b, c, α, β, γ = cell_parameters(mat)
    return Cell(cell.latticesystem, cell.spacegroup, cell.tablenumber, a, b, c, α, β, γ, cell.equivalents)
end
function Base.show(io::IO, cell::Cell)
    a, b, c, α, β, γ = Float64.(cell_parameters(cell))
    print(io, Cell, "(\"$(cell.spacegroup)\", ($a, $b, $c), ($α, $β, $γ))")
end

"""
    CIF

Representation of a .cif file.
"""
struct CIF
    cifinfo::Dict{String, Union{String, Vector{String}}}
    cell::Cell
    ids::Vector{Int}
    types::Vector{Symbol}
    pos::Matrix{Float64}
    bonds::Matrix{Bool}
end

function keep_atoms(cif::CIF, kept)
    kept_ids = sort!([cif.ids[i] for i in kept])
    unique!(kept_ids)
    idmap = Vector{Int}(undef, length(cif.types)) # upper bound on maximum(kept_ids)
    for (i,x) in enumerate(kept_ids)
        idmap[x] = i
    end
    return CIF(cif.cifinfo, cif.cell, [idmap[cif.ids[i]] for i in kept],
               cif.types[kept_ids], cif.pos[:, kept], cif.bonds[kept, kept])
end

function periodic_distance(u, v)
    dst = 0.0
    #=@inbounds=# for i in 1:3
        x = abs2(u[i] - v[i])
        if x > 0.25
            x = (1 - sqrt(x))^2
        end
        dst += x
    end
    return sqrt(dst)
end

"""
    periodic_distance(u, v, mat)

Distance between points u and v, given as triplets of fractional coordinates, in
a repeating unit cell of matrix mat.
The distance is the shortest between all equivalents of u and v.
"""
function periodic_distance(u, v, mat)
    dst = 0.0
    x = similar(u)
    uu = copy(u)
    vv = copy(v)
    #=@inbounds=# for i in 1:3
        diff = u[i] - v[i]
        if diff > 0.5
            x[i] = diff - 1
        elseif - diff > 0.5
            x[i] = diff + 1
        else
            x[i] = diff
        end
    end
    return norm(mat*x)
end

"""
    remove_partial_occupancy(::CIF)

Only keep one atom per atom site.
"""
function remove_partial_occupancy(cif::CIF)
    points::Vector{SVector{3,Float64}} = collect(eachcol(cif.pos))
    perm = sortperm(points)
    n = length(perm)
    last_triplet = points[perm[1]] .+ (1,1,1)
    last_position = 0
    toremove = Int[]
    for i in 1:n
        j = perm[i]
        thispoint = points[j]
        if norm(thispoint .- last_triplet) < 4e-4
            push!(toremove, max(last_position, j)) # keep the first that appears
        else
            last_position = j
            last_triplet = thispoint
        end
    end
    if isempty(toremove) # Nothing to do, we simply alias the CIF to help the compiler (it may be useless)
        return CIF(cif.cifinfo, cif.cell, cif.ids, cif.types, cif.pos, cif.bonds)
    end
    @ifwarn @warn "This CIF file contains a site with multiple atoms. Only one atom will be kept per atom site."
    sort!(toremove)
    ids = deleteat!(copy(cif.ids), toremove)
    tokeep = deleteat!(collect(1:n), toremove)
    bonds = cif.bonds[tokeep, tokeep]
    pos = cif.pos[:,tokeep]
    @assert allunique(collect(eachcol(pos)))
    return CIF(cif.cifinfo, cif.cell, ids, cif.types, pos, bonds)
end

"""
    prune_collisions(::CIF)

Remove all atoms that are suspiciously close to one another. This arises when all
the possible positions of at atom are superposed in the CIF file, typically for
a solvent which should be disregarded anyway.
"""
function prune_collisions(cif::CIF)
    toremove = Int[]
    points::Vector{SVector{3,Float64}} = collect(eachcol(cif.pos))
    n = length(points)
    for i in 1:n, j in (i+1):n
        if periodic_distance(points[i], points[j], Float64.(cif.cell.mat)) < 0.55
            push!(toremove, i)
            push!(toremove, j)
        end
    end
    if isempty(toremove) # Nothing to do, we simply alias the CIF to help the compiler (it may be useless)
        return CIF(cif.cifinfo, cif.cell, cif.ids, cif.types, cif.pos, cif.bonds)
    end
    @ifwarn @warn "This CIF file contains multiple colliding atoms. All colliding atoms will be removed."
    unique!(sort!(toremove))
    ids = deleteat!(copy(cif.ids), toremove)
    tokeep = deleteat!(collect(1:n), toremove)
    bonds = cif.bonds[tokeep, tokeep]
    pos = cif.pos[:,tokeep]
    return CIF(cif.cifinfo, cif.cell, ids, cif.types, pos, bonds)
end

"""
    expand_symmetry(::CIF)

Applies all the symmetry operations listed in the CIF file to the atoms and the bonds.
"""
function expand_symmetry(c::CIF)
    cif::CIF = remove_partial_occupancy(c)
    n = length(cif.ids)
    oldbonds = [(i,j) for i in 1:n for j in i:n if cif.bonds[i,j]]
    newbonds = Set{Tuple{Int,Int}}(oldbonds)
    newids::Vector{Int} = copy(cif.ids)
    newpos::Vector{Vector{Float64}} = collect(eachcol(cif.pos))
    ret = Vector{Vector{Int}}
    #=@inbounds=# for equiv in cif.cell.equivalents
        image = zeros(Int, n)
        for i in 1:length(cif.ids)
            v = newpos[i]
            p = Vector(equiv.mat*v .+ equiv.ofs)
            @. p = p - floor(p)
            for j in 1:length(newpos)
                if periodic_distance(newpos[j], p, Float64.(cif.cell.mat)) < 0.5
                    image[i] = j
                    break
                end
            end
            if image[i] == 0
                push!(newpos, p)
                push!(newids, cif.ids[i])
                image[i] = length(newpos)
            end
        end
        for (i,j) in oldbonds
            push!(newbonds, minmax(image[i], image[j]))
        end
    end
    m = length(newids)
    bonds = zeros(Bool, m, m)
    for (i,j) in newbonds
        bonds[i,j] = true
        bonds[j,i] = true
    end

    newcif = CIF(cif.cifinfo, deepcopy(cif.cell), newids, copy(cif.types), reduce(hcat, newpos), bonds)

    return prune_collisions(newcif)
end

# function set_unique_bond_type!(cif::CIF, bond_length, bonded_atoms::Tuple{Symbol, Symbol}, onlykeep, tol)
#     @assert iszero(cif.bonds)
#     indices = [i for i in 1:length(cif.ids) if cif.types[cif.ids[i]] ∈ onlykeep]
#     #=@inbounds=# for _i in 1:length(indices)
#     i = indices[_i]
#     cif.bonds[i,i] = false
#         Threads.@threads for _j in i+1:length(indices)
#             j = indices[_j]
#             if minmax(cif.types[cif.ids[i]], cif.types[cif.ids[j]]) == bonded_atoms
#                 bonded = abs2(periodic_distance(cif.pos[:,i], cif.pos[:,j], cif.cell.mat) - bond_length) <= tol
#                 cif.bonds[i,j] = bonded
#                 cif.bonds[j,i] = bonded
#             end
#         end
#     end
#     nothing
# end

"""
    edges_from_bonds(bonds, mat, pos)

Given an n×n adjacency matrix `bonds`, the 3×3 matrix of the cell `mat` and the
3×n matrix `pos` whose columns correspond to the positions of the atoms, extract
the list of PeriodicEdge3D corresponding to the bonds.
Since the adjacency matrix wraps bonds across the boundaries of the cell, the edges
are extracted so that the closest representatives are chosen to form bonds.
"""
function edges_from_bonds(bonds, mat, pos)
    @assert size(pos)[1] == 3
    n = size(pos)[2]
    edges = PeriodicEdge3D[]
    ref_dst = norm(mat*[1, 1, 1])
    for i in 1:n, k in (i+1):n
        bonds[k,i] || continue
        offset::Vector{SVector{3, Int}} = []
        old_dst = ref_dst
        for ofsx in -1:1, ofsy in -1:1, ofsz in -1:1
            # dst = norm(mat*(pos[:,i] .- (pos[:,k] .+ [ofsx, ofsy, ofsz])))
            dst = norm(pos[:,i] .- (pos[:,k] .+ (mat * [ofsx, ofsy, ofsz])))
            if abs2(dst - old_dst) < 1e-3
                push!(offset, (ofsx, ofsy, ofsz))
                old_dst = (dst + old_dst)/2
            elseif dst < old_dst
                empty!(offset)
                push!(offset, (ofsx, ofsy, ofsz))
                old_dst = dst
            end
        end
        for ofs in offset
            push!(edges, (i, k, ofs))
        end
    end
    return edges
end

"""
    Clusters

Classification of the atoms of a crystalline framework in different clusters.
For simple crystals, every atom is its own cluster.
For a MOF, a cluster is a SBU, which can be either organic or inorganic.
"""
struct Clusters
    sbus::Vector{Vector{PeriodicVertex3D}}
    classes::Vector{Int}
    attributions::Vector{Int}
    offsets::Vector{SVector{3,Int}}
end

function Clusters(n)
    sbus = [[PeriodicVertex3D(i)] for i in 1:n]
    classes = collect(1:n)
    attributions = collect(1:n)
    offsets = [zero(SVector{3,Int}) for _ in 1:n]
    return Clusters(sbus, classes, attributions, offsets)
end

Base.isempty(c::Clusters) = c.attributions == 1:length(c.attributions)

"""
    Crystal

Intermediate representation of a crystal, retaining information on the cell, and the
exact placement of the atoms and their type, as well as the residues which will be used as
vertices for the computation of the underlying topology.
"""
struct Crystal{T<:Union{Nothing,Clusters}}
    cell::Cell
    types::Vector{Symbol}
    clusters::T
    pos::Matrix{Float64}
    graph::PeriodicGraph3D
end

function ==(c1::Crystal{T}, c2::Crystal{T}) where T
    c1.cell == c2.cell && c1.types == c2.types && c1.clusters == c2.clusters &&
    c1.pos == c2.pos && c1.graph == c2.graph
end

Crystal{Nothing}(c::Crystal{Nothing}) = c
function Crystal{Nothing}(c::Crystal{Clusters})
    Crystal{Nothing}(c.cell, c.types, nothing, c.pos, c.graph)
end
function Crystal{Clusters}(c::Crystal, clusters::Clusters)
    Crystal{Clusters}(c.cell, c.types, clusters, c.pos, c.graph)
end

# function Crystal(cif::CIF)
#     c::CIF = expand_symmetry(cif)
#     if iszero(cif.bonds)
#         # The CIF files does not contain any bond, so we have to guess
#         atoms::Vector{Symbol} = unique!(sort(c1[].types))
#         if atoms == [:O, :Si] || atoms == [:Si] # zeolite
#             @info "Interpreting .cif file as a zeolite"
#             set_unique_bond_type!(c, 3.1, (:Si, :Si), (:Si,), 0.1)
#         elseif atoms == [:C] || atoms == [:C, :H]
#             set_unique_bond_type!(c, 1.54, (:C, :C), (:C,), 0.3)
#         elseif atoms == [:Cl, :Na] # toy example
#             set_unique_bond_type!(c, 2.82, (:Cl, :Na), (:Cl, :Na), 0.2)
#         else
#             error("Missing bonds on CIF object $(cif.cifinfo["data"])")
#         end
#     end
#     @assert !iszero(c.bonds)
#     bondedatoms::Vector{Int} = [i for i in 1:length(c.ids) if count(c.bonds[:,i]) > 1]
#     c = keep_atoms(c, bondedatoms) # Only keep atoms that have at least one bond
#     n = length(c.ids)
#     clusters = Clusters(n)
#     edges = edges_from_bonds(c.bonds, c.cell.mat, c.pos)
#     @assert !isempty(edges)
#     Crystal(c.cell, [c.types[x] for x in c.ids], clusters, c.pos, PeriodicGraph3D(n, edges))
# end

"""
    equilibrium(g::PeriodicGraph)

Return an equilibrium placement for the vertices of the graph, defined as a set
of positions such that each vertex is at the barycentre of its neighbors.
The returned equilibrium placement is such that the first vertex of the graph
is at the origin of the space.
"""
function equilibrium(g::PeriodicGraph{N}) where N
    n = nv(g)
    iszero(n) && return Matrix{Rational{Int128}}(undef, N, 0)
    Y = Matrix{Int}(undef, n, N)
    A = spzeros(Int, n, n)
    neigh = Vector{Int}(undef, n)
    offset = SizedVector{N,Int}(undef)
    for i in 1:n
        neigh .= 0
        offset .= 0
        count = 0
        for k in neighbors(g, i)
            k.v == i && continue
            count += 1
            neigh[k.v] += 1
            offset .-= k.ofs
        end
        Y[i,:] .= offset
        A[i,:] .= neigh
        A[i,i] = -count
    end

    return dixon_solve(Val(N), A[2:end,2:end], Y[2:end,:])
end

"""
    trim_topology(graph::PeriodicGraph)

Return a pair `(vmap, newgraph)` extracted from the input by removing vertices
of valence lower or equal to 1, and by replacing vertices of valence 2 by edges.

`vmap` maps the vertices of `newgraph` to their counterpart in `graph`.
"""
function trim_topology(graph::PeriodicGraph{N}) where N
    newgraph = PeriodicGraph{N}(nv(graph), collect(edges(graph)))
    remove_idx = Int[]
    flag = any(<=(2), degree(newgraph))
    vmap = collect(1:nv(graph))
    while flag # we alternate cycles when we remove valence 1 and valence 2
        # until no such vertex remains
        flag = any(<=(1), degree(newgraph))
        while flag
            n = nv(newgraph)
            for i in 1:n
                if degree(newgraph, i) <= 1
                    push!(remove_idx, i)
                end
            end
            map = rem_vertices!(newgraph, remove_idx)
            vmap = vmap[map]
            #pos = pos[vmap]
            empty!(remove_idx)
            flag = any(isone, degree(newgraph))
        end
        flag = any(==(2), degree(newgraph))

        while flag
            n = nv(newgraph)
            for i in 1:n
                if degree(newgraph, i) == 2
                    neigh1, neigh2 = neighbors(newgraph, i)
                    add_edge!(newgraph, PeriodicEdge{N}(neigh1.v, neigh2.v, neigh2.ofs .- neigh1.ofs))
                    push!(remove_idx, i)
                end
            end
            map = rem_vertices!(newgraph, remove_idx)
            vmap = vmap[map]
            #pos = pos[vmap]
            empty!(remove_idx)
            flag = any(==(2), degree(newgraph))
        end
        flag = any(<=(2), degree(newgraph))
    end
    return vmap, newgraph
end

trimmed_crystal(c::Crystal{Clusters}) = trimmed_crystal(coalesce_sbus(c))
function trimmed_crystal(c::Crystal{Nothing})
    vmap, graph = trim_topology(c.graph)
    return Crystal(c.cell, c.types[vmap], nothing, c.pos[:,vmap], graph)
end

"""
    CrystalNet{T<:Real}

Representation of a net as a topological abstraction of a crystal.

`T` is the numeric type used to store the exact coordinates of each vertex at the
equilibrium placement.
"""
struct CrystalNet{T<:Real}
    cell::Cell
    types::Vector{Symbol}
    pos::Vector{SVector{3,T}}
    graph::PeriodicGraph3D
end

function CrystalNet{T}(cell::Cell, types::AbstractVector{Symbol}, graph::PeriodicGraph3D,
                       placement::AbstractMatrix{T}) where T<:Real
    n = nv(graph)
    pos = Vector{SVector{3,T}}(undef, n)
    offsets = Vector{SVector{3,Int}}(undef, n)
    for (i, x) in enumerate(eachcol(placement))
        offsets[i] = floor.(Int, x)
        pos[i] = x .- offsets[i]
    end
    s = sortperm(pos)
    pos = pos[s]
    types = Symbol[types[s[i]] for i in 1:n]
    cell = Cell(cell, EquivalentPosition[])
    graph = offset_representatives!(graph, .-offsets)[s]
    # @assert all(pos[i] == mean(pos[x.v] .+ x.ofs for x in neighbors(graph, i)) for i in 1:length(pos))
    return CrystalNet(cell, types, pos, graph)
end

macro tryinttype(inttype, m, M, cell, types, graph, placement)
    quote
        if ((typemin($(esc(inttype))) <= $(esc(m))) & ($(esc(M)) <= typemax($(esc(inttype)))))
            return CrystalNet{Rational{$(esc(inttype))}}(
                        $(esc(cell)), $(esc(types)), $(esc(graph)),
                        Rational{$(esc(inttype))}.($(esc(placement))))
        end
    end
end

function trim_crystalnet!(graph, types, l, keep)
    tohandle::Vector{Int} = reduce(vcat, l)
    sort!(tohandle)
    toremove = keep ? deleteat!(collect(1:length(types)), tohandle) : tohandle
    vmap = rem_vertices!(graph, toremove)
    return vmap
    nothing
end

function CrystalNet(cell::Cell, types::AbstractVector{Symbol}, graph::PeriodicGraph3D)
    vmap, graph = trim_topology(graph)
    types = types[vmap]
    dimensions = PeriodicGraphs.dimensionality(graph)
    if haskey(dimensions, 0)
        @ifwarn @warn "Removing complex structure of dimension 0, possibly solvent residues."
        vmap = trim_crystalnet!(graph, types, dimensions[0], false)
        types = types[vmap]
        dimensions = PeriodicGraphs.dimensionality(graph)
    end

    if isempty(dimensions)
        return CrystalNet{Rational{Bool}}(cell, types, graph, Matrix{Rational{Bool}}(undef, 3, 0))
    end
    if !haskey(dimensions, 3)
        error("The input does not contain a 3D crystal.")
    end
    if length(dimensions) > 1 || only(keys(dimensions)) != 3
        @ifwarn @warn "Presence of periodic structures with a periodicity different than 3. At the moment, all substructures with periodicity other than 3 will be ignored."
        vmap = trim_crystalnet!(graph, types, get(dimensions, 3, Vector{Int}[]), true)
        types = types[vmap]
        dimensions = PeriodicGraphs.dimensionality(graph)
    end
    @assert !isempty(get(dimensions, 3, Vector{Int}[]))
    dim3 = dimensions[3]
    if length(dim3) > 1
        error("Multiple entertwined tridimensional structures. This is not currently handled.")
    end
    @assert sort!(only(dim3)) == collect(1:nv(graph))
    placement = equilibrium(graph)
    m = min(minimum(numerator.(placement)), minimum(denominator.(placement)))
    M = max(maximum(numerator.(placement)), maximum(denominator.(placement)))
    @tryinttype(Int8,   m, M, cell, types, graph, placement)
    @tryinttype(Int16,  m, M, cell, types, graph, placement)
    @tryinttype(Int32,  m, M, cell, types, graph, placement)
    @tryinttype(Int64,  m, M, cell, types, graph, placement)
    @tryinttype(Int128, m, M, cell, types, graph, placement)
    return CrystalNet{Rational{BigInt}}(cell, types, graph, placement)
    # Type-unstable function, but yields better performance than always falling back to Int128
end

function CrystalNet(c::Crystal{T}) where T
    if T === Clusters # the only alternative is c.clusters === nothing
        c = coalesce_sbus(c)
    end
    if ne(c.graph) == 0
        error("Empty graph. This probably means that the bonds have not been correctly attributed. Please switch to an input file containing explicit bonds.")
    end
    CrystalNet(c.cell, c.types, c.graph)
end

"""
    @enum ClusteringMode

Selection mode for the clustering of vertices. The choices are:
-   `InputClustering`: use the input residues as clusters. Fail if some atom does
    not belong to a residue.
-   `EachVertexClustering`: each vertex is its own cluster.
-   `MOFClustering`: discard the input residues and consider the input as a MOF. Identify
    organic and inorganic clusters using a simple heuristic based on the atom types.
-   `GuessClustering`: discard the input residues and try to identify the clusters as in `MOFClustering`.
    If it fails, use `EachVertexClustering`.
-   `AutomaticClustering`: if the input assigns each atom to a residue, use these residues
    as clusters. Otherwise, try to guess the clusters as in `GuessClustering`.
"""
@enum ClusteringMode begin
    InputClustering
    EachVertexClustering
    MOFClustering
    GuessClustering
    AutomaticClustering
end

function do_clustering(c::Crystal{T}, mode::ClusteringMode)::Tuple{Clusters, CrystalNet} where T
    n = length(c.types)
    if mode == InputClustering
        if T === Nothing
            throw(ArgumentError("Cannot use input residues as clusters: the input does not have residues."))
        else
            return Clusters(n), CrystalNet(c)
        end
    elseif mode == EachVertexClustering
        return Clusters(n), CrystalNet(Crystal{Nothing}(c))
    elseif mode == MOFClustering
        clusters = find_sbus(c)
        if length(clusters.sbus) <= 1
            throw(MissingAtomInformation("MOFClustering leads to a single cluster, choose a different clustering mode."))
        end
        return clusters, CrystalNet(Crystal{Clusters}(c, clusters))
    elseif mode == GuessClustering
        crystal = Crystal{Nothing}(c)
        try
            clusters, net = do_clustering(crystal, MOFClustering)
            if nv(net.graph) > 1
                return clusters, net
            end
        catch e
            if !(e isa MissingAtomInformation)
                rethrow()
            end
        end
        return do_clustering(crystal, EachVertexClustering)
    elseif mode == AutomaticClustering
        if T === Clusters
            return Clusters(n), CrystalNet(c)
        else
            return do_clustering(c, GuessClustering)
        end
    end
end

function CrystalNet(g::Union{PeriodicGraph3D,AbstractString,AbstractVector{PeriodicEdge3D}})
    graph = PeriodicGraph3D(g)
    cell = Cell()
    n = nv(g)
    types = [Symbol("") for _ in 1:n]
    return CrystalNet(cell, types, graph)
end

# function CrystalNet(cif::CIF)
#     @assert isempty(cif.cell.equivalents) # FIXME we only handle P1 representations
#     n = length(cif.ids)
#     @assert size(cif.pos)[1] == 3
#     @assert size(cif.pos)[2] == n
#
#     edges = edges_from_bonds(cif.bonds, cif.cell.mat, cif.pos)
#     pos = SVector{3,Rational{Int}}[rationalize.(Int,x) for x in eachcol(cif.pos)]
#     return CrystalNet(cif.cell, [cif.types[x] for x in cif.ids], pos, PeriodicGraph3D(edges))
# end


text:
function normalize_theta!(scales::AbstractArray, θ::AbstractArray)
    @assert length(scales) == size(θ, 1)
    @inbounds for (i, ti) in enumerate(eachrow(θ))
        scales[i] = norm(ti, 2)
        normalize!(ti, 2)
    end
    return
end

function rescale_xi!(Ξ::AbstractArray, scales::AbstractArray)
    @assert length(scales) == size(Ξ, 1)
    @inbounds for (si, ti) in zip(scales, eachrow(Ξ))
        ti .= ti ./ si
    end
    return
end

function rescale_theta!(θ::AbstractArray, scales::AbstractArray)
    @assert length(scales) == size(θ, 1)
    @inbounds for (i, ti) in enumerate(eachrow(θ))
        ti .= ti .* scales[i]
    end
    return
end

"""
    sparse_regression(X, Y, basis, p, t, maxiter, opt, denoise, normalize, convergence_error)
    sparse_regression!(Xi, X, Y, basis, p, t, maxiter, opt, denoise, normalize, convergence_error)
    spares_regression!(Xi, Theta, Y, maxiter, opt, denoise, normalize, convergence_error)

Performs a sparse regression via the algorithm `opt <: AbstractOptimizer`. `maxiter` specifies the upper bound of the iterations
of the optimizer, `convergence_error` the breaking condition due to convergence.
`denoise` defines if the matrix holding candidate trajectories should be thresholded via the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
`normalize` normalizes the matrix holding candidate trajectories via the L2-Norm over each function.

If the data matrices `X`, `Y` are given with a `Basis` `basis` and the additional information for parameters `p` and time points of the
measurements `t`, it returns the coefficient matrix `Xi` and the iterations taken.
This function is also available in place, which returns just the iterations.

If `Xi`, `Theta` and `Y` are given, the sparse regression will find the coefficients `Xi`, which minimize the objective and return the iterations needed.

# Example

```julia
opt = STRRidge()
maxiter = 10
c_error = 1e-3

Xi, iters = sparse_regression(X, Y, basis, [], [], maxiter, opt, false, false, c_error)

iters = sparse_regression!(Xi,X, Y, basis, [], [], maxiter, opt, false, false, c_error)

Xi2 = zeros(size(Y, 1), size(X, 1))
iters = sparse_regression!(Xi2, X, Y, maxiter, opt, false, false, c_error)
```
"""
function sparse_regression(X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray, t::AbstractVector , maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)

    Ξ = zeros(eltype(X), length(Ψ), ny)
    scales = ones(eltype(X), length(Ψ))
    θ = zeros(eltype(X), length(Ψ), nm)
    Ψ(θ, X, p, t)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing

    return Ξ, iters
end

function sparse_regression!(Ξ::AbstractArray, X::AbstractArray, Ẋ::AbstractArray, Ψ::Basis, p::AbstractArray , t::AbstractVector, maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer
    @assert size(X)[end] == size(Ẋ)[end]
    nx, nm = size(X)
    ny, nm = size(Ẋ)
    @assert size(Ξ) == (length(Ψ), ny)

    scales = ones(eltype(X), length(Ψ))
    θ = zeros(eltype(X), length(Ψ), nm)
    Ψ(θ, X, p, t)

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing

    return iters
end

# For pareto
function sparse_regression!(Ξ::AbstractArray, θ::AbstractArray, Ẋ::AbstractArray, maxiter::Int64 , opt::T, denoise::Bool, normalize::Bool, convergence_error) where T <: Optimize.AbstractOptimizer

    scales = ones(eltype(Ξ), size(θ, 1))

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    Optimize.init!(Ξ, opt, θ', Ẋ')'
    iters = Optimize.fit!(Ξ, θ', Ẋ', opt, maxiter = maxiter, convergence_error = convergence_error)

    normalize ? rescale_xi!(Ξ, scales) : nothing
    normalize ? rescale_theta!(θ, scales) : nothing

    return iters
end

# One Variable on multiple derivatives
function SINDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis, opt::T = STRRidge(); kwargs...) where {T <: Optimize.AbstractOptimizer, S <: Number}
    return SINDy(X', Ẋ, Ψ, opt; kwargs...)
end

# Multiple on one
function SINDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis, opt::T = STRRidge(); kwargs...) where {T <: Optimize.AbstractOptimizer, S <: Number}
    return SINDy(X, Ẋ', Ψ, opt; kwargs...)
end

# General
"""
    SINDy(X, Y, basis, opt = STRRidge(); p, t, maxiter, convergence_error, denoise, normalize)
    SINDy(X, Y, basis, lambdas, opt = STRRidge(); f, g, p, t, opt, maxiter, convergence_error, denoise, normalize)

Performs Sparse Identification of Nonlinear Dynamics given the data matrices `X` and `Y` via the `AbstractBasis` `basis.`
Keyworded arguments include the parameter (values) of the basis `p` and the timepoints `t` which are passed in optionally.
`opt` is an `AbstractOptimizer` useable for sparse regression, `maxiter` the maximum iterations to perform and `convergence_error` the
bound which causes the optimizer to stop.
`denoise` defines if the matrix holding candidate trajectories should be thresholded via the [optimal threshold for singular values](http://arxiv.org/abs/1305.5870).
`normalize` normalizes the matrix holding candidate trajectories via the L2-Norm over each function.

Typically `X` represent the state measurements and `Y` the measurements of the differential state. Since only the number of measurements (column dimension of the matrices) have to be equal, it is possible to augment `X` with additional
data, e.g. external forcing or inputs.

If `SINDy` is called with an additional array of thresholds contained in `lambdas`, it performs a multi objective optimization over all thresholds.
The best candidate is determined via the mapping onto a feature space `f` and an (scalar, positive definite) evaluation `g`.
The signature of should be `f(xi, theta, yi)` where `xi` are the coefficients of the sparse optimization,`theta` is the evaluated candidate library and `yi` are the rows of the matrix `Y`.
Returns a `SINDyResult`. If the pareto optimization is used, the result combines the best candidate for each row of `Y`.
"""
function SINDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 2}, Ψ::Basis, opt::T = STRRidge(); p::AbstractArray = [], t::AbstractVector = [], maxiter::Int64 = 10, denoise::Bool = false, normalize::Bool = true, convergence_error = eps()) where {T <: Optimize.AbstractOptimizer, S <: Number}
    Ξ, iters = sparse_regression(X, Ẋ, Ψ, p, t, maxiter, opt, denoise, normalize, convergence_error)
    convergence = iters < maxiter
    SparseIdentificationResult(Ξ, Ψ, iters, opt, convergence, Ẋ, X, p = p)
end

function SINDy(X::AbstractArray{S, 1}, Ẋ::AbstractArray, Ψ::Basis, thresholds::AbstractArray, opt::T = STRRidge(); kwargs...) where {T <: Optimize.AbstractOptimizer, S <: Number}
    return SINDy(X', Ẋ, Ψ, thresholds, opt; kwargs...)
end

function SINDy(X::AbstractArray{S, 2}, Ẋ::AbstractArray{S, 1}, Ψ::Basis, thresholds::AbstractArray, opt::T = STRRidge()::AbstractArray; kwargs...) where {T <: Optimize.AbstractOptimizer, S <: Number}
    return SINDy(X, Ẋ', Ψ, thresholds, opt; kwargs...)
end

function SINDy(X, DX, Ψ::Basis, thresholds::AbstractArray, opt::T = STRRidge(); f::Function = (xi, theta, dx)->[norm(xi, 0); norm(dx .- theta'*xi, 2)], g::Function = x->norm(x), p::AbstractArray = [], t::AbstractVector = [], maxiter::Int64 = 10, denoise::Bool = false, normalize::Bool = true, convergence_error = eps()) where {T <: Optimize.AbstractOptimizer}
    @assert size(X)[end] == size(DX)[end]
    nx, nm = size(X)
    ny, nm = size(DX)

    θ = Ψ(X, p, t)

    scales = ones(eltype(X), length(Ψ))

    ξ_tmp = zeros(eltype(X), length(Ψ), ny)
    Ξ_opt = zeros(eltype(X), length(Ψ), ny)

    iter_ = Inf
    _iter = Inf
    _thresh = Inf

    denoise ? optimal_shrinkage!(θ') : nothing
    normalize ? normalize_theta!(scales, θ) : nothing

    fg(xi, theta, dx) = (g∘f)(xi, theta, dx)

    @inbounds for j in 1:length(thresholds)
        set_threshold!(opt, thresholds[j])
        iter_ = sparse_regression!(view(ξ_tmp, :, :), view(θ, :, :), view(DX, :, :), maxiter, opt, false, false, convergence_error)
        normalize ? rescale_xi!(ξ_tmp, scales) : nothing
        for i in 1:ny
            if j == 1
                Ξ_opt .= ξ_tmp
                _iter = iter_
                _thresh = thresholds[j]
            else
                if evaluate_pareto!(view(Ξ_opt, :, i), view(ξ_tmp, :, i), fg, view(θ, :, :), view(DX, i, :))
                    _iter = iter_
                    _thresh = thresholds[j]
                end
            end
        end
    end

    set_threshold!(opt, _thresh)
    return SparseIdentificationResult(Ξ_opt, Ψ, _iter, opt, _iter < maxiter, DX, X, p = p)
end

ext: jl


text:
@doc raw"""
```julia
abstract type Abstract_Map end
```

This type is an abstraction of a
```math
f:X ↦ Y=f_{\hat{\theta}}(X)
```
where ``\hat{\theta}`` is map parameter vector.

Such transformations can be used in various contexts:

- calibration task: we have a reference ``\hat{X}`` to a calibrated
  ``X``, see [`Model2Fit_Recalibration`](@ref).

- parameter transformation: we map a parameter vector ``\hat{θ}`` to a new one ``θ``,
  see [`Model2Fit_Transformed_Parameters`](@ref).

# Interface

- [`parameter_size`](@ref)
- [`eval_map`](@ref)

"""
abstract type Abstract_Map end

@doc raw"""
```julia
parameter_size(::Abstract_Map) -> Int
```

Return ``\hat{θ}`` parameter vector length

# See also
- [`Abstract_Map`](@ref)

"""
parameter_size(::Abstract_Map) = @assert(false,"To implement")

@doc raw"""
```julia
eval_map(m::Abstract_Map,X_hat::AbstractVector,hat_θ::AbstractVector) -> X::AbstractVector
```

Compute ``X=X(\hat{X})``.

# See also
- [`Abstract_Map`](@ref)

"""
eval_map(m::Abstract_Map,X_hat::AbstractVector,hat_θ::AbstractVector) = @assert(false,"To implement")


text:
using LinearAlgebra

struct CountedVector{T<:Real}
    v::Vector{T}
    Av::Vector{T}
    id::Int

    function CountedVector{T}(v::Vector{T}, A, id) where {T}
        canonicalize!(v)
        new{T}(v, A*v, id)
    end
end

CountedVector(v::Vector{T}, A, id) where {T} = CountedVector{T}(v, A, id)

Base.vec(c::CountedVector) = c.v

mutable struct DoubleDescription{T<:Real}
    A::Matrix{T}
    R::Vector{CountedVector{T}}
    K::Set{Int}
    adj::Dict{Tuple{Int,Int},Bool}
    num_rays::Int
end

function double_description(ine::LiftedHRepresentation{T}) where {T}
    if !isempty(ine.linset)
        error("Linearity currently unsupported by ConvexHull.")
    end
    # FIXME add support for linearity
    dd = double_description(ine.A)
    R = Matrix{T}(undef, length(dd.R), fulldim(ine)+1)
    for (i,r) in enumerate(dd.R)
        R[i,:] = r.v
    end
    return LiftedVRepresentation{T}(R)
end

function double_description(ext::LiftedVRepresentation{T}) where {T}
    if !isempty(ext.linset)
        error("Linearity currently unsupported by ConvexHull.")
    end
    # FIXME add support for linearity
    dd = double_description(ext.R)
    A = Matrix{T}(undef, length(dd.R), fulldim(ext)+1)
    for (i,r) in enumerate(dd.R)
        A[i,:] = r.v
    end
    return LiftedHRepresentation(A)
end

function initial_description(A::Matrix{T}) where {T<:Real}
    m,n = size(A)
    B = rref(Matrix(A'))
    # find pivots
    r = 1
    K = Set{Int}()
    for i in 1:m
        comp = zeros(n)
        comp[r] = 1
        if B[:,i] == comp
            r += 1
            push!(K, i)
        end
        r > n && break
    end
    cK = sort!(collect(K))
    Ak = A[cK,:]
    if eltype(Ak) <: AbstractFloat
        R = Ak \ Matrix(LinearAlgebra.I, n, n)
    else
        # Ak \ eye(n,n) creates BigFloat from Rational{BigInt} while inv keeps Rational{BigInt}
        R = inv(Ak)
    end
    Rk = [CountedVector{T}(R[:,i],A,i) for i in 1:n]
    dd = DoubleDescription(A,Rk,K,Dict{Tuple{Int,Int},Bool}(),n)
    for i in 1:n
        # Ar = Rk[i].Av[cK]
        for j in (i+1):n
            # As = Rk[j].Av[cK]
            id = extrema([Rk[i].id,Rk[j].id])
            cache_adjacency!(dd, n, Rk[i].Av[cK], Rk[j].Av[cK], id)
        end
    end
    return dd
end

function double_description(A::Matrix{T}) where {T<:Real}
    A = [zeros(T,1,size(A,2)); A]
    A[1,1] = one(T)
    m, n = size(A)
    dd = initial_description(A)
    Kᶜ = setdiff(1:m, dd.K)
    while !isempty(Kᶜ)
        i = pop!(Kᶜ)
        # length(dd.K) == 31 && break
        update!(dd, i)
    end
    return dd
end

function is_approx_included(haystack, needle)
    n = length(needle)
    diff = zeros(n)
    for h in haystack
        diff = norm(vec(h)-needle)
        diff < n*ε && return true
    end
    return false
end

function canonicalize!(v::Vector{T}) where {T<:Real}
    n = length(v)
    val = abs(v[1])
    if val < ε
        val = abs(v[findfirst(abs.(v) .> ε)])
    end
    for i in 1:n
        v[i] = v[i] / val
    end
    return v
end

# use Lemma 8 from Fukuda (1996) to update the double description
function update!(dd::DoubleDescription{T}, i) where {T<:Real}
    m, n = size(dd.A)
    Aᵢ = reshape(dd.A[i,:], (n,))
    Rⁿᵉʷ = CountedVector{T}[]
    R⁺, R⁰, R⁻ = partition_rays(dd.R, Aᵢ)
    for r in R⁺, s in R⁻
        if isadjacent(dd,r,s)
            w = dot(Aᵢ,vec(r))*vec(s) - dot(Aᵢ,vec(s))*vec(r)
            dd.num_rays += 1
            v = CountedVector(w,dd.A,dd.num_rays)
            if sum(abs.(w)) > n*ε &&
               !is_approx_included(R⁰,   vec(w)) &&
               !is_approx_included(Rⁿᵉʷ, vec(w))
                dd.num_rays += 1
                push!(Rⁿᵉʷ, v)
            end
        end
    end
    dd.R = vcat(R⁺, R⁰, Rⁿᵉʷ)
    push!(dd.K, i)
    cK = sort(collect(dd.K))
    Ak = dd.A[cK,:]
    # should really add a test right about here to ensure
    # that old rays do not become adjacent...I think this
    # can only happen if both v,w ∈ R⁰

    # If eltype(Ak) is Rational{BigInt}, LinearAlgebra.jl is needed
    d = rank(Ak)
    for s in Rⁿᵉʷ
        # As = s.Av[cK]#Ak*vec(s)
        for r in dd.R
            r.id == s.id && continue
            # Ar = r.Av[cK]#Ak*vec(r)
            id = extrema([r.id, s.id])
            cache_adjacency!(dd, d, r.Av[cK], s.Av[cK], id)
        end
    end
    nothing
end

function partition_rays(R::Vector{CountedVector{T}}, a::Vector{T}) where T<:Real
    R⁺, R⁰, R⁻ = CountedVector{T}[], CountedVector{T}[], CountedVector{T}[]
    n = length(a)
    for r in R
        if sum(abs.(vec(r))) < n*ε
            println("have a zero vector!")
            continue
        end
        val = dot(a,vec(r))
        if val > ε
            push!(R⁺, r)
        elseif val < -ε
            push!(R⁻, r)
        else
            push!(R⁰, r)
        end
    end
    return R⁺, R⁰, R⁻
end

isadjacent(dd, r, s) = dd.adj[extrema([r.id,s.id])]

function cache_adjacency!(dd, d, Ar, As, id::Tuple{Int,Int})
    Z = active_sets(dd, Ar, As)
    if length(Z) < d - 2
        val = false
    elseif length(intersect(Z,dd.K)) ≥ d - 2
        val = true
    else
        val = (rank(dd.A[sort(collect(dd.K)),:][Z,:]) == d - 2)
    end
    dd.adj[id] = val
    return val
end

function active_sets(dd, Ar, As)
    Z = Int[]
    m = length(Ar)
    sizehint!(Z,m)
    for i in 1:m
        if abs(Ar[i]) ≤ ε && abs(As[i]) ≤ ε
            push!(Z, i)
        end
    end
    return Z
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 659840fe-502b-11eb-1c22-b971c3425be2
md"""
# Stacks & Queues
"""

# ╔═╡ 06ac35da-502b-11eb-01c1-1974550a834f
md"""
## Stacks

Implemented with:
- An Array where the last element is top of the Stack following LIFO (Last In First Out).
- A singly linked list can also be used to create a stack as we'll only have to track the front of the list, and keep adding element on top of it.

Operations:
- Push on top of the stack
- Pop from the top of the stack

"""

# ╔═╡ 18dd6bc2-502b-11eb-38e5-e900758f4b61
# Eg: Balanced Paranthesis

# ╔═╡ 7607b956-502b-11eb-01d4-ff83de1d9425
# Eg: Reverse Polish Notation

# ╔═╡ 7f79655c-502b-11eb-31f0-0588189afd4b
# Reversing a Stack

# ╔═╡ 872861fe-502b-11eb-3560-8f03be0bd627
# Minimum Bracket Reversals

# ╔═╡ 8d21363a-502b-11eb-1ee2-01b8cd790c77
md"""
## Queue

Has both Head and Tail

Implementation with
- Linked List : save references for both head and tail
- Array:
- Using a Stack:

Operations:
- Peek the head: no removal
- Dequeue from the head: get the element and remove the element
- Enqueue to the tail : add to element to the tail of the queue

"""

# ╔═╡ f41e2c08-502b-11eb-059e-57e5a1649134

# ╔═╡ f6019938-502b-11eb-267d-0b8e833e7bc5
# Reversing a Queue

# ╔═╡ c87d9b24-502b-11eb-08bb-63a08de5950d
md"""
### Deques: Double Ended Queues

Double Ended queue that goes both ways, i.e. enque or deque from both ends.

Like a Generalized version of both stack and queue and can be treated as stack or queue.

"""

# ╔═╡ f3086ec8-502b-11eb-3041-fde296acea0f

# ╔═╡ ca1eb40e-502b-11eb-0d3c-81d1bf5ee4b7
md"""
### Priority Queue

Elements are enqueued with a priority and are dequeued with the highest priority first
"""

# ╔═╡ f2624340-502b-11eb-3b70-ad34edaf5a99

# ╔═╡ Cell order:
# ╟─659840fe-502b-11eb-1c22-b971c3425be2
# ╟─06ac35da-502b-11eb-01c1-1974550a834f
# ╠═18dd6bc2-502b-11eb-38e5-e900758f4b61
# ╠═7607b956-502b-11eb-01d4-ff83de1d9425
# ╠═7f79655c-502b-11eb-31f0-0588189afd4b
# ╠═872861fe-502b-11eb-3560-8f03be0bd627
# ╠═8d21363a-502b-11eb-1ee2-01b8cd790c77
# ╠═f41e2c08-502b-11eb-059e-57e5a1649134
# ╠═f6019938-502b-11eb-267d-0b8e833e7bc5
# ╠═c87d9b24-502b-11eb-08bb-63a08de5950d
# ╠═f3086ec8-502b-11eb-3041-fde296acea0f
# ╠═ca1eb40e-502b-11eb-0d3c-81d1bf5ee4b7
# ╠═f2624340-502b-11eb-3b70-ad34edaf5a99


text:
"""
    ngcd(p,q, [k]; kwargs...)

Find numerical GCD of polynomials `p` and `q`. Refer to [`NGCD.ngcd`](@ref) for details.

In the case `degree(p) ≫ degree(q)`,  a heuristic is employed to first call one step of the Euclidean gcd approach, and then call `ngcd` with relaxed tolerances.

"""
function ngcd(p::P, q::Q, args...;kwargs...) where {T, S, P<:StandardBasisPolynomial{T}, Q <: StandardBasisPolynomial{S}}

    degree(p) < 0  && return (u=q,      v=p, w=one(q),  θ=NaN, κ=NaN)
    degree(p) == 0 && return (u=one(q), v=p, w=zero(q), θ=NaN, κ=NaN)
    degree(q) < 0  && return (u=one(q), v=p, w=zero(q), θ=NaN, κ=NaN)
    degree(q) == 0 && return (u=one(p), v=p, w=q,       θ=NaN, κ=NaN)
    assert_same_variable(p,q)

    p′,q′ = promote(p,q)

    if degree(p′) > 5*degree(q′)
        return  ngcd′(p′, q′, args...; kwargs...)
    end

    ps, qs = [p′[i] for i in eachindex(p′)], [q′[i] for i in eachindex(q′)] # need vectors, want copy
    u,v,w,Θ,κ = NGCD.ngcd(ps, qs, args...; kwargs...)
    P′ = ⟒(typeof(p′))
    (u=P′(u, indeterminate(p)), v=P′(v, indeterminate(p)), w = P′(w, indeterminate(p)), θ=Θ, κ=κ)

end

"""
    ngcd′(p,q)

When degree(p) ≫ degree(q), this uses a early call to `divrem` to bring about commensurate degrees
before calling `ngcd`.
"""
function ngcd′(p::StandardBasisPolynomial{T}, q::StandardBasisPolynomial{T};
               atol = eps(real(float(T))),
               rtol = atol,
               satol= atol,
               srtol= rtol,
               kwargs...
               ) where {T}

    a, b = divrem(p,q)

    # check if a=u (p,q) ≈ (aq,q)
    if isapprox(p, a*q, atol=atol, rtol=rtol)
        return (u=a, v=p, w=q, θ=NaN, κ=NaN)
    else
        ngcd(q, b; atol=100atol, rtol=100rtol, kwargs...)
    end
end

module NGCD
using LinearAlgebra

"""
    ngcd(ps::Vector{T}, qs::Vector{T}, [k::Int]; scale::Bool=true, atol=eps(T), rtol=eps(T), satol=atol, srtol=rtol)

Return `u, v, w, Θ, κ` where ``u⋅v ≈ ps`` and ``u⋅w ≈ qs`` (polynomial multiplication); `Θ` (`\\Theta[tab]`) is the residual error (``‖ [u⋅v,u⋅w] - [ps,qs] ‖``); and `κ` (`\\kappa[tab]`) is the numerical gcd condition number estimate. When `scale=true` (the default), ``u⋅v ≈ ps/‖ps‖`` and ``u⋅w ≈ qs/‖qs‖``

The numerical GCD problem is defined in [1] (5.4). Let ``(p,q)`` be a
polynomial pair with degree m,n. Let Ρmn be set of all such pairs. Any
given pair of polynomials has an exact greatest common divisor, ``u``, of
degree ``k``, defined up to constant factors. Let ``Ρᵏmn`` be the manifold of
all such ``(p,q)`` pairs with exact gcd of degree ``k``. A given pair ``(p,q)`` with exact gcd of degree ``j`` will
have some distance ``Θᵏ`` from ``Pᵏ``.  For a given threshold ``ϵ > 0`` a numerical GCD
of ``(p,q)`` within ``ϵ`` is an exact GCD of a pair ``(p̂,q̂)`` in ``Ρᵏ`` with

``‖ (p,q) - (p̂,q̂) ‖ <= Θᵏ``, where ``k`` is the largest value for which ``Θᵏ < ϵ``.

(In the ``ϵ → 0`` limit, this would be the exact GCD.)

Suppose ``(p,q)`` is an ``ϵ`` pertubation from ``(p̂,q̂)`` where ``(p̂,q̂)`` has an exact gcd of degree ``k``, then ``degree(gcdₑ(p,q)) = k``; as ``ϵ → 0``, ``gcdₑ(p,q) → gcd(p̂, q̂)``, and

``limsup_{(p,q)→(p̂,q̂)} inf{ ‖ (u,v,w) - (û,v̂,ŵ) ‖} / ‖ (p,q) - (p̂,q̂) ‖ < κₑ(p,q)``.

``κ`` is called the numerical GCD condition number.

The Zeng algorithm proposes a degree for ``u`` and *if* a triple ``(u,v,w)`` with ``u`` of degree ``k`` and ``(u⋅v, u⋅w)`` in ``Ρᵏmn`` can be found satisfying ``‖ (u⋅v, u⋅w) - (p,q) ‖ < ϵ`` then ``(u,v,w)`` is returned; otherwise the proposed degree is reduced and the process repeats. If not terminated, at degree ``0`` a constant gcd is returned.

The initial proposed degree is the first ``j``,  ``j=min(m,n):-1:1``, where ``Sⱼ`` is believed to have a singular value of ``0`` (``Sⱼ`` being related to the Sylvester matrix of `ps` and `qs`). The verification of the proposed degree is done using a Gauss-Newton iteration scheme holding the degree of ``u`` constant.

## Scaling:

If `scale=true` (the default when a polynomial norm is large), the gcd of ``p/‖p‖`` and ``q/‖q‖`` is identified. Scaling can reduce the condition numbers significantly.

## Tolerances:

There are two places where tolerances are utilized:

* in the identification of the rank of ``Sⱼ`` a value ``σ₁ = ‖Rx‖`` is identified. To test if this is zero a tolerance of `max(satol, ‖R‖ₒₒ ⋅ srtol)` is used.

* to test if ``(u ⋅ v, u ⋅ w) ≈ (p,q)`` a tolerance of `max(atol, λ⋅rtol)` is used with `λ` chosen to be  ``(‖(p,q)‖⋅n⋅m)⋅κ′⋅‖A‖ₒₒ`` to reflect the scale of ``p`` and ``q`` and an estimate for the condition number of ``A`` (a Jacobian involved in the solution).

This seems to work well for a reasaonable range of polynomials, however there can be issues: when the degree of ``p`` is much larger than the degree of ``q``, these choices can fail; when a higher rank is proposed, then too large a tolerance for `rtol` or `atol` can lead to a false verification; when a tolerance for `atol` or `rtol` is too strict, then a degree may not be verified.

!!! note:
    These tolerances are adjusted from those proposed in [1].

## Specified degree:

When `k` is specified, a value for ``(u,v,w)`` is identified with ``degree(u)=k``. No tolerances are utilized in computing ``Θᵏ``.

Output:

The function outputs a named tuple with names (`u`, `v`, `w`, `Θ`, `κ`). The components `u`,`v`,`w` estimate the gcd and give the divisors. The value `Θ` estimates ``Θᵏ`` and `κ` estimates the numerical condition number.

Example:

```
using Polynomials
x = variable(Polynomial{Float64})
p = (x+10)*(x^9 + x^8/3 + 1)
q = (x+10)*(x^9 + x^8/7 - 6/7)
gcd(p,q) # u a constant
gcd(p,q, method=:numerical)  # u a degree 1 polynomial
Polynomials.NGCD.ngcd(coeffs(p), coeffs(q), verbose=true) # to see some computations
```

Reference:

[1] The Numerical Greatest Common Divisor of Univariate Polynomials
by Zhonggang Zeng;
[url](http://homepages.neiu.edu/~zzeng/uvgcd.pdf);
[doi](https://doi.org/10.1090/conm/556/11014)

Note: Based on work by Andreas Varga; Requires `VERSION >= v"1.2"`.

"""
function ngcd(ps::Vector{T},
              qs::Vector{T};
              scale::Bool=false, #norm(ps) > 1e6, #false,
              atol = eps(real(T)),
              rtol = Base.rtoldefault(real(T)),
              satol = eps(real(T))^(5/6),
              srtol = eps(real(T)),
              verbose=false,
              minⱼ = -1
              ) where {T <: AbstractFloat}

    m, n = _degree.((ps, qs))

    if m < n
        out = ngcd(qs, ps; atol=atol, rtol=rtol, satol=satol, srtol=srtol, verbose=verbose)
        return (u=out.u, v=out.w, w=out.v, Θ=out.Θ, κ=out.κ)
    end

    if scale
        ps = ps./norm(ps)
        qs = qs./norm(qs)
    end

    # storage
    A0 = zeros(T, m+1, 2)
    A0[:,1] = ps
    A0[end-n:end,2] = qs

    # pre-allocate storage for Sylvester Matrices, S₁, S₂...
    Q = zeros(T, m + n, m + n)
    R = zeros(T, m + n, m + n)
    Sₓ = hcat(convmtx(ps,1),  convmtx(qs, m-n+1))
    local x::Vector{T}

    j = n  # We count down Sn, S_{n-1}, ..., S₂, S₁

    F = qr(Sₓ)
    nr, nc = size(Sₓ) # m+1, m-n+2
    Q[1:nr, 1:nr] = F.Q
    R[1:nc, 1:nc] = F.R

    while true

        V = view(R, 1:nc, 1:nc)
        flag, σ, x = smallest_singular_value(V, satol *  sqrt(1 + m - j), srtol)
        verbose && println("------ degree $j ----- σ₁: $σ  --- $flag")

        if (flag == :iszero || flag == :ispossible)

            u, v, w = initial_uvw(Val(flag), j, ps, qs, x)
            flag, ρ₁, σ₂, ρ = refine_uvw!(u,v,w, ps, qs, atol, rtol)

            verbose && println("   --- Θᵏ: $ρ₁ --- $flag (ρ=$(ρ))")

            if flag == :convergence
                return (u=u, v=v, w=w, Θ=ρ₁, κ=σ₂) # (u,v,w) verified
            end
        end

        # reduce possible degree of u and try again with Sⱼ₋₁
        # unless we hit specified minimum, in which case return it
        if j == minⱼ
            u, v, w = initial_uvw(Val(:ispossible), j, ps, qs, x)
            flag, ρ₁, σ₂, ρ = refine_uvw!(u,v,w, ps, qs, atol, rtol)
            return (u=u, v=v, w=w, Θ=ρ₁, κ=σ₂)
        end

        j -= 1
        nr += 1
        nc += 2
        nc > nr && break

        extend_QR!(Q,R, nr, nc, A0) # before Q⋅R = Sⱼ, now Q⋅R = Sⱼ₋₁

    end

    # u is a constant
    verbose && println("------ GCD is constant ------")
    u, v, w = initial_uvw(Val(:constant), j, ps, qs, x)
    flag, ρ₁, κ, ρ = refine_uvw!(u,v,w, ps, qs, atol, rtol)
    return (u=u, v=v, w=w, Θ=ρ₁, κ=κ)

end

function ngcd(ps::Vector{T}, qs::Vector{S}; kwargs...) where {T, S}
    ps′,qs′ = promote(float.(ps), float(qs))
    ngcd(ps′, qs′; kwargs...)
end

# fix the degree, k
function ngcd(ps::Vector{T},
              qs::Vector{T},
              k::Int;
              kwargs...
              ) where {T <: AbstractFloat}

    m, n = _degree.((ps,qs))

    if m < n
        out = ngcd(qs, ps, k, atol=Inf, rtol=Inf)
        return (u=out.u, v=out.w, w=out.v, Θ=out.Θ, κ=out.κ)
    end

    #    u,v,w = initial_uvw(Val(:iszero), k, ps, qs, nothing)
    Sⱼ = [convmtx(ps, n-k+1) convmtx(qs, m-k+1)]
    F = qr(Sⱼ)
    flag, σ, x = smallest_singular_value(F.R, eps(T) *  sqrt(1 + m - k), eps(T))
    if flag != :iszero
        w, v = x[1:n-k+1], -x[n-k+2:end]
        u  = [convmtx(v,k+1); convmtx(w, k+1)] \ [ps; qs]
    else
        u,v,w = initial_uvw(Val(flag), k, ps, qs, nothing)
    end
    flag, ρ₁, κ, ρ = refine_uvw!(u,v,w, ps, qs, Inf, Inf)
    return (u=u, v=v, w=w, Θ=ρ₁, κ=κ)

end

function ngcd(ps::Vector{T}, qs::Vector{S}, k::Int; kwargs...) where {T, S}
    ps′,qs′ = promote(float.(ps), float(qs))
    ngcd(ps′, qs′, k; kwargs...)
end

## -----

# return guess at smallest singular value and right sinuglar value, x
# for an upper triangular matrix, V
function smallest_singular_value(V::AbstractArray{T,2},
                                 atol=eps(T),
                                 rtol=zero(T)) where {T}

    R = UpperTriangular(V)
    k = size(R)[1]/2
    if iszero(det(R))
        return (:iszero, zero(T), T[])
    end

    m,n = size(R)

    # we are testing if ‖Ax‖ ≈ 0
    # If x is a perfect 0, but x is approximated by x' due to round off
    # then ‖A(x-x')‖ <= ‖A‖⋅‖x - x'‖ so we use ‖A‖ as scale factor
    δ = max(atol,  norm(R,Inf) * rtol)

    x = ones(T, n)
    y = zeros(T, m)
    σ₀ = σ₁ = Inf*one(real(T))
    steps, minᵢ = 1, 5

    while true
        y .= R' \ x # use iteration to converge on singular value
        x .= R  \ y
        x ./= norm(x,2)
        σ₁ = norm(R * x, 2)

        if (steps <= 50) && (steps <= minᵢ || σ₁ < 0.05*σ₀) # decreasing, keep going
            σ₀ = σ₁
        else
            break
        end

        steps += 1
    end

    if σ₁ < δ
        return (:ispossible, σ₁, x)
    else
        return (:constant, σ₁, x)
    end

end

# extend QR to next size
# Q gets a 1 in nc,nc, 0s should be elswhere
function extend_QR!(Q,R, nr, nc, A0)

    #old Q is m x m, old R is n x n; we add to these
    n = nc-2
    m = nr - 1
    k,l = size(A0)

    # add two columns to R
    # need to apply Q to top part of new columns
    R[nr-k+1:nr, (nc-1):nc] = A0
    R[1:nr-1, (nc-1):nc] = (view(Q, 1:nr-1, 1:nr-1))' *  R[1:nr-1, (nc-1):nc]

    # extend Q with row and column with identity
    Q[nr,nr] = 1

    # Make R upper triangular using Givens rotations
    for j in nr-1:-1:nc-1
        gj,_ = givens(R[j,nc-1], R[j+1,nc-1], j, j+1)
        rmul!(Q, gj')
        lmul!(gj, R)
    end

    for j in nr-1:-1:nc
        gj,_ = givens(R[j,nc], R[j+1,nc], j, j+1)
        rmul!(Q, gj')
        lmul!(gj, R)
    end

    return nothing

end

## --------------------------------------------------
## Refine u,v,w

## Find u₀,v₀,w₀ from right singular vector
function initial_uvw(::Val{:ispossible}, j, ps, qs, x)

    # Sk*[w;-v] = 0, so pick out v,w after applying permuation
    m,n = _degree.((ps, qs))
    vᵢ = vcat(2:m-n+2, m-n+4:2:length(x))
    wᵢ = vcat(1, (m-n+3):2:length(x))
    v = -x[vᵢ]
    w = x[wᵢ]
    # p194 3.9 C_k(v) u = p or Ck(w) u = q; this uses 10.2
    u = [convmtx(v,j+1); convmtx(w, j+1)] \ [ps; qs]

    return u, v, w

end

function initial_uvw(::Val{:iszero}, j, ps::Vector{T}, qs, x) where {T}
    m,n = _degree.((ps,qs))
    S = [convmtx(ps, n-j+1) convmtx(qs, m-j+1)]

    F = qr(S)
    R = UpperTriangular(F.R)

    if iszero(det(R))
        x = eigvals(R)[:,1]
    else
        x = ones(T, size(R,2))
        x .= R \ (R' \ (x/norm(x)))
        x ./= norm(x)
    end

    w = x[1:n-j+1]
    v = -x[(n-j+2):end]

    u = [convmtx(v, j+1); convmtx(w, j+1)] \ [ps;qs]

    (u,v,w)
end

function initial_uvw(::Val{:constant}, j, ps::Vector{T}, qs, x) where {T}
    u = ones(T,1)
    w = qs
    v = ps
    u,v,w
end

# find estimate for σ₂, used in a condition number (κ = 1/σ)
function σ₂(J)
    F = qr(J)
    flag, σ, x = smallest_singular_value(F.R)
    σ
end

## attempt to refine u,v,w
## check that [u ⊗ v; u ⊗ w] ≈ [p; q]
function refine_uvw!(u::Vector{T}, v, w, p, q, atol, rtol) where {T}

    m, n, l =  _degree.((u, v, w))

    uv = u ⊗ v; uw = u ⊗ w
    ρ₀, ρ₁ = one(T), residual_error(p,q,uv,uw)

    # storage
    A = zeros(T, JF_size(u, v, w)...)
    b = zeros(T, 1 + length(p) + length(q))
    Δf = zeros(T, m + n + l + 3)
    steps = 0

    h, β =  u, norm(u)^2
    minᵢ, Maxᵢ = 5, 20
    κ = NaN

    JF!(A, h, u, v, w)
    Fmp!(b,  h'*u - β, p, q, uv, uw)

    while ρ₁ > 0.0

        # update A, b, then solve A\b
        qrsolve!(Δf, A, b)

        # m + n = degree(p)
        # m + l = degree(q)
        # b has length degree(p)+degree(q) + 3
        Δv, Δw, Δu = Δf[1:(n+1)], Δf[(n+1) .+ (1:l+1)], Δf[n+1+l+1 + 1:end]

        uv, uw =  (u-Δu) ⊗ (v-Δv),  (u-Δu) ⊗ (w-Δw)
        ρ₀, ρ′ = ρ₁, residual_error(p, q, uv, uw)

        # don't worry about first few, but aftewards each step must be productive
        # though we can have really bad first steps, which we cap
        if  (steps <= Maxᵢ) && (steps <= minᵢ || ρ′ < 0.95 * ρ₀) && (  ρ′ < 100*ρ₁ )
            ρ₁ = ρ′
            u .-= Δu; v .-= Δv; w .-= Δw
            steps += 1
        else
            break
        end

        # update A,b for next iteration
        JF!(A, h, u, v, w)
        Fmp!(b,  h'*u - β, p, q, uv, uw)

    end

    # this is a heuristic
    # sensitivity is Δu / Δp <= ‖ A+ ‖ = κ
    # we use an estimate for ‖(p,q)‖ error times ‖A⁺‖⋅‖A‖ₒₒ
    κ = 1/σ₂(A) # ≈ ‖A⁺‖
    λ = norm((norm(p), norm(q))) * (m * n) * min(1, κ) * norm(A, Inf)
    ρ = max(atol, rtol * λ)

    if ρ₁ <= ρ
        return :convergence, ρ₁, κ, ρ
    else
        return :no_convergence, ρ₁, κ, ρ
    end

end

function qrsolve!(y::Vector{T}, A, b) where {T}
    y .= qr(A) \ b
end

# # Fast least-squares solver for full column rank Hessenberg-like matrices
# # By Andreas Varga
function qrsolve!(y::Vector{T}, A, b) where {T <: Float64}
    Base.require_one_based_indexing(A)
    m, n = size(A)
    m < n && error("Column dimension exceeds row dimension")
    _, τ = LinearAlgebra.LAPACK.geqrf!(A)
    T <: Complex ? tran = 'C' : tran = 'T'
    LinearAlgebra.LAPACK.ormqr!('L', tran, A, τ, view(b,:,1:1))
    y .= UpperTriangular(triu(A[1:n,:]))\b[1:n]
end

## Jacobian F(u,v,w) = [p,p'] is J(u,v,w)
function JF_size(u, v, w)

    m, k, j = _degree(u), _degree(v), _degree(w)
    n, l = m + k, m + j

    ai, aj = convmtx_size(v, m + 1)
    bi, bj = convmtx_size(u, k + 1)
    di, dj = convmtx_size(w, m + 1)
    fi, fj = convmtx_size(u, j + 1)
    ci, cj = ai, fj
    ei, ej = di, bj

    (1 + ai + di, aj + bj + cj)
end

# Jacobian matrix
function JF(u::Vector{U}, v::Vector{V}, w::Vector{W}) where {U,V,W}
    R = promote_type(U,V, W)
    M = zeros(R, JF_size(u, v, w)...)
    JF!(M, u, v, w)
    M
end

function JF!(M, h,  u::Vector{T}, v, w) where {T}

    du, dv, dw = _degree(u), _degree(v), _degree(w)
    m, n = du + dv, du + dw

    # JF_size should return these
    r11,c11 = convmtx_size(u, dv+1)
    r13,c13 = convmtx_size(v, du+1)
    r22,c22 = convmtx_size(u, dw+1)
    r23,c23 = convmtx_size(w, du+1)

    J11 = view(M, 1:r11, 1:c11)
    J13 = view(M, 1:r13, c11 + c22 .+ (1:c23))
    J22 = view(M, r11 .+ (1:r22), c11 .+ (1:c22))
    J23 = view(M, r13 .+ (1:r23), (c11 + c22) .+ (1:c23))

    convmtx!(J11, u, dv+1)
    convmtx!(J13, v, du+1)
    convmtx!(J22, u, dw+1)
    convmtx!(J23, w, du+1)
    M[end, end-du:end] = h'

    return nothing
end

## compute F(u,v,w) - [p, p'] = [u*v, u*w] - [p, p']
function Fmp!(b, γ, p, q, uv, uw)
    b[end] = γ
    for i in 1:1+length(p)-1
        j = i
        b[i] = uv[j] - p[j]
    end
    for i in 1+length(p):length(b)-1
        j = i - length(p)
        b[i] = uw[j] - q[j]
    end
    return nothing
end

# norm( [u ⊗ v; u ⊗ w] .- [p; q], 2)
function residual_error(p::Vector{T}, q, uv, uw) where {T}
    tot = zero(T)
    for (pᵢ, uvᵢ) in zip(p,uv)
        tot += (pᵢ-uvᵢ)^2
    end
    for (qᵢ, uwᵢ) in zip(q, uw)
        tot += (qᵢ-uwᵢ)^2
    end
    sqrt(tot)
end

## --------------------------------------------------
## utils

_degree(p::Vector) = length(p) - 1

"""
    conv(p,q)

Convolution of two vectors (or product of two polynomials).
Also `⊗`
"""
function conv(p::Vector{T}, q::Vector{T}) where {T}
    R = promote_type(T,eltype(q))
    n = length(p)-1
    m = length(q)-1
    a = zeros(R, m+n+1)

    @inbounds for i = 1:n+1
        for j = 1:m+1
            a[i+j-1] += p[i] * q[j]
        end
    end
    a
end
⊗(p::Vector{T}, q::Vector{T}) where {T} = conv(p,q)

"""
    convmtx(v, n::Int)
    convmtx!(M,v,n)
    convmtx_size(v,n)

Convolution matrix.
C = convmtx(v,n) returns the convolution matrix C for a vector v.
If q is a column vector of length n, then C*q is the same as conv(v,q).

"""

function convmtx!(C, v::AbstractVector{T}, n::Int) where T

    #   Form C as the Toeplitz matrix
    #   C = Toeplitz([v; zeros(n-1)],[v[1]; zeros(n-1));  put Toeplitz code inline

    nv = length(v)-1

    @inbounds for j = 1:n
        C[j:j+nv,j] = v
    end

    nothing

end

convmtx_size(v::AbstractVector, n) = (n + length(v) - 1, n)
function convmtx(v::AbstractVector{T}, n::Int) where {T}
    C = zeros(T, convmtx_size(v, n)...)
    convmtx!(C, v, n)
    C
end

end


text:
using vOptGeneric, MOBC, PyPlot, CPLEX

function parseSSCFLP(fname, solver = CplexSolver(CPX_PARAM_SCRIND = 0))
    line_to_vec(f) = parse.(Int, split(chomp(readline(f))))
    file = open(fname)
    nbfac, nbclient = line_to_vec(file)
    I = 1:nbfac
    J = 1:nbclient
    readline(file)
    s, f = Int[], Int[]
    for i = 1:nbfac
        si, fi = line_to_vec(file)
        push!(s, si)
        push!(f, fi)
    end
    readline(file)
    d = line_to_vec(file)
    readline(file)
    c = readdlm(file, Int)
    close(file)

    m = vModel(solver=solver)
    @variable(m, x[1:nbfac, 1:nbclient], Bin)
    @variable(m, y[1:nbfac], Bin)
    @addobjective(m, Min, vecdot(c, x))
    @addobjective(m, Min, dot(f, y))
    @constraint(m, [j=J], sum(x[i,j] for i ∈ I) == 1)
    @constraint(m, [i=I], sum(d[j]*x[i,j] for j ∈ J) <= s[i]*y[i])

    m
end

for folder in ["size_5_10", "size_10_20"]#, "size_15_30"]
    cd(joinpath(Pkg.dir("MOBC"), "test", "Instances_SSCFLP", folder)) do
        for file in readdir()
            m = parseSSCFLP(file)
            println(file)
            @time solve_stidsen(m, use_nsga=true, global_branch=false, docovercuts=true, lift_covers=true)
            @time solve_parragh(m, use_nsga=true, global_branch=false, docovercuts=true, lift_covers=true)
            @time solve(m, method=:epsilon, round_results=true, verbose=false, suppress_warnings=true)
        end
    end
end


text:
include("constants.jl")
include("database.jl")
include("http.jl")

if !isinteractive()
    println("spider is running...")

    database = SQLiteInitialize()
    config = JSON.parsefile("config.json")
    const ENTRY = get(config, "entry", "")
    const TOKEN = get(config, "token", "")
    request("https://$HOST/users/$ENTRY", TypeUserinfo, TOKEN, database)
    while true
        logins = SQLiteGet(database)
        for login in logins
            request("https://$HOST/users/$login/followers", TypeFollowing, TOKEN, database)
            request("https://$HOST/users/$login/following", TypeFollowers, TOKEN, database)
            sleep(5)
        end
    end
end

max_stars_repo_name: tosone/julia-spider

max_issues_repo_name: tosone/julia-spider

max_forks_repo_name: tosone/julia-spider


text:
#=

likelihood functions tribe model

Ignacio Quintero Mächler

t(-_-t)

May 15 2017

=#

"""
    Eδx(μ::Float64, ωx::Float64, δt::Float64)

Return the expected value given the weighted average with sympatric species.
"""
Eδx(μ::Float64, ωx::Float64, δt::Float64) = ωx * μ * δt

"""
    f_λ(λ::Float64, ω::Float64, δx::Float64)

Estimate rates for area colonization based
on the difference between lineage traits and area averages.
"""
function f_λ(λ::Float64, ω::Float64, δx::Float64)
  if iszero(δx)
    return λ
  else
    return λ * (1.0 + exp(-δx))^ω
  end
end

"""
    makellf(δt::Vector{Float64}, Y::Array{Int64, 3}, ntip::Int64, narea::Int64)

Make likelihood and likelihood ratio functions
for all trait matrix and biogeography history.
"""
function makellf(δt   ::Array{Float64,1},
                 Y    ::Array{Int64,3},
                 ntip ::Int64,
                 narea::Int64,
                 m    ::Int64,
                 nedge::Int64)

  # get initial range
  wf23 = Int64[]
  for j = Base.OneTo(ntip)
    push!(wf23, findfirst(Y[:,j,1] .!= 23))
  end

  # number of normal evaluations
  n = 0
  for i = wf23
    n += length(i:(m-1))
  end

  # normal constant
  normC = -0.5*log(2.0π)*n

  function llf(X     ::Array{Float64,2},
               Y     ::Array{Int64,3},
               LA    ::Array{Float64,2},
               LD    ::Array{Float64,3},
               ωx    ::Float64,
               ω1    ::Float64,
               ω0    ::Float64,
               λ1    ::Float64,
               λ0    ::Float64,
               stemev::Vector{Vector{Float64}},
               brs   ::Array{Int64,3},
               σ²    ::Float64)

    ll::Float64 = normC

    @inbounds begin

      # trait likelihood
      for j = Base.OneTo(ntip)
        @simd for i = wf23[j]:(m-1)
          ll += logdnorm_tc(X[(i+1),j],
                            X[i,j] + Eδx(LA[i,j], ωx, δt[i]),
                            δt[i]*σ²)::Float64
        end
      end

      # biogeograhic likelihood
      for k = Base.OneTo(narea)
        ll += brll(stemev[k], λ1, λ0, brs[nedge,1,k])::Float64
        for j = Base.OneTo(ntip)
          ll += bitvectorll(Y, λ1, λ0, ω1, ω0, LD, δt,
                            j, k, wf23[j], m)::Float64
        end
      end

    end

    return ll::Float64
  end

  function llrf(Xc     ::Array{Float64,2},
                Xp     ::Array{Float64,2},
                Yc     ::Array{Int64,3},
                Yp     ::Array{Int64,3},
                LAc    ::Array{Float64,2},
                LAp    ::Array{Float64,2},
                LDc    ::Array{Float64,3},
                LDp    ::Array{Float64,3},
                ωx     ::Float64,
                ω1     ::Float64,
                ω0     ::Float64,
                λ1     ::Float64,
                λ0     ::Float64,
                stemevc::Vector{Vector{Float64}},
                stemevp::Vector{Vector{Float64}},
                brs    ::Array{Int64,3},
                brsp   ::Array{Int64,3},
                σ²     ::Float64)

    llr::Float64 = 0.0

    @inbounds begin

      # trait likelihood
      for j = Base.OneTo(ntip)
        @simd for i = wf23[j]:(m-1)
          llr += llrdnorm_xμ(Xp[(i+1),j], Xc[(i+1),j],
                             Xp[i,j] + Eδx(LAp[i,j], ωx, δt[i]),
                             Xc[i,j] + Eδx(LAc[i,j], ωx, δt[i]),
                             δt[i]*σ²)::Float64
        end
      end

      # biogeograhic likelihood
      for k = Base.OneTo(narea)
        llr += brll(stemevp[k], λ1, λ0, brsp[nedge,1,k]) -
               brll(stemevc[k], λ1, λ0,  brs[nedge,1,k])
        for j = Base.OneTo(ntip)
          llr += bitvectorll(Yp, λ1, λ0, ω1, ω0, LDp, δt, j, k, wf23[j], m) -
                 bitvectorll(Yc, λ1, λ0, ω1, ω0, LDc, δt, j, k, wf23[j], m)
        end
      end

    end

    return llr::Float64
  end

  return llf, llrf
end

"""
    llr_bm(Xc ::Array{Float64,2},
           Xp ::Array{Float64,2},
           δt ::Array{Float64,1},
           σ² ::Float64,
           idx::UnitRange)

Estimate the proposal probability of a given path according to Brownian Motion.
"""
function llr_bm(Xc ::Array{Float64,2},
                Xp ::Array{Float64,2},
                idx::UnitRange,
                δt ::Array{Float64,1},
                σ² ::Float64)

  @inbounds begin
    llr::Float64 = 0.0

    @simd for i in Base.OneTo(length(idx)-1)
      llr += llrdnorm_xμ(Xc[idx[i+1]], Xp[idx[i+1]],
                         Xc[idx[i]],   Xp[idx[i]],
                         (δt[i+1] - δt[i])*σ²)
    end
  end

  return llr::Float64
end

"""
    bitvectorll(y ::Array{Int64,1},
                λ1::Float64,
                λ0::Float64,
                ω1::Float64,
                ω0::Float64,
                Δx::Array{Float64,1},
                δt::Array{Float64,1})

Return the likelihood for a
bit vector (composed of 0s and 1s).
"""
function bitvectorll(Y ::Array{Int64,3},
                     λ1::Float64,
                     λ0::Float64,
                     ω1::Float64,
                     ω0::Float64,
                     LD::Array{Float64,3},
                     δt::Array{Float64,1},
                     j ::Int64,
                     k ::Int64,
                     yi::Int64,
                     m ::Int64)

  ll::Float64 = 0.0
  i ::Int64   = yi

  @inbounds begin

    if iszero(Y[yi,j,k])

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          ll += nell(δt[i], f_λ(λ1, ω1, LD[i,j,k]))::Float64
          i += 1
        end
        i == m && break

        ll += evll(δt[i], f_λ(λ1, ω1, LD[i,j,k]))::Float64
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          ll += nell(δt[i], f_λ(λ0, ω0, LD[i,j,k]))::Float64
          i += 1
        end
        i == m && break

        ll += evll(δt[i], f_λ(λ0, ω0, LD[i,j,k]))::Float64
        i += 1
      end

    else

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          ll += nell(δt[i], f_λ(λ0, ω0, LD[i,j,k]))::Float64
          i += 1
        end
        i == m && break

        ll += evll(δt[i], f_λ(λ0, ω0, LD[i,j,k]))::Float64
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          ll += nell(δt[i], f_λ(λ1, ω1, LD[i,j,k]))::Float64
          i += 1
        end
        i == m && break

        ll += evll(δt[i], f_λ(λ1, ω1, LD[i,j,k]))::Float64
        i += 1
      end
    end

  end

  return ll::Float64
end

"""
    bitvectorllr_ω1(Y  ::Array{Int64,3},
                    λ0 ::Float64,
                    ω0c::Float64,
                    ω0p::Float64,
                    LD ::Array{Float64,3},
                    δt ::Array{Float64,1},
                    j  ::Int64,
                    k  ::Int64,
                    yi ::Int64,
                    m  ::Int64)

Return the likelihood ratio for a
bit vector when updating `ω1`.
"""
function bitvectorllr_ω1(Y  ::Array{Int64,3},
                         λ1 ::Float64,
                         ω1c::Float64,
                         ω1p::Float64,
                         LD ::Array{Float64,3},
                         δt ::Array{Float64,1},
                         j  ::Int64,
                         k  ::Int64,
                         yi ::Int64,
                         m  ::Int64)

  llr::Float64 = 0.0
  i  ::Int64   = yi

  @inbounds begin

    if iszero(Y[yi,j,k])

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_ω(δt[i], ω1c, ω1p, λ1, LD[i,j,k])
          i += 1
        end
        i >= m && break

        llr += evllr_ω(δt[i], ω1c, ω1p, λ1, LD[i,j,k])::Float64
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1
      end

    else

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_ω(δt[i], ω1c, ω1p, λ1, LD[i,j,k])::Float64
          i += 1
        end
        i >= m && break

        llr += evllr_ω(δt[i], ω1c, ω1p, λ1, LD[i,j,k])::Float64
        i += 1
      end
    end

  end

  return llr::Float64
end

"""
    bitvectorllr_ω0(Y  ::Array{Int64,3},
                    λ0 ::Float64,
                    ω0c::Float64,
                    ω0p::Float64,
                    LD ::Array{Float64,3},
                    δt ::Array{Float64,1},
                    j  ::Int64,
                    k  ::Int64,
                    yi ::Int64,
                    m  ::Int64)

Return the likelihood ratio for a
bit vector when updating `ω0`.
"""
function bitvectorllr_ω0(Y  ::Array{Int64,3},
                         λ0 ::Float64,
                         ω0c::Float64,
                         ω0p::Float64,
                         LD ::Array{Float64,3},
                         δt ::Array{Float64,1},
                         j  ::Int64,
                         k  ::Int64,
                         yi ::Int64,
                         m  ::Int64)

  llr::Float64 = 0.0
  i  ::Int64   = yi

  @inbounds begin

    if iszero(Y[yi,j,k])

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end

        i >= m && break
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_ω(δt[i], ω0c, ω0p, λ0, LD[i,j,k])
          i += 1
        end
        i >= m && break

        llr += evllr_ω(δt[i], ω0c, ω0p, λ0, LD[i,j,k])::Float64
        i += 1
      end

    else

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_ω(δt[i], ω0c, ω0p, λ0, LD[i,j,k])
          i += 1
        end
        i >= m && break

        llr += evllr_ω(δt[i], ω0c, ω0p, λ0, LD[i,j,k])::Float64
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end

        i >= m && break
        i += 1
      end
    end

  end

  return llr::Float64
end

"""
    bitvectorllr_λ1(Y  ::Array{Int64,3},
                    λ0 ::Float64,
                    ω0c::Float64,
                    ω0p::Float64,
                    LD ::Array{Float64,3},
                    δt ::Array{Float64,1},
                    j  ::Int64,
                    k  ::Int64,
                    yi ::Int64,
                    m  ::Int64)

Return the likelihood ratio for a
bit vector when updating `ω1`.
"""
function bitvectorllr_λ1(Y  ::Array{Int64,3},
                         λ1c::Float64,
                         λ1p::Float64,
                         ω1 ::Float64,
                         LD ::Array{Float64,3},
                         δt ::Array{Float64,1},
                         j  ::Int64,
                         k  ::Int64,
                         yi ::Int64,
                         m  ::Int64)

  llr = 0.0
  i   = yi

  @inbounds begin

    if iszero(Y[yi,j,k])

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_λ(δt[i], ω1, λ1c, λ1p, LD[i,j,k])
          i   += 1
        end
        i >= m && break

        llr += evllr_λ(δt[i], ω1, λ1c, λ1p, LD[i,j,k])::Float64
        i   += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1
      end

    else

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_λ(δt[i], ω1, λ1c, λ1p, LD[i,j,k])
          i += 1
        end
        i >= m && break

        llr += evllr_λ(δt[i], ω1, λ1c, λ1p, LD[i,j,k])::Float64
        i   += 1
      end
    end

  end

  return llr::Float64
end

"""
    bitvectorllr_λ0(Y  ::Array{Int64,3},
                    λ0c::Float64,
                    λ0p::Float64,
                    ω0 ::Float64,
                    LD ::Array{Float64,3},
                    δt ::Array{Float64,1},
                    j  ::Int64,
                    k  ::Int64,
                    yi ::Int64,
                    m  ::Int64)

Return the likelihood ratio for a
bit vector when updating `ω1`.
"""
function bitvectorllr_λ0(Y  ::Array{Int64,3},
                         λ0c::Float64,
                         λ0p::Float64,
                         ω0 ::Float64,
                         LD ::Array{Float64,3},
                         δt ::Array{Float64,1},
                         j  ::Int64,
                         k  ::Int64,
                         yi ::Int64,
                         m  ::Int64)

  llr::Float64 = 0.0
  i ::Int64   = yi

  @inbounds begin

    if iszero(Y[yi,j,k])

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_λ(δt[i], ω0, λ0c, λ0p, LD[i,j,k])
          i += 1
        end
        i >= m && break
        llr += evllr_λ(δt[i], ω0, λ0c, λ0p, LD[i,j,k])::Float64
        i += 1
      end

    else

      while i < m
        while i < m && Y[i,j,k] == Y[i+1,j,k]
          llr += nellr_λ(δt[i], ω0, λ0c, λ0p, LD[i,j,k])
          i += 1
        end
        i >= m && break
        llr += evllr_λ(δt[i], ω0, λ0c, λ0p, LD[i,j,k])::Float64
        i += 1

        while i < m && Y[i,j,k] == Y[i+1,j,k]
          i += 1
        end
        i >= m && break
        i += 1
      end
    end

  end

  return llr::Float64
end

"""
    bitbitll(y1::Int64, y2::Int64, λ1::Float64, λ0::Float64, ω1::Float64, ω0::Float64, Δx::Float64, δt::Float64)

Return the likelihood for a
bit vector of length 2 (composed of 0s and 1s).
"""
function bitbitll(y1::Int64,
                  y2::Int64,
                  λ1::Float64,
                  λ0::Float64,
                  ω1::Float64,
                  ω0::Float64,
                  Δx::Float64,
                  δt::Float64)

  # event or non-event
  if y1 == y2
    return nell(δt, y1 == 0 ? f_λ(λ1, ω1, Δx) : f_λ(λ0, ω0, Δx))::Float64
  else
    return evll(δt, y1 == 0 ? f_λ(λ1, ω1, Δx) : f_λ(λ0, ω0, Δx))::Float64
  end
end

"""
    makellr_biogeo(Y    ::Array{Int64,3},
                   δt   ::Vector{Float64},
                   narea::Int64,
                   ntip ::Int64,
                   m    ::Int64)

Make likelihood function for when updating ω1 & ω0.
"""
function makellr_biogeo(Y    ::Array{Int64,3},
                        δt   ::Vector{Float64},
                        narea::Int64,
                        ntip ::Int64,
                        m    ::Int64,
                        nedge::Int64)

  # which is 23 (23 = NaN) in each column
  wf23 = Int64[]
  for j = Base.OneTo(ntip)
    push!(wf23, findfirst(Y[:,j,1] .!= 23))
  end

  function fω1(Y  ::Array{Int64,3},
               λ1 ::Float64,
               ω1c::Float64,
               ω1p::Float64,
               LD ::Array{Float64,3})

    llr::Float64 = 0.0

    @inbounds begin

      for k = Base.OneTo(narea), j = Base.OneTo(ntip)
        llr +=  bitvectorllr_ω1(Y, λ1, ω1c, ω1p, LD, δt,
                                j, k, wf23[j], m)
      end
    end

    return llr
  end

  function fω0(Y  ::Array{Int64,3},
               λ0 ::Float64,
               ω0c::Float64,
               ω0p::Float64,
               LD ::Array{Float64,3})

    llr::Float64 = 0.0

    @inbounds begin

      for k = Base.OneTo(narea), j = Base.OneTo(ntip)
        llr += bitvectorllr_ω0(Y, λ0, ω0c, ω0p, LD, δt,
                               j, k, wf23[j], m)
      end
    end

    return llr
  end

  function fλ1(Y      ::Array{Int64,3},
               λ1c    ::Float64,
               λ1p    ::Float64,
               λ0     ::Float64,
               ω1     ::Float64,
               LD     ::Array{Float64,3},
               stemevc::Vector{Vector{Float64}},
               brs    ::Array{Int64,3})

    llr::Float64 = 0.0

    @inbounds begin

      for k = Base.OneTo(narea)
        llr += brll(stemevc[k], λ1p, λ0, brs[nedge,1,k])::Float64 -
               brll(stemevc[k], λ1c, λ0, brs[nedge,1,k])::Float64
        for j = Base.OneTo(ntip)
          llr += bitvectorllr_λ1(Y, λ1c, λ1p, ω1, LD, δt,
                                 j, k, wf23[j], m)
        end
      end
    end

    return llr
  end

  function fλ0(Y      ::Array{Int64,3},
               λ1     ::Float64,
               λ0c    ::Float64,
               λ0p    ::Float64,
               ω0     ::Float64,
               LD     ::Array{Float64,3},
               stemevc::Vector{Vector{Float64}},
               brs    ::Array{Int64,3})

    llr::Float64 = 0.0

    @inbounds begin

      for k = Base.OneTo(narea)
        llr += brll(stemevc[k], λ1, λ0p, brs[nedge,1,k])::Float64 -
               brll(stemevc[k], λ1, λ0c, brs[nedge,1,k])::Float64
        for j = Base.OneTo(ntip)
          llr += bitvectorllr_λ0(Y, λ0c, λ0p, ω0, LD, δt,
                                 j, k, wf23[j], m)
        end
      end
    end

    return llr
  end

  return fω1, fω0, fλ1, fλ0
end

"""
    makellf_bgiid(bridx_a::Array{Array{UnitRange{Int64},1},1},
                  δt     ::Array{Float64,1},
                  narea  ::Int64,
                  nedge  ::Int64,
                  m      ::Int64)

Make triad likelihood function for the mutual
independence model (iid), the proposal density
for data augmented biogeographic histories.
"""
function makellf_bgiid(bridx_a::Array{Array{UnitRange{Int64},1},1},
                       δt     ::Array{Float64,1},
                       narea  ::Int64,
                       nedge  ::Int64,
                       m      ::Int64)

  # prepare δts
  δtA = Array{Float64,1}[]

  for j=bridx_a[1][1:(nedge-1)]
    inds = zeros(Int64,length(j) - 1)
    for i = eachindex(inds)
      inds[i] = rowind(j[i], m)
    end
    push!(δtA, δt[inds])
  end

  function fiid(Y      ::Array{Int64,3},
                stemev::Array{Array{Float64,1},1},
                brs    ::Array{Int64,3},
                triad  ::Array{Int64,1},
                λϕ1    ::Float64,
                λϕ0    ::Float64)

    ll::Float64 = 0.0

    @inbounds begin

      pr, d1, d2 = triad::Array{Int64,1}

      if pr < nedge
        for k = Base.OneTo(narea)
          ll += bitvectorll_iid(Y, bridx_a[k][pr], λϕ1, λϕ0, δtA[pr]) +
                bitvectorll_iid(Y, bridx_a[k][d1], λϕ1, λϕ0, δtA[d1]) +
                bitvectorll_iid(Y, bridx_a[k][d2], λϕ1, λϕ0, δtA[d2])
        end
      else
        for k = Base.OneTo(narea)
          ll += bitvectorll_iid(Y, bridx_a[k][d1], λϕ1, λϕ0, δtA[d1]) +
                bitvectorll_iid(Y, bridx_a[k][d2], λϕ1, λϕ0, δtA[d2]) +
                brll(stemev[k], λϕ1, λϕ0, brs[nedge,1,k])
        end
      end

    end

    return ll::Float64
  end

  function fiidbr(Y      ::Array{Int64,3},
             stemevc::Array{Array{Float64,1},1},
             brs    ::Array{Int64,3},
             br     ::Int64,
             λϕ1    ::Float64,
             λϕ0    ::Float64)

    ll::Float64 = 0.0

    @inbounds begin
      if br == nedge
        for k = Base.OneTo(narea)
          ll += brll(stemevc[k], λϕ1, λϕ0, brs[nedge,1,k])
        end
      else
        for k = Base.OneTo(narea)
          ll += bitvectorll_iid(Y, bridx_a[k][br], λϕ1, λϕ0,  δtA[br])
        end
      end
    end

    return ll::Float64
  end

  return fiid, fiidbr
end

"""
    stem_llr(λ1     ::Float64,
             λ0     ::Float64,
             stemc  ::Array{Int64,1},
             stemp  ::Array{Int64,1},
             stemevc::Array{Array{Float64,1},1},
             stemevp::Array{Array{Float64,1},1},
             narea  ::Int64)

Estimate likelihood ratio for stem branch.
"""
function stem_llr(λ1     ::Float64,
                  λ0     ::Float64,
                  brs    ::Array{Int64,3},
                  brsp   ::Array{Int64,3},
                  stemevc::Array{Array{Float64,1},1},
                  stemevp::Array{Array{Float64,1},1},
                  narea  ::Int64,
                  nedge  ::Int64)

  @inbounds begin

    ll::Float64 = 0.0
    for k in Base.OneTo(narea)
      ll += brll(stemevp[k], λ1, λ0, brsp[nedge,1,k]) -
            brll(stemevc[k], λ1, λ0, brs[ nedge,1,k])
    end
  end

  return ll::Float64
end

"""
    stemiid_propr(λϕ1    ::Float64,
                  λϕ0    ::Float64,
                  stemc  ::Array{Int64,1},
                  stemp  ::Array{Int64,1},
                  stemevc::Array{Array{Float64,1},1},
                  stemevp::Array{Array{Float64,1},1},
                  narea  ::Int64)

Proposal ratio for stem.
"""
function stemiid_propr(λϕ1    ::Float64,
                       λϕ0    ::Float64,
                       brs    ::Array{Int64,3,},
                       brsp   ::Array{Int64,3},
                       stemevc::Array{Array{Float64,1},1},
                       stemevp::Array{Array{Float64,1},1},
                       narea  ::Int64,
                       nedge  ::Int64)

  @inbounds begin

    ll::Float64 = 0.0
    for k in Base.OneTo(narea)
      ll += brll(stemevc[k], λϕ1, λϕ0, brs[ nedge,1,k]) -
            brll(stemevp[k], λϕ1, λϕ0, brsp[nedge,1,k])
    end
  end

  return ll::Float64
end

"""
    bitvectorll_iid(y::Array{Int64,1}, λ1::Float64, λ0::Float64, δt::Array{Float64,1})

Return likelihood under the independence model
for a bit vector.
"""
function bitvectorll_iid(Y  ::Array{Int64,3},
                         idx::UnitRange{Int64},
                         λ1 ::Float64,
                         λ0 ::Float64,
                         δt ::Array{Float64,1})

  @inbounds begin

    ll::Float64  = 0.0

    cur_s::Int64   = Y[idx[1]]
    cur_λ::Float64 = cur_s == 0 ? λ1 : λ0

    for i = Base.OneTo(length(δt))
      if Y[idx[i]] == Y[idx[i+1]]
        ll += nell(δt[i], cur_λ)::Float64
      else
        ll += evll(δt[i], cur_λ)::Float64
        cur_s = 1 - cur_s
        cur_λ = cur_s == 0 ? λ1 : λ0
      end
    end

  end

  return ll::Float64
end

"""
    makellr_ωxupd(δt::Vector{Float64}, Y::Array{Int64, 3}, ntip::Int64)

Make likelihood ratio function when updating `ωx`.
"""
function makellr_ωxσupd(δt  ::Vector{Float64},
                        Y   ::Array{Int64,3},
                        ntip::Int64)

  # which is not 23 (i.e., NaN) in each column
  w23 = UnitRange{Int64}[]
  for i = Base.OneTo(ntip)
    non23 = findall(!isequal(23), Y[:,i,1])
    push!(w23, non23[1]:non23[end-1])
  end

  function fωx(X  ::Array{Float64,2},
               LAp::Array{Float64,2},
               LAn::Array{Float64,2},
               ωxc::Float64,
               ωxp::Float64,
               σ² ::Float64)

    llr::Float64 = 0.0

    @inbounds begin
      # trait likelihood
      if ωxp >= 0.0
        if ωxc >= 0.0
          # if ωxp >= 0.0 & ωxc >= 0.0
          for j = Base.OneTo(ntip)
            @simd for i = w23[j]
              llr += llrdnorm_ωx(X[(i+1),j], X[i,j],
                                 Eδx(LAp[i,j], ωxp, δt[i]),
                                 Eδx(LAp[i,j], ωxc, δt[i]),
                                 δt[i]*σ²)
            end
          end
        else
          # if ωxp >= 0.0 & ωxc < 0.0
          for j = Base.OneTo(ntip)
            @simd for i = w23[j]
              llr += llrdnorm_ωx(X[(i+1),j], X[i,j],
                                 Eδx(LAp[i,j], ωxp, δt[i]),
                                 Eδx(LAn[i,j], ωxc, δt[i]),
                                 δt[i]*σ²)
            end
          end
        end
      else
        if ωxc >= 0.0
          # if ωxp < 0.0 & ωxc >= 0.0
          for j = Base.OneTo(ntip)
            @simd for i = w23[j]
              llr += llrdnorm_ωx(X[(i+1),j], X[i,j],
                                 Eδx(LAn[i,j], ωxp, δt[i]),
                                 Eδx(LAp[i,j], ωxc, δt[i]),
                                 δt[i]*σ²)
            end
          end
        else
          # if ωxp < 0.0 & ωxc < 0.0
          for j = Base.OneTo(ntip)
            @simd for i = w23[j]
              llr += llrdnorm_ωx(X[(i+1),j], X[i,j],
                                 Eδx(LAn[i,j], ωxp, δt[i]),
                                 Eδx(LAn[i,j], ωxc, δt[i]),
                                 δt[i]*σ²)
            end
          end
        end
      end
    end

    return llr::Float64
  end

  function fσ(X  ::Array{Float64,2},
              LA ::Array{Float64,2},
              ωx ::Float64,
              σ²c::Float64,
              σ²p::Float64)

    llr::Float64 = 0.0

    @inbounds begin
      for j = Base.OneTo(ntip)
        @simd for i = w23[j]
          llr += llrdnorm_σ²(X[(i+1),j],
                             X[i,j] + Eδx(LA[i,j], ωx, δt[i]),
                             δt[i]*σ²p, δt[i]*σ²c)
        end
      end
    end

    return llr::Float64
  end

  return fωx, fσ
end

"""
    makellr_Xupd(δt::Vector{Float64}, narea::Int64)

Make likelihood function for an internal node update in `X`.
"""
function makellr_XRupd(δt   ::Vector{Float64},
                       narea::Int64,
                       wcol ::Array{Array{Int64,1},1})

  δt1 = δt[1]
  wci = wcol[1]

  function fx(xi  ::Int64,
              xpi ::Array{Float64,1},
              X   ::Array{Float64,2},
              lapi::Array{Float64,1},
              ldpi::Array{Float64,2},
              LA  ::Array{Float64,2},
              LD  ::Array{Float64,3},
              Y   ::Array{Int64,3},
              ωx  ::Float64,
              ω1  ::Float64,
              ω0  ::Float64,
              λ1  ::Float64,
              λ0  ::Float64,
              σ²  ::Float64)

    # normal likelihoods
    llr::Float64 = 0.0

    @inbounds begin

      # loop for parent nodes
      δxim1 = δt[xi-1]
      for j = wcol[xi-1]
        llr += llrdnorm_x(xpi[j], X[xi,j],
                          X[xi-1,j] + Eδx(LA[xi-1,j], ωx, δxim1),
                          δxim1*σ²)
      end

      # loop for daughter nodes
      δxi = δt[xi]
      for j = wcol[xi]
        llr += llrdnorm_μ(X[xi+1, j],
                          xpi[j]  + Eδx(lapi[j],  ωx, δxi),
                          X[xi,j] + Eδx(LA[xi,j], ωx, δxi),
                          δxi*σ²)

        for k = Base.OneTo(narea)
          llr += bitbitll(Y[xi,j,k], Y[xi+1,j,k],
                          λ1, λ0, ω1, ω0, ldpi[j,k], δxi)::Float64 -
                 bitbitll(Y[xi,j,k], Y[xi+1,j,k],
                          λ1, λ0, ω1, ω0, LD[xi,j,k], δxi)::Float64
        end
      end
    end

    return llr::Float64
  end

  function fr(xpi ::Array{Float64,1},
             X   ::Array{Float64,2},
             σ²  ::Float64)

    llr::Float64 = 0.0

    @inbounds begin

      # loop for daughter nodes
      for j = wci
        llr += llrdnorm_μ(X[2,j], xpi[j], X[1,j], δt1*σ²)
      end
    end

    return llr::Float64
  end

  return fx, fr
end

"""
    brll(brevs::Array{Float64,1}, λ1::Float64, λ0::Float64, si::Int64)

Return likelihood for a branch in continuous time.
"""
function brll(brevs::Array{Float64,1}, λ1::Float64, λ0::Float64, si::Int64)

  ll::Float64 = 0.0

  if lastindex(brevs) > 1
    for i = Base.OneTo(lastindex(brevs)-1)
      ll += evll(brevs[i], iszero(si) ? λ1 : λ0)::Float64
      si  = 1 - si
    end
  end

  ll += nell(brevs[end], iszero(si) ? λ1 : λ0)::Float64

  return ll::Float64
end

"""
    evll(t::Float64, λ::Float64)

Return log-likelihood for events.
"""
evll(t::Float64, λ::Float64) = (log(λ) - (λ * t))::Float64

"""
    evllr_ω(t  ::Float64,
             ω1c::Float64,
             ω1p::Float64,
             λ1 ::Float64,
             δx ::Float64)

Return log-likelihood ratio for events when updating `ω`.
"""
function evllr_ω(t  ::Float64,
                 ωc::Float64,
                 ωp::Float64,
                 λ ::Float64,
                 δx ::Float64)
  if iszero(δx)
    return 0.0
  else
    eϕp1 = 1.0 + exp(-δx)
    return (ωp - ωc) * log(eϕp1) +
            t * λ * (eϕp1^ωc - eϕp1^ωp)
  end
end

"""
    evllr_λ(t  ::Float64,
            ω::Float64,
            λc::Float64,
            λp::Float64,
            δx ::Float64)

Return log-likelihood ratio for non-events when updating `λ`.
"""
function evllr_λ(t ::Float64,
                 ω ::Float64,
                 λc::Float64,
                 λp::Float64,
                 δx::Float64)
  if iszero(δx)
    return log(λp/λc) +
           t * (λc - λp)
  else
    return log(λp/λc) +
           t * (1.0 + exp(-δx))^ω * (λc - λp)
  end
end

"""
    nell(t::Float64, λ::Float64)

Return log-likelihood for nonevents.
"""
nell(t::Float64, λ::Float64) = -λ*t

"""
    nellr_ω(t ::Float64,
            ωc::Float64,
            ωp::Float64,
            λ ::Float64,
            δx::Float64)

Return log-likelihood ratio for non-events when updating `ω`.
"""
function nellr_ω(t ::Float64,
                 ωc::Float64,
                 ωp::Float64,
                 λ ::Float64,
                 δx::Float64)
  if iszero(δx)
    return 0.0
  else
    eϕp1 = 1.0 + exp(-δx)
    return t * λ * (eϕp1^ωc - eϕp1^ωp)
  end
end

"""
    nellr_λ(t  ::Float64,
             ω::Float64,
             λc::Float64,
             λp::Float64,
             δx ::Float64)

Return log-likelihood ratio for non-events when updating `λ`.
"""
function nellr_λ(t  ::Float64,
                  ω::Float64,
                  λc::Float64,
                  λp::Float64,
                  δx ::Float64)
  if iszero(δx)
    return t * (λc - λp)
  else
    return t * (1.0 + exp(-δx))^ω * (λc - λp)
  end
end

"""
    allλpr(λc::Array{Float64,2}, λprior::Float64)

Return log-prior for all areas
"""
function allλpr(λ1    ::Float64,
                λ0    ::Float64,
                λprior::Float64)
  2.0*log(λprior) - λprior * (λ1 + λ0)
end

max_stars_repo_name: hmorlon/JPANDA

max_issues_repo_name: hmorlon/JPANDA

max_forks_repo_name: hmorlon/JPANDA


text:
module BulletsMaggots
using Random, ProgressMeter

include("tools.jl")
include("solver.jl")
include("build_table.jl")
include("play.jl")
include("jonah_solver.jl")

end # module


text:
using Test
using SemanticModels

tests = [
         "epidemics.jl",
         "petricospans.jl",
         ]

for test in tests
  @testset "Running $test" begin
    include(test)
  end
end


text:
struct PoissBinAliasSampler <: Sampleable{Univariate,Discrete}
    table::AliasTable
end

PoissBinAliasSampler(d::PoissonBinomial) = PoissBinAliasSampler(AliasTable(d.pmf))

rand(rng::AbstractRNG, s::PoissBinAliasSampler) = rand(rng, s.table) - 1


text:
using Pkg

using ContinuousBenchmarks
using ContinuousBenchmarks.Runner
using ContinuousBenchmarks.Config

project_root = (@__DIR__) |> dirname
Config.set_config("target.project_dir", project_root)

for bm_file in get_benchmark_files()
    Runner.run_benchmarks([bm_file], save_path="."; ignore_error=false)
end


text:
function get_dataset()
  [
    {"method" => "echo", "arg" => nothing},
    {"method" => "echo", "arg" => true},
    {"method" => "echo", "arg" => false},
    {"method" => "echo", "arg" => 1<<31},
    {"method" => "echo", "arg" => ""},
    {"method" => "echo", "arg" => "String"},
    {"method" => "echo", "arg" => []},
    {"method" => "echo", "arg" => {}},
    {"method" => "echo", "arg" => [10, 20, 30]},
    {"method" => "echo", "arg" => [10, "20", 30]},
    {"method" => "echo", "arg" => {10, 20, 30}},
    {"method" => "echo", "arg" => {{10, 20, 30}, {40, 50}, {60}} },
    {"method" => "echo", "arg" => {"key" => "value"}},
    {"method" => "echo", "arg" => {"k1" => "v1", "k2" => "v2"}},
    {"method" => "echo", "arg" => {"k1"=>nothing, "k2"=>false, "k3"=>"", "k4"=>{}, "k5"=>{"k5.1"=>false, "k5.2"=>{"k5.2.1"=>1<<31} } } },
  #  {"method" => "echo", "arg" => π },
    {"method" => "echo", "arg" => (), "expect" => {} },
    {"method" => "echo", "arg" => ((10,20,30,),(40,50,),(60,),), "expect" => {{10,20,30},{40,50},{60}} },
  ]
end

function get_dataset_for_arguments()
  [
    {"method" => "echo", "arg1" => 1, "arg2" => 2, "arg3" => 3, "arg4" => 4, "expect" => {1,2,3,4}}
  ]
end

max_stars_repo_name: keithseahus/msgpack-rpc-julia

max_issues_repo_name: keithseahus/msgpack-rpc-julia

max_forks_repo_name: keithseahus/msgpack-rpc-julia