[Scipy-svn] r4591 - in branches/Interpolate1D: . docs

scipy-svn@scip... scipy-svn@scip...
Thu Jul 31 18:58:47 CDT 2008


Author: fcady
Date: 2008-07-31 18:58:42 -0500 (Thu, 31 Jul 2008)
New Revision: 4591

Added:
   branches/Interpolate1D/_fitpack.pyf
Removed:
   branches/Interpolate1D/fitpack.pyf
Modified:
   branches/Interpolate1D/TODO.txt
   branches/Interpolate1D/docs/tutorial.rst
   branches/Interpolate1D/example_script.py
   branches/Interpolate1D/fitpack_wrapper.py
   branches/Interpolate1D/interpolate1d.py
   branches/Interpolate1D/interpolate_wrapper.py
   branches/Interpolate1D/setup.py
Log:
fitpack.pyf to _fitpack.pyf to make it more in the background, added 'nearest' interpolate method for more functionality

Modified: branches/Interpolate1D/TODO.txt
===================================================================
--- branches/Interpolate1D/TODO.txt	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/TODO.txt	2008-07-31 23:58:42 UTC (rev 4591)
@@ -15,11 +15,11 @@
     In fitpack_wrapper the smoothing parameter is s.  It now defaults
     to 0.0 (exact interpolation).  Zero smoothing and moderate (s ~ 1)
     are fast, but small non-zero s makes the algorithm VERY slow.
+    
+    This appears resolved, in that Spline can do smoothing (but defaults
+    to not), and user-defined classes are allowed to smooth.
 
-    Currently the user does see s in Interpolate1d; no smoothing 
-    is available.  The Spline class, however, has it available (default = 0.0).
 
-
 **pick best spline
     Under-the-hood machinery currently comes from _interpolate.cpp
     (used in enthought.interpolate) and FITPACK (Fortran, used in 
@@ -41,6 +41,21 @@
     that stuff needs to be cleaned up, and the
     rest needs to be either removed or set aside.
     
+**allow y to be 2-dimensional?
+    It is not decided whether this feature should be supported, but
+    I don't think it should; I think there should be another class with
+    that functionality which wraps Interpolate1d.  The main reasons
+    are:
+    1) interpolation should probably be along columns, so the user
+    can enter data as enter the data as y=array([obs1, obs2, ...]).  But
+    the interpolate_wrapper and fitpack_wrapper interpolate along
+    rows; this means there would be messy transposes in the code.
+    2) FITPACK doesn't support 2D y arrays, which means we would
+    have to store the y columns separately
+    3) If we remove bad data, bearing in mind (2), we would also
+    have to store copies of x, since the bad data in the y cols aren't
+    at the same locations, meaning that different data points should be
+    deleted for each column.
 
 **better handling of variable types
     Currently everything is cast to a float64 if it is not already
@@ -56,7 +71,7 @@
     
     
 
-*********** DOCUENTATION-TYPE AND NON-URGENT TASKS *******
+*********** DOCUMENTATION-TYPE AND NON-URGENT TASKS *******
 
 **improve regression tests
     desired for fitpack_wrapper and _interpolate_wrapper
@@ -64,10 +79,7 @@
     really basic.
 
 
-**improve unit tests
-    interpolate1d.py has its own unit tests, which also call
-    the unit tests of fitpack_wrapper and _interpolate_wrapper.
-    But they could surely be made better.
+**improve unit tests in tests directory
 
 
 **comment all files
@@ -102,7 +114,7 @@
 ********* LONGER TERM ************
 
 **update for 2D and ND
-    This will probably take the form of two additional
+    This will take the form of two additional
     classes both modelled after interpolate1d.  Thus it probably
     shouldn't be done until interpolate1d is more settled.
     
@@ -147,20 +159,3 @@
     differentiation of functionality (one for splines, one for simple
     operations, etc), rather than being a holdover of where they
     were stolen from.
-
-
-**allow y to be 2-dimensional?
-    It is not decided whether this feature should be supported, but
-    I don't think it should; I think there should be another class with
-    that functionality which wraps Interpolate1d.  The main reasons
-    are:
-    1) interpolation should probably be along columns, so the user
-    can enter data as enter the data as y=array([obs1, obs2, ...]).  But
-    the interpolate_wrapper and fitpack_wrapper interpolate along
-    rows; this means there would be messy transposes in the code.
-    2) FITPACK doesn't support 2D y arrays, which means we would
-    have to store the y columns separately
-    3) If we remove bad data, bearing in mind (2), we would also
-    have to store copies of x, since the bad data in the y cols aren't
-    at the same locations, meaning that different data points should be
-    deleted for each column.

Copied: branches/Interpolate1D/_fitpack.pyf (from rev 4587, branches/Interpolate1D/fitpack.pyf)

Modified: branches/Interpolate1D/docs/tutorial.rst
===================================================================
--- branches/Interpolate1D/docs/tutorial.rst	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/docs/tutorial.rst	2008-07-31 23:58:42 UTC (rev 4591)
@@ -481,6 +481,9 @@
 #) s 
     If s is zero, the interpolation is exact.  If s is not 0, the curve is smoothe subject to
     the constraint that sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
+    
+    BEWARE : in the current implementation of the code, if s is small but not zero,
+            instantiating Spline can become painfully slow.
 
 At calling:
 

Modified: branches/Interpolate1D/example_script.py
===================================================================
--- branches/Interpolate1D/example_script.py	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/example_script.py	2008-07-31 23:58:42 UTC (rev 4591)
@@ -82,10 +82,6 @@
     interp = I.Interpolate1d(x, y, fitpack_wrapper.Spline(k=2))
     y_quad2 = interp(newx)
 
-    # 3rd order spline with additional keyword arguments
-    interp = I.Interpolate1d(x, y, fitpack_wrapper.Spline, kindkw = {'k':3})
-    y_cubic2 = interp(newx)
-
     # 4th order spline
     interp = I.Interpolate1d(x, y, 'Quartic')
     y_quartic2 = interp(newx)
@@ -100,7 +96,6 @@
     P.plot(newx, y_block2, 'g')
     P.plot(newx, y_linear2, 'b')
     P.plot(newx, y_quad2, 'r')
-    P.plot(newx, y_cubic2, 'm')
     P.plot(newx, y_quartic2, 'y')
     P.plot(newx, y_quintic2, 'y')
     P.title( "same data through different interface" )

Deleted: branches/Interpolate1D/fitpack.pyf
===================================================================
--- branches/Interpolate1D/fitpack.pyf	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/fitpack.pyf	2008-07-31 23:58:42 UTC (rev 4591)
@@ -1,479 +0,0 @@
-!    -*- f90 -*-
-! Author: Pearu Peterson <pearu@cens.ioc.ee>
-!
-python module _dfitpack ! in
-
-  usercode '''
-
-static double dmax(double* seq,int len) {
-  double val;
-  int i;
-  if (len<1)
-    return -1e308;
-  val = seq[0];
-  for(i=1;i<len;++i)
-    if (seq[i]>val) val = seq[i];
-  return val;
-}
-static double dmin(double* seq,int len) {
-  double val;
-  int i;
-  if (len<1)
-    return 1e308;
-  val = seq[0];
-  for(i=1;i<len;++i)
-    if (seq[i]<val) val = seq[i];
-  return val;
-}
-static double calc_b(double* x,int m,double* tx,int nx) {
-  double val1 = dmin(x,m);
-  double val2 = dmin(tx,nx);
-  if (val2>val1) return val1;
-  val1 = dmax(tx,nx);
-  return val2 - (val1-val2)/nx;
-}
-static double calc_e(double* x,int m,double* tx,int nx) {
-  double val1 = dmax(x,m);
-  double val2 = dmax(tx,nx);
-  if (val2<val1) return val1;
-  val1 = dmin(tx,nx);
-  return val2 + (val2-val1)/nx;
-}
-static int imax(int i1,int i2) {
-  return MAX(i1,i2);
-}
-
-static int calc_surfit_lwrk1(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int km = MAX(kx,ky)+1;
- int ne = MAX(nxest,nyest);
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b1,b2;
- if (bx<=by) {b1=bx;b2=bx+v-ky;}
- else {b1=by;b2=by+u-kx;}
- return u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1;
-}
-static int calc_surfit_lwrk2(int m, int kx, int ky, int nxest, int nyest) {
- int u = nxest-kx-1;
- int v = nyest-ky-1;
- int bx = kx*v+ky+1;
- int by = ky*u+kx+1;
- int b2 = (bx<=by?bx+v-ky:by+u-kx);
- return u*v*(b2+1)+b2;
-}
-
-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;
-}
-'''
-
-  interface
-
-     !!!!!!!!!! Univariate spline !!!!!!!!!!!
-
-     subroutine splev(t,n,c,k,x,y,m,ier)
-       ! y = splev(t,c,k,x)
-       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
-       integer :: k
-       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 intent(hide) :: ier
-     end subroutine splev
-
-     subroutine splder(t,n,c,k,nu,x,y,m,wrk,ier)
-       ! dy = splder(t,c,k,x,[nu])
-       real*8 dimension(n) :: t
-       integer depend(t),intent(hide) :: n=len(t)
-       real*8 dimension(n),depend(n,k),check(len(c)==n),intent(in) :: c
-       integer :: k
-       integer depend(k),check(0<=nu && nu<=k) :: nu = 1
-       real*8 dimension(m) :: x
-       real*8 dimension(m),depend(m),intent(out) :: y
-       integer depend(x),intent(hide) :: m=len(x)
-       real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
-       integer intent(hide) :: ier
-     end subroutine splder
-
-     function splint(t,n,c,k,a,b,wrk)
-       ! iy = splint(t,c,k,a,b)
-       real*8 dimension(n),intent(in) :: t
-       integer intent(hide),depend(t) :: n=len(t)
-       real*8 dimension(n),depend(n),check(len(c)==n) :: c
-       integer intent(in) :: k
-       real*8 intent(in) :: a
-       real*8 intent(in) :: b
-       real*8 dimension(n),depend(n),intent(cache,hide) :: wrk
-       real*8 :: splint
-     end function splint
-
-     subroutine sproot(t,n,c,zero,mest,m,ier)
-       ! zero,m,ier = sproot(t,c[,mest])
-       real*8 dimension(n),intent(in) :: t
-       integer intent(hide),depend(t),check(n>=8) :: n=len(t)
-       real*8 dimension(n),depend(n),check(len(c)==n) :: c
-       real*8 dimension(mest),intent(out),depend(mest) :: zero
-       integer optional,intent(in),depend(n) :: mest=3*(n-7)
-       integer intent(out) :: m
-       integer intent(out) :: ier
-     end subroutine sproot
-
-     subroutine spalde(t,n,c,k,x,d,ier)
-       ! d,ier = spalde(t,c,k,x)
-
-       callprotoargument double*,int*,double*,int*,double*,double*,int*
-       callstatement {int k1=k+1; (*f2py_func)(t,&n,c,&k1,&x,d,&ier); }
-
-       real*8 dimension(n) :: t
-       integer intent(hide),depend(t) :: n=len(t)
-       real*8 dimension(n),depend(n),check(len(c)==n) :: c
-       integer intent(in) :: k
-       real*8 intent(in) :: x
-       real*8 dimension(k+1),intent(out),depend(k) :: d
-       integer intent(out) :: ier
-     end subroutine spalde
-
-     subroutine curfit(iopt,m,x,y,w,xb,xe,k,s,nest,n,t,c,fp,wrk,lwrk,iwrk,ier)
-       ! in  curfit.f
-       integer :: iopt
-       integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
-       real*8 dimension(m) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m) :: y
-       real*8 dimension(m),depend(m),check(len(w)==m) :: w
-       real*8 optional,depend(x),check(xb<=x[0]) :: xb = x[0]
-       real*8 optional,depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
-       integer optional,check(1<=k && k <=5),intent(in) :: k=3
-       real*8 optional,check(s>=0.0) :: s = 0.0
-       integer intent(hide),depend(t) :: nest=len(t)
-       integer intent(out), depend(nest) :: n=nest
-       real*8 dimension(nest),intent(inout) :: t
-       real*8 dimension(n),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk),intent(inout) :: wrk
-       integer intent(hide),depend(wrk) :: lwrk=len(wrk)
-       integer dimension(nest),intent(inout) :: iwrk
-       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) 
-       ! in percur.f
-       integer :: iopt
-       integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
-       real*8 dimension(m) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m) :: y
-       real*8 dimension(m),depend(m),check(len(w)==m) :: w
-       integer optional,check(1<=k && k <=5),intent(in) :: k=3
-       real*8 optional,check(s>=0.0) :: s = 0.0
-       integer intent(hide),depend(t) :: nest=len(t)
-       integer intent(out), depend(nest) :: n=nest
-       real*8 dimension(nest),intent(inout) :: t
-       real*8 dimension(n),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk),intent(inout) :: wrk
-       integer intent(hide),depend(wrk) :: lwrk=len(wrk)
-       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) 
-       ! in parcur.f
-       integer check(iopt>=-1 && iopt <= 1):: iopt
-       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(inout) :: 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
-       real*8 :: ub
-       real*8 :: ue
-       integer optional, check(1<=k && k<=5) :: k=3.0
-       real*8 optional, check(s>=0.0) :: s = 0.0
-       integer intent(hide), depend(t) :: nest=len(t)
-       integer intent(out), depend(nest) :: n=nest
-       real*8 dimension(nest), intent(inout) :: t
-       integer intent(hide), depend(c,nest,idim), check(nc>=idim*nest) :: nc=len(c)
-       real*8 dimension(nc), intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk), intent(inout) :: wrk
-       integer intent(hide),depend(wrk) :: lwrk=len(wrk)
-       integer dimension(nest), intent(inout) :: iwrk
-       integer intent(out) :: ier
-     end subroutine parcur
-
-
-     subroutine fpcurf0(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 = \
-       !   fpcurf0(x,y,k,[w,xb,xe,s,nest])
-
-       fortranname fpcurf
-       callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
-       callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
-       integer intent(hide) :: iopt = 0
-       real*8 dimension(m),intent(in,out) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
-       real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
-       integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
-       real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
-       real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
-       integer check(1<=k && k<=5),intent(in,out) :: k
-       real*8 check(s>=0.0),depend(m),intent(in,out) :: s = m
-       integer intent(in),depend(m,s,k,k1),check(nest>=2*k1) :: nest = (s==0.0?m+k+1:MAX(m/2,2*k1))
-       real*8 intent(hide) :: tol = 0.001
-       integer intent(hide) :: maxit = 20
-       integer intent(hide),depend(k) :: k1=k+1
-       integer intent(hide),depend(k) :: k2=k+2
-       integer intent(out) :: n
-       real*8 dimension(nest),intent(out),depend(nest) :: t
-       real*8 dimension(nest),depend(nest),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(nest),depend(nest),intent(out,cache)  :: fpint
-       real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
-       integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
-       integer intent(out) :: ier
-     end subroutine fpcurf0
-
-     subroutine fpcurf1(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 = \
-       !   fpcurf1(x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier)
-
-       fortranname fpcurf
-       callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
-       callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
-       integer intent(hide) :: iopt = 1
-       real*8 dimension(m),intent(in,out,overwrite) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out,overwrite) :: y
-       real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out,overwrite) :: w
-       integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
-       real*8 intent(in,out) :: xb
-       real*8 intent(in,out) :: xe
-       integer check(1<=k && k<=5),intent(in,out) :: k
-       real*8 check(s>=0.0),intent(in,out) :: s
-       integer intent(hide),depend(t) :: nest = len(t)
-       real*8 intent(hide) :: tol = 0.001
-       integer intent(hide) :: maxit = 20
-       integer intent(hide),depend(k) :: k1=k+1
-       integer intent(hide),depend(k) :: k2=k+2
-       integer intent(in,out) :: n
-       real*8 dimension(nest),intent(in,out,overwrite) :: t
-       real*8 dimension(nest),depend(nest),check(len(c)==nest),intent(in,out,overwrite) :: c
-       real*8 intent(in,out) :: fp
-       real*8 dimension(nest),depend(nest),check(len(fpint)==nest),intent(in,out,cache,overwrite)  :: fpint
-       real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
-       integer dimension(nest),depend(nest),check(len(nrdata)==nest),intent(in,out,cache,overwrite) :: nrdata
-       integer intent(in,out) :: ier
-     end subroutine fpcurf1
-
-     subroutine fpcurfm1(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 = \
-       !   fpcurfm1(x,y,k,t,[w,xb,xe])
-
-       fortranname fpcurf
-       callprotoargument int*,double*,double*,double*,int*,double*,double*,int*,double*,int*,double*,int*,int*,int*,int*,double*,double*,double*,double*,double*,double*,double*,double*,double*,int*,int*
-       callstatement (*f2py_func)(&iopt,x,y,w,&m,&xb,&xe,&k,&s,&nest,&tol,&maxit,&k1,&k2,&n,t,c,&fp,fpint,wrk,wrk+nest,wrk+nest*k2,wrk+nest*2*k2,wrk+nest*3*k2,nrdata,&ier)
-
-       integer intent(hide) :: iopt = -1
-       real*8 dimension(m),intent(in,out) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m),intent(in,out) :: y
-       real*8 dimension(m),depend(m),check(len(w)==m),intent(in,out) :: w = 1.0
-       integer intent(hide),depend(x),check(m>k),depend(k) :: m=len(x)
-       real*8 intent(in,out),depend(x),check(xb<=x[0]) :: xb = x[0]
-       real*8 intent(in,out),depend(x,m),check(xe>=x[m-1]) :: xe = x[m-1]
-       integer check(1<=k && k<=5),intent(in,out) :: k
-       real*8 intent(out) :: s = -1
-       integer intent(hide),depend(n) :: nest = n
-       real*8 intent(hide) :: tol = 0.001
-       integer intent(hide) :: maxit = 20
-       integer intent(hide),depend(k) :: k1=k+1
-       integer intent(hide),depend(k) :: k2=k+2
-       integer intent(out),depend(t) :: n = len(t)
-       real*8 dimension(n),intent(in,out,overwrite) :: t
-       real*8 dimension(nest),depend(nest),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(nest),depend(nest),intent(out,cache)  :: fpint
-       real*8 dimension(nest*3*k2+m*k1),intent(cache,hide),depend(nest,k1,k2,m) :: wrk
-       integer dimension(nest),depend(nest),intent(out,cache) :: nrdata
-       integer intent(out) :: ier
-     end subroutine fpcurfm1
-
-     !!!!!!!!!! Bivariate spline !!!!!!!!!!!
-
-     subroutine bispev(tx,nx,ty,ny,c,kx,ky,x,mx,y,my,z,wrk,lwrk,iwrk,kwrk,ier)
-       ! z,ier = bispev(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(mx) :: x
-       integer intent(hide),depend(x) :: mx=len(x)
-       real*8 intent(in),dimension(my) :: y
-       integer intent(hide),depend(y) :: my=len(y)
-       real*8 dimension(mx,my),depend(mx,my),intent(out,c) :: z
-       real*8 dimension(lwrk),depend(lwrk),intent(hide,cache) :: wrk
-       integer intent(hide),depend(mx,kx,my,ky) :: lwrk=mx*(kx+1)+my*(ky+1)
-       integer dimension(kwrk),depend(kwrk),intent(hide,cache) :: iwrk
-       integer intent(hide),depend(mx,my) :: kwrk=mx+my
-       integer intent(out) :: ier
-     end subroutine bispev
-
-     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)
-       !  nx,tx,ny,ty,c,fp,ier = surfit_smth(x,y,z,[w,xb,xe,yb,ye,kx,ky,s,eps,lwrk2])
-
-       fortranname surfit
-
-       integer intent(hide) :: iopt=0
-       integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
-            :: m=len(x)
-       real*8 dimension(m) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m) :: y
-       real*8 dimension(m),depend(m),check(len(z)==m) :: z
-       real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
-       real*8 optional,depend(x,m) :: xb=dmin(x,m)
-       real*8 optional,depend(x,m) :: xe=dmax(x,m)
-       real*8 optional,depend(y,m) :: yb=dmin(y,m)
-       real*8 optional,depend(y,m) :: ye=dmax(y,m)
-       integer check(1<=kx && kx<=5) :: kx = 3
-       integer check(1<=ky && ky<=5) :: ky = 3
-       real*8 optional,check(0.0<=s) :: s = m
-       integer optional,depend(kx,m),check(nxest>=2*(kx+1)) &
-            :: nxest = imax(kx+1+sqrt(m/2),2*(kx+1))
-       integer optional,depend(ky,m),check(nyest>=2*(ky+1)) &
-            :: nyest = imax(ky+1+sqrt(m/2),2*(ky+1))
-       integer intent(hide),depend(nxest,nyest) :: nmax=MAX(nxest,nyest)
-       real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
-       integer intent(out) :: nx
-       real*8 dimension(nmax),intent(out),depend(nmax) :: tx
-       integer intent(out) :: ny
-       real*8 dimension(nmax),intent(out),depend(nmax) :: ty
-       real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
-            depend(kx,ky,nxest,nyest),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk1),intent(cache,out),depend(lwrk1) :: wrk1
-       integer intent(hide),depend(m,kx,ky,nxest,nyest) &
-            :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
-       real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
-       integer optional,intent(in),depend(kx,ky,nxest,nyest) &
-            :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
-       integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
-       integer intent(hide),depend(m,nxest,nyest,kx,ky) &
-            :: kwrk=m+(nxest-2*kx-1)*(nyest-2*ky-1)
-       integer intent(out) :: ier
-     end subroutine surfit_smth
-
-     subroutine surfit_lsq(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)
-       ! tx,ty,c,fp,ier = surfit_lsq(x,y,z,tx,ty,[w,xb,xe,yb,ye,kx,ky,eps,lwrk2])
-
-       fortranname surfit
-
-       integer intent(hide) :: iopt=-1
-       integer intent(hide),depend(x,kx,ky),check(m>=(kx+1)*(ky+1)) &
-            :: m=len(x)
-       real*8 dimension(m) :: x
-       real*8 dimension(m),depend(m),check(len(y)==m) :: y
-       real*8 dimension(m),depend(m),check(len(z)==m) :: z
-       real*8 optional,dimension(m),depend(m),check(len(w)==m) :: w = 1.0
-       real*8 optional,depend(x,tx,m,nx) :: xb=calc_b(x,m,tx,nx)
-       real*8 optional,depend(x,tx,m,nx) :: xe=calc_e(x,m,tx,nx)
-       real*8 optional,depend(y,ty,m,ny) :: yb=calc_b(y,m,ty,ny)
-       real*8 optional,depend(y,ty,m,ny) :: ye=calc_e(y,m,ty,ny)
-       integer check(1<=kx && kx<=5) :: kx = 3
-       integer check(1<=ky && ky<=5) :: ky = 3
-       real*8 intent(hide) :: s = 0.0
-       integer intent(hide),depend(nx) :: nxest = nx
-       integer intent(hide),depend(ny) :: nyest = ny
-       integer intent(hide),depend(nx,ny) :: nmax=MAX(nx,ny)
-       real*8 optional,check(0.0<eps && eps<1.0) :: eps = 1e-16
-       integer intent(hide),depend(tx,kx),check(2*kx+2<=nx) :: nx = len(tx)
-       real*8 dimension(nx),intent(in,out,overwrite) :: tx
-       integer intent(hide),depend(ty,ky),check(2*ky+2<=ny) :: ny = len(ty)
-       real*8 dimension(ny),intent(in,out,overwrite) :: ty
-       real*8 dimension((nx-kx-1)*(ny-ky-1)),depend(kx,ky,nx,ny),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk1),intent(cache,hide),depend(lwrk1) :: wrk1
-       integer intent(hide),depend(m,kx,ky,nxest,nyest) &
-            :: lwrk1=calc_surfit_lwrk1(m,kx,ky,nxest,nyest)
-       real*8 dimension(lwrk2),intent(cache,hide),depend(lwrk2) :: wrk2
-       integer optional,intent(in),depend(kx,ky,nxest,nyest) &
-            :: lwrk2=calc_surfit_lwrk2(m,kx,ky,nxest,nyest)
-       integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
-       integer intent(hide),depend(m,nx,ny,kx,ky) &
-            :: 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])
-
-       fortranname regrid
-
-       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) 
-       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)
-       real*8 optional,depend(x,mx) :: xe=dmax(x,mx)
-       real*8 optional,depend(y,my) :: yb=dmin(y,my)
-       real*8 optional,depend(y,my) :: ye=dmax(y,my)
-       integer optional,check(1<=kx && kx<=5) :: kx = 3
-       integer optional,check(1<=ky && ky<=5) :: ky = 3
-       real*8 optional,check(0.0<=s) :: s = 0.0
-       integer intent(hide),depend(kx,mx),check(nxest>=2*(kx+1)) &
-            :: nxest = mx+kx+1
-       integer intent(hide),depend(ky,my),check(nyest>=2*(ky+1)) &
-            :: nyest = my+ky+1
-       integer intent(out) :: nx
-       real*8 dimension(nxest),intent(out),depend(nxest) :: tx
-       integer intent(out) :: ny
-       real*8 dimension(nyest),intent(out),depend(nyest) :: ty
-       real*8 dimension((nxest-kx-1)*(nyest-ky-1)), &
-            depend(kx,ky,nxest,nyest),intent(out) :: c
-       real*8 intent(out) :: fp
-       real*8 dimension(lwrk),intent(cache,hide),depend(lwrk) :: wrk
-       integer intent(hide),depend(mx,my,kx,ky,nxest,nyest) &
-            :: lwrk=calc_regrid_lwrk(mx,my,kx,ky,nxest,nyest)
-       integer dimension(kwrk),depend(kwrk),intent(cache,hide) :: iwrk
-       integer intent(hide),depend(mx,my,nxest,nyest) &
-            :: 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
-       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) :: xb
-       real*8 intent(in) :: xe
-       real*8 intent(in) :: yb
-       real*8 intent(in) :: ye
-       real*8 dimension(nx+ny-kx-ky-2),depend(nx,ny,kx,ky),intent(cache,hide) :: wrk
-       real*8 :: dblint
-     end function dblint
-  end interface
-end python module _dfitpack
-

Modified: branches/Interpolate1D/fitpack_wrapper.py
===================================================================
--- branches/Interpolate1D/fitpack_wrapper.py	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/fitpack_wrapper.py	2008-07-31 23:58:42 UTC (rev 4591)
@@ -1,8 +1,13 @@
 """
 This module is used for spline interpolation, and functions
 as a wrapper around the FITPACK Fortran interpolation
-package.
+package.  Its functionality is contained in the Spline
+class.
 
+Spline is primarily meant to be called by Interpolate1d
+or inter1d, but is a stand-alone class in its own right.
+
+
 The code has been modified from an older version of
 scipy.interpolate, where it was directly called by the
 user.  As such, it includes functionality not available through
@@ -29,10 +34,6 @@
 
     Can include least-squares fitting.
 
-    See also:
-
-    splrep, splev, sproot, spint, spalde - an older wrapping of FITPACK
-    BivariateSpline - a similar class for bivariate spline interpolation
     """
 
     def __init__(self, x=None, y=None, w=None, bbox = [None]*2, k=3, s=0.0):
@@ -43,16 +44,14 @@
 
         Optional input:
           w          - positive 1-d sequence of weights
-          bbox       - 2-sequence specifying the boundary of
+          bbox     - 2-sequence specifying the boundary of
                        the approximation interval.
                        By default, bbox=[x[0],x[-1]]
-          k=3        - degree of the univariate spline.
+          k=3     - degree of the univariate spline.
           s          - positive smoothing factor defined for
                        estimation condition:
                          sum((w[i]*( y[i]-s(x[i]) ))**2,axis=0) <= s
-                       Default s=len(w) which should be a good value
-                       if 1/w[i] is an estimate of the standard
-                       deviation of y[i].
+                        Default s=0
         """
         
         self._k = k

Modified: branches/Interpolate1D/interpolate1d.py
===================================================================
--- branches/Interpolate1D/interpolate1d.py	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/interpolate1d.py	2008-07-31 23:58:42 UTC (rev 4591)
@@ -1,18 +1,22 @@
 # FIXME: information strings giving mathematical descriptions of the actions
 #     of the functions.
 
-from interpolate_wrapper import linear, logarithmic, block, \
-                block_average_above, atleast_1d_and_contiguous
+from interpolate_wrapper import atleast_1d_and_contiguous, \
+                linear, logarithmic, block, block_average_above, nearest
 from fitpack_wrapper import Spline
 import numpy as np
 from numpy import array, arange, empty, float64, NaN
 
 # dictionary of interpolation functions/classes/objects
 method_register = \
-                { 'linear' : linear, 
+                { # functions
+                    'linear' : linear, 
                     'logarithmic' : logarithmic, 
                     'block' : block, 
                     'block_average_above' : block_average_above, 
+                    'nearest' : nearest,
+                    
+                    # Splines
                     'Spline' : Spline, 'spline' : Spline,
                     'Quadratic' : Spline(k=2), 'quadratic' : Spline(k=2),
                     'Quad' : Spline(k=2), 'quad' : Spline(k=2),

Modified: branches/Interpolate1D/interpolate_wrapper.py
===================================================================
--- branches/Interpolate1D/interpolate_wrapper.py	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/interpolate_wrapper.py	2008-07-31 23:58:42 UTC (rev 4591)
@@ -6,9 +6,27 @@
 import sys
 import _interpolate # C extension.  Does all the real work.
 
-def atleast_1d_and_contiguous(ary, typecode = np.float64):
-    return np.atleast_1d( np.ascontiguousarray(ary, typecode) )
+def atleast_1d_and_contiguous(ary, dtype = np.float64):
+    return np.atleast_1d( np.ascontiguousarray(ary, dtype) )
 
+def nearest(x, y, new_x):
+    """ Rounds each new_x[i] to the closest value in x
+        and returns corresponding y.
+    """
+    shifted_x = np.concatenate(( np.array([x[0]-1]) , x[0:-1] ))
+    
+    midpoints_of_x = atleast_1d_and_contiguous( .5*(x + shifted_x) )
+    new_x = atleast_1d_and_contiguous(new_x)
+    
+    TINY = 1e-10
+    indices = np.searchsorted(midpoints_of_x, new_x+TINY)-1
+    indices = np.atleast_1d(np.clip(indices, 0, np.Inf).astype(np.int))
+    new_y = np.take(y, indices, axis=-1)
+    
+    return new_y
+    
+    
+
 def linear(x, y, new_x):
     """ Linearly interpolates values in new_x based on the values in x and y
 
@@ -106,7 +124,6 @@
             For each new_x[i], finds largest j such that
             x[j] < new_x[j], and returns y[j].
         """
-
         # find index of values in x that preceed values in x
         # This code is a little strange -- we really want a routine that
         # returns the index of values where x[j] < x[index]

Modified: branches/Interpolate1D/setup.py
===================================================================
--- branches/Interpolate1D/setup.py	2008-07-31 23:46:55 UTC (rev 4590)
+++ branches/Interpolate1D/setup.py	2008-07-31 23:58:42 UTC (rev 4591)
@@ -22,7 +22,7 @@
 
     # Fortran routines (collectively "FITPACK" for spline interpolation)
     config.add_extension('_dfitpack',
-                         sources=['fitpack.pyf'],
+                         sources=['_fitpack.pyf'],
                          libraries=['_fitpack'],
                         )
                         



More information about the Scipy-svn mailing list