added to functions

4 files changed, 301 insertions, 1 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
new file mode 100755
index 0000000..df8dc89
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/functions/IndicatorBox.py
@@ -0,0 +1,65 @@
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018-2019 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+from ccpi.optimisation.functions import Function
+import numpy
+
+class IndicatorBox(Function):
+    '''Box constraints indicator function. 
+    
+    Calling returns 0 if argument is within the box. The prox operator is projection onto the box. 
+    Only implements one scalar lower and one upper as constraint on all elements. Should generalise
+    to vectors to allow different constraints one elements.
+'''
+    
+    def __init__(self,lower=-numpy.inf,upper=numpy.inf):
+        # Do nothing
+        super(IndicatorBox, self).__init__()
+        self.lower = lower
+        self.upper = upper
+        
+    
+    def __call__(self,x):
+        
+        if (numpy.all(x.array>=self.lower) and 
+            numpy.all(x.array <= self.upper) ):
+            val = 0
+        else:
+            val = numpy.inf
+        return val
+    
+    def prox(self,x,tau=None):
+        return  (x.maximum(self.lower)).minimum(self.upper)
+    
+    def proximal(self, x, tau, out=None):
+        if out is None:
+            return self.prox(x, tau)
+        else:
+            if not x.shape == out.shape:
+                raise ValueError('Norm1 Incompatible size:',
+                                 x.shape, out.shape)
+            #(x.abs() - tau*self.gamma).maximum(0) * x.sign()
+            x.abs(out = out)
+            out.__isub__(tau*self.gamma)
+            out.maximum(0, out=out)
+            if self.sign_x is None or not x.shape == self.sign_x.shape:
+                self.sign_x = x.sign()
+            else:
+                x.sign(out=self.sign_x)
+                
+            out.__imul__( self.sign_x )
diff --git a/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
new file mode 100755
index 0000000..1c51236
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/functions/MixedL21Norm.py
@@ -0,0 +1,136 @@
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018-2019 Evangelos Papoutsellis and Edoardo Pasca
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import numpy as np
+from ccpi.optimisation.functions import Function, ScaledFunction
+from ccpi.framework import DataContainer, ImageData, \
+                           ImageGeometry, BlockDataContainer 
+
+############################   mixed_L1,2NORM FUNCTIONS   #####################
+class MixedL21Norm(Function):
+    
+    def __init__(self, **kwargs):
+
+        super(MixedL21Norm, self).__init__()                      
+        self.SymTensor = kwargs.get('SymTensor',False)
+        
+    def __call__(self, x, out=None):
+        
+        ''' Evaluates L1,2Norm at point x
+            
+            :param: x is a BlockDataContainer
+                                
+        '''                        
+        if self.SymTensor:
+            
+            param = [1]*x.shape[0]
+            param[-1] = 2
+            tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])]
+            res = sum(tmp).sqrt().sum()           
+        else:
+            
+#            tmp = [ x[i]**2 for i in range(x.shape[0])]
+            tmp = [ el**2 for el in x.containers ]
+            
+#            print(x.containers)
+#            print(tmp)
+#            print(type(sum(tmp)))
+#            print(type(tmp))
+            res = sum(tmp).sqrt().sum()
+#            print(res)
+        return res           
+                            
+    def gradient(self, x, out=None):
+        return ValueError('Not Differentiable')
+                            
+    def convex_conjugate(self,x):
+        
+        ''' This is the Indicator function of ||\cdot||_{2, \infty}
+            which is either 0 if ||x||_{2, \infty} or \infty        
+        '''
+        return 0.0
+    
+    def proximal(self, x, tau, out=None):
+        
+        '''
+            For this we need to define a MixedL2,2 norm acting on BDC,
+            different form L2NormSquared which acts on DC
+        
+        '''
+        
+        pass
+    
+    def proximal_conjugate(self, x, tau, out=None): 
+        
+        if self.SymTensor:
+            
+            param = [1]*x.shape[0]
+            param[-1] = 2
+            tmp = [param[i]*(x[i] ** 2) for i in range(x.shape[0])]
+            frac = [x[i]/(sum(tmp).sqrt()).maximum(1.0) for i in range(x.shape[0])]
+            res = BlockDataContainer(*frac) 
+            
+            return res
+            
+#                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:
+            
+                tmp = [ el*el for el in x]
+                res = (sum(tmp).sqrt()).maximum(1.0) 
+                frac = [x[i]/res for i in range(x.shape[0])]
+                res = BlockDataContainer(*frac)
+                                                   
+        return res 
+    
+    def __rmul__(self, scalar):
+        return ScaledFunction(self, scalar) 
+
+#class MixedL21Norm_tensor(Function):
+#    
+#    def __init__(self):
+#        print("feerf")
+#    
+#
+if __name__ == '__main__':
+    
+    M, N, K = 2,3,5
+    ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
+    u1 = ig.allocate('random_int')
+    u2 = ig.allocate('random_int')
+    
+    U = BlockDataContainer(u1, u2, shape=(2,1))
+    
+    # Define no scale and scaled
+    f_no_scaled = MixedL21Norm() 
+    f_scaled = 0.5 * MixedL21Norm()  
+    
+    # call
+    
+    a1 = f_no_scaled(U)
+    a2 = f_scaled(U)    
+    
+    z = f_no_scaled.proximal_conjugate(U, 1)
+    
+    f_no_scaled = MixedL21Norm()
+    
+    tmp = [el*el for el in U]
+    
+
+    
diff --git a/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
new file mode 100755
index 0000000..b553d7c
--- /dev/null
+++ b/Wrappers/Python/ccpi/optimisation/functions/Norm2Sq.py
@@ -0,0 +1,98 @@
+# -*- coding: utf-8 -*-
+#   This work is part of the Core Imaging Library developed by
+#   Visual Analytics and Imaging System Group of the Science Technology
+#   Facilities Council, STFC
+
+#   Copyright 2018-2019 Jakob Jorgensen, Daniil Kazantsev and Edoardo Pasca
+
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+from ccpi.optimisation.functions import Function
+import numpy
+import warnings
+
+# Define a class for squared 2-norm
+class Norm2sq(Function):
+    '''
+    f(x) = c*||A*x-b||_2^2
+    
+    which has 
+    
+    grad[f](x) = 2*c*A^T*(A*x-b)
+    
+    and Lipschitz constant
+    
+    L = 2*c*||A||_2^2 = 2*s1(A)^2
+    
+    where s1(A) is the largest singular value of A.
+    
+    '''
+    
+    def __init__(self,A,b,c=1.0,memopt=False):
+        super(Norm2sq, self).__init__()
+    
+        self.A = A  # Should be an operator, default identity
+        self.b = b  # Default zero DataSet?
+        self.c = c  # Default 1.
+        if memopt:
+            try:
+                self.range_tmp = A.range_geometry().allocate()
+                self.domain_tmp = A.domain_geometry().allocate()
+                self.memopt = True
+            except NameError as ne:
+                warnings.warn(str(ne))
+                self.memopt = False
+            except NotImplementedError as nie:
+                print (nie)
+                warnings.warn(str(nie))
+                self.memopt = False
+        else:
+            self.memopt = False
+        
+        # Compute the Lipschitz parameter from the operator if possible
+        # Leave it initialised to None otherwise
+        try:
+            self.L = 2.0*self.c*(self.A.norm()**2)
+        except AttributeError as ae:
+            pass
+        except NotImplementedError as noe:
+            pass
+        
+    #def grad(self,x):
+    #    return self.gradient(x, out=None)
+
+    def __call__(self,x):
+        #return self.c* np.sum(np.square((self.A.direct(x) - self.b).ravel()))
+        #if out is None:
+        #    return self.c*( ( (self.A.direct(x)-self.b)**2).sum() )
+        #else:
+        y = self.A.direct(x)
+        y.__isub__(self.b)
+        #y.__imul__(y)
+        #return y.sum() * self.c
+        try:
+            return y.squared_norm() * self.c
+        except AttributeError as ae:
+            # added for compatibility with SIRF 
+            return (y.norm()**2) * self.c
+    
+    def gradient(self, x, out = None):
+        if self.memopt:
+            #return 2.0*self.c*self.A.adjoint( self.A.direct(x) - self.b )
+            
+            self.A.direct(x, out=self.range_tmp)
+            self.range_tmp -= self.b 
+            self.A.adjoint(self.range_tmp, out=out)
+            #self.direct_placehold.multiply(2.0*self.c, out=out)
+            out *= (self.c * 2.0)
+        else:
+            return (2.0*self.c)*self.A.adjoint( self.A.direct(x) - self.b )
diff --git a/Wrappers/Python/ccpi/optimisation/functions/__init__.py b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
index d6edd03..2ed36f5 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/__init__.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/__init__.py
@@ -9,4 +9,5 @@ from .BlockFunction import BlockFunction
 from .ScaledFunction import ScaledFunction
 from .FunctionOperatorComposition import FunctionOperatorComposition
 from .MixedL21Norm import MixedL21Norm
-
+from .IndicatorBox import IndicatorBox
+from .Norm2Sq import Norm2sq