docs: def fit(self, X, y=None):
          """Fits the GraphLasso covariance model to X.
  
          Closely follows sklearn.covariance.graph_lasso.GraphLassoCV.
  
          Parameters
          ----------
          X : ndarray, shape (n_samples, n_features)
              Data from which to compute the covariance estimate
          """
          # quic-specific outputs
          self.opt_ = None
          self.cputime_ = None
          self.iters_ = None
          self.duality_gap_ = None
  
          # these must be updated upon self.fit()
          self.sample_covariance_ = None
          self.lam_scale_ = None
          self.is_fitted_ = False
  
          # initialize
          X = check_array(X, ensure_min_features=2, estimator=self)
          X = as_float_array(X, copy=False, force_all_finite=False)
  
          if self.cv is None:
              cv = (3, 10)
          elif isinstance(self.cv, int):
              cv = (self.cv, 10)  # upgrade with default number of trials
          elif isinstance(self.cv, tuple):
              cv = self.cv
  
          cv = RepeatedKFold(n_splits=cv[0], n_repeats=cv[1])
  
          self.init_coefs(X)
  
          # get path
          if isinstance(self.lams, int):
              n_refinements = self.n_refinements
              lam_1 = self.lam_scale_
              lam_0 = 1e-2 * lam_1
              path = np.logspace(np.log10(lam_0), np.log10(lam_1), self.lams)[::-1]
          else:
              path = self.lams
              n_refinements = 1
  
          # run this thing a bunch
          results = list()
          t0 = time.time()
          for rr in range(n_refinements):
              if self.sc is None:
                  # parallel version
                  this_result = Parallel(
                      n_jobs=self.n_jobs, verbose=self.verbose, backend=self.backend
                  )(
                      delayed(_quic_path)(
                          X[train],
                          path,
                          X_test=X[test],
                          lam=self.lam,
                          tol=self.tol,
                          max_iter=self.max_iter,
                          Theta0=self.Theta0,
                          Sigma0=self.Sigma0,
                          method=self.method,
                          verbose=self.verbose,
                          score_metric=self.score_metric,
                          init_method=self.init_method,
                      )
                      for train, test in cv.split(X)
                  )
              else:
                  # parallel via spark
                  train_test_grid = [(train, test) for (train, test) in cv.split(X)]
                  indexed_param_grid = list(
                      zip(range(len(train_test_grid)), train_test_grid)
                  )
                  par_param_grid = self.sc.parallelize(indexed_param_grid)
                  X_bc = self.sc.broadcast(X)
  
                  # wrap function parameters so we dont pick whole self object
                  quic_path = partial(
                      _quic_path,
                      path=path,
                      lam=self.lam,
                      tol=self.tol,
                      max_iter=self.max_iter,
                      Theta0=self.Theta0,
                      Sigma0=self.Sigma0,
                      method=self.method,
                      verbose=self.verbose,
                      score_metric=self.score_metric,
                      init_method=self.init_method,
                  )
  
                  indexed_results = dict(
                      par_param_grid.map(
                          partial(_quic_path_spark, quic_path=quic_path, X_bc=X_bc)
                      ).collect()
                  )
                  this_result = [
                      indexed_results[idx] for idx in range(len(train_test_grid))
                  ]
                  X_bc.unpersist()
  
              # Little dance to transform the list in what we need
              covs, _, scores = zip(*this_result)
              covs = zip(*covs)
              scores = zip(*scores)
              results.extend(zip(path, scores, covs))
              results = sorted(results, key=operator.itemgetter(0), reverse=True)
  
              # Find the maximum (avoid using built in 'max' function to
              # have a fully-reproducible selection of the smallest alpha
              # in case of equality)
              best_score = -np.inf
              last_finite_idx = 0
              best_index = 0
              for index, (lam, scores, _) in enumerate(results):
                  # sometimes we get -np.inf in the result (in kl-loss)
                  scores = [s for s in scores if not np.isinf(s)]
                  if len(scores) == 0:
                      this_score = -np.inf
                  else:
                      this_score = np.mean(scores)
  
                  if this_score >= .1 / np.finfo(np.float64).eps:
                      this_score = np.nan
  
                  if np.isfinite(this_score):
                      last_finite_idx = index
  
                  if this_score >= best_score:
                      best_score = this_score
                      best_index = index
  
              # Refine the grid
              if best_index == 0:
                  # We do not need to go back: we have chosen
                  # the highest value of lambda for which there are
                  # non-zero coefficients
                  lam_1 = results[0][0]
                  lam_0 = results[1][0]
  
              elif best_index == last_finite_idx and not best_index == len(results) - 1:
                  # We have non-converged models on the upper bound of the
                  # grid, we need to refine the grid there
                  lam_1 = results[best_index][0]
                  lam_0 = results[best_index + 1][0]
  
              elif best_index == len(results) - 1:
                  lam_1 = results[best_index][0]
                  lam_0 = 0.01 * results[best_index][0]
  
              else:
                  lam_1 = results[best_index - 1][0]
                  lam_0 = results[best_index + 1][0]
  
              if isinstance(self.lams, int):
                  path = np.logspace(np.log10(lam_1), np.log10(lam_0), self.lams + 2)
                  path = path[1:-1]
  
              if self.verbose and n_refinements > 1:
                  print(
                      "[GraphLassoCV] Done refinement % 2i out of %i: % 3is"
                      % (rr + 1, n_refinements, time.time() - t0)
                  )
  
          results = list(zip(*results))
          grid_scores_ = list(results[1])
          lams = list(results[0])
  
          # Finally, compute the score with lambda = 0
          lams.append(0)
          grid_scores_.append(
              cross_val_score(EmpiricalCovariance(), X, cv=cv, n_jobs=self.n_jobs)
          )
          self.grid_scores_ = np.array(grid_scores_)
          self.lam_ = self.lam * lams[best_index]
          self.cv_lams_ = [self.lam * l for l in lams]
  
          # Finally fit the model with the selected lambda
          if self.method == "quic":
              (
                  self.precision_,
                  self.covariance_,
                  self.opt_,
                  self.cputime_,
                  self.iters_,
                  self.duality_gap_,
              ) = quic(
                  self.sample_covariance_,
                  self.lam_,
                  mode="default",
                  tol=self.tol,
                  max_iter=self.max_iter,
                  Theta0=self.Theta0,
                  Sigma0=self.Sigma0,
                  path=None,
                  msg=self.verbose,
              )
          else:
              raise NotImplementedError("Only method='quic' has been implemented.")
  
          self.is_fitted_ = True
          return self