diff options
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py | 69 | ||||
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/algorithms/__init__.py | 1 | ||||
| -rw-r--r-- | Wrappers/Python/wip/test_pdhg_gap.py | 129 |
3 files changed, 199 insertions, 0 deletions
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index 043fe38..8600e07 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -13,6 +13,9 @@ import time from ccpi.optimisation.operators import BlockOperator from ccpi.framework import BlockDataContainer + +import matplotlib.pyplot as plt + class PDHG(Algorithm): '''Primal Dual Hybrid Gradient''' @@ -80,3 +83,69 @@ class PDHG(Algorithm): ]) + +def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): + + # algorithmic parameters + if opt is None: + opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \ + 'memopt': False} + + if sigma is None and tau is None: + raise ValueError('Need sigma*tau||K||^2<1') + + niter = opt['niter'] if 'niter' in opt.keys() else 1000 + tol = opt['tol'] if 'tol' in opt.keys() else 1e-4 + memopt = opt['memopt'] if 'memopt' in opt.keys() else False + show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False + stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False + + + x_old = operator.domain_geometry().allocate() + y_old = operator.range_geometry().allocate() + + + xbar = x_old + x_tmp = x_old + x = x_old + + y_tmp = y_old + y = y_tmp + + # relaxation parameter + theta = 1 + + t = time.time() + + objective = [] + + for i in range(niter): + + # Gradient descent, Dual problem solution + y_tmp = y_old + sigma * operator.direct(xbar) + y = f.proximal_conjugate(y_tmp, sigma) + + # Gradient ascent, Primal problem solution + x_tmp = x_old - tau * operator.adjoint(y) + x = g.proximal(x_tmp, tau) + + #Update + xbar = x + theta * (x - x_old) + + x_old = x + y_old = y + + if i%100==0: + + primal = f(operator.direct(x)) + g(x) + dual = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y))) + print( i, primal, dual) + + plt.imshow(x.as_array()) + plt.show() +# print(f(operator.direct(x)) + g(x), i) + + t_end = time.time() + + return x, t_end - t, objective + diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py index 443bc78..a28c0bf 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py @@ -28,3 +28,4 @@ from .GradientDescent import GradientDescent from .FISTA import FISTA from .FBPD import FBPD from .PDHG import PDHG +from .PDHG import PDHG_old diff --git a/Wrappers/Python/wip/test_pdhg_gap.py b/Wrappers/Python/wip/test_pdhg_gap.py new file mode 100644 index 0000000..b196e36 --- /dev/null +++ b/Wrappers/Python/wip/test_pdhg_gap.py @@ -0,0 +1,129 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +Created on Tue Apr 2 12:26:24 2019 + +@author: vaggelis +""" + + +from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData + +import numpy as np +import matplotlib.pyplot as plt + +from ccpi.optimisation.algorithms import PDHG, PDHG_old + +from ccpi.optimisation.operators import BlockOperator, Identity, Gradient +from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \ + MixedL21Norm, BlockFunction, ScaledFunction + +from ccpi.astra.ops import AstraProjectorSimple +from skimage.util import random_noise + + +#%%############################################################################### +# Create phantom for TV tomography + +#import os +#import tomophantom +#from tomophantom import TomoP2D +#from tomophantom.supp.qualitymetrics import QualityTools + +#model = 1 # select a model number from the library +#N = 150 # set dimension of the phantom +## one can specify an exact path to the parameters file +## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' +#path = os.path.dirname(tomophantom.__file__) +#path_library2D = os.path.join(path, "Phantom2DLibrary.dat") +##This will generate a N_size x N_size phantom (2D) +#phantom_2D = TomoP2D.Model(model, N, path_library2D) +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) +#data = ImageData(phantom_2D, geometry=ig) + +N = 150 +x = np.zeros((N,N)) +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)] = 1 + +data = ImageData(x) +ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + + +detectors = 150 +angles = np.linspace(0,np.pi,100) + +ag = AcquisitionGeometry('parallel','2D',angles, detectors) +Aop = AstraProjectorSimple(ig, ag, 'cpu') +sin = Aop.direct(data) + +plt.imshow(sin.as_array()) +plt.title('Sinogram') +plt.colorbar() +plt.show() + +# Add Gaussian noise to the sinogram data +np.random.seed(10) +n1 = np.random.random(sin.shape) + +noisy_data = sin + ImageData(5*n1) + +plt.imshow(noisy_data.as_array()) +plt.title('Noisy Sinogram') +plt.colorbar() +plt.show() + + +#%% Works only with Composite Operator Structure of PDHG + +#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) + +# Create operators +op1 = Gradient(ig) +op2 = Aop + +# Form Composite Operator +operator = BlockOperator(op1, op2, shape=(2,1) ) + +alpha = 50 +f = BlockFunction( alpha * MixedL21Norm(), \ + 0.5 * L2NormSquared(b = noisy_data) ) +g = ZeroFun() + +# Compute operator Norm +normK = operator.norm() + +## Primal & dual stepsizes + +sigma = 10 +tau = 1/(sigma*normK**2) + +pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) +pdhg.max_iteration = 2000 +pdhg.update_objective_interval = 100 + +#pdhg.run(5000) + +opt = {'niter':2000} + +res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) + +#%% +#sol = pdhg.get_output().as_array() +#fig = plt.figure() +#plt.subplot(1,2,1) +#plt.imshow(noisy_data.as_array()) +##plt.colorbar() +#plt.subplot(1,2,2) +#plt.imshow(sol) +##plt.colorbar() +#plt.show() +# +# +##%% +#plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth') +#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon') +#plt.legend() +#plt.show() + + |