fix new implementation of PowerMethod

fixes different values found in #224

1 files changed, 11 insertions, 28 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index bc18f9b..f304cc6 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -24,59 +24,42 @@ class LinearOperator(Operator):
         raise NotImplementedError

     

     @staticmethod

-    def PowerMethod(operator, iterations, x0=None):

+    def PowerMethod(operator, iterations, x_init=None):

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

-        # Initialise random

         

-        if x0 is None:

-            #x0 = op.create_image_data()

+        # Initialise random

+        if x_init is None:

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

-        

+        else:

+            x0 = x_init.copy()

+            

         x1 = operator.domain_geometry().allocate()

         y_tmp = operator.range_geometry().allocate()

         s = numpy.zeros(iterations)

         # Loop

         for it in numpy.arange(iterations):

-            #x1 = operator.adjoint(operator.direct(x0))

             operator.direct(x0,out=y_tmp)

             operator.adjoint(y_tmp,out=x1)

             x1norm = x1.norm()

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

             s[it] = x1.dot(x0) / x0.squared_norm()

-            #x0 = (1.0/x1norm)*x1

-            #x1 *= (1.0 / x1norm)

-            #x0.fill(x1)

             x1.multiply((1.0/x1norm), out=x0)

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

 

     @staticmethod

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

-        # Initialise random

-        # Jakob's

-        # inputsize , outputsize = op.size()

-        #x0 = ImageContainer(numpy.random.randn(*inputsize)

-        # Edo's

-        #vg = ImageGeometry(voxel_num_x=inputsize[0],

-        #                   voxel_num_y=inputsize[1], 

-        #                   voxel_num_z=inputsize[2])

-        #

-        #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x'])

-        #print (x0)

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

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

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

         

-        if x0 is None:

-            #x0 = op.create_image_data()

+        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 = numpy.sqrt((x1*x1).sum())

             x1norm = x1.norm()

-            #print ("x0 **********" ,x0)

-            #print ("x1 **********" ,x1)

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

             x0 = (1.0/x1norm)*x1

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