[Scipysvn] r4560  branches/Interpolate1D
scipysvn@scip...
scipysvn@scip...
Thu Jul 24 14:43:55 CDT 2008
Author: fcady
Date: 20080724 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 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/Interpolate1D.py 20080724 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 inrange, belowrange, 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 : order2 splines work on linear data
make sure : order2 splines work on nonlinear 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 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/TODO.txt 20080724 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 nonzero 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
+ Underthehood 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, commentedout 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 2dimensional
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.
+**highlevel road map
+ when the module is more established, there should be a page on
+ the wiki which describes the bigpicture of the module; what
+ the capabilities are and which should be added, largescale
+ architecture of the module, etc.
**pick best spline
Underthehood 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.
**highlevel road map
when the module is more established, there should be a page on
the wiki which describes the bigpicture of the module; what
the capabilities are and which should be added, largescale
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 2dimensional
+ 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 underthehood interpolation is performed on contiguous
+ arrays.
+
+ Whatever is done must be very well commented.
Modified: branches/Interpolate1D/__init__.py
===================================================================
 branches/Interpolate1D/__init__.py 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/__init__.py 20080724 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 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/fitpack_wrapper.py 20080724 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 # lowerlevel 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  1d 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(fps)/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 leastsquares
 # 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 nuth derivative) at positions x.
Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
 branches/Interpolate1D/interpolate1d.py 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/interpolate1d.py 20080724 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 inrange, belowrange, 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 : order2 splines work on linear data
make sure : order2 splines work on nonlinear 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 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/interpolate_wrapper.py 20080724 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 20080723 20:49:37 UTC (rev 4559)
+++ branches/Interpolate1D/setup.py 20080724 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 Scipysvn
mailing list