[Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)
On Friday 09 May 2008 17:13:02 Eric Firing wrote: > Anne Archibald wrote: > > 2008/5/9 Eric Firing <[EMAIL PROTECTED]>: > >> 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: > 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. Interestingly, I notice that MaskedArray misses a __pow__ method: it uses the ndarray.__pow__ method instead, that may outputs NaNs. In other terms, I gonna have to code MaskedArray.__pow__, following Eric' example, with int64 instead of int. But, why don't we compare: abs(np.array(b).astype(int)-b).max()http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)
Pierre GM wrote: > On Friday 09 May 2008 17:13:02 Eric Firing wrote: >> Anne Archibald wrote: >>> 2008/5/9 Eric Firing <[EMAIL PROTECTED]>: 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: > >> 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. > > Interestingly, I notice that MaskedArray misses a __pow__ method: it uses the > ndarray.__pow__ method instead, that may outputs NaNs. In other terms, I > gonna have to code MaskedArray.__pow__, following Eric' example, with int64 > instead of int. > But, why don't we compare: > abs(np.array(b).astype(int)-b).max() instead ? At which point will b be considered sufficiently close to an > integer > that x**b won't return NaN ? I don't think the .max() part of that is right; the test needs to be element-wise, and turned into a mask. It is also not clear to me that the test would actually catch all the cases where x**b would return NaN. It seems like some strategic re-thinking may be needed in the long run, if not immediately. There is a wide range of combinations of arguments that will trigger invalid results, whether Inf or NaN. The only way to trap and mask all of these is to use masked_invalid after the calculation, and this only works if the user has not disabled nan output. I have not checked recently, but based on earlier strategy discussions, I suspect that numpy.ma is already strongly depending on the availability of nan and inf output to prevent exceptions being raised upon invalid calculations. Maybe this should simply be considered a requirement for the use of ma. The point of my original suggestion was that it would make power work as a user might reasonably expect under sane conditions that are likely to arise in practice. Under extreme conditions, it would leave unmasked nans or infs, or would raise an exception if that were the specified handling mode for invalid calculations. Eric ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)
On Friday 09 May 2008 18:45:33 Eric Firing wrote: > I don't think the .max() part of that is right; the test needs to be > element-wise, and turned into a mask. Quite right. I was being overzealous... > It is also not clear to me that the test would actually catch all the > cases where x**b would return NaN. Oh, probably not, but it's close enough: raise an exception if you have a negative number and an exponent that is significantly different froman integer. > It seems like some strategic re-thinking may be needed in the long run, > if not immediately. There is a wide range of combinations of arguments > that will trigger invalid results, whether Inf or NaN. Mmh, I forgot about the zero case with negative integers: right now, inf is returned. Should be easy enough to make a (x The only way to > trap and mask all of these is to use masked_invalid after the > calculation, and this only works if the user has not disabled nan > output. We'll agree that's a rather quick-and-dirty patch, not a real fix... > I have not checked recently, but based on earlier strategy > discussions, I suspect that numpy.ma is already strongly depending on > the availability of nan and inf output to prevent exceptions being > raised upon invalid calculations. Maybe this should simply be > considered a requirement for the use of ma. I wouldn't say strongly. In most cases, potential NaNs/Infs are trapped beforehand. Paul, Sasha and the other original developers had introduced DomainedOperation classes that are quite useful for that. masked_invalid may be used locally in scipy.stats.mstats as a quick fix... In our case, we need for a**b: - trap the case [a==0] >>> (a>> (a<0) & (abs(b-b.astype(int))http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)
Pierre GM wrote: > On Friday 09 May 2008 18:45:33 Eric Firing wrote: >> I don't think the .max() part of that is right; the test needs to be >> element-wise, and turned into a mask. > > Quite right. I was being overzealous... > >> It is also not clear to me that the test would actually catch all the >> cases where x**b would return NaN. > > Oh, probably not, but it's close enough: raise an exception if you have a > negative number and an exponent that is significantly different froman > integer. > >> It seems like some strategic re-thinking may be needed in the long run, >> if not immediately. There is a wide range of combinations of arguments >> that will trigger invalid results, whether Inf or NaN. > > Mmh, I forgot about the zero case with negative integers: right now, inf is > returned. Should be easy enough to make a (x >> The only way to >> trap and mask all of these is to use masked_invalid after the >> calculation, and this only works if the user has not disabled nan >> output. > > We'll agree that's a rather quick-and-dirty patch, not a real fix... > >> I have not checked recently, but based on earlier strategy >> discussions, I suspect that numpy.ma is already strongly depending on >> the availability of nan and inf output to prevent exceptions being >> raised upon invalid calculations. Maybe this should simply be >> considered a requirement for the use of ma. > > I wouldn't say strongly. In most cases, potential NaNs/Infs are trapped > beforehand. Paul, Sasha and the other original developers had introduced > DomainedOperation classes that are quite useful for that. masked_invalid may > be used locally in scipy.stats.mstats as a quick fix... > > In our case, we need for a**b: > - trap the case [a==0] (a - trap the case [(a<0) and real exponent too different from an integer] (a<0) & (abs(b-b.astype(int)) > Am I missing anything else ? > ___ > 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] Power domain (was Re: bug in oldnumeric.ma)
2008/5/9 Eric Firing <[EMAIL PROTECTED]>: > It seems like some strategic re-thinking may be needed in the long run, > if not immediately. There is a wide range of combinations of arguments > that will trigger invalid results, whether Inf or NaN. The only way to > trap and mask all of these is to use masked_invalid after the > calculation, and this only works if the user has not disabled nan > output. I have not checked recently, but based on earlier strategy > discussions, I suspect that numpy.ma is already strongly depending on > the availability of nan and inf output to prevent exceptions being > raised upon invalid calculations. Maybe this should simply be > considered a requirement for the use of ma. I think in principle the right answer is to simply run whatever underlying function, and mask any NaNs or Infs in the output. This may be a problem when it comes to seterr - is the behaviour of seterr() even defined in a multithreaded context? Can it be made thread-local? Surely hardware must support thread-local floating-point error flags. If seterr() works this way, then surely the right answer in ma is to use a try/finally block to turn off exceptions and warnings, and clean up NaNs after the fact. Anne ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
Re: [Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)
Anne Archibald wrote: > 2008/5/9 Eric Firing <[EMAIL PROTECTED]>: > > >> It seems like some strategic re-thinking may be needed in the long run, >> if not immediately. > I think in principle the right answer is to simply run whatever > underlying function, and mask any NaNs or Infs in the output. Did you already rule out promoting to complex based on the type information of the power args? I see that already: >> power( -3.0 , arange(2.8, 3.2, 0.1) ) Warning: invalid value encountered in power array([ NaN, NaN, -27., NaN, NaN]) >>> power( -3.0*(1+0j) , arange(2.8, 3.2, 0.1) ) array([-17.53465227+12.73967059j, -23.00689255 +7.47539254j, -27. +0.j, -28.66039788 -9.31232777j, -27.21107252-19.77000142j]) This gives a continuous function. So, multiply by (1+0j) and mask according to the .imag parts of the answers. It would have made more sense to have power(real, real) return complex, but of course that may break some code. Jon ___ Numpy-discussion mailing list Numpy-discussion@scipy.org http://projects.scipy.org/mailman/listinfo/numpy-discussion
