final demos

3 files changed, 221 insertions, 236 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 4d6da00..e9bad7e 100755
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_TGV_Denoising.py
@@ -23,23 +23,21 @@
 Total Generalised Variation (TGV) Denoising using PDHG algorithm:
 
 
-Problem:     min_{x} \alpha * ||\nabla x - w||_{2,1} +
-                     \beta *  || E w ||_{2,1} +
-                     fidelity
-
-            where fidelity can be as follows depending on the noise characteristics
-            of the data:
-                * Norm2Squared     \frac{1}{2} * || x - g ||_{2}^{2}
-                * KullbackLeibler  \int u - g * log(u) + Id_{u>0}
-                * L1Norm           ||u - g||_{1}
-
-             \alpha: Regularization parameter
-             \beta:  Regularization parameter
-             
+Problem:     min_{u} \alpha * ||\nabla u - w||_{2,1} +
+                     \beta *  || E u ||_{2,1} +
+                     Fidelity(u, g)
+                     
              \nabla: Gradient operator 
-              E: Symmetrized Gradient operator
+             E: Symmetrized Gradient operator             
+             \alpha: Regularization parameter
+             \beta:  Regularization parameter             
              
-             g: Noisy Data with Salt & Pepper Noise
+             g: Noisy Data 
+                          
+             Fidelity =  1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian
+                         2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper
+                         3) Kullback Leibler (\int u - g * log(u) + Id_{u>0})  if Noise is Poisson
+                                
                           
              Method = 0 ( PDHG - split ) :  K = [ \nabla, - Identity
                                                   ZeroOperator, E 
@@ -48,10 +46,15 @@ Problem:     min_{x} \alpha * ||\nabla x - w||_{2,1} +
                                                                     
              Method = 1 (PDHG - explicit ):  K = [ \nabla, - Identity
                                                   ZeroOperator, E ] 
+                          
+             Default: TGV denoising
+             noise = Gaussian
+             Fidelity = L2NormSquarred 
+             method = 0             
                                                                 
 """
 
-from ccpi.framework import ImageData, ImageGeometry
+from ccpi.framework import ImageData
 
 import numpy as np 
 import numpy                          
@@ -63,14 +66,14 @@ from ccpi.optimisation.operators import BlockOperator, Identity, \
                         Gradient, SymmetrizedGradient, ZeroOperator
 from ccpi.optimisation.functions import ZeroFunction, L1Norm, \
                       MixedL21Norm, BlockFunction, KullbackLeibler, L2NormSquared
-import sys
+                      
+from ccpi.framework import TestData
+import os, sys
 if int(numpy.version.version.split('.')[1]) > 12:
     from skimage.util import random_noise
 else:
     from demoutil import random_noise
 
-
-
 # user supplied input
 if len(sys.argv) > 1:
     which_noise = int(sys.argv[1])
@@ -81,35 +84,23 @@ print ("Applying {} noise")
 if len(sys.argv) > 2:
     method = sys.argv[2]
 else:
-    method = '1'
+    method = '0'
 print ("method ", method)
-# Create phantom for TGV SaltPepper denoising
-
-N = 100
-
-x1 = np.linspace(0, int(N/2), N)
-x2 = np.linspace(int(N/2), 0., N)
-xv, yv = np.meshgrid(x1, x2)
-
-xv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1] = yv[int(N/4):int(3*N/4)-1, int(N/4):int(3*N/4)-1].T
 
-#data = xv
 
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
-data = ig.allocate()
-data.fill(xv/xv.max())
+loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
+data = loader.load(TestData.SHAPES)
+ig = data.geometry
 ag = ig
 
 # Create noisy data. 
-# Apply Salt & Pepper noise
-# gaussian
-# poisson
 noises = ['gaussian', 'poisson', 's&p']
 noise = noises[which_noise]
 if noise == 's&p':
     n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2, seed=10)
 elif noise == 'poisson':
-    n1 = random_noise(data.as_array(), mode = noise, seed = 10)
+    scale = 5
+    n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale
 elif noise == 'gaussian':
     n1 = random_noise(data.as_array(), mode = noise, seed = 10)
 else:
@@ -128,23 +119,24 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-# Regularisation Parameters
+# Regularisation Parameter depending on the noise distribution
 if noise == 's&p':
     alpha = 0.8
 elif noise == 'poisson':
-    alpha = .1
-elif noise == 'gaussian':
     alpha = .3
+elif noise == 'gaussian':
+    alpha = .2
 
-beta = numpy.sqrt(2)* alpha
+# TODO add ref why this choice
+beta = 2 * alpha
 
-# fidelity
+# Fidelity
 if noise == 's&p':
     f3 = L1Norm(b=noisy_data)
 elif noise == 'poisson':
     f3 = KullbackLeibler(noisy_data)
 elif noise == 'gaussian':
-    f3 = L2NormSquared(b=noisy_data)
+    f3 = 0.5 * L2NormSquared(b=noisy_data)
 
 if method == '0':
     
@@ -162,7 +154,6 @@ if method == '0':
         
     f1 = alpha * MixedL21Norm()
     f2 = beta * MixedL21Norm() 
-    # f3 depends on the noise characteristics
     
     f = BlockFunction(f1, f2, f3)         
     g = ZeroFunction()
@@ -183,28 +174,27 @@ else:
     f = BlockFunction(f1, f2)         
     g = BlockFunction(f3, ZeroFunction())
      
-## Compute operator Norm
+# Compute operator Norm
 normK = operator.norm()
-#
+
 # Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
 
-
 # Setup and run the PDHG algorithm
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 200
-pdhg.run(2000, verbose = True)
+pdhg.update_objective_interval = 100
+pdhg.run(2000)
 
-#%%
+# Show results
 plt.figure(figsize=(20,5))
 plt.subplot(1,4,1)
-plt.imshow(data.as_array())
+plt.imshow(data.subset(channel=0).as_array())
 plt.title('Ground Truth')
 plt.colorbar()
 plt.subplot(1,4,2)
-plt.imshow(noisy_data.as_array())
+plt.imshow(noisy_data.subset(channel=0).as_array())
 plt.title('Noisy Data')
 plt.colorbar()
 plt.subplot(1,4,3)
@@ -212,10 +202,8 @@ plt.imshow(pdhg.get_output()[0].as_array())
 plt.title('TGV Reconstruction')
 plt.colorbar()
 plt.subplot(1,4,4)
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output()[0].as_array()[int(N/2),:], label = 'TGV reconstruction')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
+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()
-
-
+plt.show()
\ No newline at end of file
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 5aa334d..8a9920c 100644
--- a/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
+++ b/Wrappers/Python/demos/PDHG_examples/GatherAll/PDHG_Tikhonov_Denoising.py
@@ -18,26 +18,37 @@
 # limitations under the License.
 #
 #=========================================================================
+
 """ 
 
-Tikhonov Denoising using PDHG algorithm:
+``Tikhonov`` Regularization Denoising using PDHG algorithm:
 
 
-Problem:     min_{x} \alpha * ||\nabla x||_{2}^{2} + \frac{1}{2} * || x - g ||_{2}^{2}
+Problem:     min_{u},   \alpha * ||\nabla u||_{2,2} + Fidelity(u, g)
 
              \alpha: Regularization parameter
              
              \nabla: Gradient operator 
              
-             g: Noisy Data with Gaussian Noise
+             g: Noisy Data 
                           
+             Fidelity =  1) L2NormSquarred ( \frac{1}{2} * || u - g ||_{2}^{2} ) if Noise is Gaussian
+                         2) L1Norm ( ||u - g||_{1} )if Noise is Salt & Pepper
+                         3) Kullback Leibler (\int u - g * log(u) + Id_{u>0})  if Noise is Poisson
+                                                       
              Method = 0 ( PDHG - split ) :  K = [ \nabla,
                                                  Identity]
                           
                                                                     
-             Method = 1 (PDHG - explicit ):  K = \nabla  
-                                                                
-"""
+             Method = 1 (PDHG - explicit ):  K = \nabla    
+             
+             
+             Default: Tikhonov denoising
+             noise = Gaussian
+             Fidelity = L2NormSquarred 
+             method = 0
+             
+"""      
 
 from ccpi.framework import ImageData, TestData
 
@@ -62,7 +73,7 @@ else:
 if len(sys.argv) > 1:
     which_noise = int(sys.argv[1])
 else:
-    which_noise = 0
+    which_noise = 2
 print ("Applying {} noise")
 
 if len(sys.argv) > 2:
@@ -71,39 +82,25 @@ else:
     method = '0'
 print ("method ", method)
 
-# Create phantom for TV Salt & Pepper denoising
-N = 100
-
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,N))
+data = loader.load(TestData.SHAPES)
 ig = data.geometry
 ag = ig
 
-# Create noisy data. Apply Salt & Pepper noise
 # Create noisy data. 
-# Apply Salt & Pepper noise
-# gaussian
-# poisson
 noises = ['gaussian', 'poisson', 's&p']
 noise = noises[which_noise]
 if noise == 's&p':
     n1 = random_noise(data.as_array(), mode = noise, salt_vs_pepper = 0.9, amount=0.2)
 elif noise == 'poisson':
-    n1 = random_noise(data.as_array(), mode = noise, seed = 10)
+    scale = 5
+    n1 = random_noise( data.as_array()/scale, mode = noise, seed = 10)*scale
 elif noise == 'gaussian':
     n1 = random_noise(data.as_array(), mode = noise, seed = 10)
 else:
     raise ValueError('Unsupported Noise ', noise)
 noisy_data = ImageData(n1)
 
-# fidelity
-if noise == 's&p':
-    f2 = L1Norm(b=noisy_data)
-elif noise == 'poisson':
-    f2 = KullbackLeibler(noisy_data)
-elif noise == 'gaussian':
-    f2 = 0.5 * L2NormSquared(b=noisy_data)
-
 # Show Ground Truth and Noisy Data
 plt.figure(figsize=(10,5))
 plt.subplot(1,2,1)
@@ -116,15 +113,21 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-
-# Regularisation Parameter
-# no edge preservation alpha is big
+# Regularisation Parameter depending on the noise distribution
+if noise == 's&p':
+    alpha = 20
+elif noise == 'poisson':
+    alpha = 10
+elif noise == 'gaussian':
+    alpha = 5
+    
+# fidelity
 if noise == 's&p':
-    alpha = 8.
+    f2 = L1Norm(b=noisy_data)
 elif noise == 'poisson':
-    alpha = 8.
+    f2 = KullbackLeibler(noisy_data)
 elif noise == 'gaussian':
-    alpha = 8.
+    f2 = 0.5 * L2NormSquared(b=noisy_data)    
 
 if method == '0':
 
@@ -135,21 +138,14 @@ if method == '0':
     # Create BlockOperator
     operator = BlockOperator(op1, op2, shape=(2,1) ) 
 
-    # Create functions
-      
+    # Create functions      
     f1 = alpha * L2NormSquared()
-    # f2 must change according to noise
-    #f2 = 0.5 * L2NormSquared(b = noisy_data)
     f = BlockFunction(f1, f2)                                       
     g = ZeroFunction()
     
 else:
     
-    # Without the "Block Framework"
     operator = Gradient(ig)
-    f =  alpha * L2NormSquared()
-    # g must change according to noise
-    #g =  0.5 * L2NormSquared(b = noisy_data)
     g = f2
         
     
@@ -159,13 +155,12 @@ normK = operator.norm()
 # Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
-opt = {'niter':2000, 'memopt': True}
 
 # Setup and run the PDHG algorithm
-pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
 pdhg.max_iteration = 2000
-pdhg.update_objective_interval = 50
-pdhg.run(1500)
+pdhg.update_objective_interval = 100
+pdhg.run(2000)
 
 
 plt.figure(figsize=(20,5))
@@ -179,78 +174,78 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.subplot(1,4,3)
 plt.imshow(pdhg.get_output().as_array())
-plt.title('Tikhonov Reconstruction')
+plt.title('TV Reconstruction')
 plt.colorbar()
 plt.subplot(1,4,4)
-plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
-plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), data.as_array()[int(ig.shape[0]/2),:], label = 'GTruth')
+plt.plot(np.linspace(0,ig.shape[1],ig.shape[1]), pdhg.get_output().as_array()[int(ig.shape[0]/2),:], label = 'Tikhonov reconstruction')
 plt.legend()
 plt.title('Middle Line Profiles')
 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))
+#    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]))
+#
+#
+#
+#
+#
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py
index 95fe2c1..e7e33cb 100644
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TV_Tomo2D_poisson.py
@@ -84,7 +84,7 @@ sin = Aop.direct(data)
 # Create noisy data. Apply Poisson noise
 scale = 0.5
 eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
+n1 = np.random.poisson(eta + sin.as_array())
 
 noisy_data = AcquisitionData(n1, ag)
 
@@ -100,6 +100,8 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
+
+#%%
 # Regularisation Parameter
 alpha = 2
 
@@ -157,72 +159,72 @@ 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)
-    
-    tmp = noisy_data.as_array().ravel()
-    
-    fidelity = sum(kl_div(tmp, ProjMat * u))    
-    
-    constraints = [q>=fidelity, u>=0]   
-    solver = SCS
-    obj =  Minimize( regulariser +  q)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)    
-             
-    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
-           
-    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), np.reshape(u.value, (N, N))[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)
+#    
+#    tmp = noisy_data.as_array().ravel()
+#    
+#    fidelity = sum(kl_div(tmp, ProjMat * u))    
+#    
+#    constraints = [q>=fidelity, u>=0]   
+#    solver = SCS
+#    obj =  Minimize( regulariser +  q)
+#    prob = Problem(obj, constraints)
+#    result = prob.solve(verbose = True, solver = solver)    
+#             
+#    diff_cvx = np.abs(pdhg.get_output().as_array() - np.reshape(u.value, (N, N)))
+#           
+#    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), np.reshape(u.value, (N, N))[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