Operator composition (#493)

* CompositionOperator and some refactoring

* Added SumOperator and CompositionOperator

added domain_geometry and optional range_geometry as parameter of Operator. These are saved in _domain_geometry and _range_geometry.
Updated all operator and tests to take notice of this change.

* fighting with Composition

* fix direct and adjoint for CompositeOperator

* fix composition Operator

* remove target OUTPUT to trigger build

* added unit tests

* add test for ZeroOperator

* fixes

* add numba

* removed hard coded path

* removed comments

* fix direct/adjoint with out parameter

* removed __rmul__

* use calculate_norm inherited from LinearOperator

* removed hard coded version tag

* removed cached string version

* use add_custom_target instead add_custom_command

14 files changed, 932 insertions, 186 deletions

diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
index 79bc912..47c6406 100644
--- a/Wrappers/Python/CMakeLists.txt
+++ b/Wrappers/Python/CMakeLists.txt
@@ -3,7 +3,7 @@ option (BUILD_PYTHON_WRAPPER "Build Python Wrapper" ON)
 if (BUILD_PYTHON_WRAPPER)
   find_package(PythonInterp REQUIRED)
   
-  set(PYTHON_DEST_DIR "" CACHE PATH "Directory of the Python wrappers")
+  #set(PYTHON_DEST_DIR "" CACHE PATH "Directory of the Python wrappers")
   if (PYTHON_DEST_DIR)
     set(PYTHON_DEST "${PYTHON_DEST_DIR}")
   else() 
@@ -34,60 +34,68 @@ if (BUILD_PYTHON_WRAPPER)
   if (PYTHONINTERP_FOUND)
       message("Python found " ${PYTHON_EXECUTABLE})
       set(SETUP_PY_IN "${CMAKE_CURRENT_SOURCE_DIR}/setup-cil.py.in")
-      set(SETUP_PY    "${CMAKE_CURRENT_BINARY_DIR}/setup.py")
+      set(SETUP_PY    "${CMAKE_CURRENT_SOURCE_DIR}/setup.py")
       #set(DEPS        "${CMAKE_CURRENT_SOURCE_DIR}/module/__init__.py")
-      set (DEPS       "${CMAKE_BINARY_DIR}/")
-      set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/build/timestamp")
-
-      #configure_file(${SETUP_PY_IN} ${SETUP_PY})
-
+      # set (DEPS       "${CMAKE_CURRENT_SOURCE_DIR}/ccpi")
+      set(OUTPUT      "${CMAKE_CURRENT_BINARY_DIR}/timestamp")
+      file(GLOB_RECURSE DEPS ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/*.py )
+      
+      # add to add_custom_command DEPENDS the list of python files of the project.
+      # as a hack I remove ${OUTPUT}. This should trigger the new build.
+      file( REMOVE ${OUTPUT} )
+      #file( REMOVE_RECURSE ${CMAKE_CURRENT_BINARY_DIR}/build/ )
+      #message(STATUS "We should have removed the build directory now")
+      #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+      #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/test DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+      #file( COPY ${CMAKE_CURRENT_SOURCE_DIR}/setup.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR} )
+      
       
       if (CONDA_BUILD)
-        add_custom_command(OUTPUT ${OUTPUT}
-                        #COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
-                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+        add_custom_target(pythonsetup ALL
                         COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
                                                         ${PYTHON_EXECUTABLE} ${SETUP_PY} -vv install
-                                                        #echo "EDO"
-                        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+      
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
                         COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
                         DEPENDS cilacc)
 
       else()
         if (WIN32)
-          add_custom_command(OUTPUT ${OUTPUT}
-                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
+          add_custom_target(pythonsetup ALL
                         COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
                                                         PREFIX=${CMAKE_SOURCE_DIR}/src/
                                                         LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include
                                                         LIBRARY_LIB=${CMAKE_BINARY_DIR}/${CMAKE_BUILD_TYPE}
                                                         ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
                         COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
                         DEPENDS cilacc)
         else()
-          add_custom_command(OUTPUT ${OUTPUT}
-                        COMMAND ${CMAKE_COMMAND} -E copy_directory ${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR}
-                        COMMAND ${CMAKE_COMMAND} -E env PREFIX=${CMAKE_SOURCE_DIR}/src/
+          add_custom_target(pythonsetup ALL  
+                        COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                        PREFIX=${CMAKE_SOURCE_DIR}/src/
                                                         LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include
                                                         LIBRARY_LIB=${CMAKE_BINARY_DIR}/
-                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py --build-lib ${CMAKE_CURRENT_BINARY_DIR}/build/lib/
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} build_py --verbose --build-lib=${CMAKE_CURRENT_BINARY_DIR}/build/lib
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+                        
                         COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
-                        DEPENDS cilacc
+                        DEPENDS cilacc 
                         )
         endif()
         #set (PYTHON_DEST ${CMAKE_INSTALL_PREFIX}/python/)
 
         install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/build/lib/ccpi
-
-                DESTINATION ${PYTHON_DEST})
+                DESTINATION ${PYTHON_DEST} )
         install(DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/data/ DESTINATION ${CMAKE_INSTALL_PREFIX}/share/ccpi)
         #file(TOUCH ${PYTHON_DEST}/edo/__init__.py)
         
       endif()
       
       
-      add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
-      
+      #add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
+      add_custom_target(PythonWrapper ALL DEPENDS pythonsetup)
+       
   endif()
 
 endif()
diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index a0d139b..310132a 100644
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -261,7 +261,7 @@ class BlockDataContainer(object):
                     kw['out'] = out.get_item(i)
                     if operation == BlockDataContainer.AXPBY:
                         kw['y'] = ot
-                        el.axpby(kw['a'], kw['b'], kw['y'], kw['out'], kw['dtype'], kw['num_threads'])
+                        el.axpby(kw['a'], kw['b'], kw['y'], kw['out'], dtype=kw['dtype'], num_threads=kw['num_threads'])
                     else:
                         op(ot, *args, **kw)
                 else:
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index 39bd983..8fb7ab7 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -26,6 +26,66 @@ from ccpi.optimisation.functions import Function
 
 import functools
 import scipy.special
+try:
+    import numba
+    from numba import jit, prange
+    has_numba = True
+    '''Some parallelisation of KL calls'''
+    @jit(nopython=True)
+    def kl_proximal(x,b, bnoise, tau, out, eta):
+            for i in prange(x.size):
+                out.flat[i] = 0.5 *  ( 
+                    ( x.flat[i] + eta.flat[i] - bnoise.flat[i] - tau ) +\
+                    numpy.sqrt( (x.flat[i] + eta.flat[i] + bnoise.flat[i] - tau)**2. + \
+                        (4. * tau * b.flat[i]) 
+                    )
+                )
+    @jit(nopython=True)
+    def kl_proximal_conjugate(x, b, bnoise, tau, out):
+        #z = x + tau * self.bnoise
+        #return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
+
+        for i in prange(x.size):
+            z = x.flat[i] + ( tau * bnoise.flat[i] )
+            out.flat[i] = 0.5 * ( 
+                (z + 1) - numpy.sqrt((z-1)*(z-1) + 4 * tau * b.flat[i])
+                )
+    @jit(nopython=True)
+    def kl_gradient(x, b, bnoise, out, eta):
+        for i in prange(x.size):
+            out.flat[i] = 1 - b.flat[i]/(x.flat[i] + bnoise.flat[i] + eta.flat[i])
+
+    @jit(nopython=True)
+    def kl_div(x, y, eta):
+        accumulator = 0.
+        for i in prange(x.size):
+            X = x.flat[i]
+            Y = y.flat[i] + eta.flat[i]
+            if X > 0 and Y > 0:
+                # out.flat[i] = X * numpy.log(X/Y) - X + Y
+                accumulator += X * numpy.log(X/Y) - X + Y
+            elif X == 0 and Y >= 0:
+                # out.flat[i] = Y
+                accumulator += Y
+            else:
+                # out.flat[i] = numpy.inf
+                return numpy.inf
+        return accumulator
+    
+    # force a jit
+    x = numpy.asarray(numpy.random.random((10,10)), dtype=numpy.float32)
+    b = numpy.asarray(numpy.random.random((10,10)), dtype=numpy.float32)
+    bnoise = numpy.zeros_like(x)
+    eta = numpy.zeros_like(x)
+    out = numpy.empty_like(x)
+    tau = 1.
+    kl_div(b,x,eta)
+    kl_gradient(x,b,bnoise,out, eta)
+    kl_proximal(x,b, bnoise, tau, out, eta)
+    kl_proximal_conjugate(x,b, bnoise, tau, out)
+    
+except ImportError as ie:
+    has_numba = False
 
 class KullbackLeibler(Function):
     
@@ -84,11 +144,14 @@ class KullbackLeibler(Function):
         r"""Returns the value of the KullbackLeibler function at :math:`(b, x + \eta)`.
         To avoid infinity values, we consider only pixels/voxels for :math:`x+\eta\geq0`.
         """
-                        
-        tmp_sum = (x + self.eta).as_array()
-        ind = tmp_sum >= 0
-        tmp = scipy.special.kl_div(self.b.as_array()[ind], tmp_sum[ind])             
-        return numpy.sum(tmp)         
+        if has_numba:
+            # tmp = numpy.empty_like(x.as_array())
+            return kl_div(self.b.as_array(), x.as_array(), self.eta.as_array())
+        else: 
+            tmp_sum = (x + self.eta).as_array()
+            ind = tmp_sum >= 0
+            tmp = scipy.special.kl_div(self.b.as_array()[ind], tmp_sum[ind])             
+            return numpy.sum(tmp)         
         
     def log(self, datacontainer):
         '''calculates the in-place log of the datacontainer'''
@@ -106,15 +169,26 @@ class KullbackLeibler(Function):
         We require the :math:`x+\eta>0` otherwise we have inf values.
         
         """     
-                                   
-        tmp_sum_array = (x + self.eta).as_array()
-        if out is None:   
-            tmp_out = x.geometry.allocate() 
-            tmp_out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0]         
-            return tmp_out        
-        else:                 
-            x.add(self.eta, out=out)
-            out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0]                                
+        if has_numba:
+            if out is None:
+                out = (x * 0.)
+                out_np = out.as_array()
+                kl_gradient(x.as_array(), self.b.as_array(), self.bnoise.as_array(), out_np, self.eta.as_array())
+                # out.fill(out_np)
+                return out
+            else:
+                out_np = out.as_array()
+                kl_gradient(x.as_array(), self.b.as_array(), self.bnoise.as_array(), out_np, self.eta.as_array())
+                # out.fill(out_np)
+        else:                           
+            tmp_sum_array = (x + self.eta).as_array()
+            if out is None:   
+                tmp_out = x.geometry.allocate() 
+                tmp_out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0]         
+                return tmp_out        
+            else:                 
+                x.add(self.eta, out=out)
+                out.as_array()[tmp_sum_array>0] = 1 - self.b.as_array()[tmp_sum_array>0]/tmp_sum_array[tmp_sum_array>0]                                
 
             
     def convex_conjugate(self, x):
@@ -142,18 +216,30 @@ class KullbackLeibler(Function):
         where :math:`z = x + \tau \eta`
                     
         """
-          
-        if out is None:        
-            return 0.5 *( (x - self.eta - tau) + ( (x + self.eta - tau)**2 + 4*tau*self.b   ) .sqrt() )        
-        else:                      
-            x.add(self.eta, out=out)
-            out -= tau
-            out *= out
-            out.add(self.b * (4 * tau), out=out)
-            out.sqrt(out=out)  
-            out.subtract(tau, out = out)
-            out.add(x, out=out)         
-            out *= 0.5            
+        if has_numba:
+            if out is None:
+                out = (x * 0.)
+                # out_np = numpy.empty_like(out.as_array(), dtype=numpy.float64)
+                out_np = out.as_array()
+                kl_proximal(x.as_array(), self.b.as_array(), self.bnoise.as_array(), tau, out_np, self.eta.as_array())
+                out.fill(out_np)
+                return out
+            else:
+                out_np = out.as_array()
+                kl_proximal(x.as_array(), self.b.as_array(), self.bnoise.as_array(), tau, out_np, self.eta.as_array())
+                # out.fill(out_np)                    
+        else:
+            if out is None:        
+                return 0.5 *( (x - self.eta - tau) + ( (x + self.eta - tau)**2 + 4*tau*self.b   ) .sqrt() )        
+            else:                      
+                x.add(self.eta, out=out)
+                out -= tau
+                out *= out
+                out.add(self.b * (4 * tau), out=out)
+                out.sqrt(out=out)  
+                out.subtract(tau, out = out)
+                out.add(x, out=out)         
+                out *= 0.5            
         
                             
     def proximal_conjugate(self, x, tau, out=None):
@@ -174,7 +260,7 @@ class KullbackLeibler(Function):
             
             self.b.multiply(4*tau, out=out)    
             
-            out.add((tmp)**2, out=out)
+            out.add(tmp.power(2), out=out)
             out.sqrt(out=out)
             out *= -1
             tmp += 2
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 1df7745..f8aef7d 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -56,8 +56,7 @@ class FiniteDiff(LinearOperator):
         :type bnd_cond: str, default :code:`Neumann`
         
         '''
-        super(FiniteDiff, self).__init__() 
-
+        
         self.gm_domain = gm_domain
         self.gm_range = gm_range
         
@@ -75,6 +74,8 @@ 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
+        super(FiniteDiff, self).__init__(domain_geometry=gm_domain, 
+                                         range_geometry=self.gm_range) 
 
 
     def direct(self, x, out=None):
@@ -350,24 +351,24 @@ class FiniteDiff(LinearOperator):
         #else:
         #    out.fill(outa)
             
-    def range_geometry(self):
+    # def range_geometry(self):
         
-        '''
+    #     '''
         
-            Returns the range_geometry of FiniteDiff
+    #         Returns the range_geometry of FiniteDiff
         
-        '''
+    #     '''
         
-        return self.gm_range
+    #     return self.gm_range
     
-    def domain_geometry(self):
+    # def domain_geometry(self):
         
-        '''
+    #     '''
         
-            Returns the domain_geometry of FiniteDiff
+    #         Returns the domain_geometry of FiniteDiff
         
-        '''        
-        return self.gm_domain
+    #     '''        
+    #     return self.gm_domain
 
 
 if __name__ == '__main__':
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index a5feca3..fd51873 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -75,24 +75,29 @@ class Gradient(LinearOperator):
     CORRELATION_SPACE = CORRELATION_SPACE
     CORRELATION_SPACECHANNEL = CORRELATION_SPACECHANNEL
 
-    def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
+    def __init__(self, domain_geometry, bnd_cond = 'Neumann', **kwargs):
         """Constructor method
         """        
-        super(Gradient, self).__init__() 
-
+        
         backend = kwargs.get('backend',C)
 
         correlation = kwargs.get('correlation',CORRELATION_SPACE)
 
-        if correlation == CORRELATION_SPACE and gm_domain.channels > 1:
+        if correlation == CORRELATION_SPACE and domain_geometry.channels > 1:
             #numpy implementation only for now
             backend = NUMPY
             warnings.warn("Warning: correlation='Space' on multi-channel dataset will use `numpy` backend")
 
         if backend == NUMPY:
-            self.operator = Gradient_numpy(gm_domain, bnd_cond=bnd_cond, **kwargs)
+            self.operator = Gradient_numpy(domain_geometry, bnd_cond=bnd_cond, **kwargs)
         else:
-            self.operator = Gradient_C(gm_domain, bnd_cond=bnd_cond, **kwargs)
+            self.operator = Gradient_C(domain_geometry, bnd_cond=bnd_cond, **kwargs)
+        
+        super(Gradient, self).__init__(domain_geometry=domain_geometry, 
+                                       range_geometry=self.operator.range_geometry()) 
+
+        self.gm_range = self.range_geometry()
+        self.gm_domain = self.domain_geometry()
 
 
     def direct(self, x, out=None):
@@ -120,16 +125,16 @@ class Gradient(LinearOperator):
         """            
         return self.operator.adjoint(x, out=out)
 
-    def domain_geometry(self):
-        '''Returns domain_geometry of Gradient'''
+    # def domain_geometry(self):
+    #     '''Returns domain_geometry of Gradient'''
         
-        return self.operator.gm_domain
+    #     return self.operator.gm_domain
     
-    def range_geometry(self):
+    # def range_geometry(self):
         
-        '''Returns range_geometry of Gradient'''
+    #     '''Returns range_geometry of Gradient'''
         
-        return self.operator.gm_range
+    #     return self.operator.gm_range
 
 class Gradient_numpy(LinearOperator):
     
@@ -143,7 +148,6 @@ class Gradient_numpy(LinearOperator):
         :param correlation: optional, :code:`SpaceChannel` or :code:`Space`
         :type correlation: str, optional, default :code:`Space`
         '''
-        super(Gradient_numpy, self).__init__() 
                 
         self.gm_domain = gm_domain # Domain of Grad Operator
         
@@ -190,6 +194,10 @@ class Gradient_numpy(LinearOperator):
         
         # Call FiniteDiff operator        
         self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond)
+
+        super(Gradient_numpy, self).__init__(domain_geometry=self.gm_domain, 
+                                             range_geometry=self.gm_range) 
+        
         print("Initialised GradientOperator with numpy backend")               
         
     def direct(self, x, out=None):
@@ -240,14 +248,6 @@ class Gradient_numpy(LinearOperator):
         
         return self.gm_range
     
-    def __rmul__(self, scalar):
-        
-        '''Multiplication of Gradient with a scalar        
-            
-            Returns: ScaledOperator
-        '''        
-        
-        return ScaledOperator(self, scalar) 
     
     ###########################################################################
     ###############  For preconditioning ######################################
@@ -348,8 +348,6 @@ class Gradient_C(LinearOperator):
 
     def __init__(self, gm_domain, gm_range=None, bnd_cond = NEUMANN, **kwargs):
 
-        super(Gradient_C, self).__init__() 
-
         self.num_threads = kwargs.get('num_threads',NUM_THREADS)
 
         self.gm_domain = gm_domain
@@ -375,6 +373,9 @@ class Gradient_C(LinearOperator):
         else:
             raise ValueError('Number of dimensions not supported, expected 2, 3 or 4, got {}'.format(len(gm_domain.shape)))
       #self.num_threads
+        # super(Gradient_C, self).__init__() 
+        super(Gradient_C, self).__init__(domain_geometry=self.gm_domain, 
+                                             range_geometry=self.gm_range) 
         print("Initialised GradientOperator with C backend running with ", cilacc.openMPtest(self.num_threads)," threads")               
 
     @staticmethod 
@@ -418,17 +419,17 @@ class Gradient_C(LinearOperator):
         if return_val is True:
             return out
 
-    def domain_geometry(self):
+    # def domain_geometry(self):
         
-        '''Returns domain_geometry of Gradient'''
+    #     '''Returns domain_geometry of Gradient'''
         
-        return self.gm_domain
+    #     return self.gm_domain
     
-    def range_geometry(self):
+    # def range_geometry(self):
         
-        '''Returns range_geometry of Gradient'''
+    #     '''Returns range_geometry of Gradient'''
         
-        return self.gm_range
+    #     return self.gm_range
 
        
 if __name__ == '__main__':
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index a4c6ae3..2c14c5d 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -35,14 +35,14 @@ class Identity(LinearOperator):
     '''    
     
     
-    def __init__(self, gm_domain, gm_range=None):
+    def __init__(self, domain_geometry, range_geometry=None):
 
-        self.gm_domain = gm_domain
-        self.gm_range = gm_range  
-        if self.gm_range is None:
-            self.gm_range = self.gm_domain
         
-        super(Identity, self).__init__()
+        if range_geometry is None:
+            range_geometry = domain_geometry
+        
+        super(Identity, self).__init__(domain_geometry=domain_geometry,
+                                       range_geometry=range_geometry)
         
     def direct(self,x,out=None):
         
@@ -68,18 +68,6 @@ class Identity(LinearOperator):
         '''Evaluates operator norm of Identity'''        
         
         return 1.0
-        
-    def domain_geometry(self): 
-        
-        '''Returns domain_geometry of Identity'''
-        
-        return self.gm_domain
-        
-    def range_geometry(self):
-        
-        '''Returns range_geometry of Identity'''         
-        
-        return self.gm_range
     
     
     ###########################################################################
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index 6f13532..95bcc92 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -39,15 +39,16 @@ class LinearOperatorMatrix(LinearOperator):
         '''
         self.A = A
         M_A, N_A = self.A.shape
-        self.gm_domain = VectorGeometry(N_A)
-        self.gm_range = VectorGeometry(M_A)
+        domain_geometry = VectorGeometry(N_A)
+        range_geometry = VectorGeometry(M_A)
         self.s1 = None   # Largest singular value, initially unknown
-        super(LinearOperatorMatrix, self).__init__()
+        super(LinearOperatorMatrix, self).__init__(domain_geometry=domain_geometry,
+                                                   range_geometry=range_geometry)
         
     def direct(self,x, out=None):        
                 
         if out is None:
-            tmp = self.gm_range.allocate()
+            tmp = self.range_geometry().allocate()
             tmp.fill(numpy.dot(self.A,x.as_array()))
             return tmp
         else:
@@ -58,7 +59,7 @@ class LinearOperatorMatrix(LinearOperator):
 
     def adjoint(self,x, out=None):
         if out is None:
-            tmp = self.gm_domain.allocate()
+            tmp = self.domain_geometry().allocate()
             tmp.fill(numpy.dot(self.A.transpose(),x.as_array()))
             return tmp
         else:            
@@ -70,8 +71,4 @@ class LinearOperatorMatrix(LinearOperator):
     def calculate_norm(self, **kwargs):
         # If unknown, compute and store. If known, simply return it.
         return svds(self.A,1,return_singular_vectors=False)[0]
-        
-    def domain_geometry(self):
-        return self.gm_domain
-    def range_geometry(self):
-        return self.gm_range
+    
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 14d53e8..d49bc1a 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -18,13 +18,23 @@
 from __future__ import absolute_import

 from __future__ import division

 from __future__ import print_function

-

-from ccpi.optimisation.operators.ScaledOperator import ScaledOperator

+from numbers import Number

+import numpy

+import functools

 

 class Operator(object):

     '''Operator that maps from a space X -> Y'''

-    def __init__(self, **kwargs):

-        self.__norm = None

+    def __init__(self, domain_geometry, **kwargs):

+        r'''

+        Creator

+

+        :param domain_geometry: domain of the operator

+        :param range_geometry: range of the operator

+        :type range_geometry: optional, default None

+        '''

+        self._norm = None

+        self._domain_geometry = domain_geometry

+        self._range_geometry = kwargs.get('range_geometry', None)

 

     def is_linear(self):

         '''Returns if the operator is linear'''

@@ -51,20 +61,425 @@ class Operator(object):
         :parameter force: forces the recalculation of the norm

         :type force: boolean, default :code:`False`

         '''

-        if self.__norm is None or kwargs.get('force', False):

-            self.__norm = self.calculate_norm(**kwargs)

-        return self.__norm

+        if self._norm is None or kwargs.get('force', False):

+            self._norm = self.calculate_norm(**kwargs)

+        return self._norm

     def calculate_norm(self, **kwargs):

         '''Calculates the norm of the Operator'''

         raise NotImplementedError

     def range_geometry(self):

         '''Returns the range of the Operator: Y space'''

-        raise NotImplementedError

+        return self._range_geometry

     def domain_geometry(self):

         '''Returns the domain of the Operator: X space'''

-        raise NotImplementedError

+        return self._domain_geometry

     def __rmul__(self, scalar):

         '''Defines the multiplication by a scalar on the left

 

         returns a ScaledOperator'''

         return ScaledOperator(self, scalar)

+    

+    def compose(self, *other, **kwargs):

+        # TODO: check equality of domain and range of operators        

+        #if self.operator2.range_geometry != self.operator1.domain_geometry:

+        #    raise ValueError('Cannot compose operators, check domain geometry of {} and range geometry of {}'.format(self.operato1,self.operator2))    

+        

+        return CompositionOperator(self, *other, **kwargs) 

+

+    def __add__(self, other):

+        return SumOperator(self, other)

+

+    def __mul__(self, scalar):

+        return self.__rmul__(scalar)    

+    

+    def __neg__(self):

+        """ Return -self """

+        return -1 * self    

+        

+    def __sub__(self, other):

+        """ Returns the subtraction of the operators."""

+        return self + (-1) * other   

+

+

+class LinearOperator(Operator):

+    '''A Linear Operator that maps from a space X <-> Y'''

+    def __init__(self, domain_geometry, **kwargs):

+        super(LinearOperator, self).__init__(domain_geometry, **kwargs)

+    def is_linear(self):

+        '''Returns if the operator is linear'''

+        return True

+    def adjoint(self,x, out=None):

+        '''returns the adjoint/inverse operation

+        

+        only available to linear operators'''

+        raise NotImplementedError

+    

+    @staticmethod

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

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

+        

+        :param operator: input operator

+        :type operator: :code:`LinearOperator`

+        :param iterations: number of iterations to run

+        :type iteration: int

+        :param x_init: starting point for the iteration in the operator domain

+        :returns: tuple with: L, list of L at each iteration, the data the iteration worked on.

+        '''

+        

+        # Initialise random

+        if x_init is None:

+            x0 = operator.domain_geometry().allocate('random')

+        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()

+            if hasattr(x0, 'squared_norm'):

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

+            else:

+                x0norm = x0.norm()

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

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

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

+

+    def calculate_norm(self, **kwargs):

+        '''Returns the norm of the LinearOperator as calculated by the PowerMethod

+        

+        :param iterations: number of iterations to run

+        :type iterations: int, optional, default = 25

+        :param x_init: starting point for the iteration in the operator domain

+        :type x_init: same type as domain, a subclass of :code:`DataContainer`, optional, None

+        :parameter force: forces the recalculation of the norm

+        :type force: boolean, default :code:`False`

+        '''

+        x0 = kwargs.get('x_init', None)

+        iterations = kwargs.get('iterations', 25)

+        s1, sall, svec = LinearOperator.PowerMethod(self, iterations, x_init=x0)

+        return s1

+

+    @staticmethod

+    def dot_test(operator, domain_init=None, range_init=None, verbose=False, **kwargs):

+        r'''Does a dot linearity test on the operator

+        

+        Evaluates if the following equivalence holds

+        

+        .. math::

+        

+          Ax\times y = y \times A^Tx

+        

+        The equivalence is tested within a user specified precision

+

+        .. code::

+        

+          abs(desired-actual) < 1.5 * 10**(-decimal)

+

+        :param operator: operator to test

+        :param range_init: optional initialisation container in the operator range 

+        :param domain_init: optional initialisation container in the operator domain 

+        :returns: boolean, True if the test is passed.

+        :param decimal: desired precision

+        :type decimal: int, optional, default 4

+        '''

+        if range_init is None:

+            y = operator.range_geometry().allocate('random_int')

+        else:

+            y = range_init

+        if domain_init is None:

+            x = operator.domain_geometry().allocate('random_int')

+        else:

+            x = domain_init

+            

+        fx = operator.direct(x)

+        by = operator.adjoint(y)

+        a = fx.dot(y)

+        b = by.dot(x)

+        if verbose:

+            print ('Left hand side  {}, \nRight hand side {}'.format(a, b))

+            print ("decimal ", kwargs.get('decimal', 4))

+        try:

+            numpy.testing.assert_almost_equal(abs((a-b)/a), 0., decimal=kwargs.get('decimal',4))

+            return True

+        except AssertionError as ae:

+            print (ae)

+            return False

+        

+        

+class ScaledOperator(Operator):

+    

+    

+    '''ScaledOperator

+

+    A class to represent the scalar multiplication of an Operator with a scalar.

+    It holds an operator and a scalar. Basically it returns the multiplication

+    of the result of direct and adjoint of the operator with the scalar.

+    For the rest it behaves like the operator it holds.

+    

+    :param operator: a Operator or LinearOperator

+    :param scalar: a scalar multiplier

+    

+    Example:

+       The scaled operator behaves like the following:

+

+    .. code-block:: python

+

+      sop = ScaledOperator(operator, scalar)

+      sop.direct(x) = scalar * operator.direct(x)

+      sop.adjoint(x) = scalar * operator.adjoint(x)

+      sop.norm() = operator.norm()

+      sop.range_geometry() = operator.range_geometry()

+      sop.domain_geometry() = operator.domain_geometry()

+

+    '''

+    

+    def __init__(self, operator, scalar, **kwargs):

+        '''creator

+

+        :param operator: a Operator or LinearOperator

+        :param scalar: a scalar multiplier

+        :type scalar: float'''

+

+        super(ScaledOperator, self).__init__(domain_geometry=operator.domain_geometry(), 

+                                             range_geometry=operator.range_geometry())

+        if not isinstance (scalar, Number):

+            raise TypeError('expected scalar: got {}'.format(type(scalar)))

+        self.scalar = scalar

+        self.operator = operator

+    def direct(self, x, out=None):

+        '''direct method'''

+        if out is None:

+            return self.scalar * self.operator.direct(x, out=out)

+        else:

+            self.operator.direct(x, out=out)

+            out *= self.scalar

+    def adjoint(self, x, out=None):

+        '''adjoint method'''

+        if self.operator.is_linear():

+            if out is None:

+                return self.scalar * self.operator.adjoint(x, out=out)

+            else:

+                self.operator.adjoint(x, out=out)

+                out *= self.scalar

+        else:

+            raise TypeError('Operator is not linear')

+    def norm(self, **kwargs):

+        '''norm of the operator'''

+        return numpy.abs(self.scalar) * self.operator.norm(**kwargs)

+    # def range_geometry(self):

+    #     '''range of the operator'''

+    #     return self.operator.range_geometry()

+    # def domain_geometry(self):

+    #     '''domain of the operator'''

+    #     return self.operator.domain_geometry()

+    def is_linear(self):

+        '''returns whether the operator is linear

+        

+        :returns: boolean '''

+        return self.operator.is_linear()

+

+

+###############################################################################

+################   SumOperator  ###########################################

+###############################################################################      

+    

+class SumOperator(Operator):

+    

+    

+    def __init__(self, operator1, operator2):

+                

+        self.operator1 = operator1

+        self.operator2 = operator2

+        

+        # if self.operator1.domain_geometry() != self.operator2.domain_geometry():

+        #     raise ValueError('Domain geometry of {} is not equal with domain geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    

+                

+        # if self.operator1.range_geometry() != self.operator2.range_geometry():

+        #     raise ValueError('Range geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    

+            

+        self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear()            

+        

+        super(SumOperator, self).__init__(domain_geometry=self.operator1.domain_geometry(),

+                                          range_geometry=self.operator1.range_geometry()) 

+                                  

+    def direct(self, x, out=None):

+        

+        if out is None:

+            return self.operator1.direct(x) + self.operator2.direct(x)

+        else:

+            #TODO check if allcating with tmp is faster            

+            self.operator1.direct(x, out=out)

+            out.add(self.operator2.direct(x), out=out)     

+

+    def adjoint(self, x, out=None):

+        

+        if self.linear_flag:        

+            if out is None:

+                return self.operator1.adjoint(x) + self.operator2.adjoint(x)

+            else:

+                #TODO check if allcating with tmp is faster            

+                self.operator1.adjoint(x, out=out)

+                out.add(self.operator2.adjoint(x), out=out) 

+        else:

+            raise ValueError('No adjoint operation with non-linear operators')

+                                        

+    def is_linear(self):

+        return self.linear_flag 

+    

+    def calculate_norm(self, **kwargs):

+        if self.is_linear():

+            return LinearOperator.calculate_norm(self, **kwargs)

+

+###############################################################################

+################   Composition  ###########################################

+###############################################################################             

+    

+class Composition2Operator(Operator):

+    

+    def __init__(self, operator1, operator2):

+        

+        self.operator1 = operator1

+        self.operator2 = operator2

+        

+        self.linear_flag = self.operator1.is_linear() and self.operator2.is_linear()        

+        

+        if self.operator2.range_geometry() != self.operator1.domain_geometry():

+            raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    

+                

+        super(Composition2Operator, self).__init__(self.operator1.domain_geometry(),

+                                          self.operator2.range_geometry()) 

+        

+    def direct(self, x, out = None):

+

+        if out is None:

+            return self.operator1.direct(self.operator2.direct(x))

+        else:

+            tmp = self.operator2.range_geometry().allocate()

+            self.operator2.direct(x, out = tmp)

+            self.operator1.direct(tmp, out = out)

+            

+    def adjoint(self, x, out = None):

+        

+        if self.linear_flag: 

+            

+            if out is None:

+                return self.operator2.adjoint(self.operator1.adjoint(x))

+            else:

+                

+                tmp = self.operator1.domain_geometry().allocate()

+                self.operator1.adjoint(x, out=tmp)

+                self.operator2.adjoint(tmp, out=out)

+        else:

+            raise ValueError('No adjoint operation with non-linear operators')

+            

+

+    def is_linear(self):

+        return self.linear_flag             

+            

+    def calculate_norm(self, **kwargs):

+        if self.is_linear():

+            return LinearOperator.calculate_norm(self, **kwargs)

+

+class CompositionOperator(Operator):

+    

+    def __init__(self, *operators, **kwargs):

+        

+        # get a reference to the operators

+        self.operators = operators

+        

+        self.linear_flag = functools.reduce(lambda x,y: x and y.is_linear(),

+                                            self.operators, True)

+        self.preallocate = kwargs.get('preallocate', False)

+        if self.preallocate:

+            self.tmp_domain = [op.domain_geometry().allocate() for op in self.operators[:-1]]

+            self.tmp_range = [op.range_geometry().allocate() for op in self.operators[1:]]

+            # pass

+        

+        # TODO address the equality of geometries

+        # if self.operator2.range_geometry() != self.operator1.domain_geometry():

+        #     raise ValueError('Domain geometry of {} is not equal with range geometry of {}'.format(self.operator1.__class__.__name__,self.operator2.__class__.__name__))    

+                

+        super(CompositionOperator, self).__init__(

+            domain_geometry=self.operators[-1].domain_geometry(),

+            range_geometry=self.operators[0].range_geometry()) 

+        

+    def direct(self, x, out = None):

+

+        if out is None:

+            #return self.operator1.direct(self.operator2.direct(x))

+            # return functools.reduce(lambda X,operator: operator.direct(X), 

+            #                        self.operators[::-1][1:],

+            #                        self.operators[-1].direct(x))

+            for i,operator in enumerate(self.operators[::-1]):

+                if i == 0:

+                    step = operator.direct(x)

+                else:

+                    step = operator.direct(step)

+            return step

+

+        else:

+            # tmp = self.operator2.range_geometry().allocate()

+            # self.operator2.direct(x, out = tmp)

+            # self.operator1.direct(tmp, out = out)

+            

+            # out.fill (

+            #     functools.reduce(lambda X,operator: operator.direct(X), 

+            #                        self.operators[::-1][1:],

+            #                        self.operators[-1].direct(x))

+            # )

+            

+            # TODO this is a bit silly but will handle the pre allocation later

+            

+            for i,operator in enumerate(self.operators[::-1]):

+                if i == 0:

+                    operator.direct(x, out=self.tmp_range[i])

+                elif i == len(self.operators) - 1:

+                    operator.direct(self.tmp_range[i-1], out=out)

+                else:

+                    operator.direct(self.tmp_range[i-1], out=self.tmp_range[i])

+            

+            

+    def adjoint(self, x, out = None):

+        

+        if self.linear_flag: 

+            

+            if out is not None:

+                # return self.operator2.adjoint(self.operator1.adjoint(x))

+                # return functools.reduce(lambda X,operator: operator.adjoint(X), 

+                #                    self.operators[1:],

+                #                    self.operators[0].adjoint(x))

+

+                for i,operator in enumerate(self.operators):

+                    if i == 0:

+                        operator.adjoint(x, out=self.tmp_domain[i])

+                    elif i == len(self.operators) - 1:

+                        step = operator.adjoint(self.tmp_domain[i-1], out=out)

+                    else:

+                        operator.adjoint(self.tmp_domain[i-1], out=self.tmp_domain[i])

+                return

+

+            else:

+                for i,operator in enumerate(self.operators):

+                    if i == 0:

+                        step = operator.adjoint(x)

+                    else:

+                        step = operator.adjoint(step)

+                

+                return step

+        else:

+            raise ValueError('No adjoint operation with non-linear operators')

+            

+

+    def is_linear(self):

+        return self.linear_flag             

+            

+    def calculate_norm(self, **kwargs):

+        if self.is_linear():

+            return LinearOperator.calculate_norm(self, **kwargs)

+

+

+

diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 747db65..310d4e8 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -29,16 +29,16 @@ class SparseFiniteDiff(object):
     '''Create Sparse Matrices for the Finite Difference Operator'''
     
     
-    def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
+    def __init__(self, domain_geometry, range_geometry=None, 
+      direction=0, bnd_cond = 'Neumann'):
         
-        super(SparseFiniteDiff, self).__init__() 
-        self.gm_domain = gm_domain
-        self.gm_range = gm_range
+        super(SparseFiniteDiff, self).__init__(domain_geometry=domain_geometry,
+                                               range_geometry=range_geometry) 
         self.direction = direction
         self.bnd_cond = bnd_cond
         
-        if self.gm_range is None:
-            self.gm_range = self.gm_domain
+        if self.range_geometry is None:
+            self.range_geometry = self.domain_geometry
             
         self.get_dims = [i for i in gm_domain.shape]  
         
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index 6fd2e86..9c2300f 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -48,33 +48,36 @@ class SymmetrizedGradient(LinearOperator):
     CORRELATION_SPACE = "Space"
     CORRELATION_SPACECHANNEL = "SpaceChannels"
     
-    def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
+    def __init__(self, domain_geometry, bnd_cond = 'Neumann', **kwargs):
         '''creator
         
-        :param gm_domain: domain of the operator
+        :param domain_geometry: domain of the operator
         :param bnd_cond: boundary condition, either :code:`Neumann` or :code:`Periodic`.
         :type bnd_cond: str, optional, default :code:`Neumann`
         :param correlation: :code:`SpaceChannel` or :code:`Channel`
         :type correlation: str, optional, default :code:`Channel`
         '''
-        super(SymmetrizedGradient, self).__init__() 
+        # super(SymmetrizedGradient, self).__init__(domain_geometry=domain_geometry) 
                 
-        self.gm_domain = gm_domain
         self.bnd_cond = bnd_cond
         self.correlation = kwargs.get('correlation',SymmetrizedGradient.CORRELATION_SPACE)
                 
-        tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries
+        tmp_gm = len(domain_geometry.geometries)*domain_geometry.geometries
         
-        self.gm_range = BlockGeometry(*tmp_gm)
         
         # Define FD operator. We need one geometry from the BlockGeometry of the domain
-        self.FD = FiniteDiff(self.gm_domain.get_item(0), direction = 0, bnd_cond = self.bnd_cond)
+        self.FD = FiniteDiff(domain_geometry.get_item(0), direction = 0, 
+                             bnd_cond = self.bnd_cond)
         
-        if self.gm_domain.shape[0]==2:
+        if domain_geometry.shape[0]==2:
             self.order_ind = [0,2,1,3]
         else:
             self.order_ind = [0,3,6,1,4,7,2,5,8]            
-                
+        
+        super(SymmetrizedGradient, self).__init__(
+                                          domain_geometry=domain_geometry, 
+                                          range_geometry=BlockGeometry(*tmp_gm))
+        
         
     def direct(self, x, out=None):
         
@@ -83,7 +86,7 @@ class SymmetrizedGradient(LinearOperator):
         if out is None:
             
             tmp = []
-            for i in range(self.gm_domain.shape[0]):
+            for i in range(self.domain_geometry().shape[0]):
                 for j in range(x.shape[0]):
                     self.FD.direction = i
                     tmp.append(self.FD.adjoint(x.get_item(j)))
@@ -97,7 +100,7 @@ class SymmetrizedGradient(LinearOperator):
         else:
             
             ind = 0
-            for i in range(self.gm_domain.shape[0]):
+            for i in range(self.domain_geometry().shape[0]):
                 for j in range(x.shape[0]):
                     self.FD.direction = i
                     self.FD.adjoint(x.get_item(j), out=out[ind])
@@ -110,12 +113,12 @@ class SymmetrizedGradient(LinearOperator):
         
         if out is None:
             
-            tmp = [None]*self.gm_domain.shape[0]
+            tmp = [None]*self.domain_geometry().shape[0]
             i = 0
             
-            for k in range(self.gm_domain.shape[0]):
+            for k in range(self.domain_geometry().shape[0]):
                 tmp1 = 0
-                for j in range(self.gm_domain.shape[0]):
+                for j in range(self.domain_geometry().shape[0]):
                     self.FD.direction = j
                     tmp1 += self.FD.direct(x[i])                    
                     i+=1
@@ -125,23 +128,18 @@ class SymmetrizedGradient(LinearOperator):
 
         else:
             
-            tmp = self.gm_domain.allocate() 
+            tmp = self.domain_geometry().allocate() 
             i = 0
-            for k in range(self.gm_domain.shape[0]):
+            for k in range(self.domain_geometry().shape[0]):
                 tmp1 = 0
-                for j in range(self.gm_domain.shape[0]):
+                for j in range(self.domain_geometry().shape[0]):
                     self.FD.direction = j
                     self.FD.direct(x[i], out=tmp[j])
                     i+=1
                     tmp1+=tmp[j]
                 out[k].fill(tmp1)
                     
-                     
-    def domain_geometry(self):
-        return self.gm_domain
     
-    def range_geometry(self):
-        return self.gm_range
 
 if __name__ == '__main__':   
     
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index 0feafd8..bc4f08e 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -41,14 +41,12 @@ class ZeroOperator(LinearOperator):
                        
      '''
     
-    def __init__(self, gm_domain, gm_range=None):
-        
-        super(ZeroOperator, self).__init__()             
+    def __init__(self, domain_geometry, range_geometry=None):
+        if range_geometry is None:
+            range_geometry = domain_geometry.clone()
+        super(ZeroOperator, self).__init__(domain_geometry=domain_geometry, 
+                                           range_geometry=range_geometry)
 
-        self.gm_domain = gm_domain
-        self.gm_range = gm_range  
-        if self.gm_range is None:
-            self.gm_range = self.gm_domain
                    
         
     def direct(self,x,out=None):
@@ -57,18 +55,18 @@ class ZeroOperator(LinearOperator):
         
         
         if out is None:
-            return self.gm_range.allocate()
+            return self.range_geometry().allocate()
         else:
-            out.fill(self.gm_range.allocate())
+            out.fill(self.range_geometry.allocate())
     
     def adjoint(self,x, out=None):
         
         '''Returns O^{*}(y)'''        
         
         if out is None:
-            return self.gm_domain.allocate()
+            return self.domain_geometry().allocate()
         else:
-            out.fill(self.gm_domain.allocate())
+            out.fill(self.domain_geometry().allocate())
         
     def calculate_norm(self, **kwargs):
         
@@ -76,15 +74,4 @@ class ZeroOperator(LinearOperator):
         
         return 0
     
-    def domain_geometry(self): 
-        
-        '''Returns domain_geometry of ZeroOperator'''
-        
-        
-        return self.gm_domain  
-        
-    def range_geometry(self):
-        
-        '''Returns domain_geometry of ZeroOperator'''
-        
-        return self.gm_range
+    
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
index 15c8932..cc25c92 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
@@ -16,11 +16,12 @@
 #   See the License for the specific language governing permissions and

 #   limitations under the License.

 

-from .Operator import Operator

-from .LinearOperator import LinearOperator

-from .ScaledOperator import ScaledOperator

+from .Operator import Operator, LinearOperator, ScaledOperator, SumOperator,\

+ CompositionOperator, Composition2Operator

+#from .LinearOperator import LinearOperator

+#from .ScaledOperator import ScaledOperator

 from .BlockOperator import BlockOperator

-from .BlockScaledOperator import BlockScaledOperator

+# from .BlockScaledOperator import BlockScaledOperator

 

 from .SparseFiniteDiff import SparseFiniteDiff

 from .ShrinkageOperator import ShrinkageOperator

diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py
index d0f568b..2f33f08 100644
--- a/Wrappers/Python/test/test_BlockOperator.py
+++ b/Wrappers/Python/test/test_BlockOperator.py
@@ -81,8 +81,8 @@ class TestBlockOperator(unittest.TestCase):
                    ImageGeometry(10,22,31) , \
                    ImageGeometry(10,20,31) ]
             x = [ g.allocate() for g in ig ]
-            ops = [ Identity(g, gm_range=r) for g,r in zip(ig, rg0) ]
-            ops += [ Identity(g, gm_range=r) for g,r in zip(ig, rg1) ]
+            ops = [ Identity(g, range_geometry=r) for g,r in zip(ig, rg0) ]
+            ops += [ Identity(g, range_geometry=r) for g,r in zip(ig, rg1) ]
 
             K = BlockOperator(*ops, shape=(2,3))
             print ("K col comp? " , K.column_wise_compatible())
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 57dd41f..8312c6d 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -23,7 +23,12 @@ import numpy
 from timeit import default_timer as timer
 from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
 from ccpi.optimisation.operators import LinearOperator, LinearOperatorMatrix
+import numpy   
+from ccpi.optimisation.operators import SumOperator, Gradient,\
+            ZeroOperator, SymmetrizedGradient, CompositionOperator
 
+from ccpi.framework import TestData
+import os
 
 def dt(steps):
     return steps[-1] - steps[-2]
@@ -221,6 +226,10 @@ class TestOperator(CCPiTestClass):
         for n in [norm, norm2, norm3, norm4, norm5]:
             print ("norm {}", format(n))
 
+    
+
+        
+
 
 class TestGradients(CCPiTestClass): 
     def setUp(self):
@@ -339,6 +348,7 @@ class TestGradients(CCPiTestClass):
 
 
 
+
 class TestBlockOperator(unittest.TestCase):
     def assertBlockDataContainerEqual(self, container1, container2):
         print ("assert Block Data Container Equal")
@@ -645,6 +655,260 @@ class TestBlockOperator(unittest.TestCase):
         u1 = B.adjoint(w)
         self.assertEqual((w * w1).sum() , (u1*u).sum())
 
+class TestOperatorCompositionSum(unittest.TestCase):
+    def setUp(self):
+        
+        self.data = TestData().load(TestData.BOAT, size=(128,128))
+        self.ig = self.data.geometry
+
+    def test_SumOperator(self):
+
+        # numpy.random.seed(1)
+        ig = self.ig
+        data = self.data
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        # z = ZeroOperator(domain_geometry=ig)
+        # sym = SymmetrizedGradient(domain_geometry=ig)
+
+        c = SumOperator(Id1,Id2)
+        out = c.direct(data)
+
+        numpy.testing.assert_array_almost_equal(out.as_array(),3 * data.as_array())
+
+
+    def test_CompositionOperator_direct1(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id2)
+
+        out1 = G.direct(data)
+        out2 = d.direct(data)
+
+
+        numpy.testing.assert_array_almost_equal(out2.get_item(0).as_array(),  out1.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(out2.get_item(1).as_array(),  out1.get_item(1).as_array())
+
+    def test_CompositionOperator_direct2(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id2)
+
+        out1 = G.direct(data)
+        
+        d_out = d.direct(data)
+
+        d1 = Id2.direct(data)
+        d2 = G.direct(d1)
+
+        numpy.testing.assert_array_almost_equal(d2.get_item(0).as_array(),
+                                                d_out.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(d2.get_item(1).as_array(),
+                                                d_out.get_item(1).as_array())
+
+    def test_CompositionOperator_direct3(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id2)
+
+        out1 = G.direct(data)
+        
+        d_out = d.direct(data)
+
+        d1 = Id2.direct(data)
+        d2 = G.direct(d1)
+
+        numpy.testing.assert_array_almost_equal(d2.get_item(0).as_array(),
+                                                d_out.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(d2.get_item(1).as_array(),
+                                                d_out.get_item(1).as_array())
+
+        G2Id = G.compose(2*Id2)
+        d2g = G2Id.direct(data)
+
+        numpy.testing.assert_array_almost_equal(d2g.get_item(0).as_array(),
+                                                2 * d_out.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(d2g.get_item(1).as_array(),
+                                                2 * d_out.get_item(1).as_array())
+
+    def test_CompositionOperator_direct4(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id1, Id2)
+
+        out1 = G.direct(data)
+        
+        d_out = d.direct(data)
+
+        numpy.testing.assert_array_almost_equal(d_out.get_item(0).as_array(),
+                                                2 * out1.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(d_out.get_item(1).as_array(),
+                                                2 * out1.get_item(1).as_array())
+
+    def test_CompositionOperator_adjoint1(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id2)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  out1.as_array())
+        
+    def test_CompositionOperator_adjoint2(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id1)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  2 * out1.as_array())
     
+    def test_CompositionOperator_adjoint3(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = G.compose(Id1)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  2 * out1.as_array())
+
+
+    def test_CompositionOperator_adjoint4(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        
+        d = G.compose(-Id1)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  - 2 * out1.as_array())
 
+    def test_CompositionOperator_adjoint5(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 3 * Identity(ig)
+        Id = Id1 - Identity(ig)
+        d = G.compose(Id)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  2 * out1.as_array())
+    
+    def test_CompositionOperator_adjoint6(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 3 * Identity(ig)
+        Id = ZeroOperator(ig)
+        d = G.compose(Id)
+        da = d.direct(data)
+        
+        out1 = G.adjoint(da)
+        out2 = d.adjoint(da)
+
+        numpy.testing.assert_array_almost_equal(out2.as_array(),  0 * out1.as_array())
+
+    def stest_CompositionOperator_direct4(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        sym = SymmetrizedGradient(domain_geometry=ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(sym, Id2)
+
+        out1 = G.direct(data)
+        out2 = d.direct(data)
+
+
+        numpy.testing.assert_array_almost_equal(out2.get_item(0).as_array(),  out1.get_item(0).as_array())
+        numpy.testing.assert_array_almost_equal(out2.get_item(1).as_array(),  out1.get_item(1).as_array())
+    
+    def test_CompositionOperator_adjoint7(self):
+        ig = self.ig
+        data = self.data
+        G = Gradient(domain_geometry=ig)
+        
+
+        Id1 = 2 * Identity(ig)
+        Id2 = Identity(ig)
+        
+        d = CompositionOperator(G, Id1, Id2)
+
+        out1 = G.direct(data)
+        out2 = G.adjoint(out1)
+        
+        d_out = d.adjoint(out1)
+
+        numpy.testing.assert_array_almost_equal(d_out.as_array(),
+                                                2 * out2.as_array())
+        numpy.testing.assert_array_almost_equal(d_out.as_array(),
+                                                2 * out2.as_array())
+
+
+