Skip to content

Commit 131a74f

Browse files
pllimpllim
authored
Merge pull request #20412 from astrofrog/lombscargle-scipy-fit-mean
Enable ``fit_mean=True`` for Lomb-Scargle periodograms
2 parents 8e1c1a4 + 855f159 commit 131a74f

7 files changed

Lines changed: 137 additions & 44 deletions

File tree

‎astropy/timeseries/periodograms/lombscargle/implementations/main.py‎

Lines changed: 27 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@
1515
from .cython_impl import lombscargle_cython
1616
from .fast_impl import lombscargle_fast
1717
from .fastchi2_impl import lombscargle_fastchi2
18-
from .scipy_impl import lombscargle_scipy
18+
from .scipy_impl import _scipy_supports_floating_mean, lombscargle_scipy
1919
from .slow_impl import lombscargle_slow
2020

2121
METHODS = {
@@ -75,7 +75,16 @@ def _get_frequency_grid(frequency, assume_regular_frequency=False):
7575
return frequency[0], frequency[1] - frequency[0], len(frequency)
7676

7777

78-
def validate_method(method, dy, fit_mean, nterms, frequency, assume_regular_frequency): # ruff: ignore[PLR0917]
78+
def validate_method( # ruff: ignore[PLR0917]
79+
method,
80+
dy,
81+
fit_mean,
82+
nterms,
83+
frequency,
84+
assume_regular_frequency,
85+
*,
86+
method_kwds=None,
87+
):
7988
"""
8089
Validate the method argument, and if method='auto'
8190
choose the appropriate method.
@@ -95,6 +104,10 @@ def validate_method(method, dy, fit_mean, nterms, frequency, assume_regular_freq
95104
assume_regular_frequency : bool
96105
If True, assume that the frequencies are regularly spaced rather than
97106
checking.
107+
method_kwds : dict, optional, keyword-only
108+
Additional keywords to pass to the implementation, or None if none
109+
were specified. The scipy method accepts no additional keywords, so
110+
it is never chosen automatically when this is not None.
98111
99112
Returns
100113
-------
@@ -105,7 +118,14 @@ def validate_method(method, dy, fit_mean, nterms, frequency, assume_regular_freq
105118
prefer_fast = len(frequency) > 200 and (
106119
assume_regular_frequency or _is_regular(frequency)
107120
)
108-
prefer_scipy = "scipy" in methods and dy is None and not fit_mean
121+
# the scipy method accepts no additional keywords, so never select it
122+
# automatically when method_kwds are passed
123+
can_use_scipy = (
124+
"scipy" in methods
125+
and dy is None
126+
and method_kwds is None
127+
and (not fit_mean or _scipy_supports_floating_mean())
128+
)
109129

110130
# automatically choose the appropriate method
111131
if method == "auto":
@@ -116,7 +136,7 @@ def validate_method(method, dy, fit_mean, nterms, frequency, assume_regular_freq
116136
method = "chi2"
117137
elif prefer_fast:
118138
method = "fast"
119-
elif prefer_scipy:
139+
elif can_use_scipy:
120140
method = "scipy"
121141
else:
122142
method = "cython"
@@ -224,12 +244,12 @@ def lombscargle( # ruff: ignore[PLR0917]
224244
nterms=nterms,
225245
frequency=frequency,
226246
assume_regular_frequency=assume_regular_frequency,
247+
method_kwds=method_kwds,
227248
)
228249

229-
# scipy doesn't support dy or fit_mean=True
250+
# scipy doesn't support a dy array; lombscargle_scipy itself checks the
251+
# scipy version requirement for fit_mean=True
230252
if method == "scipy":
231-
if kwds.pop("fit_mean"):
232-
raise ValueError("scipy method does not support fit_mean=True")
233253
if dy is not None:
234254
dy = np.ravel(np.asarray(dy))
235255
if not np.allclose(dy[0], dy):

‎astropy/timeseries/periodograms/lombscargle/implementations/scipy_impl.py‎

Lines changed: 37 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,29 @@
1+
from functools import cache
2+
13
import numpy as np
24

5+
from astropy.timeseries.periodograms.lombscargle.utils import (
6+
compute_chi2_ref,
7+
convert_normalization,
8+
)
9+
from astropy.utils import minversion
310
from astropy.utils.compat.optional_deps import HAS_SCIPY
411

512

6-
def lombscargle_scipy(t, y, frequency, normalization="standard", center_data=True):
13+
@cache
14+
def _scipy_supports_floating_mean():
15+
# Deferred so that scipy is only imported and its version only probed when
16+
# the answer is needed; passing the module object lets minversion fall
17+
# back to scipy.__version__ when the distribution metadata is missing
18+
# (e.g. in some frozen applications).
19+
import scipy
20+
21+
return minversion(scipy, "1.15.0.dev0")
22+
23+
24+
def lombscargle_scipy(
25+
t, y, frequency, normalization="standard", center_data=True, *, fit_mean=False
26+
):
727
"""Lomb-Scargle Periodogram.
828
929
This is a wrapper of ``scipy.signal.lombscargle`` for computation of the
@@ -23,6 +43,11 @@ def lombscargle_scipy(t, y, frequency, normalization="standard", center_data=Tru
2343
center_data : bool, optional
2444
If True, pre-center the data by subtracting the weighted mean
2545
of the input data.
46+
fit_mean : bool, optional, keyword-only
47+
If True, include a constant offset as part of the model at each
48+
frequency. This can lead to more accurate results, especially in the
49+
case of incomplete phase coverage. Requires Scipy 1.15 and corresponds
50+
to the ``floating_mean`` argument in `scipy.signal.lombscargle`.
2651
2752
Returns
2853
-------
@@ -41,6 +66,9 @@ def lombscargle_scipy(t, y, frequency, normalization="standard", center_data=Tru
4166

4267
from scipy import signal
4368

69+
if fit_mean and not _scipy_supports_floating_mean():
70+
raise ValueError("fit_mean=True requires scipy 1.15 or later")
71+
4472
t, y = np.broadcast_arrays(t, y)
4573

4674
# Scipy requires floating-point input
@@ -56,17 +84,13 @@ def lombscargle_scipy(t, y, frequency, normalization="standard", center_data=Tru
5684
if center_data:
5785
y = y - y.mean()
5886

59-
# Note: scipy input accepts angular frequencies
60-
p = signal.lombscargle(t, y, 2 * np.pi * frequency)
87+
kwargs = {"floating_mean": True} if fit_mean else {}
88+
89+
# Note: scipy `freqs` input is in angular frequencies
90+
p = signal.lombscargle(t, y, 2 * np.pi * frequency, **kwargs)
6191

6292
if normalization == "psd":
63-
pass
64-
elif normalization == "standard":
65-
p *= 2 / (t.size * np.mean(y**2))
66-
elif normalization == "log":
67-
p = -np.log(1 - 2 * p / (t.size * np.mean(y**2)))
68-
elif normalization == "model":
69-
p /= 0.5 * t.size * np.mean(y**2) - p
70-
else:
71-
raise ValueError(f"normalization='{normalization}' not recognized")
72-
return p
93+
return p
94+
95+
chi2_ref = compute_chi2_ref(y, center_data=center_data, fit_mean=fit_mean)
96+
return convert_normalization(p, t.size, "psd", normalization, chi2_ref=chi2_ref)

‎astropy/timeseries/periodograms/lombscargle/tests/test_lombscargle.py‎

Lines changed: 24 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,9 @@
99
from astropy.timeseries.periodograms.lombscargle._testing import (
1010
assert_not_strictly_equal,
1111
)
12+
from astropy.timeseries.periodograms.lombscargle.implementations.scipy_impl import (
13+
_scipy_supports_floating_mean,
14+
)
1215

1316
ALL_METHODS = LombScargle.available_methods
1417
ALL_METHODS_NO_AUTO = [method for method in ALL_METHODS if method != "auto"]
@@ -70,8 +73,16 @@ def test_autofrequency(
7073
def test_all_methods(
7174
data, method, center_data, fit_mean, errors, with_units, normalization
7275
):
73-
if method == "scipy" and (fit_mean or errors != "none"):
74-
return
76+
if method == "scipy":
77+
if fit_mean and not _scipy_supports_floating_mean():
78+
pytest.skip("SciPy 1.15+ required for using `fit_mean=True`")
79+
if errors == "full":
80+
pytest.skip("scipy method only supports uniform uncertainties dy")
81+
elif errors == "partial" and normalization == "psd":
82+
pytest.skip(
83+
"scipy method drops the uniform dy, whose scale the psd "
84+
"normalization depends on"
85+
)
7586

7687
t, y, dy = data
7788
frequency = 0.8 + 0.01 * np.arange(40)
@@ -100,7 +111,10 @@ def test_all_methods(
100111
fit_mean=fit_mean,
101112
normalization=normalization,
102113
)
103-
P_expected = ls.power(frequency)
114+
# Compare against the exact pure-python reference implementation, which
115+
# supports every combination exercised here
116+
reference = "cython" if method == "slow" else "slow"
117+
P_expected = ls.power(frequency, method=reference)
104118

105119
# don't use the lagrangian approximation here; we'll test this elsewhere
106120
if method in FAST_METHODS:
@@ -115,7 +129,7 @@ def test_all_methods(
115129
else:
116130
assert not hasattr(P_method, "unit")
117131

118-
assert_quantity_allclose(P_expected, P_method)
132+
assert_quantity_allclose(P_method, P_expected, rtol=1e-10)
119133

120134

121135
@pytest.mark.parametrize("assume_regular_frequency", [True, False])
@@ -164,8 +178,11 @@ def test_cython_irregular_frequency(data, fit_mean, center_data):
164178
def test_integer_inputs(
165179
data, method, center_data, fit_mean, with_errors, normalization
166180
):
167-
if method == "scipy" and (fit_mean or with_errors):
168-
return
181+
if method == "scipy":
182+
if fit_mean and not _scipy_supports_floating_mean():
183+
pytest.skip("SciPy 1.15+ required for using `fit_mean=True`")
184+
if with_errors:
185+
pytest.skip("scipy method only supports uniform uncertainties dy")
169186

170187
t, y, dy = data
171188

@@ -234,7 +251,7 @@ def test_nterms_methods(
234251

235252
P_method = ls.power(frequency, method=method, **kwds)
236253

237-
assert_allclose(P_expected, P_method, rtol=1e-7, atol=1e-25)
254+
assert_allclose(P_method, P_expected, rtol=1e-7, atol=1e-25)
238255

239256

240257
@pytest.mark.parametrize("method", FAST_METHODS)

‎astropy/timeseries/periodograms/lombscargle/tests/test_utils.py‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -43,3 +43,23 @@ def test_convert_normalization(norm_in, norm_out, data):
4343
chi2_ref=compute_chi2_ref(y, dy),
4444
)
4545
assert_allclose(power_in_converted, power_out)
46+
47+
48+
@pytest.mark.parametrize(
49+
"norm_in, norm_out, invalid",
50+
[("spam", "standard", "spam"), ("standard", "eggs", "eggs")],
51+
)
52+
def test_convert_normalization_invalid(norm_in, norm_out, invalid):
53+
with pytest.raises(ValueError, match=f"^{invalid} is not a valid normalization$"):
54+
convert_normalization(np.ones(3), 10, norm_in, norm_out)
55+
56+
57+
def test_convert_normalization_both_invalid():
58+
59+
expected_exceptions = [
60+
pytest.RaisesExc(ValueError, match="spam is not a valid normalization"),
61+
pytest.RaisesExc(ValueError, match="eggs is not a valid normalization"),
62+
]
63+
64+
with pytest.RaisesGroup(*expected_exceptions, match="invalid normalizations"):
65+
convert_normalization(np.ones(3), 10, "spam", "eggs")

‎astropy/timeseries/periodograms/lombscargle/utils.py‎

Lines changed: 9 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -62,9 +62,15 @@ def convert_normalization(Z, N, from_normalization, to_normalization, chi2_ref=N
6262
Z = np.asarray(Z)
6363
from_to = (from_normalization, to_normalization)
6464

65-
for norm in from_to:
66-
if norm not in NORMALIZATIONS:
67-
raise ValueError(f"{from_normalization} is not a valid normalization")
65+
errors = [
66+
ValueError(f"{norm} is not a valid normalization")
67+
for norm in from_to
68+
if norm not in NORMALIZATIONS
69+
]
70+
if len(errors) == 1:
71+
raise errors[0]
72+
elif errors:
73+
raise ExceptionGroup("multiple invalid normalizations", errors)
6874

6975
if from_normalization == to_normalization:
7076
return Z

‎docs/changes/timeseries/20412.feature.rst‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,5 @@
1+
The ``scipy`` method of `~astropy.timeseries.LombScargle` now supports the
2+
``fit_mean=True`` setting with newer Scipy versions (1.15+), equivalent to its
3+
``floating_mean`` option. As a result, ``method='auto'`` may now select the
4+
``scipy`` method in cases where it previously used ``cython``, which can
5+
change results at the floating-point rounding level.

‎docs/timeseries/lombscargle.rst‎

Lines changed: 15 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -742,9 +742,10 @@ scaling is approximately :math:`O[NM]` for :math:`N` data points and
742742

743743
The ``scipy`` method wraps the C implementation of the original Lomb-Scargle
744744
periodogram which is available in :func:`scipy.signal.lombscargle`. This is
745-
slightly faster than the ``slow`` method, but does not allow for errors in
746-
data or extensions such as the floating mean. The scaling is approximately
747-
:math:`O[NM]` for :math:`N` data points and :math:`M` frequencies.
745+
slightly faster than the ``slow`` method, but only allows for uniform errors
746+
in the data, and the floating mean (``fit_mean=True``) requires Scipy 1.15 or
747+
later. The scaling is approximately :math:`O[NM]` for :math:`N` data points
748+
and :math:`M` frequencies.
748749

749750
``method='fast'``
750751
-----------------
@@ -785,17 +786,17 @@ Summary
785786

786787
The following table summarizes the features of the above algorithms:
787788

788-
============== ============================ ============= =============== ========
789-
Method Computational Observational Bias Term Multiple
790-
Scaling Uncertainties (Floating Mean) Terms
791-
============== ============================ ============= =============== ========
792-
``"slow"`` :math:`O[NM]` Yes Yes No
793-
``"cython"`` :math:`O[NM]` Yes Yes No
794-
``"scipy"`` :math:`O[NM]` No No No
795-
``"fast"`` :math:`O[N\log M]` Yes Yes No
796-
``"chi2"`` :math:`O[n_fNM]` Yes Yes Yes
797-
``"fastchi2"`` :math:`O[n_f(M + N\log M)]` Yes Yes Yes
798-
============== ============================ ============= =============== ========
789+
============== ============================ ============= ================= ========
790+
Method Computational Observational Bias Term Multiple
791+
Scaling Uncertainties (Floating Mean) Terms
792+
============== ============================ ============= ================= ========
793+
``"slow"`` :math:`O[NM]` Yes Yes No
794+
``"cython"`` :math:`O[NM]` Yes Yes No
795+
``"scipy"`` :math:`O[NM]` Uniform only Yes (Scipy 1.15+) No
796+
``"fast"`` :math:`O[N\log M]` Yes Yes No
797+
``"chi2"`` :math:`O[n_fNM]` Yes Yes Yes
798+
``"fastchi2"`` :math:`O[n_f(M + N\log M)]` Yes Yes Yes
799+
============== ============================ ============= ================= ========
799800

800801
In the Computational Scaling column, :math:`N` is the number of data points,
801802
:math:`M` is the number of frequencies, and :math:`n_f` is the number of

0 commit comments

Comments
 (0)