fix second case for TGV

1 files changed, 67 insertions, 74 deletions

diff --git a/Wrappers/Python/wip/pdhg_TGV_denoising.py b/Wrappers/Python/wip/pdhg_TGV_denoising.py
index 54b7131..0b6e594 100644
--- a/Wrappers/Python/wip/pdhg_TGV_denoising.py
+++ b/Wrappers/Python/wip/pdhg_TGV_denoising.py
@@ -55,10 +55,10 @@ plt.imshow(noisy_data.as_array())
 plt.title('Noisy data')
 plt.show()
 
-alpha, beta = 0.2, 1
+alpha, beta = 0.1, 0.5
 
 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '0'
+method = '1'
 
 if method == '0':
 
@@ -66,11 +66,9 @@ if method == '0':
     op11 = Gradient(ig)
     op12 = Identity(op11.range_geometry())
     
-    op22 = SymmetrizedGradient(op11.domain_geometry())
-    
+    op22 = SymmetrizedGradient(op11.domain_geometry())    
     op21 = ZeroOperator(ig, op22.range_geometry())
-    
-    
+        
     op31 = Identity(ig, ag)
     op32 = ZeroOperator(op22.domain_geometry(), ag)
     
@@ -90,21 +88,16 @@ else:
     op11 = Gradient(ig)
     op12 = Identity(op11.range_geometry())
     op22 = SymmetrizedGradient(op11.domain_geometry())    
-    op21 = ZeroOperator(ig, op22.range_geometry())
+    op21 = ZeroOperator(ig, op22.range_geometry())    
     
-    operator = BlockOperator(op11, -1*op12, \
-                             op21, op22, \
-                             shape=(2,2) )  
+    operator = BlockOperator(op11, -1*op12, op21, op22, shape=(2,2) )      
     
     f1 = alpha * MixedL21Norm()
     f2 = beta * MixedL21Norm()     
-    f = BlockFunction(f1, f2)     
     
-    g =  0.5 * L2NormSquared(b = noisy_data)    
+    f = BlockFunction(f1, f2)         
+    g = BlockFunction(0.5 * L2NormSquared(b = noisy_data), ZeroFunction())
      
-
-
-
 ## Compute operator Norm
 normK = operator.norm()
 #
@@ -147,65 +140,65 @@ 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((N, N))
-#    w2 = Variable((N, N))
-#    
-#    # 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 = []
-#    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())    
-#    solver = MOSEK
-#
-#    obj =  Minimize( regulariser +  fidelity)
-#    prob = Problem(obj)
-#    result = prob.solve(verbose = True, solver = solver)
-#    
-#    diff_cvx = numpy.abs( res[0].as_array() - u.value )   
-#    
-#    # Show result
-#    plt.figure(figsize=(15,15))
-#    plt.subplot(3,1,1)
-#    plt.imshow(res[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,N,N), res[0].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])) 
-#    print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max()))
+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((N, N))
+    w2 = Variable((N, N))
+    
+    # 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 = []
+    fidelity = 0.5 * sum_squares(u - noisy_data.as_array())    
+    solver = MOSEK
+
+    obj =  Minimize( regulariser +  fidelity)
+    prob = Problem(obj)
+    result = prob.solve(verbose = True, solver = solver)
+    
+    diff_cvx = numpy.abs( res[0].as_array() - u.value )   
+    
+    # Show result
+    plt.figure(figsize=(15,15))
+    plt.subplot(3,1,1)
+    plt.imshow(res[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,N,N), res[0].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])) 
+    print('Min/Max of absolute difference {}/{}'.format(diff_cvx.min(), diff_cvx.max()))