Permalink
Browse files

Merge pull request #191 from pymor/apply_inverse_adjoint

[operators] add 'apply_inverse_adjoint' to OperatorInterface
  • Loading branch information...
sdrave committed Nov 23, 2015
2 parents d97fd9e + 3c0e23a commit 541406f30bc0a1e3feb7da79e15b73a34e0685a0
@@ -150,6 +150,32 @@ def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
R.append(newton(self, V.copy(i), **options)[0])
return R
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None,
least_squares=False):
from pymor.operators.constructions import FixedParameterOperator
assembled_op = self.assemble(mu)
if assembled_op != self and not isinstance(assembled_op, FixedParameterOperator):
return assembled_op.apply_inverse_adjoint(U, ind=ind, source_product=source_product,
range_product=range_product, least_squares=least_squares)
elif source_product or range_product:
if source_product:
U = source_product.apply(U, ind=ind)
ind = None
# maybe there is a better implementation for source_product == None and range_product == None
V = self.apply_inverse_adjoint(U, mu=mu, least_squares=least_squares)
if range_product:
return range_product.apply_inverse(V)
else:
return V
else:
if not self.linear:
raise NotImplementedError
# use generic solver for the adjoint operator
from pymor.operators.constructions import AdjointOperator
options = {'inverse': self.solver_options.get('inverse_adjoint') if self.solver_options else None}
adjoint_op = AdjointOperator(self.with_(solver_options=options), with_apply_inverse=False)
return adjoint_op.apply_inverse(U, ind=ind, mu=mu, least_squares=least_squares)
def as_vector(self, mu=None):
if not self.linear:
raise TypeError('This nonlinear operator does not represent a vector or linear functional.')
@@ -351,6 +351,23 @@ def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product
else:
return PrU
def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
assert V in self.range
return V.copy(ind=ind)
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None, least_squares=False):
if source_product or range_product:
return super(IdentityOperator, self).apply_inverse_adjoint(U, ind=ind, mu=mu,
source_product=source_product,
range_product=range_product,
least_squares=least_squares)
else:
assert U in self.source
return U.copy(ind=ind)
def assemble(self, mu=None):
return self
def assemble_lincomb(self, operators, coefficients, solver_options=None, name=None):
if all(isinstance(op, IdentityOperator) for op in operators):
assert all(op.source == operators[0].source for op in operators)
@@ -535,6 +552,16 @@ def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product
else:
return ATPrU
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None, least_squares=False):
if source_product or range_product:
return super(VectorArrayOperator, self).apply_inverse_adjoint(U, ind, mu=mu,
source_product=source_product,
range_product=range_product,
least_squares=least_squares)
else:
adjoint_op = VectorArrayOperator(self._array, transposed=not self.transposed, copy=False)
return adjoint_op.apply_inverse(U, ind=ind, mu=mu, least_squares=least_squares)
def assemble_lincomb(self, operators, coefficients, solver_options=None, name=None):
transposed = operators[0].transposed
@@ -668,6 +695,12 @@ def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product
def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
return self.operator.apply_inverse(V, ind=ind, mu=self.mu, least_squares=least_squares)
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None, least_squares=False):
return self.operator.apply_inverse_adjoint(U, ind=ind, mu=self.mu,
source_product=source_product,
range_product=range_product,
least_squares=least_squares)
def jacobian(self, U, mu=None):
return self.operator.jacobian(U, mu=self.mu)
@@ -689,11 +722,18 @@ class AdjointOperator(OperatorBase):
w.r.t. which to take the adjoint.
name
If not `None`, name of the operator.
with_apply_inverse
If `True`, provide own :meth:`~pymor.operators.interfaces.OperatorInterface.apply_inverse`
and :meth:`~pymor.operator.interfaces.OperatorInterface.apply_inverse_adjoint`
implementations by calling these methods on the given `operator`.
(Is set to `False` in the default implementation of
and :meth:`~pymor.operator.interfaces.OperatorInterface.apply_inverse_adjoint`.)
"""
linear = True
def __init__(self, operator, source_product=None, range_product=None, name=None):
def __init__(self, operator, source_product=None, range_product=None, name=None,
with_apply_inverse=True):
assert isinstance(operator, OperatorInterface)
assert operator.linear
self.build_parameter_type(inherits=(operator,))
@@ -703,6 +743,7 @@ def __init__(self, operator, source_product=None, range_product=None, name=None)
self.source_product = source_product
self.range_product = range_product
self.name = name or operator.name + '_adjoint'
self.with_apply_inverse=with_apply_inverse
def apply(self, U, ind=None, mu=None):
return self.operator.apply_adjoint(U, ind=ind, mu=mu,
@@ -727,6 +768,39 @@ def apply_adjoint(self, U, ind=None, mu=None, source_product=None, range_product
return U
def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
if not self.with_apply_inverse:
return super(AdjointOperator, self).apply_inverse(V, ind=ind, mu=mu, least_squares=least_squares)
return self.operator.apply_inverse_adjoint(V, ind=ind, mu=mu,
source_product=self.source_product,
range_product=self.range_product,
least_squares=least_squares)
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None, least_squares=False):
if not self.with_apply_inverse:
return super(AdjointOperator, self).apply_inverse_adjoint(U, ind=ind, mu=mu,
source_product=source_product,
range_product=range_product,
least_squares=least_squares)
assert U in self.source
if source_product and source_product != self.range_product:
U = source_product.apply(U, ind=ind)
ind = None
if self.range_product and source_product != self.range_product:
U = self.range_product.apply_inverse(U, ind=ind)
ind = None
V = self.operator.apply_inverse(U, ind=ind, mu=mu, least_squares=least_squares)
if self.source_product and self.source_product != range_product:
V = self.source_product.apply(V)
if range_product and self.source_product != range_product:
V = range_product.apply_inverse(V)
return V
def projected(self, range_basis, source_basis, product=None, name=None):
if range_basis is not None:
if product is not None:
@@ -25,11 +25,13 @@ class OperatorInterface(ImmutableInterface, Parametric):
solver_options
If not `None`, a dict which can contain the follwing keys:
:'inverse': solver options used for
:meth:`~OperatorInterface.apply_inverse`
:'jacobian': solver options for the operators returned
by :meth:`~OperatorInterface.jacobian`
(has no effect for linear operators)
:'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.
@@ -212,6 +214,48 @@ def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
"""
pass
@abstractmethod
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None,
least_squares=False):
"""Apply the inverse adjoint operator.
Parameters
----------
U
|VectorArray| of vectors to which the inverse adjoint operator is applied.
ind
The indices of the vectors in `U` to which the inverse adjoint operator shall be
applied. (See the |VectorArray| documentation for further details.)
mu
The |Parameter| for which to evaluate the inverse adjoint operator.
source_product
See :meth:`~OperatorInterface.apply_adjoint`.
range_product
See :meth:`~OperatorInterface.apply_adjoint`.
least_squares
If `True`, solve the least squares problem::
v = argmin ||A*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 appropriate
|solver_options| are not set for the operator, most operator
implementations will choose a least squares solver by default which
may not be desirable for invertible operators.
Returns
-------
|VectorArray| of the inverse adjoint operator evaluations.
Raises
------
InversionError
The operator could not be inverted.
"""
pass
@abstractmethod
def jacobian(self, U, mu=None):
"""Return the operator's Jacobian.
@@ -256,6 +256,18 @@ def apply_inverse(self, V, ind=None, mu=None, least_squares=False):
raise InversionError(msg)
raise e
def apply_inverse_adjoint(self, U, ind=None, mu=None, source_product=None, range_product=None,
least_squares=False):
if source_product or range_product:
return super(NumpyMatrixOperator, self).apply_inverse_adjoint(U, ind=ind, mu=mu,
source_product=source_product,
range_product=range_product,
least_squares=least_squares)
else:
options = {'inverse': self.solver_options.get('inverse_adjoint') if self.solver_options else None}
adjoint_op = NumpyMatrixOperator(self._matrix.T, solver_options=options)
return adjoint_op.apply_inverse(U, ind=ind, mu=mu, least_squares=least_squares)
def projected_to_subbasis(self, dim_range=None, dim_source=None, name=None):
"""Project the operator to a subbasis.
@@ -177,9 +177,34 @@ def test_apply_inverse(operator_with_arrays):
assert U in op.source
assert len(U) == V.len_ind(ind)
VV = op.apply(U, mu=mu)
assert float_cmp_all(VV.l2_norm(), V.l2_norm(ind=ind), atol=1e-10, rtol=0.5)
assert np.all(almost_equal(VV, V, V_ind=ind, atol=1e-10, rtol=1e-3))
def test_apply_inverse_adjoint(operator_with_arrays):
op, mu, U, _ = operator_with_arrays
for ind in valid_inds(U):
try:
V = op.apply_inverse_adjoint(U, mu=mu, ind=ind)
except InversionError:
return
assert V in op.range
assert len(V) == U.len_ind(ind)
UU = op.apply_adjoint(V, mu=mu)
assert np.all(almost_equal(UU, U, V_ind=ind, atol=1e-10, rtol=1e-3))
def test_apply_inverse_adjoint_with_products(operator_with_arrays_and_products):
op, mu, U, _, sp, rp = operator_with_arrays_and_products
for ind in valid_inds(U):
try:
V = op.apply_inverse_adjoint(U, mu=mu, ind=ind, source_product=sp, range_product=rp)
except InversionError:
return
assert V in op.range
assert len(V) == U.len_ind(ind)
UU = op.apply_adjoint(V, mu=mu, source_product=sp, range_product=rp)
assert np.all(almost_equal(UU, U, V_ind=ind, atol=1e-10, rtol=1e-3))
def test_projected(operator_with_arrays):
op, mu, U, V = operator_with_arrays
op_UV = op.projected(V, U)

0 comments on commit 541406f

Please sign in to comment.