[Scipy-svn] r4577 - branches/Interpolate1D
scipy-svn@scip...
scipy-svn@scip...
Tue Jul 29 12:56:53 CDT 2008
Author: fcady
Date: 2008-07-29 12:56:52 -0500 (Tue, 29 Jul 2008)
New Revision: 4577
Modified:
branches/Interpolate1D/TODO.txt
branches/Interpolate1D/__init__.py
branches/Interpolate1D/erics_notes.txt
branches/Interpolate1D/fitpack_wrapper.py
branches/Interpolate1D/info.py
branches/Interpolate1D/interpolate1d.py
branches/Interpolate1D/interpolate_wrapper.py
Log:
incorporated many of Eric's suggestions
Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/TODO.txt 2008-07-29 17:56:52 UTC (rev 4577)
@@ -103,7 +103,7 @@
**update for 2D and ND
This will probably take the form of two additional
- classes both based on interpolate1d. Thus it probably
+ classes both modelled after interpolate1d. Thus it probably
shouldn't be done until interpolate1d is more settled.
There is an interesting problem here. Most of the extensions
@@ -114,7 +114,7 @@
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
+ But firstly it's at most 2D. 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.
@@ -127,6 +127,13 @@
based on FITPACK (or something else, depending on what
happens with the smoothing parameter), and InterpolateNd
which is based on ndimage.
+
+ Another option is to have two classes: one for uniformly spaced data and
+ and one for scattered data. Regularly spaced would use NDImage, and
+ scattered would start as an inefficient wrapper around FITPACK. But longer
+ term the scattered data could be done with Delaunay triangulation,
+ perhaps something that would implicitly calculate the convex hull
+ of the points and interpolate within it.
**high-level road map
when the module is more established, there should be a page on
Modified: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/__init__.py 2008-07-29 17:56:52 UTC (rev 4577)
@@ -9,4 +9,4 @@
from fitpack_wrapper import Spline
# wrapping for all supported interpolation types.
-from interpolate1d import Interpolate1d, interp1d
\ No newline at end of file
+from interpolate1d import Interpolate1d, interp1
\ No newline at end of file
Modified: branches/Interpolate1D/erics_notes.txt
===================================================================
--- branches/Interpolate1D/erics_notes.txt 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/erics_notes.txt 2008-07-29 17:56:52 UTC (rev 4577)
@@ -1,9 +1,39 @@
-*. I am glad to see your docstrings.
+*. Function signatures:
+
+ def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
+ kindkw={}, lowkw={}, highkw={}, \
+ remove_bad_data = False, bad_data=[], interp_axis = 0):
-*. Move tests into a seprate test directory.
+ o. the trailing \ for new lines are not necessary for functions since line continuation
+ is implicit with the open/close parentheses.
+ o. [minor] I would just use NaN instead of np.NaN.
+ o. It is dangerous to initialize containers as default arguments because they are
+ effectively a singleton for the function. We can discuss this if you don't know
+ what I am talking about.
+ o. kindkw, lowkw, and highkw aren't really necessary I don't think. They should go.
+ Xo. Do you need both remove_bad_data and bad_data? If bad_data is None, then you
+ don't remove bad_data...
+ Xo. I think I would change interp_axis to just be axis. This is consistent with many
+ other functions in numpy.
+ Xo. The choice of whether axis=0 or axis=-1 by default is a reasonable question.
+ fft defaults to axis=-1. This is also the faster axis to operate across in the
+ standard case. It is, however, the opposite of how some people think about
+ things (columns vs. rows). Talk to Travis O. for his take. Mine is to use axis=-1.
+
+ I think all of this might simplify the interface to the following:
+
+ def interp1d(x, y, new_x, kind='linear', low=NaN, high=NaN, bad_data=None, axis=-1):
-*. Follow scipy/FORMAT_GUIDLINES.txt in the main scipy directory.
+*. _remove_bad_data should be left up to the interpolation method if it knows what to do.
+ otherwise, it is handled by this top level class.
+ And, we definitely don't want the list comprehension in the remove_bad_data class.
+X*. I am glad to see your docstrings.
+
+X*. Move tests into a seprate test directory.
+
+X*. Follow scipy/FORMAT_GUIDLINES.txt in the main scipy directory.
+
For example:
test_callFormat -> test_call_format
@@ -14,14 +44,14 @@
describes them "here":http://www.python.org/doc/essays/styleguide.html. A few
reminders follow:
- o Use 4 spaces for indentation levels. Do not use tabs as they can result
+ Xo Use 4 spaces for indentation levels. Do not use tabs as they can result
in indentation confusion. Most editors have a feature that will insert 4
spaces when the tab key is hit. Also, many editors will automatically
search/replace leading tabs with 4 spaces.
o Only 80 characters on a line.
- o use all lowercase function names with underscore separated words:
+ Xo use all lowercase function names with underscore separated words:
def set_some_value()
@@ -30,7 +60,7 @@
def setSomeValue()
- o use CamelCase class names:
+ Xo use CamelCase class names:
def BaseClass()
@@ -38,7 +68,7 @@
def base_class()
-*. make_array_safe
+X*. make_array_safe
I would prefer 'make_array_safe' named atleast_1d_and_contiguous(). This is more specific
and it is immediately clear to other developers what the check does. If you add other checks,
then perhaps come up with a more generic name, but this explicit names, when possible, help
@@ -47,42 +77,11 @@
Also, numpy has an ascontiguousrarray() function that would simplify the code a line or so.
Also, this function lives in multiple places, interpolate1d and interpolate_wrapper.
-
-*. Function signatures:
-
- def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
- kindkw={}, lowkw={}, highkw={}, \
- remove_bad_data = False, bad_data=[], interp_axis = 0):
- o. the trailing \ for new lines are not necessary for functions since line continuation
- is implicit with the open/close parentheses.
- o. [minor] I would just use NaN instead of np.NaN.
- o. It is dangerous to initialize containers as default arguments because they are
- effectively a singleton for the function. We can discuss this if you don't know
- what I am talking about.
- o. kindkw, lowkw, and highkw aren't really necessary I don't think. They should go.
- o. Do you need both remove_bad_data and bad_data? If bad_data is None, then you
- don't remove bad_data...
- o. I think I would change interp_axis to just be axis. This is consistent with many
- other functions in numpy.
- o. The choice of whether axis=0 or axis=-1 by default is a reasonable question.
- fft defaults to axis=-1. This is also the faster axis to operate across in the
- standard case. It is, however, the opposite of how some people think about
- things (columns vs. rows). Talk to Travis O. for his take. Mine is to use axis=-1.
- o.
-
- I think all of this might simplify the interface to the following:
-
- def interp1d(x, y, new_x, kind='linear', low=NaN, high=NaN, bad_data=None, axis=-1):
-
-*. isinstance(value, str) should be isinstance(value, basestring) so that we handle
+X*. isinstance(value, str) should be isinstance(value, basestring) so that we handle
both strings and unicode correctly.
-*. _remove_bad_data should be left up to the interpolation method if it knows what to do.
- otherwise, it is handled by this top level class.
- And, we definitely don't want the list comprehension in the remove_bad_data class.
-
-*. If the input to interp1d is a scalar, the return value should be a scalar.
+X*. If the input to interp1d is a scalar, the return value should be a scalar.
[add test to handle this.]
This fails in the following:
Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/fitpack_wrapper.py 2008-07-29 17:56:52 UTC (rev 4577)
@@ -174,74 +174,4 @@
assert ier==0,`ier`
return z[:m]
raise NotImplementedError,\
- 'finding roots unsupported for non-cubic splines'
-
-
-
-# testing
-import unittest
-import time
-from numpy import arange, allclose, ones
-
-class Test(unittest.TestCase):
-
- def assertAllclose(self, x, y):
- self.assert_(np.allclose(x, y))
-
- def test_linearInterp(self):
- """ make sure : linear interpolation (spline with order = 1, s = 0)works
- """
- N = 3000.
- x = np.arange(N)
- y = np.arange(N)
- #T1 = time.clock()
- interp_func = Spline(x, y, k=1)
- #T2 = time.clock()
- #print "time to create order 1 spline interpolation function with N = %i:" % N, T2 - T1
- new_x = np.arange(N)+0.5
- #t1 = time.clock()
- new_y = interp_func(new_x)
- #t2 = time.clock()
- #print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
- self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
-
- def test_quadInterp(self):
- """ make sure : quadratic interpolation (spline with order = 2, s = 0)works
- """
- N = 3000.
- x = np.arange(N)
- y = x**2
- interp_func = Spline(x, y, k=2)
- #print "time to create order 1 spline interpolation function with N = %i:" % N, T2 - T1
- new_x = np.arange(N)+0.5
- #t1 = time.clock()
- new_y = interp_func(x)
- #t2 = time.clock()
- #print "time for order 1 spline interpolation with N = %i:" % N, t2 - t1
- self.assertAllclose(new_y, y)
-
-
- def test_inputFormat(self):
- """ make sure : it's possible to instantiate Spline without x and y
- """
- #print "testing input format"
- N = 3000.
- x = np.arange(N)
- y = np.arange(N)
- interp_func = Spline(k=1)
- interp_func.init_xy(x, y)
- new_x = np.arange(N)+0.5
- new_y = interp_func(new_x)
- self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
-
- def runTest(self):
- test_list = [name for name in dir(self) if name.find('test_')==0]
- for test_name in test_list:
- exec("self.%s()" % test_name)
-
-
-
-if __name__ == '__main__':
- unittest.main()
-
-
\ No newline at end of file
+ 'finding roots unsupported for non-cubic splines'
\ No newline at end of file
Modified: branches/Interpolate1D/info.py
===================================================================
--- branches/Interpolate1D/info.py 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/info.py 2008-07-29 17:56:52 UTC (rev 4577)
@@ -9,7 +9,7 @@
and extrapolation of 1D (in both input and output) real-valued. The
primary function provided is:
- interp1d(x, y, new_x) : from data points (x[i], y[i]), interpolates
+ interp1(x, y, new_x) : from data points (x[i], y[i]), interpolates
values for points in new_x and
returns them as an array. x and new_x
must both be in sorted order.
Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/interpolate1d.py 2008-07-29 17:56:52 UTC (rev 4577)
@@ -1,24 +1,14 @@
-
# FIXME: information strings giving mathematical descriptions of the actions
# of the functions.
-from interpolate_wrapper import linear, logarithmic, block, block_average_above
+from interpolate_wrapper import linear, logarithmic, block, block_average_above, atleast_1d_and_contiguous
from fitpack_wrapper import Spline
import numpy as np
from numpy import array, arange, empty, float64, NaN
-
-def make_array_safe(ary, typecode=np.float64):
- """Used to make sure that inputs and outputs are
- properly formatted.
- """
- ary = np.atleast_1d(np.asarray(ary, typecode))
- if not ary.flags['CONTIGUOUS']:
- ary = ary.copy()
- return ary
-def interp1d(x, y, new_x, kind='linear', low=np.NaN, high=np.NaN, \
- kindkw={}, lowkw={}, highkw={}, \
- remove_bad_data = False, bad_data=[], interp_axis = 0):
+def interp1d(x, y, new_x, interp = 'linear', low = NaN, high = NaN,
+ interpkw = {}, lowkw={}, highkw={},
+ bad_data = None):
""" A function for interpolation of 1D data.
Parameters
@@ -90,10 +80,15 @@
>>> interp1d(x, y, new_x)
array([.2, 2.3, 5.6, NaN])
"""
- 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)
+ return Interpolate1d(x, y,
+ interp = interp,
+ low = low,
+ high = high,
+ interpkw = interpkw,
+ lowkw = lowkw,
+ highkw = highkw,
+ bad_data = bad_data
+ )(new_x)
class Interpolate1d(object):
""" A callable class for interpolation of 1D, real-valued data.
@@ -192,11 +187,12 @@
---------
>>> import numpy
- >>> from Interpolate1D import interp1d
+ >>> from interpolate1d import Interpolate1d
>>> x = range(5) # note list is permitted
>>> y = numpy.arange(5.)
>>> new_x = [.2, 2.3, 5.6]
- >>> interp1d(x, y, new_x)
+ >>> interp_func = Interpolate1d(x, y)
+ >>> interp_fuc(new_x)
array([.2, 2.3, 5.6, NaN])
"""
# FIXME: more informative descriptions of sample arguments
@@ -204,19 +200,18 @@
# 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={}, \
- remove_bad_data = False, bad_data=[]):
+ def __init__(self, x, y, interp = 'linear', low = NaN, high = NaN,
+ interpkw={}, lowkw={}, highkw={}, bad_data = None):
# FIXME: don't allow copying multiple times.
# FIXME : allow no copying, in case user has huge dataset
# remove bad data, is there is any
- if remove_bad_data:
+ if bad_data is not None:
x, y = self._remove_bad_data(x, y, bad_data)
# check acceptable size and dimensions
- x = np.array(x)
- y = np.array(y)
+ x = np.atleast_1d(x)
+ y = np.atleast_1d(y)
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"
@@ -227,7 +222,7 @@
self._init_xy(x, y)
# store interpolation functions for each range
- self.kind = self._init_interp_method(kind, kindkw)
+ self.interp = self._init_interp_method(interp, interpkw)
self.low = self._init_interp_method(low, lowkw)
self.high = self._init_interp_method(high, highkw)
@@ -236,10 +231,10 @@
# select proper dataypes and make arrays
self._xdtype = {np.float32 : np.float32}.setdefault(type(x[0]), np.float64) # unless data is float32, cast to float64
self._ydtype = {np.float32 : np.float32}.setdefault(type(y[0]), np.float64)
- self._x = make_array_safe(x, self._xdtype).copy()
- self._y = make_array_safe(y, self._ydtype).copy()
+ self._x = atleast_1d_and_contiguous(x, self._xdtype).copy()
+ self._y = atleast_1d_and_contiguous(y, self._ydtype).copy()
- def _remove_bad_data(self, x, y, bad_data = [None, np.NaN]):
+ def _remove_bad_data(self, x, y, bad_data = [None, NaN]):
""" removes data points whose x or y coordinate is
either in bad_data or is a NaN.
"""
@@ -272,7 +267,7 @@
func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
'block_average_above':block_average_above}[interp_arg]
result = lambda new_x : func(self._x, self._y, new_x, **kw)
- elif interp_arg in ['Spline', Spline, 'spline']:
+ elif interp_arg in ['Spline', 'spline']:
# use the Spline class from fitpack_wrapper
# k = 3 unless otherwise specified
result = Spline(self._x, self._y, **kw)
@@ -289,7 +284,7 @@
result = Spline(self._x, self._y, k=4)
elif interp_arg in ['Quintic', 'quintic', 'Quin', 'quin']:
result = Spline(self._x, self._y, k=5)
- elif isinstance(interp_arg, str):
+ elif isinstance(interp_arg, basestring):
raise TypeError, "input string %s not valid" % interp_arg
# secondary usage : user passes a callable class
@@ -314,9 +309,7 @@
# user passes a function to be called
# Assume function has form of f(x, y, newx, **kw)
- # FIXME : should other function forms be allowed?
elif isfunction(interp_arg):
- # assume x, y and newx are all passed to interp_arg
result = lambda new_x : interp_arg(self._x, self._y, new_x, **kw)
# default : user has passed a default value to always be returned
@@ -332,171 +325,39 @@
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.
+ # FIXME : atleast_1d_and_contiguous may also be called within the interpolation technique.
# waste of time, but ok for the time being.
- newx = make_array_safe(newx)
+ # if input is scalar or 0-dimemsional array, output will be scalar
+ input_is_scalar = np.isscalar(newx) or (isinstance(newx, type(np.array([1.0]))) and np.shape(newx) == ())
+
+ newx_array = atleast_1d_and_contiguous(newx)
+
# masks indicate which elements fall into which interpolation region
- low_mask = newx<self._x[0]
- high_mask = newx>self._x[-1]
+ low_mask = newx_array<self._x[0]
+ high_mask = newx_array>self._x[-1]
interp_mask = (~low_mask) & (~high_mask)
-
+
# use correct function for x values in each region
- if len(newx[low_mask]) == 0: new_low=np.array([]) # FIXME : remove need for if/else.
+ if len(newx_array[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(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])
+ else: new_low = self.low(newx_array[low_mask])
+ if len(newx_array[interp_mask])==0: new_interp=np.array([])
+ else: new_interp = self.interp(newx_array[interp_mask])
+ if len(newx_array[high_mask]) == 0: new_high = np.array([])
+ else: new_high = self.high(newx_array[high_mask])
- result = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
+ result_array = np.concatenate((new_low, new_interp, new_high)) # FIXME : deal with mixed datatypes
# Would be nice to say result = zeros(dtype=?)
# and fill in
- return result
+ if input_is_scalar:
+ result = float(result_array)
+ else:
+ result = result_array
-# unit testing
-import unittest, time
-class Test(unittest.TestCase):
-
- def assertAllclose(self, x, y):
- self.assert_(np.allclose(make_array_safe(x), make_array_safe(y)))
-
- def test_interpolate_wrapper(self):
- """ run unit test contained in interpolate_wrapper.py
- """
- #print "\n\nTESTING _interpolate_wrapper MODULE"
- from interpolate_wrapper import Test
- T = Test()
- T.runTest()
-
- def test_fitpack_wrapper(self):
- """ run unit test contained in fitpack_wrapper.py
- """
- #print "\n\nTESTING _fitpack_wrapper MODULE"
- from fitpack_wrapper import Test
- T = Test()
- T.runTest()
-
- def test_instantiationFormat(self):
- """ make sure : all allowed instantiation formats are supported
- """
-
- # make sure : an instance of a callable class in which
- # x and y haven't been initiated works
- N = 7 #must be > 5
- x = np.arange(N)
- y = np.arange(N)
- interp_func = Interpolate1d(x, y, kind=Spline(k=2), low=Spline(k=2), high=Spline(k=2))
- new_x = np.arange(N+1)-0.5
- new_y = interp_func(new_x)
- self.assertAllclose(new_x, new_y)
-
- def test_callFormat(self):
- """ make sure : all allowed calling formats are supported
- """
- # make sure : having no out-of-range elements in new_x is fine
- # There was a bug with this earlier.
- N = 5
- x = arange(N)
- y = arange(N)
- new_x = arange(1,N-1)+.2
- interp_func = Interpolate1d(x, y, kind='linear', low='linear', high=np.NaN)
- new_y = interp_func(new_x)
- self.assertAllclose(new_x, new_y)
-
- def test_removeBad(self):
- """make sure : interp1d works with bad data
- """
- N = 7.0 # must be >=5
- x = arange(N); x[2] = np.NaN
- y = arange(N); y[4] = None; y[0]=np.NaN
- new_x = arange(N+1)-0.5
- new_y = interp1d(x, y, new_x, kind='linear', low='linear', high='linear', \
- remove_bad_data = True, bad_data = [None])
- self.assertAllclose(new_x, new_y)
-
- def test_intper1d(self):
- """ make sure : interp1d works, at least in the linear case
- """
- N = 7
- x = arange(N)
- y = arange(N)
- new_x = arange(N+1)-0.5
- new_y = interp1d(x, y, new_x, kind='linear', low='linear', high='linear')
- self.assertAllclose(new_x, new_y)
-
- def test_spline1_defaultExt(self):
- """ make sure : spline order 1 (linear) interpolation works correctly
- make sure : default extrapolation works
- """
- #print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
- N = 7 # must be > 5
- x = np.arange(N)
- y = np.arange(N)
- interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':1}, low=None, high=599.73)
- new_x = np.arange(N+1)-0.5
- new_y = interp_func(new_x)
-
- self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
- self.assert_(new_y[0] == None)
- 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 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
- x = np.arange(N)
- y = np.arange(N)
- T1 = time.clock()
- interp_func = Interpolate1d(x, y, kind='Spline', kindkw={'k':2}, low='spline', high='spline')
- T2 = time.clock()
- print "time to create 2nd order spline interp function with N = %i: " % N, T2 - T1
- new_x = np.arange(N+1)-0.5
- t1 = time.clock()
- new_y = interp_func(new_x)
- t2 = time.clock()
- print "time to evaluate 2nd order spline interp function with N = %i: " % N, t2 - t1
- self.assertAllclose(new_x, new_y)
-
- # make sure for non-linear data
- N = 7
- x = np.arange(N)
- y = x**2
- 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)
-
-
- def test_linear(self):
- """ make sure : linear interpolation works
- make sure : linear extrapolation works
- """
- #print "\n\nTESTING LINEAR INTERPOLATION"
- N = 7
- x = arange(N)
- y = arange(N)
- new_x = arange(N+1)-0.5
- T1 = time.clock()
- interp_func = Interpolate1d(x, y, kind='linear', low='linear', high='linear')
- T2 = time.clock()
- print "time to create linear interp function with N = %i: " % N, T2 - T1
- t1 = time.clock()
- new_y = interp_func(new_x)
- t2 = time.clock()
- print "time to create linear interp function with N = %i: " % N, t2 - t1
-
- self.assertAllclose(new_x, new_y)
-
-
-if __name__ == '__main__':
- unittest.main()
\ No newline at end of file
+ return result
+
\ No newline at end of file
Modified: branches/Interpolate1D/interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/interpolate_wrapper.py 2008-07-29 14:11:52 UTC (rev 4576)
+++ branches/Interpolate1D/interpolate_wrapper.py 2008-07-29 17:56:52 UTC (rev 4577)
@@ -6,13 +6,9 @@
import sys
import _interpolate # C extension. Does all the real work.
-def make_array_safe(ary, typecode = np.float64):
- ary = np.atleast_1d(np.asarray(ary, typecode))
- if not ary.flags['CONTIGUOUS']:
- ary = ary.copy()
- return ary
+def atleast_1d_and_contiguous(ary, typecode = np.float64):
+ return np.atleast_1d( np.ascontiguousarray(ary, typecode) )
-
def linear(x, y, new_x):
""" Linearly interpolates values in new_x based on the values in x and y
@@ -25,9 +21,9 @@
new_x
1-D array
"""
- x = make_array_safe(x, np.float64)
- y = make_array_safe(y, np.float64)
- new_x = make_array_safe(new_x, np.float64)
+ x = atleast_1d_and_contiguous(x, np.float64)
+ y = atleast_1d_and_contiguous(y, np.float64)
+ new_x = atleast_1d_and_contiguous(new_x, np.float64)
assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
if len(y.shape) == 2:
@@ -52,9 +48,9 @@
new_x
1-D array
"""
- x = make_array_safe(x, np.float64)
- y = make_array_safe(y, np.float64)
- new_x = make_array_safe(new_x, np.float64)
+ x = atleast_1d_and_contiguous(x, np.float64)
+ y = atleast_1d_and_contiguous(y, np.float64)
+ new_x = atleast_1d_and_contiguous(new_x, np.float64)
assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
if len(y.shape) == 2:
@@ -80,9 +76,9 @@
1-D array
"""
bad_index = None
- x = make_array_safe(x, np.float64)
- y = make_array_safe(y, np.float64)
- new_x = make_array_safe(new_x, np.float64)
+ x = atleast_1d_and_contiguous(x, np.float64)
+ y = atleast_1d_and_contiguous(y, np.float64)
+ new_x = atleast_1d_and_contiguous(new_x, np.float64)
assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
if len(y.shape) == 2:
@@ -122,77 +118,4 @@
# take requires the index array to be an Int
indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(np.int))
new_y = np.take(y, indices, axis=-1)
- return new_y
-
-
-# Unit Test
-import unittest
-import time
-from numpy import arange, allclose, ones, NaN, isnan
-class Test(unittest.TestCase):
-
-
- def assertAllclose(self, x, y, rtol=1.0e-5):
- for i, xi in enumerate(x):
- self.assert_(allclose(xi, y[i], rtol) or (isnan(xi) and isnan(y[i])))
-
- def test_linear(self):
- N = 3000.
- x = arange(N)
- y = arange(N)
- new_x = arange(N)+0.5
- t1 = time.clock()
- new_y = linear(x, y, new_x)
- t2 = time.clock()
- #print "time for linear interpolation with N = %i:" % N, t2 - t1
-
- self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
-
- def test_block_average_above(self):
- N = 3000.
- x = arange(N)
- y = arange(N)
-
- new_x = arange(N/2)*2
- t1 = time.clock()
- new_y = block_average_above(x, y, new_x)
- t2 = time.clock()
- #print "time for block_avg_above interpolation with N = %i:" % N, t2 - t1
- self.assertAllclose(new_y[:5], [0.0, 0.5, 2.5, 4.5, 6.5])
-
- def test_linear2(self):
- N = 3000.
- x = arange(N)
- y = ones((100,N)) * arange(N)
- new_x = arange(N)+0.5
- t1 = time.clock()
- new_y = linear(x, y, new_x)
- t2 = time.clock()
- #print "time for 2D linear interpolation with N = %i:" % N, t2 - t1
- self.assertAllclose(new_y[:5,:5],
- [[ 0.5, 1.5, 2.5, 3.5, 4.5],
- [ 0.5, 1.5, 2.5, 3.5, 4.5],
- [ 0.5, 1.5, 2.5, 3.5, 4.5],
- [ 0.5, 1.5, 2.5, 3.5, 4.5],
- [ 0.5, 1.5, 2.5, 3.5, 4.5]])
-
- def test_logarithmic(self):
- N = 4000.
- x = arange(N)
- y = arange(N)
- new_x = arange(N)+0.5
- t1 = time.clock()
- new_y = logarithmic(x, y, new_x)
- t2 = time.clock()
- #print "time for logarithmic interpolation with N = %i:" % N, t2 - t1
- correct_y = [np.NaN, 1.41421356, 2.44948974, 3.46410162, 4.47213595]
- self.assertAllclose(new_y[:5], correct_y)
-
- def runTest(self):
- test_list = [name for name in dir(self) if name.find('test_')==0]
- for test_name in test_list:
- exec("self.%s()" % test_name)
-
-if __name__ == '__main__':
- unittest.main()
-
\ No newline at end of file
+ return new_y
\ No newline at end of file
More information about the Scipy-svn
mailing list