@with_kw struct Cell
    boundary::BitVector
    interior::BitVector
    size::Int64 = sum(boundary .| interior)
end

"""
    separator!(g, subgraph_nodes, max_imbalance=0.6, nsample=20)
Call flowcutter (with s and t random) on g and pick nodes at random in the cut matching the 60% imbalance
heuristic to form a separator. Then split the graph g with respect to the separator.
Parameter max_imbalance define the heurestic to choose a cut amoung all generated by
flowcutter. Parameter max_nsample define the number of call to flowcutter.

# Arguments
- `g::SimpleGraph{Int64}`: the subgraph to run flowcutter.
- `subgraph_nodes::Vector{Int64}`: a table where element at index i is corresponding label in the original root graph.
- `max_imbalance::Float64` max imbalance to select cut from flowcutter result.
- `nsample::Int64` max number of try to get a separator.
- `seed::Int64` base seed for RNG choosing input for flowcutter.
# Return
- `sep::Vector{Int64}`: the set of indices of the nodes in separator.
- `split_parts::Array{Array{Int64,1},1}` several parts of the graph resulted from separation.
"""
function separator!(
    g::SimpleGraph{Int64},
    subgraph_nodes::Vector{Int64},
    max_imbalance::Float64,
    nsample::Int64,
    seed::Int64,
)::Tuple{Vector{Int64},Vector{Vector{Int}}}

    # run flowcutter many times and collect all of these cuts
    n = length(subgraph_nodes)
    local_cuts = [Cut[] for _ = 1:Threads.nthreads()]
    rngs = [MersenneTwister(seed + i) for i = 1:Threads.nthreads()]
    Threads.@threads for i = 1:nsample
        s, t = sample(rngs[Threads.threadid()], 1:n, 2, replace=false)
        @inbounds append!(local_cuts[Threads.threadid()], flowcutter(g, s, t))
    end
    cuts = collect(Iterators.flatten(local_cuts))
    filter!(c -> c.imbalance < max_imbalance, cuts)
    # remove dominated cuts
    candidates = Dict{Int64, Cut}()
    for c in cuts
        size = c.size
        # update candidates if they are not previous candidates with its size
        # or in case the candidate is better than its predecessor i.e same size but lower imbalance
        if !(size in keys(candidates)) || c.imbalance < candidates[size].imbalance
            candidates[size] = c
        end
    end

    # select cut with min expansion   
    cuts = [c for c in values(candidates)]
    cut = cuts[findmin(map(c->c.expansion, cuts))[2]]
    # select the nodes from the larger side of the cut
    sep = unique(map(a -> a.first, cut.arcs))
    
    # split the subgraph in several parts
    # be careful with index while removing vertices from sep
    labels = collect(1:n)
    @inbounds for rm_origin_id in sep
        rm = findfirst(id -> id == rm_origin_id, labels)
        rem_vertex!(g, rm)
        labels[rm] = labels[end]
        pop!(labels)
    end
    # catch the multiple parts created from the split
    split_parts = connected_components(g)
    # be careful to return indices which make sense in the original root graph
    return  @inbounds (
        map(node -> subgraph_nodes[node], sep),
        map(
            vector -> map(node -> subgraph_nodes[labels[node]], vector),
            split_parts,
        ),
    )
end



"""
    iterative_dissection(g, best_tw = typemax(Int64), max_imbalance = 0.6, nsample = 20, seed = 42,)
Computes an approximation of the upper bound of the treewidth of the graph g and
an order for elimination using the iterative dissection and the flow cutter algorithm.


# Arguments
- `g::SimpleGraph{Int64}`: the graph to analyse.
- `best_tw::Int64 = typemax(Int64)` stop when hitting a tw bigger than this.
- `max_imbalance::Float64 = 0.6` criteria for selecting cuts that will build the separator.
- `nsample::Int64 = 20` the number of calls to flowcutter with random inputs.
- `seed::Int64 = 4242` the base seed for the RNG sampling the inputs of flowcutter and the choosen one cut for separator.

# Return
- `order::Vector{Int64}` an array of vertices index.
- `treewidth::Int64` the approximation of treewidth.
"""
function iterative_dissection(
        g::SimpleGraph{Int64}, 
        best_tw::Int64 = typemax(Int64),
        max_imbalance::Float64 = 0.6,
        nsample::Int64 = 20,
        seed::Int64 = 4242,
    )::Pair{Vector{Int64}, Int64}
    n = nv(g)
    order = zeros(Int64, n)
    treewidth = 0
    q = PriorityQueue{Cell, Int64}(Base.Order.Reverse)
    enqueue!(q, Cell(boundary = falses(n), interior = trues(n)), n)
    i = n
    while (!isempty(q)) && (treewidth < best_tw)
        cell = dequeue!(q)
        subgraph_nodes = findall(it -> it > 0, cell.interior)
        n = length(subgraph_nodes)
        if cell.size == treewidth + 1
            while !isempty(q)
                @inbounds order[(i - n + 1):i] .= subgraph_nodes
                i -= n
                subgraph_nodes = findall(it -> it > 0, dequeue!(q).interior)
                n = length(subgraph_nodes)
            end
            @inbounds order[(i - n + 1):i] .= subgraph_nodes
            i -= n
            break
        end
        
        graph = induced_subgraph(g, subgraph_nodes)[1]
        # if graph is a tree or complete we can stop
        
        nedges = ne(graph)
        if nedges == n * (n - 1) / 2
            @inbounds order[(i - n + 1):i] .= subgraph_nodes
            i -= n
            treewidth = max(cell.size - 1, treewidth)
            # @debug "bagsize = $(cell.size) tw = $treewidth"
            continue
        elseif nedges == n - 1
            @inbounds order[(i - n + 1):i] .= tree_order!(graph, subgraph_nodes)
            i -= n
            treewidth = max(1 + sum(cell.boundary), treewidth)
            continue
        end
        
        # compute separator and cut graph in several parts (graph and indices)
        sep, toqueue = separator!(graph, subgraph_nodes, max_imbalance, nsample, seed)

        # update order and treewidth
        k = length(sep) 
        @inbounds order[(i - k + 1):i] = sep
        i -= k
        treewidth = max(k + sum(cell.boundary) - 1, treewidth)
        # @debug "bagsize = $(k + sum(cell.boundary)) tw = $treewidth"
        n = nv(g)
        # add next subgraphs to the queue
        @inbounds for new_interior in toqueue
            Ic = falses(n)
            Ic[new_interior] .= true
            Bc = falses(n)
            Bc[sep] .= true
            Bound_union_Sep = findall(Bc .| cell.boundary)
            Bc[sep] .= false
            tmp_Bc = Bound_union_Sep[findall(
                node -> any(
                    node2 -> has_edge(g, node, node2), 
                    new_interior
                    ), 
                Bound_union_Sep
            )]
            Bc[tmp_Bc] .= true
            new_cell = Cell(boundary = Bc, interior = Ic)
            enqueue!(q, new_cell, new_cell.size)
        end
    end
    return order => treewidth
end
