Fix documentation (#343)

* fix documentation algs

* fix docum functions dir

*  fix docum operator dir

* docstring title

* updated math in docstring

29 files changed, 833 insertions, 552 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index 1168db3..fec37c5 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
@@ -1,20 +1,26 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+#========================================================================
+# 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.
+#
+#=========================================================================
 
-#   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.
 import time, functools
 from numbers import Integral
 
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
index 82fbb96..c6c1d4c 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/CGLS.py
@@ -1,40 +1,47 @@
 # -*- 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.
+# 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 ccpi.optimisation.algorithms import Algorithm
-from ccpi.optimisation.functions import Norm2Sq
 import numpy
 
 class CGLS(Algorithm):
 
-    '''Conjugate Gradient Least Squares algorithm
+    r'''Conjugate Gradient Least Squares algorithm 
+    
+    Problem:  
+
+    .. math::
 
-    Parameters:
-      x_init: initial guess
-      operator: operator for forward/backward projections
-      data: data to operate on
-      tolerance: tolerance to stop algorithm
+      \min || A x - b ||^2_2
+    
+    |
+
+    Parameters :
+        
+      :parameter operator : Linear operator for the inverse problem
+      :parameter x_init : Initial guess ( Default x_init = 0)
+      :parameter data : Acquired data to reconstruct       
+      :parameter tolerance: Tolerance/ Stopping Criterion to end CGLS algorithm
       
     Reference:
         https://web.stanford.edu/group/SOL/software/cgls/
-      
     '''
-
     def __init__(self, **kwargs):
         
         super(CGLS, self).__init__()
@@ -59,9 +66,7 @@ class CGLS(Algorithm):
         self.p = self.s
         self.norms0 = self.s.norm()
         
-        ##
         self.norms = self.s.norm()
-        ##
 
         self.gamma = self.norms0**2
         self.normx = self.x.norm()
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index fa1d8d8..347dacd 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -1,43 +1,60 @@
 # -*- 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.
+# 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 ccpi.optimisation.algorithms import Algorithm
 from ccpi.optimisation.functions import ZeroFunction
 import numpy
 
 class FISTA(Algorithm):
-    '''Fast Iterative Shrinkage-Thresholding Algorithm
     
-    Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding 
-    algorithm for linear inverse problems. 
-    SIAM journal on imaging sciences,2(1), pp.183-202.
+    r'''Fast Iterative Shrinkage-Thresholding Algorithm 
+    
+    Problem : 
+    
+    .. math::
+    
+      \min_{x} f(x) + g(x)
+    
+    |
     
-    Parameters:
-      x_init: initial guess
-      f: data fidelity
-      g: regularizer
-      opt: additional options 
+    Parameters :
+        
+      :parameter x_init : Initial guess ( Default x_init = 0)
+      :parameter f : Differentiable function
+      :parameter g : Convex function with " simple " proximal operator
+
+
+    Reference:
+      
+        Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding 
+        algorithm for linear inverse problems. 
+        SIAM journal on imaging sciences,2(1), pp.183-202.
     '''
     
     
     def __init__(self, **kwargs):
-        '''initialisation can be done at creation time if all 
+        
+        '''creator 
+        
+        initialisation can be done at creation time if all 
         proper variables are passed or later with set_up'''
+        
         super(FISTA, self).__init__()
         f = kwargs.get('f', None)
         g = kwargs.get('g', ZeroFunction())
@@ -67,12 +84,11 @@ class FISTA(Algorithm):
         self.f.gradient(self.y, out=self.u)
         self.u.__imul__( -self.invL )
         self.u.__iadd__( self.y )
-        # x = g.prox(u,invL)
+
         self.g.proximal(self.u, self.invL, out=self.x)
         
         self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2)))
         
-#        self.x.subtract(self.x_old, out=self.y)
         self.y = self.x - self.x_old
         self.y.__imul__ ((self.t_old-1)/self.t)
         self.y.__iadd__( self.x )
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
index 8c2b693..d060690 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/GradientDescent.py
@@ -1,24 +1,33 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+#========================================================================
+# 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.
+#
+#=========================================================================
 
-#   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.algorithms import Algorithm
 
 class GradientDescent(Algorithm):
-    '''Implementation of Gradient Descent algorithm
+    ''' 
+    
+        Gradient Descent algorithm
+        
+        
     '''
 
     def __init__(self, **kwargs):
@@ -35,7 +44,7 @@ class GradientDescent(Algorithm):
             self.set_up(x_init=x_init, objective_function=objective_function, rate=rate)
     
     def should_stop(self):
-        '''stopping cryterion, currently only based on number of iterations'''
+        '''stopping criterion, currently only based on number of iterations'''
         return self.iteration >= self.max_iteration
     
     def set_up(self, x_init, objective_function, rate):
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 8f37765..0968872 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -1,32 +1,60 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+# 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.
 
-#   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.algorithms import Algorithm
-from ccpi.framework import ImageData, DataContainer
-import numpy as np
-import numpy
-import time
-from ccpi.optimisation.operators import BlockOperator
-from ccpi.framework import BlockDataContainer
-from ccpi.optimisation.functions import FunctionOperatorComposition
+
 
 
 class PDHG(Algorithm):
-    '''Primal Dual Hybrid Gradient'''
+    r'''Primal Dual Hybrid Gradient
+    
+    Problem: 
+    
+    .. math::
+    
+      \min_{x} f(Kx) + g(x)
+    |
+
+    Parameters : 
+        
+        :parameter operator : Linear Operator = K
+        :parameter f : Convex function with "simple" proximal of its conjugate. 
+        :parameter g : Convex function with "simple" proximal 
+        :parameter sigma : Step size parameter for Primal problem
+        :parameter tau : Step size parameter for Dual problem
+        
+        Remark: Convergence is guaranted provided that
+        
+        .. math:: \tau \sigma \|K\|^{2} <1
+        
+            
+    Reference :
+        
+        
+        (a) A. Chambolle and T. Pock (2011), "A first-order primal–dual algorithm for convex
+        problems with applications to imaging", J. Math. Imaging Vision 40, 120–145.        
+        
+        
+        (b) E. Esser, X. Zhang and T. F. Chan (2010), "A general framework for a class of first
+        order primal–dual algorithms for convex optimization in imaging science",
+        SIAM J. Imaging Sci. 3, 1015–1046.
+    '''
 
     def __init__(self, **kwargs):
         super(PDHG, self).__init__(max_iteration=kwargs.get('max_iteration',0))
@@ -78,7 +106,7 @@ class PDHG(Algorithm):
 
     def update(self):
         
-        # Gradient descent, Dual problem solution
+        # Gradient ascent for the dual variable
         self.operator.direct(self.xbar, out=self.y_tmp)
         self.y_tmp *= self.sigma
         self.y_tmp += self.y_old
@@ -86,7 +114,7 @@ class PDHG(Algorithm):
         # self.y = self.f.proximal_conjugate(self.y_old, self.sigma)
         self.f.proximal_conjugate(self.y_tmp, self.sigma, out=self.y)
         
-        # Gradient ascent, Primal problem solution
+        # Gradient descent for the primal variable
         self.operator.adjoint(self.y, out=self.x_tmp)
         self.x_tmp *= -1*self.tau
         self.x_tmp += self.x_old
@@ -98,6 +126,7 @@ class PDHG(Algorithm):
         self.xbar *= self.theta
         self.xbar += self.x
 
+        
         self.x_old.fill(self.x)
         self.y_old.fill(self.y)
 
@@ -106,4 +135,4 @@ class PDHG(Algorithm):
         p1 = self.f(self.operator.direct(self.x)) + self.g(self.x)
         d1 = -(self.f.convex_conjugate(self.y) + self.g.convex_conjugate(-1*self.operator.adjoint(self.y)))
 
-        self.loss.append([p1, d1, p1-d1])
\ No newline at end of file
+        self.loss.append([p1, d1, p1-d1])
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
index ca5b084..2b49ab0 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/SIRT.py
@@ -1,33 +1,46 @@
 # -*- 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.
+#========================================================================
+# 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 ccpi.optimisation.algorithms import Algorithm
 
 class SIRT(Algorithm):
 
-    '''Simultaneous Iterative Reconstruction Technique
+    r'''Simultaneous Iterative Reconstruction Technique
+    
+    Problem: 
+    
+    .. math::  
+    
+    A x = b
+    |
 
     Parameters:
-      x_init: initial guess
-      operator: operator for forward/backward projections
-      data: data to operate on
-      constraint: Function with prox-method, for example IndicatorBox to 
-                  enforce box constraints, default is None).
+        
+      :parameter operator : Linear operator for the inverse problem
+      :parameter x_init : Initial guess
+      :parameter data : Acquired data to reconstruct       
+      :parameter constraint : Function proximal method
+                   e.g.  x\in[0, 1], IndicatorBox to enforce box constraints
+                         Default is None).
     '''
     def __init__(self, **kwargs):
         super(SIRT, self).__init__()
diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
index 6d4b4e1..a6ac66c 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
@@ -1,20 +1,24 @@
 # -*- 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.
+#========================================================================
+# 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 ccpi.optimisation.functions import Function
 from ccpi.framework import BlockDataContainer
@@ -22,29 +26,29 @@ from numbers import Number
 
 class BlockFunction(Function):
     
-    '''BlockFunction acts as a separable sum function, i.e.,
+    r'''BlockFunction acts as a separable sum function: f = [f_1,...,f_n]
     
-      f = [f_1,...,f_n]
-      
-      f([x_1,...,x_n]) = f_1(x_1) +  .... + f_n(x_n)
+      .. math::
+
+          f([x_1,...,x_n]) = f_1(x_1) +  .... + f_n(x_n)
+      |
 
     '''
+    
     def __init__(self, *functions):
                 
         super(BlockFunction, self).__init__()
         self.functions = functions      
         self.length = len(self.functions)
-        
-        
-        
-        
+                                
     def __call__(self, x):
         
-        '''Evaluates the BlockFunction at a BlockDataContainer x
+        r'''Evaluates the BlockFunction at a BlockDataContainer x
         
-        :param: x (BlockDataContainer): must have as many rows as self.length
+            :param: x (BlockDataContainer): must have as many rows as self.length
 
-        returns sum(f_i(x_i))
+            returns ..math:: sum(f_i(x_i))
+            
         '''
         
         if self.length != x.shape[0]:
@@ -56,9 +60,9 @@ class BlockFunction(Function):
     
     def convex_conjugate(self, x):
         
-        ''' Evaluate convex conjugate of BlockFunction at x
+        r'''Convex conjugate of BlockFunction at x            
         
-        returns sum(f_i^{*}(x_i))
+            .. math:: returns sum(f_i^{*}(x_i))
         
         '''       
         t = 0                
@@ -69,9 +73,9 @@ class BlockFunction(Function):
     
     def proximal_conjugate(self, x, tau, out = None):
         
-        ''' Evaluate Proximal Operator of tau * f(\cdot) at x
-        
-        prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) 
+        r'''Proximal operator of BlockFunction at x: 
+                 
+                 .. math:: prox_{tau*f}(x) = sum_{i} prox_{tau*f_{i}}(x_{i}) 
         
         
         '''
@@ -99,11 +103,9 @@ class BlockFunction(Function):
     
     def proximal(self, x, tau, out = None):
         
-        ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x
-        
-        prox_{tau*f^{*}}(x) = sum_{i} prox_{tau*f^{*}_{i}}(x_{i}) 
-        
+        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 out is None:
@@ -130,9 +132,9 @@ class BlockFunction(Function):
     
     def gradient(self,x, out=None):
         
-        ''' Evaluate gradient of f at x: f'(x) 
+        r'''Evaluates gradient of BlockFunction at x
         
-        returns: BlockDataContainer [f_{1}'(x_{1}), ... , f_{n}'(x_{n})]
+            returns: BlockDataContainer .. math:: [f_{1}'(x_{1}), ... , f_{n}'(x_{n})]
                 
         '''
         
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Function.py b/Wrappers/Python/ccpi/optimisation/functions/Function.py
index 9a6009f..7156995 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/Function.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/Function.py
@@ -1,20 +1,24 @@
 # -*- 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.
+#========================================================================
+# 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.
+#
+#=========================================================================
 
 import warnings
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction
diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
index 24ebbf0..58d4f27 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
@@ -1,30 +1,35 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+#========================================================================
+# 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.
+#
+#=========================================================================
 
-#   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
 from ccpi.optimisation.functions import ScaledFunction
 
 
 class FunctionOperatorComposition(Function):
     
-    ''' Function composition with Operator, i.e., f(Ax)
+    '''Function composition with Operator: (f o A)(x) = f(Ax)
     
-        A: operator
-        f: function
+            : parameter A: operator
+            : parameter f: function
     
     '''
     
@@ -39,17 +44,18 @@ class FunctionOperatorComposition(Function):
         
     def __call__(self, x):
         
-        ''' Evaluate FunctionOperatorComposition at x
-        
-        returns f(Ax)
-        
-        '''
+        '''Evaluates f(Ax)'''
     
         return self.function(self.operator.direct(x))  
     
     def gradient(self, x, out=None):
-#        
-        ''' Gradient takes into account the Operator'''
+        
+        '''Evaluates gradient of f(Ax):
+        
+        ..math ::  A^{T}f'(Ax)
+            
+        '''
+        
         if out is None:
             return self.operator.adjoint(self.function.gradient(self.operator.direct(x)))
         else: 
diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
index 2001a56..51d08d1 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
@@ -18,24 +18,31 @@
 
 from ccpi.optimisation.functions import Function
 import numpy
-from ccpi.framework import ImageData
 
 class IndicatorBox(Function):
-    '''Box constraints indicator function. 
     
-    Calling returns 0 if argument is within the box. The prox operator is projection onto the box. 
-    Only implements one scalar lower and one upper as constraint on all elements. Should generalise
-    to vectors to allow different constraints one elements.
-'''
+    
+    r'''Indicator function for box constraint
+            
+      .. math:: 
+         f(x) = \mathbb{I}_{[a, b]} = \begin{cases}
+            
+                                            0, if x\in[a, b]
+                                            \infty, otherwise                            
+                                    \end{cases}
+    
+    '''
     
     def __init__(self,lower=-numpy.inf,upper=numpy.inf):
-        # Do nothing
+
         super(IndicatorBox, self).__init__()
         self.lower = lower
         self.upper = upper
         
-    
+
     def __call__(self,x):
+        
+        '''Evaluates IndicatorBox at x'''
                 
         if (numpy.all(x.array>=self.lower) and 
             numpy.all(x.array <= self.upper) ):
@@ -48,12 +55,18 @@ class IndicatorBox(Function):
         return ValueError('Not Differentiable') 
     
     def convex_conjugate(self,x):
-        # support function sup <x, z>, z \in [lower, upper]
-        # ????
+        
+        '''Convex conjugate of IndicatorBox at x'''
+
         return x.maximum(0).sum()
          
     def proximal(self, x, tau, out=None):
         
+        r'''Proximal operator of IndicatorBox at x
+
+            .. math:: prox_{\tau * f}(x)
+        '''
+        
         if out is None:
             return (x.maximum(self.lower)).minimum(self.upper)        
         else:               
@@ -61,6 +74,11 @@ class IndicatorBox(Function):
             out.minimum(self.upper, out=out) 
             
     def proximal_conjugate(self, x, tau, out=None):
+        
+        r'''Proximal operator of the convex conjugate of IndicatorBox at x:
+
+          ..math:: prox_{\tau * f^{*}}(x)
+        '''
 
         if out is None:
             
diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index 2a87205..f88c339 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -1,20 +1,22 @@
 # -*- 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
+# 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.
 
-#   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.
 import numpy
 from ccpi.optimisation.functions import Function
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction 
@@ -23,13 +25,14 @@ import scipy.special
 
 class KullbackLeibler(Function):
     
-    ''' 
-    
-    KL_div(x, y + back) = int x * log(x/(y+back)) - x + (y+back)
+    r'''Kullback-Leibler divergence function
     
-    Assumption: y>=0
-                back>=0
-                
+         .. 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}
+            
     '''
     
     def __init__(self, data, **kwargs):
@@ -42,29 +45,24 @@ class KullbackLeibler(Function):
                                                     
     def __call__(self, x):
         
-        
-        '''
-        
-            x - y * log( x + bnoise) + y * log(y) - y + bnoise
-        
-        
-        '''
-        
+
+        '''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) 
-          
 
     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):
+        if not functools.reduce(lambda x,y: x and y>0, datacontainer.as_array().ravel(), True):
             raise ValueError('KullbackLeibler. Cannot calculate log of negative number')
         datacontainer.fill( numpy.log(datacontainer.as_array()) )
 
         
     def gradient(self, x, out=None):
         
+        '''Evaluates gradient of KullbackLeibler at x'''
+        
         if out is None:
             return 1 - self.b/(x + self.bnoise)
         else:
@@ -76,11 +74,19 @@ class KullbackLeibler(Function):
             
     def convex_conjugate(self, x):
         
+        '''Convex conjugate of KullbackLeibler at x'''
+        
         xlogy = - scipy.special.xlogy(self.b.as_array(), 1 - x.as_array())
         return numpy.sum(xlogy)
             
     def proximal(self, x, tau, out=None):
         
+        r'''Proximal operator of KullbackLeibler at x
+           
+           .. math::     prox_{\tau * f}(x)
+
+        '''
+
         if out is None:        
             return 0.5 *( (x - self.bnoise - tau) + ( (x + self.bnoise - tau)**2 + 4*tau*self.b   ) .sqrt() )
         else:
@@ -103,6 +109,11 @@ class KullbackLeibler(Function):
             out *= 0.5
                             
     def proximal_conjugate(self, x, tau, out=None):
+        
+        r'''Proximal operator of the convex conjugate of KullbackLeibler at x:
+           
+           .. math::     prox_{\tau * f^{*}}(x)
+        '''
 
                 
         if out is None:
@@ -110,7 +121,6 @@ class KullbackLeibler(Function):
             return 0.5*((z + 1) - ((z-1)**2 + 4 * tau * self.b).sqrt())
         else:
             
-            #tmp = x + tau * self.bnoise
             tmp = tau * self.bnoise
             tmp += x
             tmp -= 1
@@ -126,10 +136,9 @@ class KullbackLeibler(Function):
 
     def __rmul__(self, scalar):
         
-        ''' Multiplication of L2NormSquared with a scalar
-        
-        Returns: ScaledFunction
-                        
+        '''Multiplication of KullbackLeibler with a scalar        
+            
+            Returns: ScaledFunction
         '''
         
         return ScaledFunction(self, scalar) 
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
index 08c3d17..97d03b9 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
@@ -1,20 +1,21 @@
 # -*- 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.
+# 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 ccpi.optimisation.functions import Function
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction        
@@ -23,16 +24,14 @@ from ccpi.optimisation.operators import ShrinkageOperator
 
 class L1Norm(Function):
     
-    ''' 
-    
-    Class: L1Norm
-        
-    Cases:  a) f(x) = ||x||_{1}
-    
-            b) f(x) = ||x - b||_{1}
-                             
-    '''      
-    
+    r'''L1Norm function: 
+            
+            Cases considered (with/without data):            
+                a) .. math:: f(x) = ||x||_{1}
+                b) .. math:: f(x) = ||x - b||_{1}
+                                
+    '''   
+           
     def __init__(self, **kwargs):
         
         super(L1Norm, self).__init__()
@@ -40,7 +39,7 @@ class L1Norm(Function):
         
     def __call__(self, x):
         
-        ''' Evaluate L1Norm at x: f(x) '''
+        '''Evaluates L1Norm at x'''
         
         y = x
         if self.b is not None: 
@@ -48,11 +47,12 @@ class L1Norm(Function):
         return y.abs().sum()  
     
     def gradient(self,x):
-        #TODO implement subgradient???
+        
         return ValueError('Not Differentiable')   
     
     def convex_conjugate(self,x):
-        #TODO implement Indicator infty???
+        
+        '''Convex conjugate of L1Norm at x'''
 
         y = 0        
         if self.b is not None:
@@ -61,7 +61,11 @@ class L1Norm(Function):
     
     def proximal(self, x, tau, out=None):
         
-        # TODO implement shrinkage operator, we will need it later e.g SplitBregman
+        r'''Proximal operator of L1Norm at x
+           
+           ..math::     prox_{\tau * f}(x)
+                
+        ''' 
         
         if out is None:
             if self.b is not None:
@@ -76,6 +80,12 @@ class L1Norm(Function):
                                     
     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))
@@ -88,6 +98,12 @@ class L1Norm(Function):
                 out.fill(x.divide(x.abs().maximum(1.0)) )                
             
     def __rmul__(self, scalar):
+        
+        '''Multiplication of L2NormSquared with a scalar        
+            
+            Returns: ScaledFunction
+        '''
+        
         return ScaledFunction(self, scalar)
 
 
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
index 2afae25..a625f07 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
@@ -1,44 +1,47 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+# 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.
 
-#   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
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction
-from ccpi.optimisation.functions import FunctionOperatorComposition
 
 class L2NormSquared(Function):
     
-    '''
+    r'''L2NormSquared function: 
             
-    Cases:  a) f(x) = \|x\|^{2}_{2}
-    
-            b) f(x) = ||x - b||^{2}_{2}
-                             
+            Cases considered (with/without data):            
+                a) .. math:: f(x) = \|x\|^{2}_{2} 
+                b) .. math:: f(x) = \|\|x - b\|\|^{2}_{2}
+                                
     '''    
     
     def __init__(self, **kwargs):
                                 
         super(L2NormSquared, self).__init__()
-        self.b = kwargs.get('b',None)  
+        self.b = kwargs.get('b',None) 
+        
+        # Lipschitz constant
         self.L = 2
 
     def __call__(self, x):
         
-        ''' Evaluate L2NormSquared at x: f(x) '''
+        '''Evaluates L2NormSquared at x'''
             
         y = x
         if self.b is not None: 
@@ -51,7 +54,7 @@ class L2NormSquared(Function):
                 
     def gradient(self, x, out=None):        
         
-        ''' Evaluate gradient of L2NormSquared at x: f'(x) '''
+        '''Evaluates gradient of L2NormSquared at x'''
                 
         if out is not None:
             
@@ -70,7 +73,7 @@ class L2NormSquared(Function):
                                                        
     def convex_conjugate(self, x):
         
-        ''' Evaluate convex conjugate of L2NormSquared at x: f^{*}(x)'''
+        '''Convex conjugate of L2NormSquared at x'''
         
         tmp = 0
         
@@ -82,11 +85,10 @@ class L2NormSquared(Function):
 
     def proximal(self, x, tau, out = None):
 
-        ''' Evaluate Proximal Operator of tau * f(\cdot) at x:
-        
-            prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f(z)
-        
-        '''        
+        r'''Proximal operator of L2NormSquared at x
+            
+            .. math:: prox_{\tau * f}(x)
+        '''          
         
         if out is None:
             
@@ -109,11 +111,9 @@ class L2NormSquared(Function):
     
     def proximal_conjugate(self, x, tau, out=None):
         
-        ''' Evaluate Proximal Operator of tau * f^{*}(\cdot) at x (i.e., the convex conjugate of f) :
-        
-            prox_{tau*f(\cdot)}(x) = \argmin_{z} \frac{1}{2}|| z - x ||^{2}_{2} + tau * f^{*}(z)
-        
-        '''         
+        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:
@@ -129,20 +129,13 @@ class L2NormSquared(Function):
                                         
     def __rmul__(self, scalar):
         
-        ''' Multiplication of L2NormSquared with a scalar
-        
-        Returns: ScaledFunction
-                        
-        '''
+        '''Multiplication of L2NormSquared with a scalar        
+           
+           Returns: ScaledFunction'''
         
         return ScaledFunction(self, scalar)  
 
 
-    def composition(self, operator):
-        
-        return FunctionOperatorComposition(operator)
-      
-
 
 if __name__ == '__main__':
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index d7d23a5..8cbed67 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -1,20 +1,22 @@
 # -*- 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
+# 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.
 
-#   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, ScaledFunction
 from ccpi.framework import BlockDataContainer
@@ -24,8 +26,10 @@ import functools
 class MixedL21Norm(Function):
     
     
-    '''
-        f(x) = ||x||_{2,1} = \sum |x|_{2}                   
+    r'''MixedL21Norm: .. math:: f(x) = ||x||_{2,1} = \int \|x\|_{2} dx
+
+        where x is a vector/tensor vield
+                
     '''      
     
     def __init__(self, **kwargs):
@@ -35,10 +39,9 @@ class MixedL21Norm(Function):
         
     def __call__(self, x):
         
-        ''' Evaluates L2,1Norm at point x
+        '''Evaluates MixedL21Norm at point x
             
-            :param: x is a BlockDataContainer
-                                
+           :param: x: is a BlockDataContainer
         '''
         if not isinstance(x, BlockDataContainer):
             raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x))) 
@@ -53,26 +56,31 @@ class MixedL21Norm(Function):
                             
     def convex_conjugate(self,x):
         
-        ''' This is the Indicator function of ||\cdot||_{2, \infty}
-            which is either 0 if ||x||_{2, \infty} or \infty        
+        r'''Convex conjugate of of MixedL21Norm: 
+        
+        Indicator function of .. math:: ||\cdot||_{2, \infty}
+            which is either 0 if .. math:: ||x||_{2, \infty}<1 or \infty 
+            
         '''
         
         return 0.0
         
-        #tmp = [ el**2 for el in x.containers ]
-        #print(sum(tmp).sqrt().as_array().max())
-        #return sum(tmp).sqrt().as_array().max()
     
     def proximal(self, x, tau, out=None):
         
-        '''
-            For this we need to define a MixedL2,2 norm acting on BDC,
-            different form L2NormSquared which acts on DC
-        
+        r'''Proximal operator of MixedL21Norm at x:
+           
+           .. math:: prox_{\tau * f(x)
         '''
         pass
     
     def proximal_conjugate(self, x, tau, out=None): 
+        
+        r'''Proximal operator of the convex conjugate of MixedL21Norm at x:
+           
+           .. math:: prox_{\tau * f^{*}}(x)
+
+        '''           
 
 
         if out is None:                                        
@@ -80,26 +88,19 @@ class MixedL21Norm(Function):
             res = sum(tmp).sqrt().maximum(1.0) 
             frac = [el/res for el in x.containers]
             return  BlockDataContainer(*frac)   
-        
-            
-        #TODO this is slow, why???
-#                return x.divide(x.pnorm().maximum(1.0))
+                
         else:
                             
             res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 )
             res = res1.sqrt().maximum(1.0)
             x.divide(res, out=out)
-            
-#                x.divide(sum([el*el for el in x.containers]).sqrt().maximum(1.0), out=out)
-            #TODO this is slow, why ???
-#                 x.divide(x.pnorm().maximum(1.0), out=out)
-                              
+                                          
 
     def __rmul__(self, scalar):
         
-        ''' Multiplication of L2NormSquared with a scalar
-        
-        Returns: ScaledFunction
+        '''Multiplication of MixedL21Norm with a scalar        
+            
+            Returns: ScaledFunction
              
         '''         
         return ScaledFunction(self, scalar) 
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
index eea300d..c8f5e46 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
@@ -1,36 +1,37 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+# 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.
 
-#   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
-import numpy
 import warnings
 
 # Define a class for squared 2-norm
 class Norm2Sq(Function):
-    '''
-    f(x) = c*||A*x-b||_2^2
+    r'''.. math:: f(x) = c*||A*x-b||_2^2
     
     which has 
     
-    grad[f](x) = 2*c*A^T*(A*x-b)
+    .. math:: grad[f](x) = 2*c*A^T*(A*x-b)
     
     and Lipschitz constant
     
-    L = 2*c*||A||_2^2 = 2*s1(A)^2
+    .. math:: L = 2*c*||A||_2^2 = 2*s1(A)^2
     
     where s1(A) is the largest singular value of A.
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
index df008e0..3e689e6 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/ScaledFunction.py
@@ -1,20 +1,22 @@
 # -*- coding: utf-8 -*-

-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).

+# 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.

 

-#   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 numbers import Number

 import numpy

 import warnings

diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py
index fb990f8..ca52f31 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFunction.py
@@ -1,53 +1,81 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+# 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.
 
-#   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
-from ccpi.framework import BlockDataContainer 
 
 class ZeroFunction(Function):
     
-    ''' ZeroFunction: f(x) = 0
-    
-    
+    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):
         
-        ''' This is the support function sup <x, x^{*}>  which in fact is the 
-        indicator function for the set = {0}
-        So 0 if x=0, or inf if x neq 0                
+        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:
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 5dfd6ba..e3a02ec 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -1,37 +1,37 @@
 # -*- 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
+# 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.
 
-#   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.
 import numpy
-from numbers import Number
 import functools
-from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer, DataContainer
-from ccpi.optimisation.operators import Operator, LinearOperator
-from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator
+from ccpi.framework import ImageData, BlockDataContainer, DataContainer
+from ccpi.optimisation.operators import Operator
 from ccpi.framework import BlockGeometry
        
 class BlockOperator(Operator):
-    '''A Block matrix containing Operators
+    r'''A Block matrix containing Operators
 
     The Block Framework is a generic strategy to treat variational problems in the
     following form:
 
     .. math::
     
-      min Regulariser + Fidelity
+      \min Regulariser + Fidelity
 
     
     BlockOperators have a generic shape M x N, and when applied on an 
@@ -192,7 +192,6 @@ class BlockOperator(Operator):
             else:
                 return BlockDataContainer(*res, shape=shape)
         else:
-            #tmp = self.domain_geometry().allocate()
 
             for col in range(self.shape[1]):
                 for row in range(self.shape[0]):
@@ -240,7 +239,7 @@ class BlockOperator(Operator):
     def __rmul__(self, scalar):
         '''Defines the left multiplication with a scalar
 
-        Args: scalar (number or iterable containing numbers):
+        :paramer scalar: (number or iterable containing numbers):
 
         Returns: a block operator with Scaled Operators inside'''
         if isinstance (scalar, list) or isinstance(scalar, tuple) or \
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
index 30b750d..c23c23a 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
@@ -1,44 +1,51 @@
 # -*- coding: utf-8 -*-
-#  CCP in Tomographic Imaging (CCPi) Core Imaging Library (CIL).
+# 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.
 
-#   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 numbers import Number
 import numpy
 from ccpi.optimisation.operators import ScaledOperator
 import functools
 
 class BlockScaledOperator(ScaledOperator):
+    
     '''ScaledOperator
 
-    A class to represent the scalar multiplication of an Operator with a scalar.
-    It holds an operator and a scalar. Basically it returns the multiplication
-    of the result of direct and adjoint of the operator with the scalar.
-    For the rest it behaves like the operator it holds.
+        A class to represent the scalar multiplication of an Operator with a scalar.
+        It holds an operator and a scalar. Basically it returns the multiplication
+        of the result of direct and adjoint of the operator with the scalar.
+        For the rest it behaves like the operator it holds.
+    
+        Args:
+           :param operator (Operator): a Operator or LinearOperator
+           :param scalar (Number): a scalar multiplier
+        Example:
+           The scaled operator behaves like the following:
+   .. code-block:: python
 
-    Args:
-       operator (Operator): a Operator or LinearOperator
-       scalar (Number): a scalar multiplier
-    Example:
-       The scaled operator behaves like the following:
-       sop = ScaledOperator(operator, scalar)
-       sop.direct(x) = scalar * operator.direct(x)
-       sop.adjoint(x) = scalar * operator.adjoint(x)
-       sop.norm() = operator.norm()
-       sop.range_geometry() = operator.range_geometry()
-       sop.domain_geometry() = operator.domain_geometry()
+        sop = ScaledOperator(operator, scalar)
+        sop.direct(x) = scalar * operator.direct(x)
+        sop.adjoint(x) = scalar * operator.adjoint(x)
+        sop.norm() = operator.norm()
+        sop.range_geometry() = operator.range_geometry()
+        sop.domain_geometry() = operator.domain_geometry()
+                      
     '''
     def __init__(self, operator, scalar, shape=None):
         if shape is None:
@@ -75,10 +82,4 @@ class BlockScaledOperator(ScaledOperator):
     @property
     def T(self):
         '''Return the transposed of self'''
-        #print ("transpose before" , self.shape)
-        #shape = (self.shape[1], self.shape[0])
-        ##self.shape = shape
-        ##self.operator.shape = shape
-        #print ("transpose" , shape)
-        #return self
         return type(self)(self.operator.T, self.scalar)
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index cb8cb65..9b5ae24 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -16,25 +16,31 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 from ccpi.optimisation.operators import LinearOperator
-from ccpi.framework import ImageData, BlockDataContainer
 import numpy as np
 
 class FiniteDiff(LinearOperator):
     
-    # Works for Neum/Symmetric &  periodic boundary conditions
-    # TODO add central differences???
-    # TODO not very well optimised, too many conditions
-    # TODO add discretisation step, should get that from imageGeometry
-    
-    # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-    # Grad_order = ['channels', 'direction_y', 'direction_x']
-    # Grad_order = ['direction_z', 'direction_y', 'direction_x']
-    # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
+    '''Finite Difference Operator:
+            
+            Computes first-order forward/backward differences 
+                     on 2D, 3D, 4D ImageData
+                     under Neumann/Periodic boundary conditions
+
+        Order of the Gradient ( ImageGeometry may contain channels ):
+                            
+            Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
+            Grad_order = ['channels', 'direction_y', 'direction_x']
+            Grad_order = ['direction_z', 'direction_y', 'direction_x']
+            Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']  
+                                        
     
+    '''
+
+        
     def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
-        ''''''
+
         super(FiniteDiff, self).__init__() 
-        '''FIXME: domain and range should be geometries'''
+
         self.gm_domain = gm_domain
         self.gm_range = gm_range
         
@@ -44,6 +50,7 @@ class FiniteDiff(LinearOperator):
         # Domain Geometry = Range Geometry if not stated
         if self.gm_range is None:
             self.gm_range = self.gm_domain
+            
         # check direction and "length" of geometry
         if self.direction + 1 > len(self.gm_domain.shape):
             raise ValueError('Gradient directions more than geometry domain')      
@@ -63,9 +70,7 @@ class FiniteDiff(LinearOperator):
         else:
             out = out.as_array()
             out[:]=0
-
                   
-
         ######################## Direct for 2D  ###############################
         if x_sz == 2:
             
@@ -317,11 +322,22 @@ class FiniteDiff(LinearOperator):
         return type(x)(out)
             
     def range_geometry(self):
-        '''Returns the range geometry'''
+        
+        '''
+        
+            Returns the range_geometry of FiniteDiff
+        
+        '''
+        
         return self.gm_range
     
     def domain_geometry(self):
-        '''Returns the domain geometry'''
+        
+        '''
+        
+            Returns the domain_geometry of FiniteDiff
+        
+        '''        
         return self.gm_domain
 
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index b3b23bf..baebc61 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -23,12 +23,29 @@ from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 #%%
 
 class Gradient(LinearOperator):
+    
+    
+    r'''Gradient Operator: .. math:: \nabla : X -> Y           
+            
+            Computes first-order forward/backward differences 
+                     on 2D, 3D, 4D ImageData
+                     under Neumann/Periodic boundary conditions
+                                                             
+                Example (2D): u\in X, \nabla(u) = [\partial_{y} u, \partial_{x} u]
+                              u^{*}\in Y, \nabla^{*}(u^{*}) = \partial_{y} v1 + \partial_{x} v2
+        
+                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
+                Grad_order = ['channels', 'direction_y', 'direction_x']
+                Grad_order = ['direction_z', 'direction_y', 'direction_x']
+                Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
+                
+
+    '''
+    
+                    
     CORRELATION_SPACE = "Space"
     CORRELATION_SPACECHANNEL = "SpaceChannels"
-    # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
-    # Grad_order = ['channels', 'direction_y', 'direction_x']
-    # Grad_order = ['direction_z', 'direction_y', 'direction_x']
-    # Grad_order = ['channels', 'direction_z', 'direction_y', 'direction_x']
+
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
         super(Gradient, self).__init__() 
@@ -72,8 +89,7 @@ class Gradient(LinearOperator):
          
         self.bnd_cond = bnd_cond 
         
-        # Call FiniteDiff operator
-        
+        # Call FiniteDiff operator        
         self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond)
                                                          
         
@@ -114,12 +130,24 @@ class Gradient(LinearOperator):
             
     
     def domain_geometry(self):
+        
+        '''Returns domain_geometry of Gradient'''
+        
         return self.gm_domain
     
     def range_geometry(self):
+        
+        '''Returns range_geometry of Gradient'''
+        
         return self.gm_range
     
     def __rmul__(self, scalar):
+        
+        '''Multiplication of Gradient with a scalar        
+            
+            Returns: ScaledOperator
+        '''        
+        
         return ScaledOperator(self, scalar) 
     
     ###########################################################################
diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index 54dce08..d8f86a4 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -19,11 +19,18 @@
 from ccpi.optimisation.operators import LinearOperator
 import scipy.sparse as sp
 import numpy as np
-from ccpi.framework import ImageData
 
 
 class Identity(LinearOperator):
     
+    '''Identity:  Id: X -> Y,  Id(x) = x\in Y
+                       
+                   X : gm_domain
+                   Y : gm_range ( Default: Y = X )
+                                                                                                                                   
+    '''    
+    
+    
     def __init__(self, gm_domain, gm_range=None):
 
         self.gm_domain = gm_domain
@@ -34,38 +41,58 @@ class Identity(LinearOperator):
         super(Identity, self).__init__()
         
     def direct(self,x,out=None):
+        
+        '''Returns Id(x)'''
+        
         if out is None:
             return x.copy()
         else:
             out.fill(x)
     
     def adjoint(self,x, out=None):
+        
+        '''Returns Id(x)'''         
+        
+        
         if out is None:
             return x.copy()
         else:
             out.fill(x)
         
     def calculate_norm(self, **kwargs):
+        
+        '''Evaluates operator norm of Identity'''        
+        
         return 1.0
         
-    def domain_geometry(self):       
+    def domain_geometry(self): 
+        
+        '''Returns domain_geometry of Identity'''
+        
         return self.gm_domain
         
     def range_geometry(self):
+        
+        '''Returns range_geometry of Identity'''         
+        
         return self.gm_range
     
+    
+    ###########################################################################
+    ###############  For preconditioning ######################################
+    ###########################################################################                    
     def matrix(self):
         
         return sp.eye(np.prod(self.gm_domain.shape))
     
     def sum_abs_row(self):
         
-        return self.gm_range.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
- 
+        return self.gm_range.allocate(1)
+    
     def sum_abs_col(self):
         
-        return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
-            
+        return self.gm_domain.allocate(1)
+    
     
 if __name__ == '__main__':
     
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index b254d08..8514699 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -17,7 +17,6 @@
 #   limitations under the License.

 

 from ccpi.optimisation.operators import Operator

-from ccpi.framework import ImageGeometry

 import numpy

 

 

diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
index fbc3c83..bc3312d 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperatorMatrix.py
@@ -17,14 +17,12 @@
 #   limitations under the License.
 import numpy
 from scipy.sparse.linalg import svds
-from ccpi.framework import DataContainer
-from ccpi.framework import AcquisitionData
-from ccpi.framework import VectorData
 from ccpi.framework import VectorGeometry
-from ccpi.framework import AcquisitionGeometry
-from numbers import Number
 from ccpi.optimisation.operators import LinearOperator
+
 class LinearOperatorMatrix(LinearOperator):
+    '''Matrix wrapped into a LinearOperator'''
+    
     def __init__(self,A):
         self.A = A
         M_A, N_A = self.A.shape
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
index b9ebef1..c5db47d 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
@@ -19,23 +19,33 @@ from numbers import Number
 import numpy
 
 class ScaledOperator(object):
+    
+    
     '''ScaledOperator
+
     A class to represent the scalar multiplication of an Operator with a scalar.
     It holds an operator and a scalar. Basically it returns the multiplication
     of the result of direct and adjoint of the operator with the scalar.
     For the rest it behaves like the operator it holds.
+    
     Args:
-       operator (Operator): a Operator or LinearOperator
-       scalar (Number): a scalar multiplier
+       :param operator (Operator): a Operator or LinearOperator
+       :param scalar (Number): a scalar multiplier
+    
     Example:
        The scaled operator behaves like the following:
-       sop = ScaledOperator(operator, scalar)
-       sop.direct(x) = scalar * operator.direct(x)
-       sop.adjoint(x) = scalar * operator.adjoint(x)
-       sop.norm() = operator.norm()
-       sop.range_geometry() = operator.range_geometry()
-       sop.domain_geometry() = operator.domain_geometry()
+
+    .. code-block:: python
+
+      sop = ScaledOperator(operator, scalar)
+      sop.direct(x) = scalar * operator.direct(x)
+      sop.adjoint(x) = scalar * operator.adjoint(x)
+      sop.norm() = operator.norm()
+      sop.range_geometry() = operator.range_geometry()
+      sop.domain_geometry() = operator.domain_geometry()
+
     '''
+    
     def __init__(self, operator, scalar):
         super(ScaledOperator, self).__init__()
         if not isinstance (scalar, Number):
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
index 1a9f1d8..c1f7ca4 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
@@ -20,6 +20,11 @@ from ccpi.framework import DataContainer
 
 class ShrinkageOperator():
     
+    r'''Proximal Operator for .. math:: f(x) = \|\| x \|\|_{1}
+        
+            prox_{\tau * f}(x) = x.sign() * \max( |x| - \tau, 0 )
+    '''
+    
     def __init__(self):
         pass
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 6b4d822..91d5ca9 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -22,6 +22,10 @@ from ccpi.framework import ImageData
 
 class SparseFiniteDiff(object):
     
+    
+    '''Create Sparse Matrices for the Finite Difference Operator'''
+    
+    
     def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
         
         super(SparseFiniteDiff, self).__init__() 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index bed56b6..92f8f90 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -19,25 +19,34 @@
 from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator
 from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
-from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
+from ccpi.optimisation.operators import FiniteDiff
 
 
 class SymmetrizedGradient(Gradient):
     
-    ''' Symmetrized Gradient, denoted by E: V --> W
-        where V is the Range of the Gradient Operator
-        and W is the Range of the Symmetrized Gradient.                        
+    r'''Symmetrized Gradient Operator:  E: V -> W
+        
+            V : range of the Gradient Operator
+            W : range of the Symmetrized Gradient          
+        
+            Example (2D): 
+            .. math::
+               v = (v1, v2), 
+            
+                           Ev = 0.5 * ( \nabla\cdot v + (\nabla\cdot c)^{T} )
+                           
+                           \begin{matrix} 
+                               \partial_{y} v1 & 0.5 * (\partial_{x} v1 + \partial_{y} v2) \\
+                               0.5 * (\partial_{x} v1 + \partial_{y} v2) & \partial_{x} v2 
+                           \end{matrix}
+              |                                                      
     '''
     
     
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
         super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) 
-        
-        '''
-         Domain of SymGrad is the Range of Gradient
-        '''
-        
+                
         self.gm_domain = self.gm_range 
         self.bnd_cond = bnd_cond
         
@@ -57,6 +66,8 @@ class SymmetrizedGradient(Gradient):
         
     def direct(self, x, out=None):
         
+        '''Returns E(v)'''        
+        
         if out is None:
             
             tmp = []
@@ -232,12 +243,4 @@ if __name__ == '__main__':
     print(E1.norm())
     print(E2.norm())
     
-    
-    
-    
-    
-
-#    
-#    
-#    
-#   
+  
diff --git a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
index b9345da..5f1de30 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ZeroOperator.py
@@ -22,6 +22,22 @@ from ccpi.optimisation.operators import LinearOperator
 
 class ZeroOperator(LinearOperator):
     
+    r'''ZeroOperator:  O: X -> Y,  maps any element of x\in X into the zero element in Y
+                       O(x) = O_{Y}
+                       
+                       X : gm_domain
+                       Y : gm_range ( Default: Y = X )
+                       
+                       
+                       Note: 
+                       .. math::
+
+                              O^{*}: Y^{*} -> X^{*} (Adjoint)
+                       
+                              < O(x), y > = < x, O^{*}(y) >
+                       
+     '''
+    
     def __init__(self, gm_domain, gm_range=None):
         
         super(ZeroOperator, self).__init__()             
@@ -33,22 +49,39 @@ class ZeroOperator(LinearOperator):
                    
         
     def direct(self,x,out=None):
+        
+        '''Returns O(x)'''
+        
+        
         if out is None:
             return self.gm_range.allocate()
         else:
             out.fill(self.gm_range.allocate())
     
     def adjoint(self,x, out=None):
+        
+        '''Returns O^{*}(y)'''        
+        
         if out is None:
             return self.gm_domain.allocate()
         else:
             out.fill(self.gm_domain.allocate())
         
     def calculate_norm(self, **kwargs):
-        return 0.
+        
+        '''Evaluates operator norm of ZeroOperator'''
+        
+        return 0
     
-    def domain_geometry(self):       
+    def domain_geometry(self): 
+        
+        '''Returns domain_geometry of ZeroOperator'''
+        
+        
         return self.gm_domain  
         
     def range_geometry(self):
+        
+        '''Returns domain_geometry of ZeroOperator'''
+        
         return self.gm_range
\ No newline at end of file