text:
wider(::Type{UInt64}) = UInt128
wider(::Type{UInt128}) = BigInt

inttype(::Type{UInt64}) = Int64
inttype(::Type{UInt128}) = Int128
inttype(::Type{T}) where {T} = BigInt

# Base.exponent_max(T) + Base.significand_bits(T) + 3 # "-0."
maxdigits(::Type{Float64}) = 1079
maxdigits(::Type{Float32}) = 154
maxdigits(::Type{Float16}) = 29
maxdigits(::Type{BigFloat}) = typemax(Int64)

# include a non-inlined version in case of widening (otherwise, all widened cases would fully inline)
@noinline _typeparser(::Type{T}, source, pos, len, b, code, options::Options{ignorerepeated, ignoreemptylines, Q, debug, S, D, DF}, ::Type{IntType}) where {T <: SupportedFloats, ignorerepeated, ignoreemptylines, Q, debug, S, D, DF, IntType} =
    typeparser(T, source, pos, len, b, code, options, IntType)

@inline function typeparser(::Type{T}, source, pos, len, b, code, options::Options{ignorerepeated, ignoreemptylines, Q, debug, S, D, DF}, ::Type{IntType}=UInt64) where {T <: SupportedFloats, ignorerepeated, ignoreemptylines, Q, debug, S, D, DF, IntType}
    startpos = pos
    origb = b
    x = zero(T)
    digits = zero(IntType)
    ndigits = 0
    if debug
        println("float parsing")
    end
    neg = b == UInt8('-')
    if neg || b == UInt8('+')
        pos += 1
        incr!(source)
    end
    if eof(source, pos, len)
        code |= INVALID | EOF
        @goto done
    end
    b = peekbyte(source, pos)
    if b == options.decimal
        @goto decimalcheck
    end
    b -= UInt8('0')
    if debug
        println("float 1) $(Char(b + UInt8('0')))")
    end
    if b > 0x09
        # character isn't a digit, check for special values, otherwise INVALID
        b += UInt8('0')
        if b == UInt8('n') || b == UInt8('N')
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto invalid
            end
            b = peekbyte(source, pos)
            if b == UInt8('a') || b == UInt8('A')
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto invalid
                end
                b = peekbyte(source, pos)
                if b == UInt8('n') || b == UInt8('N')
                    x = T(NaN)
                    pos += 1
                    incr!(source)
                    if eof(source, pos, len)
                        code |= EOF
                    end
                    code |= OK
                    @goto done
                end
            end
        elseif b == UInt8('i') || b == UInt8('I')
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                code |= EOF
                @goto invalid
            end
            b = peekbyte(source, pos)
            if b == UInt8('n') || b == UInt8('N')
                pos += 1
                incr!(source)
                if eof(source, pos, len)
                    code |= EOF
                    @goto invalid
                end
                b = peekbyte(source, pos)
                if b == UInt8('f') || b == UInt8('F')
                    x = ifelse(neg, T(-Inf), T(Inf))
                    code |= OK
                    pos += 1
                    incr!(source)
                    if eof(source, pos, len)
                        code |= EOF
                        @goto done
                    end
                    b = peekbyte(source, pos)
                    if b == UInt8('i') || b == UInt8('I')
                        pos += 1
                        incr!(source)
                        if eof(source, pos, len)
                            code |= EOF
                            @goto done
                        end
                        b = peekbyte(source, pos)
                        if b == UInt8('n') || b == UInt8('N')
                            pos += 1
                            incr!(source)
                            if eof(source, pos, len)
                                code |= EOF
                                @goto done
                            end
                            b = peekbyte(source, pos)
                            if b == UInt8('i') || b == UInt8('I')
                                pos += 1
                                incr!(source)
                                if eof(source, pos, len)
                                    code |= EOF
                                    @goto done
                                end
                                b = peekbyte(source, pos)
                                if b == UInt8('t') || b == UInt8('T')
                                    pos += 1
                                    incr!(source)
                                    if eof(source, pos, len)
                                        code |= EOF
                                        @goto done
                                    end
                                    b = peekbyte(source, pos)
                                    if b == UInt8('y') || b == UInt8('Y')
                                        pos += 1
                                        incr!(source)
                                        if eof(source, pos, len)
                                            code |= EOF
                                        end
                                        @goto done
                                    end
                                end
                            end
                        end
                    end
                    @goto done
                end
            end
        else
        end
@label invalid
        fastseek!(source, startpos)
        code |= INVALID
        @goto done
    end
    while true
        digits = IntType(10) * digits + b
        pos += 1
        incr!(source)
        ndigits += 1
        if eof(source, pos, len)
            digits = inttype(IntType)(digits)
            x = ifelse(neg, -T(digits), T(digits))
            code |= OK | EOF
            @goto done
        end
        b = peekbyte(source, pos) - UInt8('0')
        if debug
            println("float 2) $(Char(b + UInt8('0')))")
        end
        b > 0x09 && break
        if overflows(IntType) && digits > overflowval(IntType)
            fastseek!(source, startpos)
            return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
        elseif ndigits > maxdigits(T)
            fastseek!(source, startpos)
            code |= INVALID
            @goto done
        end
    end
    b += UInt8('0')
    if debug
        println("float 3) $(Char(b))")
    end
@label decimalcheck
    if b == options.decimal
        pos += 1
        incr!(source)
        if eof(source, pos, len)
            x = ifelse(neg, -T(digits), T(digits))
            code |= ((startpos + 1) == pos ? INVALID : OK) | EOF
            @goto done
        end
        b = peekbyte(source, pos)
        if b - UInt8('0') > 0x09 && !(b == UInt8('e') || b == UInt8('E') || b == UInt8('f') || b == UInt8('F'))
            if ndigits == 0
                code |= INVALID
                @goto done
            else
                x = ifelse(neg, -T(digits), T(digits))
                code |= OK
                @goto done
            end
        end
    end
    frac = inttype(IntType)(0)
    if debug
        println("float 4) $(Char(b))")
    end
    b -= UInt8('0')
    if b < 0x0a
        while true
            digits = IntType(10) * digits + b
            pos += 1
            incr!(source)
            frac += inttype(IntType)(1)
            if eof(source, pos, len)
                x = scale(T, digits, -frac, neg)
                code |= OK | EOF
                @goto done
            end
            b = peekbyte(source, pos) - UInt8('0')
            if debug
                println("float 5) $b")
            end
            b > 0x09 && break
            if overflows(IntType) && digits > overflowval(IntType)
                fastseek!(source, startpos)
                return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
            end
        end
    end
    b += UInt8('0')
    if debug
        println("float 6) $(Char(b))")
    end
    # check for exponent notation
    if b == UInt8('e') || b == UInt8('E') || b == UInt8('f') || b == UInt8('F')
        pos += 1
        incr!(source)
        # error to have a "dangling" 'e'
        if eof(source, pos, len)
            code |= INVALID | EOF
            @goto done
        end
        b = peekbyte(source, pos)
        if debug
            println("float 7) $(Char(b))")
        end
        exp = inttype(IntType)(0)
        negexp = b == UInt8('-')
        if negexp || b == UInt8('+')
            pos += 1
            incr!(source)
        end
        b = peekbyte(source, pos) - UInt8('0')
        if debug
            println("float 8) $b")
        end
        if b > 0x09
            # invalid to have a "dangling" 'e'
            code |= INVALID
            @goto done
        end
        while true
            exp = inttype(IntType)(10) * exp + b
            pos += 1
            incr!(source)
            if eof(source, pos, len)
                x = scale(T, digits, ifelse(negexp, -exp, exp) - frac, neg)
                code |= OK | EOF
                @goto done
            end
            b = peekbyte(source, pos) - UInt8('0')
            if debug
                println("float 9) $b")
            end
            if b > 0x09
                x = scale(T, digits, ifelse(negexp, -exp, exp) - frac, neg)
                code |= OK
                @goto done
            end
            if overflows(IntType) && exp > overflowval(IntType)
                fastseek!(source, startpos)
                return _typeparser(T, source, startpos, len, origb, code, options, wider(IntType))
            end
        end
    else
        x = scale(T, digits, -frac, neg)
        code |= OK
    end

@label done
    return x, code, pos
end

using Base.MPFR, Base.GMP, Base.GMP.MPZ

maxsig(::Type{Float16}) = 2048
maxsig(::Type{Float32}) = 16777216
maxsig(::Type{Float64}) = 9007199254740992
maxsig(::Type{BigFloat}) = 9007199254740992

ceillog5(::Type{Float16}) = 5
ceillog5(::Type{Float32}) = 11
ceillog5(::Type{Float64}) = 23
ceillog5(::Type{BigFloat}) = 23

const F16_SHORT_POWERS = [exp10(Float16(x)) for x = 0:2ceillog5(Float16)-1]
const F32_SHORT_POWERS = [exp10(Float32(x)) for x = 0:2ceillog5(Float32)-1]
const F64_SHORT_POWERS = [exp10(Float64(x)) for x = 0:2ceillog5(Float64)-1]

pow10(::Type{Float16}, e) = (@inbounds v = F16_SHORT_POWERS[e+1]; return v)
pow10(::Type{Float32}, e) = (@inbounds v = F32_SHORT_POWERS[e+1]; return v)
pow10(::Type{Float64}, e) = (@inbounds v = F64_SHORT_POWERS[e+1]; return v)
pow10(::Type{BigFloat}, e) = (@inbounds v = F64_SHORT_POWERS[e+1]; return v)

const BIGEXP10 = [1 / exp10(BigInt(e)) for e = 309:326]
const BIGFLOAT = [BigFloat()]

@inline function scale(::Type{T}, v::V, exp, neg) where {T, V <: Union{UInt128, BigInt}}
    if exp > 308
        return T(neg ? -Inf : Inf)
    elseif exp < -326
        return zero(T)
    elseif exp < -308
        y = BIGEXP10[-exp - 308]
        @inbounds x = BIGFLOAT[Threads.threadid()]
        if v > 9007199254740992
            ccall((:mpfr_set_z, :libmpfr), Int32,
                (Ref{BigFloat}, Ref{BigInt}, Int32),
                x, BigInt(v), MPFR.ROUNDING_MODE[])
        else
            ccall((:mpfr_set_d, :libmpfr), Int32,
                (Ref{BigFloat}, Float64, Int32),
                x, Float64(v), MPFR.ROUNDING_MODE[])
        end
        ccall((:mpfr_mul, :libmpfr), Int32,
            (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, Int32),
            x, x, y, MPFR.ROUNDING_MODE[])
    elseif exp < 0
        # this and the next branch are pretty slow for BigInt
        # but also extremely rare (i.e. floats w/ > 39 digits)
        x = v / exp10(V(-exp))
    else
        x = v * exp10(V(exp))
    end
    return T(neg ? -x : x)
end

@inline function scale(::Type{T}, v::UInt64, exp, neg) where {T}
    ms = maxsig(T)
    cl = ceillog5(T)
    if v < ms
        # fastest path
        if 0 <= exp < cl
            x = T(v) * pow10(T, exp)
            return ifelse(neg, -x, x)
        elseif -cl < exp < 0
            x = T(v) / pow10(T, -exp)
            return ifelse(neg, -x, x)
        end
    end
    v == 0 && return zero(T)
    if exp > 308
        return T(neg ? -Inf : Inf)
    elseif exp < -326
        return zero(T)
    end
    mant, pow = pow10spl(exp + 326)
    lz = leading_zeros(v)
    newv = v << lz
    upper, lower = two_prod(newv, mant)
    if upper & 0x1FF == 0x1FF && (lower + newv < lower)
        mant_low = mant128(exp + 325)
        product_middle2, product_low = two_prod(newv, mant_low)
        product_middle1 = lower
        product_high = upper
        product_middle = product_middle1 + product_middle2
        if product_middle < product_middle1
            product_high += 1 # overflow carry
        end
        if (product_middle + 1 == 0) && (product_high & 0x1FF == 0x1FF) &&
            (product_low + v < product_low)
            return scale(T, UInt128(v), exp, neg)
        end
        upper = product_high
        lower = product_middle
    end
    upperbit = upper >> 63
    mantissa = upper >> (upperbit + 9)
    lz += xor(1, upperbit)
    if (lower == 0) && upper & 0x1FF == 0 && mantissa & 3 == 1
        return scale(T, UInt128(v), exp, neg)
    end
    mantissa += mantissa & 1
    mantissa >>= 1
    if mantissa >= UInt64(1) << 53
        mantissa = UInt64(1) << 52
        lz -= 1
    end
    mantissa &= ~(UInt64(1) << 52)
    real_exponent = pow - lz
    if real_exponent < 1 || real_exponent > 2046
        return scale(T, UInt128(v), exp, neg)
    end
    mantissa |= real_exponent << 52
    mantissa |= (UInt64(neg) << 63)
    return T(Core.bitcast(Float64, mantissa))
end

@inline function two_prod(a, b)
    x = UInt128(a) * b
    return UInt64(x >> 64), unsafe_trunc(UInt64, x)
end

pow10spl(i) = (@inbounds x = POW10SPL[i + 1]; return x)
mant128(i) = (@inbounds x = MANTISSA128[i + 1]; return x)

const POW10SPL = [
    (0xa5ced43b7e3e9188, 7),
    (0xa5ced43b7e3e9188, 7),    (0xcf42894a5dce35ea, 10),
    (0x818995ce7aa0e1b2, 14),   (0xa1ebfb4219491a1f, 17),
    (0xca66fa129f9b60a6, 20),   (0xfd00b897478238d0, 23),
    (0x9e20735e8cb16382, 27),   (0xc5a890362fddbc62, 30),
    (0xf712b443bbd52b7b, 33),   (0x9a6bb0aa55653b2d, 37),
    (0xc1069cd4eabe89f8, 40),   (0xf148440a256e2c76, 43),
    (0x96cd2a865764dbca, 47),   (0xbc807527ed3e12bc, 50),
    (0xeba09271e88d976b, 53),   (0x93445b8731587ea3, 57),
    (0xb8157268fdae9e4c, 60),   (0xe61acf033d1a45df, 63),
    (0x8fd0c16206306bab, 67),   (0xb3c4f1ba87bc8696, 70),
    (0xe0b62e2929aba83c, 73),   (0x8c71dcd9ba0b4925, 77),
    (0xaf8e5410288e1b6f, 80),   (0xdb71e91432b1a24a, 83),
    (0x892731ac9faf056e, 87),   (0xab70fe17c79ac6ca, 90),
    (0xd64d3d9db981787d, 93),   (0x85f0468293f0eb4e, 97),
    (0xa76c582338ed2621, 100),  (0xd1476e2c07286faa, 103),
    (0x82cca4db847945ca, 107),  (0xa37fce126597973c, 110),
    (0xcc5fc196fefd7d0c, 113),  (0xff77b1fcbebcdc4f, 116),
    (0x9faacf3df73609b1, 120),  (0xc795830d75038c1d, 123),
    (0xf97ae3d0d2446f25, 126),  (0x9becce62836ac577, 130),
    (0xc2e801fb244576d5, 133),  (0xf3a20279ed56d48a, 136),
    (0x9845418c345644d6, 140),  (0xbe5691ef416bd60c, 143),
    (0xedec366b11c6cb8f, 146),  (0x94b3a202eb1c3f39, 150),
    (0xb9e08a83a5e34f07, 153),  (0xe858ad248f5c22c9, 156),
    (0x91376c36d99995be, 160),  (0xb58547448ffffb2d, 163),
    (0xe2e69915b3fff9f9, 166),  (0x8dd01fad907ffc3b, 170),
    (0xb1442798f49ffb4a, 173),  (0xdd95317f31c7fa1d, 176),
    (0x8a7d3eef7f1cfc52, 180),  (0xad1c8eab5ee43b66, 183),
    (0xd863b256369d4a40, 186),  (0x873e4f75e2224e68, 190),
    (0xa90de3535aaae202, 193),  (0xd3515c2831559a83, 196),
    (0x8412d9991ed58091, 200),  (0xa5178fff668ae0b6, 203),
    (0xce5d73ff402d98e3, 206),  (0x80fa687f881c7f8e, 210),
    (0xa139029f6a239f72, 213),  (0xc987434744ac874e, 216),
    (0xfbe9141915d7a922, 219),  (0x9d71ac8fada6c9b5, 223),
    (0xc4ce17b399107c22, 226),  (0xf6019da07f549b2b, 229),
    (0x99c102844f94e0fb, 233),  (0xc0314325637a1939, 236),
    (0xf03d93eebc589f88, 239),  (0x96267c7535b763b5, 243),
    (0xbbb01b9283253ca2, 246),  (0xea9c227723ee8bcb, 249),
    (0x92a1958a7675175f, 253),  (0xb749faed14125d36, 256),
    (0xe51c79a85916f484, 259),  (0x8f31cc0937ae58d2, 263),
    (0xb2fe3f0b8599ef07, 266),  (0xdfbdcece67006ac9, 269),
    (0x8bd6a141006042bd, 273),  (0xaecc49914078536d, 276),
    (0xda7f5bf590966848, 279),  (0x888f99797a5e012d, 283),
    (0xaab37fd7d8f58178, 286),  (0xd5605fcdcf32e1d6, 289),
    (0x855c3be0a17fcd26, 293),  (0xa6b34ad8c9dfc06f, 296),
    (0xd0601d8efc57b08b, 299),  (0x823c12795db6ce57, 303),
    (0xa2cb1717b52481ed, 306),  (0xcb7ddcdda26da268, 309),
    (0xfe5d54150b090b02, 312),  (0x9efa548d26e5a6e1, 316),
    (0xc6b8e9b0709f109a, 319),  (0xf867241c8cc6d4c0, 322),
    (0x9b407691d7fc44f8, 326),  (0xc21094364dfb5636, 329),
    (0xf294b943e17a2bc4, 332),  (0x979cf3ca6cec5b5a, 336),
    (0xbd8430bd08277231, 339),  (0xece53cec4a314ebd, 342),
    (0x940f4613ae5ed136, 346),  (0xb913179899f68584, 349),
    (0xe757dd7ec07426e5, 352),  (0x9096ea6f3848984f, 356),
    (0xb4bca50b065abe63, 359),  (0xe1ebce4dc7f16dfb, 362),
    (0x8d3360f09cf6e4bd, 366),  (0xb080392cc4349dec, 369),
    (0xdca04777f541c567, 372),  (0x89e42caaf9491b60, 376),
    (0xac5d37d5b79b6239, 379),  (0xd77485cb25823ac7, 382),
    (0x86a8d39ef77164bc, 386),  (0xa8530886b54dbdeb, 389),
    (0xd267caa862a12d66, 392),  (0x8380dea93da4bc60, 396),
    (0xa46116538d0deb78, 399),  (0xcd795be870516656, 402),
    (0x806bd9714632dff6, 406),  (0xa086cfcd97bf97f3, 409),
    (0xc8a883c0fdaf7df0, 412),  (0xfad2a4b13d1b5d6c, 415),
    (0x9cc3a6eec6311a63, 419),  (0xc3f490aa77bd60fc, 422),
    (0xf4f1b4d515acb93b, 425),  (0x991711052d8bf3c5, 429),
    (0xbf5cd54678eef0b6, 432),  (0xef340a98172aace4, 435),
    (0x9580869f0e7aac0e, 439),  (0xbae0a846d2195712, 442),
    (0xe998d258869facd7, 445),  (0x91ff83775423cc06, 449),
    (0xb67f6455292cbf08, 452),  (0xe41f3d6a7377eeca, 455),
    (0x8e938662882af53e, 459),  (0xb23867fb2a35b28d, 462),
    (0xdec681f9f4c31f31, 465),  (0x8b3c113c38f9f37e, 469),
    (0xae0b158b4738705e, 472),  (0xd98ddaee19068c76, 475),
    (0x87f8a8d4cfa417c9, 479),  (0xa9f6d30a038d1dbc, 482),
    (0xd47487cc8470652b, 485),  (0x84c8d4dfd2c63f3b, 489),
    (0xa5fb0a17c777cf09, 492),  (0xcf79cc9db955c2cc, 495),
    (0x81ac1fe293d599bf, 499),  (0xa21727db38cb002f, 502),
    (0xca9cf1d206fdc03b, 505),  (0xfd442e4688bd304a, 508),
    (0x9e4a9cec15763e2e, 512),  (0xc5dd44271ad3cdba, 515),
    (0xf7549530e188c128, 518),  (0x9a94dd3e8cf578b9, 522),
    (0xc13a148e3032d6e7, 525),  (0xf18899b1bc3f8ca1, 528),
    (0x96f5600f15a7b7e5, 532),  (0xbcb2b812db11a5de, 535),
    (0xebdf661791d60f56, 538),  (0x936b9fcebb25c995, 542),
    (0xb84687c269ef3bfb, 545),  (0xe65829b3046b0afa, 548),
    (0x8ff71a0fe2c2e6dc, 552),  (0xb3f4e093db73a093, 555),
    (0xe0f218b8d25088b8, 558),  (0x8c974f7383725573, 562),
    (0xafbd2350644eeacf, 565),  (0xdbac6c247d62a583, 568),
    (0x894bc396ce5da772, 572),  (0xab9eb47c81f5114f, 575),
    (0xd686619ba27255a2, 578),  (0x8613fd0145877585, 582),
    (0xa798fc4196e952e7, 585),  (0xd17f3b51fca3a7a0, 588),
    (0x82ef85133de648c4, 592),  (0xa3ab66580d5fdaf5, 595),
    (0xcc963fee10b7d1b3, 598),  (0xffbbcfe994e5c61f, 601),
    (0x9fd561f1fd0f9bd3, 605),  (0xc7caba6e7c5382c8, 608),
    (0xf9bd690a1b68637b, 611),  (0x9c1661a651213e2d, 615),
    (0xc31bfa0fe5698db8, 618),  (0xf3e2f893dec3f126, 621),
    (0x986ddb5c6b3a76b7, 625),  (0xbe89523386091465, 628),
    (0xee2ba6c0678b597f, 631),  (0x94db483840b717ef, 635),
    (0xba121a4650e4ddeb, 638),  (0xe896a0d7e51e1566, 641),
    (0x915e2486ef32cd60, 645),  (0xb5b5ada8aaff80b8, 648),
    (0xe3231912d5bf60e6, 651),  (0x8df5efabc5979c8f, 655),
    (0xb1736b96b6fd83b3, 658),  (0xddd0467c64bce4a0, 661),
    (0x8aa22c0dbef60ee4, 665),  (0xad4ab7112eb3929d, 668),
    (0xd89d64d57a607744, 671),  (0x87625f056c7c4a8b, 675),
    (0xa93af6c6c79b5d2d, 678),  (0xd389b47879823479, 681),
    (0x843610cb4bf160cb, 685),  (0xa54394fe1eedb8fe, 688),
    (0xce947a3da6a9273e, 691),  (0x811ccc668829b887, 695),
    (0xa163ff802a3426a8, 698),  (0xc9bcff6034c13052, 701),
    (0xfc2c3f3841f17c67, 704),  (0x9d9ba7832936edc0, 708),
    (0xc5029163f384a931, 711),  (0xf64335bcf065d37d, 714),
    (0x99ea0196163fa42e, 718),  (0xc06481fb9bcf8d39, 721),
    (0xf07da27a82c37088, 724),  (0x964e858c91ba2655, 728),
    (0xbbe226efb628afea, 731),  (0xeadab0aba3b2dbe5, 734),
    (0x92c8ae6b464fc96f, 738),  (0xb77ada0617e3bbcb, 741),
    (0xe55990879ddcaabd, 744),  (0x8f57fa54c2a9eab6, 748),
    (0xb32df8e9f3546564, 751),  (0xdff9772470297ebd, 754),
    (0x8bfbea76c619ef36, 758),  (0xaefae51477a06b03, 761),
    (0xdab99e59958885c4, 764),  (0x88b402f7fd75539b, 768),
    (0xaae103b5fcd2a881, 771),  (0xd59944a37c0752a2, 774),
    (0x857fcae62d8493a5, 778),  (0xa6dfbd9fb8e5b88e, 781),
    (0xd097ad07a71f26b2, 784),  (0x825ecc24c873782f, 788),
    (0xa2f67f2dfa90563b, 791),  (0xcbb41ef979346bca, 794),
    (0xfea126b7d78186bc, 797),  (0x9f24b832e6b0f436, 801),
    (0xc6ede63fa05d3143, 804),  (0xf8a95fcf88747d94, 807),
    (0x9b69dbe1b548ce7c, 811),  (0xc24452da229b021b, 814),
    (0xf2d56790ab41c2a2, 817),  (0x97c560ba6b0919a5, 821),
    (0xbdb6b8e905cb600f, 824),  (0xed246723473e3813, 827),
    (0x9436c0760c86e30b, 831),  (0xb94470938fa89bce, 834),
    (0xe7958cb87392c2c2, 837),  (0x90bd77f3483bb9b9, 841),
    (0xb4ecd5f01a4aa828, 844),  (0xe2280b6c20dd5232, 847),
    (0x8d590723948a535f, 851),  (0xb0af48ec79ace837, 854),
    (0xdcdb1b2798182244, 857),  (0x8a08f0f8bf0f156b, 861),
    (0xac8b2d36eed2dac5, 864),  (0xd7adf884aa879177, 867),
    (0x86ccbb52ea94baea, 871),  (0xa87fea27a539e9a5, 874),
    (0xd29fe4b18e88640e, 877),  (0x83a3eeeef9153e89, 881),
    (0xa48ceaaab75a8e2b, 884),  (0xcdb02555653131b6, 887),
    (0x808e17555f3ebf11, 891),  (0xa0b19d2ab70e6ed6, 894),
    (0xc8de047564d20a8b, 897),  (0xfb158592be068d2e, 900),
    (0x9ced737bb6c4183d, 904),  (0xc428d05aa4751e4c, 907),
    (0xf53304714d9265df, 910),  (0x993fe2c6d07b7fab, 914),
    (0xbf8fdb78849a5f96, 917),  (0xef73d256a5c0f77c, 920),
    (0x95a8637627989aad, 924),  (0xbb127c53b17ec159, 927),
    (0xe9d71b689dde71af, 930),  (0x9226712162ab070d, 934),
    (0xb6b00d69bb55c8d1, 937),  (0xe45c10c42a2b3b05, 940),
    (0x8eb98a7a9a5b04e3, 944),  (0xb267ed1940f1c61c, 947),
    (0xdf01e85f912e37a3, 950),  (0x8b61313bbabce2c6, 954),
    (0xae397d8aa96c1b77, 957),  (0xd9c7dced53c72255, 960),
    (0x881cea14545c7575, 964),  (0xaa242499697392d2, 967),
    (0xd4ad2dbfc3d07787, 970),  (0x84ec3c97da624ab4, 974),
    (0xa6274bbdd0fadd61, 977),  (0xcfb11ead453994ba, 980),
    (0x81ceb32c4b43fcf4, 984),  (0xa2425ff75e14fc31, 987),
    (0xcad2f7f5359a3b3e, 990),  (0xfd87b5f28300ca0d, 993),
    (0x9e74d1b791e07e48, 997),  (0xc612062576589dda, 1000),
    (0xf79687aed3eec551, 1003), (0x9abe14cd44753b52, 1007),
    (0xc16d9a0095928a27, 1010), (0xf1c90080baf72cb1, 1013),
    (0x971da05074da7bee, 1017), (0xbce5086492111aea, 1020),
    (0xec1e4a7db69561a5, 1023), (0x9392ee8e921d5d07, 1027),
    (0xb877aa3236a4b449, 1030), (0xe69594bec44de15b, 1033),
    (0x901d7cf73ab0acd9, 1037), (0xb424dc35095cd80f, 1040),
    (0xe12e13424bb40e13, 1043), (0x8cbccc096f5088cb, 1047),
    (0xafebff0bcb24aafe, 1050), (0xdbe6fecebdedd5be, 1053),
    (0x89705f4136b4a597, 1057), (0xabcc77118461cefc, 1060),
    (0xd6bf94d5e57a42bc, 1063), (0x8637bd05af6c69b5, 1067),
    (0xa7c5ac471b478423, 1070), (0xd1b71758e219652b, 1073),
    (0x83126e978d4fdf3b, 1077), (0xa3d70a3d70a3d70a, 1080),
    (0xcccccccccccccccc, 1083), (0x8000000000000000, 1087),
    (0xa000000000000000, 1090), (0xc800000000000000, 1093),
    (0xfa00000000000000, 1096), (0x9c40000000000000, 1100),
    (0xc350000000000000, 1103), (0xf424000000000000, 1106),
    (0x9896800000000000, 1110), (0xbebc200000000000, 1113),
    (0xee6b280000000000, 1116), (0x9502f90000000000, 1120),
    (0xba43b74000000000, 1123), (0xe8d4a51000000000, 1126),
    (0x9184e72a00000000, 1130), (0xb5e620f480000000, 1133),
    (0xe35fa931a0000000, 1136), (0x8e1bc9bf04000000, 1140),
    (0xb1a2bc2ec5000000, 1143), (0xde0b6b3a76400000, 1146),
    (0x8ac7230489e80000, 1150), (0xad78ebc5ac620000, 1153),
    (0xd8d726b7177a8000, 1156), (0x878678326eac9000, 1160),
    (0xa968163f0a57b400, 1163), (0xd3c21bcecceda100, 1166),
    (0x84595161401484a0, 1170), (0xa56fa5b99019a5c8, 1173),
    (0xcecb8f27f4200f3a, 1176), (0x813f3978f8940984, 1180),
    (0xa18f07d736b90be5, 1183), (0xc9f2c9cd04674ede, 1186),
    (0xfc6f7c4045812296, 1189), (0x9dc5ada82b70b59d, 1193),
    (0xc5371912364ce305, 1196), (0xf684df56c3e01bc6, 1199),
    (0x9a130b963a6c115c, 1203), (0xc097ce7bc90715b3, 1206),
    (0xf0bdc21abb48db20, 1209), (0x96769950b50d88f4, 1213),
    (0xbc143fa4e250eb31, 1216), (0xeb194f8e1ae525fd, 1219),
    (0x92efd1b8d0cf37be, 1223), (0xb7abc627050305ad, 1226),
    (0xe596b7b0c643c719, 1229), (0x8f7e32ce7bea5c6f, 1233),
    (0xb35dbf821ae4f38b, 1236), (0xe0352f62a19e306e, 1239),
    (0x8c213d9da502de45, 1243), (0xaf298d050e4395d6, 1246),
    (0xdaf3f04651d47b4c, 1249), (0x88d8762bf324cd0f, 1253),
    (0xab0e93b6efee0053, 1256), (0xd5d238a4abe98068, 1259),
    (0x85a36366eb71f041, 1263), (0xa70c3c40a64e6c51, 1266),
    (0xd0cf4b50cfe20765, 1269), (0x82818f1281ed449f, 1273),
    (0xa321f2d7226895c7, 1276), (0xcbea6f8ceb02bb39, 1279),
    (0xfee50b7025c36a08, 1282), (0x9f4f2726179a2245, 1286),
    (0xc722f0ef9d80aad6, 1289), (0xf8ebad2b84e0d58b, 1292),
    (0x9b934c3b330c8577, 1296), (0xc2781f49ffcfa6d5, 1299),
    (0xf316271c7fc3908a, 1302), (0x97edd871cfda3a56, 1306),
    (0xbde94e8e43d0c8ec, 1309), (0xed63a231d4c4fb27, 1312),
    (0x945e455f24fb1cf8, 1316), (0xb975d6b6ee39e436, 1319),
    (0xe7d34c64a9c85d44, 1322), (0x90e40fbeea1d3a4a, 1326),
    (0xb51d13aea4a488dd, 1329), (0xe264589a4dcdab14, 1332),
    (0x8d7eb76070a08aec, 1336), (0xb0de65388cc8ada8, 1339),
    (0xdd15fe86affad912, 1342), (0x8a2dbf142dfcc7ab, 1346),
    (0xacb92ed9397bf996, 1349), (0xd7e77a8f87daf7fb, 1352),
    (0x86f0ac99b4e8dafd, 1356), (0xa8acd7c0222311bc, 1359),
    (0xd2d80db02aabd62b, 1362), (0x83c7088e1aab65db, 1366),
    (0xa4b8cab1a1563f52, 1369), (0xcde6fd5e09abcf26, 1372),
    (0x80b05e5ac60b6178, 1376), (0xa0dc75f1778e39d6, 1379),
    (0xc913936dd571c84c, 1382), (0xfb5878494ace3a5f, 1385),
    (0x9d174b2dcec0e47b, 1389), (0xc45d1df942711d9a, 1392),
    (0xf5746577930d6500, 1395), (0x9968bf6abbe85f20, 1399),
    (0xbfc2ef456ae276e8, 1402), (0xefb3ab16c59b14a2, 1405),
    (0x95d04aee3b80ece5, 1409), (0xbb445da9ca61281f, 1412),
    (0xea1575143cf97226, 1415), (0x924d692ca61be758, 1419),
    (0xb6e0c377cfa2e12e, 1422), (0xe498f455c38b997a, 1425),
    (0x8edf98b59a373fec, 1429), (0xb2977ee300c50fe7, 1432),
    (0xdf3d5e9bc0f653e1, 1435), (0x8b865b215899f46c, 1439),
    (0xae67f1e9aec07187, 1442), (0xda01ee641a708de9, 1445),
    (0x884134fe908658b2, 1449), (0xaa51823e34a7eede, 1452),
    (0xd4e5e2cdc1d1ea96, 1455), (0x850fadc09923329e, 1459),
    (0xa6539930bf6bff45, 1462), (0xcfe87f7cef46ff16, 1465),
    (0x81f14fae158c5f6e, 1469), (0xa26da3999aef7749, 1472),
    (0xcb090c8001ab551c, 1475), (0xfdcb4fa002162a63, 1478),
    (0x9e9f11c4014dda7e, 1482), (0xc646d63501a1511d, 1485),
    (0xf7d88bc24209a565, 1488), (0x9ae757596946075f, 1492),
    (0xc1a12d2fc3978937, 1495), (0xf209787bb47d6b84, 1498),
    (0x9745eb4d50ce6332, 1502), (0xbd176620a501fbff, 1505),
    (0xec5d3fa8ce427aff, 1508), (0x93ba47c980e98cdf, 1512),
    (0xb8a8d9bbe123f017, 1515), (0xe6d3102ad96cec1d, 1518),
    (0x9043ea1ac7e41392, 1522), (0xb454e4a179dd1877, 1525),
    (0xe16a1dc9d8545e94, 1528), (0x8ce2529e2734bb1d, 1532),
    (0xb01ae745b101e9e4, 1535), (0xdc21a1171d42645d, 1538),
    (0x899504ae72497eba, 1542), (0xabfa45da0edbde69, 1545),
    (0xd6f8d7509292d603, 1548), (0x865b86925b9bc5c2, 1552),
    (0xa7f26836f282b732, 1555), (0xd1ef0244af2364ff, 1558),
    (0x8335616aed761f1f, 1562), (0xa402b9c5a8d3a6e7, 1565),
    (0xcd036837130890a1, 1568), (0x802221226be55a64, 1572),
    (0xa02aa96b06deb0fd, 1575), (0xc83553c5c8965d3d, 1578),
    (0xfa42a8b73abbf48c, 1581), (0x9c69a97284b578d7, 1585),
    (0xc38413cf25e2d70d, 1588), (0xf46518c2ef5b8cd1, 1591),
    (0x98bf2f79d5993802, 1595), (0xbeeefb584aff8603, 1598),
    (0xeeaaba2e5dbf6784, 1601), (0x952ab45cfa97a0b2, 1605),
    (0xba756174393d88df, 1608), (0xe912b9d1478ceb17, 1611),
    (0x91abb422ccb812ee, 1615), (0xb616a12b7fe617aa, 1618),
    (0xe39c49765fdf9d94, 1621), (0x8e41ade9fbebc27d, 1625),
    (0xb1d219647ae6b31c, 1628), (0xde469fbd99a05fe3, 1631),
    (0x8aec23d680043bee, 1635), (0xada72ccc20054ae9, 1638),
    (0xd910f7ff28069da4, 1641), (0x87aa9aff79042286, 1645),
    (0xa99541bf57452b28, 1648), (0xd3fa922f2d1675f2, 1651),
    (0x847c9b5d7c2e09b7, 1655), (0xa59bc234db398c25, 1658),
    (0xcf02b2c21207ef2e, 1661), (0x8161afb94b44f57d, 1665),
    (0xa1ba1ba79e1632dc, 1668), (0xca28a291859bbf93, 1671),
    (0xfcb2cb35e702af78, 1674), (0x9defbf01b061adab, 1678),
    (0xc56baec21c7a1916, 1681), (0xf6c69a72a3989f5b, 1684),
    (0x9a3c2087a63f6399, 1688), (0xc0cb28a98fcf3c7f, 1691),
    (0xf0fdf2d3f3c30b9f, 1694), (0x969eb7c47859e743, 1698),
    (0xbc4665b596706114, 1701), (0xeb57ff22fc0c7959, 1704),
    (0x9316ff75dd87cbd8, 1708), (0xb7dcbf5354e9bece, 1711),
    (0xe5d3ef282a242e81, 1714), (0x8fa475791a569d10, 1718),
    (0xb38d92d760ec4455, 1721), (0xe070f78d3927556a, 1724),
    (0x8c469ab843b89562, 1728), (0xaf58416654a6babb, 1731),
    (0xdb2e51bfe9d0696a, 1734), (0x88fcf317f22241e2, 1738),
    (0xab3c2fddeeaad25a, 1741), (0xd60b3bd56a5586f1, 1744),
    (0x85c7056562757456, 1748), (0xa738c6bebb12d16c, 1751),
    (0xd106f86e69d785c7, 1754), (0x82a45b450226b39c, 1758),
    (0xa34d721642b06084, 1761), (0xcc20ce9bd35c78a5, 1764),
    (0xff290242c83396ce, 1767), (0x9f79a169bd203e41, 1771),
    (0xc75809c42c684dd1, 1774), (0xf92e0c3537826145, 1777),
    (0x9bbcc7a142b17ccb, 1781), (0xc2abf989935ddbfe, 1784),
    (0xf356f7ebf83552fe, 1787), (0x98165af37b2153de, 1791),
    (0xbe1bf1b059e9a8d6, 1794), (0xeda2ee1c7064130c, 1797),
    (0x9485d4d1c63e8be7, 1801), (0xb9a74a0637ce2ee1, 1804),
    (0xe8111c87c5c1ba99, 1807), (0x910ab1d4db9914a0, 1811),
    (0xb54d5e4a127f59c8, 1814), (0xe2a0b5dc971f303a, 1817),
    (0x8da471a9de737e24, 1821), (0xb10d8e1456105dad, 1824),
    (0xdd50f1996b947518, 1827), (0x8a5296ffe33cc92f, 1831),
    (0xace73cbfdc0bfb7b, 1834), (0xd8210befd30efa5a, 1837),
    (0x8714a775e3e95c78, 1841), (0xa8d9d1535ce3b396, 1844),
    (0xd31045a8341ca07c, 1847), (0x83ea2b892091e44d, 1851),
    (0xa4e4b66b68b65d60, 1854), (0xce1de40642e3f4b9, 1857),
    (0x80d2ae83e9ce78f3, 1861), (0xa1075a24e4421730, 1864),
    (0xc94930ae1d529cfc, 1867), (0xfb9b7cd9a4a7443c, 1870),
    (0x9d412e0806e88aa5, 1874), (0xc491798a08a2ad4e, 1877),
    (0xf5b5d7ec8acb58a2, 1880), (0x9991a6f3d6bf1765, 1884),
    (0xbff610b0cc6edd3f, 1887), (0xeff394dcff8a948e, 1890),
    (0x95f83d0a1fb69cd9, 1894), (0xbb764c4ca7a4440f, 1897),
    (0xea53df5fd18d5513, 1900), (0x92746b9be2f8552c, 1904),
    (0xb7118682dbb66a77, 1907), (0xe4d5e82392a40515, 1910),
    (0x8f05b1163ba6832d, 1914), (0xb2c71d5bca9023f8, 1917),
    (0xdf78e4b2bd342cf6, 1920), (0x8bab8eefb6409c1a, 1924),
    (0xae9672aba3d0c320, 1927), (0xda3c0f568cc4f3e8, 1930),
    (0x8865899617fb1871, 1934), (0xaa7eebfb9df9de8d, 1937),
    (0xd51ea6fa85785631, 1940), (0x8533285c936b35de, 1944),
    (0xa67ff273b8460356, 1947), (0xd01fef10a657842c, 1950),
    (0x8213f56a67f6b29b, 1954), (0xa298f2c501f45f42, 1957),
    (0xcb3f2f7642717713, 1960), (0xfe0efb53d30dd4d7, 1963),
    (0x9ec95d1463e8a506, 1967), (0xc67bb4597ce2ce48, 1970),
    (0xf81aa16fdc1b81da, 1973), (0x9b10a4e5e9913128, 1977),
    (0xc1d4ce1f63f57d72, 1980), (0xf24a01a73cf2dccf, 1983),
    (0x976e41088617ca01, 1987), (0xbd49d14aa79dbc82, 1990),
    (0xec9c459d51852ba2, 1993), (0x93e1ab8252f33b45, 1997),
    (0xb8da1662e7b00a17, 2000), (0xe7109bfba19c0c9d, 2003),
    (0x906a617d450187e2, 2007), (0xb484f9dc9641e9da, 2010),
    (0xe1a63853bbd26451, 2013), (0x8d07e33455637eb2, 2017),
    (0xb049dc016abc5e5f, 2020), (0xdc5c5301c56b75f7, 2023),
    (0x89b9b3e11b6329ba, 2027), (0xac2820d9623bf429, 2030),
    (0xd732290fbacaf133, 2033), (0x867f59a9d4bed6c0, 2037),
    (0xa81f301449ee8c70, 2040), (0xd226fc195c6a2f8c, 2043),
    (0x83585d8fd9c25db7, 2047), (0xa42e74f3d032f525, 2050),
    (0xcd3a1230c43fb26f, 2053), (0x80444b5e7aa7cf85, 2057),
    (0xa0555e361951c366, 2060), (0xc86ab5c39fa63440, 2063),
    (0xfa856334878fc150, 2066), (0x9c935e00d4b9d8d2, 2070),
    (0xc3b8358109e84f07, 2073), (0xf4a642e14c6262c8, 2076),
    (0x98e7e9cccfbd7dbd, 2080), (0xbf21e44003acdd2c, 2083),
    (0xeeea5d5004981478, 2086), (0x95527a5202df0ccb, 2090),
    (0xbaa718e68396cffd, 2093), (0xe950df20247c83fd, 2096),
    (0x91d28b7416cdd27e, 2100), (0xb6472e511c81471d, 2103),
    (0xe3d8f9e563a198e5, 2106), (0x8e679c2f5e44ff8f, 2110)]

const MANTISSA128 = [0x419ea3bd35385e2d,
    0x52064cac828675b9,
    0x7343efebd1940993,
    0x1014ebe6c5f90bf8,
    0xd41a26e077774ef6,
    0x8920b098955522b4,
    0x55b46e5f5d5535b0,
    0xeb2189f734aa831d,
    0xa5e9ec7501d523e4,
    0x47b233c92125366e,
    0x999ec0bb696e840a,
    0xc00670ea43ca250d,
    0x380406926a5e5728,
    0xc605083704f5ecf2,
    0xf7864a44c633682e,
    0x7ab3ee6afbe0211d,
    0x5960ea05bad82964,
    0x6fb92487298e33bd,
    0xa5d3b6d479f8e056,
    0x8f48a4899877186c,
    0x331acdabfe94de87,
    0x9ff0c08b7f1d0b14,
    0x7ecf0ae5ee44dd9,
    0xc9e82cd9f69d6150,
    0xbe311c083a225cd2,
    0x6dbd630a48aaf406,
    0x92cbbccdad5b108,
    0x25bbf56008c58ea5,
    0xaf2af2b80af6f24e,
    0x1af5af660db4aee1,
    0x50d98d9fc890ed4d,
    0xe50ff107bab528a0,
    0x1e53ed49a96272c8,
    0x25e8e89c13bb0f7a,
    0x77b191618c54e9ac,
    0xd59df5b9ef6a2417,
    0x4b0573286b44ad1d,
    0x4ee367f9430aec32,
    0x229c41f793cda73f,
    0x6b43527578c1110f,
    0x830a13896b78aaa9,
    0x23cc986bc656d553,
    0x2cbfbe86b7ec8aa8,
    0x7bf7d71432f3d6a9,
    0xdaf5ccd93fb0cc53,
    0xd1b3400f8f9cff68,
    0x23100809b9c21fa1,
    0xabd40a0c2832a78a,
    0x16c90c8f323f516c,
    0xae3da7d97f6792e3,
    0x99cd11cfdf41779c,
    0x40405643d711d583,
    0x482835ea666b2572,
    0xda3243650005eecf,
    0x90bed43e40076a82,
    0x5a7744a6e804a291,
    0x711515d0a205cb36,
    0xd5a5b44ca873e03,
    0xe858790afe9486c2,
    0x626e974dbe39a872,
    0xfb0a3d212dc8128f,
    0x7ce66634bc9d0b99,
    0x1c1fffc1ebc44e80,
    0xa327ffb266b56220,
    0x4bf1ff9f0062baa8,
    0x6f773fc3603db4a9,
    0xcb550fb4384d21d3,
    0x7e2a53a146606a48,
    0x2eda7444cbfc426d,
    0xfa911155fefb5308,
    0x793555ab7eba27ca,
    0x4bc1558b2f3458de,
    0x9eb1aaedfb016f16,
    0x465e15a979c1cadc,
    0xbfacd89ec191ec9,
    0xcef980ec671f667b,
    0x82b7e12780e7401a,
    0xd1b2ecb8b0908810,
    0x861fa7e6dcb4aa15,
    0x67a791e093e1d49a,
    0xe0c8bb2c5c6d24e0,
    0x58fae9f773886e18,
    0xaf39a475506a899e,
    0x6d8406c952429603,
    0xc8e5087ba6d33b83,
    0xfb1e4a9a90880a64,
    0x5cf2eea09a55067f,
    0xf42faa48c0ea481e,
    0xf13b94daf124da26,
    0x76c53d08d6b70858,
    0x54768c4b0c64ca6e,
    0xa9942f5dcf7dfd09,
    0xd3f93b35435d7c4c,
    0xc47bc5014a1a6daf,
    0x359ab6419ca1091b,
    0xc30163d203c94b62,
    0x79e0de63425dcf1d,
    0x985915fc12f542e4,
    0x3e6f5b7b17b2939d,
    0xa705992ceecf9c42,
    0x50c6ff782a838353,
    0xa4f8bf5635246428,
    0x871b7795e136be99,
    0x28e2557b59846e3f,
    0x331aeada2fe589cf,
    0x3ff0d2c85def7621,
    0xfed077a756b53a9,
    0xd3e8495912c62894,
    0x64712dd7abbbd95c,
    0xbd8d794d96aacfb3,
    0xecf0d7a0fc5583a0,
    0xf41686c49db57244,
    0x311c2875c522ced5,
    0x7d633293366b828b,
    0xae5dff9c02033197,
    0xd9f57f830283fdfc,
    0xd072df63c324fd7b,
    0x4247cb9e59f71e6d,
    0x52d9be85f074e608,
    0x67902e276c921f8b,
    0xba1cd8a3db53b6,
    0x80e8a40eccd228a4,
    0x6122cd128006b2cd,
    0x796b805720085f81,
    0xcbe3303674053bb0,
    0xbedbfc4411068a9c,
    0xee92fb5515482d44,
    0x751bdd152d4d1c4a,
    0xd262d45a78a0635d,
    0x86fb897116c87c34,
    0xd45d35e6ae3d4da0,
    0x8974836059cca109,
    0x2bd1a438703fc94b,
    0x7b6306a34627ddcf,
    0x1a3bc84c17b1d542,
    0x20caba5f1d9e4a93,
    0x547eb47b7282ee9c,
    0xe99e619a4f23aa43,
    0x6405fa00e2ec94d4,
    0xde83bc408dd3dd04,
    0x9624ab50b148d445,
    0x3badd624dd9b0957,
    0xe54ca5d70a80e5d6,
    0x5e9fcf4ccd211f4c,
    0x7647c3200069671f,
    0x29ecd9f40041e073,
    0xf468107100525890,
    0x7182148d4066eeb4,
    0xc6f14cd848405530,
    0xb8ada00e5a506a7c,
    0xa6d90811f0e4851c,
    0x908f4a166d1da663,
    0x9a598e4e043287fe,
    0x40eff1e1853f29fd,
    0xd12bee59e68ef47c,
    0x82bb74f8301958ce,
    0xe36a52363c1faf01,
    0xdc44e6c3cb279ac1,
    0x29ab103a5ef8c0b9,
    0x7415d448f6b6f0e7,
    0x111b495b3464ad21,
    0xcab10dd900beec34,
    0x3d5d514f40eea742,
    0xcb4a5a3112a5112,
    0x47f0e785eaba72ab,
    0x59ed216765690f56,
    0x306869c13ec3532c,
    0x1e414218c73a13fb,
    0xe5d1929ef90898fa,
    0xdf45f746b74abf39,
    0x6b8bba8c328eb783,
    0x66ea92f3f326564,
    0xc80a537b0efefebd,
    0xbd06742ce95f5f36,
    0x2c48113823b73704,
    0xf75a15862ca504c5,
    0x9a984d73dbe722fb,
    0xc13e60d0d2e0ebba,
    0x318df905079926a8,
    0xfdf17746497f7052,
    0xfeb6ea8bedefa633,
    0xfe64a52ee96b8fc0,
    0x3dfdce7aa3c673b0,
    0x6bea10ca65c084e,
    0x486e494fcff30a62,
    0x5a89dba3c3efccfa,
    0xf89629465a75e01c,
    0xf6bbb397f1135823,
    0x746aa07ded582e2c,
    0xa8c2a44eb4571cdc,
    0x92f34d62616ce413,
    0x77b020baf9c81d17,
    0xace1474dc1d122e,
    0xd819992132456ba,
    0x10e1fff697ed6c69,
    0xca8d3ffa1ef463c1,
    0xbd308ff8a6b17cb2,
    0xac7cb3f6d05ddbde,
    0x6bcdf07a423aa96b,
    0x86c16c98d2c953c6,
    0xe871c7bf077ba8b7,
    0x11471cd764ad4972,
    0xd598e40d3dd89bcf,
    0x4aff1d108d4ec2c3,
    0xcedf722a585139ba,
    0xc2974eb4ee658828,
    0x733d226229feea32,
    0x806357d5a3f525f,
    0xca07c2dcb0cf26f7,
    0xfc89b393dd02f0b5,
    0xbbac2078d443ace2,
    0xd54b944b84aa4c0d,
    0xa9e795e65d4df11,
    0x4d4617b5ff4a16d5,
    0x504bced1bf8e4e45,
    0xe45ec2862f71e1d6,
    0x5d767327bb4e5a4c,
    0x3a6a07f8d510f86f,
    0x890489f70a55368b,
    0x2b45ac74ccea842e,
    0x3b0b8bc90012929d,
    0x9ce6ebb40173744,
    0xcc420a6a101d0515,
    0x9fa946824a12232d,
    0x47939822dc96abf9,
    0x59787e2b93bc56f7,
    0x57eb4edb3c55b65a,
    0xede622920b6b23f1,
    0xe95fab368e45eced,
    0x11dbcb0218ebb414,
    0xd652bdc29f26a119,
    0x4be76d3346f0495f,
    0x6f70a4400c562ddb,
    0xcb4ccd500f6bb952,
    0x7e2000a41346a7a7,
    0x8ed400668c0c28c8,
    0x728900802f0f32fa,
    0x4f2b40a03ad2ffb9,
    0xe2f610c84987bfa8,
    0xdd9ca7d2df4d7c9,
    0x91503d1c79720dbb,
    0x75a44c6397ce912a,
    0xc986afbe3ee11aba,
    0xfbe85badce996168,
    0xfae27299423fb9c3,
    0xdccd879fc967d41a,
    0x5400e987bbc1c920,
    0x290123e9aab23b68,
    0xf9a0b6720aaf6521,
    0xf808e40e8d5b3e69,
    0xb60b1d1230b20e04,
    0xb1c6f22b5e6f48c2,
    0x1e38aeb6360b1af3,
    0x25c6da63c38de1b0,
    0x579c487e5a38ad0e,
    0x2d835a9df0c6d851,
    0xf8e431456cf88e65,
    0x1b8e9ecb641b58ff,
    0xe272467e3d222f3f,
    0x5b0ed81dcc6abb0f,
    0x98e947129fc2b4e9,
    0x3f2398d747b36224,
    0x8eec7f0d19a03aad,
    0x1953cf68300424ac,
    0x5fa8c3423c052dd7,
    0x3792f412cb06794d,
    0xe2bbd88bbee40bd0,
    0x5b6aceaeae9d0ec4,
    0xf245825a5a445275,
    0xeed6e2f0f0d56712,
    0x55464dd69685606b,
    0xaa97e14c3c26b886,
    0xd53dd99f4b3066a8,
    0xe546a8038efe4029,
    0xde98520472bdd033,
    0x963e66858f6d4440,
    0xdde7001379a44aa8,
    0x5560c018580d5d52,
    0xaab8f01e6e10b4a6,
    0xcab3961304ca70e8,
    0x3d607b97c5fd0d22,
    0x8cb89a7db77c506a,
    0x77f3608e92adb242,
    0x55f038b237591ed3,
    0x6b6c46dec52f6688,
    0x2323ac4b3b3da015,
    0xabec975e0a0d081a,
    0x96e7bd358c904a21,
    0x7e50d64177da2e54,
    0xdde50bd1d5d0b9e9,
    0x955e4ec64b44e864,
    0xbd5af13bef0b113e,
    0xecb1ad8aeacdd58e,
    0x67de18eda5814af2,
    0x80eacf948770ced7,
    0xa1258379a94d028d,
    0x96ee45813a04330,
    0x8bca9d6e188853fc,
    0x775ea264cf55347d,
    0x95364afe032a819d,
    0x3a83ddbd83f52204,
    0xc4926a9672793542,
    0x75b7053c0f178293,
    0x5324c68b12dd6338,
    0xd3f6fc16ebca5e03,
    0x88f4bb1ca6bcf584,
    0x2b31e9e3d06c32e5,
    0x3aff322e62439fcf,
    0x9befeb9fad487c2,
    0x4c2ebe687989a9b3,
    0xf9d37014bf60a10,
    0x538484c19ef38c94,
    0x2865a5f206b06fb9,
    0xf93f87b7442e45d3,
    0xf78f69a51539d748,
    0xb573440e5a884d1b,
    0x31680a88f8953030,
    0xfdc20d2b36ba7c3d,
    0x3d32907604691b4c,
    0xa63f9a49c2c1b10f,
    0xfcf80dc33721d53,
    0xd3c36113404ea4a8,
    0x645a1cac083126e9,
    0x3d70a3d70a3d70a3,
    0xcccccccccccccccc,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x0,
    0x4000000000000000,
    0x5000000000000000,
    0xa400000000000000,
    0x4d00000000000000,
    0xf020000000000000,
    0x6c28000000000000,
    0xc732000000000000,
    0x3c7f400000000000,
    0x4b9f100000000000,
    0x1e86d40000000000,
    0x1314448000000000,
    0x17d955a000000000,
    0x5dcfab0800000000,
    0x5aa1cae500000000,
    0xf14a3d9e40000000,
    0x6d9ccd05d0000000,
    0xe4820023a2000000,
    0xdda2802c8a800000,
    0xd50b2037ad200000,
    0x4526f422cc340000,
    0x9670b12b7f410000,
    0x3c0cdd765f114000,
    0xa5880a69fb6ac800,
    0x8eea0d047a457a00,
    0x72a4904598d6d880,
    0x47a6da2b7f864750,
    0x999090b65f67d924,
    0xfff4b4e3f741cf6d,
    0xbff8f10e7a8921a4,
    0xaff72d52192b6a0d,
    0x9bf4f8a69f764490,
    0x2f236d04753d5b4,
    0x1d762422c946590,
    0x424d3ad2b7b97ef5,
    0xd2e0898765a7deb2,
    0x63cc55f49f88eb2f,
    0x3cbf6b71c76b25fb,
    0x8bef464e3945ef7a,
    0x97758bf0e3cbb5ac,
    0x3d52eeed1cbea317,
    0x4ca7aaa863ee4bdd,
    0x8fe8caa93e74ef6a,
    0xb3e2fd538e122b44,
    0x60dbbca87196b616,
    0xbc8955e946fe31cd,
    0x6babab6398bdbe41,
    0xc696963c7eed2dd1,
    0xfc1e1de5cf543ca2,
    0x3b25a55f43294bcb,
    0x49ef0eb713f39ebe,
    0x6e3569326c784337,
    0x49c2c37f07965404,
    0xdc33745ec97be906,
    0x69a028bb3ded71a3,
    0xc40832ea0d68ce0c,
    0xf50a3fa490c30190,
    0x792667c6da79e0fa,
    0x577001b891185938,
    0xed4c0226b55e6f86,
    0x544f8158315b05b4,
    0x696361ae3db1c721,
    0x3bc3a19cd1e38e9,
    0x4ab48a04065c723,
    0x62eb0d64283f9c76,
    0x3ba5d0bd324f8394,
    0xca8f44ec7ee36479,
    0x7e998b13cf4e1ecb,
    0x9e3fedd8c321a67e,
    0xc5cfe94ef3ea101e,
    0xbba1f1d158724a12,
    0x2a8a6e45ae8edc97,
    0xf52d09d71a3293bd,
    0x593c2626705f9c56,
    0x6f8b2fb00c77836c,
    0xb6dfb9c0f956447,
    0x4724bd4189bd5eac,
    0x58edec91ec2cb657,
    0x2f2967b66737e3ed,
    0xbd79e0d20082ee74,
    0xecd8590680a3aa11,
    0xe80e6f4820cc9495,
    0x3109058d147fdcdd,
    0xbd4b46f0599fd415,
    0x6c9e18ac7007c91a,
    0x3e2cf6bc604ddb0,
    0x84db8346b786151c,
    0xe612641865679a63,
    0x4fcb7e8f3f60c07e,
    0xe3be5e330f38f09d,
    0x5cadf5bfd3072cc5,
    0x73d9732fc7c8f7f6,
    0x2867e7fddcdd9afa,
    0xb281e1fd541501b8,
    0x1f225a7ca91a4226,
    0x3375788de9b06958,
    0x52d6b1641c83ae,
    0xc0678c5dbd23a49a,
    0xf840b7ba963646e0,
    0xb650e5a93bc3d898,
    0xa3e51f138ab4cebe,
    0xc66f336c36b10137,
    0xb80b0047445d4184,
    0xa60dc059157491e5,
    0x87c89837ad68db2f,
    0x29babe4598c311fb,
    0xf4296dd6fef3d67a,
    0x1899e4a65f58660c,
    0x5ec05dcff72e7f8f,
    0x76707543f4fa1f73,
    0x6a06494a791c53a8,
    0x487db9d17636892,
    0x45a9d2845d3c42b6,
    0xb8a2392ba45a9b2,
    0x8e6cac7768d7141e,
    0x3207d795430cd926,
    0x7f44e6bd49e807b8,
    0x5f16206c9c6209a6,
    0x36dba887c37a8c0f,
    0xc2494954da2c9789,
    0xf2db9baa10b7bd6c,
    0x6f92829494e5acc7,
    0xcb772339ba1f17f9,
    0xff2a760414536efb,
    0xfef5138519684aba,
    0x7eb258665fc25d69,
    0xef2f773ffbd97a61,
    0xaafb550ffacfd8fa,
    0x95ba2a53f983cf38,
    0xdd945a747bf26183,
    0x94f971119aeef9e4,
    0x7a37cd5601aab85d,
    0xac62e055c10ab33a,
    0x577b986b314d6009,
    0xed5a7e85fda0b80b,
    0x14588f13be847307,
    0x596eb2d8ae258fc8,
    0x6fca5f8ed9aef3bb,
    0x25de7bb9480d5854,
    0xaf561aa79a10ae6a,
    0x1b2ba1518094da04,
    0x90fb44d2f05d0842,
    0x353a1607ac744a53,
    0x42889b8997915ce8,
    0x69956135febada11,
    0x43fab9837e699095,
    0x94f967e45e03f4bb,
    0x1d1be0eebac278f5,
    0x6462d92a69731732,
    0x7d7b8f7503cfdcfe,
    0x5cda735244c3d43e,
    0x3a0888136afa64a7,
    0x88aaa1845b8fdd0,
    0x8aad549e57273d45,
    0x36ac54e2f678864b,
    0x84576a1bb416a7dd,
    0x656d44a2a11c51d5,
    0x9f644ae5a4b1b325,
    0x873d5d9f0dde1fee,
    0xa90cb506d155a7ea,
    0x9a7f12442d588f2,
    0xc11ed6d538aeb2f,
    0x8f1668c8a86da5fa,
    0xf96e017d694487bc,
    0x37c981dcc395a9ac,
    0x85bbe253f47b1417,
    0x93956d7478ccec8e,
    0x387ac8d1970027b2,
    0x6997b05fcc0319e,
    0x441fece3bdf81f03,
    0xd527e81cad7626c3,
    0x8a71e223d8d3b074,
    0xf6872d5667844e49,
    0xb428f8ac016561db,
    0xe13336d701beba52,
    0xecc0024661173473,
    0x27f002d7f95d0190,
    0x31ec038df7b441f4,
    0x7e67047175a15271,
    0xf0062c6e984d386,
    0x52c07b78a3e60868,
    0xa7709a56ccdf8a82,
    0x88a66076400bb691,
    0x6acff893d00ea435,
    0x583f6b8c4124d43,
    0xc3727a337a8b704a,
    0x744f18c0592e4c5c,
    0x1162def06f79df73,
    0x8addcb5645ac2ba8,
    0x6d953e2bd7173692,
    0xc8fa8db6ccdd0437,
    0x1d9c9892400a22a2,
    0x2503beb6d00cab4b,
    0x2e44ae64840fd61d,
    0x5ceaecfed289e5d2,
    0x7425a83e872c5f47,
    0xd12f124e28f77719,
    0x82bd6b70d99aaa6f,
    0x636cc64d1001550b,
    0x3c47f7e05401aa4e,
    0x65acfaec34810a71,
    0x7f1839a741a14d0d,
    0x1ede48111209a050,
    0x934aed0aab460432,
    0xf81da84d5617853f,
    0x36251260ab9d668e,
    0xc1d72b7c6b426019,
    0xb24cf65b8612f81f,
    0xdee033f26797b627,
    0x169840ef017da3b1,
    0x8e1f289560ee864e,
    0xf1a6f2bab92a27e2,
    0xae10af696774b1db,
    0xacca6da1e0a8ef29,
    0x17fd090a58d32af3,
    0xddfc4b4cef07f5b0,
    0x4abdaf101564f98e,
    0x9d6d1ad41abe37f1,
    0x84c86189216dc5ed,
    0x32fd3cf5b4e49bb4,
    0x3fbc8c33221dc2a1,
    0xfabaf3feaa5334a,
    0x29cb4d87f2a7400e,
    0x743e20e9ef511012,
    0x914da9246b255416,
    0x1ad089b6c2f7548e,
    0xa184ac2473b529b1,
    0xc9e5d72d90a2741e,
    0x7e2fa67c7a658892,
    0xddbb901b98feeab7,
    0x552a74227f3ea565,
    0xd53a88958f87275f,
    0x8a892abaf368f137,
    0x2d2b7569b0432d85,
    0x9c3b29620e29fc73,
    0x8349f3ba91b47b8f,
    0x241c70a936219a73,
    0xed238cd383aa0110,
    0xf4363804324a40aa,
    0xb143c6053edcd0d5,
    0xdd94b7868e94050a,
    0xca7cf2b4191c8326,
    0xfd1c2f611f63a3f0,
    0xbc633b39673c8cec,
    0xd5be0503e085d813,
    0x4b2d8644d8a74e18,
    0xddf8e7d60ed1219e,
    0xcabb90e5c942b503,
    0x3d6a751f3b936243,
    0xcc512670a783ad4,
    0x27fb2b80668b24c5,
    0xb1f9f660802dedf6,
    0x5e7873f8a0396973,
    0xdb0b487b6423e1e8,
    0x91ce1a9a3d2cda62,
    0x7641a140cc7810fb,
    0xa9e904c87fcb0a9d,
    0x546345fa9fbdcd44,
    0xa97c177947ad4095,
    0x49ed8eabcccc485d,
    0x5c68f256bfff5a74,
    0x73832eec6fff3111,
    0xc831fd53c5ff7eab,
    0xba3e7ca8b77f5e55,
    0x28ce1bd2e55f35eb,
    0x7980d163cf5b81b3,
    0xd7e105bcc332621f,
    0x8dd9472bf3fefaa7,
    0xb14f98f6f0feb951,
    0x6ed1bf9a569f33d3,
    0xa862f80ec4700c8,
    0xcd27bb612758c0fa,
    0x8038d51cb897789c,
    0xe0470a63e6bd56c3,
    0x1858ccfce06cac74,
    0xf37801e0c43ebc8,
    0xd30560258f54e6ba,
    0x47c6b82ef32a2069,
    0x4cdc331d57fa5441,
    0xe0133fe4adf8e952,
    0x58180fddd97723a6,
    0x570f09eaa7ea7648
];


text:
"
abcd_tline(l(w),c(w),r=w->0,g=w->0)

abcd_tline returns a function (of frequency and distance) which, itself, returns an abcd
matrix.

abcd_tline requires as arguments:

1. l(w) -> inductance per unit length as a function of frequency (rad/sec)
1. c(w) -> capacitance per unit length as a function of frequency (rad/sec)
1. r(w) -> series resistance per unit length as a function of frequency (rad/sec)
1. g(w) -> susceptance per unit length as a function of frequency (rad/sec)

abcd_tline returns a function f(ω,d) which, when supplied an angular frequency,
ω, and transmission line length, d, returns an ABCD matrix corresponding to this
transmission line.
"
function abcd_tline(l,c,r=w->0,g=w->0)
    return function abcd(w,d)
        Z0(w) = sqrt((r(w)+im*w*l(w))/(g(w)+im*w*c(w)))
        γ(w) = sqrt((r(w)+im*w*l(w))*(g(w)+im*w*c(w)))
        return [cosh(γ(w)*d) Z0(w)*sinh(γ(w)*d);
                 sinh(γ(w)*d)/Z0(w) cosh(γ(w)*d)]
    end
end

ext: jl


text:
abstract oneDoperator
include("partialDifferentiate/partialD.jl")

abstract chooseDimensionOperator

immutable byDimensionN{TypeOperator <: oneDoperator} <: chooseDimensionOperator
  n :: Int64
  functor :: TypeOperator
end

function resize_dims(n :: Int64, dims::Tuple{Vararg{Int64}})
  (reduce(*, dims[1:pf.n -1]), dims[pf.n], reduce(*, dims[pf.n+1:end]))
end

function eval(pf :: byDimensionN, data :: Array{Float64})
  dims_data = size(pf)
  dataResize = reshape(data, dims_data)
  res_anexample = eval(pf.functor, collect(dataResize[1,:,1]))
  res = fill(0.0 , size(dataResize,1), length(res_anexample), size(dataResize,3))
  for i1 = 1:size(res,1)
    for i3 = 1:size(res,3)
      res[i1,:,i3] = eval(pf.functor, collect(dataResize[i1,:,i3]))
    end
  end
  reshape(res, (dims_data[1:n-1]..., size(res,2), dims_data[n+1:end]...))
end


text:
function test_getters()
    model = MetNets.toy_model()
    @test MetNets.mets(model) == model.mets
    @test MetNets.rxns(model) == model.rxns

    for (i, rxn) in enumerate(MetNets.rxns(model))
        @test MetNets.ub(model, rxn) == model.ub[i]
        @test MetNets.ub(model, i) == model.ub[i]
        @test MetNets.lb(model, rxn) == model.lb[i]
        @test MetNets.lb(model, i) == model.lb[i]
        @test MetNets.rxns(model, i) == model.rxns[i]
        @test MetNets.rxns(model, rxn) == rxn
    end

    for (i, met) in enumerate(MetNets.mets(model))
        @test MetNets.mets(model, i) == model.mets[i]
        @test MetNets.mets(model, met) == model.mets[i]
    end

end
test_getters()


text:
#waits for one worker to begin before calling the next so that network is not slammed trying to send huge arrays to many workers simultaneously

function sequence_workers(wk_pool, func, args...)
    comms_chan = RemoteChannel(()->Channel{Integer}(length(wk_pool)))

    for worker in wk_pool
        remote_do(func, worker, args..., comms_chan)
        wait(comms_chan); report=take!(comms_chan)
        @assert worker==report
    end
end


text:
push!(LOAD_PATH, "./")

# if you want to experiment with the "pushed" matrix creation method,
# switch comments on the two following lines:
using Sparse23
#using Sparse23push
#using Sparse23raw
using Printf
using SparseArrays, LinearAlgebra

function dotest(dim::Int64, size::Int64)
    # 1) matrix creation:
    t1 = time_ns()
    if dim==2
        M, order, nc = PreLapl2(size)
        sizeV = size^2
    else
        M, order, nc = PreLapl3(size)
        sizeV = size^3
    end

    t1 = time_ns() - t1
    V = ones(Float64, sizeV)
    W = zeros(Float64, sizeV)

    # 2) matrix x vector product.
    t2 = time_ns()
    mul!(W, M, V)
    t2 = time_ns() - t2

    return order, nc, t1, t2
end

println()
@printf "%-5s %1s %-20s\n" "size" ":" "size of the grid."
@printf "%-5s %1s %-20s\n" "order" ":" "order of the matrix."
@printf "%-5s %1s %-20s\n" "nc" ":" "number of ceofficients."
@printf "%-5s %1s %-20s\n" "T.b" ":" "time to build the matrix (s)."
@printf "%-5s %1s %-20s\n" "T.p" ":" "time for matrix x vector product (s)."

# computation starts here:
for dim=2:3
    println("\n" * string(dim) * "d:")
    fw = open("RunningOn" * gethostname() * "-" * string(dim), "w")
    if dim==2
        size = 32
        sizeM = 2048
    else
        size = 16
        sizeM = 128
    end

    @printf "%4s %8s %9s %10s %10s \n" "size" "order" "nc" "T.b" "T.p"
    while size <= sizeM
        order, nc, t1, t2 = dotest(dim, size)   #warmup
        order, nc, t1, t2 = dotest(dim, size)

        s = @sprintf "%4d %8d %9d %10.2e %10.2e \n" size order nc float(t1)*10.0^(-9) float(t2)*10.0^(-9)
        print(s)
        write(fw, s)
        size *= 2
    end

    close(fw)
end

max_stars_repo_name: danielcort/benchmarks-python-julia-c

max_issues_repo_name: danielcort/benchmarks-python-julia-c

max_forks_repo_name: danielcort/benchmarks-python-julia-c


text:
abstract type IntegrationMethod end
struct NumericalIntegrationSimpson <: IntegrationMethod end
struct CustomSimpson <: IntegrationMethod end
struct BeyondLimber <: IntegrationMethod end

"""
    ComputeCℓ!(Cℓ::AbstractCℓ, WeightFunctionA::AbstractWeightFunction,
    WeightFunctionB::AbstractWeightFunction, BackgroundQuantities::BackgroundQuantities,
    ::AbstractCosmology, CosmologicalGrid::AbstractCosmologicalGrid,
    PowerSpectrum::AbstractPowerSpectrum, ::NumericalIntegrationSimpson)

This function evaluates the Angular Coefficients for all tomographic bins and
multipole values. In order to evaluate the numerical integrals, it is used the
Simpson numerical method from
[NumericalIntegration.jl](https://github.com/dextorious/NumericalIntegration.jl) .
"""
function  ComputeCℓ!(Cℓ::AbstractCℓ, WeightFunctionA::AbstractWeightFunction,
    WeightFunctionB::AbstractWeightFunction, BackgroundQuantities::BackgroundQuantities,
    ::AbstractCosmology, CosmologicalGrid::AbstractCosmologicalGrid,
    PowerSpectrum::AbstractPowerSpectrum, ::NumericalIntegrationSimpson)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    for idx_a in 1:size(Cℓ.CℓArray, 2)
        for idx_b in idx_a:size(Cℓ.CℓArray, 3)
            for idx_l in 1:size(Cℓ.CℓArray, 1)
                Integrand = c_0 .* WeightFunctionA.WeightFunctionArray[idx_a, :] .*
                WeightFunctionB.WeightFunctionArray[idx_b, :] ./
                (BackgroundQuantities.HZArray .* BackgroundQuantities.rZArray.^2) .*
                PowerSpectrum.InterpolatedPowerSpectrum[idx_l,:]
                Cℓ.CℓArray[idx_l, idx_a, idx_b] =
                NumericalIntegration.integrate(CosmologicalGrid.ZArray,
                Integrand, SimpsonEvenFast())
                Cℓ.CℓArray[idx_l, idx_b, idx_a] = Cℓ.CℓArray[idx_l, idx_a, idx_b]
            end
        end
    end
end

"""
    ComputeCℓ!(Cℓ::AbstractCℓ, WeightFunctionA::AbstractWeightFunction,
    WeightFunctionB::AbstractWeightFunction,
    BackgroundQuantities::AbstractBackgroundQuantities,
    ::AbstractCosmology, CosmologicalGrid::AbstractCosmologicalGrid,
    PowerSpectrum::AbstractPowerSpectrum, ::CustomSimpson)

This function evaluates the Angular Coefficients for all tomographic bins and
multipole values. In order to evaluate the numerical integrals, it has been
implemented the Simpson rule. The computation is accelerated by
[LoopVectorization.jl](https://github.com/JuliaSIMD/LoopVectorization.jl) .
"""
function  ComputeCℓ!(Cℓ::AbstractCℓ, WeightFunctionA::AbstractWeightFunction,
    WeightFunctionB::AbstractWeightFunction,
    BackgroundQuantities::AbstractBackgroundQuantities,
    ::AbstractCosmology, CosmologicalGrid::AbstractCosmologicalGrid,
    PowerSpectrum::AbstractPowerSpectrum, ::CustomSimpson)
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    check = true
    while check == true
        SimpsonWeights = SimpsonWeightArray(length(CosmologicalGrid.ZArray))
        ZStep = (last(CosmologicalGrid.ZArray)-first(CosmologicalGrid.ZArray)) /
        (length(CosmologicalGrid.ZArray)-1)
        Cℓ.CℓArray = CustomCℓIntegrator!(SimpsonWeights, Cℓ.CℓArray,
        WeightFunctionA.WeightFunctionArray, WeightFunctionB.WeightFunctionArray,
        ZStep, BackgroundQuantities.HZArray, BackgroundQuantities.rZArray,
        PowerSpectrum.InterpolatedPowerSpectrum)
        if any(isnan,Cℓ.CℓArray)

        else
            check = false
        end
    end
end

function CustomCℓIntegrator!(SimpsonWeights::Array{Float64}, CℓArray::Array{Float64, 3},
    WArrayA::Array{Float64, 2}, WArrayB::Array{Float64, 2}, ZStep::Float64,
    HZArray::Array{Float64}, rZArray::Array{Float64}, InterpolatedPmm::Array{Float64, 2})
    c_0 = 2.99792458e5 #TODO: find a package containing the exact value of
                       #physical constants involved in calculations
    @avx for i ∈ axes(CℓArray,2),
        j ∈ axes(CℓArray,3),
        l ∈ axes(CℓArray,1)
        for z ∈ axes(rZArray,1)
            CℓArray[l,i,j] += c_0 * WArrayA[i, z] * WArrayB[j, z] /
            (HZArray[z] * rZArray[z]^2) * InterpolatedPmm[l,z] * SimpsonWeights[z]
        end
    end
    CℓArray .*= ZStep
end

function  ComputeCℓ!(Cℓ::AbstractCℓ, TransferFunctionA::AbstractTransferFunction,
    TransferFunctionB::AbstractTransferFunction, ::AbstractCosmology,
    CosmologicalGrid::AbstractCosmologicalGrid, PowerSpectrum::AbstractPowerSpectrum,
    ::BeyondLimber)
    Integrand = zeros(size(Cℓ.CℓArray))
    WeightsMatrix = UnevenTrapzWeightMatrix(CosmologicalGrid.KBeyondLimberArray)
    @avx for i ∈ axes(Cℓ.CℓArray,2),
        j ∈ axes(Cℓ.CℓArray,3),
        l ∈ axes(Cℓ.CℓArray,1)
        for k ∈ axes(CosmologicalGrid.KBeyondLimberArray,2)
            Integrand[l,i,j] += 2 / π * WeightsMatrix[l, k] *
            TransferFunctionA.TransferFunctionArray[i, l, k] *
            TransferFunctionB.TransferFunctionArray[j, l, k] *
            PowerSpectrum.InterpolatedPowerSpectrumBeyondLimber[l,k] *
            CosmologicalGrid.KBeyondLimberArray[l, k]^2
        end
    end
    Cℓ.CℓArray = Integrand
end


text:
include("utils.jl")

@testset "Positive Numbers" begin
    M = PositiveNumbers()
    @testset "Positive Numbers Basics" begin
        @test repr(M) == "PositiveNumbers()"
        @test repr(PositiveVectors(2)) == "PositiveVectors(2)"
        @test repr(PositiveMatrices(2, 3)) == "PositiveMatrices(2, 3)"
        @test repr(PositiveArrays(2, 3, 4)) == "PositiveArrays(2, 3, 4)"
        @test representation_size(M) == ()
        @test manifold_dimension(M) == 1
        @test !is_manifold_point(M, -1.0)
        @test_throws DomainError is_manifold_point(M, -1.0, true)
        @test is_tangent_vector(M, 1.0, 0.0; check_base_point = false)
        @test flat(M, 1.0, FVector(TangentSpace, 1.0)) == FVector(CotangentSpace, 1.0)
        @test sharp(M, 1.0, FVector(CotangentSpace, 1.0)) == FVector(TangentSpace, 1.0)
        @test vector_transport_to(M, 1.0, 3.0, 2.0, ParallelTransport()) == 6.0
        @test retract(M, 1.0, 1.0) == exp(M, 1.0, 1.0)
        @test isinf(injectivity_radius(M))
        @test isinf(injectivity_radius(M, -2.0))
        @test isinf(injectivity_radius(M, -2.0, ExponentialRetraction()))
        @test isinf(injectivity_radius(M, ExponentialRetraction()))
        @test project(M, 1.5, 1.0) == 1.0
        @test zero_tangent_vector(M, 1.0) == 0.0
        X = similar([1.0])
        zero_tangent_vector!(M, X, 1.0)
        @test X == [0.0]
    end
    types = [Float64]
    TEST_FLOAT32 && push!(types, Float32)
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [1.0, 4.0, 2.0])
            test_manifold(
                M,
                pts,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vector_spaces = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_default_vector_transport = true,
                test_vee_hat = false,
                is_mutating = false,
            )
        end
    end
    @testset "Power of Positive Numbers" begin
        M2 = PositiveVectors(2)
        for T in types
            pts2 = [convert.(Ref(T), v) for v in [[1.0, 1.1], [3.0, 3.3], [2.0, 2.2]]]
            test_manifold(
                M2,
                pts2,
                test_forward_diff = false,
                test_reverse_diff = false,
                test_vector_spaces = false,
                test_project_tangent = true,
                test_musical_isomorphisms = true,
                test_default_vector_transport = true,
                test_vee_hat = false,
            )
        end
    end
end


text:
using MyFirstPkg
using Test

@testset "MyFirstPkg.jl" begin
    @test my_f(3) == 9
end


text:
import Pkg
Pkg.instantiate() # downloads all dependencies for the current project
Pkg.add("JSON")

using JSON

###################
### Write data ####
###################

# dictionary to write
dict1 = Dict("param1" => 1, "param2" => 2,
            "dict" => Dict("d1"=>1.,"d2"=>1.,"d3"=>1.))

# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(dict1)

# write the file with the stringdata variable information
open("write_read.json", "w") do f
        write(f, stringdata)
     end

###################
### Read data #####
###################

# create variable to write the information
dict2 = Dict()
open("write_read.json", "r") do f
    global dict2
    dicttxt = read(f,String)  # file information to string
    dict2=JSON.parse(dicttxt)  # parse and transform data
end

# print both dictionaries
println(dict1)
println(dict2)


text:
module ObjectPools

# This is inspired by
#      https://github.com/tpapp/ObjectPools.jl
# but simplified and now evolved

using DataStructures: Stack

export acquire!, release!
export VectorPool, MatrixPool, ArrayPool
using Base.Threads: nthreads, threadid

# TODO:
# * General ObjectPool
# * convert to more general ArrayPool
# * Consider allowing the most common "derived" types and enabling
#   those via dispatch; e.g. Duals? Complex? All possible Floats? ....

struct VectorPool{T}
    #        tid    stack    object
    arrays::Vector{Stack{Vector{T}}}
    VectorPool{T}() where {T} = new( [ Stack{Vector{T}}() for _=1:nthreads() ] )
end

acquire!(pool::VectorPool{T}, sz::Union{Integer, NTuple{N}}, ::Type{T}) where {T, N} =
        acquire!(pool, sz)

acquire!(pool::VectorPool{T}, len::Integer, ::Type{T}) where {T} =
        acquire!(pool, len)

function acquire!(pool::VectorPool{T}, len::Integer) where {T}
    tid = threadid()
    if !isempty(pool.arrays[tid]) > 0
        x = pop!(pool.arrays[tid])
        if len != length(x)
            resize!(x, len)
        end
        return x
    else
        return Vector{T}(undef, len)
    end
end

function release!(pool::VectorPool{T}, x::Vector{T}) where {T}
    tid = threadid()
    push!(pool.arrays[tid], x)
    return nothing
end

# Vector -> Array  -> Vector

function acquire!(pool::VectorPool{T}, sz::NTuple{N}) where {T, N}
    tid = threadid()
    len = prod(sz)::Integer
    if !isempty(pool.arrays[tid])
        x = pop!(pool.arrays[tid])
        if len != length(x)
            resize!(x, len)
        end
        return reshape(x, sz)
    else
        return Array{T, N}(undef, sz)
    end
end

release!(pool::VectorPool{T}, x::Array{T}) where {T} =
        release!(pool, reshape(x, :))

# fallbacks that allocate and relase to the main GC

acquire!(pool::VectorPool{T}, len::Integer, S::Type{T1}) where {T, T1} =
        Vector{S}(undef, len)

acquire!(pool::VectorPool{T}, sz::NTuple{N}, S::Type{T1}) where {T, N, T1} =
        Array{T, N}(undef, sz)

release!(pool::VectorPool{T}, x::AbstractVector) where {T} =
    nothing

end # module


text:
export FourRooms

mutable struct FourRooms <: AbstractGridWorld
    world::GridWorldBase{Tuple{Empty,Wall,Goal}}
    agent_pos::CartesianIndex{2}
    agent::Agent
end

function FourRooms(;n=19, agent_start_pos=CartesianIndex(2,2))
    objects = (EMPTY, WALL, GOAL)
    world = GridWorldBase(objects, n,n)
    world[EMPTY, 2:n-1, 2:n-1] .= true
    world[WALL, [1,n], 1:n] .= true
    world[WALL, 1:n, [1,n]] .= true
    world[WALL, ceil(Int,n/2), vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n)] .= true
    world[EMPTY, ceil(Int,n/2), vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n)] .= false
    world[WALL, vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n), ceil(Int,n/2)] .= true
    world[EMPTY, vcat(2:ceil(Int,n/4)-1,ceil(Int,n/4)+1:ceil(Int,n/2)-1,ceil(Int,n/2):ceil(Int,3*n/4)-1,ceil(Int,3*n/4)+1:n), ceil(Int,n/2)] .= false
    world[GOAL, n-1, n-1] = true
    world[EMPTY, n-1, n-1] = false
    FourRooms(world,agent_start_pos,Agent(dir=RIGHT))
end

function (w::FourRooms)(::MoveForward)
    dir = get_dir(w.agent)
    dest = dir(w.agent_pos)
    if !w.world[WALL, dest]
        w.agent_pos = dest
    end
    w
end


text:
# DIRS
const MODEL_RAW_DIR = joinpath(RAW_DATA_DIR, PROJ_IDER)
const MODEL_PROCESSED_DATA_DIR = joinpath(PROCESSED_DATA_DIR, basename(MODEL_RAW_DIR))
const MODEL_FIGURES_DIR = joinpath(FIGURES_DATA_DIR, basename(MODEL_RAW_DIR))
const MODEL_CACHE_DIR = joinpath(MODEL_PROCESSED_DATA_DIR, "cache")

function _create_dirs()
    for dir in [MODEL_PROCESSED_DATA_DIR,
                MODEL_FIGURES_DIR,
                MODEL_CACHE_DIR
            ]
        isdir(dir) && continue
        try; mkpath(dir); catch end
    end
end

# FILES
const MODEL_RAW_MAT_FILE = joinpath(MODEL_RAW_DIR, "e_coli_core.mat")
const BASE_MODELS_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "load_model.bson")
const EXCH_MET_MAP_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "exch_met_map.bson")
const MAXENT_RES_FILE = joinpath(MODEL_PROCESSED_DATA_DIR, "maxent_fba_ep_res.bson")


text:
using Printf
# using Plots
using Statistics
using LinearAlgebra
using BenchmarkTools

using PyPlot
pygui(true);

dim=2
mult=3
n=30

i=5

tau=2*pi
omega=1

#Delayed Mathieu parameters
kappa=0.1
delta=3.2
epsilon=1
b=0.5

#Cauchy matrices
function A(t)
 [-kappa -(delta+epsilon*cos(omega*t)); 1 0]
end

B=[0 b; 0 0]

deltat=tau/n
T=2*pi/omega

function checkT(tau,T)
    if tau <= T
        println("OK")
    else
        println("Error")
    end
end

checkT(tau,T)

#initial function
function v_ini(t)
    transpose([1*sin(omega*(t+tau/2)), 2*sin(omega*(t+tau/3)), 3*sin(omega*(t+tau/4)), 4*sin(omega*(t+tau/5)), 5*sin(omega*(t+tau/6)), 6*sin(omega*(t+tau/7))])
end

vect=zeros(n+1,mult*dim)
for j = 0:n
        vect[j+1,:]=v_ini(-tau+j*deltat)
end
vect
function simulation(init)
    k=Int16(floor(i*T/deltat))
    solmatrix=zeros(k,mult*dim)
    sol=vcat(init,solmatrix)

    for m = 1:dim:(dim*mult-(dim-1))
        for j = 0:(k-1)
            Y1=sol[n+1+j,m:m+(dim-1)]
            Y2=Y1+deltat/2*(A(j*deltat)*Y1+B*sol[j+1,m:m+(dim-1)])
            Y3=Y1+deltat/2*(A((j+1/2)*deltat)*Y2+B*(sol[j+1,m:m+(dim-1)]+sol[j+2,m:m+(dim-1)])/2)
            Y4=Y1+deltat*(A((j+1/2)*deltat)*Y3+B*(sol[j+1,m:m+(dim-1)]+sol[j+2,m:m+(dim-1)])/2)
            Y=Y1+deltat/6*((A(j*deltat)*Y1+B*sol[j+1,m:m+(dim-1)])+2*(A((j+1/2)*deltat)*Y2+B*(sol[j+1,m:m+(dim-1)]+sol[j+2,m:m+(dim-1)])/2)+2*(A((j+1/2)*deltat)*Y3+B*(sol[j+1,m:m+(dim-1)]+sol[j+2,m:m+(dim-1)])/2)+(A((j+1)*deltat)*Y4+B*sol[j+2,m:m+(dim-1)]))
            sol[n+2+j,m:m+(dim-1)]=transpose(Y)
        end
    end
    tvect=transpose(zeros(1,n+1+k))
    for j = 1:(n+1+k)
        tvect[j]=-n*deltat+j*deltat
    end
    solt=hcat(tvect,sol)
    return(solt)
end

solt1=simulation(vect)
@btime simulation(vect)
@time simulation(vect)

tplot=[solt1[j,1] for j in 1:size(solt1,1)]
v1plot=[solt1[j,2] for j in 1:size(solt1,1)]
v2plot=[solt1[j,3] for j in 1:size(solt1,1)]
v3plot=[solt1[j,4] for j in 1:size(solt1,1)]
v4plot=[solt1[j,5] for j in 1:size(solt1,1)]
v5plot=[solt1[j,6] for j in 1:size(solt1,1)]
v6plot=[solt1[j,7] for j in 1:size(solt1,1)]
plot1=plot(tplot,hcat(v1plot,v2plot,v3plot,v4plot,v5plot,v6plot))

function Ha(solt)
    rat=Int16(floor(T/tau))
    ha=zeros(n,(i+1)*mult*dim)
    for j = 1:i+1
        ha[:,(j-1)*mult*dim+1:j*mult*dim]=solt[(j-1)*n+1:j*n,2:mult*dim+1]
    end
    return(ha)
end

ha=Ha(solt1)

function eigenv(haa)
    eigena=zeros(mult*dim,1)
    Hglob=zeros(mult*dim,i*mult*dim)
    for j = 1:i
        Hglob[:,(j-1)*mult*dim+1:j*mult*dim]=inv(transpose(haa[:,(j-1)*mult*dim+1:j*mult*dim])*haa[:,(j-1)*mult*dim+1:j*mult*dim])*(transpose(haa[:,(j-1)*mult*dim+1:j*mult*dim])*haa[:,j*mult*dim+1:(j+1)*mult*dim])
    end
    for j = 1:i
        eigena=hcat(eigena,eigvals(Hglob[:,(j-1)*mult*dim+1:j*mult*dim]))
    end
    eigena=eigena[:,2:i+1]
    return(eigena)
end

U=eigenv(ha)

eigRe1=[real(U[n,1]) for n in 1:size(U,1)]
eigIm1=[imag(U[n,1]) for n in 1:size(U,1)]
plot1=plot(eigRe1,eigIm1,seriestype=:scatter)
Re_circle=[cos(n) for n in 0:0.05:2*pi]
Im_circle=[sin(n) for n in 0:0.05:2*pi]
plot!(plot1,Re_circle,Im_circle)

eigRe2=[real(U[n,2]) for n in 1:size(U,1)]
eigIm2=[imag(U[n,2]) for n in 1:size(U,1)]
plot2=plot(eigRe2,eigIm2,seriestype=:scatter)
plot!(plot2,Re_circle,Im_circle)

eigRe3=[real(U[n,3]) for n in 1:size(U,1)]
eigIm3=[imag(U[n,3]) for n in 1:size(U,1)]
plot3=plot(eigRe3,eigIm3,seriestype=:scatter)
plot!(plot3,Re_circle,Im_circle)

eigRe4=[real(U[n,4]) for n in 1:size(U,1)]
eigIm4=[imag(U[n,4]) for n in 1:size(U,1)]
plot4=plot(eigRe4,eigIm4,seriestype=:scatter)
plot!(plot4,Re_circle,Im_circle)

eigRe5=[real(U[n,5]) for n in 1:size(U,1)]
eigIm5=[imag(U[n,5]) for n in 1:size(U,1)]
plot5=plot(eigRe5,eigIm5,seriestype=:scatter)
plot!(plot5,Re_circle,Im_circle)


text:
module PkgUtils

# deps
using Pkg, Pkg.TOML, Statistics, LightGraphs, MetaGraphs
using SnoopCompile

# file includes
include("dependencies.jl")
# include("treeshaker.jl")
include("snapshot.jl")
include("environments.jl")

# exports
export get_dependents, get_dependencies, get_latest_version
export snapshot!, undo!
export @project!_str, @manifest!_str

end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    SLIC(k, m; tol=1e-4, maxiter=10)

A method for clustering geospatial data into approximately `k`
clusters using Simple Linear Iterative Clustering (SLIC).
The method produces clusters of samples that are spatially
connected based on a distance `dₛ` and that, at the same
time, are similar in terms of `vars` with distance `dᵥ`.
The tradeoff is controlled with a hyperparameter parameter
`m` in an additive model `dₜ = √(dᵥ² + m²(dₛ/s)²)`.

## Parameters

* `k`       - Approximate number of clusters
* `m`       - Hyperparameter of SLIC model
* `tol`     - Tolerance of k-means algorithm (default to `1e-4`)
* `maxiter` - Maximum number of iterations (default to `10`)
* `weights` - Dictionary with weights for each attribute (default to `nothing`)

## References

* Achanta et al. 2011. [SLIC superpixels compared to state-of-the-art
  superpixel methods](https://ieeexplore.ieee.org/document/6205760)
"""
struct SLIC{W} <: ClusteringMethod
  k::Int
  m::Float64
  tol::Float64
  maxiter::Int
  weights::W
end

function SLIC(k::Int, m::Real; tol=1e-4, maxiter=10, weights=nothing)
  @assert tol > 0 "invalid tolerance"
  @assert maxiter > 0 "invalid number of iterations"
  SLIC{typeof(weights)}(k, m, tol, maxiter, weights)
end

function partition(data, method::SLIC)
  # normalize atributes
  𝒯 = TableDistances.normalize(values(data))
  Ω = georef(first(𝒯), domain(data))

  # weights for each attribute
  weights = method.weights

  # SLIC hyperparameter
  m = method.m

  # initial spacing of clusters
  s = slic_spacing(Ω, method)

  # initialize cluster centers
  c = slic_initialization(Ω, s)

  # ball neighborhood search
  searcher = BallSearch(Ω, MetricBall(maximum(s)))

  # pre-allocate memory for label and distance
  l = fill(0, nelements(Ω))
  d = fill(Inf, nelements(Ω))

  # performance parameters
  tol     = method.tol
  maxiter = method.maxiter

  # Lloyd's (a.k.a. k-means) algorithm
  err, iter = Inf, 0
  while err > tol && iter < maxiter
    o = copy(c)

    slic_assignment!(Ω, searcher, weights, m, s, c, l, d)
    slic_update!(Ω, c, l)

    err = norm(c - o) / norm(o)
    iter += 1
  end

  orphans = findall(iszero, l)
  if length(orphans) > 0
    assigned = findall(!iszero, l)
    Ω₀ = view(domain(Ω), assigned)
    csearcher = KNearestSearch(Ω₀, 1)

    for orphan in orphans
      p = centroid(Ω, orphan)
      i = search(p, csearcher)[1]
      l[orphan] = l[assigned[i]]
    end
  end

  subsets = [findall(isequal(k), l) for k in 1:length(c)]

  Partition(data, subsets)
end

slic_spacing(data, method) = slic_srecursion(method.k, sides(boundingbox(data)))

# given the desired number of clusters and the sides of the bounding box
# of the domain, returns the spacing for each dimension recursively
function slic_srecursion(k, l)
  d = length(l)

  # base case
  d == 1 && return [l[1] / k]

  # compute the spacing for the j-th dimension
  j  = argmax(l)
  kⱼ = ceil(Int, k^(1/d))
  sⱼ = l[j]/kⱼ

  # update the new k and l
  kₙ = ceil(Int, k/kⱼ)
  lₙ = [l[begin:j-1]; l[j+1:end]]

  # then recursively compute the spacing for the remaining dimensions
  s  = slic_srecursion(kₙ, lₙ)

  [s[begin:j-1]; [sⱼ]; s[j:end]]
end

function slic_initialization(data, s)
  # efficient neighbor search
  searcher = KNearestSearch(data, 1)

  # bounding box properties
  bbox = boundingbox(data)
  lo, up = coordinates.(extrema(bbox))

  # cluster centers
  clusters = Vector{Int}()
  neighbor = Vector{Int}(undef, 1)
  ranges = [(l+sᵢ/2):sᵢ:u for (l, sᵢ, u) in zip(lo, s, up)]
  for x in Iterators.product(ranges...)
    search!(neighbor, Point(x), searcher)
    push!(clusters, neighbor[1])
  end

  unique(clusters)
end

function slic_assignment!(data, searcher, weights, m, s, c, l, d)
  sₘ = maximum(s)
  for (k, cₖ) in enumerate(c)
    pₖ = centroid(data, cₖ)
    inds = search(pₖ, searcher)

    # distance between points
    X  = (coordinates(centroid(data, ind)) for ind in inds)
    xₖ = [coordinates(pₖ)]
    dₛ = pairwise(Euclidean(), X, xₖ)

    # distance between variables
    𝒮ᵢ = view(data, inds)
    𝒮ₖ = view(data, [cₖ])
    V  = values(𝒮ᵢ)
    vₖ = values(𝒮ₖ)
    td = TableDistance(normalize=false, weights=weights)
    dᵥ = pairwise(td, V, vₖ)

    # total distance
    dₜ = @. √(dᵥ^2 + m^2 * (dₛ/sₘ)^2)

    @inbounds for (i, ind) in enumerate(inds)
      if dₜ[i] < d[ind]
        d[ind] = dₜ[i]
        l[ind] = k
      end
    end
  end
end

function slic_update!(data, c, l)
  for k in 1:length(c)
    inds = findall(isequal(k), l)
    X  = (coordinates(centroid(data, ind)) for ind in inds)
    μ  = [mean(X)]
    dₛ = pairwise(Euclidean(), X, μ)
    @inbounds c[k] = inds[argmin(vec(dₛ))]
  end
end


text:
function seq_along(x)
    1:length(x[:])
end

function seq(from, to, by = 1)
    from:by:to
end

function Range(x::AbstractArray{<:Real})
    minimum(x), maximum(x)
end

function aggregate(x::AbstractArray{T, 1}, by::AbstractArray{T2, 1}, FUN::Function = mean) where {
    T<:Real, T2<:Real }

    grps = unique(by)
    # vals = ones(length(grps));
    vals = []
    @views for grp in grps
        val = FUN(x[by .== grp])
        push!(vals, val)
    end
    Dict(grps .=> vals)
end

cbind  = hcat
rbind  = vcat
# cbind! = hcat!
# rbind! = vcat!
# hcat! not defined

export cbind, rbind,
    seq_along, seq, Range, aggregate


text:
export metadata

"""
Metadata contains information about a SAS data file.

*Fields*
- `filename`: file name/path of the SAS data set
- `encoding`: file encoding e.g. "ISO8859-1"
- `endianness`: either `:LittleEndian`` or `:BigEndian`
- `compression`: could be `:RLE`, `:RDC`, or `:none`
- `pagesize`: size of each data page in bytes
- `npages`: number of pages in the file
- `nrows`: number of data rows in the file
- `ncols`: number of data columns in the file
- `columnsinfo`: vector of column symbols and their respective types
"""
struct Metadata
    filename::AbstractString
    encoding::AbstractString          # e.g. "ISO8859-1"
    endianness::Symbol                # :LittleEndian, :BigEndian
    compression::Symbol               # :RDC, :RLE
    pagesize::Int
    npages::Int
    nrows::Int
    ncols::Int
    columnsinfo::Vector{Pair{Symbol, Type}}
end

"""
Return metadata of a SAS data file.  See `SASLib.Metadata`.
"""
function metadata(fname::AbstractString)
    h = nothing
    try
        h = SASLib.open(fname, verbose_level = 0)
        rs = read(h, 1)
        _metadata(h, rs)
    finally
        h !== nothing && SASLib.close(h)
    end
end

# construct Metadata struct using handler & result set data
function _metadata(h::Handler, rs::ResultSet)
    cmp = h.compression == compression_method_rle ?
            :RLE : (h.compression == compression_method_rdc ? :RDC : :none)
    Metadata(
        h.config.filename,
        h.file_encoding,
        h.file_endianness,
        cmp,
        h.page_length,
        h.page_count,
        h.row_count,
        h.column_count,
        [Pair(x, eltype(rs[x])) for x in names(rs)]
    )
end

# pretty print
function Base.show(io::IO, md::Metadata)
    println(io, "File: ", md.filename, " (", md.nrows, " x ", md.ncols, ")")
    displaytable(io, colfmt(md); index=true)
end

# Column display format
function colfmt(md::Metadata)
    [string(first(p), "(", typesfmt(typesof(last(p))),
        ")") for p in md.columnsinfo]
end

# Compact types format
# e.g. (Date, Missing) => "Date/Missing"
function typesfmt(ty::Tuple; excludemissing = false)
    ar = excludemissing ?
        collect(Iterators.filter(x -> x != Missing, ty)) : [ty...]
    join(string.(ar), "/")
end

# Extract types from a Union
# e.g.
#   Union{Int64, Int32, Float32} => (Int64, Int32, Float32)
#   Int32 => (Int32,)
function typesof(ty::Type)
    if ty isa Union
        y = typesof(ty.b)
        y isa Tuple ? (ty.a, y...) : (ty.a, y)
    else
        (ty,)
    end
end


text:
export conv_bias_act, conv_bias_act!

function conv_bias_act(x::AbstractArray{xT,N}, w::AbstractArray{wT,N},
                cdims::ConvDims, b::AbstractArray{bT,N}, σ=identity; kwargs...) where {xT, wT, bT, N}
    y = similar(x, promote_type(xT, wT, bT), output_size(cdims)..., channels_out(cdims), size(x,N))
    conv_bias_act!(y, x, w, cdims, b, σ; kwargs...)
    return y
end

function conv_bias_act!(y::AbstractArray{yT,5}, x::AbstractArray{xT,5}, w::AbstractArray{wT,5},
                cdims::ConvDims, b::AbstractArray{bT,5}, σ=identity; kwargs...) where {yT, xT, wT, bT}
    conv!(y, x, w, cdims)
    y .= σ.(y .+ b)
    return y
end

for N in (3, 4)
    @eval begin
        function $(Symbol("conv_bias_act!"))(
                        y::AbstractArray{yT,$N}, x::AbstractArray{xT,$N},
                        w::AbstractArray{wT,$N}, cdims::ConvDims,
                        b::AbstractArray{bT,$N}, σ=identity;
                        kwargs...) where {yT, xT, wT, bT}
            $(Symbol("conv_bias_act!"))(
                insert_singleton_spatial_dimension(y, $(5 - N)),
                insert_singleton_spatial_dimension(x, $(5 - N)),
                insert_singleton_spatial_dimension(w, $(5 - N)),
                insert_singleton_spatial_dimension(cdims, $(5 - N)),
                insert_singleton_spatial_dimension(b, $(5 - N)),
                σ;
                kwargs...
            )

            # We explicitly return `y` here, because the backend call
            # itself may return a reshaped view, which we don't want.
            return y
        end
    end
end


text:
@info "Loading packages"
using CPUTime
using Plots
using Printf

include("src/grids.jl")
include("src/util.jl")

## Parameters / Configuration
# Size Distribtion
r̅ = (10.0) * 1e-6 # mode radius of initial distribution, input micron convert to m
x̅ = mass_from_r(r̅) # mean initial droplet mass (kg)
L = (1.0)*1e-3 # total water content, g/m3 convert to kg/m3

# Grid Spacing
n = 4
α = 2^(1/n) # Mass bin scaling ratio, 2^(1/n) where 'n' is the number of bins
            # between mass doubling

m = ceil(Integer, 1 + 3*log(rm / r₁)/log(α)) # number of mass bins to populate
Δy = Δlnr = log(α) / 3  # constant grid distance of logarithmic grid

# Collision kernel
# Options:
#   1) :golovin - Golovin (1963)
#   2) :hydro - basic hydrodynamic kernel with unit coal/coll efficiency
#   3) :long - Long collision efficiency
kernel = :long

# Time Integration
tmax = 3601 # seconds
Δt = 5.0 # s
Δt_plot = 10 # minutes
nt = ceil(Integer, tmax / Δt)

# Other configs
debug = false
do_plots = true

# Report configuration to terminal
println("""
1D Stochastic Collection Equation Solver
========================================

Initial Conditions
------------------
L = $(L*1e3) g/m³ - Total LWC
r̅ = $(r̅*1e6) μm - Mode radius
x̅ = $x̅ kg - Mean droplet mass

Grid Setup
----------
$m bins spanning radii ($(r₁*1e6) μm - $(rm*1e6) μm)
α factor = $α (mass doubling every $n bins)

Time Integration
----------------
Integrating to $tmax seconds by Δt=$Δt

- DEBUG is $debug
- DO_PLOTS is $do_plots

""")

## Model Setup

# Mass and Radius Grid
#=
I find it a bit easier to reason through things in droplet-radius-space. We can
adopt the same convention of constructing a mass grid with an α scaling factor
(e.g. for α=√2 every mass bin is twice as big as the bin twice before it) and
write simple, similar recurrence relationship in radius space:

  rᵢ = α rᵢ₋₁

Given r₀ as the smallest bin, we can then solve that

  rᵢ = r₁ α^{(i-1)/3} s.t. i ∈ ℤ > 0

Computing the number of grid cells...
  ln rm > ln(r1) + ln(α) (m-1)/3
  ln rm - ln r1 > ln α (m-1)/3
  3 ln(rm / r1) / (ln α) > m - 1

=#
rᵢ = r₁*(α.^((collect(1:m) .- 1)./3)) # meter
xᵢ = mass_from_r.(rᵢ) # kg

# Initial droplet distribution; g(y, t) = 3x² * n(x, t), eq. 2 from B98
gᵢ = (3 .* xᵢ.^2) .* nc.(xᵢ, L=L, x̅=x̅) # kg / m3

# Courant Numbers
#=
The final take-home here is that these courant limits are purely a function of
the mass grid discretization. On the one hand, they can be pre-computed and then
cached for retrieval as needed. This fact is used in combination with the fact
that the same holds true for the collision kernel to pre-compute all the limits
in the original code, as well as pre-estimate where collisions will land in the
actual coad subroutine
=#

println("""
COMPUTE COURANT LIMITS ON GRID""")
CPUtic()
c, ima = courant(xᵢ)
elapsed = CPUtoq()
@printf "%3.1f seconds elapsed\n" elapsed

# Collision Kernel

println("""
COMPUTE COLLISION KERNELS ON GRID""")
CPUtic()
ck = kernels(xᵢ, kernel)
elapsed = CPUtoq()
@printf "%3.1f seconds elapsed\n" elapsed

## Plot Initial Check

# Plot initial conditions and then begin the time loop
if do_plots
  println("Plotting initial conditions")
  p = plot(
    rᵢ*1e6, gᵢ*1e3, label="t = 0 min",
    xaxis=:log, xlabel="r (μm)",
    # xlim=(0.5, 5000), xticks=[1, 10, 100, 1000],
    # ylim=(0, 0.9), yticks=0:0.06:0.9,
    ylabel="g (g / m³)",
  )
  display(p)
end

## TIME LOOP

#=
All of the time looping logic here is super old-school and can be totally
re-written for simplicity.
=#
# time integration
tlmin = 1e-6
t = 0.0
lmin = 0.0
total_runtime = 0.0

@inline function coad!(g, x, c, ck, ima, Δt, Δy; debug=false, gmin=gmin)

  # Lower and Upper integration limit i0, i1
  i0, i1 = find_bounds(g, gmin=gmin)

  if debug
    @printf "DEBUG bnds_check %6d %8d %8d\n" t i0 i1
  end

  # Main collision/coalescence loop
  for i ∈ i0:i1
    for j ∈ i:i1
      k = ima[i, j]  # Get pre-computed index of coalescence bin edge
      kp = k + 1

      # PORT - handle a weird edge condition in the initialization of ima?
      if k < 1
        continue
      end

      # We did not pre-scale the collison kernels LUT by Δt and Δy in this
      # version so we have to account for them here.
      x0 = ck[i, j] * g[i] * g[j] * Δt * Δy
      x0 = min(x0, g[i] * x[j])
      if j != k # Not sure what's going on here.
        x0 = min(x0, g[j] * x[i])
      end

      gsi = x0 / x[j]
      gsj = x0 / x[i]
      gsk = gsi + gsj
      g[i] = g[i] - gsi
      g[j] = g[j] - gsj
      gk = g[k] + gsk

      if gk > gmin

        # ORIGINAL
        # x1 = log(gᵢ[kp] / gk + 1e-60)

        # MODIFIED - Apply a limiter to avoid negative args to log
        #   We note that this may not strictly obey the formulation of the flux
        #   algorithm, but this tends to work okay in practice.
        log_arg = g[kp] / gk + 1e-60
        log_arg = max(1e-60, log_arg)
        x1 = log(log_arg)

        # We did not pre-scale the cached Courant limits for this version so we
        # have to apply them when estimating the flux directly.
        flux = gsk / x1 * (exp(0.5 * x1) - exp(x1 * (0.5 - (c[i, j]/3/Δy))))
        flux = min(flux, gk)
        g[k] = gk - flux
        g[kp] = g[kp] + flux
      end

    end # j
  end # i

end

println("""
BEGIN TIME INTEGRATION
""")
CPUtic()
for i ∈ 1:nt
  global t = t + Δt
  global tlmin = tlmin + Δt

  # Collision
  # subroutine coad
  coad!(gᵢ, xᵢ, c, ck, ima, Δt, Δy, debug=debug, gmin=gmin)

  # Plotting
  if tlmin ≥ 60
    global tlmin = tlmin - 60
    global lmin = lmin + 1

    if (lmin % Δt_plot) < 1 & do_plots
      display(plot!(p, rᵢ*1e6, gᵢ*1e3, label = "t = $lmin min"))
    end

    # Mass balance? Not sure what's going on here. Maybe numerical checking?
    x0 = 0.0
    x1 = 1.0
    imax = 0
    for i ∈ 1:m
      x0 = x0 + gᵢ[i] * Δy
      x1 = max(x1, gᵢ[i])
      if abs(x1 - gᵢ[i]) < 1e-9
        imax = i
      end
    end
    @printf "    %4d mins |" lmin
    @printf " mass %10.3e  max %10.3e imax %3d" x0 x1 imax
    @printf "\n"
  end

  # Mass balance? Not sure what's going on here. Maybe numerical checking?
  x0 = 0.0
  x1 = 1.0
  imax = 0
  for i ∈ 1:m
    x0 = x0 + gᵢ[i] * Δy
    x1 = max(x1, gᵢ[i])
    if abs(x1 - gᵢ[i]) < 1e-9
      imax = i
    end
  end
  @printf "    %4d s |" t
  @printf " mass %10.3e  max %10.3e  imax %3d" x0 x1 imax
  @printf "\n"

  local elapsed = CPUtoq()
  global total_runtime += elapsed
  CPUtic()
end
@printf "Total time - %5.2f seconds\n" total_runtime

println("End; press any key to close.")
xxx = readline()


text:
#
# Computes the norm of the projected gradient
# One method is defined for each possible scenario
# of bounds (all, only lower, only upper, none). This is
# needed such that the lower and upper bounds vectors can
# be of type Nothing if those bounds are not defined
#
function pr_gradnorm(g,x,lower,upper)
  if isnothing(lower) && isnothing(upper)
    return pr_gradnorm_no_bounds(g)
  elseif isnothing(lower)
    pr_gradnorm_upper(g,x,upper)
  elseif isnothing(upper)
    return pr_gradnorm_lower(g,x,lower)
  else
    return pr_gradnorm_both_bounds(g,x,lower,upper)
  end
end
# Both bounds
function pr_gradnorm_both_bounds(g,x,lower,upper)
  gnorm = zero(eltype(x))
  for i in eachindex(x)
    z = max(lower[i], min(upper[i], x[i]-g[i])) - x[i]
    gnorm = max(gnorm, abs(z))
  end
  return gnorm
end
# Only lower
function pr_gradnorm_lower(g,x,lower)
  gnorm = zero(eltype(x))
  for i in eachindex(x)
    z = max(lower[i], x[i]-g[i]) - x[i]
    gnorm = max(gnorm, abs(z))
  end
  return gnorm
end
# Only upper
function pr_gradnorm_upper(g,x,upper)
  gnorm = zero(eltype(x))
  for i in eachindex(x)
    z = min(upper[i], x[i]-g[i]) - x[i]
    gnorm = max(gnorm, abs(z))
  end
  return gnorm
end
# No bounds
function pr_gradnorm_no_bounds(g)
  gnorm = zero(eltype(g))
  for i in eachindex(g)
    gnorm = max(gnorm, abs(g[i]))
  end
  return gnorm
end


text:
@inline metric(κ::Kernel) = κ.metric

## Allows to iterate over kernels
Base.length(::Kernel) = 1
Base.iterate(k::Kernel) = (k,nothing)
Base.iterate(k::Kernel, ::Any) = nothing

# default fallback for evaluating a kernel with two arguments (such as vectors etc)
kappa(κ::Kernel, x, y) = kappa(κ, evaluate(metric(κ), transform(κ, x), transform(κ, y)))

### Syntactic sugar for creating matrices and using kernel functions
for k in [:ExponentialKernel,:SqExponentialKernel,:GammaExponentialKernel,:MaternKernel,:Matern32Kernel,:Matern52Kernel,:LinearKernel,:PolynomialKernel,:ExponentiatedKernel,:ZeroKernel,:WhiteKernel,:ConstantKernel,:RationalQuadraticKernel,:GammaRationalQuadraticKernel]
    @eval begin
        @inline (κ::$k)(d::Real) = kappa(κ,d) #TODO Add test
        @inline (κ::$k)(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = kappa(κ, x, y)
        @inline (κ::$k)(X::AbstractMatrix{T},Y::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,Y,obsdim=obsdim)
        @inline (κ::$k)(X::AbstractMatrix{T};obsdim::Integer=defaultobs) where {T} = kernelmatrix(κ,X,obsdim=obsdim)
    end
end

### Transform generics
@inline transform(κ::Kernel) = κ.transform
@inline transform(κ::Kernel, x) = transform(transform(κ), x)
@inline transform(κ::Kernel, x, obsdim::Int) = transform(transform(κ), x, obsdim)

## Constructors for kernels without parameters
for kernel in [:ExponentialKernel,:SqExponentialKernel,:Matern32Kernel,:Matern52Kernel,:ExponentiatedKernel]
    @eval begin
        $kernel() = $kernel(IdentityTransform())
        $kernel(ρ::Real) = $kernel(ScaleTransform(ρ))
        $kernel(ρ::AbstractVector{<:Real}) = $kernel(ARDTransform(ρ))
    end
end


text:
using LinearAlgebra

import Flux
import ZygoteRules

struct CachedDense{M <: AbstractVecOrMat, V <: AbstractVecOrMat, F}
   W :: M
   b :: V
   σ :: F
   cache :: NamedTuple{(:Wx, :y, :x̄, :W̄, :b̄), Tuple{M,M,M,M,V}}
end

Flux.trainable(m::CachedDense) = (m.W, m.b)

Flux.@functor CachedDense

CachedDense(W, b, cache) = CachedDense(W, b, identity, cache)

function CachedDense(in::Integer, out::Integer, σ = identity; batchsize::Integer,
               initW = Flux.glorot_uniform, initb = Flux.zeros)
   cache = (Wx = Flux.zeros(out, batchsize),
            y  = Flux.zeros(out, batchsize),
            x̄  = Flux.zeros(in, batchsize),
            W̄  = Flux.zeros(out, in),
            b̄  = Flux.zeros(out))
   return CachedDense(initW(out, in), initb(out), σ, cache)
end

function Base.show(io::IO, m::CachedDense)
   out, in = size(m.W)
   batchsize = size(m.cache.y, 2)
   print(io, "CachedDense($in, $out")
   m.σ == identity || print(io, ", ", m.σ)
   print(io, "; batchsize=$batchsize)")
end

(m::CachedDense)(x::AbstractVecOrMat) = cached_dense(x, m.W, m.b, m)
function cached_dense(x::AbstractVecOrMat, W::AbstractVecOrMat, b::AbstractVecOrMat, m::CachedDense)
   σ, c = m.σ, m.cache
   Wx = y = @view c.y[:,axes(x,2)]
   mul!(Wx, W, x)
   @. y = σ(Wx + b)
   return y
end

# Try to avoid hitting generic matmul in some simple cases
# Base's matmul is so slow that it's worth the extra conversion to hit BLAS
(m::CachedDense{M})(x::AbstractVecOrMat{T}) where {T <: Union{Float32,Float64}, M <: AbstractVecOrMat{T}} =
   invoke(m, Tuple{AbstractVecOrMat}, x)

(m::CachedDense{M})(x::AbstractVecOrMat{<:AbstractFloat}) where {T <: Union{Float32,Float64}, M <: AbstractVecOrMat{T}} =
   m(T.(x))

ZygoteRules.@adjoint function cached_dense(x::AbstractVecOrMat, W::AbstractVecOrMat, b::AbstractVecOrMat, m::CachedDense)
   σ, c = m.σ, m.cache
   W̄, b̄, batch = c.W̄, c.b̄, axes(x,2)
   @views Wx, y, x̄ = c.Wx[:,batch], c.y[:,batch], c.x̄[:,batch]
   W̄x̄⁺b̄ = Wx⁺b = Wx
   mul!(Wx, W, x)
   @. Wx⁺b = Wx + b
   @. y = σ(Wx⁺b)
   σ′ = σ'
   function cached_dense_adjoint(ȳ)
      @. W̄x̄⁺b̄ = ȳ * σ′(Wx⁺b)
      mul!(x̄, W', W̄x̄⁺b̄)
      mul!(W̄, W̄x̄⁺b̄, x')
      sum!(reshape(b̄, Val(2)), W̄x̄⁺b̄)
      return x̄, W̄, b̄, nothing
   end
   y, cached_dense_adjoint
end

max_stars_repo_name: AzamatB/CachedFluxLayers

max_issues_repo_name: AzamatB/CachedFluxLayers

max_forks_repo_name: AzamatB/CachedFluxLayers


text:
#=
# 8. Prescribe Traction
In this notebook we will show example of prescribing traction type boundary condition by using point source.
=#
using ViscousFlow
using Plots

#=
## Problem specification
Take $Re=100$ for example and discretize:
=#
Re = 100

#-
Δx,Δt = setstepsizes(Re,gridRe=4)
xlim, ylim = (-6,6),(-2,2)

#=
## Create channel
Two plates of length 30 centered at (0,0) are used to represent the wall of the pipe. All bodies must be put into a bodylist to create the Navier Stokes system.
=#
channel = BodyList()
tl = RigidTransformList()
push!(channel,Plate(10,1.5Δx)); push!(tl,RigidTransform((0.0,1.0),0.));
push!(channel,Plate(10,1.5Δx)); push!(tl,RigidTransform((0.0,-1.0),0.));
tl(channel)
plot(channel)

#=
## Inlet/Outlet boundary
The traction line source functions will set a specified traction strength onto a sets of specified location, so here, we use two veritcal plates with length of 2 (gap of between the wall) to represent the inlet and outlet. LineSourceParams( ) is used to collect the location of a body.
=#
inlet = Plate(2.0,1.5Δx)
T = RigidTransform((-5.0,0.0),π/2)
T(inlet)
τin = LineSourceParams(inlet)
TractionStrength_in = VectorData(collect(inlet))
fill!(TractionStrength_in.u,-1.0)
fill!(TractionStrength_in.v,0.0)
nothing #hide
#-

outlet = Plate(2.0,1.5Δx)
T = RigidTransform((5.0,0.0),π/2)
T(outlet)
τout = LineSourceParams(outlet)
TractionStrength_out = VectorData(collect(outlet))
fill!(TractionStrength_out.u,0.0)
fill!(TractionStrength_out.v,0.0)
nothing #hide
#-

plot!(inlet)
plot!(outlet)

#=
## Create N-S system
Note the order of $\tau$ in $\tau$_bc is important for prescribing their strength
=#
sys = NavierStokes(Re,Δx,xlim,ylim,Δt,channel,τ_bc=[τin,τout]);

#=
Initialize and Integrate
=#
u0 = newstate(sys)
tspan = (0.0,12.0)
integrator = init(u0,tspan,sys);

#=
## Set traction strength for inlet and outlet
Outlet is skipped because it is 0 by default.
=#
set_linesource_strength!(sys.τ_bc[1],TractionStrength_in);

#=
## Solve
=#
step!(integrator,12.0)

#=
## Examine
Plot the vorticity and streamlines
=#

plot(vorticity(integrator),sys,title="Vorticity",clim=(-10,10),levels=range(-10,10,length=50),color=:turbo,linewidth=1.5,ylim=ylim,fillrange=nothing)

#-

plot(streamfunction(integrator),sys,title="Streamlines",linewidth=1.5,levels=15,ylim=ylim,color=:Black,fillrange=nothing,linecolor=:black)

#=
Make a movie:
=#
sol = integrator.sol;
@gif for (u,t) in zip(sol.u,sol.t)
    plot(vorticity(u,sys,t),sys,clim=(-10,10),levels=range(-10,10,length=30),color=:turbo)
end every 5

max_stars_repo_path: test/literate/8.-Prescribe traction.jl

max_issues_repo_path: test/literate/8.-Prescribe traction.jl

max_forks_repo_path: test/literate/8.-Prescribe traction.jl