[Numpy-discussion] Re: scipy.stats.itemfreq: overflow with add.reduce

[Tim C. <tc...@op...>]

Hans Georg Krauthaeuser wrote:
> Hans Georg Krauthaeuser schrieb:
> 
>>Hi All,
>>
>>I was playing with scipy.stats.itemfreq when I observed the following 
>>overflow:
>>
>>In [119]:for i in [254,255,256,257,258]:
>>   .....:    l=[0]*i
>>   .....:    print i, stats.itemfreq(l), l.count(0)
>>   .....:
>>254 [ [  0 254]] 254
>>255 [ [  0 255]] 255
>>256 [ [0 0]] 256
>>257 [ [0 1]] 257
>>258 [ [0 2]] 258
>>
>>itemfreq is pretty small (in stats.py):
>>
>>----------------------------------------------------------------------
>>def itemfreq(a):
>>    """
>>Returns a 2D array of item frequencies.  Column 1 contains item values,
>>column 2 contains their respective counts.  Assumes a 1D array is passed.
>>
>>Returns: a 2D frequency table (col [0:n-1]=scores, col n=frequencies)
>>"""
>>    scores = _support.unique(a)
>>    scores = sort(scores)
>>    freq = zeros(len(scores))
>>    for i in range(len(scores)):
>>        freq[i] = add.reduce(equal(a,scores[i]))
>>    return array(_support.abut(scores, freq))
>>----------------------------------------------------------------------
>>
>>It seems that add.reduce is the source for the overflow:
>>
>>In [116]:from scipy import *
>>
>>In [117]:for i in [254,255,256,257,258]:
>>   .....:    l=[0]*i
>>   .....:    print i, add.reduce(equal(l,0))
>>   .....:
>>254 254
>>255 255
>>256 0
>>257 1
>>258 2
>>
>>Is there any possibility to avoid the overflow?

Apropos the preceding, herewith a thread on the Numpy list from a more
than a few months ago. The take-home message is that for integer arrays,
add.reduce is very fast at producing results which fall into two
categories: a) correct or b) incorrect due to overflow. Unfortunately
there is no equally quick method of determining into which of these two
categories any specific result returned by add.reduce falls.

Tim C

From: Tim Churches <tc...@op...>
To: Todd Miller <jm...@st...>
Date: Apr 12 2005 - 7:51am

Todd Miller wrote:
> On Sun, 2005-04-10 at 10:23 +1000, Tim Churches wrote:
>
>>I just got caught by code equivalent to this (with NumPy 23.8 on 32 bit
>>Linux):
>>
>> >>> import Numeric as N
>> >>> a = N.array((2000000000,1000000000),typecode=N.Int32)
>> >>> N.add.reduce(a)
>>-1294967296
>>
>>OK, it is an elementary mistake, but the silent overflow caught me
>>unawares. casting the array to Float64 before summing it avoids the
>>error, but in my instance the actual data is a rank-1 array of 21
>>million integers with a mean value of about 140 (which adds up more than
>>sys.maxint), and casting to Float64 will use quite a lot of memory (as
>>well as taking some time).
>>
>>Any advice for catching or avoiding such overflow without always
>>incurring a performance and memory hit by always casting to Float64?
>
>
> Here's what numarray does:
>
>
>>>>import numarray as N
>>>>a = N.array((2000000000,1000000000),typecode=N.Int32)
>>>>N.add.reduce(a)
>
> -1294967296
>
> So basic reductions in numarray have the same "careful while you're
> shaving" behavior as Numeric; it's fast but easy to screw up.

Sure, but how does one be careful? It seems that for any array of two
integers or more which could sum to more than sys.maxint or less than
-sys.maxint, add.reduce() in both NumPy and Numeric will give either a)
the correct answer or b) the incorrect answer, and short of adding up
the array using a safer but much slower method, there is no way of
determining if the answer provided (quickly) by add.reduce is right or
wrong? Which seems to make it fast but useless (for integer arrays, at
least? Is that an unfair summary? Can anyone point me towards a method
for using add.reduce() on small arrays of large integers with values in
the billions, or on large arrays of fairly small integer values, which
will not suddenly and without warning give the wrong answer?

>
> But:
>
>
>>>>a.sum()
>
> 3000000000L
>
>>>>a.sum(type='d')
>
> 3000000000.0
>
> a.sum() blockwise upcasts to the largest type of kind on the fly, in
> this case, Int64. This avoids the storage overhead of typecasting the
> entire array.

That's on a 64-bit platform, right? The same method could be used to
cast the accumulator to a Float64 on a 32-bit platform to avoid casting
the entire array?

> A better name for the method would have been sumall() since it sums all
> elements of a multi-dimensional array. The flattening process reduces
> on one dimension before flattening preventing a full copy of a
> discontiguous array. It could be smarter about choosing the dimension
> of the initial reduction.

OK, thanks. Unfortunately it is not possible for us to port our
application to numarray at the moment. But the insight is most helpful.

Tim C