From 217a67893a99e88372239d8150334c466c396c9d Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Wed, 22 Aug 2012 00:03:27 +0000 Subject: [PATCH] Reordered methods --- pywafo/src/wafo/interpolate.py | 213 ++++++++++++++++----------------- 1 file changed, 106 insertions(+), 107 deletions(-) diff --git a/pywafo/src/wafo/interpolate.py b/pywafo/src/wafo/interpolate.py index f0f56d7..2d4c977 100644 --- a/pywafo/src/wafo/interpolate.py +++ b/pywafo/src/wafo/interpolate.py @@ -692,13 +692,99 @@ def slopes(x,y, method='parabola', tension=0, monotone=False): yp[ii+1] = tau[ii] * beta[ii] * dydx[ii] return yp - -class StinemanInterp2(PiecewisePolynomial): - def __init__(self, x, y, yp=None, method='parabola', monotone=False): - if yp is None: - yp = slopes(x, y, method, monotone=monotone) - super(StinemanInterp,self).__init__(x, zip(y,yp)) +def stineman_interp(xi, x, y, yp=None): + """ + Given data vectors *x* and *y*, the slope vector *yp* and a new + abscissa vector *xi*, the function :func:`stineman_interp` uses + Stineman interpolation to calculate a vector *yi* corresponding to + *xi*. + + Here's an example that generates a coarse sine curve, then + interpolates over a finer abscissa:: + + x = linspace(0,2*pi,20); y = sin(x); yp = cos(x) + xi = linspace(0,2*pi,40); + yi = stineman_interp(xi,x,y,yp); + plot(x,y,'o',xi,yi) + + The interpolation method is described in the article A + CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell + W. Stineman. The article appeared in the July 1980 issue of + Creative Computing with a note from the editor stating that while + they were: + + not an academic journal but once in a while something serious + and original comes in adding that this was + "apparently a real solution" to a well known problem. + + For *yp* = *None*, the routine automatically determines the slopes + using the :func:`slopes` routine. + + *x* is assumed to be sorted in increasing order. + + For values ``xi[j] < x[0]`` or ``xi[j] > x[-1]``, the routine + tries an extrapolation. The relevance of the data obtained from + this, of course, is questionable... + + Original implementation by Halldor Bjornsson, Icelandic + Meteorolocial Office, March 2006 halldor at vedur.is + + Completely reworked and optimized for Python by Norbert Nemec, + Institute of Theoretical Physics, University or Regensburg, April + 2006 Norbert.Nemec at physik.uni-regensburg.de + """ + + # Cast key variables as float. + x = np.asarray(x, np.float_) + y = np.asarray(y, np.float_) + assert x.shape == y.shape + #N = len(y) + + if yp is None: + yp = slopes(x, y) + else: + yp = np.asarray(yp, np.float_) + + xi = np.asarray(xi, np.float_) + #yi = np.zeros(xi.shape, np.float_) + + # calculate linear slopes + dx = x[1:] - x[:-1] + dy = y[1:] - y[:-1] + s = dy / dx #note length of s is N-1 so last element is #N-2 + + # find the segment each xi is in + # this line actually is the key to the efficiency of this implementation + idx = np.searchsorted(x[1:-1], xi) + + # now we have generally: x[idx[j]] <= xi[j] <= x[idx[j]+1] + # except at the boundaries, where it may be that xi[j] < x[0] or xi[j] > x[-1] + + # the y-values that would come out from a linear interpolation: + sidx = s.take(idx) + xidx = x.take(idx) + yidx = y.take(idx) + xidxp1 = x.take(idx + 1) + yo = yidx + sidx * (xi - xidx) + + # the difference that comes when using the slopes given in yp + dy1 = (yp.take(idx) - sidx) * (xi - xidx) # using the yp slope of the left point + dy2 = (yp.take(idx + 1) - sidx) * (xi - xidxp1) # using the yp slope of the right point + + dy1dy2 = dy1 * dy2 + # The following is optimized for Python. The solution actually + # does more calculations than necessary but exploiting the power + # of numpy, this is far more efficient than coding a loop by hand + # in Python + dy1mdy2 = np.where(dy1dy2,dy1-dy2,np.inf) + dy1pdy2 = np.where(dy1dy2,dy1+dy2,np.inf) + yi = yo + dy1dy2 * np.choose(np.array(np.sign(dy1dy2), np.int32) + 1, + ((2 * xi - xidx - xidxp1) / ((dy1mdy2) * (xidxp1 - xidx)), + 0.0, + 1 / (dy1pdy2))) + return yi + class StinemanInterp(object): ''' Returns the values of an interpolating function that runs through a set of points according to the algorithm of Stineman (1980). @@ -811,100 +897,13 @@ class StinemanInterp(object): 0.0, 1 / (dy1pdy2))) return yi - + +class StinemanInterp2(PiecewisePolynomial): + def __init__(self, x, y, yp=None, method='parabola', monotone=False): + if yp is None: + yp = slopes(x, y, method, monotone=monotone) + super(StinemanInterp2,self).__init__(x, zip(y,yp)) -def stineman_interp(xi, x, y, yp=None): - """ - Given data vectors *x* and *y*, the slope vector *yp* and a new - abscissa vector *xi*, the function :func:`stineman_interp` uses - Stineman interpolation to calculate a vector *yi* corresponding to - *xi*. - - Here's an example that generates a coarse sine curve, then - interpolates over a finer abscissa:: - - x = linspace(0,2*pi,20); y = sin(x); yp = cos(x) - xi = linspace(0,2*pi,40); - yi = stineman_interp(xi,x,y,yp); - plot(x,y,'o',xi,yi) - - The interpolation method is described in the article A - CONSISTENTLY WELL BEHAVED METHOD OF INTERPOLATION by Russell - W. Stineman. The article appeared in the July 1980 issue of - Creative Computing with a note from the editor stating that while - they were: - - not an academic journal but once in a while something serious - and original comes in adding that this was - "apparently a real solution" to a well known problem. - - For *yp* = *None*, the routine automatically determines the slopes - using the :func:`slopes` routine. - - *x* is assumed to be sorted in increasing order. - - For values ``xi[j] < x[0]`` or ``xi[j] > x[-1]``, the routine - tries an extrapolation. The relevance of the data obtained from - this, of course, is questionable... - - Original implementation by Halldor Bjornsson, Icelandic - Meteorolocial Office, March 2006 halldor at vedur.is - - Completely reworked and optimized for Python by Norbert Nemec, - Institute of Theoretical Physics, University or Regensburg, April - 2006 Norbert.Nemec at physik.uni-regensburg.de - """ - - # Cast key variables as float. - x = np.asarray(x, np.float_) - y = np.asarray(y, np.float_) - assert x.shape == y.shape - #N = len(y) - - if yp is None: - yp = slopes(x, y) - else: - yp = np.asarray(yp, np.float_) - - xi = np.asarray(xi, np.float_) - #yi = np.zeros(xi.shape, np.float_) - - # calculate linear slopes - dx = x[1:] - x[:-1] - dy = y[1:] - y[:-1] - s = dy / dx #note length of s is N-1 so last element is #N-2 - - # find the segment each xi is in - # this line actually is the key to the efficiency of this implementation - idx = np.searchsorted(x[1:-1], xi) - - # now we have generally: x[idx[j]] <= xi[j] <= x[idx[j]+1] - # except at the boundaries, where it may be that xi[j] < x[0] or xi[j] > x[-1] - - # the y-values that would come out from a linear interpolation: - sidx = s.take(idx) - xidx = x.take(idx) - yidx = y.take(idx) - xidxp1 = x.take(idx + 1) - yo = yidx + sidx * (xi - xidx) - - # the difference that comes when using the slopes given in yp - dy1 = (yp.take(idx) - sidx) * (xi - xidx) # using the yp slope of the left point - dy2 = (yp.take(idx + 1) - sidx) * (xi - xidxp1) # using the yp slope of the right point - - dy1dy2 = dy1 * dy2 - # The following is optimized for Python. The solution actually - # does more calculations than necessary but exploiting the power - # of numpy, this is far more efficient than coding a loop by hand - # in Python - dy1mdy2 = np.where(dy1dy2,dy1-dy2,np.inf) - dy1pdy2 = np.where(dy1dy2,dy1+dy2,np.inf) - yi = yo + dy1dy2 * np.choose(np.array(np.sign(dy1dy2), np.int32) + 1, - ((2 * xi - xidx - xidxp1) / ((dy1mdy2) * (xidxp1 - xidx)), - 0.0, - 1 / (dy1pdy2))) - return yi - class CubicHermiteSpline(PiecewisePolynomial): ''' Piecewise Cubic Hermite Interpolation using Catmull-Rom @@ -1008,10 +1007,10 @@ def compare_methods(): ############################################################ # Sine wave test ############################################################ - + fun = np.sin # Create a example vector containing a sine wave. x = np.arange(30.0)/10. - y = np.sin(x) + y = fun(x) # Interpolate the data above to the grid defined by "xvec" xvec = np.arange(250.)/100. @@ -1032,12 +1031,12 @@ def compare_methods(): import matplotlib.pyplot as plt plt.figure() - plt.plot(x,y, 'r.') + plt.plot(x,y, 'ro', xvec, fun(xvec),'r') plt.title("pchip() Sin test code") # Plot the interpolated points - plt.plot(xvec, yvec, xvec, yvec1, xvec, yvec2,'go',xvec, yvec3) - plt.legend(['true','parbola_monoton','parabola','catmul','pchip'], frameon=False, loc=0) + plt.plot(xvec, yvec, xvec, yvec1, xvec, yvec2,'g.',xvec, yvec3) + plt.legend(['true','true','parbola_monoton','parabola','catmul','pchip'], frameon=False, loc=0) plt.ioff() plt.show() @@ -1150,8 +1149,8 @@ def test_docstrings(): if __name__ == '__main__': - test_docstrings() +# test_docstrings() #main() #test_smoothing_spline() -# compare_methods() - # demo_monoticity() \ No newline at end of file + #compare_methods() + demo_monoticity() \ No newline at end of file