code: def _cholesky_omp(
      X, y, n_nonzero_coefs, tol=None, copy_X=True, return_path=False
  ):
      """Orthogonal Matching Pursuit step using the Cholesky decomposition.
  
      Parameters
      ----------
      X : ndarray of shape (n_samples, n_features)
          Input dictionary. Columns are assumed to have unit norm.
  
      y : ndarray of shape (n_samples,)
          Input targets.
  
      n_nonzero_coefs : int
          Targeted number of non-zero elements.
  
      tol : float, default=None
          Targeted squared error, if not None overrides n_nonzero_coefs.
  
      copy_X : bool, default=True
          Whether the design matrix X must be copied by the algorithm. A false
          value is only helpful if X is already Fortran-ordered, otherwise a
          copy is made anyway.
  
      return_path : bool, default=False
          Whether to return every value of the nonzero coefficients along the
          forward path. Useful for cross-validation.
  
      Returns
      -------
      gamma : ndarray of shape (n_nonzero_coefs,)
          Non-zero elements of the solution.
  
      idx : ndarray of shape (n_nonzero_coefs,)
          Indices of the positions of the elements in gamma within the solution
          vector.
  
      coef : ndarray of shape (n_features, n_nonzero_coefs)
          The first k values of column k correspond to the coefficient value
          for the active features at that step. The lower left triangle contains
          garbage. Only returned if ``return_path=True``.
  
      n_active : int
          Number of active features at convergence.
      """
      if copy_X:
          X = X.copy("F")
      else:  # even if we are allowed to overwrite, still copy it if bad order
          X = np.asfortranarray(X)
  
      min_float = np.finfo(X.dtype).eps
      nrm2, swap = linalg.get_blas_funcs(("nrm2", "swap"), (X,))
      (potrs,) = get_lapack_funcs(("potrs",), (X,))
  
      alpha = np.dot(X.T, y)
      residual = y
      gamma = np.empty(0)
      n_active = 0
      indices = np.arange(X.shape[1])  # keeping track of swapping
  
      max_features = X.shape[1] if tol is not None else n_nonzero_coefs
  
      L = np.empty((max_features, max_features), dtype=X.dtype)
  
      if return_path:
          coefs = np.empty_like(L)
  
      while True:
          lam = np.argmax(np.abs(np.dot(X.T, residual)))
          if lam < n_active or alpha[lam] ** 2 < min_float:
              # atom already selected or inner product too small
              warnings.warn(premature, RuntimeWarning, stacklevel=2)
              break
  
          if n_active > 0:
              # Updates the Cholesky decomposition of X' X
              L[n_active, :n_active] = np.dot(
                  X[:, :n_active].T, X[:, lam]
              )
              linalg.solve_triangular(
                  L[:n_active, :n_active],
                  L[n_active, :n_active],
                  trans=0,
                  lower=1,
                  overwrite_b=True,
                  check_finite=False,
              )
              v = nrm2(L[n_active, :n_active]) ** 2
              Lkk = linalg.norm(X[:, lam]) ** 2 - v
              if Lkk <= min_float:  # selected atoms are dependent
                  warnings.warn(premature, RuntimeWarning, stacklevel=2)
                  break
              L[n_active, n_active] = sqrt(Lkk)
          else:
              L[0, 0] = linalg.norm(X[:, lam])
  
          X.T[n_active], X.T[lam] = swap(X.T[n_active], X.T[lam])
          alpha[n_active], alpha[lam] = alpha[lam], alpha[n_active]
          indices[n_active], indices[lam] = (
              indices[lam],
              indices[n_active],
          )
          n_active += 1
  
          # solves LL'x = X'y as a composition of two triangular systems
          gamma, _ = potrs(
              L[:n_active, :n_active],
              alpha[:n_active],
              lower=True,
              overwrite_b=False,
          )
  
          if return_path:
              coefs[:n_active, n_active - 1] = gamma
          residual = y - np.dot(X[:, :n_active], gamma)
          if tol is not None and nrm2(residual) ** 2 <= tol:
              break
          elif n_active == max_features:
              break
  
      if return_path:
          return (
              gamma,
              indices[:n_active],
              coefs[:, :n_active],
              n_active,
          )
      else:
          return gamma, indices[:n_active], n_active
  
imports: import warnings
  from math import sqrt
  from numbers import Integral, Real
  import numpy
  from scipy import linalg
  from scipy.linalg.lapack import get_lapack_funcs
  from ..base import MultiOutputMixin, RegressorMixin, _fit_context
  from ..model_selection import check_cv
  from ..utils import Bunch, as_float_array, check_array
  from ..utils._param_validation import Interval, StrOptions, validate_params
  from ..utils.metadata_routing import MetadataRouter, MethodMapping, _raise_for_params, _routing_enabled, process_routing
  from ..utils.parallel import Parallel, delayed
  from ._base import LinearModel, _pre_fit
label_desc: Use image node_id 1 for calling a global function with example usage: _cholesky_omp(X, y, n_nonzero_coefs, tol, copy_X, return_path) and returns: gamma,  indices,  coefs,  n_active, gamma,  indices,  n_active