[Scipy-svn] r2290 - in trunk/Lib/sandbox: . spline

scipy-svn at scipy.org scipy-svn at scipy.org
Wed Oct 18 14:28:09 CDT 2006


Author: jtravs
Date: 2006-10-18 14:27:50 -0500 (Wed, 18 Oct 2006)
New Revision: 2290

Added:
   trunk/Lib/sandbox/spline/
   trunk/Lib/sandbox/spline/README.txt
   trunk/Lib/sandbox/spline/__init__.py
   trunk/Lib/sandbox/spline/fitpack.py
   trunk/Lib/sandbox/spline/fitpack.pyf
   trunk/Lib/sandbox/spline/fitpack/
   trunk/Lib/sandbox/spline/info.py
   trunk/Lib/sandbox/spline/setup.py
   trunk/Lib/sandbox/spline/spline.py
   trunk/Lib/sandbox/spline/tests/
Modified:
   trunk/Lib/sandbox/setup.py
Log:
First check in of spline module (based on interpolate)


Modified: trunk/Lib/sandbox/setup.py
===================================================================
--- trunk/Lib/sandbox/setup.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/setup.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -76,6 +76,9 @@
     # Package for Gaussian Mixture Models
     #config.add_subpackage('pyem')
 
+    # New spline package (based on scipy.interpolate)
+    #config.add_subpackage('spline')
+
     return config
 
 if __name__ == '__main__':

Added: trunk/Lib/sandbox/spline/README.txt
===================================================================
--- trunk/Lib/sandbox/spline/README.txt	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/README.txt	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,16 @@
+This module is meant to replace the spline functionality of scipy.interpolate.
+At the moment the code base is identical with a few minor cosmetic changes.
+It does not compile yet!
+
+Planned changes are:
+
+1. remove dependence of f2c generated interface (use only f2py)
+2. cleanup!
+3. rationalise the interface (particularly the new spline/fitpack2 interface)
+4. build comprehensive unit test suite
+
+Hopefully this module will end up in scipy, with scipy.interpolate only then
+containing the interpolation functions interp1d, interp2d etc. which may depend
+on the spline module and the new delaunay module (for scattered data).
+
+John Travers

Copied: trunk/Lib/sandbox/spline/__init__.py (from rev 2289, trunk/Lib/interpolate/__init__.py)
===================================================================
--- trunk/Lib/interpolate/__init__.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/__init__.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,14 @@
+#
+# spline - Spline Tools
+#
+
+from info import __doc__
+
+from fitpack import *
+
+# New interface to fitpack library:
+from spline import *
+
+__all__ = filter(lambda s:not s.startswith('_'),dir())
+from numpy.testing import ScipyTest
+test = ScipyTest().test

Copied: trunk/Lib/sandbox/spline/fitpack (from rev 2289, trunk/Lib/interpolate/fitpack)

Copied: trunk/Lib/sandbox/spline/fitpack.py (from rev 2289, trunk/Lib/interpolate/fitpack.py)
===================================================================
--- trunk/Lib/interpolate/fitpack.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/fitpack.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,927 @@
+#!/usr/bin/env python
+"""
+fitpack (dierckx in netlib) --- A Python-C wrapper to FITPACK (by P. Dierckx).
+        FITPACK is a collection of FORTRAN programs for curve and surface
+        fitting with splines and tensor product splines.
+
+See
+ http://www.cs.kuleuven.ac.be/cwis/research/nalag/research/topics/fitpack.html
+or
+ http://www.netlib.org/dierckx/index.html
+
+Copyright 2002 Pearu Peterson all rights reserved,
+Pearu Peterson <pearu at cens.ioc.ee>
+Permission to use, modify, and distribute this software is given under the
+terms of the SciPy (BSD style) license.  See LICENSE.txt that came with
+this distribution for specifics.
+
+NO WARRANTY IS EXPRESSED OR IMPLIED.  USE AT YOUR OWN RISK.
+
+Pearu Peterson
+Modified by John Travers <jtravs at gmail.com>
+
+Running test programs:
+    $ python fitpack.py 1 3    # run test programs 1, and 3
+    $ python fitpack.py        # run all available test programs
+
+TODO: Make interfaces to the following fitpack functions:
+    For univariate splines: cocosp, concon, fourco, insert
+    For bivariate splines: profil, regrid, parsur, surev
+"""
+
+__all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
+    'bisplrep', 'bisplev']
+__version__ = "$Revision$"[10:-1]
+import _fitpack
+from numpy import atleast_1d, array, ones, zeros, sqrt, ravel, transpose, \
+     dot, sin, cos, pi, arange, empty, int32
+myasarray = atleast_1d
+
+# Try to replace _fitpack interface with
+#  f2py-generated version
+import dfitpack
+
+_iermess = {0:["""\
+    The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
+               -1:["""\
+    The spline is an interpolating spline (fp=0)""",None],
+               -2:["""\
+    The spline is weighted least-squares polynomial of degree k.
+    fp gives the upper bound fp0 for the smoothing factor s""",None],
+               1:["""\
+    The required storage space exceeds the available storage space.
+    Probable causes: data (x,y) size is too small or smoothing parameter s is too small (fp>s).""",ValueError],
+               2:["""\
+    A theoretically impossible results when finding a smoothin spline
+    with fp = s. Probably causes: s too small. (abs(fp-s)/s>0.001)""",ValueError],
+               3:["""\
+    The maximal number of iterations (20) allowed for finding smoothing
+    spline with fp=s has been reached. Probably causes: s too small.
+    (abs(fp-s)/s>0.001)""",ValueError],
+               10:["""\
+    Error on input data""",ValueError],
+               'unknown':["""\
+    An error occured""",TypeError]}
+
+_iermess2 = {0:["""\
+    The spline has a residual sum of squares fp such that abs(fp-s)/s<=0.001""",None],
+            -1:["""\
+    The spline is an interpolating spline (fp=0)""",None],
+            -2:["""\
+    The spline is weighted least-squares polynomial of degree kx and ky.
+    fp gives the upper bound fp0 for the smoothing factor s""",None],
+            -3:["""\
+    Warning. The coefficients of the spline have been computed as the minimal
+    norm least-squares solution of a rank deficient system.""",None],
+            1:["""\
+    The required storage space exceeds the available storage space.
+    Probably causes: nxest or nyest too small or s is too small. (fp>s)""",ValueError],
+            2:["""\
+    A theoretically impossible results when finding a smoothin spline
+    with fp = s. Probably causes: s too small or badly chosen eps.
+    (abs(fp-s)/s>0.001)""",ValueError],
+            3:["""\
+    The maximal number of iterations (20) allowed for finding smoothing
+    spline with fp=s has been reached. Probably causes: s too small.
+    (abs(fp-s)/s>0.001)""",ValueError],
+            4:["""\
+    No more knots can be added because the number of B-spline coefficients
+    already exceeds the number of data points m. Probably causes: either
+    s or m too small. (fp>s)""",ValueError],
+            5:["""\
+    No more knots can be added because the additional knot would coincide
+    with an old one. Probably cause: s too small or too large a weight
+    to an inaccurate data point. (fp>s)""",ValueError],
+            10:["""\
+    Error on input data""",ValueError],
+            11:["""\
+    rwrk2 too small, i.e. there is not enough workspace for computing
+    the minimal least-squares solution of a rank deficient system of linear
+    equations.""",ValueError],
+            'unknown':["""\
+    An error occured""",TypeError]}
+
+_parcur_cache = {'t': array([],float), 'wrk': array([],float),
+                 'iwrk':array([],int32), 'u': array([],float),'ub':0,'ue':1}
+
+def splprep(x,w=None,u=None,ub=None,ue=None,k=3,task=0,s=None,t=None,
+            full_output=0,nest=None,per=0,quiet=1):
+    """Find the B-spline representation of an N-dimensional curve.
+
+    Description:
+
+      Given a list of N rank-1 arrays, x, which represent a curve in N-dimensional
+      space parametrized by u, find a smooth approximating spline curve g(u).
+      Uses the FORTRAN routine parcur from FITPACK
+
+    Inputs:
+
+      x -- A list of sample vector arrays representing the curve.
+      u -- An array of parameter values.  If not given, these values are
+           calculated automatically as (M = len(x[0])):
+           v[0] = 0
+           v[i] = v[i-1] + distance(x[i],x[i-1])
+           u[i] = v[i] / v[M-1]
+      ub, ue -- The end-points of the parameters interval.  Defaults to
+                u[0] and u[-1].
+      k -- Degree of the spline.  Cubic splines are recommended.  Even values of
+           k should be avoided especially with a small s-value.
+           1 <= k <= 5.
+      task -- If task==0 find t and c for a given smoothing factor, s.
+              If task==1 find t and c for another value of the smoothing factor,
+                s. There must have been a previous call with task=0 or task=1
+                for the same set of data.
+              If task=-1 find the weighted least square spline for a given set of
+                knots, t.
+      s -- A smoothing condition.  The amount of smoothness is determined by
+           satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
+           g(x) is the smoothed interpolation of (x,y).  The user can use s to
+           control the tradeoff between closeness and smoothness of fit.  Larger
+           s means more smoothing while smaller values of s indicate less
+           smoothing. Recommended values of s depend on the weights, w.  If the
+           weights represent the inverse of the standard-deviation of y, then a
+           good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
+           where m is the number of datapoints in x, y, and w.
+      t -- The knots needed for task=-1.
+      full_output -- If non-zero, then return optional outputs.
+      nest -- An over-estimate of the total number of knots of the spline to
+              help in determining the storage space.  By default nest=m/2.
+              Always large enough is nest=m+k+1.
+      per -- If non-zero, data points are considered periodic with period
+             x[m-1] - x[0] and a smooth periodic spline approximation is returned.
+             Values of y[m-1] and w[m-1] are not used.
+      quiet -- Non-zero to suppress messages.
+
+    Outputs: (tck, u, {fp, ier, msg})
+
+      tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
+             coefficients, and the degree of the spline.
+      u -- An array of the values of the parameter.
+
+      fp -- The weighted sum of squared residuals of the spline approximation.
+      ier -- An integer flag about splrep success.  Success is indicated
+             if ier<=0. If ier in [1,2,3] an error occurred but was not raised.
+             Otherwise an error is raised.
+      msg -- A message corresponding to the integer flag, ier.
+
+    Remarks:
+
+      SEE splev for evaluation of the spline and its derivatives.
+
+    See also:
+      splrep, splev, sproot, spalde, splint - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    if task<=0:
+        _parcur_cache = {'t': array([],float), 'wrk': array([],float),
+                         'iwrk':array([],int32),'u': array([],float),
+                         'ub':0,'ue':1}
+    x=myasarray(x)
+    idim,m=x.shape
+    if per:
+        for i in range(idim):
+            if x[i][0]!=x[i][-1]:
+                if quiet<2:print 'Warning: Setting x[%d][%d]=x[%d][0]'%(i,m,i)
+                x[i][-1]=x[i][0]
+    if not 0<idim<11: raise TypeError,'0<idim<11 must hold'
+    if w is None: w=ones(m,float)
+    else: w=myasarray(w)
+    ipar=(u is not None)
+    if ipar:
+        _parcur_cache['u']=u
+        if ub is None: _parcur_cache['ub']=u[0]
+        else: _parcur_cache['ub']=ub
+        if ue is None: _parcur_cache['ue']=u[-1]
+        else: _parcur_cache['ue']=ue
+    else: _parcur_cache['u']=zeros(m,float)
+    if not (1<=k<=5): raise TypeError, '1<=k=%d<=5 must hold'%(k)
+    if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
+    if (not len(w)==m) or (ipar==1 and (not len(u)==m)):
+        raise TypeError,'Mismatch of input dimensions'
+    if s is None: s=m-sqrt(2*m)
+    if t is None and task==-1: raise TypeError, 'Knots must be given for task=-1'
+    if t is not None:
+        _parcur_cache['t']=myasarray(t)
+    n=len(_parcur_cache['t'])
+    if task==-1 and n<2*k+2:
+        raise TypeError, 'There must be at least 2*k+2 knots for task=-1'
+    if m<=k: raise TypeError, 'm>k must hold'
+    if nest is None: nest=m+2*k
+
+    if (task>=0 and s==0) or (nest<0):
+        if per: nest=m+2*k
+        else: nest=m+k+1
+    nest=max(nest,2*k+3)
+    u=_parcur_cache['u']
+    ub=_parcur_cache['ub']
+    ue=_parcur_cache['ue']
+    t=_parcur_cache['t']
+    wrk=_parcur_cache['wrk']
+    iwrk=_parcur_cache['iwrk']
+    t,c,o=_fitpack._parcur(ravel(transpose(x)),w,u,ub,ue,k,task,ipar,s,t,
+                             nest,wrk,iwrk,per)
+    _parcur_cache['u']=o['u']
+    _parcur_cache['ub']=o['ub']
+    _parcur_cache['ue']=o['ue']
+    _parcur_cache['t']=t
+    _parcur_cache['wrk']=o['wrk']
+    _parcur_cache['iwrk']=o['iwrk']
+    ier,fp,n=o['ier'],o['fp'],len(t)
+    u=o['u']
+    c.shape=idim,n-k-1
+    tcku = [t,list(c),k],u
+    if ier<=0 and not quiet:
+        print _iermess[ier][0]
+        print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
+    if ier>0 and not full_output:
+        if ier in [1,2,3]:
+            print "Warning: "+_iermess[ier][0]
+        else:
+            try:
+                raise _iermess[ier][1],_iermess[ier][0]
+            except KeyError:
+                raise _iermess['unknown'][1],_iermess['unknown'][0]
+    if full_output:
+        try:
+            return tcku,fp,ier,_iermess[ier][0]
+        except KeyError:
+            return tcku,fp,ier,_iermess['unknown'][0]
+    else:
+        return tcku
+
+_curfit_cache = {'t': array([],float), 'wrk': array([],float),
+                 'iwrk':array([],int32)}
+def splrep(x,y,w=None,xb=None,xe=None,k=3,task=0,s=1e-3,t=None,
+           full_output=0,per=0,quiet=1):
+    """Find the B-spline representation of 1-D curve.
+
+    Description:
+
+      Given the set of data points (x[i], y[i]) determine a smooth spline
+      approximation of degree k on the interval xb <= x <= xe.  The coefficients,
+      c, and the knot points, t, are returned.  Uses the FORTRAN routine
+      curfit from FITPACK.
+
+    Inputs:
+
+      x, y -- The data points defining a curve y = f(x).
+      w -- Strictly positive rank-1 array of weights the same length as x and y.
+           The weights are used in computing the weighted least-squares spline
+           fit. If the errors in the y values have standard-deviation given by the
+           vector d, then w should be 1/d. Default is ones(len(x)).
+      xb, xe -- The interval to fit.  If None, these default to x[0] and x[-1]
+                respectively.
+      k -- The order of the spline fit.  It is recommended to use cubic splines.
+           Even order splines should be avoided especially with small s values.
+           1 <= k <= 5
+      task -- If task==0 find t and c for a given smoothing factor, s.
+              If task==1 find t and c for another value of the
+                smoothing factor, s. There must have been a previous
+                call with task=0 or task=1 for the same set of data
+                (t will be stored an used internally)
+              If task=-1 find the weighted least square spline for
+                a given set of knots, t.  These should be interior knots
+                as knots on the ends will be added automatically.
+      s -- A smoothing condition.  The amount of smoothness is determined by
+           satisfying the conditions: sum((w * (y - g))**2,axis=0) <= s where
+           g(x) is the smoothed interpolation of (x,y).  The user can use s to
+           control the tradeoff between closeness and smoothness of fit.  Larger
+           s means more smoothing while smaller values of s indicate less
+           smoothing. Recommended values of s depend on the weights, w.  If the
+           weights represent the inverse of the standard-deviation of y, then a
+           good s value should be found in the range (m-sqrt(2*m),m+sqrt(2*m))
+           where m is the number of datapoints in x, y, and w.
+           default : s=m-sqrt(2*m)
+      t -- The knots needed for task=-1.  If given then task is automatically
+           set to -1.
+      full_output -- If non-zero, then return optional outputs.
+      per -- If non-zero, data points are considered periodic with period
+             x[m-1] - x[0] and a smooth periodic spline approximation is returned.
+             Values of y[m-1] and w[m-1] are not used.
+      quiet -- Non-zero to suppress messages.
+
+    Outputs: (tck, {fp, ier, msg})
+
+      tck -- (t,c,k) a tuple containing the vector of knots, the B-spline
+             coefficients, and the degree of the spline.
+
+      fp -- The weighted sum of squared residuals of the spline approximation.
+      ier -- An integer flag about splrep success.  Success is indicated if
+             ier<=0. If ier in [1,2,3] an error occurred but was not raised.
+             Otherwise an error is raised.
+      msg -- A message corresponding to the integer flag, ier.
+
+    Remarks:
+
+      See splev for evaluation of the spline and its derivatives.
+      
+    Example:
+        
+      x = linspace(0, 10, 10)
+      y = sin(x)
+      tck = splrep(x, y)
+      x2 = linspace(0, 10, 200)
+      y2 = splev(x2, tck)
+      plot(x, y, 'o', x2, y2)
+      
+    See also:
+      splprep, splev, sproot, spalde, splint - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    if task<=0:
+        _curfit_cache = {}
+    x,y=map(myasarray,[x,y])
+    m=len(x)
+    if w is None: w=ones(m,float)
+    else: w=myasarray(w)
+    if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
+    if (m != len(y)) or (m != len(w)):
+        raise TypeError, 'Lengths of the first three arguments (x,y,w) must be equal'
+    if not (1<=k<=5):
+        raise TypeError, 'Given degree of the spline (k=%d) is not supported. (1<=k<=5)'%(k)
+    if m<=k: raise TypeError, 'm>k must hold'     
+    if xb is None: xb=x[0]
+    if xe is None: xe=x[-1]
+    if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
+    if s is None: s = m-sqrt(2*m)
+    if t is not None:
+        task = -1
+    if task == -1:
+        if t is None: raise TypeError, 'Knots must be given for task=-1'
+        numknots = len(t)
+        _curfit_cache['t'] = empty((numknots + 2*k+2,),float)
+        _curfit_cache['t'][k+1:-k-1] = t
+        nest = len(_curfit_cache['t'])
+    elif task == 0:
+        if per:
+            nest = max(m+2*k,2*k+3)
+        else:
+            nest = max(m+k+1,2*k+3)
+        t = empty((nest,),float)
+        _curfit_cache['t'] = t
+    if task <= 0:
+        _curfit_cache['wrk'] = empty((m*(k+1)+nest*(7+3*k),),float)
+        _curfit_cache['iwrk'] = empty((nest,),int32)
+    try:
+        t=_curfit_cache['t']
+        wrk=_curfit_cache['wrk']
+        iwrk=_curfit_cache['iwrk']
+    except KeyError:
+        raise TypeError, "must call with task=1 only after"\
+              " call with task=0,-1"
+    if not per:
+        n,c,fp,ier = dfitpack.curfit(task, x, y, w, t, wrk, iwrk, xb, xe, k, s)
+    else:
+        n,c,fp,ier = dfitpack.percur(task, x, y, w, t, wrk, iwrk, k, s)
+    tck = (t[:n],c[:n-k-1],k)
+    if ier<=0 and not quiet:
+        print _iermess[ier][0]
+        print "\tk=%d n=%d m=%d fp=%f s=%f"%(k,len(t),m,fp,s)
+    if ier>0 and not full_output:
+        if ier in [1,2,3]:
+            print "Warning: "+_iermess[ier][0]
+        else:
+            try:
+                raise _iermess[ier][1],_iermess[ier][0]
+            except KeyError:
+                raise _iermess['unknown'][1],_iermess['unknown'][0]
+    if full_output:
+        try:
+            return tck,fp,ier,_iermess[ier][0]
+        except KeyError:
+            return tck,fp,ier,_iermess['unknown'][0]
+    else:
+        return tck
+
+def _ntlist(l): # return non-trivial list
+    return l
+    #if len(l)>1: return l
+    #return l[0]
+
+def splev(x,tck,der=0):
+    """Evaulate a B-spline and its derivatives.
+
+    Description:
+
+      Given the knots and coefficients of a B-spline representation, evaluate
+      the value of the smoothing polynomial and it's derivatives.
+      This is a wrapper around the FORTRAN routines splev and splder of FITPACK.
+
+    Inputs:
+
+      x (u) -- a 1-D array of points at which to return the value of the
+               smoothed spline or its derivatives.  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.
+      der -- The order of derivative of the spline to compute (must be less than
+             or equal to k).
+
+    Outputs: (y, )
+
+      y -- an array of values representing the spline function or curve.
+           If tck was returned from splrep, then this is a list of arrays
+           representing the curve in N-dimensional space.
+
+    See also:
+      splprep, splrep, sproot, spalde, splint - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    t,c,k=tck
+    try:
+        c[0][0]
+        return map(lambda c,x=x,t=t,k=k,der=der:splev(x,[t,c,k],der),c)
+    except: pass
+    if not (0<=der<=k):
+        raise ValueError,"0<=der=%d<=k=%d must hold"%(der,k)
+    x=myasarray(x)
+    y,ier=_fitpack._spl_(x,der,t,c,k)
+    if ier==10: raise ValueError,"Invalid input data"
+    if ier: raise TypeError,"An error occurred"
+    if len(y)>1: return y
+    return y[0]
+
+def splint(a,b,tck,full_output=0):
+    """Evaluate the definite integral of a B-spline.
+
+    Description:
+
+      Given the knots and coefficients of a B-spline, evaluate the definite
+      integral of the smoothing polynomial between two given points.
+
+    Inputs:
+
+      a, b -- The end-points of the integration interval.
+      tck -- A length 3 sequence describing the given spline (See splev).
+      full_output -- Non-zero to return optional output.
+
+    Outputs: (integral, {wrk})
+
+      integral -- The resulting integral.
+      wrk -- An array containing the integrals of the normalized B-splines defined
+             on the set of knots.
+
+
+    See also:
+      splprep, splrep, sproot, spalde, splev - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    t,c,k=tck
+    try: c[0][0];return _ntlist(map(lambda c,a=a,b=b,t=t,k=k:splint(a,b,[t,c,k]),c))
+    except: pass
+    aint,wrk=_fitpack._splint(t,c,k,a,b)
+    if full_output: return aint,wrk
+    else: return aint
+
+def sproot(tck,mest=10):
+    """Find the roots of a cubic B-spline.
+
+    Description:
+
+      Given the knots (>=8) and coefficients of a cubic B-spline return the
+      roots of the spline.
+
+    Inputs:
+
+      tck -- A length 3 sequence describing the given spline (See splev).
+             The number of knots must be >= 8.  The knots must be a montonically
+             increasing sequence.
+      mest -- An estimate of the number of zeros (Default is 10).
+
+    Outputs: (zeros, )
+
+      zeros -- An array giving the roots of the spline.
+
+    See also:
+      splprep, splrep, splint, spalde, splev - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    t,c,k=tck
+    if k==4: t=t[1:-1]
+    if k==5: t=t[2:-2]
+    try: c[0][0];return _ntlist(map(lambda c,t=t,k=k,mest=mest:sproot([t,c,k],mest),c))
+    except: pass
+    if len(t)<8:
+        raise TypeError,"The number of knots %d>=8"%(len(t))
+    z,ier=_fitpack._sproot(t,c,k,mest)
+    if ier==10:
+        raise TypeError,"Invalid input data. t1<=..<=t4<t5<..<tn-3<=..<=tn must hold."
+    if ier==0: return z
+    if ier==1:
+        print "Warning: the number of zeros exceeds mest"
+        return z
+    raise TypeError,"Unknown error"
+
+def spalde(x,tck):
+    """Evaluate all derivatives of a B-spline.
+
+    Description:
+
+      Given the knots and coefficients of a cubic B-spline compute all
+      derivatives up to order k at a point (or set of points).
+
+    Inputs:
+
+      tck -- A length 3 sequence describing the given spline (See splev).
+      x -- A point or a set of points at which to evaluate the derivatives.
+           Note that t(k) <= x <= t(n-k+1) must hold for each x.
+
+    Outputs: (results, )
+
+      results -- An array (or a list of arrays) containing all derivatives
+                 up to order k inclusive for each point x.
+
+    See also:
+      splprep, splrep, splint, sproot, splev - evaluation, roots, integral
+      bisplrep, bisplev - bivariate splines
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    t,c,k=tck
+    try:
+        c[0][0]
+        return _ntlist(map(lambda c,x=x,t=t,k=k:spalde(x,[t,c,k]),c))
+    except: pass
+    try: x=x.tolist()
+    except:
+        try: x=list(x)
+        except: x=[x]
+    if len(x)>1:
+        return map(lambda x,tck=tck:spalde(x,tck),x)
+    d,ier=_fitpack._spalde(t,c,k,x[0])
+    if ier==0: return d
+    if ier==10:
+        raise TypeError,"Invalid input data. t(k)<=x<=t(n-k+1) must hold."
+    raise TypeError,"Unknown error"
+
+#def _curfit(x,y,w=None,xb=None,xe=None,k=3,task=0,s=None,t=None,
+#           full_output=0,nest=None,per=0,quiet=1):
+
+_surfit_cache = {'tx': array([],float),'ty': array([],float),
+                 'wrk': array([],float), 'iwrk':array([],int32)}
+def bisplrep(x,y,z,w=None,xb=None,xe=None,yb=None,ye=None,kx=3,ky=3,task=0,
+             s=None,eps=1e-16,tx=None,ty=None,full_output=0,
+             nxest=None,nyest=None,quiet=1):
+    """Find a bivariate B-spline representation of a surface.
+
+    Description:
+
+      Given a set of data points (x[i], y[i], z[i]) representing a surface
+      z=f(x,y), compute a B-spline representation of the surface.
+
+    Inputs:
+
+      x, y, z -- Rank-1 arrays of data points.
+      w -- Rank-1 array of weights. By default w=ones(len(x)).
+      xb, xe -- End points of approximation interval in x.
+      yb, ye -- End points of approximation interval in y.
+                By default xb, xe, yb, ye = x.min(), x.max(), y.min(), y.max()
+      kx, ky -- The degrees of the spline (1 <= kx, ky <= 5).  Third order
+                (kx=ky=3) is recommended.
+      task -- If task=0, find knots in x and y and coefficients for a given
+                smoothing factor, s.
+              If task=1, find knots and coefficients for another value of the
+                smoothing factor, s.  bisplrep must have been previously called
+                with task=0 or task=1.
+              If task=-1, find coefficients for a given set of knots tx, ty.
+      s -- A non-negative smoothing factor.  If weights correspond
+           to the inverse of the standard-deviation of the errors in z,
+           then a good s-value should be found in the range
+           (m-sqrt(2*m),m+sqrt(2*m)) where m=len(x)
+      eps -- A threshold for determining the effective rank of an
+             over-determined linear system of equations (0 < eps < 1)
+             --- not likely to need changing.
+      tx, ty -- Rank-1 arrays of the knots of the spline for task=-1
+      full_output -- Non-zero to return optional outputs.
+      nxest, nyest -- Over-estimates of the total number of knots.
+                      If None then nxest = max(kx+sqrt(m/2),2*kx+3),
+                                   nyest = max(ky+sqrt(m/2),2*ky+3)
+      quiet -- Non-zero to suppress printing of messages.
+
+    Outputs: (tck, {fp, ier, msg})
+
+      tck -- A list [tx, ty, c, kx, ky] containing the knots (tx, ty) and
+             coefficients (c) of the bivariate B-spline representation of the
+             surface along with the degree of the spline.
+
+      fp -- The weighted sum of squared residuals of the spline approximation.
+      ier -- An integer flag about splrep success.  Success is indicated if
+             ier<=0. If ier in [1,2,3] an error occurred but was not raised.
+             Otherwise an error is raised.
+      msg -- A message corresponding to the integer flag, ier.
+
+    Remarks:
+
+      SEE bisplev to evaluate the value of the B-spline given its tck
+      representation.
+
+    See also:
+      splprep, splrep, splint, sproot, splev - evaluation, roots, integral
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    x,y,z=map(myasarray,[x,y,z])
+    x,y,z=map(ravel,[x,y,z])  # ensure 1-d arrays.
+    m=len(x)
+    if not (m==len(y)==len(z)): raise TypeError, 'len(x)==len(y)==len(z) must hold.'
+    if w is None: w=ones(m,float)
+    else: w=myasarray(w)
+    if not len(w) == m: raise TypeError,' len(w)=%d is not equal to m=%d'%(len(w),m)
+    if xb is None: xb=x.min()
+    if xe is None: xe=x.max()
+    if yb is None: yb=y.min()
+    if ye is None: ye=y.max()
+    if not (-1<=task<=1): raise TypeError, 'task must be either -1,0, or 1'
+    if s is None: s=m-sqrt(2*m)
+    if tx is None and task==-1: raise TypeError, 'Knots_x must be given for task=-1'
+    if tx is not None: _surfit_cache['tx']=myasarray(tx)
+    nx=len(_surfit_cache['tx'])
+    if ty is None and task==-1: raise TypeError, 'Knots_y must be given for task=-1'
+    if ty is not None: _surfit_cache['ty']=myasarray(ty)
+    ny=len(_surfit_cache['ty'])
+    if task==-1 and nx<2*kx+2:
+        raise TypeError, 'There must be at least 2*kx+2 knots_x for task=-1'
+    if task==-1 and ny<2*ky+2:
+        raise TypeError, 'There must be at least 2*ky+2 knots_x for task=-1'
+    if not ((1<=kx<=5) and (1<=ky<=5)):
+        raise TypeError, 'Given degree of the spline (kx,ky=%d,%d) is not supported. (1<=k<=5)'%(kx,ky)
+    if m<(kx+1)*(ky+1): raise TypeError, 'm>=(kx+1)(ky+1) must hold'
+    if nxest is None: nxest=kx+sqrt(m/2)
+    if nyest is None: nyest=ky+sqrt(m/2)
+    nxest,nyest=max(nxest,2*kx+3),max(nyest,2*ky+3)
+    if task>=0 and s==0:
+        nxest=int(kx+sqrt(3*m))
+        nyest=int(ky+sqrt(3*m))
+    if task==-1:
+        _surfit_cache['tx']=myasarray(tx)
+        _surfit_cache['ty']=myasarray(ty)
+    tx,ty=_surfit_cache['tx'],_surfit_cache['ty']
+    wrk=_surfit_cache['wrk']
+    iwrk=_surfit_cache['iwrk']
+    u,v,km,ne=nxest-kx-1,nyest-ky-1,max(kx,ky)+1,max(nxest,nyest)
+    bx,by=kx*v+ky+1,ky*u+kx+1
+    b1,b2=bx,bx+v-ky
+    if bx>by: b1,b2=by,by+u-kx
+    lwrk1=u*v*(2+b1+b2)+2*(u+v+km*(m+ne)+ne-kx-ky)+b2+1
+    lwrk2=u*v*(b2+1)+b2
+    tx,ty,c,o = _fitpack._surfit(x,y,z,w,xb,xe,yb,ye,kx,ky,task,s,eps,
+                                   tx,ty,nxest,nyest,wrk,lwrk1,lwrk2)
+    _curfit_cache['tx']=tx
+    _curfit_cache['ty']=ty
+    _curfit_cache['wrk']=o['wrk']
+    ier,fp=o['ier'],o['fp']
+    tck=[tx,ty,c,kx,ky]
+    if ier<=0 and not quiet:
+        print _iermess2[ier][0]
+        print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
+                                                           len(ty),m,fp,s)
+    ierm=min(11,max(-3,ier))
+    if ierm>0 and not full_output:
+        if ier in [1,2,3,4,5]:
+            print "Warning: "+_iermess2[ierm][0]
+            print "\tkx,ky=%d,%d nx,ny=%d,%d m=%d fp=%f s=%f"%(kx,ky,len(tx),
+                                                           len(ty),m,fp,s)
+        else:
+            try:
+                raise _iermess2[ierm][1],_iermess2[ierm][0]
+            except KeyError:
+                raise _iermess2['unknown'][1],_iermess2['unknown'][0]
+    if full_output:
+        try:
+            return tck,fp,ier,_iermess2[ierm][0]
+        except KeyError:
+            return tck,fp,ier,_iermess2['unknown'][0]
+    else:
+        return tck
+
+def bisplev(x,y,tck,dx=0,dy=0):
+    """Evaluate a bivariate B-spline and its derivatives.
+
+    Description:
+
+      Return a rank-2 array of spline function values (or spline derivative
+      values) at points given by the cross-product of the rank-1 arrays x and y.
+      In special cases, return an array or just a float if either x or y or
+      both are floats.
+
+    Inputs:
+
+      x, y -- Rank-1 arrays specifying the domain over which to evaluate the
+              spline or its derivative.
+      tck -- A sequence of length 5 returned by bisplrep containing the knot
+             locations, the coefficients, and the degree of the spline:
+             [tx, ty, c, kx, ky].
+      dx, dy -- The orders of the partial derivatives in x and y respectively.
+
+    Outputs: (vals, )
+
+      vals -- The B-pline or its derivative evaluated over the set formed by
+              the cross-product of x and y.
+
+    Remarks:
+
+      SEE bisprep to generate the tck representation.
+
+    See also:
+      splprep, splrep, splint, sproot, splev - evaluation, roots, integral
+      UnivariateSpline, BivariateSpline - an alternative wrapping 
+              of the FITPACK functions
+    """
+    tx,ty,c,kx,ky=tck
+    if not (0<=dx<kx): raise ValueError,"0<=dx=%d<kx=%d must hold"%(dx,kx)
+    if not (0<=dy<ky): raise ValueError,"0<=dy=%d<ky=%d must hold"%(dy,ky)
+    x,y=map(myasarray,[x,y])
+    if (len(x.shape) != 1) or (len(y.shape) != 1):
+        raise ValueError, "First two entries should be rank-1 arrays."
+    z,ier=_fitpack._bispev(tx,ty,c,kx,ky,x,y,dx,dy)
+    if ier==10: raise ValueError,"Invalid input data"
+    if ier: raise TypeError,"An error occurred"
+    z.shape=len(x),len(y)
+    if len(z)>1: return z
+    if len(z[0])>1: return z[0]
+    return z[0][0]
+
+
+if __name__ == "__main__":
+    import sys,string
+    runtest=range(10)
+    if len(sys.argv[1:])>0:
+        runtest=map(string.atoi,sys.argv[1:])
+    put=sys.stdout.write
+    def norm2(x):
+        return dot(transpose(x),x)
+    def f1(x,d=0):
+        if d is None: return "sin"
+        if x is None: return "sin(x)"
+        if d%4 == 0: return sin(x)
+        if d%4 == 1: return cos(x)
+        if d%4 == 2: return -sin(x)
+        if d%4 == 3: return -cos(x)
+    def f2(x,y=0,dx=0,dy=0):
+        if x is None: return "sin(x+y)"
+        d=dx+dy
+        if d%4 == 0: return sin(x+y)
+        if d%4 == 1: return cos(x+y)
+        if d%4 == 2: return -sin(x+y)
+        if d%4 == 3: return -cos(x+y)
+    def test1(f=f1,per=0,s=0,a=0,b=2*pi,N=20,at=0,xb=None,xe=None):
+        if xb is None: xb=a
+        if xe is None: xe=b
+        x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
+        x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
+        v,v1=f(x),f(x1)
+        nk=[]
+        for k in range(1,6):
+            tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
+            if at:t=tck[0][k:-k]
+            else: t=x1
+            nd=[]
+            for d in range(k+1):
+                nd.append(norm2(f(t,d)-splev(t,tck,d)))
+            nk.append(nd)
+        print "\nf = %s  s=S_k(x;t,c)  x in [%s, %s] > [%s, %s]"%(f(None),
+                                                        `round(xb,3)`,`round(xe,3)`,
+                                                          `round(a,3)`,`round(b,3)`)
+        if at: str="at knots"
+        else: str="at the middle of nodes"
+        print " per=%d s=%s Evaluation %s"%(per,`s`,str)
+        print " k :  |f-s|^2  |f'-s'| |f''-.. |f'''-. |f''''- |f'''''"
+        k=1
+        for l in nk:
+            put(' %d : '%k)
+            for r in l:
+                put(' %.1e'%r)
+            put('\n')
+            k=k+1
+    def test2(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
+              ia=0,ib=2*pi,dx=0.2*pi):
+        if xb is None: xb=a
+        if xe is None: xe=b
+        x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
+        v=f(x)
+        nk=[]
+        for k in range(1,6):
+            tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
+            nk.append([splint(ia,ib,tck),spalde(dx,tck)])
+        print "\nf = %s  s=S_k(x;t,c)  x in [%s, %s] > [%s, %s]"%(f(None),
+                                                   `round(xb,3)`,`round(xe,3)`,
+                                                    `round(a,3)`,`round(b,3)`)
+        print " per=%d s=%s N=%d [a, b] = [%s, %s]  dx=%s"%(per,`s`,N,`round(ia,3)`,`round(ib,3)`,`round(dx,3)`)
+        print " k :  int(s,[a,b]) Int.Error   Rel. error of s^(d)(dx) d = 0, .., k"
+        k=1
+        for r in nk:
+            if r[0]<0: sr='-'
+            else: sr=' '
+            put(" %d   %s%.8f   %.1e "%(k,sr,abs(r[0]),
+                                         abs(r[0]-(f(ib,-1)-f(ia,-1)))))
+            d=0
+            for dr in r[1]:
+                put(" %.1e "%(abs(1-dr/f(dx,d))))
+                d=d+1
+            put("\n")
+            k=k+1
+    def test3(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
+              ia=0,ib=2*pi,dx=0.2*pi):
+        if xb is None: xb=a
+        if xe is None: xe=b
+        x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
+        v=f(x)
+        nk=[]
+        print "  k  :     Roots of s(x) approx %s  x in [%s,%s]:"%\
+              (f(None),`round(a,3)`,`round(b,3)`)
+        for k in range(1,6):
+            tck=splrep(x,v,s=s,per=per,k=k,nest=-1,xe=xe)
+            print '  %d  : %s'%(k,`sproot(tck).tolist()`)
+    def test4(f=f1,per=0,s=0,a=0,b=2*pi,N=20,xb=None,xe=None,
+              ia=0,ib=2*pi,dx=0.2*pi):
+        if xb is None: xb=a
+        if xe is None: xe=b
+        x=a+(b-a)*arange(N+1,dtype=float)/float(N)    # nodes
+        x1=a+(b-a)*arange(1,N,dtype=float)/float(N-1) # middle points of the nodes
+        v,v1=f(x),f(x1)
+        nk=[]
+        print " u = %s   N = %d"%(`round(dx,3)`,N)
+        print "  k  :  [x(u), %s(x(u))]  Error of splprep  Error of splrep "%(f(0,None))
+        for k in range(1,6):
+            tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
+            tck=splrep(x,v,s=s,per=per,k=k,nest=-1)
+            uv=splev(dx,tckp)
+            print "  %d  :  %s    %.1e           %.1e"%\
+                  (k,`map(lambda x:round(x,3),uv)`,
+                   abs(uv[1]-f(uv[0])),
+                   abs(splev(uv[0],tck)-f(uv[0])))
+        print "Derivatives of parametric cubic spline at u (first function):"
+        k=3
+        tckp,u=splprep([x,v],s=s,per=per,k=k,nest=-1)
+        for d in range(1,k+1):
+            uv=splev(dx,tckp,d)
+            put(" %s "%(`uv[0]`))
+        print
+    def makepairs(x,y):
+        x,y=map(myasarray,[x,y])
+        xy=array(map(lambda x,y:map(None,len(y)*[x],y),x,len(x)*[y]))
+        sh=xy.shape
+        xy.shape=sh[0]*sh[1],sh[2]
+        return transpose(xy)
+    def test5(f=f2,kx=3,ky=3,xb=0,xe=2*pi,yb=0,ye=2*pi,Nx=20,Ny=20,s=0):
+        x=xb+(xe-xb)*arange(Nx+1,dtype=float)/float(Nx)
+        y=yb+(ye-yb)*arange(Ny+1,dtype=float)/float(Ny)
+        xy=makepairs(x,y)
+        tck=bisplrep(xy[0],xy[1],f(xy[0],xy[1]),s=s,kx=kx,ky=ky)
+        tt=[tck[0][kx:-kx],tck[1][ky:-ky]]
+        t2=makepairs(tt[0],tt[1])
+        v1=bisplev(tt[0],tt[1],tck)
+        v2=f2(t2[0],t2[1])
+        v2.shape=len(tt[0]),len(tt[1])
+        print norm2(ravel(v1-v2))
+    if 1 in runtest:
+        print """\
+******************************************
+\tTests of splrep and splev
+******************************************"""
+        test1(s=1e-6)
+        test1()
+        test1(at=1)
+        test1(per=1)
+        test1(per=1,at=1)
+        test1(b=1.5*pi)
+        test1(b=1.5*pi,xe=2*pi,per=1,s=1e-1)
+    if 2 in runtest:
+        print """\
+******************************************
+\tTests of splint and spalde
+******************************************"""
+        test2()
+        test2(per=1)
+        test2(ia=0.2*pi,ib=pi)
+        test2(ia=0.2*pi,ib=pi,N=50)
+    if 3 in runtest:
+        print """\
+******************************************
+\tTests of sproot
+******************************************"""
+        test3(a=0,b=15)
+        print "Note that if k is not 3, some roots are missed or incorrect"
+    if 4 in runtest:
+        print """\
+******************************************
+\tTests of splprep, splrep, and splev
+******************************************"""
+        test4()
+        test4(N=50)
+    if 5 in runtest:
+        print """\
+******************************************
+\tTests of bisplrep, bisplev
+******************************************"""
+        test5()

Copied: trunk/Lib/sandbox/spline/fitpack.pyf (from rev 2289, trunk/Lib/interpolate/fitpack.pyf)
===================================================================
--- trunk/Lib/interpolate/fitpack.pyf	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/fitpack.pyf	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,463 @@
+!    -*- f90 -*-
+! Author: Pearu Peterson <pearu at cens.ioc.ee>
+! Modified by: John Travers <jtravs at gmail.com>
+!
+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
+
+  end interface
+end python module dfitpack
+

Copied: trunk/Lib/sandbox/spline/info.py (from rev 2289, trunk/Lib/interpolate/info.py)
===================================================================
--- trunk/Lib/interpolate/info.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/info.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,28 @@
+"""
+Spline Tools
+===================
+
+Wrappers around FITPACK functions:
+
+  splrep    -- find smoothing spline given (x,y) points on curve.
+  splprep   -- find smoothing spline given parametrically defined curve.
+  splev     -- evaluate the spline or its derivatives.
+  splint    -- compute definite integral of a spline.
+  sproot    -- find the roots of a cubic spline.
+  spalde    -- compute all derivatives of a spline at given points.
+  bisplrep   -- find bivariate smoothing spline representation.
+  bisplev   -- evaluate bivariate smoothing spline.
+
+  UnivariateSpline             -- A more recent, object-oriented wrapper;
+                                  finds a (possibly smoothed) interpolating
+				  spline.
+  InterpolatedUnivariateSpline
+  LSQUnivariateSpline
+  BivariateSpline              -- A more recent, object-oriented wrapper;
+                                  finds a interpolating spline for a 
+				  bivariate function.
+
+  SmoothBivariateSpline
+"""
+
+postpone_import = 1

Copied: trunk/Lib/sandbox/spline/setup.py (from rev 2289, trunk/Lib/interpolate/setup.py)
===================================================================
--- trunk/Lib/interpolate/setup.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/setup.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,25 @@
+#!/usr/bin/env python
+
+import os
+
+def configuration(parent_package='',top_path=None):
+    from numpy.distutils.misc_util import Configuration
+
+    config = Configuration('spline', parent_package, top_path)
+
+    config.add_library('fitpack',
+                       sources=[os.path.join('fitpack', '*.f')],
+                      )
+
+    config.add_extension('dfitpack',
+                         sources=['fitpack.pyf'],
+                         libraries=['fitpack'],
+                        )
+
+    config.add_data_dir('tests')
+
+    return config
+
+if __name__ == '__main__':
+    from numpy.distutils.core import setup
+    setup(**configuration(top_path='').todict())

Copied: trunk/Lib/sandbox/spline/spline.py (from rev 2289, trunk/Lib/interpolate/fitpack2.py)
===================================================================
--- trunk/Lib/interpolate/fitpack2.py	2006-10-16 17:38:02 UTC (rev 2289)
+++ trunk/Lib/sandbox/spline/spline.py	2006-10-18 19:27:50 UTC (rev 2290)
@@ -0,0 +1,521 @@
+"""
+spline --- curve and surface fitting with splines
+
+spline is based on a collection of Fortran routines DIERCKX
+by P. Dierckx (see http://www.netlib.org/dierckx/) transformed
+to double routines by Pearu Peterson.
+"""
+# Created by Pearu Peterson, June,August 2003
+# Modified by John Travers, October 2006
+
+__all__ = [
+    'UnivariateSpline',
+    'InterpolatedUnivariateSpline',
+    'LSQUnivariateSpline',
+
+    'LSQBivariateSpline',
+    'SmoothBivariateSpline',
+    'RectBivariateSpline']
+
+import warnings
+from numpy import zeros, concatenate, alltrue, ravel, all, diff
+
+import dfitpack
+
+################ Univariate spline ####################
+
+_curfit_messages = {1:"""
+The required storage space exceeds the available storage space, as
+specified by the parameter nest: nest too small. If nest is already
+large (say nest > m/2), it may also indicate that s is too small.
+The approximation returned is the weighted least-squares spline
+according to the knots t[0],t[1],...,t[n-1]. (n=nest) the parameter fp
+gives the corresponding weighted sum of squared residuals (fp>s).
+""",
+                    2:"""
+A theoretically impossible result was found during the iteration
+proces for finding a smoothing spline with fp = s: s too small.
+There is an approximation returned but the corresponding weighted sum
+of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
+                    3:"""
+The maximal number of iterations maxit (set to 20 by the program)
+allowed for finding a smoothing spline with fp=s has been reached: s
+too small.
+There is an approximation returned but the corresponding weighted sum
+of squared residuals does not satisfy the condition abs(fp-s)/s < tol.""",
+                    10:"""
+Error on entry, no approximation returned. The following conditions
+must hold:
+xb<=x[0]<x[1]<...<x[m-1]<=xe, w[i]>0, i=0..m-1
+if iopt=-1:
+  xb<t[k+1]<t[k+2]<...<t[n-k-2]<xe"""
+                    }
+
+class UnivariateSpline(object):
+    """ Univariate spline s(x) of degree k on the interval
+    [xb,xe] calculated from a given set of data points
+    (x,y).
+
+    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, y, w=None, bbox = [None]*2, k=3, s=None):
+        """
+        Input:
+          x,y   - 1-d sequences of data points (x must be
+                  in strictly ascending order)
+
+        Optional input:
+          w          - positive 1-d sequence of weights
+          bbox       - 2-sequence specifying the boundary of
+                       the approximation interval.
+                       By default, bbox=[x[0],x[-1]]
+          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].
+        """
+        #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
+        data = dfitpack.fpcurf0(x,y,k,w=w,
+                                xb=bbox[0],xe=bbox[1],s=s)
+        if data[-1]==1:
+            # nest too small, setting to maximum bound
+            data = self._reset_nest(data)
+        self._data = data
+        self._reset_class()
+
+    def _reset_class(self):
+        data = self._data
+        n,t,c,k,ier = data[7],data[8],data[9],data[5],data[-1]
+        self._eval_args = t[:n],c[:n],k
+        if ier==0:
+            # the spline returned has a residual sum of squares fp
+            # such that abs(fp-s)/s <= tol with tol a relative
+            # tolerance set to 0.001 by the program
+            pass
+        elif ier==-1:
+            # the spline returned is an interpolating spline
+            self.__class__ = InterpolatedUnivariateSpline
+        elif ier==-2:
+            # the spline returned is the weighted least-squares
+            # polynomial of degree k. In this extreme case fp gives
+            # the upper bound fp0 for the smoothing factor s.
+            self.__class__ = LSQUnivariateSpline
+        else:
+            # error
+            if ier==1:
+                self.__class__ = LSQUnivariateSpline
+            message = _curfit_messages.get(ier,'ier=%s' % (ier))
+            warnings.warn(message)
+
+    def _reset_nest(self, data, nest=None):
+        n = data[10]
+        if nest is None:
+            k,m = data[5],len(data[0])
+            nest = m+k+1 # this is the maximum bound for nest
+        else:
+            assert n<=nest,"nest can only be increased"
+        t,c,fpint,nrdata = data[8].copy(),data[9].copy(),\
+                           data[11].copy(),data[12].copy()
+        t.resize(nest)
+        c.resize(nest)
+        fpint.resize(nest)
+        nrdata.resize(nest)
+        args = data[:8] + (t,c,n,fpint,nrdata,data[13])
+        data = dfitpack.fpcurf1(*args)
+        return data
+
+    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:
+            warnings.warn('smoothing factor unchanged for'
+                          'LSQ spline with fixed knots')
+            return
+        args = data[:6] + (s,) + data[7:]
+        data = dfitpack.fpcurf1(*args)
+        if data[-1]==1:
+            # nest too small, setting to maximum bound
+            data = self._reset_nest(data)
+        self._data = data
+        self._reset_class()
+
+    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.
+        
+        """
+        if nu is None:
+            return dfitpack.splev(*(self._eval_args+(x,)))
+        return dfitpack.splder(nu=nu,*(self._eval_args+(x,)))
+
+    def get_knots(self):
+        """ Return the positions of (boundary and interior)
+        knots of the spline.
+        """
+        data = self._data
+        k,n = data[5],data[7]
+        return data[8][k:n-k]
+
+    def get_coeffs(self):
+        """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
+        approximation: sum ((w[i]*(y[i]-s(x[i])))**2,axis=0)
+        
+        """
+        return self._data[10]
+
+    def integral(self, a, b):
+        """ Return definite integral of the spline between two
+        given points.
+        """
+        return dfitpack.splint(*(self._eval_args+(a,b)))
+
+    def derivatives(self, x):
+        """ Return all derivatives of the spline at the point x."""
+        d,ier = dfitpack.spalde(*(self._eval_args+(x,)))
+        assert ier==0,`ier`
+        return d
+
+    def roots(self):
+        """ Return the zeros of the spline.
+
+        Restriction: only cubic splines are supported by fitpack.
+        """
+        k = self._data[5]
+        if k==3:
+            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'
+
+class InterpolatedUnivariateSpline(UnivariateSpline):
+    """ Interpolated univariate spline approximation. Identical to
+    UnivariateSpline with less error checking.
+
+    """
+
+    def __init__(self, x, y, w=None, bbox = [None]*2, k=3):
+        """
+        Input:
+          x,y   - 1-d sequences of data points (x must be
+                  in strictly ascending order)
+
+        Optional input:
+          w          - positive 1-d sequence of weights
+          bbox       - 2-sequence specifying the boundary of
+                       the approximation interval.
+                       By default, bbox=[x[0],x[-1]]
+          k=3        - degree of the univariate spline.
+        """
+        #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
+        self._data = dfitpack.fpcurf0(x,y,k,w=w,
+                                      xb=bbox[0],xe=bbox[1],s=0)
+        self._reset_class()
+
+class LSQUnivariateSpline(UnivariateSpline):
+    """ 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:
+          x,y   - 1-d sequences of data points (x must be
+                  in strictly ascending order)
+          t     - 1-d sequence of the positions of user-defined
+                  interior knots of the spline (t must be in strictly
+                  ascending order and bbox[0]<t[0]<...<t[-1]<bbox[-1])
+
+        Optional input:
+          w          - positive 1-d sequence of weights
+          bbox       - 2-sequence specifying the boundary of
+                       the approximation interval.
+                       By default, bbox=[x[0],x[-1]]
+          k=3        - degree of the univariate spline.
+        """
+        #_data == x,y,w,xb,xe,k,s,n,t,c,fp,fpint,nrdata,ier
+        xb=bbox[0]
+        xe=bbox[1]
+        if xb is None: xb = x[0]
+        if xe is None: xe = x[-1]
+        t = concatenate(([xb]*(k+1),t,[xe]*(k+1)))
+        n = len(t)
+        if not alltrue(t[k+1:n-k]-t[k:n-k-1] > 0,axis=0):
+            raise ValueError,\
+                  'Interior knots t must satisfy Schoenberg-Whitney conditions'
+        data = dfitpack.fpcurfm1(x,y,k,t,w=w,xb=xb,xe=xe)
+        self._data = data[:-3] + (None,None,data[-1])
+        self._reset_class()
+
+
+################ Bivariate spline ####################
+
+_surfit_messages = {1:"""
+The required storage space exceeds the available storage space: nxest
+or nyest too small, or s too small.
+The weighted least-squares spline corresponds to the current set of
+knots.""",
+                    2:"""
+A theoretically impossible result was found during the iteration
+process for finding a smoothing spline with fp = s: s too small or
+badly chosen eps.
+Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
+                    3:"""
+the maximal number of iterations maxit (set to 20 by the program)
+allowed for finding a smoothing spline with fp=s has been reached:
+s too small.
+Weighted sum of squared residuals does not satisfy abs(fp-s)/s < tol.""",
+                    4:"""
+No more knots can be added because the number of b-spline coefficients
+(nx-kx-1)*(ny-ky-1) already exceeds the number of data points m:
+either s or m too small.
+The weighted least-squares spline corresponds to the current set of
+knots.""",
+                    5:"""
+No more knots can be added because the additional knot would (quasi)
+coincide with an old one: s too small or too large a weight to an
+inaccurate data point.
+The weighted least-squares spline corresponds to the current set of
+knots.""",
+                    10:"""
+Error on entry, no approximation returned. The following conditions
+must hold:
+xb<=x[i]<=xe, yb<=y[i]<=ye, w[i]>0, i=0..m-1
+If iopt==-1, then
+  xb<tx[kx+1]<tx[kx+2]<...<tx[nx-kx-2]<xe
+  yb<ty[ky+1]<ty[ky+2]<...<ty[ny-ky-2]<ye""",
+                    -3:"""
+The coefficients of the spline returned have been computed as the
+minimal norm least-squares solution of a (numerically) rank deficient
+system (deficiency=%i). If deficiency is large, the results may be
+inaccurate. Deficiency may strongly depend on the value of eps."""
+                    }
+
+
+class BivariateSpline(object):
+    """ Bivariate spline s(x,y) of degrees kx and ky on the rectangle
+    [xb,xe] x [yb, ye] calculated from a given set of data points
+    (x,y,z).
+
+    See also:
+
+    bisplrep, bisplev - an older wrapping of FITPACK
+    UnivariateSpline - a similar class for univariate spline interpolation
+    SmoothUnivariateSpline - to create a BivariateSpline through the 
+                             given points
+    LSQUnivariateSpline - to create a BivariateSpline using weighted
+                          least-squares fitting
+    """
+
+    def get_residual(self):
+        """ Return weighted sum of squared residuals of the spline
+        approximation: sum ((w[i]*(z[i]-s(x[i],y[i])))**2,axis=0)
+        """
+        return self.fp
+    def get_knots(self):
+        """ Return a tuple (tx,ty) where tx,ty contain knots positions
+        of the spline with respect to x-, y-variable, respectively.
+        The position of interior and additional knots are given as
+          t[k+1:-k-1] and t[:k+1]=b, t[-k-1:]=e, respectively.
+        """
+        return self.tck[:2]
+    def get_coeffs(self):
+        """ Return spline coefficients."""
+        return self.tck[2]
+    def __call__(self,x,y,mth='array'):
+        """ Evaluate spline at positions x,y."""
+        if mth=='array':
+            tx,ty,c = self.tck[:3]
+            kx,ky = self.degrees
+            z,ier = dfitpack.bispev(tx,ty,c,kx,ky,x,y)
+            assert ier==0,'Invalid input: ier='+`ier`
+            return z
+        raise NotImplementedError
+
+class SmoothBivariateSpline(BivariateSpline):
+    """ Smooth bivariate spline approximation.
+
+    See also:
+
+    bisplrep, bisplev - an older wrapping of FITPACK
+    UnivariateSpline - a similar class for univariate spline interpolation
+    LSQUnivariateSpline - to create a BivariateSpline using weighted
+                          least-squares fitting
+    """
+
+    def __init__(self, x, y, z, w=None,
+                 bbox = [None]*4, kx=3, ky=3, s=None, eps=None):
+        """
+        Input:
+          x,y,z  - 1-d sequences of data points (order is not
+                   important)
+        Optional input:
+          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)]
+          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=len(w) which should be a good value
+                       if 1/w[i] is an estimate of the standard
+                       deviation of z[i].
+          eps        - a threshold for determining the effective rank
+                       of an over-determined linear system of
+                       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)
+        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
+
+class LSQBivariateSpline(BivariateSpline):
+    """ Weighted least-squares spline approximation.
+    See also:
+
+    bisplrep, bisplev - an older wrapping of FITPACK
+    UnivariateSpline - a similar class for univariate spline interpolation
+    SmoothUnivariateSpline - to create a BivariateSpline through the 
+                             given points
+    """
+
+    def __init__(self, x, y, z, tx, ty, w=None,
+                 bbox = [None]*4,
+                 kx=3, ky=3, eps=None):
+        """
+        Input:
+          x,y,z  - 1-d sequences of data points (order is not
+                   important)
+          tx,ty  - strictly ordered 1-d sequences of knots
+                   coordinates.
+        Optional input:
+          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)]
+          kx,ky=3,3  - degrees of the bivariate spline.
+          eps        - a threshold for determining the effective rank
+                       of an over-determined linear system of
+                       equations. 0 < eps < 1, default is 1e-16.
+        """
+        nx = 2*kx+2+len(tx)
+        ny = 2*ky+2+len(ty)
+        tx1 = zeros((nx,),float)
+        ty1 = zeros((ny,),float)
+        tx1[kx+1:nx-kx-1] = tx
+        ty1[ky+1:ny-ky-1] = ty
+
+        xb,xe,yb,ye = bbox
+        tx1,ty1,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx1,ty1,w,\
+                                               xb,xe,yb,ye,\
+                                               kx,ky,eps,lwrk2=1)
+        if ier>10:
+            tx1,ty1,c,fp,ier = dfitpack.surfit_lsq(x,y,z,tx1,ty1,w,\
+                                                   xb,xe,yb,ye,\
+                                                   kx,ky,eps,lwrk2=ier)
+        if ier in [0,-1,-2]: # normal return
+            pass
+        else:
+            if ier<-2:
+                deficiency = (nx-kx-1)*(ny-ky-1)+ier
+                message = _surfit_messages.get(-3) % (deficiency)
+            else:
+                message = _surfit_messages.get(ier,'ier=%s' % (ier))
+            warnings.warn(message)
+        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
+

Copied: trunk/Lib/sandbox/spline/tests (from rev 2289, trunk/Lib/interpolate/tests)



More information about the Scipy-svn mailing list