[Scipy-svn] r6622 - in trunk/scipy/interpolate: . fitpack src
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Jul 18 15:21:26 EDT 2010
Author: charris
Date: 2010-07-18 14:21:26 -0500 (Sun, 18 Jul 2010)
New Revision: 6622
Modified:
trunk/scipy/interpolate/fitpack.py
trunk/scipy/interpolate/fitpack/splev.f
trunk/scipy/interpolate/src/__fitpack.h
trunk/scipy/interpolate/src/fitpack.pyf
Log:
ENH: Add extrapolate keyword to splev to control the its behavior
outside of the interval of definition. If extrapolate=1, then the
current behavior is preserved, if extrapolate=0, then the spline
evaluates to zero for points outside of the interval of definition.
The default value is 1.
Modified: trunk/scipy/interpolate/fitpack/splev.f
===================================================================
--- trunk/scipy/interpolate/fitpack/splev.f 2010-07-18 19:21:22 UTC (rev 6621)
+++ trunk/scipy/interpolate/fitpack/splev.f 2010-07-18 19:21:26 UTC (rev 6622)
@@ -1,9 +1,9 @@
- subroutine splev(t,n,c,k,x,y,m,ier)
+ subroutine splev(t,n,c,k,x,y,m,e,ier)
c subroutine splev evaluates in a number of points x(i),i=1,2,...,m
c a spline s(x) of degree k, given in its b-spline representation.
c
c calling sequence:
-c call splev(t,n,c,k,x,y,m,ier)
+c call splev(t,n,c,k,x,y,m,e,ier)
c
c input parameters:
c t : array,length n, which contains the position of the knots.
@@ -14,6 +14,9 @@
c be evaluated.
c m : integer, giving the number of points where s(x) must be
c evaluated.
+c e : integer, if != 0 the spline is extrapolated from the end
+c spans for points not in the support, if == 0 the spline
+c evaluates to zero for those points.
c
c output parameter:
c y : array,length m, giving the value of s(x) at the different
@@ -50,15 +53,15 @@
c++ - fixed initialization of sp to double precision value
c
c ..scalar arguments..
- integer n,k,m,ier
+ integer n, k, m, e, ier
c ..array arguments..
- real*8 t(n),c(n),x(m),y(m)
+ real*8 t(n), c(n), x(m), y(m)
c ..local scalars..
- integer i,j,k1,l,ll,l1,nk1
+ integer i, j, k1, l, ll, l1, nk1
c++..
integer k2
c..++
- real*8 arg,sp,tb,te
+ real*8 arg, sp, tb, te
c ..local array..
real*8 h(20)
c ..
@@ -67,47 +70,49 @@
ier = 10
c-- if(m-1) 100,30,10
c++..
- if(m.lt.1) go to 100
+ if (m .lt. 1) go to 100
c..++
c-- 10 do 20 i=2,m
c-- if(x(i).lt.x(i-1)) go to 100
c-- 20 continue
30 ier = 0
c fetch tb and te, the boundaries of the approximation interval.
- k1 = k+1
+ k1 = k + 1
c++..
- k2 = k1+1
+ k2 = k1 + 1
c..++
- nk1 = n-k1
+ nk1 = n - k1
tb = t(k1)
- te = t(nk1+1)
+ te = t(nk1 + 1)
l = k1
- l1 = l+1
+ l1 = l + 1
c main loop for the different points.
- do 80 i=1,m
+ do 80 i = 1, m
c fetch a new x-value arg.
arg = x(i)
-c-- if(arg.lt.tb) arg = tb
-c-- if(arg.gt.te) arg = te
+c check if arg is in the support
+ if (e .eq. 1 .or. arg .ge. tb .and. arg .le. te) go to 35
+ y(i) = 0
+ goto 80
c search for knot interval t(l) <= arg < t(l+1)
c++..
- 35 if(arg.ge.t(l) .or. l1.eq.k2) go to 40
+ 35 if (arg .ge. t(l) .or. l1 .eq. k2) go to 40
l1 = l
- l = l-1
+ l = l - 1
go to 35
c..++
- 40 if(arg.lt.t(l1) .or. l.eq.nk1) go to 50
+ 40 if(arg .lt. t(l1) .or. l .eq. nk1) go to 50
l = l1
- l1 = l+1
+ l1 = l + 1
go to 40
c evaluate the non-zero b-splines at arg.
- 50 call fpbspl(t,n,k,arg,l,h)
+ 50 call fpbspl(t, n, k, arg, l, h)
c find the value of s(x) at x=arg.
sp = 0.0d0
- ll = l-k1
- do 60 j=1,k1
- ll = ll+1
- sp = sp+c(ll)*h(j)
+ ll = l - k1
+ do 60 j = 1, k1
+ ll = ll + 1
+ sp = sp + c(ll)*h(j)
60 continue
y(i) = sp
80 continue
Modified: trunk/scipy/interpolate/fitpack.py
===================================================================
--- trunk/scipy/interpolate/fitpack.py 2010-07-18 19:21:22 UTC (rev 6621)
+++ trunk/scipy/interpolate/fitpack.py 2010-07-18 19:21:26 UTC (rev 6622)
@@ -431,29 +431,37 @@
#if len(l)>1: return l
#return l[0]
-def splev(x,tck,der=0):
- """
- Evaluate a B-spline and its derivatives.
+def splev(x, tck, der=0, extrapolate=1):
+ """Evaluate a B-spline or its derivatives.
Given the knots and coefficients of a B-spline representation, evaluate
- the value of the smoothing polynomial and it's derivatives.
- This is a wrapper around the FORTRAN routines splev and splder of FITPACK.
+ the value of the smoothing polynomial and it's derivatives. This is a
+ wrapper around the FORTRAN routines splev and splder of FITPACK.
Parameters
----------
- x (u) -- a 1-D array of points at which to return the value of the
- smoothed spline or its derivatives. 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.
- der -- The order of derivative of the spline to compute (must be less than
- or equal to k).
+ x : array_like
+ A 1-D array of points at which to return the value of the smoothed
+ spline or its derivatives. If tck was returned from splprep, then
+ the parameter values, u should be given.
+ tck : tuple
+ A sequence of length 3 returned by splrep or splprep containg the
+ knots, coefficients, and degree of the spline.
+ der : int
+ The order of derivative of the spline to compute (must be less than
+ or equal to k).
+ extrapolate : int
+ If zero, points outside of the interval defined by the knot
+ sequence evaluate to zero, otherwise the piecewise polynomials in
+ the end spans of the knot sequence are extrapolated. The default
+ value is 1.
Returns
-------
- y -- an array of values representing the spline function or curve.
- If tck was returned from splrep, then this is a list of arrays
- representing the curve in N-dimensional space.
+ y : ndarray or list of ndarrays
+ An array of values representing the spline function evaluated at
+ the points in ``x``. If tck was returned from splrep, then this is
+ a list of arrays representing the curve in N-dimensional space.
See Also
--------
@@ -472,22 +480,25 @@
on Numerical Analysis, Oxford University Press, 1993.
"""
- t,c,k=tck
+ t,c,k = tck
try:
c[0][0]
parametric = True
except:
parametric = False
if parametric:
- return map(lambda c,x=x,t=t,k=k,der=der:splev(x,[t,c,k],der),c)
+ return map(lambda c, x=x, t=t, k=k, der=der : splev(x, [t,c,k], der), c)
else:
- if not (0<=der<=k):
- raise ValueError,"0<=der=%d<=k=%d must hold"%(der,k)
- x=myasarray(x)
- y,ier=_fitpack._spl_(x,der,t,c,k)
- if ier==10: raise ValueError,"Invalid input data"
- if ier: raise TypeError,"An error occurred"
- if len(y)>1: return y
+ if not (0 <= der <= k):
+ raise ValueError, "0<=der=%d<=k=%d must hold"%(der,k)
+ x = myasarray(x)
+ y, ier =_fitpack._spl_(x, der, t, c, k, extrapolate)
+ if ier == 10:
+ raise ValueError,"Invalid input data"
+ if ier:
+ raise TypeError,"An error occurred"
+ if len(y) > 1:
+ return y
return y[0]
def splint(a,b,tck,full_output=0):
Modified: trunk/scipy/interpolate/src/__fitpack.h
===================================================================
--- trunk/scipy/interpolate/src/__fitpack.h 2010-07-18 19:21:22 UTC (rev 6621)
+++ trunk/scipy/interpolate/src/__fitpack.h 2010-07-18 19:21:26 UTC (rev 6622)
@@ -78,7 +78,7 @@
void PERCUR(int*,int*,double*,double*,double*,int*,double*,int*,int*,double*,double*,double*,double*,int*,int*,int*);
void SPALDE(double*,int*,double*,int*,double*,double*,int*);
void SPLDER(double*,int*,double*,int*,int*,double*,double*,int*,double*,int*);
-void SPLEV(double*,int*,double*,int*,double*,double*,int*,int*);
+void SPLEV(double*,int*,double*,int*,double*,double*,int*,int*,int*);
double SPLINT(double*,int*,double*,int*,double*,double*,double*);
void SPROOT(double*,int*,double*,double*,int*,int*,int*);
void PARCUR(int*,int*,int*,int*,double*,int*,double*,double*,double*,double*,int*,double*,int*,int*,double*,int*,double*,double*,double*,int*,int*,int*);
@@ -446,16 +446,16 @@
return NULL;
}
-static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k )";
+static char doc_spl_[] = " [y,ier] = _spl_(x,nu,t,c,k,e)";
static PyObject *fitpack_spl_(PyObject *dummy, PyObject *args)
{
- int n, nu, ier, k;
+ int n, nu, ier, k, e=1;
npy_intp m;
double *x, *y, *t, *c, *wrk = NULL;
PyArrayObject *ap_x = NULL, *ap_y = NULL, *ap_t = NULL, *ap_c = NULL;
PyObject *x_py = NULL, *t_py = NULL, *c_py = NULL;
- if (!PyArg_ParseTuple(args, "OiOOi",&x_py,&nu,&t_py,&c_py,&k)) {
+ if (!PyArg_ParseTuple(args, "OiOOii", &x_py, &nu, &t_py, &c_py, &k, &e)) {
return NULL;
}
ap_x = (PyArrayObject *)PyArray_ContiguousFromObject(x_py, PyArray_DOUBLE, 0, 1);
@@ -482,7 +482,7 @@
SPLDER(t, &n, c, &k, &nu, x, y, &m, wrk, &ier);
}
else {
- SPLEV(t, &n, c, &k, x, y, &m, &ier);
+ SPLEV(t, &n, c, &k, x, y, &m, &e, &ier);
}
if (wrk) {
free(wrk);
Modified: trunk/scipy/interpolate/src/fitpack.pyf
===================================================================
--- trunk/scipy/interpolate/src/fitpack.pyf 2010-07-18 19:21:22 UTC (rev 6621)
+++ trunk/scipy/interpolate/src/fitpack.pyf 2010-07-18 19:21:26 UTC (rev 6622)
@@ -64,7 +64,7 @@
return u*v*(b2+1)+b2;
}
-static int calc_regrid_lwrk(int mx, int my, int kx, int ky,
+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;
@@ -75,8 +75,8 @@
!!!!!!!!!! Univariate spline !!!!!!!!!!!
- subroutine splev(t,n,c,k,x,y,m,ier)
- ! y = splev(t,c,k,x)
+ subroutine splev(t,n,c,k,x,y,m,e,ier)
+ ! y = splev(t,c,k,x,e)
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
@@ -84,6 +84,7 @@
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 :: e=1
integer intent(hide) :: ier
end subroutine splev
@@ -161,7 +162,7 @@
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)
+ 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)
@@ -180,9 +181,9 @@
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)
+
+ 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
@@ -436,7 +437,7 @@
:: 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])
@@ -446,7 +447,7 @@
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)
+ 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)
@@ -475,7 +476,7 @@
:: 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
More information about the Scipy-svn
mailing list