# -*- coding: utf-8 -*-
# This file is part of the pyMOR project (http://www.pymor.org).
# Copyright 2013-2017 pyMOR developers and contributors. All rights reserved.
# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
from numbers import Number
import numpy as np
from pymor.algorithms import genericsolvers
from pymor.core.exceptions import InversionError
from pymor.operators.interfaces import OperatorInterface
from pymor.vectorarrays.interfaces import VectorArrayInterface
from pymor.vectorarrays.numpy import NumpyVectorSpace
[docs]class OperatorBase(OperatorInterface):
"""Base class for |Operators| providing some default implementations.
When implementing a new operator, it is usually advisable to derive
from this class.
"""
def apply2(self, V, U, mu=None):
mu = self.parse_parameter(mu)
assert isinstance(V, VectorArrayInterface)
assert isinstance(U, VectorArrayInterface)
AU = self.apply(U, mu=mu)
return V.dot(AU)
def pairwise_apply2(self, V, U, mu=None):
mu = self.parse_parameter(mu)
assert isinstance(V, VectorArrayInterface)
assert isinstance(U, VectorArrayInterface)
assert len(U) == len(V)
AU = self.apply(U, mu=mu)
return V.pairwise_dot(AU)
def jacobian(self, U, mu=None):
if self.linear:
if self.parametric:
return self.assemble(mu)
else:
return self
else:
raise NotImplementedError
def assemble(self, mu=None):
if self.parametric:
from pymor.operators.constructions import FixedParameterOperator
return FixedParameterOperator(self, mu=mu, name=self.name + '_assembled')
else:
return self
def __sub__(self, other):
if isinstance(other, Number):
assert other == 0.
return self
from pymor.operators.constructions import LincombOperator
return LincombOperator([self, other], [1, -1])
def __add__(self, other):
if isinstance(other, Number):
assert other == 0.
return self
from pymor.operators.constructions import LincombOperator
return LincombOperator([self, other], [1, 1])
__radd__ = __add__
def __mul__(self, other):
assert isinstance(other, Number)
from pymor.operators.constructions import LincombOperator
return LincombOperator([self], [other])
def __str__(self):
return '{}: R^{} --> R^{} (parameter type: {}, class: {})'.format(
self.name, self.source.dim, self.range.dim, self.parameter_type,
self.__class__.__name__)
def apply_transpose(self, V, mu=None):
if self.linear:
raise NotImplementedError
else:
raise ValueError('Trying to apply transpose of nonlinear operator.')
def apply_inverse(self, V, mu=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(V, least_squares=least_squares)
elif self.linear:
options = self.solver_options.get('inverse') if self.solver_options else None
return genericsolvers.apply_inverse(assembled_op, V, options=options, least_squares=least_squares)
else:
from pymor.algorithms.newton import newton
from pymor.core.exceptions import NewtonError
options = self.solver_options.get('inverse') if self.solver_options else None
if options:
if isinstance(options, str):
assert options == 'newton'
options = {}
else:
assert options['type'] == 'newton'
options = options.copy()
options.pop('type')
else:
options = {}
options['least_squares'] = least_squares
R = V.empty(reserve=len(V))
for i in range(len(V)):
try:
R.append(newton(self, V[i], mu=mu, **options)[0])
except NewtonError as e:
raise InversionError(e)
return R
def apply_inverse_transpose(self, U, mu=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_transpose(U, least_squares=least_squares)
else:
if not self.linear:
raise NotImplementedError
# use generic solver for the transpose operator
from pymor.operators.constructions import AdjointOperator
options = {'inverse': self.solver_options.get('inverse_transpose') if self.solver_options else None}
transpose_op = AdjointOperator(self, with_apply_inverse=False, solver_options=options)
return transpose_op.apply_inverse(U, mu=mu, least_squares=least_squares)
def as_range_array(self, mu=None):
assert isinstance(self.source, NumpyVectorSpace)
return self.apply(self.source.make_array(np.eye(self.source.dim)), mu=mu)
def as_source_array(self, mu=None):
assert isinstance(self.range, NumpyVectorSpace)
return self.apply_transpose(self.range.make_array(np.eye(self.range.dim)), mu=mu)
[docs]class ProjectedOperator(OperatorBase):
"""Generic |Operator| representing the projection of an |Operator| to a subspace.
This operator is implemented as the concatenation of the linear combination with
`source_basis`, application of the original `operator` and projection onto
`range_basis`. As such, this operator can be used to obtain a reduced basis
projection of any given |Operator|. However, no offline/online decomposition is
performed, so this operator is mainly useful for testing before implementing
offline/online decomposition for a specific application.
This operator is instantiated in :func:`pymor.algorithms.projection.project`
as a default implementation for parametric or nonlinear operators.
Parameters
----------
operator
The |Operator| to project.
source_basis
See :func:`pymor.algorithms.projection.project`.
range_basis
See :func:`pymor.algorithms.projection.project`.
product
See :func:`pymor.algorithms.projection.project`.
"""
linear = False
def __init__(self, operator, range_basis, source_basis, product=None, solver_options=None):
assert isinstance(operator, OperatorInterface)
assert source_basis is None or source_basis in operator.source
assert range_basis is None or range_basis in operator.range
assert product is None \
or (isinstance(product, OperatorInterface)
and range_basis is not None
and operator.range == product.source
and product.range == product.source)
self.build_parameter_type(operator)
self.source = NumpyVectorSpace(len(source_basis), operator.source.id) if source_basis is not None else operator.source
self.range = NumpyVectorSpace(len(range_basis), operator.range.id) if range_basis is not None else operator.range
self.solver_options = solver_options
self.name = operator.name
self.operator = operator
self.source_basis = source_basis.copy() if source_basis is not None else None
self.range_basis = range_basis.copy() if range_basis is not None else None
self.linear = operator.linear
self.product = product
@property
def T(self):
if self.product:
return super().T
else:
options = {'inverse': self.solver_options.get('inverse_transpose'),
'inverse_transpose': self.solver_options.get('inverse')} if self.solver_options else None
return ProjectedOperator(self.operator.T, self.source_basis, self.range_basis, solver_options=options)
def apply(self, U, mu=None):
mu = self.parse_parameter(mu)
if self.source_basis is None:
if self.range_basis is None:
return self.operator.apply(U, mu=mu)
elif self.product is None:
return self.range.make_array(self.operator.apply2(self.range_basis, U, mu=mu).T)
else:
V = self.operator.apply(U, mu=mu)
return self.range.make_array(self.product.apply2(V, self.range_basis))
else:
UU = self.source_basis.lincomb(U.data)
if self.range_basis is None:
return self.operator.apply(UU, mu=mu)
elif self.product is None:
return self.range.make_array(self.operator.apply2(self.range_basis, UU, mu=mu).T)
else:
V = self.operator.apply(UU, mu=mu)
return self.range.make_array(self.product.apply2(V, self.range_basis))
def jacobian(self, U, mu=None):
if self.linear:
return self.assemble(mu)
assert len(U) == 1
mu = self.parse_parameter(mu)
if self.source_basis is None:
J = self.operator.jacobian(U, mu=mu)
else:
J = self.operator.jacobian(self.source_basis.lincomb(U.data), mu=mu)
from pymor.algorithms.projection import project
pop = project(J, range_basis=self.range_basis, source_basis=self.source_basis,
product=self.product, name=self.name + '_jacobian')
if self.solver_options:
options = self.solver_options.get('jacobian')
if options:
pop = pop.with_(solver_options=options)
return pop
def assemble(self, mu=None):
op = self.operator.assemble(mu=mu)
if op == self.operator: # avoid infinite recursion in apply_inverse default impl
return self
from pymor.algorithms.projection import project
pop = project(op, range_basis=self.range_basis, source_basis=self.source_basis,
product=self.product, name=self.name + '_assembled')
if self.solver_options:
pop = pop.with_(solver_options=self.solver_options)
return pop