[Scipy-svn] r2902 - in trunk/Lib/sandbox: rbf rbf/tests spline spline/tests

scipy-svn@scip... scipy-svn@scip...
Tue Apr 10 17:06:18 CDT 2007


Author: jtravs
Date: 2007-04-10 17:04:33 -0500 (Tue, 10 Apr 2007)
New Revision: 2902

Added:
   trunk/Lib/sandbox/rbf/tests/example1.py
   trunk/Lib/sandbox/rbf/tests/example2.py
Removed:
   trunk/Lib/sandbox/rbf/tests/example.py
Modified:
   trunk/Lib/sandbox/rbf/rbf.py
   trunk/Lib/sandbox/rbf/tests/test_rbf.py
   trunk/Lib/sandbox/spline/fitpack.py
   trunk/Lib/sandbox/spline/fitpack.pyf
   trunk/Lib/sandbox/spline/tests/dierckx_test_data.py
   trunk/Lib/sandbox/spline/tests/test_fitpack.py
Log:
Changes to spline/rbf sandbox packages. New unit tests.


Modified: trunk/Lib/sandbox/rbf/rbf.py
===================================================================
--- trunk/Lib/sandbox/rbf/rbf.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/rbf/rbf.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -1,9 +1,9 @@
 #!/usr/bin/env python
-"""
-rbf - Radial basis functions for interpolation/smoothing scattered Nd data.
+"""rbf - Radial basis functions for interpolation/smoothing scattered Nd data.
 
 Written by John Travers <jtravs@gmail.com>, February 2007
 Based closely on Matlab code by Alex Chirokov
+Additional, large, improvements by Robert Hetland
 
 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
@@ -11,116 +11,141 @@
 
 NO WARRANTY IS EXPRESSED OR IMPLIED.  USE AT YOUR OWN RISK.
 
+Copyright (c) 2006-2007, Robert Hetland <hetland@tamu.edu>
+Copyright (c) 2007, John Travers <jtravs@gmail.com>
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+       notice, this list of conditions and the following disclaimer.
+
+    * Redistributions in binary form must reproduce the above
+       copyright notice, this list of conditions and the following
+       disclaimer in the documentation and/or other materials provided
+       with the distribution.
+
+    * Neither the name of Robert Hetland nor the names of any
+       contributors may be used to endorse or promote products derived
+       from this software without specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 """
 
-import scipy as s
+from numpy import sqrt, log, asarray, newaxis, all, dot, float64, eye
 import scipy.linalg
 
 class Rbf(object):
     """ A class for radial basis function approximation/interpolation of
         n-dimensional scattered data.
     """
-    def __init__(self,x,y, function='multiquadrics', constant=None, smooth=0):
+    
+    def _euclidean_norm(self, x1, x2):
+        return sqrt( ((x1 - x2)**2).sum(axis=0) )
+    
+    def _function(self, r):
+        if self.function.lower() == 'multiquadric':
+            return sqrt((self.epsilon*r)**2 + 1)
+        elif self.function.lower() == 'inverse multiquadric':
+            return 1.0/sqrt((self.epsilon*r)**2 + 1)
+        elif self.function.lower() == 'gausian':
+            return exp(-(self.epsilon*r)**2)
+        elif self.function.lower() == 'cubic':
+            return r**3
+        elif self.function.lower() == 'quintic':
+            return r**5
+        elif self.function.lower() == 'thin-plate':
+            return r**2 * log(r)
+        else:
+            raise ValueError, 'Invalid basis function name'
+    
+    def __init__(self, *args, **kwargs):
         """ Constructor for Rbf class.
-
+            
             Inputs:
-                x   (dim, n) array of coordinates for the nodes
-                y   (n,) array of values at the nodes
-                function    the radial basis function
-                            'linear', 'cubic' 'thinplate', 'multiquadrics'
-                            or 'gaussian', default is 'multiquadrics'
-                constant    adjustable constant for gaussian or multiquadrics
+                x, y, z, ..., d    
+                            Where x, y, z, ... are the coordinates of the nodes
+                            and d is the array of values at the nodes
+                
+                function    the radial basis function, based on the radius, r, given 
+                            by the norm (defult is Euclidean distance); the default
+                            is 'multiquadratic'.
+                            
+                            'multiquadric': sqrt((self.epsilon*r)**2 + 1)
+                            'inverse multiquadric': 1.0/sqrt((self.epsilon*r)**2 + 1)
+                            'gausian': exp(-(self.epsilon*r)**2)
+                            'cubic': r**3
+                            'quintic': r**5
+                            'thin-plate': r**2 * log(r)
+                            
+                epsilon     adjustable constant for gaussian or multiquadrics
                             functions - defaults to approximate average distance
                             between nodes (which is a good start)
+                            
                 smooth      values greater than zero increase the smoothness
                             of the approximation. 
                             0 is for interpolation (default), the function will
                             always go through the nodal points in this case.
+                            
+                norm        A function that returns the 'distance' between two points,
+                            with inputs as arrays of positions (x, y, z, ...), and an
+                            output as an array of distance.  E.g, the default is
+                            
+                                def euclidean_norm(self, x1, x2):
+                                    return sqrt( ((x1 - x2)**2).sum(axis=0) )
+                            
+                            which is called with x1 = x1[ndims, newaxis, :] and 
+                            x2 = x2[ndims, :, newaxis] such that the result is a 
+                            symetric, square matrix of the distances between each point
+                            to each other point.
+            
+            Outputs: 
+                Interpolator object rbfi that returns interpolated values at new positions:
+                >>> rbfi = Rbf(x, y, z, d)      # radial basis function interpolator instance
+                >>> di = rbfi(xi, yi, zi)       # interpolated values
+        """        
+        self.xi = asarray([asarray(a, dtype=float64).flatten() for a in args[:-1]])
+        self.N = self.xi.shape[-1]
+        self.di = asarray(args[-1], dtype=float64).flatten()
+        
+        assert [x.size==self.di.size for x in self.xi], \
+               'All arrays must be equal length'
+        
+        self.norm = kwargs.pop('norm', self._euclidean_norm)
+        r = self._call_norm(self.xi, self.xi)
+        self.epsilon = kwargs.pop('epsilon', r.mean())
+        self.function = kwargs.pop('function', 'multiquadric')
+        self.smooth = kwargs.pop('smooth', 0.0)
+        
+        self.A = self._function(r) - eye(self.N)*self.smooth
+        self.nodes = scipy.linalg.solve(self.A, self.di)
+    
+    def _call_norm(self, x1, x2):
+        if len(x1.shape) == 1:
+            x1 = x1[newaxis, :]
+        if len(x2.shape) == 1:
+            x2 = x2[newaxis, :]
+        x1 = x1[..., :, newaxis]
+        x2 = x2[..., newaxis, :]
+        return self.norm(x1, x2)
+    
+    def __call__(self, *args):
+        assert all([x.shape == y.shape \
+                    for x in args \
+                    for y in args]), 'Array lengths must be equal'
+        shp = args[0].shape
+        self.xa = asarray([a.flatten() for a in args], dtype=float64)
+        r = self._call_norm(self.xa, self.xi)
+        return dot(self._function(r), self.nodes).reshape(shp)
 
-            Outputs: None
-        """
-        if len(x.shape) == 1:
-            nxdim = 1
-            nx = x.shape[0]
-        else:
-            (nxdim, nx)=x.shape
-        if len(y.shape) == 1:
-            nydim = 1
-            ny = y.shape[0]
-        else:
-            (nydim, ny)=y.shape
-        x.shape = (nxdim, nx)
-        y.shape = (nydim, ny)
-        if nx != ny:
-            raise ValueError, 'x and y should have the same number of points'
-        if nydim != 1:
-            raise ValueError, 'y should be a length n vector'
-        self.x = x
-        self.y = y
-        self.function = function
-        if (constant==None 
-            and ((function == 'multiquadrics') or (function == 'gaussian'))):
-            # approx. average distance between the nodes
-            constant = (s.product(x.T.max(0)-x.T.min(0),axis=0)/nx)**(1/nxdim)
-        self.constant = constant
-        self.smooth = smooth
-        if self.function == 'linear':
-            self.phi = lambda r: r
-        elif self.function == 'cubic':
-            self.phi = lambda r: r*r*r
-        elif self.function == 'multiquadrics':
-            self.phi = lambda r: s.sqrt(1.0+r*r/(self.constant*self.constant))
-        elif self.function == 'thinplate':
-            self.phi = lambda r: r*r*s.log(r+1)
-        elif self.function == 'gaussian':
-            self.phi = lambda r: s.exp(-0.5*r*r/(self.rbfconst*self.constant))
-        else:
-            raise ValueError, 'unkown function'
-        A = self._rbf_assemble()
-        b=s.r_[y.T, s.zeros((nxdim+1, 1), float)]
-        self.coeff = s.linalg.solve(A,b)
-
-    def __call__(self, xi):
-        """ Evaluate the radial basis function approximation at points xi.
-
-            Inputs:
-                xi  (dim, n) array of coordinates for the points to evaluate at
-
-            Outputs:
-                y   (n,) array of values at the points xi
-        """
-        if len(xi.shape) == 1:
-            nxidim = 1
-            nxi = xi.shape[0]
-        else:
-            (nxidim, nxi)=xi.shape
-        xi.shape = (nxidim, nxi)
-        (nxdim, nx) = self.x.shape
-        if nxdim != nxidim:
-            raise ValueError, 'xi should have the same number of rows as an' \
-                              ' array used to create RBF interpolation'
-        f = s.zeros(nxi, float)
-        r = s.zeros(nx, float)
-        for i in range(nxi):
-            st=0.0
-            r = s.dot(xi[:,i,s.newaxis],s.ones((1,nx))) - self.x
-            r = s.sqrt(sum(r*r))
-            st = self.coeff[nx,:] + s.sum(self.coeff[0:nx,:].flatten()*self.phi(r))
-            for k in range(nxdim):
-                st=st+self.coeff[k+nx+1,:]*xi[k,i]
-            f[i] = st
-        return f
-
-    def _rbf_assemble(self):
-        (nxdim, nx)=self.x.shape
-        A=s.zeros((nx,nx), float)
-        for i in range(nx):
-            for j in range(i+1):
-                r=s.linalg.norm(self.x[:,i]-self.x[:,j])
-                temp=self.phi(r)
-                A[i,j]=temp
-                A[j,i]=temp
-            A[i,i] = A[i,i] - self.smooth
-        P = s.c_[s.ones((nx,1), float), self.x.T]
-        A = s.r_[s.c_[A, P], s.c_[P.T, s.zeros((nxdim+1,nxdim+1), float)]]
-        return A

Deleted: trunk/Lib/sandbox/rbf/tests/example.py
===================================================================
--- trunk/Lib/sandbox/rbf/tests/example.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/rbf/tests/example.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -1,52 +0,0 @@
-import scipy as s
-import scipy.interpolate
-
-from scipy.sandbox.rbf import Rbf
-
-import matplotlib
-matplotlib.use('Agg')
-import pylab as p
-
-# 1d tests - setup data
-x = s.linspace(0,10,9)
-y = s.sin(x) 
-xi = s.linspace(0,10,101)
-
-# use interpolate methods
-ius = s.interpolate.InterpolatedUnivariateSpline(x,y)
-yi = ius(xi)
-p.subplot(2,1,1)
-p.plot(x,y,'o',xi,yi, xi, s.sin(xi),'r')
-p.title('Interpolation using current scipy fitpack2')
-
-# use RBF method
-rbf = Rbf(x, y)
-fi = rbf(xi)
-p.subplot(2,1,2)
-p.plot(x,y,'bo',xi.flatten(),fi.flatten(),'g',xi.flatten(),
-                                                    s.sin(xi.flatten()),'r')
-p.title('RBF interpolation - multiquadrics')
-p.savefig('rbf1dtest.png')
-p.close()
-
-# 2-d tests - setup scattered data
-x = s.rand(50,1)*4-2
-y = s.rand(50,1)*4-2
-z = x*s.exp(-x**2-y**2)
-ti = s.linspace(-2.0,2.0,81)
-(XI,YI) = s.meshgrid(ti,ti)
-
-# use RBF
-rbf = Rbf(s.c_[x.flatten(),y.flatten()].T,z.T,constant=2)
-ZI = rbf(s.c_[XI.flatten(), YI.flatten()].T)
-ZI.shape = XI.shape
-
-# plot the result
-from enthought.tvtk.tools import mlab
-f=mlab.figure(browser=False)
-su=mlab.Surf(XI,YI,ZI,ZI,scalar_visibility=True)
-f.add(su)
-su.lut_type='blue-red'
-f.objects[0].axis.z_label='value'
-pp = mlab.Spheres(s.c_[x.flatten(), y.flatten(), z.flatten()],radius=0.03)
-f.add(pp)
\ No newline at end of file

Copied: trunk/Lib/sandbox/rbf/tests/example1.py (from rev 2897, trunk/Lib/sandbox/rbf/tests/example.py)
===================================================================
--- trunk/Lib/sandbox/rbf/tests/example.py	2007-04-04 17:27:42 UTC (rev 2897)
+++ trunk/Lib/sandbox/rbf/tests/example1.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -0,0 +1,50 @@
+import scipy as s
+import scipy.interpolate
+
+from scipy.sandbox.rbf import Rbf
+
+import matplotlib
+matplotlib.use('Agg')
+import pylab as p
+
+# 1d tests - setup data
+x = s.linspace(0,10,9)
+y = s.sin(x) 
+xi = s.linspace(0,10,101)
+
+# use interpolate methods
+ius = s.interpolate.InterpolatedUnivariateSpline(x,y)
+yi = ius(xi)
+p.subplot(2,1,1)
+p.plot(x,y,'o',xi,yi, xi, s.sin(xi),'r')
+p.title('Interpolation using current scipy fitpack2')
+
+# use RBF method
+rbf = Rbf(x, y)
+fi = rbf(xi)
+p.subplot(2,1,2)
+p.plot(x,y,'bo',xi,fi,'g',xi, s.sin(xi),'r')
+p.title('RBF interpolation - multiquadrics')
+p.show()
+
+# 2-d tests - setup scattered data
+x = s.rand(50,1)*4-2
+y = s.rand(50,1)*4-2
+z = x*s.exp(-x**2-y**2)
+ti = s.linspace(-2.0,2.0,81)
+(XI,YI) = s.meshgrid(ti,ti)
+
+# use RBF
+rbf = Rbf(x.flatten(),y.flatten(),z.flatten(),eps=2)
+ZI = rbf(XI.flatten(), YI.flatten())
+ZI.shape = XI.shape
+
+# plot the result
+from enthought.tvtk.tools import mlab
+f=mlab.figure(browser=False)
+su=mlab.Surf(XI,YI,ZI,ZI,scalar_visibility=True)
+f.add(su)
+su.lut_type='blue-red'
+f.objects[0].axis.z_label='value'
+pp = mlab.Spheres(s.c_[x.flatten(), y.flatten(), z.flatten()],radius=0.03)
+f.add(pp)

Added: trunk/Lib/sandbox/rbf/tests/example2.py
===================================================================
--- trunk/Lib/sandbox/rbf/tests/example2.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/rbf/tests/example2.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -0,0 +1,47 @@
+from numpy import sin, asarray, exp, random, mgrid, pi, cos, sqrt, ones
+from scipy.sandbox.rbf import Rbf
+import pylab as pl
+
+def truth_2d(x,y,w=2*pi):
+    "moguls"
+    return sin(w*x)*cos(w*y)
+
+def truth_nd(*args):
+    "a gausian sphere"
+    x = asarray(list(args), 'float64')
+    return exp( -sqrt((x**2).sum(axis=0)) )
+
+# 2D example
+N = 300
+xi = random.rand(N)
+yi = random.rand(N)
+di = truth_2d(xi, yi)
+xa, ya = mgrid[0:1:50j, 0:1:50j]
+s = Rbf(xi, yi, di)
+da = s(xa, ya)
+pl.figure()
+n = pl.normalize(-1., 1.)
+pl.pcolor(xa, ya, da, norm=n, cmap=pl.cm.jet)
+pl.scatter(xi, yi, 100, di, norm=n, cmap=pl.cm.jet)
+pl.axis([0., 1., 0., 1.])
+pl.colorbar()
+pl.draw()
+# 3d example
+N = 300
+xi = 2.*random.randn(N)
+yi = 2.*random.randn(N)
+zi = 2.*random.randn(N)
+di = truth_nd(xi, yi, zi)
+zas = [-0.25, 0.0, 0.25, 0.75]
+xa, ya = mgrid[-1:1:50j, -1:1:50j]
+s = Rbf(xi, yi, zi, di)
+fig = pl.figure(figsize=(12, 3))
+for idx, za in enumerate(zas):
+    da = s(xa, ya, za*ones(xa.shape, 'f'))
+    ax = fig.add_subplot(1,4,idx+1)
+    ax.pcolor(xa, ya, da, norm=pl.normalize(0, 1), \
+	  shading='flat', cmap=pl.cm.jet)
+    ax.set_aspect('equal')
+
+pl.show()
+

Modified: trunk/Lib/sandbox/rbf/tests/test_rbf.py
===================================================================
--- trunk/Lib/sandbox/rbf/tests/test_rbf.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/rbf/tests/test_rbf.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -1,9 +1,9 @@
 #!/usr/bin/env python
-# Created by John Travers, February 2007
+# Created by John Travers, Robert Hetland, 2007
 """ Test functions for rbf module """
 
 from numpy.testing import *
-import numpy as n
+from numpy import linspace, sin, random, exp 
 
 set_package_path()
 from rbf.rbf import Rbf
@@ -11,21 +11,32 @@
 
 class test_Rbf1D(NumpyTestCase):
     def check_multiquadrics(self):
-        x = n.linspace(0,10,9)
-        y = n.sin(x) 
+        x = linspace(0,10,9)
+        y = sin(x) 
         rbf = Rbf(x, y)
         yi = rbf(x)
-        assert_array_almost_equal(y.flatten(), yi)
+        assert_array_almost_equal(y, yi)
 
 class test_Rbf2D(NumpyTestCase):
     def check_multiquadrics(self):
-        x = n.random.rand(50,1)*4-2
-        y = n.random.rand(50,1)*4-2
-        z = x*n.exp(-x**2-y**2)
-        rbf = Rbf(n.c_[x.flatten(),y.flatten()].T,z.T,constant=2)
-        zi = rbf(n.c_[x.flatten(), y.flatten()].T)
+        x = random.rand(50,1)*4-2
+        y = random.rand(50,1)*4-2
+        z = x*exp(-x**2-y**2)
+        rbf = Rbf(x, y, z ,epsilon=2)
+        zi = rbf(x, y)
         zi.shape = x.shape
         assert_array_almost_equal(z, zi)
 
+class test_Rbf3D(NumpyTestCase):
+    def check_multiquadrics(self):
+        x = random.rand(50,1)*4-2
+        y = random.rand(50,1)*4-2
+        z = random.rand(50,1)*4-2
+        d = x*exp(-x**2-y**2)
+        rbf = Rbf(x, y, z, d ,epsilon=2)
+        di = rbf(x, y, z)
+        di.shape = x.shape
+        assert_array_almost_equal(di, d)
+    
 if __name__ == "__main__":
-    NumpyTest().run()
\ No newline at end of file
+    NumpyTest().run()

Modified: trunk/Lib/sandbox/spline/fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/spline/fitpack.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -25,8 +25,8 @@
     $ 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
+    For univariate splines: cocosp, concon, fourco
+    For bivariate splines: profil, parsur, surev
 """
 
 __all__ = ['splrep', 'splprep', 'splev', 'splint', 'sproot', 'spalde',
@@ -219,27 +219,36 @@
     if task==1:
         try: 
             u=_parcur_cache['u']
-            ub=_parcur_cache['ub']
-            ue=_parcur_cache['ue']
             t=_parcur_cache['t']
             wrk=_parcur_cache['wrk']
             iwrk=_parcur_cache['iwrk']
             n=_parcur_cache['n']
+            if not per:
+                ub=_parcur_cache['ub']
+                ue=_parcur_cache['ue']
         except KeyError:
             raise ValueError, 'task=1 can only be called after task=0'
-        u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth1(ipar,idim,u,x,w,
-                                             ub,ue,nest,n,t,wrk,iwrk,k=k,s=s)
+        if per:
+            u,n,t,c,fp,wrk,iwrk,ier=dfitpack.clocur_smth1(ipar,idim,u,x,
+                                             w,n,t,wrk,iwrk,k=k,s=s)
+        else:
+            u,ub,ue,n,t,c,fp,wrk,iwrk,ier=dfitpack.parcur_smth1(ipar,idim,u,x,
+                                             w,ub,ue,n,t,wrk,iwrk,k=k,s=s)
     if task==-1:
-        u,ub,ue,n,t,c,fp,ier=dfitpack.parcur_lsq(ipar,idim,u,x,w,ub,ue,
-                                                                nest,n,t,k=k)
+        if per:
+            u,n,t,c,fp,ier=dfitpack.clocur_lsq(ipar,idim,u,x,w,n,t,k=k)
+        else:
+            u,ub,ue,n,t,c,fp,ier=dfitpack.parcur_lsq(ipar,idim,u,x,w,ub,ue,
+                                                                n,t,k=k)
     if task>=0:
         _parcur_cache['n']=n
         _parcur_cache['u']=u
-        _parcur_cache['ub']=ub
-        _parcur_cache['ue']=ue
         _parcur_cache['t']=t
         _parcur_cache['wrk']=wrk
         _parcur_cache['iwrk']=iwrk
+        if not per:
+            _parcur_cache['ub']=ub
+            _parcur_cache['ue']=ue
     c = c[:n*idim]
     c.shape=idim,n
     c = c[:,:n-k-1]
@@ -379,7 +388,7 @@
             n,t,c,fp,wrk,iwrk,ier = dfitpack.percur_smth1(x,y,w,n,t,wrk,iwrk,
                                                                     k=k,s=s)
         elif task==-1:
-            n,t,c,fp,ier = dfitpack.percur_lsq(x,y,w,n,t,k=k)
+            n,t,c,fp,ier = dfitpack.percur_lsq(x,y,w,t,k=k)
         if task>=0:
             _percur_cache['t']=t
             _percur_cache['wrk']=wrk

Modified: trunk/Lib/sandbox/spline/fitpack.pyf
===================================================================
--- trunk/Lib/sandbox/spline/fitpack.pyf	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/spline/fitpack.pyf	2007-04-10 22:04:33 UTC (rev 2902)
@@ -147,7 +147,7 @@
        integer optional,check(1<=k && k <=5),intent(in) :: k=3
        real*8 intent(hide),check(s>=0.0) :: s = 0.0
        integer intent(hide),depend(m,k) :: nest=m+2*k
-       integer intent(in,out) :: n
+       integer intent(out),depend(t) :: n=len(t)
        real*8 dimension(n),intent(in,out) :: t
        real*8 dimension(n),intent(out) :: c
        real*8 intent(out) :: fp
@@ -167,10 +167,10 @@
        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(m,k) :: nest=m+2*k
-       integer intent(in,out),depend(nest) :: n=nest
+       integer intent(hide),depend(t) :: nest=len(t)
+       integer intent(in,out) :: n
        real*8 dimension(nest),intent(in,out) :: t
-       real*8 dimension(n),intent(out) :: c
+       real*8 dimension(nest),intent(out) :: c
        real*8 intent(out) :: fp
        real*8 dimension(lwrk),intent(in,out) :: wrk
        integer intent(hide),depend(m,k,nest) :: lwrk=m*(k+1)+nest*(8+5*k)
@@ -189,9 +189,9 @@
        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(m,k) :: nest=m+2*k
-       integer intent(out),depend(nest) :: n=nest
+       integer intent(out) :: n
        real*8 dimension(nest),intent(out) :: t
-       real*8 dimension(n),intent(out) :: c
+       real*8 dimension(nest),intent(out) :: c
        real*8 intent(out) :: fp
        real*8 dimension(lwrk),intent(out) :: wrk
        integer intent(hide),depend(m,k,nest) :: lwrk=m*(k+1)+nest*(8+5*k)
@@ -215,15 +215,15 @@
        real*8 intent(in,out) :: ue
        integer check(1<=k && k<=5) :: k=3.0
        real*8 intent(hide),check(s>=0.0) :: s = 0.0
-       integer intent(in) :: nest
+       integer intent(hide) :: nest=n
        integer intent(in,out) :: n
        real*8 dimension(nest), intent(in,out) :: t
        integer intent(hide), depend(nest,idim) :: nc=idim*nest
        real*8 dimension(nc), intent(out) :: c
        real*8 intent(out) :: fp
-       real*8 dimension(lwrk), intent(cache) :: wrk
+       real*8 dimension(lwrk), intent(hide,cache) :: wrk
        integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
-       integer dimension(nest), intent(cache) :: iwrk
+       integer dimension(nest), intent(hide,cache) :: iwrk
        integer intent(out) :: ier
      end subroutine parcur_lsq
 
@@ -256,6 +256,59 @@
        integer intent(out) :: ier
      end subroutine parcur_smth1
 
+     subroutine parcur_smth0(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,&
+                                                        c,fp,wrk,lwrk,iwrk,ier)
+       fortranname parcur
+       integer intent(hide) :: iopt = 0
+       integer check(ipar == 1 || ipar == 0) :: ipar
+       integer check(idim > 0 && idim < 11) :: idim
+       integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+       real*8 dimension(m), intent(in,out) :: u
+       integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+       real*8 dimension(mx) :: x
+       real*8 dimension(m) :: w
+       real*8 intent(in,out) :: ub
+       real*8 intent(in,out) :: ue
+       integer check(1<=k && k<=5) :: k=3.0
+       real*8 check(s>=0.0) :: s = 0.0
+       integer intent(in) :: nest
+       integer intent(out) :: n
+       real*8 dimension(nest), intent(out) :: t
+       integer intent(hide), depend(nest,idim) :: nc=idim*nest
+       real*8 dimension(nc), intent(out) :: c
+       real*8 intent(out) :: fp
+       real*8 dimension(lwrk), intent(out) :: wrk
+       integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
+       integer dimension(nest), intent(out) :: iwrk
+       integer intent(out) :: ier
+     end subroutine parcur_smth0
+
+     subroutine clocur_lsq(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,&
+                                                             wrk,lwrk,iwrk,ier)
+       !u,n,t,c,fp,wrk,iwrk,ier=clocur_lsq(ipar,idim,u,x,w,nest,n,t,[k,s])
+       fortranname clocur
+       integer intent(hide) :: iopt = -1
+       integer check(ipar == 1 || ipar == 0) :: ipar
+       integer check(idim > 0 && idim < 11) :: idim
+       integer intent(hide),depend(u,k),check(m>k) :: m=len(u)
+       real*8 dimension(m), intent(in,out) :: u
+       integer intent(hide),depend(x,idim,m),check(mx>=idim*m) :: mx=len(x)
+       real*8 dimension(mx) :: x
+       real*8 dimension(m) :: w
+       integer check(1<=k && k<=5) :: k=3.0
+       real*8 intent(hide) :: s = 0.0
+       integer intent(hide) :: nest=n
+       integer intent(in,out) :: n
+       real*8 dimension(nest), intent(in,out) :: t
+       integer intent(hide), depend(nest,idim) :: nc=idim*nest
+       real*8 dimension(nc), intent(out) :: c
+       real*8 intent(out) :: fp
+       real*8 dimension(lwrk), intent(cache,hide) :: wrk
+       integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(7+idim+5*k)
+       integer dimension(nest), intent(cache,hide) :: iwrk
+       integer intent(out) :: ier
+     end subroutine clocur_lsq
+
      subroutine clocur_smth0(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,&
                                                              wrk,lwrk,iwrk,ier)
        !u,n,t,c,fp,wrk,iwrk,ier=clocur_smth0(ipar,idim,u,x,w,nest,[k,s])
@@ -282,10 +335,12 @@
        integer intent(out) :: ier
      end subroutine clocur_smth0
 
-     subroutine parcur_smth0(iopt,ipar,idim,m,u,mx,x,w,ub,ue,k,s,nest,n,t,nc,&
-                                                        c,fp,wrk,lwrk,iwrk,ier)
-       fortranname parcur
-       integer intent(hide) :: iopt = 0
+     subroutine clocur_smth1(iopt,ipar,idim,m,u,mx,x,w,k,s,nest,n,t,nc,c,fp,&
+                                                             wrk,lwrk,iwrk,ier)
+       !u,n,t,c,fp,wrk,iwrk,ier=clocur_smth1(ipar,idim,u,x,w,nest,
+       !                                                    n,t,wrk,iwrk,[k,s])
+       fortranname clocur
+       integer intent(hide) :: iopt = 1
        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)
@@ -293,21 +348,19 @@
        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 intent(in,out) :: ub
-       real*8 intent(in,out) :: ue
        integer check(1<=k && k<=5) :: k=3.0
        real*8 check(s>=0.0) :: s = 0.0
        integer intent(in) :: nest
-       integer intent(out) :: n
-       real*8 dimension(nest), intent(out) :: t
+       integer intent(in,out) :: n
+       real*8 dimension(nest), intent(in,out) :: t
        integer intent(hide), depend(nest,idim) :: nc=idim*nest
        real*8 dimension(nc), intent(out) :: c
        real*8 intent(out) :: fp
-       real*8 dimension(lwrk), intent(out) :: wrk
-       integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(6+idim+3*k)
-       integer dimension(nest), intent(out) :: iwrk
+       real*8 dimension(lwrk), intent(in,out) :: wrk
+       integer intent(hide),depend(m,k,nest,idim) :: lwrk=m*(k+1)+nest*(7+idim+5*k)
+       integer dimension(nest), intent(in,out) :: iwrk
        integer intent(out) :: ier
-     end subroutine parcur_smth0
+     end subroutine clocur_smth1
 
      subroutine insert(iopt,t,n,c,k,x,tt,nn,cc,nest,ier)
        ! tt, nn, cc, ier = insert(per, t, c, k, x, nest)

Modified: trunk/Lib/sandbox/spline/tests/dierckx_test_data.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/dierckx_test_data.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/spline/tests/dierckx_test_data.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -92,6 +92,460 @@
                3.0,  3.0,  2.4,  4.2,  3.5])]
 }
 
+percur_test = {
+'x' : array([0.0,3.922,7.843,11.765,15.686,
+     19.608,23.509,27.451,31.373,35.294,39.216,43.137,47.059,50.980,
+     54.902,58.824,62.745,66.667,70.588,74.510,78.431,82.353,86.275,
+     90.196,94.118,98.039, 100.0]),
+'y' : array([10.099,14.835,21.453,25.022,22.427,
+     22.315,22.070,19.673,16.754,13.983,11.973,12.286,16.129,21.560,
+     28.041,39.205,59.489,72.559,75.960,79.137,75.925,68.809,55.758,
+     39.915,22.006,12.076, 10.099]),
+'k' : [3, 3, 3, 3, 3, 3, 3, 5, 5, 5, 5, 5, 5], 
+'iopt' : [0, 1, 1, 1, 0, 0, -1, 0, 1, 1, 1, 0, 0, -1], 
+'s' : [65000.0, 500.0, 5.0, 20.0, 20.0, 0.0, -1.0, 65000.0, 500.0, 5.0, 
+                                                    20.0, 20.0, 0.0, -1],
+'res' : [0.140221E+05,0.500033E+03,0.499936E+01,0.200010E+02,0.199981E+02,
+         0.000000E+00,0.638916E+02,0.140221E+05,0.499538E+03,0.500125E+01,
+         0.200032E+02, 0.199998E+02, 0.000000E+00, 0.627153E+02],  
+'err' : [-2, 0, 0, 0, 0, -1, 0, -2, 0, 0, 0, 0, -1, 0],
+'knots' : [[-300.000, -200.000, -100.000,    0.000,  100.000,  200.000,
+             300.000,  400.000], 
+         [-49.020,  -33.333,  -21.569,    0.000,   15.686,   27.451,
+           50.980,   66.667, 78.431,  100.000,  115.686,  127.451,  150.980], 
+      [-21.569,   -9.804,   -1.961,    0.000,    7.843,   15.686,   27.451,
+      39.216, 50.980,   58.824,   62.745,   66.667,   70.588,   74.510,
+      78.431,   90.196, 98.039,  100.000,  107.843,  115.686,  127.451], 
+      [-21.569,   -9.804,   -1.961,    0.000,    7.843,   15.686,   27.451,
+      39.216, 50.980,   58.824,   62.745,   66.667,   70.588,   74.510,
+      78.431,   90.196,  98.039,  100.000,  107.843,  115.686,  127.451], 
+      [-33.333,  -21.569,   -9.804,    0.000,    7.843,   15.686,   27.451,
+      39.216, 50.980,   58.824,   62.745,   66.667,   78.431,   90.196,
+      100.000,  107.843, 115.686,  127.451], 
+      [-9.804,   -5.882,   -1.961,    0.000,    3.922,    7.843,   11.765,
+      15.686, 19.608,   23.509,   27.451,   31.373,   35.294,   39.216,
+      43.137,   47.059, 50.980,   54.902,   58.824,   62.745,   66.667, 
+      70.588,   74.510,   78.431, 82.353,   86.275,   90.196,   94.118, 
+      98.039,  100.000,  103.922,  107.843,  111.765], 
+      [-30.000,  -20.000,  -10.000,    0.000,   10.000,   20.000,   30.000,
+      40.000, 50.000,   60.000,   70.000,   80.000,   90.000,  100.000,
+      110.000,  120.000, 130.000], 
+     [-500.000, -400.000, -300.000, -200.000, -100.000,    0.000,  100.000,
+     200.000,  300.000,  400.000,  500.000,  600.000], 
+     [-100.000,  -72.549,  -49.020,  -33.333,  -21.569,    0.000,   27.451,
+     50.980,  66.667,   78.431,  100.000,  127.451,  150.980,  166.667,
+     178.431,  200.000], 
+      [-33.333,  -29.412,  -25.490,  -21.569,   -9.804,    0.000,    7.843,
+      15.686,  27.451,   39.216,   50.980,   54.902,   58.824,   62.745, 
+      66.667,   70.588, 74.510,   78.431,   90.196,  100.000,  107.843,
+      115.686,  127.451,  139.216,  150.980], 
+      [-33.333,  -29.412,  -25.490,  -21.569,   -9.804,    0.000,    7.843, 
+      15.686, 27.451,   39.216,   50.980,   54.902,   58.824,   62.745,  
+      66.667,   70.588,  74.510,   78.431,   90.196,  100.000,  107.843,
+      115.686,  127.451,  139.216,  150.980], 
+      [-41.176,  -37.255,  -33.333,  -21.569,   -9.804,    0.000,    7.843,
+      15.686,  27.451,   39.216,   50.980,   58.824,   62.745,   66.667,
+      78.431,   90.196,  100.000,  107.843,  115.686,  127.451,  139.216,
+      150.980], 
+      [-17.647,  -13.725,   -9.804,   -5.882,   -1.961,    0.000,    3.922,
+      7.843,  11.765,  15.686,   19.608,   23.509,   27.451,   31.373,
+      35.294,   39.216,  43.137,   47.059,   50.980,   54.902,   58.824,
+      62.745,   66.667,   70.588,   74.510,   78.431,   82.353,   86.275,
+      90.196,   94.118,   98.039,  100.000,  103.922,  107.843,  111.765,
+      115.686,  119.608], 
+      [-50.000,  -40.000,  -30.000,  -20.000,  -10.000,    0.000,   10.000,
+      20.000,   30.000,   40.000,   50.000,   60.000,   70.000,   80.000,
+      90.000,  100.000,  110.000,  120.000,  130.000,  140.000,  150.000]], 
+'coef' : [[33.8253,  33.8253,  33.8253,  33.8253], 
+        [82.4859,   1.2603,  22.4815,  18.9314,   2.8969,  71.0974,  82.4859,
+        1.2603,  22.4815], 
+        [12.0792,   8.4592,  26.7900,  21.5705,  23.0754,   6.6585,  19.0962,
+        30.3922,  60.6512,  74.6763,  75.6664,  79.8148,  75.9320,  47.5211, 
+        12.0792,   8.4592,  26.7900], 
+       [13.3494,   7.8883,  25.8400,  23.0248,  21.3100,   8.9493,  15.7822,
+       34.2918,  58.3060,  72.5457,  78.4441,  79.6678,  75.9053,  46.4422, 
+       13.3494,   7.8883,  25.8400], 
+      [45.0851,  -2.1518,  27.3155,  21.7801,  22.3697,   7.9200,  16.9739,
+      33.1615,  59.0061,  80.3683,  83.5295,  45.0851,  -2.1518,  27.3155], 
+      [12.6861,   8.5744,  14.9941,  21.7397,  26.7676,  21.3221,  22.5031,
+      22.5483,  19.6945,  16.7365,  13.8898,  11.6012,  11.5439,  15.9402,
+      21.4673,  27.5536,  36.5639,  61.4126,  74.7277,  75.0275,  80.9237,
+      76.0996,  70.2252,  55.8543,  40.9132,  19.9752,  12.6861,   8.5744,
+      14.9941], 
+      [41.9234,  -3.2368,  32.0187,  19.5668,  20.2811,   9.2545,  14.9009, 
+      44.2629, 85.8093,  78.0205,  41.9234,  -3.2368,  32.0187], 
+      [33.8253,  33.8253,  33.8253,  33.8253,  33.8253,  33.8253], 
+      [93.0061,  86.0013, -54.6819,  92.8215, -49.7800,  93.0061,  86.0013,
+      -54.6819,  92.8215, -49.7800], 
+      [71.6700,  46.2404, -16.3198,  38.4125,  13.7980,  30.8425,   1.5483,
+      15.9317,  23.2525,  35.9459,  61.8592,  75.7543,  75.6340,  81.9501,
+      71.6700,  46.2404,  -16.3198,  38.4125,  13.7980], 
+      [72.7220,  44.8956, -14.3723,  36.0327,  16.7261,  27.3015,   4.5277, 
+      14.8756,  21.4604,  40.3602,  58.6143,  73.1526,  79.2995,  79.9220,
+      72.7220,  44.8956, -14.3723,  36.0327,  16.7261], 
+      [76.7467,  60.6977, -17.3393,  38.5959,  14.2538,  29.7441,   0.6054,
+      18.3639,  23.376,  64.2515,  86.7496,  76.7467,  60.6977, -17.3393, 
+      38.5959,  14.2538], 
+      [19.6912,  13.4825,   6.2608,  17.1249,  20.6067,  29.0827,  19.5270,
+      23.2109,  22.6667,  19.6866,  16.7474,  13.8240,  11.5230,  10.9217, 
+      16.1154,  21.0846,  28.1517,  33.8236,  63.1731,  76.4907,  73.0356, 
+      83.3897,  74.8104,  72.2190,  54.5815,  43.0444,  19.6912,  13.4825, 
+      6.2608,  17.1249,  20.6067], 
+      [76.0273,  49.0816, -19.3982,  45.0213,  11.9440,  25.8971,   4.7660,
+      15.1631,  39.6956,  94.5539,  76.0273,  49.0816, -19.3982,  45.0213,
+      11.9440]],
+'sp' : [[ 33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,
+33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  
+33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,
+33.825,  33.825,  33.825,  33.825],
+[ 17.610,  14.653,  15.155,  17.402,  19.674,  20.603,  20.196,  18.790,
+16.816,  14.957,  13.966,  14.600,  17.612,  23.756,  33.347,  44.926,
+56.587,  66.436,  72.800,  74.944,  72.364,  64.952,  54.201,  42.007,
+30.257,  20.848,  17.610],
+[ 10.154,  14.944,  21.971,  24.044,  23.063,  22.333,  21.684,  20.057,
+16.896,  13.521,  11.727,  12.846,  16.349,  21.241,  27.607,  39.853,
+59.205,  72.504,  76.193,  78.735,  76.496,  68.443,  55.748,  39.710, 
+22.700,  11.437,  10.154],
+[ 10.229,  14.365,  21.449,  24.103,  23.577,  22.632,  21.458,  19.571, 
+16.798,  13.995,  12.290,  12.650,  15.388,  20.651,  28.832,  41.270, 
+57.676,  71.156,  77.665,  79.088,  76.371,  68.148,  55.324,  39.540,
+23.223,  12.010,  10.229],
+[  9.743,  14.254,  21.991,  24.364,  23.262,  22.295,  21.465,  19.851,
+16.958,  13.860,  12.026,  12.636,  15.698,  20.924,  28.550,  40.879, 
+57.911,  70.869,  77.996,  79.646,  76.068,  67.642,  55.253,  39.918,
+23.566,  11.833,   9.743],
+[ 10.099,  14.835,  21.453,  25.022,  22.427,  22.315,  22.070,  19.673,
+16.754,  13.983,  11.973,  12.286,  16.129,  21.560,  28.041,  39.205,
+59.489,  72.559,  75.960,  79.137,  75.925,  68.809,  55.758,  39.915, 
+22.006,  12.076,  10.099],
+[ 10.166,  13.120,  20.713,  25.393,  24.704,  22.001,  20.333,  19.326,
+17.518,  14.652,  12.293,  12.046,  14.701,  20.681,  30.057,  42.212, 
+56.273,  69.600,  78.495,  80.503,  76.383,  67.397,  54.752,  39.615,
+23.955,  12.571,  10.166],
+[ 33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,
+33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825, 
+33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,  33.825,
+33.825,  33.825,  33.825],
+[ 17.190,  13.224,  12.783,  14.954,  18.402,  21.679,  23.510,  23.130,
+20.528,  16.633,  13.010,  11.437,  13.452,  19.920,  30.599,  43.944,
+57.555,  68.874,  75.847,  77.450,  73.684,  65.403,  54.103,  41.641, 
+29.918,  20.589,  17.189],
+[ 10.339,  14.586,  21.670,  24.495,  23.365,  21.980,  21.373,  20.137,
+17.207,  13.625,  11.619,  12.650,  16.287,  21.320,  27.802,  39.806,
+59.078,  72.423,  76.539,  78.499,  76.312,  68.654,  55.998,  39.403,
+22.514,  11.778,  10.339],
+[ 10.522,  14.367,  21.239,  24.418,  23.743,  22.313,  21.259,  19.717,
+16.995,  13.900,  12.112,  12.761,  15.686,  20.527,  28.450,  41.495,
+57.716,  71.117,  77.873,  78.930,  75.929,  68.453,  55.899,  39.343, 
+22.654,  12.043,  10.522],
+[ 10.165,  14.367,  21.641,  24.642,  23.519,  21.989,  21.225,  19.973,
+17.211,  13.835,  11.858,  12.622,  15.837,  20.810,  28.467,  41.222, 
+57.621,  71.020,  78.150,  79.291,  75.669,  68.002,  55.951,  39.793, 
+22.835,  11.746,  10.165],
+[ 10.099,  14.835,  21.453,  25.022,  22.427,  22.315,  22.070,  19.673,
+16.754,  13.983,  11.973,  12.286,  16.129,  21.560,  28.041,  39.205,
+59.489,  72.559,  75.960,  79.137,  75.925,  68.809,  55.758,  39.915, 
+22.006,  12.076,  10.099],
+[ 10.453,  13.154,  20.451,  25.336,  25.008,  22.093,  20.038,  19.155, 
+17.686,  14.895,  12.260,  11.866,  14.719,  20.778,  29.917,  42.121,
+56.387,  69.779,  78.508,  80.427,  76.147,  67.406,  55.065,  39.607,
+23.626,  12.578,  10.453]]
+}
+
+parcur_test = {
+'u' : array([120.,128.,133.,136.,138.,141.,144.,146.,149.,151.,154.,
+     161.,170.,180.,190.,200.,210.,220.,230.,240.,250.,262.,269.,
+     273.,278.,282.,287.,291.,295.,299.,305.,315.]),
+'xa' : [-1.5141,-2.0906,-1.9253,-0.8724,-0.3074,-0.5534,0.0192,
+     1.2298,2.5479,2.4710,1.7063,1.1183,0.5534,0.4727,0.3574,0.1998,
+     0.2882,0.2613,0.2652,0.2805,0.4112,0.9377,1.3527,1.5564,1.6141,
+     1.6333,1.1567,0.8109,0.2498,-0.2306,-0.7571,-1.1222],
+'xo' : [0.5150,1.3412,2.6094,3.2358,2.7401,2.7823,3.5932,
+     3.8353,2.5863,1.3105,0.6841,0.2575,0.2460,0.3689,0.2460,0.2998,
+     0.3651,0.3343,0.3881,0.4573,0.5918,0.7110,0.4035,0.0769,-0.3920,
+     -0.8570,-1.3412,-1.5641,-1.7409,-1.7178,-1.2989,-0.5572],
+'ub':120,
+'ue':320,
+'k' : [3, 3, 3, 3, 3, 3, 5, 5, 5],
+'s' : [100.0, 1.0, 0.05, 0.25, 0.25, 0.25, 0.25, 0.0, -1], 
+'ipar' : [1, 1, 1, 1, 1, 0, 0, 0, 0],
+'iopt' : [0, 1, 1, 1, 0, 0, 0, 0, -1],
+'res' : [0.559278E+02, 0.100021E+01, 0.499966E-01, 0.250216E+00, 0.249842E+00,
+                     0.250003E+00, 0.249993E+00, 0.000000E+00, 0.943099E+00],
+'err' : [-2, 0, 0, 0, 0, 0, 0, -1, 0],
+'knots' : [[120., 120., 120., 120., 320., 320., 320., 320.], 
+         [120., 120., 120., 120., 133., 138., 141., 144., 149., 154., 
+  170., 210., 250., 278., 295., 320., 320., 320., 320.], 
+         [120., 120., 120., 120., 133., 136., 138., 141., 144., 149., 
+  151., 154., 161., 170., 180., 190., 210., 250., 269., 278., 
+  287., 295., 320., 320., 320., 320.], 
+         [120., 120., 120., 120., 133., 136., 138., 141., 144., 149., 
+  151., 154., 161., 170., 180., 190., 210., 250., 269., 278., 
+  287., 295., 320., 320., 320., 320.], 
+         [120., 120., 120., 120., 133., 136., 138., 141., 144., 149., 
+  151., 154., 170., 210., 250., 278., 287., 295., 320., 320., 
+  320., 320.], 
+         [0.0000, 0.0000, 0.0000, 0.0000, 0.1197, 0.1839, 0.2232, 0.2883, 
+  0.4480, 0.6343, 0.6653, 0.6838, 0.7840, 1.0000, 1.0000, 1.0000, 
+  1.0000], 
+       [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.1197, 0.1839, 
+  0.2232, 0.4480, 0.6343, 0.6653, 0.6838, 0.7840, 1.0000, 1.0000, 
+  1.0000, 1.0000, 1.0000, 1.0000], 
+       [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.1839, 0.2232, 
+  0.2363, 0.2883, 0.3529, 0.4480, 0.5149, 0.5667, 0.6047, 0.6343, 
+  0.6420, 0.6508, 0.6596, 0.6653, 0.6674, 0.6703, 0.6740, 0.6838, 
+  0.7121, 0.7391, 0.7593, 0.7840, 0.8084, 0.8440, 0.8655, 0.8963, 
+  1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000], 
+       [0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.1250, 0.2500, 
+  0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000, 1.0000, 1.0000, 
+  1.0000, 1.0000, 1.0000]],
+'sx' : [[-1.1548, 2.7543, 0.7908,  -0.1259], 
+      [-1.5121,  -1.2336,  -3.5336,  -0.1351,  -0.6165, 0.0950, 3.1911, 
+      0.9120,   0.2207, 0.5287,  -0.2985, 2.1710, 0.4299,  -2.2076,  -0.3992], 
+      [-1.5141,  -0.8415,  -3.5384,  -0.9552,  -0.0336,  -0.8011,  -0.0719, 
+      2.8821, 2.4805, 1.1626, 1.2566, 0.3134, 0.6040, 0.1359, 0.4505,  
+      -0.0664,   1.1284, 1.8397, 1.3619,  -0.5533,  -1.4796,  -0.9382], 
+      [-1.5126,  -0.8460,  -3.6051,  -0.8411,  -0.2731,  -0.6593, 0.0111, 
+      2.6652, 2.6017, 1.1042, 1.2620, 0.3178, 0.6014, 0.1370, 0.4496,  
+      -0.0659,   1.1283, 1.8398, 1.3618,  -0.5531,  -1.4799,  -0.9379], 
+      [-1.5137,  -0.6761,  -3.7166,  -0.8442,  -0.1873,  -0.7240,  -0.0382, 
+      2.8148, 2.3950, 0.9819, 0.2063, 0.4696,  -0.1694, 1.7368, 1.5210, 
+      -0.7724, -1.2116,  -1.1237], 
+      [-1.5126,  -2.1295,  -2.5681,  -0.8519,  -0.3894, 0.0129, 4.2840,
+      1.2306,  0.2041, 0.6633, 3.2860,  -1.1886,  -1.0872], 
+      [-1.5140,  -0.2538,  -3.9579,  -2.0558,  -0.1376,  -0.7093, 3.9083,
+      3.4195,  -0.6167, 1.7279, 3.0307,  -0.5555,  -0.7361,  -1.1278], 
+      [-1.5141,  -4.8893, 3.0277,  -6.9243, 2.4051,  -1.9430, 1.8101, 0.4452, 
+  3.0970, 2.7686, 2.0549, 1.1258, 1.0620, 0.2999, 0.5931,  -0.1113, 
+  0.4195, 0.2149, 0.3078, 0.1882, 0.9922, 0.7559, 1.4062, 1.6328, 
+  1.5426, 1.8382, 1.1504, 0.6591,  -0.3452,  -0.7910,  -1.0000,  -1.1222], 
+      [-1.5140,  -1.1705,  -3.7619,  -0.4161,  -0.5726, 0.7574, 5.4849, 
+      -2.0698,  2.1028, 2.6719,  -1.5004,  -0.4675,  -1.1283]], 
+'sy' : [[2.5710, 0.4640, 0.3999,  -1.6208], 
+      [0.5161,  -1.0848, 3.3776, 2.8139, 2.7532, 4.3508, 2.1979,  -0.3601, 
+  0.7465,  -0.0769, 0.7883, 0.5575,  -3.0389,  -0.1368,  -0.6145], 
+      [0.5151, 0.5826, 1.3907, 3.5299, 2.5462, 2.5185, 4.2164, 3.5608, 
+  0.8278, 0.6202, 0.0979, 0.2647, 0.3987, 0.2091, 0.4205, 0.3339, 
+  0.9355,  -0.4440,  -1.3297,  -2.2927,  -0.4057,  -0.4095], 
+      [0.5160,  -0.4483, 2.3159, 3.1579, 2.7441, 2.6293, 4.2337, 3.2004, 
+  1.2412, 0.0951, 0.3634, 0.1744, 0.4209, 0.2015, 0.4258, 0.3310, 
+  0.9364,  -0.4443,  -1.3297,  -2.2925,  -0.4061,  -0.4092], 
+      [0.5157,  -0.0216, 1.9374, 3.3003, 2.7152, 2.4783, 4.3210, 3.3174, 
+  1.0242, 0.0591, 0.4637, 0.1370, 0.6160, 0.7280,  -1.5065,  -2.0501, 
+ -0.7021,  -0.2044], 
+       [0.5173, 0.8949, 2.5415, 3.4316, 2.4271, 4.5986, 2.6667,  -0.1322, 
+  0.3411, 0.9389,  -1.1971,  -2.5963,  -0.5047], 
+       [0.5160,  11.3627, -10.2865, 7.5268,  -0.4078, 6.4015, 2.8602, 0.5783, 
+ -0.1220, 1.5742,  -2.4089,  -1.9223,  -1.6493,  -0.5481], 
+       [0.5150, 3.0202,  -2.1849, 6.2691, 2.1376, 2.9157, 4.2898, 3.9404, 
+  2.7047, 2.0088, 0.0609, 1.2120,  -0.6203, 0.6126, 0.2469, 0.0576, 
+  0.5412, 0.2445, 0.5039, 0.3663, 1.2118, 0.4929, 0.5866, 0.0256, 
+ -0.4126,  -1.0048,  -1.2798,  -1.8261,  -1.8191,  -1.7347,  -0.6137,  
+ -0.5572], 
+       [0.5157, 0.1690, 2.5142, 3.5678, 1.7000, 5.7700, 0.8859,  -0.8397, 
+  2.4599,  -3.0704,  -1.5557,  -1.2391,  -0.5543]], 
+'sp' : [[  -1.1548,  2.5710, -0.7135,  2.3277, -0.4651,  2.1850, -0.3258,
+2.1025, -0.2370,  2.0489, -0.1097,  1.9703,  0.0107,  1.8938,  0.0872, 
+1.8441,  0.1963,  1.7711,  0.2654,  1.7236,  0.3637,  1.6539,  0.5685, 
+1.4983,  0.7840,  1.3113,  0.9646,  1.1183,  1.0882,  0.9380,  1.1601,  
+0.7672,  1.1854,  0.6032,  1.1693,  0.4427,  1.1171,  0.2830,  1.0338, 
+0.1208,  0.9247, -0.0466,  0.7670, -0.2589,  0.6649, -0.3903,  0.6040,
+-0.4685,  0.5261, -0.5696,  0.4626, -0.6535,  0.3823, -0.7624,  0.3178, 
+-0.8530,  0.2534, -0.9470,  0.1895, -1.0444,  0.0951, -1.1975, -0.0553,
+-1.4729],
+ [  -1.5121,  0.5161, -2.1158,  1.3155, -1.8367,  2.7813, -0.8954,  2.9112, 
+ -0.4539,  2.8268, -0.4538,  2.9794,  0.1608,  3.7106,  1.1126,  3.4893,  
+ 2.3766, 2.3660,  2.4950,  1.6528,  1.9671,  0.8609,  0.9699,  0.0639, 
+ 0.5592,  0.1890, 0.4118,  0.3782,  0.3610,  0.3838,  0.3472,  0.3084,
+ 0.3113,  0.2549,  0.2207, 0.2994,  0.1499,  0.4137,  0.2004,  0.5428,
+ 0.4734,  0.6320,  1.1011,  0.5682, 1.4267,  0.3319,  1.5365,  0.0978, 
+ 1.5511, -0.3202,  1.4407, -0.7531,  1.1510, -1.3044,  0.8062, -1.6447, 
+ 0.3679, -1.7989, -0.1420, -1.7090, -0.8573, -1.2727, -1.1104, -0.5605],
+ [  -1.5141,  0.5151, -2.0904,  1.3399, -1.9278,  2.6176, -0.8627,  3.2135, 
+ -0.3177,  2.7648, -0.5578,  2.7552,  0.0533,  3.6603,  1.1877,  3.7583, 
+ 2.5699,  2.6222,  2.4649,  1.3059,  1.7052,  0.6802,  1.1204,  0.2567, 
+ 0.5470,  0.2506,  0.4954,  0.3519,  0.3069,  0.2824,  0.2436,  0.2790, 
+ 0.2870,  0.3327,  0.2711,  0.3674,  0.2337,  0.3967,  0.2638,  0.4558, 
+ 0.4501,  0.5805,  0.9479,  0.6914,  1.3254,  0.4292,  1.5373,  0.0937, 
+ 1.6705, -0.4295,  1.5788, -0.8456,  1.2139, -1.3107,  0.7606, -1.6031, 
+ 0.2559, -1.7475, -0.2110, -1.6852, -0.7669, -1.3134, -1.1211, -0.5557],
+ [  -1.5126,  0.5160, -2.1059,  1.3246, -1.8862,  2.7128, -0.8361,  3.0274,
+ -0.4395,  2.8257, -0.4955,  2.8696,  0.1192,  3.6617,  1.1486,  3.6310, 
+ 2.4908,  2.5597,  2.5133,  1.5174,  1.7232,  0.5758,  1.1087,  0.2696, 
+ 0.5506,  0.2489,  0.4945,  0.3509,  0.3067,  0.2853,  0.2440,  0.2764, 
+ 0.2870,  0.3323,  0.2709,  0.3685,  0.2336,  0.3974,  0.2638,  0.4555, 
+ 0.4503,  0.5798,  0.9479,  0.6915,  1.3254,  0.4295,  1.5373,  0.0938,  
+ 1.6705, -0.4296,  1.5788, -0.8458,  1.2138, -1.3107,  0.7606, -1.6031, 
+ 0.2560, -1.7475, -0.2110, -1.6852, -0.7669, -1.3134, -1.1211, -0.5556],
+ [  -1.5137,  0.5157, -2.0941,  1.3313, -1.9242,  2.6683, -0.8241,  3.1080,
+ -0.3888,  2.8231, -0.5299,  2.7740,  0.0879,  3.6846,  1.1787,  3.7093, 
+ 2.5018,  2.5558,  2.4167,  1.4002,  1.8823,  0.6816,  1.0126,  0.2367, 
+ 0.5944,  0.2522,  0.4116,  0.3221,  0.3401,  0.3241,  0.3214,  0.3006, 
+ 0.2964,  0.2937,  0.2303,  0.3351,  0.1832,  0.4139,  0.2389,  0.5090, 
+ 0.4813,  0.5991,  1.0356,  0.5942,  1.3548,  0.3768,  1.4908,  0.1334, 
+ 1.5810, -0.3273,  1.5562, -0.8151,  1.2789, -1.3812,  0.7977, -1.6444, 
+ 0.2360, -1.7257, -0.2406, -1.6526, -0.7496, -1.3325, -1.1231, -0.5535],
+ [  -1.5126,  0.5173, -2.0962,  1.3284, -1.9154,  2.6519, -0.8579,  3.1350,
+ -0.4780,  2.7906, -0.3927,  2.8295,  0.0356,  3.5770,  1.1661,  3.8095, 
+ 2.6447,  2.5908,  2.4131,  1.3715,  1.6818,  0.5972,  1.0403,  0.2693, 
+ 0.5711,  0.2212,  0.4632,  0.2442,  0.3622,  0.2891,  0.3027,  0.3526, 
+ 0.2937,  0.4048,  0.2976,  0.4261,  0.3085,  0.4553,  0.3314,  0.4940, 
+ 0.4264,  0.5827,  0.8211,  0.5959,  1.2421,  0.3329,  1.5049,  0.0288, 
+ 1.6804, -0.3960,  1.6292, -0.8022,  1.2238, -1.3065,  0.8512, -1.5331, 
+ 0.2420, -1.7153, -0.2534, -1.7101, -0.8192, -1.4173, -1.0872, -0.5047],
+ [  -1.5140,  0.5160, -2.0913,  1.3343, -1.9224,  2.6454, -0.8413,  3.1098,
+ -0.5097,  2.8015, -0.4106,  2.8886,  0.1128,  3.5175,  1.1376,  3.8232, 
+ 2.5443,  2.6009,  2.5286,  1.3516,  1.7687,  0.6085,  0.9935,  0.3006, 
+ 0.4845,  0.2434,  0.3955,  0.2589,  0.3267,  0.2931,  0.2992,  0.3436, 
+ 0.3045,  0.3845,  0.3113,  0.4009,  0.3240,  0.4231,  0.3471,  0.4529, 
+ 0.4388,  0.5286,  0.8568,  0.6175,  1.2911,  0.4212,  1.5361,  0.1131, 
+ 1.6758, -0.3655,  1.6119, -0.8419,  1.1950, -1.3926,  0.8147, -1.6001,  
+ 0.2143, -1.7155, -0.2427, -1.6562, -0.7323, -1.3483, -1.1278, -0.5481],
+ [  -1.5141,  0.5150, -2.0906,  1.3412, -1.9253,  2.6094, -0.8724,  3.2358,
+ -0.3074,  2.7401, -0.5534,  2.7823,  0.0192,  3.5932,  1.2298,  3.8353, 
+ 2.5479,  2.5863,  2.4710,  1.3105,  1.7063,  0.6841,  1.1183,  0.2575, 
+ 0.5534,  0.2460,  0.4727,  0.3689,  0.3574,  0.2460,  0.1998,  0.2998, 
+ 0.2882,  0.3651,  0.2613,  0.3343,  0.2652,  0.3881,  0.2805,  0.4573,  
+ 0.4112,  0.5918,  0.9377,  0.7110,  1.3527,  0.4035,  1.5564,  0.0769, 
+ 1.6141, -0.3920,  1.6333, -0.8570,  1.1567, -1.3412,  0.8109, -1.5641,  
+ 0.2498, -1.7409, -0.2306, -1.7178, -0.7571, -1.2989, -1.1222, -0.5572],
+ [  -1.5140,  0.5157, -2.0927,  1.3305, -1.9113,  2.6869, -0.8937,  2.9204,
+ -0.4755,  2.9632, -0.3647,  3.0067,  0.1005,  3.3658,  1.0401,  3.7905, 
+ 2.6920,  2.7554,  2.5788,  1.2773,  1.5805,  0.4841,  0.8142,  0.2697, 
+ 0.4514,  0.2967,  0.4047,  0.3224,  0.3764,  0.3570,  0.3746,  0.3937, 
+ 0.3870,  0.4178,  0.3943,  0.4265,  0.4060,  0.4376,  0.4249,  0.4515,  
+ 0.4932,  0.4825,  0.8026,  0.4927,  1.1681,  0.3386,  1.4146,  0.1027, 
+ 1.6128, -0.3052,  1.6351, -0.7716,  1.2950, -1.3882,  0.8918, -1.6327, 
+ 0.1869, -1.7454, -0.3242, -1.6509, -0.7005, -1.3275, -1.1283, -0.5543]]
+}
+
+clocur_test = {
+'xa' : [-4.7,-7.048,-6.894,-3.75,-1.042,0.938,2.5,3.524,4.511,5.0,4.886,
+           3.524,3.2,1.302,-1.424,-3.0,-3.064,-3.665, -4.7],
+'xo' : [0.0,2.565,5.785,6.495,5.909,5.318,4.33,2.957,1.642,0.0,-1.779,
+            -2.957,-5.543,-7.386,-8.075,-5.196,-2.571,-1.334, 0.0],
+'u' : arange(19.0)*20.0,
+'k' : [3, 3, 3, 3, 3, 3, 5, 5, 5],
+'ipar' : [1, 1, 1, 1, 1, 0, 0, 0, 0],
+'iopt' : [0, 1, 1, 1, 0, 0, 0, 0, -1],
+'s' : [900., 10., 0.1, 0.5, 0.5, 0.5, 0.5, 0.0, -1],
+'res' : [0.653304E+03, 0.100003E+02, 0.999982E-01, 0.499976E+00, 0.499666E+00,
+       0.499997E+00, 0.500067E+00, 0.000000E+00, 0.207778E+01],
+'err' : [-2, 0, 0, 0, 0, 0, 0, -1, 0],
+'knots' : [[ -1080.,  -720.,  -360.,  0., 360., 720.,  1080.,  1440. ], 
+        [-180., -80., -40.,   0.,  40.,  60., 100., 180., 280., 320., 
+          360., 400., 420., 460.], 
+        [-60., -40., -20.,   0.,  20.,  40.,  60., 100., 180., 200., 
+          220., 240., 260., 280., 300., 320., 340., 360., 380., 400., 
+          420.], 
+        [-60., -40., -20.,   0.,  20.,  40.,  60., 100., 180., 200., 
+         220., 240., 260., 280., 300., 320., 340., 360., 380., 400., 
+         420.], 
+        [-120., -80., -40.,   0.,  40.,  60., 100., 180., 200., 220., 
+         240., 280., 320., 360., 400., 420., 460.], 
+        [-0.3411, -0.2121, -0.0724,  0.0000,  0.0822,  0.1584,  0.2346, 
+        0.3490, 0.5125,  0.5547,  0.5973,  0.6589,  0.7879,  0.9276,  1.0000,
+        1.0822, 1.1584,  1.2346], 
+        [-0.4453, -0.4027, -0.3411, -0.2121, -0.0724,  0.0000,  0.1584,  
+        0.2346,  0.3490,  0.5125,  0.5547,  0.5973,  0.6589,  0.7879, 
+        0.9276,  1.0000,  1.1584,  1.2346,  1.3490,  1.5125,  1.5547], 
+        [-0.2786, -0.2121, -0.1345, -0.0724, -0.0399,  0.0000,  0.0822, 
+        0.1584,  0.2346,  0.3001,  0.3490,  0.3927,  0.4332,  0.4720, 
+        0.5125,  0.5547, 0.5973,  0.6589,  0.7214,  0.7879,  0.8655, 
+        0.9276,  0.9601,  1.0000, 1.0822,  1.1584,  1.2346,  1.3001, 
+        1.3490], 
+        [-0.6250, -0.5000, -0.3750, -0.2500, -0.1250,  0.0000,  0.1250,
+        0.2500,  0.3750,  0.5000,  0.6250,  0.7500,  0.8750,  1.0000, 
+        1.1250,  1.2500, 1.3750,  1.5000,  1.6250]], 
+'sx' : [[ -0.2890, -0.2890, -0.2890, -0.2890], 
+     [-3.6414, -5.1410, -7.1471, -3.6077,  2.9569,  7.0497,  0.8157, -3.6414, 
+      -5.1410, -7.1471], 
+     [-3.6393, -4.2795, -7.6104, -7.5242, -2.2356,  2.5034,  4.8626,  5.2229, 
+      3.1970,  3.5381,  1.5313, -1.7123, -3.3029, -2.9197, -3.6393, -4.2795, 
+      -7.6104], 
+     [-3.3998, -4.4768, -7.5699, -7.4089, -2.3288,  2.4959,  5.0209,  4.9720, 
+      3.6189,  3.2870,  1.5156, -1.6335, -3.2388, -3.1205, -3.3998, -4.4768, 
+      -7.5699], 
+     [-3.1406, -3.8616, -9.3043, -2.0703,  2.3649,  4.9458,  5.1866,  3.1343, 
+      3.9923, -2.7579, -3.1406, -3.8616, -9.3043], 
+     [-2.7754, -4.3421, -7.6213, -7.4935, -2.9934,  1.5161,  5.3675,  4.7346, 
+      3.5203,  3.6277, -2.8775, -2.7754, -4.3421, -7.6213], 
+     [-5.1336, -1.1870, -4.9707,  -10.2635, -1.8832,  1.1682,  4.6332,  
+     5.9054, 2.4087,  4.9227, -5.1336, -1.1870, -4.9707,  -10.2635, -1.8832], 
+     [-2.3401, -3.8401, -5.2545, -7.7951, -8.2051, -3.8222, -1.6092,  0.6087, 
+      2.7131,  3.2464,  4.9215,  4.8089,  5.8715,  1.6678,  4.7786,  0.6271, 
+      -1.5750, -3.8446, -2.3401, -3.8401, -5.2545, -7.7951, -8.2051], 
+     [1.0044, -5.4019, -0.9630,  -12.7083, -1.1542,  1.1652,  6.9462, 
+     3.4423, 1.0044, -5.4019, -0.9630,  -12.7083, -1.1542]], 
+'sy' : [[ 0.0089,  0.0089,  0.0089,  0.0089],
+     [-2.8481, -0.2362,  5.1680,  7.0915,  4.1090,  1.8202,  -10.7190, 
+     -2.8481,  -0.2362,  5.1680], 
+     [-1.3506, -0.3010,  2.3785,  6.4303,  6.5434,  4.9464,  1.8911, -1.9032, 
+      -2.6544, -5.6451, -7.5794, -8.8289, -5.1023, -2.2260, -1.3506, -0.3010, 
+      2.3785], 
+     [-1.2300, -0.5075,  2.6980,  6.1354,  6.7093,  4.8744,  1.8243, -1.7559, 
+      -2.9185, -5.3807, -7.8932, -8.4817, -5.3102, -2.2288, -1.2300, -0.5075, 
+      2.6980], 
+     [-0.6661, -1.9633,  5.9976,  6.6889,  4.8408,  1.9295, -1.8934, -2.7464, 
+      -6.3406,  -10.2291, -0.6661, -1.9633,  5.9976], 
+     [-2.8663,  0.0279,  2.2083,  6.4113,  6.5499,  5.6375,  1.6714, -1.6993, 
+      -3.2769, -6.4064,  -10.0517, -2.8663,  0.0279,  2.2083], 
+     [-11.4885, -0.4907, -1.3396,  6.4019,  7.1875,  5.2912,  3.0296, 
+     -1.2442,-3.8350, -7.0871,  -11.4885, -0.4907, -1.3396,  6.4019,  
+     7.1875], 
+     [-3.3543, -0.9872,  0.9217,  1.8210,  7.0937,  6.6891,  6.0142,  5.5460, 
+      4.6549,  2.8892,  1.7125,  0.0930, -2.5059, -2.4784, -6.6540, -7.0488, 
+      -9.7840, -5.4089, -3.3543, -0.9872,  0.9217,  1.8210,  7.0937], 
+     [-12.0940, -3.6528, -0.7300,  5.5358,  7.2803,  5.8203, -0.2992, 
+     -2.7102, -12.0940, -3.6528, -0.7300,  5.5358,  7.2803]] ,
+'sp' : [[  -0.2890,  0.0089, -0.2890,  0.0089, -0.2890,  0.0089, -0.2890,  
+        0.0089, -0.2890,  0.0089, -0.2890,  0.0089, -0.2890,  0.0089, 
+        -0.2890, 0.0089, -0.2890,  0.0089, -0.2890,  0.0089, -0.2890, 
+        0.0089, -0.2890, 0.0089, -0.2890,  0.0089, -0.2890,  0.0089, 
+        -0.2890,  0.0089, -0.2890, 0.0089, -0.2890,  0.0089, -0.2890,  0.0089],
+     [  -5.2923,  0.4093, -6.1787,  2.9710, -6.0695,  5.3207, -4.2389,  
+     6.4366,  -1.7409,  6.1224,  0.7042,  5.1065,  2.6665,  4.1203,  4.0480,
+     3.1423,  4.8338,  1.9627,  5.0087,  0.3716,  4.5831, -1.7274,  3.6691, 
+     -3.9781,      2.4042, -5.9109,  0.9256, -7.0561, -0.6291, -6.9443,
+     -2.1199, -5.3828,    -3.3961, -3.2874, -4.3795, -1.5164],
+     [  -4.7279, -0.0293, -7.0409,  2.6072, -6.8775,  5.7691, -3.7728, 
+     6.4296, -1.0175,  6.0270,  0.9227,  5.2638,  2.4441,  4.2658,  3.6274, 
+     3.0527,  4.4745,  1.6179,  4.9874, -0.0452,  4.8553, -1.7122,  3.5915,
+     -3.0276,  3.1468, -5.4690,  1.3252, -7.4653, -1.4368, -7.9995, -2.9739,
+     -5.2440,  -3.1035, -2.5595, -3.6261, -1.3216],
+     [  -4.8128, -0.0937, -7.0276,  2.7367, -6.8007,  5.6343, -3.7924,  
+     6.4306,  -1.0745,  6.1030,  0.8984,  5.2829,  2.4686,  4.2353,  3.6817,
+     3.0092,  4.5092,  1.5998,  4.9227,  0.0020,  4.7506, -1.6513,  3.7891,
+     -3.1351,  3.0471, -5.3891,  1.2860, -7.5725, -1.3762, -7.8550, -2.9516,
+     -5.3252, -3.1868, -2.5759, -3.5327, -1.2761],
+     [  -4.8300, -0.1549, -6.8257,  2.7970, -7.0124,  5.6512, -3.7882,  
+     6.4165, -0.9607,  6.0812,  0.9211,  5.2675,  2.4192,  4.2429,  3.6200, 
+     3.0380,  4.4912,  1.6245,  5.0005, -0.0262,  4.8245, -1.7170,  3.5836, 
+     -3.0535,  3.2563, -5.4018,  1.1922, -7.5794, -1.4717, -7.8576, -2.7955,
+     -5.3774, -3.1970, -2.4761, -3.6292, -1.3149],
+     [  -4.5913, -0.0815, -7.0773,  2.6262, -6.8714,  5.7487, -3.7630, 
+     6.4357, -1.0052,  6.0385,  0.8565,  5.3025,  2.4386,  4.2481,  3.6968,  
+     2.9886, 4.5577,  1.5926,  4.9321,  0.0128,  4.6151, -1.6349,  3.8307, 
+     -3.1175,  3.1953, -5.4544,  1.1602, -7.4783, -1.4290, -7.9455, -2.7722,
+     -5.3452,  -3.2724, -2.4829, -3.7031, -1.2936],
+     [  -4.6787, -0.2073, -7.0474,  2.8389, -6.8782,  5.5855, -3.8179,  
+     6.5057, -0.9262,  6.0491,  0.8978,  5.2475,  2.3835,  4.2394,  3.6211,  
+     3.0351,  4.5499,  1.6308,  5.0208, -0.0120,  4.6762, -1.6795,  3.7408, 
+     -3.1379,  3.0946, -5.3511,  1.3669, -7.5345, -1.5010, -8.0181, -2.8910,
+     -5.2226,      -3.1657, -2.4427, -3.6475, -1.3664],
+     [  -4.7000,  0.0000, -7.0480,  2.5650, -6.8940,  5.7850, -3.7500,  
+     6.4950,  -1.0420,  5.9090,  0.9380,  5.3180,  2.5000,  4.3300,  3.5240, 
+     2.9570,  4.5110,  1.6420,  5.0000,  0.0000,  4.8860, -1.7790,  3.5240,
+     -2.9570,  3.2000, -5.5430,  1.3020, -7.3860, -1.4240, -8.0750, -3.0000,
+     -5.1960,  -3.0640, -2.5710, -3.6650, -1.3340],
+     [  -4.4547, -0.0336, -6.9471,  3.0553, -7.0308,  5.3645, -3.8468, 
+     6.4116, -0.8360,  6.2133,  0.9573,  5.3963,  2.3828,  4.2137,  3.5767, 
+     2.8525,  4.4535,  1.4584,  4.8942,  0.0339,  4.7651, -1.4260,  4.1475,
+     -3.0152,  2.7862, -5.6052,  1.0706, -7.6954, -1.0414, -7.7248, -2.9141,
+     -5.1502, -3.4430, -2.6947, -3.7218, -1.4943]]
+}
+        
+
 def f1(x,d=0):
     if d is None: return "sin"
     if x is None: return "sin(x)"

Modified: trunk/Lib/sandbox/spline/tests/test_fitpack.py
===================================================================
--- trunk/Lib/sandbox/spline/tests/test_fitpack.py	2007-04-08 15:27:44 UTC (rev 2901)
+++ trunk/Lib/sandbox/spline/tests/test_fitpack.py	2007-04-10 22:04:33 UTC (rev 2902)
@@ -13,10 +13,10 @@
 
 import sys
 from numpy.testing import *
-from numpy import array, arange, around, pi, sin, ravel
+from numpy import array, arange, around, pi, sin, ravel, zeros, asarray
 
 set_package_path()
-from spline.fitpack import splrep, splev, sproot, splint, spalde
+from spline.fitpack import splprep, splrep, splev, sproot, splint, spalde
 from spline.fitpack import bisplev, bisplrep, splprep
 restore_path()
 
@@ -24,7 +24,7 @@
 from dierckx_test_data import *
 restore_path()
 
-class test_splrep_slev(NumpyTestCase):
+class test_splrep_splev(NumpyTestCase):
     def check_curfit_against_dierckx_smth(self):
         x,y = curfit_test['x'],curfit_test['y']
         k,s = curfit_test_smth['k'],curfit_test_smth['s']
@@ -60,6 +60,101 @@
             assert_array_almost_equal(around(sp,1),
                                       curfit_test_lsq['sp'][i])
 
+    def check_percur_against_dierckx(self):
+            x,y = percur_test['x'], percur_test['y']
+            k,s = percur_test['k'], percur_test['s']
+            iopt, res = percur_test['iopt'], percur_test['res']
+            err = percur_test['err']
+            coef, knots = percur_test['coef'], percur_test['knots']
+            sp = percur_test['sp']
+            for i in range(len(k)):
+                if iopt[i] != -1:
+                    tck,fp,ier,msg = splrep(x,y,k=k[i],task=iopt[i],s=s[i],
+                                                    per=True,full_output=True)
+                else:
+                    tck,fp,ier,msg = splrep(x,y,t=knots[i],k=k[i],task=iopt[i],
+                                                     per=True,full_output=True)
+                tt,cc,kk = tck
+                tt,cc = asarray(tt), asarray(cc)
+                assert_almost_equal(ier,err[i])
+                assert_almost_equal(fp,res[i],decimal=1)
+                assert_array_almost_equal(tt,knots[i], decimal=3)
+                assert_array_almost_equal(cc,coef[i], decimal=3)
+                yy = asarray(splev(x,tck))
+                assert_array_almost_equal(yy,sp[i], decimal=3)
+            
+class test_splprep_splev(NumpyTestCase):
+    def check_parcur_against_dierckx(self):
+        xa,xo = parcur_test['xa'], parcur_test['xo']
+        k,s = parcur_test['k'], parcur_test['s']
+        u = parcur_test['u']
+        ub,ue = parcur_test['ub'], parcur_test['ue']
+        iopt, res = parcur_test['iopt'], parcur_test['res']
+        err, ipar = parcur_test['err'], parcur_test['ipar']
+        knots = parcur_test['knots']
+        sx, sy = parcur_test['sx'], parcur_test['sy']
+        sp = parcur_test['sp']
+        x = array([xa, xo])
+        for i in range(len(k)):
+            if iopt[i] != -1:
+                if ipar[i] == 1:
+                    tcku,fp,ier,msg = splprep(x,u=u,ub=ub,ue=ue,k=k[i],
+                                   task=iopt[i],s=s[i],full_output=True)
+                else:
+                    tcku,fp,ier,msg = splprep(x,ub=ub,ue=ue,k=k[i],
+                                   task=iopt[i],s=s[i],full_output=True)
+            else:
+                tcku,fp,ier,msg = splprep(x,ub=ub,ue=ue,t=knots[i],
+                                   k=k[i],task=iopt[i],full_output=True)
+            tck,u = tcku
+            tt,cc,kk = tck
+            tt,cc = asarray(tt), asarray(cc)
+            assert_almost_equal(ier,err[i])
+            assert_almost_equal(fp,res[i],decimal=3)
+            assert_array_almost_equal(tt,knots[i], decimal=3)
+            assert_array_almost_equal(cc[0],sx[i], decimal=3)
+            assert_array_almost_equal(cc[1],sy[i], decimal=3)
+            y = asarray(splev(u,tck))
+            yy = zeros(64, 'float')
+            yy[0:-1:2] = y[0]
+            yy[1::2] = y[1]
+            assert_array_almost_equal(yy,sp[i], decimal=3)
+    
+    def check_clocur_against_dierckx(self):
+        xa,xo = clocur_test['xa'], clocur_test['xo']
+        k,s = clocur_test['k'], clocur_test['s']
+        u = clocur_test['u']
+        iopt, res = clocur_test['iopt'], clocur_test['res']
+        err, ipar = clocur_test['err'], clocur_test['ipar']
+        knots = clocur_test['knots']
+        sx, sy = clocur_test['sx'], clocur_test['sy']
+        sp = clocur_test['sp']
+        x = array([xa, xo]) 
+        for i in range(len(k)):
+            if iopt[i] != -1:
+                if ipar[i] == 1:
+                    tcku,fp,ier,msg = splprep(x,u=u,k=k[i],task=iopt[i],
+                                        s=s[i],per=True,full_output=True)
+                else:
+                    tcku,fp,ier,msg = splprep(x,k=k[i],task=iopt[i],
+                                        s=s[i],per=True,full_output=True)
+            else:
+                tcku,fp,ier,msg = splprep(x,t=knots[i],k=k[i],task=iopt[i],
+                                                 per=True,full_output=True)
+            tck,u = tcku
+            tt,cc,kk = tck
+            tt,cc = asarray(tt), asarray(cc)
+            assert_almost_equal(ier,err[i])
+            assert_almost_equal(fp,res[i],decimal=3)
+            assert_array_almost_equal(tt,knots[i], decimal=3)
+            assert_array_almost_equal(cc[0],sx[i], decimal=3)
+            assert_array_almost_equal(cc[1],sy[i], decimal=3)
+            y = asarray(splev(u,tck))
+            yy = zeros(36, 'float')
+            yy[0:-1:2] = y[0,:-1]
+            yy[1::2] = y[1,:-1]
+            assert_array_almost_equal(yy,sp[i], decimal=3)
+
 class test_splint_spalde(NumpyTestCase):
     def check_splint_spalde(self):
         per = [0, 1, 0]
@@ -139,9 +234,9 @@
                 tck=splrep(x,v,s=s,per=per,k=k)
                 uv=splev(dx,tckp)
                 assert_almost_equal(0.0, around(abs(uv[1]-f(uv[0])),2), 
-                                                                     decimal=1)
+                                                                    decimal=1)
                 assert_almost_equal(0.0, 
-                            around(abs(splev(uv[0],tck)-f(uv[0])),2),decimal=1)
+                        around(abs(splev(uv[0],tck)-f(uv[0])),2),decimal=1)
 
 if __name__ == "__main__":
-    NumpyTest().run()
\ No newline at end of file
+    NumpyTest().run()



More information about the Scipy-svn mailing list