From 14d2ab57424362a32827d572d55f24a083493932 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sat, 21 Oct 2017 23:13:12 +0200 Subject: [PATCH 1/3] BUG: fix in-out aliasing bug in OperatorSum and OperatorPointwiseProduct Closes #1181 --- odl/operator/operator.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/odl/operator/operator.py b/odl/operator/operator.py index 3c0e01fa00d..f6bc6a3edbb 100644 --- a/odl/operator/operator.py +++ b/odl/operator/operator.py @@ -1154,8 +1154,10 @@ def _call(self, x, out=None): else: tmp = (self.__tmp_ran if self.__tmp_ran is not None else self.range.element()) - self.left(x, out=out) - self.right(x, out=tmp) + # Write to `tmp` first, otherwise aliased `x` and `out` lead + # to wrong result + self.left(x, out=tmp) + self.right(x, out=out) out += tmp def derivative(self, x): @@ -1492,8 +1494,10 @@ def _call(self, x, out=None): return self.left(x) * self.right(x) else: tmp = self.right.range.element() - self.left(x, out=out) - self.right(x, out=tmp) + # Write to `tmp` first, otherwise aliased `x` and `out` lead + # to wrong result + self.left(x, out=tmp) + self.right(x, out=out) out *= tmp def derivative(self, x): From 11521f26e771c7ef078e5146b3ffe4f0b8ee1fc7 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sat, 11 Nov 2017 02:15:21 +0100 Subject: [PATCH 2/3] TST: test operators with aliased in and out also --- odl/test/operator/operator_test.py | 695 +++++++++++++++-------------- 1 file changed, 370 insertions(+), 325 deletions(-) diff --git a/odl/test/operator/operator_test.py b/odl/test/operator/operator_test.py index 5d498ea4672..0b19b6362bb 100644 --- a/odl/test/operator/operator_test.py +++ b/odl/test/operator/operator_test.py @@ -18,11 +18,24 @@ MatrixOperator, OperatorLeftVectorMult, OpTypeError, OpDomainError, OpRangeError) from odl.operator.operator import _signature_from_spec, _dispatch_call_args -from odl.util.testutils import almost_equal, all_almost_equal, noise_element +from odl.util.testutils import ( + almost_equal, all_almost_equal, noise_element, noise_elements, + simple_fixture) + + +# --- Fixtures --- # + + +# This fixture is intended to activate testing of operator evaluation +# with aliased input and output, which is only possible if domain == range. +dom_eq_ran = simple_fixture('dom_eq_ran', [True, False]) + + +# --- Auxilliary --- # class MultiplyAndSquareOp(Operator): - """Example of a nonlinear operator, x --> (A*x)**2.""" + """Example of a nonlinear operator, x --> (mat*x)**2.""" def __init__(self, matrix, domain=None, range=None): dom = (odl.rn(matrix.shape[1]) @@ -33,12 +46,11 @@ def __init__(self, matrix, domain=None, range=None): super(MultiplyAndSquareOp, self).__init__(dom, ran) self.matrix = matrix - def _call(self, rhs, out=None): + def _call(self, x, out=None): if out is None: - return np.dot(self.matrix, rhs.data) ** 2 - else: - np.dot(self.matrix, rhs.data, out=out.data) - out **= 2 + out = self.range.element() + out[:] = np.dot(self.matrix, x.data) + out **= 2 def derivative(self, x): return 2 * odl.MatrixOperator(self.matrix) @@ -47,33 +59,46 @@ def __str__(self): return "MaS: " + str(self.matrix) + " ** 2" -def mult_sq_np(A, x): - # The same as MultiplyAndSquareOp but only using numpy - return np.dot(A, x) ** 2 +def mult_sq_np(mat, x): + """NumPy reference implementation of MultiplyAndSquareOp.""" + return np.dot(mat, x) ** 2 -def test_nonlinear_op(): - # Verify that the operator does indeed work as expected - A = np.random.rand(4, 3) - x = np.random.rand(3) - Aop = MultiplyAndSquareOp(A) - xvec = Aop.domain.element(x) +def check_call(operator, x, expected): + """Assert that operator(point) == expected.""" + # Out-of-place check + assert all_almost_equal(operator(x), expected) - assert all_almost_equal(Aop(xvec), mult_sq_np(A, x)) + # In-place check, no aliasing + out = operator.range.element() + operator(x, out=out) + assert all_almost_equal(out, expected) + # In-place check, aliased + if operator.domain == operator.range: + y = x.copy() + operator(y, out=y) + assert all_almost_equal(y, expected) -def check_call(operator, point, expected): - """Assert that operator(point) == expected.""" - assert all_almost_equal(operator(point), expected) +# --- Unit tests --- # - out = operator.range.element() - operator(point, out=out) - assert all_almost_equal(out, expected) +def test_operator_call(dom_eq_ran): + """Check operator evaluation against NumPy reference.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + op = MultiplyAndSquareOp(mat) + assert op.domain == op.range + else: + mat = np.random.rand(4, 3) + op = MultiplyAndSquareOp(mat) + xarr, x = noise_elements(op.domain) + assert all_almost_equal(op(x), mult_sq_np(mat, xarr)) -def test_call_in_place_wrong_return(): + +def test_operator_call_in_place_wrong_return(): """Test that operator with out parameter actually returns out.""" class BadInplaceOperator(odl.Operator): def _call(self, x, out): @@ -92,386 +117,406 @@ def _call(self, x, out): op(space.zero(), out=out) -def test_nonlinear_addition(): - # Test operator addition - A = np.random.rand(4, 3) - B = np.random.rand(4, 3) - x = np.random.rand(3) +def test_operator_sum(dom_eq_ran): + """Check operator sum against NumPy reference.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(4, 3) + mat2 = np.random.rand(4, 3) - Aop = MultiplyAndSquareOp(A) - Bop = MultiplyAndSquareOp(B) - xvec = Aop.domain.element(x) + op1 = MultiplyAndSquareOp(mat1) + op2 = MultiplyAndSquareOp(mat2) + xarr, x = noise_elements(op1.domain) # Explicit instantiation - C = OperatorSum(Aop, Bop) - - assert not C.is_linear - - check_call(C, xvec, mult_sq_np(A, x) + mult_sq_np(B, x)) + sum_op = OperatorSum(op1, op2) + assert not sum_op.is_linear + check_call(sum_op, x, mult_sq_np(mat1, xarr) + mult_sq_np(mat2, xarr)) # Using operator overloading - check_call(Aop + Bop, xvec, mult_sq_np(A, x) + mult_sq_np(B, x)) + check_call(op1 + op2, x, mult_sq_np(mat1, xarr) + mult_sq_np(mat2, xarr)) - # Verify that unmatched operators domains fail - C = np.random.rand(4, 4) - Cop = MultiplyAndSquareOp(C) + # Verify that unmatched operator domains fail + op_wrong_dom = MultiplyAndSquareOp(mat1[:, :-1]) + with pytest.raises(OpTypeError): + OperatorSum(op1, op_wrong_dom) + # Verify that unmatched operator ranges fail + op_wrong_ran = MultiplyAndSquareOp(mat1[:-1, :]) with pytest.raises(OpTypeError): - C = OperatorSum(Aop, Cop) + OperatorSum(op1, op_wrong_ran) -def test_nonlinear_scale(): - A = np.random.rand(4, 3) - x = np.random.rand(3) +def test_operator_scaling(dom_eq_ran): + """Check operator scaling against NumPy reference.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) - Aop = MultiplyAndSquareOp(A) - xvec = Aop.domain.element(x) + op = MultiplyAndSquareOp(mat) + xarr, x = noise_elements(op.domain) # Test a range of scalars (scalar multiplication could implement # optimizations for (-1, 0, 1)). scalars = [-1.432, -1, 0, 1, 3.14] - for scale in scalars: - lscaled = OperatorLeftScalarMult(Aop, scale) - rscaled = OperatorRightScalarMult(Aop, scale) + for scalar in scalars: + lscaled = OperatorLeftScalarMult(op, scalar) + rscaled = OperatorRightScalarMult(op, scalar) assert not lscaled.is_linear assert not rscaled.is_linear # Explicit - check_call(lscaled, xvec, scale * mult_sq_np(A, x)) - check_call(rscaled, xvec, mult_sq_np(A, scale * x)) + check_call(lscaled, x, scalar * mult_sq_np(mat, xarr)) + check_call(rscaled, x, mult_sq_np(mat, scalar * xarr)) # Using operator overloading - check_call(scale * Aop, xvec, scale * mult_sq_np(A, x)) - check_call(Aop * scale, xvec, mult_sq_np(A, scale * x)) + check_call(scalar * op, x, scalar * mult_sq_np(mat, xarr)) + check_call(op * scalar, x, mult_sq_np(mat, scalar * xarr)) - # Fail when scaling by wrong scalar type (A complex number) + # Fail when scaling by wrong scalar type (complex number) wrongscalars = [1j, [1, 2], (1, 2)] for wrongscalar in wrongscalars: with pytest.raises(TypeError): - OperatorLeftScalarMult(Aop, wrongscalar) + OperatorLeftScalarMult(op, wrongscalar) with pytest.raises(TypeError): - OperatorRightScalarMult(Aop, wrongscalar) + OperatorRightScalarMult(op, wrongscalar) with pytest.raises(TypeError): - Aop * wrongscalar + op * wrongscalar with pytest.raises(TypeError): - wrongscalar * Aop - + wrongscalar * op -def test_nonlinear_vector_mult(): - A = np.random.rand(4, 3) - Aop = MultiplyAndSquareOp(A) - rvec = Aop.domain.element([1, 2, 3]) - lvec = Aop.range.element([1, 2, 3, 4]) - x = Aop.domain.element([4, 5, 6]) +def test_operator_vector_mult(dom_eq_ran): + """Check operator-vector multiplication against NumPy reference.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) - # Test a range of scalars (scalar multiplication could implement - # optimizations for (-1, 0, 1). - rmult = OperatorRightVectorMult(Aop, rvec) - lmult = OperatorLeftVectorMult(Aop, lvec) + op = MultiplyAndSquareOp(mat) + right = op.domain.element(np.arange(op.domain.size)) + left = op.range.element(np.arange(op.range.size)) + xarr, x = noise_elements(op.domain) - assert not rmult.is_linear - assert not lmult.is_linear + rmult_op = OperatorRightVectorMult(op, right) + lmult_op = OperatorLeftVectorMult(op, left) + assert not rmult_op.is_linear + assert not lmult_op.is_linear - check_call(rmult, x, mult_sq_np(A, rvec * x)) - check_call(lmult, x, lvec * mult_sq_np(A, x)) + check_call(rmult_op, x, mult_sq_np(mat, right * xarr)) + check_call(lmult_op, x, left * mult_sq_np(mat, xarr)) # Using operator overloading - check_call(Aop * rvec, x, mult_sq_np(A, rvec * x)) - check_call(lvec * Aop, x, lvec * mult_sq_np(A, x)) + check_call(op * right, x, mult_sq_np(mat, right * xarr)) + check_call(left * op, x, left * mult_sq_np(mat, xarr)) -def test_nonlinear_composition(): - A = np.random.rand(5, 4) - B = np.random.rand(4, 3) - x = np.random.rand(3) - - Aop = MultiplyAndSquareOp(A) - Bop = MultiplyAndSquareOp(B) - xvec = Bop.domain.element(x) +def test_operator_composition(dom_eq_ran): + """Check operator composition against NumPy reference.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(5, 4) + mat2 = np.random.rand(4, 3) - C = OperatorComp(Aop, Bop) + op1 = MultiplyAndSquareOp(mat1) + op2 = MultiplyAndSquareOp(mat2) + xarr, x = noise_elements(op2.domain) - assert not C.is_linear - - check_call(C, xvec, mult_sq_np(A, mult_sq_np(B, x))) + comp_op = OperatorComp(op1, op2) + assert not comp_op.is_linear + check_call(comp_op, x, mult_sq_np(mat1, mult_sq_np(mat2, xarr))) # Verify that incorrect order fails - with pytest.raises(OpTypeError): - C = OperatorComp(Bop, Aop) - - -def test_linear_Op(): - # Verify that the multiply op does indeed work as expected - - A = np.random.rand(3, 3) - x = np.random.rand(3) - out = np.random.rand(3) - - Aop = MatrixOperator(A) - xvec = Aop.domain.element(x) - outvec = Aop.range.element() - - # Using out parameter - Aop(xvec, outvec) - np.dot(A, x, out) - assert all_almost_equal(out, outvec) - - # Using return value - assert all_almost_equal(Aop(xvec), np.dot(A, x)) - - -def test_linear_op_nonsquare(): - # Verify that the multiply op does indeed work as expected - A = np.random.rand(4, 3) - x = np.random.rand(3) - out = np.random.rand(4) - - Aop = MatrixOperator(A) + if not dom_eq_ran: + with pytest.raises(OpTypeError): + OperatorComp(op2, op1) - xvec = Aop.domain.element(x) - outvec = Aop.range.element() - # Using out parameter - Aop(xvec, outvec) - np.dot(A, x, out) - assert all_almost_equal(out, outvec) - - # Using return value - assert all_almost_equal(Aop(xvec), np.dot(A, x)) +def test_linear_operator_call(dom_eq_ran): + """Check call of a linear operator against NumPy, and ``is_linear``.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) + op = MatrixOperator(mat) + assert op.is_linear -def test_linear_adjoint(): - A = np.random.rand(4, 3) - x = np.random.rand(4) - out = np.random.rand(3) + xarr, x = noise_elements(op.domain) + check_call(op, x, np.dot(mat, xarr)) - Aop = MatrixOperator(A) - xvec = Aop.range.element(x) - outvec = Aop.domain.element() - # Using in-place adjoint - Aop.adjoint(xvec, outvec) - np.dot(A.T, x, out) - assert all_almost_equal(out, outvec) +def test_linear_operator_adjoint(dom_eq_ran): + """Check adjoint of a linear operator against NumPy.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) - # Using out-of-place method - assert all_almost_equal(Aop.adjoint(xvec), np.dot(A.T, x)) + op = MatrixOperator(mat) + xarr, x = noise_elements(op.range) + check_call(op.adjoint, x, np.dot(mat.T, xarr)) -def test_linear_addition(): - A = np.random.rand(4, 3) - B = np.random.rand(4, 3) - x = np.random.rand(3) - y = np.random.rand(4) +def test_linear_operator_addition(dom_eq_ran): + """Check call and adjoint of a sum of linear operators.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(4, 3) + mat2 = np.random.rand(4, 3) - Aop = MatrixOperator(A) - Bop = MatrixOperator(B) - xvec = Aop.domain.element(x) - yvec = Aop.range.element(y) + op1 = MatrixOperator(mat1) + op2 = MatrixOperator(mat2) + xarr, x = noise_elements(op1.domain) + yarr, y = noise_elements(op1.range) # Explicit instantiation - C = OperatorSum(Aop, Bop) - - assert C.is_linear - assert C.adjoint.is_linear - - assert all_almost_equal(C(xvec), np.dot(A, x) + np.dot(B, x)) - assert all_almost_equal(C.adjoint(yvec), np.dot(A.T, y) + np.dot(B.T, y)) + sum_op = OperatorSum(op1, op2) + assert sum_op.is_linear + assert sum_op.adjoint.is_linear + check_call(sum_op, x, np.dot(mat1, xarr) + np.dot(mat2, xarr)) + check_call(sum_op.adjoint, y, np.dot(mat1.T, yarr) + np.dot(mat2.T, yarr)) # Using operator overloading - assert all_almost_equal((Aop + Bop)(xvec), - np.dot(A, x) + np.dot(B, x)) - assert all_almost_equal((Aop + Bop).adjoint(yvec), - np.dot(A.T, y) + np.dot(B.T, y)) + check_call(op1 + op2, x, np.dot(mat1, xarr) + np.dot(mat2, xarr)) + check_call((op1 + op2).adjoint, + y, np.dot(mat1.T, yarr) + np.dot(mat2.T, yarr)) -def test_linear_scale(): - A = np.random.rand(4, 3) - x = np.random.rand(3) - y = np.random.rand(4) +def test_linear_operator_scaling(dom_eq_ran): + """Check call and adjoint of a scaled linear operator.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) - Aop = MatrixOperator(A) - xvec = Aop.domain.element(x) - yvec = Aop.range.element(y) + op = MatrixOperator(mat) + xarr, x = noise_elements(op.domain) + yarr, y = noise_elements(op.range) # Test a range of scalars (scalar multiplication could implement # optimizations for (-1, 0, 1). scalars = [-1.432, -1, 0, 1, 3.14] - for scale in scalars: - C = OperatorRightScalarMult(Aop, scale) + for scalar in scalars: + # Explicit instantiation + scaled_op = OperatorRightScalarMult(op, scalar) + assert scaled_op.is_linear + assert scaled_op.adjoint.is_linear + check_call(scaled_op, x, scalar * np.dot(mat, xarr)) + check_call(scaled_op.adjoint, y, scalar * np.dot(mat.T, yarr)) - assert C.is_linear - assert C.adjoint.is_linear + # Using operator overloading + check_call(scalar * op, x, scalar * np.dot(mat, xarr)) + check_call(op * scalar, x, scalar * np.dot(mat, xarr)) + check_call((scalar * op).adjoint, y, scalar * np.dot(mat.T, yarr)) + check_call((op * scalar).adjoint, y, scalar * np.dot(mat.T, yarr)) - assert all_almost_equal(C(xvec), scale * np.dot(A, x)) - assert all_almost_equal(C.adjoint(yvec), scale * np.dot(A.T, y)) - # Using operator overloading - assert all_almost_equal((scale * Aop)(xvec), - scale * np.dot(A, x)) - assert all_almost_equal((Aop * scale)(xvec), - np.dot(A, scale * x)) - assert all_almost_equal((scale * Aop).adjoint(yvec), - scale * np.dot(A.T, y)) - assert all_almost_equal((Aop * scale).adjoint(yvec), - np.dot(A.T, scale * y)) +def test_linear_right_vector_mult(dom_eq_ran): + """Check call and adjoint of linear operator x vector.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) + op = MatrixOperator(mat) + (xarr, mul_arr), (x, mul) = noise_elements(op.domain, n=2) + yarr, y = noise_elements(op.range) -def test_linear_right_vector_mult(): - A = np.random.rand(4, 3) + # Explicit instantiation + rmult_op = OperatorRightVectorMult(op, mul) + assert rmult_op.is_linear + assert rmult_op.adjoint.is_linear + check_call(rmult_op, x, np.dot(mat, mul_arr * xarr)) + check_call(rmult_op.adjoint, y, mul_arr * np.dot(mat.T, yarr)) - Aop = MatrixOperator(A) - vec = Aop.domain.element([1, 2, 3]) - x = Aop.domain.element([4, 5, 6]) - y = Aop.range.element([5, 6, 7, 8]) + # Using operator overloading + check_call(op * mul, x, np.dot(mat, mul_arr * xarr)) + check_call((op * mul).adjoint, y, mul_arr * np.dot(mat.T, yarr)) - # Test a range of scalars (scalar multiplication could implement - # optimizations for (-1, 0, 1). - C = OperatorRightVectorMult(Aop, vec) - assert C.is_linear - assert C.adjoint.is_linear +def test_linear_left_vector_mult(dom_eq_ran): + """Check call and adjoint of vector x linear operator.""" + if dom_eq_ran: + mat = np.random.rand(3, 3) + else: + mat = np.random.rand(4, 3) - assert all_almost_equal(C(x), np.dot(A, vec * x)) - assert all_almost_equal(C.adjoint(y), vec * np.dot(A.T, y)) - assert all_almost_equal(C.adjoint.adjoint(x), C(x)) + op = MatrixOperator(mat) + xarr, x = noise_elements(op.domain) + (yarr, mul_arr), (y, mul) = noise_elements(op.range, n=2) - # Using operator overloading - assert all_almost_equal((Aop * vec)(x), - np.dot(A, vec * x)) - assert all_almost_equal((Aop * vec).adjoint(y), - vec * np.dot(A.T, y)) + # Explicit instantiation + lmult_op = OperatorLeftVectorMult(op, mul) + assert lmult_op.is_linear + assert lmult_op.adjoint.is_linear + check_call(lmult_op, x, mul_arr * np.dot(mat, xarr)) + check_call(lmult_op.adjoint, y, np.dot(mat.T, mul_arr * yarr)) + # Using operator overloading + check_call(mul * op, x, mul_arr * np.dot(mat, xarr)) + check_call((mul * op).adjoint, y, np.dot(mat.T, mul_arr * yarr)) -def test_linear_composition(): - A = np.random.rand(5, 4) - B = np.random.rand(4, 3) - x = np.random.rand(3) - y = np.random.rand(5) - Aop = MatrixOperator(A) - Bop = MatrixOperator(B) - xvec = Bop.domain.element(x) - yvec = Aop.range.element(y) +def test_linear_operator_composition(dom_eq_ran): + """Check call and adjoint of linear operator composition.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(4, 3) + mat2 = np.random.rand(3, 4) - C = OperatorComp(Aop, Bop) + op1 = MatrixOperator(mat1) + op2 = MatrixOperator(mat2) + xarr, x = noise_elements(op2.domain) + yarr, y = noise_elements(op1.range) - assert C.is_linear - assert C.adjoint.is_linear + # Explicit instantiation + comp_op = OperatorComp(op1, op2) + assert comp_op.is_linear + assert comp_op.adjoint.is_linear + check_call(comp_op, x, np.dot(mat1, np.dot(mat2, xarr))) + check_call(comp_op.adjoint, y, np.dot(mat2.T, np.dot(mat1.T, yarr))) - assert all_almost_equal(C(xvec), np.dot(A, np.dot(B, x))) - assert all_almost_equal(C.adjoint(yvec), np.dot(B.T, np.dot(A.T, y))) + # Using operator overloading + check_call(op1 * op2, x, np.dot(mat1, np.dot(mat2, xarr))) + check_call((op1 * op2).adjoint, y, np.dot(mat2.T, np.dot(mat1.T, yarr))) def test_type_errors(): r3 = odl.rn(3) r4 = odl.rn(4) - Aop = MatrixOperator(np.random.rand(3, 3)) - r3Vec1 = r3.zero() - r3Vec2 = r3.zero() - r4Vec1 = r4.zero() - r4Vec2 = r4.zero() + op = MatrixOperator(np.random.rand(3, 3)) + r3_elem1 = r3.zero() + r3_elem2 = r3.zero() + r4_elem1 = r4.zero() + r4_elem2 = r4.zero() # Verify that correct usage works - Aop(r3Vec1, r3Vec2) - Aop.adjoint(r3Vec1, r3Vec2) + op(r3_elem1, r3_elem2) + op.adjoint(r3_elem1, r3_elem2) # Test that erroneous usage raises with pytest.raises(OpDomainError): - Aop(r4Vec1) + op(r4_elem1) with pytest.raises(OpDomainError): - Aop.adjoint(r4Vec1) + op.adjoint(r4_elem1) with pytest.raises(OpRangeError): - Aop(r3Vec1, r4Vec1) + op(r3_elem1, r4_elem1) with pytest.raises(OpRangeError): - Aop.adjoint(r3Vec1, r4Vec1) + op.adjoint(r3_elem1, r4_elem1) with pytest.raises(OpDomainError): - Aop(r4Vec1, r3Vec1) + op(r4_elem1, r3_elem1) with pytest.raises(OpDomainError): - Aop.adjoint(r4Vec1, r3Vec1) + op.adjoint(r4_elem1, r3_elem1) with pytest.raises(OpDomainError): - Aop(r4Vec1, r4Vec2) + op(r4_elem1, r4_elem2) with pytest.raises(OpDomainError): - Aop.adjoint(r4Vec1, r4Vec2) + op.adjoint(r4_elem1, r4_elem2) -def test_arithmetic(): +def test_arithmetic(dom_eq_ran): """Test that all standard arithmetic works.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + mat3 = np.random.rand(3, 3) + mat4 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(4, 3) + mat2 = np.random.rand(4, 3) + mat3 = np.random.rand(3, 3) + mat4 = np.random.rand(4, 4) + + op = MultiplyAndSquareOp(mat1) + op2 = MultiplyAndSquareOp(mat2) + op3 = MultiplyAndSquareOp(mat3) + op4 = MultiplyAndSquareOp(mat4) # Create elements needed for later - operator = MultiplyAndSquareOp(np.random.rand(4, 3)) - operator2 = MultiplyAndSquareOp(np.random.rand(4, 3)) - operator3 = MultiplyAndSquareOp(np.random.rand(3, 3)) - operator4 = MultiplyAndSquareOp(np.random.rand(4, 4)) - x = noise_element(operator.domain) - y = noise_element(operator.domain) - z = noise_element(operator.range) + x = noise_element(op.domain) + y = noise_element(op.domain) + z = noise_element(op.range) scalar = np.pi # Simple tests here, more in depth comes later - check_call(+operator, x, operator(x)) - check_call(-operator, x, -operator(x)) - check_call(scalar * operator, x, scalar * operator(x)) - check_call(scalar * (scalar * operator), x, scalar**2 * operator(x)) - check_call(operator * scalar, x, operator(scalar * x)) - check_call((operator * scalar) * scalar, x, operator(scalar**2 * x)) - check_call(operator + operator2, x, operator(x) + operator2(x)) - check_call(operator - operator2, x, operator(x) - operator2(x)) - check_call(operator * operator3, x, operator(operator3(x))) - check_call(operator4 * operator, x, operator4(operator(x))) - check_call(z * operator, x, z * operator(x)) - check_call(z * (z * operator), x, (z * z) * operator(x)) - check_call(operator * y, x, operator(x * y)) - check_call((operator * y) * y, x, operator((y * y) * x)) - check_call(operator + z, x, operator(x) + z) - check_call(operator - z, x, operator(x) - z) - check_call(z + operator, x, z + operator(x)) - check_call(z - operator, x, z - operator(x)) - check_call(operator + scalar, x, operator(x) + scalar) - check_call(operator - scalar, x, operator(x) - scalar) - check_call(scalar + operator, x, scalar + operator(x)) - check_call(scalar - operator, x, scalar - operator(x)) + check_call(+op, x, op(x)) + check_call(-op, x, -op(x)) + check_call(scalar * op, x, scalar * op(x)) + check_call(scalar * (scalar * op), x, scalar**2 * op(x)) + check_call(op * scalar, x, op(scalar * x)) + check_call((op * scalar) * scalar, x, op(scalar**2 * x)) + check_call(op + op2, x, op(x) + op2(x)) + check_call(op - op2, x, op(x) - op2(x)) + check_call(op * op3, x, op(op3(x))) + check_call(op4 * op, x, op4(op(x))) + check_call(z * op, x, z * op(x)) + check_call(z * (z * op), x, (z * z) * op(x)) + check_call(op * y, x, op(x * y)) + check_call((op * y) * y, x, op((y * y) * x)) + check_call(op + z, x, op(x) + z) + check_call(op - z, x, op(x) - z) + check_call(z + op, x, z + op(x)) + check_call(z - op, x, z - op(x)) + check_call(op + scalar, x, op(x) + scalar) + check_call(op - scalar, x, op(x) - scalar) + check_call(scalar + op, x, scalar + op(x)) + check_call(scalar - op, x, scalar - op(x)) def test_operator_pointwise_product(): - """Test OperatorPointwiseProduct.""" - Aop = MultiplyAndSquareOp(np.random.rand(4, 3)) - Bop = MultiplyAndSquareOp(np.random.rand(4, 3)) - x = noise_element(Aop.domain) + """Check call and adjoint of operator pointwise multiplication.""" + if dom_eq_ran: + mat1 = np.random.rand(3, 3) + mat2 = np.random.rand(3, 3) + else: + mat1 = np.random.rand(4, 3) + mat2 = np.random.rand(4, 3) - prod = odl.OperatorPointwiseProduct(Aop, Bop) + op1 = MultiplyAndSquareOp(mat1) + op2 = MultiplyAndSquareOp(mat2) + x = noise_element(op1.domain) + + prod_op = odl.OperatorPointwiseProduct(op1, op2) # Evaluate - expected = Aop(x) * Bop(x) - check_call(prod, x, expected) + expected = op1(x) * op2(x) + check_call(prod_op, x, expected) # Derivative - y = noise_element(Aop.domain) - expected = (Aop.derivative(x)(y) * Bop(x) + - Bop.derivative(x)(y) * Aop(x)) - derivative = prod.derivative(x) - assert derivative.is_linear - check_call(derivative, y, expected) - - # Adjoint - z = noise_element(Aop.range) - expected = (Aop.derivative(x).adjoint(z * Bop(x)) + - Bop.derivative(x).adjoint(z * Aop(x))) - adjoint = derivative.adjoint - assert adjoint.is_linear - check_call(adjoint, z, expected) + y = noise_element(op1.domain) + expected = (op1.derivative(x)(y) * op2(x) + + op2.derivative(x)(y) * op1(x)) + prod_deriv_op = prod_op.derivative(x) + assert prod_deriv_op.is_linear + check_call(prod_deriv_op, y, expected) + + # Derivative Adjoint + z = noise_element(op1.range) + expected = (op1.derivative(x).adjoint(z * op2(x)) + + op2.derivative(x).adjoint(z * op1(x))) + prod_deriv_adj_op = prod_deriv_op.adjoint + assert prod_deriv_adj_op.is_linear + check_call(prod_deriv_adj_op, z, expected) # FUNCTIONAL TEST @@ -540,13 +585,13 @@ def test_functional_adjoint(): def test_functional_addition(): r3 = odl.rn(3) - Aop = SumFunctional(r3) - Bop = SumFunctional(r3) + op = SumFunctional(r3) + op2 = SumFunctional(r3) x = r3.element([1, 2, 3]) y = 1 # Explicit instantiation - C = OperatorSum(Aop, Bop) + C = OperatorSum(op, op2) assert C.is_linear assert C.adjoint.is_linear @@ -558,50 +603,50 @@ def test_functional_addition(): assert all_almost_equal(C.adjoint.adjoint(x), C(x)) # Using operator overloading - assert (Aop + Bop)(x) == 2 * np.sum(x) - assert all_almost_equal((Aop + Bop).adjoint(y), y * 2 * np.ones(3)) + assert (op + op2)(x) == 2 * np.sum(x) + assert all_almost_equal((op + op2).adjoint(y), y * 2 * np.ones(3)) def test_functional_scale(): r3 = odl.rn(3) - Aop = SumFunctional(r3) + op = SumFunctional(r3) x = r3.element([1, 2, 3]) y = 1 # Test a range of scalars (scalar multiplication could implement # optimizations for (-1, 0, 1). scalars = [-1.432, -1, 0, 1, 3.14] - for scale in scalars: - C = OperatorRightScalarMult(Aop, scale) + for scalar in scalars: + C = OperatorRightScalarMult(op, scalar) assert C.is_linear assert C.adjoint.is_linear - assert C(x) == scale * np.sum(x) - assert all_almost_equal(C.adjoint(y), scale * y * np.ones(3)) + assert C(x) == scalar * np.sum(x) + assert all_almost_equal(C.adjoint(y), scalar * y * np.ones(3)) assert all_almost_equal(C.adjoint.adjoint(x), C(x)) # Using operator overloading - assert (scale * Aop)(x) == scale * np.sum(x) - assert (Aop * scale)(x) == scale * np.sum(x) - assert all_almost_equal((scale * Aop).adjoint(y), - scale * y * np.ones(3)) - assert all_almost_equal((Aop * scale).adjoint(y), - scale * y * np.ones(3)) + assert (scalar * op)(x) == scalar * np.sum(x) + assert (op * scalar)(x) == scalar * np.sum(x) + assert all_almost_equal((scalar * op).adjoint(y), + scalar * y * np.ones(3)) + assert all_almost_equal((op * scalar).adjoint(y), + scalar * y * np.ones(3)) def test_functional_left_vector_mult(): r3 = odl.rn(3) r4 = odl.rn(4) - Aop = SumFunctional(r3) + op = SumFunctional(r3) x = r3.element([1, 2, 3]) y = r4.element([3, 2, 1, 5]) # Test a range of scalars (scalar multiplication could implement # optimizations for (-1, 0, 1). - C = FunctionalLeftVectorMult(Aop, y) + C = FunctionalLeftVectorMult(op, y) assert C.is_linear assert C.adjoint.is_linear @@ -611,23 +656,23 @@ def test_functional_left_vector_mult(): assert all_almost_equal(C.adjoint.adjoint(x), C(x)) # Using operator overloading - assert all_almost_equal((y * Aop)(x), + assert all_almost_equal((y * op)(x), y * np.sum(x)) - assert all_almost_equal((y * Aop).adjoint(y), + assert all_almost_equal((y * op).adjoint(y), y.inner(y) * np.ones(3)) def test_functional_right_vector_mult(): r3 = odl.rn(3) - Aop = SumFunctional(r3) + op = SumFunctional(r3) vec = r3.element([1, 2, 3]) x = r3.element([4, 5, 6]) y = 2.0 # Test a range of scalars (scalar multiplication could implement # optimizations for (-1, 0, 1). - C = OperatorRightVectorMult(Aop, vec) + C = OperatorRightVectorMult(op, vec) assert C.is_linear assert C.adjoint.is_linear @@ -637,21 +682,21 @@ def test_functional_right_vector_mult(): assert all_almost_equal(C.adjoint.adjoint(x), C(x)) # Using operator overloading - assert all_almost_equal((Aop * vec)(x), + assert all_almost_equal((op * vec)(x), np.sum(vec * x)) - assert all_almost_equal((Aop * vec).adjoint(y), + assert all_almost_equal((op * vec).adjoint(y), vec * y) def test_functional_composition(): r3 = odl.rn(3) - Aop = SumFunctional(r3) - Bop = ConstantVector(r3) + op = SumFunctional(r3) + op2 = ConstantVector(r3) x = r3.element([1, 2, 3]) y = 1 - C = OperatorComp(Bop, Aop) + C = OperatorComp(op2, op) assert C.is_linear assert C.adjoint.is_linear @@ -661,11 +706,11 @@ def test_functional_composition(): assert all_almost_equal(C.adjoint.adjoint(x), C(x)) # Using operator overloading - assert (Aop * Bop)(y) == y * 3 - assert (Aop * Bop).adjoint(y) == y * 3 - assert all_almost_equal((Bop * Aop)(x), + assert (op * op2)(y) == y * 3 + assert (op * op2).adjoint(y) == y * 3 + assert all_almost_equal((op2 * op)(x), np.sum(x) * np.ones(3)) - assert all_almost_equal((Bop * Aop).adjoint(x), + assert all_almost_equal((op2 * op).adjoint(x), np.sum(x) * np.ones(3)) @@ -705,38 +750,38 @@ def test_nonlinear_functional_operators(): r3 = odl.rn(3) x = r3.element([1, 2, 3]) - A = SumSquaredFunctional(r3) - B = SumFunctional(r3) + mat = SumSquaredFunctional(r3) + mat2 = SumFunctional(r3) # Sum - C = A + B + C = mat + mat2 assert not C.is_linear - assert almost_equal(C(x), A(x) + B(x)) + assert almost_equal(C(x), mat(x) + mat2(x)) # Minus - C = A - B + C = mat - mat2 assert not C.is_linear - assert almost_equal(C(x), A(x) - B(x)) + assert almost_equal(C(x), mat(x) - mat2(x)) # left mul - C = 2.0 * A + C = 2.0 * mat assert not C.is_linear - assert almost_equal(C(x), 2.0 * A(x)) + assert almost_equal(C(x), 2.0 * mat(x)) # right mul - C = A * 2.0 + C = mat * 2.0 assert not C.is_linear - assert almost_equal(C(x), A(x * 2.0)) + assert almost_equal(C(x), mat(x * 2.0)) # right divide - C = A / 2.0 + C = mat / 2.0 assert not C.is_linear - assert almost_equal(C(x), A(x / 2.0)) + assert almost_equal(C(x), mat(x / 2.0)) # test functions to dispatch From 9b5e1b9cf1473699b54024a8d108ed891411c87d Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sat, 11 Nov 2017 03:10:18 +0100 Subject: [PATCH 3/3] MAINT: add workaround for Numpy bug in MatrixOperator --- odl/operator/tensor_ops.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/odl/operator/tensor_ops.py b/odl/operator/tensor_ops.py index 0abb50aa44e..3fa34103f99 100644 --- a/odl/operator/tensor_ops.py +++ b/odl/operator/tensor_ops.py @@ -858,6 +858,8 @@ def inverse(self): def _call(self, x, out=None): """Return ``self(x[, out])``.""" + from pkg_resources import parse_version + if out is None: return self.range.element(self.matrix.dot(x)) else: @@ -866,8 +868,14 @@ def _call(self, x, out=None): # sparse matrices out[:] = self.matrix.dot(x) else: - with writable_array(out) as out_arr: - self.matrix.dot(x, out=out_arr) + if (parse_version(np.__version__) < parse_version('1.13.0') and + x is out): + # Workaround for bug in Numpy < 1.13 with aliased in and + # out in np.dot + out[:] = self.matrix.dot(x) + else: + with writable_array(out) as out_arr: + self.matrix.dot(x, out=out_arr) def __repr__(self): """Return ``repr(self)``."""