[Scipysvn] r4542  branches/interpolate2
scipysvn@scip...
scipysvn@scip...
Mon Jul 14 09:28:52 CDT 2008
Author: fcady
Date: 20080714 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 20080714 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/Interpolate1D.py 20080714 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 doublereferenced. 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 20080714 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/interpolate_helper.py 20080714 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 = 1e10
 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
1D 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 = 1e10
 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 loglinear 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
1D 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 20080714 04:14:10 UTC (rev 4541)
+++ branches/interpolate2/journal.txt 20080714 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 straightup 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 reimport 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 nonstandard 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 nonfunction
+thing of block_avg_above is unique to me.
+
+Window_Average has a lowlevel 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 purefunction 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 purefunctional, 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 Scipysvn
mailing list