diff options
| -rw-r--r-- | Wrappers/Python/ccpi/optimisation/functions/KullbackLeibler.py | 126 |
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: + + + |