MAINT: signal.lombscargle: fix numerical backward-compat - #21697

Thanks for your patience on this. Perhaps the signal regulars will have a few suggestions given the ecosystem discussion in the matching issue.

No worries. I sympathize with the large ecosystem that makes use of scipy. But it also feels unfortunate not to default to the preferred/better option, and hope the user discovers it.

I've changed the default, and also tried to suggest in the docstring that the user might want to consider using floating_mean=True. I'm open to any further suggested modifications.

I still propose that the default should be changed eventually—with sufficient notice. I understand there are many downstream users. Because scipy is the most popular toolbox that has this functionality, and is geared towards generalists, I believe providing for a "floating-mean" (which is generally preferred and recommended in literature) is the better choice. On balance, the current status quo is more likely to produce erroneous results.

FWIW I checked that this patch resolves failures in astropy's test suite. Thank you !

thanks @adammj @DietBru @neutrinoceros

This still doesn't restore backwards compatibility for me in a library using lombscargle to identify the primary frequency component in a routine for fitting potentially randomly sampled and noisy/decaying sinusoidal data (https://github.com/OxfordIonTrapGroup/oitg/blob/master/oitg/fitting/sinusoid.py / https://github.com/OxfordIonTrapGroup/oitg/blob/master/test/fitting/test_sinusoid.py). Example:

The code might be a bit janky in terms of both numerics and implementation, but worked fine for many SciPy releases with precenter=True. To restore functionality on 1.15+, I need to manually do y - np.mean(y). floating_mean=True or False doesn't make a difference. Is this intentional?

@dnadlinger I'll take a look this weekend. Can you tell me which of your tests fail?

@adammj: Oh, great, thanks a lot! I've pushed a branch without the workaround and with a Poetry lock file at https://github.com/OxfordIonTrapGroup/oitg/tree/scipy-lombscargle-repro. You should be able to reproduce the error by running poetry run python -m unittest test/fitting/test_sinusoid.py. If you don't use Poetry, just running the command with SciPy available in a checkout should probably also work; I don't think there are external dependencies. The tests use random data, so might pass by accident from time to time (it would be prudent to use a fixed seed…), but by and large, they all fail on 1.15+, but all pass on 1.14 (or 1.15 after manually subtracting np.mean).

The code has some plumbing (FitBase) to pipe additional parameter initialisers, etc. around, but in a nutshell, https://github.com/OxfordIonTrapGroup/oitg/blob/ce05f2780d050b19784abf0db46121a570ac8823/oitg/fitting/sinusoid.py#L55-L58 previously gave an estimate of the frequency of the sinusoid close enough for curve_fit to later converge, while the new implementation can be significantly off. Just copying out the lombscarge() call and applying it to the t/y values from the unit tests should show the difference between the old and new implementations.

It might be useful to open an issue summarizing some of this to reduce the chances of it falling off our radar.

I think we're aiming for our next 1.16.x series release around Sept. 10th.

@tylerjereddy I don't think this is an issue with lombscargle.

@dnadlinger the issue arises in sinusoid.py. You're storing the mean of y (p_dict["c"] = np.mean(y) on line 59) -after- you call lombscargle on line 56. However, because you've set precenter=True, the following is run on line 249 of _spectral_py.py y -= y.mean(), which modifies y in place. I would recommend that you move p_dict["c"] = np.mean(y) to before the call to lombscargle. When you do this, all of your tests pass*.

I added an * there because you're using np.random without a seed. So, I can occasionally get 1 test (test_random_data) to fail when I run the unittests enough times. I'd recommend you add a seed to all of your unittests.

@adammj: Thanks a lot for investigating this even without a reduced test case! Nevertheless, isn't this still a breaking API change? lombscarge(t, y, precenter=True) now modifies y in place, while it didn't before. This also isn't documented in the API docs.

Technically, the revamped doc comment does state that it performs y -= y.mean(), which is in place, but still, it's a breaking API change that should probably have been avoided (especially since precenter is supposed to be a legacy API) or at least documented as such.

I disagree with the characterization that "it's a breaking API change", because the previous API never said what it did with y. I'm in the camp that just because something is "misused" doesn't make that misuse an unwritten requirement for the API. I also don't speak for scipy's maintainers.

My only other thought is that changing that line from y -= y.mean() to y = y - y.mean() would unnecessarily increase the memory requirements.

@tylerjereddy: I created #23474 to track further discussions.

@adammj: Thanks – I responded to your points in the description of #23474 as well. My view as a library developer is that silent breaking changes to "reasonable" user progams (i.e. even in the absence of explicit guarantees, what a user might reasonably expect) are the worst kind of breakage, especially in minor versions, as they are harder to track down than any acutal breaking API change. I'd rather see the flag removed outright, because then there would at least be an obvious error message (practially, of course, we'd want to deprecate it first).