# 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 type`str`

or`np.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`

or`False`

) identifying if the operator can be inverted by a direct solver or requires an iterative solver. This member is`True`

if the operator has also a member`A`

that contains the matrix to be inverted like for example in the`pylops.MatrixMult`

operator, and it will be`False`

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 the`pylops`

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 a`Notes`

section providing a mathematical explanation of the operatora new test has been added to an existing

`test_*.py`

file within the`pytests`

folder. The test should verify that the new operator passes the`pylops.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`\`

from`pylops.LinearOperator`

.the new operator is used within at least one

*example*(in`examples`

directory) or one*tutorial*(in`tutorials`

directory).