Norm refactor (#299)

* Fix upper/lower case inconsistency to allow astra demo to run

* fix import and calls in tests

* add storage in Operator Class of the operator norm

scaled operator is not an Operator

* add test for norm

* reduced iteration in PowerMethod

9 files changed, 36 insertions, 12 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 1d77510..10fe34c 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -104,7 +104,7 @@ class BlockOperator(Operator):
         index = row*self.shape[1]+col
         return self.operators[index]
     
-    def norm(self):
+    def calculate_norm(self):
         norm = [op.norm()**2 for op in self.operators]
         return numpy.sqrt(sum(norm))    
     
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
index aeb6c53..dc8cbf7 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
@@ -49,7 +49,7 @@ class BlockScaledOperator(ScaledOperator):
             return self.scalar * self.operator.adjoint(x, out=out)
         else:
             raise TypeError('Operator is not linear')
-    def norm(self):
+    def calculate_norm(self):
         return numpy.abs(self.scalar) * self.operator.norm()
     def range_geometry(self):
         return self.operator.range_geometry()
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 2c42532..89f3549 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -315,7 +315,7 @@ class FiniteDiff(LinearOperator):
         '''Returns the domain geometry'''
         return self.gm_domain
        
-    def norm(self):
+    def calculate_norm(self):
         x0 = self.gm_domain.allocate('random_int')
         self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0)
         return self.s1
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index e0b8a32..475022e 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -86,7 +86,7 @@ class Gradient(LinearOperator):
     def range_geometry(self):
         return self.gm_range
                                    
-    def norm(self):
+    def calculate_norm(self):
 
         x0 = self.gm_domain.allocate('random')
         self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index a58a296..61428ae 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -35,7 +35,7 @@ class Identity(LinearOperator):
         else:
             out.fill(x)
         
-    def norm(self):
+    def calculate_norm(self):
         return 1.0
         
     def domain_geometry(self):       
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 2d2089b..0cac586 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -8,6 +8,9 @@ from ccpi.optimisation.operators.ScaledOperator import ScaledOperator
 

 class Operator(object):

     '''Operator that maps from a space X -> Y'''

+    def __init__(self, **kwargs):

+        self.__norm = None

+

     def is_linear(self):

         '''Returns if the operator is linear'''

         return False

@@ -16,6 +19,11 @@ class Operator(object):
         raise NotImplementedError

     def norm(self):

         '''Returns the norm of the Operator'''

+        if self.__norm is None:

+            self.__norm = self.calculate_norm()

+        return self.__norm

+    def calculate_norm(self):

+        '''Calculates the norm of the Operator'''

         raise NotImplementedError

     def range_geometry(self):

         '''Returns the range of the Operator: Y space'''

diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index 3fc43b8..da73de6 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -102,7 +102,7 @@ class SymmetrizedGradient(Gradient):
     def range_dim(self):
         return self.gm_range
                                    
-    def norm(self):
+    def calculate_norm(self):
 #        return np.sqrt(4*len(self.domainDim()))        
         #TODO this takes time for big ImageData
         # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12)        
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index a7c5f09..2331856 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -29,7 +29,7 @@ class ZeroOp(Operator):
         else:
             return ImageData(np.zeros(self.gm_domain))
         
-    def norm(self):
+    def calculate_norm(self):
         return 0
     
     def domain_dim(self):       
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 5c97545..2fc45a5 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -51,6 +51,7 @@ class CCPiTestClass(unittest.TestCase):
 
 class TestOperator(CCPiTestClass):
     def test_ScaledOperator(self):
+        print ("test_ScaledOperator")
         ig = ImageGeometry(10,20,30)
         img = ig.allocate()
         scalar = 0.5
@@ -58,7 +59,8 @@ class TestOperator(CCPiTestClass):
         numpy.testing.assert_array_equal(scalar * img.as_array(), sid.direct(img).as_array())
         
 
-    def test_TomoIdentity(self):
+    def test_Identity(self):
+        print ("test_Identity")
         ig = ImageGeometry(10,20,30)
         img = ig.allocate()
         self.assertTrue(img.shape == (30,20,10))
@@ -107,9 +109,10 @@ class TestOperator(CCPiTestClass):
         self.assertBlockDataContainerEqual(res, w)
         
     def test_PowerMethod(self):
+        print ("test_BlockOperator")
         
         N, M = 200, 300
-        niter = 1000
+        niter = 10
         ig = ImageGeometry(N, M)
         Id = Identity(ig)
         
@@ -134,7 +137,20 @@ class TestOperator(CCPiTestClass):
         self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=2)
         #self.assertAlmostEqual(a[0], b[0])
         
+    def test_Norm(self):
+        print ("test_BlockOperator")
+        ##
+        N, M = 200, 300
 
+        ig = ImageGeometry(N, M)
+        G = Gradient(ig)
+        t0 = timer()
+        norm = G.norm()
+        t1 = timer()
+        norm2 = G.norm()
+        t2 = timer()
+        print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1))
+        self.assertLess(t2-t1, t1-t0)
 
 
 
@@ -175,7 +191,7 @@ class TestBlockOperator(unittest.TestCase):
         self.assertTrue(res)
 
     def test_BlockOperator(self):
-        
+        print ("test_BlockOperator")
         
         M, N  = 3, 4
         ig = ImageGeometry(M, N)
@@ -358,7 +374,7 @@ class TestBlockOperator(unittest.TestCase):
             print("Z1", Z1[0][1].as_array())
             print("RES1", RES1[0][1].as_array())
     def test_timedifference(self):
-        
+        print ("test_timedifference")
         M, N ,W = 100, 512, 512
         ig = ImageGeometry(M, N, W)
         arr = ig.allocate('random_int')  
@@ -427,7 +443,7 @@ class TestBlockOperator(unittest.TestCase):
         self.assertGreater(t1,t2)
     
     def test_BlockOperatorLinearValidity(self):
-        
+        print ("test_BlockOperatorLinearValidity")
         
         M, N  = 3, 4
         ig = ImageGeometry(M, N)