text:
# A SortedMultiDict is a wrapper around balancedTree.
## Unlike SortedDict, a key in SortedMultiDict can
## refer to multiple data entries.

@compat type SortedMultiDict{K, D, Ord <: Ordering}
    bt::BalancedTree23{K,D,Ord}

## Zero-argument constructor, or possibly one argument to specify order.

    function (::Type{SortedMultiDict{K,D,Ord}}){K, D, Ord <: Ordering}(o::Ord=Forward)
        bt1 = BalancedTree23{K,D,Ord}(o)
        new{K,D,Ord}(bt1)
    end
end

const SMDSemiToken = IntSemiToken

const SMDToken = Tuple{SortedMultiDict, IntSemiToken}

## This constructor takes two arrays an ordering object which defaults
## to Forward

function SortedMultiDict{K,D, Ord <: Ordering}(kk::AbstractArray{K,1},
                                               dd::AbstractArray{D,1},
                                               o::Ord=Forward)
    if length(kk) != length(dd)
        throw(ArgumentError("SortedMultiDict K and D constructor array arguments must be the same length"))
    end
    h = SortedMultiDict{K,D,Ord}(o)
    for i = 1 : length(kk)
        insert!(h, kk[i], dd[i])
    end
    h
end

## Take pairs and infer argument
## types.  Note:  this works only for the Forward ordering.

function SortedMultiDict{K,D}(p1::Pair{K,D}, ps::Pair{K,D}...)
    h = SortedMultiDict{K,D,ForwardOrdering}()
    insert!(h, p1.first, p1.second)
    for p in ps
        insert!(h, p.first, p.second)
    end
    h
end

## Take pairs and infer argument
## types.  Ordering parameter must be explicit first argument.

function SortedMultiDict{K,D, Ord <: Ordering}(o::Ord, ps::Pair{K,D}...)
    h = SortedMultiDict{K,D,Ord}(o)
    for p in ps
        insert!(h, p.first, p.second)
    end
    h
end

## This one takes an iterable; ordering type is optional.

SortedMultiDict{Ord <: Ordering}(kv, o::Ord=Forward) =
sortedmultidict_with_eltype(kv, eltype(kv), o)

function sortedmultidict_with_eltype{K,D,Ord}(kv, ::Type{Pair{K,D}}, o::Ord)
    h = SortedMultiDict{K,D,Ord}(o)
    for (k,v) in kv
        insert!(h, k, v)
    end
    h
end

## This function inserts an item into the tree.
## It returns a token that
## points to the newly inserted item.

@inline function insert!{K, D, Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}, k_, d_)
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), true)
    IntSemiToken(i)
end

## push! is an alternative to insert!; it returns the container.

@inline function push!{K, D, Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}, pr::Pair{K,D})
    insert!(m.bt, convert(K,pr[1]), convert(D,pr[2]), true)
    m
end

## First and last return the first and last (key,data) pairs
## in the SortedMultiDict.  It is an error to invoke them on an
## empty SortedMultiDict.

@inline function first(m::SortedMultiDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

@inline function last(m::SortedMultiDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return Pair(m.bt.data[i].k, m.bt.data[i].d)
end

function searchequalrange(m::SortedMultiDict, k_)
    k = convert(keytype(m),k_)
    i1 = findkeyless(m.bt, k)
    i2, exactfound = findkey(m.bt, k)
    if exactfound
        i1a = nextloc0(m.bt, i1)
        return IntSemiToken(i1a), IntSemiToken(i2)
    else
        return IntSemiToken(2), IntSemiToken(1)
    end
end

## '(k,d) in m' checks whether a key-data pair is in
## a sorted multidict.  This requires a loop over
## all data items whose key is equal to k.

function in_(k_, d_, m::SortedMultiDict)
    k = convert(keytype(m), k_)
    d = convert(valtype(m), d_)
    i1 = findkeyless(m.bt, k)
    i2,exactfound = findkey(m.bt,k)
    !exactfound && return false
    ord = m.bt.ord
    while true
        i1 = nextloc0(m.bt, i1)
        @assert(eq(ord, m.bt.data[i1].k, k))
        m.bt.data[i1].d == d && return true
        i1 == i2 && return false
    end
end

@inline eltype{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}) =  Pair{K,D}
@inline eltype{K,D,Ord <: Ordering}(::Type{SortedMultiDict{K,D,Ord}}) =  Pair{K,D}
@inline in(pr::Pair, m::SortedMultiDict) =
    in_(pr[1], pr[2], m)
@inline in(::Tuple{Any,Any}, ::SortedMultiDict) =
    throw(ArgumentError("'(k,v) in sortedmultidict' not supported in Julia 0.4 or 0.5.  See documentation"))

@inline keytype{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}) = K
@inline keytype{K,D,Ord <: Ordering}(::Type{SortedMultiDict{K,D,Ord}}) = K
@inline valtype{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}) = D
@inline valtype{K,D,Ord <: Ordering}(::Type{SortedMultiDict{K,D,Ord}}) = D
@inline ordtype{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord}) = Ord
@inline ordtype{K,D,Ord <: Ordering}(::Type{SortedMultiDict{K,D,Ord}}) = Ord
@inline orderobject(m::SortedMultiDict) = m.bt.ord

@inline function haskey(m::SortedMultiDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    exactfound
end

## Check if two SortedMultiDicts are equal in the sense of containing
## the same (K,D) pairs in the same order.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedMultiDict, m2::SortedMultiDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        throw(ArgumentError("Cannot use isequal for two SortedMultiDicts unless their element types and ordering objects are equal"))
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1,d1 = deref((m1,p1))
        k2,d2 = deref((m2,p2))
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

const SDorAssociative = Union{Associative,SortedMultiDict}

function mergetwo!{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord},
                                        m2::SDorAssociative)
    for (k,v) in m2
        insert!(m.bt, convert(K,k), convert(D,v), true)
    end
end

function packcopy{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord})
    w = SortedMultiDict((K)[], (D)[], orderobject(m))
    mergetwo!(w,m)
    w
end

function packdeepcopy{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord})
    w = SortedMultiDict((K)[], (D)[], orderobject(m))
    for (k,v) in m
        insert!(w.bt, deepcopy(k), deepcopy(v), true)
    end
    w
end

function merge!{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord},
                                     others::SDorAssociative...)
    for o in others
        mergetwo!(m,o)
    end
end

function merge{K,D,Ord <: Ordering}(m::SortedMultiDict{K,D,Ord},
                                    others::SDorAssociative...)
    result = packcopy(m)
    merge!(result, others...)
    result
end

function Base.show{K,D,Ord <: Ordering}(io::IO, m::SortedMultiDict{K,D,Ord})
    print(io, "SortedMultiDict(")
    print(io, orderobject(m), ",")
    l = length(m)
    for (count,(k,v)) in enumerate(m)
        print(io, k, " => ", v)
        if count < l
            print(io, ", ")
        end
    end
    print(io, ")")
end

similar{K,D,Ord<:Ordering}(m::SortedMultiDict{K,D,Ord}) =
   SortedMultiDict(K[], D[], orderobject(m))


text:
JULIA_PREFIX = abspath(joinpath(Base.source_path(), "..", "..", ".."))

if !("JULIA_PKGDIR" in keys(ENV))
    ENV["JULIA_PKGDIR"] = joinpath(JULIA_PREFIX, "share", "julia", "site")
    Base.LOAD_CACHE_PATH[1] = joinpath(ENV["JULIA_PKGDIR"], "lib", string("v", join(split(string(VERSION), ".")[1:2], ".")))
end

ext: jl

max_stars_repo_name: mariusvniekerk/julia-feedstock

max_issues_repo_name: mariusvniekerk/julia-feedstock

max_forks_repo_name: mariusvniekerk/julia-feedstock


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

#### Pretty printing for MVTSeries

function Base.summary(io::IO, x::MVTSeries)
    et = eltype(x) === Float64 ? "" : ",$(eltype(x))"
    typestr = "MVTSeries{$(frequencyof(x))$(et)}"

    cols = axes(x,2)
    ncols = length(cols)
    if ncols == 0
        vstr = "no variables"
    else
        vstr="variables ($(cols[1])"
        for i = 2:ncols
            cstr = ",$(cols[i])"
            if length(vstr) + length(cstr) > 20
                vstr *= ",…"
                break
            else
                vstr *= cstr
            end
        end
        vstr *= ")"
    end

    print(io, size(x,1), "×", size(x, 2), " ", typestr, " with range ", rangeof(x), " and ", vstr)
    nothing
end

Base.show(io::IO, ::MIME"text/plain", x::MVTSeries) = show(io, x)
function Base.show(io::IO, x::MVTSeries)
    summary(io, x)
    isempty(x) && return
    print(io, ":")
    limit = get(io, :limit, true)
    nval, nsym = size(x)

    from = firstdate(x)
    dheight, dwidth = displaysize(io)
    if get(io, :compact, nothing) === nothing
        io = IOContext(io, :compact => true)
    end
    dwidth -= 11

    names = collect(colnames(x))
    vals = _vals(x)

    names_str = String[]
    for n in names
        sn = string(n)
        push!(names_str, length(sn) < 10 ? sn : sn[1:10])
    end
    sd_with_names = [reshape(names_str, 1, :); x.values]
    A = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)

    all_cols = true
    if length(A) ≠ nsym
        dwidth = div(dwidth - 1, 2)
        AL = Base.alignment(io, sd_with_names, axes(sd_with_names, 1), 1:nsym, dwidth, dwidth, 2)
        AR = reverse(Base.alignment(io, sd_with_names, axes(sd_with_names, 1), reverse(1:nsym), dwidth, dwidth, 2))
        Linds = [1:length(AL)...]
        Rinds = [nsym - length(AR) + 1:nsym...]
        all_cols = false
    end

    local vdots = "\u22ee"
    local hdots = " \u2026 "
    local ddots = " \u22f1 "

    print_aligned_val(io, v, (al, ar), showsep=true; sep=showsep ? "  " : "") = begin
        sv = sprint(print, v, context=io, sizehint=0)
        if v isa Number
            vl, vr = Base.alignment(io, v)
        else
            if length(sv) > al + ar
                sv = sv[1:al + ar - 1] * '…'
            end
            vl, vr = al, length(sv) - al
        end
        print(io, repeat(" ", al - vl), sv, repeat(" ", ar - vr), sep)
    end

    print_colnames(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
            print_aligned_val(io, names[col], align, i < nLcols)
        end
        nRcols == 0 && return
        print(io, hdots)
        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
            print_aligned_val(io, names[col], align, i < nRcols)
        end
    end

    print_rows(io, rows, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for row in rows
            mit = from + (row - 1)
            print(io, '\n', lpad(mit, 8), " : ")
            for (i, (val, align)) in enumerate(zip(vals[row, Lcols], LAligns))
                print_aligned_val(io, val, align, i < nLcols)
            end
            nRcols == 0 && continue
            print(io, hdots)
            for (i, (val, align)) in enumerate(zip(vals[row, Rcols], RAligns))
                print_aligned_val(io, val, align, i < nRcols)
            end
        end
    end

    print_vdots(io, Lcols, LAligns, Rcols=[], RAligns=[]) = begin
        print(io, '\n', repeat(" ", 11))
        local nLcols = length(Lcols)
        local nRcols = length(Rcols)
        for (i, (col, align)) in enumerate(zip(Lcols, LAligns))
            print_aligned_val(io, vdots, align, i < nLcols)
        end
        nRcols == 0 && return
        print(io, ddots)
        for (i, (col, align)) in enumerate(zip(Rcols, RAligns))
            print_aligned_val(io, vdots, align, i < nRcols)
        end
    end

    if !limit
        print(io, "\n", repeat(" ", 11))
        print_colnames(io, 1:nsym, A)
        print_rows(io, 1:nval, 1:nsym, A)
    elseif nval > dheight - 6 # all rows don't fit
                # unable to show all rows
        if all_cols
            print(io, "\n", repeat(" ", 11))
            print_colnames(io, 1:nsym, A)
            top = div(dheight - 6, 2)
            print_rows(io, 1:top, 1:nsym, A)
            print_vdots(io, 1:nsym, A)
            bot = nval - dheight + 7 + top
            print_rows(io, bot:nval, 1:nsym, A)
        else # not all_cols
            print(io, "\n", repeat(" ", 11))
            print_colnames(io, Linds, AL, Rinds, AR)
            top = div(dheight - 6, 2)
            print_rows(io, 1:top, Linds, AL, Rinds, AR)
            print_vdots(io, Linds, AL, Rinds, AR)
            bot = nval - dheight + 7 + top
            print_rows(io, bot:nval, Linds, AL, Rinds, AR)
        end # all_cols
    else # all rows fit
        if all_cols
            print(io, '\n', repeat(" ", 11))
            print_colnames(io, 1:nsym, A)
            print_rows(io, 1:nval, 1:nsym, A)
        else
            print(io, '\n', repeat(" ", 11))
            print_colnames(io, Linds, AL, Rinds, AR)
            print_rows(io, 1:nval, Linds, AL, Rinds, AR)
        end
    end
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

version = v"3.3"

# Collection of sources required to build glfw
sources = [
    "https://github.com/glfw/glfw.git" =>
    "b0796109629931b6fa6e449c15a177845256a407",

]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd glfw/
cmake -DCMAKE_INSTALL_PREFIX=$prefix -DCMAKE_TOOLCHAIN_FILE=/opt/$target/$target.toolchain -DBUILD_SHARED_LIBS=ON -DGLFW_BUILD_EXAMPLES=false -DGLFW_BUILD_TESTS=false
make -j${ncore}
make install
exit

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    MacOS(:x86_64),
    Windows(:x86_64),
    Windows(:i686),
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libglfw", :libglfw)
]

# Dependencies that must be installed before this package can be built
dependencies = []

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, "glfw", version, sources, script, platforms, products, dependencies)


text:
import Base.isless

"""
Species has fields id::String and na::Int64
and represents a metabolite/species
"""
struct Species
    id::String
    na::Int64 # number of atoms
end

function Species(id::String)
    return Species(id,0)
end
function isless(x::Species,y::Species)
    return x.id < y.id
end

"""
Flux has fields:
  id::String , E::LC{String} , P::LC{String} , transition::Dict{RAS,PAS}
and represents a flux/reaction
"""
struct Flux
    id::String
    E::LC{String}
    P::LC{String}
    transition::Dict{RAS,PAS}

    function Flux(id::String,E::LC{String},P::LC{String},transition::Dict{RAS,PAS})
        E2 = E
        P2 = P
        if( any( collect(values(E)).<0) || any( collect(values(P)).<0) )
            error("Flux Constructor : Stoich. with neg. entries..")
        end
        if( isempty(E) && isempty(P) )
            error("Flux Constructor : Stoich. without any entries..")
        end
        if( isempty(E) )
            E2 = LC{String}( "*" => 1.0)
        end
        if( isempty(P) )
            P2 = LC{String}( "*" => 1.0)
        end
        new(id,E2,P2,transition)
    end
end

function Flux(id::String,stoichiometry::LC{String})
    # auto split stoich. into positive and negtive parts:
    lce = LC{String}()
    lcp = LC{String}()
    for ki in keys(stoichiometry)
        if(stoichiometry[ki]>0)
            lcp[ki] = stoichiometry[ki]
        else
            lce[ki] = -stoichiometry[ki]
        end
    end
    return Flux(id,lce,lcp,Dict{RAS,PAS}())
end

function isless(x::Flux,y::Flux)
    return x.id < y.id
end
function isioflux(f::Flux)
    return ( isinflux(f) || isoutflux(f) )
end
function isinflux(f::Flux)
    return haskey(f.E, "*")
end
function isoutflux(f::Flux)
    return haskey(f.P, "*")
end

function create_reverse_flux( f::Flux )
    name_bw_flux = string(f.id,"_bw")
    return create_reverse_flux( f , name_bw_flux )
end
function create_reverse_flux( f::Flux , name_bw_flux::String)
    tmap_reverse = Dict{RAS,PAS}()
    for ri in keys(f.transition)
        pi = f.transition[ri]
        tmap_reverse[ RAS(name_bw_flux , pi.sas , pi.ppos) ] = PAS(name_bw_flux,ri.sas,ri.rpos)
    end
    return Flux(name_bw_flux,f.P,f.E,tmap_reverse)
end

"""
Network has fields:
  fluxes::Vector{String} , species::Vector{String}
  fluxes_db::Dict{String,Flux} , species_db::Dict{String,Species}
and represents a flux/reaction
"""
struct Network
    fluxes::Vector{String}
    species::Vector{String}
    fluxes_db::Dict{String,Flux}
    species_db::Dict{String,Species}
end

function Network()
    Network( Vector{Flux}() )
end

"""
sas(id_metabolite::String,n::Network)
returns Vector{SAS}
computes all SAS for a given metabolite
"""
function sas(ss::String,n::Network)
    vsas::Vector{SAS} = Vector{SAS}()
    s::Species = n.species_db[ss]
    print(s)
    for i in 1:s.na
        push!( vsas , SAS(ss,i) )
    end
    return vsas
end

"""
ras(id_flux::String,n::Network)
returns Vector{RAS}
computes all RAS for a given flux
"""
function ras(fs::String,n::Network)
    f::Flux = n.fluxes_db[fs]

    if(isinflux(f)); return Vector{RAS}(); end

    e::Vector{String} = collect( keys( f.E ) )
    sort(e)
    vras::Vector{RAS} = Vector{RAS}()
    for ei::String in e
        if(n.species_db[ei].na>0)
            nrf  = f.E[ei]
            nr::Int64 = Int64(round(nrf))
            if( abs(nrf-nr) > (1e-8) )
                error("Isotopomer Species with non-integer number of reactands..")
            end
            for ri in 1:nr
                for rsas in sas(ei,n)
                    push!( vras , RAS( fs , rsas , ri ) )
                end
            end
        end
    end
    return vras
end

"""
pas(id_flux::String,n::Network)
returns Vector{PAS}
computes all PAS for a given flux
"""
function pas(fs::String,n::Network)
    f::Flux = n.fluxes_db[fs]

    if(isoutflux(f)); return Vector{PAS}(); end

    p::Vector{String} = collect( keys( f.P ) )
    sort(p)
    vpas::Vector{PAS} = Vector{PAS}()
    for pi::String in p
        if(n.species_db[pi].na>0)
            nrf  = f.P[pi]
            nr::Int64 = Int64(round(nrf))
            if( abs(nrf-nr) > (1e-8) )
                error("Isotopomer Species with non-integer number of reactands..")
            end
            for ri in 1:nr
                for rsas in sas(pi,n)
                    push!( vpas , PAS( fs , rsas , ri ) )
                end
            end
        end
    end
    return vpas
end

"""
ras(n::Network)
returns Vector{SAS}
computes all SAS of all species of the network
"""
function sas(n::Network)
    all_sas::Set{SAS} = Set{SAS}()
    for si in n.species
        all_sas = union( all_sas , Set(sas(si,n)) )
    end
    return all_sas
end

function sas(csemu::ConfiguredSEMU)
    return Set( csemu.semu )
end

"""
ras(n::Network)
returns Vector{RAS}
computes all RAS of all fluxes of the network
"""
function ras(n::Network)
    all_ras::Set{RAS} = Set{RAS}()
    for fi in n.fluxes
        all_ras = union( all_ras , Set(ras(fi,n)) )
    end
    return all_ras
end

"""
pas(n::Network)
returns Vector{PAS}
computes all PAS of all fluxes of the network
"""
function pas(n::Network)
    all_pas::Set{PAS} = Set{PAS}()
    for fi in n.fluxes
        all_pas = union( all_pas , Set(pas(fi,n)) )
    end
    return all_pas
end

"""
  get_N( n::Network ; sparse=false)
returns Array{Float64,2}
computes the stoichiometric matrix of the network
"""
function get_N(net::Network ; sparse_n=false)
    (n1,n2,n3) = get_full_N(net ; sparse_n=sparse_n)
    return n1
end

"""
  get_full_N( n::Network ; sparse=false)
returns Tuple{ Array{Float64,2} , Array{Float64,2}  , Array{Float64,2} }
computes the stoichiometric matrix of the network result[1], and the positive / negative
entries of the stoichiometric matrix (result[2]/result[3])
"""
function get_full_N(net::Network ; sparse_n=false)

    n_plus::Array{Float64,2}  = zeros(length(net.species),length(net.fluxes))
    n_minus::Array{Float64,2} = zeros(length(net.species),length(net.fluxes))

    d_s::Dict{String,Int64}   = Dict{String,Int64}()
    for zi=1:length(net.species); d_s[net.species[zi]] = zi; end

    for zi=1:length(net.fluxes)
        fi = net.fluxes_db[ net.fluxes[zi] ]
        for ki in keys(fi.E); if(ki=="*"); continue; end ; n_minus[d_s[ki],zi] = fi.E[ki]; end
        for ki in keys(fi.P); if(ki=="*"); continue; end ; n_plus[d_s[ki],zi]  = fi.P[ki]; end
    end

    n = (n_plus - n_minus)

    if(!sparse_n)
        return n,n_plus,n_minus
    else
        n_sp = SparseArrays.sparse(n); n_sp_p = SparseArrays.sparse(n_plus); n_sp_m = SparseArrays.sparse(n_minus);
        return n_sp , n_sp_p , n_sp_m
    end
end

"""
is_sx_flux(f::Flux)
returns Bool

checks if flux is sx flux, which is by convention equivalent of having the
"*" id as a key in reactand lc or product lc
"""
function is_sx_flux(f::Flux)
    return haskey( f.E , "*" ) || haskey( f.P , "*" )
end

"""
get_sx_fluxes(net::Network)
returns Vector{Flux}
returns sx fluxes of the network
"""
function get_sx_fluxes(net::Network)
    return map( x -> net.fluxes[x] , get_idx_sx_fluxes(net) )
end
function get_idx_sx_fluxes(net::Network)
    fsx = Vector{Int64}()
    fi_cnt = 0
    for fi in net.fluxes
        fi_cnt += 1
        if( is_sx_flux( net.fluxes_db[fi] ) ) ; push!(fsx,fi_cnt) ; end
    end
    return fsx
end
function get_non_sx_fluxes(net::Network)
    return map( x -> net.fluxes[x] , get_idx_non_sx_fluxes(net) )
end
function get_idx_non_sx_fluxes(net::Network)
    fsx = Vector{Int64}()
    fi_cnt = 0
    for fi in net.fluxes
        fi_cnt += 1
        if( !is_sx_flux( net.fluxes_db[fi] ) ) ; push!(fsx,fi_cnt) ; end
    end
    return fsx
end

"""
get_flux_polyhedron(net::Network , constraints::Vector{LinConstr} ; max_val=Inf)
returns Polyhedron

creates the flux polyhedron
CONVENTION:
System exchange fluxes are reversible, all other fluxes have 0 as lower bound
"""
function get_flux_polyhedron(net::Network , constraints::Vector{LinConstr} ; max_val=Inf)
     p_0 = get_flux_polyhedron_vanilla(net::Network , constraints::Vector{LinConstr} )

     d_s::Dict{String,Int64}   = Dict{String,Int64}()
     for zi=1:length(net.fluxes); d_s[net.fluxes[zi]] = zi; end

     # identify non-system-exchange fluxes
     fnsx = get_idx_non_sx_fluxes(net)
     #print(fnsx)

     lb =  Inf * -ones(length(net.fluxes))
     ub =  Inf *  ones(length(net.fluxes))
     lb[fnsx] = 0

     # if max_val, create bouding box. [0,max_val] for nsx fluxes, AND [-max_val,max_val] for sx fluxes
     if(max_val<Inf)
         ub[ub.>max_val] = max_val
         lb[ get_idx_sx_fluxes(net) ] = -max_val
     end

     return Polyhedron( p_0 , bounding_box(lb,ub) )
end

"""
get_flux_polyhedron_vanilla(net::Network , constraints::Vector{LinConstr} )
returns Polyhedron

creates the flux polyhedron, without >0 constraint on non-sx-fluxes

This function does not add the greater than zero constraints to the non-system-
exchange fluxes
"""
function get_flux_polyhedron_vanilla(net::Network , constraints::Vector{LinConstr} )
    m_n = get_N(net)
    p_c = lc2polyhedron(net , constraints)
    return Polyhedron( p_c.G , p_c.h , [m_n;p_c.A] , [zeros(size(m_n,1),1);p_c.b]  )
end

"""
lc2polyhedron(net::Network , lcs::Vector{LinConstr})
returns Polyhedron
creates polyhedron from the given constraints
"""
function lc2polyhedron(net::Network , lcs::Vector{LinConstr} )
    n = length(net.fluxes)
    d_s::Dict{String,Int64}   = Dict{String,Int64}()
    for zi=1:length(net.fluxes); d_s[net.fluxes[zi]] = zi; end
    d_s["1"] = n+1 # to collect the "1" data

    G::Array{ Float64 , 2} = Array{ Float64 }(0,n)
    h::Array{ Float64 , 2} = Array{ Float64 }(0,1)
    A::Array{ Float64 , 2} = Array{ Float64 }(0,n)
    b::Array{ Float64 , 2} = Array{ Float64 }(0,1)

    v_G =  Array{ RowVector{Float64,Array{Float64,1}},1 }()
    v_A =  Array{ RowVector{Float64,Array{Float64,1}},1 }()
    v_h =  Array{ Float64 , 1 }()
    v_b =  Array{ Float64 , 1 }()

    G_s = spzeros(0,n)
    A_s = spzeros(0,n)

    for lci in lcs
        row_lhs = resolve_lc( n+1 , d_s , lci.lhs )
        row_lhs_one = row_lhs[n+1]
        row_lhs = row_lhs[1:n]
        row_rhs = resolve_lc( n+1 , d_s , lci.rhs )
        row_rhs_one = row_rhs[n+1]
        row_rhs = row_rhs[1:n]

        # copy all to rhs, except for "one" part:
        p_row_lhs = row_lhs - 1.0 * row_rhs
        p_row_rhs = row_rhs_one - row_lhs_one

        if(lci.constr=='=')
            #A = [A;p_row_lhs']; b = [b;p_row_rhs]; continue;
            #push!( v_A , p_row_lhs'); push!( v_b, p_row_rhs); continue;
            A_s = vcat( A_s , p_row_lhs' ); push!( v_b, p_row_rhs); continue;
        end
        if(lci.constr=='<')
            #G = [G;p_row_lhs']; h = [h;p_row_rhs]; continue;
            #push!( v_G , p_row_lhs'); push!(v_h,p_row_rhs); continue;
            G_s = vcat( G_s , p_row_lhs' ); push!(v_h,p_row_rhs); continue;
        end
        if(lci.constr=='>')
            #G = [G;-p_row_lhs']; h = [h;-p_row_rhs]; continue;
            #vG = push!(-p_row_lhs'); push!(v_h,-p_row_rhs); continue;
            G_s = vcat( G_s , -p_row_lhs' ); push!(v_h,-p_row_rhs); continue;
        end
    end
    #print("\nDONE!\n")
    #println("\n\n $(G) \n\n $(h) \n\n $(A) \n\n $(b) \n\n ")
    #G = reduce( (x,y) -> [x;y] , v_G )
    #A = reduce( (x,y) -> [x;y] , v_A )
    h = [ h ; v_h]
    b = [ b ; v_b]
    G = full(G_s)
    A = full(A_s)
    return Polyhedron(G,h,A,b)
end


text:
__precompile__()

module SpecialFunctions

using Compat

if VERSION >= v"0.6.0-dev.2767"
    if isdefined(Base, :airyai)
        import Base: airyai, airyaix, airyaiprime, airyaiprimex,
                     airybi, airybix, airybiprime, airybiprimex,
                     besselh, besselhx, besseli, besselix, besselj, besselj0, besselj1,
                     besseljx, besselk, besselkx, bessely, bessely0, bessely1, besselyx,
                     hankelh1, hankelh1x, hankelh2, hankelh2x,
                     dawson, erf, erfc, erfcinv, erfcx, erfi, erfinv,
                     eta, digamma, invdigamma, polygamma, trigamma, zeta,
                     # deprecated
                     airy, airyx, airyprime
    else
        export
            airyai,
            airyaiprime,
            airybi,
            airybiprime,
            airyaix,
            airyaiprimex,
            airybix,
            airybiprimex,
            besselh,
            besselhx,
            besseli,
            besselix,
            besselj,
            besselj0,
            besselj1,
            besseljx,
            besselk,
            besselkx,
            bessely,
            bessely0,
            bessely1,
            besselyx,
            dawson,
            erf,
            erfc,
            erfcinv,
            erfcx,
            erfi,
            erfinv,
            eta,
            digamma,
            invdigamma,
            polygamma,
            trigamma,
            hankelh1,
            hankelh1x,
            hankelh2,
            hankelh2x,
            zeta
    end
end

if isdefined(Base.Math, :openspecfun)
    const openspecfun = Base.Math.openspecfun
else
    const openspecfun = "libopenspecfun"
end

include("bessel.jl")
include("erf.jl")
include("gamma.jl")
include("deprecated.jl")

end # module


text:
function applytourniquet!(system::CVSystem,n::Int64)
    if system.hemo.Vloss >= system.hemo.totalloss
        for i = 1:length(system.branches.ID)
            tb = system.branches.beta[i][end];
            ta = system.branches.A0[i][end];
            tc = system.branches.c0[i][end]
            for j = 1:length(system.hemo.tID)
                if i == system.hemo.tID[j]
                    if system.branches.beta[i][end] < system.hemo.bmax[j]
                        tb = (system.branches.beta[i][end]+
                            (system.hemo.bmax[j]-system.hemo.bmin[j])/system.hemo.ttn*
                            system.solverparams.h);
                    end
                    if system.branches.A0[i][end] > system.hemo.Amin[j]
                        ta = (system.branches.A0[i][end]+
                            (system.hemo.Amin[j]-system.hemo.Amax[j])/system.hemo.ttn*
                            system.solverparams.h);
                        tc = sqrt(system.branches.beta[i][end]/
                            (2*system.solverparams.rho))*system.branches.A0[i][end]^0.25;
                    end
                end
            end
            push!(system.branches.beta[i],tb)
            push!(system.branches.A0[i],ta);
            push!(system.branches.c0[i],tc);
        end
    end
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
"""
	fcr, geo = facecurrents(coeffs, basis)

Compute the value of the function with the given collection of coeffient in the provided basis in all the centroids of the mesh underlying the basis. The mesh is returned together with the currents.
"""
function facecurrents(coeffs, basis)

	T = eltype(coeffs)
	RT = real(T)

	refs = refspace(basis)
	numrefs = numfunctions(refs)

	cells, tad, a2g = assemblydata(basis)

	mesh = geometry(basis)
	D = dimension(mesh)
	# U = D+1
	U = 3

	# TODO: remove ugliness
	vals = refs(center(first(cells)))
	PT = typeof(first(coeffs)*vals[1][1])

	fcr = zeros(PT, numcells(mesh))

	for (t,cell) in enumerate(cells)

		mps = neighborhood(cell, ones(RT,D)/(D+1))
        vals = refs(mps)

		# assemble in the right hand side vector
		for i in 1 : numrefs
			fx = vals[i][1]
            for (m,a) in tad[t,i]
				# fcr[t] += coeffs[m] * a * fx
				fcr[a2g[t]] += coeffs[m] * a * fx
			end
		end
	end

	return fcr, geometry(basis)
end

function facecurrents(coeffs, basis::SpaceTimeBasis)

	space_basis = basis.space
	time_basis = basis.time

	Nt = numfunctions(time_basis)
	Δt = timestep(time_basis)

	refs = refspace(space_basis)
	trefs = refspace(time_basis)
	numrefs = numfunctions(refs)
	tnumrefs = numfunctions(trefs)

	cells, ad = assemblydata(space_basis)
	tcells, tad = assemblydata(time_basis)

	mesh = geometry(space_basis)
	T = eltype(coeffs)
	D = dimension(mesh)
	U = D+1

	# TODO: express relative to input types
	PT = SVector{U, T}
	fcr = zeros(PT, numcells(mesh), Nt)

	for (k,tcell) in enumerate(tcells)
		tmps = neighborhood(tcell,1)
		tvals = trefs(tmps)
		for (p,cell) in enumerate(cells)
			mps = center(cell)
	        vals = refs(mps)

			# assemble
			for i in 1:numrefs
				fx = vals[i][1]
				for j in 1:tnumrefs
					tfx = tvals[j]
	    			for (m,a) in ad[p,i]
						for (n,b) in tad[k,j]
							fcr[p,k] += (coeffs[m,n] * tfx * b) * a * fx
						end
					end
				end
			end

		end
	end

	return fcr, geometry(space_basis)
end

function facecurrents(u, X1, Xs...)

    offset = 0

    n = numfunctions(X1)
    fcrs = [ facecurrents(u[offset + (1:n)], X1)[1] ]
    offset += n

    for i in 1:length(Xs)
        n = numfunctions(Xs[i])
        push!(fcrs, facecurrents(u[offset + (1:n)], Xs[i])[1])
        offset += n
    end

    fcrs
end

function facecurrents(u, X::DirectProductSpace)

	if length(X.factors) == 1
		return facecurrents(u, X.factors[1]), geometry(X.factors[1])[1]
	end

	fcrs = facecurrents(u, X.factors...)
	#fcr = append!(fcrs...)
	fcr = vcat(fcrs...)
	m = weld([geometry(x) for x in X.factors]...)

	fcr, m
end

function potential(op, points, coeffs, basis; type=SVector{3,ComplexF64})
	# T = SVector{3,ComplexF64}
	T = type
	ff = zeros(T, size(points))
	store(v,m,n) = (ff[m] += v*coeffs[n])
	potential!(store, op, points, basis, type=T)
	return ff
end

function potential(op, points,coeffs, basis::SpaceTimeBasis)
	T = SVector{3,ComplexF64}
	ff = zeros(T, length(points), size(coeffs)[2])
	store(v,m,n,k,o) = (ff[m,k] += v*coeffs[n,o])
	potential!(store, op, points, basis)
	return ff
end

function potential(op, points, coeffs, space::DirectProductSpace; type=SVector{3,ComplexF64})
	# T = SVector{3,ComplexF64}
	T = type
	ff = zeros(T, size(points))
	@show size(ff)

	@assert length(coeffs) == numfunctions(space)

	offset = 0
	for fct in space.factors
		store(v,m,n) = (ff[m] += v*coeffs[offset+n])
		potential!(store, op, points, fct, type=T)
		offset += numfunctions(fct)
	end

	ff
end

function potential!(store, op, points, basis; type=SVector{3,ComplexF64})

	# T = SVector{3,ComplexF64}
	T = type
	# @show T
	z = zeros(T,length(points))

	els, ad = assemblydata(basis)
	rs = refspace(basis)

	zlocal = Array{T}(undef,numfunctions(rs))
	qdata = quaddata(op,rs,els)

	#println("Computing nearfield.")
	print("dots out of 10: ")

	todo, done, pctg = length(points), 0, 0

	for (p,y) in enumerate(points)
		for (q,el) in enumerate(els)

			fill!(zlocal,zero(T))
			qr = quadrule(op,rs,p,y,q,el,qdata)
			farfieldlocal!(zlocal,op,rs,y,el,qr)

			# assemble from local contributions
			for (r,z) in enumerate(zlocal)
                for (n,b) in ad[q,r]
					store(z*b,p,n)
				end
			end
		end

		done += 1
		new_pctg = round(Int, done / todo * 100)
		if new_pctg > pctg + 9
				#println(todo," ",done," ",new_pctg)
				print(".")
				pctg = new_pctg
		end
	end

	println("")

end

function potential!(store, op, points, basis::SpaceTimeBasis)

    space_basis = basis.space
	tb = time_basis = basis.time

	Nt = numfunctions(time_basis)
	Δt = timestep(time_basis)

	T = SVector{3,ComplexF64}

	refs = refspace(space_basis)
	trefs = refspace(time_basis)

	numrefs = numfunctions(refs)
	tnumrefs = numfunctions(trefs)

	els, ad = assemblydata(space_basis)
	tels, tad = assemblydata(time_basis)

	zlocal = Array{T}(numrefs,Nt)
	qdata = quaddata(op,refs,trefs,els,tels)

	#println("Computing nearfield.")
	print("dots out of 10: ")

	todo, done, pctg = Nt, 0, 0
    for (k,tel) in enumerate(tels)
    	for (p,y) in enumerate(points)
    		for (q,el) in enumerate(els)

    			fill!(zlocal,zero(T))
    			qr = quadrule(op,refs,trefs,p,y,q,el,k,tel,qdata)
	    		farfieldlocal!(zlocal,op,refs,trefs,p,y,q,el,k,tel,tb,qr)

    			# assemble from local contributions
    			 for r in 1:numrefs
					 for s in 1:Nt
	                    for (n,b) in ad[q,r]
	        				store(zlocal[r,s]*b,p,n,k,s)
							# end
	    				end
					end
    			end
            end
    	end

		done += 1
		new_pctg = round(Int, done / todo * 100)
		if new_pctg > pctg + 9
				#println(todo," ",done," ",new_pctg)
				print(".")
				pctg = new_pctg
		end
    end

	println("")

end

function farfieldlocal!(zlocal,op,refspace,y,el,qr)

    for q in qr
        x = q.point
        F = q.value
        dx = q.weight

        krn = kernelvals(op, y, x)
        for r in 1 : length(zlocal)
            zlocal[r] += integrand(op,krn,y,F[r],x) * dx
        end

    end

end

function farfieldlocal!(zlocal,op,trialrefs, timerefs,
        p, testel, q, trialel, k, timeel, tb,quadrule)

	M,N = size(zlocal)
	Δt = tb.timestep

    for qr in quadrule[1]
        x = qr.point
        F = qr.value
        dx = qr.weight

        krn = kernelvals(op, testel, x)

        for r = 1:M
			for l = 1:N
				t = (k - l)*Δt
            	zlocal[r,l] += integrand(op,krn,testel,F[r], t, tb) * dx
			end
		end

    end
end


text:
"""
Function fits an unsupervised model. Will be deprecated in next version.
"""
function fit_unsupervised!(M::Model, X, k; beta=0.8, alpha=0.5, init=nothing, t_init=1.0, max_iters=5000, verbose=false, tol=1e-4)
    n, d = size(X)
    if M.intercept
        X = [X ones(n)]
        d += 1
    end

    M.param["theta"] = Dict("X"=> zeros(n, k), "Y"=>zeros(d, k))

    opt = EmpiricalRiskMinimization.minimize_unsupervised(M.loss, M.reg, X, k, beta, alpha, init, t_init, max_iters=max_iters, verbose=verbose, tol=tol)
    if opt == -1
        println("Model did not converge")
        M.status = "Failed"
        return
    end

    theta_X, theta_Y, losses = opt

    M.param["theta_history"] = Dict("X"=> theta_X, "Y"=>theta_Y)
    M.param["theta"] = Dict("X"=>M.param["theta_history"]["X"][end],
                            "Y"=>M.param["theta_history"]["Y"][end])
    M.param["loss_history"] = losses

    M.status = "Converged"
end


text:
using Clang.Generators
using Yosys_jll

hapi =

cd(@__DIR__)

include_dir = joinpath(Yosys_jll.find_artifact_dir(), "share/yosys/include") |> normpath

options = load_options(joinpath(@__DIR__, "generator.toml"))

# add compiler flags, e.g. "-DXXXXXXXXX"
args = get_default_args()
push!(args, "-I$include_dir")

headers = [joinpath(include_dir, "backends/cxxrtl/cxxrtl_capi.h")]

# create context
ctx = create_context(headers, args, options)

# run generator
build!(ctx)

max_stars_repo_name: JuliaComputing/CxxrtlCosim

max_issues_repo_name: JuliaComputing/CxxrtlCosim

max_forks_repo_name: JuliaComputing/CxxrtlCosim


text:
using Test, Pkg
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    using CuArrays
    ArrayType = CuArray
else
    ArrayType = Array
end

MPI.Init()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

dst = mod(rank+1, size)
src = mod(rank-1, size)

N = 32

send_mesg = ArrayType{Float64}(undef,N)
recv_mesg = ArrayType{Float64}(undef,N)
recv_mesg_expected = ArrayType{Float64}(undef,N)
fill!(send_mesg, Float64(rank))
fill!(recv_mesg_expected, Float64(src))

rreq = MPI.Irecv!(recv_mesg, src,  src+32, comm)
sreq = MPI.Isend(send_mesg, dst, rank+32, comm)

reqs = [sreq,rreq]

(inds,stats) = MPI.Waitsome!(reqs)
@test !isempty(inds)
for ind in inds
    (onedone,stat) = MPI.Test!(reqs[ind])
    @test onedone
    @test stat == MPI.STATUS_EMPTY
end

(done, ind, stats) = MPI.Testany!(reqs)
if done && ind != 0
    (onedone,stat) = MPI.Test!(reqs[ind])
    @test onedone
    @test stat == MPI.STATUS_EMPTY
end

MPI.Waitall!(reqs)

(inds, stats) = MPI.Waitsome!(reqs)
@test isempty(inds)
@test isempty(stats)
(inds, stats) = MPI.Testsome!(reqs)
@test isempty(inds)
@test isempty(stats)

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


text:
type StreamingStats{T<:Union{Number,Vector}}
    iteration::Int
    runningmean::T
    runningvar::T
    tempdiff::T
end

StreamingStats{T<:AbstractFloat}(::Type{T}) = StreamingStats(0,zero(T),zero(T),zero(T))
StreamingStats{T<:AbstractFloat}(s::Integer,::Type{T}) = StreamingStats(0,zeros(T,s),zeros(T,s),zeros(T,s))

@inline _varupdatefactor(i::Integer,v::Number) = i==1?zero(typeof(v)):convert(typeof(v),(i-2)/(i-1))
@inline _varupdatefactor(i::Integer,v::Vector) = i==1?zero(eltype(v)):convert(eltype(v),(i-2)/(i-1))

@inline _updatetempdiff!(s::StreamingStats,v::Number) = (s.tempdiff = v - s.runningmean ; s)
@inline _updatetempdiff!(s::StreamingStats,v::Vector) = (@simd for i=1:length(v) @inbounds s.tempdiff[i] = v[i] - s.runningmean[i] end ; s)

@inline _updatemean!(s::StreamingStats,v::Number) = (s.runningmean = ((s.iteration-1)*s.runningmean + v)/s.iteration ; s)
@inline _updatemean!{T<:AbstractFloat}(s::StreamingStats,v::Vector{T}) = (@simd for i=1:length(v) @inbounds s.runningmean[i] = ((s.iteration-1)*s.runningmean[i] + v[i])/s.iteration end ; s)

@inline _updatevar!{T<:AbstractFloat}(s::StreamingStats,v::T,f::T) = (s.runningvar = f*s.runningvar + s.tempdiff*s.tempdiff/s.iteration ; s)
@inline _updatevar!{T<:AbstractFloat}(s::StreamingStats,v::Vector{T},f::T) = (@simd for i=1:length(v) @inbounds s.runningvar[i] = f*s.runningvar[i] + s.tempdiff[i]*s.tempdiff[i]/s.iteration end ; s)

function update!(s::StreamingStats,value)
    s.iteration += 1
    _updatetempdiff!(s,value)
    _updatemean!(s,value)
    _updatevar!(s,value,_varupdatefactor(s.iteration,value))
    s
end

function streamingstats(x::Vector)
    s = StreamingStats(eltype(x))
    m = zeros(x)
    v = zeros(x)
    for i=1:length(x)
        update!(s,x[i])
        m[i] = mean(s)
        v[i] = var(s)
    end
    m,v
end

function streamingstats(a::Matrix)
    s = StreamingStats(size(a,1),eltype(a))
    m = zeros(a)
    v = zeros(a)
    for i=1:size(a,2)
        update!(s,a[:,i])
        m[:,i] = mean(s)
        v[:,i] = var(s)
    end
    m,v
end

mean(s::StreamingStats) = s.runningmean
var(s::StreamingStats) = s.runningvar*(s.iteration-1)/(s.iteration-1)
std(s::StreamingStats) = sqrt(var(s))

function show(io::IO,s::StreamingStats)
    println(io,"iterations = ",s.iteration)
    println(io,"mean = ",round(mean(s),2))
    println(io,"var = ",round(var(s),2))
    println(io,"std = ",round(std(s),2))
end


text:
using EconModels, Base.Test

include("settings.jl")
include("model.jl")


text:
struct StructReinterpretArray{StructType,N,ArrayType<:AbstractArray} <: AbstractArray{StructType,N}
    parent::ArrayType
    function StructReinterpretArray(::Type{ST},data::AbstractArray{T,N}) where {ST,T,N}
        NF = check_parent_and_struct_match(ST,T)
        dim = size(data,1)
        dim % NF == 0 || throw(ArgumentError("""
        cannot reinterpret an `$(T)` array to `$(ST)` whose first dimension has size `$(dim)`.
        The resulting array would have non-integral first dimension.
        """))
        return new{ST,N,typeof(data)}(data)
    end
end

struct StructReinterpretDenseArray{StructType,N,ArrayType<:DenseArray} <: DenseArray{StructType,N}
    parent::ArrayType
    function StructReinterpretDenseArray(::Type{ST},data::DenseArray{T,N}) where {ST,T,N}
        NF = check_parent_and_struct_match(ST,T)
        dim = size(data,1)
        dim % NF == 0 || throw(ArgumentError("""
        cannot reinterpret an `$(T)` array to `$(ST)` whose first dimension has size `$(dim)`.
        The resulting array would have non-integral first dimension.
        """))
        return new{ST,N,typeof(data)}(data)
    end
end

const AbstractReinterpretArray{ST,N,AT} = Union{<:StructReinterpretArray{ST,N,AT},<:StructReinterpretDenseArray{ST,N,AT}}


text:
using Core.Compiler: CodeInfo, NewvarNode, GotoNode
using Base.Meta: isexpr
using JuliaInterpreter: is_global_ref

const SSAValues = Union{Core.Compiler.SSAValue, JuliaInterpreter.SSAValue}

const structheads = VERSION >= v"1.5.0-DEV.702" ? () : (:struct_type, :abstract_type, :primitive_type)
const trackedheads = (:method, structheads...)

isssa(stmt) = isa(stmt, Core.Compiler.SSAValue) | isa(stmt, JuliaInterpreter.SSAValue)
isslotnum(stmt) = isa(stmt, Core.Compiler.SlotNumber) | isa(stmt, JuliaInterpreter.SlotNumber)

"""
`SlotDep` is an internal type used for holding information about dependencies of a SlotNumber
variable. If `sd` is a `SlotDep`, then
- `sd.lineassigned` is the statement number on which the SlotNumber was most recently assigned
- `sd.linedeps` is the set of previous SSAValues upon which this assigment depends (encoded as Int)
"""
mutable struct SlotDep
    lineassigned::Int
    linedeps::Vector{Int}  # we use Int to handle Core.Compiler.SSAValue & JuliaInterpreter.SSAValue
end
function SlotDep(i::Int, stmt, slotdeps)
    deps = add_deps!(Int[], stmt, slotdeps)
    SlotDep(i, deps)
end
function add_deps!(linedeps, stmt, slotdeps)
    if isssa(stmt)
        push!(linedeps, stmt.id)
    elseif isslotnum(stmt)
        append!(linedeps, slotdeps[stmt.id].linedeps)
        push!(linedeps, slotdeps[stmt.id].lineassigned)
    elseif isa(stmt, Expr)
        for a in stmt.args
            if isssa(a) || isslotnum(a) || isa(a, Expr)
                add_deps!(linedeps, a, slotdeps)
            end
        end
    end
    return linedeps
end

# See docs below for `BackEdges(ci::CodeInfo)`
mutable struct Uses
    ids::Vector{Int}        # lines where it is used
    defined::Vector{Int}    # lines where it is defined
end
Uses() = Uses(Int[], Int[])
Base.push!(uses::Uses, loc::Int) = push!(uses.ids, loc)

struct BackEdges
    byid::Vector{Vector{Int}}
    byname::Dict{Union{GlobalRef,Symbol},Uses}
end
BackEdges(n::Integer) = BackEdges([Int[] for i = 1:n], Dict{Union{GlobalRef,Symbol},Uses}())

id_loc(pr::Pair) = pr.first, pr.second
id_loc(pr::Pair{<:SSAValues,Int}) = pr.first.id, pr.second
id_loc(pr::Pair{Bool,Int}) = pr.second - 1, pr.second # handles statements like `if true f() = 1 else f() = 2 end`
function Base.push!(be::BackEdges, pr::Pair{<:Union{Integer,SSAValues},Int})
    @noinline errorder(y, x) = throw(ArgumentError("SSA form requires that dependencies come after the statement, got $y and $x"))

    id, loc = id_loc(pr)
    loc > id || errorder(loc, id)
    push!(be.byid[id], loc)
    return be
end
function Base.push!(be::BackEdges, pr::Pair{named,Int}) where named <: Union{Symbol,GlobalRef}
    id, loc = pr
    if isa(id, GlobalRef)
        obj = getfield(id.mod, id.name)
        isa(obj, Core.Builtin) && return be
    end
    ref = get(be.byname, id, nothing)
    if ref === nothing
        be.byname[id] = ref = Uses()
    end
    push!(ref, loc)
    return be
end

function add_to_backedges!(backedges::BackEdges, slotdeps, loc, stmt)
    if isssa(stmt)
        push!(backedges, stmt=>loc)
    elseif isslotnum(stmt)
        isassigned(slotdeps, stmt.id) || return backedges # issue #428 (it's inside an `if false...end`)
        sd = slotdeps[stmt.id]
        if sd.lineassigned != 0
            push!(backedges, sd.lineassigned=>loc)
        end
        for id in sd.linedeps
            push!(backedges, id=>loc)
        end
    elseif stmt isa GlobalRef
        push!(backedges, stmt=>loc)
    elseif stmt isa Symbol
        push!(backedges, stmt=>loc)
    elseif stmt isa Expr
        head = stmt.head
        if head === :call && !(isssa(stmt.args[1]) || isslotnum(stmt.args[1]))
            for a in Iterators.drop(stmt.args, 1)  # don't track the callee
                add_to_backedges!(backedges, slotdeps, loc, a)
            end
        elseif head ∈ trackedheads && (name = stmt.args[1]; isa(name, Symbol) | isa(name, GlobalRef))
            push!(backedges, name=>loc)
            push!(backedges.byname[name].defined, loc)
            for a in Iterators.drop(stmt.args, 1)
                add_to_backedges!(backedges, slotdeps, loc, a)
            end
        elseif head === :(=) && (lhs = stmt.args[1]; isa(lhs, Symbol) || isa(lhs, GlobalRef))
            add_to_backedges!(backedges, slotdeps, loc, lhs)
            push!(backedges.byname[lhs].defined, loc)
            add_to_backedges!(backedges, slotdeps, loc, stmt.args[2])
        else
            for a in stmt.args
                add_to_backedges!(backedges, slotdeps, loc, a)
            end
        end
    elseif stmt isa Pair
        add_to_backedges!(backedges, slotdeps, loc, stmt.first)
        add_to_backedges!(backedges, slotdeps, loc, stmt.second)
    elseif stmt isa Tuple
        for item in stmt
            add_to_backedges!(backedges, slotdeps, loc, item)
        end
    # elseif isa(stmt, QuoteNode) || isa(stmt, NewvarNode) || isa(stmt, GotoNode) ||
    #        isa(stmt, Real) || isa(stmt, CodeInfo) || isa(stmt, Nothing) ||
    #        isa(stmt, Module) || isa(stmt, String) || isa(stmt, Char) || isa(stmt, Type) ||
    #        isa(stmt, LineNumberNode)
    # else
    #     error("unhandled stmt ", stmt, " of type ", typeof(stmt), " at ", loc)
    end
    return backedges
end

function toplevel_blocks(bbs::Core.Compiler.CFG)
    istoplevel = falses(length(bbs.blocks))
    next = 1
    for (i, block) in enumerate(bbs.blocks)
        if i == 1
            istoplevel[i] = true
        elseif i < next
        else
            istoplevel[i] = sum(block.preds .< i) == 2
        end
        if istoplevel[i] && !isempty(block.succs)
            next = maximum(block.succs)
        end
    end
    return istoplevel
end

function add_block_dependents!(backedges::BackEdges, bbs, istoplevel, i, bbidx)
    for s in bbs.blocks[bbidx].succs
        s > bbidx || continue  # follow only in the forward direction
        istoplevel[s] && continue
        r = bbs.blocks[s].stmts
        for j = Core.Compiler.first(r):Core.Compiler.last(r)
            push!(backedges, i=>j)
        end
        add_block_dependents!(backedges, bbs, istoplevel, i, s)
    end
    return backedges
end

"""
    backedges = BackEdges(code::CodeInfo)

Analyze `code` and determine the chain of dependencies.
`backedges.byid[i]` lists the succeeding lines that depend on `i`.
(In addition to SSAValue dependencies, this includes basic-block control flow
dependencies.)
`backedges.byname[sym]` lists the lines that depend on a particular symbol.
"""
function BackEdges(ci::CodeInfo)
    ci.inferred && error("supply lowered but not inferred code")
    bbs = Core.Compiler.compute_basic_blocks(ci.code)
    istoplevel = toplevel_blocks(bbs)
    codelocs = ci.codelocs
    n = length(ci.code)
    backedges = BackEdges(n)
    slotdeps = Vector{SlotDep}(undef, length(ci.slotnames))
    slotassign = zeros(Int, length(ci.slotnames))
    i, n = 1, length(ci.code)
    while i < n
        stmt = ci.code[i]
        if isa(stmt, Expr)
            head = stmt.head
            if head === :(=) && isslotnum(stmt.args[1])
                id = stmt.args[1].id
                slotdeps[id] = SlotDep(i, stmt.args[2], slotdeps)
                i += 1
            elseif head === :gotoifnot
                dep, _ = stmt.args
                add_to_backedges!(backedges, slotdeps, i, dep)
                # Add the non-toplevel successor basic blocks as dependents of this line
                bbidx = searchsortedlast(bbs.index, i) + 1 # bb index of this line
                add_block_dependents!(backedges, bbs, istoplevel, i, bbidx)
                i += 1
            elseif head === :thunk
                ci_thunk = stmt.args[1]
                if isa(ci_thunk, CodeInfo)
                    bethunk = BackEdges(ci_thunk)
                    for (name, uses) in bethunk.byname
                        add_to_backedges!(backedges, slotdeps, i, name)
                        if !isempty(uses.defined)
                            push!(backedges.byname[name].defined, i)
                        end
                    end
                end
                i += 1
            else
                add_to_backedges!(backedges, slotdeps, i, stmt)
                i += 1
            end
        else
            add_to_backedges!(backedges, slotdeps, i, stmt)
            i += 1
        end
    end
    return backedges
end

## Now that we can construct BackEdges, let's use them for analysis

function toplevel_chunks(backedges::BackEdges)
    be = [copy(usedby) for usedby in backedges.byid]
    n = length(be)
    for (name, uses) in backedges.byname
        isempty(uses.defined) && continue
        id = minimum(uses.ids)
        append!(be[id], uses.ids)
    end
    chunkid = collect(1:n)
    for i = n:-1:1
        if !isempty(be[i])
            chunkid[i] = maximum(chunkid[be[i]])
        end
    end
    chunks = UnitRange{Int}[]
    i = 1
    while i <= n
        push!(chunks, i:chunkid[i])
        i = chunkid[i]+1
    end
    return chunks
end

"""
    hastrackedexpr(src, chunk=1:length(src.code))

Detect whether the specified `chunk` of `src` contains a definition tracked
by Revise. By default these are methods and type definitions (the latter because
they might contain constructors).
"""
function hastrackedexpr(code::CodeInfo, chunk::AbstractUnitRange=axes(code.code, 1); heads=trackedheads)
    for stmtidx in chunk
        stmt = code.code[stmtidx]
        if isa(stmt, Expr)
            stmt.head == :call && is_global_ref(stmt.args[1], Core, :_typebody!) && return true
            stmt.head ∈ heads && return true
            if stmt.head == :thunk
                hastrackedexpr(stmt.args[1]; heads=heads) && return true
            end
        end
    end
    return false
end


text:
using YaoLang
using YaoLang.IBMQ
using IBMQClient
using YaoCompiler
using YaoCompiler.Intrinsics

@device function test(theta)
    1 => X
    2 => X
    @ctrl (1, 2) 3 => X
    @ctrl 1 2 => Rx(theta)
    c = @measure 1:3
    return c
end

# token = "e773394070269e3deace4372ed915c99610ee5a0e3be7b2f821e6889f4f4fe93cafdebcce46009e0ec9dd3ff8dca3ad3eb126b3bc59617ee837f2e120e99f268"

job = @ibmq nshots = 2048 begin
    for i in 1:10
        test()
    end
end


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.1.32/11.3_06072021/cudnn-11.3-linux-x64-v8.2.1.32.tgz",
                      "39412acd9ef5dd27954b6b9f5df75bd381c5d7ceb7979af6c743a7f4521f9c77")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.1.32/11.3_06072021/cudnn-11.3-linux-ppc64le-v8.2.1.32.tgz",
                      "4ee4f2afeaae34fdb06da8d4942a6802aae94ecc51f307292c45966eecbe5fb9")],
    Platform("aarch64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.1.32/11.3_06072021/cudnn-11.3-linux-aarch64sbsa-v8.2.1.32.tgz",
                      "e3a0e570cb8ba01d5d45e6eb1ebe29ff22fd5fb8ad45bfe7a448f4f95065ec1e")],
    Platform("x86_64", "windows") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/cudnn/secure/8.2.1.32/11.3_06072021/cudnn-11.3-windows-x64-v8.2.1.32.zip",
                      "5b9bf2dc4670fb1519ef55e13da5123f0b6b39fac5e6138e31388b269808d5f2")],
)


text:
export fftderivative, loc_max_amp, create_corr, get_corrs, mov_avg, get_node_files,
        plot_NS_dist, plot_tseries, plot_simple_nodes, wrap_plotter

# extracting and processing correlations
function fftderivative(A::CorrData)
    """ Computes discrete derivative via FFT"""
    C = deepcopy(A)
    detrend!(C)
    taper!(C)
    F = fft(C.corr,1)
    F .*= fftfreq(length(C.corr)).* 1im .* 2π
    C.corr = real.(ifft(F,1))/(C.fs^2)
    return C
end
# gets location and maximum amplitude for corrs in array
function loc_max_amp(corr::CorrData, filt::Array{Float64,1}=[0.1,0.5])
    """ Computes essential correlation metrics:
    latitude, longitude, max amplitude, max amplitude for anticausal side,
     max amplitude for causal side, SNR (signal-to-noise ratio), SNR (anticausal side),
        SNR (causal), correlation name
    """
    lat, lon = corr.loc.lat, corr.loc.lon
    bpass_corr = bandpass(fftderivative(fftderivative(corr)),filt[1],filt[2])
    half_len = convert(Int64,round(length(bpass_corr.corr)/2))
    max_ampl = maximum(abs.(bpass_corr.corr))
    max_ampl_anti = maximum(abs.(bpass_corr.corr[1:half_len]))
    max_ampl_causal = maximum(abs.(bpass_corr.corr[half_len:end]))
    snr = maximum(abs.(bpass_corr.corr))/std(bpass_corr.corr)
    snr_anti = maximum(abs.(bpass_corr.corr[1:half_len]))/std(bpass_corr.corr[1:half_len])
    snr_causal = maximum(abs.(bpass_corr.corr[half_len:end]))/std(bpass_corr.corr[half_len:end])
    return (lat, lon, max_ampl, max_ampl_anti, max_ampl_causal, snr, snr_anti, snr_causal, corr.name)
end
function create_corr(file::HDF5File, station::String, stacktype::String, component::String)
    """ Returns correlation with metadata from stacked HDF5 file"""
    try
        S = CorrData()
        # get corr data
        corr = file[station][component][stacktype][1:end]
        S.corr = reshape(corr,size(corr,1),1)
        S.name = station
        S.dist = read(file[station]["meta"]["dist"])
        S.loc.lat = read(file[station]["meta"]["lat"])
        S.loc.lon = read(file[station]["meta"]["lon"])
        S.loc.el = read(file[station]["meta"]["el"])
        S.azi = read(file[station]["meta"]["azi"])
        S.baz = read(file[station]["meta"]["baz"])
        S.fs = 20. # For the 20 Hz product
        return S
    catch
        return nothing
    end
end
# get information for a set of filtered stations, default to all stations
function get_corrs(file::HDF5File, stacktype::String="linear",
    component::String="ZZ", filter::String="")
    """Returns array of correlations from specified filters"""
    filtered = [name for name in names(file) if occursin(filter, name)]
    corrs = [create_corr(file, station, stacktype, component) for station in filtered]
    corrs = [corr for corr in corrs if !isnothing(corr)] # strip nothings - catches meta file etc which errors
    return corrs
end

# plotting stuff
function mov_avg(df::DataFrame, range::Int64=30)
    """ Returns dataframe with `moving_avg` column box function filter"""
    df_new = df[:,["date","max_amp"]]
    amps = df_new.max_amp
    avg, len = similar(amps), length(amps)
    for i =1:len
        first = max(1, convert(Int64, round(i-range/2)))
        last = min(len, convert(Int64, round(i+range/2)))
        window = amps[first:last]
        sd, mn = std(window), mean(window)
        avg[i] = mean(filter(x -> (x < mn + 3*sd) && (x > mn-3*sd), window))
    end
    df_new.moving_avg = avg
    return df_new
end
function get_node_files(source::String, deployment_filt::String, int_range=Array{Int64, 1})
    return filter(x -> occursin(source, x) & occursin(deployment_filt,x) & (parse(Int64, split(x, ".")[end-2][end-1:end]) in collect(int_range[1]:int_range[2])), nodes)
end

function plot_NS_dist(corrs::Array{CorrData,1}, freqs::Array{Array{Float64,1},1}, plot_sorttype::String,
    attr::Array{String,1}=["","","",""])
    lw, scale = 0.5, 3
    sort_method, p = [], []
    if plot_sorttype == "Latitude"
    sort_method = sort([corr.loc.lat for corr in corrs])
    p = sortperm([corr.loc.lat for corr in corrs]) # sort based on latitude
    elseif plot_sorttype == "Longitude"
    sort_method = sort([corr.loc.lon for corr in corrs])
    p = sortperm([corr.loc.lon for corr in corrs]) # sort based on longitude
    else # sort by distance
    plot_sorttype == "distance"
    println("Distance is not currently a supported stacktype")
    end

    lmin, lmax = minimum(sort_method), maximum(sort_method)
    new_corrs, n_corrs = corrs[p], length(corrs)
    y_labels = (round.(collect(range(0, 5*n_corrs, length=10)),digits=3),round.(collect(range(lmin, lmax, length=10)),digits=3))
    for freq_pair in freqs
        NS_plot = plot(xlims = (0,100), ylims = (-5,n_corrs*5+5),
                yticks=y_labels,xlabel = "Time (s)", ylabel = "Station Nodes by $plot_sorttype",
                xtickfontsize=5,ytickfontsize=5,fontsize=5, xguidefontsize = 10, yguidefontsize = 10,
                legendfontsize = 10, title="$(attr[3]) $(attr[2]) from $(attr[1]): $(freq_pair[1]) to $(freq_pair[2])",
                titlefontsize= 10)
        # bandpass correlations
        corrs_processed, max_amp = [], []
        for corr in new_corrs
        # bandpass, scale, shift, and add to plot
        processed = bandpass(corr, freq_pair[1], freq_pair[2]).corr ./ scale
        push!(max_amp, maximum(processed[1:convert(Int64, round(length(processed)/2))]))
        if plot_sorttype == "Latitude"
            push!(corrs_processed, processed .+ ((corr.loc.lat -lmin)/(lmax-lmin)*n_corrs*5))
        else
            push!(corrs_processed, processed .+ ((corr.loc.lon -lmin)/(lmax-lmin)*n_corrs*5))
        end
    end

    # plotting
    T = collect(-300.:1/new_corrs[1].fs:300.)
    plot!(NS_plot, -T, corrs_processed, color=:cm_maxamp, colorbar_title="Normalized Maximum Amplitude",
        line_z=max_amp', fmt = :png, linewidth = 0.5, reuse = false, legend = false)
    plot!(size=(250,400),dpi=500)
    filepath = joinpath(rootdir,
    "stack_plots/$(n_deriv_names[n_derivatives+1])_$(attr[4])/$(attr[1])/$(attr[3])_from_$(attr[1])_$(attr[2])_$(freq_pair[1])to$(freq_pair[2]).png")
    # ensure filepath is valid
    DIR = dirname(filepath)
    if !isdir(DIR)
    mkpath(DIR)
    end
    png(NS_plot,filepath)
    end
end

function plot_tseries(perm_station::String, node_stations_::Array{String,1}, source_, dir::String)
    test = DataFrame(CSV.File(perm_station))
    #get moving averages and naming schemes
    times, amplitudes = test.date, test.max_amp
    times_avg, amps_avg = mov_avg(test,90).date, mov_avg(test,90).moving_avg
    perm_name = join(split(split(perm_station,"/")[end],".")[3:4],".")
    node_range = join([split(node_stations_[1],".")[end-2],split(node_stations_[end],".")[end-2][end-1:end]],"-")
    name = "$perm_name and $node_range Peak Amplitude Time Series"

    # Plot permanent station moving average
    tseries = plot(times_avg, amps_avg, title = name, label="$perm_name Mov Avg", color="gray",
        linewidth=2, xtickfontsize=10, xlabel="Time", ylabel="Correlation Velocity Max Amp - filt: 0.1-0.5 Hz",
        legend = :outertopright, size=(900, 400), dpi=300, leftmargin= 7mm, bottommargin=5mm, linealpha=.6)

    # plot nodes
    for (ind, file) in enumerate(node_stations_)
        node_name, color = join(split(split(file,"/")[end],".")[3:4],"."), ""
        if occursin("NO.B1", node_name) # color B1 nodes green and intersecting nodes blue
            color = "green"
        else
            color = "blue"
        end
        test2 = DataFrame(CSV.File(file))
        test3 = mov_avg(test2, 20)
        times2, amplitudes2 = test2.date, test2.max_amp
        times3, amplitudes3 = test3.date, test3.moving_avg
        plot!(tseries, times3, amplitudes3, color=color, linealpha=.6, linewidth=2, label="$node_name Mov Avg")
    end
    save_path = joinpath(plot_root, "tseries_combined/$dir/CI.$(source_)_$(perm_name)_to_nodes.png")
    if !isdir(dirname(save_path))
        mkpath(dirname(save_path))
    end
    png(save_path)
end
function plot_simple_nodes(stations::Array{String,1}, source::String, vert_name::String, dir::String)
    name = "CI.$source Peak Amplitude Time Series at Intersect - $vert_name and B1"
    tseries = plot(title = name,linewidth=3, xtickfontsize=10, xlabel="Time",
        ylabel="Correlation Velocity Max Amp - filt: 0.1-0.5 Hz",
        legend = :outertopright, size=(900, 400), dpi=300, leftmargin= 7mm, bottommargin=5mm)
    for (ind, file) in enumerate(stations)
        rec_name = join(split(split(file,"/")[end],".")[3:4],".")
        series = DataFrame(CSV.File(file))
        series_m = mov_avg(series, 5)
        times, amplitudes = series.date, series.max_amp
        times_m, amplitudes_m = series_m.date, series_m.moving_avg
        color=""
        if occursin("B1", rec_name)
            color = "green"
        else
            color = "blue"
        end
        println(rec_name)
        plot!(tseries, 1:length(amplitudes_m), amplitudes_m, color=color, linealpha=.6,
                linewidth=2, label="$rec_name Mov Avg")
    end
    save_path = joinpath(plot_root,"tseries_combined/$dir/CI.$(source)_B1_$(vert_name)_intersect.png")
    if !isdir(dirname(save_path))
        mkpath(dirname(save_path))
    end
    png(save_path)
end
function wrap_plotter(vnodes::Array{String,1}, bnodes::Array{String,1}, reciever::String, vert_name, dir)
    nodes = vcat(vnodes, bnodes)
    for source_ in ["CHN","CJM","DEV","IPT","LPC","SNO","SVD","TA2"]
        try
            nodes_source = filter(x -> occursin(source_, split(x, ".")[2]), readdir("Tseries_nodes"))
            node_intersect = filter(x-> any(occursin.(join(split(x,".")[3:4],"."), nodes)), nodes_source)

            plot_simple_nodes(joinpath.("Tseries_nodes",node_intersect), source_, vert_name, dir)
            rec = "Tseries/CI.$(source_).$reciever.csv"
            println(rec)
            nodes = ["Tseries_nodes/CI.$(source_).$(node).jld2.csv" for node in nodes]
            println(nodes)
            plot_tseries(rec, joinpath.("Tseries_nodes", node_intersect), source_, dir)
        catch e
            println(e)
        end
    end
end


text:
# automatically generated by the FlatBuffers compiler, do not modify

UnionVector.eval(quote

FlatBuffers.@STRUCT struct BookReader
    books_read::Int32
end
FlatBuffers.@ALIGN(BookReader, 4)

BookReader(buf::AbstractVector{UInt8}) = FlatBuffers.read(BookReader, buf)
BookReader(io::IO) = FlatBuffers.deserialize(io, BookReader)

end)

max_stars_repo_name: rjkat/flatbuffers-julia

max_issues_repo_name: rjkat/flatbuffers-julia

max_forks_repo_name: rjkat/flatbuffers-julia


text:
export Recommender
export isbuilt, check_build_status, build!, recommend, predict, ranking

abstract type Recommender end

function check_build_status(recommender::Recommender)
    if !isbuilt(recommender)
        error("Recommender $(typeof(recommender)) is not built before making recommendation")
    end
end

isbuilt(recommender::Recommender) = true

function build!(recommender::Recommender; kwargs...)
    error("build! is not implemented for recommender type $(typeof(recommender))")
end

function recommend(recommender::Recommender, u::Int, k::Int, candidates::Array{Int})
    d = Dict{Int,Float64}()
    for candidate in candidates
        score = ranking(recommender, u, candidate)
        d[candidate] = score
    end
    sort(collect(d), by=tuple->last(tuple), rev=true)[1:k]
end

function predict(recommender::Recommender, u::Int, i::Int)
    error("predict is not implemented for recommender type $(typeof(recommender))")
end

# Return a ranking score of item i for user u
function ranking(recommender::Recommender, u::Int, i::Int)
    check_build_status(recommender)
    predict(recommender, u, i)
end


text:
module Types

using ValidatedNumerics;

immutable IntervalVector{T}
    intervals :: Array{Interval{T}, 1}
    sti :: Array{T, 1}
end

immutable Configuration{T}
  a :: Array{Interval{T}, 2}
  b :: IntervalVector{T}
  size :: Int
  sti_size :: Int
end

type Solver{T}
    current :: IntervalVector{T}
    previous :: IntervalVector{T}
    initial :: IntervalVector{T}
    system :: Configuration{T}
    roots :: Array{IntervalVector{T}, 1}
    iteration :: Int
end

immutable StopConditions
    precision :: AbstractFloat
    max_iterations :: Int
end

end


text:
@testset "Endgame" begin
    @testset "Cyclic 7" begin
        res = track.(total_degree(cyclic(7))...)

        @test count(is_success, res) == 924
        @test 4000 ≤ count(is_at_infinity, res) ≤ 4116
    end

    @testset "Hyperbolic - 6,6" begin
        # 2 solutions with multiplicity 6, projective
        @var x z
        y = 1
        # This has two roots of multiplicity 6 at the hyperplane z=0
        # each root has winding number 3
        F = [
            0.75 * x^4 + 1.5 * x^2 * y^2 - 2.5 * x^2 * z^2 + 0.75 * y^4 - 2.5 * y^2 * z^2 + 0.75 * z^4
            10 * x^2 * z + 10 * y^2 * z - 6 * z^3
        ]
        res = track.(total_degree(System(F))...)
        @test count(r -> r.winding_number == 3, res) == 12
    end

    @testset "Singular 1" begin
        # 1 singular solution with multiplicity 3
        @var x y
        z = 1
        F = [
            x^2 + 2 * y^2 + 2 * im * y * z,
            (18 + 3 * im) * x * y + 7 * im * y^2 - (3 - 18 * im) * x * z - 14 * y * z -
            7 * im * z^2,
        ]
        result = solve(System(F); start_system = :total_degree)
        @test nsingular(result) == 3
        @test nresults(result; only_nonsingular = true) == 1
    end

    @testset "Wilkinson $d" for d in [12]
        @var x
        f = System([expand(prod(x - i for i = 1:d))])
        res = track.(total_degree(f)...)
        @test all(is_success, res)
        @test round.(Int, real.(sort(first.(solution.(res)); by = abs))) == 1:d
        @test maximum(abs.(imag.(first.(solution.(res))))) < 1e-8
        @test count(r -> isnothing(r.winding_number), res) == d
    end

    @testset "(x-10)^$d" for d in [2, 8, 12]
        @var x
        f = System([(x - 10)^d])
        res = track.(total_degree(f)...)
        @test count(r -> r.winding_number == d, res) == d
    end

    @testset "Beyond Polyhedral Homotopy Example" begin
        @var x y
        f = [2.3 * x^2 + 1.2 * y^2 + 3x - 2y + 3, 2.3 * x^2 + 1.2 * y^2 + 5x + 2y - 5]
        res = track.(total_degree(System(f))...)
        @test count(is_success, res) == 2
        @test count(is_at_infinity, res) == 2
    end

    @testset "Winding Number Family d=$d" for d = 2:2:6
        @var x y
        a = [0.257, -0.139, -1.73, -0.199, 1.79, -1.32]
        f1 = (a[1] * x^d + a[2] * y) * (a[3] * x + a[4] * y) + 1
        f2 = (a[1] * x^d + a[2] * y) * (a[5] * x + a[6] * y) + 1
        res = track.(total_degree(System([f1, f2]))...)
        @test count(is_success, res) == d + 1
        @test count(is_at_infinity, res) == (d + 1)^2 - d - 1
    end

    @testset "Bacillus Subtilis" begin
        res = track.(total_degree(bacillus())...)
        @test count(is_success, res) == 44
    end

    @testset "Mohab" begin
        # Communicated by Mohab Safey El Din
        @var x y z
        F = [
            -9091098778555951517 * x^3 * y^4 * z^2 +
            5958442613080401626 * y^2 * z^7 +
            17596733865548170996 * x^2 * z^6 - 17979170986378486474 * x * y * z^6 -
            2382961149475678300 * x^4 * y^3 - 15412758154771986214 * x * y^3 * z^3 +
            133,
            -10798198881812549632 * x^6 * y^3 * z - 11318272225454111450 * x * y^9 -
            14291416869306766841 * y^9 * z - 5851790090514210599 * y^2 * z^8 +
            15067068695242799727 * x^2 * y^3 * z^4 +
            7716112995720175148 * x^3 * y * z^3 +
            171,
            13005416239846485183 * x^7 * y^3 + 4144861898662531651 * x^5 * z^4 -
            8026818640767362673 * x^6 - 6882178109031199747 * x^2 * y^4 +
            7240929562177127812 * x^2 * y^3 * z +
            5384944853425480296 * x * y * z^4 +
            88,
        ]
        @time res = track.(total_degree(System(F, [x, z, y]),)...)
        @test count(is_success, res) == 693

        # test for γ value where path jumping happened with to loose default config
        @time res =
            track.(total_degree(
                System(F, [x, z, y]),
                γ = -0.9132549847010242 + 0.4073884300256109im,
            )...)
        @test count(is_success, res) == 693
    end

    @testset "Ill-conditioned solution - look's almost diverging" begin
        sys, q₀ = fano_quintic()
        solver, starts = solver_startsolutions(
            sys;
            target_parameters = q₀,
            start_system = :polyhedral,
            seed = 0x56db29a3,
        )
        S = collect(starts)
        res = track(solver, S[1785])
        @test is_success(res)
    end
end


text:
export common_prefix_distance, generic_levenshtein, hamming_distance, levenshtein_distance, lcs_distance

"""
sim_common_prefix computes the length of the common prefix among
two strings represented as arrays
"""
function sim_common_prefix(a::T, b::T)::Int where {T <: Any}
    len_a::Int = length(a)
    len_b::Int = length(b)
    i::Int = 1
    min_len::Int = min(len_a, len_b)
    @inbounds while i <= min_len && a[i] == b[i]
    	i += 1
    end

    return i - 1
end

mutable struct CommonPrefixDistance
	calls::Int
	CommonPrefixDistance() = new(0)
end

function common_prefix_distance(a::T, b::T) where T
	return sim_common_prefix(a, b) / max(length(a), length(b))
end

"""
generic_levenshtein computes the edit distance between two strings,
this is a low level function. Please use dist_lev
"""
function generic_levenshtein(a::T, b::T, icost::Int, dcost::Int, rcost::Int)::Int where T
    if length(a) < length(b)
        a, b = b, a
    end

    alen::Int = length(a)
    blen::Int = length(b)

    alen == 0 && return blen
    blen == 0 && return alen

    C::Vector{Int} = Array(0:blen)

    prevA::Int = 0
    @inbounds for i in 1:alen
        prevA = i
        prevC::Int = C[1]
        j::Int = 1

        while j <= blen
            cost::Int = rcost
            if a[i] == b[j]
               cost = 0
            end
            C[j] = prevA
            j += 1
            prevA = min(C[j]+dcost, prevA+icost, prevC+cost)
            prevC = C[j]
        end

        C[j] = prevA
    end

    return prevA
end

"""
levenstein_distance
computes edit distance over two strings or arrays of comparable objects
"""
function levenshtein_distance(a::T, b::T)::Float64 where T
    generic_levenshtein(a, b, 1, 1, 1)
end

"""
lcs_distance
computes lcs distance over associated to the longest common subsequence; it accepts two strings or arrays of comparable objects
"""
function lcs_distance(a::T, b::T)::Float64 where T
    generic_levenshtein(a, b, 1, 1, 2)
end

# function kerrormatch(a::T1, b::T2, errors::Int)::Bool where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return alen == blen
#     blen == 0 && return true

#     C::Vector{Int} = Vector{Int}(0:blen)

#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
# 	    end

#         C[j] = prevA
#         if prevA <= errors
#             return true
#         end
#     end

#     return false
# end

# function best_match_levenshtein(a::T1, b::T2)::Int where {T1 <: Any,T2 <: Any}
#     # if length(a) < length(b)
#     #     a, b = b, a
#     # end

#     alen::Int = length(a)
#     blen::Int = length(b)

#     alen == 0 && return blen
#     blen == 0 && return alen

#     C::Vector{Int} = 1:blen |> collect

#     mindist = alen
#     @inbounds for i in 1:alen
#         prevA::Int = 0
#         prevC::Int = C[1]
#         j::Int = 1

#         while j <= blen
#             cost::Int = 1
#             if a[i] == b[j]
#                 cost = 0
#             end
#             C[j] = prevA
#             j += 1
#             prevA = min(C[j]+1, prevA+1, prevC+cost)
#             prevC = C[j]
#         end

#         C[j] = prevA
#         if prevA < mindist
#             mindist = prevA
#         end
#     end

#     return mindist
# end

"""
hamming_distance computes the hamming distance between two slices of integers
"""

function hamming_distance(a::T, b::T)::Float64 where T
    d::Int = 0

    @inbounds for i = 1:length(a)
        #if a[i] != b[i]
        #    d += 1
        #end
        d += Int(a[i] != b[i])
    end

    return d
end


text:
DIR = @__DIR__

if DIR != nothing
    const LIBRDKAFKA = DIR * "/../deps/librdkafka.so"
else
    const LIBRDKAFKA = Pkg.dir("RDKafka", "deps", "librdkafka.so")
end


text:
"""
# Description
Predict profile across all time steps for the true check.
    - if the problem is sequential, predict profiles from start to finish without the training, using only the initial profile as the initial condition.
    - if the problem is residual, predict profiles at each timestep using model-predicted difference between truth and physics-based model (KPP or TKE) prediction

Returns an n-length array of D-length vectors, where n is the number of training points and D is the
# Arguments
- `ℳ` (GP or NN). The model you will use to approximate the mapping.
- `𝒟` (ProfileData). The ProfileData object whose data 𝒢 will "test on", not necessarily the object that was used to train 𝒢.
                      If the problem is sequential, 𝒟 is the ProfileData object whose starting profile will be evolved forward using 𝒢.
# Keyword Arguments
- `postprocessed` (bool or String). If false, return whatever the model predicts directly (direct model output)
                                    If true, return the full predicted temperature profile calculated from the model output.
                                    If "both", return both.
"""
function predict(ℳ, 𝒟::ProfileData; postprocessed=true)
    D = length(𝒟.y[1])

    if typeof(𝒟.problem) <: Union{Sequential_KPP, Sequential_TKE}
        # Predict temperature profile from start to finish without the training data.

        gpr_prediction=Array{Array{Float64,1},1}()
        postprocessed_prediction=Array{Array{Float64,1},1}()
        t=1 # time index
        for (problem, n_x) in 𝒟.all_problems # n_x: number of predictors (i.e. time steps) for that problem

            post_pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            gpr__pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            post_pred_chunk[1] = unscale(𝒟.x[t], problem.scaling) # should = unscale(𝒟.x[i], problem.scaling) = postprocess_prediction(𝒟.x[i], 𝒟.y[i], problem)
            gpr__pred_chunk[1] = 𝒟.y[t] # residual -- should be zeros for this initial time step.

            for i in 1:n_x-1
                kpp_pred = scale(problem.evolve_physics_model_fn(post_pred_chunk[i]), problem.scaling)
                gpr__pred_chunk[i+1] = model_output(scale(post_pred_chunk[i], problem.scaling), t, ℳ, 𝒟)
                post_pred_chunk[i+1] = postprocess_prediction(kpp_pred, gpr__pred_chunk[i+1], problem)
                𝒟.convective_adjust(post_pred_chunk[i+1])
            end

            gpr_prediction = vcat(gpr_prediction, gpr__pred_chunk)
            postprocessed_prediction = vcat(postprocessed_prediction, post_pred_chunk)
            t += n_x
        end

        # THIS WORKS TOO BUT SLOWER:
        # gpr_prediction=Array{Array{Float64,1},1}(UndefInitializer(), 𝒟.Nt)
        # postprocessed_prediction=Array{Array{Float64,1},1}(UndefInitializer(), 𝒟.Nt)
        # kpp_pred = Array{Float64,1}(UndefInitializer(), D) # cache
        # t=1 # time index
        # for (problem, n_x) in 𝒟.all_problems # n_x: number of predictors (i.e. time steps) for that problem
        #
        #     postprocessed_prediction[t] = unscale(𝒟.x[t], problem.scaling) # should = unscale(𝒟.x[i], problem.scaling) = postprocess_prediction(𝒟.x[i], 𝒟.y[i], problem)
        #     gpr_prediction[t] = 𝒟.y[t] # residual -- should be zeros for this initial time step. Good sanity check.
        #
        #     for i in t:t+n_x-2
        #         kpp_pred .= scale(problem.evolve_physics_model_fn(postprocessed_prediction[i]), problem.scaling)
        #         gpr_prediction[i+1] = model_output(scale(postprocessed_prediction[i], problem.scaling), i, ℳ, 𝒟)
        #         postprocessed_prediction[i+1] = postprocess_prediction(kpp_pred, gpr_prediction[i+1], problem)
        #         𝒟.convective_adjust(postprocessed_prediction[i+1])
        #     end
        #
        #     t += n_x
        # end

    elseif typeof(𝒟.problem) <: ResidualProblem
        # Predict temperature profile from start to finish without the training data.

        gpr_prediction=Array{Array{Float64,1},1}()
        postprocessed_prediction=Array{Array{Float64,1},1}()
        t=1 # time index
        for (problem, n_x) in 𝒟.all_problems # n_x: number of predictors (i.e. time steps) for that problem

            post_pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            gpr__pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            post_pred_chunk[1] = unscale(𝒟.x[t], problem.scaling) # should = unscale(𝒟.x[i], problem.scaling) = postprocess_prediction(𝒟.x[i], 𝒟.y[i], problem)
            gpr__pred_chunk[1] = 𝒟.y[t] # residual -- should be zeros for this initial time step. Good sanity check.

            for i in 1:n_x-1
                kpp_pred = scale(problem.evolve_physics_model_fn(post_pred_chunk[i]), problem.scaling)
                gpr__pred_chunk[i+1] = model_output(kpp_pred, t, ℳ, 𝒟)
                post_pred_chunk[i+1] = postprocess_prediction(kpp_pred, gpr__pred_chunk[i+1], problem)
                𝒟.convective_adjust(post_pred_chunk[i+1])
            end

            gpr_prediction = vcat(gpr_prediction, gpr__pred_chunk)
            postprocessed_prediction = vcat(postprocessed_prediction, post_pred_chunk)
            t += n_x
        end

    elseif typeof(𝒟.problem) <: SequentialProblem

        # Predict temperature profile from start to finish without the training data.

        gpr_prediction=Array{Array{Float64,1},1}()
        postprocessed_prediction=Array{Array{Float64,1},1}()
        t=1 # time index
        for (problem, n_x) in 𝒟.all_problems # n_x: number of predictors (i.e. time steps) for that problem

            post_pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            gpr__pred_chunk = Array{Array{Float64,1},1}(UndefInitializer(), n_x)
            post_pred_chunk[1] = unscale(𝒟.x[t], problem.scaling) # should = unscale(𝒟.x[i], problem.scaling) = postprocess_prediction(𝒟.x[i], 𝒟.y[i], problem)
            gpr__pred_chunk[1] = 𝒟.y[t] # residual -- should be zeros for this initial time step. Good sanity check.

            for i in 1:n_x-1
                predictor = scale(post_pred_chunk[i], problem.scaling)
                # println(model_output(predictor, t, ℳ, 𝒟))
                gpr__pred_chunk[i+1] = model_output(predictor, t, ℳ, 𝒟)
                post_pred_chunk[i+1] = postprocess_prediction(predictor, gpr__pred_chunk[i+1], problem)
                𝒟.convective_adjust(post_pred_chunk[i+1])
            end

            gpr_prediction = vcat(gpr_prediction, gpr__pred_chunk)
            postprocessed_prediction = vcat(postprocessed_prediction, post_pred_chunk)
            t += n_x
        end

    elseif typeof(𝒟.problem) <: SlackProblem
        # Predict temperature profile at each timestep using model-predicted difference between truth and physics-based model (KPP or TKE) prediction
        gpr_prediction = [model_output(𝒟.x[i], i, ℳ, 𝒟) for i in 1:(𝒟.Nt)]
        postprocessed_prediction = get_postprocessed_predictions(𝒟.x, gpr_prediction, 𝒟.all_problems)
        ##
        # animate_physics_profile(ℳ, 𝒟)

    else; throw(error)
    end

    if postprocessed == "both";
        return (gpr_prediction, postprocessed_prediction)
    end

    if postprocessed
        return postprocessed_prediction
    end

    return gpr_prediction
end

function get_postprocessed_predictions(x, gpr_prediction, all_problems)

    result=Array{Array{Float64,1},1}()
    i=1
    for (problem, n_x) in all_problems # n_x: number of predictors for that problem
        result = vcat(result, GaussianProcess.postprocess_prediction(x[i : i+n_x-1], gpr_prediction[i : i+n_x-1], problem))
        i += n_x
    end

    return result
end

# function animate_physics_profile(ℳ, 𝒟)
#
#     variable = 𝒟.problem.variable # "T" or "wT"
#     xlims = x_lims[variable]
#
#     physics_data=Array{Array{Float64,1},1}()
#     for (problem, n_x) in 𝒟.all_problems # n_x: number of predictors for that problem
#         physics_data = vcat(physics_data, problem.physics_data)
#     end
#
#     # Compute error
#     n = length(physics_data)
#     total_error = 0.0
#     for i in 1:n
#         total_error += euclidean_distance(𝒟.vavg[i], physics_data[i])
#     end
#     println("ERROR $(total_error / n)")
#     println(total_error / n)
#
#     # Animate
#     animation_set = 1:30:(length(physics_data)-2)
#     anim = @animate for i in animation_set
#         day_string = string(floor(Int, 𝒟.t[i]/86400))
#         scatter(physics_data[i], 𝒟.zavg, label = "Physics") ####
#         plot!(𝒟.v[:,i], 𝒟.z, legend = :topleft, label = "LES", xlabel = "$(long_name[variable])", ylabel = "Depth [m]", title = "day " * day_string, xlims=xlims)
#     end
#
#     gif(anim, "physics_plot.gif")
# end

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
using Distributions

srand(1)

# make sure randcor works
@test sum(abs(inv(FactorAnalysis.randcor(40, 0.2))) .>= 1e-8) < 40*40/2


text:
@testset "Code replica" begin

    code = zeros(Int8, 2502)
    gpsl1 = GPSL1()
    Tracking.gen_code_replica!(
        code,
        gpsl1,
        1023e3,
        2.5e6,
        2.0,
        11,
        2480,
        SVector(-1, 0, 1),
        1
    )

    @test code[11:2492] == get_code.(gpsl1, (-1:2480) * 1023e3 / 2.5e6 .+ 2.0, 1)

    @testset "More than 1ms" begin
        code = zeros(Int8, 6502)
        gpsl1 = GPSL1()
        Tracking.gen_code_replica!(
            code,
            gpsl1,
            1023e3,
            2.5e6,
            2.0,
            11,
            6480,
            SVector(-1, 0, 1),
            1
        )

        @test code[11:6492] == get_code.(gpsl1, (-1:6480) * 1023e3 / 2.5e6 .+ 2.0, 1)
    end

    @testset "code_length is less than 1ms" begin
        code = zeros(Int8, 2502)
        gpsl1 = GPSL1()
        Tracking.gen_code_replica!(
            code,
            gpsl1,
            1023e3 * 3,
            7.5e6,
            2.0,
            11,
            2480,
            SVector(-1, 0, 1),
            1
        )

        @test code[11:2492] == get_code.(gpsl1, (-1:2480) * 1023e3 * 3 / 7.5e6 .+ 2.0, 1)
    end
end

@testset "Update code phase" begin
    gpsl1 = GPSL1()
    code_phase = 10
    code_frequency = 10Hz
    sampling_frequency = 100Hz
    num_samples = 2000
    bit_found = true
    phase = @inferred Tracking.update_code_phase(
        gpsl1,
        num_samples,
        code_frequency,
        sampling_frequency,
        code_phase,
        bit_found
    )
    @test phase ≈ mod(10 + 0.1 * 2000, 1023)
end


text:
using StatsBase: autocor
export estimate_delay, exponential_decay_fit, autocor

#####################################################################################
#                               Estimate Delay Times                                #
#####################################################################################
"""
    estimate_delay(s, method::String [, τs = 1:2:100]; kwargs...) -> τ

Estimate an optimal delay to be used in [`reconstruct`](@ref) or [`embed`](@ref).
The `method` can be one of the following:

* `"ac_zero"` : first delay at which the auto-correlation function becomes <0.
* `"ac_min"` : delay of first minimum of the auto-correlation function.
* `"mi_min"` : delay of first minimum of mutual information of `s` with itself
  (shifted for various `τs`).
  Keywords `nbins, binwidth` are propagated into [`mutualinformation`](@ref).
* `"exp_decay"` : [`exponential_decay_fit`](@ref) of the correlation function rounded
   to an integer (uses least squares on `c(t) = exp(-t/τ)` to find `τ`).
* `"exp_extrema"` : same as above but the exponential fit is done to the
  absolute value of the local extrema of the correlation function.

Both the mutual information and correlation function (`autocor`) are computed _only_
for delays `τs`. This means that the `min` methods can never return the first value
of `τs`!

The method `mi_min` is significantly more accurate than the others and also returns
good results for most timeseries. It is however the slowest method (but still quite fast!).
"""
function estimate_delay(x::AbstractVector, method::String,
    τs = 1:2:min(100, length(x)); kwargs...)

    issorted(τs) || error("`τs` must be sorted")

    if method=="ac_zero"
        c = autocor(x, τs; demean=true)
        i = 1
        while c[i] > 0 # Find 0 crossing
            i += 1
            if i == length(c)
                @warn "Did not cross 0 value, returning last `τ`."
                return τs[end]
            end
        end
        return τs[i]
    elseif method=="ac_min"
        c = autocor(x, τs, demean=true)
        return mincrossing(c, τs)
    elseif method=="mi_min"
        c = mutualinformation(x, τs; kwargs...)
        return mincrossing(c, τs)
    elseif method=="exp_decay"
        c = autocor(x, τs; demean=true)
        if any(x -> x ≤ 0, c)
            error("The correlation function has elements that are ≤ 0. "*
            "We can't fit an exponential to it. Please choose another method.")
        end
        τ = exponential_decay_fit(τs, c)
        return round(Int,τ)
    elseif method=="exp_extrema"
        c = autocor(x, τs; demean=true)
        max_ind, min_ind = localextrema(c)
        idxs = sort!(append!(max_ind, min_ind))
        ca = abs.(c[idxs])
        τa = τs[idxs]
        if any(x -> x ≤ 0, ca)
            error("The absolute correlation function has elements that are = 0. "*
            "We can't fit an exponential to it. Please choose another method.")
        end
        τ = exponential_decay_fit(τa, ca)
        return round(Int,τ)
    else
        throw(ArgumentError("Unknown method for `estimate_delay`."))
    end
end

function mincrossing(c, τs)
    i = 1
    while c[i+1] < c[i]
        i+= 1
        if i == length(c)-1
            @warn "Did not encounter a minimum, returning last `τ`."
            return τs[end]
        end
    end
    return τs[i]
end

"""
    exponential_decay_fit(x, y, weight = :equal) -> τ
Perform a least square fit of the form `y = exp(-x/τ)` and return `τ`.
Taken from:  http://mathworld.wolfram.com/LeastSquaresFittingExponential.html.
Assumes equal lengths of `x, y` and that `y ≥ 0`.

To use the method that gives more weight to small values of `y`, use `weight = :small`.
"""
function exponential_decay_fit(X, Y, weight = :equal)
    @inbounds begin
        L = length(X)
        b = if weight == :equal
            sy = sum(Y); sx = sum(X)
            a1 = sum(X[i]*Y[i]*log(Y[i]) for i in 1:L)
            a2 = sum(X[i]*Y[i] for i in 1:L)
            a3 = sum(Y[i]*log(Y[i]) for i in 1:L)
            a4 = sum(X[i]*X[i]*Y[i] for i in 1:L)
            b = (sy*a1 - a2*a3) / (sy*a4 - a2*a2)
        elseif weight == :small
            sx = sum(X)
            c1 = sum(X[i]*log(Y[i]) for i in 1:L)
            c2 = sum(log(y) for y in Y)
            c3 = sum(x*x for x in X)
            b = (L*c1 - sx*c2)/(L*c3 - sx*sx)
        end
        return -1.0/b
    end
end

"""
    localextrema(y) -> max_ind, min_ind
Find the local extrema of given array `y`, by scanning point-by-point. Return the
indices of the maxima (`max_ind`) and the indices of the minima (`min_ind`).
"""
function localextrema(y)
    @inbounds begin
        l = length(y)
        i = 1
        maxargs = Int[]
        minargs = Int[]
        if y[1] > y[2]
            push!(maxargs, 1)
        elseif y[1] < y[2]
            push!(minargs, 1)
        end

        for i in 2:l-1
            left = i-1
            right = i+1
            if  y[left] < y[i] > y[right]
                push!(maxargs, i)
            elseif y[left] > y[i] < y[right]
                push!(minargs, i)
            end
        end

        if y[l] > y[l-1]
            push!(maxargs, l)
        elseif y[l] < y[l-1]
            push!(minargs, l)
        end
        return maxargs, minargs
    end
end

#####################################################################################
#                               Estimate Delay Times                                #
#####################################################################################
export mutualinformation
"""
    mutualinformation(s, τs[; nbins, binwidth])

Calculate the mutual information between the time series `s` and its images
delayed by `τ` points for `τ` ∈ `τs`, using an _improvement_ of the method
outlined by Fraser & Swinney in [1].

## Description

The joint space of `s` and its `τ`-delayed image (`sτ`) is partitioned as a
rectangular grid, and the mutual information is computed from the joint and
marginal frequencies of `s` and `sτ` in the grid as defined in [1].
The mutual information values are returned in a vector of the same length
as `τs`.

If any of the optional keyword parameters is given, the grid will be a
homogeneous partition of the space where `s` and `sτ` are defined.
The margins of that partition will be divided in a number of bins equal
to `nbins`, such that the width of each bin will be `binwidth`, and the range
of nonzero values of `s` will be in the centre. If only of those two parameters
is given, the other will be automatically calculated to adjust the size of the
grid to the area where `s` and `sτ` are nonzero.

If no parameter is given, the space will be partitioned by a recursive bisection
algorithm based on the method given in [1].

Notice that the recursive method of [1] evaluates the joint frequencies of `s`
and `sτ` in each cell resulting from a partition, and stops when the data points
are uniformly distributed across the sub-partitions of the following levels.
For performance and stability reasons, the automatic partition method implemented
in this function is only used to divide the axes of the grid, using the marginal
frequencies of `s`.

## References
[1]: Fraser A.M. & Swinney H.L. "Independent coordinates for strange attractors
from mutual information" *Phys. Rev. A 33*(2), 1986, 1134:1140.
"""
function mutualinformation(s::AbstractVector{T}, τs::AbstractVector{Int};
    kwargs...) where {T}
    n = length(s)
    nτ = n-maximum(τs)
    perm = sortperm(s[1:nτ])
    # Choose partition method
    (bins, edges) = isempty(kwargs) ? _bisect(s[perm]) : _equalbins(s[perm]; kwargs...)
    f = zeros(length(bins))
    # Calculate the MI for each `τ`
    mi_values = zeros(length(τs))
    for (i, τ) ∈ enumerate(τs)
        sτ = view(s, τ+1:n)[perm] # delayed and reordered time series
        mi_values[i] = _mutualinfo!(f, sτ, bins, edges)
    end
    return mi_values
end

"""
    _mutualinfo!(f, sτ, edges, bins0)

Calculate the mutual information between the distribution of the delayed time
series `sτ` and its original image.

The two series are partitioned in a joint histogram with axes divided by the
points given in `edges`; the distribution of the original image is given by
`bins0`.
The values of `sτ` must be arranged such that all the points of the bin `(1,j)`
are contained in the first `bins0[1]` positions, the points of the bin `(2,j)
are contained in the following `bins[2]` positions, etc.

The vector `f` is used as a placeholder to pre-allocate the histogram.
"""
function _mutualinfo!(f::AbstractVector, sτ::AbstractVector,
    bins0::AbstractVector{<:Integer}, edges::AbstractVector)

    # Initialize values
    mi = 0.0                    # `h` will contain the result
    processed = 0               # number of points of `sτ` that have been used
    # Go through `bins0`
    n = length(sτ)
    logn = log2(n)
    for b in bins0
        frag = view(sτ, processed.+(1:b)) # fragment of `sτ` in the `b`-th bin
        processed += b
        _frequencies!(f, edges, sort(frag)) # store the values of the histogram in `f`
        for (i,g) in enumerate(f)
            (g != 0) && (mi += g/n*(log2(g/b/bins0[i])+logn))
        end
    end
    return mi
end

"""
    _frequencies!(f, s, edges)

Calculate a histogram of values of `s` along the bins defined by `edges`.
Both `s` and `edges` must be sorted ascendingly. The frequencies (counts)
of `s` in each bin will be stored in the pre-allocated vector `f`.
"""
function _frequencies!(f::AbstractVector{T}, edges::AbstractVector,
    s::AbstractVector) where {T}
    # Initialize the array of frequencies to zero
    fill!(f, zero(T))
    n = length(s)
    nbins = length(edges) - 1
    b = 1 # start in the first bin
    for i = 1:n
        if b ≥ nbins # the last bin is filled with the rest of points
            f[end] += n-i+1
            break
        end
        # when `s[i]` goes after the upper limit of the current bin,
        # move to the next bin and repeat until `s[i]` is found
        while s[i] > edges[b+1] && b < nbins
            b += 1
        end
        f[b] += 1 # add point to the current bin
    end
end

### Functions for creating histograms ###

# For type stability, all must return the same type of tuple
Histogram{T} = Tuple{Vector{<:Integer}, Vector{T}} where {T}

"""
    _equalbins(s[; nbins, binwidth])

Create a histogram of the sorted series `s` with bins of the same width.
Either the number of bins (`nbins`) or their width (`binwidth`) must be
given as keyword argument (**but not both**).
"""
function _equalbins(s::AbstractVector{T}; kwargs...)::Histogram{T} where {T}
    # only one of `nbins` or `binwidth` can be passed
    if length(kwargs) > 1
        throw(ArgumentError("the keyword argument can only be either `nbins` "*
        "or `binwidth`"))
    elseif haskey(kwargs, :nbins)     # fixed number of bins
        nbins = Int(kwargs[:nbins])
        r = range(s[1], stop=s[end], length=nbins+1)
    elseif haskey(kwargs, :binwidth)
        binwidth = T(kwargs[:binwidth])
        nbins = Int(div(s[end]-s[1],binwidth)+1)
        start = (s[1]+s[end]-nbins*binwidth)/2
        r = range(start, step=binwidth, length=nbins+1)
    else
        throw(ArgumentError("`nbins` or `binwidth` keyword argument required"))
    end
    bins = zeros(Int, nbins)
    _frequencies!(bins, r, s)
    edges = T[minimum(s); s[cumsum(bins)]]
    return (bins, edges)
end

"""
    _bisect(s)

Create a partition histogram of the sorted series `s` with a partition of its
space defined by a recursive bisection method. The first level partition
divides `s` in two segments with equal number of points; each
partition is divided into two further sub-pantitions, etc.,
until the distribution of the points in the highest level
subpartition is homogeneous.
"""
function _bisect(s::AbstractVector{T})::Histogram{T} where {T}
    # set up the vectors to return
    bins = Int[]
    edges = T[minimum(s)]
    _bisect!(bins, edges, s) # recursive step
    return (bins, edges)
end

function _bisect!(bins, edges, s)
    if _uniformtest(s) # no further partitions if it is uniform
        push!(bins, length(s))
        push!(edges, s[end])
    else
        n = length(s)
        half = div(n, 2)
        # first half
        _bisect!(bins, edges, view(s, 1:half))
        # second half
        _bisect!(bins, edges, view(s, half+1:n))
    end
    return
end

"""
    _uniformtest(s)

Test uniformity in the values of the sorted vector `s`.
"""
function _uniformtest(s)
    n = length(s)
    # less than 10 points is not enough to test
    # between 10 and 19 points, divide the range of values in two sub-ranges
    # for 20 or more, divide the range in four sub-ranges
    # and make a chi-squared test
    if n <  10
        return true
    else
        if n < 20
            nbins = 2
            critical_chisq = 2.706 # df=1, p=0.1
        else
            nbins = 4
            critical_chisq = 6.251 # df=3, p=0.1
        end
    end
    binwidth = (s[end]-s[1])/nbins
    expected = round(Int, n/nbins)
    chisq = 0.0
    # set up first iteration
    first_item = 1
    for i = 1:nbins
        next_first = searchsortedfirst(s, s[1]+i*binwidth; lt=≤)
        observed = next_first - first_item
        chisq += (observed-expected)^2/expected
        (next_first > n) && break # stop if arrived at the end of `s`
        first_item = next_first
    end
    return (chisq < critical_chisq)
end