code: def mixed_norm_solver(
      M,
      G,
      alpha,
      maxit=3000,
      tol=1e-8,
      verbose=None,
      active_set_size=50,
      debias=True,
      n_orient=1,
      solver="auto",
      return_gap=False,
      dgap_freq=10,
      active_set_init=None,
      X_init=None,
  ):
      """Solve L1/L2 mixed-norm inverse problem with active set strategy.
  
      See references :footcite:`GramfortEtAl2012,StrohmeierEtAl2016,
      BertrandEtAl2020`.
  
      Parameters
      ----------
      M : array, shape (n_sensors, n_times)
          The data.
      G : array, shape (n_sensors, n_dipoles)
          The gain matrix a.k.a. lead field.
      alpha : float
          The regularization parameter. It should be between 0 and 100.
          A value of 100 will lead to an empty active set (no active source).
      maxit : int
          The number of iterations.
      tol : float
          Tolerance on dual gap for convergence checking.
      %(verbose)s
      active_set_size : int
          Size of active set increase at each iteration.
      debias : bool
          Debias source estimates.
      n_orient : int
          The number of orientation (1 : fixed or 3 : free or loose).
      solver : 'cd' | 'bcd' | 'auto'
          The algorithm to use for the optimization. Block Coordinate Descent
          (BCD) uses Anderson acceleration for faster convergence.
      return_gap : bool
          Return final duality gap.
      dgap_freq : int
          The duality gap is computed every dgap_freq iterations of the solver on
          the active set.
      active_set_init : array, shape (n_dipoles,) or None
          The initial active set (boolean array) used at the first iteration.
          If None, the usual active set strategy is applied.
      X_init : array, shape (n_dipoles, n_times) or None
          The initial weight matrix used for warm starting the solver. If None,
          the weights are initialized at zero.
  
      Returns
      -------
      X : array, shape (n_active, n_times)
          The source estimates.
      active_set : array, shape (new_active_set_size,)
          The mask of active sources. Note that new_active_set_size is the size
          of the active set after convergence of the solver.
      E : list
          The value of the objective function over the iterations.
      gap : float
          Final duality gap. Returned only if return_gap is True.
  
      References
      ----------
      .. footbibliography::
      """
      n_dipoles = G.shape[1]
      n_positions = n_dipoles // n_orient
      _, n_times = M.shape
      alpha_max = norm_l2inf(np.dot(G.T, M), n_orient, copy=False)
      logger.info("-- ALPHA MAX : %s" % alpha_max)
      alpha = float(alpha)
      X = np.zeros((n_dipoles, n_times), dtype=G.dtype)
  
      has_sklearn = True
      try:
          from sklearn.linear_model import MultiTaskLasso  # noqa: F401
      except ImportError:
          has_sklearn = False
  
      _validate_type(solver, str, "solver")
      _check_option("solver", solver, ("cd", "bcd", "auto"))
      if solver == "auto":
          if has_sklearn and (n_orient == 1):
              solver = "cd"
          else:
              solver = "bcd"
  
      if solver == "cd":
          if n_orient == 1 and not has_sklearn:
              warn(
                  "Scikit-learn >= 0.12 cannot be found. Using block coordinate"
                  " descent instead of coordinate descent."
              )
              solver = "bcd"
          if n_orient > 1:
              warn(
                  "Coordinate descent is only available for fixed orientation. "
                  "Using block coordinate descent instead of coordinate "
                  "descent"
              )
              solver = "bcd"
  
      if solver == "cd":
          logger.info("Using coordinate descent")
          l21_solver = _mixed_norm_solver_cd
          lc = None
      else:
          assert solver == "bcd"
          logger.info("Using block coordinate descent")
          l21_solver = _mixed_norm_solver_bcd
          G = np.asfortranarray(G)
          if n_orient == 1:
              lc = np.sum(G * G, axis=0)
          else:
              lc = np.empty(n_positions)
              for j in range(n_positions):
                  G_tmp = G[:, (j * n_orient) : ((j + 1) * n_orient)]
                  lc[j] = np.linalg.norm(np.dot(G_tmp.T, G_tmp), ord=2)
  
      if active_set_size is not None:
          E = list()
          highest_d_obj = -np.inf
          if X_init is not None and X_init.shape != (
              n_dipoles,
              n_times,
          ):
              raise ValueError(
                  "Wrong dim for initialized coefficients."
              )
          active_set = (
              active_set_init
              if active_set_init is not None
              else np.zeros(n_dipoles, dtype=bool)
          )
          idx_large_corr = np.argsort(
              groups_norm2(np.dot(G.T, M), n_orient)
          )
          new_active_idx = idx_large_corr[-active_set_size:]
          if n_orient > 1:
              new_active_idx = (
                  n_orient * new_active_idx[:, None]
                  + np.arange(n_orient)[None, :]
              ).ravel()
          active_set[new_active_idx] = True
          as_size = np.sum(active_set)
          gap = np.inf
          for k in range(maxit):
              if solver == "bcd":
                  lc_tmp = lc[active_set[::n_orient]]
              elif solver == "cd":
                  lc_tmp = None
              else:
                  lc_tmp = (
                      1.01
                      * np.linalg.norm(G[:, active_set], ord=2) ** 2
                  )
              X, as_, _ = l21_solver(
                  M,
                  G[:, active_set],
                  alpha,
                  lc_tmp,
                  maxit=maxit,
                  tol=tol,
                  init=X_init,
                  n_orient=n_orient,
                  dgap_freq=dgap_freq,
              )
              active_set[active_set] = as_.copy()
              idx_old_active_set = np.where(active_set)[0]
  
              _, p_obj, d_obj, R = dgap_l21(
                  M, G, X, active_set, alpha, n_orient
              )
              highest_d_obj = max(d_obj, highest_d_obj)
              gap = p_obj - highest_d_obj
              E.append(p_obj)
              logger.info(
                  "Iteration %d :: p_obj %f :: dgap %f :: "
                  "n_active_start %d :: n_active_end %d"
                  % (
                      k + 1,
                      p_obj,
                      gap,
                      as_size // n_orient,
                      np.sum(active_set) // n_orient,
                  )
              )
              if gap < tol:
                  logger.info(
                      "Convergence reached ! (gap: %s < %s)"
                      % (gap, tol)
                  )
                  break
  
              # add sources if not last iteration
              if k < (maxit - 1):
                  idx_large_corr = np.argsort(
                      groups_norm2(np.dot(G.T, R), n_orient)
                  )
                  new_active_idx = idx_large_corr[-active_set_size:]
                  if n_orient > 1:
                      new_active_idx = (
                          n_orient * new_active_idx[:, None]
                          + np.arange(n_orient)[None, :]
                      )
                      new_active_idx = new_active_idx.ravel()
                  active_set[new_active_idx] = True
                  idx_active_set = np.where(active_set)[0]
                  as_size = np.sum(active_set)
                  X_init = np.zeros((as_size, n_times), dtype=X.dtype)
                  idx = np.searchsorted(
                      idx_active_set, idx_old_active_set
                  )
                  X_init[idx] = X
          else:
              warn("Did NOT converge ! (gap: %s > %s)" % (gap, tol))
      else:
          X, active_set, E = l21_solver(
              M,
              G,
              alpha,
              lc,
              maxit=maxit,
              tol=tol,
              n_orient=n_orient,
              init=None,
          )
          if return_gap:
              gap = dgap_l21(M, G, X, active_set, alpha, n_orient)[0]
  
      if np.any(active_set) and debias:
          bias = compute_bias(M, G[:, active_set], X, n_orient=n_orient)
          X *= bias[:, np.newaxis]
  
      logger.info(
          "Final active set size: %s" % (np.sum(active_set) // n_orient)
      )
  
      if return_gap:
          return X, active_set, E, gap
      else:
          return X, active_set, E
  
imports: import functools
  from math import sqrt
  import numpy
  from ..time_frequency._stft import istft, stft, stft_norm1, stft_norm2
  from ..utils import _check_option, _get_blas_funcs, _validate_type, logger, sum_squared, verbose, warn
  from .mxne_debiasing import compute_bias
label_desc: Use image node_id 10 for calling a global function with example usage: mixed_norm_solver(M, G, alpha, maxit, tol, verbose, active_set_size, debias, n_orient, solver, return_gap, dgap_freq, active_set_init, X_init) and returns: X,  active_set,  E,  gap, X,  active_set,  E