diff options
5 files changed, 164 insertions, 76 deletions
diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index cbce844..8152bff 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -16,17 +16,24 @@ import functools #from ccpi.optimisation.operators import Operator, LinearOperator
class BlockDataContainer(object):
- '''Class to hold DataContainers as blocks'''
+ '''Class to hold DataContainers as column vector'''
__array_priority__ = 1
def __init__(self, *args, **kwargs):
- '''containers must be passed row by row'''
+ '''containers must be consistent in shape'''
self.containers = args
+ for i, co in enumerate(args):
+ if i == 0:
+ shape = co.shape
+ else:
+ if shape != co.shape:
+ raise ValueError('Expected shape is {} got {}'.format(shape, co.shape))
self.index = 0
- shape = kwargs.get('shape', None)
- if shape is None:
- shape = (len(args),1)
+ #shape = kwargs.get('shape', None)
+ #if shape is None:
+ # shape = (len(args),1)
+ shape = (len(args),1)
self.shape = shape
- print (self.shape)
+ #print (self.shape)
n_elements = functools.reduce(lambda x,y: x*y, shape, 1)
if len(args) != n_elements:
raise ValueError(
@@ -65,14 +72,12 @@ class BlockDataContainer(object): elif isinstance(other, numpy.ndarray):
return self.shape == other.shape
return len(self.containers) == len(other.containers)
- def get_item(self, row, col=0):
+ def get_item(self, row):
if row > self.shape[0]:
raise ValueError('Requested row {} > max {}'.format(row, self.shape[0]))
- if col > self.shape[1]:
- raise ValueError('Requested col {} > max {}'.format(col, self.shape[1]))
-
- index = row*self.shape[1]+col
- return self.containers[index]
+ return self.containers[row]
+ def __getitem__(self, row):
+ return self.get_item(row)
def add(self, other, *args, **kwargs):
assert self.is_compatible(other)
@@ -96,6 +101,7 @@ class BlockDataContainer(object): shape=self.shape)
def multiply(self, other, *args, **kwargs):
+ print ("BlockDataContainer" , other)
self.is_compatible(other)
out = kwargs.get('out', None)
if isinstance(other, Number):
@@ -302,8 +308,8 @@ class BlockDataContainer(object): def __itruediv__(self, other):
'''Inline truedivision'''
return self.__idiv__(other)
- @property
- def T(self):
- '''return the transposed of self'''
- shape = (self.shape[1], self.shape[0])
- return type(self)(*self.containers, shape=shape)
+ #@property
+ #def T(self):
+ # '''return the transposed of self'''
+ # shape = (self.shape[1], self.shape[0])
+ # return type(self)(*self.containers, shape=shape)
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index a83dc8a..a3ba93e 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -10,7 +10,7 @@ from numbers import Number import functools from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer from ccpi.optimisation.operators import Operator, LinearOperator - +from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator class BlockOperator(Operator): @@ -18,7 +18,13 @@ class BlockOperator(Operator): Class to hold a number of Operators in a block. User may specify the shape of the block, by default is a row vector + + BlockOperators have a generic shape M x N, and when applied on an + Nx1 BlockDataContainer, will yield and Mx1 BlockDataContainer. + Notice: BlockDatacontainer are only allowed to have the shape of N x 1, with + N rows and 1 column. ''' + __array_priority__ = 1 def __init__(self, *args, **kwargs): ''' Class creator @@ -66,11 +72,9 @@ class BlockOperator(Operator): def direct(self, x, out=None): shape = self.get_output_shape(x.shape) - print ("direct output shape", shape) res = [] for row in range(self.shape[0]): for col in range(self.shape[1]): - print ("row {} col {}".format(row, col)) if col == 0: prod = self.get_item(row,col).direct(x.get_item(col)) else: @@ -79,7 +83,6 @@ class BlockOperator(Operator): return BlockDataContainer(*res, shape=shape) def get_output_shape(self, xshape, adjoint=False): - print ("get_output_shape", self.shape, xshape) sshape = self.shape[1] oshape = self.shape[0] if adjoint: @@ -104,7 +107,9 @@ class BlockOperator(Operator): scalars = [scalar for _ in self.operators] # create a list of ScaledOperator-s ops = [ v * op for v,op in zip(scalars, self.operators)] - return BlockOperator(*ops, shape=self.shape) + #return BlockScaledOperator(self, scalars ,shape=self.shape) + return type(self)(*ops, shape=self.shape) + @property def T(self): '''Return the transposed of self''' shape = (self.shape[1], self.shape[0]) diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py index a47bec2..aeb6c53 100644 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockScaledOperator.py @@ -1,22 +1,67 @@ -# -*- coding: utf-8 -*- -""" -Created on Thu Feb 14 12:36:40 2019 - -@author: ofn77899 -""" +from numbers import Number import numpy -from ccpi.optimisation.operators import BlockOperator +from ccpi.optimisation.operators import ScaledOperator +import functools + +class BlockScaledOperator(ScaledOperator): + '''ScaledOperator + A class to represent the scalar multiplication of an Operator with a scalar. + It holds an operator and a scalar. Basically it returns the multiplication + of the result of direct and adjoint of the operator with the scalar. + For the rest it behaves like the operator it holds. - -class BlockScaledOperator(BlockOperator): - def __init__(self, *args, **kwargs): - super(BlockScaledOperator, self).__init__(*args, **kwargs) - scalar = kwargs.get('scalar',1) - if isinstance (scalar, list) or isinstance(scalar, tuple) or \ - isinstance(scalar, numpy.ndarray): - if len(scalars) != len(self.operators): - raise ValueError('dimensions of scalars and operators do not match') + Args: + operator (Operator): a Operator or LinearOperator + scalar (Number): a scalar multiplier + Example: + The scaled operator behaves like the following: + sop = ScaledOperator(operator, scalar) + sop.direct(x) = scalar * operator.direct(x) + sop.adjoint(x) = scalar * operator.adjoint(x) + sop.norm() = operator.norm() + sop.range_geometry() = operator.range_geometry() + sop.domain_geometry() = operator.domain_geometry() + ''' + def __init__(self, operator, scalar, shape=None): + if shape is None: + shape = operator.shape + + if isinstance(scalar, (list, tuple, numpy.ndarray)): + size = functools.reduce(lambda x,y:x*y, shape, 1) + if len(scalar) != size: + raise ValueError('Scalar and operators size do not match: {}!={}' + .format(len(scalar), len(operator))) + self.scalar = scalar[:] + print ("BlockScaledOperator ", self.scalar) + elif isinstance (scalar, Number): + self.scalar = scalar + else: + raise TypeError('expected scalar to be a number of an iterable: got {}'.format(type(scalar))) + self.operator = operator + self.shape = shape + def direct(self, x, out=None): + print ("BlockScaledOperator self.scalar", self.scalar) + #print ("self.scalar", self.scalar[0]* x.get_item(0).as_array()) + return self.scalar * (self.operator.direct(x, out=out)) + def adjoint(self, x, out=None): + if self.operator.is_linear(): + return self.scalar * self.operator.adjoint(x, out=out) else: - scalar = [scalar for _ in self.operators] - self.operators = [v * op for op in self.operators]
\ No newline at end of file + raise TypeError('Operator is not linear') + def norm(self): + return numpy.abs(self.scalar) * self.operator.norm() + def range_geometry(self): + return self.operator.range_geometry() + def domain_geometry(self): + return self.operator.domain_geometry() + @property + def T(self): + '''Return the transposed of self''' + #print ("transpose before" , self.shape) + #shape = (self.shape[1], self.shape[0]) + ##self.shape = shape + ##self.operator.shape = shape + #print ("transpose" , shape) + #return self + return type(self)(self.operator.T, self.scalar)
\ No newline at end of file diff --git a/Wrappers/Python/test/test_BlockDataContainer.py b/Wrappers/Python/test/test_BlockDataContainer.py index c14a101..190fb21 100755 --- a/Wrappers/Python/test/test_BlockDataContainer.py +++ b/Wrappers/Python/test/test_BlockDataContainer.py @@ -19,7 +19,7 @@ from ccpi.framework import BlockDataContainer import functools
class TestBlockDataContainer(unittest.TestCase):
- def test_BlockDataContainerShape(self):
+ def skiptest_BlockDataContainerShape(self):
print ("test block data container")
ig0 = ImageGeometry(2,3,4)
ig1 = ImageGeometry(12,42,55,32)
@@ -34,7 +34,7 @@ class TestBlockDataContainer(unittest.TestCase): cp1 = BlockDataContainer(data2,data3)
transpose_shape = (cp0.shape[1], cp0.shape[0])
self.assertTrue(cp0.T.shape == transpose_shape)
- def test_BlockDataContainerShapeArithmetic(self):
+ def skiptest_BlockDataContainerShapeArithmetic(self):
print ("test block data container")
ig0 = ImageGeometry(2,3,4)
ig1 = ImageGeometry(12,42,55,32)
diff --git a/Wrappers/Python/test/test_BlockOperator.py b/Wrappers/Python/test/test_BlockOperator.py index f3b057b..1896978 100644 --- a/Wrappers/Python/test/test_BlockOperator.py +++ b/Wrappers/Python/test/test_BlockOperator.py @@ -8,64 +8,96 @@ import numpy class TestBlockOperator(unittest.TestCase): def test_BlockOperator(self): + print ("test_BlockOperator") ig = [ ImageGeometry(10,20,30) , \ - ImageGeometry(11,21,31) , \ - ImageGeometry(12,22,32) ] + ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) ] x = [ g.allocate() for g in ig ] ops = [ TomoIdentity(g) for g in ig ] K = BlockOperator(*ops) - #X = BlockDataContainer(*x).T + 1 X = BlockDataContainer(x[0]) Y = K.direct(X) - #self.assertTrue(Y.shape == X.shape) + self.assertTrue(Y.shape == K.shape) numpy.testing.assert_array_equal(Y.get_item(0).as_array(),X.get_item(0).as_array()) numpy.testing.assert_array_equal(Y.get_item(1).as_array(),X.get_item(0).as_array()) #numpy.testing.assert_array_equal(Y.get_item(2).as_array(),X.get_item(2).as_array()) - + + X = BlockDataContainer(*x) + 1 + Y = K.T.direct(X) + # K.T (1,3) X (3,1) => output shape (1,1) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),len(x)+zero) + def test_ScaledBlockOperatorSingleScalar(self): ig = [ ImageGeometry(10,20,30) , \ - ImageGeometry(11,21,31) , \ - ImageGeometry(12,22,32) ] + ImageGeometry(10,20,30) , \ + ImageGeometry(10,20,30) ] x = [ g.allocate() for g in ig ] ops = [ TomoIdentity(g) for g in ig ] + val = 1 + # test limit as non Scaled + scalar = 1 + k = BlockOperator(*ops) + K = scalar * k + X = BlockDataContainer(*x) + val + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),((scalar*xx).sum()+zero)) + scalar = 0.5 - K = 0.5 * BlockOperator(*ops) + k = BlockOperator(*ops) + K = scalar * k X = BlockDataContainer(*x) + 1 - print (X.shape) - X = BlockDataContainer(*x).T + 1 - print (X.shape, K.shape) - Y = K.direct(X) - self.assertTrue(Y.shape == X.shape) - - numpy.testing.assert_array_equal(Y.get_item(0).as_array(),scalar * X.get_item(0).as_array()) - numpy.testing.assert_array_equal(Y.get_item(1).as_array(),scalar * X.get_item(1).as_array()) - numpy.testing.assert_array_equal(Y.get_item(2).as_array(),scalar * X.get_item(2).as_array()) - + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),scalar*(len(x)+zero)) + + def test_ScaledBlockOperatorScalarList(self): - ig = [ImageGeometry(10, 20, 30), - ImageGeometry(11, 21, 31), - ImageGeometry(12, 22, 32)] - x = [g.allocate() for g in ig] - ops = [TomoIdentity(g) for g in ig] + ig = [ ImageGeometry(2,3) , \ + #ImageGeometry(10,20,30) , \ + ImageGeometry(2,3 ) ] + x = [ g.allocate() for g in ig ] + ops = [ TomoIdentity(g) for g in ig ] - scalar = [i*1.2 for i, el in enumerate(ig)] - K = scalar * BlockOperator(*ops) - X = BlockDataContainer(*x).T + 1 - Y = K.direct(X) - self.assertTrue(Y.shape == X.shape) + # test limit as non Scaled + scalar = numpy.asarray([1 for _ in x]) + k = BlockOperator(*ops) + K = scalar * k + val = 1 + X = BlockDataContainer(*x) + val + + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(),(scalar*xx).sum()+zero) + + scalar = numpy.asarray([i+1 for i,el in enumerate(x)]) + #scalar = numpy.asarray([6,0]) + k = BlockOperator(*ops) + K = scalar * k + X = BlockDataContainer(*x) + val + Y = K.T.direct(X) + self.assertTrue(Y.shape == (1,1)) + zero = numpy.zeros(X.get_item(0).shape) + xx = numpy.asarray([val for _ in x]) + numpy.testing.assert_array_equal(Y.get_item(0).as_array(), - scalar[0] * X.get_item(0).as_array()) - numpy.testing.assert_array_equal(Y.get_item(1).as_array(), - scalar[1] * X.get_item(1).as_array()) - numpy.testing.assert_array_equal(Y.get_item(2).as_array(), - scalar[2] * X.get_item(2).as_array()) - + (scalar*xx).sum()+zero) + def test_TomoIdentity(self): ig = ImageGeometry(10,20,30) |