code: def bond_percolation(ij, occupied_bonds):
      r
      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)
docs: 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.