From 3f5c27f14cf921246ee55582f42bf5322f3cd590 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 7 Nov 2018 14:11:44 +0000 Subject: Test norm (#162) * enable win build with unittest. skip cvx * unittest on funcs * fixed Norm1 and function * test chaining * a few fixes for unittest * investigating Norm2 * restored Norm2 TODO: rewrite it and test * removed comment --- Wrappers/Python/ccpi/framework.py | 96 +-- Wrappers/Python/ccpi/optimisation/funcs.py | 30 +- Wrappers/Python/ccpi/optimisation/ops.py | 6 +- Wrappers/Python/conda-recipe/meta.yaml | 5 +- Wrappers/Python/conda-recipe/run_test.py | 1012 +++++++++++++++------------- Wrappers/Python/wip/demo_compare_cvx.py | 62 +- 6 files changed, 673 insertions(+), 538 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py index 7249d64..4d74d2b 100644 --- a/Wrappers/Python/ccpi/framework.py +++ b/Wrappers/Python/ccpi/framework.py @@ -627,8 +627,8 @@ class DataContainer(object): .format(len(self.shape),len(new_order))) - def copy(self): - '''alias of clone''' + def copy(self): + '''alias of clone''' return self.clone() ## binary operations @@ -647,49 +647,51 @@ class DataContainer(object): elif issubclass(type(out), DataContainer) and issubclass(type(x2), DataContainer): if self.check_dimensions(out) and self.check_dimensions(x2): - pwop(self.as_array(), x2.as_array(), *args, out=out.as_array(), **kwargs ) - return type(self)(out.as_array(), - deep_copy=False, - dimension_labels=self.dimension_labels, - geometry=self.geometry) + pwop(self.as_array(), x2.as_array(), out=out.as_array(), *args, **kwargs ) + #return type(self)(out.as_array(), + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) + return out else: raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), DataContainer) and isinstance(x2, (int,float,complex)): if self.check_dimensions(out): - pwop(self.as_array(), x2, *args, out=out.as_array(), **kwargs ) - return type(self)(out.as_array(), - deep_copy=False, - dimension_labels=self.dimension_labels, - geometry=self.geometry) + pwop(self.as_array(), x2, out=out.as_array(), *args, **kwargs ) + #return type(self)(out.as_array(), + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) + return out else: raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), numpy.ndarray): if self.array.shape == out.shape and self.array.dtype == out.dtype: - pwop(self.as_array(), x2 , *args, out=out, **kwargs) - return type(self)(out, - deep_copy=False, - dimension_labels=self.dimension_labels, - geometry=self.geometry) + pwop(self.as_array(), x2 , out=out, *args, **kwargs) + #return type(self)(out, + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) else: raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) def add(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.add, other, *args, out=out, **kwargs) + return self.pixel_wise_binary(numpy.add, other, out=out, *args, **kwargs) def subtract(self, other, out=None , *args, **kwargs): - return self.pixel_wise_binary(numpy.subtract, other, *args, out=out, **kwargs) + return self.pixel_wise_binary(numpy.subtract, other, out=out, *args, **kwargs) def multiply(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.multiply, other, *args, out=out, **kwargs) + return self.pixel_wise_binary(numpy.multiply, other, out=out, *args, **kwargs) def divide(self, other , out=None ,*args, **kwargs): - return self.pixel_wise_binary(numpy.divide, other, *args, out=out, **kwargs) + return self.pixel_wise_binary(numpy.divide, other, out=out, *args, **kwargs) def power(self, other , out=None, *args, **kwargs): - return self.pixel_wise_binary(numpy.power, other, *args, out=out, **kwargs) + return self.pixel_wise_binary(numpy.power, other, out=out, *args, **kwargs) - def maximum(self,x2, out=None): - return self.pixel_wise_binary(numpy.maximum, x2=x2, out=out) + def maximum(self,x2, out=None, *args, **kwargs): + return self.pixel_wise_binary(numpy.maximum, x2=x2, out=out, *args, **kwargs) ## unary operations @@ -702,36 +704,36 @@ class DataContainer(object): geometry=self.geometry) elif issubclass(type(out), DataContainer): if self.check_dimensions(out): - pwop(self.as_array(), *args, out=out.as_array(), **kwargs ) - return type(self)(out.as_array(), - deep_copy=False, - dimension_labels=self.dimension_labels, - geometry=self.geometry) + pwop(self.as_array(), out=out.as_array(), *args, **kwargs ) + #return type(self)(out.as_array(), + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) else: raise ValueError(message(type(self),"Wrong size for data memory: ", out.shape,self.shape)) elif issubclass(type(out), numpy.ndarray): if self.array.shape == out.shape and self.array.dtype == out.dtype: - pwop(self.as_array(), *args, out=out, **kwargs) - return type(self)(out, - deep_copy=False, - dimension_labels=self.dimension_labels, - geometry=self.geometry) + pwop(self.as_array(), out=out, *args, **kwargs) + #return type(self)(out, + # deep_copy=False, + # dimension_labels=self.dimension_labels, + # geometry=self.geometry) else: raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) - def abs(self, out=None): - return self.pixel_wise_unary(numpy.abs, out=out) + def abs(self, out=None, *args, **kwargs): + return self.pixel_wise_unary(numpy.abs, out=out, *args, **kwargs) - def sign(self, out=None): + def sign(self, out=None, *args, **kwargs): #out = numpy.sign(self.as_array() ) #return type(self)(out, # deep_copy=False, # dimension_labels=self.dimension_labels, # geometry=self.geometry) - return self.pixel_wise_unary(numpy.sign , out=out) + return self.pixel_wise_unary(numpy.sign , out=out, *args, **kwargs) - def sqrt(self, out=None): - return self.pixel_wise_unary(numpy.sqrt, out=out) + def sqrt(self, out=None, *args, **kwargs): + return self.pixel_wise_unary(numpy.sqrt, out=out, *args, **kwargs) #def __abs__(self): # operation = FM.OPERATION.ABS @@ -739,9 +741,9 @@ class DataContainer(object): # __abs__ ## reductions - def sum(self): - return self.as_array().sum() - + def sum(self, out=None, *args, **kwargs): + return self.as_array().sum(*args, **kwargs) + #return self.pixel_wise_unary(numpy.sum, out=out, *args, **kwargs) class ImageData(DataContainer): '''DataContainer for holding 2D or 3D DataContainer''' @@ -805,7 +807,7 @@ class ImageData(DataContainer): raise ValueError('Please pass either a DataContainer, ' +\ 'a numpy array or a geometry') else: - if type(array) == DataContainer: + if issubclass(type(array) , DataContainer): # if the array is a DataContainer get the info from there if not ( array.number_of_dimensions == 2 or \ array.number_of_dimensions == 3 or \ @@ -817,7 +819,7 @@ class ImageData(DataContainer): # array.dimension_labels, **kwargs) super(ImageData, self).__init__(array.as_array(), deep_copy, array.dimension_labels, **kwargs) - elif type(array) == numpy.ndarray: + elif issubclass(type(array) , numpy.ndarray): if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ): raise ValueError( 'Number of dimensions are not 2 or 3 or 4 : {0}'\ @@ -914,7 +916,7 @@ class AcquisitionData(DataContainer): dimension_labels, **kwargs) else: - if type(array) == DataContainer: + if issubclass(type(array) ,DataContainer): # if the array is a DataContainer get the info from there if not ( array.number_of_dimensions == 2 or \ array.number_of_dimensions == 3 or \ @@ -926,7 +928,7 @@ class AcquisitionData(DataContainer): # array.dimension_labels, **kwargs) super(AcquisitionData, self).__init__(array.as_array(), deep_copy, array.dimension_labels, **kwargs) - elif type(array) == numpy.ndarray: + elif issubclass(type(array) ,numpy.ndarray): if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ): raise ValueError( 'Number of dimensions are not 2 or 3 or 4 : {0}'\ diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index 0ed77ae..db00e9f 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -45,6 +45,7 @@ class Function(object): def gradient(self, x, out=None): raise NotImplementedError def proximal(self, x, tau, out=None): raise NotImplementedError + class Norm2(Function): def __init__(self, @@ -108,7 +109,6 @@ class Norm2(Function): else: raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) ) - class TV2D(Norm2): @@ -193,20 +193,21 @@ class ZeroFun(Function): def prox(self,x,tau): return x.copy() - def proximal(self, x, tau, out=None): if out is None: - return self.prox(x,tau) + return self.prox(x, tau) else: if isSizeCorrect(out, x): - # check dimensionality - if issubclass(type(out), DataContainer): - out.fill(x) - - elif issubclass(type(out) , numpy.ndarray): - out[:] = x - else: - raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) ) + # check dimensionality + if issubclass(type(out), DataContainer): + out.fill(x) + + elif issubclass(type(out) , numpy.ndarray): + out[:] = x + else: + raise ValueError ('Wrong size: x{0} out{1}' + .format(x.shape,out.shape) ) + # A more interesting example, least squares plus 1-norm minimization. # Define class to represent 1-norm including prox function class Norm1(Function): @@ -229,11 +230,12 @@ class Norm1(Function): def prox(self,x,tau): return (x.abs() - tau*self.gamma).maximum(0) * x.sign() + def proximal(self, x, tau, out=None): if out is None: - return self.prox(x,tau) + return self.prox(x, tau) else: - if isSizeCorrect(out, x): + if isSizeCorrect(x,out): # check dimensionality if issubclass(type(out), DataContainer): v = (x.abs() - tau*self.gamma).maximum(0) @@ -248,7 +250,6 @@ class Norm1(Function): else: raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) ) - # Box constraints indicator function. Calling returns 0 if argument is within # the box. The prox operator is projection onto the box. Only implements one # scalar lower and one upper as constraint on all elements. Should generalise @@ -290,4 +291,3 @@ class IndicatorBox(Function): x.sign(out=self.sign_x) out.__imul__( self.sign_x ) - diff --git a/Wrappers/Python/ccpi/optimisation/ops.py b/Wrappers/Python/ccpi/optimisation/ops.py index c6775a8..96f85d8 100755 --- a/Wrappers/Python/ccpi/optimisation/ops.py +++ b/Wrappers/Python/ccpi/optimisation/ops.py @@ -47,6 +47,7 @@ class Operator(object): class Identity(Operator): def __init__(self): self.s1 = 1.0 + self.L = 1 super(Identity, self).__init__() def direct(self,x,out=None): @@ -110,18 +111,19 @@ class FiniteDiff2D(Operator): 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_x') + Ncols = x.get_dimension_size('horizontal_y') Nchannels = 1 if len(x.shape) == 4: Nchannels = x.get_dimension_size('channel') diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/conda-recipe/meta.yaml index 4b19a68..327cdcb 100644 --- a/Wrappers/Python/conda-recipe/meta.yaml +++ b/Wrappers/Python/conda-recipe/meta.yaml @@ -1,6 +1,6 @@ package: name: ccpi-framework - version: 0.11.0 + version: 0.11.1 build: @@ -12,7 +12,6 @@ build: requirements: build: - python - - numpy - setuptools run: @@ -23,7 +22,7 @@ requirements: run: - python - numpy - - cvxpy + - cvxpy # [not win] - scipy - matplotlib diff --git a/Wrappers/Python/conda-recipe/run_test.py b/Wrappers/Python/conda-recipe/run_test.py index c1dc59b..d6bc340 100755 --- a/Wrappers/Python/conda-recipe/run_test.py +++ b/Wrappers/Python/conda-recipe/run_test.py @@ -14,12 +14,17 @@ from ccpi.optimisation.funcs import Norm2sq from ccpi.optimisation.funcs import ZeroFun from ccpi.optimisation.funcs import Norm1 from ccpi.optimisation.funcs import TV2D +from ccpi.optimisation.funcs import Norm2 from ccpi.optimisation.ops import LinearOperatorMatrix from ccpi.optimisation.ops import TomoIdentity from ccpi.optimisation.ops import Identity -from cvxpy import * +try: + from cvxpy import * + cvx_not_installable = False +except ImportError: + cvx_not_installable = True def aid(x): @@ -27,76 +32,99 @@ def aid(x): # block address of an array. return x.__array_interface__['data'][0] -def dt (steps): + +def dt(steps): return steps[-1] - steps[-2] + class TestDataContainer(unittest.TestCase): - + def create_simple_ImageData(self): + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + return (ig, Phantom) + + def create_DataContainer(self, X,Y,Z, value=1): + steps = [timer()] + a = value * numpy.ones((X, Y, Z), dtype='float32') + steps.append(timer()) + t0 = dt(steps) + #print("a refcount " , sys.getrefcount(a)) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + return ds + def test_creation_nocopy(self): - shape = (2,3,4,5) + shape = (2, 3, 4, 5) size = shape[0] for i in range(1, len(shape)): size = size * shape[i] #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range( size )]) + a = numpy.asarray([i for i in range(size)]) #print("a refcount " , sys.getrefcount(a)) a = numpy.reshape(a, shape) #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z' ,'W']) + ds = DataContainer(a, False, ['X', 'Y', 'Z', 'W']) #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a),3) - self.assertEqual(ds.dimension_labels , {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) - + self.assertEqual(sys.getrefcount(a), 3) + self.assertEqual(ds.dimension_labels, {0: 'X', 1: 'Y', 2: 'Z', 3: 'W'}) + def testGb_creation_nocopy(self): - X,Y,Z = 512,512,512 - X,Y,Z = 256,512,512 + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) print("test clone") #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a),3) + self.assertEqual(sys.getrefcount(a), 3) ds1 = ds.copy() self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) ds1 = ds.clone() self.assertNotEqual(aid(ds.as_array()), aid(ds1.as_array())) - - + def testInlineAlgebra(self): - print ("Test Inline Algebra") - X,Y,Z = 1024,512,512 - X,Y,Z = 256,512,512 + print("Test Inline Algebra") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) print(t0) #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) #ds.__iadd__( 2 ) ds += 2 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],3.) - #ds.__isub__( 2 ) + self.assertEqual(ds.as_array()[0][0][0], 3.) + #ds.__isub__( 2 ) ds -= 2 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],1.) + self.assertEqual(ds.as_array()[0][0][0], 1.) #ds.__imul__( 2 ) ds *= 2 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],2.) + self.assertEqual(ds.as_array()[0][0][0], 2.) #ds.__idiv__( 2 ) ds /= 2 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],1.) - + self.assertEqual(ds.as_array()[0][0][0], 1.) + ds1 = ds.copy() #ds1.__iadd__( 1 ) ds1 += 1 @@ -104,367 +132,371 @@ class TestDataContainer(unittest.TestCase): ds += ds1 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],3.) + self.assertEqual(ds.as_array()[0][0][0], 3.) #ds.__isub__( ds1 ) ds -= ds1 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],1.) + self.assertEqual(ds.as_array()[0][0][0], 1.) #ds.__imul__( ds1 ) ds *= ds1 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],2.) + self.assertEqual(ds.as_array()[0][0][0], 2.) #ds.__idiv__( ds1 ) ds /= ds1 steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],1.) - - + self.assertEqual(ds.as_array()[0][0][0], 1.) + def test_unary_operations(self): - print ("Test unary operations") - X,Y,Z = 1024,512,512 - X,Y,Z = 256,512,512 + print("Test unary operations") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = -numpy.ones((X,Y,Z), dtype='float32') + a = -numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) print(t0) #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) - + ds = DataContainer(a, False, ['X', 'Y', 'Z']) + ds.sign(out=ds) steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],-1.) - + self.assertEqual(ds.as_array()[0][0][0], -1.) + ds.abs(out=ds) steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],1.) - - ds.__imul__( 2 ) + self.assertEqual(ds.as_array()[0][0][0], 1.) + + ds.__imul__(2) ds.sqrt(out=ds) steps.append(timer()) print(dt(steps)) - self.assertEqual(ds.as_array()[0][0][0],numpy.sqrt(2., dtype='float32')) - - - + self.assertEqual(ds.as_array()[0][0][0], + numpy.sqrt(2., dtype='float32')) + def test_binary_operations(self): self.binary_add() self.binary_subtract() self.binary_multiply() self.binary_divide() - + def binary_add(self): - print ("Test binary add") - X,Y,Z = 512,512,512 - X,Y,Z = 256,512,512 + print("Test binary add") + X, Y, Z = 512, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) - + #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) ds1 = ds.copy() - + steps.append(timer()) ds.add(ds1, out=ds) steps.append(timer()) t1 = dt(steps) - print("ds.add(ds1, out=ds)",dt(steps)) + print("ds.add(ds1, out=ds)", dt(steps)) steps.append(timer()) ds2 = ds.add(ds1) steps.append(timer()) t2 = dt(steps) - print("ds2 = ds.add(ds1)",dt(steps)) - - self.assertLess(t1,t2) - self.assertEqual(ds.as_array()[0][0][0] , 2.) - + print("ds2 = ds.add(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 2.) + ds0 = ds - ds0.add(2,out=ds0) + ds0.add(2, out=ds0) steps.append(timer()) - print ("ds0.add(2,out=ds0)", dt(steps), 3 , ds0.as_array()[0][0][0]) + print("ds0.add(2,out=ds0)", dt(steps), 3, ds0.as_array()[0][0][0]) self.assertEqual(4., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) + + dt1 = dt(steps) ds3 = ds0.add(2) steps.append(timer()) - print ("ds3 = ds0.add(2)", dt(steps), 5 , ds3.as_array()[0][0][0]) + print("ds3 = ds0.add(2)", dt(steps), 5, ds3.as_array()[0][0][0]) dt2 = dt(steps) - self.assertLess(dt1,dt2) - + self.assertLess(dt1, dt2) + def binary_subtract(self): - print ("Test binary subtract") - X,Y,Z = 512,512,512 + print("Test binary subtract") + X, Y, Z = 512, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) - + #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) ds1 = ds.copy() - + steps.append(timer()) ds.subtract(ds1, out=ds) steps.append(timer()) t1 = dt(steps) - print("ds.subtract(ds1, out=ds)",dt(steps)) + print("ds.subtract(ds1, out=ds)", dt(steps)) self.assertEqual(0., ds.as_array()[0][0][0]) - + steps.append(timer()) ds2 = ds.subtract(ds1) self.assertEqual(-1., ds2.as_array()[0][0][0]) - + steps.append(timer()) t2 = dt(steps) - print("ds2 = ds.subtract(ds1)",dt(steps)) - - self.assertLess(t1,t2) - + print("ds2 = ds.subtract(ds1)", dt(steps)) + + self.assertLess(t1, t2) + del ds1 ds0 = ds.copy() steps.append(timer()) - ds0.subtract(2,out=ds0) + ds0.subtract(2, out=ds0) #ds0.__isub__( 2 ) steps.append(timer()) - print ("ds0.subtract(2,out=ds0)", dt(steps), -2. , ds0.as_array()[0][0][0]) + print("ds0.subtract(2,out=ds0)", dt( + steps), -2., ds0.as_array()[0][0][0]) self.assertEqual(-2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) + + dt1 = dt(steps) ds3 = ds0.subtract(2) steps.append(timer()) - print ("ds3 = ds0.subtract(2)", dt(steps), 0. , ds3.as_array()[0][0][0]) + print("ds3 = ds0.subtract(2)", dt(steps), 0., ds3.as_array()[0][0][0]) dt2 = dt(steps) - self.assertLess(dt1,dt2) + self.assertLess(dt1, dt2) self.assertEqual(-2., ds0.as_array()[0][0][0]) self.assertEqual(-4., ds3.as_array()[0][0][0]) - + def binary_multiply(self): - print ("Test binary multiply") - X,Y,Z = 1024,512,512 - X,Y,Z = 256,512,512 + print("Test binary multiply") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) - + #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) ds1 = ds.copy() - + steps.append(timer()) ds.multiply(ds1, out=ds) steps.append(timer()) t1 = dt(steps) - print("ds.multiply(ds1, out=ds)",dt(steps)) + print("ds.multiply(ds1, out=ds)", dt(steps)) steps.append(timer()) ds2 = ds.multiply(ds1) steps.append(timer()) t2 = dt(steps) - print("ds2 = ds.multiply(ds1)",dt(steps)) - - self.assertLess(t1,t2) - + print("ds2 = ds.multiply(ds1)", dt(steps)) + + self.assertLess(t1, t2) + ds0 = ds - ds0.multiply(2,out=ds0) + ds0.multiply(2, out=ds0) steps.append(timer()) - print ("ds0.multiply(2,out=ds0)", dt(steps), 2. , ds0.as_array()[0][0][0]) + print("ds0.multiply(2,out=ds0)", dt( + steps), 2., ds0.as_array()[0][0][0]) self.assertEqual(2., ds0.as_array()[0][0][0]) - - dt1 = dt(steps) + + dt1 = dt(steps) ds3 = ds0.multiply(2) steps.append(timer()) - print ("ds3 = ds0.multiply(2)", dt(steps), 4. , ds3.as_array()[0][0][0]) + print("ds3 = ds0.multiply(2)", dt(steps), 4., ds3.as_array()[0][0][0]) dt2 = dt(steps) - self.assertLess(dt1,dt2) + self.assertLess(dt1, dt2) self.assertEqual(4., ds3.as_array()[0][0][0]) self.assertEqual(2., ds.as_array()[0][0][0]) - + def binary_divide(self): - print ("Test binary divide") - X,Y,Z = 1024,512,512 - X,Y,Z = 256,512,512 + print("Test binary divide") + X, Y, Z = 1024, 512, 512 + X, Y, Z = 256, 512, 512 steps = [timer()] - a = numpy.ones((X,Y,Z), dtype='float32') + a = numpy.ones((X, Y, Z), dtype='float32') steps.append(timer()) t0 = dt(steps) - + #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, False, ['X', 'Y','Z']) + ds = DataContainer(a, False, ['X', 'Y', 'Z']) ds1 = ds.copy() - + steps.append(timer()) ds.divide(ds1, out=ds) steps.append(timer()) t1 = dt(steps) - print("ds.divide(ds1, out=ds)",dt(steps)) + print("ds.divide(ds1, out=ds)", dt(steps)) steps.append(timer()) ds2 = ds.divide(ds1) steps.append(timer()) t2 = dt(steps) - print("ds2 = ds.divide(ds1)",dt(steps)) - - self.assertLess(t1,t2) - self.assertEqual(ds.as_array()[0][0][0] , 1.) - + print("ds2 = ds.divide(ds1)", dt(steps)) + + self.assertLess(t1, t2) + self.assertEqual(ds.as_array()[0][0][0], 1.) + ds0 = ds - ds0.divide(2,out=ds0) + ds0.divide(2, out=ds0) steps.append(timer()) - print ("ds0.divide(2,out=ds0)", dt(steps), 0.5 , ds0.as_array()[0][0][0]) + print("ds0.divide(2,out=ds0)", dt(steps), 0.5, ds0.as_array()[0][0][0]) self.assertEqual(0.5, ds0.as_array()[0][0][0]) - - dt1 = dt(steps) + + dt1 = dt(steps) ds3 = ds0.divide(2) steps.append(timer()) - print ("ds3 = ds0.divide(2)", dt(steps), 0.25 , ds3.as_array()[0][0][0]) + print("ds3 = ds0.divide(2)", dt(steps), 0.25, ds3.as_array()[0][0][0]) dt2 = dt(steps) - self.assertLess(dt1,dt2) + self.assertLess(dt1, dt2) self.assertEqual(.25, ds3.as_array()[0][0][0]) self.assertEqual(.5, ds.as_array()[0][0][0]) - - + def test_creation_copy(self): - shape = (2,3,4,5) + shape = (2, 3, 4, 5) size = shape[0] for i in range(1, len(shape)): size = size * shape[i] #print("a refcount " , sys.getrefcount(a)) - a = numpy.asarray([i for i in range( size )]) + a = numpy.asarray([i for i in range(size)]) #print("a refcount " , sys.getrefcount(a)) a = numpy.reshape(a, shape) #print("a refcount " , sys.getrefcount(a)) - ds = DataContainer(a, True, ['X', 'Y','Z' ,'W']) + ds = DataContainer(a, True, ['X', 'Y', 'Z', 'W']) #print("a refcount " , sys.getrefcount(a)) - self.assertEqual(sys.getrefcount(a),2) - + self.assertEqual(sys.getrefcount(a), 2) + def test_subset(self): - shape = (4,3,2) + shape = (4, 3, 2) a = [i for i in range(2*3*4)] a = numpy.asarray(a) a = numpy.reshape(a, shape) - ds = DataContainer(a, True, ['X', 'Y','Z']) + ds = DataContainer(a, True, ['X', 'Y', 'Z']) sub = ds.subset(['X']) res = True try: numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([0,6,12,18])) + numpy.asarray([0, 6, 12, 18])) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - + sub = ds.subset(['X'], Y=2, Z=0) res = True try: numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([4,10,16,22])) + numpy.asarray([4, 10, 16, 22])) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - - + sub = ds.subset(['Y']) try: numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0,2,4])) + sub.as_array(), numpy.asarray([0, 2, 4])) res = True except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - - + sub = ds.subset(['Z']) try: numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([0,1])) + sub.as_array(), numpy.asarray([0, 1])) res = True except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) sub = ds.subset(['Z'], X=1, Y=2) try: numpy.testing.assert_array_equal( - sub.as_array(), numpy.asarray([10,11])) + sub.as_array(), numpy.asarray([10, 11])) res = True except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - + print(a) - sub = ds.subset(['X', 'Y'] , Z=1) + sub = ds.subset(['X', 'Y'], Z=1) res = True try: numpy.testing.assert_array_equal(sub.as_array(), - numpy.asarray([[ 1, 3, 5], - [ 7, 9, 11], - [13, 15, 17], - [19, 21, 23]])) + numpy.asarray([[1, 3, 5], + [7, 9, 11], + [13, 15, 17], + [19, 21, 23]])) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - + def test_ImageData(self): # create ImageData from geometry vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2) vol = ImageData(geometry=vgeometry) - self.assertEqual(vol.shape , (2,3,4)) - + self.assertEqual(vol.shape, (2, 3, 4)) + vol1 = vol + 1 self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape)) - + vol1 = vol - 1 self.assertNumpyArrayEqual(vol1.as_array(), -numpy.ones(vol.shape)) - + vol1 = 2 * (vol + 1) self.assertNumpyArrayEqual(vol1.as_array(), 2 * numpy.ones(vol.shape)) - - vol1 = (vol + 1) / 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2 ) - + + vol1 = (vol + 1) / 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) / 2) + vol1 = vol + 1 - self.assertEqual(vol1.sum() , 2*3*4) - vol1 = ( vol + 2 ) ** 2 - self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4 ) - - - + self.assertEqual(vol1.sum(), 2*3*4) + vol1 = (vol + 2) ** 2 + self.assertNumpyArrayEqual(vol1.as_array(), numpy.ones(vol.shape) * 4) + + self.assertEqual(vol.number_of_dimensions, 3) + def test_AcquisitionData(self): - sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), - geom_type='parallel', pixel_num_v=3, - pixel_num_h=5 , channels=2) + sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10), + geom_type='parallel', pixel_num_v=3, + pixel_num_h=5, channels=2) sino = AcquisitionData(geometry=sgeometry) - self.assertEqual(sino.shape , (2,10,3,5)) - - + self.assertEqual(sino.shape, (2, 10, 3, 5)) + def assertNumpyArrayEqual(self, first, second): res = True try: numpy.testing.assert_array_equal(first, second) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): res = True try: numpy.testing.assert_array_almost_equal(first, second, decimal) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) + def test_DataContainerChaining(self): + dc = self.create_DataContainer(256,256,256,1) + + dc.add(9,out=dc)\ + .subtract(1,out=dc) + self.assertEqual(1+9-1,dc.as_array().flatten()[0]) + + + class TestAlgorithms(unittest.TestCase): def assertNumpyArrayEqual(self, first, second): @@ -473,311 +505,372 @@ class TestAlgorithms(unittest.TestCase): numpy.testing.assert_array_equal(first, second) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - def assertNumpyArrayAlmostEqual(self, first, second, decimal=6): res = True try: numpy.testing.assert_array_almost_equal(first, second, decimal) except AssertionError as err: res = False - print (err) + print(err) self.assertTrue(res) - def test_FISTA(self): - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0],1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt':True} - # Create object instances with the test data A and b. - f = Norm2sq(A,b,c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n,1))) - - f.grad(x_init) - - # Run FISTA for least squares plus zero function. - x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0 , opt=opt) - - # Print solution and final objective/criterion value for comparison - print("FISTA least squares plus zero function solution and objective value:") - print(x_fista0.array) - print(criter0[-1]) - - # Compare to CVXPY - - # Construct the problem. - x0 = Variable(n) - objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]) ) - prob0 = Problem(objective0) - - # The optimal objective is returned by prob.solve(). - result0 = prob0.solve(verbose=False,solver=SCS,eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus zero function solution and objective value:") - print(x0.value) - print(objective0.value) - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista0.array),x0.value,6) - def test_FISTA_Norm1(self): - - opt = {'memopt':True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0],1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt':True} - # Create object instances with the test data A and b. - f = Norm2sq(A,b,c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n,1))) - - # Create 1-norm object instance - g1 = Norm1(lam) - - g1(x_init) - g1.prox(x_init,0.02) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1,opt=opt) - - # Print for comparison - print("FISTA least squares plus 1-norm solution and objective value:") - print(x_fista1.as_array().squeeze()) - print(criter1[-1]) - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) ) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False,solver=SCS,eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fista1.array),x1.value,6) - - def test_FBPD_Norm1(self): - - opt = {'memopt':True} - # Problem data. - m = 30 - n = 20 - np.random.seed(1) - Amat = np.random.randn(m, n) - A = LinearOperatorMatrix(Amat) - bmat = np.random.randn(m) - bmat.shape = (bmat.shape[0],1) - - # A = Identity() - # Change n to equal to m. - - b = DataContainer(bmat) - - # Regularization parameter - lam = 10 - opt = {'memopt':True} - # Create object instances with the test data A and b. - f = Norm2sq(A,b,c=0.5, memopt=True) - g0 = ZeroFun() - - # Initial guess - x_init = DataContainer(np.zeros((n,1))) - - # Create 1-norm object instance - g1 = Norm1(lam) - - - # Compare to CVXPY - - # Construct the problem. - x1 = Variable(n) - objective1 = Minimize(0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1,1) ) - prob1 = Problem(objective1) - - # The optimal objective is returned by prob.solve(). - result1 = prob1.solve(verbose=False,solver=SCS,eps=1e-9) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1.value) - print(objective1.value) - - # Now try another algorithm FBPD for same problem: - x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, - Identity(), None, f, g1) - print(x_fbpd1) - print(criterfbpd1[-1]) - - self.assertNumpyArrayAlmostEqual( - numpy.squeeze(x_fbpd1.array),x1.value,6) + + def test_FISTA_cvx(self): + if not cvx_not_installable: + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + f.grad(x_init) + + # Run FISTA for least squares plus zero function. + x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt) + + # Print solution and final objective/criterion value for comparison + print("FISTA least squares plus zero function solution and objective value:") + print(x_fista0.array) + print(criter0[-1]) + + # Compare to CVXPY + + # Construct the problem. + x0 = Variable(n) + objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0])) + prob0 = Problem(objective0) + + # The optimal objective is returned by prob.solve(). + result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus zero function solution and objective value:") + print(x0.value) + print(objective0.value) + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista0.array), x0.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def test_FISTA_Norm1_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + # Create 1-norm object instance + g1 = Norm1(lam) + + g1(x_init) + g1.prox(x_init, 0.02) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g1, opt=opt) + + # Print for comparison + print("FISTA least squares plus 1-norm solution and objective value:") + print(x_fista1.as_array().squeeze()) + print(criter1[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fista1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) + + def test_FBPD_Norm1_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + # Problem data. + m = 30 + n = 20 + np.random.seed(1) + Amat = np.random.randn(m, n) + A = LinearOperatorMatrix(Amat) + bmat = np.random.randn(m) + bmat.shape = (bmat.shape[0], 1) + + # A = Identity() + # Change n to equal to m. + + b = DataContainer(bmat) + + # Regularization parameter + lam = 10 + opt = {'memopt': True} + # Create object instances with the test data A and b. + f = Norm2sq(A, b, c=0.5, memopt=True) + g0 = ZeroFun() + + # Initial guess + x_init = DataContainer(np.zeros((n, 1))) + + # Create 1-norm object instance + g1 = Norm1(lam) + + # Compare to CVXPY + + # Construct the problem. + x1 = Variable(n) + objective1 = Minimize( + 0.5*sum_squares(Amat*x1 - bmat.T[0]) + lam*norm(x1, 1)) + prob1 = Problem(objective1) + + # The optimal objective is returned by prob.solve(). + result1 = prob1.solve(verbose=False, solver=SCS, eps=1e-9) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1.value) + print(objective1.value) + + # Now try another algorithm FBPD for same problem: + x_fbpd1, itfbpd1, timingfbpd1, criterfbpd1 = FBPD(x_init, + Identity(), None, f, g1) + print(x_fbpd1) + print(criterfbpd1[-1]) + + self.assertNumpyArrayAlmostEqual( + numpy.squeeze(x_fbpd1.array), x1.value, 6) + else: + self.assertTrue(cvx_not_installable) # Plot criterion curve to see both FISTA and FBPD converge to same value. # Note that FISTA is very efficient for 1-norm minimization so it beats - # FBPD in this test by a lot. But FBPD can handle a larger class of problems + # FBPD in this test by a lot. But FBPD can handle a larger class of problems # than FISTA can. - - + # Now try 1-norm and TV denoising with FBPD, first 1-norm. - - # Set up phantom size NxN by creating ImageGeometry, initialising the + + # Set up phantom size NxN by creating ImageGeometry, initialising the # ImageData object with this geometry and empty array and finally put some # data into its array, and display as image. - def test_FISTA_denoise(self): - opt = {'memopt':True} + def test_FISTA_denoise_cvx(self): + if not cvx_not_installable: + opt = {'memopt': True} + N = 64 + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) + Phantom = ImageData(geometry=ig) + + x = Phantom.as_array() + + x[int(round(N/4)):int(round(3*N/4)), + int(round(N/4)):int(round(3*N/4))] = 0.5 + x[int(round(N/8)):int(round(7*N/8)), + int(round(3*N/8)):int(round(5*N/8))] = 1 + + # Identity operator for denoising + I = TomoIdentity(ig) + + # Data and add noise + y = I.direct(Phantom) + y.array = y.array + 0.1*np.random.randn(N, N) + + # Data fidelity term + f_denoise = Norm2sq(I, y, c=0.5, memopt=True) + + # 1-norm regulariser + lam1_denoise = 1.0 + g1_denoise = Norm1(lam1_denoise) + + # Initial guess + x_init_denoise = ImageData(np.zeros((N, N))) + + # Combine with least squares and solve using generic FISTA implementation + x_fista1_denoise, it1_denoise, timing1_denoise, \ + criter1_denoise = \ + FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) + + print(x_fista1_denoise) + print(criter1_denoise[-1]) + + # Now denoise LS + 1-norm with FBPD + x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ + criterfbpd1_denoise = \ + FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) + print(x_fbpd1_denoise) + print(criterfbpd1_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + x1_denoise = Variable(N**2, 1) + objective1_denoise = Minimize( + 0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise, 1)) + prob1_denoise = Problem(objective1_denoise) + + # The optimal objective is returned by prob.solve(). + result1_denoise = prob1_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(x1_denoise.value) + print(objective1_denoise.value) + self.assertNumpyArrayAlmostEqual( + x_fista1_denoise.array.flatten(), x1_denoise.value, 5) + + self.assertNumpyArrayAlmostEqual( + x_fbpd1_denoise.array.flatten(), x1_denoise.value, 5) + x1_cvx = x1_denoise.value + x1_cvx.shape = (N, N) + + # Now TV with FBPD + lam_tv = 0.1 + gtv = TV2D(lam_tv) + gtv(gtv.op.direct(x_init_denoise)) + + opt_tv = {'tol': 1e-4, 'iter': 10000} + + x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ + criterfbpdtv_denoise = \ + FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv, opt=opt_tv) + print(x_fbpdtv_denoise) + print(criterfbpdtv_denoise[-1]) + + # Compare to CVXPY + + # Construct the problem. + xtv_denoise = Variable((N, N)) + objectivetv_denoise = Minimize( + 0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise)) + probtv_denoise = Problem(objectivetv_denoise) + + # The optimal objective is returned by prob.solve(). + resulttv_denoise = probtv_denoise.solve( + verbose=False, solver=SCS, eps=1e-12) + + # The optimal solution for x is stored in x.value and optimal objective value + # is in result as well as in objective.value + print("CVXPY least squares plus 1-norm solution and objective value:") + print(xtv_denoise.value) + print(objectivetv_denoise.value) + + self.assertNumpyArrayAlmostEqual( + x_fbpdtv_denoise.as_array(), xtv_denoise.value, 1) + + else: + self.assertTrue(cvx_not_installable) + + +class TestFunction(unittest.TestCase): + def assertNumpyArrayEqual(self, first, second): + res = True + try: + numpy.testing.assert_array_equal(first, second) + except AssertionError as err: + res = False + print(err) + self.assertTrue(res) + + def create_simple_ImageData(self): N = 64 - ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N) + ig = ImageGeometry(voxel_num_x=N, voxel_num_y=N) Phantom = ImageData(geometry=ig) - + x = Phantom.as_array() - + x[int(round(N/4)):int(round(3*N/4)), int(round(N/4)):int(round(3*N/4))] = 0.5 x[int(round(N/8)):int(round(7*N/8)), int(round(3*N/8)):int(round(5*N/8))] = 1 - - - # Identity operator for denoising - I = TomoIdentity(ig) - - # Data and add noise - y = I.direct(Phantom) - y.array = y.array + 0.1*np.random.randn(N, N) - - - # Data fidelity term - f_denoise = Norm2sq(I,y,c=0.5,memopt=True) - - # 1-norm regulariser - lam1_denoise = 1.0 - g1_denoise = Norm1(lam1_denoise) - - # Initial guess - x_init_denoise = ImageData(np.zeros((N,N))) - - # Combine with least squares and solve using generic FISTA implementation - x_fista1_denoise, it1_denoise, timing1_denoise, \ - criter1_denoise = \ - FISTA(x_init_denoise, f_denoise, g1_denoise, opt=opt) - - print(x_fista1_denoise) - print(criter1_denoise[-1]) - - - # Now denoise LS + 1-norm with FBPD - x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise,\ - criterfbpd1_denoise = \ - FBPD(x_init_denoise, I, None, f_denoise, g1_denoise) - print(x_fbpd1_denoise) - print(criterfbpd1_denoise[-1]) - - - # Compare to CVXPY - - # Construct the problem. - x1_denoise = Variable(N**2,1) - objective1_denoise = Minimize(0.5*sum_squares(x1_denoise - y.array.flatten()) + lam1_denoise*norm(x1_denoise,1) ) - prob1_denoise = Problem(objective1_denoise) - - # The optimal objective is returned by prob.solve(). - result1_denoise = prob1_denoise.solve(verbose=False,solver=SCS,eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(x1_denoise.value) - print(objective1_denoise.value) - self.assertNumpyArrayAlmostEqual( - x_fista1_denoise.array.flatten(),x1_denoise.value,5) - - self.assertNumpyArrayAlmostEqual( - x_fbpd1_denoise.array.flatten(),x1_denoise.value,5) - x1_cvx = x1_denoise.value - x1_cvx.shape = (N,N) - - - # Now TV with FBPD - lam_tv = 0.1 - gtv = TV2D(lam_tv) - gtv(gtv.op.direct(x_init_denoise)) - - opt_tv = {'tol': 1e-4, 'iter': 10000} - - x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise,\ - criterfbpdtv_denoise = \ - FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv,opt=opt_tv) - print(x_fbpdtv_denoise) - print(criterfbpdtv_denoise[-1]) - - - - # Compare to CVXPY - - # Construct the problem. - xtv_denoise = Variable((N,N)) - objectivetv_denoise = Minimize(0.5*sum_squares(xtv_denoise - y.array) + lam_tv*tv(xtv_denoise) ) - probtv_denoise = Problem(objectivetv_denoise) - - # The optimal objective is returned by prob.solve(). - resulttv_denoise = probtv_denoise.solve(verbose=False,solver=SCS,eps=1e-12) - - # The optimal solution for x is stored in x.value and optimal objective value - # is in result as well as in objective.value - print("CVXPY least squares plus 1-norm solution and objective value:") - print(xtv_denoise.value) - print(objectivetv_denoise.value) - - self.assertNumpyArrayAlmostEqual( - x_fbpdtv_denoise.as_array(),xtv_denoise.value,1) - + + return (ig, Phantom) + + def _test_Norm2(self): + print("test Norm2") + opt = {'memopt': True} + ig, Phantom = self.create_simple_ImageData() + x = Phantom.as_array() + print(Phantom) + print(Phantom.as_array()) + + norm2 = Norm2() + v1 = norm2(x) + v2 = norm2(Phantom) + self.assertEqual(v1, v2) + + p1 = norm2.prox(Phantom, 1) + print(p1) + p2 = norm2.prox(x, 1) + self.assertNumpyArrayEqual(p1.as_array(), p2) + + p3 = norm2.proximal(Phantom, 1) + p4 = norm2.proximal(x, 1) + self.assertNumpyArrayEqual(p3.as_array(), p2) + self.assertNumpyArrayEqual(p3.as_array(), p4) + + out = Phantom.copy() + p5 = norm2.proximal(Phantom, 1, out=out) + self.assertEqual(id(p5), id(out)) + self.assertNumpyArrayEqual(p5.as_array(), p3.as_array()) # ============================================================================= # def test_upper(self): # self.assertEqual('foo'.upper(), 'FOO') -# +# # def test_isupper(self): # self.assertTrue('FOO'.isupper()) # self.assertFalse('Foo'.isupper()) -# +# # def test_split(self): # s = 'hello world' # self.assertEqual(s.split(), ['hello', 'world']) @@ -786,5 +879,6 @@ class TestAlgorithms(unittest.TestCase): # s.split(2) # ============================================================================= + if __name__ == '__main__': unittest.main() diff --git a/Wrappers/Python/wip/demo_compare_cvx.py b/Wrappers/Python/wip/demo_compare_cvx.py index ad79fa5..27b1c97 100644 --- a/Wrappers/Python/wip/demo_compare_cvx.py +++ b/Wrappers/Python/wip/demo_compare_cvx.py @@ -1,10 +1,11 @@ from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, DataContainer from ccpi.optimisation.algs import FISTA, FBPD, CGLS -from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D +from ccpi.optimisation.funcs import Norm2sq, ZeroFun, Norm1, TV2D, Norm2 from ccpi.optimisation.ops import LinearOperatorMatrix, TomoIdentity from ccpi.optimisation.ops import Identity +from ccpi.optimisation.ops import FiniteDiff2D # Requires CVXPY, see http://www.cvxpy.org/ # CVXPY can be installed in anaconda using @@ -172,6 +173,8 @@ plt.imshow(y.array) plt.title('Noisy image') plt.show() + +################### # Data fidelity term f_denoise = Norm2sq(I,y,c=0.5,memopt=True) @@ -188,9 +191,9 @@ x_fista1_denoise, it1_denoise, timing1_denoise, criter1_denoise = FISTA(x_init_d print(x_fista1_denoise) print(criter1_denoise[-1]) -plt.imshow(x_fista1_denoise.as_array()) -plt.title('FISTA LS+1') -plt.show() +#plt.imshow(x_fista1_denoise.as_array()) +#plt.title('FISTA LS+1') +#plt.show() # Now denoise LS + 1-norm with FBPD x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, \ @@ -198,9 +201,9 @@ x_fbpd1_denoise, itfbpd1_denoise, timingfbpd1_denoise, \ print(x_fbpd1_denoise) print(criterfbpd1_denoise[-1]) -plt.imshow(x_fbpd1_denoise.as_array()) -plt.title('FBPD LS+1') -plt.show() +#plt.imshow(x_fbpd1_denoise.as_array()) +#plt.title('FBPD LS+1') +#plt.show() if use_cvxpy: # Compare to CVXPY @@ -222,25 +225,46 @@ if use_cvxpy: x1_cvx = x1_denoise.value x1_cvx.shape = (N,N) + + +#plt.imshow(x1_cvx) +#plt.title('CVX LS+1') +#plt.show() + +fig = plt.figure() +plt.subplot(1,4,1) +plt.imshow(y.array) +plt.title("LS+1") +plt.subplot(1,4,2) +plt.imshow(x_fista1_denoise.as_array()) +plt.title("fista") +plt.subplot(1,4,3) +plt.imshow(x_fbpd1_denoise.as_array()) +plt.title("fbpd") +plt.subplot(1,4,4) plt.imshow(x1_cvx) -plt.title('CVX LS+1') +plt.title("cvx") plt.show() -# Now TV with FBPD +############################################################## +# Now TV with FBPD and Norm2 lam_tv = 0.1 gtv = TV2D(lam_tv) -gtv(gtv.op.direct(x_init_denoise)) +norm2 = Norm2(lam_tv) +op = FiniteDiff2D() +#gtv(gtv.op.direct(x_init_denoise)) opt_tv = {'tol': 1e-4, 'iter': 10000} x_fbpdtv_denoise, itfbpdtv_denoise, timingfbpdtv_denoise, \ - criterfbpdtv_denoise = FBPD(x_init_denoise, gtv.op, None, f_denoise, gtv,opt=opt_tv) + criterfbpdtv_denoise = FBPD(x_init_denoise, op, None, \ + f_denoise, norm2 ,opt=opt_tv) print(x_fbpdtv_denoise) print(criterfbpdtv_denoise[-1]) plt.imshow(x_fbpdtv_denoise.as_array()) plt.title('FBPD TV') -plt.show() +#plt.show() if use_cvxpy: # Compare to CVXPY @@ -262,7 +286,21 @@ if use_cvxpy: plt.imshow(xtv_denoise.value) plt.title('CVX TV') +#plt.show() + +fig = plt.figure() +plt.subplot(1,3,1) +plt.imshow(y.array) +plt.title("TV2D") +plt.subplot(1,3,2) +plt.imshow(x_fbpdtv_denoise.as_array()) +plt.title("fbpd tv denoise") +plt.subplot(1,3,3) +plt.imshow(xtv_denoise.value) +plt.title("CVX tv") plt.show() + + plt.loglog([0,opt_tv['iter']], [objectivetv_denoise.value,objectivetv_denoise.value], label='CVX TV') plt.loglog(criterfbpdtv_denoise, label='FBPD TV') -- cgit v1.2.3