How to use the pylops.Diagonal function in pylops

smooth = np.ones(self.nsmooth) / self.nsmooth
            w = filtfilt(smooth, 1, w)

        # Create operators
        Rop = MDC(self.Rtwosided_fft, self.nt2, nv=1, dt=self.dt, dr=self.dr,
                  twosided=True, conj=False, transpose=False,
                  saveGt=self.saveRt, prescaled=self.prescaled,
                  dtype=self.dtype)
        R1op = MDC(self.Rtwosided_fft, self.nt2, nv=1, dt=self.dt, dr=self.dr,
                   twosided=True, conj=True, transpose=False,
                   saveGt=self.saveRt, prescaled=self.prescaled,
                   dtype=self.dtype)
        Rollop = Roll(self.nt2 * self.ns,
                      dims=(self.nt2, self.ns),
                      dir=0, shift=-1, dtype=self.dtype)
        Wop = Diagonal(w.T.flatten())
        Iop = Identity(self.nr * self.nt2)
        Mop = Block([[Iop, -1 * Wop * Rop],
                     [-1 * Wop * Rollop * R1op, Iop]]) * BlockDiag([Wop, Wop])
        Gop = Block([[Iop, -1 * Rop],
                     [-1 * Rollop * R1op, Iop]])

        if dottest:
            Dottest(Gop, 2 * self.ns * self.nt2,
                    2 * self.nr * self.nt2,
                    raiseerror=True, verb=True)
        if dottest:
            Dottest(Mop, 2 * self.ns * self.nt2,
                    2 * self.nr * self.nt2,
                    raiseerror=True, verb=True)

        # Create input focusing function

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as pltgs

import pylops

plt.close('all')

###############################################################################
# Let's define a diagonal operator :math:`\mathbf{d}` with increasing numbers from
# ``0`` to ``N`` and a unitary model :math:`\mathbf{x}`.
N = 10
d = np.arange(N)
x = np.ones(N)

Dop = pylops.Diagonal(d)

y = Dop*x
y1 = Dop.H*x

gs = pltgs.GridSpec(1, 6)
fig = plt.figure(figsize=(7, 3))
ax = plt.subplot(gs[0, 0:3])
im = ax.imshow(Dop.matrix(), cmap='rainbow', vmin=0, vmax=N)
ax.set_title('A', size=20, fontweight='bold')
ax.set_xticks(np.arange(N-1)+0.5)
ax.set_yticks(np.arange(N-1)+0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax.axis('tight')
ax = plt.subplot(gs[0, 3])

x = np.ones((nx, ny))
print('x =\n%s' % x)

# 1st dim
d = np.arange(nx)
Dop = pylops.Diagonal(d, dims=(nx, ny), dir=0)

y = Dop*x.flatten()
y1 = Dop.H*x.flatten()

print('1st dim: y = D*x =\n%s' % y.reshape(nx, ny))
print('1st dim: xadj = D\'*x =\n%s ' % y1.reshape(nx, ny))

# 2nd dim
d = np.arange(ny)
Dop = pylops.Diagonal(d, dims=(nx, ny), dir=1)

y = Dop*x.flatten()
y1 = Dop.H*x.flatten()

print('2nd dim: y = D*x =\n%s' % y.reshape(nx, ny))
print('2nd dim: xadj = D\'*x =\n%s ' % y1.reshape(nx, ny))

# create obliquity operator
    [Kx, F] = np.meshgrid(FFTop.f1, FFTop.f2, indexing='ij')
    k = F / vel
    Kz = np.sqrt((k ** 2 - Kx ** 2).astype(dtype))
    Kz[np.isnan(Kz)] = 0

    if composition:
        OBL = Kz / (rho * np.abs(F))
        OBL[F == 0] = 0
    else:
        OBL = rho * (np.abs(F) / Kz)
        OBL[Kz == 0] = 0

    # cut off and taper
    OBL = _filter_obliquity(OBL, F, Kx, vel, critical, ntaper)
    OBLop = Diagonal(OBL.ravel(), dtype=dtype)
    return FFTop, OBLop

# Let's start by defining a simple operator that applies element-wise
# multiplication of the model with a vector ``d`` in forward mode and
# element-wise multiplication of the data with the same vector ``d`` in
# adjoint mode. This operator is present in PyLops under the
# name of :py:class:`pylops.Diagonal` and
# its implementation is discussed in more details in the :ref:`addingoperator`
# page.
import timeit
import matplotlib.pyplot as plt
import numpy as np
import pylops

n = 10
d = np.arange(n) + 1.
x = np.ones(n)
Dop = pylops.Diagonal(d)

###############################################################################
# First of all we apply the operator in the forward mode. This can be done in
# four different ways:
#
# * ``_matvec``: directly applies the method implemented for forward mode
# * ``matvec``: performs some checks before and after applying ``_matvec``
# * ``*``: operator used to map the special method ``__matmul__`` which
#   checks whether the input ``x`` is a vector or matrix and applies ``_matvec``
#   or ``_matmul`` accordingly.
# * ``@``: operator used to map the special method ``__mul__`` which
#   performs like the ``*`` opetator
#
# We will time these 4 different executions and see how using ``_matvec``
# (or ``matvec``) will result in the faster computation. It is thus advised to
# use ``*`` (or ``@``) in examples when expressivity has priority but prefer

#    # create regularized normal eqs. and solve them simultaneously
            #    PP = np.dot(PPop.A.T, PPop.A) + epsI * np.eye(nt0*nm)
            #    datarn = PPop.A.T * datar.reshape(nt0*ntheta, nspatprod)
            #    PPop_reg = MatrixMult(PP, dims=ntheta*nspatprod)
            #    minv = lstsq(PPop_reg, datarn.flatten(), **kwargs_solver)[0]
        else:
            # solve unregularized normal equations simultaneously with lop
            minv = lsqr(PPop, datar, **kwargs_solver)[0]
    else:
        # Create Thicknov regularization
        if epsI is not None:
            if isinstance(epsI, (list, tuple)):
                if len(epsI) != nm:
                    raise ValueError('epsI must be a scalar or a list of'
                                     'size nm')
                RegI = Diagonal(np.array(epsI),
                                dims=(nt0, nm,nspatprod), dir=1)
            else:
                RegI = epsI * Identity(nt0 * nm * nspatprod)

        if epsRL1 is None:
            # L2 inversion with spatial regularization
            if dims == 1:
                Regop = SecondDerivative(nt0 * nm, dtype=PPop.dtype,
                                         dims=(nt0, nm))
            elif dims == 2:
                Regop = Laplacian((nt0, nm, nx), dirs=(0, 2), dtype=PPop.dtype)
            else:
                Regop = Laplacian((nt0, nm, nx, ny), dirs=(2, 3),
                                  dtype=PPop.dtype)
            if epsI is None:
                Regop = (Regop, )