text:
streams = Dict{UInt64, kStreamInfo}()

function multistream_listener(svalues::Dict{Int, Tuple{Bool, Any}}, handle::Int, handler::Function, indices)
    return (value) -> begin
        svalues[handle] = (true, value)
        if reduce((x,y) -> x && y, true, map(tup -> tup[1], values(svalues)))
            handler(map(tup -> tup[2], map(i->svalues[i], indices))...)
            for kv in svalues
                k,v = kv
                svalues[k] = (false, v[2])
            end
        end
    end
end
function add_multiple_streams(callsof::Array{kPC,1}, handler::Function)
    vls = Dict{Int, Tuple{Bool,Any}}([(i, (false, Void())) for i in 1:length(callsof)])
    for i=1:length(callsof)
        add_stream_raw(callsof[i], multistream_listener(vls, i, handler, 1:length(callsof)))
    end
end

function add_stream_raw(tocall::kPC, handler::Function; start::Bool=true)
    kstream = Remote.KRPC.AddStream(tocall.pcall)
    streams[kstream.id] = kStreamInfo(handler, Nullable{String}(), tocall.rty, Void())
    return krpc.schema.Stream(id=kstream.id)
end

function add_stream(tocall::kPC; start::Bool=true)
    kstream = Remote.KRPC.AddStream(tocall.pcall)
    init_value = SendMessage(Remote.conn, tocall)
    outp = kStreamValue(init_value)
    streams[kstream.id] = kStreamInfo(function (rv) outp.value = rv end, Nullable{String}(), tocall.rty, Void())
    return outp
end

function add_stream(tocall::Array{kPC,1}; start::Bool=true)
    return map(pc -> add_stream(pc, start=start), tocall)
end

function clear_streams()
    for strs in streams
        Remote.KRPC.RemoveStream(strs[1])
    end
    global streams = Dict{UInt64, kStreamInfo}()
end

function stream_listener(str_conn::TCPSocket, running::Array{Bool,1})
    while running[1]
        data = readproto(RecvRawProto(str_conn), krpc.schema.StreamUpdate())
        for result in data.results
            if haskey(streams, result.id)
                listenerinfo = streams[result.id]
                ires = result.result
                if ProtoBuf.has_field(ires, :error)
                    listenerinfo.error = Nullable(kRPC.make_error(ires.error))
                    continue
                end
                ivalue = getJuliaValue(ires.value, listenerinfo.typ)
                listenerinfo.value = ivalue
                try
                    Base.invokelatest(listenerinfo.listener, ivalue)
                catch err
                    clear_streams()
                    rethrow(err)
                end
            end
        end
    end
end


text:
using ClimateMachine.Mesh.BrickMesh
using Test
using MPI

MPI.Initialized() || MPI.Init()

@testset "Linear Parition" begin
    @test BrickMesh.linearpartition(1, 1, 1) == 1:1
    @test BrickMesh.linearpartition(20, 1, 1) == 1:20
    @test BrickMesh.linearpartition(10, 1, 2) == 1:5
    @test BrickMesh.linearpartition(10, 2, 2) == 6:10
end

@testset "Hilbert Code" begin
    @test BrickMesh.hilbertcode([0, 0], bits = 1) == [0, 0]
    @test BrickMesh.hilbertcode([0, 1], bits = 1) == [0, 1]
    @test BrickMesh.hilbertcode([1, 1], bits = 1) == [1, 0]
    @test BrickMesh.hilbertcode([1, 0], bits = 1) == [1, 1]
    @test BrickMesh.hilbertcode([0, 0], bits = 2) == [0, 0]
    @test BrickMesh.hilbertcode([1, 0], bits = 2) == [0, 1]
    @test BrickMesh.hilbertcode([1, 1], bits = 2) == [0, 2]
    @test BrickMesh.hilbertcode([0, 1], bits = 2) == [0, 3]
    @test BrickMesh.hilbertcode([0, 2], bits = 2) == [1, 0]
    @test BrickMesh.hilbertcode([0, 3], bits = 2) == [1, 1]
    @test BrickMesh.hilbertcode([1, 3], bits = 2) == [1, 2]
    @test BrickMesh.hilbertcode([1, 2], bits = 2) == [1, 3]
    @test BrickMesh.hilbertcode([2, 2], bits = 2) == [2, 0]
    @test BrickMesh.hilbertcode([2, 3], bits = 2) == [2, 1]
    @test BrickMesh.hilbertcode([3, 3], bits = 2) == [2, 2]
    @test BrickMesh.hilbertcode([3, 2], bits = 2) == [2, 3]
    @test BrickMesh.hilbertcode([3, 1], bits = 2) == [3, 0]
    @test BrickMesh.hilbertcode([2, 1], bits = 2) == [3, 1]
    @test BrickMesh.hilbertcode([2, 0], bits = 2) == [3, 2]
    @test BrickMesh.hilbertcode([3, 0], bits = 2) == [3, 3]

    @test BrickMesh.hilbertcode(UInt64.([14, 3, 4])) ==
          UInt64.([0x0, 0x0, 0xe25])
end

@testset "Mesh to Hilbert Code" begin
    let
        etc = Array{Float64}(undef, 2, 4, 6)
        etc[:, :, 1] = [2.0 3.0 2.0 3.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 2] = [3.0 4.0 3.0 4.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 3] = [4.0 5.0 4.0 5.0; 4.0 4.0 5.0 5.0]
        etc[:, :, 4] = [2.0 3.0 2.0 3.0; 5.0 5.0 6.0 6.0]
        etc[:, :, 5] = [3.0 4.0 3.0 4.0; 5.0 5.0 6.0 6.0]
        etc[:, :, 6] = [4.0 5.0 4.0 5.0; 5.0 5.0 6.0 6.0]

        code_exect = UInt64[
            0x0000000000000000 0x1555555555555555 0xffffffffffffffff 0x5555555555555555 0x6aaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa
            0x0000000000000000 0x5555555555555555 0xffffffffffffffff 0x5555555555555555 0xaaaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa
        ]

        code = centroidtocode(MPI.COMM_SELF, etc)

        @test code == code_exect
    end

    let
        nelem = 1
        d = 2

        etc = Array{Float64}(undef, d, d^2, nelem)
        etc[:, :, 1] = [2.0 3.0 2.0 3.0; 4.0 4.0 5.0 5.0]
        code = centroidtocode(MPI.COMM_SELF, etc)

        @test code == zeros(eltype(code), d, nelem)
    end
end

@testset "Vertex Ordering" begin
    @test ((1,), 1) == BrickMesh.vertsortandorder(1)

    @test ((1, 2), 1) == BrickMesh.vertsortandorder(1, 2)
    @test ((1, 2), 2) == BrickMesh.vertsortandorder(2, 1)

    @test ((1, 2, 3), 1) == BrickMesh.vertsortandorder(1, 2, 3)
    @test ((1, 2, 3), 2) == BrickMesh.vertsortandorder(3, 1, 2)
    @test ((1, 2, 3), 3) == BrickMesh.vertsortandorder(2, 3, 1)
    @test ((1, 2, 3), 4) == BrickMesh.vertsortandorder(2, 1, 3)
    @test ((1, 2, 3), 5) == BrickMesh.vertsortandorder(3, 2, 1)
    @test ((1, 2, 3), 6) == BrickMesh.vertsortandorder(1, 3, 2)

    @test_throws ErrorException BrickMesh.vertsortandorder(2, 1, 1)

    @test ((1, 2, 3, 4), 1) == BrickMesh.vertsortandorder(1, 2, 3, 4)
    @test ((1, 2, 3, 4), 2) == BrickMesh.vertsortandorder(1, 3, 2, 4)
    @test ((1, 2, 3, 4), 3) == BrickMesh.vertsortandorder(2, 1, 3, 4)
    @test ((1, 2, 3, 4), 4) == BrickMesh.vertsortandorder(2, 4, 1, 3)
    @test ((1, 2, 3, 4), 5) == BrickMesh.vertsortandorder(3, 1, 4, 2)
    @test ((1, 2, 3, 4), 6) == BrickMesh.vertsortandorder(3, 4, 1, 2)
    @test ((1, 2, 3, 4), 7) == BrickMesh.vertsortandorder(4, 2, 3, 1)
    @test ((1, 2, 3, 4), 8) == BrickMesh.vertsortandorder(4, 3, 2, 1)

    @test_throws ErrorException BrickMesh.vertsortandorder(1, 3, 3, 1)
end

@testset "Mesh" begin
    let
        (etv, etc, etb, fc) = brickmesh((4:7,), (false,))
        etv_expect = [
            1 2 3
            2 3 4
        ]
        etb_expect = [
            1 0 0
            0 0 1
        ]
        fc_expect = Array{Int64, 1}[]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [4 5]
        @test etc[:, :, 2] == [5 6]
        @test etc[:, :, 3] == [6 7]
    end

    let
        (etv, etc, etb, fc) = brickmesh((4:7,), (true,))
        etv_expect = [
            1 2 3
            2 3 4
        ]
        etb_expect = [
            0 0 0
            0 0 0
        ]
        fc_expect = Array{Int64, 1}[[3, 2, 1]]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [4 5]
        @test etc[:, :, 2] == [5 6]
        @test etc[:, :, 3] == [6 7]
    end

    let
        (etv, etc, etb, fc) = brickmesh((2:5, 4:6), (false, true))

        etv_expect = [
            1 2 5 6
            2 3 6 7
            3 4 7 8
            5 6 9 10
            6 7 10 11
            7 8 11 12
        ]'
        etb_expect = [
            1 0 0 1 0 0
            0 0 1 0 0 1
            0 0 0 0 0 0
            0 0 0 0 0 0
        ]
        fc_expect = Array{Int64, 1}[[4, 4, 1, 2], [5, 4, 2, 3], [6, 4, 3, 4]]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect
        @test etc[:, :, 1] == [
            2 3 2 3
            4 4 5 5
        ]
        @test etc[:, :, 5] == [
            3 4 3 4
            5 5 6 6
        ]
    end

    let
        (etv, etc, etb, fc) =
            brickmesh((-1:2:1, -1:2:1, -1:1:1), (true, true, true))
        etv_expect = [
            1 5
            2 6
            3 7
            4 8
            5 9
            6 10
            7 11
            8 12
        ]
        etb_expect = zeros(Int64, 6, 2)

        fc_expect = Array{Int64, 1}[
            [1, 2, 1, 3, 5, 7],
            [1, 4, 1, 2, 5, 6],
            [2, 2, 5, 7, 9, 11],
            [2, 4, 5, 6, 9, 10],
            [2, 6, 1, 2, 3, 4],
        ]

        @test etv == etv_expect
        @test etb == etb_expect
        @test fc == fc_expect

        @test etc[:, :, 1] == [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            -1 -1 -1 -1 0 0 0 0
        ]

        @test etc[:, :, 2] == [
            -1 1 -1 1 -1 1 -1 1
            -1 -1 1 1 -1 -1 1 1
            0 0 0 0 1 1 1 1
        ]
    end

    let
        (etv, etc, etb, fc) = brickmesh(
            (-1:1, -1:1, -1:1),
            (false, false, false),
            boundary = ((11, 12), (13, 14), (15, 16)),
        )

        @test etb == [
            11 0 11 0 11 0 11 0
            0 12 0 12 0 12 0 12
            13 13 0 0 13 13 0 0
            0 0 14 14 0 0 14 14
            15 15 15 15 0 0 0 0
            0 0 0 0 16 16 16 16
        ]
    end

    let
        x = (1:1000,)
        p = (false,)
        b = ((1, 2),)

        (etv, etc, etb, fc) = brickmesh(x, p, boundary = b)

        n = 50
        (etv_parts, etc_parts, etb_parts, fc_parts) =
            brickmesh(x, p, boundary = b, part = 1, numparts = n)
        for j in 2:n
            (etv_j, etc_j, etb_j, fc_j) =
                brickmesh(x, p, boundary = b, part = j, numparts = n)
            etv_parts = cat(etv_parts, etv_j; dims = 2)
            etc_parts = cat(etc_parts, etc_j; dims = 3)
            etb_parts = cat(etb_parts, etb_j; dims = 2)
        end

        @test etv == etv_parts
        @test etc == etc_parts
        @test etb == etb_parts
    end

    let
        x = (-1:2:10, -1:1:1, -4:1:1)
        p = (true, false, true)
        b = ((1, 2), (3, 4), (5, 6))

        (etv, etc, etb, fc) = brickmesh(x, p, boundary = b)

        n = 50
        (etv_parts, etc_parts, etb_parts, fc_parts) =
            brickmesh(x, p, boundary = b, part = 1, numparts = n)
        for j in 2:n
            (etv_j, etc_j, etb_j, fc_j) =
                brickmesh(x, p, boundary = b, part = j, numparts = n)
            etv_parts = cat(etv_parts, etv_j; dims = 2)
            etc_parts = cat(etc_parts, etc_j; dims = 3)
            etb_parts = cat(etb_parts, etb_j; dims = 2)
        end

        @test etv == etv_parts
        @test etc == etc_parts
        @test etb == etb_parts
    end
end

@testset "Connect" begin
    let
        comm = MPI.COMM_SELF

        mesh = connectmesh(
            comm,
            partition(comm, brickmesh((0:10,), (true,))...)[1:4]...,
        )

        nelem = 10

        @test mesh[:elemtocoord][:, :, 1] == [0 1]
        @test mesh[:elemtocoord][:, :, 2] == [1 2]
        @test mesh[:elemtocoord][:, :, 3] == [2 3]
        @test mesh[:elemtocoord][:, :, 4] == [3 4]
        @test mesh[:elemtocoord][:, :, 5] == [4 5]
        @test mesh[:elemtocoord][:, :, 6] == [5 6]
        @test mesh[:elemtocoord][:, :, 7] == [6 7]
        @test mesh[:elemtocoord][:, :, 8] == [7 8]
        @test mesh[:elemtocoord][:, :, 9] == [8 9]
        @test mesh[:elemtocoord][:, :, 10] == [9 10]

        @test mesh[:elemtoelem] == [
            10 1 2 3 4 5 6 7 8 9
            2 3 4 5 6 7 8 9 10 1
        ]

        @test mesh[:elemtoface] == [
            2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1
        ]

        @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))
        @test mesh[:elemtobndy] == zeros(Int, size(mesh[:elemtoordr]))

        @test mesh[:elems] == 1:nelem
        @test mesh[:realelems] == 1:nelem
        @test mesh[:ghostelems] == nelem .+ (1:0)

        @test length(mesh[:sendelems]) == 0

        @test mesh[:nabrtorank] == Int[]
        @test mesh[:nabrtorecv] == UnitRange{Int}[]
        @test mesh[:nabrtosend] == UnitRange{Int}[]
    end

    let
        comm = MPI.COMM_SELF
        mesh = connectmesh(
            comm,
            partition(comm, brickmesh((0:4, 5:9), (false, true))...)[1:4]...,
        )

        nelem = 16

        @test mesh[:elemtocoord][:, :, 1] == [0 1 0 1; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 2] == [1 2 1 2; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 3] == [1 2 1 2; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 4] == [0 1 0 1; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 5] == [0 1 0 1; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 6] == [0 1 0 1; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 7] == [1 2 1 2; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 8] == [1 2 1 2; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 9] == [2 3 2 3; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 10] == [2 3 2 3; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 11] == [3 4 3 4; 8 8 9 9]
        @test mesh[:elemtocoord][:, :, 12] == [3 4 3 4; 7 7 8 8]
        @test mesh[:elemtocoord][:, :, 13] == [3 4 3 4; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 14] == [2 3 2 3; 6 6 7 7]
        @test mesh[:elemtocoord][:, :, 15] == [2 3 2 3; 5 5 6 6]
        @test mesh[:elemtocoord][:, :, 16] == [3 4 3 4; 5 5 6 6]

        @test mesh[:elemtoelem] == [
            1 1 4 4 5 6 6 5 8 7 10 9 14 3 2 15
            2 15 14 3 8 7 10 9 12 11 11 12 13 13 16 16
            6 7 2 1 4 5 8 3 14 9 12 13 16 15 10 11
            4 3 8 5 6 1 2 7 10 15 16 11 12 9 14 13
        ]

        @test mesh[:elemtoface] == [
            1 2 2 1 1 1 2 2 2 2 2 2 2 2 2 2
            1 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2
            4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
            3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
        ]

        @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))

        @test mesh[:elems] == 1:nelem
        @test mesh[:realelems] == 1:nelem
        @test mesh[:ghostelems] == nelem .+ (1:0)

        @test length(mesh[:sendelems]) == 0

        @test mesh[:nabrtorank] == Int[]
        @test mesh[:nabrtorecv] == UnitRange{Int}[]
        @test mesh[:nabrtosend] == UnitRange{Int}[]
    end
end

@testset "Mappings" begin
    let
        comm = MPI.COMM_SELF
        x = (0:4,)
        mesh =
            connectmesh(comm, partition(comm, brickmesh(x, (true,))...)[1:4]...)

        N = 3
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        @test vmap⁻ == reshape([1, 4, 5, 8, 9, 12, 13, 16], Nfp, nface, nelem)
        @test vmap⁺ == reshape([16, 5, 4, 9, 8, 13, 12, 1], Nfp, nface, nelem)
    end

    let
        comm = MPI.COMM_SELF
        x = (-1:1, 0:1)
        p = (false, true)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = 2
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        @test vmap⁻ == reshape(
            [
                1,
                4,
                7,  # f=1 e=1
                3,
                6,
                9,  # f=2 e=1
                1,
                2,
                3,  # f=3 e=1
                7,
                8,
                9,  # f=4 e=1
                10,
                13,
                16,  # f=1 e=2
                12,
                15,
                18,  # f=2 e=2
                10,
                11,
                12,  # f=3 e=2
                16,
                17,
                18,
            ], # f=4 e=2
            Nfp,
            nface,
            nelem,
        )

        @test vmap⁺ == reshape(
            [
                1,
                4,
                7,  # f=1 e=1
                10,
                13,
                16,  # f=1 e=2
                7,
                8,
                9,  # f=4 e=1
                1,
                2,
                3,  # f=3 e=1
                3,
                6,
                9,  # f=2 e=1
                12,
                15,
                18,  # f=2 e=2
                16,
                17,
                18,  # f=4 e=2
                10,
                11,
                12,
            ], # f=3 e=2
            Nfp,
            nface,
            nelem,
        )
    end

    let
        comm = MPI.COMM_SELF
        x = (0:1, 0:1, -1:1)
        p = (false, true, false)
        mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

        N = 2
        d = length(x)
        nelem = prod(length.(x) .- 1)
        nface = 2d
        Np = (N + 1)^d
        Nfp = (N + 1)^(d - 1)

        vmap⁻, vmap⁺ =
            mappings(N, mesh[:elemtoelem], mesh[:elemtoface], mesh[:elemtoordr])

        fmask = [
            1 3 1 7 1 19
            4 6 2 8 2 20
            7 9 3 9 3 21
            10 12 10 16 4 22
            13 15 11 17 5 23
            16 18 12 18 6 24
            19 21 19 25 7 25
            22 24 20 26 8 26
            25 27 21 27 9 27
        ]

        @test vmap⁻ == reshape([fmask[:]; fmask[:] .+ Np], Nfp, nface, nelem)

        @test vmap⁺ == reshape(
            [
                fmask[:, 1]
                fmask[:, 2]
                fmask[:, 4]
                fmask[:, 3]
                fmask[:, 5]
                fmask[:, 5] .+ Np
                fmask[:, 1] .+ Np
                fmask[:, 2] .+ Np
                fmask[:, 4] .+ Np
                fmask[:, 3] .+ Np
                fmask[:, 6]
                fmask[:, 6] .+ Np
            ],
            Nfp,
            nface,
            nelem,
        )
    end
end

@testset "Get Partition" begin
    let
        Nelem = 150
        (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
        @test so == Nelem:-1:1
        @test ss == [1, Nelem + 1]
        @test rs == [1, Nelem + 1]
    end

    let
        Nelem = 111
        code = [ones(1, Nelem); collect(Nelem:-1:1)']
        (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
        @test so == Nelem:-1:1
        @test ss == [1, Nelem + 1]
        @test rs == [1, Nelem + 1]
    end
end

@testset "Partition" begin
    (etv, etc, etb, fc) =
        brickmesh((-1:2:1, -1:2:1, -2:1:2), (true, true, true))
    (netv, netc, netb, nfc) = partition(MPI.COMM_SELF, etv, etc, etb, fc)[1:4]
    @test etv == netv
    @test etc == netc
    @test etb == netb
    @test fc == nfc
end

@testset "Comm Mappings" begin
    let
        N = 1
        d = 2
        nface = 2d

        commelems = [1, 2, 5]
        commfaces = BitArray(undef, nface, length(commelems))
        commfaces .= false
        nabrtocomm = [1:2, 3:3]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        @test vmapC == Int[]
        @test nabrtovmapC == UnitRange{Int64}[1:0, 1:0]
    end

    let
        N = 1
        d = 2
        nface = 2d

        commelems = [1, 2, 5]
        commfaces = BitArray([
            false false false
            false true false
            false true false
            false false true
        ])
        nabrtocomm = [1:2, 3:3]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        @test vmapC == [5, 6, 8, 19, 20]
        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:5]
    end

    let
        N = 2
        d = 2
        nface = 2d

        commelems = [2, 4, 5]
        commfaces = BitArray([
            true true false
            false false false
            false true true
            false false true
        ])
        nabrtocomm = [1:1, 2:3]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        @test vmapC == [10, 13, 16, 28, 29, 30, 31, 34, 37, 38, 39, 43, 44, 45]
        @test nabrtovmapC == UnitRange{Int64}[1:3, 4:14]
    end

    let
        N = 2
        d = 3
        nface = 2d

        commelems = [3, 4, 7, 9]
        commfaces = BitArray([
            true true true false
            false true false false
            false true false false
            false true false false
            false true true true
            false true false true
        ])
        nabrtocomm = [1:1, 2:4]

        vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

        @test vmapC == [
            55,
            58,
            61,
            64,
            67,
            70,
            73,
            76,
            79,
            82,
            83,
            84,
            85,
            86,
            87,
            88,
            89,
            90,
            91,
            92,
            93,
            94,
            96,
            97,
            98,
            99,
            100,
            101,
            102,
            103,
            104,
            105,
            106,
            107,
            108,
            163,
            164,
            165,
            166,
            167,
            168,
            169,
            170,
            171,
            172,
            175,
            178,
            181,
            184,
            187,
            217,
            218,
            219,
            220,
            221,
            222,
            223,
            224,
            225,
            235,
            236,
            237,
            238,
            239,
            240,
            241,
            242,
            243,
        ]
        @test nabrtovmapC == UnitRange{Int64}[1:9, 10:68]
    end

end

ext: jl


text:
module PySerial

export list_ports, Serial, open, close, read, write,
       inWaiting, flush, flushInput, flushOutput

using PyCall

# import Base, so that existing functions are overloaded correctly
import Base: open, close, read, write, flush

@pyimport serial
@pyimport serial.tools.list_ports as lp

list_ports() = lp.comports()

Serial(args...; kwargs...) = serial.Serial(args...; kwargs...)

open(s::PyObject) = s[:open]()
close(s::PyObject) = s[:close]()

read(s::PyObject, size::Int = 1) = s[:read](size)
write(s::PyObject, data) = s[:write](data)

inWaiting(s::PyObject) = s[:inWaiting]()

flush(s::PyObject) = s[:flush]()
flushInput(s::PyObject) = s[:flushInput]()
flushOutput(s::PyObject) = s[:flushOutput]()

end # module


text:
#using Plots
using FFTW

include("Domain.jl")
include("Field.jl")
include("Distributor.jl")

dom = RGDomain(512,512)

hdom = get_codomain(dom)
pdom = PowerDomain(hdom)

f = randn(hdom)
fcorr = similar(f, hdom)
kernel = field(pdom, 1 ./ (1 .+ pdom._kvec).^1.5)

fcorr = PowerDistribute(f, kernel);

corr = real(fft(fcorr))
#heatmap(real(fft(fcorr.val)))


text:
@testset "PolyTopeVertexMesh" begin
    # Simply check that constructors work as intended
    @testset "TriangleMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_TriangleMesh(T,U)
            end
        end
    end

    @testset "QuadrialteralMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_QuadrilateralMesh(T,U)
            end
        end
    end

    @testset "MixedPolygonMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_MixedPolygonMesh(T,U)
            end
        end
    end

    @testset "QuadraticTriangleMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_QuadraticTriangleMesh(T,U)
            end
        end
    end

    @testset "QuadraticQuadrialteralMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_QuadraticQuadrilateralMesh(T,U)
            end
        end
    end

    @testset "MixedQuadraticPolygonMesh" begin
        for T in Floats
            for U in [UInt16, UInt32, UInt64]
                mesh = setup_MixedQuadraticPolygonMesh(T,U)
            end
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
#=
    Set path data
=#

path = joinpath(dirname(@__FILE__), "..", "..")
push!(Base.LOAD_PATH, path)
datapath = joinpath(dirname(@__FILE__), "data")

#=
    Load required libraries
=#

using Test
using Dates
using Random
using DelimitedFiles
using SparseArrays
using JuMP

#=
    select problem types to be tested
=#
sets_to_test = Symbol[]
push!(sets_to_test, :RANDSDP)
push!(sets_to_test, :SENSORLOC)
push!(sets_to_test, :SDPLIB)
push!(sets_to_test, :MIMO)

#=
    select solvers to be tested
=#
solvers = Tuple{String, Function}[]

#=
    ProxSDP with default parameters
=#
using ProxSDP
push!(solvers, ("ProxSDP", () -> ProxSDP.Optimizer(
    log_verbose=true,
    timer_verbose=true,
    time_limit = 5 * 60.0,
    log_freq = 1_000,
    )))

#=
    ProxSDP with defautl parameters excpet for FULL RANK decomposition
=#
# using ProxSDP
# push!(solvers, ("ProxSDPfull", () -> ProxSDP.Optimizer(
#     log_verbose=true,
#     time_limit = 900.0,
#     log_freq = 1_000,
#     full_eig_decomp = true,
#     )))

#=
    First order solvers
=#

# using SCS
# push!(solvers, ("SCS", () -> SCS.Optimizer(eps = 1e-4))) # eps = 1e-4

# using COSMO
# push!(solvers, ("COSMO", () -> COSMO.Optimizer(time_limit = 900.0, max_iter = 100_000, eps_abs = 1e-4))) # eps = 1e-4

# add SDPNAL+
# using SDPNAL
# push!(solvers, ("SDPNAL", () -> SDPNAL.Optimizer())) # eps = 1e-4

#=
    Interior point solvers
=#

# using MosekTools
# push!(solvers, ("MOSEK", () -> Mosek.Optimizer(MSK_DPAR_OPTIMIZER_MAX_TIME = 900.0)))

# using CSDP
# push!(solvers, ("CSDP", () -> CSDP.Optimizer(objtol=1e-4, maxiter=100000)))

# using SDPA
# push!(solvers, ("SDPA", () -> SDPA.Optimizer()))

#=
    Funstions to write results
=#

NOW = replace("$(now())",":"=>"_")
FILE = open(joinpath(dirname(@__FILE__),"proxsdp_bench_$(NOW).log"),"w")

println(FILE, "class, prob_ref, time, obj, rank, lin_feas, sdp_feas")

function println2(FILE, solver::String, class::String, ref::String, sol)
    println(FILE, "$solver, $class, $ref, $(sol[2]), $(sol[1]), $(sol[3]), $(sol[4]), $(sol[5]), $(sol[6])")
    flush(FILE)
end

#=
    Selection of problem instances
=#

RANDSDP_TEST_SET = 1:1
SENSORLOC_TEST_SET = [
    100,
    200,
    300,
    400,
]
MIMO_TEST_SET = [
    100,
    500,
    1000,
    # 1500,
    # 2000,
    ]
GPP_TEST_SET = [
    "gpp124-1.dat-s",
    "gpp124-1.dat-s",
    "gpp124-2.dat-s",
    "gpp124-3.dat-s",
    "gpp124-4.dat-s",

    "gpp250-1.dat-s",
    "gpp250-2.dat-s",
    "gpp250-3.dat-s",
    "gpp250-4.dat-s",
    "gpp500-1.dat-s",
    "gpp500-2.dat-s",
    "gpp500-3.dat-s",
    "gpp500-4.dat-s",

    # "equalG11.dat-s",
    # "equalG51.dat-s",
]
MAXCUT_TEST_SET = [
    "mcp124-1.dat-s",
    "mcp124-1.dat-s",
    "mcp124-2.dat-s",
    "mcp124-3.dat-s",
    "mcp124-4.dat-s",

    "mcp250-1.dat-s",
    "mcp250-2.dat-s",
    "mcp250-3.dat-s",
    "mcp250-4.dat-s",
    "mcp500-1.dat-s",
    "mcp500-2.dat-s",
    "mcp500-3.dat-s",
    "mcp500-4.dat-s",

    # "maxG11.dat-s"  ,
    # "maxG51.dat-s"  ,
    # "maxG32.dat-s"  ,
    # "maxG55.dat-s"  ,
    # "maxG60.dat-s"  ,
]

#=
    Load problem testing functions
=#

include("base_randsdp.jl")
include("jump_randsdp.jl")

include("base_mimo.jl")
include("jump_mimo.jl")

include("base_sensorloc.jl")
include("jump_sensorloc.jl")

include("base_sdplib.jl")
include("jump_sdplib.jl")

_randsdp = jump_randsdp
_mimo = jump_mimo
_sensorloc = jump_sensorloc
_sdplib = jump_sdplib

#=
    Run benchmarks
=#

for optimizer in solvers
    if :RANDSDP in sets_to_test
        println("RANDSDP")
        _randsdp(optimizer[2], 0, 5, 5)
        for i in RANDSDP_TEST_SET
            @show i
            sol = _randsdp(optimizer[2], i, 5, 5)
            println2(FILE, optimizer[1], "RANDSDP", "$i", sol)
        end
    end
    if :SENSORLOC in sets_to_test
        println("SENSORLOC")
        for n in SENSORLOC_TEST_SET
            @show n
            sol = _sensorloc(optimizer[2], 0, n)
            println2(FILE, optimizer[1], "SENSORLOC", "$n", sol)
        end
    end
    if :SDPLIB in sets_to_test
        println("gpp")
        for name in GPP_TEST_SET
            println(name)
            sol = _sdplib(optimizer[2], joinpath(datapath, name))
            println2(FILE, optimizer[1], "SDPLIB_gp", name, sol)
        end
        println("max_cut")
        for name in MAXCUT_TEST_SET
            println(name)
            sol = _sdplib(optimizer[2], joinpath(datapath, name))
            println2(FILE, optimizer[1], "SDPLIB_mc", name, sol)
        end
    end
    if :MIMO in sets_to_test
        println("MIMO")
        _mimo(optimizer[2], 0, 100)
        for n in MIMO_TEST_SET
            @show n
            sol = _mimo(optimizer[2], 0, n)
            println2(FILE, optimizer[1], "MIMO", "$n", sol)
        end
    end
end

#=
    Finish benhcmark file
=#

close(FILE)


text:
function bracket(f; Tinitial = 0.0, tol = 1e-9)

    step = 0.1

    # starting value
    rightbnd = Tinitial

    # go right until function peters out
    while abs(f(rightbnd)) > tol
        rightbnd += step
    end

    # starting value
    leftbnd = Tinitial

    # go left until function peters out
    while abs(f(leftbnd)) > tol
        leftbnd -= step
    end

    # return extrema
    max(abs(leftbnd), abs(rightbnd))

end


text:
module TinyFarm
#This generates small 2d examples for datacubes, can be used for example plotting
using SurrogateCube
macro load_defaults()
  esc(quote
    Ncomp=2
    Nvar=2
    Ntime=300
    Nlat=5
    Nlon=5
    kb=0.0
    kn=0.0
    ks=0.0
    kw=0.0
    b=[ConstantBaseline(0.0), ConstantBaseline(0.0)]
    n=[WhiteNoise(1.0),WhiteNoise(1.0)]
    dt=[CubeEvent(1,[0.6],[0.6],[0.2],[0.5],[0.5],[0.5]),EmptyEvent()]
    dnoise=WhiteNoise(0.1)
    coupling=Coupling[LinearCoupling(),LinearCoupling()]
    weight=[LaplaceWeight(2),LaplaceWeight(2)]
  end)
end

macro makeCube()
  esc(:(genDataCube(b,n,dt,dnoise,coupling,weight,Ncomp,Nvar,Ntime,Nlat,Nlon,kb,kn,ks,kw)))
end

function depShift()
  @load_defaults
  kw=1.0
  #Make sure event happens in all conmponents
  dt=[CubeEvent(1,[0.6],[0.6],[0.2],[0.5],[0.5],[0.5]),CubeEvent(1,[0.6],[0.6],[0.2],[0.5],[0.5],[0.5])]
  weight=[DisturbedLaplaceWeight(2),LaplaceWeight(2)]
  @makeCube
end

function baseShift()
  @load_defaults
  ks=5.0
  @makeCube
end

function noiseshift()
  @load_defaults
  kn=5.0
  @makeCube
end

function trendjump()
  @load_defaults
  ks=5.0
  dt[1]=TrendEvent(dt[1],false)
  @makeCube
end

function trendonset()
  @load_defaults
  ks=5.0
  dt[1]=TrendEvent(OnsetEvent(0.6,0.6,0.5,0.5,0.5),true)
  @makeCube
end

function quadraticdep()
  @load_defaults
  coupling[1]=QuadraticCoupling(1.0)
  @makeCube
end

function quadraticdepshift()
  @load_defaults
  coupling[2]=QuadraticCoupling(1.0)
  kw=1.0
  dt=[CubeEvent(1,[0.6],[0.6],[0.2],[0.5],[0.5],[0.5]),CubeEvent(1,[0.6],[0.6],[0.2],[0.5],[0.5],[0.5])]
  weight=[DisturbedLaplaceWeight([0.95,sqrt(1-0.95^2)],[-sqrt(1-0.95^2),0.95]),LaplaceWeight([0.9,sqrt(1-0.9^2)])]
  weight=[DisturbedLaplaceWeight(2),LaplaceWeight(2)]
  @makeCube
end

plotVars(c)=plot(1:300,slice(c.variables,:,3,3,:))
plotEvents(c)=plot(1:300,slice(c.eventmap,:,3,3,:))
plotComps(c)=plot(1:300,slice(c.components,:,3,3,:))
function plotDep(c)
  plot(c.variables[:,:,:,1][:],c.variables[:,:,:,2][:],".")
end

end

max_stars_repo_name: CAB-LAB/DataFarm

max_issues_repo_name: esa-esdl/DataFarm

max_forks_repo_name: esa-esdl/DataFarm


text:
module TestExamplesDGMulti2D

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "dgmulti_2d")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
isdir(outdir) && rm(outdir, recursive=true)

@testset "DGMulti 2D" begin
  @trixi_testset "elixir_euler_weakform.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0013536930300254945, 0.0014315603442106193, 0.001431560344211359, 0.0047393341007602625] ./ 2.0,
      linf = [0.001514260921466004, 0.0020623991944839215, 0.002062399194485476, 0.004897700392503701]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      approximation_type = SBP(),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0074706882014934735, 0.005306220583603261, 0.005306220583613591, 0.014724842607716771] ./ 2.0,
      linf = [0.021563604940952885, 0.01359397832530762, 0.013593978324845324, 0.03270995869587523]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      element_type = Quad(),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.00031892254415307093, 0.00033637562986771894, 0.0003363756298680649, 0.0011100259064243145] ./ 2.0,
      linf = [0.001073298211445639, 0.0013568139808282087, 0.0013568139808290969, 0.0032249020004324613]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (EC) " begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.007801417730672109, 0.00708583561714128, 0.0070858356171393, 0.015217574294198809] ./ 2.0,
      linf = [0.011572828457858897, 0.013965298735070686, 0.01396529873508534, 0.04227683691807904]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.01280067571168776, 0.010607599608273302, 0.010607599608239775, 0.026408338014056548] ./ 2.0,
      linf = [0.037983023185674814, 0.05321027922533417, 0.05321027922608157, 0.13392025411844033]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (Quadrilateral elements, SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (4, 4),
      element_type = Quad(),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      approximation_type = SBP(),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0029373718090697975, 0.0030629360605489465, 0.003062936060545615, 0.0068486089344859755] ./ 2.0,
      linf = [0.01360165305316885, 0.01267402847925303, 0.012674028479251254, 0.02210545278615017]
    )
  end

  @trixi_testset "elixir_euler_bilinear.jl (Bilinear quadrilateral elements, SBP, flux differencing)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_bilinear.jl"),
      l2 = [1.0259435706215337e-5, 9.014090233720625e-6, 9.014090233223014e-6, 2.738953587401793e-5],
      linf = [7.362609083649829e-5, 6.874188055272512e-5, 6.874188052830021e-5, 0.0001912435192696904]
    )
  end

  @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, SBP, flux differencing)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
      l2 = [1.5914065936559232e-5, 1.4757575957005155e-5, 1.4757575956508324e-5, 4.511709781155535e-5],
      linf = [0.00010525416930584619, 0.00010003778091061122, 0.00010003778085621029, 0.00036426282101720275]
    )
  end

  @trixi_testset "elixir_euler_curved.jl (Quadrilateral elements, GaussSBP, flux differencing)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_curved.jl"),
      approximation_type = GaussSBP(),
      l2 = [3.578819171269236e-6, 3.4797454435720057e-6, 3.4797454434088854e-6, 1.0812715223492364e-5],
      linf = [1.2963211813321607e-5, 1.2570312152959673e-5, 1.2570312227122571e-5, 3.9389540817946767e-5]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (convergence)" begin
    mean_convergence = convergence_test(@__MODULE__, joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"), 2)
    @test isapprox(mean_convergence[:l2], [4.243843382379403, 4.128314378833922, 4.128314378397532, 4.081366752807379], rtol=0.05)
  end

  @trixi_testset "elixir_euler_weakform_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform_periodic.jl"),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0014986508075708323, 0.001528523420746786, 0.0015285234207473158, 0.004846505183839211] ./ 2.0,
      linf = [0.0015062108658376872, 0.0019373508504645365, 0.0019373508504538783, 0.004742686826709086]
    )
  end

  @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_triangulate_pkg_mesh.jl"),
      l2 = [2.344080455438114e-6, 1.8610038753097983e-6, 2.4095165666095305e-6, 6.373308158814308e-6],
      linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
    )
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      cells_per_dimension = (32, 32), tspan = (0.0, 0.2),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.11140378947116614, 0.06598161188703612, 0.10448953167839563, 0.16023209181809595] ./ 2.0,
      linf = [0.24033843177853664, 0.1659992245272325, 0.1235468309508845, 0.26911424973147735]
    )
  end

  @trixi_testset "elixir_euler_kelvin_helmholtz_instability.jl (Quadrilateral elements, GaussSBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_kelvin_helmholtz_instability.jl"),
      cells_per_dimension = (32, 32), element_type = Quad(), approximation_type=GaussSBP(), tspan = (0.0, 0.2),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.11141270656347146, 0.06598888014584121, 0.1044902203749932, 0.16023037364774995] ./ 2.0,
      linf = [0.2414760062126462, 0.1662111846065654, 0.12344140473946856, 0.26978428189564774]
    )
  end

  @trixi_testset "elixir_euler_rayleigh_taylor_instability.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_rayleigh_taylor_instability.jl"),
      cells_per_dimension = (8, 8), tspan = (0.0, 0.2),
      l2 = [0.0709665896982514, 0.005182828752164663, 0.013832655585206478, 0.03247013800580221],
      linf = [0.4783963902824797, 0.022527207050681054, 0.040307056293369226, 0.0852365428206836]
    )
  end

  @trixi_testset "elixir_euler_brown_minion_vortex.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_brown_minion_vortex.jl"),
      cells_per_dimension = 4, tspan = (0.0, 0.1),
      l2 = [0.006680001611078062, 0.02151676347585447, 0.010696524235364626, 0.15052841129694647],
      linf = [0.01544756362800248, 0.09517304772476806, 0.021957154972646383, 0.33773439650806303]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (FD SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (2, 2),
      element_type = Quad(),
      cfl = 1.0,
      approximation_type = derivative_operator(
        SummationByPartsOperators.MattssonNordström2004(),
        derivative_order=1, accuracy_order=4,
        xmin=0.0, xmax=1.0, N=12),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0008966318978421226, 0.0011418826379110242, 0.001141882637910878, 0.0030918374335671393] ./ 2.0,
      linf = [0.0015281525343109337, 0.00162430960401716, 0.0016243096040242655, 0.004447503691245913]
    )
  end

  @trixi_testset "elixir_euler_weakform.jl (FD SBP, EC)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_weakform.jl"),
      cells_per_dimension = (2, 2),
      element_type = Quad(),
      cfl = 1.0,
      approximation_type = derivative_operator(
        SummationByPartsOperators.MattssonNordström2004(),
        derivative_order=1, accuracy_order=4,
        xmin=0.0, xmax=1.0, N=12),
      volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
      surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.0014018725496871129, 0.0015887007320868913, 0.001588700732086329, 0.003870926821031202] ./ 2.0,
      linf = [0.0029541996523780867, 0.0034520465226108854, 0.003452046522624652, 0.007677153211004928]
    )
  end

  @trixi_testset "elixir_euler_fdsbp_periodic.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
      l2 = [1.3333320340010056e-6, 2.044834627970641e-6, 2.044834627855601e-6, 5.282189803559564e-6],
      linf = [2.7000151718858945e-6, 3.988595028259212e-6, 3.9885950273710336e-6, 8.848583042286862e-6]
    )
  end

  @trixi_testset "elixir_euler_fdsbp_periodic.jl (arbitrary reference domain)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
      xmin=-200.0, xmax=100.0 #= parameters for reference interval =#,
      l2 = [1.333332034149886e-6, 2.0448346280892024e-6, 2.0448346279766305e-6, 5.282189803510037e-6],
      linf = [2.700015170553627e-6, 3.988595024262409e-6, 3.988595024928543e-6, 8.84858303740188e-6]
    )
  end

  @trixi_testset "elixir_euler_fdsbp_periodic.jl (arbitrary reference and physical domains)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
      approximation_type = periodic_derivative_operator(
        derivative_order=1, accuracy_order=4, xmin=-200.0, xmax=100.0, N=100),
      coordinates_min=(-3.0, -4.0), coordinates_max=(0.0, -1.0),
      l2 = [0.07318831033918516, 0.10039910610067465, 0.1003991061006748, 0.2642450566234564],
      linf = [0.36081081739439735, 0.5244468027020845, 0.5244468027020814, 1.2210130256735705]
    )
  end

  @trixi_testset "elixir_euler_fdsbp_periodic.jl (CGSEM)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_fdsbp_periodic.jl"),
      approximation_type = SummationByPartsOperators.couple_continuously(
        SummationByPartsOperators.legendre_derivative_operator(xmin=0.0, xmax=1.0, N=4),
        SummationByPartsOperators.UniformPeriodicMesh1D(xmin=-1.0, xmax=1.0, Nx=10)),
      l2 = [1.5440402410017893e-5, 1.4913189903083485e-5, 1.4913189902797073e-5, 2.6104615985156992e-5],
      linf = [4.16334345412217e-5, 5.067812788173143e-5, 5.067812786885284e-5, 9.887976803746312e-5]
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave.jl (Quad)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave.jl"),
      cells_per_dimension = 4,
      l2 = [0.03906769915509508, 0.04923079758984701, 0.049230797589847136, 0.02660348840973283,
            0.18054907061740028, 0.019195256934309846, 0.019195256934310016, 0.027856113419468087,
            0.0016567799774264065],
      linf = [0.16447597822733662, 0.244157345789029, 0.24415734578903472, 0.11982440036793476,
              0.7450328339751362, 0.06357382685763713, 0.0635738268576378, 0.1058830287485999,
              0.005740591170062146]
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave.jl (Tri)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave.jl"),
      cells_per_dimension = 4, element_type = Tri(),
      l2 = [0.03372468091254386, 0.03971626483409167, 0.03971626483409208, 0.021427571421535722,
            0.15079331840847413, 0.015716300366650286, 0.015716300366652128, 0.022365252076279075,
            0.0009232971979900358],
      linf = [0.16290247390873458, 0.2256891306641319, 0.2256891306641336, 0.09476017042552534,
              0.6906308908961734, 0.05349939593012487, 0.05349939593013042, 0.08830587480616725,
              0.0029551359803035027]
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave_SBP.jl (Quad)" begin
    # These setups do not pass CI reliably, see
    # https://github.com/trixi-framework/Trixi.jl/pull/880 and
    # https://github.com/trixi-framework/Trixi.jl/issues/881
    @test_skip @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave_SBP.jl"),
      cells_per_dimension = 4,
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.15825983698241494, 0.19897219694837923, 0.19784182473275247, 0.10482833997417325,
            0.7310752391255246, 0.07374056714564853, 0.07371172293240634, 0.10782032253431281,
            0.004921676235111545] ./ 2.0,
      linf = [0.1765644464978685, 0.2627803272865769, 0.26358136695848144, 0.12347681727447984,
              0.7733289736898254, 0.06695360844467957, 0.06650382120802623, 0.10885097000919097,
              0.007212567638078835]
    )
  end

  @trixi_testset "elixir_mhd_weak_blast_wave_SBP.jl (Tri)" begin
    # These setups do not pass CI reliably, see
    # https://github.com/trixi-framework/Trixi.jl/pull/880 and
    # https://github.com/trixi-framework/Trixi.jl/issues/881
    @test_skip @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_weak_blast_wave_SBP.jl"),
      cells_per_dimension = 4, element_type=Tri(), tspan = (0.0, 0.2),
      # division by 2.0 corresponds to normalization by the square root of the size of the domain
      l2 = [0.13825044764021147, 0.15472815448314997, 0.1549093274293255, 0.053103596213755405,
            0.7246162776815603, 0.07730777596615901, 0.07733438386480523, 0.109893463921706,
            0.00617678167062838] ./ 2.0,
      linf = [0.22701306227317952, 0.2905255794821543, 0.2912409425436937, 0.08051361477962096,
              1.0842974228656006, 0.07866000368926784, 0.0786646354518149, 0.1614896380292925,
              0.010358210347485542]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Quad(), approximation_type = SBP(),
      l2 = [0.0020316462913319046, 0.023669019044882247, 0.03446194752754684, 1.9333465252381796e-15],
      linf = [0.010385010095182778, 0.08750628939565086, 0.12088392994348407, 9.325873406851315e-15]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Tri, SBP)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Tri(), approximation_type = SBP(),
      l2 = [0.004180680322490383, 0.07026192411558974, 0.11815151697006446, 2.329788936151192e-15],
      linf = [0.02076003852980346, 0.29169601664914424, 0.5674183379872275, 1.1546319456101628e-14]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Tri, Polynomial)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Tri(), approximation_type = Polynomial(),
      # The last l2, linf error are the L2 projection error in approximating `b`, so they are not
      # zero for general non-collocated quadrature rules (e.g., for `element_type=Tri()`, `polydeg > 2`).
      l2 = [0.0008309356912456799, 0.01522451288799231, 0.016033969387208476, 1.2820247308150876e-5],
      linf = [0.001888045014140971, 0.05466838692127718, 0.06345885709961152, 3.3989933098554914e-5]
    )
  end

  @trixi_testset "elixir_shallowwater_source_terms.jl (Quad, Polynomial)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_shallowwater_source_terms.jl"),
      cells_per_dimension = 8, element_type = Quad(), approximation_type = Polynomial(),
      # The last l2, linf error are the L2 projection error in approximating `b`. However, this is zero
      # for `Quad()` elements with `Polynomial()` approximations because the quadrature rule defaults to
      # a `(polydeg + 1)`-point Gauss quadrature rule in each coordinate (in general, StartUpDG.jl defaults
      # to the quadrature rule with the fewest number of points which exactly integrates the mass matrix).
      l2 = [7.460461950323111e-5, 0.003685589808444905, 0.0039101604749887785, 2.0636891126652983e-15],
      linf = [0.000259995400729629, 0.0072236204211630906, 0.010364675200833062, 1.021405182655144e-14]
    )
  end

end

# Clean up afterwards: delete Trixi output directory
@test_nowarn isdir(outdir) && rm(outdir, recursive=true)

end # module


text:
# To update the reference data, run
#   julia --project=. update_reference_data.jl
using ParticleDA, HDF5, StableRNGs
include(joinpath(@__DIR__, "model", "model.jl"))
using .Model

reference_file = joinpath(@__DIR__, "reference_data.h5")
rm(reference_file; force=true)

x_true, _, _ = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_1.yaml"), BootstrapFilter())
h5write(reference_file, "integration_test_1", x_true)

x_true, _, _ = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_2.yaml"), BootstrapFilter())
h5write(reference_file, "integration_test_2", x_true)

rng = StableRNG(123)
_, x_avg, _ = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "integration_test_4.yaml"), BootstrapFilter(); rng)
h5write(reference_file, "integration_test_4", x_avg)

rng = StableRNG(123)
_, x_avg, _ = ParticleDA.run_particle_filter(Model.init, joinpath(@__DIR__, "optimal_filter_test_1.yaml"), OptimalFilter(); rng)
h5write(reference_file, "integration_test_6", x_avg)


text:
using Documenter

deploydocs(;
    repo="github.com/OutlierDetectionJL/OutlierDetection.jl",
    target = "site",
    push_preview = false,
    # mkdocs-material bundles all other required dependencies
    deps = Deps.pip("mkdocs-material==7.3.6"),
    make = () -> run(`mkdocs build`)
)


text:
# Licensed to the Apache Software Foundation (ASF) under one
# or more contributor license agreements.  See the NOTICE file
# distributed with this work for additional information
# regarding copyright ownership.  The ASF licenses this file
# to you under the Apache License, Version 2.0 (the
# "License"); you may not use this file except in compliance
# with the License.  You may obtain a copy of the License at
#
# http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

FlatBuffers.@scopedenum MetadataVersion::Int16 V1 V2 V3 V4 V5

struct Null <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Null) = ()

nullStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
nullEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Struct <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Struct) = ()

structStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
structEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct List <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::List) = ()

listStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
listEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct LargeList <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::LargeList) = ()

largeListStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
largeListEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct FixedSizeList <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::FixedSizeList) = (:listSize,)

function Base.getproperty(x::FixedSizeList, field::Symbol)
    if field === :listSize
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
        return Int32(0)
    end
    return nothing
end

fixedSizeListStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
fixedSizeListAddListSize(b::FlatBuffers.Builder, listSize::Int32) = FlatBuffers.prependslot!(b, 0, listSize, 0)
fixedSizeListEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Map <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Map) = (:keysSorted,)

function Base.getproperty(x::Map, field::Symbol)
    if field === :keysSorted
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Base.Bool)
    end
    return nothing
end

mapStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
mapAddKeysSorted(b::FlatBuffers.Builder, keyssorted::Base.Bool) = FlatBuffers.prependslot!(b, 0, keyssorted, 0)
mapEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum UnionMode::Int16 Sparse Dense

struct Union <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Union) = (:mode, :typeIds)

function Base.getproperty(x::Union, field::Symbol)
    if field === :mode
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), UnionMode)
        return UnionModes.Sparse
    elseif field === :typeIds
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.Array{Int32}(x, o)
    end
    return nothing
end

unionStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
unionAddMode(b::FlatBuffers.Builder, mode::UnionMode) = FlatBuffers.prependslot!(b, 0, mode, 0)
unionAddTypeIds(b::FlatBuffers.Builder, typeIds::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, typeIds, 0)
unionStartTypeIdsVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
unionEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Int <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Int) = (:bitWidth, :is_signed)

function Base.getproperty(x::Int, field::Symbol)
    if field === :bitWidth
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
    elseif field === :is_signed
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Base.Bool)
        return false
    end
    return nothing
end

intStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
intAddBitWidth(b::FlatBuffers.Builder, bitwidth::Int32) = FlatBuffers.prependslot!(b, 0, bitwidth, 0)
intAddIsSigned(b::FlatBuffers.Builder, issigned::Base.Bool) = FlatBuffers.prependslot!(b, 1, issigned, 0)
intEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum Precision::Int16 HALF SINGLE DOUBLE

struct FloatingPoint <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::FloatingPoint) = (:precision,)

function Base.getproperty(x::FloatingPoint, field::Symbol)
    if field === :precision
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Precision)
        return Precisions.HALF
    end
    return nothing
end

floatingPointStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
floatingPointAddPrecision(b::FlatBuffers.Builder, precision::Precision) = FlatBuffers.prependslot!(b, 0, precision, 0)
floatingPointEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Utf8 <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Utf8) = ()

utf8Start(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
utf8End(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Binary <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Binary) = ()

binaryStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
binaryEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct LargeUtf8 <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::LargeUtf8) = ()

largUtf8Start(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
largUtf8End(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct LargeBinary <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::LargeBinary) = ()

largeBinaryStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
largeBinaryEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct FixedSizeBinary <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::FixedSizeBinary) = (:byteWidth,)

function Base.getproperty(x::FixedSizeBinary, field::Symbol)
    if field === :byteWidth
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
    end
    return nothing
end

fixedSizeBinaryStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
fixedSizeBinaryAddByteWidth(b::FlatBuffers.Builder, bytewidth::Int32) = FlatBuffers.prependslot!(b, 0, bytewidth, 0)
fixedSizeBinaryEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Bool <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Bool) = ()

boolStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 0)
boolEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Decimal <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Decimal) = (:precision, :scale, :bitWidth)

function Base.getproperty(x::Decimal, field::Symbol)
    if field === :precision
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
        return Int32(0)
    elseif field === :scale
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
        return Int32(0)
    elseif field === :bitWidth
        o = FlatBuffers.offset(x, 8)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
        return Int32(128)
    end
    return nothing
end

decimalStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 3)
decimalAddPrecision(b::FlatBuffers.Builder, precision::Int32) = FlatBuffers.prependslot!(b, 0, precision, 0)
decimalAddScale(b::FlatBuffers.Builder, scale::Int32) = FlatBuffers.prependslot!(b, 1, scale, 0)
decimalAddBitWidth(b::FlatBuffers.Builder, bitWidth::Int32) = FlatBuffers.prependslot!(b, 2, bitWidth, Int32(128))
decimalEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum DateUnit::Int16 DAY MILLISECOND

struct Date <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Date) = (:unit,)

function Base.getproperty(x::Date, field::Symbol)
    if field === :unit
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), DateUnit)
        return DateUnits.MILLISECOND
    end
    return nothing
end

dateStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
dateAddUnit(b::FlatBuffers.Builder, unit::DateUnit) = FlatBuffers.prependslot!(b, 0, unit, 1)
dateEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum TimeUnit::Int16 SECOND MILLISECOND MICROSECOND NANOSECOND

struct Time <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Time) = (:unit, :bitWidth)

function Base.getproperty(x::Time, field::Symbol)
    if field === :unit
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), TimeUnit)
        return TimeUnits.MILLISECOND
    elseif field === :bitWidth
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int32)
        return 32
    end
    return nothing
end

timeStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
timeAddUnit(b::FlatBuffers.Builder, unit::TimeUnit) = FlatBuffers.prependslot!(b, 0, unit, 1)
timeAddBitWidth(b::FlatBuffers.Builder, bitwidth::Int32) = FlatBuffers.prependslot!(b, 1, bitwidth, 32)
timeEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Timestamp <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Timestamp) = (:unit, :timezone)

function Base.getproperty(x::Timestamp, field::Symbol)
    if field === :unit
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), TimeUnit)
        return TimeUnits.SECOND
    elseif field === :timezone
        o = FlatBuffers.offset(x, 6)
        o != 0 && return String(x, o + FlatBuffers.pos(x))
    end
    return nothing
end

timestampStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
timestampAddUnit(b::FlatBuffers.Builder, unit::TimeUnit) = FlatBuffers.prependslot!(b, 0, unit, 0)
timestampAddTimezone(b::FlatBuffers.Builder, timezone::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, timezone, 0)
timestampEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum IntervalUnit::Int16 YEAR_MONTH DAY_TIME

struct Interval <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Interval) = (:unit,)

function Base.getproperty(x::Interval, field::Symbol)
    if field === :unit
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), IntervalUnit)
        return IntervalUnits.YEAR_MONTH
    end
    return nothing
end

intervalStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
intervalAddUnit(b::FlatBuffers.Builder, unit::IntervalUnit) = FlatBuffers.prependslot!(b, 0, unit, 0)
intervalEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Duration <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Duration) = (:unit,)

function Base.getproperty(x::Duration, field::Symbol)
    if field === :unit
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), TimeUnit)
        return TimeUnits.MILLISECOND
    end
    return nothing
end

durationStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 1)
durationAddUnit(b::FlatBuffers.Builder, unit::TimeUnit) = FlatBuffers.prependslot!(b, 0, unit, 1)
durationEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

function Type(b::UInt8)
    b == 1 && return Null
    b == 2 && return Int
    b == 3 && return FloatingPoint
    b == 4 && return Binary
    b == 5 && return Utf8
    b == 6 && return Bool
    b == 7 && return Decimal
    b == 8 && return Date
    b == 9 && return Time
    b == 10 && return Timestamp
    b == 11 && return Interval
    b == 12 && return List
    b == 13 && return Struct
    b == 14 && return Union
    b == 15 && return FixedSizeBinary
    b == 16 && return FixedSizeList
    b == 17 && return Map
    b == 18 && return Duration
    b == 19 && return LargeBinary
    b == 20 && return LargeUtf8
    b == 21 && return LargeList
    return nothing
end

function Type(::Base.Type{T})::Int16 where {T}
    T == Null && return 1
    T == Int && return 2
    T == FloatingPoint && return 3
    T == Binary && return 4
    T == Utf8 && return 5
    T == Bool && return 6
    T == Decimal && return 7
    T == Date && return 8
    T == Time && return 9
    T == Timestamp && return 10
    T == Interval && return 11
    T == List && return 12
    T == Struct && return 13
    T == Union && return 14
    T == FixedSizeBinary && return 15
    T == FixedSizeList && return 16
    T == Map && return 17
    T == Duration && return 18
    T == LargeBinary && return 19
    T == LargeUtf8 && return 20
    T == LargeList && return 21
    return 0
end

struct KeyValue <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::KeyValue) = (:key, :value)

function Base.getproperty(x::KeyValue, field::Symbol)
    if field === :key
        o = FlatBuffers.offset(x, 4)
        o != 0 && return String(x, o + FlatBuffers.pos(x))
    elseif field === :value
        o = FlatBuffers.offset(x, 6)
        o != 0 && return String(x, o + FlatBuffers.pos(x))
    end
    return nothing
end

keyValueStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 2)
keyValueAddKey(b::FlatBuffers.Builder, key::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 0, key, 0)
keyValueAddValue(b::FlatBuffers.Builder, value::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, value, 0)
keyValueEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum DictionaryKind::Int16 DenseArray

struct DictionaryEncoding <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::DictionaryEncoding) = (:id, :indexType, :isOrdered, :dictionaryKind)

function Base.getproperty(x::DictionaryEncoding, field::Symbol)
    if field === :id
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Int64)
        return Int64(0)
    elseif field === :indexType
        o = FlatBuffers.offset(x, 6)
        if o != 0
            y = FlatBuffers.indirect(x, o + FlatBuffers.pos(x))
            return FlatBuffers.init(Int, FlatBuffers.bytes(x), y)
        end
    elseif field === :isOrdered
        o = FlatBuffers.offset(x, 8)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Base.Bool)
        return false
    elseif field === :dictionaryKind
        o = FlatBuffers.offset(x, 10)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), DictionaryKind)
    end
    return nothing
end

dictionaryEncodingStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 3)
dictionaryEncodingAddId(b::FlatBuffers.Builder, id::Int64) = FlatBuffers.prependslot!(b, 0, id, 0)
dictionaryEncodingAddIndexType(b::FlatBuffers.Builder, indextype::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, indextype, 0)
dictionaryEncodingAddIsOrdered(b::FlatBuffers.Builder, isordered::Base.Bool) = FlatBuffers.prependslot!(b, 1, isordered, 0)
dictionaryEncodingEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

struct Field <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Field) = (:name, :nullable, :type, :dictionary, :children, :custom_metadata)

function Base.getproperty(x::Field, field::Symbol)
    if field === :name
        o = FlatBuffers.offset(x, 4)
        o != 0 && return String(x, o + FlatBuffers.pos(x))
    elseif field === :nullable
        o = FlatBuffers.offset(x, 6)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Base.Bool)
        return false
    elseif field === :type
        o = FlatBuffers.offset(x, 8)
        if o != 0
            T = Type(FlatBuffers.get(x, o + FlatBuffers.pos(x), UInt8))
            o = FlatBuffers.offset(x, 10)
            pos = FlatBuffers.union(x, o)
            if o != 0
                return FlatBuffers.init(T, FlatBuffers.bytes(x), pos)
            end
        end
    elseif field === :dictionary
        o = FlatBuffers.offset(x, 12)
        if o != 0
            y = FlatBuffers.indirect(x, o + FlatBuffers.pos(x))
            return FlatBuffers.init(DictionaryEncoding, FlatBuffers.bytes(x), y)
        end
    elseif field === :children
        o = FlatBuffers.offset(x, 14)
        if o != 0
            return FlatBuffers.Array{Field}(x, o)
        end
    elseif field === :custom_metadata
        o = FlatBuffers.offset(x, 16)
        if o != 0
            return FlatBuffers.Array{KeyValue}(x, o)
        end
    end
    return nothing
end

fieldStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 7)
fieldAddName(b::FlatBuffers.Builder, name::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 0, name, 0)
fieldAddNullable(b::FlatBuffers.Builder, nullable::Base.Bool) = FlatBuffers.prependslot!(b, 1, nullable, false)
fieldAddTypeType(b::FlatBuffers.Builder, ::Core.Type{T}) where {T} = FlatBuffers.prependslot!(b, 2, Type(T), 0)
fieldAddType(b::FlatBuffers.Builder, type::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 3, type, 0)
fieldAddDictionary(b::FlatBuffers.Builder, dictionary::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 4, dictionary, 0)
fieldAddChildren(b::FlatBuffers.Builder, children::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 5, children, 0)
fieldStartChildrenVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
fieldAddCustomMetadata(b::FlatBuffers.Builder, custommetadata::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 6, custommetadata, 0)
fieldStartCustomMetadataVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
fieldEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

FlatBuffers.@scopedenum Endianness::Int16 Little Big

struct Buffer <: FlatBuffers.Struct
    bytes::Vector{UInt8}
    pos::Base.Int
end

FlatBuffers.structsizeof(::Base.Type{Buffer}) = 16

Base.propertynames(x::Buffer) = (:offset, :length)

function Base.getproperty(x::Buffer, field::Symbol)
    if field === :offset
        return FlatBuffers.get(x, FlatBuffers.pos(x), Int64)
    elseif field === :length
        return FlatBuffers.get(x, FlatBuffers.pos(x) + 8, Int64)
    end
    return nothing
end

function createBuffer(b::FlatBuffers.Builder, offset::Int64, length::Int64)
    FlatBuffers.prep!(b, 8, 16)
    prepend!(b, length)
    prepend!(b, offset)
    return FlatBuffers.offset(b)
end

struct Schema <: FlatBuffers.Table
    bytes::Vector{UInt8}
    pos::Base.Int
end

Base.propertynames(x::Schema) = (:endianness, :fields, :custom_metadata)

function Base.getproperty(x::Schema, field::Symbol)
    if field === :endianness
        o = FlatBuffers.offset(x, 4)
        o != 0 && return FlatBuffers.get(x, o + FlatBuffers.pos(x), Endianness)
    elseif field === :fields
        o = FlatBuffers.offset(x, 6)
        if o != 0
            return FlatBuffers.Array{Field}(x, o)
        end
    elseif field === :custom_metadata
        o = FlatBuffers.offset(x, 8)
        if o != 0
            return FlatBuffers.Array{KeyValue}(x, o)
        end
    end
    return nothing
end

schemaStart(b::FlatBuffers.Builder) = FlatBuffers.startobject!(b, 3)
schemaAddEndianness(b::FlatBuffers.Builder, endianness::Endianness) = FlatBuffers.prependslot!(b, 0, endianness, 0)
schemaAddFields(b::FlatBuffers.Builder, fields::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 1, fields, 0)
schemaStartFieldsVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
schemaAddCustomMetadata(b::FlatBuffers.Builder, custommetadata::FlatBuffers.UOffsetT) = FlatBuffers.prependoffsetslot!(b, 2, custommetadata, 0)
schemaStartCustomMetadataVector(b::FlatBuffers.Builder, numelems) = FlatBuffers.startvector!(b, 4, numelems, 4)
schemaEnd(b::FlatBuffers.Builder) = FlatBuffers.endobject!(b)

max_stars_repo_name: jonalm/arrow-julia

max_issues_repo_name: jonalm/arrow-julia

max_forks_repo_name: jonalm/arrow-julia


text:
# table
# column
# section

function env_column(e, _)
    md = F.content(e)
    _, kwargs = lxargs(lxproc(e), "column")
    return _column(md; kwargs...)
end

function _column(md; opts="")
    return html("
") * md * html("
")
end

function env_table(e, _)
    md = F.content(e)
    _, kwargs = lxargs(lxproc(e), "table")
    return _table(md; kwargs...)
end

function _table(md; style="", caption="", class="")
    caption = ifelse(isempty(caption), "", "
$caption
")
    parser = CommonMark.Parser()
    CommonMark.enable!(parser, CommonMark.TableRule())
    h = CommonMark.html(parser(md))
    h = replace(h, "
"=>"")
    return html("""
$caption $h""")
end

function env_section(e, _)
    md = F.content(e)
    _, kwargs = lxargs(lxproc(e), "section")
    return _section(md; kwargs...)
end

function _section(md; title="", name=title, width=F.globvar("section_width"))
    width = ifelse(isnothing(width), 10, width)
    id = F.refstring(name)
    pair = (id => name)
    if pair in F.locvar("sections")
        empty!(F.LOCAL_VARS["sections"].first)
    end
    push!(F.LOCAL_VARS["sections"].first, pair)
    counter = F.globvar("section_counter")
    F.set_var!(F.GLOBAL_VARS, "section_counter", counter+1)
    class = "scrollspy"
    if iseven(counter)
        class *= " section-bg-color"
    end
    return html("""

$title
""") * md * html("""

""")
end


text:
function (+)(a::RC, b::RC)
    if isnan(a) || isnan(b)
        return ComplexNaN
    end
    if isinf(a) && isinf(b)
        return ComplexNaN
    end
    if isinf(a) || isinf(b)
        return ComplexInf
    end
    return RC(a.val + b.val)
end

function (-)(a::RC)
    return RC(-a.val, a.nan_flag, a.inf_flag)
end

function (-)(a::RC,b::RC)
    return a + (-b)
end

function (*)(a::RC, b::RC)
    if isnan(a) || isnan(b)
        return ComplexNaN
    end
    if (isinf(a)&&iszero(b))||(iszero(a)&&isinf(b))  # 0 x Inf
        return ComplexNaN
    end
    if isinf(a) || isinf(b)
        return ComplexInf
    end
    return RC(a.val * b.val)
end

function inv(z::RC{T}) where T
    if isnan(z)
        return ComplexNaN
    end
    if iszero(z)
        return ComplexInf
    end
    if isinf(z) return
        RC(zero(T))
    end
    return RC(1/z.val)
end

function (/)(a::RC{S}, b::RC{T}) where {S,T}
    if isnan(a) || isnan(b)
        return ComplexNaN
    end

    # zero denominator cases
    if iszero(b) && !iszero(a)
        return ComplexInf
    end
    if iszero(b) && iszero(a)
        return ComplexNaN
    end

    # infinite denominator cases
    if isinf(a) && isinf(b)
        return ComplexNaN
    end
    if isinf(b)
        return RC(zero(S))
    end

    # infinite numerator, nonzero/noninf denominator
    if isinf(a)
        return a
    end

    return RC(a.val / b.val)
end


text:
#TODO: seed, copy

mutable struct MultiBinary <: AbstractSpace
    n::Int
    dtype::DataType
    shape::Tuple
    MultiBinary(n::Int) = new(n, BitArray{1}, (n, ))
end

sample(multibin_obj::MultiBinary) = (multibin_obj.dtype)(rand(0:1, multibin_obj.n))

contains(multibin_obj::MultiBinary, x) = all((x .== 0) .| (x .== 1))

Base.:(==)(multibin_obj::MultiBinary, other::MultiBinary) = multibin_obj.n == other.n


text:
using Test, AdvancedHMC
using AdvancedHMC: DualValue, PhasePoint
using LinearAlgebra: dot, diagm

include("common.jl")

@testset "PhasePoint" begin
    init_z1() = PhasePoint([NaN], [NaN], DualValue(0.,[0.]), DualValue(0.,[0.]))
    init_z2() = PhasePoint([Inf], [Inf], DualValue(0.,[0.]), DualValue(0.,[0.]))

    @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z1()
    @test_logs (:warn, "The current proposal will be rejected due to numerical error(s).") init_z2()

    z1 = init_z1()
    z2 = init_z2()

    @test z1.ℓπ.value == z2.ℓπ.value
    @test z1.ℓκ.value == z2.ℓκ.value
end

@testset "Metric" begin
    n_tests = 10

    for _ in 1:n_tests
        θ_init = randn(D)
        r_init = randn(D)

        h = Hamiltonian(UnitEuclideanMetric(D), ℓπ, ∂ℓπ∂θ)
        @test -AdvancedHMC.neg_energy(h, r_init, θ_init) == sum(abs2, r_init) / 2
        @test AdvancedHMC.∂H∂r(h, r_init) == r_init

        M⁻¹ = ones(D) + abs.(randn(D))
        h = Hamiltonian(DiagEuclideanMetric(M⁻¹), ℓπ, ∂ℓπ∂θ)
        @test -AdvancedHMC.neg_energy(h, r_init, θ_init) ≈ r_init' * diagm(0 => M⁻¹) * r_init / 2
        @test AdvancedHMC.∂H∂r(h, r_init) == M⁻¹ .* r_init

        m = randn(D, D)
        M⁻¹ = m' * m
        h = Hamiltonian(DenseEuclideanMetric(M⁻¹), ℓπ, ∂ℓπ∂θ)
        @test -AdvancedHMC.neg_energy(h, r_init, θ_init) ≈ r_init' * M⁻¹ * r_init / 2
        @test AdvancedHMC.∂H∂r(h, r_init) == M⁻¹ * r_init
    end
end


text:
VERSION < v"0.7.0-beta2.199" && __precompile__()

module Tau
export tau, τ,
       sintau, costau, modtau,
       sinτ, cosτ, modτ

# Use overridden macro definition to define conversion methods for tau
Base.@irrational τ 6.28318530717958647692 (2 * big(pi))
const tau = τ

include("trig.jl")

modtau(x) = Base.mod2pi(x)
const modτ = modtau

end


text:
using Randomorg_API
using Test
@testset "Randomorg_API.jl" begin
    generateIntegerSequences(2, 6, 5, 10,false)
    generateIntegers(5, 1, 10, false)
    generateDecimalFractions(5, 6,false)
    generateGaussians(5, 10, 9.1, 3)
    generateStrings(5, 10, "abcdefghi", false)
    generateUUIDs(5)
    generateBlobs(5, 1000, "hex")
end


text:
@testset "random" begin
    @test_skip run(sess, categorical(10))
    for shape in [[10], [10,5]]
        @test_skip run(sess, uniform(shape...))
        @test_skip run(sess, normal(shape...))
        @test_skip run(sess, choice(rand(10,2), 4))
        @test_skip run(sess, choice(rand(10,2), 4, replace=true))
        @test_skip run(sess, beta(shape...;concentration1=1.0, concentration0=1.0))
        @test_skip run(sess, bernoulli(shape...; probs=0.2))
        @test_skip run(sess, gamma(shape...;concentration=3.0, rate=2.0))
        @test_skip run(sess, truncatednormal(shape...;loc=0.0, scale=1.0, low=-1., high=1.))
        @test_skip run(sess, chi(shape...;df=5))
        @test_skip run(sess, chi2(shape...;df=5))
        @test_skip run(sess, exponential(shape...;rate=1.))
        @test_skip run(sess, gumbel(shape...;loc=0, scale=1))
        @test_skip run(sess, poisson(shape...;rate=1.0))
        @test_skip run(sess, halfcauchy(shape...;loc=0, scale=1))
        @test_skip run(sess, halfnormal(shape...;scale=1))
        @test_skip run(sess, horseshoe(shape...;scale=1.0))
        @test_skip run(sess, inversegamma(shape...;concentration=1.0, scale=1.0))
        @test_skip run(sess, inversegaussian(shape...;loc=1.0, concentration=1.0))
        @test_skip run(sess, kumaraswamy(shape...;concentration1=1.0, concentration0=1.0))
        @test_skip run(sess, pareto(shape...;scale=1.0, concentration=1.0))
        @test_skip run(sess, sinharcsinh(shape...;skewness=1.0, tailweight=1.0, loc=0.0, scale=1.0))
        @test_skip run(sess, studentt(shape...;loc=0, scale=1, df=10))
        @test_skip run(sess, vonmises(shape...;concentration=1.0, loc=1.0))
    end
end


text:
include("../client/gui/graphs/graphs.jl")

using .Graphs

function test_positive_linear()
    println("Testing positive range")

    data = [3, 10, 15, 20, 22, 26, 30]

    max = findmax(data)
    min = findmin(data)

    l = linear(Float64(min[1]), Float64(max[1]), 3.0)
    @assert l == 0.0 "Expected 0.0"

    l = linear(Float64(min[1]), Float64(max[1]), 30.0)
    @assert l == 1.0 "Expected 1.0"

    println("Passed")
end

function test_negative_linear()
    println("Testing negative range")

    data = [-10, -3, 0, 5, 20, 22, 26, 30]

    max = findmax(data)
    min = findmin(data)
    # println(min, ", ", max)

    l = linear(Float64(min[1]), Float64(max[1]), -10.0)
    @assert l == 0.0 "Expected 0.0"

    l = linear(Float64(min[1]), Float64(max[1]), 30.0)
    @assert l == 1.0 "Expected 1.0"

    println("Passed")
end

function test_lerp()
    println("Testing lerp")

    l = lerp(0.0, 200.0, 0.5)
    @assert l == 100.0 "Expected 100.0, Got: $l"

    l = lerp(50.0, 200.0, 0.5)
    @assert l == 125.0 "Expected 125.0, Got: $l"

    println("Passed")
end

test_positive_linear()
test_negative_linear()
test_lerp()


text:
export getNumPublishers, getTopic, shutdown

"""
    getNumPublishers(subscriber)

Returns the number of publishers this subscriber is connected to.
"""
getNumPublishers(subscriber::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N}) = @cxx subscriber->getNumPublishers()

"""
    getTopic(subscriber)

Returns the topic that this Subscriber is subscribed on.
"""
getTopic(subscriber::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N}) = unsafe_string(@cxx subscriber->getTopic())

"""
    shutdown(subscriber)

Unsubscribe the callback associated with this Subscriber. More...
"""
shutdown(subscriber::Union{Cxx.CxxCore.CppPtr{Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N,(false, false, false)}, Cxx.CxxCore.CppValue{Cxx.CxxCore.CxxQualType{Cxx.CxxCore.CppBaseType{Symbol("ros::Subscriber")},(false, false, false)},N} where N}) = @cxx subscriber->shutdown()


text:
using MatchCore

comm_monoid = @commutative_monoid (*) 1

fold_mul = @theory begin
	a::$Number * b::$Number |> a * b
end

t = comm_monoid ∪ fold_mul

@testset "Extraction 1 - Commutative Monoid" begin
    G = EGraph(:(3 * 4))
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, t)
    @test (12 == extract!(G, extran))

    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, t)
    extr = extract!(G, extran)
	# println(extr)

    @test extr == :(b * (a * 12))
end

fold_add = @theory begin
	a::$Number + b::$Number |> a + b
end

@testset "Extraction 2" begin
	comm_group = @abelian_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = cleanast(:((x*(a+b)) + (y*(a+b))))
    G = EGraph(ex)
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, t)
    # end

    extract!(G, extran) == :((x+y) * (b+a))
end

@testset "Lazy Extraction 2" begin
	comm_group = @abelian_group (+) 0 inv

    t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ fold_mul ∪ fold_add

    # for i ∈ 1:20
    # sleep(0.3)
    ex = cleanast(:((x*(a+b)) + (y*(a+b))))
    G = EGraph(ex)
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, t)
    # end

    extract!(G, extran) == :((x+y) * (b+a))
end

@testset "Extraction - Adding analysis after saturation" begin
    G = EGraph(:(3 * 4))
    addexpr!(G, 12)
    saturate!(G, t)
    addexpr!(G, :(a * 2))
    saturate!(G, t)

    extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, t)

    @test (12 == extract!(G, extran))

    # for i ∈ 1:100
    ex = :(a * 3 * b * 4)
    G = EGraph(cleanast(ex))
    addanalysis!(G, NumberFold)
    extran = addanalysis!(G, ExtractionAnalysis, astsize)
    saturate!(G, comm_monoid)

    extr = extract!(G, extran)
    # end

    @test extr == :((12 * a) * b)
end

comm_monoid = @commutative_monoid (*) 1

comm_group = @abelian_group (+) 0 inv

powers = @theory begin
	a * a => a^2
	a => a^1
	a^n * a^m => a^(n+m)
end
logids = @theory begin
	log(a^n) => n * log(a)
	log(x * y) => log(x) * log(y)
	log(1) => 0
	log(:e) => 1
	:e^(log(x)) => x
end

t = comm_monoid ∪ comm_group ∪ distrib(:(*), :(+)) ∪ powers ∪ logids  ∪ fold_mul ∪ fold_add

@testset "Complex Extraction" begin
	G = EGraph(:(log(e) * log(e)))
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
	saturate!(G, t; timeout=7)
	@test extract!(G, extran) == 1

	G = EGraph(:(log(e) * (log(e) * e^(log(3)))  ))
	extran = addanalysis!(G, ExtractionAnalysis, astsize)
	saturate!(G, t; timeout=7)
	@test extract!(G, extran) == 3

	@time begin
		G = EGraph(:(a^3 * a^2))
		extran = addanalysis!(G, ExtractionAnalysis, astsize)
		saturate!(G, t; timeout=7)
		ex = extract!(G, extran)
	end
	@test ex == :(a^5)

	@time begin
		G = EGraph(:(a^3 * a^2))
		extran = addanalysis!(G, ExtractionAnalysis, astsize)
		saturate!(G, t; timeout=7)
		ex = extract!(G, extran)
	end
	@test ex == :(a^5)

	@time begin
		G = EGraph(:((log(e) * log(e)) * (log(a^3 * a^2))))
		extran = addanalysis!(G, ExtractionAnalysis, astsize)
		saturate!(G, t; timeout=7)
		ex = extract!(G, extran)
	end
	@test ex == :(5*log(a))
end

# EXTRACTION BUG!

costfun(n) = 1
costfun(n::Expr) = n.args[2] == :a ? 1 : 100

moveright = @theory begin
    (:b * (:a * c)) => (:a * (:b * c))
end

expr = :(a * (a * (b * (a * b))))
res = rewrite( expr , moveright)
println(res)

g = EGraph(expr)
saturate!(g, moveright)
extractor = addanalysis!(g, ExtractionAnalysis, costfun)
resg = extract!(g, extractor)
println(resg)

@testset "Symbols in Right hand" begin
    @test resg == res == :(a * (a * (a * (b * b))))
end


text:
@testset "PSSE Parsing " begin
    files = readdir(PSSE_RAW_DIR)
    if length(files) == 0
        error("No test files in the folder")
    end

    for f in files[1:1]
        @info "Parsing $f ..."
        pm_dict = PowerSystems.parse_file(joinpath(PSSE_RAW_DIR, f))
        @info "Successfully parsed $f to PowerModels dict"
        sys = PowerSystems.pm2ps_dict(pm_dict)
        @info "Successfully parsed $f to System struct"
    end

    # Test bad input
    pm_dict = PowerSystems.parse_file(joinpath(PSSE_RAW_DIR, files[1]))
    pm_dict["bus"] = Dict{String, Any}()
    @test_throws PowerSystems.DataFormatError PowerSystems.pm2ps_dict(pm_dict)
end

@testset "PSSE PowerModel Parsing " begin
    PM_PSSE_PATH = joinpath(DATA_DIR, "pm_data","pti")
    files = readdir(PM_PSSE_PATH)
    if length(files) == 0
        error("No test files in the folder")
    end

    for f in files[1:1]
        @info "Parsing $f ..."
        pm_dict = PowerSystems.parse_file(joinpath(PM_PSSE_PATH, f))
        @info "Successfully parsed $f to PowerModels dict"

        sys = PowerSystems.pm2ps_dict(pm_dict)
        @info "Successfully parsed $f to System struct"
    end
end


text:
using JuAFEM
using Tensors
using KrylovMethods
using TimerOutputs
import ProgressMeter

const ∇ = Tensors.gradient;

structNeoHook{T}
    μ::T
    λ::T
end

ProjDir = dirname(@__FILE__)
cd(ProjDir) do

  function compute_2nd_PK(mp, E)
      I = one(E)
      C = 2E + one(E)
      invC = inv(C)
      J = sqrt(det(C))
      return mp.μ *(I - invC) + mp.λ * log(J) * invC
  end

  function constitutive_driver(mp, E)
      ∂S∂E, SPK = ∇(E -> compute_2nd_PK(mp, E), E, :all)
      return SPK, ∂S∂E
  end;

  # Loop over all cells
  function assemble{dim}(grid::Grid{dim}, dh::DofHandler, K, f, cv, fv, mp, u)
      n = ndofs_per_cell(dh)
      Ke = zeros(n, n)
      fe = zeros(n)

      assembler = start_assemble(K, f)

      # loop over all cells in the grid
      @timeit "assemble" for cell in CellIterator(dh)
          # reset
          fill!(Ke, 0)
          fill!(fe, 0)

          global_dofs = celldofs(cell)
          ue = u[global_dofs] # element dofs
          @timeit "inner assemble" assemble_element!(Ke, fe, cell, cv, fv, mp, ue)

          assemble!(assembler, global_dofs, fe, Ke)
      end

      return f, K
  end;

  # Assembles the contribution from the cell to ke and fe
  function assemble_element!(ke, fe, cell, cv, fv, mp, ue)
      b = Vec{3}((0.0, -0.5, 0.0))
      t = Vec{3}((0.1, 0.0, 0.0))
      ndofs = getnbasefunctions(cv)
      reinit!(cv, cell)
      fill!(ke, 0.0)
      fill!(fe, 0.0)
      δE = Vector{SymmetricTensor{2, 3, eltype(ue), 6}}(ndofs)

      for qp in 1:getnquadpoints(cv)
          ∇u = function_gradient(cv, qp, ue)
          dΩ = getdetJdV(cv, qp)

          # strain and stress + tangent
          F = one(∇u) + ∇u
          E = symmetric(1/2 * (F' ⋅ F - one(F)))

          S, ∂S∂E = constitutive_driver(mp, E)

          # Hoist computations of δE
          for i in 1:ndofs
              δFi = shape_gradient(cv, qp, i)
              δE[i] = symmetric(1/2*(δFi'⋅F + F'⋅δFi))
          end

          for i in 1:ndofs
              δFi = shape_gradient(cv, qp, i)
              δu = shape_value(cv, qp, i)
              fe[i] += (δE[i] ⊡ S) * dΩ
              fe[i] -= (δu ⋅ b) * dΩ
              δE∂S∂E = δE[i] ⊡ ∂S∂E
              S∇δu = S ⋅ δFi'
              for j in 1:ndofs
                  δ∇uj = shape_gradient(cv, qp, j)
                  ke[i, j] += (δE∂S∂E ⊡ δE[j] + S∇δu ⊡ δ∇uj' ) * dΩ
              end
          end
      end

      for face in 1:nfaces(cell)
          if onboundary(cell, face)
              reinit!(fv, cell, face)
              for q_point in 1:getnquadpoints(fv)
                  dΓ = getdetJdV(fv, q_point)
                  for i in 1:ndofs
                      δu = shape_value(fv, q_point, i)
                      fe[i] -= (δu ⋅ t) * dΓ
                  end
              end
          end
      end
  end;

  function solve()
      reset_timer!()

      const dim = 3

      # Generate a grid
      N = 10
      L = 1.0
      left = zero(Vec{dim})
      right = L * ones(Vec{dim})
      grid = generate_grid(Tetrahedron, ntuple(x->N, dim), left, right)

      # Node sets
      addnodeset!(grid, "clamped", x -> norm(x[1]) ≈ 1)
      addnodeset!(grid, "rotation", x -> norm(x[1]) ≈ 0)

      # Material parameters
      E = 10.0
      ν = 0.3
      μ = E / (2(1 + ν))
      λ = (E * ν) / ((1 + ν) * (1 - 2ν))
      mp = NeoHook(μ, λ)

      # finite element base
      ip = Lagrange{dim, RefTetrahedron, 1}()
      qr = QuadratureRule{dim, RefTetrahedron}(1)
      qr_face = QuadratureRule{dim-1, RefTetrahedron}(1)
      cv = CellVectorValues(qr, ip)
      fv = FaceVectorValues(qr_face, ip)

      # DofHandler
      dh = DofHandler(grid)
      push!(dh, :u, dim) # Add a displacement field
      close!(dh)

      function rotation(X, t, θ = deg2rad(60.0))
          x, y, z = X
          return t * Vec{dim}(
              (0.0,
              L/2 - y + (y-L/2)*cos(θ) - (z-L/2)*sin(θ),
              L/2 - z + (y-L/2)*sin(θ) + (z-L/2)*cos(θ)
              ))
      end

      dbc = DirichletBoundaryConditions(dh)
      # Add a homogenoush boundary condition on the "clamped" edge
      add!(dbc, :u, getnodeset(grid, "clamped"), (x,t) -> [0.0, 0.0, 0.0], collect(1:dim))
      add!(dbc, :u, getnodeset(grid, "rotation"), (x,t) -> rotation(x, t), collect(1:dim))
      close!(dbc)
      t = 0.5
      update!(dbc, t)

      println("Analysis with ", length(grid.cells), " elements")

      # pre-allocate
      _ndofs = ndofs(dh)
      un = zeros(_ndofs) # previous solution vector
      u  = zeros(_ndofs)
      Δu = zeros(_ndofs)

      apply!(un, dbc)

      K = create_sparsity_pattern(dh)
      f = zeros(_ndofs)

      newton_itr = -1
      NEWTON_TOL = 1e-8
      prog = ProgressMeter.ProgressThresh(NEWTON_TOL, "Solving:")

      while true; newton_itr += 1
          u .= un .+ Δu
          f, K = assemble(grid, dh, K, f, cv, fv, mp, u)
          normg = norm(f[JuAFEM.free_dofs(dbc)])
          apply_zero!(K, f, dbc)
          ProgressMeter.update!(prog, normg; showvalues = [(:iter, newton_itr)])

          if normg < NEWTON_TOL
              break
          end

          if newton_itr > 30
              error("Reached maximum Newton iterations, aborting")
              break
          end

          @timeit "linear solve" ΔΔu, flag, relres, iter, resvec = cg(K, f; maxIter = 1000, tol = min(1e-3, normg))
          @assert flag == 0

          apply_zero!(ΔΔu, dbc)
          Δu .-= ΔΔu
      end

      # save the solution
      @timeit "export" begin
          vtk = vtk_grid("hyperelasticity", dh, u)
          vtk_save(vtk)
      end

      print_timer(linechars = :ascii)
      return u
  end

  u = solve();

  Base.Test.@test norm(u) ≈ 4.870833706518008
  println("Hyperelasticity passed!")

end


text:
using RandomMatrices
using Random: seed!
using Test

@testset "RandomMatrices" begin
seed!(1)
include("GaussianEnsembles.jl")
include("FormalPowerSeries.jl")
include("Haar.jl")
include("StochasticProcess.jl")

@testset "densities" begin
    include("densities/Semicircle.jl")
    include("densities/TracyWidom.jl")
end
end


text:
export SimpleBoreCollisionModule, SimpleBoreCollisionModuleParams

Base.@kwdef struct SimpleBoreCollisionModuleParams <: DeviceParams
  "Diameter of scanner in the y-z plane"
  scannerDiameter::typeof(1.0u"mm")
  "Geometry of the probe, centered at (0,0,0) in scanner coordinates"
  objGeometry::Geometry

  clearance::Clearance = Clearance(1.0u"mm")

  "Define the minimum and maximum points of movement into and out of the bore axis (x)"
  minMaxBoreAxis::Vector{typeof(1.0u"mm")} = [-Inf, Inf]u"mm"
end

SimpleBoreCollisionModuleParams(dict::Dict) = params_from_dict(SimpleBoreCollisionModuleParams, dict)

Base.@kwdef mutable struct SimpleBoreCollisionModule <: AbstractCollisionModule3D
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimpleBoreCollisionModuleParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String,Union{Device,Missing}}
end

collisionModuleType(cm::SimpleBoreCollisionModule) = PositionCollisionType()

function _checkCoordinate(cm::SimpleBoreCollisionModule, pos::AbstractVector{<:Unitful.Length}; returnVerbose=false)
  validStatus = Array{Bool}(undef, 3)
  errorDiff = Array{typeof(1.0u"mm")}(undef, 3)

  validStatus[1], errorDiff[1] = _checkCoordsBoreAxis(cm, pos)
  tmp = _checkCoordsCrossSection(cm, pos)

  errorDiff[2:3] .= tmp[2:3]
  validStatus[2:3] .= tmp[1]

  return validStatus, errorDiff
end

function _checkCoordsBoreAxis(cm::SimpleBoreCollisionModule, pos::AbstractVector{<:Unitful.Length})
  posBore = pos[1]
  if posBore < minimum(cm.params.minMaxBoreAxis)
    return false, posBore - minimum(cm.params.minMaxBoreAxis)
  elseif posBore > maximum(cm.params.minMaxBoreAxis)
    return false, posBore - maximum(cm.params.minMaxBoreAxis)
  else
    return true, zero(0.0u"mm")
  end
end

_checkCoordsCrossSection(cm::SimpleBoreCollisionModule, pos::AbstractVector{<:Unitful.Length}) = checkCollisionYZCircle(cm.params.objGeometry, cm.params.scannerDiameter / 2, pos[2], pos[3], cm.params.clearance.distance)

function plotSafetyErrors(cm::SimpleBoreCollisionModule, errorIndices::Vector{Int}, coords::AbstractMatrix{<:Unitful.Length})
  geo = cm.params.objGeometry;
  scannerRad = cm.params.scannerDiameter / 2;

  t = range(0, stop=2, length=200);

  x_scanner = scannerRad * cos.(t * pi);
  y_scanner = scannerRad * sin.(t * pi);
  x_scanner2 = (scannerRad - cm.params.clearance.distance) * cos.(t * pi);
  y_scanner2 = (scannerRad - cm.params.clearance.distance) * sin.(t * pi);

  fig = plot(title="Plot results - $(geo.name) positions", xlabel="y [mm]", ylabel="z [mm]", aspect_ratio=:equal)
  plot!(ustrip.(u"mm", x_scanner), ustrip.(u"mm", y_scanner), color=:blue, label="scanner")
  plot!(ustrip.(u"mm", x_scanner2), ustrip.(u"mm", y_scanner2), color=:yellow, label="scanner, with clearance")
  for i = errorIndices
    y_i = coords[i, 2];
    z_i = coords[i, 3];

    if typeof(geo) == Circle
      x_geometry = geo.diameter / 2 * cos.(t * pi) .+ y_i;
      y_geometry = geo.diameter / 2 * sin.(t * pi) .+ z_i;

      plot!(ustrip.(u"mm", x_geometry), ustrip.(u"mm", y_geometry), color=:red)

    elseif typeof(geo) == Rectangle
      # Create rectangle corner points
      # point bottom left
      p_bl = ustrip.(u"mm", [y_i - geo.width / 2, z_i - geo.height / 2]);
      # point upper left
      p_ul = ustrip.(u"mm", [y_i - geo.width / 2, z_i + geo.height / 2]);
      # point upper right
      p_ur = ustrip.(u"mm", [y_i + geo.width / 2, z_i + geo.height / 2]);
      # point bottom right
      p_br = ustrip.(u"mm", [y_i + geo.width / 2, z_i - geo.height / 2]);

      rect = transpose([p_bl p_ul p_ur p_br p_bl])
      plot!(rect[:,1], rect[:,2], color=:red);

    elseif typeof(geo) == Triangle
      # Create triangle corner points
      # point bottom left
      p_bl = ustrip.(u"mm", [y_i - geo.width / 2, z_i - geo.height / 3]);
      # upper point
      p_u = ustrip.(u"mm", [y_i, z_i + 2 / 3 * geo.height]);
      # point bottom right
      p_br = ustrip.(u"mm", [y_i + geo.width / 2, z_i - geo.height / 3]);

      tri = [p_bl p_u p_br p_bl]
      plot!(tri[1,:], tri[2,:], color=:red, label=nothing);
    end

  end
  display(fig)
end