text:
using Ghost
using Ghost: Tape, Input, Constant, mkcall, Variable, V

struct ONNXCtx
    name2var::Dict{String, Variable}
    backends::Vector{Symbol}
    exec::Bool
end

ONNXCtx(backends; exec=true) = ONNXCtx(Dict(), backends, exec)
ONNXCtx(;exec=true) = ONNXCtx(Dict(), [:ONNX], exec)

# TODO: implement rebind_context!()

"""
    getindex(tape::Tape{ONNXCtx}, onnx_name::String)

Get operation on the tape using the name in ONNX graph
"""
Base.getindex(tape::Tape{ONNXCtx}, onnx_name::String) =
    tape[tape.c.name2var[onnx_name]]

###############################################################################
#                               Operations                                    #
###############################################################################

"""
    push_call!(tape::Tape{ONNXCtx}, fn, args...; kwargs)

Shortcut for `push!(tape, mkcall(fn, args..))` also handling
keyword arguments and respecting `ONNXCtx.exec` setting.
"""
function push_call!(tape::Tape{ONNXCtx}, fn, args...; kwargs...)
    kwargs = NamedTuple(kwargs)
    if !isempty(kwargs)
        args = (kwargs, fn, args...)
        fn = Core.kwfunc(fn)
    end
    op = tape.c.exec ? mkcall(fn, args...) : mkcall(fn, args...; val=nothing)
    return push!(tape, op)
end

# A few constants to keep function signatures concise
struct OpConfig{BE, Op} end
const VarVec = Vector{Ghost.Variable}
const AttrDict = Dict{Symbol, Any}

function load_node!(tape::Tape, nd::NodeProto, backend::Symbol)
    args = [tape.c.name2var[name] for name in nd.input]
    attrs = convert(Dict{Symbol, Any}, Dict(nd.attribute))
    conf = OpConfig{backend, Symbol(nd.op_type)}()
    try
        out = load_node!(tape, conf, args, attrs)
        ismissing(out) && return out
        if out isa Tuple
            for i=1:length(nd.output)
                tape.c.name2var[nd.output[i]] = out[i]
            end
        else
            tape.c.name2var[nd.output[1]] = out
        end
    catch
        @error "Error while loading node $nd"
        rethrow()
    end
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Gemm}, args::VarVec, attrs::AttrDict)
    if (length(args) == 2 && get(attrs, :alpha, 1) == 1 &&
        get(attrs, :transA, 0) == 0 && get(attrs, :transB, 0) == 0)
        # simplified version: just matrix multiplication
        # note: arguments are swapped to account for row-major arrays
        return push_call!(tape, *, args[2], args[1])
    else
        # complete GEMM version
        kw = rename_keys(attrs, Dict(
            :transA => :tA,
            :transB => :tB,
            :alpha => :α,
            :beta => :β
        ))
        return push_call!(tape, onnx_gemm, args...; kw...)
    end
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Conv}, args::VarVec, attrs::AttrDict)
    kw = from_onnx_conv(attrs) |> NamedTuple
    return push_call!(tape, conv, args...; kw...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :MaxPool}, args::VarVec, attrs::AttrDict)
    kw = from_onnx_conv(attrs; pooling=true) |> NamedTuple
    return push_call!(tape, maxpool, args[1]; kw...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :GlobalAveragePool}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, global_average_pool, args...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Flatten}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, onnx_flatten, args...; attrs...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Add}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, add, args...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Mul}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, mul, args...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Relu}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, relu, args[1])
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :MatMul}, args::VarVec, attrs::AttrDict)
    A_ndims = ndims(args[1]._op.val)
    B_ndims = ndims(args[2]._op.val)
    if A_ndims == 2 && B_ndims == 2
        return push_call!(tape, *, args[2], args[1])
    elseif A_ndims in (2, 3) && B_ndims in (2, 3)
        return push_call!(tape, NNlib.batched_mul, args[2], args[1])
    else
        error("MatMul with arrays of $A_ndims and $B_ndims is not implemented yet")
    end
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Sigmoid}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, Broadcast.broadcast, NNlib.sigmoid, args...)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :BatchNormalization},
        args::VarVec, attrs::AttrDict)
    kw = from_onnx_norm(attrs)
    bn = push_call!(tape, batch_norm, args...; kw...)
    if bn._op.val isa Tuple
        # usual in training mode
        # unpack tuples into calls to getfield
        y = push_call!(tape, getfield, bn, 1)
        μnext = push_call!(tape, getfield, bn, 2)
        σ²next = push_call!(tape, getfield, bn, 3)
        return y, μnext, σ²next
    else
        return bn
    end
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Shape}, args::VarVec, attrs::AttrDict)
    # TODO: handle start and end attributes
    return push_call!(tape, size_vector, args[1])
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Constant}, args::VarVec, attrs::AttrDict)
    val_attr = first(keys(attrs))
    val = if val_attr == :value
        array(attrs[val_attr])
    else
        error("Don't know how to load constant value from attribute $val_attr")
    end
    return push!(tape, Constant(val))
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Gather}, args::VarVec, attrs::AttrDict)
    axis = get(attrs, :axis, 0)
    data = tape[args[1]].val
    dim = ndims(data) - axis
    return push_call!(tape, onnx_gather, args...; dim=dim)
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Unsqueeze}, args::VarVec, attrs::AttrDict)
    if length(args) == 2
        # ONNX >= v13
        return push_call!(tape, onnx_unsqueeze, args...)
    elseif length(args) == 1
        # ONNX < v13
        axes = attrs[:axes]
        v_axes = push!(tape, Constant(axes))
        return push_call!(tape, onnx_unsqueeze, args[1], v_axes)
    else
        throw(ArgumentError("Cannot load node from Unsqueeze with $(length(args)) arguments"))
    end
end

function load_node!(tape::Tape, ::OpConfig{:ONNX, :Slice}, args::VarVec, attrs::AttrDict)
    return push_call!(tape, onnx_slice, args...)
end

###############################################################################
#                                    API                                      #
###############################################################################

"""
    load(io::IO, model_args...; backends=[:ONNX], exec::Bool=true)
    load(filename::String, model_args...; backends=[:ONNX], exec::Bool=true)

Load an ONNX model as a Ghost.Tape. The way a particular ONNX node is deserialized is
controlled by methods of [load_node!](@ref) dispatched by backend and node's op_type.

`backends` parameter can be used to customize the loading process.

`exec` parameter instructs the loader to execute every added operation just after
the addition, making the debugging easier. Default is `true`.

See also: [`save!`](@ref)
"""
function load(io::IO, args...; backends=[:ONNX], exec::Bool=true)
    onnx_model = readproto(io, ModelProto());
    g = onnx_model.graph;
    tape = Tape(ONNXCtx(backends; exec=exec))
    # create map of initializers
    init_vals = Dict{String, Any}(init.name => array(init)
        for init in g.initializer)
    # load inputs; if input has init value, take it
    # otherwise take the next available argument value
    arg_idx = 1
    used_init_names = Set([])
    for inp in g.input
        val = get(init_vals, inp.name, missing)
        v = V(0)   # will be overwritten
        if val === missing && exec == true
            @assert(
                arg_idx <= length(args),
                "Neither initializer, nor argument is provided for input $(inp.name)"
            )
            val = args[arg_idx]
            arg_idx += 1
            v = push!(tape, Input(val))
        else
            # convert inputs that also have initializers to constants
            # these are usually model parameters, but may
            v = push!(tape, Constant(val))
        end
        tape.c.name2var[inp.name] = v
        push!(used_init_names, inp.name)
    end
    # load the rest of initilizers as constants
    for init in g.initializer
        name = init.name
        if !in(name, used_init_names)
            val = init_vals[name]
            v = push!(tape, Ghost.Constant(val))
            tape.c.name2var[name] = v
        end
    end
    # load nodes
    for nd in g.node
        success = false
        for backend in tape.c.backends
          if !ismissing(load_node!(tape, nd, backend))
            success = true
            @debug "Loaded $(nd.op_type) using backend $(backend)"
            break
          end
        end
        success || error("Couldn't load node for $(nd.op_type), " *
                         "tried the following backends: $(tape.c.backends)")
    end
    if length(g.output) == 1
        tape.result = Ghost.bound(tape, V(length(tape)))
    else
        # tuple output: we expect tape to contain these outputs as vars  destructured
        # from a multi-ouput op using a sequence of `getfield()` calls
        vars = [tape.c.name2var[o.name] for o in g.output]
        @assert(all(tape[v] isa Call && tape[v].fn == getfield for v in vars),
            "Don't understand this multi-output result of the graph")
        tape.result = tape[vars[1]].args[1]
    end
    return tape
end

function load(filename::String, args...; backends=[:ONNX], exec::Bool=true)
    return open(filename) do io
        load(io, args...; backends=backends, exec=exec)
    end
end


text:
# Probabilistic Principal Component Analysis

"""
This type contains probabilistic PCA model parameters.
"""
struct PPCA{T<:Real} <: LatentVariableDimensionalityReduction
    mean::Vector{T}       # sample mean: of length d (mean can be empty, which indicates zero mean)
    W::Matrix{T}          # weight matrix: of size d x p
    σ²::T                 # residual variance
end

## properties
"""
    size(M::PPCA)

Returns a tuple with values of the input dimension ``d``, *i.e* the dimension of
the observation space, and the output dimension ``p``, *i.e* the dimension of
the principal subspace.
"""
size(M::PPCA) = size(M.W)

"""
    mean(M::PPCA)

Get the mean vector (of length ``d``).
"""
mean(M::PPCA) = fullmean(size(M)[1], M.mean)

"""
    projection(M::PPCA)

Returns the projection matrix (of size ``(d, p)``). Each column of the projection
matrix corresponds to a principal component.

The principal components are arranged in descending order of the corresponding variances.

"""
projection(M::PPCA) = svd(M.W).U # recover principle components from the weight matrix

"""
    var(M::PPCA)

Returns the total residual variance of the model `M`.
"""
var(M::PPCA) = M.σ²

"""
    loadings(M::PPCA)

Returns the factor loadings matrix (of size ``(d, p)``) of the model `M`.
"""
loadings(M::PPCA) = M.W

"""
    cov(M::PPCA)

Returns the covariance of the model `M`.
"""
cov(M::PPCA) = M.W'*M.W + M.σ²*I

## use
"""
    predict(M::PPCA, x)

Transform observations `x` into latent variables. Here, `x` can be either a vector
of length `d` or a matrix where each column is an observation.
"""
function predict(M::PPCA, x::AbstractVecOrMat{T}) where {T<:Real}
    xn = centralize(x, M.mean)
    n = size(M)[2]
    return inv(cov(M))*M.W'*xn
end

"""
    reconstruct(M::PPCA, z)

Approximately reconstruct observations from the latent variable given in `z`.
Here, `z` can be either a vector of length `p` or a matrix where each column gives
the latent variables for an observation.
"""
function reconstruct(M::PPCA, z::AbstractVecOrMat{T}) where {T<:Real}
    W  = M.W
    WTW = W'W
    n = size(M)[2]
    C  = WTW + var(M) * I
    return W*inv(WTW)*C*z .+ mean(M)
end

## show

function Base.show(io::IO, M::PPCA)
    i, o = size(M)
    print(io, "Probabilistic PCA(indim = $i, outdim = $o, σ² = $(var(M)))")
end

## core algorithms
"""
    ppcaml(Z, mean; ...)

Compute probabilistic PCA using on maximum likelihood formulation for a centralized
sample matrix `Z`.

*Parameters*:
- `Z`: a centralized samples matrix
- `mean`: The mean vector of the **original** samples, which can be a vector of
length `d`, or an empty vector indicating a zero mean.

Returns the resultant [`PPCA`](@ref) model.

**Note:** This function accepts two keyword arguments: `maxoutdim` and `tol`.
"""
function ppcaml(Z::AbstractMatrix{T}, mean::Vector{T};
                tol::Real=1.0e-6, # convergence tolerance
                maxoutdim::Int=size(Z,1)-1) where {T<:Real}

    check_pcaparams(size(Z,1), mean, maxoutdim, 1.)

    d, n = size(Z)

    # SVD decomposition
    Svd = svd(Z)
    λ = Svd.S
    ord = sortperm(λ; rev=true)
    V = abs2.(λ[ord]) ./ (n-1)

    # filter 0 eigenvalues and adjust number of latent dimensions
    idxs = findall(V .< tol)
    l = length(idxs)
    l = l == 0 ? maxoutdim : l

    # variance "loss" in the projection
    σ² = sum(V[l+1:end])/(d-l)

    @inbounds for i in 1:l
        V[i] = sqrt(V[i] - σ²)
    end

    W = Svd.U[:,ord[1:l]]*diagm(0 => V[1:l])

    return PPCA(mean, W, σ²)
end

"""
    ppcaem(S, mean, n; ...)

Compute probabilistic PCA based on expectation-maximization algorithm for a given sample covariance matrix `S`.

*Parameters*:
- `S`: The sample covariance matrix.
- `mean`: The mean vector of original samples, which can be a vector of length `d`,
or an empty vector indicating a zero mean.
- `n`: The number of observations.

Returns the resultant [`PPCA`](@ref) model.

**Note:** This function accepts three keyword arguments: `maxoutdim`, `tol`, and `maxiter`.
"""
function ppcaem(S::AbstractMatrix{T}, mean::Vector{T}, n::Int;
                maxoutdim::Int=size(S,1)-1,
                tol::Real=1.0e-6,   # convergence tolerance
                maxiter::Integer=1000) where {T<:Real}

    check_pcaparams(size(S,1), mean, maxoutdim, 1.)

    d = size(S,1)
    q = maxoutdim
    Iq = Matrix{T}(I, q, q)
    Id = Matrix{T}(I, d, d)
    W = Matrix{T}(I, d, q)
    σ² = zero(T)
    M⁻¹ = inv(W'W .+ σ² * Iq)

    i = 1
    L_old = 0.
    chg = NaN
    converged = false
    while i < maxiter
        # EM-steps
        SW = S*W
        W⁺  = SW*inv(σ²*Iq + M⁻¹*W'*SW)
        σ²⁺ = tr(S - SW*M⁻¹*(W⁺)')/d
        # new parameters
        W = W⁺
        σ² = σ²⁺
        # log likelihood
        C = W*W'.+ σ²*Id
        M⁻¹ = inv(W'*W .+ σ²*Iq)
        C⁻¹ = (Id - W*M⁻¹*W')/σ²
        L = (-n/2)*(log(det(C)) + tr(C⁻¹*S))  # (-n/2)*d*log(2π) omitted

        @debug "Likelihood" iter=i L=L ΔL=abs(L_old - L)
        chg = abs(L_old - L)
        if chg < tol
            converged = true
            break
        end
        L_old = L
        i += 1
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))

    return PPCA(mean, W, σ²)
end

"""
    bayespca(S, mean, n; ...)

Compute probabilistic PCA using a Bayesian algorithm for a given sample covariance matrix `S`.

*Parameters*:
- `S`: The sample covariance matrix.
- `mean`: The mean vector of original samples, which can be a vector of length `d`,
or an empty vector indicating a zero mean.
- `n`: The number of observations.

Returns the resultant [`PPCA`](@ref) model.

**Notes:**
- This function accepts three keyword arguments: `maxoutdim`, `tol`, and `maxiter`.
- Function uses the `maxoutdim` parameter as an upper boundary when it automatically
determines the latent space dimensionality.
"""
function bayespca(S::AbstractMatrix{T}, mean::Vector{T}, n::Int;
                 maxoutdim::Int=size(S,1)-1,
                 tol::Real=1.0e-6,   # convergence tolerance
                 maxiter::Integer=1000) where {T<:Real}

    check_pcaparams(size(S,1), mean, maxoutdim, 1.)

    d = size(S,1)
    q = maxoutdim
    Iq = Matrix{T}(I, q, q)
    Id = Matrix{T}(I, d, d)
    W = Matrix{T}(I, d, q)
    wnorm = zeros(T, q)
    σ² = zero(T)
    M = W'*W .+ σ²*Iq
    M⁻¹ = inv(M)
    α = zeros(T, q)

    i = 1
    chg = NaN
    L_old = 0.
    converged = false
    while i < maxiter
        # EM-steps
        SW = S*W
        W⁺  = SW*inv(σ²*(Iq+diagm(0=>α)*M/n) + M⁻¹*W'*SW)
        σ²⁺ = tr(S - SW*M⁻¹*(W⁺)')/d
        # new parameters
        W = W⁺
        σ² = σ²⁺
        @inbounds for j in 1:q
            wnorm[j] = norm(W[:,j])^2
            α[j] = wnorm[j] < eps() ? maxintfloat(Float64) : d/wnorm[j]
        end
        # log likelihood
        C = W*W'.+ σ²*Id
        M = W'*W .+ σ²*Iq
        M⁻¹ = inv(M)
        C⁻¹ = (Id - W*M⁻¹*W')/σ²
        L = (-n/2)*(log(det(C)) + tr(C⁻¹*S))  # (-n/2)*d*log(2π) omitted

        @debug "Likelihood" iter=i L=L ΔL=abs(L_old - L)
        chg = abs(L_old - L)
        if chg < tol
            converged = true
            break
        end
        L_old = L
        i += 1
    end
    converged || throw(ConvergenceException(maxiter, chg, oftype(chg, tol)))

    return PPCA(mean, W[:,wnorm .> 0.], σ²)
end

## interface functions
"""

    fit(PPCA, X; ...)

Perform probabilistic PCA over the data given in a matrix `X`.
Each column of `X` is an observation. This method returns an instance of [`PPCA`](@ref).

**Keyword arguments:**

Let `(d, n) = size(X)` be respectively the input dimension and the number of observations:

- `method`: The choice of methods:
    - `:ml`: use maximum likelihood version of probabilistic PCA (*default*)
    - `:em`: use EM version of probabilistic PCA
    - `:bayes`: use Bayesian PCA
- `maxoutdim`: Maximum output dimension (*default* `d-1`)
- `mean`: The mean vector, which can be either of:
    - `0`: the input data has already been centralized
    - `nothing`: this function will compute the mean (*default*)
    - a pre-computed mean vector
- `tol`: Convergence tolerance (*default* `1.0e-6`)
- `maxiter`: Maximum number of iterations (*default* `1000`)

**Notes:** This function calls [`ppcaml`](@ref), [`ppcaem`](@ref) or
[`bayespca`](@ref) internally, depending on the choice of method.
"""
function fit(::Type{PPCA}, X::AbstractMatrix{T};
             method::Symbol=:ml,
             maxoutdim::Int=size(X,1)-1,
             mean=nothing,
             tol::Real=1.0e-6,   # convergence tolerance
             maxiter::Integer=1000) where {T<:Real}

    @assert !SparseArrays.issparse(X) "Use Kernel PCA for sparse arrays"

    d, n = size(X)

    # process mean
    mv = preprocess_mean(X, mean)
    if !(isempty(mv) || length(mv) == d)
        throw(DimensionMismatch("Dimensions of weight matrix and mean are inconsistent."))
    end

    if method == :ml
        Z = centralize(X, mv)
        M = ppcaml(Z, mv, maxoutdim=maxoutdim, tol=tol)
    elseif method == :em || method == :bayes
        S = covm(X, isempty(mv) ? 0 : mv, 2)
        if method == :em
            M = ppcaem(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter)
        elseif method == :bayes
            M = bayespca(S, mv, n, maxoutdim=maxoutdim, tol=tol, maxiter=maxiter)
        end
    else
        throw(ArgumentError("Invalid method name $(method)"))
    end

    return M::PPCA
end

ext: jl


text:
using Distributions, MCMC

mcmodel = model(MvNormal([1. 0.1; 0.1 1]), init=[1., 2.])

mcchain = run(mcmodel, HMC(0.75), SerialMC(steps=10000, burnin=1000))

acceptance(mcchain)

describe(mcchain)


text:
# Random function sample generation
# Following the Random API: https://docs.julialang.org/en/v1/stdlib/Random/#Hooking-into-the-Random-API
"""
    BLRFunctionSample{Tw,Tϕ}

A function sampled from
 `b::Union{BayesianLinearRegressor,BasisFunctionRegressor}` by taking a fixed
 weight sample `w ~ p(w)`. `ϕ` is the feature mapping if sampled from a
 `BasisFunctionRegressor`, or is the identity function if sampled from a
 `BayesianLinearRegressor`.
"""
struct BLRFunctionSample{Tw<:AbstractVector,Tϕ}
    w::Tw
    ϕ::Tϕ
end

(s::BLRFunctionSample)(X::AbstractMatrix{<:Real}) = s.ϕ(X)'s.w
(s::BLRFunctionSample)(X::ColVecs) = s.ϕ(X).X's.w
(s::BLRFunctionSample)(X::RowVecs) = s.ϕ(X).X * s.w

BLRorBasisFunction = Union{BayesianLinearRegressor,BasisFunctionRegressor}

function Random.Sampler(::Type{<:AbstractRNG}, blr::BLRorBasisFunction, ::Random.Repetition)
    return blr
end

function Random.rand(rng::AbstractRNG, b::BLRorBasisFunction)
    blr, ϕ = _blr_and_mapping(b)
    w = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, size(blr.mw))
    return BLRFunctionSample(w, ϕ)
end

function Random.rand(rng::AbstractRNG, b::BLRorBasisFunction, dims::Dims)
    blr, ϕ = _blr_and_mapping(b)
    ws = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, (only(size(blr.mw)), prod(dims)))
    bs = [BLRFunctionSample(w, ϕ) for w in eachcol(ws)]
    return reshape(bs, dims)
end

function Random.rand!(
    rng::AbstractRNG, A::AbstractArray{<:BLRFunctionSample}, b::BLRorBasisFunction
)
    blr, ϕ = _blr_and_mapping(b)
    ws = blr.mw .+ _cholesky(blr.Λw).U \ randn(rng, (only(size(blr.mw)), prod(size(A))))
    for i in LinearIndices(A)
        @inbounds A[i] = BLRFunctionSample(ws[:, i], ϕ)
    end
    return A
end

_blr_and_mapping(b::BayesianLinearRegressor) = b, identity
_blr_and_mapping(b::BasisFunctionRegressor) = b.blr, b.ϕ


text:
using Documenter, ApproxFun, BandedMatrices, DomainSets, LinearAlgebra

makedocs(
			doctest = false,
			clean = true,
			format = Documenter.HTML(),
			sitename = "ApproxFun.jl",
			authors = "Sheehan Olver",
			pages = Any[
					"Home" => "index.md",
					"Usage" => Any[
						"Domains" => "usage/domains.md",
						"Spaces" => "usage/spaces.md",
						"Constructors" => "usage/constructors.md",
						"Operators" => "usage/operators.md",
						"Linear Equations" => "usage/equations.md"
					],
					"Frequently Asked Questions" => "faq.md",
					"Library" => "library.md"
				]
			)

deploydocs(
    repo   = "github.com/JuliaApproximation/ApproxFun.jl.git"
    )


text:
using BioServices.UMLS
using BioMedQuery.PubMed
using BioMedQuery.DBUtils
using DataFrames

"""
    map_mesh_to_umls_async!(db, c::Credentials; timeout, append_results, verbose)

Build (using async UMLS-API calls) and store in the given database a map from
MESH descriptors to UMLS Semantic Concepts. For large queies this function will
be faster than it's synchrounous counterpart

## Arguments

* `db`: Database. Must contain TABLE:mesh_descriptor. For each of the descriptors  in that table, search and insert the associated semantic concepts into a new (cleared) TABLE:mesh2umls
* `user` : UMLS username
* `psswd` : UMLS Password
* `append_results::Bool` : If false a NEW and EMPTY mesh2umls database table in creted
* `batch_size`: Number of
"""
function map_mesh_to_umls_async!(db, user, psswd; timeout = 5, append_results=false, verbose=false, batch_size=50)

    # Determine engine
    sql_engine = (typeof(db)== MySQL.Connection) ? MySQL : SQLite
    engine_info = (sql_engine == MySQL) ? "ENGINE=InnoDB DEFAULT CHARSET=utf8" : ""

    #if the mesh2umls relationship table doesn't esxist, create it
    sql_engine.execute!(db, "CREATE table IF NOT EXISTS mesh2umls (
                                mesh VARCHAR(255),
                                umls VARCHAR(255),
                                PRIMARY KEY(mesh, umls)
                            ) $engine_info;")

    #clear the relationship table
    if !append_results
        sql_engine.execute!(db, "DELETE FROM mesh2umls")
    end

    #select all mesh descriptors
    mq = sql_engine.query(db,"SELECT name FROM mesh_desc;")

    #get the array of terms
    mesh_terms = mq[1]
    println("----------Matching MESH to UMLS-----------")
    println(mesh_terms)

    tgt = get_tgt(username = user, password = psswd)
    errors = 200*ones(length(mesh_terms))
    times = -ones(length(mesh_terms))

    for m=1:batch_size:length(mesh_terms)
        end_loop=m+batch_size
        if end_loop > length(mesh_terms)
            end_loop = length(mesh_terms)
        end
        @sync for i=m:end_loop
            #submit umls async batch query
            @async begin
                term = mesh_terms[i]
                query = Dict("string"=>term, "searchType"=>"exact" )
                # println("term: ", term)
                for attempt=1:5
                    try
                        all_results= search_umls(tgt, query, timeout=timeout)
                        @info "Descriptor $i out of $(length(mesh_terms)): $term"
                        if length(all_results) > 0
                            cui = best_match_cui(all_results)
                            if cui == ""
                                println("Nothing!")
                                println(all_results)
                            end
                            all_concepts = get_semantic_types(tgt, cui)
                            for concept in all_concepts
                                insert_row!(db, "mesh2umls", Dict(:mesh=> term, :umls=> concept), verbose)
                            end
                        end
                        break
                    catch err
                        @warn "! failed attempt $attempt out of 5 for term $term with error " err
                    end
                end
            end
        end
    end

end

"""
    map_mesh_to_umls_async(mesh_df, user, psswd; timeout, append_results, verbose)

Build (using async UMLS-API calls) and return a map from
MESH descriptors to UMLS Semantic Concepts. For large queies this function will
be faster than it's synchrounous counterpart

## Arguments

* `mesh_df`: DataFrame countaining Mesh_Descriptors. This is the dataframe with the key `mesh_desc` that is returned from pubmed_search_and_parse.
* `user` : UMLS username
* `psswd` : UMLS Password
"""
function map_mesh_to_umls_async(mesh_df::DataFrame, user, psswd; timeout = 5, verbose=false)

    mesh = Vector{String}()
    umls = Vector{String}()

    #get the array of terms
    mesh_terms = mesh_df[:name]
    num_mesh = length(mesh_terms)

    println("----------Matching MESH to UMLS-----------")
    println(mesh_terms)

    tgt = get_tgt(username = user, password = psswd)
    errors = 200*ones(num_mesh)
    times = -ones(num_mesh)
    batch_size = 50

    for m=1:batch_size:num_mesh

        end_loop = m + batch_size > num_mesh ? num_mesh : m + batch_size

        @sync for i=m:end_loop
            #submit umls async batch query
            @async begin
                term = mesh_terms[i]
                query = Dict("string"=>term, "searchType"=>"exact" )
                # println("term: ", term)
                for attempt=1:5
                    try
                        all_results= search_umls(tgt, query, timeout=timeout)
                        @info "Descriptor $i, out of $num_mesh: $term"
                        if length(all_results) > 0
                            cui = best_match_cui(all_results)
                            if cui == ""
                                println("Nothing!")
                                println(all_results)
                            end
                            all_concepts = get_semantic_types(tgt, cui)
                            for concept in all_concepts
                                push!(mesh, term)
                                push!(umls, concept)
                            end
                        end
                        break
                    catch err
                        @warn "! failed attempt $attempt out of 5 for term $term with error "  err
                    end
                end
            end
        end
    end

    return sort(unique(DataFrame(descriptor = mesh, concept = umls)))

end


text:
nodes = [1 0;
         1 1;
         0 1;
         0 0;
         .5 .5;
         1 .5;
         .5 1;
         0 .5;
         .5 0;
         .75 .25;
         .75 .75;
         .25 .75;
         .25 .25]

elem = [1 4 5 9 13 10;
        1 2 5 6 11 10;
        2 3 5 7 12 11;
        3 4 5 8 13 12]
edge = [4 9 1;
        1 6 2;
        2 7 3;
        3 8 4]
bc = [2 -1 4 -1]'
f(x) = [[1e6 0]',[0 0]']

max_stars_repo_name: billyon/first

max_issues_repo_name: billyon/first

max_forks_repo_name: billyon/first


text:
# Types for Julia implementation of pseudolikelihood estimation (Stage2)
# Claudia August 2014
# in julia: include("types.jl")
#
# Types: Edge, Node, HybridNetwork
################################################################

# procedure to create hybrid network:
# 1) create edges defined as hybrid or not
# 2) create nodes with such edges
# 3) setNode! to add nodes into edges
# 4) create hybrid network
# 5) updateGammaz!

# -------------- EDGE -------------------------#
"""
 `ANode`

Abstract node. An object of type [`Edge`](@ref) has a `node` attribute,
which is an vector of 2 `ANode` objects.
The object of type [`Node`](@ref) is an `ANode`, and has an `edge` attribute,
which is vector of `Edge` objects.
"""
abstract type ANode end

"""
    Edge(number)

Data structure for an edge and its various attributes. Most notably:
- `number` (integer): serves as unique identifier;
  remains unchanged when the network is modified,
  with a nearest neighbor interchange for example
- `node`: a vector of [`Node`]s, normally just 2 of them
- `isChild1` (boolean): `true` if `node[1]` is the child node of the edge,
  false if `node[1]` is the parent node of the edge
- `length`: branch length
- `hybrid` (boolean): whether the edge is a tree edge or a hybrid edge
  (in which case `isChild1` is important, even if the network is semi-directed)
- `gamma`: proportion of genetic material inherited by the child node via the edge;
  1.0 for a tree edge
- `isMajor` (boolean): whether the edge is the major path to the child node;
  `true` for tree edges, since a tree edge is the only path to its child node;
  normally true if `gamma>0.5`.

and other fields, used very internally
"""
mutable struct Edge
    number::Int
    length::Float64 #default 1.0
    hybrid::Bool
    y::Float64 # exp(-t), cannot set in constructor for congruence
    z::Float64 # 1-y , cannot set in constructor for congruence
    gamma::Float64 # set to 1.0 for tree edges, hybrid?gamma:1.0
    node::Array{ANode,1}
    isChild1::Bool # used for hybrid edges to set the direction (default true)
    isMajor::Bool  # major edge treated as tree edge for network traversal
                   # true if gamma>.5, or if it is the original tree edge
    inCycle::Int # = Hybrid node number if this edge is part of a cycle created by such hybrid node
    # -1 if not part of cycle. used to add new hybrid edge. updated after edge is part of a network
    containRoot::Bool # true if this edge can contain a root given the direction of hybrid edges
                      # used to add new hybrid edge. updated after edge is part of a network
    istIdentifiable::Bool # true if the parameter t (length) for this edge is identifiable as part of a network
                          # updated after part of a network
    fromBadDiamondI::Bool # true if the edge came from deleting a bad diamond I hybridization
    # inner constructors: ensure congruence among (length, y, z) and (gamma, hybrid, isMajor), and size(node)=2
    Edge(number::Int) = new(number,1.0,false,exp(-1.0),1-exp(-1.0),1.,[],true,true,-1,true,true,false)
    Edge(number::Int, length::Float64) = new(number,length,false,exp(-length),1-exp(-length),1.,[],true,true,-1,true,true,false)
    Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64) =
        new(number,length,hybrid,exp(-length),1-exp(-length), hybrid ? gamma : 1.,[],true, !hybrid || gamma>0.5 ,-1,!hybrid,true,false)
    Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,isMajor::Bool) =
        new(number,length,hybrid,exp(-length),1-exp(-length), hybrid ? gamma : 1.,[],true,isMajor,-1,!hybrid,true,false)
    function Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,node::Array{ANode,1})
        size(node,1) != 2 ?
        error("vector of nodes must have exactly 2 values") :
        new(number,length,hybrid,exp(-length),1-exp(-length),
            hybrid ? gamma : 1., node,true, !hybrid || gamma>0.5,
            -1,!hybrid,true,false)
    end
    function Edge(number::Int, length::Float64,hybrid::Bool,gamma::Float64,node::Array{ANode,1},isChild1::Bool, inCycle::Int, containRoot::Bool, istIdentifiable::Bool)
        size(node,1) != 2 ?
        error("vector of nodes must have exactly 2 values") :
        new(number,length,hybrid,exp(-length),1-exp(-length),
            hybrid ? gamma : 1., node,isChild1, !hybrid || gamma>0.5,
            inCycle,containRoot,istIdentifiable,false)
    end
end

# warning: gammaz, inCycle, isBadTriangle/Diamond updated until the node is part of a network
"""
    Node(number, leaf)
    Node(number, leaf, hybrid)

Data structure for a node and its various attributes. Most notably:

- `number` (integer): serves as unique identifier;
  remains unchanged when the network is modified,
  with a nearest neighbor interchange for example
- `leaf` (boolean): whether the node is a leaf (with data typically) or an
  internal node (no data typically)
- `name` (string): taxon name for leaves; internal node may or may not have a name
- `edge`: vector of [`Edge`](@ref)s that the node is attached to;
  1 if the node is a leaf, 2 if the node is the root, 3 otherwise, and
  potentially more if the node has a polytomy
- `hybrid` (boolean): whether the node is a hybrid node (with 2 or more parents)
  or a tree node (with a single parent)

Other more internal attributes include:

- `isBadDiamondI` and `isBadDiamondII` (booleans): whether the node is a
  hybrid node where the reticulation forms a cycle of 4 nodes (diamond),
  and where both parents of the hybrid nodes are connected to a leaf.
  In a bad diamond of type I, the hybrid node itself is also connected
  to a leaf but the common neighbor of the 2 hybrid's parents is not connected
  to a leaf.
  In a bad diamond of type II, the hybrid node has an internal node as child,
  and the common neighbor of the 2 hybrid's parents is connected to a leaf.
- `isBadTriangle`, `isVeryBadTriangle` and `isExtBadTriangle` (booleans):
  true if the reticulation forms a cycle of 3 nodes (triangle) and
  depending on the number of leaves attached these 3 nodes. The triangle means
  that the 2 parents of the hybrid node are directly related:
  one is the child of the other. `isBadTriangle` is true if the triangle is
  "good", as per Solís-Lemus & Ané (2016), that is, if all 3 nodes in the cycle
  are not connected to any leaves (the reticulation is detectable from quartet
  concordance factors, even though all branch lengths are not identifiable).
  `isVeryBadTriangle` is true if 2 (or all) of the 3 nodes are connected to a
  leaf, in which case the reticulation is undetectable from unrooted gene tree
  topologies (thus it's best to exclude these reticulations from a search).
  `isBadTriangle` is true if exactly 1 of the 3 nodes is connected to a leaf.

For details see Solís-Lemus & Ané (2016, doi:10.1371/journal.pgen.1005896)
"""
mutable struct Node <: ANode
    number::Int
    leaf::Bool
    hybrid::Bool
    gammaz::Float64  # notes file for explanation. gammaz if tree node, gamma2z if hybrid node.
                     # updated after node is part of network with updateGammaz!
    edge::Array{Edge,1}
    hasHybEdge::Bool #is there a hybrid edge in edge? only needed when hybrid=false (tree node)
    isBadDiamondI::Bool # for hybrid node, is it bad diamond case I, update in updateGammaz!
    isBadDiamondII::Bool # for hybrid node, is it bad diamond case II, update in updateGammaz!
    isExtBadTriangle::Bool # for hybrid node, is it extremely bad triangle, udpate in updateGammaz!
    isVeryBadTriangle::Bool # for hybrid node, is it very bad triangle, udpate in updateGammaz!
    isBadTriangle::Bool # for hybrid node, is it very bad triangle, udpate in updateGammaz!
    inCycle::Int # = hybrid node if this node is part of a cycle created by such hybrid node, -1 if not part of cycle
    prev::Union{Nothing,ANode} # previous node in cycle, used in updateInCycle. set to "nothing" to begin with
    k::Int # num nodes in cycle, only stored in hybrid node, updated after node becomes part of network
           # default -1
    typeHyb::Int8 # type of hybridization (1,2,3,4, or 5), needed for quartet network only. default -1
    name::AbstractString
    # inner constructor: set hasHybEdge depending on edge
    Node() = new(-1,false,false,-1.,Edge[],false,false,false,false,false,false,-1,nothing,-1,-1,"")
    Node(number::Int, leaf::Bool) = new(number,leaf,false,-1.,[],false,false,false,false,false,false,-1.,nothing,-1,-1,"")
    Node(number::Int, leaf::Bool, hybrid::Bool) = new(number,leaf,hybrid,-1.,[],hybrid,false,false,false,false,false,-1.,nothing,-1,-1,"")
    Node(number::Int, leaf::Bool, hybrid::Bool, edge::Array{Edge,1})=new(number,leaf,hybrid,-1.,edge,!all((e->!e.hybrid),edge),false,false,false,false,false,-1.,nothing,-1,-1,"")
    Node(number::Int, leaf::Bool, hybrid::Bool,gammaz::Float64, edge::Array{Edge,1}) = new(number,leaf,hybrid,gammaz,edge,!all((e->!e.hybrid), edge),false,false,false,false,false,-1.,nothing,-1,-1,"")
end

# partition type
mutable struct Partition
    cycle::Vector{Int} #hybrid node number for cycle (or cycles)
    edges::Vector{Edge} #edges in partition
end

abstract type Network end

# warning: no attempt to make sure the direction of edges matches with the root
# warning: no check if it is network or tree, node array can have no hybrids
# warning: nodes and edges need to be defined and linked before adding to a network
"""
    HybridNetwork

Subtype of abstract `Network` type.
Explicit network or tree with the following attributes:

- numTaxa (taxa are tips, i.e. nodes attached to a single edge)
- numNodes (total number of nodes: tips and internal nodes)
- numEdges
- numHybrids (number of hybrid nodes)
- edge (array of Edges)
- node (array of Nodes)
- root (index of root in vector 'node'. May be artificial, for printing and traversal purposes only.)
- hybrid (array of Nodes: those are are hybrid nodes)
- leaf (array of Nodes: those that are leaves)
- loglik (score after fitting network to data, i.e. negative log pseudolik for SNaQ)
- isRooted (true or false)
"""
mutable struct HybridNetwork <: Network
    numTaxa::Int  # cannot set in constructor for congruence
    numNodes::Int
    numEdges::Int
    node::Array{Node,1}
    edge::Array{Edge,1}
    root::Int # node[root] is the root node, default 1
    names::Array{String,1} # translate table for taxon names --but also includes hybrid names...
    hybrid::Array{Node,1} # array of hybrid nodes in network
    numHybrids::Int # number of hybrid nodes
    cladewiseorder_nodeIndex::Vector{Int} # index in 'node' for "cladewise" preorder in main tree
    visited::Array{Bool,1} # reusable array of booleans
    edges_changed::Array{Edge,1} # reusable array of edges
    nodes_changed::Array{Node,1} # reusable array of nodes. used for preorder traversal
    leaf::Array{Node,1} # array of leaves
    ht::Vector{Float64} # vector of parameters to optimize
    numht::Vector{Int} # vector of number of the hybrid nodes and edges in ht e.g. [3,6,8,...], 2 hybrid nodes 3,6, and edge 8 is the 1st identifiable
    numBad::Int # number of bad diamond I hybrid nodes, set as 0
    hasVeryBadTriangle::Bool # true if the network has extremely/very bad triangles that should be ignored
    index::Vector{Int} #index in net.edge, net.node of elements in net.ht to make updating easy
    loglik::Float64 # value of the min -loglik after optBL
    blacklist::Vector{Int} # reusable array of integers, used in afterOptBL
    partition::Vector{Partition} # to choose edges from a partition only to avoid intersecting cycles
    cleaned::Bool # attribute to know if the network has been cleaned after readm default false
    isRooted::Bool # to know if network is rooted, e.g. after directEdges! (which updates isChild1 of each edge)
    # inner constructor
    function HybridNetwork(node::Array{Node,1},edge::Array{Edge,1})
        hybrid=Node[];
        leaf=Node[];
        for n in node
            if n.hybrid push!(hybrid,n); end
            if n.leaf   push!(leaf,  n); end
        end
        new(size(leaf,1),size(node,1),size(edge,1),node,edge,1,[],hybrid,size(hybrid,1), #numTaxa,...,numHybrids
            [],[],[],[],leaf,[],[], #cladewiseorder,...,numht
            0,false,[],0,[],[],false,false) #numBad...
    end
    HybridNetwork() = new(0,0,0,[],[],0,[],[],0, # numTaxa ... numHybrid
                          [],[],[],[],[],[],[], # cladewiseorder...
                          0,false,[],0,[],[],false,false); # numBad ...
end

# type created from a HybridNetwork only to extract a given quartet
"""
    QuartetNetwork(net::HybridNetwork)

Subtype of `Network` abstract type.
A `QuartetNetwork` object is an internal type used to calculate the
expected CFs of quartets on a given network.
Attributes of the `QuartetNetwork` objects need not be updated at a given time (see below).

The procedure to calculate expected CFs for a given network is as follows:
1. A `QuartetNetwork` object is created for each `Quartet` using
   `extractQuartet!(net,d)` for `net::HybridNetwork` and `d::DataCF`
2. The vector `d.quartet` has all the `Quartet` objects, each with a `QuartetNetwork`
   object (`q.qnet`). Attibutes in `QuartetNetwork` are not updated at this point
3. Attributes in `QuartetNetwork` are partially updated when calculating the
   expected CF (`calculateExpCFAll!`). To calculate the expected CF for this quartet,
   we need to update the attributes: `which`, `typeHyb`, `t1`, `split`, `formula`, `expCF`.
   To do this, we need to modify the `QuartetNetwork` object (i.e. merge edges,...).
   But we do not want to modify it directly because it is connected to the original
   `net` via a map of the edges and nodes, so we use a deep copy:
   `qnet=deepcopy(q.qnet)` and then `calculateExpCFAll!(qnet)`.
   Attributes that are updated on the original `QuartetNetwork` object `q.qnet` are:
    - `q.qnet.hasEdge`: array of booleans of length equal to `net.edge` that shows which identifiable edges and gammas of `net` (`net.ht`) are in `qnet` (and still identifiable). Note that the first elements of the vector correspond to the gammas.
    - `q.qnet.index`: length should match the number of trues in `qnet.hasEdge`. It has the indexes in `qnet.edge` from the edges in `qnet.hasEdge`. Note that the first elements of the vector correspond to the gammas.
    - `q.qnet.edge`: list of edges in `QuartetNetwork`. Note that external edges in `net` are collapsed when they appear in `QuartetNetwork`, so only internal edges map directly to edges in `net`
    - `q.qnet.expCF`: expected CF for this `Quartet`

Why not modify the original `QuartetNetwork`? We wanted to keep the original
`QuartetNetwork` stored in `DataCF` with all the identifiable edges, to be able
to determine if this object had been changed or not after a certain optimization.

The process is:

1. Deep copy of full network to create `q.qnet` for `Quartet q`.
   This `QuartetNetwork` object has only 4 leaves now, but does not have merged edges
   (the identifiable ones) so that we can correspond to the edges in net.
   This `QuartetNetwork` does not have other attributes updated.
2. For the current set of branch lengths and gammas, we can update the attributes
   in `q.qnet` to compute the expected CF. The functions that do this will "destroy"
   the `QuartetNetwork` object by merging edges, removing nodes, etc... So, we do
   this process in `qnet=deepcopy(q.qnet)`, and at the end, only update `q.qnet.expCF`.
3. After we optimize branch lengths in the full network, we want to update the
   branch lengths in `q.qnet`. The edges need to be there (which is why we do
   not want to modify this `QuartetNetwork` object by merging edges), and
   we do not do a deep-copy of the full network again. We only change the values
   of branch lengths and gammas in `q.qnet`, and we can re-calculate the expCF
   by creating a deep copy `qnet=deepcopy(q.qnet)` and run the other functions
   (which merge edges, etc) to get the `expCF`.

Future work: there are definitely more efficient ways to do this (without the deep copies).
In addition, currently edges that are no longer identifiable in `QuartetNetwork`
do not appear in `hasEdge` nor `index`. Need to study this.

```jldoctest
julia> net0 = readTopology("(s17:13.76,(((s3:10.98,(s4:8.99,s5:8.99)I1:1.99)I2:0.47,(((s6:2.31,s7:2.31)I3:4.02,(s8:4.97,#H24:0.0::0.279)I4:1.36)I5:3.64,((s9:8.29,((s10:2.37,s11:2.37)I6:3.02,(s12:2.67,s13:2.67)I7:2.72)I8:2.89)I9:0.21,((s14:2.83,(s15:1.06,s16:1.06)I10:1.78)I11:2.14)#H24:3.52::0.72)I12:1.47)I13:1.48)I14:1.26,(((s18:5.46,s19:5.46)I15:0.59,(s20:4.72,(s21:2.40,s22:2.40)I16:2.32)I17:1.32)I18:2.68,(s23:8.56,(s1:4.64,s2:4.64)I19:3.92)I20:0.16)I21:3.98)I22:1.05);");

julia> net = readTopologyLevel1(writeTopology(net0)) ## need level1 attributes for functions below
HybridNetwork, Un-rooted Network
46 edges
46 nodes: 23 tips, 1 hybrid nodes, 22 internal tree nodes.
tip labels: s17, s3, s4, s5, ...
(s4:8.99,s5:8.99,(s3:10.0,((((s6:2.31,s7:2.31)I3:4.02,(s8:4.97,#H24:0.0::0.279)I4:1.36)I5:3.64,((s9:8.29,((s10:2.37,s11:2.37)I6:3.02,(s12:2.67,s13:2.67)I7:2.72)I8:2.89)I9:0.21,((s14:2.83,(s15:1.06,s16:1.06)I10:1.78)I11:2.14)#H24:3.52::0.721)I12:1.47)I13:1.48,((((s18:5.46,s19:5.46)I15:0.59,(s20:4.72,(s21:2.4,s22:2.4)I16:2.32)I17:1.32)I18:2.68,(s23:8.56,(s1:4.64,s2:4.64)I19:3.92)I20:0.16)I21:3.98,s17:10.0)I22:1.26)I14:0.47)I2:1.99)I1;

julia> q1 = Quartet(1,["s1", "s16", "s18", "s23"],[0.296,0.306,0.398])
number: 1
taxon names: ["s1", "s16", "s18", "s23"]
observed CF: [0.296, 0.306, 0.398]
pseudo-deviance under last used network: 0.0 (meaningless before estimation)
expected CF under last used network: Float64[] (meaningless before estimation)

julia> qnet = PhyloNetworks.extractQuartet!(net,q1)
taxa: ["s1", "s16", "s18", "s23"]
number of hybrid nodes: 1

julia> sum([e.istIdentifiable for e in net.edge]) ## 23 identifiable edges in net
23

julia> idedges = [ee.number for ee in net.edge[[e.istIdentifiable for e in net.edge]]];

julia> print(idedges)
[5, 6, 9, 11, 12, 13, 17, 20, 21, 22, 26, 27, 28, 29, 30, 31, 34, 38, 39, 40, 44, 45, 46]

julia> length(qnet.hasEdge) ## 24 = 1 gamma + 23 identifiable edges
24

julia> sum(qnet.hasEdge) ## 8 = 1 gamma + 7 identifiable edges in qnet
8

julia> print(idedges[qnet.hasEdge[2:end]]) ## 7 id. edges: [12, 13, 29, 30, 31, 45, 46]
[12, 13, 29, 30, 31, 45, 46]

julia> qnet.edge[qnet.index[1]].number ## 11 = minor hybrid edge
11
```
"""
mutable struct QuartetNetwork <: Network
    numTaxa::Int
    numNodes::Int
    numEdges::Int
    node::Array{Node,1}
    edge::Array{Edge,1}
    hybrid::Array{Node,1} # array of hybrid nodes in network
    leaf::Array{Node,1} # array of leaves
    numHybrids::Int # number of hybrid nodes
    hasEdge::Array{Bool,1} # array of boolean with all the original identifiable edges of HybridNetwork and gammas (net.ht)
    quartetTaxon::Array{String,1} # the quartet taxa in the order it represents. Points to same array as its Quartet.taxon
    which::Int8 # 0 it tree quartet, 1 is equivalent to tree quartet and 2 if two minor CF different, default -1
    typeHyb::Array{Int8,1} #array with the type of hybridization of each hybrid node in the quartet
    t1::Float64 # length of internal edge, used when qnet.which=1, default = -1
    names::Array{String,1} # taxon and node names, same order as in network.node
    split::Array{Int8,1} # split that denotes to which side each leaf is from the split, i.e. [1,2,2,1] means that leaf1 and 4 are on the same side of the split, default -1,-1,-1,-1
    formula::Array{Int8,1} # array for qnet.which=1 that indicates if the expCf is major (1) or minor (2) at qnet.expCF[i] depending on qnet.formula[i], default -1,-1,-1
    expCF::Array{Float64,1} # three expected CF in order 12|34, 13|24, 14|23 (matching obsCF from qnet.quartet), default [0,0,0]
    indexht::Vector{Int} # index in net.ht for each edge in qnet.ht
    changed::Bool # true if the expCF would be changed with the current parameters in the optimization, to recalculate, default true
    index::Vector{Int} # index in qnet.edge (qnet.node for gammaz) of the members in qnet.indexht to know how to find quickly in qnet
    # inner constructor
    function QuartetNetwork(net::HybridNetwork)
        net2 = deepcopy(net); #fixit: maybe we dont need deepcopy of all, maybe only arrays
        new(net2.numTaxa,net2.numNodes,net2.numEdges,net2.node,net2.edge,net2.hybrid,net2.leaf,net2.numHybrids, [true for e in net2.edge],[],-1,[], -1.,net2.names,Int8[-1,-1,-1,-1],Int8[-1,-1,-1],[0,0,0],[],true,[])
        #new(sum([n.leaf?1:0 for n in net.node]),size(net.node,1),size(net.edge,1),copy(net.node),copy(net.edge),copy(net.hybrid),size(net.hybrid,1), [true for e in net2.edge],[],-1,[],-1.,net2.names,[-1,-1,-1,-1],[-1,-1,-1],[],true,[])
    end
    QuartetNetwork() = new(0,0,0,[],[],[],[],0,[],[],-1,[],-1.0,[],[],[],[],[],true,[])
end

abstract type AQuartet end

"""
    Quartet

type that saves the information on a given 4-taxon subset. It contains the following attributes:

- number: integer
- taxon: vector of taxon names, like t1 t2 t3 t4
- obsCF: vector of observed CF, in order 12|34, 13|24, 14|23
- logPseudoLik
- ngenes: number of gene trees used to compute the observed CF; -1.0 if unknown
- qnet: [`QuartetNetwork`](@ref), which saves the expCF after snaq estimation to
  emphasize that the expCF depend on a specific network, not the data

see also: [`QuartetT`](@ref) for quartet with data of user-defined type `T`,
using a mapping between quartet indices and quartet taxa.
"""
mutable struct Quartet <: AQuartet
    number::Int
    taxon::Array{String,1} # taxa 1234. qnet.quartetTaxon points to the same array.
    obsCF::Array{Float64,1} # three observed CF in order 12|34, 13|24, 14|23
    qnet::QuartetNetwork # quartet network for the current network (want to keep as if private attribute)
    logPseudoLik::Float64 # log pseudolik value for the quartet. 0.0 by default
    ngenes::Float64 # number of gene trees used to compute the obsCV, default -1.; Float in case ngenes is average
    # inner constructor: to guarantee obsCF are only three and add up to 1
    function Quartet(number::Integer,t1::AbstractString,t2::AbstractString,t3::AbstractString,t4::AbstractString,obsCF::Array{Float64,1})
        size(obsCF,1) != 3 ? error("observed CF vector should have size 3, not $(size(obsCF,1))") : nothing
        0.99 < sum(obsCF) < 1.02 || @warn "observed CF should add up to 1, not $(sum(obsCF))"
        new(number,[t1,t2,t3,t4],obsCF,QuartetNetwork(),0.0,-1.0);
    end
    function Quartet(number::Integer,t1::Array{String,1},obsCF::Array{Float64,1})
        size(obsCF,1) != 3 ? error("observed CF vector should have size 3, not $(size(obsCF,1))") : nothing
        0.99< sum(obsCF) < 1.02 || @warn "observed CF should add up to 1, not $(sum(obsCF))"
        size(t1,1) != 4 ? error("array of taxa should have size 4, not $(size(t1,1))") : nothing
        0.0 <= obsCF[1] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[1]) for $(t1)")
        0.0 <= obsCF[2] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[2]) for $(t1)")
        0.0 <= obsCF[3] <= 1.0 || error("obsCF must be between (0,1), but it is $(obsCF[3]) for $(t1)")
        new(number,t1,obsCF,QuartetNetwork(),0.0,-1.0);
    end
    Quartet() = new(0,[],[],QuartetNetwork(),0.0,-1.0)
end

"""
QuartetT{T}

Generic type for 4-taxon sets. Fields:
- `number`: rank of the 4-taxon set
- `taxonnumber`: static vector of 4 integers, assumed to be distinct and sorted
- `data`: object of type `T`

For easier look-up, a unique mapping is used between the rank (`number`) of a
4-taxon set and its 4 taxa (see [`quartetrank`](@ref) and [`nchoose1234`](@ref)):

rank-1 = (t1-1) choose 1 + (t2-1) choose 2 + (t3-1) choose 3 + (t4-1) choose 4

# examples

```jldoctest
julia> nCk = PhyloNetworks.nchoose1234(5)
6×4 Matrix{Int64}:
 0   0   0  0
 1   0   0  0
 2   1   0  0
 3   3   1  0
 4   6   4  1
 5  10  10  5

julia> PhyloNetworks.QuartetT(1,3,4,6, [.92,.04,.04, 100], nCk)
4-taxon set number 8; taxon numbers: 1,3,4,6
data: [0.92, 0.04, 0.04, 100.0]
```
"""
struct QuartetT{T} <: AQuartet where T
    number::Int
    taxonnumber::StaticArrays.SVector{4,Int}
    data::T
end
function Base.show(io::IO, obj::QuartetT{T}) where T
    disp = "4-taxon set number $(obj.number); taxon numbers: "
    disp *= join(obj.taxonnumber,",")
    disp *= "\ndata: "
    print(io, disp)
    print(io, obj.data)
end
function QuartetT(tn1::Int,tn2::Int,tn3::Int,tn4::Int, data::T, nCk::Matrix, checksorted=true::Bool) where T
    if checksorted
        (tn1<tn2 && tn2<tn3 && tn3<tn4) || error("taxon numbers must be sorted")
    end
    QuartetT{T}(quartetrank(tn1,tn2,tn3,tn4,nCk), SVector(tn1,tn2,tn3,tn4), data)
end

"""
    quartetrank(t1,t2,t3,t4, nCk::Matrix)
    quartetrank([t1,t2,t3,t4], nCk)

Return the rank of a four-taxon set with taxon numbers `t1,t2,t3,t4`,
assuming that `ti`s are positive integers such that t1<t2, t2<t3 and t3<t4
(assumptions not checked!).
`nCk` should be a matrix of "n choose k" binomial coefficients:
see [`nchoose1234`](@ref).

# examples

```jldoctest
julia> nCk = PhyloNetworks.nchoose1234(5)
6×4 Matrix{Int64}:
 0   0   0  0
 1   0   0  0
 2   1   0  0
 3   3   1  0
 4   6   4  1
 5  10  10  5

julia> PhyloNetworks.quartetrank([1,2,3,4], nCk)
1

julia> PhyloNetworks.quartetrank([3,4,5,6], nCk)
15
```
"""
@inline function quartetrank(tnum::AbstractVector, nCk::Matrix)
    quartetrank(tnum..., nCk)
end
@inline function quartetrank(t1::Int, t2::Int, t3::Int, t4::Int, nCk::Matrix)
    # rank-1 = t1-1 choose 1 + t2-1 choose 2 + t3-1 choose 3 + t4-1 choose 4
    return nCk[t1,1] + nCk[t2,2] + nCk[t3,3] + nCk[t4,4] + 1
end

"""
    nchoose1234(nmax)

`nmax+1 x 4` matrix containing the binomial coefficient
"n choose k" in row `n+1` and column `k`. In other words,
`M[i,k]` gives "i-1 choose k". It is useful to store these
values and look them up to rank (a large number of) 4-taxon sets:
see [`quartetrank`](@ref).
"""
function nchoose1234(nmax::Int)
    # compute nC1, nC2, nC3, nC4 for n in [0, nmax]: used for ranking quartets
    M = Matrix{Int}(undef, nmax+1, 4)
    for i in 1:(nmax+1)
        M[i,1] = i-1 # n choose 1 = n. row i is for n=i-1
    end
    M[1,2:4] .= 0 # 0 choose 2,3,4 = 0
    for i in 2:(nmax+1)
        for k in 2:4 # to choose k items in 1..n: the largest could be n, else <= n-1
            M[i,k] = M[i-1,k-1] + M[i-1,k]
        end
    end
    return M
end

# Data on quartet concordance factors -------

"""
    DataCF

type that contains the following attributes:

- quartet (vector of Quartets)
- numQuartets
- tree (vector of trees: empty if a table of CF was input instead of list of trees)
- numTrees (-1 if a table CF was input instead of list of trees)
- repSpecies (taxon names that were repeated in table of CF or input gene trees: used inside snaq for multiple alleles case)

The list of Quartet may be accessed with the attribute .quartet.
If the input was a list of trees, the HybridNetwork's can be accessed with the attribute .tree.
For example, if the DataCF object is named d, d.quartet[1] will show the first quartet
and d.tree[1] will print the first input tree.
"""
mutable struct DataCF # fixit
    quartet::Array{Quartet,1} # array of quartets read from CF output table or list of quartets in file
    numQuartets::Integer # number of quartets
    tree::Vector{HybridNetwork} #array of input gene trees
    numTrees::Integer # number of gene trees
    repSpecies::Vector{String} #repeated species in the case of multiple alleles
    DataCF(quartet::Array{Quartet,1}) = new(quartet,length(quartet),[],-1,[])
    DataCF(quartet::Array{Quartet,1},trees::Vector{HybridNetwork}) = new(quartet,length(quartet),trees,length(trees),[])
    DataCF() = new([],0,[],-1,[])
end

# aux type for the updateBL function
mutable struct EdgeParts
    edgenum::Int
    part1::Vector{Node}
    part2::Vector{Node}
    part3::Vector{Node}
    part4::Vector{Node}
end

struct RootMismatch <: Exception
    msg::String
end
RootMismatch() = RootMismatch("")
Base.showerror(io::IO, e::RootMismatch) = print(io, "RootMismatch: ", e.msg);


text:
@info "Loading CaesarZMQ RoMEPlotting specific utilities."

# FIXME should be a @requires
# global isLoadingPlottingEndpoints = false
# try
#     getfield(Main, :RoMEPlotting)
#     # If that works, then we can enable plotting
#     global isLoadingPlottingEndpoints = true
# catch ex
# end

# if isLoadingPlottingEndpoints
@info "[ZMQ Server] Plotting endpoints are enabled!"

function plotKDE(configDict, fg, requestDict)::Dict{String, Any}
    error("Not implemented yet!")
end

function plotPose(configDict, fg, requestDict)::Dict{String, Any}
    error("Not implemented yet!")
end

function drawPoses(configDict, fg, requestDict)::Dict{String, Any}
    plotRequest = Unmarshal.unmarshal(PlotRequest, requestDict["payload"])
    pl = RoMEPlotting.drawPoses(fg);
    io = IOBuffer();
    Gadfly.draw(PNG(io, (plotRequest.widthPx)px, (plotRequest.heightPx)px), pl);
    return Dict{String, Any}("status" => "ok", "mimetype" => "image/png", "encoding" => "none", "data" => io.data)
end

function drawPosesLandms(configDict, fg, requestDict)::Dict{String, Any}
    plotRequest = Unmarshal.unmarshal(PlotRequest, requestDict["payload"])
    pl = RoMEPlotting.drawPosesLandms(fg);
    io = IOBuffer();
    Gadfly.draw(PNG(io, (plotRequest.widthPx)px, (plotRequest.heightPx)px), pl);
    return Dict{String, Any}("status" => "ok", "mimetype" => "image/png", "encoding" => "none", "data" => io.data)
end

    # else
#     @info "[ZMQ Server] Plotting endpoints are disabled in ZMQ server. Please call 'using RoMEPlotting' before 'using Caesar' if you would like to enable them."
# end


text:
using LULESH
using Test
using MPI

function run_example(name, args; nranks=nothing)
    example = joinpath(@__DIR__, "..", "examples", name)

    @testset "$(basename(example)) $args $nranks" begin
        cmd = `$(Base.julia_cmd()) --startup-file=no $(example) $(args)`
        if nranks !== nothing
            mpiexecjl = joinpath(dirname(pathof(MPI)), "..", "bin", "mpiexecjl")
            cmd = `$(mpiexecjl) -n $(nranks) $(cmd)`
        end
        @debug "Testing $example" cmd
        @test success(pipeline(cmd, stderr=stderr))
    end
end

@testset "examples" begin
    @testset "benchmark.jl" begin
        run_example("benchmark.jl", `-s 45`)
        run_example("benchmark.jl", `-s 45`; nranks=1)
        if Sys.CPU_THREADS >= 8
            run_example("benchmark.jl", `-s 45`; nranks=8)
        end
    end
end


text:
@with_kw struct OptimizerParams
    iterations_limit::Int64 = 1e3
    time_limit::Float64 = 30.0
    show_trace::Bool = true
    trace_length::Int64 = 20
    cvcheckers::Set{ConvergenceChecker} = Set{ConvergenceChecker}([FirstOrderOptimality()])
end


text:
@testset "StepRange" begin
    for R in (StepMRange, StepSRange)
        @testset "$R" begin
            r = R(1, 1, 10)
            b = StepRange(1, 1, 10)
            @test R(1:10) == 1:1:10
            @test eltype(R(UInt(1), UInt(1), UInt(10))) == UInt
            @test R{Int,Int}(r) === r
            @test eltype(R{UInt,UInt}(r)) == UInt
            @test first(r) == r.start
            @test last(r) == r.stop

            @test intersect(r, r[2]) == intersect(b, b[2])
            @test intersect(r, r[2]) == intersect(b, b[2])

            @test_throws ErrorException r.notfield
            if R == StepMRange
                setproperty!(r, :step, 2)
                @test r == StepMRange(1, 2, 9)

                setproperty!(r, :start, 2)
                @test r == StepMRange(2, 2, 8)

                setproperty!(r, :stop, 10)
                @test r == StepMRange(2, 2, 10)

                @test_throws ErrorException setproperty!(r, :anything, 3)
            end
        end
    end
end


text:
using OlfactoryTopographicBehavior
using DataFrames
using Unitful: Hz, s
using IntervalSets
using PyPlot: PyPlot, plt
using Test

@testset "OlfactoryTopographicBehavior.jl" begin
    # TrialData extraction
    df = DataFrame("soundlick"     => [0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0],
                   "lickometer"    => [1,-1, 1,-1, 1, 0, 0, 0, 1,-1, 1, 0, 1,-1, 1],
                   "odortiming"    => [0, 0, 0, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 0, 0],
                   "odordirection" => [1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1])
    fs = 2Hz
    tds = trialdatas(df; lohi_sound = (0.25, 0.75), fs, soundduration=2s, fcarrierlick=1Hz, minlickduration=1s)
    td = only(tds)
    @test td.trialrange == 2:length(df.soundlick)
    @test td.soundinterval == 1s .. 3s
    @test only(td.lickonsets) == 2s
    @test td.odordirection == 2
    @test td.odorinterval == 3.5s .. 4.5s
    tds = trialdatas(df; lohi_sound = (0.25, 0.75), fs, soundduration=2s, fcarrierlick=1Hz, minlickduration=0s)
    td = only(tds)
    @test td.lickonsets == [2s, 5s]
    @test_throws ErrorException("peak frequency 1.0 Hz does not match expected 60 Hz") trialdatas(df; lohi_sound = (0.25, 0.75), fs, soundduration=3s)
    @test_throws ErrorException("soundlick pulse of duration 2.0 s starting at 4 not recognized") trialdatas(df; lohi_sound = (0.25, 0.75), fs, soundduration=3s, fcarrierlick=1Hz, minlickduration=1s)

    # Display
    io = IOBuffer()
    print(io, td)
    str = String(take!(io))
    @test startswith(str, "Trial over indices 2:15")

    # Construction from a matrix
    mat = [df.soundlick df.lickometer df.odortiming df.odordirection]
    dfmat = session_frame(mat; fs, soundlick=1, lickometer=2, odortiming=3, odordirection=4)
    @test dfmat.df == df
    print(io, dfmat)
    str = String(take!(io))
    @test endswith(str, "Recorded at 2.0f0 Hz")
    @test dfmat.lickometer == df.lickometer
    @test trialdatas(dfmat;  lohi_sound = (0.25, 0.75), soundduration=2s, fcarrierlick=1Hz, minlickduration=1s) ==
          trialdatas(df; fs, lohi_sound = (0.25, 0.75), soundduration=2s, fcarrierlick=1Hz, minlickduration=1s)

    # lick detection
    td = TrialData(1:0, 1s .. 2s, 1s, 1.5s .. 3.0s, -1, [0.8s])
    @test !is_lick(td)
    push!(td.lickonsets, 6s)
    @test !is_lick(td)
    @test  is_lick(td; deadline=10s)
    empty!(td.lickonsets)
    @test !is_lick(td)
    push!(td.lickonsets, 1.2s)
    @test  is_lick(td)
    empty!(td.lickonsets)
    push!(td.lickonsets, 0.8s, 1.2s)
    @test  is_lick(td)

    # trial-type aggregation
    tds = [TrialData(1:0, 1s .. 2s, 2.5s, 2s .. 3s, 1, [0.8s]),
           TrialData(1:0, 1s .. 2s, 2.5s, 2s .. 3s, 1, [2.7s]),
           TrialData(1:0, 1s .. 2s, 2.5s, 2s .. 3s, 2, [0.8s]),
           TrialData(1:0, 1s .. 2s, 2.5s, 2s .. 3s, 2, [2.7s])
    ]
    agg = collect_by_trialtype(td->isodd(td.odordirection), [1,2,3,4], tds)
    @test agg["TP"] == [2]
    @test agg["TN"] == [3]
    @test agg["FP"] == [4]
    @test agg["FN"] == [1]

    # Trial-interval selection
    licks = repeat(Bool[0, 1], 20)
    @test trialrange(licks, 0.5) == firstindex(licks):lastindex(licks)
    @test trialrange(licks, 0.85) != firstindex(licks):lastindex(licks)
    @test trialrange(licks, 0.15) != firstindex(licks):lastindex(licks)
    n = lastindex(licks)
    push!(licks, false)
    @test trialrange(licks, 0.5) == firstindex(licks):lastindex(licks)
    append!(licks, fill(false, 40))
    idx = trialrange(licks, 0.5)
    @test first(idx) == firstindex(licks)
    @test n <= last(idx)+1 < lastindex(licks)
    # A fairly naturalistic pattern: an overabundance of licks early in the session,
    # followed by near-ideal behavior, followed by underlicking due to satiation.
    n1, n2, n3 = 50, 100, 50
    licks = [repeat(Bool[0, 1, 1, 1], n1); repeat(Bool[0, 1], n2); repeat(Bool[0, 0, 0, 1], n3)]
    idx = trialrange(licks, 0.5)
    @test 4*n1 .+ (1:2*n2-1) ⊆ idx
    @test first(idx) > firstindex(licks)
    @test last(idx) < lastindex(licks)
    n1, n2, n3 = 100, 100, 100
    for pval in (0.05, 0.4)
        nfail = 0
        for _ = 1:100
            empty!(licks)
            for _ = 1:n1
                push!(licks, rand() < 0.85)
            end
            for _ = 1:n2
                push!(licks, rand() < 0.5)
            end
            for _ = 1:n3
                push!(licks, rand() < 0.15)
            end
            idx = trialrange(licks, 0.5, pval)
            nfail += !(n1+ceil(Int, 0.1*n2):n1+floor(Int, 0.9*n2) ⊆ idx)
            nfail += first(idx) <= firstindex(licks)
            nfail += last(idx)  >= lastindex(licks)
            # @test n1+1:n1+n2 ⊆ idx
            # @test first(idx) > firstindex(licks)
            # @test last(idx)  < lastindex(licks)
        end
        @test 30*pval< nfail < 300*pval
    end
    licks = repeat(Bool[0, 1, 1, 1], 10)
    @test trialrange(licks, 0.5) != firstindex(licks):lastindex(licks)
    @test trialrange(licks, 0.75) == firstindex(licks):lastindex(licks)

    # Plotting
    dfplt = copy(df)
    dfplt.sniff = randn(size(dfplt, 1))
    dfrplt = OlfactoryTopographicBehavior.DataFrameRate(dfplt, fs)
    tds = trialdatas(dfrplt; lohi_sound = (0.25, 0.75), soundduration=2s, fcarrierlick=1Hz, minlickduration=1s)
    td = only(tds)
    fig = plt.figure()
    axsniff, axtiming = plottrial(dfrplt, td, fig)
    show(io, axsniff)
    str = String(take!(io))
    @test occursin("AxesSubplot:xlabel='Time (s)', ylabel='Chest expansion", str)
    plt.close(fig)

    # Sniff analysis
    fill!(dfrplt.sniff, 0)
    dfrplt.sniff[end÷2+1:end] .= 1
    ss = OlfactoryTopographicBehavior.sniffslope(dfrplt, 0.5s)
    @test all(<(0.01), ss[1:3])
    @test all(<(0.01), ss[end-2:end])
    val, idx = findmax(ss)
    @test val > 0.3
    @test idx ∈ (7, 8)
    @test ss[7] ≈ ss[8]
    time, freq, sgs = OlfactoryTopographicBehavior.sniffgrams(dfrplt, tds)
    @test only(time) == 1s
    @test freq == [0Hz, 0.5Hz, 1Hz]
    @test only(sgs) == reshape([2, 0, 0], (3, 1))
    @test OlfactoryTopographicBehavior.sniffgrams_power(freq, sgs; finterval=0.5Hz .. 1Hz) == [0]
end


text:
using Statistics
export p_dim, hvIndicator,p_MADS, BiMADS_status, paretoCoverage,plot_Bpoint
#TODO BiMADS
logocolors = Colors.JULIA_LOGO_COLORS
"""
points for BiMADS Algo in objective space
"""
mutable struct B_points
    cost::Vector{Float64}  #f1 and f2 value in objective space
    x_now::Vector{Float64}
    weight::Int            #weight for each undominated point to determine the δ
    function B_points(cost::Vector{Float64},x::Vector{Float64}, w::Int=0)
        p = new()
        p.cost = cost
        p.x_now=x
        p.weight = w
        return p
    end
end
"""
Status for BiMADS main iteration, also contain stop conditions

"""
mutable struct BiMADS_status
    iteration::Int64                #Total iteration spent
    func_evaluation::Int64          #Total function evaluation spent
    total_time::Float64             #Totlal running time
    hypervolume::Float64            #Hypervolume for current Pareto set
    opt_status::OptimizationStatus
    opt_string::String
    start_time::Float64
    function BiMADS_status()
        s=new()
        s.iteration=0
        s.func_evaluation=0
        s.total_time=0.0
        s.hypervolume=0.0
        s.opt_status=Unoptimized
        s.opt_string="Unoptimized"
        s.start_time=time()
        return s
    end
end

"""
    plot_Bpoint(points::Vector{B_points})

Plot the Pareto front
"""
function plot_Bpoint(points::Vector{B_points})
    fig = Plots.scatter()
    for i = 1:length(points)
        fig = plot!([points[i].cost[1]],[points[i].cost[2]],seriestype = :scatter,aspect_ratio=1,legend = false,color=logocolors.red,show = true)
    end
    plot!(fig,xlabel="f1 cost",ylabel="f2 cost")
return fig
end

function get_adj_dis(points::Vector{B_points})
    println("-----------------")
    for i=1:length(points)-1
        display(LinearAlgebra.norm(points[i+1].cost - points[i].cost)^2)
    end
    println("-----------------")
end

function paretoCoverage(paretoSet::Vector{B_points})::Tuple{Float64,Float64}
    points=Vector{Float64}()
    for i=1:length(paretoSet)-1
        push!(points,LinearAlgebra.norm(paretoSet[i+1].cost - paretoSet[i].cost)^2)
    end
    return mean(points), std(points)
end

"""
    hvIndicator(paretoSet::Vector{B_points},factor=1.1)::Float64

Hyper-Volume indicator for evalating the Pareto Front
factor is used for choosing the referencepoint for hypervolume calculation
"""
function hvIndicator(paretoSet::Vector{B_points},factor=1.1)::Float64
    points=Vector{Vector{Float64}}()
    for i=1:length(paretoSet)-1
        push!(points,paretoSet[i].cost)
    end
    length(points)<2 && return 0.
     ref=factor.*[last(points)[1],first(points)[2]]
     # normalize_factor=(ref[1]-first(points)[1]).*(ref[2]-last(points)[2])./2
     normalize_factor=(last(points)[1]-first(points)[1]).*(first(points)[2]-last(points)[2])./2
     hv_volume=0.
     area(x::Vector{Float64},y::Vector{Float64})=(abs(x[1]-y[1])).*(abs(x[2]-y[2]))
     for p in points
         hv_volume+=area(ref,p)
         ref[2]=p[2]
     end
     return hv_volume/normalize_factor
end

"""
    p_split(p::DSProblem)::Tuple{DSProblem,DSProblem}

To split Bi-obj problem into two single obj problems
"""
function p_split(p::DSProblem)::Tuple{DSProblem,DSProblem}
    x = zeros(Float64, p.N)
    f1, f2 = p.objective(x)
    p1::DSProblem = deepcopy(p::DSProblem)
    p2::DSProblem = deepcopy(p::DSProblem)
    SetObjective(p1, f1)
    SetObjective(p2, f2)
    return p1, p2
end

"""
    get_Bpoints(p1::DSProblem,p2::DSProblem,flag::Int)

Get the points to determine the reference point and Pareto front
flag for identifying optimize p1 or p2
"""
function get_Bpoints(p1::DSProblem, p2::DSProblem, flag::Int)
    temp_Bpoints = Vector{B_points}()
    Setup(p1)
    # push B_point for initial point
    flag == 1 &&
        push!(temp_Bpoints, B_points([p1.objective(p1.x), p2.objective(p1.x)],p1.x))

    while _check_stoppingconditions(p1)
        p1.full_output && OutputIterationDetails(p1)
        if OptimizeLoop(p1) == Dominating
            flag == 1 &&
                push!(temp_Bpoints, B_points([p1.x_cost, p2.objective(p1.x)],p1.x))
            flag == 2 &&
                push!(temp_Bpoints, B_points([p2.objective(p1.x), p1.x_cost],p1.x))
        end
    end

    Finish(p1)

    if p1.status.optimization_status != MeshPrecisionLimit && p1.status.optimization_status != PollPrecisionLimit && p1.status.optimization_status != KeyInterrupt
        error("The $(p1.status.optimization_status_string) is not enough for BiMADS, please increase it")
    end
    return temp_Bpoints
end

"""
    pareto_front(points::Vector{B_points}, minimize = true)

find the Pareto points, only for this BiMADS.jl, not general methods
"""
function pareto_front(points::Vector{B_points}, minimize = true)
    cmp_strict = minimize ? (<) : (>)
    cmp = minimize ? (<=) : (>=)
    dims = 2 #number of objective functions, needs modification for multi-MADS
    strictly_dominates(u::B_points, v::B_points) =
        all(cmp(u.cost[i], v.cost[i]) for i = 1:dims) &&
        any(cmp_strict(u.cost[i], v.cost[i]) for i = 1:dims)
    undominated(p::B_points) =
        !any(
            strictly_dominates(p2::B_points, p::B_points) for
            p2 ∈ points if p2.cost != p.cost
        )
    filter!(undominated, points)
    sort!(points, by = v -> v.cost, rev = false)
    return points
end

"""
    initial_X_L(p1, p2)::Vector{Float64}

get the initial undominated points from p1 and p2
"""
function initial_X_L(p1, p2)::Vector{B_points}
    x_L = Vector{B_points}()
    # optmization for p1
    append!(x_L, get_Bpoints(p1, p2, 1))
    x_L = reverse(x_L)
    # optmization for p2
    append!(x_L, get_Bpoints(p2, p1, 2))
    return x_L
end

"""
    get_ref(points::Vector{B_points})::Int

get the reference point for L>2
"""
function get_ref(points::Vector{B_points})::Int
    dis(x::Int) = LinearAlgebra.norm(points[x].cost - points[x-1].cost)^2+LinearAlgebra.norm(points[x].cost - points[x+1].cost)^2
    max_index = 2 #point with largest distance to its adjacent points
    max_distance=0.

    for i = 2:length(points)-1
        if dis(i)/(points[i].weight+1) >= max_distance
            max_distance = dis(i)/(points[i].weight+1)
            max_index = i
        end
    end
    return max_index
end

"""
    ReferencePointDetermination(undominated_points::Vector{B_points})::Tuple{Int,Float64,Vector{Float64}}

get the reference point for any case
"""
function ReferencePointDetermination(undominated_points::Vector{B_points})::Tuple{Int,Float64,Vector{Float64}}
    L=length(undominated_points)
    if L==1
        j=1
        ref_point=undominated_points[1].cost
        δ= 1.
    elseif L==2
        j=1 #different definition in paper(j=2) and books(j=1)
        ref_point=[undominated_points[2].cost[1],undominated_points[1].cost[2]]
        δ=(LinearAlgebra.norm(undominated_points[1].cost - undominated_points[2].cost)^2)/(undominated_points[2].weight+1)
    else
        j = get_ref(undominated_points)
        ref_point=[undominated_points[j+1].cost[1],undominated_points[j-1].cost[2]]
        δ=(LinearAlgebra.norm(undominated_points[j].cost - undominated_points[j-1].cost)^2
        +LinearAlgebra.norm(undominated_points[j].cost - undominated_points[j+1].cost)^2)/(undominated_points[j].weight+1)
    end
    undominated_points[j].weight+=1
    return j, δ, ref_point
end

"""
    function phi(f1::Function,f2::Function, r::Vector{Float64}, x::Vector{Float64})

Auxiliary Function from BiMADS. Aim to refomulate the Bi-obj problem into
single-obj problem. Return the cost of the refomulated objective function

x is the point to be evaluated, r is the reference point
"""
function phi(f1::Function,f2::Function, r::Vector{Float64}, x::Vector{Float64})
    cost1=f1(x)
    cost2=f2(x)
    if (cost1 <= r[1]) && (cost2 <= r[2]) #if x dominant r
        return -(r[1] - cost1)^2 * (r[2] - cost2)^2
    elseif (cost1 == r[1]) || (cost2 == r[2])
        return 0
    else
        return (max(cost1 - r[1], 0))^2 + (max(cost2 - r[2], 0))^2
    end
end

"""
    update_p(p::DSProblem)
reset the reformulated problem for new MADS iteration
Including reset the Mesh/Poll size, IterationLimit, FunctionEvaluationLimit, RuntimeLimit, etc.

"""
# update p before the main iteration
function update_p(p,p1::DSProblem,p2::DSProblem,status::BiMADS_status)
    for c in p.stoppingconditions
        i=_get_conditionindexes(p,typeof(c))[1]
        if c isa IterationStoppingCondition
            status.iteration=p1.status.iteration+p2.status.iteration
            p.stoppingconditions[i].limit-=status.iteration
            p.stoppingconditions[i].limit<=1 && error("The Iteration limit is not enough for BiMADS, please increase it")
        elseif c isa FunctionEvaluationStoppingCondition
            status.func_evaluation=p1.status.function_evaluations+p2.status.function_evaluations
            p.stoppingconditions[i].limit-=status.func_evaluation
            p.stoppingconditions[i].limit<=1 && error("The Function evaluation limit is not enough for BiMADS, please increase it")
        elseif c isa RuntimeStoppingCondition
            status.total_time=p1.status.runtime_total+p2.status.runtime_total
            p.stoppingconditions[i].limit -=status.total_time
            p.stoppingconditions[i].limit<=0 && error("The Runtime limit is not enough for BiMADS, please increase it")
        end
    end
end

# update p during the main iteration
function update_p(p_init::DSProblem,p_reform::DSProblem,status::BiMADS_status)
    status.iteration+=p_reform.status.iteration
    status.func_evaluation+=p_reform.status.function_evaluations
    status.total_time=p_reform.status.runtime_total
    p_reform=deepcopy(p_init)

    for c in p_reform.stoppingconditions
        i=_get_conditionindexes(p_reform,typeof(c))[1]
        if c isa IterationStoppingCondition
            p_reform.stoppingconditions[i].limit=p_init.stoppingconditions[i].limit-status.iteration
        elseif c isa FunctionEvaluationStoppingCondition
            p_reform.stoppingconditions[i].limit=p_init.stoppingconditions[i].limit-status.func_evaluation
        end
    end
    return p_reform
end
# update p after the main iteration
# function update_p(p_reform::DSProblem,status::BiMADS_status)
#     for c in p.stoppingconditions
#         i=_get_conditionindexes(p,typeof(c))[1]
#         if c isa IterationStoppingCondition
#             status.iteration+=p_reform.status.iteration
#         elseif c isa FunctionEvaluationStoppingCondition
#             status.func_evaluation+=p_reform.status.function_evaluations
#         end
#     end
# end

"""
    checkKeyInterrupt(status::BiMADS_status)::Bool
Check the input from REPL and stop the optimization

"""
# for BiMADS
function checkKeyInterrupt(status::BiMADS_status)::Bool
    bb = bytesavailable(stdin)
    data = read(stdin, bb)
    if bb>0 && data==UInt8[0x71] #q for quit
        println("quit")
        status.opt_status=KeyInterrupt
        return false
    end
    return true
end

"""
    checkBiMADSStopping(p::DSProblem,status::BiMADS_status,undominated_points::Vector{B_points})::Bool
Check the stopping condition for BiMADS

"""
function checkBiMADSStopping(p::DSProblem,status::BiMADS_status,undominated_points::Vector{B_points})::Bool
    p.status.runtime_total = time() - status.start_time
    if p.status.optimization_status!=PollPrecisionLimit && p.status.optimization_status!=MeshPrecisionLimit
        # add the numbers from last iteration
        status.opt_status=p.status.optimization_status
        status.iteration+=p.status.iteration
        status.func_evaluation+=p.status.function_evaluations
        return false
    end

    for c in p.stoppingconditions
        if c isa RuntimeStoppingCondition
            status.opt_status=RuntimeLimit
            i=_get_conditionindexes(p,RuntimeStoppingCondition)[1]
            p.status.runtime_total>p.stoppingconditions[i].limit && return false
        end
        if c isa HypervolumeStoppingCondition
            status.opt_status=HypervolumeLimit
            i=_get_conditionindexes(p,HypervolumeStoppingCondition)[1]
            diff=hvIndicator(undominated_points)-status.hypervolume
            status.hypervolume=hvIndicator(undominated_points)
            status.hypervolume>=p.stoppingconditions[i].limit && return false #if the HV smaller than a certain value
            # diff<=p.stoppingconditions[i].limit && return false #if the change of HV smaller than threshold
        end
    end

    !checkKeyInterrupt(status) && return false
    return true
end

"""
    Optimize_Bi!(p::DSProblem)

Run the BiMADS algorithm on problem `p`.

`p` must have had its initial point and objective function set. If extreme
barrier constraints have been set then the initial point must be value for
those constraints.
"""
function Optimize_Bi!(p::DSProblem)
    println("BiMADS")
    status=BiMADS_status()
    p1::DSProblem, p2::DSProblem = p_split(p)
    f1 = p1.objective
    f2 = p2.objective

# SetIterationLimit(p1,10)
# SetIterationLimit(p2,2)

    # Initialization
    undominated_points = pareto_front(initial_X_L(p1, p2))
# display(undominated_points)
# display(plot_Bpoint(undominated_points))

    p_reform=deepcopy(p)
    update_p(p_reform,p1,p2,status)

    iteration_count=0
    ref_point=Vector{Float64}()
    j = 0     #initial point for the refomulated sigle-obj problem
    δ = 0.0 # quantify the quality of the coverage of the undominated points
            # here is the coverage of the point with largest distance to its adjacent points
# fig1=plot_Bpoint(undominated_points)
# fig2=plot_Bpoint(undominated_points)
    # Main iteration
    while true
        count=0
# count2=0
# eva=500
# tti=10000
# length(undominated_points)>200 && break
        #Reference point determination
        j, δ, ref_point=ReferencePointDetermination(undominated_points)

        #Single-objective formulation minimization
        f_reform(x)=phi(f1,f2, ref_point, x)
        SetObjective(p_reform, f_reform)
        SetInitialPoint(p_reform,undominated_points[j].x_now)
# SetIterationLimit(p_reform,10)
# fig1=plot_Bpoint(undominated_points)
# p_reform.full_output=true
        #run MADS for refomulated problem and add new undominated points
        Setup(p_reform)

        while _check_stoppingconditions(p_reform)
            p_reform.full_output && OutputIterationDetails(p_reform)
            if OptimizeLoop(p_reform) == Dominating
                    push!(undominated_points, B_points([p1.objective(p_reform.x), p2.objective(p_reform.x)],p_reform.x))
                    count+=1
# count2+=1

# p_reform.status.function_evaluations>eva && break
# p_reform.status.optimization_status=PollPrecisionLimit
# count>3000 && break
# plot!(fig2,[p1.objective(p_reform.x)],[p2.objective(p_reform.x)],seriestype = :scatter,color=logocolors.blue)
            end
        end
        # p_reform.status.runtime_total = time() - status.start_time
# fig2=plot_Bpoint(undominated_points)

# display(fig2)

        #update X_L
        pareto_front(undominated_points)
        println("===============================")
        println("Start of iteration ",iteration_count)
        println("j=",j,"  ref=",ref_point)
        println("Add new points:",count)
        println("Hyper-Volume:",hvIndicator(undominated_points))
        println("Total undominated points: ", length(undominated_points))
        # println("===============================")
        iteration_count+=1
# iteration_count>tti && break
# @show p_reform.status.optimization_status
# p_reform=update_p(p,p_reform,status)
        !checkBiMADSStopping(p_reform,status,undominated_points) && break
        p_reform=update_p(p,p_reform,status)
# length(undominated_points)>=474 && break
        # n_iteration>=p.stoppingconditions[1].limit ? (n_iteration+=1) : (n_iteration+=1)
        # display(undominated_points)
        # fig3=plot_Bpoint(undominated_points)
        # figall=plot(fig2,fig3,aspect_ratio=1)
        # savefig(figall, "/Users/zyy/Desktop/XJTLU/MSc_Project/Julia/test_julia/Results/DS_result_$(iteration_count-1).pdf");

    end
    # update_p(p_reform,status)
    println("===============================")
    println("Total undominated points:", length(undominated_points))
    println("Total Iterations: ",status.iteration)
    println("Total Function Evaluations: ",status.func_evaluation)
    println("Total Run Time: ",p_reform.status.runtime_total)
    println("Hyper-Volume: ",hvIndicator(undominated_points))
    println("Optimization Status: ",status.opt_status)
    println("===============================")
    return undominated_points
end

"""
    p_dim(obj::Function,p.N)::Int

The dimension of the given objective function
If the dimension==1, reset the objective function for MADS
"""
function p_dim(p::DSProblem)::Int
    dim=length(p.objective(ones(p.N)))
    if dim==1
        x = zeros(Float64, p.N)
        f= p.objective(x)
        SetObjective(p, f[1])
        return dim
    end
    return dim
end


text:
using SPICE
using Statistics: mean

file = "pck00010.tpc"
furnsh(file)

ids = Set([])
re = r"BODY(?<id>[0-9]+)_PM\s+="
lines = open(readlines, file)
for line in lines
    s = string(line)
    if occursin(re, s)
        m = match(re, s)
        push!(ids, m["id"])
    end
end

open("pck.jl", "w") do f
    write(f, "# Automatically generated by $(basename(@__FILE__)), do not edit!\n\n")
    for id in ids
        body = replace(titlecase(bodc2n(parse(Int, id))), r"\s"=>"")
        body = body == "Moon" ? "Luna" : body

        # Radii
        try
            r = bodvrd(id, "RADII")
            subplan, along, polar = r
            write(f, "subplanetary_radius(::Type{Float64}, ::$body) = $subplan\n")
            write(f, "subplanetary_radius(::$body) = subplanetary_radius(Float64, $body())\n")
            write(f, "along_orbit_radius(::Type{Float64}, ::$body) = $along\n")
            write(f, "along_orbit_radius(::$body) = along_orbit_radius(Float64, $body())\n")
            if subplan == along
                write(f, "equatorial_radius(::Type{Float64}, ::$body) = $subplan\n")
                write(f, "equatorial_radius(::$body) = equatorial_radius(Float64, $body())\n")
            end
            write(f, "polar_radius(::Type{Float64}, ::$body) = $polar\n")
            write(f, "polar_radius(::$body) = polar_radius(Float64, $body())\n")
            write(f, "mean_radius(::Type{Float64}, ::$body) = $(mean(r))\n")
            write(f, "mean_radius(::$body) = mean_radius(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end

        # Right ascension
        try
            ra = deg2rad.(bodvrd(id, "POLE_RA"))
            n = length(ra)
            ra1 = zeros(3)
            ra1[1:n] .+= ra[1:n]
            alpha0, alpha1, alpha2 = ra1
            write(f, "alpha0(::Type{Float64}, ::$body) = $alpha0\n")
            write(f, "alpha0(::$body) = alpha0(Float64, $body())\n")
            write(f, "alpha1(::Type{Float64}, ::$body) = $alpha1\n")
            write(f, "alpha1(::$body) = alpha1(Float64, $body())\n")
            write(f, "alpha2(::Type{Float64}, ::$body) = $alpha2\n")
            write(f, "alpha2(::$body) = alpha2(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end
        try
            n = length(bodvrd(id[1:1], "NUT_PREC_ANGLES")) ÷ 2
            ra = zeros(n)
            ra′ = deg2rad.(bodvrd(id, "NUT_PREC_RA"))
            ra[1:length(ra′)] .+= ra′
            write(f, "alpha(::Type{Float64}, ::$body) = [$(join(ra, ", "))]\n")
            write(f, "alpha(::$body) = alpha(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end

        # Declination
        try
            dec = deg2rad.(bodvrd(id, "POLE_DEC"))
            n = length(dec)
            dec1 = zeros(3)
            dec1[1:n] .+= dec[1:n]
            delta0, delta1, delta2 = dec1
            write(f, "delta0(::Type{Float64}, ::$body) = $delta0\n")
            write(f, "delta0(::$body) = delta0(Float64, $body())\n")
            write(f, "delta1(::Type{Float64}, ::$body) = $delta1\n")
            write(f, "delta1(::$body) = delta1(Float64, $body())\n")
            write(f, "delta2(::Type{Float64}, ::$body) = $delta2\n")
            write(f, "delta2(::$body) = delta2(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end
        try
            n = length(bodvrd(id[1:1], "NUT_PREC_ANGLES")) ÷ 2
            dec = zeros(n)
            dec′ = deg2rad.(bodvrd(id, "NUT_PREC_DEC"))
            dec[1:length(dec′)] .+= dec′
            write(f, "delta(::Type{Float64}, ::$body) = [$(join(dec, ", "))]\n")
            write(f, "delta(::$body) = delta(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end

        # Rotation
        try
            pm = deg2rad.(bodvrd(id, "PM"))
            n = length(pm)
            pm1 = zeros(3)
            pm1[1:n] .+= pm[1:n]
            omega0, omega1, omega2 = pm1
            write(f, "omega0(::Type{Float64}, ::$body) = $omega0\n")
            write(f, "omega0(::$body) = omega0(Float64, $body())\n")
            write(f, "omega1(::Type{Float64}, ::$body) = $omega1\n")
            write(f, "omega1(::$body) = omega1(Float64, $body())\n")
            write(f, "omega2(::Type{Float64}, ::$body) = $omega2\n")
            write(f, "omega2(::$body) = omega2(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end
        try
            n = length(bodvrd(id[1:1], "NUT_PREC_ANGLES")) ÷ 2
            pm = zeros(n)
            pm′ = deg2rad.(bodvrd(id, "NUT_PREC_PM"))
            pm[1:length(pm′)] .+= pm′
            write(f, "omega(::Type{Float64}, ::$body) = [$(join(pm, ", "))]\n")
            write(f, "omega(::$body) = omega(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end

        # Nutation / Precession
        try
            length(id) != 3 && continue
            np = deg2rad.(reshape(bodvrd(id[1:1], "NUT_PREC_ANGLES"), 2, :))
            write(f, "theta0(::Type{Float64}, ::$body) = [$(join(np[1,:], ", "))]\n")
            write(f, "theta0(::$body) = theta0(Float64, $body())\n")
            write(f, "theta1(::Type{Float64}, ::$body) = [$(join(np[2,:], ", "))]\n")
            write(f, "theta1(::$body) = theta1(Float64, $body())\n")
        catch err
            err isa SpiceError || rethrow(err)
        end
    end
end


text:
# ---
# title: JuliaRL\_IQN\_CartPole
# cover: assets/JuliaRL_IQN_CartPole.png
# description: IQN applied to CartPole
# date: 2021-05-22
# author: "[Jun Tian](https://github.com/findmyway)"
# ---

#+ tangle=true
using ReinforcementLearning
using StableRNGs
using Flux
using Flux.Losses

function RL.Experiment(
    ::Val{:JuliaRL},
    ::Val{:IQN},
    ::Val{:CartPole},
    ::Nothing;
    seed = 123,
)
    rng = StableRNG(seed)
    env = CartPoleEnv(; T = Float32, rng = rng)
    ns, na = length(state(env)), length(action_space(env))
    init = glorot_uniform(rng)
    Nₑₘ = 16
    n_hidden = 64
    κ = 1.0f0

    nn_creator() =
        ImplicitQuantileNet(
            ψ = Dense(ns, n_hidden, relu; init = init),
            ϕ = Dense(Nₑₘ, n_hidden, relu; init = init),
            header = Dense(n_hidden, na; init = init),
        ) |> cpu

    agent = Agent(
        policy = QBasedPolicy(
            learner = IQNLearner(
                approximator = NeuralNetworkApproximator(
                    model = nn_creator(),
                    optimizer = ADAM(0.001),
                ),
                target_approximator = NeuralNetworkApproximator(model = nn_creator()),
                κ = κ,
                N = 8,
                N′ = 8,
                Nₑₘ = Nₑₘ,
                K = 32,
                γ = 0.99f0,
                stack_size = nothing,
                batch_size = 32,
                update_horizon = 1,
                min_replay_history = 100,
                update_freq = 1,
                target_update_freq = 100,
                default_priority = 1.0f2,
                rng = rng,
                device_rng = rng,
            ),
            explorer = EpsilonGreedyExplorer(
                kind = :exp,
                ϵ_stable = 0.01,
                decay_steps = 500,
                rng = rng,
            ),
        ),
        trajectory = CircularArrayPSARTTrajectory(
            capacity = 1000,
            state = Vector{Float32} => (ns,),
        ),
    )

    stop_condition = StopAfterStep(10_000, is_show_progress=!haskey(ENV, "CI"))
    hook = TotalRewardPerEpisode()
    Experiment(agent, env, stop_condition, hook, "")
end

#+ tangle=false
using Plots
pyplot() #hide
ex = E`JuliaRL_IQN_CartPole`
run(ex)
plot(ex.hook.rewards)
savefig("assets/JuliaRL_IQN_CartPole.png") #hide

# ![](assets/JuliaRL_IQN_CartPole.png)


text:
module clj_julia_helper
"Ruthlessly stolen from pyjulia"

"""
    num_utf8_trailing(d::Vector{UInt8})

If `d` ends with an incomplete UTF8-encoded character, return the number of trailing incomplete bytes.
Otherwise, return `0`.

Taken from IJulia.jl.
"""
function num_utf8_trailing(d::Vector{UInt8})
    i = length(d)
    # find last non-continuation byte in d:
    while i >= 1 && ((d[i] & 0xc0) == 0x80)
        i -= 1
    end
    i < 1 && return 0
    c = d[i]
    # compute number of expected UTF-8 bytes starting at i:
    n = c <= 0x7f ? 1 : c < 0xe0 ? 2 : c < 0xf0 ? 3 : 4
    nend = length(d) + 1 - i # num bytes from i to end
    return nend == n ? 0 : nend
end

function pipe_stream(sender::IO, receiver, buf::IO = IOBuffer())
    receiver("begin pipe stream")
    try
        while !eof(sender)
            receiver("in loop")
            nb = bytesavailable(sender)
            write(buf, read(sender, nb))

            # Taken from IJulia.send_stream:
            d = take!(buf)
            n = num_utf8_trailing(d)
            dextra = d[end-(n-1):end]
            resize!(d, length(d) - n)
            s = String(copy(d))

            write(buf, dextra)
            receiver(s)  # check isvalid(String, s)?
        end
    catch e
        receiver("error detected")
        if !isa(e, InterruptException)
            rethrow()
        end
        pipe_stream(sender, receiver, buf)
    end
end

const read_stdout = Ref{Base.PipeEndpoint}()
const read_stderr = Ref{Base.PipeEndpoint}()

function pipe_std_outputs(out_receiver, err_receiver)
    global readout_task
    global readerr_task
    global out_rec = out_receiver
    global err_rec = err_receiver
    read_stdout[], = redirect_stdout()
    readout_task = @task pipe_stream(read_stdout[], out_receiver)
    read_stderr[], = redirect_stderr()
    readerr_task = @task pipe_stream(read_stderr[], err_receiver)
    (readout_task,readerr_task)
end

end

max_stars_repo_name: cnuernber/libjulia-clj

max_issues_repo_name: cnuernber/julia-clj

max_forks_repo_name: cnuernber/julia-clj


text:
# ## Ghost
set_camera!(vis, cam_pos=[-1,1,0], zoom=1)

z_sim = get_maximal_state(storage)
timesteps = [5, 10, 15]# .+ 150

for t in timesteps
    name = Symbol("robot_$t")
    build_robot(mech, vis=vis, name=name, color= magenta_light)
    z = z_sim[t]
    set_robot(vis, mech, z, name=name)
end

z_sim[220][1] -= 0.1
z_sim[220][13 + 1] -= 0.05
set_robot(vis, mech, z_sim[220])

z_sim[104][1] -= 0.025
z_sim[104][13 + 1] -= 0.05
set_robot(vis, mech, z_sim[104])


text:
struct StrangSplitting{T,R,V} <: AbstractIntegrator{T}
    uv::Tuple{V,V}
    E::Vector{<:AbstractExponentiator{T}}
    function StrangSplitting(u::V, E::AbstractExponentiator{T}...) where {T,V}
        foreach(e -> checkorder(e,StrangSplitting), E)
        new{T,real(T),V}((u,similar(u)), Vector{AbstractExponentiator{T}}(E...))
    end
end

order(::StrangSplitting) = 2
order(::Type{StrangSplitting}) = 2

function step!(strang::StrangSplitting{T,R,V}, t::R, Δt::R) where {T,R,V}
    # Integrate all operators over time-step, to second order since
    # the Strang splitting itself introduces a second-order splitting
    # error.
    foreach(e -> midpoint!(operator(e), t, Δt), strang.E)

    a,b = 1,2

    # When the operator is L = A + B + C + D + ..., the Strang
    # splitting is ... (D/2) (C/2) (B/2) A (B/2) (C/2) (D/2) ...

    # ... (D/2) (C/2) (B/2)
    for E in reverse(strang.E[2:end])
        expv!(strang.uv[b], E, Δt/2, strang.uv[a], -1)
        inplace(E) || @swap! a b
    end

    # A
    expv!(strang.uv[b], E₁, Δt, strang.uv[a], 0)
    inplace(strang.E[1]) || @swap! a b

    # (B/2) (C/2) (D/2) ...
    for E in strang.E[2:end]
        expv!(strang.uv[b], E, Δt/2, strang.uv[a], 1)
        inplace(E) || @swap! a b
    end

    # If the propagated solution is in the v vector, copy it to the u
    # vector.
    a == 2 && copyto!(strang.uv[1], strang.uv[2])
end


text:
module Citations

using Documenter, Markdown

import Documenter: Expanders, Builder, Documents, Anchors
import Documenter.Utilities: Selectors
import Documenter.Utilities.DOM: @tags

function read_bib_file(filename)
    entries = open(filename) do io
        contents = read(filename, String)

        parser = btp[:bparser][:BibTexParser](common_strings = true)
        parser[:ignore_nonstandard_types] = true
        parser[:interpolate_strings] = true

        parsed = btp[:loads](contents)
        parsed[:entries]
    end

    citations = Dict{String, Any}()
    for entry in entries
        citations[entry["ID"]] = entry
    end

    return(citations)
end

include("BibliographyBlocks.jl")
include("nodes.jl")
include("GenerateBibliography.jl")
include("HTMLWriter.jl")
include("RegisterCitations.jl")
include("walk_and_replace.jl")

end


text:
# ------------------------------------------------------------------------------

function check_primality_zerodim(J::Singular.sideal{Singular.spoly{Singular.n_Q}})
    J = Singular.std(J)
    basis = gens(Singular.kbase(J))
    dim = length(basis)
    S = Nemo.MatrixSpace(Nemo.QQ, dim, dim)
    matrices = []

    for v in gens(base_ring(J))
        M = zero(S)
        for (i, vec) in enumerate(basis)
            image = Singular.reduce(v * vec, J)
            for (j, base_vec) in enumerate(basis)
                M[i, j] = Nemo.QQ(coeff(image, base_vec))
            end
        end
        push!(matrices, M)
    end
    generic_multiplication = sum(rand(1:100) * M for M in matrices)
    return isirreducible(Nemo.charpoly(generic_multiplication))
end

#------------------------------------------------------------------------------
"""
    check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly}, extra_relations::Array{fmpq_mpoly, 1})

The function checks if the ideal generated by the polynomials and saturated at the leading coefficient with respect to the corresponding variables is prime over rationals.

The `extra_relations` allows to add more polynomials to the generators (not affecting the saturation).
"""
function check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly}, extra_relations::Array{fmpq_mpoly, 1})
    leaders = collect(keys(polys))
    ring = parent(leaders[1])

    Rsing, vsing = Singular.PolynomialRing(Singular.QQ, [var_to_str(l) for l in leaders])
    eval_point = [v in keys(polys) ? v : ring(rand(1:100)) for v in gens(ring)]
    all_polys = vcat(collect(values(polys)), extra_relations)
    zerodim_ideal = Singular.Ideal(Rsing, map(p -> parent_ring_change(evaluate(p, eval_point), Rsing), all_polys))

    return check_primality_zerodim(zerodim_ideal)
end

#------------------------------------------------------------------------------
"""
    check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly})

The function checks if the ideal generated by the polynomials and saturated at the leading coefficient with respect to the corresponding variables is prime over rationals.
"""
function check_primality(polys::Dict{fmpq_mpoly, fmpq_mpoly})
    return check_primality(polys, Array{fmpq_mpoly, 1}())
end

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