fix with cvxpy tomo2D recon

6 files changed, 116 insertions, 34 deletions

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]))
+    
+