text:
function combine_union(g1::AbstractGraph, g2::AbstractGraph; kwargs...)
    g3 = union(g1, g2)
    cliques = maximal_cliques(cg)
    cliquetree, nb_lc = get_cliquetree(cliques)
    return cliques, cliquetree, g3, nb_lc, Dict()
end

function combine_intersect(g1, g2; dec_type::Symbol=:cholesky, kwargs...)
    g3 = intersect(g1, g2)
    if is_chordal(g3)
        cliques, cliquetree, nb_lc = get_decomposition_data(g3)
        return cliques, cliquetree, g3, nb_lc, Dict()
    else
        return decompose(g3, dec_type)
    end
end

function combine_intersect_dist(g1, g2)
    g3 = intersect(g1, g2)
    if is_chordal(g3)
        cliques, cliquetree, nb_lc = get_decomposition_data(g3)
        return cliques, cliquetree, g3, nb_lc, Dict()
    else
        miss_edges = union(setdiff(edges(g1), edges(g3)), setdiff(edges(g2), edges(g3)))
        while !isempty(miss_edges) || is_chordal(g3)
            @info "$(length(miss_edges))"
            dists = [length(a_star(g3, e.src, e.dst)) for e in miss_edges]
            ei = sortperm(dists)[1]
            add_edge!(g3, miss_edges[ei])
            deleteat!(miss_edges, ei)
        end
    end
end

function combine_intersect_percent(g1, g2; percent=0.1, dec_type::Symbol=:cholesky, kwargs...)
    g3 = intersect(g1, g2)
    if is_chordal(g3)
        cliques, cliquetree, nb_lc = get_decomposition_data(g3)
        return cliques, cliquetree, g3, nb_lc, Dict()
    else
        miss_edges_g1 = setdiff(edges(g1), edges(g3))
        miss_edges_g2 = setdiff(edges(g2), edges(g3))
        nb_edges_to_add_g1 = floor(Int, percent * length(miss_edges_g1))
        nb_edges_to_add_g2 = floor(Int, percent * length(miss_edges_g2))
        for e in sample(miss_edges_g1, nb_edges_to_add_g1, replace=false)
            add_edge!(g3, e)
        end
        for e in sample(miss_edges_g2, nb_edges_to_add_g2, replace=false)
            add_edge!(g3, e)
        end
        return decompose(g3, dec_type)
    end
end

# TODO: use Molzahn heuristic to choose which cliques to select
#function molzahn_heuristic_combination(cliques, clique, cliquetree_graph, dstmx)
#    add_vertex!(cliquetree_graph)
#    index = nv(cliquetree_graph)
#    set_prop!(cliquetree_graph, nv(cliquetree_graph), :clique, clique)
#    dstmx = hcat(dstmx, zeros(Int, size(dstmx)[1]))
#    dstmx = vcat(dstmx, zeros(Int, size(dstmx)[2]))
#    for v in vertices(cliquetree_graph)[begin:end - 1]
#        length_intersect = length(intersect(clique, get_prop(cliquetree_graph, v, :clique)))
#        if length_intersect > 0
#            add_edge!(cliquetree_graph, v, index)
#            dstmx[v, index] = length_intersect
#            dstmx[index, v] = length_intersect
#        end
#    end
#end

#function get_neighbor(cliques, clique)
#    neighbor = cliques[1]
#    neighbor_similarity = length(intersect(cliques[1], clique))
#    for c in cliques[2:end]
#        similarity = length(intersect(c, clique))
#        if similarity > neighbor_similarity
#            neighbor = c
#            neighbor_similarity = similarity
#        end
#    end
#    return neighbor
#end

function combine_cliques_intersect(g1, g2; instance_path, dec_type=:cholesky, kwargs...)
    cliques1 = sort(map(sort, maximal_cliques(g1)))
    cliques2 = sort(map(sort, maximal_cliques(g2)))
    cliques3 = intersect(cliques1, cliques2)
    isempty(cliques3) && return nothing
    g3 = construct_graph_from_cliques(cliques3)
    g3 = union(g3, construct_network_graph(instance_path))
    return decompose(g3, dec_type)
end

const combine_map = Dict(
    :union => combine_union,
    :intersect => combine_intersect,
    :intersect_percent => combine_intersect_percent,
    :intersect_cliques => combine_cliques_intersect
)

function combine_decomposition(g1::AbstractGraph, g2::AbstractGraph; combine_type::Symbol=:intersect, kwargs...)
    return combine_map[combine_type](g1, g2; kwargs...)
end

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
##############################################################################
##
## We use SubDataFrame's to maintain a reference to a subset of a DataFrame
## without making copies.
##
##############################################################################

"""
A view of row subsets of an AbstractDataFrame

A `SubDataFrame` is meant to be constructed with `view`.  A
SubDataFrame is used frequently in split/apply sorts of operations.

```julia
view(d::AbstractDataFrame, rows)
```

### Arguments

* `d` : an AbstractDataFrame
* `rows` : any indexing type for rows, typically an Int,
  AbstractVector{Int}, AbstractVector{Bool}, or a Range

### Notes

A `SubDataFrame` is an AbstractDataFrame, so expect that most
DataFrame functions should work. Such methods include `describe`,
`dump`, `nrow`, `size`, `by`, `stack`, and `join`. Indexing is just
like a DataFrame; copies are returned.

To subset along columns, use standard column indexing as that creates
a view to the columns by default. To subset along rows and columns,
use column-based indexing with `view`.

### Examples

```julia
df = DataFrame(a = repeat([1, 2, 3, 4], outer=[2]),
               b = repeat([2, 1], outer=[4]),
               c = randn(8))
sdf1 = view(df, 1:6)
sdf2 = view(df, df[:a] .> 1)
sdf3 = view(df[[1,3]], df[:a] .> 1)  # row and column subsetting
sdf4 = groupby(df, :a)[1]  # indexing a GroupedDataFrame returns a SubDataFrame
sdf5 = view(sdf1, 1:3)
sdf1[:,[:a,:b]]
```

"""
immutable SubDataFrame{T <: AbstractVector{Int}} <: AbstractDataFrame
    parent::DataFrame
    rows::T # maps from subdf row indexes to parent row indexes

    function (::Type{SubDataFrame{T}}){T}(parent::DataFrame, rows::T)
        if length(rows) > 0
            rmin, rmax = extrema(rows)
            if rmin < 1 || rmax > size(parent, 1)
                throw(BoundsError())
            end
        end
        new{T}(parent, rows)
    end
end

function SubDataFrame{T <: AbstractVector{Int}}(parent::DataFrame, rows::T)
    return SubDataFrame{T}(parent, rows)
end

function SubDataFrame(parent::DataFrame, row::Integer)
    return SubDataFrame(parent, [row])
end

function SubDataFrame{S <: Integer}(parent::DataFrame, rows::AbstractVector{S})
    return view(parent, Int(rows))
end

function Base.view{S <: Real}(df::DataFrame, rowinds::AbstractVector{S})
    return SubDataFrame(df, rowinds)
end

function Base.view{S <: Real}(sdf::SubDataFrame, rowinds::AbstractVector{S})
    return SubDataFrame(sdf.parent, sdf.rows[rowinds])
end

function Base.view(df::DataFrame, rowinds::AbstractVector{Bool})
    return view(df, getindex(SimpleIndex(size(df, 1)), rowinds))
end

function Base.view(sdf::SubDataFrame, rowinds::AbstractVector{Bool})
    return view(sdf, getindex(SimpleIndex(size(sdf, 1)), rowinds))
end

function Base.view(adf::AbstractDataFrame, rowinds::Integer)
    return SubDataFrame(adf, Int[rowinds])
end

function Base.view(adf::AbstractDataFrame, rowinds::Any)
    return view(adf, getindex(SimpleIndex(size(adf, 1)), rowinds))
end

function Base.view(adf::AbstractDataFrame, rowinds::Any, colinds::Any)
    return view(adf[[colinds]], rowinds)
end

##############################################################################
##
## AbstractDataFrame interface
##
##############################################################################

index(sdf::SubDataFrame) = index(sdf.parent)

# TODO: Remove these
nrow(sdf::SubDataFrame) = ncol(sdf) > 0 ? length(sdf.rows)::Int : 0
ncol(sdf::SubDataFrame) = length(index(sdf))

function Base.getindex(sdf::SubDataFrame, colinds::Any)
    return sdf.parent[sdf.rows, colinds]
end

function Base.getindex(sdf::SubDataFrame, rowinds::Any, colinds::Any)
    return sdf.parent[sdf.rows[rowinds], colinds]
end

function Base.setindex!(sdf::SubDataFrame, val::Any, colinds::Any)
    sdf.parent[sdf.rows, colinds] = val
    return sdf
end

function Base.setindex!(sdf::SubDataFrame, val::Any, rowinds::Any, colinds::Any)
    sdf.parent[sdf.rows[rowinds], colinds] = val
    return sdf
end

##############################################################################
##
## Miscellaneous
##
##############################################################################

Base.map(f::Function, sdf::SubDataFrame) = f(sdf) # TODO: deprecate

function Base.delete!(sdf::SubDataFrame, c::Any) # TODO: deprecate?
    return SubDataFrame(delete!(sdf.parent, c), sdf.rows)
end

without(sdf::SubDataFrame, c::Vector{Int}) = view(without(sdf.parent, c), sdf.rows)
without(sdf::SubDataFrame, c::Int) = view(without(sdf.parent, c), sdf.rows)
without(sdf::SubDataFrame, c::Any) = view(without(sdf.parent, c), sdf.rows)

ext: jl


text:
module SOH
include("boundary_conditions.jl")
include("flux.jl")
include("init.jl")
include("plot_save.jl")
include("run.jl")
include("scheme.jl")
include("toolbox.jl")

export boundary_conditions_xy!, boundary_conditions_ρuv!
export flux_x
export random_init, quadratic_potential_force, flat_quadratic_potential_force
export plot_rhoUV, update_plot!, save_data!, radial_density
export run!
export scheme_iter!
export make_new_dir, coefficients_Vicsek, nice_float2string

end     #module


text:
using UnivariateFunctions
using StochasticIntegrals
using Dates
using LinearAlgebra
using Distributions: Normal, MersenneTwister, quantile
using Statistics: var, mean, cov
using Sobol
using Random
const tol = 10*eps()

USD_mean_reversion  = 0.15
USD_vol = 0.01
GBP_mean_reversion = 0.1
GBP_vol = 0.001
GBP_FX_Vol = 0.2

today    = Date(2016,12,1)
tommorow = Date(2016,12,2)

Barclays_vol = PE_Function(0.1, 0.0, 0.0, 0)
CS_vol       = PE_Function(0.2, 0.0, 0.0, 0)

brownian_corr_matrix = Symmetric([1.0 0.75;
                                  0.75 1.0])
brownian_ids = [ :CS, :BARC]
BARC_ito    = ItoIntegral(:BARC, Barclays_vol)
CS_ito    = ItoIntegral(:CS, CS_vol)
ito_integrals = Dict([:CS, :BARC] .=> [CS_ito, BARC_ito])

ito_set_ = ItoSet(brownian_corr_matrix, brownian_ids, ito_integrals)
covar = CovarianceAtDate(ito_set_, years_from_global_base(today), years_from_global_base(tommorow))

T = Float64
conditioning_draws = Dict{Symbol,T}(:CS => 0.02)
conditional_mean, conditional_covar = generate_conditioned_distribution(covar, conditioning_draws)
conditional_mean[1] > 100*eps()


text:
###Importing dependencies
using Pkg
Pkg.activate("./")
using DataFrames, MClassification, MLJ, Queryverse, Plots
@load KNNClassifier

###Loading the dataset
dataset_path = "./datasets/sinthetic/"
dataset = Queryverse.load(dataset_path * "UG_2C_5D.csv", header_exists=false) |> DataFrame

###Setting the data
X = dataset[:,1:end-1]
y = dataset[:,end]

training_set_size = 150
y_predicted = [Array{Any, 1}() for i in 1:2]

###KNN
knn_classifier = MLJ.machine(KNNClassifier(K=3), X, coerce(y, autotype(y)))
train, test = (collect(1:training_set_size), collect(training_set_size+1:length(y)))
MLJ.fit!(knn_classifier, rows=train)
y_predicted[1] = predict_mode(knn_classifier, rows=test)

###MClassification
m_classifier = MClassification.fit(X[1:150, :], y[1:150], 0.1)
y_predicted[2] = MClassification.predict(m_classifier, convert(Array{Float64, 2}, X[151:end, :]))

n_instances = 0
hates = [Array{Float64, 1}() for i in 1:2]
hits = zeros(Int64, 2)
for i in 1:length(y_predicted[2])
    global n_instances += 1

    if y[i+150] == y_predicted[1][i]
        hits[1] += 1
    end

    if y[i+150] == y_predicted[2][i]
        hits[2] += 1
    end

    if n_instances % 2000 == 0
            append!(hates[1], hits[1] / 2000)
            append!(hates[2], hits[2] / 2000)
            hits[1] = 0
            hits[2] = 0
    end
end
if(hits[1] > 0 || hits[2] > 0)
    append!(hates[1], hits[1] / (length(y_predicted[2])%2000))
    append!(hates[2], hits[2] / (length(y_predicted[2])%2000))
end
plot(hates, ylims = (0, 1), labels = permutedims(["Static-KNN", "MClassification"]), legend=:outertopright)

y_test = y[151:end]

accuracy =
if y_test[i]==
y_predicted
function Accuracy(y::Array{T, N}, y_predicted::Array{T, N}, steps=1) where {T == T}
{
    if length(y_predicted) != length(y)
        return catch
    end
    hits = 0
    for i in 1:length(y)
        if y[i] == y_predicted[i]
            hits += 1
        end
    end
}


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   This file contains functions reltated to input handling. This code was
#   adapted from the on in TextUserInterfaces.jl.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

################################################################################
#                                  Constants
################################################################################

include("keycodes.jl")

################################################################################
#                                  Functions
################################################################################

"""
    _jlgetch(stream::IO)

Wait for an keystroke in the stream `stream` and return it (see
[`Keystroke`](@ref)).

"""
function _jlgetch((@nospecialize stream::IO))
    c_raw = read(stream, UInt8)::UInt8

    c::UInt32 = UInt32(c_raw)
    nc::Int32 = 0

    if c == 27

        s = string(Char(c))

        # Read the entire sequence limited to 10 characters.
        for i = 1:10
            stream.buffer.size == i && break
            nc = read(stream, Char)::Char
            s *= string(Char(nc))
            haskey(keycodes, s) && break
        end

        if length(s) == 1
            return Keystroke(raw = s, value = :esc)
        elseif haskey(keycodes, s)
            aux = keycodes[s]
            return Keystroke(
                raw = s,
                value = aux.value,
                alt = aux.alt,
                ctrl = aux.ctrl,
                shift = aux.shift
            )
        else
            return Keystroke(raw = s, value = :undefined)
        end
    elseif c == nocharval
        return Keystroke(raw = string(c), value = :undefined)
    elseif 192 <= c <= 223 # utf8 based logic starts here
        bs1 = UInt8(c)
        bs2 = read(stream, UInt8)::UInt8
        return Keystroke(
            raw = string(bs1) * ", " * string(bs2),
            value = String([bs1, bs2])
        )
    elseif c < 192 || c > 253
        if c == 4
            return Keystroke(raw = string(c), value = :eot)
        elseif c == 9
            return Keystroke(raw = string(c), value = :tab)
        elseif c == 10
            return Keystroke(raw = string(c), value = :enter)
        elseif c == 13
            return Keystroke(raw = string(c), value = :enter)
        elseif c == 21
            return Keystroke(raw = string(c), value = :shiftin)
        elseif c == 127
            return Keystroke(raw = string(c), value = :backspace)
        elseif c == 410
            return Keystroke(raw = string(c), value = :resize)
        else
            return Keystroke(raw = string(c), value = string(Char(c)))
        end
    elseif  224 <= c <= 239
        bs1 = UInt8(c)
        bs2 = read(stream, UInt8)::UInt8
        bs3 = read(stream, UInt8)::UInt8
        return Keystroke(
            raw = string(bs1) * ", " *
                  string(bs2) * ", " *
                  string(bs3),
            value = String([bs1, bs2, bs3])
        )
    elseif  240 <= c <= 247
        bs1 = UInt8(c)
        bs2 = read(stream, UInt8)::UInt8
        bs3 = read(stream, UInt8)::UInt8
        bs4 = read(stream, UInt8)::UInt8
        return Keystroke(
            raw = string(bs1) * ", " *
                  string(bs2) * ", " *
                  string(bs3) * ", " *
                  string(bs4),
            value = String([bs1, bs2, bs3, bs4])
        )
    elseif  248 <= c <= 251
        bs1 = UInt8(c)
        bs2 = read(stream, UInt8)::UInt8
        bs3 = read(stream, UInt8)::UInt8
        bs4 = read(stream, UInt8)::UInt8
        bs5 = read(stream, UInt8)::UInt8
        return Keystroke(
            raw = string(bs1) * ", " *
                  string(bs2) * ", " *
                  string(bs3) * ", " *
                  string(bs4) * ", " *
                  string(bs5),
            value = String([bs1, bs2, bs3, bs4, bs5])
        )
    elseif  252 <= c <= 253
        bs1 = UInt8(c)
        bs2 = read(stream, UInt8)::UInt8
        bs3 = read(stream, UInt8)::UInt8
        bs4 = read(stream, UInt8)::UInt8
        bs5 = read(stream, UInt8)::UInt8
        bs6 = read(stream, UInt8)::UInt8
        return Keystroke(
            raw = string(bs1) * ", " *
                  string(bs2) * ", " *
                  string(bs3) * ", " *
                  string(bs4) * ", " *
                  string(bs5) * ", " *
                  string(bs6),
            value = String([bs1, bs2, bs3, bs4, bs5, bs6])
        )
    end

    return Keystroke(value = :undefined)
end


text:
using Soss, DifferentialEquations
using Plots, StatsPlots

function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, γ, δ  = p
  du[1] = (α - β*y)x # dx =
  du[2] = (δ*x - γ)y # dy =
end
p = [1.5, 1.0, 3.0, 1.0]
u0 = [1.0,1.0]
prob1 = ODEProblem(lotka_volterra,u0,(0.0,10.0),p)

sol1 = solve(prob1,Tsit5(),saveat=0.1)
odedata = rand(For(Poisson, Array(sol1)))

plot(sol1, legend = false); scatter!(sol1.t, odedata')

m = @model ode begin
    α ~ Uniform()
    β ~ Uniform()
    γ ~ Uniform()
    δ ~ Uniform()

    # p = [2α+0.5,2β,3γ+1,2δ]
    p = [20α+0.5,20β,30γ+1,20δ]
    prob = remake(ode, p=p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

    data ~ For(Array(predicted)) do λ Poisson(λ) end
end

post = dynamicHMC(m(ode=prob1) | (data=odedata,))

plt = plot(sol1, legend = false)

for x in post
    p = [2 * x.α + 0.5,2 * x.β,3 * x.γ + 1,2 * x.δ]

    thisprob = remake(prob1, p=p)
    thissol = solve(thisprob,Tsit5(),saveat=0.1)
    plot!(plt, thissol, alpha=0.02, linewidth=2, color = :black, legend = false)
end
# display(pl)
plot!(sol1, w=1, legend = false)
scatter!(sol1.t, odedata')


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

# TODO: radical

const SMARTS_ATOM_COND_SYMBOL = Dict(
    'X' => :connectivity,
    'D' => :nodedegree,
    'v' => :valence,
    'H' => :hydrogenconnected,
    'r' => :smallestsssr,
    'R' => :sssrcount
)

const SMARTS_CHARGE_SIGN = Dict(
    '+' => 1,
    '-' => -1
)

"""
    atomsymbol!(state::SmartsParserState) -> Union{Pair,Nothing}

Atomsymbol <- Br / Cl / [AaBCcNnOoPpSsFI*]
"""
function atomsymbol!(state::SmartsParserState)
    sym1 = read(state)
    sym2 = lookahead(state, 1)
    if sym1 == 'C' && sym2 == 'l'
        forward!(state, 2)
        return QueryFormula(
            :and, Set([QueryFormula(:atomsymbol, :Cl), QueryFormula(:isaromatic, false)]))
    elseif sym1 == 'B' && sym2 == 'r'
        forward!(state, 2)
        return QueryFormula(
            :and, Set([QueryFormula(:atomsymbol, :Br), QueryFormula(:isaromatic, false)]))
    elseif sym1 in "BCNOPSFI"
        forward!(state)
        return QueryFormula(
            :and, Set([QueryFormula(:atomsymbol, Symbol(sym1)), QueryFormula(:isaromatic, false)]))
    elseif sym1 in "bcnops"
        forward!(state)
        return QueryFormula(
            :and, Set([
                QueryFormula(:atomsymbol, Symbol(uppercase(sym1))),
                QueryFormula(:isaromatic, true)
            ]))
    elseif sym1 == 'A'
        forward!(state)
        return QueryFormula(:isaromatic, false)
    elseif sym1 == 'a'
        forward!(state)
        return QueryFormula(:isaromatic, true)
    elseif sym1 == '*'
        forward!(state)
        return QueryFormula(:any, true)
    end
    # return nothing
end

"""
    atomprop!(state::SmartsParserState) -> Union{Pair,Nothing}

AtomProp <- '\$(' RecursiveQuery ')' / Mass / Symbol / AtomNum / Stereo / CHG / [DHRrvX]
"""
function atomprop!(state::SmartsParserState)
    sym1 = read(state)
    sym2 = lookahead(state, 1)
    if haskey(ATOMSYMBOLMAP, string(uppercase(sym1), sym2))
        # Two letter atoms
        forward!(state, 2)
        if string(uppercase(sym1), sym2) in ("As", "Se")
            return QueryFormula(
                :and, Set([
                    QueryFormula(:atomsymbol, Symbol(uppercase(sym1), sym2)),
                    QueryFormula(:isaromatic, islowercase(sym1))
                ]))
        end
        @assert isuppercase(sym1)
        return QueryFormula(:atomsymbol, Symbol(sym1, sym2))
    elseif haskey(SMARTS_ATOM_COND_SYMBOL, sym1)
        # Neighbor and ring conditions
        forward!(state)
        if isdigit(sym2)
            num = parse(Int, sym2)
            forward!(state)
        else
            sym1 in ('r', 'R') && return QueryFormula(:not, QueryFormula(:sssrcount, 0))
            num = 1
        end
        return QueryFormula(SMARTS_ATOM_COND_SYMBOL[sym1], num)
    elseif sym1 in ('A', 'a', '*')
        forward!(state)
        sym1 == 'A' && return QueryFormula(:isaromatic, false)
        sym1 == 'a' && return QueryFormula(:isaromatic, true)
        sym1 == '*' && return QueryFormula(:any, true)
    elseif haskey(ATOMSYMBOLMAP, string(uppercase(sym1)))
        # Single letter atoms
        forward!(state)
        if uppercase(sym1) in "BCNOPS"
            return QueryFormula(
                :and, Set([
                    QueryFormula(:atomsymbol, Symbol(uppercase(sym1))),
                    QueryFormula(:isaromatic, islowercase(sym1))
                ]))
        end
        @assert isuppercase(sym1)
        return QueryFormula(:atomsymbol, Symbol(sym1))
    elseif sym1 == '#'
        # Atomic number
        forward!(state)
        start = state.pos
        while isdigit(lookahead(state, 1))
            forward!(state)
        end
        num = parse(Int, SubString(state.input, start, state.pos))
        forward!(state)
        return QueryFormula(:atomsymbol, atomsymbol(num))
    elseif sym1 in keys(SMARTS_CHARGE_SIGN)
        # Charge
        forward!(state)
        if isdigit(sym2)
            chg = parse(Int, sym2)
            forward!(state)
        else
            chg = 1
            while read(state) == sym1
                forward!(state)
                chg += 1
            end
        end
        return QueryFormula(:charge, chg * SMARTS_CHARGE_SIGN[sym1])
    elseif sym1 == '@'
        # Stereo
        # @ -> anticlockwise, @@ -> clockwise, ? -> or not specified
        forward!(state)
        cw = false
        if read(state) == '@'
            forward!(state)
            cw = true
        end
        if read(state) == '?'
            forward!(state)
            return QueryFormula(:not, QueryFormula(:stereo, cw ? :anticlockwise : :clockwise))
        else
            return QueryFormula(:stereo, cw ? :clockwise : :anticlockwise)
        end
    elseif isdigit(sym1)
        # Isotope
        start = state.pos
        while isdigit(lookahead(state, 1))
            forward!(state)
        end
        num = SubString(state.input, start, state.pos)
        forward!(state)
        return QueryFormula(:mass, parse(Int, num))
    elseif sym1 == '$' && sym2 == '('
        # Recursive
        forward!(state, 2)
        start = state.pos
        toclose = 1
        while true
            if read(state) == ')'
                if toclose == 1
                    break
                else
                    toclose -= 1
                end
            elseif read(state) == '('
                toclose += 1
            end
            forward!(state)
        end
        q = SubString(state.input, start, state.pos - 1)
        forward!(state)
        return QueryFormula(:recursive, q)
    end
    # return nothing
end

"""
    atom!(state::SmilesParser) -> Vector{SmilesAtom}

Atom <- '[' AtomProp+ ']' / AtomSymbol
"""
function atom!(state::SmilesParser)
    sym1 = read(state)
    if sym1 == '['
        forward!(state)
        fml = lghighand!(state, atomprop!)
        fml = tidyformula(fml)
        symcls = read(state)
        @assert symcls == ']' "(atom!) unexpected token: $(symcls)"
        forward!(state)
        @assert findformula(fml, :atomsymbol) !== nothing || findformula(fml, :hydrogenconnected) >= 1
        atoms = [SmilesAtom(
            something(findformula(fml, :atomsymbol), :H),
            something(findformula(fml, :charge), 0),
            1,  # TODO: multiplicity
            findformula(fml, :mass),
            something(findformula(fml, :isaromatic), false),
            something(findformula(fml, :stereo), :unspecified)
        )]
        hcnt = something(findformula(fml, :hydrogenconnected), 0)
        if findformula(fml, :atomsymbol) === nothing
            hcnt -= 1
        end
        for h in 1:hcnt
            push!(atoms, SmilesAtom(:H))
        end
        return atoms
    else
        fml = atomsymbol!(state)
        if fml === nothing
            return SmilesAtom[]
        else
            asym = findformula(fml, :atomsymbol)
            arom = findformula(fml, :isaromatic)
            return [SmilesAtom(asym, 0, 1, nothing, arom, :unspecified)]
        end
    end
end

"""
    atom!(state::SmartsParser) -> Vector{SmartsAtom}

Atom <- '[' (AtomProp / LogicalOperator)+ ']' / AtomSymbol
"""
function atom!(state::SmartsParser)
    sym1 = read(state)
    if sym1 == '['
        forward!(state)
        fml = lglowand!(state, atomprop!)
        @assert fml !== nothing "(atom!) empty atomprop"
        fml = tidyformula(fml)
        symcls = read(state)
        @assert symcls == ']' "(atom!) unexpected token: $(symcls)"
        forward!(state)
        return [SmartsAtom(fml)]
    else
        fml = atomsymbol!(state)
        if fml === nothing
            return SmartsAtom[]
        else
            return [SmartsAtom(fml)]
        end
    end
end


text:
module MSTParser

include("tarjan.jl")

end


text:
using Sialia
using Documenter

DocMeta.setdocmeta!(Sialia, :DocTestSetup, :(using Sialia); recursive=true)

makedocs(;
    modules=[Sialia],
    authors="Mark Wardle <mark@potient.com> and contributors",
    repo="https://github.com/wardlem/Sialia.jl/blob/{commit}{path}#{line}",
    sitename="Sialia.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://wardlem.github.io/Sialia.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/wardlem/Sialia.jl",
)


text:
# Autogenerated wrapper script for Calcium_jll for x86_64-w64-mingw32
export libcalcium

using FLINT_jll
using Arb_jll
using Antic_jll
using GMP_jll
using MPFR_jll
JLLWrappers.@generate_wrapper_header("Calcium")
JLLWrappers.@declare_library_product(libcalcium, "libcalcium-0.dll")
function __init__()
    JLLWrappers.@generate_init_header(FLINT_jll, Arb_jll, Antic_jll, GMP_jll, MPFR_jll)
    JLLWrappers.@init_library_product(
        libcalcium,
        "bin\\libcalcium-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# -----------------------------------------------------------------------
# Codes used during the experiments
# -----------------------------------------------------------------------
# Required packages
using PyCall, Dates, JLD
include("cross_entropy.jl")
include("ga_decimal.jl")

# Including path to scanivalve
push!(pyimport("sys")["path"], "path to scanivalve folder")

#---------------------------------------------------------------------------------------------------------------
# SCANIVALVE CONNECTION AND CONFIGURING
#---------------------------------------------------------------------------------------------------------------

"""
    scani_open(;ip="191.30.80.131")

## Description
Open connection with Scanivalve DSA 3217
## Arguments
- ip: IP device
## Example
```jldoctest
julia> dev = open_scani();

# Example of configuration
julia> dev.config(AVG=1, PERIOD=1000, FPS=1000) # for 16.0 s acquisition
```
"""
function scani_open(;ip="191.30.80.131")
    scani = pyimport("scanivalve")
    s = scani.Scanivalve(ip=ip)
    return s
end

"""
    scani_acquire(s; sname=missing)

## Description
Function to acquire and save measured data
## Arguments
- s: scanivalve connection already configured
- sname: name of the file to be saved. If missing, then only returns the pressure, the frequency and the DateTime of the measurements
## Example
```jldoctest
julia> s = scani_open()

julia> s.config(AVG=1, PERIOD=1000, FPS=1000)

julia> scani_acquire(s, "teste.jld")
```
"""
function scani_acquire(s; sname=missing)
    p, f = s.acquire()
    h = now()

    if ismissing(sname)
        return p, f, h
    else
        save(sname, "pressure", p, "freq", f, "h", h)
        return p, f, h
    end
end

#---------------------------------------------------------------------------------------------------------------
# AUTOMATIC TEST WITH PREDEFINED positions
#---------------------------------------------------------------------------------------------------------------

"""
    testmatrix_onecylinder(;xlims::AbstractVector, ylim::AbstractVector, dx, dy)

## Description
Create a matrix with the measurements points for a test with one cylinder (no GA)
## Arguments
- xlims: vector with the low and the up boundaries for x direction
- ylims: vector with the low and the up boundaries for y direction
- dx, dy: x and y steps
## Example
```jldoctest
julia> testmatrix_onecylinder(xlims=[0,5],ylims=[-1,1], dx=5, dy=2)
4-element Array{OneCylType,1}:
 OneCylType(0.0, -1.0)
 OneCylType(0.0, 1.0)
 OneCylType(5.0, -1.0)
 OneCylType(5.0, 1.0)
```
"""
function testmatrix_onecylinder(;xlims::AbstractVector, ylims::AbstractVector, dx, dy)

    x = range(xlims[1], stop=xlims[2], step=dx)
    y = range(ylims[1], stop=ylims[2], step=dy)

    p = [OneCylType(i,j) for i in x for j in y]

    return p
end

"""
    autotest_ocylinder(s, robot, pos_vec::Vector{OneCylType}; sname="testdata.jld")

## Description
Perform multiple measurements in a given grid of positions (for 1 cylinder only).
## Arguments
- s: scanivalve connection configured
- robot: robot connection configured
- pos_vec: vector of OneCylType - grid of positions
- returntozero: if true then the cylinder returns back to the origin before new position
- sname: file name to be saved
## Example
```jldoctest
julia> s = scani_open();

julia> s.config(AVG=1, PERIOD=1000, FPS=1000);

julia> dev = stepper_init(1)

julia> grid = testmatrix_onecylinder(xlims=[0.0,90.0], ylims=[-30.0,20.0], dx=15.0, dy=5.0);

julia> d,f = autotest_1cylinder(s, dev, grid)
```
"""
function autotest_onecylinder(s, robot, pos_vec::Vector{OneCylType}; returntozero=false, sname="testdata.jld")

    coords = convert.(Array,pos_vec)
    pressure = Array{Array{Float64,2},1}()
    freq = Array{Float64,1}()
    testdate = Array{DateTime,1}()

    for pos in pos_vec
        println("Moving to position: x = ", pos.x, ", y = ", pos.y, "\n")
        if returntozero
            movecyl!(robot, OneCylType(0.0,0.0))
        end
        movecyl!(robot, pos)
        println("Acquiring data...\n")
        p, f, h = scani_acquire(s)
        push!(pressure, p)
        push!(freq, f)
        push!(testdate, h)
        println("Saving on file...\n")
        save(sname, "coords", coords, "pressure", pressure, "freq", freq, "h", testdate)
    end

    println("End of test")

    return pressure, freq, testdate
end

#---------------------------------------------------------------------------------------------------------------
# CROSS ENTROPY METHOD
#---------------------------------------------------------------------------------------------------------------

"""
    crossentropy(s, robot, pos_vec::AbstractVector, nbest::Int; mean_tol::AbstractVector, std_tol::AbstractVector, maxiter=10, dx=5.0, dy=0.0)

## Description
Cross entropy test for one cylinder
## Arguments
- sname: name of the file to be saved
- s: scanivalve connection configured
- robot: robot connection configured
- N: number of positions for the second iteration and on (for the first operation pos_vec is used)
- pos_vec: initial positions vector
- nbest: number of values which should be selected for the CE
- mean_tol: tolerance criteria for mean variation (for each direction)
- std_tol: tolerance criteria for standard deviation (for each direction)
- maxiter: maximum number of iterations
- dx, dy: horizontal and vertical gap between cylinders at the origin
## Output
Array with all the positions and mean step base pressure as well as the best position found
"""
function crossentropy(sname, s, robot, N::Int, pos_vec::AbstractVector, nbest::Int; mean_tol::AbstractVector,
                        std_tol::AbstractVector, maxiter=10)

    iter = 1
    df = []
    data = convert(DataFrame, pos_vec)
    data[!,:pm] = repeat([0.0], nrow(data))

    best = DataFrame()

    while iter <= maxiter
        println("Iteration ", iter)
        pm = []

        for pos in pos_vec
            println("Moving to position ", convert(Array, pos), "\n")
            strategy = beststrategy(robot, pos)
            movecyl!(robot, pos, strategy=strategy)
            println("Acquiring data...\n")
            p,_,_ = scani_acquire(s)
            push!(pm, mean(p[:,1:9]))
        end

        data.pm = pm
        push!(df, deepcopy(data))
        append!(best, data)
        best = best[best.pm .== maximum(best.pm),:]

        println(data)

        lims, converge = crossentropy(df, nbest; mean_tol=mean_tol, std_tol=std_tol)
        println(lims)
        println("converged? ", converge)

        if converge
            return df, best
            break
        else
            pos_vec = randpos(N, lims)
            data = convert(DataFrame, pos_vec)
            data[!,:pm] = repeat([0.0], nrow(data))
        end

        save(sname, "df", df, "best", best)
        iter += 1
    end

    return df, best
end

#---------------------------------------------------------------------------------------------------------------
# GENETIC ALGORITHM
#---------------------------------------------------------------------------------------------------------------
#
# DECIMAL-CODED
#
"""
    gadecimal(s, robot, pop_init::AbstractVector; Tm::AbstractVector, α0::AbstractVector, maxgen::Int, ucb::Int, lcb::Int, digits=0, boundaries::AbstractVector, mean_tol=AbstractVector, std_tol=AbstractVector, dx=5.0, dy=0.0, measureall=false)

## Description
Decimal-coded, mutation based genetic ALGORITHM
## Arguments
- s: scanivalve connection configured
- robot: robot connection configured
- pos_vec: initial position vector (CylType)
- Tm: mutation rate for each direction
- α0: shift for the first generation
- maxgen: total number of generations
- ucb: upper class boundary (elitism)
- lcb: lower class boundary
- digits: decimal places
- boundaries: domain limits for the direction
- mean_tol: vector with the tolerance for the mean values
- std_tol: vector with the tolerance for the std values
- dx, dy: horizontal and vertical gap between cylinders at the origin
- measureall: if true then all the positions are measured each iteration (even repeated ones)
"""
function gadecimal(sname, s, robot, pos_vec::AbstractVector; Tm::AbstractVector, α0::AbstractVector, maxgen::Int, ucb::Int,
        lcb::Int, digits=0, boundaries::AbstractVector, mean_tol::AbstractVector, std_tol::AbstractVector, measureall=false)

    gen = 1
    N = length(pos_vec)
    df = []
    tomeasure = 1:N

    df = []
    data = convert(DataFrame, pos_vec)
    data[!,:fitness] = repeat([0.0], N)

    while gen <= maxgen
        println("Generation ", gen)
        pm = data.fitness

        for k in tomeasure
            println("Moving to position ", convert(Array, pos_vec[k]), "\n")
            strategy = beststrategy(robot, pos_vec[k])
            movecyl!(robot, pos_vec[k], strategy=strategy)
            println("Acquiring data...\n")
            p,_,_ = scani_acquire(s)
            pm[k] = mean(p[:,1:9])
        end

        data.fitness = pm
        push!(df, deepcopy(data))

        println(data)

        converge = gadecimal_criteria(df, ucb=ucb, mean_tol=mean_tol, std_tol=std_tol)
        if converge
            return df
            break
        end

        data,tomeasure = gadecimal_step(data, Tm=Tm, α0=α0, gen=gen, maxgen=maxgen, ucb=ucb, lcb=lcb,
                        digits=digits, boundaries=boundaries)

        if measureall
            tomeasure = 1:N
        end

        pos_vec = convert(eltype(pos_vec), data[!,Not(:fitness)])

        save(sname, "df", df)
        gen += 1
    end

    return df
end


text:
const libGR_syms = [
:gr_setcallback
:gr_initgr
:gr_opengks
:gr_closegks
:gr_inqdspsize
:gr_openws
:gr_closews
:gr_activatews
:gr_deactivatews
:gr_clearws
:gr_updatews
:gr_polyline
:gr_polymarker
:gr_text
:gr_inqtext
:gr_fillarea
:gr_cellarray
:gr_nonuniformcellarray
:gr_polarcellarray
:gr_nonuniformpolarcellarray
:gr_gdp
:gr_path
:gr_spline
:gr_gridit
:gr_setlinetype
:gr_setlinewidth
:gr_setlinecolorind
:gr_setmarkertype
:gr_setmarkersize
:gr_setmarkercolorind
:gr_settextfontprec
:gr_setcharexpan
:gr_setcharspace
:gr_settextcolorind
:gr_setcharheight
:gr_inqcharheight
:gr_setcharup
:gr_settextpath
:gr_settextalign
:gr_setfillintstyle
:gr_setfillstyle
:gr_setfillcolorind
:gr_setcolorrep
:gr_setscale
:gr_inqscale
:gr_setwindow
:gr_inqwindow
:gr_setviewport
:gr_inqviewport
:gr_selntran
:gr_setclip
:gr_setwswindow
:gr_setwsviewport
:gr_createseg
:gr_copysegws
:gr_redrawsegws
:gr_setsegtran
:gr_closeseg
:gr_emergencyclosegks
:gr_updategks
:gr_setspace
:gr_textext
:gr_inqtextext
:gr_axes
:gr_axeslbl
:gr_grid
:gr_grid3d
:gr_verrorbars
:gr_herrorbars
:gr_polyline3d
:gr_polymarker3d
:gr_axes3d
:gr_titles3d
:gr_surface
:gr_contour
:gr_contourf
:gr_hexbin
:gr_setcolormap
:gr_colorbar
:gr_inqcolor
:gr_inqcolorfromrgb
:gr_hsvtorgb
:gr_tick
:gr_validaterange
:gr_adjustlimits
:gr_adjustrange
:gr_beginprint
:gr_beginprintext
:gr_endprint
:gr_ndctowc
:gr_wctondc
:gr_wc3towc
:gr_drawrect
:gr_fillrect
:gr_drawarc
:gr_fillarc
:gr_drawpath
:gr_setarrowstyle
:gr_setarrowsize
:gr_drawarrow
:gr_readimage
:gr_drawimage
:gr_importgraphics
:gr_setshadow
:gr_settransparency
:gr_setcoordxform
:gr_begingraphics
:gr_endgraphics
:gr_getgraphics
:gr_drawgraphics
:gr_mathtex
:gr_inqmathtex
:gr_setregenflags
:gr_inqregenflags
:gr_savestate
:gr_restorestate
:gr_selectcontext
:gr_destroycontext
:gr_uselinespec
:gr_delaunay
:gr_interp2
:gr_trisurface
:gr_tricontour
:gr_gradient
:gr_quiver
:gr_reducepoints
:gr_version
:gr_shadepoints
:gr_shadelines
:gr_setcolormapfromrgb
:gr_panzoom
:gr_setborderwidth
:gr_setbordercolorind
:gr_setprojectiontype
:gr_setperspectiveprojection
:gr_setorthographicprojection
:gr_settransformationparameters
:gr_setresamplemethod
:gr_setwindow3d
:gr_setspace3d
:gr_text3d
:gr_inqtext3d
:gr_settextencoding
:gr_inqtextencoding
:gr_loadfont
:gr_inqvpsize
:gr_inqtransformationparameters
:gr_polygonmesh3d
]


text:
abstract type RotationType{T} end

"""
    Rotation(Quaternion(1., 0., 0., 0.))

Rotation around the origin implemented from a rotation type `R`.
`R` can be any kind of rotation which implements a method for
`rotate(x, r::R) where {R<:RotationType}`. In 3D, this can be
`Quaternion` or `Euler`, for example.
"""
struct Rotation{T,R<:RotationType{T}} <: Transform
  rot::R
end
(r::Rotation)(p::AbstractVector) = r.rot * p

(≈)(x::Rotation, y::Rotation) = x.rot ≈ y.rot

inv(r::Rotation) = typeof(r)(-r.rot)

const Rotated{O,R<:Rotation} = Transformed{O,R}
Rotated(obj::O, transf::R) where {O,R} = Rotated{O,R}(obj, transf)

struct Quaternion{T} <: RotationType{T}
  coords::SVector{4,T}
end

Quaternion(coords::AbstractArray) = Quaternion(SVector{length(coords),eltype(coords)}(coords))
Quaternion(coords::Vararg{T}) where {T} = Quaternion(SVector{length(coords),T}(coords))


text:
function MOI.supports_constraint(
    ::Optimizer,
    ::Type{MOI.VectorAffineFunction{Float64}},
    ::Type{<:MOI.IndicatorSet{A, S}}
) where {
    A, S <: Union{MOI.LessThan{Float64}, MOI.GreaterThan{Float64}, MOI.EqualTo{Float64}}
}
    return true
end

function MOI.is_valid(
    model::Optimizer,
    c::MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, S}
) where {
    S <: MOI.IndicatorSet
}
    info_tuple = get(model.indicator_constraint_info, c.value, nothing)
    if info_tuple === nothing
        return false
    end
    return typeof(info_tuple[1].set) == S
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintSet,
    c::MOI.ConstraintIndex{<:MOI.VectorAffineFunction, <:MOI.IndicatorSet}
)
    MOI.throw_if_not_valid(model, c)
    info = model.indicator_constraint_info[c.value][1]
    return info.set
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintFunction,
    c::MOI.ConstraintIndex{<:MOI.VectorAffineFunction, <:MOI.IndicatorSet}
)
    MOI.throw_if_not_valid(model, c)
    func = model.indicator_constraint_info[c.value][2]
    return func
end

function MOI.add_constraint(
    model::Optimizer,
    func::MOI.VectorAffineFunction{Float64},
    s::MOI.IndicatorSet{A, S}
) where {
    A, S <: Union{MOI.GreaterThan{Float64}, MOI.LessThan{Float64}, MOI.EqualTo{Float64}}
}
    first_index_terms  = [v.scalar_term for v in func.terms if v.output_index == 1]
    scalar_index_terms = [v.scalar_term for v in func.terms if v.output_index != 1]
    if length(first_index_terms) != 1
        throw(ArgumentError("There should be exactly one term in output_index 1, found $(length(first_index_terms))"))
    end
    # getting the indicator variable
    sense_code = _get_indicator_sense(s.set)
    active_bool = A == MOI.ACTIVATE_ON_ONE ? Cint(0) : Cint(1)
    var_inner = _info(model, first_index_terms[1].variable_index).column

    # keeping track of the constraints
    model.last_constraint_index += 1
    info = ConstraintInfo(length(model.indicator_constraint_info) + 1, s)
    model.indicator_constraint_info[model.last_constraint_index] = (info, func)

    linear_idx = Vector{Cint}(undef, length(scalar_index_terms))
    linear_coefficients = Vector{Cdouble}(undef, length(scalar_index_terms))
    for (i, term) in enumerate(scalar_index_terms)
        linear_idx[i] = _info(model, term.variable_index).column
        linear_coefficients[i] = term.coefficient
    end
    # switching constant on other side
    # a^T x + b <= c <==> a^T x <= c - b
    fcons = MOI.constant(func)
    rhs = Cdouble(MOI.constant(s.set) - fcons[2])
    add_indicator_constraint(model.inner, linear_idx, linear_coefficients, sense_code, rhs, var_inner, active_bool)
    return MOI.ConstraintIndex{typeof(func), typeof(s)}(model.last_constraint_index)
end

_get_indicator_sense(::MOI.GreaterThan) = Cint('G')
_get_indicator_sense(::MOI.LessThan) = Cint('L')
_get_indicator_sense(::MOI.EqualTo) = Cint('E')

function MOI.get(
    model::Optimizer,
    ::MOI.ListOfConstraintIndices{<:MOI.VectorAffineFunction, S},
) where {
    S <: MOI.IndicatorSet
}
    indices = MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, S}[
        MOI.ConstraintIndex{MOI.VectorAffineFunction{Float64}, S}(key)
        for (key, info_func) in model.indicator_constraint_info
        if typeof(info_func[1].set) === S
    ]
    return sort!(indices, by = x -> x.value)
end

function MOI.get(
    model::Optimizer,
    ::MOI.NumberOfConstraints{MOI.VectorAffineFunction{Float64}, S}
) where {
    S <: MOI.IndicatorSet
}
    nindices = 0
    for (key, info_func) in model.indicator_constraint_info
        info = info_func[1]
        if typeof(info.set) === S
            nindices += 1
        end
    end
    return nindices
end

function MOI.get(
    model::Optimizer,
    ::MOI.ConstraintName,
    c::MOI.ConstraintIndex{<:MOI.VectorAffineFunction, <:MOI.IndicatorSet}
)
    MOI.throw_if_not_valid(model, c)
    info = model.indicator_constraint_info[c.value][1]
    return info.name
end

function MOI.set(
    model::Optimizer,
    ::MOI.ConstraintName,
    c::MOI.ConstraintIndex{<:MOI.VectorAffineFunction, S},
    name::String,
) where {
    S <: MOI.IndicatorSet
}
    MOI.throw_if_not_valid(model, c)
    info = model.indicator_constraint_info[c.value][1]
    info.name = name
    return
end


text:
mutable struct Point
    x::Float64 # Coordonnée x
    y::Float64 # Coordonnée x
    id::String # Identifiant
end

mutable struct Solution
    # Structure définissant une solution (TSP/TSP-D)
    ordreVisiteCamion::Vector{Int64}
    operationsDrone::Vector{Vector{Int64}}
    tempsParcours::Float64 # Temps de parcours total
    tempsAttente::Float64 # Temps d'attente total (drone + camion)
end

mutable struct Instance
    vitesseCamion::Float64
    vitesseDrone::Float64
    nbPoints::Int64
    listePoints::Vector{Point}
    D::Array{Float64,2} # Distancier
    solution::Solution
    codeResolution::Int64 # Permet de savoir quel est l'état courant de la résolution
    # 0 = Initialisation (càd. pas de résolution)
    # 1 = TSP : solution exacte
    # 2 = TSP : solution approchée LK
    # 3 = TSP-D : AEP mono-objectif
    # 4 = TSP-D : EP mono-objectif
    # 5 = TSP-D : AEP bi-objectif
end

max_stars_repo_name: xgandibleux/tspd-solver

max_issues_repo_name: xgandibleux/tspd-solver

max_forks_repo_name: xgandibleux/tspd-solver


text:
export chemical_formula

"""
Returns the chemical formula of an AbstractSystem as a string.
"""
function chemical_formula(symbols::AbstractVector{Symbol})
    parts = map(collect(Set(symbols))) do sym
        sym_count = count(isequal(sym), symbols)
        sym_count < 2 && return string(sym)

        str_count = string(sym_count)
        for i in 0:9
            str_count = replace(str_count, ('0' + i) => ('₀' + i))  # Make subscripts
        end
        string(sym) * str_count
    end
    join(sort(parts))
end
chemical_formula(system) = chemical_formula(atomic_symbol(system))

function show_system(io::IO, system::AbstractSystem{D}) where {D}
    print(io, "($(chemical_formula(system)), ")
    bc = boundary_conditions(system)
    if all(isequal(bc[1]), bc)
        print(io, typeof(bc[1]), ", ")
    end
    box = bounding_box(system)
    if box != infinite_box(D)
        box_str = ["[" * join(ustrip.(bvector), ", ") * "]" for bvector in box]
        print(io, "box=[", join(box_str, ", "), "]u\"$(unit(box[1][1]))\"")
    else
        print(io, "box=infinite")
    end
    print(io, ")")
end


text:
using VisualParticipationAnalytics
using Clustering
using DataFrames
using Dates
using Distances
using Embeddings
using Languages
using LinearAlgebra
using MLJ
using MultivariateStats
using NearestNeighbors
using PGFPlotsX
using Random
using SQLite
using Statistics
using StatsBase
using TextAnalysis

Random.seed!(1)

CONFIG = Dict(:out_dir=>"out", :do_lsa=>false, :do_lda=>false, :do_cluster=>true)

"""
"""
function initplots()
    empty!(PGFPlotsX.CUSTOM_PREAMBLE)
    push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"""
        \usepgfplotslibrary{colorbrewer}
        \pgfplotsset{
            colormap/Paired-12,
            cycle list/Paired-12,
            cycle multiindex* list={
                mark list*\nextlist
                Paired-12\nextlist
            }
        }
        \usepackage{libertine}
        \usepackage{unicode-math}
        \setmathfont[Scale=MatchUppercase]{libertinusmath-regular.otf}
        \newcommand{\mlv}[1]{\mathit{#1}}
    """)
end

"""
"""
function creategraph(dbpath, tablename)
    # Load data
    db = SQLite.DB(dbpath)
    df = DBInterface.execute(db, """
        SELECT * FROM $tablename;
    """) |> DataFrame

    for i = 1:5
        println(df[i])
    end
end

"""
"""
function preparedata(dbpath, tablename)
    # Load data
    df = readdb(dbpath, tablename)
    display(describe(df))

    # Pre-process data
    # Textual data
    lang = Languages.German()
    authors = df[:, ["author"]]
    #categorical!(authors)
    #ohe = OneHotEncoder()
    #mach = MLJ.fit!(machine(ohe, authors))
    #authors_ohe = MLJ.transform(mach, authors)
    #df = [df authors_ohe]

    titles = df[:, "title"]
    titles = preprocess.(titles, [lang])
    tfidf_titles = tf_idf(titles)
    crps_titles = Corpus(StringDocument.(titles))
    update_lexicon!(crps_titles)
    dtm_titles = DocumentTermMatrix(crps_titles)

    contents = df[:, "content"]
    contents = preprocess.(contents, [lang])
    tfidf_contents = tf_idf(contents)
    crps_contents = Corpus(StringDocument.(contents))
    update_lexicon!(crps_contents)
    dtm_contents = DocumentTermMatrix(crps_contents)

    # Temporal data
    # created_at

    # Other data
    #
    #category = df[:, ["category"]]
    #categorical!(category)
    #ohe = OneHotEncoder()
    #mach = MLJ.fit!(machine(ohe, category))
    #category_ohe = MLJ.transform(mach, category)
    #df = [df category_ohe]

    return df, crps_titles, crps_contents, dtm_titles, dtm_contents, tfidf_titles, tfidf_contents
end

"""
"""
function assignmentplot(assignments, x, y, xlabel, ylabel)
    t = @pgf PGFPlotsX.Table({ meta="cluster" }, x=x, y=y, cluster=assignments)
    p = @pgf Plot({
        scatter,
        "only marks",
        scatter_src="explicit",
        mark_size="1pt"
    }, t)
    return @pgf Axis({ xlabel=xlabel, ylabel=ylabel }, p)
end

"""
"""
function clusterkmeans(data, ks, dist, x, y, name)
    clusterings = []
    Js = []
    silhouette_coefficients = []

    distances = pairwise(dist, data')
    for k in ks
        clustering = kmeans(data', k, distance=dist)
        push!(clusterings, clustering)
        push!(Js, clustering.totalcost)
        push!(silhouette_coefficients, mean(silhouettes(clustering, distances)))
    end

    ## Assignment plots
    for k in ks
        plt = assignmentplot(assignments(clusterings[k - 1]), x, y, "longitude",
                "latitude")
        pgfsave(CONFIG[:out_dir] * "/kmeans_$(name)_assignments_k_" * (k < 10 ? "0$k" : "$k") * ".pdf", plt)
    end

    ## Elbow method plot
    c = @pgf Coordinates(zip(ks, Js))
    p = @pgf Plot({ color="red", mark="x" }, c)
    plt = @pgf Axis({ xlabel=raw"\(k\)", ylabel=raw"\(J\)" }, p)
    pgfsave(CONFIG[:out_dir] * "/kmeans_$(name)_elbow.pdf", plt)

    ## Silhouette coefficient plot
    c = @pgf Coordinates(zip(ks, silhouette_coefficients))
    p = @pgf Plot({ color="red", mark="x" }, c)
    plt = @pgf Axis({ xlabel=raw"\(k\)", ylabel="silhouette coefficient" }, p)
    pgfsave(CONFIG[:out_dir] * "/kmeans_$(name)_silhouette.pdf", plt)
end

"""
"""
function clusterdbscan(data, minpoints, epsilons, dist, optimize, x, y, name)
    eps = 0.5825
    #silhouette_coefficients = []
    for minpts in minpoints
        if optimize
            tree = BallTree(data', dist)
            _, knn_distances = knn(tree, data', minpts + 1, true)
            mean_knn_distances = mean.(knn_distances)
            sort!(mean_knn_distances)

            n = length(mean_knn_distances)
            ds = []
            first_to_last = [n, mean_knn_distances[n]] - [1, mean_knn_distances[1]]
            normalize!(first_to_last)
            for i = 1:n
                first_to_i = [i, mean_knn_distances[i]] - [1, mean_knn_distances[1]]
                d = norm(first_to_i - dot(first_to_i, first_to_last) * first_to_last)
                push!(ds, d)
            end
            eps = mean_knn_distances[argmax(ds)]

            c = @pgf Coordinates(zip(1:n, mean_knn_distances))
            p = @pgf Plot({ color="red" }, c)
            plt = @pgf Axis({ xlabel="instance", ylabel="$minpts-nn distance" }, p,
                    HLine({ dotted }, eps))
            pgfsave(CONFIG[:out_dir] * "/dbscan_$(name)_eps_min_pts_" * (minpts < 10 ? "0" : "") * "$minpts.pdf", plt)
        end

        distances = pairwise(dist, data')
        for ϵ in epsilons
            clustering = dbscan(distances, ϵ, minpts)
            assigns = assignments(clustering)
            #filter!(!iszero, assigns)
            assigns = assigns .+ 1
            #push!(silhouette_coefficients, mean(silhouettes(assigns, distances)))

            save(clustering, CONFIG[:out_dir] * "/" * (minpts < 10 ? "0" : "") * "$(minpts)_$(ϵ).json")

            ## Assignment plots
            plt = assignmentplot(assignments(clustering), x, y, "longitude", "latitude")
            pgfsave(CONFIG[:out_dir] * "/dbscan_$(name)_assignments_min_pts_" * (minpts < 10 ? "0" : "")
                    * "$(minpts)" * "_eps_" * string(floor(ϵ, digits=3)) * ".pdf" , plt)
        end
    end
end

"""
"""
function process(dbpath, tablename)
    CONFIG[:out_dir] = CONFIG[:out_dir] * "/" * split(dbpath, "/")[end]
    !ispath(CONFIG[:out_dir]) && mkpath(CONFIG[:out_dir])

    df, crps_titles, crps_contents, dtm_titles, dtm_contents, tfidf_titles, tfidf_contents = preparedata(dbpath, tablename)
    display(dtm_titles)
    display(dtm_contents)
    display(Array(tfidf_titles))
    display(Array(tfidf_contents))

    #M = fit(PCA, Array(tfidf_contents)'; maxoutdim=1000)
    #tfidf_contents = MultivariateStats.transform(M, tfidf_contents')'

    if CONFIG[:do_lsa]
        lsa(tfidf_titles)
        lsa(tfidf_contents)
        M = fit(UnitRangeTransform, tfidf_contents, dims=2)
        tfidf_contents = StatsBase.transform(M, tfidf_contents)
    end
    if CONFIG[:do_lda]
        lda_titles = lda(dtm_titles, 4, 1000, 0.1, 0.1)
        lda_contents = lda(dtm_contents, 4, 1000, 0.1, 0.1)
        for i = 1:4
            println("TOPIC $i")
            display(topkwords(lda_titles[1], i, crps_titles, 10))
            display(topkwords(lda_contents[1], i, crps_contents, 10))
        end
    end

    # Cluster data
    if CONFIG[:do_cluster]
        longitude = df[:, "long"]
        latitude = df[:, "lat"]
        clusterdbscan([longitude latitude], 2:12, 5:5:200, Haversine(), false, longitude, latitude, "long_lat")
        #clusterkmeans(tfidf_contents, 2:15, Euclidean(), longitude, latitude, "tfidf_contents_euclidean")
        #clusterkmeans(tfidf_contents, 2:15, CosineDist(), longitude, latitude, "tfidf_contents_cosine")
        #clusterkmedoids(tfidf_contents, 2:12, CosineDist(), longitude, latitude, "tfidf_contents_cosine")
    end
end

"""
"""
function main()
    initplots()
    #process("~/Datasets/participation/databases/liqd_laermorte_melden.sqlite", "contribution")
    process("~/Datasets/participation/databases/liqd_mauerpark.sqlite", "contribution")
    #process("~/Datasets/participation/databases/liqd_blankenburger_sueden.sqlite", "comment_a")
    #process("~/Datasets/participation/databases/liqd_blankenburger_sueden.sqlite", "comment_b")
    #process("~/Datasets/participation/databases/liqd_blankenburger_sueden.sqlite", "comment_c")
end

"""
"""
function distances()
    !ispath(CONFIG[:out_dir]) && mkpath(CONFIG[:out_dir])

    df = readdb("~/Datasets/participation/databases/liqd_mauerpark.sqlite", "contribution")
    longitude = df[:, "long"]
    latitude = df[:, "lat"]
    distances = pairwise(Haversine(), [longitude latitude]')
    save(distances, CONFIG[:out_dir] * "/long_lat_haversine.json")
    sorteddistances = createrankings(distances)
    save(Dict(sorteddistances), CONFIG[:out_dir] * "/long_lat_haversine_sorted.json")

    createdat = df[:, "created_at"]
    createdat = map(t -> split(t, ".")[1], createdat)
    format = DateFormat("Y-m-dTH:M:S", "english")
    timestamps = DateTime.(createdat, format)
    n = length(timestamps)
    distances = zeros(Millisecond, n, n)
    for i = 1:n
        for j = 1:n
            distances[i, j] = timestamps[i] - timestamps[j]
        end
    end
    save(string.(canonicalize.(Dates.CompoundPeriod.(distances))), CONFIG[:out_dir] * "/timestamps_en.json")
    save(map(t -> t.value, round.(distances, Second)), CONFIG[:out_dir] * "/timestamps_seconds.json")
    save(map(t -> t.value, round.(distances, Minute)), CONFIG[:out_dir] * "/timestamps_minutes.json")
    save(map(t -> t.value, round.(distances, Hour)), CONFIG[:out_dir] * "/timestamps_hours.json")
    save(map(t -> t.value, round.(distances, Day)), CONFIG[:out_dir] * "/timestamps_days.json")

    authors = df[:, "author"]
    n = length(authors)
    distances = zeros(Int, n, n)
    for i = 1:n
        for j = 1:n
            distances[i, j] = authors[i] == authors[j] ? 0 : 1
        end
    end
    save(distances, CONFIG[:out_dir] * "/authors.json")
    sorteddistances = createrankings(distances)
    save(Dict(sorteddistances), CONFIG[:out_dir] * "/authors_sorted.json")

    titles = df[:, "title"]
    lang = Languages.German()
    titles = preprocess.(titles, [lang], dostem=false)
    embeddingtable = load_embeddings(FastText_Text{:de})
    wordindextable = Dict(word => index for (index, word) in enumerate(embeddingtable.vocab))
    embeddingsize = size(embeddingtable.embeddings, 1)
    titleembeddings = Matrix{Float64}(undef, embeddingsize, length(titles))
    for (index, title) in enumerate(titles)
        words = split(title, " ")
        titleembedding = zeros(embeddingsize)
        numunknownwords = 0
        for word in words
            if !haskey(wordindextable, word)
                numunknownwords += 1
                continue
            end
            wordindex = wordindextable[word]
            embedding = embeddingtable.embeddings[:, wordindex]
            titleembedding += embedding
        end
        numknownwords = length(words) - numunknownwords
        if numknownwords > 0
            titleembedding /= numknownwords
        end
        titleembeddings[:, index] = titleembedding
    end
    distances = pairwise(Euclidean(), titleembeddings)
    save(distances, CONFIG[:out_dir] * "/titles_euclidean.json")
    sorteddistances = createrankings(distances)
    save(Dict(sorteddistances), CONFIG[:out_dir] * "/titles_euclidean_sorted.json")
end


text:
using Colors

"""
    Dataset

A supertype for all loadable datasets
"""
abstract type Dataset end

"""
    VideoGraphicDataset{C <: Color} <: Dataset

A supertype for all videographic datasets
"""
abstract type VideoGraphicDataset{C <: Color} <: Dataset end


text:
# MAPI_acoustic_BZ.jl

push!(LOAD_PATH,"../../src") # Temporary versions of modules in PWD
using JuliaPhonons

P=read_POSCAR(open("../POSCAR"),expansion=[4,4,4])
eigenvectors,eigenmodes=read_meshyaml(open("../mesh.yaml"),P)

output_conflated_xyz(P,0,(eigenvectors[1],eigenvectors[1]),
                    (1,1.0),
                    steps=128, repeats=1, q=[0.25,0.25,0.25]) # generates files anim_{count}.xyz

output_conflated_xyz(P,1,(eigenvectors[2],eigenvectors[2]),
                    (1,1.0),
                    steps=128, repeats=1, q=[0.25,0.25,0.25]) # generates files anim_{count}.xyz

output_conflated_xyz(P,2,(eigenvectors[3],eigenvectors[3]),
                    (1,1.0),
                    steps=128, repeats=1, q=[0.25,0.25,0.25]) # generates files anim_{count}.xyz

end
output_conflated_xyz(P,0,(eigenvectors[1], eigenvectors[2],eigenvectors[3]),
                    (1,1,1),
                    steps=64, repeats=1, q=[0.25,0.25,0.25]) # generates files anim_{count}.xyz

max_stars_repo_name: jarvist/Julia-Phonons

max_issues_repo_name: lucydot/Julia-Phonons

max_forks_repo_name: lucydot/Julia-Phonons


text:
function draw_wireframe_sphere(pos, r; npts=11, color=(255,255,255), name=nothing, opacity=1.)
    ϕs = range(0, stop=2π, length=npts)
    θs = range(-π/2, stop=π/2, length=npts)
    # vertical
    v_xs = [[r * cos(θ) * sin(ϕ) + pos[1] for θ in θs] for ϕ in ϕs[1:end-1]]
    v_ys = [[r * cos(θ) * cos(ϕ) + pos[2] for θ in θs] for ϕ in ϕs[1:end-1]]
    v_zs = [[r * sin(θ) + pos[3] for θ in θs] for ϕ in ϕs[1:end-1]]
    # horizontal
    h_xs = [[r * cos(θ) * sin(ϕ) + pos[1] for ϕ in ϕs] for θ in θs]
    h_ys = [[r * cos(θ) * cos(ϕ) + pos[2] for ϕ in ϕs] for θ in θs]
    h_zs = [[r * sin(θ) + pos[3] for ϕ in ϕs] for θ in θs]

    hover = isnothing(name) ? "skip" : nothing

    xs, ys, zs = Float64[], Float64[], Float64[]
    for i=1:length(v_xs)
        append!(xs, v_xs[i])
        push!(xs, NaN)
        append!(ys, v_ys[i])
        push!(ys, NaN)
        append!(zs, v_zs[i])
        push!(zs, NaN)
    end
    for j=1:length(h_xs)
        append!(xs, h_xs[j])
        push!(xs, NaN)
        append!(ys, h_ys[j])
        push!(ys, NaN)
        append!(zs, h_zs[j])
        push!(zs, NaN)
    end

    return scatter3d(;x=xs, y=ys, z=zs,
                      mode="lines",
                      line=attr(color="rgb$color"),
                      opacity=opacity,
                      showlegend=false, showscale=false, name=name,
                      hovertemplate=name, hoverinfo=hover,
            )
end

draw_central_body(bd::Union{CelestialBody,Star}; kwargs...
    ) = draw_wireframe_sphere(SVector{3,Float64}(0.,0.,0.), bd.eqradius; color=bd.color, name=bd.name, kwargs...)

draw_orbiting_body(bd::CelestialBody, t; kwargs...
    ) = draw_wireframe_sphere(time_orbital_position(t, bd.orbit), bd.eqradius; color=bd.color, name=bd.name, kwargs...)

draw_soi(bd::CelestialBody, t; kwargs...
    ) = draw_wireframe_sphere(time_orbital_position(t, bd.orbit), bd.SoI; color=bd.color, name=bd.name, opacity=0.25, kwargs...)


text:
using Test
# using Base.Test
# using Base.@__FILE__

include("../src/mlst_download_functions.jl")
include("../src/build_db_functions.jl")
include("../src/calling_functions.jl")
include("../src/db_graph.jl")

TEST_DIR = (dirname(@__FILE__))
TMPDIR = mktempdir()
VERSION = "testing"
# encapsulate in a big test to get a summary of all in the end; if a testset fails, it will show the specific testset results,
# otherwise will show just the total # of test passes.
@testset "mentalist_tests" begin

  c_jejuni_pubmlst_dir = string(TMPDIR, "/", "c_jejuni_pubmlst")
  c_jejuni_pubmlst_db_file = string(TMPDIR, "/", "c_jejuni_pubmlst_31.db")
  c_jejuni_pubmlst_copy_dir = string(TMPDIR, "/", "c_jejuni_pubmlst_copy")
  c_jejuni_fastq_file = string("../data/SRR5824107_small.fastq.gz")
  l_pneumophila_cgmlst_dir = string(TMPDIR, "/", "l_pneumophila_cgmlst")
  l_pneumophila_cgmlst_db_file = string(TMPDIR,"/", "l_pneumophila_cgmlst_31.db")
  l_pneumophila_fastq_file = string("../data/SRR801793.fastq.gz")

  c_jejuni_pubmlst_loci_files = ""
  @testset "download_pubmlst_scheme" begin
    mkdir(c_jejuni_pubmlst_dir)
    c_jejuni_pubmlst_loci_files, c_jejuni_profile_file = download_pubmlst_scheme("Campylobacter jejuni", c_jejuni_pubmlst_dir)
    c_jejuni_allele_filenames = open(string(TEST_DIR, "/", "test_data", "/", "c_jejuni_pubmlst_allele_filenames.txt"))
    for c_jejuni_allele_filename in eachline(c_jejuni_allele_filenames)
       @test isfile(string(c_jejuni_pubmlst_dir, "/", chomp(c_jejuni_allele_filename)))
    end
  end

  l_pneumophila_cgmlst_loci_files = ""
  @testset "download_cgmlst_scheme" begin
    mkdir(l_pneumophila_cgmlst_dir)
    l_pneumophila_cgmlst_loci_files = download_cgmlst_scheme("Legionella pneumophila", l_pneumophila_cgmlst_dir)
    l_pneumophila_allele_filenames = open(string(TEST_DIR, "/", "test_data", "/", "l_pneumophila_cgmlst_allele_filenames.txt"))
    for l_pneumophila_allele_filename in eachline(l_pneumophila_allele_filenames)
      @test isfile(string(l_pneumophila_cgmlst_dir, "/", chomp(l_pneumophila_allele_filename)))
    end
  end

  @testset "build_l_pneumophila_db" begin
    K::Int8 = 31
    results, loci = kmer_class_for_each_locus(DNAKmer{K}, l_pneumophila_cgmlst_loci_files, 1)
    @test typeof(results) == Vector{Tuple{Dict{UInt64,Vector{Int16}},Vector{Int16},Dict{UInt64,Int64},Vector{Int64}}}
    @test typeof(loci) == Vector{String}

    kmer_classification = combine_loci_classification(DNAKmer{K}, results, loci)
    loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus = kmer_classification
    @test typeof(loci_list) == Vector{Int32}
    @test typeof(weight_list) == Vector{Int16}
    @test typeof(alleles_list) == Vector{Int16}
    @test typeof(kmer_list) == Vector{String}
    @test typeof(allele_ids_per_locus) == Vector{Int64}
    @test length(weight_list) == length(kmer_list)
    @test all([x != 0 for x in weight_list])
    @test all([length(kmer) == K for kmer in kmer_list])

    profile = nothing
    args = Dict("k" => K, "fasta_files" => l_pneumophila_cgmlst_loci_files)
    save_db(DNAKmer{K}, kmer_classification, loci, l_pneumophila_cgmlst_db_file, profile, args, VERSION)
    @test isfile(l_pneumophila_cgmlst_db_file)
  end

  @testset "build_c_jejuni_db" begin
    K::Int8 = 31
    results, loci = kmer_class_for_each_locus(DNAKmer{K}, c_jejuni_pubmlst_loci_files, 1)
    @test typeof(results) == Vector{Tuple{Dict{UInt64,Vector{Int16}},Vector{Int16},Dict{UInt64,Int64},Vector{Int64}}}
    @test typeof(loci) == Vector{String}

    kmer_classification = combine_loci_classification(DNAKmer{K}, results, loci)
    loci_list, weight_list, alleles_list, kmer_list, allele_ids_per_locus = kmer_classification
    @test typeof(loci_list) == Vector{Int32}
    @test typeof(weight_list) == Vector{Int16}
    @test typeof(alleles_list) == Vector{Int16}
    @test typeof(kmer_list) == Vector{String}
    @test typeof(allele_ids_per_locus) == Vector{Int64}
    @test length(weight_list) == length(kmer_list)
    @test all([x != 0 for x in weight_list])
    @test all([length(kmer) == K for kmer in kmer_list])

    profile = nothing
    args = Dict("k" => K, "fasta_files" => c_jejuni_pubmlst_loci_files)
    save_db(DNAKmer{K}, kmer_classification, loci, c_jejuni_pubmlst_db_file, profile, args, VERSION)
    @test isfile(c_jejuni_pubmlst_db_file)
  end

  @testset "call_l_pneumophila" begin
    # open the db:
    kmer_db, loci, loci2alleles, coverage, k, profile, build_args = open_db(l_pneumophila_cgmlst_db_file)
    @test typeof(kmer_db) == Dict{Kmer{DNA,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [l_pneumophila_fastq_file])
    # test kmer counts:
    @test kmer_count[DNAKmer{k}("AACTAAAACACTGGCTTTACTGCCCCAAATG")] == 19
    @test kmer_count[DNAKmer{k}("TTATGCCTATGAAAATATGATTTTATCTGAA")] == 24
    @test kmer_count[DNAKmer{k}("AAAGGAATTACAATGAATATCTAAATAATTC")] == 14
    @test kmer_count[DNAKmer{k}("AGGAAGTTTTGGTTATATGCTATATGATGCG")] == 21
    @test kmer_count[DNAKmer{k}("AAAAGAGCCAATTTGATGAGCCACCACATCA")] == 20

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles, coverage)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DefaultDict{Int16,Int64,Int64}

    # some parameters:
    kmer_thr, max_mutations, output_votes = 1, 5, true
    # call:
    allele_calls, voting_result = call_alleles(DNAKmer{k}, kmer_count, votes, loci_votes, loci, loci2alleles, l_pneumophila_cgmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test allele_calls[1].allele == "1"
    @test allele_calls[1].depth == 6
    @test allele_calls[2].allele == "1"
    @test allele_calls[2].depth == 9
    @test allele_calls[132].allele == "N" # novel_allele
    @test allele_calls[662].allele == "N" # novel_allele
    @test allele_calls[536].allele == "1"
    @test allele_calls[536].flag == "+" # multiple alleles possible
    @test allele_calls[652].allele == "1"
    @test allele_calls[652].flag == "+" # multiple alleles possible
  end

  @testset "call_c_jejuni" begin
    # open the db:
    kmer_db, loci, loci2alleles, coverage, k, profile, build_args = open_db(c_jejuni_pubmlst_db_file)
    @test typeof(kmer_db) == Dict{Kmer{DNA,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [c_jejuni_fastq_file])

    # test kmer counts:
    @test kmer_count[DNAKmer{k}("TCATTTAAGGACTTTTCAGTGATTAAAATCA")] == 2
    @test kmer_count[DNAKmer{k}("CACTCCAATTTTTTCAAATAAAGTAGCTAAG")] == 0
    @test kmer_count[DNAKmer{k}("ATTCTTTTACTCCTATTATCGGTTATACTAA")] == 0
    @test kmer_count[DNAKmer{k}("GAAAAAAGTAATCCAAGGTGCGCAAAAAGCA")] == 4
    @test kmer_count[DNAKmer{k}("AAATATAGTCAATAAATTATAAAAAAAACTT")] == 0

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles, coverage)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DefaultDict{Int16,Int64,Int64}

    # some parameters:
    kmer_thr, max_mutations, output_votes = 2, 5, true
    # call:
    allele_calls, voting_result = call_alleles(DNAKmer{k}, kmer_count, votes, loci_votes, loci, loci2alleles, c_jejuni_pubmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test [ac.allele for ac in allele_calls] == ["2", "17", "2", "3", "2", "1", "5"]
  end

  @testset "call_c_jejuni_after_moving_db" begin
    mkdir(c_jejuni_pubmlst_copy_dir)
    cp(c_jejuni_pubmlst_db_file, string(c_jejuni_pubmlst_copy_dir, "/", basename(c_jejuni_pubmlst_db_file)))
    cp(c_jejuni_pubmlst_dir, string(c_jejuni_pubmlst_copy_dir, "/", basename(c_jejuni_pubmlst_dir)))
    # open the db:
    kmer_db, loci, loci2alleles, coverage, k, profile, build_args = open_db(string(dirname(c_jejuni_pubmlst_copy_dir), "/", basename(c_jejuni_pubmlst_db_file)))
    @test typeof(kmer_db) == Dict{Kmer{DNA,31},Vector{Tuple{Int16,Int16,Vector{Int16}}}}
    @test typeof(loci) == Vector{String}
    @test typeof(loci2alleles) == Dict{Int16,Vector{Int16}}

    # count kmers:
    kmer_count = count_kmers(DNAKmer{k}, [c_jejuni_fastq_file])

    # test kmer counts:
    @test kmer_count[DNAKmer{k}("TCATTTAAGGACTTTTCAGTGATTAAAATCA")] == 2
    @test kmer_count[DNAKmer{k}("CACTCCAATTTTTTCAAATAAAGTAGCTAAG")] == 0
    @test kmer_count[DNAKmer{k}("ATTCTTTTACTCCTATTATCGGTTATACTAA")] == 0
    @test kmer_count[DNAKmer{k}("GAAAAAAGTAATCCAAGGTGCGCAAAAAGCA")] == 4
    @test kmer_count[DNAKmer{k}("AAATATAGTCAATAAATTATAAAAAAAACTT")] == 0

    # call:
    votes, loci_votes = count_votes(kmer_count, kmer_db, loci2alleles, coverage)
    @test typeof(votes) == Dict{Int16,Dict{Int16,Int64}}
    @test typeof(loci_votes) == DefaultDict{Int16,Int64,Int64}

    # some parameters:
    kmer_thr, max_mutations, output_votes = 2, 5, true
    # call:
    allele_calls, voting_result = call_alleles(DNAKmer{k}, kmer_count, votes, loci_votes, loci, loci2alleles, c_jejuni_pubmlst_loci_files, kmer_thr, max_mutations, output_votes)
    @test typeof(allele_calls) == Vector{AlleleCall}
    @test [ac.allele for ac in allele_calls] == ["2", "17", "2", "3", "2", "1", "5"]
  end

  @testset "complement_alleles" begin
    @test complement_alleles([], 0) == Int16[]
    @test complement_alleles([], 1) == [Int16(1)]
    @test complement_alleles([-1], 0) == Int16[]
    @test complement_alleles([-1], 1) == [Int16(1)]
    @test complement_alleles([0], 0) == Int16[]
    @test complement_alleles([0], 1) == [Int16(1)]
    @test complement_alleles([1], 0) == Int16[]
    @test complement_alleles([1], 1) == Int16[]
    @test complement_alleles([0,1,3,5,7,9], 10) == Int16[2,4,6,8,10]
    @test complement_alleles([0,2,4,6,8,10], 10) == Int16[1,3,5,7,9]
  end

  @testset "twin" begin
    @test twin(DNAKmer("A")) == DNAKmer("T")
    @test twin(DNAKmer("C")) == DNAKmer("G")
    @test twin(DNAKmer("G")) == DNAKmer("C")
    @test twin(DNAKmer("T")) == DNAKmer("A")
    @test twin(DNAKmer("ACTG")) == DNAKmer("CAGT")
  end

  # rm(TMPDIR, recursive=true)
end

max_stars_repo_name: liberuser/MentaLiST

max_issues_repo_name: liberuser/MentaLiST

max_forks_repo_name: liberuser/MentaLiST


text:
# Temperature and pressure coupling

export
    NoCoupling,
    apply_coupling!,
    AndersenThermostat,
    RescaleThermostat,
    BerendsenThermostat

"""
    NoCoupling()

Placeholder coupler that does nothing.
"""
struct NoCoupling end

"""
    apply_coupling!(system, simulator, coupling)

Apply a coupler to modify a simulation.
Custom couplers should implement this function.
"""
apply_coupling!(sys, sim, ::NoCoupling) = sys

"""
    AndersenThermostat(temperature, coupling_const)

Rescale random velocities according to the Andersen thermostat.
"""
struct AndersenThermostat{T, C}
    temperature::T
    coupling_const::C
end

function apply_coupling!(sys, sim, thermostat::AndersenThermostat)
    for i in 1:length(sys)
        if rand() < (sim.dt / thermostat.coupling_const)
            mass = sys.atoms[i].mass
            sys.velocities[i] = velocity(mass, thermostat.temperature;
                                            dims=n_dimensions(sys))
        end
    end
    return sys
end

"""
    RescaleThermostat(temperature)

The velocity rescaling thermostat.
This thermostat should not be used in general as it can lead to simulation artifacts.
"""
struct RescaleThermostat{T}
    temperature::T
end

function apply_coupling!(sys, sim, thermostat::RescaleThermostat)
    sys.velocities *= sqrt(thermostat.temperature / temperature(sys))
    return sys
end

"""
    BerendsenThermostat(temperature, coupling_const)

The Berendsen thermostat.
This thermostat should not be used in general as it can lead to simulation artifacts.
"""
struct BerendsenThermostat{T, C}
    temperature::T
    coupling_const::C
end

function apply_coupling!(sys, sim, thermostat::BerendsenThermostat)
    λ2 = 1 + (sim.dt / thermostat.coupling_const) * ((thermostat.temperature / temperature(sys)) - 1)
    sys.velocities *= sqrt(λ2)
    return sys
end


text:
module EMLflow_measurement
	using DanaTypes
	using DotPlusInheritance
	using Reexport
	@reexport using ...streams.EMLstreams
	import EMLtypes.length
	include("flow_measurement/flow_measurement.jl")
end


text:
# module for drawing moving robots
# add bias error of movement speed

module AnimeMoveSpeedBias
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))

  function animate_per_time(time_interval, world, bias_robot,
                            no_bias_robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(bias_robot)
    draw!(no_bias_robot)
  end

  function main(time_interval=0.1; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    init_pose = [0.0, 0.0, 0.0]
    circling = Agent(0.2, 10.0/180*pi)

    br = DifferentialWheeledRobot(init_pose, 0.2, "gray",
                                  circling, time_interval,
                                  bias_rate_stds=[0.2, 0.2])

    nbr = DifferentialWheeledRobot(init_pose, 0.2, "red",
                                   circling, time_interval,
                                   bias_rate_stds=[0.0, 0.0])

    if is_test == false
      anime = @animate for t in 0:time_interval:30
        animate_per_time(t, w, br, nbr)
      end
      path = "src/model/uncertainty/movement/speed_bias/anime_move_speed_bias.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
    else
      for t in 0:time_interval:10
        animate_per_time(t, w, br, nbr)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
#TODO change location in API docs
"""
    alleles(locus::T; miss::Bool = false) where T<:GenoArray
Return an array of all the non-missing alleles of a locus. Use
`miss = true` to include missing values.
"""
@inline function alleles(locus::T; miss::Bool = false) where T<:GenoArray
    if miss==true
        Base.Iterators.flatten(locus) |> collect
    else
        @inbounds Base.Iterators.flatten(locus[.!ismissing.(locus)]) |> collect
    end
end

"""
    unique_alleles(locus::T) where T<:GenoArray
Return an array of all the unique non-missing alleles of a locus.
"""
@inline function unique_alleles(locus::GenoArray)
    unique(alleles(locus))
end

"""
    convert_coord(coordinate::String)
Takes non-decimal-degree format as a `String` and returns it as a decimal degree
`Float32`. Can be broadcasted over an array of coordinate strings to convert them.
## Formatting requirements
- Coordinates as a `String` separated by spaces (`"11 43 41"`) or colons (`"11:43:41"`)
- Must use negative sign (`"-11 43.52"`) or single-letter cardinal direction (`"11 43.52W"`)
- Missing data should be coded as the string `"missing"` (can be accomplished with `replace!()`)
- Can mix colons and spaces (although it's bad practice)
### Example
```
julia> convert_coord("-41 31.52")
-41.5253f0

julia> convert_coord.(["-41 31.52", "25 11:54S"])
2-element Array{Float32,1}:
-41.5253
-25.1983
```
"""
function convert_coord(coordinate::String)
    lowercase(coordinate) == "missing" && return missing
    coord_strip = replace(uppercase(coordinate), r"[NSEW]" => "")
    split_coord = parse.(Float32, split(coord_strip, r"\s|:"))
    split_coord[2] /=60.0
    if length(split_coord) == 3
        split_coord[3] /=3600.0
    end
    conv = sum(abs.(split_coord))
    # N + E are positive | S + W are negative
    if split_coord[1] < 0 || occursin(r"[SW]", uppercase(coordinate))
        # negative
        return round(conv * -1, digits = 4)
    else
        # positive
        return round(conv, digits = 4)
    end
end

function Base.copy(data::PopData)
    PopData(copy.([data.meta,data.loci])...)
end

"""
    motivational_quote()
Randomly returns one of ~50 _somewhat_ motivational quotes.
"""
function motivational_quote()
# quotes retrieved on 06/11/2020 from https://www.coburgbanks.co.uk/blog/friday-funnies/50-funny-motivational-quotes/
    quotes =  [
    "\"The elevator to success is out of order. You’ll have to use the stairs, one step at a time.\" Joe Girard",
    "\"People often say that motivation doesn’t last. Well, neither does bathing – that’s why we recommend it daily.\" Zig Ziglar",
    "\"I always wanted to be somebody, but now I realise I should have been more specific.\" Lily Tomlin",
    "\"I am so clever that sometimes I don’t understand a single word of what I am saying.\" Oscar Wilde",
    "\"People say nothing is impossible, but I do nothing every day.\" Winnie the Pooh",
    "\"Life is like a sewer… what you get out of it depends on what you put into it.\" Tom Lehrer",
    "\"You can’t have everything. Where would you put it?\" Steven Wright",
    "\"Work until your bank account looks like a phone number.\" Unknown ",
    "\"Change is not a four letter word… but often your reaction to it is!\" Jeffrey Gitomer",
    "\"If you think you are too small to make a difference, try sleeping with a mosquito.\" Dalai Lama",
    "\"Bad decisions make good stories.\" Ellis Vidler",
    "\"I’ll probably never fully become what I wanted to be when I grew up, but that’s probably because I wanted to be a ninja princess.\" Cassandra Duffy",
    "\"When life gives you lemons, squirt someone in the eye.\" Cathy Guisewite",
    "\"A clear conscience is a sure sign of a bad memory.\" Mark Twain",
    "\"Well-behaved women seldom make history.\" Laurel Thatcher Ulrich",
    "\"I didn’t fail the test. I just found 100 ways to do it wrong.\" Benjamin Franklin",
    "\"I used to think I was indecisive, but now I’m not so sure.\" Unknown",
    "\"Don’t worry about the world coming to an end today. It’s already tomorrow in Australia.\" Charles Schulz",
    "\"Think like a proton. Always positive.\" Unknown",
    "\"Be happy – it drives people crazy.\" Unknown",
    "\"Optimist: someone who figures that taking a step backward after taking a step forward is not a disaster, it’s more like a cha-cha.\" Robert Brault",
    "\"The question isn’t who is going to let me, it’s who is going to stop me.\" Ayn Rand.",
    "\"You’re only given a little spark of madness. You mustn’t lose it.\" Robin Williams",
    "\"I am an early bird and a night owl… so I am wise and I have worms\" Michael Scott",
    "\"If you let your head get too big, it’ll break your neck.\" Elvis Presley",
    "\"The road to success is dotted with many tempting parking spaces.\" Will Rogers",
    "\"Live each day like it’s your second to the last. That way you can fall asleep at night.\" Jason Love",
    "\"Even a stopped clock is right twice every day. After some years, it can boast of a long series of successes.\" Marie Von Ebner-Eschenbach",
    "\"Honest criticism is hard to take, particularly from a relative, a friend, an acquaintance, or a stranger.\" Franklin P. Jones",
    "\"Opportunity is missed by most people because it is dressed in overalls and looks like work.\" Thomas Eddison",
    "\"A diamond is merely a lump of coal that did well under pressure.\" Unknown",
    "\"Nothing is impossible, the word itself says “I’m possible!\" Audrey Hepburn",
    "\"Friendship is like peeing on yourself: everyone can see it, but only you get the warm feeling that it brings.\" Robert Bloch",
    "\"Women who seek to be equal with men lack ambition.\" Marilyn Monroe",
    "\"By working faithfully eight hours a day you may eventually get to be boss and work twelve hours a day.\" Robert Frost",
    "\"The trouble with having an open mind, of course, is that people will insist on coming along and trying to put things in it.\" Terry Pratchett",
    "\"Age is of no importance unless you’re a cheese.\" Billie Burke",
    "\"When tempted to fight fire with fire, remember that the Fire Department usually uses water.\" Unknown",
    "\"Trying is the first step toward failure.\" Homer Simpson",
    "\"Happiness is just sadness that hasn’t happened yet.\" Unknown",
    "\"The best things in life are actually really expensive.\" Unknown",
    "\"A few harmless flakes working together can unleash an avalanche of destruction.\" Justin Sewell",
    "\"It could be that your purpose in life is to serve as a warning to others.\"  Ashleigh Brilliant",
    "\"If the world didn’t suck we’d all fly into space.\" Unknown",
    "\"Always remember that you are unique – just like everybody else.\" Unknown"
    ]
    return quotes[rand(1:length(quotes))]
end

"""
    multitest_missing(pvals::Vector{T}, correction::String) where T <: Union{Missing, <:AbstractFloat}
Modification to `MultipleTesting.adjust` to include `missing` values in the
returned array. Missing values are first removed from the array, the appropriate
correction made, then missing values are re-added to the array at their original
positions. See MultipleTesting.jl docs for full more detailed information.
#### example
`multitest_missing([0.1, 0.01, 0.005, 0.3], "bh")`

### `correction` methods (case insensitive)
- `"bonferroni"` : Bonferroni adjustment
- `"holm"` : Holm adjustment
- `"hochberg"` : Hochberg adjustment
- `"bh"` : Benjamini-Hochberg adjustment
- `"by"` : Benjamini-Yekutieli adjustment
- `"bl"` : Benjamini-Liu adjustment
- `"hommel"` : Hommel adjustment
- `"sidak"` : Šidák adjustment
- `"forwardstop"` or `"fs"` : Forward-Stop adjustment
- `"bc"` : Barber-Candès adjustment
"""
@inline function multitest_missing(pvals::Vector{T}, correction::String) where T <: Union{Missing, <:AbstractFloat}
    # get indices of where original missing are
    miss_idx = findall(i -> i === missing, pvals)

    # make seperate array for non-missing P vals
    p_no_miss = x[.!ismissing.(x)]

    # make a dict of all possible tests and their respective functions
    d = Dict(
        "bonferroni" => Bonferroni(),
        "holm" => Holm(),
        "hochberg" => Hochberg(),
        "bh" => BenjaminiHochberg(),
        "by" => BenjaminiYekutieli(),
        "bl" => BenjaminiLiu(),
        "hommel" => Hommel(),
        "sidak" => Sidak(),
        "forwardstop" => ForwardStop(),
        "fs" => ForwardStop(),
        "bc" => BarberCandes(),
    )

    correct = adjust(p_no_miss, d[lowercase(correction)]) |> Vector{Union{Missing, Float64}}

    # re-add missing to original positions
    @inbounds for i in miss_idx
        @inbounds insert!(correct, i, missing)
    end
    return correct
end

"""
    nonmissing(vec::T) where T<:AbstractArray
Convenience function to count the number of non-`missing` values
in a vector.
"""
function nonmissing(vec::T) where T<:AbstractArray
    count(!ismissing, vec)
end

"""
    quickstart()
Prints helpful text of how to get started using PopGen.
"""
function quickstart()
    printstyled("\n        Quickstart for PopGen\n\n", bold = true)
    println("Documentation: https://pdimens.github.io/PopGen.jl/")
    println("Motivational(?) quote: ", motivational_quote())
    println("\nA few things things you can do to get started:")

    printstyled("\nLoad in data\n\n", color = :magenta)
    println("- read_from(filename; kwargs...)")
    println("- genepop(infile; kwargs...)  or similar file-specific importer")
    println("- use available gulfsharks() or nancycats() datasets")

    printstyled("\nExplore PopData\n\n", color = :blue)
    println("- populations(PopData) to view population information")
    println("- loci(PopData) to view locus names")
    println("- samples(PopData) to view sample names")
    println("- missing(PopData, by = ...) to view missing information")

    printstyled("\nManipulate PopData\n\n", color = :light_red)
    println("- populations!(PopData, ...) to rename populations")
    println("- locations!(PopData, ...) to add geographical coordinates")
    println("- exclude!(PopData, kwargs...) to selectively remove data")

    printstyled("\nAnalyses\n\n", color = :green)
    println("- richness(PopData) to calculate allelic richness")
    println("- allele_avg(PopData) to calculate average # of alleles")
    println("- summary(PopData) to calculate F-statistics, heterozygosity, etc.")
    println("- hwe_test(PopData) to test for Hardy-Weinberg Equilibrium")

    return
end

"""
    reciprocal(num::T) where T <: Signed
Returns the reciprocal (1/number) of a number. Will return `0` when
the number is `0` instead of returning `Inf`.
"""
function reciprocal(num::T) where T <: Real
    !iszero(num) ? 1.0/float(num) : 0.0
end

#TODO add to API docs
"""
    strict_shuffle(x::T) where T <: AbstractArray
Shuffle only the non-missing values of a Vector and return a copy of the vector,
keeping the `missing` values at their original locations.
Use `strict_shuffle!` to edit in-place instead of returning a copy.
"""
@inline function strict_shuffle(x::T) where T <: AbstractArray
    # get indices of where original missing are
    miss_idx = findall(i -> i === missing, x)
    out_vec = shuffle(x[.!ismissing.(x)])

    insert!.(Ref(out_vec), miss_idx, missing)
    return out_vec
end

"""
    strict_shuffle!(x::T)! where T <: AbstractArray
Shuffle only the non-missing values of a Vector, keeping the
`missing` values at their original locations. Use `strict_shuffle`
to return a copy instead of editing in-place.
"""
function strict_shuffle!(x::T) where T <: AbstractArray
    @inbounds shuffle!(@view x[.!ismissing.(x)])
    return x
end


text:
mutable struct FVIntegrator{T1,mType,F,T2,tType}
  alg::T1
  mesh::mType
  Flux::F
  numvars::Int
  fluxes :: T2
  dt :: tType
  use_threads :: Bool
end

mutable struct FVDiffIntegrator{T1,mType,F,B,T2,tType}
  alg::T1
  mesh::mType
  Flux::F
  DiffMat::B
  M::Int
  fluxes::T2
  dt :: tType
  use_threads :: Bool
end

function compute_du!(du, fluxes, mesh, noscalar::Bool, ::Type{Val{true}})
    Threads.@threads for cell in cell_indices(mesh)
        @inbounds left = left_node_idx(cell, mesh)
        @inbounds right = right_node_idx(cell, mesh)
        if noscalar
            du[:,cell] = -(fluxes[:,right] - fluxes[:,left] ) / cell_volume(mesh.mesh, cell)
        else
            du[cell] = -(fluxes[right] - fluxes[left] ) / cell_volume(mesh.mesh, cell)
        end
    end
end

function compute_du!(du, fluxes, mesh, noscalar::Bool, ::Type{Val{false}})
    for cell in cell_indices(mesh)
        @inbounds left = left_node_idx(cell, mesh)
        @inbounds right = right_node_idx(cell, mesh)
        if noscalar
            du[:,cell] = -(fluxes[:,right] - fluxes[:,left] ) / cell_volume(mesh.mesh, cell)
        else
            du[cell] = -(fluxes[right] - fluxes[left] ) / cell_volume(mesh.mesh, cell)
        end
    end
end

function (fv::FVIntegrator)(du::AbstractArray{T}, u::AbstractArray{T}, p, t) where {T}
  mesh = fv.mesh; alg = fv.alg; Flux = fv.Flux; numvars = fv.numvars
  fluxes = fv.fluxes; dt = fv.dt; use_threads = fv.use_threads
  compute_fluxes!(fluxes, Flux, u, mesh, dt, alg, numvars > 1, Val{use_threads})
  apply_bc_in_fluxes!(fluxes, mesh)
  compute_du!(du, fluxes, mesh, numvars > 1, Val{use_threads})
  apply_bc_in_du!(du, mesh)
  nothing
end

function (fv::FVDiffIntegrator)(du::AbstractArray{T}, u::AbstractArray{T}, p, t) where {T}
  mesh = fv.mesh; alg = fv.alg; Flux = fv.Flux; M = fv.numvars
  fluxes = fv.fluxes; DiffMat = fv.DiffMat; dt = fv.dt; use_threads = fv.use_threads
  compute_Dfluxes!(fluxes, Flux, DiffMat, u, mesh, dt, alg, M>1, Val{use_threads})
  apply_bc_in_fluxes!(fluxes, mesh)
  compute_du!(du, fluxes, mesh, numvars > 1, Val{use_threads})
  apply_bc_in_du!(du, mesh)
  nothing
end


text:
function evolve_boundary!(
    KS::SolverSet,
    ctr::T1,
    a1face::T2,
    a2face::T2,
    dt,
    mode,
    bc,
) where {T1<:AbstractArray{ControlVolume2D2F,2},T2<:AbstractArray{Interface2D2F,2}}

    nx, ny, dx, dy = begin
        if KS.ps isa CSpace2D
            KS.ps.nr, KS.ps.nθ, KS.ps.dr, KS.ps.darc
        else
            KS.ps.nx, KS.ps.ny, KS.ps.dx, KS.ps.dy
        end
    end

    bcs = ifelse(bc isa Symbol, [bc, bc, bc, bc], bc)

    if bcs[1] == :maxwell
        @inbounds Threads.@threads for j = 1:ny
            vn = KS.vSpace.u .* a1face[1, j].n[1] .+ KS.vSpace.v .* a1face[1, j].n[2]
            vt = KS.vSpace.v .* a1face[1, j].n[1] .- KS.vSpace.u .* a1face[1, j].n[2]
            xc = (KS.ps.vertices[1, j, 1, 1] + KS.ps.vertices[1, j, 4, 1]) / 2
            yc = (KS.ps.vertices[1, j, 1, 2] + KS.ps.vertices[1, j, 4, 2]) / 2
            bcL = local_frame(KS.ib.bc(xc, yc), a1face[1, j].n[1], a1face[1, j].n[2])
            flux_boundary_maxwell!(
                a1face[1, j].fw,
                a1face[1, j].fh,
                a1face[1, j].fb,
                bcL, # left
                ctr[1, j].h,
                ctr[1, j].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dy[1, j],
                1,
            )
            a1face[1, j].fw .=
                global_frame(a1face[1, j].fw, a1face[1, j].n[1], a1face[1, j].n[2])
        end
    end

    if bcs[2] == :maxwell
        @inbounds Threads.@threads for j = 1:ny
            vn = KS.vSpace.u .* a1face[nx+1, j].n[1] .+ KS.vSpace.v .* a1face[nx+1, j].n[2]
            vt = KS.vSpace.v .* a1face[nx+1, j].n[1] .- KS.vSpace.u .* a1face[nx+1, j].n[2]
            xc = (KS.ps.vertices[nx, j, 2, 1] + KS.ps.vertices[nx, j, 3, 1]) / 2
            yc = (KS.ps.vertices[nx, j, 2, 2] + KS.ps.vertices[nx, j, 3, 2]) / 2
            bcR = local_frame(KS.ib.bc(xc, yc), a1face[nx+1, j].n[1], a1face[nx+1, j].n[2])
            flux_boundary_maxwell!(
                a1face[nx+1, j].fw,
                a1face[nx+1, j].fh,
                a1face[nx+1, j].fb,
                bcR, # right
                ctr[nx, j].h,
                ctr[nx, j].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dy[nx, j],
                -1,
            )
            a1face[nx+1, j].fw .=
                global_frame(a1face[nx+1, j].fw, a1face[nx+1, j].n[1], a1face[nx+1, j].n[2])
        end
    end

    if bcs[3] == :maxwell
        @inbounds Threads.@threads for i = 1:nx
            vn = KS.vSpace.u .* a2face[i, 1].n[1] .+ KS.vSpace.v .* a2face[i, 1].n[2]
            vt = KS.vSpace.v .* a2face[i, 1].n[1] .- KS.vSpace.u .* a2face[i, 1].n[2]
            xc = (KS.ps.vertices[i, 1, 1, 1] + KS.ps.vertices[i, 1, 2, 1]) / 2
            yc = (KS.ps.vertices[i, 1, 1, 2] + KS.ps.vertices[i, 1, 2, 2]) / 2
            bcD = local_frame(KS.ib.bc(xc, yc), a2face[i, 1].n[1], a2face[i, 1].n[2])

            flux_boundary_maxwell!(
                a2face[i, 1].fw,
                a2face[i, 1].fh,
                a2face[i, 1].fb,
                bcD, # left
                ctr[i, 1].h,
                ctr[i, 1].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dx[i, 1],
                1,
            )
            a2face[i, 1].fw .=
                global_frame(a2face[i, 1].fw, a2face[i, 1].n[1], a2face[i, 1].n[2])
        end
    end

    if bcs[4] == :maxwell
        @inbounds Threads.@threads for i = 1:nx
            vn = KS.vSpace.u .* a2face[i, ny+1].n[1] .+ KS.vSpace.v .* a2face[i, ny+1].n[2]
            vt = KS.vSpace.v .* a2face[i, ny+1].n[1] .- KS.vSpace.u .* a2face[i, ny+1].n[2]
            xc = (KS.ps.vertices[i, ny, 3, 1] + KS.ps.vertices[i, ny, 4, 1]) / 2
            yc = (KS.ps.vertices[i, ny, 3, 2] + KS.ps.vertices[i, ny, 4, 2]) / 2
            bcU = local_frame(KS.ib.bc(xc, yc), a2face[i, ny+1].n[1], a2face[i, ny+1].n[2])

            flux_boundary_maxwell!(
                a2face[i, ny+1].fw,
                a2face[i, ny+1].fh,
                a2face[i, ny+1].fb,
                bcU, # right
                ctr[i, ny].h,
                ctr[i, ny].b,
                vn,
                vt,
                KS.vSpace.weights,
                KS.gas.K,
                dt,
                dx[i, ny],
                -1,
            )
            a2face[i, ny+1].fw .=
                global_frame(a2face[i, ny+1].fw, a2face[i, ny+1].n[1], a2face[i, ny+1].n[2])
        end
    end

end


text:
using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

const optimizer = Powersense.Optimizer()
MOI.set(optimizer, MOI.RawParameter("external_optimizer"), GLPK.Optimizer)
MOI.set(optimizer, MOI.RawParameter("max_iter"), 3000)
MOI.set(optimizer, MOI.RawParameter("tol_residual"), 1.e-2)
MOI.set(optimizer, MOI.RawParameter("tol_infeas"), 1.e-2)
MOI.set(optimizer, MOI.RawParameter("OutputFlag"), 0)

const config_no_duals = MOIT.TestConfig(atol=1e-2, rtol=1e-2, duals=false, optimal_status=MOI.LOCALLY_SOLVED)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "Powersense"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_default_copy_to(optimizer, false)
    @test !MOIU.supports_default_copy_to(optimizer, true)
end

@testset "Unit ($algo)" for algo in ["Line Search", "Trust Region"]
# @testset "Unit ($algo)" for algo in ["Line Search"]
    MOI.set(optimizer, MOI.RawParameter("algorithm"), algo)
    bridged = MOIB.full_bridge_optimizer(optimizer, Float64)
    # A number of test cases are excluded because loadfromstring! works only
    # if the solver supports variable and constraint names.
    exclude = [
        "delete_variable", # Deleting not supported.
        "delete_variables", # Deleting not supported.
        "getvariable", # Variable names not supported.
        "solve_zero_one_with_bounds_1", # Variable names not supported.
        "solve_zero_one_with_bounds_2", # Variable names not supported.
        "solve_zero_one_with_bounds_3", # Variable names not supported.
        "getconstraint", # Constraint names not suported.
        "variablenames", # Variable names not supported.
        "solve_with_upperbound", # loadfromstring!
        "solve_with_lowerbound", # loadfromstring!
        "solve_integer_edge_cases", # loadfromstring!
        "solve_affine_lessthan", # loadfromstring!
        "solve_affine_greaterthan", # loadfromstring!
        "solve_affine_equalto", # loadfromstring!
        "solve_affine_interval", # loadfromstring!
        "get_objective_function", # Function getters not supported.
        "solve_constant_obj",  # loadfromstring!
        "solve_blank_obj", # loadfromstring!
        "solve_singlevariable_obj", # loadfromstring!
        "solve_objbound_edge_cases", # ObjectiveBound not supported.
        "solve_affine_deletion_edge_cases", # Deleting not supported.
        "solve_unbounded_model", # `NORM_LIMIT`
        "number_threads", # NumberOfThreads not supported
        "delete_nonnegative_variables", # get ConstraintFunction n/a.
        "update_dimension_nonnegative_variables", # get ConstraintFunction n/a.
        "delete_soc_variables", # VectorOfVar. in SOC not supported
        "solve_result_index", # DualObjectiveValue not supported
        "solve_farkas_interval_lower",
        "solve_farkas_lessthan",
        "solve_farkas_interval_upper",
        "solve_farkas_greaterthan",
        "solve_farkas_variable_lessthan_max",
        "solve_farkas_variable_lessthan",
        "solve_farkas_equalto_upper",
        "solve_farkas_equalto_lower",
    ]
    MOIT.unittest(bridged, config_no_duals, exclude)
    MOI.empty!(optimizer)
end

@testset "MOI Linear tests ($algo)" for algo in ["Line Search", "Trust Region"]
    MOI.set(optimizer, MOI.RawParameter("algorithm"), algo)
    exclude = ["linear8a", # Behavior in infeasible case doesn't match test.
               "linear12", # Same as above.
               "linear8b", # Behavior in unbounded case doesn't match test.
               "linear8c", # Same as above.
               "linear7",  # VectorAffineFunction not supported.
               "linear15", # VectorAffineFunction not supported.
               ]
    model_for_Powersense = MOIU.UniversalFallback(MOIU.Model{Float64}())
    linear_optimizer = MOI.Bridges.Constraint.SplitInterval{Float64}(
                         MOIU.CachingOptimizer(model_for_Powersense, optimizer))
    MOIT.contlineartest(linear_optimizer, config_no_duals, exclude)
    # Tests setting bounds of `SingleVariable` constraint
    # MOIT.linear1test(linear_optimizer, config_no_duals)
    MOI.empty!(optimizer)
end

# FIXME: LP subproblems are numerically instable in the trust region method.
@testset "MOI QP tests ($algo)" for algo in ["Line Search"]
    MOI.set(optimizer, MOI.RawParameter("algorithm"), algo)
    qp_optimizer = MOIU.CachingOptimizer(MOIU.Model{Float64}(), optimizer)
    MOIT.qptest(qp_optimizer, config_no_duals)
    MOI.empty!(optimizer)
end

@testset "MOI QCQP tests ($algo)" for algo in ["Line Search", "Trust Region"]
    MOI.set(optimizer, MOI.RawParameter("algorithm"), algo)
    qp_optimizer = MOIU.CachingOptimizer(MOIU.Model{Float64}(), optimizer)
    exclude = ["qcp1"] # VectorAffineFunction not supported.
    MOIT.qcptest(qp_optimizer, config_no_duals, exclude)
    MOI.empty!(optimizer)
end

# FIXME: LP subproblems are numerically instable in the trust region method.
@testset "MOI NLP tests ($algo)" for algo in ["Line Search"]
    MOI.set(optimizer, MOI.RawParameter("algorithm"), algo)
    MOIT.nlptest(optimizer, config_no_duals)
    # MOIT.hs071_test(optimizer, config_no_duals)
    MOI.empty!(optimizer)
end

#=@testset "Testing getters" begin
    MOI.Test.copytest(MOI.instantiate(Powersense.Optimizer, with_bridge_type=Float64), MOIU.Model{Float64}())
end=#

@testset "Bounds set twice" begin
    MOI.Test.set_lower_bound_twice(optimizer, Float64)
    MOI.Test.set_upper_bound_twice(optimizer, Float64)
end