diff options
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py | 272 |
1 files changed, 59 insertions, 213 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index 4bbc536..b2af8fc 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -27,11 +27,14 @@ class BlockOperator(Operator): Do not include the `self` parameter in the ``Args`` section. Args: - vararg (Operator): Operators in the block. varargs are passed in a - row-by-row fashion. - shape (:obj:`tuple`, optional): If passed the Operators listed - in the vararg are laid out as described. Shape and number - of Operators must match. + vararg (Operator): Operators in the block. + shape (:obj:`tuple`, optional): If shape is passed the Operators in + vararg are considered input in a row-by-row fashion. + Shape and number of Operators must match. + + Example: + BlockOperator(op0,op1) results in a row block + BlockOperator(op0,op1,shape=(1,2)) results in a column block ''' self.operators = args shape = kwargs.get('shape', None) @@ -95,33 +98,63 @@ class BlockOperator(Operator): raise ValueError('Incompatible shapes {} {}'.format(self.shape, xshape)) return (oshape, xshape[-1]) -''' - def direct(self, x, out=None): - - out = [None]*self.dimension[0] - for i in range(self.dimension[0]): - z1 = ImageData(np.zeros(self.compMat[i][0].range_dim())) - for j in range(self.dimension[1]): - z1 += self.compMat[i][j].direct(x[j]) - out[i] = z1 - - return out - - def adjoint(self, x, out=None): + +class BlockLinearOperator(BlockOperator): + '''Class to hold a block 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 + ''' + def __init__(self, *args, **kwargs): + ''' + Class creator + + Note: + Do not include the `self` parameter in the ``Args`` section. + + Args: + vararg (Operator): LinearOperators in the block. + shape (:obj:`tuple`, optional): If shape is passed the Operators in + vararg are considered input in a row-by-row fashion. + Shape and number of Operators must match. + + Example: + BlockLinearOperator(op0,op1) results in a row block + BlockLinearOperator(op0,op1,shape=(1,2)) results in a column block + ''' + for i,op in enumerate(args): + if not op.is_linear(): + raise ValueError('Operator {} must be LinearOperator'.format(i)) + super(BlockLinearOperator, self).__init__(*args, **kwargs) + + def adjoint(self, x, out=None): + '''Adjoint operation for the BlockOperator - out = [None]*self.dimension[1] - for i in range(self.dimension[1]): - z2 = ImageData(np.zeros(self.compMat[0][i].domain_dim())) - for j in range(self.dimension[0]): - z2 += self.compMat[j][i].adjoint(x[j]) - out[i] = z2 -''' -from ccpi.optimisation.algorithms import CGLS + only available on BlockLinearOperator + ''' + shape = self.get_output_shape(x.shape, adjoint=True) + res = [] + for row in range(self.shape[1]): + for col in range(self.shape[0]): + if col == 0: + prod = self.get_item(col,row).adjoint(x.get_item(row)) + else: + prod += self.get_item(col,row).adjoint(x.get_item(row)) + res.append(prod) + return BlockDataContainer(*res, shape=shape) + + + + + + if __name__ == '__main__': #from ccpi.optimisation.Algorithms import GradientDescent + from ccpi.optimisation.algorithms import CGLS + from ccpi.plugins.ops import CCPiProjectorSimple from ccpi.optimisation.ops import PowerMethodNonsquare from ccpi.optimisation.ops import TomoIdentity @@ -131,193 +164,6 @@ if __name__ == '__main__': #from ccpi.optimisation.Algorithms import CGLS import matplotlib.pyplot as plt - ig0 = ImageGeometry(2,3,4) - ig1 = ImageGeometry(12,42,55,32) - - data0 = ImageData(geometry=ig0) - data1 = ImageData(geometry=ig1) + 1 - - data2 = ImageData(geometry=ig0) + 2 - data3 = ImageData(geometry=ig1) + 3 - - cp0 = BlockDataContainer(data0,data1) - cp1 = BlockDataContainer(data2,data3) -# - a = [ (el, ot) for el,ot in zip(cp0.containers,cp1.containers)] - print (a[0][0].shape) - #cp2 = BlockDataContainer(*a) - cp2 = cp0.add(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.) - - cp2 = cp0 + cp1 - assert (cp2.get_item(0,0).as_array()[0][0][0] == 2.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == 4.) - cp2 = cp0 + 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - cp2 = cp0 + [1 ,2] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 3., decimal = 5) - cp2 += cp1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +3. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5) - - cp2 += 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , +4. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +7., decimal = 5) - - cp2 += [-2,-1] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 2. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 6., decimal = 5) - - - cp2 = cp0.subtract(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.) - cp2 = cp0 - cp1 - assert (cp2.get_item(0,0).as_array()[0][0][0] == -2.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == -2.) - - cp2 = cp0 - 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0, decimal = 5) - cp2 = cp0 - [1 ,2] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -1. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -1., decimal = 5) - - cp2 -= cp1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -3. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5) - - cp2 -= 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -4. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -5., decimal = 5) - - cp2 -= [-2,-1] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -2. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -4., decimal = 5) - - - cp2 = cp0.multiply(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.) - cp2 = cp0 * cp1 - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - assert (cp2.get_item(1,0).as_array()[0][0][0] == 3.) - - cp2 = cp0 * 2 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5) - cp2 = 2 * cp0 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2, decimal = 5) - cp2 = cp0 * [3 ,2] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - cp2 = cp0 * numpy.asarray([3 ,2]) - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - - cp2 = [3,2] * cp0 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - cp2 = numpy.asarray([3,2]) * cp0 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - cp2 = [3,2,3] * cp0 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 2., decimal = 5) - - cp2 *= cp1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0 , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5) - - cp2 *= 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , +6., decimal = 5) - - cp2 *= [-2,-1] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -6., decimal = 5) - - - cp2 = cp0.divide(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4) - cp2 = cp0/cp1 - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1./3., decimal=4) - - cp2 = cp0 / 2 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5) - cp2 = cp0 / [3 ,2] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5) - cp2 = cp0 / numpy.asarray([3 ,2]) - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5) - cp3 = numpy.asarray([3 ,2]) / (cp0+1) - numpy.testing.assert_almost_equal(cp3.get_item(0,0).as_array()[0][0][0] , 3. , decimal=5) - numpy.testing.assert_almost_equal(cp3.get_item(1,0).as_array()[0][0][0] , 1, decimal = 5) - - cp2 += 1 - cp2 /= cp1 - # TODO fix inplace division - - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 1./2 , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1.5/3., decimal = 5) - - cp2 /= 1 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0.5 , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 0.5, decimal = 5) - - cp2 /= [-2,-1] - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , -0.5/2. , decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , -0.5, decimal = 5) - #### - - cp2 = cp0.power(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4) - cp2 = cp0**cp1 - assert (cp2.get_item(0,0).as_array()[0][0][0] == 0.) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4) - - cp2 = cp0 ** 2 - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0] , 0., decimal=5) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0] , 1., decimal = 5) - - cp2 = cp0.maximum(cp1) - assert (cp2.get_item(0,0).as_array()[0][0][0] == cp1.get_item(0,0).as_array()[0][0][0]) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], cp2.get_item(1,0).as_array()[0][0][0], decimal=4) - - - cp2 = cp0.abs() - numpy.testing.assert_almost_equal(cp2.get_item(0,0).as_array()[0][0][0], 0., decimal=4) - numpy.testing.assert_almost_equal(cp2.get_item(1,0).as_array()[0][0][0], 1., decimal=4) - - cp2 = cp0.subtract(cp1) - s = cp2.sign() - numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], -1., decimal=4) - numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], -1., decimal=4) - - cp2 = cp0.add(cp1) - s = cp2.sqrt() - numpy.testing.assert_almost_equal(s.get_item(0,0).as_array()[0][0][0], numpy.sqrt(2), decimal=4) - numpy.testing.assert_almost_equal(s.get_item(1,0).as_array()[0][0][0], numpy.sqrt(4), decimal=4) - - s = cp0.sum() - numpy.testing.assert_almost_equal(s[0], 0, decimal=4) - s0 = 1 - s1 = 1 - for i in cp0.get_item(0,0).shape: - s0 *= i - for i in cp0.get_item(1,0).shape: - s1 *= i - - numpy.testing.assert_almost_equal(s[1], cp0.get_item(0,0).as_array()[0][0][0]*s0 +cp0.get_item(1,0).as_array()[0][0][0]*s1, decimal=4) # Set up phantom size N x N x vert by creating ImageGeometry, initialising the # ImageData object with this geometry and empty array and finally put some |