The Operator
class and its subclasses are function factories. It means that in order to perform computations using an operator class, we first need to create an instance of it:
A = FFTOperator(1024)
and then use the instance as a function. The object A
is a Python callable and takes numpy ’s Ndimensional array object as inputs and outputs:
>>> pyoperators.memory.verbose = True
>>> x = np.arange(1024, dtype=complex)
>>> y = A(x)
Info: Allocating (1024,) complex128 = 0.015625 MiB in FFTOperator.
When using the operator in such a way, a nosideeffect policy is enforced: it is garanteed that the input array will not be modified and as a consequence, a new buffer will be allocated every time the operator is called.
Similarly to ufuncs, such behaviour can be changed by providing the output array:
>>> pyoperators.memory.verbose = True
>>> y = np.empty(1024, dtype=complex)
>>> A(x, out=y) # or simply A(x, y)
Inplace operations can be performed in a similar way by providing the input array as the output.
In this documentation, the same term operator will be used for the Operator
class, subclasses and instances, though the context will disambiguate the meaning.
Operators are easy to manipulate: they can be multiplied by a scalar, combined by addition, elementwise multiplication, composition or partition. The synthax uses the usual arithmetic signs, but be careful that the *
sign depends on the operands’ linear
flag: if both operands are linear, it stands for composition otherwise it stands for elementwise multiplication.
>>> A = DiagonalOperator([1., 2.])
>>> B = DiagonalOperator([1., 1.])
>>> A.flags.linear, B.flags.linear
(True, True)
>>> x = [2, 2]
>>> np.allclose((3 * A)(x), 3 * A(x))
True
>>> np.allclose((A + B)(x), A(x) + B(x))
True
>>> np.allclose(A(B)(x), A(B(x)))
True
>>> np.allclose((A * B)(x), A(B(x))) # because A and B are linear
True
>>> np.allclose(MultiplicationOperator([A, B])(x), A(x) * B(x))
True
Unless an algebraic simplification is possible, the resulting operator is a composite operator:
addition  AdditionOperator 
\( A(x) + B(x) \) 
elementwise multiplication  MultiplicationOperator 
\( A(x) \times B(x) \) 
composition  CompositionOperator 
\( A(B(x)) \) 
partition  BlockColumnOperator 
\( \begin{bmatrix} A \\ B \end{bmatrix} \) 
BlockDiagonalOperator 
\( \begin{bmatrix} A & 0 \\ 0 & B \end{bmatrix} \) 

BlockRowOperator 
\( \begin{bmatrix} A & B \end{bmatrix} \) 
The conjugate, transpose, adjoint and inverse operators, when defined, can be obtained by using the following attributes:
conjugate  .C 
transpose  .T 
adjoint  .H 
inverse  .I 
For example:
>>> o = Operator()
>>> o.T.C is o.H
True
>>> o = Operator(flags='real,symmetric')
>>> o.T is o.H is o
True
Operators can be created in two ways. The first one is to define a direct function which will replace the usual matrixvector operation and to instantiate the Operator
class by passing the function as the first argument:
>>> def f(x, out):
... out[...] = 2 * x
>>> P = Operator(f)
Transforming a singleargument ufunc into an Operator
instance is straightforward:
>>> sqrt = Operator(np.sqrt)
The alternative way is more flexible, it consists in subclassing the Operator
class:
>>> class MyOperator(Operator):
... def __init__(self, coef):
... self.coef = coef
... Operator.__init__(self)
... def direct(self, x, out):
... out[...] = self.coef * x
...
>>> P = MyOperator(2)
The operator’s data type is a data type object (an instance of numpy.dtype class). It can be set using the dtype
keyword.
>>> a = Operator(dtype=np.float32)
>>> a.dtype
dtype('float32')
It is used to determine the output data type, as being the common type of the operator and input data types, following the standard coercion rules. For example:
>>> a = DiagonalOperator([1, 2])
>>> a.dtype
dtype('int64')
>>> a([1, 1j]).dtype
dtype('complex128')
>>> a(np.array([2, 2], np.uint8)).dtype
dtype('int64')
The operator’s data type can be equal to None
(the default), in which case the output data type is always the input data type.
An operator P can have the following algebraic properties, accessible with the flags
attribute:
linear  
real  conj P = P 
symmetric  P ^{T} = P 
hermitian  P ^{H} = P 
idempotent  P P = P 
involutary  P P = I 
orthogonal  P ^{T} P = I 
unitary  P ^{H} P = I 
They can be set using the flags
keyword:
>>> def f(x, out):
... out[...] = x
>>> P = Operator(direct=f, flags='symmetric,involutary')
or by using operator decorators:
>>> @pyoperators.flags.symmetric
... @pyoperators.flags.involutary
... class MyOperator(Operator):
... def direct(self, x, out):
... out[...] = x
...
>>> P = MyOperator()
These flags are used to simplify expressions:
>>> C = Operator(flags='idempotent')
>>> C * C is C
True
>>> D = Operator(flags='involutary')
>>> D * D
IdentityOperator()
Other properties are described by the flags
attribute:
square  input and output have the same shape 
separable  can be separated if applied by or applying on a block operator 
inplace  handle inplace operations (more info) 
update_output  can update the operator output (more info) 
alignment_input  only handle aligned input 
alignment_output  only handle aligned output 
contiguous_input  only handle contiguous input 
contiguous_output  only handle contiguous output 
shape_input  explicit, implicit or unconstrained (more info) 
shape_output  explicit, implicit or unconstrained (more info) 
The Operator
class has a flexible way to deal with inputs and outputs. Its output shape can be:
explicit: the operator only returns arrays of a specified shape.
implicit: the output shape is derived from the input shape through the method reshapein
.
unconstrained: the output shape does not depend on the input shape.
Likewise, its input shape can be:
explicit: the operator only applies over arrays of a specified shape.
implicit: the operator can only apply over arrays whose shape is derived from the output shape through the method reshapeout
.
unconstrained: any input shape can be handled by the operator.
The operator’s flags record this information in the shape_input
and shape_output
attributes. If the output (input) shape is
shapeout
(in
) is equal to that tuple. The reshapein
(out
) method should not be overridden.shapeout
(in
) is None
and the method reshapein
(out
) applied over the input (output) shape must return the output (input) shape.shapeout
(in
) is None
. The reshapein
(out
) should not be overridden.The explicit shapes are specified using the shapein
and shapeout
keywords. The following example appends a zero to a fixedsize input.
>>> class Op(Operator):
... def __init__(self, value):
... Operator.__init__(self, shapein=3, shapeout=4)
... self.value = value
... def direct(self, input, output):
... output[0:3] = input
... output[3] = self.value
...
>>> op = Op(0.)
>>> op.shapein, op.shapeout
((3,), (4,))
>>> op([3., 2. ,1.])
array([3., 2., 1., 0.])
This case does not happen because the input shape can be obtained from the explicit output shape and the reshapeout
method.
Such operator is obtained by only setting the shapeout
keyword. The following example returns an array of size 2 whose elements are the sum and product of the input elements.
>>> class Op(Operator):
... def __init__(self):
... Operator.__init__(self, shapeout=2)
... def direct(self, input, output):
... output[0] = np.sum(input)
... output[1] = np.product(input)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, (2,))
>>> op([1., 2., 3., 4.])
array([ 10., 24.])
This case does not happen because the output shape can be obtained from the explicit input shape and the reshapein
method.
This operator defines both reshapein
and reshapeout
methods. Square
decorated operators are examples of this kind of operators: their reshapein
and reshapeout
methods simply are lambda x:x
by default. The interest for an implicit output shape operator of also handling implicit input shapes arises when a reverse operator such as .T
or .H
is defined because in which the reshapein
and reshapeout
methods are swapped. The operator in the following example adds a zero to its input.
>>> @pyoperators.flags.linear
... class Op(Operator):
... def direct(self, input, output):
... output[:1] = input
... output[1] = 0
... def transpose(self, input, output):
... output[...] = input[:1]
... def reshapein(self, shapein):
... return (shapein[0]+1,)
... def reshapeout(self, shapeout):
... return (shapeout[0]1,)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op([1,2,3])
array([1, 2, 3, 0])
>>> op.T([1, 2, 3, 0])
array([1, 2, 3])
>>> np.array_equal(op.T.todense(4), op.todense(3).T)
True
This operator can handle arbitrary inputs and its output dimensions are derived from those of the input. To obtain such operator, one simply has to define a reshapein
method, which takes the input shape as input and returns the output shape. The following example stretches the input by a factor 2 along the first dimension.
>>> class Op(Operator):
... def direct(self, input, output):
... output[::2,...] = input
... output[1::2,...] = input
... def reshapein(self, shapein):
... return (2*shapein[0],) + shapein[1:]
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op.reshapein((2,3))
(4,3)
>>> op([[1, 2, 3],
... [2, 3, 4]])
array([[1, 2, 3],
[1, 2, 3],
[2, 3, 4],
[2, 3, 4]])
Such operator is obtained by only setting the shapein
keyword. The following example fills the output with an arithmetic progression whose coefficients are given by the two elements of the input.
>>> class Op(Operator):
... def __init__(self):
... Operator.__init__(self, shapein=2)
... def direct(self, input, output):
... output[...] = input[0] + input[1] * np.arange(output.size).reshape(output.shape)
...
>>> op = Op()
>>> op.shapein, op.shapeout
((2,), None)
>>> y = np.empty((2,2), int)
>>> op([1, 2], out=y)
array([[1, 3],
[5, 7]])
This kind of operator is of limited interest and is shown here for completeness. Since there is no way to compute the output shape from the input, the output argument must be provided. The operator in this example returns the input from which the last element has been removed.
>>> class Op(Operator):
... def direct(self, input, output):
... output[...] = input[:1]
... def reshapeout(self, shapeout):
... return (shapeout[0] + 1,)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> op([1,2,3])
ValueError: The output shape of an implicit input shape and unconstrained output shape operator cannot be inferred.
>>> y = np.empty(4)
>>> op([1.,2.,3.,4.,5.], out=y)
array([ 1., 2., 3., 4.])
>>> op(np.arange(4.), out=y)
ValueError: The input has an invalid shape '(4,)'. Expected shape is '(5,)'.
This kind of operator is the default: this behaviour is obtained by not setting the shapein
and shapeout
keywords and not defining the reshapein
and reshapeout
methods. ConstantOperator
is an example of such an operator. If the output is not provided as argument, it is assumed that the output shape is that of the input. The following example fills the output with the value of the sum of the elements of the input.
>>> class Op(Operator):
... def __init__(self):
... Operator.__init__(self)
... def direct(self, input, output):
... output[...] = np.sum(input)
...
>>> op = Op()
>>> op.shapein, op.shapeout
(None, None)
>>> y = np.empty((2,2), int)
>>> op([2, 1], out=y)
array([[3, 3],
[3, 3]])
It is possible to validate the shapes of the input or output of an operator by defining the validatein
or validateout
methods, which take the shapes of the input or output and raise a ValueError
exception to signal a validation failure. There are two advantages of using these methods instead of performing the validation in the direct
method:
The following example checks the number of axes and the dimension along the axes.
>>> class Op(Operator):
... def direct(self, input, output):
... output[...] = 1
... def validatein(self, shapein):
... if len(shapein) != 2:
... raise ValueError('The number of axes of the input is not 2.')
... if shapein[0] % 16 != 0 or shapein[1] % 16 != 0:
... raise ValueError('The input dimensions are not a multiple of 16.')
...
>>> op = Op()
>>> op(np.ones(10))
ValueError: The number of axes of the input is not 2.
>>> op(np.ones(16,17))
ValueError: The input dimensions are not a multiple of 16.
An operator propagates the class of the input array, if the latter subclasses numpy.ndarray
:
>>> class ndarray2(np.ndarray):
... pass
>>> type(I(ndarray2((2,2))))
__main__.ndarray2
It is possible to change the subclass too, by using the classout
keyword:
>>> class ndarray3(np.ndarray):
... pass
>>> I2 = IdentityOperator(classout=ndarray3)
>>> type(I2(ndarray((2,2))))
__main__.ndarray3
>>> type(I2(ndarray2((2,2))))
__main__.ndarray3
An operator also propagates attributes:
>>> class ndarray2(np.ndarray):
... pass
...
>>> x = ndarray2((2,2))
>>> x.foo = 'bar'
>>> I(x).foo
'bar'
It is also possible to add or change an attribute, by setting the attrout
keyword to a dictionary whose keys are the attribute names and whose values are the attribute values:
>>> I2 = IdentityOperator(attrout={'foo':'new_bar'})
>>> I2(x).foo
'new_bar'
>>> x.foo
'bar'
>>> I2(np.ones((2,2))).foo
'new_bar'
More flexibility is possible by passing a function, instead of a dictionary, to the attrout
keyword. This function expects a dictionary whose keys are the attribute names and whose values are the attribute values (it is the output’s __dict__
):
>>> class I2Operator(IdentityOperator):
... def __init__(self):
... IdentityOperator.__init__(self, attrout=self.add_history)
... def add_history(self, attr):
... from time import ctime
... if 'history' not in attr:
... attr['history'] = []
... attr['history'] += [ctime() + ' : ' + self.__class__.__name__]
...
>>> I2 = I2Operator()
>>> I2(I2(0)).history
The set of operators that can be obtained by using any number of times the '.C'
, '.T'
, '.H'
and '.I'
attributes over a given operator is finite. The set is made of the operator itself, its conjugate, transpose, adjoint, inverse, inverseconjugate, inversetranspose and inverseadjoint.
It is possible to define such operators by using a unary rule of the form ‘subject’ is ‘predicate’. The subject being
'C' 
operator’s conjugate 
'T' 
operator’s transpose 
'H' 
operator’s adjoint 
'I' 
operator’s inverse 
'IC' 
operator’s inverseconjugate 
'IT' 
operator’s inversetranspose 
'IH' 
operator’s inverseadjoint 
and the predicate being
'.' 
the operator itself  
function 
callback function 
For instance, the rule 'C'
→ '.'
means that the operator’s conjugate is the operator itself. This is how the real
decorator is translated internally. Likewise, the involutary
decorator is represented by the rule 'I'
→ '.'
.
Unary rules are attached to an operator by using the set_rule
method:
>>> @pyoperators.flags.linear
... class Op1(Operator):
... def __init__(self, **keywords):
... Operator.__init__(self, **keywords)
... self.set_rule('T', lambda s: ReverseOperatorFactory(Op2, s))
... def direct(self, input, output):
... output[:1] = input
... output[1] = 0
... def reshapein(self, shapein):
... return (shapein[0]+1,)
...
>>> @pyoperators.flags.linear
... class Op2(Operator):
... def __init__(self, **keywords):
... Operator.__init__(self, **keywords)
... self.set_rule('T', lambda s: ReverseOperatorFactory(Op1, s))
... def direct(self, input, output):
... output[...] = input[:1]
... def reshapein(self, shapein):
... return (shapein[0]1,)
...
>>> op1 = Op1()
>>> op2 = op1.T
>>> type(op2)
__main__.Op2
>>> array_equal(op1.todense(4), op2.todense(5).T)
True
The ReverseOperatorFactory
is a convenience factory which returns an operator of the type its first argument and by swapping the attributes attrin
(out
), classin
(out
), shapein
(out
), reshapein
(out
), toshapein
(out
), validatein
(out
) from the second argument.
PyOperators has the ability to reduce arithmetic expressions involving operators. This is achieved by means of binary rules which are applied to pairs of operators. Such a rule is specific to an operation (addition, multiplication or composition) and has the form ‘subject1 \& subject2 → predicate’, in which subject1 and subject2 are properties that will be matched against the pair of actual operators. One of them must be ‘.’, and stands for the reference operator to which the rule is attached. For instance and for a given operator P
, the subject ('.', '.')
will match the pair (P
, P
), the subject ('T', '.')
the pair (P.T
, P
) and ('.', 'H')
the pair (P
, P.H
).
In addition to the subjects that can be used for unary rules, it is possible to specify an Operator
subclass for subclass matching. For example the subject (‘.’, Operator) will match any Operator instance on the right of the reference operator.
When the pair of operators is matched by the subject, it is replaced by the predicate:
'.' 
reference operator 
'C' 
reference operator’s conjugate 
'T' 
reference operator’s transpose 
'H' 
reference operator’s ajoint 
'I' 
reference operator’s inverse 
'1' 
identity operator 
function 
callback function 
For a noncommutative operation (such as composition), the rules attached to the lefthand side operator applying on a righthandside operator are combined to the rules of the righthand side operators applying on a lefthandside operator:
Rules attached to an operator are also prioritised:
Let’s examine some examples for the composition operation:
('T', '.')
→ '1'
: if the expression P.T*P
is matched, it will be replaced by the identity. This is how the orthogonal property is translated into a rule.('.', 'I')
→ '1'
: if the expression P*P.I
is matched, it will be replaced by the identity.('.', Operator)
→ myfunc1
: every time P
appears as the lefthandside operator of a pair, the callback function myfunc
will be called and the pair of operators will be replaced by the result of this function, unless it is None
, in which case the pair of operators is unchanged.The rules attached to an operator are stored in the rules
attribute. It is a dictionary whose keys can be AdditionOperator
, CompositionOperator
or MultiplicationOperator
. For the noncommutative composition operation, rules are further split according to whether the operator is matched on the left or right hand side of the pair.
Rules are added by using the set_rule
method, which has 3 arguments:
AdditionOperator
, CompositionOperator
, MultiplicationOperator
)The predicate can be a function, in which case it must be a static method with two arguments for the pair of operators that matches the subject. Rules, such as the ones set in the superclass, can be deleted with the del_rule
method. In the following example, we construct an operator Power
that raises its input to a specified power. We add a rule to reduce the composition of two instances of the Power
class.
>>> @pyoperators.flags.square
... class Power(Operator):
... def __init__(self, exponent):
... Operator.__init__(self)
... self.exponent = exponent
... self.set_rule(('.', Power), self.rule_power, CompositionOperator)
... def direct(input, output):
... output[...] = input ** self.exponent
... @staticmethod
... def rule_power(p1, p2):
... return Power(p1.exponent * p2.exponent)
...
>>> p1 = Power(2)
>>> p2 = Power(3)
>>> p = p1(p2)
>>> type(p)
__main__.Power
>>> p.exponent
6
inplace operator: operator whose direct
method can handle input and output arguments pointing to the same memory location.
outofplace operator: operator whose direct
method cannot.
Such a property is set by the inplace
decorator. For instance, the operator DiagonalOperator
is an inplace operator and the following computation
>>> d = DiagonalOperator([1, 2])
>>> x = np.array([0, 1])
>>> d(x, x)
array([0, 2])
is done inplace and is equivalent to:
>>> x = np.array([0, 1])
>>> x *= [1, 2]
An inplace operator also has to handle outofplace operations. In the following example, the operation is performed outofplace to satisfy the nosideeffect policy:
>>> d = DiagonalOperator([1, 2])
>>> x = np.array([0, 1])
>>> y = d(x)
>>> y
array([0, 2])
which is equivalent to:
>>> x = np.array([0, 1])
>>> y = x * [1, 2]
The inplace property is useful to avoid intermediate variables. To minimise the memorycache transfers during a composition, an algorithm has been put in place to determine the intermediate variables to be extracted from the memory manager, by maximising the temporal locality (i.e. by reusing as much as possible the same memory area). This algorithm depends on:
As an example, let’s consider the composition of an inplace operator IN
by an outofplace operator OUT
. The outofplace composition
>>> (IN * OUT)(x, out=y)
requires an intermediate variable only if the size of OUT
‘s output is larger than that of the y
variable. Otherwise, y
’s buffer is used as output for the OUT
operator and the application of the operator IN
is performed inplace on y
.
Concerning the inplace composition:
>>> (IN * OUT)(x, out=x)
it is not possible to avoid the use of a temporary variable, since a buffer different from the x
variable is required for the OUT
operator. The application of the operator IN
is performed outofplace because its output is the x
variable.
Another possibility to avoid temporaries in the addition or elementwise multiplication of operators is to give an operator the possibility to directly update the running output argument in the direct
method. It is particularly interesting for sums of matrices with few nonzero rows, since it gives a mechanism to skip the zero ones.
Let’s first see the default behaviour, in which the operands do not update their output argument:
>>> (o1 + o2)(x, out=y)
The following steps are performed:
y
is used as o1
’s output,o2
’s output,y
buffer.If we now assume that the operator o2
can update its output argument, the intermediate variable is not required anymore:
y
is used as o1
’s outputo2
’s output argument to be updated inplace.One can enable such property by adding the flag update_output
and the keyword operation
to the operator’s direct method.
>>> import operator
>>> import pyoperators
>>> from pyoperators import Operator, operation_assignment
>>> pyoperators.memory.verbose=True
>>> @pyoperators.flags.linear
... @pyoperators.flags.update_output
... class P(Operator):
... def __init__(self, index):
... self.index = index
... Operator.__init__(self, shapein=1, shapeout=11)
... def direct(self, input, output, operation=operation_assignment):
... if operation is operation_assignment:
... output[...] = 0
... elif operation is not operator.iadd:
... raise NotImplementedError()
... output[self.index] += input
...
>>> Q = sum(P(i) for i in [0, 5, 10])
>>> y = np.empty(11)
>>> Q([2.], out=y)
array([ 2., 0., 0., 0., 0., 2., 0., 0., 0., 0., 2.])
No intermediate variable is needed and all operations are performed inplace. The AdditionOperator
Q
calls the direct
method of its first operand P
with the operation
keyword set to the function operation_assignment
, which means that P
‘s output should be assigned to the output argument. In the subsequent calls, the operation
keyword is set to the iadd
function from Python’s standard operator
module (If Q
were a MultiplicationOperator
, it would be imul
), which results in adding the operand’s output to the output argument inplace. This example shows the interest of the flag update_output
when only few rows are nonzero: the output only need to be updated for these rows, which avoids thrashing the cache with zeros everytime an operand is called.
Currently, there are two ways to partition operators, depending on whether the partitioning is done in an already existing dimension or not.
Stack partition:  outputs are stacked along a new dimension 
Chunk partition:  outputs are concatenated along an existing dimension 
A stack partition operator is an instance of BlockRowOperator
initialised with the new_axisin
keyword, BlockDiagonalOperator
with the keywords new_axisin
or new_axisout
or BlockColumnOperator
initialised with the new_axisout
keyword. There is a strong constraint on the input and output shapes of the blocks, since they must be the same for all blocks. The partition along the specified axis is always explicit and is equal to a tuple of ones with as many elements as the number of blocks. In the following example is shown a stack partition block column operator:
>>> C = BlockColumnOperator([I, 2*I, 3*I], new_axisout=0)
>>> C(np.ones(2))
array([[ 1., 1.],
[ 2., 2.],
[ 3., 3.]])
The outputs of the blocks are stacked along the dimension specified by new_axisout
keyword, i.e the first one. The following example shows a stack partition block diagonal operator:
>>> D = BlockDiagonalOperator([I, 2*I, 3*I], new_axisin=1)
>>> x = np.arange(2*3).reshape((2,3))
>>> x
array([[0, 1, 2],
[3, 4, 5]])
>>> D(x)
array([[ 0, 2, 6],
[ 3, 8, 15]])
It can be seen that the first block I
is applied over x[:,0]
, the second block 2*I
over x[:,1]
, and the third block 3*I
over x[:,2]
.
A chunk partition operator can be obtained as an instance of BlockRowOperator
initialised with the axisin
keyword, BlockDiagonalOperator
with the keywords axisin
or axisout
or BlockColumnOperator
initialised with the axisout
keyword. There is a lesser constraint on the shapes of the partitioned input and output of the blocks: all dimensions must be same except along the partitioned dimension(s). The partition along the specified axis can be implicit or explicit, which is reflected in the keywords partitionin
and partitionout
. The following example shows a block diagonal chunk partition operator with an explicit partition
>>> D = BlockDiagonalOperator([I, 2*I, 3*I], axisin=1, partitionin=(2,3,2))
>>> D(np.ones(7))
array([ 1., 1., 2., 2., 2., 3., 3.])
It can be seen that the first block I
is applied over the first two elements, the second block 2*I
over the three following elements, and the third block 3*I
over the last two elements. Block row and block column chunk partition operators can handle implicit partitions:
>>> C = BlockColumnOperator([I, 2*I, 3*I], axisout=1)
>>> C(np.ones((2,2)))
array([[ 1., 1., 2., 2., 3., 3.],
[ 1., 1., 2., 2., 3., 3.]])
Likewise, the outputs of the blocks are concatenated along the dimension specified by the axisout
keyword, i.e the last one in this case.
More operators can be found in the project PySimulators, projection, discrete differences, compression, downsampling, MPI operators, etc.
There is an ongoing work consisting in migrating the most generic ones into the pyoperators package.