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

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))