del sym_grad

22 files changed, 595 insertions, 494 deletions

diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
index 386cd87..166014b 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/BlockDataContainer.py
@@ -270,6 +270,7 @@ class BlockDataContainer(object):
     def squared_norm(self):
         y = numpy.asarray([el.squared_norm() for el in self.containers])
         return y.sum() 
+        
     
     def norm(self):
         return numpy.sqrt(self.squared_norm())   
@@ -442,9 +443,12 @@ class BlockDataContainer(object):
         '''Inline truedivision'''
         return self.__idiv__(other)
     
+    def dot(self, other):
+#        
+        tmp = [ self.containers[i].dot(other.containers[i]) for i in range(self.shape[0])]
+        return sum(tmp)
     
-
-    
+       
     
 if __name__ == '__main__':
     
@@ -456,6 +460,7 @@ if __name__ == '__main__':
     BG = BlockGeometry(ig, ig)
     
     U = BG.allocate('random_int')
+    V = BG.allocate('random_int')
     
     
     print ("test sum BDC " )
@@ -468,10 +473,10 @@ if __name__ == '__main__':
     z1 = sum(U**2).sqrt().as_array()    
     numpy.testing.assert_array_equal(z, z1)   
     
-    
-
     z2 = U.pnorm(2)
     
+    zzz = U.dot(V)
+    
    
 
 
diff --git a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
index 5dd6750..ed44d99 100644
--- a/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
+++ b/Wrappers/Python/build/lib/ccpi/framework/BlockGeometry.py
@@ -31,7 +31,50 @@ class BlockGeometry(object):
         '''returns the Geometry in the BlockGeometry located at position index'''

         return self.geometries[index]            

 

-    def allocate(self, value=0, dimension_labels=None):

+    def allocate(self, value=0, dimension_labels=None, **kwargs):

+        

+        symmetry = kwargs.get('symmetry',False)        

         containers = [geom.allocate(value) for geom in self.geometries]

+        

+        if symmetry == True:

+                        

+            # for 2x2       

+            # [ ig11, ig12\

+            #   ig21, ig22]

+            

+            # Row-wise Order

+            

+            if len(containers)==4:

+                containers[1]=containers[2]

+            

+            # for 3x3  

+            # [ ig11, ig12, ig13\

+            #   ig21, ig22, ig23\

+            #   ig31, ig32, ig33]            

+                      

+            elif len(containers)==9:

+                containers[1]=containers[3]

+                containers[2]=containers[6]

+                containers[5]=containers[7]

+            

+            # for 4x4  

+            # [ ig11, ig12, ig13, ig14\

+            #   ig21, ig22, ig23, ig24\

+            #   ig31, ig32, ig33, ig34

+            #   ig41, ig42, ig43, ig44]   

+            

+            elif len(containers) == 16:

+                containers[1]=containers[4]

+                containers[2]=containers[8]

+                containers[3]=containers[12]

+                containers[6]=containers[9]

+                containers[7]=containers[10]

+                containers[11]=containers[15]

+                

+                

+                

+        

         return BlockDataContainer(*containers)

+    

+    

          

diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py
index ed95c3f..12cbabc 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/Algorithm.py
@@ -145,13 +145,40 @@ class Algorithm(object):
         if self.should_stop():
             print ("Stop cryterion has been reached.")
         i = 0
+        
+#        print("Iteration {:<5} Primal {:<5} Dual {:<5} PDgap".format('','',''))
         for _ in self:
-            if verbose and self.iteration % self.update_objective_interval == 0:
-                print ("Iteration {}/{}, objective {}".format(self.iteration, 
-                       self.max_iteration, self.get_last_objective()) )
-            else:
+            
+            
+            if self.iteration % self.update_objective_interval == 0:
+            
                 if callback is not None:
-                    callback(self.iteration, self.get_last_objective())
+                    callback(self.iteration, self.get_last_objective(), self.x)
+            
+                else:
+                    
+                    if verbose:
+            
+#            if verbose and self.iteration % self.update_objective_interval == 0:
+                #pass
+                # \t for tab
+#                print( "{:04}/{:04} {:<5} {:.4f} {:<5} {:.4f} {:<5} {:.4f}".\
+#                      format(self.iteration, self.max_iteration,'', \
+#                             self.get_last_objective()[0],'',\
+#                             self.get_last_objective()[1],'',\
+#                             self.get_last_objective()[2]))
+                
+                
+                        print ("Iteration {}/{}, {}".format(self.iteration, 
+                               self.max_iteration, self.get_last_objective()) )                
+                
+                #print ("Iteration {}/{}, Primal, Dual, PDgap = {}".format(self.iteration, 
+                #       self.max_iteration, self.get_last_objective()) )
+                
+                
+#                else:
+#                    if callback is not None:
+#                        callback(self.iteration, self.get_last_objective(), self.x)
             i += 1
             if i == iterations:
                 break
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py
index 8ea2b6c..ee51049 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/FISTA.py
@@ -65,10 +65,10 @@ class FISTA(Algorithm):
             
         # initialization
         if memopt:
-            self.y = x_init.clone()
-            self.x_old = x_init.clone()
-            self.x = x_init.clone()
-            self.u = x_init.clone()
+            self.y = x_init.copy()
+            self.x_old = x_init.copy()
+            self.x = x_init.copy()
+            self.u = x_init.copy()            
         else:
             self.x_old = x_init.copy()
             self.y = x_init.copy()
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
index a41bd48..0f5e8ef 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/algorithms/PDHG.py
@@ -13,6 +13,7 @@ import time
 from ccpi.optimisation.operators import BlockOperator
 from ccpi.framework import BlockDataContainer
 from ccpi.optimisation.functions import FunctionOperatorComposition
+import matplotlib.pyplot as plt
 
 class PDHG(Algorithm):
     '''Primal Dual Hybrid Gradient'''
@@ -46,34 +47,43 @@ class PDHG(Algorithm):
         self.y_old = self.operator.range_geometry().allocate()
         
         self.xbar = self.x_old.copy()
-        
+        self.x_tmp = self.x_old.copy()
         self.x = self.x_old.copy()
-        self.y = self.y_old.copy()
-        if self.memopt:
-            self.y_tmp = self.y_old.copy()
-            self.x_tmp = self.x_old.copy()
-        #y = y_tmp
+        
+        self.y_tmp = self.y_old.copy()        
+        self.y = self.y_tmp.copy()
+        
+        
+        
+        #self.y = self.y_old.copy()
+        
+        
+        #if self.memopt:
+        #    self.y_tmp = self.y_old.copy()
+        #    self.x_tmp = self.x_old.copy()
+
             
         # relaxation parameter
         self.theta = 1
 
     def update(self):
+        
         if self.memopt:
             # Gradient descent, Dual problem solution
             # self.y_old += self.sigma * self.operator.direct(self.xbar)
             self.operator.direct(self.xbar, out=self.y_tmp)
             self.y_tmp *= self.sigma
-            self.y_old += self.y_tmp
+            self.y_tmp += self.y_old
 
             #self.y = self.f.proximal_conjugate(self.y_old, self.sigma)
-            self.f.proximal_conjugate(self.y_old, self.sigma, out=self.y)
+            self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y)
 
             # Gradient ascent, Primal problem solution
             self.operator.adjoint(self.y, out=self.x_tmp)
-            self.x_tmp *= self.tau
-            self.x_old -= self.x_tmp
+            self.x_tmp *= -1*self.tau
+            self.x_tmp += self.x_old
 
-            self.g.proximal(self.x_old, self.tau, out=self.x)
+            self.g.proximal(self.x_tmp, self.tau, out=self.x)
 
             #Update
             self.x.subtract(self.x_old, out=self.xbar)
@@ -82,7 +92,8 @@ class PDHG(Algorithm):
 
             self.x_old.fill(self.x)
             self.y_old.fill(self.y)
-
+            
+         
         else:
             # Gradient descent, Dual problem solution
             self.y_old += self.sigma * self.operator.direct(self.xbar)
@@ -93,18 +104,28 @@ class PDHG(Algorithm):
             self.x = self.g.proximal(self.x_old, self.tau)
 
             #Update
-            #xbar = x + theta * (x - x_old)
-            self.xbar.fill(self.x)
-            self.xbar -= self.x_old 
+            self.x.subtract(self.x_old, out=self.xbar)
             self.xbar *= self.theta
             self.xbar += self.x
 
-            self.x_old = self.x
-            self.y_old = self.y
+            self.x_old.fill(self.x)
+            self.y_old.fill(self.y)
+            
+            #xbar = x + theta * (x - x_old)
+#            self.xbar.fill(self.x)
+#            self.xbar -= self.x_old 
+#            self.xbar *= self.theta
+#            self.xbar += self.x
+
+#            self.x_old.fill(self.x)
+#            self.y_old.fill(self.y)
+            
+          
 
     def update_objective(self):
+        
         p1 = self.f(self.operator.direct(self.x)) + self.g(self.x)
-        d1 = -(self.f.convex_conjugate(self.y) + self.g(-1*self.operator.adjoint(self.y)))
+        d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y)))
 
         self.loss.append([p1,d1,p1-d1])
 
@@ -151,8 +172,8 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
 
         
         if not memopt:
-            
-            y_tmp = y_old +  sigma * operator.direct(xbar)            
+                   
+            y_tmp = y_old +  sigma * operator.direct(xbar)                        
             y = f.proximal_conjugate(y_tmp, sigma)
                 
             x_tmp = x_old - tau*operator.adjoint(y)
@@ -161,20 +182,19 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
             x.subtract(x_old, out=xbar)
             xbar *= theta
             xbar += x
+            
+            if i%50==0:
+                
+                p1 = f(operator.direct(x)) + g(x)
+                d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) )            
+                primal.append(p1)
+                dual.append(d1)
+                pdgap.append(p1-d1) 
+                print(p1, d1, p1-d1)            
                                                   
             x_old.fill(x)
             y_old.fill(y)
-            
-            
-            if i%10==0:
-#                
-                p1 = f(operator.direct(x)) + g(x)
-                print(p1)
-#                d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
-#                primal.append(p1)
-#                dual.append(d1)
-#                pdgap.append(p1-d1)      
-#                print(p1, d1, p1-d1)
+                        
             
         else:
             
@@ -184,7 +204,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
             f.proximal_conjugate(y_tmp, sigma, out=y)
 
             operator.adjoint(y, out = x_tmp)   
-            x_tmp *= -tau
+            x_tmp *= -1*tau
             x_tmp += x_old
 
             g.proximal(x_tmp, tau, out = x)
@@ -192,19 +212,20 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
             x.subtract(x_old, out=xbar)
             xbar *= theta
             xbar += x
-                                                
-            x_old.fill(x)
-            y_old.fill(y)
             
-            if i%10==0:
-#                
+            if i%50==0:
+                
                 p1 = f(operator.direct(x)) + g(x)
-                print(p1)
-#                d1 = - ( f.convex_conjugate(y) + g(-1*operator.adjoint(y)) )
-#                primal.append(p1)
-#                dual.append(d1)
-#                pdgap.append(p1-d1)      
-#                print(p1, d1, p1-d1)
+                d1 = - ( f.convex_conjugate(y) + g.convex_conjugate(-1*operator.adjoint(y)) )  
+                primal.append(p1)
+                dual.append(d1)
+                pdgap.append(p1-d1)  
+                print(p1, d1, p1-d1)            
+                                    
+            x_old.fill(x)
+            y_old.fill(y)
+                    
+
             
                          
     t_end = time.time()        
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py b/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py
index efc465c..b2b9791 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/funcs.py
@@ -17,7 +17,6 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
-from ccpi.optimisation.ops import Identity, FiniteDiff2D
 import numpy
 from ccpi.framework import DataContainer
 import warnings
@@ -99,12 +98,6 @@ class Norm2(Function):
                 raise ValueError ('Wrong size: x{0} out{1}'.format(x.shape,out.shape) )
         
 
-class TV2D(Norm2):
-    
-    def __init__(self, gamma):
-        super(TV2D,self).__init__(gamma, 0)
-        self.op = FiniteDiff2D()
-        self.L = self.op.get_max_sing_val()
         
 
 # Define a class for squared 2-norm
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
index bf627a5..0836324 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/BlockFunction.py
@@ -93,16 +93,28 @@ class BlockFunction(Function):
         
         
         '''
+        
+        if out is None:
                 
-        out = [None]*self.length
-        if isinstance(tau, Number):
-            for i in range(self.length):
-                out[i] = self.functions[i].proximal(x.get_item(i), tau)
+            out = [None]*self.length
+            if isinstance(tau, Number):
+                for i in range(self.length):
+                    out[i] = self.functions[i].proximal(x.get_item(i), tau)
+            else:
+                for i in range(self.length):
+                    out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i))                        
+            
+            return BlockDataContainer(*out)  
+        
         else:
-            for i in range(self.length):
-                out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(i))                        
-
-        return BlockDataContainer(*out)     
+            if isinstance(tau, Number):
+                for i in range(self.length):
+                    self.functions[i].proximal(x.get_item(i), tau, out[i])
+            else:
+                for i in range(self.length):
+                    self.functions[i].proximal(x.get_item(i), tau.get_item(i), out[i])
+            
+            
     
     def gradient(self,x, out=None):
         
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py
index 70511bb..8895f3d 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/FunctionOperatorComposition.py
@@ -19,16 +19,13 @@ class FunctionOperatorComposition(Function):
     
     '''
     
-    def __init__(self, operator, function):
+    def __init__(self, function, operator):
         
         super(FunctionOperatorComposition, self).__init__()
+        
         self.function = function     
         self.operator = operator
-        alpha = 1
-        
-        if isinstance (function, ScaledFunction):
-            alpha = function.scalar
-        self.L = 2 * alpha * operator.norm()**2
+        self.L = function.L * operator.norm()**2 
         
         
     def __call__(self, x):
@@ -39,47 +36,57 @@ class FunctionOperatorComposition(Function):
         
         '''
     
-        return self.function(self.operator.direct(x))   
+        return self.function(self.operator.direct(x))  
+    
+    def gradient(self, x, out=None):
+#        
+        ''' Gradient takes into account the Operator'''
+        if out is None:
+            return self.operator.adjoint(self.function.gradient(self.operator.direct(x)))
+        else: 
+            tmp = self.operator.range_geometry().allocate()
+            self.operator.direct(x, out=tmp)
+            self.function.gradient(tmp, out=tmp)
+            self.operator.adjoint(tmp, out=out)
 
-    #TODO do not know if we need it
-    def call_adjoint(self, x):
     
-        return self.function(self.operator.adjoint(x))  
+    
 
+    #TODO do not know if we need it
+    #def call_adjoint(self, x):
+    #
+    #    return self.function(self.operator.adjoint(x))  
 
-    def convex_conjugate(self, x):
-        
-        ''' convex_conjugate does not take into account the Operator'''
-        return self.function.convex_conjugate(x)
 
-    def proximal(self, x, tau, out=None):
-        
-        '''proximal does not take into account the Operator'''                
-        if out is None:
-            return self.function.proximal(x, tau)
-        else:
-            self.function.proximal(x, tau, out=out)
-            
+    #def convex_conjugate(self, x):
+    #    
+    #        ''' convex_conjugate does not take into account the Operator'''
+    #    return self.function.convex_conjugate(x)
 
-    def proximal_conjugate(self, x, tau, out=None):    
+    
 
-        ''' proximal conjugate does not take into account the Operator'''
-        if out is None:
-            return self.function.proximal_conjugate(x, tau)
-        else:
-            self.function.proximal_conjugate(x, tau, out=out) 
 
-    def gradient(self, x, out=None):
-        
-        ''' Gradient takes into account the Operator'''
-        if out is None:
-            return self.operator.adjoint(
-                self.function.gradient(self.operator.direct(x))
-                )
-        else:
-            self.operator.adjoint(
-                self.function.gradient(self.operator.direct(x), 
-                out=out)
-            )
+                
+if __name__ == '__main__':   
+
+    from ccpi.framework import ImageGeometry
+    from ccpi.optimisation.operators import Gradient
+    from ccpi.optimisation.functions import L2NormSquared
+    
+    M, N, K = 2,3
+    ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
+    
+    G = Gradient(ig)
+    alpha = 0.5
+    
+    f = L2NormSquared()    
+    f_comp = FunctionOperatorComposition(G, alpha * f)
+    x = ig.allocate('random_int')
+    print(f_comp.gradient(x).shape
+          
+          )
+    
+
+             
         
                        
\ No newline at end of file
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py
index 7889cad..b53f669 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/KullbackLeibler.py
@@ -20,7 +20,8 @@
 import numpy
 from ccpi.optimisation.functions import Function
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction 
-from ccpi.framework import ImageData
+from ccpi.framework import ImageData, ImageGeometry
+import functools
 
 class KullbackLeibler(Function):
     
@@ -37,16 +38,23 @@ class KullbackLeibler(Function):
 
                                                 
     def __call__(self, x):
-        
-        # TODO check
                 
-        tmp = x + self.bnoise  
-        ind = tmp.as_array()>0
+        res_tmp = numpy.zeros(x.shape)
+        
+        tmp = x + self.bnoise 
+        ind = x.as_array()>0
+        
+        res_tmp[ind] = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind])
+        
+        return res_tmp.sum()
+
+    def log(self, datacontainer):
+        '''calculates the in-place log of the datacontainer'''
+        if not functools.reduce(lambda x,y: x and y>0,
+                                datacontainer.as_array().ravel(), True):
+            raise ValueError('KullbackLeibler. Cannot calculate log of negative number')
+        datacontainer.fill( numpy.log(datacontainer.as_array()) )
 
-        res = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind])
-            
-        return sum(res)
-    
         
     def gradient(self, x, out=None):
         
@@ -54,27 +62,42 @@ class KullbackLeibler(Function):
         if out is None:
             return 1 - self.b/(x + self.bnoise)
         else:
-            self.b.divide(x+self.bnoise, out=out)
+            x.add(self.bnoise, out=out)
+            self.b.divide(out, out=out)
             out.subtract(1, out=out)
-    
+            out *= -1
+            
     def convex_conjugate(self, x):
         
-        tmp = self.b/( 1 - x ) 
+        tmp = self.b/(1-x)
         ind = tmp.as_array()>0
         
-        sel
-        
-        return (self.b * ( ImageData( numpy.log(tmp) ) - 1 ) - self.bnoise * (x - 1)).sum()
+        return (self.b.as_array()[ind] * (numpy.log(tmp.as_array()[ind])-1)).sum()
+
     
     def proximal(self, x, tau, out=None):
         
         if out is None:        
             return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt() )
         else:
-            tmp =  0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt() )
-            out.fill(tmp)
             
-    
+            tmp =  0.5 *( (x - self.bnoise - tau) + 
+                        ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt()
+                        )
+            x.add(self.bnoise, out=out)
+            out -= tau
+            out *= out
+            tmp = self.b * (4 * tau)
+            out.add(tmp, out=out)
+            out.sqrt(out=out)
+            
+            x.subtract(self.bnoise, out=tmp)
+            tmp -= tau
+            
+            out += tmp
+            
+            out *= 0.5
+                            
     def proximal_conjugate(self, x, tau, out=None):
 
                 
@@ -82,10 +105,21 @@ class KullbackLeibler(Function):
             z = x + tau * self.bnoise
             return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
         else:
-            z = x + tau * self.bnoise
-            res = 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
-            out.fill(res)
+#            z = x + tau * self.bnoise
+#            out.fill( 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt()) )
+                        
+            tmp1 = x + tau * self.bnoise - 1
+            tmp2 = tmp1 + 2
             
+            self.b.multiply(4*tau, out=out)       
+            tmp1.multiply(tmp1, out=tmp1)
+            out += tmp1
+            out.sqrt(out=out)
+                        
+            out *= -1
+            out += tmp2
+            out *= 0.5
+
         
     
     def __rmul__(self, scalar):
@@ -106,6 +140,7 @@ if __name__ == '__main__':
     
     from ccpi.framework import ImageGeometry
     import numpy
+    
     N, M = 2,3
     ig  = ImageGeometry(N, M)
     data = ImageData(numpy.random.randint(-10, 10, size=(M, N)))
@@ -137,4 +172,4 @@ if __name__ == '__main__':
 #    
 
     
-        
\ No newline at end of file
+        
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py
index 2004e5f..e8f6da4 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/MixedL21Norm.py
@@ -43,19 +43,9 @@ class MixedL21Norm(Function):
         '''
         if not isinstance(x, BlockDataContainer):
             raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) 
-                     
-        if self.SymTensor:
-            
-            #TODO fix this case
-            param = [1]*x.shape[0]
-            param[-1] = 2
-            tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])]
-            res = sum(tmp).sqrt().sum()   
-            
-        else:
-                        
-            tmp = [ el**2 for el in x.containers ]
-            res = sum(tmp).sqrt().sum()
+                                         
+        tmp = [ el**2 for el in x.containers ]
+        res = sum(tmp).sqrt().sum()
 
         return res
                             
@@ -67,7 +57,12 @@ class MixedL21Norm(Function):
         ''' This is the Indicator function of ||\cdot||_{2, \infty}
             which is either 0 if ||x||_{2, \infty} or \infty        
         '''
+        
         return 0.0
+        
+        #tmp = [ el**2 for el in x.containers ]
+        #print(sum(tmp).sqrt().as_array().max())
+        #return sum(tmp).sqrt().as_array().max()
     
     def proximal(self, x, tau, out=None):
         
@@ -80,35 +75,24 @@ class MixedL21Norm(Function):
     
     def proximal_conjugate(self, x, tau, out=None): 
 
-        if self.SymTensor:
-            
-            param = [1]*x.shape[0]
-            param[-1] = 2
-            tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])]
-            frac = [x[i]/(sum(tmp).sqrt()).maximum(1.0) for i in range(x.shape[0])]
-            res = BlockDataContainer(*frac) 
-            
-            return res
-        
-        else:
 
-            if out is None:                                        
-                tmp = [ el*el for el in x.containers]
-                res = sum(tmp).sqrt().maximum(1.0) 
-                frac = [el/res for el in x.containers]
-                return  BlockDataContainer(*frac)   
+        if out is None:                                        
+            tmp = [ el*el for el in x.containers]
+            res = sum(tmp).sqrt().maximum(1.0) 
+            frac = [el/res for el in x.containers]
+            return  BlockDataContainer(*frac)   
+        
             
-                
-            #TODO this is slow, why???
+        #TODO this is slow, why???
 #                return x.divide(x.pnorm().maximum(1.0))
-            else:
-                                
-                res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
-                res = res1.sqrt().maximum(1.0)
-                x.divide(res, out=out)
-                
+        else:
+                            
+            res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
+            res = res1.sqrt().maximum(1.0)
+            x.divide(res, out=out)
+            
 #                x.divide(sum([el*el for el in x.containers]).sqrt().maximum(1.0), out=out)
-                #TODO this is slow, why ???
+            #TODO this is slow, why ???
 #                 x.divide(x.pnorm().maximum(1.0), out=out)
                               
 
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py
index 7caeab2..1db223b 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ScaledFunction.py
@@ -59,7 +59,8 @@ class ScaledFunction(object):
         if out is None:            

             return self.scalar * self.function.gradient(x)    

         else:

-            out.fill( self.scalar * self.function.gradient(x) )

+            self.function.gradient(x, out=out)

+            out *= self.scalar

 

     def proximal(self, x, tau, out=None):

         '''This returns the proximal operator for the function at x, tau

diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py
index cce519a..a019815 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/functions/ZeroFunction.py
@@ -39,14 +39,8 @@ class ZeroFunction(Function):
         indicator function for the set = {0}
         So 0 if x=0, or inf if x neq 0                
         '''
+        return x.maximum(0).sum()
         
-        if x.shape[0]==1:
-            return x.maximum(0).sum()
-        else:
-            if isinstance(x, BlockDataContainer):
-                return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum()
-            else:
-                return x.maximum(0).sum() + x.maximum(0).sum()
     
     def proximal(self, x, tau, out=None):
         if out is None:
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py
index 1d77510..c8bd546 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/BlockOperator.py
@@ -266,15 +266,30 @@ class BlockOperator(Operator):
             # column BlockOperator
             return self.get_item(0,0).domain_geometry()
         else:
-            shape = (self.shape[0], 1)
-            return BlockGeometry(*[el.domain_geometry() for el in self.operators],
-                    shape=shape)
+            # get the geometries column wise
+            # we need only the geometries from the first row
+            # since it is compatible from __init__
+            tmp = []
+            for i in range(self.shape[1]):
+                tmp.append(self.get_item(0,i).domain_geometry())
+            return BlockGeometry(*tmp)                
+                                    
+            #shape = (self.shape[0], 1)
+            #return BlockGeometry(*[el.domain_geometry() for el in self.operators],
+            #        shape=self.shape)
 
     def range_geometry(self):
         '''returns the range of the BlockOperator'''
-        shape = (self.shape[1], 1)
-        return BlockGeometry(*[el.range_geometry() for el in self.operators],
-                    shape=shape)
+        
+        tmp = []
+        for i in range(self.shape[0]):
+            tmp.append(self.get_item(i,0).range_geometry())
+        return BlockGeometry(*tmp)            
+        
+        
+        #shape = (self.shape[1], 1)
+        #return BlockGeometry(*[el.range_geometry() for el in self.operators],
+        #            shape=shape)
         
     def sum_abs_row(self):
         
@@ -312,7 +327,8 @@ class BlockOperator(Operator):
 if __name__ == '__main__':
     
     from ccpi.framework import ImageGeometry
-    from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
+    from ccpi.optimisation.operators import Gradient, Identity, \
+                            SparseFiniteDiff, SymmetrizedGradient, ZeroOperator
 
         
     M, N = 4, 3
@@ -363,4 +379,39 @@ if __name__ == '__main__':
     
     
     
+    ###########################################################################
+    # Block Operator for TGV reconstruction
+    
+    M, N = 2,3
+    ig = ImageGeometry(M, N)
+    ag = ig
+    
+    op11 = Gradient(ig)
+    op12 = Identity(op11.range_geometry())
+    
+    op22 = SymmetrizedGradient(op11.domain_geometry())
+    
+    op21 = ZeroOperator(ig, op22.range_geometry())
+    
+    
+    op31 = Identity(ig, ag)
+    op32 = ZeroOperator(op22.domain_geometry(), ag)
+    
+    operator = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) )    
+    
+    z1 = operator.domain_geometry()
+    z2 = operator.range_geometry()
+    
+    print(z1.shape)
+    print(z2.shape)
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
     
\ No newline at end of file
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 835f96d..f459334 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -7,7 +7,6 @@ Created on Fri Mar  1 22:51:17 2019
 """
 
 from ccpi.optimisation.operators import LinearOperator
-from ccpi.optimisation.ops import PowerMethodNonsquare
 from ccpi.framework import ImageData, BlockDataContainer
 import numpy as np
 
@@ -42,7 +41,7 @@ class FiniteDiff(LinearOperator):
         
         #self.voxel_size = kwargs.get('voxel_size',1)
         # this wrongly assumes a homogeneous voxel size
-        self.voxel_size = self.gm_domain.voxel_size_x
+#        self.voxel_size = self.gm_domain.voxel_size_x
 
 
     def direct(self, x, out=None):
@@ -318,7 +317,7 @@ class FiniteDiff(LinearOperator):
        
     def norm(self):
         x0 = self.gm_domain.allocate('random_int')
-        self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0)
+        self.s1, sall, svec = LinearOperator.PowerMethod(self, 25, x0)
         return self.s1
     
     
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
index 4935435..e0b8a32 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/GradientOperator.py
@@ -7,7 +7,6 @@ Created on Fri Mar  1 22:50:04 2019
 """
 
 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator
-from ccpi.optimisation.ops import PowerMethodNonsquare
 from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
 from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
@@ -90,7 +89,7 @@ class Gradient(LinearOperator):
     def norm(self):
 
         x0 = self.gm_domain.allocate('random')
-        self.s1, sall, svec = PowerMethodNonsquare(self, 10, x0)
+        self.s1, sall, svec = LinearOperator.PowerMethod(self, 10, x0)
         return self.s1
     
     def __rmul__(self, scalar):
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py
index e19304f..f304cc6 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/LinearOperator.py
@@ -6,6 +6,8 @@ Created on Tue Mar  5 15:57:52 2019
 """

 

 from ccpi.optimisation.operators import Operator

+from ccpi.framework import ImageGeometry

+import numpy

 

 

 class LinearOperator(Operator):

@@ -20,3 +22,46 @@ class LinearOperator(Operator):
         

         only available to linear operators'''

         raise NotImplementedError

+    

+    @staticmethod

+    def PowerMethod(operator, iterations, x_init=None):

+        '''Power method to calculate iteratively the Lipschitz constant'''

+        

+        # Initialise random

+        if x_init is None:

+            x0 = operator.domain_geometry().allocate(ImageGeometry.RANDOM_INT)

+        else:

+            x0 = x_init.copy()

+            

+        x1 = operator.domain_geometry().allocate()

+        y_tmp = operator.range_geometry().allocate()

+        s = numpy.zeros(iterations)

+        # Loop

+        for it in numpy.arange(iterations):

+            operator.direct(x0,out=y_tmp)

+            operator.adjoint(y_tmp,out=x1)

+            x1norm = x1.norm()

+            s[it] = x1.dot(x0) / x0.squared_norm()

+            x1.multiply((1.0/x1norm), out=x0)

+        return numpy.sqrt(s[-1]), numpy.sqrt(s), x0

+

+    @staticmethod

+    def PowerMethodNonsquare(op,numiters , x_init=None):

+        '''Power method to calculate iteratively the Lipschitz constant'''

+        

+        if x_init is None:

+            x0 = op.allocate_direct()

+            x0.fill(numpy.random.randn(*x0.shape))

+        else:

+            x0 = x_init.copy()

+        

+        s = numpy.zeros(numiters)

+        # Loop

+        for it in numpy.arange(numiters):

+            x1 = op.adjoint(op.direct(x0))

+            x1norm = x1.norm()

+            s[it] = (x1*x0).sum() / (x0.squared_norm())

+            x0 = (1.0/x1norm)*x1

+        return numpy.sqrt(s[-1]), numpy.sqrt(s), x0

+

+

diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index c38458d..243809b 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -7,7 +7,6 @@ Created on Fri Mar  1 22:53:55 2019
 """
 
 from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator
-from ccpi.optimisation.ops import PowerMethodNonsquare
 from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
 from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
@@ -15,6 +14,12 @@ from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 
 class SymmetrizedGradient(Gradient):
     
+    ''' Symmetrized Gradient, denoted by E: V --> W
+        where V is the Range of the Gradient Operator
+        and W is the Range of the Symmetrized Gradient.                        
+    '''
+    
+    
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
         super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) 
@@ -22,165 +27,218 @@ class SymmetrizedGradient(Gradient):
         '''
          Domain of SymGrad is the Range of Gradient
         '''
+        
         self.gm_domain = self.gm_range 
         self.bnd_cond = bnd_cond
         
         self.channels = self.gm_range.get_item(0).channels
         
-        if self.correlation=='Space':
-            if self.channels>1:
-                pass
-            else: 
-#                # 2D image ---> Dx v1, Dyv2, Dx
-                tmp = self.gm_domain.geometries + (self.gm_domain.get_item(0),)
-                self.gm_range = BlockGeometry(*tmp )
-                self.ind1 = range(self.gm_domain.get_item(0).length)
-                self.ind2 = range(self.gm_domain.get_item(0).length-1, -1, -1)
-#                self.order = myorder = [0,1,2 3]
-                
-        elif self.correlation=='SpaceChannels':
-            if self.channels>1:
-                pass
-            else:
-                raise ValueError('No channels to correlate')        
-                                                     
-        
-    def direct(self, x, out=None):
+        tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries
         
-#        tmp = numpy.zeros(self.gm_range)
-#        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0])
-#        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1])
-#        tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) +
-#                        FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) )
-#        
-#        return type(x)(tmp)
-
-        tmp = [[None]*2]*2
-        for i in range(2):
-            for j in range(2):  
-                tmp[i][j]=FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).adjoint(x.get_item(j))
-        tmp = numpy.array(tmp)
-        z = 0.5 * (tmp.T + tmp)
-        z = z.to
+        self.gm_range = BlockGeometry(*tmp_gm)
         
-        return BlockDataContainer(*z.tolist())
-
-    
-    def adjoint(self, x, out=None):
-        pass
+        self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond)
         
-        res = []
-        for i in range(2):
-            tmp = ImageData(np.zeros(x.get_item(0)))
-            for j in range(2):                
-                tmp += FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).direct(x.get_item(j))
-            res.append(tmp)   
-        return res            
-                
+        if self.gm_domain.shape[0]==2:
+            self.order_ind = [0,2,1,3]
+        else:
+            self.order_ind = [0,3,6,1,4,7,2,5,8]            
                 
         
-#        for 
+    def direct(self, x, out=None):
+        
+        if out is None:
+            
+            tmp = []
+            for i in range(self.gm_domain.shape[0]):
+                for j in range(x.shape[0]):
+                    self.FD.direction = i
+                    tmp.append(self.FD.adjoint(x.get_item(j)))
+                    
+            tmp1 = [tmp[i] for i in self.order_ind]        
+                    
+            res = [0.5 * sum(x) for x in zip(tmp, tmp1)]   
+                    
+            return BlockDataContainer(*res) 
+    
+        else:
+            
+            ind = 0
+            for i in range(self.gm_domain.shape[0]):
+                for j in range(x.shape[0]):
+                    self.FD.direction = i
+                    self.FD.adjoint(x.get_item(j), out=out[ind])
+                    ind+=1                    
+            out1 = BlockDataContainer(*[out[i] for i in self.order_ind])          
+            out.fill( 0.5 * (out + out1) )
+            
+                                               
+    def adjoint(self, x, out=None):
         
-#        tmp = numpy.zeros(self.gm_domain)
-#        
-#        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) +  \
-#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2])
-#                 
-#        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) +  \
-#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1])                 
-#
-#        return type(x)(tmp)          
+        if out is None:
             
-    def alloc_domain_dim(self):
-        return ImageData(numpy.zeros(self.gm_domain))
-    
-    def alloc_range_dim(self):
-        return ImageData(numpy.zeros(self.range_dim))
-    
-    def domain_dim(self):
+            tmp = [None]*self.gm_domain.shape[0]
+            i = 0
+            
+            for k in range(self.gm_domain.shape[0]):
+                tmp1 = 0
+                for j in range(self.gm_domain.shape[0]):
+                    self.FD.direction = j
+                    tmp1 += self.FD.direct(x[i])                    
+                    i+=1
+                tmp[k] = tmp1  
+            return BlockDataContainer(*tmp)
+            
+
+        else:
+            
+            tmp = self.gm_domain.allocate() 
+            i = 0
+            for k in range(self.gm_domain.shape[0]):
+                tmp1 = 0
+                for j in range(self.gm_domain.shape[0]):
+                    self.FD.direction = j
+                    self.FD.direct(x[i], out=tmp[j])
+                    i+=1
+                    tmp1+=tmp[j]
+                out[k].fill(tmp1)
+#            tmp = self.adjoint(x)
+#            out.fill(tmp)
+                    
+                     
+    def domain_geometry(self):
         return self.gm_domain
     
-    def range_dim(self):
+    def range_geometry(self):
         return self.gm_range
                                    
     def norm(self):
-#        return np.sqrt(4*len(self.domainDim()))        
-        #TODO this takes time for big ImageData
-        # for 2D ||grad|| = sqrt(8), 3D ||grad|| = sqrt(12)        
-        x0 = ImageData(np.random.random_sample(self.domain_dim()))
-        self.s1, sall, svec = PowerMethodNonsquare(self, 25, x0)
-        return self.s1 
+
+#        TODO need dot method for BlockDataContainer
+#        return numpy.sqrt(4*self.gm_domain.shape[0])
+    
+#        x0 = self.gm_domain.allocate('random')
+        self.s1, sall, svec = LinearOperator.PowerMethod(self, 50)        
+        return self.s1
     
 
 
 if __name__ == '__main__':   
     
     ###########################################################################  
-    ## Symmetrized Gradient
+    ## Symmetrized Gradient Tests
     from ccpi.framework import DataContainer
     from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
     import numpy as np
     
     N, M = 2, 3
     K = 2
+    C = 2
     
     ig1 = ImageGeometry(N, M)
-    ig2 = ImageGeometry(N, M, channels=K)
+    ig2 = ImageGeometry(N, M, channels=C)
     
     E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann')
-    E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Periodic')
     
-    print(E1.domain_geometry().shape)
-    print(E2.domain_geometry().shape)
+    try:
+        E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann')
+    except:
+        print("No Channels to correlate")
+        
+    E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Neumann')  
+
+    print(E1.domain_geometry().shape, E1.range_geometry().shape)
+    print(E2.domain_geometry().shape, E2.range_geometry().shape)   
+    
+    #check Linear operator property
+    
+    u1 = E1.domain_geometry().allocate('random_int')
+    u2 = E2.domain_geometry().allocate('random_int')
+    
+    # Need to allocate random_int at the Range of SymGradient
+    
+    #a1 = ig1.allocate('random_int')
+    #a2 = ig1.allocate('random_int')
+    #a3 = ig1.allocate('random_int')
+    
+    #a4 = ig1.allocate('random_int')
+    #a5 = ig1.allocate('random_int')    
+    #a6 = ig1.allocate('random_int')
+    
+    # TODO allocate has to create this symmetry by default!!!!!
+    #w1 = BlockDataContainer(*[a1, a2, \
+    #                           a2, a3]) 
+    w1 = E1.range_geometry().allocate('random_int',symmetry=True)
+    
+    LHS = (E1.direct(u1) * w1).sum()
+    RHS = (u1 * E1.adjoint(w1)).sum()
+    
+    numpy.testing.assert_equal(LHS, RHS)     
     
-    u1 = E1.gm_domain.allocate('random_int')
     u2 = E2.gm_domain.allocate('random_int')
-        
     
-    res = E1.direct(u1) 
+    #aa1 = ig2.allocate('random_int')
+    #aa2 = ig2.allocate('random_int')
+    #aa3 = ig2.allocate('random_int')
+    #aa4 = ig2.allocate('random_int')
+    #aa5 = ig2.allocate('random_int')
+    #aa6 = ig2.allocate('random_int')  
     
-    res1 = E1.adjoint(res)
+    #w2 = BlockDataContainer(*[aa1, aa2, aa3, \
+    #                          aa2, aa4, aa5, \
+    #                          aa3, aa5, aa6])     
+    w2 = E2.range_geometry().allocate('random_int',symmetry=True)
+ 
     
-#    Dx = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann')
-#    Dy = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann')
-#    
-#    B = BlockOperator(Dy, Dx)
-#    V = BlockDataContainer(u1,u2)
-#    
-#    res = B.adjoint(V)
+    LHS1 = (E2.direct(u2) * w2).sum()
+    RHS1 = (u2 * E2.adjoint(w2)).sum()
     
-#    ig = (N,M)
-#    ig2 = (2,) + ig
-#    ig3 = (3,) + ig
-#    u1 = ig.allocate('random_int')
-#    w1 = E.gm_range.allocate('random_int')
-#    DataContainer(np.random.randint(10, size=ig3))
+    numpy.testing.assert_equal(LHS1, RHS1)      
     
+    out = E1.range_geometry().allocate()
+    E1.direct(u1, out=out)
+    a1 = E1.direct(u1)
+    numpy.testing.assert_array_equal(a1[0].as_array(), out[0].as_array()) 
+    numpy.testing.assert_array_equal(a1[1].as_array(), out[1].as_array()) 
+    numpy.testing.assert_array_equal(a1[2].as_array(), out[2].as_array()) 
+    numpy.testing.assert_array_equal(a1[3].as_array(), out[3].as_array()) 
     
     
-#    d1 = E.direct(u1)
-#    d2 = E.adjoint(w1)
+    out1 = E1.domain_geometry().allocate()
+    E1.adjoint(w1, out=out1)
+    b1 = E1.adjoint(w1)    
     
-#    LHS = (d1.as_array()[0]*w1.as_array()[0] + \
-#           d1.as_array()[1]*w1.as_array()[1] + \
-#           2*d1.as_array()[2]*w1.as_array()[2]).sum()
-#    
-#    RHS = (u1.as_array()[0]*d2.as_array()[0] + \
-#           u1.as_array()[1]*d2.as_array()[1]).sum()   
-#    
-#    
-#    print(LHS, RHS, E.norm())
-#    
+    LHS_out = (out * w1).sum()
+    RHS_out = (u1 * out1).sum()
+    print(LHS_out, RHS_out)
     
-#    
     
+    out2 = E2.range_geometry().allocate()
+    E2.direct(u2, out=out2)
+    a2 = E2.direct(u2)
+    
+    out21 = E2.domain_geometry().allocate()
+    E2.adjoint(w2, out=out21)
+    b2 = E2.adjoint(w2)    
     
+    LHS_out = (out2 * w2).sum()
+    RHS_out = (u2 * out21).sum()
+    print(LHS_out, RHS_out)    
     
     
+    out = E1.range_geometry().allocate()
+    E1.direct(u1, out=out)
+    E1.adjoint(out, out=out1)
     
+    print(E1.norm())
+    print(E2.norm())
     
     
     
     
     
-    
\ No newline at end of file
+
+#    
+#    
+#    
+#   
\ No newline at end of file
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py
index a7c5f09..8168f0b 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/ZeroOperator.py
@@ -8,32 +8,37 @@ Created on Wed Mar  6 19:25:53 2019
 
 import numpy as np
 from ccpi.framework import ImageData
-from ccpi.optimisation.operators import Operator
+from ccpi.optimisation.operators import LinearOperator
 
-class ZeroOp(Operator):
+class ZeroOperator(LinearOperator):
     
-    def __init__(self, gm_domain, gm_range):
+    def __init__(self, gm_domain, gm_range=None):
+        
+        super(ZeroOperator, self).__init__()             
+
         self.gm_domain = gm_domain
-        self.gm_range = gm_range
-        super(ZeroOp, self).__init__()
+        self.gm_range = gm_range  
+        if self.gm_range is None:
+            self.gm_range = self.gm_domain
+                   
         
     def direct(self,x,out=None):
         if out is None:
-            return ImageData(np.zeros(self.gm_range))
+            return self.gm_range.allocate()
         else:
-            return ImageData(np.zeros(self.gm_range))
+            out.fill(self.gm_range.allocate())
     
     def adjoint(self,x, out=None):
         if out is None:
-            return ImageData(np.zeros(self.gm_domain))
+            return self.gm_domain.allocate()
         else:
-            return ImageData(np.zeros(self.gm_domain))
+            out.fill(self.gm_domain.allocate())
         
     def norm(self):
         return 0
     
-    def domain_dim(self):       
+    def domain_geometry(self):       
         return self.gm_domain  
         
-    def range_dim(self):
+    def range_geometry(self):
         return self.gm_range
\ No newline at end of file
diff --git a/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py
index 7040d3a..23222d4 100644
--- a/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py
+++ b/Wrappers/Python/build/lib/ccpi/optimisation/operators/__init__.py
@@ -18,6 +18,6 @@ from .FiniteDifferenceOperator import FiniteDiff
 from .GradientOperator import Gradient

 from .SymmetrizedGradientOperator import SymmetrizedGradient

 from .IdentityOperator import Identity

-from .ZeroOperator import ZeroOp

-

+from .ZeroOperator import ZeroOperator

+from .LinearOperatorMatrix import LinearOperatorMatrix

 

diff --git a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py
index ef7a9ec..10d15fa 100644
--- a/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py
+++ b/Wrappers/Python/wip/Demos/IMAT_Reconstruction/TV_WhiteBeam_reconstruction.py
@@ -52,11 +52,11 @@ pixv = sinogram_wb.shape[1] # detectors
 igWB = ImageGeometry(voxel_num_x = pixh, voxel_num_y = pixv)
 
 # Load Golden angles
-with open("/media/newhd/vaggelis/CCPi/CCPi-Block/CCPi-Framework/Wrappers/Python/wip/Demos/IMAT_Reconstruction/golden_angles.txt") as f:
+with open("golden_angles.txt") as f:
     angles_string = [line.rstrip() for line in f]
     angles = numpy.array(angles_string).astype(float)
 agWB = AcquisitionGeometry('parallel', '2D',  angles * numpy.pi / 180, pixh)
-op_WB = AstraProjectorSimple(igWB, agWB, 'gpu')
+op_WB = AstraProjectorSimple(igWB, agWB, 'cpu')
 sinogram_aqdata = AcquisitionData(sinogram_wb, agWB)
 
 # BackProjection
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
index 4903c44..58978ae 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
@@ -48,7 +48,7 @@ noisy_data = ImageData(n1)
 # Regularisation Parameter
 alpha = 2
 
-method = '1'
+method = '0'
 
 if method == '0':
 
diff --git a/Wrappers/Python/wip/test_symGrad_method1.py b/Wrappers/Python/wip/test_symGrad_method1.py
deleted file mode 100644
index 36adee1..0000000
--- a/Wrappers/Python/wip/test_symGrad_method1.py
+++ /dev/null
@@ -1,178 +0,0 @@
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 22 14:53:03 2019
-
-@author: evangelos
-"""
-
-from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer
-
-import numpy as np 
-import numpy                          
-import matplotlib.pyplot as plt
-
-from ccpi.optimisation.algorithms import PDHG, PDHG_old
-
-from ccpi.optimisation.operators import BlockOperator, Identity, \
-                        Gradient, SymmetrizedGradient, ZeroOperator
-from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
-                      MixedL21Norm, BlockFunction
-
-from skimage.util import random_noise
-
-from timeit import default_timer as timer
-#def dt(steps):
-#    return steps[-1] - steps[-2]
-
-# Create phantom for TGV Gaussian denoising
-
-N = 3
-
-data = np.zeros((N,N))
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
-
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
-
-data = xv
-data = data/data.max()
-
-plt.imshow(data)
-plt.show()
-
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-ag = ig
-
-# Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.005, seed=10)
-noisy_data = ImageData(n1)
-
-
-plt.imshow(noisy_data.as_array())
-plt.title('Noisy data')
-plt.show()
-
-alpha, beta = 0.2, 1
-
-#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '1'
-
-
-# Create operators
-op11 = Gradient(ig)
-op12 = Identity(op11.range_geometry())
-
-op22 = SymmetrizedGradient(op11.domain_geometry())
-
-op21 = ZeroOperator(ig, op22.range_geometry())
-
-
-op31 = Identity(ig, ag)
-op32 = ZeroOperator(op22.domain_geometry(), ag)
-
-operator1 = BlockOperator(op11, -1*op12, op21, op22, op31, op32, shape=(3,2) ) 
-    
-    
-f1 = alpha * MixedL21Norm()
-f2 = beta * MixedL21Norm() 
-f3 = ZeroFunction()   
-f_B3 = BlockFunction(f1, f2, f3)         
-g_B3 = ZeroFunction()
-    
-    
-    
-# Create operators
-op11 = Gradient(ig)
-op12 = Identity(op11.range_geometry())
-
-op22 = SymmetrizedGradient(op11.domain_geometry())
-
-op21 = ZeroOperator(ig, op22.range_geometry())
-
-operator2 = BlockOperator(op11, -1*op12, \
-                         op21, op22, \
-                         shape=(2,2) )  
-
-#f1 = alpha * MixedL21Norm()
-#f2 = beta * MixedL21Norm()     
-f_B2 = BlockFunction(f1, f2)     
-g_B2 =  0.5 * L2NormSquared(b = noisy_data)    
-     
-
-#%%
-
-x_old1 = operator1.domain_geometry().allocate('random_int')
-y_old1 = operator1.range_geometry().allocate()       
-            
-xbar1 = x_old1.copy()
-x_tmp1 = x_old1.copy()
-x1 = x_old1.copy()
-    
-y_tmp1 = y_old1.copy()
-y1 = y_tmp1.copy()
-
-x_old2 = x_old1.copy()
-y_old2 = operator2.range_geometry().allocate()       
-            
-xbar2 = x_old2.copy()
-x_tmp2 = x_old2.copy()
-x2 = x_old2.copy()
-    
-y_tmp2 = y_old2.copy()
-y2 = y_tmp2.copy()
-
-sigma = 0.4
-tau = 0.4
-
-y_tmp1 = y_old1 +  sigma * operator1.direct(xbar1) 
-y_tmp2 = y_old2 +  sigma * operator2.direct(xbar2)  
-
-numpy.testing.assert_array_equal(y_tmp1[0][0].as_array(), y_tmp2[0][0].as_array()) 
-numpy.testing.assert_array_equal(y_tmp1[0][1].as_array(), y_tmp2[0][1].as_array()) 
-numpy.testing.assert_array_equal(y_tmp1[1][0].as_array(), y_tmp2[1][0].as_array()) 
-numpy.testing.assert_array_equal(y_tmp1[1][1].as_array(), y_tmp2[1][1].as_array()) 
-
-
-y1 = f_B3.proximal_conjugate(y_tmp1, sigma)
-y2 = f_B2.proximal_conjugate(y_tmp2, sigma)
-
-numpy.testing.assert_array_equal(y1[0][0].as_array(), y2[0][0].as_array()) 
-numpy.testing.assert_array_equal(y1[0][1].as_array(), y2[0][1].as_array()) 
-numpy.testing.assert_array_equal(y1[1][0].as_array(), y2[1][0].as_array()) 
-numpy.testing.assert_array_equal(y1[1][1].as_array(), y2[1][1].as_array()) 
-
-
-x_tmp1 = x_old1 - tau * operator1.adjoint(y1)
-x_tmp2 = x_old2 - tau * operator2.adjoint(y2)
-
-numpy.testing.assert_array_equal(x_tmp1[0].as_array(), x_tmp2[0].as_array()) 
-
-
-
-                      
-
-
-
-
-
-
-
-##############################################################################
-#x_1 = operator1.domain_geometry().allocate('random_int')
-#
-#x_2 = BlockDataContainer(x_1[0], x_1[1])
-#
-#res1 = operator1.direct(x_1)
-#res2 = operator2.direct(x_2)
-#
-#print(res1[0][0].as_array(), res2[0][0].as_array())
-#print(res1[0][1].as_array(), res2[0][1].as_array())
-#
-#print(res1[1][0].as_array(), res2[1][0].as_array())
-#print(res1[1][1].as_array(), res2[1][1].as_array())
-#
-##res1 = op11.direct(x1[0]) - op12.direct(x1[1])
-##res2 = op21.direct(x1[0]) - op22.direct(x1[1])