From df538e0d55513274ba27285764783359ef7b5c6c Mon Sep 17 00:00:00 2001 From: epapoutsellis Date: Wed, 17 Apr 2019 17:53:31 +0100 Subject: fix with cvxpy tomo2D recon --- .../Python/ccpi/optimisation/algorithms/PDHG.py | 2 +- .../ccpi/optimisation/functions/L2NormSquared.py | 12 +-- Wrappers/Python/wip/pdhg_TV_denoising.py | 17 +++- .../Python/wip/pdhg_TV_denoising_salt_pepper.py | 9 +- Wrappers/Python/wip/pdhg_TV_tomography2D.py | 102 +++++++++++++++++---- Wrappers/Python/wip/pdhg_tv_denoising_poisson.py | 8 ++ 6 files changed, 116 insertions(+), 34 deletions(-) (limited to 'Wrappers/Python') diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py index c0c2c10..6360ac1 100644 --- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py +++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py @@ -182,7 +182,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs): f.proximal_conjugate(y_tmp, sigma, out=y) operator.adjoint(y, out = x_tmp) - x_tmp *= -tau + x_tmp *= -1*tau x_tmp += x_old g.proximal(x_tmp, tau, out = x) diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py index 8a16c28..ca6e7a7 100644 --- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py +++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py @@ -94,12 +94,12 @@ class L2NormSquared(Function): return x/(1+2*tau) else: - tmp = x.subtract(self.b) - #tmp -= self.b - tmp /= (1+2*tau) - tmp += self.b - return tmp -# return (x-self.b)/(1+2*tau) + self.b +# tmp = x.subtract(self.b) +# tmp -= self.b +# tmp /= (1+2*tau) +# tmp += self.b +# return tmp + return (x-self.b)/(1+2*tau) + self.b # if self.b is not None: # out=x diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py index 0e7f628..f3980c4 100755 --- a/Wrappers/Python/wip/pdhg_TV_denoising.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising.py @@ -45,7 +45,7 @@ plt.title('Noisy data') plt.show() # Regularisation Parameter -alpha = 2 +alpha = 0.5 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ") method = '1' @@ -83,8 +83,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) @@ -149,7 +149,8 @@ if cvx_not_installable: if 'MOSEK' in installed_solvers(): solver = MOSEK else: - solver = SCS + solver = SCS + obj = Minimize( regulariser + fidelity) prob = Problem(obj) result = prob.solve(verbose = True, solver = solver) @@ -172,6 +173,14 @@ if cvx_not_installable: 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 cec9770..364d879 100644 --- a/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py +++ b/Wrappers/Python/wip/pdhg_TV_denoising_salt_pepper.py @@ -54,13 +54,8 @@ if method == '0': op1 = Gradient(ig) op2 = Identity(ig, ag) - # Form Composite Operator operator = BlockOperator(op1, op2, shape=(2,1) ) - #### Create functions -# f = FunctionComposition_new(operator, mixed_L12Norm(alpha), \ -# L2NormSq(0.5, b = noisy_data) ) - f1 = alpha * MixedL21Norm() f2 = L1Norm(b = noisy_data) @@ -73,12 +68,12 @@ else: # No Composite # ########################################################################### operator = Gradient(ig) - f = alpha * FunctionOperatorComposition(operator, MixedL21Norm()) + f = alpha * MixedL21Norm() g = L1Norm(b = noisy_data) ########################################################################### #%% -diag_precon = True +diag_precon = False if diag_precon: diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py index d0ef097..1d4023b 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 = 75 +N = 50 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 = 150 -angles = np.linspace(0,np.pi,100) +detectors = 50 +angles = np.linspace(0,np.pi,50) ag = AcquisitionGeometry('parallel','2D',angles, detectors) Aop = AstraProjectorSimple(ig, ag, 'cpu') @@ -75,11 +75,6 @@ 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 @@ -109,19 +104,18 @@ if diag_precon: tau, sigma = tau_sigma_precond(operator) -else: +else: + normK = operator.norm() sigma = 1 tau = 1/(sigma*normK**2) # Compute operator Norm -normK = operator.norm() + # Primal & dual stepsizes -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) @@ -132,9 +126,7 @@ t3 = timer() res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) t4 = timer() # -print ("No memopt in {}s, memopt in {}s ".format(t2-t1, t4 -t3)) -#%% plt.figure(figsize=(15,15)) plt.subplot(3,1,1) plt.imshow(res.as_array()) @@ -160,3 +152,81 @@ diff = (res1 - res).abs().as_array().max() # 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 diff --git a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py index ec745a2..82384ef 100644 --- a/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py +++ b/Wrappers/Python/wip/pdhg_tv_denoising_poisson.py @@ -169,9 +169,17 @@ if cvx_not_installable: 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), u1.value[int(N/2),:], label = 'CVX') + plt.legend() + print('Primal Objective (CVX) {} '.format(obj.value)) print('Primal Objective (PDHG) {} '.format(primal[-1])) + + -- cgit v1.2.3