[Scipy-svn] r4560 - branches/Interpolate1D
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Jul 24 15:43:55 EDT 2008
Author: fcady
Date: 2008-07-24 14:43:52 -0500 (Thu, 24 Jul 2008)
New Revision: 4560
Modified:
branches/Interpolate1D/Interpolate1D.py
branches/Interpolate1D/TODO.txt
branches/Interpolate1D/__init__.py
branches/Interpolate1D/fitpack_wrapper.py
branches/Interpolate1D/interpolate1d.py
branches/Interpolate1D/interpolate_wrapper.py
branches/Interpolate1D/setup.py
Log:
better init, fitpack_wrapper trimmer down, and TODO.txt more clear and informative
Modified: branches/Interpolate1D/Interpolate1D.py
===================================================================
--- branches/Interpolate1D/Interpolate1D.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/Interpolate1D.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -11,7 +11,7 @@
Classes provided include:
- interpolate1d : an object for interpolation of
+ Interpolate1d : an object for interpolation of
various kinds. interp1d is a wrapper
around this class.
@@ -117,7 +117,7 @@
>>> interp1d(x, y, new_x)
array([.2, 2.3, 5.6, NaN])
"""
- return Interpolate1D(x, y, kind=kind, low=low, high=high, \
+ return Interpolate1d(x, y, kind=kind, low=low, high=high, \
kindkw=kindkw, lowkw=lowkw, highkw=highkw, \
remove_bad_data = remove_bad_data, bad_data=bad_data)(new_x)
@@ -255,7 +255,6 @@
from inspect import isclass, isfunction
- # FIXME : more string options available ('cubic', etc)
if interp_arg in ['linear', 'logarithmic', 'block', 'block_average_above']:
# string used to indicate interpolation method, Select appropriate function
func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
@@ -264,6 +263,13 @@
elif interp_arg in ['Spline', Spline, 'spline']:
# spline is a special case of above
result = Spline(self._x, self._y, **kw)
+ elif interp_arg in ['cubic', 'Cubic', 'Quadratic', \
+ 'quadratic', 'Quad', 'quad']:
+ # specify certain kinds of splines
+ if interp_arg in ['cubic', 'Cubic']:
+ result = Spline(self._x, self._y, k=3)
+ elif interp_arg in ['Quadratic', 'quadratic', 'Quad', 'quad']:
+ result = Spline(self._x, self._y, k=2)
elif isfunction(interp_arg):
# assume user has passed a function
result = lambda new_x : interp_arg(new_x, **kw)
@@ -279,7 +285,8 @@
Breaks x into pieces in-range, below-range, and above range.
Performs appropriate operation on each and concatenates results.
"""
-
+ # FIXME : make_array_safe may also be called within the interpolation technique.
+ # waste of time, but ok for the time being.
x = make_array_safe(x)
# masks indicate which elements fall into which interpolation region
@@ -349,6 +356,8 @@
"""
make sure : order-2 splines work on linear data
make sure : order-2 splines work on non-linear data
+ make sure : 'cubic' and 'quad' as arguments yield
+ the desired spline
"""
print "\n\nTESTING 2nd ORDER SPLINE"
N = 7 #must be > 5
@@ -369,7 +378,7 @@
N = 7
x = np.arange(N)
y = x**2
- interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':2}, low='spline', high='spline')
+ interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':2}, low='quad', high='cubic')
new_x = np.arange(N+1)-0.5
new_y = interp_func(new_x)
self.assertAllclose(new_x**2, new_y)
Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/TODO.txt 2008-07-24 19:43:52 UTC (rev 4560)
@@ -5,81 +5,121 @@
at appropriate places in the code.
+**handle smoothing in fitpack_wrapper
+ The default smoothing parameter (s) means we aren't actually
+ doing interpolation. Also, small but non-zero s makes the
+ algorithm very slow. Figure out how to handle this, prob by
+ not letting the user see s and setting it to 0.
+
+ This may be best handled by modifying fitpack_wrapper.py.
+
+
+**pick best spline
+ Under-the-hood machinery currently comes from _interpolate.cpp
+ (used in enthought.interpolate) and FITPACK (Fortran, used in
+ scipy.interpolate). This isn't necessarily the best (for example,
+ speed of FITPACK is highly sensitive to parameter s). Other code
+ is used in scipy.ndimage and scipy.signal. There is surely other
+ code out there too. Figure out what is best and incorporate it.
+
+ Signal code is slower than FITPACK, and NDImage really only does
+ filtering. I'm inclined to stay with FITPACK.
+
+
+**clean up fitpack_wrapper.py
+ Currently is all hacked together from scipy.interpolate.
+ There is unused functionality, commented-out functionality,
+ and other undesirables. Once it has been established what
+ we will include, that stuff needs to be cleaned up, and the
+ rest needs to be either removed or set aside.
+
+
**comment interpolate1d
-There's comments there already, but they should be
-made better.
+ There's comments there already, but they should be
+ made better.
**doc strings for interpolate1d and its members
-There's docstrings there already, but they should be
-made better. In particular, it must be ensured that
-they are of the proper format and include examples.
+ There's docstrings there already, but they should be
+ made better. In particular, it must be ensured that
+ they are of the proper format and include examples.
-The doc strings for __init__.py, interpolate1d.py,
-Interpolate1d, and interp1d are virtually identical
-and very long; perhaps a master string can be stored
-somewhere that they all reference. This would make
-updates of documentation easier.
+ The doc strings for __init__.py, interpolate1d.py,
+ Interpolate1d, and interp1d are virtually identical
+ and very long; perhaps a master string can be stored
+ somewhere that they all reference. This would make
+ updates of documentation easier.
-**more strings user can pass ('cubic', etc)
-User can specify interpolation type as a string argument
-to interpolate1d at initialization. More strings should work.
-
-
**figure out NumPy version stuff with vectorize.
-In function interpolate1d._format_array.
-It would be nice to remove the hack I used.
-I believe vectorize is supposed to handle arrays of
-length 0, but it's not working on my computer.
+ In function interpolate1d._format_array.
+ It would be nice to remove the hack I used.
+ I believe vectorize is supposed to handle arrays of
+ length 0, but it's not working on my computer.
**better handling of variable types
-Currently everything is cast to a float64 if it is not already
-a float32. Is this the best way to do it?
+ Currently everything is cast to a float64 if it is not already
+ a float32. Is this the best way to do it?
-Also, for the future, code should be added for record arrays,
-which mix real values with strings. This is, I believe already
-largely supported, but that's not because the code was written
-with that in mind. I haven't thought through the details.
+ Also, for the future, code should be added for record arrays,
+ which mix real values with strings. This is, I believe already
+ largely supported, but that's not because the code was written
+ with that in mind. I haven't thought through the details.
-Perhaps this should be done as another function/class which
-wraps interpolate1d.
+ Perhaps this should be done as another function/class which
+ wraps interpolate1d.
-**allow y to be 2-dimensional
-That way the interpolated function is from R1 -> Rn, and
-not just R1 -> R1. This requires some thinking about axes.
+**improve regression tests
+ desired for fitpack_wrapper and _interpolate_wrapper
+ as well as interpolate1d.
-**improve regression tests
-desired for fitpack_wrapper and _interpolate_wrapper
-as well as interpolate1d.
+**high-level road map
+ when the module is more established, there should be a page on
+ the wiki which describes the big-picture of the module; what
+ the capabilities are and which should be added, large-scale
+ architecture of the module, etc.
-**pick best spline
-Under-the-hood machinery currently comes from _interpolate.cpp
-(used in enthought.interpolate) and FITPACK (Fortran, used in
-scipy.interpolate). This isn't necessarily the best (for example,
-speed of FITPACK is highly sensitive to parameter k). Other code
-is used in scipy.ndimage and scipy.signal. There is surely other
-code out there too. Figure out what is best and incorporate it.
+ It might note which underlying C/Fortran modules can or should
+ be modified or merged. It would be great if either 1) there were
+ only 1 extension module, or 2) the modules showed natural
+ differentiation of functionality (one for splines, one for simple
+ operations, etc), rather than being a holdover of where they
+ were stolen from.
-**high-level road map
-when the module is more established, there should be a page on
-the wiki which describes the big-picture of the module; what
-the capabilities are and which should be added, large-scale
-architecture of the module, etc.
+**update for 2D and ND
+ This will probably take the form of two additional
+ classes both based on interpolate1d. Thus it probably
+ shouldn't be done until interpolate1d is more settled.
-It might note which underlying C/Fortran modules can or should
-be modified or merged. It would be great if either 1) there were
-only 1 extension module, or 2) the modules showed natural
-differentiation of functionality (one for splines, one for simple
-operations, etc), rather than being a holdover of where they
-were stolen from.
+**more convenient way to enter kw arguments
+ currently users have to pass in dictionaries of additional
+ keyword arguments for specific interpolation types.
+ This is kind of awkward, and it would be nice to let
+ them pass in a single argument. But this is low priority.
-**update for 2D and ND
-This will probably take the form of two additional
-classes both based on interpolate1d. Thus it probably
-shouldn't be done until interpolate1d is more settled.
+
+**allow y to be 2-dimensional
+ That way the interpolated function is from R1 -> Rn, and
+ not just R1 -> R1. This requires some thinking about axes.
+
+ The interpolation should, I think, be along columns, so each
+ row of y is one complete observation. This is because 1) rows
+ are separated more from each other when an array is displayed,
+ and 2) the user may want to enter the data as y=array([obs1, obs2, ...]).
+
+ There are two big problem. First, interpolate_wrapper interpolates
+ along rows, because each row is contiguous. Second,
+ FITPACK doesn't support 2D y arrays.
+
+ The solution is to have the user indicate which axis to interpolate
+ along, and take the transpose if need be (also copy, so that
+ the columns become contiguous). Then when the user enters newx,
+ all under-the-hood interpolation is performed on contiguous
+ arrays.
+
+ Whatever is done must be very well commented.
Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/__init__.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -1,32 +1,12 @@
-#FIXME : better docstring
-"""
- Interpolation of 1D data
- This module provides several functions and classes for interpolation
- and extrapolation of 1D data (1D in both input and output). The
- primary function provided is:
+# information about included functions
+from info import __doc__
- interp1d(x, y, new_x) : from data points x and y, interpolates
- values for points in new_x and
- returns them as an array.
-
- Classes provided include:
-
- Interpolate1d : an object for interpolation of
- various kinds. interp1d is a wrapper
- around this class.
-
- Spline : an object for spline interpolation
-
- Functions provided include:
-
- linear : linear interpolation
- logarithmic : logarithmic interpolation
- block : block interpolation
- block_average_above : block average above interpolation
-
-"""
-
+# basic interpolation routines
from interpolate_wrapper import linear, logarithmic, block, block_average_above
+
+#support for spline interpolation
from fitpack_wrapper import Spline
+
+# wrapping for all supported interpolation types.
from interpolate1d import Interpolate1d, interp1d
\ No newline at end of file
Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/fitpack_wrapper.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -2,19 +2,19 @@
This module is used for spline interpolation, and functions
as a wrapper around the FITPACK Fortran interpolation
package. It is not intended to be directly accessed by
-the user, but rather through the class Interpolate1D.
+the user, but rather through the class Interpolate1d.
The code has been modified from an older version of
scipy.interpolate, where it was directly called by the
user. As such, it includes functionality not available through
-Interpolate1D. For this reason, users may wish to get
+Interpolate1d. For this reason, users may wish to get
under the hood.
"""
import numpy as np
-import dfitpack
+import dfitpack # lower-level wrapper around FITPACK
class Spline(object):
@@ -30,7 +30,7 @@
BivariateSpline - a similar class for bivariate spline interpolation
"""
- def __init__(self, x, y, w=None, bbox = [None]*2, k=3, s=None):
+ def __init__(self, x, y, w=None, bbox = [None]*2, k=3, s=0.0):
"""
Input:
x,y - 1-d sequences of data points (x must be
@@ -44,46 +44,21 @@
k=3 - degree of the univariate spline.
s - positive smoothing factor defined for
estimation condition:
- sum((w[i]*(y[i]-s(x[i])))**2,axis=0) <= s
+ sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
Default s=len(w) which should be a good value
if 1/w[i] is an estimate of the standard
deviation of y[i].
"""
#_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
- data = dfitpack.fpcurf0(x,y,k,w=w,
- xb=bbox[0],xe=bbox[1],s=s)
+ data = dfitpack.fpcurf0(x, y, k, w=w,
+ xb=bbox[0], xe=bbox[1], s=s)
if data[-1]==1:
# nest too small, setting to maximum bound
data = self._reset_nest(data)
self._data = data
- self._reset_class()
-
- def _reset_class(self):
- data = self._data
+ # the relevant part of self._reset_class()
n,t,c,k,ier = data[7],data[8],data[9],data[5],data[-1]
self._eval_args = t[:n],c[:n],k
- if ier==0:
- # the spline returned has a residual sum of squares fp
- # such that abs(fp-s)/s <= tol with tol a relative
- # tolerance set to 0.001 by the program
- pass
- elif ier==-1:
- # the spline returned is an interpolating spline
- #self.__class__ = InterpolatedUnivariateSpline
- pass
- elif ier==-2:
- # the spline returned is the weighted least-squares
- # polynomial of degree k. In this extreme case fp gives
- # the upper bound fp0 for the smoothing factor s.
- #self.__class__ = LSQUnivariateSpline
- pass
- else:
- # error
- #if ier==1:
- # self.__class__ = LSQUnivariateSpline
- #message = _curfit_messages.get(ier,'ier=%s' % (ier))
- #warnings.warn(message)
- pass
def _reset_nest(self, data, nest=None):
n = data[10]
@@ -118,7 +93,7 @@
# nest too small, setting to maximum bound
data = self._reset_nest(data)
self._data = data
- self._reset_class()
+ #self._reset_class()
def __call__(self, x, nu=None):
""" Evaluate spline (or its nu-th derivative) at positions x.
Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/interpolate1d.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -11,7 +11,7 @@
Classes provided include:
- interpolate1d : an object for interpolation of
+ Interpolate1d : an object for interpolation of
various kinds. interp1d is a wrapper
around this class.
@@ -117,7 +117,7 @@
>>> interp1d(x, y, new_x)
array([.2, 2.3, 5.6, NaN])
"""
- return Interpolate1D(x, y, kind=kind, low=low, high=high, \
+ return Interpolate1d(x, y, kind=kind, low=low, high=high, \
kindkw=kindkw, lowkw=lowkw, highkw=highkw, \
remove_bad_data = remove_bad_data, bad_data=bad_data)(new_x)
@@ -255,7 +255,6 @@
from inspect import isclass, isfunction
- # FIXME : more string options available ('cubic', etc)
if interp_arg in ['linear', 'logarithmic', 'block', 'block_average_above']:
# string used to indicate interpolation method, Select appropriate function
func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
@@ -264,6 +263,13 @@
elif interp_arg in ['Spline', Spline, 'spline']:
# spline is a special case of above
result = Spline(self._x, self._y, **kw)
+ elif interp_arg in ['cubic', 'Cubic', 'Quadratic', \
+ 'quadratic', 'Quad', 'quad']:
+ # specify certain kinds of splines
+ if interp_arg in ['cubic', 'Cubic']:
+ result = Spline(self._x, self._y, k=3)
+ elif interp_arg in ['Quadratic', 'quadratic', 'Quad', 'quad']:
+ result = Spline(self._x, self._y, k=2)
elif isfunction(interp_arg):
# assume user has passed a function
result = lambda new_x : interp_arg(new_x, **kw)
@@ -279,7 +285,8 @@
Breaks x into pieces in-range, below-range, and above range.
Performs appropriate operation on each and concatenates results.
"""
-
+ # FIXME : make_array_safe may also be called within the interpolation technique.
+ # waste of time, but ok for the time being.
x = make_array_safe(x)
# masks indicate which elements fall into which interpolation region
@@ -349,6 +356,8 @@
"""
make sure : order-2 splines work on linear data
make sure : order-2 splines work on non-linear data
+ make sure : 'cubic' and 'quad' as arguments yield
+ the desired spline
"""
print "\n\nTESTING 2nd ORDER SPLINE"
N = 7 #must be > 5
@@ -369,7 +378,7 @@
N = 7
x = np.arange(N)
y = x**2
- interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':2}, low='spline', high='spline')
+ interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':2}, low='quad', high='cubic')
new_x = np.arange(N+1)-0.5
new_y = interp_func(new_x)
self.assertAllclose(new_x**2, new_y)
Modified: branches/Interpolate1D/interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/interpolate_wrapper.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/interpolate_wrapper.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -31,7 +31,7 @@
assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
if len(y.shape) == 2:
new_y = np.zeros((y.shape[0], len(new_x)), np.float64)
- for i in range(len(new_y)):
+ for i in range(len(new_y)): # for each row
_interpolate.linear_dddd(x, y[i], new_x, new_y[i])
else:
new_y = np.zeros(len(new_x), np.float64)
Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py 2008-07-23 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/setup.py 2008-07-24 19:43:52 UTC (rev 4560)
@@ -6,7 +6,7 @@
def configuration(parent_package='',top_path=None):
from numpy.distutils.misc_util import Configuration
- config = Configuration('', parent_package, top_path) #first arg was 'interpolate'
+ config = Configuration('', parent_package, top_path)
config.add_extension('_interpolate',
More information about the Scipy-svn
mailing list