From 4290663849f7f396fa3a6f8b4de8a312372e2556 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Fri, 8 Mar 2019 07:47:20 -0500 Subject: working implementation of CGLS Tikhonov example CGLS with Tikhonov regularisation with Block structures --- .../Python/ccpi/framework/BlockDataContainer.py | 19 +++------ .../ccpi/optimisation/operators/BlockOperator.py | 39 +++++++++++-------- Wrappers/Python/wip/CGLS_tikhonov.py | 45 ++++++++++++++-------- 3 files changed, 58 insertions(+), 45 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 8152bff..d509d25 100755 --- a/Wrappers/Python/ccpi/framework/BlockDataContainer.py +++ b/Wrappers/Python/ccpi/framework/BlockDataContainer.py @@ -19,14 +19,8 @@ class BlockDataContainer(object): '''Class to hold DataContainers as column vector''' __array_priority__ = 1 def __init__(self, *args, **kwargs): - '''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: @@ -38,7 +32,7 @@ class BlockDataContainer(object): if len(args) != n_elements: raise ValueError( 'Dimension and size do not match: expected {} got {}' - .format(n_elements,len(args))) + .format(n_elements, len(args))) def __iter__(self): @@ -60,7 +54,6 @@ class BlockDataContainer(object): if isinstance(other, Number): return True elif isinstance(other, list): - # TODO look elements should be numbers for ot in other: if not isinstance(ot, (Number,\ numpy.int, numpy.int8, numpy.int16, numpy.int32, numpy.int64,\ @@ -72,10 +65,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): 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) @@ -308,8 +303,4 @@ class BlockDataContainer(object): def __itruediv__(self, other): '''Inline truedivision''' return self.__idiv__(other) - #@property - #def T(self): - # '''return the transposed of self''' - # shape = (self.shape[1], self.shape[0]) - # return type(self)(*self.containers, shape=shape) + diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index c9bf794..f102f1e 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -81,6 +81,28 @@ class BlockOperator(Operator): prod += self.get_item(row,col).direct(x.get_item(col)) res.append(prod) return BlockDataContainer(*res, shape=shape) + + def adjoint(self, x, out=None): + '''Adjoint operation for the BlockOperator + + BlockOperator may contain both LinearOperator and Operator + This method exists in BlockOperator as it is not known what type of + Operator it will contain. + + Raises: ValueError if the contained Operators are not linear + ''' + if not functools.reduce(lambda x,y: x and y, self.operators.is_linear(), True): + raise ValueError('Not all operators in Block are linear.') + shape = self.get_output_shape(x.shape, adjoint=True) + res = [] + for row in range(self.shape[1]): + for col in range(self.shape[0]): + if col == 0: + prod = self.get_item(col,row).adjoint(x.get_item(row)) + else: + prod += self.get_item(col,row).adjoint(x.get_item(row)) + res.append(prod) + return BlockDataContainer(*res, shape=shape) def get_output_shape(self, xshape, adjoint=False): sshape = self.shape[1] @@ -142,22 +164,7 @@ class BlockLinearOperator(BlockOperator): 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 - - 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) + diff --git a/Wrappers/Python/wip/CGLS_tikhonov.py b/Wrappers/Python/wip/CGLS_tikhonov.py index 7178510..f247896 100644 --- a/Wrappers/Python/wip/CGLS_tikhonov.py +++ b/Wrappers/Python/wip/CGLS_tikhonov.py @@ -105,46 +105,57 @@ B = BlockDataContainer(b, # setup a tomo identity Ibig = 1e5 * TomoIdentity(geometry=ig) Ismall = 1e-5 * TomoIdentity(geometry=ig) +Iok = 1e1 * TomoIdentity(geometry=ig) # composite operator Kbig = BlockOperator(A, Ibig, shape=(2,1)) Ksmall = BlockOperator(A, Ismall, shape=(2,1)) - +Kok = BlockOperator(A, Iok, shape=(2,1)) + #out = K.direct(X_init) f = Norm2sq(Kbig,B) f.L = 0.00003 fsmall = Norm2sq(Ksmall,B) -f.L = 0.00003 - +fsmall.L = 0.00003 + +fok = Norm2sq(Kok,B) +fok.L = 0.00003 + simplef = Norm2sq(A, b) simplef.L = 0.00003 gd = GradientDescent( x_init=x_init, objective_function=simplef, rate=simplef.L) gd.max_iteration = 10 - + +Kbig.direct(X_init) +Kbig.adjoint(B) cg = CGLS() cg.set_up(X_init, Kbig, B ) -cg.max_iteration = 1 +cg.max_iteration = 5 cgsmall = CGLS() cgsmall.set_up(X_init, Ksmall, B ) -cgsmall.max_iteration = 1 +cgsmall.max_iteration = 5 cgs = CGLS() cgs.set_up(x_init, A, b ) cgs.max_iteration = 6 -# # + +cgok = CGLS() +cgok.set_up(X_init, Kok, B ) +cgok.max_iteration = 6 +# # #out.__isub__(B) #out2 = K.adjoint(out) #(2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b ) for _ in gd: - print ("iteration {} {}".format(gd.iteration, gd.get_current_loss())) + print ("iteration {} {}".format(gd.iteration, gd.get_last_loss())) cg.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) @@ -152,6 +163,7 @@ cgs.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) +cgok.run(10, verbose=True) # # for _ in cg: # print ("iteration {} {}".format(cg.iteration, cg.get_current_loss())) # # @@ -164,19 +176,22 @@ cgsmall.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val) # print ("iteration {} {}".format(cgs.iteration, cgs.get_current_loss())) # # fig = plt.figure() -plt.subplot(1,5,1) +plt.subplot(1,6,1) plt.imshow(Phantom.subset(vertical=0).as_array()) plt.title('Simulated Phantom') -plt.subplot(1,5,2) +plt.subplot(1,6,2) plt.imshow(gd.get_output().subset(vertical=0).as_array()) plt.title('Simple Gradient Descent') -plt.subplot(1,5,3) +plt.subplot(1,6,3) plt.imshow(cgs.get_output().subset(vertical=0).as_array()) plt.title('Simple CGLS') -plt.subplot(1,5,4) -plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) +plt.subplot(1,6,4) +plt.imshow(cg.get_output().get_item(0).subset(vertical=0).as_array()) plt.title('Composite CGLS\nbig lambda') -plt.subplot(1,5,5) -plt.imshow(cgsmall.get_output().get_item(0,0).subset(vertical=0).as_array()) +plt.subplot(1,6,5) +plt.imshow(cgsmall.get_output().get_item(0).subset(vertical=0).as_array()) plt.title('Composite CGLS\nsmall lambda') +plt.subplot(1,6,6) +plt.imshow(cgok.get_output().get_item(0).subset(vertical=0).as_array()) +plt.title('Composite CGLS\nok lambda') plt.show() -- cgit v1.2.3