comparison algs

2 files changed, 141 insertions, 16 deletions

diff --git a/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py b/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py
index df3c153..2dcaa89 100644
--- a/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py
+++ b/Wrappers/Python/wip/Demos/FISTA_vs_CGLS.py
@@ -44,17 +44,10 @@ detectors = N
 angles = np.linspace(0, np.pi, N, dtype=np.float32)
 
 ag = AcquisitionGeometry('parallel','2D',angles, detectors)
-Aop = AstraProjectorSimple(ig, ag, 'gpu')
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
 sin = Aop.direct(data)
 
-# Create noisy data. Apply Gaussian noise
-
-np.random.seed(10)
-noisy_data = sin + AcquisitionData(np.random.normal(0, 3, sin.shape))
-
-# Regularisation Parameter
-
-operator = Gradient(ig)
+noisy_data = sin 
 
 fidelity = FunctionOperatorComposition(L2NormSquared(b=noisy_data), Aop)
 regularizer = ZeroFunction()
@@ -64,16 +57,21 @@ x_init = ig.allocate()
 ## Setup and run the FISTA algorithm
 opt = {'tol': 1e-4, 'memopt':True}
 fista = FISTA(x_init=x_init , f=fidelity, g=regularizer, opt=opt)
-fista.max_iteration = 20
-fista.update_objective_interval = 1
-fista.run(20, verbose=True)
+fista.max_iteration = 500
+fista.update_objective_interval = 50
+fista.run(500, verbose=True)
 
 ## Setup and run the CGLS algorithm        
 cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data)
-cgls.max_iteration = 20
-cgls.run(20, verbose=True)
+cgls.max_iteration = 500
+cgls.update_objective_interval = 50
+cgls.run(500, verbose=True)
 
-diff = np.abs(fista.get_output().as_array() - cgls.get_output().as_array())
+diff = fista.get_output() - cgls.get_output()
+
+
+#%%
+print( diff.norm())
 
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
@@ -85,7 +83,7 @@ plt.imshow(cgls.get_output().as_array())
 plt.title('CGLS reconstruction')
 plt.colorbar()
 plt.subplot(3,1,3)
-plt.imshow(diff)
+plt.imshow(diff.abs().as_array())
 plt.title('Difference reconstruction')
 plt.colorbar()
 plt.show()
diff --git a/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py b/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py
new file mode 100644
index 0000000..3155654
--- /dev/null
+++ b/Wrappers/Python/wip/Demos/PDHG_vs_CGLS.py
@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from ccpi.framework import ImageData, ImageGeometry, AcquisitionGeometry, AcquisitionData
+
+import numpy as np 
+import numpy                          
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG, CGLS
+
+from ccpi.optimisation.operators import Gradient
+from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, FunctionOperatorComposition
+from skimage.util import random_noise
+from ccpi.astra.ops import AstraProjectorSimple
+            
+#%%
+
+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)
+
+detectors = N
+angles = np.linspace(0, np.pi, N, dtype=np.float32)
+
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
+sin = Aop.direct(data)
+
+noisy_data = sin 
+
+x_init = ig.allocate()
+
+## Setup and run the CGLS algorithm        
+cgls = CGLS(x_init=x_init, operator=Aop, data=noisy_data)
+cgls.max_iteration = 500
+cgls.update_objective_interval = 50
+cgls.run(500, verbose=True)
+
+# Create BlockOperator
+operator = Aop
+f = 0.5 * L2NormSquared(b = noisy_data)
+g = ZeroFunction()
+  
+## Compute operator Norm
+normK = operator.norm()
+
+## Primal & dual stepsizes
+sigma = 0.1
+tau = 1/(sigma*normK**2)
+#
+#
+## Setup and run the PDHG algorithm
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+pdhg.max_iteration = 2000
+pdhg.update_objective_interval = 50
+pdhg.run(2000)
+
+#%%
+
+diff = pdhg.get_output() - cgls.get_output()
+print( diff.norm())
+#
+plt.figure(figsize=(15,15))
+plt.subplot(3,1,1)
+plt.imshow(pdhg.get_output().as_array())
+plt.title('PDHG reconstruction')
+plt.colorbar()
+plt.subplot(3,1,2)
+plt.imshow(cgls.get_output().as_array())
+plt.title('CGLS reconstruction')
+plt.colorbar()
+plt.subplot(3,1,3)
+plt.imshow(diff.abs().as_array())
+plt.title('Difference reconstruction')
+plt.colorbar()
+plt.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+#
+#
+#
+#
+#
+#
+#
+#