Merge norm ref (#301)

* 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

* simplified calls to norm for LinearOperators, added default in base class

* add pillow as dependency

* resolve conflict

14 files changed, 33 insertions, 98 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 10fe34c..1577f06 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -104,8 +104,8 @@ class BlockOperator(Operator):
         index = row*self.shape[1]+col
         return self.operators[index]
     
-    def calculate_norm(self):
-        norm = [op.norm()**2 for op in self.operators]
+    def calculate_norm(self, **kwargs):
+        norm = [op.norm(**kwargs)**2 for op in self.operators]
         return numpy.sqrt(sum(norm))    
     
     def direct(self, x, out=None):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
index dc8cbf7..74ba9e4 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
@@ -49,8 +49,8 @@ class BlockScaledOperator(ScaledOperator):
             return self.scalar * self.operator.adjoint(x, out=out)
         else:
             raise TypeError('Operator is not linear')
-    def calculate_norm(self):
-        return numpy.abs(self.scalar) * self.operator.norm()
+    def norm(self, **kwargs):
+        return numpy.abs(self.scalar) * self.operator.norm(**kwargs)
     def range_geometry(self):
         return self.operator.range_geometry()
     def domain_geometry(self):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 89f3549..05278ef 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -314,13 +314,8 @@ class FiniteDiff(LinearOperator):
     def domain_geometry(self):
         '''Returns the domain geometry'''
         return self.gm_domain
-       
-    def calculate_norm(self):
-        x0 = self.gm_domain.allocate('random_int')
-        self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0)
-        return self.s1
-    
-    
+
+
 if __name__ == '__main__':
     
     from ccpi.framework import ImageGeometry
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 475022e..8aa7f5b 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -85,12 +85,7 @@ class Gradient(LinearOperator):
     
     def range_geometry(self):
         return self.gm_range
-                                   
-    def calculate_norm(self):
 
-        x0 = self.gm_domain.allocate('random')
-        self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0)
-        return self.s1
     
     def __rmul__(self, scalar):
         return ScaledOperator(self, scalar) 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index 61428ae..f084504 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 calculate_norm(self):
+    def calculate_norm(self, **kwargs):
         return 1.0
         
     def domain_geometry(self):       
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index f304cc6..55eb692 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -29,7 +29,7 @@ class LinearOperator(Operator):
         

         # Initialise random

         if x_init is None:

-            x0 = operator.domain_geometry().allocate(ImageGeometry.RANDOM_INT)

+            x0 = operator.domain_geometry().allocate(type(operator.domain_geometry()).RANDOM_INT)

         else:

             x0 = x_init.copy()

             

@@ -45,23 +45,11 @@ class LinearOperator(Operator):
             x1.multiply((1.0/x1norm), out=x0)

         return numpy.sqrt(s[-1]), numpy.sqrt(s), x0

 

-    @staticmethod

-    def PowerMethodNonsquare(op,numiters , x_init=None):

-        '''Power method to calculate iteratively the Lipschitz constant'''

-        

-        if x_init is None:

-            x0 = op.allocate_direct()

-            x0.fill(numpy.random.randn(*x0.shape))

-        else:

-            x0 = x_init.copy()

-        

-        s = numpy.zeros(numiters)

-        # Loop

-        for it in numpy.arange(numiters):

-            x1 = op.adjoint(op.direct(x0))

-            x1norm = x1.norm()

-            s[it] = (x1*x0).sum() / (x0.squared_norm())

-            x0 = (1.0/x1norm)*x1

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

+        iterations = kwargs.get('iterations', 25)

+        s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0)

+        return s1

 

 

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index 62e22e0..8cc4c71 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -30,22 +30,7 @@ class LinearOperatorMatrix(LinearOperator):
     def size(self):
         return self.A.shape
     
-    def get_max_sing_val(self):
+    def calculate_norm(self, **kwargs):
         # If unknown, compute and store. If known, simply return it.
-        if self.s1 is None:
-            self.s1 = svds(self.A,1,return_singular_vectors=False)[0]
-            return self.s1
-        else:
-            return self.s1
-    def allocate_direct(self):
-        '''allocates the memory to hold the result of adjoint'''
-        #numpy.dot(self.A.transpose(),x.as_array())
-        M_A, N_A = self.A.shape
-        out = numpy.zeros((N_A,1))
-        return DataContainer(out)
-    def allocate_adjoint(self):
-        '''allocate the memory to hold the result of direct'''
-        #numpy.dot(self.A.transpose(),x.as_array())
-        M_A, N_A = self.A.shape
-        out = numpy.zeros((M_A,1))
-        return DataContainer(out)
+        return svds(self.A,1,return_singular_vectors=False)[0]
+        
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 0cac586..c1adf62 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -17,12 +17,12 @@ class Operator(object):
     def direct(self,x, out=None):

         '''Returns the application of the Operator on x'''

         raise NotImplementedError

-    def norm(self):

+    def norm(self, **kwargs):

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

         if self.__norm is None:

-            self.__norm = self.calculate_norm()

+            self.__norm = self.calculate_norm(**kwargs)

         return self.__norm

-    def calculate_norm(self):

+    def calculate_norm(self, **kwargs):

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

         raise NotImplementedError

     def range_geometry(self):

diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
index ba0049e..f7be5de 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
@@ -40,8 +40,8 @@ class ScaledOperator(object):
                 out *= self.scalar
         else:
             raise TypeError('Operator is not linear')
-    def norm(self):
-        return numpy.abs(self.scalar) * self.operator.norm()
+    def norm(self, **kwargs):
+        return numpy.abs(self.scalar) * self.operator.norm(**kwargs)
     def range_geometry(self):
         return self.operator.range_geometry()
     def domain_geometry(self):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 5e318ff..eb51cc3 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -10,7 +10,7 @@ import scipy.sparse as sp
 import numpy as np
 from ccpi.framework import ImageData
 
-class SparseFiniteDiff():
+class SparseFiniteDiff(object):
     
     def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
         
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index da73de6..1897040 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -75,42 +75,12 @@ class SymmetrizedGradient(Gradient):
                 tmp += FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).direct(x.get_item(j))
             res.append(tmp)   
         return res            
-                
-                
-        
-#        for 
-        
-#        tmp = numpy.zeros(self.gm_domain)
-#        
-#        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) +  \
-#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2])
-#                 
-#        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) +  \
-#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1])                 
-#
-#        return type(x)(tmp)          
-            
-    def alloc_domain_dim(self):
-        return ImageData(numpy.zeros(self.gm_domain))
-    
-    def alloc_range_dim(self):
-        return ImageData(numpy.zeros(self.range_dim))
-    
+
     def domain_dim(self):
         return self.gm_domain
     
     def range_dim(self):
         return self.gm_range
-                                   
-    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)        
-        x0 = ImageData(np.random.random_sample(self.domain_dim()))
-        self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0)
-        return self.s1 
-    
-
 
 if __name__ == '__main__':   
     
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index 2331856..a2a94b7 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -29,8 +29,8 @@ class ZeroOp(Operator):
         else:
             return ImageData(np.zeros(self.gm_domain))
         
-    def calculate_norm(self):
-        return 0
+    def calculate_norm(self, **kwargs):
+        return 0.
     
     def domain_dim(self):       
         return self.gm_domain  
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml
index dd3238e..163d326 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/conda-recipe/meta.yaml
@@ -36,6 +36,7 @@ requirements:
     - scipy
     #- matplotlib
     - h5py
+    - pillow
   
 about:
   home: http://www.ccpi.ac.uk
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 2fc45a5..d890e46 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -121,19 +121,20 @@ class TestOperator(CCPiTestClass):
         uid = Id.domain_geometry().allocate(ImageGeometry.RANDOM_INT, seed=1)
         
         a = LinearOperator.PowerMethod(Id, niter, uid)
-        b = LinearOperator.PowerMethodNonsquare(Id, niter, uid)
-        
+        #b = LinearOperator.PowerMethodNonsquare(Id, niter, uid)
+        b = LinearOperator.PowerMethod(Id, niter)
         print ("Edo impl", a[0])
-        print ("old impl", b[0])
+        print ("None impl", b[0])
         
         #self.assertAlmostEqual(a[0], b[0])
         self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=6)
         
         a = LinearOperator.PowerMethod(G, niter, uid)
-        b = LinearOperator.PowerMethodNonsquare(G, niter, uid)
+        b = LinearOperator.PowerMethod(G, niter)
+        #b = LinearOperator.PowerMethodNonsquare(G, niter, uid)
         
         print ("Edo impl", a[0])
-        print ("old impl", b[0])
+        #print ("old impl", b[0])
         self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=2)
         #self.assertAlmostEqual(a[0], b[0])