[SciPy-User] Seg. fault from scipy.interpolate or numpy

[Bruce Southey]

On Fri, Mar 25, 2011 at 3:17 PM, James Turner <jturner at gemini.edu> wrote:
> I'm getting a segmentation fault with no traceback from my Python
> script, even after upgrading to the latest stable NumPy 1.5.1 and
> SciPy 0.9.0. I have managed to whittle the code down to 30 lines
> of Python in a single file that still reproduces the error, with
> only NumPy, SciPy and PyFITS as dependencies (previously there was
> a whole chain of PyRAF scripts). However, my test script still
> requires a 109Mb input file.
>
> Would someone mind pointing me in the right direction either to
> report this properly or troubleshoot it further, please? Should I
> make a Trac ticket on the SciPy developers' page? I can put the
> input file on an sftp server if that helps (and perhaps reduce it
> to 40Mb). I have figured out how to run gdb on my script (thanks
> to STScI) and am pasting the output from "run" and "where" below,
> but my familiarity with the debugger is otherwise very limited.
>
> As you can see from gdb, the crash occurs in numpy, but I suspect
> the cause of the problem may be the array that is produced at the
> previous line by interp.splev(). Printing **des in gdb shows some
> funny characters and an apparently-huge reference count, but I
> don't really know what I'm doing / looking for.
>
> Thanks a lot!
>
> James.
>
[snip]
>
> import pyfits
> import numpy as np
> import scipy.interpolate as interp
>
> input = pyfits.open("brgN20080607S0485.fits")[1].data
>
> gain = 1.0
> readn = 3.5
>
> xdim = input.shape[1]
> ydim = input.shape[0]
> xorder = 2*xdim/4.3
> yorder = 2*ydim/4.5
> imgx = range(xdim)
> imgy = range(ydim)
> xstep = (xdim-1) / xorder
> ystep = (ydim-1) / yorder
> knots = np.arange(xstep, xdim-1-xstep, xstep)
>
> oldoutput = input.copy()
> objfit = oldoutput.copy()
>
> for yc in range(ydim):
>     tck = interp.splrep(imgx, oldoutput[yc], k=3, t=knots)
>     objfit[yc] = interp.splev(imgx, tck=tck)
>
> noise = np.sqrt(objfit.clip(min=0.0)*gain + readn*readn) / gain
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user
>
[I did see the ticket comments and it doesn't crash for me on Linux
and Python2.7.]

Why do you use 'input' for as the variable name?
This could be a problem because input is a builtin function. If the
error goes away then it is most likely a Python bug.

Second, is there a reason for using copy() for 'objfit' instead of one of these?
    zeros_like : Return an array of zeros with shape and type of input.
    ones_like : Return an array of ones with shape and type of input.
    empty_like : Return an empty array with shape and type of input.
    ones : Return a new array setting values to one.
    empty : Return a new uninitialized array.

So can you determine exactly which line it crashes on (especially the
value yc) and the shapes of the arrays?

Bruce