code: def _gamma_map_opt(
      M,
      G,
      alpha,
      maxit=10000,
      tol=1e-6,
      update_mode=1,
      group_size=1,
      gammas=None,
      verbose=None,
  ):
      """Hierarchical Bayes (Gamma-MAP).
  
      Parameters
      ----------
      M : array, shape=(n_sensors, n_times)
          Observation.
      G : array, shape=(n_sensors, n_sources)
          Forward operator.
      alpha : float
          Regularization parameter (noise variance).
      maxit : int
          Maximum number of iterations.
      tol : float
          Tolerance parameter for convergence.
      group_size : int
          Number of consecutive sources which use the same gamma.
      update_mode : int
          Update mode, 1: MacKay update (default), 3: Modified MacKay update.
      gammas : array, shape=(n_sources,)
          Initial values for posterior variances (gammas). If None, a
          variance of 1.0 is used.
      %(verbose)s
  
      Returns
      -------
      X : array, shape=(n_active, n_times)
          Estimated source time courses.
      active_set : array, shape=(n_active,)
          Indices of active sources.
      """
      G = G.copy()
      M = M.copy()
  
      if gammas is None:
          gammas = np.ones(G.shape[1], dtype=np.float64)
  
      eps = np.finfo(float).eps
  
      n_sources = G.shape[1]
      n_sensors, n_times = M.shape
  
      # apply normalization so the numerical values are sane
      M_normalize_constant = np.linalg.norm(np.dot(M, M.T), ord="fro")
      M /= np.sqrt(M_normalize_constant)
      alpha /= M_normalize_constant
      G_normalize_constant = np.linalg.norm(G, ord=np.inf)
      G /= G_normalize_constant
  
      if n_sources % group_size != 0:
          raise ValueError(
              "Number of sources has to be evenly dividable by the "
              "group size"
          )
  
      n_active = n_sources
      active_set = np.arange(n_sources)
  
      gammas_full_old = gammas.copy()
  
      if update_mode == 2:
          denom_fun = np.sqrt
      else:
          # do nothing
          def denom_fun(x):
              return x
  
      last_size = -1
      for itno in range(maxit):
          gammas[np.isnan(gammas)] = 0.0
  
          gidx = np.abs(gammas) > eps
          active_set = active_set[gidx]
          gammas = gammas[gidx]
  
          # update only active gammas (once set to zero it stays at zero)
          if n_active > len(active_set):
              n_active = active_set.size
              G = G[:, gidx]
  
          CM = np.dot(G * gammas[np.newaxis, :], G.T)
          CM.flat[:: n_sensors + 1] += alpha
          # Invert CM keeping symmetry
          U, S, _ = _safe_svd(CM, full_matrices=False)
          S = S[np.newaxis, :]
          del CM
          CMinv = np.dot(U / (S + eps), U.T)
          CMinvG = np.dot(CMinv, G)
          A = np.dot(
              CMinvG.T, M
          )  # mult. w. Diag(gamma) in gamma update
  
          if update_mode == 1:
              # MacKay fixed point update (10) in [1]
              numer = gammas**2 * np.mean((A * A.conj()).real, axis=1)
              denom = gammas * np.sum(G * CMinvG, axis=0)
          elif update_mode == 2:
              # modified MacKay fixed point update (11) in [1]
              numer = gammas * np.sqrt(
                  np.mean((A * A.conj()).real, axis=1)
              )
              denom = np.sum(
                  G * CMinvG, axis=0
              )  # sqrt is applied below
          else:
              raise ValueError("Invalid value for update_mode")
  
          if group_size == 1:
              if denom is None:
                  gammas = numer
              else:
                  gammas = numer / np.maximum(
                      denom_fun(denom), np.finfo("float").eps
                  )
          else:
              numer_comb = np.sum(numer.reshape(-1, group_size), axis=1)
              if denom is None:
                  gammas_comb = numer_comb
              else:
                  denom_comb = np.sum(
                      denom.reshape(-1, group_size), axis=1
                  )
                  gammas_comb = numer_comb / denom_fun(denom_comb)
  
              gammas = np.repeat(gammas_comb / group_size, group_size)
  
          # compute convergence criterion
          gammas_full = np.zeros(n_sources, dtype=np.float64)
          gammas_full[active_set] = gammas
  
          err = np.sum(np.abs(gammas_full - gammas_full_old)) / np.sum(
              np.abs(gammas_full_old)
          )
  
          gammas_full_old = gammas_full
  
          breaking = err < tol or n_active == 0
          if len(gammas) != last_size or breaking:
              logger.info(
                  "Iteration: %d\t active set size: %d\t convergence: "
                  "%0.3e" % (itno, len(gammas), err)
              )
              last_size = len(gammas)
  
          if breaking:
              break
  
      if itno < maxit - 1:
          logger.info("\nConvergence reached !\n")
      else:
          warn("\nConvergence NOT reached !\n")
  
      # undo normalization and compute final posterior mean
      n_const = np.sqrt(M_normalize_constant) / G_normalize_constant
      x_active = n_const * gammas[:, None] * A
  
      return x_active, active_set
  
imports: import numpy
  from ..fixes import _safe_svd
  from ..forward import is_fixed_orient
  from ..minimum_norm.inverse import _check_reference, _log_exp_var
  from ..utils import logger, verbose, warn
  from .mxne_inverse import _check_ori, _compute_residual, _make_dipoles_sparse, _make_sparse_stc, _prepare_gain, _reapply_source_weighting
label_desc: Use image node_id 1 for calling a global function with example usage: _gamma_map_opt(M, G, alpha, maxit, tol, update_mode, group_size, gammas, verbose) and returns: x_active,  active_set, x