From c2cb31c0432dadf9f0257c6ad4a87a669a6fa6b2 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 16 Jan 2020 14:57:15 +0000 Subject: add force recalculation of norm (#487) * add force recalculation of norm closes #394 better handling and documentation of Operator norm method. * added docstring and test for dot_test --- .../ccpi/optimisation/operators/LinearOperator.py | 25 ++++++++++++--- .../Python/ccpi/optimisation/operators/Operator.py | 21 +++++++++++-- Wrappers/Python/test/test_Operator.py | 36 ++++++++++++++++++++-- 3 files changed, 74 insertions(+), 8 deletions(-) diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py index 5ff5b01..60a4c2a 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py @@ -72,14 +72,22 @@ class LinearOperator(Operator): return numpy.sqrt(s[-1]), numpy.sqrt(s), x0 def calculate_norm(self, **kwargs): - '''Returns the norm of the LinearOperator as calculated by the PowerMethod''' - x0 = kwargs.get('x0', None) + '''Returns the norm of the LinearOperator as calculated by the PowerMethod + + :param iterations: number of iterations to run + :type iterations: int, optional, default = 25 + :param x_init: starting point for the iteration in the operator domain + :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, None + :parameter force: forces the recalculation of the norm + :type force: boolean, default :code:`False` + ''' + x0 = kwargs.get('x_init', None) iterations = kwargs.get('iterations', 25) s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0) return s1 @staticmethod - def dot_test(operator, domain_init=None, range_init=None, verbose=False): + def dot_test(operator, domain_init=None, range_init=None, verbose=False, **kwargs): r'''Does a dot linearity test on the operator Evaluates if the following equivalence holds @@ -88,10 +96,18 @@ class LinearOperator(Operator): Ax\times y = y \times A^Tx + The equivalence is tested within a user specified precision + + .. code:: + + abs(desired-actual) < 1.5 * 10**(-decimal) + :param operator: operator to test :param range_init: optional initialisation container in the operator range :param domain_init: optional initialisation container in the operator domain :returns: boolean, True if the test is passed. + :param decimal: desired precision + :type decimal: int, optional, default 4 ''' if range_init is None: y = operator.range_geometry().allocate('random_int') @@ -108,8 +124,9 @@ class LinearOperator(Operator): b = by.dot(x) if verbose: print ('Left hand side {}, \nRight hand side {}'.format(a, b)) + print ("decimal ", kwargs.get('decimal', 4)) try: - numpy.testing.assert_almost_equal(abs((a-b)/a), 0, decimal=4) + numpy.testing.assert_almost_equal(abs((a-b)/a), 0., decimal=kwargs.get('decimal',4)) return True except AssertionError as ae: print (ae) diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py index 3cd0899..14d53e8 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py @@ -33,8 +33,25 @@ class Operator(object): '''Returns the application of the Operator on x''' raise NotImplementedError def norm(self, **kwargs): - '''Returns the norm of the Operator''' - if self.__norm is None: + '''Returns the norm of the Operator + + Calling norm triggers the calculation of the norm of the operator. Normally this + is a computationally expensive task, therefore we store the result of norm into + a member of the class. If the calculation has already run, following calls to + norm just return the saved member. + It is possible to force recalculation by setting the optional force parameter. Notice that + norm doesn't take notice of how many iterations or of the initialisation of the PowerMethod, + so in case you want to recalculate by setting a higher number of iterations or changing the + starting point or both you need to set :code:`force=True` + + :param iterations: number of iterations to run + :type iterations: int, optional, default = 25 + :param x_init: starting point for the iteration in the operator domain + :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, default None + :parameter force: forces the recalculation of the norm + :type force: boolean, default :code:`False` + ''' + if self.__norm is None or kwargs.get('force', False): self.__norm = self.calculate_norm(**kwargs) return self.__norm def calculate_norm(self, **kwargs): diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py index 5c659a4..57dd41f 100644 --- a/Wrappers/Python/test/test_Operator.py +++ b/Wrappers/Python/test/test_Operator.py @@ -191,6 +191,7 @@ class TestOperator(CCPiTestClass): def test_Norm(self): print ("test_BlockOperator") ## + numpy.random.seed(1) N, M = 200, 300 ig = ImageGeometry(N, M) @@ -200,8 +201,26 @@ class TestOperator(CCPiTestClass): t1 = timer() norm2 = G.norm() t2 = timer() - print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1)) + norm3 = G.norm(force=True) + t3 = timer() + print ("Norm dT1 {} dT2 {} dT3 {}".format(t1-t0,t2-t1, t3-t2)) self.assertLess(t2-t1, t1-t0) + self.assertLess(t2-t1, t3-t2) + + numpy.random.seed(1) + t4 = timer() + norm4 = G.norm(iterations=50, force=True) + t5 = timer() + self.assertLess(t2-t1, t5-t4) + + numpy.random.seed(1) + t4 = timer() + norm5 = G.norm(x_init=ig.allocate('random'), iterations=50, force=True) + t5 = timer() + self.assertLess(t2-t1, t5-t4) + for n in [norm, norm2, norm3, norm4, norm5]: + print ("norm {}", format(n)) + class TestGradients(CCPiTestClass): def setUp(self): @@ -301,9 +320,22 @@ class TestGradients(CCPiTestClass): lhs3 = E3.direct(u3).dot(w3) rhs3 = u3.dot(E3.adjoint(w3)) numpy.testing.assert_almost_equal(lhs3, rhs3) + self.assertAlmostEqual(lhs3, rhs3) + print (lhs3, rhs3, abs((rhs3-lhs3)/rhs3) , 1.5 * 10**(-4), abs((rhs3-lhs3)/rhs3) < 1.5 * 10**(-4)) + self.assertTrue( LinearOperator.dot_test(E3, range_init = w3, domain_init=u3) ) + def test_dot_test(self): + Grad3 = Gradient(self.ig3, correlation = 'Space', backend='numpy') + # self.assertAlmostEqual(lhs3, rhs3) - # self.assertTrue( LinearOperator.dot_test(E3) ) + self.assertTrue( LinearOperator.dot_test(Grad3 , verbose=True)) + self.assertTrue( LinearOperator.dot_test(Grad3 , decimal=6, verbose=True)) + def test_dot_test2(self): + Grad3 = Gradient(self.ig3, correlation = 'SpaceChannel', backend='c') + + # self.assertAlmostEqual(lhs3, rhs3) + self.assertTrue( LinearOperator.dot_test(Grad3 , verbose=True)) + self.assertTrue( LinearOperator.dot_test(Grad3 , decimal=6, verbose=True)) -- cgit v1.2.3