Fix sym grad (#393)

* fix func style

* fix symGrad

* recover from master

* fix L1

* fix L1

* Update KullbackLeibler.py

* correct adjoint FD

* fix symgrad

* no change in adjoint FD

* add comment on FD

* added unittest

* changed syntax for assert_almost_equal

* split test of sym gradient

* fixed test

6 files changed, 160 insertions, 112 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index d71f597..a6c9545 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -240,4 +240,4 @@ if __name__ == '__main__':
         
 
     
-        
+    
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
index 09e550e..1c2c43f 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
@@ -177,4 +177,4 @@ if __name__ == '__main__':
         
         
         
-      
+      
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
index 92e0116..f5108ba 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
@@ -280,4 +280,4 @@ if __name__ == '__main__':
     
     numpy.testing.assert_array_almost_equal(res1.as_array(), \
                                             res2.as_array(), decimal=4)
-                                            
+                                            
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
index 6e2f49a..0dd7d4b 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/FiniteDifferenceOperator.py
@@ -194,7 +194,7 @@ class FiniteDiff(LinearOperator):
 
                     
     def adjoint(self, x, out=None):
-        
+                      
         x_asarr = x.as_array()
         x_sz = len(x.shape)
         outnone = False 
diff --git a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
index 8d14cf8..c85abfa 100644
--- a/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
+++ b/Wrappers/Python/ccpi/optimisation/operators/SymmetrizedGradientOperator.py
@@ -21,15 +21,12 @@ from __future__ import print_function
 from __future__ import unicode_literals
 
 
-from ccpi.optimisation.operators import Gradient, Operator, LinearOperator,\
-        ScaledOperator
-from ccpi.framework import ImageData, ImageGeometry, BlockGeometry, \
-        BlockDataContainer
-import numpy 
+from ccpi.optimisation.operators import LinearOperator
+from ccpi.framework import BlockGeometry, BlockDataContainer
 from ccpi.optimisation.operators import FiniteDiff
 
 
-class SymmetrizedGradient(Gradient):
+class SymmetrizedGradient(LinearOperator):
     
     r'''Symmetrized Gradient Operator:  E: V -> W
         
@@ -49,21 +46,23 @@ class SymmetrizedGradient(Gradient):
               |                                                      
     '''
     
+    CORRELATION_SPACE = "Space"
+    CORRELATION_SPACECHANNEL = "SpaceChannels"
     
     def __init__(self, gm_domain, bnd_cond = 'Neumann', **kwargs):
         
-        super(SymmetrizedGradient, self).__init__(gm_domain, bnd_cond, **kwargs) 
+        super(SymmetrizedGradient, self).__init__() 
                 
-        self.gm_domain = self.gm_range 
+        self.gm_domain = gm_domain
         self.bnd_cond = bnd_cond
-        
-        self.channels = self.gm_range.get_item(0).channels
-        
+        self.correlation = kwargs.get('correlation',SymmetrizedGradient.CORRELATION_SPACE)
+                
         tmp_gm = len(self.gm_domain.geometries)*self.gm_domain.geometries
         
         self.gm_range = BlockGeometry(*tmp_gm)
         
-        self.FD = FiniteDiff(self.gm_domain, direction = 0, bnd_cond = self.bnd_cond)
+        # Define FD operator. We need one geometry from the BlockGeometry of the domain
+        self.FD = FiniteDiff(self.gm_domain.get_item(0), direction = 0, bnd_cond = self.bnd_cond)
         
         if self.gm_domain.shape[0]==2:
             self.order_ind = [0,2,1,3]
@@ -130,8 +129,6 @@ class SymmetrizedGradient(Gradient):
                     i+=1
                     tmp1+=tmp[j]
                 out[k].fill(tmp1)
-#            tmp = self.adjoint(x)
-#            out.fill(tmp)
                     
                      
     def domain_geometry(self):
@@ -144,110 +141,59 @@ if __name__ == '__main__':
     
     ###########################################################################  
     ## Symmetrized Gradient Tests
-    from ccpi.framework import DataContainer
-    from ccpi.optimisation.operators import Gradient, BlockOperator, FiniteDiff
+    from ccpi.framework import ImageGeometry
+    from ccpi.optimisation.operators import Gradient
     import numpy as np
     
     N, M = 2, 3
     K = 2
     C = 2
     
-    ig1 = ImageGeometry(N, M)
-    ig2 = ImageGeometry(N, M, channels=C)
+    ###########################################################################
+    # 2D geometry no channels
+    ig = ImageGeometry(N, M)
+    Grad = Gradient(ig)
     
-    E1 = SymmetrizedGradient(ig1, correlation = 'Space', bnd_cond='Neumann')
+    E1 = SymmetrizedGradient(Grad.range_geometry())
+    np.testing.assert_almost_equal(E1.norm(), np.sqrt(8), 1e-5)
     
-    try:
-        E1 = SymmetrizedGradient(ig1, correlation = 'SpaceChannels', bnd_cond='Neumann')
-    except:
-        print("No Channels to correlate")
-        
-    E2 = SymmetrizedGradient(ig2, correlation = 'SpaceChannels', bnd_cond='Neumann')  
-
     print(E1.domain_geometry().shape, E1.range_geometry().shape)
-    print(E2.domain_geometry().shape, E2.range_geometry().shape)   
-    
-    #check Linear operator property
-    
     u1 = E1.domain_geometry().allocate('random_int')
-    u2 = E2.domain_geometry().allocate('random_int')
-    
-    # Need to allocate random_int at the Range of SymGradient
-    
-    #a1 = ig1.allocate('random_int')
-    #a2 = ig1.allocate('random_int')
-    #a3 = ig1.allocate('random_int')
-    
-    #a4 = ig1.allocate('random_int')
-    #a5 = ig1.allocate('random_int')    
-    #a6 = ig1.allocate('random_int')
-    
-    # TODO allocate has to create this symmetry by default!!!!!
-    #w1 = BlockDataContainer(*[a1, a2, \
-    #                           a2, a3]) 
-    w1 = E1.range_geometry().allocate('random_int',symmetry=True)
-    
-    LHS = (E1.direct(u1) * w1).sum()
-    RHS = (u1 * E1.adjoint(w1)).sum()
-    
-    numpy.testing.assert_equal(LHS, RHS)     
+    w1 = E1.range_geometry().allocate('random_int', symmetry = True)
+       
     
-    u2 = E2.gm_domain.allocate('random_int')
+    lhs = E1.direct(u1).dot(w1)
+    rhs = u1.dot(E1.adjoint(w1))
+    np.testing.assert_almost_equal(lhs, rhs)
     
-    #aa1 = ig2.allocate('random_int')
-    #aa2 = ig2.allocate('random_int')
-    #aa3 = ig2.allocate('random_int')
-    #aa4 = ig2.allocate('random_int')
-    #aa5 = ig2.allocate('random_int')
-    #aa6 = ig2.allocate('random_int')  
+    ###########################################################################
+    # 2D geometry with channels
+    ig2 = ImageGeometry(N, M, channels = C)
+    Grad2 = Gradient(ig2, correlation = 'Space')
     
-    #w2 = BlockDataContainer(*[aa1, aa2, aa3, \
-    #                          aa2, aa4, aa5, \
-    #                          aa3, aa5, aa6])     
-    w2 = E2.range_geometry().allocate('random_int',symmetry=True)
- 
+    E2 = SymmetrizedGradient(Grad2.range_geometry())
+    np.testing.assert_almost_equal(E2.norm(), np.sqrt(12), 1e-6)
     
-    LHS1 = (E2.direct(u2) * w2).sum()
-    RHS1 = (u2 * E2.adjoint(w2)).sum()
-    
-    numpy.testing.assert_equal(LHS1, RHS1)      
-    
-    out = E1.range_geometry().allocate()
-    E1.direct(u1, out=out)
-    a1 = E1.direct(u1)
-    numpy.testing.assert_array_equal(a1[0].as_array(), out[0].as_array()) 
-    numpy.testing.assert_array_equal(a1[1].as_array(), out[1].as_array()) 
-    numpy.testing.assert_array_equal(a1[2].as_array(), out[2].as_array()) 
-    numpy.testing.assert_array_equal(a1[3].as_array(), out[3].as_array()) 
-    
-    
-    out1 = E1.domain_geometry().allocate()
-    E1.adjoint(w1, out=out1)
-    b1 = E1.adjoint(w1)    
-    
-    LHS_out = (out * w1).sum()
-    RHS_out = (u1 * out1).sum()
-    print(LHS_out, RHS_out)
-    
-    
-    out2 = E2.range_geometry().allocate()
-    E2.direct(u2, out=out2)
-    a2 = E2.direct(u2)
-    
-    out21 = E2.domain_geometry().allocate()
-    E2.adjoint(w2, out=out21)
-    b2 = E2.adjoint(w2)    
-    
-    LHS_out = (out2 * w2).sum()
-    RHS_out = (u2 * out21).sum()
-    print(LHS_out, RHS_out)    
-    
-    
-    out = E1.range_geometry().allocate()
-    E1.direct(u1, out=out)
-    E1.adjoint(out, out=out1)
-    
-    print(E1.norm())
-    print(E2.norm())
-    
-  
+    print(E2.domain_geometry().shape, E2.range_geometry().shape)
+    u2 = E2.domain_geometry().allocate('random_int')
+    w2 = E2.range_geometry().allocate('random_int', symmetry = True)
+#    
+    lhs2 = E2.direct(u2).dot(w2)
+    rhs2 = u2.dot(E2.adjoint(w2))
+    np.testing.assert_almost_equal(lhs2, rhs2)
+    
+    ###########################################################################
+    # 3D geometry no channels
+    ig3 = ImageGeometry(N, M, K)
+    Grad3 = Gradient(ig3, correlation = 'Space')
+    
+    E3 = SymmetrizedGradient(Grad3.range_geometry())
+    np.testing.assert_almost_equal(E3.norm(), np.sqrt(12), 1e-6)
+    
+    print(E3.domain_geometry().shape, E3.range_geometry().shape)
+    u3 = E3.domain_geometry().allocate('random_int')
+    w3 = E3.range_geometry().allocate('random_int', symmetry = True)
+#    
+    lhs3 = E3.direct(u3).dot(w3)
+    rhs3 = u3.dot(E3.adjoint(w3))
+    np.testing.assert_almost_equal(lhs3, rhs3)  
\ No newline at end of file
diff --git a/Wrappers/Python/test/test_Operator.py b/Wrappers/Python/test/test_Operator.py
index b26bb5d..3372b9b 100644
--- a/Wrappers/Python/test/test_Operator.py
+++ b/Wrappers/Python/test/test_Operator.py
@@ -18,13 +18,14 @@
 import unittest
 from ccpi.framework import ImageGeometry, ImageData, BlockDataContainer, DataContainer
 from ccpi.optimisation.operators import BlockOperator, BlockScaledOperator,\
-    FiniteDiff
+    FiniteDiff, SymmetrizedGradient
 import numpy
 from timeit import default_timer as timer
 from ccpi.framework import ImageGeometry
 from ccpi.optimisation.operators import Gradient, Identity, SparseFiniteDiff
 from ccpi.optimisation.operators import LinearOperator
 
+
 def dt(steps):
     return steps[-1] - steps[-2]
 
@@ -170,6 +171,107 @@ class TestOperator(CCPiTestClass):
         print ("Norm dT1 {} dT2 {}".format(t1-t0,t2-t1))
         self.assertLess(t2-t1, t1-t0)
 
+class TestGradients(CCPiTestClass): 
+    def setUp(self):
+        N, M = 20, 30
+        K = 20
+        C = 20
+        self.decimal = 1
+        self.iterations = 50
+        ###########################################################################
+        # 2D geometry no channels
+        self.ig = ImageGeometry(N, M)
+        self.ig2 = ImageGeometry(N, M, channels = C)
+        self.ig3 = ImageGeometry(N, M, K)
+
+    def test_SymmetrizedGradient1a(self):
+        ###########################################################################  
+        ## Symmetrized Gradient Tests
+        print ("Test SymmetrizedGradient")
+        ###########################################################################
+        # 2D geometry no channels
+        # ig = ImageGeometry(N, M)
+        Grad = Gradient(self.ig)
+        
+        E1 = SymmetrizedGradient(Grad.range_geometry())
+        numpy.testing.assert_almost_equal(E1.norm(iterations=self.iterations), numpy.sqrt(8), decimal = self.decimal)
+        
+    def test_SymmetrizedGradient1b(self):
+        ###########################################################################  
+        ## Symmetrized Gradient Tests
+        print ("Test SymmetrizedGradient")
+        ###########################################################################
+        # 2D geometry no channels
+        # ig = ImageGeometry(N, M)
+        Grad = Gradient(self.ig)
+        
+        E1 = SymmetrizedGradient(Grad.range_geometry())
+        u1 = E1.domain_geometry().allocate('random_int')
+        w1 = E1.range_geometry().allocate('random_int', symmetry = True)
+        
+        
+        lhs = E1.direct(u1).dot(w1)
+        rhs = u1.dot(E1.adjoint(w1))
+        # self.assertAlmostEqual(lhs, rhs)
+        numpy.testing.assert_almost_equal(lhs, rhs)
+            
+    def test_SymmetrizedGradient2(self):        
+        ###########################################################################
+        # 2D geometry with channels
+        # ig2 = ImageGeometry(N, M, channels = C)
+        Grad2 = Gradient(self.ig2, correlation = 'Space')
+        
+        E2 = SymmetrizedGradient(Grad2.range_geometry())
+        
+        u2 = E2.domain_geometry().allocate('random_int')
+        w2 = E2.range_geometry().allocate('random_int', symmetry = True)
+    #    
+        lhs2 = E2.direct(u2).dot(w2)
+        rhs2 = u2.dot(E2.adjoint(w2))
+            
+        numpy.testing.assert_almost_equal(lhs2, rhs2)
+        
+    def test_SymmetrizedGradient2a(self):        
+        ###########################################################################
+        # 2D geometry with channels
+        # ig2 = ImageGeometry(N, M, channels = C)
+        Grad2 = Gradient(self.ig2, correlation = 'Space')
+        
+        E2 = SymmetrizedGradient(Grad2.range_geometry())
+        numpy.testing.assert_almost_equal(E2.norm(iterations=self.iterations), 
+           numpy.sqrt(8), decimal = self.decimal)
+        
+    
+    def test_SymmetrizedGradient3a(self):
+        ###########################################################################
+        # 3D geometry no channels
+        #ig3 = ImageGeometry(N, M, K)
+        Grad3 = Gradient(self.ig3, correlation = 'Space')
+        
+        E3 = SymmetrizedGradient(Grad3.range_geometry())
+
+        norm1 = E3.norm()
+        norm2 = E3.calculate_norm(iterations=100)
+        print (norm1,norm2)
+        numpy.testing.assert_almost_equal(norm2, numpy.sqrt(12), decimal = self.decimal)
+        
+    def test_SymmetrizedGradient3b(self):
+        ###########################################################################
+        # 3D geometry no channels
+        #ig3 = ImageGeometry(N, M, K)
+        Grad3 = Gradient(self.ig3, correlation = 'Space')
+        
+        E3 = SymmetrizedGradient(Grad3.range_geometry())
+        
+        u3 = E3.domain_geometry().allocate('random_int')
+        w3 = E3.range_geometry().allocate('random_int', symmetry = True)
+    #    
+        lhs3 = E3.direct(u3).dot(w3)
+        rhs3 = u3.dot(E3.adjoint(w3))
+        numpy.testing.assert_almost_equal(lhs3, rhs3)  
+        # self.assertAlmostEqual(lhs3, rhs3)
+        # self.assertTrue( LinearOperator.dot_test(E3) )
+