From 485ecfcd2fa1d88f27e0430c4685b4a463756d32 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Fri, 2 Dec 2011 16:34:22 +0000 Subject: [PATCH] Small updates --- pywafo/src/wafo/fig.py | 108 +++++++------ pywafo/src/wafo/kdetools.py | 296 ++++++++++++++++++++++++++++++++++-- pywafo/src/wafo/wafodata.py | 2 +- 3 files changed, 337 insertions(+), 69 deletions(-) diff --git a/pywafo/src/wafo/fig.py b/pywafo/src/wafo/fig.py index 77bcbb7..5dacb32 100644 --- a/pywafo/src/wafo/fig.py +++ b/pywafo/src/wafo/fig.py @@ -41,10 +41,10 @@ import win32con import msvcrt import numpy -__all__ = ['close','cycle','hide','keep','maximize','minimize','pile','restore','stack','tile'] +__all__ = ['close', 'cycle', 'hide', 'keep', 'maximize', 'minimize', 'pile', 'restore', 'stack', 'tile'] # Figure format strings to recognize in window title -_FIG_FORMATS = ('Figure','TVTK Scene','Chaco Plot Window: Figure') +_FIG_FORMATS = ('Figure', 'TVTK Scene', 'Chaco Plot Window: Figure') _SCREENSIZE = None def _getScreenSize(wnds): @@ -98,7 +98,7 @@ def _findTopWindows(wantedTitle=None): if wantedTitle == None: return topWindows else: - return [(hwnd,windowTxt) for hwnd, windowTxt in topWindows if windowTxt.startswith(wantedTitle)] + return [(hwnd, windowTxt) for hwnd, windowTxt in topWindows if windowTxt.startswith(wantedTitle)] def findallfigs(): ''' @@ -120,10 +120,10 @@ def findallfigs(): def _figparse(*args): figs = [] for arg in args: - if isinstance(arg,(list,tuple,set)): + if isinstance(arg, (list, tuple, set)): for val in arg: figs.append(int(val)) - elif isinstance(arg,int): + elif isinstance(arg, int): figs.append(arg) elif arg == 'all': figs = 'all' @@ -132,7 +132,7 @@ def _figparse(*args): raise TypeError('Only integers arguments accepted!') #raise TypeError('Unrecognized argument type (%s)!'%type(arg)) - if len(figs)==0 or figs =='all': + if len(figs) == 0 or figs == 'all': figs = findallfigs() return figs @@ -145,12 +145,12 @@ def _fig2wnd(figs): for fig in figs: for format_ in _FIG_FORMATS: winTitle = format_ + ' %d' % fig - hwnd = FindWindow(None,winTitle) + hwnd = FindWindow(None, winTitle) if not hwnd == 0: wnd_handles.append(hwnd) return wnd_handles -def _show_figure(figs,cmdshow): +def _show_figure(figs, cmdshow): ''' sets the specified figure's show state. @param figs: vector for figure numbers @@ -186,13 +186,13 @@ def _show_figure(figs,cmdshow): for fig in figs: for format_ in _FIG_FORMATS: winTitle = format_ + ' %d' % fig - hwnd = FindWindow(None,winTitle) + hwnd = FindWindow(None, winTitle) if not hwnd == 0: #ShowWindow(hwnd,cmdshow) BringWindowToTop(hwnd) - ShowWindow(hwnd,cmdshow) + ShowWindow(hwnd, cmdshow) -def _show_windows(wnds,cmdshow): +def _show_windows(wnds, cmdshow): ''' sets the specified window's show state. @param wnds: list of window handles numbers @@ -228,7 +228,7 @@ def _show_windows(wnds,cmdshow): if not hwnd == 0: #ShowWindow(hwnd,cmdshow) BringWindowToTop(hwnd) - ShowWindow(hwnd,cmdshow) + ShowWindow(hwnd, cmdshow) def keep(*figs): ''' Keeps figure windows of your choice and closes the rest. @@ -255,15 +255,15 @@ def keep(*figs): ''' figs2keep = [] for fig in figs: - if isinstance(fig,(list,tuple,set)): + if isinstance(fig, (list, tuple, set)): for val in fig: figs2keep.append(int(val)) - elif isinstance(fig,int): + elif isinstance(fig, int): figs2keep.append(fig) else: raise TypeError('Only integers arguments accepted!') - if len(figs2keep)>0: + if len(figs2keep) > 0: allfigs = set(findallfigs()) # Remove figure handles in the "keep" list @@ -476,16 +476,16 @@ def pile(*figs): figlist = _figparse(*figs) wnds = _fig2wnd(figlist) numfigs = len(wnds) - if numfigs>0: + if numfigs > 0: pos = _getScreenSize(wnds) - pos[3] = int(pos[3]/2) - pos[2] = int(pos[2]/2.5) - pos[1] = int(pos[3]/2) - pos[0] = int(pos[2]/2) + pos[3] = int(pos[3] / 2) + pos[2] = int(pos[2] / 2.5) + pos[1] = int(pos[3] / 2) + pos[0] = int(pos[2] / 2) BringWindowToTop = win32gui.BringWindowToTop MoveWindow = win32gui.MoveWindow for wnd in wnds: - MoveWindow(wnd, pos[0], pos[1], pos[2], pos[3],1) + MoveWindow(wnd, pos[0], pos[1], pos[2], pos[3], 1) BringWindowToTop(wnd) @@ -521,17 +521,17 @@ def stack(*figs): figlist = _figparse(*figs) wnds = _fig2wnd(figlist) numfigs = len(wnds) - if numfigs>0: + if numfigs > 0: screenpos = _getScreenSize(wnds) - maxfigs = numpy.fix(screenpos[3]/20) + maxfigs = numpy.fix(screenpos[3] / 20) - if (numfigs>maxfigs): # figure limit check - print(' More than %d requested '% maxfigs) + if (numfigs > maxfigs): # figure limit check + print(' More than %d requested ' % maxfigs) return BringWindowToTop = win32gui.BringWindowToTop MoveWindow = win32gui.MoveWindow - GetWindowRect = win32gui.GetWindowRect + GetWindowRect = win32gui.GetWindowRect # # #% tile figures by postiion @@ -543,13 +543,13 @@ def stack(*figs): pos[3] -= pos[1] pos[2] -= pos[0] #print('[x, y, w, h] = ', pos) - ypos = screenpos[1]+iy*20 #int(screenpos[3] - iy*20 -pos[3] -70) # figure location (row) - xpos = int(iy*5 + 15+screenpos[0]) # figure location (column) - MoveWindow(wnds[iy],xpos,ypos, pos[2],pos[3],1) + ypos = screenpos[1] + iy * 20 #int(screenpos[3] - iy*20 -pos[3] -70) # figure location (row) + xpos = int(iy * 5 + 15 + screenpos[0]) # figure location (column) + MoveWindow(wnds[iy], xpos, ypos, pos[2], pos[3], 1) BringWindowToTop(wnds[iy]) -def tile(*figs,**kwds): +def tile(*figs, **kwds): ''' Tile figure windows. Parameters @@ -587,17 +587,17 @@ def tile(*figs,**kwds): nfigs = len(wnds); # Number of windows. - if nfigs>0: - nfigspertile = kwds.get('pairs',nfigs) + if nfigs > 0: + nfigspertile = kwds.get('pairs', nfigs) ceil = numpy.ceil sqrt = numpy.sqrt maximum = numpy.maximum - nlayers = int(ceil(nfigs/nfigspertile)) + nlayers = int(ceil(nfigs / nfigspertile)) nh = int(ceil(sqrt(nfigspertile))) # Number of figures horisontally. - nv = int(ceil(nfigspertile/nh)); # Number of figures vertically. + nv = int(ceil(nfigspertile / nh)); # Number of figures vertically. nh = maximum(nh, 2); @@ -622,15 +622,15 @@ def tile(*figs,**kwds): #% 4 - Window vertical size #%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - hspc = 10 # Horisontal space. + hspc = 10 # Horisontal space. topspc = 20; # Space above top figure. medspc = 10; # Space between figures. botspc = 20; # Space below bottom figure. #print('scrwid = %d' % scrwid) - figwid = ( scrwid - (nh+1)*hspc )/nh + figwid = (scrwid - (nh + 1) * hspc) / nh #print('figwid = %d' % figwid) - fighgt = ( scrhgt - (topspc+botspc) - (nv-1)*medspc )/nv; + fighgt = (scrhgt - (topspc + botspc) - (nv - 1) * medspc) / nv; figwid = int(numpy.round(figwid)) fighgt = int(numpy.round(fighgt)) @@ -642,19 +642,19 @@ def tile(*figs,**kwds): for unused_ix in xrange(nlayers): for row in xrange(nv): for col in xrange(nh): - if (row)*nh + col < nfigspertile : - if idx0: - maximize = kwds.get('maximize',False) - pairs = kwds.get('pairs',1) + if numfigs > 0: + maximize = kwds.get('maximize', False) + pairs = kwds.get('pairs', 1) if maximize or pairs == None: nfigspercycle = 1 @@ -713,18 +713,18 @@ def cycle(*figs,**kwds): i = 0; escape_key = chr(27) backspace_key = chr(8) - while 0<=i and i 1 else In + + kw = np.zeros((inc,)*d) + for i in range(d): + kw += exp(-Inc[i]) + y = kwds.get('y', 1.0) + d, n = self.dataset.shape + # Find the binned kernel weights, c. + c = gridcount(self.dataset, X, y=y)/n + # Perform the convolution. + at = dctn(c) * kw + z = idctn(at)*at.size/np.prod(R) + return z*(z>0.0) + #ix = (slice(0, inc),)*d + #return z[ix] * (z[ix] > 0.0) + + def _eval_grid_fun(self, eval_grd, *args, **kwds): + output = kwds.pop('output', 'value') + f = eval_grd(*args, **kwds) + if output == 'value': + return f + else: + titlestr = 'Kernel density estimate (%s)' % self.kernel.name + kwds2 = dict(title=titlestr) + kwds2['plot_kwds'] = dict(plotflag=1) + kwds2.update(**kwds) + args = self.args + if self.d == 1: + args = args[0] + wdata = WafoData(f, args, **kwds2) + if self.d > 1: + PL = np.r_[10:90:20, 95, 99, 99.9] + try: + ql = qlevels(f, p=PL) + wdata.clevels = ql + wdata.plevels = PL + except: + pass + return wdata + def _check_shape(self, points): + points = atleast_2d(points) + d, m = points.shape + if d != self.d: + if d == 1 and m == self.d: + # points was passed in as a row vector + points = np.reshape(points, (self.d, 1)) + else: + msg = "points have dimension %s, dataset has dimension %s" % (d, + self.d) + raise ValueError(msg) + return points + def eval_points(self, points, **kwds): + """Evaluate the estimated pdf on a set of points. + + Parameters + ---------- + points : (# of dimensions, # of points)-array + Alternatively, a (# of dimensions,) vector can be passed in and + treated as a single point. + + Returns + ------- + values : (# of points,)-array + The values at each point. + + Raises + ------ + ValueError if the dimensionality of the input points is different than + the dimensionality of the KDE. + """ + + points = self._check_shape(points) + return self._eval_points(points, **kwds) + + def _eval_points(self, points, **kwds): + pass + + __call__ = eval_grid_fast class _KDE(object): """ Kernel-Density Estimator base class. @@ -341,7 +554,7 @@ class TKDE(_KDE): def __init__(self, data, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128, L2=None): self.L2 = L2 - _KDE.__init__(self, data, hs, kernel, alpha, xmin, xmax, inc) + super(TKDE, self).__init__(data, hs, kernel, alpha, xmin, xmax, inc) def _initialize(self): self._check_xmin() @@ -574,7 +787,7 @@ class KDE(_KDE): """ def __init__(self, data, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128): - _KDE.__init__(self, data, hs, kernel, alpha, xmin, xmax, inc) + super(KDE, self).__init__(data, hs, kernel, alpha, xmin, xmax, inc) def _initialize(self): self._compute_smoothing() @@ -638,11 +851,14 @@ class KDE(_KDE): Xn = np.dot(self.inv_hs, np.vstack(Xnc)) # Obtain the kernel weights. + kw = self.kernel(Xn) + norm_fact = (kw.sum()*dx.prod()*self.n) + norm_fact2 = (self._norm_factor * self.kernel.norm_factor(d, self.n)) + + kw = kw/norm_fact r = kwds.get('r', 0) - if r==0: - kw = self.kernel(Xn) / (self._norm_factor * self.kernel.norm_factor(d, self.n)) - else: - kw = np.vstack(Xnc) ** r * self.kernel(Xn) / (self._norm_factor * self.kernel.norm_factor(d, self.n)) + if r!=0: + kw *= np.vstack(Xnc) ** r kw.shape = shape0 kw = np.fft.ifftshift(kw) fftn = np.fft.fftn @@ -1434,7 +1650,8 @@ class Kernel(object): c = gridcount(A[dim], xa) N = len(set(A[dim])) - a = dct(c/c.sum(), norm=None) + #a = dct(c/c.sum(), norm=None) + a = dct(c/len(A[dim]), norm=None) #% now compute the optimal bandwidth^2 using the referenced method I = np.asfarray(np.arange(1, inc))**2 @@ -2443,7 +2660,7 @@ def gridcount(data, X, y=1): >>> import numpy as np >>> import wafo.kdetools as wk >>> import pylab as plb - >>> N = 20; + >>> N = 20; >>> data = np.random.rayleigh(1,N) >>> x = np.linspace(0,max(data)+1,50) >>> dx = x[1]-x[0] @@ -2551,8 +2768,8 @@ def kde_demo1(): import scipy.stats as st x = np.linspace(-4, 4, 101) x0 = x / 2.0 - data = np.random.normal(loc=0, scale=1.0, size=7) #rndnorm(0,1,7,1); - kernel = Kernel('gaus') + data = np.random.normal(loc=0, scale=1.0, size=7) + kernel = Kernel('gauss') hs = kernel.hns(data) hVec = [hs / 2, hs, 2 * hs] @@ -2742,7 +2959,59 @@ def kreg_demo1(hs=None, fast=False, fun='hisj'): print(kreg.tkde.tkde.inv_hs) print(kreg.tkde.tkde.hs) - + +def kde_gauss_demo(n=50000): + ''' + KDEDEMO Demonstrate the KDEgauss + + KDEDEMO1 shows the true density (dotted) compared to KDE based on 7 + observations (solid) and their individual kernels (dashed) for 3 + different values of the smoothing parameter, hs. + ''' + + import scipy.stats as st + #x = np.linspace(-4, 4, 101) + #data = np.random.normal(loc=0, scale=1.0, size=n) + #data = np.random.exponential(scale=1.0, size=n) +# n1 = 128 +# I = (np.arange(n1)*pi)**2 *0.01*0.5 +# kw = exp(-I) +# pylab.plot(idctn(kw)) +# return + dist = st.norm + dist = st.expon + data = dist.rvs(loc=0, scale=1.0, size=n) + d, N = np.atleast_2d(data).shape + + if d==1: + plot_options = [dict(color='red'), dict(color='green'), dict(color='black')] + else: + plot_options = [dict(colors='red'), dict(colors='green'), dict(colors='black')] + + pylab.figure(1) + kde0 = KDE(data, kernel=Kernel('gauss', 'hisj')) + f0 = kde0.eval_grid_fast(output='plot', ylab='Density') + f0.plot(**plot_options[0]) + + kde1 = TKDE(data, kernel=Kernel('gauss', 'hisj'), L2=0) + f1 = kde1.eval_grid_fast(output='plot', ylab='Density') + f1.plot(**plot_options[1]) + + kde2 = KDEgauss(data) + f2 = kde2(output='plot', ylab='Density') + x = f2.args + f2.plot(**plot_options[2]) + + fmax = dist.pdf(x, 0, 1).max() + if d==1: + pylab.plot(x, dist.pdf(x, 0, 1), 'k:') + pylab.axis([x.min(), x.max(), 0, fmax]) + pylab.show() + print(fmax/f2.data.max()) + format_ = ''.join(('%g, ')*d) + format_ = 'hs0=%s hs2=%s' % (format_,format_) + print(format_ % tuple(kde0.hs.tolist()+kde2.hs.tolist())) + def test_docstrings(): import doctest doctest.testmod() @@ -2750,4 +3019,5 @@ def test_docstrings(): if __name__ == '__main__': #test_docstrings() #kde_demo2() - kreg_demo1() \ No newline at end of file + #kreg_demo1(fast=True) + kde_gauss_demo() \ No newline at end of file diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py index 062cec1..2440903 100644 --- a/pywafo/src/wafo/wafodata.py +++ b/pywafo/src/wafo/wafodata.py @@ -190,7 +190,7 @@ class Plotter_1d(object): def plot(self, wdata, *args, **kwds): plotflag = kwds.pop('plotflag', False) if plotflag: - h1 = self._plot(plotflag, wdata, **kwds) + h1 = self._plot(plotflag, wdata, *args, **kwds) else: if isinstance(wdata.args, (list, tuple)): args1 = tuple((wdata.args)) + (wdata.data,) + args