fix algebra BDC

5 files changed, 580 insertions, 391 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
index f1d6d9a..529a1ce 100755
--- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py
+++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py
@@ -1,356 +1,470 @@
-    # -*- coding: utf-8 -*-

-"""

-Created on Tue Mar  5 16:04:45 2019

-

-@author: ofn77899

-"""

-from __future__ import absolute_import

-from __future__ import division

-from __future__ import print_function

-from __future__ import unicode_literals

-

-import numpy

-from numbers import Number

-import functools

-from ccpi.framework import DataContainer

-#from ccpi.framework import AcquisitionData, ImageData

-#from ccpi.optimisation.operators import Operator, LinearOperator

- 

-class BlockDataContainer(object):

-    '''Class to hold DataContainers as column vector'''

-    __array_priority__ = 1

-    def __init__(self, *args, **kwargs):

-        ''''''

-        self.containers = args

-        self.index = 0

-        shape = kwargs.get('shape', None)

-        if shape is None:

-           shape = (len(args),1)

-#        shape = (len(args),1)

-        self.shape = shape

-

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

-        if len(args) != n_elements:

-            raise ValueError(

-                    'Dimension and size do not match: expected {} got {}'

-                    .format(n_elements, len(args)))

-

-        

-    def __iter__(self):

-        '''BlockDataContainer is Iterable'''

-        return self

-    def next(self):

-        '''python2 backwards compatibility'''

-        return self.__next__()

-    def __next__(self):

-        try:

-            out = self[self.index]

-        except IndexError as ie:

-            raise StopIteration()

-        self.index+=1

-        return out

-    

-    def is_compatible(self, other):

-        '''basic check if the size of the 2 objects fit'''

-        

-#        for i in range(len(self.containers)):

-#            if type(self.containers[i])==type(self):

-#                self = self.containers[i]

-        

-        if isinstance(other, Number):

-            return True   

-        elif isinstance(other, list):

-            for ot in other:

-                if not isinstance(ot, (Number,\

-                                 numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\

-                                 numpy.float, numpy.float16, numpy.float32, numpy.float64, \

-                                 numpy.complex)):

-                    raise ValueError('List/ numpy array can only contain numbers {}'\

-                                     .format(type(ot)))

-            return len(self.containers) == len(other)

-        elif isinstance(other, numpy.ndarray):

-            return len(self.containers) == len(other)

-        

-        elif issubclass(other.__class__, DataContainer):

-            ret = True

-            for i, el in enumerate(self.containers):

-                if isinstance(el, BlockDataContainer):

-                    a = el.is_compatible(other)

-                else:

-                    a = el.shape == other.shape

-#                print ("current element" , el.shape, "other ", other.shape, "same shape" , a)

-                ret = ret and a

-            return ret

-        

-        

-#        elif issubclass(other.__class__, DataContainer):

-#            return self.get_item(0).shape == other.shape

-        return len(self.containers) == len(other.containers)

-

-    def get_item(self, row):

-        if row > self.shape[0]:

-            raise ValueError('Requested row {} > max {}'.format(row, self.shape[0]))

-        return self.containers[row]

-

-    def __getitem__(self, row):

-        return self.get_item(row)

-                

-    def add(self, other, *args, **kwargs):

-        if not self.is_compatible(other):

-            raise ValueError('Incompatible for add')

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

-        #print ("args" , *args)

-        if isinstance(other, Number):

-            return type(self)(*[ el.add(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

-        elif issubclass(other.__class__, DataContainer):

-            # try to do algebra with one DataContainer. Will raise error if not compatible

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

-        

-        return type(self)(

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

-            shape=self.shape)

-        

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

-        if not self.is_compatible(other):

-            raise ValueError('Incompatible for subtract')

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

-        if isinstance(other, Number):

-            return type(self)(*[ el.subtract(other,  *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,  *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

-        elif issubclass(other.__class__, DataContainer):

-            # try to do algebra with one DataContainer. Will raise error if not compatible

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

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

-                          shape=self.shape)

-

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

-        if not self.is_compatible(other):

-            raise ValueError('{} Incompatible for multiply'.format(other))

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

-        if isinstance(other, Number):

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

-        elif isinstance(other, list):

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

-        elif isinstance(other, numpy.ndarray):           

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

-        elif issubclass(other.__class__, DataContainer):

-            # try to do algebra with one DataContainer. Will raise error if not compatible

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

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

-                          shape=self.shape)

-    

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

-        if not self.is_compatible(other):

-            raise ValueError('Incompatible for divide')

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

-        if isinstance(other, Number):

-            return type(self)(*[ el.divide(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

-        elif issubclass(other.__class__, DataContainer):

-            # try to do algebra with one DataContainer. Will raise error if not compatible

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

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

-                          shape=self.shape)

-    

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

-        if not self.is_compatible(other):

-            raise ValueError('Incompatible for power')

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

-        if isinstance(other, Number):

-            return type(self)(*[ el.power(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

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

-    

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

-        if not self.is_compatible(other):

-            raise ValueError('Incompatible for maximum')

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

-        if isinstance(other, Number):

-            return type(self)(*[ el.maximum(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)

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

-    

-    ## unary operations    

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

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

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

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

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

-        return type(self)(*[ el.sqrt(*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], shape=self.shape)

-    

-    ## reductions

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

-        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() 

-    def norm(self):

-        return numpy.sqrt(self.squared_norm())    

-    def copy(self):

-        '''alias of clone'''    

-        return self.clone()

-    def clone(self):

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

-

-    def fill(self, other):

-        if isinstance (other, BlockDataContainer):

-            if not self.is_compatible(other):

-                raise ValueError('Incompatible containers')

-            for el,ot in zip(self.containers, other.containers):

-                el.fill(ot)

-        else:

-            return ValueError('Cannot fill with object provided {}'.format(type(other)))

-    

-    def __add__(self, other):

-        return self.add( other )

-    # __radd__

-    

-    def __sub__(self, other):

-        return self.subtract( other )

-    # __rsub__

-    

-    def __mul__(self, other):

-        return self.multiply(other)

-    # __rmul__

-    

-    def __div__(self, other):

-        return self.divide(other)

-    # __rdiv__

-    def __truediv__(self, other):

-        return self.divide(other)

-    

-    def __pow__(self, other):

-        return self.power(other)

-    # reverse operand

-    def __radd__(self, other):

-        '''Reverse addition

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return self + other

-    # __radd__

-    

-    def __rsub__(self, other):

-        '''Reverse subtraction

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return (-1 * self) + other

-    # __rsub__

-    

-    def __rmul__(self, other):

-        '''Reverse multiplication

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return self * other

-    # __rmul__

-    

-    def __rdiv__(self, other):

-        '''Reverse division

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return pow(self / other, -1)

-    # __rdiv__

-    def __rtruediv__(self, other):

-        '''Reverse truedivision

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return self.__rdiv__(other)

-    

-    def __rpow__(self, other):

-        '''Reverse power

-        

-        to make sure that this method is called rather than the __mul__ of a numpy array

-        the class constant __array_priority__ must be set > 0

-        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__

-        '''

-        return other.power(self)

-    

-    def __iadd__(self, other):

-        '''Inline addition'''

-        if isinstance (other, BlockDataContainer):

-            for el,ot in zip(self.containers, other.containers):

-                el += ot

-        elif isinstance(other, Number):

-            for el in self.containers:

-                el += other

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

-            if not self.is_compatible(other):

-                raise ValueError('Incompatible for __iadd__')

-            for el,ot in zip(self.containers, other):

-                el += ot

-        return self

-    # __iadd__

-    

-    def __isub__(self, other):

-        '''Inline subtraction'''

-        if isinstance (other, BlockDataContainer):

-            for el,ot in zip(self.containers, other.containers):

-                el -= ot

-        elif isinstance(other, Number):

-            for el in self.containers:

-                el -= other

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

-            if not self.is_compatible(other):

-                raise ValueError('Incompatible for __isub__')

-            for el,ot in zip(self.containers, other):

-                el -= ot

-        return self

-    # __isub__

-    

-    def __imul__(self, other):

-        '''Inline multiplication'''

-        if isinstance (other, BlockDataContainer):

-            for el,ot in zip(self.containers, other.containers):

-                el *= ot

-        elif isinstance(other, Number):

-            for el in self.containers:

-                el *= other

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

-            if not self.is_compatible(other):

-                raise ValueError('Incompatible for __imul__')

-            for el,ot in zip(self.containers, other):

-                el *= ot

-        return self

-    # __imul__

-    

-    def __idiv__(self, other):

-        '''Inline division'''

-        if isinstance (other, BlockDataContainer):

-            for el,ot in zip(self.containers, other.containers):

-                el /= ot

-        elif isinstance(other, Number):

-            for el in self.containers:

-                el /= other

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

-            if not self.is_compatible(other):

-                raise ValueError('Incompatible for __idiv__')

-            for el,ot in zip(self.containers, other):

-                el /= ot

-        return self

-    # __rdiv__

-    def __itruediv__(self, other):

-        '''Inline truedivision'''

-        return self.__idiv__(other)

-

+    # -*- coding: utf-8 -*-
+"""
+Created on Tue Mar  5 16:04:45 2019
+
+@author: ofn77899
+"""
+from __future__ import absolute_import
+from __future__ import division
+from __future__ import print_function
+from __future__ import unicode_literals
+
+import numpy
+from numbers import Number
+import functools
+from ccpi.framework import DataContainer
+#from ccpi.framework import AcquisitionData, ImageData
+#from ccpi.optimisation.operators import Operator, LinearOperator
+ 
+class BlockDataContainer(object):
+    '''Class to hold DataContainers as column vector
+    
+    Provides basic algebra between BlockDataContainer's, DataContainer's and
+    subclasses and Numbers
+
+    1) algebra between `BlockDataContainer`s will be element-wise, only if 
+       the shape of the 2 `BlockDataContainer`s is the same, otherwise it 
+       will fail 
+    2) algebra between `BlockDataContainer`s and `list` or `numpy array` will 
+       work as long as the number of `rows` and element of the arrays match,
+       indipendently on the fact that the `BlockDataContainer` could be nested
+    3) algebra between `BlockDataContainer` and one `DataContainer` is possible.
+       It will require that all the `DataContainers` in the block to be
+       compatible with the `DataContainer` we want to algebra with. Should we 
+       require that the `DataContainer` is the same type? Like `ImageData` or `AcquisitionData`?
+    4) algebra between `BlockDataContainer` and a `Number` is possible and it
+       will be done with each element of the `BlockDataContainer` even if nested
+
+    A = [ [B,C] , D] 
+    A * 3 = [ 3 * [B,C] , 3* D] = [ [ 3*B, 3*C]  , 3*D ]
+    
+    '''
+    ADD = 'add'
+    SUBTRACT = 'subtract'
+    MULTIPLY = 'multiply'
+    DIVIDE = 'divide'
+    POWER = 'power'
+    __array_priority__ = 1
+    def __init__(self, *args, **kwargs):
+        ''''''
+        self.containers = args
+        self.index = 0
+        shape = kwargs.get('shape', None)
+        if shape is None:
+           shape = (len(args),1)
+#        shape = (len(args),1)
+        self.shape = shape
+
+        n_elements = functools.reduce(lambda x,y: x*y, shape, 1)
+        if len(args) != n_elements:
+            raise ValueError(
+                    'Dimension and size do not match: expected {} got {}'
+                    .format(n_elements, len(args)))
+
+        
+    def __iter__(self):
+        '''BlockDataContainer is Iterable'''
+        return self
+    def next(self):
+        '''python2 backwards compatibility'''
+        return self.__next__()
+    def __next__(self):
+        try:
+            out = self[self.index]
+        except IndexError as ie:
+            raise StopIteration()
+        self.index+=1
+        return out
+    
+    def is_compatible(self, other):
+        '''basic check if the size of the 2 objects fit'''
+
+        if isinstance(other, Number):
+            return True   
+        elif isinstance(other, (list, numpy.ndarray)) :
+            for ot in other:
+                if not isinstance(ot, (Number,\
+                                 numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\
+                                 numpy.float, numpy.float16, numpy.float32, numpy.float64, \
+                                 numpy.complex)):
+                    raise ValueError('List/ numpy array can only contain numbers {}'\
+                                     .format(type(ot)))
+            return len(self.containers) == len(other)
+        elif issubclass(other.__class__, DataContainer):
+            ret = True
+            for i, el in enumerate(self.containers):
+                if isinstance(el, BlockDataContainer):
+                    a = el.is_compatible(other)
+                else:
+                    a = el.shape == other.shape
+                print ("current element" , el.shape, "other ", other.shape, "same shape" , a)
+                ret = ret and a
+            return ret
+            #return self.get_item(0).shape == other.shape
+        return len(self.containers) == len(other.containers)
+
+    def get_item(self, row):
+        if row > self.shape[0]:
+            raise ValueError('Requested row {} > max {}'.format(row, self.shape[0]))
+        return self.containers[row]
+
+    def __getitem__(self, row):
+        return self.get_item(row)
+                
+#    def add(self, other, *args, **kwargs):
+#        if not self.is_compatible(other):
+#            raise ValueError('Incompatible for add')
+#        out = kwargs.get('out', None)
+#        #print ("args" , *args)
+#        if isinstance(other, Number):
+#            return type(self)(*[ el.add(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+#        elif issubclass(other.__class__, DataContainer):
+#            # try to do algebra with one DataContainer. Will raise error if not compatible
+#            return type(self)(*[ el.add(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+#        
+#        return type(self)(
+#            *[ el.add(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+#            shape=self.shape)
+#        
+#    def subtract(self, other, *args, **kwargs):
+#        if not self.is_compatible(other):
+#            raise ValueError('Incompatible for subtract')
+#        out = kwargs.get('out', None)
+#        if isinstance(other, Number):
+#            return type(self)(*[ el.subtract(other,  *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,  *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+#        elif issubclass(other.__class__, DataContainer):
+#            # try to do algebra with one DataContainer. Will raise error if not compatible
+#            return type(self)(*[ el.subtract(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+#        return type(self)(*[ el.subtract(ot,  *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+#                          shape=self.shape)
+#
+#    def multiply(self, other, *args, **kwargs):
+#        if not self.is_compatible(other):
+#            raise ValueError('{} Incompatible for multiply'.format(other))
+#        out = kwargs.get('out', None)
+#        if isinstance(other, Number):
+#            return type(self)(*[ el.multiply(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+#        elif isinstance(other, list):
+#            return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+#        elif isinstance(other, numpy.ndarray):           
+#            return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+#        elif issubclass(other.__class__, DataContainer):
+#            # try to do algebra with one DataContainer. Will raise error if not compatible
+#            return type(self)(*[ el.multiply(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+#        return type(self)(*[ el.multiply(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+#                          shape=self.shape)
+#    
+#    def divide_old(self, other, *args, **kwargs):
+#        if not self.is_compatible(other):
+#            raise ValueError('Incompatible for divide')
+#        out = kwargs.get('out', None)
+#        if isinstance(other, Number):
+#            return type(self)(*[ el.divide(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+#        elif issubclass(other.__class__, DataContainer):
+#            # try to do algebra with one DataContainer. Will raise error if not compatible
+#            if out is not None:
+#                kw = kwargs.copy()
+#                for i,el in enumerate(self.containers):
+#                    kw['out'] = out.get_item(i)
+#                    el.divide(other, *args, **kw)
+#                return
+#            else:
+#                return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape)
+#        return type(self)(*[ el.divide(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)],
+#                          shape=self.shape)
+    def add(self, other, *args, **kwargs):
+        out = kwargs.get('out', None)
+        if out is not None:
+            self.binary_operations(BlockDataContainer.ADD, other, *args, **kwargs)
+        else:
+            return self.binary_operations(BlockDataContainer.ADD, other, *args, **kwargs)
+    def subtract(self, other, *args, **kwargs):
+        out = kwargs.get('out', None)
+        if out is not None:
+            self.binary_operations(BlockDataContainer.SUBTRACT, other, *args, **kwargs)
+        else:
+            return self.binary_operations(BlockDataContainer.SUBTRACT, other, *args, **kwargs)
+    def multiply(self, other, *args, **kwargs):
+        out = kwargs.get('out', None)
+        if out is not None:
+            self.binary_operations(BlockDataContainer.MULTIPLY, other, *args, **kwargs)
+        else:
+            return self.binary_operations(BlockDataContainer.MULTIPLY, other, *args, **kwargs)
+    def divide(self, other, *args, **kwargs):
+        out = kwargs.get('out', None)
+        if out is not None:
+            self.binary_operations(BlockDataContainer.DIVIDE, other, *args, **kwargs)
+        else:
+            return self.binary_operations(BlockDataContainer.DIVIDE, other, *args, **kwargs)
+    
+
+    def binary_operations(self, operation, other, *args, **kwargs):
+        if not self.is_compatible(other):
+            raise ValueError('Incompatible for divide')
+        out = kwargs.get('out', None)
+        if isinstance(other, Number) or issubclass(other.__class__, DataContainer):
+            # try to do algebra with one DataContainer. Will raise error if not compatible
+            kw = kwargs.copy()
+            res = []
+            for i,el in enumerate(self.containers):
+                if operation == BlockDataContainer.ADD:
+                    op = el.add
+                elif operation == BlockDataContainer.SUBTRACT:
+                    op = el.subtract
+                elif operation == BlockDataContainer.MULTIPLY:
+                    op = el.multiply
+                elif operation == BlockDataContainer.DIVIDE:
+                    op = el.divide
+                elif operation == BlockDataContainer.POWER:
+                    op = el.power
+                else:
+                    raise ValueError('Unsupported operation', operation)
+                if out is not None:
+                    kw['out'] = out.get_item(i)
+                    op(other, *args, **kw)
+                else:
+                    res.append(op(other, *args, **kw))
+            if out is not None:
+                return
+            else:
+                return type(self)(*res, shape=self.shape)
+        elif isinstance(other, (list, numpy.ndarray, BlockDataContainer)):
+            # try to do algebra with one DataContainer. Will raise error if not compatible
+            kw = kwargs.copy()
+            res = []
+            if isinstance(other, BlockDataContainer):
+                the_other = other.containers
+            else:
+                the_other = other
+            for i,zel in enumerate(zip ( self.containers, the_other) ):
+                el = zel[0]
+                ot = zel[1]
+                if operation == BlockDataContainer.ADD:
+                    op = el.add
+                elif operation == BlockDataContainer.SUBTRACT:
+                    op = el.subtract
+                elif operation == BlockDataContainer.MULTIPLY:
+                    op = el.multiply
+                elif operation == BlockDataContainer.DIVIDE:
+                    op = el.divide
+                elif operation == BlockDataContainer.POWER:
+                    op = el.power
+                else:
+                    raise ValueError('Unsupported operation', operation)
+                if out is not None:
+                    kw['out'] = out.get_item(i)
+                    op(ot, *args, **kw)
+                else:
+                    res.append(op(ot, *args, **kw))
+            if out is not None:
+                return
+            else:
+                return type(self)(*res, shape=self.shape)
+            return type(self)(*[ operation(ot, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+        else:
+            raise ValueError('Incompatible type {}'.format(type(other)))
+    
+
+    def power(self, other, *args, **kwargs):
+        if not self.is_compatible(other):
+            raise ValueError('Incompatible for power')
+        out = kwargs.get('out', None)
+        if isinstance(other, Number):
+            return type(self)(*[ el.power(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+        return type(self)(*[ el.power(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)
+    
+    def maximum(self,other, *args, **kwargs):
+        if not self.is_compatible(other):
+            raise ValueError('Incompatible for maximum')
+        out = kwargs.get('out', None)
+        if isinstance(other, Number):
+            return type(self)(*[ el.maximum(other, *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, *args, **kwargs) for el,ot in zip(self.containers,other)], shape=self.shape)
+        return type(self)(*[ el.maximum(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], shape=self.shape)
+    
+    ## unary operations    
+    def abs(self, *args,  **kwargs):
+        return type(self)(*[ el.abs(*args, **kwargs) for el in self.containers], shape=self.shape)
+    def sign(self, *args,  **kwargs):
+        return type(self)(*[ el.sign(*args, **kwargs) for el in self.containers], shape=self.shape)
+    def sqrt(self, *args,  **kwargs):
+        return type(self)(*[ el.sqrt(*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], shape=self.shape)
+    
+    ## reductions
+    def sum(self, *args, **kwargs):
+        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() 
+    def norm(self):
+        return numpy.sqrt(self.squared_norm())    
+    def copy(self):
+        '''alias of clone'''    
+        return self.clone()
+    def clone(self):
+        return type(self)(*[el.copy() for el in self.containers], shape=self.shape)
+    def fill(self, other):
+        if isinstance (other, BlockDataContainer):
+            if not self.is_compatible(other):
+                raise ValueError('Incompatible containers')
+            for el,ot in zip(self.containers, other.containers):
+                el.fill(ot)
+        else:
+            return ValueError('Cannot fill with object provided {}'.format(type(other)))
+    
+    def __add__(self, other):
+        return self.add( other )
+    # __radd__
+    
+    def __sub__(self, other):
+        return self.subtract( other )
+    # __rsub__
+    
+    def __mul__(self, other):
+        return self.multiply(other)
+    # __rmul__
+    
+    def __div__(self, other):
+        return self.divide(other)
+    # __rdiv__
+    def __truediv__(self, other):
+        return self.divide(other)
+    
+    def __pow__(self, other):
+        return self.power(other)
+    # reverse operand
+    def __radd__(self, other):
+        '''Reverse addition
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return self + other
+    # __radd__
+    
+    def __rsub__(self, other):
+        '''Reverse subtraction
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return (-1 * self) + other
+    # __rsub__
+    
+    def __rmul__(self, other):
+        '''Reverse multiplication
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return self * other
+    # __rmul__
+    
+    def __rdiv__(self, other):
+        '''Reverse division
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return pow(self / other, -1)
+    # __rdiv__
+    def __rtruediv__(self, other):
+        '''Reverse truedivision
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return self.__rdiv__(other)
+    
+    def __rpow__(self, other):
+        '''Reverse power
+        
+        to make sure that this method is called rather than the __mul__ of a numpy array
+        the class constant __array_priority__ must be set > 0
+        https://docs.scipy.org/doc/numpy-1.15.1/reference/arrays.classes.html#numpy.class.__array_priority__
+        '''
+        return other.power(self)
+    
+    def __iadd__(self, other):
+        '''Inline addition'''
+        if isinstance (other, BlockDataContainer):
+            for el,ot in zip(self.containers, other.containers):
+                el += ot
+        elif isinstance(other, Number):
+            for el in self.containers:
+                el += other
+        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+            if not self.is_compatible(other):
+                raise ValueError('Incompatible for __iadd__')
+            for el,ot in zip(self.containers, other):
+                el += ot
+        return self
+    # __iadd__
+    
+    def __isub__(self, other):
+        '''Inline subtraction'''
+        if isinstance (other, BlockDataContainer):
+            for el,ot in zip(self.containers, other.containers):
+                el -= ot
+        elif isinstance(other, Number):
+            for el in self.containers:
+                el -= other
+        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+            if not self.is_compatible(other):
+                raise ValueError('Incompatible for __isub__')
+            for el,ot in zip(self.containers, other):
+                el -= ot
+        return self
+    # __isub__
+    
+    def __imul__(self, other):
+        '''Inline multiplication'''
+        if isinstance (other, BlockDataContainer):
+            for el,ot in zip(self.containers, other.containers):
+                el *= ot
+        elif isinstance(other, Number):
+            for el in self.containers:
+                el *= other
+        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+            if not self.is_compatible(other):
+                raise ValueError('Incompatible for __imul__')
+            for el,ot in zip(self.containers, other):
+                el *= ot
+        return self
+    # __imul__
+    
+    def __idiv__(self, other):
+        '''Inline division'''
+        if isinstance (other, BlockDataContainer):
+            for el,ot in zip(self.containers, other.containers):
+                el /= ot
+        elif isinstance(other, Number):
+            for el in self.containers:
+                el /= other
+        elif isinstance(other, list) or isinstance(other, numpy.ndarray):
+            if not self.is_compatible(other):
+                raise ValueError('Incompatible for __idiv__')
+            for el,ot in zip(self.containers, other):
+                el /= ot
+        return self
+    # __rdiv__
+    def __itruediv__(self, other):
+        '''Inline truedivision'''
+        return self.__idiv__(other)
+
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 5323e76..439149c 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -67,24 +67,22 @@ class PDHG(Algorithm):
 
             #self.y = self.f.proximal_conjugate(self.y_old, self.sigma)
             self.f.proximal_conjugate(self.y_old, self.sigma, out=self.y)
-            
+
             # Gradient ascent, Primal problem solution
             self.operator.adjoint(self.y, out=self.x_tmp)
             self.x_tmp *= self.tau
             self.x_old -= self.x_tmp
-            
+
             self.g.proximal(self.x_old, self.tau, out=self.x)
-            
+
             #Update
             self.x.subtract(self.x_old, out=self.xbar)
-            #self.xbar -= self.x_old 
             self.xbar *= self.theta
             self.xbar += self.x
-                            
+
             self.x_old.fill(self.x)
             self.y_old.fill(self.y)
-            #self.y_old = self.y.copy()
-            #self.x_old = self.x.copy()
+
         else:
             # Gradient descent, Dual problem solution
             self.y_old += self.sigma * self.operator.direct(self.xbar)
@@ -93,19 +91,16 @@ class PDHG(Algorithm):
             # Gradient ascent, Primal problem solution
             self.x_old -= self.tau * self.operator.adjoint(self.y)
             self.x = self.g.proximal(self.x_old, self.tau)
-            
+
             #Update
             #xbar = x + theta * (x - x_old)
             self.xbar.fill(self.x)
             self.xbar -= self.x_old 
             self.xbar *= self.theta
             self.xbar += self.x
-                            
-            self.x_old.fill(self.x)
-            self.y_old.fill(self.y)
-            #self.y_old = self.y.copy()
-            #self.x_old = self.x.copy()
-            #self.y = self.y_old
+
+            self.x_old = self.x
+            self.y_old = self.y
 
     def update_objective(self):
         p1 = self.f(self.operator.direct(self.x)) + self.g(self.x)
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
index 0c7eb85..24c47f4 100755
--- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -21,6 +21,7 @@ import numpy as np
 from ccpi.optimisation.functions import Function, ScaledFunction
 from ccpi.framework import DataContainer, ImageData, \
                            ImageGeometry, BlockDataContainer 
+import functools
 
 ############################   mixed_L1,2NORM FUNCTIONS   #####################
 class MixedL21Norm(Function):
@@ -36,7 +37,9 @@ class MixedL21Norm(Function):
             
             :param: x is a BlockDataContainer
                                 
-        '''                        
+        '''
+        if not isinstance(x, BlockDataContainer):
+            raise ValueError('__call__ expected BlockDataContainer, got {}'.format(type(x)))                      
         if self.SymTensor:
             
             param = [1]*x.shape[0]
@@ -73,7 +76,6 @@ class MixedL21Norm(Function):
             different form L2NormSquared which acts on DC
         
         '''
-        
         pass
     
     def proximal_conjugate(self, x, tau, out=None): 
@@ -90,17 +92,18 @@ class MixedL21Norm(Function):
         
         else:
 
-                              
-         if out is None:
-            
-            tmp = [ el*el for el in x.containers]
-            res = (sum(tmp).sqrt()).maximum(1.0)                                                                
-            return x/res
-        
-         else:
-                         
-            res = (sum(x**2).sqrt()).maximum(1.0)             
-            out.fill(x/res)
+            if out is None:                                        
+                tmp = [ el*el for el in x.containers]
+                res = sum(tmp).sqrt().maximum(1.0) 
+                frac = [el/res for el in x.containers]
+                res = BlockDataContainer(*frac)    
+                return res
+            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)
+                #for i,el in enumerate(x.containers):
+                #    el.divide(res, out=out.get_item(i))
 
     def __rmul__(self, scalar):
         return ScaledFunction(self, scalar) 
diff --git a/Wrappers/Python/test/test_BlockDataContainer.py b/Wrappers/Python/test/test_BlockDataContainer.py
index 2fca23c..0dd0657 100755
--- a/Wrappers/Python/test/test_BlockDataContainer.py
+++ b/Wrappers/Python/test/test_BlockDataContainer.py
@@ -136,12 +136,12 @@ class TestBlockDataContainer(unittest.TestCase):
         print  (a[0][0].shape)

         #cp2 = BlockDataContainer(*a)

         cp2 = cp0.add(cp1)

-        assert (cp2.get_item(0).as_array()[0][0][0] == 2.)

-        assert (cp2.get_item(1).as_array()[0][0][0] == 4.)

+        self.assertEqual (cp2.get_item(0).as_array()[0][0][0] , 2.)

+        self.assertEqual (cp2.get_item(1).as_array()[0][0][0] , 4.)

         

         cp2 = cp0 + cp1 

-        assert (cp2.get_item(0).as_array()[0][0][0] == 2.)

-        assert (cp2.get_item(1).as_array()[0][0][0] == 4.)

+        self.assertTrue (cp2.get_item(0).as_array()[0][0][0] == 2.)

+        self.assertTrue (cp2.get_item(1).as_array()[0][0][0] == 4.)

         cp2 = cp0 + 1 

         numpy.testing.assert_almost_equal(cp2.get_item(0).as_array()[0][0][0] , 1. , decimal=5)

         numpy.testing.assert_almost_equal(cp2.get_item(1).as_array()[0][0][0] , 2., decimal = 5)

@@ -427,6 +427,21 @@ class TestBlockDataContainer(unittest.TestCase):
         cp0.fill(cp2)

         self.assertBlockDataContainerEqual(cp0, cp2)

 

+    def test_NestedBlockDataContainer(self):

+        ig0 = ImageGeometry(2,3,4)

+        ig1 = ImageGeometry(2,3,5)

+        

+        data0 = ig0.allocate(0)

+        data2 = ig0.allocate(1)

+        

+        cp0 = BlockDataContainer(data0,data2)

+        #cp1 = BlockDataContainer(data2,data3)

+

+        nested = BlockDataContainer(cp0, data2, data2)

+        out = BlockDataContainer(BlockDataContainer(data0 , data0), data0, data0)

+        nested.divide(data2,out=out)

+        self.assertBlockDataContainerEqual(out, nested)

+

 

     def assertBlockDataContainerEqual(self, container1, container2):

         print ("assert Block Data Container Equal")

diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py
index 1891afd..bc1f034 100644
--- a/Wrappers/Python/test/test_functions.py
+++ b/Wrappers/Python/test/test_functions.py
@@ -20,7 +20,7 @@ from ccpi.optimisation.operators import Gradient
 #from ccpi.optimisation.functions import SimpleL2NormSq
 from ccpi.optimisation.functions import L2NormSquared
 #from ccpi.optimisation.functions import SimpleL1Norm
-from ccpi.optimisation.functions import L1Norm
+from ccpi.optimisation.functions import L1Norm, MixedL21Norm
 
 from ccpi.optimisation.funcs import Norm2sq
 # from ccpi.optimisation.functions.L2NormSquared import SimpleL2NormSq, L2NormSq
@@ -29,13 +29,46 @@ from ccpi.optimisation.funcs import Norm2sq
 from ccpi.optimisation.functions import ZeroFun
 
 from ccpi.optimisation.functions import FunctionOperatorComposition
-import unittest  
-import numpy      
+import unittest
+import numpy
 
 #    
 
                     
 class TestFunction(unittest.TestCase):
+    def assertBlockDataContainerEqual(self, container1, container2):
+        print ("assert Block Data Container Equal")
+        self.assertTrue(issubclass(container1.__class__, container2.__class__))
+        for col in range(container1.shape[0]):
+            if issubclass(container1.get_item(col).__class__, DataContainer):
+                print ("Checking col ", col)
+                self.assertNumpyArrayEqual(
+                    container1.get_item(col).as_array(), 
+                    container2.get_item(col).as_array()
+                    )
+            else:
+                self.assertBlockDataContainerEqual(container1.get_item(col),container2.get_item(col))
+    
+    def assertNumpyArrayEqual(self, first, second):
+        res = True
+        try:
+            numpy.testing.assert_array_equal(first, second)
+        except AssertionError as err:
+            res = False
+            print(err)
+        self.assertTrue(res)
+
+    def assertNumpyArrayAlmostEqual(self, first, second, decimal=6):
+        res = True
+        try:
+            numpy.testing.assert_array_almost_equal(first, second, decimal)
+        except AssertionError as err:
+            res = False
+            print(err)
+            print("expected " , second)
+            print("actual " , first)
+
+        self.assertTrue(res)
     def test_Function(self):
     
     
@@ -280,8 +313,37 @@ class TestFunction(unittest.TestCase):
         ynew = new_chisq.gradient(u)
         numpy.testing.assert_array_equal(yold.as_array(), ynew.as_array())
         
+    def test_mixedL12Norm(self):
+        M, N, K = 2,3,5
+        ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
+        u1 = ig.allocate('random_int')
+        u2 = ig.allocate('random_int')
+        
+        U = BlockDataContainer(u1, u2, shape=(2,1))
+        
+        # Define no scale and scaled
+        f_no_scaled = MixedL21Norm() 
+        #f_scaled = 0.5 * MixedL21Norm()  
+        
+        # call
+        
+        # a1 = f_no_scaled(U)
+        # a2 = f_scaled(U)
+        # self.assertBlockDataContainerEqual(a1,a2)
+        tmp = [ el**2 for el in U.containers ]
+        self.assertBlockDataContainerEqual(BlockDataContainer(*tmp),
+                                           U.power(2))
+        
+        z1 = f_no_scaled.proximal_conjugate(U, 1)
+        u3 = ig.allocate('random_int')
+        u4 = ig.allocate('random_int')
+        
+        z3 = BlockDataContainer(u3, u4, shape=(2,1))
+        
+        
+        f_no_scaled.proximal_conjugate(U, 1, out=z3)
+        self.assertBlockDataContainerEqual(z3,z1)
 
-    
 #    
 #    f1 = L2NormSq(alpha=1, b=noisy_data)
 #    print(f1(noisy_data))