text:
const MYID = myid()
const OTHERIDS = filter(id-> id != MYID, procs())[rand(1:(nprocs()-1))]

function check_leaks()
    if length(DistributedArrays.refs) > 0
        sleep(0.1)  # allow time for any cleanup to complete and test again
        length(DistributedArrays.refs) > 0 && warn("Probable leak of ", length(DistributedArrays.refs), " darrays")
    end
end

@testset "test distribute" begin
    A = rand(1:100, (100,100))

    @testset "test default distribute" begin
        DA = distribute(A)
        @test length(procs(DA)) == nworkers()
        @test sum(DA) == sum(A)
        close(DA)
    end

    @testset "test distribute with procs arguments" begin
        DA = distribute(A, procs = procs())
        @test length(procs(DA)) == nprocs()
        @test sum(DA) == sum(A)
        close(DA)
    end

    @testset "test distribute with procs and dist arguments" begin
        DA = distribute(A, procs = [1, 2], dist = [1,2])
        @test size(procs(DA)) == (1,2)
        @test sum(DA) == sum(A)
        close(DA)
    end

    @testset "Create darray with unconventional distribution and distibute like it" begin
        block = 10
        Y = nworkers() * block
        X = nworkers() * block
        remote_parts = map(workers()) do wid
            remotecall(rand, wid, block, Y)
        end
        DA1 = DArray(reshape(remote_parts, (length(remote_parts), 1)))
        A = rand(X, Y)
        DA2 = distribute(A, DA1)

        @test size(DA1) == size(DA2)

        close(DA1)
        close(DA2)
    end
end

check_leaks()

@testset "test DArray equality" begin
    D = drand((200,200), [MYID, OTHERIDS])
    DC = copy(D)

    @testset "test isequal(::DArray, ::DArray)" begin
        @test D == DC
    end

    @testset "test copy(::DArray) does a copy of each localpart" begin
        @spawnat OTHERIDS localpart(DC)[1] = 0
        @test fetch(@spawnat OTHERIDS localpart(D)[1] != 0)
    end

    close(D)
    close(DC)
end

check_leaks()

@testset "test DArray similar" begin
    D = drand((200,200), [MYID, OTHERIDS])
    DS = similar(D,Float16)

    @testset "test eltype of a similar" begin
        @test eltype(DS) == Float16
    end

    @testset "test dims of a similar" begin
        @test size(D) == size(DS)
    end
    close(D)
    close(DS)
end

check_leaks()

@testset "test DArray reshape" begin
    D = drand((200,200), [MYID, OTHERIDS])

    @testset "Test error-throwing in reshape" begin
        @test_throws DimensionMismatch reshape(D,(100,100))
    end

    DR = reshape(D,(100,400))
    @testset "Test reshape" begin
        @test size(DR) == (100,400)
    end
    close(D)
end

check_leaks()

@testset "test @DArray comprehension constructor" begin

    @testset "test valid use of @DArray" begin
        D = @DArray [i+j for i=1:10, j=1:10]
        @test D == [i+j for i=1:10, j=1:10]
        close(D)
    end

    @testset "test invalid use of @DArray" begin
        @test_throws ArgumentError eval(:((@DArray [1,2,3,4])))
    end
end

check_leaks()

@testset "test DArray / Array conversion" begin
    D = drand((200,200), [MYID, OTHERIDS])

    @testset "test convert(::Array, ::(Sub)DArray)" begin
        S = convert(Matrix{Float64}, D[1:150, 1:150])
        A = convert(Matrix{Float64}, D)

        @test A[1:150,1:150] == S
        D2 = convert(DArray{Float64,2,Matrix{Float64}}, A)
        @test D2 == D
        @test fetch(@spawnat MYID localpart(D)[1,1]) == D[1,1]
        @test fetch(@spawnat OTHERIDS localpart(D)[1,1]) == D[1,101]
        close(D2)

        S2 = convert(Vector{Float64}, D[4, 23:176])
        @test A[4, 23:176] == S2

        S3 = convert(Vector{Float64}, D[23:176, 197])
        @test A[23:176, 197] == S3

        S4 = zeros(4)
        setindex!(S4, D[3:4, 99:100], :)
        @test S4 == vec(D[3:4, 99:100])
        @test S4 == vec(A[3:4, 99:100])

        S5 = zeros(2,2)
        setindex!(S5, D[1,1:4], :, 1:2)
        @test vec(S5) == D[1, 1:4]
        @test vec(S5) == A[1, 1:4]
    end
    close(D)
end

check_leaks()

@testset "copy!" begin
    D1 = dzeros((10,10))
    r1 = remotecall_wait(() -> randn(3,10), workers()[1])
    r2 = remotecall_wait(() -> randn(7,10), workers()[2])
    D2 = DArray(reshape([r1; r2], 2, 1))
    copy!(D2, D1)
    @test D1 == D2
    close(D1)
    close(D2)
end

check_leaks()

@testset "test DArray reduce" begin
    D = DArray(id->fill(myid(), map(length,id)), (10,10), [MYID, OTHERIDS])

    @testset "test reduce" begin
        @test reduce(+, D) == ((50*MYID) + (50*OTHERIDS))
    end

    @testset "test map / reduce" begin
        D2 = map(x->1, D)
        @test reduce(+, D2) == 100
        close(D2)
    end

    @testset "test map! / reduce" begin
        map!(x->1, D)
        @test reduce(+, D) == 100
    end
    close(D)
end

check_leaks()

@testset "test scale" begin
    A = randn(100,100)
    DA = distribute(A)
    @test scale!(DA, 2) == scale!(A, 2)
    close(DA)
end

check_leaks()

@testset "test scale!(b, A)" begin
    A = randn(100, 100)
    b = randn(100)
    DA = distribute(A)
    @test scale!(b, A) == scale!(b, DA)
    close(DA)
    A = randn(100, 100)
    b = randn(100)
    DA = distribute(A)
    @test scale!(A, b) == scale!(DA, b)
    close(DA)
end

check_leaks()

@testset "test mapreduce on DArrays" begin
    for _ = 1:25, f = [x -> Int128(2x), x -> Int128(x^2), x -> Int128(x^2 + 2x - 1)], opt = [+, *]
        A = rand(1:5, rand(2:30))
        DA = distribute(A)
        @test mapreduce(f, opt, DA) - mapreduce(f, opt, A) == 0
        close(DA)
    end
end

check_leaks()

@testset "test mapreducedim on DArrays" begin
    D = DArray(I->fill(myid(), map(length,I)), (73,73), [MYID, OTHERIDS])
    D2 = map(x->1, D)
    @test mapreducedim(t -> t*t, +, D2, 1) == mapreducedim(t -> t*t, +, convert(Array, D2), 1)
    @test mapreducedim(t -> t*t, +, D2, 2) == mapreducedim(t -> t*t, +, convert(Array, D2), 2)
    @test mapreducedim(t -> t*t, +, D2, (1,2)) == mapreducedim(t -> t*t, +, convert(Array, D2), (1,2))

    # Test non-regularly chunked DArrays
    r1 = DistributedArrays.remotecall(() -> sprandn(3, 10, 0.1), workers()[1])
    r2 = DistributedArrays.remotecall(() -> sprandn(7, 10, 0.1), workers()[2])
    D = DArray(reshape([r1; r2], (2,1)))
    @test Array(sum(D, 2)) == sum(Array(D), 2)

    # close(D)
    # close(D2)
    darray_closeall()   # temp created by the mapreduce above
end

check_leaks()

@testset "test mapreducdim, reducedim on DArrays" begin
    dims = (20,20,20)
    DA = drandn(dims)
    A = convert(Array, DA)

    @testset "dimension $dms" for dms in (1, 2, 3, (1,2), (1,3), (2,3), (1,2,3))
        @test mapreducedim(t -> t*t, +, A, dms) ≈ mapreducedim(t -> t*t, +, DA, dms)
        @test mapreducedim(t -> t*t, +, A, dms, 1.0) ≈ mapreducedim(t -> t*t, +, DA, dms, 1.0)
        @test reducedim(*, A, dms) ≈ reducedim(*, DA, dms)
        @test reducedim(*, A, dms, 2.0) ≈ reducedim(*, DA, dms, 2.0)
    end
    close(DA)
    darray_closeall()   # temp created by the mapreduce above
end

check_leaks()

@testset "test statistical functions on DArrays" begin
    dims = (20,20,20)
    DA = drandn(dims)
    A = convert(Array, DA)

    @testset "test $f for dimension $dms" for f in (mean, ), dms in (1, 2, 3, (1,2), (1,3), (2,3), (1,2,3))
        # std is pending implementation
        @test f(DA,dms) ≈ f(A,dms)
    end

    close(DA)
    darray_closeall()   # temporaries created above
end

check_leaks()

@testset "test sum on DArrays" begin
    A = randn(100,100)
    DA = distribute(A)

    # sum either throws an ArgumentError or a CompositeException of ArgumentErrors
    try
        sum(DA, -1)
    catch err
        if isa(err, CompositeException)
            @test !isempty(err.exceptions)
            for excep in err.exceptions
                # Unpack the remote exception
                orig_err = excep.ex.captured.ex
                @test isa(orig_err, ArgumentError)
            end
        else
            @test isa(err, ArgumentError)
        end
    end
    try
        sum(DA, 0)
    catch err
        if isa(err, CompositeException)
            @test !isempty(err.exceptions)
            for excep in err.exceptions
                # Unpack the remote exception
                orig_err = excep.ex.captured.ex
                @test isa(orig_err, ArgumentError)
            end
        else
            @test isa(err, ArgumentError)
        end
    end

    @test sum(DA) ≈ sum(A)
    @test sum(DA,1) ≈ sum(A,1)
    @test sum(DA,2) ≈ sum(A,2)
    @test sum(DA,3) ≈ sum(A,3)
    close(DA)
    darray_closeall()   # temporaries created above
end

check_leaks()

@testset "test size on DArrays" begin

    A = randn(100,100)
    DA = distribute(A)

    @test_throws BoundsError size(DA, 0)
    @test size(DA,1) == size(A,1)
    @test size(DA,2) == size(A,2)
    @test size(DA,3) == size(A,3)
    close(DA)
end

check_leaks()

# test length / endof
@testset "test collections API" begin
    A = randn(23,23)
    DA = distribute(A)

    @testset "test length" begin
        @test length(DA) == length(A)
    end

    @testset "test endof" begin
        @test endof(DA) == endof(A)
    end
    close(DA)
end

check_leaks()

@testset "test max / min / sum" begin
    a = map(x->Int(round(rand() * 100)) - 50, Array(Int, 100,1000))
    d = distribute(a)

    @test sum(d) == sum(a)
    @test maximum(d) == maximum(a)
    @test minimum(d) == minimum(a)
    @test maxabs(d) == maxabs(a)
    @test minabs(d) == minabs(a)
    @test sumabs(d) == sumabs(a)
    @test sumabs2(d) == sumabs2(a)
    close(d)
end

check_leaks()

@testset "test all / any" begin
    a = map(x->Int(round(rand() * 100)) - 50, Array(Int, 100,1000))
    a = [true for i in 1:100]
    d = distribute(a)

    @test all(d)
    @test any(d)

    close(d)

    a[50] = false
    d = distribute(a)
    @test !all(d)
    @test any(d)

    close(d)

    a = [false for i in 1:100]
    d = distribute(a)
    @test !all(d)
    @test !any(d)

    close(d)

    d = dones(10,10)
    @test !all(x-> x>1.0, d)
    @test all(x-> x>0.0, d)

    close(d)

    a = ones(10,10)
    a[10] = 2.0
    d = distribute(a)
    @test any(x-> x == 1.0, d)
    @test any(x-> x == 2.0, d)
    @test !any(x-> x == 3.0, d)

    close(d)
end

check_leaks()

@testset "test count"  begin
    a = ones(10,10)
    a[10] = 2.0
    d = distribute(a)

    @test count(x-> x == 2.0, d) == 1
    @test count(x-> x == 1.0, d) == 99
    @test count(x-> x == 0.0, d) == 0

    close(d)
end

check_leaks()

@testset "test prod" begin
    a = fill(2, 10);
    d = distribute(a);
    @test prod(d) == 2^10

    close(d)
end

check_leaks()

@testset "test zeros" begin
    @testset "1D dzeros default element type" begin
        A = dzeros(10)
        @test A == zeros(10)
        @test eltype(A) == Float64
        @test size(A) == (10,)
        close(A)
    end

    @testset "1D dzeros with specified element type" begin
        A = dzeros(Int, 10)
        @test A == zeros(10)
        @test eltype(A) == Int
        @test size(A) == (10,)
        close(A)
    end

    @testset "2D dzeros default element type, Dims constuctor" begin
        A = dzeros((10,10))
        @test A == zeros((10,10))
        @test eltype(A) == Float64
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dzeros specified element type, Dims constructor" begin
        A = dzeros(Int, (10,10))
        @test A == zeros(Int, (10,10))
        @test eltype(A) == Int
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dzeros, default element type" begin
        A = dzeros(10,10)
        @test A == zeros(10,10)
        @test eltype(A) == Float64
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dzeros, specified element type" begin
        A = dzeros(Int, 10, 10)
        @test A == zeros(Int, 10, 10)
        @test eltype(A) == Int
        @test size(A) == (10,10)
        close(A)
    end
end

check_leaks()

@testset "test dones" begin
    @testset "1D dones default element type" begin
        A = dones(10)
        @test A == ones(10)
        @test eltype(A) == Float64
        @test size(A) == (10,)
        close(A)
    end

    @testset "1D dones with specified element type" begin
        A = dones(Int, 10)
        @test eltype(A) == Int
        @test size(A) == (10,)
        close(A)
    end

    @testset "2D dones default element type, Dims constuctor" begin
        A = dones((10,10))
        @test A == ones((10,10))
        @test eltype(A) == Float64
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dones specified element type, Dims constructor" begin
        A = dones(Int, (10,10))
        @test A == ones(Int, (10,10))
        @test eltype(A) == Int
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dones, default element type" begin
        A = dones(10,10)
        @test A == ones(10,10)
        @test eltype(A) == Float64
        @test size(A) == (10,10)
        close(A)
    end

    @testset "2D dones, specified element type" begin
        A = dones(Int, 10, 10)
        @test A == ones(Int, 10, 10)
        @test eltype(A) == Int
        @test size(A) == (10,10)
        close(A)
    end
end

check_leaks()

@testset "test drand" begin
    @testset "1D drand" begin
        A = drand(100)
        @test eltype(A) == Float64
        @test size(A) == (100,)
        @test all(x-> x >= 0.0 && x <= 1.0, A)
        close(A)
    end

    @testset "1D drand, specified element type" begin
        A = drand(Int, 100)
        @test eltype(A) == Int
        @test size(A) == (100,)
        close(A)
    end

    @testset "2D drand, Dims constructor" begin
        A = drand((50,50))
        @test eltype(A) == Float64
        @test size(A) == (50,50)
        @test all(x-> x >= 0.0 && x <= 1.0, A)
        close(A)
    end

    @testset "2D drand" begin
        A = drand(100,100)
        @test eltype(A) == Float64
        @test size(A) == (100,100)
        @test all(x-> x >= 0.0 && x <= 1.0, A)
        close(A)
    end

    @testset "2D drand, Dims constructor, specified element type" begin
        A = drand(Int, (100,100))
        @test eltype(A) == Int
        @test size(A) == (100,100)
        close(A)
    end

    @testset "2D drand, specified element type" begin
        A = drand(Int, 100, 100)
        @test eltype(A) == Int
        @test size(A) == (100,100)
        close(A)
    end
end

check_leaks()

@testset "test randn" begin
    @testset "1D drandn" begin
        A = drandn(100)
        @test eltype(A) == Float64
        @test size(A) == (100,)
        close(A)
    end

    @testset "2D drandn, Dims constructor" begin
        A = drandn((50,50))
        @test eltype(A) == Float64
        @test size(A) == (50,50)
        close(A)
    end

    @testset "2D drandn" begin
        A = drandn(100,100)
        @test eltype(A) == Float64
        @test size(A) == (100,100)
        close(A)
    end
end

check_leaks()

@testset "test c/transpose" begin
    @testset "test ctranspose real" begin
        A = drand(Float64, 100, 200)
        @test A' == Array(A)'
        close(A)
    end
    @testset "test ctranspose complex" begin
        A = drand(Complex128, 200, 100)
        @test A' == Array(A)'
        close(A)
    end
    @testset "test transpose real" begin
        A = drand(Float64, 200, 100)
        @test A.' == Array(A).'
        close(A)
    end
    @testset "test ctranspose complex" begin
        A = drand(Complex128, 100, 200)
        @test A.' == Array(A).'
        close(A)
    end

    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test convert from subdarray" begin
    a = drand(20, 20);

    s = view(a, 1:5, 5:8)
    @test isa(s, SubDArray)
    @test s == convert(DArray, s)

    s = view(a, 6:5, 5:8)
    @test isa(s, SubDArray)
    @test s == convert(DArray, s)
    close(a)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test scalar math" begin
    a = drand(20, 20);
    b = convert(Array, a)
    @testset "$f" for f in (-, abs, abs2, acos, acosd, acot,
              acotd, acsch, angle, asech, asin,
              asind, asinh, atan, atand, atanh,
              big, cbrt, ceil, cis, complex, conj,
              cos, cosc, cosd, cosh, cospi, cot,
              cotd, coth, csc, cscd, csch, dawson,
              deg2rad, digamma, erf, erfc, erfcinv,
              erfcx, erfi, erfinv, exp, exp10, exp2,
              expm1, exponent, float, floor, gamma, imag,
              invdigamma, isfinite, isinf, isnan, lfact,
              lgamma, log, log10, log1p, log2, rad2deg, real,
              sec, secd, sech, sign, sin, sinc, sind,
              sinh, sinpi, sqrt, tan, tand, tanh, trigamma)
        @test f.(a) == f.(b)
    end
    a = a + 1
    b = b + 1
    @testset "$f" for f in (asec, asecd, acosh, acsc, acscd, acoth)
        @test f.(a) == f.(b)
    end
    close(a)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test mapslices" begin
    A = randn(5,5,5)
    D = distribute(A, procs = workers(), dist = [1, 1, min(nworkers(), 5)])
    @test mapslices(svdvals, D, (1,2)) ≈ mapslices(svdvals, A, (1,2))
    @test mapslices(svdvals, D, (1,3)) ≈ mapslices(svdvals, A, (1,3))
    @test mapslices(svdvals, D, (2,3)) ≈ mapslices(svdvals, A, (2,3))
    @test mapslices(sort, D, (1,)) ≈ mapslices(sort, A, (1,))
    @test mapslices(sort, D, (2,)) ≈ mapslices(sort, A, (2,))
    @test mapslices(sort, D, (3,)) ≈ mapslices(sort, A, (3,))

    # issue #3613
    B = mapslices(sum, dones(Float64, (2,3,4), workers(), [1,1,min(nworkers(),4)]), [1,2])
    @test size(B) == (1,1,4)
    @test all(B.==6)

    # issue #5141
    C1 = mapslices(x-> maximum(-x), D, [])
    @test C1 == -D

    # issue #5177
    c = dones(Float64, (2,3,4,5), workers(), [1,1,1,min(nworkers(),5)])
    m1 = mapslices(x-> ones(2,3), c, [1,2])
    m2 = mapslices(x-> ones(2,4), c, [1,3])
    m3 = mapslices(x-> ones(3,4), c, [2,3])
    @test size(m1) == size(m2) == size(m3) == size(c)

    n1 = mapslices(x-> ones(6), c, [1,2])
    n2 = mapslices(x-> ones(6), c, [1,3])
    n3 = mapslices(x-> ones(6), c, [2,3])
    n1a = mapslices(x-> ones(1,6), c, [1,2])
    n2a = mapslices(x-> ones(1,6), c, [1,3])
    n3a = mapslices(x-> ones(1,6), c, [2,3])
    @test (size(n1a) == (1,6,4,5) && size(n2a) == (1,3,6,5) && size(n3a) == (2,1,6,5))
    @test (size(n1) == (6,1,4,5) && size(n2) == (6,3,1,5) && size(n3) == (2,6,1,5))
    close(D)
    close(c)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test scalar ops" begin
    a = drand(20,20)
    b = convert(Array, a)
    c = drand(20,20)
    d = convert(Array, c)

    @testset "$f" for f in (+, -, .+, .-, .*, ./, .%, div, mod)
        x = rand()
        @test f(a, x) == f(b, x)
        @test f(x, a) == f(x, b)
        @test f(a, c) == f(b, d)
    end

    close(a)
    close(c)

    a = dones(Int, 20, 20)
    b = convert(Array, a)
    @testset "$f" for f in (.<<, .>>)
        @test f(a, 2) == f(b, 2)
        @test f(2, a) == f(2, b)
        @test f(a, a) == f(b, b)
    end

    @testset "$f" for f in (rem,)
        x = rand()
        @test f(a, x) == f(b, x)
    end
    close(a)
    close(c)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test broadcast ops" begin
    wrkrs = workers()
    nwrkrs = length(wrkrs)
    nrows = 20 * nwrkrs
    ncols = 10 * nwrkrs
    a = drand((nrows,ncols), wrkrs, (1, nwrkrs))
    m = mean(a, 1)
    c = a .- m
    d = convert(Array, a) .- convert(Array, m)
    @test c == d
    darray_closeall()
end

check_leaks()

@testset "test matrix multiplication" begin
    A = drandn(20,20)
    b = drandn(20)
    B = drandn(20,20)

    @test norm(convert(Array, A*b) - convert(Array, A)*convert(Array, b), Inf) < sqrt(eps())
    @test norm(convert(Array, A*B) - convert(Array, A)*convert(Array, B), Inf) < sqrt(eps())
    @test norm(convert(Array, A'*b) - convert(Array, A)'*convert(Array, b), Inf) < sqrt(eps())
    @test norm(convert(Array, A'*B) - convert(Array, A)'*convert(Array, B), Inf) < sqrt(eps())
    close(A)
    close(b)
    close(B)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test norm" begin
    x = drandn(20)

    @test abs(norm(x) - norm(convert(Array, x))) < sqrt(eps())
    @test abs(norm(x, 1) - norm(convert(Array, x), 1)) < sqrt(eps())
    @test abs(norm(x, 2) - norm(convert(Array, x), 2)) < sqrt(eps())
    @test abs(norm(x, Inf) - norm(convert(Array, x), Inf)) < sqrt(eps())
    close(x)
end

check_leaks()

@testset "test axpy!" begin
    x = drandn(20)
    y = drandn(20)

    @test norm(convert(Array, LinAlg.axpy!(2.0, x, copy(y))) - LinAlg.axpy!(2.0, convert(Array, x), convert(Array, y))) < sqrt(eps())
    @test_throws DimensionMismatch LinAlg.axpy!(2.0, x, zeros(length(x) + 1))
    close(x)
    close(y)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test ppeval" begin
    A = drandn((10, 10, nworkers()), workers(), [1, 1, nworkers()])
    B = drandn((10, nworkers()), workers(), [1, nworkers()])

    R = zeros(10, nworkers())
    for i = 1:nworkers()
        R[:, i] = convert(Array, A)[:, :, i]*convert(Array, B)[:, i]
    end
    @test convert(Array, ppeval(*, A, B)) ≈ R
    @test sum(ppeval(eigvals, A)) ≈ sum(ppeval(eigvals, A, eye(10, 10)))
    close(A)
    close(B)
    darray_closeall()  # close the temporaries created above
end

check_leaks()

@testset "test nnz" begin
    A = sprandn(10, 10, 0.5)
    @test nnz(distribute(A)) == nnz(A)
end

@testset "test matmatmul" begin
    A = drandn(30, 30)
    B = drandn(30, 20)
    a = convert(Array, A)
    b = convert(Array, B)

    AB = A * B
    AtB = A.' * B
    AcB = A' * B

    ab = a * b
    atb = a.' * b
    acb = a' * b

    @test AB ≈ ab
    @test AtB ≈ atb
    @test AcB ≈ acb
    darray_closeall()  # close the temporaries created above
end

@testset "sort, T = $T" for i in 0:6, T in [Int, Float64]
    d = DistributedArrays.drand(T, 10^i)
    @testset "sample = $sample" for sample in Any[true, false, (minimum(d),maximum(d)), rand(T, 10^i>512 ? 512 : 10^i)]
        d2 = DistributedArrays.sort(d; sample=sample)

        @test length(d) == length(d2)
        @test sort(convert(Array, d)) == convert(Array, d2)
    end
    darray_closeall()  # close the temporaries created above
end

check_leaks()

darray_closeall()

@testset "test for any leaks" begin
    sleep(1.0)     # allow time for any cleanup to complete
    allrefszero = Bool[remotecall_fetch(()->length(DistributedArrays.refs) == 0, p) for p in procs()]
    @test all(allrefszero)

    allregistrieszero = Bool[remotecall_fetch(()->length(DistributedArrays.registry) == 0, p) for p in procs()]
    @test all(allregistrieszero)
end


text:
#### Land sources
export PhaseChange

function heaviside(x::FT) where {FT}
    if x >= FT(0)
        output = FT(1)
    else
        output = FT(0)
    end
    return output
end

"""
    PhaseChange <: TendencyDef{Source}
The function which computes the freeze/thaw source term for Richard's equation.
"""
Base.@kwdef struct PhaseChange{FT} <: TendencyDef{Source}
    "Typical resolution in the vertical"
    Δz::FT = FT(NaN)
end

prognostic_vars(::PhaseChange) =
    (VolumetricLiquidFraction(), VolumetricIceFraction())

function precompute(land::LandModel, args, tt::Source)
    dtup = DispatchedSet(map(land.source) do s
        (s, precompute(s, land, args, tt))
    end)
    return (; dtup)
end

function precompute(source_type::PhaseChange, land::LandModel, args, tt::Source)
    @unpack state, diffusive, aux, t, direction = args

    FT = eltype(state)

    param_set = parameter_set(land)
    _ρliq = FT(ρ_cloud_liq(param_set))
    _ρice = FT(ρ_cloud_ice(param_set))
    _Tfreeze = FT(T_freeze(param_set))
    _LH_f0 = FT(LH_f0(param_set))
    _g = FT(grav(param_set))

    ϑ_l, θ_i = get_water_content(land.soil.water, aux, state, t)
    ν = land.soil.param_functions.porosity
    eff_porosity = ν - θ_i
    θ_l = volumetric_liquid_fraction(ϑ_l, eff_porosity)
    T = get_temperature(land.soil.heat, aux, t)
    κ_dry = k_dry(land.param_set, land.soil.param_functions)
    S_r = relative_saturation(θ_l, θ_i, ν)
    kersten = kersten_number(θ_i, S_r, land.soil.param_functions)
    κ_sat = saturated_thermal_conductivity(
        θ_l,
        θ_i,
        land.soil.param_functions.κ_sat_unfrozen,
        land.soil.param_functions.κ_sat_frozen,
    )
    κ = thermal_conductivity(κ_dry, kersten, κ_sat)
    ρc_ds = land.soil.param_functions.ρc_ds
    ρc_s = volumetric_heat_capacity(θ_l, θ_i, ρc_ds, land.param_set)
    θ_r = land.soil.param_functions.θ_r

    hydraulics = land.soil.water.hydraulics
    θ_m = min(_ρice * θ_i / _ρliq + θ_l, ν)
    ψ0 = matric_potential(hydraulics, θ_m)
    ψT = _LH_f0 / _g / _Tfreeze * (T - _Tfreeze)
    if T < _Tfreeze
        θstar = θ_r + (ν - θ_r) * inverse_matric_potential(hydraulics, ψ0 + ψT)
    else
        θstar = θ_l
    end
    Δz = source_type.Δz
    τLTE = FT(ρc_s * Δz^2 / κ)
    ΔT = norm(diffusive.soil.heat.κ∇T) / κ * Δz
    ρ_w = FT(0.5) * (_ρliq + _ρice)

    τpt = τLTE * (ρ_w * _LH_f0 * (ν - θ_r)) / (ρc_s * ΔT)

    τft = max(τLTE, τpt)
    freeze_thaw =
        FT(1) / τft * (
            _ρliq *
            (θ_l - θstar) *
            heaviside(_Tfreeze - T) *
            heaviside(θ_l - θstar) - _ρice * θ_i * heaviside(T - _Tfreeze)
        )
    return (; freeze_thaw)
end

function source(
    ::VolumetricLiquidFraction,
    s::PhaseChange,
    land::LandModel,
    args,
)
    @unpack state = args
    @unpack freeze_thaw = args.precomputed.dtup[s]
    FT = eltype(state)
    _ρliq = FT(ρ_cloud_liq(land.param_set))
    return -freeze_thaw / _ρliq
end

function source(::VolumetricIceFraction, s::PhaseChange, land::LandModel, args)
    @unpack state = args
    @unpack freeze_thaw = args.precomputed.dtup[s]
    FT = eltype(state)
    _ρice = FT(ρ_cloud_ice(land.param_set))
    return freeze_thaw / _ρice
end

ext: jl


text:
module ExampleFs
using ..FuseApi

using CStructures: CStruct, CStructGuarded, CVector, Cserialize

import Base: open, read, write, unlink, rename, readdir

using Base.Filesystem
import Base.Filesystem: S_IFDIR, S_IFREG, S_IFCHR, S_IFBLK, S_IFIFO, S_IFLNK, S_IFSOCK, S_IFMT
import Base.Filesystem: S_IRWXU, S_IRWXG, S_IRWXO
import Base.Filesystem: S_IRUSR, S_IRGRP, S_IROTH
import Base.Filesystem: S_IWUSR, S_IWGRP, S_IWOTH
import Base.Filesystem: S_IXUSR, S_IXGRP, S_IXOTH
import Base.Filesystem: S_ISUID, S_ISGID, S_ISVTX

import Base: UV_ENOENT, UV_ENOTDIR, UV_EEXIST, UV_ENOTEMPTY, UV_ENAMETOOLONG, UV_ERANGE
const UV_ENODATA = -61 # missing in Base

const X_UGO = S_IRWXU | S_IRWXG | S_IRWXO
const X_NOX = X_UGO & ~((S_IXUSR | S_IXGRP | S_IXOTH) & X_UGO)

const DIR1 = "."
const DIR2 = ".."

const PATH_MAX = 4096

const DEFAULT_MODE = S_IFREG | X_UGO

const T_ATIME = 0x0001
const T_MTIME = 0x0002
const T_CTIME = 0x0004

struct Direntry
    ino::FuseIno
    name::String
end

mutable struct Inode
    ino::FuseIno
    mode::FuseMode
    nlink::Int64
    size::UInt64
    uid::UInt32
    gid::UInt32
    atime::Timespec
    mtime::Timespec
    ctime::Timespec
    lookup::UInt64
    function Inode(ino, mode, nlink, size, uid, gid)
        t = timespec_now()
        new(ino, mode, nlink, size, uid, gid, t, t, t, 0)
    end
end

is_directory(inode::Inode) = inode.mode & S_IFMT == S_IFDIR
is_regular(inode::Inode) = inode.mode & S_IFMT == S_IFREG

struct FilesystemData
    inodes::Dict{FuseIno, Inode}
    data::Dict{FuseIno, Vector{UInt8}}
    data_dir::Dict{FuseIno, Vector{Direntry}}
    xattr::Dict{FuseIno, Vector{Tuple{String,Vector{UInt8}}}}
    FilesystemData(root) = new(Dict(FUSE_INO_ROOT => root), Dict(), Dict(), Dict())
end

reset_inodes(fs::FilesystemData) = for inode in values(fs.inodes); inode.lookup = 0; end
get_inode(fs::FilesystemData, ino::FuseIno) = get(fs.inodes, ino, nothing)
has_inode(fs::FilesystemData, ino::FuseIno) = haskey(fs.inodes, ino)
next_ino(fs::FilesystemData) = maximum(keys(fs.inodes)) + 1
push_inode!(fs::FilesystemData, inode::Pair{FuseIno,Inode}) = push!(fs.inodes, inode)

get_data(fs::FilesystemData, ino::FuseIno) = get(fs.data, ino) do; UInt8[] end
get_data!(fs::FilesystemData, ino::FuseIno) = get!(fs.data, ino) do; UInt8[] end
set_data!(fs::FilesystemData, ino::FuseIno, v) = setindex!(fs.data, v, ino)
has_data(fs::FilesystemData, ino::FuseIno) = haskey(fs.data, ino)
push_data!(fs::FilesystemData, inode::Pair{FuseIno}) = push!(fs.data, inode)
delete_data!(fs::FilesystemData, inode::FuseIno) = delete!(fs.data, inode)

function lookup_ino(fs::FilesystemData, parent::Integer, name::AbstractString)
    haskey(fs.data_dir, parent) || return 0
    direntries = fs.data_dir[parent]
    k = findfirst(de -> de.name == name, direntries)
    k === nothing ? 0 : direntries[k].ino
end
get_direntries(fs::FilesystemData, ino::Integer) = get(fs.data_dir, ino) do; Direntry[] end
get_direntries!(fs::FilesystemData, ino::Integer) = get!(fs.data_dir, ino) do; Direntry[] end

get_xattr(fs::FilesystemData, ino::Integer) = get(fs.xattr, ino) do; emptyxvec() end
set_xattr!(fs::FilesystemData, ino::Integer, vec) = fs.xattr[ino] = vec

function destroy_ino!(fs::FilesystemData, ino::FuseIno)
    delete!(fs.inodes, ino)
    delete!(fs.data, ino)
    delete!(fs.data_dir, ino)
    delete!(fs.xattr, ino)
    nothing
end

status(inode::Inode) = CStructGuarded{Cstat}(Cserialize(Cstat, inode))
status(fs::FilesystemData, ino::FuseIno) = status(get_inode(fs, ino))

function example(fs::FilesystemData, parent::Integer, name::String, data=nothing)
    inode = do_create(fs, FuseIno(parent), name, S_IFREG | X_UGO, 0, 0)
    if data !== nothing
        d =  convert(Vector{UInt8}, data)
        inode.size = length(d)
        push_data!(fs, inode.ino => d)
    end
end

# implementation of callback functions

export init
function init(fs::FilesystemData, conn::CStruct{FuseConnInfo})
    #println("init userdata = '$fs'")
    #println("init conninfo = $conn")
    reset_inodes(fs)
end

export destroy
function destroy(fs::FilesystemData)
    #println("destroy userdata = '$fs'")
end

export lookup
function lookup(fs::FilesystemData, req::FuseReq, parent::FuseIno, name::String)
    entry = do_lookup(fs, req, parent, name)
    entry isa Number && return entry
    inode = get_inode(fs, entry.ino)
    inode.lookup += 1
    fuse_reply_entry(req, entry)
end

function find_direntry(fs::FilesystemData, parent::Integer, name::AbstractString)
    parent = FuseIno(parent)
    has_inode(fs, parent) || return UV_ENOENT
    is_directory(fs.inodes[parent]) || return UV_ENOTDIR
    ino = lookup_ino(fs, parent, name)
    ino == 0 ? UV_ENOENT : ino
end

function do_lookup(fs::FilesystemData, req::FuseReq, parent::FuseIno, name::String)
    ino = find_direntry(fs, parent, name)
    ino <= 0 && return ino
    sta = status(fs, ino)
    entry = CStructGuarded(FuseEntryParam, (ino=ino, attr=sta))
    entry.attr_timeout = 0.0
    entry.entry_timeout = 1e19
    entry
end

export forget
function forget(fs::Any, req::FuseReq, ino::FuseIno, count::Integer)
    println("forget($ino, $count)")
    inode = get_inode(fs, ino)::Inode
    nl = inode.lookup
    nl = nl >= count ? nl - count : 0
    inode.lookup = nl
    if nl <= 0
        destroy_ino!(fs, ino)
    end
    fuse_reply_none(req)
end

export getattr
function getattr(fs::Any, req::FuseReq, ino::FuseIno, ::CStruct{FuseFileInfo})
    has_inode(fs, ino) || return UV_ENOENT
    attr = status(fs, ino)
    attr_timeout = 10.0
    fuse_reply_attr(req, attr, attr_timeout)
end

export setattr
function setattr(fs::Any, req::FuseReq, ino::FuseIno, st::CStruct{Cstat}, to_set::Integer, fi::CStruct{FuseFileInfo})
    has_inode(fs, ino) || return UV_ENOENT
    inode = get_inode(fs, ino)::Inode
    resetsuid = false
    if to_set & FUSE_SET_ATTR_MODE != 0
        inode.mode = st.mode
    end
    if to_set & FUSE_SET_ATTR_UID != 0
        resetsuid |= inode.uid != st.uid
        inode.uid = st.uid
    end
    if to_set & FUSE_SET_ATTR_GID != 0
        inode.gid = st.gid
    end
    if to_set & FUSE_SET_ATTR_SIZE != 0
        resetsuid |= inode.size != st.size
        inode.size = st.size
    end
    if to_set & FUSE_SET_ATTR_ATIME != 0
        inode.atime = st.atime
    end
    if to_set & FUSE_SET_ATTR_MTIME != 0
        inode.mtime = st.mtime
    end
    if to_set & (FUSE_SET_ATTR_ATIME_NOW | FUSE_SET_ATTR_MTIME_NOW) != 0
        t = timespec_now()
        if to_set & FUSE_SET_ATTR_ATIME_NOW != 0
            inode.atime = t
        end
        if to_set & FUSE_SET_ATTR_MTIME_NOW != 0
            inode.mtime = t
        end
    end
    if to_set & FUSE_SET_ATTR_CTIME != 0
        inode.ctime = st.ctime
    end
    if resetsuid
        inode.mode &= ~(S_ISUID | S_ISGID)
    end
    getattr(fs, req, ino, fi)
end

export opendir
function opendir(fs::Any, req::FuseReq, ino::FuseIno, fi::CStruct{FuseFileInfo})
    fuse_reply_open(req, fi)
end

export releasedir
function releasedir(fs::Any, req::FuseReq, ino::FuseIno, fi::CStruct{FuseFileInfo})
    fuse_reply_err(req, 0)
end

export readdir
function readdir(fs::Any, req::FuseReq, ino::FuseIno, size::Integer, off::Integer, ::CStruct{FuseFileInfo})
    buf = Vector{UInt8}(undef, size)
    dir = get_inode(fs, ino)::Inode
    is_directory(dir) || return UV_ENOTDIR
    p = 1
    des = get_direntries(fs, ino)
    n = length(des)
    rem = size
    i = off + 1
    while i <= n + 2
        if i <= 2
            st = status(dir)
            name = i == 1 ? DIR1 : DIR2
        else
            de = des[i-2]
            st = status(fs, de.ino)
            name = de.name
        end
        entsize = fuse_add_direntry(req, view(buf, p:size), name, st, i)
        if entsize > rem
            break
        end
        p += entsize
        rem -= entsize
        i += 1
    end
    if p > 1
        touch(dir, T_ATIME)
    end
    fuse_reply_buf(req, buf, size - rem)
end

function do_create(fs::FilesystemData, parent::FuseIno, name::String, mode::Integer, uid::Integer, gid::Integer)
    dir = get_inode(fs, parent)::Inode
    is_directory(dir) || return UV_ENOTDIR
    lookup_ino(fs, parent, name) == 0 || return UV_EEXIST
    ino = next_ino(fs)
    inode = Inode(ino, mode, 1, 0, uid, gid)
    push_inode!(fs, ino => inode)
    direntries = get_direntries!(fs, parent)
    push!(direntries, Direntry(ino, name))
    touchdir(dir, direntries, T_MTIME | T_CTIME)
    inode
end

export create
function create(fs::Any, req::FuseReq, parent::FuseIno, name::String, mode::FuseMode, fi::CStruct{FuseFileInfo})
    ctx = fuse_req_ctx(req)
    uid = ctx.uid
    gid = ctx.gid
    umask = ctx.umask
    mode = mode & ~umask | mode & ~X_UGO
    do_create(fs, parent, name, mode, uid, gid)
    entry = do_lookup(fs, req, parent, name)
    inode = get_inode(fs, entry.ino)
    inode.lookup += 1
    fuse_reply_create(req, entry, fi)
end

export open
function open(fs::Any, req::FuseReq, ino::FuseIno, fi::CStruct{FuseFileInfo})
    inode = get_inode(fs, ino)::Inode
    if is_regular(inode) && fi.flags & JL_O_TRUNC != 0
        inode.size = 0
        delete_data!(fs, ino)
    end
    fuse_reply_open(req, fi)
end

export read
function read(fs::Any, req::FuseReq, ino::FuseIno, size::Integer, off::Integer, ::CStruct{FuseFileInfo})
    has_inode(fs, ino) || return UV_ENOENT
    inode = get_inode(fs, ino)::Inode
    is_directory(inode) && return UV_ENOENT
    data = get_data(fs, ino)
    sz = max(min(length(data) - off, size), 0)
    bufv = CStructGuarded{FuseBufvec}(Cserialize(FuseBufvec, (count=1, buf=[(size=sz, mem = data)])))
    if sz > 0
        touch(inode, T_ATIME)
    end
    fuse_reply_data(req, bufv, FuseBufCopyFlags(0))
end

export write
function write(fs::Any, req::FuseReq, ino::FuseIno, cvec::CVector{UInt8}, size::Integer, off::Integer, ::CStruct{FuseFileInfo})
    inode = get_inode(fs, ino)::Inode
    is_directory(inode) && return UV_ENOENT
    data = get_data!(fs, ino)
    total = off + size
    if total > length(data)
        resize!(data, total)
        for i = length(data)+1:off
            data[i] = 0
        end
    end
    for i = 1:size
        data[off+i] = cvec[i]
    end
    inode.size = length(data)
    if size > 0
        touch(inode, T_MTIME | T_CTIME)
        inode.mode &= ~(S_ISUID | S_ISGID)
    end
    fuse_reply_write(req, size)
end

export link
function link(fs::Any, req::FuseReq, ino::FuseIno, newparent::FuseIno, name::String)
    inode = get_inode(fs, ino)::Inode
    is_directory(inode) && return UV_EPERM
    dir = get_inode(fs, newparent)::Inode
    is_directory(dir) || return UV_ENOTDIR
    lookup_ino(fs, newparent, name) == 0 || return UV_EEXIST
    direntries = get_direntries!(fs, newparent)
    push!(direntries, Direntry(ino, name))
    touchdir(dir, direntries, T_MTIME | T_CTIME)
    entry = do_lookup(fs, req, newparent, name)
    inode.nlink += 1
    inode.lookup += 1
    fuse_reply_entry(req, entry)
end

export unlink
function unlink(fs::Any, req::FuseReq, parent::FuseIno, name::String)
    ino = find_direntry(fs, parent, name)
    ino <= 0 && return ino
    direntries = get_direntries(fs, parent)
    deleteat!(direntries, findall(de -> de.name == name, direntries))
    dir = get_inode(fs, parent)::Inode
    touchdir(dir, direntries, T_MTIME | T_CTIME)
    inode = get_inode(fs, ino)::Inode
    touch(inode, T_CTIME)
    return fuse_reply_err(req, 0)
end

export rename
function rename(fs::Any, req::FuseReq, parent::FuseIno, name::String, newparent::FuseIno, newname::String, flags::Integer)
    # println("rename($parent/$name => $newparent/$newname)")
    if parent == newparent && name == newname
        return fuse_reply_err(req, 0)
    end
    dir = get_inode(fs, parent)::Inode
    dirnew = parent == newparent ? dir : get_inode(fs, newparent)::Inode
    ino = lookup_ino(fs, parent, name)
    ino == 0 && return UV_ENOENT
    newino = lookup_ino(fs, newparent, newname)
    newino == 0 && flags == RENAME_EXCHANGE && return UV_ENOENT
    newino != 0 && flags == RENAME_NOREPLACE && return UV_EEXIST
    des = get_direntries(fs, parent)
    desnew = parent == newparent ? des : get_direntries!(fs, newparent)
    k = findfirst(de -> de.name == name, des)::Integer
    de = des[k]
    knew = newino == 0 ? 0 : findfirst(de -> de.name == newname, desnew)::Integer

    if flags == RENAME_EXCHANGE
        den = desnew[knew]
        des[k] = Direntry(de.ino, den.name)
        desnew[knew] = Direntry(den.ino, de.name)
    else
        if parent == newparent
            des[k] = Direntry(de.ino, newname)
        else
            push!(desnew, name == newname ? des[k] : Direntry(de.ino, newname))
            deleteat!(des, k)
        end
        if knew != 0 # && flags != RENAME_NOREPLACE
            deleteat!(desnew, knew)
        end
    end

    touchdir(dir, des, T_MTIME | T_CTIME)
    if parent != newparent
        touchdir(dirnew, desnew, T_MTIME | T_CTIME)
    end
    fuse_reply_err(req, 0)
end

export mkdir
function mkdir(fs::Any, req::FuseReq, parent::FuseIno, name::String, mode::FuseMode)
    ctx = fuse_req_ctx(req)
    uid = ctx.uid
    gid = ctx.gid
    umask = ctx.umask
    mode = mode & X_UGO & ~umask | mode & (S_ISUID | S_ISGID | S_ISVTX) | S_IFDIR
    inode = do_create(fs, parent, name, mode, uid, gid)
    entry = do_lookup(fs, req, parent, name)
    inode.lookup += 1
    fuse_reply_entry(req, entry)
end

export rmdir
function rmdir(fs::Any, req::FuseReq, parent::FuseIno, name::String)
    ino = find_direntry(fs, parent, name)
    ino <= 0 && return ino
    !isempty(get_direntries(fs, ino)) && return UV_ENOTEMPTY
    direntries = get_direntries(fs, parent)
    deleteat!(direntries, findall(de -> de.name == name, direntries))
    dir = get_inode(fs, parent)::Inode
    touchdir(dir, direntries, T_MTIME | T_CTIME)
    return fuse_reply_err(req, 0)
end

export symlink
function symlink(fs::Any, req::FuseReq, link::String, parent::FuseIno, name::String)
    # println("symlink($parent/$name => $link)")
    n = length(link)
    0 < n < PATH_MAX || return UV_ENAMETOOLONG
    ctx = fuse_req_ctx(req)
    uid = ctx.uid
    gid = ctx.gid
    mode = X_UGO | S_IFLNK
    inode = do_create(fs, parent, name, mode, uid, gid)
    inode.size = n
    vlink = collect(codeunits(link))
    set_data!(fs, inode.ino, push!(vlink, 0x00))
    entry = do_lookup(fs, req, parent, name)
    inode.lookup += 1
    fuse_reply_entry(req, entry)
end

export readlink
function readlink(fs::Any, req::FuseReq, ino::FuseIno)
    link = get_data(fs, ino)
    fuse_reply_readlink(req, link)
end

export setxattr
function setxattr(fs::Any, req::FuseReq, ino::FuseIno, name::String, value::Vector{UInt8}, flags::Integer)
    vec = get_xattr(fs, ino)
    k = findfirst(x->x[1] == name, vec)
    hk = k === nothing
    flags == XATTR_CREATE && hk && return UV_EEXIST
    flags == XATTR_REPLACE && !hk && return UV_ENODATA
    if k === nothing
        set_xattr!(fs, ino, push!(vec, (name, value)))
    else
        vec[k] = (name, value)
    end
    fuse_reply_err(req, 0)
end

export getxattr
function getxattr(fs::Any, req::FuseReq, ino::FuseIno, name::String, size::Integer)
    vec = get_xattr(fs, ino)
    k = findfirst(x->x[1] == name, vec)
    k === nothing && return UV_ENODATA
    val = vec[k][2]
    #println("getxattr($name) = $val")
    bsize = sizeof(val)
    if size == 0
        fuse_reply_xattr(req, bsize)
    else
        bsize > size && return UV_ERANGE
        fuse_reply_buf(req, val, bsize)
    end
end

export removexattr
function removexattr(fs::Any, req::FuseReq, ino::FuseIno, name::String)
    vec = get_xattr(fs, ino)
    k = findfirst(x->x[1] == name, vec)
    if k !== nothing
        deleteat!(vec, k)
    end
    fuse_reply_err(req, 0)
end

export listxattr
function listxattr(fs::Any, req::FuseReq, ino::FuseIno, size::Integer)
    str = join([t[1] for t in get_xattr(fs, ino)], '\0')
    val = collect(codeunits(str))
    !isempty(val) && push!(val, 0x0)
    #println("listxattr($ino, $size); val = $val")
    bsize = sizeof(val)
    if size == 0
        fuse_reply_xattr(req, bsize)
    else
        bsize > size && return UV_ERANGE
        fuse_reply_buf(req, val, bsize)
    end
end

# utility functions

function touchdir(dir::Inode, direntries, mode)
    dir.size = length(direntries)
    touch(dir, mode)
end
emptyxvec() = Tuple{String,Vector{UInt8}}[]

function Base.touch(inode::Inode, mode::UInt16)
    t = timespec_now()
    if mode & T_ATIME != 0
        inode.atime = t
    end
    if mode & T_MTIME != 0
        inode.mtime = t
    end
    if mode & T_CTIME != 0
        inode.ctime = t
    end
    t
end

args = ["mountpoint", "-s", "-d", "-h", "-V", "-f" , "-o", "clone_fd", "-o", "max_idle_threads=5"]

const FS = FilesystemData(Inode(FUSE_INO_ROOT, S_IFDIR | X_UGO, 1, 0, 0, 0))
example(FS, 1, "xxx", codeunits("hello world 4711!\n"))

using Base.Threads

bg() = @spawn main_loop($args, @__MODULE__, FS)
fg() = main_loop(args, @__MODULE__, FS)

end # module ExampleFs


text:
# boundary condition functions

@doc """
Now actually we are integrating
  ∫ G∇ϕ:[q] dΓ

Input:
  mesh
  sbp
  eqn
  opts
Output:

"""->
function evalBoundaryIntegrals_vector(mesh::AbstractMesh{Tmsh},
                                      sbp::AbstractOperator,
                                      eqn::NSData{Tsol, Tres, Tdim},
                                      opts) where {Tmsh, Tsol, Tres, Tdim}

  sbpface = mesh.sbpface
  Dx = Array{Tmsh}((mesh.numNodesPerElement, mesh.numNodesPerElement, Tdim))
  R = sview(sbpface.interp[:,:])
  w = sview(sbpface.wface, :)
  res = sview(eqn.res, :,:,:)

  numNodes_elem = mesh.numNodesPerElement
  numNodes_face = mesh.numNodesPerFace
  stencilsize   = sbpface.stencilsize
  q_bnd = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  dq    = Array{Tsol}(mesh.numDofPerNode, mesh.numNodesPerFace)
  Gt = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  RDx = zeros(Tmsh, mesh.numNodesPerFace, mesh.numNodesPerElement, Tdim)
  GtRDx = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace, mesh.numNodesPerElement)
  nrm     = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  nrm1 = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  area = Array{Tmsh}(mesh.numNodesPerFace)

  # loop over all the boundaries
  for bc = 1:mesh.numBC
    indx0 = mesh.bndry_offsets[bc]
    indx1 = mesh.bndry_offsets[bc+1] - 1

    for f = indx0:indx1
      bndry = mesh.bndryfaces[f]
      elem = bndry.element
      face = bndry.face
      p = sview(sbpface.perm, :, face)

      # compute RDx
      calcDx(mesh, sbp, elem, Dx)

      for i = 1 : length(RDx)
        RDx[i] = 0.0
      end

      for d =  1 : Tdim
        for row = 1 : numNodes_face
          for col = 1 : numNodes_elem
            for s = 1 : stencilsize
              RDx[row, col, d] +=  R[s, row] * Dx[p[s], col, d]
            end
          end
        end
      end

      vecflux = sview(eqn.vecflux_bndry, :,:,:,f)
      for i = 1 : numNodes_elem
        for j = 1 : numNodes_face
          for iDof = 2 : mesh.numDofPerNode
            # res[iDof, i, elem] +=  ( RDx[j, i, 1] * vecflux[1, iDof, j]
            # + RDx[j, i, 2] * vecflux[2, iDof, j] ) * w[j]
            tmp = 0.0
            for iDim = 1 : Tdim
              tmp += RDx[j, i, iDim] * vecflux[iDim, iDof, j]
            end
            res[iDof, i, elem] += tmp * w[j]
          end
        end
      end
    end
  end

  return nothing
end  # end evalBoundaryIntegrals_vector

function calcViscousFlux_boundary(mesh::AbstractMesh{Tmsh},
                                  sbp::AbstractOperator,
                                  eqn::NSData{Tsol, Tres, Tdim},
                                  opts) where {Tmsh, Tsol, Tres, Tdim}
  # freestream info
  Ma = eqn.params.euler_params.Ma
  Re = eqn.params.Re
  gamma_1 = eqn.params.euler_params.gamma_1
  Pr = 0.72
  coef_nondim = Ma/Re

  p = opts["order"]
  sat_type = opts["SAT_type"]
  const_tii = (p + 1.0)*(p + Tdim)/Tdim
  sbpface = mesh.sbpface
  dq    = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  dqn   = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  q_bnd = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
  pMat  = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, mesh.numNodesPerFace)
  Gt    = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  Gt_bnd  = zeros(Tsol, mesh.numDofPerNode, mesh.numDofPerNode, Tdim, Tdim, mesh.numNodesPerFace)
  Fv_face = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  Fv_bnd  = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)
  vecflux = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace)

  nrm1 = Array{Tmsh}(Tdim, mesh.numNodesPerFace)
  area = Array{Tmsh}(mesh.numNodesPerFace)
  area_sum = sview(eqn.area_sum, :)

  # sigma = calcTraceInverseInequalityConst(sbp, sbpface)
  dqdx_elem = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerElement )
  dqdx_face = zeros(Tsol, Tdim, mesh.numDofPerNode, mesh.numNodesPerFace )
  for iBC = 1 : mesh.numBC
    indx0 = mesh.bndry_offsets[iBC]
    indx1 = mesh.bndry_offsets[iBC+1] - 1

    # specify boundary value function
    # TODO: Put it into a function
    local bnd_functor::AbstractBoundaryValueType
    key_i = string("BC", iBC, "_name")
    val = opts[key_i]
    Gt_functor = calcDiffusionTensor
    if val == "FreeStreamBC"
      bnd_functor = Farfield()
    elseif val == "ExactChannelBC"
      bnd_functor = ExactChannel()
    elseif val == "nonslipBC"
      bnd_functor = AdiabaticWall()
      Gt_functor = calcDiffusionTensorOnAdiabaticWall
    elseif val == "noPenetrationBC"
      continue
    elseif val == "zeroPressGradientBC"
      bnd_functor = Farfield()
    else
      error("iBC = ", iBC, ", Only 'FreeStreamBC' and 'nonslipBC' available")
    end

    for f = indx0 : indx1

      flux  = zeros(Tsol, mesh.numDofPerNode, mesh.numNodesPerFace)
      bndry = mesh.bndryfaces[f]
      elem = bndry.element
      face = bndry.face
      perm = sview(sbpface.perm, :, face)
      xy = sview(mesh.coords_bndry, :, :, f)

      # Compute geometric info on face
      nrm_xy = ro_sview(mesh.nrm_bndry, :, :, f)
      for n = 1 : mesh.numNodesPerFace
        area[n] = norm(ro_sview(nrm_xy, :, n))

        for i = 1 : Tdim
          nrm1[i,n] = nrm_xy[i,n] / area[n]
        end
      end

      # We need viscous flux and diffusion tensor on interfaces, and there
      # are different ways to compute them. For viscous flux:
      # 1) since we have Q on face nodes, if 𝛻Q is available on interface, then we are done.
      # 2) we can comoute viscous flux on volume nodes and then interpolate to interface node.
      # It's logically simple but computationally expensive.

      # Compute boundary viscous flux, F(q_b, ∇q) = G(q_b)∇q.
      # so we need viscousity tensor G, and derivative of q.
      q_face = sview(eqn.q_bndry, :, :, f)
      bnd_functor(q_face, xy, nrm1, eqn.params, q_bnd)

      # diffusion matrix used in penalty term should be computed from q_face rather than q_bnd
      if val == "nonslipBC"
        Gt_functor(eqn.params, q_bnd, nrm1, Gt)
      else
        Gt_functor(eqn.params, q_bnd, Gt)
      end
      q_elem = sview(eqn.q, :, :, elem)
      calcGradient(mesh, sbp, elem, q_elem, dqdx_elem)

      #
      # TODO: we can consider the first 2 dimension as a single dimension,
      # then we will not need slice here any more.
      #
      for d = 1 : Tdim
        q_x_node = sview(dqdx_elem, d, :, :)
        q_x_face = sview(dqdx_face, d, :, :)
        boundaryinterpolate(sbpface, bndry, q_x_node, q_x_face)
      end

      calcFvis(eqn.params, Gt, dqdx_face, Fv_face)

      # compute penalty matrix
      cmptBPMat(mesh, sbp, eqn, opts, f, Gt, pMat)

      # Start to compute fluxes.  We have 3 terms on interfaces:
      # 1) -{Fv}⋅[ϕ]
      # 2) -{G^T ∇ϕ}⋅[q]
      # 3) +δ{G}[q]:[ϕ]
      for n = 1 : mesh.numNodesPerFace
        for iDof = 1 : mesh.numDofPerNode
          dq[iDof, n] = q_face[iDof, n] - q_bnd[iDof, n]
          # dqn[1, iDof, n] = -dq[iDof, n]*nrm_xy[1, n]
          # dqn[2, iDof, n] = -dq[iDof, n]*nrm_xy[2, n]
        end
      end
      # DEBUG BEGIN
      # if maximum(abs.(real(dq))) > 1.e-11
        # println(real(dq))
      # end
      # DEBUG END

      # -----------------------------------------------
      # This part computes the contribution of
      # ∫ {G^T∇ϕ}:[q] dΓ = ∫ ∇ϕ⋅F dΓ ,
      # where
      # [q] = (q+ - q-) ⊗ n,
      # G = G(q_b) depends on boudanry value.
      # Then we can consider Δq⊗n as ∇q and F as viscous flux.
      # -----------------------------------------------

      # calcFvis(Gt, dqn, vecflux)
      fill!(vecflux, 0.0)
      for n = 1 : mesh.numNodesPerFace
        for iDof = 1 : mesh.numDofPerNode
          for iDim = 1 : Tdim
            for jDim = 1 : Tdim
              tmp = 0.0
              for jDof = 1 : mesh.numDofPerNode
                tmp += Gt[iDof, jDof, iDim, jDim, n]
              end
              vecflux[iDim, iDof, n] += tmp * nrm_xy[jDim, n]
            end
            vecflux[iDim,iDof,n] *=  dq[iDof,n]
          end
        end
      end

      for n = 1 : mesh.numNodesPerFace
        for iDof = 2 : mesh.numDofPerNode
          for iDim = 1 : Tdim
            flux[iDof, n] -= Fv_face[iDim, iDof, n]*nrm_xy[iDim,n]
          end
        end
      end

      for n = 1 : mesh.numNodesPerFace
        for iDof = 2 : mesh.numDofPerNode
          for jDof = 1: mesh.numDofPerNode
            flux[iDof, n] +=  pMat[iDof, jDof, n]*dq[jDof, n]
          end
        end
      end

      # accumulate fluxes
      for n = 1 : mesh.numNodesPerFace
        for iDof = 2 : Tdim+2
          for iDim = 1 : Tdim
            eqn.vecflux_bndry[iDim, iDof, n, f] -=  vecflux[iDim, iDof, n]*coef_nondim
          end
          eqn.bndryflux[iDof, n, f] += flux[iDof, n]*coef_nondim
        end
      end
    end # loop over faces of one BC
  end # loop over BCs
  return nothing
end

#------------------------------------------------------------------------------
# functions for computing qR

abstract type AbstractBoundaryValueType end

@doc """

compute the adiabatic wall boundary value

Input:
q_in    :: conservative variables on boundary face
norm    :: outward unit normal of the boundary face
Output:
q_bnd    :: the boundary value

"""->
mutable struct AdiabaticWall <: AbstractBoundaryValueType
end
function (obj::AdiabaticWall)(
              q_in::AbstractArray{Tsol, 2},
              xy::AbstractArray{Tmsh, 2},
              norm::AbstractArray{Tmsh, 2},
              params::ParamType{Tdim},
              q_bnd::AbstractArray{Tsol, 2}) where {Tsol, Tmsh, Tdim}
  @assert( size(q_in, 1) == size(q_bnd,  1))
  @assert( size(q_in, 2) == size(q_bnd,  2))
  @assert( size(q_in, 2) == size(norm,  2))

  dim = size(norm, 1)
  numNodes = size(q_in, 2)

  for n = 1 : numNodes
    q_bnd[1, n] = q_in[1, n]
    # noslip condition
    q_bnd[2:dim+1, n] = 0.0
    q_bnd[dim+2, n] = q_in[dim+2, n]
  end

  return nothing
end

mutable struct ExactChannel<: AbstractBoundaryValueType
end
function (obj::ExactChannel)(
                          q_in::AbstractArray{Tsol, 2},
                          xyz::AbstractArray{Tmsh, 2},
                          norm::AbstractArray{Tmsh, 2},
                          params::ParamType{3},
                          qg::AbstractArray{Tsol, 2}) where {Tsol, Tmsh}

  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1
  sigma = 0.01
  aoa = params.euler_params.aoa
  beta = params.sideslip_angle
  rhoInf = 1.0
  uInf = params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = params.euler_params.Ma * sin(beta) * -1
  wInf = params.euler_params.Ma * cos(beta) * sin(aoa)
  TInf = 1.0
  numNodes = size(qg, 2)
  for n = 1 : numNodes
    x = xyz[1,n]
    y = xyz[2,n]
    z = xyz[3,n]

    rho = rhoInf * (1 + sigma*x*y*z)
    ux = sin(pi*x) + 1
    uy = sin(pi*y) + 1
    uz = sin(pi*z) + 1
    u  = (1 + sigma*ux * uy * uz )* uInf
    vx = sin(pi*x) + 1
    vy = sin(pi*y) + 1
    vz = sin(pi*z) + 1
    v  = (1 + sigma*vx * vy * vz )* vInf
    wx = sin(pi*x) + 1
    wy = sin(pi*y) + 1
    wz = sin(pi*z) + 1
    w  = (1 + sigma*wx * wy * wz) * wInf
    T  = TInf

    if !params.isViscous
      u += 0.2 * uInf
      v += 0.2 * vInf
      w += 0.2 * wInf
    end

    qg[1,n] = rho
    qg[2,n] = rho*u
    qg[3,n] = rho*v
    qg[4,n] = rho*w
    qg[5,n] = T/(gamma * gamma_1) + 0.5 * (u*u + v*v + w*w)
    qg[5,n] *= rho
  end

  return nothing
end

function (obj::ExactChannel)(
                          q_in::AbstractArray{Tsol, 2},
                          xy::AbstractArray{Tmsh, 2},
                          norm::AbstractArray{Tmsh, 2},
                          params::ParamType{2},
                          q_bnd::AbstractArray{Tsol, 2}) where {Tsol, Tmsh}
  @assert( size(q_in, 1) == size(q_bnd,  1))
  @assert( size(q_in, 2) == size(q_bnd,  2))
  @assert( size(q_in, 2) == size(norm,  2))
  dim = size(norm, 1)
  numNodes = size(q_in, 2)
  sigma = 0.1

  gamma = params.euler_params.gamma
  gamma_1 = gamma - 1

  aoa = params.euler_params.aoa
  qRef = zeros(Float64, dim+2)
  qRef[1] = 1.0
  qRef[2] = params.euler_params.Ma*cos(aoa)
  qRef[3] = params.euler_params.Ma*sin(aoa)
  qRef[4] = 1.0

  for n = 1 : numNodes
    x = xy[1, n]
    y = xy[2, n]
    rho = qRef[1] * (sigma*exp(sin(0.5*pi*(x+y))) +  1.0)
    ux  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[2]
    uy  = exp(y) * sin(pi*y)
    u   = ux * uy
    vx  = (exp(x) * sin(pi*x) * sigma + 1) * qRef[3]
    vy  = exp(y) * sin(pi*y)
    v   = vx * vy
    T   = (1 + sigma*exp(0.1*x+0.1*y)) * qRef[4]
    # T   = qRef[4]
    if !params.isViscous
      u += 0.2 * qRef[2]
    end
    q_bnd[1, n] = rho
    q_bnd[2, n] = rho * u
    q_bnd[3, n] = rho * v
    q_bnd[4, n] = T/(gamma*gamma_1) + 0.5 * (u*u + v*v)
    q_bnd[4, n] *= rho
  end

  return nothing
end

@doc """

compute the farfield boundary value

Input:
q_in    :: conservative variables on boundary face
norm    :: outward unit normal of the boundary face
Output:
q_bnd    :: the boundary value

"""->
mutable struct Farfield <: AbstractBoundaryValueType
end

function (obj::Farfield)(
              q_in::AbstractArray{Tsol, 2},
              xy::AbstractArray{Tmsh, 2},
              norm::AbstractArray{Tmsh, 2},
              params::ParamType{2},
              q_bnd::AbstractArray{Tsol, 2}) where {Tsol, Tmsh}
  @assert( size(q_in, 1) == size(q_bnd,  1))
  @assert( size(q_in, 2) == size(q_bnd,  2))
  @assert( size(q_in, 2) == size(norm,  2))

  dim      = size(norm, 1)
  numNodes = size(q_in, 2)
  lambda   = zeros(Float64, 3)
  gamma    = params.euler_params.gamma
  gamma_1  = params.euler_params.gamma - 1
  aoa      = params.euler_params.aoa
  gg_1     = (gamma*gamma_1)
  MaInf    = params.euler_params.Ma

  # freestream values
  qInf = zeros(Float64, dim + 2)
  qInf[1] = 1.0
  qInf[2] = qInf[1]*MaInf*cos(aoa)
  qInf[3] = qInf[1]*MaInf*sin(aoa)
  qInf[4] = qInf[1]*(1.0/gg_1 + 0.5*MaInf*MaInf)

  for n = 1 : numNodes
    #
    # contravariant velocity
    #
    u = q_in[2,n] / q_in[1,n]
    v = q_in[3,n] / q_in[1,n]
    vn = u*norm[1, n] + v*norm[2, n]
    v2 = u*u + v*v
    T = gg_1*(q_in[4,n]/q_in[1,n] - 0.5*v2)
    a = sqrt(T)
    #
    # eigenvalues
    #
    lambda[1] = real(vn)  # TODO: taking the real part is not complex safe!
    lambda[2] = real(vn + a)
    lambda[3] = real(vn - a)

    if lambda[2] <= 0.0            # supersonic inflow
      q_bnd[:, n] = qInf[:]
    elseif lambda[3] >= 0.0        # supersonic outflow
      q_bnd[:, n] = q_in[:, n]
    elseif lambda[1] <= 0.0        # subsonic inflow
      p = q_in[1,n]*T/gamma
      q_bnd[1, n] = qInf[1]
      q_bnd[2, n] = qInf[2]
      q_bnd[3, n] = qInf[3]
      q_bnd[4, n] = p/gamma_1 + 0.5*MaInf*MaInf*qInf[1]
    else                        # subsonic outflow
      pInf = 1.0/gamma
      q_bnd[1, n] = q_in[1, n]
      q_bnd[2, n] = q_in[2, n]
      q_bnd[3, n] = q_in[3, n]
      q_bnd[4, n] = pInf/gamma_1 + 0.5*q_in[1,n] * v2
    end

    # DEBUG ONLY
    # q_bnd[:,n] = qInf[:]
    # DEBUG END
  end

  return nothing
end

#
# TODO: combine this with the 2D version. The only difference
# between 2D and 3D is the freestream conditions.
#
function (obj::Farfield)(
              q_in::AbstractArray{Tsol, 2},
              xy::AbstractArray{Tmsh, 2},
              norm::AbstractArray{Tmsh, 2},
              params::ParamType{3},
              q_bnd::AbstractArray{Tsol, 2}) where {Tsol, Tmsh}
  @assert( size(q_in, 1) == size(q_bnd,  1))
  @assert( size(q_in, 2) == size(q_bnd,  2))
  @assert( size(q_in, 2) == size(norm,  2))
  @assert( size(norm, 1) == 3 )

  dim      = size(norm, 1)
  numNodes = size(q_in, 2)
  lambda   = zeros(Float64, 3)
  gamma    = params.euler_params.gamma
  gamma_1  = params.euler_params.gamma_1
  aoa      = params.euler_params.aoa
  beta     = params.sideslip_angle
  gg_1     = (gamma*gamma_1)
  MaInf    = params.euler_params.Ma

  uInf = params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = params.euler_params.Ma * sin(beta) * -1
  wInf = params.euler_params.Ma * cos(beta) * sin(aoa)
  # freestream values
  qInf = zeros(Float64, dim + 2)
  qInf[1] = 1.0
  qInf[2] = qInf[1]*uInf
  qInf[3] = qInf[1]*vInf
  qInf[4] = qInf[1]*wInf
  qInf[5] = qInf[1]*(1.0/gg_1 + 0.5*MaInf*MaInf)

  for n = 1 : numNodes
    #
    # contravariant velocity
    #
    vn = q_in[2, n]*norm[1, n] + q_in[3, n]*norm[2, n] + q_in[4,n] *norm[3,n]
    vn = vn/q_in[1,n]
    v2 = q_in[2,n]*q_in[2,n] + q_in[3,n]*q_in[3,n] + q_in[4,n]*q_in[4,n]
    v2 /= q_in[1,n]*q_in[1,n]
    T = gg_1*(q_in[5,n]/q_in[1,n] - 0.5*v2)
    a = sqrt(T)
    #
    # eigenvalues
    #
    lambda[1] = real(vn)
    lambda[2] = real(vn + a)
    lambda[3] = real(vn - a)

    # depending on eigenvalue, set state boundary value to exterior or interior (depending on field)
    if lambda[2] <= 0.0            # supersonic inflow
      q_bnd[:, n] = qInf[:]
    elseif lambda[3] >= 0.0        # supersonic outflow
      q_bnd[:, n] = q_in[:, n]
    elseif lambda[1] <= 0.0        # subsonic inflow
      p = q_in[1,n]*T/gamma
      q_bnd[1:dim+1, n] = qInf[1:dim+1]
      q_bnd[dim+2, n] = p/gamma_1 + 0.5*MaInf*MaInf*qInf[1]
    else                        # subsonic outflow
      pInf = 1.0/gamma
      q_bnd[1:dim+1, n] = q_in[1:dim+1, n]
      rhoV2 = q_in[1, n] * v2
      q_bnd[dim+2, n] = pInf/gamma_1 + 0.5*rhoV2
    end

    # debug only
    q_bnd[:,n] = qInf[:]        # for MMS in which q at boundary equals to qInf
    # end
  end

  return nothing
end

#------------------------------------------------------------------------------
# Boundary condition functors
mutable struct nonslipBC <: BCType
end
# low level function
function (obj::nonslipBC)(
              params::ParamType,
              q::AbstractArray{Tsol,1},
              aux_vars::AbstractArray{Tres, 1},
              x::AbstractArray{Tmsh,1},
              nrm_xy::AbstractArray{Tmsh,1},
              bndryflux::AbstractArray{Tres, 1}) where {Tmsh, Tsol, Tres}

  dim = length(nrm_xy)
	qg = params.qg
  # adiabatic wall
	qg[1] = q[1]
	qg[2:dim+1] = 0.0
	qg[dim+2] = q[dim+2]
  # isothermal wall
	# qg[1] = q[1]
	# rhoV2 = (q[2]*q[2] + q[3]*q[3])/q[1]
	# qg[2:dim+1] = 0.0
	# qg[dim+2] = q[4] - 0.5*rhoV2

	v_vals = params.v_vals
	convertFromNaturalToWorkingVars(params, qg, v_vals)
	# this is a problem: q is in conservative variables even if
	# params says we are using entropy variables
	calcEulerFlux(params, v_vals, aux_vars, nrm_xy, bndryflux)

	return nothing
end

mutable struct ExactChannelBC <: BCType
end
# low level function
function (obj::ExactChannelBC)(
              params::ParamType{3},
              q::AbstractArray{Tsol,1},
              aux_vars::AbstractArray{Tres, 1},
              xyz::AbstractArray{Tmsh,1},
              nrm_xy::AbstractArray{Tmsh,1},
              bndryflux::AbstractArray{Tres, 1}) where {Tmsh, Tsol, Tres}

  sigma = 0.01
  gamma = params.euler_params.gamma
  gamma_1 = params.euler_params.gamma - 1
  aoa = params.euler_params.aoa
  beta = params.sideslip_angle
  rhoInf = 1.0
  uInf = params.euler_params.Ma * cos(beta) * cos(aoa)
  vInf = params.euler_params.Ma * sin(beta) * -1
  wInf = params.euler_params.Ma * cos(beta) * sin(aoa)
  TInf = 1.0
  x = xyz[1]
  y = xyz[2]
  z = xyz[3]

  rho = rhoInf * (1 + sigma*x*y*z)
  ux = sin(pi*x) + 1
  uy = sin(pi*y) + 1
  uz = sin(pi*z) + 1
  u  = (1 + sigma*ux * uy * uz )* uInf
  vx = sin(pi*x) + 1
  vy = sin(pi*y) + 1
  vz = sin(pi*z) + 1
  v  = (1 + sigma*vx * vy * vz )* vInf
  wx = sin(pi*x) + 1
  wy = sin(pi*y) + 1
  wz = sin(pi*z) + 1
  w  = (1 + sigma*wx * wy * wz) * wInf
  T  = TInf

  if !params.isViscous
    u += 0.2 * uInf
    v += 0.2 * vInf
    w += 0.2 * wInf
  end

  qg = Array{Tsol}(5)
	qg[1] = rho
	qg[2] = rho*u
	qg[3] = rho*v
	qg[4] = rho*w
  qg[5] = T/(gamma * gamma_1) + 0.5 * (u*u + v*v + w*w)
  qg[5] *= rho

	v_vals = params.v_vals
	convertFromNaturalToWorkingVars(params, qg, v_vals)
	# this is a problem: q is in conservative variables even if
	# params says we are using entropy variables
	# calcEulerFlux(params, v_vals, aux_vars, [nx2, ny2], bndryflux)
  RoeSolver(params, q, qg, aux_vars, nrm_xy, bndryflux)

	return nothing
end

# low level function
function (obj::ExactChannelBC)(
              params::ParamType{2},
              q::AbstractArray{Tsol,1},
              aux_vars::AbstractArray{Tres, 1},
              x::AbstractArray{Tmsh,1},
              nrm_xy::AbstractArray{Tmsh,1},
              bndryflux::AbstractArray{Tres, 1}) where {Tmsh, Tsol, Tres}

  # functor ExactChannel takes varibales on multiple nodes, so we need to reshape some variables
  xy = reshape(x, length(x), 1)
  norm = reshape(nrm_xy, length(nrm_xy), 1)
  q_in = reshape(q, length(q), 1)
  q_bnd = zeros(Tsol, 4, 1)
  bnd_functor = ExactChannel()
  bnd_functor(q_in, xy, norm, params, q_bnd)
  qg = reshape(q_bnd, length(q_bnd))

	v_vals = params.v_vals
	convertFromNaturalToWorkingVars(params, qg, v_vals)
  RoeSolver(params, q, qg, aux_vars, nrm_xy, bndryflux)

	return nothing
end

mutable struct zeroPressGradientBC <: BCType
end

# low level function
function (obj::zeroPressGradientBC)(
                                params::ParamType,
                                q::AbstractArray{Tsol,1},
                                aux_vars::AbstractArray{Tres, 1},
                                x::AbstractArray{Tmsh,1},
                                nrm_xy::AbstractArray{Tmsh,1},
                                bndryflux::AbstractArray{Tres, 1}) where {Tmsh, Tsol, Tres}

  dim = length(nrm_xy)

	gamma = params.euler_params.gamma
	gamma_1 = params.euler_params.gamma_1
	qg = params.qg
	dim = 2
  rhoV2 = (norm(view(q, 2:dim+1))) / q[1]
	# rhoV2 = (q[2]*q[2] + q[3]*q[3]) / q[1]
	pinf = 1./gamma
	qg[1:dim+1] = q[1:dim+1]
	qg[dim+2] = pinf/gamma_1 + 0.5*rhoV2

	v_vals = params.v_vals
	convertFromNaturalToWorkingVars(params, qg, v_vals)
	# this is a problem: q is in conservative variables even if
	# params says we are using entropy variables
	calcEulerFlux(params, v_vals, aux_vars, nrm_xy, bndryflux)

	return nothing
end

global const BCDict = Dict{String, BCType}(
  "nonslipBC" => nonslipBC(),
  "ExactChannelBC" => ExactChannelBC(),
  "zeroPressGradientBC" => zeroPressGradientBC(),
)

"""
### NavierStokesnMod.getBCFunctors

  This function uses the opts dictionary to populate mesh.bndry_funcs with
  the functors

    func(params::ParamType,
         q::AbstractArray{Tsol,1},
         aux_vars::AbstractArray{Tres, 1},  coords::AbstractArray{Tmsh,1},
         nrm_xy::AbstractArray{Tmsh,1},
         bndryflux::AbstractArray{Tres, 1},
         bndry::BoundaryNode=NullBoundaryNode)

  This is a high level function.
"""
# use this function to populate access the needed values in BCDict
function getBCFunctors(mesh::AbstractMesh, sbp::AbstractOperator, eqn::NSData, opts)

  # the boundary conditions should switch over to the same model as the
  # other physics modules, but for now `calcViscousFlux_boundary` internally
  # picks the right functor
  # this will likely require each `eqn` object having its own `bndry_funcs`
  # array, rather than the mesh having a single one for all physics
  # In the meantime, `calcViscousFlux_boundary` is getting its functors
  # internally, and the invisicd functors are stored in the mesh.
  #
  # As of this writing, I believe the only boundary condition that works is
  # the FreeStreamBC, because it exists for both viscous and inviscid.
  # Eventually the inviscid part of every viscous boundary condition will
  # need to exist and both boundary conditions applied (related to the
  # comment above about each eqn object having its own BC functor array.
  error("getBCFunctors is not (currently) used for $PhysicsName")
  for i=1:mesh.numBC
    key_i = string("BC", i, "_name")
    val = opts[key_i]
    mesh.bndry_funcs[i] = BCDict[val]
  end

  return nothing
end # end function getBCFunctors


text:
module Matter

include("./geometry/vector.jl")
include("./geometry/vertex.jl")
export Vector, add!, sub!, mult!, div!, magnitude, normalize!, dot, det, cross, rotate!
export Vertex, Contact, vertices, center, mean, area, inertia, translate!, rotate!

end # module


text:
@testset "Dual Type" begin
    h = LPHRepresentation(spzeros(Int, 2, 2), [1, 2], [3, 4], [4, 5], [6, 7])
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Rational{BigInt}, LPHRepresentation{Rational{BigInt}, SparseMatrixCSC{Rational{BigInt},Int}}, Polyhedra.Hull{Rational{BigInt}, Vector{Rational{BigInt}}, Int}}
    h = hrep(zeros(2, 2), zeros(2))
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}
    h = hrep(spzeros(2, 2), zeros(2))
    p = polyhedron(h)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, SparseMatrixCSC{Float64, Int}}, MixedMatVRep{Float64, Matrix{Float64}}}
    v = vrep(zeros(2, 3))
    p = polyhedron(v)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, Matrix{Float64}}}
    v = vrep(spzeros(2, 3))
    p = polyhedron(v)
    @test p isa Polyhedra.DefaultPolyhedron{Float64, MixedMatHRep{Float64, Matrix{Float64}}, MixedMatVRep{Float64, SparseMatrixCSC{Float64, Int}}}
end

struct BadPoly{T} <: Polyhedra.Polyhedron{T}
end

@testset "Unimplemented methods" begin
    p = BadPoly{Int}()
    @test_throws ErrorException detectvlinearity!(p)
    @test_throws ErrorException detecthlinearity!(p)
    @test_throws ErrorException eliminate(p, [1], FourierMotzkin())
    @test_throws ErrorException eliminate(p, [1], BlockElimination())
    @test_throws ErrorException hrep(p)
    @test_throws ErrorException vrep(p)
    @test_throws ErrorException loadpolyhedron!(p, "a", "b")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ine")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ext")
end

@testset "Polyhedra.DefaultPolyhedron constructor with nothing" begin
    vr = convexhull([-1, 0], [0, -1]) + conichull([1, 1], [-1, 1])
    hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
    p = Polyhedra.DefaultPolyhedron{Int, typeof(hr), typeof(vr)}(hr, nothing, lp_solver)
    @test hrep(p) === hr
    @test !vrepiscomputed(p)
    @test p.solver === lp_solver
    p = Polyhedra.DefaultPolyhedron{Int, typeof(hr), typeof(vr)}(nothing, vr, lp_solver)
    @test !hrepiscomputed(p)
    @test vrep(p) === vr
    @test p.solver === lp_solver
end


text:
srand(1)

n = 3
k = 1
a = 2
# [Locus 1, Locus 2, Locus 3]
links = [0b101, 0b011, 0b110]
contribs = NK.AlleleContribs([
  Dict{NK.AlleleString,Float64}(
    # Locus 1
    0b000 => 0.1,
    0b100 => 0.2,
    0b001 => 0.3,
    0b101 => 0.4
  ),
  Dict{NK.AlleleString,Float64}(
    # Locus 2
    0b000 => 0.4,
    0b001 => 0.3,
    0b010 => 0.2,
    0b011 => 0.1,
  ),
  Dict{NK.AlleleString,Float64}(
    # Locus 3
    0b000 => 0.2,
    0b010 => 0.1,
    0b100 => 0.4,
    0b110 => 0.3
  )
])

ls = NK.NKLandscape(n, k, a, links, contribs)
pop = NK.Population([NK.Genotype(alleles, ls) for alleles = 0:7])

fits = [
  (0.1 + 0.4 + 0.2) / 3.0,
  (0.3 + 0.3 + 0.2) / 3.0,
  (0.1 + 0.2 + 0.1) / 3.0,
  (0.3 + 0.1 + 0.1) / 3.0,
  (0.2 + 0.4 + 0.4) / 3.0,
  (0.4 + 0.3 + 0.4) / 3.0,
  (0.2 + 0.2 + 0.3) / 3.0,
  (0.4 + 0.1 + 0.3) / 3.0
]

context("Fast functional tests") do
  context("NK.fitness(...)") do
    function testfitness(g, f0)
      f1 = NK.fitness(g)
      @fact f0 --> roughly(f1) "Expected $(f0), found $(f1), for $(g)"
    end

    for i = 1:8
      testfitness(pop.genotypes[i], fits[i])
    end
  end

  context("NK.popfits(...)") do
    fitsa = NK.popfits(pop)

    for i = 1:8
      @fact fits[i] --> roughly(fitsa[i]) "Expected $(fits[i]), found $(fitsa[i]) for $(i)th"
    end
  end

  context("NK.propsel(...)") do
    srand(0)
    newpop = NK.propsel(pop)

    counts = zeros(Int64, 8)
    for i = 1:8
      for j = 1:8
        if newpop.genotypes[i] == pop.genotypes[j]
          counts[j] += 1
        end
      end
    end
    @fact any(c -> c > 1, counts) --> true "Expected a different population after selection"
    oldmean = mean(fits)
    newmean = mean(NK.popfits(newpop))
    @fact oldmean --> less_than(newmean) "Expected greater mean fitness after selection"
  end

  context("NK.tournsel(...)") do
    srand(0)
    newpop = NK.tournsel(pop, 8)

    counts = zeros(Int64, 8)
    for i = 1:8
      for j = 1:8
        if newpop.genotypes[i] == pop.genotypes[j]
          counts[j] += 1
        end
      end
    end
    @fact any(c -> c > 1, counts) --> true "Expected a different population after selection"
    oldmean = mean(fits)
    newmean = mean(NK.popfits(newpop))
    @fact oldmean --> less_than(newmean) "Expected greater mean fitness after selection"
  end

  context("NK.moransel(...)") do
    srand(0)
    newpop = NK.moransel(pop, 16)

    counts = zeros(Int64, 8)
    for i = 1:8
      for j = 1:8
        if newpop.genotypes[i] == pop.genotypes[j]
          counts[j] += 1
        end
      end
    end
    @fact any(c -> c > 1, counts) --> true "Expected a different population after selection"
    oldmean = mean(fits)
    newmean = mean(NK.popfits(newpop))
    @fact oldmean --> less_than(newmean) "Expected greater mean fitness after selection"
  end

  context("NK.elitesel(...)") do
    srand(0)
    newpop = NK.elitesel(pop, 2)

    counts = zeros(Int64, 8)
    for i = 1:8
      for j = 1:8
        if newpop.genotypes[i] == pop.genotypes[j]
          counts[j] += 1
        end
      end
    end
    @fact any(c -> c > 1, counts) --> true "Expected a different population after selection"
    @fact newpop.genotypes |> unique |> length --> 2 "Expected two genotypes in new population"
    oldmean = mean(fits)
    newmean = mean(NK.popfits(newpop))
    @fact oldmean --> less_than(newmean) "Expected greater mean fitness after selection"
    elites = sort(pop.genotypes, by=(g) -> NK.fitness(g))[end-1:end]
    for g = newpop.genotypes
      @fact g --> anyof(elites...) "Expected members of new population to be elites"
    end
    for g = elites
      @fact g --> anyof(newpop.genotypes...) "Expected all elites to be present in new population"
    end
  end

  context("NK.bwmutate(...)") do
    srand(0)
    trials = 1000
    mutprob = 0.1
    total = 0.0
    for _ = 1:trials
      newpop = NK.bwmutate(pop, mutprob)
      for i = 1:NK.popsize(pop)
        total += (pop.genotypes[i].alleles $ newpop.genotypes[i].alleles) |> count_ones
      end
    end
    meanmuts = total / (trials * NK.popsize(pop) * n)
    @fact mutprob --> roughly(meanmuts; atol=mutprob / 10) "Expected $(mutprob) mutrate, found $(meanmuts)"
  end

  context("NK.bsmutate(...)") do
    srand(0)
    newpop = NK.bsmutate(pop, 1.0)

    for i = 1:NK.popsize(pop)
      delta = abs(count_ones(newpop.genotypes[i].alleles) - count_ones(pop.genotypes[i].alleles))
      @fact delta --> 1 "Expected 1 difference in $(pop.genotypes[i]), $(newpop.genotypes[i]), found $(delta)"
    end
  end
end


text:
using Documenter, NeuralPDE

include("pages.jl")

makedocs(
    sitename="NeuralPDE.jl",
    authors="#",
    clean=true,
    doctest=false,
    modules=[NeuralPDE],

    format=Documenter.HTML(  analytics = "UA-90474609-3",
                             assets=["assets/favicon.ico"],
                             canonical="https://neuralpde.sciml.ai/stable/"),
    pages=pages
)

deploydocs(
   repo="github.com/SciML/NeuralPDE.jl.git";
   push_preview=true
)


text:
export
    LocalSimulation,
    simulate
"""
    LocalSimulation

Describe a Local-BPS sampling
"""
struct LocalSimulation
    fg::FactorGraph
    x0::Vector{AllowedVarType}  # Starting point
    v0::Vector{AllowedVarType}  # Starting velocity
    T::Real                     # Simulation time
    maxnevents::Int             # Maximum number of events
    lambdaref::Real             # Refreshment rate
    # -- implicit
    nvars::Int                  # Number of nodes
    nfactors::Int               # Number of factors
    function LocalSimulation(fg, x0, v0, T, maxnevents, lambdaref)
        # check none of the default named arguments went through
        cond =  !(  fg == :undefined ||
                    x0 == :undefined ||
                    v0 == :undefined ||
                    T  == :undefined )
        @assert cond "One or more essential argument undefined"
        cond = (length(x0) == length(v0) == fg.structure.nvars > 0)
        @assert cond "inconsistent dims"
        cond = (T > 0.0 && maxnevents > 0 && lambdaref > 0.0)
        @assert cond "inconsistent params"
        # creating the object
        new( fg, x0, v0, T,
             maxnevents,
             lambdaref,
             length(x0),
             length(fg.factors) )
    end
end
# Constructor with named arguments
function LocalSimulation(;
            factorgraph = :undefined,
            x0          = :undefined,
            v0          = :undefined,
            T           = :undefined,
            maxnevents  = 50000,
            lambdaref   = 0.5)::LocalSimulation
    # calling the unnamed constructor
    LocalSimulation( factorgraph,
                     x0, v0, T,
                     maxnevents,
                     lambdaref )
end

"""
    simulate(sim)

Run a local BPS simulation following the specifications of `sim`.
"""
function simulate(sim::LocalSimulation)::Tuple{AllEventList, Dict}

    (start, all_evlist, pq, tref) = ls_init(sim)

    # counters for the first branch (bounce) and the second branch (refresh)
    ev_firstbranch  = 0
    ev_secondbranch = 0
    nevents         = 0 # sum of above two
    # global clock
    globalclock     = 0.0

    # global iteration
    prog = Progress(sim.maxnevents, 1)

    for evnum in 1:sim.maxnevents
        # get bounce and dequeue
        (fidx, tbounce) = peek(pq)
        dequeue!(pq)
        t = tbounce
        globalclock += t
        # ------------
        # FIRST BRANCH (t<tref => )
        if t < tref
            ev_firstbranch += 1
            ls_firstbranch!(sim.fg, fidx, all_evlist, pq, t)
        # -------------
        # SECOND BRANCH (t>=tref => refreshment)
        else
            ev_secondbranch += 1
            pq    = ls_refreshment(sim.fg, tref, all_evlist)
            tref += randexp()/sim.lambdaref
        end
        # if use while loop
        # nevents += 1
        next!(prog)
    end

    details = Dict(
        "clocktime"       => time()-start,
        "nevents"         => nevents,
        "ev_firstbranch"  => ev_firstbranch,
        "ev_secondbranch" => ev_secondbranch,
        "globalclock"     => globalclock
    )

    (all_evlist, details)
end

# =============================================================================
# NOTE: The functions with a `ls_` prefix are helper functions exclusively
# designed for the local simulation. They aren't expected to be used elsewhere.

"""
    ls_init(sim)

Initialise a local simulation: store the starting clock time, create a list of
all eventlists and push the original event, create the priorityqueue and fill
it with the initial bounce times
"""
function ls_init(sim::LocalSimulation)
    # instantiate wall clock
    start = time()
    # initialisation of the eventlists for every node in the graph
    all_evlist = AllEventList([typeof(sim.x0[i]) for i in 1:sim.nvars])
    for i in 1:sim.nvars
        evi = Event(sim.x0[i], sim.v0[i], 0.0)
        pushevent!(all_evlist.evl[i], evi)
    end
    # initialisation of the priority queue and the refreshment time
    pq   = PriorityQueue{Int, Real}()
    tref = randexp()/sim.lambdaref
    # filling of the priority queue with initial position
    for fidx in 1:sim.fg.structure.nfactors
        (xf, vf, g) = ls_retrieve(sim.fg, fidx, all_evlist, 0.0)
        ls_updatepq!(pq, sim.fg, fidx, xf, vf, g, 0.0)
    end
    return (start, all_evlist, pq, tref)
end

"""
    ls_firstbranch(fg, fidx, all_evlist, pq, t)

Operation corresponding to the first branch of events in simulate.
"""
function ls_firstbranch!(fg::FactorGraph, fidx::Int, all_evlist::AllEventList,
                         pq::PriorityQueue, t::Real)
    # retrieve xf, vf corresponding to factor
    (xf, vf, g, vars) = ls_retrieve(fg, fidx, all_evlist, t, true)
    ls_saveupdate!(all_evlist, vars, xf, vf, t)
    ls_updatepq!(pq, fg, fidx, xf, vf, g, t)
    # same story for linked factors (fp)
    for fpidx in linkedfactors(fg, fidx)
        # we don't need to retrieve `vars` here
        (xfp, vfp, gp) = ls_retrieve(fg, fpidx, all_evlist, t)
        ls_updatepq!(pq, fg, fpidx, xfp, vfp, gp, t)
    end
    return (all_evlist, pq)
end

"""
    ls_reshape(u, v)

(Local Simulation, helper function) Reshape a single vector `u` to a vector of
AllowedVarType following the structure of `v`. For example, let
`v = [0., [0., 0.]]` and `u=[1., 2., 3.]` then `ls_reshape(u,v)` will return
`[1., [2., 3.]]`.
"""
function ls_reshape(u::Vector{<:Real}, v::V) where V <: Vector{AllowedVarType}
    w = similar(v)
    # offsets to know where to look for next block of information
    offset = 0
    for i ∈ 1:length(v)
        tmpi = length(v[i])
        # if length is 1, return only single point, otherwise vector
        @inbounds w[i] = (tmpi==1) ? u[offset+1] : u[offset .+ (1:tmpi)]
        offset += tmpi
    end
    return w
end

"""
    ls_retrieve(fg, fidx, all_evlist, t, doreflect)

(Local simulation, helper function) Retrieve the structure corresponding to all
positions of the nodes attached to a factor (`xf`) and the velocities (`vf`) as
well as the gradient at `xf` and the list of indices of the attached variables.
"""
function ls_retrieve(fg::FactorGraph, fidx::Int,
                    all_evlist::AllEventList,
                    t::Real, doreflect::Bool=false)
    # indices of the variable associated with factor fidx
    vars = assocvariables(fg, fidx)
    # allocate xf, vf, note the different variables
    # don't necessarily have the same types
    vf = Vector{AllowedVarType}(undef, length(vars))
    xf = similar(vf)
    # shortcut to initial event
    if t == 0.0
        for (i, k) ∈ enumerate(vars)
            ev = getevent(all_evlist.evl[k], 1)
            xf[i], vf[i] = ev.x, ev.v
        end
    # retrieving later than initial event (standard case)
    else
        # retrieve vf, xf(t)
        for (i, k) ∈ enumerate(vars)
            # get the eventlist corresponding to variable k
            ev = getlastevent(all_evlist.evl[k])
            # store the components and extrapolate from ev.x as needed
            vf[i] = ev.v
            xf[i] = ev.x + (t-ev.t)*ev.v
        end
    end
    # compute the gradient of the factor's log likelihood at xf
    g  = fg.factors[fidx].gll(vcat(xf...))
    # compute the reflection (if required)
    vf = doreflect ? ls_reshape(reflect_bps!(g, vcat(vf...)), vf) : vf
    (xf, vf, g, vars)
end

"""
    ls_saveupdate!(all_evlist, vars, xf, vf, t)

(Local simulation, helper function) Append the information contained in `xf`,
`vf` in the appropriate EventLists corresponding to the relevant variables
indexed in `vars` (variables associated with the factor `f`, `xf` and `vf` are
the recovered positions and velocities obtained when considering factor `f`).
"""
function ls_saveupdate!(all_evlist::AllEventList, vars::Vector{Int},
                        xf::Vector{AllowedVarType}, vf::Vector{AllowedVarType},
                        t::Real)
    # Add the new sample to the list
    for (i, k) ∈ enumerate(vars)
        pushevent!(all_evlist.evl[k], Event(xf[i], vf[i], t))
    end
    return all_evlist
end

"""
    ls_updatepq!(pq, fg, fidx, xf, vf, g, t)

(Local simulation, helper function) Compute the bounce time (possibly via
thinning depending on the factor's implementation) associated with a factor and
push it to the PriorityQueue `pq`.
"""
function ls_updatepq!(pq::PriorityQueue, fg::FactorGraph, fidx::Int,
                      xf::Vector{AllowedVarType}, vf::Vector{AllowedVarType},
                      g::Vector{<:Real}, t::Real)
    # useful temporary variables
    vcxf, vcvf = vcat(xf...), vcat(vf...)
    # Update time in Priority Queue for the current factor
    acc  = false
    tauf = 0.0
    while !acc
        bounce = fg.factors[fidx].nextevent(vcxf, vcvf)
        acc, tauf = bounce.dobounce(g, vcvf), bounce.tau
    end
    # add the time to current priorityQueue
    pq[fidx] = t+tauf
    return pq
end

"""
    ls_random(evl)

(Local simulation, helper function) Generate N(0,1) random numbers of dimension
corresponding to the passed `EventList`.
"""
function ls_random(evl::EventList)
    l = length(evl.xs[end])
    return l > 1 ? randn(l) : randn()
end

"""
    ls_refreshment(fg, t, all_evlist)

(Local simulation, helper function) Generate a new PriorityQueue by drawing
new velocities for each node and storing it (cf. NewQueue algorithm in the BPS
paper).
"""
function ls_refreshment(fg::FactorGraph, t::Real, all_evlist::AllEventList)
    # draw a new bunch of velocities for each node
    v = Vector{AllowedVarType}(undef, fg.structure.nvars)
    for i ∈ 1:fg.structure.nvars
        @inbounds v[i] = ls_random(all_evlist.evl[i])
    end
    # Instantiate a new priority queue
    pq = PriorityQueue{Int, Real}()
    for fidx ∈ 1:fg.structure.nfactors
        # retrieve xf, vf corresponding to factor
        (xf, vf, g, vars) = ls_retrieve(fg, fidx, all_evlist, t)
        # update the priority queue, note v[vars] and not vf
        # since we are using the refreshed velocity
        ls_updatepq!(pq, fg, fidx, xf, v[vars], g, t)
    end
    return pq
end


text:
module NoNewline

"f"
f(x) = x

end # EOF at end of this line.


text:
nodes_test= [Bus(1 , #number
                 "Bus 1", #Name
                 "REF" , #BusType (REF, PV, PQ)
                 0, #Angle in radians
                 1.02, #Voltage in pu
                 (min=0.94, max=1.06), #Voltage limits in pu
                 69), #Base voltage in kV
             Bus(2 , "Bus 2"  , "PV" , 0 , 1.0 , (min=0.94, max=1.06), 69)]

nodes_test2= [Bus(1 , #number
              "Bus 1", #Name
              "REF" , #BusType (REF, PV, PQ)
              0, #Angle in radians
              1.02, #Voltage in pu
              (min=0.94, max=1.06), #Voltage limits in pu
              69), #Base voltage in kV
          Bus(2 , "Bus 2"  , "PV" , 0 , 1.0 , (min=0.94, max=1.06), 69),
          Bus(3 , "Bus 3"  , "PQ" , 0 , 1.0 , (min=0.94, max=1.06), 69)]

branch_test = [Line("Line1", #name
                    true, #available
                    0.0, #active power flow initial condition (from-to)
                    0.0, #reactive power flow initial condition (from-to)
                    Arc(from=nodes_test[1], to=nodes_test[2]), #Connection between buses
                    0.01, #resistance in pu
                    0.05, #reactance in pu
                    (from=0.0, to=0.0), #susceptance in pu
                    18.046, #rate in MW
                    1.04)]  #angle limits (-min and max)

branch_test_fault = [Line("Line1", #name
                    true, #available
                    0.0, #active power flow initial condition (from-to)
                    0.0, #reactive power flow initial condition (from-to)
                    Arc(from=nodes_test[1], to=nodes_test[2]), #Connection between buses
                    0.02, #resistance in pu
                    0.1, #reactance in pu
                    (from=0.0, to=0.0), #susceptance in pu
                    18.046, #rate in MW
                    1.04)]  #angle limits (-min and max)

####

branch_test2 = [Line("Line1", #name
                    true, #available
                    0.0, #active power flow initial condition (from-to)
                    0.0, #reactive power flow initial condition (from-to)
                    Arc(from=nodes_test2[1], to=nodes_test2[2]), #Connection between buses
                    0.01, #resistance in pu
                    0.05, #reactance in pu
                    (from=0.0, to=0.0), #susceptance in pu
                    18.046, #rate in MW
                    1.04),
                    Line("Line2", #name
                            true, #available
                            0.0, #active power flow initial condition (from-to)
                            0.0, #reactive power flow initial condition (from-to)
                            Arc(from=nodes_test2[2], to=nodes_test2[3]), #Connection between buses
                            0.02, #resistance in pu
                            0.1, #reactance in pu
                            (from=0.0, to=0.0), #susceptance in pu
                            18.046, #rate in MW
                            1.04)]  #angle limits (-min and max)

branch_test2_fault = [Line("Line1", #name
                true, #available
                0.0, #active power flow initial condition (from-to)
                0.0, #reactive power flow initial condition (from-to)
                Arc(from=nodes_test2[1], to=nodes_test2[2]), #Connection between buses
                0.02, #resistance in pu
                0.1, #reactance in pu
                (from=0.0, to=0.0), #susceptance in pu
                18.046, #rate in MW
                1.04),
                Line("Line2", #name
                        true, #available
                        0.0, #active power flow initial condition (from-to)
                        0.0, #reactive power flow initial condition (from-to)
                        Arc(from=nodes_test2[2], to=nodes_test2[3]), #Connection between buses
                        0.02, #resistance in pu
                        0.1, #reactance in pu
                        (from=0.0, to=0.0), #susceptance in pu
                        18.046, #rate in MW
                        1.04)]  #angle limits (-min and max)


text:
module SoapySDR

using soapysdr_jll
const lib = soapysdr_jll.libsoapysdr

using Intervals
using Unitful
using Unitful.DefaultSymbols
const dB = u"dB"
const GC = Base.GC

export @u_str

include("lowlevel/Constants.jl") # Done
include("lowlevel/Errors.jl")    # Done
include("lowlevel/Formats.jl")   # Done
include("lowlevel/Types.jl")     # Done
include("lowlevel/Device.jl")
include("highlevel.jl")

const SDRStream = Stream
export SDRStream

end


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license
#
# =========================================================================== #
# !!!! AUTO-GENERATED FILE !!!! Do not edit. See the bindgen repository for
# the generation code (https://github.com/chemfiles/bindgen).
# This file contains Julia interface to the C API
# =========================================================================== #

# === Manually translated from the header
const Cbool = Cuchar
const chfl_vector3d = Array{Cdouble, 1}

struct chfl_match
    size    ::UInt64
    atoms_1 ::UInt64
    atoms_2 ::UInt64
    atoms_3 ::UInt64
    atoms_4 ::UInt64
end

struct chfl_format_metadata
    name :: Ptr{Cchar}
    extension :: Ptr{Cchar}
    description :: Ptr{Cchar}
    reference :: Ptr{Cchar}
    read :: Cbool
    write :: Cbool
    memory :: Cbool
    positions :: Cbool
    velocities :: Cbool
    unit_cell :: Cbool
    atoms :: Cbool
    bonds :: Cbool
    residues :: Cbool
end
# === End of manual type defintion

struct CHFL_TRAJECTORY end

struct CHFL_CELL end

struct CHFL_ATOM end

struct CHFL_FRAME end

struct CHFL_TOPOLOGY end

struct CHFL_SELECTION end

struct CHFL_RESIDUE end

struct CHFL_PROPERTY end

# enum chfl_status
const chfl_status = UInt32
const CHFL_SUCCESS = chfl_status(0)
const CHFL_MEMORY_ERROR = chfl_status(1)
const CHFL_FILE_ERROR = chfl_status(2)
const CHFL_FORMAT_ERROR = chfl_status(3)
const CHFL_SELECTION_ERROR = chfl_status(4)
const CHFL_CONFIGURATION_ERROR = chfl_status(5)
const CHFL_OUT_OF_BOUNDS = chfl_status(6)
const CHFL_PROPERTY_ERROR = chfl_status(7)
const CHFL_GENERIC_ERROR = chfl_status(254)
const CHFL_CXX_ERROR = chfl_status(255)

# enum chfl_bond_order
const chfl_bond_order = UInt32
const CHFL_BOND_UNKNOWN = chfl_bond_order(0)
const CHFL_BOND_SINGLE = chfl_bond_order(1)
const CHFL_BOND_DOUBLE = chfl_bond_order(2)
const CHFL_BOND_TRIPLE = chfl_bond_order(3)
const CHFL_BOND_QUADRUPLE = chfl_bond_order(4)
const CHFL_BOND_QUINTUPLET = chfl_bond_order(5)
const CHFL_BOND_AMIDE = chfl_bond_order(254)
const CHFL_BOND_AROMATIC = chfl_bond_order(255)

# enum chfl_property_kind
const chfl_property_kind = UInt32
const CHFL_PROPERTY_BOOL = chfl_property_kind(0)
const CHFL_PROPERTY_DOUBLE = chfl_property_kind(1)
const CHFL_PROPERTY_STRING = chfl_property_kind(2)
const CHFL_PROPERTY_VECTOR3D = chfl_property_kind(3)

# enum chfl_cellshape
const chfl_cellshape = UInt32
const CHFL_CELL_ORTHORHOMBIC = chfl_cellshape(0)
const CHFL_CELL_TRICLINIC = chfl_cellshape(1)
const CHFL_CELL_INFINITE = chfl_cellshape(2)


text:
export TxDAQControllerParams, TxDAQController, controlTx

Base.@kwdef mutable struct TxDAQControllerParams <: DeviceParams
  phaseAccuracy::Float64
  amplitudeAccuracy::Float64
  controlPause::Float64
  maxControlSteps::Int64 = 20
  correctCrossCoupling::Bool = false
end
TxDAQControllerParams(dict::Dict) = params_from_dict(TxDAQControllerParams, dict)

struct ControlledChannel
  seqChannel::PeriodicElectricalChannel
  daqChannel::TxChannelParams
end

Base.@kwdef mutable struct TxDAQController <: VirtualDevice
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::TxDAQControllerParams
  "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}}

  currTx::Union{Matrix{ComplexF64}, Nothing} = nothing
  desiredTx::Union{Matrix{ComplexF64}, Nothing} = nothing
  sinLUT::Union{Matrix{Float64}, Nothing} = nothing
  cosLUT::Union{Matrix{Float64}, Nothing} = nothing
  controlledChannels::Vector{ControlledChannel} = []
end

function init(tx::TxDAQController)
  @info "Initializing TxDAQController with ID `$(tx.deviceID)`."
  tx.present = true
end

neededDependencies(::TxDAQController) = [AbstractDAQ]
optionalDependencies(::TxDAQController) = []

function controlTx(txCont::TxDAQController, seq::Sequence, initTx::Union{Matrix{ComplexF64}, Nothing} = nothing)
  # Prepare and check channel under control
  daq = dependency(txCont, AbstractDAQ)
  seqControlledChannel = getControlledChannel(seq)
  missingControlDef = []
  txCont.controlledChannels = []
  for seqChannel in seqControlledChannel
    name = id(seqChannel)
    daqChannel = get(daq.params.channels, name, nothing)
    if isnothing(daqChannel) || isnothing(daqChannel.feedback) || !in(daqChannel.feedback.channelID, daq.refChanIDs)
      push!(missingControlDef, name)
    else
      push!(txCont.controlledChannels, ControlledChannel(seqChannel, daqChannel))
    end
  end
  if length(missingControlDef) > 0
    message = "The sequence requires control for the following channel " * join(string.(missingControlDef), ", ", " and") * ", but either the channel was not defined or had no defined feedback channel."
    throw(IllegalStateException(message))
  end

  # Check that channels only have one component
  if any(x -> length(x.components) > 1, seqControlledChannel)
    throw(IllegalStateException("Sequence has channel with more than one component. Such a channel cannot be controlled by this controller"))
  end

  if !isnothing(initTx)
    s = size(initTx)
    # Not square or does not fit controlled channel matrix
    if !(length(s) == 0 || all(isequal(s[1]), s))
      throw(IllegalStateException("Given initTx for control tx has dimenions $s that is either not square or does not match the amount of controlled channel"))
    end
  end

  # Prepare values
  if isnothing(initTx)
    txCont.currTx = convert(Matrix{ComplexF64}, diagm(ustrip.(u"V", [channel.daqChannel.limitPeak/10 for channel in txCont.controlledChannels])))
  else
    txCont.currTx = initTx
  end
  sinLUT, cosLUt = createLUTs(seqControlledChannel, seq::Sequence)
  txCont.sinLUT = sinLUT
  txCont.cosLUT = cosLUt
  Ω = calcDesiredField(seqControlledChannel)

  # Start Tx
  setTxParams(daq, txFromMatrix(txCont, txCont.currTx)...)
  startTx(daq)

  controlPhaseDone = false
  i = 1
  while !controlPhaseDone && i <= txCont.params.maxControlSteps
    @info "CONTROL STEP $i"
    period = currentPeriod(daq)
    uMeas, uRef = readDataPeriods(daq, 1, period + 1, acqNumAverages(seq))

    controlPhaseDone = doControlStep(txCont, seq, uRef, Ω)

    sleep(txCont.params.controlPause)
    i += 1
  end
  stopTx(daq)
  setTxParams(daq, txFromMatrix(txCont, txCont.currTx)...)
  return txFromMatrix(txCont, txCont.currTx)
end

getControlledChannel(seq::Sequence) = [channel for field in seq.fields if field.control for channel in field.channels if typeof(channel) <: PeriodicElectricalChannel]

function txFromMatrix(txCont::TxDAQController, Γ::Matrix{ComplexF64})
  amplitudes = Dict{String, Vector{Union{Float32, Nothing}}}()
  phases = Dict{String, Vector{Union{Float32, Nothing}}}()
  for (d, channel) in enumerate(txCont.controlledChannels)
    amps = []
    phs = []
    for (e, channel) in enumerate(txCont.controlledChannels)
      push!(amps, abs(Γ[d, e]))
      push!(phs, angle(Γ[d, e]))
    end
    amplitudes[id(channel.seqChannel)] = amps
    phases[id(channel.seqChannel)] = phs
  end
  return amplitudes, phases
end

function createLUTs(seqChannel::Vector{PeriodicElectricalChannel}, seq::Sequence)
  N = rxNumSamplingPoints(seq)
  D = length(seqChannel)
  cycle = ustrip(dfCycle(seq))
  base = ustrip(dfBaseFrequency(seq))
  dfFreq = [base/x.components[1].divider for x in seqChannel]

  sinLUT = zeros(N,D)
  cosLUT = zeros(N,D)
  for d=1:D
    Y = round(Int64, cycle*dfFreq[d] )
    for n=1:N
      sinLUT[n,d] = sin(2 * pi * (n-1) * Y / N) / N #sqrt(N)*2
      cosLUT[n,d] = cos(2 * pi * (n-1) * Y / N) / N #sqrt(N)*2
    end
  end
  return sinLUT, cosLUT
end

function doControlStep(txCont::TxDAQController, seq::Sequence, uRef, Ω::Matrix)

  Γ = calcFieldFromRef(txCont,seq, uRef)
  daq = dependency(txCont, AbstractDAQ)

  @info "reference Γ=" Γ

  if controlStepSuccessful(Γ, Ω, txCont)
    @info "Could control"
    return true
  else
    newTx = newDFValues(Γ, Ω, txCont)
    oldTx = txCont.currTx
    @info "oldTx=" oldTx
    @info "newTx=" newTx

    if checkDFValues(newTx, oldTx, Γ,txCont)
      txCont.currTx[:] = newTx
      setTxParams(daq, txFromMatrix(txCont, txCont.currTx)...)
    else
      @warn "New values are above voltage limits or are different than expected!"
    end

    @info "Could not control !"
    return false
  end
end

function calcFieldFromRef(txCont::TxDAQController, seq::Sequence, uRef)
  len = length(txCont.controlledChannels)
  Γ = zeros(ComplexF64, len, len)

  for d=1:len
    for e=1:len
      c = ustrip(u"T/V", txCont.controlledChannels[d].daqChannel.feedback.calibration)

      uVolt = float(uRef[1:rxNumSamplingPoints(seq),d,1])

      a = 2*sum(uVolt.*txCont.cosLUT[:,e])
      b = 2*sum(uVolt.*txCont.sinLUT[:,e])

      Γ[d,e] = c*(b+im*a)
    end
  end
  return Γ
end

function calcDesiredField(seqChannel::Vector{PeriodicElectricalChannel})
  temp = [ustrip(amplitude(ch.components[1])) * exp(im*ustrip(phase(ch.components[1]))) for ch in seqChannel]
  return convert(Matrix{ComplexF64}, diagm(temp))
end

function controlStepSuccessful(Γ::Matrix, Ω::Matrix, txCont::TxDAQController)

  if txCont.params.correctCrossCoupling
    diff = Ω - Γ
  else
    diff = diagm(diag(Ω)) - diagm(diag(Γ))
  end
  deviation = maximum(abs.(diff)) / maximum(abs.(Ω))
  @info "Ω = " Ω
  @info "Γ = " Γ
  @info "Ω - Γ = " diff
  @info "deviation = $(deviation)   allowed= $(txCont.params.amplitudeAccuracy)"
  return deviation < txCont.params.amplitudeAccuracy
end

function newDFValues(Γ::Matrix, Ω::Matrix, txCont::TxDAQController)

  κ = txCont.currTx
  if txCont.params.correctCrossCoupling
    β = Γ*inv(κ)
  else
    @show size(Γ), size(κ)
    β = diagm(diag(Γ))*inv(diagm(diag(κ)))
  end
  newTx = inv(β)*Ω

  @warn "here are the values"
  @show κ
  @show Γ
  @show Ω

  return newTx
end

function checkDFValues(newTx, oldTx, Γ, txCont::TxDAQController)

  calibFieldToVoltEstimate = [ustrip(u"V/T", ch.daqChannel.calibration) for ch in txCont.controlledChannels]
  calibFieldToVoltMeasured = abs.(diag(oldTx) ./ diag(Γ))

  @info "" calibFieldToVoltEstimate[1] calibFieldToVoltMeasured[1]

  deviation = abs.(1.0 .- calibFieldToVoltMeasured./calibFieldToVoltEstimate)

  @info "We expected $(calibFieldToVoltEstimate) and got $(calibFieldToVoltMeasured), deviation: $deviation"

  return all( abs.(newTx) .<  ustrip.(u"V", [channel.daqChannel.limitPeak for channel in txCont.controlledChannels]) ) && maximum( deviation ) < 0.2
end