precond test

3 files changed, 143 insertions, 8 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
index 752fd21..484dc61 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/BlockOperator.py
@@ -8,7 +8,7 @@ Created on Thu Feb 14 12:36:40 2019
 import numpy
 from numbers import Number
 import functools
-from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer
+from ccpi.framework import AcquisitionData, ImageData, BlockDataContainer, DataContainer
 from ccpi.optimisation.operators import Operator, LinearOperator
 from ccpi.optimisation.operators.BlockScaledOperator import BlockScaledOperator
 from ccpi.framework import BlockGeometry
@@ -224,7 +224,7 @@ class BlockOperator(Operator):
         
         res = []
         for row in range(self.shape[0]):
-            for col in range(self.shape[1]):
+            for col in range(self.shape[1]):                            
                 if col == 0:
                     prod = self.get_item(row,col).sum_abs_row()
                 else:
@@ -269,8 +269,8 @@ if __name__ == '__main__':
     
     B = BlockOperator(G, Id)
     
-    print(B.sum_abs_row().as_array())
-    
+    print(B.sum_abs_row())
+#    
     Gx = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
     Gy = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
     
@@ -282,17 +282,17 @@ if __name__ == '__main__':
     d_res = numpy.reshape(d1 + d2 + d3, ig.shape, 'F')
     
     print(d_res)
-    
+#    
     z1 = abs(Gx.matrix()).toarray().sum(axis=1)
     z2 = abs(Gy.matrix()).toarray().sum(axis=1)
     z3 = abs(Id.matrix()).toarray().sum(axis=1)
-
+#
     z_res = BlockDataContainer(BlockDataContainer(ImageData(numpy.reshape(z2, ig.shape, 'F')),\
                                                   ImageData(numpy.reshape(z1, ig.shape, 'F'))),\
                                                   ImageData(numpy.reshape(z3, ig.shape, 'F')))
-
+#
     ttt = B.sum_abs_col()
-    
+#    
     numpy.testing.assert_array_almost_equal(z_res[0][0].as_array(), ttt[0][0].as_array(), decimal=4)    
     numpy.testing.assert_array_almost_equal(z_res[0][1].as_array(), ttt[0][1].as_array(), decimal=4)    
     numpy.testing.assert_array_almost_equal(z_res[1].as_array(), ttt[1].as_array(), decimal=4)    
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index 87723f0..cd65ee4 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -122,6 +122,7 @@ if __name__ == '__main__':
     d2 = G_neum.sum_abs_col()
     print(d2)    
     
+    d1 * d2
     
     
     
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
new file mode 100644
index 0000000..518ead2
--- /dev/null
+++ b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
@@ -0,0 +1,134 @@
+#!/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 ZeroFun, L2NormSquared, \
+                      MixedL21Norm, FunctionOperatorComposition, BlockFunction, ScaledFunction
+
+from skimage.util import random_noise
+
+
+
+# ############################################################################
+# Create phantom for TV denoising
+
+N = 100
+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', seed=10)
+noisy_data = ImageData(n1)
+
+
+#%%
+
+# 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
+#    f = FunctionComposition_new(operator, mixed_L12Norm(alpha), \
+#                                    L2NormSq(0.5, b = noisy_data) )    
+    
+    f1 = alpha * MixedL21Norm()
+    f2 = 0.5 * L2NormSquared(b = noisy_data)
+    
+    f = BlockFunction(f1, f2 )                                        
+    g = ZeroFun()
+    
+else:
+    
+    ###########################################################################
+    #         No Composite #
+    ###########################################################################
+    operator = Gradient(ig)
+    f = alpha * FunctionOperatorComposition(operator, MixedL21Norm())
+    g = 0.5 * L2NormSquared(b = noisy_data)
+    ###########################################################################
+#%%
+
+diag_precon = True 
+
+if diag_precon:
+    tmp_tau = operator.sum_abs_row()
+    tmp_sigma = operator.sum_abs_col()
+    
+    tmp_sigma[0][0].as_array()[tmp_sigma[0][0].as_array()==0]=1
+    tmp_sigma[0][1].as_array()[tmp_sigma[0][1].as_array()==0]=1
+    tmp_sigma[1].as_array()[tmp_sigma[1].as_array()==0]=1
+    
+    tau = 1/tmp_tau
+    sigma = 1/tmp_sigma
+    
+else:
+    # Compute operator Norm
+    normK = operator.norm()
+    print ("normK", normK)
+    # Primal & dual stepsizes
+    sigma = 1
+    tau = 1/(sigma*normK**2)
+
+#%%
+    
+#opt = {'niter':2000}
+
+#res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt)    
+    
+    
+#pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
+#pdhg.max_iteration = 2000
+#pdhg.update_objective_interval = 10
+#
+#pdhg.run(2000)
+#
+#    
+#
+#sol = pdhg.get_output().as_array()
+##sol = result.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[int(N/2),:], label = 'GTruth')
+#plt.plot(np.linspace(0,N,N), sol[int(N/2),:], label = 'Recon')
+#plt.legend()
+#plt.show()
+#
+
+#%%
+#