scipy.cluster.hierarchy.single segmentation fault with 2**16 datapoints

[timbalam]

Hi,

Running the following code in the interactive prompt executes successfully:

>>> import scipy.cluster, numpy.random
>>> scipy.cluster.hierarchy.single(numpy.random.random_sample((2**15, 10)))

The following code results in a segmentation fault:

>>> import scipy.cluster, numpy.random
>>> scipy.cluster.hierarchy.single(numpy.random.random_sample((2**16, 10)))

Thanks,
Tim

gdb output:

GNU gdb (Ubuntu/Linaro 7.4-2012.04-0ubuntu2.1) 7.4-2012.04
Copyright (C) 2012 Free Software Foundation, Inc.
License GPLv3+: GNU GPL version 3 or later <http://gnu.org/licenses/gpl.html>
This is free software: you are free to change and redistribute it.
There is NO WARRANTY, to the extent permitted by law.  Type "show copying"
and "show warranty" for details.
This GDB was configured as "x86_64-linux-gnu".
For bug reporting instructions, please see:
<http://bugs.launchpad.net/gdb-linaro/>...
Reading symbols from /srv/sw/python/2.7.4/bin/python...done.
(gdb) run -c "import scipy.cluster, numpy.random; scipy.cluster.hierarchy.single(num
py.random.random_sample((2**16,10)))"
Starting program: /srv/sw/python/2.7.4/bin/python -c "import scipy.cluster, numpy.ra
ndom; scipy.cluster.hierarchy.single(numpy.random.random_sample((2**16,10)))"
warning: no loadable sections found in added symbol-file system-supplied DSO at 0x7f
fff7ffa000
[Thread debugging using libthread_db enabled]
Using host libthread_db library "/lib/x86_64-linux-gnu/libthread_db.so.1".

Program received signal SIGSEGV, Segmentation fault.
__pyx_pf_5scipy_7cluster_10_hierarchy_26slink (__pyx_v_n=65536, __pyx_v_Z=...,
    __pyx_v_dists=..., __pyx_self=<optimised out>)
    at scipy/cluster/_hierarchy.c:10392
10392   scipy/cluster/_hierarchy.c: No such file or directory.

gdb backtrace:

(gdb) bt
#0  __pyx_pf_5scipy_7cluster_10_hierarchy_26slink (__pyx_v_n=65536, __pyx_v_Z=...,
    __pyx_v_dists=..., __pyx_self=<optimised out>)
    at scipy/cluster/_hierarchy.c:10392
#1  __pyx_pw_5scipy_7cluster_10_hierarchy_27slink (__pyx_self=<optimised out>,
    __pyx_args=<optimised out>, __pyx_kwds=<optimised out>)
    at scipy/cluster/_hierarchy.c:10160
#2  0x000000000049ec5d in call_function (oparg=<optimised out>,
    pp_stack=0x7fffffffd530) at Python/ceval.c:4021
#3  PyEval_EvalFrameEx (f=<optimised out>, throwflag=<optimised out>)
    at Python/ceval.c:2666
#4  0x00000000004a0806 in PyEval_EvalCodeEx (co=<optimised out>,
    globals=<optimised out>, locals=<optimised out>, args=<optimised out>,
    argcount=1, kws=0xd1dc40, kwcount=2, defs=0xf525c0, defcount=2, closure=0x0)
    at Python/ceval.c:3253
#5  0x000000000049e6e5 in fast_function (nk=<optimised out>, na=1,
    n=<optimised out>, pp_stack=0x7fffffffd720, func=0xfc2cf8)
    at Python/ceval.c:4117
#6  call_function (oparg=<optimised out>, pp_stack=0x7fffffffd720)
    at Python/ceval.c:4042
#7  PyEval_EvalFrameEx (f=<optimised out>, throwflag=<optimised out>)
    at Python/ceval.c:2666
#8  0x00000000004a0806 in PyEval_EvalCodeEx (co=<optimised out>,
    globals=<optimised out>, locals=<optimised out>, args=<optimised out>,
    argcount=1, kws=0x85c568, kwcount=0, defs=0x0, defcount=0, closure=0x0)
    at Python/ceval.c:3253
#9  0x000000000049e6e5 in fast_function (nk=<optimised out>, na=1,
    n=<optimised out>, pp_stack=0x7fffffffd910, func=0xfc29b0)
    at Python/ceval.c:4117
#10 call_function (oparg=<optimised out>, pp_stack=0x7fffffffd910)
    at Python/ceval.c:4042
#11 PyEval_EvalFrameEx (f=<optimised out>, throwflag=<optimised out>)
    at Python/ceval.c:2666
#12 0x00000000004a0806 in PyEval_EvalCodeEx (co=<optimised out>,
    globals=<optimised out>, locals=<optimised out>, args=<optimised out>,
    argcount=0, kws=0x0, kwcount=0, defs=0x0, defcount=0, closure=0x0)
    at Python/ceval.c:3253
#13 0x00000000004a0932 in PyEval_EvalCode (co=<optimised out>,
    globals=<optimised out>, locals=<optimised out>) at Python/ceval.c:667
#14 0x00000000004c5b2a in run_mod (arena=0x85ddb0, flags=<optimised out>,
    locals=0x7f65f0, globals=0x7f65f0, filename=0x51e326 "<string>",
    mod=<optimised out>) at Python/pythonrun.c:1365
#15 PyRun_StringFlags (flags=<optimised out>, locals=0x7f65f0, globals=0x7f65f0,
    start=257, str=<optimised out>) at Python/pythonrun.c:1328
#16 PyRun_SimpleStringFlags (command=<optimised out>, flags=<optimised out>)
    at Python/pythonrun.c:969
#17 0x00000000004152d5 in Py_Main (argc=3, argv=<optimised out>)
    at Modules/main.c:584
#18 0x00007ffff711d76d in __libc_start_main (main=0x4149e0 <main>, argc=3,
    ubp_av=0x7fffffffdc98, init=<optimised out>, fini=<optimised out>,
    rtld_fini=<optimised out>, stack_end=0x7fffffffdc88) at libc-start.c:226
#19 0x0000000000414a11 in _start ()

[argriffing]

This segfaults for me too. Is it related to #2089?

[timbalam]

Hi argriffing,

The segfault seems specific to the single linkage specialised code.

Running the same code with average linkage for me gives:

>>> scipy.cluster.hierarchy.average(numpy.random.random_sample((2**16,10)))
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File
"/srv/sw/scipy/0.17.0/lib/python2.7/site-packages/scipy/cluster/hierarchy.py",
line 308, in average
    return linkage(y, method='average', metric='euclidean')
  File
"/srv/sw/scipy/0.17.0/lib/python2.7/site-packages/scipy/cluster/hierarchy.py",
line 660, in linkage
    int(_cpy_non_euclid_methods[method]))
  File "scipy/cluster/_hierarchy.pyx", line 733, in
scipy.cluster._hierarchy.linkage (scipy/cluster/_hierarchy.c:9058)
ValueError: negative dimensions are not allowed

Line 733 of _hierarchy.pyx is:

cdef double[:] D = np.ndarray(n * (n - 1) / 2, dtype=np.double)

I guess there is an overflow when computing the array size n * (n - 1) /
2 here.

In the single linkage code, I would guess there is an overflow when
computing condensed_index on line 874:

    M[j] = dists[condensed_index(n, i, j)]

since the condensed_index function returns an int.

On 28/04/2016 11:32 AM, argriffing wrote:

This segfaults for me too. Is it related to #2089
#2089?

—
You are receiving this because you authored the thread.
Reply to this email directly or view it on GitHub
#6107 (comment)

[rgommers]

I guess there is an overflow when computing the array size n * (n - 1) / 2 here.

That sounds right, addressed in gh-6638. I don't have enough memory right now to test this fix ((2**16, 10) takes 16 GB). Can someone check that the problem is gone?

In the single linkage code, I would guess there is an overflow when computing condensed_index on line 874:

That code disappeared with the rewrite of single linkage in gh-6495.

[mortonjt]

Hi @rgommers . We've tested this on 100k points using the fix that you suggested, and it doesn't look like the problem is going away.

When we first ran this without the fix, we got the following message.

ValueError: got differing extents in dimension 0 (got 75923957 and
4370891253)

Note that 4370891253 - 2**32 = 75923957, so this is likely due to integer overflow.

After compiling with the fix that you suggested, we get the following error

More recent error message:

/var/spool/torque/mom_priv/jobs/122263.barnacle.ucsd.edu.SC: line 11:
17554 Segmentation fault      (core dumped) python
/home/jlefler/balance_analysis/scripts/tree.py

We are running on 100GB of memory, on a CentOS machine.
Are there any other diagnostics that you'd recommend to further probe into this problem?
Let us know, we'd also love to see this happen!

cc @wasade

[rgommers]

Thanks @mortonjt, looks like there's another issue on top of the fixed overflow then. I'll get a machine with more RAM next Friday, so should be able to reproduce this myself. Probably easier than some back and forth debug instructions. But good that you can test, I'll ping you if I make more changes in the PR.

[mortonjt]

Hi @rgommers . @jlefty has followed up with this issue. There were a couple of modifications on top of your fix that could make this work for clustering ~93,000 elements.

1. It looks like there is another issue coming up from condensed_index If all of the int declarations (both parameters and return types) are replaced with np.int64_t, that will resolve this problem
2. i_start needs to be declared as np.int64_t.

Does this make sense? We can even push in a PR to resolve this issue.

[rgommers]

@mortonjt that does make sense. condensed_index multiplies two large ints, so will indeed overflow in the current implementation.

I'll get around to updating my PR in the weekend. If you want to send me a PR (to my cluster-overflow branch to get it into gh-6638), that would be very welcome of course.

[Josh-Lefler]

I just pushed a PR here. Not sure what is going on with the previous commits ...