[Scipy-svn] r2634 - trunk/Lib/interpolate

scipy-svn at scipy.org scipy-svn at scipy.org
Tue Jan 30 06:16:36 CST 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