text:
using LinearAlgebra, SparseArrays, Test
include("inconsistentvrep.jl")

using JuMP, StaticArrays

function change_fulldim_test(Rep)
    T = Int
    RepT = Rep{T, Matrix{T}}
    @testset "Change Polyhedra.FullDim with $RepT" begin
        @test fulldim(RepT) == -1
        @test Polyhedra.coefficient_type(RepT) == T
        changedrep = Polyhedra.similar_type(RepT, 10)
        @test fulldim(changedrep) == -1
        @test Polyhedra.FullDim(changedrep) == -1
        #@test (@inferred Polyhedra.FullDim(changedrep)) == 10
        @test Polyhedra.coefficient_type(changedrep) == T
    end
end

function isempty_vtest(vr::VRepresentation, nr, np)
    @testset "isempty not working correctly for iterators #17 with V-rep" begin
        hasr = nr > 0
        hasp = np > 0
        @test npoints(vr) == length(points(vr)) == np
        @test nallrays(vr) == length(allrays(vr)) == nr
        @test nlines(vr) == length(lines(vr)) == 0
        @test nrays(vr) == length(rays(vr)) == nr
        @test haspoints(vr) == !isempty(points(vr)) == hasp
        @test hasallrays(vr) == !isempty(allrays(vr)) == hasr
        @test haslines(vr) == !isempty(lines(vr)) == false
        @test hasrays(vr) == !isempty(rays(vr)) == hasr
    end
end

function isempty_htest(hr::HRepresentation, ne, ni)
    @testset "isempty not working correctly for iterators #17 with H-rep" begin
        hase = ne > 0
        hasi = ni > 0
        @test nallhalfspaces(hr) == length(allhalfspaces(hr)) == 2ne + ni
        @test nhyperplanes(hr) == length(hyperplanes(hr)) == ne
        @test nhalfspaces(hr) == length(halfspaces(hr)) == ni
        @test hasallhalfspaces(hr) == !isempty(allhalfspaces(hr)) == hase || hasi
        @test hashyperplanes(hr) == !isempty(hyperplanes(hr)) == hase
        @test hashalfspaces(hr) == !isempty(halfspaces(hr)) == hasi
    end
end

function htest(hr::Polyhedra.HRepresentation, AT::Type{<:AbstractVector})
    @test (@inferred Polyhedra.coefficient_type(hr)) == Float64
    @test                                               (@inferred eltype(allhalfspaces(hr)))  == HalfSpace{Float64, AT}
    @test                                               (@inferred collect(allhalfspaces(hr))) isa Vector{HalfSpace{Float64, AT}}
    @test isempty(allhalfspaces(hr)) == iszero(nallhalfspaces(hr))
    @test (@inferred Polyhedra.halfspacetype(hr))    == (@inferred eltype(halfspaces(hr)))     == HalfSpace{Float64, AT}
    @test                                               (@inferred collect(halfspaces(hr)))    isa Vector{HalfSpace{Float64, AT}}
    @test isempty(halfspaces(hr)) == iszero(nhalfspaces(hr))
    @test (@inferred Polyhedra.hyperplanetype(hr))   == (@inferred eltype(hyperplanes(hr)))    == HyperPlane{Float64, AT}
    @test                                               (@inferred collect(hyperplanes(hr)))   isa Vector{HyperPlane{Float64, AT}}
    @test isempty(hyperplanes(hr)) == iszero(nhyperplanes(hr))
    @test_throws DimensionMismatch ones(2, 3) \ hr
end

function vtest(vr::VRepresentation, AT::Type{<:AbstractVector})
    @test (@inferred Polyhedra.coefficient_type(vr)) == Int
    @test (@inferred Polyhedra.pointtype(vr))    == (@inferred eltype(points(vr)))     == AT
    @test                                           (@inferred collect(points(vr)))    isa Vector{AT}
    @test isempty(points(vr)) == iszero(npoints(vr))
    @test                                           (@inferred eltype(allrays(vr)))    == Ray{Int, AT}
    @test                                           (@inferred collect(allrays(vr)))   isa Vector{Ray{Int, AT}}
    @test isempty(allrays(vr)) == iszero(nallrays(vr))
    @test (@inferred Polyhedra.linetype(vr))     == (@inferred eltype(lines(vr)))      == Line{Int, AT}
    @test                                           (@inferred collect(lines(vr)))     isa Vector{Line{Int, AT}}
    @test isempty(lines(vr)) == iszero(nlines(vr))
    @test (@inferred Polyhedra.raytype(vr))      == (@inferred eltype(rays(vr)))       == Ray{Int, AT}
    @test                                           (@inferred collect(rays(vr)))      isa Vector{Ray{Int, AT}}
    @test isempty(rays(vr)) == iszero(nrays(vr))
    @test_throws DimensionMismatch ones(2, 1) * vr
end

function eltype_incorrect()
    hps = [HyperPlane([1, 2, 3], 7.)]
    shps = [@inferred HyperPlane((@SVector [1, 2, 3]), 7.)]
    @test eltype(shps) == HyperPlane{Float64, SVector{3, Float64}}
    hss = [HalfSpace([4, 5, 6.], 8)]
    shss = [@inferred HalfSpace((@SVector [4., 5., 6.]), 8)]
    @test eltype(shss) == HalfSpace{Float64, SVector{3, Float64}}
    #htest((@inferred hrep(hps)), Vector{Float64})
    htest(hrep(hps), Vector{Float64})
    #htest((@inferred hrep(shps)), SVector{3, Float64})
    htest((hrep(shps)), SVector{3, Float64})
    #htest((@inferred hrep(hss)), Vector{Float64})
    htest(hrep(hss), Vector{Float64})
    #htest((@inferred hrep(shss)), SVector{3, Float64})
    htest(hrep(shss), SVector{3, Float64})
    #htest((@inferred hrep(hps, hss)), Vector{Float64})
    htest(hrep(hps, hss), Vector{Float64})
    #htest((@inferred hrep(shps, shss)), SVector{3, Float64})
    htest(hrep(shps, shss), SVector{3, Float64})
    htest(hrep([1 2 3; 4 5 6], [7., 8], BitSet([1])), Vector{Float64})
    htest(hrep(@SMatrix([1 2 3; 4 5 6]), @SVector([7., 8]), BitSet([1])), SVector{3, Float64})
    htest(hrep(spzeros(2, 3), [7., 8], BitSet([1])), SparseVector{Float64, Int})
    htest(hrep([1 2 3; 4 5 6], [7., 8]), Vector{Float64})
    htest(hrep(@SMatrix([1 2 3; 4 5 6]), @SVector([7., 8])), SVector{3, Float64})
    ps = [[1, 2], [3, 4]]
    sps = [(@SVector [1, 2]), (@SVector [3, 4])]
    rs = [Ray([0, 1])]
    srs = [Ray(@SVector [0, 1])]
    ls = [Line([0, 1])]
    sls = [Line(@SVector [0, 1])]
    @test eltype(sps) == SVector{2, Int}
    vtest(vrep(ps), Vector{Int})
    #vtest((@inferred vrep(sps)), SVector{2, Int})
    vtest(vrep(sps), SVector{2, Int})
    vtest(convexhull(ps...), Vector{Int})
    #vtest((@inferred convexhull(sps...)), SVector{2, Int})
    vtest(convexhull(sps...), SVector{2, Int})
    #vtest((@inferred vrep(ls)), Vector{Int})
    vtest(vrep(ls), Vector{Int})
#        vtest((@inferred vrep(sls)), SVector{2, Int})
#        vtest((@inferred vrep(rs)), Vector{Int})
#        vtest((@inferred vrep(srs)), SVector{2, Int})
#        vtest((@inferred vrep(ls, rs)), Vector{Int})
#        vtest((@inferred vrep(sls, srs)), SVector{2, Int})
#        vtest((@inferred vrep(ps, ls, rs)), Vector{Int})
#        vtest((@inferred vrep(sps, sls, srs)), SVector{2, Int})
    vtest(vrep(sls), SVector{2, Int})
    vtest(vrep(rs), Vector{Int})
    vtest(vrep(srs), SVector{2, Int})
    vtest(vrep(ls, rs), Vector{Int})
    vtest(vrep(sls, srs), SVector{2, Int})
    vtest(vrep(ps, ls, rs), Vector{Int})
    vtest(vrep(sps, sls, srs), SVector{2, Int})
    vtest(vrep([1 2; 3 4]), Vector{Int})
    vtest(vrep(spzeros(Int, 2, 2)), SparseVector{Int, Int})
    vtest(vrep([1 2; 3 4], zeros(Int, 0, 0), BitSet()), Vector{Int})
    vtest(vrep([1 2; 3 4], zeros(Int, 0, 0)), Vector{Int})
    vtest(vrep([1 2; 3 4]), Vector{Int})
end

function polar_test()
    v = convexhull([-1, 2]) + conichull([1, 0], Line([1, 0]))
    for p in [v, polyhedron(v)]
        h = polar(p)
        @test collect(hyperplanes(h)) == [HyperPlane([1, 0], 0)]
        @test collect(halfspaces(h)) == [HalfSpace([-1, 2], 1), HalfSpace([1, 0], 0)]
    end
    h = intersect(HalfSpace([1, 0], 1e-13), HyperPlane([1, 1], -1e-13))
    v = polar(h)
    @test collect(points(v)) == [[0.0, 0.0]]
    @test collect(lines(v)) == [Line([1.0, 1.0])]
    @test collect(rays(v)) == [Ray([1.0, 0.0])]
    h = intersect(HalfSpace([1, 0], 2), HalfSpace([1, 1], 0), HyperPlane([0, 1], 0))
    v = polar(h)
    @test collect(points(v)) == [[0.5, 0.0]]
    @test collect(lines(v)) == [Line([0.0, 1.0])]
    @test collect(rays(v)) == [Ray([1.0, 1.0])]
    h = intersect(HalfSpace([1, 0], 2), HalfSpace([1, 1], 0))
    for p in [h, polyhedron(h)]
        v = polar(p)
        @test collect(points(v)) == [[0.5, 0.0]]
        @test !haslines(v)
        @test collect(rays(v)) == [Ray([1.0, 1.0])]
    end
end

@testset "Representation tests" begin

    include("vecrep.jl")
    include("lphrep.jl")
    include("matrep.jl")
    include("liftedrep.jl")

    @testset "eltype for some iterators is incorrect #7" begin
        eltype_incorrect()
    end

    @testset "Iterating over halfspaces of a MixedMatHRep broken #9" begin
        A = [1 2; 3 4; 5 6]
        b = [1, 2, 3]
        halfspace = [1, 3]
        hyperplane = [2]
        linset = BitSet(2)
        hr = hrep(A, b, linset)
        Aall = [3 4; -3 -4; 1 2; 5 6]
        ball = [2, -2, 1, 3]
        for (i, h) in enumerate(allhalfspaces(hr))
            @test h.a == Aall[i, :]
            @test h.β == ball[i]
            @test isa(h, HalfSpace{Int})
        end
        for (i, h) in enumerate(halfspaces(hr))
            @test h.a == A[halfspace[i], :]
            @test h.β == b[halfspace[i]]
            @test isa(h, HalfSpace{Int})
        end
        for (i, h) in enumerate(hyperplanes(hr))
            @test h.a == A[hyperplane[i], :]
            @test h.β == b[hyperplane[i]]
            @test isa(h, HyperPlane{Int})
        end
    end

    @testset "Building rep with different type" begin
        @test Polyhedra.coefficient_type(MixedMatHRep{Float64}([1 2; 3 4], [1, 2], BitSet())) == Float64
        @test Polyhedra.coefficient_type(MixedMatVRep{Float64}([1 2; 3 4], [1 2; 3 4], BitSet())) == Float64
        @test Polyhedra.coefficient_type(LiftedHRepresentation{Float64}([1 2; 3 4])) == Float64
        @test Polyhedra.coefficient_type(LiftedVRepresentation{Float64}([1 2; 3 4])) == Float64
    end

    @testset "Radius of largest inscribed ball at a given center" begin
        p = hrep(Matrix(1I, 2, 2), zeros(2))
        r = Polyhedra.maximum_radius_with_center(p, zeros(2))
        @test r == 0

        p = hrep([0 0], [0])
        @test_throws ErrorException Polyhedra.maximum_radius_with_center(p, zeros(2))
        p = hrep([1 0], [-1])
        @test_throws ErrorException Polyhedra.maximum_radius_with_center(p, zeros(2))

        A = [ 2  1
              2 -1
             -1  2
             -1 -2]
        b = ones(4)
        p = hrep(A, b)
        r = Polyhedra.maximum_radius_with_center(p, zeros(2))
        @test r ≈ inv(√5) atol=1e-6

        p = polyhedron(vrep([ 0  0
                              1  0
                              1  1
                              0  1]))
        @test Polyhedra.maximum_radius_with_center(p, [0.00, 0.00]) == 0.00
        @test Polyhedra.maximum_radius_with_center(p, [1.00, 1.00]) == 0.00
        @test Polyhedra.maximum_radius_with_center(p, [0.50, 0.50]) == 0.50
        @test Polyhedra.maximum_radius_with_center(p, [0.50, 0.25]) == 0.25
        @test Polyhedra.maximum_radius_with_center(p, [0.25, 0.50]) == 0.25
        @test Polyhedra.maximum_radius_with_center(p, [0.25, 0.25]) == 0.25
    end

    @testset "Chebyshev center" begin
        p = hrep(Matrix(1I, 2, 2), zeros(2))
        @test_throws ErrorException chebyshevcenter(p, lp_solver) # unbounded

        p = hrep([1 1; -1 -1], [0, -1])
        @test_throws ErrorException chebyshevcenter(p, lp_solver) # empty

        # examples/chebyshevcenter.ipynb
        A = [ 2  1
              2 -1
             -1  2
             -1 -2]
        b = ones(4)
        p = hrep(A, b)
        c, r = chebyshevcenter(p, lp_solver)
        @test c ≈ [0, 0] atol=1e-6
        @test r ≈ 0.4472135955 atol=1e-6

        _interval(c, r) = HalfSpace([1], c + r) ∩ HalfSpace([-1], -c + r)
        c_exp = [-4, 3, 8]
        for rs in [[1, 2, 3], [3, 1, 2], [3, 2, 1], [2, 1, 2]]
            h = _interval(c_exp[1], rs[1]) * _interval(c_exp[2], rs[2]) * _interval(c_exp[3], rs[3])
            c, r = hchebyshevcenter(h, lp_solver, verbose=1)
            @test c ≈ c_exp atol=1e-6
            @test r ≈ 1 atol=1e-6

            c, r = hchebyshevcenter(h, lp_solver, proper=false)
            for i in 1:3
                @test c_exp[i] - rs[i] + 1 - 1e-6 <= c[i]
                @test c[i] <= c_exp[i] + rs[i] - 1 + 1e-6
            end
            @test r ≈ 1 atol=1e-6
        end

        h = _interval(0, 0) * _interval(0, 1)
        c, r = hchebyshevcenter(h, lp_solver, linearity_detected=true, verbose=1)
        @test c[1] ≈ 0 atol=1e-6
        @test -1 - 1e-6 <= c[2] <= 1 + 1e-6
        @test r ≈ 0 atol=1e-6

        p = convexhull([0, 0], [0, 1], [1, 0])
        @test_throws ErrorException chebyshevcenter(p) # Not yet implemented
    end

    @testset "V-consistency with iterator constructor" begin
        T = Int
        AT = Vector{Int}
        for VRepType in (Polyhedra.LiftedVRepresentation{T, Matrix{T}},
                         Polyhedra.MixedMatVRep{T, Matrix{T}},
                         Polyhedra.Hull{T, AT, Int})
            @test_throws ErrorException VRepType(2, AT[], [Line([1, 2])])
            @test_throws ErrorException VRepType(2, AT[], Line{T, AT}[], [Ray([1, 2])])
            @test_throws ErrorException VRepType(2, AT[], [Line([1, 2])], [Ray([1, 2])])
            v = VRepType(2, [Line([1, 2])])
            @test collect(points(v)) == [[0, 0]]
            @test collect(lines(v)) == [Line([1, 2])]
            @test !hasrays(v)
        end
        for vinc in (InconsistentVRep{T, AT, Int}(2, AT[], Line{T, AT}[], [Ray([1, 2])]),
                     InconsistentVRep{T, AT, Int}(2, AT[], [Line([1, 2])], Ray{T, AT}[]),
                     InconsistentVRep{T, AT, Int}(2, AT[], [Line([1, 2])], [Ray([1, 2])]))
            @test_throws ErrorException Polyhedra.checkvconsistency(vinc)
            pinc = polyhedron(vinc)
            @test_throws ErrorException Polyhedra.checkvconsistency(pinc)
        end
        let
            AT = StaticArrays.SVector{1, Int}
            VRepType = Polyhedra.Hull{T, AT, Int}
            @test isempty(VRepType(2, AT[], Line{T, AT}[], Ray{T, AT}[]))
            @test isempty(VRepType(2, Line{T, AT}[], Ray{T, AT}[]))
        end
    end
    @testset "Combination of different coefficient type" begin
        @testset "V-representation" begin
            generator_fulltest(convexhull([1, 0], Line([0, 1.])), convexhull(Line([0, 1]), [1, 0.]))
            @test conichull(convexhull([1, 0.]), conichull([0, 1])) isa Polyhedra.RaysHull{Float64}
            @test convexhull(convexhull([1, 0]), Line([0, 1.])) isa Polyhedra.Hull{Float64}
            @test convexhull(Line([0, 1.]), convexhull([1, 0])) isa Polyhedra.Hull{Float64}
            @test convexhull(convexhull(Line([0, 1.])), [1, 0]) isa Polyhedra.Hull{Float64}
            @test convexhull(convexhull(Line([1, 0])), Line([0, 1.])) isa Polyhedra.LinesHull{Float64}
            @test convexhull(conichull([1, 0.]), Line([0, 1])) isa Polyhedra.RaysHull{Float64}
            @test convexhull(conichull([1, 0.]), [0, 1]) isa Polyhedra.Hull{Float64}
        end
        @testset "H-representation" begin
            #FIXME inference fails @test (@inferred (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1., 0.], 1))) isa Polyhedra.HyperPlanesIntersection{Float64}
            @test ((HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1., 0.], 1))) isa Polyhedra.HyperPlanesIntersection{Float64}
            @test HyperPlane([1, 1], 0) ∩ HalfSpace([1., 0.], 1.) isa Polyhedra.Intersection{Float64}
            #FIXME inference fails @test (@inferred (HalfSpace([1., 0.], 1.) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)))) isa Polyhedra.Intersection{Float64}
            @test ((HalfSpace([1., 0.], 1.) ∩ (HyperPlane([1, 1], 0) ∩ HyperPlane([1, 0], 1)))) isa Polyhedra.Intersection{Float64}
        end
    end
    @testset "Conversion with different array type" begin
        @testset "V-representation" begin
            vv = convexhull(@SVector [0, 1]) + conichull((@SVector [1, 1]), Line(@SVector [1, 0]))
            mv = vrep([0 1], [1 1; 1 0], BitSet([2]))
            generator_fulltest(vv, mv)
            mvv = @inferred convert(typeof(mv), vv)
            vmv = @inferred convert(typeof(vv), mv)
            generator_fulltest(vv, mvv)
            generator_fulltest(vmv, vv)
        end
        @testset "H-representation" begin
            vh = HalfSpace((@SVector [1, 0]), 0) ∩ HyperPlane((@SVector [0, 1]), 1)
            mh = hrep([0 1; 1 0], [1, 0], BitSet(1))
            inequality_fulltest(vh, mh)
            mvh = @inferred convert(typeof(mh), vh)
            vmh = @inferred convert(typeof(vh), mh)
            inequality_fulltest(vh, mvh)
            inequality_fulltest(vmh, vh)
        end
    end
    @testset "Copy test" begin
        @testset "V-representation" begin
            vr = convexhull(@SVector [0, 1])
            vc = copy(vr)
            @test vc !== vr
            @test vc.points !== vr.points
            generator_fulltest(vc, vr)
            vr = conichull(Line(@SVector [1, 1]), Line(@SVector [1, 0]))
            vc = copy(vr)
            @test vc !== vr
            @test vc.lines !== vr.lines
            generator_fulltest(vc, vr)
            vr = conichull(Line(@SVector [1, 1]), Line(@SVector [1, 0]), @SVector [1, 1])
            vc = copy(vr)
            @test vc !== vr
            @test vc.rays !== vr.rays
            @test vc.lines !== vr.lines
            @test vc.lines.lines !== vr.lines.lines
            generator_fulltest(vc, vr)
            vr = convexhull(@SVector [0, 1]) + conichull(Line(@SVector [1, 1]), (@SVector [1, 1]), Line(@SVector [1, 0]))
            vc = copy(vr)
            @test vc !== vr
            @test vc.points !== vr.points
            @test vc.points.points !== vr.points.points
            @test vc.rays !== vr.rays
            @test vc.rays.lines !== vr.rays.lines
            @test vc.rays.lines.lines !== vr.rays.lines.lines
            @test vc.rays.rays !== vr.rays.rays
            generator_fulltest(vc, vr)
        end
        @testset "H-representation" begin
            hr = HyperPlane([1, 0], 1) ∩ HyperPlane([0, 1], -1)
            hc = copy(hr)
            @test hc !== hr
            @test hc.hyperplanes !== hr.hyperplanes
            inequality_fulltest(hc, hr)
            hr = HyperPlane([1, 0], 1) ∩ HyperPlane([1, 1], 0) ∩ HalfSpace([0, 1], -1)
            hc = copy(hr)
            @test hc !== hr
            @test hc.hyperplanes !== hr.hyperplanes
            @test hc.hyperplanes.hyperplanes !== hr.hyperplanes.hyperplanes
            @test hc.halfspaces !== hr.halfspaces
            inequality_fulltest(hc, hr)
        end
    end
    @testset "Preserving sparsity" begin
        for h in (HalfSpace(sparsevec([1], [1], 2), 1) ∩ HyperPlane(sparsevec([2], [-1]), 3), hrep(sparse([1, 2], [1, 2], [1, -1]), [1, 3]))
            @test Polyhedra.Intersection(h) isa Polyhedra.Intersection{Int, SparseVector{Int, Int}}
            @test LPHRep(h) isa LPHRep{Int}
            @test MixedMatHRep(h) isa MixedMatHRep{Int, SparseMatrixCSC{Int, Int}}
            @test LiftedHRepresentation(h) isa LiftedHRepresentation{Int, SparseMatrixCSC{Int, Int}}
            @test Polyhedra.Intersection{Float64}(h) isa Polyhedra.Intersection{Float64, SparseVector{Float64, Int}}
            @test LPHRep{Float64}(h) isa LPHRep{Float64}
            @test MixedMatHRep{Float64}(h) isa MixedMatHRep{Float64, SparseMatrixCSC{Float64, Int}}
            @test LiftedHRepresentation{Float64}(h) isa LiftedHRepresentation{Float64, SparseMatrixCSC{Float64, Int}}
        end
        for v in (convexhull(sparsevec([1], [1], 2), sparsevec([2], [-1])), vrep(sparse([1, 2], [1, 2], [1, -1])))
            @test Polyhedra.Hull(v) isa Polyhedra.Hull{Int, SparseVector{Int, Int}}
            @test MixedMatVRep(v) isa MixedMatVRep{Int, SparseMatrixCSC{Int, Int}}
            @test LiftedVRepresentation(v) isa LiftedVRepresentation{Int, SparseMatrixCSC{Int, Int}}
            @test Polyhedra.Hull{Float64}(v) isa Polyhedra.Hull{Float64, SparseVector{Float64, Int}}
            @test MixedMatVRep{Float64}(v) isa MixedMatVRep{Float64, SparseMatrixCSC{Float64, Int}}
            @test LiftedVRepresentation{Float64}(v) isa LiftedVRepresentation{Float64, SparseMatrixCSC{Float64, Int}}
        end
    end

    @testset "Polar tests" begin
        polar_test()
    end
end


text:
########## Definitions of various stochastic process ##########

abstract type AbstractRandomField end
abstract type AbstractRealRandomField <: AbstractRandomField end

"""
Generic zero-mean real-valued stochastic process.
"""
abstract type StochasticProcess{T<:TimeStyle} <: AbstractRealRandomField end
const ContinuousTimeStochasticProcess = StochasticProcess{ContinuousTime}
const DiscreteTimeStochasticProcess = StochasticProcess{DiscreteTime}

abstract type StationaryProcess{T} <: StochasticProcess{T} end
const ContinuousTimeStationaryProcess = StationaryProcess{ContinuousTime}
const DiscreteTimeStationaryProcess = StationaryProcess{DiscreteTime}

abstract type SelfSimilarProcess <: ContinuousTimeStochasticProcess end
abstract type SSSIProcess <: SelfSimilarProcess end  #

@doc raw"""
Stationary process resulting from filtration (without DC) of another process.

The causal version is defined as:
    \sum_{n=0}^{N-1} a[n+1] X(t-nδ)

and the anti-causal version:
    \sum_{n=0}^{N-1} a[n+1] X(t+nδ)

# Note
- Under assumption of stationarity, the causality has no effect on the auto-covariance function. The implementation of concrete type should follow the anti-causal convention.
- The aimed concrete type of the abstract type `IncrementProcess` is `FractionalGaussianNoise` (fGn) which is stationary. However due to the lack of multiple inheritance in Julia it is very hard to make `IncrementProcess` a subtype of `StationaryProcess`. Possible solutions to this problem include 1) SimpleTrait.jl,  2) copy functions defined for `StationaryProcess`,  3) force `FilteredProcess` to be a subtype of `StationaryProcess`. We adopt the solution 3) here.
"""
abstract type FilteredProcess{T<:TimeStyle, P<:StochasticProcess{>:T}} <: StationaryProcess{T} end

"""
Differential process `X(t±δ) - X(t)`
"""
abstract type DifferentialProcess{T<:TimeStyle, P<:StochasticProcess{>:T}} <: FilteredProcess{T, P} end  # Process as first order finite difference of another process.

"""
Discrete time differential process `X((n±l)δ) - X(nδ)` with `l` being the lag.
"""
abstract type IncrementProcess{P<:StochasticProcess} <: DifferentialProcess{DiscreteTime, P} end

#### Process specific functions ####
"""
Return the self-similar exponent of the process.
"""
ss_exponent(X::SelfSimilarProcess) = throw(NotImplementedError())

"""
Return the filter of the process.
"""
filter(X::FilteredProcess)::AbstractVector = throw(NotImplementedError())

"""
Test whether the filtration is causal:
"""
iscausal(X::FilteredProcess) = throw(NotImplementedError())

"""
Return the parent process of the filtered process.
"""
parent_process(X::FilteredProcess) = throw(NotImplementedError())

"""
Return the step `δ` of the filtered process.
"""
step(X::FilteredProcess)::Real = throw(NotImplementedError())

"""
Filter of differential process, defined as
- the causal case: X(t) - X(t-lδ)
- the anti-causal case: X(t+lδ) - X(t)
with `l = lag(X)`.
"""
function filter(X::DifferentialProcess)
    filt = vcat(1, zeros(lag(X)-1), -1)
    return iscausal(X) ? filt : reverse(filt)
end

"""
Return the lag of differential process.
"""
lag(X::DifferentialProcess) = throw(NotImplementedError())

#### Generic identifiers ####
"""
Test whether a process is time discrete (a time series).
"""
iscontinuoustime(X::StochasticProcess{ContinuousTime}) = true
iscontinuoustime(X::StochasticProcess{DiscreteTime}) = false

"""
Test whether a process is multivariate.
"""
ismultivariate(X::StochasticProcess) = false

"""
Test whether a process has stationary increments.
"""
isincrementstationary(X::StochasticProcess) = false
isincrementstationary(X::SSSIProcess) = true
isincrementstationary(X::StationaryProcess) = true

"""
Test whether a process is stationary.
"""
isstationary(X::StochasticProcess) = false
isstationary(X::StationaryProcess) = true
# isstationary(X::FilteredProcess{T, P}) where {T, P<:StationaryProcess} = true
# isstationary(X::IncrementProcess{P}) where {P<:SSSIProcess} = true

"""
Determine whether a grid has the constant step.
"""
function isregulargrid(G::AbstractVector)
    return if 1 <= length(G) <= 2
        true
    elseif length(G) > 2  # second order difference should be close to 0
        isapprox(maximum(abs.(diff(diff(G)))), 0.0; atol=1e-10)
        # ddG = unique(diff(diff(G)))
        # length(ddG) == 1 && ddG[1] == 0.
    else
        false
    end
end

#### Statistics for stochastic process ####
"""
Auto-covariance function of a stochastic process.
"""
autocov(X::StochasticProcess{T}, t::T, s::T) where T<:TimeStyle = throw(NotImplementedError())

autocov(X::StationaryProcess{T}, t::T) where T<:TimeStyle = throw(NotImplementedError())
autocov(X::StationaryProcess{T}, t::T, s::T) where T<:TimeStyle = autocov(X, t-s)

"""
Compute the auto-covariance matrix of a stochastic process on a sampling grid.
"""
function autocov!(C::Matrix{<:Real}, X::StochasticProcess{T}, G::AbstractVector{<:T}) where T<:TimeStyle
    @assert size(C, 1) == size(C, 2) == length(G)

    # construct the covariance matrix (a symmetric matrix)
    N = size(C, 1)  # dimension of the auto-covariance matrix
    for c = 1:N, r = 1:c
        C[r,c] = autocov(X, G[r], G[c])
    end
    for c = 1:N, r = (c+1):N
        C[r,c] = C[c,r]
    end
    return Symmetric(C)
end

function autocov!(C::Matrix{<:Real}, X::StochasticProcess{T}, G1::AbstractVector{<:T}, G2::AbstractVector{<:T}) where T<:TimeStyle
    @assert size(C, 1) == length(G1) && size(C, 2) == length(G2)

    # construct the covariance matrix (a symmetric matrix)
    N,M = size(C)  # dimension of the auto-covariance matrix
    for c = 1:M, r = 1:N
        C[r,c] = autocov(X, G1[r], G2[c])
    end
    return C
end

"""
Compute the auto-covariance sequence of a stationary process on a regular grid.
"""
function autocov!(C::AbstractVector{<:Real}, X::StationaryProcess{T}, G::AbstractVector{<:T}) where T<:TimeStyle
    # @assert isregulargrid(G)
    @assert length(C) == length(G)  # check dimension

    # construct the auto-covariance kernel
    for n = 1:length(C)
        C[n] = autocov(X, G[n]-G[1])
    end
    return C
end

"""
    covseq(X::StationaryProcess, G::RegularGrid)

Return the auto-covariance sequence of a stationary process on a regular grid.
"""
covseq(X::StationaryProcess{T}, G::AbstractVector{<:T}) where T<:TimeStyle = autocov!(zeros(length(G)), X, G)
covseq(X::StationaryProcess, N::Integer) = covseq(X, 1:N)

function autocov!(C::Matrix{<:Real}, X::StationaryProcess{T}, G::AbstractVector{<:T}) where T<:TimeStyle
    # check dimension
    @assert size(C, 1) == size(C, 2) == length(G)
    # println("autocov! of StationaryProcess, $(ss_exponent(X))")

    return if isregulargrid(G)
        # construct the covariance matrix (a Toeplitz matrix)
        covmat!(C, covseq(X,G))
    else
        # if G is not regular the `covseq` can not be applied, invoke the function of `StochasticProcess`.
        invoke(autocov!, Tuple{Matrix{<:Real}, StochasticProcess{S}, AbstractVector{<:S}} where S<:TimeStyle, C, X, G)
    end
end

"""
Construct the covariance matrix from the covariance sequence.
"""
function covmat!(C::Matrix{T}, S::AbstractVector{T}) where {T<:Real}
    for c = 1:length(S), r = 1:length(S)
        C[r,c] = S[abs(r-c)+1]
    end
    return C
end

covmat(S::AbstractVector{T}) where {T<:Real} = covmat!(zeros(T, length(S), length(S)), S)

"""
Return the auto-covariance matrix of a stochastic process on a sampling grid.
"""
covmat(X::StochasticProcess, G::AbstractVector) = autocov!(zeros(length(G), length(G)), X, G)

"""
Return the auto-covariance matrix of a stochastic process between two sampling grids.

The `(i,j)`-th coefficient in the matrix is `autocov(G1[i], G2[j])`.
"""
covmat(X::StochasticProcess, G1::AbstractVector, G2::AbstractVector) = autocov!(zeros(length(G1), length(G2)), X, G1, G2)

"""
Return the auto-covariance matrix on an integer sampling grid `1:N`.
"""
covmat(X::StochasticProcess, N::Integer) = covmat(X, 1:N)
covmat(X::StochasticProcess, N::Integer, M::Integer) = covmat(X, 1:N, 1:M)

"""
Return the partial correlation function of a discrete time stationary process.
"""
partcorr(X::DiscreteTimeStationaryProcess, n::DiscreteTime) = throw(NotImplementedError())

"""
Compute the partial correlation sequence of a discrete time stationary process on a regular integer sampling grid.
"""
function partcorr!(C::Vector{<:Real}, X::DiscreteTimeStationaryProcess, G::AbstractVector{<:DiscreteTime})
    # check dimension
    @assert length(C) == length(G)
    @assert isregulargrid(G)

    for n = 1:length(C)
        C[n] = partcorr(X, G[n]-G[1]+1)
    end
    return C
end

"""
Return the partial correlation function of a time discrete stationary process.

# Args
- method: if set to `:LevinsonDurbin` it will use the Levinson-Durbin method which needs only the autocovariance expression of the process.
"""
function partcorr(X::DiscreteTimeStationaryProcess, G::AbstractVector{<:DiscreteTime}, method::Symbol=:None)
    if method == :LevinsonDurbin # use Levinson-Durbin
        cseq = covseq(X, G)
        pseq, sseq, rseq = LevinsonDurbin(cseq)
        return rseq
    else
        return partcorr!(zeros(length(G)), X, G)
    end
end

#### Statistical inference on stochastic process ####

"""
    cond_mean_cov(P::StochasticProcess{T}, Gx::AbstractVector{<:T}, Gy::AbstractVector{<:T}, Y::AbstractVector{<:Real}) where T<:TimeStyle

Conditional mean and covariance of a zero-mean Gaussian process `P` on the position `Gx` given the value `Y` on the position `Gy`.
"""
function cond_mean_cov(P::StochasticProcess{T}, Gx::AbstractVector{<:T}, Gy::AbstractVector{<:T}, Y::AbstractVector{<:Real}) where T<:TimeStyle
    @assert length(Gy) == length(Y)

    Σxx = covmat(P, Gx)
    Σxy = covmat(P, Gx, Gy)
    Σyy = covmat(P, Gy)
    # inverse of Σyy, method 1: pseudo inverse
    iΣyy = pinv(Matrix(Σyy))
    # # method 2: LU
    # iL = inv(cholesky(Σyy).L)
    # iΣyy = iL' * iL

    μc = Σxy * iΣyy * Y

    # μc = Σxy * iΣyy * Y
    # μc = Σxy * (Σyy\Y)
    Σc = Σxx - Σxy * iΣyy * Σxy'

    return (μ=μc, Σ=Σc, C=Σxy * iΣyy)
end

cond_mean_cov(P::StochasticProcess, gx::ContinuousTime, Gy::AbstractVector, Y::AbstractVector) = cond_mean_cov(P, [gx], Gy, Y)

"""
Conditional mean and covariance on regular grid.
"""
cond_mean_cov(P::StochasticProcess, n::DiscreteTime, Y::AbstractVector) = cond_mean_cov(P, (1:n).+length(Y), 1:length(Y), Y)

"""
EXPERIMENTAL: linear coefficient of recursive conditional mean

Predict `k` steps into the future.
"""
function cond_mean_coeff(P::StochasticProcess{T}, k::Integer, l::Integer; mode::Symbol=:recursive) where T<:TimeStyle
    @assert k>0

    Σyy = covmat(P, 0:l-1)

    if mode == :recursive
        Σxy = covmat(P, [l], 0:l-1)
        cv = Σxy * pinv(Matrix(Σyy))
        M = vcat(diagm(1 => ones(l-1))[1:l-1,:], cv)
        Cv = [cv]

        for i=2:k
            push!(Cv, Cv[end]*M)
        end

        return vcat(Cv...)
    else
        Cv = covmat(P, l:l+k-1, 0:l-1) * pinv(Matrix(Σyy))
        return Cv
    end
end

include("FBM.jl")  # Fractional Brownian Motion related
include("MFBM.jl")  # Multi-Fractional Brownian Motion related
include("FARIMA.jl")  # Fractional ARIMA related

ext: jl

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
import Base

"""
    TrenchObject
An abstract type that is the root of an object as implemented in
Specifically, see also `subtypes(foo).`

```julia-repl
julia> subtypes(TrenchObject)
2-element Vector{Any}:
 DiffableFunction
 Variable
```
"""
abstract type TrenchObject end

"""
    DiffableFunction
An abstract type that is the parent type of differentiable functions;

all function in JITrench must be children of this type.

# Examples

```julia-repl
julia> subtypes(DiffableFunction)
24-element Vector{Any}:
 JITrench.Add
 JITrench.BroadcastTo
 JITrench.Broadcasting
 JITrench.Cos
 JITrench.Div
 JITrench.Exp
 JITrench.Flatten
 JITrench.GetIndex
 JITrench.GetIndexGrad
 JITrench.Log
 ⋮
 JITrench.Reshape
 JITrench.Sigmoid
 JITrench.Sin
 JITrench.Sub
 JITrench.Sum
 JITrench.SumTo
 JITrench.Tan
 JITrench.Transpose
 JITrench._Linear
```
"""
abstract type DiffableFunction  <: TrenchObject end

"""
    Variable(values, [creator, grad, generation, name])

A type that is a variable in JITrench.

For a detailed explanation of this, see the documentation(here).

Only real numbers or arrays consisting of real numbers can be stored.

# Examples
```julia-repl
julia> Variable(2)
name: nothing
values: 2
creator: User-Defined(nothing)

julia> Variable([1.2, 3.5])
name: nothing
values: [1.2, 3.5]
creator: User-Defined(nothing)
```

"""
mutable struct Variable{T} <: TrenchObject
    values::T
    creator
    grad::Union{Variable,Nothing}
    generation::Int
    name::Union{Nothing,String}
    function Variable(values::T; creator=nothing, grad=nothing, generation=0, name=nothing) where {T <: Union{<:Real,AbstractArray{<:Real}}}
        new{T}(values, creator, grad, generation, name)
    end
end

Base.promote_rule(::Type{<:Real}, ::Type{<:Variable}) = Variable

Base.convert(::Type{Variable}, x::AbstractArray) = Variable(x)

Base.convert(::Type{Variable}, x::Real) = Variable(x)

Base.convert(::Type{Variable}, x::Variable) = x

Base.size(x::Variable) = size(x.values)

function get_output_str(var::Variable)
    output = ""
    output *= "name: $(var.name) \n"
    output *= "values: $(var.values)\n"
    if (var.grad !== nothing)
        output *= "grad: $(var.grad)\n"
    end
    if (var.creator !== nothing)
        output *= "creator: $(typeof(var.creator))"
    else
        output *= "creator: User-Defined(nothing)"
    end
    return output
end

# print()
function Base.show(io::IO, var::Variable)
    print(io, "Variable($(var.values))")
end

# REPL
function Base.show(io::IO, ::MIME"text/plain", var::Variable)
    print(io, get_output_str(var))
end


text:
using Arena.Amphitheatre
#Extending the functionallity of Amphitheatre requires 2 things
# 1. A type definition of the object to draw inhereting from AbstractAmphitheatre, eg. DrawSphereAmphi <: AbstractAmphitheatre
# 2. A visualize! function of the form `function visualize!(vis::Visualizer, object_to_draw::DrawSphereAmphi)::Nothing` to do the drawing
import Arena.Amphitheatre: visualize!

# using some usefull drawing functions
using MeshCat, GeometryTypes, CoordinateTransformations
# ============================================================
# --------------------DrawSphereInMeshCat--------------------
# ============================================================

# 1. the object to draw
struct DrawSphereAmphi <: AbstractAmphitheatre
	robotId::String
	sessionId::String
	spheres::Dict{Symbol, Vector{Float64}}
	isdrawn::Dict{Symbol, Bool}
	sphereScale::Float64
end
# a helpfull constructor
DrawSphereAmphi(robotId::String, sessionId::String, spheres::Dict{Symbol, Vector{Float64}}, sphereScale::Float64 = 0.1) =
	DrawSphereAmphi(robotId, sessionId, spheres, Dict{Symbol, Bool}(), sphereScale)

# 2. The visualize! function that does the drawing
function visualize!(vis::Visualizer, sa::DrawSphereAmphi)::Nothing

	for key in keys(sa.spheres)

		if !get(sa.isdrawn, key, false)
			haskey(sa.isdrawn, key) ? sa.isdrawn[key] = true : push!(sa.isdrawn, key=>true)

			sphere = HyperSphere(Point(0., 0, 0), sa.sphereScale)
			matcolor = MeshPhongMaterial(color=RGBA(0,1,0,0.5))
			setobject!(vis[sa.robotId][sa.sessionId][key], sphere, matcolor)
		end
		tform = Translation(sa.spheres[key][1:3]...)
		settransform!(vis[sa.robotId][sa.sessionId][key], tform)
	end
end

## thats it, now lets test or new visualizer

#create the object
balVis = DrawSphereAmphi("robot", "session1", Dict(:s0=>[0.,0,0], :s1=>[1.,0,0]))

# add it to the object to draw array of type AbstractAmphitheatre
visdatasets = AbstractAmphitheatre[balVis]

#start the visualizer
vis, vistask = visualize(visdatasets, start_browser=true)

#now push in a few more spheres and see them pop up
push!(balVis.spheres, :s2=>[2.,0,0])
push!(balVis.spheres, :s3=>[2.,1,1])
push!(balVis.spheres, :s4=>[1.,1,1])

# you can easily create a second one
push!(visdatasets, DrawSphereAmphi("robot", "session2", Dict(:s0=>[0.,-1,0], :s1=>[1.,-1,0]), 0.2))

## finally run stopAmphiVis!() to stop the visualizer
stopAmphiVis!()


text:
function solve()
    n, m = [parse(Int, x) for x in split(readline())]
    as = [parse(Int, x) for x in split(readline())]
    edges = Dict{Int, Set{Int}}()
    for i=1:n
        edges[i] = Set{Int}()
    end
    for _=1:m
        x, y = [parse(Int, x) for x in split(readline())]
        push!(edges[x], y)
    end

    dp = ones(Int, n) * 10^10
    for i=1:n
        for j in edges[i]
            dp[j] = min(dp[j], dp[i], as[i])
        end
    end
    maximum((as - dp)[2:end])
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text:
using DataStructures, SCS, Mosek, OPFInstances
using MathProgComplex

using JuMP, KNITRO

!isdefined(:MPC) && (const MPC = MathProgComplex)

function main()
    problem_c, pt = MPC.import_from_dat(getinstancepath("Matpower", "QCQP", "WB5"))
    problem = MPC.pb_cplx2real(problem_c)

    workpath = joinpath("Mosek_runs", "worksdp")
    ispath(workpath) && rm(workpath, recursive=true)
    mkpath(workpath)

    # cstobj = problem.objective[MPC.Exponent()]
    # cstobj = 0
    # problem.objective += -cstobj
    # @assert !haskey(problem.objective, MPC.Exponent())

    di = Dict("LOAD_2_Im" => 1,
            "LOAD_2_Re" => 2,
            "LOAD_3_Im" => 1,
            "LOAD_3_Re" => 2,
            "LOAD_4_Im" => 1,
            "LOAD_4_Re" => 1,
            "UNIT_1_Im" => 1,
            "UNIT_1_Re" => 1,
            "UNIT_5_Im" => 1,
            "UNIT_5_Re" => 1,
            "VOLTM_1_Re" => 1,
            "VOLTM_2_Re" => 1,
            "VOLTM_3_Re" => 1,
            "VOLTM_4_Re" => 1,
            "VOLTM_5_Re" => 1)

    relax_ctx = MPC.set_relaxation(problem; hierarchykind=:Real,
                                            issparse = false,
                                            # d = 2,
                                            di = di,
                                            params = Dict(:opt_outlev=>1,
                                                          :opt_outmode=>0,
                                                          :opt_pbsolved=>:MomentRelaxation,
                                                          :opt_solver=>:MosekCAPI))

    primal = SortedDict{Tuple{String,String,String}, Float64}()
    dual = SortedDict{Tuple{String, String, String}, Float64}()

    primobj_JuMP, dualobj_JuMP = MPC.run_hierarchy(problem, relax_ctx, indentedprint=true, save_pbs=true,
                                                                        primsol = primal,
                                                                        dualsol = dual)

    varname2ind = Dict([varname=>ind for (ind, varname) in enumerate(keys(problem.variables))])
    ind2varname = Dict([ind=>varname for (varname, ind) in varname2ind])

    n = length(varname2ind)

    mat = Matrix{Float64}(n, n)

    for ((clique, var1, var2), val) in primal
        if haskey(varname2ind, var1) && haskey(varname2ind, var2)
            mat[varname2ind[var1], varname2ind[var2]] = val
            mat[varname2ind[var2], varname2ind[var1]] = val
        end
    end

    eigenvals, eigenvecs = eig(mat)

    # mat2 = zeros(n, n)
    # for i=1:n
    #     mat2 += eigenvals[i] * (eigenvecs[: , i] * transpose(eigenvecs[: , i]))
    # end

    # println("\nEigenvals:")
    # display(eigenvals)

    # println("\nWorking with:")
    maxvec = eigenvecs[:, n]
    display(maxvec)

    ## Deal with sign
    if count(x->x<=0, maxvec) > n/2
        maxvec *= -1
    end

    ## Find the closest corresponding measure
    point_maxvec = Point()
    for i=1:n
        varname = ind2varname[i]
        vartype = problem.variables[varname]
        setindex!(point_maxvec, maxvec[i], Variable(varname, vartype))
    end

    localpb = Problem()
    localpb.variables = problem.variables

    obj = Polynomial()
    for ((clique, var1, var2), val) in primal
        if clique == "clique1" && var2 <= var1
            add!(obj, (parsevar(var1) * parsevar(var2) - val)^2)
        end
    end

    set_objective!(localpb, obj)

    solver = KnitroSolver()
    m, JuMPvar = get_JuMP_cartesian_model(localpb, solver)
    solve(m)

    pointknitro = Point()
    for (varname, jumpvar) in JuMPvar
        println("$jumpvar       ", getvalue(jumpvar))
        setindex!(pointknitro, getvalue(jumpvar), Variable(varname, Real))
    end

    println("--------------------------------------")
    println("\nslacks points - maxvec:\n", get_slacks(problem, point_maxvec))
    println("\nslacks points - knitro:\n", get_slacks(problem, pointknitro))

    return nothing
end

function parsevar(var::String)
    if var == "1"
        return 1.0
    elseif ismatch(r"\AVOLT_\d_(Re|Im)\Z", var)
        return Variable(var, Real)
    elseif ismatch(r"\AVOLT_\d_(Re|Im)\*VOLT_\d_(Re|Im)\Z", var)
        var1, var2 = matchall(r"VOLT_\d_(Re|Im)", var)
        return Variable(var1, Real) * Variable(var2, Real)
    elseif ismatch(r"\AVOLT_\d_(Re|Im)\^2\Z", var)
        var1 = matchall(r"VOLT_\d_(Re|Im)", var)
        return Variable(string(var1), Real)^2
    else
        @show var
        error()
    end
end

main()


text:
module LULESH

using CEnum
using CUDA
using MPI
using Parameters
using Printf
using Random

abstract type AbstractDomain end
# Device vector types
const VD{T} = Union{Vector{T}, CuVector{T}} where T

# Index type
const IndexT = Int
const MAX_FIELDS_PER_MPI_COMM = 6

struct LuleshProblem
    num_iters::Int
    structured::Bool
    nx::Int
    nr::Int
    balance::Int
    cost::Int
    col::Int
    row::Int
    plane::Int
    side::Int
    devicetype
    floattype
    comm::Union{MPI.Comm, Nothing}
    function LuleshProblem(num_iters, structured, nx, nr, balance, cost, devicetype, floattype, comm)
        !MPI.Initialized() && MPI.Init()
        col, row, plane, side = InitMeshDecomp(comm)
        return new(num_iters, structured, nx, nr, balance, cost, col, row, plane, side, devicetype, floattype, comm)
    end
end

include("bc.jl")
include("types.jl")
include("domain.jl")
include("mpi.jl")
include("utils.jl")

export printUsage, IndexT, Domain, LuleshProblem, getMyRank, getNumRanks, getWtime, lagrangeLeapFrog, comm_max

function InitMeshDecomp(comm)
    # Assume cube processor layout for now
    numRanks = getNumRanks(comm)
    myRank = getMyRank(comm)
    testProcs = floor(cbrt(numRanks+0.5))
    @show testProcs, numRanks
    if (testProcs*testProcs*testProcs != numRanks)
        error("Num processors must be a cube of an integer (1, 8, 27, ...)")
    end
    @show typeof(testProcs)
    # TODO This is not good for the padding
    if (MAX_FIELDS_PER_MPI_COMM > getCacheCoherencePad(testProcs))
        error("corner element comm buffers too small. MAX_FIELDS_PER_MPI_COMM > CACHE_COHERENCE_PAD_REAL ($MAX_FIELDS_PER_MPI_COMM > $(getCacheCoherencePad(testProcs))")
    end
    testProcs = convert(Int, testProcs)
    dx = convert(Int, testProcs)
    dy = convert(Int, testProcs)
    dz = convert(Int, testProcs)

    # temporary test
    if dx*dy*dz != numRanks
        error("error -- must have as many domains as procs\n") ;
    end
    remainder = (dx*dy*dz) % numRanks
    myDom = Int(0)
    if (myRank < remainder)
        myDom = myRank*( 1+ (div(dx*dy*dz, numRanks)))
    else
      myDom = remainder*( 1+ (div(dx*dy*dz, numRanks))) +
         (myRank - remainder)*(div(dx*dy*dz, numRanks))
    end
   @show typeof(myDom)
   col = myDom % dx
   row = div(myDom, dx) % dy
   plane = div(myDom, (dx*dy))
   side = testProcs

   return col, row, plane, side
end
export InitMeshDecomp
end # module


text:
module Err
export PatternCompilationError,
    InternalException, SyntaxError, UnknownExtension, PatternUnsolvedException

PatternUnsolvedException = ErrorException

struct PatternCompilationError <: Exception
    line::Union{LineNumberNode, Nothing}
    ex::Exception
end

struct InternalException <: Exception
    msg::String
end

struct SyntaxError <: Exception
    msg::String
end

struct UnknownExtension <: Exception
    ext::Union{String, Symbol}
end

end


text:
@testset "precision" begin
    precdefault = 256
    prec = 64

    @testset "precision: $T" for T in
                                 (Arf, Arb, Acb, ArbPoly, ArbSeries, AcbPoly, AcbSeries)
        @test precision(T()) == precdefault
        @test precision(T(prec = prec)) == prec

        @test precision(Arblib.cstruct(T())) == precdefault
        @test precision(Arblib.cstruct(T(prec = prec))) == precdefault

        @test precision(T) == precdefault
        @test precision(Arblib.cstructtype(T)) == precdefault
        @test precision(Ptr{Arblib.cstructtype(T)}) == precdefault

        if T in [ArbPoly, ArbSeries, AcbPoly, AcbSeries]
            @test precision(T()[0]) == precdefault
            @test precision(T(prec = prec)[0]) == prec
        end
    end

    @testset "setprecision: $T" for T in
                                    (Arf, Arb, Acb, ArbPoly, ArbSeries, AcbPoly, AcbSeries)
        x = T()
        @test precision(x) == precdefault

        y = setprecision(x, prec)
        @test precision(y) == prec
        @test isequal(x, y)

        if T in [Arf, Arb, Acb]
            Arblib.set!(y, 2)
            @test !isequal(x, y)
        elseif T in [ArbPoly, ArbSeries, AcbPoly, AcbSeries]
            @test precision(y[0]) == prec
            y[0] = eltype(T)(2)
            @test !isequal(x, y)
        end
    end

    @testset "setprecision 2" begin
        x = Arb("0.1")
        @test precision(x) == 256
        @test precision(Arb) == 256
        @test string(x) ==
              "[0.1000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 1.95e-78]"

        setprecision(Arb, 64)
        @test precision(x) == 256
        @test precision(Arb) == 64
        y = Arb("0.1")
        @test precision(y) == 64
        @test string(y) == "[0.100000000000000000 +/- 1.22e-20]"

        setprecision(Arb, 256)
    end

    @testset "precision change: $T" for T in
                                        [ArbMatrix, ArbRefMatrix, AcbMatrix, AcbRefMatrix]
        A = T(3, 3; prec = 96)
        @test precision(A) == 96
        @test precision(A[1, 1]) == 96
        B = setprecision(A, 128)
        @test precision(B) == 128
        @test precision(B[1, 1]) == 128
        A[1, 1][] = 1
        @test A[1, 1] == B[1, 1]
    end

    @testset "_precision" begin
        let _precision = Arblib._precision
            @test _precision(Arf(prec = 64)) == 64
            @test _precision(Arb(prec = 64)) == 64
            @test _precision(Acb(prec = 64)) == 64
            @test _precision(BigFloat(precision = 64)) == 64
            @test _precision(Arf(prec = 64) + im * Arf(prec = 80)) == 80
            @test _precision(Arb(prec = 80) + im * Arb(prec = 64)) == 80
            @test _precision(Complex(BigFloat(precision = 64), BigFloat(precision = 80))) ==
                  80

            @test _precision(1) == precision(Arb)
            @test _precision(1.0) == precision(Arb)
            @test _precision(1 // 1) == precision(Arb)
            @test _precision(1 + im) == precision(Arb)
            @test _precision(1.0 + im) == precision(Arb)
            @test _precision(1 // 1 + im) == precision(Arb)

            @test _precision((Arf(prec = 64), Arb(prec = 80))) == 80
            @test _precision((Arf(prec = 80), Arb(prec = 64))) == 80
            @test _precision((BigFloat(precision = 80), Acb(prec = 64))) == 80
            @test _precision((BigFloat(precision = 64), Acb(prec = 80))) == 80
            @test _precision((im * Arb(prec = 80), Acb(prec = 64))) == 80
            @test _precision((im * Arb(prec = 64), Acb(prec = 80))) == 80
            @test _precision((
                Complex(BigFloat(precision = 64), BigFloat(precision = 64)),
                Arf(prec = 80),
            )) == 80
            @test _precision((
                Arf(prec = 80),
                Complex(BigFloat(precision = 64), BigFloat(precision = 64)),
            )) == 80

            @test _precision((Arf(prec = 64), 1)) == 64
            @test _precision((1, Arf(prec = 64))) == 64
            @test _precision((BigFloat(precision = 64), 1)) == 64
            @test _precision((1, BigFloat(precision = 64))) == 64
            @test _precision((im * Arf(prec = 64), 1)) == 64
            @test _precision((1, im * Arf(prec = 64))) == 64

            @test _precision((1.0, 2)) == precision(Arb)
        end
    end
end


text:
#------------------------------------------------------------------
# helper methods for the executor
#------------------------------------------------------------------
@timetrack function wait_for_executors(runenv::RunEnv)
    @sync for node in runenv.nodes
        @async begin
            empty!(node.last_task_stat)
            if node.broker_task !== nothing
                push!(node.last_task_stat, (node.brokerid, node.broker_task))
            end
            append!(node.last_task_stat, zip(node.executorids, node.executor_tasks))
            empty!(node.executor_tasks)
            node.broker_task = nothing
            for (pid,task) in node.last_task_stat
                isready(task) || wait(task)
            end
        end
    end

    delete_meta(runenv.masterid, runenv.nodes, runenv.rootpath)
    remotetrack_end(runenv.remotetrack)
    profile_end(runenv.profile, "master_$(runenv.masterid)")
    nothing
end

function print_stats(runenv::RunEnv)
    info("execution stats:")
    for node in runenv.nodes
        map((x)->begin
            pid = x[1]
            result = fetch(x[2])
            role = (pid === node.brokerid) ? "broker" : "executor"
            if isa(result, Tuple)
                if role == "executor"
                    stole, shared, executed = fetch(x[2])
                    info("  executor $pid stole $stole, shared $shared, executed $executed tasks")
                else
                    stole, shared = fetch(x[2])
                    info("broker $pid stole $stole, shared $shared tasks")
                end
            else
                info("$role $pid: $result")
            end
        end, node.last_task_stat)
    end
    nothing
end

function handle_failure(runenv::RunEnv, env::ExecutionCtx, pid, cause)
    meta = env.meta
    isa(cause, ProcessExitedException) && (cause = :killed)
    info("PID $pid failed with $cause")

    for node in runenv.nodes
        executors = node.executorids
        if (pid == node.brokerid) || (pid in executors)
            # remove node from runenv
            filter!((x)->(x!==node), runenv.nodes)
            tasks_to_reschedule = detach(meta, node.brokerid)
            info("detached brokerid ", node.brokerid, " with ", length(tasks_to_reschedule), " tasks to reschedule")
            cleanupnode = ((cause === :killed) && (pid === node.brokerid)) ? first(executors) : node.brokerid

            # identify a cleanup node and terminate the others
            nodestoterminate = convert(Vector{Int}, copy(executors))
            push!(nodestoterminate, node.brokerid)
            (cause === :killed) && filter!((x)->x!=pid, nodestoterminate)
            filter!((x)->x!=cleanupnode, nodestoterminate)
            info("terminating nodes ", nodestoterminate)
            for tpid in nodestoterminate
                try
                    rmprocs(tpid)
                catch ex
                    warn("error $ex terminating $tpid")
                end
            end

            # cleanup and terminate the last (cleanup) node
            info("cleaning up and terminating last node ", cleanupnode)
            try
                remotecall_wait(cleanup_meta, cleanupnode, META_IMPL[:node], runenv.rootpath, node.brokerid)
                rmprocs(Int(cleanupnode))
            catch ex
                warn("error $ex terminating $cleanupnode")
            end
            isempty(runenv.nodes) && error("PID $pid failed with $cause")

            # reschedule tasks being executed by the terminated node
            for (tid,annotatedtid) in tasks_to_reschedule
                info("rescheduling task $tid annotated as $annotatedtid")
                enqueue(env, annotatedtid, false, true)
            end
            break
        end
    end

    nothing
end

function check_failures(runenv::RunEnv, env::ExecutionCtx)
    nfailures = 0

    # check for process failure
    procsnow = sort!(procs())
    if hash(procsnow) !== runenv.nodehash
        usedprocs = [runenv.masterid]
        for node in runenv.nodes
            push!(usedprocs, node.brokerid)
            append!(usedprocs, node.executorids)
        end
        deadprocs = setdiff(usedprocs, procsnow)

        for pid in deadprocs
            handle_failure(runenv, env, pid, :killed)
            nfailures += 1
        end
    end

    @sync for node in runenv.nodes
        @async begin
            tasklist = []
            if node.broker_task !== nothing
                push!(tasklist, (node.brokerid, node.broker_task))
            end
            append!(tasklist, zip(node.executorids, node.executor_tasks))
            for (pid,task) in tasklist
                try
                    if isready(task)
                        result = fetch(task)
                        if !isa(result, Tuple)
                            handle_failure(runenv, env, pid, result)
                            nfailures += 1
                        end
                    end
                catch ex
                    handle_failure(runenv, env, pid, ex)
                    nfailures += 1
                end
            end
        end
    end
    nfailures
end

function delete_meta(masterid::UInt64, nodes::Vector{NodeEnv}, rootpath::String)
    @sync begin
        @async delete_meta(META_IMPL[:cluster], rootpath, masterid)
        for node in nodes
            if node.brokerid !== masterid
                @async remotecall_wait(delete_meta, node.brokerid, META_IMPL[:node], rootpath, node.brokerid)
            end
        end
    end
    nothing
end

function delete_meta(metaimpl::String, rootpath::String, brokerid::Integer)
    broker_rootpath = joinpath(rootpath, string(brokerid))
    M = metastore(metaimpl, broker_rootpath, 0)
    M.brokerid = myid()
    delete!(M)
    nothing
end

function cleanup_meta(masterid::UInt64, nodes::Vector{NodeEnv}, rootpath::String)
    @sync begin
        @async cleanup_meta(META_IMPL[:cluster], rootpath, masterid)
        for node in nodes
            if node.brokerid !== masterid
                @async remotecall_wait(cleanup_meta, node.brokerid, META_IMPL[:node], rootpath, node.brokerid)
            end
        end
    end
    nothing
end

function cleanup_meta(metaimpl::String, rootpath::String, brokerid::Integer)
    broker_rootpath = joinpath(rootpath, string(brokerid))
    M = metastore(metaimpl, broker_rootpath, 0)
    M.brokerid = brokerid
    DagScheduler.cleanup(M)
    nothing
end

#------------------------------------------------------------------
# per process execution engine context
#------------------------------------------------------------------
const genv = Ref{Union{ExecutionCtx,Void}}(nothing)
const upstream_genv = Ref{Union{ExecutionCtx,Void}}(nothing)
const join_count = Ref(0)
const start_cond = Channel{Void}(1)

#------------------------------------------------------------------
# helper methods for the broker
#------------------------------------------------------------------
function join_cluster()
    join_count[] -= 1
    if join_count[] == 0
        put!(start_cond, nothing)
    end
end

function upstream_result_to_local(env, task, res)
    try
        if !has_result(env.meta, task)
            t = get_executable(env, task)
            set_result(env.meta, task, res; processlocal=false)
        end
    catch ex
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    end
    nothing
end

function local_result_to_upstream(upenv, task, res)
    try
        if ((task === taskid(upenv.dag_root)) || was_stolen(upenv, task)) && !has_result(upenv.meta, task)
            t = get_executable(upenv, task)
            set_result(upenv.meta, task, res; processlocal=false)
        end
    catch ex
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    end
    nothing
end

function unify_trigger(env::ExecutionCtx, root::TaskIdType, unified_trigger::Channel{Void}, do_trigger::Ref{Bool})
    while !has_result(env.meta, root)
        wait_trigger(env.meta)
        isready(unified_trigger) || put!(unified_trigger, nothing)
    end
    do_trigger[] = false
    put!(unified_trigger, nothing)
    nothing
end

function broker_task_selector(tasks::Vector{UInt64}, brokerid::UInt64, rev::Bool)
    env = (genv[])::ExecutionCtx
    costs = env.costs
    del_annotation = env.meta.del_annotation
    affs = map((tid)->affinity(del_annotation(tid), brokerid, costs), tasks)
    _val, idx = rev ? findmin(affs) : findmax(affs)
    idx
end

function broker_tasks(upenv::ExecutionCtx, env::ExecutionCtx, unified_trigger::Channel{Void}, do_trigger::Ref{Bool})
    try
        costs = env.costs
        hascosts = (costs !== nothing) && !isempty(costs)
        brokerid = env.brokerid
        broker_from_upstream_task_selector = (tasks) -> broker_task_selector(tasks, brokerid, false)  # get task with highest affinity for the broker
        broker_from_downstream_task_selector = (tasks) -> broker_task_selector(tasks, brokerid, true) # get task with least affinity for the broker
        last_upstream = last_downstream = 0.0

        while do_trigger[]
            #tasklog(env, "doing broker tasks")
            tnow = time()

            upstat, downstat = upenv.meta.sharemode, env.meta.sharemode
            needed_upstream, needed_downstream = should_share(upenv), should_share(env)
            have_upstream = (upstat.ncreated > 0) ? (upstat.nshared > 0) : true
            have_downstream = (downstat.ncreated > 0) ? (downstat.nshared > 0) : true
            from_upstream = needed_downstream && have_upstream && ((tnow - last_downstream) > 0.5)
            from_downstream = needed_upstream && !needed_downstream && have_downstream && ((tnow - last_upstream) > 0.5)

            tasklog(env, "up,down | shared(", upstat.nshared, ",", downstat.nshared, ") created(", upstat.ncreated, ",", downstat.ncreated, ") del(", upstat.ndeleted, ",", downstat.ndeleted, ") share(", needed_upstream, ",", needed_downstream, ")")

            if from_upstream
                tasklog(env, "bring new upstream tasks into local node")
                # bring new upstream tasks into local node
                task = hascosts ? steal(upenv, broker_from_upstream_task_selector) : steal(upenv)
                tasklog(env, "stole ", (task === NoTask) ? "NoTask" : string(task))
                if task !== NoTask
                    enqueue(env, task, false, true)
                    from_downstream = false
                    last_upstream = tnow
                end
            end

            if from_downstream
                tasklog(env, "export shared node tasks to upstream brokers")
                # export shared node tasks to upstream brokers
                task = hascosts ? steal(env, broker_from_downstream_task_selector) : steal(env)
                tasklog(env, "stole ", (task === NoTask) ? "NoTask" : string(task))
                if task !== NoTask
                    enqueue(upenv, task, false, true)
                    last_downstream = tnow
                end
            end

            tasklog(env, "waiting on trigger")
            # wait for trigger
            take!(unified_trigger)
            while isready(unified_trigger)
                take!(unified_trigger)
            end
        end
    catch ex
        close(unified_trigger)
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    end
    nothing
end

#-------------------------------------------------------------------
# broker provides the initial tasks and coordinates among executors
# by stealing spare tasks from all peers and letting peers steal
#--------------------------------------------------------------------
function runbroker(rootpath::String, id::UInt64, upstream_brokerid::UInt64, root_t::Thunk, costs::Vector;
        debug::Bool=false, remotetrack::Bool=false, profile::Bool=false, upstream_help_threshold::Int=typemax(Int), downstream_help_threshold::Int=typemax(Int))
    profile_init(profile, "broker_$(id)")
    remotetrack_init(remotetrack)
    if genv[] === nothing
        env = genv[] = ExecutionCtx(META_IMPL[:node], rootpath, id, id, :broker, downstream_help_threshold; debug=debug)
    else
        env = (genv[])::ExecutionCtx
    end

    if upstream_genv[] === nothing
        upenv = upstream_genv[] = ExecutionCtx(META_IMPL[:cluster], rootpath, id, upstream_brokerid, :executor, upstream_help_threshold; debug=debug)
    else
        upenv = (upstream_genv[])::ExecutionCtx
    end

    env.debug = debug
    env.remotetrack = remotetrack
    upenv.debug = debug
    upenv.remotetrack = remotetrack
    init(env, root_t; result_callback=(task,res)->local_result_to_upstream(upenv,task,res))
    init(upenv, root_t; result_callback=(task,res)->upstream_result_to_local(env,task,res))
    env.costs = upenv.costs = costs
    remotecall_wait(join_cluster, Int(upstream_brokerid))
    tasklog(env, "invoked")

    try
        root = taskid(root_t)
        unified_trigger = Channel{Void}(10)
        do_trigger = Ref(true)

        @sync begin
            # monitor the global and node queues
            @async unify_trigger(upenv, root, unified_trigger, do_trigger)
            @async unify_trigger(env, root, unified_trigger, do_trigger)

            # broker tasks between this node and other nodes
            @async broker_tasks(upenv, env, unified_trigger, do_trigger)
        end

        tasklog(env, "stole ", env.nstolen, " shared ", env.nshared, " tasks")
        #info("broker stole $(env.nstolen), shared $(env.nshared) tasks")

        return (env.nstolen, env.nshared)
    catch ex
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    finally
        remotetrack_end(remotetrack)
        profile_end(profile, "broker_$(id)")
        reset(env)
        reset(upenv)
    end
end

#-------------------------------------------------------------------
# broker provides the initial tasks and coordinates among executors
# by stealing spare tasks from all peers and letting peers steal
#--------------------------------------------------------------------
function master_schedule(env, execstages, scheduled, completed)
    tasklog(env, "scheduling ", length(execstages), " execstages entries")
    # filter out fully scheduled stages
    isempty(execstages) || filter!((task)->!(taskid(task) in scheduled), execstages)

    if !isempty(execstages)
        schedulable_depth = Vector{Int}()
        schedulable_tids = Vector{TaskIdType}()
        for task in execstages
            # filter out child stages that are completed
            filter!(x->(isa(x, Thunk) && !(taskid(x) in completed)), task)
            # pick child stages that have no inputs or all inputs ready (indicated by the same condition)
            walk_dag(task, false) do x,d
                if isa(x, Thunk)
                    tid = taskid(x)
                    if !(tid in scheduled) && isempty(x.inputs)
                        push!(schedulable_depth, d)
                        push!(schedulable_tids, tid)
                    end
                end
                nothing
            end
        end
        tasklog(env, "found ", length(schedulable_tids), " schedulable tasks")
        # schedule for execution, in decreasing order of node depth
        sp = sortperm(schedulable_depth; rev=true)
        for idx in 1:length(sp)
            tid = schedulable_tids[sp[idx]]
            if !(tid in scheduled)
                keep(env, tid, 0, false)
                push!(scheduled, tid)
            end
        end
    end
    nothing
end

function runmaster(runenv::RunEnv, root_t::Thunk, execstages, costs::Vector; debug::Bool=false, remotetrack::Bool=false)
    rootpath = runenv.rootpath
    id = UInt64(myid())
    brokerid = runenv.masterid
    if genv[] === nothing
        env = genv[] = ExecutionCtx(META_IMPL[:cluster], rootpath, id, brokerid, :master, typemax(Int); debug=debug)
    else
        env = (genv[])::ExecutionCtx
    end
    env.debug = debug
    env.remotetrack = remotetrack
    env.costs = costs

    completed = Vector{TaskIdType}()
    scheduled = Vector{TaskIdType}()

    init(env, root_t; result_callback=(task,res)->push!(completed,task))
    tasklog(env, "invoked")

    remotecall_wait(join_cluster, Int(id))
    started = false
    @schedule while !started && !isready(start_cond)
        check_failures(runenv, env)
        sleep(0.01)
    end
    take!(start_cond)
    started = true
    tasklog(env, "all brokers and executors joined cluster")

    try
        root = taskid(root_t)
        tasklog(env, "started with ", root)

        last_ncompleted = 0
        while !has_result(env.meta, root)
            # while critical path items exist, process them first
            if isempty(scheduled) || (length(completed) > last_ncompleted)
                last_ncompleted = length(completed)
                master_schedule(env, execstages, scheduled, completed)
            end
            wait_trigger(env.meta)
            check_failures(runenv, env)
        end

        tasklog(env, "stole ", env.nstolen, " shared ", env.nshared, " tasks")
        #info("broker stole $(env.nstolen), shared $(env.nshared) tasks")
        return get_result(env.meta, root)
    catch ex
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    finally
        reset(env)
    end
end

function runexecutor(rootpath::String, id::UInt64, brokerid::UInt64, masterid::UInt64, root_t::Thunk;
        debug::Bool=false, remotetrack::Bool=false, profile::Bool=false, help_threshold::Int=typemax(Int))
    profile_init(profile, "executor_$(brokerid)_$(id)")
    remotetrack_init(remotetrack)
    if genv[] === nothing
        env = genv[] = ExecutionCtx(META_IMPL[:node], rootpath, id, brokerid, :executor, help_threshold; debug=debug)
    else
        env = (genv[])::ExecutionCtx
    end
    env.debug = debug
    env.remotetrack = remotetrack
    tasklog(env, "starting")
    init(env, root_t)
    remotecall_wait(join_cluster, Int(masterid))

    try
        root = taskid(root_t)
        task = NoTask

        while !has_result(env.meta, root)
            tasklog(env, "trying to do tasks, remaining ", length(env.reserved))

            # check if we should share some of our reserved tasks out
            if should_share(env, length(env.reserved))
                reserve_to_share(env)
            end

            task = reserve(env)

            # if no tasks in own queue, steal tasks from peers
            if task === NoTask
                tasklog(env, "got notask or lasttask, will steal")
                task = steal(env)

                if task !== NoTask
                    keep(env, task, 0, true)
                    tasklog(env, "stole ", task, ", remaining ", length(env.reserved))
                end
            end

            if task !== NoTask
                # execute task
                if !(complete = has_result(env.meta, task))
                    complete = keep(env, task, 1, true)
                    if !complete
                        if runnable(env, task)
                            complete = exec(env, task)
                        end
                    end
                end
                release(env, task, complete)
                tasklog(env, task, " is complete: ", complete, ", remaining ", length(env.reserved))
            else
                wait_trigger(env.meta)
            end
        end
        tasklog(env, "stole ", env.nstolen, ", shared ", env.nshared, ", completed ", env.nexecuted, " tasks")
        return (env.nstolen, env.nshared, env.nexecuted)
    catch ex
        taskexception(env, ex, catch_backtrace())
        rethrow(ex)
    finally
        remotetrack_end(remotetrack)
        profile_end(profile, "executor_$(brokerid)_$(id)")
        reset(env)
    end
end

function cleanup(runenv::RunEnv)
    if nothing !== runenv.reset_task
        wait(runenv.reset_task)
        runenv.reset_task = nothing
    end
    delete_meta(runenv.masterid, runenv.nodes, runenv.rootpath)
    for node in runenv.nodes
        empty!(node.executor_tasks)
    end
    @everywhere DagScheduler.genv[] = nothing
    @everywhere DagScheduler.upstream_genv[] = nothing
    cleanup_meta(runenv.masterid, runenv.nodes, runenv.rootpath)
    nothing
end

function rundag(runenv::RunEnv, dag::Thunk)
    profile_init(runenv.profile, "master_$(runenv.masterid)")
    remotetrack_init(runenv.remotetrack)
    #=
    dag, _elapsedtime, _bytes, _gctime, _memallocs = @timed begin
        dref_to_fref(dag)
    end
    =#
    #info("dag preparation time: $_elapsedtime")

    # determine critical path
    execstages, costs = schedule(runenv, dag)

    if nothing !== runenv.reset_task
        wait(runenv.reset_task)
        runenv.reset_task = nothing
    end

    upstream_help_threshold = length(runenv.nodes) - 1
    executor_help_threshold = sum(length(node.executorids) for node in runenv.nodes)
    join_count[] = 1 # start with 1 for the master
    while isready(start_cond)
        take!(start_cond)
    end

    for node in runenv.nodes
        downstream_help_threshold = length(node.executorids) - 1

        for idx in 1:length(node.executorids)
            executorid = node.executorids[idx]
            #info("spawning executor $executorid")
            join_count[] += 1
            executor_task = @spawnat executorid runexecutor(runenv.rootpath, executorid, node.brokerid, runenv.masterid, dag; debug=runenv.debug, remotetrack=runenv.remotetrack, profile=runenv.profile, help_threshold=executor_help_threshold)
            push!(node.executor_tasks, executor_task)
        end

        if node.brokerid !== runenv.masterid
            #info("spawning broker $(node.brokerid)")
            join_count[] += 1
            node.broker_task = @spawnat node.brokerid runbroker(runenv.rootpath, node.brokerid, runenv.masterid, dag, costs;
                debug=runenv.debug, remotetrack=runenv.remotetrack, profile=runenv.profile,
                upstream_help_threshold=upstream_help_threshold,
                downstream_help_threshold=downstream_help_threshold)
        end
    end

    #info("spawning master broker")
    res = runmaster(runenv, dag, execstages, costs; debug=runenv.debug, remotetrack=runenv.remotetrack)

    runenv.reset_task = @schedule wait_for_executors(runenv)
    res
end


text:
using RecurrenceAnalysis, SparseArrays
using Test

cells = [(1,1),(2,2),(3,2),(6,2),(7,2),(8,2),(6,3),(7,3),(10,3),(2,4),(4,4),
    (1,5),(2,5),(3,5),(5,5),(9,5),(10,5),(2,6),(3,6),(6,6),(7,6),(10,6),
    (3,7),(8,7),(4,8),(1,9),(5,9),(7,9),(9,9),(9,10),(10,10),(3,11),(7,11)]
i = [x[1] for x in cells]
j = [x[2] for x in cells]
rmat = CrossRecurrenceMatrix(sparse(i,j,trues(length(i))))
### Graph of the recurrence matrix:
## · : non-recurrent (white) points
## x : recurrent (black) points
## o : recurrent points within Theiler window = 2
#
#    1 2 3 4 5 6 7 8 9 A B
#  1 o · · · x · · · x · ·
#  2 · o · x x x · · · · ·
#  3 · o · · x x x · · · x
#  4 · · · o · · · x · · ·
#  5 · · · · o · · · x · ·
#  6 · x x · · o · · · · ·
#  7 · x x · · o · · x · x
#  8 · x · · · · o · · · ·
#  9 · · · · x · · · o o ·
#  A · · x · x x · · · o ·
#
###

@testset "Recurrence structures" begin
    @testset "Default parameters" begin
        histograms = recurrencestructures(rmat)
        true_histograms = Dict("diagonal" => [11,7,1,0,1],
                               "vertical" => [15,6,2],
                               "recurrencetimes" => [3,1,4,4])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=0, lmin=1, border=1)
        @test rqa_params[:RR] == 33/110
        @test rqa_params[:DET] == 1.0
        @test rqa_params[:L] == 33/20
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.996 atol=0.001
        @test rqa_params[:TREND] ≈ -33.8 atol=0.1
        @test rqa_params[:LAM] == 1.0
        @test rqa_params[:TT] == 33/23
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 33/12
        @test rqa_params[:RTE] ≈ 1.286 atol=0.001
        @test rqa_params[:NMPRT] == 4
    end
    @testset "With Theiler window" begin
        histograms = recurrencestructures(rmat, theiler=2)
        true_histograms = Dict("diagonal" => [9,4,0,0,1],
                               "vertical" => [10,3,2],
                               "recurrencetimes" => [1,0,3,0,0,0,2])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=2, lmin=1, border=1)
        @test rqa_params[:RR] == 22/110
        @test rqa_params[:DET] == 1.0
        @test rqa_params[:L] == 22/14
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.830 atol=0.001
        @test rqa_params[:TREND] ≈ -12.2 atol=0.1
        @test rqa_params[:LAM] == 1.0
        @test rqa_params[:TT] == 22/15
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 24/6
        @test rqa_params[:RTE] ≈ 1.011 atol=0.001
        @test rqa_params[:NMPRT] == 3
    end
    @testset "With minimum line" begin
        histograms = recurrencestructures(rmat, lmin=2)
        true_histograms = Dict("diagonal" => [11,7,1,0,1],
                               "vertical" => [15,6,2],
                               "recurrencetimes" => [3,1,4,4])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=0, lmin=2, border=1)
        @test rqa_params[:RR] == 33/110
        @test rqa_params[:DET] == 22/33
        @test rqa_params[:L] == 22/9
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.684 atol=0.001
        @test rqa_params[:TREND] ≈ -33.8 atol=0.1
        @test rqa_params[:LAM] == 18/33
        @test rqa_params[:TT] == 18/8
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 33/12
        @test rqa_params[:RTE] ≈ 1.286 atol=0.001
        @test rqa_params[:NMPRT] == 4
    end
    @testset "Theiler and minimum line" begin
        histograms = recurrencestructures(rmat, theiler=2, lmin=2)
        true_histograms = Dict("diagonal" => [9,4,0,0,1],
                               "vertical" => [10,3,2],
                               "recurrencetimes" => [1,0,3,0,0,0,2])
        for k in keys(histograms)
            @test histograms[k] == true_histograms[k]
        end
        rqa_params = rqa(rmat, theiler=2, lmin=2, border=1)
        @test rqa_params[:RR] == 22/110
        @test rqa_params[:DET] == 13/22
        @test rqa_params[:L] == 13/5
        @test rqa_params[:Lmax] == 5
        @test rqa_params[:DIV] == 0.2
        @test rqa_params[:ENTR] ≈ 0.500 atol=0.001
        @test rqa_params[:TREND] ≈ -12.2 atol=0.1
        @test rqa_params[:LAM] == 12/22
        @test rqa_params[:TT] == 12/5
        @test rqa_params[:Vmax] == 3
        @test rqa_params[:MRT] == 24/6
        @test rqa_params[:RTE] ≈ 1.011 atol=0.001
        @test rqa_params[:NMPRT] == 3
    end
end


text:
export compute_and_plot

# function round_values!(df::DataFrame, round::Function)
#     for i in 1:ncol(df)
#         if eltype(df[:, i]) <: Real
#             df[:, i] = round.(df[:, i])
#         end
#     end
# end

function compute_and_plot(df_cases::DataFrame, case_name::String, k_gen::Int, ylabel_R::String, ylabel_N::String, databasis::String, pop_days::Int=0)
    label_4_days = "RKI Nowcast 4 days"
    label_7_days = "RKI Nowcast 7 days"
    label_projected_7_days = "Projected 7 days"

    N_4_days, R_4_days = build_R(df_cases, k_gen - 1, 0, k_gen)
    N_7_days, R_7_days = build_R(df_cases, 5, 1, k_gen)
    N_projected_7_days, R_projected_7_days = build_R_acausal(df_cases, k_gen - 1, 3, k_gen)

    last_day = last(R_projected_7_days).days - Day(pop_days)

    Ns = [N_4_days, N_7_days, N_projected_7_days]
    Rs = [R_4_days, R_7_days, R_projected_7_days]
    labels = [label_4_days, label_7_days, label_projected_7_days]

    # rename columns
    for (label, i) in zip(labels, 1:length(Ns))
        Ns[i] = rename(Ns[i], :cases => label)
        Rs[i] = rename(Rs[i], :R => label)
    end

    # join columns based on days
    N = outerjoin(Ns...; on = :days)
    R = outerjoin(Rs...; on = :days)

    # replace missing values
    N = coalesce.(N, NaN)
    R = coalesce.(R, NaN)

    # restrict to meaningful days, i.e. most recent day - Day(pop_days)
    N = N[N.days .<= last_day, :]
    R = R[R.days .<= last_day, :]

    # # round values
    # N = round_values(N, n -> !isnan(n) ? round(Int, n) : n)
    # R = round_values(R, r -> round(r, digits=1))

    # write to CSV
    CSV.write("results-N-"*case_name*".csv", N, delim = ",")
    CSV.write("results-R-"*case_name*".csv", R, delim = ",")

    # plot
    last_R_proj_7 = round(R[end, label_projected_7_days], digits=2)
    plot_df(R, "days", ylabel_R, "R based on $(databasis). R = $last_R_proj_7. (Last updated: $(today()); valid for $last_day)", "reproduction-numbers-"*case_name*".png", :topright)
    plot_df(N, "days", ylabel_N, "Number of cases based on $(databasis). (Last updated: $(today()); valid for $last_day)", "cases-"*case_name*".png", :topleft)

    # return values
    N, R
end

function compute_and_plot(d::Dict, case_name::String, k_gen::Int, ylabel_R::String, ylabel_N::String, databasis::String, pop_days::Int=0)
    N = Vector{DataFrame}()
    R = Vector{DataFrame}()
    for (name, df) in d
        display("processing $name")
        n, r = compute_and_plot(df, case_name*"-"*name, k_gen, ylabel_R*"-"*name, ylabel_N*"-"*name, databasis, pop_days)
        push!(N, n)
        push!(R, r)
    end
    N, R
end

function cut_to_day!(df::DataFrame, day::Date)
    df = df[df.days .<= day, :]
end

function plot_df(df, day_col::String, ylabel::String, title::String, file_name::String, position=:topright)
    days = df[!, day_col]
    data = df[!, .!(names(df) .== day_col)]

    closeall()
    gr()
    for (i, name) in enumerate(names(data))
        plot_fun = i==1 ? plot : plot!
        plot_fun(days, data[!, name], label=name, ylabel=ylabel, titlefontsize=7, margin=7mm, legend=position, title=title)
    end
    savefig(file_name)
end


text:
ALL_TASK_OFFSETS[("windows-x86_64", v"1.5.3")] = Dict(
    :END => 704,
    :backtrace => 56,
    :bufsz => 376,
    :copy_stack => 384,
    :ctx => 96,
    :donenotify => 32,
    :eh => 392,
    :exception => 48,
    :excstack => 408,
    :gcstack => 400,
    :locks => 432,
    :logstate => 64,
    :next => 0,
    :prio => 426,
    :queue => 8,
    :result => 40,
    :sizeof_ctx => 272,
    :start => 72,
    :state => 24,
    :sticky => 80,
    :stkbuf => 368,
    :tid => 424,
    :timing_stack => 688,
    :tls => 16,
    :tls_base_context => 6608,
    :world_age => 416,
)


text:
using NPZ

# We expects the data to be in npy format, return a dict of {group: items}, each file is a different group
function load_data(path::String; prefix::String="", swapDimension::Bool = true)
    groups_dict = Dict()
    arr = npzread(path * prefix * ".npy")
    for (index, value) in enumerate(arr)
        if isnan(value)
            arr[index] = 0.0
        end
    end
    return swapDimension ? transpose(arr) : arr
end

# #Note that we expect the log_likelihood_array to be in rows (samples) x columns (clusters) , this is due to making it more efficent that way.
function sample_log_cat_array!(labels::AbstractArray{Int64,1}, log_likelihood_array::AbstractArray{Float32,2})
    # println("lsample log cat" * string(log_likelihood_array))
    log_likelihood_array[isnan.(log_likelihood_array)] .= -Inf #Numerical errors arent fun
    max_log_prob_arr = maximum(log_likelihood_array, dims = 2)
    log_likelihood_array .-= max_log_prob_arr
    map!(exp,log_likelihood_array,log_likelihood_array)
    # println("lsample log cat2" * string(log_likelihood_array))
    sum_prob_arr = sum(log_likelihood_array, dims =[2])
    log_likelihood_array ./=  sum_prob_arr
    for i=1:length(labels)
        labels[i] = sample(1:size(log_likelihood_array,2), ProbabilityWeights(log_likelihood_array[i,:]))
    end
end

function create_sufficient_statistics(dist::distribution_hyper_params, pts::Array{Any,1})
    return create_sufficient_statistics(dist,dist, Array{Float32}(undef, 0, 0))
end

function get_labels_histogram(labels)
    hist_dict = Dict()
    for v in labels
        if haskey(hist_dict,v) == false
            hist_dict[v] = 0
        end
        hist_dict[v] += 1
    end
    return sort(collect(hist_dict), by=x->x[1])
end

function get_node_leaders_dict()
    leader_dict = Dict()
    cur_leader = (nworkers()== 0 ? procs() : workers())[1]
    leader_dict[cur_leader] = []
    for i in (nworkers()== 0 ? procs() : workers())
        if i in procs(cur_leader)
            push!(leader_dict[cur_leader], i)
        else
            cur_leader = i
            leader_dict[cur_leader] = [i]
        end
    end
    return leader_dict
end

function log_multivariate_gamma(x::Number, D::Number)
    res::Float64 = D*(D-1)/4*log(pi)
    for d = 1:D
        res += logabsgamma(x+(1-d)/2)[1]
    end
    return res
end

function dcolwise_dot!(r::AbstractArray, a::AbstractMatrix, b::AbstractMatrix)
    n = length(r)
    @simd for j = 1:n
        v = zero(promote_type(eltype(a), eltype(b)))
        @simd for i = 1:size(a, 1)
            @fastmath @inbounds v += a[i, j]*b[i, j]
        end
        r[j] = v
    end
end


text:
@testset "BernoulliLikelihood" begin
    lik = BernoulliLikelihood()
    test_interface(lik, SqExponentialKernel(), rand(10))
end


text:
module ABMExamples

include("SchellingsSegregation.jl")
include("ForestFire.jl")
include("HKOpinionDynamics.jl")
include("Flocking.jl")
include("PredatorPrey.jl")

import .SchellingExample: run_schelling_model!
import .ForestFireExample: run_forest_fire_model!
import .HKOpinionDynamicsExample: hk_model_run_and_plot!
import .FlockingExample: run_flocking_example!
import .PredatorPreyExample: run_predatorprey_model!

export run_schelling_model!
export run_forest_fire_model!
export hk_model_run_and_plot!
export run_flocking_example!
export run_predatorprey_model!

end


text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
using LinearAlgebra
using RecursiveArrayTools
using Distributions
using Distances
using Random
using Profile
using JLD
using MAT
using IterativeSolvers
using BlockArrays
Random.seed!(123)

# Creating uniform random data
NTrain = 30000
d =  2
NTest = 50
N = NTest + NTrain
x = rand(d, N)
setSplit = randperm(N)
trainInds = setSplit[1:NTrain]
testInds = setSplit[(NTrain + 1) : end]
xTrain = x[ :, trainInds ]
xTest = x[:, testInds ]
x = hcat(xTrain, xTest)

# # Creating gridded data
# n = 20
# h = 1 / (n - 1)
# xTrain = Matrix(VectorOfArray([[i * h; j * h] for i in 1:n, j in 1:n][:]))
# NTrain = size(xTrain,2)
# NTest = 50
# xTest = rand(2, NTest)
# N = NTrain + NTest
# x = hcat(xTrain, xTest) # .+ 0.5 * (1/sqrt(N)) * randn(2, N)

ρ = 3.0
λ = 1.5

# reordering points and forming skeletons
@time colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
x = x[:, P]

LTest = SparseMatrixCSC(N, N, colptr, rowval, ones(length(rowval)))
LTest22 = LTest[(NTest + 1) : end, (NTest + 1) : end]

xTest = x[:, 1:NTest]
xTrain = x[:, (NTest + 1) : end]
skeletons = construct_skeletons(colptr, rowval, distances, λ)

l = 1.1
σ = vcat(zeros(NTest), fill(1.0, NTrain))
# σ = vcat(zeros(NTest), 0.1 * exp.(randn(NTrain)) .+ 1e-4)
σ = PseudoBlockArray(σ, [NTest, NTrain])
# covfunc = r -> exp(-r/l)
# covfunc = r -> matern12!(r,l)
function covfunc!(r)
  matern12!(r,l)
end

# Creating
Θ = covfunc!(pairwise(Euclidean(), x, dims=2))
L_Θ = cholesky(Matrix(Θ))
y = L_Θ.L * randn(N)
Θ = Θ + Diagonal(σ)

# Creating data
Θ = PseudoBlockArray(Θ, [NTest, NTrain], [NTest, NTrain])
y = PseudoBlockArray(y, [NTest, NTrain])
y[Block(2)] .+= Diagonal(sqrt.(σ[Block(2)])) * randn(NTrain)

L = assembleL(skeletons, x, covfunc!)
LNaive = assembleL(skeletons, x, covfunc!, σ)
L22 = L[(NTest + 1) : end, (NTest + 1) : end]
A22 = L22 * L22' + Diagonal(1 ./ max.(σ[Block(2)], 1e-30))
# U22 = sparse(triu(A22))
U22 = squareSparse(L22)+ Diagonal(1 ./ max.(σ[Block(2)], 1e-30))

icholU_high_level!(U22)
# U22 = UpperTriangular(U22)
U22Exct = cholesky(Matrix(A22)).U

# The exact posterior mean
μFull = Θ[Block(1,2)] * (Θ[Block(2,2)] \ y[Block(2)])
# display(μFull)

# The vecchia posterior mean
# z = (Diagonal(σ[Block(2)]) \ (  U22 \ (U22' \ (L22 * (L22' * y[Block(2)])))))
noisecov = NoiseCov(L22, U22, σ[Block(2)])
# z = (L22 * (L22' * y[Block(2)]))
# z = cg(A22 , z, Pl = IChol(U22',U22), verbose=true, maxiter=10, tol=1e-16)
# z = (Diagonal(σ[Block(2)]))\ z

z = noisecov \ (y[Block(2)])

μVechia = (L' \ (L \ (vcat(zero(y[Block(1)]), z))))[1:NTest]
# display(μVechia)
@show mean(abs.(μFull - μVechia))
@show median(abs.(μFull - μVechia))

# The vecchia posterior mean
# z = (Diagonal(σ[Block(2)]) \ ( U22Exct \ (U22Exct' \ (L22 * (L22' * y[Block(2)])))))
z = (L22 * (L22' * y[Block(2)]))
z = cg(A22 , z, Pl = IChol(U22Exct',U22Exct), verbose=true, maxiter=10)
z = (Diagonal(σ[Block(2)]))\ z

μVechiaExct = (L' \ (L \ (vcat(zero(y[Block(1)]), z))))[1:NTest]
# display(μVechiaExct)
@show mean(abs.(μFull - μVechiaExct))
@show median(abs.(μFull - μVechiaExct))

# The naive application of vecchia
μNaive = -LNaive[1:NTest, 1:NTest]' \ (LNaive[(NTest + 1) : end, 1 : NTest]' * y[Block(2)])
# display(μNaive)
@show mean(abs.(μFull - μNaive))
@show median(abs.(μFull - μNaive))

# # Direct application of Vecchia
# # μ = Θ[Block(1,2)] * (L * (L' * y[Block(2)]))
#
# # Version with Cholesky factor
# μ = Θ[Block(1,2)] * (σ \ (  U \ (U' \ (L * (L' * y[Block(2)])))))
#
# display(μTrue)
# display(μ)
# display(y[Block(1)])
#
# @show norm(μTrue - μ) / norm(μTrue)
#
#


text:
# Have all been moved to ModelConstructors.jl

#=using DSGE, ModelConstructors
using Test, Distributions, InteractiveUtils, Nullables

@testset "Ensure transformations to the real line/model space are valid" begin
    for T in subtypes(Transform)
        global u = parameter(:σ_pist, 2.5230, (1e-8, 5.), (1e-8, 5.), T(), fixed=false)
        @test ( transform_to_real_line(u) |> x -> transform_to_model_space(u,x) ) == u.value

        if !isa(T,Type{ModelConstructors.Untransformed})
            # check transform_to_real_line and transform_to_model_space to different things if T is not ModelConstructors.Untransformed
            @test transform_to_real_line(u,u.value) != transform_to_model_space(u,u.value)
        end
    end
end

# probability
N = 10^2
u = parameter(:bloop, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.SquareRoot(); fixed = true)
v = parameter(:cat, 2.5230, (1e-8, 5.), (1e-8, 5.), ModelConstructors.Exponential(), Gamma(2.00, 0.1))

pvec =  ParameterVector{Float64}(undef, N)
for i in 1:length(pvec)
	pvec[i] = (i%2 == 0) ? u : v
end
@testset "Check logpdf/pdf function approximations" begin
    @test logpdf(pvec) ≈ 50*logpdf(v)
    @test pdf(pvec) ≈ exp(50*logpdf(v))
end

updated = update(pvec, ones(length(pvec)))
ModelConstructors.update!(pvec, ones(length(pvec)))

@testset "Check if update! preserves dimensions and values" begin
    @test all(updated .== pvec)
    @test logpdf(pvec) == logpdf(updated)
end

# test we only update non-fixed parameters
@testset "Ensure only non-fixed parameters are updated" begin
    for p in pvec
        if p.fixed
            @test p.value == 2.5230
        elseif isa(p, Parameter)
            @test p.value == one(Float64)
        end
    end
end

# vector of new values must be the same length
@testset "Ensure update! enforces the same length of the parameter vector being updated" begin
    @test_throws AssertionError ModelConstructors.update!(pvec, ones(length(pvec)-1))
end

@testset "Ensure parameters being updated are of the same type." begin
    for w in [parameter(:moop, 3.0, fixed=false), parameter(:moop, 3.0; scaling = log, fixed=false)]
        # new values must be of the same type
        @test_throws MethodError parameter(w, one(Int))

        # new value is out of bounds
        @test_throws ParamBoundsError parameter(w, -1.)
    end
end

# subspecs
function sstest(m::AnSchorfheide)

    # Change all the fields of an unfixed parameter
    m <= parameter(:ι_w, 0.000, (0.0, .9999), (0.0,0.9999), ModelConstructors.Untransformed(), Normal(0.0,1.0), fixed=false,
                   description="ι_w: A new parameter.",
                   tex_label="\\iota_w")

    # Change an unfixed parameter to be fixed
    m <= parameter(:ι_p, 0.000, fixed=true,
                   description= "ι_p: Another new parameter",
                   tex_label="\\iota_p")

    # Change a fixed parameter
    m <= parameter(:δ, 0.02,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    # Overwrite a fixed parameter with an unfixed parameter
    m <= parameter(:ϵ_p, 0.750, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),        fixed=false,  scaling = x -> 1 + x/100,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    steadystate!(m)
end

m = AnSchorfheide()
sstest(m)

@testset "Test steady-state parameters" begin
    @test m[:ι_w].value == 0.0
    @test m[:ι_w].valuebounds == (0.0, .9999)
    @test m[:ι_w].transform == ModelConstructors.Untransformed()
    @test m[:ι_w].transform_parameterization == (0.0,0.9999)
    @test isa(m[:ι_w].prior.value, Normal)

    @test m[:ι_p].value == 0.0
    @test m[:ι_p].valuebounds == (0.0, 0.0)
    @test isnull(m[:ι_p].prior)
    @test m[:ι_p].fixed == true

    @test m[:δ].value == 0.02
    @test m[:δ].valuebounds == (0.02, 0.02)
    @test isnull(m[:δ].prior)
    @test m[:δ].fixed == true

    @test m[:ϵ_p].value == 0.750
    @test m[:ϵ_p].transform == ModelConstructors.Exponential()
    @test isa(m[:ϵ_p].prior.value, Gamma)
    @test m[:ϵ_p].fixed==false
end

nothing
=#