add also cvx demos

3 files changed, 211 insertions, 120 deletions

diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
index e9bad7e..8453d20 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
@@ -206,4 +206,77 @@ plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]
 plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'TGV reconstruction')
 plt.legend()
 plt.title('Middle Line Profiles')
-plt.show()
\ No newline at end of file
+plt.show()
+
+#%% 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:    
+    
+    u = Variable(ig.shape)
+    w1 = Variable(ig.shape)
+    w2 = Variable(ig.shape)
+    
+    # create TGV regulariser
+    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) - vec(w1), \
+                                           DY.matrix() * vec(u) - vec(w2)]), 2, axis = 0)) + \
+                  beta * sum(norm(vstack([ DX.matrix().transpose() * vec(w1), DY.matrix().transpose() * vec(w2), \
+                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ), \
+                                      0.5 * ( DX.matrix().transpose() * vec(w2) + DY.matrix().transpose() * vec(w1) ) ]), 2, axis = 0  ) )  
+    
+    constraints = []
+    
+    # choose solver
+    if 'MOSEK' in installed_solvers():
+        solver = MOSEK
+    else:
+        solver = SCS      
+    
+    # fidelity
+    if noise == 's&p':
+        fidelity = pnorm( u - noisy_data.as_array(),1)
+    elif noise == 'poisson':
+        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
+        solver = SCS
+    elif noise == 'gaussian':
+        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
+        
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = solver)
+    
+    diff_cvx = numpy.abs( pdhg.get_output()[0].as_array() - u.value )
+        
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(pdhg.get_output()[0].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,ig.shape[1],ig.shape[1]), pdhg.get_output()[0].as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
index c472f36..74e7901 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TV_Denoising.py
@@ -72,7 +72,7 @@ else:
 if len(sys.argv) > 1:
     which_noise = int(sys.argv[1])
 else:
-    which_noise = 1
+    which_noise = 0
 print ("Applying {} noise")
 
 if len(sys.argv) > 2:
@@ -83,7 +83,7 @@ print ("method ", method)
 
 
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SHAPES)
+data = loader.load(TestData.SHAPES, size=(50,50))
 ig = data.geometry
 ag = ig
 
@@ -211,61 +211,71 @@ else:
     plt.show()
 
 
-###%% 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)
-#    
-#    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
-#        
-#    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(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), 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]))
+##%% 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      
+    
+    # fidelity
+    if noise == 's&p':
+        fidelity = pnorm( u - noisy_data.as_array(),1)
+    elif noise == 'poisson':
+        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
+        solver = SCS
+    elif noise == 'gaussian':
+        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
+    
+            
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = solver)
+    
+    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
+        
+    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(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,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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]))
diff --git a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
index 8a9920c..e16c5a6 100644
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
@@ -185,65 +185,73 @@ plt.show()
 
 
 
-###%% 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_squares(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( pdhg.get_output().as_array() - u.value )
-#        
-#    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(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), 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]))
+##%% 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_squares(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+    
+    # choose solver
+    if 'MOSEK' in installed_solvers():
+        solver = MOSEK
+    else:
+        solver = SCS      
+    
+    # fidelity
+    if noise == 's&p':
+        fidelity = pnorm( u - noisy_data.as_array(),1)
+    elif noise == 'poisson':
+        fidelity = sum(kl_div(noisy_data.as_array(), u)) 
+        solver = SCS
+    elif noise == 'gaussian':
+        fidelity = 0.5 * sum_squares(noisy_data.as_array() - u)
+        
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = solver)
+    
+    diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
+        
+    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(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,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'PDHG')
+    plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), u.value[int(ig.shape[0]/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]))
 #
 #
 #