docs: def poisson_cluster(data, k, init=None, max_iters=100):
      """
      Performs Poisson hard EM on the given data.
  
      Args:
          data (array): A 2d array- genes x cells. Can be dense or sparse; for best performance, sparse matrices should be in CSC format.
          k (int): Number of clusters
          init (array, optional): Initial centers - genes x k array. Default: None, use kmeans++
          max_iters (int, optional): Maximum number of iterations. Default: 100
  
      Returns:
          a tuple of two arrays: a cells x 1 vector of cluster assignments,
          and a genes x k array of cluster means.
      """
      # TODO: be able to use a combination of fixed and unknown starting points
      # e.g., have init values only for certain genes, have a row of all
      # zeros indicating that kmeans++ should be used for that row.
      genes, cells = data.shape
      #print 'starting: ', centers
      if sparse.issparse(data) and not sparse.isspmatrix_csc(data):
          data = sparse.csc_matrix(data)
      init, assignments = kmeans_pp(data, k, centers=init)
      centers = np.copy(init)
      assignments = np.zeros(cells)
      for it in range(max_iters):
          lls = poisson_ll(data, centers)
          #cluster_dists = np.zeros((cells, k))
          new_assignments = np.argmax(lls, 1)
          if np.equal(assignments, new_assignments).all():
              #print 'ending: ', centers
              return new_assignments, centers
          for c in range(k):
              if sparse.issparse(data):
                  if data[:,new_assignments==c].shape[0]==0:
                      # re-initialize centers?
                      new_c, _ = kmeans_pp(data, k, centers[:,:c])
                      centers[:,c] = new_c[:,c]
                  else:
                      centers[:,c] = np.asarray(data[:,new_assignments==c].mean(1)).flatten()
              else:
                  if len(data[:,new_assignments==c])==0:
                      new_c, _ = kmeans_pp(data, k, centers[:,:c])
                      centers[:,c] = new_c[:,c]
                  else:
                      centers[:,c] = np.mean(data[:,new_assignments==c], 1)
          assignments = new_assignments
      return assignments, centers