Small updates

pywafo/src/wafo/kdetools.py

@@ -10,10 +10,12 @@
 #-------------------------------------------------------------------------------
 #!/usr/bin/env python
 from __future__ import division
+import copy
 import warnings
 import numpy as np
 from numpy import pi, sqrt, atleast_2d, exp, newaxis, array #@UnresolvedImport
 import scipy
+from scipy import interpolate
 from scipy import linalg
 from scipy.special import gamma
 from misc import tranproc, trangood
@@ -112,6 +114,14 @@ class TKDE(object):
     array([ 1.03982714,  0.45839018,  0.39514782,  0.32860602,  0.26433318,
             0.20717946,  0.15907684,  0.1201074 ,  0.08941027,  0.06574882])
             
+    >>> kde.eval_grid(x)
+    array([ 1.03982714,  0.45839018,  0.39514782,  0.32860602,  0.26433318,
+            0.20717946,  0.15907684,  0.1201074 ,  0.08941027,  0.06574882])
+            
+    >>> kde.eval_grid_fast(x)
+    array([ 0.        ,  1.16200356,  0.99256178,  0.81930973,  0.65479862,
+            0.51021576,  0.3896221 ,  0.29266142,  0.        ,  0.        ])
+            
     import pylab as plb          
     h1 = plb.plot(x, f) #  1D probability density plot
     t = np.trapz(f, x)   
@@ -125,11 +135,13 @@ class TKDE(object):
         self.alpha = alpha
         self.xmin = xmin
         self.xmax = xmax
+        self.inc = inc
         self.L2 = L2
         self.d, self.n = self.dataset.shape
         self.initialize()
     
     def initialize(self):
+        self._set_xlimits()
         tdataset = self._dat2gaus(self.dataset)
         xmin = self.xmin
         if xmin is not None:
@@ -137,8 +149,26 @@ class TKDE(object):
         xmax = self.xmax
         if xmax is not None:
             xmax = self._dat2gaus(xmax)
-        self.kde = KDE(tdataset, self.hs, self.kernel, self.alpha, xmin, xmax,
+        self.tkde = KDE(tdataset, self.hs, self.kernel, self.alpha, xmin, xmax,
                        self.inc)
+    def _set_xlimits(self):
+        amin = self.dataset.min(axis=-1)
+        amax = self.dataset.max(axis=-1)
+        xyzrange = amax-amin
+        offset = xyzrange/4.0
+        if self.xmin is None:
+            self.xmin = amin - offset
+        else:
+            self.xmin = self.xmin * np.ones(self.d)
+        if self.xmax is None:
+            self.xmax = amax + offset
+        else:
+            self.xmax = self.xmax * np.ones(self.d)
+        
+        if self.L2 is not None:
+            L2 = np.atleast_1d(self.L2) * np.ones(self.d) # default no transformation
+            self.xmin = np.where(L2!=1, np.maximum(self.xmin, amin/2.0), self.xmin)
+            
         
     def _check_shape(self, points):
         points = atleast_2d(points)
@@ -160,11 +190,22 @@ class TKDE(object):
         
         L2 = np.atleast_1d(self.L2) * np.ones(self.d) # default no transformation
         
-        tpoints = points.copy()
+        tpoints = copy.copy(points)
         for i, v2 in enumerate(L2.tolist()):
-            tpoints[i] = np.where(v2 == 0, np.log(points[i]), points[i] ** v2)
+            tpoints[i] = np.log(points[i]) if v2 == 0 else points[i] ** v2
         return tpoints
     
+    def _gaus2dat(self, tpoints):
+        if self.L2 is None:
+            return tpoints # default no transformation
+        
+        L2 = np.atleast_1d(self.L2) * np.ones(self.d) # default no transformation
+        
+        points = copy.copy(tpoints)
+        for i, v2 in enumerate(L2.tolist()):
+            points[i] = np.exp(tpoints[i]) if v2 == 0 else tpoints[i] ** (1.0/v2)
+        return points
+    
     def _scale_pdf(self, pdf, points):
         if self.L2 is None:
             return pdf
@@ -177,7 +218,70 @@ class TKDE(object):
                     transformation. Check the KDE for spurious spikes'''
             warnings.warn(msg)
         return pdf
+    def eval_grid_fast(self, *args):
+        """Evaluate the estimated pdf on a grid.
 
+        Parameters
+        ----------
+        arg_0,arg_1,... arg_d-1 : vectors
+            Alternatively, if no vectors is passed in then
+             arg_i = linspace(self.xmin[i], self.xmax[i], self.inc)
+
+        Returns
+        -------
+        values : array-like
+            The values evaluated at meshgrid(*args).
+
+        """
+        
+        if self.L2 is None:
+            f = self.tkde.eval_grid_fast(*args)
+            self.args = self.tkde.args
+            return f
+        #targs = self._dat2gaus(list(args)) if len(args) else args
+        tf = self.tkde.eval_grid_fast()
+        self.args = self._gaus2dat(list(self.tkde.args))
+        points = meshgrid(*self.args) if self.d>1 else self.args
+        f = self._scale_pdf(tf, points)
+        if len(args):
+            if self.d==1:
+                pdf = interpolate.interp1d(points[0],f, bounds_error=False, fill_value=0.0)
+            elif self.d==2:
+                pdf = interpolate.interp2d(points[0],points[1], f, bounds_error=False, fill_value=0.0)
+            ipoints = meshgrid(*args) if self.d>1 else args
+            fi = pdf(*ipoints)
+            fi.shape = ipoints[0].shape
+            return fi
+        return f
+    def eval_grid(self, *args):
+        """Evaluate the estimated pdf on a grid.
+
+        Parameters
+        ----------
+        arg_0,arg_1,... arg_d-1 : vectors
+            Alternatively, if no vectors is passed in then
+             arg_i = linspace(self.xmin[i], self.xmax[i], self.inc)
+
+        Returns
+        -------
+        values : array-like
+            The values evaluated at meshgrid(*args).
+
+        """
+        if len(args)==0:
+            args = []
+            for i in range(self.d):
+                args.append(np.linspace(self.xmin[i], self.xmax[i], self.inc))
+        self.args = args
+        if self.L2 is None:
+            return self.tkde.eval_grid(*args)
+        targs = self._dat2gaus(list(args))
+        tf = self.tkde.eval_grid(*targs)
+        points = meshgrid(*args) if self.d>1 else self.args
+        f = self._scale_pdf(tf, points)
+        return f
+        
+        return self.tkde.eval_grid(*args)
     def evaluate(self, points):
         """Evaluate the estimated pdf on a set of points.
 
@@ -198,10 +302,10 @@ class TKDE(object):
         the dimensionality of the KDE.
         """
         if self.L2 is None:
-            return self.kde(points)
+            return self.tkde(points)
         points = self._check_shape(points)
         tpoints = self._dat2gaus(points)
-        tf = self.kde(tpoints)
+        tf = self.tkde(tpoints)
         f = self._scale_pdf(tf, points)
         return f
     
@@ -258,6 +362,24 @@ class KDE(object):
     array([ 0.17252055,  0.41014271,  0.61349072,  0.57023834,  0.37198073,
             0.21409279,  0.12738463,  0.07460326,  0.03956191,  0.01887164])
     
+    >>> kde.eval_grid(x)
+    array([ 0.17252055,  0.41014271,  0.61349072,  0.57023834,  0.37198073,
+            0.21409279,  0.12738463,  0.07460326,  0.03956191,  0.01887164])
+    
+    >>> kde0 = wk.KDE(data, hs=0.5, alpha=0.0)
+    >>> kde0.evaluate(x)
+    array([ 0.2039735 ,  0.40252503,  0.54595078,  0.52219649,  0.3906213 ,
+            0.26381501,  0.16407362,  0.08270612,  0.02991145,  0.00720821])
+    
+    >>> kde0.eval_grid(x)
+    array([ 0.2039735 ,  0.40252503,  0.54595078,  0.52219649,  0.3906213 ,
+            0.26381501,  0.16407362,  0.08270612,  0.02991145,  0.00720821])       
+    
+    >>> f = kde0.eval_grid_fast()
+    >>> np.interp(x, kde0.args[0], f)
+    array([ 0.54344   ,  1.04793706,  1.38047458,  1.28324876,  0.943131  ,
+            0.63570524,  0.39404219,  0.19450807,  0.08505344,  0.08505344])
+            
     import pylab as plb          
     h1 = plb.plot(x, f) #  1D probability density plot
     t = np.trapz(f, x)   
@@ -345,8 +467,10 @@ class KDE(object):
 
         """
         if len(args)==0:
+            args = []
             for i in range(self.d):
                 args.append(np.linspace(self.xmin[i], self.xmax[i], self.inc))
+        self.args = args
         return self._eval_grid_fast(*args)
     def _eval_grid_fast(self, *args):
         X = np.vstack(args)
@@ -360,16 +484,16 @@ class KDE(object):
         for i in range(d):
             Xn.append(x0*dx[i])
         
-        Xnc = meshgrid(*Xn)
+        Xnc = meshgrid(*Xn) if d>1 else Xn
         
         shape0 = Xnc[0].shape
         for i in range(d):
             Xnc[i].shape = (-1,)
  
-        Xn = np.dot(np.vstack(Xnc), self.inv_hs)
+        Xn = np.dot(self.inv_hs, np.vstack(Xnc))
         
         # Obtain the kernel weights.        
-        kw = self.kernel(Xn)
+        kw = self.kernel(Xn)/self._norm_factor
         kw.shape = shape0
         kw = np.fft.ifftshift(kw)
         fftn = np.fft.fftn
@@ -401,13 +525,15 @@ class KDE(object):
         """
 
         if len(args)==0:
+            args = []
             for i in range(self.d):
                 args.append(np.linspace(self.xmin[i], self.xmax[i], self.inc))
+        self.args = args
         return self._eval_grid(*args)
     
     def _eval_grid(self, *args):
         
-        grd = meshgrid(*args)
+        grd = meshgrid(*args) if len(args)>1 else list(args)
         shape0 = grd[0].shape
         d = len(grd)
         for i in range(d):

pywafo/src/wafo/test/test_kdetools.py

@@ -166,7 +166,6 @@ def test_gridcount_2D():
            [ 0.02063536,  0.31054405,  0.71865964,  0.13486633,  0.        ],
            [ 0.        ,  0.        ,  0.        ,  0.        ,  0.        ]])
     
-    
     h = plb.plot(x, c, '.')   # 1D histogram
    
     h1 = plb.plot(x, c / dx / N) #  1D probability density plot