From f04fb14dc64d0759a37ddd6cd06948c1c2574886 Mon Sep 17 00:00:00 2001 From: Per A Brodtkorb Date: Sun, 28 Aug 2016 19:06:19 +0200 Subject: [PATCH] Simplified code. --- wafo/containers.py | 209 ++++++++++++++++++++++++--------------------- wafo/misc.py | 85 +++++++++--------- wafo/sg_filter.py | 58 ++++++------- 3 files changed, 178 insertions(+), 174 deletions(-) diff --git a/wafo/containers.py b/wafo/containers.py index 3fd55bd..fe3d65c 100644 --- a/wafo/containers.py +++ b/wafo/containers.py @@ -62,12 +62,11 @@ class PlotData(object): >>> h1 = d2() # Plot with confidence interval - >>> d3 = PlotData(np.sin(x), x) - >>> d3.children = [PlotData(np.vstack([np.sin(x)*0.9, - ... np.sin(x)*1.2]).T, x)] - >>> d3.plot_args_children = [':r'] - - >>> h = d3.plot() + >>> ci = PlotData(np.vstack([np.sin(x)*0.9, np.sin(x)*1.2]).T, x, + ... plot_args=[':r']) + >>> d3 = PlotData(np.sin(x), x, children=[ci]) + >>> h = d3.plot() # plot data, CI red dotted line + >>> h = d3.plot(plot_args_children=['b--']) # CI with blue dashed line ''' @@ -165,9 +164,10 @@ class PlotData(object): def integrate(self, a=None, b=None, **kwds): ''' >>> x = np.linspace(0,5,60) - >>> d = PlotData(np.sin(x), x) - >>> d.dataCI = np.vstack((d.data*.9,d.data*1.1)).T - >>> d.integrate(0,np.pi/2, return_ci=True) + >>> y = np.sin(x) + >>> ci = PlotData(np.vstack((y*.9, y*1.1)).T, x) + >>> d = PlotData(y, x, children=[ci]) + >>> d.integrate(0, np.pi/2, return_ci=True) array([ 0.99940055, 0.85543644, 1.04553343]) >>> np.allclose(d.integrate(0, 5, return_ci=True), ... d.integrate(return_ci=True)) @@ -187,12 +187,18 @@ class PlotData(object): b = x[-1] ix = np.flatnonzero((a < x) & (x < b)) xi = np.hstack((a, x.take(ix), b)) - fi = np.hstack((self.eval_points(a), self.data.take(ix), - self.eval_points(b))) + if self.data.ndim > 1: + fi = np.vstack((self.eval_points(a), + self.data[ix, :], + self.eval_points(b))).T + else: + fi = np.hstack((self.eval_points(a), self.data.take(ix), + self.eval_points(b))) res = fun(fi, xi, **kwds) if return_ci: - return np.hstack( - (res, fun(self.dataCI[ix, :].T, xi[1:-1], **kwds))) + res_ci = [child.integrate(a, b, method=method) + for child in self.children] + return np.hstack((res, res_ci)) return res def plot(self, *args, **kwds): @@ -205,15 +211,13 @@ class PlotData(object): if not plotflag and self.children is not None: axis.hold('on') tmp = [] - child_args = kwds.pop( - 'plot_args_children', - tuple( - self.plot_args_children)) + child_args = kwds.pop('plot_args_children', + tuple(self.plot_args_children)) child_kwds = dict(self.plot_kwds_children).copy() child_kwds.update(kwds.pop('plot_kwds_children', {})) child_kwds['axis'] = axis for child in self.children: - tmp1 = child(*child_args, **child_kwds) + tmp1 = child.plot(*child_args, **child_kwds) if tmp1 is not None: tmp.append(tmp1) if len(tmp) == 0: @@ -267,8 +271,8 @@ class AxisLabels: return self.__str__() def __str__(self): - return '%s\n%s\n%s\n%s\n' % ( - self.title, self.xlab, self.ylab, self.zlab) + txt = 'AxisLabels(title={}, xlab={}, ylab={}, zlab={})' + return txt.format(self.title, self.xlab, self.ylab, self.zlab) def copy(self): newcopy = empty_copy(self) @@ -318,11 +322,6 @@ class Plotter_1d(object): if plotmethod is None: plotmethod = 'plot' self.plotmethod = plotmethod - # self.plotbackend = plotbackend -# try: -# self.plotfun = getattr(plotbackend, plotmethod) -# except: -# pass def show(self, *args, **kwds): plt.show(*args, **kwds) @@ -350,72 +349,65 @@ class Plotter_1d(object): def _plot(self, axis, plotflag, wdata, *args, **kwds): x = wdata.args - data = transformdata(x, wdata.data, plotflag) + data = transformdata_1d(x, wdata.data, plotflag) dataCI = getattr(wdata, 'dataCI', ()) + if dataCI: + dataCI = transformdata_1d(x, dataCI, plotflag) h1 = plot1d(axis, x, data, dataCI, plotflag, *args, **kwds) return h1 __call__ = plot -def plot1d(axis, args, data, dataCI, plotflag, *varargin, **kwds): - - plottype = np.mod(plotflag, 10) - if plottype == 0: # No plotting - return [] - elif plottype == 1: - H = axis.plot(args, data, *varargin, **kwds) - elif plottype == 2: - H = axis.step(args, data, *varargin, **kwds) - elif plottype == 3: - H = axis.stem(args, data, *varargin, **kwds) - elif plottype == 4: - H = axis.errorbar(args, data, - yerr=[dataCI[:, 0] - data, - dataCI[:, 1] - data], - *varargin, **kwds) - elif plottype == 5: - H = axis.bar(args, data, *varargin, **kwds) - elif plottype == 6: - level = 0 - if np.isfinite(level): - H = axis.fill_between(args, data, level, *varargin, **kwds) - else: - H = axis.fill_between(args, data, *varargin, **kwds) - elif plottype == 7: - H = axis.plot(args, data, *varargin, **kwds) - H = axis.fill_between(args, dataCI[:, 0], dataCI[:, 1], - alpha=0.2, color='r') - +def set_plot_scale(axis, f_max, plotflag): scale = plotscale(plotflag) - logXscale = 'x' in scale - logYscale = 'y' in scale - logZscale = 'z' in scale - - if logXscale: + log_x_scale = 'x' in scale + log_y_scale = 'y' in scale + log_z_scale = 'z' in scale + if log_x_scale: axis.set(xscale='log') - if logYscale: + if log_y_scale: axis.set(yscale='log') - if logZscale: + if log_z_scale: axis.set(zscale='log') - - transFlag = np.mod(plotflag // 10, 10) - logScale = logXscale or logYscale or logZscale - if logScale or (transFlag == 5 and not logScale): + trans_flag = np.mod(plotflag // 10, 10) + log_scale = log_x_scale or log_y_scale or log_z_scale + if log_scale or (trans_flag == 5 and not log_scale): ax = list(axis.axis()) - fmax1 = data.max() - if transFlag == 5 and not logScale: - ax[3] = 11 * np.log10(fmax1) + if trans_flag == 8 and not log_scale: + ax[3] = 11 * np.log10(f_max) ax[2] = ax[3] - 40 else: - ax[3] = 1.15 * fmax1 + ax[3] = 1.15 * f_max ax[2] = ax[3] * 1e-4 - axis.axis(ax) - if np.any(dataCI) and plottype < 3: - axis.hold(True) - plot1d(axis, args, dataCI, (), plotflag, 'r--') - return H + +def plot1d(axis, args, data, dataCI, plotflag, *varargin, **kwds): + h = [] + plottype = np.mod(plotflag, 10) + if plottype == 0: # No plotting + return h + fun = {1: 'plot', 2: 'step', 3: 'stem', 5: 'bar'}.get(plottype) + if fun is not None: + plotfun = getattr(axis, fun) + h.extend(plotfun(args, data, *varargin, **kwds)) + if np.any(dataCI) and plottype < 3: + axis.hold(True) + h.extend(plotfun(args, dataCI, 'r--')) + elif plottype == 4: + h = axis.errorbar(args, data, + yerr=[dataCI[:, 0] - data, + dataCI[:, 1] - data], + *varargin, **kwds) + elif plottype == 6: + h = axis.fill_between(args, data, 0, *varargin, **kwds) + elif plottype == 7: + h = axis.plot(args, data, *varargin, **kwds) + h.extend(axis.fill_between(args, dataCI[:, 0], dataCI[:, 1], + alpha=0.2, color='r')) + fmax1 = data.max() + set_plot_scale(axis, fmax1, plotflag) + return h def plotscale(plotflag): @@ -462,6 +454,10 @@ def plotscale(plotflag): 'ylog' >>> plotscale(10000) 'zlog' + >>> plotscale(1100) + 'xylog' + >>> plotscale(11100) + 'xyzlog' See also --------- @@ -480,25 +476,41 @@ def plotscale(plotflag): return scales[scaleId] -def transformdata(x, f, plotflag): - transFlag = np.mod(plotflag // 10, 10) - if transFlag == 0: - data = f - elif transFlag == 1: - data = 1 - f - elif transFlag == 2: - data = cumtrapz(f, x) - elif transFlag == 3: - data = 1 - cumtrapz(f, x) - if transFlag in (4, 5): - if transFlag == 4: - data = -np.log1p(-cumtrapz(f, x)) - else: - if any(f < 0): - raise ValueError('Invalid plotflag: Data or dataCI is ' + - 'negative, but must be positive') - data = 10 * np.log10(f) - return data +def plotflag2plottype_1d(plotflag): + plottype = np.mod(plotflag, 10) + return ['', 'plot', 'step', 'stem', 'errorbar', 'bar'][plottype] + + +def plotflag2transform_id(plotflag): + transform_id = np.mod(plotflag // 10, 10) + return ['f', '1-f', + 'cumtrapz(f)', '1-cumtrapz(f)', + 'log(f)', 'log(1-f)' + 'log(cumtrapz(f))', 'log(cumtrapz(f))', + 'log10(f)'][transform_id] + + +def transform_id2plotflag2(transform_id): + return {'': 0, 'None': 0, 'f': 0, '1-f': 1, + 'cumtrapz(f)': 2, '1-cumtrapz(f)': 3, + 'log(f)': 4, 'log(1-f)': 5, + 'log(cumtrapz(f))': 6, 'log(1-cumtrapz(f))': 7, + '10log10(f)': 8}[transform_id] * 10 + + +def transformdata_1d(x, f, plotflag): + transform_id = np.mod(plotflag // 10, 10) + transform = [lambda f, x: f, + lambda f, x: 1 - f, + lambda f, x: cumtrapz(f, x), + lambda f, x: 1 - cumtrapz(f, x), + lambda f, x: np.log(f), + lambda f, x: np.log1p(-f), + lambda f, x: np.log(cumtrapz(f, x)), + lambda f, x: np.log1p(-cumtrapz(f, x)), + lambda f, x: 10*np.log10(f) + ][transform_id] + return transform(f, x) class Plotter_2d(Plotter_1d): @@ -589,14 +601,15 @@ def test_plotdata(): x = np.linspace(0, np.pi, 9) xi = np.linspace(0, np.pi, 4*9) - d = PlotData(np.sin(x)/2, x, xlab='x', ylab='sin', title='sinus', + d = PlotData(np.sin(x)/2, x, dataCI=[], xlab='x', ylab='sin', title='sinus', plot_args=['r.']) di = PlotData(d.eval_points(xi, method='cubic'), xi) unused_hi = di.plot() unused_h = d.plot() f = di.to_cdf() + for i in range(4): - _ = f.plot(plotflag=i) + _ = di.plot(plotflag=i) d.show('hold') @@ -607,5 +620,5 @@ def test_docstrings(): if __name__ == '__main__': - test_docstrings() - # test_plotdata() + # test_docstrings() + test_plotdata() diff --git a/wafo/misc.py b/wafo/misc.py index a29c884..0ade928 100644 --- a/wafo/misc.py +++ b/wafo/misc.py @@ -7,7 +7,7 @@ import fractions import numpy as np from numpy import ( meshgrid, - abs, amax, any, logical_and, arange, linspace, atleast_1d, + amax, logical_and, arange, linspace, atleast_1d, asarray, ceil, floor, frexp, hypot, sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, finfo, inf, pi, interp, isscalar, zeros, ones, linalg, @@ -39,7 +39,7 @@ __all__ = ['now', 'spaceline', 'narg_smallest', 'args_flat', 'is_numlike', 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart', 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', - 'lazyselect' + 'lazyselect', 'piecewise', 'valarray', 'check_random_state'] @@ -269,8 +269,8 @@ def lazyselect(condlist, choicelist, arrays, default=0): arrays = np.broadcast_arrays(*arrays) tcode = np.mintypecode([a.dtype.char for a in arrays]) out = valarray(np.shape(arrays[0]), value=default, typecode=tcode) - for index in range(len(condlist)): - func, cond = choicelist[index], condlist[index] + for index, cond in enumerate(condlist): + func = choicelist[index] if np.all(cond is False): continue cond, _ = np.broadcast_arrays(cond, arrays[0]) @@ -669,8 +669,8 @@ class Bunch(object): self.__dict__.update(kwargs) -def printf(format, *args): # @ReservedAssignment - sys.stdout.write(format % args) +def printf(format_, *args): + sys.stdout.write(format_ % args) def sub_dict_select(somedict, somekeys): @@ -1612,8 +1612,8 @@ def mctp2rfc(fmM, fMm=None): # end # end fx = 0.0 - if (max(abs(e)) > 1e-6 and - max(abs(NN)) > 1e-6 * max(MA[0], mA[0])): + if (max(np.abs(e)) > 1e-6 and + max(np.abs(NN)) > 1e-6 * max(MA[0], mA[0])): PMm = _get_PMm(AA1, MA, nA) A = PMm @@ -1722,7 +1722,7 @@ def _findrfc2(y, h, method=0): t1, y1 = (t0, y0) if y0 > yi else (ti, yi) # Update y if y0 is a turning point - if abs(z0 - z1) == 2: + if np.abs(z0 - z1) == 2: j += 1 t[j] = t0 @@ -1731,7 +1731,7 @@ def _findrfc2(y, h, method=0): # end # Update y if last y0 is greater than (or equal) threshold - if cmpfun1(h, abs(y0 - y[t[j]])): + if cmpfun1(h, np.abs(y0 - y[t[j]])): j += 1 t[j] = t0 return t[:j + 1] @@ -1965,28 +1965,20 @@ def findtc(x_in, v=None, kind=None): first_is_down_crossing = (x[v_ind[0]] > x[v_ind[0] + 1]) if first_is_down_crossing: - for i in range(n_tc): - # trough - j = 2 * i - ind[j] = x[v_ind[j] + 1:v_ind[j + 1] + 1].argmin() - # crest - ind[j + 1] = x[v_ind[j + 1] + 1:v_ind[j + 2] + 1].argmax() - - if (2 * n_tc + 1 < n_c) and (kind in (None, 'tw')): - # trough - ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1].argmin() - - else: # the first is a up-crossing - for i in range(n_tc): - # crest - j = 2 * i - ind[j] = x[v_ind[j] + 1:v_ind[j + 1] + 1].argmax() - # trough - ind[j + 1] = x[v_ind[j + 1] + 1:v_ind[j + 2] + 1].argmin() - - if (2 * n_tc + 1 < n_c) and (kind in (None, 'cw')): - # crest - ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1].argmax() + f1, f2 = np.argmin, np.argmax + else: + f1, f2 = np.argmax, np.argmin + + for i in range(n_tc): + # trough or crest + j = 2 * i + ind[j] = f1(x[v_ind[j] + 1:v_ind[j + 1] + 1]) + # crest or trough + ind[j + 1] = f2(x[v_ind[j + 1] + 1:v_ind[j + 2] + 1]) + + if (2 * n_tc + 1 < n_c) and (kind in (None, 'tw', 'cw')): + # trough or crest + ind[n_c - 2] = f1(x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1]) return v_ind[:n_c - 1] + ind + 1, v_ind @@ -2081,7 +2073,7 @@ def findoutliers(x, zcrit=0.0, dcrit=None, ddcrit=None, verbose=False): def _find_consecutive_equal_values(dxn, zcrit): - mask_small = (abs(dxn) <= zcrit) + mask_small = (np.abs(dxn) <= zcrit) i_small = np.flatnonzero(mask_small) if verbose: if zcrit == 0.: @@ -2130,7 +2122,7 @@ def findoutliers(x, zcrit=0.0, dcrit=None, ddcrit=None, verbose=False): indg, = nonzero(indg) if verbose: - print('Found the total of %d spurious points' % ind.size) + print('Found the total of %d spurious points' % np.size(ind)) return ind, indg @@ -2174,8 +2166,7 @@ def common_shape(*args, ** kwds): ''' shape = kwds.get('shape') x0 = 1 if shape is None else np.ones(shape) - x1 = np.broadcast(x0, *args) - return tuple(x1.shape) + return tuple(np.broadcast(x0, *args).shape) def argsreduce(condition, * args): @@ -2275,15 +2266,15 @@ def stirlerr(n): n500 = 500 < n1 y[n500] = (S0 - S1 / nn[n500]) / n1[n500] n80 = logical_and(80 < n1, n1 <= 500) - if any(n80): + if np.any(n80): y[n80] = (S0 - (S1 - S2 / nn[n80]) / nn[n80]) / n1[n80] n35 = logical_and(35 < n1, n1 <= 80) - if any(n35): + if np.any(n35): nn35 = nn[n35] y[n35] = (S0 - (S1 - (S2 - S3 / nn35) / nn35) / nn35) / n1[n35] n15 = logical_and(15 < n1, n1 <= 35) - if any(n15): + if np.any(n15): nn15 = nn[n15] y[n15] = ( S0 - (S1 - (S2 - (S3 - S4 / nn15) / nn15) / nn15) / nn15) / n1[n15] @@ -2426,7 +2417,7 @@ def binomln(z, w): Example ------- - >>> abs(binomln(3,2)- 1.09861229)<1e-7 + >>> np.abs(binomln(3,2)- 1.09861229)<1e-7 array([ True], dtype=bool) See also @@ -2463,7 +2454,7 @@ def betaloge(z, w): Example ------- >>> import wafo.misc as wm - >>> abs(wm.betaloge(3,2)+2.48490665)<1e-7 + >>> np.abs(wm.betaloge(3,2)+2.48490665)<1e-7 array([ True], dtype=bool) See also @@ -2545,7 +2536,7 @@ def nextpow2(x): if (t > 1): f, n = frexp(t) else: - f, n = frexp(abs(x)) + f, n = frexp(np.abs(x)) if (f == 0.5): n = n - 1 @@ -2615,7 +2606,7 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5): x = linspace(a, b, n) y = fun(x) y00 = interp(x, x0, y0) - err = 0.5 * amax(abs((y00 - y) / (abs(y00 + y) + _TINY))) + err = 0.5 * amax(np.abs((y00 - y) / (np.abs(y00 + y) + _TINY))) return x, y @@ -2644,7 +2635,7 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5): fy = fun(y) fy0 = interp(y, x, fx) - erri = 0.5 * (abs((fy0 - fy) / (abs(fy0 + fy) + _TINY))) + erri = 0.5 * (np.abs((fy0 - fy) / (np.abs(fy0 + fy) + _TINY))) err = erri.max() @@ -2794,7 +2785,7 @@ def trangood(x, f, min_n=None, min_x=None, max_x=None, max_n=inf): xn = xo[-1] x0 = xo[0] L = float(xn - x0) - if ((nf < min_n) or (max_n < nf) or any(abs(ddx) > 10 * _EPS * (L))): + if ((nf < min_n) or (max_n < nf) or np.any(np.abs(ddx) > 10 * _EPS * (L))): # pab 07.01.2001: Always choose the stepsize df so that # it is an exactly representable number. # This is important when calculating numerical derivatives and is @@ -3145,10 +3136,10 @@ def num2pistr(x, n=3, numerator_max=10, denominator_max=10): True ''' def _denominator_text(den): - return '' if abs(den) == 1 else '/%d' % den + return '' if np.abs(den) == 1 else '/%d' % den def _numerator_text(num): - if abs(num) == 1: + if np.abs(num) == 1: return '-' if num == -1 else '' return '{:d}'.format(num) frac = fractions.Fraction.from_float(x / pi).limit_denominator(int(1e+13)) diff --git a/wafo/sg_filter.py b/wafo/sg_filter.py index ad05428..22a6533 100644 --- a/wafo/sg_filter.py +++ b/wafo/sg_filter.py @@ -2,7 +2,7 @@ from __future__ import absolute_import, division import numpy as np # from math import pow # from numpy import zeros,dot -from numpy import (pi, abs, size, convolve, linalg, concatenate, sqrt) +from numpy import (pi, convolve, linalg, concatenate, sqrt) from scipy.sparse import spdiags from scipy.sparse.linalg import spsolve, expm from scipy.signal import medfilt @@ -132,14 +132,14 @@ class SavitzkyGolay(object): def smooth_last(self, signal, k=0): coeff = self._coeff - n = size(coeff - 1) // 2 + n = np.size((coeff - 1) // 2 y = np.squeeze(signal) if n == 0: return y if y.ndim > 1: coeff.shape = (-1, 1) - first_vals = y[0] - abs(y[n:0:-1] - y[0]) - last_vals = y[-1] + abs(y[-2:-n - 2:-1] - y[-1]) + first_vals = y[0] - np.abs(y[n:0:-1] - y[0]) + last_vals = y[-1] + np.abs(y[-2:-n - 2:-1] - y[-1]) y = concatenate((first_vals, y, last_vals)) return (y[-2 * n - 1 - k:-k] * coeff).sum(axis=0) @@ -147,9 +147,8 @@ class SavitzkyGolay(object): return self.smooth(signal) def smooth(self, signal): - x = np.asarray(signal) - if x.dtype != np.float64 and x.dtype != np.float32: - x = x.astype(np.float64) + dtype = np.result_type(signal, np.float) + x = np.asarray(signal, dtype=dtype) coeffs = self._coeff mode, axis = self.mode, self.axis @@ -180,13 +179,13 @@ class SavitzkyGolay(object): the smoothed signal (or it's n-th derivative). """ coeff = self._coeff - n = size(coeff - 1) // 2 + n = np.size(coeff - 1) // 2 y = np.squeeze(signal) if n == 0: return y if pad: - first_vals = y[0] - abs(y[n:0:-1] - y[0]) - last_vals = y[-1] + abs(y[-2:-n - 2:-1] - y[-1]) + first_vals = y[0] - np.abs(y[n:0:-1] - y[0]) + last_vals = y[-1] + np.abs(y[-2:-n - 2:-1] - y[-1]) y = concatenate((first_vals, y, last_vals)) n *= 2 d = y.ndim @@ -424,8 +423,8 @@ class _Filter(object): @staticmethod def _studentized_residuals(r, I, h): - median_abs_deviation = np.median(abs(r[I] - np.median(r[I]))) - return abs(r / (1.4826 * median_abs_deviation) / sqrt(1 - h)) + median_abs_deviation = np.median(np.abs(r[I] - np.median(r[I]))) + return np.abs(r / (1.4826 * median_abs_deviation) / sqrt(1 - h)) def robust_weights(self, r, I, h): """Return weights for robust smoothing.""" @@ -458,10 +457,10 @@ class _Filter(object): def check_smooth_parameter(self, s): if self.auto_smooth: - if abs(np.log10(s) - np.log10(self.s_min)) < self.errp: + if np.abs(np.log10(s) - np.log10(self.s_min)) < self.errp: warnings.warn('''s = %g: the lower bound for s has been reached. Put s as an input variable if required.''' % s) - elif abs(np.log10(s) - np.log10(self.s_max)) < self.errp: + elif np.abs(np.log10(s) - np.log10(self.s_max)) < self.errp: warnings.warn('''s = %g: the Upper bound for s has been reached. Put s as an input variable if required.''' % s) @@ -723,8 +722,8 @@ def test_smoothn_1d(): 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), - unused_h = plt.plot(x, y, 'r.', x, z, 'k', linewidth=2) + _h0 = 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) @@ -735,18 +734,18 @@ def test_smoothn_1d(): def test_smoothn_2d(): # import mayavi.mlab as plt - xp = np.r_[0:1:.02] + xp = np.r_[0:1:0.02] [x, y] = np.meshgrid(xp, xp) f = np.exp(x + y) + np.sin((x - 2 * y) * 3) fn = f + np.random.randn(*f.shape) * 0.5 - fs, s = smoothn(fn, fulloutput=True) # @UnusedVariable + _fs, s = smoothn(fn, fulloutput=True) fs2 = smoothn(fn, s=2 * s) - plt.subplot(131), - plt.contourf(xp, xp, fn) - plt.subplot(132), - plt.contourf(xp, xp, fs2) - plt.subplot(133), - plt.contourf(xp, xp, f) + _h = plt.subplot(131), + _h = plt.contourf(xp, xp, fn) + _h = plt.subplot(132), + _h = plt.contourf(xp, xp, fs2) + _h = plt.subplot(133), + _h = plt.contourf(xp, xp, f) plt.show('hold') @@ -1048,7 +1047,7 @@ class Kalman(object): def _update_covariance(self, P, K): return P - np.dot(K, np.dot(self.H, P)) - return np.dot(np.eye(len(P)) - K * self.H, P) + # return np.dot(np.eye(len(P)) - K * self.H, P) def _filter_main(self, z, u): ''' This is the code which implements the discrete Kalman filter: @@ -1187,7 +1186,7 @@ def test_kalman_sine(): x = np.zeros((n, m)) for i, zi in enumerate(z): - x[i] = filt(zi, u=0).ravel() + x[i] = np.ravel(filt(zi, u=0)) _hz = plt.plot(z, 'r.', label='observations') # a-posteriori state estimates: @@ -1321,7 +1320,8 @@ class HampelFilter(object): self.adaptive = adaptive self.fulloutput = fulloutput - def _check(self, dx): + @staticmethod + def _check(dx): if not np.isscalar(dx): raise ValueError('DX must be a scalar.') if dx < 0: @@ -1357,8 +1357,8 @@ class HampelFilter(object): while S0Rel > self.adaptive: Y0Tmp[j], S0Tmp[j] = localwindow(X, Y, DXTmp[j], i) if j > 0: - S0Rel = abs((S0Tmp[j - 1] - S0Tmp[j]) / - (S0Tmp[j - 1] + S0Tmp[j]) / 2) + S0Rel = np.abs((S0Tmp[j - 1] - S0Tmp[j]) / + (S0Tmp[j - 1] + S0Tmp[j]) / 2) j += 1 Y0[i] = Y0Tmp[j - 2]