From 9d03383a707a7be94b42ac339d52b991cf5d5bc1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca 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 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