lots of changes and added tests

12 files changed, 257 insertions, 66 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index 1bfc98c..cbce844 100755
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -26,7 +26,7 @@ class BlockDataContainer(object):
         if shape is None:

             shape = (len(args),1)

         self.shape = shape

-        

+        print (self.shape)

         n_elements = functools.reduce(lambda x,y: x*y, shape, 1)

         if len(args) != n_elements:

             raise ValueError(

@@ -78,74 +78,79 @@ class BlockDataContainer(object):
         assert self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.add(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.add(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list) or isinstance(other, numpy.ndarray):

-            return type(self)(*[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(

+            *[ el.add(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],

+            shape=self.shape)

         

     def subtract(self, other, *args, **kwargs):

         assert self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.subtract(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.subtract(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list) or isinstance(other, numpy.ndarray):

-            return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(*[ el.subtract(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],

+                          shape=self.shape)

 

     def multiply(self, other, *args, **kwargs):

         self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.multiply(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.multiply(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list):

-            return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

+            return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

         elif isinstance(other, numpy.ndarray):           

-            return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(*[ el.multiply(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],

+                          shape=self.shape)

     

     def divide(self, other, *args, **kwargs):

         self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.divide(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.divide(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list) or isinstance(other, numpy.ndarray):

-            return type(self)(*[ el.divide(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.divide(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.divide(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(*[ el.divide(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],

+                          shape=self.shape)

     

     def power(self, other, *args, **kwargs):

         assert self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.power(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.power(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list) or isinstance(other, numpy.ndarray):

-            return type(self)(*[ el.power(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.power(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.power(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(*[ el.power(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)

     

     def maximum(self,other, *args, **kwargs):

         assert self.is_compatible(other)

         out = kwargs.get('out', None)

         if isinstance(other, Number):

-            return type(self)(*[ el.maximum(other, out, *args, **kwargs) for el in self.containers])

+            return type(self)(*[ el.maximum(other, out, *args, **kwargs) for el in self.containers], shape=self.shape)

         elif isinstance(other, list) or isinstance(other, numpy.ndarray):

-            return type(self)(*[ el.maximum(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)])

-        return type(self)(*[ el.maximum(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)])

+            return type(self)(*[ el.maximum(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

+        return type(self)(*[ el.maximum(ot, out, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)

     

     ## unary operations    

     def abs(self, *args,  **kwargs):

         out = kwargs.get('out', None)

-        return type(self)(*[ el.abs(out, *args, **kwargs) for el in self.containers]) 

+        return type(self)(*[ el.abs(out, *args, **kwargs) for el in self.containers], shape=self.shape)

     def sign(self, *args,  **kwargs):

         out = kwargs.get('out', None)

-        return type(self)(*[ el.sign(out, *args, **kwargs) for el in self.containers])

+        return type(self)(*[ el.sign(out, *args, **kwargs) for el in self.containers], shape=self.shape)

     def sqrt(self, *args,  **kwargs):

         out = kwargs.get('out', None)

-        return type(self)(*[ el.sqrt(out, *args, **kwargs) for el in self.containers])

+        return type(self)(*[ el.sqrt(out, *args, **kwargs) for el in self.containers], shape=self.shape)

     def conjugate(self, out=None):

-        return type(self)(*[el.conjugate() for el in self.containers])

+        return type(self)(*[el.conjugate() for el in self.containers], shape=self.shape)

     

     ## reductions

     def sum(self, *args, **kwargs):

-        return numpy.asarray([ el.sum(*args, **kwargs) for el in self.containers])

+        return numpy.sum([ el.sum(*args, **kwargs) for el in self.containers])

     def squared_norm(self):

         y = numpy.asarray([el.squared_norm() for el in self.containers])

         return y.sum() 

@@ -155,7 +160,7 @@ class BlockDataContainer(object):
         '''alias of clone'''    

         return self.clone()

     def clone(self):

-        return type(self)(*[el.copy() for el in self.containers])

+        return type(self)(*[el.copy() for el in self.containers], shape=self.shape)

     

     def __add__(self, other):

         return self.add( other )

@@ -297,4 +302,8 @@ class BlockDataContainer(object):
     def __itruediv__(self, other):

         '''Inline truedivision'''

         return self.__idiv__(other)

-    

+    @property

+    def T(self):

+       '''return the transposed of self'''

+       shape = (self.shape[1], self.shape[0])

+       return type(self)(*self.containers, shape=shape)

diff --git a/Wrappers/Python/ccpi/framework/__init__.py b/Wrappers/Python/ccpi/framework/__init__.py
index 083f547..4683c21 100755
--- a/Wrappers/Python/ccpi/framework/__init__.py
+++ b/Wrappers/Python/ccpi/framework/__init__.py
@@ -20,5 +20,5 @@ from .framework import ImageData, AcquisitionData
 from .framework import ImageGeometry, AcquisitionGeometry

 from .framework import find_key, message

 from .framework import DataProcessor

-from .framework import AX, PixelByPixelDataProcessor

+from .framework import AX, PixelByPixelDataProcessor, CastDataContainer

 from .BlockDataContainer import BlockDataContainer

diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index 3159cc7..1413e21 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -734,7 +734,9 @@ class DataContainer(object):
     
     def sqrt(self, *args,  **kwargs):
         return self.pixel_wise_unary(numpy.sqrt, *args,  **kwargs)
-    
+
+    def conjugate(self, *args,  **kwargs):
+        return self.pixel_wise_unary(numpy.conjugate, *args,  **kwargs)
     #def __abs__(self):
     #    operation = FM.OPERATION.ABS
     #    return self.callFieldMath(operation, None, self.mask, self.maskOnValue)
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
index 903bc30..7e500e8 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
@@ -27,3 +27,4 @@ from .CGLS import CGLS
 from .GradientDescent import GradientDescent
 from .FISTA import FISTA
 from .FBPD import FBPD
+
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index b2af8fc..a83dc8a 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -66,9 +66,11 @@ class BlockOperator(Operator):
     
     def direct(self, x, out=None):
         shape = self.get_output_shape(x.shape)
+        print ("direct output shape", shape)
         res = []
         for row in range(self.shape[0]):
             for col in range(self.shape[1]):
+                print ("row {} col {}".format(row, col))
                 if col == 0:
                     prod = self.get_item(row,col).direct(x.get_item(col))
                 else:
@@ -76,19 +78,8 @@ class BlockOperator(Operator):
             res.append(prod)
         return BlockDataContainer(*res, shape=shape)
     
-    def adjoint(self, x, out=None):
-        shape = self.get_output_shape(x.shape, adjoint=True)
-        res = []
-        for row in range(self.shape[1]):
-            for col in range(self.shape[0]):
-                if col == 0:
-                    prod = self.get_item(col,row).adjoint(x.get_item(row))
-                else:
-                    prod += self.get_item(col,row).adjoint(x.get_item(row))
-            res.append(prod)
-        return BlockDataContainer(*res, shape=shape)
-    
     def get_output_shape(self, xshape, adjoint=False):
+        print ("get_output_shape", self.shape, xshape)
         sshape = self.shape[1]
         oshape = self.shape[0]
         if adjoint:
@@ -98,8 +89,27 @@ class BlockOperator(Operator):
             raise ValueError('Incompatible shapes {} {}'.format(self.shape, xshape))
         return (oshape, xshape[-1])
     
-    
-    
+    def __rmul__(self, scalar):
+        '''Defines the left multiplication with a scalar
+
+        Args: scalar (number or iterable containing numbers):
+
+        Returns: a block operator with Scaled Operators inside'''
+        if isinstance (scalar, list) or isinstance(scalar, tuple) or \
+                isinstance(scalar, numpy.ndarray):
+            if len(scalar) != len(self.operators):
+                raise ValueError('dimensions of scalars and operators do not match')
+            scalars = scalar
+        else:
+            scalars = [scalar for _ in self.operators]
+        # create a list of ScaledOperator-s
+        ops = [ v * op for v,op in zip(scalars, self.operators)]
+        return BlockOperator(*ops, shape=self.shape)
+    def T(self):
+        '''Return the transposed of self'''
+        shape = (self.shape[1], self.shape[0])
+        return type(self)(*self.operators, shape=shape)
+
 class BlockLinearOperator(BlockOperator):
     '''Class to hold a block operator
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
index 29dacb8..a47bec2 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py
@@ -4,12 +4,8 @@ Created on Thu Feb 14 12:36:40 2019
 
 @author: ofn77899
 """
-#from ccpi.optimisation.ops import Operator
 import numpy
-from numbers import Number
-import functools
-from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer
-from ccpi.optimisation.operators import Operator, LinearOperator
+from ccpi.optimisation.operators import BlockOperator
 
 
        
diff --git a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
index d0e7804..e19304f 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/LinearOperator.py
@@ -7,8 +7,11 @@ Created on Tue Mar  5 15:57:52 2019
 

 from ccpi.optimisation.operators import Operator

 

+

 class LinearOperator(Operator):

-    '''Operator that maps from a space X -> Y'''

+    '''A Linear Operator that maps from a space X <-> Y'''

+    def __init__(self):

+        super(LinearOperator, self).__init__()

     def is_linear(self):

         '''Returns if the operator is linear'''

         return True

diff --git a/Wrappers/Python/ccpi/optimisation/operators/Operator.py b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
index 11b9b87..95082f4 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/Operator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/Operator.py
@@ -4,20 +4,15 @@ Created on Tue Mar  5 15:55:56 2019
 

 @author: ofn77899

 """

-from ccpi.framework,operators import ScaledOperator

+from ccpi.optimisation.operators import ScaledOperator

 

 class Operator(object):

     '''Operator that maps from a space X -> Y'''

-    def __init__(self, **kwargs):

-        self.scalar = 1

     def is_linear(self):

         '''Returns if the operator is linear'''

         return False

     def direct(self,x, out=None):

         raise NotImplementedError

-    def size(self):

-        # To be defined for specific class

-        raise NotImplementedError

     def norm(self):

         raise NotImplementedError

     def range_geometry(self):

@@ -26,4 +21,3 @@ class Operator(object):
         raise NotImplementedError

     def __rmul__(self, scalar):

         return ScaledOperator(self, scalar)

-    

diff --git a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
index 1181604..adcc6d9 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/ScaledOperator.py
@@ -1,10 +1,10 @@
-from ccpi.optimisation.operators import LinearOperator
 from numbers import Number
+import numpy
 
-class ScaledOperator(LinearOperator):
+class ScaledOperator(object):
     '''ScaledOperator
 
-    A class to represent the scalar multiplication of an Operator with a scalar. 
+    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.
@@ -22,8 +22,9 @@ class ScaledOperator(LinearOperator):
        sop.domain_geometry() = operator.domain_geometry()
     '''
     def __init__(self, operator, scalar):
+        super(ScaledOperator, self).__init__()
         if not isinstance (scalar, Number):
-            raise TypeError('expected scalar: got {}'.format(type(scalar))
+            raise TypeError('expected scalar: got {}'.format(type(scalar)))
         self.scalar = scalar
         self.operator = operator
     def direct(self, x, out=None):
@@ -31,12 +32,11 @@ class ScaledOperator(LinearOperator):
     def adjoint(self, x, out=None):
         if self.operator.is_linear():
             return self.scalar * self.operator.adjoint(x, out=out)
-    def size(self):
-        return self.operator.size()
+        else:
+            raise TypeError('Operator is not linear')
     def norm(self):
-        return self.operator.norm()
+        return numpy.abs(self.scalar) * self.operator.norm()
     def range_geometry(self):
         return self.operator.range_geometry()
     def domain_geometry(self):
         return self.operator.domain_geometry()
-        
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_BlockDataContainer.py b/Wrappers/Python/test/test_BlockDataContainer.py
index 824abf6..c14a101 100755
--- a/Wrappers/Python/test/test_BlockDataContainer.py
+++ b/Wrappers/Python/test/test_BlockDataContainer.py
@@ -16,8 +16,80 @@ from ccpi.framework import ImageData, AcquisitionData
 #from ccpi.optimisation.algorithms import GradientDescent

 from ccpi.framework import BlockDataContainer

 #from ccpi.optimisation.Algorithms import CGLS

+import functools

 

 class TestBlockDataContainer(unittest.TestCase):

+    def test_BlockDataContainerShape(self):

+        print ("test block data container")

+        ig0 = ImageGeometry(2,3,4)

+        ig1 = ImageGeometry(12,42,55,32)

+

+        data0 = ImageData(geometry=ig0)

+        data1 = ImageData(geometry=ig1) + 1

+

+        data2 = ImageData(geometry=ig0) + 2

+        data3 = ImageData(geometry=ig1) + 3

+

+        cp0 = BlockDataContainer(data0,data1)

+        cp1 = BlockDataContainer(data2,data3)

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp0.T.shape == transpose_shape)

+    def test_BlockDataContainerShapeArithmetic(self):

+        print ("test block data container")

+        ig0 = ImageGeometry(2,3,4)

+        ig1 = ImageGeometry(12,42,55,32)

+

+        data0 = ImageData(geometry=ig0)

+        data1 = ImageData(geometry=ig1) + 1

+

+        data2 = ImageData(geometry=ig0) + 2

+        data3 = ImageData(geometry=ig1) + 3

+

+        cp0 = BlockDataContainer(data0,data1)

+        #cp1 = BlockDataContainer(data2,data3)

+        cp1 = cp0 + 1

+        self.assertTrue(cp1.shape == cp0.shape)

+        cp1 = cp0.T + 1

+

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T - 1

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = (cp0.T + 1)*2

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = (cp0.T + 1)/2

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.power(2.2)

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.maximum(3)

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.abs()

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.sign()

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.sqrt()

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

+        cp1 = cp0.T.conjugate()

+        transpose_shape = (cp0.shape[1], cp0.shape[0])

+        self.assertTrue(cp1.shape == transpose_shape)

+

     def test_BlockDataContainer(self):

         print ("test block data container")

         ig0 = ImageGeometry(2,3,4)

@@ -198,7 +270,9 @@ class TestBlockDataContainer(unittest.TestCase):
         numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], numpy.sqrt(4), decimal=4)

         

         s = cp0.sum()

-        numpy.testing.assert_almost_equal(s[0], 0, decimal=4)

+        size = functools.reduce(lambda x,y: x*y, data1.shape, 1)

+        print ("size" , size)

+        numpy.testing.assert_almost_equal(s, 0 + size, decimal=4)

         s0 = 1

         s1 = 1

         for i in cp0.get_item(0,0).shape:

@@ -206,5 +280,5 @@ class TestBlockDataContainer(unittest.TestCase):
         for i in cp0.get_item(1,0).shape:

             s1 *= i

             

-        numpy.testing.assert_almost_equal(s[1], cp0.get_item(0,0).as_array()[0][0][0]*s0 +cp0.get_item(1,0).as_array()[0][0][0]*s1, decimal=4)

+        #numpy.testing.assert_almost_equal(s[1], cp0.get_item(0,0).as_array()[0][0][0]*s0 +cp0.get_item(1,0).as_array()[0][0][0]*s1, decimal=4)

     
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py
new file mode 100644
index 0000000..f3b057b
--- /dev/null
+++ b/Wrappers/Python/test/test_BlockOperator.py
@@ -0,0 +1,78 @@
+import unittest
+from ccpi.optimisation.operators import BlockOperator
+from ccpi.framework import BlockDataContainer
+from ccpi.optimisation.ops import TomoIdentity
+from ccpi.framework import ImageGeometry, ImageData
+import numpy
+
+class TestBlockOperator(unittest.TestCase):
+
+    def test_BlockOperator(self):
+        ig = [ ImageGeometry(10,20,30) , \
+               ImageGeometry(11,21,31) , \
+               ImageGeometry(12,22,32) ]
+        x = [ g.allocate() for g in ig ]
+        ops = [ TomoIdentity(g) for g in ig ]
+
+        K = BlockOperator(*ops)
+        #X = BlockDataContainer(*x).T + 1
+        X = BlockDataContainer(x[0])
+        Y = K.direct(X)
+        #self.assertTrue(Y.shape == X.shape)
+
+        numpy.testing.assert_array_equal(Y.get_item(0).as_array(),X.get_item(0).as_array())
+        numpy.testing.assert_array_equal(Y.get_item(1).as_array(),X.get_item(0).as_array())
+        #numpy.testing.assert_array_equal(Y.get_item(2).as_array(),X.get_item(2).as_array())
+
+
+    def test_ScaledBlockOperatorSingleScalar(self):
+        ig = [ ImageGeometry(10,20,30) , \
+               ImageGeometry(11,21,31) , \
+               ImageGeometry(12,22,32) ]
+        x = [ g.allocate() for g in ig ]
+        ops = [ TomoIdentity(g) for g in ig ]
+
+        scalar = 0.5
+        K = 0.5 * BlockOperator(*ops)
+        X = BlockDataContainer(*x) + 1
+        print (X.shape)
+        X = BlockDataContainer(*x).T + 1
+        print (X.shape, K.shape)
+        Y = K.direct(X)
+        self.assertTrue(Y.shape == X.shape)
+
+        numpy.testing.assert_array_equal(Y.get_item(0).as_array(),scalar * X.get_item(0).as_array())
+        numpy.testing.assert_array_equal(Y.get_item(1).as_array(),scalar * X.get_item(1).as_array())
+        numpy.testing.assert_array_equal(Y.get_item(2).as_array(),scalar * X.get_item(2).as_array())
+
+    def test_ScaledBlockOperatorScalarList(self):
+        ig = [ImageGeometry(10, 20, 30),
+              ImageGeometry(11, 21, 31),
+              ImageGeometry(12, 22, 32)]
+        x = [g.allocate() for g in ig]
+        ops = [TomoIdentity(g) for g in ig]
+
+        scalar = [i*1.2 for i, el in enumerate(ig)]
+
+        K = scalar * BlockOperator(*ops)
+        X = BlockDataContainer(*x).T + 1
+        Y = K.direct(X)
+        self.assertTrue(Y.shape == X.shape)
+
+        numpy.testing.assert_array_equal(Y.get_item(0).as_array(),
+                                         scalar[0] * X.get_item(0).as_array())
+        numpy.testing.assert_array_equal(Y.get_item(1).as_array(),
+                                         scalar[1] * X.get_item(1).as_array())
+        numpy.testing.assert_array_equal(Y.get_item(2).as_array(),
+                                         scalar[2] * X.get_item(2).as_array())
+
+
+    def test_TomoIdentity(self):
+        ig = ImageGeometry(10,20,30)
+        img = ig.allocate()
+        self.assertTrue(img.shape == (30,20,10))
+        self.assertEqual(img.sum(), 0)
+        Id = TomoIdentity(ig)
+        y = Id.direct(img)
+        numpy.testing.assert_array_equal(y.as_array(), img.as_array())
+
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
new file mode 100644
index 0000000..46e8c7c
--- /dev/null
+++ b/Wrappers/Python/test/test_Operator.py
@@ -0,0 +1,24 @@
+import unittest
+#from ccpi.optimisation.operators import Operator
+from ccpi.optimisation.ops import TomoIdentity
+from ccpi.framework import ImageGeometry, ImageData
+import numpy
+
+class TestOperator(unittest.TestCase):
+    def test_ScaledOperator(self):
+        ig = ImageGeometry(10,20,30)
+        img = ig.allocate()
+        scalar = 0.5
+        sid = scalar * TomoIdentity(ig)
+        numpy.testing.assert_array_equal(scalar * img.as_array(), sid.direct(img).as_array())
+        
+
+    def test_TomoIdentity(self):
+        ig = ImageGeometry(10,20,30)
+        img = ig.allocate()
+        self.assertTrue(img.shape == (30,20,10))
+        self.assertEqual(img.sum(), 0)
+        Id = TomoIdentity(ig)
+        y = Id.direct(img)
+        numpy.testing.assert_array_equal(y.as_array(), img.as_array())
+