code: def iterative_tf_mixed_norm_solver(
      M,
      G,
      alpha_space,
      alpha_time,
      n_tfmxne_iter,
      wsize=64,
      tstep=4,
      maxit=3000,
      tol=1e-8,
      debias=True,
      n_orient=1,
      dgap_freq=10,
      verbose=None,
  ):
      """Solve TF L0.5/L1 + L0.5 inverse problem with BCD + active set approach.
  
      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_space: float
          The spatial regularization parameter. The higher it is the less there
          will be active sources.
      alpha_time : float
          The temporal regularization parameter. The higher it is the smoother
          will be the estimated time series. 0 means no temporal regularization,
          a.k.a. irMxNE.
      n_tfmxne_iter : int
          Number of TF-MxNE iterations. If > 1, iterative reweighting is applied.
      wsize : int or array-like
          Length of the STFT window in samples (must be a multiple of 4).
          If an array is passed, multiple TF dictionaries are used (each having
          its own wsize and tstep) and each entry of wsize must be a multiple
          of 4.
      tstep : int or array-like
          Step between successive windows in samples (must be a multiple of 2,
          a divider of wsize and smaller than wsize/2) (default: wsize/2).
          If an array is passed, multiple TF dictionaries are used (each having
          its own wsize and tstep), and each entry of tstep must be a multiple
          of 2 and divide the corresponding entry of wsize.
      maxit : int
          The maximum number of iterations for each TF-MxNE problem.
      tol : float
          If absolute difference between estimates at 2 successive iterations
          is lower than tol, the convergence is reached. Also used as criterion
          on duality gap for each TF-MxNE problem.
      debias : bool
          Debias source estimates.
      n_orient : int
          The number of orientation (1 : fixed or 3 : free or loose).
      dgap_freq : int or np.inf
          The duality gap is evaluated every dgap_freq iterations.
      %(verbose)s
  
      Returns
      -------
      X : array, shape (n_active, n_times)
          The source estimates.
      active_set : array
          The mask of active sources.
      E : list
          The value of the objective function over iterations.
      """
      n_sensors, n_times = M.shape
      n_sources = G.shape[1]
      n_positions = n_sources // n_orient
  
      tstep = np.atleast_1d(tstep)
      wsize = np.atleast_1d(wsize)
      if len(tstep) != len(wsize):
          raise ValueError(
              "The same number of window sizes and steps must be "
              "passed. Got tstep = %s and wsize = %s" % (tstep, wsize)
          )
  
      n_steps = np.ceil(n_times / tstep.astype(float)).astype(int)
      n_freqs = wsize // 2 + 1
      n_coefs = n_steps * n_freqs
      phi = _Phi(wsize, tstep, n_coefs, n_times)
      phiT = _PhiT(tstep, n_freqs, n_steps, n_times)
  
      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)
  
      # space and time penalties, and inverse of their derivatives:
      def g_space(Z):
          return np.sqrt(
              np.sqrt(
                  phi.norm(Z, ord=2).reshape(-1, n_orient).sum(axis=1)
              )
          )
  
      def g_space_prime_inv(Z):
          return 2.0 * g_space(Z)
  
      def g_time(Z):
          return np.sqrt(
              np.sqrt(
                  np.sum(
                      (np.abs(Z) ** 2.0).reshape(
                          (n_orient, -1), order="F"
                      ),
                      axis=0,
                  )
              ).reshape((-1, Z.shape[1]), order="F")
          )
  
      def g_time_prime_inv(Z):
          return 2.0 * g_time(Z)
  
      E = list()
  
      active_set = np.ones(n_sources, dtype=bool)
      Z = np.zeros((n_sources, phi.n_coefs.sum()), dtype=np.complex128)
  
      for k in range(n_tfmxne_iter):
          active_set_0 = active_set.copy()
          Z0 = Z.copy()
  
          if k == 0:
              w_space = None
              w_time = None
          else:
              w_space = 1.0 / g_space_prime_inv(Z)
              w_time = g_time_prime_inv(Z)
              w_time[w_time == 0.0] = -1.0
              w_time = 1.0 / w_time
              w_time[w_time < 0.0] = 0.0
  
          (
              X,
              Z,
              active_set_,
              E_,
              _,
          ) = _tf_mixed_norm_solver_bcd_active_set(
              M,
              G[:, active_set],
              alpha_space,
              alpha_time,
              lc[active_set[::n_orient]],
              phi,
              phiT,
              Z_init=Z,
              w_space=w_space,
              w_time=w_time,
              n_orient=n_orient,
              maxit=maxit,
              tol=tol,
              dgap_freq=dgap_freq,
              verbose=None,
          )
  
          active_set[active_set] = active_set_
  
          if active_set.sum() > 0:
              l21_penalty = np.sum(g_space(Z.copy()))
              l1_penalty = phi.norm(g_time(Z.copy()), ord=1).sum()
  
              p_obj = (
                  0.5
                  * np.linalg.norm(
                      M - np.dot(G[:, active_set], X), "fro"
                  )
                  ** 2.0
                  + alpha_space * l21_penalty
                  + alpha_time * l1_penalty
              )
              E.append(p_obj)
  
              logger.info(
                  "Iteration %d: active set size=%d, E=%f"
                  % (k + 1, active_set.sum() / n_orient, p_obj)
              )
  
              # Check convergence
              if np.array_equal(active_set, active_set_0):
                  max_diff = np.amax(np.abs(Z - Z0))
                  if max_diff < tol:
                      print(
                          "Convergence reached after %d reweightings!"
                          % k
                      )
                      break
          else:
              p_obj = 0.5 * np.linalg.norm(M) ** 2.0
              E.append(p_obj)
              logger.info(
                  "Iteration %d: as_size=%d, E=%f"
                  % (k + 1, active_set.sum() / n_orient, p_obj)
              )
              break
  
      if debias:
          if active_set.sum() > 0:
              bias = compute_bias(
                  M, G[:, active_set], X, n_orient=n_orient
              )
              X *= bias[:, np.newaxis]
  
      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 23 for calling a global function with example usage: iterative_tf_mixed_norm_solver(M, G, alpha_space, alpha_time, n_tfmxne_iter, wsize, tstep, maxit, tol, debias, n_orient, dgap_freq, verbose) and returns: X,  active_set,  E, np, unknown, np, unknown