text:
println("------------------------------------")
println("|     Fields and vector spaces     |")
println("------------------------------------")
ti = time()
@timedtestset "Fields" begin
    @test isa(ℝ, Field)
    @test isa(ℂ, Field)
    @test eval(Meta.parse(sprint(show, ℝ))) == ℝ
    @test eval(Meta.parse(sprint(show, ℂ))) == ℂ
    @test ℝ ⊆ ℝ
    @test ℝ ⊆ ℂ
    @test ℂ ⊆ ℂ
    @test !(ℂ ⊆ ℝ)

    for T in (Int8, Int16, Int32, Int64, BigInt)
        @test one(T) ∈ ℝ
        @test one(Rational{T}) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test !(one(Complex{Rational{T}}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Rational{T}) ∈ ℂ
        @test one(Complex{T}) ∈ ℂ
        @test one(Complex{Rational{T}} ∈ ℂ)

        @test T ⊆ ℝ
        @test Rational{T} ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test !(Complex{Rational{T}} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Rational{T} ⊆ ℂ
        @test Complex{T} ⊆ ℂ
        @test Complex{Rational{T}} ⊆ ℂ
    end
    for T in (Float32, Float64, BigFloat)
        @test one(T) ∈ ℝ
        @test !(one(Complex{T}) ∈ ℝ)
        @test one(T) ∈ ℂ
        @test one(Complex{T} ∈ ℂ)

        @test T ⊆ ℝ
        @test !(Complex{T} ⊆ ℝ)
        @test T ⊆ ℂ
        @test Complex{T} ⊆ ℂ
    end
end

@timedtestset "ElementarySpace: CartesianSpace" begin
    d = 2
    V = ℝ^d
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, CartesianSpace)
    @test !isdual(V)
    @test !isdual(V')
    @test V == CartesianSpace(Trivial() => d) == CartesianSpace(Dict(Trivial() => d))
    @test @constinferred(hash(V)) == hash(deepcopy(V))
    @test V == @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V))
    @test field(V) == ℝ
    @test @constinferred(sectortype(V)) == Trivial
    @test ((@constinferred sectors(V))...,) == (Trivial(),)
    @test length(sectors(V)) == 1
    @test @constinferred(TensorLabXD.hassector(V, Trivial()))
    @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial()))
    @test dim(@constinferred(typeof(V)())) == 0
    @test (sectors(typeof(V)())...,) == ()
    @test @constinferred(TensorLabXD.axes(V)) == Base.OneTo(d)
    @test ℝ^d == ℝ[](d) == CartesianSpace(d) == typeof(V)(d)
    W = @constinferred ℝ^1
    @test @constinferred(oneunit(V)) == W == oneunit(typeof(V))
    @test @constinferred(⊕(V,V)) == ℝ^(2d)
    @test @constinferred(⊕(V,oneunit(V))) == ℝ^(d+1)
    @test @constinferred(⊕(V,V,V,V)) == ℝ^(4d)
    @test @constinferred(fuse(V,V)) == ℝ^(d^2)
    @test @constinferred(fuse(V,V',V,V')) == ℝ^(d^4)
    @test @constinferred(flip(V)) == V'
    @test flip(V) ≅ V
    @test flip(V) ≾ V
    @test flip(V) ≿ V
    @test V ≺ ⊕(V,V)
    @test !(V ≻ ⊕(V,V))
    @test @constinferred(infimum(V, ℝ^3)) == V
    @test @constinferred(supremum(V', ℝ^3)) == ℝ^3
end

@timedtestset "ElementarySpace: ComplexSpace" begin
    d = 2
    V = ℂ^d
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, V'))) == V'
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, ComplexSpace)
    @test !isdual(V)
    @test isdual(V')
    @test V == ComplexSpace(Trivial() => d) == ComplexSpace(Dict(Trivial() => d))
    @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V)) != V
    @test @constinferred(field(V)) == ℂ
    @test @constinferred(sectortype(V)) == Trivial
    @test @constinferred(sectortype(V)) == Trivial
    @test ((@constinferred sectors(V))...,) == (Trivial(),)
    @test length(sectors(V)) == 1
    @test @constinferred(TensorLabXD.hassector(V, Trivial()))
    @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial()))
    @test dim(@constinferred(typeof(V)())) == 0
    @test (sectors(typeof(V)())...,) == ()
    @test @constinferred(TensorLabXD.axes(V)) == Base.OneTo(d)
    @test ℂ^d == Vect[Trivial](d) == Vect[](Trivial()=>d) == ℂ[](d) == typeof(V)(d)
    W = @constinferred ℂ^1
    @test @constinferred(oneunit(V)) == W == oneunit(typeof(V))
    @test @constinferred(⊕(V, V)) == ℂ^(2d)
    @test_throws SpaceMismatch (⊕(V, V'))
    @test_throws MethodError (⊕(ℝ^d, ℂ^d))
    @test_throws MethodError (⊗(ℝ^d, ℂ^d))
    @test @constinferred(⊕(V,V)) == ℂ^(2d)
    @test @constinferred(⊕(V,oneunit(V))) == ℂ^(d+1)
    @test @constinferred(⊕(V,V,V,V)) == ℂ^(4d)
    @test @constinferred(fuse(V,V)) == ℂ^(d^2)
    @test @constinferred(fuse(V,V',V,V')) == ℂ^(d^4)
    @test @constinferred(flip(V)) == V'
    @test flip(V) ≅ V
    @test flip(V) ≾ V
    @test flip(V) ≿ V
    @test V ≺ ⊕(V,V)
    @test !(V ≻ ⊕(V,V))
    @test @constinferred(infimum(V, ℂ^3)) == V
    @test @constinferred(supremum(V', (ℂ^3)')) == dual(ℂ^3) == conj(ℂ^3)
end

@timedtestset "ElementarySpace: GeneralSpace" begin
    d = 2
    V = GeneralSpace{ℂ}(d)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, dual(V)))) == dual(V)
    @test eval(Meta.parse(sprint(show, conj(V)))) == conj(V)
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    @test !isdual(V)
    @test isdual(V')
    @test !isdual(conj(V))
    @test isdual(conj(V'))
    @test !TensorLabXD.isconj(V)
    @test !TensorLabXD.isconj(V')
    @test TensorLabXD.isconj(conj(V))
    @test TensorLabXD.isconj(conj(V'))
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test !isa(V, InnerProductSpace)
    @test !isa(V, EuclideanSpace)
    @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test @constinferred(dual(V)) != @constinferred(conj(V)) != V
    @test @constinferred(field(V)) == ℂ
    @test @constinferred(sectortype(V)) == Trivial
    @test @constinferred(TensorLabXD.hassector(V, Trivial()))
    @test @constinferred(dim(V)) == d == @constinferred(dim(V, Trivial()))
    @test @constinferred(TensorLabXD.axes(V)) == Base.OneTo(d)
end

@timedtestset "ElementarySpace: $(TensorLabXD.type_repr(Vect[I]))" for I in sectorlist
    if Base.IteratorSize(values(I)) === Base.IsInfinite()
        set = unique(vcat(one(I), [randsector(I) for k = 1:10]))
        gen = (c=>2 for c in set)
    else
        gen = (values(I)[k]=>(k+1) for k in 1:length(values(I)))
    end
    V = GradedSpace(gen)
    @test eval(Meta.parse(TensorLabXD.type_repr(typeof(V)))) == typeof(V)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, V'))) == V'
    @test V' == GradedSpace(gen; dual = true)
    @test V == @constinferred GradedSpace(gen...)
    @test V' == @constinferred GradedSpace(gen...; dual = true)
    @test V == @constinferred GradedSpace(tuple(gen...))
    @test V' == @constinferred GradedSpace(tuple(gen...); dual = true)
    @test V == @constinferred GradedSpace(Dict(gen))
    @test V' == @constinferred GradedSpace(Dict(gen); dual = true)
    @test V == @inferred Vect[I](gen)
    @test V' == @constinferred Vect[I](gen; dual = true)
    @test V == @constinferred Vect[I](gen...)
    @test V' == @constinferred Vect[I](gen...; dual = true)
    @test V == @constinferred Vect[I](Dict(gen))
    @test V' == @constinferred Vect[I](Dict(gen); dual = true)
    @test V == @constinferred typeof(V)(c=>dim(V,c) for c in sectors(V))
    if I isa ZNIrrep
        @test V == @constinferred typeof(V)(V.dims)
        @test V' == @constinferred typeof(V)(V.dims; dual = true)
    end
    @test @constinferred(hash(V)) == hash(deepcopy(V)) != hash(V')
    @test V == GradedSpace(reverse(collect(gen))...)
    @test eval(Meta.parse(sprint(show, V))) == V
    @test eval(Meta.parse(sprint(show, typeof(V)))) == typeof(V)
    # space with no sectors
    @test dim(@constinferred(typeof(V)())) == 0
    # space with a single sector
    W = @constinferred GradedSpace(one(I)=>1)
    @test W == GradedSpace(one(I)=>1, randsector(I) => 0)
    @test @constinferred(oneunit(V)) == W == oneunit(typeof(V))
    # randsector never returns trivial sector, so this cannot error
    @test_throws ArgumentError GradedSpace(one(I)=>1, randsector(I) => 0, one(I)=>3)
    @test eval(Meta.parse(sprint(show, W))) == W
    @test isa(V, VectorSpace)
    @test isa(V, ElementarySpace)
    @test isa(V, InnerProductSpace)
    @test isa(V, EuclideanSpace)
    @test isa(V, GradedSpace)
    @test isa(V, GradedSpace{I})
    @test @constinferred(dual(V)) == @constinferred(conj(V)) == @constinferred(adjoint(V)) != V
    @test @constinferred(field(V)) == ℂ
    @test @constinferred(sectortype(V)) == I
    slist = @constinferred sectors(V)
    @test @constinferred(TensorLabXD.hassector(V, first(slist)))
    @test @constinferred(dim(V)) == sum(dim(s)*dim(V, s) for s in slist)
    @constinferred dim(V, first(slist))
    if hasfusiontensor(I)
        @test @constinferred(TensorLabXD.axes(V)) == Base.OneTo(dim(V))
    end
    @test @constinferred(⊕(V,V)) == Vect[I](c=>2dim(V,c) for c in sectors(V))
    @test @constinferred(⊕(V,V,V,V)) == Vect[I](c=>4dim(V,c) for c in sectors(V))
    @test @constinferred(⊕(V,oneunit(V))) ==
            Vect[I](c=>isone(c)+dim(V,c) for c in sectors(V))
    @test @constinferred(fuse(V,oneunit(V))) == V
    d = Dict{I,Int}()
    for a in sectors(V), b in sectors(V)
        for c in a ⊗ b
            d[c] = get(d, c, 0) + dim(V, a)*dim(V, b)*Nsymbol(a,b,c)
        end
    end
    @test @constinferred(fuse(V,V)) == GradedSpace(d)
    @test @constinferred(flip(V)) ==
            Vect[I](conj(c)=>dim(V,c) for c in sectors(V))'
    @test flip(V) ≅ V
    @test flip(V) ≾ V
    @test flip(V) ≿ V
    @test @constinferred(⊕(V, V)) == @constinferred supremum(V, ⊕(V, V))
    @test V == @constinferred infimum(V, ⊕(V, V))
    @test V ≺ ⊕(V,V)
    @test !(V ≻ ⊕(V,V))
    @test infimum(V, GradedSpace(one(I)=>3)) == GradedSpace(one(I)=>2)
    @test_throws SpaceMismatch (⊕(V, V'))
end

@timedtestset "ProductSpace{ℂ}" begin
    V1, V2, V3, V4 = ℂ^1, ℂ^2, ℂ^3, ℂ^4
    P = @constinferred ProductSpace(V1, V2, V3, V4)
    @test eval(Meta.parse(sprint(show, P))) == P
    @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P)
    @test isa(P, VectorSpace)
    @test isa(P, CompositeSpace)
    @test spacetype(P) == ComplexSpace
    @test sectortype(P) == Trivial
    @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P')
    @test P == deepcopy(P)
    @test P == typeof(P)(P...)
    @constinferred (x->tuple(x...))(P)
    @test @constinferred(dual(P)) == P'
    @test @constinferred(field(P)) == ℂ
    @test @constinferred(*(V1, V2, V3, V4)) == P
    @test @constinferred(⊗(V1, V2, V3, V4)) == P
    @test @constinferred(⊗(V1, V2⊗V3⊗V4)) == P
    @test @constinferred(⊗(V1⊗V2, V3⊗V4)) == P
    @test @constinferred(⊗(V1, V2, V3⊗V4)) == P
    @test @constinferred(⊗(V1, V2⊗V3, V4)) == P
    @test @constinferred(insertunit(P, 3)) == V1 * V2 * oneunit(V1) * V3 * V4
    @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3
    @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3
    @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3
    @test fuse(flip(V1), V2, flip(V3)) ≅ V1 ⊗ V2 ⊗ V3
    @test @constinferred(⊗(P)) == P
    @test @constinferred(⊗(V1)) == ProductSpace(V1)
    @test eval(Meta.parse(sprint(show, ⊗(V1)))) == ⊗(V1)
    @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) ==
                @constinferred(one(P)) == @constinferred(one(typeof(P))) ==
                ProductSpace{ComplexSpace}(())
    @test eval(Meta.parse(sprint(show, one(P)))) == one(P)
    @test @constinferred(⊗(one(P), P)) == P
    @test @constinferred(⊗(P, one(P))) == P
    @test @constinferred(⊗(one(P), one(P))) == one(P)
    @test @constinferred(adjoint(P)) == dual(P) == V4' ⊗ V3' ⊗ V2' ⊗ V1'
    @test @constinferred(dims(P)) == map(dim, (V1, V2, V3, V4))
    @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3, V4))
    @test @constinferred(dim(P, 2)) == dim(V2)
    @test @constinferred(sectors(P)) ==
            (mapreduce(sectors, (a, b)->tuple(a..., b...), (V1, V2, V3, V4)),)
    cube(x) = x^3
    @test @constinferred(cube(V1)) == V1 ⊗ V1 ⊗ V1
    N = 3
    @test V1^N == V1 ⊗ V1 ⊗ V1
    @test P^2 == P ⊗ P
    @test @constinferred(dims(P, first(sectors(P)))) == dims(P)
    @test ((@constinferred blocksectors(P))...,) == (Trivial(),)
    @test (blocksectors(P ⊗ ℂ^0)...,) == ()
    @test @constinferred(blockdim(P, first(blocksectors(P)))) == dim(P)
    @test Base.IteratorEltype(P) == Base.IteratorEltype(typeof(P)) ==
                                    Base.IteratorEltype(P.spaces)
    @test Base.IteratorSize(P) == Base.IteratorSize(typeof(P)) ==
                                    Base.IteratorSize(P.spaces)
    @test Base.eltype(P) == Base.eltype(typeof(P)) == typeof(V1)
    @test eltype(collect(P)) == typeof(V1)
    @test collect(P) == [V1, V2, V3, V4]
end

@timedtestset "ProductSpace{SU₂Space}" begin
    V1, V2, V3 = SU₂Space(0=>3, 1//2=>1), SU₂Space(0=>2, 1=>1), SU₂Space(1//2=>1, 1=>1)'
    P = @constinferred ProductSpace(V1, V2, V3)
    @test eval(Meta.parse(sprint(show, P))) == P
    @test eval(Meta.parse(sprint(show, typeof(P)))) == typeof(P)
    @test isa(P, VectorSpace)
    @test isa(P, CompositeSpace)
    @test spacetype(P) == SU₂Space
    @test sectortype(P) == Irrep[SU₂] == SU2Irrep
    @test @constinferred(hash(P)) == hash(deepcopy(P)) != hash(P')
    @test @constinferred(dual(P)) == P'
    @test @constinferred(field(P)) == ℂ
    @test @constinferred(*(V1, V2, V3)) == P
    @test @constinferred(⊗(V1, V2, V3)) == P
    @test @constinferred(adjoint(P)) == dual(P) == V3' ⊗ V2' ⊗ V1'
    @test @constinferred(insertunit(P, 3; conj = true)) == V1 * V2 * oneunit(V1)' * V3
    @test fuse(V1, V2', V3) ≅ V1 ⊗ V2' ⊗ V3
    @test fuse(V1, V2', V3) ≾ V1 ⊗ V2' ⊗ V3 ≾ fuse(V1 ⊗ V2' ⊗ V3)
    @test fuse(V1, V2') ⊗ V3 ≾ V1 ⊗ V2' ⊗ V3
    @test fuse(V1, V2', V3) ≿ V1 ⊗ V2' ⊗ V3 ≿ fuse(V1 ⊗ V2' ⊗ V3)
    @test V1 ⊗ fuse(V2', V3) ≿ V1 ⊗ V2' ⊗ V3
    @test fuse(flip(V1) ⊗ V2) ⊗ flip(V3) ≅ V1 ⊗ V2 ⊗ V3
    @test @constinferred(⊗(V1)) == ProductSpace(V1)
    @test @constinferred(one(V1)) == @constinferred(one(typeof(V1))) ==
                @constinferred(one(P)) == @constinferred(one(typeof(P))) ==
                ProductSpace{ComplexSpace}(())
    @test @constinferred(dims(P)) == map(dim, (V1, V2, V3))
    @test @constinferred(dim(P)) == prod(dim, (V1, V2, V3))
    for s in @constinferred(sectors(P))
        @test hassector(P, s)
        @test @constinferred(dims(P, s)) == dim.((V1, V2, V3), s)
    end
    @test sum(dim(c)*blockdim(P, c) for c in @constinferred(blocksectors(P))) == dim(P)
end

@timedtestset "Deligne tensor product of spaces" begin
    V1 = SU₂Space(0=>3, 1//2=>1)
    V2 = SU₂Space(0=>2, 1=>1)'
    V3 = ℤ₃Space(0=>3, 1=>2, 2=>1)
    V4 = ℂ^3

    for W1 in (V1, V2, V3, V4)
        for W2 in (V1, V2, V3, V4)
            for W3 in (V1, V2, V3, V4)
                for W4 in (V1, V2, V3, V4)
                    Ws = @constinferred(W1 ⊠ W2 ⊠ W3 ⊠ W4)
                    @test Ws == @constinferred((W1 ⊠ W2) ⊠ (W3 ⊠ W4)) ==
                                @constinferred(((W1 ⊠ W2) ⊠ W3) ⊠ W4) ==
                                @constinferred((W1 ⊠ (W2 ⊠ W3)) ⊠ W4) ==
                                @constinferred(W1 ⊠ ((W2 ⊠ W3)) ⊠ W4) ==
                                @constinferred(W1 ⊠ (W2 ⊠ (W3 ⊠ W4)))
                    I1, I2, I3, I4 = map(sectortype, (W1, W2, W3, W4))
                    I = sectortype(Ws)
                    @test I == @constinferred((I1 ⊠ I2) ⊠ (I3 ⊠ I4)) ==
                                @constinferred(((I1 ⊠ I2) ⊠ I3) ⊠ I4) ==
                                @constinferred((I1 ⊠ (I2 ⊠ I3)) ⊠ I4) ==
                                @constinferred(I1 ⊠ ((I2 ⊠ I3)) ⊠ I4) ==
                                @constinferred(I1 ⊠ (I2 ⊠ (I3 ⊠ I4)))
                    @test dim(Ws) == dim(W1) * dim(W2) * dim(W3) * dim(W4)
                end
            end
        end
    end
    @test sectortype(@constinferred((V1 ⊗ V2) ⊠ V3)) == @constinferred(Irrep[SU₂ × ℤ₃])
    @test dim((V1 ⊗ V2) ⊠ V3) == dim(V1) * dim(V2) * dim(V3)
    @test sectortype((V1 ⊗ V2) ⊠ V3 ⊠ V4) == Irrep[SU₂ × ℤ₃]
    @test dim((V1 ⊗ V2) ⊠ V3 ⊠ V4) == dim(V1) * dim(V2) * dim(V3) * dim(V4)
    @test fuse(V2 ⊠ V4) == fuse(V4 ⊠ V2) == SU₂Space(0=>6, 1=>3)
    @test fuse(V3 ⊠ V4) == fuse(V4 ⊠ V3) == ℤ₃Space(0=>9, 1=>6, 2=>3)
end

@timedtestset "HomSpace" begin
    V1, V2, V3, V4, V5 = SU₂Space(0=>3, 1//2=>1), SU₂Space(0=>2, 1=>1),
                            SU₂Space(1//2=>1, 1=>1)', SU₂Space(0=>2, 1//2=>2),
                            SU₂Space(0=>1, 1//2=>1, 3//2=>1)'
    W = TensorLabXD.HomSpace(V1 ⊗ V2, V3 ⊗ V4 ⊗ V5)
    @test W == (V3 ⊗ V4 ⊗ V5 → V1 ⊗ V2)
    @test W == (V1 ⊗ V2 ← V3 ⊗ V4 ⊗ V5)
    @test W' == (V1 ⊗ V2 → V3 ⊗ V4 ⊗ V5)
    @test eval(Meta.parse(sprint(show, W))) == W
    @test eval(Meta.parse(sprint(show, typeof(W)))) == typeof(W)
    @test spacetype(W) == SU₂Space
    @test sectortype(W) == Irrep[SU₂]
    @test W[1] == V1
    @test W[2] == V2
    @test W[3] == V3'
    @test W[4] == V4'
    @test W[5] == V5'
    @test @constinferred(hash(W)) == hash(deepcopy(W)) != hash(W')
    @test W == deepcopy(W)
end

tf = time()
printstyled("Finished vector space tests in ",
            string(round(tf-ti; sigdigits=3)),
            " seconds."; bold = true, color = Base.info_color())
println()


text:
module Day23

using AdventOfCode2019

struct Computer
    task::Union{Task,Nothing}
    input::Union{Channel{Int},Nothing}
    output::Union{Channel{Int},Nothing}
    waiting::Union{Channel{Bool},Nothing}
    address::Int
end

function all_waiting(network::Dict{Int,Computer})
    for (k, v) in network
        k == 255 && continue
        if !isready(v.waiting) || isready(v.input)
            return false
        end
    end
    return true
end

function has_unprocessed_output(network::Dict{Int,Computer})
    for (k, v) in network
        k == 255 && continue
        if isready(v.output)
            return true
        end
    end
    return false
end

function initialize_network(program::Array{Int,1})
    network = Dict{Int,Computer}()
    for i = 0:49
        inp = Channel{Int}(100)
        out = Channel{Int}(100)
        waiting = Channel{Bool}(1)
        address = i
        t = @async AdventOfCode2019.IntCode.run_program!(copy(program), inp, out, nothing; waitingForInput = waiting)
        network[i] = Computer(t, inp, out, waiting, address)

        # Initialize computer by assigning its network address
        put!(network[i].input, i)
        put!(network[i].input, -1)
    end
    network[255] = Computer(nothing, Channel{Int}(100), Channel{Int}(100), Channel{Bool}(1), 255)
    return network
end

function shutdown(network::Dict{Int,Computer})
    for (k, v) in network
        close(v.input)
        close(v.output)
        close(v.waiting)
    end
end

function day23(input::String = readInput(joinpath(@__DIR__, "..", "data", "day23.txt")))
    program = parse.(Int, split(input, ','))
    return solve(program)
end

function solve(program::Array{Int,1})
    network = initialize_network(program)
    firstYChannel = Channel{Int}(1)
    lastYChannel = Channel{Int}(1)
    currentNATPacket = [-1, -1]
    @label resume
    while true
        read = false
        for (k, v) in network
            k == 255 && continue
            if isready(v.output)
                read = true
                a = take!(v.output)
                x = take!(v.output)
                y = take!(v.output)
                put!(network[a].input, x)
                put!(network[a].input, y)
            end
        end
        if !read
            for (k, v) in network
                k == 255 && continue
                put!(v.input, -1)
            end
        end
        while true
            if has_unprocessed_output(network)
                break
            elseif all_waiting(network)
                while isready(network[255].input)
                    currentNATPacket[1] = take!(network[255].input)
                    currentNATPacket[2] = take!(network[255].input)
                    if !isready(firstYChannel)
                        put!(firstYChannel, currentNATPacket[2])
                    end
                end
                put!(network[0].input, currentNATPacket[1])
                put!(network[0].input, currentNATPacket[2])
                if isready(lastYChannel)
                    if take!(lastYChannel) == currentNATPacket[2]
                        part1 = take!(firstYChannel)
                        close(firstYChannel)
                        close(lastYChannel)
                        shutdown(network)
                        return [part1, currentNATPacket[2]]
                    end
                end
                put!(lastYChannel, currentNATPacket[2])
                break
            end
            yield()
        end
    end
end

end # module

ext: jl


text:
### helpful utils
# Used for removing LineNumberNodes from generated expressions for equality comparison
# purposes.
using MacroTools: MacroTools

@testset "Utils" begin
    TiledArrays = OneDNN.TiledArrays
    ### cdiv
    @test TiledArrays.cdiv(Int32(10), Int64(9)) == 2
    @test TiledArrays.cdiv(10, 10) == 1
    @test TiledArrays.cdiv(10, 11) == 1

    ### Test the "head" utility function.
    @test TiledArrays.head(1) == ()
    @test TiledArrays.head(1, 2) == (1,)
    @test TiledArrays.head(1, 2, 3) == (1, 2)
    @test TiledArrays.head(1, 2, 3, 4) == (1, 2, 3)
    @test TiledArrays.head(1, 2, 3, 4, 5) == (1, 2, 3, 4)
    @test TiledArrays.head(1, 2, 3, 4, 5, 6) == (1, 2, 3, 4, 5)

    ### strip_leading_tiles
    Tile = OneDNN.TiledArrays.Tile
    strip_leading_tiles = OneDNN.TiledArrays.strip_leading_tiles
    x = (1, 2, 3)
    @test strip_leading_tiles(x) == x

    x = (Tile(1, 2), 3, 4, 5)
    @test strip_leading_tiles(x) == (3, 4, 5)

    x = (Tile(1, 2), Tile(3, 4), 5, 6)
    @test strip_leading_tiles(x) == (5, 6)

    x = (Tile(1, 2), Tile(3, 4), 5, Tile(6, 7), 8)
    @test strip_leading_tiles(x) == (5, Tile(6, 7), 8)

    ### adjust_for_padding
    adjust_for_padding = OneDNN.TiledArrays.adjust_for_padding
    size = (13, 17, 512)
    padded_size = (16, 32, 512)
    @test adjust_for_padding(size, padded_size, (0, 0, 0)) == (0, 0, 0)
    @test adjust_for_padding(size, padded_size, (1, 1, 1)) == (1, 1, 1)
    @test adjust_for_padding(size, padded_size, (12, 17, 10)) == (12, 32, 10)
    @test adjust_for_padding(size, padded_size, (13, 17, 10)) == (16, 32, 10)
    @test adjust_for_padding(size, padded_size, (14, 16, 512)) == (16 + 1, 16, 512)
    @test adjust_for_padding(size, padded_size, (27, 18, 513)) == (2 * 16 + 1, 32 + 1, 513)
end

@testset "Validation" begin
    TiledArrays = OneDNN.TiledArrays

    # First, a normal layout should work fine.
    validate = TiledArrays.validate
    Tile = TiledArrays.Tile
    @test validate((1,))
    @test validate((1, 2))
    @test validate((1, 2, 3))
    @test validate((1, 2, 3, 4))

    # Permutations should be fine
    @test validate((2, 1))
    @test validate((2, 3, 1))
    @test validate((4, 2, 3, 1))

    # Negative indices should throw an error.
    @test_throws Exception validate((-1,))
    @test_throws Exception validate((1, 2, -3))

    # Skipping an index should be an error
    @test_throws Exception validate((2,))
    @test_throws Exception validate((1, 3))

    # Correct usage of tiling.
    @test validate((Tile(3, 8), 2, 1, 3))
    # Nested Tiling
    @test validate((Tile(3, 8), 2, Tile(3, 2), 1, 3))
    # Multiple Tiling
    @test validate((Tile(3, 8), Tile(1, 8), 2, 1, 3))
    # Multiple Tiling with Nesting
    @test validate((Tile(3, 8), Tile(1, 8), 2, Tile(3, 8), 1, 3))

    # Incorrect usage.
    # Tile after int.
    @test_throws Exception validate((Tile(3, 8), 2, 1, 3, 4, Tile(3, 2)))
    # Multiple Ints
    @test_throws Exception validate((1, 1))
    @test_throws Exception validate((1, Tile(2, 2), 1, 2))
    # Tile with no int
    @test_throws Exception validate((Tile(3, 8), 1, 2))
    # Tile with negative index
    @test_throws Exception validate((Tile(-2, 8), 1, 2))
end

@testset "Testing Splitting" begin
    TiledArrays = OneDNN.TiledArrays
    Tile = TiledArrays.Tile
    _splitindex(vlayout, I) = TiledArrays.splitindex(vlayout, I .- one(Int))

    # Now that'd we've tested the code generation logic ... lets test the ACTUAL logic.
    vlayout = Val((2, 1, 3))
    @test _splitindex(vlayout, (10, 20, 30)) == (20, 10, 30) .- 1

    vlayout = Val((Tile(3, 8), 2, 1, 3))
    @test _splitindex(vlayout, (1, 1, 1)) == (0, 0, 0, 0)
    @test _splitindex(vlayout, (1, 1, 2)) == (1, 0, 0, 0)
    @test _splitindex(vlayout, (1, 1, 8)) == (7, 0, 0, 0)
    # Hit the tile barrier
    @test _splitindex(vlayout, (1, 1, 9)) == (0, 0, 0, 1)
    @test _splitindex(vlayout, (1, 1, 16)) == (7, 0, 0, 1)
    @test _splitindex(vlayout, (1, 1, 17)) == (0, 0, 0, 2)

    @test _splitindex(vlayout, (2, 3, 8)) == (7, 3 - 1, 2 - 1, 0)
    @test _splitindex(vlayout, (2, 3, 9)) == (0, 3 - 1, 2 - 1, 1)

    # Multiple levels of tiling
    vlayout = Val((Tile(3, 8), 2, Tile(3, 2), 1, 3))
    @test _splitindex(vlayout, (1, 1, 1)) == (0, 0, 0, 0, 0)
    @test _splitindex(vlayout, (1, 1, 2)) == (1, 0, 0, 0, 0)
    @test _splitindex(vlayout, (1, 1, 8)) == (7, 0, 0, 0, 0)
    @test _splitindex(vlayout, (1, 1, 9)) == (0, 0, 1, 0, 0)
    @test _splitindex(vlayout, (1, 1, 16)) == (7, 0, 1, 0, 0)
    @test _splitindex(vlayout, (1, 1, 17)) == (0, 0, 0, 0, 1)
    @test _splitindex(vlayout, (1, 1, 24)) == (7, 0, 0, 0, 1)
    @test _splitindex(vlayout, (1, 1, 25)) == (0, 0, 1, 0, 1)
end

@testset "Testing Layout Conversion" begin
    # The goal here is to test that we're doing the correct thing when converting from
    # OneDNN layouts for ALL of the layout types in OneDNN.
    #
    # So, we need to:
    #
    # 1. Generate a candidate set of arrays from 1 to 12 dimensions.
    #    Note that we need to have kind of weird dimensionality on these arrays to make
    #    sure we deal with OneDNN's padding correctly.
    #
    # 2. Iterate through all the unique format tags available, converting the arrays of
    #    the correct dimension to that format.
    #
    # 3. Use our handy-dandy format conversion tool to pull out the appropriate indices,
    #    collect the Memory's parent using a view with these indices, and check if a
    #    reshaped version of the view matches the original array.

    # Make the candidate set of arrays.
    # Construct 5 arrays of each dimension.
    arrays_per_dim = 5
    cdiv(a, b) = ceil(Int, a / b)

    # Characters for determining the number of dimensions.
    # Just go up to 4 dims to keep this somewhat reasonable.
    # Keep the length of `chars` to 5 so we can detect if there are 6 or more dimensions
    # and bail.
    chars = ['a', 'b', 'c', 'd', 'e']
    # Drop the first 5 characters to avoid the leading "dnnl_".
    getndims(x::Symbol) = something(findlast(in(lowercase(string(x))[5:end]), chars))

    function rand_in_range(lower, upper, ndims)
        return randn(
            Float32, (rand(cdiv(lower, 4 * j):cdiv(upper, 4 * j)) for j in 1:ndims)...
        )
    end

    dimension_to_array = Dict(
        map(1:length(chars)) do i
            # Cut down upper bound to try to conserve memory consumption.
            return i => [rand_in_range(13, 512, i) for _ in 1:arrays_per_dim]
        end,
    )
    @show Base.summarysize(dimension_to_array)

    # Create initial "Memory" objects so we don't always have to be recreating these.
    dimension_to_memory = Dict(k => OneDNN.Memory.(v) for (k, v) in dimension_to_array)

    # Formats to ignore
    ignore = [:dnnl_format_tag_undef, :dnnl_format_tag_any]
    exit_on = :dnnl_format_tag_last
    seen = Set{UInt32}()

    # Formats that have failed the round-trip test.
    broken_formats = Set([:dnnl_ABc4a4b, :dnnl_ABcd4a4b])

    failed_formats = Set{Symbol}()
    translation_failed = Set{Symbol}()
    meter = ProgressMeter.Progress(Int(OneDNN.Lib.dnnl_format_tag_last), 1)
    for (name, value) in CEnum.name_value_pairs(OneDNN.Lib.dnnl_format_tag_t)
        ProgressMeter.next!(meter)
        # Bookkeeping ...
        in(name, ignore) && continue
        name == exit_on && break
        in(value, seen) && error("Already seen: $name => $(Int(value))")
        push!(seen, value)

        ndims = getndims(name)
        ndims === lastindex(chars) && continue

        @test in(ndims, 1:length(chars))

        arrays = dimension_to_array[ndims]
        memories = dimension_to_memory[ndims]

        for (array, memory) in zip(arrays, memories)
            new_md = OneDNN.memorydesc(
                Float32, size(memory), OneDNN.Lib.dnnl_format_tag_t(value)
            )
            reordered = OneDNN.reorder(new_md, memory)

            # Can we at least round-trip `materialize`.
            passed = (array == OneDNN.materialize(reordered))

            # Some formats were broken when these tests were written.
            # Here, we put in a test to see if they remain broken in the future.
            if !passed
                push!(failed_formats, name)
                continue
            end

            # Now, we assume that everything is working correctly.
            @test passed

            indices = OneDNN.generate_linear_indices(reordered)
            formatted = reshape(view(parent(reordered), indices), size(array))
            translation_passed = (array == formatted)
            @test translation_passed
            if !translation_passed
                push!(translation_failed, name)
            end
        end
    end
    ProgressMeter.finish!(meter)
    @test failed_formats == broken_formats
    if !isempty(translation_failed)
        @show translation_failed
    end
end


text:
module NTupleKmers

export
    Kmer,
    DNAKmer,
    DNAKmer27,
    DNAKmer31,
    DNAKmer63,
    kmertype,
    capacity,
    n_unused

using BioSequences

import BioSequences: twobitnucs, BioSequence, Alphabet, decode

###
### Type definition
###

# In BioSequences.jl kmers are called Mer, not Kmer so there is no name clash.
"""
    Kmer{A<:NucleicAcidAlphabet{2},K,N} <: BioSequence{A}

A parametric, immutable, bitstype for representing Kmers - short sequences.

Given the number of Kmers generated from raw sequencing reads, avoiding
repetetive memory allocation and triggering of garbage collection is important,
as is the ability to effectively pack Kmers into arrays and similar collections.

In julia this means an immutable bitstype must represent such shorter Kmer
sequences. Thankfully this is not much of a limitation - kmers are rarely
manipulated and so by and large don't have to be mutable like `LongSequence`s.

Excepting their immutability, they fulfill the rest of the API and behaviours
expected from a concrete `BioSequence` type.

!!! warning
    Given their immutability, `setindex` and mutating sequence transformations
    are not implemented for kmers e.g. `reverse_complement!`.

!!! tip
    Note that some sequence transformations that are not mutating are
    available, since they can return a new kmer value as a result e.g.
    `reverse_complement`.
"""
struct Kmer{A<:NucleicAcidAlphabet{2},K,N}
    data::NTuple{N,UInt64}
end

###
### Shortcuts
###

"Shortcut for the type `Kmer{DNAAlphabet{2},K,N}`"
const DNAKmer{K,N} = Kmer{DNAAlphabet{2},K,N}

"Shortcut for the type `DNAKmer{27,1}`"
const DNAKmer27 = DNAKmer{27,1}

"Shortcut for the type `DNAKmer{31,1}`"
const DNAKmer31 = DNAKmer{31,1}

"Shortcut for the type `DNAKmer{63,2}`"
const DNAKmer63 = DNAKmer{63,2}

# TODO: Delete when Kmer is made to <: BioSequence.
Alphabet(kmer::Kmer{A,K,N}) where {A,K,N} = A()

###
### Base Functions
###

"""
    kmertype(::Type{Kmer{A,K}}) where {A,K}

Fully resolve and incomplete kmer typing, computing the N parameter of
`Kmer{A,K,N}`, given only `Kmer{A,K}`.

## Example

```julia
julia> DNAKmer{63}
Kmer{DNAAlphabet{2},63,N} where N

julia> kmertype(DNAKmer{63})
Kmer{DNAAlphabet{2},63,2}

```
"""
@inline function kmertype(::Type{Kmer{A,K}}) where {A,K}
    return Kmer{A,K,ifelse(rem(2K, 64) != 0, div(2K, 64) + 1, div(2K, 64))}
end

@inline capacity(::Type{Kmer{A,K,N}}) where {A,K,N} = div(64N, 2)
@inline n_unused(::Type{Kmer{A,K,N}}) where {A,K,N} = capacity(Kmer{A,K,N}) - K
@inline n_unused(seq::Kmer) = n_unused(typeof(seq))

@inline function checkmer(::Type{Kmer{A,K,N}}) where {A,K,N}
    c = capacity(Kmer{A,K,N})
    if !(1 ≤ K ≤ c)
        throw(ArgumentError("K must be within 1..$c"))
    end
end

@inline Base.eltype(::Type{Kmer{A,K,N}}) where {A,K,N} = eltype(A)
@inline Base.length(x::Kmer{A,K,N}) where {A,K,N} = K
@inline Base.summary(x::Kmer{DNAAlphabet{2},K,N}) where {K,N} = string("DNA ", K, "-mer")
@inline Base.summary(x::Kmer{RNAAlphabet{2},K,N}) where {K,N} = string("RNA ", K, "-mer")

function Base.typemin(::Type{Kmer{A,K,N}}) where {A,K,N}
    checkmer(T)
    return T(ntuple(i -> zero(UInt64), N))
end

function Base.typemax(::Type{Kmer{A,K,N}}) where {A,K,N}
    checkmer(Kmer{A,K,N})
    return Kmer{A,K,N}((typemax(UInt64) >> (64N - 2K), ntuple(i -> typemax(UInt64),N - 1)...))
end

function Base.rand(::Type{Kmer{A,K,N}}) where {A,K,N}
    return Kmer{A,K,N}(ntuple(i -> rand(UInt64), N))
end

###
### Constructors
###

# Create a Mer from a sequence.
function (::Type{Kmer{A,K,N}})(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet{2},K,N}
    seqlen = length(seq)
    if seqlen != K
        throw(ArgumentError("seq does not contain the correct number of nucleotides ($seqlen ≠ $K)"))
    end
    if seqlen > capacity(Kmer{A,K,N})
        throw(ArgumentError("Cannot build a mer longer than $(capacity(Kmer{A,K,N}))bp long"))
    end

    # Construct the head.
    bases_in_head = div(64 - (64N - 2K), 2)
    head = zero(UInt64)
    @inbounds for i in 1:bases_in_head
        nt = convert(eltype(typeof(seq)), seq[i])
        head = (head << 2) | UInt64(twobitnucs[reinterpret(UInt8, nt) + 0x01])
    end

    # And the rest of the sequence
    idx = Ref(bases_in_head + 1)

    tail = ntuple(Val{N - 1}()) do i
        Base.@_inline_meta
        body = zero(UInt64)
        @inbounds for i in 1:32
            nt = convert(eltype(typeof(seq)), seq[idx[]])
            body = (body << 2) | UInt64(twobitnucs[reinterpret(UInt8, nt) + 0x01])
            idx[] += 1
        end
        return body
    end

    return Kmer{A,K,N}((head, tail...))
end

function (::Type{BigMer{A,K}})(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet{2},K}
    seqlen = length(seq)
    if seqlen != K
        throw(ArgumentError("seq does not contain the correct number of nucleotides ($seqlen ≠ $K)"))
    end
    if seqlen > BioSequences.capacity(BigMer{A,K})
        throw(ArgumentError("Cannot build a mer longer than $(BioSequences.capacity(BigMer{A,K}))bp long"))
    end

    x = zero(BioSequences.encoded_data_type(BigMer{A,K}))
    for c in seq
        nt = convert(eltype(BigMer{A,K}), c)
        x = (x << 2) | BioSequences.encoded_data_type(BigMer{A,K})(twobitnucs[reinterpret(UInt8, nt) + 0x01])
    end

    return BigMer{A,K}(x)
end

include("bitindex.jl")
include("indexing.jl")
include("predicates.jl")
include("transformations.jl")

@inline function choptail(x::NTuple{N,UInt64}) where {N}
    ntuple(Val{N - 1}()) do i
        Base.@_inline_meta
        return @inbounds x[i]
    end
end

@inline function setlast(x::NTuple{N,UInt64}, nt::DNA) where {N}
    @inbounds begin
        bits = UInt64(twobitnucs[reinterpret(UInt8, nt) + 0x01])
        tail = (x[N] & (typemax(UInt64) - UInt64(3))) | bits
    end
    return (choptail(x)..., tail)
end

"""
    shiftright

It is important to be able to efficiently shift the all the nucleotides in a kmer
one space to the left or right, as it is a key operation in iterating through
de bruijn graph neighbours or in building kmers a nucleotide at a time.
"""
function shiftright end

@inline shiftright(x::BigDNAMer{K}) where {K} = BigDNAMer{K}(reinterpret(UInt128, x) >> 2)

@inline function shiftright(x::Kmer{A,K,N}) where {A,K,N}
    return Kmer{A,K,N}(_shiftright(zero(UInt64), x.data...))
end

@inline function _shiftright(carry::UInt64, head::UInt64, tail...)
    return ((head >> 2) | carry, _shiftright((head & UInt64(3)) << 62, tail...)...)
end
@inline _shiftright(carry::UInt64) = ()

"""
    shiftleft

It is important to be able to efficiently shift the all the nucleotides in a kmer
one space to the left or right, as it is a key operation in iterating through
de bruijn graph neighbours or in building kmers a nucleotide at a time.
"""
function shiftleft end

@inline shiftleft(x::BigDNAMer{K}) where {K} = BigDNAMer{K}(reinterpret(UInt128, x) << 2)

@inline function shiftleft(x::Kmer{A,K,N}) where {A,K,N}
    _, newbits = _shiftleft(x.data...)
    # TODO: The line below is a workaround for julia issues #29114 and #3608
    newbits′ = newbits isa UInt64 ? (newbits,) : newbits
    return Kmer{A,K,N}(_cliphead(64N - 2K, newbits′...))
end

@inline function shiftleft(x::NTuple{N,UInt64}) where {N}
    _, newbits = _shiftleft(x...)
    # TODO: The line below is a workaround for julia issues #29114 and #3608
    return newbits isa UInt64 ? (newbits,) : newbits
end

@inline function _cliphead(by::Integer, head::UInt64, tail...)
    return (head & (typemax(UInt64) >> by), tail...)
end

@inline function _shiftleft(head::UInt64, tail...)
    carry, newtail = _shiftleft(tail...)
    # TODO: The line below is a workaround for julia issues #29114 and #36087
    newtail′ = newtail isa UInt64 ? (newtail,) : newtail
    return head >> 62, ((head << 2) | carry, newtail′...)
end

@inline _shiftleft(head::UInt64) = (head & 0xC000000000000000) >> 62, head << 2

#=
@inline function shiftleft2(x::Kmer{A,K,N}) where {A,K,N}
    return _cliphead(64N - 2K, _shiftleft2(x.data...)...)
end

@inline function _shiftleft2(head::UInt64, next::UInt64, tail...)
    return (head << 2 | (next >> 62), _shiftleft2(next, tail...)...)
end
@inline _shiftleft2(head::UInt64) = head << 2
=#

end # module


text:
using TraitWrappers
using Test

using Base:IteratorEltype, IteratorSize

@testset "IterableTraitWrapper" begin
    a = IterableTraitWrapper(1:8)
    @test IteratorEltype(a) == IteratorEltype(1:8)
    @test IteratorSize(a) == IteratorSize(1:8)
    @test iterate(a) == iterate(1:8)
    @test iterate(a, 2) == iterate(1:8, 2)
    @test eltype(a) == eltype(1:8)
    @test length(a) == length(1:8)
    @test size(a) == size(1:8)
end


text:
# Main parser function
# Constructs an AST from the code
function parse(src)
    toks = lparse(src, regs)
    ast = pparse(toks)

    return ast
end


text:
using Pkg
Pkg.build()
Pkg.test("CoverageCore"; coverage=true)


text:
function rayleigh(v, x, len, E0, E, P)
   w = v / sqrt(dot(v, P, v))
   return 2.0 * (E(x+len/2*w) - 2.0 * E0 + E(x-len/2*v)) / len^2
end

localmerit(x, x0, v0, len, g0, λ0, E) = (
   0.5 * ( E(x+len/2*v0) + E(x-len/2*v0) )
   - 2.0 * dot(v0, g0) * dot(v0, x-x0)
   - λ0 * dot(v0, x-x0)^2  )

function run!{T}(method::StaticDimer, E, dE, x0::Vector{T}, v0::Vector{T})

   # read all the parameters
   @unpack a_trans, a_rot, tol_trans, tol_rot, maxnumdE, len,
            precon_prep!, verbose, precon_rot, rescale_v = method
   P0=method.precon
   # initialise variables
   x, v = copy(x0), copy(v0)
   nit = 0
   numdE, numE = 0, 0
   log = DimerLog()
   # and just start looping
   if verbose >= 2
      @printf(" nit |  |∇E|_∞    |∇R|_∞     λ   \n")
      @printf("-----|-----------------------------\n")
   end
   for nit = 0:maxnumdE  # there will be at most this many evaluations
      # normalise v
      P0 = precon_prep!(P0, x)
      v /= sqrt(dot(v, P0, v))
      # evaluate gradients, and more stuff
      dE0 = dE(x)
      dEv = dE(x + len * v)
      numdE += 2
      Hv = (dEv - dE0) / len

      # NEWTON TYPE RESCALING IN v DIRECTION
      if rescale_v
         P = PreconSMW(P0, v, abs(dot(Hv, v)) - 1.0)
         v /= sqrt(dot(v, P, v))
      else
         P = P0
      end

      # translation and rotation residual, store history
      res_trans = vecnorm(dE0, Inf)
      q_rot = - Hv + dot(v, Hv) * (P * v)
      res_rot = vecnorm(q_rot, Inf)
      push!(log, numE, numdE, res_trans, res_rot)
      if verbose >= 2
         @printf("%4d | %1.2e  %1.2e  %4.2f  \n", nit, res_trans, res_rot, dot(v, Hv))
      end
      if res_trans <= tol_trans && res_rot <= tol_rot
         if verbose >= 1
            println("StaticDimer terminates succesfully after $(nit) iterations")
         end
         return x, v, log
      end
      if numdE > maxnumdE
         if verbose >= 1
            println("StaticDimer terminates unsuccesfully due to numdE >= $(maxnumdE)")
         end
         return x, v, log
      end

      # translation step
      p_trans = - (P \ dE0) + 2.0 * dot(v, dE0) * v
      x += a_trans * p_trans
      # rotation step
      if precon_rot
         p_rot = - (P \ Hv) + dot(Hv, v) * v
      else
         p_rot = q_rot
      end
      v += a_rot * p_rot
   end

   error("why am I here?")
end

function run!{T}(method::BBDimer, E, dE, x0::Vector{T}, v0::Vector{T})

   # read all the parameters
   @unpack a0_trans, a0_rot, tol_trans, tol_rot, maxnumdE, len,
            precon_prep!, verbose, precon_rot, rescale_v, ls = method
   P0=method.precon
   # initialise variables
   x, v = copy(x0), copy(v0)
   β, γ = -1.0, -1.0   # this tells the loop that they haven't been initialised
   numdE, numE = 0, 0
   log = DimerLog()
   p_trans_old = zeros(length(x0))
   p_rot_old = zeros(length(x0))
   # and just start looping
   if verbose >= 2
      @printf(" nit |  |∇E|_∞    |∇R|_∞        λ         β         γ \n")
      @printf("-----|------------------------------------------------\n")
   end
   nit = 0
   while true
      nit += 1

      if any(isnan, v) || any(isnan, x)
         error("""BBDimer method has encountered NANs. This can happen with
                  BB type step size selection, which can be fast but is sometimes
                  unstable. To prevent this, use try different initial conditions
                  or use a different saddle-search method.""")
      end

      # normalise v
      P0 = precon_prep!(P0, x)
      v /= sqrt(dot(v, P0, v))
      # evaluate gradients, and more stuff
      dEm = dE(x - len/2 * v)
      dEp = dE(x + len/2 * v)
      numdE += 2
      Hv = (dEp - dEm) / len
      dE0 = 0.5 * (dEp + dEm)

      # NEWTON TYPE RESCALING IN v DIRECTION
      if rescale_v
         P = PreconSMW(P0, v, abs(dot(Hv, v)) - 1.0)
         v /= sqrt(dot(v, P, v))
      else
         P = P0
      end

      # translation and rotation residual, store history
      res_trans = vecnorm(dE0, Inf)
      λ = dot(v, Hv)
      q_rot = - Hv + λ * (P * v)
      res_rot = vecnorm(q_rot, Inf)
      push!(log, numE, numdE, res_trans, res_rot)
      if verbose >= 2
         @printf("%4d | %1.2e  %1.2e  %1.2e  %1.2e  %1.2e \n",
               nit, res_trans, res_rot, λ, β, γ)
      end

      # check whether to terminate
      if res_trans <= tol_trans && res_rot <= tol_rot
         if verbose >= 1
            println("BBDimer terminates succesfully after $(nit) iterations")
         end
         return x, v, log
      end
      if numdE > maxnumdE
         if verbose >= 1
            println("BBDimer terminates unsuccesfully due to numdE >= $(maxnumdE)")
         end
         return x, v, log
      end

      # compute the two "search" directions
      p_trans = - (P \ dE0) + 2.0 * dot(v, dE0) * v
      if precon_rot
         p_rot = - (P \ Hv - λ * v)
      else
         p_rot = q_rot
      end

      # initial step-sizes guess
      if nit == 1     # first iteration
         β = a0_trans
         γ = a0_rot
      else             # subsequent iterations: BB step
         # Δx, Δv has already been computed (see end of loop)
         Δg = p_trans - p_trans_old
         Δd = p_rot - p_rot_old
         β = abs( dot(Δx, P, Δg) / dot(Δg, P, Δg) )
         γ = abs( dot(Δv, P, Δd) / dot(Δd, P, Δd) )
      end

      # perform linesearch on the translation step
      F_trans = xx -> localmerit(xx, x, v, len, dE0, λ, E)
      β, numE_ls, _ = linesearch!(ls, F_trans, F_trans(x), - dot(p_trans, P, p_trans),
                                    x, p_trans, β)
      numE += numE_ls

      # perform line-search on rotation step
      E0 = E(x);  numE += 1
      F_rot = vv -> rayleigh(vv, x, len, E0, E, P)
      γ, numE_ls, _ = linesearch!(ls, F_rot, F_rot(v), -dot(p_rot, P, p_rot),
                                    v, p_rot, γ)
      numE += numE_ls

      if isnan(β) || isnan(γ)
         if verbose >= 1
            println("BBDimer terminates unsuccesfully due to unsuccesful linesearch")
         end
         return x, v, log
      end

      x_old, v_old = x, v
      # translation step
      x += β * p_trans
      # rotation step
      v += γ * p_rot

      # remember the change in x and v, also remember the old p_trans and p_rot
      Δx = β * p_trans
      Δv = γ * p_rot
      copy!(p_trans_old, p_trans)
      copy!(p_rot_old, p_rot)
   end
   error("why am I here?")
end

function dimer_ode(z, dE, P, precon_prep!, len)
   n = length(z) ÷ 2
   x, v = z[1:n], z[n+1:end]
   P = precon_prep!(P, x)
   # evaluate gradients, and more stuff
   dE0 = dE(x)
   dEv = dE(x + len * v)
   Hv = (dEv - dE0) / len

   p_trans = - (P \ dE0) + 2.0 * dot(v, dE0) * v
   p_rot = - (P \ Hv) + dot(Hv, v) * v
   F = [p_trans; p_rot]
   return F, norm(F, Inf), 2, dot
end

function dimer_project(z, P, precon_prep!)
   n = length(z) ÷ 2
   x, v = z[1:n], z[n+1:end]
   P = precon_prep!(P, x)
   v /= sqrt(dot(v, P, v))
   znew = copy(z)
   znew[n+1:end] = v
   return znew
end

function run!(method::ODEDimer, E, dE, x0::Vector, v0::Vector)

   # read all the parameters
   @unpack tol_trans, tol_rot, maxnumdE, len,
            precon_prep!, verbose, precon_rot, rescale_v, ode = method
   P0=method.precon

   # initial condition
   n = length(x0)
   z0 = [x0; v0]
   # nonlinear system
   F = (z, P, nit) -> dimer_ode(z, dE, P0, precon_prep!, len)
   # projection (normalisation) step
   G = (z, P) -> dimer_project(z, P0, precon_prep!)
   # initialise a log
   log = PathLog()

   # run the ODE solver
   zout, log = odesolve(ode,
                        F,
                        z0,   # initial condition
                        log,  # store iteration information in this log
                        g = G,
                        maxnit = maxnumdE,
                        verbose = verbose,
                        tol = min(tol_trans, tol_rot))
   z = zout[end]
   return z[1:n], z[n+1:end], log
end

function dimer_jacobian(z, dE, ddE, P, precon_prep!, len)

   n = length(z) ÷ 2
   x, v = z[1:n], z[n+1:end]
   P = precon_prep!(P, x)

   # evaluate gradients, etc
   dE0 = dE(x)
   dEv = dE(x + len * v)
   Hv = (dEv - dE0) / len
   ddE0 = ddE(x)
   ddEv = ddE(x + len * v)

   Fxx = ∂ₓFˣ(x, v, ddE0)
   Fxv = ∂ᵥFˣ(x, v, dE0)
   Fvx = ∂ₓFᵛ(x, v, len, ddE0, ddEv)
   Fvv = ∂ᵥFᵛ(x, v, Hv, ddEv)

   dF = [Fxx Fxv; Fvx Fvv]

   return dF
end

function run!(method::AccelDimer, E, dE, ddE, x0::Vector, v0::Vector)

   # read all the parameters
   @unpack tol_trans, tol_rot, maxnumdE, len,
            precon_prep!, verbose, precon_rot, rescale_v,
            a0, b, fd_scheme, redistrib, = method
   accel = momentum_descent(h = a0, b = b , fd_scheme = fd_scheme,
                              redistrib = redistrib)
   P0=method.precon

   # initial condition
   n = length(x0)
   z0 = [x0; v0]
   # nonlinear system
   F = (z, P, nit) -> dimer_ode(z, dE, P0, precon_prep!, len)
   dF = (z, P) -> dimer_jacobian(z, dE, ddE, P0, precon_prep!, len)
   # projection (normalisation) step
   G = (z, P) -> dimer_project(z, P0, precon_prep!)
   # initialise a log
   log = PathLog()

   # run the ODE solver
   zout, log = odesolve(accel,
                        F,
                        dF,
                        z0,   # initial condition
                        log,  # store iteration information in this log
                        verbose = verbose,
                        g = G,
                        maxnit = maxnumdE,
                        tol = min(tol_trans, tol_rot))
   z = zout[end]
   return z[1:n], z[n+1:end], log
end


text:
include("../src/robin_dict.jl")

@testset "Constructors" begin
    h1 = RobinDict()
    @test length(h1) == 0
    @test isempty(h1) == true
    @test h1.idxfloor == 0
    @test length(h1.keys) == 16
    @test length(h1.vals) == 16
    @test length(h1.hashes) == 16
    @test eltype(h1) == Pair{Any, Any}
    @test keytype(h1) == Any
    @test valtype(h1) == Any
end

@testset "RobinDict" begin
    h = RobinDict()
    for i=1:10000
        h[i] = i+1
    end
    for i=1:10000
        @test (h[i] == i+1)
    end
    for i=1:2:10000
        delete!(h, i)
    end
    for i=1:2:10000
        h[i] = i+1
    end
    for i=1:10000
        @test (h[i] == i+1)
    end
    for i=1:10000
        delete!(h, i)
    end
    @test isempty(h)
    h[77] = 100
    @test h[77] == 100
    for i=1:10000
        h[i] = i+1
    end
    for i=1:2:10000
        delete!(h, i)
    end
    for i=10001:20000
        h[i] = i+1
    end
    for i=2:2:10000
        @test h[i] == i+1
    end
    for i=10000:20000
        @test h[i] == i+1
    end
    h = RobinDict{Any,Any}("a" => 3)
    @test h["a"] == 3
    h["a","b"] = 4
    @test h["a","b"] == h[("a","b")] == 4
    h["a","b","c"] = 4
    @test h["a","b","c"] == h[("a","b","c")] == 4

    @testset "eltype, keytype and valtype" begin
        @test eltype(h) == Pair{Any,Any}
        @test keytype(h) == Any
        @test valtype(h) == Any

        td = RobinDict{AbstractString,Float64}()
        @test eltype(td) == Pair{AbstractString,Float64}
        @test keytype(td) == AbstractString
        @test valtype(td) == Float64
        @test keytype(Dict{AbstractString,Float64}) === AbstractString
        @test valtype(Dict{AbstractString,Float64}) === Float64
    end
    # test rethrow of error in ctor
    @test_throws DomainError RobinDict((sqrt(p[1]), sqrt(p[2])) for p in zip(-1:2, -1:2))
end

@testset "RobinDict on pairs" begin
	let x = RobinDict(3=>3, 5=>5, 8=>8, 6=>6)
	    pop!(x, 5)
	    for k in keys(x)
	        RobinDict{Int,Int}(x)
	        @test k in [3, 8, 6]
	    end
	end

    let y = RobinDict{Any, Int}(3=>3, 5=>5, "8"=>8, 6=>6)
        pop!(y, "8")
        for k in keys(y)
            RobinDict{Int,Int}(y)
            @test k in [3, 5, 6]
        end
    end

    d = @inferred RobinDict(Pair(1,1), Pair(2,2), Pair(3,3))
    @test isa(d, RobinDict)
    @test d == RobinDict(1=>1, 2=>2, 3=>3)
    @test eltype(d) == Pair{Int,Int}
end

@testset "KeyError" begin
    let z = RobinDict()
        get_KeyError = false
        try
            z["a"]
        catch _e123_
            get_KeyError = isa(_e123_,KeyError)
        end
        @test get_KeyError
    end
end

@testset "Filter function" begin
    _d = RobinDict("a"=>0)
    @test isa([k for k in filter(x->length(x)==1, collect(keys(_d)))], Vector{String})

    h = RobinDict{Int, Int}()
    for i in 1:100
        h[i] = i
    end
    filter!(x->isodd(x.first), h)
    for i in 1:2:100
        @test h[i] == i
    end
end

@testset "typeof" begin
    d = RobinDict(((1, 2), (3, 4)))
    @test d[1] === 2
    @test d[3] === 4
    d2 = RobinDict(1 => 2, 3 => 4)
    d3 = RobinDict((1 => 2, 3 => 4))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == RobinDict{Int,Int}

    d = RobinDict(((1, 2), (3, "b")))
    @test d[1] === 2
    @test d[3] == "b"
    d2 = RobinDict(1 => 2, 3 => "b")
    d3 = RobinDict((1 => 2, 3 => "b"))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == RobinDict{Int,Any}

    d = RobinDict(((1, 2), ("a", 4)))
    @test d[1] === 2
    @test d["a"] === 4
    d2 = RobinDict(1 => 2, "a" => 4)
    d3 = RobinDict((1 => 2, "a" => 4))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == RobinDict{Any,Int}

    d = RobinDict(((1, 2), ("a", "b")))
    @test d[1] === 2
    @test d["a"] == "b"
    d2 = RobinDict(1 => 2, "a" => "b")
    d3 = RobinDict((1 => 2, "a" => "b"))
    @test d == d2 == d3
    @test typeof(d) == typeof(d2) == typeof(d3) == RobinDict{Any,Any}
end

@testset "type of RobinDict constructed from varargs of Pairs" begin
    @test RobinDict(1=>1, 2=>2.0) isa RobinDict{Int,Real}
    @test RobinDict(1=>1, 2.0=>2) isa RobinDict{Real,Int}
    @test RobinDict(1=>1.0, 2.0=>2) isa RobinDict{Real,Real}

    for T in (Nothing, Missing)
        @test RobinDict(1=>1, 2=>T()) isa RobinDict{Int,Union{Int,T}}
        @test RobinDict(1=>T(), 2=>2) isa RobinDict{Int,Union{Int,T}}
        @test RobinDict(1=>1, T()=>2) isa RobinDict{Union{Int,T},Int}
        @test RobinDict(T()=>1, 2=>2) isa RobinDict{Union{Int,T},Int}
    end
end

@testset "equality" for eq in (isequal, ==)
    @test  eq(RobinDict(), RobinDict())
    @test  eq(RobinDict(1 => 1), RobinDict(1 => 1))
    @test !eq(RobinDict(1 => 1), RobinDict())
    @test !eq(RobinDict(1 => 1), RobinDict(1 => 2))
    @test !eq(RobinDict(1 => 1), RobinDict(2 => 1))

    # Generate some data to populate dicts to be compared
    data_in = [ (rand(1:1000), randstring(2)) for _ in 1:1001 ]

    # Populate the first dict
    d1 = RobinDict{Int, AbstractString}()
    for (k, v) in data_in
        d1[k] = v
    end
    data_in = collect(d1)
    # shuffle the data
    for i in 1:length(data_in)
        j = rand(1:length(data_in))
        data_in[i], data_in[j] = data_in[j], data_in[i]
    end
    # Inserting data in different (shuffled) order should result in
    # equivalent dict.
    d2 = RobinDict{Int, AbstractString}()
    for (k, v) in data_in
        d2[k] = v
    end

    @test eq(d1, d2)
    d3 = copy(d2)
    d4 = copy(d2)
    # Removing an item gives different dict
    delete!(d1, data_in[rand(1:length(data_in))][1])
    @test !eq(d1, d2)
    # Changing a value gives different dict
    d3[data_in[rand(1:length(data_in))][1]] = randstring(3)
    !eq(d1, d3)
    # Adding a pair gives different dict
    d4[1001] = randstring(3)
    @test !eq(d1, d4)

    @test eq(RobinDict(), sizehint!(RobinDict(),96))

    # Dictionaries of different types
    @test_throws MethodError eq(RobinDict(1 => 2), RobinDict("dog" => "bone"))
    @test eq(RobinDict{Int,Int}(), RobinDict{AbstractString,AbstractString}())
end

@testset "equality special cases" begin
    @test RobinDict(1=>0.0) == RobinDict(1=>-0.0)
    @test !isequal(RobinDict(1=>0.0), RobinDict(1=>-0.0))

    @test RobinDict(0.0=>1) != RobinDict(-0.0=>1)
    @test !isequal(RobinDict(0.0=>1), RobinDict(-0.0=>1))

    @test RobinDict(1=>NaN) != RobinDict(1=>NaN)
    @test isequal(RobinDict(1=>NaN), RobinDict(1=>NaN))

    @test RobinDict(NaN=>1) == RobinDict(NaN=>1)
    @test isequal(RobinDict(NaN=>1), RobinDict(NaN=>1))

    @test ismissing(RobinDict(1=>missing) == RobinDict(1=>missing))
    @test isequal(RobinDict(1=>missing), RobinDict(1=>missing))

    @test RobinDict(missing=>1) == RobinDict(missing=>1)
    @test isequal(RobinDict(missing=>1), RobinDict(missing=>1))
end

@testset "get!" begin
    f(x) = x^2
    d = RobinDict(8=>19)
    @test get!(d, 8, 5) == 19
    @test get!(d, 19, 2) == 2

    @test get!(d, 42) do  # d is updated with f(2)
        f(2)
    end == 4

    @test get!(d, 42) do  # d is not updated
        f(200)
    end == 4

    @test get(d, 13) do   # d is not updated
        f(4)
    end == 16

    @test d == RobinDict(8=>19, 19=>2, 42=>4)
end

@testset "push!" begin
    d = RobinDict()
    @test push!(d, 'a' => 1) === d
    @test d['a'] == 1
    @test push!(d, 'b' => 2, 'c' => 3) === d
    @test d['b'] == 2
    @test d['c'] == 3
    @test push!(d, 'd' => 4, 'e' => 5, 'f' => 6) === d
    @test d['d'] == 4
    @test d['e'] == 5
    @test d['f'] == 6
    @test length(d) == 6
end

@testset "pop!" begin
    d = RobinDict(1=>2, 3=>4)
    @test pop!(d, 1) == 2
    @test_throws KeyError pop!(d, 1)
    @test pop!(d, 1, 0) == 0
    @test pop!(d) == (3=>4)
    @test_throws ArgumentError pop!(d)
end

@testset "keys as a set" begin
    d = RobinDict(1=>2, 3=>4)
    @test keys(d) isa AbstractSet
    @test empty(keys(d)) isa AbstractSet
    let i = keys(d) ∩ Set([1,2])
        @test i isa AbstractSet
        @test i == Set([1])
    end
    @test Set(string(k) for k in keys(d)) == Set(["1","3"])
end

@testset "find" begin
    @test findall(isequal(1), RobinDict(:a=>1, :b=>2)) == [:a]
    @test sort(findall(isequal(1), RobinDict(:a=>1, :b=>1))) == [:a, :b]
    @test isempty(findall(isequal(1), RobinDict()))
    @test isempty(findall(isequal(1), RobinDict(:a=>2, :b=>3)))

    @test findfirst(isequal(1), RobinDict(:a=>1, :b=>2)) == :a
    @test findfirst(isequal(1), RobinDict(:a=>1, :b=>1, :c=>3)) in (:a, :b)
    @test findfirst(isequal(1), RobinDict()) === nothing
    @test findfirst(isequal(1), RobinDict(:a=>2, :b=>3)) === nothing
end

@testset "haskey" begin
    h = RobinDict(1=>2, 2=>3)
    @test haskey(h, 1) == true
    @test haskey(h, 2) == true
    @test haskey(h, 3) == false
    @test_throws MethodError haskey(h, "1")
end

@testset "getkey" begin
    h = RobinDict(1=>2, 3 => 6, 5=>10)
    @test getkey(h, 1, 7) == 1
    @test getkey(h, 4, 6) == 6
    @test_throws MethodError getkey(h, "1", 8)
end

@testset "empty" begin
    h = RobinDict()
    for i=1:1000
        h[i] = i+1
    end
    length0 = length(h.hashes)
    empty!(h)
    @test h.count == 0
    @test h.idxfloor == 0
    @test length(h.hashes) == length(h.keys) == length(h.vals) == length0
end

@testset "ArgumentError" begin
    @test_throws ArgumentError RobinDict(0)
    @test_throws ArgumentError RobinDict([1])
    @test_throws ArgumentError RobinDict([(1,2),0])
end

@testset "empty tuple" begin
    h = RobinDict(())
    @test length(h) == 0
end

@testset "get_idxfloor" begin
    h = RobinDict()
    @test get_idxfloor(h) == 0

    h["a"] = 1
    h[2] = "b"
    @test h.idxfloor == get_idxfloor(h)
    pop!(h)
    @test h.idxfloor == get_idxfloor(h)
    pop!(h)
    @test h.idxfloor == get_idxfloor(h) == 0
end

@testset "invariants" begin
    h1 = RobinDict{Int, Int}()
    for i in 1:300
        h1[i] = i
    end

    for i in 1:length(h1.keys)
        if isslotfilled(h1, i)
            @test hash_key(h1.keys[i]) == h1.hashes[i]
        end
    end

    h2 = RobinDict{Float64, Float64}()
    for i in 1:300
        h2[rand()] = rand()
    end

    for i in 1:length(h2.keys)
        if isslotfilled(h2, i)
            @test hash_key(h2.keys[i]) == h2.hashes[i]
        end
    end

    h3 = RobinDict{String, Int}()
    for i in 1:300
        h3[randstring()] = i
    end

    for i in 1:length(h3.keys)
        if isslotfilled(h3, i)
            @test hash_key(h3.keys[i]) == h3.hashes[i]
        end
    end

    max_disp = 0
    function check_invariants(h::RobinDict)
        cnt = 0
        min_idx = 0
        sz = length(h.keys)
        for i=1:length(h.keys)
            isslotfilled(h, i) || continue
            (min_idx == 0) && (min_idx = i)
            @assert hash_key(h.keys[i]) == h.hashes[i]
            @assert (h.hashes[i] & 0x80000000) != 0
            cnt += 1
            @assert typeof(h.hashes[i]) == UInt32
            des_ind = desired_index(h.hashes[i], sz)
            pos_diff = 0
            if (i >= des_ind)
                pos_diff = i - des_ind
            else
                pos_diff = sz - des_ind + i
            end
            dist = calculate_distance(h, i)
            @assert dist == pos_diff
            max_disp = max(max_disp, dist)
            distlast = (i != 1) ? isslotfilled(h, i-1) ? calculate_distance(h, i-1) : 0 : isslotfilled(h, sz) ? calculate_distance(h, sz) : 0
            @assert dist <= distlast + 1
        end
        @assert h.idxfloor == min_idx
        @assert cnt == length(h)
    end

    h = RobinDict()
    for i = 1:10000
        h[i] = i+1
    end
    check_invariants(h)
end


text:
using NetCDF, Plots, Printf

filename = "convection_day"

base_dir = pwd()*"/les_data/"*filename
filename = base_dir*"/"*filename*".nc"

filename = pwd() * sub_directory * "_middepth.nc"
ncinfo(filename)
x = Array(NetCDF.open(filename, "yC"))
y = Array(NetCDF.open(filename, "yC"))
z = Array(NetCDF.open(filename, "zC"))
t = Array(NetCDF.open(filename, "time"))
sim_day = t ./ 86400
z = (z[2:end] + z[1:end-1]) / 2
b = NetCDF.open(filename, "b")
u = NetCDF.open(filename, "u")
v = NetCDF.open(filename, "v")
w = NetCDF.open(filename, "w")
##

# T = b./(g*α)
# @. T /= 9.81 * 2e-4
T = b./(9.81 * 2e-4)

cmax = maximum(Array(b[2:end-1,2:end-1, 1, 1:10:end]))
cmin = minimum(Array(b[:,:, 1, 1:10:end]))
clims = (cmin, cmax)
anim = @animate for i in 1:1:length(b[1, 1, 1, :])
    b_array = Array(b[:, :, 1, i])
    # p1 = contourf(y, z, b_array, fill = true, linewidth = 0, color = :ocean, clim = (-0.0013871555898403098, -3.3776441214941526e-6))
    day_label = @sprintf("%.2f ", sim_day[i])
    p1 = contourf(x, y, b_array',
    color = :ocean, title = "Buoyancy at Middepth at day " * day_label ,
     xlabel = "Meridional [m]", ylabel = "Depth [m]"
     , clims = clims)
end
gif(anim, filename *".gif", fps = 15)

max_stars_repo_name: adelinehillier/LearnConvection

max_issues_repo_name: adelinehillier/LearnConvection

max_forks_repo_name: adelinehillier/LearnConvection


text:
import Base: show, +, *, ==, size, getindex, setindex!, show, intersect, symdiff, union, push!, append!, isless

using Statistics

# vectors with elements that are sorted and unique
include("suvector.jl")
# pauli group elements
include("pauli.jl")
# hamiltonian and rg
include("rghamiltonian.jl")

max_stars_repo_name: aormorningstar/SBRG

max_issues_repo_name: aormorningstar/SBRG

max_forks_repo_name: aormorningstar/SBRG


text:
using SparseArrays, ENN

pn=TPN(
    [:is_trig1, :is_on1, :is_off1, :is_trig2, :is_on2, :is_off2, :counter],
    [:in1, :in2, :on1, :off1, :on2, :off2],
    sparse([
        0; 0; 0;  0; 0; 0;  0 ;; # fire first input
        0; 0; 0;  0; 0; 0;  0 ;; # fire second input
        1; 0; 1;  0; 0; 0;  0 ;; # turn 1 on
        0; 1; 0;  0; 0; 0;  1 ;; # turn 1 off
        0; 0; 0;  1; 0; 1;  0 ;; # turn 2 on
        0; 0; 0;  0; 1; 0;  1    # turn 2 off
    ]),
    sparse([
        1; 0; 0;  0; 0; 0;  0 ;; # fire first input
        0; 0; 0;  1; 0; 0;  0 ;; # fire second input
       -1; 1;-1;  0; 0; 0;  1 ;; # turn 1 on and increase counter
        0;-1; 1;  0; 0; 0; -1 ;; # turn 1 off and decrease counter
        0; 0; 0; -1; 1;-1;  1 ;; # turn 2 on and increase counter
        0; 0; 0;  0;-1; 1; -1    # turn 2 off and decrease counter
    ]),
    [0.0, 0.0, 0.0, 1.0, 0.0, 1.0],
    [Inf, Inf, 0.0, 1.0, 0.0, 1.0],
    [0, 0, 1, 0, 0, 1, 0]
)

f,ax,p = graphplot(pn); f


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*--------------------------------------------------------------------
#* Model of basic streams
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt and Rafael de P. Soares
#* $Id$
#*---------------------------------------------------------------------
type simple_sink2
	simple_sink2()=begin
		new(
			stream(Dict{Symbol,Any}(
				:Brief=>"Inlet Stream",
				:PosX=>1,
				:PosY=>0.5308,
				:Protected=>true,
				:Symbol=>"_{in}"
			)),
			[:Inlet,]
		)
	end
	Inlet::stream
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export simple_sink2
function atributes(in::simple_sink2,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Sink2"
	fields[:Brief]="Simple material stream sink"
	fields[:Info]="
	This model should be used for boundary streams when no additional
	information about the stream is desired.
	"
	drive!(fields,_)
	return fields
end
simple_sink2(_::Dict{Symbol,Any})=begin
	newModel=simple_sink2()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
# Julia wrapper for header: /Users/wookyoung/arrayfire/forge/include/forge.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function fg_update_vertex_buffer(pBufferId::UInt32,pBufferSize::UInt32,pBufferData)
    ccall((:fg_update_vertex_buffer,libforge),fg_err,(UInt32,Csize_t,Ptr{Void}),pBufferId,pBufferSize,pBufferData)
end

function fg_update_pixel_buffer(pBufferId::UInt32,pBufferSize::UInt32,pBufferData)
    ccall((:fg_update_pixel_buffer,libforge),fg_err,(UInt32,Csize_t,Ptr{Void}),pBufferId,pBufferSize,pBufferData)
end

function fg_finish()
    ccall((:fg_finish,libforge),fg_err,())
end

function fg_get_last_error(msg,len)
    ccall((:fg_get_last_error,libforge),Void,(Ptr{Cstring},Ptr{Cint}),msg,len)
end

function fg_err_to_string(err::fg_err)
    ccall((:fg_err_to_string,libforge),Cstring,(fg_err,),err)
end

function fg_create_font(pFont)
    ccall((:fg_create_font,libforge),fg_err,(Ptr{fg_font},),pFont)
end

function fg_destroy_font(pFont::fg_font)
    ccall((:fg_destroy_font,libforge),fg_err,(fg_font,),pFont)
end

function fg_load_font_file(pFont::fg_font,pFileFullPath)
    ccall((:fg_load_font_file,libforge),fg_err,(fg_font,Cstring),pFont,pFileFullPath)
end

function fg_load_system_font(pFont::fg_font,pFontName)
    ccall((:fg_load_system_font,libforge),fg_err,(fg_font,Cstring),pFont,pFontName)
end

function fg_create_image(pImage,pWidth::UInt32,pHeight::UInt32,pFormat::fg_channel_format,pType::fg_dtype)
    ccall((:fg_create_image,libforge),fg_err,(Ptr{fg_image},UInt32,UInt32,fg_channel_format,fg_dtype),pImage,pWidth,pHeight,pFormat,pType)
end

function fg_destroy_image(pImage::fg_image)
    ccall((:fg_destroy_image,libforge),fg_err,(fg_image,),pImage)
end

function fg_set_image_alpha(pImage::fg_image,pAlpha::Cfloat)
    ccall((:fg_set_image_alpha,libforge),fg_err,(fg_image,Cfloat),pImage,pAlpha)
end

function fg_set_image_aspect_ratio(pImage::fg_image,pKeepRatio::Bool)
    ccall((:fg_set_image_aspect_ratio,libforge),fg_err,(fg_image,Bool),pImage,pKeepRatio)
end

function fg_get_image_width(pOut,pImage::fg_image)
    ccall((:fg_get_image_width,libforge),fg_err,(Ptr{UInt32},fg_image),pOut,pImage)
end

function fg_get_image_height(pOut,pImage::fg_image)
    ccall((:fg_get_image_height,libforge),fg_err,(Ptr{UInt32},fg_image),pOut,pImage)
end

function fg_get_image_pixelformat(pOut,pImage::fg_image)
    ccall((:fg_get_image_pixelformat,libforge),fg_err,(Ptr{fg_channel_format},fg_image),pOut,pImage)
end

function fg_get_image_type(pOut,pImage::fg_image)
    ccall((:fg_get_image_type,libforge),fg_err,(Ptr{fg_dtype},fg_image),pOut,pImage)
end

function fg_get_pixel_buffer(pOut,pImage::fg_image)
    ccall((:fg_get_pixel_buffer,libforge),fg_err,(Ptr{UInt32},fg_image),pOut,pImage)
end

function fg_get_image_size(pOut,pImage::fg_image)
    ccall((:fg_get_image_size,libforge),fg_err,(Ptr{UInt32},fg_image),pOut,pImage)
end

function fg_render_image(pWindow::fg_window,pImage::fg_image,pX::Cint,pY::Cint,pWidth::Cint,pHeight::Cint)
    ccall((:fg_render_image,libforge),fg_err,(fg_window,fg_image,Cint,Cint,Cint,Cint),pWindow,pImage,pX,pY,pWidth,pHeight)
end

function fg_create_plot(pPlot,pNPoints::UInt32,pType::fg_dtype,pChartType::fg_chart_type,pPlotType::fg_plot_type,pMarkerType::fg_marker_type)
    ccall((:fg_create_plot,libforge),fg_err,(Ptr{fg_plot},UInt32,fg_dtype,fg_chart_type,fg_plot_type,fg_marker_type),pPlot,pNPoints,pType,pChartType,pPlotType,pMarkerType)
end

function fg_destroy_plot(pPlot::fg_plot)
    ccall((:fg_destroy_plot,libforge),fg_err,(fg_plot,),pPlot)
end

function fg_set_plot_color(pPlot::fg_plot,pRed::Cfloat,pGreen::Cfloat,pBlue::Cfloat,pAlpha::Cfloat)
    ccall((:fg_set_plot_color,libforge),fg_err,(fg_plot,Cfloat,Cfloat,Cfloat,Cfloat),pPlot,pRed,pGreen,pBlue,pAlpha)
end

function fg_set_plot_legend(pPlot::fg_plot,pLegend)
    ccall((:fg_set_plot_legend,libforge),fg_err,(fg_plot,Cstring),pPlot,pLegend)
end

function fg_set_plot_marker_size(pPlot::fg_plot,pMarkerSize::Cfloat)
    ccall((:fg_set_plot_marker_size,libforge),fg_err,(fg_plot,Cfloat),pPlot,pMarkerSize)
end

function fg_get_plot_vertex_buffer(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_vertex_buffer,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_color_buffer(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_color_buffer,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_alpha_buffer(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_alpha_buffer,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_radii_buffer(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_radii_buffer,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_vertex_buffer_size(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_vertex_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_color_buffer_size(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_color_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_alpha_buffer_size(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_alpha_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_get_plot_radii_buffer_size(pOut,pPlot::fg_plot)
    ccall((:fg_get_plot_radii_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_plot),pOut,pPlot)
end

function fg_create_surface(pSurface,pXPoints::UInt32,pYPoints::UInt32,pType::fg_dtype,pPlotType::fg_plot_type,pMarkerType::fg_marker_type)
    ccall((:fg_create_surface,libforge),fg_err,(Ptr{fg_surface},UInt32,UInt32,fg_dtype,fg_plot_type,fg_marker_type),pSurface,pXPoints,pYPoints,pType,pPlotType,pMarkerType)
end

function fg_destroy_surface(pSurface::fg_surface)
    ccall((:fg_destroy_surface,libforge),fg_err,(fg_surface,),pSurface)
end

function fg_set_surface_color(pSurface::fg_surface,pRed::Cfloat,pGreen::Cfloat,pBlue::Cfloat,pAlpha::Cfloat)
    ccall((:fg_set_surface_color,libforge),fg_err,(fg_surface,Cfloat,Cfloat,Cfloat,Cfloat),pSurface,pRed,pGreen,pBlue,pAlpha)
end

function fg_set_surface_legend(pSurface::fg_surface,pLegend)
    ccall((:fg_set_surface_legend,libforge),fg_err,(fg_surface,Cstring),pSurface,pLegend)
end

function fg_get_surface_vertex_buffer(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_vertex_buffer,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_get_surface_color_buffer(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_color_buffer,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_get_surface_alpha_buffer(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_alpha_buffer,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_get_surface_vertex_buffer_size(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_vertex_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_get_surface_color_buffer_size(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_color_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_get_surface_alpha_buffer_size(pOut,pSurface::fg_surface)
    ccall((:fg_get_surface_alpha_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_surface),pOut,pSurface)
end

function fg_create_vector_field(pField,pNPoints::UInt32,pType::fg_dtype,pChartType::fg_chart_type)
    ccall((:fg_create_vector_field,libforge),fg_err,(Ptr{fg_vector_field},UInt32,fg_dtype,fg_chart_type),pField,pNPoints,pType,pChartType)
end

function fg_destroy_vector_field(pField::fg_vector_field)
    ccall((:fg_destroy_vector_field,libforge),fg_err,(fg_vector_field,),pField)
end

function fg_set_vector_field_color(pField::fg_vector_field,pRed::Cfloat,pGreen::Cfloat,pBlue::Cfloat,pAlpha::Cfloat)
    ccall((:fg_set_vector_field_color,libforge),fg_err,(fg_vector_field,Cfloat,Cfloat,Cfloat,Cfloat),pField,pRed,pGreen,pBlue,pAlpha)
end

function fg_set_vector_field_legend(pField::fg_vector_field,pLegend)
    ccall((:fg_set_vector_field_legend,libforge),fg_err,(fg_vector_field,Cstring),pField,pLegend)
end

function fg_get_vector_field_vertex_buffer(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_vertex_buffer,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_color_buffer(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_color_buffer,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_alpha_buffer(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_alpha_buffer,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_direction_buffer(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_direction_buffer,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_vertex_buffer_size(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_vertex_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_color_buffer_size(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_color_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_alpha_buffer_size(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_alpha_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_get_vector_field_direction_buffer_size(pOut,pField::fg_vector_field)
    ccall((:fg_get_vector_field_direction_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_vector_field),pOut,pField)
end

function fg_create_histogram(pHistogram,pNBins::UInt32,pDataType::fg_dtype)
    ccall((:fg_create_histogram,libforge),fg_err,(Ptr{fg_histogram},UInt32,fg_dtype),pHistogram,pNBins,pDataType)
end

function fg_destroy_histogram(pHistogram::fg_histogram)
    ccall((:fg_destroy_histogram,libforge),fg_err,(fg_histogram,),pHistogram)
end

function fg_set_histogram_color(pHistogram::fg_histogram,pRed::Cfloat,pGreen::Cfloat,pBlue::Cfloat,pAlpha::Cfloat)
    ccall((:fg_set_histogram_color,libforge),fg_err,(fg_histogram,Cfloat,Cfloat,Cfloat,Cfloat),pHistogram,pRed,pGreen,pBlue,pAlpha)
end

function fg_set_histogram_legend(pHistogram::fg_histogram,pLegend)
    ccall((:fg_set_histogram_legend,libforge),fg_err,(fg_histogram,Cstring),pHistogram,pLegend)
end

function fg_get_histogram_vertex_buffer(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_vertex_buffer,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_get_histogram_color_buffer(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_color_buffer,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_get_histogram_alpha_buffer(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_alpha_buffer,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_get_histogram_vertex_buffer_size(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_vertex_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_get_histogram_color_buffer_size(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_color_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_get_histogram_alpha_buffer_size(pOut,pHistogram::fg_histogram)
    ccall((:fg_get_histogram_alpha_buffer_size,libforge),fg_err,(Ptr{UInt32},fg_histogram),pOut,pHistogram)
end

function fg_create_chart(pHandle,pChartType::fg_chart_type)
    ccall((:fg_create_chart,libforge),fg_err,(Ptr{fg_chart},fg_chart_type),pHandle,pChartType)
end

function fg_destroy_chart(pHandle::fg_chart)
    ccall((:fg_destroy_chart,libforge),fg_err,(fg_chart,),pHandle)
end

function fg_set_chart_axes_titles(pHandle::fg_chart,pX,pY,pZ)
    ccall((:fg_set_chart_axes_titles,libforge),fg_err,(fg_chart,Cstring,Cstring,Cstring),pHandle,pX,pY,pZ)
end

function fg_set_chart_axes_limits(pHandle::fg_chart,pXmin::Cfloat,pXmax::Cfloat,pYmin::Cfloat,pYmax::Cfloat,pZmin::Cfloat,pZmax::Cfloat)
    ccall((:fg_set_chart_axes_limits,libforge),fg_err,(fg_chart,Cfloat,Cfloat,Cfloat,Cfloat,Cfloat,Cfloat),pHandle,pXmin,pXmax,pYmin,pYmax,pZmin,pZmax)
end

function fg_get_chart_axes_limits(pXmin,pXmax,pYmin,pYmax,pZmin,pZmax,pHandle::fg_chart)
    ccall((:fg_get_chart_axes_limits,libforge),fg_err,(Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat},Ptr{Cfloat},fg_chart),pXmin,pXmax,pYmin,pYmax,pZmin,pZmax,pHandle)
end

function fg_set_chart_legend_position(pHandle::fg_chart,pX::Cfloat,pY::Cfloat)
    ccall((:fg_set_chart_legend_position,libforge),fg_err,(fg_chart,Cfloat,Cfloat),pHandle,pX,pY)
end

function fg_add_image_to_chart(pImage,pHandle::fg_chart,pWidth::UInt32,pHeight::UInt32,pFormat::fg_channel_format,pType::fg_dtype)
    ccall((:fg_add_image_to_chart,libforge),fg_err,(Ptr{fg_image},fg_chart,UInt32,UInt32,fg_channel_format,fg_dtype),pImage,pHandle,pWidth,pHeight,pFormat,pType)
end

function fg_add_histogram_to_chart(pHistogram,pHandle::fg_chart,pNBins::UInt32,pType::fg_dtype)
    ccall((:fg_add_histogram_to_chart,libforge),fg_err,(Ptr{fg_histogram},fg_chart,UInt32,fg_dtype),pHistogram,pHandle,pNBins,pType)
end

function fg_add_plot_to_chart(pPlot,pHandle::fg_chart,pNPoints::UInt32,pType::fg_dtype,pPlotType::fg_plot_type,pMarkerType::fg_marker_type)
    ccall((:fg_add_plot_to_chart,libforge),fg_err,(Ptr{fg_plot},fg_chart,UInt32,fg_dtype,fg_plot_type,fg_marker_type),pPlot,pHandle,pNPoints,pType,pPlotType,pMarkerType)
end

function fg_add_surface_to_chart(pSurface,pHandle::fg_chart,pXPoints::UInt32,pYPoints::UInt32,pType::fg_dtype,pPlotType::fg_plot_type,pMarkerType::fg_marker_type)
    ccall((:fg_add_surface_to_chart,libforge),fg_err,(Ptr{fg_surface},fg_chart,UInt32,UInt32,fg_dtype,fg_plot_type,fg_marker_type),pSurface,pHandle,pXPoints,pYPoints,pType,pPlotType,pMarkerType)
end

function fg_add_vector_field_to_chart(pField,pHandle::fg_chart,pNPoints::UInt32,pType::fg_dtype)
    ccall((:fg_add_vector_field_to_chart,libforge),fg_err,(Ptr{fg_vector_field},fg_chart,UInt32,fg_dtype),pField,pHandle,pNPoints,pType)
end

function fg_render_chart(pWindow::fg_window,pChart::fg_chart,pX::Cint,pY::Cint,pWidth::Cint,pHeight::Cint)
    ccall((:fg_render_chart,libforge),fg_err,(fg_window,fg_chart,Cint,Cint,Cint,Cint),pWindow,pChart,pX,pY,pWidth,pHeight)
end

function fg_get_chart_type(pChartType,pChart::fg_chart)
    ccall((:fg_get_chart_type,libforge),fg_err,(Ptr{fg_chart_type},fg_chart),pChartType,pChart)
end

function fg_create_window(pWindow,pWidth::Cint,pHeight::Cint,pTitle,pShareWindow::fg_window,pInvisible::Bool)
    ccall((:fg_create_window,libforge),fg_err,(Ptr{fg_window},Cint,Cint,Cstring,fg_window,Bool),pWindow,pWidth,pHeight,pTitle,pShareWindow,pInvisible)
end

function fg_destroy_window(pWindow::fg_window)
    ccall((:fg_destroy_window,libforge),fg_err,(fg_window,),pWindow)
end

function fg_set_window_font(pWindow::fg_window,pFont::fg_font)
    ccall((:fg_set_window_font,libforge),fg_err,(fg_window,fg_font),pWindow,pFont)
end

function fg_set_window_title(pWindow::fg_window,pTitle)
    ccall((:fg_set_window_title,libforge),fg_err,(fg_window,Cstring),pWindow,pTitle)
end

function fg_set_window_position(pWindow::fg_window,pX::Cint,pY::Cint)
    ccall((:fg_set_window_position,libforge),fg_err,(fg_window,Cint,Cint),pWindow,pX,pY)
end

function fg_set_window_size(pWindow::fg_window,pWidth::UInt32,pHeight::UInt32)
    ccall((:fg_set_window_size,libforge),fg_err,(fg_window,UInt32,UInt32),pWindow,pWidth,pHeight)
end

function fg_set_window_colormap(pWindow::fg_window,pColorMap::fg_color_map)
    ccall((:fg_set_window_colormap,libforge),fg_err,(fg_window,fg_color_map),pWindow,pColorMap)
end

function fg_get_window_context_handle(pContext,pWindow::fg_window)
    ccall((:fg_get_window_context_handle,libforge),fg_err,(Ptr{Clonglong},fg_window),pContext,pWindow)
end

function fg_get_window_display_handle(pDisplay,pWindow::fg_window)
    ccall((:fg_get_window_display_handle,libforge),fg_err,(Ptr{Clonglong},fg_window),pDisplay,pWindow)
end

function fg_get_window_width(pWidth,pWindow::fg_window)
    ccall((:fg_get_window_width,libforge),fg_err,(Ptr{Cint},fg_window),pWidth,pWindow)
end

function fg_get_window_height(pHeight,pWindow::fg_window)
    ccall((:fg_get_window_height,libforge),fg_err,(Ptr{Cint},fg_window),pHeight,pWindow)
end

function fg_make_window_current(pWindow::fg_window)
    ccall((:fg_make_window_current,libforge),fg_err,(fg_window,),pWindow)
end

function fg_get_window_grid(pRows,pCols,pWindow::fg_window)
    ccall((:fg_get_window_grid,libforge),fg_err,(Ptr{Cint},Ptr{Cint},fg_window),pRows,pCols,pWindow)
end

function fg_hide_window(pWindow::fg_window)
    ccall((:fg_hide_window,libforge),fg_err,(fg_window,),pWindow)
end

function fg_show_window(pWindow::fg_window)
    ccall((:fg_show_window,libforge),fg_err,(fg_window,),pWindow)
end

function fg_close_window(pIsClosed,pWindow::fg_window)
    ccall((:fg_close_window,libforge),fg_err,(Ptr{Bool},fg_window),pIsClosed,pWindow)
end

function fg_draw_image(pWindow::fg_window,pImage::fg_image,pKeepAspectRatio::Bool)
    ccall((:fg_draw_image,libforge),fg_err,(fg_window,fg_image,Bool),pWindow,pImage,pKeepAspectRatio)
end

function fg_draw_chart(pWindow::fg_window,pChart::fg_chart)
    ccall((:fg_draw_chart,libforge),fg_err,(fg_window,fg_chart),pWindow,pChart)
end

function fg_setup_window_grid(pRows::Cint,pCols::Cint,pWindow::fg_window)
    ccall((:fg_setup_window_grid,libforge),fg_err,(Cint,Cint,fg_window),pRows,pCols,pWindow)
end

function fg_draw_image_to_cell(pWindow::fg_window,pRowId::Cint,pColId::Cint,pImage::fg_image,pTitle,pKeepAspectRatio::Bool)
    ccall((:fg_draw_image_to_cell,libforge),fg_err,(fg_window,Cint,Cint,fg_image,Cstring,Bool),pWindow,pRowId,pColId,pImage,pTitle,pKeepAspectRatio)
end

function fg_draw_chart_to_cell(pWindow::fg_window,pRowId::Cint,pColId::Cint,pChart::fg_chart,pTitle)
    ccall((:fg_draw_chart_to_cell,libforge),fg_err,(fg_window,Cint,Cint,fg_chart,Cstring),pWindow,pRowId,pColId,pChart,pTitle)
end

function fg_swap_window_buffers(pWindow::fg_window)
    ccall((:fg_swap_window_buffers,libforge),fg_err,(fg_window,),pWindow)
end

function fg_save_window_framebuffer(pFullPath,pWindow::fg_window)
    ccall((:fg_save_window_framebuffer,libforge),fg_err,(Cstring,fg_window),pFullPath,pWindow)
end


text:
@inline function logdensity(
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractVector{T},
    sptr::StackPointer = STACK_POINTER_REF[]
) where {D,T}
    @boundscheck length(θ) == D || PaddedMatrices.ThrowBoundsError()
    GC.@preserve θ begin
        θptr = PtrVector{D,T,D}(pointer(θ))
        lp = logdensity(ℓ, θptr, sptr)
    end
    lp
end
@inline function logdensity(
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractMutableFixedSizeVector{D,T},
    sptr::StackPointer = STACK_POINTER_REF[]
) where {D,T}
    GC.@preserve θ begin
        θptr = PtrVector{D,T,D}(pointer(θ))
        lp = logdensity(ℓ, θptr, sptr)
    end
    lp
end
@inline function logdensity_and_gradient!(
    ∇::AbstractVector{T},
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractVector{T},
    sptr::StackPointer = STACK_POINTER_REF[]
) where {D,T}
    @boundscheck min(length(∇),length(θ)) < D && PaddedMatrices.ThrowBoundsError()
    GC.@preserve ∇ θ begin
        ∇ptr = PtrVector{D,T,D}(pointer(∇))
        θptr = PtrVector{D,T,D}(pointer(θ))
        lp = logdensity_and_gradient!(∇ptr, ℓ, θptr, sptr)
    end
    lp
end
@inline function logdensity_and_gradient!(
    ∇::AbstractMutableFixedSizeVector{D,T},
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractMutableFixedSizeVector{D,T},
    sptr::StackPointer = STACK_POINTER_REF[]
) where {D,T}
    GC.@preserve ∇ θ begin
        ∇ptr = PtrVector{D,T,D}(pointer(∇));
        θptr = PtrVector{D,T,D}(pointer(θ));
        lp = logdensity_and_gradient!(∇ptr, ℓ, θptr, sptr)
    end
    lp
end
@inline function logdensity_and_gradient(
    sp::StackPointer,
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractMutableFixedSizeVector{D,T}
) where {D,T}
    ∇ = PtrVector{D,T,D}(pointer(sp,T))
    sp += VectorizationBase.align(D*sizeof(T))
    sp, (logdensity_and_gradient!(∇, ℓ, θ, sp), ∇)
end
@inline function logdensity_and_gradient(
    sp::StackPointer,
    ℓ::AbstractProbabilityModel{D},
    θ::AbstractVector{T}
) where {D,T}
    @boundscheck length(θ) == D || PaddedMatrices.ThrowBoundsError()
    ∇ = PtrVector{D,T,D}(pointer(sp,T))
    sp += VectorizationBase.align(D*sizeof(T))
    sp, (logdensity_and_gradient!(∇, ℓ, θ, sp), ∇)
end
@inline function logdensity_and_gradient(
    l::AbstractProbabilityModel{D}, θ, sptr::StackPointer = STACK_POINTER_REF[]
) where {D}
    ∇ = PaddedMatrices.mutable_similar(θ)
    logdensity_and_gradient!(∇, l, θ, sptr), ∇
end


text:
chars = Set{Char}()
for line in eachline(STDIN)
    m = match(r"^\d+\t(.+?)\t", line)
    if m == nothing; continue; end
    w = m.captures[1]
    for c in w
        push!(chars, c)
    end
end
println(length(chars))