From f013cc87a6be499300a7a1b2a195679a8a7564e1 Mon Sep 17 00:00:00 2001 From: pupperemeritus Date: Mon, 21 Aug 2023 02:53:21 +0530 Subject: [PATCH] Documentation, Algorithm Fixes, API Change --- stingray/fourier.py | 150 ++++++++++++++++++++++++---------------- stingray/lombscargle.py | 92 +++++++++++++++++------- 2 files changed, 156 insertions(+), 86 deletions(-) diff --git a/stingray/fourier.py b/stingray/fourier.py index 7931c3b58..2eca42daf 100644 --- a/stingray/fourier.py +++ b/stingray/fourier.py @@ -1932,11 +1932,12 @@ def lsft_fast( t: npt.ArrayLike, freqs: npt.ArrayLike, sign: Optional[int] = 1, - fullspec: Optional[bool] = False, oversampling: Optional[int] = 5, -): +) -> npt.NDArray: """ Calculates the Lomb-Scargle Fourier transform of a light curve. + Only considers non-negative frequencies. + Subtracts mean from data as it is required for the working of the algorithm. Parameters ---------- @@ -1963,46 +1964,42 @@ def lsft_fast( ft_res : numpy.ndarray An array of Fourier transformed data. """ - if sign not in [1, -1]: - raise ValueError("sign must be 1 or -1") - + y_ = copy.deepcopy(y) - np.mean(y) + freqs = freqs[freqs >= 0] # Constants initialization - sum_xx = np.sum(y) - num_xt = len(y) + sum_xx = np.sum(y_) + num_xt = len(y_) num_ww = len(freqs) - const1 = np.sqrt(0.5) * np.sqrt(num_ww) - const2 = const1 # Arrays initialization ft_real = ft_imag = np.zeros((num_ww)) f0, df, N = freqs[0], freqs[1] - freqs[0], len(freqs) # Sum (y_i * cos(wt - wtau)) - Sh, Ch = trig_sum(t, y, df, N, f0, oversampling=oversampling) + Sh, Ch = trig_sum(t, y_, df, N, f0, oversampling=oversampling) - # Summation of (cos(wt - wtau))^2 and (sin(wt - wtau))^2 - _, scos2x = trig_sum( - t, - np.ones_like(y) / len(y), - df, - N, - f0, - freq_factor=2, - oversampling=oversampling, - ) - # cos^2(x) = (1 + cos(2x))/2 - # sin^2(x) = (1 - cos(2x))/2 - C2, S2 = (1 + scos2x) / 2, (1 - scos2x) / 2 # Angular Frequency w = freqs * 2 * np.pi # Summation of cos(2wt) and sin(2wt) csum, ssum = trig_sum( - t, np.ones_like(y) / len(y), df, N, f0, freq_factor=2, oversampling=oversampling + t, np.ones_like(y_) / len(y_), df, N, f0, freq_factor=2, oversampling=oversampling ) wtau = 0.5 * np.arctan2(ssum, csum) + S2, C2 = trig_sum( + t, + np.ones_like(y_), + df, + N, + f0, + freq_factor=2, + oversampling=oversampling, + ) + + const = np.sqrt(0.5) * np.sqrt(num_ww) + coswtau = np.cos(wtau) sinwtau = np.sin(wtau) @@ -2015,21 +2012,17 @@ def lsft_fast( scos2 = 0.5 * (N + C2 * cos2wtau + S2 * sin2wtau) ssin2 = 0.5 * (N - C2 * cos2wtau - S2 * sin2wtau) - # Calculating the real and imaginary components of the Fourier transform - ft_real = const1 * sumr / np.sqrt(scos2) - ft_imag = const2 * sumi / np.sqrt(ssin2) + ft_real = const * sumr / np.sqrt(scos2) + ft_imag = const * np.sign(sign) * sumi / np.sqrt(ssin2) - # Looping through all the frequencies - ft_real[0] = sum_xx / np.sqrt(num_xt) - ft_imag[0] = 0 + ft_real[freqs == 0] = sum_xx / np.sqrt(num_xt) + ft_imag[freqs == 0] = 0 - # Resultant fourier transform for current angular frequency - ft_res = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * w * t[0]) + phase = wtau - w * t[0] - if sign == -1: - ft_res = np.conjugate(ft_res) + ft_res = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * phase) - return ft_res if fullspec else ft_res[freqs > 0] + return ft_res def lsft_slow( @@ -2037,10 +2030,11 @@ def lsft_slow( t: npt.ArrayLike, freqs: npt.ArrayLike, sign: Optional[int] = 1, - fullspec: Optional[bool] = False, -): +) -> npt.NDArray: """ Calculates the Lomb-Scargle Fourier transform of a light curve. + Only considers non-negative frequencies. + Subtracts mean from data as it is required for the working of the algorithm. Parameters ---------- @@ -2064,22 +2058,27 @@ def lsft_slow( ft_res : numpy.ndarray An array of Fourier transformed data. """ - if sign not in (1, -1): - raise ValueError("sign must be 1 or -1") - - # Constants initialization - const1 = np.sqrt(0.5) * np.sqrt(len(y)) - const2 = const1 + y_ = copy.deepcopy(y) - np.mean(y) + freqs = freqs[freqs >= 0] ft_real = np.zeros_like(freqs) ft_imag = np.zeros_like(freqs) ft_res = np.zeros_like(freqs, dtype=np.complex128) - for i, freq in enumerate(freqs): - wrun = freq * 2 * np.pi - if i == 0: - ft_real[i] = sum(y) / np.sqrt(len(y)) - ft_imag[i] = 0 + num_y = len(y_) + num_freqs = len(freqs) + sum_y = np.sum(y_) + const1 = np.sqrt(0.5) * np.sqrt(num_y) + const2 = const1 * np.sign(sign) + ft_real = ft_imag = np.zeros(num_freqs) + ft_res = np.zeros(num_freqs, dtype=np.complex128) + + for i in range(num_freqs): + wrun = freqs[i] * 2 * np.pi + if wrun == 0: + ft_real = sum_y / np.sqrt(num_y) + ft_imag = 0 + phase_this = 0 else: csum = np.sum(np.cos(2.0 * wrun * t)) ssum = np.sum(np.sin(2.0 * wrun * t)) @@ -2087,20 +2086,51 @@ def lsft_slow( watan = np.arctan2(ssum, csum) wtau = 0.5 * watan - cos_wt_wtau = np.cos(wrun * t - wtau) - sin_wt_wtau = np.sin(wrun * t - wtau) + sumr = np.sum(y_ * np.cos(wrun * t - wtau)) + sumi = np.sum(y_ * np.sin(wrun * t - wtau)) - sumr = np.sum(y * cos_wt_wtau) - sumi = np.sum(y * sin_wt_wtau) + scos2 = np.sum(np.power(np.cos(wrun * t - wtau), 2)) + ssin2 = np.sum(np.power(np.sin(wrun * t - wtau), 2)) - scos2 = np.sum(np.power(cos_wt_wtau, 2)) - ssin2 = np.sum(np.power(sin_wt_wtau, 2)) + ft_real = const1 * sumr / np.sqrt(scos2) + ft_imag = const2 * sumi / np.sqrt(ssin2) + phase_this = wtau - wrun * t[0] + ft_res[i] = np.complex128(ft_real + (1j * ft_imag)) * np.exp(-1j * phase_this) + return ft_res - ft_real[i] = const1 * sumr / np.sqrt(scos2) - ft_imag[i] = const2 * sumi / np.sqrt(ssin2) - ft_res[i] = np.complex128(ft_real[i] + (1j * ft_imag[i])) * np.exp(-1j * wrun * t[0]) - if sign == -1: - ft_res = np.conjugate(ft_res) +def impose_symmetry_lsft( + ft_res: npt.ArrayLike, + sum_y: float, + len_y: int, + freqs: npt.ArrayLike, +) -> npt.ArrayLike: + """ + Impose symmetry on the input fourier transform. - return ft_res if fullspec else ft_res[freqs > 0] + Parameters + ---------- + ft_res : np.array + The Fourier transform of the signal. + sum_y : float + The sum of the values of the signal. + len_y : int + The length of the signal. + freqs : np.array + An array of frequencies at which the transform is sampled. + + Returns + ------- + lsft_res : np.array + The Fourier transform of the signal with symmetry imposed. + freqs_new : np.array + The new frequencies + """ + zero_present = np.amy(freqs == 0) + if not zero_present: + ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res)), 0, ft_res]) + freqs_new = np.concatenate([-np.flip(freqs), 0, freqs]) + else: + ft_res_new = np.concatenate([np.conjugate(np.flip(ft_res))[1:], ft_res]) + freqs_new = np.concatenate([-np.flip(freqs)[1:], freqs]) + return ft_res_new, freqs_new diff --git a/stingray/lombscargle.py b/stingray/lombscargle.py index 60d6e4498..e285e190a 100644 --- a/stingray/lombscargle.py +++ b/stingray/lombscargle.py @@ -8,7 +8,7 @@ from .crossspectrum import Crossspectrum from .events import EventList from .exceptions import StingrayError -from .fourier import lsft_fast, lsft_slow +from .fourier import lsft_fast, lsft_slow, impose_symmetry_lsft from .lightcurve import Lightcurve from .utils import simon @@ -33,7 +33,7 @@ class LombScargleCrossspectrum(Crossspectrum): norm: {``frac``, ``abs``, ``leahy``, ``none``}, string, optional, default ``none`` The normalization of the cross spectrum. - power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``real`` + power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``all`` Parameter to choose among complete, real part and magnitude of the cross spectrum. fullspec: boolean, optional, default ``False`` @@ -115,7 +115,7 @@ def __init__( data1: Optional[Union[EventList, Lightcurve]] = None, data2: Optional[Union[EventList, Lightcurve]] = None, norm: Optional[str] = "none", - power_type: Optional[str] = "real", + power_type: Optional[str] = "all", dt: Optional[float] = None, fullspec: Optional[bool] = False, skip_checks: bool = False, @@ -453,7 +453,9 @@ def _ls_cross(self, lc1, lc2, freq=None, fullspec=False, method="fast", oversamp fit_mean=False, center_data=False, normalization="psd", - ).autofrequency(minimum_frequency=self.min_freq, maximum_frequency=self.max_freq), + ).autofrequency( + minimum_frequency=max(self.min_freq, 0), maximum_frequency=self.max_freq + ), )[0] freqs2 = ( LombScargle( @@ -462,34 +464,43 @@ def _ls_cross(self, lc1, lc2, freq=None, fullspec=False, method="fast", oversamp fit_mean=False, center_data=False, normalization="psd", - ).autofrequency(minimum_frequency=self.min_freq, maximum_frequency=self.max_freq), + ).autofrequency( + minimum_frequency=max(self.min_freq, 0), maximum_frequency=self.max_freq + ), )[0] if max(freqs2) > max(freq): freq = freqs2 if method == "slow": - lsft1 = lsft_slow(lc1.counts, lc1.time, freq, sign=-1, fullspec=fullspec) - lsft2 = lsft_slow(lc2.counts, lc2.time, freq, sign=1, fullspec=fullspec) + lsft1 = lsft_slow(lc1.counts, lc1.time, freq) + lsft2 = lsft_slow(lc2.counts, lc2.time, freq) elif method == "fast": lsft1 = lsft_fast( lc1.counts, lc1.time, freq, - sign=-1, - fullspec=fullspec, oversampling=oversampling, ) lsft2 = lsft_fast( lc2.counts, lc2.time, freq, - fullspec=fullspec, oversampling=oversampling, ) - cross = np.multiply(lsft1, lsft2) - if not fullspec: - freq = freq[freq > 0] - cross = cross[freq > 0] + if fullspec: + lsft1, _ = impose_symmetry_lsft( + lsft1, + np.sum((lc1.counts)), + lc1.n, + freq, + ) + lsft2, freq = impose_symmetry_lsft( + lsft2, + np.sum(lc2.counts), + lc2.n, + freq, + ) + cross = np.multiply(lsft1, np.conjugate(lsft2)) return freq, cross def _initialize_empty(self): @@ -515,22 +526,51 @@ def _initialize_empty(self): self.variance2 = None return + def phase_lag(self): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'phase_lag' ! Not applicable for unevenly sampled data" + ) + 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 + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'time_lag' ! Not applicable for unevenly sampled data" + ) - ..math:: + def classical_significances(self, threshold=1, trial_correction=False): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'classical_significances' ! Not applicable for unevenly sampled data" + ) - \tau = \frac{\phi}{\two pi \nu} + def from_time_array(): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'from_time_array' ! Not applicable for unevenly sampled data" + ) - where :math:`\nu` is the center of the frequency bins. - """ - if self.__class__ == LombScargleCrossspectrum: - ph_lag = self.phase_lag() + def from_events(): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'from_events' ! Not applicable for unevenly sampled data" + ) - return ph_lag / (2 * np.pi * self.freq) - else: - raise AttributeError("Object has no attribute named 'time_lag' !") + def from_lightcurve(): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'from_lightcurve' ! Not applicable for unevenly sampled data" + ) + + def from_lc_iterable(): + """Not applicable for unevenly sampled data""" + raise AttributeError( + "Object has no attribute named 'from_lc_iterable' ! Not applicable for unevenly sampled data" + ) + + def _initialize_from_any_input(): + """Not required for unevenly sampled data""" + raise AttributeError("Object has no attribute named '_initialize_from_any_input' !") class LombScarglePowerspectrum(LombScargleCrossspectrum): @@ -552,7 +592,7 @@ class LombScarglePowerspectrum(LombScargleCrossspectrum): norm: {``frac``, ``abs``, ``leahy``, ``none``}, string, optional, default ``none`` The normalization of the cross spectrum. - power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``real`` + power_type: {``real``, ``absolute``, ``all`}, string, optional, default ``all`` Parameter to choose among complete, real part and magnitude of the cross spectrum. fullspec: boolean, optional, default ``False``