import numpy as np
from sklearn.utils import check_random_state
try:
# astropy.timeseries is new in v3.2
from astropy.timeseries import LombScargle
except ImportError:
from astropy.stats import LombScargle
from ..utils.decorators import deprecated
from ..utils.exceptions import AstroMLDeprecationWarning
[docs]@deprecated('0.4', alternative='astropy.stats.LombScargle',
warning_type=AstroMLDeprecationWarning)
def lomb_scargle(t, y, dy, omega, generalized=True,
subtract_mean=True, significance=None):
"""
(Generalized) Lomb-Scargle Periodogram with Floating Mean
Parameters
----------
t : array_like
sequence of times
y : array_like
sequence of observations
dy : array_like
sequence of observational errors
omega : array_like
frequencies at which to evaluate p(omega)
generalized : bool
if True (default) use generalized lomb-scargle method
otherwise, use classic lomb-scargle.
subtract_mean : bool
if True (default) subtract the sample mean from the data before
computing the periodogram. Only referenced if generalized is False
significance : None or float or ndarray
if specified, then this is a list of significances to compute
for the results.
Returns
-------
p : array_like
Lomb-Scargle power associated with each frequency omega
z : array_like
if significance is specified, this gives the levels corresponding
to the desired significance (using the Scargle 1982 formalism)
Notes
-----
The algorithm is based on reference [1]_. The result for generalized=False
is given by equation 4 of this work, while the result for generalized=True
is given by equation 20.
Note that the normalization used in this reference is different from that
used in other places in the literature (e.g. [2]_). For a discussion of
normalization and false-alarm probability, see [1]_.
To recover the normalization used in Scargle [3]_, the results should
be multiplied by (N - 1) / 2 where N is the number of data points.
References
----------
.. [1] M. Zechmeister and M. Kurster, A&A 496, 577-584 (2009)
.. [2] W. Press et al, Numerical Recipes in C (2002)
.. [3] Scargle, J.D. 1982, ApJ 263:835-853
"""
# delegate to astropy's Lomb-Scargle
ls = LombScargle(t, y, dy, fit_mean=generalized, center_data=subtract_mean)
frequency = np.asarray(omega) / (2 * np.pi)
p_omega = ls.power(frequency, method='cython')
if significance is not None:
N = t.size
M = 2 * N
z = (-2.0 / (N - 1.)
* np.log(1 - (1 - np.asarray(significance)) ** (1. / M)))
return p_omega, z
else:
return p_omega
[docs]@deprecated('0.4', alternative='astropy.stats.LombScargle.false_alarm_probability',
warning_type=AstroMLDeprecationWarning)
def lomb_scargle_bootstrap(t, y, dy, omega,
generalized=True, subtract_mean=True,
N_bootstraps=100, random_state=None):
"""Use a bootstrap analysis to compute Lomb-Scargle significance
Parameters
----------
The first set of parameters are passed to the lomb_scargle algorithm
t : array_like
sequence of times
y : array_like
sequence of observations
dy : array_like
sequence of observational errors
omega : array_like
frequencies at which to evaluate p(omega)
generalized : bool
if True (default) use generalized lomb-scargle method
otherwise, use classic lomb-scargle.
subtract_mean : bool
if True (default) subtract the sample mean from the data before
computing the periodogram. Only referenced if generalized is False
Remaining parameters control the bootstrap
N_bootstraps : int
number of bootstraps
random_state : None, int, or RandomState object
random seed, or random number generator
Returns
-------
D : ndarray
distribution of the height of the highest peak
"""
random_state = check_random_state(random_state)
t = np.asarray(t)
y = np.asarray(y)
dy = np.asarray(dy) + np.zeros_like(y)
D = np.zeros(N_bootstraps)
for i in range(N_bootstraps):
ind = random_state.randint(0, len(y), len(y))
p = lomb_scargle(t, y[ind], dy[ind], omega,
generalized=generalized, subtract_mean=subtract_mean)
D[i] = p.max()
return D
def lomb_scargle_AIC(P, y, dy, n_harmonics=1):
"""Compute the AIC for a Lomb-Scargle Periodogram
Parameters
----------
P : array_like
lomb-scargle power
y : array_like
observations
dy : array_like
errors
n_harmonics : int (optional)
the number of harmonics used in the Lomb-Scargle fit. Default is 1
Returns
-------
AIC : ndarray
AIC value corresponding to values in P
"""
P, y, dy = map(np.asarray(P, y, dy))
w = 1. / dy ** 2
mu = np.dot(w, y) / w.sum()
return np.sum(((y - mu) / dy) ** 2) * P - (2 * n_harmonics + 1) * 2
def lomb_scargle_BIC(P, y, dy, n_harmonics=1):
"""Compute the BIC for a Lomb-Scargle Periodogram
Parameters
----------
P : array_like
lomb-scargle power
y : array_like
observations
dy : array_like
errors
n_harmonics : int (optional)
the number of harmonics used in the Lomb-Scargle fit. Default is 1
Returns
-------
BIC : ndarray
BIC value corresponding to values in P
"""
P, y, dy = map(np.asarray, (P, y, dy))
w = 1. / dy ** 2
mu = np.dot(w, y) / w.sum()
N = len(y)
return np.sum(((y - mu) / dy) ** 2) * P - (2 * n_harmonics + 1) * np.log(N)
[docs]@deprecated('0.4', alternative='astropy.stats.LombScargle',
warning_type=AstroMLDeprecationWarning)
def multiterm_periodogram(t, y, dy, omega, n_terms=3):
"""Perform a multiterm periodogram at each omega
This calculates the chi2 for the best-fit least-squares solution
for each frequency omega.
Parameters
----------
t : array_like
sequence of times
y : array_like
sequence of observations
dy : array_like
sequence of observational errors
omega : float or array_like
frequencies at which to evaluate p(omega)
Returns
-------
power : ndarray
P = 1. - chi2 / chi2_0
where chi2_0 is the chi-square for a simple mean fit to the data
"""
# TODO: deprecate this
ls = LombScargle(t, y, dy, nterms=n_terms)
frequency = np.asarray(omega) / (2 * np.pi)
return ls.power(frequency)
[docs]def search_frequencies(t, y, dy,
LS_func=lomb_scargle,
LS_kwargs=None,
initial_guess=25,
limit_fractions=[0.04, 0.3, 0.9, 0.99],
n_eval=10000,
n_retry=5,
n_save=50):
"""Utility Routine to find the best frequencies
To find the best frequency with a Lomb-Scargle periodogram requires
searching a large range of frequencies at a very fine resolution.
This is an iterative routine that searches progressively finer
grids to narrow-in on the best result.
Parameters
----------
t: array_like
observed times
y: array_like
observed fluxes or magnitudes
dy: array_like
observed errors on y
Other Parameters
----------------
LS_func : function
Function used to perform Lomb-Scargle periodogram. The call signature
should be LS_func(t, y, dy, omega, **kwargs)
(Default is astroML.periodogram.lomb_scargle)
LS_kwargs : dict
dictionary of keyword arguments to pass to LS_func in addition to
(t, y, dy, omega)
initial_guess : float
the initial guess of the best period
limit_fractions : array_like
the list of fractions to use when zooming in on peak possibilities.
On the i^th iteration, with f_i = limit_fractions[i], the range
probed around each candidate will be
(candidate * f_i, candidate / f_i).
n_eval : integer or list
The number of point to evaluate in the range on each iteration.
If n_eval is a list, it should have the same length as limit_fractions.
n_retry : integer or list
Number of top points to search on each iteration. If n_retry is a list,
it should have the same length as limit_fractions.
n_save : integer or list
Number of evaluations to save on each iteration.
If n_save is a list, it should have the same length as limit_fractions.
Returns
-------
omega_top, power_top: ndarrays
The saved values of omega and power. These will have size
1 + n_save * (1 + n_retry * len(limit_fractions))
as long as n_save > n_retry
"""
if LS_kwargs is None:
LS_kwargs = dict()
omega_best = [initial_guess]
power_best = LS_func(t, y, dy, omega_best, **LS_kwargs)
for (Ne, Nr, Ns, frac) in np.broadcast(n_eval, n_retry,
n_save, limit_fractions):
# make sure we explore differing regions
log_ob = np.log(omega_best)
width = 0.1 * np.log(frac)
log_ob = np.floor(-log_ob / width).astype(int)
indices = np.arange(len(log_ob))
for i in range(Nr):
if len(indices) == 0:
break
omega_try = omega_best[indices[-1]]
non_duplicates = (log_ob != log_ob[-1])
log_ob = log_ob[non_duplicates]
indices = indices[non_duplicates]
omega = np.linspace(omega_try * frac, omega_try / frac, Ne)
power = LS_func(t, y, dy, omega, **LS_kwargs)
i = np.argsort(power)[-Ns:]
power_best = np.concatenate([power_best, power[i]])
omega_best = np.concatenate([omega_best, omega[i]])
i = np.argsort(power_best)
power_best = power_best[i]
omega_best = omega_best[i]
i = np.argsort(omega_best)
return omega_best[i], power_best[i]
[docs]class MultiTermFit:
"""Multi-term Fourier fit to a light curve
Parameters
----------
omega : float
angular frequency of the fundamental mode
n_terms : int
the number of Fourier modes to use in the fit
"""
[docs] def __init__(self, omega, n_terms):
self.omega = omega
self.n_terms = n_terms
def _make_X(self, t):
t = np.asarray(t)
k = np.arange(1, self.n_terms + 1)
X = np.hstack([np.ones(t[:, None].shape),
np.sin(k * self.omega * t[:, None]),
np.cos(k * self.omega * t[:, None])])
return X
def fit(self, t, y, dy):
"""Fit multiple Fourier terms to the data
Parameters
----------
t: array_like
observed times
y: array_like
observed fluxes or magnitudes
dy: array_like
observed errors on y
Returns
-------
self :
The MultiTermFit object is returned
"""
t = np.asarray(t)
y = np.asarray(y)
dy = np.asarray(dy)
X_scaled = self._make_X(t) / dy[:, None]
y_scaled = y / dy
self.t_ = t
self.w_ = np.linalg.solve(np.dot(X_scaled.T, X_scaled),
np.dot(X_scaled.T, y_scaled))
return self
def predict(self, Nphase, return_phased_times=False, adjust_offset=True):
"""Compute the phased fit, and optionally return phased times
Parameters
----------
Nphase : int
Number of terms to use in the phased fit
return_phased_times : bool
If True, then return a phased version of the input times
adjust_offset : bool
If true, then shift results so that the minimum value is at phase 0
Returns
-------
phase, y_fit : ndarrays
The phase and y value of the best-fit light curve
phased_times : ndarray
The phased version of the training times. Returned if
return_phased_times is set to True.
"""
phase_fit = np.linspace(0, 1, Nphase + 1)[:-1]
X_fit = self._make_X(2 * np.pi * phase_fit / self.omega)
y_fit = np.dot(X_fit, self.w_)
i_offset = np.argmin(y_fit)
if adjust_offset:
y_fit = np.concatenate([y_fit[i_offset:], y_fit[:i_offset]])
if return_phased_times:
if adjust_offset:
offset = phase_fit[i_offset]
else:
offset = 0
phased_times = (self.t_ * self.omega * 0.5 / np.pi - offset) % 1
return phase_fit, y_fit, phased_times
else:
return phase_fit, y_fit