[Scipy-svn] r5232 - in trunk/scipy/interpolate: . fitpack src tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sun Dec 7 12:45:14 EST 2008
Author: ptvirtan
Date: 2008-12-07 11:44:59 -0600 (Sun, 07 Dec 2008)
New Revision: 5232
Added:
trunk/scipy/interpolate/fitpack/bispeu.f
Modified:
trunk/scipy/interpolate/fitpack2.py
trunk/scipy/interpolate/src/fitpack.pyf
trunk/scipy/interpolate/tests/test_fitpack.py
Log:
interpolate: evaluating bivariate splines at arbitrary points
Added: trunk/scipy/interpolate/fitpack/bispeu.f
===================================================================
--- trunk/scipy/interpolate/fitpack/bispeu.f 2008-12-07 17:30:50 UTC (rev 5231)
+++ trunk/scipy/interpolate/fitpack/bispeu.f 2008-12-07 17:44:59 UTC (rev 5232)
@@ -0,0 +1,64 @@
+ subroutine bispeu(tx,nx,ty,ny,c,kx,ky,x,y,z,m,wrk,lwrk, ier)
+c subroutine bispeu evaluates on a set of points (x(i),y(i)),i=1,...,m
+c a bivariate spline s(x,y) of degrees kx and ky, given in the
+c b-spline representation.
+c
+c calling sequence:
+c call bispeu(tx,nx,ty,ny,c,kx,ky,x,y,z,m,wrk,lwrk,
+c * iwrk,kwrk,ier)
+c
+c input parameters:
+c tx : real array, length nx, which contains the position of the
+c knots in the x-direction.
+c nx : integer, giving the total number of knots in the x-direction
+c ty : real array, length ny, which contains the position of the
+c knots in the y-direction.
+c ny : integer, giving the total number of knots in the y-direction
+c c : real array, length (nx-kx-1)*(ny-ky-1), which contains the
+c b-spline coefficients.
+c kx,ky : integer values, giving the degrees of the spline.
+c x : real array of dimension (mx).
+c y : real array of dimension (my).
+c m : on entry m must specify the number points. m >= 1.
+c wrk : real array of dimension lwrk. used as workspace.
+c lwrk : integer, specifying the dimension of wrk.
+c lwrk >= kx+ky+2
+c
+c output parameters:
+c z : real array of dimension m.
+c on succesful exit z(i) contains the value of s(x,y)
+c at the point (x(i),y(i)), i=1,...,m.
+c ier : integer error flag
+c ier=0 : normal return
+c ier=10: invalid input data (see restrictions)
+c
+c restrictions:
+c m >=1, lwrk>=mx*(kx+1)+my*(ky+1), kwrk>=mx+my
+c tx(kx+1) <= x(i-1) <= x(i) <= tx(nx-kx), i=2,...,mx
+c ty(ky+1) <= y(j-1) <= y(j) <= ty(ny-ky), j=2,...,my
+c
+c other subroutines required:
+c fpbisp,fpbspl
+c
+c ..scalar arguments..
+ integer nx,ny,kx,ky,m,lwrk,kwrk,ier
+c ..array arguments..
+ real*8 tx(nx),ty(ny),c((nx-kx-1)*(ny-ky-1)),x(m),y(m),z(m),
+ * wrk(lwrk)
+c ..local scalars..
+ integer iwrk(2)
+ integer i,iw,lwest
+c ..
+c before starting computations a data check is made. if the input data
+c are invalid control is immediately repassed to the calling program.
+ ier = 10
+ lwest = kx+ky+2
+ if (lwrk.lt.lwest) go to 100
+ if (m.lt.1) go to 100
+ ier = 0
+ do 10 i=1,m
+ call fpbisp(tx,nx,ty,ny,c,kx,ky,x(i),1,y(i),1,z(i),wrk(1),
+ * wrk(kx+2),iwrk(1),iwrk(2))
+ 10 continue
+ 100 return
+ end
Modified: trunk/scipy/interpolate/fitpack2.py
===================================================================
--- trunk/scipy/interpolate/fitpack2.py 2008-12-07 17:30:50 UTC (rev 5231)
+++ trunk/scipy/interpolate/fitpack2.py 2008-12-07 17:44:59 UTC (rev 5232)
@@ -362,6 +362,16 @@
return z
raise NotImplementedError
+ def ev(self, xi, yi):
+ """
+ Evaluate spline at points (x[i], y[i]), i=0,...,len(x)-1
+ """
+ tx,ty,c = self.tck[:3]
+ kx,ky = self.degrees
+ zi,ier = dfitpack.bispeu(tx,ty,c,kx,ky,xi,yi)
+ assert ier==0, 'Invalid input: ier='+`ier`
+ return zi
+
def integral(self, xa, xb, ya, yb):
"""
Evaluate the integral of the spline over area [xa,xb] x [ya,yb].
Modified: trunk/scipy/interpolate/src/fitpack.pyf
===================================================================
--- trunk/scipy/interpolate/src/fitpack.pyf 2008-12-07 17:30:50 UTC (rev 5231)
+++ trunk/scipy/interpolate/src/fitpack.pyf 2008-12-07 17:44:59 UTC (rev 5232)
@@ -329,6 +329,25 @@
integer intent(out) :: ier
end subroutine bispev
+ subroutine bispeu(tx,nx,ty,ny,c,kx,ky,x,y,z,m,wrk,lwrk,ier)
+ ! z,ier = bispeu(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(m) :: x
+ real*8 intent(in),dimension(m) :: y
+ integer intent(hide),depend(x) :: m=len(x)
+ real*8 dimension(m),depend(m),intent(out,c) :: z
+ real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
+ integer intent(hide),depend(kx,ky) :: lwrk=kx+ky+2
+ integer intent(out) :: ier
+ end subroutine bispeu
+
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)
Modified: trunk/scipy/interpolate/tests/test_fitpack.py
===================================================================
--- trunk/scipy/interpolate/tests/test_fitpack.py 2008-12-07 17:30:50 UTC (rev 5231)
+++ trunk/scipy/interpolate/tests/test_fitpack.py 2008-12-07 17:44:59 UTC (rev 5232)
@@ -150,5 +150,18 @@
lut = RectBivariateSpline(x,y,z)
assert_array_almost_equal(lut(x,y),z)
+ def test_evaluate(self):
+ x = array([1,2,3,4,5])
+ y = array([1,2,3,4,5])
+ z = array([[1,2,1,2,1],[1,2,1,2,1],[1,2,3,2,1],[1,2,2,2,1],[1,2,1,2,1]])
+ lut = RectBivariateSpline(x,y,z)
+
+ xi = [1, 2.3, 5.3, 0.5, 3.3, 1.2, 3]
+ yi = [1, 3.3, 1.2, 4.0, 5.0, 1.0, 3]
+ zi = lut.ev(xi, yi)
+ zi2 = array([lut(xp, yp)[0,0] for xp, yp in zip(xi, yi)])
+
+ assert_almost_equal(zi, zi2)
+
if __name__ == "__main__":
run_module_suite()
More information about the Scipy-svn
mailing list