Source code for pymor.operators.numpy

# 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)

"""This module provides the following |NumPy| based |Operators|:

  - |NumpyMatrixOperator| wraps a 2D |NumPy array| as an |Operator|.
  - |NumpyMatrixBasedOperator| should be used as base class for all |Operators|
    which assemble into a |NumpyMatrixOperator|.
  - |NumpyGenericOperator| wraps an arbitrary Python function between
    |NumPy arrays| as an |Operator|.
"""

from functools import reduce

import numpy as np
from scipy.io import mmwrite, savemat
from scipy.linalg import solve
import scipy.sparse
from scipy.sparse import issparse

from pymor.core.config import config
from pymor.core.defaults import defaults
from pymor.core.exceptions import InversionError
from pymor.core.interfaces import abstractmethod
from pymor.core.logger import getLogger
from pymor.operators.basic import OperatorBase
from pymor.vectorarrays.numpy import NumpyVectorSpace


[docs]class NumpyGenericOperator(OperatorBase): """Wraps an arbitrary Python function between |NumPy arrays| as an |Operator|. Parameters ---------- mapping The function to wrap. If `parameter_type` is `None`, the function is of the form `mapping(U)` and is expected to be vectorized. In particular:: mapping(U).shape == U.shape[:-1] + (dim_range,). If `parameter_type` is not `None`, the function has to have the signature `mapping(U, mu)`. adjoint_mapping The adjoint function to wrap. If `parameter_type` is `None`, the function is of the form `adjoint_mapping(U)` and is expected to be vectorized. In particular:: adjoint_mapping(U).shape == U.shape[:-1] + (dim_source,). If `parameter_type` is not `None`, the function has to have the signature `adjoint_mapping(U, mu)`. dim_source Dimension of the operator's source. dim_range Dimension of the operator's range. linear Set to `True` if the provided `mapping` and `adjoint_mapping` are linear. parameter_type The |ParameterType| of the |Parameters| the mapping accepts. solver_options The |solver_options| for the operator. name Name of the operator. """ def __init__(self, mapping, adjoint_mapping=None, dim_source=1, dim_range=1, linear=False, parameter_type=None, source_id=None, range_id=None, solver_options=None, name=None): self.__auto_init(locals()) self.source = NumpyVectorSpace(dim_source, source_id) self.range = NumpyVectorSpace(dim_range, range_id)
[docs] def apply(self, U, mu=None): assert U in self.source if self.parametric: mu = self.parse_parameter(mu) return self.range.make_array(self.mapping(U.to_numpy(), mu=mu)) else: return self.range.make_array(self.mapping(U.to_numpy()))
[docs] def apply_adjoint(self, V, mu=None): if self.adjoint_mapping is None: raise ValueError('NumpyGenericOperator: adjoint mapping was not defined.') assert V in self.range V = V.to_numpy() if self.parametric: mu = self.parse_parameter(mu) return self.source.make_array(self.adjoint_mapping(V, mu=mu)) else: return self.source.make_array(self.adjoint_mapping(V))
[docs]class NumpyMatrixBasedOperator(OperatorBase): """Base class for operators which assemble into a |NumpyMatrixOperator|. Attributes ---------- sparse `True` if the operator assembles into a sparse matrix, `False` if the operator assembles into a dense matrix, `None` if unknown. """ linear = True sparse = None @property def H(self): if not self.parametric: return self.assemble().H else: return super().H @abstractmethod def _assemble(self, mu=None): pass
[docs] def assemble(self, mu=None): return NumpyMatrixOperator(self._assemble(self.parse_parameter(mu)), source_id=self.source.id, range_id=self.range.id, solver_options=self.solver_options, name=self.name)
[docs] def apply(self, U, mu=None): return self.assemble(mu).apply(U)
[docs] def apply_adjoint(self, V, mu=None): return self.assemble(mu).apply_adjoint(V)
[docs] def as_range_array(self, mu=None): return self.assemble(mu).as_range_array()
[docs] def as_source_array(self, mu=None): return self.assemble(mu).as_source_array()
[docs] def apply_inverse(self, V, mu=None, least_squares=False): return self.assemble(mu).apply_inverse(V, least_squares=least_squares)
[docs] def export_matrix(self, filename, matrix_name=None, output_format='matlab', mu=None): """Save the matrix of the operator to a file. Parameters ---------- filename Name of output file. matrix_name The name, the output matrix is given. (Comment field is used in case of Matrix Market output_format.) If `None`, the |Operator|'s `name` is used. output_format Output file format. Either `matlab` or `matrixmarket`. mu The |Parameter| to assemble the to be exported matrix for. """ assert output_format in {'matlab', 'matrixmarket'} matrix = self.assemble(mu).matrix matrix_name = matrix_name or self.name if output_format is 'matlab': savemat(filename, {matrix_name: matrix}) else: mmwrite(filename, matrix, comment=matrix_name)
[docs]class NumpyMatrixOperator(NumpyMatrixBasedOperator): """Wraps a 2D |NumPy Array| as an |Operator|. Parameters ---------- matrix The |NumPy array| which is to be wrapped. source_id The id of the operator's `source` |VectorSpace|. range_id The id of the operator's `range` |VectorSpace|. solver_options The |solver_options| for the operator. name Name of the operator. """ def __init__(self, matrix, source_id=None, range_id=None, solver_options=None, name=None): assert matrix.ndim <= 2 if matrix.ndim == 1: matrix = np.reshape(matrix, (1, -1)) try: matrix.setflags(write=False) # make numpy arrays read-only except AttributeError: pass self.__auto_init(locals()) self.source = NumpyVectorSpace(matrix.shape[1], source_id) self.range = NumpyVectorSpace(matrix.shape[0], range_id) self.sparse = issparse(matrix) @classmethod def from_file(cls, path, key=None, source_id=None, range_id=None, solver_options=None, name=None): from pymor.tools.io import load_matrix matrix = load_matrix(path, key=key) return cls(matrix, solver_options=solver_options, source_id=source_id, range_id=range_id, name=name or key or path) @property def H(self): options = {'inverse': self.solver_options.get('inverse_adjoint'), 'inverse_adjoint': self.solver_options.get('inverse')} if self.solver_options else None if self.sparse: adjoint_matrix = self.matrix.transpose(copy=False).conj(copy=False) elif np.isrealobj(self.matrix): adjoint_matrix = self.matrix.T else: adjoint_matrix = self.matrix.T.conj() return self.with_(matrix=adjoint_matrix, source_id=self.range_id, range_id=self.source_id, solver_options=options, name=self.name + '_adjoint') def _assemble(self, mu=None): pass
[docs] def assemble(self, mu=None): return self
[docs] def as_range_array(self, mu=None): return self.range.make_array(self.matrix.T.copy())
[docs] def as_source_array(self, mu=None): return self.source.make_array(self.matrix.copy()).conj()
[docs] def apply(self, U, mu=None): assert U in self.source return self.range.make_array(self.matrix.dot(U.to_numpy().T).T)
[docs] def apply_adjoint(self, V, mu=None): assert V in self.range return self.H.apply(V, mu=mu)
[docs] @defaults('check_finite', 'default_sparse_solver_backend') def apply_inverse(self, V, mu=None, least_squares=False, check_finite=True, default_sparse_solver_backend='scipy'): """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. check_finite Test if solution only contains finite values. default_sparse_solver_backend Default sparse solver backend to use (scipy, pyamg, generic). Returns ------- |VectorArray| of the inverse operator evaluations. Raises ------ InversionError The operator could not be inverted. """ assert V in self.range if V.dim == 0: if self.source.dim == 0 or least_squares: return self.source.make_array(np.zeros((len(V), self.source.dim))) else: raise InversionError if self.source.dim != self.range.dim and not least_squares: raise InversionError options = self.solver_options.get('inverse') if self.solver_options else None assert self.sparse or not options if self.sparse: if options: solver = options if isinstance(options, str) else options['type'] backend = solver.split('_')[0] else: backend = default_sparse_solver_backend if backend == 'scipy': from pymor.bindings.scipy import apply_inverse as apply_inverse_impl elif backend == 'pyamg': if not config.HAVE_PYAMG: raise RuntimeError('PyAMG support not enabled.') from pymor.bindings.pyamg import apply_inverse as apply_inverse_impl elif backend == 'generic': logger = getLogger('pymor.bindings.scipy.scipy_apply_inverse') logger.warning('You have selected a (potentially slow) generic solver for a NumPy matrix operator!') from pymor.algorithms.genericsolvers import apply_inverse as apply_inverse_impl else: raise NotImplementedError return apply_inverse_impl(self, V, options=options, least_squares=least_squares, check_finite=check_finite) else: if least_squares: try: R, _, _, _ = np.linalg.lstsq(self.matrix, V.to_numpy().T) except np.linalg.LinAlgError as e: raise InversionError(f'{str(type(e))}: {str(e)}') R = R.T else: try: R = solve(self.matrix, V.to_numpy().T).T except np.linalg.LinAlgError as e: raise InversionError(f'{str(type(e))}: {str(e)}') if check_finite: if not np.isfinite(np.sum(R)): raise InversionError('Result contains non-finite values') return self.source.make_array(R)
[docs] def apply_inverse_adjoint(self, U, mu=None, least_squares=False): return self.H.apply_inverse(U, mu=mu, least_squares=least_squares)
def _assemble_lincomb(self, operators, coefficients, identity_shift=0., solver_options=None, name=None): if not all(isinstance(op, NumpyMatrixOperator) for op in operators): return None common_mat_dtype = reduce(np.promote_types, (op.matrix.dtype for op in operators if hasattr(op, 'matrix'))) common_coef_dtype = reduce(np.promote_types, (type(c) for c in coefficients + [identity_shift])) common_dtype = np.promote_types(common_mat_dtype, common_coef_dtype) if coefficients[0] == 1: matrix = operators[0].matrix.astype(common_dtype) else: matrix = operators[0].matrix * coefficients[0] if matrix.dtype != common_dtype: matrix = matrix.astype(common_dtype) for op, c in zip(operators[1:], coefficients[1:]): if c == 1: try: matrix += op.matrix except NotImplementedError: matrix = matrix + op.matrix elif c == -1: try: matrix -= op.matrix except NotImplementedError: matrix = matrix - op.matrix else: try: matrix += (op.matrix * c) except NotImplementedError: matrix = matrix + (op.matrix * c) if identity_shift != 0: if identity_shift.imag == 0: identity_shift = identity_shift.real if operators[0].sparse: try: matrix += (scipy.sparse.eye(matrix.shape[0]) * identity_shift) except NotImplementedError: matrix = matrix + (scipy.sparse.eye(matrix.shape[0]) * identity_shift) else: matrix += (np.eye(matrix.shape[0]) * identity_shift) return NumpyMatrixOperator(matrix, source_id=self.source.id, range_id=self.range.id, solver_options=solver_options) def __getstate__(self): if hasattr(self.matrix, 'factorization'): # remove unplicklable SuperLU factorization del self.matrix.factorization return self.__dict__ def _format_repr(self, max_width, verbosity): if self.sparse: matrix_repr = f'<{self.range.dim}x{self.source.dim} sparse, {self.matrix.nnz} nnz>' else: matrix_repr = f'<{self.range.dim}x{self.source.dim} dense>' return super()._format_repr(max_width, verbosity, override={'matrix': matrix_repr})