From 4a07f1c9c36ecd47e3fce1738c87d42615d53f8c Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:12:59 +0100 Subject: [PATCH 01/26] DEV: Add a function for matrix representation. See issue #49. --- odl/operator/utility.py | 77 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 77 insertions(+) create mode 100644 odl/operator/utility.py diff --git a/odl/operator/utility.py b/odl/operator/utility.py new file mode 100644 index 00000000000..64403caf453 --- /dev/null +++ b/odl/operator/utility.py @@ -0,0 +1,77 @@ +# Copyright 2014, 2015 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() + +# External +import numpy as np + +# Internal +from odl.space.base_ntuples import FnBase +from odl.set.pspace import ProductSpace + + +def matrix_representation(op): + """Returns a matrix representation of a linear operator. + + Parameters + ---------- + op : :class:`~odl.Operator` + The linear operator of which one wants a matrix representation. + + Returns + ---------- + matrix : `numpy.ndarray` + The matrix representation of the operator. + + Notes + ---------- + The algorithm works by letting the operator act on all unit vectors, and + stacking the output as a matrix. + + """ + + if not op.is_linear: + print('WARNING: The operator is not linear; cannot produce matrix', + 'representation of it.') + return + + if not isinstance(op.domain, FnBase) or isinstance(op.domain, ProductSpace): + print('WARNING: The operator domain is not discreate or produc space;', + 'cannot produce matrix representation of it.') + return + + if not isinstance(op.range, FnBase): + print('WARNING: The operator range is not discreate; cannot produce', + 'matrix representation of it.') + return + + n = op.range.size + m = op.domain.size + + matrix = np.zeros([n, m]) + v = op.domain.element() + tmp = op.range.element() + for i in range(m): + v.set_zero() + v[i] = 1.0 + matrix[:,i] = op(v, out=tmp) + + return matrix \ No newline at end of file From 8fddde260af6ea1e6de8491dd99dca671634327c Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:28:09 +0100 Subject: [PATCH 02/26] TST: Add test for the matrix representation function. Add a small test in R3 for the matrix representation function, to see that one gets the same matrix back. --- test/operator/utility_test.py | 70 +++++++++++++++++++++++++++++++++++ 1 file changed, 70 insertions(+) create mode 100644 test/operator/utility_test.py diff --git a/test/operator/utility_test.py b/test/operator/utility_test.py new file mode 100644 index 00000000000..0757ee22e0b --- /dev/null +++ b/test/operator/utility_test.py @@ -0,0 +1,70 @@ +# Copyright 2014, 2015 The ODL development group +# +# This file is part of ODL. +# +# ODL is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ODL is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with ODL. If not, see . + + +# Imports for common Python 2/3 codebase +from __future__ import print_function, division, absolute_import +from future import standard_library +standard_library.install_aliases() +from builtins import str, super + +# External module imports +import pytest +import numpy as np + +# ODL imports +import odl +from odl.operator.utility import matrix_representation +from odl.util.testutils import almost_equal + + +class MultiplyOp(odl.Operator): + + """Multiply with matrix. + """ + + def __init__(self, matrix, domain=None, range=None): + domain = (odl.Rn(matrix.shape[1]) + if domain is None else domain) + range = (odl.Rn(matrix.shape[0]) + if range is None else range) + self.matrix = matrix + + super().__init__(domain, range, linear=True) + + def _apply(self, rhs, out): + np.dot(self.matrix, rhs.data, out=out.data) + + @property + def adjoint(self): + return MultiplyOp(self.matrix.T, self.range, self.domain) + + +def test_matrix_representation(): + # Verify that the matrix representation function returns the correct matrix + + A = np.random.rand(3, 3) + + Aop = MultiplyOp(A) + + the_matrix = matrix_representation(Aop) + + assert almost_equal(np.sum(np.abs(A - the_matrix)), 1e-6) + + +if __name__ == '__main__': + pytest.main(str(__file__.replace('\\', '/')) + ' -v') From 76646d405a94a8b4bc51a74195206a3f2558661f Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:29:50 +0100 Subject: [PATCH 03/26] STY: pep8 changes. --- odl/operator/utility.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/operator/utility.py b/odl/operator/utility.py index 64403caf453..817d0a2909d 100644 --- a/odl/operator/utility.py +++ b/odl/operator/utility.py @@ -72,6 +72,6 @@ def matrix_representation(op): for i in range(m): v.set_zero() v[i] = 1.0 - matrix[:,i] = op(v, out=tmp) + matrix[:, i] = op(v, out=tmp) - return matrix \ No newline at end of file + return matrix From 39a44003ace3ea6acda5a75ccf9988dd9ac38182 Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:35:31 +0100 Subject: [PATCH 04/26] MAINT: Move to discr. See #49. --- odl/{operator => discr}/utility.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename odl/{operator => discr}/utility.py (100%) diff --git a/odl/operator/utility.py b/odl/discr/utility.py similarity index 100% rename from odl/operator/utility.py rename to odl/discr/utility.py From beeb9cd774cf8eaeff12f25b60d3f23cd2c2dd2e Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:37:40 +0100 Subject: [PATCH 05/26] MAINT: Move to discr. See #49. --- test/{operator => discr}/utility_test.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename test/{operator => discr}/utility_test.py (100%) diff --git a/test/operator/utility_test.py b/test/discr/utility_test.py similarity index 100% rename from test/operator/utility_test.py rename to test/discr/utility_test.py From 0c53b52ca04fd76793e5f80f4bc6e609cafc5eee Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:38:35 +0100 Subject: [PATCH 06/26] BUG: Fix bug due to the move. See #49. --- test/discr/utility_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/discr/utility_test.py b/test/discr/utility_test.py index 0757ee22e0b..dc76418fafb 100644 --- a/test/discr/utility_test.py +++ b/test/discr/utility_test.py @@ -28,7 +28,7 @@ # ODL imports import odl -from odl.operator.utility import matrix_representation +from odl.discr.utility import matrix_representation from odl.util.testutils import almost_equal From 9e7f91dba200422ad56ff1dadd52445fe70ad7fc Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 10:56:51 +0100 Subject: [PATCH 07/26] BUG/STY: Change to pep8, and solved a bug on the same line. See #49. --- odl/discr/utility.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/odl/discr/utility.py b/odl/discr/utility.py index 817d0a2909d..431da6e7c5a 100644 --- a/odl/discr/utility.py +++ b/odl/discr/utility.py @@ -53,7 +53,8 @@ def matrix_representation(op): 'representation of it.') return - if not isinstance(op.domain, FnBase) or isinstance(op.domain, ProductSpace): + if not (isinstance(op.domain, FnBase) or + isinstance(op.domain, ProductSpace)): print('WARNING: The operator domain is not discreate or produc space;', 'cannot produce matrix representation of it.') return From b4dc31a3e8ef2289993cba4f4c69fa730465432b Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 15:23:05 +0100 Subject: [PATCH 08/26] API/STY: Move matrix representation to new file to remove errors when building. See 49. --- odl/discr/{utility.py => operator_utilities.py} | 0 test/discr/{utility_test.py => operator_utilities_test.py} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename odl/discr/{utility.py => operator_utilities.py} (100%) rename test/discr/{utility_test.py => operator_utilities_test.py} (100%) diff --git a/odl/discr/utility.py b/odl/discr/operator_utilities.py similarity index 100% rename from odl/discr/utility.py rename to odl/discr/operator_utilities.py diff --git a/test/discr/utility_test.py b/test/discr/operator_utilities_test.py similarity index 100% rename from test/discr/utility_test.py rename to test/discr/operator_utilities_test.py From 56b520798d7c7fd31b1d1a66a9c5ec985ec71920 Mon Sep 17 00:00:00 2001 From: aringh Date: Mon, 23 Nov 2015 15:24:10 +0100 Subject: [PATCH 09/26] TST: Change test to call correct function. In agreement with previous commit (which renamed things). See #49. --- test/discr/operator_utilities_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/discr/operator_utilities_test.py b/test/discr/operator_utilities_test.py index dc76418fafb..b0ed47be969 100644 --- a/test/discr/operator_utilities_test.py +++ b/test/discr/operator_utilities_test.py @@ -28,7 +28,7 @@ # ODL imports import odl -from odl.discr.utility import matrix_representation +from odl.discr.operator_utilities import matrix_representation from odl.util.testutils import almost_equal From d7b4ff22d6d1df4d349603927609a7f2a19e13c9 Mon Sep 17 00:00:00 2001 From: aringh Date: Tue, 24 Nov 2015 10:26:53 +0100 Subject: [PATCH 10/26] ENH: Update according to comments. See #49 and #73. --- odl/discr/operator_utilities.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/odl/discr/operator_utilities.py b/odl/discr/operator_utilities.py index 431da6e7c5a..756e8dda88c 100644 --- a/odl/discr/operator_utilities.py +++ b/odl/discr/operator_utilities.py @@ -15,6 +15,10 @@ # You should have received a copy of the GNU General Public License # along with ODL. If not, see . +"""Usefull utility functions on discreate spaces (i.e., either Rn/Cn or +discretized function spaces), for example obtaining a matrix representation of +an operator. """ + # Imports for common Python 2/3 codebase from __future__ import print_function, division, absolute_import from future import standard_library @@ -53,8 +57,7 @@ def matrix_representation(op): 'representation of it.') return - if not (isinstance(op.domain, FnBase) or - isinstance(op.domain, ProductSpace)): + if not isinstance(op.domain, FnBase): print('WARNING: The operator domain is not discreate or produc space;', 'cannot produce matrix representation of it.') return From ab5e8d0a51832424fb097510efea3cafa307aebd Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 09:38:00 +0100 Subject: [PATCH 11/26] MAINT: Move files to new location. --- odl/{discr => operator}/operator_utilities.py | 0 test/{discr => operator}/operator_utilities_test.py | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename odl/{discr => operator}/operator_utilities.py (100%) rename test/{discr => operator}/operator_utilities_test.py (100%) diff --git a/odl/discr/operator_utilities.py b/odl/operator/operator_utilities.py similarity index 100% rename from odl/discr/operator_utilities.py rename to odl/operator/operator_utilities.py diff --git a/test/discr/operator_utilities_test.py b/test/operator/operator_utilities_test.py similarity index 100% rename from test/discr/operator_utilities_test.py rename to test/operator/operator_utilities_test.py From 5a443f32fb5a65cb37a3a81d409a61ec3b36553f Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 09:42:56 +0100 Subject: [PATCH 12/26] TST: Update test according to comments and commit d167e34. --- test/operator/operator_utilities_test.py | 32 +++++------------------- 1 file changed, 6 insertions(+), 26 deletions(-) diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index b0ed47be969..1d5ba4e92e3 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -20,7 +20,6 @@ from __future__ import print_function, division, absolute_import from future import standard_library standard_library.install_aliases() -from builtins import str, super # External module imports import pytest @@ -28,38 +27,19 @@ # ODL imports import odl -from odl.discr.operator_utilities import matrix_representation +from odl.operator.operator_utilities import matrix_representation +from odl.space.ntuples import MatVecOperator from odl.util.testutils import almost_equal -class MultiplyOp(odl.Operator): - - """Multiply with matrix. - """ - - def __init__(self, matrix, domain=None, range=None): - domain = (odl.Rn(matrix.shape[1]) - if domain is None else domain) - range = (odl.Rn(matrix.shape[0]) - if range is None else range) - self.matrix = matrix - - super().__init__(domain, range, linear=True) - - def _apply(self, rhs, out): - np.dot(self.matrix, rhs.data, out=out.data) - - @property - def adjoint(self): - return MultiplyOp(self.matrix.T, self.range, self.domain) - - def test_matrix_representation(): # Verify that the matrix representation function returns the correct matrix - A = np.random.rand(3, 3) + n = 3 + rn = odl.Rn(n) + A = np.random.rand(n, n) - Aop = MultiplyOp(A) + Aop = MatVecOperator(rn, rn, A) #MultiplyOp(A) the_matrix = matrix_representation(Aop) From 76b0cead3a7b08e6aeba6d51960461ea686620f3 Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 09:51:10 +0100 Subject: [PATCH 13/26] TST: Updated the operator tests to use MatVecOperator. See the comment in #73 (which actually concerns something else but where it was found). --- test/operator/operator_test.py | 80 ++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 32 deletions(-) diff --git a/test/operator/operator_test.py b/test/operator/operator_test.py index 8738a5a286e..119cc87b4e1 100644 --- a/test/operator/operator_test.py +++ b/test/operator/operator_test.py @@ -30,8 +30,8 @@ import odl from odl import (Operator, OperatorSum, OperatorComp, OperatorLeftScalarMult, OperatorRightScalarMult, - FunctionalLeftVectorMult, - OperatorRightVectorMult) + FunctionalLeftVectorMult, OperatorRightVectorMult, + MatVecOperator) from odl.util.testutils import almost_equal, all_almost_equal @@ -184,26 +184,26 @@ def test_nonlinear_composition(): C = OperatorComp(Bop, Aop) -class MultiplyOp(Operator): - - """Multiply with matrix. - """ - - def __init__(self, matrix, domain=None, range=None): - domain = (odl.Rn(matrix.shape[1]) - if domain is None else domain) - range = (odl.Rn(matrix.shape[0]) - if range is None else range) - self.matrix = matrix - - super().__init__(domain, range, linear=True) - - def _apply(self, rhs, out): - np.dot(self.matrix, rhs.data, out=out.data) - - @property - def adjoint(self): - return MultiplyOp(self.matrix.T, self.range, self.domain) +#class MultiplyOp(Operator): +# +# """Multiply with matrix. +# """ +# +# def __init__(self, matrix, domain=None, range=None): +# domain = (odl.Rn(matrix.shape[1]) +# if domain is None else domain) +# range = (odl.Rn(matrix.shape[0]) +# if range is None else range) +# self.matrix = matrix +# +# super().__init__(domain, range, linear=True) +# +# def _apply(self, rhs, out): +# np.dot(self.matrix, rhs.data, out=out.data) +# +# @property +# def adjoint(self): +# return MultiplyOp(self.matrix.T, self.range, self.domain) def test_linear_Op(): @@ -213,7 +213,9 @@ def test_linear_Op(): x = np.random.rand(3) out = np.random.rand(3) - Aop = MultiplyOp(A) + r3 = odl.Rn(3) + + Aop = MatVecOperator(r3, r3, A) xvec = Aop.domain.element(x) outvec = Aop.range.element() @@ -232,7 +234,10 @@ def test_linear_op_nonsquare(): x = np.random.rand(3) out = np.random.rand(4) - Aop = MultiplyOp(A) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + Aop = MatVecOperator(r3, r4, A) + xvec = Aop.domain.element(x) outvec = Aop.range.element() @@ -250,7 +255,9 @@ def test_linear_adjoint(): x = np.random.rand(4) out = np.random.rand(3) - Aop = MultiplyOp(A) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + Aop = MatVecOperator(r3, r4, A) xvec = Aop.range.element(x) outvec = Aop.domain.element() @@ -269,8 +276,10 @@ def test_linear_addition(): x = np.random.rand(3) y = np.random.rand(4) - Aop = MultiplyOp(A) - Bop = MultiplyOp(B) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + Aop = MatVecOperator(r3, r4, A) + Bop = MatVecOperator(r3, r4, B) xvec = Aop.domain.element(x) yvec = Aop.range.element(y) @@ -295,7 +304,9 @@ def test_linear_scale(): x = np.random.rand(3) y = np.random.rand(4) - Aop = MultiplyOp(A) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + Aop = MatVecOperator(r3, r4, A) xvec = Aop.domain.element(x) yvec = Aop.range.element(y) @@ -325,7 +336,9 @@ def test_linear_scale(): def test_linear_right_vector_mult(): A = np.random.rand(4, 3) - Aop = MultiplyOp(A) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + Aop = MatVecOperator(r3, r4, A) vec = Aop.domain.element([1, 2, 3]) x = Aop.domain.element([4, 5, 6]) y = Aop.range.element([5, 6, 7, 8]) @@ -354,8 +367,11 @@ def test_linear_composition(): x = np.random.rand(3) y = np.random.rand(5) - Aop = MultiplyOp(A) - Bop = MultiplyOp(B) + r3 = odl.Rn(3) + r4 = odl.Rn(4) + r5 = odl.Rn(5) + Aop = MatVecOperator(r4, r5, A) + Bop = MatVecOperator(r3, r4, B) xvec = Bop.domain.element(x) yvec = Aop.range.element(y) @@ -372,7 +388,7 @@ def test_type_errors(): r3 = odl.Rn(3) r4 = odl.Rn(4) - Aop = MultiplyOp(np.random.rand(3, 3)) + Aop = MatVecOperator(r3, r3, np.random.rand(3, 3)) r3Vec1 = r3.zero() r3Vec2 = r3.zero() r4Vec1 = r4.zero() From 96115b8a4ec85d8d003efb86c1537466955e5508 Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 09:54:03 +0100 Subject: [PATCH 14/26] TST: Remove the oldmatrix multiplication operator. See #73, and compare commit 792b359 where it was just commented. --- test/operator/operator_test.py | 22 ---------------------- 1 file changed, 22 deletions(-) diff --git a/test/operator/operator_test.py b/test/operator/operator_test.py index 119cc87b4e1..3a7e5a274df 100644 --- a/test/operator/operator_test.py +++ b/test/operator/operator_test.py @@ -184,28 +184,6 @@ def test_nonlinear_composition(): C = OperatorComp(Bop, Aop) -#class MultiplyOp(Operator): -# -# """Multiply with matrix. -# """ -# -# def __init__(self, matrix, domain=None, range=None): -# domain = (odl.Rn(matrix.shape[1]) -# if domain is None else domain) -# range = (odl.Rn(matrix.shape[0]) -# if range is None else range) -# self.matrix = matrix -# -# super().__init__(domain, range, linear=True) -# -# def _apply(self, rhs, out): -# np.dot(self.matrix, rhs.data, out=out.data) -# -# @property -# def adjoint(self): -# return MultiplyOp(self.matrix.T, self.range, self.domain) - - def test_linear_Op(): # Verify that the multiply op does indeed work as expected From 410246f9cf73df5c23b066d4191ee9a9ae669a9e Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 09:58:21 +0100 Subject: [PATCH 15/26] MAINT: Minor last changes. --- odl/operator/operator_utilities.py | 7 +++---- test/operator/operator_utilities_test.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/odl/operator/operator_utilities.py b/odl/operator/operator_utilities.py index 756e8dda88c..4c301c0f8ee 100644 --- a/odl/operator/operator_utilities.py +++ b/odl/operator/operator_utilities.py @@ -15,7 +15,7 @@ # You should have received a copy of the GNU General Public License # along with ODL. If not, see . -"""Usefull utility functions on discreate spaces (i.e., either Rn/Cn or +"""Usefull utility functions on discrete spaces (i.e., either Rn/Cn or discretized function spaces), for example obtaining a matrix representation of an operator. """ @@ -29,7 +29,6 @@ # Internal from odl.space.base_ntuples import FnBase -from odl.set.pspace import ProductSpace def matrix_representation(op): @@ -58,12 +57,12 @@ def matrix_representation(op): return if not isinstance(op.domain, FnBase): - print('WARNING: The operator domain is not discreate or produc space;', + print('WARNING: The operator domain is not discrete or produc space;', 'cannot produce matrix representation of it.') return if not isinstance(op.range, FnBase): - print('WARNING: The operator range is not discreate; cannot produce', + print('WARNING: The operator range is not discrete; cannot produce', 'matrix representation of it.') return diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index 1d5ba4e92e3..b710bf44f00 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -39,7 +39,7 @@ def test_matrix_representation(): rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) #MultiplyOp(A) + Aop = MatVecOperator(rn, rn, A) the_matrix = matrix_representation(Aop) From d41acdec23c70b8028a1e34a55adc5540196f893 Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 17:24:16 +0100 Subject: [PATCH 16/26] ENH: Update so that matrix representation can handle product spaces. This version of the matrix representation function can handle operators on product spaces, both in range and domain. See #49 --- odl/operator/operator_utilities.py | 60 +++++++++++++++++++++++------- 1 file changed, 47 insertions(+), 13 deletions(-) diff --git a/odl/operator/operator_utilities.py b/odl/operator/operator_utilities.py index 4c301c0f8ee..942954322bc 100644 --- a/odl/operator/operator_utilities.py +++ b/odl/operator/operator_utilities.py @@ -29,6 +29,7 @@ # Internal from odl.space.base_ntuples import FnBase +from odl.set.pspace import ProductSpace def matrix_representation(op): @@ -56,25 +57,58 @@ def matrix_representation(op): 'representation of it.') return - if not isinstance(op.domain, FnBase): + if not (isinstance(op.domain, FnBase) or + isinstance(op.domain, ProductSpace)): print('WARNING: The operator domain is not discrete or produc space;', 'cannot produce matrix representation of it.') return - if not isinstance(op.range, FnBase): + if not (isinstance(op.range, FnBase) or + isinstance(op.range, ProductSpace)): print('WARNING: The operator range is not discrete; cannot produce', 'matrix representation of it.') - return - - n = op.range.size - m = op.domain.size - matrix = np.zeros([n, m]) - v = op.domain.element() - tmp = op.range.element() - for i in range(m): - v.set_zero() - v[i] = 1.0 - matrix[:, i] = op(v, out=tmp) + # Get the size of the range, and handle ProductSpace + op_ran = op.range + op_ran_is_prod_space = isinstance(op_ran, ProductSpace) + if op_ran_is_prod_space: + num_ran = op_ran.size + n = [op_ran[i].size for i in range(num_ran)] + else: + num_ran = 1 + n = [op_ran.size] + + # Get the size of the domain, and handle ProductSpace + op_dom = op.domain + op_dom_is_prod_space = isinstance(op_dom, ProductSpace) + if op_dom_is_prod_space: + num_dom = op_dom.size + m = [op_dom[i].size for i in range(num_dom)] + else: + num_dom = 1 + m = [op_dom.size] + + # Generate the matrix + matrix = np.zeros([np.sum(n), np.sum(m)]) + tmp_1 = op_ran.element() + v = op_dom.element() + index = 0 + + for i in range(num_dom): + for j in range(m[i]): + v.set_zero() + if op_dom_is_prod_space: + v[i][j] = 1.0 + else: + v[j] = 1.0 + tmp_2 = op(v, out=tmp_1) + tmp_3 = [] + if op_ran_is_prod_space: + for k in range(num_ran): + tmp_3 = np.concatenate((tmp_3, tmp_2[k].asarray())) + else: + tmp_3 = tmp_2.asarray() + matrix[:, index] = tmp_3 + index += 1 return matrix From 281ce2cb75a2735f0ca7d1e21a7fbfc811e74b00 Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 17:24:51 +0100 Subject: [PATCH 17/26] TST: Update of test for matrix representation. See #49. --- test/operator/operator_utilities_test.py | 96 ++++++++++++++++++++++++ 1 file changed, 96 insertions(+) diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index b710bf44f00..e39450ae65f 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -28,6 +28,9 @@ # ODL imports import odl from odl.operator.operator_utilities import matrix_representation +from odl.set.pspace import ProductSpace +from odl.operator.pspace_ops import ProductSpaceOperator + from odl.space.ntuples import MatVecOperator from odl.util.testutils import almost_equal @@ -46,5 +49,98 @@ def test_matrix_representation(): assert almost_equal(np.sum(np.abs(A - the_matrix)), 1e-6) +def test_three_matrix_representation_product_to_lin_space(): + # Verify that the matrix representation function returns the correct matrix + + n = 3 + rn = odl.Rn(n) + A = np.random.rand(n, n) + Aop = MatVecOperator(rn, rn, A) + + m = 2 + rm = odl.Rn(m) + B = np.random.rand(n, m) + Bop = MatVecOperator(rm, rn, B) + + dom = ProductSpace(rn, rm) + ran = ProductSpace(rn, 1) + + AB_matrix = np.hstack([A, B]) + ABop = ProductSpaceOperator([Aop, Bop], dom, ran) + + the_matrix = matrix_representation(ABop) + + assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) + + +def test_matrix_representation_lin_space_to_product(): + # Verify that the matrix representation function returns the correct matrix + + n = 3 + rn = odl.Rn(n) + A = np.random.rand(n, n) + Aop = MatVecOperator(rn, rn, A) + + m = 2 + rm = odl.Rn(m) + B = np.random.rand(m, n) + Bop = MatVecOperator(rn, rm, B) + + dom = ProductSpace(rn, 1) + ran = ProductSpace(rn, rm) + + AB_matrix = np.vstack([A, B]) + ABop = ProductSpaceOperator([[Aop], [Bop]], dom, ran) + + the_matrix = matrix_representation(ABop) + + assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) + + +def test_matrix_representation_product_to_product(): + # Verify that the matrix representation function returns the correct matrix + + n = 3 + rn = odl.Rn(n) + A = np.random.rand(n, n) + Aop = MatVecOperator(rn, rn, A) + + m = 2 + rm = odl.Rn(m) + B = np.random.rand(m, m) + Bop = MatVecOperator(rm, rm, B) + + ran_and_dom = ProductSpace(rn, rm) + + AB_matrix = np.vstack([np.hstack([A, np.zeros((n, m))]), + np.hstack([np.zeros((m, n)), B])]) + ABop = ProductSpaceOperator([[Aop, None], [None, Bop]], + ran_and_dom, ran_and_dom) + the_matrix = matrix_representation(ABop) + + assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) + + +def test_five_matrix_representation_product_to_product_two(): + # Verify that the matrix representation function returns the correct matrix + + n = 3 + rn = odl.Rn(n) + A = np.random.rand(n, n) + Aop = MatVecOperator(rn, rn, A) + + B = np.random.rand(n, n) + Bop = MatVecOperator(rn, rn, B) + + ran_and_dom = ProductSpace(rn, 2) + + AB_matrix = np.vstack([np.hstack([A, np.zeros((n, n))]), + np.hstack([np.zeros((n, n)), B])]) + ABop = ProductSpaceOperator([[Aop, None], [None, Bop]], + ran_and_dom, ran_and_dom) + the_matrix = matrix_representation(ABop) + + assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) + if __name__ == '__main__': pytest.main(str(__file__.replace('\\', '/')) + ' -v') From 751cc9ceae52c5bac029bb95f71344ad85816412 Mon Sep 17 00:00:00 2001 From: aringh Date: Wed, 25 Nov 2015 17:35:07 +0100 Subject: [PATCH 18/26] MAINT: Correct typos. --- test/operator/operator_utilities_test.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index e39450ae65f..5225bc1d6f1 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -49,7 +49,7 @@ def test_matrix_representation(): assert almost_equal(np.sum(np.abs(A - the_matrix)), 1e-6) -def test_three_matrix_representation_product_to_lin_space(): +def test_matrix_representation_product_to_lin_space(): # Verify that the matrix representation function returns the correct matrix n = 3 @@ -121,7 +121,7 @@ def test_matrix_representation_product_to_product(): assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) -def test_five_matrix_representation_product_to_product_two(): +def test_matrix_representation_product_to_product_two(): # Verify that the matrix representation function returns the correct matrix n = 3 From dc77b712c9d3f5e7e0b2bc84d43e5c1bca9efc09 Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 09:44:09 +0100 Subject: [PATCH 19/26] ENH: Update according to comments. Comments found in #73. See also #49 for the issue. --- odl/operator/operator_utilities.py | 62 +++++++++++++++--------------- 1 file changed, 32 insertions(+), 30 deletions(-) diff --git a/odl/operator/operator_utilities.py b/odl/operator/operator_utilities.py index 942954322bc..e8964716054 100644 --- a/odl/operator/operator_utilities.py +++ b/odl/operator/operator_utilities.py @@ -53,62 +53,64 @@ def matrix_representation(op): """ if not op.is_linear: - print('WARNING: The operator is not linear; cannot produce matrix', - 'representation of it.') - return + raise ValueError('The operator is not linear') if not (isinstance(op.domain, FnBase) or - isinstance(op.domain, ProductSpace)): - print('WARNING: The operator domain is not discrete or produc space;', - 'cannot produce matrix representation of it.') - return + (isinstance(op.domain, ProductSpace) and + all(isinstance(spc, FnBase) for spc in op.domain))): + raise TypeError('Operator domain {} is not FnBase, nor ProductSpace' + 'with only FnBase components'.format(op.domain)) if not (isinstance(op.range, FnBase) or - isinstance(op.range, ProductSpace)): - print('WARNING: The operator range is not discrete; cannot produce', - 'matrix representation of it.') + (isinstance(op.range, ProductSpace) and + all(isinstance(spc, FnBase) for spc in op.range))): + raise TypeError('Operator range {} is not FnBase, nor ProductSpace' + 'with only FnBase components'.format(op.range)) # Get the size of the range, and handle ProductSpace - op_ran = op.range - op_ran_is_prod_space = isinstance(op_ran, ProductSpace) + # Store for reuse in loop + op_ran_is_prod_space = isinstance(op.range, ProductSpace) if op_ran_is_prod_space: - num_ran = op_ran.size - n = [op_ran[i].size for i in range(num_ran)] + num_ran = op.range.size + n = [ran.size for ran in op.range] else: num_ran = 1 - n = [op_ran.size] + n = [op.range.size] # Get the size of the domain, and handle ProductSpace - op_dom = op.domain - op_dom_is_prod_space = isinstance(op_dom, ProductSpace) + # Store for reuse in loop + op_dom_is_prod_space = isinstance(op.domain, ProductSpace) if op_dom_is_prod_space: - num_dom = op_dom.size - m = [op_dom[i].size for i in range(num_dom)] + num_dom = op.domain.size + m = [dom.size for dom in op.domain] else: num_dom = 1 - m = [op_dom.size] + m = [op.domain.size] # Generate the matrix matrix = np.zeros([np.sum(n), np.sum(m)]) - tmp_1 = op_ran.element() - v = op_dom.element() + tmp_ran = op.range.element() # Store for reuse in loop + tmp_dom = op.domain.element() # Store for reuse in loop index = 0 for i in range(num_dom): for j in range(m[i]): - v.set_zero() + tmp_dom.set_zero() if op_dom_is_prod_space: - v[i][j] = 1.0 + tmp_dom[i][j] = 1.0 else: - v[j] = 1.0 - tmp_2 = op(v, out=tmp_1) - tmp_3 = [] + tmp_dom[j] = 1.0 + op(tmp_dom, out=tmp_ran) if op_ran_is_prod_space: + tmp_result = np.empty(np.sum(n)) + tmp_idx = 0 for k in range(num_ran): - tmp_3 = np.concatenate((tmp_3, tmp_2[k].asarray())) + tmp_result[tmp_idx: tmp_idx + op.range[k].size] = \ + tmp_ran[k].asarray() + tmp_idx += op.range[k].size else: - tmp_3 = tmp_2.asarray() - matrix[:, index] = tmp_3 + tmp_result = tmp_ran.asarray() + matrix[:, index] = tmp_result index += 1 return matrix From 06ba7de8f66203fe6d33ed48a971737d7c8e5d46 Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 09:45:14 +0100 Subject: [PATCH 20/26] TST; Update of tests according to previous commit. New tests checking that the code raises correct errors. See #73 and #49. --- test/operator/operator_utilities_test.py | 53 ++++++++++++++++++++++++ 1 file changed, 53 insertions(+) diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index 5225bc1d6f1..f0fc323be2a 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -142,5 +142,58 @@ def test_matrix_representation_product_to_product_two(): assert almost_equal(np.sum(np.abs(AB_matrix - the_matrix)), 1e-6) + +def test_matrix_representation_not_linear_op(): + # Verify that the matrix representation function gives correct error + class small_nonlin_op(odl.Operator): + """Small nonlinear test operator""" + def __init__(self): + super().__init__(domain=odl.Rn(3), range=odl.Rn(4), linear=False) + + def _apply(self, x, out): + return odl.Rn(np.random.rand(4)) + + nonlin_op = small_nonlin_op() + with pytest.raises(ValueError): + matrix_representation(nonlin_op) + + +def test_matrix_representation_wrong_domain(): + # Verify that the matrix representation function gives correct error + class small_op(odl.Operator): + """Small nonlinear test operator""" + def __init__(self): + super().__init__(domain=ProductSpace(odl.Rn(3), + ProductSpace(odl.Rn(3), + odl.Rn(3))), + range=odl.Rn(4), linear=True) + + def _apply(self, x, out): + return odl.Rn(np.random.rand(4)) + + nonlin_op = small_op() + with pytest.raises(TypeError): + matrix_representation(nonlin_op) + + +def test_matrix_representation_wrong_range(): + # Verify that the matrix representation function gives correct error + class small_op(odl.Operator): + """Small nonlinear test operator""" + def __init__(self): + super().__init__(domain=odl.Rn(3), + range=ProductSpace(odl.Rn(3), + ProductSpace(odl.Rn(3), + odl.Rn(3))), + linear=True) + + def _apply(self, x, out): + return odl.Rn(np.random.rand(4)) + + nonlin_op = small_op() + with pytest.raises(TypeError): + matrix_representation(nonlin_op) + + if __name__ == '__main__': pytest.main(str(__file__.replace('\\', '/')) + ' -v') From 607c4d4873cdf02d5dffc778de00cbeacde3c035 Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 10:01:37 +0100 Subject: [PATCH 21/26] TST: Import super from builtins. Done to not break the build in python 2, see http://odl.readthedocs.org/guide/faq.html#errors-related-to-python-2-3. Issues #49 and #73. --- test/operator/operator_utilities_test.py | 1 + 1 file changed, 1 insertion(+) diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index f0fc323be2a..b1d0bd5bc16 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -20,6 +20,7 @@ from __future__ import print_function, division, absolute_import from future import standard_library standard_library.install_aliases() +from builtins import super # External module imports import pytest From 872df47157d3157020bdc136df958268a5ddd8f1 Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 13:23:24 +0100 Subject: [PATCH 22/26] TST: Rebase and use MatVecOperator according to suggestion. The rebase was needed since commit e98cdc0 changed the interface of it. Moreover, explicit spaces where removed from the call, since this is easily inferred from the matrix dimensions. --- test/operator/operator_test.py | 35 +++++++----------------- test/operator/operator_utilities_test.py | 25 +++++++++-------- 2 files changed, 23 insertions(+), 37 deletions(-) diff --git a/test/operator/operator_test.py b/test/operator/operator_test.py index 3a7e5a274df..102f4b5f35b 100644 --- a/test/operator/operator_test.py +++ b/test/operator/operator_test.py @@ -191,9 +191,7 @@ def test_linear_Op(): x = np.random.rand(3) out = np.random.rand(3) - r3 = odl.Rn(3) - - Aop = MatVecOperator(r3, r3, A) + Aop = MatVecOperator(A) xvec = Aop.domain.element(x) outvec = Aop.range.element() @@ -212,9 +210,7 @@ def test_linear_op_nonsquare(): x = np.random.rand(3) out = np.random.rand(4) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r4, A) + Aop = MatVecOperator(A) xvec = Aop.domain.element(x) outvec = Aop.range.element() @@ -233,9 +229,7 @@ def test_linear_adjoint(): x = np.random.rand(4) out = np.random.rand(3) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r4, A) + Aop = MatVecOperator(A) xvec = Aop.range.element(x) outvec = Aop.domain.element() @@ -254,10 +248,8 @@ def test_linear_addition(): x = np.random.rand(3) y = np.random.rand(4) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r4, A) - Bop = MatVecOperator(r3, r4, B) + Aop = MatVecOperator(A) + Bop = MatVecOperator(B) xvec = Aop.domain.element(x) yvec = Aop.range.element(y) @@ -282,9 +274,7 @@ def test_linear_scale(): x = np.random.rand(3) y = np.random.rand(4) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r4, A) + Aop = MatVecOperator(A) xvec = Aop.domain.element(x) yvec = Aop.range.element(y) @@ -314,9 +304,7 @@ def test_linear_scale(): def test_linear_right_vector_mult(): A = np.random.rand(4, 3) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r4, A) + Aop = MatVecOperator(A) vec = Aop.domain.element([1, 2, 3]) x = Aop.domain.element([4, 5, 6]) y = Aop.range.element([5, 6, 7, 8]) @@ -345,11 +333,8 @@ def test_linear_composition(): x = np.random.rand(3) y = np.random.rand(5) - r3 = odl.Rn(3) - r4 = odl.Rn(4) - r5 = odl.Rn(5) - Aop = MatVecOperator(r4, r5, A) - Bop = MatVecOperator(r3, r4, B) + Aop = MatVecOperator(A) + Bop = MatVecOperator(B) xvec = Bop.domain.element(x) yvec = Aop.range.element(y) @@ -366,7 +351,7 @@ def test_type_errors(): r3 = odl.Rn(3) r4 = odl.Rn(4) - Aop = MatVecOperator(r3, r3, np.random.rand(3, 3)) + Aop = MatVecOperator(np.random.rand(3, 3)) r3Vec1 = r3.zero() r3Vec2 = r3.zero() r4Vec1 = r4.zero() diff --git a/test/operator/operator_utilities_test.py b/test/operator/operator_utilities_test.py index b1d0bd5bc16..5603d4e0a44 100644 --- a/test/operator/operator_utilities_test.py +++ b/test/operator/operator_utilities_test.py @@ -40,10 +40,9 @@ def test_matrix_representation(): # Verify that the matrix representation function returns the correct matrix n = 3 - rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) + Aop = MatVecOperator(A) the_matrix = matrix_representation(Aop) @@ -56,12 +55,12 @@ def test_matrix_representation_product_to_lin_space(): n = 3 rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) + Aop = MatVecOperator(A) m = 2 rm = odl.Rn(m) B = np.random.rand(n, m) - Bop = MatVecOperator(rm, rn, B) + Bop = MatVecOperator(B) dom = ProductSpace(rn, rm) ran = ProductSpace(rn, 1) @@ -80,12 +79,12 @@ def test_matrix_representation_lin_space_to_product(): n = 3 rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) + Aop = MatVecOperator(A) m = 2 rm = odl.Rn(m) B = np.random.rand(m, n) - Bop = MatVecOperator(rn, rm, B) + Bop = MatVecOperator(B) dom = ProductSpace(rn, 1) ran = ProductSpace(rn, rm) @@ -104,18 +103,19 @@ def test_matrix_representation_product_to_product(): n = 3 rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) + Aop = MatVecOperator(A) m = 2 rm = odl.Rn(m) B = np.random.rand(m, m) - Bop = MatVecOperator(rm, rm, B) + Bop = MatVecOperator(B) ran_and_dom = ProductSpace(rn, rm) AB_matrix = np.vstack([np.hstack([A, np.zeros((n, m))]), np.hstack([np.zeros((m, n)), B])]) - ABop = ProductSpaceOperator([[Aop, None], [None, Bop]], + ABop = ProductSpaceOperator([[Aop, 0], + [0, Bop]], ran_and_dom, ran_and_dom) the_matrix = matrix_representation(ABop) @@ -128,16 +128,17 @@ def test_matrix_representation_product_to_product_two(): n = 3 rn = odl.Rn(n) A = np.random.rand(n, n) - Aop = MatVecOperator(rn, rn, A) + Aop = MatVecOperator(A) B = np.random.rand(n, n) - Bop = MatVecOperator(rn, rn, B) + Bop = MatVecOperator(B) ran_and_dom = ProductSpace(rn, 2) AB_matrix = np.vstack([np.hstack([A, np.zeros((n, n))]), np.hstack([np.zeros((n, n)), B])]) - ABop = ProductSpaceOperator([[Aop, None], [None, Bop]], + ABop = ProductSpaceOperator([[Aop, 0], + [0, Bop]], ran_and_dom, ran_and_dom) the_matrix = matrix_representation(ABop) From 9ffebf7e8abec6f742d3697f25b4f08c7a806bbe Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 14:35:28 +0100 Subject: [PATCH 23/26] ENH: Minor change in print when raising exception. See #49, #73. --- odl/operator/operator_utilities.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/odl/operator/operator_utilities.py b/odl/operator/operator_utilities.py index e8964716054..0c43b951e91 100644 --- a/odl/operator/operator_utilities.py +++ b/odl/operator/operator_utilities.py @@ -58,13 +58,13 @@ def matrix_representation(op): if not (isinstance(op.domain, FnBase) or (isinstance(op.domain, ProductSpace) and all(isinstance(spc, FnBase) for spc in op.domain))): - raise TypeError('Operator domain {} is not FnBase, nor ProductSpace' + raise TypeError('Operator domain {} is not FnBase, nor ProductSpace ' 'with only FnBase components'.format(op.domain)) if not (isinstance(op.range, FnBase) or (isinstance(op.range, ProductSpace) and all(isinstance(spc, FnBase) for spc in op.range))): - raise TypeError('Operator range {} is not FnBase, nor ProductSpace' + raise TypeError('Operator range {} is not FnBase, nor ProductSpace ' 'with only FnBase components'.format(op.range)) # Get the size of the range, and handle ProductSpace From da1e5d1050ca79ec354062e4b7f3049f4f2783fd Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 14:40:54 +0100 Subject: [PATCH 24/26] ENH: Change name on files. See #49, #73. --- odl/operator/{operator_utilities.py => oputils.py} | 0 test/operator/{operator_utilities_test.py => oputils_test.py} | 0 2 files changed, 0 insertions(+), 0 deletions(-) rename odl/operator/{operator_utilities.py => oputils.py} (100%) rename test/operator/{operator_utilities_test.py => oputils_test.py} (100%) diff --git a/odl/operator/operator_utilities.py b/odl/operator/oputils.py similarity index 100% rename from odl/operator/operator_utilities.py rename to odl/operator/oputils.py diff --git a/test/operator/operator_utilities_test.py b/test/operator/oputils_test.py similarity index 100% rename from test/operator/operator_utilities_test.py rename to test/operator/oputils_test.py From 1c678f02d24ffdd4700c446537ef21cf1d5a9c25 Mon Sep 17 00:00:00 2001 From: aringh Date: Thu, 26 Nov 2015 14:42:50 +0100 Subject: [PATCH 25/26] TST: Change test according to new name in commit da1e5d1. See #49, #73. --- test/operator/oputils_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/operator/oputils_test.py b/test/operator/oputils_test.py index 5603d4e0a44..cba3bb451b7 100644 --- a/test/operator/oputils_test.py +++ b/test/operator/oputils_test.py @@ -28,7 +28,7 @@ # ODL imports import odl -from odl.operator.operator_utilities import matrix_representation +from odl.operator.oputils import matrix_representation from odl.set.pspace import ProductSpace from odl.operator.pspace_ops import ProductSpaceOperator From 38ba1a610eb4d77b1a727bcafe3155db930b6b40 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 26 Nov 2015 15:13:36 +0100 Subject: [PATCH 26/26] ENH: optimize usage of tmps in matrix_representation --- odl/operator/oputils.py | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/odl/operator/oputils.py b/odl/operator/oputils.py index 0c43b951e91..78019585b88 100644 --- a/odl/operator/oputils.py +++ b/odl/operator/oputils.py @@ -90,27 +90,29 @@ def matrix_representation(op): # Generate the matrix matrix = np.zeros([np.sum(n), np.sum(m)]) tmp_ran = op.range.element() # Store for reuse in loop - tmp_dom = op.domain.element() # Store for reuse in loop + tmp_dom = op.domain.zero() # Store for reuse in loop index = 0 + last_i = last_j = 0 for i in range(num_dom): for j in range(m[i]): - tmp_dom.set_zero() if op_dom_is_prod_space: + tmp_dom[last_i][last_j] = 0.0 tmp_dom[i][j] = 1.0 else: + tmp_dom[last_j] = 0.0 tmp_dom[j] = 1.0 op(tmp_dom, out=tmp_ran) if op_ran_is_prod_space: - tmp_result = np.empty(np.sum(n)) tmp_idx = 0 for k in range(num_ran): - tmp_result[tmp_idx: tmp_idx + op.range[k].size] = \ - tmp_ran[k].asarray() + matrix[tmp_idx: tmp_idx + op.range[k].size, index] = ( + tmp_ran[k]) tmp_idx += op.range[k].size else: - tmp_result = tmp_ran.asarray() - matrix[:, index] = tmp_result + matrix[:, index] = tmp_ran.asarray() index += 1 + last_j = j + last_i = i return matrix