code: def _compute_mxne_sure(
      M,
      gain,
      alpha_grid,
      sigma,
      n_mxne_iter,
      maxit,
      tol,
      n_orient,
      active_set_size,
      debias,
      solver,
      dgap_freq,
      random_state,
      verbose,
  ):
      """Stein Unbiased Risk Estimator (SURE).
  
      Implements the finite-difference Monte-Carlo approximation
      of the SURE for Multi-Task LASSO.
  
      See reference :footcite:`DeledalleEtAl2014`.
  
      Parameters
      ----------
      M : array, shape (n_sensors, n_times)
          The data.
      gain : array, shape (n_sensors, n_dipoles)
          The gain matrix a.k.a. lead field.
      alpha_grid : array, shape (n_alphas,)
          The grid of alphas used to evaluate the SURE.
      sigma : float
          The true or estimated noise level in the data. Usually 1 if the data
          has been previously whitened using MNE whitener.
      n_mxne_iter : int
          The number of MxNE iterations. If > 1, iterative reweighting is
          applied.
      maxit : int
          Maximum number of iterations.
      tol : float
          Tolerance parameter.
      n_orient : int
          The number of orientation (1 : fixed or 3 : free or loose).
      active_set_size : int
          Size of active set increase at each iteration.
      debias : bool
          Debias source estimates.
      solver : 'cd' | 'bcd' | 'auto'
          The algorithm to use for the optimization.
      dgap_freq : int or np.inf
          The duality gap is evaluated every dgap_freq iterations.
      random_state : int | None
          The random state used in a random number generator for delta and
          epsilon used for the SURE computation.
  
      Returns
      -------
      X : array, shape (n_active, n_times)
          Coefficient matrix.
      active_set : array, shape (n_dipoles,)
          Array of indices of non-zero coefficients.
      best_alpha_ : float
          Alpha that minimizes the SURE.
  
      References
      ----------
      .. footbibliography::
      """
  
      def g(w):
          return np.sqrt(np.sqrt(groups_norm2(w.copy(), n_orient)))
  
      def gprime(w):
          return 2.0 * np.repeat(g(w), n_orient).ravel()
  
      def _run_solver(
          alpha, M, n_mxne_iter, as_init=None, X_init=None, w_init=None
      ):
          if n_mxne_iter == 1:
              X, active_set, _ = mixed_norm_solver(
                  M,
                  gain,
                  alpha,
                  maxit=maxit,
                  tol=tol,
                  active_set_size=active_set_size,
                  n_orient=n_orient,
                  debias=debias,
                  solver=solver,
                  dgap_freq=dgap_freq,
                  active_set_init=as_init,
                  X_init=X_init,
                  verbose=False,
              )
          else:
              X, active_set, _ = iterative_mixed_norm_solver(
                  M,
                  gain,
                  alpha,
                  n_mxne_iter,
                  maxit=maxit,
                  tol=tol,
                  n_orient=n_orient,
                  active_set_size=active_set_size,
                  debias=debias,
                  solver=solver,
                  dgap_freq=dgap_freq,
                  weight_init=w_init,
                  verbose=False,
              )
          return X, active_set
  
      def _fit_on_grid(gain, M, eps, delta):
          coefs_grid_1_0 = np.zeros(
              (len(alpha_grid), gain.shape[1], M.shape[1])
          )
          coefs_grid_2_0 = np.zeros(
              (len(alpha_grid), gain.shape[1], M.shape[1])
          )
          active_sets, active_sets_eps = [], []
          M_eps = M + eps * delta
          # warm start - first iteration (leverages convexity)
          logger.info("Warm starting...")
          for j, alpha in enumerate(alpha_grid):
              logger.info("alpha: %s" % alpha)
              X, a_set = _run_solver(alpha, M, 1)
              X_eps, a_set_eps = _run_solver(alpha, M_eps, 1)
              coefs_grid_1_0[j][a_set, :] = X
              coefs_grid_2_0[j][a_set_eps, :] = X_eps
              active_sets.append(a_set)
              active_sets_eps.append(a_set_eps)
          # next iterations
          if n_mxne_iter == 1:
              return coefs_grid_1_0, coefs_grid_2_0, active_sets
          else:
              coefs_grid_1 = coefs_grid_1_0.copy()
              coefs_grid_2 = coefs_grid_2_0.copy()
              logger.info("Fitting SURE on grid.")
              for j, alpha in enumerate(alpha_grid):
                  logger.info("alpha: %s" % alpha)
                  if active_sets[j].sum() > 0:
                      w = gprime(coefs_grid_1[j])
                      X, a_set = _run_solver(
                          alpha, M, n_mxne_iter - 1, w_init=w
                      )
                      coefs_grid_1[j][a_set, :] = X
                      active_sets[j] = a_set
                  if active_sets_eps[j].sum() > 0:
                      w_eps = gprime(coefs_grid_2[j])
                      X_eps, a_set_eps = _run_solver(
                          alpha, M_eps, n_mxne_iter - 1, w_init=w_eps
                      )
                      coefs_grid_2[j][a_set_eps, :] = X_eps
                      active_sets_eps[j] = a_set_eps
  
              return coefs_grid_1, coefs_grid_2, active_sets
  
      def _compute_sure_val(coef1, coef2, gain, M, sigma, delta, eps):
          n_sensors, n_times = gain.shape[0], M.shape[1]
          dof = (gain @ (coef2 - coef1) * delta).sum() / eps
          df_term = np.linalg.norm(M - gain @ coef1) ** 2
          sure = df_term - n_sensors * n_times * sigma**2
          sure += 2 * dof * sigma**2
          return sure
  
      sure_path = np.empty(len(alpha_grid))
  
      rng = check_random_state(random_state)
      # See Deledalle et al. 20214 Sec. 5.1
      eps = 2 * sigma / (M.shape[0] ** 0.3)
      delta = rng.randn(*M.shape)
  
      coefs_grid_1, coefs_grid_2, active_sets = _fit_on_grid(
          gain, M, eps, delta
      )
  
      logger.info("Computing SURE values on grid.")
      for i, (coef1, coef2) in enumerate(
          zip(coefs_grid_1, coefs_grid_2)
      ):
          sure_path[i] = _compute_sure_val(
              coef1, coef2, gain, M, sigma, delta, eps
          )
          if verbose:
              logger.info(
                  "alpha %s :: sure %s" % (alpha_grid[i], sure_path[i])
              )
      best_alpha_ = alpha_grid[np.argmin(sure_path)]
  
      X = coefs_grid_1[np.argmin(sure_path)]
      active_set = active_sets[np.argmin(sure_path)]
  
      X = X[active_set, :]
  
      return X, active_set, best_alpha_
  
imports: import numpy
  from .._fiff.proj import deactivate_proj
  from ..dipole import Dipole
  from ..fixes import _safe_svd
  from ..forward import is_fixed_orient
  from ..minimum_norm.inverse import _check_reference, _log_exp_var, _prepare_forward, combine_xyz
  from ..source_estimate import SourceEstimate, _BaseSourceEstimate, _make_stc
  from ..utils import _check_depth, _check_option, _validate_type, check_random_state, logger, sum_squared, verbose, warn
  from .mxne_optim import _Phi, groups_norm2, iterative_mixed_norm_solver, iterative_tf_mixed_norm_solver, mixed_norm_solver, norm_epsilon_inf, norm_l2inf, tf_mixed_norm_solver
label_desc: Use image node_id 13 for calling a global function with example usage: _compute_mxne_sure(M, gain, alpha_grid, sigma, n_mxne_iter, maxit, tol, n_orient, active_set_size, debias, solver, dgap_freq, random_state, verbose) and returns: X,  active_set,  best_alpha_, np, unknown, X,  active_set, sure, coefs_grid_1_0,  coefs_grid_2_0,  active_sets, coefs_grid_1,  coefs_grid_2,  active_sets