text:
using XInts
using XInts: SLimb, slimbmin, slimbmax, Limb, ClongMax, CulongMax, CdoubleMax, BITS_PER_LIMB,
             add!, sub!, mul!, com!, xor!, lshift!, rshift!, and!, ior!, SLimbW, LimbW, vec, is_short
using BitIntegers
using Random

function validate(x::XInt)
    @test x isa XInt
    if XInts.is_short(x)
        @test x.x != slimbmin
    else
        xv = x.v
        xl = abs(x.x)
        @test length(xv) >= xl >= 1
        @test !iszero(xv[xl])
        if xl == 1
            @test xv[1] >= XInts.limb1min
        end
    end
    x
end

validate(x) = @test x isa XInt

vint(x) = validate(XInt(x))

@testset "constructor" begin
    function test(a)
        x = validate(XInt(a))
        @test x isa XInt
        if a isa XInt
            @test x === a
        else
            @test x == a
        end
        if !(a isa BigInt)
            @test x == BigInt(a)
        end
    end

    @testset "XInt(::XInt)" begin
        test(XInt(2))
        test(XInt(big(2)^200))
    end
    @testset "XInt(::Bool)" begin
        test(true)
        test(false)
    end

    @testset "XInt(::$T)" for T in [Base.BitInteger_types...,
                                    Int256, UInt256, Int512, UInt512, Int1024, UInt1024]
        for x = rand(T, 20)
            test(x)
        end
        for d = T[0, 1, 2]
            test(d)
            if T <: Signed
                test(-d)
            end
            test(typemin(T) + d)
            test(typemax(T) - d)
            test((typemax(T) >>> 1) + d)
            test((typemax(T) >>> 1) - d)
        end
    end

    @testset "XInt(::BigInt)" begin
        for T = (Bool, Int8, Int, Int128)
            for x = BigInt.(rand(T, 5))
                test(x)
            end
        end
    end

    @testset "XInt(::$F)" for F in [Float16, Float32, Float64]
        for T = (Bool, Int8, Int32)
            for x = F.(rand(T, 5))
                isinteger(x) && test(x)
            end
        end
        @test XInt(F(2.0)) === XInt(2)
        @test_throws InexactError XInt(F(2.2))
    end
end

@testset "conversions" begin
    @testset "rem(::XInt, ::Type{Bool})" begin
        for x = Any[rand(Limb), big(rand(Int128)), rand(-big(2)^200:big(2)^200)]
            @test vint(x) % Bool === x % Bool
        end
    end
    @testset "rem(::XInt, ::Type{$T})" for T in Base.BitInteger_types
        xs = T[0, 2, rand(T, 10)..., typemax(T)]
        if T <: Signed
            push!(xs, -1, -2, typemin(T))
        end
        for x = xs
            @test vint(x) % T === x
            @test vint(big(rand(T)) << (8*sizeof(T)) +
                       big(x)) % T === x
            @test vint(big(rand(T)) << (16*sizeof(T)) +
                       big(rand(T)) << (8*sizeof(T)) +
                       big(x)) % T === x
        end
    end

    @testset "Bool(::XInt)" begin
        @test Bool(XInt(0)) == false
        @test Bool(XInt(1)) == true
        for x = Any[rand(Int8, 10); rand(Int128, 10)]
            (x == 0 || x == 1) && continue
            @test_throws InexactError Bool(XInt(x))
        end
    end
    @testset "$T(::XInt)" for T in Base.BitInteger_types
        for x = rand(T, 10)
            @test T(vint(x)) === x
        end
        if sizeof(T) > sizeof(SLimb)
            for x = T.(rand(T <: Signed ? Int32 : UInt32, 10))
                @test T(vint(x)) === x
            end
        end
        @test T(vint(typemax(T))) === typemax(T)
        @test T(vint(typemin(T))) === typemin(T)
        @test_throws InexactError T(XInt(typemax(T))+XInt(1))
        @test_throws InexactError T(XInt(typemax(T))+XInt(+10))
        @test_throws InexactError T(XInt(typemin(T))+XInt(-1))
        @test_throws InexactError T(XInt(typemin(T))+XInt(-10))
        @test_throws InexactError T(XInt(2)^200) # for sizeof(T) < sizeof(Limb)
        # test un-reduced XInt
        @test T(XInts._XInt(1, Limb[1])) === T(1)
    end

    @testset "BigInt(::XInt)" begin
        for T = (Bool, Int8, Int, Int128)
            for x = rand(T, 5)
                z = BigInt(x)
                y = vint(x)
                @test z == BigInt(y) isa BigInt
                @test z == big(y) isa BigInt
                @test z == y % BigInt isa BigInt
            end
        end
    end

    @testset "$F(::XInt, $R)" for
        F in (Float16, Float32, Float64, BigFloat, AbstractFloat, float),
        R in ((F===float) ? ("",) : (RoundToZero, RoundDown, RoundUp, RoundNearest, ""))

        for T = Base.BitInteger_types
            for x = rand(T, 10)
                if R == ""
                    @test F(vint(x)) == F(BigInt(x))
                    # have to use BigInt because of a Julia bug
                else
                    @test F(vint(x), R) == F(BigInt(x), R)
                end
            end
        end
    end
end

@testset "math" begin
    @testset "isodd/iseven" begin
        for x = Any[rand(Int8, 5)...,
                    rand((-1, 1)) .* rand(0:big(2)^200, 5)...]
            y = vint(x)
            @test iseven(y) == iseven(x) != isodd(y)
        end
    end

    @testset "iszero/isone" begin
        @test iszero(XInt(0))
        @test zero(XInt) === XInt(0)
        @test isone(XInt(1))
        @test one(XInt) === XInt(1)
        x = rand(Int8)
        @test iszero(x) == iszero(XInt(x))
        @test isone(x) == isone(XInt(x))
    end

    @testset "sign" begin
        @test sign(zero(XInt)) === XInt(0)
        @test signbit(zero(XInt)) === false
        @test signbit(one(XInt)) === false
        @test signbit(-one(XInt)) === true

        for x in rand(Int8, 10)
            y = XInt(x)
            @test sign(y) isa XInt
            @test sign(y) == sign(x)
            @test signbit(y) == signbit(x)
        end
    end

    @testset "digits" begin
        for r = big(2).^[5, 8, 64, 200]
            for b = rand(-r:r, 5)
                pad = rand(0:2)
                for base = [rand(8:12, 3); 100] # 100 > 62, cf. implementation of digits!
                    @test ndigits(b, pad=pad, base=base) == ndigits(XInt(b), pad=pad, base=base)
                    @test digits(b, pad=pad, base=base) == digits(XInt(b), pad=pad, base=base)
                end
            end
        end
    end

    @testset "prevpow & nextpow" begin
        # TODO: deduplicate this xs definition (from binary ops)
        xs = Any[1, 2, rand(3:1000, 10)..., slimbmax, slimbmax-1, slimbmax-2,
                 rand(Int128(2)^65:Int128(2)^100, 2)..., rand(big(3):big(2)^200, 2)...,
                 1024, big(2)^100, big(2)^200] # popcount == 1

        for x = xs, b = Int[2, 3, 4, 5, 9, 10, rand(11:100, 6)...]
            y = validate(nextpow(XInt(b), XInt(x)))
            @test y isa XInt
            @test y == nextpow(big(b), big(x))
            if b == 2
                # accepts other integer types for base
                @test y == nextpow(2, XInt(x))
            end
            y = validate(prevpow(XInt(b), XInt(x)))
            @test y isa XInt
            @test y == prevpow(big(b), big(x))
            if b == 2
                # accepts other integer types for base
                @test y == validate(prevpow(2, XInt(x)))
            end
        end
        for b = XInt[1, 0, -1, rand(-100:-2, 5)...]
            @test_throws DomainError prevpow(b, XInt(rand(1:100)))
            @test_throws DomainError nextpow(b, XInt(rand(1:100)))
        end
        for x = [0, -1, rand(-100:-2, 5)...]
            @test_throws DomainError prevpow(XInt(rand(2:20)), vint(x))
            @test_throws DomainError nextpow(XInt(rand(2:20)), vint(x))
        end

        # add some x which lead to increased size (in limbs)
        append!(xs, big.(rand(UInt, 4)) .<< (BITS_PER_LIMB .* [0, 1, 2, 3])')
        xs = [xs; .-(xs)]
        push!(xs, slimbmin)
        for x=xs
            y = validate(Base._nextpow2(XInt(x)))
            @test y isa XInt
            @test y == Base._nextpow2(big(x))
            y = validate(Base._prevpow2(XInt(x)))
            @test y isa XInt
            @test y == Base._prevpow2(big(x))
        end
    end

    @testset "sum" begin
        function testsum(bs)
            bs0 = deepcopy(bs)
            xs = XInt.(bs)
            xs0 = deepcopy(xs)
            sumbs = sum(bs)
            sumxs = sum(xs)
            # we test that sum didn't mess up with internal mutation
            @test is_short(sumxs) || all(x -> x !== sumxs, xs)
            @test xs == xs0 == bs0
            @test sumbs == sumxs
        end
        xs = BigInt.(rand(Int8, rand(1:10)))
        testsum(xs)
        append!(xs, rand(SLimb, rand(1:10)))
        testsum(shuffle!(xs))
        append!(xs, rand(SLimbW, rand(1:10)))
        testsum(shuffle!(xs))
        append!(xs, rand(-big(2)^200:big(2)^200, rand(1:10)))
        testsum(shuffle!(xs))
        # the first number x is a non-short, so when added with the initial
        # 0 in sum, it will return x itself: we must not mutate it later on
        xs = rand(-big(2)^65:big(2)^65, 10)
        testsum(xs)
        # test only big positive or negative numbers
        axs = abs.(xs)
        testsum(axs)
        testsum(.-axs)
        # with some shorts
        append!(axs, rand(slimbmin+1:slimbmax, 10))
        testsum(axs)
        testsum(.-axs)

        # tuples
        z1, z2, z3, z4 = zs = XInt.(xs[1:4])
        @test sum(zs) == sum(Tuple(zs))
        @test sum(zs) == z1 + z2 + z3 + z4
    end
end

@testset "misc" begin
    @testset "[deep]copy, signed, Signed, widen, hastypemax" begin
        # automatic implementation from XInt <: Signed
        for x = XInt[2, big(2)^70]
            @test copy(x) === x
            let y = validate(deepcopy(x))
                @test y == x
                @test x == 2 ? x === y : y !== x
            end
            @test signed(x) === x
            @test Signed(x) === x
            @test widen(x) === x
        end
        @test widen(XInt) == XInt
        @test !Base.hastypemax(XInt)
        if VERSION >= v"1.5"
            @test signed(XInt) == XInt
        end
    end
    @testset "parse, string, show" begin
        @test validate(parse(XInt, "123")) === XInt(123)

        for x = BigInt[-10:10;
                       rand(Int8, 10);
                       rand(Int128, 10);
                       rand(-big(2)^200:big(2)^200, 10)]
            @test string(x) == string(XInt(x))
            @test sprint(show, x) == sprint(show, XInt(x))
            for base = [-62:-2; 2:62], pad = rand(0:220, 30)
                @test string(x, base=base, pad=pad) == string(XInt(x), base=base, pad=pad)
            end
        end
    end
    @testset "float(XInt)" begin
        @test float(XInt) === BigFloat
    end
end

@testset "trunc" begin
    for T in [#=Float16,=# Float32, Float64] # TODO: add Float16, and also for BigInt
        @test validate(trunc(XInt, T(2.2))) === XInt(2)
        @test validate(trunc(XInt, T(-2.2))) === XInt(-2)
        @test_throws InexactError trunc(XInt, T(NaN))
        @test_throws InexactError trunc(XInt, T(Inf))
        @test validate(unsafe_trunc(XInt, T(2.2))) === XInt(2)
        @test validate(unsafe_trunc(XInt, T(-2.2))) === XInt(-2)
    end
end

@testset "comparisons" begin
    @testset "== (integers)" begin
        xs = BigInt[slimbmin, 0, 1, 2, slimbmax-1, slimbmax, slimbmax+1, big(2)^100]
        xs = [xs..., (-).(xs)...]
        for (a, b) in Iterators.product(xs, xs)
            x, y = XInt(a), XInt(b)
            @test x == x
            @test (a == b) == (x == y)
            @test a == x && x == a
            @test b == y && y == b
            for T in Base.BitInteger_types
                if typemin(T) <= b <= typemax(T)
                    z = T(b)
                    @test (x == z) == (a == z) == (z == x)
                end
            end
        end
    end
    @testset "XInt with floats" begin
        for T = Base.uniontypes(CdoubleMax)
            b = Int(maxintfloat(T))
            for x = rand(0:b, 10)
                @test vint(x) == T(x)
                @test T(x) == vint(x)
                @test vint(x-1) < T(x)
                @test T(x) < vint(x+1)
                @test T(x) <= vint(x) <= T(x)

            end
        end
    end
end

opmap(x) = x
opmap(::typeof(add!)) = +
opmap(::typeof(sub!)) = -
opmap(::typeof(mul!)) = *
opmap(::typeof(com!)) = ~
opmap(::typeof(lshift!)) = <<
opmap(::typeof(rshift!)) = >>
opmap(::typeof(and!)) = &
opmap(::typeof(ior!)) = |
opmap(::typeof(xor!)) = xor

function testop2(op, x, y)
    for a = (-x, x), b = (-y, y)
        s = if op === invmod
            try
                op(big(a), big(b))
            catch
                @test_throws DomainError op(XInt(a), XInt(b))
                continue
            end
        else
            opmap(op)(big(a), big(b))
        end
        @assert !(a isa XInt) # so that if op is mutating, it doesn't mutate a
        r = op(vint(a), vint(b))
        if op === divrem
            @test r isa Tuple{XInt,XInt}
            validate.(r)
        elseif op === gcdx
            @test r isa Tuple{XInt,XInt,XInt}
            validate.(r)
        elseif !(s isa BigInt)
            @test typeof(s) == typeof(r)
        else
            validate(r)
        end
        @test s == r
        if op === div
            for R = (RoundToZero, RoundDown, RoundUp)
                @test div(big(a), big(b), R) ==
                    validate(div(XInt(a), XInt(b), R)) isa XInt
            end
        end
        if opmap(op) != op
            # mutating
            tmp = XInts._XInt(0, Limb[])
            aa, bb = vint(a), vint(b)
            rtmp = validate(op(tmp, aa, bb))
            @test rtmp == r
            if !is_short(r)
                @test vec(rtmp) === vec(tmp) || !is_short(aa) && vec(rtmp) === vec(aa) ||
                    vec(rtmp) === vec(bb)
                aa2 = XInts._copy(aa)
                # test aliasing works between the source(s) and destination
                if aa != bb
                    rtmp = validate(op(aa2, aa2, bb))
                else
                    rtmp = validate(op(aa2, aa2, aa2))
                end
                @test rtmp == r
            end
            tmp = XInts._XInt(0)
            @test validate(op(tmp, vint(a), vint(b))) == r
        end
    end
end

function test_alloc(op, x, y, r=nothing)
    if VERSION >= v"1.5"
        @test 0 == @allocated op(x, y)
    end
    z = op(x, y)
    @test validate(z) == r
end

@testset "operations" begin

    xs = BigInt[0, 1, 2, rand(0:1000, 5)..., slimbmax, slimbmax-1, slimbmax-2,
                rand(SLimb, 5)..., rand(Limb, 5)...,
                Int128(slimbmax)+1, Int128(slimbmax)+2,
                (LimbW.(rand(Limb, 3)) .<< BITS_PER_LIMB)...,
                typemax(LimbW), LimbW(typemax(Limb)) << BITS_PER_LIMB, # and!
                (LimbW(typemax(Limb)) << BITS_PER_LIMB) + 1, # and!
                (let h = rand(Limb) # ior!
                     [big(h) << 2BITS_PER_LIMB | big(rand(LimbW)),
                      big(~h) << 2BITS_PER_LIMB | big(rand(LimbW)),
                      ]
                 end)...,
                (Int128(slimbmax).+rand(UInt8, 5))...,
                rand(Int128(2)^65:Int128(2)^100, 2)..., rand(big(0):big(2)^200, 2)...,
                slimbmin] # TODO: slimbmin and 0 are used twice when negated in test(...)

    # add!
    let
        x4 = rand(1:typemax(Limb))
        y4 = x4-1
        z = Limb(1)<<63
        y1 = rand(z : typemax(Limb))
        x1 = rand(Limb(0) : y1-z) # y1-z >= 0
        # y1-x1 >= z, so x1-y1 doesn't have top bit set
        x2 = rand(1:typemax(Limb))
        y2 = x2-1 # 2nd limb must absorb carry, and become 0
        x = XInt(x4) << 3BITS_PER_LIMB | XInt(x2) << BITS_PER_LIMB | x1
        y = XInt(y4) << 3BITS_PER_LIMB | XInt(y2) << BITS_PER_LIMB | y1
        u2 = rand(0:x2-1)
        u = XInt(y4) << 3BITS_PER_LIMB | XInt(u2) << BITS_PER_LIMB | y1
        # w, v: similar x and y, but with a 3rd limb, to have carry propagate
        w3 = rand(1:typemax(Limb))
        w = x | XInt(w3) << 2BITS_PER_LIMB
        v2 = rand(x2:typemax(Limb))
        v3 = rand(0:w3-1)
        v = XInt(y4) << 3BITS_PER_LIMB | XInt(v3) << 2BITS_PER_LIMB | XInt(v2) << BITS_PER_LIMB | y1
        push!(xs, x, y, u, v, w)
    end
    let # xor!
        # x < 0 with 3 limbs, y > 0 with 2 limbs, result with null 3rd limb
        # formula: x ⊻ y = -((-x-1) ⊻ y + 1), so find such result via
        # propagating the carry of +1 till 3rd limb
        x1 = rand(1:typemax(Limb))
        y1 = ~(x1-1)
        y2 = rand(1:typemax(Limb)) # y2 must not be null
        x2 = ~y2
        x3 = typemax(Limb)
        x = XInt(x3) << 2BITS_PER_LIMB | XInt(x2) << BITS_PER_LIMB | XInt(x1)
        y = XInt(y2) << BITS_PER_LIMB | XInt(y1)
        push!(xs, x, y)
    end

    @testset "$op(::XInt, ::XInt)" for op = (+, -, *, mod, rem, gcd, gcdx, lcm, &, |, xor,
                                             /, div, divrem, fld, cld, invmod,
                                             cmp, <, <=, >, >=, ==, flipsign, binomial,
                                             add!, sub!, mul!, and!, ior!, xor!)
        for x=xs, y=xs
            iszero(y) && op ∈ (/, mod, rem, div, divrem, fld, cld, invmod) &&
                continue
            op === binomial && !all(0 .<= [x, y] .<= 1000) && continue
            testop2(op, x, y)
        end
        if op === +
            x = XInt(2)^128
            y = validate(-x+1)
            @test validate(x+y) == big(x) + big(y)
        end

        # special cases
        if op === &
            # special case where the result takes one more limb than the operands
            x = -XInt(XInt(typemax(UInt))<<64 + 1)
            y = -XInt(2)^64
            z = validate(x & y)
            @test z == big(x) & big(y)
        elseif op === |
            x = -(XInt(2)^128)
            y = XInt(2)^128-1
            test_alloc(|, x, y, -1)
        end
    end

    function make_values(T, n=5)
        S = T === BigInt ? Int128 : T
        as = T[0, 1, 2, 3, rand(S, 2n)..., typemax(S), typemax(S)-1, typemax(S)-2]
        xs = BigInt.(filter(isinteger, as))
        if T == Limb # xor!
            x1 = rand(1:typemax(T))
            y = ~(x1-1)
            x = BigInt(typemax(T)) << BITS_PER_LIMB | x1
            push!(as, y)
            push!(xs, x)
        end
        push!(xs, BigInt(2)^BITS_PER_LIMB) # also for xor!
        if T <: Signed
            append!(as, (-).(as))
            push!(as, typemin(S))
        end
        if T <: Integer
            push!(xs, typemin(S))
            append!(xs, rand(S, 2n))
        end
        append!(xs, rand(big(2)^65:big(2)^100, n))
        append!(xs, rand(big(0):big(2)^200, n))
        as, xs
    end

    @testset "$op(::XInt, ::$T) / $op(::$T, ::XInt)" for
        op in (+, -, *, /, cmp, <, <=, >, >=, ==, !=, flipsign, gcd, gcdx, binomial,
               &, |, xor),
        T in [Base.BitInteger_types...; BigInt; Int256; UInt256;
              if op ∈ (gcd, gcdx, binomial, &, |, xor)
                  []
              else
                  [Base.uniontypes(CdoubleMax); BigFloat]
              end]
        as, xs = make_values(T)
        for a = as, y = xs, z = (-y, y), x = (vint(z),)
            op === binomial && !all(-1000 .<= [a, z] .<= 1000) && continue
            for (r, s) = ((op == binomial) ?
                          [(op(x, a), op(z, a))] : # TODO: remove special case
                          Any[(op(a, x), op(a, z)),
                              (op(x, a), op(z, a))])

                if op === flipsign
                    if r isa XInt
                        @test s isa BigInt
                    else
                        @test typeof(s) == typeof(r)
                    end
                elseif s isa BigInt
                    @test validate(r) isa XInt
                elseif op === gcdx
                    @test r isa Tuple{XInt,XInt,XInt}
                    validate.(r)
                else
                    @test typeof(r) == typeof(s)
                end
                @test isequal(r, s) || iszero(r) && iszero(s) && T == BigFloat
                # not `==` for NaN
                # TODO: fix the 0.0 vs -0.0 discrepancy with BigFloat?
            end
        end
    end

    @testset "powermod(::$T/XInt, ::$T/XInt, ::XInt)" for
            T in [Base.uniontypes(CulongMax); Base.uniontypes(ClongMax); [BigInt]]
        # TODO: add more combinations for first two parameters
        as1, xs1 = make_values(T, 2)
        as2, xs2 = make_values(T, 2)
        for a=as1, b=filter(>=(0), as2), x=filter(!iszero, xs1)
            s = powermod(a, b, x)
            r = powermod(a, b, vint(x))
            @test s == r
            @test validate(r) isa XInt
        end
        for a=as1, x=filter(>=(0), xs1), y=filter(!iszero, xs2)
            s = powermod(a, x, y)
            r = powermod(a, vint(x), vint(y))
            @test s == r
            @test validate(r) isa XInt
            r = powermod(vint(a), XInt(x), XInt(y))
            @test s == r
            @test validate(r) isa XInt
        end
    end

    @testset "checked operations" begin
        xs = XInt.([rand(1:1000, 4); rand(1:big(2)^65, 5)])
        xs = [xs; .-(xs)]
        fns = [getfield(Base, f[1]) => f[2] for f in
               [:checked_add => +, :checked_sub => -, :checked_mul => *,
                :checked_div => div, :checked_rem => rem, :checked_fld => fld,
                :checked_mod => mod, :checked_cld => cld]]

        for x=xs
            @test Base.checked_abs(x) == abs(x)
            @test Base.checked_neg(x) == -x
            for y=xs, (c, f) = fns
                @test f(x, y) == c(x, y)
            end
            for y=xs
                @test Base.add_with_overflow(x, y) == (x+y, false)
                @test Base.sub_with_overflow(x, y) == (x-y, false)
                @test Base.mul_with_overflow(x, y) == (x*y, false)
            end
        end
    end
end

@testset "bit and unary ops, etc." begin
    xs = BigInt[1, 2, 3, rand(UInt8, 5)..., rand(UInt, 5)..., rand(SLimb, 10)...,
                typemax(UInt), typemax(UInt)-1, typemax(UInt)-2,
                typemax(SLimb), typemax(SLimb)-1, typemax(SLimb)-2, typemax(SLimb)-3,
                big(2)^192-1] # specifically for rshift!

    append!(xs, rand(big(2)^65:big(2)^100, 5))
    append!(xs, rand(big(0):big(2)^200, 5))
    append!(xs, (-).(xs))
    push!(xs, 0, typemin(SLimb), typemin(SLimb) >> 1, big(typemin(SLimb))-1)

    @testset "$op(::XInt)" for op = (-, ~, isqrt, trailing_zeros, trailing_ones, count_ones,
                                     abs, factorial, com!)

        @assert xs isa Vector{BigInt} # so that mutating ops don't mutate xs themselves
        for x = xs
            op === factorial && !(0 <= x <= 1000) && continue
            s = try
                opmap(op)(x)
            catch
                @test_throws Union{InexactError,DomainError,DivideError} op(XInt(x))
                continue
            end
            r = op(vint(x))
            @test r == s
            if s isa BigInt
                @test validate(r) isa XInt
            else
                @assert s isa Int
                @test r isa Int
            end
        end
    end

    cs = [0, 1, 2, 3, 4, BITS_PER_LIMB, 2*BITS_PER_LIMB, rand(5:4*BITS_PER_LIMB, 20)...]
    cs2 = [cs; (-).(cs)]
    @testset "$op(::XInt, c)" for op = (<<, >>, >>>, (^), lshift!, rshift!)
        for T = [Base.BitInteger_types..., BigInt, Bool, XInt]
            for c = (T <: Signed && op !== (^) ? cs2 :
                     T === Bool ?                [true, false] :
                                                 cs)
                !Base.hastypemax(T) || typemin(T) <= c <= typemax(T) || continue
                for x = xs
                    s = opmap(op)(x, T(c))
                    @assert !(x isa XInt) # for mutation
                    r = op(vint(x), T(c))
                    @test s == r
                    validate(r)
                    if opmap(op) != op
                        # mutating
                        tmp = XInts._XInt(0, Limb[])
                        xx = vint(x)
                        rtmp = validate(op(tmp, xx, T(c)))
                        @test rtmp == r
                        if !is_short(r)
                            @test vec(rtmp) === vec(tmp) || vec(rtmp) === vec(xx)
                        end
                        tmp = XInts._XInt(0)
                        @test validate(op(tmp, vint(x), T(c))) == r
                    end
                end
            end
            if op === (^) && !Base.hastypemax(T) || sizeof(T) > sizeof(Culong)
                c = T(typemax(Culong)) + rand(1:1000)
                z = XInt(0)
                @test z^c === z
                @test_throws OverflowError XInt(2)^c
            end
        end
    end
end

@testset "hashing" begin
    xs = BigInt[rand(UInt, 10); rand(UInt128, 10); rand(1:big(2)^300, 5)]
    xs = [xs; .-(xs)]
    for x=xs
        h = rand(UInt)
        a = XInt(x)
        b = big(x)

        @test hash(b) == hash(a)
        @test hash(b, h) == hash(a, h)
        @test Base.hash_integer(b, h) == Base.hash_integer(a, h)

        for y=xs
            iszero(y) && continue
            a2 = XInt(y)
            b2 = big(y)
            @test hash(b//b2) == hash(a//a2)
            @test hash(b//b2, h) == hash(a//a2, h)
        end
    end
end

@testset "rand" begin
    xs = BigInt[[0, 1, 2]; rand(Int8, 5); rand(SLimb, 10); rand(Int128, 5);
                rand(big(1):big(2)^300, 5)]
    as = [xs; .-(xs)]
    rng = Random.MersenneTwister()

    for a=as, x=xs, b = (x, a+x)
        # tests (a big, b small) and vice-versa, and length(r) small even if a, b bigs
        r = XInt(a):XInt(b)
        if isempty(r)
            @test_throws ArgumentError rand(r)
        else
            @test validate(rand(r)) ∈ r
            sp = Random.Sampler(Random.MersenneTwister, r)
            z = XInt(0)
            @test validate(Random.rand!(rng, z, sp)) ∈ r
        end
    end
end


text:
using MATLAB
using Base.Test

include("engine.jl")
include("matfile.jl")
include("matstr.jl")
include("mstatements.jl")
include("mxarray.jl")

ext: jl


text:
using Test

using Interact

@test gettheme() == Interact.Bulma()
settheme!(:nativehtml)
@test gettheme() == InteractBase.NativeHTML()
settheme!(:bulma)
@test gettheme() == Interact.Bulma()

@test isfile(Interact.bulma_css())
@test isfile(Interact.bulma_confined_css())

@test first(eachline(Interact.bulma_css()))[1:20] == "/*! bulma.io v0.7.4 "
@test first(eachline(Interact.bulma_confined_css()))[1:20] == ".interact-widget{/*!"

@test length(InteractBase.libraries(Interact.Bulma())) == 3
@test all(isfile, InteractBase.libraries(Interact.Bulma()))


text:
export Newton

"""
A globalized Newton algorithm with Line Search
Inputs:
    - An AbstractNLPModel (our problem)
    - An NLPStopping (stopping fonctions/criterion of our problem)
    - (opt) hessian_rep, different way to implement the Hessian. The default
        way is a LinearOperator
Outputs:
    - An NLPAtX, so all the information of the last iterate
    - A boolean true if we reahced an optimal solution, false otherwise
"""
function Newton(nlp             :: AbstractNLPModel,
                 nlp_stop       :: NLPStopping;
                 linesearch     :: Function = armijo_ls,
                 Nwtdirection   :: Function = NwtdirectionLDLT,
                 hessian_rep    :: Function = hessian_dense,
                 kwargs...)

    nlp_at_x = nlp_stop.current_state

    n = nlp.meta.nvar

    xt = copy(nlp.meta.x0)
	T = eltype(xt)
    ∇ft = Array{Float64}(undef, n)

    f = obj(nlp, xt)
    ∇f = grad(nlp, xt)
	nlp_stop.meta.optimality0 = norm(copy(∇f))
    OK = update_and_start!(nlp_stop, x = xt, fx = f, gx = ∇f)#, g0 = ∇f)
    # ∇fNorm = BLAS.nrm2(n, nlp_at_x.gx, 1)
    ∇fNorm = norm(nlp_at_x.gx)
    !OK && Stopping.update!(nlp_at_x, Hx = hessian_rep(nlp, xt))

    iter = 0

    @info log_header([:iter, :f, :nrm_g, :α], [Int64, T, T, T])
	@info log_row(Any[iter, nlp_at_x.fx, ∇fNorm, 0.0])
    β = 0.0
    # d = zeros(nlp_at_x.gx)
    d = zero(nlp_at_x.gx)

    h = LineSearch.LineModel(nlp, xt, d)

    while !OK
        d = Nwtdirection(nlp_at_x.Hx, nlp_at_x.gx)
        # slope = BLAS.dot(n, d, 1, nlp_at_x.gx, 1)
        slope = d' * nlp_at_x.gx

        h = LineSearch.redirect!(h, xt, d)

        ls_at_t = LSAtT(0.0, h₀ = nlp_at_x.fx, g₀ = slope)
        stop_ls = LS_Stopping(h, (x, y) -> armijo(x, y, τ₀ = 1e-09), ls_at_t)
        ls_at_t, good_step_size = linesearch(h, stop_ls; kwargs...)
        # good_step_size || (nlp_stop.meta.stalled_linesearch = true)

        xt = nlp_at_x.x + ls_at_t.x * d
        ft = obj(nlp, xt); ∇ft = grad(nlp, xt)

        # BLAS.blascopy!(n, nlp_at_x.x, 1, xt, 1)
        # BLAS.axpy!(n, ls_at_t.x, d, 1, xt, 1) #BLAS.axpy!(n, t, d, 1, xt, 1)
        ∇ft = grad!(nlp, xt, ∇ft)

        # Move on.
        OK = update_and_stop!(nlp_stop, x = xt, fx = ft, gx = ∇ft, Hx = hessian_rep(nlp, xt))

        # norm(∇f) bug: https://github.com/JuliaLang/julia/issues/11788
        # ∇fNorm = BLAS.nrm2(n, nlp_at_x.gx, 1)
        ∇fNorm = norm(nlp_at_x.gx)
        iter = iter + 1
		@show nlp_at_x.x[1]
		@show nlp_at_x.x[2]
		@info log_row(Any[iter, nlp_at_x.fx, ∇fNorm, ls_at_t.x])
    end

    return nlp_stop, nlp_stop.meta.optimal
end


text:
"""
    radiality_centrality(g)

Calculate the [radiality centrality](http://www.cbmc.it/fastcent/doc/Radiality.htm)
of a graph `g` across all vertices. Return a vector representing the centrality
calculated for each node in `g`.

The radiality centrality ``R_u`` of a vertex ``u`` is defined as
``
R_u = \\frac{D_g + 1 - \\frac{\\sum_{v∈V}d_{u,v}}{|V|-1}}{D_g}
``

where ``D_g`` is the diameter of the graph and ``d_{u,v}`` is the
length of the shortest path from ``u`` to ``v``.

### References
- Brandes, U.: A faster algorithm for betweenness centrality. J Math Sociol 25 (2001) 163-177

# Examples
```jldoctest
julia> using LightGraphs

julia> radiality_centrality(StarGraph(4))
4-element Array{Float64,1}:
 1.0
 0.6666666666666666
 0.6666666666666666
 0.6666666666666666

julia> radiality_centrality(PathGraph(3))
3-element Array{Float64,1}:
 0.75
 1.0
 0.75
```
"""
function radiality_centrality(g::AbstractGraph)::Vector{Float64}
    n_v = nv(g)
    vs = vertices(g)

    meandists = zeros(Float64, n_v)
    dmtr = 0.0
    for v in vs
        d = dijkstra_shortest_paths(g, v)
        dmtr = max(dmtr, maximum(d.dists))
        meandists[v] = sum(d.dists) / (n_v - 1) # ignore the source vx
    end
    meandists = (dmtr + 1) .- (meandists)
    return meandists ./ dmtr
end


text:
# This script gets run once, on a developer's machine, to generate artifacts
using Pkg.Artifacts

artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")

rm(artifacts_toml)

plotschema_url = "https://api.plot.ly/v2/plot-schema?sha1"
plotschema_hash = create_artifact() do artifact_dir
    download(plotschema_url, joinpath(artifact_dir, "plotschema.json"))
end
bind_artifact!(artifacts_toml, "plotschema", plotschema_hash; download_info = [
    (plotschema_url, plotschema_hash)
])

plotly_url = "https://cdn.plot.ly/plotly-latest.min.js"
plotly_hash = create_artifact() do artifact_dir
    download(plotly_url, joinpath(artifact_dir, "plotly-latest.min.js"))
end
bind_artifact!(artifacts_toml, "plotly", plotly_hash; download_info = [
    (plotly_url, plotly_hash)
])


text:
using ProgressMeter, KitBase
import KitML

begin
    # space
    x0 = -1.5
    x1 = 1.5
    y0 = -1.5
    y1 = 1.5
    nx = 50#100
    ny = 50#100
    dx = (x1 - x0) / nx
    dy = (y1 - y0) / ny

    pspace = KitBase.PSpace2D(x0, x1, nx, y0, y1, ny)

    # time
    tEnd = 1.0
    cfl = 0.5
    dt = cfl / 2 * (dx * dy) / (dx + dy)

    # quadrature
    quadratureorder = 5
    points, triangulation = KitBase.octa_quadrature(quadratureorder)
    weights = KitBase.quadrature_weights(points, triangulation)
    nq = size(points, 1)

    # particle
    SigmaS = 1 * ones(ny + 4, nx + 4)
    SigmaA = 0 * ones(ny + 4, nx + 4)
    SigmaT = SigmaS + SigmaA

    # moments
    L = 1
    ne = (L + 1)^2
    phi = zeros(ne, nx, ny)
    α = zeros(ne, nx, ny)
    m = KitBase.eval_spherharmonic(points, L)
end

begin
    s2 = 0.03^2
    flr = 1e-4
    init_field(x, y) = max(flr, 1.0 / (4.0 * pi * s2) * exp(-(x^2 + y^2) / 4.0 / s2))
    for j = 1:nx
        for i = 1:ny
            y = y0 + (i - 0.5) * dy
            x = x0 + (j - 0.5) * dx
            # only zeroth order moment is non-zero
            phi[1, i, j] = init_field(x, y)
        end
    end
end

global t = 0.0
flux1 = zeros(ne, nx + 1, ny)
flux2 = zeros(ne, nx, ny + 1)

@showprogress for iter = 1:20
    # regularization
    Threads.@threads for j = 1:ny
        for i = 1:nx
            res = KitBase.optimize_closure(
                α[:, i, j],
                m,
                weights,
                phi[:, i, j],
                KitBase.maxwell_boltzmann_dual,
            )
            α[:, i, j] .= res.minimizer

            phi[:, i, j] .= KitBase.realizable_reconstruct(
                res.minimizer,
                m,
                weights,
                KitBase.maxwell_boltzmann_dual_prime,
            )
        end
    end

    # flux
    fη1 = zeros(nq)
    for j = 1:ny
        for i = 2:nx
            KitBase.flux_kfvs!(
                fη1,
                KitBase.maxwell_boltzmann_dual.(α[:, i-1, j]' * m)[:],
                KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:],
                points[:, 1],
                dt,
            )

            for k in axes(flux1, 1)
                flux1[k, i, j] = sum(m[k, :] .* weights .* fη1)
            end
        end
    end

    fη2 = zeros(nq)
    for i = 1:nx
        for j = 2:ny
            KitBase.flux_kfvs!(
                fη2,
                KitBase.maxwell_boltzmann_dual.(α[:, i, j-1]' * m)[:],
                KitBase.maxwell_boltzmann_dual.(α[:, i, j]' * m)[:],
                points[:, 2],
                dt,
            )

            for k in axes(flux2, 1)
                flux2[k, i, j] = sum(m[k, :] .* (weights .* fη2))
            end
        end
    end

    # update
    for j = 2:ny-1
        for i = 2:nx-1
            for q = 1:ne
                phi[q, i, j] =
                    phi[q, i, j] +
                    (flux1[q, i, j] - flux1[q, i+1, j]) / dx +
                    (flux2[q, i, j] - flux2[q, i, j+1]) / dy #+
                #(integral - phi[q, i, j]) * dt
            end
        end
    end

    global t += dt
end

X = zeros(Float32, 4, nx * ny)
Y = zeros(Float32, 4, nx * ny)
for j = 1:ny
    for i = 1:nx
        res = KitBase.optimize_closure(
            α[:, i, j],
            m,
            weights,
            phi[:, i, j],
            KitBase.maxwell_boltzmann_dual,
        )
        α[:, i, j] .= res.minimizer
        phi[:, i, j] .= KitBase.realizable_reconstruct(
            res.minimizer,
            m,
            weights,
            KitBase.maxwell_boltzmann_dual_prime,
        )

        X[:, (j-1)*nx+i] .= phi[:, i, j]
        Y[:, (j-1)*nx+i] .= α[:, i, j]
    end
end

using Flux
#data = [(X, Y)]
data = Flux.Data.DataLoader(X, Y, batchsize = 1, shuffle = true)

using Flux: @epochs
```nn = Chain(Dense(4, 16, tanh), Dense(16, 16, tanh), Dense(16, 16, tanh), Dense(16, 16, tanh), Dense(16, 4))```
nn = Chain(
    Dense(4, 16, relu),
    SkipConnection(Chain(Dense(16, 16, relu), Dense(16, 16, relu)), +),
    Dense(16, 4, relu),
)
ps = params(nn)

function loss(x, y)
    sum(abs2, nn(x) - y)
end
cb() = @show(loss(X, Y))
@epochs 10 Flux.train!(loss, ps, data, ADAM(), cb = Flux.throttle(cb, 1))

struct Shortcut{T}
    f::T
    σ::Function
end
Flux.@functor Shortcut
(m::Shortcut)(x) = m.σ.(m.f(x) .+ x)

nn = Chain(
    Dense(4, 16, tanh),
    Shortcut(Chain(Dense(16, 16, relu), Dense(16, 16, relu)), tanh),
    Dense(16, 4, tanh),
)
ps = params(nn)

using DiffEqFlux
ann = FastChain(
    FastDense(4, 16, tanh),
    FastDense(16, 16, tanh),
    FastDense(16, 16, tanh),
    FastDense(16, 4),
)
p = initial_params(ann)
function loss2(θ)
    sum(abs2, ann(X, θ) - Y)
end
callback = function (p, l)
    println("loss: $l")
    return false
end

res = DiffEqFlux.sciml_train(loss2, p, ADAM(), cb = callback, maxiters = 200)
res = DiffEqFlux.sciml_train(
    loss2,
    res.minimizer,
    ADAM(),
    cb = Flux.throttle(callback, 1),
    maxiters = 5000,
)

p, re = Flux.destructure(nn)
function loss3(θ)
    sum(abs2, re(θ)(X) - Y)
end
res = DiffEqFlux.sciml_train(loss3, p, ADAM(), cb = callback, maxiters = 200)


text:
to_func_name(x::Symbol) = string(x) |> lowercase |> Symbol

"""
     default_plot_signatures(funcname, PlotType)
Creates all the different overloads for `funcname` that need to be supported for the plotting frontend!
Since we add all these signatures to different functions, we make it reusable with this function
"""
function default_plot_signatures(funcname, funcname!, PlotType)
    quote
        """
            $($(funcname))(args...; attributes...)

        Command works on plot args 1:N and accepts keyword arguments to style the plot. Creates a new scene!
        """
        ($funcname)(args...; attributes...) = plot!(Scene(), $PlotType, Attributes(attributes), args...)

        """
            $($(funcname!))(args...; attributes...)

        Command works on plot args 1:N and accepts keyword arguments to style the plot. Adds new plot to `current_scene()`
        """
        ($funcname!)(args...; attributes...) = plot!(current_scene(), $PlotType, Attributes(attributes), args...)

        """
            $($(funcname!))(scene::SceneLike, args...; attributes...)

        Command works on plot args 1:N and accepts keyword arguments to style the plot. Adds new plot to `scene`!
        """
        ($funcname!)(scene::SceneLike, args...; attributes...) = plot!(scene, $PlotType, Attributes(attributes), args...)

        """
            $($(funcname))(attributes::Attributes, args...; attributes...)

        Like $($(funcname))(args...; attributes...) but accepts a theme as first argument. Creates a new scene!
        """
        ($funcname)(attributes::Attributes, args...; kw_attributes...) = plot!(Scene(), $PlotType, merge!(Attributes(kw_attributes), attributes), args...)

        """
            $($(funcname!))(attributes::Attributes, args...; attributes...)

        Like $($(funcname!))(args...; attributes...) but accepts a theme as first argument. Adds new plot to `current_scene()`!
        """
        ($funcname!)(attributes::Attributes, args...; kw_attributes...) = plot!(current_scene(), $PlotType, merge!(Attributes(kw_attributes), attributes), args...)

        """
            $($(funcname!))(attributes::Attributes, args...; attributes...)

        Like $($(funcname!))(scene, args...; attributes...) but accepts a theme as second argument. Adds new plot to `scene`!
        """
        ($funcname!)(scene::SceneLike, attributes::Attributes, args...; kw_attributes...) = plot!(scene, $PlotType, merge!(Attributes(kw_attributes), attributes), args...)
    end
end

"""

# type recipe are really simple and just overload the argument conversion pipeline.
convert_arguments(P::Type{<:AbstractPlot}, x::MyType) = convert_arguments(P, rand(10, 10))
# only apply this for a certain plot type:
convert_arguments(P::Type{<:Scatter}, x::MyType) = convert_arguments(P, rand(10, 10))
# optionally define the plotting type, when someone just says plot(x::MyType)
# - will fall back to whatever the standard plot type for what you return from argument_convert is!
plottype(::MyType) = Surface

# full recipes:

# (x, y, z) && themes are optional
@recipe(MyPlot, x, y, z) do scene
    Theme(
        plot_color => :red
    )
end
# this macro defines the following:

const Myplot{ArgTypes} = Combined{myplot, ArgTypes}
# the type parameter contains the function instead of e.g. a symbol. This way
# the mapping from Myplot to myplot is safer and simpler.
# Downside, we always need a function. Would there be rejection against that?

# all the signatures to make it nice to use:
myplot(args...; kw_args...) = ...
myplot!(scene, args...; kw_args...) = ...
myplot(kw_args::Dict, args...) = ...
myplot!(scene, kw_args::Dict, args...) = ...
#etc (not 100% settled what signatures there will be)

# the next one is optional, but it will allow to say plot_object[:x] to fetch the first argument
# from the call `plot_object = myplot(rand(10), rand(10), rand(10))`.
# If you leave out the (x, y, z)
# in the recipe macro, it will default to plot_object[1] (always works), and plot_object[:arg1]
argument_names(::Type{<: MyPlot}) = (:x, :y, :z)

# this function will insert the theme into any scene that plots Myplot.
function default_theme(scene, ::Myplot)
    Theme(
        plot_color => :red
    )
end
#-----------------------------
# Implement the recipe
function plot!(plot::MyPlot)
    # normal plotting code, building up on any previously defined recipes
    # or atomic plotting operations, and adds it to the combined `plot`:
    lines!(plot, rand(10), color = plot[:plot_color])
    plot!(plot, plot[:x], plot[:y])
    plot
end

# Add specialization
const MyVolume = MyPlot{Tuple{ArgTypes <: AbstractArray{<: AbstractFloat, 3}}}
argument_names(::Type{<: MyVolume}) = (:volume,) # again, optional
function plot!(plot::MyVolume)
    # plot a volume with a colormap going from fully transparent to plot_color
    volume!(plot, plot[:volume], colormap = :transparent => plot[:plot_color])
    plot
end
"""
macro recipe(theme_func, Tsym::Symbol, args::Symbol...)
    funcname_sym = to_func_name(Tsym)
    funcname! = esc(Symbol("$(funcname_sym)!"))
    PlotType = esc(Tsym)
    funcname = esc(funcname_sym)
    expr = quote
        $(funcname)() = not_implemented_for($funcname)
        const $(PlotType){$(esc(:ArgType))} = Combined{$funcname, $(esc(:ArgType))}
        Base.show(io::IO, ::Type{<: $PlotType}) = print(io, $(string(Tsym)), "{...}")
        $(default_plot_signatures(funcname, funcname!, PlotType))
        Base.@__doc__($funcname)
        AbstractPlotting.default_theme(scene, ::Type{<: $PlotType}) = $(esc(theme_func))(scene)
        export $PlotType, $funcname, $funcname!
    end
    if !isempty(args)
        push!(expr.args, :($(esc(:argument_names))(::Type{<: $PlotType}, len::Integer) = $args))
    end
    expr
end


text:
#=
Created on 05/01/2021 20:10:45
Last update: -

@author: Michiel Stock
michielfmstock@gmail.com

Example of some interesting matrices
- Strang matrix
=#

# STRANG MATRIX

struct Strang <: AbstractMatrix{Int}
    n::Int
end

Base.size(S::Strang) = (S.n, S.n)
Base.getindex(S::Strang, i, j) = i==j ? 2 : (abs(i - j) == 1 ?  -1 : 0)

S = Strang(1000)

@time sum(S)

Base.sum(S::Strang) = 2

@time sum(S)

v = randn(1000)

@time S * v

function Base.:*(S::Strang, v::Vector)
    n = length(v)
    @assert size(S, 2) == n
    x = similar(v)
    for i in 1:n
        x[i] = v[i]
        i > 1 && (x[i] += v[i-1])
        i < n && (x[i] += v[i+1])
    end
    return x
end

@time S * v

# Vandermonde matrix

struct Vandermonde{T,VT} <: AbstractMatrix{T}
    α::VT
    m::Int
    Vandermonde(α::AbstractVector{T}, m) where {T} = new{T,typeof(α)}(α,m)
end

Vandermonde(α::Vector{<:Number}) = Vandermonde(α, length(α))

Base.size(V::Vandermonde) = (length(V.α), V.m)
Base.getindex(V::Vandermonde, i, j) = V.α[i]^(j-1)

α = [1, 2, 3, 4]

V = Vandermonde(α, 4)

using LinearAlgebra

LinearAlgebra.det(V::Vandermonde) = ((xi - xj) for (i, xi) in enumerate(V.α), (j, xj) in enumerate(V.α) if i < j) |> prod


text:
__precompile__()
module FiniteVolumeMesh

export HomogeneousMesh, Vec, Vec2D, TriangleCell, read_input, CellProblemAdvecTemp, PCG

include("utils.jl")

abstract type AbstractVec{T} <: DenseArray{T,1} end
abstract type AbstractTen{T} <: DenseArray{T,2} end

abstract type AbstractCell2Node{NN,NF} <: NTupleWrap{NN,UInt} end
abstract type AbstractCell2Node2D{NN,NF} <: AbstractCell2Node{NN,NF} end

abstract type AbstractFace2Node{N} <: NTupleWrap{N,UInt} end
const AbstractFace2Node1D = AbstractFace2Node{1}
const AbstractFace2Node2D = AbstractFace2Node{2}

abstract type AbstractMesh end

include("KrylovCtrl.jl")
include("pcg.jl")
include("vec.jl")
include("vecarray.jl")
include("tensors.jl")
include("tensorsarray.jl")
include("fieldtypes.jl")
include("geometry.jl")
include("cells.jl")
include("faces.jl")
include("boundaryconditions.jl")
include("mesh.jl")
include("MatrixLike.jl")
include("advection.jl")
include("input.jl")
include("readneutral.jl")
include("writevtk.jl")
include("derivatives.jl")
include("solvers.jl")
include("laplacian.jl")

end # module


text:
using CompScienceMeshes
using Test

for U in [Float32, Float64]
    V = [
        point(U,0,0,0),
        point(U,1,0,0),
        point(U,1,1,0),
        point(U,0,1,0)]

    F = [
        index(1,2,3),
        index(1,3,4)]

    m = Mesh(V,F)

    f = CompScienceMeshes.FlippedMesh(m)
    g = -f
    @test g === m

    T = [index(2,1,3), index(3,1,4)]
    for (i,_c) in enumerate(cells(f))
        @test _c == T[i]
    end
end
# @test f.faces[1] == index(2,1,3)
# @test f.faces[2] == index(3,1,4)

max_stars_repo_name: krcools/CompScienceMeshes

max_issues_repo_name: krcools/CompScienceMeshes

max_forks_repo_name: krcools/CompScienceMeshes


text:
using Tk

window = Toplevel("Hello World", 200, 100, false)
pack_stop_propagate(window)

fr = Frame(window)
pack(fr, expand=true, fill="both")

txt = Label(fr, "Hello World")
pack(txt, expand=true)

set_visible(window, true)

# sleep(7)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using LinearAlgebra
using MAT
using AdFem
using PyPlot
using SparseArrays

# geometry setup in domain [0,1]^2
solid_left = 0.45
solid_right = 0.55
solid_top = 0.5
solid_bottom = 0.52

chip_left = 0.48
chip_right = 0.52
chip_top = 0.5
chip_bottom = 0.505

k_mold = 0.014531
k_chip = Variable(1.0) #2.60475
k_air = 0.64357
nu = 0.47893 # equal to 1/Re
power_source = 82.46295 #82.46295 = 1.0e6 divide by air rho cp   #0.0619 = 1.0e6 divide by chip die rho cp
buoyance_coef = 299102.83

u_std = 0.001
p_std = 0.000001225
T_infty = 300

m = 200
n = 200
h = 1/n
NT = 5    # number of iterations for Newton's method, 8 is good for m=400

# compute solid indices and chip indices
solid_fem_idx = Array{Int64, 1}([])
solid_fvm_idx = Array{Int64, 1}([])
chip_fem_idx = Array{Int64, 1}([])
chip_fvm_idx = Array{Int64, 1}([])

for i = 1:(m+1)
    for j = 1:(n+1)
        if (i-1)*h >= solid_left-1e-9 && (i-1)*h <= solid_right+1e-9 && (j-1)*h >= solid_top-1e-9 && (j-1)*h <= solid_bottom+1e-9
            # print(i, j)
            global solid_fem_idx = [solid_fem_idx; (j-1)*(m+1)+i]
            if (i-1)*h >= chip_left-1e-9 && (i-1)*h <= chip_right+1e-9 && (j-1)*h >= chip_top-1e-9 && (j-1)*h <= chip_bottom+1e-9
                global chip_fem_idx = [chip_fem_idx; (j-1)*(m+1)+i]
            end
        end
    end
end

for i = 1:m
    for j = 1:n
        if (i-1)*h + h/2 >= solid_left-1e-9 && (i-1)*h + h/2 <= solid_right+1e-9 &&
            (j-1)*h + h/2 >= solid_top-1e-9 && (j-1)*h + h/2 <= solid_bottom+1e-9
            global solid_fvm_idx = [solid_fvm_idx; (j-1)*m+i]
            if (i-1)*h + h/2 >= chip_left-1e-9 && (i-1)*h + h/2 <= chip_right+1e-9 && (j-1)*h + h/2 >= chip_top-1e-9 && (j-1)*h + h/2<= chip_bottom+1e-9
                global chip_fvm_idx = [chip_fvm_idx; (j-1)*m+i]
            end
        end
    end
end

# initialize space varying k and heat source
k_fem = k_air * constant(ones((m+1)*(n+1)))
k_fem = scatter_update(k_fem, solid_fem_idx, k_mold * ones(length(solid_fem_idx)))
k_fem = scatter_update(k_fem, chip_fem_idx, k_chip * ones(length(chip_fem_idx)))
kgauss = fem_to_gauss_points(k_fem, m, n, h)

heat_source_fem = zeros((m+1)*(n+1))
heat_source_fem[chip_fem_idx] .= power_source #/ h^2
heat_source_gauss = fem_to_gauss_points(heat_source_fem, m, n, h)

# chip_gauss_idx = [ 4 .* chip_fvm_idx; 4 .* chip_fvm_idx .- 1; 4 .* chip_fvm_idx .- 2; 4 .* chip_fvm_idx .- 3]
# heat_source_gauss = zeros(4*m*n)
# heat_source_gauss[chip_gauss_idx] .= power_source

# F1 = compute_fem_source_term1(eval_f_on_gauss_pts(ffunc_, m, n, h), m, n, h)
# F2 = compute_fem_source_term1(eval_f_on_gauss_pts(gfunc_, m, n, h), m, n, h)
# H = h^2*eval_f_on_fvm_pts(hfunc_, m, n, h)
B = constant(compute_interaction_matrix(m, n, h))

# compute F
Laplace = nu * constant(compute_fem_laplace_matrix1(m, n, h))
heat_source = constant(compute_fem_source_term1(heat_source_gauss, m, n, h))

LaplaceK = constant(compute_fem_laplace_matrix1(kgauss, m, n, h))
# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k = @. k_nn(x, y); k=stack(k)
# kgauss = fem_to_gauss_points(k, m, n, h)
# LaplaceK = compute_fem_laplace_matrix1(kgauss, m, n, h)

bd = bcnode("all", m, n, h)

# only apply Dirichlet to velocity; set left bottom two points to zero to fix rank deficient problem for pressure

bd = [bd; bd .+ (m+1)*(n+1);
     2*(m+1)*(n+1)+1; 2*(m+1)*(n+1)+m;
    #  (2*(m+1)*(n+1)+m*n )+1:(2*(m+1)*(n+1)+m*n )+m+1]
     bd .+ (2*(m+1)*(n+1)+m*n )]

# add solid region into boundary condition for u, v, p
bd = [bd; solid_fem_idx; solid_fem_idx .+ (m+1)*(n+1); solid_fvm_idx .+ 2(m+1)*(n+1)]

function compute_residual(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]
    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    interaction = compute_interaction_term(p, m, n, h) # julia kernel needed
    f1 = compute_fem_source_term1(ugauss.*ux, m, n, h)
    f2 = compute_fem_source_term1(vgauss.*uy, m, n, h)
    f3 = -interaction[1:(m+1)*(n+1)]
    f4 = Laplace*u
    # f5 = -F1
    F = f1 + f2 + f3 + f4 #+ f5

    g1 = compute_fem_source_term1(ugauss.*vx, m, n, h)
    g2 = compute_fem_source_term1(vgauss.*vy, m, n, h)
    g3 = -interaction[(m+1)*(n+1)+1:end]
    g4 = Laplace*v
    # g5 = -F2
    T_gauss = fem_to_gauss_points(T, m, n, h)
    buoyance_term = - buoyance_coef * compute_fem_source_term1(T_gauss, m, n, h)

    G = g1 + g2 + g3 + g4 + buoyance_term #+ g5

    H0 = -B * [u;v] # + H

    T0 = LaplaceK * T + compute_fem_advection_matrix1(ugauss,vgauss, m, n, h) * T - heat_source
    R = [F;G;H0;T0]
    return R
end

function compute_jacobian(S)
    u, v, p, T = S[1:(m+1)*(n+1)],
        S[(m+1)*(n+1)+1:2(m+1)*(n+1)],
        S[2(m+1)*(n+1)+1:2(m+1)*(n+1)+m*n],
        S[2(m+1)*(n+1)+m*n+1:end]

    G = eval_grad_on_gauss_pts([u;v], m, n, h)
    ugauss = fem_to_gauss_points(u, m, n, h)
    vgauss = fem_to_gauss_points(v, m, n, h)
    ux, uy, vx, vy = G[:,1,1], G[:,1,2], G[:,2,1], G[:,2,2]

    M1 = constant(compute_fem_mass_matrix1(ux, m, n, h))
    M2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h)) # a julia kernel needed
    M3 = Laplace
    Fu = M1 + M2 + M3

    Fv = constant(compute_fem_mass_matrix1(uy, m, n, h))

    N1 = constant(compute_fem_mass_matrix1(vy, m, n, h))
    N2 = constant(compute_fem_advection_matrix1(constant(ugauss), constant(vgauss), m, n, h))
    N3 = Laplace
    Gv = N1 + N2 + N3

    Gu = constant(compute_fem_mass_matrix1(vx, m, n, h))

    M = LaplaceK + constant(compute_fem_advection_matrix1(ugauss,vgauss, m, n, h))

    gradT = eval_grad_on_gauss_pts1(T, m, n, h)
    Tx, Ty = gradT[:,1], gradT[:,2]
    DU_TX = constant(compute_fem_mass_matrix1(Tx, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)
    DV_TY = constant(compute_fem_mass_matrix1(Ty, m, n, h))       # (m+1)*(n+1), (m+1)*(n+1)

    T_mat = -buoyance_coef * constant(compute_fem_mass_matrix1(m, n, h))
    T_mat = [SparseTensor(spzeros((m+1)*(n+1), (m+1)*(n+1))); T_mat]

    J0 = [Fu Fv
          Gu Gv]

    J1 = [J0 -B' T_mat
        -B spdiag(zeros(size(B,1))) SparseTensor(spzeros(m*n, (m+1)*(n+1)))]

    N = 2*(m+1)*(n+1) + m*n
    J = [J1
        [DU_TX DV_TY SparseTensor(spzeros((m+1)*(n+1), m*n)) M]]
end

function solve_steady_cavityflow_one_step(S)
    residual = compute_residual(S)
    J = compute_jacobian(S)

    J, _ = fem_impose_Dirichlet_boundary_condition1(J, bd, m, n, h)
    residual = scatter_update(residual, bd, zeros(length(bd)))    # residual[bd] .= 0.0 in Tensorflow syntax

    d = J\residual
    residual_norm = norm(residual)
    op = tf.print("residual norm", residual_norm)
    d = bind(d, op)
    S_new = S - d
    return S_new
end

function condition(i, S_arr)
    i <= NT + 1
end

function body(i, S_arr)
    S = read(S_arr, i-1)
    op = tf.print("i=",i)
    i = bind(i, op)
    S_new = solve_steady_cavityflow_one_step(S)
    S_arr = write(S_arr, i, S_new)
    return i+1, S_arr
end

S_arr = TensorArray(NT+1)
S_arr = write(S_arr, 1, zeros(m*n+3*(m+1)*(n+1)))

i = constant(2, dtype=Int32)

_, S = while_loop(condition, body, [i, S_arr])
S = set_shape(stack(S), (NT+1, 2*(m+1)*(n+1)+m*n+(m+1)*(n+1)))

# sess = Session(); init(sess)
# output = run(sess, S)
V_computed = S[end, :]

V_data = matread("xzchipv0_data.mat")["V"]

sample_size = 10
idx = rand(1:(m+1)*(n+1), sample_size)
idx = [idx; (m+1)*(n+1) .+ idx; 2*(m+1)*(n+1)+m*n .+ idx] # observe velocity and temperature
observed_data = V_data[idx]

noise = false
noise_level = 0.05
if noise
    noise_ratio = (1 - noise_level) .+ 2 * noise_level * rand(Float64, size(observed_data)) # uniform on (1-noise_level, 1+noise_level)
    observed_data = observed_data .* noise_ratio
end

loss = mean((V_computed[idx] .- observed_data)^2)
loss = loss * 1e10
# ---------------------------------------------------
# create a session and run
max_iter = 10
sess = Session(); init(sess)
loss_ = BFGS!(sess, loss, max_iter)
figure(); semilogy(loss_); savefig("xzchipv0_loss.png")

# xy = fem_nodes(m, n, h)
# x, y = xy[:, 1], xy[:, 2]
# k0 = @. k_exact(x,y)

# figure(figsize=(14,4));
# subplot(131)
# visualize_scalar_on_fem_points(k0, m, n, h); title("conductivity exact")
# subplot(132)
# visualize_scalar_on_fem_points(run(sess, k), m, n, h); title("conductivity prediction")
# subplot(133)
# visualize_scalar_on_fem_points(k0.-run(sess, k), m, n, h); title("conductivity difference")
# savefig("xzchipv0_k.png")


text:
## Network for representing continous functions (value or policy)
mutable struct LatentConditionedNetwork <: NetworkPolicy
    policy
    z
end

device(π::LatentConditionedNetwork) = device(π.policy)

Flux.@functor LatentConditionedNetwork

Flux.trainable(π::LatentConditionedNetwork) = Flux.trainable(π.policy)

layers(π::LatentConditionedNetwork) = layers(π.policy)

POMDPs.action(π::LatentConditionedNetwork, s; z=π.z) = value(π, s; z=z)

function POMDPs.value(π::LatentConditionedNetwork, s; z=π.z)
    if size(z, 2) != size(s)[end]
        z = repeat(z, 1, ndims(s) == 1 ? 1 : size(s)[end])
    end
    value(π.policy, vcat(z,s))
end

function POMDPs.value(π::LatentConditionedNetwork, s, a; z=π.z)
    if size(z, 2) != size(s)[end]
        z = repeat(z, 1, ndims(s) == 1 ? 1 : size(s)[end])
    end
    value(π.policy, vcat(z,s), a)
end

action_space(π::LatentConditionedNetwork) = action_space(π.policy)

@with_kw mutable struct OffPolicyLatentSolver <: Solver
    π # Policy
    S::AbstractSpace # State space
    A::AbstractSpace = action_space(π) # Action space
    N::Int = 1000 # Number of environment interactions
    ΔN::Int = 4 # Number of interactions between updates
    max_steps::Int = 100 # Maximum number of steps per episode
    log::Union{Nothing, LoggerParams} = nothing # The logging parameters
    i::Int = 0 # The current number of environment interactions
    a_opt::Union{Nothing, TrainingParams} = nothing # Training parameters for the actor
    c_opt::TrainingParams # Training parameters for the critic
    post_batch_callback = (𝒟) -> nothing

    # Off-policy-specific parameters
    π⁻ = deepcopy(π)
    π_explore::Policy # exploration noise
    target_update = (π⁻, π; kwargs...) -> polyak_average!(π⁻, π, 0.005f0) # Function for updating the target network
    target_fn # Target for critic regression with input signature (π⁻, 𝒟, γ; i)
    buffer_size = 1000 # Size of the buffer
    required_columns = Symbol[]
    buffer::ExperienceBuffer = ExperienceBuffer(S, A, buffer_size, required_columns) # The replay buffer
    buffer_init::Int = max(c_opt.batch_size, 200) # Number of observations to initialize the buffer with
    extra_buffers = [] # extra buffers (i.e. for experience replay in continual learning)
    buffer_fractions = [1.0] # Fraction of the minibatch devoted to each buffer
    z_dists = MvNormal[]
end

TIER(;π::ActorCritic, ΔN=50, π_explore=GaussianNoiseExplorationPolicy(0.1f0),  a_opt::NamedTuple=(;), c_opt::NamedTuple=(;), log::NamedTuple=(;), π_smooth::Policy=GaussianNoiseExplorationPolicy(0.1f0, ϵ_min=-0.5f0, ϵ_max=0.5f0), kwargs...) =
    OffPolicyLatentSolver(;
        π=π,
        ΔN=ΔN,
        log=LoggerParams(;dir = "log/ddpg", log...),
        a_opt=TrainingParams(;loss=TD3_actor_loss_w_latent, name="actor_", a_opt...),
        c_opt=TrainingParams(;loss=double_Q_loss_w_latent, name="critic_", epochs=ΔN, c_opt...),
        π_explore=π_explore,
        target_fn=TD3_target_w_latent(π_smooth),
        kwargs...)

# function optimize_latent(loss, latent_dim)
#     rng = MersenneTwister(0)
#     z_prospective = [Float32.(rand(rng, Uniform(-1,1), latent_dim)) for i=1:100]
#     vals = [loss(z) for z in z_prospective]
#     z_prospective[argmin(vals)]
# end

function cross_entropy_optimization(f, P, latent_dim)
    m = max(floor(Int, 1000 * norm(P.Σ)), 5*(latent_dim+3))
    m_elite = floor(Int, 0.2*m)
    # P = MvNormal(zeros(latent_dim), 0.5*I)
    for i=1:1
    # mi = max(floor(Int, m / length(Ps)), 1)
    # samples = clamp.(hcat([rand(P,mi) for P in Ps]...), -1f0, 1f0)
        # P = MvNormal(P.μ, P.Σ + 1f-5*I)
        samples = clamp.(rand(P, m), -1f0, 1f0)
        samples .+= rand(Uniform(-1f-5, 1f-5), size(samples)...)
        order = sortperm([f(samples[:,i]) for i=1:m])
        P = fit(MvNormal, Float64.(samples[:, order[1:m_elite]]))
    end
    P
end

function find_latent!(loss, P, latent_dim)
    # 𝒟[:μ_z] .= optimize_latent(loss, latent_dim)
    # Ps = [MvNormal(𝒟[:μ_z][:,i], 𝒟[:Σ_z][:,:,i] + 1f-5 * I ) for i=1:length(𝒟)]
    P = cross_entropy_optimization(loss, P, latent_dim)
    # 𝒟[:μ_z] .= P.μ
    # P
    # 𝒟[:Σ_z] .= P.Σ
end

function latent_loss(π⁻, π, γ, 𝒟)
    (z) -> begin
        y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* value(critic(π⁻), 𝒟[:sp], action(actor(π⁻), 𝒟[:sp], z=z), z=z)
        Q = value(critic(π), 𝒟[:s], 𝒟[:a], z=z)
        Flux.mae(Q, y)
    end
end

function TD3_latent_loss(π⁻, π, γ, 𝒟)
    (z) -> begin
        y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* value(critic(π).N1, 𝒟[:sp], action(actor(π), 𝒟[:sp], z=z), z=z)
        Q = value(critic(π).N1, 𝒟[:s], 𝒟[:a], z=z)
        Flux.mae(Q, y)
    end
end

function latent_target(π, 𝒟, γ; z=𝒟[:μ_z], kwargs...)
    𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* value(critic(π), 𝒟[:sp], action(actor(π), 𝒟[:sp], z=z), z=z)
end

function critic_loss_w_latent(π, 𝒟, y; loss=Flux.mae, weighted=false, name=:Qavg, info=Dict())
    Q = value(critic(π), 𝒟[:s], 𝒟[:a]; z=𝒟[:μ_z])

    # Store useful information
    ignore() do
        info[name] = mean(Q)
    end

    loss(Q, y, agg = weighted ? weighted_mean(𝒟[:weight]) : mean)
end

function double_Q_loss_w_latent(π, 𝒟, y; info=Dict(), weighted=false)
    q1loss = critic_loss_w_latent(π.C.N1, 𝒟, y, info=info, name=:Q1avg, weighted=weighted)
    q2loss = critic_loss_w_latent(π.C.N2, 𝒟, y, info=info, name=:Q2avg, weighted=weighted)
    q1loss + q2loss
end

function actor_loss_w_latent(π, 𝒟; info=Dict())
    -mean(value(critic(π), 𝒟[:s], action(actor(π), 𝒟[:s], z=𝒟[:μ_z]), z=𝒟[:μ_z]))
end

TD3_actor_loss_w_latent(π, 𝒟; info = Dict()) = -mean(value(critic(π).N1, 𝒟[:s], action(actor(π), 𝒟[:s],z=𝒟[:μ_z]), z=𝒟[:μ_z]))

function TD3_target_w_latent(π_smooth)
    (π, 𝒟, γ::Float32; i,  z=𝒟[:μ_z]) -> begin
        ap, _ = exploration(π_smooth, 𝒟[:sp], π_on=π, i=i)
        y = 𝒟[:r] .+ γ .* (1.f0 .- 𝒟[:done]) .* min.(value(critic(π).N1, 𝒟[:sp], ap, z=z), value(critic(π).N2, 𝒟[:sp], ap, z=z))
    end
end

action_regularization_tier(π, 𝒟s) = length(𝒟s) == 0 ? 0 : mean([Flux.mse(action(actor(π), 𝒟[:s], z=𝒟[:μ_z]), 𝒟[:a]) for 𝒟 in 𝒟s])
action_value_regularization_tier(π, 𝒟s) = length(𝒟s) == 0 ? 0 : mean([Flux.mse(value(critic(π).N1, 𝒟[:s], 𝒟[:a], z=𝒟[:μ_z]), 𝒟[:value]) for 𝒟 in 𝒟s]) +  mean([Flux.mse(value(critic(π).N2, 𝒟[:s], 𝒟[:a], z=𝒟[:μ_z]), 𝒟[:value]) for 𝒟 in 𝒟s])

function POMDPs.solve(𝒮::OffPolicyLatentSolver, mdp)
    # Compute the latent dimension
    latent_dim = length(actor(𝒮.π).z)

    # Add data for normal distributions
    # C = capacity(𝒮.buffer)
    # 𝒮.buffer.data[:μ_z] = zeros(Float32, latent_dim, C)
    # 𝒮.buffer.data[:Σ_z] = 0.5f0*repeat(Array{Float32, 2}(I, latent_dim, latent_dim), outer=[1,1,C])

    # Create minibatch buffers for each buffer
    allbuffs = [𝒮.extra_buffers..., 𝒮.buffer]
    push!(𝒮.z_dists, MvNormal(zeros(latent_dim), 0.5*I))
    #
    # for b in allbuffs
    #     @assert haskey(b, :μ_z) && haskey(b, :Σ_z)
    # end

    batches = split_batches(𝒮.c_opt.batch_size, 𝒮.buffer_fractions)
    𝒟s = [buffer_like(b, capacity=batchsize, device=device(𝒮.π)) for (b, batchsize) in zip(allbuffs, batches)]
    # Add latent dimension
    # last_z = [zeros(Float32, latent_dim) for _ in 𝒟s]

    for 𝒟 in 𝒟s
        𝒟.data[:μ_z] = zeros(Float32, latent_dim, capacity(𝒟))
        # 𝒟.data[:Σ_z] = 0.5f0*repeat(Array{Float32, 2}(I, latent_dim, latent_dim), outer=[1,1,capacity(𝒟)])
        𝒟.data[:value] = zeros(Float32, 1, capacity(𝒟))
    end

    γ = Float32(discount(mdp))
    s = Sampler(mdp, 𝒮.π, S=𝒮.S, A=𝒮.A, max_steps=𝒮.max_steps, π_explore=𝒮.π_explore, required_columns=extra_columns(𝒮.buffer))
    isnothing(𝒮.log.sampler) && (𝒮.log.sampler = s)

    # Log the pre-train performance
    log(𝒮.log, 𝒮.i)

    # Fill the buffer with initial observations before training
    𝒮.i += fillto!(𝒮.buffer, s, 𝒮.buffer_init, i=𝒮.i, explore=true)

    # Loop over the desired number of environment interactions
    for 𝒮.i in range(𝒮.i, stop=𝒮.i + 𝒮.N - 𝒮.ΔN, step=𝒮.ΔN)
        # Sample transitions into the replay buffer
        data = steps!(s, Nsteps=𝒮.ΔN, explore=true, i=𝒮.i)
        # data[:μ_z] = repeat(actor(𝒮.π).z, outer=[1,𝒮.ΔN])
        # data[:Σ_z] = repeat(Array{Float32, 2}(I, latent_dim, latent_dim), outer=[1,1,𝒮.ΔN])
        push!(𝒮.buffer, data)

        # callback for potentially updating the buffer
        𝒮.post_batch_callback(𝒮.buffer)

        # Determine the latent variables
        info_z = Dict()
        for (b, index, Dz) in zip(allbuffs, 1:length(𝒟s), 𝒮.z_dists)
            D = minibatch(b, rand(1:length(b), 1000))
            z_before = 𝒮.z_dists[index].μ
            𝒮.z_dists[index] = find_latent!(TD3_latent_loss(𝒮.π⁻, 𝒮.π, γ, D), Dz, latent_dim)
            # b[:μ_z][:,𝒟.indices] .= 𝒟[:μ_z]
            # b[:Σ_z][:,:,𝒟.indices] .= 𝒟[:Σ_z]
            z_after = 𝒮.z_dists[index].μ
            info_z["Δz$index"] = norm(z_after .- z_before) # Store the change in task identifier
            info_z["z$index"] = norm(z_after) # Store the change in task identifier
            info_z["Σ$index"] = norm(𝒮.z_dists[index].Σ) # Store the norm of the covariance
            info_z["td_error$index"] = mean(abs.(value(critic(𝒮.π).N1, D[:s], D[:a], z=z_after)  .- 𝒮.target_fn(𝒮.π⁻, D, γ, i=𝒮.i, z=z_after))) # store td error on the ith task
        end

        infos = []
        # Loop over the desired number of training steps
        for epoch in 1:𝒮.c_opt.epochs
            # Sample a random minibatch of 𝑁 transitions (sᵢ, aᵢ, rᵢ, sᵢ₊₁) from 𝒟
            for (𝒟, b, index) in zip(𝒟s, allbuffs, 1:length(𝒟s))
                rand!(𝒟, b, i=𝒮.i) # sample a batch
                𝒟[:μ_z] .= 𝒮.z_dists[index].μ # set the parameters
            end

            # Set the latent variable of the current actor and critic
            z = 𝒮.z_dists[end].μ
            actor(𝒮.π).z .= z
            critic(𝒮.π).N1.z .= z
            critic(𝒮.π).N2.z .= z

            # concatenate the minibatch buffers
            𝒟 = hcat(𝒟s...)

            # Compute target
            y = 𝒮.target_fn(𝒮.π⁻, 𝒟, γ, i=𝒮.i)

            # Train the critic
            info = train!(critic(𝒮.π), (;kwargs...) -> 𝒮.c_opt.loss(𝒮.π, 𝒟, y; kwargs...) + action_value_regularization_tier(𝒮.π, 𝒟s[1:end-1]), 𝒮.c_opt)

            # Train the actor
            if !isnothing(𝒮.a_opt) && ((epoch-1) % 𝒮.a_opt.update_every) == 0
                info_a = train!(actor(𝒮.π), (;kwargs...) -> 𝒮.a_opt.loss(𝒮.π, 𝒟; kwargs...) + action_regularization_tier(𝒮.π, 𝒟s[1:end-1]), 𝒮.a_opt)
                info = merge(info, info_a)

                # Update the target network
                𝒮.target_update(𝒮.π⁻, 𝒮.π)
            end

            # Store the training information
            push!(infos, info)

        end
        # If not using a separate actor, update target networks after critic training
        isnothing(𝒮.a_opt) && 𝒮.target_update(𝒮.π⁻, 𝒮.π, i=𝒮.i + 1:𝒮.i + 𝒮.ΔN)

        # Log the results
        log(𝒮.log, 𝒮.i + 1:𝒮.i + 𝒮.ΔN, aggregate_info(infos), info_z)
    end
    𝒮.i += 𝒮.ΔN
    𝒮.π
end

# function optimize_latent(π, loss, i, kmax=10)
#     z, y = π.z, loss()
#     z_best, y_best = π.z, y
#     scale = Float32((1/5000)*i + 1)
#     for k in 1:kmax
#         π.z = min.(max.(z .+ Float32.(randn(size(z)...)) ./ scale, -1f0), 1f0)
#         y′ = loss()
#         Δy = y′ - y
#         if Δy ≤ 0 || rand() < exp(-Δy*i)
#             z, y = π.z, y′
#         end
#         if y′ < y_best
#             z_best, y_best = π.z, y′
#         end
#     end
#     π.z = z_best
# end


text:
module TestJoin

using Base.Test
using TimeData
using DataArrays
using Dates

using DataFrames

println("Running join tests")

##########################
## find dates functions ##
##########################

idx1 = [5, 6, 8, 9, 10]
idx2 = [2, 3, 4, 5, 6, 7, 9]

## with numeric indices
equalIds = TimeData.findAllOccurringDates(idx1, idx2)
ids = TimeData.getIndices(equalIds, idx1, idx2)
@test ids == sort(unique([idx1, idx2]))

## with date indices
dats1 = [Date(2000,1, ii) for ii=idx1]
dats2 = [Date(2000,1, ii) for ii=idx2]
equalIds = TimeData.findAllOccurringDates(idx1, idx2)
ids = TimeData.getIndices(equalIds, dats1, dats2)
@test ids == sort(unique([dats1, dats2]))

## left join with numeric indices
equalIds = TimeData.findDatesOfLeftInRight(idx1, idx2)

vals = TimeData.getIndices(equalIds, idx2) |>
ids -> DataFrame(x = ids) |>
complete_cases! |>
array
@test vals[:] == sort(unique(intersect(idx1, idx2)))

## left join with dates
equalIds = TimeData.findDatesOfLeftInRight(dats1, dats2)

vals = TimeData.getIndices(equalIds, dats2) |>
inds -> DataFrame(x = inds) |>
complete_cases! |>
array
@test vals[:] == sort(unique(intersect(dats1, dats2)))

## inner join with numeric indices
equalIds = TimeData.findCommonDates(idx1, idx2)
vals = TimeData.getIndices(equalIds, idx1, idx2)
@test vals == sort(unique(intersect(idx1, idx2)))

## join for numeric index
tm1 = Timematr(rand(length(idx1), 2), idx1)
tm2 = Timematr(rand(length(idx2), 3), idx2)

joinSortedIdx_left(tm1, tm2)
joinSortedIdx_right(tm1, tm2)
joinSortedIdx_inner(tm1, tm2)
joinSortedIdx_outer(tm1, tm2)

## join for date index
tm1 = Timematr(rand(length(idx1), 2), dats1)
tm2 = Timematr(rand(length(idx2), 3), dats2)

joinSortedIdx_left(tm1, tm2)
joinSortedIdx_right(tm1, tm2)
joinSortedIdx_inner(tm1, tm2)
joinSortedIdx_outer(tm1, tm2)

## outer join old bug
##-------------------

tm1 = Timenum(DataFrame(Adj_Close = NA), [Date(2012,1,3)])
tm2 = Timenum(DataFrame(Adj_Close2 = [3, 4, 5, 2, 1]),
              [Date(2012,1,1):Date(2012,1,5)])
actOut = joinSortedIdx_outer(tm1, tm2)

expOut = Timedata(DataFrame(Adj_Close = DataArray(NAtype, 5),
                           Adj_Close2 = [3, 4, 5, 2, 1]),
              [Date(2012,1,1):Date(2012,1,5)])
@test isequal(actOut, expOut)

end


text:
 using Base

mutable struct LimitOrder
    price::UInt128
    size::UInt128
    trader_id::UInt128
    side::Bool
    order_id::UInt128
end

function Base.size(o::LimitOrder)::UInt128
    return o.size
end

function trader(o::LimitOrder)::UInt128
    return o.trader_id
end

function side(o::LimitOrder)::Bool
    return o.side
end

function id(o::LimitOrder)::UInt128
    return o.order_id
end


text:
using Random, Test
using CollegeEntry

mdl = CollegeEntry;

function max_choice_test()
    rng = MersenneTwister(12);
    vWork = 1.0;
    vCollege_cV = [0.5, 3.0];
    v, ic = mdl.max_choice(vWork, vCollege_cV, [true, true]);
    @test v == 3.0  &&  ic == 2
    v, ic = mdl.max_choice(vWork, vCollege_cV, [true, false]);
    @test v == vWork  &&  ic == 0
    v, ic = mdl.max_choice(vWork, vCollege_cV, [false, false]);
    @test v == vWork  &&  ic == 0

    J = 5;
    nc = 7;
    vWorkV = randn(rng, J);
    vCollege_jcM = rand(rng, J, nc);
    admitV = rand(rng, Bool, nc);
    vV, icV = mdl.max_choices(vWorkV, vCollege_jcM, admitV);
    for j = 1 : J
        v, ic = mdl.max_choice(vWorkV[j], vec(vCollege_jcM[j,:]), admitV);
        @test isapprox(vV[j], v)  &&  isequal(icV[j], ic)
    end
end

function helpers_test()
    @testset "Helpers" begin
        x = rand(4,3,2) ./ 5.0;
        x[2 : 2 : 8] .= -0.000000001;
        x ./= sum(x);
        x .+= 0.000000001;
        mdl.make_valid_probs!(x);
        @test sum(x) <= 1.0
        @test all(x .>= 0.0)
	end
end

function find_class_test()
    @testset "find_class" begin
        cutoffV = [0.1, 0.5, 0.9];
        @test mdl.find_class(0.09, cutoffV) == 1
        @test mdl.find_class(0.11, cutoffV) == 2
        @test mdl.find_class(0.49, cutoffV) == 2
        @test mdl.find_class(0.51, cutoffV) == 3
        @test mdl.find_class(0.94, cutoffV) == 4
    end
end

@testset "Helpers" begin
    max_choice_test()
    helpers_test();
    find_class_test();
end

# -------------


text:
# global constants, to be also imported by parallel processes

# timeout = 1000000
timeout = 240


text:
#
# Date created: 2022-03-09
# Author: aradclif
#
#
############################################################################################

function aminusb_exp_sumbody(N::Int, D)
    body = Expr(:block)
    params = D.parameters
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    e = Expr(:(=), c, Expr(:call, :+, c, Expr(:call, :exp, Expr(:call, :-, a, b))))
    push!(body.args, e)
    body
end

function aminusb_exp_sum_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusb_exp_sumbody(N, D)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return C
    end
end

@generated function aminusb_exp_sum!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    aminusb_exp_sum_quote(N, D)
end

function logself_plusb!(C::AbstractArray{Tₒ, N}, B::AbstractArray{T, N}) where {Tₒ, T, N}
    @turbo for i ∈ eachindex(C)
        C[i] = B[i] + log(C[i])
    end
    C
end

function lvlogsumexp(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        # return hvncat(ntuple(i -> 1, Val(N)), true, lvlogsumexp1(A))
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvlogsumexp1(A))
    else
        B = lvmaximum(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        C = zeros(Base.promote_op(exp, T), Dᴮ)
        aminusb_exp_sum!(C, A, B, Dᴮ′)
        logself_plusb!(C, B)
        # return C
    end
    return C
end

lvlogsumexp(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvlogsumexp(A, (dims,))
lvlogsumexp(A::AbstractArray{T, N}; dims=:) where {T, N} = lvlogsumexp(A, dims)
lvlogsumexp(A::AbstractArray{T, N}) where {T, N} = lvlogsumexp1(A)
lvlogsumexp(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvlogsumexp1(A)

function lvlogsumexp1(A::AbstractArray{T, N}) where {T, N}
    α = typemin(T)
    s = zero(promote_type(T, Float64))
    @turbo for i ∈ eachindex(A)
        α = max(A[i], α)
    end
    @turbo for i ∈ eachindex(A)
        s += exp(A[i] - α)
    end
    α + log(s)
end

################ threaded version
function taminusb_exp_sum_quote(N::Int, D)
    ls = loopgen(N)
    body = aminusb_exp_sumbody(N, D)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return B
    end
end

@generated function taminusb_exp_sum!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N}, B::AbstractArray{T, N}, dims::D) where {Tₒ, T, N, D}
    taminusb_exp_sum_quote(N, D)
end

function tlogself_plusb!(C::AbstractArray{Tₒ, N}, B::AbstractArray{T, N}) where {Tₒ, T, N}
    @tturbo for i ∈ eachindex(C)
        C[i] = B[i] + log(C[i])
    end
    C
end

function lvtlogsumexp(A::AbstractArray{T, N}, dims::NTuple{M, Int}) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvtlogsumexp1(A))
    else
        B = lvtmaximum(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        C = zeros(Base.promote_op(exp, T), Dᴮ)
        taminusb_exp_sum!(C, A, B, Dᴮ′)
        tlogself_plusb!(C, B)
        C
    end
    return C
end

lvtlogsumexp(A::AbstractArray{T, N}, dims::Int) where {T, N} = lvtlogsumexp(A, (dims,))
lvtlogsumexp(A::AbstractArray{T, N}; dims=:) where {T, N} = lvtlogsumexp(A, dims)
lvtlogsumexp(A::AbstractArray{T, N}) where {T, N} = lvtlogsumexp1(A)
lvtlogsumexp(A::AbstractArray{T, N}, ::Colon) where {T, N} = lvtlogsumexp1(A)

function lvtlogsumexp1(A::AbstractArray{T, N}) where {T, N}
    α = typemin(T)
    s = zero(promote_type(T, Float64))
    @tturbo for i ∈ eachindex(A)
        α = max(A[i], α)
    end
    @tturbo for i ∈ eachindex(A)
        s += exp(A[i] - α)
    end
    α + log(s)
end


text:
"""
Constant Velocity Motion Model
"""
mutable struct MotionModel
    prev_time::Float64
    """
    Element of SE(3) group.
    """
    prev_wc::SMatrix{4, 4, Float64}
    """
    Element of se(3) algebra.
    """
    log_rel_t::SMatrix{4, 4, Float64}
end

function MotionModel(;
    prev_time::Float64 = -1.0,
    prev_wc::SMatrix{4, 4, Float64} = SMatrix{4, 4, Float64}(I),
    log_rel_t::SMatrix{4, 4, Float64} = zeros(SMatrix{4, 4, Float64}),
)
    MotionModel(prev_time, prev_wc, log_rel_t)
end

function reset!(m::MotionModel)
    m.prev_time = -1.0
    m.log_rel_t = zeros(SMatrix{4, 4, Float64})
end

"""
Apply Motion Model to a given `wc` transformation.
"""
function (m::MotionModel)(wc::SMatrix{4, 4, Float64}, time)
    m.prev_time < 0 && return wc
    # `wc` and `m.prev_wc` should be equal here,
    # since prev_wc is updated right after pose computation.
    # If not - update to stay consistent.
    δ = log_lie(SE3, wc * inv(SE3, m.prev_wc))
    all(isapprox.(δ, 0; atol=1e-5)) || (m.prev_wc = wc;)

    δt = time - m.prev_time
    wc * exp_lie(SE3, m.log_rel_t .* δt)
end

function update!(m::MotionModel, wc::SMatrix{4, 4, Float64}, time)
    if m.prev_time < 0
        m.prev_time = time
        m.prev_wc = wc
        return
    end

    δt = time - m.prev_time
    m.prev_time = time
    δt < 0 && error(
        "Got older than previous image! " *
        "Previous time $(m.prev_time) vs time $time.")

    new_rel = inv(SE3, m.prev_wc) * wc
    m.log_rel_t = log_lie(SE3, new_rel) ./ δt
    m.prev_wc = wc
end


text:
abstract type AbstractKeyHandler end
# required get methods: get_name, get_key, get_api_url
# required set methods: set_name!, set_key!, set_api_url!, load_key!

"""
`get_name(kh::AbstractKeyHandler)`

Returns the `name` property of the key handler.
"""
get_name(kh::AbstractKeyHandler) = getproperty(kh,:name)

"""
`get_api_url(kh::AbstractKeyHandler)`

Returns the `api_url` property of the key handler.
"""
get_api_url(kh::AbstractKeyHandler) = getproperty(kh,:api_url)

"""
`get_key(kh::AbstractKeyHandler)`

Returns the `key` property of the key handler.
"""
get_key(kh::AbstractKeyHandler) = getproperty(kh,:key)

"""
`set_name!(kh::AbstractKeyHandler)`

Returns the `name` property of the key handler.
"""
set_name!(kh::AbstractKeyHandler,newname::AbstractString) = setproperty!(kh,:name,newname)

"""
`set_api_url!(kh::AbstractKeyHandler)`

Returns the `api_url` property of the key handler.
"""
set_api_url!(kh::AbstractKeyHandler,newurl::AbstractString) = setproperty!(kh,:api_url,newurl)

"""
`set_key!(kh::AbstractKeyHandler)`

Returns the `key` property of the key handler.
"""
set_key!(kh::AbstractKeyHandler,newkey::AbstractString) = setproperty!(kh,:key,newkey)

function Base.show(io::IO, kh::AbstractKeyHandler)
    @printf io "API Key Handler\n"
    @printf io "\tname: %s\n" get_name(kh)
    @printf io "\turl: %s\n" get_api_url(kh)
    @printf io "\tkey: %s\n" get_key(kh)
end

"""
`load_key!(kh::AbstractKeyHandler; KEY_ENV_NAME="", KEY_FILE_NAME="", verbose=false)`

Attempts to load the `key` from one of two stored locations. It first looks in `ENV[KEY_ENV_NAME]`
and then looks the file `KEY_FILE_NAME`, which should be a path to the file.

`KEY_FILE_NAME` can take the form of any of configurations available in `ConfParser.jl` provided that
 1. For `syntax=keyonly": the key by itself on a single line, with nothing else in the file.
 2. For all others: your key is indexed by the word "key", outside of a named block.

If found, the key is set using `set_key!(kh,key)`.
"""
function load_key!(kh::AbstractKeyHandler;KEY_ENV_NAME::AbstractString="", KEY_FILE_NAME::AbstractString="", verbose::Bool=false)
    key = ""
    if !isempty(KEY_ENV_NAME) & haskey(ENV,KEY_ENV_NAME)
        verbose && println("Key loaded from ENV variable $(KEY_ENV_NAME)")
        key = ENV[KEY_ENV_NAME]
        set_key!(kh,key)
    elseif !isempty(KEY_FILE_NAME) & isfile(KEY_FILE_NAME)
        conf = ConfParse(KEY_FILE_NAME)
        parse_conf!(conf)
        key = retrieve(conf,"key")
        set_key!(kh,key)
        verbose && println("Key loaded from file $(KEY_FILE_NAME)")
    else
        error("Key not detected.")
    end
    isempty(key) && error("Key not detected.")
    nothing
end

"""
`load!(kh::AbstractKeyHandler,filename::AbstractString;verbose::Bool=false)`

Load the key handler data from a configuration file. `filename` can take the form of any of configurations available in `ConfParser.jl` provided that
 1. The information sits in the first block of the file.
 2. The information is indexed by `name`, `api_url`, and `key`.
"""
function load!(kh::AbstractKeyHandler,filename::AbstractString;verbose::Bool=false)
    isempty(filename) && error("Invalid filename: $(filename)")
    !isfile(filename) && error("Invalid filename (does not exist): $(filename)")

    conf = ConfParse(filename)
    parse_conf!(conf)
    verbose && println("Configuration file found with syntax $(conf._syntax).")
    if !isequal(conf._syntax,"ini")
        nm = retrieve(conf,"name")
        url = retrieve(conf,"api_url")
        key = retrieve(conf,"key")
    else
        blockname = first(keys(conf._data))
        nm = retrieve(conf,blockname,"name")
        url = retrieve(conf,blockname,"api_url")
        key = retrieve(conf,blockname,"key")
    end

    set_name!(kh,nm)
    set_api_url!(kh,url)
    set_key!(kh,key)

    nothing
end

"""
`save!(kh::AbstractKeyHandler,filename::AbstractString;verbose::Bool=false)`

Save the key handler data to a configuration file. `filename` can take the form of any of configurations available in `ConfParser.jl`. If `filename` exists, the same configuration is used but the data overwritten. If it does not exist, `defaultsyntax` is used.
"""
function ConfParser.save!(kh::AbstractKeyHandler,filename::AbstractString;defaultsyntax::AbstractString="ini")
    isempty(filename) && error("Invalid filename: $(filename)")
    if isfile(filename)
        conf = ConfParse(filename)
        parse_conf!(conf)
        if !isequal(conf._syntax,"ini")
            commit!(conf,"name",get_name(kh))
            commit!(conf,"api_url",get_api_url(kh))
            commit!(conf,"key",get_key(kh))
        else
            blockname = first(keys(conf._data))
            commit!(conf,blockname,"name",get_name(kh))
            commit!(conf,blockname,"api_url",get_api_url(kh))
            commit!(conf,blockname,"key",get_key(kh))

        end

        save!(conf,filename)
    else
        fh = open(filename,"w")
        if isequal(defaultsyntax,"keyonly")
            print(fh,get_key(kh))
            close(fh)
        else
            sep = ""
            if isequal(defaultsyntax,"ini")
                sep = "="
            elseif isequal(defaultsyntax,"http")
                sep = ":"
            elseif isequal(defaultsyntax,"simple")
                sep = " "
            else
                error("Syntax $(defaultsyntax) not found.")
            end
            println(fh,"name",sep,get_name(kh))
            println(fh,"api_url",sep,get_api_url(kh))
            print(fh,"key",sep,get_key(kh))
            close(fh)
        end
    end
end


text:
# Test whether the accelerator is restarted when rho is adapted (the ADMM operator changed)

using COSMO, SparseArrays, LinearAlgebra, Test

@testset "Rho adaption triggers restarts" begin
    q = [1; 1.];
    P = sparse([4. 1; 1 2]);
    A = [1. 1; 1 0; 0 1];
    l = [1.; 0; 0];
    u = [1; 0.7; 0.7];

    Aa = [-A; A]
    ba = [u; -l]
    constraint1 = COSMO.Constraint(Aa, ba, COSMO.Nonnegatives);

    accelerator = with_options(AndersonAccelerator{Float64, Type2{NormalEquations}, RollingMemory, NoRegularizer}, mem = 5, activate_logging = true)
    settings = COSMO.Settings(adaptive_rho_interval = 23, rho = 1e-4, accelerator = accelerator, safeguard = false);
    model = COSMO.Model();
    assemble!(model, P, q, constraint1, settings = settings);
    res = COSMO.optimize!(model);

    @test length(filter(x -> x[2] == :rho_adapted, model.accelerator.acceleration_status)) == num_rho_adaptions(model)
end


text:
function PrintStats(x::LangOptions)
    @check_ptrs x
    return clang_LangOptions_PrintStats(x)
end


text:
using STCA
using Base.Test

@test STCA.tostr(0x3031_3233)=="0123"
@test STCA.tostr(0x3031_3233_3435_3637)=="01234567"
@test STCA.toval32("0123")==0x3031_3233
@test STCA.toval64("01234567")==0x3031_3233_3435_3637

@test STCA.rotate(0x11223344)==0x44112233
@test STCA.rotate(0x11223344_55667788)==0x44112233_88556677

@test STCA.north(0x11223344)==0x11
@test STCA.east(0x11223344)==0x22
@test STCA.south(0x11223344)==0x33
@test STCA.west(0x11223344)==0x44

d = [
0xaaaaaaaa 0xffff33ff 0xbbbbbbbb;
0xff22ffff 0xeeeeeeee 0xffffff44;
0xcccccccc 0x11ffffff 0xdddddddd
]
sd = view(d,1:3,1:3)
@test STCA.qn(sd)==0x33
@test STCA.qe(sd)==0x44
@test STCA.qs(sd)==0x11
@test STCA.qw(sd)==0x22
@test STCA.get_target_state(sd)==0xeeeeeeee_33441122
STCA.set_target_state(sd,0xffffffff_11223344)
@test STCA.get_target_state(sd)==0xffffffff_11223344

STCA.qn!(sd,0x0000_00ff)
@test d[4]==0xffff_ffff
STCA.qe!(sd,0x0000_00ff)
@test d[8]==0xffff_ffff
STCA.qs!(sd,0x0000_00ff)
@test d[6]==0xffff_ffff
STCA.qw!(sd,0x0000_00ff)
@test d[2]==0xffff_ffff


text:
# Triangle, pentagonal, and hexagonal numbers are generated by the following
# formulae:
#
# Triangle        Tn=n(n+1)/2     1, 3, 6, 10, 15, ...
# Pentagonal      Pn=n(3n−1)/2    1, 5, 12, 22, 35, ...
# Hexagonal       Hn=n(2n−1)      1, 6, 15, 28, 45, ...
#
# It can be verified that T285 = P165 = H143 = 40755.
#
# Find the next triangle number that is also pentagonal and hexagonal.

using ProjectEulerSolutions

# Calculate different pentagonal and hexagonal numbers, since every hexagonal
# number is also a triangular number.  Then find the intersection
function p045solution_set(n::Integer=166)::Integer

    pentvals = collect(165:n)
    hexvals = collect(143:n)
    pentset = Set(fld.(pentvals .* (3 .* pentvals .- 1), 2))
    hexset = Set(hexvals .* (2 .* hexvals .- 1))

    return maximum(intersect(pentset, hexset))
end

# Ignoring triangular numbers, solving for when the pentagonal and hexagonal
# equations are equal gives a "Diophantine equation" of the form
# axx + bxy + cyy + dx + ey + f = 0.  Plugging in, we can get a recursive
# solution from https://www.alpertron.com.ar/QUAD.HTM and iterate through a
# single step.  This idea came from https://projecteuler.net/thread=45
function p045solution_diophantine(n::Integer=0)::Integer
    # x0 = 1
    # y0 = 1
    x1 = 165
    y1 = 143
    x2 = 97 * x1 + 112 * y1 - 44
    y2 = 84 * x1 + 97 * y1 - 38

    return y2 * (2 * y2 - 1) # convert back to hexagonal number from index
end

p045 = Problems.Problem(Dict("Set" => p045solution_set,
                             "Diophantine" => p045solution_diophantine))

Problems.benchmark(p045, 50_000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
using NURBS
using Base.Test

@testset "Bezier Curves" begin

    @testset "Bezier Curve" begin
        b = [
        	1.0 2 4 3;
        	1.0 3 3 1;
        	0.0 0 0 0
        ]
        ans1 = [
    		1.0 1.56481 2.18519 2.75 3.14815 3.26852 3.0;
    		1.0 1.83333 2.33333 2.5  2.33333 1.83333 1.0;
    		0.0 0.0     0.0     0.0  0.0     0.0     0.0
    	]
        ans2 = [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]
        @test typeof(bezier(4, b, 7)[1]) == Array{Float64,2}
        @test typeof(bezier(4, b, 7)[2]) == Array{Array{Int64,1},1}
        @test isapprox(bezier(4, b, 7)[1], ans1, atol=1e-5)
        @test bezier(4, b, 7)[2] == ans2
    end

    @testset "Bezier Curves with derivatives" begin
    	b = [
    		1.0 2 4 3;
    		1.0 3 3 1;
    		0.0 0 0 0
    	]
    	ansP = [
    		1.0 1.56481 2.18519 2.75 3.14815 3.26852 3.0;
    		1.0 1.83333 2.33333 2.5  2.33333 1.83333 1.0;
    		0.0 0.0     0.0     0.0  0.0     0.0     0.0
    	]
    	ansD1 = [
    		3.66667  3.66667  3.0   1.66667  3.26852
			4.0      2.0      0.0  -2.0      1.83333
			0.0      0.0      0.0   0.0      0.0
		]
    	ansD2 = [
    		  2.0   -2.0   -6.0  -10.0  3.26852
			-12.0  -12.0  -12.0  -12.0  1.83333
			  0.0    0.0    0.0    0.0  0.0
    	]
    	ansEV = [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6], [6, 7]]
    	ansDEV = [[1, 2], [2, 3], [3, 4], [4, 5]]

    	@test typeof(dbezier(4, b, 7)[1]) == Array{Float64,2}
    	@test typeof(dbezier(4, b, 7)[2]) == Array{Float64,2}
    	@test typeof(dbezier(4, b, 7)[3]) == Array{Float64,2}
        @test typeof(dbezier(4, b, 7)[4]) == Array{Array{Int64,1},1}
        @test typeof(dbezier(4, b, 7)[5]) == Array{Array{Int64,1},1}

        @test isapprox(dbezier(4, b, 7)[1], ansP, atol=1e-5)
        @test isapprox(dbezier(4, b, 7)[2], ansD1, atol=1e-5)
        @test isapprox(dbezier(4, b, 7)[3], ansD2, atol=1e-5)
        @test dbezier(4, b, 7)[4] == ansEV
        @test dbezier(4, b, 7)[5] == ansDEV
    end

    @testset "Bernstein Basis" begin
        ans1 = [
            [1.0, 0.614, 0.275, 0.125, 0.042, 0.003, 0.0],
            [0.0, 0.325, 0.444, 0.375, 0.239, 0.058, 0.0],
            [0.0, 0.058, 0.239, 0.375, 0.444, 0.325, 0.0],
            [0.0, 0.003, 0.042, 0.125, 0.275, 0.614, 1.0],
        ]
        @test typeof(bern_basis(1, 1, 0.5)) == Float64
        @test isapprox(bern_basis(3, 2, 0.5), 0.375, atol = 1e-3)
        @test isapprox(map(y -> map(x -> bern_basis(3, y, x), [0, 0.15, 0.35, 0.5, 0.65, 0.85, 1]),[0, 1, 2, 3]), ans1, atol=1e-2)
    end
end


text:
module ExprTree
using AbstractTrees, LinearAlgebra, StaticArrays
using ..DiagTree
# using Unrolled
# using InteractiveUtils

using Printf, PyCall

const ADD, MUL = 1, 2
export ADD, MUL

include("common.jl")

# struct Cache and struct Pool
include("pool.jl")

include("tree.jl")
export ExpressionTree

# IO operations
include("io.jl")

# diagram evaluation
include("eval.jl")

include("build.jl")

end


text:
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.0.0 (2018-08-08)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia> isequal(NaN,NaN)
true

julia> NaN == NaN
false

# complex number
julia> im
im

julia> (1 + 2im)*(2 - 3im)
8 + 1im

julia> sqrt(9)
3.0

julia> sqrt(-1)
ERROR: DomainError with -1.0:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
 [1] throw_complex_domainerror(::Symbol, ::Float64) at .\math.jl:31
 [2] sqrt at .\math.jl:479 [inlined]
 [3] sqrt(::Int64) at .\math.jl:505
 [4] top-level scope at none:0

julia> sqrt(im)
0.7071067811865476 + 0.7071067811865475im

julia> real(2 -3im)
2

julia> imag(2 -3im)
-3

julia> a = 1; b = 2; a + b*im
1 + 2im

julia> complex(a,b)
1 + 2im

julia> 6//2
3//1

julia> sin(0)
0.0

julia> cos(0)
1.0

julia> log(0)
-Inf

julia> log(1)
0.0

julia> log(10)
2.302585092994046

julia>


text:
""" Serialize abstract wiring diagrams as JSON.

JSON data formats are convenient when programming for the web. Unfortunately, no
standard for JSON graph formats has gained any kind of widespread adoption. We
adopt a format compatible with that used by the KEILER project and its successor
ELK (Eclipse Layout Kernel). This format is roughly feature compatible with
GraphML, supporting nested graphs and ports. It also supports layout information
like node position and size.

References:

- KEILER's JSON graph format:
  https://rtsys.informatik.uni-kiel.de/confluence/display/KIELER/JSON+Graph+Format
- ELK's JSON graph format:
  https://www.eclipse.org/elk/documentation/tooldevelopers/graphdatastructure/jsonformat.html
"""
module JSONWiringDiagrams
export read_json_graph, parse_json_graph, write_json_graph, generate_json_graph,
  convert_from_json_graph_data, convert_to_json_graph_data

using DataStructures: OrderedDict
import JSON

using ..DirectedWiringDiagrams, ..WiringDiagramSerialization
import ..DirectedWiringDiagrams: PortData

const JSONObject = OrderedDict{String,Any}

# Serialization
###############

""" Write a wiring diagram to a file as JSON.
"""
function write_json_graph(diagram::WiringDiagram, filename::String;
                          indent::Union{Int,Nothing}=nothing)
  open(filename, "w") do io
    JSON.print(io, generate_json_graph(diagram), indent)
  end
end

""" Generate a JSON dict representing a wiring diagram.
"""
function generate_json_graph(diagram::WiringDiagram)::AbstractDict
  generate_json_box(diagram, Int[])
end

function generate_json_box(diagram::WiringDiagram, path::Vector{Int})
  json_object_with_value(diagram.value,
    "id" => box_id(path),
    "ports" => generate_json_ports(diagram),
    "children" => [
      generate_json_box(box(diagram, v), [path; v]) for v in box_ids(diagram)
    ],
    "edges" => [
      json_object_with_value(wire.value,
        "id" => wire_id(path, i),
        "source" => box_id(diagram, [path; wire.source.box]),
        "sourcePort" => port_name(diagram, wire.source),
        "target" => box_id(diagram, [path; wire.target.box]),
        "targetPort" => port_name(diagram, wire.target),
      )
      for (i, wire) in enumerate(wires(diagram))
    ],
  )
end

function generate_json_box(box::Box, path::Vector{Int})
  json_object_with_value(box.value,
    "id" => box_id(path),
    "ports" => generate_json_ports(box),
  )
end

function generate_json_ports(box::AbstractBox)::AbstractArray
  [
    [
      json_object_with_value(port,
        "id" => port_name(InputPort, i),
        "portkind" => "input",
      )
      for (i, port) in enumerate(input_ports(box))
    ];
    [
      json_object_with_value(port,
        "id" => port_name(OutputPort, i),
        "portkind" => "output",
      )
      for (i, port) in enumerate(output_ports(box))
    ];
  ]
end

function json_object_with_value(value, args...)
  obj = JSONObject(args...)
  props = convert_to_json_graph_data(value)
  if !isempty(props)
    obj["properties"] = props
  end
  obj
end

convert_to_json_graph_data(x) = convert_to_graph_data(x)

# Deserialization
#################

""" Read a wiring diagram from a JSON file.
"""
function read_json_graph(
    BoxValue::Type, PortValue::Type, WireValue::Type, filename::String)
  parse_json_graph(BoxValue, PortValue, WireValue, JSON.parsefile(filename))
end

""" Parse a wiring diagram from a JSON string or dict.
"""
function parse_json_graph(
    BoxValue::Type, PortValue::Type, WireValue::Type, s::Union{AbstractString,IO})
  parse_json_graph(BoxValue, PortValue, WireValue, JSON.parse(s))
end
function parse_json_graph(
    ::Type{BoxValue}, ::Type{PortValue}, ::Type{WireValue},
    node::AbstractDict)::WiringDiagram where {BoxValue, PortValue, WireValue}
  diagram, ports = parse_json_box(BoxValue, PortValue, WireValue, node)
  diagram
end

function parse_json_box(
    BoxValue::Type, PortValue::Type, WireValue::Type, node::AbstractDict)
  # Read the ports of the box.
  ports, input_ports, output_ports = parse_json_ports(PortValue, node)

  # Dispense with case when node is an atomic box.
  if !haskey(node, "children")
    value = parse_json_graph_data(BoxValue, node)
    return (Box(value, input_ports, output_ports), ports)
  end

  # If we get here, we're reading a wiring diagram.
  # FIXME: We should not assume that diagram data has same type as box data.
  value = haskey(node, "properties") ? parse_json_graph_data(BoxValue, node) : nothing
  diagram = WiringDiagram(value, input_ports, output_ports)
  all_ports = Dict{Tuple{String,String},Port}()
  for (key, port_data) in ports
    all_ports[key] = port_data.kind == InputPort ?
      Port(input_id(diagram), OutputPort, port_data.port) :
      Port(output_id(diagram), InputPort, port_data.port)
  end

  # Read the nodes recursively.
  for subnode in node["children"]
    box, subports = parse_json_box(BoxValue, PortValue, WireValue, subnode)
    v = add_box!(diagram, box)
    for (key, port_data) in subports
      all_ports[key] = Port(v, port_data.kind, port_data.port)
    end
  end

  # Read the edges.
  for edge in node["edges"]
    value = parse_json_graph_data(WireValue, edge)
    source = all_ports[(edge["source"], edge["sourcePort"])]
    target = all_ports[(edge["target"], edge["targetPort"])]
    add_wire!(diagram, Wire(value, source, target))
  end

  (diagram, ports)
end

function parse_json_ports(PortValue::Type, node::AbstractDict)
  ports = Dict{Tuple{String,String},PortData}()
  input_ports, output_ports = PortValue[], PortValue[]
  node_id = node["id"]
  for port in node["ports"]
    port_id = port["id"]
    port_kind = port["portkind"]
    value = parse_json_graph_data(PortValue, port)
    if port_kind == "input"
      push!(input_ports, value)
      ports[(node_id, port_id)] = PortData(InputPort, length(input_ports))
    elseif port_kind == "output"
      push!(output_ports, value)
      ports[(node_id, port_id)] = PortData(OutputPort, length(output_ports))
    else
      error("Invalid port kind: $portkind")
    end
  end
  (ports, input_ports, output_ports)
end

function parse_json_graph_data(Value::Type, obj::AbstractDict)
  convert_from_json_graph_data(Value, get(obj, "properties", Dict()))
end

convert_from_json_graph_data(T::Type, data) = convert_from_graph_data(T, data)

end


text:
module LatestPkg

export LP

include(joinpath(@__DIR__, "..", "deps", "deps.jl"))
include(pkgpath)
using .Pkg

const LP = Pkg

end # module


text:
using Arpack, Test, TestSetExtensions, LinearOperatorsRevamp

include("test_aux.jl")

include("test_linop.jl")
include("test_linop_allocs.jl")
include("test_adjtrans.jl")
include("test_cat.jl")
include("test_lbfgs.jl")
include("test_lsr1.jl")
include("test_kron.jl")
include("test_callable.jl")
include("test_deprecated.jl")


text:
include("physics.jl")

function centre_of_mass(contraption::Contraption)
    CoM = (0, 0)

    for i in eachindex(contraption.mass)
        CoM += contraption.mass[i]  * contraption.position[:, i]
    end

    CoM /= sum(contraption.mass)

    return CoM
end

function moment_of_inertia(contraption::Contraption)
    CoM = centre_of_mass(contraption)
    I = 0

    for i in eachindex(contraption.mass)
        r = Distances.norm(contraption.position[:, i] - CoM)
        I += contraption.mass[i] * r^2
    end

    return I
end

function momentum(contraption::Contraption)
    momentum = (0, 0)

    for i in eachindex(contraption.mass)
        momentum += contraption.velocity[:, i] * contraption.mass[i]
    end

    return momentum
end

# FIXME: requires 3D coordinates everywhere in order to work
function angular_momentum(contraption::Contraption)
    CoM = centre_of_mass(contraption)
    L = (0, 0, 0)

    for i in eachindex(contraption.mass)
        r = contraption.position[:, i] - CoM
        p = contraption.velocity[:, i] * contraption.mass[i]
        L += Base.LinAlg.cross(r, p)
    end

    return L
end

# Computes average angular velocity derived from angular momentum and moment of inertia
function angular_velocity(contraption::Contraption)
    ω = angular_momentum(contraption) / moment_of_inertia(contraption)

    return ω
end

function kinetic_energy(contraption::Contraption)
    E = 0

    for i in eachindex(contraption.mass)
        E += 0.5 * contraption.mass[i] * contraption.velocity[:, i] ^ 2
    end

    return E
end

function elastic_energy(contraption::Contraption)
    E = 0

    for i in eachindex(contraption.springs)
        if contraption.springs != 0
            Δs = Distances.norm(contraption.position[:, i[1]] - contraption.position[:, i[2]])

            E += 0.5 * contraption.springs[i] * Δs ^ 2
        end
    end

    return E
end

# Gravitational potential energy, computed relative to y = 0
function gravitational_energy(contraption::Contraption, g::Float64)

    E = 0
    for i in eachindex(contraption.mass)
        E += contraption.mass[i] * g * contraption.position[2, i]
    end

    return E
end

function total_energy(contraption::Contraption)

    E = sum(kinetic_energy(contraption),
            elastic_energy(contraption),
            gravitational_energy(contraption))

    return E
end

function detect_collision(velocity_series)

end

function detect_rest(velocity_series)

end

function detect_oscillation(velocity_series)

end


text:
using Documenter, PseudoArcLengthContinuation

makedocs(doctest = false,
	sitename = "Pseudo Arc Length Continuation in Julia",
	format = Documenter.HTML(collapselevel = 1),
	# format = DocumenterLaTeX.LaTeX(),
	authors = "Romain Veltz",
	pages = Any[
		"Home" => "index.md",
		"Overview" => "guidelines.md",
		"Tutorials" => "tutorials.md",
		"Functionalities" => [
			"Plotting" => "plotting.md",
			"Linear / Eigen Solvers" => "linearsolver.md",
			"Bordered linear solvers" => "borderedlinearsolver.md",
			"Bifurcation detection" => "detectionBifurcation.md",
			"Fold / Hopf Continuation (codim 2)" => "codim2Continuation.md",
			"Normal form" =>[
				"Simple branch point" => "simplebp.md",
				"Non-simple branch point" => "nonsimplebp.md",
				"Simple Hopf point" => "simplehopf.md",
			],
			"Branch switching" => "branchswitching.md",
			"Deflated problem" => "deflatedproblem.md",
			"Constrained problem" => "constrainedproblem.md",
			"Periodic Orbits" => [
				"Introduction" => "periodicOrbit.md",
				"Finite Differences" => "periodicOrbitFD.md",
				"Shooting" => "periodicOrbitShooting.md",
				],
			"DiffEq wrapper" => "diffeq.md",
			"Bordered arrays" => "Borderedarrays.md",
			"Iterator Interface" => "iterator.md",
		],
		"Frequently Asked Questions" => "faq.md",
		"Library" => "library.md"
	]
	)

deploydocs(
	repo = "github.com/rveltz/PseudoArcLengthContinuation.jl.git",
)