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

[Eric Firing]

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<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)

You need abs(a) there, otherwise you will knock out all negative a values.


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

Did you mean int64 here? This is also assuming the output is float64, 
correct?  Then I think it will fail if the inputs are float32, because 
the underlying ufunc will need to return a float32 result, but the 
criterion was based on float64.

There are still all the cases where the combined (large) sizes of a and 
b will make a**b = Inf; I don't think that anything short of letting the 
inf be generated and then using masked_invalid will work for that.

Eric

> 
> Am I missing anything else ?
> _______________________________________________
> Numpy-discussion mailing list
> Numpy-discussion@scipy.org
> http://projects.scipy.org/mailman/listinfo/numpy-discussion