text:
const Poly = Union{ArbPoly,AcbPoly}
const Series = Union{ArbSeries,AcbSeries}

##
## Length and degree
##
Base.length(p::Union{arb_poly_struct,ArbPoly}) =
    ccall(@libarb(arb_poly_length), Int, (Ref{arb_poly_struct},), p)
Base.length(p::Union{acb_poly_struct,AcbPoly}) =
    ccall(@libarb(acb_poly_length), Int, (Ref{acb_poly_struct},), p)
Base.length(p::Series) = p.degree + 1

# Only define size for Poly. Series inherits the size from being
# Numbers
Base.size(p::Poly) = (length(p),)
Base.size(p::Poly, d::Integer) = d < 1 ? throw(BoundsError()) : d == 1 ? length(p) : 1

# We have to define these separately so that it overrides the default
# correctly
degree(p::ArbSeries) = p.degree
degree(p::AcbSeries) = p.degree

_degree(p::Series, q::Series) = min(degree(p), degree(q))

##
## Get and set coefficients
##

Base.eltype(::Type{<:Union{ArbPoly,ArbSeries}}) = Arb
Base.eltype(::Type{<:Union{AcbPoly,AcbSeries}}) = Acb

coeffs(p::Union{Poly,Series}) = [@inbounds p[i] for i = 0:degree(p)]

Base.firstindex(::Union{Poly,Series}) = 0
Base.lastindex(p::Union{Poly,Series}) = degree(p)

Base.checkbounds(::Type{Bool}, ::Poly, i::Integer) = i >= 0
Base.checkbounds(::Type{Bool}, p::Series, i::Integer) = 0 <= i <= degree(p)
Base.checkbounds(p::Union{Poly,Series}, i::Integer) =
    checkbounds(Bool, p, i) || throw(BoundsError(p, i))

Base.@propagate_inbounds function Base.getindex(
    p::T,
    i::Integer,
) where {T<:Union{Poly,Series}}
    @boundscheck checkbounds(p, i)
    return get_coeff!(eltype(T)(prec = precision(p)), p, i)
end

Base.getindex(p::Union{Poly,Series}, ::Colon) = coeffs(p)
Base.getindex(p::Union{Poly,Series}, I::AbstractRange{<:Integer}) = [p[i] for i in I]

Base.@propagate_inbounds function Base.setindex!(
    p::Union{ArbPoly,ArbSeries},
    x::ArbLike,
    i::Integer,
)
    @boundscheck checkbounds(p, i)
    set_coeff!(p, i, x)
    return x
end

Base.@propagate_inbounds function Base.setindex!(
    p::Union{AcbPoly,AcbSeries},
    x::AcbLike,
    i::Integer,
)
    @boundscheck checkbounds(p, i)
    set_coeff!(p, i, x)
    return x
end

Base.@propagate_inbounds function Base.setindex!(p::Union{Poly,Series}, x, i::Integer)
    @boundscheck checkbounds(p, i)
    set_coeff!(p, i, eltype(p)(x, prec = precision(p)))
    return x
end

Base.@propagate_inbounds function ref(p::Union{ArbPoly,ArbSeries}, i::Integer)
    @boundscheck 0 <= i <= degree(p) || throw(BoundsError(p, i))
    ptr = cstruct(p).coeffs + i * sizeof(arb_struct)
    return ArbRef(ptr, precision(p), cstruct(p))
end

Base.@propagate_inbounds function ref(p::Union{AcbPoly,AcbSeries}, i::Integer)
    @boundscheck 0 <= i <= degree(p) || throw(BoundsError(p, i))
    ptr = cstruct(p).coeffs + i * sizeof(acb_struct)
    return AcbRef(ptr, precision(p), cstruct(p))
end

##
## Constructors
##

for (TPoly, TSeries) in [(:ArbPoly, :ArbSeries), (:AcbPoly, :AcbSeries)]
    @eval $TPoly(p::cstructtype($TPoly); prec::Integer = DEFAULT_PRECISION[]) =
        set!($TPoly(prec = prec), p)

    @eval function $TPoly(coeff; prec::Integer = _precision(coeff))
        p = $TPoly(prec = prec)
        p[0] = coeff
        return p
    end

    @eval function $TPoly(coeffs::AbstractVector; prec::Integer = _precision(first(coeffs)))
        p = fit_length!($TPoly(prec = prec), length(coeffs))
        @inbounds for i = 1:length(coeffs)
            p[i-1] = coeffs[i]
        end
        return p
    end

    @eval $TPoly(p::Union{$TPoly,$TSeries}; prec::Integer = precision(p)) =
        set!($TPoly(prec = prec), p)
end
function AcbPoly(p::Union{ArbPoly,ArbSeries}; prec = precision(p))
    res = fit_length!(AcbPoly(prec = prec), length(p))
    @inbounds for i = 0:degree(p)
        res[i] = p[i]
    end
    return res
end

for (TSeries, TPoly) in [(:ArbSeries, :ArbPoly), (:AcbSeries, :AcbPoly)]
    @eval $TSeries(
        p::cstructtype($TSeries);
        degree::Integer = degree(p),
        prec::Integer = DEFAULT_PRECISION[],
    ) = set!($TSeries(degree = degree, prec = prec), p)

    @eval function $TSeries(coeff; degree::Integer = 0, prec::Integer = _precision(coeff))
        p = $TSeries(degree = degree, prec = prec)
        p[0] = coeff
        return p
    end

    @eval function $TSeries(
        coeffs::AbstractVector;
        degree::Integer = length(coeffs) - 1,
        prec::Integer = _precision(first(coeffs)),
    )
        p = fit_length!($TSeries(degree = degree, prec = prec), degree + 1)
        @inbounds for i = 1:min(length(coeffs), degree + 1)
            p[i-1] = coeffs[i]
        end
        return p
    end

    @eval $TSeries(p::Union{$TPoly,$TSeries}; prec::Integer = precision(p)) =
        set!($TSeries(degree = degree(p), prec = prec), p)
end
function AcbSeries(p::Union{ArbPoly,ArbSeries}; degree = degree(p), prec = precision(p))
    res = AcbSeries(degree = degree, prec = prec)
    @inbounds for i = 0:min(Arblib.degree(p), degree)
        res[i] = p[i]
    end
    return res
end

Base.zero(p::T) where {T<:Poly} = T(prec = precision(p))
Base.one(p::Poly) = one!(zero(p))

Base.zero(p::T) where {T<:Series} = T(degree = degree(p), prec = precision(p))
Base.one(p::Series) = one!(zero(p))

Base.zero(::Type{T}) where {T<:Union{Poly,Series}} = T()
Base.one(::Type{T}) where {T<:Union{Poly,Series}} = one!(zero(T))

fromroots(::Type{ArbPoly}, roots::ArbVector; prec::Integer = DEFAULT_PRECISION[]) =
    product_roots!(ArbPoly(prec = prec), roots)
fromroots(::Type{ArbPoly}, roots::AbstractVector; prec::Integer = DEFAULT_PRECISION[]) =
    fromroots(ArbPoly, ArbVector(roots, prec = prec), prec = prec)

fromroots(
    ::Type{ArbPoly},
    real_roots::ArbVector,
    complex_roots::AcbVector;
    prec::Integer = DEFAULT_PRECISION[],
) = product_roots_complex!(
    ArbPoly(prec = prec),
    real_roots,
    length(real_roots),
    complex_roots,
    length(complex_roots),
)
fromroots(
    ::Type{ArbPoly},
    real_roots::AbstractVector,
    complex_roots::AbstractVector;
    prec::Integer = DEFAULT_PRECISION[],
) = fromroots(
    ArbPoly,
    ArbVector(real_roots, prec = prec),
    AcbVector(complex_roots, prec = prec),
    prec = prec,
)

fromroots(::Type{AcbPoly}, roots::AcbVector; prec::Integer = DEFAULT_PRECISION[]) =
    product_roots!(AcbPoly(prec = prec), roots)
fromroots(::Type{AcbPoly}, roots::AbstractVector; prec::Integer = DEFAULT_PRECISION[]) =
    fromroots(AcbPoly, AcbVector(roots, prec = prec), prec = prec)

Base.copy(p::Union{Poly,Series}) = set!(zero(p), p)

##
## Arithmetic
##

Base.:-(p::Union{Poly,Series}) = neg!(zero(p), p)

Base.:+(p::T, q::T) where {T<:Poly} = add!(T(prec = _precision((p, q))), p, q)
Base.:-(p::T, q::T) where {T<:Poly} = sub!(T(prec = _precision((p, q))), p, q)
Base.:*(p::T, q::T) where {T<:Poly} = mul!(T(prec = _precision((p, q))), p, q)

function Base.:+(p::AcbPoly, q::ArbPoly)
    res = AcbPoly(q, prec = _precision((p, q)))
    return add!(res, p, res)
end
Base.:+(p::ArbPoly, q::AcbPoly) = q + p

function Base.:-(p::AcbPoly, q::ArbPoly)
    res = AcbPoly(q, prec = _precision((p, q)))
    return sub!(res, p, res)
end
function Base.:-(p::ArbPoly, q::AcbPoly)
    res = AcbPoly(p, prec = _precision((p, q)))
    return sub!(res, res, q)
end

function Base.:*(p::AcbPoly, q::ArbPoly)
    res = AcbPoly(q, prec = _precision((p, q)))
    return mul!(res, p, res)
end
Base.:*(p::ArbPoly, q::AcbPoly) = q * p

# We can't define these as `(p::T, q::T) where {T <: Series}` due to
# method ambiguity issues.
for T in [ArbSeries, AcbSeries]
    @eval function Base.:+(p::$T, q::$T)
        deg = _degree(p, q)
        return add_series!($T(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
    end
    @eval function Base.:-(p::$T, q::$T)
        deg = _degree(p, q)
        return sub_series!($T(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
    end
    @eval function Base.:*(p::$T, q::$T)
        deg = _degree(p, q)
        return mullow!($T(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
    end
    @eval function Base.:/(p::$T, q::$T)
        deg = _degree(p, q)
        return div_series!($T(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
    end
end

function Base.:+(p::AcbSeries, q::ArbSeries)
    deg = _degree(p, q)
    res = AcbSeries(q, degree = deg, prec = _precision((p, q)))
    return add_series!(res, p, res, deg + 1)
end
Base.:+(p::ArbSeries, q::AcbSeries) = q + p

function Base.:-(p::AcbSeries, q::ArbSeries)
    deg = _degree(p, q)
    res = AcbSeries(q, degree = deg, prec = _precision((p, q)))
    return sub_series!(res, p, res, deg + 1)
end
function Base.:-(p::ArbSeries, q::AcbSeries)
    deg = _degree(p, q)
    res = AcbSeries(p, degree = deg, prec = _precision((p, q)))
    return sub_series!(res, res, q, deg + 1)
end

function Base.:*(p::AcbSeries, q::ArbSeries)
    deg = _degree(p, q)
    res = AcbSeries(q, degree = deg, prec = _precision((p, q)))
    return mullow!(res, p, res, deg + 1)
end
Base.:*(p::ArbSeries, q::AcbSeries) = q * p

function Base.:/(p::AcbSeries, q::ArbSeries)
    deg = _degree(p, q)
    res = AcbSeries(q, degree = deg, prec = _precision((p, q)))
    return div_series!(res, p, res, deg + 1)
end
function Base.:/(p::ArbSeries, q::AcbSeries)
    deg = _degree(p, q)
    res = AcbSeries(p, degree = deg, prec = _precision((p, q)))
    return div_series!(res, res, q, deg + 1)
end

Base.inv(p::Series) = inv_series!(zero(p), p, degree(p) + 1)

# TODO: Implement separate ÷ and rem as well? They are only
# implemented for vectors in Arb so we would have to call those
# functions manually.
function Base.divrem(p::T, q::T) where {T<:Poly}
    quotient = T(prec = _precision((p, q)))
    remainder = T(prec = _precision((p, q)))
    divrem!(quotient, remainder, p, q)
    return quotient, remainder
end

##
## Scalar arithmetic
##

for (T, Tel, Tel_inplace) in [
    (Union{ArbPoly,ArbSeries}, Real, Union{ArbOrRef,ArfLike,Unsigned,Integer}),
    (Union{AcbPoly,AcbSeries}, Number, Union{AcbOrRef,ArbOrRef,Unsigned,Integer}),
]
    # Since we use inplace methods we need to manually normalise the
    # polynomials after the operation. This is to make sure that it
    # correctly recognizes the possibly new degree of the polynomial.
    # For example for iszero(ArbPoly(1) - 1) to work.

    @eval function Base.:+(p::$T, c::$Tel_inplace)
        res = copy(p)
        # Handle p being the zero polynomial
        iszero(p) && return set_coeff!(res, 0, c)
        res0 = ref(res, 0)
        add!(res0, res0, c)
        return normalise!(res)
    end
    @eval Base.:+(p::$T, c::$Tel) = p + convert(eltype(p), c)

    @eval function Base.:-(p::$T, c::$Tel_inplace)
        res = copy(p)
        # Handle p being the zero polynomial
        if iszero(p)
            set_coeff!(res, 0, c)
            res0 = ref(res, 0)
            neg!(res0, res0)
            return res
        end
        res0 = ref(res, 0)
        sub!(res0, res0, c)
        return normalise!(res)
    end
    @eval Base.:-(p::$T, c::$Tel) = p - convert(eltype(p), c)
    @eval function Base.:-(c::$Tel_inplace, p::$T)
        res = -p
        # Handle p being the zero polynomial
        iszero(p) && return set_coeff!(res, 0, c)
        res0 = ref(res, 0)
        add!(res0, res0, c)
        return normalise!(res)
    end
    @eval Base.:-(c::$Tel, p::$T) = convert(eltype(p), c) - p

    @eval Base.:*(p::$T, c::$Tel) = mul!(zero(p), p, convert(eltype(p), c))

    @eval Base.:/(p::$T, c::$Tel) = div!(zero(p), p, convert(eltype(p), c))
end

# Avoid conversion in these cases
@eval Base.:*(p::Union{ArbPoly,ArbSeries}, c::ArbOrRef) = mul!(zero(p), p, c)
@eval Base.:*(p::Union{AcbPoly,AcbSeries}, c::AcbOrRef) = mul!(zero(p), p, c)

@eval Base.:/(p::Union{ArbPoly,ArbSeries}, c::ArbOrRef) = div!(zero(p), p, c)
@eval Base.:/(p::Union{AcbPoly,AcbSeries}, c::AcbOrRef) = div!(zero(p), p, c)

# Promotion to complex
for (T, complexT) in [(ArbPoly, AcbPoly), (ArbSeries, AcbSeries)]
    @eval function Base.:+(p::$T, c::AcbOrRef)
        res = $complexT(p)
        res0 = ref(res, 0)
        add!(res0, res0, c)
        return res
    end

    @eval function Base.:-(p::$T, c::AcbOrRef)
        res = $complexT(p)
        res0 = ref(res, 0)
        sub!(res0, res0, c)
        return res
    end
    @eval function Base.:-(c::AcbOrRef, p::$T)
        res = $complexT(p)
        neg!(res, res)
        res0 = ref(res, 0)
        add!(res0, res0, c)
        return res
    end

    @eval function Base.:*(p::$T, c::AcbOrRef)
        res = $complexT(p)
        return mul!(res, res, c)
    end

    @eval function Base.:/(p::$T, c::AcbOrRef)
        res = $complexT(p)
        return div!(res, res, c)
    end
end
Base.:+(p::Union{ArbPoly,ArbSeries}, c::Complex) = p + convert(Acb, c)
Base.:-(p::Union{ArbPoly,ArbSeries}, c::Complex) = p - convert(Acb, c)
Base.:-(c::Complex, p::Union{ArbPoly,ArbSeries}) = convert(Acb, c) - p
Base.:*(p::Union{ArbPoly,ArbSeries}, c::Complex) = p * convert(Acb, c)
Base.:/(p::Union{ArbPoly,ArbSeries}, c::Complex) = p / convert(Acb, c)

Base.:+(c::Number, p::Union{Poly,Series}) = p + c
Base.:*(c::Number, p::Union{Poly,Series}) = p * c

function Base.:/(c::Real, p::ArbSeries)
    res = inv(p)
    return mul!(res, res, convert(Arb, c))
end
function Base.:/(c::Union{AcbOrRef,Complex}, p::ArbSeries)
    res = AcbSeries(inv(p))
    return mul!(res, res, convert(Acb, c))
end
function Base.:/(c::Number, p::AcbSeries)
    res = inv(p)
    return mul!(res, res, convert(Acb, c))
end

##
## Composition
##

taylor_shift(p::Union{Poly,Series}, c) = taylor_shift!(zero(p), p, convert(eltype(p), c))

compose(p::T, q::T) where {T<:Poly} = compose!(T(prec = _precision((p, q))), p, q)
function compose(p::T, q::T) where {T<:Series}
    iszero(ref(q, 0)) ||
        throw(ArgumentError("constant term of q must be zero, got q[0] = $(q[0])"))
    deg = _degree(p, q)
    res = T(degree = deg, prec = _precision((p, q)))
    return compose_series!(res, p, q, deg + 1)
end

function revert(p::Series)
    degree(p) >= 1 || throw(ArgumentError("p must have degree at least 1"))
    iszero(ref(p, 0)) ||
        throw(ArgumentError("constant term of p must be zero, got p[0] = $(p[0])"))
    !iszero(ref(p, 1)) ||
        throw(ArgumentError("linear term of p must be non-zero, got p[0] = $(p[0])"))
    return revert_series!(zero(p), p, degree(p) + 1)
end

##
## Evaluation
##

(p::Union{ArbPoly,ArbSeries})(x::ArbOrRef) = evaluate!(Arb(prec = precision(p)), p, x)

(p::Union{ArbPoly,ArbSeries})(x::Real) =
    evaluate!(Arb(prec = precision(p)), p, convert(Arb, x))

(p::Union{Poly,Series})(x::AcbOrRef) = evaluate!(Acb(prec = precision(p)), p, x)

(p::Union{Poly,Series})(x) = evaluate!(Acb(prec = precision(p)), p, convert(Acb, x))

function evaluate2(p::Union{ArbPoly,ArbSeries}, x::ArbOrRef)
    res1, res2 = Arb(prec = precision(p)), Arb(prec = precision(p))
    evaluate2!(res1, res2, p, x)
    return (res1, res2)
end

function evaluate2(p::Union{ArbPoly,ArbSeries}, x::Real)
    res1, res2 = Arb(prec = precision(p)), Arb(prec = precision(p))
    evaluate2!(res1, res2, p, convert(Arb, x))
    return (res1, res2)
end

function evaluate2(p::Union{Poly,Series}, x::AcbOrRef)
    res1, res2 = Acb(prec = precision(p)), Acb(prec = precision(p))
    evaluate2!(res1, res2, p, x)
    return (res1, res2)
end

evaluate2(p::Union{Poly,Series}, x::T) where {T} = evaluate2(p, convert(Acb, x))

##
## Differentiation and integration
##

derivative(p::Poly) = derivative!(zero(p), p)
derivative(p::T) where {T<:Series} =
    derivative!(T(degree = degree(p) - 1, prec = precision(p)), p)

integral(p::Poly) = integral!(zero(p), p)
integral(p::T) where {T<:Series} =
    integral!(T(degree = degree(p) + 1, prec = precision(p)), p)

##
## Power methods
##

Base.:^(p::Poly, e::Integer) = pow!(zero(p), p, convert(UInt, e))

function Base.:^(p::ArbSeries, q::ArbSeries)
    deg = _degree(p, q)
    return pow_series!(ArbSeries(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
end
function Base.:^(p::AcbSeries, q::AcbSeries)
    deg = _degree(p, q)
    return pow_series!(AcbSeries(degree = deg, prec = _precision((p, q))), p, q, deg + 1)
end

Base.:^(p::ArbSeries, e::Real) = pow_arb_series!(zero(p), p, convert(Arb, e), length(p))
function Base.:^(p::ArbSeries, e::Number)
    res = AcbSeries(p)
    return pow_acb_series!(res, res, convert(Acb, e), length(p))
end
function Base.:^(x::Real, p::ArbSeries)
    res = ArbSeries(x, degree = degree(p), prec = precision(p))
    return pow_series!(res, res, p, length(p))
end
function Base.:^(x::Number, p::ArbSeries)
    res = AcbSeries(x, degree = degree(p), prec = precision(p))
    return pow_series!(res, res, AcbSeries(p), length(p))
end

Base.:^(p::AcbSeries, e::Number) = pow_acb_series!(zero(p), p, convert(Acb, e), length(p))
function Base.:^(x::Number, p::AcbSeries)
    res = AcbSeries(x, degree = degree(p), prec = precision(p))
    return pow_series!(res, res, p, length(p))
end

# Disambiguation
Base.:^(p::ArbSeries, e::Integer) = pow_arb_series!(zero(p), p, convert(Arb, e), length(p))
Base.:^(p::AcbSeries, e::Integer) = pow_acb_series!(zero(p), p, convert(Acb, e), length(p))
Base.:^(p::ArbSeries, e::Rational) = pow_arb_series!(zero(p), p, convert(Arb, e), length(p))
Base.:^(p::AcbSeries, e::Rational) = pow_acb_series!(zero(p), p, convert(Acb, e), length(p))

##
## Series methods
##

for f in [:sqrt, :log, :log1p, :exp, :sin, :cos, :tan, :atan, :sinh, :cosh]
    @eval Base.$f(p::Series) = $(Symbol(f, :_series!))(zero(p), p, length(p))
end

Base.asin(p::ArbSeries) = asin_series!(zero(p), p, length(p))
Base.acos(p::ArbSeries) = acos_series!(zero(p), p, length(p))

rsqrt(p::Series) = rsqrt_series!(zero(p), p, length(p))

Base.sinpi(p::Series) = sin_pi_series!(zero(p), p, length(p))
Base.cospi(p::Series) = cos_pi_series!(zero(p), p, length(p))
cotpi(p::Series) = cot_pi_series!(zero(p), p, length(p))
# Julias definition of sinc is equivalent to Arbs definition of sincpi
Base.sinc(p::ArbSeries) = sinc_pi_series!(zero(p), p, length(p))

function Base.sincos(p::Series)
    s, c = zero(p), zero(p)
    sin_cos_series!(s, c, p, length(p))
    return (s, c)
end
function sincospi(p::Series)
    s, c = zero(p), zero(p)
    sin_cos_pi_series!(s, c, p, length(p))
    return (s, c)
end
function sinhcosh(p::Series)
    s, c = zero(p), zero(p)
    sinh_cosh_series!(s, c, p, length(p))
    return (s, c)
end


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# common definitions

const Float = Float64
# const Int = Int64
# const UInt = UInt64

# float and int types for storing all-pairs shortest paths data
# can reduce precision in order to reduce memory used
const FloatSpTime = Float64 # precision for storing shortest path travel times
const FloatSpDist = Float64 # precision for storing shortest path travel distances
const IntRNode = Int16 # precision for storing number of nodes in rGraph
const IntFadj = Int8 # precision for storing maximum number of nodes adjacent to any node in rGraph (= length of longest vector in network.rGraph.fadjList)

const sourceDir = @__DIR__
const sourcePath = sourceDir # for compat
const jemssDir = realpath(joinpath(sourceDir, ".."))

const Deployment = Vector{Int} # deployment[i] gives index of station to deploy ambulance i to
const CompTable = Array{Int,2} # compTable[i,j] gives number of ambulances to place at station j when there are i free ambulances
const MultiCompTable = Vector{Vector{Vector{Int}}} # multiCompTable[i][j][k] is jth option for number of ambs to assign to station k with i free ambs
const NestedCompTable = Vector{Int} # nested compliance table can be represented as a list of station indices; nestedCompTable[1:i] gives indices of stations that i free ambulances should be placed at
const PriorityList = Vector{Int} # a priority list gives the order of preference for which ambulances should be redeployed (moved up) to stations; priorityList[1] has the index of the station with highest priority

const maxNumThreads = nthreads() # value set by environment variable JULIA_NUM_THREADS

# run modes
const checkMode = true # for data checking, e.g. assertions that are checked frequently

# file chars
const delimiter = ','
const newline = "\r\n"

# misc globals
configFileDir = "" # hacky; to allow using "$configFileDir" in config files
genConfigFileDir = "" # hacky; to allow using "$genConfigFileDir" in gen config files

# misc null values
const nullIndex = -1
const nullX = -200 # outside lat/lon range
const nullY = -200
const nullTime = -1.0
const nullDist = -1.0
const nullHist = fit(Histogram, [], 0:0)
NullHist() = deepcopy(nullHist)
nullFunction() = nothing

# call/travel priorities
@enum Priority nullPriority=0 highPriority=1 medPriority=2 lowPriority=3
const priorities = setdiff([instances(Priority)...], [nullPriority])
const numPriorities = length(priorities)

# ambulance classes
@enum AmbClass nullAmbClass=0 als=1 bls=2 # als = advanced life support, bls = basic life support

# event types/forms
@enum EventForm nullEvent ambGoesToSleep ambWakesUp callArrives considerDispatch ambDispatched ambMobilised ambReachesCall ambGoesToHospital ambReachesHospital ambBecomesFree ambReturnsToStation ambReachesStation considerMoveUp ambMoveUpToStation

# ambulance statuses
@enum AmbStatus ambNullStatus ambSleeping ambIdleAtStation ambMobilising ambGoingToCall ambAtCall ambGoingToHospital ambAtHospital ambFreeAfterCall ambReturningToStation ambMovingUpToStation
@enum AmbStatusSet ambWorking ambBusy ambFree ambTravelling ambGoingToStation

# call statuses
@enum CallStatus callNullStatus callScreening callQueued callWaitingForAmb callOnSceneTreatment callGoingToHospital callAtHospital callProcessed

# route statuses
@enum RouteStatus routeNullStatus routeBeforeStartNode routeOnPath routeAfterEndNode

@enum MoveUpModule nullMoveUpModule compTableModule ddsmModule dmexclpModule multiCompTableModule priorityListModule zhangIpModule temp0Module temp1Module temp2Module

ext: jl


text:
using Jecco.KG_3_1

grid = SpecCartGrid3D(
    x_min            = -5.0,
    x_max            =  5.0,
    x_nodes          =  64,
    y_min            = -5.0,
    y_max            =  5.0,
    y_nodes          =  64,
    u_max            =  1.001,
    u_domains        =  1,
    u_nodes          =  48,
)

potential = SquarePotential()

evoleq = AffineNull(
    potential      = potential,
)

id = Sine2D(
    kx   = 1,
    ky   = 2,
    Lx   = grid.x_max - grid.x_min,
    Ly   = grid.y_max - grid.y_min,
)

io = InOut(
    out_bulk_every_t            = 0.04,
    checkpoint_every_walltime_hours = 1,
    remove_existing   = true,
)

integration = Integration(
    dt              = 0.002,
    ODE_method      = KG_3_1.RK4(),
    adaptive        = false,
    tmax            = 4.0,
)

run_model(grid, id, evoleq, integration, io)


text:
module KoalaPlots

export recordfig

import Plots.savefig

"""
    recordfig(filename, models...)

Save the current `Plots` figure as a PNG file called `filename` (which
must include ".png" extension) and generate a markdown report
containing an embedded version of the figure and the output of
`showall(model)` for each `model` in `models` (ususually metadata
associated with the plot).

For example, if the filename is "assets/myplot.png" then two files
"assets/myplot.png" and "assets/myplot.md" are created.

"""
function recordfig(figure_filename::String, models...)
    local_figure_filename = match(r"([^\/]*)$", figure_filename)[1]
    title = match(r"(.*)\..*$", local_figure_filename)[1]
    md_filename = match(r"(.*)\..*$", figure_filename)[1] * ".md"
    savefig(figure_filename)
    open(md_filename,"w") do md_file
        write(md_file, "# $title\n\n")
        write(md_file, "![$local_figure_filename]($local_figure_filename)\n\n")
        for model in models
            showall(md_file, model)
            println(md_file, "\n\n")
        end
    end
    info("Recorded figure and metadata at $md_filename.")
end

end # of module


text:
@cfgdef struct HGFRobertaConfig <: HGFConfig
  vocab_size::Int = 30522
  hidden_size::Int = 768
  num_hidden_layers::Int = 12
  num_attention_heads::Int = 12
  intermediate_size::Int = 3072
  hidden_act::String = "gelu"
  hidden_dropout_prob::Float64 = 0.1
  attention_probs_dropout_prob::Float64 = 0.1
  max_position_embeddings::Int = 512
  type_vocab_size::Int = 2
  initializer_range::Float32 = 0.02
  layer_norm_eps::Float32 = 1e-12
  pad_token_id::Int = 1
  bos_token_id::Int = 0
  eos_token_id::Int = 2
end

config_type(::Val{:roberta}) = HGFRobertaConfig


text:
using Revise

using Devito

configuration!("log-level", "DEBUG")
configuration!("language", "openmp")
configuration!("mpi", false)

x = SpaceDimension(name="x", spacing=Constant(name="h_x", value=5.0))
z = SpaceDimension(name="z", spacing=Constant(name="h_z", value=5.0))

function ricker!(src, f, _t, t₀)
    t = reshape(_t, size(src))
    src .= (1.0 .- 2.0 * (pi * f * (t .- t₀)).^2) .* exp.(-(pi * f * (t .- t₀)).^2)
end

grid = Grid(
    dimensions = (z,x),
    shape = (501,251),
    origin = (0.0,0.0),
    extent = (2500.0,1250.0),
    dtype = Float32)

b = Devito.Function(name="b", grid=grid, space_order=8)
v = Devito.Function(name="v", grid=grid, space_order=8)
q = Devito.Function(name="woverq", grid=grid, space_order=8)

b_data = data(b)
v_data = data(v)
q_data = data(q)

b_data .= 1
v_data .= 1.5
q_data .= 1/1000

time_range = 0.0f0:0.5f0:1000.0f0

p = TimeFunction(name="p", grid=grid, time_order=2, space_order=8)
z,x,t = dimensions(p)

src = SparseTimeFunction(name="src", grid=grid, npoint=1, nt=length(time_range))
src_coords = coordinates(src)
src_coords .= [625.0; 5.0]
src_data = data(src)
ricker!(src_data, 0.01, collect(time_range), 125)
src_term = inject(src; field=forward(p), expr=src*spacing(t)^2*v^2/b)

nz,nx,δz,δx = size(grid)...,spacing(grid)...
rec = SparseTimeFunction(name="rec", grid=grid, npoint=nx, nt=length(time_range))
rec_coords = coordinates(rec)
rec_coords[1,:] .= δx*(0:nx-1)
rec_coords[2,:] .= 10.0

rec_term = interpolate(rec, expr=p)

g1(field) = dx(field,x0=x+spacing(x)/2)
g3(field) = dz(field,x0=z+spacing(z)/2)
g1_tilde(field) = dx(field,x0=x-spacing(x)/2)
g3_tilde(field) = dz(field,x0=z-spacing(z)/2)

update_p = spacing(t)^2 * v^2 / b *
    (g1_tilde(b * g1(p)) + g3_tilde(b * g3(p))) +
    (2 - spacing(t) * q) * p +
    (spacing(t) * q - 1) * backward(p)

stencil_p = Eq(forward(p), update_p)

dt = step(time_range)
smap = spacing_map(grid)
smap[spacing(t)] = dt

op = Operator([stencil_p, src_term, rec_term], subs=smap, name="OpExampleIso")

summary = apply(op)
@show summary

_p = data(p)
d = data(rec)

using PyPlot

_p = data(p)

d = data(rec)

figure();imshow(_p[:,:,1])
display(gcf())

figure();imshow(d, aspect="auto",cmap="gray",clim=.1*[-1,1]*maximum(abs,d))
display(gcf())


text:
using Dates
using Test

using Stonks
using Stonks: JSONContent, APIResponseError, ContentParserError
using Stonks.Models
using Stonks.Parsers

include("test_utils.jl")

@testset "Alphavantage parsers" begin
  price_content = get_test_data("data/alphavantage_prices.json")
  exchange_content = get_test_data("data/alphavantage_exchange.json")
  is_content = get_test_data("data/alphavantage_income_statement.json")
  bs_content = get_test_data("data/alphavantage_balance_sheet.json")
  cf_content = get_test_data("data/alphavantage_cashflow.json")
  eps_content = get_test_data("data/alphavantage_earnings.json")
  price_parser = Parsers.AlphavantagePriceParser
  info_parser = Parsers.AlphavantageInfoParser
  exchange_parser = Parsers.AlphavantageExchangeRateParser
  is_parser = Parsers.AlphavantageIncomeStatementParser
  bs_parser = Parsers.AlphavantageBalanceSheetParser
  cf_parser = Parsers.AlphavantageCashflowParser
  eps_parser = Parsers.AlphavantageEarningsParser

  @testset "Succesful price response" begin
    data = parse_content(price_parser, price_content)
    dates = map(x -> x.date, data)
    @test isa(data, Vector{AssetPrice})
    @test minimum(dates) == Date("2022-02-15")
    @test maximum(dates) == Date("2022-02-18")
  end

  @testset "Succesful price response with from date limit" begin
    from = Date("2022-02-16")
    data = parse_content(price_parser, price_content; from=from)
    dates = map(x -> x.date, data)
    @test minimum(dates) == from
    @test maximum(dates) == Date("2022-02-18")
  end

  @testset "Succesful price response with from and to date limits" begin
    data = parse_content(
      price_parser, price_content; from=Date("2022-02-16"), to=Date("2022-02-17")
    )
    dates = map(x -> x.date, data)
    @test minimum(dates) == Date("2022-02-16")
    @test maximum(dates) == Date("2022-02-17")
  end

  @testset "Error response when API key is missing or invalid" begin
    content = """{
      "Error Message": "the parameter apikey is invalid or missing.
       Please claim your free API key on (https://www.alphavantage.co/support/#api-key).
       It should take less than 20 seconds."
    }""" |> c-> replace(c, "\n" => "")
    data = parse_content(price_parser, content)
    @test isa(data, APIResponseError)
  end

  @testset "Succesful overview (info) response" begin
    overview_content = get_test_data("data/alphavantage_overview.json")
    data = parse_content(info_parser, overview_content)
    @test length(data) == 1
    @test first(data).symbol == "IBM"
  end

  @testset "Successful exchange rate response" begin
    data = parse_content(exchange_parser, exchange_content)
    @test unique(map(x -> (x.base, x.target), data)) == [("EUR", "USD")]
    @test isa(data, Vector{ExchangeRate})
    @test length(data) == 5
  end

  @testset "Successful exchange rate response with from limit" begin
    from = Date("2022-02-18")
    data = parse_content(exchange_parser, exchange_content; from=from)
    @test minimum(map(x -> x.date, data)) == from
  end

  @testset "Successful exchange rate response with from and to limits" begin
    from, to = Date("2022-02-18"), Date("2022-02-21")
    data = parse_content(exchange_parser, exchange_content; from=from, to=to)
    dates = map(x -> x.date, data)
    @test minimum(dates) == from
    @test maximum(dates) == to
  end

  @testset "Successful income statement response" begin
    data = parse_content(is_parser, is_content)
    dates = map(x -> x.date, data)
    @test isa(data, Vector{Models.IncomeStatement})
    @test length(data) == 25
    @test first(data).symbol == "IBM"
  end

  @testset "Successful income statement response with yearly frequency" begin
    data = parse_content(is_parser, is_content; frequency="yearly")
    @test length(data) == 5
  end

  @testset "Successful income statement response with quarterly frequency" begin
    data = parse_content(is_parser, is_content; frequency="quarterly")
    @test length(data) == 20
  end

  @testset "Successful income statement response with date filters" begin
    date_min, date_max = Date("2020-01-01"), Date("2020-12-31")
    data = parse_content(is_parser, is_content; from=date_min, to=date_max)
    dates = map(x -> x.date, data)
    @test minimum(dates) >= date_min && maximum(dates) <= date_max
  end

  @testset "Empty income statement response" begin
    res = parse_content(is_parser, "{}")
    @test isa(res, ContentParserError)
  end

  @testset "Successful balance sheet response" begin
    data = parse_content(bs_parser, bs_content)
    dates = map(x -> x.date, data)
    @test isa(data, Vector{Models.BalanceSheet})
    @test length(data) == 25
    @test first(data).symbol == "IBM"
  end

  @testset "Successful balance sheet response with quarterly frequency" begin
    data = parse_content(bs_parser, bs_content; frequency="quarterly")
    @test length(data) == 20
  end

  @testset "Successful balance sheet response with date filters" begin
    date_min, date_max = Date("2020-01-01"), Date("2020-12-31")
    data = parse_content(bs_parser, bs_content; from=date_min, to=date_max)
    dates = map(x -> x.date, data)
    @test minimum(dates) >= date_min && maximum(dates) <= date_max
  end

  @testset "Empty balance sheet response" begin
    res = parse_content(bs_parser, "{}")
    @test isa(res, ContentParserError)
  end

  @testset "Successful earnings response" begin
    data = parse_content(eps_parser, eps_content)
    dates = map(x -> x.date, data)
    @test isa(data, Vector{Models.Earnings})
    @test length(data) == 15
    @test first(data).symbol == "IBM"
  end

  @testset "Successful earnings response with quarterly frequency" begin
    data = parse_content(eps_parser, eps_content; frequency="quarterly")
    @test length(data) == 12
    @test unique(map(x -> x.frequency, data)) == ["quarterly"]
  end

  @testset "Successful earnings response with date filters" begin
    date_min, date_max = Date("2020-01-01"), Date("2020-12-31")
    data = parse_content(eps_parser, eps_content; from=date_min, to=date_max)
    dates = map(x -> x.date, data)
    @test minimum(dates) >= date_min && maximum(dates) <= date_max
  end

  @testset "Successful cashflow statement response" begin
    data = parse_content(cf_parser, cf_content)
    dates = map(x -> x.date, data)
    @test isa(data, Vector{Models.CashflowStatement})
    @test length(data) == 11
    @test first(data).symbol == "IBM"
  end

  @testset "Successful cash flow response with quarterly frequency" begin
    data = parse_content(cf_parser, cf_content; frequency="quarterly")
    @test length(data) == 8
    @test unique(map(x -> x.frequency, data)) == ["quarterly"]
  end

  @testset "Successful cash flow response with date filters" begin
    date_min, date_max = Date("2020-01-01"), Date("2020-12-31")
    data = parse_content(cf_parser, cf_content; from=date_min, to=date_max)
    dates = map(x -> x.date, data)
    @test minimum(dates) >= date_min && maximum(dates) <= date_max
  end

end


text:
using StochasticDiffEq, DiffEqNoiseProcess

f(du,u,p,t) = (du .= u)
g(du,u,p,t) = (du .= u)
u0 = rand(4,2)

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRIW1())
sol = solve(prob,SRIW1())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRI())
sol = solve(prob,SRI())

W = WienerProcess(0.0,0.0,0.0,reset=false)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRI())
@test_throws ErrorException sol = solve(prob,SRI())

g(du,u,p,t) = (du .= 1)

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRA1())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRA2())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRA3())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SOSRA())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SOSRA2())

W = WienerProcess(0.0,0.0,0.0)
prob = SDEProblem(f,g,u0,(0.0,1.0),noise=W)
sol = solve(prob,SRA())


text:
abstract type AbstractPacket end
abstract type AbstractRegister end
mutable struct Packet <: AbstractPacket
    val::Float64
end

function __repr__(self::AbstractPacket)::String
    return AbstractPacket(self.val)
end

function __eq__(self::AbstractPacket, other::AbstractPacket)::Bool
    return __key(self) == __key(other)
end

function __key(self::AbstractPacket)
    (self.val)
end

mutable struct Register <: AbstractRegister
    PACKET::AbstractPacket
    VALUE::Int64
end

if abspath(PROGRAM_FILE) == @__FILE__
    a = VALUE(Register, 10)
    @assert(is_value(a))
    value(a)
    b = PACKET(Register, Packet(1.3))
    @assert(is_packet(b))
    packet(b)
    println("OK")
end


text:
function deduce_mapping_from_counts(input_line, pattern_dict, possible_mappings)
    # if the input line contains the patterns for
    # "1", "2", "3", only, we should expect the following
    # counts of the true segments:
    # a: 8 b: 6 c: 8 d: 7 e: 4 f: 9 g: 7
    for encoded_segment in collect("abcdefg")
        if count(x->(x==encoded_segment), input_line) == 4
            possible_mappings[encoded_segment] = "e"
        elseif count(x->(x==encoded_segment), input_line) == 6
            possible_mappings[encoded_segment] = "b"
        elseif count(x->(x==encoded_segment), input_line) == 7
            possible_mappings[encoded_segment] = "dg"
        elseif count(x->(x==encoded_segment), input_line) == 8
            possible_mappings[encoded_segment] = "ac"
        elseif count(x->(x==encoded_segment), input_line) == 9
            possible_mappings[encoded_segment] = "f"
        end
    end
    return possible_mappings
end

function find_intersections(digit, true_pattern, possible_mappings)
    # modify a dictionary such that the value for a given key ('digit') is the
    # intersection between the current value and a new string
    for d in digit
        possible_mappings[d] = join(intersect(possible_mappings[d], true_pattern))
    end
    return possible_mappings
end

function deduce_mapping_from_combinations(input_line, pattern_dict, possible_mappings)
    # just look at 1, 4, 7 and narrow down possible mappings according to what segments
    # are in each of them.
    for digit in split(input_line)
        # sort it alphabetically
        digit = join(sort(collect(digit)))
        if length(digit) == 2  # it's a "1"
            possible_mappings = find_intersections(digit, pattern_dict["1"], possible_mappings)
        elseif length(digit) == 3 # it's a "7"
            possible_mappings = find_intersections(digit, pattern_dict["7"], possible_mappings)
        elseif length(digit) == 4 # it's a "4"
            possible_mappings = find_intersections(digit, pattern_dict["4"], possible_mappings)
        end
    end
    return possible_mappings
end

function cleanup_mappings(possible_mappings)
    ## if we already have enough information to unambiguously decode, but
    ## haven't yet propagated it through the map, do that here.
    for k in keys(possible_mappings)
        if length(possible_mappings[k]) == 1
            continue
        end
        newstring = ""
        for c in possible_mappings[k]
            if (!).(string(c) in values(possible_mappings))
                newstring *= c
            end
        end
        possible_mappings[k] = newstring
    end
    return possible_mappings
end

function deduce_mapping(input_line, pattern_dict)
    # create an initial dict of encoded_segment -> possible_plaintexts
    # that we will then try to disambiguate using two different methods
    possible_mappings_init = Dict('a'=>"abcdefg",
                                  'b'=>"abcdefg",
                                  'c'=>"abcdefg",
                                  'd'=>"abcdefg",
                                  'e'=>"abcdefg",
                                  'f'=>"abcdefg",
                                  'g'=>"abcdefg")
    possible_mappings = deduce_mapping_from_counts(input_line,
                                                   pattern_dict,
                                                   possible_mappings_init)
    possible_mappings = deduce_mapping_from_combinations(input_line,
                                                         pattern_dict,
                                                         possible_mappings)
    # should have all the info there, but need to propagate through the dict
    possible_mappings = cleanup_mappings(possible_mappings)
    return possible_mappings
end

function get_value_string(encoded_digits, reverse_patterns_dict, mapping)
    # use the mapping, and the reverse digit->segments dict, to decode
    # the four numbers in an encoded output.
    output_string = ""
    for encoded_digit in split(encoded_digits)
        pattern_key = ""
        for d in encoded_digit
            pattern_key *= mapping[d]
        end
        # sort into alphabetical order
        pattern_key = join(sort(collect(pattern_key)))
        number = reverse_patterns_dict[pattern_key]
        println("found number ",number)
        output_string *= number
    end
    return output_string
end

function process_row(signal_pattern, encoded_output, pattern_dict)
    # put everything together for a single signal pattern and output
    mapping = deduce_mapping(signal_pattern, pattern_dict)
    patterns_reverse = Dict(value=>key for (key, value) in pattern_dict)
    output_string = get_value_string(encoded_output, patterns_reverse, mapping)
    return parse(Int, output_string)
end

patterns = Dict("0"=>"abcefg",
                "1"=>"cf",
                "2"=>"acdeg",
                "3"=>"acdfg",
                "4"=>"bcdf",
                "5"=>"abdfg",
                "6"=>"abdefg",
                "7"=>"acf",
                "8"=>"abcdefg",
                "9"=>"abcdfg")

# parse the input into strings
input_lines = readlines("input.txt")

# now go through each line (and pray for the best!)
total = 0
for i in 1:length(input_lines)
    signal_pattern, output = split(input_lines[i],"|")
    global total += process_row(signal_pattern, output, patterns)
end

println("total is ", total)


text:
using RobotZoo
using RobotDynamics
using Test
using StaticArrays
using BenchmarkTools
using Rotations
using Random

@testset "Allocations" begin
    include("models_test.jl")
end

@testset "LinearModels" begin
    include("linear_models_test.jl")
    include("linear_models_examples_test.jl")
end


text:
function compare_assert(a::Int64, b::Int64)
    @assert(a == b)
    @assert(!(0 == 1))
end

function main()
    @assert(true)
    @assert(!(false))
    compare_assert(1, 1)
    @assert(true)
    @assert(true)
    println(join(["OK"], " "))
end

main()


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

const line_width = 60

const alu = string(
   "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG",
   "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA",
   "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT",
   "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA",
   "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG",
   "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC",
   "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")

const iub1 = b"acgtBDHKMNRSVWY"
const iub2 = [0.27, 0.12, 0.12, 0.27, 0.02,0.02, 0.02, 0.02, 0.02, 0.02,0.02, 0.02, 0.02, 0.02, 0.02]

const homosapiens1 = b"acgt"
const homosapiens2 = [0.3029549426680, 0.1979883004921,0.1975473066391, 0.3015094502008]

const IM   = Int32(139968)
const IA   = Int32(3877)
const IC   = Int32(29573)
const state = Ref(Int32(42))
gen_random() = (state[] = ((state[] * IA) + IC) % IM)

function repeat_fasta(io, src, n)
    k = length(src)
    s = string(src, src, src[1:(n % k)])
    I = Iterators.cycle(src)
    col = 1
    count = 1
    c, state = iterate(I)
    write(io, c % UInt8)
    while count < n
        col += 1
        c, state = iterate(I, state)
        write(io, c % UInt8)
        if col == line_width
            write(io, '\n')
            col = 0
        end
        count += 1
    end
    write(io, '\n')
    return
end

function choose_char(cs)
    k = length(cs)
    r = gen_random() / IM
    r < cs[1] && return 1
    a = 1
    b = k
    while b > a + 1
        c = fld(a + b, 2)
        if r < cs[c]
            b = c
        else
            a = c
        end
    end
    return b
end

function random_fasta(io, symb, pr, n)
    cs = cumsum(pr)
    k = n
    while k > 0
        m = min(k, line_width)
        @inbounds for i = 1:m
            write(io, symb[choose_char(cs)])
        end
        write(io, '\n'%UInt8)
        k -= line_width
    end
    return
end

function perf_fasta(_n=25000000, io = stdout)
  write(io, ">ONE Homo sapiens alu\n")
  repeat_fasta(io, alu, 2n)

  write(io, ">TWO IUB ambiguity codes\n")
  random_fasta(io, iub1, iub2, 3n)
  write(io, ">THREE Homo sapiens frequency\n")
  random_fasta(io, homosapiens1, homosapiens2, 5n)
end

n = parse(Int,ARGS[1])
perf_fasta(n)


text:
#
# Record field.
#

FieldQuery(dom::Domain, pos::Int) =
    let fs = fields(dom)
        @assert 1 <= pos <= length(fs)
        Query(
            FieldSig(pos),
            Input(dom),
            fs[pos])
    end

FieldQuery(oty::Output, pos) =
    FieldQuery(domain(oty), pos)

immutable FieldSig <: AbstractPrimitive
    pos::Int
end

describe(io::IO, sig::FieldSig) = print(io, "_$(sig.pos)")

ev(sig::FieldSig, dv::DataVector) =
    column(dv, sig.pos)

ev(sig::FieldSig, vs::AbstractVector) =
    [v[sig.pos] for v in vs]

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
# News Boy Problem solved by L-Shaped Benders Decomposition
# Author: Edward J. Xu, edxu96@outlook.com
# Date: April 5th, 2019
# 0. Set working path, install packages and use module -----------------------------------------------------------------
push!(LOAD_PATH, "$(homedir())/Desktop/News Boy Problem, Stochastic Programming")
cd("$(homedir())/Desktop/News Boy Problem, Stochastic Programming")
using BendersLshapedStochasMilp_EDXU
using JuMP
using GLPKMathProgInterface
# 1. Date Input --------------------------------------------------------------------------------------------------------
demand = [12, 14, 16, 18, 20, 22, 24, 26, 28, 30] # Demand of newspapers in each scenario
lengthS = length(demand)
vec_prob = [0.05, 0.10, 0.10, 0.10, 0.15, 0.15, 0.10, 0.10, 0.10, 0.05 ] # probability of scenario
c = 20 # purchase price
p = 70 # selling price
h = 10 # scrap value
y = 20
# 2. Data Transfer -----------------------------------------------------------------------------------------------------
n_x = 10  # Number of all the x
n_y = 1
vec_min_y = hcat([0])
vec_max_y = hcat([30])
vec_f = hcat(c - 10)
vec_pi = hcat(vec_prob)
mat_c = zeros(10, 1, 1)
for i = 1:10
    mat_c[i, :, :] = hcat([- 70 + 10])
end
mat_h = zeros(10, 2, 1)
for i = 1: 10
    mat_h[i, 1, :] = hcat(- demand[i])
    mat_h[i, 2, :] = hcat(0)
end
mat3_t = zeros(10, 2, 1)
for i = 1:10
    mat3_t[i, 1, :] = hcat(0)
    mat3_t[i, 2, :] = hcat(1)
end
mat3_w = - ones(10, 2, 1)
epsilon = 0.001
timesIterationMax = 500
# 3. Begin Optimization ------------------------------------------------------------------------------------------------
BendersLshaped(; n_x = n_x, n_y = n_y, vec_min_y = vec_min_y, vec_max_y = vec_max_y, vec_f = vec_f,
    vec_pi = vec_pi, mat_c = mat_c, mat_h = mat_h, mat3_t = mat3_t, mat3_w = mat3_w,
    epsilon = epsilon, timesIterationMax = timesIterationMax)
# answer: obj = - 976


text:
module TestBenchDictHistogram

using ConcurrentCollections
using ConcurrentCollectionsBenchmarks.BenchDictHistogram:
    default_ntasks_list, generate, hist_parallel!, hist_parallel_dac, hist_seq!
using Test

function diffvalues(actual, expected)
    diffs = []
    for (key, ve) in expected
        va = actual[key]
        if ve != va
            push!(diffs, (; key, actual = va, expected = ve))
        end
    end
    return diffs
end

function test()
    datasize_list = [10, 2^5, 2^10, 2^20]
    fulldata = generate(datasize = datasize_list[end])
    @testset for datasize in datasize_list
        data = view(fulldata, 1:datasize)
        test(data)
    end
end

function test(data)
    dbase = hist_seq!(Dict{String,Int}(), data)
    @testset "seq" begin
        cdseq = hist_seq!(ConcurrentDict{String,Int}(), data)
        @test sort(collect(setdiff(keys(dbase), keys(cdseq)))) == []
        @test sort(collect(setdiff(keys(cdseq), keys(dbase)))) == []
        @test diffvalues(cdseq, dbase) == []
        @test Dict(cdseq) == dbase
    end
    @testset for ntasks in default_ntasks_list()
        cdpar = hist_parallel!(ConcurrentDict{String,Int}(), data; ntasks = ntasks)
        #=
        if Dict(cdpar) != dbase
            global FAILED = (; cdpar, dbase)
        end
        =#
        @test sort(collect(setdiff(keys(dbase), keys(cdpar)))) == []
        @test sort(collect(setdiff(keys(cdpar), keys(dbase)))) == []
        @test diffvalues(cdpar, dbase) == []
        @test Dict(cdpar) == dbase
    end
    @testset "dac" begin
        @testset for ntasks in default_ntasks_list()
            dpar = hist_parallel_dac(data; ntasks = ntasks)
            @test sort(collect(setdiff(keys(dbase), keys(dpar)))) == []
            @test sort(collect(setdiff(keys(dpar), keys(dbase)))) == []
            @test diffvalues(dpar, dbase) == []
            @test dpar == dbase
        end
    end
end

end  # module


text:
s1 = StreamingStats(Float64)
s2 = StreamingStats(10,Float64)

ntestelements = 1000
x1 = randn(ntestelements)
x2 = randn(10,ntestelements)
m1,v1 = streamingstats(x1)
m2,v2 = streamingstats(x2)

for i=1:ntestelements
    update!(s1,x1[i])
    update!(s2,x2[:,i])
    #test the results of the update function against ground truth
    @test_approx_eq_eps mean(s1) mean(x1[1:i]) 1e-10
    @test_approx_eq_eps var(s1) var(x1[1:i]) 1e-10
    @test_approx_eq_eps std(s1) std(x1[1:i]) 1e-10
    @test_approx_eq_eps mean(s2) mean(x2[:,1:i],2) 1e-10
    @test_approx_eq_eps var(s2) var(x2[:,1:i],2) 1e-10
    @test_approx_eq_eps std(s2) std(x2[:,1:i],2) 1e-10
    #test the results of the streamstats function
    @test_approx_eq_eps m1[i] mean(s1) 1e-10
    @test_approx_eq_eps v1[i] var(s1) 1e-10
    @test_approx_eq_eps m2[:,i] mean(s2) 1e-10
    @test_approx_eq_eps v2[:,i] var(s2) 1e-10
end

println("==================")
println("Test show function")
println("==================")
show(s1)
show(s2)
println("=============")
println("End test show")
println("=============")


text:
module Simulate_YouBot

using  Modia3D
import Modia3D.ModiaMath
using  Modia3D.ModiaMath.Unitful

vmat1 = Modia3D.Material(color="LightBlue" , transparency=0.0)    # material of SolidFileMesh
vmat2 = Modia3D.Material(color="Grey"      , transparency=0.3)    # material of work piece
vmat3 = Modia3D.Material(color="LightBlue" , transparency=0.3)    # material of table

workpiece_Lx = 0.04
workpiece_Ly = 0.04
workpiece_Lz = 0.08

xPosTable = 0.71
tableX    = 0.4
tableY    = 0.3
tableZ    = 0.01
heigthLeg = 0.35
widthLeg  = 0.02
@assembly Table(world) begin
    solidPlate = Solid(SolidBox(tableX, tableY, tableZ, rsmall=0.0) , "DryWood", vmat1, contactMaterial="DryWood")
    plate = Object3D(world, solidPlate, fixed=true, r=[xPosTable , 0.0, heigthLeg], visualizeFrame=false)
    leg   = Solid(SolidBox(widthLeg, widthLeg, heigthLeg) , "DryWood", vmat3)
    leg1  = Object3D(plate, leg, fixed=true, r=[tableX/2 - widthLeg/2 , tableY/2 - widthLeg/2, -heigthLeg/2])
    leg2  = Object3D(plate, leg, fixed=true, r=[-tableX/2 + widthLeg/2 , tableY/2 - widthLeg/2, -heigthLeg/2])
    leg3  = Object3D(plate, leg, fixed=true, r=[tableX/2 - widthLeg/2 , -tableY/2 + widthLeg/2, -heigthLeg/2])
    leg4  = Object3D(plate, leg, fixed=true, r=[-tableX/2 + widthLeg/2 , -tableY/2 + widthLeg/2, -heigthLeg/2])
end

@assembly YouBot begin
    world           = Object3D(visualizeFrame=true)

    # Table + workpiece
    table     = Table(world)
    workpiece = Modia3D.ContactBox(table.plate, scale=[workpiece_Lx, workpiece_Ly, workpiece_Lz], massProperties="DryWood",
                                   material=vmat2, contactMaterial="DryWood",
                                   r=[-tableX/2+1.2*workpiece_Lx/2, 0.0, workpiece_Lz/2+tableZ/2], visualizeFrame=false, fixed=false)
end

gravField = UniformGravityField(g=9.81, n=[0,0,-1])
youBot    = YouBot(sceneOptions=SceneOptions(gravityField=gravField, visualizeFrames=false,
                                             defaultFrameLength=0.1, enableContactDetection=true, elasticContactReductionFactor=1e-4))

#Modia3D.visualizeAssembly!(youBot)
model  = Modia3D.SimulationModel( youBot )
result = ModiaMath.simulate!(model, stopTime=1.0, log=false, interval=0.001, tolerance=1e-5)

ModiaMath.plot(result, "workpiece.box.r[3]", figure=1)

println("... success of Simulate_YouBotBoxOnTable.jl!")

end


text:
using CUDA: CUDA, @cuda
using DataFrames: DataFrame
using LoopVectorization: @avx

export satisfies, satisfies_all, satisfies_flows_down, satisfies_flows,
count_downflow, downflow_all

#####################
# Circuit logical evaluation
#####################

# evaluate a circuit as a function
function (root::LogicCircuit)(data...)
    satisfies(root, data...)
end

"Evaluate satisfaction of the circuit bottom-up for a given input"
satisfies(root::LogicCircuit, data::Real...) =
    satisfies(root, collect(data))

satisfies(root::LogicCircuit, data::AbstractVector{Bool}) =
    satisfies(root, DataFrame(reshape(BitVector(data), 1, :)))[1]

satisfies(root::LogicCircuit, data::AbstractVector{<:AbstractFloat}) =
    satisfies(root, DataFrame(reshape(data, 1, :)))[1]

satisfies(circuit::LogicCircuit, data::DataFrame) =
    satisfies(same_device(BitCircuit(circuit, data), data) , data)

function satisfies(circuit::BitCircuit, data::DataFrame)::AbstractVector
    output = satisfies_all(circuit,data)[:,end]
    if isbinarydata(data)
        return AbstractBitVector(output, num_examples(data))
    else
        @assert isfpdata(data) "Only floating point and binary flows are supported"
        return output
    end
end

#####################
# Circuit evaluation of *all* nodes in circuit
#####################

"Evaluate the circuit bottom-up for a given input and return the value of all nodes"
function satisfies_all(circuit::BitCircuit, data, reuse=nothing)
    @assert num_features(data) == num_features(circuit)
    @assert iscomplete(data)
    values = init_satisfies(data, reuse, num_nodes(circuit))
    satisfies_layers(circuit, values)
    return values
end

"Initialize values from the data (data frames)"
function init_satisfies(data, reuse, num_nodes)
    if isbinarydata(data)
        flowtype = isgpu(data) ? CuMatrix{UInt} : Matrix{UInt}
        values = similar!(reuse, flowtype, num_chunks(data), num_nodes)

        @views values[:,TRUE_BITS] .= typemax(UInt)
        @views values[:,FALSE_BITS] .= typemin(UInt)
        nfeatures = num_features(data)
        for i=1:nfeatures
            @views values[:,2+i] .= chunks(data,i)
            @views values[:,2+nfeatures+i] .= .~ chunks(data,i)
        end
        # warning: there are now some 1 bits beyond the BitVector mask.
        # therefore, reset those to 0
        @views values[end,1:2+2*num_features(data)] .&= _msk_end(data)

    else
        @assert isfpdata(data) "Only floating point and binary flows are supported"
        flowtype = isgpu(data) ? CuMatrix{eltype(data)} : Matrix{eltype(data)}
        values = similar!(reuse, flowtype, num_examples(data), num_nodes)
        @views values[:,TRUE_BITS] .= one(eltype(values))
        @views values[:,FALSE_BITS] .= zero(eltype(values))
        nfeatures = num_features(data)
        for i=1:nfeatures
            @views values[:,2+i] .= feature_values(data,i)
            @views values[:,2+nfeatures+i] .= one(eltype(values)) .- feature_values(data,i)
        end

    end
    return values
end

# upward pass helpers on CPU

"Evaluate the layers of a bit circuit on the CPU (SIMD & multi-threaded)"
function satisfies_layers(circuit::BitCircuit, values::Matrix)
    els = circuit.elements
    for layer in circuit.layers[2:end]
        Threads.@threads for dec_id in layer
            j = @inbounds circuit.nodes[1,dec_id]
            els_end = @inbounds circuit.nodes[2,dec_id]
            if j == els_end
                assign_value(values, dec_id, els[2,j], els[3,j])
                j += 1
            else
                assign_value(values, dec_id, els[2,j], els[3,j], els[2,j+1], els[3,j+1])
                j += 2
            end
            while j <= els_end
                if j == els_end
                    accum_value(values, dec_id, els[2,j], els[3,j])
                    j += 1
                else
                    accum_value(values, dec_id, els[2,j], els[3,j], els[2,j+1], els[3,j+1])
                    j += 2
                end
            end
        end
    end
end

assign_value(v::Matrix{<:AbstractFloat}, i, e1p, e1s) =
    @views @. @avx v[:,i] = v[:,e1p] * v[:,e1s]

accum_value(v::Matrix{<:AbstractFloat}, i, e1p, e1s) =
    @views @. v[:,i] += v[:,e1p] * v[:,e1s] # adding @avx crashes macro

assign_value(v::Matrix{<:AbstractFloat}, i, e1p, e1s, e2p, e2s) =
    @views @. @avx v[:,i] = v[:,e1p] * v[:,e1s] + v[:,e2p] * v[:,e2s]

accum_value(v::Matrix{<:AbstractFloat}, i, e1p, e1s, e2p, e2s) =
    @views @. v[:,i] += v[:,e1p] * v[:,e1s] + v[:,e2p] * v[:,e2s] # adding @avx crashes macro

assign_value(v::Matrix{<:Unsigned}, i, e1p, e1s) =
    @views @. @avx v[:,i] = v[:,e1p] & v[:,e1s]

accum_value(v::Matrix{<:Unsigned}, i, e1p, e1s) =
    @views @. v[:,i] |= v[:,e1p] & v[:,e1s] # adding @avx crashes macro

assign_value(v::Matrix{<:Unsigned}, i, e1p, e1s, e2p, e2s) =
    @views @. @avx v[:,i] = v[:,e1p] & v[:,e1s] | v[:,e2p] & v[:,e2s]

accum_value(v::Matrix{<:Unsigned}, i, e1p, e1s, e2p, e2s) =
    @views @. v[:,i] |= v[:,e1p] & v[:,e1s] | v[:,e2p] & v[:,e2s] # adding @avx crashes macro

# upward pass helpers on GPU

"Evaluate the layers of a bit circuit on the GPU"
function satisfies_layers(circuit::BitCircuit, values::CuMatrix;  dec_per_thread = 8, log2_threads_per_block = 8)
    CUDA.@sync for layer in circuit.layers[2:end]
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))
        @cuda threads=num_threads blocks=num_blocks satisfies_layers_cuda(layer, circuit.nodes, circuit.elements, values)
    end
end

"assign threads to examples and decision nodes so that everything is a power of 2"
function balance_threads(num_examples, num_decisions, total_log2)
    ratio = num_examples / num_decisions
    k = ceil(Int, (log2(ratio) + total_log2)/2)
    k = min(max(0, k), total_log2)
    l = total_log2-k
    (2^k, 2^l)
end

"CUDA kernel for circuit evaluation"
function satisfies_layers_cuda(layer, nodes, elements, values)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for j = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            decision_id = @inbounds layer[i]
            k = @inbounds nodes[1,decision_id]
            els_end = @inbounds nodes[2,decision_id]
            @inbounds values[j, decision_id] =
                el_value(values[j, elements[2,k]], values[j, elements[3,k]])
            while k < els_end
                k += 1
                @inbounds accum_el_value(values, j, decision_id, values[j,
                            elements[2,k]], values[j, elements[3,k]])
            end # would loop unrolling help here as on CPU? probably not?
        end
    end
    return nothing
end

el_value(p::AbstractFloat, s) = p * s
el_value(p::Unsigned, s) = p & s

accum_el_value(values, j, decision_id, p::AbstractFloat, s) =
    @inbounds values[j, decision_id] += el_value(p, s)
accum_el_value(values, j, decision_id, p::Unsigned, s) =
    @inbounds values[j, decision_id] |= el_value(p, s)

#####################
# Bit circuit values and flows (up and downward pass)
#####################

"Compute the value and flow of each node"
function satisfies_flows(circuit::LogicCircuit, data,
                         reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop,
                         weights=nothing)
    bc = same_device(BitCircuit(circuit, data), data)
    satisfies_flows(bc, data, reuse_values, reuse_flows; on_node, on_edge, weights)
end

function satisfies_flows(circuit::BitCircuit, data,
                         reuse_values=nothing, reuse_flows=nothing; on_node=noop, on_edge=noop,
                         weights=nothing)
    @assert isgpu(data) == isgpu(circuit) "BitCircuit and data need to be on the same device"
    values = satisfies_all(circuit, data, reuse_values)
    flows = satisfies_flows_down(circuit, values, reuse_flows; on_node, on_edge, weights)
    return values, flows
end

#####################
# Bit circuit flows downward pass
#####################

"When values of nodes have already been computed, do a downward pass computing the flows at each node"
function satisfies_flows_down(circuit::BitCircuit, values, reuse=nothing; on_node=noop, on_edge=noop, weights=nothing)
    flows = similar!(reuse, typeof(values), size(values)...)
    if !(weights isa Nothing) && values isa CuMatrix
        weights = convert(CuVector{Float32}, weights)
    end
    satisfies_flows_down_layers(circuit, flows, values, on_node, on_edge; weights = weights)
    return flows
end

# downward pass helpers on CPU

"Evaluate the layers of a bit circuit on the CPU (SIMD & multi-threaded)"
function satisfies_flows_down_layers(circuit::BitCircuit, flows::Matrix, values::Matrix, on_node, on_edge; weights = nothing)
    els = circuit.elements
    for layer in Iterators.reverse(circuit.layers)
        Threads.@threads for dec_id in layer
            par_start = @inbounds circuit.nodes[3,dec_id]
            if iszero(par_start)
                if dec_id == num_nodes(circuit)
                    # populate root flow from values
                    @inbounds @views @. flows[:, dec_id] = values[:, dec_id]
                end
                # no parents, ignore (can happen for false/true node and root)
            else
                par_end = @inbounds circuit.nodes[4,dec_id]
                for j = par_start:par_end
                    par = @inbounds circuit.parents[j]
                    grandpa = @inbounds els[1,par]
                    sib_id = sibling(els, par, dec_id)
                    # TODO: check if things would speed up if we also allowed the first parent to write flows in lower-down layers
                    single_child = has_single_child(circuit.nodes, grandpa)
                    if single_child
                        if j == par_start
                            @inbounds @views @. flows[:, dec_id] = flows[:, grandpa]
                        else
                            accum_flow(flows, dec_id, grandpa)
                        end
                    else
                        if j == par_start
                            assign_flow(flows, values, dec_id, grandpa, sib_id)
                        else
                            accum_flow(flows, values, dec_id, grandpa, sib_id)
                        end
                    end
                    # report edge flow only once:
                    sib_id > dec_id && on_edge(flows, values, dec_id, sib_id, par, grandpa, single_child, weights)
                end
            end
            on_node(flows, values, dec_id, weights)
        end
    end
end

assign_flow(f::Matrix{<:Unsigned}, v, d, g, s) =
    @inbounds @views @. f[:, d] = v[:, s] & v[:, d] & f[:, g]

function assign_flow(f::Matrix{<:AbstractFloat}, v, d, g, s)
    @avx for j in 1:size(f,1)
        edge_flow = v[j, s] * v[j, d] / v[j, g] * f[j, g]
        edge_flow = vifelse(isfinite(edge_flow), edge_flow, zero(eltype(f)))
        f[j, d] = edge_flow
    end
end

accum_flow(f::Matrix{<:Unsigned}, d, g) =
    @inbounds @views @. f[:, d] |= f[:, g]

accum_flow(f::Matrix{<:AbstractFloat}, d, g) =
    @inbounds @views @. f[:, d] += f[:, g]

accum_flow(f::Matrix{<:Unsigned}, v, d, g, s) =
    @inbounds @views @. f[:, d] |= v[:, s] & v[:, d] & f[:, g]

function accum_flow(f::Matrix{<:AbstractFloat}, v, d, g, s)
    @avx for j in 1:size(f,1)
        edge_flow = v[j, s] * v[j, d] / v[j, g] * f[j, g]
        edge_flow = vifelse(isfinite(edge_flow), edge_flow, zero(eltype(f)))
        f[j, d] += edge_flow
    end
end

# # downward pass helpers on GPU

"Pass flows down the layers of a bit circuit on the GPU"
function satisfies_flows_down_layers(circuit::BitCircuit, flows::CuMatrix, values::CuMatrix, on_node, on_edge;
            dec_per_thread = 8, log2_threads_per_block = 7, weights = nothing)
    CUDA.@sync for layer in Iterators.reverse(circuit.layers)
        num_examples = size(values, 1)
        num_decision_sets = length(layer)/dec_per_thread
        num_threads =  balance_threads(num_examples, num_decision_sets, log2_threads_per_block)
        num_blocks = (ceil(Int, num_examples/num_threads[1]),
                      ceil(Int, num_decision_sets/num_threads[2]))
        @cuda threads=num_threads blocks=num_blocks satisfies_flows_down_layers_cuda(layer, circuit.nodes, circuit.elements, circuit.parents, flows, values, on_node, on_edge, weights)
    end
end

"CUDA kernel for passing flows down circuit"
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values, on_node, on_edge, weights::Nothing)
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, nothing)
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, nothing)
        end
    end
    return nothing
end
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values, on_node, on_edge, weights::CUDA.CuDeviceArray{<:AbstractFloat,1,1})
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        start_idx = 64 * k - 63
        end_idx = CUDA.min(64 * k, CUDA.length(weights))
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            if dec_id == prime
                                for bit_idx::UInt32 = 0 : (end_idx - start_idx)
                                    weight::Float32 = @inbounds weights[start_idx + bit_idx]
                                    on_edge(flows, values, prime, sub, par, grandpa, bit_idx, edge_flow, single_child, weight)
                                end
                            end
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            for bit_idx::UInt32 = 0 : (end_idx - start_idx)
                weight::Float32 = @inbounds weights[start_idx + bit_idx]
                on_node(flows, values, dec_id, bit_idx, flow, weight)
            end
        end
    end
    return nothing
end
function satisfies_flows_down_layers_cuda(layer, nodes, elements, parents, flows, values::CUDA.CuDeviceArray{<:AbstractFloat,2,1}, on_node, on_edge, weights::CUDA.CuDeviceArray{<:AbstractFloat,1,1})
    index_x = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    index_y = (blockIdx().y - 1) * blockDim().y + threadIdx().y
    stride_x = blockDim().x * gridDim().x
    stride_y = blockDim().y * gridDim().y
    for k = index_x:stride_x:size(values,1)
        for i = index_y:stride_y:length(layer)
            dec_id = @inbounds layer[i]
            if dec_id == size(nodes,2)
                # populate root flow from values
                flow = values[k, dec_id]
            else
                par_start = @inbounds nodes[3,dec_id]
                flow = zero(eltype(flows))
                if !iszero(par_start)
                    par_end = @inbounds nodes[4,dec_id]
                    for j = par_start:par_end
                        par = @inbounds parents[j]
                        grandpa = @inbounds elements[1,par]
                        v_gp = @inbounds values[k, grandpa]
                        prime = elements[2,par]
                        sub = elements[3,par]
                        if !iszero(v_gp) # edge flow only gets reported when non-zero
                            f_gp = @inbounds flows[k, grandpa]
                            single_child = has_single_child(nodes, grandpa)
                            if single_child
                                edge_flow = f_gp
                            else
                                v_prime = @inbounds values[k, prime]
                                v_sub = @inbounds values[k, sub]
                                edge_flow = compute_edge_flow( v_prime, v_sub, v_gp, f_gp)
                            end
                            flow = sum_flows(flow, edge_flow)
                            # report edge flow only once:
                            dec_id == prime && on_edge(flows, values, prime, sub, par, grandpa, k, edge_flow, single_child, weights[k])
                        end
                    end
                end
            end
            @inbounds flows[k, dec_id] = flow
            on_node(flows, values, dec_id, k, flow, weights[k])
        end
    end
    return nothing
end

@inline sum_flows(a,b::AbstractFloat) = a + b
@inline sum_flows(a,b::Unsigned) = a | b

@inline compute_edge_flow(p_up::AbstractFloat, s_up, n_up, n_down) = p_up * s_up / n_up * n_down
@inline compute_edge_flow(p_up::Unsigned, s_up, n_up, n_down) = p_up & s_up & n_down

# API to get flows

function count_downflow(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit)
    dec_id = n.data.node_id
    sum(1:size(flows,1)) do i
        count_ones(flows[i, dec_id])
    end
end

function downflow_all(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit)
    dec_id = n.data.node_id
    indices = map(1:size(flows,1)) do i
        digits(Bool, flows[i, dec_id], base=2, pad=64)
    end
    BitArray(vcat(indices...)[1:N])
end

function count_downflow(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, c::LogicCircuit)
    grandpa = n.data.node_id
    if isleafgate(c)
        par = c.data.node_id
        return sum(1:size(flows,1)) do i
            count_ones(values[i, par] & flows[i, grandpa])
        end
    else
        ids = [x.data.node_id for x in children(c)]
        return sum(1:size(flows,1)) do i
            indices = flows[i, grandpa]
            for id in ids
                indices &= values[i, id]
            end
            count_ones(indices)
        end
    end
end

function downflow_all(values::Matrix{UInt64}, flows::Matrix{UInt64}, N, n::LogicCircuit, c::LogicCircuit)
    grandpa = n.data.node_id
    if isleafgate(c)
        par = c.data.node_id
        edge = map(1:size(flows,1)) do i
            digits(Bool, values[i, par] & flows[i, grandpa], base=2, pad=64)
        end
    else
        ids = [x.data.node_id for x in children(c)]
        edge = map(1:size(flows,1)) do i
            indices = flows[i, grandpa]
            for id in ids
                indices = indices & values[i, id]
            end
            digits(Bool, indices, base=2, pad=64)
        end
    end
    BitArray(vcat(edge...)[1:N])
end


text:
import FFTW
import Luna: SFA, Tools, Grid, Maths, PhysData, Ionisation
import Luna.PhysData: wlfreq
import PyPlot: plt, pygui
pygui(true)

τfwhm = 30e-15
λ0 = 800e-9
Ipeak = 3e18
gas = :Ar

grid = Grid.RealGrid(0, λ0, (10e-9, 1500e-9), 4*τfwhm)

It = Maths.gauss.(grid.t; fwhm=τfwhm) .* Ipeak
Et = Tools.intensity_to_field.(It) .* cos.(grid.t .* wlfreq(λ0))

irf! = Ionisation.ionrate_fun!_PPTcached(gas, λ0)

D = SFA.sfa_dipole(grid.t, Et, gas, λ0; gate=true, depletion=true, irf!)
eV = PhysData.ħ .* grid.ω ./ PhysData.electron
Dω = FFTW.rfft(D)
tg = collect(range(extrema(grid.t)...; length=2048))
gab = Maths.gabor(grid.t, D, tg, 300e-18)
##
IpeV = PhysData.ionisation_potential(gas; unit=:eV)
# Ponderomotive energy:
Up = PhysData.electron^2 * maximum(Et)^2/(4*PhysData.m_e * wlfreq(λ0)^2)
UpeV = Up/PhysData.electron
# Cutoff law:
cutoffeV = IpeV + 3.17UpeV

eV0 = PhysData.ħ*wlfreq(λ0)/PhysData.electron
harmonics = eV./eV0

##
plt.figure()
plt.semilogy(harmonics, abs2.(Dω), label="Harmonic spectrum")
plt.axvline(cutoffeV/eV0, linestyle="--", color="k", label="IP + 3.17 Up")
plt.xlim(0, 1.2cutoffeV/eV0)
plt.xlabel("Harmonic order")
plt.ylabel("Spectral energy density")
plt.legend(frameon=false)

##
plt.figure()
plt.pcolormesh(tg*1e15, eV[1:4:end], Maths.log10_norm(abs2.(gab[1:4:end, :])))
plt.axhline(cutoffeV, linestyle="--", color="w")
plt.clim(-10, -4)
plt.ylim(20, 1.2*cutoffeV)
plt.xlim(-τfwhm*1e15, τfwhm*1e15)
cb = plt.colorbar()
cb.set_label("Log10 SED")
plt.xlabel("Time (fs)")
plt.ylabel("Photon energy (eV)")


text:
function init_state(model::AbstractGPModel)
    state = init_local_vars((;), model)
    state = init_opt_state(state, model)
    state = init_hyperopt_state(state, model)
    if n_output(model) > 1
        state = init_state_A(state, model)
    end
    return state
end

@traitfn function init_local_vars(
    state, model::TGP
) where {TGP <: AbstractGPModel; !IsMultiOutput{TGP}}
    if inference(model) isa Union{Analytic,AnalyticVI,GibbsSampling}
        local_vars = init_local_vars(likelihood(model), batchsize(model))
        return merge(state, (; local_vars))
    else
        return state
    end
end

@traitfn function init_local_vars(
    state, model::TGP
) where {TGP <: AbstractGPModel; IsMultiOutput{TGP}}
    if inference(model) isa Union{Analytic,AnalyticVI,GibbsSampling}
        local_vars = init_local_vars.(likelihood(model), batchsize(model))
        return merge(state, (; local_vars))
    else
        return state
    end
end

function init_opt_state(state, model::AbstractGPModel)
    if inference(model) isa VariationalInference
        opt_state = map(model.f) do gp
            init_opt_state(gp, inference(model))
        end
        return merge(state, (; opt_state))
    else
        return state
    end
end

function init_opt_state(
    ::Union{VarLatent{T},TVarLatent{T}}, vi::VariationalInference
) where {T}
    return (;)
end

function init_opt_state(gp::Union{VarLatent{T},TVarLatent{T}}, vi::NumericalVI) where {T}
    return (;
        ν=zeros(T, batchsize(vi)), # Derivative -<dv/dx>_qn
        λ=zeros(T, batchsize(vi)), # Derivative  <d²V/dx²>_qm
        state_μ=Optimisers.state(opt(vi).optimiser, mean(gp)),
        state_Σ=Optimisers.state(opt(vi), cov(gp).data),
        ∇η₁=zero(mean(gp)),
        ∇η₂=zero(cov(gp).data),
    )
end

function init_opt_state(gp::SparseVarLatent{T}, vi::VariationalInference) where {T}
    state = (; ∇η₁=zero(mean(gp)), ∇η₂=zero(cov(gp).data))
    if vi isa AnalyticVI && is_stochastic(vi)
        state = merge(
            state,
            (;
                state_η₁=Optimisers.state(opt(vi).optimiser, nat1(gp)),
                state_η₂=Optimisers.state(opt(vi), nat2(gp).data),
            ),
        )
    end
    if vi isa NumericalVI
        state = merge(
            state,
            (;
                ν=zeros(T, batchsize(vi)), # Derivative -<dv/dx>_qn
                λ=zeros(T, batchsize(vi)), # Derivative  <d²V/dx²>_qm
                state_μ=Optimisers.state(opt(vi).optimiser, mean(gp)),
                state_Σ=Optimisers.state(opt(vi), cov(gp).data),
            ),
        )
    end
    return state
end

function init_opt_state(gp::OnlineVarLatent{T}, vi::VariationalInference) where {T}
    state = (; ∇η₁=zero(mean(gp)), ∇η₂=zero(cov(gp).data))
    if is_stochastic(vi)
        state_η₁ = state(opt(vi), nat1(gp))
        state_η₂ = state(opt(vi), nat2(gp).data)
        state = merge(state, (; state_η₁, state_η₂))
    end
    k = dim(gp)
    prev𝓛ₐ = zero(T)
    invDₐ = Symmetric(Matrix{T}(I(k)))
    prevη₁ = zeros(T, k)
    return merge(state, (; previous_gp=(; prev𝓛ₐ, invDₐ, prevη₁)))
end

function init_state_A(state, model::AbstractGPModel)
    A_state = [
        [Optimisers.init(model.A_opt, model.A[i][j]) for j in 1:model.nf_per_task[i]] for
        i in 1:n_output(model)
    ]
    return merge(state, (; A_state))
end

function init_hyperopt_state(state, model::GP)
    hyperopt_state = init_hyperopt_state(model.f)
    return merge(state, (; hyperopt_state))
end

function init_hyperopt_state(state, model::AbstractGPModel)
    hyperopt_state = map(model.f) do gp
        init_hyperopt_state(gp)
    end
    return merge(state, (; hyperopt_state))
end

@traitfn function init_hyperopt_state(gp::TGP) where {TGP <: AbstractLatent; IsFull{TGP}}
    hyperopt_state = (;)
    if !isnothing(opt(gp))
        k = kernel(gp)
        state_k = Optimisers.state(opt(gp), k)
        hyperopt_state = merge(hyperopt_state, (; state_k))
    end
    hyperopt_state = init_priormean_state(hyperopt_state, pr_mean(gp))
    return hyperopt_state
end

@traitfn function init_hyperopt_state(gp::TGP) where {TGP <: AbstractLatent; !IsFull{TGP}}
    hyperopt_state = (;)
    if !isnothing(opt(gp))
        k = kernel(gp)
        state_k = Optimisers.state(opt(gp), k)
        hyperopt_state = merge(hyperopt_state, (; state_k))
    end
    if !isnothing(Zopt(gp))
        Z = Zview(gp)
        state_Z = Optimisers.state(opt(gp), Z)
        hyperopt_state = merge(hyperopt_state, (; state_Z))
    end
    hyperopt_state = init_priormean_state(hyperopt_state, pr_mean(gp))
    return hyperopt_state
end

function Optimisers.state(opt, Z::Union{ColVecs,RowVecs})
    return Optimisers.state(opt, Z.X)
end


text:
module LightOSM

using Parameters
using DataStructures: DefaultDict, OrderedDict, MutableLinkedList, PriorityQueue, dequeue!, dequeue_pair!
using Statistics: mean
using SparseArrays: SparseMatrixCSC, sparse
using Graphs: AbstractGraph, DiGraph, nv, outneighbors, weakly_connected_components, vertices, all_neighbors, indegree, outdegree, add_edge!
using StaticGraphs: StaticDiGraph
using SimpleWeightedGraphs: SimpleWeightedDiGraph
using MetaGraphs: MetaDiGraph
using NearestNeighbors: KDTree, knn
using HTTP
using JSON
using LightXML
using DataFrames
using GeoInterface
using RecipesBase

export GeoLocation,
       AbstractOSMGraph,
       OSMGraph,
       SimplifiedOSMGraph,
       Node,
       Way,
       Restriction,
       Building,
       distance,
       heading,
       calculate_location,
       download_osm_network,
       graph_from_object,
       graph_from_download,
       graph_from_file,
       shortest_path,
       weights_from_path,
       total_path_weight,
       nearest_node,
       download_osm_buildings,
       buildings_from_object,
       buildings_from_download,
       buildings_from_file,
       simplify_graph,
       node_gdf,
       edge_gdf,
       highway_gdf

include("types.jl")
include("constants.jl")
include("utilities.jl")
include("geometry.jl")
include("download.jl")
include("parse.jl")
include("graph.jl")
include("traversal.jl")
include("shortest_path.jl")
include("nearest_node.jl")
include("buildings.jl")
include("simplification.jl")
include("geodataframes.jl")
include("plotrecipes.jl")

end # module


text:
using Test

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

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

bridged_mock = MOIB.Variable.RSOCtoPSD{Float64}(mock)

# Can come from a SOC of dimension 2 which makes more sense
# FIXME should it be moved to contconic or is RSOC of dimension 2 too exotic ?
@testset "RSOC of dimension 2" begin
    MOI.empty!(bridged_mock)
    xy, cxy = MOI.add_constrained_variables(bridged_mock, MOI.RotatedSecondOrderCone(2))
    x, y = xy
    fx = MOI.SingleVariable(x)
    fy = MOI.SingleVariable(y)
    c = MOI.add_constraint(bridged_mock, 1.0fx + 1.0fy, MOI.LessThan(1.0))
    obj = 1.0fy
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}(), obj)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0.0, 2.0],
        #(MOI.VectorOfVariables) => [0.0, 1.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0])
    MOI.optimize!(bridged_mock)
    @test MOI.get(bridged_mock, MOI.ObjectiveValue()) == 1.0
    @test MOI.get(bridged_mock, MOI.DualObjectiveValue()) == 1.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), xy) == [0.0, 1.0]
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cxy) == [0.0, 1.0]
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cxy) == [1.0, 0.0]
    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), c) == 1.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c) == -1.0

    MOI.set(mock, MOI.ConstraintName(), c, "c")
    @testset "Test mock model" begin
        var_names = ["a", "b"]
        MOI.set(mock, MOI.VariableName(),
                MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
        nonneg = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.Nonnegatives}())
        @test length(nonneg) == 1
        MOI.set(mock, MOI.ConstraintName(), nonneg[1], "cab")

        s = """
        variables: a, b
        cab: [a, b] in MathOptInterface.Nonnegatives(2)
        c: a + 0.5b <= 1.0
        maxobjective: 0.5b
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, var_names, ["cab", "c"])
    end

    @testset "Test bridged model" begin
        var_names = ["x", "y"]
        MOI.set(bridged_mock, MOI.VariableName(),
                MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
        MOI.set(bridged_mock, MOI.ConstraintName(), cxy, "cxy")

        s = """
        variables: x, y
        cxy: [x, y] in MathOptInterface.RotatedSecondOrderCone(2)
        c: x + y <= 1.0
        maxobjective: 1.0y
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["cxy", "c"])
    end

    test_delete_bridged_variables(bridged_mock, xy, MOI.RotatedSecondOrderCone, 2, (
        (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
        (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle, 0),
        (MOI.SingleVariable, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0),
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 1)
    ))
end

@testset "RSOC4" begin
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1.0, 1.0, 2.0, 1.0, 0.0, 2.0],
        (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64})  => [0.25],
        (MOI.SingleVariable, MOI.EqualTo{Float64})  => [-0.5],
        (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
        (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle) => [[1.0, -0.5, 0.25, -0.5, 0.25, 0.25]])
    mock.eval_variable_constraint_dual = false
    MOIT.rotatedsoc4test(bridged_mock, config)
    mock.eval_variable_constraint_dual = true

    @testset "Test mock model" begin
        var_names = ["Q$i$j" for j in 1:3 for i in 1:j]
        MOI.set(
            mock, MOI.VariableName(),
            MOI.get(mock, MOI.ListOfVariableIndices()), var_names)
        psd = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle}())
        @test length(psd) == 1
        MOI.set(mock, MOI.ConstraintName(), psd[1], "psd")
        off_diag = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.SingleVariable, MOI.EqualTo{Float64}}())
        @test length(off_diag) == 1
        MOI.set(mock, MOI.ConstraintName(), off_diag[1], "off_diag23")
        diag = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}())
        @test length(diag) == 1
        MOI.set(mock, MOI.ConstraintName(), diag[1], "diag33")
        c = MOI.get(mock, MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
        @test length(c) == 1
        MOI.set(mock, MOI.ConstraintName(), c[1], "c")

        s = """
        variables: Q11, Q12, Q13, Q22, Q23, Q33
        psd: [Q11, Q12, Q22, Q13, Q23, Q33] in MathOptInterface.PositiveSemidefiniteConeTriangle(3)
        off_diag23: Q23 == 0.0
        diag33: Q22 + -1.0Q33 == 0.0
        c: Q11 + 0.5Q22 <= 2.0
        maxobjective: Q12 + Q13
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(mock, model, var_names, ["psd", "off_diag23", "diag33", "c"])
    end

    @testset "Test bridged model" begin
        var_names = ["t", "u", "x", "y"]
        MOI.set(
            bridged_mock, MOI.VariableName(),
            MOI.get(bridged_mock, MOI.ListOfVariableIndices()), var_names)
        rsoc = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{
            MOI.VectorOfVariables, MOI.RotatedSecondOrderCone}())
        @test length(rsoc) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), rsoc[1], "rsoc")
        c = MOI.get(bridged_mock, MOI.ListOfConstraintIndices{
            MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}}())
        @test length(c) == 1
        MOI.set(bridged_mock, MOI.ConstraintName(), c[1], "c")

        s = """
        variables: t, u, x, y
        rsoc: [t, u, x, y] in MathOptInterface.RotatedSecondOrderCone(4)
        c: t + u <= 2.0
        maxobjective: x + y
        """
        model = MOIU.Model{Float64}()
        MOIU.loadfromstring!(model, s)
        MOIU.test_models_equal(bridged_mock, model, var_names, ["rsoc", "c"])
    end

    @testset "Delete" begin
        v = MOI.get(bridged_mock, MOI.ListOfVariableIndices())
        @test length(v) == 4

        message = string("Cannot delete variable as it is constrained with other",
                         " variables in a `MOI.VectorOfVariables`.")
        for i in 1:4
            err = MOI.DeleteNotAllowed(v[i], message)
            @test_throws err MOI.delete(bridged_mock, v[i])
        end

        test_delete_bridged_variables(bridged_mock, v, MOI.RotatedSecondOrderCone, 4, (
            (MOI.VectorOfVariables, MOI.Nonnegatives, 0),
            (MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle, 0),
            (MOI.SingleVariable, MOI.EqualTo{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}, 0)
        ))
    end
end


text:
function get_nnd(s)
    nnd = []
    clusters = unique(s)
    for c in clusters
        g_c = findall(x-> x==c, s)
        forwd = diff(g_c)
        backwd = map(x->abs(x),diff(reverse(g_c)))
        nnd_s = map((x,y)-> minimum([x,y]),forwd,backwd)
        push!(nnd,nnd_s)
    end
    return vcat(nnd...)
end
function get_H(s::Array{Int64,1})
    T = Dict{Any,Float64}()
    n = length(s)
    for i = 1:length(s)
        T[s[i]] = get(T, s[i],0) + 1/n
    end
    p = collect(values(T))
    return mapreduce(x-> -x*log(x),+,p)
end
function get_freq(s)
    T = Dict{Int64,Int64}()
    n = length(s)
    for i = 1:n
        T[s[i]] = get(T, s[i],0) + 1
    end
    return T
end
function get_nnd_distro(s,dmax::Int64)
    P = zeros(dmax)
    n = length(s)
    for i in 1:dmax
        ndata = findall(x-> x==i, s)
        if !isempty(ndata)
            P[i] = (length(ndata) + 1/dmax)/(n + 1/dmax)
        else
            P[i] = (1/dmax) / (n+1/dmax)
        end
    end
    return P
end
KS_d(x,y) = maximum(abs.(cumsum(x) - cumsum(y)))
function get_gene_cor(gene_exp)
    pcor = Float64[]
    g_di = Float64[]
    lg = size(gene_exp)[1]
    for i in 1:lg-1
        for j in i+1:lg
            #try
            push!(pcor, cor(gene_exp[i,2:end-1], gene_exp[j,2:end-1]))
            push!(g_di, abs(gene_exp[i,end]- gene_exp[j,end]))
        end
    end
    g_ord = sortperm(g_di)
    return [g_di[g_ord] pcor[g_ord]]
end
function get_gene_null(gene_exp)
    pcor = Float64[]
    g_di = Float64[]
    lg = size(gene_exp)[1]
    for i in 1:lg-1
        for j in i+1:lg
            #try
            push!(pcor, cor(shuffle(gene_exp[i,2:end-1]), gene_exp[j,2:end-1]))
            push!(g_di, abs(gene_exp[i,end]- gene_exp[j,end]))
        end
    end
    g_ord = sortperm(g_di)
    return [g_di[g_ord] pcor[g_ord]]
end
function get_smooth_cor(gene_cor; num_w=1000, q_1=0.25, q_2=0.75)
    w_s = round((gene_cor[:,1][end] - gene_cor[:,1][1] )/ num_w)
    med_dis = Float64[]
    med_cor = Float64[]
    q_cor = []
    for nw in 0:num_w-1
        ini = nw*w_s
        fin = nw*w_s + w_s
        ix_dw = findall(x-> ini<x<fin, gene_cor[:,1])
        if isempty(ix_dw); continue; end
        push!(med_dis, median(gene_cor[ix_dw,1]))
        m_cor = median(gene_cor[ix_dw,2])
        push!(med_cor, m_cor)
        qs = quantile(gene_cor[ix_dw,2],[q_1,q_2])
        push!(q_cor,[m_cor-qs[1] qs[2]-m_cor])
    end
    q_cor =vcat(q_cor...)
    return med_dis, med_cor, q_cor
end
moving_average(vs,n) = [sum(@view vs[i:(i+n-1)])/n for i in 1:(length(vs)-(n-1))]
function get_plot_cor(gene_cor,gene_null)
    dis_cor, val_cor, q_cor = get_smooth_cor(gene_cor)
    dis_null, val_null, q_null = get_smooth_cor(gene_null)
    mv_w = 11
    dis_cor_mv = moving_average(dis_cor,mv_w)
    val_cor_mv = moving_average(val_cor,mv_w)
    q_cor_mv = [moving_average(q_cor[:,1],mv_w) moving_average(q_cor[:,2],mv_w)]
    dis_null_mv = moving_average(dis_null,mv_w)
    val_null_mv = moving_average(val_null,mv_w)
    q_null_mv = [moving_average(q_null[:,1],mv_w) moving_average(q_null[:,2],mv_w)]
    pt = plot(leg=false, #title="$gtitle",
        size=(1100,1000)
    )
    scatter!(pt,gene_cor[:,1]/10^8,gene_cor[:,2],
        alpha=0.08,
        ylims=(-0.3,1.0),
        color="#36983B"
    )
    scatter!(pt,gene_null[:,1]/10^8,gene_null[:,2],
        alpha=0.08,
        color="#EA5108"
    )
    plot!(pt,dis_cor_mv/10^8, val_cor_mv,ribbon=q_cor_mv, color="#36983B",linewidth=2)
    plot!(pt,dis_null_mv/10^8, val_null_mv,ribbon=q_null_mv, color="#EA5108",linewidth=2)
    return pt
end
function get_genes_ch(gene_tags, genes_data)
    g_ch = []
    g_st = []
    ix_x = []
    genes= []
    for g in gene_tags
        ix = findfirst(x-> x==g,genes_data[:,3])
        if isnothing(ix)
            #println("gene $g not found")
            continue
        end
        push!(ix_x, ix)
        push!(genes, g)
        push!(g_ch, genes_data[ix,5])
        push!(g_st, genes_data[ix,4])
    end
    return [genes g_st g_ch]
end

max_stars_repo_name: josemaz/gene-matrices

max_issues_repo_name: josemaz/gene-matrices

max_forks_repo_name: josemaz/gene-matrices