On Mon, Jan 23, 2012 at 1:33 PM, Travis Oliphant <teoliphant@gmail.com>wrote:

Can you determine where the problem is, precisely. In other words, can you verify that c is not getting filled in correctly?

You are no doubt going to get overflow in the summation as you have a uint8 parameter. But, having that overflow be exactly '0' would be surprising.

Can you verify that a and b are getting created correctly? Also, 'c' should be a 2-d array, can you verify that? Can you take the sum along the -1 axis and the 0 axis separately:

print a.shape print b.shape print c.shape

c[1000000:].sum(axis=0) d = c[1000000:].sum(axis=-1) print d[:100] print d[-100:]

I am getting the same results as David. It looks like c just "stopped filling in" partway through the array. I don't think there is any overflow issue, since the result of sum() is up-promoted to uint64 when I do that. Travis, here are the outputs at my end - I cut out many zeros for brevity: In [7]: print a.shape (5000000, 972) In [8]: print b.shape (4993210,) In [9]: print c.shape (4993210, 972) In [10]: c[1000000:].sum(axis=0) Out[10]: array([0, 0, 0, .... , 0]) In [11]: d = c[1000000:].sum(axis=-1) In [12]: print d[:100] [0 0 0 ... 0 0] In [13]: print d[-100:] [0 0 0 ... 0 0 0] I looked at sparse subsamples with matplotlib - specifically, imshow(a[::1000, :]) - and the a array looks correct (random values everywhere), but c is zero past a certain row number. In fact, it looks like it becomes zero at row 575419 - I think for all rows in c beyond row 574519, the values will be zero. For lower row numbers, I think they are correctly filled (at least, by the sparse view in matplotlib). In [15]: a[b[574519], 350:360] Out[15]: array([143, 155, 11, 30, 212, 149, 110, 164, 165, 120], dtype=uint8) In [16]: c[574519, 350:360] Out[16]: array([143, 155, 11, 30, 212, 149, 0, 0, 0, 0], dtype=uint8) I'm using EPD 7.1, numpy 1.6.1, Linux installation (I don't know the kernel details) HTH, Aronne

On Mon, Jan 23, 2012 at 7:55 PM, David Warde-Farley <wardefar@iro.umontreal.ca> wrote:

I've reproduced this (rather serious) bug myself and confirmed that it exists in master, and as far back as 1.4.1.

I'd really appreciate if someone could reproduce and confirm on another machine, as so far all my testing has been on our single high-memory machine.

I see the same behaviour on a Winodows machine with numpy 1.6.1. But I don't think it is an indexing problem - rather something with the random number creation. a itself is already zeros for high indexes.  In [8]: b[1000000:1000010] Out[8]: array([3429029, 1251819, 4292918, 2249483, 757620, 3977130, 3455449, 2005054, 2565207, 3114930]) In [9]: a[b[1000000:1000010]] Out[9]: array([[0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], ..., [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0]], dtype=uint8) In [41]: a[581350:,0].sum() Out[41]: 0 Cheers Robin

Thanks, David

On Mon, Jan 23, 2012 at 05:23:28AM -0500, David Warde-Farley wrote:

...

A colleague has run into this weird behaviour with NumPy 1.6.1, EPD 7.1-2, on Linux (Fedora Core 14) 64-bit:

...

a = numpy.array(numpy.random.randint(256,size=(5000000,972)),dtype='uint8') b = numpy.random.randint(5000000,size=(4993210,)) c = a[b]

It seems c is not getting filled in full, namely:

...

In [14]: c[1000000:].sum() Out[14]: 0

I haven't been able to reproduce this quite yet, I'll try to find a machine with sufficient memory tomorrow. But does anyone have any insight in the mean time? It smells like some kind of integer overflow bug.

Thanks,

David _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

---

NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion

On 1/23/2012 12:33 PM, David Warde-Farley wrote:

On Mon, Jan 23, 2012 at 08:38:44PM +0100, Robin wrote:

...

On Mon, Jan 23, 2012 at 7:55 PM, David Warde-Farley <wardefar@iro.umontreal.ca> wrote:

...

I've reproduced this (rather serious) bug myself and confirmed that it exists in master, and as far back as 1.4.1.

I'd really appreciate if someone could reproduce and confirm on another machine, as so far all my testing has been on our single high-memory machine.

I see the same behaviour on a Winodows machine with numpy 1.6.1. But I don't think it is an indexing problem - rather something with the random number creation. a itself is already zeros for high indexes.  In [8]: b[1000000:1000010] Out[8]: array([3429029, 1251819, 4292918, 2249483, 757620, 3977130, 3455449, 2005054, 2565207, 3114930])

In [9]: a[b[1000000:1000010]] Out[9]: array([[0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], ..., [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0], [0, 0, 0, ..., 0, 0, 0]], dtype=uint8)

In [41]: a[581350:,0].sum() Out[41]: 0

Hmm, this seems like a separate bug to mine. In mine, 'a' is indeed being filled in -- the problem arises with c alone.

So, another Windows-specific bug to add to the pile, perhaps? :(

David

Maybe this explains the win-amd64 behavior: There are a couple of places in mtrand where array indices and sizes are C long instead of npy_intp, for example in the randint function: <https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx#L8...> Christoph

Den 23.01.2012 22:08, skrev Christoph Gohlke:

Maybe this explains the win-amd64 behavior: There are a couple of places in mtrand where array indices and sizes are C long instead of npy_intp, for example in the randint function:

<https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx#L8...>

Both i and length could overflow here. It should overflow on allocation of more than 2 GB. There is also a lot of C longs in the internal state (line 55-105), as well as the other functions. Producing 2 GB of random ints twice fails:

...

...

import numpy as np np.random.randint(5000000,size=(2*1024**3,)) array([0, 0, 0, ..., 0, 0, 0]) np.random.randint(5000000,size=(2*1024**3,))

Traceback (most recent call last): File "<pyshell#3>", line 1, in <module> np.random.randint(5000000,size=(2*1024**3,)) File "mtrand.pyx", line 881, in mtrand.RandomState.randint (numpy\random\mtrand\mtrand.c:6040) MemoryError

...

...

Sturla

On 24.01.2012 06:32, Sturla Molden wrote:

The use of C long affects all the C and Pyrex source code in mtrand module, not just mtrand.pyx. All of it is fubar on Win64.

randomkit.c handles C long correctly, I think. There are different codes for 32 and 64 bit C long, and buffer sizes are size_t. Sturla

On Tue, Jan 24, 2012 at 08:37, Sturla Molden <sturla@molden.no> wrote:

On 24.01.2012 09:21, Sturla Molden wrote:

...

randomkit.c handles C long correctly, I think. There are different codes for 32 and 64 bit C long, and buffer sizes are size_t.

distributions.c take C longs as parameters e.g. for the binomial distribution. mtrand.pyx correctly handles this, but it can give an unexpected overflow error on 64-bit Windows:

In [1]: np.random.binomial(2**31, .5) --------------------------------------------------------------------------- OverflowError                             Traceback (most recent call last) C:\Windows\system32\<ipython-input-1-000aa0626c42> in <module>() ----> 1 np.random.binomial(2**31, .5)

C:\Python27\lib\site-packages\numpy\random\mtrand.pyd in mtrand.RandomState.binomial (numpy\random\mtrand\mtrand.c:13770)()

OverflowError: Python int too large to convert to C long

On systems where C longs are 64 bit, this is likely not to produce an error.

This begs the question if also randomkit.c and districutions.c should be changed to use npy_intp for consistency across all platforms.

There are two different uses of long that you need to distinguish. One is for sizes, and one is for parameters and values. The sizes should definitely be upgraded to npy_intp. The latter shouldn't; these should remain as the default integer type of Python and numpy, a C long. The reason longs are used for sizes is that I wrote mtrand for Numeric and Python 2.4 before numpy was even announced (and I don't think we had npy_intp at the time I merged it into numpy, but I could be wrong). Using longs for sizes was the order of the day. I don't think I had even touched a 64-bit machine that wasn't a DEC Alpha at the time, so I knew very little about the issues. So yes, please, fix whatever you can. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth."   -- Umberto Eco

On Tue, Jan 24, 2012 at 09:15:01AM +0000, Robert Kern wrote:

On Tue, Jan 24, 2012 at 08:37, Sturla Molden <sturla@molden.no> wrote:

...

On 24.01.2012 09:21, Sturla Molden wrote:

...

randomkit.c handles C long correctly, I think. There are different codes for 32 and 64 bit C long, and buffer sizes are size_t.

distributions.c take C longs as parameters e.g. for the binomial distribution. mtrand.pyx correctly handles this, but it can give an unexpected overflow error on 64-bit Windows:

In [1]: np.random.binomial(2**31, .5) --------------------------------------------------------------------------- OverflowError                             Traceback (most recent call last) C:\Windows\system32\<ipython-input-1-000aa0626c42> in <module>() ----> 1 np.random.binomial(2**31, .5)

C:\Python27\lib\site-packages\numpy\random\mtrand.pyd in mtrand.RandomState.binomial (numpy\random\mtrand\mtrand.c:13770)()

OverflowError: Python int too large to convert to C long

On systems where C longs are 64 bit, this is likely not to produce an error.

This begs the question if also randomkit.c and districutions.c should be changed to use npy_intp for consistency across all platforms.

There are two different uses of long that you need to distinguish. One is for sizes, and one is for parameters and values. The sizes should definitely be upgraded to npy_intp. The latter shouldn't; these should remain as the default integer type of Python and numpy, a C long.

Hmm. Seeing as the width of a C long is inconsistent, does this imply that the random number generator will produce different results on different platforms? Or do the state dynamics prevent it from ever growing in magnitude to the point where that's an issue? David

On Thu, Jan 26, 2012 at 17:39, Sturla Molden <sturla@molden.no> wrote:

Den 24.01.2012 17:19, skrev David Warde-Farley:

...

Hmm. Seeing as the width of a C long is inconsistent, does this imply that the random number generator will produce different results on different platforms?

If it does, it is a C programming mistake. C code should never depend on the exact size of a long, only it's minimum size.  ISO C defines other datatypes if an exact integer size is needed (include stdint.h), but ANSI C used for NumPy does not.

I think you're subtly misunderstanding his question. He's not asking if the code is written such that it semantically requires a long to have one specific size or another (and indeed, it is not). However, it is true that the code may behave differently for the same inputs on different machines with different long sizes. Namely, some part of the computation may overflow on 32-bit longs while giving an accurate answer with 64-bit longs. They just have different domains of accuracy over their inputs. It is not necessarily a mistake to take advantage of the extra room when it is available. That is the reason that Python ints are C longs and why numpy follows suit. But unfortunately, it is true that at least some of the distributions do have different behavior when using 64-bit longs than when using 32-bit longs. Here is an example of drawing from a binomial distribution with a large N on a 32-bit process and comparing it with results from a 64-bit process: [~]$ $PYBIN/python Enthought Python Distribution -- www.enthought.com Version: 7.1-2 (32-bit) Python 2.7.2 |EPD 7.1-2 (32-bit)| (default, Jul 27 2011, 13:29:32) [GCC 4.0.1 (Apple Inc. build 5493)] on darwin Type "packages", "demo" or "enthought" for more information.

...

...

import numpy as np prng = np.random.RandomState(1234567890) x32 = prng.binomial(500000000, 0.5, size=10000) x64 = np.load('x64.npy') np.save('x32.npy', x32) bad = (x64 != x32) bad.sum() 9449 bad[-9449:] array([ True, True, True, ..., True, True, True], dtype=bool) bad[-9449:].sum() 9449

binomial() is using a rejection algorithm for this set of inputs. For each draw, it is going to generate random numbers until they meet a certain condition. Both the 32-bit process and the 64-bit process draw the same exact numbers until the 552nd draw. Then, I suspect, there is an integer overflow in the 32-bit process causing the rejection algorithm to terminate either earlier or later than it otherwise should. Since the two processes have consumed different amounts of random numbers, the underlying uniform PRNG is no longer in the same state, so all of the numbers thereafter will be different. It's not clear to me how problematic this is. I haven't seen any difference when using reasonable input values (N=500000000 is a ridiculous number to be using with a binomial distribution). If I'm right that there is an overflow when using the 32-bit longs, then the results should not be trusted anyways, so there is no point in comparing them to the 64-bit results. It's just that the domain of validity with a 32-bit long is a bit smaller than when using a 64-bit long. The deviation of x32[551] from the mean is larger than the maximum deviation from the 64-bit results, so it is reasonably likely that the draw is just bogus.

...

...

np.max(abs(x64 - 250000000)) 44519 x32[551] - 250000000 47368

Often, the acceptance criterion is something of the form (X < something) while expecting X to be positive. An integer overflow would introduce a negative value somewhere in the computation and could easily "pass" this acceptance criterion when it really shouldn't have if the intermediate computations had been done without overflow. If anyone wants to debug this more thoroughly, this bit of code will get the PRNG into exactly the right state to see the difference on the next binomial() draw:

...

...

import numpy as np prng = np.random.RandomState(1234567890) blah = prng.binomial(500000000, 0.5, size=551)

If you run python under gdb, you can then set a breakpoint in rk_binomial_btpe() in distributions.c to step through the next call to prng.binomial(). Sometimes you can fix these issues in a rejection algorithm by checking for overflow and rejecting those cases. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth."   -- Umberto Eco

On Tue, Jan 24, 2012 at 08:47, Sturla Molden <sturla@molden.no> wrote:

The coding is also inconsistent, compare for example:

https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx#L1...

https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx#L2...

I'm sorry, what are you demonstrating there? -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth."   -- Umberto Eco

On Tue, Jan 24, 2012 at 09:19, Sturla Molden <sturla@molden.no> wrote:

On 24.01.2012 10:16, Robert Kern wrote:

...

I'm sorry, what are you demonstrating there?

Both npy_intp and C long are used for sizes and indexing.

Ah, yes. I think Travis added the multiiter code to cont1_array(), which does broadcasting, so he used npy_intp as is proper (and necessary to pass into the multiiter API). The other functions don't do broadcasting, so he didn't touch them. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth."   -- Umberto Eco

On Tue, Jan 24, 2012 at 6:24 PM, David Warde-Farley <wardefar@iro.umontreal.ca> wrote:

On Tue, Jan 24, 2012 at 06:00:05AM +0100, Sturla Molden wrote:

...

Den 23.01.2012 22:08, skrev Christoph Gohlke:

...

Maybe this explains the win-amd64 behavior: There are a couple of places in mtrand where array indices and sizes are C long instead of npy_intp, for example in the randint function:

<https://github.com/numpy/numpy/blob/master/numpy/random/mtrand/mtrand.pyx#L8...>

Both i and length could overflow here. It should overflow on allocation of more than 2 GB.

There is also a lot of C longs in the internal state (line 55-105), as well as the other functions.

Producing 2 GB of random ints twice fails:

Sturla, since you seem to have access to Win64 machines, do you suppose you could try this code:

...

...

...

a = numpy.ones((1, 972)) b = numpy.zeros((4993210,), dtype=int) c = a[b]

and verify that there's a whole lot of 0s in the matrix, specifically,

...

...

...

c[574519:].sum() 356.0 c[574520:].sum() 0.0

is the case on Linux 64-bit; is it the case on Windows 64?

Yes - I get exactly the same numbers in 64 bit windows with 1.6.1. Cheers Robin

On Tue, Jan 24, 2012 at 01:02:44PM -0500, David Warde-Farley wrote:

On Tue, Jan 24, 2012 at 06:37:12PM +0100, Robin wrote:

...

Yes - I get exactly the same numbers in 64 bit windows with 1.6.1.

Alright, so that rules out platform specific effects.

I'll try and hunt the bug down when I have some time, if someone more familiar with the indexing code doesn't beat me to it.

I've figured it out. In numpy/core/src/multiarray/mapping.c, PyArray_GetMap is using an int for a counter variable where it should be using an npy_intp. I've filed a pull request at https://github.com/numpy/numpy/pull/188 with a regression test. David

On Tue, Jan 24, 2012 at 3:30 PM, David Warde-Farley < wardefar@iro.umontreal.ca> wrote:

On Tue, Jan 24, 2012 at 01:02:44PM -0500, David Warde-Farley wrote:

...

On Tue, Jan 24, 2012 at 06:37:12PM +0100, Robin wrote:

...

Yes - I get exactly the same numbers in 64 bit windows with 1.6.1.

Alright, so that rules out platform specific effects.

I'll try and hunt the bug down when I have some time, if someone more familiar with the indexing code doesn't beat me to it.

I've figured it out. In numpy/core/src/multiarray/mapping.c, PyArray_GetMap is using an int for a counter variable where it should be using an npy_intp.

I've filed a pull request at https://github.com/numpy/numpy/pull/188 with a regression test.

I think this bug, or one like it, was reported a couple of years ago. But I don't recall if there was ever a ticket opened. Chuck