pdhg examples

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)