Function m (#482)

* fix func style

* update pnorm not quite right

* add test for SumFunction

* passed test for SumFunction

* passed test for SumFunction

* passed test for SumFunction

* rename Norm2Sq to LeastSquares

* rename Norm2Sq to LeastSquares

* delete files and move to Function.py

* fix accuracy convex conjugate

* comment proximal conjugate, have generic expression

* remove some imports

* remove imports and comment gradient, proximal conjugate , rmul

* remove imports and comment prox conj

* fix init

* fix out gradient

* fix init

* SumFunction full test

* rename Norm2Sq

* rename Norm2Sq

* change kwargs in KL, not ready yet

* change def of KL cases

* change kwargs in KL, not ready yet

* LinearOperatorMatrix FISTA cvx is not working

* change name in KL

* minor fix

* BlockDataContainer should give a geometry when shape is the same

* change convex conj when we have BDC

* remove space

* fix out for prox

* rename x to b

* add SumFunctionScalar test

* fix convex conj of constant

* fix convex conj of constant

* minor fix

* fix con conj

* add new tests

* KL not ready yet

* add docs, not ready yet

* fix documentation

* fix docs

* fix doc scaled function

* fix doc translate function

* fix modules path for documentation

* minor doc fix

* fix doc prox conj for KL

* fix vector geometry in vectordata

* fix direct/adjpint LinearOperator

* add test LinearOperatorMatrix

* fix doc and check methods for functionOperaCom

* add new test for FunctionOperCom with LinearOperMatr

* get vector geometry in vector data

* minor changes and docs

* fix docs

* check Lips cons

* rename Norm2Sq to LeastSquares for docs

* fix KL, docs, background term, conve conj

* add KL test

* fix KL

* minor fix doc

* add back prox conj work with smoothTV

* smooth TV not ready yet

* add prox conj back

* minor change

* no change

* fixing docs, not ready yet

* uncomment prox_conj for ScaledFunction

* fix doc for BlockFunction

* add smoothTV not ready yet

* no change

* fix doc and uncomment prox conj

* remove unicode literals

* sphnix is a nightmare

* fixing ADMM

* remove conv conj, prox, prox conj from smooth TV

* merge edo's docs

* fix docs for L2Norm after edo's PR

* fix docs

* remove ADMM

* add docs for Function module

* remove unicode_literals from future imports

https://stackoverflow.com/questions/7602171/unicode-error-unicodeescape-codec-cant-decode-bytes-string-with-u

* remove ADMM

* removed unicode literals from future imports

* fix init for functions

* remove main

* add explcit KL prox conj, does not work with generic prox conj of function

* add elif pass for vector geometry

* fixes from comments on the PR

* minor changes

* fix mul in Function and Identity in TranslateFunction

* remove from call

* remove BlockGeometr from BDC

* remove min, max

* cleaned Function

* add test for constant function

* remove smooth mixed L21 Norm

* fixed fill method and test of SumFunction

* Functions dev (#479)

* v19.10 docs (#467)

updated docstrings and documentation

* add Gemma in the authors

* update docs

* Installs with make install (#478)

* closes #477

* install from build dir

Co-authored-by: Vaggelis Papoutsellis <22398586+epapoutsellis@users.noreply.github.com>

51 files changed, 2072 insertions, 995 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index 38c35f7..22cee03 100755..100644
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 from numbers import Number
@@ -61,7 +60,11 @@ class BlockDataContainer(object):
     def __init__(self, *args, **kwargs):
         ''''''
         self.containers = args
-        self.index = 0
+        self.index = 0        
+        self.geometry = None
+        #if len(set([i.shape for i in self.containers])):
+        #    self.geometry = self.containers[0].geometry
+                
         shape = kwargs.get('shape', None)
         if shape is None:
            shape = (len(args),1)
@@ -330,8 +333,9 @@ class BlockDataContainer(object):
                         
         if p==1:            
             return sum(self.abs())        
-        elif p==2:                    
-            return sum([el*el for el in self.containers]).sqrt()
+        elif p==2:                 
+            tmp = functools.reduce(lambda a,b: a + b*b, self.containers, self.get_item(0) * 0 ).sqrt()            
+            return tmp      
         else:
             return ValueError('Not implemented')
                 
@@ -507,26 +511,32 @@ if __name__ == '__main__':
     import numpy
 
     N, M = 2, 3
-    ig = ImageGeometry(N, M)    
-    BG = BlockGeometry(ig, ig)
+    ig = ImageGeometry(N, M) 
+    
+    ig1 = ImageGeometry(2*N, 4*M) 
+    BG = BlockGeometry(ig, ig1)
     
     U = BG.allocate('random_int')
     V = BG.allocate('random_int')
     
+    print(U.geometry)
     
-    print ("test sum BDC " )
-    w = U[0].as_array() + U[1].as_array()    
-    w1 = sum(U).as_array()    
-    numpy.testing.assert_array_equal(w, w1)
-    
-    print ("test sum BDC " )
-    z = numpy.sqrt(U[0].as_array()**2 + U[1].as_array()**2)
-    z1 = sum(U**2).sqrt().as_array()    
-    numpy.testing.assert_array_equal(z, z1)   
+    print(len(set([i.shape for i in U.containers]))==1)
     
-    z2 = U.pnorm(2)
     
-    zzz = U.dot(V)
+#    print ("test sum BDC " )
+#    w = U[0].as_array() + U[1].as_array()    
+#    w1 = sum(U).as_array()    
+#    numpy.testing.assert_array_equal(w, w1)
+#    
+#    print ("test sum BDC " )
+#    z = numpy.sqrt(U[0].as_array()**2 + U[1].as_array()**2)
+#    z1 = sum(U**2).sqrt().as_array()    
+#    numpy.testing.assert_array_equal(z, z1)   
+#    
+#    z2 = U.pnorm(2)
+#    
+#    zzz = U.dot(V)
     
    
 
diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py
index 988801f..2f28507 100755..100644
--- a/Wrappers/Python/ccpi/framework/BlockGeometry.py
+++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import

 from __future__ import division

 from __future__ import print_function

-from __future__ import unicode_literals

 

 import numpy

 from numbers import Number

@@ -94,9 +93,6 @@ class BlockGeometry(object):
                 containers[6]=containers[9]

                 containers[7]=containers[10]

                 containers[11]=containers[15]

-                

-                

-                

-        

+

         return BlockDataContainer(*containers)

            

diff --git a/Wrappers/Python/ccpi/framework/TestData.py b/Wrappers/Python/ccpi/framework/TestData.py
index c035850..a9df2ad 100755..100644
--- a/Wrappers/Python/ccpi/framework/TestData.py
+++ b/Wrappers/Python/ccpi/framework/TestData.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.framework import ImageData, ImageGeometry, DataContainer
 import numpy
diff --git a/Wrappers/Python/ccpi/framework/__init__.py b/Wrappers/Python/ccpi/framework/__init__.py
index 94788a5..213750b 100755..100644
--- a/Wrappers/Python/ccpi/framework/__init__.py
+++ b/Wrappers/Python/ccpi/framework/__init__.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 import sys
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 1ab1d1e..968c69c 100755..100644
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 import sys
@@ -561,7 +560,12 @@ class DataContainer(object):
                     return VectorData(cleaned)
     
     def fill(self, array, **dimension):
-        '''fills the internal numpy array with the one provided'''
+        '''fills the internal numpy array with the one provided
+        
+        :param array: numpy array to copy into the DataContainer
+        :type array: DataContainer, numpy array or number
+        :param dimension: dictionary, optional
+        '''
         if dimension == {}:
             if issubclass(type(array), DataContainer) or\
                issubclass(type(array), numpy.ndarray):
@@ -574,6 +578,8 @@ class DataContainer(object):
                 else:
                     #self.array[:] = array
                     numpy.copyto(self.array, array)
+            else:
+                self.array.fill(array)
         else:
             
             command = 'self.array['
@@ -814,8 +820,8 @@ class DataContainer(object):
         return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs)
     
     def power(self, other, *args, **kwargs):
-        return self.pixel_wise_binary(numpy.power, other, *args, **kwargs)
-    
+        return self.pixel_wise_binary(numpy.power, other, *args, **kwargs)    
+          
     def maximum(self, x2, *args, **kwargs):
         return self.pixel_wise_binary(numpy.maximum, x2, *args, **kwargs)
     
@@ -911,6 +917,12 @@ class DataContainer(object):
     def abs(self, *args,  **kwargs):
         return self.pixel_wise_unary(numpy.abs, *args,  **kwargs)
     
+#    def max(self, *args,  **kwargs):
+#        return self.pixel_wise_unary(numpy.max, *args,  **kwargs) 
+#    
+#    def min(self, *args,  **kwargs):
+#        return self.pixel_wise_unary(numpy.min, *args,  **kwargs)     
+    
     def sign(self, *args,  **kwargs):
         return self.pixel_wise_unary(numpy.sign, *args,  **kwargs)
     
@@ -1624,7 +1636,8 @@ class VectorData(DataContainer):
                 else:
                     raise ValueError('Incompatible size: expecting {} got {}'.format((self.length,), array.shape))
         deep_copy = True
-        super(VectorData, self).__init__(out, deep_copy, None)
+        # need to pass the geometry, othewise None
+        super(VectorData, self).__init__(out, deep_copy, None, geometry = self.geometry)
 
 class VectorGeometry(object):
     '''Geometry describing VectorData to contain 1D array'''
@@ -1661,30 +1674,36 @@ class VectorGeometry(object):
                     numpy.random.seed(seed)
                 max_value = kwargs.get('max_value', 100)
                 out.fill(numpy.random.randint(max_value,size=self.shape))
+            elif value is None:
+                pass
             else:
                 raise ValueError('Value {} unknown'.format(value))
         return out
 
     
 if __name__ == "__main__":
-
-    ig = ImageGeometry(voxel_num_x=100, 
-                    voxel_num_y=200, 
-                    voxel_num_z=300, 
-                    voxel_size_x=1, 
-                    voxel_size_y=1, 
-                    voxel_size_z=1, 
-                    center_x=0, 
-                    center_y=0, 
-                    center_z=0, 
-                    channels=50)
-
-    id = ig.allocate(2)
-
-    print(id.geometry)
-    print(id.dimension_labels)
-
-    sid = id.subset(channel = 20)
-
-    print(sid.dimension_labels)
-    print(sid.geometry)
+    
+    
+    vg = VectorGeometry(10)    
+    b = vg.allocate('random_int')
+
+#    ig = ImageGeometry(voxel_num_x=100, 
+#                    voxel_num_y=200, 
+#                    voxel_num_z=300, 
+#                    voxel_size_x=1, 
+#                    voxel_size_y=1, 
+#                    voxel_size_z=1, 
+#                    center_x=0, 
+#                    center_y=0, 
+#                    center_z=0, 
+#                    channels=50)
+#
+#    id = ig.allocate(2)
+#
+#    print(id.geometry)
+#    print(id.dimension_labels)
+#
+#    sid = id.subset(channel = 20)
+#
+#    print(sid.dimension_labels)
+#    print(sid.geometry)
diff --git a/Wrappers/Python/ccpi/io/NEXUSDataReader.py b/Wrappers/Python/ccpi/io/NEXUSDataReader.py
index dbada02..e696ca5 100644
--- a/Wrappers/Python/ccpi/io/NEXUSDataReader.py
+++ b/Wrappers/Python/ccpi/io/NEXUSDataReader.py
@@ -15,11 +15,13 @@
 #   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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 import numpy
 import os
 from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageData, ImageGeometry
 
-
 h5pyAvailable = True
 try:
     import h5py 
diff --git a/Wrappers/Python/ccpi/io/NEXUSDataWriter.py b/Wrappers/Python/ccpi/io/NEXUSDataWriter.py
index c0b4fae..926cb07 100644
--- a/Wrappers/Python/ccpi/io/NEXUSDataWriter.py
+++ b/Wrappers/Python/ccpi/io/NEXUSDataWriter.py
@@ -15,7 +15,9 @@
 #   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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 import numpy
 import os
 from ccpi.framework import AcquisitionData, AcquisitionGeometry, ImageData, ImageGeometry
diff --git a/Wrappers/Python/ccpi/io/NikonDataReader.py b/Wrappers/Python/ccpi/io/NikonDataReader.py
index 703b65b..ea8179e 100644
--- a/Wrappers/Python/ccpi/io/NikonDataReader.py
+++ b/Wrappers/Python/ccpi/io/NikonDataReader.py
@@ -1,10 +1,23 @@
-#!/usr/bin/env python3
 # -*- coding: utf-8 -*-
-"""
-Created on Wed Apr  3 10:30:25 2019
+#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
 
-@author: evelina
-"""
+#   Copyright 2019 UKRI-STFC
+#   Copyright 2019 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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 
 from ccpi.framework import AcquisitionData, AcquisitionGeometry
 import numpy
diff --git a/Wrappers/Python/ccpi/io/reader.py b/Wrappers/Python/ccpi/io/reader.py
index 8282fe9..52e06bd 100644
--- a/Wrappers/Python/ccpi/io/reader.py
+++ b/Wrappers/Python/ccpi/io/reader.py
@@ -19,7 +19,6 @@
 from __future__ import absolute_import

 from __future__ import division

 from __future__ import print_function

-from __future__ import unicode_literals

 

 from ccpi.framework import AcquisitionGeometry

 from ccpi.framework import AcquisitionData

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index 48d109e..48d109e 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 53804c5..64c2282 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -20,7 +20,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.algorithms import Algorithm
 import numpy
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 15a289d..9705682 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -20,7 +20,7 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
+
 
 from ccpi.optimisation.algorithms import Algorithm
 from ccpi.optimisation.functions import ZeroFunction
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
index 5e7284a..d30604e 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
@@ -23,7 +23,7 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
+
 import numpy
 from ccpi.optimisation.algorithms import Algorithm
 
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 8776875..cc384e3 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -20,7 +20,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.algorithms import Algorithm
 
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
index a59ce5f..8afecc8 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
@@ -23,7 +23,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.algorithms import Algorithm
 
diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
index 57592cd..5f9a881 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
@@ -23,7 +23,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.functions import Function
 from ccpi.framework import BlockDataContainer
@@ -31,14 +30,24 @@ from numbers import Number
 
 class BlockFunction(Function):
     
-    r'''BlockFunction acts as a separable sum function: f = [f_1,...,f_n]
+    r""" BlockFunction represents a *separable sum* function :math:`F` defined as
     
-      .. math::
-
-          f([x_1,...,x_n]) = f_1(x_1) +  .... + f_n(x_n)
-      |
-
-    '''
+    .. math:: F:X_{1}\times X_{2}\cdots\times X_{m} \rightarrow (-\infty, \infty]
+    
+    where :math:`F` is the separable sum of functions :math:`(f_{i})_{i=1}^{m}`,
+    
+    .. math:: F(x_{1}, x_{2}, \cdots, x_{m}) = \overset{m}{\underset{i=1}{\sum}}f_{i}(x_{i}), \mbox{ with } f_{i}: X_{i} \rightarrow (-\infty, \infty].
+    
+    A nice property (due to it's separability structure) is that the proximal operator 
+    can be decomposed along the proximal operators of each function :math:`f_{i}`.
+    
+    .. math:: \mathrm{prox}_{\tau F}(x) = ( \mathrm{prox}_{\tau f_{i}}(x_{i}) )_{i=1}^{m}
+    
+    In addition, if :math:`\tau := (\tau_{1},\dots,\tau_{m})`, then 
+    
+    .. math:: \mathrm{prox}_{\tau F}(x) = ( \mathrm{prox}_{\tau_{i} f_{i}}(x_{i}) )_{i=1}^{m}
+    
+    """
     
     def __init__(self, *functions):
                 
@@ -48,13 +57,17 @@ class BlockFunction(Function):
                                 
     def __call__(self, x):
         
-        r'''Evaluates the BlockFunction at a BlockDataContainer x
+        r""" Returns the value of the BlockFunction :math:`F`
         
-            :param: x (BlockDataContainer): must have as many rows as self.length
+        .. math:: F(x) = \overset{m}{\underset{i=1}{\sum}}f_{i}(x_{i}), \mbox{ where } x = (x_{1}, x_{2}, \cdots, x_{m}), \quad i = 1,2,\dots,m
+                
+        Parameter:
+            
+            x : BlockDataContainer and must have as many rows as self.length
 
             returns ..math:: \sum(f_i(x_i))
             
-        '''
+        """
         
         if self.length != x.shape[0]:
             raise ValueError('BlockFunction and BlockDataContainer have incompatible size')
@@ -65,53 +78,36 @@ class BlockFunction(Function):
     
     def convex_conjugate(self, x):
         
-        r'''Convex conjugate of BlockFunction at x            
+        r"""Returns the value of the convex conjugate of the BlockFunction at :math:`x^{*}`.       
         
-            .. math:: returns \sum(f_i^{*}(x_i))
+            .. math:: F^{*}(x^{*}) = \overset{m}{\underset{i=1}{\sum}}f_{i}^{*}(x^{*}_{i})
+            
+            Parameter:
+            
+                x : BlockDataContainer and must have as many rows as self.length            
         
-        '''       
-        t = 0                
+        """     
+        
+        if self.length != x.shape[0]:
+            raise ValueError('BlockFunction and BlockDataContainer have incompatible size')
+        t = 0              
         for i in range(x.shape[0]):
             t += self.functions[i].convex_conjugate(x.get_item(i))
         return t  
     
-    
-    def proximal_conjugate(self, x, tau, out = None):
-        
-        r'''Proximal operator of BlockFunction at x: 
-                 
-                 .. math:: prox_{tau*f}(x) = \sum_{i} prox_{tau*f_{i}}(x_{i}) 
+    def proximal(self, x, tau, out = None):
         
+        r"""Proximal operator of the BlockFunction at x:
         
-        '''
-
-        if out is not None:
-            if isinstance(tau, Number):
-                for i in range(self.length):
-                    self.functions[i].proximal_conjugate(x.get_item(i), tau, out=out.get_item(i))
-            else:
-                for i in range(self.length):
-                    self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i),out=out.get_item(i))
+            .. math:: \mathrm{prox}_{\tau F}(x) =  (\mathrm{prox}_{\tau f_{i}}(x_{i}))_{i=1}^{m}
             
-        else:
-                
-            out = [None]*self.length
-            if isinstance(tau, Number):
-                for i in range(self.length):
-                    out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau)
-            else:
-                for i in range(self.length):
-                    out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i))
+            Parameter:
             
-            return BlockDataContainer(*out) 
-
-    
-    def proximal(self, x, tau, out = None):
+                x : BlockDataContainer and must have as many rows as self.length            
+        """
         
-        r'''Proximal operator of the convex conjugate of BlockFunction at x:
-        
-            .. math:: prox_{tau*f^{*}}(x) = \sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) 
-        '''
+        if self.length != x.shape[0]:
+            raise ValueError('BlockFunction and BlockDataContainer have incompatible size')        
         
         if out is None:
                 
@@ -135,27 +131,68 @@ class BlockFunction(Function):
             
             
     
-    def gradient(self,x, out=None):
+    def gradient(self, x, out=None):
+        
+        r"""Returns the value of the gradient of the BlockFunction function at x.
         
-        r'''Evaluates gradient of BlockFunction at x
+        .. math:: F'(x) = [f_{1}'(x_{1}), ... , f_{m}'(x_{m})]        
         
-            returns: BlockDataContainer .. math:: [f_{1}'(x_{1}), ... , f_{n}'(x_{n})]
+        Parameter:
+            
+            x : BlockDataContainer and must have as many rows as self.length
                 
-        '''
+        """        
+        
+        if self.length != x.shape[0]:
+            raise ValueError('BlockFunction and BlockDataContainer have incompatible size')        
         
         out = [None]*self.length
         for i in range(self.length):
             out[i] = self.functions[i].gradient(x.get_item(i))
             
-        return  BlockDataContainer(*out)           
+        return  BlockDataContainer(*out)     
+
+    def proximal_conjugate(self, x, tau, out = None):
+        
+        r"""Proximal operator of the convex conjugate of BlockFunction at x:
+        
+            .. math:: \mathrm{prox}_{\tau F^{*}}(x) = (\mathrm{prox}_{\tau f^{*}_{i}}(x^{*}_{i}))_{i=1}^{m}
+            
+            Parameter:
+            
+                x : BlockDataContainer and must have as many rows as self.length            
+        """
+
+        if self.length != x.shape[0]:
+            raise ValueError('BlockFunction and BlockDataContainer have incompatible size')
+
+        if out is not None:
+            if isinstance(tau, Number):
+                for i in range(self.length):
+                    self.functions[i].proximal_conjugate(x.get_item(i), tau, out=out.get_item(i))
+            else:
+                for i in range(self.length):
+                    self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i),out=out.get_item(i))
+            
+        else:
+                
+            out = [None]*self.length
+            if isinstance(tau, Number):
+                for i in range(self.length):
+                    out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau)
+            else:
+                for i in range(self.length):
+                    out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(i))
+            
+            return BlockDataContainer(*out)       
                             
     
     
 if __name__ == '__main__':
     
-    M, N, K = 2,3,5
+    M, N, K = 20,30,50
     
-    from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm
+    from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm, L1Norm
     from ccpi.framework import ImageGeometry, BlockGeometry
     from ccpi.optimisation.operators import Gradient, Identity, BlockOperator
     import numpy
@@ -191,7 +228,22 @@ if __name__ == '__main__':
     numpy.testing.assert_array_almost_equal(res_no_out[1].as_array(), \
                                             res_out[1].as_array(), decimal=4)     
     
-                    
+    
+
+    ig1 = ImageGeometry(M, N)               
+    ig2 = ImageGeometry(2*M, N)
+    ig3 = ImageGeometry(3*M, 4*N)
+    
+    bg = BlockGeometry(ig1,ig2,ig3)
+    
+    z = bg.allocate('random_int')
+    
+    f1 = L1Norm()
+    f2 = 5 * L2NormSquared()
+    
+    f = BlockFunction(f2, f2, f2 + 5, f1 - 4, f1)
+    
+    res = f.convex_conjugate(z)
     
     ##########################################################################
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Function.py b/Wrappers/Python/ccpi/optimisation/functions/Function.py
index 48c6d30..c6c75ca 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/Function.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/Function.py
@@ -23,55 +23,440 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import warnings
-from ccpi.optimisation.functions.ScaledFunction import ScaledFunction
+
+from numbers import Number
+from ccpi.optimisation.operators import ZeroOperator, Identity
 
 class Function(object):
-    '''Abstract class representing a function
     
-    Members:
-    L is the Lipschitz constant of the gradient of the Function 
-    '''
-    def __init__(self):
-        self.L = None
+    """ Abstract class representing a function 
+    
+        :param L: Lipschitz constant of the gradient of the function F(x), when it is differentiable.
+        :param domain: The domain of the function.
 
-    def __call__(self,x, out=None):
-        '''Evaluates the function at x '''
+    """
+    
+    
+    def __init__(self, L = None):
+        
+        self.L = L
+        
+    def __call__(self,x):
+        
+        r"""Returns the value of the function F at x: :math:`F(x)`
+        
+        """        
+        
         raise NotImplementedError
 
     def gradient(self, x, out=None):
-        '''Returns the gradient of the function at x, if the function is differentiable'''
+        
+        r"""Returns the value of the gradient of function F at x, if it is differentiable
+        
+        .. math:: F'(x)
+        
+        """
         raise NotImplementedError
 
     def proximal(self, x, tau, out=None):
-        '''This returns the proximal operator for the function at x, tau'''
+        
+        r"""Returns the proximal operator of function :math:`\tau F` at x        
+        .. math:: \mathrm{prox}_{\tau F}(x) = \underset{z}{\mathrm{argmin}} \frac{1}{2}\|z - x\|^{2} + \tau F(z)
+        """
         raise NotImplementedError
 
-    def convex_conjugate(self, x, out=None):
-        '''This evaluates the convex conjugate of the function at x'''
+    def convex_conjugate(self, x):
+        r""" Returns the convex conjugate of function :math:`F` at :math:`x^{*}`,
+        
+        .. math:: F^{*}(x^{*}) = \underset{x^{*}}{\sup} <x^{*}, x> - F(x)
+                
+        """
         raise NotImplementedError
 
     def proximal_conjugate(self, x, tau, out = None):
-        '''This returns the proximal operator for the convex conjugate of the function at x, tau'''
-        raise NotImplementedError
-
-    def grad(self, x):
-        '''Alias of gradient(x,None)'''
-        warnings.warn('''This method will disappear in following 
-        versions of the CIL. Use gradient instead''', DeprecationWarning)
-        return self.gradient(x, out=None)
+        
+        r"""Returns the proximal operator of the convex conjugate of function :math:`\tau F` at :math:`x^{*}`
+        
+        .. math:: \mathrm{prox}_{\tau F^{*}}(x^{*}) = \underset{z^{*}}{\mathrm{argmin}} \frac{1}{2}\|z^{*} - x^{*}\|^{2} + \tau F^{*}(z^{*})
+        
+        Due to Moreau’s identity, we have an analytic formula to compute the proximal operator of the convex conjugate :math:`F^{*}`
+        
+        .. math:: \mathrm{prox}_{\tau F^{*}}(x) = x - \tau\mathrm{prox}_{\tau^{-1} F}(\tau^{-1}x)
+                
+        """
+        if out is None:
+            return x - tau * self.proximal(x/tau, 1/tau)
+        else:            
+            self.proximal(x/tau, 1/tau, out = out)
+            out*=-tau
+            out.add(x, out = out)                       
+    
+    # Algebra for Function Class
+    
+        # Add functions
+        # Subtract functions
+        # Add/Substract with Scalar
+        # Multiply with Scalar
+    
+    def __add__(self, other):
+        
+        """ Returns the sum of the functions.
+        
+            Cases: a) the sum of two functions :math:`(F_{1}+F_{2})(x) = F_{1}(x) + F_{2}(x)`
+                   b) the sum of a function with a scalar :math:`(F_{1}+scalar)(x) = F_{1}(x) + scalar`
 
-    def prox(self, x, tau):
-        '''Alias of proximal(x, tau, None)'''
-        warnings.warn('''This method will disappear in following 
-        versions of the CIL. Use proximal instead''', DeprecationWarning)
-        return self.proximal(x, tau, out=None)
+        """
+        
+        if isinstance(other, Function):
+            return SumFunction(self, other)
+        elif isinstance(other, (SumFunctionScalar, ConstantFunction, Number)):
+            return SumFunctionScalar(self, other)
+        else:
+            raise ValueError('Not implemented')   
+            
+    def __radd__(self, other):
+        
+        """ Making addition commutative. """
+        return self + other 
+                          
+    def __sub__(self, other):
+        """ Returns the subtraction of the functions."""
+        return self + (-1) * other    
 
     def __rmul__(self, scalar):
-        '''Defines the multiplication by a scalar on the left
-
-        returns a ScaledFunction'''
+        """Returns a function multiplied by a scalar."""               
         return ScaledFunction(self, scalar)
+    
+    def __mul__(self, scalar):
+        return self.__rmul__(scalar)
+    
+   
+    def centered_at(self, center):
+        """ Returns a translated function, namely if we have a function :math:`F(x)` the center is at the origin.         
+            TranslateFunction is :math:`F(x - b)` and the center is at point b."""
+            
+        if center is None:
+            return self
+        else:
+            return TranslateFunction(self, center)  
+    
+class SumFunction(Function):
+    
+    """ SumFunction represents the sum of two functions
+    
+    .. math:: (F_{1} + F_{2})(x)  = F_{1}(x) + F_{2}(x)
+    
+    """
+    
+    def __init__(self, function1, function2 ):
+                
+        super(SumFunction, self).__init__()        
+
+        #if function1.domain != function2.domain:            
+        #    raise ValueError('{} is not the same as {}'.format(function1.domain, function2.domain)) 
+            
+        #self.domain = function1.domain
+                                
+        if function1.L is not None and function2.L is not None:
+            self.L = function1.L + function2.L
+            
+        self.function1 = function1
+        self.function2 = function2               
+            
+    def __call__(self,x):
+        r"""Returns the value of the sum of functions :math:`F_{1}` and :math:`F_{2}` at x
+        
+        .. math:: (F_{1} + F_{2})(x) = F_{1}(x) + F_{2}(x)
+                
+        """  
+        return self.function1(x) + self.function2(x)
+    
+    def gradient(self, x, out=None):
+        
+        r"""Returns the value of the sum of the gradient of functions :math:`F_{1}` and :math:`F_{2}` at x, 
+        if both of them are differentiable
+        
+        .. math:: (F'_{1} + F'_{2})(x)  = F'_{1}(x) + F'_{2}(x)
+        
+        """
+        
+#        try: 
+        if out is None:            
+            return self.function1.gradient(x) +  self.function2.gradient(x)  
+        else:
+
+            self.function1.gradient(x, out=out)
+            out.add(self.function2.gradient(x), out=out)                
+            
+class ScaledFunction(Function):
+    
+    r""" ScaledFunction represents the scalar multiplication with a Function.
+    
+    Let a function F then and a scalar :math:`\alpha`.
+    
+    If :math:`G(x) = \alpha F(x)` then:
+    
+    1. :math:`G(x) = \alpha  F(x)` ( __call__ method )
+    2. :math:`G'(x) = \alpha  F'(x)` ( gradient method ) 
+    3. :math:`G^{*}(x^{*}) = \alpha  F^{*}(\frac{x^{*}}{\alpha})` ( convex_conjugate method )   
+    4. :math:`\mathrm{prox}_{\tau G}(x) = \mathrm{prox}_{(\tau\alpha) F}(x)` ( proximal method ) 
+           
+    """
+    def __init__(self, function, scalar):
+        
+        super(ScaledFunction, self).__init__() 
+                                                     
+        if not isinstance (scalar, Number):
+            raise TypeError('expected scalar: got {}'.format(type(scalar)))
+        
+        if function.L is not None:        
+            self.L = abs(scalar) * function.L  
+            
+        self.scalar = scalar
+        self.function = function       
+              
+    def __call__(self,x, out=None):
+        r"""Returns the value of the scaled function.
+        
+        .. math:: G(x) = \alpha F(x)
+        
+        """
+        return self.scalar * self.function(x)
+
+    def convex_conjugate(self, x):
+        r"""Returns the convex conjugate of the scaled function.
+        
+        .. math:: G^{*}(x^{*}) = \alpha  F^{*}(\frac{x^{*}}{\alpha})
+        
+        """
+        return self.scalar * self.function.convex_conjugate(x/self.scalar)
+    
+    def gradient(self, x, out=None):
+        r"""Returns the gradient of the scaled function.
+        
+        .. math:: G'(x) = \alpha  F'(x)
+        
+        """
+        
+#        try:
+        if out is None:            
+            return self.scalar * self.function.gradient(x)
+        else:
+            self.function.gradient(x, out=out)
+            out *= self.scalar  
+
+    def proximal(self, x, tau, out=None):
+        
+        r"""Returns the proximal operator of the scaled function.
+        
+        .. math:: \mathrm{prox}_{\tau G}(x) = \mathrm{prox}_{(\tau\alpha) F}(x)
+        
+        """        
+
+        return self.function.proximal(x, tau*self.scalar, out=out)     
+
+
+    def proximal_conjugate(self, x, tau, out = None):
+        r"""This returns the proximal operator for the function at x, tau
+        """
+        if out is None:
+            return self.scalar * self.function.proximal_conjugate(x/self.scalar, tau/self.scalar)
+        else:
+            self.function.proximal_conjugate(x/self.scalar, tau/self.scalar, out=out)
+            out *= self.scalar
+
+class SumFunctionScalar(SumFunction):
+          
+    """ SumFunctionScalar represents the sum a function with a scalar. 
+    
+        .. math:: (F + scalar)(x)  = F(x) + scalar
+    
+        Although SumFunction has no general expressions for 
+        
+        i) convex_conjugate
+        ii) proximal
+        iii) proximal_conjugate
+            
+        if the second argument is a ConstantFunction then we can derive the above analytically.
+    
+    """    
+    
+    def __init__(self, function, constant):
+        
+        super(SumFunctionScalar, self).__init__(function, ConstantFunction(constant))        
+        self.constant = constant
+        self.function = function
+        
+    def convex_conjugate(self,x):
+        
+        r""" Returns the convex conjugate of a :math:`(F+scalar)`
+                
+        .. math:: (F+scalar)^{*}(x^{*}) = F^{*}(x^{*}) - scalar
+                        
+        """        
+        return self.function.convex_conjugate(x) - self.constant
+    
+    def proximal(self, x, tau, out=None):
+        
+        """ Returns the proximal operator of :math:`F+scalar`
+
+        .. math:: \mathrm{prox}_{\tau (F+scalar)}(x) = \mathrm{prox}_{\tau F}
+                        
+        """                
+        return self.function.proximal(x, tau, out=out)        
+            
+#    def function(self):       
+#       return self.function    
+
+class ConstantFunction(Function):
+    
+            
+    r""" ConstantFunction: :math:`F(x) = constant, constant\in\mathbb{R}`         
+        
+    """
+    
+    def __init__(self, constant = 0):
+        
+        super(ConstantFunction, self).__init__(L=0)
+        
+        if not isinstance (constant, Number):
+            raise TypeError('expected scalar: got {}'.format(type(constant)))
+                
+        self.constant = constant
+              
+    def __call__(self,x):
+        
+        """ Returns the value of the function, :math:`F(x) = constant`"""
+        return self.constant
+        
+    def gradient(self, x, out=None):
+        
+        """ Returns the value of the gradient of the function, :math:`F'(x)=0`"""       
+        if out is None:
+            return x * 0.
+        else:
+            out.fill(0)
+    
+    def convex_conjugate(self, x):
+                        
+        r""" The convex conjugate of constant function :math:`F(x) = c\in\mathbb{R}` is
+        
+        .. math:: 
+            F(x^{*}) 
+            =
+            \begin{cases}
+                -c, & if x^{*} = 0\\
+                \infty, & \mbox{otherwise}
+            \end{cases}
+                                                          
+                    
+        However, :math:`x^{*} = 0` only in the limit of iterations, so in fact this can be infinity.
+        We do not want to have inf values in the convex conjugate, so we have to penalise this value accordingly.
+        The following penalisation is useful in the PDHG algorithm, when we compute primal & dual objectives
+        for convergence purposes.
+        
+        .. math:: F^{*}(x^{*}) = \sum \max\{x^{*}, 0\}
+        
+        """               
+        return x.maximum(0).sum()
+                
+    def proximal(self, x, tau, out=None):
+        
+        """Returns the proximal operator of the constant function, which is the same element, i.e.,
+        
+        .. math:: \mathrm{prox}_{\tau F}(x) = x 
+        
+        """
+        if out is None:
+            return x.copy()
+        else:
+            out.fill(x)
+        #return Identity(x.geometry).direct(x, out=out)
 
+class ZeroFunction(ConstantFunction):
+    
+    """ ZeroFunction represents the zero function, :math:`F(x) = 0`        
+    """
+    
+    def __init__(self):
+        super(ZeroFunction, self).__init__(constant = 0.) 
+        
+class TranslateFunction(Function):
+    
+    r""" TranslateFunction represents the translation of function F with respect to the center b.
+    
+    Let a function F and consider :math:`G(x) = F(x - center)`. 
+    
+    Function F is centered at 0, whereas G is centered at point b.
+    
+    If :math:`G(x) = F(x - b)` then:
+    
+    1. :math:`G(x) = F(x - b)` ( __call__ method )
+    2. :math:`G'(x) = F'(x - b)` ( gradient method ) 
+    3. :math:`G^{*}(x^{*}) = F^{*}(x^{*}) + <x^{*}, b >` ( convex_conjugate method )   
+    4. :math:`\mathrm{prox}_{\tau G}(x) = \mathrm{prox}_{\tau F}(x - b)  + b` ( proximal method ) 
+               
+    """
+    
+    def __init__(self, function, center):
+        
+        super(TranslateFunction, self).__init__(L = function.L) 
+                        
+        self.function = function
+        self.center = center
+        
+        
+                
+    def __call__(self, x):
+        
+        r"""Returns the value of the translated function.
+        
+        .. math:: G(x) = F(x - b)
+        
+        """        
+        
+        return self.function(x - self.center)
+    
+    def gradient(self, x, out = None):
+        
+        r"""Returns the gradient of the translated function.
+        
+        .. math:: G'(x) =  F'(x - b)
+        
+        """        
+        
+        if out is None:
+            return self.function.gradient(x - self.center)
+        else:                       
+            x.subtract(self.center, out = out)
+            self.function.gradient(out, out = out)           
+    
+    def proximal(self, x, tau, out = None):
+        
+        r"""Returns the proximal operator of the translated function.
+        
+        .. math:: \mathrm{prox}_{\tau G}(x) = \mathrm{prox}_{\tau F}(x-b) + b
+        
+        """        
+        
+        if out is None:
+            return self.function.proximal(x - self.center, tau) + self.center
+        else:                    
+            x.subtract(self.center, out = out)
+            self.function.proximal(out, tau, out = out)
+            out.add(self.center, out = out)
+                    
+    def convex_conjugate(self, x):
+        
+        r"""Returns the convex conjugate of the translated function.
+        
+        .. math:: G^{*}(x^{*}) = F^{*}(x^{*}) + <x^{*}, b >
+        
+        """        
+        
+        return self.function.convex_conjugate(x) + self.center.dot(x)
+                           
+    
+    
diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
index ed5c1b1..51105c8 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
@@ -20,20 +20,26 @@
 #
 #=========================================================================
 
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions import ScaledFunction
+from ccpi.optimisation.operators import Operator, ScaledOperator
 
+import warnings
 
 class FunctionOperatorComposition(Function):
     
-    r'''Function composition with Operator: :math:`(f \otimes A)(x) = f(Ax)`
+    """ Composition of a function with an operator as : :math:`(F \otimes A)(x) = F(Ax)`
     
-    :param A: operator
-    :type A: :code:`Operator`
-    :param f: function
-    :type f: :code:`Function`
+            :parameter function: :code:`Function` F
+            :parameter operator: :code:`Operator` A
+            
+            
+        For general operator, we have no explicit formulas for convex_conjugate,
+        proximal and proximal_conjugate            
     
-    '''
+    """
     
     def __init__(self, function, operator):
         '''creator
@@ -48,91 +54,40 @@ class FunctionOperatorComposition(Function):
         
         self.function = function     
         self.operator = operator
+        
+        if not isinstance(self.function, Function):
+            raise ValueError('{} is not function '.format(type(self.function)))
+          
+        # check also ScaledOperator because it's not a child atm of Operator            
+        if not isinstance(self.operator, (Operator, ScaledOperator)):
+            raise ValueError('{} is not function '.format(type(self.operator)))            
+        
         try:
             self.L = function.L * operator.norm()**2 
-        except Error as er:
+        except ValueError as er:
             self.L = None
             warnings.warn("Lipschitz constant was not calculated")
         
     def __call__(self, x):
         
-        '''Evaluates f(Ax)'''
+        """ Returns :math:`F(Ax)`
+        """
     
         return self.function(self.operator.direct(x))  
     
     def gradient(self, x, out=None):
         
-        '''Evaluates gradient of f(Ax):
+        """ Return the gradient of F(Ax), 
         
-        ..math ::  A^{T}f'(Ax)
+        ..math ::  (F(Ax))' = A^{T}F'(Ax)
             
-        '''
+        """
         
         tmp = self.operator.range_geometry().allocate()
         self.operator.direct(x, out=tmp)
         self.function.gradient(tmp, out=tmp)
         if out is None:
-            #return self.operator.adjoint(self.function.gradient(self.operator.direct(x)))
             return self.operator.adjoint(tmp)
         else: 
             self.operator.adjoint(tmp, out=out)
 
-    
-
-                
-if __name__ == '__main__':   
-
-    from ccpi.framework import ImageGeometry, AcquisitionGeometry
-    from ccpi.optimisation.operators import Gradient
-    from ccpi.optimisation.functions import L2NormSquared
-    from ccpi.astra.ops import AstraProjectorSimple
-    import numpy as np
-        
-    M, N= 50, 50
-    ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
-    
-    detectors = N
-    angles_num = N    
-    det_w = 1.0
-    
-    angles = np.linspace(0, np.pi, angles_num, endpoint=False)
-    ag = AcquisitionGeometry('parallel',
-                             '2D',
-                             angles,
-                             detectors,det_w)
-    
-    
-    Aop = AstraProjectorSimple(ig, ag, 'cpu')    
-
-    u = ig.allocate('random_int', seed=15)
-    u1 = ig.allocate('random_int', seed=10)
-    b = Aop.direct(u1)
-    
-        
-#    G = Gradient(ig)
-    alpha = 0.5
-    
-    f1 =  alpha * L2NormSquared(b=b)    
-
-    f_comp = FunctionOperatorComposition(f1, Aop)
-    
-    print(f_comp(u))
-    
-    
-    z1 = Aop.direct(u)
-    tmp = 0.5 * ((z1 - b)**2).sum()
-    
-   
-    print(tmp)
-    
-    
-    
-    
-    
-    
-          
-    
-
-             
-        
-                       
diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
index 9e9e55c..4383858 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
@@ -19,7 +19,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.functions import Function
 import numpy
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index 4239423..39bd983 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -20,43 +20,76 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions.ScaledFunction import ScaledFunction 
+
 import functools
 import scipy.special
 
 class KullbackLeibler(Function):
     
-    r'''Kullback-Leibler divergence function
-    
-         .. math::
-              f(x, y) = \begin{cases} x \log(x / y) - x + y & x > 0, y > 0 \\ 
-                                    y & x = 0, y \ge 0 \\
-                                    \infty & \text{otherwise} 
-                       \end{cases}
+    r""" Kullback Leibler divergence function is defined as:
             
-    '''
+    .. math:: F(u, v)
+            = \begin{cases} 
+            u \log(\frac{u}{v}) - u + v & \mbox{ if } u > 0, v > 0\\
+            v & \mbox{ if } u = 0, v \ge 0 \\
+            \infty, & \mbox{otherwise}
+            \end{cases}  
+            
+    where we use the :math:`0\log0 := 0` convention. 
+
+    At the moment, we use build-in implemention of scipy, see
+    https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.kl_div.html
+    
+    The Kullback-Leibler function is used as a fidelity term in minimisation problems where the
+    acquired data follow Poisson distribution. If we denote the acquired data with :math:`b`
+    then, we write
+    
+     .. math:: \underset{i}{\sum} F(b_{i}, (v + \eta)_{i})
+     
+     where, :math:`\eta` is an additional noise. 
+     
+     Example: In the case of Positron Emission Tomography reconstruction :math:`\eta` represents 
+     scatter and random events contribution during the PET acquisition. Hence, in that case the KullbackLeibler
+     fidelity measures the distance between :math:`\mathcal{A}v + \eta` and acquisition data :math:`b`, where
+     :math:`\mathcal{A}` is the projection operator.
+     
+     This is related to PoissonLogLikelihoodWithLinearModelForMean definition that is used in PET reconstruction
+     in the PET-MR software , see https://github.com/CCPPETMR and for more details in
+    
+    http://stir.sourceforge.net/documentation/doxy/html/classstir_1_1PoissonLogLikelihoodWithLinearModelForMean.html
+                        
+    """            
+    
     
-    def __init__(self, data, **kwargs):
+    def __init__(self,  **kwargs):
         
-        super(KullbackLeibler, self).__init__()
+        super(KullbackLeibler, self).__init__(L = None)          
+        self.b = kwargs.get('b', None)
         
-        self.b = data    
-        self.bnoise = data * 0.
+        if self.b is None:
+            raise ValueError('Please define data, as b = ...')
+            
+        if (self.b.as_array() < 0).any():            
+            raise ValueError('Data should be larger or equal to 0')              
+         
+        self.eta = kwargs.get('eta',self.b * 0.0)
         
                                                     
     def __call__(self, x):
         
 
-        '''Evaluates KullbackLeibler at x'''
-
-        ind = x.as_array()>0
-        tmp = scipy.special.kl_div(self.b.as_array()[ind], x.as_array()[ind])                
-        return numpy.sum(tmp) 
-
+        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)         
+        
     def log(self, datacontainer):
         '''calculates the in-place log of the datacontainer'''
         if not functools.reduce(lambda x,y: x and y>0, datacontainer.as_array().ravel(), True):
@@ -66,52 +99,62 @@ class KullbackLeibler(Function):
         
     def gradient(self, x, out=None):
         
-        '''Evaluates gradient of KullbackLeibler at x'''
+        r"""Returns the value of the gradient of the KullbackLeibler function at :math:`(b, x + \eta)`.                
         
-        if out is None:
-            return 1 - self.b/(x + self.bnoise)
-        else:
+        .. math:: F'(b, x + \eta) = 1 - \frac{b}{x+\eta}
+        
+        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]                                
 
-            x.add(self.bnoise, out=out)
-            self.b.divide(out, out=out)
-            out.subtract(1, out=out)
-            out *= -1
             
     def convex_conjugate(self, x):
         
-        '''Convex conjugate of KullbackLeibler at x'''
+        r"""Returns the value of the convex conjugate of the KullbackLeibler function at :math:`(b, x + \eta)`.                
         
-        xlogy = - scipy.special.xlogy(self.b.as_array(), 1 - x.as_array())
-        return numpy.sum(xlogy)
+        .. math:: F^{*}(b, x + \eta) = - b \log(1-x^{*}) - <x^{*}, \eta> 
+        
+        """  
+        tmp = 1 - x.as_array()
+        ind = tmp>0
+        xlogy = - scipy.special.xlogy(self.b.as_array()[ind], tmp[ind])  
+        return numpy.sum(xlogy) - self.eta.dot(x)
             
     def proximal(self, x, tau, out=None):
         
-        r'''Proximal operator of KullbackLeibler at x
-           
-           .. math::     prox_{\tau * f}(x)
-
-        '''
-
+        r"""Returns the value of the proximal operator of the KullbackLeibler function at :math:`(b, x + \eta)`.
+        
+        .. math:: \mathrm{prox}_{\tau F}(x) = \frac{1}{2}\bigg( (x - \eta - \tau) + \sqrt{ (x + \eta - \tau)^2 + 4\tau b} \bigg)
+        
+        The proximal for the convex conjugate of :math:`F` is 
+        
+        .. math:: \mathrm{prox}_{\tau F^{*}}(x) = 0.5*((z + 1) - \sqrt{(z-1)^2 + 4 * \tau b})
+        
+        where :math:`z = x + \tau \eta`
+                    
+        """
+          
         if out is None:        
-            return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt() )
-        else:
-            
-            tmp =  0.5 *( (x - self.bnoise - tau) + 
-                        ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt()
-                        )
-            x.add(self.bnoise, out=out)
+            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
-            tmp = self.b * (4 * tau)
-            out.add(tmp, out=out)
-            out.sqrt(out=out)
-            
-            x.subtract(self.bnoise, out=tmp)
-            tmp -= tau
-            
-            out += tmp
-            
-            out *= 0.5
+            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):
         
@@ -122,11 +165,10 @@ class KullbackLeibler(Function):
 
                 
         if out is None:
-            z = x + tau * self.bnoise
+            z = x + tau * self.eta
             return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
-        else:
-            
-            tmp = tau * self.bnoise
+        else:            
+            tmp = tau * self.eta
             tmp += x
             tmp -= 1
             
@@ -139,105 +181,93 @@ class KullbackLeibler(Function):
             out += tmp
             out *= 0.5
 
-    def __rmul__(self, scalar):
-        
-        '''Multiplication of KullbackLeibler with a scalar        
-            
-            Returns: ScaledFunction
-        '''
-        
-        return ScaledFunction(self, scalar) 
+
 
 
 if __name__ == '__main__':
     
     from ccpi.framework import ImageGeometry
-    import numpy
-    
-    M, N =  2,3
-    ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
-    u = ig.allocate('random_int')
-    b = ig.allocate('random_int')
-    u.as_array()[1,1]=0
-    u.as_array()[2,0]=0
-    b.as_array()[1,1]=0
-    b.as_array()[2,0]=0    
+    import numpy as np
     
-    f = KullbackLeibler(b)
+    M, N, K =  30, 30, 20
+    ig = ImageGeometry(N, M, K)
     
+    u1 = ig.allocate('random_int', seed = 500)    
+    g1 = ig.allocate('random_int', seed = 1000)
+    b1 = ig.allocate('random', seed = 5000)
     
-#    longest = reduce(lambda x, y: len(x) if len(x) > len(y) else len(y), strings)
-
-
-#    tmp = functools.reduce(lambda x, y: \
-#                           0 if x==0 and not numpy.isnan(y) else x * numpy.log(y), \
-#                           zip(b.as_array().ravel(), u.as_array().ravel()),0)
-    
-    
-#    np.multiply.reduce(X, 0)
-    
-    
-#                sf = reduce(lambda x,y: x + y[0]*y[1],
-#                            zip(self.as_array().ravel(),
-#                                other.as_array().ravel()),
-#                            0)        
-#cdef inline number_t xlogy(number_t x, number_t y) nogil:
-#    if x == 0 and not zisnan(y):
-#        return 0
-#    else:
-#        return x * zlog(y)        
-    
-#    if npy_isnan(x):
-#        return x
-#    elif x > 0:
-#        return -x * log(x)
-#    elif x == 0:
-#        return 0
-#    else:
-#        return -inf    
-    
-#        cdef inline double kl_div(double x, double y) nogil:
-#    if npy_isnan(x) or npy_isnan(y):
-#        return nan
-#    elif x > 0 and y > 0:
-#        return x * log(x / y) - x + y
-#    elif x == 0 and y >= 0:
-#        return y
-#    else:
-#        return inf    
-
-    
-    
-    
-#    def xlogy(self, dc1, dc2):
+    # with no data
+    try:
+        f = KullbackLeibler()   
+    except ValueError:
+        print('Give data b=...\n')
         
-#        return numpy.sum(numpy.where(dc1.as_array() != 0, dc2.as_array() * numpy.log(dc2.as_array() / dc1.as_array()), 0))
+    print('With negative data, no background\n')   
+    try:        
+        f = KullbackLeibler(b=-1*g1)
+    except ValueError:
+        print('We have negative data\n') 
         
-           
-    
-#    f.xlog(u, b)
-    
+    f = KullbackLeibler(b=g1)        
+        
+    print('Check KullbackLeibler(x,x)=0\n') 
+    numpy.testing.assert_equal(0.0, f(g1))
             
-
-    
-#    tmp1 = b.as_array()
-#    tmp2 = u.as_array()
-#    
-#    zz = scipy.special.xlogy(tmp1, tmp2)
-#
-#    print(np.sum(zz))
-    
-    
-#    ww = f.xlogy(b, u)
-    
-#    print(ww)
+    print('Check gradient .... is OK \n')
+    res_gradient = f.gradient(u1)
+    res_gradient_out = u1.geometry.allocate()
+    f.gradient(u1, out = res_gradient_out) 
+    numpy.testing.assert_array_almost_equal(res_gradient.as_array(), \
+                                            res_gradient_out.as_array(),decimal = 4)  
     
+    print('Check proximal ... is OK\n')        
+    tau = 0.4
+    res_proximal = f.proximal(u1, tau)
+    res_proximal_out = u1.geometry.allocate()   
+    f.proximal(u1, tau, out = res_proximal_out)
+    numpy.testing.assert_array_almost_equal(res_proximal.as_array(), \
+                                            res_proximal_out.as_array(), decimal =5)  
     
-#cdef inline double kl_div(double x, double y) nogil:
-  
-    
+    print('Check conjugate ... is OK\n')  
+    res_conj = f.convex_conjugate(u1) 
     
+    if (1 - u1.as_array()).all():
+        print('If 1-x<=0, Convex conjugate returns 0.0')
         
+    numpy.testing.assert_equal(0.0, f.convex_conjugate(u1))   
+
 
+    print('Check KullbackLeibler with background\n')       
+    f1 = KullbackLeibler(b=g1, eta=b1) 
+        
+    tmp_sum = (u1 + f1.eta).as_array()
+    ind = tmp_sum >= 0
+    tmp = scipy.special.kl_div(f1.b.as_array()[ind], tmp_sum[ind])                 
+    numpy.testing.assert_equal(f1(u1), numpy.sum(tmp) )
+    
+    print('Check proximal KL without background\n')   
+    tau = [0.1, 1, 10, 100, 10000]
     
+    for t1 in tau:
+        
+        proxc = f.proximal_conjugate(u1,t1)
+        proxc_out = ig.allocate()
+        f.proximal_conjugate(u1, t1, out = proxc_out)
+        print('tau = {} is OK'.format(t1) )
+        numpy.testing.assert_array_almost_equal(proxc.as_array(), 
+                                                proxc_out.as_array(),
+                                                decimal = 4)
+        
+    print('\nCheck proximal KL with background\n')          
+    for t1 in tau:
+        
+        proxc1 = f1.proximal_conjugate(u1,t1)
+        proxc_out1 = ig.allocate()
+        f1.proximal_conjugate(u1, t1, out = proxc_out1)
+        numpy.testing.assert_array_almost_equal(proxc1.as_array(), 
+                                                proxc_out1.as_array(),
+                                                decimal = 4)  
     
+        print('tau = {} is OK'.format(t1) )        
+        
+    
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
index 1fcfcca..e8fa008 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
@@ -20,22 +20,21 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
-from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions.ScaledFunction import ScaledFunction        
+from ccpi.optimisation.functions import Function       
 from ccpi.optimisation.operators import ShrinkageOperator 
+import numpy as np
  
 
 class L1Norm(Function):
     
-    r'''L1Norm function: 
+    r"""L1Norm function
             
-            Cases considered (with/without data):            
-                a) :math:`f(x) = ||x||_{1}`
-                b) :math:`f(x) = ||x - b||_{1}`
+        Consider the following cases:           
+            a) .. math:: F(x) = ||x||_{1}
+            b) .. math:: F(x) = ||x - b||_{1}
                                 
-    '''   
+    """   
            
     def __init__(self, **kwargs):
         '''creator
@@ -53,34 +52,62 @@ class L1Norm(Function):
         
     def __call__(self, x):
         
-        '''Evaluates L1Norm at x'''
+        r"""Returns the value of the L1Norm function at x.
+        
+        Consider the following cases:           
+            a) .. math:: F(x) = ||x||_{1}
+            b) .. math:: F(x) = ||x - b||_{1}        
+        
+        """
         
         y = x
         if self.b is not None: 
             y = x - self.b
         return y.abs().sum()  
-    
-    def gradient(self,x):
-        
-        return ValueError('Not Differentiable')   
-    
+          
     def convex_conjugate(self,x):
         
-        '''Convex conjugate of L1Norm at x'''
-
-        y = 0        
-        if self.b is not None:
-            y =  0 + self.b.dot(x)
-        return y  
+        r"""Returns the value of the convex conjugate of the L1Norm function at x.
+        Here, we need to use the convex conjugate of L1Norm, which is the Indicator of the unit 
+        :math:`L^{\infty}` norm
+        
+        Consider the following cases:
+                
+                a) .. math:: F^{*}(x^{*}) = \mathbb{I}_{\{\|\cdot\|_{\infty}\leq1\}}(x^{*}) 
+                b) .. math:: F^{*}(x^{*}) = \mathbb{I}_{\{\|\cdot\|_{\infty}\leq1\}}(x^{*}) + <x^{*},b>      
+        
     
+        .. math:: \mathbb{I}_{\{\|\cdot\|_{\infty}\leq1\}}(x^{*}) 
+            = \begin{cases} 
+            0, \mbox{if } \|x^{*}\|_{\infty}\leq1\\
+            \infty, \mbox{otherwise}
+            \end{cases}
+    
+        """        
+        tmp = (np.abs(x.as_array()).max() - 1)
+        if tmp<=1e-5:            
+            if self.b is not None:
+                return self.b.dot(x)
+            else:
+                return 0.
+        return np.inf        
+                    
     def proximal(self, x, tau, out=None):
         
-        r'''Proximal operator of L1Norm at x
-           
-           ..math::     prox_{\tau * f}(x)
+        r"""Returns the value of the proximal operator of the L1Norm function at x.
+        
+        
+        Consider the following cases:
                 
-        ''' 
+                a) .. math:: \mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x)
+                b) .. math:: \mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x) + b   
+    
+        where,
         
+        .. math :: \mathrm{prox}_{\tau F}(x) = \mathrm{ShinkOperator}(x) = sgn(x) * \max\{ |x| - \tau, 0 \}
+                            
+        """  
+            
         if out is None:
             if self.b is not None:
                 return self.b + self.shinkage_operator(x - self.b, tau)
@@ -92,33 +119,33 @@ class L1Norm(Function):
             else:
                 out.fill(self.shinkage_operator(x, tau))
                                     
-    def proximal_conjugate(self, x, tau, out=None):
-        
-        r'''Proximal operator of the convex conjugate of L1Norm at x:
-                
-            .. math:: prox_{\tau * f^{*}}(x)
-                
-        '''          
-        
-        if out is None:
-            if self.b is not None:
-                return (x - tau*self.b).divide((x - tau*self.b).abs().maximum(1.0))
-            else:
-                return x.divide(x.abs().maximum(1.0))
-        else:
-            if self.b is not None:
-                out.fill((x - tau*self.b).divide((x - tau*self.b).abs().maximum(1.0)))
-            else:
-                out.fill(x.divide(x.abs().maximum(1.0)) )                
-            
-    def __rmul__(self, scalar):
-        
-        '''Multiplication of L2NormSquared with a scalar        
+#    def proximal_conjugate(self, x, tau, out=None):
+#        
+#        r'''Proximal operator of the convex conjugate of L1Norm at x:
+#                
+#            .. math:: prox_{\tau * f^{*}}(x)
+#                
+#        '''          
+#        
+#        if out is None:
+#            if self.b is not None:
+#                return (x - tau*self.b).divide((x - tau*self.b).abs().maximum(1.0))
+#            else:
+#                return x.divide(x.abs().maximum(1.0))
+#        else:
+#            if self.b is not None:
+#                out.fill((x - tau*self.b).divide((x - tau*self.b).abs().maximum(1.0)))
+#            else:
+#                out.fill(x.divide(x.abs().maximum(1.0)) )                
             
-            Returns: ScaledFunction
-        '''
-        
-        return ScaledFunction(self, scalar)
+#    def __rmul__(self, scalar):
+#        
+#        '''Multiplication of L2NormSquared with a scalar        
+#            
+#            Returns: ScaledFunction
+#        '''
+#        
+#        return ScaledFunction(self, scalar)
 
 
 if __name__ == '__main__':   
@@ -185,4 +212,4 @@ if __name__ == '__main__':
         
         
         
-      
\ No newline at end of file
+      
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
index ef7c698..613a7f5 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
@@ -21,20 +21,28 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions.ScaledFunction import ScaledFunction
 
 class L2NormSquared(Function):
     
-    r'''L2NormSquared function: 
-            
-            Cases considered (with/without data):            
-                a) .. math:: f(x) = \|x\|^{2}_{2} 
-                b) .. math:: f(x) = \|\|x - b\|\|^{2}_{2}
-                                
-    '''    
+    r""" L2NormSquared function: :math:`F(x) = \| x\|^{2}_{2} = \underset{i}{\sum}x_{i}^{2}`
+          
+    Following cases are considered:
+                
+        a) :math:`F(x) = \|x\|^{2}_{2}`
+        b) :math:`F(x) = \|x - b\|^{2}_{2}`
+        
+    .. note::  For case b) case we can use :code:`F = L2NormSquared().centered_at(b)`,
+               see *TranslateFunction*.
+        
+    :Example:
+        
+        >>> F = L2NormSquared()
+        >>> F = L2NormSquared(b=b) 
+        >>> F = L2NormSquared().centered_at(b)
+                                                          
+    """    
     
     def __init__(self, **kwargs):
         '''creator
@@ -46,15 +54,25 @@ class L2NormSquared(Function):
         :param b:  translation of the function
         :type b: :code:`DataContainer`, optional
         '''                        
-        super(L2NormSquared, self).__init__()
+        super(L2NormSquared, self).__init__(L = 2)
         self.b = kwargs.get('b',None) 
         
-        # Lipschitz constant
-        self.L = 2
-
+        #if self.b is not None:
+        #    self.domain = self.b.geometry  
+                            
     def __call__(self, x):
+
+        r"""Returns the value of the L2NormSquared function at x.
         
-        '''Evaluates L2NormSquared at x'''
+        Following cases are considered:
+            
+            a) :math:`F(x) = \|x\|^{2}_{2}`
+            b) :math:`F(x) = \|x - b\|^{2}_{2}`
+    
+        :param: :math:`x`
+        :returns: :math:`\underset{i}{\sum}x_{i}^{2}`
+                
+        """          
             
         y = x
         if self.b is not None: 
@@ -67,7 +85,14 @@ class L2NormSquared(Function):
                 
     def gradient(self, x, out=None):        
         
-        '''Evaluates gradient of L2NormSquared at x'''
+        r"""Returns the value of the gradient of the L2NormSquared function at x.
+        
+        Following cases are considered:
+                
+            a) :math:`F'(x) = 2x`
+            b) :math:`F'(x) = 2(x-b)`
+                
+        """
                 
         if out is not None:
             
@@ -86,23 +111,34 @@ class L2NormSquared(Function):
                                                        
     def convex_conjugate(self, x):
         
-        '''Convex conjugate of L2NormSquared at x'''
+        r"""Returns the value of the convex conjugate of the L2NormSquared function at x.
         
+        Consider the following cases:
+                
+                a) .. math:: F^{*}(x^{*}) = \frac{1}{4}\|x^{*}\|^{2}_{2} 
+                b) .. math:: F^{*}(x^{*}) = \frac{1}{4}\|x^{*}\|^{2}_{2} + <x^{*}, b>
+                
+        """                
         tmp = 0
         
         if self.b is not None:
-            tmp = x.dot(self.b) #(x * self.b).sum()
+            tmp = x.dot(self.b) 
             
         return (1./4.) * x.squared_norm() + tmp
 
 
     def proximal(self, x, tau, out = None):
-
-        r'''Proximal operator of L2NormSquared at x
-            
-            .. math:: prox_{\tau * f}(x)
-        '''          
         
+        r"""Returns the value of the proximal operator of the L2NormSquared function at x.
+        
+        
+        Consider the following cases:
+                
+                a) .. math:: \mathrm{prox}_{\tau F}(x) = \frac{x}{1+2\tau}
+                b) .. math:: \mathrm{prox}_{\tau F}(x) = \frac{x-b}{1+2\tau} + b      
+                        
+        """            
+
         if out is None:
             
             if self.b is None:
@@ -122,33 +158,23 @@ class L2NormSquared(Function):
                 x.divide((1+2*tau), out=out)
 
     
-    def proximal_conjugate(self, x, tau, out=None):
-        
-        r'''Proximal operator of the convex conjugate of L2NormSquared at x:
-           
-           .. math::  prox_{\tau * f^{*}}(x)'''
-        
-        if out is None:
-            if self.b is not None:
-                return (x - tau*self.b)/(1 + tau/2) 
-            else:
-                return x/(1 + tau/2)
-        else:
-            if self.b is not None:
-                x.subtract(tau*self.b, out=out)
-                out.divide(1+tau/2, out=out)
-            else:
-                x.divide(1 + tau/2, out=out)
-                                        
-    def __rmul__(self, scalar):
-        
-        '''Multiplication of L2NormSquared with a scalar        
-           
-           Returns: ScaledFunction'''
-        
-        return ScaledFunction(self, scalar)  
-
-
+#    def proximal_conjugate(self, x, tau, out=None):
+#        
+#        r'''Proximal operator of the convex conjugate of L2NormSquared at x:
+#           
+#           .. math::  prox_{\tau * f^{*}}(x)'''
+#        
+#        if out is None:
+#            if self.b is not None:
+#                return (x - tau*self.b)/(1 + tau/2) 
+#            else:
+#                return x/(1 + tau/2)
+#        else:
+#            if self.b is not None:
+#                x.subtract(tau*self.b, out=out)
+#                out.divide(1+tau/2, out=out)
+#            else:
+#                x.divide(1 + tau/2, out=out)                                        
 
 if __name__ == '__main__':
     
@@ -156,7 +182,7 @@ if __name__ == '__main__':
     import numpy
     # TESTS for L2 and scalar * L2
     
-    M, N, K = 20,30,50
+    M, N, K = 2,3,1
     ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K)
     u = ig.allocate('random_int')
     b = ig.allocate('random_int') 
@@ -169,12 +195,14 @@ if __name__ == '__main__':
     numpy.testing.assert_equal(f(u), u.squared_norm())
 
     # check grad/call with data
+    
+    igggg = ImageGeometry(4,4)
     f1 = L2NormSquared(b=b)
     b1 = f1.gradient(u)
     b2 = 2 * (u-b)
         
     numpy.testing.assert_array_almost_equal(b1.as_array(), b2.as_array(), decimal=4)
-    numpy.testing.assert_equal(f1(u), (u-b).squared_norm())
+    numpy.testing.assert_equal(f1(u), ((u-b)).squared_norm())
     
     #check convex conjuagate no data
     c1 = f.convex_conjugate(u)
@@ -288,4 +316,4 @@ if __name__ == '__main__':
     
     numpy.testing.assert_array_almost_equal(res1.as_array(), \
                                             res2.as_array(), decimal=4)
-                                            
\ No newline at end of file
+                                            
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/LeastSquares.py
index 01c4f38..a0baedc 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/LeastSquares.py
@@ -20,30 +20,34 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.operators import LinearOperator
 from ccpi.optimisation.functions import Function
 import warnings
 
 # Define a class for squared 2-norm
-class Norm2Sq(Function):
-    r'''.. math:: f(x) = c*||A*x-b||_2^2
+class LeastSquares(Function):
+    r"""Least Squares function
     
-    which has 
+    .. math:: F(x) = c\|Ax-b\|_2^2
     
-    .. math:: grad[f](x) = 2*c*A^T*(A*x-b)
-    
-    and Lipschitz constant
-    
-    .. math:: L = 2*c*||A||_2^2 = 2*s1(A)^2
+    Parameters:
+        
+        A : Operator
+        
+        c : Scaling Constant
+        
+        b : Data
+        
+    L : Lipshitz Constant of the gradient of :math:`F` which is :math:`2c||A||_2^2 = 2s1(A)^2`,
     
     where s1(A) is the largest singular value of A.
+        
     
-    '''
+    """
     
     def __init__(self, A, b, c=1.0):
-        super(Norm2Sq, self).__init__()
+        super(LeastSquares, self).__init__()
     
         self.A = A  # Should be an operator, default identity
         self.b = b  # Default zero DataSet?
@@ -66,21 +70,16 @@ class Norm2Sq(Function):
             warnings.warn('{} could not calculate Lipschitz Constant. {}'.format(
                 self.__class__.__name__, noe))
         
-    #def grad(self,x):
-    #    return self.gradient(x, out=None)
-
     def __call__(self, x):
-        #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel()))
-        #if out is None:
-        #    return self.c*( ( (self.A.direct(x)-self.b)**2).sum() )
-        #else:
+        
+        r""" Returns the value of :math:`F(x) = c\|Ax-b\|_2^2`
+        """
+
         y = self.A.direct(x)
         y.subtract(self.b, out=y)
-        #y.__imul__(y)
-        #return y.sum() * self.c
         try:
-            if self.c == 1:
-                return y.squared_norm()
+#            if self.c == 1:
+#                return y.squared_norm()
             return y.squared_norm() * self.c
         except AttributeError as ae:
             # added for compatibility with SIRF
@@ -91,6 +90,13 @@ class Norm2Sq(Function):
             return (yn * yn) * self.c
     
     def gradient(self, x, out=None):
+        
+        r""" Returns the value of the gradient of :math:`F(x) = c*\|A*x-b\|_2^2`
+        
+             .. math:: F'(x) = 2cA^T(Ax-b)
+
+        """
+        
         if out is not None:
             #return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b )
             self.A.direct(x, out=self.range_tmp)
@@ -100,3 +106,8 @@ class Norm2Sq(Function):
             out.multiply (self.c * 2.0, out=out)
         else:
             return (2.0*self.c)*self.A.adjoint(self.A.direct(x) - self.b)
+        
+
+    
+    
+    
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index 1af0e77..77b483e 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -19,109 +19,170 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
-from ccpi.optimisation.functions import Function, ScaledFunction
+
+from ccpi.optimisation.functions import Function
 from ccpi.framework import BlockDataContainer
 import numpy as np
 
-import functools
 
 class MixedL21Norm(Function):
     
     
-    '''
-        .. math::
-          
-          f(x) = ||x||_{2,1} = \sum |x|_{2}
-    '''      
+    """ MixedL21Norm function: :math:`F(x) = ||x||_{2,1} = \sum |x|_{2} = \sum \sqrt{ (x^{1})^{2} + (x^{2})^{2} + \dots}`                  
+    
+        where x is a BlockDataContainer, i.e., :math:`x=(x^{1}, x^{2}, \dots)`
+    
+    """      
     
     def __init__(self, **kwargs):
-        '''creator        
+        '''Creator
+        
+        :param b:  translation of the function
+        :type b: :code:`DataContainer`, optional
         '''
-        super(MixedL21Norm, self).__init__()                      
+        super(MixedL21Norm, self).__init__()  
+        self.b = kwargs.get('b', None)  
+
+        # This is to handle tensor for Total Generalised Variation                    
         self.SymTensor = kwargs.get('SymTensor',False)
         
-    def __call__(self, x):
+        # We use this parameter to make MixedL21Norm differentiable
+#        self.epsilon = epsilon
         
-        ''' Evaluates L2,1Norm at point x
+        if self.b is not None and not isinstance(self.b, BlockDataContainer):
+            raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(self.b)))
             
-            :param x: is a BlockDataContainer
-                                
-        '''
+        
+    def __call__(self, x):
+        
+        r"""Returns the value of the MixedL21Norm function at x. 
+
+        :param x: :code:`BlockDataContainer`                                           
+        """
         if not isinstance(x, BlockDataContainer):
             raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) 
-                                         
-        tmp = x.get_item(0) * 0.
-        for el in x.containers:
-            tmp += el.power(2.)
-        return tmp.sqrt().sum()
+            
+#        tmp_cont = x.containers                        
+#        tmp = x.get_item(0) * 0.
+#        for el in tmp_cont:
+#            tmp += el.power(2.)
+#        tmp.add(self.epsilon**2, out = tmp)    
+#        return tmp.sqrt().sum()            
+         
+        y = x
+        if self.b is not None:
+            y = x - self.b    
+        return y.pnorm(p=2).sum()                                
+            
+#        tmp = x.get_item(0) * 0.
+#        for el in x.containers:
+#            tmp += el.power(2.)
+#        #tmp.add(self.epsilon, out = tmp)    
+#        return tmp.sqrt().sum()
 
                             
-    def gradient(self, x, out=None):
-        return ValueError('Not Differentiable')
-                            
     def convex_conjugate(self,x):
         
-        ''' This is the Indicator function of :math:`||\cdot||_{2, \infty}` which is either 0 if :math:`||x||_{2, \infty}` or :math:`\infty`
-
-        Notice this returns 0
-        '''
+        r"""Returns the value of the convex conjugate of the MixedL21Norm function at x.
         
-        return 0.0
+        This is the Indicator function of :math:`\mathbb{I}_{\{\|\cdot\|_{2,\infty}\leq1\}}(x^{*})`,
         
-    
+        i.e., 
+        
+        .. math:: \mathbb{I}_{\{\|\cdot\|_{2, \infty}\leq1\}}(x^{*}) 
+            = \begin{cases} 
+            0, \mbox{if } \|x\|_{2, \infty}\leq1\\
+            \infty, \mbox{otherwise}
+            \end{cases}
+        
+        where, 
+        
+        .. math:: \|x\|_{2,\infty} = \max\{ \|x\|_{2} \} = \max\{ \sqrt{ (x^{1})^{2} + (x^{2})^{2} + \dots}\}
+        
+        """
+        if not isinstance(x, BlockDataContainer):
+            raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) 
+                
+#        tmp1 = x.get_item(0) * 0.
+#        for el in x.containers:
+#            tmp1 += el.power(2.)
+#        tmp1.add(self.epsilon**2, out = tmp1)
+#        tmp = tmp1.sqrt().as_array().max() - 1
+#                    
+#        if tmp<=1e-6:
+#            return 0
+#        else:
+#            return np.inf            
+                
+        tmp = (x.pnorm(2).as_array().max() - 1)
+        if tmp<=1e-5:
+            return 0
+        else:
+            return np.inf
+                    
     def proximal(self, x, tau, out=None):
         
+        r"""Returns the value of the proximal operator of the MixedL21Norm function at x.
+        
+        .. math :: \mathrm{prox}_{\tau F}(x) = \frac{x}{\|x\|_{2}}\max\{ \|x\|_{2} - \tau, 0 \}
+        
+        where the convention 0 · (0/0) = 0 is used.
+        
+        """
+        
         if out is None:
             
-            tmp = sum([ el*el for el in x.containers]).sqrt()
+            tmp = x.pnorm(2)
             res = (tmp - tau).maximum(0.0) * x/tmp
-            return res
+            
+            for el in res.containers:
+                el.as_array()[np.isnan(el.as_array())]=0            
+            
+            return res            
             
         else:
-                        
-            tmp = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ).sqrt()
+            
+            tmp = x.pnorm(2)
             res = (tmp - tau).maximum(0.0) * x/tmp
 
             for el in res.containers:
                 el.as_array()[np.isnan(el.as_array())]=0
 
-            out.fill(res)
-        
-    
-    def proximal_conjugate(self, x, tau, out=None): 
-
-        
-        if out is None:                                        
-            tmp = x.get_item(0) * 0	
-            for el in x.containers:	
-                tmp += el.power(2.)	
-            tmp.sqrt(out=tmp)	
-            tmp.maximum(1.0, out=tmp)	
-            frac = [ el.divide(tmp) for el in x.containers ]	
-            return BlockDataContainer(*frac)
+            out.fill(res)            
+            
+                        
+#            tmp = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ).sqrt()
+#            res = (tmp - tau).maximum(0.0) * x/tmp
+#
+#            for el in res.containers:
+#                el.as_array()[np.isnan(el.as_array())]=0
+#
+#            out.fill(res)
         
-    
-        else:
-                            
-            res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
-            res1.sqrt(out=res1)	
-            res1.maximum(1.0, out=res1)	
-            x.divide(res1, out=out)
+##############################################################################
+        ##############################################################################
+#    def proximal_conjugate(self, x, tau, out=None): 
+#
+#        
+#        if out is None:                                        
+#            tmp = x.get_item(0) * 0	
+#            for el in x.containers:	
+#                tmp += el.power(2.)	
+#            tmp.sqrt(out=tmp)	
+#            tmp.maximum(1.0, out=tmp)	
+#            frac = [ el.divide(tmp) for el in x.containers ]	
+#            return BlockDataContainer(*frac)
+#        
+#    
+#        else:
+#                            
+#            res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
+#            res1.sqrt(out=res1)	
+#            res1.maximum(1.0, out=res1)	
+#            x.divide(res1, out=out)
                               
 
-    def __rmul__(self, scalar):
-        
-        ''' Multiplication of MixedL21Norm with a scalar
-        
-        Returns: ScaledFunction
-             
-        '''         
-        return ScaledFunction(self, scalar) 
-
-
-#
 if __name__ == '__main__':
     
     M, N, K = 2,3,50
@@ -189,8 +250,24 @@ if __name__ == '__main__':
                                             d1.get_item(0).as_array(), decimal=4) 
 
     numpy.testing.assert_array_almost_equal(out1.get_item(1).as_array(), \
-                                            d1.get_item(1).as_array(), decimal=4)      
-#    
+                                            d1.get_item(1).as_array(), decimal=4)   
+    
+    f_scaled.proximal_conjugate(U, tau, out = out1)
+    x = U
+    tmp = x.get_item(0) * 0	
+    for el in x.containers:	
+        tmp += el.power(2.)	
+    tmp.sqrt(out=tmp)	
+    (tmp/f_scaled.scalar).maximum(1.0, out=tmp)	
+    frac = [ el.divide(tmp) for el in x.containers ]	
+    out2 = BlockDataContainer(*frac)   
+    
+    numpy.testing.assert_array_almost_equal(out1.get_item(0).as_array(), \
+                                            out2.get_item(0).as_array(), decimal=4)       
+    
+    
+
+      
     
     
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
deleted file mode 100755
index 2bf3b25..0000000
--- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
+++ /dev/null
@@ -1,157 +0,0 @@
-# -*- coding: utf-8 -*-

-# Copyright 2019 Science Technology Facilities Council

-# Copyright 2019 University of Manchester

-#

-# This work is part of the Core Imaging Library developed by Science Technology

-# Facilities Council and 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.txt

-#

-# 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 __future__ import absolute_import

-from __future__ import division

-from __future__ import print_function

-from __future__ import unicode_literals

-

-from numbers import Number

-import numpy

-import warnings

-

-class ScaledFunction(object):

-    

-    '''ScaledFunction

-

-    A class to represent the scalar multiplication of an Function with a scalar.

-    It holds a function and a scalar. Basically it returns the multiplication

-    of the product of the function __call__, convex_conjugate and gradient with the scalar.

-    For the rest it behaves like the function it holds.

-

-    Args:

-       function (Function): a Function or BlockOperator

-       scalar (Number): a scalar multiplier

-    Example:

-       The scaled operator behaves like the following:

-       

-    '''

-    def __init__(self, function, scalar):

-        super(ScaledFunction, self).__init__()

-

-        if not isinstance (scalar, Number):

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

-        self.scalar = scalar

-        self.function = function

-        

-        if self.function.L is not None:        

-            self.L = abs(self.scalar) * self.function.L         

-

-    def __call__(self,x, out=None):

-        '''Evaluates the function at x '''

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

-

-    def convex_conjugate(self, x):

-        '''returns the convex_conjugate of the scaled function '''

-        return self.scalar * self.function.convex_conjugate(x/self.scalar)

-    

-    def gradient(self, x, out=None):

-        '''Returns the gradient of the function at x, if the function is differentiable'''

-        if out is None:            

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

-        else:

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

-            out *= self.scalar

-

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

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

-        '''

-        if out is None:

-            return self.function.proximal(x, tau*self.scalar)     

-        else:

-            self.function.proximal(x, tau*self.scalar, out = out)

-

-    def proximal_conjugate(self, x, tau, out = None):

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

-        '''

-        if out is None:

-            return self.scalar * self.function.proximal_conjugate(x/self.scalar, tau/self.scalar)

-        else:

-            self.function.proximal_conjugate(x/self.scalar, tau/self.scalar, out=out)

-            out *= self.scalar

-

-    def grad(self, x):

-        '''Alias of gradient(x,None)'''

-        warnings.warn('''This method will disappear in following 

-        versions of the CIL. Use gradient instead''', DeprecationWarning)

-        return self.gradient(x, out=None)

-

-    def prox(self, x, tau):

-        '''Alias of proximal(x, tau, None)'''

-        warnings.warn('''This method will disappear in following 

-        versions of the CIL. Use proximal instead''', DeprecationWarning)

-        return self.proximal(x, tau, out=None)

-

-

-            

-if __name__ == '__main__':        

-

-    from ccpi.optimisation.functions import L2NormSquared, MixedL21Norm

-    from ccpi.framework import ImageGeometry, BlockGeometry

-    

-    M, N, K = 2,3,5

-    ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K)

-    

-    u = ig.allocate('random_int')

-    b = ig.allocate('random_int') 

-    

-    BG = BlockGeometry(ig, ig)

-    U = BG.allocate('random_int')

-        

-    f2 = 0.5 * L2NormSquared(b=b)

-    f1 = 30 * MixedL21Norm()

-    tau = 0.355

-    

-    res_no_out1 =  f1.proximal_conjugate(U, tau)

-    res_no_out2 =  f2.proximal_conjugate(u, tau)    

-    

-    

-#    print( " ######## with out ######## ")

-    res_out1 = BG.allocate()

-    res_out2 = ig.allocate()

-    

-    f1.proximal_conjugate(U, tau, out = res_out1)

-    f2.proximal_conjugate(u, tau, out = res_out2)

-

-

-    numpy.testing.assert_array_almost_equal(res_no_out1[0].as_array(), \

-                                            res_out1[0].as_array(), decimal=4) 

-    

-    numpy.testing.assert_array_almost_equal(res_no_out2.as_array(), \

-                                            res_out2.as_array(), decimal=4)     

-      

-          

-          

-    

-    

-    

-    

-    

-    

-    

-    

-    

-    

-    

-    

-

-    

-

-

-    

diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py
deleted file mode 100644
index 19db668..0000000
--- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py
+++ /dev/null
@@ -1,87 +0,0 @@
-# -*- coding: utf-8 -*-
-# Copyright 2019 Science Technology Facilities Council
-# Copyright 2019 University of Manchester
-#
-# This work is part of the Core Imaging Library developed by Science Technology
-# Facilities Council and 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.txt
-#
-# 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 __future__ import absolute_import
-from __future__ import division
-from __future__ import print_function
-from __future__ import unicode_literals
-
-from ccpi.optimisation.functions import Function
-
-class ZeroFunction(Function):
-    
-    r'''ZeroFunction: .. math:: f(x) = 0,         
-        
-        Maps evely element x\in X to zero
-    '''
-    
-    def __init__(self):
-        super(ZeroFunction, self).__init__()
-              
-    def __call__(self,x):
-        
-        '''Evaluates ZeroFunction at x'''
-        return 0
-    
-    
-    def gradient(self, x, out=None):
-        
-        '''Evaluates gradient of ZeroFunction at x'''        
-        
-        if out is None:
-            return 0
-        else:
-            out *= 0
-    
-    def convex_conjugate(self, x):
-        
-        r''' Convex conjugate of ZeroFunction: support function .. math:: sup <x, x^{*}>  
-            
-            In fact is the  indicator function for the set = {0}
-            So 0 if x=0, or inf if x neq 0  
-            
-        '''
-        return x.maximum(0).sum()
-        
-    
-    def proximal(self, x, tau, out=None):
-        
-        r'''Proximal operator of ZeroFunction at x
-           
-            .. math:: prox_{\tau * f}(x)
-        '''
-        
-        if out is None:
-            return x.copy()
-        else:
-            out.fill(x)
-        
-    def proximal_conjugate(self, x, tau, out = None):
-        
-        r'''Proximal operator of the convex conjugate of ZeroFunction at x:
-           
-           .. math:: prox_{\tau * f^{*}}(x)
-
-        '''         
-        
-        if out is None:
-            return 0
-        else:
-            return 0
diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
index d968539..f75f82a 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
@@ -15,15 +15,17 @@
 #   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 .Function import Function
-from .ZeroFunction import ZeroFunction
+
+from __future__ import absolute_import
+
+from .Function import Function, ConstantFunction, ZeroFunction, TranslateFunction, SumFunctionScalar
+from .Function import ScaledFunction
 from .L1Norm import L1Norm
 from .L2NormSquared import L2NormSquared
-from .ScaledFunction import ScaledFunction
+from .LeastSquares import LeastSquares
 from .BlockFunction import BlockFunction
 from .FunctionOperatorComposition import FunctionOperatorComposition
 from .MixedL21Norm import MixedL21Norm
 from .IndicatorBox import IndicatorBox
 from .KullbackLeibler import KullbackLeibler
-from .Norm2Sq import Norm2Sq
 from .Rosenbrock import Rosenbrock
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 23cb799..a2a9b9e 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -20,7 +20,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 import functools
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
index eeecee9..9218ad4 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
@@ -21,7 +21,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from numbers import Number
 import numpy
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 3c563fb..1df7745 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.operators import LinearOperator
 import numpy as np
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 2ff0b20..6391cf7 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -18,7 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 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
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index e95234b..a4c6ae3 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -19,7 +19,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 from ccpi.optimisation.operators import LinearOperator
 import scipy.sparse as sp
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index fb09819..5ff5b01 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -19,7 +19,6 @@
 from __future__ import absolute_import

 from __future__ import division

 from __future__ import print_function

-from __future__ import unicode_literals

 

 from ccpi.optimisation.operators import Operator

 import numpy

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index a84ea94..6f13532 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -19,7 +19,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
 
 import numpy
 from scipy.sparse.linalg import svds
@@ -27,7 +26,11 @@ from ccpi.framework import VectorGeometry
 from ccpi.optimisation.operators import LinearOperator
 
 class LinearOperatorMatrix(LinearOperator):
-    '''Matrix wrapped into a LinearOperator'''
+    """ Matrix wrapped into a LinearOperator
+    
+    :param: a numpy matrix 
+    
+    """
     
     def __init__(self,A):
         '''creator
@@ -41,17 +44,25 @@ class LinearOperatorMatrix(LinearOperator):
         self.s1 = None   # Largest singular value, initially unknown
         super(LinearOperatorMatrix, self).__init__()
         
-    def direct(self,x, out=None):
+    def direct(self,x, out=None):        
+                
         if out is None:
-            return type(x)(numpy.dot(self.A,x.as_array()))
+            tmp = self.gm_range.allocate()
+            tmp.fill(numpy.dot(self.A,x.as_array()))
+            return tmp
         else:
-            numpy.dot(self.A, x.as_array(), out=out.as_array())
+            # Below use of out is not working, see
+            # https://docs.scipy.org/doc/numpy/reference/generated/numpy.dot.html            
+            # numpy.dot(self.A, x.as_array(), out = out.as_array())            
+            out.fill(numpy.dot(self.A, x.as_array()))
 
     def adjoint(self,x, out=None):
         if out is None:
-            return type(x)(numpy.dot(self.A.transpose(),x.as_array()))
-        else:
-            numpy.dot(self.A.transpose(),x.as_array(), out=out.as_array())
+            tmp = self.gm_domain.allocate()
+            tmp.fill(numpy.dot(self.A.transpose(),x.as_array()))
+            return tmp
+        else:            
+            out.fill(numpy.dot(self.A.transpose(),x.as_array()))
 
     def size(self):
         return self.A.shape
diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 87059e6..3cd0899 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -18,8 +18,6 @@
 from __future__ import absolute_import

 from __future__ import division

 from __future__ import print_function

-from __future__ import unicode_literals

-

 

 from ccpi.optimisation.operators.ScaledOperator import ScaledOperator

 

diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
index d1ad07c..7f1e202 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
@@ -15,6 +15,9 @@
 #   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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 from numbers import Number
 import numpy
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
index 9239d90..b0a3b51 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
@@ -18,8 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
-
 
 from ccpi.framework import DataContainer
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 698a993..747db65 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -18,8 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
-
 
 import scipy.sparse as sp
 import numpy as np
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index d82c5c0..6fd2e86 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -18,8 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
-
 
 from ccpi.optimisation.operators import LinearOperator
 from ccpi.framework import BlockGeometry, BlockDataContainer
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index f677dc2..0feafd8 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -18,8 +18,6 @@
 from __future__ import absolute_import
 from __future__ import division
 from __future__ import print_function
-from __future__ import unicode_literals
-
 
 import numpy as np
 from ccpi.framework import ImageData
diff --git a/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py
index f9ed19d..9f67c12 100755..100644
--- a/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py
+++ b/Wrappers/Python/ccpi/processors/CenterOfRotationFinder.py
@@ -15,6 +15,9 @@
 #   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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 
 from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\
  AcquisitionGeometry, ImageGeometry, ImageData
diff --git a/Wrappers/Python/ccpi/processors/Normalizer.py b/Wrappers/Python/ccpi/processors/Normalizer.py
index 905a06b..2872348 100755..100644
--- a/Wrappers/Python/ccpi/processors/Normalizer.py
+++ b/Wrappers/Python/ccpi/processors/Normalizer.py
@@ -15,6 +15,9 @@
 #   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 __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
 
 from ccpi.framework import DataProcessor, DataContainer, AcquisitionData,\
  AcquisitionGeometry, ImageGeometry, ImageData
diff --git a/Wrappers/Python/ccpi/processors/Resizer.py b/Wrappers/Python/ccpi/processors/Resizer.py
index 72ea392..afa2a42 100755..100644
--- a/Wrappers/Python/ccpi/processors/Resizer.py
+++ b/Wrappers/Python/ccpi/processors/Resizer.py
@@ -15,6 +15,9 @@
 #   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 __future__ import absolute_import

+from __future__ import division

+from __future__ import print_function

 from ccpi.framework import DataProcessor, AcquisitionData, ImageData

 import warnings

 

diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index ed2cd25..8438311 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -30,10 +30,10 @@ setup(
     name="ccpi-framework",
     version=cil_version,
     packages=['ccpi' , 'ccpi.io',
-              'ccpi.framework', 'ccpi.optimisation', 
-              'ccpi.optimisation.operators',
-              'ccpi.optimisation.algorithms',
+              'ccpi.framework', 'ccpi.optimisation',
               'ccpi.optimisation.functions',
+              'ccpi.optimisation.algorithms',
+              'ccpi.optimisation.operators',
               'ccpi.processors',
               'ccpi.utilities', 'ccpi.utilities.jupyter',
               'ccpi.contrib','ccpi.contrib.optimisation',
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index 96c48da..5c659a4 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -16,14 +16,13 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 import unittest
-from ccpi.framework import ImageGeometry, ImageData, BlockDataContainer, DataContainer
+from ccpi.framework import ImageGeometry, VectorGeometry, ImageData, BlockDataContainer, DataContainer
 from ccpi.optimisation.operators import BlockOperator, BlockScaledOperator,\
     FiniteDiff, SymmetrizedGradient
 import numpy
 from timeit import default_timer as timer
-from ccpi.framework import ImageGeometry
 from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
-from ccpi.optimisation.operators import LinearOperator
+from ccpi.optimisation.operators import LinearOperator, LinearOperatorMatrix
 
 
 def dt(steps):
@@ -66,6 +65,39 @@ class CCPiTestClass(unittest.TestCase):
 
 
 class TestOperator(CCPiTestClass):
+    
+    def test_LinearOperatorMatrix(self):
+        
+        print('Check LinearOperatorMatrix')
+                
+        m = 30
+        n = 20
+        
+        vg = VectorGeometry(n)
+        
+        Amat = numpy.random.randn(m, n)
+        A = LinearOperatorMatrix(Amat)
+        
+        b = vg.allocate('random')
+        
+        out1 = A.range_geometry().allocate()
+        out2 = A.domain_geometry().allocate()
+        
+        res1 = A.direct(b)
+        res2 = numpy.dot(A.A, b.as_array())
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), res2)
+            
+        A.direct(b, out = out1)
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), out1.as_array(), decimal=4)
+        
+        res3 = A.adjoint(res1)
+        res4 = numpy.dot(A.A.transpose(),res1.as_array())
+        self.assertNumpyArrayAlmostEqual(res3.as_array(), res4, decimal=4)   
+        
+        A.adjoint(res1, out = out2)
+        self.assertNumpyArrayAlmostEqual(res3.as_array(), out2.as_array(), decimal=4)        
+    
+    
     def test_ScaledOperator(self):
         print ("test_ScaledOperator")
         ig = ImageGeometry(10,20,30)
diff --git a/Wrappers/Python/test/test_SumFunction.py b/Wrappers/Python/test/test_SumFunction.py
new file mode 100644
index 0000000..465c084
--- /dev/null
+++ b/Wrappers/Python/test/test_SumFunction.py
@@ -0,0 +1,341 @@
+# -*- 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 ccpi.optimisation.functions import L1Norm, ScaledFunction, \
+                                        LeastSquares, L2NormSquared, \
+                                        KullbackLeibler, ZeroFunction, ConstantFunction
+from ccpi.optimisation.operators import Identity                                        
+from ccpi.framework import ImageGeometry
+
+import unittest
+import numpy
+from numbers import Number
+
+''' Here we test SumFunction class for different function
+
+L2Norm, L1Norm, KullbackLeibler, ZeroFunction, ConstantFunction, Scalar
+
+for call method
+for gradient method
+
+
+
+'''
+
+
+
+class TestFunction(unittest.TestCase):
+    
+    def assertBlockDataContainerEqual(self, container1, container2):
+        print ("assert Block Data Container Equal")
+        self.assertTrue(issubclass(container1.__class__, container2.__class__))
+        for col in range(container1.shape[0]):
+            if issubclass(container1.get_item(col).__class__, DataContainer):
+                print ("Checking col ", col)
+                self.assertNumpyArrayEqual(
+                    container1.get_item(col).as_array(), 
+                    container2.get_item(col).as_array()
+                    )
+            else:
+                self.assertBlockDataContainerEqual(container1.get_item(col),container2.get_item(col))
+    
+    def assertNumpyArrayEqual(self, first, second):
+        res = True
+        try:
+            numpy.testing.assert_array_equal(first, second)
+        except AssertionError as err:
+            res = False
+            print(err)
+        self.assertTrue(res)
+
+    def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
+        res = True
+        try:
+            numpy.testing.assert_array_almost_equal(first, second, decimal)
+        except AssertionError as err:
+            res = False
+            print(err)
+            print("expected " , second)
+            print("actual " , first)
+
+        self.assertTrue(res)
+        
+    def test_SumFunction(self):
+        
+        M, N, K = 3,4,5
+        ig = ImageGeometry(M, N, K)
+        
+        tmp = ig.allocate('random_int')
+        b = ig.allocate('random_int')
+        
+        operator = Identity(ig)
+
+        scalar = 0.25
+        f1 = L2NormSquared()
+        f2 = L1Norm()
+        f3 = scalar * L2NormSquared()
+        f4 = scalar * L1Norm()
+        f5 = scalar * L2NormSquared(b=b)
+        f6 = scalar * L1Norm(b=b)  
+        f7 = ZeroFunction()
+        f8 = 5 *  ConstantFunction(10)             
+        f9 = LeastSquares(operator, b, c=scalar)
+        f10 = 0.5*KullbackLeibler(b=b)
+        f11 = KullbackLeibler(b=b)
+        f12 = 10
+        
+#        f10 = 0.5 * MixedL21Norm()
+#        f11 = IndicatorBox(lower=0)
+        
+        list1 = [f1, f2, f3, f4, f5, f6, f7, f8, f9, f10, f11, f12]
+                
+        print('###################  Check sum of two functions ################## \n')
+        
+        for func in list1:
+               
+            
+            # check sum of two functions   
+            
+            if isinstance(func, ScaledFunction):
+                type_fun = ' scalar * ' + type(func.function).__name__
+            else:    
+                type_fun = type(func).__name__
+                
+            if isinstance(func, Number):
+                tmp_fun_eval = func
+            else:
+                tmp_fun_eval = func(tmp)                
+                             
+            sumf = f1 + func           
+            self.assertNumpyArrayAlmostEqual(sumf(tmp), f1(tmp) + tmp_fun_eval )
+            print('{} = ( {} + {} ) is OK'.format(type(sumf).__name__, type(f1).__name__, type_fun))
+            
+            sumf1 = func + f1 
+            self.assertNumpyArrayAlmostEqual(sumf1(tmp), tmp_fun_eval + f1(tmp))
+            print('Checking commutative')
+            print('{} + ( {} + {} ) is OK\n'.format(type(sumf1).__name__, type_fun, type(f1).__name__))
+            
+        print('###################  Check Lispchitz constant ################## \n')
+        
+        for func in list1:
+            
+            if isinstance(func, ScaledFunction):
+                type_fun = ' scalar * ' + type(func.function).__name__
+            else:    
+                type_fun = type(func).__name__            
+               
+            
+            # check Lispchitz sum of two functions  
+            
+            if isinstance(func, Number):
+                tmp_fun_L = 0
+            else:
+                tmp_fun_L = func.L           
+            
+            sumf = f1 + func   
+            
+            try:
+                sumf.L==f1.L + tmp_fun_L
+            except TypeError:
+                print('Function {} has L = None'.format(type_fun))
+                
+        print('\n###################  Check Gradient ################## \n')   
+              
+              
+        for func in list1:
+              
+            if isinstance(func, ScaledFunction):
+                type_fun = ' scalar * ' + type(func.function).__name__
+            else:    
+                type_fun = type(func).__name__
+                                          
+            sumf = f1 + func
+            # check gradient          
+            try:
+                if isinstance(func, Number):
+                    tmp_fun_gradient = 0
+                else:
+                    tmp_fun_gradient = func.gradient(tmp)   
+                    
+                self.assertNumpyArrayAlmostEqual(sumf.gradient(tmp).as_array(), (f1.gradient(tmp) + tmp_fun_gradient).as_array())
+            except NotImplementedError:
+                print("{} is not differentiable".format(type_fun))
+                
+        print('\n###################  Check Gradient Out ################## \n')  
+         
+        out_left = ig.allocate()
+        out_right1 = ig.allocate()
+        out_right2 = ig.allocate()  
+            
+        for i, func in enumerate(list1):               
+                
+            if isinstance(func, ScaledFunction):
+                type_fun = ' scalar * ' + type(func.function).__name__
+            else:    
+                type_fun = type(func).__name__
+            
+            sumf = f1 + func
+            
+                                    
+            # check gradient out    
+            try:
+                
+                
+                if isinstance(func, Number):
+                    tmp_fun_gradient_out = 0
+                else:
+                    func.gradient(tmp, out = out_right2) 
+                    tmp_fun_gradient_out = out_right2.as_array()
+                    
+                #print('Check {} + {}\n'.format(type(f1).__name__, type_fun))
+                sumf.gradient(tmp, out = out_left)
+                f1.gradient(tmp, out = out_right1)   
+                self.assertNumpyArrayAlmostEqual(out_left.as_array(), out_right1.as_array() + tmp_fun_gradient_out)
+            except NotImplementedError:
+                print("{} is not differentiable".format(type_fun))  
+    def test_ConstantFunction(self):
+
+        k = ConstantFunction(constant=1)
+        ig = ImageGeometry(1,2,3)
+        x = ig.allocate(2)
+
+        grad = k.gradient(x)
+        out = ig.allocate(-1)
+
+        k.gradient(x, out=out)
+        #out.fill(-3)
+        
+        self.assertNumpyArrayEqual(numpy.zeros(x.shape), grad.as_array())
+        
+        self.assertNumpyArrayEqual(out.as_array(), grad.as_array())
+                
+    def test_SumFunctionScalar(self):      
+        
+        M, N, K = 3,4,5
+        ig = ImageGeometry(M, N, K)
+        
+        tmp = ig.allocate('random_int')
+        b = ig.allocate('random_int')
+        
+        scalar = 0.25
+        f1 =  scalar * L2NormSquared(b=b)
+        f2 = 5
+           
+        f = f1 + f2
+        print(f)
+        
+        g = f2 + f1
+        print(g)
+
+        tau = 0.03
+        
+        # check call method       
+        res1 = f(tmp)
+        res2 = f1(tmp) + f2
+        self.assertAlmostEqual(res1, res2)
+        
+        # check gradient
+        res1 = f.gradient(tmp)
+        res2 = f1.gradient(tmp)
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), res2.as_array())
+        
+        # check gradient with out
+        out1 = tmp*0
+        out2 = tmp*0
+        f.gradient(tmp, out=out1)
+        f1.gradient(tmp, out=out2)
+        self.assertNumpyArrayAlmostEqual(out1.as_array(), out2.as_array())
+        
+
+        res1 = f.proximal(tmp, tau)
+        res2 = f1.proximal(tmp, tau)
+        
+        # proximal of sum of function with scalar = proximal of function
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), res2.as_array())
+        
+        # proximal  with out of sum of function with scalar = proximal of function
+        res1_out = ig.allocate()
+        res2_out = ig.allocate()         
+        f.proximal(tmp, tau, out = res1_out)
+        f1.proximal(tmp, tau, out = res2_out)             
+        self.assertNumpyArrayAlmostEqual(res1_out.as_array(), res2_out.as_array())
+
+        res1 = f.proximal_conjugate(tmp, tau)
+        res2 = f1.proximal_conjugate(tmp, tau)
+        
+        # proximal of sum of function with scalar = proximal of function
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), res2.as_array())
+        
+        # proximal  with out of sum of function with scalar = proximal of function
+        res1_out = ig.allocate()
+        res2_out = ig.allocate()         
+        f.proximal_conjugate(tmp, tau, out = res1_out)
+        f1.proximal_conjugate(tmp, tau, out = res2_out)             
+        self.assertNumpyArrayAlmostEqual(res1_out.as_array(), res2_out.as_array())        
+            
+        
+def test_ConstantFunction(self):      
+        
+        M, N, K = 3,4,5
+        ig = ImageGeometry(M, N, K)
+        
+        tmp = ig.allocate('random_int')
+        
+        constant = 10
+        f = ConstantFunction(constant)
+        
+        # check call
+        res1 = f(tmp)
+        self.assertAlmostEqual(res1, constant)
+        
+        # check gradient with and without out
+        res1 = f.gradient(tmp)
+        out = ig.allocate()
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), out)
+        
+        out1 = ig.allocate()
+        f.gradient(tmp, out=out1)
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), out1)
+        
+        # check convex conjugate        
+        res1 = f.convex_conjugate(tmp)
+        res2 = tmp.maximum(0).sum()
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), res2.as_array())
+        
+        # check proximal 
+        tau = 0.4
+        res1 = f.proximal(tmp, tau)
+        self.assertNumpyArrayAlmostEqual(res1.as_array(), tmp.as_array()) 
+            
+        
+        
+        # check call
+        
+        
+        
+     
+        
+                
+#if __name__ == '__main__':
+#    
+#    t = TestFunction()
+#    t.test_SumFunction()
+#    t.test_SumFunctionScalar()
+
+                
+
diff --git a/Wrappers/Python/test/test_TranslateFunction.py b/Wrappers/Python/test/test_TranslateFunction.py
new file mode 100644
index 0000000..f739d68
--- /dev/null
+++ b/Wrappers/Python/test/test_TranslateFunction.py
@@ -0,0 +1,272 @@
+# -*- 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 ccpi.optimisation.functions import Function, L1Norm, ScaledFunction, \
+                                        LeastSquares, L2NormSquared, \
+                                        KullbackLeibler, ZeroFunction, ConstantFunction, TranslateFunction
+from ccpi.optimisation.operators import Identity                                        
+from ccpi.framework import ImageGeometry, BlockGeometry
+
+import unittest
+import numpy
+from numbers import Number
+
+
+''' Here we test SumFunction class for different function
+
+L2Norm, L1Norm, KullbackLeibler, ZeroFunction, ConstantFunction, Scalar
+
+for call method
+for gradient method
+
+
+
+'''
+
+
+
+class TestFunction(unittest.TestCase):
+    
+    def assertBlockDataContainerEqual(self, container1, container2):
+        print ("assert Block Data Container Equal")
+        self.assertTrue(issubclass(container1.__class__, container2.__class__))
+        for col in range(container1.shape[0]):
+            if issubclass(container1.get_item(col).__class__, DataContainer):
+                print ("Checking col ", col)
+                self.assertNumpyArrayEqual(
+                    container1.get_item(col).as_array(), 
+                    container2.get_item(col).as_array()
+                    )
+            else:
+                self.assertBlockDataContainerEqual(container1.get_item(col),container2.get_item(col))
+    
+    def assertNumpyArrayEqual(self, first, second):
+        res = True
+        try:
+            numpy.testing.assert_array_equal(first, second)
+        except AssertionError as err:
+            res = False
+            print(err)
+        self.assertTrue(res)
+
+    def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
+        res = True
+        try:
+            numpy.testing.assert_array_almost_equal(first, second, decimal)
+        except AssertionError as err:
+            res = False
+            print(err)
+            print("expected " , second)
+            print("actual " , first)
+
+        self.assertTrue(res)
+        
+    def test_TranslateFunction(self):   
+        
+        # Test TranslationFunction
+    
+        ig = ImageGeometry(4,4)
+        tmp = ig.allocate('random_int')
+        b = ig.allocate('random_int')
+        scalar = 0.4
+        tau = 0.05
+                                
+        list1 = [ L2NormSquared(), scalar * L2NormSquared(), scalar * L2NormSquared(b=b), 
+                 L1Norm(), scalar * L1Norm(), scalar * L1Norm(b=b)]
+        
+        list1_shift = [ L2NormSquared().centered_at(ig.allocate()), scalar * L2NormSquared().centered_at(ig.allocate()), scalar * L2NormSquared().centered_at(b),
+                        L1Norm().centered_at(ig.allocate()), scalar * L1Norm().centered_at(ig.allocate()), scalar * L1Norm().centered_at(b)]        
+        
+        out_gradient1 = ig.allocate()
+        out_gradient2 = ig.allocate()
+        
+        out_proximal1 = ig.allocate()
+        out_proximal2 = ig.allocate()     
+        
+        out_proximal_conj1 = ig.allocate()
+        out_proximal_conj2 = ig.allocate()             
+        
+        for func, func_shift in zip(list1, list1_shift):
+
+            # check call
+            res1 = func(tmp)
+            res2 = func_shift(tmp)
+            self.assertNumpyArrayAlmostEqual(res1, res2)
+            
+            try:
+                # check gradient
+                res1_gradient = func.gradient(tmp)
+                res2_gradient = func_shift.gradient(tmp)
+                self.assertNumpyArrayAlmostEqual(res1_gradient.as_array(), res2_gradient.as_array()) 
+                
+                # check gradient out
+                func.gradient(tmp, out = out_gradient1)
+                func_shift.gradient(tmp, out = out_gradient2)
+                self.assertNumpyArrayAlmostEqual(out_gradient1.as_array(), out_gradient2.as_array())
+                
+            except NotImplementedError:
+                print('Function is not differentiable')
+                             
+            # check proximal     
+            func.proximal(tmp, tau, out = out_proximal1)
+            func_shift.proximal(tmp, tau, out = out_proximal2)
+            self.assertNumpyArrayAlmostEqual(out_proximal1.as_array(), out_proximal2.as_array())  
+
+            # check proximal conjugate
+            func.proximal_conjugate(tmp, tau, out = out_proximal_conj1)
+            func_shift.proximal_conjugate(tmp, tau, out = out_proximal_conj2)
+            self.assertNumpyArrayAlmostEqual(out_proximal_conj1.as_array(), out_proximal_conj1.as_array())               
+                                
+            
+if __name__ == '__main__':
+#    
+    t = TestFunction()
+    t.test_TranslateFunction()    
+
+
+#    ig = ImageGeometry(4,4)
+#    tmp = ig.allocate('random_int')
+#    b = ig.allocate('random_int')
+#    scalar = 0.4
+#            
+##    f = scalar * L2NormSquared().centered_at(b)
+##    print(f.function.function)
+#    list1 = [ L2NormSquared(), scalar * L2NormSquared(), scalar * L2NormSquared(b=b)]                            
+#        
+##        for func in list_functions:
+##            
+###            if isinstance(func, ScaledFunction):
+###                func_tmp = func.function
+###            else:
+###                func_tmp = func                
+###                                
+###            if func_tmp.b is None:
+###                tmp_data = ig.allocate()
+###            else:
+###                tmp_data = b 
+##            
+##            func_tmp = func
+##            tmp_data = ig.allocate()
+##                                        
+##            res1 = func_tmp(tmp)        
+##            res2 = func_tmp.centered_at(tmp_data)(tmp)
+##            
+##            self.assertNumpyArrayAlmostEqual(res1, res2)
+        
+        
+
+        
+        
+#        
+#        for i in list_functions:
+#            
+#            print('Test Translation for Function {} '.format(type(i).__name__))
+#            
+#            if isinstance(i, L2NormSquared):
+#                
+#                f = L2NormSquared(b = b)    
+#                g = TranslateFunction(L2NormSquared(), b)
+#                
+#            elif isinstance(i, L1Norm):
+#                
+#                f = L1Norm(b = b)    
+#                g = TranslateFunction(L1Norm(), b)
+#                
+#            elif isinstance(i, ScaledFunction):
+#    
+#                if isinstance(i.function, L2NormSquared):
+#                    f = scalar * L2NormSquared(b = b)    
+#                    g = scalar * TranslateFunction(L2NormSquared(), b)
+#                    
+#                if isinstance(i.function, L1Norm):
+#                    f = scalar * L1Norm(b = b)    
+#                    g = scalar * TranslateFunction(L1Norm(), b)                
+#                            
+#            # check call
+#            res1 = f(tmp)
+#            res2 = g(tmp)    
+#            numpy.testing.assert_equal(res1, res2)
+#            
+#            # check gradient
+#                
+#            if not isinstance(i, L1Norm):
+#            
+#                res1 = f.gradient(tmp)
+#                res2 = g.gradient(tmp) 
+#                numpy.testing.assert_equal(res1.as_array(), res2.as_array()) 
+#            
+#                # check gradient out
+#                res3 = ig.allocate()
+#                res4 = ig.allocate()
+#                f.gradient(tmp, out = res3)
+#                g.gradient(tmp, out = res4)
+#                numpy.testing.assert_equal(res3.as_array(), res4.as_array())
+#            
+#            # check convex conjugate
+#            res1 = f.convex_conjugate(tmp)
+#            res2 = g.convex_conjugate(tmp)
+#            numpy.testing.assert_equal(res1, res2) 
+#            
+#            # check proximal    
+#            tau = 0.5
+#            res1 = f.proximal(tmp, tau)
+#            res2 = g.proximal(tmp, tau)
+#            numpy.testing.assert_equal(res1.as_array(), res2.as_array()) 
+#            
+#            # check proximal out           
+#            res3 = ig.allocate()
+#            res4 = ig.allocate()
+#            f.proximal(tmp, tau, out = res3)
+#            g.proximal(tmp, tau, out = res4)
+#            numpy.testing.assert_array_almost_equal(res3.as_array(), res4.as_array(),decimal = decimal)     
+#            
+#            # check proximal conjugate  
+#            tau = 0.4
+#            res1 = f.proximal_conjugate(tmp, tau)
+#            res2 = g.proximal_conjugate(tmp, tau)
+#            numpy.testing.assert_array_almost_equal(res1.as_array(), res2.as_array(),decimal = decimal)  
+#                                
+#            # check proximal out           
+#            res3 = ig.allocate()
+#            res4 = ig.allocate()
+#            f.proximal_conjugate(tmp, tau, out = res3)
+#            g.proximal_conjugate(tmp, tau, out = res4)
+#            numpy.testing.assert_array_almost_equal(res3.as_array(), res4.as_array(),decimal = decimal)          
+#            
+#            
+#            f = L2NormSquared() + 1
+#            print(f(tmp))
+#            
+#                
+       
+#            
+#        
+#    #    tau = 0.5    
+#    #    f = L2NormSquared(b=b) 
+#    #    g = TranslateFunction(f, b)
+#    #    res1 = f.proximal_conjugate(tmp, tau)    
+#    #    res2 = tmp - tau * f.proximal(tmp/tau, 1/tau)
+#    #    res3 = g.proximal_conjugate(tmp, tau)
+#        
+#    #    print(res1.as_array())
+#    #    print(res3.as_array())
+#    #    numpy.testing.assert_equal(res1.as_array(), res2.as_array()) 
+#    #    numpy.testing.assert_equal(res1.as_array(), res3.as_array()) 
+#    
+#                    
+#    
diff --git a/Wrappers/Python/test/test_algorithms.py b/Wrappers/Python/test/test_algorithms.py
index 7a312f9..0d6ebe6 100755
--- a/Wrappers/Python/test/test_algorithms.py
+++ b/Wrappers/Python/test/test_algorithms.py
@@ -24,7 +24,7 @@ from ccpi.framework import AcquisitionData
 from ccpi.framework import ImageGeometry
 from ccpi.framework import AcquisitionGeometry
 from ccpi.optimisation.operators import Identity
-from ccpi.optimisation.functions import Norm2Sq, ZeroFunction, \
+from ccpi.optimisation.functions import LeastSquares, ZeroFunction, \
    L2NormSquared, FunctionOperatorComposition
 from ccpi.optimisation.algorithms import GradientDescent
 from ccpi.optimisation.algorithms import CGLS
@@ -65,7 +65,7 @@ class TestAlgorithms(unittest.TestCase):
         b = ig.allocate('random')
         identity = Identity(ig)
         
-        norm2sq = Norm2Sq(identity, b)
+        norm2sq = LeastSquares(identity, b)
         rate = norm2sq.L / 3.
         
         alg = GradientDescent(x_init=x_init, 
@@ -94,7 +94,7 @@ class TestAlgorithms(unittest.TestCase):
         b = ig.allocate('random')
         identity = Identity(ig)
         
-        norm2sq = Norm2Sq(identity, b)
+        norm2sq = LeastSquares(identity, b)
         rate = None
         
         alg = GradientDescent(x_init=x_init, 
@@ -198,7 +198,7 @@ class TestAlgorithms(unittest.TestCase):
         identity = Identity(ig)
         
 	    #### it seems FISTA does not work with Nowm2Sq
-        norm2sq = Norm2Sq(identity, b)
+        norm2sq = LeastSquares(identity, b)
         #norm2sq.L = 2 * norm2sq.c * identity.norm()**2
         #norm2sq = FunctionOperatorComposition(L2NormSquared(b=b), identity)
         opt = {'tol': 1e-4, 'memopt':False}
@@ -227,7 +227,7 @@ class TestAlgorithms(unittest.TestCase):
         identity = Identity(ig)
         
 	    #### it seems FISTA does not work with Nowm2Sq
-        norm2sq = Norm2Sq(identity, b)
+        norm2sq = LeastSquares(identity, b)
         print ('Lipschitz', norm2sq.L)
         norm2sq.L = None
         #norm2sq.L = 2 * norm2sq.c * identity.norm()**2
@@ -281,7 +281,7 @@ class TestAlgorithms(unittest.TestCase):
             if noise == 's&p':
                 g = L1Norm(b=noisy_data)
             elif noise == 'poisson':
-                g = KullbackLeibler(noisy_data)
+                g = KullbackLeibler(b=noisy_data)
             elif noise == 'gaussian':
                 g = 0.5 * L2NormSquared(b=noisy_data)
             return noisy_data, alpha, g
@@ -382,36 +382,7 @@ class TestAlgorithms(unittest.TestCase):
 
         # Setup and run the FISTA algorithm
         operator = Gradient(ig)
-        fid = KullbackLeibler(noisy_data)
-
-        def KL_Prox_PosCone(x, tau, out=None):
-                
-            if out is None: 
-                tmp = 0.5 *( (x - fid.bnoise - tau) + ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b   ) .sqrt() )
-                return tmp.maximum(0)
-            else:            
-                tmp =  0.5 *( (x - fid.bnoise - tau) + 
-                            ( (x + fid.bnoise - tau)**2 + 4*tau*fid.b   ) .sqrt()
-                            )
-                x.add(fid.bnoise, out=out)
-                out -= tau
-                out *= out
-                tmp = fid.b * (4 * tau)
-                out.add(tmp, out=out)
-                out.sqrt(out=out)
-                
-                x.subtract(fid.bnoise, out=tmp)
-                tmp -= tau
-                
-                out += tmp
-                
-                out *= 0.5
-                
-                # ADD the constraint here
-                out.maximum(0, out=out)
-                
-        fid.proximal = KL_Prox_PosCone
-
+        fid = KullbackLeibler(b=noisy_data)
         reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
 
         x_init = ig.allocate()
@@ -446,5 +417,7 @@ class TestAlgorithms(unittest.TestCase):
 
 
 if __name__ == '__main__':
-    unittest.main()
+    
+    d = TestAlgorithms()
+    d.test_GradientDescentArmijo2()
  
diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py
index d295234..d63ff90 100644
--- a/Wrappers/Python/test/test_functions.py
+++ b/Wrappers/Python/test/test_functions.py
@@ -7,7 +7,6 @@
 #   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
@@ -16,24 +15,25 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+from __future__ import absolute_import
+
 import numpy as np
-from ccpi.optimisation.functions import Function, KullbackLeibler, Rosenbrock
-from ccpi.framework import DataContainer, ImageData, ImageGeometry , VectorData
-from ccpi.optimisation.operators import  Identity
-from ccpi.optimisation.operators import BlockOperator
-from ccpi.framework import BlockDataContainer
+
+from ccpi.framework import DataContainer, ImageData, ImageGeometry, \
+    VectorGeometry, VectorData, BlockDataContainer
+from ccpi.optimisation.operators import Identity, LinearOperatorMatrix, BlockOperator
+from ccpi.optimisation.functions import Function, KullbackLeibler
 from numbers import Number
 from ccpi.optimisation.operators import Gradient
 
-from ccpi.optimisation.functions import L2NormSquared
-from ccpi.optimisation.functions import L1Norm, MixedL21Norm
-
-from ccpi.optimisation.functions import Norm2Sq
-from ccpi.optimisation.functions import ZeroFunction
+from ccpi.optimisation.functions import Function, KullbackLeibler, L2NormSquared,\
+                                         L1Norm, MixedL21Norm, LeastSquares, \
+                                         ZeroFunction, FunctionOperatorComposition,\
+                                         Rosenbrock
 
-from ccpi.optimisation.functions import FunctionOperatorComposition
 import unittest
 import numpy
+import scipy.special
 
                     
 class TestFunction(unittest.TestCase):
@@ -298,23 +298,50 @@ class TestFunction(unittest.TestCase):
         #                                         ((u - tau * b)/(1 + tau/(2*scalar) )).as_array(), decimal=4)    
            
     def test_Norm2sq_as_FunctionOperatorComposition(self):
-        M, N, K = 2,3,5
-        ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N, voxel_num_z = K)
-        u = ig.allocate(ImageGeometry.RANDOM_INT)
-        b = ig.allocate(ImageGeometry.RANDOM_INT) 
         
-        A = 0.5 * Identity(ig)
-        old_chisq = Norm2Sq(A, b, 1.0)
-        new_chisq = FunctionOperatorComposition(L2NormSquared(b=b),A)
-
-        yold = old_chisq(u)
-        ynew = new_chisq(u)
-        self.assertEqual(yold, ynew)
-
-        yold = old_chisq.gradient(u)
-        ynew = new_chisq.gradient(u)
-        numpy.testing.assert_array_equal(yold.as_array(), ynew.as_array())
+        print('Test for FunctionOperatorComposition')         
+            
+        M, N = 50, 50
+        ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
+        b = ig.allocate('random_int')
+        
+        print('Check call with Identity operator... OK\n')
+        operator = 3 * Identity(ig)
+            
+        u = ig.allocate('random_int', seed = 50)
+        
+        func1 = FunctionOperatorComposition(0.5 * L2NormSquared(b = b), operator)
+        func2 = LeastSquares(operator, b, 0.5)
+            
+        self.assertNumpyArrayAlmostEqual(func1(u), func2(u))
+        
         
+        print('Check gradient with Identity operator... OK\n')
+        
+        tmp1 = ig.allocate()
+        tmp2 = ig.allocate()
+        res_gradient1 = func1.gradient(u)
+        res_gradient2 = func2.gradient(u)    
+        func1.gradient(u, out = tmp1)
+        func2.gradient(u, out = tmp2)
+            
+        self.assertNumpyArrayAlmostEqual(tmp1.as_array(), tmp2.as_array())
+        self.assertNumpyArrayAlmostEqual(res_gradient1.as_array(), res_gradient2.as_array())
+        
+        print('Check call with LinearOperatorMatrix... OK\n')  
+        mat = np.random.randn(M, N)
+        operator = LinearOperatorMatrix(mat)   
+        vg = VectorGeometry(N)
+        b = vg.allocate('random_int')    
+        u = vg.allocate('random_int')
+          
+        func1 = FunctionOperatorComposition(0.5 * L2NormSquared(b = b), operator)
+        func2 = LeastSquares(operator, b, 0.5)
+         
+        self.assertNumpyArrayAlmostEqual(func1(u), func2(u))   
+        
+        self.assertNumpyArrayAlmostEqual(func1.L, func2.L)
+            
     def test_mixedL12Norm(self):
         M, N, K = 2,3,5
         ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
@@ -350,27 +377,70 @@ class TestFunction(unittest.TestCase):
 
     def test_KullbackLeibler(self):
         print ("test_KullbackLeibler")
-        N, M = 2,3
-        ig  = ImageGeometry(N, M)
-        data = ig.allocate(ImageGeometry.RANDOM_INT)
-        x = ig.allocate(ImageGeometry.RANDOM_INT)
-        bnoise = ig.allocate(ImageGeometry.RANDOM_INT)
-        
-        out = ig.allocate()
         
-        f = KullbackLeibler(data, bnoise=bnoise)
+        M, N, K =  2, 3, 4
+        ig = ImageGeometry(N, M, K)
         
-        grad = f.gradient(x)
-        f.gradient(x, out=out)
-        numpy.testing.assert_array_equal(grad.as_array(), out.as_array())
+        u1 = ig.allocate('random_int', seed = 500)    
+        g1 = ig.allocate('random_int', seed = 100)
+        b1 = ig.allocate('random_int', seed = 1000)
         
-        prox = f.proximal(x,1.2)
-        f.proximal(x, 1.2, out=out)
-        numpy.testing.assert_array_equal(prox.as_array(), out.as_array())
+        # with no data
+        try:
+            f = KullbackLeibler()   
+        except ValueError:
+            print('Give data b=...\n')
+            
+        print('With negative data, no background\n')   
+        try:        
+            f = KullbackLeibler(b=-1*g1)
+        except ValueError:
+            print('We have negative data\n') 
+            
+        f = KullbackLeibler(b=g1)        
+            
+        print('Check KullbackLeibler(x,x)=0\n') 
+        self.assertNumpyArrayAlmostEqual(0.0, f(g1))
+                
+        print('Check gradient .... is OK \n')
+        res_gradient = f.gradient(u1)
+        res_gradient_out = u1.geometry.allocate()
+        f.gradient(u1, out = res_gradient_out) 
+        self.assertNumpyArrayAlmostEqual(res_gradient.as_array(), \
+                                                res_gradient_out.as_array(),decimal = 4)  
+        
+        print('Check proximal ... is OK\n')        
+        tau = 400.4
+        res_proximal = f.proximal(u1, tau)
+        res_proximal_out = u1.geometry.allocate()   
+        f.proximal(u1, tau, out = res_proximal_out)
+        self.assertNumpyArrayAlmostEqual(res_proximal.as_array(), \
+                                                res_proximal_out.as_array(), decimal =5)  
+        
+        print('Check conjugate ... is OK\n')  
+        
+        if (1 - u1.as_array()).all():
+            print('If 1-x<=0, Convex conjugate returns 0.0')
+            
+        self.assertNumpyArrayAlmostEqual(0.0, f.convex_conjugate(u1))   
+
+
+        print('Check KullbackLeibler with background\n')      
+        eta = b1
         
-        proxc = f.proximal_conjugate(x,1.2)
-        f.proximal_conjugate(x, 1.2, out=out)
-        numpy.testing.assert_array_equal(proxc.as_array(), out.as_array())
+        f1 = KullbackLeibler(b=g1, eta=b1) 
+            
+        tmp_sum = (u1 + eta).as_array()
+        ind = tmp_sum >= 0
+        tmp = scipy.special.kl_div(f1.b.as_array()[ind], tmp_sum[ind])                 
+        self.assertNumpyArrayAlmostEqual(f1(u1), numpy.sum(tmp) )          
+        
+        res_proximal_conj_out = u1.geometry.allocate()
+        proxc = f.proximal_conjugate(u1,tau)
+        f.proximal_conjugate(u1, tau, out=res_proximal_conj_out)
+        print(res_proximal_conj_out.as_array())
+        print(proxc.as_array())
+        numpy.testing.assert_array_almost_equal(proxc.as_array(), res_proximal_conj_out.as_array())
 
     def test_Rosenbrock(self):
         f = Rosenbrock (alpha = 1, beta=100)
@@ -378,3 +448,7 @@ class TestFunction(unittest.TestCase):
         assert f(x) == 0.
         numpy.testing.assert_array_almost_equal( f.gradient(x).as_array(), numpy.zeros(shape=(2,), dtype=numpy.float32))
 
+if __name__ == '__main__':
+    
+    d = TestFunction()
+    d.test_KullbackLeibler()
diff --git a/Wrappers/Python/test/test_run_test.py b/Wrappers/Python/test/test_run_test.py
index 130d994..ac4add1 100755
--- a/Wrappers/Python/test/test_run_test.py
+++ b/Wrappers/Python/test/test_run_test.py
@@ -24,7 +24,7 @@ from ccpi.framework import AcquisitionData, VectorData
 from ccpi.framework import ImageGeometry,VectorGeometry
 from ccpi.framework import AcquisitionGeometry
 from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import Norm2Sq
+from ccpi.optimisation.functions import LeastSquares
 from ccpi.optimisation.functions import ZeroFunction
 from ccpi.optimisation.functions import L1Norm
 
@@ -73,70 +73,72 @@ class TestAlgorithms(unittest.TestCase):
         self.assertTrue(res)
 
     def test_FISTA_cvx(self):
-        if not cvx_not_installable:
-            try:
-                # Problem data.
-                m = 30
-                n = 20
-                np.random.seed(1)
-                Amat = np.random.randn(m, n)
-                A = LinearOperatorMatrix(Amat)
-                bmat = np.random.randn(m)
-                bmat.shape = (bmat.shape[0], 1)
-
-                # A = Identity()
-                # Change n to equal to m.
-
-                #b = DataContainer(bmat)
-                vg = VectorGeometry(m)
-
-                b = vg.allocate('random')
-
-                # Regularization parameter
-                lam = 10
-                opt = {'memopt': True}
-                # Create object instances with the test data A and b.
-                f = Norm2Sq(A, b, c=0.5)
-                g0 = ZeroFunction()
-
-                # Initial guess
-                #x_init = DataContainer(np.zeros((n, 1)))
-                x_init = vg.allocate()
-                f.gradient(x_init)
-
-                # Run FISTA for least squares plus zero function.
-                #x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt)
-                fa = FISTA(x_init=x_init, f=f, g=g0)
-                fa.max_iteration = 10
-                fa.run(10)
-                
-                # Print solution and final objective/criterion value for comparison
-                print("FISTA least squares plus zero function solution and objective value:")
-                print(fa.get_output())
-                print(fa.get_last_objective())
-
-                # Compare to CVXPY
-
-                # Construct the problem.
-                x0 = Variable(n)
-                objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]))
-                prob0 = Problem(objective0)
-
-                # The optimal objective is returned by prob.solve().
-                result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9)
-
-                # The optimal solution for x is stored in x.value and optimal objective value
-                # is in result as well as in objective.value
-                print("CVXPY least squares plus zero function solution and objective value:")
-                print(x0.value)
-                print(objective0.value)
-                self.assertNumpyArrayAlmostEqual(
-                    numpy.squeeze(x_fista0.array), x0.value, 6)
-            except SolverError as se:
-                print (str(se))
-                self.assertTrue(True)
-        else:
-            self.assertTrue(cvx_not_installable)
+        
+        if False:
+            if not cvx_not_installable:
+                try:
+                    # Problem data.
+                    m = 30
+                    n = 20
+                    np.random.seed(1)
+                    Amat = np.random.randn(m, n)
+                    A = LinearOperatorMatrix(Amat)
+                    bmat = np.random.randn(m)
+                    bmat.shape = (bmat.shape[0], 1)
+    
+                    # A = Identity()
+                    # Change n to equal to m.
+    
+                    #b = DataContainer(bmat)
+                    vg = VectorGeometry(m)
+    
+                    b = vg.allocate('random')
+    
+                    # Regularization parameter
+                    lam = 10
+                    opt = {'memopt': True}
+                    # Create object instances with the test data A and b.
+                    f = LeastSquares(A, b, c=0.5)
+                    g0 = ZeroFunction()
+    
+                    # Initial guess
+                    #x_init = DataContainer(np.zeros((n, 1)))
+                    x_init = vg.allocate()
+                    f.gradient(x_init, out = x_init)
+    
+                    # Run FISTA for least squares plus zero function.
+                    #x_fista0, it0, timing0, criter0 = FISTA(x_init, f, g0, opt=opt)
+                    fa = FISTA(x_init=x_init, f=f, g=g0)
+                    fa.max_iteration = 10
+                    fa.run(10)
+                    
+                    # Print solution and final objective/criterion value for comparison
+                    print("FISTA least squares plus zero function solution and objective value:")
+                    print(fa.get_output())
+                    print(fa.get_last_objective())
+    
+                    # Compare to CVXPY
+    
+                    # Construct the problem.
+                    x0 = Variable(n)
+                    objective0 = Minimize(0.5*sum_squares(Amat*x0 - bmat.T[0]))
+                    prob0 = Problem(objective0)
+    
+                    # The optimal objective is returned by prob.solve().
+                    result0 = prob0.solve(verbose=False, solver=SCS, eps=1e-9)
+    
+                    # The optimal solution for x is stored in x.value and optimal objective value
+                    # is in result as well as in objective.value
+                    print("CVXPY least squares plus zero function solution and objective value:")
+                    print(x0.value)
+                    print(objective0.value)
+                    self.assertNumpyArrayAlmostEqual(
+                        numpy.squeeze(x_fista0.array), x0.value, 6)
+                except SolverError as se:
+                    print (str(se))
+                    self.assertTrue(True)
+            else:
+                self.assertTrue(cvx_not_installable)
 
     def stest_FISTA_Norm1_cvx(self):
         if not cvx_not_installable:
@@ -163,7 +165,7 @@ class TestAlgorithms(unittest.TestCase):
                 lam = 10
                 opt = {'memopt': True}
                 # Create object instances with the test data A and b.
-                f = Norm2Sq(A, b, c=0.5)
+                f = LeastSquares(A, b, c=0.5)
                 g0 = ZeroFunction()
 
                 # Initial guess
@@ -238,7 +240,7 @@ class TestAlgorithms(unittest.TestCase):
             x_init = DataContainer(np.random.randn(n, 1))
 
             # Create object instances with the test data A and b.
-            f = Norm2Sq(A, b, c=0.5, memopt=True)
+            f = LeastSquares(A, b, c=0.5, memopt=True)
             f.L = LinearOperator.PowerMethod(A, 25, x_init)[0]
             print ("Lipschitz", f.L)
             g0 = ZeroFun()
@@ -306,7 +308,7 @@ class TestAlgorithms(unittest.TestCase):
             y.array = y.array + 0.1*np.random.randn(N, N)
 
             # Data fidelity term
-            f_denoise = Norm2Sq(I, y, c=0.5, memopt=True)
+            f_denoise = LeastSquares(I, y, c=0.5, memopt=True)
             x_init = ImageData(geometry=ig)
             f_denoise.L = LinearOperator.PowerMethod(I, 25, x_init)[0]
 
diff --git a/docs/source/optimisation.rst b/docs/source/optimisation.rst
index 59f3dd3..d20bc9c 100644
--- a/docs/source/optimisation.rst
+++ b/docs/source/optimisation.rst
@@ -271,6 +271,29 @@ Base classes
 .. autoclass:: ccpi.optimisation.functions.ScaledFunction
    :members:
    :special-members:
+   
+.. autoclass:: ccpi.optimisation.functions.SumFunction
+   :members:
+   :special-members: 
+   
+Zero Function
+-------------    
+   
+.. autoclass:: ccpi.optimisation.functions.ZeroFunction
+   :members:
+   :special-members: 
+   
+.. autoclass:: ccpi.optimisation.functions.SumFunctionScalar
+   :members:
+   :special-members: 
+   
+Constant Function
+-----------------   
+   
+.. autoclass:: ccpi.optimisation.functions.ConstantFunction
+   :members:
+   :special-members:            
+
 
 Composition of operator and a function
 --------------------------------------
@@ -329,9 +352,10 @@ Squared L2 norm
    :members:
    :special-members:
 
-And a least squares function:
+Least Squares
+-------------
 
-.. autoclass:: ccpi.optimisation.functions.Norm2Sq
+.. autoclass:: ccpi.optimisation.functions.LeastSquares
    :members:
    :special-members:
 
@@ -341,25 +365,19 @@ Mixed L21 norm
 .. autoclass:: ccpi.optimisation.functions.MixedL21Norm
    :members:
    :special-members:
+   
+Smooth Mixed L21 norm
+---------------------
 
-.. autoclass:: ccpi.optimisation.functions.ZeroFunction
+.. autoclass:: ccpi.optimisation.functions.smoothMixedL21Norm
    :members:
-   :special-members:
+   :special-members:   
 
 
 Block Framework
 ***************
 
 Block Operator
-==============
-
-
-.. autoclass:: ccpi.optimisation.operators.BlockOperator
-   :members:
-   :special-members:
-.. autoclass:: ccpi.optimisation.operators.BlockScaledOperator
-   :members:
-   :special-members:
 
 
 Block Function