[Scipy-svn] r4547 - branches/Interpolate1D
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Jul 16 18:10:41 EDT 2008
Author: fcady
Date: 2008-07-16 17:10:40 -0500 (Wed, 16 Jul 2008)
New Revision: 4547
Added:
branches/Interpolate1D/Interpolate1D.py
branches/Interpolate1D/Interpolate1D.pyc
branches/Interpolate1D/__init__.py
branches/Interpolate1D/__init__fit.py
branches/Interpolate1D/_interpolate.cpp
branches/Interpolate1D/_interpolate.pyd
branches/Interpolate1D/fitpack.pyf
branches/Interpolate1D/fitpack_wrapper.py
branches/Interpolate1D/fitpack_wrapper.pyc
branches/Interpolate1D/info_fit.py
branches/Interpolate1D/interpolate.h
branches/Interpolate1D/interpolate_wrapper.py
branches/Interpolate1D/interpolate_wrapper.pyc
branches/Interpolate1D/setup.py
Log:
a few files were not added last time
Added: branches/Interpolate1D/Interpolate1D.py
===================================================================
--- branches/Interpolate1D/Interpolate1D.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/Interpolate1D.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,130 @@
+""" A module for intepolation
+"""
+
+# fixme: information strings giving mathematical descriptions of the actions
+# of the functions.
+
+from interpolate_wrapper import linear, logarithmic, block, block_average_above
+from fitpack_wrapper import Spline
+import numpy as np
+from numpy import array, arange, empty, float64, NaN
+
+# fixme: use this to ensure proper type of all inputs and outputs in Interpolate1D
+def make_array_safe(ary, typecode=np.float64):
+ ary = np.atleast_1d(np.asarray(ary, typecode))
+ if not ary.flags['CONTIGUOUS']:
+ ary = ary.copy()
+ return ary
+
+
+class Interpolate1D(object):
+ # see enthought.interpolate
+
+ # fixme: Handle other data types.
+
+ def __init__(self, x, y, k=1, kind='linear', low=None, high=None):
+
+ # fixme: Handle checking if they are the correct size.
+ self._x = make_array_safe(x).copy()
+ self._y = make_array_safe(y).copy()
+ self._xdtype = type(self._x[0])
+ self._ydtype = type(self._y[0])
+
+ assert( len(x) == len(y) , "x and y must be of the same length" )
+ assert( x.ndim == 1 , "x must be one-dimensional" )
+ assert( y.ndim == 1 , "y must be one-dimensional" )
+ # fixme: let y be 2-dimensional. Involves reworking of Interpolate1D.__call__
+ # because Spline enumerates y along the last, rather then first, axis,
+ # while concatenate works along first axis
+
+ self.kind = self._init_interp_method(self._x, self._y, k, kind)
+ self.low = self._init_interp_method(self._x, self._y, k, low)
+ self.high = self._init_interp_method(self._x, self._y, k, high)
+
+ def _init_interp_method(self, x, y, k, interp_arg):
+ from inspect import isclass, isfunction
+
+ if interp_arg in ['linear', 'logarithmic', 'block', 'block_average_above']:
+ func = {'linear':linear, 'logarithmic':logarithmic, 'block':block, \
+ 'block_average_above':block_average_above}[interp_arg]
+ result = lambda new_x : func(self._x, self._y, new_x)
+ elif interp_arg in ['Spline', Spline, 'spline']:
+ result = Spline(self._x, self._y, k=k)
+ elif isfunction(interp_arg):
+ result = interp_arg
+ elif isclass(interp_arg):
+ result = interp_arg(x, y)
+ else:
+ print "warning: defaulting on extrapolation"
+ result = np.vectorize(lambda new_x : interp_arg)
+ return result
+
+ def __call__(self, x):
+
+ x = make_array_safe(x)
+ low_mask = x<self._x[0]
+ high_mask = x>self._x[-1]
+ interp_mask = (~low_mask) & (~high_mask)
+
+ # hack, since getting an error when self.low or self.high gets 0-length array
+ # and they return None or NaN
+ if len(x[low_mask]) == 0: new_low=np.array([])
+ else: new_low = self.low(x[low_mask])
+ if len(x[interp_mask])==0: new_interp=np.array([])
+ else: new_interp = self.kind(x[interp_mask])
+ if len(x[high_mask]) == 0: new_high = np.array([])
+ else: new_high = self.high(x[high_mask])
+
+ result = np.concatenate((new_low, new_interp, new_high))
+
+ return result
+
+# unit testing
+import unittest, time
+class Test(unittest.TestCase):
+
+ def assertAllclose(self, x, y):
+ self.assert_(np.allclose(make_array_safe(x), make_array_safe(y)))
+
+ # fixme: run the test contained in the wrapper modules
+
+ def test_Interp_linearSpl(self):
+ #return
+ print "\n\nTESTING LINEAR (1st ORDER) SPLINE"
+ N = 7
+ x = np.arange(N)
+ y = np.arange(N)
+ T1 = time.clock()
+ interp_func = Interpolate1D(x, y, k=1, kind='Spline', low=None, high=None)
+ T2 = time.clock()
+ print 'time to create linear interp function: ', T2 - T1
+ new_x = np.arange(N)-0.5
+ t1 = time.clock()
+ new_y = interp_func(new_x)
+ t2 = time.clock()
+ print '1d interp (sec):', t2 - t1
+
+ print "new_y: ", new_y
+ self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
+ self.assert_(new_y[0] == None)
+
+ def test_linear(self):
+ print "\n\nTESTING LINEAR INTERPOLATION"
+ N = 7
+ x = arange(N)
+ y = arange(N)
+ new_x = arange(N+1)-0.5
+ T1 = time.clock()
+ interp_func = Interpolate1D(x, y, kind='linear', low=None, high=None)
+ T2 = time.clock()
+ print 'time to create linear interp function: ', T2 - T1
+ t1 = time.clock()
+ new_y = interp_func(new_x)
+ t2 = time.clock()
+ print '1d interp (sec):', t2 - t1
+
+ self.assertAllclose(new_y[1:5], [0.5, 1.5, 2.5, 3.5])
+ self.assert_(new_y[0] == None)
+
+if __name__ == '__main__':
+ unittest.main()
\ No newline at end of file
Added: branches/Interpolate1D/Interpolate1D.pyc
===================================================================
(Binary files differ)
Property changes on: branches/Interpolate1D/Interpolate1D.pyc
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: branches/Interpolate1D/__init__.py
===================================================================
--- branches/Interpolate1D/__init__.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/__init__.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,3 @@
+from interpolate_wrapper import linear, logarithmicm, block_average_above
+
+from Interpolate1D import Interpolate1D
\ No newline at end of file
Added: branches/Interpolate1D/__init__fit.py
===================================================================
--- branches/Interpolate1D/__init__fit.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/__init__fit.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,15 @@
+#
+# interpolate - Interpolation Tools
+#
+
+from info import __doc__
+
+#from interpolate import *
+#from fitpack import *
+
+# New interface to fitpack library:
+from fitpack2 import *
+
+__all__ = filter(lambda s:not s.startswith('_'),dir())
+from numpy.testing import NumpyTest
+test = NumpyTest().test
Added: branches/Interpolate1D/_interpolate.cpp
===================================================================
--- branches/Interpolate1D/_interpolate.cpp 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/_interpolate.cpp 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,236 @@
+#include "Python.h"
+#include <stdlib.h>
+
+#include "interpolate.h"
+#include "numpy/arrayobject.h"
+
+using namespace std;
+
+extern "C" {
+
+static PyObject* linear_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:linear_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ linear((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* loginterp_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ loginterp((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* window_average_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+ double width;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOOd:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y, &width))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ window_average((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0), width);
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyObject* block_average_above_method(PyObject*self, PyObject* args, PyObject* kywds)
+{
+ static char *kwlist[] = {"x","y","new_x","new_y", NULL};
+ PyObject *py_x, *py_y, *py_new_x, *py_new_y;
+ py_x = py_y = py_new_x = py_new_y = NULL;
+ PyObject *arr_x, *arr_y, *arr_new_x, *arr_new_y;
+ arr_x = arr_y = arr_new_x = arr_new_y = NULL;
+
+ if(!PyArg_ParseTupleAndKeywords(args,kywds,"OOOO:loginterp_dddd",kwlist,&py_x, &py_y, &py_new_x, &py_new_y))
+ return NULL;
+ arr_x = PyArray_FROMANY(py_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_x) {
+ PyErr_SetString(PyExc_ValueError, "x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_y = PyArray_FROMANY(py_y, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_y) {
+ PyErr_SetString(PyExc_ValueError, "y must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_x = PyArray_FROMANY(py_new_x, PyArray_DOUBLE, 1, 1, NPY_IN_ARRAY);
+ if (!arr_new_x) {
+ PyErr_SetString(PyExc_ValueError, "new_x must be a 1-D array of floats");
+ goto fail;
+ }
+ arr_new_y = PyArray_FROMANY(py_new_y, PyArray_DOUBLE, 1, 1, NPY_INOUT_ARRAY);
+ if (!arr_new_y) {
+ PyErr_SetString(PyExc_ValueError, "new_y must be a 1-D array of floats");
+ goto fail;
+ }
+
+ block_average_above((double*)PyArray_DATA(arr_x), (double*)PyArray_DATA(arr_y),
+ PyArray_DIM(arr_x,0), (double*)PyArray_DATA(arr_new_x),
+ (double*)PyArray_DATA(arr_new_y), PyArray_DIM(arr_new_x,0));
+
+ Py_DECREF(arr_x);
+ Py_DECREF(arr_y);
+ Py_DECREF(arr_new_x);
+ Py_DECREF(arr_new_y);
+
+ Py_RETURN_NONE;
+
+fail:
+ Py_XDECREF(arr_x);
+ Py_XDECREF(arr_y);
+ Py_XDECREF(arr_new_x);
+ Py_XDECREF(arr_new_y);
+ return NULL;
+}
+
+static PyMethodDef interpolate_methods[] = {
+ {"linear_dddd", (PyCFunction)linear_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"loginterp_dddd", (PyCFunction)loginterp_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"window_average_ddddd", (PyCFunction)window_average_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {"block_average_above_dddd", (PyCFunction)block_average_above_method, METH_VARARGS|METH_KEYWORDS,
+ ""},
+ {NULL, NULL, 0, NULL}
+};
+
+
+PyMODINIT_FUNC init_interpolate(void)
+{
+ PyObject* m;
+ m = Py_InitModule3("_interpolate", interpolate_methods,
+ "A few interpolation routines.\n"
+ );
+ if (m == NULL)
+ return;
+ import_array();
+}
+
+} // extern "C"
Added: branches/Interpolate1D/_interpolate.pyd
===================================================================
(Binary files differ)
Property changes on: branches/Interpolate1D/_interpolate.pyd
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: branches/Interpolate1D/fitpack.pyf
===================================================================
--- branches/Interpolate1D/fitpack.pyf 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/fitpack.pyf 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,479 @@
+! -*- f90 -*-
+! Author: Pearu Peterson <pearu at cens.ioc.ee>
+!
+python module dfitpack ! in
+
+ usercode '''
+
+static double dmax(double* seq,int len) {
+ double val;
+ int i;
+ if (len<1)
+ return -1e308;
+ val = seq[0];
+ for(i=1;i<len;++i)
+ if (seq[i]>val) val = seq[i];
+ return val;
+}
+static double dmin(double* seq,int len) {
+ double val;
+ int i;
+ if (len<1)
+ return 1e308;
+ val = seq[0];
+ for(i=1;i<len;++i)
+ if (seq[i]<val) val = seq[i];
+ return val;
+}
+static double calc_b(double* x,int m,double* tx,int nx) {
+ double val1 = dmin(x,m);
+ double val2 = dmin(tx,nx);
+ if (val2>val1) return val1;
+ val1 = dmax(tx,nx);
+ return val2 - (val1-val2)/nx;
+}
+static double calc_e(double* x,int m,double* tx,int nx) {
+ double val1 = dmax(x,m);
+ double val2 = dmax(tx,nx);
+ if (val2<val1) return val1;
+ val1 = dmin(tx,nx);
+ return val2 + (val2-val1)/nx;
+}
+static int imax(int i1,int i2) {
+ return MAX(i1,i2);
+}
+
+static int calc_surfit_lwrk1(int m, int kx, int ky, int nxest, int nyest) {
+ int u = nxest-kx-1;
+ int v = nyest-ky-1;
+ int km = MAX(kx,ky)+1;
+ int ne = MAX(nxest,nyest);
+ int bx = kx*v+ky+1;
+ int by = ky*u+kx+1;
+ int b1,b2;
+ if (bx<=by) {b1=bx;b2=bx+v-ky;}
+ else {b1=by;b2=by+u-kx;}
+ return u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1;
+}
+static int calc_surfit_lwrk2(int m, int kx, int ky, int nxest, int nyest) {
+ int u = nxest-kx-1;
+ int v = nyest-ky-1;
+ int bx = kx*v+ky+1;
+ int by = ky*u+kx+1;
+ int b2 = (bx<=by?bx+v-ky:by+u-kx);
+ return u*v*(b2+1)+b2;
+}
+
+static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
+ int nxest, int nyest) {
+ int u = MAX(my, nxest);
+ return 4+nxest*(my+2*kx+5)+nyest*(2*ky+5)+mx*(kx+1)+my*(ky+1)+u;
+}
+'''
+
+ interface
+
+ !!!!!!!!!! Univariate spline !!!!!!!!!!!
+
+ subroutine splev(t,n,c,k,x,y,m,ier)
+ ! y = splev(t,c,k,x)
+ real*8 dimension(n),intent(in) :: t
+ integer intent(hide),depend(t) :: n=len(t)
+ real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
+ integer :: k
+ real*8 dimension(m),intent(in) :: x
+ real*8 dimension(m),depend(m),intent(out) :: y
+ integer intent(hide),depend(x) :: m=len(x)
+ integer intent(hide) :: ier
+ end subroutine splev
+
+ subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
+ ! dy = splder(t,c,k,x,[nu])
+ real*8 dimension(n) :: t
+ integer depend(t),intent(hide) :: n=len(t)
+ real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
+ integer :: k
+ integer depend(k),check(0<=nu && nu<=k) :: nu = 1
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),intent(out) :: y
+ integer depend(x),intent(hide) :: m=len(x)
+ real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
+ integer intent(hide) :: ier
+ end subroutine splder
+
+ function splint(t,n,c,k,a,b,wrk)
+ ! iy = splint(t,c,k,a,b)
+ real*8 dimension(n),intent(in) :: t
+ integer intent(hide),depend(t) :: n=len(t)
+ real*8 dimension(n),depend(n),check(len(c)==n) :: c
+ integer intent(in) :: k
+ real*8 intent(in) :: a
+ real*8 intent(in) :: b
+ real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
+ real*8 :: splint
+ end function splint
+
+ subroutine sproot(t,n,c,zero,mest,m,ier)
+ ! zero,m,ier = sproot(t,c[,mest])
+ real*8 dimension(n),intent(in) :: t
+ integer intent(hide),depend(t),check(n>=8) :: n=len(t)
+ real*8 dimension(n),depend(n),check(len(c)==n) :: c
+ real*8 dimension(mest),intent(out),depend(mest) :: zero
+ integer optional,intent(in),depend(n) :: mest=3*(n-7)
+ integer intent(out) :: m
+ integer intent(out) :: ier
+ end subroutine sproot
+
+ subroutine spalde(t,n,c,k,x,d,ier)
+ ! d,ier = spalde(t,c,k,x)
+
+ callprotoargument double*,int*,double*,int*,double*,double*,int*
+ callstatement {int k1=k+1; (*f2py_func)(t,&n,c,&k1,&x,d,&ier); }
+
+ real*8 dimension(n) :: t
+ integer intent(hide),depend(t) :: n=len(t)
+ real*8 dimension(n),depend(n),check(len(c)==n) :: c
+ integer intent(in) :: k
+ real*8 intent(in) :: x
+ real*8 dimension(k+1),intent(out),depend(k) :: d
+ integer intent(out) :: ier
+ end subroutine spalde
+
+ subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
+ ! in curfit.f
+ integer :: iopt
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m) :: w
+ real*8 optional,depend(x),check(xb<=x[0]) :: xb = x[0]
+ real*8 optional,depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
+ integer optional,check(1<=k && k <=5),intent(in) :: k=3
+ real*8 optional,check(s>=0.0) :: s = 0.0
+ integer intent(hide),depend(t) :: nest=len(t)
+ integer intent(out), depend(nest) :: n=nest
+ real*8 dimension(nest),intent(inout) :: t
+ real*8 dimension(n),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk),intent(inout) :: wrk
+ integer intent(hide),depend(wrk) :: lwrk=len(wrk)
+ integer dimension(nest),intent(inout) :: iwrk
+ integer intent(out) :: ier
+ end subroutine curfit
+
+ subroutine percur(iopt,m,x,y,w,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
+ ! in percur.f
+ integer :: iopt
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m) :: w
+ integer optional,check(1<=k && k <=5),intent(in) :: k=3
+ real*8 optional,check(s>=0.0) :: s = 0.0
+ integer intent(hide),depend(t) :: nest=len(t)
+ integer intent(out), depend(nest) :: n=nest
+ real*8 dimension(nest),intent(inout) :: t
+ real*8 dimension(n),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk),intent(inout) :: wrk
+ integer intent(hide),depend(wrk) :: lwrk=len(wrk)
+ integer dimension(nest),intent(inout) :: iwrk
+ integer intent(out) :: ier
+ end subroutine percur
+
+
+ subroutine parcur(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,c,fp,wrk,lwrk,iwrk,ier)
+ ! in parcur.f
+ integer check(iopt>=-1 && iopt <= 1):: iopt
+ integer check(ipar == 1 || ipar == 0) :: ipar
+ integer check(idim > 0 && idim < 11) :: idim
+ integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+ real*8 dimension(m), intent(inout) :: u
+ integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+ real*8 dimension(mx) :: x
+ real*8 dimension(m) :: w
+ real*8 :: ub
+ real*8 :: ue
+ integer optional, check(1<=k && k<=5) :: k=3.0
+ real*8 optional, check(s>=0.0) :: s = 0.0
+ integer intent(hide), depend(t) :: nest=len(t)
+ integer intent(out), depend(nest) :: n=nest
+ real*8 dimension(nest), intent(inout) :: t
+ integer intent(hide), depend(c,nest,idim), check(nc>=idim*nest) :: nc=len(c)
+ real*8 dimension(nc), intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk), intent(inout) :: wrk
+ integer intent(hide),depend(wrk) :: lwrk=len(wrk)
+ integer dimension(nest), intent(inout) :: iwrk
+ integer intent(out) :: ier
+ end subroutine parcur
+
+
+ subroutine fpcurf0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
+ ! fpcurf0(x,y,k,[w,xb,xe,s,nest])
+
+ fortranname fpcurf
+ callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
+ callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
+
+ integer intent(hide) :: iopt = 0
+ real*8 dimension(m),intent(in,out) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
+ real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
+ integer check(1<=k && k<=5),intent(in,out) :: k
+ real*8 check(s>=0.0),depend(m),intent(in,out) :: s = m
+ integer intent(in),depend(m,s,k,k1),check(nest>=2*k1) :: nest = (s==0.0?m+k+1:MAX(m/2,2*k1))
+ real*8 intent(hide) :: tol = 0.001
+ integer intent(hide) :: maxit = 20
+ integer intent(hide),depend(k) :: k1=k+1
+ integer intent(hide),depend(k) :: k2=k+2
+ integer intent(out) :: n
+ real*8 dimension(nest),intent(out),depend(nest) :: t
+ real*8 dimension(nest),depend(nest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
+ real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
+ integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
+ integer intent(out) :: ier
+ end subroutine fpcurf0
+
+ subroutine fpcurf1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
+ ! fpcurf1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
+
+ fortranname fpcurf
+ callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
+ callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
+
+ integer intent(hide) :: iopt = 1
+ real*8 dimension(m),intent(in,out,overwrite) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out,overwrite) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out,overwrite) :: w
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 intent(in,out) :: xb
+ real*8 intent(in,out) :: xe
+ integer check(1<=k && k<=5),intent(in,out) :: k
+ real*8 check(s>=0.0),intent(in,out) :: s
+ integer intent(hide),depend(t) :: nest = len(t)
+ real*8 intent(hide) :: tol = 0.001
+ integer intent(hide) :: maxit = 20
+ integer intent(hide),depend(k) :: k1=k+1
+ integer intent(hide),depend(k) :: k2=k+2
+ integer intent(in,out) :: n
+ real*8 dimension(nest),intent(in,out,overwrite) :: t
+ real*8 dimension(nest),depend(nest),check(len(c)==nest),intent(in,out,overwrite) :: c
+ real*8 intent(in,out) :: fp
+ real*8 dimension(nest),depend(nest),check(len(fpint)==nest),intent(in,out,cache,overwrite) :: fpint
+ real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
+ integer dimension(nest),depend(nest),check(len(nrdata)==nest),intent(in,out,cache,overwrite) :: nrdata
+ integer intent(in,out) :: ier
+ end subroutine fpcurf1
+
+ subroutine fpcurfm1(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,c,fp,fpint,wrk,nrdata,ier)
+ ! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
+ ! fpcurfm1(x,y,k,t,[w,xb,xe])
+
+ fortranname fpcurf
+ callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
+ callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
+
+ integer intent(hide) :: iopt = -1
+ real*8 dimension(m),intent(in,out) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
+ real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
+ integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
+ real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
+ real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
+ integer check(1<=k && k<=5),intent(in,out) :: k
+ real*8 intent(out) :: s = -1
+ integer intent(hide),depend(n) :: nest = n
+ real*8 intent(hide) :: tol = 0.001
+ integer intent(hide) :: maxit = 20
+ integer intent(hide),depend(k) :: k1=k+1
+ integer intent(hide),depend(k) :: k2=k+2
+ integer intent(out),depend(t) :: n = len(t)
+ real*8 dimension(n),intent(in,out,overwrite) :: t
+ real*8 dimension(nest),depend(nest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(nest),depend(nest),intent(out,cache) :: fpint
+ real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
+ integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
+ integer intent(out) :: ier
+ end subroutine fpcurfm1
+
+ !!!!!!!!!! Bivariate spline !!!!!!!!!!!
+
+ subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,iwrk,kwrk,ier)
+ ! z,ier = bispev(tx,ty,c,kx,ky,x,y)
+ real*8 dimension(nx),intent(in) :: tx
+ integer intent(hide),depend(tx) :: nx=len(tx)
+ real*8 dimension(ny),intent(in) :: ty
+ integer intent(hide),depend(ty) :: ny=len(ty)
+ real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
+ check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
+ integer :: kx
+ integer :: ky
+ real*8 intent(in),dimension(mx) :: x
+ integer intent(hide),depend(x) :: mx=len(x)
+ real*8 intent(in),dimension(my) :: y
+ integer intent(hide),depend(y) :: my=len(y)
+ real*8 dimension(mx,my),depend(mx,my),intent(out,c) :: z
+ real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
+ integer intent(hide),depend(mx,kx,my,ky) :: lwrk=mx*(kx+1)+my*(ky+1)
+ integer dimension(kwrk),depend(kwrk),intent(hide,cache) :: iwrk
+ integer intent(hide),depend(mx,my) :: kwrk=mx+my
+ integer intent(out) :: ier
+ end subroutine bispev
+
+ subroutine surfit_smth(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
+ nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
+ iwrk,kwrk,ier)
+ ! nx,tx,ny,ty,c,fp,ier = surfit_smth(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
+
+ fortranname surfit
+
+ integer intent(hide) :: iopt=0
+ integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
+ :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(z)==m) :: z
+ real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
+ real*8 optional,depend(x,m) :: xb=dmin(x,m)
+ real*8 optional,depend(x,m) :: xe=dmax(x,m)
+ real*8 optional,depend(y,m) :: yb=dmin(y,m)
+ real*8 optional,depend(y,m) :: ye=dmax(y,m)
+ integer check(1<=kx && kx<=5) :: kx = 3
+ integer check(1<=ky && ky<=5) :: ky = 3
+ real*8 optional,check(0.0<=s) :: s = m
+ integer optional,depend(kx,m),check(nxest>=2*(kx+1)) &
+ :: nxest = imax(kx+1+sqrt(m/2),2*(kx+1))
+ integer optional,depend(ky,m),check(nyest>=2*(ky+1)) &
+ :: nyest = imax(ky+1+sqrt(m/2),2*(ky+1))
+ integer intent(hide),depend(nxest,nyest) :: nmax=MAX(nxest,nyest)
+ real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
+ integer intent(out) :: nx
+ real*8 dimension(nmax),intent(out),depend(nmax) :: tx
+ integer intent(out) :: ny
+ real*8 dimension(nmax),intent(out),depend(nmax) :: ty
+ real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
+ depend(kx,ky,nxest,nyest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk1),intent(cache,out),depend(lwrk1) :: wrk1
+ integer intent(hide),depend(m,kx,ky,nxest,nyest) &
+ :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
+ real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
+ integer optional,intent(in),depend(kx,ky,nxest,nyest) &
+ :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
+ integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
+ integer intent(hide),depend(m,nxest,nyest,kx,ky) &
+ :: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
+ integer intent(out) :: ier
+ end subroutine surfit_smth
+
+ subroutine surfit_lsq(iopt,m,x,y,z,w,xb,xe,yb,ye,kx,ky,s,nxest,nyest,&
+ nmax,eps,nx,tx,ny,ty,c,fp,wrk1,lwrk1,wrk2,lwrk2,&
+ iwrk,kwrk,ier)
+ ! tx,ty,c,fp,ier = surfit_lsq(x,y,z,tx,ty,[w,xb,xe,yb,ye,kx,ky,eps,lwrk2])
+
+ fortranname surfit
+
+ integer intent(hide) :: iopt=-1
+ integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
+ :: m=len(x)
+ real*8 dimension(m) :: x
+ real*8 dimension(m),depend(m),check(len(y)==m) :: y
+ real*8 dimension(m),depend(m),check(len(z)==m) :: z
+ real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
+ real*8 optional,depend(x,tx,m,nx) :: xb=calc_b(x,m,tx,nx)
+ real*8 optional,depend(x,tx,m,nx) :: xe=calc_e(x,m,tx,nx)
+ real*8 optional,depend(y,ty,m,ny) :: yb=calc_b(y,m,ty,ny)
+ real*8 optional,depend(y,ty,m,ny) :: ye=calc_e(y,m,ty,ny)
+ integer check(1<=kx && kx<=5) :: kx = 3
+ integer check(1<=ky && ky<=5) :: ky = 3
+ real*8 intent(hide) :: s = 0.0
+ integer intent(hide),depend(nx) :: nxest = nx
+ integer intent(hide),depend(ny) :: nyest = ny
+ integer intent(hide),depend(nx,ny) :: nmax=MAX(nx,ny)
+ real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
+ integer intent(hide),depend(tx,kx),check(2*kx+2<=nx) :: nx = len(tx)
+ real*8 dimension(nx),intent(in,out,overwrite) :: tx
+ integer intent(hide),depend(ty,ky),check(2*ky+2<=ny) :: ny = len(ty)
+ real*8 dimension(ny),intent(in,out,overwrite) :: ty
+ real*8 dimension((nx-kx-1)*(ny-ky-1)),depend(kx,ky,nx,ny),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk1),intent(cache,hide),depend(lwrk1) :: wrk1
+ integer intent(hide),depend(m,kx,ky,nxest,nyest) &
+ :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
+ real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
+ integer optional,intent(in),depend(kx,ky,nxest,nyest) &
+ :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
+ integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
+ integer intent(hide),depend(m,nx,ny,kx,ky) &
+ :: kwrk=m+(nx-2*kx-1)*(ny-2*ky-1)
+ integer intent(out) :: ier
+ end subroutine surfit_lsq
+
+ subroutine regrid_smth(iopt,mx,x,my,y,z,xb,xe,yb,ye,kx,ky,s,&
+ nxest,nyest,nx,tx,ny,ty,c,fp,wrk,lwrk,iwrk,kwrk,ier)
+ ! nx,tx,ny,ty,c,fp,ier = regrid_smth(x,y,z,[xb,xe,yb,ye,kx,ky,s])
+
+ fortranname regrid
+
+ integer intent(hide) :: iopt=0
+ integer intent(hide),depend(x,kx),check(mx>kx) :: mx=len(x)
+ real*8 dimension(mx) :: x
+ integer intent(hide),depend(y,ky),check(my>ky) :: my=len(y)
+ real*8 dimension(my) :: y
+ real*8 dimension(mx*my),depend(mx,my),check(len(z)==mx*my) :: z
+ real*8 optional,depend(x,mx) :: xb=dmin(x,mx)
+ real*8 optional,depend(x,mx) :: xe=dmax(x,mx)
+ real*8 optional,depend(y,my) :: yb=dmin(y,my)
+ real*8 optional,depend(y,my) :: ye=dmax(y,my)
+ integer optional,check(1<=kx && kx<=5) :: kx = 3
+ integer optional,check(1<=ky && ky<=5) :: ky = 3
+ real*8 optional,check(0.0<=s) :: s = 0.0
+ integer intent(hide),depend(kx,mx),check(nxest>=2*(kx+1)) &
+ :: nxest = mx+kx+1
+ integer intent(hide),depend(ky,my),check(nyest>=2*(ky+1)) &
+ :: nyest = my+ky+1
+ integer intent(out) :: nx
+ real*8 dimension(nxest),intent(out),depend(nxest) :: tx
+ integer intent(out) :: ny
+ real*8 dimension(nyest),intent(out),depend(nyest) :: ty
+ real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
+ depend(kx,ky,nxest,nyest),intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk),intent(cache,hide),depend(lwrk) :: wrk
+ integer intent(hide),depend(mx,my,kx,ky,nxest,nyest) &
+ :: lwrk=calc_regrid_lwrk(mx,my,kx,ky,nxest,nyest)
+ integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
+ integer intent(hide),depend(mx,my,nxest,nyest) &
+ :: kwrk=3+mx+my+nxest+nyest
+ integer intent(out) :: ier
+ end subroutine regrid_smth
+
+ function dblint(tx,nx,ty,ny,c,kx,ky,xb,xe,yb,ye,wrk)
+ ! iy = dblint(tx,ty,c,kx,ky,xb,xe,yb,ye)
+ real*8 dimension(nx),intent(in) :: tx
+ integer intent(hide),depend(tx) :: nx=len(tx)
+ real*8 dimension(ny),intent(in) :: ty
+ integer intent(hide),depend(ty) :: ny=len(ty)
+ real*8 intent(in),dimension((nx-kx-1)*(ny-ky-1)),depend(nx,ny,kx,ky),&
+ check(len(c)==(nx-kx-1)*(ny-ky-1)):: c
+ integer :: kx
+ integer :: ky
+ real*8 intent(in) :: xb
+ real*8 intent(in) :: xe
+ real*8 intent(in) :: yb
+ real*8 intent(in) :: ye
+ real*8 dimension(nx+ny-kx-ky-2),depend(nx,ny,kx,ky),intent(cache,hide) :: wrk
+ real*8 :: dblint
+ end function dblint
+ end interface
+end python module dfitpack
+
Added: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/fitpack_wrapper.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,204 @@
+""" fit_helper.py
+
+mimics the functionality of enthought.interpolate that is
+contained in the module fitting.py
+
+"""
+
+import numpy as np
+
+import dfitpack #fixme: rename module fitpack_wrapper
+
+
+class Spline(object):
+ """ Univariate spline s(x) of degree k on the interval
+ [xb,xe] calculated from a given set of data points
+ (x,y).
+
+ Can include least-squares fitting.
+
+ See also:
+
+ splrep, splev, sproot, spint, spalde - an older wrapping of FITPACK
+ BivariateSpline - a similar class for bivariate spline interpolation
+ """
+
+ def __init__(self, x, y, w=None, bbox = [None]*2, k=3, s=None):
+ """
+ Input:
+ x,y - 1-d sequences of data points (x must be
+ in strictly ascending order)
+
+ Optional input:
+ w - positive 1-d sequence of weights
+ bbox - 2-sequence specifying the boundary of
+ the approximation interval.
+ By default, bbox=[x[0],x[-1]]
+ k=3 - degree of the univariate spline.
+ s - positive smoothing factor defined for
+ estimation condition:
+ sum((w[i]*(y[i]-s(x[i])))**2,axis=0) <= s
+ Default s=len(w) which should be a good value
+ if 1/w[i] is an estimate of the standard
+ deviation of y[i].
+ """
+ #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
+ data = dfitpack.fpcurf0(x,y,k,w=w,
+ xb=bbox[0],xe=bbox[1],s=s)
+ if data[-1]==1:
+ # nest too small, setting to maximum bound
+ data = self._reset_nest(data)
+ self._data = data
+ self._reset_class()
+
+ def _reset_class(self):
+ data = self._data
+ n,t,c,k,ier = data[7],data[8],data[9],data[5],data[-1]
+ self._eval_args = t[:n],c[:n],k
+ if ier==0:
+ # the spline returned has a residual sum of squares fp
+ # such that abs(fp-s)/s <= tol with tol a relative
+ # tolerance set to 0.001 by the program
+ pass
+ elif ier==-1:
+ # the spline returned is an interpolating spline
+ #self.__class__ = InterpolatedUnivariateSpline
+ pass
+ elif ier==-2:
+ # the spline returned is the weighted least-squares
+ # polynomial of degree k. In this extreme case fp gives
+ # the upper bound fp0 for the smoothing factor s.
+ #self.__class__ = LSQUnivariateSpline
+ pass
+ else:
+ # error
+ #if ier==1:
+ # self.__class__ = LSQUnivariateSpline
+ #message = _curfit_messages.get(ier,'ier=%s' % (ier))
+ #warnings.warn(message)
+ pass
+
+ def _reset_nest(self, data, nest=None):
+ n = data[10]
+ if nest is None:
+ k,m = data[5],len(data[0])
+ nest = m+k+1 # this is the maximum bound for nest
+ else:
+ assert n<=nest,"nest can only be increased"
+ t,c,fpint,nrdata = data[8].copy(),data[9].copy(),\
+ data[11].copy(),data[12].copy()
+ t.resize(nest)
+ c.resize(nest)
+ fpint.resize(nest)
+ nrdata.resize(nest)
+ args = data[:8] + (t,c,n,fpint,nrdata,data[13])
+ data = dfitpack.fpcurf1(*args)
+ return data
+
+ def set_smoothing_factor(self, s):
+ """ Continue spline computation with the given smoothing
+ factor s and with the knots found at the last call.
+
+ """
+ data = self._data
+ if data[6]==-1:
+ warnings.warn('smoothing factor unchanged for'
+ 'LSQ spline with fixed knots')
+ return
+ args = data[:6] + (s,) + data[7:]
+ data = dfitpack.fpcurf1(*args)
+ if data[-1]==1:
+ # nest too small, setting to maximum bound
+ data = self._reset_nest(data)
+ self._data = data
+ self._reset_class()
+
+ def __call__(self, x, nu=None):
+ """ Evaluate spline (or its nu-th derivative) at positions x.
+ Note: x can be unordered but the evaluation is more efficient
+ if x is (partially) ordered.
+
+ """
+ print 'length of x: ', len(x)
+ if len(x) == 0: return np.array([]) #hack to cope with shape (0,)
+ if nu is None:
+ return dfitpack.splev(*(self._eval_args+(x,)))
+ return dfitpack.splder(nu=nu,*(self._eval_args+(x,)))
+
+ def get_knots(self):
+ """ Return the positions of (boundary and interior)
+ knots of the spline.
+ """
+ data = self._data
+ k,n = data[5],data[7]
+ return data[8][k:n-k]
+
+ def get_coeffs(self):
+ """Return spline coefficients."""
+ data = self._data
+ k,n = data[5],data[7]
+ return data[9][:n-k-1]
+
+ def get_residual(self):
+ """Return weighted sum of squared residuals of the spline
+ approximation: sum ((w[i]*(y[i]-s(x[i])))**2,axis=0)
+
+ """
+ return self._data[10]
+
+ def integral(self, a, b):
+ """ Return definite integral of the spline between two
+ given points.
+ """
+ return dfitpack.splint(*(self._eval_args+(a,b)))
+
+ def derivatives(self, x):
+ """ Return all derivatives of the spline at the point x."""
+ d,ier = dfitpack.spalde(*(self._eval_args+(x,)))
+ assert ier==0,`ier`
+ return d
+
+ def roots(self):
+ """ Return the zeros of the spline.
+
+ Restriction: only cubic splines are supported by fitpack.
+ """
+ k = self._data[5]
+ if k==3:
+ z,m,ier = dfitpack.sproot(*self._eval_args[:2])
+ assert ier==0,`ier`
+ return z[:m]
+ raise NotImplementedError,\
+ 'finding roots unsupported for non-cubic splines'
+
+
+
+# testing
+import unittest
+import time
+from numpy import arange, allclose, ones
+
+class Test(unittest.TestCase):
+
+ def assertAllclose(self, x, y):
+ self.assert_(np.allclose(x, y))
+
+ def test_linearSpl(self):
+ N = 3000.
+ x = np.arange(N)
+ y = np.arange(N)
+ T1 = time.clock()
+ interp_func = Spline(x, y, k=1)
+ T2 = time.clock()
+ print 'time to create linear interp function: ', T2 - T1
+ new_x = np.arange(N)+0.5
+ t1 = time.clock()
+ new_y = interp_func(new_x)
+ t2 = time.clock()
+ print '1d interp (sec):', t2 - t1
+ self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
+
+if __name__ == '__main__':
+ unittest.main()
+
+
\ No newline at end of file
Added: branches/Interpolate1D/fitpack_wrapper.pyc
===================================================================
(Binary files differ)
Property changes on: branches/Interpolate1D/fitpack_wrapper.pyc
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: branches/Interpolate1D/info_fit.py
===================================================================
--- branches/Interpolate1D/info_fit.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/info_fit.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,12 @@
+"""
+Interpolation Tools
+===================
+
+Primary interpolation agent:
+
+ Interpolate1D -- callable class, initialized with x and y to interpolate
+ from and indicator of how to deal with
+ extrapolation
+"""
+
+postpone_import = 1
Added: branches/Interpolate1D/interpolate.h
===================================================================
--- branches/Interpolate1D/interpolate.h 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/interpolate.h 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,206 @@
+#include <time.h>
+#include <math.h>
+#include <iostream>
+#include <algorithm>
+
+template <class T>
+void linear(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ for (int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ int index;
+ if (new_x <= x_vec[0])
+ index = 0;
+ else if (new_x >=x_vec[len-1])
+ index = len-2;
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ index = which - x_vec-1;
+ }
+
+ if(new_x == x_vec[index])
+ {
+ // exact value
+ new_y_vec[i] = y_vec[index];
+ }
+ else
+ {
+ //interpolate
+ double x_lo = x_vec[index];
+ double x_hi = x_vec[index+1];
+ double y_lo = y_vec[index];
+ double y_hi = y_vec[index+1];
+ double slope = (y_hi-y_lo)/(x_hi-x_lo);
+ new_y_vec[i] = slope * (new_x-x_lo) + y_lo;
+ }
+ }
+}
+
+template <class T>
+void loginterp(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ for (int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ int index;
+ if (new_x <= x_vec[0])
+ index = 0;
+ else if (new_x >=x_vec[len-1])
+ index = len-2;
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ index = which - x_vec-1;
+ }
+
+ if(new_x == x_vec[index])
+ {
+ // exact value
+ new_y_vec[i] = y_vec[index];
+ }
+ else
+ {
+ //interpolate
+ double x_lo = x_vec[index];
+ double x_hi = x_vec[index+1];
+ double y_lo = log10(y_vec[index]);
+ double y_hi = log10(y_vec[index+1]);
+ double slope = (y_hi-y_lo)/(x_hi-x_lo);
+ new_y_vec[i] = pow(10.0, (slope * (new_x-x_lo) + y_lo));
+ }
+ }
+}
+
+template <class T>
+int block_average_above(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len)
+{
+ int bad_index = -1;
+ int start_index = 0;
+ T last_y = 0.0;
+ T thickness = 0.0;
+
+ for(int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ if ((new_x < x_vec[0]) || (new_x > x_vec[len-1]))
+ {
+ bad_index = i;
+ break;
+ }
+ else if (new_x == x_vec[0])
+ {
+ // for the first sample, just return the cooresponding y value
+ new_y_vec[i] = y_vec[0];
+ }
+ else
+ {
+ T* which = std::lower_bound(x_vec, x_vec+len, new_x);
+ int index = which - x_vec-1;
+
+ // calculate weighted average
+
+ // Start off with "residue" from last interval in case last x
+ // was between to samples.
+ T weighted_y_sum = last_y * thickness;
+ T thickness_sum = thickness;
+ for(int j=start_index; j<=index; j++)
+ {
+ T next_x;
+ if (x_vec[j+1] < new_x)
+ thickness = x_vec[j+1] - x_vec[j];
+ else
+ thickness = new_x -x_vec[j];
+ weighted_y_sum += y_vec[j] * thickness;
+ thickness_sum += thickness;
+ }
+ new_y_vec[i] = weighted_y_sum/thickness_sum;
+
+ // Store the thickness between the x value and the next sample
+ // to add to the next weighted average.
+ last_y = y_vec[index];
+ thickness = x_vec[index+1] - new_x;
+
+ // start next weighted average at next sample
+ start_index =index+1;
+ }
+ }
+ return bad_index;
+}
+
+template <class T>
+int window_average(T* x_vec, T* y_vec, int len,
+ T* new_x_vec, T* new_y_vec, int new_len,
+ T width)
+{
+ for(int i=0;i<new_len;i++)
+ {
+ T new_x = new_x_vec[i];
+ T bottom = new_x - width/2;
+ T top = new_x + width/2;
+
+ T* which = std::lower_bound(x_vec, x_vec+len, bottom);
+ int bottom_index = which - x_vec;
+ if (bottom_index < 0)
+ {
+ //bottom = x_vec[0];
+ bottom_index = 0;
+ }
+
+ which = std::lower_bound(x_vec, x_vec+len, top);
+ int top_index = which - x_vec;
+ if (top_index >= len)
+ {
+ //top = x_vec[len-1];
+ top_index = len-1;
+ }
+ //std::cout << std::endl;
+ //std::cout << bottom_index << " " << top_index << std::endl;
+ //std::cout << bottom << " " << top << std::endl;
+ // calculate weighted average
+ T thickness =0.0;
+ T thickness_sum =0.0;
+ T weighted_y_sum =0.0;
+ for(int j=bottom_index; j < top_index; j++)
+ {
+ thickness = x_vec[j+1] - bottom;
+ weighted_y_sum += y_vec[j] * thickness;
+ thickness_sum += thickness;
+ bottom = x_vec[j+1];
+ /*
+ std::cout << "iter: " << j - bottom_index << " " <<
+ "index: " << j << " " <<
+ "bottom: " << bottom << " " <<
+ "x+1: " << x_vec[j+1] << " " <<
+ "x: " << x_vec[j] << " " <<
+ "y: " << y_vec[j] << " " <<
+ "weighted_sum: " << weighted_y_sum <<
+ "thickness: " << thickness << " " <<
+ "thickness_sum: " << thickness_sum << std::endl;
+ */
+ //std::cout << x_vec[j] << " ";
+ //std::cout << thickness << " ";
+ }
+
+ // last element
+ thickness = top - bottom;
+ weighted_y_sum += y_vec[top_index] * thickness;
+ thickness_sum += thickness;
+ /*
+ std::cout << "iter: last" << " " <<
+ "index: " << top_index << " " <<
+ "x: " << x_vec[top_index] << " " <<
+ "y: " << y_vec[top_index] << " " <<
+ "weighted_sum: " << weighted_y_sum <<
+ "thickness: " << thickness << " " <<
+ "thickness_sum: " << thickness_sum << std::endl;
+ */
+ //std::cout << x_vec[top_index] << " " << thickness_sum << std::endl;
+ new_y_vec[i] = weighted_y_sum/thickness_sum;
+ }
+ return -1;
+}
Added: branches/Interpolate1D/interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/interpolate_wrapper.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/interpolate_wrapper.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,259 @@
+""" helper_funcs.py
+"""
+
+import numpy
+import sys; sys.path.append('C:\home\python\branches\interpolate2\entinterp')
+import _interpolate
+
+def make_array_safe(ary, typecode = numpy.float64):
+ ary = numpy.atleast_1d(numpy.asarray(ary, typecode))
+ if not ary.flags['CONTIGUOUS']:
+ ary = ary.copy()
+ return ary
+
+
+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
+
+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
+
+def block_average_above(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
+ """
+ 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)
+
+ 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)):
+ 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)
+ 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 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(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))
+ new_y = numpy.take(y, indices, axis=-1)
+ 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]
+
+
+# Below is stuff from scipy.interpolate and fitpack
+
+# Unit Test
+import unittest
+import time
+from numpy import arange, allclose, ones, NaN, isnan
+class Test(unittest.TestCase):
+
+ #def assertAllclose(self, x, y, rtol = 1.0e-5):
+ # self.assert_(numpy.allclose(x, y, rtol = rtol))
+
+ def assertAllclose(self, x, y, rtol=1.0e-5):
+ for i, xi in enumerate(x):
+ self.assert_(allclose(xi, y[i], rtol) or (isnan(xi) and isnan(y[i])))
+
+ def test_linear(self):
+ 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
+
+ self.assertAllclose(new_y[:5], [0.5, 1.5, 2.5, 3.5, 4.5])
+
+ def test_block_average_above(self):
+ 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
+ self.assertAllclose(new_y[:5], [0.0, 0.5, 2.5, 4.5, 6.5])
+
+ def test_linear2(self):
+ 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
+ self.assertAllclose(new_y[:5,:5],
+ [[ 0.5, 1.5, 2.5, 3.5, 4.5],
+ [ 0.5, 1.5, 2.5, 3.5, 4.5],
+ [ 0.5, 1.5, 2.5, 3.5, 4.5],
+ [ 0.5, 1.5, 2.5, 3.5, 4.5],
+ [ 0.5, 1.5, 2.5, 3.5, 4.5]])
+
+ def test_logarithmic(self):
+ N = 3000.
+ x = arange(N)
+ y = arange(N)
+ new_x = arange(N)+0.5
+ t1 = time.clock()
+ new_y = logarithmic(x, y, new_x)
+ t2 = time.clock()
+ print 'logarithmic interp (sec):', t2 - t1
+ correct_y = [numpy.NaN, 1.41421356, 2.44948974, 3.46410162, 4.47213595]
+ self.assertAllclose(new_y[:5], correct_y)
+ print "logo"
+
+if __name__ == '__main__':
+ unittest.main()
+
\ No newline at end of file
Added: branches/Interpolate1D/interpolate_wrapper.pyc
===================================================================
(Binary files differ)
Property changes on: branches/Interpolate1D/interpolate_wrapper.pyc
___________________________________________________________________
Name: svn:mime-type
+ application/octet-stream
Added: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py 2008-07-16 22:08:51 UTC (rev 4546)
+++ branches/Interpolate1D/setup.py 2008-07-16 22:10:40 UTC (rev 4547)
@@ -0,0 +1,24 @@
+#!/usr/bin/env python
+
+from os.path import join
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.misc_util import Configuration
+
+ config = Configuration('', parent_package, top_path)
+
+ config.add_extension('_interpolate',
+ ['_interpolate.cpp'],
+ include_dirs = ['.'],
+ depends = ['interpolate.h'])
+
+ config.add_extension('dfitpack',
+ sources=['fitpack.pyf'],
+ libraries=['fitpack'],
+ )
+
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(**configuration(top_path='').todict() )
\ No newline at end of file
More information about the Scipy-svn
mailing list