BUG: fix integer roundoff corner case in Lomb-Scargle - #4950

Hopefully this will fix the allowed test failures on the numpy-dev build.

Test failure looks unrelated to this PR: appears a VM died mid-way through the test run.

The good news is that the previously failing numpy-dev test is now passing!

@astrofrog I'm happy for this to be merged but I am currently traveling and it will need to be merged onto the 1.2 branch as well.

I'm still puzzled why it raised errors before. It seemed in all the tests dy=1, so what would go wrong?
...
...
Hmm, this is clearly a bug in numpy (now at numpy/numpy#7649)

np.array(1, dtype=int) ** -2
# 0
np.array(1, dtype=float) ** -2
# 1.

but this also suggests a simpler fix for now: just set dy = 1. (i.e., add the dot) for the default dy=None case. This would seem good practice anyway.

I was thinking that this was a bit weird, so glad to hear it is a Numpy bug!

just set dy = 1. (i.e., add the dot) for the default dy=None case.

That would do nothing in the case of a user passing an integer. The solution in this PR should work in all cases, and is much cleaner than trying to enforce dy to be a float everywhere.

Oh, and I opened a numpy issue for this last night: numpy/numpy#7648 😄

good point. though I might then do 1./dy**2. But really this is a numpy bug, so I labelled it as such. If numpy is quick about, I'd use their fix...

But perhaps I'm being silly: it certainly doesn't hurt to merge this as is either.

I agree with @mhvk . Adding a ".0" will not hurt performance or anything, right? So, might as well.

No, the .0 doesn't matter, though just doing 1/a**2 is faster:

a = np.ones((100000,), dtype=int)
%timeit a ** -2
# 100 loops, best of 3: 7.21 ms per loop
%timeit a ** -2.0
# 100 loops, best of 3: 7.29 ms per loop
%timeit 1./a**2
1000 loops, best of 3: 1.1 ms per loop

But I doubt any of these differences matter compared to the overall runtime.

This is now fixed upstream, so I'm just closing this; it is better if code just works without having to worry about whether a power is float or integer!

No, please don't close this! This needs to be merged. The upstream fix will not affect old numpy releases, and there is still the possibility for this to break things here. Note that there are two separate issues happening: the broken behavior in dev on np.array(1) ** -2, and the broken behavior on all current releases with, e.g. np.int32(2) ** -2 (which is platform dependent). The fix in this PR is relevant to the latter case, which is not fixed upstream.

Again, this PR needs to be merged, and I can't think of any good reason not to merge it.

Is this not an issue only with the developer version of Numpy though?

Not necessarily, the problem is that np.array(integer) ** -2 can be zero because numpy doesn't follow the coercion rules that int ** int could be a float. It's like python2 behaviour with int / int = int. 😄 But if I understand @jakevdp correctly this may be platform-dependent.

I believe there are platforms on which failures would occur due to numpy's int-casting issues. I don't think any of the Travis platforms happen to have the right combination of features to show that failure.

Really, though, this PR changes an int to a float in a case where we already expect a cast to happen anyway... I'm not sure this volume of discussion is necessary before merge 😄

I pushed another PR that tests all methods to make sure integer and float inputs yield the same results: #4965.

@jakevdp - there are no older versions of numpy that have this bug (the PR that caused it has not been in any released version), and numpy itself now has a test to ensure this doesn't recur.

@MSeifert04 - ah, yes, you're right: we're not safe for dy = 2 in user input.

@mhvk – there are two different bugs. The one you're talking about was only in dev, the one I'm talking about is in all current releases. For example:

>>> np.__version__
'1.10.4'
>>> np.int64(2) ** -2
0.25
>>> np.int32(2) ** -2
0

This Bug really messes with my head...

NumPy integer don't always follow the normal python3 coercion rules (with power-operations) and NumPy raises now Exceptions when dividing by 0. This means the Bug wasn't even cleanly recognizable before because it just issued a Warning... and even the fix still has some problems with special cases?

I tried the example @jakevdp gave and got different results (same version on windows 64 bit python3)

np.__version__    # 1.10.4
np.array(2) ** -2 # 0
np.array(2).dtype # int32
np.int32(2) ** -2 # 0.25
np.int64(2) ** -2 # 0.25
2 ** -2           # 0.25

I'm not sure if that's relevant for this discussion but @mhvk: did they fix this difference?

@MSeifert04 – it's a subtle issue having to do with the casting rules used in binary ufuncs, and the result is highly platform dependent. See the discussion in the numpy PR with the fix: numpy/numpy#7459

OK, while I would have done the cast to real differently, that is really nitpicking, and it seems this may still avoid real issues. So, let me just merge this.

@astrofrog - not quite sure about changelog since this hasn't gone in anywhere -- this should get backported to 1.2 though!