diff options
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/funcs.py | 2 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py | 220 |
2 files changed, 4 insertions, 218 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/funcs.py b/Wrappers/Python/ccpi/optimisation/funcs.py index 9b9fc36..99af275 100755 --- a/Wrappers/Python/ccpi/optimisation/funcs.py +++ b/Wrappers/Python/ccpi/optimisation/funcs.py @@ -154,6 +154,8 @@ class Norm2sq(Function): self.L = 2.0*self.c*(self.A.get_max_sing_val()**2) except AttributeError as ae: pass + except NotImplementedError as noe: + pass def grad(self,x): #return 2*self.c*self.A.adjoint( self.A.direct(x) - self.b ) diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py index f102f1e..8298c03 100755 --- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py +++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py @@ -98,9 +98,9 @@ class BlockOperator(Operator): 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)) + prod = self.get_item(row, col).adjoint(x.get_item(col)) else: - prod += self.get_item(col,row).adjoint(x.get_item(row)) + prod += self.get_item(row, col).adjoint(x.get_item(col)) res.append(prod) return BlockDataContainer(*res, shape=shape) @@ -137,221 +137,5 @@ class BlockOperator(Operator): shape = (self.shape[1], self.shape[0]) return type(self)(*self.operators, shape=shape) -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) - - - - - - - - - if __name__ == '__main__': pass - #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 -# from ccpi.optimisation.funcs import Norm2sq, Norm1 -# from ccpi.framework import ImageGeometry, AcquisitionGeometry -# from ccpi.optimisation.Algorithms import GradientDescent -# #from ccpi.optimisation.Algorithms import CGLS -# import matplotlib.pyplot as plt - - -# # 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 -# # data into its array, and display one slice as image. - -# # Image parameters -# N = 128 -# vert = 4 - -# # Set up image geometry -# ig = ImageGeometry(voxel_num_x=N, -# voxel_num_y=N, -# voxel_num_z=vert) - -# # Set up empty image data -# Phantom = ImageData(geometry=ig, -# dimension_labels=['horizontal_x', -# 'horizontal_y', -# 'vertical']) -# Phantom += 0.05 -# # Populate image data by looping over and filling slices -# i = 0 -# while i < vert: -# if vert > 1: -# x = Phantom.subset(vertical=i).array -# else: -# x = Phantom.array -# x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 -# x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.94 -# if vert > 1 : -# Phantom.fill(x, vertical=i) -# i += 1 - - -# perc = 0.02 -# # Set up empty image data -# noise = ImageData(numpy.random.normal(loc = 0.04 , -# scale = perc , -# size = Phantom.shape), geometry=ig, -# dimension_labels=['horizontal_x', -# 'horizontal_y', -# 'vertical']) -# Phantom += noise - -# # Set up AcquisitionGeometry object to hold the parameters of the measurement -# # setup geometry: # Number of angles, the actual angles from 0 to -# # pi for parallel beam, set the width of a detector -# # pixel relative to an object pixe and the number of detector pixels. -# angles_num = 20 -# det_w = 1.0 -# det_num = N - -# angles = numpy.linspace(0,numpy.pi,angles_num,endpoint=False,dtype=numpy.float32)*\ -# 180/numpy.pi - -# # Inputs: Geometry, 2D or 3D, angles, horz detector pixel count, -# # horz detector pixel size, vert detector pixel count, -# # vert detector pixel size. -# ag = AcquisitionGeometry('parallel', -# '3D', -# angles, -# N, -# det_w, -# vert, -# det_w) - -# # Set up Operator object combining the ImageGeometry and AcquisitionGeometry -# # wrapping calls to CCPi projector. -# A = CCPiProjectorSimple(ig, ag) - -# # Forward and backprojection are available as methods direct and adjoint. Here -# # generate test data b and some noise - -# b = A.direct(Phantom) - - -# #z = A.adjoint(b) - - -# # Using the test data b, different reconstruction methods can now be set up as -# # demonstrated in the rest of this file. In general all methods need an initial -# # guess and some algorithm options to be set. Note that 100 iterations for -# # some of the methods is a very low number and 1000 or 10000 iterations may be -# # needed if one wants to obtain a converged solution. -# x_init = ImageData(geometry=ig, -# dimension_labels=['horizontal_x','horizontal_y','vertical']) -# X_init = BlockDataContainer(x_init) -# B = BlockDataContainer(b, -# ImageData(geometry=ig, dimension_labels=['horizontal_x','horizontal_y','vertical'])) - -# # setup a tomo identity -# Ibig = 1e5 * TomoIdentity(geometry=ig) -# Ismall = 1e-5 * TomoIdentity(geometry=ig) - -# # composite operator -# Kbig = BlockOperator(A, Ibig, shape=(2,1)) -# Ksmall = BlockOperator(A, Ismall, shape=(2,1)) - -# #out = K.direct(X_init) - -# f = Norm2sq(Kbig,B) -# f.L = 0.00003 - -# fsmall = Norm2sq(Ksmall,B) -# f.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 - -# cg = CGLS() -# cg.set_up(X_init, Kbig, B ) -# cg.max_iteration = 1 - -# cgsmall = CGLS() -# cgsmall.set_up(X_init, Ksmall, B ) -# cgsmall.max_iteration = 1 - - -# cgs = CGLS() -# cgs.set_up(x_init, A, b ) -# cgs.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())) - -# cg.run(10, lambda it,val: print ("iteration {} objective {}".format(it,val))) - -# 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))) -# # for _ in cg: -# # print ("iteration {} {}".format(cg.iteration, cg.get_current_loss())) -# # -# # fig = plt.figure() -# # plt.imshow(cg.get_output().get_item(0,0).subset(vertical=0).as_array()) -# # plt.title('Composite CGLS') -# # plt.show() -# # -# # for _ in cgs: -# # print ("iteration {} {}".format(cgs.iteration, cgs.get_current_loss())) -# # -# fig = plt.figure() -# plt.subplot(1,5,1) -# plt.imshow(Phantom.subset(vertical=0).as_array()) -# plt.title('Simulated Phantom') -# plt.subplot(1,5,2) -# plt.imshow(gd.get_output().subset(vertical=0).as_array()) -# plt.title('Simple Gradient Descent') -# plt.subplot(1,5,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.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.title('Composite CGLS\nsmall lambda') -# plt.show() |