docs: def arcte_with_lazy_pagerank(adjacency_matrix, rho, epsilon, number_of_threads=None):
      """
      Extracts local community features for all graph nodes based on the partitioning of node-centric similarity vectors.
  
      Inputs:  - A in R^(nxn): Adjacency matrix of an undirected network represented as a SciPy Sparse COOrdinate matrix.
               - rho: Restart probability
               - epsilon: Approximation threshold
  
      Outputs: - X in R^(nxC_n): The latent space embedding represented as a SciPy Sparse COOrdinate matrix.
      """
      adjacency_matrix = sparse.csr_matrix(adjacency_matrix)
      number_of_nodes = adjacency_matrix.shape[0]
  
      if number_of_threads is None:
          number_of_threads = get_threads_number()
      if number_of_threads == 1:
          # Calculate natural random walk transition probability matrix.
          rw_transition, out_degree, in_degree = get_natural_random_walk_matrix(adjacency_matrix, make_shared=False)
  
          a = adjacency_matrix.copy()
          a.data = np.ones_like(a.data)
          edge_count_vector = np.squeeze(np.asarray(a.sum(axis=0), dtype=np.int64))
  
          iterate_nodes = np.where(edge_count_vector != 0)[0]
          argsort_indices = np.argsort(edge_count_vector[iterate_nodes])
          iterate_nodes = iterate_nodes[argsort_indices][::-1]
          iterate_nodes = iterate_nodes[np.where(edge_count_vector[iterate_nodes] > 1.0)[0]]
  
          # iterate_nodes = np.where(out_degree != 0)[0]
          # argsort_indices = np.argsort(out_degree[iterate_nodes])
          # iterate_nodes = iterate_nodes[argsort_indices][::-1]
          # iterate_nodes = iterate_nodes[np.where(out_degree[iterate_nodes] > 1.0)[0]]
  
          local_features = arcte_with_lazy_pagerank_worker(iterate_nodes,
                                                           rw_transition.indices,
                                                           rw_transition.indptr,
                                                           rw_transition.data,
                                                           out_degree,
                                                           in_degree,
                                                           rho,
                                                           epsilon)
      else:
          # Calculate natural random walk transition probability matrix.
          rw_transition, out_degree, in_degree = get_natural_random_walk_matrix(adjacency_matrix, make_shared=True)
  
          a = adjacency_matrix.copy()
          a.data = np.ones_like(a.data)
          edge_count_vector = np.squeeze(np.asarray(a.sum(axis=0), dtype=np.int64))
  
          iterate_nodes = np.where(edge_count_vector != 0)[0]
          argsort_indices = np.argsort(edge_count_vector[iterate_nodes])
          iterate_nodes = iterate_nodes[argsort_indices][::-1]
          iterate_nodes = iterate_nodes[np.where(edge_count_vector[iterate_nodes] > 1.0)[0]]
  
          # iterate_nodes = np.where(out_degree != 0)[0]
          # argsort_indices = np.argsort(out_degree[iterate_nodes])
          # iterate_nodes = iterate_nodes[argsort_indices][::-1]
          # iterate_nodes = iterate_nodes[np.where(out_degree[iterate_nodes] > 1.0)[0]]
  
          pool = mp.Pool(number_of_threads)
          node_chunks = list(parallel_chunks(iterate_nodes, number_of_threads))
          node_count = 0
          for chunk in node_chunks:
              node_count += len(list(chunk))
          results = list()
          for chunk_no in range(len(pool._pool)):
              pool.apply_async(arcte_with_lazy_pagerank_worker,
                               args=(node_chunks[chunk_no],
                                     rw_transition.indices,
                                     rw_transition.indptr,
                                     rw_transition.data,
                                     out_degree,
                                     in_degree,
                                     rho,
                                     epsilon),
                               callback=results.append)
          pool.close()
          pool.join()
          # local_features = sparse.hstack(results)
          local_features = results[0]
          for additive_features in results[1:]:
              local_features += additive_features
          local_features = sparse.csr_matrix(local_features)
  
      # Form base community feature matrix.
      identity_matrix = sparse.csr_matrix(sparse.eye(number_of_nodes, number_of_nodes, dtype=np.float64))
      adjacency_matrix_ones = adjacency_matrix
      adjacency_matrix_ones.data = np.ones_like(adjacency_matrix.data)
      base_community_features = identity_matrix + adjacency_matrix_ones
  
      # Stack horizontally matrices to form feature matrix.
      try:
          features = sparse.hstack([base_community_features, local_features]).tocsr()
      except ValueError as e:
          print("Failure with horizontal feature stacking.")
          features = base_community_features
  
      return features