add docstring

5 files changed, 185 insertions, 115 deletions

diff --git a/Wrappers/Python/ccpi/framework/framework.py b/Wrappers/Python/ccpi/framework/framework.py
index a8a0ab4..387b5c1 100755
--- a/Wrappers/Python/ccpi/framework/framework.py
+++ b/Wrappers/Python/ccpi/framework/framework.py
@@ -760,13 +760,6 @@ class DataContainer(object):
         return numpy.sqrt(self.squared_norm())
     
     
-#    def dot(self, other, *args, **kwargs):
-#        '''return the inner product of 2 DataContainers viewed as vectors'''
-#        if self.shape == other.shape:
-#            return numpy.dot(self.as_array().ravel(), other.as_array().ravel())
-#        else:
-#            raise ValueError('Shapes are not aligned: {} != {}'.format(self.shape, other.shape))
-
     def dot(self, other, *args, **kwargs):
         '''return the inner product of 2 DataContainers viewed as vectors'''
         method = kwargs.get('method', 'reduce')
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index e0b8a32..6ffaf70 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -111,7 +111,7 @@ class Gradient(LinearOperator):
         return BlockDataContainer(*mat)    
 
 
-    def sum_abs_row(self):
+    def sum_abs_col(self):
         
         tmp = self.gm_range.allocate()
         res = self.gm_domain.allocate()
@@ -120,7 +120,7 @@ class Gradient(LinearOperator):
             res += spMat.sum_abs_row()
         return res
     
-    def sum_abs_col(self):
+    def sum_abs_row(self):
         
         tmp = self.gm_range.allocate()
         res = []
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
index 4903c44..3c295f5 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_Poisson.py
@@ -17,6 +17,27 @@
 #   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, x>0  \alpha * ||\nabla x||_{1} + \int x - g * log(x)
+
+             \nabla: Gradient operator 
+             g: Noisy Data with Poisson Noise
+             \alpha: Regularization parameter
+             
+             Method = 0:  K = [ \nabla,
+                                 Identity]
+                                                                    
+             Method = 1:  K = \nabla    
+             
+             
+"""
+
 from ccpi.framework import ImageData, ImageGeometry
 
 import numpy as np 
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
index 4189acb..f5d4ce4 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Denoising_SaltPepper.py
@@ -17,6 +17,27 @@
 #   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, x>0  \alpha * ||\nabla x||_{1} + ||x-g||_{1}
+
+             \nabla: Gradient operator 
+             g: Noisy Data with Salt & Pepper Noise
+             \alpha: Regularization parameter
+             
+             Method = 0:  K = [ \nabla,
+                                 Identity]
+                                                                    
+             Method = 1:  K = \nabla    
+             
+             
+"""
+
 from ccpi.framework import ImageData, ImageGeometry
 
 import numpy as np 
diff --git a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
index 0711e91..75286e5 100644
--- a/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
+++ b/Wrappers/Python/wip/Demos/PDHG_TV_Tomo2D.py
@@ -31,6 +31,25 @@ from ccpi.optimisation.functions import ZeroFunction, KullbackLeibler, \
 
 from ccpi.astra.ops import AstraProjectorSimple
 
+""" 
+
+Total Variation Denoising using PDHG algorithm:
+
+             min_{x} max_{y} < K x, y > + g(x) - f^{*}(y) 
+
+
+Problem:     min_x, x>0  \alpha * ||\nabla x||_{1} + int A x -g log(Ax + \eta)
+
+             \nabla: Gradient operator 
+             
+             A: Projection Matrix
+             g: Noisy sinogram corrupted with Poisson Noise
+             
+             \eta: Background Noise
+             \alpha: Regularization parameter
+ 
+"""
+
 # Create phantom for TV 2D tomography 
 N = 75
 x = np.zeros((N,N))
@@ -70,12 +89,26 @@ f = BlockFunction(f1, f2)
                                       
 g = ZeroFunction()
     
-# Compute operator Norm
-normK = operator.norm()
+diag_precon =  True
+
+if diag_precon:
+    
+    def tau_sigma_precond(operator):
+        
+        tau = 1/operator.sum_abs_row()
+        sigma = 1/ operator.sum_abs_col()
+               
+        return tau, sigma
+
+    tau, sigma = tau_sigma_precond(operator)
+             
+else:
+    # Compute operator Norm
+    normK = operator.norm()
+    # Primal & dual stepsizes
+    sigma = 10
+    tau = 1/(sigma*normK**2)
 
-# Primal & dual stepsizes
-sigma = 1
-tau = 1/(sigma*normK**2)
 
 
 # Setup and run the PDHG algorithm
@@ -84,6 +117,8 @@ pdhg.max_iteration = 2000
 pdhg.update_objective_interval = 50
 pdhg.run(2000)
 
+
+#%%
 plt.figure(figsize=(15,15))
 plt.subplot(3,1,1)
 plt.imshow(data.as_array())
@@ -108,104 +143,104 @@ plt.show()
 
 #%% Check with CVX solution
 
-from ccpi.optimisation.operators import SparseFiniteDiff
-import astra
-import numpy
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-    
-
-    ##Construct problem    
-    u = Variable(N*N)
-    #q = Variable()
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    
-    # create matrix representation for Astra operator
-
-    vol_geom = astra.create_vol_geom(N, N)
-    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
-
-    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
-
-    matrix_id = astra.projector.matrix(proj_id)
-
-    ProjMat = astra.matrix.get(matrix_id)
-    
-    fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) 
-    #constraints = [q>= fidelity, u>=0]
-    constraints = [u>=0]
-        
-    solver = SCS
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj, constraints)
-    result = prob.solve(verbose = True, solver = solver)    
-     
-
-##%% Check with CVX solution
-
-from ccpi.optimisation.operators import SparseFiniteDiff
-
-try:
-    from cvxpy import *
-    cvx_not_installable = True
-except ImportError:
-    cvx_not_installable = False
-
-
-if cvx_not_installable:
-
-    ##Construct problem    
-    u = Variable(ig.shape)
-    
-    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
-    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
-    
-    # Define Total Variation as a regulariser
-    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
-    fidelity = pnorm( u - noisy_data.as_array(),1)
-    
-    # choose solver
-    if 'MOSEK' in installed_solvers():
-        solver = MOSEK
-    else:
-        solver = SCS  
-        
-    obj =  Minimize( regulariser +  fidelity)
-    prob = Problem(obj)
-    result = prob.solve(verbose = True, solver = solver)
-    
-           
-    plt.figure(figsize=(15,15))
-    plt.subplot(3,1,1)
-    plt.imshow(pdhg.get_output().as_array())
-    plt.title('PDHG solution')
-    plt.colorbar()
-    plt.subplot(3,1,2)
-    plt.imshow(np.reshape(u.value, (N, N)))
-    plt.title('CVX solution')
-    plt.colorbar()
-    plt.subplot(3,1,3)
-    plt.imshow(diff_cvx)
-    plt.title('Difference')
-    plt.colorbar()
-    plt.show()    
-    
-    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.legend()
-    plt.title('Middle Line Profiles')
-    plt.show()
-            
-    print('Primal Objective (CVX) {} '.format(obj.value))
-    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file
+#from ccpi.optimisation.operators import SparseFiniteDiff
+#import astra
+#import numpy
+#
+#try:
+#    from cvxpy import *
+#    cvx_not_installable = True
+#except ImportError:
+#    cvx_not_installable = False
+#
+#
+#if cvx_not_installable:
+#    
+#
+#    ##Construct problem    
+#    u = Variable(N*N)
+#    #q = Variable()
+#    
+#    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+#    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+#
+#    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+#    
+#    # create matrix representation for Astra operator
+#
+#    vol_geom = astra.create_vol_geom(N, N)
+#    proj_geom = astra.create_proj_geom('parallel', 1.0, detectors, angles)
+#
+#    proj_id = astra.create_projector('strip', proj_geom, vol_geom)
+#
+#    matrix_id = astra.projector.matrix(proj_id)
+#
+#    ProjMat = astra.matrix.get(matrix_id)
+#    
+#    fidelity = sum( ProjMat * u - noisy_data.as_array().ravel() * log(ProjMat * u)) 
+#    #constraints = [q>= fidelity, u>=0]
+#    constraints = [u>=0]
+#        
+#    solver = SCS
+#    obj =  Minimize( regulariser +  fidelity)
+#    prob = Problem(obj, constraints)
+#    result = prob.solve(verbose = True, solver = solver)    
+#     
+#
+###%% Check with CVX solution
+#
+#from ccpi.optimisation.operators import SparseFiniteDiff
+#
+#try:
+#    from cvxpy import *
+#    cvx_not_installable = True
+#except ImportError:
+#    cvx_not_installable = False
+#
+#
+#if cvx_not_installable:
+#
+#    ##Construct problem    
+#    u = Variable(ig.shape)
+#    
+#    DY = SparseFiniteDiff(ig, direction=0, bnd_cond='Neumann')
+#    DX = SparseFiniteDiff(ig, direction=1, bnd_cond='Neumann')
+#    
+#    # Define Total Variation as a regulariser
+#    regulariser = alpha * sum(norm(vstack([DX.matrix() * vec(u), DY.matrix() * vec(u)]), 2, axis = 0))
+#    fidelity = pnorm( u - noisy_data.as_array(),1)
+#    
+#    # choose solver
+#    if 'MOSEK' in installed_solvers():
+#        solver = MOSEK
+#    else:
+#        solver = SCS  
+#        
+#    obj =  Minimize( regulariser +  fidelity)
+#    prob = Problem(obj)
+#    result = prob.solve(verbose = True, solver = solver)
+#    
+#           
+#    plt.figure(figsize=(15,15))
+#    plt.subplot(3,1,1)
+#    plt.imshow(pdhg.get_output().as_array())
+#    plt.title('PDHG solution')
+#    plt.colorbar()
+#    plt.subplot(3,1,2)
+#    plt.imshow(np.reshape(u.value, (N, N)))
+#    plt.title('CVX solution')
+#    plt.colorbar()
+#    plt.subplot(3,1,3)
+#    plt.imshow(diff_cvx)
+#    plt.title('Difference')
+#    plt.colorbar()
+#    plt.show()    
+#    
+#    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.legend()
+#    plt.title('Middle Line Profiles')
+#    plt.show()
+#            
+#    print('Primal Objective (CVX) {} '.format(obj.value))
+#    print('Primal Objective (PDHG) {} '.format(pdhg.objective[-1][0]))
\ No newline at end of file