[Numpy-discussion] Power domain (was Re: bug in oldnumeric.ma)

[Pierre GM]

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<np.finfo(float).tiny) trap... 

> 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<np.finfo(float).tiny)
- trap the case [(a<0) and real exponent too different from an integer]
>>> (a<0) & (abs(b-b.astype(int))<np.finfo(float).precision)

Am I missing anything else ?