From 9d03383a707a7be94b42ac339d52b991cf5d5bc1 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 16 Apr 2019 14:50:51 +0100
Subject: updated tests

---
 Wrappers/Python/test/test_DataContainer.py | 64 +++++++++++++++++++++++-------
 Wrappers/Python/test/test_Operator.py      | 10 ++++-
 2 files changed, 59 insertions(+), 15 deletions(-)

diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 40cd244..8e8ab87 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -174,7 +174,7 @@ class TestDataContainer(unittest.TestCase):
     def binary_add(self):
         print("Test binary add")
         X, Y, Z = 512, 512, 512
-        X, Y, Z = 1024, 512, 512
+        #X, Y, Z = 1024, 512, 512
         steps = [timer()]
         a = numpy.ones((X, Y, Z), dtype='float32')
         steps.append(timer())
@@ -183,9 +183,10 @@ class TestDataContainer(unittest.TestCase):
         #print("a refcount " , sys.getrefcount(a))
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
         ds1 = ds.copy()
+        out = ds.copy()
 
         steps.append(timer())
-        ds.add(ds1, out=ds)
+        ds.add(ds1, out=out)
         steps.append(timer())
         t1 = dt(steps)
         print("ds.add(ds1, out=ds)", dt(steps))
@@ -196,20 +197,25 @@ class TestDataContainer(unittest.TestCase):
         print("ds2 = ds.add(ds1)", dt(steps))
 
         self.assertLess(t1, t2)
-        self.assertEqual(ds.as_array()[0][0][0], 2.)
-
+        self.assertEqual(out.as_array()[0][0][0], 2.)
+        self.assertNumpyArrayEqual(out.as_array(), ds2.as_array())
+        
         ds0 = ds
-        ds0.add(2, out=ds0)
         steps.append(timer())
-        print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0])
-        self.assertEqual(4., ds0.as_array()[0][0][0])
+        ds0.add(2, out=out)
+        steps.append(timer())
+        print("ds0.add(2,out=out)", dt(steps), 3, ds0.as_array()[0][0][0])
+        self.assertEqual(3., out.as_array()[0][0][0])
 
         dt1 = dt(steps)
+        steps.append(timer())
         ds3 = ds0.add(2)
         steps.append(timer())
         print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0])
         dt2 = dt(steps)
         self.assertLess(dt1, dt2)
+        
+        self.assertNumpyArrayEqual(out.as_array(), ds3.as_array())
 
     def binary_subtract(self):
         print("Test binary subtract")
@@ -222,16 +228,17 @@ class TestDataContainer(unittest.TestCase):
         #print("a refcount " , sys.getrefcount(a))
         ds = DataContainer(a, False, ['X', 'Y', 'Z'])
         ds1 = ds.copy()
+        out = ds.copy()
 
         steps.append(timer())
-        ds.subtract(ds1, out=ds)
+        ds.subtract(ds1, out=out)
         steps.append(timer())
         t1 = dt(steps)
         print("ds.subtract(ds1, out=ds)", dt(steps))
-        self.assertEqual(0., ds.as_array()[0][0][0])
+        self.assertEqual(0., out.as_array()[0][0][0])
 
         steps.append(timer())
-        ds2 = ds.subtract(ds1)
+        ds2 = out.subtract(ds1)
         self.assertEqual(-1., ds2.as_array()[0][0][0])
 
         steps.append(timer())
@@ -247,8 +254,8 @@ class TestDataContainer(unittest.TestCase):
         #ds0.__isub__( 2 )
         steps.append(timer())
         print("ds0.subtract(2,out=ds0)", dt(
-            steps), -2., ds0.as_array()[0][0][0])
-        self.assertEqual(-2., ds0.as_array()[0][0][0])
+            steps), -1., ds0.as_array()[0][0][0])
+        self.assertEqual(-1., ds0.as_array()[0][0][0])
 
         dt1 = dt(steps)
         ds3 = ds0.subtract(2)
@@ -256,8 +263,8 @@ class TestDataContainer(unittest.TestCase):
         print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0])
         dt2 = dt(steps)
         self.assertLess(dt1, dt2)
-        self.assertEqual(-2., ds0.as_array()[0][0][0])
-        self.assertEqual(-4., ds3.as_array()[0][0][0])
+        self.assertEqual(-1., ds0.as_array()[0][0][0])
+        self.assertEqual(-3., ds3.as_array()[0][0][0])
 
     def binary_multiply(self):
         print("Test binary multiply")
@@ -300,6 +307,9 @@ class TestDataContainer(unittest.TestCase):
         self.assertLess(dt1, dt2)
         self.assertEqual(4., ds3.as_array()[0][0][0])
         self.assertEqual(2., ds.as_array()[0][0][0])
+        
+        ds.multiply(2.5, out=ds0)
+        self.assertEqual(2.5*2., ds0.as_array()[0][0][0])
 
     def binary_divide(self):
         print("Test binary divide")
@@ -575,6 +585,32 @@ class TestDataContainer(unittest.TestCase):
         norm = dc.norm()
         self.assertEqual(sqnorm, 8.0)
         numpy.testing.assert_almost_equal(norm, numpy.sqrt(8.0), decimal=7)
+        
+    def test_multiply_out(self):
+        print ("test multiply_out")
+        import functools
+        ig = ImageGeometry(10,11,12)
+        u = ig.allocate()
+        a = numpy.ones(u.shape)
+        
+        u.fill(a)
+        
+        numpy.testing.assert_array_equal(a, u.as_array())
+        
+        #u = ig.allocate(ImageGeometry.RANDOM_INT, seed=1)
+        l = functools.reduce(lambda x,y: x*y, (10,11,12), 1)
+        
+        a = numpy.zeros((l, ), dtype=numpy.float32)
+        for i in range(l):
+            a[i] = numpy.sin(2 * i* 3.1415/l)
+        b = numpy.reshape(a, u.shape)
+        u.fill(b)
+        numpy.testing.assert_array_equal(b, u.as_array())
+        
+        u.multiply(2, out=u)
+        c = b * 2
+        numpy.testing.assert_array_equal(u.as_array(), c)
+        
 
 
 
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 0d6cb8b..5c97545 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -68,6 +68,7 @@ class TestOperator(CCPiTestClass):
         numpy.testing.assert_array_equal(y.as_array(), img.as_array())
 
     def test_FiniteDifference(self):
+        print ("test FiniteDifference")
         ##
         N, M = 2, 3
 
@@ -94,10 +95,17 @@ class TestOperator(CCPiTestClass):
         G = Gradient(ig)
 
         u = G.range_geometry().allocate(ImageGeometry.RANDOM_INT)
-        res = G.domain_geometry().allocate(ImageGeometry.RANDOM_INT)
+        res = G.domain_geometry().allocate()
         G.adjoint(u, out=res)
         w = G.adjoint(u)
         self.assertNumpyArrayEqual(res.as_array(), w.as_array())
+        
+        u = G.domain_geometry().allocate(ImageGeometry.RANDOM_INT)
+        res = G.range_geometry().allocate()
+        G.direct(u, out=res)
+        w = G.direct(u)
+        self.assertBlockDataContainerEqual(res, w)
+        
     def test_PowerMethod(self):
         
         N, M = 200, 300
-- 
cgit v1.2.3


From 9a8c03aff07d63f2c635c75fade0da835eb1fb98 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
Date: Tue, 16 Apr 2019 15:28:55 +0100
Subject: fix new implementation of PowerMethod

fixes different values found in #224
---
 .../ccpi/optimisation/operators/LinearOperator.py  | 39 ++++++----------------
 1 file 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
-- 
cgit v1.2.3