How to use the pylops.utils.dottest function in pylops

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'],

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]:
        Rop = Roll(par['ny'] * par['nx'] * par['nx'],
                   dims=(par['ny'], par['nx'], par['nx']),
                   dir=dir, shift=3, dtype=par['dtype'])
        y = Rop * x[str(dir)].flatten()
        assert dottest(Rop, par['ny'] * par['nx'] * par['nx'],
                       par['ny'] * par['nx'] * par['nx'])

        xinv = Rop.H * y
        assert_array_almost_equal(x[str(dir)].ravel(), xinv, decimal=3)

def test_DWT_2dsignal(par):
    """Dot-test and inversion for DWT operator for 2d signal
    """
    for dir in [0, 1]:
        DWTop = DWT(dims=(par['nt'], par['nx']),
                    dir=dir, wavelet='haar', level=3)
        x = np.random.normal(0., 1., (par['nt'], par['nx'])) + \
            par['imag'] * np.random.normal(0., 1., (par['nt'], par['nx']))

        assert dottest(DWTop, DWTop.shape[0], DWTop.shape[1],
                       complexflag=0 if par['imag'] == 0 else 3)

        y = DWTop * x.ravel()
        xadj = DWTop.H * y  # adjoint is same as inverse for dwt
        xinv = lsqr(DWTop, y, damp=1e-10, iter_lim=10, show=0)[0]

        assert_array_almost_equal(x.ravel(), xadj, decimal=8)
        assert_array_almost_equal(x.ravel(), xinv, decimal=8)

def test_Flip2D(par):
    """Dot-test, forward and adjoint for Flip operator on 2d signal
    """
    np.random.seed(10)
    x = {}
    x['0'] = np.outer(np.arange(par['ny']), np.ones(par['nx'])) + \
             par['imag'] * np.outer(np.arange(par['ny']), np.ones(par['nx']))
    x['1'] = np.outer(np.ones(par['ny']), np.arange(par['nx'])) + \
             par['imag'] * np.outer(np.ones(par['ny']), np.arange(par['nx']))

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

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

def test_Sum2D(par):
    """Dot-test for Sum operator on 2d signal
    """
    for dir in [0, 1]:
        dim_d = [par['ny'], par['nx']]
        dim_d.pop(dir)
        Sop = Sum(dims=(par['ny'], par['nx']),
                  dir=dir,  dtype=par['dtype'])
        assert dottest(Sop, np.prod(dim_d), par['ny'] * par['nx'])

par['imag'] * np.random.normal(0, 1, (par['nx'], par['nt']))

    # fixed indeces
    iava = np.vstack((np.arange(0, 10),
                      np.arange(0, 10)))
    Iop = Bilinear(iava, dims=(par['nx'], par['nt']), dtype=par['dtype'])
    assert dottest(Iop, 10, par['nx'] * par['nt'],
                   complexflag=0 if par['imag'] == 0 else 3)

    # decimal indeces
    Nsub = int(np.round(par['nx'] * par['nt'] * perc_subsampling))
    iavadec = np.vstack((np.random.uniform(0, par['nx'] - 1, Nsub),
                      np.random.uniform(0, par['nt'] - 1, Nsub)))
    Idecop = Bilinear(iavadec, dims=(par['nx'], par['nt']),
                      dtype=par['dtype'])
    assert dottest(Idecop, Nsub, par['nx'] * par['nt'],
                   complexflag=0 if par['imag'] == 0 else 3)

    # repeated indeces
    with pytest.raises(ValueError):
        iava_rep = iava.copy()
        iava_rep[-2] = [0, 0]
        iava_rep[-1] = [0, 0]
        _, _ = Bilinear(iava_rep, dims=(par['nx'], par['nt']),
                        dtype=par['dtype'])

    y = (Iop * x.ravel())
    assert_array_almost_equal(y, x[iava[0], iava[1]])

axs[1].imshow(drec1.T, cmap='gray', vmin=-clip, vmax=clip)
axs[1].set_title('Rec. from Seislet (%.1f %% coeffs.)' %
                 (100 * (nx - levels_cum[1]) / nx))
axs[1].axis('tight')
axs[2].imshow(d.T - drec1.T, cmap='gray', vmin=-clip, vmax=clip)
axs[2].set_title('Rec. error from Seislet')
axs[2].axis('tight')

############################################
# To conclude it is worth noting that the Seislet transform, opposite to the
# Wavelet transform, is not an orthogonal transformation: in other words,
# its adjoint and inverse are not equivalent. While we have used it the forward
# and inverse transformations, when used as part of a linear operator to be
# inverted the Seislet transform requires the forward-adjoint pair that is
# implemented in PyLops and passes the dot-test as shown below
pylops.utils.dottest(Sop, nt*nx, nt*nx, verb=True)