docs: def _single_mp_run(x, Phi, bound, max_iter, verbose=False, pad=0,
                     random_state=None, memory=Memory(None)):
      """ run of the RSSMP algorithm """
  
      rng = check_random_state(random_state)
      pad = int(pad)
      x = np.concatenate((np.zeros(pad), x, np.zeros(pad)))
  
      n = x.size
      m = Phi.doth(x).size
      err_mse = []
  
      # Initialisation
      residual = np.concatenate((x.copy(), np.zeros(max(Phi.sizes) // 2)))
  
      s = np.zeros(m)
      x_est = np.zeros(n)
      # Main algorithm
      coeffs = np.zeros(m)
      it_number = 0
      current_lambda = 1
      err_mse.append(linalg.norm(residual))
  
      # 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 scale_idx, size in enumerate(Phi.sizes):
              shift = rng.randint(low=0, high=size // 4)
              coeffs[scale_idx * n:(scale_idx + 1) * n] = mdct(
                  residual[shift:shift + n], size).ravel()
              rndshifts.append(shift)
  
          # Select a new element
          idx = np.argmax(np.abs(coeffs))
  
          # Update coefficients
          s[idx] += coeffs[idx]
  
          # Only one method now : local update via a cached waveform
          # find scale and frequency bin of selected atom
          mdct_wf = memory.cache(mdct_waveform)
  
          scale_idx = idx // n
          size = Phi.sizes[scale_idx]
          F = n // (size // 2)
          frame = (idx - (scale_idx * n)) % F
          freq_bin = ((idx - (scale_idx * n))) // F
          pos = (frame * size // 2) - size // 4 + rndshifts[scale_idx]
          residual[pos:pos + size] -= coeffs[idx] * mdct_wf(size, freq_bin)
  
          # also add it to the reconstruction
          x_est[pos:pos + size] += coeffs[idx] * mdct_wf(size, freq_bin)
  
          # error computation (err_mse)
          err_mse.append(linalg.norm(residual))
  
          current_lambda = np.sqrt(1 - err_mse[-1] / err_mse[-2])
          if current_lambda <= bound:
              x_est[pos:pos + size] -= coeffs[idx] * mdct_wf(size, freq_bin)
          if verbose:
              print("Iteration %d : Current lambda of %1.4f" % (
                    it_number, current_lambda))
          it_number += 1
  
      return x_est, err_mse