# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2019 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from pymor.core.interfaces import ImmutableInterface, abstractmethod
from pymor.parameters.base import Parametric
from pymor.vectorarrays.numpy import NumpyVectorSpace
[docs]class OperatorInterface(ImmutableInterface, Parametric):
"""Interface for |Parameter| dependent discrete operators.
An operator in pyMOR is simply a mapping which for any given
|Parameter| maps vectors from its `source` |VectorSpace|
to vectors in its `range` |VectorSpace|.
Note that there is no special distinction between functionals
and operators in pyMOR. A functional is simply an operator with
|NumpyVectorSpace| `(1)` as its `range` |VectorSpace|.
Attributes
----------
solver_options
If not `None`, a dict which can contain the following keys:
:'inverse': solver options used for
:meth:`~OperatorInterface.apply_inverse`
:'inverse_adjoint': solver options used for
:meth:`~OperatorInterface.apply_inverse_adjoint`
:'jacobian': solver options for the operators returned
by :meth:`~OperatorInterface.jacobian`
(has no effect for linear operators)
If `solver_options` is `None` or a dict entry is missing
or `None`, default options are used.
The interpretation of the given solver options is up to
the operator at hand. In general, values in `solver_options`
should either be strings (indicating a solver type) or
dicts of options, usually with an entry `'type'` which
specifies the solver type to use and further items which
configure this solver.
linear
`True` if the operator is linear.
source
The source |VectorSpace|.
range
The range |VectorSpace|.
H
The adjoint operator, i.e. ::
self.H.apply(V, mu) == self.apply_adjoint(V, mu)
for all V, mu.
"""
solver_options = None
@property
def H(self):
from pymor.operators.constructions import AdjointOperator
return AdjointOperator(self)
[docs] @abstractmethod
def apply(self, U, mu=None):
"""Apply the operator to a |VectorArray|.
Parameters
----------
U
|VectorArray| of vectors to which the operator is applied.
mu
The |Parameter| for which to evaluate the operator.
Returns
-------
|VectorArray| of the operator evaluations.
"""
pass
[docs] @abstractmethod
def apply2(self, V, U, mu=None):
"""Treat the operator as a 2-form by computing ``V.dot(self.apply(U))``.
If the operator is a linear operator given by multiplication with a matrix
M, then `apply2` is given as::
op.apply2(V, U) = V^T*M*U.
In the case of complex numbers, note that `apply2` is anti-linear in the
first variable by definition of `dot`.
Parameters
----------
V
|VectorArray| of the left arguments V.
U
|VectorArray| of the right right arguments U.
mu
The |Parameter| for which to evaluate the operator.
Returns
-------
A |NumPy array| with shape `(len(V), len(U))` containing the 2-form
evaluations.
"""
pass
[docs] @abstractmethod
def pairwise_apply2(self, V, U, mu=None):
"""Treat the operator as a 2-form by computing ``V.dot(self.apply(U))``.
Same as :meth:`OperatorInterface.apply2`, except that vectors from `V`
and `U` are applied in pairs.
Parameters
----------
V
|VectorArray| of the left arguments V.
U
|VectorArray| of the right right arguments U.
mu
The |Parameter| for which to evaluate the operator.
Returns
-------
A |NumPy array| with shape `(len(V),) == (len(U),)` containing
the 2-form evaluations.
"""
pass
[docs] @abstractmethod
def apply_adjoint(self, V, mu=None):
"""Apply the adjoint operator.
For any given linear |Operator| `op`, |Parameter| `mu` and
|VectorArrays| `U`, `V` in the :attr:`~OperatorInterface.source`
resp. :attr:`~OperatorInterface.range` we have::
op.apply_adjoint(V, mu).dot(U) == V.dot(op.apply(U, mu))
Thus, when `op` is represented by a matrix `M`, `apply_adjoint` is
given by left-multplication of (the complex conjugate of) `M` with `V`.
Parameters
----------
V
|VectorArray| of vectors to which the adjoint operator is applied.
mu
The |Parameter| for which to apply the adjoint operator.
Returns
-------
|VectorArray| of the adjoint operator evaluations.
"""
pass
[docs] @abstractmethod
def apply_inverse(self, V, mu=None, least_squares=False):
"""Apply the inverse operator.
Parameters
----------
V
|VectorArray| of vectors to which the inverse operator is applied.
mu
The |Parameter| for which to evaluate the inverse operator.
least_squares
If `True`, solve the least squares problem::
u = argmin ||op(u) - v||_2.
Since for an invertible operator the least squares solution agrees
with the result of the application of the inverse operator,
setting this option should, in general, have no effect on the result
for those operators. However, note that when no appropriate
|solver_options| are set for the operator, most implementations
will choose a least squares solver by default which may be
undesirable.
Returns
-------
|VectorArray| of the inverse operator evaluations.
Raises
------
InversionError
The operator could not be inverted.
"""
pass
[docs] @abstractmethod
def apply_inverse_adjoint(self, U, mu=None, least_squares=False):
"""Apply the inverse adjoint operator.
Parameters
----------
U
|VectorArray| of vectors to which the inverse adjoint operator is applied.
mu
The |Parameter| for which to evaluate the inverse adjoint operator.
least_squares
If `True`, solve the least squares problem::
v = argmin ||op^*(v) - u||_2.
Since for an invertible operator the least squares solution agrees
with the result of the application of the inverse operator,
setting this option should, in general, have no effect on the result
for those operators. However, note that when no appropriate
|solver_options| are set for the operator, most operator
implementations will choose a least squares solver by default which
may be undesirable.
Returns
-------
|VectorArray| of the inverse adjoint operator evaluations.
Raises
------
InversionError
The operator could not be inverted.
"""
pass
[docs] @abstractmethod
def jacobian(self, U, mu=None):
"""Return the operator's Jacobian as a new |Operator|.
Parameters
----------
U
Length 1 |VectorArray| containing the vector for which to compute
the Jacobian.
mu
The |Parameter| for which to compute the Jacobian.
Returns
-------
Linear |Operator| representing the Jacobian.
"""
pass
[docs] @abstractmethod
def d_mu(self, component, index=()):
"""Return the operator's derivative with respect to an index of a parameter component.
Parameters
----------
component
Parameter component
index
index in the parameter component
Returns
-------
New |Operator| representing the partial derivative.
"""
pass
[docs] def as_range_array(self, mu=None):
"""Return a |VectorArray| representation of the operator in its range space.
In the case of a linear operator with |NumpyVectorSpace| as
:attr:`~OperatorInterface.source`, this method returns for every |Parameter|
`mu` a |VectorArray| `V` in the operator's :attr:`~OperatorInterface.range`,
such that ::
V.lincomb(U.to_numpy()) == self.apply(U, mu)
for all |VectorArrays| `U`.
Parameters
----------
mu
The |Parameter| for which to return the |VectorArray| representation.
Returns
-------
V
The |VectorArray| defined above.
"""
assert isinstance(self.source, NumpyVectorSpace) and self.linear
raise NotImplementedError
[docs] def as_source_array(self, mu=None):
"""Return a |VectorArray| representation of the operator in its source space.
In the case of a linear operator with |NumpyVectorSpace| as
:attr:`~OperatorInterface.range`, this method returns for every |Parameter|
`mu` a |VectorArray| `V` in the operator's :attr:`~OperatorInterface.source`,
such that ::
self.range.make_array(V.dot(U).T) == self.apply(U, mu)
for all |VectorArrays| `U`.
Parameters
----------
mu
The |Parameter| for which to return the |VectorArray| representation.
Returns
-------
V
The |VectorArray| defined above.
"""
assert isinstance(self.range, NumpyVectorSpace) and self.linear
raise NotImplementedError
[docs] def as_vector(self, mu=None):
"""Return a vector representation of a linear functional or vector operator.
Depending on the operator's :attr:`~OperatorInterface.source` and
:attr:`~OperatorInterface.range`, this method is equivalent to calling
:meth:`~OperatorInterface.as_range_array` or :meth:`~OperatorInterface.as_source_array`
respectively. The resulting |VectorArray| is required to have length 1.
Parameters
----------
mu
The |Parameter| for which to return the vector representation.
Returns
-------
V
|VectorArray| of length 1 containing the vector representation.
"""
if not self.linear:
raise TypeError('This nonlinear operator does not represent a vector or linear functional.')
if self.source.is_scalar:
return self.as_range_array(mu)
elif self.range.is_scalar:
return self.as_source_array(mu)
else:
raise TypeError('This operator does not represent a vector or linear functional.')
[docs] @abstractmethod
def assemble(self, mu=None):
"""Assemble the operator for a given parameter.
The result of the method strongly depends on the given operator.
For instance, a matrix-based operator will assemble its matrix, a |LincombOperator|
will try to form the linear combination of its operators, whereas an arbitrary
operator might simply return a :class:`~pymor.operators.constructions.FixedParameterOperator`.
The only assured property of the assembled operator is that it no longer
depends on a |Parameter|.
Parameters
----------
mu
The |Parameter| for which to assemble the operator.
Returns
-------
Parameter-independent, assembled |Operator|.
"""
pass
def _assemble_lincomb(self, operators, coefficients, identity_shift=0., solver_options=None, name=None):
"""Try to assemble a linear combination of the given operators.
Returns a new |Operator| which represents the sum ::
c_1*O_1 + ... + c_N*O_N + s*I
where `O_i` are |Operators|, `c_i`, `s` scalar coefficients and `I` the identity.
This method is called in the :meth:`assemble` method of |LincombOperator| on
the first of its operators. If an assembly of the given linear combination
is possible, e.g. the linear combination of the system matrices of the
operators can be formed, then the assembled operator is returned.
Otherwise, the method returns `None` to indicate that assembly is not possible.
Parameters
----------
operators
List of |Operators| `O_i` whose linear combination is formed.
coefficients
List of the corresponding linear coefficients `c_i`.
identity_shift
The coefficient `s`.
solver_options
|solver_options| for the assembled operator.
name
Name of the assembled operator.
Returns
-------
The assembled |Operator| if assembly is possible, otherwise `None`.
"""
return None
[docs] def restricted(self, dofs):
"""Restrict the operator range to a given set of degrees of freedom.
This method returns a restricted version `restricted_op` of the
operator along with an array `source_dofs` such that for any
|VectorArray| `U` in `self.source` the following is true::
self.apply(U, mu).dofs(dofs)
== restricted_op.apply(NumpyVectorArray(U.dofs(source_dofs)), mu))
Such an operator is mainly useful for
:class:`empirical interpolation <pymor.operators.ei.EmpiricalInterpolatedOperator>`
where the evaluation of the original operator only needs to be known
for few selected degrees of freedom. If the operator has a small
stencil, only few `source_dofs` will be needed to evaluate the
restricted operator which can make its evaluation very fast
compared to evaluating the original operator.
Parameters
----------
dofs
One-dimensional |NumPy array| of degrees of freedom in the operator
:attr:`~OperatorInterface.range` to which to restrict.
Returns
-------
restricted_op
The restricted operator as defined above. The operator will have
|NumpyVectorSpace| `(len(source_dofs))` as :attr:`~OperatorInterface.source`
and |NumpyVectorSpace| `(len(dofs))` as :attr:`~OperatorInterface.range`.
source_dofs
One-dimensional |NumPy array| of source degrees of freedom as
defined above.
"""
raise NotImplementedError
[docs] @abstractmethod
def __add__(self, other):
"""Sum of two operators."""
pass
def __sub__(self, other):
return self + (- other)
[docs] @abstractmethod
def __mul__(self, other):
"""Product of operator by a scalar."""
pass
def __rmul__(self, other):
return self * other
[docs] @abstractmethod
def __matmul__(self, other):
"""Concatenation of two operators."""
pass
def __neg__(self):
return self * (-1.)