Implementation of the Cuthill-McKee algorithm.

    Permute a symmetric binary matrix into a band matrix form with a small bandwidth.

    Parameters
    ----------
    matrix : ndarray, dtype=bool, shape = [n, n]
        The matrix is internally converted to a symmetric matrix by setting each element [i,j] to True if either
        [i,j] or [j,i] evaluates to true.

    Returns
    -------
    order : list of int
        Permutation intices

    Examples
    --------
    >>> A = np.array([[0,0,1,1], [0,0,0,0], [1,0,1,0], [1,0,0,0]])
    >>> p = cuthill_mckee(A)
    >>> A
    array([[0, 0, 1, 1],
           [0, 0, 0, 0],
           [1, 0, 1, 0],
           [1, 0, 0, 0]])
    >>> A[p,:][:,p]
    array([[0, 0, 0, 0],
           [0, 0, 1, 0],
           [0, 1, 0, 1],
           [0, 0, 1, 1]])

def cuthill_mckee(matrix):
    """Implementation of the Cuthill-McKee algorithm.

    Permute a symmetric binary matrix into a band matrix form with a small bandwidth.

    Parameters
    ----------
    matrix : ndarray, dtype=bool, shape = [n, n]
        The matrix is internally converted to a symmetric matrix by setting each element [i,j] to True if either
        [i,j] or [j,i] evaluates to true.

    Returns
    -------
    order : list of int
        Permutation intices

    Examples
    --------
    >>> A = np.array([[0,0,1,1], [0,0,0,0], [1,0,1,0], [1,0,0,0]])
    >>> p = cuthill_mckee(A)
    >>> A
    array([[0, 0, 1, 1],
           [0, 0, 0, 0],
           [1, 0, 1, 0],
           [1, 0, 0, 0]])
    >>> A[p,:][:,p]
    array([[0, 0, 0, 0],
           [0, 0, 1, 0],
           [0, 1, 0, 1],
           [0, 0, 1, 1]])
    """
    matrix = np.atleast_2d(matrix)
    n, m = matrix.shape
    assert(n == m)

    # make sure the matrix is really symmetric. This is equivalent to
    # converting a directed adjacency matrix into a undirected adjacency matrix.
    matrix = np.logical_or(matrix, matrix.T)

    degree = np.sum(matrix, 0)
    order = [np.argmin(degree)]

    for i in range(n):
        adj = np.nonzero(matrix[order[i]])[0]
        adj = [a for a in adj if a not in order]
        if not adj:
            idx = [i for i in range(n) if i not in order]
            order.append(idx[np.argmin(degree[idx])])
        else:
            if len(adj) == 1:
                order.append(adj[0])
            else:
                adj = np.asarray(adj)
                i = adj[np.argsort(degree[adj])]
                order.extend(i.tolist())
        if len(order) == n:
            break

    return order