diff options
Diffstat (limited to 'Wrappers')
| -rwxr-xr-x | Wrappers/Python/wip/pdhg_TV_denoising.py | 162 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py | 167 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_TV_tomography2D.py | 162 | ||||
| -rw-r--r-- | Wrappers/Python/wip/pdhg_tv_denoising_poisson.py | 8 |
4 files changed, 304 insertions, 195 deletions
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py index f3980c4..24fe216 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -26,7 +26,7 @@ from timeit import default_timer as timer # Create phantom for TV Gaussian denoising -N = 100 +N = 500 data = np.zeros((N,N)) data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5 @@ -48,7 +48,7 @@ plt.show() alpha = 0.5 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '1' +method = '0' if method == '0': @@ -73,29 +73,43 @@ else: # No Composite # ########################################################################### operator = Gradient(ig) - f = alpha * MixedL21Norm() - g = 0.5 * L2NormSquared(b = noisy_data) + f = alpha * MixedL21Norm() + g = 0.5 * L2NormSquared(b = noisy_data) -# Compute operator Norm -normK = operator.norm() + +diag_precon = False + +if diag_precon: + + def tau_sigma_precond(operator): + + tau = 1/operator.sum_abs_row() + sigma = 1/ operator.sum_abs_col() + + return tau, sigma -# Primal & dual stepsizes -sigma = 1 -tau = 1/(sigma*normK**2) + tau, sigma = tau_sigma_precond(operator) + +else: + # Compute operator Norm + normK = operator.norm() + + # Primal & dual stepsizes + sigma = 1 + tau = 1/(sigma*normK**2) -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} + +opt = {'niter':2000} +opt1 = {'niter':2000, 'memopt': True} t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) t2 = timer() - t3 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) t4 = timer() -# -##%% + plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(res.as_array()) @@ -124,66 +138,66 @@ print(" Max of abs difference is {}".format(diff)) #%% Check with CVX solution -from ccpi.optimisation.operators import SparseFiniteDiff - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - ##Construct problem - u = Variable(ig.shape) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - # Define Total Variation as a regulariser - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - - diff_cvx = numpy.abs( res.as_array() - u.value ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u.value) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.legend() - - - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1])) +#from ccpi.optimisation.operators import SparseFiniteDiff +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +# +#if cvx_not_installable: +# +# ##Construct problem +# u = Variable(ig.shape) +# +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# # Define Total Variation as a regulariser +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +# fidelity = 0.5 * sum_squares(u - noisy_data.as_array()) +# +# # choose solver +# if 'MOSEK' in installed_solvers(): +# solver = MOSEK +# else: +# solver = SCS +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +# diff_cvx = numpy.abs( res.as_array() - u.value ) +# +# # Show result +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(res.as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# +# plt.subplot(3,1,2) +# plt.imshow(u.value) +# plt.title('CVX solution') +# plt.colorbar() +# +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.legend() +# +# +# +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(primal[-1])) diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py index 364d879..bd330fc 100644 --- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py @@ -8,18 +8,20 @@ Created on Fri Feb 22 14:53:03 2019 from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer -import numpy as np +import numpy as np +import numpy 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, L1Norm, \ - MixedL21Norm, FunctionOperatorComposition, BlockFunction - +from ccpi.optimisation.functions import ZeroFunction, L1Norm, \ + MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction from skimage.util import random_noise +from timeit import default_timer as timer + # ############################################################################ @@ -60,7 +62,7 @@ if method == '0': f2 = L1Norm(b = noisy_data) f = BlockFunction(f1, f2 ) - g = ZeroFun() + g = ZeroFunction() else: @@ -89,20 +91,48 @@ if diag_precon: else: # Compute operator Norm normK = operator.norm() - print ("normK", normK) + # Primal & dual stepsizes - sigma = 1/normK - tau = 1/normK -# tau = 1/(sigma*normK**2) + sigma = 1 + tau = 1/(sigma*normK**2) + -opt = {'niter':2000} +opt = {'niter':5000} +opt1 = {'niter':5000, 'memopt': True} +t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) - -plt.figure(figsize=(5,5)) +t2 = timer() + + +t3 = timer() +res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) +t4 = timer() + +plt.figure(figsize=(15,15)) +plt.subplot(3,1,1) plt.imshow(res.as_array()) +plt.title('no memopt') plt.colorbar() +plt.subplot(3,1,2) +plt.imshow(res1.as_array()) +plt.title('memopt') +plt.colorbar() +plt.subplot(3,1,3) +plt.imshow((res1 - res).abs().as_array()) +plt.title('diff') +plt.colorbar() +plt.show() +# +plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'memopt') +plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'no memopt') +plt.legend() plt.show() +# +print ("Time: No memopt in {}s, \n Time: Memopt in {}s ".format(t2-t1, t4 -t3)) +diff = (res1 - res).abs().as_array().max() +# +print(" Max of abs difference is {}".format(diff)) #pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma) #pdhg.max_iteration = 2000 @@ -132,43 +162,106 @@ plt.show() #plt.show() -#%% Compare with cvx - -try_cvx = input("Do you want CVX comparison (0/1)") +#%% Check with CVX solution -if try_cvx=='0': +from ccpi.optimisation.operators import SparseFiniteDiff +try: from cvxpy import * - import sys - sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts') - from cvx_functions import TV_cvx + cvx_not_installable = True +except ImportError: + cvx_not_installable = False - u = Variable((N, N)) + +if cvx_not_installable: + + u = Variable(ig.shape) + + DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') + DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') + + # Define Total Variation as a regulariser + regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) + fidelity = pnorm( u - noisy_data.as_array(),1) - regulariser = alpha * TV_cvx(u) - solver = MOSEK + + # choose solver + if 'MOSEK' in installed_solvers(): + solver = MOSEK + else: + solver = SCS + obj = Minimize( regulariser + fidelity) - constraints = [] - prob = Problem(obj, constraints) - - # Choose solver (SCS is fast but less accurate than MOSEK) + prob = Problem(obj) result = prob.solve(verbose = True, solver = solver) - - print('Objective value is {} '.format(obj.value)) - - diff_pdhg_cvx = np.abs(u.value - res.as_array()) - plt.imshow(diff_pdhg_cvx) + + diff_cvx = numpy.abs( res.as_array() - u.value ) + +# Show result + plt.figure(figsize=(15,15)) + plt.subplot(3,1,1) + plt.imshow(res.as_array()) + plt.title('PDHG solution') + plt.colorbar() + + plt.subplot(3,1,2) + plt.imshow(u.value) + plt.title('CVX solution') + plt.colorbar() + + plt.subplot(3,1,3) + plt.imshow(diff_cvx) + plt.title('Difference') plt.colorbar() - plt.title('|CVX-PDHG|') plt.show() - + + plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') - plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG') plt.legend() - plt.show() + -else: - print('No CVX solution available') + + + print('Primal Objective (CVX) {} '.format(obj.value)) + print('Primal Objective (PDHG) {} '.format(primal[-1])) + + + +#try_cvx = input("Do you want CVX comparison (0/1)") +# +#if try_cvx=='0': +# +# from cvxpy import * +# import sys +# sys.path.insert(0,'/Users/evangelos/Desktop/Projects/CCPi/CCPi-Framework/Wrappers/Python/ccpi/optimisation/cvx_scripts') +# from cvx_functions import TV_cvx +# +# u = Variable((N, N)) +# fidelity = pnorm( u - noisy_data.as_array(),1) +# regulariser = alpha * TV_cvx(u) +# solver = MOSEK +# obj = Minimize( regulariser + fidelity) +# constraints = [] +# prob = Problem(obj, constraints) +# +# # Choose solver (SCS is fast but less accurate than MOSEK) +# result = prob.solve(verbose = True, solver = solver) +# +# print('Objective value is {} '.format(obj.value)) +# +# diff_pdhg_cvx = np.abs(u.value - res.as_array()) +# plt.imshow(diff_pdhg_cvx) +# plt.colorbar() +# plt.title('|CVX-PDHG|') +# plt.show() +# +# plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX') +# plt.plot(np.linspace(0,N,N), res.as_array()[int(N/2),:], label = 'PDHG') +# plt.legend() +# plt.show() +# +#else: +# print('No CVX solution available') diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index 1d4023b..2713cfd 100644 --- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py +++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py @@ -43,7 +43,7 @@ from timeit import default_timer as timer #ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) #data = ImageData(phantom_2D, geometry=ig) -N = 50 +N = 75 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 @@ -52,8 +52,8 @@ data = ImageData(x) ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N) -detectors = 50 -angles = np.linspace(0,np.pi,50) +detectors = N +angles = np.linspace(0,np.pi,N) ag = AcquisitionGeometry('parallel','2D',angles, detectors) Aop = AstraProjectorSimple(ig, ag, 'cpu') @@ -82,7 +82,7 @@ op2 = Aop # Form Composite Operator operator = BlockOperator(op1, op2, shape=(2,1) ) -alpha = 50 +alpha = 10 f = BlockFunction( alpha * MixedL21Norm(), \ 0.5 * L2NormSquared(b = noisy_data) ) g = ZeroFunction() @@ -114,8 +114,8 @@ else: # Primal & dual stepsizes -opt = {'niter':5000} -opt1 = {'niter':5000, 'memopt': True} +opt = {'niter':2000} +opt1 = {'niter':2000, 'memopt': True} t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) @@ -155,78 +155,78 @@ print(" Max of abs difference is {}".format(diff)) #%% Check with CVX solution -from ccpi.optimisation.operators import SparseFiniteDiff -import astra -import numpy - -try: - from cvxpy import * - cvx_not_installable = True -except ImportError: - cvx_not_installable = False - - -if cvx_not_installable: - - - u = Variable( N*N) - - DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') - DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') - - regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) - - # create matrix representation for Astra operator - - vol_geom = astra.create_vol_geom(N, N) - proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) - - proj_id = astra.create_projector('strip', proj_geom, vol_geom) - - matrix_id = astra.projector.matrix(proj_id) - - ProjMat = astra.matrix.get(matrix_id) - - fidelity = 0.5 * sum_squares(ProjMat * u - noisy_data.as_array().ravel()) - - # choose solver - if 'MOSEK' in installed_solvers(): - solver = MOSEK - else: - solver = SCS - - obj = Minimize( regulariser + fidelity) - prob = Problem(obj) - result = prob.solve(verbose = True, solver = solver) - -#%% - - u_rs = np.reshape(u.value, (N,N)) - - diff_cvx = numpy.abs( res.as_array() - u_rs ) - - # Show result - plt.figure(figsize=(15,15)) - plt.subplot(3,1,1) - plt.imshow(res.as_array()) - plt.title('PDHG solution') - plt.colorbar() - - plt.subplot(3,1,2) - plt.imshow(u_rs) - plt.title('CVX solution') - plt.colorbar() - - plt.subplot(3,1,3) - plt.imshow(diff_cvx) - plt.title('Difference') - plt.colorbar() - plt.show() - - plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') - plt.plot(np.linspace(0,N,N), u_rs[int(N/2),:], label = 'CVX') - plt.legend() - - - print('Primal Objective (CVX) {} '.format(obj.value)) - print('Primal Objective (PDHG) {} '.format(primal[-1]))
\ No newline at end of file +#from ccpi.optimisation.operators import SparseFiniteDiff +#import astra +#import numpy +# +#try: +# from cvxpy import * +# cvx_not_installable = True +#except ImportError: +# cvx_not_installable = False +# +# +#if cvx_not_installable: +# +# +# u = Variable( N*N) +# +# DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann') +# DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann') +# +# regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0)) +# +# # create matrix representation for Astra operator +# +# vol_geom = astra.create_vol_geom(N, N) +# proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles) +# +# proj_id = astra.create_projector('strip', proj_geom, vol_geom) +# +# matrix_id = astra.projector.matrix(proj_id) +# +# ProjMat = astra.matrix.get(matrix_id) +# +# fidelity = 0.5 * sum_squares(ProjMat * u - noisy_data.as_array().ravel()) +# +# # choose solver +# if 'MOSEK' in installed_solvers(): +# solver = MOSEK +# else: +# solver = SCS +# +# obj = Minimize( regulariser + fidelity) +# prob = Problem(obj) +# result = prob.solve(verbose = True, solver = solver) +# +##%% +# +# u_rs = np.reshape(u.value, (N,N)) +# +# diff_cvx = numpy.abs( res.as_array() - u_rs ) +# +# # Show result +# plt.figure(figsize=(15,15)) +# plt.subplot(3,1,1) +# plt.imshow(res.as_array()) +# plt.title('PDHG solution') +# plt.colorbar() +# +# plt.subplot(3,1,2) +# plt.imshow(u_rs) +# plt.title('CVX solution') +# plt.colorbar() +# +# plt.subplot(3,1,3) +# plt.imshow(diff_cvx) +# plt.title('Difference') +# plt.colorbar() +# plt.show() +# +# plt.plot(np.linspace(0,N,N), res1.as_array()[int(N/2),:], label = 'PDHG') +# plt.plot(np.linspace(0,N,N), u_rs[int(N/2),:], label = 'CVX') +# plt.legend() +# +# +# print('Primal Objective (CVX) {} '.format(obj.value)) +# print('Primal Objective (PDHG) {} '.format(primal[-1]))
\ No newline at end of file diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py index 82384ef..2c7e145 100644 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py @@ -48,7 +48,7 @@ alpha = 2 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") -method = '0' +method = '1' if method == '0': # Create operators @@ -81,8 +81,8 @@ normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) -opt = {'niter':2000} -opt1 = {'niter':2000, 'memopt': True} +opt = {'niter':5000} +opt1 = {'niter':5000, 'memopt': True} t1 = timer() res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) @@ -92,6 +92,8 @@ t3 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) t4 = timer() +print(pdgap[-1]) + plt.figure(figsize=(15,15)) plt.subplot(3,1,1) |