[Scipy-svn] r4542 - branches/interpolate2

scipy-svn@scip... scipy-svn@scip...
Mon Jul 14 09:28:52 CDT 2008


Author: fcady
Date: 2008-07-14 09:28:50 -0500 (Mon, 14 Jul 2008)
New Revision: 4542

Added:
   branches/interpolate2/journal.txt
Modified:
   branches/interpolate2/Interpolate1D.py
   branches/interpolate2/interpolate_helper.py
Log:
commiting so that I can access from office computer

Modified: branches/interpolate2/Interpolate1D.py
===================================================================
--- branches/interpolate2/Interpolate1D.py	2008-07-14 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/Interpolate1D.py	2008-07-14 14:28:50 UTC (rev 4542)
@@ -2,7 +2,8 @@
 """
 
 #from enthought.interpolate import DataFit, Linear, Block
-from interpolate_helper import block, linear, block, Linear, Block, logarithmic
+from interpolate_helper import * #block, linear, block, Linear, \
+    #Block, logarithmic, Block_Average_Above, Window_Average
 from numpy import array, arange, empty, float64
 
 
@@ -50,25 +51,25 @@
         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)
+        self.kind = self._init_interp_method(self._x, self._y, kind)
+        self.low = self._init_interp_method(self._x, self._y, low)
+        self.high = self._init_interp_method(self._x, self._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) )
+                #fixme: import the object
+                exec('dummy_var = %s' % kind)
             except:
-                raise
+                raise "Could not import indicated class or function"
         
         if isclass(kind):
             result = kind(x, y)
         elif isfunction(kind):
             result = kind
         else:
-            print 'WTF??'
             raise
         return result
     
@@ -89,12 +90,19 @@
         result[interp_mask] = new_interp
         
         return result
+        
+    def __clone__(self):
+        # fixme: objects are currently double-referenced.  Make them copied over.
+        clone_x = self._x.clone()
+        clone_y = self._y.clone()
+        return Interpolate1D(self._x, self._y, low=self.low, kind=self.kind, high=self.high)
 
 if __name__ == '__main__':
+    print "hey world"
     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) )
+    interp = Interpolate1D(a, b, low = 'Block', \
+        kind='Window_Average', high=lambda x: block(a,b,x) )
     print 'c equals: ', c
     print 'interp(c) equals: ', interp(c)                   
\ No newline at end of file

Modified: branches/interpolate2/interpolate_helper.py
===================================================================
--- branches/interpolate2/interpolate_helper.py	2008-07-14 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/interpolate_helper.py	2008-07-14 14:28:50 UTC (rev 4542)
@@ -3,7 +3,6 @@
 
 import numpy
 import _interpolate
-print "this is the module"
 
 def make_array_safe(ary, typecode):
     ary = numpy.atleast_1d(numpy.asarray(ary, typecode))
@@ -12,28 +11,7 @@
     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():
+def linear(x, y, new_x):
     """ Linearly interpolates values in new_x based on the values in x and y
 
         Parameters
@@ -45,42 +23,23 @@
         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)
+    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"
-        
-    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
+    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)
 
-def block(x, y, new_x):
-    """ Used when only one element is available in the log.
-    """
+    return new_y
 
-    # 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
+def logarithmic(x, y, new_x):
+    """ Linearly interpolates values in new_x based in the log space of y.
 
         Parameters
         ----------
@@ -99,37 +58,16 @@
     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])
+            _interpolate.loginterp_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)
+        _interpolate.loginterp_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)
+def block_average_above(x, y, new_x):
+    """ Linearly interpolates values in new_x based on the values in x and 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
@@ -139,6 +77,7 @@
         new_x
             1-D array
     """
+    bad_index = None
     x = make_array_safe(x, numpy.float64)
     y = make_array_safe(y, numpy.float64)
     new_x = make_array_safe(new_x, numpy.float64)
@@ -147,10 +86,93 @@
     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])
+            bad_index = _interpolate.block_averave_above_dddd(x, y[i], 
+                                                            new_x, new_y[i])
+            if bad_index is not None:
+                break                                                
     else:
         new_y = numpy.zeros(len(new_x), numpy.float64)
-        _interpolate.loginterp_dddd(x, y, new_x, new_y)
+        bad_index = _interpolate.block_average_above_dddd(x, y, new_x, new_y)
 
+    if bad_index is not None:
+        msg = "block_average_above cannot extrapolate and new_x[%d]=%f "\
+              "is out of the x range (%f, %f)" % \
+              (bad_index, new_x[bad_index], x[0], x[-1])
+        raise ValueError, msg
+              
     return new_y
+
+
+
+
+def test_helper():
+    """ use numpy.allclose to test
+    """
     
+    print "now testing accuracy of interpolation of linear data"
+    
+    x = numpy.arange(10.)
+    y = 2.0*x
+    c = numpy.array([-1.0, 2.3, 10.5])
+    
+    assert numpy.allclose( linear(x, y, c) , [-2.0, 4.6, 21.0] ), "problem in linear"
+    assert numpy.allclose( logarithmic(x, y, c) , [0. , 4.51738774 , 21.47836848] ), \
+                    "problem with logarithmic"
+    assert numpy.allclose( block_average_above(x, y, c) , [0., 2., 4.] ), \
+                    "problem with block_average_above"
+
+def compare_runtimes():
+    from scipy import arange, ones
+    import time
+    
+    # basic linear interp
+    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 '1d interp (sec):', t2 - t1
+    print new_y[:5]
+
+    # basic block_average_above
+    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 '1d block_average_above (sec):', t2 - t1
+    print new_y[:5]
+    
+    # y linear with y having multiple params
+    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 'fast interpolate (sec):', t2 - t1
+    print new_y[:5,:5]
+
+    # scipy with multiple params
+    import scipy
+    N = 3000.
+    x = arange(N)
+    y = ones((100, N)) * arange(N)
+    new_x = arange(N)
+    t1 = time.clock()
+    interp = scipy.interpolate.interp1d(x, y)
+    new_y = interp(new_x)
+    t2 = time.clock()
+    print 'scipy interp1d (sec):', t2 - t1
+    print new_y[:5,:5]
+
+if __name__ == '__main__':
+    compare_runtimes()
+    test_helper()
+
+# Below is stuff from scipy.interpolate and fitpack
\ No newline at end of file

Added: branches/interpolate2/journal.txt
===================================================================
--- branches/interpolate2/journal.txt	2008-07-14 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/journal.txt	2008-07-14 14:28:50 UTC (rev 4542)
@@ -0,0 +1,83 @@
+# Log in making a new Interpolate module
+
+# The goal
+The current interpolation software is wanting.  My goal is to create a module to 
+replace scipy.interpolate.  Currently functionality is mostly in scipy.interpolate
+and enthought.interpolate.  A few other functions are also in NDImage and 
+signalprocessing.  
+
+First off, I just want to have a good 1D interpolation callable class, along
+with a straight-up function.  It should have functionality for linear, logarithmic,
+block and spline interpolation.  There should also be extrapolation high and low.
+
+# starting off
+The core function was derived from the skeleton written by Eric Jones while we 
+were talking.  It was basically a callable class with __init__ and __call__ methods.
+__init__ would set x and y, and also decide which method (linear, block, etc) to
+use.
+
+I started off scavenging the enthought.interpolate.py functions.  I turned the
+functions into other callable classes, and created functions which are wrappers
+around the classes.  The Interpolate1D.__init__ method should call a subroutine
+which will appropriately deal with the interpolation argument.  It could be a class,
+a function, or a string.  The string was difficult, as I included an instruction to
+import the named class/function.  But this caused a namespace problem, as iPython
+wouldn''t re-import it after editing and calling Interpolate1D again.
+
+# stuff incorporated
+Now I've incorporatd the functionality form interpolate.py as callable classes.
+I have several directions to go.
+1) create tests for all of the functionality incorporated so far.  This would be boring
+    but good practice and really must be done.
+2) add handlers for NaN and non-standard data types.  This also must be done.
+    In particular, there's lists, arrays, and naked atomic types (ints, floats, etc).
+    I think everything should be cast to a numpy array.  But should it be a numpy
+    array of float64?  
+3) incorporate the functionality of fitting.py.  This uses splines and calls
+    scipy.interpolate.  It also uses traits, which is taboo.  I should really revamp
+    the whole thing; it should directly call the fortran subroutines, rather than calling
+    them through scipy.interpolate.  Thus, I should really be scavenging scipy but
+    putting a good face on it.
+    
+There's a question of whether interpolate_helper should include all of these functions
+or should be split into a couple modules.  At least preliminarily, I think several modules
+is a good thing because it makes testing easier.
+
+But for the time being, I'll write a test function for the stuff currently in 
+interpolate_helper.
+
+# July 12 2:18 PM
+the block function is strange; wrong, it appears.  The answers depend on other 
+elements in the argument array.  Block doesn't call C code, so it's more suspect.  
+Perhaps I just don't understand the function
+well enough.  Linear works erfectly, logarithmic seems good too, though it's zero 
+for negative values.  block_average above looks good.
+Window_width has a problem with the C call.  So window_width and block need
+more thought.
+
+Block shows same problem in mine and enthought.interpolate, but the non-function
+thing of block_avg_above is unique to me.
+
+Window_Average has a low-level bug.  I'm scrapping it for now, since it's rarely used
+anyway.  Block comes not from interpolate, but from fitting.  Thus I've gotten all
+interpolate.py functionality included.  Yay.
+
+
+# July 12 5 PM
+moving on to fitting.py.  It has a parent DataFit class whose children call functions in interpolate.py.
+That makes for more pure-function stuff, thanks to interpolate.py.  This has objects
+wrapped around funcs, whereas I have funcs wrapped around objects.  Objects around
+funcs has less overhead (I think) and is more pure-functional, so I'll redo things to
+follow that.
+
+Redoing it amounts to copying interpolate.py over everything I've done.  Damn it.
+But the main function 1) needs assert statement, and 2) needlessly defines x
+and y a lot (though that makes one per function)
+
+
+# July 14
+Looking at DataFit in fitting.py I'm debating how to replicate it.  In
+particular, I don't need set_xy and such.  However, those functions
+allow for more flexibility in future code (setting may be more complicated
+for some descendent class) so they're good.  However, I'm preliminarily
+making that function private; x and y must be initialized.
\ No newline at end of file



More information about the Scipy-svn mailing list