add methods for precond pdhg

3 files changed, 195 insertions, 4 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index ec14b8f..87723f0 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -8,9 +8,9 @@ Created on Fri Mar  1 22:50:04 2019
 
 from ccpi.optimisation.operators import Operator, LinearOperator, ScaledOperator
 from ccpi.optimisation.ops import PowerMethodNonsquare
-from ccpi.framework import ImageData, ImageGeometry, BlockGeometry
+from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
 import numpy 
-from ccpi.optimisation.operators import FiniteDiff
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 
 #%%
 
@@ -45,7 +45,6 @@ class Gradient(LinearOperator):
         
         tmp = self.gm_range.allocate()
         
-        
         for i in range(tmp.shape[0]):
             tmp.get_item(i).fill(FiniteDiff(self.gm_domain, direction = self.ind[i], bnd_cond = self.bnd_cond).direct(x))
         return tmp    
@@ -73,6 +72,56 @@ class Gradient(LinearOperator):
     def __rmul__(self, scalar):
         return ScaledOperator(self, scalar) 
     
+    
+    def matrix(self):
+        
+        tmp = self.gm_range.allocate()
+        
+        mat = []
+        for i in range(tmp.shape[0]):
+            
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            mat.append(spMat.matrix())
+    
+        return BlockDataContainer(*mat)    
+
+
+    def sum_abs_row(self):
+        
+        tmp = self.gm_range.allocate()
+        res = self.gm_domain.allocate()
+        for i in range(tmp.shape[0]):
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            res += spMat.sum_abs_row()
+        return res
+    
+    def sum_abs_col(self):
+        
+        tmp = self.gm_range.allocate()
+        res = []
+        for i in range(tmp.shape[0]):
+            spMat = SparseFiniteDiff(self.gm_domain, direction=self.ind[i], bnd_cond=self.bnd_cond)
+            res.append(spMat.sum_abs_col())
+        return BlockDataContainer(*res)
+   
+       
 if __name__ == '__main__':
     
-    pass
+    M, N = 2, 3
+    ig = ImageGeometry(M, N)
+    arr = ig.allocate('random_int' )
+    
+    G_neum = Gradient(ig)
+    
+    d = G_neum.matrix()
+    print(d[1])
+    
+    d1 = G_neum.sum_abs_row()
+    print(d1.as_array())
+    
+    d2 = G_neum.sum_abs_col()
+    print(d2)    
+    
+    
+    
+    
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
new file mode 100644
index 0000000..0fb5efb
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -0,0 +1,140 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+"""
+Created on Tue Apr  2 14:06:15 2019
+
+@author: vaggelis
+"""
+
+import scipy.sparse as sp
+import numpy as np
+from ccpi.framework import ImageData
+
+class SparseFiniteDiff():
+    
+    def __init__(self, gm_domain, gm_range=None, direction=0, bnd_cond = 'Neumann'):
+        
+        super(SparseFiniteDiff, self).__init__() 
+        self.gm_domain = gm_domain
+        self.gm_range = gm_range
+        self.direction = direction
+        self.bnd_cond = bnd_cond
+        
+        if self.gm_range is None:
+            self.gm_range = self.gm_domain
+            
+        self.get_dims = [i for i in gm_domain.shape]  
+        
+        if self.direction + 1 > len(self.gm_domain.shape):
+            raise ValueError('Gradient directions more than geometry domain')         
+            
+    def matrix(self):    
+        
+            i = self.direction 
+            
+            mat = sp.spdiags(np.vstack([-np.ones((1,self.get_dims[i])),np.ones((1,self.get_dims[i]))]), [0,1], self.get_dims[i], self.get_dims[i], format = 'lil')
+
+            if self.bnd_cond == 'Neumann':
+                mat[-1,:] = 0
+            elif self.bnd_cond == 'Periodic':
+                mat[-1,0] = 1    
+                
+            tmpGrad = mat if i == 0 else sp.eye(self.get_dims[0])
+            
+            for j in range(1, self.gm_domain.length):
+
+                tmpGrad = sp.kron(mat, tmpGrad ) if j == i else sp.kron(sp.eye(self.get_dims[j]), tmpGrad ) 
+                
+            return tmpGrad    
+                         
+    def T(self):        
+        return self.matrix().T
+     
+    def direct(self, x):
+        
+        x_asarr = x.as_array()
+        res = np.reshape( self.matrix() * x_asarr.flatten('F'), self.gm_domain.shape, 'F')
+        return type(x)(res)
+    
+    def adjoint(self, x):
+        
+        x_asarr = x.as_array()
+        res = np.reshape( self.matrix().T * x_asarr.flatten('F'), self.gm_domain.shape, 'F')
+        return type(x)(res) 
+    
+    def sum_abs_row(self):
+        
+        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=0), self.gm_domain.shape, 'F')))
+    
+    def sum_abs_col(self):
+        
+        return ImageData(np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F')))
+        
+if __name__ == '__main__':
+    
+    from ccpi.framework import ImageGeometry
+    from ccpi.optimisation.operators import FiniteDiff
+    
+    # 2D
+    M, N= 2, 3
+    ig = ImageGeometry(M, N)
+    arr = ig.allocate('random_int')
+    
+    for i in [0,1]:
+    
+        # Neumann
+        sFD_neum = SparseFiniteDiff(ig, direction=i, bnd_cond='Neumann')
+        G_neum = FiniteDiff(ig, direction=i, bnd_cond='Neumann')
+        
+        # Periodic
+        sFD_per = SparseFiniteDiff(ig, direction=i, bnd_cond='Periodic')
+        G_per = FiniteDiff(ig, direction=i, bnd_cond='Periodic')    
+    
+        u_neum_direct = G_neum.direct(arr)
+        u_neum_sp_direct = sFD_neum.direct(arr)
+        np.testing.assert_array_almost_equal(u_neum_direct.as_array(), u_neum_sp_direct.as_array(), decimal=4)
+        
+        u_neum_adjoint = G_neum.adjoint(arr)
+        u_neum_sp_adjoint = sFD_neum.adjoint(arr)
+        np.testing.assert_array_almost_equal(u_neum_adjoint.as_array(), u_neum_sp_adjoint.as_array(), decimal=4)    
+        
+        u_per_direct = G_neum.direct(arr)
+        u_per_sp_direct = sFD_neum.direct(arr)  
+        np.testing.assert_array_almost_equal(u_per_direct.as_array(), u_per_sp_direct.as_array(), decimal=4)
+        
+        u_per_adjoint = G_per.adjoint(arr)
+        u_per_sp_adjoint = sFD_per.adjoint(arr)
+        np.testing.assert_array_almost_equal(u_per_adjoint.as_array(), u_per_sp_adjoint.as_array(), decimal=4)      
+    
+    # 3D
+    M, N, K = 2, 3, 4
+    ig3D = ImageGeometry(M, N, K)
+    arr3D = ig3D.allocate('random_int')
+    
+    for i in [0,1,2]:
+    
+        # Neumann
+        sFD_neum3D = SparseFiniteDiff(ig3D, direction=i, bnd_cond='Neumann')
+        G_neum3D = FiniteDiff(ig3D, direction=i, bnd_cond='Neumann')
+        
+        # Periodic
+        sFD_per3D = SparseFiniteDiff(ig3D, direction=i, bnd_cond='Periodic')
+        G_per3D = FiniteDiff(ig3D, direction=i, bnd_cond='Periodic')    
+    
+        u_neum_direct3D = G_neum3D.direct(arr3D)
+        u_neum_sp_direct3D = sFD_neum3D.direct(arr3D)
+        np.testing.assert_array_almost_equal(u_neum_direct3D.as_array(), u_neum_sp_direct3D.as_array(), decimal=4)
+        
+        u_neum_adjoint3D = G_neum3D.adjoint(arr3D)
+        u_neum_sp_adjoint3D = sFD_neum3D.adjoint(arr3D)
+        np.testing.assert_array_almost_equal(u_neum_adjoint3D.as_array(), u_neum_sp_adjoint3D.as_array(), decimal=4)    
+        
+        u_per_direct3D = G_neum3D.direct(arr3D)
+        u_per_sp_direct3D = sFD_neum3D.direct(arr3D)  
+        np.testing.assert_array_almost_equal(u_per_direct3D.as_array(), u_per_sp_direct3D.as_array(), decimal=4)
+        
+        u_per_adjoint3D = G_per3D.adjoint(arr3D)
+        u_per_sp_adjoint3D = sFD_per3D.adjoint(arr3D)
+        np.testing.assert_array_almost_equal(u_per_adjoint3D.as_array(), u_per_sp_adjoint3D.as_array(), decimal=4)      
+    
+    
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/operators/__init__.py b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
index 1c09faf..57d89ad 100755
--- a/Wrappers/Python/ccpi/optimisation/operators/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/__init__.py
@@ -17,3 +17,5 @@ from .GradientOperator import Gradient
 from .SymmetrizedGradientOperator import SymmetrizedGradient

 from .IdentityOperator import Identity

 from .ZeroOperator import ZeroOp

+

+from .SparseFiniteDiff import SparseFiniteDiff