text:
"""
branchAndBound(
    model::SequencingModel,
    maxWidth::U;

    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}

Solve a sequencing problem to optimality using a branch and bound scheme
Returns bestKnownSolution (Vector{Int}), bestKnownValue (Real), bound (Real)

# Arguments
- `model::SequencingModel`: The sequencing problem
- `maxWidth::Int`: The max width allowed in a decision diagram layer
- `orderingType::OrderingOption`:Which heuristic to use to sort the edges, defaults to 'edge'
- `bestKnownSolution::Union{Vector{Int},Nothing}`: An optional parameter that gives the DD a starting solution
- `bestKnownValue::Union{T,Nothing}`: An optional parameter, The value of the best known feasible solution
- `rrbSetting::Union{RRBSetting,Nothing}`: An optional parameter that chooses an rough relaxed bounding method
- `rrbPacket`: Optional data to be passed to the rrbFunction
- `timeLimit::Union{Nothing, W}`: An optional parameter that adds a time limit
- `timeElapsed::X`: An optional parameter denoting how much time passed already
- `loggingOn::Bool`: An optional parameter that can be used to turn off logging to the console and local file
- `debugOn::Bool`: An optional parameter that can be used to turn on debug statements in other functions
- `fileName::Union{Nothing,String}`: The file to write results to
"""
function branchAndBound(
    model::SequencingModel,
    maxWidth::U;
    orderingType::OrderingOption=edge,
    bestKnownSolution::Union{Array{U},Nothing}=nothing,
    bestKnownValue::Union{T,Nothing}=nothing,
    rrbSetting::Union{RRBSetting,Nothing}=nothing,
    rrbPacket=nothing,
    timeLimit::Union{Nothing, W}=nothing,
    timeElapsed::X=0.0,
    loggingOn::Bool=true,
    debugOn::Bool=false,
    fileName::Union{Nothing,String}=nothing
)where{T<:Real,U<:Int, W<:Real, X<:Real}
    startTime = time()
    bestBound = nothing
    timedOut() = isTimedOut(timeLimit,startTime, timeElapsed)
    cleanRestrictedDD(node) = restrictedSequencingDD(model, maxWidth,rootNode = node, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket)
    cleanRestrictedDD(node, domain) = restrictedSequencingDD(model, maxWidth,rootNode = node, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket, rootDomain=domain)
    cleanRelaxedDD(dd) = relaxedSequencingDD(model, maxWidth,splitOrder, dd = dd, debugOn=debugOn,bestKnownSolution=bestKnownSolution, bestKnownValue=bestKnownValue, rrbSetting=rrbSetting, rrbPacket=rrbPacket, timeLimit=timeLimit, timeElapsed=timeElapsed, startTime=startTime, widthOne=true)
    splitOrder = findSplitOrdering(model,orderingType)

     #get initial bounds
     if isnothing(bestKnownSolution)
        bestKnownSolution, bestKnownValue, exact = restrictedSequencingDD(model,maxWidth,debugOn=debugOn)
        if exact
            if loggingOn
                finishBranchAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
    end
    if timedOut()
        return bestKnownSolution, bestKnownValue, 0
    end
    #improve initial bounds using rrb
    if !isnothing(rrbSetting)
        bestKnownSolution, bestKnownValue, exact = cleanRestrictedDD(nothing)
        #if solution is found finish
        if exact
            if loggingOn
                finishBranchAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
            end
            return bestKnownSolution, bestKnownValue, bestKnownValue
        end
        if timedOut()
            return bestKnownSolution, bestKnownValue, 0
        end
    end
    if loggingOn
        logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Initial Solution:",fileName=fileName)
        logRuntime(startTime, timeElapsed, fileName=fileName)
    end
    #generate initial relaxed DD
    relaxedDD, relaxedBound, relaxedSolution, bestKnownSolution, bestKnownValue, exact = cleanRelaxedDD(nothing)

    if exact || relaxedBound == bestKnownValue
        if loggingOn
            finishBranchAndBound(bestKnownSolution, bestKnownValue,startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
        end
        return bestKnownSolution, bestKnownValue, bestKnownValue
    elseif timedOut()
        return bestKnownSolution, bestKnownValue, 0
    else
        bestBound = relaxedBound
        if loggingOn
            logOptimalityGap(bestBound, bestKnownValue, fileName=fileName)
            logRuntime(startTime, timeElapsed, fileName=fileName)
        end
    end

    #initialize queue
    queue = Vector{Tuple{RelaxedStateSequencingNode, BitVector}}()
    for node in getLastExactLayer(relaxedDD)
        filterArcs!(node.arcsOut,model,bestKnownValue = bestKnownValue,rrbSetting=rrbSetting, rrbPacket=rrbPacket)
        if !isempty(node.arcsOut)
            addNodeToBnBQueue!(queue, node,model)
        end
    end
    relaxedDD = nothing
    if loggingOn
        logQueueLength(length(queue),fileName=fileName)
    end
    #process queue
    while !isempty(queue)
        currentValue = bestKnownValue
        currentBlock = pop!(queue)
        if debugOn
            if isempty(queue)
                println("In Queue: ",length(queue)+1, "  From: ",currentBlock[1].lengthToTerminal+currentBlock[1].lengthToRoot, " - ",currentBlock[1].lengthToTerminal+currentBlock[1].lengthToRoot, "\n")
            else
                println("In Queue: ",length(queue)+1, "  From: ",currentBlock[1].lengthToTerminal+currentBlock[1].lengthToRoot, " - ",queue[1][1].lengthToTerminal+queue[1][1].lengthToRoot, "\n")
            end
        end

        bestKnownSolution, bestKnownValue, restrictedExact = cleanRestrictedDD(currentBlock[1], currentBlock[2])
        #time limit block
        if timedOut()
            return bestKnownSolution, bestKnownValue, bestBound
        end

        if !restrictedExact
            #relax the node
            relaxedDD, localBound, emptyPath, bestKnownSolution, bestKnownValue, localExactness = cleanRelaxedDD(createWidthOneDDFromNode(model,rootNode=currentBlock[1],masterDomain=currentBlock[2]))
            #time limit block
            if timedOut()
                return bestKnownSolution, bestKnownValue, bestBound
            end
            if !localExactness
                #add the DDs to the queue
                for node in getLastExactLayer(relaxedDD)
                    if isBetterSolutionValue(model, node.lengthToRoot+node.lengthToTerminal, bestKnownValue)
                        filterArcs!(node.arcsOut,model,bestKnownValue = bestKnownValue,rrbSetting=rrbSetting, rrbPacket=rrbPacket)
                        if !isempty(node.arcsOut)
                            addNodeToBnBQueue!(queue, node,model)
                        end
                    end
                end
            end
        end

        while !isempty(queue) && isBetterSolutionValue(model, bestKnownValue, queue[1][1].lengthToRoot + queue[1][1].lengthToTerminal)
            deleteat!(queue, 1)
        end

        if !isempty(queue)
            currentBound = last(queue)[1].lengthToRoot + last(queue)[1].lengthToTerminal
            if loggingOn && isBetterSolutionValue(model, bestKnownValue, currentValue)
                logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Improved Solution:", fileName=fileName)
                logOptimalityGap(currentBound, bestKnownValue, fileName=fileName)
                logRuntime(startTime, timeElapsed, fileName=fileName)
                logQueueLength(length(queue),fileName=fileName)
            elseif loggingOn && isBetterSolutionValue(model, bestBound, currentBound)
                logOptimalityGap(currentBound, bestKnownValue, fileName=fileName)
                logRuntime(startTime, timeElapsed, fileName=fileName)
                logQueueLength(length(queue),fileName=fileName)
            end
            bestBound = currentBound
        end

        #time limit block
        if timedOut()
            return bestKnownSolution, bestKnownValue, bestBound
        end
    end

    if loggingOn
        finishBranchAndBound(bestKnownSolution, bestKnownValue, startTime=startTime, timeElapsed=timeElapsed, fileName=fileName)
    end
    return bestKnownSolution, bestKnownValue, bestKnownValue
end

"""
    addNodeToBnBQueue!(queue::Vector{Tuple{RelaxedStateSequencingNode, BitVector}}, node::RelaxedStateSequencingNode,model::SequencingModel)

Adds a tuple containing a realxed node and its domain

# Arguments
- `queue::Vector{Tuple{RelaxedStateSequencingNode, BitVector}}`: The queue
- `node::RelaxedStateSequencingNode`: The block to add to the queue
- `model::SequencingModel`: The sequencing problem

"""
function addNodeToBnBQueue!(queue::Vector{Tuple{RelaxedStateSequencingNode, BitVector}}, node::RelaxedStateSequencingNode,model::SequencingModel)
    #get domain then delete ouut arcs
    domain = zeros(length(node.state))
    for arc in node.arcsOut
        domain[arc.assignment] = true
    end
    deleteArcs!(node.arcsOut)
    #delete in arcs
    deleteArcs!(node.arcsIn)

    block = (node, domain)
    #put worst bound at back of queue
    if isempty(queue)
        push!(queue, block)
    else
        insert!(
            queue,
            searchsorted(queue, block, by = x -> x[1].lengthToTerminal+x[1].lengthToRoot,rev= (model.objective!=minimization)).start,
            block
        )
    end
end

"""
    finishBranchAndBound(bestKnownSolution::Vector{Int}, bestKnownValue::T;startTime::Union{U,Nothing}=nothing, timeElapsed::V=0,fileName::Union{Nothing,String}=nothing)where{T<:Real, U<:Real, V<:Real}

Print the optimal solution

# Arguments
- `bestKnownSolution::Vector{Int}`: The solution as an ordered list of nodes
- `bestKnownValue::T`: The cost of the solution
- `startTime::U`: the time the program started (in seconds since the epoch)
- `additionalTimeElapsed::V`: additional time used (in seconds)
- `fileName::Union{Nothing,String}`: An optional file to write the log statement to
"""
function finishBranchAndBound(bestKnownSolution::Vector{Int}, bestKnownValue::T;startTime::Union{U,Nothing}=nothing, timeElapsed::V=0,fileName::Union{Nothing,String}=nothing)where{T<:Real, U<:Real, V<:Real}
    logSequencingSolution(bestKnownSolution, bestKnownValue,firstLine="Optimal Solution Found:",fileName=fileName)
    if !isnothing(startTime)
        logRuntime(startTime, timeElapsed,fileName=fileName)
    end
    println()
end


text:
@testset "nobs" begin
    @test_throws MethodError nobs(X,X)
    @test_throws MethodError nobs(X,y)

    @testset "Array, SparseArray, and Tuple" begin
        @test_throws DimensionMismatch nobs((X,XX,rand(100)))
        @test_throws DimensionMismatch nobs((X,X'))
        @test_throws DimensionMismatch nobs((X,XX), obsdim = :first)
        for var in (Xs, ys, vars...)
            @test @inferred(nobs(var, obsdim = :last)) === 150
            @test @inferred(nobs(var, obsdim = 100)) === 1
            @test @inferred(nobs(var)) === 150
            @test @inferred(nobs(var, ObsDim.Last())) === 150
        end
        @test @inferred(nobs(())) === 0
        @test @inferred(nobs((), obsdim = :first)) === 0
        @test @inferred(nobs((), obsdim = :last)) === 0
        @test @inferred(nobs((), obsdim = 3)) === 0
    end

    @testset "0-dim SubArray" begin
        v = view([3], 1)
        @test @inferred(nobs(v)) === 1
        @test @inferred(getobs(v)) === 3
        @test @inferred(getobs(v, 1)) === 3
        @test_throws BoundsError getobs(v, 2)
        @test_throws ArgumentError getobs(v, 2:3)
    end

    @testset "SubArray" begin
        @test @inferred(nobs(view(X,:,:))) === 150
        @test @inferred(nobs(view(X,:,:))) === 150
        @test @inferred(nobs(view(XX,:,:,:))) === 150
        @test @inferred(nobs(view(XXX,:,:,:,:))) === 150
        @test @inferred(nobs(view(y,:))) === 150
        @test @inferred(nobs(view(Y,:,:))) === 150
        @test @inferred(nobs(view(X,:,:), obsdim = :last)) === 150
        @test @inferred(nobs(view(XX,:,:,:), obsdim = :last)) === 150
        @test @inferred(nobs(view(XXX,:,:,:,:), obsdim = :last)) === 150
        @test @inferred(nobs(view(y,:), obsdim = :last)) === 150
        @test @inferred(nobs(view(Y,:,:), obsdim = :last)) === 150
    end

    @testset "various obsdim" begin
        @test_throws ArgumentError nobs(X, obsdim = 1.0)
        @test_throws ArgumentError nobs(X, obsdim = :one)
        @test_throws MethodError nobs(X, obsdim = ObsDim.Undefined())
        @test_throws DimensionMismatch nobs((X',X), (ObsDim.First(),ObsDim.Last(),ObsDim.Last()))
        @test_throws DimensionMismatch nobs((X',X), (ObsDim.First(),))
        @test_throws DimensionMismatch nobs((X',X), obsdim=(1,2,2))
        @test @inferred(nobs(X, ObsDim.Undefined())) === 150 # fallback
        @test @inferred(nobs(Xs, obsdim = 1)) === 10
        @test @inferred(nobs(Xs, obsdim = :first)) === 10
        @test @inferred(nobs(XXX, obsdim = 1)) === 3
        @test @inferred(nobs(XXX, ObsDim.First())) === 3
        @test @inferred(nobs(XXX, obsdim = :first)) === 3
        @test @inferred(nobs(XXX, obsdim = 2)) === 20
        @test @inferred(nobs(XXX, obsdim = 3)) === 30
        @test @inferred(nobs(XXX, ObsDim.Constant(3))) === 30
        @test @inferred(nobs(XXX, obsdim = 4)) === 150
        @test @inferred(nobs((X,y), obsdim = :last)) === 150
        @test @inferred(nobs((X',y), obsdim = :first)) === 150
        @test @inferred(nobs((X',X'), obsdim = :first)) === 150
        @test @inferred(nobs((X',X), obsdim = (:first,:last))) === 150
        @test @inferred(nobs((X',X), obsdim = (1,2))) === 150
        @test @inferred(nobs((X',X,X), obsdim = (1,2,2))) === 150
        @test @inferred(nobs((X',X,X), obsdim = (1,2,:last))) === 150
        @test @inferred(nobs((X',X), (ObsDim.First(),ObsDim.Last()))) === 150
        @test @inferred(nobs((X,X), obsdim = :first)) === 4
    end

    @testset "custom types" begin
        # test that fallback bouncing doesn't cause stackoverflow
        @test_throws MethodError nobs(EmptyType())
        @test_throws MethodError nobs(EmptyType(), ObsDim.Undefined())
        @test_throws MethodError nobs(EmptyType(), ObsDim.Last())
        @test_throws MethodError nobs(EmptyType(), (ObsDim.Last(),ObsDim.Last()))
        @test_throws MethodError nobs(EmptyType(), obsdim = 1)
        @test_throws MethodError nobs(EmptyType(), obsdim = (1,1))
        @test_throws MethodError nobs(EmptyType(), obsdim = :last)
        # test types that don't use the obsdim
        @test_throws MethodError nobs(CustomType(), obsdim = 1)
        @test_throws MethodError nobs(CustomType(), ObsDim.Last())
        @test nobs(CustomType()) === 100
    end
end

@testset "getobs" begin
    @testset "Array and Subarray" begin
        # access outside nobs bounds
        @test_throws BoundsError getobs(X, -1)
        @test_throws BoundsError getobs(X, 0)
        @test_throws BoundsError getobs(X, 0, obsdim = 1)
        @test_throws BoundsError getobs(X, 151)
        @test_throws BoundsError getobs(X, 151, obsdim = 2)
        @test_throws BoundsError getobs(X, 151, obsdim = 1)
        @test_throws BoundsError getobs(X, 5, obsdim = 1)
        @test typeof(@inferred(getobs(Xv))) <: Array
        @test typeof(@inferred(getobs(yv))) <: Array
        @test all(getobs(Xv) .== X)
        @test all(getobs(yv) .== y)
        @test getobs(X, ObsDim.Undefined()) === X
        @test getobs(X, ObsDim.Constant(1)) === X
        @test_throws MethodError getobs(X, obsdim = ObsDim.Undefined())
        @test getobs(X, obsdim = ObsDim.Constant(1)) === X
        @test @inferred(getobs(X))   === X
        @test @inferred(getobs(XX))  === XX
        @test @inferred(getobs(XXX)) === XXX
        @test @inferred(getobs(y))   === y
        @test @inferred(getobs(X, 45)) == getobs(X', 45, obsdim = 1)
        @test @inferred(getobs(X, 3:10)) == getobs(X', 3:10, obsdim = 1)'
        for i in (2, 2:20, [2,1,4])
            @test_throws ErrorException @inferred(getobs(XX, i, obsdim = 1))
            @test @inferred(getobs(XX, i, obsdim = :last) == getindex(XX,:,:,i))
            @test @inferred(getobs(XX, i, ObsDim.First())) == getindex(XX,i,:,:)
            @test @inferred(getobs(XX, i, ObsDim.Constant(2))) == getindex(XX,:,i,:)
            @test @inferred(getobs(XX, i, ObsDim.Last())) == getindex(XX,:,:,i)
            @test getobs(XX, i, obsdim = 1) == XX[i,:,:]
            @test getobs(XX, i, obsdim = :first) == XX[i,:,:]
            @test getobs(XX, i, obsdim = 2) == XX[:,i,:]
            @test getobs(XX, i, obsdim = :last) == XX[:,:,i]
        end
        for i in (2, 1:150, 2:10, [2,5,7], [2,1])
            @test_throws MethodError getobs(X, i, ObsDim.Undefined())
            @test_throws DimensionMismatch getobs(X, i, obsdim = 12)
            @test typeof(getobs(Xv, i)) <: Array
            @test typeof(getobs(yv, i)) <: ((typeof(i) <: Int) ? String : Array)
            @test all(getobs(Xv,i) .== X[:,i])
            @test @inferred(getobs(Xv,i))  == X[:,i]
            @test @inferred(getobs(X,i))   == X[:,i]
            @test @inferred(getobs(XX,i))  == XX[:,:,i]
            @test @inferred(getobs(XXX,i)) == XXX[:,:,:,i]
            @test @inferred(getobs(y,i))   == ((typeof(i) <: Int) ? y[i] : y[i])
            @test @inferred(getobs(yv,i))  == ((typeof(i) <: Int) ? y[i] : y[i])
            @test @inferred(getobs(Y,i))   == Y[:,i]
        end
    end

    @testset "SparseArray" begin
        @test @inferred(getobs(Xs)) === Xs
        @test @inferred(getobs(ys)) === ys
        @test @inferred(getobs(Xs, 45)) == getobs(Xs', 45, obsdim = 1)
        @test @inferred(getobs(Xs, 3:9)) == getobs(Xs', 3:9, obsdim = 1)'
        @test typeof(getobs(Xs,2)) <: SparseVector
        @test typeof(getobs(Xs,1:5)) <: SparseMatrixCSC
        @test typeof(getobs(ys,2)) <: Float64
        @test typeof(getobs(ys,1:5)) <: SparseVector
        for i in (2, 2:10, [2,1,4])
            @test @inferred(getobs(Xs, i, ObsDim.First())) == Xs[i,:]
            @test_throws ErrorException @inferred(getobs(Xs, i, obsdim = 1))
            @test getobs(Xs, i, obsdim = 1) == Xs[i,:]
            @test getobs(Xs, i, obsdim = :first) == Xs[i,:]
            @test getobs(Xs, i, obsdim = 2) == Xs[:,i]
            @test getobs(Xs, i, obsdim = :last) == Xs[:,i]
        end
        for i in (2, 1:150, 2:10, [2,5,7], [2,1])
            @test_throws MethodError getobs(Xs, i, ObsDim.Undefined())
            @test_throws DimensionMismatch getobs(Xs, i, obsdim = 12)
            @test @inferred(getobs(Xs,i)) == Xs[:,i]
            @test @inferred(getobs(ys,i)) == ys[i]
            @test @inferred(getobs(ys, i, obsdim = :last)) == ys[i]
            @test getobs(ys, i, obsdim = :last) == ys[i]
            @test getobs(ys, i, obsdim = :first) == ys[i]
        end
    end

    @testset "Tuple" begin
        @test_throws DimensionMismatch getobs((X,yv), obsdim=2)
        # bounds checking correctly
        @test_throws BoundsError getobs((X,y), 151)
        # special case empty tuple
        @test @inferred(getobs((), 10, obsdim = 1)) === ()
        @test getobs((), obsdim=2) === ()
        @test @inferred(getobs(())) === ()
        @test @inferred(getobs((), ObsDim.Last())) === ()
        @test @inferred(getobs((), 10)) === ()
        @test getobs((), 10, obsdim = 1) === ()
        @test @inferred(getobs((X,y))) === (X,y)
        @test @inferred(getobs((X,yv))) == (X,y)
        @test @inferred(getobs((Xv,y))) == (X,y)
        @test @inferred(getobs((Xv,yv))) == (X,y)
        @test_throws DimensionMismatch @inferred(getobs((X',y))) == (X',y)
        @test @inferred(getobs((XX,X,y))) === (XX,X,y)
        @test @inferred(getobs((XXX,XX,X,y))) === (XXX,XX,X,y)
        tx, ty = getobs((Xv,yv))
        @test typeof(tx) <: Array
        @test typeof(ty) <: Array
        tx, ty = getobs((Xv,yv), 10:50)
        @test typeof(tx) <: Array
        @test typeof(ty) <: Array
        for i in (1:150, 2:10, [2,5,7], [2,1])
            @test_throws DimensionMismatch getobs((X',y), i)
            @test_throws DimensionMismatch getobs((X,y),  i, obsdim=2)
            @test_throws DimensionMismatch getobs((X',y), i, obsdim=2)
            @test_throws DimensionMismatch getobs((X,y), i, obsdim=(1,2))
            @test_throws DimensionMismatch getobs((X,y), i, obsdim=(2,1,1))
            @test_throws DimensionMismatch getobs((XX,X,y), i, obsdim=(2,2,1))
            @test_throws DimensionMismatch getobs((XX,X,y), i, obsdim=(3,2))
            @test @inferred(getobs((X,y), i))  == (X[:,i], y[i])
            @test @inferred(getobs((X,yv), i)) == (X[:,i], y[i])
            @test @inferred(getobs((Xv,y), i)) == (X[:,i], y[i])
            @test @inferred(getobs((X,Y), i))  == (X[:,i], Y[:,i])
            @test @inferred(getobs((X,yt), i)) == (X[:,i], yt[:,i])
            @test @inferred(getobs((XX,X,y), i)) == (XX[:,:,i], X[:,i], y[i])
            @test_throws ErrorException @inferred(getobs((XX,X,y), i, obsdim=(3,2,1)))
            @test getobs((XX,X,y), i, obsdim=(3,2,1)) == (XX[:,:,i], X[:,i], y[i])
            @test getobs((X, y), i, obsdim=:last)  == (X[:,i], y[i])
            @test getobs((X',y), i, obsdim=:first) == (X'[i,:], y[i])
            @test getobs((X,yv), i, obsdim=:last)  == (X[:,i], y[i])
            @test getobs((Xv,y), i, obsdim=:last)  == (X[:,i], y[i])
            @test getobs((X, Y), i, obsdim=:last)  == (X[:,i], Y[:,i])
            @test getobs((X',y), i, obsdim=:first)  == (X'[i,:], y[i])
            @test getobs((X, y), i, obsdim=(:last,:last))  == (X[:,i], y[i])
            @test getobs((X, y), i, obsdim=(2,1))  == (X[:,i], y[i])
            @test getobs((X',y), i, obsdim=(1,1))  == (X'[i,:], y[i])
            @test getobs((X',yt), i, obsdim=(1,2))  == (X'[i,:], yt[:,i])
            @test getobs((X',yt), i, obsdim=(:first,:last))  == (X'[i,:], yt[:,i])
            @test getobs((XX,X,y), i, obsdim=:last) == (XX[:,:,i], X[:,i], y[i])
            # compare if obs match in tuple
            x1, y1 = getobs((X1,Y1), i)
            @test all(x1' .== y1)
            x1, y1, z1 = getobs((X1,Y1,sparse(X1)), i)
            @test all(x1' .== y1)
            @test all(x1 .== z1)
        end
        @test @inferred(getobs((X,y), 2)) == (X[:,2], y[2])
        @test @inferred(getobs((X,y), 2, obsdim=:last)) == (X[:,2], y[2])
        @test_throws ErrorException @inferred(getobs((X,y), 2, obsdim=(:last,:first)))
        @test getobs((X,y), 2, obsdim=:last) == (X[:,2], y[2])
        @test getobs((X,y), 2, obsdim=(:last,:first)) == (X[:,2], y[2])
        @test getobs((X,y), 2, obsdim=(2,:first)) == (X[:,2], y[2])
        @test @inferred(getobs((Xv,y), 2)) == (X[:,2], y[2])
        @test @inferred(getobs((X,yv), 2)) == (X[:,2], y[2])
        @test @inferred(getobs((X,Y), 2)) == (X[:,2], Y[:,2])
        @test @inferred(getobs((XX,X,y), 2)) == (XX[:,:,2], X[:,2], y[2])
    end

    @testset "type without getobs support" begin
        @test_throws MethodError getobs(EmptyType())
        @test_throws MethodError getobs(EmptyType(), obsdim=1)
        # test that fallback bouncing doesn't cause stackoverflow
        @test_throws MethodError getobs(EmptyType(), 1, ObsDim.Last())
        @test_throws MethodError getobs(EmptyType(), 1)
        @test_throws MethodError getobs(EmptyType(), 1:10)
        @test_throws MethodError getobs(EmptyType(), 1, obsdim=:first)
        @test_throws MethodError getobs(EmptyType(), 1, obsdim=10)
        @test_throws MethodError getobs(EmptyType(), 1:10, obsdim=:last)
    end

    @testset "custom type with getobs support" begin
        @test_throws MethodError getobs(CustomType(), obsdim=1)
        @test_throws MethodError getobs(CustomType(), ObsDim.Last())
        @test_throws MethodError getobs(CustomType(), 4:40, obsdim=1)
        @test @inferred(getobs(CustomType(), 11)) === 11
        @test @inferred(getobs(CustomType(), 4:40)) == collect(4:40)
        # defaults to 1:nobs
        @test @inferred(getobs(CustomType())) == collect(1:100)
        # No bounds checking here
        @test @inferred(getobs(CustomType(), 200)) === 200
        @test @inferred(getobs(CustomType(), [2,200,1])) == [2,200,1]
    end
end

@testset "getobs!" begin
    @testset "Array and Subarray" begin
        Xbuf = similar(X)
        # interpreted as idx
        @test_throws Exception getobs!(Xbuf, X, ObsDim.Undefined())
        @test_throws Exception getobs!(Xbuf, X, ObsDim.Constant(1))
        # obsdim not defined without some idx
        @test_throws MethodError getobs!(Xbuf, X, obsdim = ObsDim.Undefined())
        @test_throws MethodError getobs!(Xbuf, X, obsdim = ObsDim.Constant(1))
        # access outside nobs bounds
        @test_throws BoundsError getobs!(Xbuf, X, -1)
        @test_throws BoundsError getobs!(Xbuf, X, 0)
        @test_throws BoundsError getobs!(Xbuf, X, 0, obsdim = 1)
        @test_throws BoundsError getobs!(Xbuf, X, 151)
        @test_throws BoundsError getobs!(Xbuf, X, 151, obsdim = 2)
        @test_throws BoundsError getobs!(Xbuf, X, 151, obsdim = 1)
        @test_throws BoundsError getobs!(Xbuf, X, 5, obsdim = 1)
        @test @inferred(getobs!(Xbuf, X)) === Xbuf
        @test Xbuf == X
        @test all(getobs!(similar(Xv), Xv) .== X)
        @test all(getobs!(similar(yv), yv) .== y)
        @test @inferred(getobs!(similar(XX), XX))   == XX
        @test @inferred(getobs!(similar(XXX), XXX)) == XXX
        @test @inferred(getobs!(similar(y), y))     == y
        xbuf1 = zeros(4)
        xbuf2 = zeros(4)
        @test @inferred(getobs!(xbuf1, X, 45)) == getobs!(xbuf2, X', 45, obsdim = 1)
        Xbuf1 = zeros(4,8)
        Xbuf2 = zeros(8,4)
        @test @inferred(getobs!(Xbuf1, X, 3:10)) == getobs!(Xbuf2, X', 3:10, obsdim = 1)'
        # obsdim = 2
        Xbuf1 = zeros(20,150)
        @test_throws ErrorException @inferred(getobs!(Xbuf1, XX, 2, obsdim = 2))
        @test @inferred(getobs!(Xbuf1, XX, 5, ObsDim.Constant(2))) == XX[:,5,:]
        @test getobs!(Xbuf1, XX, 11, obsdim = 2) == XX[:,11,:]
        Xbuf2 = zeros(20,5,150)
        @test_throws ErrorException @inferred(getobs!(Xbuf2, XX, 6:10, obsdim = 2))
        @test @inferred(getobs!(Xbuf2, XX, 6:10, ObsDim.Constant(2))) == XX[:,6:10,:]
        @test getobs!(Xbuf2, XX, 11:15, obsdim = 2) == XX[:,11:15,:]
        # string vector
        @test getobs!("setosa", y, 1) == "setosa"
        @test getobs!(nothing, y, 1) == "setosa"
    end

    @testset "SparseArray" begin
        # Sparse Arrays opt-out of buffer usage
        @test @inferred(getobs!(nothing, Xs)) === getobs(Xs)
        @test @inferred(getobs!(nothing, Xs, 1)) == getobs(Xs, 1)
        @test @inferred(getobs!(nothing, Xs, 5:10)) == getobs(Xs, 5:10)
        @test @inferred(getobs!(nothing, Xs, 2, ObsDim.First())) == getobs(Xs, 2, obsdim=1)
        @test getobs!(nothing, Xs, 2, obsdim = 1) == getobs(Xs, 2, obsdim=1)
        @test @inferred(getobs!(nothing, ys)) === getobs(ys)
        @test @inferred(getobs!(nothing, ys, 1)) === getobs(ys, 1)
        @test @inferred(getobs!(nothing, ys, 5:10)) == getobs(ys, 5:10)
        @test @inferred(getobs!(nothing, ys, 5:10, ObsDim.First())) == getobs(ys, 5:10)
        @test getobs!(nothing, ys, 5:10, obsdim=1) == getobs(ys, 5:10)
    end

    @testset "Tuple" begin
        @test_throws MethodError getobs!((nothing,nothing), (X,y))
        @test_throws MethodError getobs!((nothing,nothing), (X,y), 1:5)
        @test_throws DimensionMismatch getobs!((nothing,nothing,nothing), (X,y))
        xbuf = zeros(4,2)
        ybuf = ["foo", "bar"]
        @test_throws DimensionMismatch getobs!((xbuf,), (X,y))
        @test_throws DimensionMismatch getobs!((xbuf,ybuf,ybuf), (X,y))
        @test_throws DimensionMismatch getobs!((xbuf,), (X,y), 1:5)
        @test_throws DimensionMismatch getobs!((xbuf,ybuf,ybuf), (X,y), 1:5)
        @test @inferred(getobs!((xbuf,ybuf),(X,y), 2:3)) === (xbuf,ybuf)
        @test xbuf == getobs(X, 2:3)
        @test ybuf == getobs(y, 2:3)
        @test @inferred(getobs!((xbuf,ybuf),(X,y), [50,150])) === (xbuf,ybuf)
        @test xbuf == getobs(X, [50,150])
        @test ybuf == getobs(y, [50,150])

        xbuf2 = zeros(2,4)
        @test @inferred(getobs!((xbuf2,ybuf),(X',y), 4:5, ObsDim.First())) === (xbuf2,ybuf)
        @test xbuf2 == getobs(X', 4:5, obsdim=1)
        @test ybuf  == getobs(y, 2:3)

        @test @inferred(getobs!((xbuf2,ybuf,xbuf),(X',y,X), 99:100, (ObsDim.First(),ObsDim.Last(),ObsDim.Last()))) === (xbuf2,ybuf,xbuf)
        @test xbuf2 == getobs(X', 99:100, obsdim=1)
        @test ybuf  == getobs(y, 99:100)
        @test xbuf == getobs(X, 99:100)

        @test getobs!((xbuf2,ybuf,xbuf),(X',y,X), 9:10, obsdim=(1,1,2)) === (xbuf2,ybuf,xbuf)
        @test xbuf2 == getobs(X', 9:10, obsdim=1)
        @test ybuf  == getobs(y, 9:10)
        @test xbuf == getobs(X, 9:10)

        @test getobs!((nothing,xbuf),(Xs,X), 3:4) == (getobs(Xs,3:4),xbuf)
        @test xbuf == getobs(X,3:4)
    end

    @testset "type without getobs support" begin
        # buffer is ignored if getobs! is not defined
        @test_throws MethodError getobs!(nothing, EmptyType())
        @test_throws MethodError getobs!(nothing, EmptyType(), 1)
        @test_throws MethodError getobs!(nothing, EmptyType(), obsdim=1)
        @test_throws MethodError getobs!(nothing, EmptyType(), ObsDim.Last())
        @test_throws MethodError getobs!(nothing, CustomType(), obsdim=1)
        @test_throws MethodError getobs!(nothing, CustomType(), ObsDim.Last())
        @test_throws MethodError getobs!(nothing, CustomType(), 4:40, obsdim=1)
    end

    @testset "custom type with getobs support" begin
        # No-op unless defined
        @test @inferred(getobs!(nothing, CustomType())) == collect(1:100)
        @test @inferred(getobs!(nothing, CustomType(), 11)) === 11
        @test @inferred(getobs!(nothing, CustomType(), 4:40)) == collect(4:40)
        # No bounds checking here
        @test @inferred(getobs!(nothing, CustomType(), 200)) === 200
        @test @inferred(getobs!(nothing, CustomType(), [2,200,1])) == [2,200,1]
    end
end

ext: jl


text:
"""
    point_type_type(::BifurcationProblem) :: Type
"""
function point_type_type end  # TODO: rename

"""
    eigvals_prototype(::BifurcationProblem, ::ContinuationCache)
"""
function eigvals_prototype end

"""
    regular_point(::Type{P <: AbstractSpecialPoint}) :: P
"""
function regular_point end

"""
    analyze!(cache::BifurcationCache, opts)
"""
function analyze! end

"""
    re_analyze!(solver::BifurcationSolver, u::AbstractVector)
"""
function re_analyze! end

#=
"""
    curves(ctx [, by]) :: Vector{<: Curve}

Return a list of curves in `ctx` (sweep, solution, or solver)
optionally grouped (chunked) by a key `by` (e.g., `:stability`).
"""
function curves end
=#

"""
    measure(curve, key)
"""
function measure end
# Better names?  I thought maybe `values` but it is used for
# associative collections.  How about `observe`?  `gauge`?
# `evaluate`?  `extract`?  Maybe just define `Base.getindex`?

#=
"""
    measure(curve, keys::Tuple{T₁, ..., Tₙ}) :: Tuple{V₁, ..., Vₙ}
"""
measure(c::Any, keys::Tuple) = measure.((c,), keys)

# This makes plotting functions easy but I'm not sure if this is fine,
# especially if I were to expose it as an API.
measure(::Any, ::Nothing) = nothing
=#


text:
# Automatically generated by gen_gm.jl, do not edit!

grav_param(::Float64, ::SolarSystemBarycenter) = 1.3271244004193938e11
grav_param(::SolarSystemBarycenter) = grav_param(Float64, SolarSystemBarycenter())
grav_param(::Type{Float64}, ::MercuryBarycenter) = 22031.78000000002
grav_param(::MercuryBarycenter) = grav_param(Float64, MercuryBarycenter())
grav_param(::Type{Float64}, ::VenusBarycenter) = 324858.59200000006
grav_param(::VenusBarycenter) = grav_param(Float64, VenusBarycenter())
grav_param(::Type{Float64}, ::EarthBarycenter) = 403503.2355022598
grav_param(::EarthBarycenter) = grav_param(Float64, EarthBarycenter())
grav_param(::Type{Float64}, ::MarsBarycenter) = 42828.37521400002
grav_param(::MarsBarycenter) = grav_param(Float64, MarsBarycenter())
grav_param(::Type{Float64}, ::JupiterBarycenter) = 1.267127648000002e8
grav_param(::JupiterBarycenter) = grav_param(Float64, JupiterBarycenter())
grav_param(::Type{Float64}, ::SaturnBarycenter) = 3.79405852e7
grav_param(::SaturnBarycenter) = grav_param(Float64, SaturnBarycenter())
grav_param(::Type{Float64}, ::UranusBarycenter) = 5.794548600000008e6
grav_param(::UranusBarycenter) = grav_param(Float64, UranusBarycenter())
grav_param(::Type{Float64}, ::NeptuneBarycenter) = 6.836527100580024e6
grav_param(::NeptuneBarycenter) = grav_param(Float64, NeptuneBarycenter())
grav_param(::Type{Float64}, ::PlutoBarycenter) = 977.0000000000007
grav_param(::PlutoBarycenter) = grav_param(Float64, PlutoBarycenter())
grav_param(::Type{Float64}, ::Sun) = 1.3271244004193938e11
grav_param(::Sun) = grav_param(Float64, Sun())
grav_param(::Type{Float64}, ::Mercury) = 22031.78000000002
grav_param(::Mercury) = grav_param(Float64, Mercury())
grav_param(::Type{Float64}, ::Venus) = 324858.59200000006
grav_param(::Venus) = grav_param(Float64, Venus())
grav_param(::Type{Float64}, ::Earth) = 398600.435436096
grav_param(::Earth) = grav_param(Float64, Earth())
grav_param(::Type{Float64}, ::Mars) = 42828.37362069909
grav_param(::Mars) = grav_param(Float64, Mars())
grav_param(::Type{Float64}, ::Jupiter) = 1.266865349218008e8
grav_param(::Jupiter) = grav_param(Float64, Jupiter())
grav_param(::Type{Float64}, ::Saturn) = 3.793120749865224e7
grav_param(::Saturn) = grav_param(Float64, Saturn())
grav_param(::Type{Float64}, ::Uranus) = 5.793951322279009e6
grav_param(::Uranus) = grav_param(Float64, Uranus())
grav_param(::Type{Float64}, ::Neptune) = 6.835099502439672e6
grav_param(::Neptune) = grav_param(Float64, Neptune())
grav_param(::Type{Float64}, ::Pluto) = 869.6138177608748
grav_param(::Pluto) = grav_param(Float64, Pluto())
grav_param(::Type{Float64}, ::Luna) = 4902.800066163796
grav_param(::Luna) = grav_param(Float64, Luna())
grav_param(::Type{Float64}, ::Phobos) = 0.0007087546066894452
grav_param(::Phobos) = grav_param(Float64, Phobos())
grav_param(::Type{Float64}, ::Deimos) = 9.615569648120313e-5
grav_param(::Deimos) = grav_param(Float64, Deimos())
grav_param(::Type{Float64}, ::Io) = 5959.916033410404
grav_param(::Io) = grav_param(Float64, Io())
grav_param(::Type{Float64}, ::Europa) = 3202.738774922892
grav_param(::Europa) = grav_param(Float64, Europa())
grav_param(::Type{Float64}, ::Ganymede) = 9887.834453334144
grav_param(::Ganymede) = grav_param(Float64, Ganymede())
grav_param(::Type{Float64}, ::Callisto) = 7179.28936139727
grav_param(::Callisto) = grav_param(Float64, Callisto())
grav_param(::Type{Float64}, ::Amalthea) = 0.1378480571202615
grav_param(::Amalthea) = grav_param(Float64, Amalthea())
grav_param(::Type{Float64}, ::Mimas) = 2.503522884661795
grav_param(::Mimas) = grav_param(Float64, Mimas())
grav_param(::Type{Float64}, ::Enceladus) = 7.211292085479989
grav_param(::Enceladus) = grav_param(Float64, Enceladus())
grav_param(::Type{Float64}, ::Tethys) = 41.21117207701302
grav_param(::Tethys) = grav_param(Float64, Tethys())
grav_param(::Type{Float64}, ::Dione) = 73.11635322923193
grav_param(::Dione) = grav_param(Float64, Dione())
grav_param(::Type{Float64}, ::Rhea) = 153.9422045545342
grav_param(::Rhea) = grav_param(Float64, Rhea())
grav_param(::Type{Float64}, ::Titan) = 8978.138845307376
grav_param(::Titan) = grav_param(Float64, Titan())
grav_param(::Type{Float64}, ::Hyperion) = 0.3718791714191668
grav_param(::Hyperion) = grav_param(Float64, Hyperion())
grav_param(::Type{Float64}, ::Iapetus) = 120.5134781724041
grav_param(::Iapetus) = grav_param(Float64, Iapetus())
grav_param(::Type{Float64}, ::Phoebe) = 0.5531110414633374
grav_param(::Phoebe) = grav_param(Float64, Phoebe())
grav_param(::Type{Float64}, ::Janus) = 0.1266231296945636
grav_param(::Janus) = grav_param(Float64, Janus())
grav_param(::Type{Float64}, ::Epimetheus) = 0.03513977490568457
grav_param(::Epimetheus) = grav_param(Float64, Epimetheus())
grav_param(::Type{Float64}, ::Atlas) = 0.0003759718886965353
grav_param(::Atlas) = grav_param(Float64, Atlas())
grav_param(::Type{Float64}, ::Prometheus) = 0.01066368426666134
grav_param(::Prometheus) = grav_param(Float64, Prometheus())
grav_param(::Type{Float64}, ::Pandora) = 0.0091037683110543
grav_param(::Pandora) = grav_param(Float64, Pandora())
grav_param(::Type{Float64}, ::Ariel) = 83.46344431770477
grav_param(::Ariel) = grav_param(Float64, Ariel())
grav_param(::Type{Float64}, ::Umbriel) = 85.09338094489388
grav_param(::Umbriel) = grav_param(Float64, Umbriel())
grav_param(::Type{Float64}, ::Titania) = 226.9437003741248
grav_param(::Titania) = grav_param(Float64, Titania())
grav_param(::Type{Float64}, ::Oberon) = 205.3234302535623
grav_param(::Oberon) = grav_param(Float64, Oberon())
grav_param(::Type{Float64}, ::Miranda) = 4.3195168992321
grav_param(::Miranda) = grav_param(Float64, Miranda())
grav_param(::Type{Float64}, ::Triton) = 1427.598140725034
grav_param(::Triton) = grav_param(Float64, Triton())
grav_param(::Type{Float64}, ::Charon) = 105.8799888601881
grav_param(::Charon) = grav_param(Float64, Charon())
grav_param(::Type{Float64}, ::Nix) = 0.00304817564816976
grav_param(::Nix) = grav_param(Float64, Nix())
grav_param(::Type{Float64}, ::Hydra) = 0.003211039206155255
grav_param(::Hydra) = grav_param(Float64, Hydra())
grav_param(::Type{Float64}, ::Kerberos) = 0.001110040850536676
grav_param(::Kerberos) = grav_param(Float64, Kerberos())
grav_param(::Type{Float64}, ::Ceres) = 63.13
grav_param(::Ceres) = grav_param(Float64, Ceres())
grav_param(::Type{Float64}, ::Pallas) = 13.73
grav_param(::Pallas) = grav_param(Float64, Pallas())
grav_param(::Type{Float64}, ::Vesta) = 17.29
grav_param(::Vesta) = grav_param(Float64, Vesta())
grav_param(::Type{Float64}, ::Psyche) = 1.81
grav_param(::Psyche) = grav_param(Float64, Psyche())
grav_param(::Type{Float64}, ::Eros) = 0.0004463
grav_param(::Eros) = grav_param(Float64, Eros())
grav_param(::Type{Float64}, ::Davida) = 2.26
grav_param(::Davida) = grav_param(Float64, Davida())


text:
module VB

using Distributions
using Flux, Flux.Tracker

import Random # shuffle, seed
import Dates
ShowTime() = Dates.format(Dates.now(), "yyyy-mm-dd HH:MM:SS")

import Cytof5.Model: solveBeta, gen_beta_est

include("ADVI/ADVI.jl")
include("Priors.jl")
include("Constants.jl")
include("vae.jl")
include("State.jl")

include("compute_Z.jl")
include("prob_miss.jl")
include("loglike.jl")
include("logprior.jl")
include("logq.jl")
include("compute_elbo.jl")

include("fit.jl")

end # VB


text:
using QuasiArrays, Test

include("test_abstractquasiarray.jl")
include("test_quasibroadcast.jl")
include("test_arrayops.jl")
include("test_quasisubarray.jl")

include("test_dense.jl")
include("test_quasiadjtrans.jl")
include("test_continuous.jl")
include("test_matmul.jl")

include("test_quasiconcat.jl")

include("test_ldiv.jl")


text:
using LinearAlgebra
using Random
# using CairoMakie
using Distances
# CairoMakie.activate!(type = "svg")
using StatsBase

# function Merge(y1,y2)
#     N1 = length(y1)
#     N2 = length(y2)
#     N  = N1 + N2
#     i,j,k,b = 1,1,1,0
#     y = zeros(N)
#     while i ≤ N1 && j ≤ N2
#         if y1[i] ≤ y2[j]
#             y[k]  = y1[i]
#             b = b+1
#             i = i+1
#         else
#             y[k] = y2[j]
#             z[y2[j]] = z[y2[j]] + b
#             j = j+1
#         end
#         k = k+1
#     end
#     if i ≤ N1
#         y[k:end] .= y2[k:end]
#     end
# end

function Sort(y)
    N = length(y)
    if N == 1
        return y
    else
        m = floor(Int,N/2)
        y1= Sort(y[1:m])
        y2= Sort(y[m+1:N])
#         y = Merge(y1,y2)
        return y
    end
end

y = [2,4,5,3,1]
Sort(y)
z = zeros(length(y))


text:
include("test_firinterp.jl")
include("test_firfilter.jl")
include("test_firdecim.jl")

max_stars_repo_name: JayKickliter/LiquidDSP

max_issues_repo_name: JayKickliter/LiquidDSP

max_forks_repo_name: JayKickliter/LiquidDSP


text:
# ---
# title: 350. Intersection of Two Arrays II
# id: problem350
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Hash Table, Two Pointers, Binary Search, Sort
# link: <https://leetcode.com/problems/intersection-of-two-arrays-ii/description/>
# hidden: true
# ---
#
# Given two arrays, write a function to compute their intersection.
#
# **Example 1:**
#
#
#
#     Input: nums1 = [1,2,2,1], nums2 = [2,2]
#     Output: [2,2]
#
#
# **Example 2:**
#
#
#
#     Input: nums1 = [4,9,5], nums2 = [9,4,9,8,4]
#     Output: [4,9]
#
# **Note:**
#
#   * Each element in the result should appear as many times as it shows in both arrays.
#   * The result can be in any order.
#
# **Follow up:**
#
#   * What if the given array is already sorted? How would you optimize your algorithm?
#   * What if _nums1_ 's size is small compared to _nums2_ 's size? Which algorithm is better?
#   * What if elements of _nums2_ are stored on disk, and the memory is limited such that you cannot load all elements into the memory at once?
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using ScikitLearnBase
using Random

## 1. Generate synthetic data from two distinct Gaussians: n_samples_A and
##    n_samples_B data points
## 2. Fit the GMM to it
## 3. Count how many points were classified into each Gaussian, test that
##    it's either n_samples_A or n_samples_B

n_samples_A = 300
n_samples_B = 600

# generate spherical data centered on (20, 20)
Random.seed!(42)
shifted_gaussian = randn(n_samples_A, 2) .+ [20, 20]'

# generate twice as many points from zero centered stretched Gaussian data
C = [0. -0.7
     3.5 .7]
stretched_gaussian = randn(n_samples_B, 2) * C

# concatenate the two datasets into the final training set
X_train = vcat(shifted_gaussian, stretched_gaussian)

# fit a Gaussian Mixture Model with two components
gmm = fit!(GMM(n_components=2, kind=:full), X_train)

# Check that the training points are correctly classified
@assert sum(predict(gmm, X_train) .== 1) in [n_samples_A, n_samples_B]


text:
# dual only versus dual only no scaling

using JLD

include("perfplots.jl")

timeshift = 10.0

# produced by process_csv.jl in perfprofile mode
jldsource = joinpath(dirname(@__FILE__),"solvercomparison_perf.jld")

d = load(jldsource)
time_table = d["time_table"] + timeshift
solvers = d["solvers"]
@assert solvers == ["CPLEX_MISOCP","PAJ_CPLEX_MOSEK_msd","PAJ_CPLEX_MOSEK","SCIP_MISOCP","BONMIN","PAJ_CBC_ECOS","PAJ_GLPK_ECOS"]
instance_names = d["instance_names"]

Plots.scalefontsizes(2.1)

performance_profile(time_table[:,1:2], ["CPLEX","MSD CPLEX"],logscale=false, ymax=1.0, xmax=15, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"commercialtime.eps"))

performance_profile(time_table[:,5:7], ["BONMIN","iter. Cbc","iter. GLPK"],logscale=false, ymax=1.0, xmax=15, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"opensourcetime.eps"))

performance_profile(time_table[:,[4,6]], ["iter. Cbc","SCIP_MISOCP"],logscale=false, ymax=1.0, xmax=15, linewidth=3)
Plots.savefig(joinpath(dirname(@__FILE__),"itercbcsciptime.eps"))

max_stars_repo_name: mlubin/PajaritoSupplement

max_issues_repo_name: mlubin/PajaritoSupplement

max_forks_repo_name: mlubin/PajaritoSupplement


text:
# Orchestra module.
module Orchestra

# Load source files
include("system.jl")
include("types.jl")
include("util.jl")
include("transformers.jl")

end # module


text:
module PerformanceRepr

using MIDI, DataStructures
using ..NoteSequences

# Performance representation related constants
const NOTE_ON = 1
const NOTE_OFF = 2
const TIME_SHIFT = 3
const VELOCITY = 4
const DEFAULT_MAX_SHIFT_STEPS = 100

export NOTE_ON, NOTE_OFF, TIME_SHIFT, VELOCITY

include("performance.jl")
include("encoder_decoder.jl")

end


text:
function Base.summary(adf::AbstractDataFrame) # -> UTF8String
    nrows, ncols = size(adf)
    return @sprintf "%dx%d %s" nrows ncols typeof(adf)
end

begin
    local io = IOBuffer(Array(Uint8, 80), true, true)
    global ourstrwidth
    function ourstrwidth(x::Any)
        truncate(io, 0)
        ourshowcompact(io, x)
        return position(io)
    end
end

ourshowcompact(io::IO, x::Any) = showcompact(io, x)
ourshowcompact(io::IO, x::String) = print(io, x)

# Determine the maximum string length of any entry in each DataFrame column
function getmaxwidths(adf::AbstractDataFrame,
                      rowindices1::AbstractVector{Int},
                      rowindices2::AbstractVector{Int},
                      rowlabel::String) # -> Vector{Int}
    ncols = size(adf, 2)
    cnames = names(adf)
    maxwidths = Array(Int, ncols + 1)
    for j in 1:ncols
        # (1) Consider length of column name
        maxwidths[j] = ourstrwidth(cnames[j])

        # (2) Consider length of longest entry in that column
        for i in rowindices1
            maxwidths[j] = max(maxwidths[j], ourstrwidth(adf[i, j]))
        end
        for i in rowindices2
            maxwidths[j] = max(maxwidths[j], ourstrwidth(adf[i, j]))
        end
    end
    rowmaxwidth1 = isempty(rowindices1) ? 0 : ndigits(maximum(rowindices1))
    rowmaxwidth2 = isempty(rowindices2) ? 0 : ndigits(maximum(rowindices2))
    maxwidths[ncols + 1] = max(max(rowmaxwidth1, rowmaxwidth2),
                               ourstrwidth(rowlabel))
    return maxwidths
end

# Determine width of printing DataFrame in a single chunk
function getprintedwidth(maxwidths::Vector{Int}) # -> Int
    n = length(maxwidths)
    # Length of line-initial |
    width = 1
    for i in 1:n
        # Length of field + 2 spaces + trailing |
        width += maxwidths[i] + 3
    end
    return width
end

# Given available space, split the columns into a set of chunks
# each of which should fit on the screen at once
function getchunkbounds(maxwidths::Vector{Int},
                        splitchunks::Bool) # -> Vector{Int}
    ncols = length(maxwidths) - 1
    rowmaxwidth = maxwidths[ncols + 1]
    availablewidth = Base.tty_cols()
    if splitchunks
        chunkbounds = [0]
        # Include 2 spaces + 2 | characters for row/col label
        totalwidth = rowmaxwidth + 4
        for j in 1:ncols
            # Include 2 spaces + | character in per-column character count
            totalwidth += maxwidths[j] + 3
            if totalwidth > availablewidth
                push!(chunkbounds, j - 1)
                totalwidth = rowmaxwidth + 4 + maxwidths[j] + 3
            end
        end
        push!(chunkbounds, ncols)
    else
        chunkbounds = [0, ncols]
    end
    return chunkbounds
end

function showrowindices(io::IO,
                        adf::AbstractDataFrame,
                        rowindices::AbstractVector{Int},
                        rowmaxwidth::Int,
                        maxwidths::Vector{Int},
                        leftcol::Int,
                        rightcol::Int) # -> Nothing
    for i in rowindices
        # Print row ID
        @printf io "| %d" i
        padding = rowmaxwidth - ndigits(i)
        for itr in 1:padding
            write(io, ' ')
        end
        print(io, " | ")
        # Print DataFrame entry
        for j in leftcol:rightcol
            strlen = ourstrwidth(adf[i, j])
            ourshowcompact(io, adf[i, j])
            padding = maxwidths[j] - strlen
            for itr in 1:padding
                write(io, ' ')
            end
            if j == rightcol
                if i == rowindices[end]
                    print(io, " |")
                else
                    print(io, " |\n")
                end
            else
                print(io, " | ")
            end
        end
    end
    return
end

function showrows(io::IO,
                  adf::AbstractDataFrame,
                  rowindices1::AbstractVector{Int},
                  rowindices2::AbstractVector{Int},
                  maxwidths::Vector{Int},
                  splitchunks::Bool = false,
                  rowlabel::String = "Row #",
                  displaysummary::Bool = true) # -> Nothing
    ncols = size(adf, 2)
    cnames = names(adf)

    if displaysummary
        println(io, summary(adf))
    end

    if isempty(rowindices1)
        return
    end

    rowmaxwidth = maxwidths[ncols + 1]
    chunkbounds = getchunkbounds(maxwidths, splitchunks)
    nchunks = length(chunkbounds) - 1

    for chunkindex in 1:nchunks
        leftcol = chunkbounds[chunkindex] + 1
        rightcol = chunkbounds[chunkindex + 1]

        # Print table bounding line
        write(io, '|')
        for itr in 1:(rowmaxwidth + 2)
            write(io, '-')
        end
        write(io, '|')
        for j in leftcol:rightcol
            for itr in 1:(maxwidths[j] + 2)
                write(io, '-')
            end
            write(io, '|')
        end
        write(io, '\n')

        # Print column names
        @printf io "| %s" rowlabel
        padding = rowmaxwidth - ourstrwidth(rowlabel)
        for itr in 1:padding
            write(io, ' ')
        end
        @printf io " | "
        for j in leftcol:rightcol
            s = cnames[j]
            ourshowcompact(io, s)
            padding = maxwidths[j] - ourstrwidth(s)
            for itr in 1:padding
                write(io, ' ')
            end
            if j == rightcol
                print(io, " |\n")
            else
                print(io, " | ")
            end
        end

        # Print main table body, potentially in two abbreviated sections
        showrowindices(io,
                       adf,
                       rowindices1,
                       rowmaxwidth,
                       maxwidths,
                       leftcol,
                       rightcol)
        if !isempty(rowindices2)
            print(io, "\n⋮\n")
            showrowindices(io,
                           adf,
                           rowindices2,
                           rowmaxwidth,
                           maxwidths,
                           leftcol,
                           rightcol)
        end

        # Print newlines to separate chunks
        if chunkindex < nchunks
            print(io, "\n\n")
        end
    end

    return
end

function Base.show(io::IO,
                   adf::AbstractDataFrame,
                   splitchunks::Bool = false,
                   rowlabel::String = "Row #",
                   displaysummary::Bool = true) # -> Nothing
    nrows = size(adf, 1)
    availableheight = Base.tty_rows() - 5
    nrowssubset = fld(availableheight, 2)
    bound = min(nrowssubset - 1, nrows)
    if nrows <= availableheight
        rowindices1 = 1:nrows
        rowindices2 = 1:0
    else
        rowindices1 = 1:bound
        rowindices2 = max(bound + 1, nrows - nrowssubset + 1):nrows
    end
    maxwidths = getmaxwidths(adf, rowindices1, rowindices2, rowlabel)
    width = getprintedwidth(maxwidths)
    if width > Base.tty_cols() && !splitchunks
        column_summary(io, adf)
    else
        showrows(io,
                 adf,
                 rowindices1,
                 rowindices2,
                 maxwidths,
                 splitchunks,
                 rowlabel,
                 displaysummary)
    end
    return
end

function Base.show(adf::AbstractDataFrame,
                   splitchunks::Bool = false) # -> Nothing
    show(STDOUT, adf, splitchunks)
end

function Base.showall(io::IO,
                      adf::AbstractDataFrame,
                      splitchunks::Bool = false,
                      rowlabel::String = "Row #",
                      displaysummary::Bool = true) # -> Nothing
    rowindices1 = 1:size(adf, 1)
    rowindices2 = 1:0
    maxwidths = getmaxwidths(adf, rowindices1, rowindices2, rowlabel)
    width = getprintedwidth(maxwidths)
    showrows(io,
             adf,
             rowindices1,
             rowindices2,
             maxwidths,
             splitchunks,
             rowlabel,
             displaysummary)
    return
end

function Base.showall(adf::AbstractDataFrame,
                      splitchunks::Bool = false) # -> Nothing
    showall(STDOUT, adf, splitchunks)
    return
end

function colmissing(adf::AbstractDataFrame) # -> Vector{Int}
    nrows, ncols = size(adf)
    missing = zeros(Int, ncols)
    for j in 1:ncols
        for i in 1:nrows
            missing[j] += int(isna(adf[i, j]))
        end
    end
    return missing
end

function column_summary(io::IO, adf::AbstractDataFrame) # -> Nothing
    println(io, summary(adf))
    metadata = DataFrame(Name = names(adf),
                         Type = types(adf),
                         Missing = colmissing(adf))
    showall(io, metadata, true, "Col #", false)
    return
end

function Base.print(io::IO, adf::AbstractDataFrame)
    show(io, adf)
    print(io, '\n')
end


text:
using Test, KongYiji, Pkg, JLD2, FileIO, ProgressMeter, DataFrames
#using Knet

include("icwb.jl")

#include("ChTreebank.jl")
include("HMM.jl")
#include("ZhuXian.jl")
#include("LM.jl")
#include("ooHMM.jl")
#include("LMHMM.jl")
#include("NGramHMM.jl")

#=
@testset "Test hmm-lstm..." begin
        #KongYiji.trhmmlstm()
        ctb = ChTreebank()
        #ctb = [[[("pa", "wa")]]]
        m = HMM2(ctb; embedsz=64, hiddensz=64, epochs=1, batchsz=256, seqlen=100, file=joinpath(pathof(KongYiji), "..", "..", "data", "hmm2.jld2"))
        tks = m("一个脱离了低级趣味的人")
        @show tks
end
=#

#=
@testset "Extract wrong cases on icwb..." begin
        home = joinpath("d:\\", "icwb2-data")
        tk = Kong()
        @showprogress 1 "" for source in ["pku", "msr", "cityu", "as"]
                test_file = joinpath(home, "testing", "$(source)_test.utf8")
                truth_file = joinpath(home, "gold", source == "as" ? "as_testing_gold.utf8" : "$(source)_test_gold.utf8")
                dict_file = joinpath(home, "gold", "$(source)_training_words.utf8")
                #tk = Kong(; user_dict_path=dict_file)
                output = tk([line for line in eachline(test_file) if length(line) > 0])
                truth = [split(line) for line in eachline(truth_file) if length(line) > 0]
                debug_file = joinpath(pathof(KongYiji), "..", "..", "test", "wrong.cases.$(source).2.txt")
                open(debug_file, "w") do io
                        for i in 1:length(truth)
                                o, s = output[i], truth[i]
                                ms = eachmatch(o, s)
                                for i = 1:length(ms) - 1
                                        lo, ls = ms[i] .+ 1
                                        ro, rs = ms[i + 1] .- 1
                                        if lo <= ro
                                                print(io, "S : "); for j in ls:rs print(io, s[j], '\t') end; println(io)
                                                print(io, "O : "); for j in lo:ro print(io, o[j], '\t') end; println(io)
                                                println(io, "-------------------------")
                                        end
                                end
                                println(io, "=================================")
                        end
                end
        end
end
=#

#=
@testset "Do statistics on ctb.." begin
        @time ctb_home = KongYiji.unzip7(joinpath(pathof(KongYiji), "..", "..", "data", "ctb.jld2.7z"))
        @time ctb = load(ctb_home)["ctb"]
        df = DataFrame(cw=String[], cpos=String[], lw=String[], rw=String[], lpos=String[], rpos=String[])
        for doc in ctb, sent in doc
                nw = length(sent)
                for i in 1:nw
                        cpos, cw = sent[i]
                        if word == cw | pos == cpos
                                #l = i == 1  ? "^" : sent[i - 1][2][end:end]
                                #r = i == nw ? "\$" : sent[i + 1][2][end:end]
                                lpos, lw = i == 1  ? ("^", "") : sent[i - 1]
                                rpos, rw = i == nw ? ("\$", "") : sent[i + 1]
                                push!(df, (cw, cpos, lw, rw, lpos, rpos))
                        end
                end
        end
        println(size(df))
end
=#

#=
@testset "Testing on icwb..." begin
        home = joinpath("d:\\", "icwb2-data")
        @showprogress 1 "" for source in ["pku", "msr", "cityu", "as"]
                test_file = joinpath(home, "testing", "$(source)_test.utf8")
                truth_file = joinpath(home, "gold", source == "as" ? "as_testing_gold.utf8" : "$(source)_test_gold.utf8")
                dict_file = joinpath(home, "gold", "$(source)_training_words.utf8")
                tk = Kong(; user_dict_path=dict_file)
                output = tk([line for line in eachline(test_file) if length(line) > 0])
                truth = [split(line) for line in eachline(truth_file) if length(line) > 0]
                debug_file = joinpath(pathof(KongYiji), "..", "..", "test", "debug.$(source).txt")
                open(debug_file, "w") do io
                       for i in 1:length(truth)
                               print(io, "S : ")
                               for w in truth[i] print(io, w, "\t") end
                               println(io)
                               print(io, "O : ")
                               for w in output[i] print(io, w, "\t") end
                               println(io)
                               println(io)
                       end
                end
                #=
                dict = Set([word for word in eachline(dict_file) if length(word) > 0])
                println(KongYiji.IcwbScoreTable(dict, truth, output, source))
                println("==================================")
                =#
        end
end
=#

#=
@testset "Testing on icwb..." begin
        home = joinpath("d:\\", "icwb2-data")
        scorer = joinpath(home, "scripts", "score.pl")
        tk = Kong()
        for source in ["as"] #["pku", "msr", "cityu", "as"]
                test = joinpath(home, "testing", "$(source)_test.utf8")
                test_out = joinpath(pathof(KongYiji), "..", "..", "test", "$(source)_kongyiji_tokens.txt")
                open(test_out, "w") do io
                        for line in tk(collect(eachline(test)))
                                for token in line print(io, token, ' ') end
                                println(io)
                        end
                end
                dict = joinpath(home, "gold", "$(source)_training_words.utf8")
                truth = joinpath(home, "gold", source == "as" ? "as_testing_gold.utf8" : "$(source)_test_gold.utf8")
                cmd = "perl -w $(scorer) $(dict) $(truth) $(test_out)"
                run(pipeline(ifelse(Sys.iswindows(), `cmd /c $cmd`, `sh -c $cmd`);
                                stdout=joinpath(pathof(KongYiji), "..", "..", "test", "$(source)_kongyiji_score.txt")))
        end
end
=#

#=
@testset "Generating REQUIRE file..." begin
        println(Pkg.METADATA_compatible_uuid("KongYiji"))
        PT = Pkg.Types
        Pkg.activate("..")             # current directory as the project
        ctx = PT.Context()
        pkg = ctx.env.pkg
        if pkg ≡ nothing
            @error "Not in a package, I won't generate REQUIRE."
            exit(1)
        else
            @info "found package" pkg = pkg
        end

        deps = PT.get_deps(ctx)
        non_std_deps = sort(collect(setdiff(keys(deps), values(ctx.stdlibs))))

        open(joinpath("..", "REQUIRE"), "w") do io
            println(io, "julia 0.7")
            for d in non_std_deps
                println(io, d)
                @info "listing $d"
            end
        end
end
=#

#=
@testset "Generating HMM model file of CTB..." begin
        @time ctb_home = KongYiji.unzip7(joinpath(pathof(KongYiji), "..", "..", "data", "ctb.jld2.7z"))
        @time ctb = load(ctb_home)["ctb"]
        @time hmm = KongYiji.HMM(ctb)
        home = joinpath(pathof(KongYiji), "..", "..", "data", "hmm.jld2")
        mkpath(dirname(home))
        @time @save home hmm
        @time zhome = KongYiji.zip7(home)
        rm(home)
        @time home2 = KongYiji.unzip7(zhome)
        @assert home == home2
        @time hmm2 = load(home2)["hmm"]
        @test hmm == hmm2
end
=#

#=
@testset "Test KongYiji(1) with Hand written examples..." begin
        tk = Kong()
        input = "一个脱离了低级趣味的人"
        output = tk(input)
        @show output

        input = "一/个/脱离/了/低级/趣味/的/人"
        tk(input, "/")

        inputs = [
                "他/说/的/确实/在理",
                "这/事/的确/定/不/下来",
                "费孝通/向/人大/常委会/提交/书面/报告",
                "邓颖超/生前/使用/过/的/物品",
                "停电/范围/包括/沙坪坝区/的/犀牛屙屎/和/犀牛屙屎抽水",
        ]
        println("Input :")
        for input in inputs
                println(input)
        end

        println("raw output :")
        for input in inputs
                println(tk(filter(c -> c != '/', input)))
        end

        tk2 = Kong(; user_dict_array=[("VV", "定"),
                                      ("VA", "在理"),
                                       "邓颖超",
                                       "沙坪坝区",
                                       "犀牛屙屎",
                                       "犀牛屙屎抽水",
                                     ]
        )
        println("output with user dict supplied :")
        for input in inputs
                println(tk2(filter(c -> c != '/', input)))
        end
end

=#


text:
struct STL{Dc,Dp,T} # <: GridTopology
  vertex_to_coordinates::Vector{Point{Dp,T}}
  n_m_to_nface_to_mfaces::Matrix{Table{Int32,Vector{Int32},Vector{Int32}}}
  polytope::Polytope{Dc}
  facedims::Vector{Int8}
  offsets::Vector{Int32}
end

function STL(model::DiscreteModel{Dc,Dp}) where {Dc,Dp}
  Dc == Dp-1 || error("STL must be a surface")
  topo = get_grid_topology(model)
  v_coords = get_vertex_coordinates(topo)
  polytope = only(get_polytopes(topo))
  facedims = get_facedims(topo)
  offsets = Int32.(get_offsets(topo))
  T = Table{Int32,Vector{Int32},Vector{Int32}}
  n_m_to_nf_to_mf = Matrix{T}(undef,Dc+1,Dc+1)
  for n in 0:Dc, m in 0:Dc
    n_m_to_nf_to_mf[n+1,m+1] = get_faces(topo,n,m)
  end
  STL(
    v_coords,
    n_m_to_nf_to_mf,
    polytope,
    facedims,
    offsets)
end

num_dims(stl::STL{Dc}) where Dc = Dc

num_point_dims(stl::STL{Dc,Dp}) where {Dc,Dp} = Dp

get_faces(stl::STL,n,m) = stl.n_m_to_nface_to_mfaces[n+1,m+1]

get_vertex_coordinates(stl::STL) = stl.vertex_to_coordinates

get_polytope(stl::STL) = stl.polytope

get_offsets(stl::STL) = stl.offsets

get_offset(stl::STL,d::Integer) = get_offsets(stl)[d+1]

get_facedims(stl::STL) = stl.facedims

get_face_vertices(stl::STL,d::Integer) = get_faces(stl,d,0)

get_cell_vertices(stl::STL{Dc}) where Dc = get_face_vertices(stl,Dc)

num_faces(stl::STL,d::Integer) = length(get_faces(stl,d,d))

num_faces(stl::STL) = length(get_facedims(stl))

num_vertices(stl::STL) = num_faces(stl,0)

num_cells(stl::STL{Dc}) where Dc = num_faces(stl,Dc)

function get_face_vertices(stl::STL{D}) where D
  @notimplementedif D ≤ 1
  f_to_v = lazy_append(get_face_vertices(stl,0),get_face_vertices(stl,1))
  for d in 2:D
    f_to_v = lazy_append(f_to_v,get_face_vertices(stl,d))
  end
  f_to_v
end

function get_cell!(c,stl::STL{Dc},i::Integer) where Dc
  get_dface!(c,stl,i,Val{Dc}())
end

function get_dface!(c,stl::STL,i::Integer,::Val{d}) where d
  T = get_face_vertices(stl,d)
  X = get_vertex_coordinates(stl)
  p = get_polytope(stl)
  get_dface!(c,T,X,p,i,Val{d}())
end

function get_cell(stl::STL,i::Integer)
  c = get_cell_cache(stl)
  get_cell!(c,stl,i)
end

function get_cell_cache(stl::STL{Dc}) where Dc
  get_dface_cache(stl,Dc)
end

function get_dface_cache(stl::STL,d::Integer)
  T = get_face_vertices(stl,d)
  array_cache(T)
end

function get_cell!(cache,T,X,p::Polytope{D},i::Integer) where D
  get_dface!(cache,T,X,p,i,Val{D}())
end

function get_dface!(cache,T,X,p::Polytope,i::Integer,::Val{d}) where d
  @notimplementedif !is_simplex(p)
  get_simplex_dface!(cache,T,X,i,Val{d}())
end

function get_simplex_dface!(cache,T,X,i::Integer,::Val{d}) where d
  nodes = getindex!(cache,T,i)
  vertices = ntuple( i -> X[nodes[i]], Val{d+1}() )
  simplex_face( vertices )
end

function is_open_surface(stl::STL)
  Dc = num_dims(stl)
  e_to_f = get_faces(stl,Dc-1,Dc)
  1 ≤ maximum(length,e_to_f) ≤ 2 || return false
  1 ≤ minimum(length,e_to_f) ≤ 2 || return false
  true
end

# Read STL from file

function read_stl(filename::String)
  F = _file_format(filename)
  mesh = load(File{F}(filename))
  vertex_to_coordinates = MeshIO.decompose(MeshIO.Point3,mesh)
  facet_to_vertices = MeshIO.decompose(MeshIO.TriangleFace,mesh)
  vertex_to_normals = MeshIO.decompose_normals(mesh)
  vertex_to_coordinates = Vector{Point{3,Float64}}( vertex_to_coordinates )
  facet_to_vertices = Table( Vector{Vector{Int}}( facet_to_vertices ) )
  vertex_to_normals = Vector{VectorValue{3,Float64}}( vertex_to_normals )
  facet_to_normals = vertex_to_normals[1:3:end]
  vertex_to_coordinates,facet_to_vertices,facet_to_normals
end

function _file_format(filename)
  ext =  lowercase(last(splitext(filename)))
  stl,obj = ".stl",".obj"
  ext == stl && !_is_binary(filename) && return format"STL_ASCII"
  ext == stl && _is_binary(filename) && return format"STL_BINARY"
  ext == obj && return format"OBJ"
  msg = "Invalid file format: $filename"
  @unreachable msg
end

function _get_encoding(filename)
  cmd = `file --mime-encoding $filename`
  read(cmd,String)
end

_is_binary(filename) = contains(_get_encoding(filename),"binary")

function Base.convert(::Type{Point{D,T}},x::MeshIO.Point{D}) where {D,T}
  Point{D,T}(x.data)
end

function Base.convert(::Type{Vector{T}},x::MeshIO.TriangleFace) where T
  convert.(T,collect(x))
end

# Save STL to file

function save_as_stl(stl::DiscreteModel{Dc,Dp},filename) where {Dc,Dp}
  (Dc == 2 && Dp == 3 ) || error("Geometry incompatible with STL format")
  filename *= ".stl"
  file = File{format"STL_BINARY"}(filename)
  facetype = MeshIO.GLTriangleFace
  pointtype = MeshIO.Point3f0
  normaltype = MeshIO.Vec3f0
  faces = Array{facetype}(undef, num_cells(stl))
  vertices = Array{pointtype}(undef, num_cells(stl) * 3)
  normals = Array{normaltype}(undef, num_cells(stl) * 3)
  for i in 1:num_cells(stl)
    faces[i] = facetype( (i-1)*3+1,(i-1)*3+2, (i-1)*3+3 )
    n = normal(get_cell(stl,i))
    normals[ (i-1)*3+1 ] = normaltype(Tuple(n)...)
    normals[ (i-1)*3+2 ] = normals[ (i-1)*3+1 ]
    normals[ (i-1)*3+3 ] = normals[ (i-1)*3+1 ]
    coords = get_cell_coordinates(stl)[i]
    vertices[ (i-1)*3+1 ] = pointtype(Tuple(coords[1])...)
    vertices[ (i-1)*3+2 ] = pointtype(Tuple(coords[2])...)
    vertices[ (i-1)*3+3 ] = pointtype(Tuple(coords[3])...)
  end
  mesh = MeshIO.Mesh(MeshIO.meta(vertices; normals=normals), faces)
  save(file,mesh)
  filename
end

function save_as_stl(stls::Vector{<:DiscreteModel},filename)
  files = []
  for (i,stl) in enumerate(stls)
    f = save_as_stl(stl,filename*"_$i")
    push!(files,f)
  end
  files
end

# STL as Grid and DiscreteModel

function get_cell_cache(grid::Grid)
  T = get_cell_node_ids(grid)
  array_cache(T)
end

function get_cell!(c,grid::Grid,i)
  T = get_cell_node_ids(grid)
  X = get_node_coordinates(grid)
  p = get_polytope( only(get_reffes(grid)) )
  get_cell!(c,T,X,p,i)
end
function get_facet_cache(model::DiscreteModel)
  Dc = num_dims(model)
  get_dface_cache(model,Dc-1)
end

function get_dface_cache(model::DiscreteModel,d)
  topo = get_grid_topology(model)
  T = get_face_vertices(topo,d)
  array_cache(T)
end

function get_facet!(c,model::DiscreteModel{Dc},i) where Dc
  get_dface!(c,model,i,Val{Dc-1}())
end

function get_dface!(c,model::DiscreteModel,i,::Val{d}) where d
  topo = get_grid_topology(model)
  T = get_face_vertices(topo,d)
  X = get_vertex_coordinates(topo)
  p = get_polytope( only(get_reffes(model)) )
  get_dface!(c,T,X,p,i,Val{d}())
end

function get_dface(model::DiscreteModel,i,::Val{d}) where d
  c = get_dface_cache(model,d)
  get_dface!(c,model,i,Val{d}())
end

function get_vertex(stl::DiscreteModel,vertex::Integer)
  get_vertex_coordinates(get_grid_topology(stl))[vertex]
end

function get_edge(stl::DiscreteModel,edge::Integer)
  get_dface(stl,edge,Val{1}())
end

function get_facet(stl::DiscreteModel{Dc},facet::Integer) where Dc
  get_dface(stl,facet,Val{Dc-1}())
end

function get_cell(stl::DiscreteModel{Dc},cell::Integer) where Dc
  get_dface(stl,cell,Val{Dc}())
end

function Base.eps(T::Type{<:AbstractFloat},grid::Grid)
  pmin,pmax = get_bounding_box(grid)
  vmax = max(abs.(Tuple(pmin))...,abs.(Tuple(pmax))...)
  eps(T(vmax))
end

Base.eps(grid::Grid) = eps(Float64,grid)

Base.eps(model::DiscreteModel) = eps(get_grid(model))

function compute_stl_model(
  cell_to_vertices::Table,
  vertex_to_coordinates::Vector{<:Point{D}}) where D

  p = Polytope( tfill(TET_AXIS,Val{D-1}()) )
  compute_model(cell_to_vertices,vertex_to_coordinates,p)
end

function compute_model(
  cell_to_vertices::Table,
  vertex_to_coordinates::Vector{<:Point},
  p::Polytope)

  grid = compute_grid(cell_to_vertices,vertex_to_coordinates,p)
  UnstructuredDiscreteModel(grid)
end

function compute_grid(
  cell_to_nodes::AbstractArray,
  node_to_coordinates::Vector{<:Point},
  p::Polytope)

  T = Table(cell_to_nodes)
  X = node_to_coordinates
  reffes = [LagrangianRefFE(p)]
  cell_types = fill(1,length(T))
  UnstructuredGrid(X,T,reffes,cell_types)
end

function merge_and_collapse(stl::DiscreteModel;atol=eps(Float32,stl))
  stl = merge_nodes(stl;atol)
  collapse_small_facets!(stl;atol)
  stl = merge_nodes(stl;atol)
  stl = delete_duplicated_faces(stl)
  preprocess_small_facets(stl;atol)
end

function collapse_small_facets!(stl::DiscreteModel;atol)
  D = num_dims(stl)
  stl_topo = get_grid_topology(stl)
  coords = get_vertex_coordinates(stl_topo)
  f_to_v = get_faces(stl_topo,D,0)
  f_to_e = get_faces(stl_topo,D,D-1)
  e_to_v = get_faces(stl_topo,D-1,0)
  fv_c = array_cache(f_to_v)
  fe_c = array_cache(f_to_e)
  ev_c = array_cache(e_to_v)
  c = get_facet_cache(stl)
  for f in 1:num_cells(stl)
    for v in getindex!(fv_c,f_to_v,f)
      vertex = coords[v]
      for e in getindex!(fe_c,f_to_e,f)
        vertices = getindex!(ev_c,e_to_v,e)
        v ∉ vertices || continue
        edge = get_facet!(c,stl,e)
        if distance(vertex,edge) < atol
          _vertex = projection(vertex,edge)
          coords[v] = _vertex
        end
      end
    end
  end
  stl
end

function merge_nodes(stl::DiscreteModel;atol=eps(Float32,stl))
  X,T = delete_repeated_vertices(stl;atol)
  compute_stl_model(T,X)
end

function delete_repeated_vertices(stl::DiscreteModel;atol)
  vertex_to_equal = _map_equal_vertices(stl;atol)
  topo = get_grid_topology(stl)
  X = get_vertex_coordinates(topo)
  T = get_cell_vertices(topo)
  X,T = _delete_vertices(X,T,vertex_to_equal)
  _delete_empty_cells!(T)
  T = Table(T)
  X,T
end

function _map_equal_vertices(stl::DiscreteModel;atol)
  _map_equal_vertices_from_cloud(stl;atol)
end

function _map_equal_vertices_from_cloud(stlmodel::DiscreteModel;atol)
  stl = get_grid(stlmodel)
  group_to_vertices =  _group_vertices(stl;atol=atol*10)
  vertices_map = collect(1:num_nodes(stl))
  for _vertices in group_to_vertices
    for i in 1:length(_vertices), j in i+1:length(_vertices)
      if _vertices[i] != _vertices[j]
        vi = get_node_coordinates(stl)[_vertices[i]]
        vj = get_node_coordinates(stl)[_vertices[j]]
        if vertices_map[_vertices[i]] == _vertices[i] &&
           vertices_map[_vertices[j]] == _vertices[j]

          if distance(vi,vj) < atol
            vertices_map[_vertices[j]] = vertices_map[_vertices[i]]
          end
        end
      end
    end
  end
  for i in 1:length(vertices_map)
    nmax = 10
    n = 0
    v = i
    while vertices_map[v] != v
      v = vertices_map[v]
      n += 1
      @assert n < nmax
    end
    if vertices_map[i] != v
      vertices_map[i] = v
    end
  end
  vertices_map
end

function _delete_vertices(X,T,vertex_to_equal_vertex)
  u = unique(vertex_to_equal_vertex)
  m = vertex_to_equal_vertex .== 1:length(vertex_to_equal_vertex)
  m = cumsum(m) .* m
  vertices_map = m[vertex_to_equal_vertex]
  @assert count(iszero,vertices_map) == 0
  X = X[u]
  T = map( i -> vertices_map[i], T )
  X,T
end

function _delete_empty_cells!(T)
  filter!( f -> length(f) == _num_uniques(f) , T )
end

function _num_uniques(a::AbstractVector)
  c = 0
  for i in 1:length(a)
    if a[i] ∉ view(a,1:i-1)
      c += 1
    end
  end
  c
end

function _group_vertices(stl::Grid{Dc,D};atol) where {Dc,D}
  min_length = _compute_min_length(stl;atol)
  digits = -Int(floor(log10(min_length)))
  pmin, pmax = get_bounding_box(stl)
  cells = NTuple{D,Int}[]
  vertices = Int[]
  m = CartesianIndices( tfill(2,Val{D}()) )
  for (iv,v) in enumerate(get_node_coordinates(stl))
    p = v-pmin
    for i in m
      f = d -> i.I[d] == 1 ? floor(p[d];digits) : ceil(p[d];digits)
      c = ntuple( f, Val{D}() )
      c = Int.(round.(exp10(digits).*c)).+1
      push!(cells,c)
      push!(vertices,iv)
    end
  end
  u = unique(cells)
  dict = Dict( zip(u,1:length(u)) )
  cells = map( i -> dict[i], cells )
  cell_to_vertices = [ Int[] for _ in 1:maximum(cells) ]
  for (cell,vertex) in zip(cells,vertices)
    push!(cell_to_vertices[cell],vertex)
  end
  cell_to_vertices
end

function _compute_min_length(mesh::Grid;atol)
  min_length = Inf
  for nodes in get_cell_node_ids(mesh)
    for i in 1:length(nodes), j in i+1:length(nodes)
      vi = get_node_coordinates(mesh)[nodes[i]]
      vj = get_node_coordinates(mesh)[nodes[j]]
      dist = distance(vi,vj)
      if dist < min_length && dist > atol
        min_length = dist
      end
    end
  end
  min_length
end

function preprocess_small_facets(stl::DiscreteModel;atol)
  @notimplementedif !is_edge_manifold(stl) "The geometry is not edge manifold"
  @notimplementedif !is_water_tight(stl) "The geometry is not watter tight"
  max_iters = 100
  for i in 1:max_iters
    stl,incomplete = _preprocess_small_facets(stl;atol)
    incomplete || return stl
  end
  @warn "Unable to fix small facets in $max_iters iterations"
  println("@warn: Unable to fix small facets in $max_iters iterations")
  stl
end

function _preprocess_small_facets(stl::DiscreteModel{Dc};atol) where Dc
  @notimplementedif Dc ≠ 2
  hang_v,cut_e = get_hanging_vertices_and_edges(stl;atol)
  _stl = stl
  incomplete = false
  if !isempty(hang_v)
    incomplete = true
    touched = falses(num_cells(stl))
    topo = get_grid_topology(stl)
    f_to_v = get_faces(topo,Dc,0)
    e_to_v = get_faces(topo,Dc-1,0)
    e_to_f = get_faces(topo,Dc-1,Dc)
    fv = array_cache(f_to_v)
    ev = array_cache(e_to_v)
    ef = array_cache(e_to_f)
    c = get_facet_cache(stl)
    X = get_vertex_coordinates(topo)
    T = Vector( f_to_v )
    Tnew = Vector{Int32}[]
    for (v,e) in zip(hang_v,cut_e)
      edge = get_facet!(c,stl,e)
      p = projection(X[v],edge)
      X[v] = p
      faces_around = getindex!(ef,e_to_f,e)
      if any(f->touched[f],faces_around)
        continue
      end
      for f in faces_around
        touched[f] = true
        ks = split_face!(fv,ev,f_to_v,e_to_v,f,e,v)
        push!(Tnew,ks...)
      end
    end
    deleteat!(T,touched)
    append!(T,Tnew)
    _delete_empty_cells!(T)
    _stl = compute_stl_model(Table(T),X)
    _stl = merge_nodes(_stl;atol)
    _stl = delete_duplicated_faces(_stl)
    if !is_water_tight(_stl)
      _stl = stl
      incomplete = false
    end
  end
  _stl,incomplete
end

function delete_duplicated_faces(stl::DiscreteModel)
  topo = get_grid_topology(stl)
  D = num_dims(stl)
  v_to_f = get_faces(topo,0,D)
  f_to_v = get_faces(topo,D,0)
  vf = array_cache(v_to_f)
  fvi = array_cache(f_to_v)
  fvj = array_cache(f_to_v)
  is_duplicated = falses(num_cells(stl))
  for v in 1:num_vertices(topo)
    facets = getindex!(vf,v_to_f,v)
    for i in 1:length(facets), j in i+1:length(facets)
      fi = facets[i]
      fj = facets[j]
      fi ≠ fj || continue
      vi = getindex!(fvi,f_to_v,fi)
      vj = getindex!(fvj,f_to_v,fj)
      if vi ⊆ vj
        is_duplicated[fi] = is_duplicated[fj] = true
      end
    end
  end
  facets = findall(d->!d,is_duplicated)
  T = f_to_v[facets]
  X = get_vertex_coordinates(topo)
  compute_stl_model(T,X)
end

function split_face!(fv,ev,f_to_v,e_to_v,f,e,v)
  facet = getindex!(fv,f_to_v,f)
  edge = getindex!(ev,e_to_v,e)
  replace( facet, edge[1] => v ),
  replace( facet, edge[2] => v )
end

function get_hanging_vertices_and_edges(stl::DiscreteModel{Dc};atol) where Dc
  @notimplementedif Dc ≠ 2
  topo = get_grid_topology(stl)
  vertex_coordinates = get_vertex_coordinates(topo)
  f_to_e = get_faces(topo,Dc,Dc-1)
  f_to_v = get_faces(topo,Dc,0)
  e_to_v = get_faces(topo,Dc-1,0)
  ce = array_cache(f_to_e)
  cv = array_cache(f_to_v)
  _cv = array_cache(e_to_v)
  c = get_facet_cache(stl)
  hanging_vertices = Int32[]
  cut_edges = Int32[]
  for f in 1:num_cells(stl)
    min_dist = Inf
    hang_v = UNSET
    cut_e = UNSET
    for e in getindex!(ce,f_to_e,f)
      for v in getindex!(cv,f_to_v,f)
        _vertices = getindex!(_cv,e_to_v,e)
        v ∉ _vertices || continue
        edge = get_facet!(c,stl,e)
        point = vertex_coordinates[v]
        dist = distance(edge,point)
        if dist < min_dist
          min_dist = dist
          hang_v = v
          cut_e = e
        end
      end
    end
    if min_dist < atol
      push!(hanging_vertices,hang_v)
      push!(cut_edges,cut_e)
    end
  end
  hanging_vertices,cut_edges
end

function split_disconnected_parts(stl::DiscreteModel)
  Dc = num_dims(stl)
  f_to_v = get_faces(get_grid_topology(stl),Dc,0)
  v_to_f = get_faces(get_grid_topology(stl),0,Dc)
  v_to_part = fill(UNSET,num_vertices(stl))
  f_to_part = fill(UNSET,num_cells(stl))
  num_parts = 0
  stack = Int[]
  for f in 1:num_cells(stl)
    f_to_part[f] == UNSET || continue
    num_parts += 1
    f_to_part[f] = num_parts
    empty!(stack)
    push!(stack,f)
    while !isempty(stack)
      f_curr = pop!(stack)
      for v in f_to_v[f_curr]
        v_to_part[v] = num_parts
        for f_neig in v_to_f[v]
          f_curr ≠ f_neig || continue
          f_to_part[f_neig] == UNSET || continue
          f_to_part[f_neig] = num_parts
          push!(stack,f_neig)
        end
      end
    end
  end
  @assert all(!isequal(UNSET),v_to_part)
  @assert all(!isequal(UNSET),f_to_part)
  coords = get_vertex_coordinates(get_grid_topology(stl))
  stls = typeof(stl)[]
  for part in 1:num_parts
    faces = findall(isequal(part),f_to_part)
    vertices = findall(isequal(part),v_to_part)
    v_to_part_v = fill(UNSET,length(v_to_part))
    v_to_part_v[vertices] = 1:length(vertices)
    _f_to_v = Table(map(i->v_to_part_v[i],f_to_v[faces]))
    _coords = coords[vertices]
    _stl = compute_stl_model(_f_to_v,_coords)
    push!(stls,_stl)
  end
  stls
end

function check_requisites(stl::DiscreteModel,bgmodel::DiscreteModel;
  verbose=true,max_num_facets=10000)

  !verbose || println(join(fill('-',40)))
  fulfill = true
  if !is_surface(stl)
    !verbose || println("Is not a surface")
    fulfill = false
  end
  if !is_vertex_manifold(stl)
    !verbose || println("Is not vertex manifold")
    fulfill = false
  end
  if !is_edge_manifold(stl)
    !verbose || println("Is not edge manifold")
    fulfill = false
  end
  if !is_water_tight(stl)
    !verbose || println("Is not water tight")
    fulfill = false
  end
  if has_sharp_edges(stl)
    !verbose || println("Has sharp edges")
    fulfill = false
  end
  max_fv = max_num_facets_per_vertex(stl)
  max_fbc = max_num_facets_per_bgcell(stl,bgmodel)
  println("Maximum num facets per vertex: $max_fv")
  println("Maximum num facets per bgcell: $max_fbc")
  if max_fbc > max_num_facets
    fulfill = false
    if max_fv > max_num_facets
      !verbose || println("Unable to run geometry $max_num_facets")
    else
      !verbose || println("Please refine your mesh")
    end
  end
  fulfill
end

function is_vertex_manifold(stlmodel::DiscreteModel{2,3})
  stl = get_grid_topology(stlmodel)
  D = num_dims(stl)
  v_to_e = get_faces(stl,0,1)
  v_to_f = get_faces(stl,0,D)
  f_to_v = get_faces(stl,D,0)
  fc = array_cache(v_to_f)
  vc = array_cache(f_to_v)
  for v in 1:num_vertices(stl)
    vfacets = getindex!(fc,v_to_f,v)
    !isempty(vfacets) || continue
    f0 = vfacets[1]
    fnext = f0
    nf = 0
    while true
      nf += 1
      i = findfirst(isequal(v), getindex!(vc,f_to_v,fnext) )
      inext = i == D+1 ? 1 : i+1
      vnext = getindex!(vc,f_to_v,fnext)[inext]
      faces = getindex!(fc,v_to_f,vnext)
      i = findfirst( f -> f ≠ fnext && v ∈ getindex!(vc,f_to_v,f), faces )
      fnext = isnothing(i) ? UNSET : faces[i]
      fnext ≠ UNSET || break
      fnext ≠ f0 || break
    end
    if nf ≠ length(vfacets)
      return false
    end
  end
  true
end

function is_edge_manifold(stlmodel::DiscreteModel{2,3})
  stl = get_grid_topology(stlmodel)
  e_to_f = get_faces(stl,1,2)
  c = array_cache(e_to_f)
  maximum(e->length(getindex!(c,e_to_f,e)),1:length(e_to_f)) ≤ 2
end

function max_num_facets_per_bgcell(stlmodel,bgmodel)
  stl = STL(stlmodel)
  grid = get_grid(bgmodel)
  c_to_stlf =compute_cell_to_facets(grid,stl)
  maximum(length,c_to_stlf)
end

function max_num_facets_per_vertex(stlmodel)
  stl = get_grid_topology(stlmodel)
  v_to_f = get_faces(stl,0,2)
  c = array_cache(v_to_f)
  maximum(v->length(getindex!(c,v_to_f,v)),1:length(v_to_f))
end

function has_sharp_edges(stlmodel)
  stl = STL(stlmodel)
  Π = get_facet_planes(stl)
  e_to_f = get_faces(stl,1,2)
  c = array_cache(e_to_f)
  for e in 1:length(e_to_f)
    facets = getindex!(c,e_to_f,e)
    Π1 = Π[facets[1]]
    Π2 = Π[facets[2]]
    n1 = normal(Π1)
    n2 = normal(Π2)
    n1 ⋅ n2 ≉ -1 || return true
  end
  false
end

function get_bounding_box(mesh::Grid)
  pmin = get_node_coordinates(mesh)[1]
  pmax = get_node_coordinates(mesh)[1]
  for vertex in get_node_coordinates(mesh)
    pmin = Point(min.(Tuple(pmin),Tuple(vertex)))
    pmax = Point(max.(Tuple(pmax),Tuple(vertex)))
  end
  pmin,pmax
end

function get_bounding_box(model::DiscreteModel)
  get_bounding_box(get_grid_topology(model))
end

function get_bounding_box(msh::T) where T<:Union{GridTopology,STL}
  vertices = get_vertex_coordinates(msh)
  pmin = pmax = vertices[1]
  for vertex in vertices
    pmin = Point( min.( Tuple(pmin), Tuple(vertex) ) )
    pmax = Point( max.( Tuple(pmax), Tuple(vertex) ) )
  end
  pmin,pmax
end

is_water_tight(model::DiscreteModel) = is_water_tight(get_grid_topology(model))

function is_water_tight(top::GridTopology{Dc,Dp}) where {Dc,Dp}
  Dc == Dp - 1 || error("The Geometry is not a surface")
  c_to_f = get_faces(top,Dc-1,Dc)
  maximum(length,c_to_f) == minimum(length,c_to_f) == 2
end

is_surface(model::DiscreteModel) = is_surface(get_grid_topology(model))

function is_surface(::GridTopology{Dc,Dp}) where {Dc,Dp}
  Dc == Dp-1
end

function measure(a::Grid,mask)
  m = BigFloat(0.0,precision=128)
  p = get_polytope(only(get_reffes(a)))
  T = get_cell_node_ids(a)
  X = get_node_coordinates(a)
  c = array_cache(T)
  for i in 1:num_cells(a)
    if mask[i]
      f = get_cell!(c,T,X,p,i)
      m += measure(f)
    end
  end
  Float64(m)
end

function measure(a::Grid,cell_to_inoutcut,inoutcut::Symbol)
  inout_dict = Dict{Symbol,Int8}(
    :in => FACE_IN,
    :out => FACE_OUT,
    :cut => FACE_CUT )
  measure(a,cell_to_inoutcut,inout_dict[inoutcut])
end

function measure(a::Grid,cell_to_val,val)
  measure(a,lazy_map(isequal(val),cell_to_val))
end

function measure(a::Grid)
  measure(a,lazy_map(i->true,1:num_cells(a)))
end

volume(a::Grid{D,D},args...) where D = measure(a,args...)

function surface(a::Grid{Df,Dp},args...) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  measure(a,args...)
end

surface(a::DiscreteModel) = surface(get_grid(a))

function measures(a::Grid,cell_to_inoutcut)
  measure(a,cell_to_inoutcut,FACE_IN),
  measure(a,cell_to_inoutcut,FACE_OUT),
  measure(a,cell_to_inoutcut,FACE_CUT)
end

volumes(a::Grid{D,D},args...) where D = measures(a,args...)

function surfaces(a::Grid{Df,Dp},args...) where {Df,Dp}
  @notimplementedif Df ≠ Dp-1
  measures(a,args...)
end

function min_height(grid::Grid)
  min_h = Inf
  c = get_cell_cache(grid)
  for i in 1:num_cells(grid)
    cell = get_cell!(c,grid,i)
    h = min_height(cell)
    if h < min_h
      min_h = h
    end
  end
  min_h
end

min_height(model::DiscreteModel) = min_height(get_grid(model))

# Cartesian Grid

function get_bounding_box(grid::CartesianGrid)
  desc = get_cartesian_descriptor(grid)
  pmin = desc.origin
  pmax = desc.origin + VectorValue(desc.partition .* desc.sizes)
  pmin,pmax
end

function measure(a::CartesianGrid)
  pmin,pmax = get_bounding_box(a)
  measure(pmin,pmax)
end

function measure(a::CartesianGrid{D,T,typeof(identity)},mask) where {D,T}
  cell_vol = prod( get_cartesian_descriptor(a).sizes )
  cell_vol * count(mask)
end

function compute_cell_to_facets(grid::CartesianGrid,stl::STL)
  CELL_EXPANSION_FACTOR = 1e-3
  desc = get_cartesian_descriptor(grid)
  @assert length(get_reffes(grid)) == 1
  p = get_polytope(get_cell_reffe(grid)[1])
  @notimplementedif desc.map !== identity
  cell_to_stl_facets = [ Int32[] for _ in 1:num_cells(grid) ]
  n = Threads.nthreads()
  thread_to_cells = [ Int32[] for _ in 1:n ]
  thread_to_stl_facets = [ Int32[] for _ in 1:n ]
  coords = get_node_coordinates(grid)
  cell_to_nodes = get_cell_node_ids(grid)
  c = [ ( get_cell_cache(stl), array_cache(cell_to_nodes) ) for _ in 1:n ]
  δ = CELL_EXPANSION_FACTOR
  Threads.@threads for stl_facet in 1:num_cells(stl)
    i = Threads.threadid()
    cc,nc = c[i]
    f = get_cell!(cc,stl,stl_facet)
    pmin,pmax = get_bounding_box(f)
    for cid in get_cells_around(desc,pmin,pmax)
      cell = LinearIndices(desc.partition)[cid.I...]
      nodes = getindex!(nc,cell_to_nodes,cell)
      _pmin = coords[nodes[1]]
      _pmax = coords[nodes[end]]
      Δ = (_pmax - _pmin) * δ
      _pmin = _pmin - Δ
      _pmax = _pmax + Δ
      if voxel_intersection(f,_pmin,_pmax,p)
        push!(thread_to_cells[i],cell)
        push!(thread_to_stl_facets[i],stl_facet)
      end
    end
  end
  cell_to_stl_facets = [ Int32[] for _ in 1:num_cells(grid) ]
  for (cells,stl_facets) in zip(thread_to_cells,thread_to_stl_facets)
    for (cell,stl_facet) in zip(cells,stl_facets)
      push!(cell_to_stl_facets[cell],stl_facet)
    end
  end
  cell_to_stl_facets
end

function get_cells_around(desc::CartesianDescriptor{D},pmin::Point,pmax::Point) where D
  cmin,_ = get_cell_bounds(desc,pmin)
  _,cmax = get_cell_bounds(desc,pmax)
  cmin = CartesianIndices(desc.partition)[cmin]
  cmax = CartesianIndices(desc.partition)[cmax]
  ranges = ntuple( i -> cmin.I[i] : cmax.I[i], Val{D}() )
  CartesianIndices( ranges )
end

function get_cell_bounds(desc::CartesianDescriptor,p::Point)
  function _get_cell(cell)
    cell = Int.(cell)
    cell = max.(cell,1)
    cell = min.(cell,desc.partition)
    LinearIndices(desc.partition)[cell...]
  end
  tol = 0.1
  coords = Tuple(p-desc.origin)./desc.sizes
  cell = floor.(coords).+1
  cell_min = cell .- ( (coords.-(floor.(coords))) .< tol )
  cell_max = cell .+ ( (coords.-(floor.(coords))) .> (1-tol) )
  _get_cell(cell_min),_get_cell(cell_max)
end

# Bisectors

function bisector_plane!(
   cache,
   stl::STL{Dc,Dp},
   d::Integer,
   dface::Integer,
   Πf::AbstractArray) where {Dc,Dp}

  @notimplementedif Dc ≠ Dp-1
  @notimplementedif d ≠ Dc-1
  c,fc = cache
  e_to_f = get_faces(stl,d,Dc)
  facets = getindex!(c,e_to_f,dface)
  length(facets) == 2 || return Πf[ only(facets) ]
  edge = get_dface!(fc,stl,dface,Val{Dc-1}())
  Π1 = Πf[ facets[1] ]
  Π2 = Πf[ facets[2] ]
  bisector_plane(edge,Π1,Π2)
end

function bisector_plane(
   stl::STL,
   d::Integer,
   dface::Integer,
   Πf::AbstractArray)

  c = bisector_plane_cache(stl,d)
  bisector_plane(c,stl,d,dface,Πf)
end

function bisector_plane_cache(stl::STL,d::Integer)
  Dc = num_dims(stl)
  e_to_f = get_faces(stl,d,Dc)
  c = array_cache(e_to_f)
  fc = get_dface_cache(stl,d)
  c,fc
end


text:
#' ---
#' title: Intro to Julia
#' ---

#' **Originally Contributed by**: Juan Pablo Vielma

#' ## Introduction

#' Since JuMP is embedded in Julia, knowing some basic Julia is important
#' for learning JuMP. This notebook is designed to provide a minimalist
#' crash course in the basics of Julia. You can find resources that provide
#' a more comprehensive introduction to Julia [here](https://julialang.org/learning/).

#' ### How to Print
#' In Julia, we usually use println() to print

println("Hello, World!")

#' ### Basic Data Types
#' Integers

typeof(1 + -2)

#' Floating point numbers

typeof(1.2 - 2.3)

#' There are also some cool things like an irrational representation of π. To make π
#' (and most other greek letters), type \pi and then press [TAB].

π

#+

typeof(π)

#' Julia has native support for complex numbers

typeof(2 + 3im)

#' Double quotes are used for strings

typeof("This is Julia")

#' Unicode is fine in strings

typeof("π is about 3.1415")

#'Julia symbols provide a way to make human readable unique identifiers.

:my_id
typeof(:my_id)

#' ### Arithmetic and Equality Testing

#' Julia is great for math

1 + 1

#'Even math involving complex numbers

(2 + 1im) * (1 - 2im)

#' We can also write things like the following using √ (\sqrt)

sin(2π/3) == √3/2

#' Wait. What???

sin(2π/3) - √3/2

#' Let's try again using ≈ (\approx).

sin(2π/3) ≈ √3/2

#' Note that this time we used ≈ instead of ==. That is because computers don't use
#' real numbers. They use a discrete representation called floating point. If you aren't
#' careful, this can throw up all manner of issues. For example:

1 + 1e-16 == 1

#' It even turns out that floating point numbers aren't associative!

(1 + 1e-16) - 1e-16 == 1 + (1e-16 - 1e-16)

#' ### Vectors, Matrices and Arrays
#' Similar to Matlab, Julia has native support for vectors, matrices and tensors; all of which are represented by arrays of different dimensions.
#' Vectors are constructed by comma-separated elements surrounded by square brackets:

b = [5, 6]

#' Matrices can by constructed with spaces separating the columns, and semicolons separating the rows:

A = [1 2; 3 4]

#' We can do linear algebra:

x = A \ b

#+

A * x

#+

A * x == b

#' Note that when multiplying vectors and matrices, dimensions matter. For example, you can't multiply a vector by a vector:
#+ tangle = false

b * b

#' But multiplying transposes works:

@show b' * b
@show b * b';

#' ### Tuples
#' Julia makes extensive use of a simple data structure called Tuples.  Tuples are immutable collections of values.
#' For example,

t = ("hello", 1.2, :foo)

#+

typeof(t)

#' Tuples can be accessed by index, similar to arrays,

t[2]

#' And can be "unpacked" like so,

a, b, c = t
b

#' The values can also be given names, which is a convenient way of making light-weight data structures.

t = (word="hello", num=1.2, sym=:foo)

#' Then values can be accessed using a dot syntax,

t.word

#' ### Dictionaries
#' Similar to Python, Julia has native support for dictionaries.  Dictionaries provide a very generic way of mapping keys to values.  For example, a map of integers to strings,

d1 = Dict(1 => "A", 2 => "B", 4 => "D")

#' Looking up a values uses the bracket syntax,

d1[2]

#' Dictionaries support non-integer keys and can mix data types,

Dict("A" => 1, "B" => 2.5, "D" => 2 - 3im)

#' Dictionaries can be nested

d2 = Dict("A" => 1, "B" => 2, "D" => Dict(:foo => 3, :bar => 4))

#+

d2["B"]

#+

d2["D"][:foo]

#' ### For-Each Loops
#' Julia has native support for for-each style loops with the syntax `for <value> in <collection> end`.

for i in 1:5
    println(i)
end

#+

for i in [1.2, 2.3, 3.4, 4.5, 5.6]
    println(i)
end

#' This for-each loop also works with dictionaries.

for (key, value) in Dict("A" => 1, "B" => 2.5, "D" => 2 - 3im)
    println("$key: $value")
end

#' Note that in contrast to vector languages like Matlab and R, loops do not result in a significant performance degradation in Julia.

#' ### Control Flow
#' Julia control flow is similar to Matlab, using the keywords `if-elseif-else-end`, and the logical operators `||` and `&&` for *or* and *and* respectively.

i = 10
for i in 0:3:15
    if i < 5
        println("$(i) is less than 5")
    elseif i < 10
        println("$(i) is less than 10")
    else
        if i == 10
            println("the value is 10")
        else
            println("$(i) is bigger than 10")
        end
    end
end

#' ### Comprehensions
#' Similar to languages like Haskell and Python, Julia supports the use of simple loops in the construction of arrays and dictionaries, called comprehenions.
#'
#' A list of increasing integers,

[i for i in 1:5]

#' Matrices can be built by including multiple indices,

[i*j for i in 1:5, j in 5:10]

#' Conditional statements can be used to filter out some values,

[i for i in 1:10 if i%2 == 1]

#' A similar syntax can be used for building dictionaries

Dict("$i" => i for i in 1:10 if i%2 == 1)

#' ### Functions
#' A simple function is defined as follows,

function print_hello()
    println("hello")
end
print_hello()

#' Arguments can be added to a function,

function print_it(x)
    println(x)
end
print_it("hello")
print_it(1.234)
print_it(:my_id)

#' Optional keyword arguments are also possible

function print_it(x; prefix="value:")
    println("$(prefix) $x")
end
print_it(1.234)
print_it(1.234, prefix="val:")

#' The keyword `return` is used to specify the return values of a function.

function mult(x; y=2.0)
    return x * y
end
mult(4.0)

#+

mult(4.0, y=5.0)

#' ### Other notes on types
#' Usually, specifing types is not required to use Julia.  However, it can be helpful to understand the basics of Julia types for debugging.
#' For example this list has a type of `Array{Int64,1}` indicating that it is a one dimensional array of integer values.

[1, 5, -2, 7]

#' In this example, the decimal values lead to a one dimensional array of floating point values, i.e. `Array{Float64,1}`.  Notice that the integer `7` is promoted to a `Float64`, because all elements in the array need share a common type.

[1.0, 5.2, -2.1, 7]

#' ### Mutable vs immutable objects
#' Some types in Julia are *mutable*, which means you can change the values inside them. A good example is an array. You can modify the contents of an array without having to make a new array.
#' In contrast, types like `Float64` are *immutable*. You can't modify the contents of a `Float64`.
#' This is something to be aware of when passing types into functions. For example:

function mutability_example(mutable_type::Vector{Int}, immutable_type::Int)
    mutable_type[1] += 1
    immutable_type += 1
    return
end

mutable_type = [1, 2, 3]
immutable_type = 1

mutability_example(mutable_type, immutable_type)

println("mutable_type: $(mutable_type)")
println("immutable_type: $(immutable_type)")

#' Because `Vector{Int}` is a mutable type, modifying the variable inside the function changed the value outside of the function. In constrast, the change to `immutable_type` didn't modify the value outside the function.
#' You can check mutability with the `isimmutable` function.

@show isimmutable([1, 2, 3])
@show isimmutable(1);

#' ### Using Packages and the Package Manager
#' No matter how wonderful Julia's base language is, at some point you will want to use an extension package.  Some of these are built-in, for example random number generation is available in the `Random` package in the standard library. These packages are loaded with the commands `using` and `import`.

using Random
Random.seed!(33);

[rand() for i in 1:10]

#' The Package Manager is used to install packages that are not part of Julia's standard library.
#' For example the following can be used to install JuMP,
#+ tangle = false

using Pkg
Pkg.add("JuMP")

#' For a complete list of registed Julia packages see the package listing at https://pkg.julialang.org/.
#' From time to you may wish to use a Julia package that is not registered.  In this case a git repository URL can be used to install the package.
#+ tangle = false

using Pkg
Pkg.add("https://github.com/user-name/MyPackage.jl.git")

#' Note that for clarity this example uses the package manager `Pkg`.  Julia includes an interactive package manager that can be accessed using `]`.  [This video](https://youtu.be/76KL8aSz0Sg) gives an overview of using the interactive package manager environment.
#' The state of installed packages can also be saved in two files: `Project.toml` and `Manifest.toml`. If these files are stored in the same directory than a notebook, the state of the packages can be recovered by running
#+ tangle = false
import Pkg
Pkg.activate(@__DIR__)
Pkg.instantiate()

#' ### HELP!
#' Julia includes a help mode that can be accessed using `?`.  Entering any object (e.g. function, type, struct, ...) into the help mode will show its documentation, if any is available.

#' ### Some Common Gotchas

#' #### MethodError
#' A common error in Julia is `MethodError`, which indicates that the function is not defined for the given value.  For example, by default the `ceil` function is not defined for complex numbers.  The "closest candidates" list suggest some Julia types that the function is defined for.
#+ tangle = false

ceil(1.2 + 2.3im)