Repository navigation
Enable fit_mean=True for Lomb-Scargle periodograms - #20412
Conversation
|
Thank you for your contribution to Astropy! 🌌 This checklist is meant to remind the package maintainers who will review this pull request of some common things to look for.
|
|
This is not ready to look at yet, will finalize tomorrow |
astrofrog
changed the title
fit_mean=True for Lomb-Scargle periodograms
b042094 to
853be33
Compare
| if fit_mean and SCIPY_LT_1_15: | ||
| pytest.skip("SciPy 1.15+ required for using `fit_mean=True`") | ||
| elif errors == "full" or (errors == "partial" and normalization == "psd"): | ||
| pytest.skip("scipy method only supports uniform uncertainties dy") |
There was a problem hiding this comment.
@jakevdp suggested in https://github.com/astropy/astropy/pull/17170/changes#r1798477154 that we could pass through non-uniform errors, which is true in principle but this still fails due to a pre-existing issue on main, so I'll address in a follow-up PR to keep this one focused.
|
@dhomeier @neutrinoceros - let me know if you are happy with this approach and the result, or if there are things from your PRs that you think should be included here. |
853be33 to
ac6b4ee
Compare
neutrinoceros
left a comment
neutrinoceros
left a comment
There was a problem hiding this comment.
Thanks for reviving this Tom.
Generally this looks good to me, but I do have a couple suggestions I would like addressed before I formally approve.
the new fit_mean argument of
lombscargle_scipyis keyword-only, so existing positional callers are unaffected
this doesn't make sense to me. Did you mean "has a default value" instead of "is keyword-only" ? To be clear:
- keyword-only parameters usually have default values
- parameters with default values are usually expected to be passed as keywords
- even though the overlap is considerable in practice, none of these set is a subset to the other; one can define a keyword-only parameter without a default value, and default values can be added to any parameter (including positional-only ones), as long as every parameter that follow also has one.
(note that I'm also distinguishing "parameters", which are part of a function signature from "arguments", which are values actually passed to functions at call sites)
| # importing scipy itself; if scipy is installed without metadata (e.g. in | ||
| # some frozen applications) assume it is recent, and if it is not installed | ||
| # at all this value is never consulted since HAS_SCIPY gates all uses. | ||
| try: |
There was a problem hiding this comment.
this should use HAS_SCIPY instead of a try/except
There was a problem hiding this comment.
alternatively... why not use minversion here ?
| P_float = LombScargle(t, y, dy, **kwds).power(frequency, method=method) | ||
| P_int = LombScargle(t_int, y_int, dy_int, **kwds).power(frequency, method=method) | ||
| assert_allclose(P_float, P_int) | ||
| assert_allclose(P_int, P_float) |
There was a problem hiding this comment.
if I'm following correctly you're switching arguments to better fit the convention that the know-good value (expected) is on the second position. In this specific instance it looks surprising to me; I expect int-based math to produce more reliable results so I would expect it to serve as the reference. What is your reasoning ?
There was a problem hiding this comment.
This was done in the original PRs I combined here - but I agree it's not necessary, so reverted in 1ea6fde
There was a problem hiding this comment.
In my PR I was indeed trying to match the calling pattern of actual, desired. In this case, t_int, y_int etc. are cast (and rounded) from the original float values, so the latter should really be the reference. But I don't think it makes any difference for assert_allclose internally.
There was a problem hiding this comment.
internally, indeed it doesn't. But switching the arguments has a visible effect of changing the error message raised in case of a mismatch, and then it's preferable to have the reference on the right side
There was a problem hiding this comment.
Yes, that was probably the reason I switched a number of those, plus it's easier to grasp when looking at the code.
There was a problem hiding this comment.
For context, the whole test was added in #4965 to check integer inputs (which in the case of scipy are simply cast back to float before passing to signal.lombscargle).
There was a problem hiding this comment.
Looks like I dropped the ball on this, but just noting that afaics the other methods are also just internally converting integer inputs to float again, when multiplying in weights at the latest. A native integer implementation of a periodogram does not seem to make much sense to me anyway.
| for norm in from_to: | ||
| if norm not in NORMALIZATIONS: | ||
| raise ValueError(f"{from_normalization} is not a valid normalization") | ||
| raise ValueError(f"{norm} is not a valid normalization") |
There was a problem hiding this comment.
there's an opportunity to improve this validation here by:
- collecting all violations first (building a
list[ValueError]by incremental appends) - if the list ends up empty, do nothing
- otherwise:
- if there's exactly one exception, raise it
- otherwise: raise all exceptions under an
ExceptionGroup
There was a problem hiding this comment.
TBC, I'm not requiring this to be implemented here, but it should be easy enough that maybe we might as well
There was a problem hiding this comment.
I'd prefer not to do this, as otherwise downstream users will have to catch either ValueError or ExceptionGroup depending on the number of issues.
There was a problem hiding this comment.
not really, they'd just have to use except* ValueError instead of except ValueError
>>> try:
... raise ExceptionGroup('a group', [ValueError('1'), ValueError('2')])
... except ValueError:
... print('gotcha')
...
+ Exception Group Traceback (most recent call last):
| File "<python-input-2>", line 2, in <module>
| raise ExceptionGroup('a group', [ValueError('1'), ValueError('2')])
| ExceptionGroup: a group (2 sub-exceptions)
+-+---------------- 1 ----------------
| ValueError: 1
+---------------- 2 ----------------
| ValueError: 2
+------------------------------------
>>> try:
... raise ExceptionGroup('a group', [ValueError('1'), ValueError('2')])
... except* ValueError:
... print('gotcha')
...
gotcha
>>> try:
... raise ExceptionGroup('a group', [ValueError])
... except* ValueError:
... print('gotcha')
...
gotcha
>>> try:
... raise ValueError
... except* ValueError:
... print('gotcha')
...
gotcha
Yes I meant has a default value here |
… versions Co-authored-by: P. L. Lim <[email protected]>
…g-mean model, where the reference chi-squared is the variance about the mean rather than about zero (extracted from astropy#17600)
… now that the normalizations handle uncentered data
…ribution metadata, let lombscargle_scipy own the version guard, avoid auto-selecting scipy when method_kwds are passed, and reuse the shared normalization helpers
…ods instead of re-deriving the auto selection
Co-authored-by: Clément Robert <[email protected]>
… via minversion instead of an import-time metadata probe
…handled dy scale rather than a uniformity restriction
…n supported scipy versions and skip instead of returning
…in test_integer_inputs independent of the version skip, restore the original assert_allclose argument order, and reformat the validate_method signature
… method_kwds in validate_method to satisfy numpydoc validation
34eac7a to
83f0c59
Compare
… otherwise unchanged by this PR
|
@neutrinoceros - I think I've addressed all your review comments, would you mind taking another look? |
|
thanks. I would like us to converge on #20412 (comment) but then I can switch my review status |
…an ExceptionGroup when both are invalid
|
@neutrinoceros - I've done this now (implemented the ExceptionGroup) |
| with pytest.raises(ExceptionGroup, match="invalid normalizations") as excinfo: | ||
| convert_normalization(np.ones(3), 10, "spam", "eggs") | ||
| assert [str(exc) for exc in excinfo.value.exceptions] == [ | ||
| "spam is not a valid normalization", | ||
| "eggs is not a valid normalization", | ||
| ] |
There was a problem hiding this comment.
Much preferable to use pytest.RaisesGroup and pytest.RaisesExc here. Note that we'll need pytest 8.4 as the minimal requirement but that should be perfectly okay to bump now. I'm actually planning to require 9.0+ in January
If you need an example of how to use these see https://github.com/neutrinoceros/gpgi/blob/54f5174a43ac3d9445aa468e8643f8db5f5cdda8/tests/test_api.py#L337
neutrinoceros
left a comment
neutrinoceros
left a comment
There was a problem hiding this comment.
the test needs a small rewrite
Co-authored-by: Clément Robert <[email protected]>
|
@neutrinoceros - I've now re-written the test |
neutrinoceros
left a comment
neutrinoceros
left a comment
There was a problem hiding this comment.
thanks ! that's a go for me
|
So you want to keep treating the cast to int version as test reference? |
|
I don't feel strongly about it |
|
Yes, it's been backwards all the time without apparently creating much confusion. |
|
Just to be clear... is there a follow-up issue to be opened here? This PR is already merged. |
|
There was an unresolved discussion on the tests, but it's probably not worth opening a new PR. |
|
Maybe could keep it in mind when addressing #20412 (comment) though |
Description
With my astropy.timeseries maintainer hat on, and in the interest of unblocking two competing (and now stale) PRs to enable
fit_mean=Truein the Lomb-Scargle scipy mode, I have made this PR which is a combination of the two PRs:fit_mean=Trueinlombscargle_scipyfor Scipy 1.15+ #17170 by @dhomeierfit_mean=Truetolombscargle_scipy#17600 by @neutrinocerosand am closing the two PRs. We should now focus all discussions here.
I have used all the commits from the first except the precenter one (since scipy will deprecate this parameter), as well as the normalization fix from the second, and then applied some additional fixes on top of this. Compared to
main, the changes are as follows:fit_mean=Truefor the scipy Lomb-Scargle method by passing scipy 1.15+'sfloating_meanoption through toscipy.signal.lombscargle. On older scipy versions a clearValueErroris raisedmethod='auto'can now select the scipy backend for the default configuration (fit_mean=True, nody) - previously it was excluded whenever the mean was fitted - but never whenmethod_kwdsare passed (scipy accepts none)dyarrays for the scipy method instead of rejecting all uncertaintiesstandard/log/modelnormalizations for the floating-mean model with uncentered data, now computed via the package's sharedcompute_chi2_ref()/convert_normalization()helpers; results match theslow/cythonimplementations (to ~1e-12)fit_meanargument oflombscargle_scipyis keyword-only, so existing positional callers are unaffectedfit_mean=Truestays the default everywhere, and plainLombScargle(t, y).power(freq)returns the same results as beforefit_mean=Trueand uniform-dycombinations, comparing every method against the exactslowreference; docs (method table and narrative) updated to reflect the new capabilitiesNote that no deprecation is needed because the old scipy method raised an error for anything except an explicit
fit_mean=False, which still behaves exactly as before, so there are no existing users whose unspecifieddefault could silently change meaning.
AI Disclosure
If AI tools were used to develop this pull request, describe the tools including specific model and version, how they were used, and what content is AI generated. Otherwise enter "N/A".
I used Claude Fable to help review the changes from the two PRs and understand how to best combine them, and help identify and fix some of the issues that came up.
Merge method