Fits event lists

CHANGELOG.rst

@@ -1,3 +1,31 @@
+v2.2 (2024-10-22)
+-----------------
+
+New Features
+^^^^^^^^^^^^
+
+- Add a compute_rms function to LombScarglePowerspectrum (`#828 <https://github.com/StingraySoftware/stingray/pull/828>`__)
+- Introduced FITSReader class for lazy-loading of event lists (`#834 <https://github.com/StingraySoftware/stingray/pull/834>`__)
+- implementation of the shift-and-add technique for QPOs and other varying power spectral features (`#849 <https://github.com/StingraySoftware/stingray/pull/849>`__)
+
+
+Bug Fixes
+^^^^^^^^^
+
+- The ``fold_events`` function now checks if the keyword arguments (`kwargs`) are in the list of optional parameters.
+  If any unidentified keys are present, it raises a `ValueError`.
+  This fix ensures that the function only accepts valid optional parameters and provides a clear error message for unsupported keys. (`#837 <https://github.com/StingraySoftware/stingray/pull/837>`__)
+
+
+Internal Changes
+^^^^^^^^^^^^^^^^
+
+- Eliminated runtime dependency on setuptools (`#852 <https://github.com/StingraySoftware/stingray/pull/852>`__)
+- Moved configuration to pyproject.toml as recommended by PEP 621 (`#842 <https://github.com/StingraySoftware/stingray/pull/842>`__)
+- Added pre-commit hooks in ``pre-commit-config.yaml`` (`#847 <https://github.com/StingraySoftware/stingray/pull/847>`__)
+- Improved main page of the documentation (`#748 <https://github.com/StingraySoftware/stingray/pull/748>`__)
+
+
 v2.1 (2024-05-29)
 -----------------
 

README.rst

@@ -7,7 +7,7 @@ Stingray
    :widths: 50, 50, 50, 50
 
    |Python version|, |GitHub release|, |Build Status Master|, |Slack|
-   |Docs|, |joss|, |Repo status|,  " "
+   |Docs|, |joss|, |Repo status|,  |pyOpenSci Peer-Reviewed|
    |License|, |doi|, |Coverage Status Master|, " "
 
 
@@ -98,6 +98,8 @@ The code is distributed under the MIT license; see `LICENSE.rst <LICENSE.rst>`_
    :target: https://www.repostatus.org/#active
 .. |License| image:: https://img.shields.io/badge/License-MIT-yellow.svg
    :target: https://opensource.org/licenses/MIT
+.. |pyOpenSci Peer-Reviewed| image:: https://pyopensci.org/badges/peer-reviewed.svg 
+   :target: https://github.com/pyOpenSci/software-review/issues/201
 
 .. _Astropy: https://www.github.com/astropy/astropy
 .. _Issues: https://www.github.com/stingraysoftware/stingray/issues

docs/citing.rst

@@ -12,7 +12,7 @@ If possible, we ask that you cite a DOI corresponding to the specific version of
 
 .. include:: _zenodo.rst
 
-If this isn't possible — for example, because you worked with an unreleased version of the code — you can cite Stingray's `concept DOI <https://zenodo.org/help/versioning>`__, `10.5281/zenodo.1490116 <https://zenodo.org/record/1490116>`__ (`BibTeX <https://zenodo.org/record/1490116/export/hx>`__), which will always resolve to the latest release.
+If this isn't possible — for example, because you worked with an unreleased version of the code — you can cite Stingray's `concept DOI <https://zenodo.org/help/versioning>`__, `10.5281/zenodo.1490116 <https://zenodo.org/records/1490116>`__ (`BibTeX <https://zenodo.org/records/1490116/export/hx>`__), which will always resolve to the latest release.
 
 Papers
 ======

docs/conf.py

@@ -158,8 +158,13 @@
 # (source start file, name, description, authors, manual section).
 man_pages = [("index", project.lower(), project + " Documentation", [author], 1)]
 
-# Trust the links from doi.org, even if they might have Client errors or other minor issues
-linkcheck_ignore = [r"https://doi.org/"]
+# Trust the links from these sites, even if they might have Client errors or other minor issues
+linkcheck_ignore = [
+    r"https://doi.org/",
+    r"https://arxiv.org/",
+    r"https://.*adsabs.harvard.edu/",
+    r"https://zenodo.org/",
+]
 
 # -- Options for the edit_on_github extension ---------------------------------
 
@@ -197,7 +202,7 @@ class Release(object):
 
     @property
     def zenodo_url(self):
-        return f"https://zenodo.org/record/{self.doi.split('.')[-1]}"
+        return f"https://zenodo.org/records/{self.doi.split('.')[-1]}"
 
     @property
     def github_url(self):

docs/index.rst

@@ -308,6 +308,7 @@ Advanced
    :maxdepth: 2
 
    timeseries
+   largedata
    api
 
 Additional information

docs/largedata.rst

@@ -0,0 +1,7 @@
+Working with large data sets
+****************************
+
+.. toctree::
+   :maxdepth: 2
+
+   notebooks/Performance/Dealing with large data files.ipynb

docs/timeseries.rst

@@ -6,3 +6,4 @@ Working with more generic time series
 
    notebooks/StingrayTimeseries/StingrayTimeseries Tutorial.ipynb
    notebooks/StingrayTimeseries/Working with weights and polarization.ipynb
+   notebooks/Modeling/GP_Modeling/GP_modeling_tutorial.ipynb

pyproject.toml

@@ -34,8 +34,6 @@ keywords = [
 dependencies = [
     "numpy>=1.17",
     "astropy>=4.0",
-    "setuptools",
-    "setuptools_scm",
     "scipy>=1.1.0",
     "matplotlib>=3.0,!=3.4.00",
 ]
@@ -152,7 +150,7 @@ filterwarnings = [
     "ignore:.*Number of segments used in averaging.*:UserWarning",
     "ignore:.*:astropy.units.core.UnitsWarning",
     "ignore:.*cannot be added to FITS Header:astropy.utils.exceptions.AstropyUserWarning",
-    "ignore:Numba not installed:UserWarning",
+    "ignore:The recommended numba package is not installed:UserWarning",
     "ignore:More than 20 figures have been opened.:RuntimeWarning",
     "ignore:This platform does not support:RuntimeWarning",
     "ignore:Some error bars in the Averaged:UserWarning",

stingray/crossspectrum.py

@@ -2132,7 +2132,22 @@ class DynamicalCrossspectrum(AveragedCrossspectrum):
         The number of averaged powers in each spectral bin (initially 1, it changes after rebinning).
     """
 
-    def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_time=None):
+    def __init__(
+        self, data1=None, data2=None, segment_size=None, norm="frac", gti=None, sample_time=None
+    ):
+        self.segment_size = segment_size
+        self.sample_time = sample_time
+        self.gti = gti
+        self.norm = norm
+
+        if segment_size is None and data1 is None and data2 is None:
+            self._initialize_empty()
+            self.dyn_ps = None
+            return
+
+        if segment_size is None or data1 is None or data2 is None:
+            raise TypeError("data1, data2, and segment_size must all be specified")
+
         if isinstance(data1, EventList) and sample_time is None:
             raise ValueError("To pass input event lists, please specify sample_time")
         elif isinstance(data1, Lightcurve):
@@ -2145,11 +2160,6 @@ def __init__(self, data1, data2, segment_size, norm="frac", gti=None, sample_tim
         if segment_size < 2 * sample_time:
             raise ValueError("Length of the segment is too short to form a light curve!")
 
-        self.segment_size = segment_size
-        self.sample_time = sample_time
-        self.gti = gti
-        self.norm = norm
-
         self._make_matrix(data1, data2)
 
     def _make_matrix(self, data1, data2):
@@ -2370,6 +2380,74 @@ def trace_maximum(self, min_freq=None, max_freq=None):
 
         return np.array(max_positions)
 
+    def shift_and_add(self, f0_list, nbins=100, output_obj_type=AveragedCrossspectrum, rebin=None):
+        """Shift and add the dynamical cross spectrum.
+
+        This is the basic operation for the shift-and-add operation used to track
+        kHz QPOs in X-ray binaries (e.g. Méndez et al. 1998, ApJ, 494, 65).
+
+        Parameters
+        ----------
+        freqs : np.array
+            Array of frequencies, the same for all powers. Must be sorted and on a uniform
+            grid.
+        power_list : list of np.array
+            List of power spectra. Each power spectrum must have the same length
+            as the frequency array.
+        f0_list : list of float
+            List of central frequencies
+
+        Other parameters
+        ----------------
+        nbins : int, default 100
+            Number of bins to extract
+        rebin : int, default None
+            Rebin the final spectrum by this factor. At the moment, the rebinning
+            is linear.
+        output_obj_type : class, default :class:`AveragedCrossspectrum`
+            The type of the output object. Can be, e.g. :class:`AveragedCrossspectrum` or
+            :class:`AveragedPowerspectrum`.
+
+        Returns
+        -------
+        output: :class:`AveragedPowerspectrum` or :class:`AveragedCrossspectrum`
+            The final averaged power spectrum.
+
+        Examples
+        --------
+        >>> power_list = [[2, 5, 2, 2, 2], [1, 1, 5, 1, 1], [3, 3, 3, 5, 3]]
+        >>> power_list = np.array(power_list).T
+        >>> freqs = np.arange(5) * 0.1
+        >>> f0_list = [0.1, 0.2, 0.3, 0.4]
+        >>> dps = DynamicalCrossspectrum()
+        >>> dps.dyn_ps = power_list
+        >>> dps.freq = freqs
+        >>> dps.df = 0.1
+        >>> dps.m = 1
+        >>> output = dps.shift_and_add(f0_list, nbins=5)
+        >>> assert np.array_equal(output.m, [2, 3, 3, 3, 2])
+        >>> assert np.array_equal(output.power, [2. , 2. , 5. , 2. , 1.5])
+        >>> assert np.allclose(output.freq, [0.05, 0.15, 0.25, 0.35, 0.45])
+        """
+        from .fourier import shift_and_add
+
+        final_freqs, final_powers, count = shift_and_add(
+            self.freq, self.dyn_ps.T, f0_list, nbins=nbins, M=self.m, df=self.df, rebin=rebin
+        )
+
+        output = output_obj_type()
+        good = ~np.isnan(final_powers)
+
+        output.freq = final_freqs[good]
+        output.power = final_powers[good]
+        output.m = count[good]
+        output.df = self.df
+        output.norm = self.norm
+        output.gti = self.gti
+        output.segment_size = self.segment_size
+
+        return output
+
     def power_colors(
         self,
         freq_edges=[1 / 256, 1 / 32, 0.25, 2.0, 16.0],

stingray/events.py

@@ -14,7 +14,8 @@
 from .base import StingrayTimeseries
 from .filters import get_deadtime_mask
 from .gti import generate_indices_of_boundaries
-from .io import load_events_and_gtis, pi_to_energy
+from .io import pi_to_energy, get_file_extension
+from .io import FITSTimeseriesReader
 from .lightcurve import Lightcurve
 from .utils import simon, njit
 from .utils import histogram
@@ -620,28 +621,18 @@ def read(cls, filename, fmt=None, rmf_file=None, **kwargs):
         ev: :class:`EventList` object
             The :class:`EventList` object reconstructed from file
         """
-
+        if fmt is None:
+            for fits_ext in ["fits", "evt"]:
+                if fits_ext in get_file_extension(filename).lower():
+                    fmt = "hea"
+                    break
         if fmt is not None and fmt.lower() in ("hea", "ogip"):
-            evtdata = load_events_and_gtis(filename, **kwargs)
-
-            evt = EventList(
-                time=evtdata.ev_list,
-                gti=evtdata.gti_list,
-                pi=evtdata.pi_list,
-                energy=evtdata.energy_list,
-                mjdref=evtdata.mjdref,
-                instr=evtdata.instr,
-                mission=evtdata.mission,
-                header=evtdata.header,
-                detector_id=evtdata.detector_id,
-                ephem=evtdata.ephem,
-                timeref=evtdata.timeref,
-                timesys=evtdata.timesys,
-            )
-            if "additional_columns" in kwargs:
-                for key in evtdata.additional_data:
-                    if not hasattr(evt, key.lower()):
-                        setattr(evt, key.lower(), evtdata.additional_data[key])
+            additional_columns = kwargs.pop("additional_columns", None)
+
+            evt = FITSTimeseriesReader(
+                filename, output_class=EventList, additional_columns=additional_columns
+            )[:]
+
             if rmf_file is not None:
                 evt.convert_pi_to_energy(rmf_file)
             return evt