text:
## 1-dimensional ranges ##

typealias Dims (Int...)

abstract Range{T} <: AbstractArray{T,1}

## ordinal ranges

abstract OrdinalRange{T,S} <: Range{T}

immutable StepRange{T,S} <: OrdinalRange{T,S}
    start::T
    step::S
    stop::T

    function StepRange(start::T, step::S, stop::T)
        if T<:FloatingPoint || S<:FloatingPoint
            error("StepRange should not be used with floating point")
        end
        z = zero(S)
        step == z && error("step cannot be zero")
        step != step && error("step cannot be NaN")

        if stop == start
            last = stop
        else
            if (step > z) != (stop > start)
                # empty range has a special representation where stop = start-1
                # this is needed to avoid the wrap-around that can happen computing
                # start - step, which leads to a range that looks very large instead
                # of empty.
                if step > z
                    last = start - one(step)
                else
                    last = start + one(step)
                end
            else
                diff = stop - start
                if T<:Signed && (diff > zero(diff)) != (stop > start)
                    # handle overflowed subtraction with unsigned rem
                    if diff > zero(diff)
                        remain = -oftype(T, unsigned(-diff) % step)
                    else
                        remain = oftype(T, unsigned(diff) % step)
                    end
                else
                    remain = diff % step
                end
                last = stop - remain
            end
        end

        new(start, step, last)
    end
end

StepRange{T,S}(start::T, step::S, stop::T) = StepRange{T,S}(start, step, stop)

immutable UnitRange{T<:Real} <: OrdinalRange{T,Int}
    start::T
    stop::T

    UnitRange(start, stop) = new(start, ifelse(stop >= start, stop, start-1))
end
UnitRange{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop)

colon(a::Real, b::Real) = colon(promote(a,b)...)

colon{T<:Real}(start::T, stop::T) = UnitRange{T}(start, stop)

range(a::Real, len::Integer) = UnitRange{typeof(a)}(a, a+len-1)

colon{T}(start::T, stop::T) = StepRange(start, one(stop-start), stop)

range{T}(a::T, len::Integer) =
    StepRange{T, typeof(a-a)}(a, one(a-a), a+oftype(a-a,(len-1)))

# first promote start and stop, leaving step alone
# this is for non-numeric ranges where step can be quite different
colon{A<:Real,C<:Real}(a::A, b, c::C) = colon(convert(promote_type(A,C),a), b, convert(promote_type(A,C),c))

colon{T<:Real}(start::T, step, stop::T) = StepRange(start, step, stop)

colon{T}(start::T, step, stop::T) = StepRange(start, step, stop)

range{T,S}(a::T, step::S, len::Integer) = StepRange{T,S}(a, step, a+step*(len-1))

## floating point ranges

immutable FloatRange{T<:FloatingPoint} <: Range{T}
    start::T
    step::T
    len::T
    divisor::T
end
FloatRange(a::FloatingPoint, s::FloatingPoint, l::Real, d::FloatingPoint) =
    FloatRange{promote_type(typeof(a),typeof(s),typeof(d))}(a,s,l,d)

# float rationalization helper
function rat(x)
    y = x
    a = d = 1
    b = c = 0
    m = maxintfloat(Float32)
    while abs(y) <= m
        f = itrunc(y)
        y -= f
        a, c = f*a + c, a
        b, d = f*b + d, b
        max(abs(a),abs(b)) <= convert(Int,m) || return c, d
        oftype(x,a)/oftype(x,b) == x && break
        y = inv(y)
    end
    return a, b
end

function colon{T<:FloatingPoint}(start::T, step::T, stop::T)
    step == 0                    && error("range step cannot be zero")
    start == stop                && return FloatRange{T}(start,step,1,1)
    (0 < step) != (start < stop) && return FloatRange{T}(start,step,0,1)

    # float range "lifting"
    r = (stop-start)/step
    n = round(r)
    lo = prevfloat((prevfloat(stop)-nextfloat(start))/n)
    hi = nextfloat((nextfloat(stop)-prevfloat(start))/n)
    if lo <= step <= hi
        a0, b = rat(start)
        a = convert(T,a0)
        if a/convert(T,b) == start
            c0, d = rat(step)
            c = convert(T,c0)
            if c/convert(T,d) == step
                e = lcm(b,d)
                a *= div(e,b)
                c *= div(e,d)
                eT = convert(T,e)
                if (a+n*c)/eT == stop
                    return FloatRange{T}(a, c, n+1, eT)
                end
            end
        end
    end
    FloatRange{T}(start, step, floor(r)+1, one(step))
end

colon{T<:FloatingPoint}(a::T, b::T) = colon(a, one(a), b)

colon{T<:Real}(a::T, b::FloatingPoint, c::T) = colon(promote(a,b,c)...)
colon{T<:FloatingPoint}(a::T, b::FloatingPoint, c::T) = colon(promote(a,b,c)...)
colon{T<:FloatingPoint}(a::T, b::Real, c::T) = colon(promote(a,b,c)...)

range(a::FloatingPoint, len::Integer) = FloatRange(a,one(a),len,one(a))
range(a::FloatingPoint, st::FloatingPoint, len::Integer) = FloatRange(a,st,len,one(a))
range(a::Real, st::FloatingPoint, len::Integer) = FloatRange(float(a), st, len, one(st))
range(a::FloatingPoint, st::Real, len::Integer) = FloatRange(a, float(st), len, one(a))

## interface implementations

similar(r::Range, T::Type, dims::Dims) = Array(T, dims)

size(r::Range) = (length(r),)

isempty(r::StepRange) =
    (r.start != r.stop) && ((r.step > zero(r.step)) != (r.stop > r.start))
isempty(r::UnitRange) = r.start > r.stop
isempty(r::FloatRange) = length(r)==0

step(r::StepRange) = r.step
step(r::UnitRange) = 1
step(r::FloatRange) = r.step/r.divisor

function length(r::StepRange)
    n = integer(div(r.stop+r.step - r.start, r.step))
    isempty(r) ? zero(n) : n
end
length(r::UnitRange) = integer(r.stop - r.start + 1)
length(r::FloatRange) = integer(r.len)

function length{T<:Union(Int,Uint,Int64,Uint64)}(r::StepRange{T})
    isempty(r) && return zero(T)
    if r.step > 1
        return checked_add(oftype(T, div(unsigned(r.stop - r.start), r.step)), one(T))
    elseif r.step < -1
        return checked_add(oftype(T, div(unsigned(r.start - r.stop), -r.step)), one(T))
    else
        checked_add(div(checked_sub(r.stop, r.start), r.step), one(T))
    end
end

length{T<:Union(Int,Uint,Int64,Uint64)}(r::UnitRange{T}) =
    checked_add(checked_sub(r.stop, r.start), one(T))

first{T}(r::OrdinalRange{T}) = oftype(T, r.start)
first(r::FloatRange) = r.start/r.divisor

last{T}(r::StepRange{T}) = r.stop
last(r::UnitRange) = r.stop
last{T}(r::FloatRange{T}) = oftype(T, (r.start + (r.len-1)*r.step)/r.divisor)

minimum(r::UnitRange) = isempty(r) ? error("range must be non-empty") : first(r)
maximum(r::UnitRange) = isempty(r) ? error("range must be non-empty") : last(r)
minimum(r::Range)  = isempty(r) ? error("range must be non-empty") : min(first(r), last(r))
maximum(r::Range)  = isempty(r) ? error("range must be non-empty") : max(first(r), last(r))

ctranspose(r::Range) = [x for _=1, x=r]
transpose(r::Range) = r'

# Ranges are immutable
copy(r::Range) = r

## iteration

start(r::FloatRange) = 0
next{T}(r::FloatRange{T}, i) = (oftype(T, (r.start + i*r.step)/r.divisor), i+1)
done(r::FloatRange, i) = (length(r) <= i)

# NOTE: For ordinal ranges, we assume start+step might be from a
# lifted domain (e.g. Int8+Int8 => Int); use that for iterating.
start(r::StepRange) = convert(typeof(r.start+r.step), r.start)
next{T}(r::StepRange{T}, i) = (oftype(T,i), i+r.step)
done{T,S}(r::StepRange{T,S}, i) = (i!=r.stop) & ((r.step>zero(S))==(i>r.stop))

start(r::UnitRange) = oftype(r.start+1, r.start)
next{T}(r::UnitRange{T}, i) = (oftype(T,i), i+1)
done(r::UnitRange, i) = i==oftype(i,r.stop)+1

## indexing

getindex(r::Range, i::Real) = getindex(r, to_index(i))

function getindex{T}(r::Range{T}, i::Integer)
    1 <= i <= length(r) || error(BoundsError)
    oftype(T, first(r) + (i-1)*step(r))
end
function getindex{T}(r::FloatRange{T}, i::Integer)
    1 <= i <= length(r) || error(BoundsError)
    oftype(T, (r.start + (i-1)*r.step)/r.divisor)
end

function getindex(r::UnitRange, s::UnitRange{Int})
    sl = length(s)
    if sl > 0
        if !(1 <= last(s) <= length(r))
            throw(BoundsError())
        end
        st = r[s.start]
    else
        st = oftype(r.start, r.start + s.start-1)
    end
    range(st, sl)
end

function getindex(r::StepRange, s::Range{Int})
    sl = length(s)
    if sl > 0
        if !(1 <= last(s) <= length(r))
            throw(BoundsError())
        end
        st = r[first(s)]
    else
        st = oftype(r.start, r.start + (first(s)-1)*step(r))
    end
    range(st, step(r)*step(s), sl)
end

getindex(r::FloatRange, s::UnitRange) = r[first(s)]:step(r):r[last(s)]

function show(io::IO, r::Range)
    print(io, repr(first(r)), ':', repr(step(r)), ':', repr(last(r)))
end
show(io::IO, r::UnitRange) = print(io, repr(first(r)), ':', repr(last(r)))

=={T<:Range}(r::T, s::T) = (first(r) == first(s)) & (step(r) == step(s)) & (last(r) == last(s))

function ==(r::Range, s::Range)
    lr = length(r)
    if lr != length(s)
        return false
    end
    u, v = start(r), start(s)
    while !done(r, u)
        x, u = next(r, u)
        y, v = next(s, v)
        if x != y
            return false
        end
    end
    return true
end

intersect{T1<:Integer, T2<:Integer}(r::UnitRange{T1}, s::UnitRange{T2}) = max(r.start,s.start):min(last(r),last(s))

intersect{T<:Integer}(i::Integer, r::UnitRange{T}) =
    i < first(r) ? (first(r):i) :
    i > last(r)  ? (i:last(r))  : (i:i)

intersect{T<:Integer}(r::UnitRange{T}, i::Integer) = intersect(i, r)

function intersect{T1<:Integer, T2<:Integer}(r::UnitRange{T1}, s::StepRange{T2})
    if length(s) == 0
        range(first(r), 0)
    elseif step(s) == 0
        intersect(first(s), r)
    elseif step(s) < 0
        intersect(r, reverse(s))
    else
        sta = first(s)
        ste = step(s)
        sto = last(s)
        lo = first(r)
        hi = last(r)
        i0 = max(sta, lo + mod(sta - lo, ste))
        i1 = min(sto, hi - mod(hi - sta, ste))
        i0:ste:i1
    end
end

function intersect{T1<:Integer, T2<:Integer}(r::StepRange{T1}, s::UnitRange{T2})
    if step(r) < 0
        reverse(intersect(s, reverse(r)))
    else
        intersect(s, r)
    end
end

function intersect{T1<:Integer, T2<:Integer}(r::StepRange{T1}, s::StepRange{T2})
    if length(r) == 0 || length(s) == 0
        return range(first(r), step(r), 0)
    elseif step(s) < 0
        return intersect(r, reverse(s))
    elseif step(r) < 0
        return reverse(intersect(reverse(r), s))
    end

    start1 = first(r)
    step1 = step(r)
    stop1 = last(r)
    start2 = first(s)
    step2 = step(s)
    stop2 = last(s)
    a = lcm(step1, step2)

    # if a == 0
    #     # One or both ranges have step 0.
    #     if step1 == 0 && step2 == 0
    #         return start1 == start2 ? r : Range(start1, 0, 0)
    #     elseif step1 == 0
    #         return start2 <= start1 <= stop2 && rem(start1 - start2, step2) == 0 ? r : Range(start1, 0, 0)
    #     else
    #         return start1 <= start2 <= stop1 && rem(start2 - start1, step1) == 0 ? (start2:step1:start2) : Range(start1, step1, 0)
    #     end
    # end

    g, x, y = gcdx(step1, step2)

    if rem(start1 - start2, g) != 0
        # Unaligned, no overlap possible.
        return range(start1, a, 0)
    end

    z = div(start1 - start2, g)
    b = start1 - x * z * step1
    # Possible points of the intersection of r and s are
    # ..., b-2a, b-a, b, b+a, b+2a, ...
    # Determine where in the sequence to start and stop.
    m = max(start1 + mod(b - start1, a), start2 + mod(b - start2, a))
    n = min(stop1 - mod(stop1 - b, a), stop2 - mod(stop2 - b, a))
    m:a:n
end

function intersect(r::Range, s::Range...)
    i = r
    for t in s
        i = intersect(i, t)
    end
    i
end

# findin (the index of intersection)
function _findin{T1<:Integer, T2<:Integer}(r::Range{T1}, span::UnitRange{T2})
    local ifirst
    local ilast
    fspan = first(span)
    lspan = last(span)
    fr = first(r)
    lr = last(r)
    sr = step(r)
    if sr > 0
        ifirst = fr >= fspan ? 1 : iceil((fspan-fr)/sr)+1
        ilast = lr <= lspan ? length(r) : length(r) - iceil((lr-lspan)/sr)
    elseif sr < 0
        ifirst = fr <= lspan ? 1 : iceil((lspan-fr)/sr)+1
        ilast = lr >= fspan ? length(r) : length(r) - iceil((lr-fspan)/sr)
    else
        ifirst = fr >= fspan ? 1 : length(r)+1
        ilast = fr <= lspan ? length(r) : 0
    end
    ifirst, ilast
end

function findin{T1<:Integer, T2<:Integer}(r::UnitRange{T1}, span::UnitRange{T2})
    ifirst, ilast = _findin(r, span)
    ifirst:ilast
end

function findin{T1<:Integer, T2<:Integer}(r::Range{T1}, span::UnitRange{T2})
    ifirst, ilast = _findin(r, span)
    ifirst:1:ilast
end

## linear operations on ranges ##

-(r::OrdinalRange) = range(-r.start, -step(r), length(r))
-(r::FloatRange)   = FloatRange(-r.start, -r.step, r.len, r.divisor)

+(x::Real, r::UnitRange)  = range(x + r.start, length(r))
+(x::Real, r::Range) = (x+first(r)):step(r):(x+last(r))
#+(x::Real, r::StepRange)  = range(x + r.start, r.step, length(r))
+(x::Real, r::FloatRange) = FloatRange(r.divisor*x + r.start, r.step, r.len, r.divisor)
+(r::Range, x::Real)      = x + r
#+(r::FloatRange, x::Real) = x + r

-(x::Real, r::Range)      = (x-first(r)):-step(r):(x-last(r))
-(x::Real, r::FloatRange) = FloatRange(r.divisor*x - r.start, -r.step, r.len, r.divisor)
-(r::UnitRange, x::Real)  = range(r.start-x, length(r))
-(r::StepRange , x::Real) = range(r.start-x, r.step, length(r))
-(r::FloatRange, x::Real) = FloatRange(r.start - r.divisor*x, r.step, r.len, r.divisor)

.*(x::Real, r::OrdinalRange) = range(x*r.start, x*step(r), length(r))
.*(x::Real, r::FloatRange)   = FloatRange(x*r.start, x*r.step, r.len, r.divisor)
.*(r::Range, x::Real)        = x .* r
.*(r::FloatRange, x::Real)   = x .* r

./(r::OrdinalRange, x::Real) = range(r.start/x, step(r)/x, length(r))
./(r::FloatRange, x::Real)   = FloatRange(r.start/x, r.step/x, r.len, r.divisor)

promote_rule{T1,T2}(::Type{FloatRange{T1}},::Type{FloatRange{T2}}) =
    FloatRange{promote_type(T1,T2)}
convert{T}(::Type{FloatRange{T}}, r::FloatRange) =
    FloatRange{T}(r.start,r.step,r.len,r.divisor)

promote_rule{F,OR<:OrdinalRange}(::Type{FloatRange{F}}, ::Type{OR}) =
    FloatRange{promote_type(F,eltype(OR))}
convert{T}(::Type{FloatRange{T}}, r::OrdinalRange) =
    FloatRange{T}(start(r), step(r), length(r), one(T))

# +/- of ranges is defined in operators.jl (to be able to use @eval etc.)

## non-linear operations on ranges ##

./(x::Number, r::Range) = [ x/y for y=r ]
./(r::Range, y::Number) = [ x/y for x=r ]

.^(x::Number, r::Range) = [ x^y for y=r ]
.^(r::Range, y::Number) = [ x^y for x=r ]

## concatenation ##

function vcat{T}(r::Range{T})
    n = length(r)
    a = Array(T,n)
    i = 1
    for x in r
        @inbounds a[i] = x
        i += 1
    end
    return a
end

convert{T}(::Type{Array{T,1}}, r::Range{T}) = vcat(r)

function vcat{T}(rs::Range{T}...)
    n = sum(length,rs)::Int
    a = Array(T,n)
    i = 1
    for r in rs
        for x in r
            @inbounds a[i] = x
            i += 1
        end
    end
    return a
end

reverse(r::OrdinalRange) = range(last(r), -step(r), length(r))
reverse(r::FloatRange)   = FloatRange(last(r), -r.step, r.len, r.divisor)

## sorting ##

issorted(r::UnitRange) = true
issorted(r::Range) = step(r) >= 0

sort(r::UnitRange) = r
sort!(r::UnitRange) = r

sort{T<:Real}(r::Range{T}) = issorted(r) ? r : reverse(r)

sortperm(r::UnitRange) = 1:length(r)
sortperm{T<:Real}(r::Range{T}) = issorted(r) ? (1:1:length(r)) : (length(r):-1:1)

function sum{T<:Real}(r::Range{T})
    l = length(r)
    # note that a little care is required to avoid overflow in l*(l-1)/2
    return l * first(r) + (iseven(l) ? (step(r) * (l-1)) * (l>>1)
                                     : (step(r) * l) * ((l-1)>>1))
end

function map!(f::Callable, dest, r::Range)
    i = 1
    for ri in r dest[i] = f(ri); i+=1; end
    dest
end

function map_range_to!(f::Callable, first, dest, r::Range, state)
    dest[1] = first
    i = 2
    while !done(r, state)
        ri, state = next(r, state)
        dest[i] = f(ri)
        i += 1
    end
    dest
end

function map(f::Callable, r::Range)
    if isempty(r); return {}; end
    state = start(r)
    (ri, state) = next(r, state)
    first = f(ri)
    map_range_to!(f, first, Array(typeof(first), length(r)), r, state)
end

function in(x, r::Range)
    n = step(r) == 0 ? 1 : iround((x-first(r))/step(r))+1
    n >= 1 && n <= length(r) && r[n] == x
end

in{T<:Integer}(x, r::Range{T}) = isinteger(x) && !isempty(r) && x>=minimum(r) && x<=maximum(r) && (step(r)==0 || mod(int(x)-first(r),step(r))==0)

max_stars_repo_name: aplund/julia

max_issues_repo_name: aplund/julia

max_forks_repo_name: aplund/julia


text:
# test/runtests.jl
using Reproduce, Test, FileIO

include("parse.jl")
include("search.jl")
include("args_iter.jl")

function tests()

    PARSE_TESTS.@testparse
    SEARCH_TESTS.@testsearch
    ARGSITER_TESTS.@tests
end

tests()

ext: jl


text:
export minimumjerk, minimumjerkf, functionform

"""
Retorna os coeficientes do polinômio de trajetória com base no princípio do minimum jerk considerando o tempo inicial sendo 0(zero).
Entradas
x0 -> posição inicial
v0 -> velocidade inicial
a0 -> aceleração inicial
t0 -> tempo inicial
xf -> posição final
vf -> velocidade final
af -> aceleração final
tf -> tempo para chegar ao ponto final
"""
function minimumjerk(x0::T, v0::T, a0::T, t0::T, xf::T, vf::T, af::T, tf::T)  where {T<:AbstractFloat}
    b=[ 1 t0 t0^2 t0^3 t0^4 t0^5;
    0 1 2*t0 3*(t0^2) 4*(t0^3) 5*(t0^4);
    0 0 2 6*t0 12*(t0^2) 20*(t0^3);
    1 tf tf^2 tf^3 tf^4 tf^5;
    0 1 2*tf 3*(tf^2) 4*(tf^3) 5*(tf^4);
    0 0 2 6*tf 12*(tf^2) 20*(tf^3)]
    c = [x0, v0, a0, xf, vf, af]
    a = inv(b)*c
    return a
end

"""
Retorna os coeficientes do polinômio de trajetória com base no princípio do minimum jerk considerando o tempo inicial sendo 0(zero).
Entradas
x0 -> posição inicial
v0 -> velocidade inicial
a0 -> aceleração inicial
xf -> posição final
vf -> velocidade final
af -> aceleração final
T -> tempo para chegar ao ponto final
"""
function minimumjerk(x0::T, v0::T, a0::T, xf::T, vf::T, af::T, tf::T) where T<:AbstractFloat
    t0 = 0
    minimumjerk(x0,v0,a0,t0,xf,vf,af,tf)
end

function minimumjerkf(x0::T, v0::T, a0::T, t0::T, xf::T, vf::T, af::T, tf::T) where T<:AbstractFloat
    functionform(minimumjerk(x0,v0,a0,t0,xf,vf,af,tf))
end

function minimumjerkf(x0::T, v0::T, a0::T, xf::T, vf::T, af::T, tf::T) where T<:AbstractFloat
    t0 = 0.0
    functionform(minimumjerk(x0,v0,a0,t0,xf,vf,af,tf))
end

# Retorna as funções de posição, velocidade, aceleração e jerk a parti dos coeficientes do polinômio do mínimo jerk
# Entradas
# a -> vetor com os coeficientes do polinômio
function functionform(a::Vector)
    @eval out1 = $(Meta.parse("t -> $(a[1]) + $(a[2])*t + $(a[3])*(t^2) + $(a[4])*(t^3) + $(a[5])*(t^4) + $(a[6])*(t^5)"))
    @eval out2 = $(Meta.parse("t -> $(a[2]) + 2*$(a[3])*t + 3*$(a[4])*(t^2) + 4*$(a[5])*(t^3) + 5*$(a[6])*(t^4)"))
    @eval out3 = $(Meta.parse("t ->  2*$(a[3]) + 6*$(a[4])*t + 12*$(a[5])*(t^2) + 20*$(a[6])*(t^3)"))
    @eval out4 = $(Meta.parse("t -> 6*$(a[4]) + 24*$(a[5])*t + 60*$(a[6])*(t^2)"))
    return out1, out2, out3, out4
end

# Retorna as funções de posição, velocidade, aceleração e jerk
# partindo de um vetor de posições finais de cada junta e do tempo final
# essa função considera que as condições iniciais (psição, velocidade e aceleração) são
# nulas em cada junta
function minimumjerk(final_positions::Vector{T},tend::T) where T<:AbstractFloat
    xr = Function[]
    vr = Function[]
    ar = Function[]
    jr = Function[]
    for i in final_positions
        x, v, a, j = minimumjerkf(0.0, 0.0, 0.0,0.0, i, 0.0, 0.0, tend)
        push!(xr,x)
        push!(vr,v)
        push!(ar,a)
        push!(jr,j)
    end
    xr, vr, ar, jr
end


text:
using KLIEPInference
using JLD

scratch_dir = ARGS[1]

p = 0.05:0.05:0.95

for m in [25, 50, 100]
    for sgn in [-1, 0, 1]
        coverage = zeros(length(p), 6)

        nrep = 0
        for rep in 1:1000
            try
                file = jldopen("$(scratch_dir)/res_$(m)_$(sgn)_$(rep).jld", "r")
                boot1 = read(file, "boot1")
                boot2 = read(file, "boot2")
                close(file)

                q = boot_quantile(boot1, p)
                covered1 = maximum(abs.(boot1.θhat)) .<= q

                q = boot_quantile(boot2, p)
                covered2 = maximum(abs.(boot2.θhat)) .<= q

                coverage[:, 1] .+= covered1
                coverage[:, 2] .+= covered2
                coverage[:, 3] .+= covered1 .== covered2

                q, w = boot_quantile_studentized(boot1, p)
                covered1 = maximum(abs.(boot1.θhat ./ w)) .<= q

                q, w = boot_quantile_studentized(boot2, p)
                covered2 = maximum(abs.(boot2.θhat ./ w)) .<= q

                coverage[:, 4] .+= covered1
                coverage[:, 5] .+= covered2
                coverage[:, 6] .+= covered1 .== covered2

                nrep += 1
            catch
                @warn "could not open $(m) $(sgn) $(rep)"
            end
        end

        coverage = broadcast(/, coverage, nrep)

        @save "$(scratch_dir)/coverage_$(m)_$(sgn).jld" coverage
    end
end


text:
include("./relaxedSequencingNode.jl")
include("./relaxedStateSequencingNode.jl")
include("./relaxedDDSubRoutines.jl")
include("./relaxedDD.jl")
include("./roughRelaxedBounding.jl")


text:
module HRP

using Statistics, LinearAlgebra, Clustering

export hrp, ivp

dist(x, y) = norm(collect(skipmissing(x .- y)), 2)

function covar(x, y)
    id = findall(z -> !ismissing(z), x .- y)
    return cov(x[id], y[id])
end

function ret2dist(r)
    n = size(r, 2)
    Σ = fill(0.0, n, n)
    D = fill(0.0, n, n)
    for i = 1:n
        for j = i:n
            D[i, j] = dist(r[:, i], r[:, j])
            D[j, i] = D[i, j]
            Σ[i, j] = covar(r[:, i], r[:, j])
            Σ[j, i] = Σ[i, j]
        end
    end
    return Σ, D
end

function ivp(C)
    iv = 1.0 ./ diag(C)
    return iv ./ sum(iv)
end

function clustervar(C, cluster)
    V = C[cluster, cluster]
    u = ivp(V)
    return dot(u' * V, u)
end

function bisect(items)
    tmp = Vector{Int}[]
    for item in items
        n1 = length(item)
        if n1 > 1
            n2 = floor(Int, n1 / 2)
            push!(tmp, item[1:n2])
            push!(tmp, item[n2+1:end])
        end
    end
    return tmp
end

function update!(w, items, Σ)
    for k in 1:2:length(items)
        c1, c2 = items[k], items[k+1]
        v1, v2 = clustervar(Σ, c1), clustervar(Σ, c2)
        α = v1 / (v1 + v2)
        w[c1] .= w[c1] .* α
        w[c2] .= w[c2] .* (1.0 - α)
    end
end

function hrp(x; input_type="ret")
    Σ = fill(0.0, size(x, 2), size(x, 2))
    D = similar(Σ)
    if input_type == "ret"
        Σ, D = ret2dist(x)
    elseif input_type == "cov"
        Σ, D = x, cov2dist(x)
    end
    h = hclust(D).order
    n = length(h)
    w = fill(1.0, n)
    items = [h]
    while !isempty(items)
        items = bisect(items)
        update!(w, items, Σ)
    end
    return ivp(Σ), w
end

#= function hrp(Σ)
    D = sqrt.(2.0*(1.0.-cov2corr(Σ)))
    id = 1:size(D,2)
    w = fill(1.0, size(D,2))
    while length(id) > 1
        id, k = split(D, id)
        v1 = clustervar(Σ, id)
        v2 = Σ[k,k]
        α = ivp([v1,v2])[1]
        w[id] .= w[id] * α
        w[k] = w[k] * (1. - α)
    end
    return w
end =#

function split(D::Matrix{Float64}, id::Union{UnitRange{Int64},Vector{Int64}})
    if length(id) == 1
        return nothing, id
    end
    tmp = []
    for i in id
        push!(tmp, minimum([D[i, k] for k in setdiff(id, i)]))
    end
    k = id[argmax(tmp)]
    id = setdiff(id, k)
    return id, k
end

function hcluster_order(D)
    n = size(D, 2)
    id = 1:n
    out = []
    for i = 1:n
        id, k = split(D, id)
        push!(out, k)
    end
    return reverse(collect(Iterators.flatten(out)))
end

function cov2corr(C)
    σinv = diagm(1.0 ./ sqrt.(diag(C)))
    return σinv * C * σinv
end

function cov2dist(C)
    R = cov2corr(C)
    for i ∈ 1:size(R, 2)
        R[i, i] = 1.0
    end
    D = sqrt.(0.5 * (1.0 .- R))
    for i ∈ 1:size(D, 2)
        for j ∈ i+1:size(D, 2)
            D[i, j] = D[j, i]
        end
    end
    return D
end

end

max_stars_repo_name: banachtech/hrp

max_issues_repo_name: banachtech/hrp

max_forks_repo_name: banachtech/hrp


text:
@system VaporPressure begin
    # Campbell and Norman (1998), p 41 Saturation vapor pressure in kPa
    a => 0.611 ~ preserve(u"kPa", parameter)
    b => 17.502 ~ preserve(parameter)
    c => 240.97 ~ preserve(parameter) # °C

    es(a, b, c; T(u"°C")): saturation => (t = Cropbox.deunitfy(T); a*exp((b*t)/(c+t))) ~ call(u"kPa")
    ea(es; T(u"°C"), RH(u"percent")): ambient => es(T) * RH ~ call(u"kPa")
    D(es; T(u"°C"), RH(u"percent")): deficit => es(T) * (1 - RH) ~ call(u"kPa")
    RH(es; T(u"°C"), VPD(u"kPa")): relative_humidity => 1 - VPD / es(T) ~ call(u"NoUnits")

    # slope of the sat vapor pressure curve: first order derivative of Es with respect to T
    Δ(es, b, c; T(u"°C")): saturation_slope_delta => (e = es(T); t = Cropbox.deunitfy(T); e*(b*c)/(c+t)^2 / u"K") ~ call(u"kPa/K")
    s(Δ; T(u"°C"), P(u"kPa")): saturation_slope => Δ(T) / P ~ call(u"K^-1")
end


text:
using ITensors
using Test
import Random: seed!

seed!(12345)

@testset "ITensor contraction sequence optimization" begin
  d = 100
  i = Index(d, "i")
  A = randomITensor(i', dag(i))

  @test !ITensors.using_contraction_sequence_optimization()

  A2 = A' * A
  @test hassameinds(A2, (i'', i))

  ITensors.enable_contraction_sequence_optimization()
  @test ITensors.using_contraction_sequence_optimization()

  @test A' * A ≈ A2

  ITensors.disable_contraction_sequence_optimization()
  @test !ITensors.using_contraction_sequence_optimization()

  A3 = A'' * A' * A
  @test hassameinds(A3, (i''', i))
  @test contract([A'', A', A]) ≈ A3
  @test contract([A'', A', A]; sequence = "automatic") ≈ A3
  @test contract([A'', A', A]; sequence = "left_associative") ≈ A3
  @test contract([A'', A', A]; sequence = "right_associative") ≈ A3
  @test contract([A'', A', A]; sequence = [[1, 2], 3]) ≈ A3
  @test contract([A'', A', A]; sequence = [[2, 3], 1]) ≈ A3
  # A bad sequence
  @test contract([A'', A', A]; sequence = [[1, 3], 2]) ≈ A3

  ITensors.enable_contraction_sequence_optimization()
  @test ITensors.using_contraction_sequence_optimization()

  @test A'' * A' * A ≈ A3
  @test A * A'' * A' ≈ A3
  @test contract([A'', A', A]) ≈ A3
  @test contract([A, A'', A']) ≈ A3
  @test contract([A'', A', A]; sequence = "automatic") ≈ A3
  @test contract([A'', A', A]; sequence = "left_associative") ≈ A3
  @test contract([A'', A', A]; sequence = "right_associative") ≈ A3
  @test contract([A'', A', A]; sequence = [[1, 2], 3]) ≈ A3
  @test contract([A'', A', A]; sequence = [[2, 3], 1]) ≈ A3
  @test contract([A'', A', A]; sequence = [[1, 3], 2]) ≈ A3

  ITensors.disable_contraction_sequence_optimization()
  @test !ITensors.using_contraction_sequence_optimization()

  # This is not the only sequence
  @test ITensors.optimized_contraction_sequence([A, A'', A']) == Any[1, Any[2, 3]]

  time_without_opt = @elapsed A * A'' * A'

  ITensors.enable_contraction_sequence_optimization()
  @test ITensors.using_contraction_sequence_optimization()

  time_with_opt = @elapsed A * A'' * A'

  @test time_with_opt < time_without_opt

  ITensors.disable_contraction_sequence_optimization()
  @test !ITensors.using_contraction_sequence_optimization()

  A4 = A''' * A'' * A' * A
  @test hassameinds(A4, (i'''', i))
  @test contract([A''', A'', A', A]; sequence = [[[1, 2], 3], 4]) ≈ A4

  ITensors.enable_contraction_sequence_optimization()
  @test ITensors.using_contraction_sequence_optimization()

  @test A'' * A * A''' * A' ≈ A4

  ITensors.disable_contraction_sequence_optimization()
  @test !ITensors.using_contraction_sequence_optimization()

  # This is not the only sequence
  @test ITensors.optimized_contraction_sequence([A, A'', A', A''']) == Any[Any[1, 3], Any[2, 4]]

  time_without_opt = @elapsed A * A'' * A' * A'''

  ITensors.enable_contraction_sequence_optimization()
  @test ITensors.using_contraction_sequence_optimization()

  time_with_opt = @elapsed A * A'' * A' * A'''

  @test time_with_opt < time_without_opt

  ITensors.disable_contraction_sequence_optimization()
  @test !ITensors.using_contraction_sequence_optimization()
end


text:
using CairoMakie, MLBase, Random
Random.seed!(123)
# fake data
gt = rand(1:10,1000)
pred = rand(1:10,1000)
data = confusmat(10, gt, pred)
ticks = string.(sort(unique(gt)))
# plot
fig = Figure(resolution = (600, 400))
ax = Axis(fig[1,1]; aspect = 1, ylabel = "Ground truth",
    xlabel = "Prediction")
hmap = heatmap!(ax, data, colormap = :magma)
for i in 1:10, j in 1:10
    txtcolor = data[i,j] < 10 ? :white : :black
    text!(ax, "$(round(data[i,j], digits = 2))", position = (i,j),
        color = txtcolor, align = (:center, :center), textsize = 12)
end
cbar = Colorbar(fig[1,2], hmap, label = "Counts", width = 25,
    halign = :left, ticksize = 25, tickalign = 1,
    height = Relative(3.55/4))
ax.xticks = (1:10, ticks)
ax.yticks = (1:10, ticks)
save("confusionMatrix.png", fig, px_per_unit = 2)
fig


text:
@testset "dlanczos" begin
println("Testing `dlanczos`")

# create A = -Δ , B = diag(V)   with V positive
d = 30
x = range(-1,1,length=d)
A = d^2 * SymTridiagonal(2*ones(d), -ones(d-1))
B = Diagonal(exp.(- x.^2) .+ 0.1)

# construct an index-1 system
μ = 30.0
H = A - μ * B
σ = sort(eigvals(Matrix(H)))
@assert σ[1] < 0 && σ[2] > 0

b = ones(d) + 0.1 * sin.(x)

# define a linear function
fl = x -> H * x

# then try the same with a quadratic nonlinearity mixed in
#   (hessian remains the same!)
q = x -> [ [x[i]*x[i+1] for i = 1:d-1]; x[end] * x[1] ]
fq = x -> H * x + q(x)

# testing unpreconditioned lanczos
for (f, V0, msg) in [ (fl, reshape(b, d, 1), "Linear - Basic Lanczos"),
                      (fl, [b rand(d)], "Linear - Block Lanczos"),
                      (fq, reshape(b, d, 1), "Nonlinear - Basic Lanczos"),
                      (fq, [b rand(d)], "Nonlinear - Block Lanczos")
                     ]
   @info("Testing: $msg")
   x, G, numf, success = dlanczos(zeros(d), f, zeros(d), b, 1e-6, d; V0 = V0, debug = false)
   λ, v = G[1]
   print_tf(@test success)
   print_tf(@test norm(H * x - b) < 1e-6)
   print_tf(@test norm(H * v - λ * v, Inf) < 1e-6)
   print_tf(@test abs(λ - σ[1]) < 1e-7)
   println()
end

# next we add some preconditioning
P = 0.9 * sparse(A) + μ/2 * I
λP = minimum(eigvals(Matrix(H), Matrix(P)))
Random.seed!(12345)
vrand = rand(d)

for (V0, msg, numfo) in [ (reshape(P\b, d, 1), "Preconditioned Lanczos", 9),
                           ([P\b P\vrand], "Preconditioned Block-Lanczos", 14) ]
   @info("Testing: $msg")
   x, G, numf, success = dlanczos(zeros(d), fq, zeros(d), b, 1e-6, d;
                                P = P, debug = false, V0=V0)
   print_tf(@test numf <= numfo)
   print_tf(@test success)
   λ, v = G[1]
   print_tf(@test norm(H * x - b) < 1e-6)
   print_tf(@test norm(H * v - λ * P * v, Inf) < 1e-6)
   print_tf(@test abs(λ - λP) < 1e-6)
   if numf != numfo
      warn("numf was $numfo in original tests; now it is $numf")
   end
   println()
end

end  # @testset "blocklanczos"


text:
module TestNoEffectTempSize
    using BioEnergeticFoodWebs
    using Test

    # WET MASS (default)
    EC = [0 0 1 ; 0 0 1 ; 0 0 0]
    LFC = [0 1 0 ; 0 0 1 ; 0 0 0]
    bsize_ratio = 10.0 #consumers are 10 times larger than their resource
    p_ec = model_parameters(EC, Z = bsize_ratio)
    p_lfc = model_parameters(LFC, Z = bsize_ratio)
    #body sizes are CR ratio ^ (trophic rank - 1)
    expected_mass_ec = [10.0, 10.0, 1.0]
    @test p_ec[:bodymass] == expected_mass_ec
    expected_mass_lfc = [100.0, 10.0, 1.0]
    @test p_lfc[:bodymass] == expected_mass_lfc

    #if body masses are provided, it overwrites Z
    bm = [12.3, 2.3, 1.0]
    p_bm = model_parameters(EC, Z = bsize_ratio, bodymass = bm)
    @test p_bm[:bodymass] == bm

    # DRY MASS
    dm = [6.3, 5.2, 1.2]
    p_dm = model_parameters(EC, TSR_type = :no_response, dry_mass_293 = dm)
    expected_mass_dm = dm .* 6.5
    @test p_dm[:bodymass] == expected_mass_dm
end

module TestEffectTempSize
    using BioEnergeticFoodWebs
    using Statistics
    using Test

    EC = [0 0 1 ; 0 0 1 ; 0 0 0]
    dm = [1.0, 0.8, 0.2]
    wm = dm .* 6.5
    temp = 295.0
    temp_c = temp-273.15
    temp2 = 285.0
    temp2_c = temp2-273.15

    # MEAN AQUATIC
    # dry mass
    p_aqua_1 = model_parameters(EC, TSR_type = :mean_aquatic, dry_mass_293 = dm, T = temp)
    pcm_aqua_1 = -3.90 .- 0.53 .* log10.(dm)
    TSr_aqua_1 = log.(pcm_aqua_1 ./ 100 .+ 1)
    expected_bm_aqua = (dm .* 6.5) .* exp.(TSr_aqua_1 .* (temp_c-20))
    @test p_aqua_1[:bodymass] == expected_bm_aqua
    p_aqua_2 = model_parameters(EC, TSR_type = :mean_aquatic, dry_mass_293 = dm, T = temp2)
    @test p_aqua_2[:bodymass] > p_aqua_1[:bodymass]
    # wet mass
    p_aquawm = model_parameters(EC, TSR_type = :mean_aquatic, bodymass = wm, T = temp)
    pcm_aquawm = -3.90 .- 0.53 .* log10.(dm)
    TSr_aquawm = log.(pcm_aquawm ./ 100 .+ 1)
    expected_bm_aquawm = wm .* exp.(TSr_aquawm .* (temp_c-20))
    @test p_aquawm[:bodymass] == expected_bm_aqua == expected_bm_aquawm
    # Z
    p_aqua_z = model_parameters(EC, TSR_type = :mean_aquatic, Z = 10.0, T = temp)
    pcm_aqua_z = -3.90 .- 0.53 .* log10.([10.0, 10.0, 1.0] ./ 6.5)
    TSr_aqua_z = log.(pcm_aqua_z ./ 100 .+ 1)
    expected_bm_aqua_z = [10.0, 10.0, 1.0] .* exp.(TSr_aqua_z .* (temp_c-20))
    @test p_aqua_z[:bodymass] == expected_bm_aqua_z

    # MEAN TERRESTRIAL
    # dry mass
    p_terr_1 = model_parameters(EC, TSR_type = :mean_terrestrial, dry_mass_293 = dm, T = temp)
    pcm_terr_1 = -1.72 .+ 0.54 .* log10.(dm)
    TSr_terr_1 = log.(pcm_terr_1 ./ 100 .+ 1)
    expected_bm_terr = (dm .* 6.5) .* exp.(TSr_terr_1 .* (temp_c-20))
    @test p_terr_1[:bodymass] == expected_bm_terr
    p_terr_2 = model_parameters(EC, TSR_type = :mean_terrestrial, dry_mass_293 = dm, T = temp2)
    @test p_terr_2[:bodymass] > p_terr_1[:bodymass]
    # wet mass
    p_terrwm = model_parameters(EC, TSR_type = :mean_terrestrial, bodymass = wm, T = temp)
    pcm_terrwm = -1.72 .+ 0.54 .* log10.(dm)
    TSr_terrwm = log.(pcm_terrwm ./ 100 .+ 1)
    expected_bm_terrwm = wm .* exp.(TSr_terrwm .* (temp_c-20))
    @test p_terrwm[:bodymass] == expected_bm_terr == expected_bm_terrwm
    # Z
    p_terr_z = model_parameters(EC, TSR_type = :mean_terrestrial, Z = 10.0, T = temp)
    pcm_terr_z = -1.72 .+ 0.54 .* log10.([10.0, 10.0, 1.0] ./ 6.5)
    TSr_terr_z = log.(pcm_terr_z ./ 100 .+ 1)
    expected_bm_terr_z = [10.0, 10.0, 1.0] .* exp.(TSr_terr_z .* (temp_c-20))
    @test p_terr_z[:bodymass] == expected_bm_terr_z

    # MAXIMUM
    # dry mass
    p_max_1 = model_parameters(EC, TSR_type = :maximum, dry_mass_293 = dm, T = temp)
    pcm_max_1 = -8
    TSr_max_1 = log.(pcm_max_1 ./ 100 .+ 1)
    expected_bm_max = (dm .* 6.5) .* exp.(TSr_max_1 .* (temp_c-20))
    @test p_max_1[:bodymass] == expected_bm_max
    p_max_2 = model_parameters(EC, TSR_type = :maximum, dry_mass_293 = dm, T = temp2)
    @test p_max_2[:bodymass] > p_max_1[:bodymass]
    # wet mass
    p_maxwm = model_parameters(EC, TSR_type = :maximum, bodymass = wm, T = temp)
    pcm_maxwm = -8
    TSr_maxwm = log.(pcm_maxwm ./ 100 .+ 1)
    expected_bm_maxwm = wm .* exp.(TSr_maxwm .* (temp_c-20))
    @test p_maxwm[:bodymass] == expected_bm_max == expected_bm_maxwm
    # Z
    p_max_z = model_parameters(EC, TSR_type = :maximum, Z = 10.0, T = temp)
    pcm_max_z = -8
    TSr_max_z = log.(pcm_max_z ./ 100 .+ 1)
    expected_bm_max_z = [10.0, 10.0, 1.0] .* exp.(TSr_max_z .* (temp_c-20))
    @test p_max_z[:bodymass] == expected_bm_max_z

    # REVERSE
    # dry mass
    p_rev_1 = model_parameters(EC, TSR_type = :reverse, dry_mass_293 = dm, T = temp)
    pcm_rev_1 = 4
    TSr_rev_1 = log.(pcm_rev_1 ./ 100 .+ 1)
    expected_bm_rev = (dm .* 6.5) .* exp.(TSr_rev_1 .* (temp_c-20))
    @test p_rev_1[:bodymass] == expected_bm_rev
    p_rev_2 = model_parameters(EC, TSR_type = :reverse, dry_mass_293 = dm, T = temp2)
    @test p_rev_2[:bodymass] < p_rev_1[:bodymass]
    # wet mass
    p_revwm = model_parameters(EC, TSR_type = :reverse, bodymass = wm, T = temp)
    pcm_revwm = 4
    TSr_revwm = log.(pcm_revwm ./ 100 .+ 1)
    expected_bm_revwm = wm .* exp.(TSr_revwm .* (temp_c-20))
    @test p_revwm[:bodymass] == expected_bm_rev == expected_bm_revwm
    # Z
    p_rev_z = model_parameters(EC, TSR_type = :reverse, Z = 10.0, T = temp)
    pcm_rev_z = 4
    TSr_rev_z = log.(pcm_rev_z ./ 100 .+ 1)
    expected_bm_rev_z = [10.0, 10.0, 1.0] .* exp.(TSr_rev_z .* (temp_c-20))
    @test p_rev_z[:bodymass] == expected_bm_rev_z

end

max_stars_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_issues_repo_name: SimmonsBI/BioEnergeticFoodWebs

max_forks_repo_name: SimmonsBI/BioEnergeticFoodWebs


text:
@testset "Static graphs" begin
    @testset "DiGraphs" begin
        g = @inferred(complete_digraph(5))
        @test nv(g) == 5 && ne(g) == 20
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = complete_digraph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = complete_digraph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = complete_digraph(2)
        @test nv(g) == 2 && ne(g) == 2
        g = @inferred complete_digraph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127 * (127 - 1)
        @test eltype(g) == Int8
    end

    @testset "Graphs" begin
        g = @inferred(complete_graph(5))
        @test nv(g) == 5 && ne(g) == 10
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = complete_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = complete_graph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = complete_graph(2)
        @test nv(g) == 2 && ne(g) == 1
        g = @inferred complete_graph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127 * (127 - 1) ÷ 2
        @test eltype(g) == Int8
    end

    @testset "Bipartite Graphs" begin
        g = @inferred(complete_bipartite_graph(5, 8))
        @test nv(g) == 13 && ne(g) == 40
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = complete_bipartite_graph(0,0)
        @test nv(g) == 0 && ne(g) == 0
        g = complete_bipartite_graph(5,0)
        @test nv(g) == 5 && ne(g) == 0
        g = complete_bipartite_graph(0,5)
        @test nv(g) == 5 && ne(g) == 0
        g = @inferred complete_bipartite_graph(Int8(100), Int8(27))
        @test nv(g) == 127 && ne(g) == 100 * 27
        @test eltype(g) == Int8
        g = complete_bipartite_graph(Int8(127), Int8(0))
        @test nv(g) == 127 && ne(g) == 0
        g = complete_bipartite_graph(Int8(0), Int8(127))
        @test nv(g) == 127 && ne(g) == 0
    end

    function iscompletemultipartite(g, partitions)
      sum(partitions) != nv(g) && return false
      n = nv(g)

      edges = 0
      for p in partitions
        edges += p*(Int(n)-p)
      end
      edges = div(edges, 2)

      edges != ne(g) && return false

      cur = 1
      for p in partitions
        currange = cur:(cur+p-1)
        lowerrange = 1:(cur-1)
        upperrange = (cur+p):n
        for u in currange
          for v in currange # check that no vertices are connected to vertices in the same partition
            has_edge(g, u, v) && return false
          end

          for v in lowerrange # check all lower partition vertices
            !has_edge(g, u, v) && return false
          end

          for v in upperrange # check all higher partition vertices
            !has_edge(g, u, v) && return false
          end
        end
        cur += p
      end
      return true
    end

    @testset "Multipartite Graphs" begin
        g = @inferred(complete_multipartite_graph([5, 8, 3]))
        @test nv(g) == 16 && ne(g) == 79
        @test isvalid_simplegraph(g)
        @test iscompletemultipartite(g, [5, 8, 3])
        # tests for extreme values
        g = complete_multipartite_graph([0, 0, 0])
        @test nv(g) == 0 && ne(g) == 0
        g = complete_multipartite_graph(Int[])
        @test nv(g) == 0 && ne(g) == 0
        g = complete_multipartite_graph([5, 0, 3])
        @test nv(g) == 8 && ne(g) == 15
        @test iscompletemultipartite(g, [5, 0, 3])
        g = complete_multipartite_graph(Int8[100, 25, 2])
        @test nv(g) == 127 && ne(g) == 2750
        @test eltype(g) == Int8
        @test iscompletemultipartite(g, [100, 25, 2])
        g = complete_multipartite_graph([0, 10, 0, 2, 0, 3, 0, 0, 6, 0])
        @test nv(g) == 21 && ne(g) == 146
        @test iscompletemultipartite(g, [0, 10, 0, 2, 0, 3, 0, 0, 6, 0])
        g = complete_multipartite_graph([5, 3, -1])
        @test nv(g) == 0 && ne(g) == 0
        g = complete_multipartite_graph([-5, 0, -11])
        @test nv(g) == 0 && ne(g) == 0
        g = complete_multipartite_graph([5, 3])
        h = complete_multipartite_graph([5, 0, 3])
        j = complete_multipartite_graph([0, 5, 0, 3, 0, 0])
        @test g == h == j
        @test_throws InexactError complete_multipartite_graph(Int8[100, 20, 20])
    end

    function retrievepartitions(n, r)
      partitions = partitions = Vector{Int}(undef, r)
      c = cld(n,r)
      f = fld(n,r)
      for i in 1:(n%r)
        partitions[i] = c
      end
      for i in ((n%r)+1):r
        partitions[i] = f
      end
      return partitions
    end

    @testset "Turan Graphs" begin
        g = @inferred(turan_graph(13, 4))
        @test nv(g) == 13 && ne(g) == 63
        @test isvalid_simplegraph(g)
        @test iscompletemultipartite(g, retrievepartitions(13, 4))
        # testing smaller int types
        g = turan_graph(Int8(11), Int8(6))
        @test nv(g) == 11 && ne(g) == 50
        @test iscompletemultipartite(g, retrievepartitions(11, 6))
        g = turan_graph(Int16(35), Int16(17))
        @test nv(g) == 35 && ne(g) == 576
        @test iscompletemultipartite(g, retrievepartitions(35, 17))
        # tests for extreme values
        g = turan_graph(15,15)
        @test nv(g) == 15 && ne(g) == 105
        @test iscompletemultipartite(g, retrievepartitions(15, 15))
        g = turan_graph(10, 1)
        @test nv(g) == 10 && ne(g) == 0
        @test iscompletemultipartite(g, retrievepartitions(10, 1))
        @test_throws DomainError turan_graph(3,0)
        @test_throws DomainError turan_graph(0,4)
        @test_throws DomainError turan_graph(3,4)
        @test_throws DomainError turan_graph(-1,5)
        @test_throws DomainError turan_graph(3,-6)
    end

    @testset "Star Graphs" begin
        g = @inferred(star_graph(5))
        @test nv(g) == 5 && ne(g) == 4
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = star_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = star_graph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = star_graph(2)
        @test nv(g) == 2 && ne(g) == 1
        g = @inferred star_graph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127 - 1
        @test eltype(g) == Int8
    end

    @testset "Star DiGraphs" begin
        g = @inferred(star_digraph(5))
        @test nv(g) == 5 && ne(g) == 4
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = star_digraph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = star_digraph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = star_digraph(2)
        @test nv(g) == 2 && ne(g) == 1
        @test first(edges(g)) == Edge(1, 2) # edges should point outwards from vertex 1
        g = @inferred star_digraph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127 - 1
        @test eltype(g) == Int8

    end

    @testset "Path DiGraphs" begin
        g = @inferred(path_digraph(5))
        @test nv(g) == 5 && ne(g) == 4
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = path_digraph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = path_digraph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = @inferred path_digraph(Int8(127))
        @test nv(g) == 127 && ne(g) == 126
        @test eltype(g) == Int8

    end

    @testset "Path Graphs" begin
        g = @inferred(path_graph(5))
        @test nv(g) == 5 && ne(g) == 4
        @test isvalid_simplegraph(g)
        # tests for extreme values
        g = path_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = path_graph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = @inferred path_graph(Int8(127))
        @test nv(g) == 127 && ne(g) == 126
        @test eltype(g) == Int8
    end

    @testset "Cycle DiGraphs" begin
        g = @inferred(cycle_digraph(5))
        @test nv(g) == 5 && ne(g) == 5
        @test isvalid_simplegraph(g)
         # tests for extreme values
        g = cycle_digraph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = cycle_digraph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = cycle_digraph(2)
        @test nv(g) == 2 && ne(g) == 2
        g = @inferred cycle_digraph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127
        @test eltype(g) == Int8
    end

    @testset "Cycle Graphs" begin
        g = @inferred(cycle_graph(5))
        @test nv(g) == 5 && ne(g) == 5
        @test isvalid_simplegraph(g)
         # tests for extreme values
        g = cycle_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = cycle_graph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = cycle_graph(2)
        @test nv(g) == 2 && ne(g) == 1
        g = @inferred cycle_graph(Int8(127))
        @test nv(g) == 127 && ne(g) == 127
        @test eltype(g) == Int8
    end

    @testset "Wheel DiGraphs" begin
        g = @inferred(wheel_digraph(5))
        @test nv(g) == 5 && ne(g) == 8
        @test isvalid_simplegraph(g)
          # tests for extreme values
        g = wheel_digraph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = wheel_digraph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = wheel_digraph(2)
        @test nv(g) == 2 && ne(g) == 1
        g = wheel_digraph(3)
        @test nv(g) == 3 && ne(g) == 4
        g = @inferred wheel_digraph(Int8(127))
        @test nv(g) == 127 && ne(g) == 2 * 126
        @test eltype(g) == Int8
    end

    @testset "Wheel Graphs" begin
        g = @inferred(wheel_graph(5))
        @test nv(g) == 5 && ne(g) == 8
        @test isvalid_simplegraph(g)
          # tests for extreme values
        g = wheel_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = wheel_graph(1)
        @test nv(g) == 1 && ne(g) == 0
        g = wheel_graph(2)
        @test nv(g) == 2 && ne(g) == 1
        g = wheel_graph(3)
        @test nv(g) == 3 && ne(g) == 3
        g = @inferred wheel_graph(Int8(127))
        @test nv(g) == 127 && ne(g) == 2 * 126
        @test eltype(g) == Int8
    end

    @testset "grids" begin
        g = @inferred(grid([3, 3, 4]))
        @test nv(g) == 3 * 3 * 4
        @test ne(g) == 75
        @test Δ(g) == 6
        @test δ(g) == 3
        @test isvalid_simplegraph(g)

        g2 = @inferred(grid((3, 3, 4)))
        @test g2 == g

        g = @inferred(grid([3, 3, 4], periodic=true))
        @test nv(g) == 3 * 3 * 4
        @test ne(g) == 108
        @test Δ(g) == 6
        @test δ(g) == 6
        @test isvalid_simplegraph(g)

        g2 = @inferred(grid((3, 3, 4), periodic=true))
        @test g2 == g
    end

    @testset "Clique Graphs" begin
        g = @inferred(clique_graph(3, 5))
        @test nv(g) == 15 && ne(g) == 20
        @test g[1:3] == complete_graph(3)
        @test isvalid_simplegraph(g)
    end

    @testset "Binary Trees" begin
        g = @inferred(crosspath(3, binary_tree(2)))
        # f = Vector{Vector{Int}}[[2 3 4];
        # [1 5];
        # [1 6];
        # [1 5 6 7];
        # [2 4 8];
        # [3 4 9];
        # [4 8 9];
        # [5 7];
        # [6 7]
        # ]
        I = [1, 1, 1, 2, 2, 3, 3, 4, 4, 4, 4, 5, 5, 5, 6, 6, 6, 7, 7, 7, 8, 8, 9, 9]
        J = [2, 3, 4, 1, 5, 1, 6, 1, 5, 6, 7, 2, 4, 8, 3, 4, 9, 4, 8, 9, 5, 7, 6, 7]
        V = ones(Int, length(I))
        Adj = sparse(I, J, V)
        @test Adj == sparse(g)
        @test isvalid_simplegraph(g)
        @test_throws DomainError binary_tree(Int8(8))

        g = @inferred(double_binary_tree(3))
        # [[3, 2, 8]
        # [4, 1, 5]
        # [1, 6, 7]
        # [2]
        # [2]
        # [3]
        # [3]
        # [10, 9, 1]
        # [11, 8, 12]
        # [8, 13, 14]
        # [9]
        # [9]
        # [10]
        # [10]]
        I = [1, 1, 1, 2, 2, 2, 3, 3, 3, 4, 5, 6, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 12, 13, 14]
        J = [3, 2, 8, 4, 1, 5, 1, 6, 7, 2, 2, 3, 3, 10, 9, 1, 11, 8, 12, 8, 13, 14, 9, 9, 10, 10]
        V = ones(Int, length(I))
        Adj = sparse(I, J, V)
        @test Adj == sparse(g)
        @test isvalid_simplegraph(g)
    end

    @testset "Roach Graphs" begin
        rg3 = @inferred(roach_graph(3))
        # [3]
        # [4]
        # [1, 5]
        # [2, 6]
        # [3, 7]
        # [4, 8]
        # [9, 8, 5]
        # [10, 7, 6]
        # [11, 10, 7]
        # [8, 9, 12]
        # [9, 12]
        # [10, 11]
        I = [1, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7, 7, 8, 8, 8, 9, 9, 9, 10, 10, 10, 11, 11, 12, 12]
        J = [3, 4, 1, 5, 2, 6, 3, 7, 4, 8, 9, 8, 5, 10, 7, 6, 11, 10, 7, 8, 9, 12, 9, 12, 10, 11]
        V = ones(Int, length(I))
        Adj = sparse(I, J, V)
        @test Adj == sparse(rg3)
        @test isvalid_simplegraph(rg3)
    end

    function isladdergraph(g)
      n = nv(g)÷2
      !(degree(g, 1) == degree(g, n) == degree(g, n+1) == degree(g, 2*n) == 2) && return false
      !(has_edge(g, 1, n+1) && has_edge(g, n, 2*n)) && return false
      !(has_edge(g, 1, 2) && has_edge(g, n, n-1) && has_edge(g, n+1, n+2) && has_edge(g, 2*n, 2*n-1) ) && return false
      for i in 2:(n-1)
        !(degree(g, i) == 3 && degree(g, n+i) == 3) && return false
        !(has_edge(g, i, i%n +1) && has_edge(g, i, i+n)) && return false
        !(has_edge(g, n+i,n+(i%n +1)) && has_edge(g, n+i, i)) && return false
      end
      return true
    end

    @testset "Ladder Graphs" begin
        g = @inferred(ladder_graph(5))
        @test nv(g) == 10 && ne(g) == 13
        @test isvalid_simplegraph(g)
        @test isladdergraph(g)
        # tests for extreme values
        g = ladder_graph(0)
        @test nv(g) == 0 && ne(g) == 0
        g = ladder_graph(1)
        @test nv(g) == 2 && ne(g) == 1
        g = @inferred ladder_graph(-1)
        @test nv(g) == 0 && ne(g) == 0
        g = ladder_graph(Int8(63))
        @test nv(g) == 126 && ne(g) == 187
        @test eltype(g) == Int8
        @test isladdergraph(g)
        # tests for errors
        @test_throws InexactError ladder_graph(Int8(64))
    end

    function iscircularladdergraph(g)
      n = nv(g)÷2
      for i in 1:n
        !(degree(g, i) == 3 && degree(g, n+i) == 3) && return false
        !(has_edge(g, i, i%n +1) && has_edge(g, i, i+n)) && return false
        !(has_edge(g, n+i,n+(i%n +1)) && has_edge(g, n+i, i)) && return false
      end
      return true
    end

    @testset "Circular Ladder Graphs" begin
        g = @inferred(circular_ladder_graph(5))
        @test nv(g) == 10 && ne(g) == 15
        @test isvalid_simplegraph(g)
        @test iscircularladdergraph(g)
        # tests for extreme values
        g = circular_ladder_graph(3)
        @test nv(g) == 6 && ne(g) == 9
        @test iscircularladdergraph(g)
        g = circular_ladder_graph(Int8(63))
        @test nv(g) == 126 && ne(g) == 189
        @test eltype(g) == Int8
        @test iscircularladdergraph(g)
        # tests for errors
        @test_throws InexactError circular_ladder_graph(Int8(64))
        @test_throws DomainError circular_ladder_graph(-1)
        @test_throws DomainError circular_ladder_graph(0)
        @test_throws DomainError circular_ladder_graph(1)
        @test_throws DomainError circular_ladder_graph(2)
    end

    # checking that the nodes are organized correctly
    # see the docstring implementation notes for lollipop_graph
    function isbarbellgraph(g, n1, n2)
      nv(g) != n1+n2 && return false
      ne(g) != n1*(n1-1)÷2+n2*(n2-1)÷2 +1 && return false
      for i in 1:n1
        for j in (i+1):n1
          !has_edge(g, i, j) && return false
        end
      end

      for i in n1 .+ 1:n2
        for j in (i+1):(n1+n2)
          !has_edge(g, i, j) && return false
        end
      end

      !has_edge(g, n1, n1+1) && return false
      return true
    end

    @testset "Barbell Graphs" begin
        g = @inferred(barbell_graph(5, 6))
        @test nv(g) == 11 && ne(g) == 26
        @test isvalid_simplegraph(g)
        @test isbarbellgraph(g, 5, 6)
        g = barbell_graph(Int8(5), Int8(6))
        @test nv(g) == 11 && ne(g) == 26
        @test isvalid_simplegraph(g)
        @test isbarbellgraph(g, 5, 6)
        # extreme values
        g = barbell_graph(1, 5)
        @test nv(g) == 6 && ne(g) == 11
        @test isbarbellgraph(g, 1, 5)
        g = barbell_graph(5, 1)
        @test nv(g) == 6 && ne(g) == 11
        @test isbarbellgraph(g, 5, 1)
        g = barbell_graph(1, 1)
        @test nv(g) == 2 && ne(g) == 1
        @test isbarbellgraph(g, 1, 1)
        @test_throws OverflowError barbell_graph(Int8(100), Int8(50))
        @test_throws DomainError barbell_graph(1, 0)
        @test_throws DomainError barbell_graph(0, 1)
        @test_throws DomainError barbell_graph(-1, -1)
    end

    # checking that the nodes are organized correctly
    # see the docstring implementation notes for lollipop_graph
    function islollipopgraph(g, n1, n2)
      nv(g) != n1+n2 && return false
      ne(g) != n1*(n1-1)÷2+n2 && return false
      for i in 1:n1
        for j in (i+1):n1
          !has_edge(g, i, j) && return false
        end
      end

      for i in n1 .+ 1:(n2-1)
        !has_edge(g, i, i+1) && return false
      end

      !has_edge(g, n1, n1+1) && return false

      return true
    end

    @testset "Lollipop Graphs" begin
        g = @inferred(lollipop_graph(3, 5))
        @test nv(g) == 8 && ne(g) == 8
        @test isvalid_simplegraph(g)
        @test islollipopgraph(g, 3, 5)
        g = lollipop_graph(Int8(7), Int8(6))
        @test nv(g) == 13 && ne(g) == 27
        @test isvalid_simplegraph(g)
        @test islollipopgraph(g, 7, 6)
        # extreme values
        g = lollipop_graph(1, 3)
        @test nv(g) == 4 && ne(g) == 3
        @test islollipopgraph(g, 1, 3)
        g = lollipop_graph(3, 1)
        @test nv(g) == 4 && ne(g) == 4
        @test islollipopgraph(g, 3, 1)
        g = lollipop_graph(1, 1)
        @test nv(g) == 2 && ne(g) == 1
        @test islollipopgraph(g, 1, 1)
        @test_throws OverflowError lollipop_graph(Int8(100), Int8(50))
        @test_throws DomainError lollipop_graph(1, 0)
        @test_throws DomainError lollipop_graph(0, 1)
        @test_throws DomainError lollipop_graph(-1, -1)
    end
end


text:
"""This module implements functions to access solution and make prediction."""
# Author: Henri Gérard <hgerard.proy@gmail.com>
# License: MIT

function getsolution(x, ambiguity, nbfeatures)
    if ambiguity == "KLdivergence"
        return x[3:3+nbfeatures-1]
    elseif ambiguity ==  "wasserstein"
        return x[2:2+nbfeatures-1]
    elseif ambiguity ==  "entropic"
        return x[2:2+nbfeatures-1]
    end
end

function pred(data,θ, regression::RegressionModel)
    error("'error_loss' not defined for $(typedof(s))")
end

function pred(data, θ, regression::LinearRegression)
    pred = zeros(size(data)[1])
    for i in 1:length(pred)
        pred[i] = θ'*data[i,1:end-1]
    end
    return pred
end

function pred(data, θ, regression::LogisticRegression)
    pred = zeros(size(data)[1])
    for i in 1:length(data[:,end])
        if (1/(1+exp(θ'*data[i,1:end-1])) > 0.5)
            pred[i] = -1
        end
        if (1/(1+exp(-θ'*data[i,1:end-1])) > 0.5)
            pred[i] = 1
        end
    end
    return pred
end

function positive_rate(x, data)
    aux = zeros(size(data)[1])
    for i in 1:size(data)[1]
        aux[i] = 1/(1+exp(-x'*data[i,1:end-1]))
    end
    return aux
end

max_stars_repo_name: Henri-Gerard/robox

max_issues_repo_name: Henri-Gerard/robox

max_forks_repo_name: Henri-Gerard/robox


text:
using Distributions
using StatsBase
using Random
using DataFrames
using StatsPlots
using BenchmarkTools;

function rate_to_proportion(r::Float64,t::Float64)
    1-exp(-r*t)
end;

@enum InfectionStatus Susceptible Infected Recovered

function sir_abm(u,p,t)
    du = deepcopy(u)
    (β,c,γ,δt) = p
    N = length(u)
    for i in 1:N # loop through agents
        # If recovered
        if u[i]==Recovered continue
        # If susceptible
        elseif u[i]==Susceptible
            ncontacts = rand(Poisson(c*δt))
            du[i]=Susceptible
            ncontacts = rand(Poisson(c*δt))
            while ncontacts > 0
                j = sample(1:N)
                if j==i
                    continue
                end
                a = u[j]
                if a==Infected && rand() < β
                    du[i] = Infected
                    break
                end
                ncontacts -= 1
            end
        # If infected
    else u[i]==Infected
            if rand() < γ
                du[i] = Recovered
            end
        end
    end
    du
end;

function sir_abm!(du,u,p,t)
    (β,c,γ,δt) = p
    N = length(u)
    # Initialize du to u
    for i in 1:N
        du[i] = u[i]
    end
    for i in 1:N # loop through agents
        # If recovered
        if u[i]==Recovered
            continue
        # If susceptible
        elseif u[i]==Susceptible
            ncontacts = rand(Poisson(c*δt))
            while ncontacts > 0
                j = sample(1:N)
                if j==i
                    continue
                end
                a = u[j]
                if a==Infected && rand() < β
                    du[i] = Infected
                    break
                end
                ncontacts -= 1
            end
        # If infected
        else u[i]==Infected
            if rand() < γ
                du[i] = Recovered
            end
        end
    end
    nothing
end;

δt = 0.1
nsteps = 400
tf = nsteps*δt
tspan = (0.0,nsteps)
t = 0:δt:tf;

β = 0.05
c = 10.0
γ = rate_to_proportion(0.25,δt)
p = [β,c,γ,δt];

N = 1000
I0 = 10
u0 = Array{InfectionStatus}(undef,N)
for i in 1:N
    if i <= I0
        s = Infected
    else
        s = Susceptible
    end
    u0[i] = s
end;

Random.seed!(1234);

susceptible(x) = count(i == Susceptible for i in x)
infected(x) = count(i == Infected for i in x)
recovered(x) = count(i == Recovered for i in x);

function sim(u0,nsteps,dt)
    u = copy(u0)
    t = 0.0
    ta = []
    Sa = []
    Ia = []
    Ra =[]
    push!(ta,t)
    push!(Sa,susceptible(u))
    push!(Ia,infected(u))
    push!(Ra,recovered(u))
    for i in 1:nsteps
        u=sir_abm(u,p,t)
        t = t + dt
        push!(ta,t)
        push!(Sa,susceptible(u))
        push!(Ia,infected(u))
        push!(Ra,recovered(u))
    end
    DataFrame(t=ta,S=Sa,I=Ia,R=Ra)
end;

function sim!(u0,nsteps,dt)
    u = copy(u0)
    du = copy(u0)
    t = 0.0
    ta = []
    Sa = []
    Ia = []
    Ra =[]
    push!(ta,t)
    push!(Sa,susceptible(u))
    push!(Ia,infected(u))
    push!(Ra,recovered(u))
    for i in 1:nsteps
        sir_abm!(du,u,p,t)
        u,du = du,u
        t = t + dt
        push!(ta,t)
        push!(Sa,susceptible(u))
        push!(Ia,infected(u))
        push!(Ra,recovered(u))
    end
    DataFrame(t=ta,S=Sa,I=Ia,R=Ra)
end;

df_abm = sim(u0,nsteps,δt);

df_abm! = sim!(u0,nsteps,δt);

@df df_abm plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number",
    title="New state")

@df df_abm! plot(:t,
    [:S :I :R],
    label=["S" "I" "R"],
    xlabel="Time",
    ylabel="Number",
    title="In-place")

@benchmark sim(u0,nsteps,δt)

@benchmark sim!(u0,nsteps,δt)


text:
module PPLM

using Flux
using Zygote: @adjoint
using Flux: onehotbatch
using Flux.Data: DataLoader
using Flux.Optimise: update!
using DataFrames: DataFrame
using Transformers.HuggingFace
using DataDeps
using BSON
using CUDA
using Parameters: @with_kw
using MLDataUtils
using ProgressMeter
using Pkg.Artifacts

const ARTIFACTS_TOML = joinpath(@__DIR__, "Artifacts.toml")

const SMALL_CONST = Float32(1.0e-15)

abstract type PretrainedTokenizer end

include("artifacts.jl")
include("bow.jl")
include("discriminator.jl")
include("data_preprocess.jl")
include("discrim_train.jl")
include("tokenizer_datadeps.jl")
include("tokenizer.jl")
include("generate.jl")
include("utils.jl")
include("pplm.jl")

export load_pretrained_tokenizer, GPT2, encode, decode, detokenize, tokenize, load_data, load_cached_data, data_preprocess, ClassifierHead, DisciminatorV1, DiscriminatorV2

"""
    get_gpt2()

Function to load gpt2 lmheadmodel along with the `tokenizer`.
"""
function get_gpt2()
    model = hgf"gpt2:lmheadmodel"
    tokenizer = load_pretrained_tokenizer(GPT2)
    return tokenizer, model
end

"""
    get_gpt2_medium()

Function to load gpt2-medium lmhead model along with the `tokenizer`.

**Note**: In case this function gives error of permission denied, try changing the file permissions for the Artifacts.toml file of Transformers.jl package (as it is read only by default) under the `src/huggingface` folder.
"""
function get_gpt2_medium()
    model = hgf"gpt2-medium:lmheadmodel"
    tokenizer = load_pretrained_tokenizer(GPT2)
    return tokenizer, model
end

device_id=-1
device=Flux.cpu

"""
    set_device(d_id=0)

Function to set cuda device if available and also to disallow scalar operations
"""
function set_device(d_id=0)
    global device_id, device
    if CUDA.has_cuda()       # Check if CUDA is available
        device = Flux.gpu
        device_id = d_id
        CUDA.device!(d_id)
        @info "CUDA is on"
        CUDA.allowscalar(false)
    end
end

function __init__()
    init_tokenizer_datadeps()
    set_device()
end

end


text:
module ClassicalSystems

    export integrate,lyapunov_spectrum,lyapunov_exponent,ClassicalSystem

    using DiffEqCallbacks,OrdinaryDiffEq,DiffEqBase
    using ProgressMeter
    using RecursiveArrayTools
    using LinearAlgebra
    using Statistics

    """
    ```julia
    abstract struct ClassicalSystem
    ```
    This abstract object represents a classical system that may be passed to multiple functions
    in this module. To generate a system, use [`ClassicalDicke.ClassicalDickeSystem`](@ref Main.ClassicalDicke.ClassicalDickeSystem)
    or use [`ClassicalLMG.ClassicalLMGSystem`](@ref Main.ClassicalLMG.ClassicalLMGSystem).
    """
    abstract type ClassicalSystem end

    function parameters(system::ClassicalSystem)
        error("Please define a method ClassicalSystems.parameters(system::$(typeof(system)) that returns the list
        of parameters for this system")
    end
    function step(system::ClassicalSystem)
        error("Please define a method ClassicalSystems.step(system::$(typeof(system)) that returns ode system
        of equations")
    end
    function out_of_bounds(system::ClassicalSystem)
        error("Please define a method ClassicalSystems.out_of_bounds(system::$(typeof(system)) that returns a function
        out_of_bounds(u,_,t) that determines if u is out of bounds")
    end
    function varnames(system::ClassicalSystem)
        error("Please define a method ClassicalSystems.varnames(system::$(typeof(system)) that returns a function
        varnames(u,_,t) that reurns the coordinate names in order")
    end
    function nvars(system::ClassicalSystem,vars...)

        vs=[]
        d=Dict(zip(varnames(system),1:length(varnames(system))))
        for v in vars
            push!(vs,d[v])
        end
        return vs
    end
    eye(n)=Matrix(Diagonal(ones(n)))

    """
    ```julia
    function integrate(system::ClassicalSystem,
                       u₀::AbstractVector{<:Real},
                       t::Real;
                       t₀=0.0::Real,
                       tol=1e-12::Real,
                       get_fundamental_matrix::Bool=false,
                       integrator_alg=TsitPap8(),
                       use_big_numbers::Bool=false,
                       integate_backwards::Bool=false,
                       fundamental_matrix_initial=nothing,
                       kargs...)
    ```
    This function integrates initial condition `u₀` from `t₀` to `t` under the Hamiltonian system determined by `system`, returning an instance of [`OrdinaryDiffEq.ODESolution`](https://diffeq.sciml.ai/dev/basics/solution/).

    # Arguments
    - `system` is an instance of [`ClassicalSystems.ClassicalSystem`](@ref).
    - `u₀` is an array which codifies the initial condition `[Q,q,P,p]` for `ClassicalDicke` and `[Q,P]` for `ClassicalLMG`.
    - `t` is the start time of the integration.
    # Keyword arguments
    - `t₀` is the start of the integration (defaults to `t₀ = 0.0`) [`ClassicalSystems.ClassicalSystem`](@ref).
    - `tol` is the tolerance for the integration, which determines both `abstol` and `reltol` in [`OrdinaryDiffEq.solve`](https://diffeq.sciml.ai/stable/basics/common_solver_opts/)
    - `get_fundamental_matrix` determines whether to also compute the fundametal matrix of the system. If `true`, the result at each time is an [`ArrayPartition(x,Φ)`](https://diffeq.sciml.ai/stable/features/diffeq_arrays/#ArrayPartitions), so that `x=result.x[2]` retrieves the coordinate and `Ψ=result.x[2]` retrieves the fundamental matrix. Default is `false`. Note that the integration is consideribly slowed down if this parameter is set to `true`.
    - `integrator_alg` is the integration algorithm to use. Defaults to `TsitPap8` (Tsitouras-Papakostas 8/7 Runge-Kutta method). See [the `DifferentialEquations` documentation](https://diffeq.sciml.ai/stable/solvers/ode_solve/#Full-List-of-Methods) for other options.
    - `use_big_numbers` forces the integration to be performed with `BigFloat` instead of `Float`, allowing for infinite numerical precision, but hindering speed substantially. Defaults to `false`.
    - `integate_backwards` tells the integrator to integrate back in time, from `-t₀` to `-t`. Defaults to  `false`.
    - `fundamental_matrix_initial` is the initial value if `get_fundamental_matrix` is `true`. The default value is the identity matrix.
    - Additional `kargs` are passed to [`OrdinaryDiffEq.solve`](https://diffeq.sciml.ai/stable/basics/common_solver_opts/).
    """
    function integrate(system::ClassicalSystem,
                            u₀::AbstractVector{<:Real},
                            t::Real;
                            t₀=0.0::Real,
                            tol=1e-12::Real,
                            get_fundamental_matrix::Bool=false,
                            integrator_alg=TsitPap8(),
                            use_big_numbers::Bool=false,
                            integate_backwards::Bool=false,
                            fundamental_matrix_initial::Union{Nothing,AbstractMatrix{<:Real}}=nothing,
                            kargs...)
        if t<t₀
           error("Use integate_backwards=true instead of negative time")
        end
        if integate_backwards==true
            integate_backwards=-1
        else
            integate_backwards=1
        end
        if length(u₀)!=length(varnames(system))
            error("u₀ should have length $(length(varnames(system))), which is the dimension of the phase space of system")
        end
        if out_of_bounds(system)(u₀,nothing,0)
            error("The initial condition u₀=$u₀ is out of bounds")
        end

        u₀=u₀
        if use_big_numbers
            u₀=big.(u₀)
        end
        sys=nothing

        if get_fundamental_matrix
            n=length(u₀)

            if fundamental_matrix_initial==nothing
                fundamental_matrix_initial=eye(n)
            end
            cacheJ=eye(n)
            if use_big_numbers
                fundamental_matrix_initial=big.(fundamental_matrix_initial)
                cacheJ=big.(cacheJ)
            end

            u₀ = ArrayPartition(u₀,fundamental_matrix_initial)
            function sistemaconmatriz(du,u,p,t)
                step(system)(du.x[1],u.x[1],p,t)
                step(system).jac(cacheJ,u.x[1],p,t)
                mul!(du.x[2], cacheJ, u.x[2])
            end
            sys = ODEFunction(sistemaconmatriz)

        else
            sys=step(system)
        end

        tspan = (t₀,t)
        params=[parameters(system);integate_backwards]
        prob = ODEProblem(sys,u₀,tspan,params)
        s= solve(prob;alg=integrator_alg,isoutofdomain=out_of_bounds(system),progress=false,dtmin=tol,force_dtmin=true,abstol=tol,reltol=tol,kargs...)
        return s
    end

    """
    ```julia
    function lyapunov_exponent(system::ClassicalSystem,
        x::AbstractVector{<:Real},
        t::Real = 100; kargs...)
    ```
    Returns the maximal [Lyapunov exponent](https://en.wikipedia.org/wiki/Lyapunov_exponent) for `x` by
    calling `maximum(lyapunov_spectrum(system, x, t; kargs...))`. Same arguments as [`lyapunov_spectrum`](@ref) apply.
    """
    function lyapunov_exponent(system::ClassicalSystem,
        x::AbstractVector{<:Real},
        t::Real = 100;kargs...)
        return maximum(lyapunov_spectrum(system,x,t;kargs...))

     end
    """
    ```julia
    function lyapunov_spectrum(system::ClassicalSystem,
        x::AbstractVector{<:Real},
        t::Real = 100;
        λtol::Real = 1e-5,
        λreltol::Real = 1e-3,
        maxiters::Integer = 100,
        verbose::Bool = true,
        kargs...)
    ```
    Calculates the [Lyapunov spectrum](https://en.wikipedia.org/wiki/Lyapunov_exponent#Definition_of_the_Lyapunov_spectrum) at
    point `x`. The algorithm used is given in Fig. 3.7 of Ref. [Parker1989Book](@cite), which basically integrates the fundamental
    matrix and renormalizes in intervals of length `t`, until the oscillations of the Lyapunov exponets go bellow `λtol + λ*λreltol`.
    # Arguments
    - `system` is an instance of [`ClassicalSystems.ClassicalSystem`](@ref).
    - `x` is an array which codifies the initial condition `[Q,q,P,p]` for `ClassicalDicke` and `[Q,P]` for `ClassicalLMG`.
    - `t` is the time intervals for which to integrate. This value must not be too big, or the fundamental matrix may
      get too large and errors will occur.  Default value is `100`, which is fine most of the time.
    # Keyword arguments
    - `λtol` is the absolute tolerance (``E_a`` in Fig. 3.7 of Ref. [Parker1989Book](@cite)). Default is `1e-5`.
    - `λreltol` is the relative tolerance (``E_r`` in Fig. 3.7 of Ref. [Parker1989Book](@cite)). Default is `1e-3`.
    - `maxiters` is the maximum number of iterations. Default is `100`.
    - `verbose` is a boolean toggling the progress bar (default is `true`, which enables the progress bar).
    - `kargs` are redirected to [`integrate`](@ref integrate(::ClassicalSystem,::AbstractVector{<:Real},::Real)).
    """
    function lyapunov_spectrum(system::ClassicalSystem, x::AbstractVector{<:Real},
        t::Real = 100;
        λtol::Real = 1e-4,
        λreltol::Real = 1e-3,
        maxiters::Integer = 100,
        verbose::Bool = true,
        kargs...)
        n = length(varnames(system))
        u=eye(n)
        λ=zeros(n)
        λ_old=zeros(n)
        sums=zeros(n)
        v=zeros(n,n)
        if verbose
            prog=ProgressUnknown("Integrating")
        end

        for k in 1:maxiters
            λ_old .= λ
            x,δx=integrate(system,x,t;kargs...,fundamental_matrix_initial=u,save_everystep=false,get_fundamental_matrix=true,verbose=verbose).u[end].x

            for i in 1:n
                v[:,i] .= @view δx[:,i]
                for j in 1:(i-1)
                    v[:,i] .-= dot((@view v[:,i]),(@view u[:,j])) .* (@view u[:,j])
                end
                nrmv = norm(@view v[:,i])
                u[:,i] .= (@view v[:,i]) ./ nrmv
                sums[i] += log(nrmv)
                λ[i] = sums[i]/(k*t)
                err= norm(λ_old - λ)
                if err < λreltol*norm(λ) + λtol
                    return λ
                end
                if verbose
                    prog.desc="max_λ_estimate = $(round(maximum(λ),digits=6)), error = 10^$(round(log10(err),digits=1)), iterations ="
                    next!(prog)
                end
            end
        end
        error("Maximum iterations")
    end
end


text:
#=
LineSearch.jl

    Provides line search methods, such as backtracking, to determine the optimal
    step size in optimization and root-finding algorithms.

@author: Quint Wiersma <q.wiersma@vu.nl>

@date: 2022/02/02
=#

module LineSearch

using LinearAlgebra

export
    # Structs
    BackTrack,
    # Backtracking
    backtrack!

# Include programs
include("backtrack.jl")

end


text:
function ConModel(A,b,C,d,XLB,XUB,y;Cuts = [])
    m,n = size(A);
    nObj = size(C)[1];
    MyModel = Model();
    @variable(MyModel, XLB[i] <= X[i=1:n] <= XUB[i])
    @variable(MyModel, Z[i=1:nObj]>=0)
    @objective(MyModel,Min,sum(1*(sum(C[j,i]*X[i] for i = 1:n) + d[j]) for j=1:nObj))
    @constraintref Ucons[1:m];
    for i = 1:m
        Ucons[i] = @constraint(MyModel,sum(A[i,j]*X[j] for j = 1:n) >= b[i] )
    end
    for j = 1:nObj
        @constraint(MyModel, Z[j] == sum(C[j,i]*X[i] for i = 1:n) + d[j] )
    end
    for j = 2:nObj
        @constraint(MyModel, y[j]*(sum(C[1,i]*X[i] for i = 1:n) + d[1]) == y[1]*(sum(C[j,i]*X[i] for i = 1:n) + d[j]) )
        # @constraint(MyModel, y[j]*Z[1] == y[1]*Z[j] )
    end
    for c = 1:length(Cuts)
        @constraint(MyModel,sum(Cuts[c].Coeff[i]*(sum(C[i,j]*X[j] for j = 1:n )+d[i])  for i = 1:nObj) >= Cuts[c].RHS)
    end
    return MyModel;

end


text:
using KNITRO
using Base.Test
using Compat

 ## Solve a small QCQP (quadratically constrained quadratic programming)
 #  test problem.
 #
 #  min   1000 - x0^2 - 2 x1^2 - x2^2 - x0 x1 - x0 x2
 #  s.t.  8 x0 + 14 x1 + 7 x2 - 56 = 0
 #        x0^2 + x1^2 + x2^2 - 25 >= 0
 #        x0 >= 0, x1 >= 0, x2 >= 0
 #
 #  The start point (2, 2, 2) converges to the minimum at (0, 0, 8),
 #  with final objective = 936.0.  From a different start point,
 #  KNITRO may converge to an alternate local solution at (7, 0, 0),
 #  with objective = 951.0.
 ##

function eval_f(x::Vector{Float64})
    1000.0 - x[1]^2 - 2.0*x[2]^2 - x[3]^2 - x[1]^2 - x[1]*x[3]
end

function eval_g(x::Vector{Float64}, cons::Vector{Float64})
    cons[1] = 8.0*x[1] + 14.0*x[2] + 7.0*x[3] - 56.0
    cons[2] = x[1]^2 + x[2]^2 + x[3]^2 - 25.0
end

function eval_grad_f(x::Vector{Float64}, grad::Vector{Float64})
    grad[1] = -2.0*x[1] - x[2] - x[3]
    grad[2] = -4.0*x[2] - x[1]
    grad[3] = -2.0*x[3] - x[1]
end

function eval_jac_g(x::Vector{Float64}, jac::Vector{Float64})
    #---- GRADIENT OF THE FIRST CONSTRAINT, c[0].
    jac[1] =  8.0
    jac[2] = 14.0
    jac[3] =  7.0

    #---- GRADIENT OF THE SECOND CONSTRAINT, c[1].
    jac[4] = 2.0*x[1]
    jac[5] = 2.0*x[2]
    jac[6] = 2.0*x[3]
end

function eval_h(x::Vector{Float64}, lambda::Vector{Float64},
                sigma::Float64, hess::Vector{Float64})
    hess[1] = -2.0*sigma + 2.0*lambda[2]
    hess[2] = -1.0*sigma
    hess[3] = -1.0*sigma
    hess[4] = -4.0*sigma + 2.0*lambda[2]
    hess[5] = -2.0*sigma + 2.0*lambda[2]
end

function eval_hv(x::Vector{Float64}, lambda::Vector{Float64},
                 sigma::Float64, hv::Vector{Float64})
    hv1 = (-2.0*sigma + 2.0*lambda[2])*hv[1] - sigma*hv[2] - sigma*hv[3]
    hv2 = -sigma*hv[1] + (-4.0*sigma + 2.0*lambda[2])*hv[2]
    hv3 = -sigma*hv[1] + (-2.0*sigma + 2.0*lambda[2])*hv[3]
    hv[1] = hv1
    hv[2] = hv2
    hv[3] = hv3
end

objGoal = KTR_OBJGOAL_MINIMIZE
objType = KTR_OBJTYPE_QUADRATIC

n = 3
x_L = zeros(3)
x_U = fill(KTR_INFBOUND, 3)
x = [-2.0,1.0] # initial guess

m = 2
c_Type = [KTR_CONTYPE_LINEAR, KTR_CONTYPE_QUADRATIC]
c_L = zeros(2)
c_U = [0.0, KTR_INFBOUND]

jac_con = Int32[0,0,0,1,1,1]
jac_var = Int32[0,1,2,0,1,2]
hess_row = Int32[0,0,0,1,2]
hess_col = Int32[0,1,2,1,2]

x = [2.0,2.0,2.0] # initial guess

kp = createProblem()
@fact applicationReturnStatus(kp) --> :Uninitialized
setOption(kp, "outlev", "all")
setOption(kp, "hessopt", 1)
setOption(kp, "hessian_no_f", 1)
setOption(kp, "feastol", 1.0e-10)

@fact applicationReturnStatus(kp) --> :Uninitialized
@fact kp.eval_status --> @compat(Int32(0))

initializeProblem(kp, objGoal, objType, x_L, x_U, c_Type, c_L, c_U,
                  jac_var, jac_con, hess_row, hess_col; initial_x = x)
@fact applicationReturnStatus(kp) --> :Initialized
@fact kp.eval_status --> @compat(Int32(0))

#---- ALLOCATE ARRAYS FOR REVERSE COMMUNICATIONS OPERATION.
cons = Array(Float64, m)
objGrad = Array(Float64, n)
jac = Array(Float64, length(jac_con))
hess = Array(Float64, length(hess_row))
hessVector = Array(Float64, n)

#---- SOLVE THE PROBLEM.  IN REVERSE COMMUNICATIONS MODE, KNITRO
#---- RETURNS WHENEVER IT NEEDS MORE PROBLEM INFORMATION.  THE CALLING
#---- PROGRAM MUST INTERPRET KNITRO'S RETURN STATUS AND CONTINUE
#---- SUPPLYING PROBLEM INFORMATION UNTIL KNITRO IS COMPLETE.
while kp.status > 0
    @fact applicationReturnStatus(kp) --> anyof(:Initialized, :ReverseComms)
    solveProblem(kp, cons, objGrad, jac, hess, hessVector)
    if kp.status == KTR_RC_EVALFC
        #---- KNITRO WANTS obj AND c EVALUATED AT THE POINT x.
        kp.obj_val[1] = eval_f(kp.x)
        eval_g(kp.x,cons)
    elseif kp.status == KTR_RC_EVALGA
        #---- KNITRO WANTS objGrad AND jac EVALUATED AT THE POINT x.
        eval_grad_f(kp.x, objGrad)
        eval_jac_g(kp.x, jac)
    elseif kp.status == KTR_RC_EVALH
        #---- KNITRO WANTS hess EVALUATED AT THE POINT x.
        eval_h(kp.x, kp.lambda, 1.0, hess)
    elseif kp.status == KTR_RC_EVALH_NO_F
        #---- KNITRO WANTS hess EVALUATED AT THE POINT x
        #---- WITHOUT OBJECTIVE COMPONENT.
        eval_h(kp.x, kp.lambda, 0.0, hess)
    elseif kp.status == KTR_RC_EVALHV
        #---- KNITRO WANTS hessVector EVALUATED AT THE POINT x.
        eval_hv(kp.x, kp.lambda, 1.0, hessVector)
    elseif kp.status == KTR_RC_EVALHV_NO_F
        #---- KNITRO WANTS hessVector EVALUATED AT THE POINT x
        #---- WITHOUT OBJECTIVE COMPONENT.
        eval_hv(kp.x, kp.lambda, 0.0, hessVector)
    end
    #/*---- ASSUME THAT PROBLEM EVALUATION IS ALWAYS SUCCESSFUL.
    #*---- IF A FUNCTION OR ITS DERIVATIVE COULD NOT BE EVALUATED
    #*---- AT THE GIVEN (x, lambda), THEN SET evalStatus = 1 BEFORE
    #*---- CALLING KTR_solve AGAIN. */
    @fact applicationReturnStatus(kp) --> anyof(:ReverseComms, :Optimal)
    @fact kp.eval_status --> @compat(Int32(0))
end

# --- test optimal solutions ---
facts("Test optimal solutions") do
    @fact applicationReturnStatus(kp) --> :Optimal
    @fact kp.x --> roughly(
        [0.0, 0.0, 8.0], 1e-5)
    @fact kp.obj_val[1] --> roughly(936.0, 1e-5)
end


text:
#===============================================================================

The most naive way of computing n15 requires fourteen multiplications:

n × n × ... × n = n15

But using a "binary" method you can compute it in six multiplications:

n × n = n^2
n^2 × n^2 = n^4
n^4 × n^4 = n^8
n^8 × n^4 = n^12
n^12 × n^2 = n^14
n^14 × n = n^15

However it is yet possible to compute it in only five multiplications:

n × n = n^2
n^2 × n = n^3
n^3 × n^3 = n^6
n^6 × n^6 = n^12
n^12 × n^3 = n^15

We shall define m(k) to be the minimum number of multiplications to compute n^k;
for example m(15) = 5.

For 1 ≤ k ≤ 200, find ∑ m(k).

--------------------------------------------------------------------------------

To start with, this problem can be reformulated by asking how many additions are
necessary to reach a number starting from 1 and using only additions.

The general approach is bottom up. At every iteration we build a new set of sets
of numbers. An inner set contains all the numbers that can be together obtained
by 'iter' calculations.

An inner set is generated from a previous inner set by taking a copy of the
later, then by selecting one of its elements, by adding it to the initial
elements and adding the numbers so obtained to the new set. This must be done
for all elements of the initial set.

At the next iteration, it is only necessary to consider the new elements that
have been added when selecting the element to add to the other.

Because a same inner set can be obtained from different paths, each new inner
set is associated with the numbers selected to produce it.

Each time a new number is generated, the corresponding iteration is recorded
with it.

===============================================================================#

using Base.Test

function compute_nbr_steps(limit)
    nbr_steps = Array{Int}(limit)
    # nbr_steps[i] is the number of addition needed to compute i (i.e. it is
    # m(i)).

    fill!(nbr_steps, -1)
    sets = Dict{Set{Int}, Set{Int}}()
    # Key is all numbers in the set, Value is the set of the numbers added
    # during the previous iteration.

    const empty_value = Set{Int}()
    # empty_value is used as an indicator

    completion = 1
    iter = 0
    nbr_steps[1] = 0
    sets[Set{Int}(1)] = Set{Int}(1)

    while completion<limit
        iter += 1
        updated_sets = typeof(sets)()

        for (set, added_numbers) in sets, l in set, m in added_numbers
            n=l+m
            if n ≤ limit && n ∉ set
                updated_set = copy(set)
                push!(updated_set, n)

                nums = get(updated_sets, updated_set, empty_value)

                if nums == empty_value
                    nums = Set{Int}(n)
                else
                    push!(nums, n)
                end

                updated_sets[updated_set] = nums

                if nbr_steps[n]<0
                    nbr_steps[n] = iter
                    completion += 1
                    println("\t completion = $(completion)")
                    if completion == limit
                        return nbr_steps
                    end
                end
            end
        end
        sets = updated_sets
    end
end

result = compute_nbr_steps(200)

@test result[15] == 5
@test length(result) == 200

print("Euler 122: $(sum(result))")

max_stars_repo_name: dpieroux/euler

max_issues_repo_name: dpieroux/euler

max_forks_repo_name: dpieroux/euler


text:
export ActivationFunction
export Neurons

export forward, backward

############################################################
# Each activation function / neuron should be a sub type
# of ActivationFunction and should implement the following
# methods.
#
# - forward(backend :: Backend, neuron :: MyActivationFunction, output :: Blob)
#   Update the output blob by the activation function. Computation should be done in-place.
#   Note we call is "output" because an activation function is applied to the "output" of
#   a layer.
#
# - backward(backend :: Backend, neuron :: MyActivationFunction, output :: Blob, gradient :: Blob)
#   The gradient blob contains the gradient with respect to the output of the activation
#   function, update (IN-PLACE) the gradient to be with respect to the output of the activation
#   function.
############################################################

@compat abstract type ActivationFunction end

############################################################
# A module to hold built-in activation functions to avoid
# messy namespace
############################################################
module Neurons
using ..Mocha.ActivationFunction
# Identity
type Identity <: ActivationFunction
end

# Rectified-Linear: ReLU(eps)(x) = max(x,eps)
type ReLU <: ActivationFunction
  epsilon::Float64 # optional floor value, default zero
end
ReLU() = ReLU(0.0)

# Exponential: Exponential(x) = exp(x)
type Exponential <: ActivationFunction
end

# Leaky Rectified-Linear: LReLU(x) = x > 0 ? x : 0.01x
type LReLU <: ActivationFunction
end

# Sigmoid: Sigmoid(x) = 1 / (1 + exp(-x))
type Sigmoid <: ActivationFunction
end

# Sigmoid: Tanh(x) = (1 + exp(-2x)) / (1 + exp(-2x))
type Tanh <: ActivationFunction
end
end # module Neurons

############################################################
# Identity
############################################################
function forward(backend :: Backend, neuron :: Neurons.Identity, output :: Blob)
  # do nothing
end
function backward(backend :: Backend, neuron :: Neurons.Identity, output :: Blob, gradient :: Blob)
  # do nothing
end

############################################################
# Rectified-Linear
############################################################
function forward(backend :: CPUBackend, neuron :: Neurons.ReLU, output :: Blob)
  function _forward(output::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds output[i] = max(neuron.epsilon, output[i])
    end
  end
  _forward(output.data)
end
function backward(backend :: CPUBackend, neuron :: Neurons.ReLU, output :: Blob, gradient :: Blob)
  function _backward(output::AbstractArray, gradient::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds gradient[i] *= (output[i] > neuron.epsilon)
    end
  end
  _backward(output.data, gradient.data)
end

############################################################
# Leaky Rectified-Linear
############################################################
function forward(backend :: CPUBackend, neuron :: Neurons.LReLU, output :: Blob)
  function _forward(output::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds output[i] = output[i] > 0 ? output[i] : 0.01 * output[i]
    end
  end
  _forward(output.data)
end
function backward(backend :: CPUBackend, neuron :: Neurons.LReLU, output :: Blob, gradient :: Blob)
  function _backward(output::AbstractArray, gradient::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds gradient[i] *= ((output[i] > 0) + 0.01 * (output[i] <= 0))
    end
  end
  _backward(output.data, gradient.data)
end

############################################################
# Sigmoid
############################################################
function forward(backend :: CPUBackend, neuron :: Neurons.Sigmoid, output :: Blob)
  function _forward(output::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds output[i] = 1 / (1 + exp(-output[i]))
    end
  end
  _forward(output.data)
end
function backward(backend :: CPUBackend, neuron :: Neurons.Sigmoid, output :: Blob, gradient :: Blob)
  function _backward(output::AbstractArray, gradient::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds gradient[i] *= output[i] * (1-output[i])
    end
  end
  _backward(output.data, gradient.data)
end

############################################################
# Tanh
############################################################
function forward(backend :: CPUBackend, neuron :: Neurons.Tanh, output :: Blob)
  function _forward(output::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds output[i] = tanh(output[i])
    end
  end
  _forward(output.data)
end
function backward(backend :: CPUBackend, neuron :: Neurons.Tanh, output :: Blob, gradient :: Blob)
  function _backward(output::AbstractArray, gradient::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds gradient[i] *= (1 - output[i] * output[i])
    end
  end
  _backward(output.data, gradient.data)
end

############################################################
# Exponential
############################################################
function forward(backend :: CPUBackend, neuron :: Neurons.Exponential, output :: Blob)
  function _forward(output::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds output[i] = exp(output[i])
    end
  end
  _forward(output.data)
end
function backward(backend :: CPUBackend, neuron :: Neurons.Exponential, output :: Blob, gradient :: Blob)
  function _backward(output::AbstractArray, gradient::AbstractArray)
    @simd for i = 1:length(output)
      @inbounds gradient[i] *= output[i]
    end
  end
  _backward(output.data, gradient.data)
end