code: def dgap_l21l1(
      M,
      G,
      Z,
      active_set,
      alpha_space,
      alpha_time,
      phi,
      phiT,
      n_orient,
      highest_d_obj,
      w_space=None,
      w_time=None,
  ):
      """Duality gap for the time-frequency mixed norm inverse problem.
  
      See :footcite:`GramfortEtAl2012,NdiayeEtAl2016`
  
      Parameters
      ----------
      M : array, shape (n_sensors, n_times)
          The data.
      G : array, shape (n_sensors, n_sources)
          Gain matrix a.k.a. lead field.
      Z : array, shape (n_active, n_coefs)
          Sources in TF domain.
      active_set : array of bool, shape (n_sources, )
          Mask of active sources.
      alpha_space : float
          The spatial regularization parameter.
      alpha_time : float
          The temporal regularization parameter. The higher it is the smoother
          will be the estimated time series.
      phi : instance of _Phi
          The TF operator.
      phiT : instance of _PhiT
          The transpose of the TF operator.
      n_orient : int
          Number of dipoles per locations (typically 1 or 3).
      highest_d_obj : float
          The highest value of the dual objective so far.
      w_space : array, shape (n_positions, )
          Array of spatial weights.
      w_time : array, shape (n_positions, n_coefs)
          Array of TF weights.
  
      Returns
      -------
      gap : float
          Dual gap
      p_obj : float
          Primal objective
      d_obj : float
          Dual objective. gap = p_obj - d_obj
      R : array, shape (n_sensors, n_times)
          Current residual (M - G * X)
  
      References
      ----------
      .. footbibliography::
      """
      X = phiT(Z)
      GX = np.dot(G[:, active_set], X)
      R = M - GX
  
      # some functions need w_time only on active_set, other need it completely
      if w_time is not None:
          w_time_as = w_time[active_set[::n_orient]]
      else:
          w_time_as = None
      if w_space is not None:
          w_space_as = w_space[active_set[::n_orient]]
      else:
          w_space_as = None
  
      penaltyl1 = norm_l1_tf(Z, phi, n_orient, w_time_as)
      penaltyl21 = norm_l21_tf(Z, phi, n_orient, w_space_as)
      nR2 = sum_squared(R)
      p_obj = (
          0.5 * nR2 + alpha_space * penaltyl21 + alpha_time * penaltyl1
      )
  
      l1_ratio = alpha_time / (alpha_space + alpha_time)
      dual_norm = norm_epsilon_inf(
          G, R, phi, l1_ratio, n_orient, w_space=w_space, w_time=w_time
      )
      scaling = min(1.0, (alpha_space + alpha_time) / dual_norm)
  
      d_obj = (scaling - 0.5 * (scaling**2)) * nR2 + scaling * np.sum(
          R * GX
      )
      d_obj = max(d_obj, highest_d_obj)
  
      gap = p_obj - d_obj
      return gap, p_obj, d_obj, R
  
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 19 for calling a global function with example usage: dgap_l21l1(M, G, Z, active_set, alpha_space, alpha_time, phi, phiT, n_orient, highest_d_obj, w_space, w_time) and returns: gap,  p_obj,  d_obj,  R