danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Simplified code.

Browse Source
master
Per A Brodtkorb 8 years ago
parent 29e1aed0ca
commit f04fb14dc6
3 changed files with 178 additions and 174 deletions
Show Stats Download Patch File Download Diff File

199
wafo/containers.py
Unescape Escape View File

@ -62,12 +62,11 @@ class PlotData(object):
>>> h1 = d2() >>> h1 = d2()
# Plot with confidence interval # Plot with confidence interval
>>> d3 = PlotData(np.sin(x), x) >>> ci = PlotData(np.vstack([np.sin(x)*0.9, np.sin(x)*1.2]).T, x,
>>> d3.children = [PlotData(np.vstack([np.sin(x)*0.9, ... plot_args=[':r'])
... np.sin(x)*1.2]).T, x)] >>> d3 = PlotData(np.sin(x), x, children=[ci])
>>> d3.plot_args_children = [':r'] >>> h = d3.plot() # plot data, CI red dotted line
>>> h = d3.plot(plot_args_children=['b--']) # CI with blue dashed line
>>> h = d3.plot()
''' '''
@ -165,8 +164,9 @@ class PlotData(object):
def integrate(self, a=None, b=None, **kwds): def integrate(self, a=None, b=None, **kwds):
''' '''
>>> x = np.linspace(0,5,60) >>> x = np.linspace(0,5,60)
>>> d = PlotData(np.sin(x), x) >>> y = np.sin(x)
>>> d.dataCI = np.vstack((d.data*.9,d.data*1.1)).T >>> 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) >>> d.integrate(0, np.pi/2, return_ci=True)
array([ 0.99940055, 0.85543644, 1.04553343]) array([ 0.99940055, 0.85543644, 1.04553343])
>>> np.allclose(d.integrate(0, 5, return_ci=True), >>> np.allclose(d.integrate(0, 5, return_ci=True),
@ -187,12 +187,18 @@ class PlotData(object):
b = x[-1] b = x[-1]
ix = np.flatnonzero((a < x) & (x < b)) ix = np.flatnonzero((a < x) & (x < b))
xi = np.hstack((a, x.take(ix), b)) xi = np.hstack((a, x.take(ix), 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), fi = np.hstack((self.eval_points(a), self.data.take(ix),
self.eval_points(b))) self.eval_points(b)))
res = fun(fi, xi, **kwds) res = fun(fi, xi, **kwds)
if return_ci: if return_ci:
return np.hstack( res_ci = [child.integrate(a, b, method=method)
(res, fun(self.dataCI[ix, :].T, xi[1:-1], **kwds))) for child in self.children]
return np.hstack((res, res_ci))
return res return res
def plot(self, *args, **kwds): def plot(self, *args, **kwds):
@ -205,15 +211,13 @@ class PlotData(object):
if not plotflag and self.children is not None: if not plotflag and self.children is not None:
axis.hold('on') axis.hold('on')
tmp = [] tmp = []
child_args = kwds.pop( child_args = kwds.pop('plot_args_children',
'plot_args_children', tuple(self.plot_args_children))
tuple(
self.plot_args_children))
child_kwds = dict(self.plot_kwds_children).copy() child_kwds = dict(self.plot_kwds_children).copy()
child_kwds.update(kwds.pop('plot_kwds_children', {})) child_kwds.update(kwds.pop('plot_kwds_children', {}))
child_kwds['axis'] = axis child_kwds['axis'] = axis
for child in self.children: for child in self.children:
tmp1 = child(*child_args, **child_kwds) tmp1 = child.plot(*child_args, **child_kwds)
if tmp1 is not None: if tmp1 is not None:
tmp.append(tmp1) tmp.append(tmp1)
if len(tmp) == 0: if len(tmp) == 0:
@ -267,8 +271,8 @@ class AxisLabels:
return self.__str__() return self.__str__()
def __str__(self): def __str__(self):
return '%s\n%s\n%s\n%s\n' % ( txt = 'AxisLabels(title={}, xlab={}, ylab={}, zlab={})'
self.title, self.xlab, self.ylab, self.zlab) return txt.format(self.title, self.xlab, self.ylab, self.zlab)
def copy(self): def copy(self):
newcopy = empty_copy(self) newcopy = empty_copy(self)
@ -318,11 +322,6 @@ class Plotter_1d(object):
if plotmethod is None: if plotmethod is None:
plotmethod = 'plot' plotmethod = 'plot'
self.plotmethod = plotmethod self.plotmethod = plotmethod
# self.plotbackend = plotbackend
# try:
# self.plotfun = getattr(plotbackend, plotmethod)
# except:
# pass
def show(self, *args, **kwds): def show(self, *args, **kwds):
plt.show(*args, **kwds) plt.show(*args, **kwds)
@ -350,72 +349,65 @@ class Plotter_1d(object):
def _plot(self, axis, plotflag, wdata, *args, **kwds): def _plot(self, axis, plotflag, wdata, *args, **kwds):
x = wdata.args x = wdata.args
data = transformdata(x, wdata.data, plotflag) data = transformdata_1d(x, wdata.data, plotflag)
dataCI = getattr(wdata, 'dataCI', ()) dataCI = getattr(wdata, 'dataCI', ())
if dataCI:
dataCI = transformdata_1d(x, dataCI, plotflag)
h1 = plot1d(axis, x, data, dataCI, plotflag, *args, **kwds) h1 = plot1d(axis, x, data, dataCI, plotflag, *args, **kwds)
return h1 return h1
__call__ = plot __call__ = plot
def plot1d(axis, args, data, dataCI, plotflag, *varargin, **kwds): def set_plot_scale(axis, f_max, plotflag):
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')
scale = plotscale(plotflag) scale = plotscale(plotflag)
logXscale = 'x' in scale log_x_scale = 'x' in scale
logYscale = 'y' in scale log_y_scale = 'y' in scale
logZscale = 'z' in scale log_z_scale = 'z' in scale
if log_x_scale:
if logXscale:
axis.set(xscale='log') axis.set(xscale='log')
if logYscale: if log_y_scale:
axis.set(yscale='log') axis.set(yscale='log')
if logZscale: if log_z_scale:
axis.set(zscale='log') axis.set(zscale='log')
trans_flag = np.mod(plotflag // 10, 10)
transFlag = np.mod(plotflag // 10, 10) log_scale = log_x_scale or log_y_scale or log_z_scale
logScale = logXscale or logYscale or logZscale if log_scale or (trans_flag == 5 and not log_scale):
if logScale or (transFlag == 5 and not logScale):
ax = list(axis.axis()) ax = list(axis.axis())
fmax1 = data.max() if trans_flag == 8 and not log_scale:
if transFlag == 5 and not logScale: ax[3] = 11 * np.log10(f_max)
ax[3] = 11 * np.log10(fmax1)
ax[2] = ax[3] - 40 ax[2] = ax[3] - 40
else: else:
ax[3] = 1.15 * fmax1 ax[3] = 1.15 * f_max
ax[2] = ax[3] * 1e-4 ax[2] = ax[3] * 1e-4
axis.axis(ax) axis.axis(ax)
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: if np.any(dataCI) and plottype < 3:
axis.hold(True) axis.hold(True)
plot1d(axis, args, dataCI, (), plotflag, 'r--') h.extend(plotfun(args, dataCI, 'r--'))
return H 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): def plotscale(plotflag):
@ -462,6 +454,10 @@ def plotscale(plotflag):
'ylog' 'ylog'
>>> plotscale(10000) >>> plotscale(10000)
'zlog' 'zlog'
>>> plotscale(1100)
'xylog'
>>> plotscale(11100)
'xyzlog'
See also See also
--------- ---------
@ -480,25 +476,41 @@ def plotscale(plotflag):
return scales[scaleId] return scales[scaleId]
def transformdata(x, f, plotflag): def plotflag2plottype_1d(plotflag):
transFlag = np.mod(plotflag // 10, 10) plottype = np.mod(plotflag, 10)
if transFlag == 0: return ['', 'plot', 'step', 'stem', 'errorbar', 'bar'][plottype]
data = f
elif transFlag == 1:
data = 1 - f def plotflag2transform_id(plotflag):
elif transFlag == 2: transform_id = np.mod(plotflag // 10, 10)
data = cumtrapz(f, x) return ['f', '1-f',
elif transFlag == 3: 'cumtrapz(f)', '1-cumtrapz(f)',
data = 1 - cumtrapz(f, x) 'log(f)', 'log(1-f)'
if transFlag in (4, 5): 'log(cumtrapz(f))', 'log(cumtrapz(f))',
if transFlag == 4: 'log10(f)'][transform_id]
data = -np.log1p(-cumtrapz(f, x))
else:
if any(f < 0): def transform_id2plotflag2(transform_id):
raise ValueError('Invalid plotflag: Data or dataCI is ' + return {'': 0, 'None': 0, 'f': 0, '1-f': 1,
'negative, but must be positive') 'cumtrapz(f)': 2, '1-cumtrapz(f)': 3,
data = 10 * np.log10(f) 'log(f)': 4, 'log(1-f)': 5,
return data '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): class Plotter_2d(Plotter_1d):
@ -589,14 +601,15 @@ def test_plotdata():
x = np.linspace(0, np.pi, 9) x = np.linspace(0, np.pi, 9)
xi = np.linspace(0, np.pi, 4*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.']) plot_args=['r.'])
di = PlotData(d.eval_points(xi, method='cubic'), xi) di = PlotData(d.eval_points(xi, method='cubic'), xi)
unused_hi = di.plot() unused_hi = di.plot()
unused_h = d.plot() unused_h = d.plot()
f = di.to_cdf() f = di.to_cdf()
for i in range(4): for i in range(4):
_ = f.plot(plotflag=i) _ = di.plot(plotflag=i)
d.show('hold') d.show('hold')
@ -607,5 +620,5 @@ def test_docstrings():
if __name__ == '__main__': if __name__ == '__main__':
test_docstrings() # test_docstrings()
# test_plotdata() test_plotdata()

77
wafo/misc.py
Unescape Escape View File

@ -7,7 +7,7 @@ import fractions
import numpy as np import numpy as np
from numpy import ( from numpy import (
meshgrid, meshgrid,
abs, amax, any, logical_and, arange, linspace, atleast_1d, amax, logical_and, arange, linspace, atleast_1d,
asarray, ceil, floor, frexp, hypot, asarray, ceil, floor, frexp, hypot,
sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff, sqrt, arctan2, sin, cos, exp, log, log1p, mod, diff,
finfo, inf, pi, interp, isscalar, zeros, ones, linalg, 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', 'betaloge', 'gravity', 'nextpow2', 'discretize', 'polar2cart',
'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc', 'cart2polar', 'meshgrid', 'ndgrid', 'trangood', 'tranproc',
'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere', 'plot_histgrm', 'num2pistr', 'test_docstrings', 'lazywhere',
'lazyselect' 'lazyselect',
'piecewise', 'piecewise',
'valarray', 'check_random_state'] 'valarray', 'check_random_state']
@ -269,8 +269,8 @@ def lazyselect(condlist, choicelist, arrays, default=0):
arrays = np.broadcast_arrays(*arrays) arrays = np.broadcast_arrays(*arrays)
tcode = np.mintypecode([a.dtype.char for a in arrays]) tcode = np.mintypecode([a.dtype.char for a in arrays])
out = valarray(np.shape(arrays[0]), value=default, typecode=tcode) out = valarray(np.shape(arrays[0]), value=default, typecode=tcode)
for index in range(len(condlist)): for index, cond in enumerate(condlist):
func, cond = choicelist[index], condlist[index] func = choicelist[index]
if np.all(cond is False): if np.all(cond is False):
continue continue
cond, _ = np.broadcast_arrays(cond, arrays[0]) cond, _ = np.broadcast_arrays(cond, arrays[0])
@ -669,8 +669,8 @@ class Bunch(object):
self.__dict__.update(kwargs) self.__dict__.update(kwargs)
def printf(format, *args): # @ReservedAssignment def printf(format_, *args):
sys.stdout.write(format % args) sys.stdout.write(format_ % args)
def sub_dict_select(somedict, somekeys): def sub_dict_select(somedict, somekeys):
@ -1612,8 +1612,8 @@ def mctp2rfc(fmM, fMm=None):
# end # end
# end # end
fx = 0.0 fx = 0.0
if (max(abs(e)) > 1e-6 and if (max(np.abs(e)) > 1e-6 and
max(abs(NN)) > 1e-6 * max(MA[0], mA[0])): max(np.abs(NN)) > 1e-6 * max(MA[0], mA[0])):
PMm = _get_PMm(AA1, MA, nA) PMm = _get_PMm(AA1, MA, nA)
A = PMm A = PMm
@ -1722,7 +1722,7 @@ def _findrfc2(y, h, method=0):
t1, y1 = (t0, y0) if y0 > yi else (ti, yi) t1, y1 = (t0, y0) if y0 > yi else (ti, yi)
# Update y if y0 is a turning point # Update y if y0 is a turning point
if abs(z0 - z1) == 2: if np.abs(z0 - z1) == 2:
j += 1 j += 1
t[j] = t0 t[j] = t0
@ -1731,7 +1731,7 @@ def _findrfc2(y, h, method=0):
# end # end
# Update y if last y0 is greater than (or equal) threshold # 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 j += 1
t[j] = t0 t[j] = t0
return t[:j + 1] 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]) first_is_down_crossing = (x[v_ind[0]] > x[v_ind[0] + 1])
if first_is_down_crossing: if first_is_down_crossing:
for i in range(n_tc): f1, f2 = np.argmin, np.argmax
# trough else:
j = 2 * i f1, f2 = np.argmax, np.argmin
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): for i in range(n_tc):
# crest # trough or crest
j = 2 * i j = 2 * i
ind[j] = x[v_ind[j] + 1:v_ind[j + 1] + 1].argmax() ind[j] = f1(x[v_ind[j] + 1:v_ind[j + 1] + 1])
# trough # crest or trough
ind[j + 1] = x[v_ind[j + 1] + 1:v_ind[j + 2] + 1].argmin() 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, 'cw')): if (2 * n_tc + 1 < n_c) and (kind in (None, 'tw', 'cw')):
# crest # trough or crest
ind[n_c - 2] = x[v_ind[n_c - 2] + 1:v_ind[n_c - 1]+1].argmax() 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 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): def _find_consecutive_equal_values(dxn, zcrit):
mask_small = (abs(dxn) <= zcrit) mask_small = (np.abs(dxn) <= zcrit)
i_small = np.flatnonzero(mask_small) i_small = np.flatnonzero(mask_small)
if verbose: if verbose:
if zcrit == 0.: if zcrit == 0.:
@ -2130,7 +2122,7 @@ def findoutliers(x, zcrit=0.0, dcrit=None, ddcrit=None, verbose=False):
indg, = nonzero(indg) indg, = nonzero(indg)
if verbose: 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 return ind, indg
@ -2174,8 +2166,7 @@ def common_shape(*args, ** kwds):
''' '''
shape = kwds.get('shape') shape = kwds.get('shape')
x0 = 1 if shape is None else np.ones(shape) x0 = 1 if shape is None else np.ones(shape)
x1 = np.broadcast(x0, *args) return tuple(np.broadcast(x0, *args).shape)
return tuple(x1.shape)
def argsreduce(condition, * args): def argsreduce(condition, * args):
@ -2275,15 +2266,15 @@ def stirlerr(n):
n500 = 500 < n1 n500 = 500 < n1
y[n500] = (S0 - S1 / nn[n500]) / n1[n500] y[n500] = (S0 - S1 / nn[n500]) / n1[n500]
n80 = logical_and(80 < n1, n1 <= 500) 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] y[n80] = (S0 - (S1 - S2 / nn[n80]) / nn[n80]) / n1[n80]
n35 = logical_and(35 < n1, n1 <= 80) n35 = logical_and(35 < n1, n1 <= 80)
if any(n35): if np.any(n35):
nn35 = nn[n35] nn35 = nn[n35]
y[n35] = (S0 - (S1 - (S2 - S3 / nn35) / nn35) / nn35) / n1[n35] y[n35] = (S0 - (S1 - (S2 - S3 / nn35) / nn35) / nn35) / n1[n35]
n15 = logical_and(15 < n1, n1 <= 35) n15 = logical_and(15 < n1, n1 <= 35)
if any(n15): if np.any(n15):
nn15 = nn[n15] nn15 = nn[n15]
y[n15] = ( y[n15] = (
S0 - (S1 - (S2 - (S3 - S4 / nn15) / nn15) / nn15) / nn15) / n1[n15] S0 - (S1 - (S2 - (S3 - S4 / nn15) / nn15) / nn15) / nn15) / n1[n15]
@ -2426,7 +2417,7 @@ def binomln(z, w):
Example Example
------- -------
>>> abs(binomln(3,2)- 1.09861229)<1e-7 >>> np.abs(binomln(3,2)- 1.09861229)<1e-7
array([ True], dtype=bool) array([ True], dtype=bool)
See also See also
@ -2463,7 +2454,7 @@ def betaloge(z, w):
Example Example
------- -------
>>> import wafo.misc as wm >>> 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) array([ True], dtype=bool)
See also See also
@ -2545,7 +2536,7 @@ def nextpow2(x):
if (t > 1): if (t > 1):
f, n = frexp(t) f, n = frexp(t)
else: else:
f, n = frexp(abs(x)) f, n = frexp(np.abs(x))
if (f == 0.5): if (f == 0.5):
n = n - 1 n = n - 1
@ -2615,7 +2606,7 @@ def _discretize_linear(fun, a, b, tol=0.005, n=5):
x = linspace(a, b, n) x = linspace(a, b, n)
y = fun(x) y = fun(x)
y00 = interp(x, x0, y0) 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 return x, y
@ -2644,7 +2635,7 @@ def _discretize_adaptive(fun, a, b, tol=0.005, n=5):
fy = fun(y) fy = fun(y)
fy0 = interp(y, x, fx) 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() 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] xn = xo[-1]
x0 = xo[0] x0 = xo[0]
L = float(xn - x0) 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 # pab 07.01.2001: Always choose the stepsize df so that
# it is an exactly representable number. # it is an exactly representable number.
# This is important when calculating numerical derivatives and is # 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 True
''' '''
def _denominator_text(den): 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): def _numerator_text(num):
if abs(num) == 1: if np.abs(num) == 1:
return '-' if num == -1 else '' return '-' if num == -1 else ''
return '{:d}'.format(num) return '{:d}'.format(num)
frac = fractions.Fraction.from_float(x / pi).limit_denominator(int(1e+13)) frac = fractions.Fraction.from_float(x / pi).limit_denominator(int(1e+13))

56
wafo/sg_filter.py
Unescape Escape View File

@ -2,7 +2,7 @@ from __future__ import absolute_import, division
import numpy as np import numpy as np
# from math import pow # from math import pow
# from numpy import zeros,dot # 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 import spdiags
from scipy.sparse.linalg import spsolve, expm from scipy.sparse.linalg import spsolve, expm
from scipy.signal import medfilt from scipy.signal import medfilt
@ -132,14 +132,14 @@ class SavitzkyGolay(object):
def smooth_last(self, signal, k=0): def smooth_last(self, signal, k=0):
coeff = self._coeff coeff = self._coeff
n = size(coeff - 1) // 2 n = np.size((coeff - 1) // 2
y = np.squeeze(signal) y = np.squeeze(signal)
if n == 0: if n == 0:
return y return y
if y.ndim > 1: if y.ndim > 1:
coeff.shape = (-1, 1) coeff.shape = (-1, 1)
first_vals = y[0] - abs(y[n:0:-1] - y[0]) first_vals = y[0] - np.abs(y[n:0:-1] - y[0])
last_vals = y[-1] + abs(y[-2:-n - 2:-1] - y[-1]) last_vals = y[-1] + np.abs(y[-2:-n - 2:-1] - y[-1])
y = concatenate((first_vals, y, last_vals)) y = concatenate((first_vals, y, last_vals))
return (y[-2 * n - 1 - k:-k] * coeff).sum(axis=0) return (y[-2 * n - 1 - k:-k] * coeff).sum(axis=0)
@ -147,9 +147,8 @@ class SavitzkyGolay(object):
return self.smooth(signal) return self.smooth(signal)
def smooth(self, signal): def smooth(self, signal):
x = np.asarray(signal) dtype = np.result_type(signal, np.float)
if x.dtype != np.float64 and x.dtype != np.float32: x = np.asarray(signal, dtype=dtype)
x = x.astype(np.float64)
coeffs = self._coeff coeffs = self._coeff
mode, axis = self.mode, self.axis mode, axis = self.mode, self.axis
@ -180,13 +179,13 @@ class SavitzkyGolay(object):
the smoothed signal (or it's n-th derivative). the smoothed signal (or it's n-th derivative).
""" """
coeff = self._coeff coeff = self._coeff
n = size(coeff - 1) // 2 n = np.size(coeff - 1) // 2
y = np.squeeze(signal) y = np.squeeze(signal)
if n == 0: if n == 0:
return y return y
if pad: if pad:
first_vals = y[0] - abs(y[n:0:-1] - y[0]) first_vals = y[0] - np.abs(y[n:0:-1] - y[0])
last_vals = y[-1] + abs(y[-2:-n - 2:-1] - y[-1]) last_vals = y[-1] + np.abs(y[-2:-n - 2:-1] - y[-1])
y = concatenate((first_vals, y, last_vals)) y = concatenate((first_vals, y, last_vals))
n *= 2 n *= 2
d = y.ndim d = y.ndim
@ -424,8 +423,8 @@ class _Filter(object):
@staticmethod @staticmethod
def _studentized_residuals(r, I, h): def _studentized_residuals(r, I, h):
median_abs_deviation = np.median(abs(r[I] - np.median(r[I]))) median_abs_deviation = np.median(np.abs(r[I] - np.median(r[I])))
return abs(r / (1.4826 * median_abs_deviation) / sqrt(1 - h)) return np.abs(r / (1.4826 * median_abs_deviation) / sqrt(1 - h))
def robust_weights(self, r, I, h): def robust_weights(self, r, I, h):
"""Return weights for robust smoothing.""" """Return weights for robust smoothing."""
@ -458,10 +457,10 @@ class _Filter(object):
def check_smooth_parameter(self, s): def check_smooth_parameter(self, s):
if self.auto_smooth: 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. warnings.warn('''s = %g: the lower bound for s has been reached.
Put s as an input variable if required.''' % s) 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. warnings.warn('''s = %g: the Upper bound for s has been reached.
Put s as an input variable if required.''' % s) 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]) y[np.r_[70, 75, 80]] = np.array([5.5, 5, 6])
z = smoothn(y) # Regular smoothing z = smoothn(y) # Regular smoothing
zr = smoothn(y, robust=True) # Robust smoothing zr = smoothn(y, robust=True) # Robust smoothing
plt.subplot(121), _h0 = plt.subplot(121),
unused_h = plt.plot(x, y, 'r.', x, z, 'k', linewidth=2) _h = plt.plot(x, y, 'r.', x, z, 'k', linewidth=2)
plt.title('Regular smoothing') plt.title('Regular smoothing')
plt.subplot(122) plt.subplot(122)
plt.plot(x, y, 'r.', x, zr, 'k', linewidth=2) plt.plot(x, y, 'r.', x, zr, 'k', linewidth=2)
@ -735,18 +734,18 @@ def test_smoothn_1d():
def test_smoothn_2d(): def test_smoothn_2d():
# import mayavi.mlab as plt # import mayavi.mlab as plt
xp = np.r_[0:1:.02] xp = np.r_[0:1:0.02]
[x, y] = np.meshgrid(xp, xp) [x, y] = np.meshgrid(xp, xp)
f = np.exp(x + y) + np.sin((x - 2 * y) * 3) f = np.exp(x + y) + np.sin((x - 2 * y) * 3)
fn = f + np.random.randn(*f.shape) * 0.5 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) fs2 = smoothn(fn, s=2 * s)
plt.subplot(131), _h = plt.subplot(131),
plt.contourf(xp, xp, fn) _h = plt.contourf(xp, xp, fn)
plt.subplot(132), _h = plt.subplot(132),
plt.contourf(xp, xp, fs2) _h = plt.contourf(xp, xp, fs2)
plt.subplot(133), _h = plt.subplot(133),
plt.contourf(xp, xp, f) _h = plt.contourf(xp, xp, f)
plt.show('hold') plt.show('hold')
@ -1048,7 +1047,7 @@ class Kalman(object):
def _update_covariance(self, P, K): def _update_covariance(self, P, K):
return P - np.dot(K, np.dot(self.H, P)) 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): def _filter_main(self, z, u):
''' This is the code which implements the discrete Kalman filter: ''' This is the code which implements the discrete Kalman filter:
@ -1187,7 +1186,7 @@ def test_kalman_sine():
x = np.zeros((n, m)) x = np.zeros((n, m))
for i, zi in enumerate(z): 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') _hz = plt.plot(z, 'r.', label='observations')
# a-posteriori state estimates: # a-posteriori state estimates:
@ -1321,7 +1320,8 @@ class HampelFilter(object):
self.adaptive = adaptive self.adaptive = adaptive
self.fulloutput = fulloutput self.fulloutput = fulloutput
def _check(self, dx): @staticmethod
def _check(dx):
if not np.isscalar(dx): if not np.isscalar(dx):
raise ValueError('DX must be a scalar.') raise ValueError('DX must be a scalar.')
if dx < 0: if dx < 0:
@ -1357,7 +1357,7 @@ class HampelFilter(object):
while S0Rel > self.adaptive: while S0Rel > self.adaptive:
Y0Tmp[j], S0Tmp[j] = localwindow(X, Y, DXTmp[j], i) Y0Tmp[j], S0Tmp[j] = localwindow(X, Y, DXTmp[j], i)
if j > 0: if j > 0:
S0Rel = abs((S0Tmp[j - 1] - S0Tmp[j]) / S0Rel = np.abs((S0Tmp[j - 1] - S0Tmp[j]) /
(S0Tmp[j - 1] + S0Tmp[j]) / 2) (S0Tmp[j - 1] + S0Tmp[j]) / 2)
j += 1 j += 1

Write Preview
Loading…
Cancel
Save