How to use the pylops.Identity function in pylops

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
        if G0 is None:

w = filtfilt(smooth, 1, w)

        # Create operators
        Rop = MDC(self.Rtwosided_fft, self.nt2, nv=nvs,
                  dt=self.dt, dr=self.dr, twosided=True,
                  conj=False, transpose=False, prescaled=self.prescaled,
                  dtype=self.dtype)
        R1op = MDC(self.Rtwosided_fft, self.nt2, nv=nvs,
                   dt=self.dt, dr=self.dr, twosided=True,
                   conj=True, transpose=False, prescaled=self.prescaled,
                   dtype=self.dtype)
        Rollop = Roll(self.ns * nvs * self.nt2,
                      dims=(self.nt2, self.ns, nvs),
                      dir=0, shift=-1, dtype=self.dtype)
        Wop = Diagonal(w.transpose(2, 0, 1).flatten())
        Iop = Identity(self.nr * nvs * 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.nr * nvs * self.nt2,
                    2 * self.nr * nvs * self.nt2,
                    raiseerror=True, verb=True)
        if dottest:
            Dottest(Mop, 2 * self.ns * nvs * self.nt2,
                    2 * self.nr * nvs * self.nt2,
                    raiseerror=True, verb=True)

        # Create input focusing function
        if G0 is None:

"""
    nffts = (int(nffts[0]) if nffts[0] is not None else nr[0],
             int(nffts[1]) if nffts[1] is not None else nr[1],
             int(nffts[2]) if nffts[2] is not None else nt)

    # create obliquity operator
    FFTop, OBLop = \
        _obliquity3D(nt, nr, dt, dr, rho, vel,
                     nffts=nffts,
                     critical=critical, ntaper=ntaper,
                     composition=True, dtype=dtype)

    # create up-down modelling operator
    UDop = (BlockDiag([FFTop.H, scaling * FFTop.H]) * \
            Block([[Identity(nffts[0] * nffts[1] * nffts[2], dtype=dtype),
                    Identity(nffts[0] * nffts[1] * nffts[2], dtype=dtype)],
                   [OBLop, -OBLop]]) * \
            BlockDiag([FFTop, FFTop]))
    return UDop

# Let's define now the sampling operator as well as create our covariance
# matrices in terms of linear operators. This may not be strictly necessary
# here but shows how even Bayesian-type of inversion can very easily scale to
# large model and data spaces.

# Sampling operator
perc_subsampling = 0.2
ntsub = int(np.round(nt*perc_subsampling))
iava = np.sort(np.random.permutation(np.arange(nt))[:ntsub])
iava[-1] = nt-1 # assume we have the last sample to avoid instability
Rop = pylops.Restriction(nt, iava, dtype='float64')

# Covariance operators
Cm_op = \
    pylops.signalprocessing.Convolve1D(nt, diag_ave, offset=N)
Cd_op = sigmad**2 * pylops.Identity(ntsub)

###############################################################################
# We model now our data and add noise that respects our prior definition
n = np.random.normal(0, sigmad, nt)
y = Rop * x
yn = Rop * (x + n)
ymask = Rop.mask(x)
ynmask = Rop.mask(x + n)

###############################################################################
# First we apply the Bayesian inversion equation
xbayes = x0 + Cm_op * Rop.H * (lsqr(Rop * Cm_op * Rop.H + Cd_op,
                                    yn - Rop*x0, iter_lim=400)[0])

# Visualize
fig, ax = plt.subplots(1, 1, figsize=(12, 5))

ax.axis('off')
ax = plt.subplot(gs[0, 5])
ax.imshow(y[:, np.newaxis], cmap='rainbow')
ax.set_title('y', size=20, fontweight='bold')
ax.set_xticks([])
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color='white')
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
fig.colorbar(im, ax=ax, ticks=[0, 1], pad=0.3, shrink=0.7)

###############################################################################
# Similarly we can consider the case with data bigger than model
N, M = 10, 5
x = np.arange(M)
Iop = pylops.Identity(N, M, dtype='int')

y = Iop*x
xadj = Iop.H*y

print('x = %s ' % x)
print('I*x = %s ' % y)
print('I\'*y = %s ' % xadj)

###############################################################################
# and model bigger than data
N, M = 5, 10
x = np.arange(M)
Iop = pylops.Identity(N, M, dtype='int')

y = Iop*x
xadj = Iop.H*y

N, M = 10, 5
x = np.arange(M)
Iop = pylops.Identity(N, M, dtype='int')

y = Iop*x
xadj = Iop.H*y

print('x = %s ' % x)
print('I*x = %s ' % y)
print('I\'*y = %s ' % xadj)

###############################################################################
# and model bigger than data
N, M = 5, 10
x = np.arange(M)
Iop = pylops.Identity(N, M, dtype='int')

y = Iop*x
xadj = Iop.H*y

print('x = %s ' % x)
print('I*x = %s ' % y)
print('I\'*y = %s ' % xadj)

def _MDC(G, nt, nv, dt=1., dr=1., twosided=True, fast=None, dtype=None,
         transpose=True, saveGt=True, conj=False, prescaled=False,
         _Identity=Identity, _Transpose=Transpose, _FFT=FFT,
         _Fredholm1=Fredholm1, args_Identity={}, args_Transpose={},
         args_FFT={}, args_Identity1={}, args_Transpose1={},
         args_FFT1={}, args_Fredholm1={}):
    r"""Multi-dimensional convolution.

    Used to be able to provide operators from different libraries to
    MDC. It operates in the same way as public method
    (PoststackLinearModelling) but has additional input parameters allowing
    passing a different operator and additional arguments to be passed to such
    operator.

    """
    warnings.warn('A new implementation of MDC is provided in v1.5.0. This '
                  'currently affects only the inner working of the operator, '
                  'end-users can continue using the operator in the same way. '
                  'Nevertheless, it is now recommended to start using the '

into data and viceversa.
"""
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as pltgs

import pylops

plt.close('all')

###############################################################################
# Let's define an identity operator :math:`\mathbf{Iop}` with same number of
# elements for data and model (:math:`N=M`).
N, M = 5, 5
x = np.arange(M)
Iop = pylops.Identity(M, dtype='int')

y = Iop*x
xadj = Iop.H*y

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