From Operator to NumPy operator

docs/source/substitutions.py

@@ -113,6 +113,12 @@
 .. |array| replace:: :class:`NumPy array <numpy.ndarray>`
 .. |Array| replace:: :class:`NumPy array <numpy.ndarray>`
 
+.. |SciPy| replace:: :mod:`SciPy <scipy>`
+.. |SciPy spmatrix| replace:: :class:`SciPy spmatrix <scipy.sparse.spmatrix>`
+.. |SciPy spmatrices| replace:: :class:`SciPy spmatrices <scipy.sparse.spmatrix>`
+.. |Scipy spmatrix| replace:: :class:`SciPy spmatrix <scipy.sparse.spmatrix>`
+.. |Scipy spmatrices| replace:: :class:`SciPy spmatrices <scipy.sparse.spmatrix>`
+
 .. |OrderedDict| replace:: :class:`~collections.OrderedDict`
 
 .. |solver_options| replace:: :attr:`~pymor.operators.interfaces.OperatorInterface.solver_options`

src/pymor/algorithms/to_matrix.py

@@ -0,0 +1,118 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright 2013-2016 pyMOR developers and contributors. All rights reserved.
+# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+import scipy.linalg as spla
+import scipy.sparse as sps
+import scipy.sparse.linalg as spsla
+
+from pymor.operators.block import BlockOperator
+from pymor.operators.constructions import (AdjointOperator, ComponentProjection, Concatenation, IdentityOperator,
+                                           LincombOperator, VectorArrayOperator, ZeroOperator)
+from pymor.operators.numpy import NumpyMatrixOperator
+
+
+def to_matrix(op, format=None, mu=None):
+    """Transfrom construction of NumpyMatrixOperators to NumPy or SciPy array
+
+    Parameters
+    ----------
+    op
+        Operator.
+    format
+        Format of the resulting |SciPy spmatrix|.
+        If `None`, a dense format is used.
+    mu
+        |Parameter|.
+
+    Returns
+    -------
+    res
+        Equivalent matrix.
+    """
+    assert format is None or format in ('bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil')
+    op = op.assemble(mu)
+    mapping = {
+        'bsr': sps.bsr_matrix,
+        'coo': sps.coo_matrix,
+        'csc': sps.csc_matrix,
+        'csr': sps.csr_matrix,
+        'dia': sps.dia_matrix,
+        'dok': sps.dok_matrix,
+        'lil': sps.lil_matrix
+    }
+    return _to_matrix(op, format, mapping, mu)
+
+
+def _to_matrix(op, format, mapping, mu):
+    if isinstance(op, NumpyMatrixOperator):
+        if format is None:
+            if not op.sparse:
+                res = op._matrix
+            else:
+                res = op._matrix.toarray()
+        else:
+            res = mapping[format](op._matrix)
+    elif isinstance(op, BlockOperator):
+        op_blocks = op._blocks
+        mat_blocks = [[] for i in range(op.num_range_blocks)]
+        for i in range(op.num_range_blocks):
+            for j in range(op.num_source_blocks):
+                if op_blocks[i, j] is None:
+                    if format is None:
+                        mat_blocks[i].append(np.zeros((op.range.subtype[i].dim, op.source.subtype[j].dim)))
+                    else:
+                        mat_blocks[i].append(None)
+                else:
+                    mat_blocks[i].append(_to_matrix(op_blocks[i, j], format, mapping, mu))
+        if format is None:
+            res = np.bmat(mat_blocks)
+        else:
+            res = sps.bmat(mat_blocks, format=format)
+    elif isinstance(op, AdjointOperator):
+        res = _to_matrix(op.operator, format, mapping, mu).T
+        if op.range_product is not None:
+            res = res.dot(_to_matrix(op.range_product, format, mapping, mu))
+        if op.source_product is not None:
+            if format is None:
+                res = spla.solve(_to_matrix(op.source_product, format, mapping, mu), res)
+            else:
+                res = spsla.spsolve(_to_matrix(op.source_product, format, mapping, mu), res)
+    elif isinstance(op, ComponentProjection):
+        if format is None:
+            res = np.zeros((op.range.dim, op.source.dim))
+            for i, j in enumerate(op.components):
+                res[i, j] = 1
+        else:
+            data = np.ones((op.range.dim,))
+            i = np.arange(op.range.dim)
+            j = op.components
+            res = sps.coo_matrix((data, (i, j)), shape=(op.range.dim, op.source.dim))
+            res = res.asformat(format)
+    elif isinstance(op, Concatenation):
+        res = _to_matrix(op.second, format, mapping, mu).dot(_to_matrix(op.first, format, mapping, mu))
+    elif isinstance(op, IdentityOperator):
+        if format is None:
+            res = np.eye(op.source.dim)
+        else:
+            res = sps.eye(op.source.dim, format=format)
+    elif isinstance(op, LincombOperator):
+        op_coefficients = op.evaluate_coefficients(mu)
+        res = op_coefficients[0] * _to_matrix(op.operators[0], format, mapping, mu)
+        for i in range(1, len(op.operators)):
+            res = res + op_coefficients[i] * _to_matrix(op.operators[i], format, mapping, mu)
+    elif isinstance(op, VectorArrayOperator):
+        res = op._array.data if op.transposed else op._array.data.T
+        if format is not None:
+            res = mapping[format](res)
+    elif isinstance(op, ZeroOperator):
+        if format is None:
+            res = np.zeros((op.range.dim, op.source.dim))
+        else:
+            res = mapping[format]((op.range.dim, op.source.dim))
+    else:
+        raise ValueError('Encountered unsupported operator type {}'.format(type(op)))
+    return res

src/pymortests/to_matrix.py

@@ -0,0 +1,44 @@
+# This file is part of the pyMOR project (http://www.pymor.org).
+# Copyright 2013-2016 pyMOR developers and contributors. All rights reserved.
+# License: BSD 2-Clause License (http://opensource.org/licenses/BSD-2-Clause)
+
+from __future__ import absolute_import, division, print_function
+
+import numpy as np
+import scipy.sparse as sps
+
+from pymor.algorithms.to_matrix import to_matrix
+from pymor.operators.block import BlockDiagonalOperator
+from pymor.operators.constructions import (AdjointOperator, Concatenation, IdentityOperator, LincombOperator,
+                                           VectorArrayOperator)
+from pymor.operators.numpy import NumpyMatrixOperator
+from pymor.vectorarrays.numpy import NumpyVectorArray, NumpyVectorSpace
+
+
+def test_to_matrix():
+    np.random.seed(0)
+    A = np.random.randn(2, 2)
+    B = np.random.randn(3, 3)
+    C = np.random.randn(3, 3)
+
+    X = np.bmat([[np.eye(2) + A, np.zeros((2, 3))], [np.zeros((3, 2)), B.dot(C.T)]])
+
+    C = sps.csc_matrix(C)
+
+    Aop = NumpyMatrixOperator(A)
+    Bop = NumpyMatrixOperator(B)
+    Cop = NumpyMatrixOperator(C)
+
+    Xop = BlockDiagonalOperator([LincombOperator([IdentityOperator(NumpyVectorSpace(2)), Aop],
+                                                 [1, 1]), Concatenation(Bop, AdjointOperator(Cop))])
+
+    assert np.allclose(X, to_matrix(Xop))
+    assert np.allclose(X, to_matrix(Xop, format='csr').toarray())
+
+    np.random.seed(0)
+    V = np.random.randn(10, 2)
+    Vva = NumpyVectorArray(V.T)
+    Vop = VectorArrayOperator(Vva)
+    assert np.allclose(V, to_matrix(Vop))
+    Vop = VectorArrayOperator(Vva, transposed=True)
+    assert np.allclose(V, to_matrix(Vop).T)