fista trials

5 files changed, 61 insertions, 79 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
index 04e7c38..419958a 100755
--- a/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/FISTA.py
@@ -67,7 +67,8 @@ class FISTA(Algorithm):
         
         self.t = 0.5*(1 + numpy.sqrt(1 + 4*(self.t_old**2)))
         
-        self.x.subtract(self.x_old, out=self.y)
+#        self.x.subtract(self.x_old, out=self.y)
+        self.y = self.x - self.x_old
         self.y.__imul__ ((self.t_old-1)/self.t)
         self.y.__iadd__( self.x )
         
diff --git a/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py b/Wrappers/Python/demos/CompareAlgorithms/CGLS_FISTA_PDHG_LeastSquares.py
index 0875c20..672d4bc 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, ImageGeometry, AcquisitionGeometry
+from ccpi.framework import ImageData, TestData, AcquisitionGeometry
 
 import numpy as np 
 import numpy                          
@@ -42,17 +42,17 @@ 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 astra   
+import os, sys
 
-# Create Ground truth phantom and Sinogram
 
-N = 128
-x = np.zeros((N,N))
-x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
-x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
-
-data = ImageData(x)
-ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+# 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
 
 detectors = N
 angles = np.linspace(0, np.pi, N, dtype=np.float32)
@@ -69,13 +69,14 @@ sin = Aop.direct(data)
 
 noisy_data = sin 
 
+###############################################################################
 # 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('strip', proj_geom, vol_geom)
+proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
 # Create a sinogram from a phantom
-sinogram_id, sinogram = astra.create_sino(x, proj_id)
+sinogram_id, sinogram = astra.create_sino(data.as_array(), proj_id)
 
 # Create a data object for the reconstruction
 rec_id = astra.data2d.create('-vol', vol_geom)
@@ -92,6 +93,7 @@ astra.algorithm.run(alg_id, 1000)
 
 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)
@@ -99,12 +101,15 @@ cgls.max_iteration = 1000
 cgls.update_objective_interval = 200
 cgls.run(1000, verbose=False)
 
+###############################################################################
 # Setup and run the PDHG algorithm 
 operator = Aop
 f = L2NormSquared(b = noisy_data)
 g = ZeroFunction()
+
 ## Compute operator Norm
 normK = operator.norm()
+
 ## Primal & dual stepsizes
 sigma = 1
 tau = 1/(sigma*normK**2)
@@ -112,17 +117,17 @@ 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=False)
+pdhg.run(1000, verbose=True)
 
+###############################################################################
 # Setup and run the FISTA algorithm 
 fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop)
 regularizer = ZeroFunction()
 
-opt = {'memopt':True}
-fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt)
+fista = FISTA(x_init=x_init , f=fidelity, g=regularizer)
 fista.max_iteration = 1000
 fista.update_objective_interval = 200
-fista.run(1000, verbose=False)
+fista.run(1000, verbose=True)
 
 #%% Show results
 
@@ -131,18 +136,22 @@ plt.suptitle('Reconstructions ', fontsize=16)
 
 plt.subplot(2,2,1)
 plt.imshow(cgls.get_output().as_array())
+plt.colorbar()
 plt.title('CGLS reconstruction')
 
 plt.subplot(2,2,2)
 plt.imshow(fista.get_output().as_array())
+plt.colorbar()
 plt.title('FISTA reconstruction')
 
 plt.subplot(2,2,3)
 plt.imshow(pdhg.get_output().as_array())
+plt.colorbar()
 plt.title('PDHG reconstruction')
 
 plt.subplot(2,2,4)
 plt.imshow(recon_cgls_astra)
+plt.colorbar()
 plt.title('CGLS astra')
 
 diff1 = pdhg.get_output() - cgls.get_output()
@@ -167,28 +176,14 @@ plt.title('Diff CLGS astra vs CGLS')
 plt.colorbar()
 
 
+#%%
 
+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]))
 
 
+true_obj = (Aop.direct(cglsd.get_output())-noisy_data).squared_norm()
+print('True objective {}'.format(true_obj))
 
 
-
-
-
-
-
-
-
-
-
-
-
-
-#
-#
-#
-#
-#
-#
-#
-#
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
index 68fb649..1a96a2d 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_CGLS.py
@@ -21,14 +21,15 @@
 
 from ccpi.framework import AcquisitionGeometry
 from ccpi.optimisation.algorithms import FISTA
-from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, L2NormSquared, FunctionOperatorComposition
-from ccpi.astra.ops import AstraProjectorSimple
+from ccpi.optimisation.functions import IndicatorBox, ZeroFunction, \
+                         L2NormSquared, FunctionOperatorComposition
+from ccpi.astra.operators import AstraProjectorSimple
 
 import numpy as np
 import matplotlib.pyplot as plt
 from ccpi.framework import TestData
 import os, sys
-
+from ccpi.optimisation.funcs import Norm2sq
 
 # Load Data 
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
@@ -117,16 +118,15 @@ plt.show()
 # demonstrated in the rest of this file. In general all methods need an initial 
 # guess and some algorithm options to be set:
 
-f = FunctionOperatorComposition(0.5 * L2NormSquared(b=sin), Aop)
+f = FunctionOperatorComposition(L2NormSquared(b=sin), Aop)
+#f = Norm2sq(Aop, sin, c=0.5, memopt=True)
 g = ZeroFunction()
 
 x_init = ig.allocate()
-opt = {'memopt':True}
-
-fista = FISTA(x_init=x_init, f=f, g=g, opt=opt)
-fista.max_iteration = 100
-fista.update_objective_interval = 1
-fista.run(100, verbose=True)
+fista = FISTA(x_init=x_init, f=f, g=g)
+fista.max_iteration = 1000
+fista.update_objective_interval = 100
+fista.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista.get_output().as_array())
@@ -134,18 +134,12 @@ plt.title('FISTA unconstrained')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista.objective)
-plt.title('FISTA objective')
-plt.show()
-
-#%%
 
 # Run FISTA for least squares with lower/upper bound 
-fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1), opt=opt)
-fista0.max_iteration = 100
-fista0.update_objective_interval = 1
-fista0.run(100, verbose=True)
+fista0 = FISTA(x_init=x_init, f=f, g=IndicatorBox(lower=0,upper=1))
+fista0.max_iteration = 1000
+fista0.update_objective_interval = 100
+fista0.run(1000, verbose=True)
 
 plt.figure()
 plt.imshow(fista0.get_output().as_array())
@@ -153,10 +147,10 @@ plt.title('FISTA constrained in [0,1]')
 plt.colorbar()
 plt.show()
 
-plt.figure()
-plt.semilogy(fista0.objective)
-plt.title('FISTA constrained in [0,1]')
-plt.show()
+#plt.figure()
+#plt.semilogy(fista0.objective)
+#plt.title('FISTA constrained in [0,1]')
+#plt.show()
 
 #%% Check with CVX solution
 
@@ -178,10 +172,10 @@ if cvx_not_installable:
     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)
+    proj_id = astra.create_projector('linear', proj_geom, vol_geom)
 
     matrix_id = astra.projector.matrix(proj_id)
-
+    
     ProjMat = astra.matrix.get(matrix_id)
     
     tmp = sin.as_array().ravel()
@@ -216,4 +210,6 @@ if cvx_not_installable:
     plt.legend()
     plt.title('Middle Line Profiles')
     plt.show()
-            
\ No newline at end of file
+            
+    print('Primal Objective (CVX) {} '.format(obj.value))
+    print('Primal Objective (FISTA) {} '.format(fista0.loss[1]))    
\ No newline at end of file
diff --git a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
index 41219f4..6007990 100644
--- a/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
+++ b/Wrappers/Python/demos/FISTA_examples/FISTA_Tikhonov_Poisson_Denoising.py
@@ -21,7 +21,7 @@
 
 """ 
 
-Tikhonov for Poisson denoising using FISTA algorithm:
+"Tikhonov regularization" for Poisson denoising using FISTA algorithm:
 
 Problem:     min_x, x>0  \alpha * ||\nabla x||_{2}^{2} + \int x - g * log(x) 
 
@@ -52,14 +52,14 @@ from skimage.util import random_noise
 loader = TestData(data_dir=os.path.join(sys.prefix, 'share','ccpi'))
 
 # Load Data                      
-N = 150
-M = 150
+N = 100
+M = 100
 data = loader.load(TestData.SIMPLE_PHANTOM_2D, size=(N,M), scale=(0,1))
 
 ig = data.geometry
 ag = ig
 
-# Create Noisy data. Add Gaussian noise
+# Create Noisy data with Poisson noise
 n1 = random_noise(data.as_array(), mode = 'poisson', seed = 10)
 noisy_data = ImageData(n1)
 
@@ -75,7 +75,6 @@ plt.title('Noisy Data')
 plt.colorbar()
 plt.show()
 
-#%%
 # Regularisation Parameter
 alpha = 10
 
@@ -114,8 +113,7 @@ fid.proximal = KL_Prox_PosCone
 reg = FunctionOperatorComposition(alpha * L2NormSquared(), operator)
 
 x_init = ig.allocate()
-opt = {'memopt':True}
-fista = FISTA(x_init=x_init , f=reg, g=fid, opt=opt)
+fista = FISTA(x_init=x_init , f=reg, g=fid)
 fista.max_iteration = 2000
 fista.update_objective_interval = 500
 fista.run(2000, verbose=True)
@@ -196,7 +194,6 @@ if cvx_not_installable:
     plt.title('Middle Line Profiles')
     plt.show()
             
-    #TODO what is the output of fista.objective, fista.loss
     print('Primal Objective (CVX) {} '.format(obj.value))
     print('Primal Objective (FISTA) {} '.format(fista.loss[1]))
 
diff --git a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
index bc0081f..422deb9 100644
--- a/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
+++ b/Wrappers/Python/demos/PDHG_examples/Tomo/PDHG_TGV_Tomo2D.py
@@ -85,15 +85,8 @@ data = ImageData(data/data.max())
 ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
 ag = ig
 
-# Create noisy data. Add Gaussian noise
-scale = 0.5
-eta = 0 
-n1 = scale * np.random.poisson(eta + sin.as_array()/scale)
-
-noisy_data = AcquisitionData(n1, ag)
 
 #Create Acquisition Data and apply poisson noise
-
 detectors = N
 angles = np.linspace(0, np.pi, N)