[SciPy-dev] Periodic spline interpolation bug / memory error?

[Zachary Pincus]

Hi folks,

I'm trying to estimate smooth contours with scipy's parametric  
splines, given a set of (x, y) points known to lie along the contour.

Now, because contours are closed, special care must be taken to  
ensure that the interpolated value at the last point be the same as  
the interpolated value at the first point.

I assume that the 'per' option to scipy.interpolate.splprep is  
designed to allow for this, as per the documentation:
> per -- If non-zero, data points are considered periodic with period
>              x[m-1] - x[0] and a smooth periodic spline  
> approximation is returned.
>              Values of y[m-1] and w[m-1] are not used.

However, I cannot get this to work on my computer. Setting 'per =  
true' always results in memory errors or other problems.

Here is a simple example to reproduce the problem:
# first make x and y points along the unit circle, from zero to just  
below two pi
In [1]: import numpy, scipy, scipy.interpolate
In [2]: twopi = numpy.arange(0, 2 * numpy.pi, 0.1)
In [3]: xs = numpy.cos(twopi)
In [4]: ys = numpy.sin(twopi)
In [5]: tck, uout = scipy.interpolate.splprep([xs, ys], u = twopi,  
per = True)
Warning: Setting x[0][63]=x[0][0]
Warning: Setting x[1][63]=x[1][0]
...[here my machine grinds for 2-3 minutes]...
Warning:     The required storage space exceeds the available  
strorage space.
     Probably causes: nest to small or s is too small. (fp>s)


At this point, the returned tck arrays are just all zeros. Sometimes  
I get other malloc errors printed to stdout and memory error  
exceptions, e.g.:
Python(6820,0xa000ed68) malloc: *** vm_allocate(size=3716243456)  
failed (error code=3)
Python(6820,0xa000ed68) malloc: *** error: can't allocate region
Python(6820,0xa000ed68) malloc: *** set a breakpoint in szone_error  
to debug

During the time that my machine is grinding, python is using very  
little CPU -- the grind is all because python is allocating huge  
amounts of memory, causing the pager to go nuts.

If I explicitly make the last value of the x and y input arrays equal  
to the first value (as the warnings say that the function is doing),  
I get the same problem:
In [6]: xs[-1] = xs[0]
In [7]: ys[-1] = ys[0]
In [8]: tck, uout = scipy.interpolate.splprep([xs, ys], u = twopi,  
per = True)
#same thing

Any thoughts?

Zach Pincus

Program in Biomedical Informatics and Department of Biochemistry
Stanford University School of Medicine