text:
"A `FunctionLevelSet` is a set that derives from the levels of a function."
abstract type FunctionLevelSet{T} <: Domain{T} end

# Convenience: assume the function and the level are stored fields
levelfun(d::FunctionLevelSet) = d.f
levelfun(d::FunctionLevelSet, x) = d.f(x)
level(d::FunctionLevelSet) = d.level

"Supertype of level set domains of the form `f(x)=C`."
abstract type AbstractLevelSet{T} <: FunctionLevelSet{T} end

indomain(x, d::AbstractLevelSet) = levelfun(d, x) == level(d)

show(io::IO, d::AbstractLevelSet) =
    print(io, "level set f(x) = $(level(d)) with f = $(levelfun(d))")

==(d1::AbstractLevelSet, d2::AbstractLevelSet) = levelfun(d1)==levelfun(d2) &&
    level(d1)==level(d2)
hash(d::AbstractLevelSet, h::UInt) = hashrec("AbstractLevelSet", levelfun(d), level(d), h)

"The domain defined by `f(x)=0` for a given function `f`."
struct ZeroSet{T,F} <: AbstractLevelSet{T}
    f   ::  F
end

ZeroSet(f) = ZeroSet{Float64}(f)
ZeroSet{T}(f::F) where {T,F} = ZeroSet{T,F}(f)

level(d::ZeroSet) = 0

similardomain(d::ZeroSet, ::Type{T}) where {T} = ZeroSet{T}(levelfun(d))

"The domain defined by `f(x)=C` for a given function `f` and constant `C`."
struct LevelSet{T,F,S} <: AbstractLevelSet{T}
    f       ::  F
    level   ::  S
end

LevelSet(f, level) = LevelSet{typeof(level)}(f, level)
LevelSet{T}(f::F, level::S) where {T,F,S} = LevelSet{T,F,S}(f, level)

similardomain(d::LevelSet, ::Type{T}) where {T} = LevelSet{T}(levelfun(d), level(d))

convert(::Type{LevelSet}, d::ZeroSet{T}) where {T} = LevelSet{T}(levelfun(d), level(d))
convert(::Type{LevelSet{T}}, d::ZeroSet) where {T} = LevelSet{T}(levelfun(d), level(d))

"Supertype of sublevel set domains."
abstract type AbstractSublevelSet{T,C} <: FunctionLevelSet{T} end

indomain(x, d::AbstractSublevelSet{T,:closed}) where {T} = levelfun(d, x) <= level(d)
indomain(x, d::AbstractSublevelSet{T,:open}) where {T} = levelfun(d, x) < level(d)

show(io::IO, d::AbstractSublevelSet{T,:closed}) where {T} =
    print(io, "sublevel set f(x) <= $(level(d)) with f = $(levelfun(d))")
show(io::IO, d::AbstractSublevelSet{T,:open}) where {T} =
    print(io, "sublevel set f(x) < $(level(d)) with f = $(levelfun(d))")

==(d1::AbstractSublevelSet, d2::AbstractSublevelSet) = levelfun(d1)==levelfun(d2) &&
    level(d1)==level(d2)
hash(d::AbstractSublevelSet, h::UInt) =
    hashrec("AbstractSublevelSet", levelfun(d), level(d), h)

"The domain where `f(x) <= 0` (or `f(x) < 0`)."
struct SubzeroSet{T,C,F} <: AbstractSublevelSet{T,C}
    f   ::  F
end

SubzeroSet(f) = SubzeroSet{Float64}(f)
SubzeroSet{T}(f) where {T} = SubzeroSet{T,:closed}(f)
SubzeroSet{T,C}(f::F) where {T,C,F} = SubzeroSet{T,C,F}(f)

level(d::SubzeroSet) = 0

similardomain(d::SubzeroSet{S,C}, ::Type{T}) where {S,C,T} =
    SubzeroSet{T,C}(levelfun(d))

interior(d::SubzeroSet{T}) where {T} = SubzeroSet{T,:open}(d.f)
closure(d::SubzeroSet{T}) where {T} = SubzeroSet{T,:closed}(d.f)

boundary(d::SubzeroSet{T}) where {T} = ZeroSet{T}(d.f)

"The domain defined by `f(x) <= C` (or `f(x) < C`) for a given function `f` and constant `C`."
struct SublevelSet{T,C,F,S} <: AbstractSublevelSet{T,C}
    f       ::  F
    level   ::  S
end

SublevelSet(f, level) = SublevelSet{typeof(level)}(f, level)
SublevelSet{T}(f, level) where {T} = SublevelSet{T,:closed}(f, level)
SublevelSet{T,C}(f::F, level::S) where {T,C,F,S} = SublevelSet{T,C,F,S}(f, level)

similardomain(d::SublevelSet{S,C}, ::Type{T}) where {S,C,T} =
    SublevelSet{T,C}(levelfun(d), level(d))

interior(d::SublevelSet{T}) where {T} = SublevelSet{T,:open}(d.f, d.level)
closure(d::SublevelSet{T}) where {T} = SublevelSet{T,:closed}(d.f, d.level)

boundary(d::SublevelSet{T}) where {T} = LevelSet{T}(d.f, d.level)

"Supertype of superlevel set domains."
abstract type AbstractSuperlevelSet{T,C} <: FunctionLevelSet{T} end

indomain(x, d::AbstractSuperlevelSet{T,:closed}) where {T} = levelfun(d, x) >= level(d)
indomain(x, d::AbstractSuperlevelSet{T,:open}) where {T} = levelfun(d, x) > level(d)

show(io::IO, d::AbstractSuperlevelSet{T,:closed}) where {T} =
    print(io, "superlevel set f(x) >= $(level(d)) with f = $(levelfun(d))")
show(io::IO, d::AbstractSuperlevelSet{T,:open}) where {T} =
    print(io, "superlevel set f(x) > $(level(d)) with f = $(levelfun(d))")

==(d1::AbstractSuperlevelSet, d2::AbstractSuperlevelSet) =
    levelfun(d1)==levelfun(d2) && level(d1)==level(d2)
hash(d::AbstractSuperlevelSet, h::UInt) =
    hashrec("AbstractSuperlevelSet", levelfun(d), level(d), h)

"The domain where `f(x) >= 0` (or `f(x) > 0`)."
struct SuperzeroSet{T,C,F} <: AbstractSuperlevelSet{T,C}
    f   ::  F
end

SuperzeroSet(f) = SuperzeroSet{Float64}(f)
SuperzeroSet{T}(f) where {T} = SuperzeroSet{T,:closed}(f)
SuperzeroSet{T,C}(f::F) where {T,C,F} = SuperzeroSet{T,C,F}(f)

level(d::SuperzeroSet) = 0

similardomain(d::SuperzeroSet{S,C}, ::Type{T}) where {S,C,T} =
    SuperzeroSet{T,C}(levelfun(d))

interior(d::SuperzeroSet{T}) where {T} = SuperzeroSet{T,:open}(d.f)
closure(d::SuperzeroSet{T}) where {T} = SuperzeroSet{T,:closed}(d.f)

boundary(d::SuperzeroSet{T}) where {T} = ZeroSet{T}(d.f)

"The domain defined by `f(x) >= C` (or `f(x) > C`) for a given function `f` and constant `C`."
struct SuperlevelSet{T,C,F,S} <: AbstractSuperlevelSet{T,C}
    f       ::  F
    level   ::  S
end

SuperlevelSet(f, level) = SuperlevelSet{typeof(level)}(f, level)
SuperlevelSet{T}(f, level) where {T} = SuperlevelSet{T,:closed}(f, level)
SuperlevelSet{T,C}(f::F, level::S) where {T,C,F,S} = SuperlevelSet{T,C,F,S}(f, level)

similardomain(d::SuperlevelSet{S,C}, ::Type{T}) where {S,C,T} =
    SuperlevelSet{T,C}(levelfun(d), level(d))

interior(d::SuperlevelSet{T}) where {T} = SuperlevelSet{T,:open}(d.f, d.level)
closure(d::SuperlevelSet{T}) where {T} = SuperlevelSet{T,:closed}(d.f, d.level)

boundary(d::SuperlevelSet{T}) where {T} = LevelSet{T}(d.f, d.level)

## Additional functionality

pseudolevel(d::AbstractLevelSet, epsilon) =
    _pseudolevel(d, epsilon, levelfun(d), level(d))
_pseudolevel(d::AbstractLevelSet{T}, epsilon, fun, C) where {T} =
    SublevelSet{T,:open}(x -> norm(fun(x)-C), epsilon)

pseudolevel(d::Domain, epsilon) = pseudolevel(convert(LevelSet, d), epsilon)


text:
using Statistics
using Random

struct ReactionSystem
    N::Integer
    M::Integer
    c::Vector{T} where T <: Real
    h::Vector{Function}
    update_X::Vector{Vector{Tuple{Integer,Integer}}}
    update_a::Vector{Vector{Integer}}
end

function gillespie_algorithm(rs::ReactionSystem, Xi; tmax = 100.0, nmax = 500)
    t = 0.0
    n = 1
    ts = Vector{typeof(t)}(undef, nmax)
    ts[1] = t
    X = copy(Xi)
    Xs = Array{typeof(Xi[1]),2}(undef, nmax, rs.N)
    Xs[1,:] = X
    a = [rs.c[ν]*(rs.h[ν](X)) for ν in 1:rs.M]
    a0 = sum(a)
    while t < tmax && n < nmax && a0 > 0
        τ = log(1/rand())/a0
        r = rand()
        μ = 1
        s = a[μ]
        while r*a0 > s
            μ += 1
            s += a[μ]
        end
        # update X
        for (i,j) in rs.update_X[μ]
            X[i] += j
        end
        # update a
        for ν in rs.update_a[μ]
            a[ν] = rs.c[ν]*(rs.h[ν](X))
        end
        a0 = sum(a) # recompute the entire sum to avoid rounding errors
        ts[n+1] = t
        Xs[n+1,:] = X
        t += τ
        n += 1
    end
    return ts[1:n], Xs[1:n,:]
end

function SI(α::Real)
    rs = ReactionSystem(2, 1, [α],
    [X -> X[1]*X[2]],
    [[(1,-1),(2,1)]],
    [[1]]
    )
end

function LotkaVolterra(c)
    rs = ReactionSystem(3,3,c,
    [X -> X[1]*X[3],
     X -> X[1]*X[2],
     X -> X[2]],
    [[(1,1)],[(1,-1),(2,1)],[(2,-1)]],
    [[1],[1,2,3],[2,3]]
    )
end

function ClosedLotkaVolterra(c)
    rs = ReactionSystem(3,3,c,
    [X -> X[1]*X[3],
     X -> X[1]*X[2],
     X -> X[2]],
    [[(1,1),(3,-1)],[(1,-1),(2,1)],[(2,-1)]],
    [[1],[1,2,3],[2,3]]
    )
end

function Brusselator(c)
    rs = ReactionSystem(4,4,c,
    [X -> X[3],
     X -> X[4]*X[1],
     X -> X[2]*X[1]*(X[1]-1)/2,
     X -> X[1]],
     [[(1,1)], [(1,-1),(2,1)], [(1,1),(2,-1)], [(1,-1)]],
     [[2,3,4], [2,3,4], [2,3,4], [2,3,4]]
    )
end

using Plots

si = SI(0.5)
ts, Xs = gillespie_algorithm(si, [100,1])
plot(ts, Xs[:,1], line=:steppre)
scatter!(ts, Xs[:,1], markersize=1)
plot!(ts, 100*exp.(-0.5*ts))

X0 = 100000
lv = LotkaVolterra([10/X0,0.01,10])

ts, Xs = gillespie_algorithm(lv, [1000,1000,X0], nmax = 10^6)

plot(ts, Xs[:,1], line=:steppre)
plot!(ts, Xs[:,2], line=:steppre)

plot(Xs[:,1],Xs[:,2], linewidth=0.1, label="")

X1 = 100000
X2 = 100000
br = Brusselator([5000/X0, 50/X2, 0.00005, 5])

ts, Xs = gillespie_algorithm(br, [1000,2000,X1,X2], tmax=10000.0, nmax=2*10^6)

plot(ts, Xs[:,1], line=:steppre)
plot!(ts, Xs[:,2], line=:steppre)

plot(Xs[:,1],Xs[:,2], linewidth=0.1, label="")

ext: jl

max_stars_repo_name: csimal/metaplex-networks

max_issues_repo_name: csimal/metaplex-networks

max_forks_repo_name: csimal/metaplex-networks


text:
function dgpnorm(n::Int, x::Array{Float64}, xl::Array{Float64},
                 xu::Array{Float64}, g::Array{Float64})
    inf_norm = 0.0
    for i=1:n
        if xl[i] != xu[i]
            if x[i] == xl[i]
                v = (min(g[i], 0.0))^2
            elseif x[i] == xu[i]
                v = (max(g[i], 0.0))^2
            else
                v = g[i]^2
            end

            v = sqrt(v)
            inf_norm = (inf_norm > v) ? inf_norm : v
        end
    end

    return inf_norm
end


text:
export StraightLineHomotopy, StraightLineHomotopyCache, gamma

"""
    StraightLineHomotopy(G, F; gamma=exp(i * 2π*rand()))

Construct the homotopy ``H(x, t) = γtG(x) + (1-t)F(x)``.
"""
struct StraightLineHomotopy{S<:AbstractSystem,T<:AbstractSystem} <: AbstractHomotopy
    start::S
    target::T
    gamma::Complex{Float64}

end
function StraightLineHomotopy(start::AbstractSystem, target::AbstractSystem; gamma=randomish_gamma())
    StraightLineHomotopy(start, target, gamma)
end
(H::StraightLineHomotopy)(x, t, c=cache(H, x, t)) = evaluate(H, x, t, c)

"""
    StraightLineHomotopyCache

An simple cache for `StartTargetHomotopy`s consisting of the caches for the
start and target system as well as a `Vector` and a `Matrix`.
"""
struct StraightLineHomotopyCache{SC, TC, T1, T2} <: AbstractHomotopyCache
    start::SC
    target::TC

    u::Vector{T1} # for evaluation
    U::Matrix{T2} # for Jacobian
end

function cache(H::StraightLineHomotopy, x, t)
    start_cache = cache(H.start, x)
    target_cache = cache(H.target, x)
    u = complex.(float.(evaluate(H.start, x, start_cache)))
    U = complex.(float.(jacobian(H.start, x, start_cache)))
    StraightLineHomotopyCache(start_cache, target_cache, u, U)
end

Base.size(H::StraightLineHomotopy) = size(H.start)

"""
    gamma(H::StraightLineHomotopy)

Obtain the gamma used in the StraightLineHomotopy.
"""
gamma(H::StraightLineHomotopy) = H.gamma

"""
    γ(H)

Obtain the gamma used in the StraightLineHomotopy.
"""
γ(H::StraightLineHomotopy) = gamma(H)

function evaluate!(u, H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    @inbounds evaluate!(c.u, H.start, x, c.start)
    @inbounds evaluate!(u, H.target, x, c.target)

    @inbounds for i in eachindex(u)
        u[i] = (γ(H) * t) * c.u[i] + (1.0 - t) * u[i]
    end

    u
end

function dt!(u, H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    @inbounds evaluate!(c.u, H.start, x, c.start)
    @inbounds evaluate!(u, H.target, x, c.target)

    @inbounds for i in eachindex(u)
        u[i] = γ(H) * c.u[i] - u[i]
    end
    u
end

function jacobian!(U, H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    @inbounds jacobian!(c.U, H.start, x, c.start)
    @inbounds jacobian!(U, H.target, x, c.target)

    @inbounds for i in eachindex(U)
        U[i] = (γ(H) * t) * c.U[i] + (1.0 - t) * U[i]
    end

    U
end

function evaluate_and_jacobian!(u, U, H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    @inbounds evaluate_and_jacobian!(c.u, c.U, H.start, x, c.start)
    @inbounds evaluate_and_jacobian!(u, U, H.target, x, c.target)

    @inbounds for i in eachindex(u)
        u[i] = (γ(H) * t) * c.u[i] + (1.0 - t) * u[i]
    end
    @inbounds for i in eachindex(U)
        U[i] = (γ(H) * t) * c.U[i] + (1.0 - t) * U[i]
    end

    nothing
end

function jacobian_and_dt!(U, u, H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    @inbounds evaluate_and_jacobian!(c.u, c.U, H.start, x, c.start)
    @inbounds evaluate_and_jacobian!(u, U, H.target, x, c.target)

    @inbounds for i in eachindex(U)
        U[i] = (γ(H) * t) * c.U[i] + (1.0 - t) * U[i]
    end
    @inbounds for i in eachindex(u)
        u[i] = γ(H) * c.u[i] - u[i]
    end

    nothing
end

function evaluate(H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    M, N = size(H)
    evaluate!(similar(c.u, M), H, x, t, c)
end
function jacobian(H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    M, N = size(H)
    jacobian!(similar(c.U, M, N), H, x, t, c)
end
function dt(H::StraightLineHomotopy, x, t, c::StraightLineHomotopyCache)
    M, N = size(H)
    dt!(similar(c.u, M), H, x, t, c)
end


text:
export GitIndex, GitIndexEntry, add_bypath!, write_tree!, write!, reload!, clear!,
       remove!, removedir!, read_tree!, add_all!, update_all!,
       remove_all!, has_conflicts, add!, is_fully_merged, conflicts

type GitIndex
    ptr::Ptr{Void}

    function GitIndex(ptr::Ptr{Void})
        @assert ptr != C_NULL
        idx = new(ptr)
        finalizer(idx, free!)
        return idx
    end
end

free!(idx::GitIndex) = begin
    if idx.ptr != C_NULL
        ccall((:git_index_free, libgit2), Void, (Ptr{Void},), idx.ptr)
        idx.ptr = C_NULL
    end
end

Base.convert(::Type{Ptr{Void}}, idx::GitIndex) = idx.ptr

GitIndex(path::String) = begin
    idxptr = Ptr{Void}[0]
    @check ccall((:git_index_open, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Uint8}), idxptr, path)
    return GitIndex(idxptr[1])
end

GitIndex(r::GitRepo) = begin
    idxptr = Ptr{Void}[0]
    @check ccall((:git_repository_index, libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}), idxptr, r)
    return GitIndex(idxptr[1])
end

has_conflicts(idx::GitIndex) =
    bool(ccall((:git_index_has_conflicts, libgit2), Cint, (Ptr{Void},), idx))

has_conflicts(idx::GitRepo) = has_conflicts(GitIndex(repo))

is_fully_merged(idx::GitIndex) = !has_conflicts(idx)
is_fully_merged(repo::GitRepo) = is_fully_merged(GitIndex(repo))

function clear!(idx::GitIndex)
    ccall((:git_index_clear, libgit2), Void, (Ptr{Void},), idx)
    return idx
end
clear!(repo::GitRepo) = clear!(GitIndex(repo))

function reload!(idx::GitIndex)
    @check ccall((:git_index_read, libgit2), Cint, (Ptr{Void}, Cint), idx, zero(Cint))
    return idx
end
reload!(repo::GitRepo) = reload!(GitIndex(repo))

function write!(idx::GitIndex)
    @check ccall((:git_index_write, libgit2), Cint, (Ptr{Void},), idx)
    return idx
end
write!(repo::GitRepo) = write!(GitIndex(repo))

function remove!(idx::GitIndex, path::String, stage::Integer=0)
    @check ccall((:git_index_remove, libgit2), Cint, (Ptr{Void}, Ptr{Uint8}, Cint),
                 idx, path, stage)
    return idx
end

function removedir!(idx::GitIndex, path::String, stage::Integer=0)
    @check ccall((:git_index_remove_directory, libgit2), Cint, (Ptr{Void}, Ptr{Uint8}, Cint),
                  idx, path, stage)
    return idx
end

Base.length(idx::GitIndex) =
    int(ccall((:git_index_entrycount, libgit2), Cint, (Ptr{Void},), idx))

function write_tree!(idx::GitIndex)
    id_ptr = [Oid()]
    @check ccall((:git_index_write_tree, libgit2), Cint,
                 (Ptr{Oid}, Ptr{Void}), id_ptr, idx)
    return id_ptr[1]
end

type GitIndexEntry
    path::ByteString
    id::Oid
    ctime::Float64
    mtime::Float64
    file_size::Int
    dev::Int
    ino::Int
    mode::Int
    uid::Int
    gid::Int
    valid::Bool
    stage::Int
end

immutable IndexTimeStruct
    seconds::GitTimeT
    nanoseconds::Cuint
end

IndexTimeStruct() = IndexTimeStruct(0, 0)

immutable IndexEntryStruct
    ctime::IndexTimeStruct
    mtime::IndexTimeStruct
    dev::Cuint
    ino::Cuint
    mode::Cuint
    uid::Cuint
    gid::Cuint
    file_size::GitOffT
    id::Oid
    flags::Cushort
    flags_extended::Cushort
    path::Ptr{Uint8}
end

IndexEntryStruct() = IndexEntryStruct(IndexTimeStruct(),
                                      IndexTimeStruct(),
                                      zero(Cuint),
                                      zero(Cuint),
                                      zero(Cuint),
                                      zero(Cuint),
                                      zero(Cuint),
                                      zero(GitOffT),
                                      Oid(),
                                      zero(Cushort),
                                      zero(Cushort),
                                      zero(Ptr{Uint8}))
Oid(entry::GitIndexEntry) = entry.id

IndexEntryStruct(entry::GitIndexEntry) = begin
    flags = uint16(0x0)
    flags &= ~GitConst.IDXENTRY_STAGEMASK
    flags |= (uint16(entry.stage) << GitConst.IDXENTRY_STAGESHIFT) & GitConst.IDXENTRY_STAGEMASK
    flags &= ~GitConst.IDXENTRY_VALID
    entry.valid && (flags |= GitConst.IDXENTRY_VALID)
    #TODO: rounding
    ctime = IndexTimeStruct(ifloor(entry.ctime), ifloor(entry.ctime - floor(entry.ctime) * 1e3))
    mtime = IndexTimeStruct(ifloor(entry.mtime), ifloor(entry.mtime - floor(entry.mtime) * 1e3))
    return IndexEntryStruct(ctime,
                            ctime,
                            entry.dev,
                            entry.ino,
                            entry.mode,
                            entry.uid,
                            entry.gid,
                            entry.file_size,
                            entry.id,
                            flags,
                            0,
                            convert(Ptr{Uint8}, entry.path))
end

function GitIndexEntry(ptr::Ptr{IndexEntryStruct})
    @assert ptr != C_NULL
    entry = unsafe_load(ptr)
    path  = bytestring(entry.path)
    ctime = entry.ctime.seconds + (entry.ctime.nanoseconds / 1e3)
    mtime = entry.mtime.seconds + (entry.mtime.nanoseconds / 1e3)
    dev   = int(entry.dev)
    ino   = int(entry.ino)
    mode  = int(entry.mode)
    uid   = int(entry.uid)
    gid   = int(entry.gid)
    valid = bool(entry.flags & GitConst.IDXENTRY_VALID)
    stage = int((entry.flags & GitConst.IDXENTRY_STAGEMASK) >> GitConst.IDXENTRY_STAGESHIFT)
    file_size = int(entry.file_size)
    id = Oid(hex(entry.id))
    return GitIndexEntry(path, id, ctime, mtime, file_size,
                         dev, ino, mode, uid, gid, valid, stage)
end

Base.push!(idx::GitIndex, entry::GitIndexEntry) = begin
    estruct = IndexEntryStruct(entry)
    @check ccall((:git_index_add, libgit2), Cint,
                 (Ptr{Void}, Ptr{IndexEntryStruct}), idx, &estruct)
    return idx
end

Base.push!(idx::GitIndex, path::String) = begin
    @check ccall((:git_index_add_bypath, libgit2), Cint,
                 (Ptr{Void}, Ptr{Uint8}), idx, path)
    return idx
end

Base.getindex(idx::GitIndex, i::Integer) = begin
    i <= 0 && throw(BoundsError())
    entryptr = ccall((:git_index_get_byindex, libgit2), Ptr{IndexEntryStruct},
                     (Ptr{Void}, Csize_t), idx, i-1)
    entryptr == C_NULL && return nothing
    return GitIndexEntry(entryptr)
end

Base.getindex(idx::GitIndex, path::String, stage::Integer=0) = begin
    entryptr = ccall((:git_index_get_bypath, libgit2), Ptr{IndexEntryStruct},
                      (Ptr{Void}, Ptr{Uint8}, Cint), idx, path, stage)
    entryptr == C_NULL && return nothing
    return GitIndexEntry(entryptr)
end

Base.(:(==))(e1::GitIndexEntry, e2::GitIndexEntry) = begin
    return e1.path  == e2.path &&
           e1.id == e2.id &&
           e1.ctime == e2.ctime &&
           e1.mtime == e2.mtime &&
           e1.file_size == e2.file_size &&
           e1.dev  == e2.dev &&
           e1.ino  == e2.ino &&
           e1.mode == e2.mode &&
           e1.uid  == e2.uid &&
           e1.gid  == e2.gid &&
           e1.valid == e2.valid &&
           e1.stage == e2.stage
end
Base.isequal(e1::GitIndexEntry, e2::GitIndexEntry) = e1 == e2

Base.start(idx::GitIndex) = begin
    if length(idx) == 0
        return (0, nothing)
    end
    return (1, idx[1])
end

Base.done(idx::GitIndex, state) = begin
    i = state[1]
    return i == 0 || i > length(idx)
end

Base.next(idx::GitIndex, state) = begin
    i = state[1] + 1
    return (state[2], (i, idx[i]))
end

function read_tree!(idx::GitIndex, tree::GitTree)
    @check ccall((:git_index_read_tree, libgit2), Cint,
                 (Ptr{Void}, Ptr{Void}), idx, tree)
    return idx
end

function add_all!(idx::GitIndex, pathspec::String;
                 force::Bool=false,
                 disable_pathsepc_match::Bool=false,
                 check_pathspec::Bool=false)
    return add_all!(idx, [pathspec],
                    force=force,
                    disable_pathsepc_match=disable_pathsepc_match,
                    check_pathspec=check_pathspec)
end

function add_all!{T<:String}(idx::GitIndex, pathspecs::Vector{T};
                             force::Bool=false,
                             disable_pathsepc_match::Bool=false,
                             check_pathspec::Bool=false)
    flags = GitConst.INDEX_ADD_DEFAULT
    if force
        flags |= GitConst.INDEX_ADD_FORCE
    end
    if disable_pathsepc_match
        flags |= GitConst.INDEX_ADD_DISABLE_PATHSPEC_MATCH
    end
    if check_pathspec
        flags |= GitConst.INDEX_ADD_CHECK_PATHSPEC
    end
    sa_ptr = [StrArrayStruct(pathspecs)]
    ex_ptr = Cint[0]
    @check ccall((:git_index_add_all, libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Cuint, Ptr{Void}, Ptr{Cint}),
                 idx, sa_ptr, flags, C_NULL, ex_ptr)
    ccall((:git_strarray_free, libgit2), Void, (Ptr{StrArrayStruct},), sa_ptr)
    if ex_ptr[1] != GitErrorConst.GIT_OK
        throw(LibGitError(ex_ptr[1]))
    end
    return idx
end

function update_all!{T<:String}(idx::GitIndex, pathspecs::Vector{T})
    sa_ptr = [StrArrayStruct(pathspecs)]
    ex_ptr = Cint[0]
    @check ccall((:git_index_update_all, libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Cint}),
                 idx, sa_ptr, C_NULL, ex_ptr)
    ccall((:git_strarray_free, libgit2), Void, (Ptr{StrArrayStruct},), sa_ptr)
    if ex_ptr[1] != GitErrorConst.GIT_OK
        throw(LibGitError(ex_ptr[1]))
    end
    return idx
end
update_all!(idx::GitIndex, pathspec::String) = update_all!(idx, [pathspec])
update_all!(idx::GitIndex) = update_all!(idx, [""])

function remove_all!{T<:String}(idx::GitIndex, pathspecs::Vector{T})
    sa_ptr = [StrArrayStruct(pathspecs)]
    ex_ptr = Cint[0]
    @check ccall((:git_index_remove_all, libgit2), Cint,
                 (Ptr{Void}, Ptr{StrArrayStruct}, Ptr{Void}, Ptr{Cint}),
                 idx, sa_ptr, C_NULL, ex_ptr)
    ccall((:git_strarray_free, libgit2), Void, (Ptr{StrArrayStruct},), sa_ptr)
    if ex_ptr[1] != GitErrorConst.GIT_OK
        throw(LibGitError(ex_ptr[1]))
    end
    return idx
end
remove_all!(idx::GitIndex, pathspec::String) = remove_all!(idx, [pathspec])

conflicts(r::GitRepo, idx::GitIndex) = begin
    ancestor, ours, theirs = nothing, nothing, nothing
    current_path = nothing
    res = {}
    for entry in idx
        if entry.stage == 0 # staged
            continue
        end
        if current_path != nothing && entry.path != current_path
            push!(res, (ancestor, ours, theirs))
            ancestor, ours, theirs = nothing, nothing, nothing
        end
        current_path = entry.path
        if entry.stage == 1 # ancestor
            ancestor = entry
        elseif  entry.stage == 2 # ours
            ours = entry
        elseif entry.stage == 3 # theirs
            theirs = entry
        else
            error("entry $entry has unexpected stage level $(entry.stage)")
        end
    end
    if current_path != nothing
        push!(res, (ancestor, ours, theirs))
    end
    return res
end


text:
# Color Definitions
# to define or change the RGB based on a 288 scale, divide by 288 for each rgb and set A = 1
import Plots
NUCLEAR = Colors.RGBA(0.615, 0.108, 0.125, 1)
COAL = Colors.RGBA(0.118, 0.118, 0.118, 1)
HYDRO = Colors.RGBA(0.083, 0.441, 0.514, 1)
GAS_CC = Colors.RGBA(0.285, 0.462, 0.372, 1)
GAS_CT = Colors.RGBA(0.493, 0.406, 0.563, 1)
STORAGE = Colors.RGBA(0.128, 0.528, 0.556, 1)
OIL_ST = Colors.RGBA(0.462, 0.212, 0.351, 1) # petroleum
OIL_CT = Colors.RGBA(0.462, 0.212, 0.351, 1) # petroleum
SYNC_COND = Colors.RGBA(0.462, 0.212, 0.351, 1) # petroleum
WIND = Colors.RGBA(0.000, 0.632, 0.830, 1)
SOLAR = Colors.RGBA(0.885, 0.594, 0.007, 1)
CSP = Colors.RGBA(0.875, 0.410, 0.094, 1)
CURTAILMENT = Colors.RGBA(0.847, 0.219, 0.295, 1)

# Out of a 288 rgba scale
NUCLEAR_288 = "rgba(177, 31, 36, 1)"
COAL_288 = "rgba(34, 34, 34, 1)"
HYDRO_288 = "rgba(24, 127, 148, 1)"
GAS_CC_288 = "rgba(82, 133, 107, 1)"
GAS_CT_288 = "rgba(142, 117, 162, 1)"
STORAGE_288 = "rgba(37, 152, 160, 1)"
OIL_ST_288 = "rgba(133, 61, 101, 1)" # petroleum
OIL_CT_288 = "rgba(133, 61, 101, 1)" # petroleum
SYNC_COND_288 = "rgba(133, 61, 101, 1)" # petroleum
WIND_288 = "rgba(0, 182, 239, 1)"
SOLAR_288 = "rgba(255, 171, 2, 1)"
CSP_288 = "rgba(252, 118, 27, 1)"
CURTAILMENT_288 = "rgba(244, 63, 85, 1)"

GR_DEFAULT = hcat(
    NUCLEAR,
    COAL,
    HYDRO,
    GAS_CC,
    GAS_CT,
    STORAGE,
    OIL_ST,
    OIL_CT,
    SYNC_COND,
    WIND,
    SOLAR,
    CSP,
    CURTAILMENT,
)

FUEL_DEFAULT = vcat(
    NUCLEAR,
    COAL,
    HYDRO,
    GAS_CC,
    GAS_CT,
    STORAGE,
    OIL_ST,
    OIL_CT,
    SYNC_COND,
    WIND,
    SOLAR,
    CSP,
    CURTAILMENT,
)

PLOTLY_DEFAULT = vcat(
    NUCLEAR_288,
    COAL_288,
    HYDRO_288,
    GAS_CC_288,
    GAS_CT_288,
    STORAGE_288,
    OIL_ST_288,
    OIL_CT_288,
    SYNC_COND_288,
    WIND_288,
    SOLAR_288,
    CSP_288,
    CURTAILMENT_288,
)

VARIABLE_TYPES = ["P", "Spin", "Reg", "Flex"]

function _filter_variables(results::IS.Results; kwargs...)
    filter_results = Dict()
    reserves = get(kwargs, :reserves, false)
    if reserves
        for (key, var) in results.variables
            start = split("$key", "_")[1]
            if in(start, VARIABLE_TYPES)
                filter_results[key] = var
            end
        end
    else
        for (key, var) in results.variables
            start = split("$key", "_")[1]
            if start == "P"
                filter_results[key] = var
            end
        end
    end
    results = Results(
        filter_results,
        results.total_cost,
        results.optimizer_log,
        results.time_stamp,
    )
    return results
end

function match_fuel_colors(
    stack::StackedGeneration,
    bar::BarGeneration,
    backend::Any,
    default::Array,
)
    if backend == Plots.PlotlyJSBackend()
        color_range = PLOTLY_DEFAULT
    else
        color_range = FUEL_DEFAULT
    end
    fuels = [
        "Nuclear",
        "Coal",
        "Hydro",
        "Gas_CC",
        "Gas_CT",
        "Storage",
        "Oil_ST",
        "Oil_CT",
        "Sync_Cond",
        "Wind",
        "Solar",
        "CSP",
        "curtailment",
    ]
    color_fuel = DataFrames.DataFrame(fuels = fuels, colors = color_range)
    default =
        [(color_fuel[findall(in(["$(bar.labels[1])"]), color_fuel.fuels), :][:, :colors])[1]]
    for i in 2:length(bar.labels)
        specific_color =
            (color_fuel[findall(in(["$(bar.labels[i])"]), color_fuel.fuels), :][
                :,
                :colors,
            ])[1]
        default = hcat(default, specific_color)
    end
    return default
end
"""
    fuel_plot(results, system)

This function makes a stack plot of the results by fuel type
and assigns each fuel type a specific color.

# Arguments

- `res::Results = results`: results to be plotted
- `system::PSY.System`: The power systems system

# Example

```julia
res = solve_op_problem!(OpProblem)
fuel_plot(res, sys)
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""
function fuel_plot(res::IS.Results, sys::PSY.System; kwargs...)
    ref = make_fuel_dictionary(sys, res)
    fuel_plot(res, ref; kwargs...)
end

function fuel_plot(res::Array{IS.Results}, sys::PSY.System; kwargs...)
    ref = make_fuel_dictionary(sys, res)
    fuel_plot(res, ref; kwargs...)
end
"""
    fuel_plot(results::IS.Results, generators)

This function makes a stack plot of the results by fuel type
and assigns each fuel type a specific color.

# Arguments

- `res::IS.Results = results`: results to be plotted
- `generators::Dict`: the dictionary of fuel type and an array
 of the generators per fuel type, or some other specified category order

# Example

```julia
res = solve_op_problem!(OpProblem)
generator_dict = make_fuel_dictionary(sys, res)
fuel_plot(res, generator_dict)
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function fuel_plot(res::IS.Results, generator_dict::Dict; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    res = _filter_variables(res; kwargs...)
    stack = get_stacked_aggregation_data(res, generator_dict)
    bar = get_bar_aggregation_data(res, generator_dict)
    backend = Plots.backend()
    default_colors = match_fuel_colors(stack, bar, backend, FUEL_DEFAULT)
    seriescolor = get(kwargs, :seriescolor, default_colors)
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _fuel_plot_internal(stack, bar, seriescolor, backend, save_fig, set_display; kwargs...)
end

function fuel_plot(results::Array{}, generator_dict::Dict; kwargs...)
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    new_res = _filter_variables(results[1]; kwargs...)
    stack = get_stacked_aggregation_data(new_res, generator_dict)
    bar = get_bar_aggregation_data(new_res, generator_dict)
    for i in 2:length(results)
        new_res = _filter_variables(results[i]; kwargs...)
        stack = hcat(stack, get_stacked_aggregation_data(new_res, generator_dict))
        bar = hcat(bar, get_bar_aggregation_data(new_res, generator_dict))
    end
    backend = Plots.backend()
    default_colors = match_fuel_colors(stack[1], bar[1], backend, FUEL_DEFAULT)
    seriescolor = get(kwargs, :seriescolor, default_colors)
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _fuel_plot_internal(stack, bar, seriescolor, backend, save_fig, set_display; kwargs...)
end

function _fuel_plot_internal(
    stack::Any,
    bar::Any,
    seriescolor::Array,
    backend::Plots.PlotlyJSBackend,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    plotly_stack_gen(stack, seriescolor; kwargs...)
    plotly_bar_gen(bar, seriescolor; kwargs...)
end

function _fuel_plot_internal(
    stack::Any,
    bar::Any,
    seriescolor::Array,
    backend::Any,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    p1 = RecipesBase.plot(stack, seriescolor)
    p2 = RecipesBase.plot(bar, seriescolor)
    set_display && display(p1)
    set_display && display(p2)
    if !isnothing(save_fig)
        Plots.savefig(p1, joinpath(save_fig, "Fuel_Stack.png"))
        Plots.savefig(p2, joinpath(save_fig, "Fuel_Bar.png"))
    end
end

"""
   bar_plot(results::IS.Results)

This function plots a bar plot for the generators in each variable within
the results variables dictionary, and makes a bar plot for all of the variables.

# Arguments
- `res::IS.Results = results`: results to be plotted

# Example

```julia
results = solve_op_problem!(OpProblem)
bar_plot(results)
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function bar_plot(res::IS.Results; kwargs...)
    backend = Plots.backend()
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    res = _filter_variables(res; kwargs...)
    bar_gen = get_bar_gen_data(res)
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _bar_plot_internal(res, bar_gen, backend, save_fig, set_display; kwargs...)
end

"""
   bar_plot(results::Array{IS.Results})

This function plots a subplot for each result. Each subplot has a bar plot for the generators in each variable within
the results variables dictionary, and makes a bar plot for all of the variables per result object.

# Arguments
- `res::Array{IS.Results} = [results1; results2]`: results to be plotted

# Example

```julia
results1 = solve_op_problem!(OpProblem1)
results2 = solve_op_problem!(OpProblem2)
bar_plot([results1; results2])
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function bar_plot(results::Array{}; kwargs...)
    backend = Plots.backend()
    set_display = get(kwargs, :display, true)
    save_fig = get(kwargs, :save, nothing)
    res = _filter_variables(results[1]; kwargs...)
    bar_gen = get_bar_gen_data(res)
    for i in 2:size(results, 1)
        filter = _filter_variables(results[i]; kwargs...)
        res = hcat(res, filter)
        bar_gen = hcat(bar_gen, get_bar_gen_data(filter))
    end
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _bar_plot_internal(res, bar_gen, backend, save_fig, set_display; kwargs...)
end

function bar_plot(res::IS.Results, variables::Array; kwargs...)
    res_var = Dict()
    for variable in variables
        res_var[variable] = res.variables[variable]
    end
    results = Results(res_var, res.total_cost, res.optimizer_log, res.time_stamp)
    bar_plot(results; kwargs...)
end

function bar_plot(results::Array{IS.Results}, variables::Array; kwargs...)
    new_results = []
    for res in results
        res_var = Dict()
        for variable in variables
            res_var[variable] = res.variables[variable]
        end
        new_res = Results(res_var, res.total_cost, res.optimizer_log, res.time_stamp)
        new_results = vcat(new_results, new_res)
    end
    bar_plot(new_results; kwargs...)
end

function _bar_plot_internal(
    res::IS.Results,
    bar_gen::BarGeneration,
    backend::Plots.PlotlyJSBackend,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, PLOTLY_DEFAULT)
    plotly_bar_plots(res, seriescolor; kwargs...)
    plotly_bar_gen(bar_gen, seriescolor; kwargs...)
end

function _bar_plot_internal(
    res::Array,
    bar_gen::Array,
    backend::Plots.PlotlyJSBackend,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, PLOTLY_DEFAULT)
    plotly_bar_plots(res, seriescolor; kwargs...)
    plotly_bar_gen(bar_gen, seriescolor; kwargs...)
end

function _bar_plot_internal(
    res::IS.Results,
    bar_gen::BarGeneration,
    backend::Any,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, GR_DEFAULT)
    for name in string.(keys(res.variables))
        variable_bar = get_bar_plot_data(res, name)
        p = RecipesBase.plot(variable_bar, name, seriescolor)
        set_display && display(p)
        if !isnothing(save_fig)
            Plots.savefig(p, joinpath(save_fig, "$(name)_Bar.png"))
        end
    end
    p2 = RecipesBase.plot(bar_gen, seriescolor)
    set_display && display(p2)
    if !isnothing(save_fig)
        Plots.savefig(p2, joinpath(save_fig, "Bar_Generation.png"))
    end
end

function _bar_plot_internal(
    results::Array,
    bar_gen::Array,
    backend::Any,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, GR_DEFAULT)
    for name in string.(keys(results[1, 1].variables))
        variable_bar = get_bar_plot_data(results[1, 1], name)
        for i in 2:length(results)
            variable_bar = hcat(variable_bar, get_bar_plot_data(results[i], name))
        end
        p = RecipesBase.plot(variable_bar, name, seriescolor)
        set_display && display(p)
        if !isnothing(save_fig)
            Plots.savefig(p, joinpath(save_fig, "$(name)_Bar.png"))
        end
    end
    p2 = RecipesBase.plot(bar_gen, seriescolor)
    set_display && display(p2)
    if !isnothing(save_fig)
        Plots.savefig(p2, joinpath(save_fig, "Bar_Generation.png"))
    end
end

"""
     stack_plot(results::IS.Results)

This function plots a stack plot for the generators in each variable within
the results variables dictionary, and makes a stack plot for all of the variables.

# Arguments
- `res::IS.Results = results`: results to be plotted

# Examples

```julia
results = solve_op_problem!(OpProblem)
stack_plot(results)
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function stack_plot(res::IS.Results; kwargs...)
    set_display = get(kwargs, :set_display, true)
    backend = Plots.backend()
    save_fig = get(kwargs, :save, nothing)
    _filter_variables(res; kwargs...)
    stacked_gen = get_stacked_generation_data(res)
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _stack_plot_internal(res, stacked_gen, backend, save_fig, set_display; kwargs...)
end

"""
     stack_plot(results::Array{IS.Results})

This function plots a subplot for each result object. Each subplot stacks the generators in each variable within
results variables dictionary, and makes a stack plot for all of the variables per result object.

# Arguments
- `res::Array{IS.Results} = [results1, results2]`: results to be plotted

# Examples

```julia
results1 = solve_op_problem!(OpProblem1)
results2 = solve_op_problem!(OpProblem2)
stack_plot([results1; results2])
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function stack_plot(results::Array{}; kwargs...)
    set_display = get(kwargs, :set_display, true)
    backend = Plots.backend()
    save_fig = get(kwargs, :save, nothing)
    new_results = _filter_variables(results[1]; kwargs...)
    stacked_gen = get_stacked_generation_data(new_results)
    for res in 2:length(results)
        filtered = _filter_variables(results[res]; kwargs...)
        new_results = hcat(new_results, filtered)
        stacked_gen = hcat(stacked_gen, get_stacked_generation_data(filtered))
    end
    if isnothing(backend)
        throw(IS.ConflictingInputsError("No backend detected. Type gr() to set a backend."))
    end
    _stack_plot_internal(
        new_results,
        stacked_gen,
        backend,
        save_fig,
        set_display;
        kwargs...,
    )
end

"""
     stack_plot(results::IS.Results, variables::Array)

This function plots a stack plot for the generators in each variable within
the results variables dictionary, and makes a stack plot for all of the variables in the array.

# Arguments
- `res::IS.Results = results`: results to be plotted
- `variables::Array`: list of variables to be plotted in the results

# Examples

```julia
results = solve_op_problem!(OpProblem)
variables = [:var1, :var2, :var3]
stack_plot(results, variables)
```

# Accepted Key Words
- `display::Bool`: set to false to prevent the plots from displaying
- `save::String = "file_path"`: set a file path to save the plots
- `seriescolor::Array`: Set different colors for the plots
- `reserves::Bool`: if reserves = true, the researves will be plotted with the active power
"""

function stack_plot(res::IS.Results, variables::Array; kwargs...)
    res_var = Dict()
    for variable in variables
        res_var[variable] = res.variables[variable]
    end
    results = Results(res_var, res.total_cost, res.optimizer_log, res.time_stamp)
    stack_plot(results; kwargs...)
end

function stack_plot(results::Array{IS.Results}, variables::Array; kwargs...)
    new_results = []
    for res in results
        res_var = Dict()
        for variable in variables
            res_var[variable] = res.variables[variable]
        end
        new_res = Results(res_var, res.total_cost, res.optimizer_log, res.time_stamp)
        new_results = vcat(new_results, new_res)
    end
    stack_plot(new_results; kwargs...)
end

function _stack_plot_internal(
    res::IS.Results,
    stack::StackedGeneration,
    backend::Plots.PlotlyJSBackend,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, PLOTLY_DEFAULT)
    plotly_stack_plots(res, seriescolor; kwargs...)
    plotly_stack_gen(stack, seriescolor; kwargs...)
end

function _stack_plot_internal(
    res::Array{},
    stack::Array{},
    backend::Plots.PlotlyJSBackend,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, PLOTLY_DEFAULT)
    plotly_stack_plots(res, seriescolor; kwargs...)
    plotly_stack_gen(stack, seriescolor; kwargs...)
end

function _stack_plot_internal(
    res::IS.Results,
    stack::StackedGeneration,
    backend::Any,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, GR_DEFAULT)
    for name in string.(keys(res.variables))
        variable_stack = get_stacked_plot_data(res, name)
        p = RecipesBase.plot(variable_stack, name, seriescolor)
        set_display && display(p)
        if !isnothing(save_fig)
            Plots.savefig(p, joinpath(save_fig, "$(name)_Stack.png"))
        end
    end
    p2 = RecipesBase.plot(stack, seriescolor)
    set_display && display(p2)
    if !isnothing(save_fig)
        Plots.savefig(p2, joinpath(save_fig, "Stack_Generation.png"))
    end
end

function _stack_plot_internal(
    results::Any,
    stack::Any,
    backend::Any,
    save_fig::Any,
    set_display::Bool;
    kwargs...,
)
    seriescolor = get(kwargs, :seriescolor, GR_DEFAULT)

    for name in string.(keys(results[1].variables))
        variable_stack = get_stacked_plot_data(results[1], name)
        for res in 2:length(results)
            variable_stack = hcat(variable_stack, get_stacked_plot_data(results[res], name))
        end
        p = RecipesBase.plot(variable_stack, name, seriescolor)
        set_display && display(p)
        if !isnothing(save_fig)
            Plots.savefig(p, joinpath(save_fig, "$(name)_Stack.png"))
        end
    end
    p2 = RecipesBase.plot(stack, seriescolor)
    set_display && display(p2)
    if !isnothing(save_fig)
        Plots.savefig(p2, joinpath(save_fig, "Stack_Generation.png"))
    end
end


text:
@testset "Index: $IndexType" for IndexType in [NaiveIndex, BruteTreeIndex,
                                    KDTreeIndex, HNSWIndex, IVFIndex, NoopIndex]
    T = eltype(1.0)
    data = T[0 0 0 5 5 5; 0 1 2 10 11 12]
    spdata = sparse(data)
    point = T[5.1, 10]
    true_length = size(data, 2)

    if IndexType === IVFIndex
        _indexfunc = d->IVFIndex(d; kc=4, k=2, m=1)
        index = _indexfunc(data)
        spindex = _indexfunc(spdata)
    else
        index = IndexType(data)
        spindex = IndexType(data)
    end
    @test index isa IndexType
    idxs, scores = Garamond.knn_search(index, point, 10; w=4)
    @test idxs isa Vector{Int}
    @test scores isa Vector{T}
    if IndexType != NoopIndex
        @test all(i in idxs for i in 1:true_length)
    else
        @test isempty(idxs)
        @test isempty(scores)
    end
    @test length(index) == length(spindex) == true_length

    # Test push!, pop! interface
    idx = 1
    point = data[:, idx]
    if IndexType in [NaiveIndex, BruteTreeIndex, IVFIndex]
        # push!
        @test try push!(index, point); true; catch; false; end
        @test length(index) == true_length + 1

        # pop!
        if !isa(index, IVFIndex)
            @test pop!(index) == point
        else
            @test pop!(index) isa Vector{T} # ivfadc quantizes the vector, it will not be the same
        end
        @test length(index) == true_length

        # pushfirst!
        @test try pushfirst!(index, point); true; catch; false; end
        @test length(index) == true_length + 1

        # popfirst!
        if !isa(index, IVFIndex)
            @test popfirst!(index) == point
        else
            @test popfirst!(index) isa Vector{T}
        end
        @test length(index) == true_length

        # delete_from_index!
        @test Garamond.delete_from_index!(index, [idx]) === nothing
        @test length(index) == true_length - 1
    elseif IndexType == NoopIndex
        # push!
        @test try push!(index, point); true; catch; false; end
        @test length(index) == true_length  # nothing pushed

        # pop!
        @test pop!(index) === nothing
        @test length(index) == true_length

        # pushfirst!
        @test try pushfirst!(index, point); true; catch; false; end
        @test length(index) == true_length  # nothing pushed

        # popfirst!
        @test popfirst!(index) === nothing
        @test length(index) == true_length

        # delete_from_index!
        @test Garamond.delete_from_index!(index, [idx]) === nothing
        @test length(index) == true_length

    elseif IndexType in [KDTreeIndex, HNSWIndex]
        @test_throws MethodError push!(idx, point)
        @test_throws MethodError pushfirst!(idx, point)
        @test_throws MethodError pop!(idx)
        @test_throws MethodError popfirst!(idx)
        @test_throws MethodError Garamond.delete_from_index!(idx, [1,2])
    end
end


text:
"""
Example of how to obtain the bid which maximises the profit of the Generator at the slack
bus using the Package MaximumGeneratorProfitTCRDD. The Algorithm is based on the
Transmission Constrained Residual Demand Derivative (TCRDD) [1].

[1] L. Xu, R. Baldick and Y. Sutjandra, "Bidding Into Electricity Markets:
    A Transmission-Constrained Residual Demand Derivative Approach," in IEEE
    Transactions on Power Systems, vol. 26, no. 3, pp. 1380-1388, Aug. 2011,
    doi: 10.1109/TPWRS.2010.2083702.
"""
# Packages
using DataFrames
using Dates
using DataStructures
using Ipopt
using MaximumGeneratorProfitTCRDD
using PowerSystems
using PowerSimulations
using Test
using TimeSeries

# Build Test Case
include("data/example_case4/buildup_case4.jl")
include("data/example_case4/c_sys_case4.jl")
sys = c_sys_case4()

# Settings
# Solver for OPF
solver = optimizer_with_attributes(Ipopt.Optimizer)
# Network representation
network = StandardPTDFModel
# Tolerance to identify binding constraints
dual_lines_tol = 1e-1
dual_gen_tol = 1e-1
# Segments to evaluate approximated profit function
segm_bid_argmax_profit = 50000
# Maximum iterations for Screening and Bisection Loops
maxit_scr = 3
maxit_bi = 5
# Convergence Tolerance between bids in [pu]
epsilon = 0.01 #tolerance bid_lo - bid_hi < epsilon
# Flags
print_results = true
print_progress = true
print_plots = true

# Initial Bid
bid0 = 0.3

# Maximize Generator Profit
bid_opt_found, bid_opt, bid_lo, bid_hi, iter_scr, iter_bi = maxGenProfit_tcrdd(
    sys,
    bid0;
    dual_lines_tol,
    dual_gen_tol,
    segm_bid_argmax_profit,
    maxit_scr,
    maxit_bi,
    epsilon,
    print_results,
    print_progress,
    print_plots,
    network,
    solver
)

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
using VDPTag
using POMDPs
using Base.Test
using POMDPToolbox
# using Plots
using MCTS

# include("discretization.jl")

pomdp = VDPTagPOMDP()
gen = NextMLFirst(mdp(pomdp), MersenneTwister(31))
s = TagState([1.0, 1.0], [-1.0, -1.0])

struct MyNode end
MCTS.n_children(::MyNode) = rand(1:10)

@inferred next_action(gen, pomdp, s, MyNode())
@inferred next_action(gen, pomdp, initial_state_distribution(pomdp), MyNode())

dpomdp = ADiscreteVDPTagPOMDP(pomdp, 10)
@show actions(dpomdp)
for sao in stepthrough(dpomdp, RandomPolicy(dpomdp), "sao", max_steps=10)
    @show sao
end


text:
############################################################
# joBlockDiag ##############################################
############################################################

export joBlockDiag, joBlockDiagException

struct joBlockDiagException <: Exception
    msg :: String
end

############################################################
## outer constructors

"""
    julia> op = joBlockDiag(op1[,op2][,...];[weights=...,][name=...])

Block-diagonal operator composed from different square JOLI operators

# Signature

    joBlockDiag(ops::joAbstractLinearOperator...;
        weights::LocalVector{WDT}=zeros(0),name::String="joBlockDiag")
            where {WDT<:Number}

# Arguments

- `op#`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `weights`: vector of waights for each operator
    - `name`: custom name

# Notes

- all operators must be square (M(i)==N(i))
- all given operators must have same domain/range types
- the domain/range types of joBlockDiag are equal to domain/range types of the given operators

# Example

define operators

    a=rand(ComplexF64,4,4);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")
    b=rand(ComplexF64,8,8);
    B=joMatrix(b;DDT=ComplexF32,RDT=ComplexF64,name="B")
    c=rand(ComplexF64,6,6);
    C=joMatrix(c;DDT=ComplexF32,RDT=ComplexF64,name="C")

define weights if needed

    w=rand(ComplexF64,3)

basic block diagonal

    BD=joBlockDiag(A,B,C)

weighted block diagonal

    BD=joBlockDiag(A,B,C;weights=w)

"""
function joBlockDiag(ops::joAbstractLinearOperator...;
           weights::LocalVector{WDT}=zeros(0),name::String="joBlockDiag") where {WDT<:Number}

    isempty(ops) && throw(joBlockDiagException("empty argument list"))
    l=length(ops)
    for i=1:l
        ops[i].m==ops[i].n || throw(joBlockDiagException("non-square $i operator"))
        deltype(ops[i])==deltype(ops[1]) || throw(joBlockDiagException("domain type mismatch for $i operator"))
        reltype(ops[i])==reltype(ops[1]) || throw(joBlockDiagException("range type mismatch for $i operator"))
    end
    (length(weights)==l || length(weights)==0) || throw(joBlockDiagException("lenght of weights vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    mo=zeros(Int,l)
    no=zeros(Int,l)
    for i=1:l
        ms[i]=ops[i].m
        ns[i]=ops[i].n
    end
    for i=2:l
        mo[i]=mo[i-1]+ms[i-1]
        no[i]=no[i-1]+ns[i-1]
    end
    m=sum(ms)
    n=sum(ns)
    weighted=(length(ws)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*ops[i])
            push!(fops_T,ws[i]*transpose(ops[i]))
            push!(fops_A,conj(ws[i])*adjoint(ops[i]))
            push!(fops_C,conj(ws[i])*conj(ops[i]))
        else
            push!(fops,ops[i])
            push!(fops_T,transpose(ops[i]))
            push!(fops_A,adjoint(ops[i]))
            push!(fops_C,conj(ops[i]))
        end
    end
    return joCoreBlock{deltype(fops[1]),reltype(fops[1])}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end
"""
    julia> op = joBlockDiag(l,op;[weights=...,][name=...])

Block-diagonal operator composed from l-times replicated square JOLI operator

# Signature

    joBlockDiag(l::Integer,op::joAbstractLinearOperator;
        weights::LocalVector{WDT}=zeros(0),name::String="joBlockDiag")
            where {WDT<:Number}

# Arguments

- `l`: # of replcated blocks
- `op`: JOLI operators (subtypes of joAbstractLinearOperator)
- keywords
    - `weights`: vector of waights for each operator
    - `name`: custom name

# Notes

- all given operators must have same domain/range types
- the domain/range types of joBlockDiag are equal to domain/range types of the given operators

# Example

define operator

    a=rand(ComplexF64,4,4);
    A=joMatrix(a;DDT=ComplexF32,RDT=ComplexF64,name="A")

define weights if needed

    w=rand(ComplexF64,3)

basic block diagonal

    BD=joBlockDiag(3,A)

weighted block diagonal

    BD=joBlockDiag(3,A;weights=w)

"""
function joBlockDiag(l::Integer,op::joAbstractLinearOperator;
           weights::LocalVector{WDT}=zeros(0),name::String="joBlockDiag") where {WDT<:Number}

    op.m==op.n || throw(joBlockDiagException("non-square operator"))
    (length(weights)==l || length(weights)==0) || throw(joBlockDiagException("lenght of weights vector does not match number of operators"))
    ws=Base.deepcopy(weights)
    ms=zeros(Int,l)
    ns=zeros(Int,l)
    mo=zeros(Int,l)
    no=zeros(Int,l)
    for i=1:l
        ms[i]=op.m
        ns[i]=op.n
    end
    for i=2:l
        mo[i]=mo[i-1]+ms[i-1]
        no[i]=no[i-1]+ns[i-1]
    end
    m=l*op.m
    n=l*op.n
    weighted=(length(weights)==l)
    fops=Vector{joAbstractLinearOperator}(undef,0)
    fops_T=Vector{joAbstractLinearOperator}(undef,0)
    fops_A=Vector{joAbstractLinearOperator}(undef,0)
    fops_C=Vector{joAbstractLinearOperator}(undef,0)
    for i=1:l
        if weighted
            push!(fops,ws[i]*op)
            push!(fops_T,ws[i]*transpose(op))
            push!(fops_A,conj(ws[i])*adjoint(op))
            push!(fops_C,conj(ws[i])*conj(op))
        else
            push!(fops,op)
            push!(fops_T,transpose(op))
            push!(fops_A,adjoint(op))
            push!(fops_C,conj(op))
        end
    end
    return joCoreBlock{deltype(op),reltype(op)}(name*"($l)",m,n,l,ms,ns,mo,no,ws,
                      fops,fops_T,fops_A,fops_C,@joNF,@joNF,@joNF,@joNF)
end


text:
using BAT
using ValueShapes
using ParallelProcessingTools
using DensityInterface

export get_prior, get_likelihood
export plot_model_space, plot_data_space

"""
    get_prior(pdf_params)

Get a suitable prior for a given PDF parametrisation.
"""
function get_prior end

function get_prior(pdf_params::ValencePDFParams)

    prior = NamedTupleDist(
        θ_tmp = Dirichlet(pdf_params.weights),
        λ_u = Uniform(0, 1),
        K_u = Uniform(2, 10),
        λ_d = Uniform(0, 1),
        K_d = Uniform(2, 10),
        λ_g1 = Uniform(0, 1),
        λ_g2 = Uniform(-1, 0),
        K_g =  Uniform(2, 10),
        λ_q = Uniform(-1, 0),
    )

    return prior
end

function get_prior(pdf_params::DirichletPDFParams)

    prior = NamedTupleDist(
        θ = Dirichlet(pdf_params.weights),
        K_u = Uniform(2, 10),
        K_d = Uniform(2, 10),
        λ_g1 = Uniform(0, 1),
        λ_g2 = Uniform(-1, 0),
        K_g =  Uniform(2, 10),
        λ_q = Uniform(-1, 0),
    )

    return prior
end

"""
    get_likelihood(pdf_params)

Get a suitable likelihood for a given PDF parametrisation.
"""
function get_likelihood end

function get_likelihood(pdf_params::ValencePDFParams, sim_data::Dict{String, Any},
                        qcdnum_params::QCDNUMParameters, splint_params::SPLINTParameters,
                        quark_coeffs::QuarkCoefficients)

    likelihood = let d = sim_data

        counts_obs_ep = Int.(d["counts_obs_ep"])
        counts_obs_em = Int.(d["counts_obs_em"])
        nbins = d["nbins"]

        logfuncdensity(function (params)

            θ = get_scaled_θ(params.λ_u, params.K_u, params.λ_d,
                             params.K_d, Vector(params.θ_tmp))

            pdf_params = ValencePDFParams(λ_u=params.λ_u, K_u=params.K_u, λ_d=params.λ_d,
                                          K_d=params.K_d, λ_g1=params.λ_g1, λ_g2=params.λ_g2,
                                          K_g=params.K_g, λ_q=params.λ_q, θ=θ)

            counts_pred_ep, counts_pred_em = @critical forward_model(pdf_params, qcdnum_params,
                                                                     splint_params, quark_coeffs)

            ll_value = 0.0
            for i in 1:nbins

                if counts_pred_ep[i] < 0
                    @debug "counts_pred_ep[i] < 0, setting to 0" i counts_pred_ep[i]
                    counts_pred_ep[i] = 0
                end

                if counts_pred_em[i] < 0
                    @debug "counts_pred_em[i] < 0, setting to 0" i counts_pred_em[i]
                    counts_pred_em[i] = 0
                end

                ll_value += logpdf(Poisson(counts_pred_ep[i]), counts_obs_ep[i])
                ll_value += logpdf(Poisson(counts_pred_em[i]), counts_obs_em[i])
            end

            return ll_value
        end)

    end

    return likelihood
end

function get_likelihood(pdf_params::DirichletPDFParams, sim_data::Dict{String, Any},
                        qcdnum_params::QCDNUMParameters, splint_params::SPLINTParameters,
                        quark_coeffs::QuarkCoefficients)

    likelihood = let d = sim_data

        counts_obs_ep = Int.(d["counts_obs_ep"])
        counts_obs_em = Int.(d["counts_obs_em"])
        nbins = d["nbins"]

        logfuncdensity(function (params)

            pdf_params = DirichletPDFParams(K_u=params.K_u, K_d=params.K_d,
                                            λ_g1=params.λ_g1, λ_g2=params.λ_g2,
                                            K_g=params.K_g, λ_q=params.λ_q,
                                            θ=Vector(params.θ))

            counts_pred_ep, counts_pred_em = @critical forward_model(pdf_params, qcdnum_params,
                                                                     splint_params, quark_coeffs)

            ll_value = 0.0
            for i in 1:nbins

                if counts_pred_ep[i] < 0
                    @debug "counts_pred_ep[i] < 0, setting to 0" i counts_pred_ep[i]
                    counts_pred_ep[i] = 0
                end

                if counts_pred_em[i] < 0
                    @debug "counts_pred_em[i] < 0, setting to 0" i counts_pred_em[i]
                    counts_pred_em[i] = 0
                end

                ll_value += logpdf(Poisson(counts_pred_ep[i]), counts_obs_ep[i])
                ll_value += logpdf(Poisson(counts_pred_em[i]), counts_obs_em[i])
            end

            return ll_value
        end)

    end

    return likelihood
end

"""
    plot_model_space(pdf_params, samples)

Compare truth and posterior samples in the model space.
"""
function plot_model_space end

function plot_model_space(pdf_params::AbstractPDFParams, samples; xmin::Float64=1e-3,
                          xmax::Float64=1.0, nx::Integer=50, nsamples::Integer=200,
                          truth_color=:black, sample_color=:skyblue3)

    x_grid = range(xmin, stop=xmax, length=nx)

    p = plot(x_grid, [xtotx(x, pdf_params) for x in x_grid], color=truth_color,
             lw=3, label="Truth")

    sub_samples = bat_sample(samples, OrderedResampling(nsamples=nsamples)).result

    p = plot_model_space_impl(x_grid, pdf_params, sub_samples, p, color=sample_color)

    p = plot!(xaxis=:log, yaxis=:log, xlabel="x", ylabel="xtotx")
    p = ylims!(1e-5, 50.0)

    return p
end

function plot_model_space_impl(x_grid::StepRangeLen{Float64}, pdf_params::ValencePDFParams, samples, p; color=:skyblue3)

    for i in eachindex(samples)

        θ_i = get_scaled_θ(samples.v.λ_u[i], samples.v.K_u[i], samples.v.λ_d[i],
                           samples.v.K_d[i], Vector(samples.v.θ_tmp[i]))

        pdf_params_i = ValencePDFParams(λ_u=samples.v.λ_u[i], K_u=samples.v.K_u[i],
                                        λ_d=samples.v.λ_d[i], K_d=samples.v.K_d[i],
                                        λ_g1=samples.v.λ_g1[i], λ_g2=samples.v.λ_g2[i],
                                        K_g=samples.v.K_g[i], λ_q=samples.v.λ_q[i],
                                        θ=θ_i)
        p = plot!(x_grid, [xtotx(x, pdf_params_i) for x in x_grid], color=color, lw=3,
                  alpha=0.01, label="")

    end

    return p
end

function plot_model_space_impl(x_grid::StepRangeLen{Float64}, pdf_params::DirichletPDFParams, samples, p; color=:skyblue3)

    for i in eachindex(sub_samples)

        pdf_params_i = DirichletPDFParams(λ_u=sub_samples.v.λ_u[i], K_u=sub_samples.v.K_u[i],
                                          λ_d=sub_samples.v.λ_d[i], K_d=sub_samples.v.K_d[i],
                                          λ_g1=sub_samples.v.λ_g1[i], λ_g2=sub_samples.v.λ_g2[i],
                                          K_g=sub_samples.v.K_g[i], λ_q=sub_samples.v.λ_q[i],
                                          θ=Vector(sub_samples.v.θ[i]))
        p = plot!(x_grid, [xtotx(x, pdf_params_i) for x in x_grid], color=color, lw=3,
                  alpha=0.01, label="")

    end

    return p
end

"""
    plot_data_space(pdf_params, sim_data, samples, qcdnum_grid,
                    qcdnum_params, splint_params, quark_coeffs)

Compare truth and posterior samples in the data space.
"""
function plot_data_space end

function plot_data_space(pdf_params::AbstractPDFParams, sim_data::Dict{String, Any}, samples,
                         qcdnum_grid::QCDNUMGrid, qcdnum_params::QCDNUMParameters,
                         splint_params::SPLINTParameters, quark_coeffs::QuarkCoefficients;
                         ep_color=:firebrick, em_color=:teal, nsamples::Integer=100, plot_size=(1000, 500))

    forward_model_init(qcdnum_grid, qcdnum_params, splint_params)

    nbins = sim_data["nbins"]

    p1 = scatter(1:nbins, sim_data["counts_obs_ep"], label="Observed counts (eP)",
                 color=ep_color, lw=3)
    p2 = scatter(1:nbins, sim_data["counts_obs_em"], label="Observed counts (eM)",
                 color=em_color, lw=3)

    sub_samples = bat_sample(samples, OrderedResampling(nsamples=nsamples)).result

    p1, p2 = plot_data_space_impl(pdf_params, sub_samples, qcdnum_params,
                                  splint_params, quark_coeffs, p1, p2, nbins)

    plot(p1, p2, size=plot_size, xlabel="Bin number")
end

function plot_data_space_impl(pdf_params::ValencePDFParams, samples, qcdnum_params::QCDNUMParameters,
                              splint_params::SPLINTParameters, quark_coeffs::QuarkCoefficients,
                              p1, p2, nbins::Integer; ep_color=:firebrick, em_color=:teal)

    for i in eachindex(samples)

        counts_obs_ep_i = zeros(UInt64, nbins)
        counts_obs_em_i = zeros(UInt64, nbins)

        θ_i = get_scaled_θ(samples.v.λ_u[i], samples.v.K_u[i], samples.v.λ_d[i],
                           samples.v.K_d[i], Vector(samples.v.θ_tmp[i]))

        pdf_params_i = ValencePDFParams(λ_u=samples.v.λ_u[i], K_u=samples.v.K_u[i],
                                        λ_d=samples.v.λ_d[i], K_d=samples.v.K_d[i],
                                        λ_g1=samples.v.λ_g1[i], λ_g2=samples.v.λ_g2[i],
                                        K_g=samples.v.K_g[i], λ_q=samples.v.λ_q[i],
                                        θ=θ_i)

        counts_pred_ep_i, counts_pred_em_i = forward_model(pdf_params_i, qcdnum_params,
                                                           splint_params, quark_coeffs)

        for j in 1:nbins

            if counts_pred_ep_i[j] < 0

                @warn "Predicted counts (eP) found to be < 0, setting to 0" i j counts_pred_ep_i[j]
                counts_pred_ep_i[j] = 0

            end

            if counts_pred_em_i[j] < 0

                @warn "Predicted counts (eM) found to be < 0, setting to 0" i j counts_pred_em_i[j]
                counts_pred_em_i[j] = 0

            end

            counts_obs_ep_i[j] = rand(Poisson(counts_pred_ep_i[j]))
            counts_obs_em_i[j] = rand(Poisson(counts_pred_em_i[j]))

        end

        p1 = scatter!(p1, 1:nbins, counts_obs_ep_i, label="", color=ep_color,
                      lw=3, alpha=0.01)
        p2 = scatter!(p2, 1:nbins, counts_obs_em_i, label="", color=em_color,
                      lw=3, alpha=0.01)

    end

    return p1, p2
end

function plot_data_space_impl(pdf_params::DirichletPDFParams, samples, qcdnum_params::QCDNUMParameters,
                              splint_params::SPLINTParameters, quark_coeffs::QuarkCoefficients,
                              p1, p2; ep_color=:firebrick, em_color=:teal)

    for i in eachindex(samples)

        counts_obs_ep_i = zeros(UInt64, nbins)
        counts_obs_em_i = zeros(UInt64, nbins)

        pdf_params_i = DirichletPDFParams(λ_u=samples.v.λ_u[i], K_u=samples.v.K_u[i],
                                          λ_d=samples.v.λ_d[i], K_d=samples.v.K_d[i],
                                          λ_g1=samples.v.λ_g1[i], λ_g2=samples.v.λ_g2[i],
                                          K_g=samples.v.K_g[i], λ_q=samples.v.λ_q[i],
                                          θ=Vector(samples.v.θ))

        counts_pred_ep_i, counts_pred_em_i = forward_model(pdf_params_i, qcdnum_params,
                                                           splint_params, quark_coeffs)

        for j in 1:nbins

            if counts_pred_ep_i[j] < 0

                @warn "Predicted counts (eP) found to be < 0, setting to 0" i j counts_pred_ep_i[j]
                counts_pred_ep_i[j] = 0

            end

            if counts_pred_em_i[j] < 0

                @warn "Predicted counts (eM) found to be < 0, setting to 0" i j counts_pred_em_i[j]
                counts_pred_em_i[j] = 0

            end

            counts_obs_ep_i[j] = rand(Poisson(counts_pred_ep_i[j]))
            counts_obs_em_i[j] = rand(Poisson(counts_pred_em_i[j]))

        end

        p1 = scatter!(p1, 1:nbins, counts_obs_ep_i, label="", color=ep_color,
                      lw=3, alpha=0.01)
        p2 = scatter!(p2, 1:nbins, counts_obs_em_i, label="", color=em_color,
                      lw=3, alpha=0.01)

    end

    return p1, p2
end


text:
module Densratio

export
    densratio,
    compute_density_ratio,
    search_sigma_and_lambda,
    compute_kernel_Gaussian,
    kernel_Gaussian

using PyCall
const densratio_py = PyNULL()

function __init__()
    copy!(densratio_py, pyimport_conda("densratio", "densratio"))
end

function densratio(x, y)
    densratio.densratio(x, y)
end

end # module


text:
# This file includes an example for two dimensional interpolation over triangular interpolation domain.

using FractalTools
using Triangulation
using PyPlot
using Statistics

close("all")

datapath = joinpath(@__DIR__,  "data/gmsh_data")
filename = "t3.msh"

# Read mesh file.
mesh = read_mesh(joinpath(datapath, filename))
scale_factor = 1 / 10  # Rescale the mesh.
for node in mesh.nodes
    node.x /= scale_factor
    node.y /= scale_factor
end
meshes = triangular_partition(mesh)
regions = gmsh2matplotlibtri.(meshes)
# region = regions[1]

# Compute z-values on mesh points.
func(x, y) = [
    x^2 + y^2,
    x^2 - y^2
]
number_of_regions = length(regions)
z = Vector{Vector}(undef, number_of_regions)
t = Vector{Vector}(undef, number_of_regions)
for k = 1 : number_of_regions
    z[k] = getindex.(func.(regions[k].x, regions[k].y), 1)
    t[k] = getindex.(func.(regions[k].x, regions[k].y), 2)
end

# Refine region to evaluate interpolated function
subdiv = 3
x, y = refine(gmsh2matplotlibtri(mesh), subdiv, true)

# Compute initial and real values.
func0(xi, yi) = [0., 0.]

# Perform interpolation.
α7, α8, α11, α12 = fill(0.001, 4)
interpolant, coeffs = hiddenfis(regions, z, t, α7=α7, α8=α8, α11=α11, α12=α12, func0=func0, num_iter=10, get_coeffs=true)

# Interpolation and error
real_values = func.(x, y)
real_values_z = getindex.(real_values, 1)
real_values_t = getindex.(real_values, 2)
interpolated_values = interpolant.(x, y)
interpolated_values_z = getindex.(interpolated_values, 1)
interpolated_values_t = getindex.(interpolated_values, 2)
absolute_error_z = abs.(real_values_z - interpolated_values_z)
absolute_error_t = abs.(real_values_t - interpolated_values_t)
relative_error_z = absolute_error_z ./ (abs.(real_values_z .+ eps())) * 100
relative_error_t = absolute_error_t ./ (abs.(real_values_t .+ eps())) * 100
frerr_z = mean(absolute_error_z.^2)
frerr_t = mean(absolute_error_t.^2)

# Plot the results.
xd = xcoords(mesh)
yd = ycoords(mesh)
ztd = func.(xd, yd)
zd = getindex.(ztd, 1)
td = getindex.(ztd, 2)
minx, meanx, maxx = minimum(x), mean(x), maximum(x)
miny, meany, maxy = minimum(y), mean(y), maximum(y)

# Plot the results
figure("Real and Interpolated Values z")
plot_trisurf(x, y, real_values_z)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(real_values_z), maximum(real_values_z), length=5))
tight_layout()

plot_trisurf(x, y, interpolated_values_z)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(interpolated_values_z), maximum(interpolated_values_z), length=5))
tight_layout()
display(gcf())

figure("Absolute Values z")
plot_trisurf(x, y, absolute_error_z)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(absolute_error_z), maximum(absolute_error_z), length=5))
tight_layout()
display(gcf())

figure("Relative Error z")
plot_trisurf(x, y, relative_error_z)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(relative_error_z), maximum(relative_error_z), length=5))
tight_layout()
display(PyPlot.gcf())

figure("Real and Interpolated Values t")
plot_trisurf(x, y, real_values_t)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(real_values_t), maximum(real_values_t), length=5))
tight_layout()

plot_trisurf(x, y, interpolated_values_t)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(interpolated_values_t), maximum(interpolated_values_t), length=5))
tight_layout()
display(gcf())

figure("Absolute Values t")
plot_trisurf(x, y, absolute_error_t)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(absolute_error_t), maximum(absolute_error_t), length=5))
tight_layout()
display(gcf())

figure("Relative Error t")
plot_trisurf(x, y, relative_error_t)
colors = ["black", "red", "blue", "green", "purple", "magenta"]
for (color, region) in zip(colors, regions)
    triplot(region, color=color)
    plot(region.x, region.y, ".", color="orange")
end
xlabel("x", labelpad=10, fontsize=12)
ylabel("y", labelpad=10, fontsize=12)
zlabel("z", labelpad=10, fontsize=12)
xticks(range(minx, maxx, length=5))
yticks(range(miny, maxy, length=5))
# zticks(range(minimum(relative_error_t), maximum(relative_error_t), length=5))
tight_layout()
display(PyPlot.gcf())


text:
# ---
# title: 58. Length of Last Word
# id: problem58
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/length-of-last-word/description/>
# hidden: true
# ---
#
# Given a string `s` consists of some words separated by spaces, return _the
# length of the last word  in the string. If the last word does not exist,
# return _`0`.
#
# A **word** is a maximal substring consisting of non-space characters only.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "Hello World"
#     Output: 5
#
#
# **Example 2:**
#
#
#
#     Input: s = " "
#     Output: 0
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 104`
#   * `s` consists of only English letters and spaces `' '`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module Hangul

# Hangul Compatibility Jamo
export Jamo
export Letter, Consonant, Vowel
export Syllable
export ComposeError, DecomposeError
export to_char, to_syllable
include("jamo.jl")

# Hangul Jamo
export YetJamo
export Initial, Medial, Final
export YetSyllable
export to_string
include("yetjamo.jl")

end # module


text:
#=
# Plotting Graphs with `GraphMakie.jl`
## The `graphplot` Command
Plotting your first `AbstractGraph` from [`Graphs.jl`](https://github.com/JuliaGraphs/Graphs.jl)
is as simple as
=#
using CairoMakie
CairoMakie.activate!(type="png") # hide
set_theme!(resolution=(800, 400)) #hide
CairoMakie.inline!(true) # hide
using GraphMakie
using Graphs
import Random; Random.seed!(2) # hide

using GraphMakie.NetworkLayout

g = SimpleGraph(5)
add_edge!(g, 1, 2); add_edge!(g, 2, 4);
add_edge!(g, 4, 3); add_edge!(g, 3, 2);
add_edge!(g, 2, 5); add_edge!(g, 5, 4);
add_edge!(g, 4, 1); add_edge!(g, 1, 5);

offsets = 0.15 * (p[:node_pos][] .- p[:node_pos][][1])
offsets[1] = Point2f(0, 0.3)
p.nlabels_offset[] = offsets
autolimits!(ax)
f # hide

# ## Adding Edge Labels
Random.seed!(42)
g = barabasi_albert(6, 2)

labels =  repr.(1:ne(g))

f, ax, p = graphplot(g, elabels=labels,
                     elabels_color=[:black for i in 1:ne(g)],
                     edge_color=[:black for i in 1:ne(g)])
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
## Indicate Edge Direction

It is possible to put arrows on the edges using the `arrow_show` parameter. This parameter
is `true` for `SimpleDiGraph` by default. The position and size of each arrowhead can be
change using the `arrow_shift` and `arrow_size` parameters.
=#
g = wheel_digraph(10)
arrow_size = [10+i for i in 1:ne(g)]
arrow_shift = range(0.1, 0.8, length=ne(g))
f, ax, p = graphplot(g; arrow_size, arrow_shift)
hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

set_theme!(resolution=(800, 800)) #hide
g = SimpleGraph(8); add_edge!(g, 1, 2); add_edge!(g, 3, 4); add_edge!(g, 5, 6); add_edge!(g, 7, 8)

waypoints = Dict(1 => [(.25,  0.25), (.75, -0.25)],
                 2 => [(.25, -0.25), (.75, -0.75)],
                 3 => [(.25, -0.75), (.75, -1.25)],
                 4 => [(.25, -1.25), (.75, -1.75)])
waypoint_radius = Dict(1 => nothing,
                       2 => 0,
                       3 => 0.05,
                       4 => 0.15)

f = Figure(); f[1,1] = ax = Axis(f)
using Makie.Colors # hide
for i in 3:4 #hide
    poly!(ax, Circle(Point2f(waypoints[i][1]), waypoint_radius[i]), color=RGBA(0.0,0.44705883,0.69803923,0.2)) #hide
    poly!(ax, Circle(Point2f(waypoints[i][2]), waypoint_radius[i]), color=RGBA(0.0,0.44705883,0.69803923,0.2)) #hide
end #hide

p = graphplot!(ax, g; layout=SquareGrid(cols=2, dy=-0.5),
               waypoints, waypoint_radius,
               nlabels=["","r = nothing (equals :spline)",
                        "","r = 0 (straight lines)",
                        "","r = 0.05 (in data space)",
                        "","r = 0.1"],
               nlabels_distance=30, nlabels_align=(:left,:center))

for i in 1:4 #hide
    scatter!(ax, waypoints[i], color=RGBA(0.0,0.44705883,0.69803923,1.0)) #hide
end #hide
xlims!(ax, (-0.1, 2.25)), hidedecorations!(ax); hidespines!(ax); ax.aspect = DataAspect()
f # hide

#=
## Plot Graphs in 3D
If the layout returns points in 3 dimensions, the plot will be in 3D. However this is a bit
experimental. Feel free to file an issue if there are any problems.
=#
set_theme!(resolution=(800, 800)) #hide
g = smallgraph(:cubical)
elabels_shift = [0.5 for i in 1:ne(g)]
elabels_shift[[2,7,8,9]] .= 0.3
elabels_shift[10] = 0.25
graphplot(g; layout=Spring(dim=3, seed=5),
          elabels="Edge ".*repr.(1:ne(g)),
          elabels_textsize=12,
          elabels_opposite=[3,5,7,8,12],
          elabels_shift,
          elabels_distance=3,
          arrow_show=true,
          arrow_shift=0.9,
          arrow_size=15)

#=
Using [`JSServe.jl`](https://github.com/SimonDanisch/JSServe.jl) and [`WGLMakie.jl`](https://github.com/JuliaPlots/WGLMakie.jl)
we can also add some interactivity:
=#
using JSServe
Page(exportable=true, offline=true)
#
using WGLMakie
WGLMakie.activate!()
set_theme!(resolution=(800, 600))
g = smallgraph(:dodecahedral)
graphplot(g, layout=Spring(dim=3), node_size=100)


text:
struct UnitDict{R,T}
    keys::R
    values::Vector{T}

    UnitDict(unitrange, values::Vector{T}) where T = new{typeof(unitrange),T}(unitrange, values)
    UnitDict(values::Vector) = UnitDict(Base.OneTo(length(values)), values)
    function UnitDict(dict::Dict{T1,T2}) where {T1,T2}
        if isempty(dict)
            UnitDict(zero(T1):zero(T1), T2[])
        else
            a = minimum(keys(dict))
            b = maximum(keys(dict))

            a == 1 ? (range = Base.OneTo(b)) : (range = a:b)
            values = T2[]
            for k in range
                push!(values, dict[k])
            end

            UnitDict(range, values)
        end
    end
end

function Base.show(io::IO, mime::MIME{Symbol("text/plain")}, dict::UnitDict{R,T}) where {R,T}
    summary(io, dict)
    println(io,"")
    println(io, " keys:   "*string(dict.keys))
    println(io, " values: "*string(dict.values))
end

@inline Base.length(dict::UnitDict) = length(dict.values)
@inline Base.firstindex(dict::UnitDict) = first(dict.keys)
@inline Base.lastindex(dict::UnitDict) = last(dict.keys)

@inline Base.getindex(dict::UnitDict{Base.OneTo{K},T}, key::K) where {T,K} = dict.values[key]
@inline Base.getindex(dict::UnitDict{UnitRange{K},T}, key::K) where {T,K} = dict.values[key - first(dict.keys) + 1]
@inline Base.getindex(dict::UnitDict{Base.OneTo{K},T}, key::AbstractVector{K}) where {T,K} = dict.values[key]
@inline Base.getindex(dict::UnitDict{UnitRange{K},T}, key::AbstractVector{K}) where {T,K} = dict.values[key .- (first(dict.keys) + 1)]
@inline Base.getindex(dict::UnitDict{Base.OneTo{K},T}, key::UnitRange{K}) where {T,K} = dict.values[key]
@inline Base.getindex(dict::UnitDict{UnitRange{K},T}, key::UnitRange{K}) where {T,K} = dict.values[key .- (first(dict.keys) + 1)]
@inline Base.setindex!(dict::UnitDict{Base.OneTo{K},T}, value::T, key::K) where {T,K} = setindex!(dict.values, value, key)
@inline Base.setindex!(dict::UnitDict{UnitRange{K},T}, value::T, key::K) where {T,K} = setindex!(dict.values, value, key - first(dict.keys) + 1)

@inline Base.iterate(d::UnitDict, i = 1) = i > length(d) ? nothing : (d.values[i], i + 1)
@inline Base.pairs(d::UnitDict, i = 1) = Base.Generator(=>, d.keys, d.values)

@inline Base.foreach(f, itr::UnitDict, arg...) = (for x in itr; f(x, arg...); end; return)
@inline function foreachactive(f, itr::UnitDict, arg...)
    for x in itr
        isinactive(x) && continue
        f(x, arg...)
    end
    return
end

@inline Base.haskey(d::UnitDict, key) = Base.ht_keyindex(d, key)

function Base.ht_keyindex(h::UnitDict, key)
    for el in h.keys
        if el == key
            return true
        end
    end
    return false
end


text:
include("Graph.jl")

function LRtest(G; needEmbedding = false)
    if (G.n > 2) && (G.m > 3*G.n-6)
        return nothing
    end
    INF = 131072

    g = Array{Array{Int, 1}, 1}(undef, G.n)
    gID = Array{Array{Int, 1}, 1}(undef, G.n)
    for i = 1 : G.n
        g[i] = []
        gID[i] = []
    end
    for (ID, u, v, w) in G.E
        push!(g[u], v)
        push!(gID[u], ID)
        push!(g[v], u)
        push!(gID[v], ID)
    end

    Roots = Array{Int, 1}()
    height = zeros(Int, G.n)
    lowpt = zeros(Int, G.m)
    lowpt2 = zeros(Int, G.m)
    nestingDepth = zeros(Int, G.m)
    fill!(height, INF)
    parentEdge = zeros(Int, G.n)
    ori = zeros(Int, G.m)

    ref = zeros(Int, G.m)
    side = ones(Int, G.m)
    S = zeros(Int, 2*G.m, 4)
    top = 0
    stackBottom = zeros(Int, G.m)
    lowptEdge = zeros(Int, G.m)

    # orientation
    DFS1(v) = begin
        e = parentEdge[v]
        for i = 1 : size(g[v], 1)
            w = g[v][i]
            wid = gID[v][i]
            if ori[wid] == 0
                ori[wid] = (v == G.E[wid][2]) ? 1 : -1
                lowpt[wid] = height[v]
                lowpt2[wid] = height[v]
                if height[w] == INF
                    parentEdge[w] = wid
                    height[w] = height[v] + 1
                    DFS1(w)
                else
                    lowpt[wid] = height[w]
                end
                # determine nesting depth
                nestingDepth[wid] = 2*lowpt[wid]
                if lowpt2[wid] < height[v]
                    nestingDepth[wid] += 1
                end
                # update lowpoints of parent edge e
                if e != 0
                    if lowpt[wid] < lowpt[e]
                        lowpt2[e] = min(lowpt[e], lowpt2[wid])
                        lowpt[e] = lowpt[wid]
                    elseif lowpt[wid] > lowpt[e]
                        lowpt2[e] = min(lowpt2[e], lowpt[wid])
                    else
                        lowpt2[e] = min(lowpt2[e], lowpt2[wid])
                    end
                end
            end
        end
    end

    for i = 1 : G.n
        if height[i] == INF
            push!(Roots, i)
            height[i] = 0
            DFS1(i)
        end
    end

    #println("Lowpt :")
    #foreach(i -> print(lowpt[i], " "), 1 : G.m)
    #println()
    #println("Lowpt2 :")
    #foreach(i -> print(lowpt2[i], " "), 1 : G.m)
    #println()
    #println("Height :")
    #foreach(i -> print(height[i], " "), 1 : G.n)
    #println()
    #for i = 1 : G.m
    #    if ori[i] == 1
    #        println(G.E[i][2], " -> ", G.E[i][3])
    #    elseif ori[i] == -1
    #        println(G.E[i][3], " -> ", G.E[i][2])
    #    end
    #end

    # testing

    ep = Array{Array{Int, 1}, 1}(undef, G.n)
    eptar = Array{Array{Int, 1}, 1}(undef, G.n)
    foreach(i -> ep[i] = [], 1 : G.n)
    foreach(i -> eptar[i] = [], 1 : G.n)
    for (ID, u, v, w) in G.E
        if ori[ID] == 1
            push!(ep[u], ID)
        else
            push!(ep[v], ID)
        end
    end
    for i = 1 : G.n
        sort!(ep[i], by = x->nestingDepth[x])
        for j = 1 : size(ep[i], 1)
            if ori[ep[i][j]] == 1
                push!(eptar[i], G.E[ep[i][j]][3])
            else
                push!(eptar[i], G.E[ep[i][j]][2])
            end
        end
    end

    #for i = 1 : G.n
    #    for j = 1 : size(ep[i], 1)
    #        print(eptar[i][j], "[", ep[i][j], "](", nestingDepth[ep[i][j]], ") ")
    #    end
    #    println()
    #end

    conflicting(I, b) = begin
        return ((I != (0, 0)) && (lowpt[I[2]] > lowpt[b]))
    end
    lowest(P) = begin
        if (P[1] == 0) && (P[2] == 0)
            return lowpt[P[3]]
        elseif (P[3] == 0) && (P[4] == 0)
            return lowpt[P[1]]
        end
        return min(lowpt[P[1]], lowpt[P[3]])
    end
    target(x) = begin
        if ori[x] == 1
            return G.E[x][3]
        else
            return G.E[x][2]
        end
    end
    swapLR(Q) = begin
        return [Q[3], Q[4], Q[1], Q[2]]
    end

    DFS2(v) = begin
        e = parentEdge[v]
        for i = 1 : size(ep[v], 1)
            eiid = ep[v][i]
            eitar = eptar[v][i]
            stackBottom[eiid] = top
            if eiid == parentEdge[eitar]
                if !DFS2(eitar)
                    return false
                end
            else
                lowptEdge[eiid] = eiid
                top += 1
                S[top, :] = [0, 0, eiid, eiid]
            end
            # integrate new return edges
            if lowpt[eiid] < height[v]
                if i == 1
                    lowptEdge[e] = lowptEdge[eiid]
                else
                    # add constraints of ei
                    P = zeros(Int, 4)
                    # merge return edges of ei into P.R
                    while true
                        Q = copy(S[top, :])
                        top -= 1
                        ((Q[1] != 0) || (Q[2] != 0)) ? (Q = copy(swapLR(Q))) : nothing
                        if (Q[1] != 0) || (Q[2] != 0)
                            return false
                        else
                            if lowpt[Q[3]] > lowpt[e]
                                if (P[3] == 0) && (P[4] == 0)
                                    P[4] = Q[4]
                                else
                                    ref[P[3]] = Q[4]
                                end
                                P[3] = Q[3]
                            else
                                ref[Q[3]] = lowptEdge[e]
                            end
                        end
                        if (top == stackBottom[eiid])
                            break
                        end
                    end
                    # merge confliting return edges of e1, ..., ei-1 into P.L
                    #println(eiid, " ", S[top, 1], " ", S[top, 2], " ", S[top, 3], " ", S[top, 4])
                    while ((conflicting((S[top, 1], S[top, 2]), eiid)) || (conflicting((S[top, 3], S[top, 4]), eiid)))
                        Q = copy(S[top, :])
                        top -= 1
                        conflicting((Q[3], Q[4]), eiid) ? (Q = copy(swapLR(Q))) : nothing
                        if conflicting((Q[3], Q[4]), eiid)
                            return false
                        else
                            ref[P[3]] = Q[4]
                            (Q[3] != 0) ? (P[3] = Q[3]) : nothing
                        end
                        if (P[1] == 0) && (P[2] == 0)
                            P[2] = Q[2]
                        else
                            ref[P[1]] = Q[2]
                        end
                        P[1] = Q[1]
                    end
                    if sum(P) != 0
                        top += 1
                        S[top, :] = copy(P)
                    end
                end
            end
        end
        # remove back edges returning to parent
        if e != 0
            u = (G.E[e][2] != v) ? (G.E[e][2]) : (G.E[e][3])
            # trim back edges ending at parent u
            ## drop entire conflict pairs
            while (top > 0) && (lowest(S[top, :]) == height[u])
                P = copy(S[top, :])
                top -= 1
                (P[1] != 0) ? (side[P[1]] = -1) : nothing
            end
            if top > 0
                P = copy(S[top, :])
                top -= 1
                # trim left interval
                while (P[2] != 0) && (target(P[2]) == u)
                    P[2] = ref[P[2]]
                end
                if (P[2] == 0) && (P[1] != 0)
                    ref[P[1]] = P[3]
                    side[P[1]] = -1
                    P[1] = 0
                end
                # trim right interval
                while (P[4] != 0) && (target(P[4]) == u)
                    P[4] = ref[P[4]]
                end
                if (P[4] == 0) && (P[3] != 0)
                    ref[P[3]] = P[1]
                    side[P[3]] = -1
                    P[3] = 0
                end
                top += 1
                S[top, :] = copy(P)
            end
            # side of e is side of a highest return edge
            if lowpt[e] < height[u]
                hl = S[top, 2]
                hr = S[top, 4]
                if (hl != 0) && ((hr == 0) || (lowpt[hl] > lowpt[hr]))
                    ref[e] = hl
                else
                    ref[e] = hr
                end
            end
        end
        return true
    end

    for s in Roots
        if !DFS2(s)
            return nothing
        end
    end

    #foreach(i -> print(side[i], " "), 1 : G.m)
    #println()

    Emb = Array{Array{Int, 1}, 1}(undef, G.n)
    foreach(i -> Emb[i] = [], 1 : G.n)

    if (needEmbedding == false)
        return Emb
    end

    # get embedding

    sign(eid) = begin
        if ref[eid] != 0
            side[eid] = side[eid] * sign(ref[eid])
            ref[eid] = 0
        end
        return side[eid]
    end

    foreach(i -> nestingDepth[i] = sign(i) * nestingDepth[i], 1 : G.m)

    foreach(i -> ep[i] = [], 1 : G.n)
    foreach(i -> eptar[i] = [], 1 : G.n)
    for (ID, u, v, w) in G.E
        if ori[ID] == 1
            push!(ep[u], ID)
        else
            push!(ep[v], ID)
        end
    end
    for i = 1 : G.n
        sort!(ep[i], by = x->nestingDepth[x])
        for j = 1 : size(ep[i], 1)
            if ori[ep[i][j]] == 1
                push!(eptar[i], G.E[ep[i][j]][3])
            else
                push!(eptar[i], G.E[ep[i][j]][2])
            end
        end
    end

    sp = 0
    LL = zeros(Int, 3*G.m)
    RR = zeros(Int, 3*G.m)
    Val = zeros(Int, 3*G.m)
    Fst = zeros(Int, G.n)
    rightRef = zeros(Int, G.n)
    leftRef = zeros(Int, G.n)
    lkid = Array{Array{Int, 1}, 1}(undef, G.n)
    foreach(i -> lkid[i] = [], 1 : G.n)

    AddToFirst(v, x) = begin
        sp += 1
        RR[sp] = Fst[v]
        (Fst[v] != 0) ? (LL[Fst[v]] = sp) : nothing
        Val[sp] = x
        Fst[v] = sp
    end

    AddBefore(v, id, x) = begin
        sp += 1
        RR[sp] = id
        LL[sp] = LL[id]
        (LL[id] != 0) ? (RR[LL[id]] = sp) : nothing
        (Fst[v] == id) ? (Fst[v] = sp) : nothing
        LL[id] = sp
        Val[sp] = x
    end

    AddAfter(id, x) = begin
        sp += 1
        RR[sp] = RR[id]
        LL[sp] = id
        (RR[id] != 0) ? (LL[RR[id]] = sp) : nothing
        RR[id] = sp
        Val[sp] = x
    end

    for i = 1 : G.n
        if size(ep[i], 1) > 0
            AddToFirst(i, ep[i][1])
            push!(lkid[i], sp)
            for j = 2 : size(ep[i], 1)
                AddAfter(sp, ep[i][j])
                push!(lkid[i], sp)
            end
        end
    end

    DFS3(v) = begin
        for i = 1 : size(ep[v], 1)
            eiid = ep[v][i]
            w = eptar[v][i]
            eilkid = lkid[v][i]
            if eiid == parentEdge[w]
                AddToFirst(w, eiid)
                leftRef[v] = eilkid
                rightRef[v] = eilkid
                DFS3(w)
            else
                if side[eiid] == 1
                    AddAfter(rightRef[w], eiid)
                else
                    AddBefore(w, leftRef[w], eiid)
                    leftRef[w] = sp
                end
            end
        end
    end

    for s in Roots
        DFS3(s)
    end

    for i = 1 : G.n
        x = Fst[i]
        while x != 0
            push!(Emb[i], Val[x])
            x = RR[x]
        end
    end

    return Emb
end


text:
function calc_prewhitening_filter(Rxx)
    F = eigen(Rxx)
    F.vectors * Diagonal(1 ./ sqrt.(F.values)) * F.vectors' .* sqrt(mean_power(F.values))
end

function calc_amplitude_filter(Rxx)
    F = eigen(Rxx)
    F.vectors * Diagonal(1 ./ F.values) * F.vectors' .* mean_power(F.values)
end

function calc_eigen_beamformer(Rxx)
    eigvecs(Rxx)[:,end]
end

function calc_variance_covariance(signal)
    Hermitian(signal' * signal) / size(signal, 1)
end

function mean_power(eigenvalues)
    1 / (sum(1 ./ eigenvalues) / length(eigenvalues))
end