[Scipy-svn] r4538 - in branches: . interpolate2
scipy-svn at scipy.org
scipy-svn at scipy.org
Thu Jul 10 19:29:26 EDT 2008
Author: fcady
Date: 2008-07-10 18:29:12 -0500 (Thu, 10 Jul 2008)
New Revision: 4538
Added:
branches/interpolate2/
branches/interpolate2/Interpolate1D.py
branches/interpolate2/_interpolate.pyd
branches/interpolate2/interpolate_helper.py
Log:
added rudimentary but working new interpolation module
Added: branches/interpolate2/Interpolate1D.py
===================================================================
--- branches/interpolate2/Interpolate1D.py 2008-07-10 03:32:10 UTC (rev 4537)
+++ branches/interpolate2/Interpolate1D.py 2008-07-10 23:29:12 UTC (rev 4538)
@@ -0,0 +1,100 @@
+""" A module for intepolation
+"""
+
+#from enthought.interpolate import DataFit, Linear, Block
+from interpolate_helper import block, linear, block, Linear, Block, logarithmic
+from numpy import array, arange, empty, float64
+
+
+
+
+class Interpolate1D(object):
+ """ Class for interpolation and extrapolation of functions
+ from 1D data.
+
+ Parameters
+ ----------
+ x : array_like
+ x values to interpolate from. Pass in as NumPy array
+ or Python list.
+ y : array_like
+ y values to interpolate from. Pass in as NumPy array
+ or Python list.x and y must be the same length.
+ kind : optional
+ low : optional
+ high : optional
+
+
+ Example
+ -------
+ # fixme: This should also demonstrate extrapolation.
+
+ >>> from numpy import arange, array
+ >>> a = arange(10.)
+ >>> b = arange(10.) * 2.0
+ >>> new_a = array([1.0, 2.5, 8.0, 1.0, -1, 12])
+ >>> interp = Interpolate1D(a, b)
+ >>> interp(new_a)
+ array([ 2., 5., 16., 2., 0., 18.])
+
+ """
+
+ def __init__(self, x, y, kind='Linear', low='Block', high='Block'):
+ """ Class constructor
+
+ x : array-like
+
+ """
+
+ # fixme: Handle checking if they are the correct size.
+ self._x = array(x)
+ self._y = array(y)
+
+ self.kind = self._init_interp_method(x, y, kind)
+ self.low = self._init_interp_method(x, y, low)
+ self.high = self._init_interp_method(x, y, high)
+
+ def _init_interp_method(self, x, y, kind):
+ from inspect import isclass, isfunction
+
+ if isinstance(kind, str):
+ try:
+ exec('from interpolate_helper import %s; kind = %s' % (kind, kind) )
+ except:
+ raise
+
+ if isclass(kind):
+ result = kind(x, y)
+ elif isfunction(kind):
+ result = kind
+ else:
+ print 'WTF??'
+ raise
+ return result
+
+
+ def __call__(self, x):
+ low_mask = x<self._x[0]
+ high_mask = x>self._x[-1]
+ interp_mask = (~low_mask) & (~high_mask)
+
+ new_interp = self.kind(x[interp_mask])
+ new_low = self.low(x[low_mask])
+ new_high = self.high(x[high_mask])
+
+ # fixme: Handle other data types.
+ result = empty(x.shape, dtype=float64)
+ result[low_mask] = new_low
+ result[high_mask] = new_high
+ result[interp_mask] = new_interp
+
+ return result
+
+if __name__ == '__main__':
+ a = arange(10.)
+ b = 2*a
+ c = array([-1, 4.5, 19])
+ interp = Interpolate1D(a, b, kind=lambda x: linear(a,b,x), \
+ low='Block', high=lambda x: block(a,b,x) )
+ print 'c equals: ', c
+ print 'interp(c) equals: ', interp(c)
\ No newline at end of file
Added: branches/interpolate2/_interpolate.pyd
===================================================================
(Binary files differ)
Property changes on: branches/interpolate2/_interpolate.pyd
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: branches/interpolate2/interpolate_helper.py
===================================================================
--- branches/interpolate2/interpolate_helper.py 2008-07-10 03:32:10 UTC (rev 4537)
+++ branches/interpolate2/interpolate_helper.py 2008-07-10 23:29:12 UTC (rev 4538)
@@ -0,0 +1,156 @@
+""" helper_funcs.py
+"""
+
+import numpy
+import _interpolate
+print "this is the module"
+
+def make_array_safe(ary, typecode):
+ ary = numpy.atleast_1d(numpy.asarray(ary, typecode))
+ if not ary.flags['CONTIGUOUS']:
+ ary = ary.copy()
+ return ary
+
+
+class Block():
+ """ Used when only one element is available in the log.
+ """
+ def __init__(self, x, y):
+ self._x = x
+ self._y = y
+
+ def __call__(self, new_x):
+ # find index of values in x that preceed values in x
+ # This code is a little strange -- we really want a routine that
+ # returns the index of values where x[j] < x[index]
+ TINY = 1e-10
+ indices = numpy.searchsorted(new_x, new_x+TINY)-1
+
+ # If the value is at the front of the list, it'll have -1.
+ # In this case, we will use the first (0), element in the array.
+ # take requires the index array to be an Int
+ indices = numpy.atleast_1d(numpy.clip(indices, 0, numpy.Inf).astype(numpy.int))
+ result = numpy.take(self._y, indices, axis=-1)
+ return result
+
+class Linear():
+ """ Linearly interpolates values in new_x based on the values in x and y
+
+ Parameters
+ ----------
+ x
+ 1-D array
+ y
+ 1-D or 2-D array
+ new_x
+ 1-D array
+ """
+ def __init__(self, x, y):
+ self._x = make_array_safe(x, numpy.float64)
+ self._y = make_array_safe(y, numpy.float64)
+
+ assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
+
+ def __call__(self, new_x):
+ new_x = make_array_safe(new_x, numpy.float64)
+ if len(self._y.shape) == 2:
+ new_y = numpy.zeros((y.shape[0], len(new_x)), numpy.float64)
+ for i in range(len(new_y)):
+ _interpolate.linear_dddd(x, y[i], new_x, new_y[i])
+ else:
+ new_y = numpy.zeros(len(new_x), numpy.float64)
+ _interpolate.linear_dddd(x, y, new_x, new_y)
+ return new_y
+
+def block(x, y, new_x):
+ """ Used when only one element is available in the log.
+ """
+
+ # find index of values in x that preceed values in x
+ # This code is a little strange -- we really want a routine that
+ # returns the index of values where x[j] < x[index]
+ TINY = 1e-10
+ indices = numpy.searchsorted(new_x, new_x+TINY)-1
+
+ # If the value is at the front of the list, it'll have -1.
+ # In this case, we will use the first (0), element in the array.
+ # take requires the index array to be an Int
+ indices = numpy.atleast_1d(numpy.clip(indices, 0, numpy.Inf).astype(numpy.int))
+ result = numpy.take(y, indices, axis=-1)
+ return result
+
+def linear(x, y, new_x):
+ """ Linearly interpolates values in new_x based on the values in x and y
+
+ Parameters
+ ----------
+ x
+ 1-D array
+ y
+ 1-D or 2-D array
+ new_x
+ 1-D array
+ """
+ x = make_array_safe(x, numpy.float64)
+ y = make_array_safe(y, numpy.float64)
+ new_x = make_array_safe(new_x, numpy.float64)
+
+ assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
+ if len(y.shape) == 2:
+ new_y = numpy.zeros((y.shape[0], len(new_x)), numpy.float64)
+ for i in range(len(new_y)):
+ _interpolate.linear_dddd(x, y[i], new_x, new_y[i])
+ else:
+ new_y = numpy.zeros(len(new_x), numpy.float64)
+ _interpolate.linear_dddd(x, y, new_x, new_y)
+
+ return new_y
+
+class Logarithmic:
+ """ For log-linear interpolation
+ """
+
+ def __init__(self, x, y):
+ self._x = make_array_safe(x, numpy.float64)
+ self._y = make_array_safe(y, numpy.float64)
+
+ def __call__(new_x):
+ new_x = make_array_safe(new_x, numpy.float64)
+ assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
+ if len(y.shape) == 2:
+ new_y = numpy.zeros((y.shape[0], len(new_x)), numpy.float64)
+ for i in range(len(new_y)):
+ _interpolate.loginterp_dddd(x, y[i], new_x, new_y[i])
+ else:
+ new_y = numpy.zeros(len(new_x), numpy.float64)
+ _interpolate.loginterp_dddd(x, y, new_x, new_y)
+
+ return new_y
+
+def logarithmic(x, y, new_x):
+ """ Linearly interpolates values in new_x based in the log space of y.
+
+ Parameters
+ ----------
+ x
+ 1-D array
+ y
+ 1-D or 2-D array
+ new_x
+ 1-D array
+ """
+ x = make_array_safe(x, numpy.float64)
+ y = make_array_safe(y, numpy.float64)
+ new_x = make_array_safe(new_x, numpy.float64)
+
+ assert len(y.shape) < 3, "function only works with 1D or 2D arrays"
+ if len(y.shape) == 2:
+ new_y = numpy.zeros((y.shape[0], len(new_x)), numpy.float64)
+ for i in range(len(new_y)):
+ _interpolate.loginterp_dddd(x, y[i], new_x, new_y[i])
+ else:
+ new_y = numpy.zeros(len(new_x), numpy.float64)
+ _interpolate.loginterp_dddd(x, y, new_x, new_y)
+
+ return new_y
+
More information about the Scipy-svn
mailing list