How to use the pylops.basicoperators.MatrixMult function in pylops

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.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(par['nx']) + par['imag']*np.ones(par['nx'])

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

    xlsqr = lsqr(Vop, Vop * 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
    V1op = VStack([G1, MatrixMult(G2, dtype=par['dtype'])],
                  dtype=par['dtype'])
    assert dottest(V1op, 2 * par['ny'], 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')
    V2op = VStack([G1, MatrixMult(G2, dtype=par['dtype'])],
                  dtype=par['dtype'])
    assert dottest(V2op, 2 * par['ny'], par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

xlsqr = lsqr(Bop, Bop * 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
    B1op = Block([[G11,
                   MatrixMult(G12, dtype=par['dtype'])],
                  [MatrixMult(G21, dtype=par['dtype']),
                   G22]], dtype=par['dtype'])
    assert dottest(B1op, 2 * par['ny'], 2 * par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

    # use scipy matrix directly in the definition of the operator
    G11 = sp_random(par['ny'], par['nx'], density=0.4).astype('float32')
    B2op = Block([[G11,
                   MatrixMult(G12, dtype=par['dtype'])],
                  [MatrixMult(G21, dtype=par['dtype']),
                   G22]], dtype=par['dtype'])
    assert dottest(B2op, 2 * par['ny'], 2 * par['nx'],
                   complexflag=0 if par['imag'] == 0 else 3)

def test_Sliding2D(par):
    """Dot-test and inverse for Sliding2D operator
    """
    Op = MatrixMult(np.ones((par['nwiny'] * par['nt'], par['ny'] * par['nt'])))

    Slid = Sliding2D(Op, dims=(par['ny']*par['winsy'], par['nt']),
                     dimsd=(par['npy'], par['nt']),
                     nwin=par['nwiny'], nover=par['novery'],
                     tapertype=par['tapertype'])
    assert dottest(Slid, par['npy']*par['nt'],
                   par['ny']*par['nt']*par['winsy'])
    x = np.ones((par['ny']*par['winsy'], par['nt']))
    y = Slid * x.flatten()

    xinv = LinearOperator(Slid) / y
    assert_array_almost_equal(x.flatten(), xinv)

def test_VStack_incosistent_columns(par):
    """Check error is raised if operators with different number of columns
    are passed to VStack
    """
    G1 = np.random.normal(0, 10, (par['ny'], par['nx'])).astype(par['dtype'])
    G2 = np.random.normal(0, 10, (par['ny'], par['nx'] + 1)).astype(par['dtype'])
    with pytest.raises(ValueError):
        VStack([MatrixMult(G1, dtype=par['dtype']),
                MatrixMult(G2, dtype=par['dtype'])],
               dtype=par['dtype'])

def test_NormalEquationsInversion(par):
    """Solve normal equations in least squares sense
    """
    np.random.seed(10)
    G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
        par['imag']*np.random.normal(0, 10,
                                     (par['ny'], par['nx'])).astype('float32')
    Gop = MatrixMult(G, dtype=par['dtype'])

    Reg = MatrixMult(np.eye(par['nx']), dtype=par['dtype'])
    NReg = MatrixMult(np.eye(par['nx']), dtype=par['dtype'])
    Weigth = Diagonal(np.ones(par['ny']), dtype=par['dtype'])
    x = np.ones(par['nx']) + par['imag']*np.ones(par['nx'])
    x0 = np.random.normal(0, 10, par['nx']) + \
         par['imag']*np.random.normal(0, 10, par['nx']) if par['x0'] else None
    y = Gop * x

    # normal equations with regularization
    xinv = NormalEquationsInversion(Gop, [Reg], y, epsI=0,
                                    epsRs=[1e-8], x0=x0,
                                    returninfo=False,
                                    **dict(maxiter=200, tol=1e-10))
    assert_array_almost_equal(x, xinv, decimal=3)
    # normal equations with weight

def test_NormalEquationsInversion(par):
    """Solve normal equations in least squares sense
    """
    np.random.seed(10)
    G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
        par['imag']*np.random.normal(0, 10,
                                     (par['ny'], par['nx'])).astype('float32')
    Gop = MatrixMult(G, dtype=par['dtype'])

    Reg = MatrixMult(np.eye(par['nx']), dtype=par['dtype'])
    NReg = MatrixMult(np.eye(par['nx']), dtype=par['dtype'])
    Weigth = Diagonal(np.ones(par['ny']), dtype=par['dtype'])
    x = np.ones(par['nx']) + par['imag']*np.ones(par['nx'])
    x0 = np.random.normal(0, 10, par['nx']) + \
         par['imag']*np.random.normal(0, 10, par['nx']) if par['x0'] else None
    y = Gop * x

    # normal equations with regularization
    xinv = NormalEquationsInversion(Gop, [Reg], y, epsI=0,
                                    epsRs=[1e-8], x0=x0,
                                    returninfo=False,
                                    **dict(maxiter=200, tol=1e-10))
    assert_array_almost_equal(x, xinv, decimal=3)
    # normal equations with weight
    xinv = NormalEquationsInversion(Gop, None, y, Weight=Weigth, epsI=0,
                                    x0=x0, returninfo=False,
                                    **dict(maxiter=200, tol=1e-10))

def test_RegularizedInversion(par):
    """Solve regularized inversion in least squares sense
    """
    np.random.seed(10)
    G = np.random.normal(0, 10, (par['ny'], par['nx'])).astype('float32') + \
        par['imag']*np.random.normal(0, 10,
                                     (par['ny'], par['nx'])).astype('float32')
    Gop = MatrixMult(G, dtype=par['dtype'])
    Reg = MatrixMult(np.eye(par['nx']), dtype=par['dtype'])
    Weigth = Diagonal(np.ones(par['ny']), dtype=par['dtype'])
    x = np.ones(par['nx']) + par['imag']*np.ones(par['nx'])
    x0 = np.random.normal(0, 10, par['nx']) + \
         par['imag']*np.random.normal(0, 10, par['nx']) if par['x0'] else None
    y = Gop*x

    # regularized inversion with regularization
    xinv = RegularizedInversion(Gop, [Reg], y, epsRs=[1e-8], x0=x0,
                                returninfo=False,
                                **dict(damp=0, iter_lim=200, show=0))
    assert_array_almost_equal(x, xinv, decimal=3)
    # regularized inversion with weight
    xinv = RegularizedInversion(Gop, None, y, Weight=Weigth,
                                x0=x0,
                                returninfo=False,

def __init__(self, ops, dtype=None):
        self.ops = ops
        nops = np.zeros(len(self.ops), dtype=np.int)
        for iop, oper in enumerate(ops):
            if not isinstance(oper, (LinearOperator, spLinearOperator)):
                self.ops[iop] = MatrixMult(oper, dtype=oper.dtype)
            nops[iop] = self.ops[iop].shape[0]
        self.nops = nops.sum()
        mops = [oper.shape[1] for oper in self.ops]
        if len(set(mops)) > 1:
            raise ValueError('operators have different number of columns')
        self.mops = mops[0]
        self.nnops = np.insert(np.cumsum(nops), 0, 0)
        self.shape = (self.nops, self.mops)
        if dtype is None:
            self.dtype = _get_dtype(self.ops)
        else:
            self.dtype = np.dtype(dtype)
        self.explicit = False

# create sinc interpolation matrix
    nreg = M if dims is None else dims[dir]
    ireg = np.arange(nreg)
    sinc = np.tile(iava[:, np.newaxis], (1, nreg)) - \
           np.tile(ireg, (len(iava), 1))
    sinc = np.sinc(sinc)

    # identify additional dimensions and create MatrixMult operator
    otherdims = None
    if dims is not None:
        otherdims = np.array(dims)
        otherdims = np.roll(otherdims, -dir)
        otherdims = otherdims[1:]
        print('dims', dims)
        print('otherdims', otherdims)
    Op = MatrixMult(sinc, dims=otherdims, dtype=dtype)

    # create Transpose operator that brings dir to first dimension
    if dir > 0:
        axes = np.arange(len(dims), dtype=np.int)
        axes = np.roll(axes, -dir)
        dimsd =  list(dims)
        dimsd[dir] = len(iava)
        Top = Transpose(dims, axes=axes, dtype=dtype)
        T1op = Transpose(dimsd, axes=axes, dtype=dtype)
        Op = T1op.H * Op * Top
    return Op