Documentation, Algorithm Fixes, API Change

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,32 +2012,29 @@ 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(
     y: npt.ArrayLike,
     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,43 +2058,79 @@ 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))
 
             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

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``