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

def test_eigs(par):
    """Eigenvalues and condition number estimate with ARPACK
    """
    # explicit=True
    diag = np.arange(par['nx'], 0, -1) +\
           par['imag'] * np.arange(par['nx'], 0, -1)
    Op = MatrixMult(np.vstack((np.diag(diag),
                               np.zeros((par['ny'] - par['nx'], par['nx'])))))
    eigs = Op.eigs()
    assert_array_almost_equal(diag[:eigs.size], eigs, decimal=3)

    cond = Op.cond()
    assert_array_almost_equal(np.real(cond), par['nx'], decimal=3)

    # explicit=False
    Op = Diagonal(diag, dtype=par['dtype'])
    if par['ny'] > par['nx']:
        Op = VStack([Op, Zero(par['ny'] - par['nx'], par['nx'])])
    eigs = Op.eigs()
    assert_array_almost_equal(diag[:eigs.size], eigs, decimal=3)

    # uselobpcg cannot be used for square non-symmetric complex matrices
    if np.iscomplex(Op):
        eigs1 = Op.eigs(uselobpcg=True)
        assert_array_almost_equal(eigs, eigs1, decimal=3)

    cond = Op.cond()
    assert_array_almost_equal(np.real(cond), par['nx'], decimal=3)

    # uselobpcg cannot be used for square non-symmetric complex matrices
    if np.iscomplex(Op):
        cond1 = Op.cond(uselobpcg=True, niter=100)

def test_sparse(par):
    """Sparse matrix representation
    """
    diag = np.arange(par['nx']) +\
           par['imag'] * np.arange(par['nx'])
    D = np.diag(diag)
    Dop = Diagonal(diag, dtype=par['dtype'])
    S = Dop.tosparse()
    assert_array_equal(S.A, D)

def test_overloads(par):
    """Apply various overloaded operators (.H, -, +, *) and ensure that the
    returned operator is still of pylops LinearOperator type
    """
    diag = np.arange(par['nx']) +\
           par['imag'] * np.arange(par['nx'])
    Dop = Diagonal(diag, dtype=par['dtype'])

    # .H
    assert isinstance(Dop.H, LinearOperator)
    # negate
    assert isinstance(-Dop, LinearOperator)
    # multiply by scalar
    assert isinstance(2*Dop, LinearOperator)
    # +
    assert isinstance(Dop + Dop, LinearOperator)
    # -
    assert isinstance(Dop - 2*Dop, LinearOperator)
    # *
    assert isinstance(Dop * Dop, LinearOperator)
    # **
    assert isinstance(Dop **2, LinearOperator)

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,
                                **dict(damp=0, iter_lim=200, show=0))
    assert_array_almost_equal(x, xinv, decimal=3)

str(mwin_ins), str(mwin_ends))
        logging.warning('data wins - start:%s, end:%s',
                        str(dwin_ins), str(dwin_ends))
    if nwins*Op.shape[1]//dims[1] != dims[0]:
        raise ValueError('Model shape (dims=%s) is not consistent with chosen '
                         'number of windows. Choose dims[0]=%d for the '
                         'operator to work with estimated number of windows, '
                         'or create the operator with design=True to find '
                         'out the optimal number of windows for the current '
                         'model size...'
                         % (str(dims), nwins*Op.shape[1]//dims[1]))
    # transform to apply
    if tapertype is None:
        OOp = BlockDiag([Op for _ in range(nwins)])
    else:
        OOp = BlockDiag([Diagonal(taps[itap].flatten()) * Op
                         for itap in range(nwins)])

    combining = HStack([Restriction(np.prod(dimsd), range(win_in, win_end),
                                    dims=dimsd).H
                        for win_in, win_end in zip(dwin_ins, dwin_ends)])
    Sop = combining * OOp
    return Sop

:math:`\mathbf{v}(x, y)`:

    .. math::
        df_\mathbf{v}(x,y) =
            \nabla f(x,y) \mathbf{v}(x,y)

    where we have here considered the 2-dimensional case.

    This operator can be easily implemented as the concatenation of the
    :py:class:`pylops.Gradient` operator and the :py:class:`pylops.Diagonal`
    operator with :math:`\mathbf{v}` along the main diagonal.

    """
    Gop = Gradient(dims, sampling=sampling, edge=edge, kind=kind, dtype=dtype)
    if v.ndim == 1:
        Dop = Diagonal(v, dims=[len(dims)]+list(dims), dir=0, dtype=dtype)
    else:
        Dop = Diagonal(v.ravel(), dtype=dtype)
    Sop = Sum(dims=[len(dims)]+list(dims), dir=0, dtype=dtype)
    ddop = Sop * Dop * Gop
    return LinearOperator(ddop)

.. math::
        df_\mathbf{v}(x,y) =
            \nabla f(x,y) \mathbf{v}(x,y)

    where we have here considered the 2-dimensional case.

    This operator can be easily implemented as the concatenation of the
    :py:class:`pylops.Gradient` operator and the :py:class:`pylops.Diagonal`
    operator with :math:`\mathbf{v}` along the main diagonal.

    """
    Gop = Gradient(dims, sampling=sampling, edge=edge, kind=kind, dtype=dtype)
    if v.ndim == 1:
        Dop = Diagonal(v, dims=[len(dims)]+list(dims), dir=0, dtype=dtype)
    else:
        Dop = Diagonal(v.ravel(), dtype=dtype)
    Sop = Sum(dims=[len(dims)]+list(dims), dir=0, dtype=dtype)
    ddop = Sop * Dop * Gop
    return LinearOperator(ddop)

if nwins*Op.shape[1]//dims[2] != dims[0]*dims[1]:
        raise ValueError('Model shape (dims=%s) is not consistent with chosen '
                         'number of windows. Choose dims[0]=%d and '
                         'dims[1]=%d for the operator to work with '
                         'estimated number of windows, or create '
                         'the operator with design=True to find out the'
                         'optimal number of windows for the current '
                         'model size...'
                         % (str(dims), nwins0*Op.shape[1]//(nop[1]*dims[2]),
                            nwins1 * Op.shape[1]//(nop[0]*dims[2])))
    # transform to apply
    if tapertype is None:
        OOp = BlockDiag([Op for _ in range(nwins)])
    else:
        OOp = BlockDiag([Diagonal(tap.flatten()) * Op
                         for _ in range(nwins)])

    hstack = HStack([Restriction(dimsd[1] * dimsd[2] * nwin[0],
                                 range(win_in, win_end),
                                 dims=(nwin[0], dimsd[1], dimsd[2]),
                                 dir=1).H
                     for win_in, win_end in zip(dwin1_ins,
                                                dwin1_ends)])

    combining1 = BlockDiag([hstack]*nwins0)
    combining0 = HStack([Restriction(np.prod(dimsd),
                                     range(win_in, win_end),
                                     dims=dimsd, dir=0).H
                         for win_in, win_end in zip(dwin0_ins, dwin0_ends)])
    Sop = combining0 * combining1 * OOp
    return Sop

if epsRs is None and Regs is not None:
        epsRs = [1] * len(Regs)

    # Normal equations
    if Weight is not None:
        y_normal = OpH * Weight * data
    else:
        y_normal = OpH * data
    if Weight is not None:
        Op_normal = OpH * Weight * Op
    else:
        Op_normal = OpH * Op

    # Add regularization terms
    if epsI > 0:
        Op_normal += epsI ** 2 * Diagonal(np.ones(Op.shape[1]))

    if Regs is not None:
        for epsR, Reg, datareg in zip(epsRs, Regs, dataregs):
            RegH = Reg.H
            y_normal += epsR ** 2 * RegH * datareg
            Op_normal += epsR ** 2 * RegH * Reg

    if NRegs is not None:
        for epsNR, NReg in zip(epsNRs, NRegs):
            Op_normal += epsNR ** 2 * NReg

    # CG solver
    if x0 is not None:
        y_normal = y_normal - Op_normal*x0
    xinv, istop = cg(Op_normal, y_normal, **kwargs_cg)
    if x0 is not None:

outside = (iava >= lastsample - 1)
    if sum(outside) > 0:
        logging.warning('at least one value is beyond penultimate sample, '
                        'forced to be at penultimate sample')
    iava[outside] = lastsample - 1 - 1e-10
    _checkunique(iava)

    # find indices and weights
    iva_l = np.floor(iava).astype(np.int)
    iva_r = iva_l + 1
    weights = iava - iva_l

    # create operators
    Op = Diagonal(1 - weights, dims=dimsd, dir=dir, dtype=dtype) * \
         Restriction(M, iva_l, dims=dims, dir=dir, dtype=dtype) + \
         Diagonal(weights, dims=dimsd, dir=dir, dtype=dtype) * \
         Restriction(M, iva_r, dims=dims, dir=dir, dtype=dtype)
    return Op, iava