From 41b536a3f2a33d5e527d8b7ada63b47b1abbbca8 Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Thu, 11 Apr 2019 16:55:20 +0100 Subject: changes for memopt option --- .../Python/ccpi/framework/BlockDataContainer.py | 35 +++++++++++++++++-- Wrappers/Python/ccpi/framework/framework.py | 2 ++ .../Python/ccpi/optimisation/algorithms/PDHG.py | 22 ++++++------ .../ccpi/optimisation/functions/MixedL21Norm.py | 40 ++++++---------------- Wrappers/Python/wip/pdhg_TV_denoising.py | 13 ++++--- Wrappers/Python/wip/pdhg_TV_tomography2D.py | 31 ++++++++++++++--- 6 files changed, 93 insertions(+), 50 deletions(-) diff --git a/Wrappers/Python/ccpi/framework/BlockDataContainer.py b/Wrappers/Python/ccpi/framework/BlockDataContainer.py index 69b6931..8934f49 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 @@ -341,3 +341,32 @@ class BlockDataContainer(object): '''Inline truedivision''' return self.__idiv__(other) +if __name__ == '__main__': + + M, N, K = 2,3,5 + from ccpi.framework import ImageGeometry, BlockGeometry, BlockDataContainer + + ig = ImageGeometry(N, M) + u = ig.allocate('random_int') + + BG = BlockGeometry(ig, ig) + U = BG.allocate('random_int') + + U_nested = BlockDataContainer(BlockDataContainer(u, u), u) + + + res1 = U + u + res2 = U_nested + u +# res2 = u + U + + + + + + + + + + + + \ No newline at end of file diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py index 07c2ead..e03a29c 100755 --- a/Wrappers/Python/ccpi/framework/framework.py +++ b/Wrappers/Python/ccpi/framework/framework.py @@ -1109,6 +1109,7 @@ class AcquisitionData(DataContainer): class DataProcessor(object): + '''Defines a generic DataContainer processor accepts DataContainer as inputs and @@ -1154,6 +1155,7 @@ class DataProcessor(object): raise NotImplementedError('Implement basic checks for input DataContainer') def get_output(self, out=None): + for k,v in self.__dict__.items(): if v is None and k != 'output': raise ValueError('Key {0} is None'.format(k)) diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 2a69857..5323e76 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -161,18 +161,20 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): if not memopt: - y_tmp = y_old + sigma * operator.direct(xbar) - y = f.proximal_conjugate(y_tmp, sigma) + y_old += sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_old, sigma) - x_tmp = x_old - tau * operator.adjoint(y) - x = g.proximal(x_tmp, tau) - - xbar = x + theta * (x - x_old) + x_old -= tau*operator.adjoint(y) + x = g.proximal(x_old, tau) + + xbar.fill(x) + xbar -= x_old + xbar *= theta + xbar += x - x_old = x - y_old = y - - + x_old.fill(x) + y_old.fill(y) + else: operator.direct(xbar, out = y_tmp) diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py index 0ce7d8a..0c658a4 100755 --- a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py +++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py @@ -77,9 +77,7 @@ class MixedL21Norm(Function): pass def proximal_conjugate(self, x, tau, out=None): - - - + if self.SymTensor: param = [1]*x.shape[0] @@ -91,38 +89,22 @@ class MixedL21Norm(Function): return res else: -# pass - -# # 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 + tmp = [ el*el for el in x.containers] + res = (sum(tmp).sqrt()).maximum(1.0) + return x/res else: - - tmp = [ el*el for el in x] - res = (sum(tmp).sqrt()).maximum(1.0) - out.fill(x/res) - - - - # 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 ) + + res = (sum(x**2).sqrt()).maximum(1.0) + out.fill(x/res) + + + - def __rmul__(self, scalar): return ScaledFunction(self, scalar) diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py index 598acb0..22fee90 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -23,11 +23,13 @@ from timeit import default_timer as timer def dt(steps): return steps[-1] - steps[-2] +#%% + # ############################################################################ # Create phantom for TV denoising -N = 100 +N = 200 data = np.zeros((N,N)) data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -62,8 +64,8 @@ if method == '0': #### Create functions - f1 = MixedL21Norm() - f2 = L2NormSquared(b = noisy_data) + f1 = alpha * MixedL21Norm() + f2 = 0.5 * L2NormSquared(b = noisy_data) f = BlockFunction(f1, f2) g = ZeroFun() @@ -88,15 +90,18 @@ print ("normK", normK) sigma = 1 tau = 1/(sigma*normK**2) -<<<<<<< HEAD opt = {'niter':100} opt1 = {'niter':100, 'memopt': True} +t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +print(timer()-t1) print("with memopt \n") +t2 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +print(timer()-t2) plt.figure(figsize=(5,5)) plt.imshow(res.as_array()) diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index f06f166..159f2ea 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -21,6 +21,7 @@ from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \ from ccpi.astra.ops import AstraProjectorSimple from skimage.util import random_noise +from timeit import default_timer as timer #%%############################################################################### @@ -118,15 +119,37 @@ else: # #pdhg.run(5000) -opt = {'niter':2000} -# +opt = {'niter':300} +opt1 = {'niter':300, 'memopt': True} + + +t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) +print(timer()-t1) plt.figure(figsize=(5,5)) plt.imshow(res.as_array()) plt.colorbar() -plt.show() - +plt.show() + +#%% +print("with memopt \n") +# +t2 = timer() +res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +#print(timer()-t2) +# +# +plt.figure(figsize=(5,5)) +plt.imshow(res1.as_array()) +plt.colorbar() +plt.show() +# +#%% +plt.figure(figsize=(5,5)) +plt.imshow(np.abs(res1.as_array()-res.as_array())) +plt.colorbar() +plt.show() #%% #sol = pdhg.get_output().as_array() #fig = plt.figure() -- cgit v1.2.3