fix precond

3 files changed, 38 insertions, 22 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
index a58a296..a853b8d 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/IdentityOperator.py
@@ -50,7 +50,7 @@ class Identity(LinearOperator):
     
     def sum_abs_row(self):
         
-        return self.gm_domain.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
+        return self.gm_range.allocate(1)#ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
  
     def sum_abs_col(self):
         
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 5e318ff..c5c2ec9 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -65,13 +65,13 @@ class SparseFiniteDiff():
     def sum_abs_row(self):
         
         res = np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F'))
-        res[res==0]=1
+        #res[res==0]=0
         return ImageData(res)
     
     def sum_abs_col(self):
         
         res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F') )
-        res[res==0]=1
+        #res[res==0]=0
         return ImageData(res)
         
 if __name__ == '__main__':
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
index 39bbb2c..860e76e 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Gaussian.py
@@ -17,6 +17,28 @@
 #   See the License for the specific language governing permissions and
 #   limitations under the License.
 
+
+""" 
+
+Total Variation Denoising using PDHG algorithm:
+
+             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
+
+
+Problem:     min_x \alpha * ||\nabla x||_{1} + || x - g ||_{2}^{2}
+
+             \nabla: Gradient operator 
+             g: Noisy Data with Gaussian Noise
+             \alpha: Regularization parameter
+             
+             Method = 0:  K = [ \nabla,
+                                 Identity]
+                                                                    
+             Method = 1:  K = \nabla    
+             
+             
+"""
+
 from ccpi.framework import ImageData, ImageGeometry
 
 import numpy as np 
@@ -29,10 +51,9 @@ from ccpi.optimisation.operators import BlockOperator, Identity, Gradient
 from ccpi.optimisation.functions import ZeroFunction, L2NormSquared, \
                       MixedL21Norm, BlockFunction
 
-from skimage.util import random_noise
 
 # Create phantom for TV Gaussian denoising
-N = 100
+N = 200
 
 data = np.zeros((N,N))
 data[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
@@ -42,13 +63,13 @@ ig = ImageGeometry(voxel_num_x = N, voxel_num_y = N)
 ag = ig
 
 # Create noisy data. Add Gaussian noise
-n1 = random_noise(data.as_array(), mode = 'gaussian', mean=0, var = 0.05, seed=10)
-noisy_data = ImageData(n1)
+np.random.seed(10)
+noisy_data = ImageData( data.as_array() + np.random.normal(0, 0.05, size=ig.shape) )
 
 # Regularisation Parameter
 alpha = 2
 
-method = '1'
+method = '0'
 
 if method == '0':
 
@@ -70,11 +91,10 @@ if method == '0':
 else:
     
     # Without the "Block Framework"
-    operator = Gradient(ig)
+    operator =  Gradient(ig)
     f =  alpha * MixedL21Norm()
     g =  0.5 * L2NormSquared(b = noisy_data)
         
-    
 # Compute operator Norm
 normK = operator.norm()
 
@@ -82,13 +102,11 @@ normK = operator.norm()
 sigma = 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)
-
+pdhg.max_iteration = 3000
+pdhg.update_objective_interval = 200
+pdhg.run(3000, verbose=False)
 
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
@@ -104,7 +122,7 @@ plt.imshow(pdhg.get_output().as_array())
 plt.title('TV Reconstruction')
 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), pdhg.get_output().as_array()[int(N/2),:], label = 'TV reconstruction')
 plt.legend()
@@ -112,7 +130,7 @@ plt.title('Middle Line Profiles')
 plt.show()
 
 
-##%% Check with CVX solution
+#%% Check with CVX solution
 
 from ccpi.optimisation.operators import SparseFiniteDiff
 
@@ -143,7 +161,7 @@ if cvx_not_installable:
         
     obj =  Minimize( regulariser +  fidelity)
     prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
+    result = prob.solve(verbose = True, solver = MOSEK)
     
     diff_cvx = numpy.abs( pdhg.get_output().as_array() - u.value )
         
@@ -164,6 +182,8 @@ if cvx_not_installable:
     
     plt.plot(np.linspace(0,N,N), pdhg.get_output().as_array()[int(N/2),:], label = 'PDHG')
     plt.plot(np.linspace(0,N,N), u.value[int(N/2),:], label = 'CVX')
+    plt.plot(np.linspace(0,N,N), data.as_array()[int(N/2),:], label = 'Truth')
+    
     plt.legend()
     plt.title('Middle Line Profiles')
     plt.show()
@@ -171,7 +191,3 @@ if cvx_not_installable:
     print('Primal Objective (CVX) {} '.format(obj.value))
     print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
 
-
-
-
-