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

[scipy-svn@scip...]

Author: oliphant
Date: 2007-06-01 07:09:10 -0500 (Fri, 01 Jun 2007)
New Revision: 3068

Modified:
   trunk/Lib/interpolate/__fitpack.h
   trunk/Lib/interpolate/interpolate.py
Log:
Add Lagrange interpolating polynomial

Modified: trunk/Lib/interpolate/__fitpack.h
===================================================================
--- trunk/Lib/interpolate/__fitpack.h	2007-06-01 08:47:46 UTC (rev 3067)
+++ trunk/Lib/interpolate/__fitpack.h	2007-06-01 12:09:10 UTC (rev 3068)
@@ -829,7 +829,7 @@
 "integer-spaced, or cardinal spline matrix a bit faster.";
 static PyObject *_bsplmat(PyObject *dummy, PyObject *args) {
     int k,N,i,numbytes,j, equal;
-    int dims[2];
+    npy_intp dims[2];
     PyObject *x_i_py=NULL;
     PyArrayObject *x_i=NULL, *BB=NULL;
     double *t=NULL, *h=NULL, *ptr;
@@ -970,7 +970,7 @@
 "then it produces the result as if the sample distance were dx";
 static PyObject *_bspldismat(PyObject *dummy, PyObject *args) {
     int k,N,i,j, equal, m;
-    int dims[2];
+    npy_intp dims[2];
     PyObject *x_i_py=NULL;
     PyArrayObject *x_i=NULL, *BB=NULL;
     double *t=NULL, *h=NULL, *ptr, *dptr;

Modified: trunk/Lib/interpolate/interpolate.py
===================================================================
--- trunk/Lib/interpolate/interpolate.py	2007-06-01 08:47:46 UTC (rev 3067)
+++ trunk/Lib/interpolate/interpolate.py	2007-06-01 12:09:10 UTC (rev 3068)
@@ -4,7 +4,7 @@
 """
 
 __all__ = ['interp1d', 'interp2d', 'spline', 'spleval', 'splmake', 'spltopp',
-           'ppform']
+           'ppform', 'lagrange']
 
 from numpy import shape, sometrue, rank, array, transpose, \
      swapaxes, searchsorted, clip, take, ones, putmask, less, greater, \
@@ -23,6 +23,21 @@
         all = sometrue(all,axis=0)
     return all
 
+def lagrange(x, w):
+    """Return the Lagrange interpolating polynomial of the data-points (x,w)
+    """
+    M = len(x)
+    p = poly1d(0.0)
+    for j in xrange(M):
+        pt = poly1d(w[j])
+        for k in xrange(M):
+            if k == j: continue
+            fac = x[j]-x[k]
+            pt *= poly1d([1.0,-x[k]])/fac
+        p += pt
+    return p
+
+
 # !! Need to find argument for keeping initialize.  If it isn't
 # !! found, get rid of it!