changes for Symmetrized Gradient

8 files changed, 155 insertions, 82 deletions

diff --git a/Wrappers/Python/ccpi/framework/BlockGeometry.py b/Wrappers/Python/ccpi/framework/BlockGeometry.py
index 0f43155..5dd6750 100755
--- a/Wrappers/Python/ccpi/framework/BlockGeometry.py
+++ b/Wrappers/Python/ccpi/framework/BlockGeometry.py
@@ -16,17 +16,16 @@ class BlockGeometry(object):
         ''''''

         self.geometries = args

         self.index = 0

-        #shape = kwargs.get('shape', None)

-        #if shape is None:

-        #   shape = (len(args),1)

+

         shape = (len(args),1)

         self.shape = shape

-        #print (self.shape)

+

         n_elements = functools.reduce(lambda x,y: x*y, shape, 1)

         if len(args) != n_elements:

             raise ValueError(

                     'Dimension and size do not match: expected {} got {}'

                     .format(n_elements, len(args)))

+                      

             

     def get_item(self, index):

         '''returns the Geometry in the BlockGeometry located at position index'''

diff --git a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
index d0e27ae..7c6bc8a 100644
--- a/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
+++ b/Wrappers/Python/ccpi/optimisation/algorithms/PDHG.py
@@ -104,8 +104,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
 
     x_old = operator.domain_geometry().allocate()
     y_old = operator.range_geometry().allocate()       
-        
-    
+            
     xbar = x_old
     x_tmp = x_old
     x = x_old
@@ -118,7 +117,9 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
     
     t = time.time()
     
-    objective = []
+    primal = []
+    dual = []
+    pdgap = []
     
     
     for i in range(niter):
@@ -137,11 +138,18 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
         x_old = x
         y_old = y   
                 
-        if i%100==0:
+#        if i%100==0:
+            
+        p1 = f(operator.direct(x)) + g(x)
+        d1 = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y)))
+        pd1 = p1 - d1
+        
+        primal.append(p1)
+        dual.append(d1)
+        pdgap.append(pd1)
+        
             
-            primal = f(operator.direct(x)) + g(x)
-            dual = -(f.convex_conjugate(y) + g(-1*operator.adjoint(y)))
-            print( i, primal, dual, primal-dual)
+#            print( i, primal, dual, primal-dual)
             
 #            plt.imshow(x.as_array())
 #            plt.show()
@@ -149,7 +157,7 @@ def PDHG_old(f, g, operator, tau = None, sigma = None, opt = None, **kwargs):
                          
     t_end = time.time()        
         
-    return x, t_end - t, objective
+    return x, t_end - t, primal, dual, pdgap
 
 
 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 24c4e4b..0faba22 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -6,12 +6,12 @@ Created on Fri Mar  1 22:51:17 2019
 @author: evangelos
 """
 
-from ccpi.optimisation.operators import Operator
+from ccpi.optimisation.operators import LinearOperator
 from ccpi.optimisation.ops import PowerMethodNonsquare
 from ccpi.framework import ImageData, BlockDataContainer
 import numpy as np
 
-class FiniteDiff(Operator):
+class FiniteDiff(LinearOperator):
     
     # Works for Neum/Symmetric &  periodic boundary conditions
     # TODO add central differences???
diff --git a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
index e00de0c..0c267fc 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/GradientOperator.py
@@ -43,8 +43,7 @@ class Gradient(LinearOperator):
         
     def direct(self, x, out=None):
         
-        tmp = self.gm_range.allocate()
-        
+        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    
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
index 1b88cba..d54db9b 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SparseFiniteDiff.py
@@ -70,7 +70,7 @@ class SparseFiniteDiff():
     
     def sum_abs_col(self):
         
-        res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'C'))
+        res = np.array(np.reshape(abs(self.matrix()).sum(axis=1), self.gm_domain.shape, 'F'))
         res[res==0]=1
         return ImageData(res)
         
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index d908e49..ea3ba8f 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -6,59 +6,93 @@ Created on Fri Mar  1 22:53:55 2019
 @author: evangelos
 """
 
-from ccpi.optimisation.operators import Operator
-from ccpi.optimisation.operators import FiniteDiff
+from ccpi.optimisation.operators import Gradient, Operator, LinearOperator, ScaledOperator
 from ccpi.optimisation.ops import PowerMethodNonsquare
-from ccpi.framework import ImageData, DataContainer
-import numpy as np
+from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, BlockDataContainer
+import numpy 
+from ccpi.optimisation.operators import FiniteDiff, SparseFiniteDiff
 
 
-class SymmetrizedGradient(Operator):
+class SymmetrizedGradient(Gradient):
     
-    def __init__(self, gm_domain, gm_range, bnd_cond = 'Neumann', **kwargs):
+    def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
-        super(SymmetrizedGradient, self).__init__() 
+        super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) 
         
-        self.gm_domain = gm_domain # Domain of Grad Operator
-        self.gm_range = gm_range # Range of Grad Operator
-        self.bnd_cond = bnd_cond # Boundary conditions of Finite Differences
-    
-        # Kwargs Default options            
-        self.memopt = kwargs.get('memopt',False)
-        self.correlation = kwargs.get('correlation','Space') 
+        '''
+         Domain of SymGrad is the Range of Gradient
+        '''
+        self.gm_domain = self.gm_range 
+        self.bnd_cond = bnd_cond
+        
+        self.channels = self.gm_range.get_item(0).channels
         
-        #TODO not tested yet, operator norm???
-        self.voxel_size = kwargs.get('voxel_size',[1]*len(gm_domain))  
-                                             
+        if self.correlation=='Space':
+            if self.channels>1:
+                pass
+            else: 
+#                # 2D image ---> Dx v1, Dyv2, Dx
+                tmp = self.gm_domain.geometries + (self.gm_domain.get_item(0),)
+                self.gm_range = BlockGeometry(*tmp )
+                self.ind1 = range(self.gm_domain.get_item(0).length)
+                self.ind2 = range(self.gm_domain.get_item(0).length-1, -1, -1)
+#                self.order = myorder = [0,1,2 3]
+                
+        elif self.correlation=='SpaceChannels':
+            if self.channels>1:
+                pass
+            else:
+                raise ValueError('No channels to correlate')        
+                                                     
         
     def direct(self, x, out=None):
         
-        tmp = np.zeros(self.gm_range)
-        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0])
-        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1])
-        tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) +
-                        FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) )
+#        tmp = numpy.zeros(self.gm_range)
+#        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0])
+#        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1])
+#        tmp[2] = 0.5 * (FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).adjoint(x.as_array()[0]) +
+#                        FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).adjoint(x.as_array()[1]) )
+#        
+#        return type(x)(tmp)
+
+        tmp = [[None]*2]*2
+        for i in range(2):
+            for j in range(2):  
+                tmp[i][j]=FiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).adjoint(x.get_item(j))
+        tmp = numpy.array(tmp)
+        z = 0.5 * (tmp.T + tmp)
         
-        return type(x)(tmp)
-    
+        return BlockDataContainer(z.tolist())
+
     
     def adjoint(self, x, out=None):
+        pass
         
-        tmp = np.zeros(self.gm_domain)
+        res = []
+        for i in range(2):
+            for j in range(2):
+                
+                restmpFiniteDiff(self.gm_domain.get_item(0), direction = i, bnd_cond = self.bnd_cond).direct(x.get_item(j))
+                
+                
         
-        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) +  \
-                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2])
-                 
-        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) +  \
-                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1])                 
-
-        return type(x)(tmp)          
+#        for 
+        
+#        tmp = numpy.zeros(self.gm_domain)
+#        
+#        tmp[0] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[0]) +  \
+#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[2])
+#                 
+#        tmp[1] = FiniteDiff(self.gm_domain[1:], direction = 1, bnd_cond = self.bnd_cond).direct(x.as_array()[2]) +  \
+#                 FiniteDiff(self.gm_domain[1:], direction = 0, bnd_cond = self.bnd_cond).direct(x.as_array()[1])                 
+#
+#        return type(x)(tmp)          
             
     def alloc_domain_dim(self):
-        return ImageData(np.zeros(self.gm_domain))
+        return ImageData(numpy.zeros(self.gm_domain))
     
     def alloc_range_dim(self):
-        return ImageData(np.zeros(self.range_dim))
+        return ImageData(numpy.zeros(self.range_dim))
     
     def domain_dim(self):
         return self.gm_domain
@@ -80,29 +114,58 @@ if __name__ == '__main__':
     
     ###########################################################################  
     ## Symmetrized Gradient
+    from ccpi.framework import DataContainer
+    from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
+    import numpy as np
     
     N, M = 2, 3
-    ig = (N,M)
-    ig2 = (2,) + ig
-    ig3 = (3,) + ig
-    u1 = DataContainer(np.random.randint(10, size=ig2))
-    w1 = DataContainer(np.random.randint(10, size=ig3))
+    K = 2
+    
+    ig1 = ImageGeometry(N, M)
+    ig2 = ImageGeometry(N, M, channels=K)
+    
+    E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann')
+    E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Periodic')
     
-    E = SymmetrizedGradient(ig2,ig3)
+    print(E1.domain_geometry().shape)
+    print(E2.domain_geometry().shape)
     
-    d1 = E.direct(u1)
-    d2 = E.adjoint(w1)
+    u1 = E1.gm_domain.allocate('random_int')
+    u2 = E2.gm_domain.allocate('random_int')
+        
+    
+    res = E1.direct(u1) 
+    
+    Dx = FiniteDiff(ig1, direction = 1, bnd_cond = 'Neumann')
+    Dy = FiniteDiff(ig1, direction = 0, bnd_cond = 'Neumann')
+    
+    B = BlockOperator(Dy, Dx)
+    V = BlockDataContainer(u1,u2)
+    
+    res = B.adjoint(V)
     
-    LHS = (d1.as_array()[0]*w1.as_array()[0] + \
-           d1.as_array()[1]*w1.as_array()[1] + \
-           2*d1.as_array()[2]*w1.as_array()[2]).sum()
+#    ig = (N,M)
+#    ig2 = (2,) + ig
+#    ig3 = (3,) + ig
+#    u1 = ig.allocate('random_int')
+#    w1 = E.gm_range.allocate('random_int')
+#    DataContainer(np.random.randint(10, size=ig3))
     
-    RHS = (u1.as_array()[0]*d2.as_array()[0] + \
-           u1.as_array()[1]*d2.as_array()[1]).sum()   
     
     
-    print(LHS, RHS, E.norm())
+#    d1 = E.direct(u1)
+#    d2 = E.adjoint(w1)
     
+#    LHS = (d1.as_array()[0]*w1.as_array()[0] + \
+#           d1.as_array()[1]*w1.as_array()[1] + \
+#           2*d1.as_array()[2]*w1.as_array()[2]).sum()
+#    
+#    RHS = (u1.as_array()[0]*d2.as_array()[0] + \
+#           u1.as_array()[1]*d2.as_array()[1]).sum()   
+#    
+#    
+#    print(LHS, RHS, E.norm())
+#    
     
 #    
     
diff --git a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
index d5c021d..426ce8b 100644
--- a/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
+++ b/Wrappers/Python/wip/pdhg_TV_denoising_precond.py
@@ -74,7 +74,7 @@ else:
     ###########################################################################
 #%%
 
-diag_precon = True 
+diag_precon =  False
 
 if diag_precon:
     
@@ -82,7 +82,7 @@ if diag_precon:
         
         tau = 1/operator.sum_abs_row()
         sigma = 1/ operator.sum_abs_col()
-    
+               
         return tau, sigma
 
     tau, sigma = tau_sigma_precond(operator)
@@ -99,14 +99,19 @@ else:
 #%%
     
 opt = {'niter':2000}
-#
-res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
  
-aaa = res[0].as_array()
-#    
-plt.imshow(aaa)
+plt.figure(figsize=(5,5))
+plt.imshow(res.as_array())
 plt.colorbar()
 plt.show()
+
+#aaa = res[0].as_array()
+#    
+#plt.imshow(aaa)
+#plt.colorbar()
+#plt.show()
 #c2 = aaa
 #del aaa
 #%%
@@ -114,9 +119,9 @@ plt.show()
 #c2 = aaa
 ##%%    
 #%%
-z = c1 - c2
-plt.imshow(np.abs(z[0:95,0:95]))
-plt.colorbar()
+#z = c1 - c2
+#plt.imshow(np.abs(z[0:95,0:95]))
+#plt.colorbar()
 
 #%%
 #pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
diff --git a/Wrappers/Python/wip/pdhg_TV_tomography2D.py b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
index 9feb05b..f06f166 100644
--- a/Wrappers/Python/wip/pdhg_TV_tomography2D.py
+++ b/Wrappers/Python/wip/pdhg_TV_tomography2D.py
@@ -13,7 +13,7 @@ from ccpi.framework import ImageData, ImageGeometry, BlockDataContainer, Acquisi
 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, \
@@ -109,7 +109,7 @@ if diag_precon:
     tau, sigma = tau_sigma_precond(operator)
     
 else:    
-    sigma = 10
+    sigma = 1
     tau = 1/(sigma*normK**2)
 
 #pdhg = PDHG(f=f,g=g,operator=operator, tau=tau, sigma=sigma)
@@ -120,13 +120,12 @@ else:
     
 opt = {'niter':2000}
 #
-res = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
+res, time, primal, dual, pdgap = PDHG_old(f, g, operator, tau = tau, sigma = sigma, opt = opt) 
 
-aaa = res[0].as_array()
-#    
-plt.imshow(aaa)
+plt.figure(figsize=(5,5))
+plt.imshow(res.as_array())
 plt.colorbar()
-plt.show()    
+plt.show() 
 
 #%%
 #sol = pdhg.get_output().as_array()