Test norm (#162)

* enable win build with unittest. skip cvx

* unittest on funcs

* fixed Norm1 and function

* test chaining

* a few fixes for unittest

* investigating Norm2

* restored Norm2

TODO: rewrite it and test

* removed comment

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')