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
run of the RSSMP algorithm