[Scipy-svn] r2694 - in trunk/Lib: interpolate interpolate/fitpack sandbox/spline sandbox/spline/fitpack sandbox/spline/tests
scipy-svn@scip...
scipy-svn@scip...
Thu Feb 8 14:42:36 CST 2007
Author: jtravs
Date: 2007-02-08 14:42:26 -0600 (Thu, 08 Feb 2007)
New Revision: 2694
Added:
trunk/Lib/sandbox/spline/tests/visual_check_parametric.py
Modified:
trunk/Lib/interpolate/fitpack.py
trunk/Lib/interpolate/fitpack/insert.f
trunk/Lib/sandbox/spline/fitpack.py
trunk/Lib/sandbox/spline/fitpack.pyf
trunk/Lib/sandbox/spline/fitpack/insert.f
Log:
Add patch to fitpack/insert.f for Zach Pincus. Add insert function to sandbox.spline.
Modified: trunk/Lib/interpolate/fitpack/insert.f
===================================================================
--- trunk/Lib/interpolate/fitpack/insert.f 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/interpolate/fitpack/insert.f 2007-02-08 20:42:26 UTC (rev 2694)
@@ -61,7 +61,7 @@
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
-c latest update : march 1987
+c latest update : february 2007 (second interval search added)
c
c ..scalar arguments..
integer iopt,n,k,nn,nest,ier
@@ -69,7 +69,7 @@
c ..array arguments..
real*8 t(nest),c(nest),tt(nest),cc(nest)
c ..local scalars..
- integer kk,k1,l,nk,nk1
+ integer kk,k1,l,nk
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.
@@ -79,11 +79,18 @@
nk = n-k
if(x.lt.t(k1) .or. x.gt.t(nk)) go to 40
c search for knot interval t(l) <= x < t(l+1).
- nk1 = nk-1
l = k1
- 10 if(x.lt.t(l+1) .or. l.eq.nk1) go to 20
+ 10 if(x.lt.t(l+1)) go to 20
l = l+1
+ if(l.eq.nk) go to 14
go to 10
+c if no interval found above, then reverse the search and
+c look for knot interval t(l) < x <= t(l+1).
+ 14 l = nk-1
+ 16 if(x.gt.t(l)) go to 20
+ l = l-1
+ if(l.eq.k) go to 40
+ go to 16
20 if(t(l).ge.t(l+1)) go to 40
if(iopt.eq.0) go to 30
kk = 2*k
Modified: trunk/Lib/interpolate/fitpack.py
===================================================================
--- trunk/Lib/interpolate/fitpack.py 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/interpolate/fitpack.py 2007-02-08 20:42:26 UTC (rev 2694)
@@ -818,7 +818,6 @@
if ier: raise TypeError,"An error occurred"
return (tt, cc, k)
-
if __name__ == "__main__":
import sys,string
runtest=range(10)
Modified: trunk/Lib/sandbox/spline/fitpack/insert.f
===================================================================
--- trunk/Lib/sandbox/spline/fitpack/insert.f 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/sandbox/spline/fitpack/insert.f 2007-02-08 20:42:26 UTC (rev 2694)
@@ -61,7 +61,7 @@
c celestijnenlaan 200a, b-3001 heverlee, belgium.
c e-mail : Paul.Dierckx@cs.kuleuven.ac.be
c
-c latest update : march 1987
+c latest update : february 2007 (second interval search added)
c
c ..scalar arguments..
integer iopt,n,k,nn,nest,ier
@@ -69,7 +69,7 @@
c ..array arguments..
real*8 t(nest),c(nest),tt(nest),cc(nest)
c ..local scalars..
- integer kk,k1,l,nk,nk1
+ integer kk,k1,l,nk
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.
@@ -79,11 +79,18 @@
nk = n-k
if(x.lt.t(k1) .or. x.gt.t(nk)) go to 40
c search for knot interval t(l) <= x < t(l+1).
- nk1 = nk-1
l = k1
- 10 if(x.lt.t(l+1) .or. l.eq.nk1) go to 20
+ 10 if(x.lt.t(l+1)) go to 20
l = l+1
+ if(l.eq.nk) go to 14
go to 10
+c if no interval found above, then reverse the search and
+c look for knot interval t(l) < x <= t(l+1).
+ 14 l = nk-1
+ 16 if(x.gt.t(l)) go to 20
+ l = l-1
+ if(l.eq.k) go to 40
+ go to 16
20 if(t(l).ge.t(l+1)) go to 40
if(iopt.eq.0) go to 30
kk = 2*k
Modified: trunk/Lib/sandbox/spline/fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.py 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/sandbox/spline/fitpack.py 2007-02-08 20:42:26 UTC (rev 2694)
@@ -30,11 +30,11 @@
"""
__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
- 'bisplrep', 'bisplev']
+ 'bisplrep', 'bisplev', 'insert']
__version__ = "$Revision$"[10:-1]
from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
- dot, sin, cos, pi, arange, empty, int32, where
+ dot, sin, cos, pi, arange, empty, int32, where, resize
myasarray = atleast_1d
# f2py-generated interface to fitpack
@@ -210,7 +210,11 @@
else: nest=m+k+1
nest=max(nest,2*k+3)
if task==0:
- u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth0(ipar,idim,u,x,
+ if per:
+ u,n,t,c,fp,wrk,iwrk,ier=dfitpack.clocur_smth0(ipar,idim,u,x,w,nest,
+ k=k,s=s)
+ else:
+ u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth0(ipar,idim,u,x,
w,ub,ue,nest,k=k,s=s)
if task==1:
try:
@@ -589,7 +593,63 @@
raise TypeError,"Invalid input data. t(k)<=x<=t(n-k+1) must hold."
raise TypeError,"Unknown error"
+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]
+ parametric = True
+ except:
+ parametric = False
+ if parametric:
+ 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)
+ else:
+ n = len(t)
+ nest = n+m
+ c = c.copy()
+ c.resize(nest)
+ t = t.copy()
+ t.resize(nest)
+ while(n<nest):
+ tt, nn, cc, ier = dfitpack.insert(per, t, n, c, k, x)
+ if ier==10: raise ValueError,"Invalid input data"
+ if ier: raise TypeError,"An error occurred"
+ t = tt
+ c = cc
+ n+=1
+ return (tt, cc, k)
+
_surfit_cache = {}
def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,
Modified: trunk/Lib/sandbox/spline/fitpack.pyf
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.pyf 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/sandbox/spline/fitpack.pyf 2007-02-08 20:42:26 UTC (rev 2694)
@@ -256,10 +256,34 @@
integer intent(out) :: ier
end subroutine parcur_smth1
+ subroutine clocur_smth0(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,&
+ wrk,lwrk,iwrk,ier)
+ !u,n,t,c,fp,wrk,iwrk,ier=clocur_smth0(ipar,idim,u,x,w,nest,[k,s])
+ fortranname clocur
+ integer intent(hide) :: iopt = 0
+ integer check(ipar == 1 || ipar == 0) :: ipar
+ integer check(idim > 0 && idim < 11) :: idim
+ integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+ real*8 dimension(m), intent(in,out) :: u
+ integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+ real*8 dimension(mx) :: x
+ real*8 dimension(m) :: w
+ integer check(1<=k && k<=5) :: k=3.0
+ real*8 check(s>=0.0) :: s = 0.0
+ integer intent(in) :: nest
+ integer intent(out) :: n
+ real*8 dimension(nest), intent(out) :: t
+ integer intent(hide), depend(nest,idim) :: nc=idim*nest
+ real*8 dimension(nc), intent(out) :: c
+ real*8 intent(out) :: fp
+ real*8 dimension(lwrk), intent(out) :: wrk
+ integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(7+idim+5*k)
+ integer dimension(nest), intent(out) :: iwrk
+ integer intent(out) :: ier
+ end subroutine clocur_smth0
+
subroutine parcur_smth0(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,&
c,fp,wrk,lwrk,iwrk,ier)
- !u,ub,ue,n,t,c,fp,wrk,iwrk,ier=parcur_smth0(ipar,idim,u,x,w,
- ! ub,ue,nest,[k,s])
fortranname parcur
integer intent(hide) :: iopt = 0
integer check(ipar == 1 || ipar == 0) :: ipar
@@ -285,6 +309,21 @@
integer intent(out) :: ier
end subroutine parcur_smth0
+ subroutine insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
+ ! tt, nn, cc, ier = insert(per, t, c, k, x, nest)
+ integer intent(in):: iopt
+ real*8 dimension(nest),intent(in) :: t
+ integer intent(in) :: n
+ real*8 dimension(nest),depend(nest),intent(in) :: c
+ integer intent(in) :: k
+ real*8 intent(in) :: x
+ real*8 dimension(nest),depend(nest),intent(out) :: tt
+ integer intent(out) :: nn
+ real*8 dimension(nest),depend(nest),intent(out) :: cc
+ integer intent(hide),depend(t) :: nest=len(t)
+ integer intent(out) :: ier
+ end subroutine insert
+
subroutine fpcurf_smth0(iopt,x,y,w,m,xb,xe,k,s,nest,tol,maxit,k1,k2,n,t,&
c,fp,fpint,wrk,nrdata,ier)
! x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier = \
Added: trunk/Lib/sandbox/spline/tests/visual_check_parametric.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/visual_check_parametric.py 2007-02-08 17:07:46 UTC (rev 2693)
+++ trunk/Lib/sandbox/spline/tests/visual_check_parametric.py 2007-02-08 20:42:26 UTC (rev 2694)
@@ -0,0 +1,51 @@
+from scipy import arange, cos, linspace, randn, pi, sin
+from scipy.sandbox.spline import splprep, splev
+import matplotlib
+matplotlib.use('WXAgg')
+import pylab
+
+# make ascending spiral in 3-space
+t=linspace(0,1.75*2*pi,100)
+
+x = sin(t)
+y = cos(t)
+z = t
+
+# add noise
+x+= 0.1*randn(*x.shape)
+y+= 0.1*randn(*y.shape)
+z+= 0.1*randn(*z.shape)
+
+# spline parameters
+s=3.0 # smoothness parameter
+k=2 # spline order
+nest=-1 # estimate of number of knots needed (-1 = maximal)
+
+# find the knot points
+tckp,u = splprep([x,y,z],s=s,k=k,nest=-1)
+
+# evaluate spline, including interpolated points
+xnew,ynew,znew = splev(linspace(0,1,400),tckp)
+
+pylab.subplot(2,2,1)
+data,=pylab.plot(x,y,'bo-',label='data')
+fit,=pylab.plot(xnew,ynew,'r-',label='fit')
+pylab.legend()
+pylab.xlabel('x')
+pylab.ylabel('y')
+
+pylab.subplot(2,2,2)
+data,=pylab.plot(x,z,'bo-',label='data')
+fit,=pylab.plot(xnew,znew,'r-',label='fit')
+pylab.legend()
+pylab.xlabel('x')
+pylab.ylabel('z')
+
+pylab.subplot(2,2,3)
+data,=pylab.plot(y,z,'bo-',label='data')
+fit,=pylab.plot(ynew,znew,'r-',label='fit')
+pylab.legend()
+pylab.xlabel('y')
+pylab.ylabel('z')
+
+pylab.show()
\ No newline at end of file
More information about the Scipy-svn
mailing list