fix cgls

2 files changed, 56 insertions, 46 deletions

diff --git a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
index abcb26c..e8f4063 100644
--- a/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
+++ b/Wrappers/Python/demos/CGLS_examples/CGLS_Tomography.py
@@ -50,7 +50,7 @@ import os
 
 # Load  Shepp-Logan phantom 
 model = 1 # select a model number from the library
-N = 64 # set dimension of the phantom
+N = 128 # set dimension of the phantom
 path = os.path.dirname(tomophantom.__file__)
 path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
 phantom_2D = TomoP2D.Model(model, N, path_library2D)
@@ -59,7 +59,7 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
 data = ImageData(phantom_2D)
 
 detectors =  N
-angles = np.linspace(0, np.pi, 90, dtype=np.float32)
+angles = np.linspace(0, np.pi, 180, dtype=np.float32)
 
 ag = AcquisitionGeometry('parallel','2D', angles, detectors)
 
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
index 672d4bc..568df38 100644
--- a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
+++ b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
@@ -32,7 +32,7 @@ Problem:     min_x || A x - g ||_{2}^{2}
 """
 
 
-from ccpi.framework import ImageData, TestData, AcquisitionGeometry
+from ccpi.framework import ImageData, AcquisitionData, ImageGeometry, AcquisitionGeometry
 
 import numpy as np 
 import numpy                          
@@ -43,40 +43,47 @@ from ccpi.optimisation.algorithms import PDHG, CGLS, FISTA
 from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition
 from ccpi.astra.ops import AstraProjectorSimple
 import astra   
-import os, sys
 
+import tomophantom
+from tomophantom import TomoP2D
+import os
 
-# Load Data 
-loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))  
-                     
-N = 50
-M = 50
-data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
-ig = data.geometry
+# Load  Shepp-Logan phantom 
+model = 1 # select a model number from the library
+N = 128 # set dimension of the phantom
+path = os.path.dirname(tomophantom.__file__)
+path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+phantom_2D = TomoP2D.Model(model, N, path_library2D)
 
-detectors = N
-angles = np.linspace(0, np.pi, N, dtype=np.float32)
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+data = ImageData(phantom_2D)
+
+detectors =  N
+angles = np.linspace(0, np.pi, 180, dtype=np.float32)
 
-device = input('Available device: GPU==1 / CPU==0 ')
 ag = AcquisitionGeometry('parallel','2D', angles, detectors)
-if device=='1':
+
+device = input('Available device: GPU==1 / CPU==0 ')
+
+if device =='1':
     dev = 'gpu'
 else:
     dev = 'cpu'
-    
-Aop = AstraProjectorSimple(ig, ag, dev)
+
+Aop = AstraProjectorSimple(ig, ag, dev)    
 sin = Aop.direct(data)
 
-noisy_data = sin 
+np.random.seed(10)
+noisy_data = AcquisitionData( sin.as_array() + np.random.normal(0,1,ag.shape))
 
 ###############################################################################
 # Setup and run Astra CGLS algorithm
 vol_geom = astra.create_vol_geom(N, N)
 proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-proj_id = astra.create_projector('linear', proj_geom, vol_geom)
+proj_id = astra.create_projector('line', proj_geom, vol_geom)
 
 # Create a sinogram from a phantom
-sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id)
+sinogram_id = astra.data2d.create('-sino', proj_geom, noisy_data.as_array())
 
 # Create a data object for the reconstruction
 rec_id = astra.data2d.create('-vol', vol_geom)
@@ -89,17 +96,21 @@ cgls_astra['ProjectorId'] = proj_id
 # Create the algorithm object from the configuration structure
 alg_id = astra.algorithm.create(cgls_astra)
 
-astra.algorithm.run(alg_id, 1000)
+astra.algorithm.run(alg_id, 25)
+recon_cgls_astra = ImageData(astra.data2d.get(rec_id))
 
-recon_cgls_astra = astra.data2d.get(rec_id)
 
-###############################################################################
-# Setup and run the CGLS algorithm  
-x_init = ig.allocate()             
-cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data)
-cgls.max_iteration = 1000
-cgls.update_objective_interval = 200
-cgls.run(1000, verbose=False)
+#%%
+
+# Setup and run the simple CGLS algorithm  
+x_init = ig.allocate()  
+
+cgls = CGLS(x_init = x_init, operator = Aop, data = noisy_data)
+cgls.max_iteration = 25
+cgls.update_objective_interval = 5
+cgls.run(25, verbose = True)
+
+#%%
 
 ###############################################################################
 # Setup and run the PDHG algorithm 
@@ -115,19 +126,20 @@ sigma = 1
 tau = 1/(sigma*normK**2)
 
 pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
-pdhg.max_iteration = 1000
-pdhg.update_objective_interval = 200
-pdhg.run(1000, verbose=True)
+pdhg.max_iteration = 400
+pdhg.update_objective_interval = 20
+pdhg.run(400, verbose=True)
 
+#%%
 ###############################################################################
 # Setup and run the FISTA algorithm 
 fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop)
 regularizer = ZeroFunction()
 
 fista = FISTA(x_init=x_init , f=fidelity, g=regularizer)
-fista.max_iteration = 1000
-fista.update_objective_interval = 200
-fista.run(1000, verbose=True)
+fista.max_iteration = 20
+fista.update_objective_interval = 5
+fista.run(20, verbose = True)
 
 #%% Show results
 
@@ -150,40 +162,38 @@ plt.colorbar()
 plt.title('PDHG reconstruction')
 
 plt.subplot(2,2,4)
-plt.imshow(recon_cgls_astra)
+plt.imshow(recon_cgls_astra.as_array())
 plt.colorbar()
 plt.title('CGLS astra')
 
-diff1 = pdhg.get_output() - cgls.get_output()
-diff2 = fista.get_output() - cgls.get_output()
-diff3 = ImageData(recon_cgls_astra) - cgls.get_output()
+diff1 = pdhg.get_output() - recon_cgls_astra
+diff2 = fista.get_output() - recon_cgls_astra
+diff3 = cgls.get_output() - recon_cgls_astra
 
 plt.figure(figsize=(15,15))
 
 plt.subplot(3,1,1)
 plt.imshow(diff1.abs().as_array())
-plt.title('Diff PDHG vs CGLS')
+plt.title('Diff PDHG vs CGLS astra')
 plt.colorbar()
 
 plt.subplot(3,1,2)
 plt.imshow(diff2.abs().as_array())
-plt.title('Diff FISTA vs CGLS')
+plt.title('Diff FISTA vs CGLS astra')
 plt.colorbar()
 
 plt.subplot(3,1,3)
 plt.imshow(diff3.abs().as_array())
-plt.title('Diff CLGS astra vs CGLS')
+plt.title('Diff CLGS vs CGLS astra')
 plt.colorbar()
 
 
-#%%
+cgls_astra_obj = fidelity(ImageData(recon_cgls_astra))
 
 print('Primal Objective (FISTA) {} '.format(fista.objective[-1]))
 print('Primal Objective (CGLS) {} '.format(cgls.objective[-1]))
 print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
+print('Primal Objective (CGLS_astra) {} '.format(cgls_astra_obj))
 
 
-true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm()
-print('True objective {}'.format(true_obj))
-