docs: def nmf_ensemble(data, k, n_runs=10, W_list=[], **nmf_params):
      """
      Runs an ensemble method on the list of NMF W matrices...
  
      Args:
          data: genes x cells array (should be log + cell-normalized)
          k: number of classes
          n_runs (optional): number of random initializations of state estimation
          M_list (optional): list of M arrays from state estimation
          se_params (optional): optional poisson_estimate_state params
  
      Returns:
          W_new
          H_new
      """
      nmf = NMF(k)
      if len(W_list)==0:
          W_list = []
          for i in range(n_runs):
              W = nmf.fit_transform(data)
              W_list.append(W)
      W_stacked = np.hstack(W_list)
      nmf_w = nmf.fit_transform(W_stacked)
      nmf_h = nmf.components_
      H_new = data.T.dot(nmf_w).T
      nmf2 = NMF(k, init='custom')
      nmf_w = nmf2.fit_transform(data, W=nmf_w, H=H_new)
      H_new = nmf2.components_
      #W_new = W_new/W_new.sum(0)
      # alternatively, use nmf_w and h_new as initializations for another NMF round?
      return nmf_w, H_new