[Numpy-discussion] bug in oldnumeric.ma
The following code works with numpy.ma but not numpy.oldnumeric.ma, obviously it shouldn't have missing values... Note that replacing 2. with 2 (int) works, dtype does not seem to matter import numpy.oldnumeric.ma as MA,numpy s = MA.array([ 12.16271591, 11.19478798, 10.27440453, 9.60334778, 9.2451086, 9.13312435, 9.11890984, 9.03033924, 8.7346344, 8.18558788, 7.43921041, 6.62947559, 5.90856123, 5.37711906, 5.03489971, 4.77586126, 4.4189887 , 3.76522899, 2.65106893, 0.99725384, -1.1513288, -3.58740163, -5.93836021, -7.72085428, -8.4840517, -7.98519516, -6.31602335, -3.89117384, -1.29867446, 0.90583926, 2.3683548, 3.00160909, 2.94496846, 2.44637036, 1.7288276, 0.89646715, -0.07413431, -1.25960362, -2.67903948, -4.22316313, -5.66368866, -6.74747849, -7.31541586, -7.37950087, -7.12536144, -6.82764053, -6.74864483, -7.04398584, -7.73052883, -8.70125961, -9.77933311, -10.78447914, -11.59126091, -12.1706562, -12.60001087, -13.03606987, -13.64736843 , -14.51338005, -15.52330303, -16.33303833, -16.46606064, -15.5264, -13.45486546, -10.63810158, -7.79133797, -5.66280842, -4.72991323, -5.03621674, -6.21381569, -7.6610961, -8.78764057, -9.22090816, -8.91826916, -8.14496899, -7.33739614, -6.91325617, -7.09938431, -7.84141493, -8.829319, -9.63222504, -9.88238335, -9.44395733, -8.4723177, -7.35424185, -6.54324722, -6.36360407, -6.8554, -7.74132967, -8.51024342, -8.62310696, -7.73065567, -5.82352972, -3.24726033, -0.5881027 , 1.51040995, 2.52599192, 2.19137073, 0.56473863, -1.97977543, -4.84694195, -7.40427446, -9.2083149, -10.16563416, -10.53237438, -10.75102997, -11.21533489, -12.09201908, -13.28770447, -14.54858589, -15.62684536, -16.40594101, -16.91949844, -17.27025032, -17.51157188, -17.57112694, -17.25679016, -16.32791138, -14.58930111, -11.96431351, -8.52671432, -4.49569035, -0.20288418, 3.96036053, 7.60605049, 10.41992378, 12.21882153, 12.98544502, 12.87247849],dtype='f') s2c=MA.power(s,2.) print MA.count(s) print MA.count(s2c) ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
Charles Doutriaux wrote: > The following code works with numpy.ma but not numpy.oldnumeric.ma, No, this is a bug in numpy.ma also; power is broken: In [1]:import numpy as np In [2]:x = np.ma.array([-1.1]) In [3]:x**2.0 ### This works Out[3]: masked_array(data = [1.21], mask = [False], fill_value=1e+20) In [4]:np.ma.power(x, 2.0) ### This doesn't Out[4]: masked_array(data = [--], mask = [ True], fill_value=1e+20) Here is the code in ma/core.py, which is masking the output for negative inputs: def power(a, b, third=None): """Computes a**b elementwise. Masked values are set to 1. """ if third is not None: raise MAError, "3-argument power not supported." ma = getmask(a) mb = getmask(b) m = mask_or(ma, mb) fa = getdata(a) fb = getdata(b) if fb.dtype.char in typecodes["Integer"]: return masked_array(umath.power(fa, fb), m) md = make_mask((fa < 0), shrink=True) wrong m = mask_or(m, md) if m is nomask: return masked_array(umath.power(fa, fb)) else: fa = fa.copy() fa[(fa < 0)] = 1 wrong return masked_array(umath.power(fa, fb), m) Eric ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/7 Eric Firing <[EMAIL PROTECTED]>: > Charles Doutriaux wrote: > > The following code works with numpy.ma but not numpy.oldnumeric.ma, > > No, this is a bug in numpy.ma also; power is broken: While it's tempting to just call power() and mask out any NaNs that result, that's going to be a problem if people have their environments set to raise exceptions on the production of NaNs. Is it an adequate criterion to check (a<0) & (round(b)==b)? We have to be careful: In [16]: np.array([-1.0])**(2.0**128) Warning: invalid value encountered in power Out[16]: array([ nan]) 2.0**128 cannot be distinguished from nearby non-integral values, so this is reasonable behaviour (and a weird corner case), but In [23]: np.round(2.0**128) == 2.0**128 Out[23]: True Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
All, Yes, there is a problem with ma.power: masking negative data should be restricted to the case of an exponent between -1. and 1. only, don't you think ? On Wednesday 07 May 2008 18:47:18 Anne Archibald wrote: > 2008/5/7 Eric Firing <[EMAIL PROTECTED]>: > > Charles Doutriaux wrote: > > > The following code works with numpy.ma but not numpy.oldnumeric.ma, > > > > No, this is a bug in numpy.ma also; power is broken: > > While it's tempting to just call power() and mask out any NaNs that > result, that's going to be a problem if people have their environments > set to raise exceptions on the production of NaNs. Is it an adequate > criterion to check (a<0) & (round(b)==b)? We have to be careful: > > In [16]: np.array([-1.0])**(2.0**128) > Warning: invalid value encountered in power > Out[16]: array([ nan]) > > 2.0**128 cannot be distinguished from nearby non-integral values, so > this is reasonable behaviour (and a weird corner case), but > > In [23]: np.round(2.0**128) == 2.0**128 > Out[23]: True > > Anne > ___ > Numpy-discussion mailing list > Numpy-discussion@scipy.org > http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/7 Pierre GM <[EMAIL PROTECTED]>: > All, > Yes, there is a problem with ma.power: masking negative data should be > restricted to the case of an exponent between -1. and 1. only, don't you > think ? No, there's a problem with any fractional exponent (with even denominator): x**(3/2) == (x**3)**(1/2). And of course in a floating-point world, you can't really ask whether the denominator is even or not. So any non-integer power is trouble. The draconian approach would be to simply disallow negative numbers to be raised to exponents of type float, but that's going to annoy a lot of people who have integers which happen to be represented in floats. (Representing integers in floats is exact up to some fairly large value.) So the first question is "how do we recognize floats with integer values?". Unfortunately that's not the real problem. The real problem is "how do we predict when power() is going to produce a NaN?" Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
On Wednesday 07 May 2008 20:38:22 Anne Archibald wrote: > 2008/5/7 Pierre GM <[EMAIL PROTECTED]>: > > All, > > Yes, there is a problem with ma.power: masking negative data should be > > restricted to the case of an exponent between -1. and 1. only, don't you > > think ? > > No, there's a problem with any fractional exponent (with even > denominator): x**(3/2) == (x**3)**(1/2). Argh. Good point. > The real > problem is "how do we predict when power() is going to produce a NaN?" An alternative would be to forget about it: let power() output NaNs, and fix them afterwards with fix_invalid. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
Anne Archibald wrote: > 2008/5/7 Pierre GM <[EMAIL PROTECTED]>: > >> All, >> Yes, there is a problem with ma.power: masking negative data should be >> restricted to the case of an exponent between -1. and 1. only, don't you >> think ? >> > > No, there's a problem with any fractional exponent (with even > denominator): x**(3/2) == (x**3)**(1/2). And of course in a > floating-point world, you can't really ask whether the denominator is > even or not. So any non-integer power is trouble. > > The draconian approach would be to simply disallow negative numbers to > be raised to exponents of type float, but that's going to annoy a lot > of people who have integers which happen to be represented in floats. > (Representing integers in floats is exact up to some fairly large > value.) So the first question is "how do we recognize floats with > integer values?". Unfortunately that's not the real problem. The real > problem is "how do we predict when power() is going to produce a NaN?" > > Anne > Is there a rule against squaring away the negatives? def not_your_normal_pow( x, y ): return exp( log( power( x, 2) ) * y / 2 ) Which still needs some work for x==0. Jon ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/7 Jonathan Wright <[EMAIL PROTECTED]>: > Is there a rule against squaring away the negatives? > > def not_your_normal_pow( x, y ): return exp( log( power( x, 2) ) * y / 2 ) > > Which still needs some work for x==0. Well, it means (-1.)**(3.) becomes 1., which is probably not what the user expected... Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/7 Pierre GM <[EMAIL PROTECTED]>: > On Wednesday 07 May 2008 20:38:22 Anne Archibald wrote: > > 2008/5/7 Pierre GM <[EMAIL PROTECTED]>: > > > All, > > > Yes, there is a problem with ma.power: masking negative data should be > > > restricted to the case of an exponent between -1. and 1. only, don't you > > > think ? > > > > No, there's a problem with any fractional exponent (with even > > denominator): x**(3/2) == (x**3)**(1/2). > > Argh. Good point. > > > > The real > > problem is "how do we predict when power() is going to produce a NaN?" > An alternative would be to forget about it: let power() output NaNs, and fix > them afterwards with fix_invalid. Tempting, but the user may have used seterr() to arrange that exceptions are raised when this happens, which is going to put a spanner in the works. (And temporarily changing seterr() is problematic in a multithreaded context...) Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
On Wed, May 7, 2008 at 12:12 PM, Pierre GM <[EMAIL PROTECTED]> wrote: > Yes, there is a problem with ma.power: masking negative data should be > restricted to the case of an exponent between -1. and 1. only, don't you > think ? Charles Doutriaux has suggested that 1.1.0 shouldn't be released until this bug is fixed. I am inclined to somewhat agree with this, but would like to know if others agree with this sentiment. Please respond with a +1 or -1 to making this a 1.1.0 blocker. I also not sure how closely this issue is related with this open ticket: http://projects.scipy.org/scipy/numpy/ticket/301 I would appreciate it if someone could clarify this issue. Thanks, -- Jarrod Millman Computational Infrastructure for Research Labs 10 Giannini Hall, UC Berkeley phone: 510.643.4014 http://cirl.berkeley.edu/ ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
Jarrod Millman wrote: > On Wed, May 7, 2008 at 12:12 PM, Pierre GM <[EMAIL PROTECTED]> wrote: >> Yes, there is a problem with ma.power: masking negative data should be >> restricted to the case of an exponent between -1. and 1. only, don't you >> think ? > > Charles Doutriaux has suggested that 1.1.0 shouldn't be released until > this bug is fixed. I am inclined to somewhat agree with this, but > would like to know if others agree with this sentiment. Please > respond with a +1 or -1 to making this a 1.1.0 blocker. > > I also not sure how closely this issue is related with this open > ticket: http://projects.scipy.org/scipy/numpy/ticket/301 > I would appreciate it if someone could clarify this issue. It is a completely different issue; #301 applies to numpy.power with integer arguments including a negative power; Charles Doutriaux pointed out a bug in numpy.ma.power (and in the oldnumeric.ma) that masked the result if the power was negative. Pierre GM fixed this in r5137, so that the result is masked only if the first argument is negative and the power is between -1 and 1. Eric ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/9 Eric Firing <[EMAIL PROTECTED]>: > Pierre GM fixed this in r5137, so that the result is masked only if the > first argument is negative and the power is between -1 and 1. Why should the output be masked only when the power is between -1 and 1? Other powers also produce nan's. Either way, the docstring should be updated: Computes a**b elementwise. Masked values are set to 1. Regards Stéfan ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
Stefan, (and Jarrod and Pierre) (Context for anyone new to the thread: the subject is slightly misleading, because the bug is/was present in both oldnumeric.ma and numpy.ma; the discussion of fix pertains to the latter only.) Regarding your objections to r5137: good point. I wondered about that. I think the function should look like this (although it might be possible to speed up the implementation for the most common case): def power(a, b, third=None): """Computes a**b elementwise. Where a is negative and b has a non-integer value, the result is masked. """ if third is not None: raise MAError, "3-argument power not supported." ma = getmask(a) mb = getmask(b) m = mask_or(ma, mb) fa = getdata(a) fb = getdata(b) if fb.dtype.char in typecodes["Integer"]: return masked_array(umath.power(fa, fb), m) md = make_mask((fb != fb.astype(int)) & (fa < 0), shrink=True) m = mask_or(m, md) if m is nomask: return masked_array(umath.power(fa, fb)) else: fa = fa.copy() fa[m] = 1 return masked_array(umath.power(fa, fb), m) I don't have time right now to turn this into a proper patch, complete with test case, but if no one else can do it sooner then I can probably do it in the next day or two. Here is a quick partial example of the behavior of the version above: In [1]:import numpy as np In [2]:xx = np.ma.array([-2.2, -2.2, 2.2, 2.2], mask = [True, False, False, True]) In [3]:xx Out[3]: masked_array(data = [-- -2.2 2.2 --], mask = [ True False False True], fill_value=1e+20) In [4]:np.ma.power(xx, 2.0) Out[4]: masked_array(data = [-- 4.84 4.84 --], mask = [ True False False True], fill_value=1e+20) In [5]:np.ma.power(xx, -2.0) Out[5]: masked_array(data = [-- 0.206611570248 0.206611570248 --], mask = [ True False False True], fill_value=1e+20) In [6]:np.ma.power(xx, -2.1) Out[6]: masked_array(data = [-- -- 0.190946793699 --], mask = [ True True False True], fill_value=1e+20) Eric Stéfan van der Walt wrote: > 2008/5/9 Eric Firing <[EMAIL PROTECTED]>: >> Pierre GM fixed this in r5137, so that the result is masked only if the >> first argument is negative and the power is between -1 and 1. > > Why should the output be masked only when the power is between -1 and > 1? Other powers also produce nan's. Either way, the docstring should > be updated: > > Computes a**b elementwise. > > Masked values are set to 1. > > Regards > Stéfan > ___ > Numpy-discussion mailing list > Numpy-discussion@scipy.org > http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
On Friday 09 May 2008 12:55:39 Eric Firing wrote: > Stefan, (and Jarrod and Pierre) > Regarding your objections to r5137: good point. I wondered about that. > I think the function should look like this (although it might be > possible to speed up the implementation for the most common case): OK, on it. ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/9 Eric Firing <[EMAIL PROTECTED]>: > Stefan, (and Jarrod and Pierre) > > (Context for anyone new to the thread: the subject is slightly > misleading, because the bug is/was present in both oldnumeric.ma and > numpy.ma; the discussion of fix pertains to the latter only.) > > Regarding your objections to r5137: good point. I wondered about that. > I think the function should look like this (although it might be > possible to speed up the implementation for the most common case): [...] > md = make_mask((fb != fb.astype(int)) & (fa < 0), shrink=True) Unfortunately this isn't quite the right condition: In [18]: x = 2.**35; numpy.array([-1.])**x; numpy.array(x).astype(int)==x Out[18]: array([ 1.]) Out[18]: False Switching to int64 seems to help: In [27]: x = 2.**62; numpy.array([-1.])**x; numpy.array(x).astype(numpy.int64)==x Out[27]: array([ 1.]) Out[27]: True This seems to work, but may be platform-dependent: 2**62+1 cannot be represented as an IEEE float, so whether pow() successfully deals with it may be different for machines that don't work with 80-bit floating-point internally. A suspenders-and-belt approach would check for NaNs and add them to the mask, but this still doesn't cover the case where the user has numpy set to raise exceptions any time NaNs are generated. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
Anne Archibald wrote: > 2008/5/9 Eric Firing <[EMAIL PROTECTED]>: >> Stefan, (and Jarrod and Pierre) >> >> (Context for anyone new to the thread: the subject is slightly >> misleading, because the bug is/was present in both oldnumeric.ma and >> numpy.ma; the discussion of fix pertains to the latter only.) >> >> Regarding your objections to r5137: good point. I wondered about that. >> I think the function should look like this (although it might be >> possible to speed up the implementation for the most common case): > [...] >> md = make_mask((fb != fb.astype(int)) & (fa < 0), shrink=True) > > Unfortunately this isn't quite the right condition: > > In [18]: x = 2.**35; numpy.array([-1.])**x; numpy.array(x).astype(int)==x > Out[18]: array([ 1.]) > Out[18]: False > > Switching to int64 seems to help: > > In [27]: x = 2.**62; numpy.array([-1.])**x; > numpy.array(x).astype(numpy.int64)==x > Out[27]: array([ 1.]) > Out[27]: True > > This seems to work, but may be platform-dependent: 2**62+1 cannot be > represented as an IEEE float, so whether pow() successfully deals with > it may be different for machines that don't work with 80-bit > floating-point internally. > > A suspenders-and-belt approach would check for NaNs and add them to > the mask, but this still doesn't cover the case where the user has > numpy set to raise exceptions any time NaNs are generated. There may not be a perfect solution, but I suspect your suggestion to use int64 is more than good enough to get things going for a 1.1 release. The docstring could note the limitation. If it is established that the calculation will fail for a power outside some domain, then such a domain check could be added to the mask. Eric > > Anne > ___ > Numpy-discussion mailing list > Numpy-discussion@scipy.org > http://projects.scipy.org/mailman/listinfo/numpy-discussion ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
All, I fixed the power function in numpy.ma following Anne's suggestion: compute first, mask the problems afterwards. It's a quick and dirty fix that crashes if the user has set its error system to raise an exception on invalid (np.seterr(invalid='raise')), but it works otherwise and keeps subclasses (such as TimeSeries). I will have to modify the .__pow__ method so that ma.power is called: right now, a**b calls ndarray(a).__pow__(b), which may yield NaNs and Infs. What should I do with oldnumeric.ma.power ? Try to fix it the same way, or leave the bug ? I'm not that enthusiastic to have to debug the old package, but if it's part of the job... ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
2008/5/12 Pierre GM <[EMAIL PROTECTED]>: > I fixed the power function in numpy.ma following Anne's suggestion: compute > first, mask the problems afterwards. It's a quick and dirty fix that crashes > if the user has set its error system to raise an exception on invalid > (np.seterr(invalid='raise')), but it works otherwise and keeps subclasses > (such as TimeSeries). > I will have to modify the .__pow__ method so that ma.power is called: right > now, a**b calls ndarray(a).__pow__(b), which may yield NaNs and Infs. > What should I do with oldnumeric.ma.power ? Try to fix it the same way, or > leave the bug ? I'm not that enthusiastic to have to debug the old package, > but if it's part of the job... We should leave the oldnumeric.ma package alone. Even if its `power` is broken, some packages may depend on it. We'll provide it in 1.2 for backward compatibility, and get rid of it after 1.3. Regards Stéfan ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] bug in oldnumeric.ma
On Monday 12 May 2008 12:04:24 Stéfan van der Walt wrote: > 2008/5/12 Pierre GM <[EMAIL PROTECTED]>: > > What should I do with oldnumeric.ma.power ? Try to fix it the same > > way, or leave the bug ? I'm not that enthusiastic to have to debug the > > old package, but if it's part of the job... > > We should leave the oldnumeric.ma package alone. Even if its `power` > is broken, some packages may depend on it. We'll provide it in 1.2 > for backward compatibility, and get rid of it after 1.3. OK then, I prefer that... Additional question: when raising to a power in place, NaNs/Infs can show up in the _data part: should I set those invalid data to fill_value or not ? ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
