add cvx script

1 files changed, 101 insertions, 102 deletions

diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
index 75286e5..87d5328 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
@@ -110,7 +110,6 @@ else:
     tau = 1/(sigma*normK**2)
 
 
-
 # Setup and run the PDHG algorithm
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
 pdhg.max_iteration = 2000
@@ -143,104 +142,104 @@ plt.show()
 
 #%% 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:
-#    
-#
-#    ##Construct problem    
-#    u = Variable(N*N)
-#    #q = Variable()
-#    
-#    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 = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) 
-#    #constraints = [q>= fidelity, u>=0]
-#    constraints = [u>=0]
-#        
-#    solver = SCS
-#    obj =  Minimize( regulariser +  fidelity)
-#    prob = Problem(obj, constraints)
-#    result = prob.solve(verbose = True, solver = solver)    
-#     
-#
-###%% 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 = pnorm( u - noisy_data.as_array(),1)
-#    
-#    # 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)
-#    
-#           
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(pdhg.get_output().as_array())
-#    plt.title('PDHG solution')
-#    plt.colorbar()
-#    plt.subplot(3,1,2)
-#    plt.imshow(np.reshape(u.value, (N, N)))
-#    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), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
-#    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
-#    plt.legend()
-#    plt.title('Middle Line Profiles')
-#    plt.show()
-#            
-#    print('Primal Objective (CVX) {} '.format(obj.value))
-#    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ 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:
+    
+
+    ##Construct problem    
+    u = Variable(N*N)
+    #q = Variable()
+    
+    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 = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) 
+    #constraints = [q>= fidelity, u>=0]
+    constraints = [u>=0]
+        
+    solver = SCS
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj, constraints)
+    result = prob.solve(verbose = True, solver = solver)    
+     
+
+##%% 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 = pnorm( u - noisy_data.as_array(),1)
+    
+    # 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)
+    
+           
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(pdhg.get_output().as_array())
+    plt.title('PDHG solution')
+    plt.colorbar()
+    plt.subplot(3,1,2)
+    plt.imshow(np.reshape(u.value, (N, N)))
+    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), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
+    plt.legend()
+    plt.title('Middle Line Profiles')
+    plt.show()
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file