[Scipy-svn] r2634 - trunk/Lib/interpolate
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Jan 30 07:16:36 EST 2007
Author: jtravs
Date: 2007-01-30 06:16:32 -0600 (Tue, 30 Jan 2007)
New Revision: 2634
Modified:
trunk/Lib/interpolate/__fitpack.h
trunk/Lib/interpolate/_fitpackmodule.c
trunk/Lib/interpolate/fitpack.py
Log:
Added patch for 'insert' funtion from dierckx fitpack for Zachary Pincus.
Modified: trunk/Lib/interpolate/__fitpack.h
===================================================================
--- trunk/Lib/interpolate/__fitpack.h 2007-01-30 04:08:29 UTC (rev 2633)
+++ trunk/Lib/interpolate/__fitpack.h 2007-01-30 12:16:32 UTC (rev 2634)
@@ -16,6 +16,7 @@
{"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
{"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
{"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
+ {"_insert", fitpack_insert, METH_VARARGS, doc_insert},
*/
/* link libraries: (one item per line)
ddierckx
@@ -36,6 +37,7 @@
#define SURFIT surfit
#define BISPEV bispev
#define PARDER parder
+#define INSERT insert
#else
#define CURFIT curfit_
#define PERCUR percur_
@@ -49,6 +51,7 @@
#define SURFIT surfit_
#define BISPEV bispev_
#define PARDER parder_
+#define INSERT insert_
#endif
void CURFIT(int*,int*,double*,double*,double*,double*,double*,int*,double*,int*,int*,double*,double*,double*,double*,int*,int*,int*);
void PERCUR(int*,int*,double*,double*,double*,int*,double*,int*,int*,double*,double*,double*,double*,int*,int*,int*);
@@ -62,6 +65,7 @@
void SURFIT(int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*,double*,int*,int*,int*,double*,int*,double*,int*,double*,double*,double*,double*,int*,double*,int*,int*,int*,int*);
void BISPEV(double*,int*,double*,int*,double*,int*,int*,double*,int*,double*,int*,double*,double*,int*,int*,int*,int*);
void PARDER(double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,int*,double*,int*,double*,double*,int*,int*,int*,int*);
+void INSERT(int*,double*,int*,double*,int*,double*,double*,int*,double*,int*,int*);
/* Note that curev, cualde need no interface. */
@@ -545,6 +549,43 @@
return NULL;
}
+static char doc_insert[] = " [tt,cc,ier] = _insert(iopt,t,c,k,x,m)";
+static PyObject *fitpack_insert(PyObject *dummy, PyObject*args) {
+ int iopt, n, nn, k, nest, ier, m;
+ double x;
+ double *t, *c, *tt, *cc;
+ PyArrayObject *ap_t = NULL, *ap_c = NULL, *ap_tt = NULL, *ap_cc = NULL;
+ PyObject *t_py = NULL, *c_py = NULL;
+ if (!PyArg_ParseTuple(args, "iOOidi",&iopt,&t_py,&c_py,&k, &x, &m)) return NULL;
+ ap_t = (PyArrayObject *)PyArray_ContiguousFromObject(t_py, PyArray_DOUBLE, 0, 1);
+ ap_c = (PyArrayObject *)PyArray_ContiguousFromObject(c_py, PyArray_DOUBLE, 0, 1);
+ if (ap_t == NULL || ap_c == NULL) goto fail;
+ t = (double *) ap_t->data;
+ c = (double *) ap_c->data;
+ n = ap_t->dimensions[0];
+ nest = n + m;
+ ap_tt = (PyArrayObject *)PyArray_FromDims(1,&nest,PyArray_DOUBLE);
+ ap_cc = (PyArrayObject *)PyArray_FromDims(1,&nest,PyArray_DOUBLE);
+ if (ap_tt == NULL || ap_cc == NULL) goto fail;
+ tt = (double *) ap_tt->data;
+ cc = (double *) ap_cc->data;
+ for ( ; n < nest; n++) {
+ INSERT(&iopt, t, &n, c, &k, &x, tt, &nn, cc, &nest, &ier);
+ if (ier) break;
+ t = tt;
+ c = cc;
+ }
+ Py_DECREF(ap_c);
+ Py_DECREF(ap_t);
+ PyObject* ret = Py_BuildValue("NNi",PyArray_Return(ap_tt),PyArray_Return(ap_cc),ier);
+ return ret;
+
+ fail:
+ Py_XDECREF(ap_c);
+ Py_XDECREF(ap_t);
+ return NULL;
+ }
+
Modified: trunk/Lib/interpolate/_fitpackmodule.c
===================================================================
--- trunk/Lib/interpolate/_fitpackmodule.c 2007-01-30 04:08:29 UTC (rev 2633)
+++ trunk/Lib/interpolate/_fitpackmodule.c 2007-01-30 12:16:32 UTC (rev 2634)
@@ -14,6 +14,7 @@
{"_parcur", fitpack_parcur, METH_VARARGS, doc_parcur},
{"_surfit", fitpack_surfit, METH_VARARGS, doc_surfit},
{"_bispev", fitpack_bispev, METH_VARARGS, doc_bispev},
+{"_insert", fitpack_insert, METH_VARARGS, doc_insert},
{NULL, NULL, 0, NULL}
};
PyMODINIT_FUNC init_fitpack(void) {
Modified: trunk/Lib/interpolate/fitpack.py
===================================================================
--- trunk/Lib/interpolate/fitpack.py 2007-01-30 04:08:29 UTC (rev 2633)
+++ trunk/Lib/interpolate/fitpack.py 2007-01-30 12:16:32 UTC (rev 2634)
@@ -29,7 +29,7 @@
"""
__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
- 'bisplrep', 'bisplev']
+ 'bisplrep', 'bisplev', 'insert']
__version__ = "$Revision$"[10:-1]
import _fitpack
from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
@@ -751,7 +751,50 @@
if len(z[0])>1: return z[0]
return z[0][0]
+def insert(x,tck,m=1,per=0):
+ """Insert knots into a B-spline.
+ Description:
+
+ Given the knots and coefficients of a B-spline representation, create a
+ new B-spline with a knot inserted m times at point x.
+ This is a wrapper around the FORTRAN routine insert of FITPACK.
+
+ Inputs:
+
+ x (u) -- A 1-D point at which to insert a new knot(s). If tck was returned
+ from splprep, then the parameter values, u should be given.
+ tck -- A sequence of length 3 returned by splrep or splprep containg the
+ knots, coefficients, and degree of the spline.
+ m -- The number of times to insert the given knot (its multiplicity).
+ per -- If non-zero, input spline is considered periodic.
+
+ Outputs: tck
+
+ tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
+ coefficients, and the degree of the new spline.
+
+ Requirements:
+ t(k+1) <= x <= t(n-k), where k is the degree of the spline.
+ In case of a periodic spline (per != 0) there must be
+ either at least k interior knots t(j) satisfying t(k+1)<t(j)<=x
+ or at least k interior knots t(j) satisfying x<=t(j)<t(n-k).
+ """
+ t,c,k=tck
+ try:
+ c[0][0]
+ cc = []
+ for c_vals in c:
+ tt, cc_val, kk = insert(x, [t, c_vals, k], m)
+ cc.append(cc_val)
+ return (tt, cc, kk)
+ except: pass
+ tt, cc, ier = _fitpack._insert(per, t, c, k, x, m)
+ if ier==10: raise ValueError,"Invalid input data"
+ if ier: raise TypeError,"An error occurred"
+ return (tt, cc, k)
+
+
if __name__ == "__main__":
import sys,string
runtest=range(10)
More information about the Scipy-svn
mailing list