Geometry allocation (#211)

* initial revision

* Removed class members of Algorithm class

added update_objective

* initial version. Fix inline __idiv__

* First implementation of CompositeOperator/DataContainer

* removed __getitem__ added get_item added shape

* added CGLS

* working unit test, initial tomography test

* added reverse multiplication of operator with number

* added operators directory

* fixed typo

* added unittest for CompositeDataContainer

* fix TomoIdentity with scalar

* check numerical types from numpy

* add default stop criterion and run method

* add run method

* first working implementation of CGLS with CompositeOperator/DataContainer

notice problem with _rmul_ and _mul_ methods precedence with numpy.

* removed line endings

* removed dos line ending

* Added allocate to ImageGeometry and AcquisitionGeometry

* remove composite operator

* Delete Algorithms.py

* remove setup.py from PR

* readded setup.py

* added newline at the end of the file

* added newline at EOF in setup.py

4 files changed, 176 insertions, 122 deletions

diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 0c23628..71a9f3f 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -111,8 +111,12 @@ class ImageGeometry(object):
         repres += "voxel_size : x{0},y{1},z{2}\n".format(self.voxel_size_x, self.voxel_size_y, self.voxel_size_z)
         repres += "center : x{0},y{1},z{2}\n".format(self.center_x, self.center_y, self.center_z)
         return repres
-    
-    
+    def allocate(self, value=0, dimension_labels=None):
+        '''allocates an ImageData according to the size expressed in the instance'''
+        out = ImageData(geometry=self, dimension_labels=dimension_labels)
+        if value != 0:
+            out += value
+        return out
 class AcquisitionGeometry(object):
     
     def __init__(self, 
@@ -192,9 +196,12 @@ class AcquisitionGeometry(object):
         repres += "distance center-detector: {0}\n".format(self.dist_source_center)
         repres += "number of channels: {0}\n".format(self.channels)
         return repres
-
-
-            
+    def allocate(self, value=0, dimension_labels=None):
+        '''allocates an AcquisitionData according to the size expressed in the instance'''
+        out = AcquisitionData(geometry=self, dimension_labels=dimension_labels)
+        if value != 0:
+            out += value
+        return out
 class DataContainer(object):
     '''Generic class to hold data
     
@@ -741,6 +748,7 @@ class DataContainer(object):
         return numpy.sqrt(self.squared_norm())
     
     
+    
 class ImageData(DataContainer):
     '''DataContainer for holding 2D or 3D DataContainer'''
     def __init__(self, 
@@ -903,8 +911,11 @@ class AcquisitionData(DataContainer):
                         elif dim == 'horizontal':
                             shape.append(horiz)
                     if len(shape) != len(dimension_labels):
-                        raise ValueError('Missing {0} axes'.format(
-                                len(dimension_labels) - len(shape)))
+                        raise ValueError('Missing {0} axes.\nExpected{1} got {2}}'\
+                            .format(
+                                len(dimension_labels) - len(shape),
+                                dimension_labels, shape) 
+                            )
                     shape = tuple(shape)
                     
                 array = numpy.zeros( shape , dtype=numpy.float32) 
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
index 322e9eb..798fb61 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FBPD.py
@@ -1,86 +1,86 @@
-# -*- coding: utf-8 -*-

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

-#   Visual Analytics and Imaging System Group of the Science Technology

-#   Facilities Council, STFC

-

-#   Copyright 2019 Edoardo Pasca

-

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

-"""

-Created on Thu Feb 21 11:09:03 2019

-

-@author: ofn77899

-"""

-

-from ccpi.optimisation.algorithms import Algorithm

-from ccpi.optimisation.funcs import ZeroFun

-        

-class FBPD(Algorithm):

-    '''FBPD Algorithm

-    

-    Parameters:

-      x_init: initial guess

-      f: constraint

-      g: data fidelity

-      h: regularizer

-      opt: additional algorithm 

-    '''

-    constraint = None

-    data_fidelity = None

-    regulariser = None

-    def __init__(self, **kwargs):

-        pass

-    def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\

-         regulariser=None, opt=None):

-

-        # default inputs

-        if constraint    is None: 

-            self.constraint    = ZeroFun()

-        else:

-            self.constraint = constraint

-        if data_fidelity is None:

-            data_fidelity = ZeroFun()

-        else:

-            self.data_fidelity = data_fidelity

-        if regulariser   is None:

-            self.regulariser   = ZeroFun()

-        else:

-            self.regulariser = regulariser

-        

-        # algorithmic parameters

-        

-        

-        # step-sizes

-        self.tau   = 2 / (self.data_fidelity.L + 2)

-        self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L

-        

-        self.inv_sigma = 1/self.sigma

-    

-        # initialization

-        self.x = x_init

-        self.y = operator.direct(self.x)

-        

-    

-    def update(self):

-    

-        # primal forward-backward step

-        x_old = self.x

-        self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) )

-        self.x = self.constraint.prox(self.x, self.tau);

-    

-        # dual forward-backward step

-        self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old);

-        self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma);   

-

-        # time and criterion

-        self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x))

+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2019 Edoardo Pasca
+
+#   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.
+"""
+Created on Thu Feb 21 11:09:03 2019
+
+@author: ofn77899
+"""
+
+from ccpi.optimisation.algorithms import Algorithm
+from ccpi.optimisation.funcs import ZeroFun
+        
+class FBPD(Algorithm):
+    '''FBPD Algorithm
+    
+    Parameters:
+      x_init: initial guess
+      f: constraint
+      g: data fidelity
+      h: regularizer
+      opt: additional algorithm 
+    '''
+    constraint = None
+    data_fidelity = None
+    regulariser = None
+    def __init__(self, **kwargs):
+        pass
+    def set_up(self, x_init, operator=None, constraint=None, data_fidelity=None,\
+         regulariser=None, opt=None):
+
+        # default inputs
+        if constraint    is None: 
+            self.constraint    = ZeroFun()
+        else:
+            self.constraint = constraint
+        if data_fidelity is None:
+            data_fidelity = ZeroFun()
+        else:
+            self.data_fidelity = data_fidelity
+        if regulariser   is None:
+            self.regulariser   = ZeroFun()
+        else:
+            self.regulariser = regulariser
+        
+        # algorithmic parameters
+        
+        
+        # step-sizes
+        self.tau   = 2 / (self.data_fidelity.L + 2)
+        self.sigma = (1/self.tau - self.data_fidelity.L/2) / self.regulariser.L
+        
+        self.inv_sigma = 1/self.sigma
+    
+        # initialization
+        self.x = x_init
+        self.y = operator.direct(self.x)
+        
+    
+    def update(self):
+    
+        # primal forward-backward step
+        x_old = self.x
+        self.x = self.x - self.tau * ( self.data_fidelity.grad(self.x) + self.operator.adjoint(self.y) )
+        self.x = self.constraint.prox(self.x, self.tau);
+    
+        # dual forward-backward step
+        self.y = self.y + self.sigma * self.operator.direct(2*self.x - x_old);
+        self.y = self.y - self.sigma * self.regulariser.prox(self.inv_sigma*self.y, self.inv_sigma);   
+
+        # time and criterion
+        self.loss = self.constraint(self.x) + self.data_fidelity(self.x) + self.regulariser(self.operator.direct(self.x))
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
index 52fe6d7..903bc30 100755..100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
@@ -1,29 +1,29 @@
-# -*- coding: utf-8 -*-

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

-#   Visual Analytics and Imaging System Group of the Science Technology

-#   Facilities Council, STFC

-

-#   Copyright 2019 Edoardo Pasca

-

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

-"""

-Created on Thu Feb 21 11:03:13 2019

-

-@author: ofn77899

-"""

-

-from .Algorithm import Algorithm

-from .CGLS import CGLS

-from .GradientDescent import GradientDescent

-from .FISTA import FISTA

-from .FBPD import FBPD
\ No newline at end of file
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2019 Edoardo Pasca
+
+#   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.
+"""
+Created on Thu Feb 21 11:03:13 2019
+
+@author: ofn77899
+"""
+
+from .Algorithm import Algorithm
+from .CGLS import CGLS
+from .GradientDescent import GradientDescent
+from .FISTA import FISTA
+from .FBPD import FBPD
diff --git a/Wrappers/Python/test/test_DataContainer.py b/Wrappers/Python/test/test_DataContainer.py
index 05f3fe8..3def054 100755
--- a/Wrappers/Python/test/test_DataContainer.py
+++ b/Wrappers/Python/test/test_DataContainer.py
@@ -457,6 +457,49 @@ class TestDataContainer(unittest.TestCase):
                                         pixel_num_h=5, channels=2)
         sino = AcquisitionData(geometry=sgeometry)
         self.assertEqual(sino.shape, (2, 10, 3, 5))
+    def test_ImageGeometry_allocate(self):
+        vgeometry = ImageGeometry(voxel_num_x=4, voxel_num_y=3, channels=2)
+        image = vgeometry.allocate()
+        self.assertEqual(0,image.as_array()[0][0][0])
+        image = vgeometry.allocate(1)
+        self.assertEqual(1,image.as_array()[0][0][0])
+        default_order = ['channel' , 'horizontal_y' , 'horizontal_x']
+        self.assertEqual(default_order[0], image.dimension_labels[0])
+        self.assertEqual(default_order[1], image.dimension_labels[1])
+        self.assertEqual(default_order[2], image.dimension_labels[2])
+        order = [ 'horizontal_x' , 'horizontal_y', 'channel' ]
+        image = vgeometry.allocate(0,dimension_labels=order)
+        self.assertEqual(order[0], image.dimension_labels[0])
+        self.assertEqual(order[1], image.dimension_labels[1])
+        self.assertEqual(order[2], image.dimension_labels[2])
+    def test_AcquisitionGeometry_allocate(self):
+        ageometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=10),
+                                        geom_type='parallel', pixel_num_v=3,
+                                        pixel_num_h=5, channels=2)
+        sino = ageometry.allocate()
+        shape = sino.shape
+        print ("shape", shape)
+        self.assertEqual(0,sino.as_array()[0][0][0][0])
+        self.assertEqual(0,sino.as_array()[shape[0]-1][shape[1]-1][shape[2]-1][shape[3]-1])
+        
+        sino = ageometry.allocate(1)
+        self.assertEqual(1,sino.as_array()[0][0][0][0])
+        self.assertEqual(1,sino.as_array()[shape[0]-1][shape[1]-1][shape[2]-1][shape[3]-1])
+        print (sino.dimension_labels, sino.shape, ageometry)
+        
+        default_order = ['channel' , ' angle' ,
+                                          'vertical' , 'horizontal']
+        self.assertEqual(default_order[0], sino.dimension_labels[0])
+        self.assertEqual(default_order[1], sino.dimension_labels[1])
+        self.assertEqual(default_order[2], sino.dimension_labels[2])
+        self.assertEqual(default_order[3], sino.dimension_labels[3])
+        order = ['vertical' , 'horizontal', 'channel' , 'angle' ]
+        sino = ageometry.allocate(0,dimension_labels=order)
+        print (sino.dimension_labels, sino.shape, ageometry)
+        self.assertEqual(order[0], sino.dimension_labels[0])
+        self.assertEqual(order[1], sino.dimension_labels[1])
+        self.assertEqual(order[2], sino.dimension_labels[2])
+        self.assertEqual(order[2], sino.dimension_labels[2])
 
     def assertNumpyArrayEqual(self, first, second):
         res = True