[Scipysvn] r4564  branches/Interpolate1D
scipysvn@scip...
scipysvn@scip...
Mon Jul 28 14:41:45 CDT 2008
Author: fcady
Date: 20080728 14:41:44 0500 (Mon, 28 Jul 2008)
New Revision: 4564
Added:
branches/Interpolate1D/example_script.py
Modified:
branches/Interpolate1D/Interpolate1D.py
branches/Interpolate1D/TODO.txt
branches/Interpolate1D/interpolate1d.py
Log:
added example script demonstrating use of interpolate1d.py
Modified: branches/Interpolate1D/Interpolate1D.py
===================================================================
 branches/Interpolate1D/Interpolate1D.py 20080725 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/Interpolate1D.py 20080728 19:41:44 UTC (rev 4564)
@@ 130,7 +130,7 @@
x  list or NumPy array
x includes the xvalues for the data set to
interpolate from. It must be sorted in
 ascending order
+ ascending order.
y  list or NumPy array
y includes the yvalues for the data set to
@@ 159,7 +159,7 @@
a number') for all values outside the range of x.
remove_bad_data  bool
 indicates whether to remove bad data.
+ indicates whether to remove bad data points from x and y.
bad_data  list
List of values (in x or y) which indicate unacceptable data. All points
@@ 168,7 +168,7 @@
numpy.NaN is always considered bad data.
 Acceptable Input Strings
+ Some Acceptable Input Strings

"linear"  linear interpolation : default
@@ 192,7 +192,8 @@
"""
# FIXME: more informative descriptions of sample arguments
# FIXME: examples in doc string
 # FIXME : Allow copying or not of arrays. noncopy + remove_bad_data should flash a warning (esp if we interpolate missing values), but work anyway.
+ # FIXME : Allow copying or not of arrays. noncopy + remove_bad_data should flash
+ # a warning (esp if we interpolate missing values), but work anyway.
def __init__(self, x, y, kind='linear', low=np.NaN, high=np.NaN, \
kindkw={}, lowkw={}, highkw={}, \
@@ 203,8 +204,7 @@
# check acceptable size and dimensions
x = np.array(x)
y = np.array(y)
 assert len(x) > 0 and len(y) > 0 , "Interpolate1D does not support\
 arrays of length 0"
+ assert len(x) > 0 and len(y) > 0 , "Arrays cannot be of zero length"
assert x.ndim == 1 , "x must be onedimensional"
assert y.ndim == 1 , "y must be onedimensional"
assert len(x) == len(y) , "x and y must be of the same length"
@@ 213,13 +213,14 @@
if remove_bad_data:
x, y = self._remove_bad_data(x, y, bad_data)
+ # store data
# FIXME : may be good to let x and y be initialized later, or changed afterthefact
self._init_xy(x, y)
# store interpolation functions for each range
 self.kind = self._init_interp_method(self._x, self._y, kind, kindkw)
 self.low = self._init_interp_method(self._x, self._y, low, lowkw)
 self.high = self._init_interp_method(self._x, self._y, high, highkw)
+ self.kind = self._init_interp_method(kind, kindkw)
+ self.low = self._init_interp_method(low, lowkw)
+ self.high = self._init_interp_method(high, highkw)
def _remove_bad_data(self, x, y, bad_data = [None, np.NaN]):
""" removes data points whose x or y coordinate is
@@ 242,7 +243,7 @@
self._x = make_array_safe(x, self._xdtype).copy()
self._y = make_array_safe(y, self._ydtype).copy()
 def _init_interp_method(self, x, y, interp_arg, kw):
+ def _init_interp_method(self, interp_arg, kw):
"""
User provides interp_arg and dictionary kw. _init_interp_method
returns the interpolating function from x and y specified by interp_arg,
@@ 295,7 +296,11 @@
result.set_xy(self._x, self._y)
# user passes a function to be called
 elif isfunction(interp_arg) and interp.func_code.argcount == 3:
+ # FIXME : I think there is too much flexibility allowed here; it makes
+ # there be more pathological side cases to consider. Functions
+ # should perhaps be reqired to be of the form f(x, y, newx, **kw)
+ elif isfunction(interp_arg) and interp.func_code.argcount >= 3:
+ # assume x, y and newx are all passed to interp_arg
result = lambda new_x : interp_arg(self._x, self._y, new_x, **kw)
elif isfunction(interp_arg):
result = lambda new_x : interp_arg(new_x, **kw)
@@ 306,7 +311,7 @@
return result
 def __call__(self, x):
+ def __call__(self, newx):
"""
Input x must be in sorted order.
Breaks x into pieces inrange, belowrange, and above range.
@@ 314,24 +319,24 @@
"""
# 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)
+ newx = make_array_safe(newx)
# masks indicate which elements fall into which interpolation region
 low_mask = x<self._x[0]
 high_mask = x>self._x[1]
+ low_mask = newx<self._x[0]
+ high_mask = newx>self._x[1]
interp_mask = (~low_mask) & (~high_mask)
# use correct function for x values in each region
 if len(x[low_mask]) == 0: new_low=np.array([]) # FIXME : remove need for if/else.
+ if len(newx[low_mask]) == 0: new_low=np.array([]) # FIXME : remove need for if/else.
# if/else is a hack, since vectorize is failing
# to work on lists/arrays of length 0
# on the computer where this is being
# developed
 else: new_low = self.low(x[low_mask])
 if len(x[interp_mask])==0: new_interp=np.array([])
 else: new_interp = self.kind(x[interp_mask])
 if len(x[high_mask]) == 0: new_high = np.array([])
 else: new_high = self.high(x[high_mask])
+ else: new_low = self.low(newx[low_mask])
+ if len(newx[interp_mask])==0: new_interp=np.array([])
+ else: new_interp = self.kind(newx[interp_mask])
+ if len(newx[high_mask]) == 0: new_high = np.array([])
+ else: new_high = self.high(newx[high_mask])
result = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
# Would be nice to say result = zeros(dtype=?)
@@ 349,7 +354,7 @@
def test_interpolate_wrapper(self):
""" run unit test contained in interpolate_wrapper.py
"""
 print "\n\nTESTING _interpolate_wrapper MODULE"
+ #print "\n\nTESTING _interpolate_wrapper MODULE"
from interpolate_wrapper import Test
T = Test()
T.runTest()
@@ 357,7 +362,7 @@
def test_fitpack_wrapper(self):
""" run unit test contained in fitpack_wrapper.py
"""
 print "\n\nTESTING _fitpack_wrapper MODULE"
+ #print "\n\nTESTING _fitpack_wrapper MODULE"
from fitpack_wrapper import Test
T = Test()
T.runTest()
@@ 367,8 +372,7 @@
"""
# make sure : an instance of a callable class in which
 # x and y haven't been initiated works
 print 'hello'
+ # x and y haven't been initiated works
N = 7 #must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 402,9 +406,8 @@
self.assertAllclose(new_x, new_y)
def test_intper1d(self):
+ """ make sure : interp1d works, at least in the linear case
"""
 make sure : interp1d works, at least in the linear case
 """
N = 7
x = arange(N)
y = arange(N)
@@ 413,11 +416,10 @@
self.assertAllclose(new_x, new_y)
def test_spline1_defaultExt(self):
 """
 make sure : spline order 1 (linear) interpolation works correctly
+ """ make sure : spline order 1 (linear) interpolation works correctly
make sure : default extrapolation works
"""
 print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
+ #print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
N = 7 # must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 430,13 +432,12 @@
self.assert_(new_y[1] == 599.73)
def test_spline2(self):
 """
 make sure : order2 splines work on linear data
+ """ 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"
+ #print "\n\nTESTING 2nd ORDER SPLINE"
N = 7 #must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 462,11 +463,10 @@
def test_linear(self):
 """
 make sure : linear interpolation works
+ """ make sure : linear interpolation works
make sure : linear extrapolation works
"""
 print "\n\nTESTING LINEAR INTERPOLATION"
+ #print "\n\nTESTING LINEAR INTERPOLATION"
N = 7
x = arange(N)
y = arange(N)
@@ 482,11 +482,6 @@
self.assertAllclose(new_x, new_y)





if __name__ == '__main__':
unittest.main()
\ No newline at end of file
Modified: branches/Interpolate1D/TODO.txt
===================================================================
 branches/Interpolate1D/TODO.txt 20080725 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/TODO.txt 20080728 19:41:44 UTC (rev 4564)
@@ 5,26 +5,25 @@
at appropriate places in the code.
**include license for FITPACK
 take it from scipy.interpolate.
+************ PRESSING ISSUES ***********
+**handle smoothing
+ There is a question of whether, if, and how much to allow
+ smoothing of the data. If we smooth, we're not technically
+ interpolating, but users often want to smooth data.
**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.
+ In fitpack_wrapper the smoothing parameter is s. It now defaults
+ to 0.0 (exact interpolation). Zero smoothing and moderate (s ~ 1)
+ are fast, but small nonzero s makes the algorithm VERY slow.
 Then again, users very often want to smooth data at the same
 time. FITPACK suffers from VERY slow performance when s
 is very small but not 0. This is a problem.
+ Currently the user does see s in Interpolate1d; no smoothing
+ is available. The Spline class, however, has it available (default = 0.0).
**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 smoothing parameter s). Other code
+ scipy.interpolate). This isn't necessarily the best. 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.
@@ 35,15 +34,44 @@
**clean up fitpack_wrapper.py
 Currently is all hacked together from scipy.interpolate.
+ Currently it 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
+ we will include, and what should be accessible to the user,
+ that stuff needs to be cleaned up, and the
rest needs to be either removed or set aside.
+
+**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?
**comment interpolate1d
 There's comments there already, but they should be
+ There's also the question of Y being 2dimensional and/or
+ incorporating strings for record arrays. My instinct is to
+ have Interpolate1d only interpolate functions that are from
+ R1 > R1. That's how it is currently designed. Other stuff
+ can be made as wrappers around Interpolate1d.
+
+ See also "allow y to be 2dimensional?" below
+
+
+
+*********** DOCUENTATIONTYPE AND NONURGENT TASKS *******
+
+**improve regression tests
+ desired for fitpack_wrapper and _interpolate_wrapper
+ as well as Interpolate1d. What I have now is
+ really basic.
+
+
+**improve unit tests
+ interpolate1d.py has its own unit tests, which also call
+ the unit tests of fitpack_wrapper and _interpolate_wrapper.
+ But they could surely be made better.
+
+
+**comment all files
+ There are comments there already, but they should be
made better. Plus, features are changing, which requires
updating the documentation.
@@ 61,30 +89,42 @@
**figure out NumPy stuff with vectorize.
 In function interpolate1d.__call__
+ In function Interpolate1d.__call__
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.
+
+
+
+********* LONGER TERM ************
+**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.
+
+ There is an interesting problem here. Most of the extensions
+ I have assume a regular grid. First off, this is often not general.
+ Secondly, if I DO use a regular grid, how do I deal with bad
+ data? The best way is probably a preprocessing where you
+ interpolate values for the bad points (linear would be a nice simple
+ way to do it at first, just to get it working)
+
+ We should probably use something other than FITPACK for this.
+ First off it's at most 2D. But much worse, it doesn't evaluate at
+ a set of points; it evaluates over a grid, which requires inputs
+ being in sorted order (both in x and y coordinates). This makes
+ input inconvenient and the code runs a lot slower than ndimage.
+ ND image has 2 main downsides :1) it requires x and y be uniform
+ spacing (since its really interpolating entries of an array, rather
+ than values of a function on Rn), and 2) the C code is TERRIBLY
+ documented/commented. But that can be fixed.
+
+ So we may have two separate classes: Interpolate1d which is
+ based on FITPACK (or something else, depending on what
+ happens with the smoothing parameter), and InterpolateNd
+ which is based on ndimage.
**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?

 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.


**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
@@ 99,19 +139,6 @@
were stolen from.
**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.

 There is an interesting problem here. Most of the extensions
 I have assume a regular grid. First off, this is often not general.
 Secondly, if I DO use a regular grid, how do I deal with bad
 data? The best way is probably a preprocessing where you
 interpolate values for the bad points (linear would be a nice simple
 way to do it at first, just to get it working)


**allow y to be 2dimensional?
It is not decided whether this feature should be supported, but
I don't think it should; I think there should be another class with
Added: branches/Interpolate1D/example_script.py
===================================================================
 branches/Interpolate1D/example_script.py 20080725 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/example_script.py 20080728 19:41:44 UTC (rev 4564)
@@ 0,0 +1,41 @@
+# sample operation script
+import numpy as np
+import interpolate1d as I
+import matplotlib.pyplot as P
+
+
+N = 10.0
+x = np.arange(N)
+y = np.sin(x)
+
+
+## Interpolating inrange data
+newx = np.arange(.05, N, .05)
+
+# block interpolation
+interp = I.Interpolate1d(x, y, 'block')
+y_block = interp(newx)
+
+# linear interpolation
+interp = I.Interpolate1d(x, y, 'linear')
+y_linear = interp(newx)
+
+# 2nd order spline
+interp = I.Interpolate1d(x, y, 'quad')
+y_quad = interp(newx)
+
+# 3rd order spline
+interp = I.Interpolate1d(x, y, 'cubic')
+y_cubic = interp(newx)
+
+
+# plot result
+print "plotting results"
+P.hold(True)
+P.plot(newx, y_block)
+P.plot(newx, y_linear)
+P.plot(newx, y_quad)
+P.plot(newx, y_cubic)
+P.show()
+print "plotted results"
+
Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
 branches/Interpolate1D/interpolate1d.py 20080725 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/interpolate1d.py 20080728 19:41:44 UTC (rev 4564)
@@ 130,7 +130,7 @@
x  list or NumPy array
x includes the xvalues for the data set to
interpolate from. It must be sorted in
 ascending order
+ ascending order.
y  list or NumPy array
y includes the yvalues for the data set to
@@ 159,7 +159,7 @@
a number') for all values outside the range of x.
remove_bad_data  bool
 indicates whether to remove bad data.
+ indicates whether to remove bad data points from x and y.
bad_data  list
List of values (in x or y) which indicate unacceptable data. All points
@@ 168,7 +168,7 @@
numpy.NaN is always considered bad data.
 Acceptable Input Strings
+ Some Acceptable Input Strings

"linear"  linear interpolation : default
@@ 192,7 +192,8 @@
"""
# FIXME: more informative descriptions of sample arguments
# FIXME: examples in doc string
 # FIXME : Allow copying or not of arrays. noncopy + remove_bad_data should flash a warning (esp if we interpolate missing values), but work anyway.
+ # FIXME : Allow copying or not of arrays. noncopy + remove_bad_data should flash
+ # a warning (esp if we interpolate missing values), but work anyway.
def __init__(self, x, y, kind='linear', low=np.NaN, high=np.NaN, \
kindkw={}, lowkw={}, highkw={}, \
@@ 203,8 +204,7 @@
# check acceptable size and dimensions
x = np.array(x)
y = np.array(y)
 assert len(x) > 0 and len(y) > 0 , "Interpolate1D does not support\
 arrays of length 0"
+ assert len(x) > 0 and len(y) > 0 , "Arrays cannot be of zero length"
assert x.ndim == 1 , "x must be onedimensional"
assert y.ndim == 1 , "y must be onedimensional"
assert len(x) == len(y) , "x and y must be of the same length"
@@ 213,13 +213,14 @@
if remove_bad_data:
x, y = self._remove_bad_data(x, y, bad_data)
+ # store data
# FIXME : may be good to let x and y be initialized later, or changed afterthefact
self._init_xy(x, y)
# store interpolation functions for each range
 self.kind = self._init_interp_method(self._x, self._y, kind, kindkw)
 self.low = self._init_interp_method(self._x, self._y, low, lowkw)
 self.high = self._init_interp_method(self._x, self._y, high, highkw)
+ self.kind = self._init_interp_method(kind, kindkw)
+ self.low = self._init_interp_method(low, lowkw)
+ self.high = self._init_interp_method(high, highkw)
def _remove_bad_data(self, x, y, bad_data = [None, np.NaN]):
""" removes data points whose x or y coordinate is
@@ 242,7 +243,7 @@
self._x = make_array_safe(x, self._xdtype).copy()
self._y = make_array_safe(y, self._ydtype).copy()
 def _init_interp_method(self, x, y, interp_arg, kw):
+ def _init_interp_method(self, interp_arg, kw):
"""
User provides interp_arg and dictionary kw. _init_interp_method
returns the interpolating function from x and y specified by interp_arg,
@@ 295,7 +296,11 @@
result.set_xy(self._x, self._y)
# user passes a function to be called
 elif isfunction(interp_arg) and interp.func_code.argcount == 3:
+ # FIXME : I think there is too much flexibility allowed here; it makes
+ # there be more pathological side cases to consider. Functions
+ # should perhaps be reqired to be of the form f(x, y, newx, **kw)
+ elif isfunction(interp_arg) and interp.func_code.argcount >= 3:
+ # assume x, y and newx are all passed to interp_arg
result = lambda new_x : interp_arg(self._x, self._y, new_x, **kw)
elif isfunction(interp_arg):
result = lambda new_x : interp_arg(new_x, **kw)
@@ 306,7 +311,7 @@
return result
 def __call__(self, x):
+ def __call__(self, newx):
"""
Input x must be in sorted order.
Breaks x into pieces inrange, belowrange, and above range.
@@ 314,24 +319,24 @@
"""
# 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)
+ newx = make_array_safe(newx)
# masks indicate which elements fall into which interpolation region
 low_mask = x<self._x[0]
 high_mask = x>self._x[1]
+ low_mask = newx<self._x[0]
+ high_mask = newx>self._x[1]
interp_mask = (~low_mask) & (~high_mask)
# use correct function for x values in each region
 if len(x[low_mask]) == 0: new_low=np.array([]) # FIXME : remove need for if/else.
+ if len(newx[low_mask]) == 0: new_low=np.array([]) # FIXME : remove need for if/else.
# if/else is a hack, since vectorize is failing
# to work on lists/arrays of length 0
# on the computer where this is being
# developed
 else: new_low = self.low(x[low_mask])
 if len(x[interp_mask])==0: new_interp=np.array([])
 else: new_interp = self.kind(x[interp_mask])
 if len(x[high_mask]) == 0: new_high = np.array([])
 else: new_high = self.high(x[high_mask])
+ else: new_low = self.low(newx[low_mask])
+ if len(newx[interp_mask])==0: new_interp=np.array([])
+ else: new_interp = self.kind(newx[interp_mask])
+ if len(newx[high_mask]) == 0: new_high = np.array([])
+ else: new_high = self.high(newx[high_mask])
result = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
# Would be nice to say result = zeros(dtype=?)
@@ 349,7 +354,7 @@
def test_interpolate_wrapper(self):
""" run unit test contained in interpolate_wrapper.py
"""
 print "\n\nTESTING _interpolate_wrapper MODULE"
+ #print "\n\nTESTING _interpolate_wrapper MODULE"
from interpolate_wrapper import Test
T = Test()
T.runTest()
@@ 357,7 +362,7 @@
def test_fitpack_wrapper(self):
""" run unit test contained in fitpack_wrapper.py
"""
 print "\n\nTESTING _fitpack_wrapper MODULE"
+ #print "\n\nTESTING _fitpack_wrapper MODULE"
from fitpack_wrapper import Test
T = Test()
T.runTest()
@@ 367,8 +372,7 @@
"""
# make sure : an instance of a callable class in which
 # x and y haven't been initiated works
 print 'hello'
+ # x and y haven't been initiated works
N = 7 #must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 402,9 +406,8 @@
self.assertAllclose(new_x, new_y)
def test_intper1d(self):
+ """ make sure : interp1d works, at least in the linear case
"""
 make sure : interp1d works, at least in the linear case
 """
N = 7
x = arange(N)
y = arange(N)
@@ 413,11 +416,10 @@
self.assertAllclose(new_x, new_y)
def test_spline1_defaultExt(self):
 """
 make sure : spline order 1 (linear) interpolation works correctly
+ """ make sure : spline order 1 (linear) interpolation works correctly
make sure : default extrapolation works
"""
 print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
+ #print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
N = 7 # must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 430,13 +432,12 @@
self.assert_(new_y[1] == 599.73)
def test_spline2(self):
 """
 make sure : order2 splines work on linear data
+ """ 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"
+ #print "\n\nTESTING 2nd ORDER SPLINE"
N = 7 #must be > 5
x = np.arange(N)
y = np.arange(N)
@@ 462,11 +463,10 @@
def test_linear(self):
 """
 make sure : linear interpolation works
+ """ make sure : linear interpolation works
make sure : linear extrapolation works
"""
 print "\n\nTESTING LINEAR INTERPOLATION"
+ #print "\n\nTESTING LINEAR INTERPOLATION"
N = 7
x = arange(N)
y = arange(N)
@@ 482,11 +482,6 @@
self.assertAllclose(new_x, new_y)





if __name__ == '__main__':
unittest.main()
\ No newline at end of file
More information about the Scipysvn
mailing list