update L1Norm add tests

2 files changed, 54 insertions, 25 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
new file mode 100644
index 0000000..f47c655
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/operators/ShrinkageOperator.py
@@ -0,0 +1,19 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Mar  6 19:30:51 2019
+
+@author: evangelos
+"""
+
+from ccpi.framework import DataContainer
+
+class ShrinkageOperator():
+    
+    def __init__(self):
+        pass
+
+    def __call__(self, x, tau, out=None):
+        
+        return x.sign() * (x.abs() - tau).maximum(0) 
+   
\ No newline at end of file
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising.py b/Wrappers/Python/wip/pdhg_TV_denoising.py
index a8e721f..e9787ac 100755
--- a/Wrappers/Python/wip/pdhg_TV_denoising.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising.py
@@ -11,7 +11,7 @@ from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer
 import numpy as np                           
 import matplotlib.pyplot as plt
 
-from ccpi.optimisation.algorithms import PDHG
+from ccpi.optimisation.algorithms import PDHG, PDHG_old
 
 from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
 from ccpi.optimisation.functions import ZeroFun, L2NormSquared, \
@@ -33,7 +33,7 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
 ag = ig
 
 # Create noisy data. Add Gaussian noise
-n1 = random_noise(data, mode='gaussian', seed=10)
+n1 = random_noise(data, mode = 'gaussian', seed=10)
 noisy_data = ImageData(n1)
 
 
@@ -43,7 +43,8 @@ noisy_data = ImageData(n1)
 alpha = 2
 
 #method = input("Enter structure of PDHG (0=Composite or 1=NotComposite): ")
-method = '0'
+method = '1'
+
 if method == '0':
 
     # Create operators
@@ -70,7 +71,7 @@ else:
     ###########################################################################
     operator = Gradient(ig)
     f = alpha * FunctionOperatorComposition(operator, MixedL21Norm())
-    g = 0.5 * L2NormSquared(b = noisy_data)
+    g = L2NormSquared(b = noisy_data)
     ###########################################################################
 #%%
     
@@ -81,32 +82,41 @@ print ("normK", normK)
 sigma = 1
 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 = 10
-
-pdhg.run(2000)
+opt = {'niter':2000}
 
-    
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+ 
+plt.figure(figsize=(5,5))
+plt.imshow(res.as_array())
+plt.colorbar()
+plt.show()
 
-sol = pdhg.get_output().as_array()
-#sol = result.as_array()
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 10
 #
-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()
+#pdhg.run(2000)
 #
-
+#    
+#
+#sol = pdhg.get_output().as_array()
+##sol = result.as_array()
 ##
-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()
+#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[int(N/2),:], label = 'GTruth')
+#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+#plt.legend()
+#plt.show()
 
 
 #%%