fix call conv_conj

1 files changed, 105 insertions, 21 deletions

diff --git a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
index 4dd77cf..0d3c8f5 100644
--- a/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
+++ b/Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py
@@ -20,33 +20,42 @@
 import numpy
 from ccpi.optimisation.functions import Function
 from ccpi.optimisation.functions.ScaledFunction import ScaledFunction 
-from ccpi.framework import ImageData, ImageGeometry
 import functools
+import scipy.special
 
 class KullbackLeibler(Function):
     
-    ''' Assume that data >= 0
+    ''' 
+    
+    KL_div(x, y + back) = int x * log(x/(y+back)) - x + (y+back)
+    
+    Assumption: y>=0
+                back>=0
                 
     '''
     
-    def __init__(self,data, **kwargs):
+    def __init__(self, data, **kwargs):
         
         super(KullbackLeibler, self).__init__()
         
-        self.b = data        
-        self.bnoise = kwargs.get('bnoise', 0)
-
-                                                
+        self.b = data    
+        self.bnoise = 0
+        
+                                                    
     def __call__(self, x):
-                
-        res_tmp = numpy.zeros(x.shape)
         
-        tmp = x + self.bnoise 
-        ind = x.as_array()>0
         
-        res_tmp[ind] = x.as_array()[ind] - self.b.as_array()[ind] * numpy.log(tmp.as_array()[ind])
+        '''
         
-        return res_tmp.sum()
+            x - y * log( x + bnoise) + y * log(y) - y + bnoise
+        
+        
+        '''
+        
+        # TODO avoid scipy import ????
+        tmp = scipy.special.kl_div(self.b.as_array(), x.as_array())                
+        return numpy.sum(tmp) 
+          
 
     def log(self, datacontainer):
         '''calculates the in-place log of the datacontainer'''
@@ -58,7 +67,6 @@ class KullbackLeibler(Function):
         
     def gradient(self, x, out=None):
         
-        #TODO Division check
         if out is None:
             return 1 - self.b/(x + self.bnoise)
         else:
@@ -70,12 +78,10 @@ class KullbackLeibler(Function):
             
     def convex_conjugate(self, x):
         
-        tmp = self.b/(1-x)
-        ind = tmp.as_array()>0
-        
-        return (self.b.as_array()[ind] * (numpy.log(tmp.as_array()[ind])-1)).sum()
-
-    
+        # TODO avoid scipy import ????
+        xlogy = scipy.special.xlogy(self.b.as_array(), 1 - x.as_array())
+        return numpy.sum(-xlogy)
+            
     def proximal(self, x, tau, out=None):
         
         if out is None:        
@@ -140,8 +146,86 @@ if __name__ == '__main__':
     M, N =  2,3
     ig = ImageGeometry(voxel_num_x=M, voxel_num_y = N)
     u = ig.allocate('random_int')
-    b = np.random.normal(0, 0.1, size=ig.shape)  
+    b = ig.allocate('random_int')
+    u.as_array()[1,1]=0
+    u.as_array()[2,0]=0
+    b.as_array()[1,1]=0
+    b.as_array()[2,0]=0    
+    
+    f = KullbackLeibler(b)
+    
+    
+#    longest = reduce(lambda x, y: len(x) if len(x) > len(y) else len(y), strings)
+
+
+#    tmp = functools.reduce(lambda x, y: \
+#                           0 if x==0 and not numpy.isnan(y) else x * numpy.log(y), \
+#                           zip(b.as_array().ravel(), u.as_array().ravel()),0)
+    
+    
+#    np.multiply.reduce(X, 0)
+    
+    
+#                sf = reduce(lambda x,y: x + y[0]*y[1],
+#                            zip(self.as_array().ravel(),
+#                                other.as_array().ravel()),
+#                            0)        
+#cdef inline number_t xlogy(number_t x, number_t y) nogil:
+#    if x == 0 and not zisnan(y):
+#        return 0
+#    else:
+#        return x * zlog(y)        
+    
+#    if npy_isnan(x):
+#        return x
+#    elif x > 0:
+#        return -x * log(x)
+#    elif x == 0:
+#        return 0
+#    else:
+#        return -inf    
+    
+#        cdef inline double kl_div(double x, double y) nogil:
+#    if npy_isnan(x) or npy_isnan(y):
+#        return nan
+#    elif x > 0 and y > 0:
+#        return x * log(x / y) - x + y
+#    elif x == 0 and y >= 0:
+#        return y
+#    else:
+#        return inf    
+
+    
+    
+    
+#    def xlogy(self, dc1, dc2):
+        
+#        return numpy.sum(numpy.where(dc1.as_array() != 0, dc2.as_array() * numpy.log(dc2.as_array() / dc1.as_array()), 0))
         
+           
+    
+#    f.xlog(u, b)
+    
+            
+
+    
+#    tmp1 = b.as_array()
+#    tmp2 = u.as_array()
+#    
+#    zz = scipy.special.xlogy(tmp1, tmp2)
+#
+#    print(np.sum(zz))
+    
+    
+#    ww = f.xlogy(b, u)
+    
+#    print(ww)
+    
+    
+#cdef inline double kl_div(double x, double y) nogil:
+  
+    
+