From ef856af7ae10c0f5fa0528ba916dc7c7b87e6da9 Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Thu, 22 Dec 2011 02:16:49 +0000 Subject: [PATCH] Completed smoothn. Tested on 3 examples and works OK. --- pywafo/src/wafo/kdetools.py | 443 ++++++++++++++++++---------------- pywafo/src/wafo/polynomial.py | 2 +- 2 files changed, 233 insertions(+), 212 deletions(-) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 9335815..d586c9f 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -193,7 +193,7 @@ class KDEgauss(object): def _eval_grid_fast(self, *args, **kwds): X = np.vstack(args) d, inc = X.shape - dx = X[:, 1] - X[:, 0] + #dx = X[:, 1] - X[:, 0] R = X.max(axis=-1)- X.min(axis=-1) t_star = (self.hs/R)**2 @@ -2401,24 +2401,24 @@ def _compute_qth_weighted_percentile(a, q, axis, out, method, weights, overwrite shape0 = a.shape if axis is None: - sorted = a.ravel() + sorted_ = a.ravel() else: taxes = range(a.ndim) taxes[-1], taxes[axis] = taxes[axis], taxes[-1] - sorted = np.transpose(a, taxes).reshape(-1, shape0[axis]) + sorted_ = np.transpose(a, taxes).reshape(-1, shape0[axis]) - ind = sorted.argsort(axis= -1) + ind = sorted_.argsort(axis= -1) if overwrite_input: - sorted.sort(axis= -1) + sorted_.sort(axis= -1) else: - sorted = np.sort(sorted, axis= -1) + sorted_ = np.sort(sorted_, axis= -1) w = np.atleast_1d(weights) n = len(w) w = w * n / w.sum() # Work on each column separately because of weight vector - m = sorted.shape[0] + m = sorted_.shape[0] nq = len(q) y = np.zeros((m, nq)) pk_fun = _PKDICT.get(method, 1) @@ -2426,8 +2426,8 @@ def _compute_qth_weighted_percentile(a, q, axis, out, method, weights, overwrite sortedW = w[ind[i]] # rearrange the weight according to ind k = sortedW.cumsum() # cumulative weight pk = pk_fun(k, sortedW, n) # different algorithm to compute percentile - # Interpolation between pk and sorted for given value of q - y[i] = np.interp(q, pk, sorted[i]) + # Interpolation between pk and sorted_ for given value of q + y[i] = np.interp(q, pk, sorted_[i]) if axis is None: return np.squeeze(y) else: @@ -2448,9 +2448,9 @@ _KDICT = {1:lambda p, n: p * (n - 1), 4:lambda p, n: p * (n + 1) - 1, 5:lambda p, n: p * (n + 1. / 3) - 2. / 3, 6:lambda p, n: p * (n + 1. / 4) - 5. / 8} -def _compute_qth_percentile(sorted, q, axis, out, method): +def _compute_qth_percentile(sorted_, q, axis, out, method): if not np.isscalar(q): - p = [_compute_qth_percentile(sorted, qi, axis, None, method) + p = [_compute_qth_percentile(sorted_, qi, axis, None, method) for qi in q] if out is not None: out.flat = p @@ -2460,8 +2460,8 @@ def _compute_qth_percentile(sorted, q, axis, out, method): if (q < 0) or (q > 1): raise ValueError, "percentile must be in the range [0,100]" - indexer = [slice(None)] * sorted.ndim - Nx = sorted.shape[axis] + indexer = [slice(None)] * sorted_.ndim + Nx = sorted_.shape[axis] k_fun = _KDICT.get(method, 1) index = np.clip(k_fun(q, Nx), 0, Nx - 1) i = int(index) @@ -2473,14 +2473,14 @@ def _compute_qth_percentile(sorted, q, axis, out, method): indexer[axis] = slice(i, i + 2) j = i + 1 weights1 = np.array([(j - index), (index - i)], float) - wshape = [1] * sorted.ndim + wshape = [1] * sorted_.ndim wshape[axis] = 2 weights1.shape = wshape sumval = weights1.sum() # Use add.reduce in both cases to coerce data type as well as # check and use out array. - return np.add.reduce(sorted[indexer] * weights1, axis=axis, out=out) / sumval + return np.add.reduce(sorted_[indexer] * weights1, axis=axis, out=out) / sumval def percentile(a, q, axis=None, out=None, overwrite_input=False, method=1, weights=None): """ @@ -2594,17 +2594,17 @@ def percentile(a, q, axis=None, out=None, overwrite_input=False, method=1, weigh return _compute_qth_weighted_percentile(a, q, axis, out, method, weights, overwrite_input) elif overwrite_input: if axis is None: - sorted = a.ravel() - sorted.sort() + sorted_ = a.ravel() + sorted_.sort() else: a.sort(axis=axis) - sorted = a + sorted_ = a else: - sorted = np.sort(a, axis=axis) + sorted_ = np.sort(a, axis=axis) if axis is None: axis = 0 - return _compute_qth_percentile(sorted, q, axis, out, method) + return _compute_qth_percentile(sorted_, q, axis, out, method) def iqrange(data, axis=None): ''' @@ -2736,7 +2736,7 @@ def gridcount(data, X, y=1): binx = np.asarray(np.floor((dat - xlo[:, newaxis]) / dx), dtype=int) w = dx.prod() - abs = np.abs + abs = np.abs #@ReservedAssignment if d == 1: x.shape = (-1,) c = np.asarray((acfun(binx, (x[binx + 1] - dat) * y, size=(inc, )) + @@ -2859,68 +2859,40 @@ def evar(y): return noisevar/M.mean()**2 #fun = lambda x : score_fun(x, S2, dcty2) Lopt = optimize.fminbound(score_fun, -38, 38, args=(S2, dcty2)) - M = 1-1./(1+10**Lopt*S2) + M = 1.0-1.0/(1+10**Lopt*S2) noisevar = (dcty2*M**2).mean() return noisevar -def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e-3, maxiter=100): +def smoothn(data, s=None, weight=None, robust=False, z0=None, tolz=1e-3, maxiter=100, fulloutput=False): ''' - SMOOTHN Robust spline smoothing for 1-D to N-D data. - SMOOTHN provides a fast, automatized and robust discretized smoothing - spline for data of any dimension. - - Z = SMOOTHN(Y) automatically smoothes the uniformly-sampled array Y. Y - can be any N-D noisy array (time series, images, 3D data,...). Non - finite data (NaN or Inf) are treated as missing values. - - Z = SMOOTHN(Y,S) smoothes the array Y using the smoothing parameter S. - S must be a real positive scalar. The larger S is, the smoother the - output will be. If the smoothing parameter S is omitted (see previous - option) or empty (i.e. S = []), it is automatically determined using - the generalized cross-validation (GCV) method. - - Z = SMOOTHN(Y,W) or Z = SMOOTHN(Y,W,S) specifies a weighting array W of - real positive values, that must have the same size as Y. Note that a - nil weight corresponds to a missing value. - - Robust smoothing - ---------------- - Z = SMOOTHN(...,'robust') carries out a robust smoothing that minimizes - the influence of outlying data. - - [Z,S] = SMOOTHN(...) also returns the calculated value for S so that - you can fine-tune the smoothing subsequently if needed. - - An iteration process is used in the presence of weighted and/or missing - values. Z = SMOOTHN(...,OPTION_NAME,OPTION_VALUE) smoothes with the - termination parameters specified by OPTION_NAME and OPTION_VALUE. They - can contain the following criteria: - ----------------- - TolZ: Termination tolerance on Z (default = 1e-3) - TolZ must be in ]0,1[ - MaxIter: Maximum number of iterations allowed (default = 100) - Initial: Initial value for the iterative process (default = - original data) - ----------------- - Syntax: [Z,...] = SMOOTHN(...,'MaxIter',500,'TolZ',1e-4,'Initial',Z0); - - [Z,S,EXITFLAG] = SMOOTHN(...) returns a boolean value EXITFLAG that - describes the exit condition of SMOOTHN: - 1 SMOOTHN converged. - 0 Maximum number of iterations was reached. - - Class Support - ------------- - Input array can be numeric or logical. The returned array is of class - double. + SMOOTHN fast and robust spline smoothing for 1-D to N-D data. - Notes - ----- - The N-D (inverse) discrete cosine transform functions DCTN and IDCTN are required. + + Parameters + ---------- + data : array like + uniformly-sampled data array to smooth. Non finite values (NaN or Inf) + are treated as missing values. + s : real positive scalar + smooting parameter. The larger S is, the smoother the output will be. + Default value is automatically determined using the generalized + cross-validation (GCV) method. + weight : string or array weights + weighting array of real positive values, that must have the same size as DATA. + Note that a zero weight corresponds to a missing value. + robust : bool + If true carry out a robust smoothing that minimizes the influence of outlying data. + tolz : real positive scalar + Termination tolerance on Z (default = 1e-3) + maxiter : scalar integer + Maximum number of iterations allowed (default = 100) + z0 : array-like + Initial value for the iterative process (default = original data) + + Returns + ------- + z : array like + smoothed data To be made ---------- @@ -2930,30 +2902,35 @@ def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e- --------- Garcia D, Robust smoothing of gridded data in one and higher dimensions with missing values. Computational Statistics & Data Analysis, 2010. - PDF download + http://www.biomecardio.com/pageshtm/publi/csda10.pdf Examples: -------- - 1-D example - x = linspace(0,100,2^8); - y = cos(x/10)+(x/50).^2 + randn(size(x))/10; - y([70 75 80]) = [5.5 5 6]; - z = smoothn(y); Regular smoothing - zr = smoothn(y,'robust'); Robust smoothing - subplot(121), plot(x,y,'r.',x,z,'k','LineWidth',2) - axis square, title('Regular smoothing') - subplot(122), plot(x,y,'r.',x,zr,'k','LineWidth',2) - axis square, title('Robust smoothing') + + 1-D example + >>> import matplotlib.pyplot as plt + >>> x = np.linspace(0,100,2**8) + >>> y = cos(x/10)+(x/50)**2 + np.random.randn(x.shape)/10 + >>> y[np.r_[70, 75, 80]] = np.array([5.5, 5, 6]) + >>> z = smoothn(y) # Regular smoothing + >>> zr = smoothn(y,robust=True) # Robust smoothing + >>> plt.subplot(121), + >>> h = plt.plot(x,y,'r.',x,z,'k',LineWidth=2) + >>> plt.title('Regular smoothing') + >>> plt.subplot(122) + >>> plt.plot(x,y,'r.',x,zr,'k',LineWidth=2) + >>> plt.title('Robust smoothing') 2-D example - xp = 0:.02:1; - [x,y] = meshgrid(xp); - f = exp(x+y) + sin((x-2*y)*3); - fn = f + randn(size(f))*0.5; - fs = smoothn(fn); - subplot(121), surf(xp,xp,fn), zlim([0 8]), axis square - subplot(122), surf(xp,xp,fs), zlim([0 8]), axis square + >>> xp = np.r_[0:1:.02] + >>> [x,y] = np.meshgrid(xp,xp) + >>> f = np.exp(x+y) + np.sin((x-2*y)*3); + >>> fn = f + randn(size(f))*0.5; + >>> fs = smoothn(fn); + >>> plt.subplot(121), + >>> plt.contourf(xp,xp,fn) + >>> plt.subplot(122) + >>> plt.contourf(xp,xp,fs) 2-D example with missing data n = 256; @@ -2980,7 +2957,8 @@ def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e- subplot(122), slice(x,y,z,v,xslice,yslice,zslice,'cubic') title('Smoothed data') - Cardioid + Cardioid + t = linspace(0,2*pi,1000); x = 2*cos(t).*(1-cos(t)) + randn(size(t))*0.1; y = 2*sin(t).*(1-cos(t)) + randn(size(t))*0.1; @@ -3024,7 +3002,9 @@ def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e- if weight is None: pass elif isinstance(weight, str): - weightstr = W.lower() + weightstr = weight.lower() + else: + W = weight # Weights. Zero weights are assigned to not finite values (Inf or NaN), # (Inf/NaN values = missing data). @@ -3047,10 +3027,10 @@ def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e- # penalized least squares process. d = y.ndim Lambda = np.zeros(sizy) - siz0 = (1,)*d + siz0 = [1,]*d for i in range(d): - siz0[i] = sizy(i) - Lambda = Lambda + np.cos(pi*np.arange(sizy(i))/sizy[i]).reshape(siz0) + siz0[i] = sizy[i] + Lambda = Lambda + np.cos(pi*np.arange(sizy[i])/sizy[i]).reshape(siz0) siz0[i] = 1 Lambda = -2*(d-Lambda) @@ -3066,66 +3046,70 @@ def unfinished_smoothn(data,s=None, weight=None, robust=False, z0=None, tolz=1e- N = (np.array(sizy)!=1).sum() # tensor rank of the y-array hMin = 1e-6; hMax = 0.99; - sMinBnd = (((1+sqrt(1+8*hMax**(2/N)))/4./hMax**(2/N))**2-1)/16 - sMaxBnd = (((1+sqrt(1+8*hMin**(2/N)))/4./hMin**(2/N))**2-1)/16 + sMinBnd = (((1+sqrt(1+8*hMax**(2./N)))/4./hMax**(2./N))**2-1)/16 + sMaxBnd = (((1+sqrt(1+8*hMin**(2./N)))/4./hMin**(2./N))**2-1)/16 # Initialize before iterating Wtot = W; # Initial conditions for z - if weight is None: - z = np.zeros(sizy) - else: + if isweighted: # With weighted/missing data # An initial guess is provided to ensure faster convergence. For that # purpose, a nearest neighbor interpolation followed by a coarse # smoothing are performed. if z0 is None: - z = InitialGuess(y,IsFinite); + z = InitialGuess(y,IsFinite) else: # an initial guess (z0) has been provided - z = z0 + z = z0 + else: + z = np.zeros(sizy) z0 = z y[1-IsFinite] = 0 # arbitrary values for missing y-data tol = 1 RobustIterativeProcess = True RobustStep = 1; - nit = 0; + # Error on p. Smoothness parameter s = 10^p errp = 0.1; - opt = optimset('TolX',errp); + # Relaxation factor RF: to speedup convergence RF = 1 + 0.75 if weight else 1.0 + norm = linalg.norm # Main iterative process - while RobustIterativeProcess + while RobustIterativeProcess: # "amount" of weights (see the function GCVscore) aow = Wtot.sum()/noe # 0 < aow <= 1 - - while tol>TolZ && nit tol=0 (no iteration) - tol = norm(z0(:)-z(:))/norm(z(:)) if weight else 0.0 - + tol = norm(z0.ravel()-z.ravel())/norm(z.ravel()) if isweighted else 0.0 + if tol<=tolz: + break z0 = z # re-initialization - end - exitflag = nit0.9: # aow = 1 means that all of the data are equally weighted - # very much faster: does not require any inverse DCT - RSS = norm(DCTy.ravel()*(Gamma.ravel()-1))**2 - else: - # take account of the weights to calculate RSS: - yhat = idctn(Gamma*DCTy) - RSS = norm(sqrt(Wtot[IsFinite])*(y[IsFinite]-yhat[IsFinite]))**2 - end - - TrH = Gamma.sum() - GCVscore = RSS/nof/(1.0-TrH/noe)**2 - return GCVScore - - - # Robust weights - def RobustWeights(r,I,h,wstr): - #weights for robust smoothing. - MAD = median(abs(r[I]-median(r[I]))) # median absolute deviation - u = abs(r/(1.4826*MAD)/sqrt(1-h)) # studentized residuals - if wstr == 'cauchy': - c = 2.385; - W = 1./(1+(u/c).^2); % Cauchy weights - elif wstr=='talworth': - c = 2.795; - W = u0.9: # aow = 1 means that all of the data are equally weighted + # very much faster: does not require any inverse DCT + RSS = linalg.norm(DCTy.ravel()*(Gamma.ravel()-1))**2 + else: + # take account of the weights to calculate RSS: + yhat = idctn(Gamma*DCTy) + RSS = linalg.norm(sqrt(Wtot[IsFinite])*(y[IsFinite]-yhat[IsFinite]))**2 + #end + + TrH = Gamma.sum() + GCVscore = RSS/nof/(1.0-TrH/noe)**2 + return GCVscore + + +# Robust weights +def RobustWeights(r,I,h,wstr): + #weights for robust smoothing. + MAD = np.median(abs(r[I]-np.median(r[I]))) # median absolute deviation + u = abs(r/(1.4826*MAD)/sqrt(1-h)) # studentized residuals + if wstr == 'cauchy': + c = 2.385; + W = 1./(1+(u/c)**2) # Cauchy weights + elif wstr=='talworth': + c = 2.795 + W = u