6 files changed, 225 insertions, 20 deletions
diff --git a/NOTICE.txt b/NOTICE.txt
new file mode 100644
index 0000000..c107329
--- /dev/null
+++ b/NOTICE.txt
@@ -0,0 +1,15 @@
+CCPi Core Imaging Library (CIL).
+Copyright 2017 Rutherford Appleton Laboratory STFC
+Copyright 2017 University of Manchester
+
+This software product is developed for the Collaborative Computational 
+Project in Tomographic Imaging CCPi (http://www.ccpi.ac.uk/) at RAL STFC (http://www.stfc.ac.uk), University of Manchester
+and other contributing institutions.
+
+Main contributors in alphabetical order:
+Evelina Ametova (UoM)
+Jakob Jorgensen (UoM)
+Daniil Kazantsev (Diamond Light Source)
+Srikanth Nagella (STFC)
+Edoardo Pasca (STFC)
+Evangelos Papoutsellis (UoM)
diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index ffc91ae..7236e0e 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -764,12 +764,26 @@ class DataContainer(object):
         return numpy.sqrt(self.squared_norm())
     def dot(self, other, *args, **kwargs):
         '''return the inner product of 2 DataContainers viewed as vectors'''
+        method = kwargs.get('method', 'reduce')
+        if method not in ['numpy','reduce']:
+            raise ValueError('dot: specified method not valid. Expecting numpy or reduce got {} '.format(
+                    method))
         if self.shape == other.shape:
-            return numpy.dot(self.as_array().ravel(), other.as_array().ravel())
+            # return (self*other).sum()
+            if method == 'numpy':
+                return numpy.dot(self.as_array().ravel(), other.as_array())
+            elif method == 'reduce':
+                # see https://github.com/vais-ral/CCPi-Framework/pull/273
+                # notice that Python seems to be smart enough to use
+                # the appropriate type to hold the result of the reduction
+                sf = reduce(lambda x,y: x + y[0]*y[1],
+                            zip(self.as_array().ravel(),
+                                other.as_array().ravel()),
+                            0)
+                return sf
         else:
             raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape))
-    
-    
+   
     
     
     
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
index ed95c3f..78fe196 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/Algorithm.py
@@ -16,7 +16,7 @@
 #   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
+import time, functools
 from numbers import Integral
 
 class Algorithm(object):
@@ -34,7 +34,7 @@ class Algorithm(object):
       method will stop when the stopping cryterion is met. 
    '''
 
-    def __init__(self):
+    def __init__(self, **kwargs):
         '''Constructor
         
         Set the minimal number of parameters:
@@ -48,11 +48,11 @@ class Algorithm(object):
                                        when evaluating the objective is computationally expensive.
         '''
         self.iteration = 0
-        self.__max_iteration = 0
+        self.__max_iteration = kwargs.get('max_iteration', 0)
         self.__loss = []
         self.memopt = False
         self.timing = []
-        self.update_objective_interval = 1
+        self.update_objective_interval = kwargs.get('update_objective_interval', 1)
     def set_up(self, *args, **kwargs):
         '''Set up the algorithm'''
         raise NotImplementedError()
@@ -91,9 +91,11 @@ class Algorithm(object):
             if self.iteration % self.update_objective_interval == 0:
                 self.update_objective()
             self.iteration += 1
+        
     def get_output(self):
         '''Returns the solution found'''
         return self.x
+    
     def get_last_loss(self):
         '''Returns the last stored value of the loss function
         
@@ -146,13 +148,60 @@ class Algorithm(object):
             print ("Stop cryterion has been reached.")
         i = 0
         for _ in self:
-            if verbose and self.iteration % self.update_objective_interval == 0:
-                print ("Iteration {}/{}, objective {}".format(self.iteration, 
-                       self.max_iteration, self.get_last_objective()) )
-            else:
+            if i == 0 and verbose:
+                print (self.verbose_header())
+            if (self.iteration -1) % self.update_objective_interval == 0:                
+                if verbose:
+                    #print ("Iteration {:>7} max: {:>7}, = {}".format(self.iteration-1, 
+                    #   self.max_iteration, self.get_last_objective()) )
+                    print (self.verbose_output())
                 if callback is not None:
-                    callback(self.iteration, self.get_last_objective())
+                    callback(self.iteration -1, self.get_last_objective(), self.x)
             i += 1
             if i == iterations:
                 break
-    
+    def verbose_output(self):
+        '''Creates a nice tabulated output'''
+        timing = self.timing[-self.update_objective_interval-1:-1]
+        if len (timing) == 0:
+            t = 0
+        else:
+            t = sum(timing)/len(timing)
+        el = [ self.iteration-1, 
+               self.max_iteration,
+               "{:.3f}".format(t), 
+               self.get_last_objective() ]
+        
+        string = self.objective_to_string()
+        out = "{:>9} {:>10} {:>13} {}".format(*el[:-1] , string)
+        return out
+
+    def objective_to_string(self):
+        el = self.get_last_objective()
+        if type(el) == list:
+            string = functools.reduce(lambda x,y: x+' {:>13.5e}'.format(y), el[:-1],'')
+            string += '{:>15.5e}'.format(el[-1])
+        else:
+            string = "{:>20.5e}".format(el)
+        return string
+    def verbose_header(self):
+        el = self.get_last_objective()
+        if type(el) == list:
+            out = "{:>9} {:>10} {:>13} {:>13} {:>13} {:>15}\n".format('Iter', 
+                                                      'Max Iter',
+                                                      'Time/Iter',
+                                                      'Primal' , 'Dual', 'Primal-Dual')
+            out += "{:>9} {:>10} {:>13} {:>13} {:>13} {:>15}".format('', 
+                                                      '',
+                                                      '[s]',
+                                                      'Objective' , 'Objective', 'Gap')
+        else:
+            out = "{:>9} {:>10} {:>13} {:>20}\n".format('Iter', 
+                                                      'Max Iter',
+                                                      'Time/Iter',
+                                                      'Objective')
+            out += "{:>9} {:>10} {:>13} {:>20}".format('', 
+                                                      '',
+                                                      '[s]',
+                                                      '')
+        return out
+\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index 89b5519..f5ba85e 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -20,13 +20,8 @@
 import numpy
 import time
 
-from ccpi.optimisation.functions import Function
-from ccpi.optimisation.functions import ZeroFunction
-from ccpi.framework import ImageData 
-from ccpi.framework import AcquisitionData
-from ccpi.optimisation.spdhg import spdhg 
-from ccpi.optimisation.spdhg import KullbackLeibler
-from ccpi.optimisation.spdhg import KullbackLeiblerConvexConjugate
+
+
 
 def FISTA(x_init, f=None, g=None, opt=None):
     '''Fast Iterative Shrinkage-Thresholding Algorithm
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 8e8ab87..e92d4c6 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -455,6 +455,11 @@ class TestDataContainer(unittest.TestCase):
             self.assertTrue(False)
         except ValueError as ve:
             self.assertTrue(True)
+            
+        print ("test dot numpy")
+        n0 = (ds0 * ds1).sum()
+        n1 = ds0.as_array().ravel().dot(ds1.as_array().ravel())
+        self.assertEqual(n0, n1)
         
         
 
diff --git a/Wrappers/Python/wip/compare_CGLS_algos.py b/Wrappers/Python/wip/compare_CGLS_algos.py
new file mode 100644
index 0000000..119752c
--- /dev/null
+++ b/Wrappers/Python/wip/compare_CGLS_algos.py
@@ -0,0 +1,127 @@
+# This demo illustrates how to use the SIRT algorithm without and with 
+# nonnegativity and box constraints. The ASTRA 2D projectors are used.
+
+# First make all imports
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, \
+    AcquisitionData
+from ccpi.optimisation.algs import FISTA, FBPD, CGLS, SIRT
+from ccpi.astra.operators import AstraProjectorSimple
+
+from ccpi.optimisation.algorithms import CGLS as CGLSalg
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+# Choose either a parallel-beam (1=parallel2D) or fan-beam (2=cone2D) test case
+test_case = 1
+
+# Set up phantom size NxN by creating ImageGeometry, initialising the 
+# ImageData object with this geometry and empty array and finally put some
+# data into its array, and display as image.
+N = 128
+ig = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=ig)
+
+x = Phantom.as_array()
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+#plt.figure()
+#plt.imshow(x)
+#plt.title('Phantom image')
+#plt.show()
+
+# Set up AcquisitionGeometry object to hold the parameters of the measurement
+# setup geometry: # Number of angles, the actual angles from 0 to 
+# pi for parallel beam and 0 to 2pi for fanbeam, set the width of a detector 
+# pixel relative to an object pixel, the number of detector pixels, and the 
+# source-origin and origin-detector distance (here the origin-detector distance 
+# set to 0 to simulate a "virtual detector" with same detector pixel size as
+# object pixel size).
+angles_num = 20
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+if test_case==1:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('parallel',
+                             '2D',
+                             angles,
+                             det_num,det_w)
+elif test_case==2:
+    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+    ag = AcquisitionGeometry('cone',
+                             '2D',
+                             angles,
+                             det_num,
+                             det_w,
+                             dist_source_center=SourceOrig, 
+                             dist_center_detector=OrigDetec)
+else:
+    NotImplemented
+
+# Set up Operator object combining the ImageGeometry and AcquisitionGeometry
+# wrapping calls to ASTRA as well as specifying whether to use CPU or GPU.
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
+
+# Forward and backprojection are available as methods direct and adjoint. Here 
+# generate test data b and do simple backprojection to obtain z.
+b = Aop.direct(Phantom)
+z = Aop.adjoint(b)
+
+#plt.figure()
+#plt.imshow(b.array)
+#plt.title('Simulated data')
+#plt.show()
+
+#plt.figure()
+#plt.imshow(z.array)
+#plt.title('Backprojected data')
+#plt.show()
+
+# Using the test data b, different reconstruction methods can now be set up as
+# demonstrated in the rest of this file. In general all methods need an initial 
+# guess and some algorithm options to be set:
+x_init = ImageData(np.zeros(x.shape),geometry=ig)
+opt = {'tol': 1e-4, 'iter': 7}
+
+# First a CGLS reconstruction using the function version of CGLS can be done:
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt)
+
+#plt.figure()
+#plt.imshow(x_CGLS.array)
+#plt.title('CGLS')
+#plt.colorbar()
+#plt.show()
+
+#plt.figure()
+#plt.semilogy(criter_CGLS)
+#plt.title('CGLS criterion')
+#plt.show()
+
+
+
+# Now CLGS using the algorithm class
+CGLS_alg = CGLSalg()
+CGLS_alg.set_up(x_init, Aop, b )
+CGLS_alg.max_iteration = 2000
+CGLS_alg.run(opt['iter'])
+x_CGLS_alg = CGLS_alg.get_output()
+
+#plt.figure()
+#plt.imshow(x_CGLS_alg.as_array())
+#plt.title('CGLS ALG')
+#plt.colorbar()
+#plt.show()
+
+#plt.figure()
+#plt.semilogy(CGLS_alg.objective)
+#plt.title('CGLS criterion')
+#plt.show()
+
+print(criter_CGLS)
+print(CGLS_alg.objective)
+
+print((x_CGLS - x_CGLS_alg).norm())
+\ No newline at end of file