[Scipy-svn] r2902 - in trunk/Lib/sandbox: rbf rbf/tests spline spline/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Tue Apr 10 18:06:18 EDT 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 at 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 at tamu.edu>
+Copyright (c) 2007, John Travers <jtravs at 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