refactoring add unittest

8 files changed, 67 insertions, 429 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
index bc9b62d..d6c98c4 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/BlockFunction.py
@@ -7,7 +7,8 @@ Created on Fri Mar  8 10:01:31 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import BlockDataContainer
 
 class BlockFunction(Function):
diff --git a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
index a67b884..0f3defe 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/FunctionOperatorComposition.py
@@ -7,7 +7,8 @@ Created on Fri Mar  8 09:55:36 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 
 
 class FunctionOperatorComposition(Function):
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
index ec7aa5b..f83de6f 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L1Norm.py
@@ -7,12 +7,13 @@ Created on Wed Mar  6 19:42:34 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import DataContainer, ImageData, ImageGeometry 
 
 
 ############################   L1NORM FUNCTIONS   #############################
-class SimpleL1NormEdo(Function):
+class SimpleL1Norm(Function):
     
     def __init__(self, alpha=1):
         
@@ -35,7 +36,7 @@ class SimpleL1NormEdo(Function):
     def proximal_conjugate(self, x, tau):
         return x.divide((x.abs()/self.alpha).maximum(1.0))
     
-class L1Norm(SimpleL1NormEdo):
+class L1Norm(SimpleL1Norm):
     
     def __init__(self, alpha=1, **kwargs):
         
diff --git a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
index 8b73be6..5817317 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/L2NormSquared.py
@@ -9,7 +9,8 @@ Created on Thu Feb  7 13:10:56 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import DataContainer, ImageData, ImageGeometry 
 
     
diff --git a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
index b41cc26..9def741 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/ZeroFun.py
@@ -7,7 +7,8 @@ Created on Wed Mar  6 19:44:10 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import DataContainer, ImageData
 from ccpi.framework import BlockDataContainer 
 
@@ -33,8 +34,11 @@ class ZeroFun(Function):
             else:
                 return x.maximum(0).sum() + x.maximum(0).sum()
     
-    def proximal(self,x,tau):
-        return x.copy()
+    def proximal(self,x,tau, out=None):
+        if out is None:
+            return x.copy()
+        else:
+            out.fill(x)
         
     def proximal_conjugate(self, x, tau):
         return 0
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
index 7ce617a..c4ba0a6 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
@@ -1,6 +1,6 @@
 # -*- coding: utf-8 -*-
 
-
+from .Function import Function
 from .ZeroFun import ZeroFun
 from .L1Norm import *
 from .L2NormSquared import *
diff --git a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py
index 0ce1d28..8fe8620 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/mixed_L12Norm.py
@@ -7,9 +7,10 @@ Created on Wed Mar  6 19:43:12 2019
 """
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import DataContainer, ImageData, ImageGeometry 
-
+from ccpi.optimisation.functions.FunctionOperatorComposition import FunctionOperatorComposition
 
 
 ############################   mixed_L1,2NORM FUNCTIONS   #############################
@@ -59,7 +60,7 @@ class mixed_L12Norm(Function):
     def composition_with(self, operator):
         
         if self.b is None:
-            return FunctionComposition(mixed_L12Norm(self.alpha), operator)
+            return FunctionOperatorComposition(operator, mixed_L12Norm(self.alpha))
         else:
-            return FunctionComposition(mixed_L12Norm(self.alpha, b=self.b), operator)    
+            return FunctionOperatorComposition(operator, mixed_L12Norm(self.alpha, b=self.b))
     
diff --git a/Wrappers/Python/test/test_functions.py b/Wrappers/Python/test/test_functions.py
index 54a952a..0741d1c 100644
--- a/Wrappers/Python/test/test_functions.py
+++ b/Wrappers/Python/test/test_functions.py
@@ -6,438 +6,67 @@ Created on Sat Mar  2 19:24:37 2019
 @author: evangelos
 """
 
-# -*- coding: utf-8 -*-
-
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Feb  7 13:10:56 2019
-
-@author: evangelos
-"""
 
 import numpy as np
-from ccpi.optimisation.funcs import Function
+#from ccpi.optimisation.funcs import Function
+from ccpi.optimisation.functions import Function
 from ccpi.framework import DataContainer, ImageData, ImageGeometry 
-from operators import CompositeDataContainer, Identity, CompositeOperator
+from ccpi.optimisation.operators import  Identity
+from ccpi.optimisation.operators import BlockOperator
+from ccpi.framework import BlockDataContainer
 from numbers import Number
-from GradientOperator import Gradient
-
-    
-class SimpleL2NormSq(Function):
-    
-    def __init__(self, alpha=1):
-        
-        super(SimpleL2NormSq, self).__init__()         
-        self.alpha = alpha
-        
-    def __call__(self, x):
-        return self.alpha * x.power(2).sum()
-    
-    def gradient(self,x):
-        return 2 * self.alpha * x
-    
-    def convex_conjugate(self,x):
-        return (1/(4*self.alpha)) * x.power(2).sum()
-    
-    def proximal(self, x, tau):
-        return x.divide(1+2*tau*self.alpha)
-    
-    def proximal_conjugate(self, x, tau):
-        return x.divide(1 + tau/(2*self.alpha) )    
-
-
-############################   L2NORM FUNCTIONS   #############################
-class L2NormSq(SimpleL2NormSq):
-    
-    def __init__(self, alpha, **kwargs):
-                 
-        super(L2NormSq, self).__init__(alpha)
-        self.alpha = alpha
-        self.b = kwargs.get('b',None)
-        self.L = 1        
-        
-    def __call__(self, x):
-        
-        if self.b is None:
-            return SimpleL2NormSq.__call__(self, x)
-        else:
-            return SimpleL2NormSq.__call__(self, x - self.b) 
-        
-    def gradient(self, x):
-        
-        if self.b is None:
-            return 2*self.alpha * x 
-        else:
-            return 2*self.alpha * (x - self.b) 
-                    
-    def composition_with(self, operator):        
-        
-        if self.b is None:
-            return FunctionComposition(L2NormSq(self.alpha), operator)
-        else:
-            return FunctionComposition(L2NormSq(self.alpha, b=self.b), operator)
-                                             
-    def convex_conjugate(self, x):
-        
-        ''' The convex conjugate corresponds to the simple functional i.e., 
-        f(x) = alpha * ||x - b||_{2}^{2}
-        '''
-        
-        if self.b is None:
-            return SimpleL2NormSq.convex_conjugate(self, x)
-        else:
-            return SimpleL2NormSq.convex_conjugate(self, x) + (self.b * x).sum()
-    
-    def proximal(self, x, tau):
-        
-        ''' The proximal operator corresponds to the simple functional i.e., 
-        f(x) = alpha * ||x - b||_{2}^{2}
-        
-        argmin_x { 0.5||x - u||^{2} + tau f(x) }
-        '''
-        
-        if self.b is None:
-            return SimpleL2NormSq.proximal(self, x, tau)
-        else:
-            return self.b + SimpleL2NormSq.proximal(self, x - self.b , tau)
-            
-    
-    def proximal_conjugate(self, x, tau):
-        
-        ''' The proximal operator corresponds to the simple convex conjugate 
-        functional i.e., f^{*}(x^{)        
-        argmin_x { 0.5||x - u||^{2} + tau f(x) }
-        '''
-        if self.b is None:
-            return SimpleL2NormSq.proximal_conjugate(self, x, tau)
-        else:
-            return SimpleL2NormSq.proximal_conjugate(self, x - tau * self.b, tau)
-     
-        
-############################   L1NORM FUNCTIONS   #############################
-class SimpleL1Norm(Function):
-    
-    def __init__(self, alpha=1):
-        
-        super(SimpleL1Norm, self).__init__()         
-        self.alpha = alpha
-        
-    def __call__(self, x):
-        return self.alpha * x.abs().sum()
-    
-    def gradient(self,x):
-        return ValueError('Not Differentiable')
-            
-    def convex_conjugate(self,x):
-        return 0
-    
-    def proximal(self, x, tau):
-        ''' Soft Threshold'''
-        return x.sign() * (x.abs() - tau * self.alpha).maximum(1.0)
-        
-    def proximal_conjugate(self, x, tau):
-        return x.divide((x.abs()/self.alpha).maximum(1.0))
-    
-class L1Norm(SimpleL1Norm):
-    
-    def __init__(self, alpha=1, **kwargs):
-        
-        super(L1Norm, self).__init__()         
-        self.alpha = alpha 
-        
-        self.A = kwargs.get('A',None)
-        self.b = kwargs.get('b',None)
-        
-    def __call__(self, x):
-        
-        if self.b is None:
-            return SimpleL1Norm.__call__(self, self.A.direct(x))
-        else:
-            return SimpleL1Norm.__call__(self, self.A.direct(x) - self.b)
-        
-    def eval_norm(self, x):
-        
-        return SimpleL1Norm.__call__(self, x)        
-    
-    def gradient(self, x):
-        return ValueError('Not Differentiable')
-            
-    def convex_conjugate(self,x):
-        if self.b is None:
-            return SimpleL1Norm.convex_conjugate(self, x)
-        else:
-            return SimpleL1Norm.convex_conjugate(self, x) + (self.b * x).sum()
-    
-    def proximal(self, x, tau):
-        
-        if self.b is None:
-            return SimpleL1Norm.proximal(self, x, tau)
-        else:
-            return self.b + SimpleL1Norm.proximal(self, x + self.b , tau)
-        
-    def proximal_conjugate(self, x, tau):
-        
-        if self.b is None:
-            return SimpleL1Norm.proximal_conjugate(self, x, tau)
-        else:
-            return SimpleL1Norm.proximal_conjugate(self, x - tau*self.b, tau)
-                        
-
-############################   mixed_L1,2NORM FUNCTIONS   #############################
-class mixed_L12Norm(Function):
-    
-    def __init__(self, alpha, **kwargs):
+from ccpi.optimisation.operators import Gradient
 
-        super(mixed_L12Norm, self).__init__() 
-        
-        self.alpha = alpha 
-        self.b = kwargs.get('b',None)                
-        self.sym_grad = kwargs.get('sym_grad',False)
-        
-    def __call__(self,x):
-        
-        if self.b is None:
-            tmp1 = x
-        else:
-            tmp1 = x - self.b            
-#        
-        if self.sym_grad:
-            tmp = np.sqrt(tmp1.as_array()[0]**2 +  tmp1.as_array()[1]**2 +  2*tmp1.as_array()[2]**2)
-        else:
-            tmp = ImageData(tmp1.power(2).sum(axis=0)).sqrt()
-            
-        return self.alpha*tmp.sum()          
-                            
-    def gradient(self,x):
-        return ValueError('Not Differentiable')
-                            
-    def convex_conjugate(self,x):
-        return 0
-    
-    def proximal(self, x, tau):
-        pass
-    
-    def proximal_conjugate(self, x, tau): 
-        
-        if self.sym_grad:
-                tmp2 = np.sqrt(x.as_array()[0]**2 +  x.as_array()[1]**2 +  2*x.as_array()[2]**2)/self.alpha
-                res = x.divide(ImageData(tmp2).maximum(1.0))                                
-        else:
-                res = x.divide((ImageData(x.power(2).sum(axis=0)).sqrt()/self.alpha).maximum(1.0))  
-                                                   
-        return res 
-    
-    def composition_with(self, operator):
-        
-        if self.b is None:
-            return FunctionComposition(mixed_L12Norm(self.alpha), operator)
-        else:
-            return FunctionComposition(mixed_L12Norm(self.alpha, b=self.b), operator)    
-    
+from ccpi.optimisation.functions import SimpleL2NormSq
+from ccpi.optimisation.functions import L2NormSq
+from ccpi.optimisation.functions import SimpleL1Norm
+from ccpi.optimisation.functions import L1Norm
+# from ccpi.optimisation.functions.L2NormSquared import SimpleL2NormSq, L2NormSq
+# from ccpi.optimisation.functions.L1Norm import SimpleL1Norm, L1Norm
+from ccpi.optimisation.functions import mixed_L12Norm
+from ccpi.optimisation.functions import ZeroFun
 
-#%%
-                  
-class ZeroFun(Function):
-    
-    def __init__(self):
-        super(ZeroFun, self).__init__()
-              
-    def __call__(self,x):
-        return 0
-    
-    def convex_conjugate(self, x):
-        ''' This is the support function sup <x, x^{*}>  which in fact is the 
-        indicator function for the set = {0}
-        So 0 if x=0, or inf if x neq 0                
-        '''
-        
-        if x.shape[0]==1:
-            return x.maximum(0).sum()
-        else:
-            return x.get_item(0).maximum(0).sum() + x.get_item(1).maximum(0).sum()
-    
-    def proximal(self,x,tau):
-        return x.copy()
-        
-    def proximal_conjugate(self, x, tau):
-        return 0
-            
-        
-class CompositeFunction(Function):
-    
-    def __init__(self, *functions, blockMatrix):
-
-        self.blockMatrix = blockMatrix        
-        self.functions = functions
-        self.length = len(self.functions)
-        
-    def get_item(self, ind):        
-        return self.functions[ind]        
-                
-    def __call__(self,x):
-                                
-        t = 0                
-        for i in range(x.shape[0]):
-            t += self.functions[i](x.get_item(i))               
-        return t            
-    
-
-    def convex_conjugate(self, x):
-        
-        z = 0
-        t = 0
-        for i in range(x.shape[0]):
-            t += self.functions[z].convex_conjugate(x.get_item(i))
-            z += 1        
-
-        return t 
-    
-    def proximal(self, x, tau, out = None):
-        
-        if isinstance(tau, Number):
-            tau = CompositeDataContainer(tau)
-        out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal(x.get_item(i), tau.get_item(0))
-        return CompositeDataContainer(*out) 
+from ccpi.optimisation.functions import FunctionOperatorComposition
+import unittest        
 
-                            
-    def proximal_conjugate(self, x, tau, out = None):
-        
-        if isinstance(tau, Number):
-            tau = CompositeDataContainer(tau)
-        if isinstance(x, ImageData):
-            x = CompositeDataContainer(x)
-        out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau.get_item(0))
-        return CompositeDataContainer(*out)    
-    
 #    
-class  FunctionComposition(Function):
-    
-    def __init__(self, function, operator):
-        
-        self.function = function
-        self.alpha = self.function.alpha
-        self.b  = self.function.b         
-        self.operator = operator
-        
-    
-        super(FunctionComposition, self).__init__()
-    
-    '''    overide call and gradient '''
-    def __call__(self, x):        
-        return self.function(x)
-    
-    def gradient(self,x):        
-        return self.operator.adjoint(self.function.gradient(self.operator.direct(x)))
-    
-    ''' Same as in the parent class'''
-    def proximal(self,x, tau):
-        return self.function.proximal(x, tau)
-    
-    def proximal_conjugate(self,x, tau):
-        return self.function.proximal_conjugate(x, tau)
 
-    def convex_conjugate(self,x):
-        return self.function.convex_conjugate(x)    
-    
                     
-class FunctionComposition_new(Function):
+class TestFunction(unittest.TestCase):
+    def test_Function(self):
+    
     
-     def __init__(self, operator, *functions):
+        N = 3
+        ig = (N,N)
+        ag = ig       
+        op1 = Gradient(ig)
+        op2 = Identity(ig, ag)
+
+        # Form Composite Operator
+        operator = BlockOperator((2,1), op1, op2 ) 
         
-        self.functions = functions      
-        self.operator = operator
-        self.length = len(self.functions)
+        # Create functions
+        noisy_data = ImageData(np.random.randint(10, size=ag))
         
-#        if self.length==1:
-#            self.L = self.functions[0].alpha*(self.operator.norm()**2)    
-                        
-        # length == to operator.shape[0]#    
-        super(FunctionComposition_new, self).__init__()
-                    
-     def __call__(self, x):
-         
-        if isinstance(x, ImageData):
-            x =  CompositeDataContainer(x)         
+        d = ImageData(np.random.randint(10, size=ag))
         
-        t = 0                
-        for i in range(x.shape[0]):
-            t += self.functions[i](x.get_item(i))               
-        return t 
-    
-     def convex_conjugate(self, x):
-         
-        if isinstance(x, ImageData):
-            x =  CompositeDataContainer(x)         
+        f = mixed_L12Norm(alpha = 1).composition_with(op1)
+        g = L2NormSq(alpha=0.5, b=noisy_data)
         
-        t = 0                
-        for i in range(x.shape[0]):
-            t += self.functions[i].convex_conjugate(x.get_item(i))               
-        return t     
-                        
-     def proximal_conjugate(self, x, tau, out = None):
-     
-        if isinstance(x, ImageData):
-            x =  CompositeDataContainer(x)
-            
-        out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal_conjugate(x.get_item(i), tau)
+        # Compare call of f
+        a1 = ImageData(op1.direct(d).power(2).sum(axis=0)).sqrt().sum()
+        #print(a1, f(d))
+        self.assertEqual (a1, f(d))
         
-        if self.length==1:
-            return ImageData(*out)   
-        else:
-            return CompositeDataContainer(*out)   
-
-     def proximal(self, x, tau, out = None):
-     
-        if isinstance(x, ImageData):
-            x =  CompositeDataContainer(x)
-            
-        out = [None]*self.length
-        for i in range(self.length):
-            out[i] = self.functions[i].proximal(x.get_item(i), tau)
+        # Compare call of g
+        a2 = g.alpha*(d - noisy_data).power(2).sum()
+        #print(a2, g(d)) 
+        self.assertEqual(a2, g(d))
         
-        if self.length==1:
-            return ImageData(*out)   
-        else:
-            return CompositeDataContainer(*out)            
-    
-    
-if __name__ == '__main__':    
-    
-    N = 3
-    ig = (N,N)
-    ag = ig       
-    op1 = Gradient(ig)
-    op2 = Identity(ig, ag)
-
-    # Form Composite Operator
-    operator = CompositeOperator((2,1), op1, op2 ) 
-    
-    # Create functions
-    noisy_data = ImageData(np.random.randint(10, size=ag))
-    
-    d = ImageData(np.random.randint(10, size=ag))
-    
-    f = mixed_L12Norm(alpha = 1).composition_with(op1)
-    g = L2NormSq(alpha=0.5, b=noisy_data)
-    
-    # Compare call of f
-    a1 = ImageData(op1.direct(d).power(2).sum(axis=0)).sqrt().sum()
-    print(a1, f(d))
-    
-    # Compare call of g
-    a2 = g.alpha*(d - noisy_data).power(2).sum()
-    print(a2, g(d)) 
-    
-    # Compare convex conjugate of g
-    a3 = 0.5 * d.power(2).sum() + (d*noisy_data).sum()
-    print( a3, g.convex_conjugate(d))
+        # Compare convex conjugate of g
+        a3 = 0.5 * d.power(2).sum() + (d*noisy_data).sum()
+        self.assertEqual(a3, g.convex_conjugate(d))
+        #print( a3, g.convex_conjugate(d))