How to use pylops - 10 common examples

def test_ISTA_FISTA(par):
    """Invert problem with ISTA/FISTA
    """
    np.random.seed(42)
    Aop = MatrixMult(np.random.randn(par['ny'], par['nx']))

    x = np.zeros(par['nx'])
    x[par['nx'] // 2] = 1
    x[3] = 1
    x[par['nx'] - 4] = -1
    y = Aop * x

    eps = 0.5
    perc = 30
    maxit = 2000

    # ISTA with too high alpha (check that exception is raised)
    with pytest.raises(ValueError):
        xinv, _, _ = ISTA(Aop, y, maxit, eps=eps, alpha=1e5, monitorres=True,
                          tol=0, returninfo=True)

np.abs(taxis[1] - taxis[1]))
            Pop = FFT2D(dims=dims, nffts=nffts,
                        sampling=sampling)
            Pop = Pop.H
        SIop = Rop * Pop
    elif 'radon' in kind:
        prec = True
        dotcflag = 0
        kindradon = kind.split('-')[-1]
        if ndims == 3:
            Pop = Radon3D(taxis, spataxis, spat1axis, paxis, p1axis,
                          centeredh=centeredh, kind=kindradon, engine=engine)
            dimsp = (paxis.size, p1axis.size, taxis.size)

        else:
            Pop = Radon2D(taxis, spataxis, paxis,
                          centeredh=centeredh, kind=kindradon, engine=engine)
            dimsp = (paxis.size, taxis.size)
        SIop = Rop * Pop
    elif kind == 'sliding':
        prec = True
        dotcflag = 0
        if ndims == 3:
            dspat = np.abs(spataxis[1]-spataxis[0])
            dspat1 = np.abs(spat1axis[1]-spat1axis[0])
            nspat, nspat1 = spataxis.size, spat1axis.size
            spataxis_local = np.linspace(-dspat * nwin[0] // 2,
                                         dspat * nwin[0] // 2,
                                         nwin[0])
            spat1axis_local = np.linspace(-dspat1 * nwin[1] // 2,
                                          dspat1 * nwin[1] // 2,
                                          nwin[1])

# Blurring guassian operator
nh = [15, 25]
hz = np.exp(-0.1*np.linspace(-(nh[0]//2), nh[0]//2, nh[0])**2)
hx = np.exp(-0.03*np.linspace(-(nh[1]//2), nh[1]//2, nh[1])**2)
hz /= np.trapz(hz) # normalize the integral to 1
hx /= np.trapz(hx) # normalize the integral to 1
h = hz[:, np.newaxis] * hx[np.newaxis, :]

fig, ax = plt.subplots(1, 1, figsize=(5, 3))
him = ax.imshow(h)
ax.set_title('Blurring operator')
fig.colorbar(him, ax=ax)
ax.axis('tight')

Cop = pylops.signalprocessing.Convolve2D(Nz * Nx, h=h,
                                         offset=(nh[0] // 2,
                                                 nh[1] // 2),
                                         dims=(Nz, Nx), dtype='float32')

###############################################################################
# We first apply the blurring operator to the sharp image. We then
# try to recover the sharp input image by inverting the convolution operator
# from the blurred image. Note that when we perform inversion without any
# regularization, the deblurred image will show some ringing due to the
# 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,

% linearization)

        G = [np.hstack((np.diag(G1[itheta] * np.ones(nt0)),
                        np.diag(G2[itheta] * np.ones(nt0)),
                        np.diag(G3[itheta] * np.ones(nt0))))
             for itheta in range(ntheta)]
        G = np.vstack(G).reshape(ntheta * nt0, 3 * nt0)

        # Create wavelet operator
        C = convmtx(wav, nt0)[:, len(wav) // 2:-len(wav) // 2 + 1]
        C = [C] * ntheta
        C = block_diag(*C)

        # Combine operators
        M = np.dot(C, np.dot(G, D))
        return MatrixMult(M, dims=spatdims)

    else:
        # Create wavelet operator
        Cop = Convolve1D(np.prod(np.array(dims)), h=wav,
                         offset=len(wav)//2, dir=0, dims=dims)

        # create AVO operator
        AVOop = AVOLinearModelling(theta, vsvp, spatdims=spatdims,
                                   linearization=linearization)

        # Create derivative operator
        dimsm = list(dims)
        dimsm[1] = AVOop.npars
        Dop = FirstDerivative(np.prod(np.array(dimsm)), dims=dimsm,
                              dir=0, sampling=1.)
        return Cop*AVOop*Dop

Cop = pylops.signalprocessing.Convolve2D(Nz * Nx, h=h,
                                         offset=(nh[0] // 2,
                                                 nh[1] // 2),
                                         dims=(Nz, Nx), dtype='float32')

###############################################################################
# We first apply the blurring operator to the sharp image. We then
# try to recover the sharp input image by inverting the convolution operator
# from the blurred image. Note that when we perform inversion without any
# regularization, the deblurred image will show some ringing due to the
# 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,

np.diag(G3[itheta] * np.ones(nt0))))
             for itheta in range(ntheta)]
        G = np.vstack(G).reshape(ntheta * nt0, 3 * nt0)

        # Create wavelet operator
        C = convmtx(wav, nt0)[:, len(wav) // 2:-len(wav) // 2 + 1]
        C = [C] * ntheta
        C = block_diag(*C)

        # Combine operators
        M = np.dot(C, np.dot(G, D))
        return MatrixMult(M, dims=spatdims)

    else:
        # Create wavelet operator
        Cop = Convolve1D(np.prod(np.array(dims)), h=wav,
                         offset=len(wav)//2, dir=0, dims=dims)

        # create AVO operator
        AVOop = AVOLinearModelling(theta, vsvp, spatdims=spatdims,
                                   linearization=linearization)

        # Create derivative operator
        dimsm = list(dims)
        dimsm[1] = AVOop.npars
        Dop = FirstDerivative(np.prod(np.array(dimsm)), dims=dimsm,
                              dir=0, sampling=1.)
        return Cop*AVOop*Dop

def test_PoststackLinearModelling2d(par):
    """Dot-test and inversion for PoststackLinearModelling in 2d
    """

    # Dense
    PPop_dense = PoststackLinearModelling(wav, nt0=nz, spatdims=nx,
                                          explicit=True)
    assert dottest(PPop_dense, nz * nx, nz * nx, tol=1e-4)

    # Linear operator
    PPop = PoststackLinearModelling(wav, nt0=nz, spatdims=nx,
                                    explicit=False)
    assert dottest(PPop, nz * nx, nz * nx, tol=1e-4)

    # Compare data
    d = (PPop * m2d.flatten()).reshape(nz, nx)
    d_dense = (PPop_dense * m2d.flatten()).reshape(nz, nx)
    assert_array_almost_equal(d, d_dense, decimal=4)

    # Inversion
    for explicit in [True, False]:
        if explicit and not par['simultaneous'] and par['epsR'] is None:
            dict_inv = {}
        elif explicit and not par['simultaneous'] and par['epsR'] is not None:
            dict_inv = dict(damp=0 if par['epsI'] is None else par['epsI'],
                            iter_lim=10)
        else:
            dict_inv = dict(damp=0 if par['epsI'] is None else par['epsI'],
                            iter_lim=10)

np.ones(par['nx']) + \
             par['imag'] * np.outer(np.ones(par['ny']),
                                    np.arange(par['nx']))[:, :, np.newaxis] * \
             np.ones(par['nx'])
    x['2'] = np.outer(np.ones(par['ny']),
                      np.ones(par['nx']))[:, :, np.newaxis] * \
             np.arange(par['nx']) + \
             par['imag'] * np.outer(np.ones(par['ny']),
                                    np.ones(par['nx']))[:, :, np.newaxis] * \
             np.arange(par['nx'])

    for dir in [0, 1, 2]:
        Fop = Flip(par['ny']*par['nx']*par['nx'],
                   dims=(par['ny'], par['nx'], par['nx']),
                   dir=dir, dtype=par['dtype'])
        assert dottest(Fop, par['ny']*par['nx']*par['nx'],
                       par['ny']*par['nx']*par['nx'])

        y = Fop * x[str(dir)].flatten()
        xadj = Fop.H * y.flatten()
        xadj = xadj.reshape(par['ny'], par['nx'], par['nx'])
        assert_array_equal(x[str(dir)], xadj)

complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(Hop, Hop * x, damp=1e-20, iter_lim=300, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=4)

    # use numpy matrix directly in the definition of the operator
    H1op = HStack([G1, MatrixMult(G2, dtype=par['dtype'])],
                  dtype=par['dtype'])
    assert dottest(H1op, par['ny'], 2 * par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

    # use scipy matrix directly in the definition of the operator
    G1 = sp_random(par['ny'], par['nx'], density=0.4).astype('float32')
    H2op = HStack([G1, MatrixMult(G2, dtype=par['dtype'])],
                  dtype=par['dtype'])
    assert dottest(H2op, par['ny'], 2 * par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

def test_BlockDiag(par):
    """Dot-test and inversion for BlockDiag operator
    """
    np.random.seed(0)
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    x = np.ones(2*par['nx']) + par['imag']*np.ones(2*par['nx'])

    BDop = BlockDiag([MatrixMult(G1, dtype=par['dtype']),
                      MatrixMult(G2, dtype=par['dtype'])],
                     dtype=par['dtype'])
    assert dottest(BDop, 2*par['ny'], 2*par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

    xlsqr = lsqr(BDop, BDop * x, damp=1e-20, iter_lim=500, show=0)[0]
    assert_array_almost_equal(x, xlsqr, decimal=3)

    # use numpy matrix directly in the definition of the operator
    BD1op = BlockDiag([MatrixMult(G1, dtype=par['dtype']), G2],
                      dtype=par['dtype'])
    assert dottest(BD1op, 2 * par['ny'], 2 * par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

    # use scipy matrix directly in the definition of the operator
    G2 = sp_random(par['ny'], par['nx'], density=0.4).astype('float32')
    BD2op = BlockDiag([MatrixMult(G1, dtype=par['dtype']), G2],
                      dtype=par['dtype'])
    assert dottest(BD2op, 2 * par['ny'], 2 * par['nx'],