__all__ = ["Identity"]
import numpy as np
from pylops import LinearOperator
from pylops.utils.backend import get_array_module, inplace_set
from pylops.utils.decorators import reshaped
from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray
[docs]
class Identity(LinearOperator):
r"""Identity operator.
Simply move model to data in forward model and viceversa in adjoint mode if
:math:`M = N`. If :math:`M > N` removes last :math:`M - N` elements from
model in forward and pads with :math:`0` in adjoint. If :math:`N > M`
removes last :math:`N - M` elements from data in adjoint and pads with
:math:`0` in forward.
Note that the identity operator can handle both 1d and nd arrays; in the
case of nd arrays, all elements of N must be larger or equal than those of M
(or all elements of M must be larger or equal than those of N).
Parameters
----------
N : :obj:`int` or :obj:`tuple`
Number of samples in data (and model, if ``M`` is not provided).
If a tuple is provided, this is interpreted as the data (and model)
are nd-arrays.
M : :obj:`int` or :obj:`tuple`, optional
Number of samples in model. If a tuple is provided, this is interpreted
as the model is an nd-array. Note that when `M` is a tuple, `N` must be
also a tuple with the same number of elements.
inplace : :obj:`bool`, optional
Work inplace (``True``) or make a new copy (``False``). By default,
data is a reference to the model (in forward) and model is a reference
to the data (in adjoint).
forceflat : :obj:`bool`, optional
.. versionadded:: 2.2.0
Force an array to be flattened after matvec and rmatvec. Note that this is only
required when `N` and `M` are tuples (input and output arrays are nd-arrays).
dtype : :obj:`str`, optional
Type of elements in input array.
name : :obj:`str`, optional
.. versionadded:: 2.0.0
Name of operator (to be used by :func:`pylops.utils.describe.describe`)
Attributes
----------
dims : :obj:`tuple`
Shape of the array after the adjoint, but before flattening.
For example, ``x_reshaped = (Op.H * y.ravel()).reshape(Op.dims)``.
dimsd : :obj:`tuple`
Shape of the array after the forward, but before flattening.
For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dimsd)``.
mode : :obj:`str`
Operation mode: ``same`` for :math:`M = N`,
``model`` for :math:`N > M,
and ``data`` for :math:`M > N`.
shape : :obj:`tuple`
Operator shape.
Raises
------
ValueError
- If ``M`` is a tuple with different number of elements of ``N``
- If ``N`` ``M`` are non-identical tuples and some values are largers
in ``N`` and some in ``M``
NotImplementedError
If ``N`` or ``M`` have type different from int or tuple/list
Notes
-----
For :math:`M = N`, an *Identity* operator simply moves the model
:math:`\mathbf{x}` to the data :math:`\mathbf{y}` in forward mode and
viceversa in adjoint mode:
.. math::
y_i = x_i \quad \forall i=1,2,\ldots,N
or in matrix form:
.. math::
\mathbf{y} = \mathbf{I} \mathbf{x} = \mathbf{x}
and
.. math::
\mathbf{x} = \mathbf{I} \mathbf{y} = \mathbf{y}
For :math:`M > N`, the *Identity* operator takes the first :math:`M`
elements of the model :math:`\mathbf{x}` into the data :math:`\mathbf{y}`
in forward mode
.. math::
y_i = x_i \quad \forall i=1,2,\ldots,N
and all the elements of the data :math:`\mathbf{y}` into the first
:math:`M` elements of model in adjoint mode (other elements are ``O``):
.. math::
x_i = y_i \quad \forall i=1,2,\ldots,M
x_i = 0 \quad \forall i=M+1,\ldots,N
"""
def __init__(
self,
N: int | InputDimsLike,
M: int | InputDimsLike | None = None,
inplace: bool = True,
forceflat: bool | None = None,
dtype: DTypeLike = "float64",
name: str = "I",
) -> None:
M = N if M is None else M
if isinstance(N, int) and isinstance(M, int):
# N and M are scalars (1d-arrays)
super().__init__(
dtype=np.dtype(dtype),
dims=(M,),
dimsd=(N,),
forceflat=forceflat,
name=name,
)
# identify behaviour for matvec/rmatvec: 'same' for N=M,
# 'data' for N>M, and 'model' for M>N
if N == M:
self.mode = "same"
elif N < M:
self.mode = "model"
self._sliceN = slice(0, N)
self._sliceM = slice(0, M)
else:
self.mode = "data"
self._sliceN = slice(0, N)
self._sliceM = slice(0, M)
elif isinstance(N, (tuple, list)) and isinstance(M, (tuple, list)):
# N and M are tuples (nd-arrays)
# First check that all elements in N and M are the same or that
# all elements of either N or M are bigger than the other one and
# raise error is not the case
if np.array_equal(N, M):
self.mode = "same"
elif np.array_equal(M, np.maximum(N, M)):
self.mode = "model"
self._sliceN = tuple([slice(0, n) for n in N])
self._sliceM = tuple([slice(0, m) for m in M])
elif np.array_equal(N, np.maximum(N, M)):
self.mode = "data"
self._sliceN = tuple([slice(0, n) for n in N])
self._sliceM = tuple([slice(0, m) for m in M])
else:
msg = (
"N and M are not identical, "
"and some values are larger in N and some in M"
)
raise ValueError(msg)
super().__init__(
dtype=np.dtype(dtype), dims=M, dimsd=N, forceflat=forceflat, name=name
)
else:
msg = (
"N and M must have same type and equal to "
"int, tuple, or list, instead their types"
f" are type(N)={type(N)} and type(M)={type(M)}"
)
raise NotImplementedError(msg)
self.inplace = inplace
@reshaped
def _matvec(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
if not self.inplace:
x = x.copy()
if self.mode == "same":
y = x
elif self.mode == "model":
y = x[self._sliceN]
else:
y = ncp.zeros(self.dimsd, dtype=self.dtype)
y = inplace_set(x, y, self._sliceM)
return y
@reshaped
def _rmatvec(self, x: NDArray) -> NDArray:
ncp = get_array_module(x)
if not self.inplace:
x = x.copy()
if self.mode == "same":
y = x
elif self.mode == "model":
y = ncp.zeros(self.dims, dtype=self.dtype)
y = inplace_set(x, y, self._sliceN)
else:
y = x[self._sliceM]
return y
PyLops