add old pdhg, test gap

3 files changed, 199 insertions, 0 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index 043fe38..8600e07 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -13,6 +13,9 @@ import time
 from ccpi.optimisation.operators import BlockOperator
 from ccpi.framework import BlockDataContainer
 
+
+import matplotlib.pyplot as plt
+
 class PDHG(Algorithm):
     '''Primal Dual Hybrid Gradient'''
 
@@ -80,3 +83,69 @@ class PDHG(Algorithm):
         ])
 
 
+
+def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
+        
+    # algorithmic parameters
+    if opt is None: 
+        opt = {'tol': 1e-6, 'niter': 500, 'show_iter': 100, \
+               'memopt': False} 
+        
+    if sigma is None and tau is None:
+        raise ValueError('Need sigma*tau||K||^2<1') 
+                
+    niter = opt['niter'] if 'niter' in opt.keys() else 1000
+    tol = opt['tol'] if 'tol' in opt.keys() else 1e-4
+    memopt = opt['memopt'] if 'memopt' in opt.keys() else False  
+    show_iter = opt['show_iter'] if 'show_iter' in opt.keys() else False 
+    stop_crit = opt['stop_crit'] if 'stop_crit' in opt.keys() else False 
+
+
+    x_old = operator.domain_geometry().allocate()
+    y_old = operator.range_geometry().allocate()       
+        
+    
+    xbar = x_old
+    x_tmp = x_old
+    x = x_old
+    
+    y_tmp = y_old
+    y = y_tmp
+        
+    # relaxation parameter
+    theta = 1
+    
+    t = time.time()
+    
+    objective = []
+    
+    for i in range(niter):
+        
+        # Gradient descent, Dual problem solution
+        y_tmp = y_old + sigma * operator.direct(xbar)
+        y = f.proximal_conjugate(y_tmp, sigma)
+        
+        # Gradient ascent, Primal problem solution
+        x_tmp = x_old - tau * operator.adjoint(y)
+        x = g.proximal(x_tmp, tau)
+        
+        #Update
+        xbar = x + theta * (x - x_old)
+                                
+        x_old = x
+        y_old = y   
+        
+        if i%100==0:
+            
+            primal = f(operator.direct(x)) + g(x)
+            dual = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y)))
+            print( i, primal, dual)
+            
+            plt.imshow(x.as_array())
+            plt.show()
+#            print(f(operator.direct(x)) + g(x), i)
+                         
+    t_end = time.time()        
+        
+    return x, t_end - t, objective
+
diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
index 443bc78..a28c0bf 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/__init__.py
@@ -28,3 +28,4 @@ from .GradientDescent import GradientDescent
 from .FISTA import FISTA
 from .FBPD import FBPD
 from .PDHG import PDHG
+from .PDHG import PDHG_old
diff --git a/Wrappers/Python/wip/test_pdhg_gap.py b/Wrappers/Python/wip/test_pdhg_gap.py
new file mode 100644
index 0000000..b196e36
--- /dev/null
+++ b/Wrappers/Python/wip/test_pdhg_gap.py
@@ -0,0 +1,129 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  2 12:26:24 2019
+
+@author: vaggelis
+"""
+
+
+from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, AcquisitionGeometry, AcquisitionData
+
+import numpy as np                           
+import matplotlib.pyplot as plt
+
+from ccpi.optimisation.algorithms import PDHG, PDHG_old
+
+from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
+from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \
+                      MixedL21Norm, BlockFunction, ScaledFunction
+
+from ccpi.astra.ops import AstraProjectorSimple
+from skimage.util import random_noise
+
+
+#%%###############################################################################
+# Create phantom for TV tomography
+
+#import os
+#import tomophantom
+#from tomophantom import TomoP2D
+#from tomophantom.supp.qualitymetrics import QualityTools
+
+#model = 1 # select a model number from the library
+#N = 150 # set dimension of the phantom
+## one can specify an exact path to the parameters file
+## path_library2D = '../../../PhantomLibrary/models/Phantom2DLibrary.dat'
+#path = os.path.dirname(tomophantom.__file__)
+#path_library2D = os.path.join(path, "Phantom2DLibrary.dat")
+##This will generate a N_size x N_size phantom (2D)
+#phantom_2D = TomoP2D.Model(model, N, path_library2D)
+#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+#data = ImageData(phantom_2D, geometry=ig)
+
+N = 150
+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 = 150
+angles = np.linspace(0,np.pi,100)
+
+ag = AcquisitionGeometry('parallel','2D',angles, detectors)
+Aop = AstraProjectorSimple(ig, ag, 'cpu')
+sin = Aop.direct(data)
+
+plt.imshow(sin.as_array())
+plt.title('Sinogram')
+plt.colorbar()
+plt.show()
+
+# Add Gaussian noise to the sinogram data
+np.random.seed(10)
+n1 = np.random.random(sin.shape)
+
+noisy_data = sin + ImageData(5*n1)
+
+plt.imshow(noisy_data.as_array())
+plt.title('Noisy Sinogram')
+plt.colorbar()
+plt.show()
+
+
+#%% Works only with Composite Operator Structure of PDHG
+
+#ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+
+# Create operators
+op1 = Gradient(ig)
+op2 = Aop
+
+# Form Composite Operator
+operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+alpha = 50
+f = BlockFunction( alpha * MixedL21Norm(), \
+                   0.5 * L2NormSquared(b = noisy_data) )
+g = ZeroFun()
+
+# Compute operator Norm
+normK = operator.norm()
+
+## Primal & dual stepsizes
+
+sigma = 10
+tau = 1/(sigma*normK**2)
+
+pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+pdhg.max_iteration = 2000
+pdhg.update_objective_interval = 100
+
+#pdhg.run(5000)
+
+opt = {'niter':2000}
+
+res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt)
+
+#%%
+#sol = pdhg.get_output().as_array()
+#fig = plt.figure()
+#plt.subplot(1,2,1)
+#plt.imshow(noisy_data.as_array())
+##plt.colorbar()
+#plt.subplot(1,2,2)
+#plt.imshow(sol)
+##plt.colorbar()
+#plt.show()
+#
+#
+##%%
+#plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'GTruth')
+#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+#plt.legend()
+#plt.show()
+
+