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 N-dimensional
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 no-side-effect 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)
In-place 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, element-wise 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 element-wise 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:
Composite operator | Code | Mathematical representation |
---|---|---|
AdditionOperator |
A + B |
\(x \mapsto A(x) + B(x)\) |
MultiplicationOperator |
MultiplicationOperatoor([A, B]) or A * B \(^1\) |
\(x \mapsto A(x) \odot B(x)\) |
CompositionOperator |
CompositionOperator([A, B]) or A * B \(^2\) |
\(x \mapsto A(B(x))\) |
BlockColumnOperator |
\(\begin{bmatrix} A \\ B \end{bmatrix}\) | |
BlockDiagonalOperator |
\(\begin{bmatrix} A & 0 \\ 0 & B \end{bmatrix}\) | |
BlockRowOperator |
\(\begin{bmatrix} A & B \end{bmatrix}\) |
\(^1\) if A or B are not linear, \(^2\) if A and B are linear.
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 matrix-vector 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 single-argument
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 fixed-size 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, inverse-conjugate, inverse-transpose and
inverse-adjoint.
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 inverse-conjugate |
'IT' |
operator’s inverse-transpose |
'IH' |
operator’s inverse-adjoint |
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 adjoint |
'I' |
reference operator’s inverse |
'1' |
identity operator |
function |
callback function |
For a non-commutative operation (such as composition), the rules attached to the left-hand side operator applying on a right-hand-side operator are combined to the rules of the right-hand side operators applying on a left-hand-side 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
left-hand-side 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
non-commutative 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
in-place operator: operator whose direct
method can handle input
and output arguments pointing to the same memory location.
out-of-place operator: operator whose direct
method cannot.
Such a property is set by the inplace
decorator. For instance, the
operator DiagonalOperator
is an in-place operator and the following
computation
>>> d = DiagonalOperator([1, 2])
>>> x = np.array([0, 1])
>>> d(x, x)
array([0, 2])
is done in-place and is equivalent to:
>>> x = np.array([0, 1])
>>> x *= [1, 2]
An in-place operator also has to handle out-of-place operations. In the following example, the operation is performed out-of-place to satisfy the no-side-effect 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 in-place property is useful to avoid intermediate variables. To minimise the memory-cache 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 in-place operator
IN
by an out-of-place operator OUT
. The out-of-place 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 in-place on y
.
Concerning the in-place 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 out-of-place
because its output is the x
variable.
Another possibility to avoid temporaries in the addition or element-wise
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 non-zero 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 in-place.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
in-place. 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 in-place. This example shows the interest of the flag
update_output
when only few rows are non-zero: 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 on-going work consisting in migrating the most generic ones
into the pyoperators package.