summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--Wrappers/Python/ccpi/framework.py96
-rwxr-xr-xWrappers/Python/ccpi/optimisation/funcs.py30
-rwxr-xr-xWrappers/Python/ccpi/optimisation/ops.py6
-rw-r--r--Wrappers/Python/conda-recipe/meta.yaml5
-rwxr-xr-xWrappers/Python/conda-recipe/run_test.py1012
-rw-r--r--Wrappers/Python/wip/demo_compare_cvx.py62
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')
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-14 22:13:52 +0000