Repository navigation
LombScargle occasionally returns all-NaN powers #10070
Description
Description
I want to use LombScargle to calculate the spectral window function of my observations, but sometimes I get an array of NaNs returned for powers. I'm not sure what governs the conditions for this error, but so far I've only encountered it when the input y grid is populated with a single value.
Expected behavior
Input an array t and a same-sized array y where all y-values are 1, get back some powers corresponding to the window function.
Actual behavior
Sometimes this works as expected. Other times, the output powers are an array of NaNs and the following warning is raised:
/Users/mbedell/miniconda3/envs/exoplanet/lib/python3.8/site-packages/astropy/timeseries/periodograms/lombscargle/implementations/fast_impl.py:125: RuntimeWarning: invalid value encountered in true_divide
power /= YY
Steps to Reproduce
Here is a snippet of code that reproduces the error:
from astropy.timeseries import LombScargle
np.random.seed(42)
N = 180
xs = np.arange(N) + np.random.normal(0., 0.1, N)
ys = np.ones_like(xs)
frequency, power = LombScargle(xs, ys).autopower(method='fast')
np.isfinite(power).any()(should raise a warning and return False)
Changes to the above code that do not fix the error:
- changing autopower method (I have tried 'slow' and 'chi2' with same results)
- changing the random seed/rerunning the randomization
- specifying a frequency grid as an input
- making the
xsgrid perfectly uniformly spaced (no random perturbations)
Changes to the above code that do fix the error and return sensible powers:
- changing
Nto nearly any other integer (I haven't tried an exhaustive list but for example 179 or 181 work just fine, as do other random ones I've tried. However,N = 180is definitely not the only number that can reproduce this bug -N = 256behaves similarly.) - scaling
ysby some factor - this fixes the error if and only if the factor is not a multiple of 180, as far as I can tell. That is,ys = np.ones_like(xs)/100fixes the error, butys = np.ones_like(xs)/180does not. Again, I haven't tried factors exhaustively, so take that theory with a grain of salt. - making the values of
ysdifferent from each other, e.g.:
ys = np.ones_like(xs) * np.random.normal(1., 0.01, N)
System Details
macOS-10.15.3-x86_64-i386-64bit
Python 3.8.2 | packaged by conda-forge | (default, Mar 23 2020, 17:55:48)
[Clang 9.0.1 ]
Numpy 1.18.1
astropy 4.0
Scipy 1.4.1
Activity
pllim commentedon Mar 24, 2020
@bsipocz , you are familiar with this functionality, right?
dhomeier commentedon Apr 4, 2020
@megbedell you've certainly spotted some mysterious behaviour here, but I don't think your attempts with a rescaled
ysreturn the desired result either - trying this with different choices ofmethodI found constant power of 0.0 forslow, near-constant 1.0 (!) forchi2, and a very noisy spectrum that showed indeed hints of the windowing function (peaks at integer frequencies), but diverging towards low frequencies.The
slowresult is indeed the correct one for that case, sinceLombScargleby default is using the auto-centred version of the spectrum, i.e. Eq. 2 of VanderPlas & Ivecić 2015 withy_k - <y>, i.e. all zeros in the numerator in the case of your constantys!The NaNs are a result of the default normalisation over the weighted square of the (centred) input
YY; apparently you are just getting nonzero values thanks to numerical noise for some scaling factors, e.g.You can avoid the NaNs by selecting
LombScargle(xs, ys, normalization='psd'), which does not divide by any power ofY, but what you really want is probably a nonzero output of the power spectrum that does give you the desired signal of the window function:This seems to produce consistent and sensible spectra for all methods.
I think for the other cases it would be reasonable to add a check for
YY != 0.0and either raise an exception, or a warning and return the power directly (which should be all zero then anyway).