[Scipy-svn] r4538 - in branches: . interpolate2

scipy-svn@scip... scipy-svn@scip...
Thu Jul 10 18:29:26 CDT 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