memopt test

8 files changed, 236 insertions, 188 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index 9664037..888950d 100755
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -187,7 +187,7 @@ class BlockDataContainer(object):
     def clone(self):

         return type(self)(*[el.copy() for el in self.containers], shape=self.shape)

     def fill(self, x):

-        for el,ot in zip(self.containers, x):

+        for el,ot in zip(self.containers, x.containers):

             el.fill(ot)

     

     def __add__(self, other):

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 46a1969..d07005a 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -6,8 +6,9 @@ Created on Mon Feb  4 16:18:06 2019
 @author: evangelos
 """
 from ccpi.optimisation.algorithms import Algorithm
-from ccpi.framework import ImageData
+from ccpi.framework import ImageData, DataContainer
 import numpy as np
+import numpy
 import time
 from ccpi.optimisation.operators import BlockOperator
 from ccpi.framework import BlockDataContainer
@@ -80,6 +81,29 @@ class PDHG(Algorithm):
             -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(- 1 * self.operator.adjoint(self.y)))
         ])
 
+    
+    
+def assertBlockDataContainerEqual( container1, container2):
+    print ("assert Block Data Container Equal")
+    assert issubclass(container1.__class__, container2.__class__)
+    for col in range(container1.shape[0]):
+        if issubclass(container1.get_item(col).__class__, DataContainer):
+            print ("Checking col ", col)
+            assertNumpyArrayEqual(
+                container1.get_item(col).as_array(), 
+                container2.get_item(col).as_array()
+                )
+        else:
+            assertBlockDataContainerEqual(container1.get_item(col),container2.get_item(col))    
+
+def assertNumpyArrayEqual(first, second):
+    res = True
+    try:
+        numpy.testing.assert_array_equal(first, second)
+    except AssertionError as err:
+        res = False
+        print(err)
+    assert res
 
 
 def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
@@ -102,12 +126,12 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
     x_old = operator.domain_geometry().allocate()
     y_old = operator.range_geometry().allocate()       
             
-    xbar = x_old
-    x_tmp = x_old
-    x = x_old
+    xbar = x_old.copy()
+    x_tmp = x_old.copy()
+    x = x_old.copy()
     
-    y_tmp = y_old
-    y = y_old
+    y_tmp = y_old.copy()
+    y = y_old.copy()
         
     # relaxation parameter
     theta = 1
@@ -137,37 +161,26 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
         
         else:
             
-            operator.direct(xbar, out = y_tmp) 
-            
-            y_tmp.multiply(sigma, out = y_tmp)
-            y_tmp.add(y_old, out = y_tmp)
-#            y_tmp.__imul__(sigma)
-#            y_tmp.__iadd__(y_old)
-            
-#            y_tmp *= sigma
-#            y_tmp += y_old
-            
-#            y_tmp = y_old + sigma * operator.direct(xbar)                        
+
+            operator.direct(xbar, out = y_tmp)             
+            y_tmp *= sigma
+            y_tmp += y_old                      
             f.proximal_conjugate(y_tmp, sigma, out=y)
-            
-#            x_tmp = x_old - tau * operator.adjoint(y)
-            
+
             operator.adjoint(y, out = x_tmp)   
-            x_tmp.multiply(-tau, out = x_tmp)
-            x_tmp.add(x_old, out = x_tmp)
-            
-            
-#            x_tmp *= -tau
-#            x_tmp += x_old
+            x_tmp *= -tau
+            x_tmp += x_old
 
             g.proximal(x_tmp, tau, out = x)
             
             xbar = x - x_old
             xbar *= theta
             xbar += x
-            
-            x_old.fill(x)
-            y_old.fill(y)
+                        
+                        
+            x_old = x.copy()
+            y_old = y.copy()
+
             
 #            pass
 #        
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index dd463c0..639f7bf 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -91,6 +91,8 @@ class MixedL21Norm(Function):
                 res = BlockDataContainer(*frac) 
                 return res
             else:
+                
+                
                 pass
                         
 #                tmp2 = np.sqrt(x.as_array()[0]**2 +  x.as_array()[1]**2 +  2*x.as_array()[2]**2)/self.alpha
@@ -101,19 +103,25 @@ class MixedL21Norm(Function):
                 
                 tmp = [ el*el for el in x]
                 res = (sum(tmp).sqrt()).maximum(1.0) 
-                frac = [x[i]/res for i in range(x.shape[0])]
-                res = BlockDataContainer(*frac)
                 
-                return res 
+#                res = sum(x**2).sqrt().maximum(1.0) 
+                
+                #frac = [x[i]/res for i in range(x.shape[0])]
+                #res = BlockDataContainer(*frac)
+                
+                return x/res 
             
             else:
                 
-                tmp = [ el*el for el in x]
+                tmp = x**2 #[ el*el for el in x]
                 res = (sum(tmp).sqrt()).maximum(1.0) 
-                frac = [x[i]/res for i in range(x.shape[0])]                
-#                res = (sum(x**2).sqrt()).maximum(1.0) 
-#                return x/res
-                out.fill(frac)
+                
+#                res = sum(x**2).sqrt().maximum(1.0)
+#                frac = [x[i]/res for i in range(x.shape[0])]
+#                res = BlockDataContainer(*frac)
+#                res = sum(x**2).sqrt().maximum(1.0) 
+                out.fill(x/res)
+#                return res 
 
 
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
index 88d9b64..6d21acb 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
@@ -45,14 +45,17 @@ class ZeroFun(Function):
             else:
                 return x.maximum(0).sum() + x.maximum(0).sum()
     
-    def proximal(self,x,tau, out=None):
+    def proximal(self, x, tau, out=None):
         if out is None:
             return x.copy()
         else:
             out.fill(x)
         
-    def proximal_conjugate(self, x, tau):
-        return 0
+    def proximal_conjugate(self, x, tau, out = None):
+        if out is None:
+            return 0
+        else:
+            return 0
 
     def domain_geometry(self):
         pass
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 3d2a96b..8e73ff2 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -52,38 +52,33 @@ class FiniteDiff(LinearOperator):
         
         if out is None:        
             out = np.zeros_like(x_asarr)
-            fd_arr = out
         else:
-            fd_arr = out.as_array()        
-        
-#        if out is None:        
-#            out = self.gm_domain.allocate().as_array()
-#        
-#        fd_arr = out.as_array()
-#        fd_arr = self.gm_domain.allocate().as_array()
-          
+            out = out.as_array()
+            out[:]=0
+
+                  
         ######################## Direct for 2D  ###############################
         if x_sz == 2:
             
             if self.direction == 1:
                 
-                np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,0:-1] )
+                np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = out[:,0:-1] )
                 
                 if self.bnd_cond == 'Neumann':
                     pass                                        
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,-1] )
+                    np.subtract( x_asarr[:,0], x_asarr[:,-1], out = out[:,-1] )
                 else: 
                     raise ValueError('No valid boundary conditions')
                 
             if self.direction == 0:
                 
-                np.subtract( x_asarr[1:], x_asarr[0:-1], out = fd_arr[0:-1,:] )
+                np.subtract( x_asarr[1:], x_asarr[0:-1], out = out[0:-1,:] )
 
                 if self.bnd_cond == 'Neumann':
                     pass                                        
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[-1,:] ) 
+                    np.subtract( x_asarr[0,:], x_asarr[-1,:], out = out[-1,:] ) 
                 else:    
                     raise ValueError('No valid boundary conditions') 
                     
@@ -92,35 +87,35 @@ class FiniteDiff(LinearOperator):
                     
             if self.direction == 0:  
                 
-                np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[0:-1,:,:] )
+                np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = out[0:-1,:,:] )
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[-1,:,:] ) 
+                    np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = out[-1,:,:] ) 
                 else:    
                     raise ValueError('No valid boundary conditions')                      
                                                              
             if self.direction == 1:
                 
-                np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,0:-1,:] ) 
+                np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = out[:,0:-1,:] ) 
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,-1,:] )
+                    np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = out[:,-1,:] )
                 else:    
                     raise ValueError('No valid boundary conditions')                      
                                 
              
             if self.direction == 2:
                 
-                np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,0:-1] ) 
+                np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = out[:,:,0:-1] ) 
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,-1] )
+                    np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = out[:,:,-1] )
                 else:    
                     raise ValueError('No valid boundary conditions')  
                     
@@ -128,42 +123,42 @@ class FiniteDiff(LinearOperator):
         elif x_sz == 4:
                     
             if self.direction == 0:                            
-                np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[0:-1,:,:,:] )
+                np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = out[0:-1,:,:,:] )
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[-1,:,:,:] )
+                    np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = out[-1,:,:,:] )
                 else:    
                     raise ValueError('No valid boundary conditions') 
                                                 
             if self.direction == 1:
-                np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,0:-1,:,:] ) 
+                np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = out[:,0:-1,:,:] ) 
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,-1,:,:] )
+                    np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = out[:,-1,:,:] )
                 else:    
                     raise ValueError('No valid boundary conditions')                 
                 
             if self.direction == 2:
-                np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,0:-1,:] ) 
+                np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = out[:,:,0:-1,:] ) 
                 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,-1,:] )
+                    np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = out[:,:,-1,:] )
                 else:    
                     raise ValueError('No valid boundary conditions')                   
                 
             if self.direction == 3:
-                np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,0:-1] )                 
+                np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = out[:,:,:,0:-1] )                 
 
                 if self.bnd_cond == 'Neumann':
                     pass
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,-1] )
+                    np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = out[:,:,:,-1] )
                 else:    
                     raise ValueError('No valid boundary conditions')                   
                                 
@@ -177,14 +172,18 @@ class FiniteDiff(LinearOperator):
     def adjoint(self, x, out=None):
         
         x_asarr = x.as_array()
-        #x_asarr = x
         x_sz = len(x.shape)
         
         if out is None:        
             out = np.zeros_like(x_asarr)
-            fd_arr = out
         else:
-            fd_arr = out.as_array()          
+            out = out.as_array()        
+        
+#        if out is None:        
+#            out = np.zeros_like(x_asarr)
+#            fd_arr = out
+#        else:
+#            fd_arr = out.as_array()          
         
 #        if out is None:        
 #            out = self.gm_domain.allocate().as_array()
@@ -198,28 +197,28 @@ class FiniteDiff(LinearOperator):
         
             if self.direction == 1:
                 
-                np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = fd_arr[:,1:] )
+                np.subtract( x_asarr[:,1:], x_asarr[:,0:-1], out = out[:,1:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[:,0], 0, out = fd_arr[:,0] )
-                    np.subtract( -x_asarr[:,-2], 0, out = fd_arr[:,-1] )
+                    np.subtract( x_asarr[:,0], 0, out = out[:,0] )
+                    np.subtract( -x_asarr[:,-2], 0, out = out[:,-1] )
                     
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,0], x_asarr[:,-1], out = fd_arr[:,0] )                                        
+                    np.subtract( x_asarr[:,0], x_asarr[:,-1], out = out[:,0] )                                        
                     
                 else:   
                     raise ValueError('No valid boundary conditions') 
                                     
             if self.direction == 0:
                 
-                np.subtract( x_asarr[1:,:], x_asarr[0:-1,:], out = fd_arr[1:,:] )
+                np.subtract( x_asarr[1:,:], x_asarr[0:-1,:], out = out[1:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[0,:], 0, out = fd_arr[0,:] )
-                    np.subtract( -x_asarr[-2,:], 0, out = fd_arr[-1,:] ) 
+                    np.subtract( x_asarr[0,:], 0, out = out[0,:] )
+                    np.subtract( -x_asarr[-2,:], 0, out = out[-1,:] ) 
                     
                 elif self.bnd_cond == 'Periodic':  
-                    np.subtract( x_asarr[0,:], x_asarr[-1,:], out = fd_arr[0,:] ) 
+                    np.subtract( x_asarr[0,:], x_asarr[-1,:], out = out[0,:] ) 
                     
                 else:   
                     raise ValueError('No valid boundary conditions')     
@@ -229,35 +228,35 @@ class FiniteDiff(LinearOperator):
                 
             if self.direction == 0:          
                   
-                np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = fd_arr[1:,:,:] )
+                np.subtract( x_asarr[1:,:,:], x_asarr[0:-1,:,:], out = out[1:,:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[0,:,:], 0, out = fd_arr[0,:,:] )
-                    np.subtract( -x_asarr[-2,:,:], 0, out = fd_arr[-1,:,:] )
+                    np.subtract( x_asarr[0,:,:], 0, out = out[0,:,:] )
+                    np.subtract( -x_asarr[-2,:,:], 0, out = out[-1,:,:] )
                 elif self.bnd_cond == 'Periodic':                     
-                    np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = fd_arr[0,:,:] )
+                    np.subtract( x_asarr[0,:,:], x_asarr[-1,:,:], out = out[0,:,:] )
                 else:   
                     raise ValueError('No valid boundary conditions')                     
                                     
             if self.direction == 1:
-                np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = fd_arr[:,1:,:] )
+                np.subtract( x_asarr[:,1:,:], x_asarr[:,0:-1,:], out = out[:,1:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[:,0,:], 0, out = fd_arr[:,0,:] )
-                    np.subtract( -x_asarr[:,-2,:], 0, out = fd_arr[:,-1,:] )
+                    np.subtract( x_asarr[:,0,:], 0, out = out[:,0,:] )
+                    np.subtract( -x_asarr[:,-2,:], 0, out = out[:,-1,:] )
                 elif self.bnd_cond == 'Periodic':                     
-                    np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = fd_arr[:,0,:] )
+                    np.subtract( x_asarr[:,0,:], x_asarr[:,-1,:], out = out[:,0,:] )
                 else:   
                     raise ValueError('No valid boundary conditions')                                 
                 
             if self.direction == 2:
-                np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = fd_arr[:,:,1:] )
+                np.subtract( x_asarr[:,:,1:], x_asarr[:,:,0:-1], out = out[:,:,1:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[:,:,0], 0, out = fd_arr[:,:,0] ) 
-                    np.subtract( -x_asarr[:,:,-2], 0, out = fd_arr[:,:,-1] ) 
+                    np.subtract( x_asarr[:,:,0], 0, out = out[:,:,0] ) 
+                    np.subtract( -x_asarr[:,:,-2], 0, out = out[:,:,-1] ) 
                 elif self.bnd_cond == 'Periodic':                     
-                    np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = fd_arr[:,:,0] )
+                    np.subtract( x_asarr[:,:,0], x_asarr[:,:,-1], out = out[:,:,0] )
                 else:   
                     raise ValueError('No valid boundary conditions')                                 
         
@@ -265,51 +264,51 @@ class FiniteDiff(LinearOperator):
         elif x_sz == 4:                
                 
             if self.direction == 0:                            
-                np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = fd_arr[1:,:,:,:] )
+                np.subtract( x_asarr[1:,:,:,:], x_asarr[0:-1,:,:,:], out = out[1:,:,:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[0,:,:,:], 0, out = fd_arr[0,:,:,:] )
-                    np.subtract( -x_asarr[-2,:,:,:], 0, out = fd_arr[-1,:,:,:] )
+                    np.subtract( x_asarr[0,:,:,:], 0, out = out[0,:,:,:] )
+                    np.subtract( -x_asarr[-2,:,:,:], 0, out = out[-1,:,:,:] )
                     
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = fd_arr[0,:,:,:] )
+                    np.subtract( x_asarr[0,:,:,:], x_asarr[-1,:,:,:], out = out[0,:,:,:] )
                 else:    
                     raise ValueError('No valid boundary conditions') 
                                 
             if self.direction == 1:
-                np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = fd_arr[:,1:,:,:] )
+                np.subtract( x_asarr[:,1:,:,:], x_asarr[:,0:-1,:,:], out = out[:,1:,:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                   np.subtract( x_asarr[:,0,:,:], 0, out = fd_arr[:,0,:,:] )
-                   np.subtract( -x_asarr[:,-2,:,:], 0, out = fd_arr[:,-1,:,:] )
+                   np.subtract( x_asarr[:,0,:,:], 0, out = out[:,0,:,:] )
+                   np.subtract( -x_asarr[:,-2,:,:], 0, out = out[:,-1,:,:] )
                     
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = fd_arr[:,0,:,:] )
+                    np.subtract( x_asarr[:,0,:,:], x_asarr[:,-1,:,:], out = out[:,0,:,:] )
                 else:    
                     raise ValueError('No valid boundary conditions') 
                     
                 
             if self.direction == 2:
-                np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = fd_arr[:,:,1:,:] )
+                np.subtract( x_asarr[:,:,1:,:], x_asarr[:,:,0:-1,:], out = out[:,:,1:,:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[:,:,0,:], 0, out = fd_arr[:,:,0,:] ) 
-                    np.subtract( -x_asarr[:,:,-2,:], 0, out = fd_arr[:,:,-1,:] ) 
+                    np.subtract( x_asarr[:,:,0,:], 0, out = out[:,:,0,:] ) 
+                    np.subtract( -x_asarr[:,:,-2,:], 0, out = out[:,:,-1,:] ) 
                     
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = fd_arr[:,:,0,:] )
+                    np.subtract( x_asarr[:,:,0,:], x_asarr[:,:,-1,:], out = out[:,:,0,:] )
                 else:    
                     raise ValueError('No valid boundary conditions')                 
                 
             if self.direction == 3:
-                np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = fd_arr[:,:,:,1:] )
+                np.subtract( x_asarr[:,:,:,1:], x_asarr[:,:,:,0:-1], out = out[:,:,:,1:] )
                 
                 if self.bnd_cond == 'Neumann':
-                    np.subtract( x_asarr[:,:,:,0], 0, out = fd_arr[:,:,:,0] ) 
-                    np.subtract( -x_asarr[:,:,:,-2], 0, out = fd_arr[:,:,:,-1] )   
+                    np.subtract( x_asarr[:,:,:,0], 0, out = out[:,:,:,0] ) 
+                    np.subtract( -x_asarr[:,:,:,-2], 0, out = out[:,:,:,-1] )   
                     
                 elif self.bnd_cond == 'Periodic':
-                    np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = fd_arr[:,:,:,0] )
+                    np.subtract( x_asarr[:,:,:,0], x_asarr[:,:,:,-1], out = out[:,:,:,0] )
                 else:    
                     raise ValueError('No valid boundary conditions')                  
                               
@@ -328,8 +327,7 @@ class FiniteDiff(LinearOperator):
         return self.gm_domain
        
     def norm(self):
-        x0 = self.gm_domain.allocate()
-        x0.fill( np.random.random_sample(x0.shape) )
+        x0 = self.gm_domain.allocate('random_int')
         self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0)
         return self.s1
     
@@ -337,23 +335,46 @@ class FiniteDiff(LinearOperator):
 if __name__ == '__main__':
     
     from ccpi.framework import ImageGeometry
+    import numpy
     
     N, M = 2, 3
 
     ig = ImageGeometry(N, M)
 
 
-    FD = FiniteDiff(ig, direction = 0, bnd_cond = 'Neumann')
+    FD = FiniteDiff(ig, direction = 1, bnd_cond = 'Neumann')
     u = FD.domain_geometry().allocate('random_int')
-    
-    
+        
     res = FD.domain_geometry().allocate()
+    res1 = FD.range_geometry().allocate()
     FD.direct(u, out=res)
-    print(res.as_array())
-#    z = FD.direct(u)
-    
+
+    z = FD.direct(u)    
 #    print(z.as_array(), res.as_array())
 
+    for i in range(10):
+#        
+        z1 = FD.direct(u) 
+        FD.direct(u, out=res)
+        
+        u = ig.allocate('random_int')
+        res = u
+        z1  = u
+        numpy.testing.assert_array_almost_equal(z1.as_array(), \
+                                                res.as_array(), decimal=4)
+        
+#        print(z1.as_array(), res.as_array())
+        z2 = FD.adjoint(z1) 
+        FD.adjoint(z1, out=res1)  
+        numpy.testing.assert_array_almost_equal(z2.as_array(), \
+                                                res1.as_array(), decimal=4)        
+        
+        
+        
+
+        
+        
+        
     
 #    w = G.range_geometry().allocate('random_int')
     
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 9c573cb..e535847 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -49,35 +49,32 @@ class Gradient(LinearOperator):
         
                 
         if out is not None:
+            
             for i in range(self.gm_range.shape[0]):
-                self.FD.direction=self.ind[i]
+                self.FD.direction = self.ind[i]
                 self.FD.direct(x, out = out[i])
-#                FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).direct(x, out=out[i])
-            return out
         else:
             tmp = self.gm_range.allocate()        
             for i in range(tmp.shape[0]):
                 self.FD.direction=self.ind[i]
                 tmp.get_item(i).fill(self.FD.direct(x))
-#                tmp.get_item(i).fill(FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).direct(x))
             return tmp    
         
     def adjoint(self, x, out=None):
         
         if out is not None:
           
-            tmp = self.gm_domain.allocate()            
+            tmp = self.gm_domain.allocate()   
+            
             for i in range(x.shape[0]):
                 self.FD.direction=self.ind[i] 
                 self.FD.adjoint(x.get_item(i), out = tmp)
-#                FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).adjoint(x.get_item(i), out=tmp)
                 out+=tmp               
         else:            
             tmp = self.gm_domain.allocate()
             for i in range(x.shape[0]):
                 self.FD.direction=self.ind[i]
                 tmp+=self.FD.adjoint(x.get_item(i))
-#                tmp+=FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).adjoint(x.get_item(i))
             return tmp    
             
     
@@ -96,7 +93,9 @@ class Gradient(LinearOperator):
     def __rmul__(self, scalar):
         return ScaledOperator(self, scalar) 
     
-    
+    ###########################################################################
+    ###############  For preconditioning ######################################
+    ###########################################################################
     def matrix(self):
         
         tmp = self.gm_range.allocate()
@@ -238,21 +237,4 @@ if __name__ == '__main__':
     
     
     
-#
-#
-#    res = G.range_geometry().allocate()
-##    
-#    G.direct(u, out=res)
-#    z = G.direct(u)
-##    
-#    print(res[0].as_array())
-#    print(z[0].as_array())
-#
-    
-    
-    
-    
-##    LHS = (G.direct(u)*w).sum()
-##    RHS = (u * G.adjoint(w)).sum()
-    
-#    
+
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py
index feb09ee..d0cb198 100755
--- a/Wrappers/Python/wip/pdhg_TV_denoising.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising.py
@@ -58,10 +58,10 @@ if method == '0':
 
     #### Create functions
       
-    f1 = alpha * MixedL21Norm()
-    f2 = 0.5 * L2NormSquared(b = noisy_data)
-    
-    f = BlockFunction(f1, f2 )                                        
+    f1 = MixedL21Norm()
+    f2 = L2NormSquared(b = noisy_data)    
+    f = BlockFunction(f1, f2)  
+                                      
     g = ZeroFun()
     
 else:
@@ -79,6 +79,7 @@ else:
 # Compute operator Norm
 normK = operator.norm()
 print ("normK", normK)
+
 # Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
@@ -86,9 +87,12 @@ tau = 1/(sigma*normK**2)
 opt = {'niter':100}
 opt1 = {'niter':100, 'memopt': True}
 
-res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
 res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
 
+print("with memopt \n")
+
+res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
+
 plt.figure(figsize=(5,5))
 plt.imshow(res.as_array())
 plt.colorbar()
diff --git a/Wrappers/Python/wip/test_profile.py b/Wrappers/Python/wip/test_profile.py
index a97ad8d..f14c0c3 100644
--- a/Wrappers/Python/wip/test_profile.py
+++ b/Wrappers/Python/wip/test_profile.py
@@ -10,58 +10,75 @@ Created on Mon Apr  8 13:57:46 2019
 
 from ccpi.framework import ImageGeometry
 from ccpi.optimisation.operators import Gradient, BlockOperator, Identity
+from ccpi.optimisation.functions import MixedL21Norm, L2NormSquared, BlockFunction
+import numpy
 
-N, M, K = 200, 300, 100
+N, M, K = 2, 3, 2
 
-ig = ImageGeometry(N, M, K)
+ig = ImageGeometry(N, M)
+b = ig.allocate('random_int')
 
 G = Gradient(ig)
 Id = Identity(ig)
 
-u = G.domain_geometry().allocate('random_int')
-w = G.range_geometry().allocate('random_int')
-
-
-res = G.range_geometry().allocate()
-res1 = G.domain_geometry().allocate()
-#
-#
-#LHS = (G.direct(u)*w).sum()
-#RHS = (u * G.adjoint(w)).sum()
-#
-#print(G.norm())
-#print(LHS, RHS)
-#
-##%%%re
-##
-#G.direct(u, out=res)
-#G.adjoint(w, out=res1)
-##
-#LHS1 = (res * w).sum()
-#RHS1 = (u * res1).sum()
-##
-#print(LHS1, RHS1)
-
-B = BlockOperator(2*G, 3*Id)
-uB = B.domain_geometry().allocate('random_int')
-resB = B.range_geometry().allocate()
-
-#z2 = B.direct(uB)
-#B.direct(uB, out = resB)
-
-#%%
+#operator = BlockOperator(G, Id)
+operator = G
+
+f1 = MixedL21Norm()
+f2 = L2NormSquared(b = b)
+
+f = BlockFunction( f1, f2)
+
+
+x_old = operator.domain_geometry().allocate()
+y_old = operator.range_geometry().allocate('random_int') 
+  
+
+xbar = operator.domain_geometry().allocate('random_int')
+
+x_tmp = x_old.copy()
+x = x_old.copy()
+    
+y_tmp = operator.range_geometry().allocate()
+y = y_old.copy()
+
+y1 = y.copy()
+
+sigma = 20
 
 for i in range(100):
-#    
-#    z2 = B.direct(uB)
-#    
-    B.direct(uB, out = resB)
     
-#    z1 = G.adjoint(w)
-#    z = G.direct(u)
+    operator.direct(xbar, out = y_tmp)  
+    y_tmp *= sigma
+    y_tmp += y_old
+
+    
+    y_tmp1 =  sigma * operator.direct(xbar) + y_old
+    
+    print(i)    
+    print(" y_old :", y_old[0].as_array(), "\n")
+    print(" y_tmp[0] :", y_tmp[0].as_array(),"\n")
+    print(" y_tmp1[0] :", y_tmp1[0].as_array())
+    
+    
+    numpy.testing.assert_array_equal(y_tmp[0].as_array(), \
+                                            y_tmp1[0].as_array())
+    
+    numpy.testing.assert_array_equal(y_tmp[1].as_array(), \
+                                            y_tmp1[1].as_array()) 
+    
+    
+    y1 = f.proximal_conjugate(y_tmp1, sigma)
+    f.proximal_conjugate(y_tmp, sigma, y)
+    
+
+    
+    
+ 
+    
+
+    
     
-#    G.adjoint(w, out=res1)
     
-#    G.direct(u, out=res)
     
     
\ No newline at end of file