diff --git a/stingray/lombscargle.py b/stingray/lombscargle.py index b3ed92d91..f706a066f 100644 --- a/stingray/lombscargle.py +++ b/stingray/lombscargle.py @@ -13,7 +13,6 @@ from .utils import simon -# method default will change to fast after implementation of the fast algorithm class LombScargleCrossspectrum(Crossspectrum): main_array_attr = "freq" type = "crossspectrum" @@ -64,6 +63,9 @@ class LombScargleCrossspectrum(Crossspectrum): The method to be used by the Lomb-Scargle Fourier Transformation function. `fast` and `slow` are the allowed values. Default is `fast`. fast uses the optimized Press and Rybicki O(n*log(n)) + + oversampling : float, optional, default: 5 + Interpolation Oversampling Factor (for the fast algorithm) Attributes ---------- @@ -113,7 +115,8 @@ def __init__( min_freq: float = 0, max_freq: float = None, df: float = None, - method="slow", + method: str = "fast", + oversampling: int = 5, ): self._type = None good_input = data1 is not None and data2 is not None @@ -128,8 +131,15 @@ def __init__( min_freq=min_freq, max_freq=max_freq, df=df, + method=method, + oversampling=oversampling, ) + if dt is None: + if isinstance(data1, Lightcurve): + dt = data1.dt + elif isinstance(data1, EventList): + raise ValueError("dt must be provided for EventLists") self.dt = dt norm = norm.lower() self.norm = norm @@ -138,12 +148,6 @@ def __init__( if not good_input: return self._initialize_empty() - if type(data1) not in [EventList, Lightcurve] or type(data2) not in [ - EventList, - Lightcurve, - ]: - raise TypeError("One of the arguments is not of type eventlist or lightcurve") - if isinstance(data1, EventList): self.lc1 = data1.to_lc(self.dt) else: @@ -161,7 +165,9 @@ def __init__( self.min_freq = min_freq self.max_freq = max_freq - self._make_crossspectrum(self.lc1, self.lc2, fullspec, method) + self._make_crossspectrum( + self.lc1, self.lc2, fullspec, method=method, oversampling=oversampling + ) if self.power_type == "abs": self.power = np.abs(self.power) self.power_err = np.abs(self.power_err) @@ -174,7 +180,6 @@ def __init__( self.unnorm_power_err = np.real(self.unnorm_power_err) self._make_auxil_pds(self.lc1, self.lc2) - # method default will change to fast after implementation of the fast algorithm def initial_checks( self, data1=None, @@ -186,7 +191,8 @@ def initial_checks( max_freq=None, fullspec=False, df=None, - method="slow", + method="fast", + oversampling=5, ): if not isinstance(norm, str): raise TypeError("norm must be a string") @@ -207,10 +213,28 @@ def initial_checks( if max_freq is not None: if max_freq < min_freq or max_freq < 0: raise ValueError("max_freq must be non-negative and greater than min_freq") + + if method not in ["fast", "slow"]: + raise ValueError("method must be one of ['fast','slow']") + + if not isinstance(oversampling, int): + raise TypeError("oversampling must be an integer") + + if not isinstance(df, float) and df is not None: + raise TypeError("df must be a float") + + if not isinstance(fullspec, bool): + raise TypeError("fullspec must be a boolean") + + if type(data1) not in [EventList, Lightcurve] or type(data2) not in [ + EventList, + Lightcurve, + ]: + raise TypeError("One of the arguments is not of type eventlist or lightcurve") + return True - # method default will change to fast after implementation of the fast algorithm - def _make_crossspectrum(self, lc1, lc2, fullspec=False, method="fast"): + def _make_crossspectrum(self, lc1, lc2, fullspec=False, method="fast", oversampling=5): """ Auxiliary method computing the normalized cross spectrum from two light curves. This includes checking for the presence of and @@ -272,7 +296,13 @@ def _make_crossspectrum(self, lc1, lc2, fullspec=False, method="fast"): self.m = 1 - self.freq, self.unnorm_power = self._ls_cross(self.lc1, self.lc2, fullspec=fullspec) + self.freq, self.unnorm_power = self._ls_cross( + self.lc1, + self.lc2, + fullspec=fullspec, + method=method, + oversampling=oversampling, + ) self.power = self._normalize_crossspectrum(self.unnorm_power) if lc1.err_dist.lower() != lc2.err_dist.lower(): @@ -327,8 +357,7 @@ def _make_auxil_pds(self, lc1, lc2): df=self.df, ) - # method default will change to fast after implementation of the fast algorithm - def _ls_cross(self, lc1, lc2, ww=None, fullspec=False, method="slow"): + def _ls_cross(self, lc1, lc2, ww=None, fullspec=False, method="fast", oversampling=5): """ Lomb-Scargle Fourier transform the two light curves, then compute the cross spectrum. Computed as CS = lc1 x lc2* (where lc2 is the one that gets @@ -387,8 +416,18 @@ def _ls_cross(self, lc1, lc2, ww=None, fullspec=False, method="slow"): if method == "slow": lsft1 = lsft_slow(lc1.counts, lc1.time, ww, sign=1, fullspec=fullspec) lsft2 = lsft_slow(lc2.counts, lc2.time, ww, sign=-1, fullspec=fullspec) - # elif method == "fast": - # lsft1 = lsft_fast(lc1, ww, fullspec=fullspec) + elif method == "fast": + lsft1 = lsft_fast( + lc1.counts, lc1.time, ww, fullspec=fullspec, oversampling=oversampling + ) + lsft2 = lsft_fast( + lc2.counts, + lc2.time, + ww, + fullspec=fullspec, + sign=-1, + oversampling=oversampling, + ) cross = np.multiply(lsft1, lsft2) freq = ww if not fullspec: @@ -412,6 +451,23 @@ def _initialize_empty(self): self.k = 1 return + def time_lag(self): + r"""Calculate the fourier time lag of the cross spectrum. + The time lag is calculated by taking the phase lag :math:`\phi` and + + ..math:: + + \tau = \frac{\phi}{\two pi \nu} + + where :math:`\nu` is the center of the frequency bins. + """ + if self.__class__ == LombScargleCrossspectrum: + ph_lag = self.phase_lag() + + return ph_lag / (2 * np.pi * self.freq) + else: + raise AttributeError("Object has no attribute named 'time_lag' !") + class LombScarglePowerspectrum(LombScargleCrossspectrum): type = "powerspectrum" @@ -458,6 +514,9 @@ class LombScarglePowerspectrum(LombScargleCrossspectrum): and `slow` are the allowed values. Default is `fast`. fast uses the optimized Press and Rybicki O(n*log(n)) + oversampling : float, optional, default: 5 + Interpolation Oversampling Factor (for the fast algorithm) + Attributes ---------- norm: {"leahy" | "frac" | "abs" | "none" } @@ -503,6 +562,7 @@ def __init__( max_freq: Optional[float] = None, df: Optional[float] = None, method: Optional[str] = "fast", + oversampling: Optional[int] = 5, ): self._type = None data1 = copy.deepcopy(data) @@ -520,6 +580,7 @@ def __init__( max_freq=max_freq, df=df, method=method, + oversampling=oversampling, ) if type(data) not in [EventList, Lightcurve, None]: good_input = False @@ -543,6 +604,7 @@ def __init__( max_freq=max_freq, df=df, method=method, + oversampling=oversampling, ) self.nphots = self.nphots1