Source code for lightkurve.periodogram

"""Defines the Periodogram class and associated tools."""
from __future__ import division, print_function

import copy
import logging
import math
import re
import warnings

import numpy as np
from matplotlib import pyplot as plt

import astropy
from astropy.table import Table
from astropy import units as u
from astropy.units import cds
from astropy.convolution import convolve, Box1DKernel
from astropy.time import Time

# LombScargle was moved from astropy.stats to astropy.timeseries in AstroPy v3.2
try:
    from astropy.timeseries import LombScargle
    from astropy.timeseries import implementations  # for .main._is_regular
except ImportError:
    from astropy.stats import LombScargle
    from astropy.stats.lombscargle import implementations


from . import MPLSTYLE
from .utils import LightkurveWarning, validate_method
from .lightcurve import LightCurve

log = logging.getLogger(__name__)

__all__ = ["Periodogram", "LombScarglePeriodogram", "BoxLeastSquaresPeriodogram"]


[docs]class Periodogram(object): """Generic class to represent a power spectrum (frequency vs power data). The Periodogram class represents a power spectrum, with values of frequency on the x-axis (in any frequency units) and values of power on the y-axis (in units of flux^2 / [frequency units]). Attributes ---------- frequency : `~astropy.units.Quantity` Array of frequencies as an AstroPy Quantity object. power : `~astropy.units.Quantity` Array of power-spectral-densities. The Quantity must have units of `flux^2 / freq_unit`, where freq_unit is the unit of the frequency attribute. nyquist : float The Nyquist frequency of the lightcurve. In units of freq_unit, where freq_unit is the unit of the frequency attribute. label : str Human-friendly object label, e.g. "KIC 123456789". targetid : str Identifier of the target. default_view : "frequency" or "period" Should plots be shown in frequency space or period space by default? meta : dict Free-form metadata associated with the Periodogram. """ frequency = None """The array of frequency values.""" power = None """The array of power values."""
[docs] def __init__( self, frequency, power, nyquist=None, label=None, targetid=None, default_view="frequency", meta={}, ): # Input validation if not isinstance(frequency, u.quantity.Quantity): raise ValueError("frequency must be an `astropy.units.Quantity` object.") if not isinstance(power, u.quantity.Quantity): raise ValueError("power must be an `astropy.units.Quantity` object.") # Frequency must have frequency units try: frequency.to(u.Hz) except u.UnitConversionError: raise ValueError("Frequency must be in units of 1/time.") # Frequency and power must have sensible shapes if frequency.shape[0] <= 1: raise ValueError("frequency and power must have a length greater than 1.") if frequency.shape != power.shape: raise ValueError("frequency and power must have the same length.") self.frequency = frequency self.power = power self.nyquist = nyquist self.label = label self.targetid = targetid self.default_view = self._validate_view(default_view) self.meta = meta
def _validate_view(self, view): """Verifies whether `view` is is one of {"frequency", "period"} and raises a helpful `ValueError` if not. """ if view is None and hasattr(self, "default_view"): view = self.default_view return validate_method(view, ["frequency", "period"]) def _is_evenly_spaced(self): """Returns true if the values in ``frequency`` are evenly spaced. This helper method exists because some features, such as ``smooth()``, ``estimate_numax()``, and ``estimate_deltanu()``, require a grid of evenly-spaced frequencies. """ # verify that the first differences are all equal freqdiff = np.diff(self.frequency.value) if np.allclose(freqdiff[0], freqdiff): return True return False @property def period(self): """The array of periods, i.e. 1/frequency.""" return 1.0 / self.frequency @property def max_power(self): """Power of the highest peak in the periodogram.""" return np.nanmax(self.power) @property def frequency_at_max_power(self): """Frequency value corresponding to the highest peak in the periodogram.""" return self.frequency[np.nanargmax(self.power)] @property def period_at_max_power(self): """Period value corresponding to the highest peak in the periodogram.""" return 1.0 / self.frequency_at_max_power
[docs] def bin(self, binsize=10, method="mean"): """Bins the power spectrum. Parameters ---------- binsize : int The factor by which to bin the power spectrum, in the sense that the power spectrum will be smoothed by taking the mean in bins of size N / binsize, where N is the length of the original frequency array. Defaults to 10. method : str, one of 'mean' or 'median' Method to use for binning. Default is 'mean'. Returns ------- binned_periodogram : a `Periodogram` object Returns a new `Periodogram` object which has been binned. """ # Input validation if binsize < 1: raise ValueError("binsize must be larger than or equal to 1") method = validate_method(method, ["mean", "median"]) m = int(len(self.power) / binsize) # length of the binned arrays if method == "mean": binned_freq = self.frequency[: m * binsize].reshape((m, binsize)).mean(1) binned_power = self.power[: m * binsize].reshape((m, binsize)).mean(1) elif method == "median": binned_freq = np.nanmedian( self.frequency[: m * binsize].reshape((m, binsize)), axis=1 ) binned_power = np.nanmedian( self.power[: m * binsize].reshape((m, binsize)), axis=1 ) binned_pg = self.copy() binned_pg.frequency = binned_freq binned_pg.power = binned_power return binned_pg
[docs] def smooth(self, method="boxkernel", filter_width=0.1): """Smooths the power spectrum using the 'boxkernel' or 'logmedian' method. If `method` is set to 'boxkernel', this method will smooth the power spectrum by convolving with a numpy Box1DKernel with a width of `filter_width`, where `filter width` is in units of frequency. This is best for filtering out noise while maintaining seismic mode peaks. This method requires the Periodogram to have an evenly spaced grid of frequencies. A `ValueError` exception will be raised if this is not the case. If `method` is set to 'logmedian', it smooths the power spectrum using a moving median which moves across the power spectrum in a steps of log10(x0) + 0.5 * filter_width where `filter width` is in log10(frequency) space. This is best for estimating the noise background, as it filters over the seismic peaks. Periodograms that are unsmoothed have multiplicative noise that is distributed as chi squared 2 degrees of freedom. This noise distribution has a well defined mean and median but the two are not equivalent. The mean of a chi squared 2 dof distribution is 2, but the median is 2(8/9)**3. (see https://en.wikipedia.org/wiki/Chi-squared_distribution) In order to maintain consistency between 'boxkernel' and 'logmedian' a correction factor of (8/9)**3 is applied to (i.e., the median is divided by the factor) to the median values. In addition to consistency with the 'boxkernel' method, the correction of the median values is useful when applying the periodogram flatten method. The flatten method divides the periodgram by the smoothed periodogram using the 'logmedian' method. By appyling the correction factor we follow asteroseismic convention that the signal-to-noise power has a mean value of unity. (note the signal-to-noise power is really the signal plus noise divided by the noise and hence should be unity in the absence of any signal) Parameters ---------- method : str, one of 'boxkernel' or 'logmedian' The smoothing method to use. Defaults to 'boxkernel'. filter_width : float If `method` = 'boxkernel', this is the width of the smoothing filter in units of frequency. If method = `logmedian`, this is the width of the smoothing filter in log10(frequency) space. Returns ------- smoothed_pg : `Periodogram` object Returns a new `Periodogram` object in which the power spectrum has been smoothed. """ method = validate_method(method, ["boxkernel", "logmedian"]) if method == "boxkernel": if filter_width <= 0.0: raise ValueError( "the `filter_width` parameter must be " "larger than 0 for the 'boxkernel' method." ) try: filter_width = u.Quantity(filter_width, self.frequency.unit) except u.UnitConversionError: raise ValueError( "the `filter_width` parameter must have " "frequency units." ) # Check to see if we have a grid of evenly spaced periods instead. if not self._is_evenly_spaced(): raise ValueError( "the 'boxkernel' method requires the periodogram " "to have a grid of evenly spaced frequencies." ) fs = np.mean(np.diff(self.frequency)) box_kernel = Box1DKernel(math.ceil((filter_width / fs).value)) smooth_power = convolve(self.power.value, box_kernel) smooth_pg = self.copy() smooth_pg.power = u.Quantity(smooth_power, self.power.unit) return smooth_pg if method == "logmedian": if isinstance(filter_width, astropy.units.quantity.Quantity): raise ValueError( "the 'logmedian' method requires a dimensionless " "value for `filter_width` in log10(frequency) space." ) count = np.zeros(len(self.frequency.value), dtype=int) bkg = np.zeros_like(self.frequency.value) x0 = np.log10(self.frequency[0].value) corr_factor = (8.0 / 9.0) ** 3 while x0 < np.log10(self.frequency[-1].value): m = np.abs(np.log10(self.frequency.value) - x0) < filter_width if len(bkg[m] > 0): bkg[m] += np.nanmedian(self.power[m].value) / corr_factor count[m] += 1 x0 += 0.5 * filter_width bkg /= count smooth_pg = self.copy() smooth_pg.power = u.Quantity(bkg, self.power.unit) return smooth_pg
[docs] def plot( self, scale="linear", ax=None, xlabel=None, ylabel=None, title="", style="lightkurve", view=None, unit=None, **kwargs ): """Plots the Periodogram. Parameters ---------- scale: str Set x,y axis to be "linear" or "log". Default is linear. ax : `~matplotlib.axes.Axes` A matplotlib axes object to plot into. If no axes is provided, a new one will be generated. xlabel : str Plot x axis label ylabel : str Plot y axis label title : str Plot set_title style : str Path or URL to a matplotlib style file, or name of one of matplotlib's built-in stylesheets (e.g. 'ggplot'). Lightkurve's custom stylesheet is used by default. view : str {'frequency', 'period'}. Default 'frequency'. If 'frequency', x-axis units will be frequency. If 'period', the x-axis units will be period and 'log' scale. kwargs : dict Dictionary of arguments to be passed to `matplotlib.pyplot.plot`. Returns ------- ax : `~matplotlib.axes.Axes` The matplotlib axes object. """ if isinstance(unit, u.quantity.Quantity): unit = unit.unit view = self._validate_view(view) if unit is None: unit = self.frequency.unit if view == "period": unit = self.period.unit if style is None or style == "lightkurve": style = MPLSTYLE if ylabel is None: ylabel = "Power" if self.power.unit.to_string() != "": unit_label = self.power.unit.to_string("latex") # The line below is a workaround for AstroPy bug #9218. # It can be removed once the fix for that issue is widespread. # See https://github.com/astropy/astropy/pull/9218 unit_label = re.sub( r"\^{([^}]+)}\^{([^}]+)}", r"^{\g<1>^{\g<2>}}", unit_label ) ylabel += " [{}]".format(unit_label) # This will need to be fixed with housekeeping. Self.label currently doesnt exist. if ("label" not in kwargs) and ("label" in dir(self)): kwargs["label"] = self.label with plt.style.context(style): if ax is None: fig, ax = plt.subplots() # Plot frequency and power if view == "frequency": ax.plot(self.frequency.to(unit), self.power, **kwargs) if xlabel is None: xlabel = "Frequency [{}]".format(unit.to_string("latex")) elif view == "period": ax.plot(self.period.to(unit), self.power, **kwargs) if xlabel is None: xlabel = "Period [{}]".format(unit.to_string("latex")) ax.set_xlabel(xlabel) ax.set_ylabel(ylabel) # Show the legend if labels were set legend_labels = ax.get_legend_handles_labels() if np.sum([len(a) for a in legend_labels]) != 0: ax.legend(loc="best") ax.set_yscale(scale) ax.set_xscale(scale) ax.set_title(title) return ax
[docs] def flatten(self, method="logmedian", filter_width=0.01, return_trend=False): """Estimates the Signal-To-Noise (SNR) spectrum by dividing out an estimate of the noise background. This method divides the power spectrum by a background estimated using a moving filter in log10 space by default. For details on the `method` and `filter_width` parameters, see `Periodogram.smooth()` Dividing the power through by the noise background produces a spectrum with no units of power. Since the signal is divided through by a measure of the noise, we refer to this as a `Signal-To-Noise` spectrum. Parameters ---------- method : str, one of 'boxkernel' or 'logmedian' Background estimation method passed on to `Periodogram.smooth()`. Defaults to 'logmedian'. filter_width : float If `method` = 'boxkernel', this is the width of the smoothing filter in units of frequency. If method = `logmedian`, this is the width of the smoothing filter in log10(frequency) space. return_trend : bool If True, then the background estimate, alongside the SNR spectrum, will be returned. Returns ------- snr_spectrum : `Periodogram` object Returns a periodogram object where the power is an estimate of the signal-to-noise of the spectrum, creating by dividing the powers with a simple estimate of the noise background using a smoothing filter. bkg : `Periodogram` object The estimated power spectrum of the background noise. This is only returned if `return_trend = True`. """ bkg = self.smooth(method=method, filter_width=filter_width) snr_pg = self / bkg.power snr = SNRPeriodogram( snr_pg.frequency, snr_pg.power, nyquist=self.nyquist, targetid=self.targetid, label=self.label, meta=self.meta, ) if return_trend: return snr, bkg return snr
[docs] def to_table(self): """Exports the Periodogram as an Astropy Table. Returns ------- table : `~astropy.table.Table` object An AstroPy Table with columns 'frequency', 'period', and 'power'. """ return Table( data=(self.frequency, self.period, self.power), names=("frequency", "period", "power"), meta=self.meta, )
[docs] def copy(self): """Returns a copy of the Periodogram object. This method uses the `copy.deepcopy` function to ensure that all objects stored within the Periodogram are copied. Returns ------- pg_copy : Periodogram A new `Periodogram` object which is a copy of the original. """ return copy.deepcopy(self)
def __repr__(self): return "Periodogram(ID: {})".format(self.label) def __getitem__(self, key): copy_self = self.copy() copy_self.frequency = self.frequency[key] copy_self.power = self.power[key] return copy_self def __add__(self, other): copy_self = self.copy() copy_self.power = copy_self.power + u.Quantity(other, self.power.unit) return copy_self def __radd__(self, other): return self.__add__(other) def __sub__(self, other): return self.__add__(-other) def __rsub__(self, other): copy_self = self.copy() copy_self.power = other - copy_self.power return copy_self def __mul__(self, other): copy_self = self.copy() copy_self.power = other * copy_self.power return copy_self def __rmul__(self, other): return self.__mul__(other) def __truediv__(self, other): return self.__mul__(1.0 / other) def __rtruediv__(self, other): copy_self = self.copy() copy_self.power = other / copy_self.power return copy_self def __div__(self, other): return self.__truediv__(other) def __rdiv__(self, other): return self.__rtruediv__(other)
[docs] def show_properties(self): """Prints a summary of the non-callable attributes of the Periodogram object. Prints in order of type (ints, strings, lists, arrays and others). Prints in alphabetical order. """ attrs = {} for attr in dir(self): if not attr.startswith("_"): res = getattr(self, attr) if callable(res): continue if isinstance(res, astropy.units.quantity.Quantity): unit = res.unit res = res.value attrs[attr] = {"res": res} attrs[attr]["unit"] = unit.to_string() else: attrs[attr] = {"res": res} attrs[attr]["unit"] = "" if attr == "hdu": attrs[attr] = {"res": res, "type": "list"} for idx, r in enumerate(res): if idx == 0: attrs[attr]["print"] = "{}".format(r.header["EXTNAME"]) else: attrs[attr]["print"] = "{}, {}".format( attrs[attr]["print"], "{}".format(r.header["EXTNAME"]) ) continue if isinstance(res, int): attrs[attr]["print"] = "{}".format(res) attrs[attr]["type"] = "int" elif isinstance(res, float): attrs[attr]["print"] = "{}".format(np.round(res, 4)) attrs[attr]["type"] = "float" elif isinstance(res, np.ndarray): attrs[attr]["print"] = "array {}".format(res.shape) attrs[attr]["type"] = "array" elif isinstance(res, list): attrs[attr]["print"] = "list length {}".format(len(res)) attrs[attr]["type"] = "list" elif isinstance(res, str): if res == "": attrs[attr]["print"] = "{}".format("None") else: attrs[attr]["print"] = "{}".format(res) attrs[attr]["type"] = "str" elif attr == "wcs": attrs[attr]["print"] = "astropy.wcs.wcs.WCS" attrs[attr]["type"] = "other" else: attrs[attr]["print"] = "{}".format(type(res)) attrs[attr]["type"] = "other" output = Table( names=["Attribute", "Description", "Units"], dtype=[object, object, object] ) idx = 0 types = ["int", "str", "float", "list", "array", "other"] for typ in types: for attr, dic in attrs.items(): if dic["type"] == typ: output.add_row([attr, dic["print"], dic["unit"]]) idx += 1 print("lightkurve.Periodogram properties:") output.pprint(max_lines=-1, max_width=-1)
[docs] def to_seismology(self, **kwargs): """Returns a `~lightkurve.seismology.Seismology` object to analyze the periodogram. Returns ------- seismology : `~lightkurve.seismology.Seismology` Helper object to run asteroseismology methods. """ from .seismology import Seismology return Seismology(self)
class SNRPeriodogram(Periodogram): """Defines a Signal-to-Noise Ratio (SNR) Periodogram class. This class is nearly identical to the standard :class:`Periodogram` class, but has different plotting defaults. """ def __init__(self, *args, **kwargs): super(SNRPeriodogram, self).__init__(*args, **kwargs) def __repr__(self): return "SNRPeriodogram(ID: {})".format(self.label) def plot(self, **kwargs): """Plot the SNR spectrum using matplotlib's `plot` method. See `Periodogram.plot` for details on the accepted arguments. Parameters ---------- kwargs : dict Dictionary of arguments ot be passed to `Periodogram.plot`. Returns ------- ax : `~matplotlib.axes.Axes` The matplotlib axes object. """ ax = super(SNRPeriodogram, self).plot(**kwargs) if "ylabel" not in kwargs: ax.set_ylabel("Signal to Noise Ratio (SNR)") return ax class LombScarglePeriodogram(Periodogram): """Subclass of :class:`Periodogram <lightkurve.periodogram.Periodogram>` representing a power spectrum generated using the Lomb Scargle method. """ def __init__(self, *args, **kwargs): self._LS_object = kwargs.pop("ls_obj", None) self.nterms = kwargs.pop("nterms", 1) self.ls_method = kwargs.pop("ls_method", "fastchi2") super(LombScarglePeriodogram, self).__init__(*args, **kwargs) def __repr__(self): return "LombScarglePeriodogram(ID: {})".format(self.label)
[docs] @staticmethod def from_lightcurve( lc, minimum_frequency=None, maximum_frequency=None, minimum_period=None, maximum_period=None, frequency=None, period=None, nterms=1, nyquist_factor=1, oversample_factor=None, freq_unit=None, normalization="amplitude", ls_method="fast", **kwargs ): """Creates a `Periodogram` from a LightCurve using the Lomb-Scargle method. By default, the periodogram will be created for a regular grid of frequencies from one frequency separation to the Nyquist frequency, where the frequency separation is determined as 1 / the time baseline. The min frequency and/or max frequency (or max period and/or min period) can be passed to set custom limits for the frequency grid. Alternatively, the user can provide a custom regular grid using the `frequency` parameter or a custom regular grid of periods using the `period` parameter. The sampling of the spectrum can be changed using the `oversample_factor` parameter. An oversampled spectrum (oversample_factor > 1) is useful for displaying the full details of the spectrum, allowing the frequencies and amplitudes to be measured directly from the plot itself, with no fitting required. This is recommended for most applications, with a value of 5 or 10. On the other hand, an oversample_factor of 1 means the spectrum is critically sampled, where every point in the spectrum is independent of the others. This may be used when Lorentzians are to be fitted to modes in the power spectrum, in cases where the mode lifetimes are shorter than the time-base of the data (which is sometimes the case for solar-like oscillations). An oversample_factor of 1 is suitable for these stars because the modes are usually fully resolved. That is, the power from each mode is spread over a range of frequencies due to damping. Hence, any small error from measuring mode frequencies by taking the maximum of the peak is negligible compared with the intrinsic linewidth of the modes. The `normalization` parameter will normalize the spectrum to either power spectral density ("psd") or amplitude ("amplitude"). Users doing asteroseismology on classical pulsators (e.g. delta Scutis) typically prefer `normalization="amplitude"` because "amplitude" has higher dynamic range (high and low peaks visible simultaneously), and we often want to read off amplitudes from the plot. If `normalization="amplitude"`, the default value for `oversample_factor` is set to 5 and `freq_unit` is 1/day. Alternatively, users doing asteroseismology on solar-like oscillators tend to prefer `normalization="psd"` because power density has a scaled axis that depends on the length of the observing time, and is used when we are interested in noise levels (e.g. granulation) and are looking at damped oscillations. If `normalization="psd"`, the default value for `oversample_factor` is set to 1 and `freq_unit` is set to microHz. Default values of `freq_unit` and `oversample_factor` can be overridden. See Appendix A of Kjeldsen & Bedding, 1995 for a full discussion of normalization and measurement of oscillation amplitudes (http://adsabs.harvard.edu/abs/1995A%26A...293...87K). The parameter nterms controls how many Fourier terms are used in the model. Setting the Nyquist_factor to be greater than 1 will sample the space beyond the Nyquist frequency, which may introduce aliasing. The `freq_unit` parameter allows a request for alternative units in frequency space. By default frequency is in (1/day) and power in (amplitude). Asteroseismologists for example may want frequency in (microHz) in which case they would pass `freq_unit=u.microhertz`. By default this method uses the LombScargle 'fast' method, which assumes a regular grid. If a regular grid of periods (i.e. an irregular grid of frequencies) it will use the 'slow' method. If nterms > 1 is passed, it will use the 'fastchi2' method for regular grids, and 'chi2' for irregular grids. Caution: this method assumes that the LightCurve's time (lc.time) is given in units of days. Parameters ---------- lc : LightCurve object The LightCurve from which to compute the Periodogram. minimum_frequency : float If specified, use this minimum frequency rather than one over the time baseline. maximum_frequency : float If specified, use this maximum frequency rather than nyquist_factor times the nyquist frequency. minimum_period : float If specified, use 1./minium_period as the maximum frequency rather than nyquist_factor times the nyquist frequency. maximum_period : float If specified, use 1./maximum_period as the minimum frequency rather than one over the time baseline. frequency : array-like The grid of frequencies to use. If given a unit, it is converted to units of freq_unit. If not, it is assumed to be in units of freq_unit. This over rides any set frequency limits. period : array-like The grid of periods to use (as 1/period). If given a unit, it is converted to units of freq_unit. If not, it is assumed to be in units of 1/freq_unit. This overrides any set period limits. nterms : int Default 1. Number of terms to use in the Fourier fit. nyquist_factor : int Default 1. The multiple of the average Nyquist frequency. Is overriden by maximum_frequency (or minimum period). oversample_factor : int Default: None. The frequency spacing, determined by the time baseline of the lightcurve, is divided by this factor, oversampling the frequency space. This parameter is identical to the samples_per_peak parameter in astropy.LombScargle(). If normalization='amplitude', oversample_factor will be set to 5. If normalization='psd', it will be 1. These defaults can be overridden. freq_unit : `astropy.units.core.CompositeUnit` Default: None. The desired frequency units for the Lomb Scargle periodogram. This implies that 1/freq_unit is the units for period. With default normalization ('amplitude'), the freq_unit is set to 1/day, which can be overridden. 'psd' normalization will set freq_unit to microhertz. normalization : 'psd' or 'amplitude' Default: `'amplitude'`. The desired normalization of the spectrum. Can be either power spectral density (`'psd'`) or amplitude (`'amplitude'`). ls_method : str Default: `'fast'`. Passed to the `method` keyword of `astropy.stats.LombScargle()`. kwargs : dict Keyword arguments passed to `astropy.stats.LombScargle()` Returns ------- Periodogram : `Periodogram` object Returns a Periodogram object extracted from the lightcurve. """ # Input validation normalization = validate_method(normalization, ["psd", "amplitude"]) if np.isnan(lc.flux).any(): lc = lc.remove_nans() log.debug( "Lightcurve contains NaN values." "These are removed before creating the periodogram." ) # Setting default frequency units if freq_unit is None: freq_unit = 1 / u.day if normalization == "amplitude" else u.microhertz # Default oversample factor if oversample_factor is None: oversample_factor = 5.0 if normalization == "amplitude" else 1.0 if "min_period" in kwargs: warnings.warn( "`min_period` keyword is deprecated, " "please use `minimum_period` instead.", LightkurveWarning, ) minimum_period = kwargs.pop("min_period", None) if "max_period" in kwargs: warnings.warn( "`max_period` keyword is deprecated, " "please use `maximum_period` instead.", LightkurveWarning, ) maximum_period = kwargs.pop("max_period", None) if "min_frequency" in kwargs: warnings.warn( "`min_frequency` keyword is deprecated, " "please use `minimum_frequency` instead.", LightkurveWarning, ) minimum_frequency = kwargs.pop("min_frequency", None) if "max_frequency" in kwargs: warnings.warn( "`max_frequency` keyword is deprecated, " "please use `maximum_frequency` instead.", LightkurveWarning, ) maximum_frequency = kwargs.pop("max_frequency", None) # Check if any values of period have been passed and set format accordingly if not all(b is None for b in [period, minimum_period, maximum_period]): default_view = "period" else: default_view = "frequency" # If period and frequency keywords have both been set, throw an error if (not all(b is None for b in [period, minimum_period, maximum_period])) & ( not all( b is None for b in [frequency, minimum_frequency, maximum_frequency] ) ): raise ValueError( "You have input keyword arguments for both frequency and period. " "Please only use one." ) time = lc.time.copy() # Approximate Nyquist Frequency and frequency bin width in terms of days nyquist = 0.5 * (1.0 / (np.median(np.diff(time.value)))) * (1 / cds.d) fs = (1.0 / (time[-1] - time[0])) / oversample_factor # Convert these values to requested frequency unit nyquist = nyquist.to(freq_unit) fs = fs.to(freq_unit) # Warn if there is confusing input if (frequency is not None) & ( any([a is not None for a in [minimum_frequency, maximum_frequency]]) ): log.warning( "You have passed both a grid of frequencies " "and min_frequency/maximum_frequency arguments; " "the latter will be ignored." ) if (period is not None) & ( any([a is not None for a in [minimum_period, maximum_period]]) ): log.warning( "You have passed a grid of periods " "and minimum_period/maximum_period arguments; " "the latter will be ignored." ) # Tidy up the period stuff... if maximum_period is not None: # minimum_frequency MUST be none by this point. minimum_frequency = 1.0 / maximum_period if minimum_period is not None: # maximum_frequency MUST be none by this point. maximum_frequency = 1.0 / minimum_period # If the user specified a period, copy it into the frequency. if period is not None: frequency = 1.0 / period # Do unit conversions if user input min/max frequency or period if frequency is None: if minimum_frequency is not None: minimum_frequency = u.Quantity(minimum_frequency, freq_unit) if maximum_frequency is not None: maximum_frequency = u.Quantity(maximum_frequency, freq_unit) if (minimum_frequency is not None) & (maximum_frequency is not None): if minimum_frequency > maximum_frequency: if default_view == "frequency": raise ValueError( "minimum_frequency cannot be larger than maximum_frequency" ) if default_view == "period": raise ValueError( "minimum_period cannot be larger than maximum_period" ) # If nothing has been passed in, set them to the defaults if minimum_frequency is None: minimum_frequency = fs if maximum_frequency is None: maximum_frequency = nyquist * nyquist_factor # Create frequency grid evenly spaced in frequency frequency = np.arange( minimum_frequency.value, maximum_frequency.value, fs.value ) # Convert to desired units frequency = u.Quantity(frequency, freq_unit) # Change to compatible ls method if sampling not even in frequency if not implementations.main._is_regular(frequency) and ls_method in [ "fastchi2", "fast", ]: oldmethod = ls_method ls_method = {"fastchi2": "chi2", "fast": "slow"}[ls_method] log.warning( "The requested periodogram is not evenly sampled in frequency.\n" "Method has been changed from '{}' to '{}' to allow for this.".format( oldmethod, ls_method ) ) if (nterms > 1) and (ls_method not in ["fastchi2", "chi2"]): warnings.warn( "Building a Lomb Scargle Periodogram using the `slow` method. " "`nterms` has been set to >1, however this is not supported under the `{}` method. " "To run with higher nterms, set `ls_method` to either 'fastchi2', or 'chi2'. " "Please refer to the `astropy.timeseries.periodogram.LombScargle` documentation.".format( ls_method ), LightkurveWarning, ) nterms = 1 if float(astropy.__version__[0]) >= 3: LS = LombScargle( time, lc.flux, nterms=nterms, normalization="psd", **kwargs ) power = LS.power(frequency, method=ls_method) else: LS = LombScargle(time, lc.flux, nterms=nterms, **kwargs) power = LS.power(frequency, method=ls_method, normalization="psd") if normalization == "psd": # Power spectral density # Rescale from the unnormalized power output by Astropy's # Lomb-Scargle function to units of flux_variance / [frequency unit] # that may be of more interest for asteroseismology. power *= 2.0 / (len(time) * oversample_factor * fs) elif normalization == "amplitude": power = np.sqrt(power) * np.sqrt(4.0 / len(lc.time)) # Periodogram needs properties return LombScarglePeriodogram( frequency=frequency, power=power, nyquist=nyquist, targetid=lc.meta.get("TARGETID"), label=lc.meta.get("LABEL"), default_view=default_view, ls_obj=LS, nterms=nterms, ls_method=ls_method, meta=lc.meta, )
def model(self, time, frequency=None): """Obtain the flux model for a given frequency and time Parameters ---------- time : np.ndarray Time points to evaluate model. frequency : frequency to evaluate model. Default is the frequency at max power. Returns ------- result : lightkurve.LightCurve Model object with the time and flux model """ if self._LS_object is None: raise ValueError("No `astropy` Lomb Scargle object exists.") if frequency is None: frequency = self.frequency_at_max_power f = self._LS_object.model(time, frequency) lc = LightCurve( time=time, flux=f, meta={"FREQUENCY": frequency}, label="LS Model", targetid="{} LS Model".format(self.targetid), ) return lc.normalize() class BoxLeastSquaresPeriodogram(Periodogram): """Subclass of :class:`Periodogram <lightkurve.periodogram.Periodogram>` representing a power spectrum generated using the Box Least Squares (BLS) method. """ def __init__(self, *args, **kwargs): self.duration = kwargs.pop("duration", None) self.depth = kwargs.pop("depth", None) self.snr = kwargs.pop("snr", None) self._BLS_result = kwargs.pop("bls_result", None) self._BLS_object = kwargs.pop("bls_obj", None) self.transit_time = kwargs.pop("transit_time", None) self.time = kwargs.pop("time", None) self.flux = kwargs.pop("flux", None) self.time_unit = kwargs.pop("time_unit", None) super(BoxLeastSquaresPeriodogram, self).__init__(*args, **kwargs) def __repr__(self): return "BoxLeastSquaresPeriodogram(ID: {})".format(self.label)
[docs] @staticmethod def from_lightcurve(lc, **kwargs): """Creates a `Periodogram` from a LightCurve using the Box Least Squares (BLS) method. Parameters ---------- lc : `LightCurve` object The LightCurve from which to compute the Periodogram. duration : float, array_like, or `~astropy.units.Quantity`, optional The set of durations that will be considered. Default to `[0.05, 0.10, 0.15, 0.20, 0.25, 0.33]` if not specified. period : array_like or `~astropy.units.Quantity`, optional The periods where the Periodogram should be computed. If not provided, a default will be created using `BoxLeastSquares.autoperiod() <astropy.timeseries.BoxLeastSquares.autoperiod>`. minimum_period, maximum_period : float or `~astropy.units.Quantity`, optional If ``period`` is not provided, the minimum/maximum periods to search. The defaults will be computed as described in the notes below. frequency_factor : float, optional If ``period`` is not provided, a factor to control the frequency spacing of periods to be considered. kwargs : dict Keyword arguments passed to `BoxLeastSquares.power() <astropy.timeseries.BoxLeastSquares.power>` Returns ------- Periodogram : `Periodogram` object Returns a Periodogram object extracted from the lightcurve. Notes ----- If ``period`` is not provided, the default minimum period is computed from maximum duration and the median observation time gap as .. code-block:: python minimum_period = max(median(diff(lc.time)) * 4, max(duration) + median(diff(lc.time))) The default maximum period is computed as .. code-block:: python maximum_period = (max(lc.time) - min(lc.time)) / 3 ensuring that any systems with at least 3 transits are within the range of searched periods. """ # BoxLeastSquares was added to `astropy.stats` in AstroPy v3.1 and then # moved to `astropy.timeseries` in v3.2, which makes the import below # somewhat complicated. try: from astropy.timeseries import BoxLeastSquares except ImportError: try: from astropy.stats import BoxLeastSquares except ImportError: raise ImportError("BLS requires AstroPy v3.1 or later") # Validate user input for `lc` # (BoxLeastSquares will not work if flux or flux_err contain NaNs) lc = lc.remove_nans() if np.isfinite(lc.flux_err).all(): dy = lc.flux_err else: dy = None # Validate user input for `duration` duration = kwargs.pop("duration", [0.05, 0.10, 0.15, 0.20, 0.25, 0.33]) if duration is not None and ~np.all(np.isfinite(duration)): raise ValueError( "`duration` parameter contains illegal nan or inf value(s)" ) # Validate user input for `period` period = kwargs.pop("period", None) minimum_period = kwargs.pop("minimum_period", None) maximum_period = kwargs.pop("maximum_period", None) if period is not None and ~np.all(np.isfinite(period)): raise ValueError("`period` parameter contains illegal nan or inf value(s)") if minimum_period is None: if period is None: minimum_period = np.max( [ np.median(np.diff(lc.time.value)) * 4, np.max(duration) + np.median(np.diff(lc.time.value)), ] ) else: minimum_period = np.min(period) if maximum_period is None: if period is None: maximum_period = (np.max(lc.time.value) - np.min(lc.time.value)) / 3.0 else: maximum_period = np.max(period) # Validate user input for `time_unit` time_unit = kwargs.pop("time_unit", "day") if time_unit not in dir(u): raise ValueError( "{} is not a valid value for `time_unit`".format(time_unit) ) # Validate user input for `frequency_factor` frequency_factor = kwargs.pop("frequency_factor", 10) df = ( frequency_factor * np.min(duration) / (np.max(lc.time.value) - np.min(lc.time.value)) ** 2 ) npoints = int(((1 / minimum_period) - (1 / maximum_period)) / df) if npoints > 1e7: raise ValueError( "`period` contains {} points." "Periodogram is too large to evaluate. " "Consider setting `frequency_factor` to a higher value." "".format(np.round(npoints, 4)) ) elif npoints > 1e5: log.warning( "`period` contains {} points." "Periodogram is likely to be large, and slow to evaluate. " "Consider setting `frequency_factor` to a higher value." "".format(np.round(npoints, 4)) ) # Create BLS object and run the BLS search bls = BoxLeastSquares(lc.time, lc.flux, dy) if period is None: period = bls.autoperiod( duration, minimum_period=minimum_period, maximum_period=maximum_period, frequency_factor=frequency_factor, ) result = bls.power(period, duration, **kwargs) if not isinstance(result.period, u.quantity.Quantity): result.period = u.Quantity(result.period, time_unit) if not isinstance(result.power, u.quantity.Quantity): result.power = result.power * u.dimensionless_unscaled if not isinstance(result.duration, u.quantity.Quantity): result.duration = u.Quantity(result.duration, time_unit) return BoxLeastSquaresPeriodogram( frequency=1.0 / result.period, power=result.power, default_view="period", label=lc.meta.get("LABEL"), targetid=lc.meta.get("TARGETID"), transit_time=result.transit_time, duration=result.duration, depth=result.depth, bls_result=result, snr=result.depth_snr, bls_obj=bls, time=lc.time, flux=lc.flux, time_unit=time_unit, )
def compute_stats(self, period=None, duration=None, transit_time=None): """Computes commonly used vetting statistics for a transit model. See astropy.stats.bls docs for further details. Parameters ---------- period : float or Quantity Period of the transits. Default is `period_at_max_power` duration : float or Quantity Duration of the transits. Default is `duration_at_max_power` transit_time : float or Quantity Transit midpoint of the transits. Default is `transit_time_at_max_power` Returns ------- stats : dict Dictionary of vetting statistics """ if period is None: period = self.period_at_max_power log.warning("No period specified. Using period at max power") if duration is None: duration = self.duration_at_max_power log.warning("No duration specified. Using duration at max power") if transit_time is None: transit_time = self.transit_time_at_max_power log.warning("No transit time specified. Using transit time at max power") if not isinstance(transit_time, Time): transit_time = Time( transit_time, format=self.time.format, scale=self.time.scale ) return self._BLS_object.compute_stats( u.Quantity(period, "d").value, u.Quantity(duration, "d").value, transit_time ) def get_transit_model(self, period=None, duration=None, transit_time=None): """Computes the transit model using the BLS, returns a lightkurve.LightCurve See astropy.stats.bls docs for further details. Parameters ---------- period : float or Quantity Period of the transits. Default is `period_at_max_power` duration : float or Quantity Duration of the transits. Default is `duration_at_max_power` transit_time : float or Quantity Transit midpoint of the transits. Default is `transit_time_at_max_power` Returns ------- model : lightkurve.LightCurve Model of transit """ from .lightcurve import LightCurve if period is None: period = self.period_at_max_power log.warning("No period specified. Using period at max power") if duration is None: duration = self.duration_at_max_power log.warning("No duration specified. Using duration at max power") if transit_time is None: transit_time = self.transit_time_at_max_power log.warning("No transit time specified. Using transit time at max power") if not isinstance(transit_time, Time): transit_time = Time( transit_time, format=self.time.format, scale=self.time.scale ) model_flux = self._BLS_object.model( self.time, u.Quantity(period, "d").value, u.Quantity(duration, "d").value, transit_time, ) model = LightCurve(time=self.time, flux=model_flux, label="Transit Model Flux") return model def get_transit_mask(self, period=None, duration=None, transit_time=None): """Returns a boolean array that is ``True`` during transits and ``False`` elsewhere. Parameters ---------- period : float or Quantity Period of the transits. Default is `period_at_max_power` duration : float or Quantity Duration of the transits. Default is `duration_at_max_power` transit_time : float or Quantity Transit midpoint of the transits. Default is `transit_time_at_max_power` Returns ------- transit_mask : np.array of bool Mask that flags transits. Mask is ``True`` where there are transits. """ model = self.get_transit_model( period=period, duration=duration, transit_time=transit_time ) return model.flux != np.median(model.flux) @property def transit_time_at_max_power(self): """Returns the transit time corresponding to the highest peak in the periodogram.""" return self.transit_time[np.nanargmax(self.power)] @property def duration_at_max_power(self): """Returns the duration corresponding to the highest peak in the periodogram.""" return self.duration[np.nanargmax(self.power)] @property def depth_at_max_power(self): """Returns the depth corresponding to the highest peak in the periodogram.""" return self.depth[np.nanargmax(self.power)] def plot(self, **kwargs): """Plot the BoxLeastSquaresPeriodogram spectrum using matplotlib's `plot` method. See `Periodogram.plot` for details on the accepted arguments. Parameters ---------- kwargs : dict Dictionary of arguments ot be passed to `Periodogram.plot`. Returns ------- ax : `~matplotlib.axes.Axes` The matplotlib axes object. """ ax = super(BoxLeastSquaresPeriodogram, self).plot(**kwargs) if "ylabel" not in kwargs: ax.set_ylabel("BLS Power") return ax def flatten(self, **kwargs): raise NotImplementedError( "`flatten` is not implemented for `BoxLeastSquaresPeriodogram`." ) def smooth(self, **kwargs): raise NotImplementedError( "`smooth` is not implemented for `BoxLeastSquaresPeriodogram`. " )