C lib (#458)

C library implemented with optimised axpy fucntions and gradient operator in c

11 files changed, 626 insertions, 35 deletions

diff --git a/Wrappers/Python/CMakeLists.txt b/Wrappers/Python/CMakeLists.txt
new file mode 100644
index 0000000..f325265
--- /dev/null
+++ b/Wrappers/Python/CMakeLists.txt
@@ -0,0 +1,88 @@
+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")
+  if (PYTHON_DEST_DIR)
+  set(PYTHON_DEST "${PYTHON_DEST_DIR}")
+  else() 
+  set(PYTHON_DEST "${CMAKE_INSTALL_PREFIX}/python")
+  endif()
+  message(STATUS "Python wrappers will be installed in " ${PYTHON_DEST})
+
+  message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}")
+
+  set(CMAKE_BUILD_TYPE "Release")
+  
+  find_package(PythonLibs)
+  if (PYTHONINTERP_FOUND)
+    message(STATUS "Found PYTHON_EXECUTABLE=${PYTHON_EXECUTABLE}")
+    message(STATUS "Python version ${PYTHON_VERSION_STRING}")
+  endif()
+  if (PYTHONLIBS_FOUND)
+    message(STATUS "Found PYTHON_INCLUDE_DIRS=${PYTHON_INCLUDE_DIRS}")
+    message(STATUS "Found PYTHON_LIBRARIES=${PYTHON_LIBRARIES}")
+  endif()
+
+  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(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})
+
+      message("Core binary dir " ${CMAKE_BINARY_DIR}/Core/${CMAKE_BUILD_TYPE})
+      
+      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}
+                        COMMAND ${CMAKE_COMMAND} -E env CIL_VERSION=${CIL_VERSION}
+                                                        #PREFIX=${CMAKE_SOURCE_DIR}/src/Core
+                                                        #LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/Core
+                                                        #LIBRARY_LIB=${CMAKE_BINARY_DIR}/src/Core
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} -vv install
+                                                        #echo "EDO"
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_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}
+                        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_ext --inplace
+                        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/
+                                                        LIBRARY_INC=${CMAKE_SOURCE_DIR}/src/include
+                                                        LIBRARY_LIB=${CMAKE_BINARY_DIR}/
+                                                        ${PYTHON_EXECUTABLE} ${SETUP_PY} build_ext --inplace
+                        COMMAND ${CMAKE_COMMAND} -E touch ${OUTPUT}
+                        DEPENDS cilacc
+                        )
+        endif()
+        #set (PYTHON_DEST ${CMAKE_INSTALL_PREFIX}/python/)
+        install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/
+                DESTINATION ${PYTHON_DEST})
+        #file(TOUCH ${PYTHON_DEST}/edo/__init__.py)
+        
+      endif()
+      
+      
+      add_custom_target(PythonWrapper ALL DEPENDS ${OUTPUT})
+      #install(CODE "execute_process(COMMAND ${PYTHON} ${SETUP_PY} install)")
+  endif()
+
+endif()
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 01ce7ef..1ab1d1e 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -26,6 +26,26 @@ from datetime import timedelta, datetime
 import warnings
 from functools import reduce
 from numbers import Number
+import ctypes, platform
+
+# dll = os.path.abspath(os.path.join( 
+#          os.path.abspath(os.path.dirname(__file__)),
+#          'libfdiff.dll')
+# )
+
+# check for the extension
+if platform.system() == 'Linux':
+    dll = 'libcilacc.so'
+elif platform.system() == 'Windows':
+    dll = 'cilacc.dll'
+elif platform.system() == 'Darwin':
+    dll = 'libcilacc.dylib'
+else:
+    raise ValueError('Not supported platform, ', platform.system())
+
+#print ("dll location", dll)
+cilacc = ctypes.cdll.LoadLibrary(dll)
+
 
 def find_key(dic, val):
     """return the key of dictionary dic given the value"""
@@ -802,7 +822,69 @@ class DataContainer(object):
     def minimum(self,x2, out=None, *args, **kwargs):
         return self.pixel_wise_binary(numpy.minimum, x2=x2, out=out, *args, **kwargs)
 
-    
+    @staticmethod
+    def axpby(a,x,b,y,out,dtype=numpy.float32):
+        '''performs axpby with cilacc C library
+        
+        Does the operation .. math:: a*x+b*y and stores the result in out
+
+        :param a: scalar
+        :param x: DataContainer
+        :param b: scalar
+        :param y: DataContainer
+        :param out: DataContainer to store the result
+        :param dtype: optional, data type of the DataContainers
+        '''
+
+        c_float_p = ctypes.POINTER(ctypes.c_float)
+        c_double_p = ctypes.POINTER(ctypes.c_double)
+        # get the reference to the data
+        ndx = x.as_array()
+        ndy = y.as_array()
+        ndout = out.as_array()
+
+        if ndx.dtype != dtype:
+            ndx = ndx.astype(dtype)
+        if ndy.dtype != dtype:
+            ndy = ndy.astype(dtype)
+        
+        if dtype == numpy.float32:
+            x_p = ndx.ctypes.data_as(c_float_p)
+            y_p = ndy.ctypes.data_as(c_float_p)
+            out_p = ndout.ctypes.data_as(c_float_p)
+            f = cilacc.saxpby
+
+        elif dtype == numpy.float64:
+            ndx = ndx.astype(numpy.float64)
+            b = b.astype(numpy.float64)
+            x_p = ndx.ctypes.data_as(c_double_p)
+            y_p = ndy.ctypes.data_as(c_double_p)
+            out_p = ndout.ctypes.data_as(c_double_p)
+            f = cilacc.daxpby
+        else:
+            raise TypeError('Unsupported type {}. Expecting numpy.float32 or numpy.float64'.format(dtype))
+
+        #out = numpy.empty_like(a)
+
+        
+        # int psaxpby(float * x, float * y, float * out, float a, float b, long size)
+        cilacc.saxpby.argtypes = [ctypes.POINTER(ctypes.c_float),  # pointer to the first array 
+                                  ctypes.POINTER(ctypes.c_float),  # pointer to the second array 
+                                  ctypes.POINTER(ctypes.c_float),  # pointer to the third array 
+                                  ctypes.c_float,                  # type of A (float)
+                                  ctypes.c_float,                  # type of B (float)
+                                  ctypes.c_long]                   # type of size of first array 
+        cilacc.daxpby.argtypes = [ctypes.POINTER(ctypes.c_double), # pointer to the first array 
+                                  ctypes.POINTER(ctypes.c_double), # pointer to the second array 
+                                  ctypes.POINTER(ctypes.c_double), # pointer to the third array 
+                                  ctypes.c_double,                 # type of A (c_double)
+                                  ctypes.c_double,                 # type of B (c_double)
+                                  ctypes.c_long]                   # type of size of first array 
+
+        if f(x_p, y_p, out_p, a, b, ndx.size) != 0:
+            raise RuntimeError('axpby execution failed')
+        
+
     ## unary operations
     def pixel_wise_unary(self, pwop, *args,  **kwargs):
         out = kwargs.get('out', None)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 3c32a93..8e07802 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -21,15 +21,37 @@ from __future__ import print_function
 from __future__ import unicode_literals
 
 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
-from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
+import warnings
 
-#%%
+NEUMANN = 'Neumann'
+PERIODIC = 'Periodic'
+C = 'c'
+NUMPY = 'numpy'
+CORRELATION_SPACE = "Space"
+CORRELATION_SPACECHANNEL = "SpaceChannels"
 
 class Gradient(LinearOperator):
+
+    """This is a class to compute the first-order forward/backward differences on ImageData
     
+    :param gm_domain: Set up the domain of the function
+    :type gm_domain: `ImageGeometry`
+    :param bnd_cond: Set the boundary conditions to use 'Neumann' or 'Periodic', defaults to 'Neumann'
+    :type bnd_cond: str, optional    
     
+    :param \**kwargs:
+        See below
+
+    :Keyword Arguments:
+        * *correlation* (``str``) --
+          'Space' or 'SpaceChannels', defaults to 'Space'
+        * *backend* (``str``) --
+          'c' or 'numpy', defaults to 'c' if correlation is 'SpaceChannels' or channels = 1
+    """
+
     r'''Gradient Operator: .. math:: \nabla : X -> Y           
             
             Computes first-order forward/backward differences 
@@ -39,29 +61,82 @@ class Gradient(LinearOperator):
                 Example (2D): u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u]
                               u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2
         
-                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-                Grad_order = ['channels', 'direction_y', 'direction_x']
-                Grad_order = ['direction_z', 'direction_y', 'direction_x']
-                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-                
-
     '''
+
+    #kept here for backwards compatability
+    CORRELATION_SPACE = CORRELATION_SPACE
+    CORRELATION_SPACECHANNEL = CORRELATION_SPACECHANNEL
+
+    def __init__(self, gm_domain, 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:
+            #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)
+        else:
+            self.operator = Gradient_C(gm_domain, bnd_cond=bnd_cond)
+
+
+    def direct(self, x, out=None):
+        """Computes the first-order forward differences
+
+        :param x: Image data
+        :type x: `ImageData`
+        :param out: pre-allocated output memory to store result
+        :type out: `BlockDataContainer`, optional        
+        :return: result data if not passed as parameter
+        :rtype: `BlockDataContainer`
+        """        
+        return self.operator.direct(x, out=out)
+        
+        
+    def adjoint(self, x, out=None):
+        """Computes the first-order backward differences
+
+        :param x: Gradient images for each dimension in ImageGeometry domain
+        :type x: `BlockDataContainer`
+        :param out: pre-allocated output memory to store result
+        :type out: `ImageData`, optional      
+        :return: result data if not passed as parameter
+        :rtype: `ImageData`
+        """            
+        return self.operator.adjoint(x, out=out)
+
+    def domain_geometry(self):
+        '''Returns domain_geometry of Gradient'''
+        
+        return self.operator.gm_domain
     
-                    
-    CORRELATION_SPACE = "Space"
-    CORRELATION_SPACECHANNEL = "SpaceChannels"
+    def range_geometry(self):
+        
+        '''Returns range_geometry of Gradient'''
+        
+        return self.operator.gm_range
 
+class Gradient_numpy(LinearOperator):
+    
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
-        super(Gradient, self).__init__() 
+        super(Gradient_numpy, self).__init__() 
                 
         self.gm_domain = gm_domain # Domain of Grad Operator
         
-        self.correlation = kwargs.get('correlation',Gradient.CORRELATION_SPACE)
+        self.correlation = kwargs.get('correlation',CORRELATION_SPACE)
         
-        if self.correlation==Gradient.CORRELATION_SPACE:
-            if self.gm_domain.channels>1:
+        if self.correlation==CORRELATION_SPACE:
+            if self.gm_domain.channels > 1:
                 self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length-1)] )
+
                 if self.gm_domain.length == 4:
                     # 3D + Channel
                     # expected Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
@@ -70,6 +145,7 @@ class Gradient(LinearOperator):
                     # 2D + Channel
                     # expected Grad_order = ['channels', 'direction_y', 'direction_x']
                     expected_order = [ImageGeometry.CHANNEL, ImageGeometry.HORIZONTAL_Y, ImageGeometry.HORIZONTAL_X]
+
                 order = self.gm_domain.get_order_by_label(self.gm_domain.dimension_labels, expected_order)
                 self.ind = order[1:]
                 #self.ind = numpy.arange(1,self.gm_domain.length)
@@ -83,10 +159,12 @@ class Gradient(LinearOperator):
                 else:
                     # 2D
                     expected_order = [ImageGeometry.HORIZONTAL_Y, ImageGeometry.HORIZONTAL_X]    
+
                 self.ind = self.gm_domain.get_order_by_label(self.gm_domain.dimension_labels, expected_order)
                 # self.ind = numpy.arange(self.gm_domain.length)
-        elif self.correlation==Gradient.CORRELATION_SPACECHANNEL:
-            if self.gm_domain.channels>1:
+                
+        elif self.correlation==CORRELATION_SPACECHANNEL:
+            if self.gm_domain.channels > 1:
                 self.gm_range = BlockGeometry(*[self.gm_domain for _ in range(self.gm_domain.length)])
                 self.ind = range(self.gm_domain.length)
             else:
@@ -189,6 +267,146 @@ class Gradient(LinearOperator):
             res.append(spMat.sum_abs_col())
         return BlockDataContainer(*res)
    
+import ctypes, platform
+
+# check for the extension
+if platform.system() == 'Linux':
+    dll = 'libcilacc.so'
+elif platform.system() == 'Windows':
+    dll = 'cilacc.dll'
+elif platform.system() == 'Darwin':
+    dll = 'libcilacc.dylib'
+else:
+    raise ValueError('Not supported platform, ', platform.system())
+
+#print ("dll location", dll)
+cilacc = ctypes.cdll.LoadLibrary(dll)
+
+#FD = ctypes.CDLL(dll)
+
+c_float_p = ctypes.POINTER(ctypes.c_float)
+
+cilacc.openMPtest.restypes = ctypes.c_int32
+
+cilacc.fdiff4D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+cilacc.fdiff3D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+cilacc.fdiff2D.argtypes = [ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.POINTER(ctypes.c_float),
+                       ctypes.c_long,
+                       ctypes.c_long,
+                       ctypes.c_int32,
+                       ctypes.c_int32]
+
+
+class Gradient_C(LinearOperator):
+    
+    '''Finite Difference Operator:
+            
+            Computes first-order forward/backward differences 
+                     on 2D, 3D, 4D ImageData
+                     under Neumann/Periodic boundary conditions'''
+
+    def __init__(self, gm_domain, gm_range=None, bnd_cond = NEUMANN):
+
+        super(Gradient_C, self).__init__() 
+
+        self.gm_domain = gm_domain
+        self.gm_range = gm_range
+        
+        #default is 'Neumann'
+        self.bnd_cond = 0
+        
+        if bnd_cond == PERIODIC:
+            self.bnd_cond = 1
+        
+        # Domain Geometry = Range Geometry if not stated
+        if self.gm_range is None:
+            self.gm_range = BlockGeometry(*[gm_domain for _ in range(len(gm_domain.shape))])
+        
+
+        if len(gm_domain.shape) == 4:
+            self.fd = cilacc.fdiff4D
+        elif len(gm_domain.shape) == 3:
+            self.fd = cilacc.fdiff3D
+        elif len(gm_domain.shape) == 2:
+            self.fd = cilacc.fdiff2D
+        else:
+            raise ValueError('Number of dimensions not supported, expected 2, 3 or 4, got {}'.format(len(gm_domain.shape)))
+        
+                        
+    @staticmethod 
+    def datacontainer_as_c_pointer(x):
+        ndx = x.as_array()
+        return ndx, ndx.ctypes.data_as(c_float_p)
+        
+    def direct(self, x, out=None):
+        ndx , x_p = Gradient_C.datacontainer_as_c_pointer(x)
+        
+        return_val = False
+        if out is None:
+            out = self.gm_range.allocate(None)
+            return_val = True
+
+        arg1 = [Gradient_C.datacontainer_as_c_pointer(out.get_item(i))[1] for i in range(self.gm_range.shape[0])]
+        arg2 = [el for el in x.shape]
+        args = arg1 + arg2 + [self.bnd_cond, 1]
+        self.fd(x_p, *args)
+        
+        if return_val is True:
+            return out
+
+
+    def adjoint(self, x, out=None):
+
+        return_val = False
+        if out is None:
+            out = self.gm_domain.allocate(None)
+            return_val = True
+
+        ndout, out_p = Gradient_C.datacontainer_as_c_pointer(out)
+
+        arg1 = [Gradient_C.datacontainer_as_c_pointer(x.get_item(i))[1] for i in range(self.gm_range.shape[0])]
+        arg2 = [el for el in out.shape]
+        args = arg1 + arg2 + [self.bnd_cond, 0]
+
+        self.fd(out_p, *args)
+
+        if return_val is True:
+            return out
+
+    def domain_geometry(self):
+        
+        '''Returns domain_geometry of Gradient'''
+        
+        return self.gm_domain
+    
+    def range_geometry(self):
+        
+        '''Returns range_geometry of Gradient'''
+        
+        return self.gm_range
+
        
 if __name__ == '__main__':
     
@@ -302,8 +520,5 @@ if __name__ == '__main__':
     G.direct(u, out = res)          
     G.adjoint(w, out = res1)
           
-    print( (res*w).sum(),  (u*res1).sum() )     
-    
+    print( (res*w).sum(),  (u*res1).sum() )
     
-    
-
diff --git a/Wrappers/Python/conda-recipe/bld.bat b/Wrappers/Python/recipe/bld.bat
index 97a4e62..97a4e62 100644
--- a/Wrappers/Python/conda-recipe/bld.bat
+++ b/Wrappers/Python/recipe/bld.bat
diff --git a/Wrappers/Python/conda-recipe/build.sh b/Wrappers/Python/recipe/build.sh
index 43e85d5..43e85d5 100644
--- a/Wrappers/Python/conda-recipe/build.sh
+++ b/Wrappers/Python/recipe/build.sh
diff --git a/Wrappers/Python/conda-recipe/conda_build_config.yaml b/Wrappers/Python/recipe/conda_build_config.yaml
index e5e168b..e5e168b 100644
--- a/Wrappers/Python/conda-recipe/conda_build_config.yaml
+++ b/Wrappers/Python/recipe/conda_build_config.yaml
diff --git a/Wrappers/Python/conda-recipe/meta.yaml b/Wrappers/Python/recipe/meta.yaml
index 9d03220..9d03220 100644
--- a/Wrappers/Python/conda-recipe/meta.yaml
+++ b/Wrappers/Python/recipe/meta.yaml
diff --git a/Wrappers/Python/setup-cil.py.in b/Wrappers/Python/setup-cil.py.in
new file mode 100644
index 0000000..d1cf7fc
--- /dev/null
+++ b/Wrappers/Python/setup-cil.py.in
@@ -0,0 +1,70 @@
+# -*- coding: utf-8 -*-
+#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+
+#   Copyright 2017 UKRI-STFC
+#   Copyright 2017 University of Manchester
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#   http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from distutils.core import setup
+import os
+import sys
+
+
+cil_version=os.environ['CIL_VERSION']
+if  cil_version == '':
+    print("Please set the environmental variable CIL_VERSION")
+    sys.exit(1)
+
+setup(
+    name="ccpi-framework",
+    version=cil_version,
+    packages=['ccpi' , 'ccpi.io',
+              'ccpi.framework', 'ccpi.optimisation', 
+              'ccpi.optimisation.operators',
+              'ccpi.optimisation.algorithms',
+              'ccpi.optimisation.functions',
+              'ccpi.processors',
+              'ccpi.utilities', 'ccpi.utilities.jupyter',
+              'ccpi.contrib','ccpi.contrib.optimisation',
+              'ccpi.contrib.optimisation.algorithms'],
+    data_files = [('share/ccpi', ['data/boat.tiff',
+                                  'data/peppers.tiff',
+                                  'data/camera.png',
+                                  'data/resolution_chart.tiff',
+                                  'data/shapes.png',
+                                  'data/24737_fd_normalised.nxs'])],
+
+    # Project uses reStructuredText, so ensure that the docutils get
+    # installed or upgraded on the target machine
+    #install_requires=['docutils>=0.3'],
+
+#    package_data={
+#        # If any package contains *.txt or *.rst files, include them:
+#        '': ['*.txt', '*.rst'],
+#        # And include any *.msg files found in the 'hello' package, too:
+#        'hello': ['*.msg'],
+#    },
+    # zip_safe = False,
+
+    # metadata for upload to PyPI
+    author="CCPi developers",
+    maintainer="Edoardo Pasca",
+    maintainer_email="edoardo.pasca@stfc.ac.uk",
+    description='CCPi Core Imaging Library - Python Framework Module',
+    license="Apache v2.0",
+    keywords="Python Framework",
+    url="http://www.ccpi.ac.uk/CIL",   # project home page, if any
+
+    # could also include long_description, download_url, classifiers, etc.
+)
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 9b9a8dd..d2d5b05 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -732,6 +732,19 @@ class TestDataContainer(unittest.TestCase):
         u.multiply(2, out=u)
         c = b * 2
         numpy.testing.assert_array_equal(u.as_array(), c)
+
+    def test_axpby(self):
+        print ("test axpby")
+        ig = ImageGeometry(10,10)                                               
+        d1 = ig.allocate(1)                                                     
+        d2 = ig.allocate(2)                                                     
+        out = ig.allocate(None)                                                 
+        # equals to 2 * [1] + 1 * [2] = [4]
+        DataContainer.axpby(2,d1,1,d2,out)
+        res = numpy.ones_like(d1.as_array()) * 4.
+        numpy.testing.assert_array_equal(res, out.as_array())
+
+
         
 
 
diff --git a/Wrappers/Python/test/test_Gradient.py b/Wrappers/Python/test/test_Gradient.py
index 5dc8137..5aeede0 100755
--- a/Wrappers/Python/test/test_Gradient.py
+++ b/Wrappers/Python/test/test_Gradient.py
@@ -23,6 +23,7 @@ from ccpi.framework import BlockDataContainer
 import functools
 
 from ccpi.optimisation.operators import Gradient, Identity, BlockOperator
+from ccpi.optimisation.operators import LinearOperator
 
 class TestGradient(unittest.TestCase):
     def test_Gradient(self): 
@@ -36,17 +37,17 @@ class TestGradient(unittest.TestCase):
         ig3 = ImageGeometry(voxel_num_x = M, voxel_num_y = N, channels = channels) 
         ig4 = ImageGeometry(voxel_num_x = M, voxel_num_y = N, channels = channels, voxel_num_z= K) 
         
-        G1 = Gradient(ig1, correlation = 'Space') 
+        G1 = Gradient(ig1, correlation = 'Space', backend='numpy') 
         print(G1.range_geometry().shape, '2D no channels')
             
-        G4 = Gradient(ig3, correlation = 'SpaceChannels')
-        print(G4.range_geometry().shape, '2D with channels corr')
-        G5 = Gradient(ig3, correlation = 'Space')
+        G4 = Gradient(ig3, correlation = 'SpaceChannels', backend='numpy')
+        print(G4.range_geometry().shape, '2D with channels corr')  
+        G5 = Gradient(ig3, correlation = 'Space', backend='numpy')
         print(G5.range_geometry().shape, '2D with channels no corr')
         
-        G6 = Gradient(ig4, correlation = 'Space')
+        G6 = Gradient(ig4, correlation = 'Space', backend='numpy')
         print(G6.range_geometry().shape, '3D with channels no corr')
-        G7 = Gradient(ig4, correlation = 'SpaceChannels')
+        G7 = Gradient(ig4, correlation = 'SpaceChannels', backend='numpy')
         print(G7.range_geometry().shape, '3D with channels with corr')
         
         
@@ -77,8 +78,8 @@ class TestGradient(unittest.TestCase):
         #check direct/adjoint for space/channels correlation
         
         ig_channel = ImageGeometry(voxel_num_x = 2, voxel_num_y = 3, channels = 2)
-        G_no_channel = Gradient(ig_channel, correlation = 'Space')
-        G_channel = Gradient(ig_channel, correlation = 'SpaceChannels')
+        G_no_channel = Gradient(ig_channel, correlation = 'Space', backend='numpy')
+        G_channel = Gradient(ig_channel, correlation = 'SpaceChannels', backend='numpy')
         
         u3 = ig_channel.allocate('random_int')
         res_no_channel = G_no_channel.direct(u3)
@@ -95,7 +96,7 @@ class TestGradient(unittest.TestCase):
         
         ig2D = ImageGeometry(voxel_num_x = 2, voxel_num_y = 3)
         u4 = ig2D.allocate('random_int')
-        G2D = Gradient(ig2D)
+        G2D = Gradient(ig2D, backend='numpy')
         res = G2D.direct(u4)  
         print(res[0].as_array())
         print(res[1].as_array())
@@ -105,12 +106,133 @@ class TestGradient(unittest.TestCase):
         arr = ig.allocate('random_int' )
         
         # check direct of Gradient and sparse matrix
-        G = Gradient(ig)
+        G = Gradient(ig, backend='numpy')
         norm1 = G.norm(iterations=300)
         print ("should be sqrt(8) {} {}".format(numpy.sqrt(8), norm1))
         numpy.testing.assert_almost_equal(norm1, numpy.sqrt(8), decimal=1)
         ig4 = ImageGeometry(M,N, channels=3)
-        G4 = Gradient(ig4, correlation=Gradient.CORRELATION_SPACECHANNEL)
+        G4 = Gradient(ig4, correlation="SpaceChannels", backend='numpy')
         norm4 = G4.norm(iterations=300)
-        print ("should be sqrt(12) {} {}".format(numpy.sqrt(12), norm4))
+        print("should be sqrt(12) {} {}".format(numpy.sqrt(12), norm4))
         self.assertTrue((norm4 - numpy.sqrt(12))/norm4 < 0.2)
+
+    def test_GradientOperator_4D(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        size = nc * nz * ny * nx
+        dim = [nc, nz, ny, nx]
+
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+        arr = numpy.arange(size).reshape(dim).astype(numpy.float32)**2
+
+        data = ig.allocate()
+        data.fill(arr)
+
+        #neumann
+        grad_py = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad_py.direct(data)
+        gold_adjoint = grad_py.adjoint(gold_direct)
+
+        grad_c = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        out_direct = grad_c.direct(data)
+        out_adjoint = grad_c.adjoint(out_direct)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+        #periodic
+        grad_py = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad_py.direct(data)
+        gold_adjoint = grad_py.adjoint(gold_direct)
+
+        grad_c = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='c')
+        out_direct = grad_c.direct(data)
+        out_adjoint = grad_c.adjoint(out_direct)
+
+        #print("Gradient, 4D, bnd_cond='Periodic', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("Gradient, 4D, bnd_cond='Periodic', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+    def test_Gradient_4D_allocate(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        size = nc * nz * ny * nx
+        dim = [nc, nz, ny, nx]
+
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+
+        arr = numpy.arange(size).reshape(dim).astype(numpy.float32)**2
+
+        data = ig.allocate()
+        data.fill(arr)
+
+        #numpy
+        grad1 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        gold_direct = grad1.direct(data)
+        gold_adjoint = grad1.adjoint(gold_direct)
+
+        grad2 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        out_direct = grad2.range_geometry().allocate()
+        out_adjoint = grad2.domain_geometry().allocate()
+        grad2.direct(data, out=out_direct)
+        grad2.adjoint(out_direct, out=out_adjoint)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+        #c
+        grad1 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        gold_direct = grad1.direct(data)
+        gold_adjoint = grad1.adjoint(gold_direct)
+
+        grad2 = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        out_direct = grad2.range_geometry().allocate()
+        out_adjoint = grad2.domain_geometry().allocate()
+        grad2.direct(data, out=out_direct)
+        grad2.adjoint(out_direct, out=out_adjoint)
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', direct")
+        numpy.testing.assert_array_equal(out_direct.get_item(0).as_array(), gold_direct.get_item(0).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(1).as_array(), gold_direct.get_item(1).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(2).as_array(), gold_direct.get_item(2).as_array())
+        numpy.testing.assert_array_equal(out_direct.get_item(3).as_array(), gold_direct.get_item(3).as_array())
+
+        #print("GradientOperator, 4D, bnd_cond='Neumann', adjoint")
+        numpy.testing.assert_array_equal(out_adjoint.as_array(), gold_adjoint.as_array())
+
+    def test_Gradient_linearity(self):
+
+        nc, nz, ny, nx = 3, 4, 5, 6
+        ig = ImageGeometry(voxel_num_x=nx, voxel_num_y=ny, voxel_num_z=nz, channels=nc)
+
+        grad = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='c')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='c')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Neumann', correlation='SpaceChannels', backend='numpy')
+        self.assertTrue(LinearOperator.dot_test(grad))
+
+        grad = Gradient(ig, bnd_cond='Periodic', correlation='SpaceChannels', backend='numpy')
+        self.assertTrue(LinearOperator.dot_test(grad))
+        
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 3372b9b..96c48da 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -494,12 +494,13 @@ class TestBlockOperator(unittest.TestCase):
             print("Z1", Z1[0][1].as_array())
             print("RES1", RES1[0][1].as_array())
     def test_timedifference(self):
+
         print ("test_timedifference")
         M, N ,W = 100, 512, 512
         ig = ImageGeometry(M, N, W)
         arr = ig.allocate('random_int')  
         
-        G = Gradient(ig)
+        G = Gradient(ig, backend='numpy')
         Id = Identity(ig)
         
         B = BlockOperator(G, Id)