text:
import Base: uinttype

const MANTISSA_MASK = 0x000fffffffffffff
const EXP_MASK = 0x00000000000007ff

sign_mask(::Type{Float64}) =        0x8000_0000_0000_0000
exponent_mask(::Type{Float64}) =    0x7ff0_0000_0000_0000
exponent_one(::Type{Float64}) =     0x3ff0_0000_0000_0000
exponent_half(::Type{Float64}) =    0x3fe0_0000_0000_0000
significand_mask(::Type{Float64}) = 0x000f_ffff_ffff_ffff

sign_mask(::Type{Float32}) =        0x8000_0000
exponent_mask(::Type{Float32}) =    0x7f80_0000
exponent_one(::Type{Float32}) =     0x3f80_0000
exponent_half(::Type{Float32}) =    0x3f00_0000
significand_mask(::Type{Float32}) = 0x007f_ffff

sign_mask(::Type{Float16}) =        0x8000
exponent_mask(::Type{Float16}) =    0x7c00
exponent_one(::Type{Float16}) =     0x3c00
exponent_half(::Type{Float16}) =    0x3800
significand_mask(::Type{Float16}) = 0x03ff

for T in (Float16, Float32, Float64)
    @eval significand_bits(::Type{$T}) = $(trailing_ones(significand_mask(T)))
    @eval exponent_bits(::Type{$T}) = $(sizeof(T)*8 - significand_bits(T) - 1)
    @eval exponent_bias(::Type{$T}) = $(Int(exponent_one(T) >> significand_bits(T)))
    # maximum float exponent
    @eval exponent_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)) - exponent_bias(T))
    # maximum float exponent without bias
    @eval exponent_raw_max(::Type{$T}) = $(Int(exponent_mask(T) >> significand_bits(T)))
end

neededdigits(::Type{Float64}) = 309 + 17
neededdigits(::Type{Float32}) = 39 + 9 + 2
neededdigits(::Type{Float16}) = 9 + 5 + 9

function writeshortest(x::T,
        plus::Bool=false,
        space::Bool=false,
        hash::Bool=true,
        precision::Integer=-1,
        expchar::UInt8=UInt8('e'),
        padexp::Bool=false,
        decchar::UInt8=UInt8('.'),
        typed::Bool=false,
        compact::Bool=false) where {T <: Base.IEEEFloat}
    buf = Base.StringVector(neededdigits(T))
    pos = writeshortest(buf, 1, x, plus, space, hash, precision, expchar, padexp, decchar, typed, compact)
    return String(resize!(buf, pos - 1))
end

"""
    b, e10 = reduce_shortest(f[, maxsignif])

Reduce to shortest decimal representation where `abs(f) == b * 10^e10` and `b` is an
integer. If a `maxsignif` argument is provided, then `b < maxsignif`.
"""
@inline function reduce_shortest(f::T, maxsignif=nothing) where {T}
    U = uinttype(T)
    uf = reinterpret(U, f)
    m = uf & significand_mask(T)
    e = ((uf & exponent_mask(T)) >> significand_bits(T)) % Int

    ## Step 1
    #  mf * 2^ef == f
    mf = (one(U) << significand_bits(T)) | m
    ef = e - exponent_bias(T) - significand_bits(T)
    f_isinteger = mf & ((one(U) << -ef) - one(U)) == 0

    if ef > 0 || ef < -significand_bits(T) || !f_isinteger
        # fixup subnormals
        if e == 0
            ef = 1 - exponent_bias(T) - significand_bits(T)
            mf = m
        end

        ## Step 2
        #  u * 2^e2 == (f + prevfloat(f))/2
        #  v * 2^e2 == f
        #  w * 2^e2 == (f + nextfloat(f))/2
        e2 = ef - 2
        mf_iseven = iseven(mf) # trailing bit of significand is zero

        v = U(4) * mf
        w = v + U(2)
        u_shift_half = m == 0 && e > 1 # if first element of binade, other than first normal one
        u = v - U(2) + u_shift_half

        ## Step 3
        #  a == floor(u * 2^e2 / 10^e10), exact if a_allzero
        #  b == floor(v * 2^e2 / 10^e10), exact if b_allzero
        #  c == floor(w * 2^e2 / 10^e10)
        a_allzero = false
        b_allzero = false
        b_lastdigit = 0x00
        if e2 >= 0
            q = log10pow2(e2) - (T == Float64 ? (e2 > 3) : 0)
            e10 = q
            k = pow5_inv_bitcount(T) + pow5bits(q) - 1
            i = -e2 + q + k
            a, b, c = mulshiftinvsplit(T, u, v, w, q, i)
            if T == Float32 || T == Float16
                if q != 0 && div(c - 1, 10) <= div(a, 10)
                    l = pow5_inv_bitcount(T) + pow5bits(q - 1) - 1
                    mul = pow5invsplit_lookup(T, q-1)
                    b_lastdigit = (mulshift(v, mul, -e2 + q - 1 + l) % 10) % UInt8
                end
            end
            if q <= qinvbound(T)
                if ((v % UInt32) - 5 * div(v, 5)) == 0
                    b_allzero = pow5(v, q)
                elseif mf_iseven
                    a_allzero = pow5(u, q)
                else
                    c -= pow5(w, q)
                end
            end
        else
            q = log10pow5(-e2) - (T == Float64 ? (-e2 > 1) : 0)
            e10 = q + e2
            i = -e2 - q
            k = pow5bits(i) - pow5_bitcount(T)
            j = q - k
            a, b, c = mulshiftsplit(T, u, v, w, i, j)
            if T == Float32 || T == Float16
                if q != 0 && div(c - 1, 10) <= div(a, 10)
                    j = q - 1 - (pow5bits(i + 1) - pow5_bitcount(T))
                    mul = pow5split_lookup(T, i+1)
                    b_lastdigit = (mulshift(v, mul, j) % 10) % UInt8
                end
            end
            if q <= 1
                b_allzero = true
                if mf_iseven
                    a_allzero = !u_shift_half
                else
                    c -= 1
                end
            elseif q < qbound(T)
                b_allzero = pow2(v, q - (T != Float64))
            end
        end

        ## Step 4: reduction
        if a_allzero || b_allzero
            # a) slow loop
            while true
                c_div10 = div(c, 10)
                a_div10 = div(a, 10)
                if c_div10 <= a_div10
                    break
                end
                a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                a_allzero &= a_mod10 == 0
                b_allzero &= b_lastdigit == 0
                b_lastdigit = b_mod10 % UInt8
                b = b_div10
                c = c_div10
                a = a_div10
                e10 += 1
            end
            if a_allzero
                while true
                    a_div10 = div(a, 10)
                    a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
                    if a_mod10 != 0 && (maxsignif === nothing || b < maxsignif)
                        break
                    end
                    c_div10 = div(c, 10)
                    b_div10 = div(b, 10)
                    b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                    b_allzero &= b_lastdigit == 0
                    b_lastdigit = b_mod10 % UInt8
                    b = b_div10
                    c = c_div10
                    a = a_div10
                    e10 += 1
                end
            end
            if b_allzero && b_lastdigit == 5 && iseven(b)
                b_lastdigit = UInt8(4)
            end
            roundup = (b == a && (!mf_iseven || !a_allzero)) || b_lastdigit >= 5
        else
            # b) specialized for common case (99% Float64, 96% Float32)
            roundup = b_lastdigit >= 5
            c_div100 = div(c, 100)
            a_div100 = div(a, 100)
            if c_div100 > a_div100
                b_div100 = div(b, 100)
                b_mod100 = (b % UInt32) - UInt32(100) * (b_div100 % UInt32)
                roundup = b_mod100 >= 50
                b = b_div100
                c = c_div100
                a = a_div100
                e10 += 2
            end
            while true
                c_div10 = div(c, 10)
                a_div10 = div(a, 10)
                if c_div10 <= a_div10
                    break
                end
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                roundup = b_mod10 >= 5
                b = b_div10
                c = c_div10
                a = a_div10
                e10 += 1
            end
            roundup = (b == a || roundup)
        end
        if maxsignif !== nothing && b > maxsignif
            # reduce to max significant digits
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b <= maxsignif
                    break
                end
                b = b_div10
                roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
                b_allzero &= b_mod10 == 0
                e10 += 1
            end
            b = b + roundup

            # remove trailing zeros
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b_mod10 != 0
                    break
                end
                b = b_div10
                e10 += 1
            end
        else
            b = b + roundup
        end
    else
        # c) specialized f an integer < 2^53
        b = mf >> -ef
        e10 = 0

        if maxsignif !== nothing && b > maxsignif
            b_allzero = true
            # reduce to max significant digits
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b <= maxsignif
                    break
                end
                b = b_div10
                roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
                b_allzero &= b_mod10 == 0
                e10 += 1
            end
            b = b + roundup
        end
        while true
            b_div10 = div(b, 10)
            b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
            if b_mod10 != 0
                break
            end
            b = b_div10
            e10 += 1
        end
    end
    return b, e10
end

@inline function writeshortest(buf::Vector{UInt8}, pos, x::T,
                               plus=false, space=false, hash=true,
                               precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'),
                               typed=false, compact=false) where {T}
    @assert 0 < pos <= length(buf)
    neg = signbit(x)
    # special cases
    if x == 0
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('0')
        pos += 1
        if hash
            buf[pos] = decchar
            pos += 1
        end
        if precision == -1
            buf[pos] = UInt8('0')
            pos += 1
            if typed && x isa Float32
                buf[pos] = UInt8('f')
                buf[pos + 1] = UInt8('0')
                pos += 2
            end
            return pos
        end
        while precision > 1
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
        return pos
    elseif isnan(x)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        if typed
            if x isa Float32
                buf[pos + 3] = UInt8('3')
                buf[pos + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + 3] = UInt8('1')
                buf[pos + 4] = UInt8('6')
            end
        end
        return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    elseif !isfinite(x)
        if neg
            buf[pos] = UInt8('-')
        end
        buf[pos + neg] = UInt8('I')
        buf[pos + neg + 1] = UInt8('n')
        buf[pos + neg + 2] = UInt8('f')
        if typed
            if x isa Float32
                buf[pos + neg + 3] = UInt8('3')
                buf[pos + neg + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + neg + 3] = UInt8('1')
                buf[pos + neg + 4] = UInt8('6')
            end
        end
        return pos + neg + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    end

    output, nexp = reduce_shortest(x, compact ? 999_999 : nothing)

    if typed && x isa Float16
        buf[pos] = UInt8('F')
        buf[pos + 1] = UInt8('l')
        buf[pos + 2] = UInt8('o')
        buf[pos + 3] = UInt8('a')
        buf[pos + 4] = UInt8('t')
        buf[pos + 5] = UInt8('1')
        buf[pos + 6] = UInt8('6')
        buf[pos + 7] = UInt8('(')
        pos += 8
    end
    if neg
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end

    olength = decimallength(output)
    exp_form = true
    pt = nexp + olength
    if -4 < pt <= (precision == -1 ? (T == Float16 ? 3 : 6) : precision)
        exp_form = false
        if pt <= 0
            buf[pos] = UInt8('0')
            pos += 1
            buf[pos] = decchar
            pos += 1
            for _ = 1:abs(pt)
                buf[pos] = UInt8('0')
                pos += 1
            end
            # elseif pt >= olength
            # nothing to do at this point
            # else
            # nothing to do at this point
        end
    else
        pos += 1
    end
    i = 0
    ptr = pointer(buf)
    ptr2 = pointer(DIGIT_TABLE)
    if (output >> 32) != 0
        q = output ÷ 100000000
        output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
        output = q

        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        d = output2 % UInt32(10000)
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        d0 = (d % 100) << 1
        d1 = (d ÷ 100) << 1
        memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
        memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
        memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
        i += 8
    end
    output2 = output % UInt32
    while output2 >= 10000
        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
        i += 4
    end
    if output2 >= 100
        c = (output2 % UInt32(100)) << 1
        output2 = div(output2, UInt32(100))
        memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
        i += 2
    end
    if output2 >= 10
        c = output2 << 1
        buf[pos + 1] = DIGIT_TABLE[c + 2]
        buf[pos - exp_form] = DIGIT_TABLE[c + 1]
    else
        buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
    end

    if !exp_form
        if pt <= 0
            pos += olength
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        elseif pt >= olength
            pos += olength
            precision -= olength
            for _ = 1:nexp
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
            if hash
                buf[pos] = decchar
                pos += 1
                if precision < 0
                    buf[pos] = UInt8('0')
                    pos += 1
                end
                while precision > 0
                    buf[pos] = UInt8('0')
                    pos += 1
                    precision -= 1
                end
            end
        else
            pointoff = olength - abs(nexp)
            memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
            buf[pos + pointoff] = decchar
            pos += olength + 1
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
    else
        if olength > 1 || hash
            buf[pos] = decchar
            pos += olength
            precision -= olength
        end
        if hash && olength == 1
            buf[pos] = UInt8('0')
            pos += 1
        end
        while hash && precision > 0
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end

        buf[pos] = expchar
        pos += 1
        exp2 = nexp + olength - 1
        if exp2 < 0
            buf[pos] = UInt8('-')
            pos += 1
            exp2 = -exp2
        elseif padexp
            buf[pos] = UInt8('+')
            pos += 1
        end

        if exp2 >= 100
            c = exp2 % 10
            memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
            buf[pos + 2] = UInt8('0') + (c % UInt8)
            pos += 3
        elseif exp2 >= 10
            memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
            pos += 2
        else
            if padexp
                buf[pos] = UInt8('0')
                pos += 1
            end
            buf[pos] = UInt8('0') + (exp2 % UInt8)
            pos += 1
        end
    end
    if typed && x isa Float16
        buf[pos] = UInt8(')')
        pos += 1
    end

    return pos
end

memcpy(d, doff, s, soff, n) = (ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), d + doff - 1, s + soff - 1, n); nothing)
memmove(d, doff, s, soff, n) = (ccall(:memmove, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), d + doff - 1, s + soff - 1, n); nothing)

# Note: these are smaller than the values given in Figure 4 from the paper
# see https://github.com/ulfjack/ryu/issues/119
pow5_bitcount(::Type{Float16}) = 30
pow5_bitcount(::Type{Float32}) = 61
pow5_bitcount(::Type{Float64}) = 121

pow5_inv_bitcount(::Type{Float16}) = 30
pow5_inv_bitcount(::Type{Float32}) = 59
pow5_inv_bitcount(::Type{Float64}) = 122

qinvbound(::Type{Float16}) = 4
qinvbound(::Type{Float32}) = 9
qinvbound(::Type{Float64}) = 21

qbound(::Type{Float16}) = 15
qbound(::Type{Float32}) = 31
qbound(::Type{Float64}) = 63

log10pow2(e) = (e * 78913) >> 18
log10pow5(e) = (e * 732923) >> 20
pow5bits(e) = ((e * 1217359) >> 19) + 1

@inline function mulshift(m::U, mul, j) where {U<:Unsigned}
    W = widen(U)
    nbits = 8*sizeof(U)
    return ((((W(m) * (mul % U)) >> nbits) + W(m) * (mul >> nbits)) >> (j - nbits)) % promote_type(U,UInt32)
end

indexforexp(e) = div(e + 15, 16)
pow10bitsforindex(idx) = 16 * idx + 120
lengthforindex(idx) = div(((Int64(16 * idx) * 1292913986) >> 32) + 1 + 16 + 8, 9)

@inline function pow5(x, p)
    count = 0
    while true
        q = div(x, 5)
        r = x - 5 * q
        r != 0 && return count >= p
        x = q
        count += 1
    end
end

pow2(x, p) = (x & ((Int64(1) << p) - 1)) == 0

@inline function decimallength(v)
    v >= 10000000000000000 && return 17
    v >= 1000000000000000 && return 16
    v >= 100000000000000 && return 15
    v >= 10000000000000 && return 14
    v >= 1000000000000 && return 13
    v >= 100000000000 && return 12
    v >= 10000000000 && return 11
    v >= 1000000000 && return 10
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end
@inline function decimallength(v::UInt32)
    v >= 100000000 && return 9
    v >= 10000000 && return 8
    v >= 1000000 && return 7
    v >= 100000 && return 6
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end
@inline function decimallength(v::UInt16)
    v >= 10000 && return 5
    v >= 1000 && return 4
    v >= 100 && return 3
    v >= 10 && return 2
    return 1
end

@inline function mulshiftinvsplit(::Type{T}, mv, mp, mm, i, j) where {T}
    mul = pow5invsplit_lookup(T, i)
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

@inline function mulshiftsplit(::Type{T}, mv, mp, mm, i, j) where {T}
    mul = pow5split_lookup(T, i)
    vr = mulshift(mv, mul, j)
    vp = mulshift(mp, mul, j)
    vm = mulshift(mm, mul, j)
    return vr, vp, vm
end

function pow5invsplit(::Type{T}, i) where {T<:AbstractFloat}
    W = widen(uinttype(T))
    pow = big(5)^i
    inv = div(big(1) << (ndigits(pow, base=2) - 1 + pow5_inv_bitcount(T)), pow) + 1
    return W(inv)
end

function pow5invsplit_lookup end
for T in (Float64, Float32, Float16)
    e2_max = exponent_max(T) - precision(T) - 2
    i_max = log10pow2(e2_max)
    table = [pow5invsplit(T, i) for i = 0:i_max]
    @eval pow5invsplit_lookup(::Type{$T}, i) = @inbounds($table[i+1])
end

function pow5split(::Type{T}, i) where {T<:AbstractFloat}
    W = widen(uinttype(T))
    pow = big(5)^i
    return W(pow >> (ndigits(pow, base=2) - pow5_bitcount(T)))
end

function pow5split_lookup end
for T in (Float64, Float32, Float16)
    e2_min = 1 - exponent_bias(T) - significand_bits(T) - 2
    i_max = 1 - e2_min - log10pow5(-e2_min)
    table = [pow5split(T, i) for i = 0:i_max]
    @eval pow5split_lookup(::Type{$T}, i) = @inbounds($table[i+1])
end

const DIGIT_TABLE = UInt8[
  '0','0','0','1','0','2','0','3','0','4','0','5','0','6','0','7','0','8','0','9',
  '1','0','1','1','1','2','1','3','1','4','1','5','1','6','1','7','1','8','1','9',
  '2','0','2','1','2','2','2','3','2','4','2','5','2','6','2','7','2','8','2','9',
  '3','0','3','1','3','2','3','3','3','4','3','5','3','6','3','7','3','8','3','9',
  '4','0','4','1','4','2','4','3','4','4','4','5','4','6','4','7','4','8','4','9',
  '5','0','5','1','5','2','5','3','5','4','5','5','5','6','5','7','5','8','5','9',
  '6','0','6','1','6','2','6','3','6','4','6','5','6','6','6','7','6','8','6','9',
  '7','0','7','1','7','2','7','3','7','4','7','5','7','6','7','7','7','8','7','9',
  '8','0','8','1','8','2','8','3','8','4','8','5','8','6','8','7','8','8','8','9',
  '9','0','9','1','9','2','9','3','9','4','9','5','9','6','9','7','9','8','9','9'
]


text:
using BabyNeuralODE
using Test

data = read_halfmoons("../data/halfmoons.dat", 2000);

# no middle layer
train(data, 0);

# Dense middle layer
train(data, 1);

# ODE middle layer
train(data, 2);

# using Plots
#
# t1 = []
# t2 = []
# for _ in 1:5
#     push!(t1, Essais.train(data, 1))
#     plot()
#     plot!(t1, label="Dense", color=:red)
#     plot!(t2, label="ODE", color=:blue)
#     savefig("plot.pdf")
#
#     push!(t2, Essais.train(data, 2))
#     plot()
#     plot!(t1, label="Dense", color=:red)
#     plot!(t2, label="ODE", color=:blue)
#     savefig("plot.pdf")
# end

ext: jl


text:
# Block-based deque

#######################################
#
#  DequeBlock
#
#######################################

type DequeBlock{T}
    data::Vector{T}  # only data[front:back] is valid
    capa::Int
    front::Int
    back::Int
    prev::DequeBlock{T}  # ref to previous block
    next::DequeBlock{T}  # ref to next block

    function DequeBlock(capa::Int, front::Int)
        data = Array(T, capa)
        blk = new(data, capa, front, front-1)
        blk.prev = blk
        blk.next = blk
        blk
    end
end

# block at the rear of the chain, elements towards the front
rear_deque_block{T}(ty::Type{T}, n::Int) = DequeBlock{T}(n, 1)

# block at the head of the train, elements towards the back
head_deque_block{T}(ty::Type{T}, n::Int) = DequeBlock{T}(n, n+1)

capacity(blk::DequeBlock) = blk.capa
length(blk::DequeBlock) = blk.back - blk.front + 1
isempty(blk::DequeBlock) = blk.back < blk.front

# reset the block to empty, and position

function reset!{T}(blk::DequeBlock{T}, front::Int)
    blk.front = front
    blk.back = front - 1
    blk.prev = blk
    blk.next = blk
end

function show(io::IO, blk::DequeBlock)  # avoids recursion into prev and next
    x = blk.data[blk.front:blk.back]
    print(io, "$(typeof(blk))(capa = $(blk.capa), front = $(blk.front), back = $(blk.back)): $x")
end

#######################################
#
#  Deque
#
#######################################

const DEFAULT_DEQUEUE_BLOCKSIZE = 1024

type Deque{T}
    nblocks::Int
    blksize::Int
    len::Int
    head::DequeBlock{T}
    rear::DequeBlock{T}

    function Deque(blksize::Int)
        head = rear = rear_deque_block(T, blksize)
        new(1, blksize, 0, head, rear)
    end

    Deque() = Deque{T}(DEFAULT_DEQUEUE_BLOCKSIZE)
end

deque{T}(::Type{T}) = Deque{T}()

isempty(q::Deque) = q.len == 0
length(q::Deque) = q.len
num_blocks(q::Deque) = q.nblocks

function front(q::Deque)
    isempty(q) && error("Attempted to front at an empty deque.")
    blk = q.head
    blk.data[blk.front]
end

function back(q::Deque)
    isempty(q) && error("Attempted to back at an empty deque.")
    blk = q.rear
    blk.data[blk.back]
end

# Iteration

immutable DequeIterator{T}
    is_done::Bool
    cblock::DequeBlock{T}  # current block
    i::Int
end

start{T}(q::Deque{T}) = DequeIterator{T}(isempty(q), q.head, q.head.front)

function next{T}(q::Deque{T}, s::DequeIterator{T})
    cb = s.cblock
    i::Int = s.i
    x::T = cb.data[i]

    is_done = false

    i += 1
    if i > cb.back
        cb_next = cb.next
        if is(cb, cb_next)
            is_done = true
        else
            cb = cb_next
            i = 1
        end
    end

    (x, DequeIterator{T}(is_done, cb, i))
end

done{T}(q::Deque{T}, s::DequeIterator{T}) = s.is_done

function Base.collect{T}(q::Deque{T})
    r = T[]
    for x::T in q
        push!(r, x)
    end
    return r
end

# Showing

function show(io::IO, q::Deque)
    print(io, "Deque [$(collect(q))]")
end

function dump(io::IO, q::Deque)
    println(io, "Deque (length = $(q.len), nblocks = $(q.nblocks))")
    cb::DequeBlock = q.head
    i = 1
    while true
        print(io, "block $i [$(cb.front):$(cb.back)] ==> ")
        for j = cb.front : cb.back
            print(io, cb.data[j])
            print(io, ' ')
        end
        println(io)

        cb_next::DequeBlock = cb.next
        if !is(cb, cb_next)
            cb = cb_next
            i += 1
        else
            break
        end
    end
end

# Manipulation

function empty!{T}(q::Deque{T})
    # release all blocks except the head
    if q.nblocks > 1
        cb::DequeBlock{T} = q.rear
        while cb != q.head
            empty!(cb.data)
            cb = cb.prev
        end
    end

    # clean the head block (but retain the block itself)
    reset!(q.head, 1)

    # reset queue fields
    q.nblocks = 1
    q.len = 0
    q.rear = q.head
    q
end

function push!{T}(q::Deque{T}, x)  # push back
    rear = q.rear

    if isempty(rear)
        rear.front = 1
        rear.back = 0
    end

    if rear.back < rear.capa
        @inbounds rear.data[rear.back += 1] = convert(T, x)
    else
        new_rear = rear_deque_block(T, q.blksize)
        new_rear.back = 1
        new_rear.data[1] = convert(T, x)
        new_rear.prev = rear
        q.rear = rear.next = new_rear
        q.nblocks += 1
    end
    q.len += 1
    q
end

function unshift!{T}(q::Deque{T}, x)   # push front
    head = q.head

    if isempty(head)
        n = head.capa
        head.front = n + 1
        head.back = n
    end

    if head.front > 1
        @inbounds head.data[head.front -= 1] = convert(T, x)
    else
        n::Int = q.blksize
        new_head = head_deque_block(T, n)
        new_head.front = n
        new_head.data[n] = convert(T, x)
        new_head.next = head
        q.head = head.prev = new_head
        q.nblocks += 1
    end
    q.len += 1
    q
end

function pop!{T}(q::Deque{T})   # pop back
    isempty(q) && error("Attempted to pop from an empty deque.")
    rear = q.rear
    @assert rear.back >= rear.front

    @inbounds x = rear.data[rear.back]
    rear.back -= 1
    if rear.back < rear.front
        if q.nblocks > 1
            # release and detach the rear block
            empty!(rear.data)
            q.rear = rear.prev::DequeBlock{T}
            q.rear.next = q.rear
            q.nblocks -= 1
        end
    end
    q.len -= 1
    x
end

function shift!{T}(q::Deque{T})  # pop front
    isempty(q) && error("Attempted to pop from an empty deque.")
    head = q.head
    @assert head.back >= head.front

    @inbounds x = head.data[head.front]
    head.front += 1
    if head.back < head.front
        if q.nblocks > 1
            # release and detach the head block
            empty!(head.data)
            q.head = head.next::DequeBlock{T}
            q.head.prev = q.head
            q.nblocks -= 1
        end
    end
    q.len -= 1
    x
end

max_stars_repo_name: ambujone/DIDebug

max_issues_repo_name: ambujone/DIDebug

max_forks_repo_name: ambujone/DIDebug


text:
using HybridArrays
using StaticArrays
using Test

struct ScalarTest end
Base.:(+)(x::Number, y::ScalarTest) = x
Broadcast.broadcastable(x::ScalarTest) = Ref(x)

@testset "broadcasting" begin
    Ai = rand(Int, 2, 3, 4)
    Bi = HybridArray{Tuple{2,3,StaticArrays.Dynamic()}, Int, 3}(Ai)

    Af = rand(Float64, 2, 3, 4)
    Bf = HybridArray{Tuple{2,3,StaticArrays.Dynamic()}, Float64, 3}(Af)

    Bi[1,2,:] .= [110, 111, 112, 113]
    @test Bi[1,2,:] == @SVector [110, 111, 112, 113]

    @testset "Scalar Broadcast" begin
        @test Bf == @inferred(Bf .+ ScalarTest())
        @test Bf .+ 1 == @inferred(Bf .+ Ref(1))
    end

    @testset "AbstractArray-of-HybridArray with scalar math" begin
        v = [Bf]
        @test @inferred(v .* 1.0)::typeof(v) == v
    end

    @testset "2x2 HybridMatrix with HybridVector" begin
        m = HybridMatrix{2,StaticArrays.Dynamic()}([1 2; 3 4])
        v = HybridVector{2}([1, 4])
        @test @inferred(broadcast(+, m, v)) == @SMatrix [2 3; 7 8]
        @test @inferred(m .+ v) == @SMatrix [2 3; 7 8]
        @test @inferred(v .+ m) == @SMatrix [2 3; 7 8]
        @test @inferred(m .* v) == @SMatrix [1 2; 12 16]
        @test @inferred(v .* m) == @SMatrix [1 2; 12 16]
        @test @inferred(m ./ v) == @SMatrix [1 2; 3/4 1]
        @test @inferred(v ./ m) == @SMatrix [1 1/2; 4/3 1]
        @test @inferred(m .- v) == @SMatrix [0 1; -1 0]
        @test @inferred(v .- m) == @SMatrix [0 -1; 1 0]
        @test @inferred(m .^ v) == @SMatrix [1 2; 81 256]
        @test @inferred(v .^ m) == @SMatrix [1 1; 64 256]
        # StaticArrays Issue #546
        @test @inferred(m ./ (v .* v')) == @SMatrix [1.0 0.5; 0.75 0.25]
        testinf(m, v) = m ./ (v .* v')
        @test @inferred(testinf(m, v)) == @SMatrix [1.0 0.5; 0.75 0.25]
    end

    @testset "2x2 HybridMatrix with 1x2 HybridMatrix" begin
        # StaticArrays Issues #197, #242: broadcast between SArray and row-like SMatrix
        m1 = HybridMatrix{2,StaticArrays.Dynamic()}([1 2; 3 4])
        m2 = HybridMatrix{1,StaticArrays.Dynamic()}([1 4])
        @test @inferred(broadcast(+, m1, m2)) == @SMatrix [2 6; 4 8]
        @test @inferred(m1 .+ m2) == @SMatrix [2 6; 4 8]
        @test @inferred(m2 .+ m1) == @SMatrix [2 6; 4 8]
        @test @inferred(m1 .* m2) == @SMatrix [1 8; 3 16]
        @test @inferred(m2 .* m1) == @SMatrix [1 8; 3 16]
        @test @inferred(m1 ./ m2) == @SMatrix [1 1/2; 3 1]
        @test @inferred(m2 ./ m1) == @SMatrix [1 2; 1/3 1]
        @test @inferred(m1 .- m2) == @SMatrix [0 -2; 2 0]
        @test @inferred(m2 .- m1) == @SMatrix [0 2; -2 0]
        @test @inferred(m1 .^ m2) == @SMatrix [1 16; 3 256]
    end

    @testset "1x2 HybridMatrix with SVector" begin
        # StaticArrays Issues #197, #242: broadcast between SVector and row-like SVector
        m = HybridMatrix{1,StaticArrays.Dynamic()}([1 2])
        v = SVector(1, 4)
        @test @inferred(broadcast(+, m, v)) == @SMatrix [2 3; 5 6]
        @test @inferred(m .+ v) == @SMatrix [2 3; 5 6]
        @test @inferred(v .+ m) == @SMatrix [2 3; 5 6]
        @test @inferred(m .* v) == @SMatrix [1 2; 4 8]
        @test @inferred(v .* m) == @SMatrix [1 2; 4 8]
        @test @inferred(m ./ v) == @SMatrix [1 2; 1/4 1/2]
        @test @inferred(v ./ m) == @SMatrix [1 1/2; 4 2]
        @test @inferred(m .- v) == @SMatrix [0 1; -3 -2]
        @test @inferred(v .- m) == @SMatrix [0 -1; 3 2]
        @test @inferred(m .^ v) == @SMatrix [1 2; 1 16]
        @test @inferred(v .^ m) == @SMatrix [1 1; 4 16]
    end

    @testset "HybridMatrix with HybridMatrix" begin
        m1 = HybridMatrix{2,StaticArrays.Dynamic()}([1 2; 3 4])
        m2 = HybridMatrix{2,StaticArrays.Dynamic()}([1 3; 4 5])
        @test @inferred(broadcast(+, m1, m2)) == @SMatrix [2 5; 7 9]
        @test @inferred(m1 .+ m2) == @SMatrix [2 5; 7 9]
        @test @inferred(m2 .+ m1) == @SMatrix [2 5; 7 9]
        @test @inferred(m1 .* m2) == @SMatrix [1 6; 12 20]
        @test @inferred(m2 .* m1) == @SMatrix [1 6; 12 20]
        # StaticArrays Issue #199: broadcast with empty SArray
        @test @inferred(HybridVector{1}([1]) .+ HybridVector{0,Int}([])) === SVector{0,Union{}}()
        @test_broken @inferred(HybridVector{0,Int}([]) .+ SVector(1)) === SVector{0,Union{}}()
        # StaticArrays Issue #200: broadcast with Adjoint
        @test @inferred(m1 .+ m2') == @SMatrix [2 6; 6 9]
        @test @inferred(m1 .+ transpose(m2)) == @SMatrix [2 6; 6 9]
        # StaticArrays Issue 382: infinite recursion in Base.Broadcast.broadcast_indices with Adjoint
        @test @inferred(HybridVector{2}([1,1])' .+ [1, 1]) == [2 2; 2 2]
        @test @inferred(transpose(HybridVector{2}([1,1])) .+ [1, 1]) == [2 2; 2 2]
        @test @inferred(HybridVector{StaticArrays.Dynamic()}([1,1])' .+ [1, 1]) == [2 2; 2 2]
        @test @inferred(transpose(HybridVector{StaticArrays.Dynamic()}([1,1])) .+ [1, 1]) == [2 2; 2 2]
    end

    @testset "HybridMatrix with Scalar" begin
        m = HybridMatrix{2,StaticArrays.Dynamic()}([1 2; 3 4])
        @test @inferred(broadcast(+, m, 2)) == @SMatrix [3 4; 5 6]
        @test @inferred(m .+ 2) == @SMatrix [3 4; 5 6]
        @test @inferred(2 .+ m) == @SMatrix [3 4; 5 6]
        @test @inferred(m .* 2) == @SMatrix [2 4; 6 8]
        @test @inferred(2 .* m) == @SMatrix [2 4; 6 8]
        @test @inferred(m ./ 2) == @SMatrix [1/2 1; 3/2 2]
        @test @inferred(2 ./ m) == @SMatrix [2 1; 2/3 1/2]
        @test @inferred(m .- 2) == @SMatrix [-1 0; 1 2]
        @test @inferred(2 .- m) == @SMatrix [1 0; -1 -2]
        @test @inferred(m .^ 2) == @SMatrix [1 4; 9 16]
        @test @inferred(2 .^ m) == @SMatrix [2 4; 8 16]
    end
    @testset "Empty arrays" begin
        @test @inferred(1.0 .+ HybridMatrix{2,0,Float64}(zeros(2,0))) == HybridMatrix{2,0,Float64}(zeros(2,0))
        @test @inferred(1.0 .+ HybridMatrix{0,2,Float64}(zeros(0,2))) == HybridMatrix{0,2,Float64}(zeros(0,2))
        @test @inferred(1.0 .+ HybridArray{Tuple{2,StaticArrays.Dynamic(),0},Float64}(zeros(2,3,0))) == HybridArray{Tuple{2,StaticArrays.Dynamic(),0},Float64}(zeros(2,3,0))
        @test @inferred(HybridVector{0,Float64}(zeros(0)) .+ HybridMatrix{0,2,Float64}(zeros(0,2))) == HybridMatrix{0,2,Float64}(zeros(0,2))
        m = HybridMatrix{0,2,Float64}(zeros(0,2))
        @test @inferred(broadcast!(+, m, m, HybridVector{0,Float64}(zeros(0)))) == HybridMatrix{0,2,Float64}(zeros(0,2))
    end

    @testset "Mutating broadcast!" begin
        # No setindex! error
        A = HybridMatrix{2,StaticArrays.Dynamic()}([1 0; 0 1])
        @test @inferred(broadcast!(+, A, A, SVector(1, 4))) == @MMatrix [2 1; 4 5]
        A = HybridMatrix{2,StaticArrays.Dynamic()}([1 0; 0 1])
        @test @inferred(broadcast!(+, A, A, @SMatrix([1  4]))) == @MMatrix [2 4; 1 5]
        A = HybridMatrix{1,StaticArrays.Dynamic()}([1 0])
        @test_throws DimensionMismatch broadcast!(+, A, A, SVector(1, 4))
        A = HybridMatrix{1,StaticArrays.Dynamic()}([1 0])
        @test @inferred(broadcast!(+, A, A, @SMatrix([1 4]))) == @MMatrix [2 4]
        A = HybridMatrix{1,StaticArrays.Dynamic()}([1 0])
        @test @inferred(broadcast!(+, A, A, 2)) == @MMatrix [3 2]
    end

    @testset "broadcast! with mixtures of SArray and Array" begin
        A = HybridVector{StaticArrays.Dynamic()}([0, 0])
        @test @inferred(broadcast!(+, A, [1,2])) == [1,2]
    end

    @testset "eltype after broadcast" begin
        # test cases StaticArrays issue #198
        let a = HybridVector{4, Number}(Number[2, 2.0, 4//2, 2+0im])
            @test eltype(a .+ 2) == Number
            @test eltype(a .- 2) == Number
            @test eltype(a * 2) == Number
            @test eltype(a / 2) == Number
        end
        let a = HybridVector{3, Real}(Real[2, 2.0, 4//2])
            @test eltype(a .+ 2) == Real
            @test eltype(a .- 2) == Real
            @test eltype(a * 2) == Real
            @test eltype(a / 2) == Real
        end
        let a = HybridVector{3, Real}(Real[2, 2.0, 4//2])
            @test eltype(a .+ 2.0) == Float64
            @test eltype(a .- 2.0) == Float64
            @test eltype(a * 2.0) == Float64
            @test eltype(a / 2.0) == Float64
        end
        let a = broadcast(Float32, HybridVector{3}([3, 4, 5]))
            @test eltype(a) == Float32
        end
    end

    @testset "broadcast general scalars" begin
        # StaticArrays Issue #239 - broadcast with non-numeric element types
        @eval @enum Axis aX aY aZ
        @test (HybridVector{3}([aX,aY,aZ]) .== Ref(aX)) == HybridVector{3}([true,false,false])
        mv = HybridVector{3}([aX,aY,aZ])
        @test broadcast!(identity, mv, Ref(aX)) == HybridVector{3}([aX,aX,aX])
        @test mv == HybridVector{3}([aX,aX,aX])
    end

    @testset "broadcast! with Array destination" begin
        # Issue #385
        a = zeros(3, 3)
        b = HybridMatrix{3,StaticArrays.Dynamic()}([1 2 3; 4 5 6; 7 8 9])
        a .= b
        @test a == b

        c = HybridVector{3}([1, 2, 3])
        a .= c
        @test a == [1 1 1; 2 2 2; 3 3 3]

        d = HybridVector{4}([1, 2, 3, 4])
        @test_throws DimensionMismatch a .= d
    end
end


text:
# This file is a part of ParallelProcessingTools.jl, licensed under the MIT License (MIT).

function _workpart_hi(n_items::T, n_partitions::T, i::T) where {T<:Integer}
    @assert n_items >= 0
    @assert n_partitions > 0
    @assert i >= 0 && i <= n_partitions

    a = div(n_items, n_partitions)
    b = rem(n_items, n_partitions)
    c = div(b * i, n_partitions)
    (a * i + c)::T
end

function _workpart_scheme(range::Base.OneTo{<:Integer}, n_partitions::Integer, i::Integer)
    n_items = last(range)
    n_items_T, n_partitions_T, i_T = promote(n_items, n_partitions, i)

    if !(0 < i_T <= n_partitions_T)
        one(i_T):zero(i_T)
    else
        (_workpart_hi(n_items_T, n_partitions_T, i_T - 1) + 1):_workpart_hi(n_items_T, n_partitions_T, i_T)
    end
end

_workpart_scheme(range::StepRange, n::Integer, i::Integer) =
    step(range) * (_workpart_scheme(Base.OneTo(length(range)), n, i) .- 1) .+ first(range)

_workpart_scheme(range::AbstractUnitRange, n::Integer, i::Integer) =
    _workpart_scheme(Base.OneTo(length(range)), n, i) .+ (first(range) - 1)

_workpart_scheme(A::AbstractArray, n::Integer, i::Integer) =
    view(A, _workpart_scheme(eachindex(A), n, i))

_is_sorted_and_unique(A::AbstractRange) = true

function _is_sorted_and_unique(A::AbstractVector)
    if isempty(A)
        return true
    elseif !issorted(A)
        return false
    else
        idxs = LinearIndices(A)
        prev = A[first(idxs)]
        @inbounds for i in (first(idxs)+1):last(idxs)
            current = A[i]
            if current == prev
                return false
            else
                prev = current
            end
        end
        return true
    end
end

"""
    workpart(data::AbstractArray, workersel::AbstractVector{W}, current_worker::W) where {W}

Get the part of `data` that the execution unit `current_worker` is
responsible for. Implies a partition of `data` across the workers listed in
`workersel`.

For generic `data` arrays, `workpart` will return a view. If `data` is a
`Range` (e.g. indices to be processed), a sub-range will be returned.

Type `W` will typically be `Int` and `workersel` will usually be a range/array
of thread/process IDs.

Note: `workersel` is required to be sorted in ascending order and to contain
no duplicate entries.

Examples:

```julia
using Distributed, Base.Threads
A = rand(100)
# ...
sub_A = workpart(A, workers(), myid())
# ...
idxs = workpart(eachindex(sub_A), allthreads(), threadid())
for i in idxs
    # ...
end
```
"""
function workpart end
export workpart

function workpart(data::AbstractArray, workersel::AbstractVector{W}, current_worker::W) where {W}
    _is_sorted_and_unique(workersel) || throw(ArgumentError("List of selected workers must be sorted and contain no duplicates"))

    worker_idxs = eachindex(workersel)
    n = length(worker_idxs)
    i = searchsortedfirst(workersel, current_worker)
    checkindex(Bool, worker_idxs, i) && workersel[i] == current_worker || throw(ArgumentError("Worker $current_worker is not in specified list of selected workers"))

    _workpart_scheme(data, n, i)
end

function workpart(data::AbstractArray, workersel::Integer, current_worker::Integer)
    current_worker == workersel || throw(ArgumentError("Worker $current_worker does not equal specified worker"))
    data
end

@deprecate(
    workpartition(A::AbstractArray, n::Integer, i::Integer),
    workpart(A, 1::n, i)
)

@deprecate(
    threadpartition(A::AbstractArray, n_threads::Integer = length(allthreads()), i::Integer = threadid()),
    workpart(A, 1::n_threads, i)
)


text:
using Flux
using StatsBase, Statistics, IterTools
using Flux: Data.DataLoader
using Printf, EvalMetrics

@doc """
    residual_block(f::Int64, ks::Int64)

Create residual block of `f` filters and `ks` kernel-size.
""" ->
function residual_block(f::Int64, ks::Int64)
    sigmoid_tanh = Chain(x->σ.(x) .* tanh.(x))
    return Chain(
        Dropout(0.1),
        Conv((ks,), f => f, identity; stride=ks),
        sigmoid_tanh,
        Conv((1,), f => f, identity; stride=1)
    )
end

@doc """
    build_wavenet(features::Int64, in_dim::Int64;
                     filters::Int64=5, kernel_size::Int64=2, out_dim::Int64)

Build the wavenet-nn for samples which consists of `in_dim` sub-samples and each sub sample has got `features` dimension.
The `filters` and `kernel_size` are numbers of filters and kernel size for each convolution, `out_dim` is number of target classes.
Returns the model.
""" ->
function build_wavenet(features::Int64, in_dim::Int64; filters::Int64=12, kernel_size::Int64=2, out_dim::Int64=1)

    parts = []

    stack_conv1d = Conv((1,), features=>filters, identity; stride=1)
    append!(parts, [stack_conv1d])

    res_blocks = []
    res_connetcions_mean = []
    res_connetcions_max = []
    res_lc_coefs = []
    res_blocks_cnt::Int64 = in_dim

    #=
    Create residual blocks.
    =#
    while res_blocks_cnt >= kernel_size
        res_block = residual_block(filters, kernel_size)
        append!(res_blocks, [res_block])
        res_blocks_cnt = Int(floor(res_blocks_cnt / kernel_size))

        a = 1.0
        b = 1.0
        append!(res_connetcions_mean, [Chain(AdaptiveMeanPool((res_blocks_cnt,)), x-> x.*a)])
        append!(res_connetcions_max, [Chain(AdaptiveMaxPool((res_blocks_cnt,)), x-> x.*b)])
        append!(res_lc_coefs, [[a, b]])
    end

    #=
    Create skip blocks.
    =#
    skip_blocks = []
    for i=1:length(res_blocks)
        append!(skip_blocks, [Chain(Conv((1,), filters => 1, identity; stride=1),
                                    Flux.flatten,
                                    Dense(Int(floor(in_dim/2^i)), 1, identity))])
    end
    σ_ = Chain(x->x.+0.5, x->σ.(x))
    final = Chain(Dense(length(skip_blocks)+1, out_dim, σ_))
    append!(parts, [final])

    #=
    Connect blocks.
    =#
    function model(x)
        x::Array{Float32, 3} = stack_conv1d(x)
        res_blocks_cnt = length(res_blocks)
        skip_connections = Array{Float32, 2}(zeros(1, size(x)[3]))
        for i::Int64=1:res_blocks_cnt
            x_res_block = res_blocks[i]
            x_max = res_connetcions_max[i]
            x_mean = res_connetcions_mean[i]
            x = Parallel(+, x_res_block, x_max, x_mean)(x, x, x)
            skip_connections = vcat(skip_connections, skip_blocks[i](x))
        end
        return final(skip_connections)
    end
    return model, Flux.params(
        parts,
        res_blocks,
        res_lc_coefs,
        skip_blocks
    )
end

# Compute MSE loss function for 'network'.
loss_mse(x, y, network) = Flux.mse(network(x)', y)
# Accuracy of 'network' for given treshold 't', data 'x' and targets 'y'.
acc_t(x, y, t, network) = sum((network(x) .> t) .== y')/length(y)

@doc """
    train!(ps, network, trn::Tuple{Array{Float32, 3}, Vector{Int64}},
        val::Tuple{Array{Float32, 3}, Vector{Int64}};
        loss=loss_mse, acc=acc_t, threshold::Float64=0.5, iter::Int64=100)

Train parametres 'ps' of the 'network' on train 'trn' dataset.
""" ->
function train!(ps, network, trn::Tuple{Array{Float32, 3}, Vector{Int64}},
                    val::Tuple{Array{Float32, 3}, Vector{Int64}};
                    loss=loss_mse, acc=acc_t, threshold::Float64=0.5, iter::Int64=100)
    X, Y = trn[1], trn[2]
    x, y = val[1], val[2]
    for i=1:iter
        gs = gradient(() -> loss(X, Y, network), ps)
        Flux.Optimise.update!(ADAM(), ps, gs)
        if i%10 == 0
            @printf "iter: %d/%d, trn_loss: %.4f, trn_acc: %.2f; val_loss: %.4f, val_acc: %.2f\n" i iter loss(X, Y, network) acc(X, Y, threshold, network) loss(x, y, network) acc(x, y, threshold, network)
            @printf "trn_auc %.2f, val_auc %.2f\n\n" binary_eval_report(Y, Vector{Float64}((network(X)[1,:])))["au_roccurve"] binary_eval_report(y, Vector{Float64}((network(x)[1,:])))["au_roccurve"]

        end
    end
    ps
end

@doc """
    print_conmat(ŷ::Vector{Float64}, y::Vector{Int64}, absolut_values::Bool)

Print confusion matrix for predicted probabilities 'ŷ' (with treshold 0.5) and targets 'y'.
""" ->
function print_conmat(ŷ::Vector{Float64}, y::Vector{Int64}, absolut_values::Bool)
    cm = ConfusionMatrix(y, ŷ, 0.5)
    if absolut_values
        @printf "
        +---------------------+------------------+------------------+
        |          .          | Actual positives | Actual negatives |
        +---------------------+------------------+------------------+
        | Predicted positives | %12d     | %12d     |
        | Predicted negatives | %12d     | %12d     |
        +---------------------+------------------+------------------+
        \n" cm.tp cm.fp cm.fn cm.tn
    else
        total = cm.p+cm.n
        @printf "
        +---------------------+------------------+------------------+
        |          .          | Actual positives | Actual negatives |
        +---------------------+------------------+------------------+
        | Predicted positives |       %.2f       |       %.2f       |
        | Predicted negatives |       %.2f       |       %.2f       |
        +---------------------+------------------+------------------+
        \n" cm.tp/total cm.fp/total cm.fn/total cm.tn/total
    end
end

max_stars_repo_name: kozvojtex/MMWavenet

max_issues_repo_name: kozvojtex/MMWavenet

max_forks_repo_name: kozvojtex/MMWavenet


text:
module RentCalculator

using Parameters

@with_kw struct Param
    buy_price = 800_000.0 # home purchase price (dollars)
    stay = 5.0 # time in residence (years)
    mortgage_rate = 3.5 # mortgage rate (percent)
    mortgage_down = 20.0 # mortgage down payment (percent of purchase price)
    mortgage_length = 30.0 # mortgage length (years)
    rent_monthly = 3_600.0 # monthly rent (dollars)
    rent_deposit = 3_500.0 # rental deposit (dollars)
    rent_broker = 0.0 # rent broker's fee (dollars)
    rent_insurance = 0.5 # rental insurance rate (percent of purchase price)
    price_rate = 4.0 # annual home price appreciation rate (percent)
    rent_rate = 2.0 # annual rent rate increase (percent)
    market_rate = 7.0 # annual nominal return on investment (percent)
    inflation_rate = 2.0 # annual inflation rate (percent)
    tax_property = 1.0 # annual property tax rate (percent)
    tax_marginal = 33.3 # marginal tax rate (percent)
    tax_investment = 24.3 # tax rate on investments (percent)
    tax_capitalgains = 100_000.0 # captial gains deduction (dollars)
    buy_closing = 3.0 # cost of buying house (percent of purchase price)
    buy_selling = 5.0 # cost of selling house (percent of selling price)
    buy_renovation = 0.5 # annual renovation costs (percent of purchase price)
    buy_maintenance = 0.5 # annual maintenance costs (percent of purchase price)
    buy_insurance = 0.2 # homeowners insurance (percent of purchase price)
    buy_utilities = 200.0 # additional monthly utilities compared to renting (dollars)
    buy_common = 400.0 # common fee (dollars)
    buy_deductible = 20.0 # percent of common fee that is deductible (percent)
end

"""
    interest(p, horizon)

Calculate yearly interest on a mortgage loan with details specified in Param `p` for time `horizon`.

# Examples
```julia-repl
julia> interest(Param(), 30)
[long output]
```
"""
function interest(p::Param, horizon)
    interest_cumulative = zeros(horizon)
    interest_gains = zeros(horizon)
    return true
end

end

#


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2019
#
# This is a Julia package for discrete event simulation
#

"supertype for clocks in `DiscreteEvents.jl`"
abstract type AbstractClock end

"supertype for events"
abstract type AbstractEvent end

# abstract types for clock state machines
abstract type ClockState end
abstract type ClockEvent end

"""
    Timing

Enumeration type for scheduling events and timed conditions:

- `at`: schedule an event at a given time,
- `after`: schedule an event a given time after current time,
- `every`: schedule an event every given time from now on,
- `before`: a timed condition is true before a given time,
- `until`: delay until t.
"""
@enum Timing at after every before until

"""
    Action

An action is either  a `Function` or an `Expr` or a `Tuple` of them. It can
be scheduled in an event for later execution.

!!! warning "Evaluating expressions is slow"
    Expr should be avoided in time critical parts of applications. You will get a one
    time warning if you use them. They can be replaced easily by `fun`s or function
    closures. They are evaluated at global scope in Module `Main` only. Other modules using
    `DiscreteEvents.jl` cannot use Expr in events and have to use functions.
"""
const Action = Union{Function, Expr, Tuple}

"""
    DiscreteEvent{T<:Action} <: AbstractEvent

A discrete event is an `Action` to be executed at an event time.

# Arguments, fields
- `ex::T`: a function or an expression or a tuple of them,
- `t::Float64`: event time,
- `Δt::Float64`: repeat rate with for repeating events.
"""
struct DiscreteEvent{T<:Action} <: AbstractEvent
    ex::T
    t::Float64
    Δt::Float64
end
DiscreteEvent(ex::T, t::Real, Δt::Real) where {T<:Action} =
    DiscreteEvent(ex, float(t), float(Δt))

"""
    DiscreteCond{S<:Action, T<:Action} <: AbstractEvent

A condition to be evaluated repeatedly with expressions or functions
to be executed if conditions are met.

# Arguments, fields
- `cond::S`: a conditional function or an expression or a tuple of them
    (conditions must evaluate to `Bool`),
- `ex::T`: a function or an expression or a tuple of them to be executed
    if conditions are met,
"""
struct DiscreteCond{S<:Action, T<:Action} <: AbstractEvent
    cond::S
    ex::T
end

"""
    Sample{T<:Action} <: AbstractEvent

Sampling actions are executed at sampling time.

# Arguments, fields
- `ex<:Action`: an [`Action`](@ref) to be executed at sample time.
"""
struct Sample{T<:Action} <: AbstractEvent
    ex::T
end

"""
    ClockException(ev::ClockEvent, value=nothing)

Define a ClockException, which can be thrown to processes.

# Arguments, fields
- `ev::ClockEvent`: delivers an event to the interrupted task
- `value=nothing`: deliver some other value
"""
struct ClockException <: Exception
  ev::ClockEvent
  value::Any
  ClockException(ev::ClockEvent, value=nothing) = new(ev, value)
end

"""
    Prc(id, f, arg...; kw...)

Prepare a function to run as a process (asynchronous task) in a simulation.

# Arguments, fields
- `id`: some unique identification for registration,
- `f::Function`: a function `f(clk, arg...; kw...)`, must take `clk` (a [`Clock`](@ref))
    as its first argument,
- `arg...`: further arguments to `f`
- `kw...`: keyword arguments to `f`

# Fields
- `task::Union{Task,Nothing}`: a task structure,
- `clk::Union{AbstractClock,Nothing}`: clock where the process is registered,
- `state::ClockState`: process state,

!!! note
    A function started as a Prc runs in a loop. It has to give back control
    by e.g. doing a `take!(input)` or by calling [`delay!`](@ref) or [`wait!`](@ref),
    which will `yield` it. Otherwise it will starve everything else!
"""
mutable struct Prc
    id::Any
    task::Union{Task,Nothing}
    clk::Union{AbstractClock,Nothing}
    f::Function
    arg::Tuple
    kw::Base.Iterators.Pairs

    Prc( id, f::Function, arg...; kw...) = new(id, nothing, nothing, f, arg, kw)
end

"""
    Schedule()

A Schedule contains events, conditional events and sampling functions to be
executed or evaluated on the clock's time line.

# Fields
- `events::PriorityQueue{DiscreteEvent,Float64}`: scheduled events,
- `cevents::Array{DiscreteCond,1}`: conditional events to evaluate at each tick,
- `samples::Array{Sample,1}`: sampling expressions to evaluate at each tick,
"""
mutable struct Schedule
    events::PriorityQueue{DiscreteEvent,Float64}
    cevents::Array{DiscreteCond,1}
    samples::Array{Sample,1}

    Schedule() = new(PriorityQueue{DiscreteEvent,Float64}(), DiscreteCond[], Sample[])
end

"""
    ClockChannel

Provide a channel to an active clock or a real time clock.

# Fields
- `ref::Ref{Task}`: a pointer to an active clock,
- `ch::Channel`: a communication channel to an active clock,
- `id::Int`: the thread id of the active clock.
- `done::Bool`: flag indicating if the active clock has completed its cycle.
"""
mutable struct ClockChannel{T <: ClockEvent}
    ref::Ref{Task}
    forth::Channel{T}
    back::Channel{T}
    thread::Int
    done::Bool
    load::Int
end

"""
```
Clock(Δt::T=0; t0::U=0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}
```
Create a new simulation clock.

# Arguments
- `Δt::T=0`: time increment. If no Δt is given, the simulation doesn't tick,
    but jumps from event to event. Δt can be set later with `sample_time!`.
- `t0::U=0`: start time for simulation
- `unit::FreeUnits=NoUnits`: clock time unit. Units can be set explicitely by
    setting e.g. `unit=minute` or implicitly by giving Δt as a time or else setting
    t0 to a time, e.g. `t0=60s`.

# Fields
- `id::Int`: clock ident number, 0: master clock, ≥ 1: parallel clock,
- `state::ClockState`: clock state,
- `time::Float64`: clock time,
- `unit::FreeUnits`: time unit,
- `end_time::Float64`: end time for simulation,
- `Δt::Float64`: sampling time, timestep between ticks,
- `ac::Vector{ClockChannel}`: [`channels`](@ref ClockChannel) to active clocks on parallel threads,
- `sc::Schedule`: the clock [`Schedule`](@ref) (events, cond events and sampling),
- `processes::Dict{Any, Prc}`: registered `Prc`es,
- `tn::Float64`: next timestep,
- `tev::Float64`: next event time,
- `evcount::Int`: event counter,
- `scount::Int`: sample counter

# Examples

```jldoctest
julia> using DiscreteEvents, Unitful

julia> import Unitful: s, minute, hr

julia> c = Clock()                 # create a unitless clock (standard)
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c1 = Clock(1s, unit=minute)  # create a clock with unit [minute]
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 minute, Δt=0.01667 minute, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c2 = Clock(1s)               # create a clock with implicit unit [s]
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 s, Δt=1.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c3 = Clock(t0=60s)           # another clock with implicit unit [s]
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=60.0 s, Δt=0.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0

julia> c4 = Clock(1s, t0=1hr)       # here Δt's unit [s] takes precedence
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=3600.0 s, Δt=1.0 s, prc:0
  scheduled ev:0, cev:0, sampl:0
```
"""
mutable struct Clock <: AbstractClock
    id::Int
    state::ClockState
    time::Float64
    unit::FreeUnits
    Δt::Float64
    ac::Vector{ClockChannel}
    sc::Schedule
    processes::Dict{Any, Prc}
    tn::Float64
    tev::Float64
    end_time::Float64
    evcount::Int
    scount::Int

    function Clock(Δt::T=0;
                   t0::U=0, unit::FreeUnits=NoUnits) where {T<:Number,U<:Number}
        if 1unit isa Time
            Δt = isa(Δt, Time) ? uconvert(unit, Δt).val : Δt
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
        elseif Δt isa Time
            unit = Unitful.unit(Δt)
            t0 = isa(t0, Time) ? uconvert(unit, t0).val : t0
            Δt = Δt.val
        elseif t0 isa Time
            unit = Unitful.unit(t0)
            t0 = t0.val
        else
            nothing
        end
        new(0, Undefined(), t0, unit, Δt, ClockChannel[], Schedule(),
            Dict{Any, Prc}(), t0 + Δt, t0, t0, 0, 0)
    end
end

"""
    ActiveClock(clock, master, cmd, ans)

A thread specific clock which can be operated via a channel.

# Fields
- `clock::Clock`: the thread specific clock,
- `master::Ref{Clock}`: a pointer to the master clock (on thread 1),
- `cmd::Channel{ClockEvent}`: the command channel from master,
- `ans::Channel{ClockEvent}`: the response channel to master.
- `id::Int`: the id in master's ac array,
- `thread::Int`: the thread, the active clock runs on.

!!! note
You should not setup an `ActiveClock` explicitly. Rather this is done
implicitly by [`fork!`](@ref)ing a [`Clock`](@ref) to other available threads
or directly with [`PClock`](@ref).
It then can be accessed via [`pclock`](@ref) as in the following example.

!!! note
Directly accessing the `clock`-field of parallel `ActiveClock`s in simulations is
possible but not recommended since it breaks parallel operation. The right way is
to pass `event!`s to the `ActiveClock`-variable. The communication then happens
over the channel to the `ActiveClock` as it should be.

# Example
```jldoctest; filter = r"^.*[0-9]+ ac.*"
julia> using DiscreteEvents

julia> clk = Clock()
Clock 0, thrd 1 (+ 0 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> fork!(clk)

julia> clk    #  ⬇ here you got 3 parallel active clocks
Clock 0, thrd 1 (+ 3 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.0 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> clk = PClock()
Clock 0, thrd 1 (+ 3 ac): state=DiscreteEvents.Undefined(), t=0.0 , Δt=0.01 , prc:0
  scheduled ev:0, cev:0, sampl:0

julia> ac1 = pclock(clk, 1)    # get access to the 1st active clock (on thread 2)
Active clock 1 on thrd 2: state=DiscreteEvents.Idle(), t=0.0 , Δt=0.01 , prc:0
   scheduled ev:0, cev:0, sampl:0
```
"""
mutable struct ActiveClock{T <: ClockEvent} <: AbstractClock
    clock::Clock
    master::Ref{Clock}
    forth::Channel{T}
    back ::Channel{T}
    id::Int
    thread::Int
end

"""
    RTClock{T <: ClockEvent} <: AbstractClock

Real time clocks use system time and are controlled over channels. They are independent
from each other and other clocks and run asynchronously as tasks on arbitrary threads.
Multiple real time clocks can be setup with arbitrary frequencies to do different things.

# Fields
- `clock::Clock`:
- `cmd::Channel{T}`:
- `back::Channel{T}`:
- `id::Int`:
- `thread::Int`:
- `t0::Float64`:
"""
mutable struct RTClock{T <: ClockEvent} <: AbstractClock
    clock::Clock
    cmd::Channel{T}
    back::Channel{T}
    id::Int
    thread::Int
    t0::Float64
end


text:
module Musculoskeletal

include("model.jl")

end # module


text:
# clones a table
# If the flush option is not enabled, then any data the source table currently has in memory will not exist in the clone.
# A cloned table copies the configuration of the source table.
# However the permissions of the source table are not copied to the clone.
# After a clone is created, only the user that created the clone can read and write to it.
function table_clone(session::AccumuloSession, tablename::AbstractString, newtablename::AbstractString; flush::Bool=true, set_properties::Dict=Dict(), exclude_properties::SET=())
    setprops = Dict{UTF8String,UTF8String}([utf8(string(n))=>utf8(string(v)) for (n,v) in set_properties]...)
    exclprops = Set(UTF8String[utf8(string(x)) for x in exclude_properties])
    cloneTable(client(session), handle(session), utf8(tablename), utf8(newtablename), flush, setprops, exclprops)
end

# create, delete, rename tables
#
# When a table is created, by default it is configured to use the VersioningIterator and keep one version. A table can be created without the VersioningIterator.
# Accumulo can be configured to return the top k versions, or versions later than a given date. The default is to return the one most recent version.
# The version policy can be changed by changing the VersioningIterator options with the `config` api.
#
# Using logical timestamps ensures that timestamps set by Accumulo always move forward. This helps avoid problems caused by TabletServers that have different time settings.
# The default is to use milliseconds as timestamp.
function table_create(session::AccumuloSession, tablename::AbstractString; versioning::Bool=true, logical_timestamp::Bool=false)
    time_type = logical_timestamp ? TimeType.LOGICAL : TimeType.MILLIS
    createTable(client(session), handle(session), utf8(tablename), versioning, time_type)
end
table_delete(session::AccumuloSession, tablename::AbstractString) = deleteTable(client(session), handle(session), utf8(tablename))
table_rename(session::AccumuloSession, tablename::AbstractString, new_name::AbstractString) = renameTable(client(session), handle(session), utf8(tablename), utf8(new_name))

# prints how much space, in bytes, is used by files referenced by a table.  When multiple tables are specified it prints how much space, in bytes, is used by files shared between tables, if any.
function table_du(session::AccumuloSession, tables::AbstractString...)
    tableset = Set(UTF8String[utf8(string(x)) for x in tables])
    usage = getDiskUsage(client(session), handle(session), tableset)
    (usage[1]).usage
end

table_exists(session::AccumuloSession, table::AbstractString) = tableExists(client(session), handle(session), utf8(table))

# get or set table specific properties
table_config(session::AccumuloSession, tablename::AbstractString) = getTableProperties(client(session), handle(session), utf8(tablename))
table_config!(session::AccumuloSession, tablename::AbstractString, property::AbstractString, value::AbstractString) = setTableProperty(client(session), handle(session), utf8(tablename), utf8(property), utf8(value))

# limit versions in a table
function table_versions!(session::AccumuloSession, tablename::AbstractString, nversions::Integer)
    (nversions > 0) || error("at least 1 version required")
    viter = iter("vers", "org.apache.accumulo.core.iterators.user.VersioningIterator", 10, Dict("maxVersions"=>nversions))
    add_iter(session, tablename, viter, (IteratorScope.MINC, IteratorScope.MAJC, IteratorScope.SCAN))
end

# import/export tables
#
# Accumulo supports exporting tables for the purpose of copying tables to another cluster.
# Exporting and importing tables preserves the tables configuration, splits, and logical time.
# Tables are exported and then copied via the hadoop distcp command.
# To export a table, it must be offline and stay offline while discp runs. The reason it needs to stay offline is to prevent files from being deleted.
# A table can be cloned and the clone taken offline inorder to avoid losing access to the table.
table_export(session::AccumuloSession, tablename::AbstractString, export_dir::AbstractString) = exportTable(client(session), handle(session), utf8(tablename), utf8(export_dir))
table_import(session::AccumuloSession, tablename::AbstractString, import_dir::AbstractString) = importTable(client(session), handle(session), utf8(tablename), utf8(import_dir))

# take tables online/offline
table_offline(session::AccumuloSession, tablename::AbstractString; wait::Bool=true) = offlineTable(client(session), handle(session), utf8(tablename), wait)
table_online(session::AccumuloSession, tablename::AbstractString; wait::Bool=true) = onlineTable(client(session), handle(session), utf8(tablename), wait)

# get a list of all existing tables
tables(session::AccumuloSession) = listTables(client(session), handle(session))

#=
namespaces not supported yet
# displays a list of all existing namespaces
function namespaces()
end

# creates a new namespace
function createnamespace()
end

# deletes a namespace
function deletenamespace()
end

# renames a namespace
function renamenamespace()
end
=#


text:
module Arduino

using SerialPorts

"""
    SerialPorts.Arduino.reset

Reset an Arduino board.
"""
function Base.reset(s::SerialPort)
    setDTR(s, false)
    sleep(1)
    setDTR(s, true)
    sleep(1)
    nothing
end

end #module


text:
function valuetype end
function makevalue end

struct Environment
    contacts::Vector{Contact}

    function Environment(mechanism, urdf, obstacles)
        elements = visual_elements(mechanism, URDFVisuals(urdf; tag="collision"))
        point_elements = filter(e -> e.geometry isa GeometryTypes.HyperSphere && GeometryTypes.radius(e.geometry) == 0, elements)
        points = map(point_elements) do element
            p = element.transform(SVector(GeometryTypes.origin(element.geometry)))
            Point3D(element.frame, p)
        end
        contacts = vec(map(Base.Iterators.product(points, obstacles)) do p
            point, obstacle = p
            body = body_fixed_frame_to_body(mechanism, point.frame)
            Contact(mechanism, body, point, obstacle)
        end)
        # removes collision with itself
        contacts = filter(c -> c.point.frame != c.obstacle.normal.frame, contacts)
        new(contacts)
    end
end

struct EnvironmentJacobian{T}
    contact_jacobians::Vector{ContactJacobian{T}}

    function EnvironmentJacobian{T}(env::Environment) where T
        contact_jacobians = [ContactJacobian{T}(contact) for contact in env.contacts]

        new{T}(contact_jacobians)
    end
end

struct EnvironmentJacobianCache <: RigidBodyDynamics.AbstractTypeDict
    env::Environment
    keys::Vector{Tuple{UInt64, Int}}
    values::Vector{EnvironmentJacobian}
end

function EnvironmentJacobianCache(env::Environment)
    EnvironmentJacobianCache(env, [], [])
end

Base.show(io::IO, ::EnvironmentJacobianCache) = print(io, "EnvironmentJacobianCache{…}(…)")

@inline RigidBodyDynamics.valuetype(::Type{EnvironmentJacobianCache}, ::Type{T}) where {T} = EnvironmentJacobian{T}
@inline RigidBodyDynamics.makevalue(envj_c::EnvironmentJacobianCache, ::Type{T}) where {T} = EnvironmentJacobian{T}(envj_c.env)

function contact_jacobian!(envj::EnvironmentJacobian, state::MechanismState)
    # TODO: parallel
    for cj in envj.contact_jacobians
        contact_jacobian!(cj, state)
    end
end


text:
using DrWatson, PyCall
@quickactivate "GCN_HM_GRN-Integration"

pushfirst!(PyVector(pyimport("sys")."path"), "");
run_sim = pyimport("port").run_sim

allparams = Dict(
    :k => [2, 3],
    :epochs => [500, 750],
    :cl => ["E116"]
)

dicts = dict_list(allparams)

function makesim(d::Dict)
    @unpack k, epochs, cl = d
    best_auc = run_sim(cl, epochs, k)
    fulld = copy(d)
    fulld[:auc] = best_auc
    return fulld
end

for (i, d) in enumerate(dicts)
    f = makesim(d)
    print("Finished $(i)")
    @tagsave(datadir("port_no_sep", savename(d, "bson")), f)
end


text:
using Test
using GeometricIntegrators
using GeometricProblems.LotkaVolterra2d
using GeometricProblems.LotkaVolterra2d: Δt
using SrkMethodsForDegenerateLagrangianSystems

const tableaus = (
    tableaus_vprk_glrk(),
    tableaus_vprk_lobatto_IIIA_IIIB(),
    tableaus_vprk_lobatto_IIIB_IIIA(),
    tableaus_vprk_radau(),
    tableaus_srk_glrk(),
    tableaus_firk_glrk(),
)

const iode = lotka_volterra_2d_iode()
const nt = 10

for list in tableaus
    for run in list
        tab, file = run

        if length(run) ≥ 3
            integrator = run[3]
        else
            integrator = Integrator
        end

        @test_nowarn begin
            sol = Solution(iode, Δt, nt)
            int = integrator(iode, tab, Δt)
            integrate!(int, sol)
        end
    end
end


text:
module Day15

export day15

function day15(input::Vector{Int}, upto::Int)
    sequence = Dict{Int, Tuple{Int, Int}}()
    for (i,k) in enumerate(input)
        sequence[k] = (i,i)
    end
    last = input[end]
    for i in length(input)+1:upto
        (prev,prevprev) = sequence[last]
        last = prev - prevprev
        (a,b) = get(sequence, last, (i,i))
        sequence[last] = (i,a)
    end
    last
end

end


text:
@testsuite "linear algebra" AT->begin
    @testset "adjoint and transpose" begin
        @test compare(adjoint, AT, rand(Float32, 32, 32))
        @test compare(adjoint!, AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
        @test compare(transpose, AT, rand(Float32, 32, 32))
        @test compare(transpose!, AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
        @test compare((x,y)->copyto!(x, adjoint(y)), AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
        @test compare((x,y)->copyto!(x, transpose(y)), AT, rand(Float32, 32, 32), rand(Float32, 32, 32))
        @test compare(transpose!, AT, Array{Float32}(undef, 32, 32), rand(Float32, 32, 32))
        @test compare(transpose!, AT, Array{Float32}(undef, 128, 32), rand(Float32, 32, 128))
    end

    @testset "copytri!" begin
        @testset for uplo in ('U', 'L')
            @test compare(x -> LinearAlgebra.copytri!(x, uplo), AT, rand(Float32, 128, 128))
        end
    end

    @testset "copyto! for triangular" begin
        for TR in (UpperTriangular, LowerTriangular)
            @test compare(transpose!, AT, Array{Float32}(undef, 128, 32), rand(Float32, 32, 128))

            cpu_a = Array{Float32}(undef, 128, 128)
            gpu_a = AT{Float32}(undef, 128, 128)
            gpu_b = AT{Float32}(undef, 128, 128)

            rand!(gpu_a)
            copyto!(cpu_a, TR(gpu_a))
            @test cpu_a == Array(collect(TR(gpu_a)))

            rand!(gpu_a)
            gpu_c = copyto!(gpu_b, TR(gpu_a))
            @test all(Array(gpu_b) .== TR(Array(gpu_a)))
            @test gpu_c isa AT
        end

        for TR in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular)
            gpu_a = AT{Float32}(undef, 128, 128) |> rand! |> TR
            gpu_b = AT{Float32}(undef, 128, 128) |> TR

            gpu_c = copyto!(gpu_b, gpu_a)
            @test all(Array(gpu_b) .== Array(gpu_a))
            @test all(Array(gpu_c) .== Array(gpu_a))
            @test gpu_c isa TR
        end
    end

    @testset "permutedims" begin
        @test compare(x -> permutedims(x, (2, 1)), AT, rand(Float32, 2, 3))
        @test compare(x -> permutedims(x, (2, 1, 3)), AT, rand(Float32, 4, 5, 6))
        @test compare(x -> permutedims(x, (3, 1, 2)), AT, rand(Float32, 4, 5, 6))
        @test compare(x -> permutedims(x, [2,1,4,3]), AT, randn(ComplexF64,3,4,5,1))
    end

    @testset "issymmetric/ishermitian" begin
        n = 128
        areal = randn(n,n)/2
        aimg  = randn(n,n)/2

        @testset for eltya in (Float32, Float64, ComplexF32, ComplexF64)
            a = convert(Matrix{eltya}, eltya <: Complex ? complex.(areal, aimg) : areal)
            asym = transpose(a) + a        # symmetric indefinite
            aherm = a' + a                 # Hermitian indefinite
            @test issymmetric(asym)
            @test ishermitian(aherm)
        end
    end
    @testset "Array + Diagonal" begin
        n = 128
        A = AT{Float32}(undef, n, n)
        d = AT{Float32}(undef, n)
        rand!(A)
        rand!(d)
        D = Diagonal(d)
        B = A + D
        @test collect(B) ≈ collect(A) + collect(D)
    end

    @testset "$f! with diagonal $d" for (f, f!) in ((triu, triu!), (tril, tril!)),
                                        d in -2:2
        A = randn(10, 10)
        @test f(A, d) == Array(f!(AT(A), d))
    end

    @testset "$T gemv y := $f(A) * x * a + y * b" for f in (identity, transpose, adjoint), T in supported_eltypes()
        y, A, x = rand(T, 4), rand(T, 4, 4), rand(T, 4)

        # workaround for https://github.com/JuliaLang/julia/issues/35163#issue-584248084
        T <: Integer && (y .%= T(10); A .%= T(10); x .%= T(10))

        @test compare(*, AT, f(A), x)
        @test compare(mul!, AT, y, f(A), x)
        @test compare(mul!, AT, y, f(A), x, Ref(T(4)), Ref(T(5)))
        @test typeof(AT(rand(3, 3)) * AT(rand(3))) <: AbstractVector

        if f !== identity
            @test compare(mul!, AT, rand(T, 2,2), rand(T, 2,1), f(rand(T, 2)))
        end
    end

    @testset "$T gemm C := $f(A) * $g(B) * a + C * b" for f in (identity, transpose, adjoint), g in (identity, transpose, adjoint), T in supported_eltypes()
        A, B, C = rand(T, 4, 4), rand(T, 4, 4), rand(T, 4, 4)

        # workaround for https://github.com/JuliaLang/julia/issues/35163#issue-584248084
        T <: Integer && (A .%= T(10); B .%= T(10); C .%= T(10))

        @test compare(*, AT, f(A), g(B))
        @test compare(mul!, AT, C, f(A), g(B))
        @test compare(mul!, AT, C, f(A), g(B), Ref(T(4)), Ref(T(5)))
        @test typeof(AT(rand(3, 3)) * AT(rand(3, 3))) <: AbstractMatrix
    end

    @testset "lmul! and rmul!" for (a,b) in [((3,4),(4,3)), ((3,), (1,3)), ((1,3), (3))], T in supported_eltypes()
        @test compare(rmul!, AT, rand(T, a), Ref(rand(T)))
        @test compare(lmul!, AT, Ref(rand(T)), rand(T, b))
    end
end


text:
mutable struct MyMessage
    field1::Int32
    field2::Float64
end

function encode(msg::MyMessage)
    buf = IOBuffer()
    write(buf, hton(msg.field1))
    write(buf, hton(msg.field2))
    buf.data
end

function decode(data, msg::Type{MyMessage})
    buf = IOBuffer(data)
    MyMessage(ntoh(read(buf, Int32)), ntoh(read(buf, Float64)))
end


text:
using RecipesBase
using ColorTypes

@recipe function f(w::T,cache::BasicILMCache; layers=true) where {T<:GridData}

  framestyle --> :frame
  xlims --> (-Inf,Inf)
  ylims --> (-Inf,Inf)
  aspect_ratio := 1
  legend --> :none
  grid --> false

  @series begin
    trim := 2
    w, cache.g
  end

  if layers
    if typeof(w) <: ScalarGridData
      @series begin
        cache
      end
    else
      if typeof(w) <: VectorGridData
        nsub = 2
      elseif typeof(w) <: TensorGridData
        nsub = 4
      end
      for i in 1:nsub
        @series begin
          subplot := i
          cache
        end
      end
    end
  end

end

@recipe function f(w1::T1,w2::T2,cache::BasicILMCache; layers=true) where {T1<:GridData,T2<:GridData}
    @series begin
      w1, cache.g
    end

    @series begin
      linestyle --> :dash
      w2, cache.g
    end

    if layers
      #pts = points(cache)
      @series begin
        cache
      end
    end

end

@recipe function f(w::T,cache::BasicILMCache; bodyid=1) where {T<:ScalarData}

  legend --> :none
  grid --> false
  @series begin
    view(arcs(cache),cache,bodyid), view(w,cache,bodyid)
  end
end

@recipe function f(w::T,cache::BasicILMCache; bodyid=1) where {T<:VectorData}

  @series begin
    w.u, cache
  end

  @series begin
    w.v, cache
  end
end

@recipe function f(cache::BasicILMCache)
  pts = points(cache)
  #x = [pts.u; pts.u[1]]
  #y = [pts.v; pts.v[1]]
  x = [pts.u;]
  y = [pts.v;]
  seriestype := :scatter
  markersize --> 1
  markershape --> :circle
  markercolor --> :black
  xlims --> (-Inf,Inf)
  ylims --> (-Inf,Inf)
  aspect_ratio := 1
  legend --> :none
  grid --> false
  #linecolor --> :black
  #fillrange --> 0
  #fillalpha --> 0
  #fillcolor --> :black
  x := x
  y := y
  ()
end

@recipe function f(sys::ILMSystem)
    sys.base_cache
end

@recipe function f(w::T,sys::ILMSystem) where T<:Union{GridData,ScalarData,VectorData}
    w, sys.base_cache
end

@recipe function f(w1::T1,w2::T2,sys::ILMSystem) where {T1<:GridData,T2<:GridData}
    w1, w2, sys.base_cache
end


text:
"""
`SimpleTools` is a collection of useful functions that don't have an
obvious home, but are useful for my `SimpleWorld` modules.
"""
module SimpleTools

include("dict_compose.jl")
include("make_continuous.jl")
include("flush_print.jl")

end # end of module


text:
using CxxWrap
using Test


text:
"""
    function find_solution(board::Matrix{Int})

Find all the solutions that can be obtained in the least number of moves for the
given `board`. The used algorithm avoids cycles and returns an empty vector of
solutions if no solution exists for `board`.
"""
function find_solution(board::Matrix{Int})
    solutions = Vector{Vector{Move}}()
    current_possible_moves = get_possible_moves(board)
    for move in current_possible_moves
        push!(solutions, Vector{Move}([move]))
    end
    solution_found = false
    valid_solutions = Vector{Vector{Move}}()
    begin
        red_moves = _robot_possible_moves(board, get_red_position(board))
        if CartesianIndex(3,3) in [red_moves[i][2] for i in 1:length(red_moves)]
            push!(valid_solutions, [(get_red_position(board) => CartesianIndex(3,3))])
            solution_found = true
        end
    end

    while !solution_found && !isempty(solutions)
        nb_solutions = length(solutions)
        for _ in 1:nb_solutions
            current_solution = solutions[1]
            new_board = copy(board)
            list_of_boards = Vector{Matrix{Int}}()
            push!(list_of_boards, copy(new_board))
            for move in current_solution
                move_robot!(new_board, move)
                push!(list_of_boards, copy(new_board))
            end
            current_possible_moves = get_possible_moves(new_board)
            for move in current_possible_moves
                move_robot!(new_board, move)
                if new_board in list_of_boards
                    move_robot!(new_board, reverse(move))
                    continue
                end
                push!(solutions, copy(current_solution))
                push!(solutions[end], move)
                if is_finished(new_board)
                    solution_found = true
                    push!(valid_solutions, copy(solutions[end]))
                end
                move_robot!(new_board, reverse(move))
            end
            popfirst!(solutions)
        end
    end

    return valid_solutions
end

"""
    print_solution(_board::Matrix{Int}, solution::Vector{Move})

Print a solution in the format
1:   Color - Direction
2:   Color - Direction
etc.
"""
function print_solution(_board::Matrix{Int}, solution::Vector{Move})
    board = copy(_board)
    robotcolors = Dict{Int, String}([
        1 => "Red",
        2 => "Orange",
        3 => "Green",
        4 => "Purple",
        5 => "Yellow",
        6 => "Blue"
    ])
    i = 1
    for move in solution
        print("$i:\t$(robotcolors[board[move[1]]]) - ")
        if move[1][1] != move[2][1]
            if move[1][1] < move[2][1]
                println("Down")
            else
                println("Up")
            end
        else #move[1][2] != move[2][2]
            if move[1][2] < move[2][2]
                println("Right")
            else
                println("Left")
            end
        end
        i += 1
        move_robot!(board, move)
    end
    return nothing
end

"""
    print_solutions(board::Matrix{Int}, solutions::Vector{Vector{Move}})

Print all solutions stored in `solutions`.
"""
function print_solutions(board::Matrix{Int}, solutions::Vector{Vector{Move}})
    i = 1
    for solution in solutions
        println("Solution $i:")
        print_solution(board, solution)
        println()
        i += 1
    end
    return nothing
end

max_stars_repo_name: remi-garcia/LunarLockout

max_issues_repo_name: remi-garcia/LunarLockout

max_forks_repo_name: remi-garcia/LunarLockout


text:
abstract type AbstractDynamicalSystemResampling end

"""
    NoResampling

Indicates that no resampling should be performed.
"""
struct NoResampling <: AbstractDynamicalSystemResampling end

export AbstractDynamicalSystemResampling, NoResampling


text:
using Test
using ComparativeAutograder
using ComparativeAutograder: results_dict

@testset "Test results_dict" begin
    @testset "Test results_dict with add method." begin
        add = (x, y) -> x + y
        suite = TestSuite([
            TestCase(1, 2),
            TestCase(-1, 0),
        ])
        result_sub = runtestsuite(add, suite)
        result_sol = runtestsuite(+, suite)

        output = results_dict(suite, result_sub, result_sol)
        @test output["passed"] == true
        @test length(output["testCases"]) == 2

        @test output["testCases"][1]["passed"] == true
        @test output["testCases"][1]["passed"] == true
    end
end


text:
using StaticLint
using StaticLint: quoted, unquoted, should_mark_missing_getfield_ref, hasref,
    parentof, is_in_fexpr, errorof, haserror

using CSTParser
using CSTParser: EXPR, headof, isidentifier, isnonstdid, valof, iscall

using SymbolServer
using Printf
using Dates
using Sockets

# ..:: Variables and data structures ::..

mutable struct ErrorSpan
    beg_char::Integer
    end_char::Integer
end

SL = StaticLint
SS = SymbolServer
CSTP = SL.CSTParser
T_Error = Tuple{ErrorSpan,EXPR}

error = "error"
warn = "warn"
info = "info"

const LintCodeToErrorLevel = Dict{SL.LintCodes,String}(
    SL.IncorrectCallArgs => error,
    SL.IncorrectIterSpec => error,
    SL.NothingEquality => error,
    SL.NothingNotEq => error,
    SL.ConstIfCondition => info,
    SL.EqInIfConditional => error,
    SL.PointlessOR => info,
    SL.PointlessAND => info,
    SL.UnusedBinding => info,
    SL.InvalidTypeDeclaration => error,
    SL.UnusedTypeParameter => info,
    SL.IncludeLoop => warn,
    SL.MissingFile => error,
    SL.InvalidModuleName => error,
    SL.TypePiracy => warn,
    SL.UnusedFunctionArgument => warn,
    SL.CannotDeclareConst => error,
    SL.InvalidRedefofConst => error,
    SL.NotEqDef => warn,
    SL.KwDefaultMismatch => error,
    SL.InappropriateUseOfLiteral => warn,
    SL.ShouldBeInALoop => warn,
    SL.TypeDeclOnGlobalVariable => error,
    SL.UnsupportedConstLocalVariable => warn,
    SL.UnassignedKeywordArgument => error,
    SL.CannotDefineFuncAlreadyHasValue => error)

# ..:: Functions ::..

""" Make a string holding the current time. """
print_time() = @sprintf("[%s] ", Dates.format(now(), "E U yyyy HH:MM:SS"))

"""
    shorten_string(str)

Shorten the string to the last N of its characters.

# Arguments
- `str`: original string.

# Returns
The shortened string (just its tail).
"""
function shorten_string(str::AbstractString)::String
    max_length = 25
    return (length(str)<=max_length) ? str : "..."*last(str, max_length)
end # function

"""
    print_error(msg)

Print the error message.

# Arguments
- `msg`: the error message.
"""
function print_error(msg::String)::Nothing
    splitter = ".jl:"
    msg = split(msg, splitter)
    msg[1] = shorten_string(msg[1])
    msg = string(msg[1], splitter, msg[2])
    @info print_time()*msg
    return nothing
end # function

"""
    format_error(p, pos, desc, lvl)

Print out the error.

# Arguments
- `p`: the file path.
- `pos`: the (start,end) character positions.
- `desc`: the error description
- `lvl`: the error level (error? warning? info?)

# Returns
- `msg`: the error message in the agreed format.
"""
function format_error(p::String, pos::ErrorSpan, desc::String, lvl::String)::String
    msg = @sprintf "%s:%d:%d: %s: %s\n" p pos.beg_char pos.end_char lvl desc
    return msg
end

"""
    lint_collect_hints(x, server[, missingrefs, isquoted, errs, pos])

Collect linter output about the file.

This is basically copied from the StaticLint.collect_hints function, however I
updated it slightly to output the (start, end) character positions of the
error, and incorporated the feedback from
https://github.com/julia-vscode/CSTParser.jl/issues/238 on how to make the
character positions match correctly.

Note: the second output is used internally by the function. You don't need to
look at it when using the function to collect the array of errors (the first
output).
"""
function lint_collect_hints(x::EXPR,
                            server,
                            missingrefs=:all,
                            isquoted=false,
                            errs=T_Error[],
                            pos=1)::Tuple{Array{T_Error,1},Integer}
    if isquoted && unquoted(x)
        isquoted = false
    elseif quoted(x)
        isquoted = true
    end

    if headof(x) === :errortoken
        # collect parse errors
        push!(errs, (ErrorSpan(pos, pos+x.span), x))
    elseif !isquoted
        if missingrefs != :none &&
            isidentifier(x) &&
            !hasref(x) &&
            !(valof(x) == "var" &&
              parentof(x) isa EXPR &&
              isnonstdid(parentof(x))) &&
              !((valof(x) == "stdcall" ||
                 valof(x) == "cdecl" ||
                 valof(x) == "cdecl" ||
                 valof(x) == "fastcall" ||
                 valof(x) == "thiscall" ||
                 valof(x) == "llvmcall") &&
                is_in_fexpr(x, x -> (iscall(x) &&
                                     isidentifier(x.args[1]) &&
                                     valof(x.args[1]) == "ccall")))
            push!(errs, (ErrorSpan(pos, pos+x.span), x))
        elseif haserror(x) && errorof(x) isa SL.LintCodes
            # collect lint hints
            push!(errs, (ErrorSpan(pos, pos+x.span), x))
        end
    elseif (isquoted &&
            missingrefs == :all &&
            should_mark_missing_getfield_ref(x, server))
        push!(errs, (ErrorSpan(pos, pos+x.span), x))
    end

    for a in x
        if a.args === nothing
            lint_collect_hints(a, server, missingrefs, isquoted, errs, pos)
            pos += a.fullspan
        else
            _,pos = lint_collect_hints(a, server, missingrefs, isquoted, errs, pos)
        end
    end

    return errs, pos
end

"""
    convert_pos_byte_to_char!(src, offset)

Convert byte to character count.

Convert error start and end character locations from a byte measurement to a
strict character count measurement.

# Arguments
- `src`: the source code (as a string) which this error is for.
- `offset`: the error location in bytes. Modified in place.
"""
function convert_pos_byte_to_char!(src::String, offset::ErrorSpan)
    strlen_byte = lastindex(src)
    offset.beg_char = length(src, 1, min(offset.beg_char, strlen_byte))
    offset.end_char = length(src, 1, min(offset.end_char, strlen_byte))
end

"""
    get_error_level(x)

Convert error to an error level.

# Arguments
- `err`: the error object.

# Returns
- `lvl`: the error level.
"""
function get_error_level(x::EXPR)::String
    error_code = SL.errorof(x)
    lvl = error # Default error level
    if haskey(LintCodeToErrorLevel, error_code)
        lvl = LintCodeToErrorLevel[error_code]
    end
    return lvl
end

"""
    lint_file(rootfile, server, conn)

Static lint a file.

# Arguments
- `rootfile`: the file to be linted.
- `server`: the StaticLint server.
- `conn`: the TCP connection to send errors back over.
"""
function lint_file(rootfile::String,
                   server::SL.FileServer,
                   conn::TCPSocket)::Nothing
    @info print_time()*"Parsing file: "*shorten_string(rootfile)

    empty!(server.files)

    @info print_time()*"Loading file..."

    # Load and parse the file
    f = SL.loadfile(server, rootfile)

    @info print_time()*"Semanting pass..."

    # SL's main run- finding variables, scopes, etc.
    SL.semantic_pass(f)

    @info print_time()*"Running lint checks..."

    # Run lint checks.
    hints = Dict()
    slopts = SL.LintOptions(:)
    for (path, file) in server.files
        SL.check_all(file.cst, slopts, server)
        hints[path],_ = lint_collect_hints(file.cst, server)
    end

    # Send errors back to client
    files_src = Dict(f.path=>f.source)
    for (p, hs) in hints
        for (offset, x) in hs
            if !haskey(files_src, p)
                files_src[p] = read(p, String)
            end

            if p!=rootfile
                # Skip errors from other files
                continue
            end

            if (SL.haserror(x) &&
                SL.errorof(x) isa SL.LintCodes)
                description = SL.LintCodeDescriptions[SL.errorof(x)]
                level = get_error_level(x)
                convert_pos_byte_to_char!(files_src[p], offset)
                error_msg = format_error(p, offset, description, level)
                print_error(error_msg)
                write(conn, error_msg)
            else
                # Use of undeclared variable
                description = string("Missing reference for ", CSTP.valof(x))
                level = error
                convert_pos_byte_to_char!(files_src[p], offset)
                error_msg = format_error(p, offset, description, level)
                print_error(error_msg)
                write(conn, error_msg)
            end
        end
    end

    @info print_time()*"Finished"

    return nothing
end

# ..:: Start the server loop ::..

@info print_time()*"Initializing server"

# Julia configuration
depot = first(SS.Pkg.depots())
cache = joinpath(dirname(pathof(SS)), "..", "store")
env = dirname(SS.Pkg.Types.Context().env.project_file)

# Setup server
server = SL.FileServer()
ssi = SymbolServerInstance(depot, cache)
_, server.symbolserver = SS.getstore(ssi, env)
server.symbol_extends  = SS.collect_extended_methods(server.symbolserver)

@info print_time()*"Started server"

tcpserver = listen(1111) # Start the server
while true
    # Wait for client connection
    conn = accept(tcpserver)
    # Read the file to parse
    rootfile = readline(conn)
    # Stop the server?
    if rootfile == "stop"
        @info print_time()*"Received stop command, closing"
        close(conn)
        break
    end
    # Line the file
    try
        lint_file(rootfile, server, conn)
        write(conn, "<<end_lint>>")
    catch err
        bt = catch_backtrace()
        msg = sprint(showerror, err, bt)
        @error print_time()*msg
    end
    close(conn)
end

close(tcpserver)

@info print_time()*"Shutdown server"

max_stars_repo_name: dmalyuta/julia-staticlint

max_issues_repo_name: dmalyuta/julia-staticlint

max_forks_repo_name: dmalyuta/julia-staticlint


text:
@define HeaderCoefficientsRK begin
    name::Symbol
    o::Int
    s::Int
end

@define CoefficientsRK begin
    a::Matrix{T}
    b::Vector{T}
    c::Vector{T}
end

"Holds the coefficients of a Runge-Kutta method."
struct CoefficientsRK{T} <: AbstractCoefficients{T}
    @HeaderCoefficientsRK
    @CoefficientsRK
    â::Matrix{T}
    b̂::Vector{T}
    ĉ::Vector{T}

    function CoefficientsRK{T}(name,o,s,a,b,c) where {T}
        new(name,o,s,a,b,c,zero(a),zero(b),zero(c))
    end

    function CoefficientsRK{T}(name,o,s,a,b,c,â,b̂,ĉ) where {T}
        @assert T <: Real
        @assert isa(name, Symbol)
        @assert isa(o, Integer)
        @assert isa(s, Integer)
        @assert s > 0 "Number of stages must be > 0"
        @assert s==size(a,1)==size(a,2)==length(b)==length(c)
        @assert s==size(â,1)==size(â,2)==length(b̂)==length(ĉ)

        if !get_config(:tab_compensated_summation)
            â = zero(â)
            b̂ = zero(b̂)
            ĉ = zero(ĉ)
        end

        new(name,o,s,a,b,c,â,b̂,ĉ)
    end
end

function CoefficientsRK(T::Type, name::Symbol, order::Int, a::AbstractArray{CT,2}, b::AbstractArray{CT,1}, c::AbstractArray{CT,1}) where {CT}
    a̅ = Matrix{T}(a)
    b̅ = Vector{T}(b)
    c̅ = Vector{T}(c)

    if get_config(:tab_compensated_summation)
        â = Matrix{T}(a-Matrix{eltype(a)}(a̅))
        b̂ = Vector{T}(b-Vector{eltype(b)}(b̅))
        ĉ = Vector{T}(c-Vector{eltype(c)}(c̅))
    else
        â = zero(a̅)
        b̂ = zero(b̅)
        ĉ = zero(c̅)
    end

    CoefficientsRK{T}(name, order, length(c), a̅, b̅, c̅, â, b̂, ĉ)
end

function CoefficientsRK(name::Symbol, order::Int, a::AbstractArray{T,2}, b::AbstractArray{T,1}, c::AbstractArray{T,1}) where {T}
    CoefficientsRK{T}(name, order, length(c), a, b, c)
end

function CoefficientsRK(name::Symbol, order::Int,
        a::AbstractArray{T,2}, b::AbstractArray{T,1}, c::AbstractArray{T,1},
        â::AbstractArray{T,2}, b̂::AbstractArray{T,1}, ĉ::AbstractArray{T,1}) where {T}
    CoefficientsRK{T}(name, order, length(c), a, b, c, â, b̂, ĉ)
end

Base.hash(tab::CoefficientsRK, h::UInt) = hash(tab.o, hash(tab.s, hash(tab.a, hash(tab.b, hash(tab.c, hash(:CoefficientsRK, h))))))

Base.:(==)(tab1::CoefficientsRK, tab2::CoefficientsRK) = (tab1.o == tab2.o
                                                       && tab1.s == tab2.s
                                                       && tab1.a == tab2.a
                                                       && tab1.b == tab2.b
                                                       && tab1.c == tab2.c)

Base.isequal(tab1::CoefficientsRK{T1}, tab2::CoefficientsRK{T2}) where {T1,T2} = (tab1 == tab2 && T1 == T2 && typeof(tab1) == typeof(tab2))

"Print Runge-Kutta coefficients."
function Base.show(io::IO, tab::CoefficientsRK)
    print(io, "Runge-Kutta Coefficients ", tab.name, " with ", tab.s, " stages and order ", tab.o)
    print(io, "  a = ", tab.a)
    print(io, "  b = ", tab.b)
    print(io, "  c = ", tab.c)
end

function get_symplectic_conjugate_coefficients(a::Matrix{T}, b::Vector{T}, a̅::Matrix{T}) where {T}
    @assert size(a) == size(a̅)
    @assert length(b) == size(a,1) == size(a,2)

    for i in 1:size(a̅, 1)
        for j in 1:size(a̅, 2)
            a̅[i,j] = b[j] * ( one(T) - a[j,i] / b[i] )
        end
    end
    nothing
end

function get_symplectic_conjugate_coefficients(coeff::CoefficientsRK)
    a̅ = zero(coeff.a)
    ã = zero(coeff.a)

    get_symplectic_conjugate_coefficients(coeff.a, coeff.b, a̅)

    if all(x -> x≠0, coeff.b̂)
        get_symplectic_conjugate_coefficients(coeff.â, coeff.b̂, ã)
    end

    CoefficientsRK(coeff.name, coeff.s, a̅, coeff.b, coeff.c, ã, coeff.b̂, coeff.ĉ)
end

function symplecticize(coeff::CoefficientsRK; name=nothing, T=Float64)
    name == nothing ? Symbol(string(coeff.name)*"S") : nothing
    a̅ = zero(coeff.a)
    get_symplectic_conjugate_coefficients(coeff.a, coeff.b, a̅)
    CoefficientsRK(T, name, coeff.o, 0.5*(coeff.a + a̅), coeff.b, coeff.c)
end

function compute_symplecticity_error(coeff::CoefficientsRK)
    [coeff.b[i] * coeff.a[i,j] + coeff.b[j] * coeff.a[j,i] - coeff.b[i] * coeff.b[j] for i in 1:size(coeff.a,1), j in 1:size(coeff.a,2)]
end

function check_symplecticity(coeff::CoefficientsRK{T}) where {T}
    symplectic = falses(coeff.s, coeff.s)

    for i in 1:size(coeff.a, 1)
        for j in 1:size(coeff.a, 2)
            symplectic[i,j] = isapprox(coeff.b[i] * coeff.a[i,j] + coeff.b[j] * coeff.a[j,i], coeff.b[i] * coeff.b[j], atol=16*eps(T), rtol=16*eps(T))
        end
    end

    symplectic
end

function check_symmetry(coeff::CoefficientsRK{T}) where {T}
    symmetric = falses(coeff.s, coeff.s)

    for i in 1:size(coeff.a, 1)
        for j in 1:size(coeff.a, 2)
            symmetric[i,j] = isapprox(coeff.a[coeff.s+1-i, coeff.s+1-j] + coeff.a[i,j], coeff.b[j], atol=16*eps(T), rtol=16*eps(T))
        end
    end

    symmetric
end

function check_order_conditions_B(coeff::CoefficientsRK{T}, k) where {T}
    local res::T = 0

    for i in 1:coeff.s
        res += coeff.b[i] * coeff.c[i]^(k-1)
    end

    return isapprox(res, 1/k, atol=16*eps(T), rtol=16*eps(T))
end

function check_order_conditions_C(coeff::CoefficientsRK{T}, k) where {T}
    local order  = falses(coeff.s)
    local res::T

    for i in 1:size(coeff.a, 1)
        res = 0
        for j in 1:size(coeff.a, 2)
            res += coeff.a[i,j] * coeff.c[j]^(k-1)
        end
        order[i] = isapprox(res, coeff.c[i]^k/k, atol=16*eps(T), rtol=16*eps(T))
        # println(res - coeff.c[i]^k/k)
    end

    order
end

function check_order_conditions_D(coeff::CoefficientsRK{T}, k) where {T}
    local order  = falses(coeff.s)
    local res::T

    for i in 1:size(coeff.a, 1)
        res = 0
        for j in 1:size(coeff.a, 2)
            res += coeff.b[j] * coeff.c[j]^(k-1) * coeff.a[j,i]
        end
        order[i] = isapprox(res, coeff.b[i] * (1 - coeff.c[i]^k)/k, atol=16*eps(T), rtol=16*eps(T))
        # println(res - coeff.b[i] * (1 - coeff.c[i]^k)/k)
    end

    order
end


text:
"""
    cvresults =  cv_fpw(X, r_range; numFolds = 10, seed=1)

Performs cross validation for the radius parameter.
Candidate values for the radius parameter are specified in `r_range`.
Default number of folds is 10.
Seed controls the generation of the folds.
Returns the a matrix of dimensions (number of radii candidates)×(number of folds) of log-likelihoods evaluated on left out folds.

# Example
```julia-repl
julia> using Statistics, PyPlot
julia> X = spiraldata(300)
julia> r_range = LinRange(0.01, 2.0, 100)
julia> cvresults =  cv_fpw(X, r_range)
julia> r_perf = mean(cvresults, dims=2)
julia> best_index = argmax(r_perf)
julia> r_best = r_range[best_index]
julia> mix = fpw(X, r_best)
julia> x = rand(mix, 1000)'
julia> plot(X[:,1], X[:,2], "bo", label="data")
julia> plot(x[:,1], x[:,2], ".r", label="generated", alpha=0.7)
julia> legend()
```
"""
function cv_fpw(X, r_range; numFolds = 10, seed = 1, gamma = 1e-6, randrepeats = 7)

  Random.seed!(seed)

  # dimensions of data

  N, D = size(X)

  # partition datasets into disjoint sets

  sets = collect(Kfold(N, numFolds))

  score = zeros(length(r_range), numFolds)

  progressbar = Progress(length(r_range))

  for (r_index, r) in enumerate(r_range)

    for i = 1:numFolds

      # define indices for training set
      trainInd = sets[i]
      Xtrain   = @view X[trainInd,:]

      # define indices for testing set
      testInd = setdiff(collect(1:N), trainInd)
      Xtest   = @view X[testInd,:]

      for repeat=1:randrepeats

        # get Parzen window
        mixturemodel = fpw(Xtrain, r)

        # get log-likelihood on test set
        logpdfmix(x) = logpdf(mixturemodel, x)

        for nn=1:length(testInd)
          @inbounds score[r_index, i] += logpdfmix(@view Xtest[nn,:]) / randrepeats
        end

      end

    end

    ProgressMeter.next!(progressbar; showvalues = [(:r,r)])

  end

  return score

end