fixing other method pf pdhg denoising

2 files changed, 280 insertions, 7 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/algs.py b/Wrappers/Python/ccpi/optimisation/algs.py
index 6b6ae2c..da6adc1 100755
--- a/Wrappers/Python/ccpi/optimisation/algs.py
+++ b/Wrappers/Python/ccpi/optimisation/algs.py
@@ -72,7 +72,7 @@ def FISTA(x_init, f=None, g=None, opt=None):
     
     t_old = 1
     
-    c = f(x_init) + g(x_init)
+#    c = f(x_init) + g(x_init)
 
     # algorithm loop
     for it in range(0, max_iter):
@@ -99,9 +99,9 @@ def FISTA(x_init, f=None, g=None, opt=None):
             
             
         else:
-            u = y - invL*f.grad(y)
+            u = y - invL*f.gradient(y)
             
-            x = g.prox(u,invL)
+            x = g.proximal(u,invL)
             
             t = 0.5*(1 + numpy.sqrt(1 + 4*(t_old**2)))
             
@@ -111,8 +111,8 @@ def FISTA(x_init, f=None, g=None, opt=None):
             t_old = t
         
         # time and criterion
-        timing[it] = time.time() - time0
-        criter[it] = f(x) + g(x);
+#        timing[it] = time.time() - time0
+#        criter[it] = f(x) + g(x);
         
         # stopping rule
         #if np.linalg.norm(x - x_old) < tol * np.linalg.norm(x_old) and it > 10:
@@ -121,9 +121,9 @@ def FISTA(x_init, f=None, g=None, opt=None):
         #print(it, 'out of', 10, 'iterations', end='\r');
 
     #criter = criter[0:it+1];
-    timing = numpy.cumsum(timing[0:it+1]);
+#    timing = numpy.cumsum(timing[0:it+1]);
     
-    return x, it, timing, criter
+    return x #, it, timing, criter
 
 def FBPD(x_init, operator=None, constraint=None, data_fidelity=None,\
          regulariser=None, opt=None):
diff --git a/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py b/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py
new file mode 100644
index 0000000..e142d94
--- /dev/null
+++ b/Wrappers/Python/wip/test_pdhg_profile/profile_pdhg_TV_denoising.py
@@ -0,0 +1,273 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb 22 14:53:03 2019
+
+@author: evangelos
+"""
+
+from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer
+
+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 ZeroFunction, L2NormSquared, \
+                      MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction
+
+from skimage.util import random_noise
+
+from timeit import default_timer as timer
+def dt(steps):
+    return steps[-1] - steps[-2]
+
+#%%
+
+# Create phantom for TV denoising
+
+N = 500
+
+data = np.zeros((N,N))
+data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+data[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
+ag = ig
+
+# Create noisy data. Add Gaussian noise
+n1 = random_noise(data, mode = 'gaussian', mean=0, var = 0.05, seed=10)
+noisy_data = ImageData(n1)
+
+plt.imshow(noisy_data.as_array())
+plt.show()
+
+#%%
+
+# Regularisation Parameter
+alpha = 2
+
+#method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
+method = '0'
+
+if method == '0':
+
+    # Create operators
+    op1 = Gradient(ig)
+    op2 = Identity(ig, ag)
+
+    # Form Composite Operator
+    operator = BlockOperator(op1, op2, shape=(2,1) ) 
+
+    #### Create functions
+      
+    f1 = alpha * MixedL21Norm()
+    f2 = 0.5 * L2NormSquared(b = noisy_data)    
+    f = BlockFunction(f1, f2)  
+                                      
+    g = ZeroFunction()
+    
+else:
+    
+    ###########################################################################
+    #         No Composite #
+    ###########################################################################
+    operator = Gradient(ig)
+    f = alpha * FunctionOperatorComposition(operator, MixedL21Norm())
+    g = L2NormSquared(b = noisy_data)
+    
+    ###########################################################################
+#%%
+    
+# Compute operator Norm
+normK = operator.norm()
+
+# Primal & dual stepsizes
+sigma = 1
+tau = 1/(sigma*normK**2)
+
+opt = {'niter':2000}
+opt1 = {'niter':2000, 'memopt': True}
+
+t1 = timer()
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+t2 = timer()
+
+
+t3 = timer()
+res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt1) 
+t4 = timer()
+#
+print ("No memopt in {}s, memopt in  {}s ".format(t2-t1, t4 -t3))
+
+#    
+#%%
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 100
+##
+#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+#pdhgo.max_iteration = 2000
+#pdhgo.update_objective_interval = 100
+##
+#steps = [timer()]
+#pdhgo.run(2000)
+#steps.append(timer())
+#t1 = dt(steps)
+##
+#pdhg.run(2000)
+#steps.append(timer())
+#t2 = dt(steps)
+#
+#print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1))
+#res = pdhg.get_output()
+#res1 = pdhgo.get_output()
+
+#%%
+#plt.figure(figsize=(15,15))
+#plt.subplot(3,1,1)
+#plt.imshow(res.as_array())
+#plt.title('no memopt')
+#plt.colorbar()
+#plt.subplot(3,1,2)
+#plt.imshow(res1.as_array())
+#plt.title('memopt')
+#plt.colorbar()
+#plt.subplot(3,1,3)
+#plt.imshow((res1 - res).abs().as_array())
+#plt.title('diff')
+#plt.colorbar()
+#plt.show()
+
+
+#plt.figure(figsize=(15,15))
+#plt.subplot(3,1,1)
+#plt.imshow(pdhg.get_output().as_array())
+#plt.title('no memopt class')
+#plt.colorbar()
+#plt.subplot(3,1,2)
+#plt.imshow(res.as_array())
+#plt.title('no memopt')
+#plt.colorbar()
+#plt.subplot(3,1,3)
+#plt.imshow((pdhg.get_output() - res).abs().as_array())
+#plt.title('diff')
+#plt.colorbar()
+#plt.show()
+#    
+#    
+#
+#plt.figure(figsize=(15,15))
+#plt.subplot(3,1,1)
+#plt.imshow(pdhgo.get_output().as_array())
+#plt.title('memopt class')
+#plt.colorbar()
+#plt.subplot(3,1,2)
+#plt.imshow(res1.as_array())
+#plt.title('no memopt')
+#plt.colorbar()
+#plt.subplot(3,1,3)
+#plt.imshow((pdhgo.get_output() - res1).abs().as_array())
+#plt.title('diff')
+#plt.colorbar()
+#plt.show()
+    
+
+    
+
+    
+#    print ("Time difference {}s {}s {}s Speedup {:.2f}".format(t1,t2,t2-t1, t2/t1))
+#    res = pdhg.get_output()
+#    res1 = pdhgo.get_output()
+#    
+#    diff = (res-res1)
+#    print ("diff norm {} max {}".format(diff.norm(), diff.abs().as_array().max()))
+#    print ("Sum ( abs(diff) )  {}".format(diff.abs().sum()))
+#    
+#    
+#    plt.figure(figsize=(5,5))
+#    plt.subplot(1,3,1)
+#    plt.imshow(res.as_array())
+#    plt.colorbar()
+#    #plt.show()
+#     
+#    #plt.figure(figsize=(5,5))
+#    plt.subplot(1,3,2)
+#    plt.imshow(res1.as_array())
+#    plt.colorbar()
+    
+#plt.show()
+
+
+
+#=======
+## opt = {'niter':2000, 'memopt': True}
+#
+## res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+# 
+#>>>>>>> origin/pdhg_fix
+#
+#
+## opt = {'niter':2000, 'memopt': False}
+## res1, time1, primal1, dual1, pdgap1 = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+#
+## plt.figure(figsize=(5,5))
+## plt.subplot(1,3,1)
+## plt.imshow(res.as_array())
+## plt.title('memopt')
+## plt.colorbar()
+## plt.subplot(1,3,2)
+## plt.imshow(res1.as_array())
+## plt.title('no memopt')
+## plt.colorbar()
+## plt.subplot(1,3,3)
+## plt.imshow((res1 - res).abs().as_array())
+## plt.title('diff')
+## plt.colorbar()
+## plt.show()
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 100
+#
+#
+#pdhgo = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma, memopt=True)
+#pdhgo.max_iteration = 2000
+#pdhgo.update_objective_interval = 100
+#
+#steps = [timer()]
+#pdhgo.run(200)
+#steps.append(timer())
+#t1 = dt(steps)
+#
+#pdhg.run(200)
+#steps.append(timer())
+#t2 = dt(steps)
+#
+#print ("Time difference {} {} {}".format(t1,t2,t2-t1))
+#sol = pdhg.get_output().as_array()
+##sol = result.as_array()
+##
+#fig = plt.figure()
+#plt.subplot(1,3,1)
+#plt.imshow(noisy_data.as_array())
+#plt.colorbar()
+#plt.subplot(1,3,2)
+#plt.imshow(sol)
+#plt.colorbar()
+#plt.subplot(1,3,3)
+#plt.imshow(pdhgo.get_output().as_array())
+#plt.colorbar()
+#
+#plt.show()
+###
+##
+####
+##plt.plot(np.linspace(0,N,N), data[int(N/2),:], label = 'GTruth')
+##plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+##plt.legend()
+##plt.show()
+#
+#
+##%%
+##