# Source code for pymor.algorithms.projection

```# This file is part of the pyMOR project (http://www.pymor.org).

import numpy as np

from pymor.algorithms.rules import RuleTable, match_class, match_generic
from pymor.core.exceptions import RuleNotMatchingError
from pymor.operators.basic import ProjectedOperator
from pymor.operators.constructions import (LincombOperator, Concatenation, ConstantOperator,
IdentityOperator)
from pymor.operators.ei import EmpiricalInterpolatedOperator, ProjectedEmpiciralInterpolatedOperator
from pymor.operators.interfaces import OperatorInterface
from pymor.operators.numpy import NumpyMatrixOperator
from pymor.vectorarrays.numpy import NumpyVectorSpace

[docs]def project(op, range_basis, source_basis, product=None):
"""Petrov-Galerkin projection of a given |Operator|.

Given an inner product `( ⋅, ⋅)`, source vectors `b_1, ..., b_N`
and range vectors `c_1, ..., c_M`, the projection `op_proj` of `op`
is defined by ::

[ op_proj(e_j) ]_i = ( c_i, op(b_j) )

for all i,j, where `e_j` denotes the j-th canonical basis vector of R^N.

In particular, if the `c_i` are orthonormal w.r.t. the given product,
then `op_proj` is the coordinate representation w.r.t. the `b_i/c_i` bases
of the restriction of `op` to `span(b_i)` concatenated with the
orthogonal projection onto `span(c_i)`.

From another point of view, if `op` is viewed as a bilinear form
(see :meth:`apply2`) and `( ⋅, ⋅ )` is the Euclidean inner
product, then `op_proj` represents the matrix of the bilinear form restricted
to `span(b_i) / span(c_i)` (w.r.t. the `b_i/c_i` bases).

How the projection is realized will depend on the given |Operator|.
While a projected |NumpyMatrixOperator| will
again be a |NumpyMatrixOperator|, only a generic
:class:`~pymor.operators.basic.ProjectedOperator` can be returned
in general. The exact algorithm is specified in :class:`ProjectRules`.

Parameters
----------
range_basis
The vectors `c_1, ..., c_M` as a |VectorArray|. If `None`, no
projection in the range space is performed.
source_basis
The vectors `b_1, ..., b_N` as a |VectorArray| or `None`. If `None`,
no restriction of the source space is performed.
product
An |Operator| representing the inner product.  If `None`, the
Euclidean inner product is chosen.

Returns
-------
The projected |Operator| `op_proj`.
"""
assert source_basis is None or source_basis in op.source
assert range_basis is None or range_basis in op.range
assert product is None or product.source == product.range == op.range

return ProjectRules(range_basis, source_basis, product).apply(op)

[docs]class ProjectRules(RuleTable):
"""|RuleTable| for the :func:`project` algorithm."""

def __init__(self, range_basis, source_basis, product):
super().__init__(use_caching=True)
self.range_basis, self.source_basis, self.product = \
range_basis, source_basis, product

@match_class(ZeroOperator)
def action_ZeroOperator(self, op):
range_basis, source_basis = self.range_basis, self.source_basis
if source_basis is not None and range_basis is not None:
from pymor.operators.numpy import NumpyMatrixOperator
return NumpyMatrixOperator(np.zeros((len(range_basis), len(source_basis))),
source_id=op.source.id, range_id=op.range.id,
name=op.name)
else:
new_source = (NumpyVectorSpace(len(source_basis), op.source.id) if source_basis is not None else
op.source)
new_range = (NumpyVectorSpace(len(range_basis), op.range.id) if range_basis is not None else
op.range)
return ZeroOperator(new_range, new_source, name=op.name)

@match_class(ConstantOperator)
def action_ConstantOperator(self, op):
range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
if range_basis is not None:
projected_value = NumpyVectorSpace.make_array(range_basis.inner(op._value, product).T, op.range.id)
else:
projected_value = op._value
if source_basis is None:
return ConstantOperator(projected_value, op.source, name=op.name)
else:
return ConstantOperator(projected_value, NumpyVectorSpace(len(source_basis), op.source.id),
name=op.name)

@match_generic(lambda op: op.linear and not op.parametric, 'linear and not parametric')
def action_apply_basis(self, op):
range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
if source_basis is None:
if range_basis is None:
return op
else:
try:
V = op.apply_adjoint(product.apply(range_basis) if product else range_basis)
except NotImplementedError:
if isinstance(op.source, NumpyVectorSpace):
from pymor.operators.numpy import NumpyMatrixOperator
return NumpyMatrixOperator(V.to_numpy(),
source_id=op.source.id,
range_id=op.range.id,
name=op.name)
else:
from pymor.operators.constructions import VectorArrayOperator
else:
if range_basis is None:
V = op.apply(source_basis)
if isinstance(op.range, NumpyVectorSpace):
from pymor.operators.numpy import NumpyMatrixOperator
return NumpyMatrixOperator(V.to_numpy().T,
source_id=op.source.id,
range_id=op.range.id,
name=op.name)
else:
from pymor.operators.constructions import VectorArrayOperator
elif product is None:
from pymor.operators.numpy import NumpyMatrixOperator
return NumpyMatrixOperator(op.apply2(range_basis, source_basis),
source_id=op.source.id,
range_id=op.range.id,
name=op.name)
else:
from pymor.operators.numpy import NumpyMatrixOperator
V = op.apply(source_basis)
return NumpyMatrixOperator(product.apply2(range_basis, V),
source_id=op.source.id,
range_id=op.range.id,
name=op.name)

@match_class(Concatenation)
def action_Concatenation(self, op):
if len(op.operators) == 1:
return self.apply(op.operators[0])

range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
last, first = op.operators[0], op.operators[-1]

if source_basis is not None and first.linear and not first.parametric:
V = first.apply(source_basis)
return type(self)(range_basis, V, product).apply(op.with_(operators=op.operators[:-1]))
elif range_basis is not None and last.linear and not last.parametric:
if product:
range_basis = product.apply(range_basis)
return type(self)(V, source_basis, None).apply(op.with_(operators=op.operators[1:]))
else:
projected_first = type(self)(None, source_basis, product=None).apply(first)
projected_last = type(self)(range_basis, None, product=product).apply(last)
return Concatenation((projected_last,) + op.operators[1:-1] + (projected_first,), name=op.name)

range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
if range_basis is not None:
if product is not None:
range_basis = product.apply(range_basis)
if op.source_product:
range_basis = op.source_product.apply_inverse(range_basis)

if source_basis is not None and op.range_product:
source_basis = op.range_product.apply(source_basis)

operator = type(self)(source_basis, range_basis, None).apply(op.operator)
range_product = op.range_product if source_basis is None else None
source_product = op.source_product if range_basis is None else None
name=op.name)

@match_class(EmpiricalInterpolatedOperator)
def action_EmpiricalInterpolatedOperator(self, op):
range_basis, source_basis, product = self.range_basis, self.source_basis, self.product
if len(op.interpolation_dofs) == 0:
return self.apply(ZeroOperator(op.range, op.source, op.name))
elif not hasattr(op, 'restricted_operator') or source_basis is None:
raise RuleNotMatchingError('Has no restricted operator or source_basis is None')
if range_basis is not None:
projected_collateral_basis = NumpyVectorSpace.make_array(op.collateral_basis.inner(range_basis,
product),
op.range.id)
else:
projected_collateral_basis = op.collateral_basis

return ProjectedEmpiciralInterpolatedOperator(op.restricted_operator, op.interpolation_matrix,
NumpyVectorSpace.make_array(source_basis.dofs(op.source_dofs)),
projected_collateral_basis, op.triangular,
op.source.id, None, op.name)

@match_class(AffineOperator)
def action_AffineOperator(self, op):
return self.apply(op.affine_shift + op.linear_part)

@match_class(LincombOperator)
def action_LincombOperator(self, op):
return self.replace_children(op).with_(solver_options=None)

@match_class(SelectionOperator)
def action_SelectionOperator(self, op):
return self.replace_children(op)

@match_class(OperatorInterface)
def action_generic_projection(self, op):
op.logger.warning('Using inefficient generic projection operator')
return ProjectedOperator(op, self.range_basis, self.source_basis, self.product)

[docs]def project_to_subbasis(op, dim_range=None, dim_source=None):
"""Project already projected |Operator| to a subbasis.

The purpose of this method is to further project an operator that has been
obtained through :meth:`project` to subbases of the original projection bases, i.e. ::

project_to_subbasis(project(op, r_basis, s_basis, prod), dim_range, dim_source)

should be the same as ::

project(op, r_basis[:dim_range], s_basis[:dim_source], prod)

For a |NumpyMatrixOperator| this amounts to extracting the upper-left
(dim_range, dim_source) corner of its matrix.

The subbasis projection algorithm is specified in :class:`ProjectToSubbasisRules`.

Parameters
----------
dim_range
Dimension of the range subbasis.
dim_source
Dimension of the source subbasis.

Returns
-------
The projected |Operator|.
"""
assert dim_source is None or (isinstance(op.source, NumpyVectorSpace) and dim_source <= op.source.dim)
assert dim_range is None or (isinstance(op.range, NumpyVectorSpace) and dim_range <= op.range.dim)

return ProjectToSubbasisRules(dim_range, dim_source).apply(op)

[docs]class ProjectToSubbasisRules(RuleTable):
"""|RuleTable| for the :func:`project_to_subbasis` algorithm."""

def __init__(self, dim_range, dim_source):
super().__init__(use_caching=True)
self.dim_range, self.dim_source = dim_range, dim_source

@match_class(LincombOperator)
def action_recurse(self, op):
return self.replace_children(op)

@match_class(NumpyMatrixOperator)
def action_NumpyMatrixOperator(self, op):
# copy instead of just slicing the matrix to ensure contiguous memory
return NumpyMatrixOperator(op.matrix[:self.dim_range, :self.dim_source].copy(),
source_id=op.source.id,
range_id=op.range.id,
solver_options=op.solver_options,
name=op.name)

@match_class(ConstantOperator)
def action_ConstantOperator(self, op):
dim_range, dim_source = self.dim_range, self.dim_source
source = op.source if dim_source is None else NumpyVectorSpace(dim_source, op.source.id)
value = op._value if dim_range is None else NumpyVectorSpace(op._value.to_numpy()[:, :dim_range], op.range.id)
return ConstantOperator(value, source, name=op.name)

@match_class(IdentityOperator)
def action_IdentityOperator(self, op):
dim_range, dim_source = self.dim_range, self.dim_source
if dim_range != dim_source:
raise RuleNotMatchingError('dim_range and dim_source must be equal.')
space = op.source if dim_source is None else NumpyVectorSpace(dim_source, op.source.id)
return IdentityOperator(space, name=op.name)

@match_class(ProjectedEmpiciralInterpolatedOperator)
def action_ProjectedEmpiciralInterpolatedOperator(self, op):
if not isinstance(op.projected_collateral_basis.space, NumpyVectorSpace):
raise NotImplementedError

restricted_operator = op.restricted_operator

old_pcb = op.projected_collateral_basis
projected_collateral_basis = NumpyVectorSpace.make_array(old_pcb.to_numpy()[:, :self.dim_range],
old_pcb.space.id)

old_sbd = op.source_basis_dofs
source_basis_dofs = NumpyVectorSpace.make_array(old_sbd.to_numpy()[:self.dim_source])

return ProjectedEmpiciralInterpolatedOperator(restricted_operator, op.interpolation_matrix,
source_basis_dofs, projected_collateral_basis, op.triangular,
op.source.id, solver_options=op.solver_options, name=op.name)

@match_class(ProjectedOperator)
def action_ProjectedOperator(self, op):
dim_range, dim_source = self.dim_range, self.dim_source
source_basis = op.source_basis if dim_source is None \
else op.source_basis[:dim_source]
range_basis = op.range_basis if dim_range is None \
else op.range_basis[:dim_range]
return ProjectedOperator(op.operator, range_basis, source_basis, product=None,
solver_options=op.solver_options)
```