docs: def _single_multichannel_mp_run(X, Phi, bound, selection_rule, stop_crit,
                                  max_iter, verbose=False, pad=0,
                                  random_state=None, memory=Memory(None)):
      """ run of the structured variant of the RSSMP algorithm """
      rng = check_random_state(random_state)
  
      # padding as v stak
      pad = int(pad)
      n_channels = X.shape[0]
      X = np.hstack((np.zeros((n_channels, pad)), X,
                     np.zeros((n_channels, pad))))
      n_samples = X.shape[1]
      n_projs = Phi.doth(X).shape[1]
      err_mse = {}
  
      # Initialisation
      residual = np.hstack((X.copy(), np.zeros((n_channels,
                                                max(Phi.sizes) / 2))))
  
      s_rep = np.zeros((n_channels, n_projs))
      X_est = np.zeros((n_channels, n_samples))
      # Main algorithm
      coeffs = np.zeros((n_channels, n_projs))
  
      it_number = 0
      current_lambda = 1
      for c_idx in range(n_channels):
          err_mse[c_idx] = []
          err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))
  
      # Decomposition loop: stopping criteria is either SNR or iteration number
      while (current_lambda > bound) & (it_number < max_iter):
  
          # pick a shift at random : in each size
          rndshifts = {}
          for c_idx in range(n_channels):
              rndshifts[c_idx] = []
          for s_idx, L in enumerate(Phi.sizes):
              shift = rng.randint(low=0, high=L / 4)
              for c_idx in range(n_channels):
                  coeffs[c_idx, s_idx * n_samples:(s_idx + 1) * n_samples] = \
                      mdct(residual[c_idx, shift:shift + n_samples], L).ravel()
                  rndshifts[c_idx].append(shift)
  
          # Multichannel mode : we combine projections
          combined = selection_rule(coeffs ** 2)
  
          # Select a new element
          idx = np.argmax(np.abs(combined))
          # find scale and frequency bin of selected atom
          s_idx = idx // n_samples
          L = Phi.sizes[s_idx]
          F = n_samples // (L // 2)
          frame = (idx - (s_idx * n_samples)) % F
          freq_bin = ((idx - (s_idx * n_samples))) // F
  
          mdct_wf = memory.cache(mdct_waveform)
  
          # Update coefficients and residual
          current_lambda_array = np.zeros(n_channels)
          for c_idx in range(n_channels):
              s_rep[c_idx, idx] += coeffs[c_idx, idx]
  
              # Only one method now : local update via a cached waveform
              pos = (frame * L / 2) - L / 4 + rndshifts[c_idx][s_idx]
              residual[c_idx, pos:pos + L] -= coeffs[c_idx, idx] * \
                  mdct_wf(L, freq_bin)
  
              # also add it to the reconstruction
              X_est[c_idx, pos:pos + L] += coeffs[c_idx, idx] * \
                  mdct_wf(L, freq_bin)
  
              # error computation (err_mse)
              err_mse[c_idx].append(linalg.norm(residual[c_idx, :]))
  
              current_lambda_array[c_idx] = np.sqrt(
                  1. - err_mse[c_idx][-1] / err_mse[c_idx][-2])
  
          current_lambda = stop_crit(current_lambda_array)
  
          if verbose:
              print("Iteration %d : Current lambda of %1.4f" % (
                    it_number, current_lambda))
          it_number += 1
  
      return X_est[:, pad: -pad], err_mse