Merge pull request #90 from vais-ral/doc_opti

First opti doc

3 files changed, 65 insertions, 21 deletions

diff --git a/README.md b/README.md
index 4a14723..06d30fe 100644
--- a/README.md
+++ b/README.md
@@ -41,10 +41,54 @@ In `ccpi.framework` we define a number of common classes normally used in tomogr
  
  ### `ccpi.optimisation`
  
- This package allows writing of optimisation algorithms. The main actors here are:
+ This package allows rapid prototyping of optimisation-based reconstruction problems, 
+ i.e. defining and solving different optimization problems to enforce different properties 
+ on the reconstructed image.
  
-  * `Function`
-  * `Operator`
+ Firstly, it provides an object-oriented framework for defining mathematical operators and functions 
+ as well a collection of useful example operators and functions. Both smooth and 
+ non-smooth functions can be used. 
+ 
+ Further, it provides a number of high-level generic 
+ implementations of optimisation algorithms to solve genericlly formulated 
+ optimisation problems constructed from operator and function objects. 
+ 
+ The fundamental components are:
+ 
+  * `Operator`: A class specifying a (currently linear) operator
+  * `Function`: A class specifying mathematical functions such as a least squares data fidelity.
+  * `Algorithm`: Implementation of an optimisation algorithm to solve a particular generic optimisation problem. These are currently python functions by may be changed to operators in another release.
+
+ #### `Operator`
+ 
+ The two most important methods are `direct` and `adjoint` methods that describe the result of 
+ applying the operator, and its adjoint respectively, onto a compatible `DataContainer` input. 
+ The output is another `DataContainer` object or subclass hereof. An important 
+ special case is to represent the tomographic forward and backprojection operations.
+ 
+ #### `Function`
+ 
+ A `function` represents a mathematical function of one or more inputs and is intended 
+ to accept `DataContainer`s as input as well as any additional parameters. 
+ Fixed parameters can be passed in during the creation of the `function` object.
+ The methods of the `function` reflect the properties of it, for example, 
+ if the function represented is differentiable 
+ the `function` should contain a method `grad` which should return the gradient of the function evaluated at
+ an input point. If the function is not differentiable but allows a simple proximal operator, the method 
+ `prox` should return the proxial operator evaluated at an input point. The function value 
+ is evaluated by calling the function itself, e.g. `f(x)` for a `function` 
+ `f` and input point `x`.
+ 
+ #### `Algorithm`
+ 
+ A number of generic algorithm implementations are provided including CGLS and FISTA. An algorithm 
+ is designed for a particular generic optimisation problem accepts and number of `function`s and/or 
+ `operator`s as input to define a specific instance of the generic optimisation problem to be solved.
+ 
+ #### Examples
+ 
+ Please see the demos for examples of defining and using operators, functions and algorithms 
+ to specify and solve optimisation-based reconstruction problems.
  
  
  
diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index 5942055..a45100c 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -20,7 +20,7 @@
 import numpy
 import time
 
-from ccpi.optimisation.funcs import BaseFunction
+from ccpi.optimisation.funcs import Function
 
 def FISTA(x_init, f=None, g=None, opt=None):
     '''Fast Iterative Shrinkage-Thresholding Algorithm
@@ -37,8 +37,8 @@ def FISTA(x_init, f=None, g=None, opt=None):
       opt: additional algorithm 
     '''
     # default inputs
-    if f   is None: f = BaseFunction()
-    if g   is None: g = BaseFunction()
+    if f   is None: f = Function()
+    if g   is None: g = Function()
     
     # algorithmic parameters
     if opt is None: 
@@ -84,7 +84,7 @@ def FISTA(x_init, f=None, g=None, opt=None):
         
         # time and criterion
         timing[it] = time.time() - time0
-        criter[it] = f.fun(x) + g.fun(x);
+        criter[it] = f(x) + g(x);
         
         # stopping rule
         #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
@@ -108,9 +108,9 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None):
       opt: additional algorithm 
     '''
     # default inputs
-    if f   is None: f = BaseFunction()
-    if g   is None: g = BaseFunction()
-    if h   is None: h = BaseFunction()
+    if f   is None: f = Function()
+    if g   is None: g = Function()
+    if h   is None: h = Function()
     
     # algorithmic parameters
     if opt is None: 
@@ -156,7 +156,7 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None):
 
         # time and criterion
         timing[it] = time.time() - t
-        criter[it] = f.fun(x) + g.fun(x) + h.fun(h.op.direct(x));
+        criter[it] = f(x) + g(x) + h(h.op.direct(x));
            
         # stopping rule
         #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py
index 0dbe28f..d11d6c3 100755
--- a/Wrappers/Python/ccpi/optimisation/funcs.py
+++ b/Wrappers/Python/ccpi/optimisation/funcs.py
@@ -21,14 +21,14 @@ from ccpi.optimisation.ops import Identity, FiniteDiff2D
 import numpy
 
 
-class BaseFunction(object):
+class Function(object):
     def __init__(self):
         self.op = Identity()
-    def fun(self,x):      return 0
+    def __call__(self,x):      return 0
     def grad(self,x):     return 0
     def prox(self,x,tau): return x
 
-class Norm2(BaseFunction):
+class Norm2(Function):
     
     def __init__(self, 
                  gamma=1.0, 
@@ -37,7 +37,7 @@ class Norm2(BaseFunction):
         self.gamma     = gamma;
         self.direction = direction; 
     
-    def fun(self, x):
+    def __call__(self, x):
         
         xx = numpy.sqrt(numpy.sum(numpy.square(x.as_array()), self.direction,
                                   keepdims=True))
@@ -63,7 +63,7 @@ class TV2D(Norm2):
         
 
 # Define a class for squared 2-norm
-class Norm2sq(BaseFunction):
+class Norm2sq(Function):
     '''
     f(x) = c*||A*x-b||_2^2
     
@@ -93,19 +93,19 @@ class Norm2sq(BaseFunction):
         #return 2*self.c*self.A.adjoint( self.A.direct(x) - self.b )
         return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b )
     
-    def fun(self,x):
+    def __call__(self,x):
         #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel()))
         return self.c*( ( (self.A.direct(x)-self.b)**2).sum() )
 
 
-class ZeroFun(BaseFunction):
+class ZeroFun(Function):
     
     def __init__(self,gamma=0,L=1):
         self.gamma = gamma
         self.L = L
         super(ZeroFun, self).__init__()
     
-    def fun(self,x):
+    def __call__(self,x):
         return 0
     
     def prox(self,x,tau):
@@ -113,7 +113,7 @@ class ZeroFun(BaseFunction):
 
 # A more interesting example, least squares plus 1-norm minimization.
 # Define class to represent 1-norm including prox function
-class Norm1(BaseFunction):
+class Norm1(Function):
     
     def __init__(self,gamma):
         # Do nothing
@@ -121,7 +121,7 @@ class Norm1(BaseFunction):
         self.L = 1
         super(Norm1, self).__init__()
     
-    def fun(self,x):
+    def __call__(self,x):
         return self.gamma*(x.abs().sum())
     
     def prox(self,x,tau):