From 1dec48d390df5cbc3436832cedf559b64f4651bc Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 11 Apr 2019 15:43:31 +0100 Subject: fix out call --- .../ccpi/optimisation/functions/MixedL21Norm.py | 58 +++++++++++++--------- 1 file changed, 34 insertions(+), 24 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index ed1d5e5..c6b6e95 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): @@ -88,29 +90,37 @@ class MixedL21Norm(Function): return res else: -# pass + 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) + + if False: + # works but not memory efficient as allocating a new BlockDataContainer + a = x / res + out.fill(a) + elif False: + # this leads to error +# File "ccpi\framework\BlockDataContainer.py", line 142, in divide +# return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape) +# File "ccpi\framework\BlockDataContainer.py", line 142, in +# return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape) +# File "ccpi\framework\framework.py", line 814, in divide +# return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs) +# File "ccpi\framework\framework.py", line 802, in pixel_wise_binary +# raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) +# ValueError: ImageData: incompatible class: true_divide + x.divide(res, out=out) + else: + for i,el in enumerate(x.containers): + #a = out.get_item(i) + el.divide(res, out=out.get_item(i)) - -# # tmp2 = np.sqrt(x.as_array()[0]**2 + x.as_array()[1]**2 + 2*x.as_array()[2]**2)/self.alpha -# # res = x.divide(ImageData(tmp2).maximum(1.0)) -# if out is None: - - tmp = [ el*el for el in x] - res = (sum(tmp).sqrt()).maximum(1.0) - frac = [x[i]/res for i in range(x.shape[0])] - res = BlockDataContainer(*frac) - - return res - # else: - # tmp = [ el*el for el in x] - # res = (sum(tmp).sqrt()).maximum(1.0) - # #frac = [x[i]/res for i in range(x.shape[0])] - # for i in range(x.shape[0]): - # a = out.get_item(i) - # b = x.get_item(i) - # b /= res - # a.fill( b ) - def __rmul__(self, scalar): return ScaledFunction(self, scalar) -- cgit v1.2.3 From e80f8c108871245f06dc3e570502e95a4acba64b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 11 Apr 2019 15:43:54 +0100 Subject: clean code closes #240 --- .../ccpi/optimisation/functions/MixedL21Norm.py | 23 ++-------------------- 1 file changed, 2 insertions(+), 21 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index c6b6e95..a655e03 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -99,27 +99,8 @@ class MixedL21Norm(Function): else: res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) res = res1.sqrt().maximum(1.0) - - if False: - # works but not memory efficient as allocating a new BlockDataContainer - a = x / res - out.fill(a) - elif False: - # this leads to error -# File "ccpi\framework\BlockDataContainer.py", line 142, in divide -# return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape) -# File "ccpi\framework\BlockDataContainer.py", line 142, in -# return type(self)(*[ el.divide(other, *args, **kwargs) for el in self.containers], shape=self.shape) -# File "ccpi\framework\framework.py", line 814, in divide -# return self.pixel_wise_binary(numpy.divide, other, *args, **kwargs) -# File "ccpi\framework\framework.py", line 802, in pixel_wise_binary -# raise ValueError (message(type(self), "incompatible class:" , pwop.__name__, type(out))) -# ValueError: ImageData: incompatible class: true_divide - x.divide(res, out=out) - else: - for i,el in enumerate(x.containers): - #a = out.get_item(i) - el.divide(res, out=out.get_item(i)) + for i,el in enumerate(x.containers): + el.divide(res, out=out.get_item(i)) def __rmul__(self, scalar): return ScaledFunction(self, scalar) -- cgit v1.2.3 From 350a889c38805dcda98a299315af1ab64510fa5b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 11 Apr 2019 15:45:56 +0100 Subject: add test for mixed L21 Norm --- Wrappers/Python/test/test_functions.py | 70 ++++++++++++++++++++++++++++++++-- 1 file changed, 66 insertions(+), 4 deletions(-) (limited to 'Wrappers') 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)) -- cgit v1.2.3 From 6ce64e15b13cf7c6ae55cf9bc891980679268ac4 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 11 Apr 2019 15:49:24 +0100 Subject: minor code beautification --- .../Python/ccpi/optimisation/algorithms/PDHG.py | 25 ++++++++++------------ 1 file changed, 11 insertions(+), 14 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 086e322..2ac3eba 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -66,24 +66,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) @@ -92,19 +90,18 @@ 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.fill(self.x) + #self.y_old.fill(self.y) + 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) -- cgit v1.2.3 From 10aae87e1416d291906b94927acb4aac5737a44e Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 11 Apr 2019 17:13:08 +0100 Subject: fixing algebra with nested block data containers --- .../Python/ccpi/framework/BlockDataContainer.py | 45 +++++++++++++++++++--- 1 file changed, 40 insertions(+), 5 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 13663c2..85cd05a 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -53,9 +53,9 @@ class BlockDataContainer(object): 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] + #for i in range(len(self.containers)): + # if type(self.containers[i])==type(self): + # self = self.containers[i] if isinstance(other, Number): return True @@ -71,7 +71,16 @@ class BlockDataContainer(object): elif isinstance(other, numpy.ndarray): return len(self.containers) == len(other) elif issubclass(other.__class__, DataContainer): - return self.get_item(0).shape == other.shape + 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): @@ -139,10 +148,36 @@ class BlockDataContainer(object): 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) + 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 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 + 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) + 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) 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') -- cgit v1.2.3 From 4129fe6a81d20b5f7d05554222d3a78f18f014f0 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 12 Apr 2019 11:56:51 +0100 Subject: convert eol to unix --- .../Python/ccpi/framework/BlockDataContainer.py | 827 +++++++++++---------- 1 file changed, 450 insertions(+), 377 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 85cd05a..fee0cda 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -1,377 +1,450 @@ - # -*- 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 - #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 - 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 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 - 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) - 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) - 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 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)): + # try to do algebra with one DataContainer. Will raise error if not compatible + kw = kwargs.copy() + res = [] + for i,zel in enumerate(zip ( self.containers, 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) + elif isinstance(other, BlockDataContainer): + return type(self)(*[ el.divide(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], + 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) + -- cgit v1.2.3 From 4049917fa20353da6f88ff03b7c5b11d67743a00 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 12 Apr 2019 12:09:39 +0100 Subject: fix BlockDataContainer algebra closes #242 --- .../Python/ccpi/framework/BlockDataContainer.py | 152 ++++++++++++--------- 1 file changed, 86 insertions(+), 66 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index fee0cda..529a1ce 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -111,79 +111,98 @@ class BlockDataContainer(object): 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): - 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) - + 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): - 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) - + 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): - 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) + 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): @@ -215,11 +234,15 @@ class BlockDataContainer(object): return else: return type(self)(*res, shape=self.shape) - elif isinstance(other, (list, numpy.ndarray)): + elif isinstance(other, (list, numpy.ndarray, BlockDataContainer)): # try to do algebra with one DataContainer. Will raise error if not compatible kw = kwargs.copy() res = [] - for i,zel in enumerate(zip ( self.containers, other) ): + 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: @@ -244,9 +267,6 @@ class BlockDataContainer(object): 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) - elif isinstance(other, BlockDataContainer): - return type(self)(*[ el.divide(ot, *args, **kwargs) for el,ot in zip(self.containers,other.containers)], - shape=self.shape) else: raise ValueError('Incompatible type {}'.format(type(other))) -- cgit v1.2.3 From 04c3f99648e38756e8180519db5f32ac3344ea2b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 12 Apr 2019 12:10:46 +0100 Subject: update test --- Wrappers/Python/test/test_BlockDataContainer.py | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) (limited to 'Wrappers') 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") -- cgit v1.2.3 From 60f73fe79526dc61d721b2dc76a942f2e9a082d1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 12 Apr 2019 12:11:49 +0100 Subject: uses new Algebra of BDC --- Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index a655e03..c5be084 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -99,8 +99,9 @@ class MixedL21Norm(Function): else: res1 = functools.reduce(lambda a,b: a + b*b, x.containers, x.get_item(0) * 0 ) res = res1.sqrt().maximum(1.0) - for i,el in enumerate(x.containers): - el.divide(res, out=out.get_item(i)) + 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) -- cgit v1.2.3 From 7fc291f6d2d71b0d5aa7f3fcf11966910dcea7ab Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 12 Apr 2019 12:12:25 +0100 Subject: removes copy in non optimised alg --- Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py | 2 -- 1 file changed, 2 deletions(-) (limited to 'Wrappers') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 2ac3eba..835c979 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -98,8 +98,6 @@ class PDHG(Algorithm): self.xbar *= self.theta self.xbar += self.x - #self.x_old.fill(self.x) - #self.y_old.fill(self.y) self.x_old = self.x self.y_old = self.y -- cgit v1.2.3