moved Power Method in Linear Operator

closes #196

3 files changed, 90 insertions, 294 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index e19304f..8ecefb8 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -6,6 +6,8 @@ Created on Tue Mar  5 15:57:52 2019
 """

 

 from ccpi.optimisation.operators import Operator

+from ccpi.framework import ImageGeometry

+import numpy

 

 

 class LinearOperator(Operator):

@@ -20,3 +22,62 @@ class LinearOperator(Operator):
         

         only available to linear operators'''

         raise NotImplementedError

+    

+    @staticmethod

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

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

+        # Initialise random

+        

+        if x0 is None:

+            #x0 = op.create_image_data()

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

+        

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

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

+        

+        if x0 is None:

+            #x0 = op.create_image_data()

+            x0 = op.allocate_direct()

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

+        

+        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

+

+

diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py
deleted file mode 100755
index fcd0d9e..0000000
--- a/Wrappers/Python/ccpi/optimisation/ops.py
+++ /dev/null
@@ -1,294 +0,0 @@
-# -*- coding: utf-8 -*-
-#   This work is part of the Core Imaging Library developed by
-#   Visual Analytics and Imaging System Group of the Science Technology
-#   Facilities Council, STFC
-
-#   Copyright 2018 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
-
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-
-#       http://www.apache.org/licenses/LICENSE-2.0
-
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-import numpy
-from scipy.sparse.linalg import svds
-from ccpi.framework import DataContainer
-from ccpi.framework import AcquisitionData
-from ccpi.framework import ImageData
-from ccpi.framework import ImageGeometry
-from ccpi.framework import AcquisitionGeometry
-from numbers import Number
-# Maybe operators need to know what types they take as inputs/outputs
-# to not just use generic DataContainer
-
-
-class Operator(object):
-    '''Operator that maps from a space X -> Y'''
-    def __init__(self, **kwargs):
-        self.scalar = 1
-    def is_linear(self):
-        '''Returns if the operator is linear'''
-        return False
-    def direct(self,x, out=None):
-        raise NotImplementedError
-    def size(self):
-        # To be defined for specific class
-        raise NotImplementedError
-    def norm(self):
-        raise NotImplementedError
-    def allocate_direct(self):
-        '''Allocates memory on the Y space'''
-        raise NotImplementedError
-    def allocate_adjoint(self):
-        '''Allocates memory on the X space'''
-        raise NotImplementedError
-    def range_geometry(self):
-        raise NotImplementedError
-    def domain_geometry(self):
-        raise NotImplementedError
-    def __rmul__(self, other):
-        '''reverse multiplication of Operator with number sets the variable scalar in the Operator'''
-        assert isinstance(other, Number)
-        self.scalar = other
-        return self
-
-class LinearOperator(Operator):
-    '''Operator that maps from a space X -> Y'''
-    def is_linear(self):
-        '''Returns if the operator is linear'''
-        return True
-    def adjoint(self,x, out=None):
-        raise NotImplementedError
-        
-class Identity(Operator):
-    def __init__(self):
-        self.s1 = 1.0
-        self.L = 1
-        super(Identity, self).__init__()
-        
-    def direct(self,x,out=None):
-        if out is None:
-            return x.copy()
-        else:
-            out.fill(x)
-    
-    def adjoint(self,x, out=None):
-        if out is None:
-            return x.copy()
-        else:
-            out.fill(x)
-    
-    def size(self):
-        return NotImplemented
-    
-    def get_max_sing_val(self):
-        return self.s1
-
-class TomoIdentity(Operator):
-    def __init__(self, geometry, **kwargs):
-        super(TomoIdentity, self).__init__()
-        self.s1 = 1.0
-        self.geometry = geometry
-        
-    def is_linear(self):
-        return True
-    def direct(self,x,out=None):
-        
-        if out is None:
-            if self.scalar != 1:
-                return x * self.scalar
-            return x.copy()
-        else:
-            if self.scalar != 1:
-                out.fill(x * self.scalar)
-                return
-            out.fill(x)
-            return
-    
-    def adjoint(self,x, out=None):
-        return self.direct(x, out)
-    
-    def norm(self):
-        return self.s1
-    
-    def get_max_sing_val(self):
-        return self.s1
-    def allocate_direct(self):
-        if issubclass(type(self.geometry), ImageGeometry):
-            return ImageData(geometry=self.geometry)
-        elif issubclass(type(self.geometry), AcquisitionGeometry):
-            return AcquisitionData(geometry=self.geometry)
-        else:
-            raise ValueError("Wrong geometry type: expected ImageGeometry of AcquisitionGeometry, got ", type(self.geometry))
-    def allocate_adjoint(self):
-        return self.allocate_direct()
-    def range_geometry(self):
-        return self.geometry
-    def domain_geometry(self):
-        return self.geometry
-    
-    
-
-class FiniteDiff2D(Operator):
-    def __init__(self):
-        self.s1 = 8.0
-        super(FiniteDiff2D, self).__init__()
-        
-    def direct(self,x, out=None):
-        '''Forward differences with Neumann BC.'''
-        # FIXME this seems to be working only with numpy arrays
-        
-        d1 = numpy.zeros_like(x.as_array())
-        d1[:,:-1] = x.as_array()[:,1:] - x.as_array()[:,:-1]
-        d2 = numpy.zeros_like(x.as_array())
-        d2[:-1,:] = x.as_array()[1:,:] - x.as_array()[:-1,:]
-        d = numpy.stack((d1,d2),0)
-        #x.geometry.voxel_num_z = 2
-        return type(x)(d,False,geometry=x.geometry)
-    
-    def adjoint(self,x, out=None):
-        '''Backward differences, Neumann BC.'''
-        Nrows = x.get_dimension_size('horizontal_x')
-        Ncols = x.get_dimension_size('horizontal_y')
-        Nchannels = 1
-        if len(x.shape) == 4:
-            Nchannels = x.get_dimension_size('channel')
-        zer = numpy.zeros((Nrows,1))
-        xxx = x.as_array()[0,:,:-1]
-        #
-        h = numpy.concatenate((zer,xxx), 1) 
-        h -= numpy.concatenate((xxx,zer), 1)
-        
-        zer = numpy.zeros((1,Ncols))
-        xxx = x.as_array()[1,:-1,:]
-        #
-        v  = numpy.concatenate((zer,xxx), 0) 
-        v -= numpy.concatenate((xxx,zer), 0)
-        return type(x)(h + v, False, geometry=x.geometry)
-    
-    def size(self):
-        return NotImplemented
-    
-    def get_max_sing_val(self):
-        return self.s1
-
-def PowerMethodNonsquareOld(op,numiters):
-    # Initialise random
-    # Jakob's
-    #inputsize = op.size()[1]
-    #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))
-    
-    x0 = op.create_image_data()
-    
-    s = numpy.zeros(numiters)
-    # Loop
-    for it in numpy.arange(numiters):
-        x1 = op.adjoint(op.direct(x0))
-        x1norm = numpy.sqrt((x1**2).sum())
-        #print ("x0 **********" ,x0)
-        #print ("x1 **********" ,x1)
-        s[it] = (x1*x0).sum() / (x0*x0).sum()
-        x0 = (1.0/x1norm)*x1
-    return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
-
-#def PowerMethod(op,numiters):
-#    # Initialise random
-#    x0 = np.random.randn(400)
-#    s = np.zeros(numiters)
-#    # Loop
-#    for it in np.arange(numiters):
-#        x1 = np.dot(op.transpose(),np.dot(op,x0))
-#        x1norm = np.sqrt(np.sum(np.dot(x1,x1)))
-#        s[it] = np.dot(x1,x0) / np.dot(x1,x0)
-#        x0 = (1.0/x1norm)*x1
-#    return s, x0
-    
-
-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))
-    
-    if x0 is None:
-        #x0 = op.create_image_data()
-        x0 = op.allocate_direct()
-        x0.fill(numpy.random.randn(*x0.shape))
-    
-    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
-
-class LinearOperatorMatrix(Operator):
-    def __init__(self,A):
-        self.A = A
-        self.s1 = None   # Largest singular value, initially unknown
-        super(LinearOperatorMatrix, self).__init__()
-        
-    def direct(self,x, out=None):
-        if out is None:
-            return type(x)(numpy.dot(self.A,x.as_array()))
-        else:
-            numpy.dot(self.A, x.as_array(), out=out.as_array())
-            
-    
-    def adjoint(self,x, out=None):
-        if out is None:
-            return type(x)(numpy.dot(self.A.transpose(),x.as_array()))
-        else:
-            numpy.dot(self.A.transpose(),x.as_array(), out=out.as_array())
-            
-    
-    def size(self):
-        return self.A.shape
-    
-    def get_max_sing_val(self):
-        # 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)
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 293fb43..70a1cd6 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -8,6 +8,7 @@ import numpy
 from timeit import default_timer as timer
 from ccpi.framework import ImageGeometry
 from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
+from ccpi.optimisation.operators import LinearOperator
 
 def dt(steps):
     return steps[-1] - steps[-2]
@@ -97,7 +98,35 @@ class TestOperator(CCPiTestClass):
         G.adjoint(u, out=res)
         w = G.adjoint(u)
         self.assertNumpyArrayEqual(res.as_array(), w.as_array())
+    def test_PowerMethod(self):
         
+        N, M = 2, 3
+
+        ig = ImageGeometry(N, M)
+        Id = Identity(ig)
+        
+        G = Gradient(ig)
+        
+        uid = Id.domain_geometry().allocate(ImageGeometry.RANDOM_INT)
+        
+        a = LinearOperator.PowerMethod(Id, 10, uid)
+        b = LinearOperator.PowerMethodNonsquare(Id, 10, uid)
+        
+        print ("Edo impl", a[0])
+        print ("old impl", b[0])
+        
+        #self.assertAlmostEqual(a[0], b[0])
+        self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=6)
+        
+        a = LinearOperator.PowerMethod(G, 10, uid)
+        b = LinearOperator.PowerMethodNonsquare(G, 10, uid)
+        
+        print ("Edo impl", a[0])
+        print ("old impl", b[0])
+        self.assertNumpyArrayAlmostEqual(a[0],b[0],decimal=2)
+        #self.assertAlmostEqual(a[0], b[0])
+        
+