[Scipy-svn] r2254 - in trunk/Lib/interpolate: . tests

scipy-svn at scipy.org scipy-svn at scipy.org
Wed Oct 11 18:29:20 CDT 2006


Author: stefan
Date: 2006-10-11 18:28:59 -0500 (Wed, 11 Oct 2006)
New Revision: 2254

Modified:
   trunk/Lib/interpolate/fitpack.pyf
   trunk/Lib/interpolate/fitpack2.py
   trunk/Lib/interpolate/tests/test_fitpack.py
Log:
Add RectBivariateSpline [for John Travers].  Fix long lines.


Modified: trunk/Lib/interpolate/fitpack.pyf
===================================================================
--- trunk/Lib/interpolate/fitpack.pyf	2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/fitpack.pyf	2006-10-11 23:28:59 UTC (rev 2254)
@@ -63,6 +63,12 @@
  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
@@ -411,7 +417,46 @@
             :: 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
+
   end interface
 end python module dfitpack
 

Modified: trunk/Lib/interpolate/fitpack2.py
===================================================================
--- trunk/Lib/interpolate/fitpack2.py	2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/fitpack2.py	2006-10-11 23:28:59 UTC (rev 2254)
@@ -13,10 +13,11 @@
     'LSQUnivariateSpline',
 
     'LSQBivariateSpline',
-    'SmoothBivariateSpline']
+    'SmoothBivariateSpline',
+    'RectBivariateSpline']
 
 import warnings
-from numpy import zeros, concatenate, alltrue
+from numpy import zeros, concatenate, alltrue, ravel, all, diff
 
 import dfitpack
 
@@ -134,6 +135,7 @@
     def set_smoothing_factor(self, s):
         """ Continue spline computation with the given smoothing
         factor s and with the knots found at the last call.
+        
         """
         data = self._data
         if data[6]==-1:
@@ -150,8 +152,9 @@
 
     def __call__(self, x, nu=None):
         """ Evaluate spline (or its nu-th derivative) at positions x.
-        Note: x can be unordered but the evaluation is
-        more efficient if x is (partially) ordered.
+        Note: x can be unordered but the evaluation is more efficient
+        if x is (partially) ordered.
+        
         """
         if nu is None:
             return dfitpack.splev(*(self._eval_args+(x,)))
@@ -166,14 +169,15 @@
         return data[8][k:n-k]
 
     def get_coeffs(self):
-        """ Return spline coefficients."""
+        """Return spline coefficients."""
         data = self._data
         k,n = data[5],data[7]
         return data[9][:n-k-1]
 
     def get_residual(self):
-        """ Return weighted sum of squared residuals of the spline
+        """Return weighted sum of squared residuals of the spline
         approximation: sum ((w[i]*(y[i]-s(x[i])))**2,axis=0)
+        
         """
         return self._data[10]
 
@@ -199,11 +203,15 @@
             z,m,ier = dfitpack.sproot(*self._eval_args[:2])
             assert ier==0,`ier`
             return z[:m]
-        raise NotImplementedError,'finding roots unsupported for non-cubic splines'
+        raise NotImplementedError,\
+              'finding roots unsupported for non-cubic splines'
 
 class InterpolatedUnivariateSpline(UnivariateSpline):
-    """ Interpolated univariate spline approximation. Identical to UnivariateSpline with less error checking."""
+    """ Interpolated univariate spline approximation. Identical to
+    UnivariateSpline with less error checking.
 
+    """
+
     def __init__(self, x, y, w=None, bbox = [None]*2, k=3):
         """
         Input:
@@ -223,8 +231,12 @@
         self._reset_class()
 
 class LSQUnivariateSpline(UnivariateSpline):
-    """ Weighted least-squares univariate spline approximation. Appears to be identical to UnivariateSpline with more error checking."""
+    """ Weighted least-squares univariate spline
+    approximation. Appears to be identical to UnivariateSpline with
+    more error checking.
 
+    """
+
     def __init__(self, x, y, t, w=None, bbox = [None]*2, k=3):
         """
         Input:
@@ -361,7 +373,8 @@
           w          - positive 1-d sequence of weights
           bbox       - 4-sequence specifying the boundary of
                        the rectangular approximation domain.
-                       By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
+                       By default, bbox=[min(x,tx),max(x,tx),
+                                         min(y,ty),max(y,ty)]
           kx,ky=3,3  - degrees of the bivariate spline.
           s          - positive smoothing factor defined for
                        estimation condition:
@@ -374,9 +387,10 @@
                        equations. 0 < eps < 1, default is 1e-16.
         """
         xb,xe,yb,ye = bbox
-        nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth(x,y,z,w,\
-                                                         xb,xe,yb,ye,\
-                                                         kx,ky,s=s,eps=eps,lwrk2=1)
+        nx,tx,ny,ty,c,fp,wrk1,ier = dfitpack.surfit_smth(x,y,z,w,
+                                                         xb,xe,yb,ye,
+                                                         kx,ky,s=s,
+                                                         eps=eps,lwrk2=1)
         if ier in [0,-1,-2]: # normal return
             pass
         else:
@@ -410,7 +424,8 @@
           w          - positive 1-d sequence of weights
           bbox       - 4-sequence specifying the boundary of
                        the rectangular approximation domain.
-                       By default, bbox=[min(x,tx),max(x,tx),min(y,ty),max(y,ty)]
+                       By default, bbox=[min(x,tx),max(x,tx),
+                                         min(y,ty),max(y,ty)]
           kx,ky=3,3  - degrees of the bivariate spline.
           eps        - a threshold for determining the effective rank
                        of an over-determined linear system of
@@ -443,3 +458,63 @@
         self.fp = fp
         self.tck = tx1,ty1,c
         self.degrees = kx,ky
+
+class RectBivariateSpline(BivariateSpline):
+    """ Bivariate spline approximation over a rectangular mesh.
+
+    Can be used for both smoothing or interpolating data.
+
+    See also:
+
+    SmoothBivariateSpline - a smoothing bivariate spline for scattered data
+    bisplrep, bisplev - an older wrapping of FITPACK
+    UnivariateSpline - a similar class for univariate spline interpolation
+    """
+
+    def __init__(self, x, y, z,
+                 bbox = [None]*4, kx=3, ky=3, s=0):
+        """
+        Input:
+          x,y  - 1-d sequences of coordinates in strictly ascending order
+            z  - 2-d array of data with shape (x.size,y.size)
+        Optional input:
+          bbox       - 4-sequence specifying the boundary of
+                       the rectangular approximation domain.
+                       By default, bbox=[min(x,tx),max(x,tx),
+                                         min(y,ty),max(y,ty)]
+          kx,ky=3,3  - degrees of the bivariate spline.
+          s          - positive smoothing factor defined for
+                       estimation condition:
+                         sum((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0) <= s
+                       Default s=0 which is for interpolation
+        """
+        x,y = ravel(x),ravel(y)
+        if not all(diff(x) > 0.0):
+            raise TypeError,'x must be strictly increasing'
+        if not all(diff(y) > 0.0):
+            raise TypeError,'y must be strictly increasing'
+        if not ((x.min() == x[0]) and (x.max() == x[-1])):
+            raise TypeError, 'x must be strictly ascending' 
+        if not ((y.min() == y[0]) and (y.max() == y[-1])):
+            raise TypeError, 'y must be strictly ascending' 
+        if not x.size == z.shape[0]:
+            raise TypeError,\
+                  'x dimension of z must have same number of elements as x'
+        if not y.size == z.shape[1]:
+            raise TypeError,\
+                  'y dimension of z must have same number of elements as y'
+        z = ravel(z)
+        xb,xe,yb,ye = bbox
+        nx,tx,ny,ty,c,fp,ier = dfitpack.regrid_smth(x,y,z,
+                                                    xb,xe,yb,ye,
+                                                    kx,ky,s)
+        if ier in [0,-1,-2]: # normal return
+            pass
+        else:
+            message = _surfit_messages.get(ier,'ier=%s' % (ier))
+            warnings.warn(message)
+
+        self.fp = fp
+        self.tck = tx[:nx],ty[:ny],c[:(nx-kx-1)*(ny-ky-1)]
+        self.degrees = kx,ky
+

Modified: trunk/Lib/interpolate/tests/test_fitpack.py
===================================================================
--- trunk/Lib/interpolate/tests/test_fitpack.py	2006-10-11 20:44:55 UTC (rev 2253)
+++ trunk/Lib/interpolate/tests/test_fitpack.py	2006-10-11 23:28:59 UTC (rev 2254)
@@ -14,10 +14,12 @@
 
 import sys
 from numpy.testing import *
+from numpy import array
 set_package_path()
 from interpolate.fitpack2 import UnivariateSpline,LSQUnivariateSpline,\
      InterpolatedUnivariateSpline
-from interpolate.fitpack2 import LSQBivariateSpline, SmoothBivariateSpline
+from interpolate.fitpack2 import LSQBivariateSpline, SmoothBivariateSpline,\
+     RectBivariateSpline
 restore_path()
 
 class test_UnivariateSpline(ScipyTestCase):
@@ -73,5 +75,13 @@
         assert_almost_equal(lut.get_residual(),0.0)
         assert_array_almost_equal(lut([1,1.5,2],[1,1.5]),[[0,0],[1,1],[2,2]])
 
+class test_RectBivariateSpline(ScipyTestCase):
+    def check_defaults(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)
+        assert_array_almost_equal(lut(x,y),z)
+
 if __name__ == "__main__":
     ScipyTest().run()



More information about the Scipy-svn mailing list