def bond_percolation(ij, occupied_bonds):
    r"""
    Calculates the site and bond occupancy status for a bond percolation
    process given a list of occupied bonds.

    Parameters
    ----------
    ij : array_like
        An N x 2 array of [site_A, site_B] connections.  A site is
        considered occupied if any of it's connecting bonds are occupied.

    occupied_bonds: boolean
        A list indicating whether a bond is occupied or not

    Returns
    -------
    A tuple contain a list of site and bond labels, indicating which
    cluster each belongs to.  A value of -1 indicates uninvaded.

    Notes
    -----
    The ``connected_components`` function of scipy.csgraph will give ALL
    sites a cluster number whether they are occupied or not, so this
    function essentially adjusts the cluster numbers to represent a
    percolation process.

    """
    from collections import namedtuple
    Np = sp.amax(ij) + 1
    adj_mat = sprs.csr_matrix((occupied_bonds, (ij[:, 0], ij[:, 1])),
                              shape=(Np, Np))
    adj_mat.eliminate_zeros()
    clusters = csgraph.connected_components(csgraph=adj_mat, directed=False)[1]
    valid_clusters = sp.bincount(clusters) > 1
    mapping = -sp.ones(shape=(clusters.max()+1, ), dtype=int)
    mapping[valid_clusters] = sp.arange(0, valid_clusters.sum())
    s_labels = mapping[clusters]
    b_labels = sp.amin(s_labels[ij], axis=1)
    tup = namedtuple('cluster_labels', ('sites', 'bonds'))
    return tup(s_labels, b_labels)