diff options
| -rw-r--r-- | Wrappers/Python/ccpi/framework.py | 96 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/funcs.py | 30 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/ops.py | 6 | ||||
| -rw-r--r-- | Wrappers/Python/conda-recipe/meta.yaml | 5 | ||||
| -rwxr-xr-x | Wrappers/Python/conda-recipe/run_test.py | 1012 | ||||
| -rw-r--r-- | Wrappers/Python/wip/demo_compare_cvx.py | 62 |
6 files changed, 673 insertions, 538 deletions
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') |