[Scipy-svn] r4542 - branches/interpolate2
scipy-svn at scipy.org
scipy-svn at scipy.org
Mon Jul 14 10:28:52 EDT 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