Repository: StingraySoftware/stingray (python)

Issue:
Bad tail of covariance spectra at small count rates.
### Description of the Bug

At low count rates, covariance spectra can appear like this:

<img width="562" alt="Screenshot 2025-01-30 at 22 11 43" src="https://github.com/user-attachments/assets/df572565-aff6-45fa-ad8c-8eaed811fdae" />

It turns out that this is due to a bad estimate of the Poisson level (that gets subtracted to calculate the rms!) when the counts per segment are few. In this case, the Poisson level is underestimated, and the rms appears higher. 

### Steps/Code to Replicate the Bug

TBD

### Expected Results

TBD

### Actual Results

TBD

Pull request:
Covariance spectrum fixes
Resolves #873 

Patch:
diff --git a/docs/changes/874.bugfix.rst b/docs/changes/874.bugfix.rst
new file mode 100644
index 00000000..5e7e2387
--- /dev/null
+++ b/docs/changes/874.bugfix.rst
@@ -0,0 +1,1 @@
+Set a lower limit to the number of photons in a segment of data used for certain ``VarEnergySpectrum`` subclasses. This avoids, e.g., spurious high covariance measurements in low-count data sets.
\ No newline at end of file
diff --git a/stingray/base.py b/stingray/base.py
index 57f99355..ae311077 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -319,15 +319,17 @@ class StingrayObject(object):
                 if not getattr(self, attr, None) == getattr(other_ts, attr, None):
                     return False
             else:
-                if not np.array_equal(getattr(self, attr, None), getattr(other_ts, attr, None)):
+                if not np.array_equal(
+                    getattr(self, attr, None), getattr(other_ts, attr, None), equal_nan=True
+                ):
                     return False
 
         for attr in self.array_attrs():
-            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr)):
+            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr), equal_nan=True):
                 return False
 
         for attr in self.internal_array_attrs():
-            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr)):
+            if not np.array_equal(getattr(self, attr), getattr(other_ts, attr), equal_nan=True):
                 return False
 
         return True
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index 4369d80f..4a2fb6f2 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -1,3 +1,4 @@
+from collections.abc import Iterable
 import numpy as np
 import warnings
 from stingray.base import StingrayObject
@@ -188,30 +189,50 @@ class VarEnergySpectrum(StingrayObject, metaclass=ABCMeta):
 
     def __init__(
         self,
-        events,
-        freq_interval,
-        energy_spec,
+        events=None,
+        freq_interval=None,
+        energy_spec=None,
         ref_band=None,
         bin_time=1,
         use_pi=False,
         segment_size=None,
         events2=None,
         return_complex=False,
+        min_phot_per_segment=10,
     ):
-        self.events1 = events
-        self.events2 = assign_value_if_none(events2, events)
-        self._analyze_inputs()
+        if isinstance(energy_spec, tuple):
+            energies = _decode_energy_specification(energy_spec)
+        elif isinstance(energy_spec, Iterable) and not isinstance(energy_spec, str):
+            energies = np.asanyarray(energy_spec)
+        elif events is not None:
+            raise ValueError(
+                f"Energy specification must be a tuple or a list (input: {energy_spec})"
+            )
+
         # This will be set to True in ComplexCovariance
         self.return_complex = return_complex
-
         self.freq_interval = freq_interval
         self.use_pi = use_pi
         self.bin_time = bin_time
-
-        if isinstance(energy_spec, tuple):
-            energies = _decode_energy_specification(energy_spec)
+        self.min_phot_per_segment = min_phot_per_segment
+        if energy_spec is not None:
+            self.energy_intervals = list(zip(energies[0:-1], energies[1:]))
         else:
-            energies = np.asanyarray(energy_spec)
+            self.energy_intervals = None
+        self.ref_band = np.asanyarray(assign_value_if_none(ref_band, [0, np.inf]))
+        if len(self.ref_band.shape) <= 1:
+            self.ref_band = np.asanyarray([self.ref_band])
+
+        self.segment_size = self.delta_nu = None
+
+        if events is None:
+            self.events1 = self.events2 = None
+            # This will be set to True in ComplexCovariance
+            return
+
+        self.events1 = events
+        self.events2 = assign_value_if_none(events2, events)
+        self._analyze_inputs()
 
         self.energy_intervals = list(zip(energies[0:-1], energies[1:]))
 
@@ -235,7 +256,16 @@ class VarEnergySpectrum(StingrayObject, metaclass=ABCMeta):
     @property
     def energy(self):
         """Give the centers of the energy intervals."""
-        return np.sum(self.energy_intervals, axis=1) / 2
+        if not hasattr(self, "_energy") or self._energy is None:
+            if self.energy_intervals is None or len(self.energy_intervals) == 0:
+                self._energy = []
+            else:
+                self._energy = np.sum(self.energy_intervals, axis=1) / 2
+        return self._energy
+
+    @energy.setter
+    def energy(self, value):
+        self._energy = value
 
     def _analyze_inputs(self):
         """Make some checks on the inputs and set some internal variable.
@@ -414,14 +444,27 @@ class VarEnergySpectrum(StingrayObject, metaclass=ABCMeta):
     def _spectrum_function(self):
         pass
 
-    def from_astropy_table(self, *args, **kwargs):
-        raise NotImplementedError("from_XXXX methods are not implemented for VarEnergySpectrum")
-
-    def from_xarray(self, *args, **kwargs):
-        raise NotImplementedError("from_XXXX methods are not implemented for VarEnergySpectrum")
+    def check_phot_per_segment(self, phot_per_segment, eint):
+        """Check if the number of photons per segment is enough.
 
-    def from_pandas(self, *args, **kwargs):
-        raise NotImplementedError("from_XXXX methods are not implemented for VarEnergySpectrum")
+        Parameters
+        ----------
+        phot_per_segment: float or int
+            Number of photons in the segment
+        eint: iterable
+            Energy interval
+        """
+        if phot_per_segment < self.min_phot_per_segment:
+            streint = [f"{x:g}" for x in eint]
+            warnings.warn(
+                f"Low count rate in the {streint[0]}-{streint[1]} subject band: "
+                f"{phot_per_segment:g} ct/segment (<{self.min_phot_per_segment:g}). Skipping. "
+                "If you know what you're doing, you can set the `min_phot_per_segment` "
+                "parameter to a lower value, but in general, we recommend to use fewer "
+                "subject bands",
+            )
+            return False
+        return True
 
 
 class RmsSpectrum(VarEnergySpectrum):
@@ -487,8 +530,8 @@ class RmsSpectrum(VarEnergySpectrum):
 
     def __init__(
         self,
-        events,
-        energy_spec,
+        events=None,
+        energy_spec=None,
         ref_band=None,
         freq_interval=[0, 1],
         bin_time=1,
@@ -524,6 +567,8 @@ class RmsSpectrum(VarEnergySpectrum):
             sub_events = self._get_times_from_energy_range(self.events1, eint)
             countrate_sub = get_average_ctrate(sub_events, self.gti, self.segment_size)
             sub_power_noise = poisson_level(norm="abs", meanrate=countrate_sub)
+            if not self.check_phot_per_segment(countrate_sub * self.segment_size, eint):
+                continue
 
             # If we provided the `events2` array, calculate the rms from the
             # cospectrum, otherwise from the PDS
@@ -533,6 +578,8 @@ class RmsSpectrum(VarEnergySpectrum):
                 sub_events2 = self._get_times_from_energy_range(self.events2, eint)
                 countrate_sub2 = get_average_ctrate(sub_events2, self.gti, self.segment_size)
                 sub2_power_noise = poisson_level(norm="abs", meanrate=countrate_sub2)
+                if not self.check_phot_per_segment(countrate_sub2 * self.segment_size, eint):
+                    continue
 
                 # Calculate the cross spectrum
                 results = avg_cs_from_timeseries(
@@ -783,9 +830,9 @@ class LagSpectrum(VarEnergySpectrum):
     # events, freq_interval, energy_spec, ref_band = None
     def __init__(
         self,
-        events,
-        freq_interval,
-        energy_spec,
+        events=None,
+        freq_interval=None,
+        energy_spec=None,
         ref_band=None,
         bin_time=1,
         use_pi=False,
@@ -955,8 +1002,8 @@ class ComplexCovarianceSpectrum(VarEnergySpectrum):
 
     def __init__(
         self,
-        events,
-        energy_spec,
+        events=None,
+        energy_spec=None,
         ref_band=None,
         freq_interval=[0, 1],
         bin_time=1,
@@ -965,6 +1012,7 @@ class ComplexCovarianceSpectrum(VarEnergySpectrum):
         events2=None,
         norm="frac",
         return_complex=True,
+        min_phot_per_segment=10,
     ):
         self.norm = norm
         VarEnergySpectrum.__init__(
@@ -978,6 +1026,7 @@ class ComplexCovarianceSpectrum(VarEnergySpectrum):
             segment_size=segment_size,
             events2=events2,
             return_complex=return_complex,
+            min_phot_per_segment=min_phot_per_segment,
         )
 
     def _spectrum_function(self):
@@ -1007,6 +1056,10 @@ class ComplexCovarianceSpectrum(VarEnergySpectrum):
             # Extract events from the subject band
             sub_events = self._get_times_from_energy_range(self.events1, eint)
             countrate_sub = get_average_ctrate(sub_events, self.gti, self.segment_size)
+            phot_per_segment = countrate_sub * self.segment_size
+            if not self.check_phot_per_segment(phot_per_segment, eint):
+                continue
+
             sub_power_noise = poisson_level(norm="abs", meanrate=countrate_sub)
 
             results_cross = avg_cs_from_timeseries(
@@ -1131,8 +1184,8 @@ class CovarianceSpectrum(ComplexCovarianceSpectrum):
 
     def __init__(
         self,
-        events,
-        energy_spec,
+        events=None,
+        energy_spec=None,
         ref_band=None,
         freq_interval=[0, 1],
         bin_time=1,


Test patch:
diff --git a/stingray/tests/test_varenergyspectrum.py b/stingray/tests/test_varenergyspectrum.py
index e3767894..a1e5cf44 100644
--- a/stingray/tests/test_varenergyspectrum.py
+++ b/stingray/tests/test_varenergyspectrum.py
@@ -1,3 +1,4 @@
+import copy
 from multiprocessing import Event
 import os
 import numpy as np
@@ -89,6 +90,14 @@ class TestVarEnergySpectrum(object):
         with pytest.raises(TypeError):
             ref_int = VarEnergySpectrum(self.events, [0.0, 10000], (0.5, 5, 10, "log"), [0.3, 10])
 
+    @pytest.mark.parametrize("energy_spec", [2, "a"])
+    def test_invalid_energy_spec(self, energy_spec):
+        with pytest.raises(
+            ValueError,
+            match=f"Energy specification must be a tuple or a list .input: {energy_spec}.",
+        ):
+            DummyVarEnergy(self.events, [0.0, 10000], energy_spec=energy_spec)
+
     def test_ref_band_none(self):
         events = EventList(
             [0.09, 0.21, 0.23, 0.32, 0.4, 0.54], energy=[0, 0, 0, 0, 1, 1], gti=[[0, 0.65]]
@@ -210,20 +219,16 @@ class TestRmsAndCovSpectrum(object):
         assert np.all(np.iscomplex(spec.spectrum))
 
     @pytest.mark.parametrize("cross", [True, False])
-    @pytest.mark.parametrize("kind", ["rms", "cov", "lag"])
-    def test_empty_subband_cov(self, cross, kind):
+    def test_empty_subband_lag(self, cross):
         ev2 = None
         if cross:
             ev2 = self.test_ev2_small
 
-        if kind == "rms":
-            func = RmsSpectrum
-        elif kind == "lag":
-            func = LagSpectrum
-        elif kind == "cov":
-            func = ComplexCovarianceSpectrum
-
-        spec = func(
+        # Note: energy_spec is a list, so it's actually the edges
+        # of the energy bins. So, the covariance spectrum will be
+        # calculated in two bands: 0.3-12 keV and 12-15 keV. But
+        # the 12-15 keV band is empty (see definition of test_ev1_small)
+        spec = LagSpectrum(
             self.test_ev1_small,
             freq_interval=[0.00001, 0.1],
             energy_spec=[0.3, 12, 15],
@@ -235,6 +240,51 @@ class TestRmsAndCovSpectrum(object):
         good = ~np.isnan(spec.spectrum)
         assert np.count_nonzero(good) == 1
 
+    @pytest.mark.parametrize("cross", [True, False])
+    @pytest.mark.parametrize("kind", ["rms", "cov"])
+    def test_empty_subband_cov(self, cross, kind):
+        ev2 = None
+        if cross:
+            ev2 = self.test_ev2_small
+
+        if kind == "rms":
+            func = RmsSpectrum
+        elif kind == "cov":
+            func = ComplexCovarianceSpectrum
+        # Note: energy_spec is a list, so it's actually the edges
+        # of the energy bins. So, the covariance spectrum will be
+        # calculated in two bands: 0.3-12 keV and 12-15 keV. But
+        # the 12-15 keV band is empty (see definition of test_ev1_small)
+        with pytest.warns(UserWarning, match="Low count rate in the 12-15 subject band"):
+            spec = func(
+                self.test_ev1_small,
+                freq_interval=[0.00001, 0.1],
+                energy_spec=[0.3, 12, 15],
+                ref_band=[[0.3, 12]],
+                bin_time=self.bin_time / 2,
+                segment_size=200,
+                events2=ev2,
+            )
+        good = ~np.isnan(spec.spectrum)
+        assert np.count_nonzero(good) == 1
+
+    def test_empty_subband_cov_ev2(self):
+        ev2 = copy.deepcopy(self.test_ev2_small)
+        # We empty out only the second event list above 5 keV
+        ev2.filter_energy_range([0.3, 5], inplace=True)
+
+        with pytest.warns(UserWarning, match="Low count rate in the 5-12 subject band"):
+            spec = RmsSpectrum(
+                self.test_ev1_small,
+                freq_interval=[0.00001, 0.1],
+                energy_spec=[0.3, 5, 12],
+                bin_time=self.bin_time / 2,
+                segment_size=200,
+                events2=ev2,
+            )
+        good = ~np.isnan(spec.spectrum)
+        assert np.count_nonzero(good) == 1
+
     @pytest.mark.parametrize("norm", ["frac", "abs"])
     def test_correct_rms_values_vs_cross(self, norm):
         """The rms calculated with independent event lists (from the cospectrum)
@@ -355,6 +405,132 @@ class TestRmsAndCovSpectrum(object):
         assert np.all(np.isnan(rms.spectrum_error))
 
 
+import abc
+
+
+class BaseTestIO(abc.ABC):
+    @property
+    @abc.abstractmethod
+    def variant(self):
+        pass
+
+    @classmethod
+    def setup_class(cls):
+        if cls.variant == "rms":
+            cls.func = RmsSpectrum
+        elif cls.variant == "complcov":
+            cls.func = ComplexCovarianceSpectrum
+        elif cls.variant == "cov":
+            cls.func = CovarianceSpectrum
+        elif cls.variant == "lag":
+            cls.func = LagSpectrum
+        spec = cls.func(energy_spec=[0.3, 12, 15])
+        spec.freq_interval = [0.1, 0.2]
+        spec.ref_band = [0.3, 12]
+        spec.bin_time = 0.01
+        spec.segment_size = 100
+        spec.cross = cls.variant == "complcov"
+        cls.sting_obj = spec
+
+    def test_astropy_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        ts = so.to_astropy_table()
+        new_so = self.func.from_astropy_table(ts)
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_XARRAY")
+    def test_xarray_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        ts = so.to_xarray()
+        new_so = self.func.from_xarray(ts)
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_PANDAS")
+    def test_pandas_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        ts = so.to_pandas()
+        new_so = self.func.from_pandas(ts)
+        assert so == new_so
+
+    def test_astropy_roundtrip_empty(self):
+        # Set an attribute to a DummyStingrayObj. It will *not* be saved
+        so = self.func()
+        ts = so.to_astropy_table()
+        new_so = self.func.from_astropy_table(ts)
+        assert new_so.energy == []
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_XARRAY")
+    def test_xarray_roundtrip_empty(self):
+        so = self.func()
+        ts = so.to_xarray()
+        new_so = self.func.from_xarray(ts)
+        assert new_so.energy == []
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_PANDAS")
+    def test_pandas_roundtrip_empty(self):
+        so = self.func()
+        ts = so.to_pandas()
+        new_so = self.func.from_pandas(ts)
+        assert new_so.energy == []
+        assert so == new_so
+
+    @pytest.mark.skipif("not _HAS_H5PY")
+    def test_hdf_roundtrip(self):
+        so = copy.deepcopy(self.sting_obj)
+        so.write("dummy.hdf5")
+        new_so = so.read("dummy.hdf5")
+        os.unlink("dummy.hdf5")
+
+        assert so == new_so
+
+    def test_file_roundtrip_fits(self):
+        so = copy.deepcopy(self.sting_obj)
+        with pytest.warns(
+            UserWarning, match=".* output does not serialize the metadata at the moment"
+        ):
+            so.write("dummy.fits")
+        new_so = self.func.read("dummy.fits")
+        os.unlink("dummy.fits")
+        assert so == new_so
+
+    @pytest.mark.parametrize("fmt", ["ascii", "ascii.ecsv"])
+    def test_file_roundtrip(self, fmt):
+        so = copy.deepcopy(self.sting_obj)
+        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
+            so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = self.func.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        assert so == new_so
+
+    def test_file_roundtrip_pickle(self):
+        fmt = "pickle"
+        so = copy.deepcopy(self.sting_obj)
+        so.write(f"dummy.{fmt}", fmt=fmt)
+        new_so = self.func.read(f"dummy.{fmt}", fmt=fmt)
+        os.unlink(f"dummy.{fmt}")
+
+        assert so == new_so
+
+
+class TestCovarianceIO(BaseTestIO):
+    variant = "cov"
+
+
+class TestComplexCovarianceIO(BaseTestIO):
+    variant = "complcov"
+
+
+class TestRmsIO(BaseTestIO):
+    variant = "rms"
+
+
+class TestLagIO(BaseTestIO):
+    variant = "lag"
+
+
 @pytest.mark.slow
 class TestLagEnergySpectrum(object):
     @classmethod
@@ -427,80 +603,3 @@ class TestLagEnergySpectrum(object):
 
         assert np.all(np.isnan(lag.spectrum))
         assert np.all(np.isnan(lag.spectrum_error))
-
-
-class TestRoundTrip:
-    @classmethod
-    def setup_class(cls):
-        tstart = 0.0
-        tend = 100.0
-        nphot = 1000
-        alltimes = np.random.uniform(tstart, tend, nphot)
-        alltimes.sort()
-        cls.events = EventList(
-            alltimes, energy=np.random.uniform(0.3, 12, nphot), gti=[[tstart, tend]]
-        )
-        cls.vespec = DummyVarEnergy(
-            cls.events, [0.0, 10000], (0.5, 5, 10, "lin"), [0.3, 10], bin_time=0.1
-        )
-        cls.vespec.spectrum = np.zeros_like(cls.vespec.energy)
-        cls.vespec.spectrum_error = np.zeros_like(cls.vespec.energy)
-
-    def _check_equal(self, so, table):
-        for attr in ["energy", "spectrum", "spectrum_error"]:
-            assert np.allclose(getattr(so, attr), table[attr])
-
-        if hasattr(table, "meta"):
-            for attr in ["freq_interval"]:
-                assert getattr(so, attr) == table.meta[attr]
-        if hasattr(table, "attrs"):
-            for attr in ["freq_interval"]:
-                assert getattr(so, attr) == table.attrs[attr]
-
-    def test_astropy_export(self):
-        so = self.vespec
-        ts = so.to_astropy_table()
-        self._check_equal(so, ts)
-        with pytest.raises(NotImplementedError):
-            so.from_astropy_table(ts)
-
-    @pytest.mark.skipif("not _HAS_XARRAY")
-    def test_xarray_export(self):
-        so = self.vespec
-        ts = so.to_xarray()
-        self._check_equal(so, ts)
-        with pytest.raises(NotImplementedError):
-            so.from_xarray(ts)
-
-    @pytest.mark.skipif("not _HAS_PANDAS")
-    def test_pandas_export(self):
-        so = self.vespec
-        ts = so.to_pandas()
-        self._check_equal(so, ts)
-        with pytest.raises(NotImplementedError):
-            so.from_pandas(ts)
-
-    @pytest.mark.skipif("not _HAS_H5PY")
-    def test_hdf_export(self):
-        so = self.vespec
-        so.write("dummy.hdf5")
-        new_so = Table.read("dummy.hdf5")
-        os.unlink("dummy.hdf5")
-        self._check_equal(so, new_so)
-
-    @pytest.mark.parametrize("fmt", ["ascii.ecsv", "ascii", "fits"])
-    def test_file_export(self, fmt):
-        so = self.vespec
-        with pytest.warns(UserWarning, match=".* output does not serialize the metadata"):
-            so.write("dummy", fmt=fmt)
-        new_so = Table.read("dummy", format=fmt)
-        os.unlink("dummy")
-        self._check_equal(so, new_so)
-
-    @pytest.mark.parametrize("fmt", ["pickle"])
-    def test_file_export_pickle(self, fmt):
-        so = self.vespec
-        so.write("dummy", fmt=fmt)
-        new_so = so.read("dummy", fmt=fmt)
-        os.unlink("dummy")
-        self._check_equal(so, new_so.to_astropy_table())