docs: def mroc(adjacency_matrix, alpha):
      """
      Extracts hierarchical community features using the MROC method.
  
      Introduced in: Wang, X., Tang, L., Liu, H., & Wang, L. (2013).
                     Learning with multi-resolution overlapping communities.
                     Knowledge and information systems, 36(2), 517-535.
  
      Inputs:  - A in R^(nxn): Adjacency matrix of an undirected network represented as a SciPy Sparse COOrdinate matrix.
               - alpha: A maximum community size stopping threshold.
  
      Outputs: - X in R^(nxC_n): The latent space embedding represented as a SciPy Sparse COOrdinate matrix.
      """
      # Find number of nodes
      number_of_nodes = adjacency_matrix.shape[0]
  
      ####################################################################################################################
      # Base community calculation
      ####################################################################################################################
      # Initialize empty lists
      base_list = list()
      base_row = list()
      base_col = list()
  
      # Save function handles for speed
      append_base_list = base_list.append
      append_base_row = base_row.append
      append_base_col = base_col.append
  
      # Find base communities
      adjacency_matrix = adjacency_matrix.tocsc()
      number_of_base_communities = 0
      for i in range(number_of_nodes):
          # Calculate base community
          base_community = set(adjacency_matrix.getcol(i).indices)
          base_community.add(i)
          flag = True
          for c in base_list:
              if c == base_community:
                  flag = False
                  break
          if flag:
              append_base_list(base_community)
  
              for n in base_community:
                  append_base_row(n)
                  append_base_col(number_of_base_communities)
  
              number_of_base_communities += 1
  
      # Form sparse matrices
      base_row = np.array(base_row)
      base_col = np.array(base_col)
      base_data = np.ones(base_row.size, dtype=np.float64)
      features = sparse.coo_matrix((base_data, (base_row, base_col)),
                                   shape=(number_of_nodes, number_of_base_communities))
  
      features = features.tocsr()
      base_community_number = features.shape[1]
  
      print('Base communities calculated.')
  
      reverse_index_csr = copy.copy(features)
      reverse_index_csc = reverse_index_csr.tocsc()
      reverse_index_csr = reverse_index_csr.tocsr()
  
      reverse_index_rows = np.ndarray(number_of_nodes, dtype=np.ndarray)
      reverse_index_cols = np.ndarray(number_of_nodes, dtype=np.ndarray)
      for n in range(number_of_nodes):
          reverse_index_row = reverse_index_csr.getrow(n)
          reverse_index_rows[n] = reverse_index_row.indices
  
          if n < base_community_number:
              reverse_index_col = reverse_index_csc.getcol(n)
              reverse_index_cols[n] = reverse_index_col.indices
  
      flag = True
  
      print('Start merge iterations.')
  
      iteration = 0
  
      while flag:
          level_row = list()
          level_col = list()
  
          append_level_row = level_row.append
          append_level_col = level_col.append
  
          unavailable_communities = -1*np.ones(reverse_index_csc.shape[1])
          unavailable_communities_counter = 0
          next_level_communities = list()
          append_next_level_community = next_level_communities.append
          number_of_communities = 0
          for j in range(reverse_index_csr.shape[1]):
              if j in unavailable_communities:
                  continue
              must_break = reverse_index_csr.shape[1] - unavailable_communities_counter
              print(must_break)
              if must_break < 1:
                  break
              unavailable_communities[unavailable_communities_counter] = j
              unavailable_communities_counter += 1
              c_j = reverse_index_cols[j]
  
              indices = community_neighbors(c_j, reverse_index_rows, unavailable_communities, unavailable_communities_counter)
  
              max_similarity = -1
              community_index = 0
              for jj in indices:
                  c_jj = reverse_index_cols[jj]
                  similarity = jaccard(c_j, c_jj)
                  if similarity > max_similarity:
                      max_similarity = similarity
                      community_index = jj
  
              jj = community_index
              if max_similarity > 0:
                  # Merge two communities
                  c_jj = reverse_index_cols[jj]
                  c_new = np.union1d(c_j, c_jj)
  
                  flag_1 = np.setdiff1d(c_new, c_j)
                  flag_2 = np.setdiff1d(c_new, c_jj)
                  if (flag_1.size != 0) and (flag_2.size != 0):
                      for n in c_new:
                          append_level_row(n)
                          append_level_col(number_of_communities)
  
                      if c_new.size < alpha:
                          append_next_level_community(number_of_communities)
  
                      number_of_communities += 1
                  unavailable_communities[unavailable_communities_counter] = jj
                  unavailable_communities_counter += 1
  
          level_row = np.array(level_row)
          level_col = np.array(level_col)
          level_data = np.ones(level_row.size, dtype=np.float64)
          communities = sparse.coo_matrix((level_data, (level_row, level_col)),
                                          shape=(number_of_nodes, number_of_communities))
  
          if communities.getnnz() == 0:
              break
  
          features = sparse.hstack([features, communities])
  
          reverse_index_csc = copy.copy(communities)
          reverse_index_csc = reverse_index_csc.tocsc()
          reverse_index_csc = reverse_index_csc[:, np.array(next_level_communities)]
          reverse_index_csr = reverse_index_csc.tocsr()
  
          reverse_index_rows = np.ndarray(number_of_nodes, dtype=np.ndarray)
          reverse_index_cols = np.ndarray(len(next_level_communities), dtype=np.ndarray)
          for n in range(number_of_nodes):
              reverse_index_row = reverse_index_csr.getrow(n)
              reverse_index_rows[n] = reverse_index_row.indices
  
              if n < len(next_level_communities):
                  reverse_index_col = reverse_index_csc.getcol(n)
                  reverse_index_cols[n] = reverse_index_col.indices
  
          if len(next_level_communities) > 1:
              flag = True
  
          iteration += 1
          print('Iteration: ', iteration)
          print('List length', len(next_level_communities))
  
      return features