[Scipy-svn] r4564 - branches/Interpolate1D
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Jul 28 15:41:45 EDT 2008
Author: fcady
Date: 2008-07-28 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 2008-07-25 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/Interpolate1D.py 2008-07-28 19:41:44 UTC (rev 4564)
@@ -130,7 +130,7 @@
x -- list or NumPy array
x includes the x-values for the data set to
interpolate from. It must be sorted in
- ascending order
+ ascending order.
y -- list or NumPy array
y includes the y-values 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. non-copy + remove_bad_data should flash a warning (esp if we interpolate missing values), but work anyway.
+ # FIXME : Allow copying or not of arrays. non-copy + 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 one-dimensional"
assert y.ndim == 1 , "y must be one-dimensional"
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 after-the-fact
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 in-range, below-range, 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 : order-2 splines work on linear data
+ """ 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"
+ #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 2008-07-25 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/TODO.txt 2008-07-28 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 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.
+ 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 non-zero 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
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 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, commented-out 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 2-dimensional 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 2-dimensional?" below
+
+
+
+*********** DOCUENTATION-TYPE AND NON-URGENT 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 pre-processing 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.
-
-
**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
@@ -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 pre-processing 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 2-dimensional?
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 2008-07-25 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/example_script.py 2008-07-28 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 in-range 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 2008-07-25 20:46:43 UTC (rev 4563)
+++ branches/Interpolate1D/interpolate1d.py 2008-07-28 19:41:44 UTC (rev 4564)
@@ -130,7 +130,7 @@
x -- list or NumPy array
x includes the x-values for the data set to
interpolate from. It must be sorted in
- ascending order
+ ascending order.
y -- list or NumPy array
y includes the y-values 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. non-copy + remove_bad_data should flash a warning (esp if we interpolate missing values), but work anyway.
+ # FIXME : Allow copying or not of arrays. non-copy + 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 one-dimensional"
assert y.ndim == 1 , "y must be one-dimensional"
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 after-the-fact
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 in-range, below-range, 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 : order-2 splines work on linear data
+ """ 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"
+ #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 Scipy-svn
mailing list