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

[James Turner]

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.

------

(gdb) run scripts/test.py
Starting program: /astro/iraf/i686/gempylocal/bin/python scripts/test.py
[Thread debugging using libthread_db enabled]
[New Thread 0x5555ae10 (LWP 27190)]

Program received signal SIGSEGV, Segmentation fault.
0x55a811b8 in _update_descr_and_dimensions (des=0xffecee74,
     newdims=0xffeced50, newstrides=0x0, oldnd=0, isfortran=0)
     at numpy/core/src/multiarray/ctors.c:222
222     numpy/core/src/multiarray/ctors.c: No such file or directory.
         in numpy/core/src/multiarray/ctors.c
(gdb) where
#0  0x55a811b8 in _update_descr_and_dimensions (des=0xffecee74,
     newdims=0xffeced50, newstrides=0x0, oldnd=0, isfortran=0)
     at numpy/core/src/multiarray/ctors.c:222
#1  0x55a4b869 in PyArray_NewFromDescr (subtype=0x55aa6940, descr=0x9196920,
     nd=0, dims=0x0, strides=0x0, data=0x0, flags=0, obj=0x0)
     at numpy/core/src/multiarray/ctors.c:1401
#2  0x55a81ac9 in Array_FromPyScalar (op=0x8d08fdc, typecode=0x55aa6940)
     at numpy/core/src/multiarray/ctors.c:1027
#3  0x55a4c066 in PyArray_FromAny (op=0x8d08fdc, newtype=0x9196920,
     min_depth=0, max_depth=0, flags=1281, context=0x0)
     at numpy/core/src/multiarray/ctors.c:1833
#4  0x55a5ddbf in PyArray_Clip (self=0x91144d0, min=0x8d08fdc, max=0x0,
     out=0x0) at numpy/core/src/multiarray/calculation.c:894
#5  0x55a8b9fb in array_clip (self=0x91144d0, args=0x5555d02c, kwds=0x9205934)
     at numpy/core/src/multiarray/methods.c:1970
#6  0x08113847 in PyCFunction_Call (func=0x91e762c, arg=0x5555d02c,
     kw=0x9205934) at Objects/methodobject.c:77
#7  0x0805e227 in PyObject_Call (func=0x91e762c, arg=0x5555d02c, kw=0x9205934)
     at Objects/abstract.c:1860
#8  0x080c4b67 in do_call (func=0x91e762c, pp_stack=0xffecf118, na=0, nk=1)
     at Python/ceval.c:3775
#9  0x080c45d8 in call_function (pp_stack=0xffecf118, oparg=256)
     at Python/ceval.c:3587
---Type <return> to continue, or q <return> to quit---
#10 0x080c19de in PyEval_EvalFrameEx (f=0x8d28664, throwflag=0)
     at Python/ceval.c:2267
#11 0x080c2d49 in PyEval_EvalCodeEx (co=0x555ae1d0, globals=0x55576acc,
     locals=0x55576acc, args=0x0, argcount=0, kws=0x0, kwcount=0, defs=0x0,
     defcount=0, closure=0x0) at Python/ceval.c:2831
#12 0x080be00f in PyEval_EvalCode (co=0x555ae1d0, globals=0x55576acc,
     locals=0x55576acc) at Python/ceval.c:494
#13 0x080e4a73 in run_mod (mod=0x8d2c6b0,
     filename=0xffed1495 "scripts/test.py", globals=0x55576acc,
     locals=0x55576acc, flags=0xffecf380, arena=0x8cd7a20)
     at Python/pythonrun.c:1271
#14 0x080e4a17 in PyRun_FileExFlags (fp=0x8cd2008,
     filename=0xffed1495 "scripts/test.py", start=257, globals=0x55576acc,
     locals=0x55576acc, closeit=1, flags=0xffecf380) at Python/pythonrun.c:1257
#15 0x080e3af2 in PyRun_SimpleFileExFlags (fp=0x8cd2008,
     filename=0xffed1495 "scripts/test.py", closeit=1, flags=0xffecf380)
     at Python/pythonrun.c:877
#16 0x080e3482 in PyRun_AnyFileExFlags (fp=0x8cd2008,
     filename=0xffed1495 "scripts/test.py", closeit=1, flags=0xffecf380)
     at Python/pythonrun.c:696
#17 0x08057269 in Py_Main (argc=2, argv=0xffecf474) at Modules/main.c:523
#18 0x0805666e in main (argc=2, argv=0xffecf474) at Modules/python.c:23
(gdb)

-----

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