code: def site_percolation(ij, occupied_sites):
      r
      from collections import namedtuple
      Np = sp.size(occupied_sites)
      occupied_bonds = sp.all(occupied_sites[ij], axis=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]
      clusters[~occupied_sites] = -1
      s_labels = ps.tools.make_contiguous(clusters + 1)
      if sp.any(~occupied_sites):
          s_labels -= 1
      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 site percolation
      process given a list of occupied sites.
  
      Parameters
      ----------
      ij : array_like
          An N x 2 array of [site_A, site_B] connections.  If two connected
          sites are both occupied they are part of the same cluster, as it
          the bond connecting them.
  
      occupied_sites : boolean
          A list indicating whether sites are occupied or not
  
      Returns
      -------
      A tuple containing a list of site and bond labels, indicating which
      cluster each belongs to.  A value of -1 indicates unoccupied.
  
      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.