[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