rcm.jl

function sort_inverse(p)
    out = [ i for i in 1:length(p) ]
    for i = 1:length(p)
        out[p[i]] = i
    end
    out
end

# TODO: limit allocations to improve performance
# Reverse Cuthill McKee
function rcm(A)
    n = size(A, 1)
    e = zeros(n)
    
    visited = falses(n)
    touched = falses(n)
    
    # find out which is our first node
    degrees = sum(A .!= 0, 1)[:]
    sorted = sortperm(degrees)
    sorted_inv = sort_inverse(sorted)
    
    # prepare the queue
    queue = [0 for i in 1:n]
    queue_start = 1
    queue_end = 1
    
    while queue_start < n
        # enqueue x₀ 
        x₀ = sorted[sorted_inv[!visited]][1]
        touched[x₀] = true
        queue[queue_end] = x₀
        queue_end += 1
        
        while queue_start < queue_end
            x = queue[queue_start]
            visited[x] = true
            queue_start += 1
            
            e[x] = 1
            adj = (A * e) .!= 0
            e[x] = 0
            
            adj_idx = find(adj & !touched)
            if length(adj_idx) == 0
                continue
            end
            
            sort!(adj_idx, by = adj -> sorted_inv[adj])
            touched[adj_idx] = true

            queue[queue_end : queue_end + length(adj_idx) - 1] = adj_idx
            queue_end += length(adj_idx)
        end
    end
    
    reverse(queue)
end