Merge dd44796 into 8500dc2

src/pymor/algorithms/numpy.py

@@ -0,0 +1,117 @@
+# 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_numpy_operator(op, format=None):
+    """Transfrom construction of NumpyMatrixOperators to NumpyMatrixOperator
+
+    Parameters
+    ----------
+    op
+        Operator.
+    format
+        Format of the |SciPy sparray| in the resulting NumpyMatrixOperator.
+        If `None`, a dense format is used.
+
+    Returns
+    -------
+    res
+        Equivalent NumpyMatrixOperator.
+    """
+    assert format is None or format in ('bsr', 'coo', 'csc', 'csr', 'dia', 'dok', 'lil')
+    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 NumpyMatrixOperator(_to_matrix(op, format, mapping))
+
+
+def _to_matrix(op, format, mapping):
+    if isinstance(op, NumpyMatrixOperator):
+        if format is None:
+            if not op.sparse:
+                res = op._matrix
+            else:
+                res = op._matrix.toarray()
+        else:
+            if not op.sparse:
+                res = mapping[format](op._matrix)
+            else:
+                res = op._matrix.asformat(format)
+    elif isinstance(op, BlockOperator):
+        op_blocks = op._blocks
+        mat_blocks = [[] for i in xrange(op.num_range_blocks)]
+        for i in xrange(op.num_range_blocks):
+            for j in xrange(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))
+        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).T
+        if op.range_product is not None:
+            res = res.dot(_to_matrix(op.range_product, format, mapping))
+        if op.source_product is not None:
+            if format is None:
+                res = spla.solve(_to_matrix(op.source_product, format, mapping), res)
+            else:
+                res = spsla.spsolve(_to_matrix(op.source_product, format, mapping), 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).dot(_to_matrix(op.first, format, mapping))
+    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):
+        res = op.coefficients[0] * _to_matrix(op.operators[0], format, mapping)
+        for i in xrange(1, len(op.operators)):
+            res = res + op.coefficients[i] * _to_matrix(op.operators[i], format, mapping)
+    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_numpy_operator.py

@@ -0,0 +1,43 @@
+# 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.numpy import to_numpy_operator
+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_numpy_operator():
+    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_numpy_operator(Xop)._matrix)
+
+    np.random.seed(0)
+    V = np.random.randn(10, 2)
+    Vva = NumpyVectorArray(V.T)
+    Vop = VectorArrayOperator(Vva)
+    assert np.allclose(V, to_numpy_operator(Vop)._matrix)
+    Vop = VectorArrayOperator(Vva, transposed=True)
+    assert np.allclose(V, to_numpy_operator(Vop)._matrix.T)