[Scipy-svn] r4547 - branches/Interpolate1D

scipy-svn@scip... scipy-svn@scip...
Wed Jul 16 17:10:41 CDT 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@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