Implementing new operators#
Users are welcome to create new operators and add them to the PyLops library.
In this tutorial, we will go through the key steps in the definition of an operator, using a simplified version of the
pylops.Diagonal
operator as an example. This is a very simple operator that applies a diagonal matrix
to the model in forward mode and to the data in adjoint mode.
Creating the operator#
The first thing we need to do is to create a new file with the name of the operator we would like to implement. Note that as the operator will be a class, we need to follow the UpperCaseCamelCase convention both for the class itself and for the filename.
Once we have created the file, we will start by importing the modules that will be needed by the operator.
While this varies from operator to operator, you will always need to import the pylops.LinearOperator
class,
which will be used as parent class for any of our operators:
from pylops import LinearOperator
This class was originally defined as a child of the scipy.sparse.linalg.LinearOperator
, implementing
the same methods of its parent class as well as additional methods for quick inversion (e.g., A \ y
),
eigenvalues computation, conversion to dense matrices, etc. From version v2.1
onwards, our linear operator
class has become stand-alone; however, by keeping the same naming structure of the original scipy class, we still
allow inter-operability with scipy-native linear operators and solvers.
After that we define our new object:
class Diagonal(LinearOperator):
followed by a numpydoc docstring
(starting with r"""
and ending with """
) containing the documentation of the operator. Such docstring should
contain at least a short description of the operator, a Parameters
section with a detailed description of the
input parameters and a Notes
section providing a mathematical explanation of the operator. Take a look at
some of the core operators of PyLops to get a feeling of the level of details of the mathematical explanation.
Initialization (__init__
)#
We then need to create the __init__
where the input parameters are passed and saved as members of our class.
While the input parameters change from operator to operator, it is always required to create three members:
dtype
: data type object (of typestr
ornp.dtype
) of the model and data;shape
: a tuple containing the dimensions of the operator in the data and model space;explicit
: a boolean (True
orFalse
) identifying if the operator can be inverted by a direct solver or requires an iterative solver. This member isTrue
if the operator has also a memberA
that contains the matrix to be inverted like for example in thepylops.MatrixMult
operator, and it will beFalse
otherwise.
In this specific case, we have another member called d
which is equal to the input vector containing the diagonal
elements of the matrix we want to multiply to the model and data.
def __init__(self, d, dtype=None):
self.d = d.ravel()
self.shape = (len(self.d), len(self.d))
self.dtype = np.dtype(dtype)
self.explicit = False
Alternatively, since version v2.0.0
, the recommended way of initializing operators derived from the base
pylops.LinearOperator
class is to invoke super
to assign the required attributes:
def __init__(self, d, dtype=None):
self.d = d.ravel()
super().__init__(dtype=np.dtype(dtype), shape=(len(self.d), len(self.d)))
In this case, there is no need to declare explicit
as it already defaults to False
.
Moreover, since version v2.0.0
, every pylops.LinearOperator
class is imbued with dims
,
dimsd
, and clinear
in addition to the required dtype
, shape
, and explicit
. Note that
dims
and dimsd
can be defined in spite of shape
, which will be automatically assigned within the
super
method: the main difference between dims
/dimsd
and shape
is the the former variables can be
used the define the n-dimensional nature of the input of an operator, whilst the latter variable refers to their overall
shape when the input is flattened.
See the docs of pylops.LinearOperator
for more information about what these
attributes mean.
Forward mode (_matvec
)#
We can then move onto writing the forward mode in the method _matvec
. In other words, we will need to write
the piece of code that will implement the following operation \(\mathbf{y} = \mathbf{A}\mathbf{x}\).
Such method is always composed of two inputs (the object itself self
and the input model x
).
In our case the code to be added to the forward is very simple, we will just need to apply element-wise multiplication
between the model \(\mathbf{x}\) and the elements along the diagonal contained in the array \(\mathbf{d}\).
We will finally need to return
the result of this operation:
def _matvec(self, x):
return self.d * x
Note that since version v2.0.0
, this method can be decorated by the decorator @reshaped
. As discussed in
more details in the decorator documentation, by adding such decorator the input x
is initially reshaped into
a nd-array of shape dims
, fed to the actual code in _matvec
and then flattened.
Adjoint mode (_rmatvec
)#
Finally we need to implement the adjoint mode in the method _rmatvec
. In other words, we will need to write
the piece of code that will implement the following operation \(\mathbf{x} = \mathbf{A}^H\mathbf{y}\).
Such method is also composed of two inputs (the object itself self
and the input data y
).
In our case the code to be added to the forward is the same as the one from the forward (but this will obviously be
different from operator to operator):
def _rmatvec(self, x):
return self.d * x
And that’s it, we have implemented our first linear operator!
Similar to _matvec
, since version v2.0.0
, this method can also be decorated by the decorator @reshaped
.
When doing so, the input x
is initially reshaped into
a nd-array of shape dimsd
, fed to the actual code in _rmatvec
and then flattened.
Testing the operator#
Being able to write an operator is not yet a guarantee of the fact that the operator is correct, or in other words that the adjoint code is actually the adjoint of the forward code. Luckily for us, a simple test can be performed to check the validity of forward and adjoint operators, the so called dot-test.
We can generate random vectors \(\mathbf{u}\) and \(\mathbf{v}\) and verify the the following equality within a numerical tolerance:
The method pylops.utils.dottest
implements such a test for you, all you need to do is create a new test
within an existing test_*.py
file in the pytests
folder (or in a new file).
Generally a test file will start with a number of dictionaries containing different parameters we would like to use in the testing of one or more operators. The test itself starts with a decorator that contains a list of all (or some) of dictionaries that will would like to use for our specific operator, followed by the definition of the test
@pytest.mark.parametrize("par", [(par1),(par2)])
def test_Diagonal(par):
At this point we can first of all create the operator and run the pylops.utils.dottest
preceded by the
assert
command. Moreover, the forward and adjoint methods should tested towards expected outputs or even
better, when the operator allows it (i.e., operator is invertible), a small inversion should be run and the inverted
model tested towards the input model.
"""Dot-test and inversion for diagonal operator
"""
d = np.arange(par['nx']) + 1.
Dop = Diagonal(d)
assert dottest(Dop, par['nx'], par['nx'],
complexflag=0 if par['imag'] == 1 else 3)
x = np.ones(par['nx'])
xlsqr = lsqr(Dop, Dop * x, damp=1e-20, iter_lim=300, show=0)[0]
assert_array_almost_equal(x, xlsqr, decimal=4)
Documenting the operator#
Once the operator has been created, we can add it to the documentation of PyLops. To do so, simply add the name of
the operator within the index.rst
file in docs/source/api
directory.
Moreover, in order to facilitate the user of your operator by other users, a simple example should be provided as part of the
Sphinx-gallery of the documentation of the PyLops library. The directory examples
contains several scripts that
can be used as template.
Final checklist#
Before submitting your new operator for review, use the following checklist to ensure that your code adheres to the guidelines of PyLops:
you have created a new file containing a single class (or a function when the new operator is a simple combination of existing operators - see
pylops.Laplacian
for an example of such operator) and added to a new or existing directory within thepylops
package.the new class contains at least
__init__
,_matvec
and_matvec
methods.the new class (or function) has a numpydoc docstring documenting at least the input
Parameters
and with aNotes
section providing a mathematical explanation of the operatora new test has been added to an existing
test_*.py
file within thepytests
folder. The test should verify that the new operator passes thepylops.utils.dottest
. Moreover it is advisable to create a small toy example where the operator is applied in forward mode and the resulting data is inverted using\
frompylops.LinearOperator
.the new operator is used within at least one example (in
examples
directory) or one tutorial (intutorials
directory).