How to use the pylops.optimization.sparsity.FISTA function in pylops

assert_array_almost_equal(x, xinv, decimal=1)

        # FISTA
        xinv, _, _ = FISTA(Aop, y, maxit, eps=eps, threshkind=threshkind,
                           tol=0, returninfo=True, show=False)
        assert_array_almost_equal(x, xinv, decimal=1)

    # Percentile based ISTA and FISTA
    for threshkind in ['hard-percentile', 'soft-percentile', 'half-percentile']:
        # ISTA
        xinv, _, _ = ISTA(Aop, y, maxit, perc=perc, threshkind=threshkind,
                          tol=0, returninfo=True, show=False)
        assert_array_almost_equal(x, xinv, decimal=1)

        # FISTA
        xinv, _, _ = FISTA(Aop, y, maxit, perc=perc, threshkind=threshkind,
                           tol=0, returninfo=True, show=False)
        assert_array_almost_equal(x, xinv, decimal=1)

def test_ISTA_FISTA_missing_perc():
    """Check error is raised if perc=None and threshkind is percentile based
    """
    with pytest.raises(ValueError):
        _ = ISTA(Identity(5), np.ones(5), 10, perc=None,
                 threshkind='soft-percentile')
    with pytest.raises(ValueError):
        _ = FISTA(Identity(5), np.ones(5), 10, perc=None,
                  threshkind='soft-percentile')

lsm = pylops.waveeqprocessing.LSM(z, x, t, sources, recs, v0, wav, wavc,
                                  mode='analytic')

d = lsm.Demop * refl.ravel()
d = d.reshape(ns, nr, nt)

madj = lsm.Demop.H * d.ravel()
madj = madj.reshape(nx, nz)

minv = lsm.solve(d.ravel(), solver=lsqr, **dict(iter_lim=100))
minv = minv.reshape(nx, nz)

minv_sparse = lsm.solve(d.ravel(),
                        solver=pylops.optimization.sparsity.FISTA,
                        **dict(eps=1e2, niter=100))
minv_sparse = minv_sparse.reshape(nx, nz)

# demigration
dadj = lsm.Demop * madj.ravel()
dadj = dadj.reshape(ns, nr, nt)

dinv = lsm.Demop * minv.ravel()
dinv = dinv.reshape(ns, nr, nt)

dinv_sparse = lsm.Demop * minv_sparse.ravel()
dinv_sparse = dinv_sparse.reshape(ns, nr, nt)

# sphinx_gallery_thumbnail_number = 2
fig, axs = plt.subplots(2, 2, figsize=(10, 8))
axs[0][0].imshow(refl.T, cmap='gray', vmin=-1, vmax=1)

# radon operator
npx = 61
pxmax = 5e-4
px = np.linspace(-pxmax, pxmax, npx)

Rop = pylops.signalprocessing.Radon2D(taxis, xaxis, px, kind='linear',
                                      interp='nearest',  centeredh=False,
                                      dtype='float64')

# adjoint Radon transform
xadj = Rop.H * y.flatten()
xadj = xadj.reshape(npx, par['nt'])

# sparse Radon transform
xinv, niter, cost = \
    pylops.optimization.sparsity.FISTA(Rop, y.flatten(), 15,
                                       eps=1e1, returninfo=True)
xinv = xinv.reshape(npx, par['nt'])

# filtering
xfilt = np.zeros_like(xadj)
xfilt[npx//2-3:npx//2+4] = xadj[npx//2-3:npx//2+4]

yfilt = Rop * xfilt.flatten()
yfilt = yfilt.reshape(par['nx'], par['nt'])

# filtering on sparse transform
xinvfilt = np.zeros_like(xinv)
xinvfilt[npx//2-3:npx//2+4] = xinv[npx//2-3:npx//2+4]

yinvfilt = Rop * xinvfilt.flatten()
yinvfilt = yinvfilt.reshape(par['nx'], par['nt'])

#   .. math::
#        J_1 = ||\mathbf{y} - \mathbf{R} \mathbf{F} \mathbf{p}||_2^2 +
#              \epsilon ||\mathbf{p}||_1
#
# where :math:`\mathbf{F}` is the FFT operator. We will thus use the
# :py:class:`pylops.optimization.sparsity.ISTA` and
# :py:class:`pylops.optimization.sparsity.FISTA` solvers to estimate our input
# signal.

pista, niteri, costi = \
    pylops.optimization.sparsity.ISTA(Rop*FFTop.H, y, niter=1000,
                                      eps=0.1, tol=1e-7, returninfo=True)
xista = FFTop.H*pista

pfista, niterf, costf = \
    pylops.optimization.sparsity.FISTA(Rop*FFTop.H, y, niter=1000,
                                       eps=0.1, tol=1e-7, returninfo=True)
xfista = FFTop.H*pfista

fig, axs = plt.subplots(2, 1, figsize=(12, 8))
fig.suptitle('Data reconstruction with sparsity', fontsize=14,
             fontweight='bold', y=0.9)
axs[0].plot(f, np.abs(X), 'k', lw=3)
axs[0].plot(f, np.abs(pista), '--r', lw=3)
axs[0].plot(f, np.abs(pfista), '--g', lw=3)
axs[0].set_xlim(0, 30)
axs[0].set_title('Frequency domain')
axs[1].plot(t[iava], y, '.k', ms=20, label='available samples')
axs[1].plot(t, x, 'k', lw=3, label='original')
axs[1].plot(t, xista, '--r', lw=3, label='ISTA')
axs[1].plot(t, xfista, '--g', lw=3, label='FISTA')
axs[1].set_title('Time domain')

# noisy
xls = \
    pylops.optimization.leastsquares.RegularizedInversion(Cop, [], yn,
                                                          returninfo=False,
                                                          **dict(damp=1e-1,
                                                                 atol=1e-3,
                                                                 iter_lim=100,
                                                                 show=0))

xista, niteri, costi = \
    pylops.optimization.sparsity.ISTA(Cop, yn, niter=100, eps=5e-1,
                                      tol=1e-5, returninfo=True)

xfista, niterf, costf = \
    pylops.optimization.sparsity.FISTA(Cop, yn, niter=100, eps=5e-1,
                                       tol=1e-5, returninfo=True)

fig, ax = plt.subplots(1, 1, figsize=(8, 3))
ax.plot(t, x, 'k', lw=8, label=r'$x$')
ax.plot(t, y, 'r', lw=4, label=r'$y=Ax$')
ax.plot(t, yn, '--b', lw=4, label=r'$y_n$')
ax.plot(t, xls, '--g', lw=4, label=r'$x_{LS}$')
ax.plot(t, xista, '--m', lw=4, label=r'$x_{ISTA} (niter=%d)$' % niteri)
ax.plot(t, xfista, '--y', lw=4, label=r'$x_{FISTA} (niter=%d)$' % niterf)
ax.set_title('Noisy deconvolution', fontsize=14, fontweight='bold')
ax.legend()
plt.tight_layout()

fig, ax = plt.subplots(1, 1, figsize=(8, 3))
ax.semilogy(costi, 'm', lw=2, label=r'$x_{ISTA} (niter=%d)$' % niteri)
ax.semilogy(costf, 'y', lw=2, label=r'$x_{FISTA} (niter=%d)$' % niterf)

# instabilities of the inverse process. Using a L1 solver with a DWT
# preconditioner or TV regularization allows to recover sharper contrasts.
imblur = Cop * im.flatten()

imdeblur = \
    pylops.optimization.leastsquares.NormalEquationsInversion(Cop, None,
                                                              imblur,
                                                              maxiter=50)

Wop = pylops.signalprocessing.DWT2D((Nz, Nx), wavelet='haar', level=3)
Dop = [pylops.FirstDerivative(Nz * Nx, dims=(Nz, Nx), dir=0, edge=False),
       pylops.FirstDerivative(Nz * Nx, dims=(Nz, Nx), dir=1, edge=False)]
DWop = Dop + [Wop, ]

imdeblurfista = \
    pylops.optimization.sparsity.FISTA(Cop * Wop.H, imblur, eps=1e-1,
                                       niter=100)[0]
imdeblurfista = Wop.H * imdeblurfista

imdeblurtv = \
    pylops.optimization.sparsity.SplitBregman(Cop, Dop, imblur.flatten(),
                                              niter_outer=10, niter_inner=5,
                                              mu=1.5, epsRL1s=[2e0, 2e0],
                                              tol=1e-4, tau=1., show=False,
                                              ** dict(iter_lim=5, damp=1e-4))[0]

imdeblurtv1 = \
    pylops.optimization.sparsity.SplitBregman(Cop, DWop,
                                              imblur.flatten(),
                                              niter_outer=10, niter_inner=5,
                                              mu=1.5, epsRL1s=[1e0, 1e0, 1e0],
                                              tol=1e-4, tau=1., show=False,