Repository: StingraySoftware/stingray
Pull request #781 (2023-12-07 20:06:11)

Colors (as in color-color or hardness-intensity diagrams)
~Depends on #754~

pep8speaks: Hello @matteobachetti! Thanks for opening this PR. We checked the lines you've touched for [PEP 8](https://www.python.org/dev/peps/pep-0008) issues, and found:



* In the file [`stingray/base.py`](https://github.com/StingraySoftware/stingray/blob/28e78332698a01be44682880aabd75b347fffb35/stingray/base.py):

> [Line 2453:101](https://github.com/StingraySoftware/stingray/blob/28e78332698a01be44682880aabd75b347fffb35/stingray/base.py#L2453): [E501](https://duckduckgo.com/?q=pep8%20E501) line too long (121 > 100 characters)






























matteobachetti: @dhuppenkothen the problem is that some missions do not declare a time resolution. What do you use for dt then? 0 means simply that there is no known time resolution, or that it's probably irrelevant for the purpose of the analysis.
matteobachetti: Ok, I changed the logic to eliminate some of the confusion, and eliminated method-level imports as promised ;) 
matteobachetti: Note to self: fix the return value in the docstring of get_intensity_evolution 

Patch:
diff --git a/docs/changes/781.feature.rst b/docs/changes/781.feature.rst
new file mode 100644
index 00000000..c170f426
--- /dev/null
+++ b/docs/changes/781.feature.rst
@@ -0,0 +1,1 @@
+Calculate colors and intensities on a segment-by-segment basis in event lists
diff --git a/stingray/base.py b/stingray/base.py
index 9b0cc62f..954f5ff0 100644
--- a/stingray/base.py
+++ b/stingray/base.py
@@ -7,7 +7,6 @@ from collections import OrderedDict
 import pickle
 import warnings
 import copy
-import os
 import logging
 
 import numpy as np
@@ -33,8 +32,9 @@ from .gti import (
     gti_border_bins,
     get_btis,
     merge_gtis,
-    check_separate,
-    append_gtis,
+    get_total_gti_length,
+    bin_intervals_from_gtis,
+    time_intervals_from_gtis,
 )
 
 from typing import TYPE_CHECKING, Type, TypeVar, Union
@@ -247,8 +247,8 @@ class StingrayObject(object):
 
         This is useful for debugging, and for interactive use.
 
-        Optional parameters
-        -------------------
+        Other parameters
+        ----------------
         func_to_apply : function
             A function that modifies the attributes listed in ``attrs_to_apply``.
             It must return the modified attributes and a label to be printed.
@@ -1053,6 +1053,7 @@ class StingrayObject(object):
             raise IndexError("The index must be either an integer or a slice " "object !")
 
         new_ts = type(self)()
+
         for attr in self.meta_attrs():
             setattr(new_ts, attr, copy.deepcopy(getattr(self, attr)))
 
@@ -1062,6 +1063,17 @@ class StingrayObject(object):
         return new_ts
 
 
+def _ts_sum(ts):
+    """Sum the number of values of a time series object.
+
+    If it has a ``counts`` attribute, sum the counts. Otherwise, count the number
+    of time samples. Masked values are ignored.
+    """
+    if hasattr(ts, "counts"):
+        return np.sum(ts.counts[ts.mask])
+    return np.count_nonzero(ts.mask)
+
+
 class StingrayTimeseries(StingrayObject):
     """Basic class for time series data.
 
@@ -1128,7 +1140,10 @@ class StingrayTimeseries(StingrayObject):
 
     dt: float
         The time resolution of the measurements. Can be a scalar or an array attribute (useful
-        for non-evenly sampled data or events from different instruments)
+        for non-evenly sampled data or events from different instruments). It can also be 0, which
+        means that the time series is not evenly sampled and the effects of the time resolution are
+        considered negligible for the analysis. This is sometimes the case for events from
+        high-energy telescopes.
 
     mjdref : float
         The MJD used as a reference for the time array.
@@ -1216,8 +1231,12 @@ class StingrayTimeseries(StingrayObject):
 
     @property
     def mask(self):
-        if self._mask is None:
-            self._mask = create_gti_mask(self.time, self.gti, dt=self.dt)
+        if self._mask is not None:
+            return self._mask
+        if self._gti is not None:
+            self._mask = create_gti_mask(self.time, self._gti, dt=self.dt)
+        else:
+            self._mask = np.ones_like(self.time, dtype=bool)
         return self._mask
 
     @property
@@ -1291,6 +1310,19 @@ class StingrayTimeseries(StingrayObject):
             )
         return value
 
+    @property
+    def exposure(self):
+        """
+        Return the total exposure of the time series, i.e. the sum of the GTIs.
+
+        Returns
+        -------
+        total_exposure : float
+            The total exposure of the time series, in seconds.
+        """
+
+        return get_total_gti_length(self.gti)
+
     def __eq__(self, other_ts):
         return super().__eq__(other_ts)
 
@@ -1704,12 +1736,10 @@ class StingrayTimeseries(StingrayObject):
         >>> ts = StingrayTimeseries(time, array_attrs={"counts": count}, dt=1)
         >>> ts_new = ts.truncate(start=2, stop=8)
         >>> assert np.allclose(ts_new.counts, [30, 40, 50, 60, 70, 80])
-        >>> ts_new.time
-        array([3, 4, 5, 6, 7, 8])
+        >>> assert np.allclose(ts_new.time, [3, 4, 5, 6, 7, 8])
         >>> # Truncation can also be done by time values
         >>> ts_new = ts.truncate(start=6, method='time')
-        >>> ts_new.time
-        array([6, 7, 8, 9])
+        >>> assert np.allclose(ts_new.time, [6, 7, 8, 9])
         >>> assert np.allclose(ts_new.counts, [60, 70, 80, 90])
         """
 
@@ -2405,6 +2435,225 @@ class StingrayTimeseries(StingrayObject):
                 )
         return ax
 
+    def estimate_segment_size(self, min_counts=None, min_samples=None, even_sampling=None):
+        """Estimate a reasonable segment length for segment-by-segment analysis.
+
+        The user has to specify a criterion based on a minimum number of counts (if
+        the time series has a ``counts`` attribute) or a minimum number of time samples.
+        At least one between ``min_counts`` and ``min_samples`` must be specified.
+        In the special case of a time series with ``dt=0`` (event list-like, where each time
+        stamp correspond to a single count), the two definitions are equivalent.
+
+        Other Parameters
+        ----------------
+        min_counts : int
+            Minimum number of counts for each chunk. Optional (but needs ``min_samples``
+            if left unspecified). Only makes sense if the series has a ``counts`` attribute and
+            it is evenly sampled.
+        min_samples : int
+            Minimum number of time bins. Optional (but needs ``min_counts`` if left unspecified).
+        even_sampling : bool
+            Force the treatment of the data as evenly sampled or not. If None, the data are
+            considered evenly sampled if ``self.dt`` is larger than zero and the median
+            separation between subsequent times is within 1% of ``self.dt``.
+
+        Returns
+        -------
+        segment_size : float
+            The length of the light curve chunks that satisfies the conditions
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> time = np.arange(150)
+        >>> counts = np.zeros_like(time) + 3
+        >>> ts = StingrayTimeseries(time, counts=counts, dt=1)
+        >>> assert np.isclose(ts.estimate_segment_size(min_counts=10, min_samples=3), 4.0)
+        >>> assert np.isclose(ts.estimate_segment_size(min_counts=10, min_samples=5), 5.0)
+        >>> counts[2:4] = 1
+        >>> ts = StingrayTimeseries(time, counts=counts, dt=1)
+        >>> assert np.isclose(ts.estimate_segment_size(min_counts=3, min_samples=1), 3.0)
+        >>> # A slightly more complex example
+        >>> dt=0.2
+        >>> time = np.arange(0, 1000, dt)
+        >>> counts = np.random.poisson(100, size=len(time))
+        >>> ts = StingrayTimeseries(time, counts=counts, dt=dt)
+        >>> assert np.isclose(ts.estimate_segment_size(100, 2), 0.4)
+        >>> min_total_bins = 40
+        >>> assert np.isclose(ts.estimate_segment_size(100, 40), 8.0)
+        """
+        if min_counts is None and min_samples is None:
+            raise ValueError("You have to specify at least one of min_counts or min_samples")
+
+        mean_data_separation = np.median(np.diff(self.time))
+
+        if even_sampling is None:
+            # The time series is considered evenly sampled if the median separation between
+            # subsequent times is within 1% of the time resolution
+            even_sampling = False
+            if (
+                self.dt is not None
+                and self.dt > 0
+                and np.isclose(mean_data_separation, self.dt, rtol=0.01)
+            ):
+                even_sampling = True
+            logging.info(f"Data are {'not' if not even_sampling else ''} evenly sampled")
+
+        if min_counts is None:
+            if even_sampling and hasattr(self, "counts"):
+                min_counts = 0
+            else:
+                min_counts = min_samples
+
+        mean_ctrate = _ts_sum(self) / self.exposure
+
+        rough_estimate = np.ceil(min_counts / mean_ctrate)
+
+        # If data are evenly sampled, even sampling make the segment an integer multiple of dt.
+        # Otherwise, just use steps of 1 second.
+        if even_sampling:
+            step = self.dt
+        else:
+            step = 1.0
+
+        rough_estimate = np.ceil(min_counts / mean_ctrate / step) * step
+
+        segment_size = np.max([rough_estimate, min_samples * step])
+
+        keep_searching = True
+
+        while keep_searching:
+            start_times, stop_times, results = self.analyze_segments(_ts_sum, segment_size)
+            mincounts = np.min(results)
+            if mincounts >= min_counts:
+                keep_searching = False
+            else:
+                segment_size += step
+
+        return segment_size
+
+    def analyze_segments(self, func, segment_size, fraction_step=1, **kwargs):
+        """Analyze segments of the light curve with any function.
+
+        Parameters
+        ----------
+        func : function
+            Function accepting a :class:`StingrayTimeseries` object as single argument, plus
+            possible additional keyword arguments, and returning a number or a
+            tuple - e.g., ``(result, error)`` where both ``result`` and ``error`` are
+            numbers.
+        segment_size : float
+            Length in seconds of the light curve segments. If None, the full GTIs are considered
+            instead as segments.
+
+        Other parameters
+        ----------------
+        fraction_step : float
+            If the step is not a full ``segment_size`` but less (e.g. a moving window),
+            this indicates the ratio between step step and ``segment_size`` (e.g.
+            0.5 means that the window shifts of half ``segment_size``)
+        kwargs : keyword arguments
+            These additional keyword arguments, if present, they will be passed
+            to ``func``
+
+        Returns
+        -------
+        start_times : array
+            Lower time boundaries of all time segments.
+        stop_times : array
+            upper time boundaries of all segments.
+        result : array of N elements
+            The result of ``func`` for each segment of the light curve
+
+        Examples
+        --------
+        >>> import numpy as np
+        >>> time = np.arange(0, 10, 0.1)
+        >>> counts = np.zeros_like(time) + 10
+        >>> ts = StingrayTimeseries(time, counts=counts, dt=0.1)
+        >>> # Define a function that calculates the mean
+        >>> mean_func = lambda ts: np.mean(ts.counts)
+        >>> # Calculate the mean in segments of 5 seconds
+        >>> start, stop, res = ts.analyze_segments(mean_func, 5)
+        >>> len(res) == 2
+        True
+        >>> np.allclose(res, 10)
+        True
+        """
+
+        if segment_size is None:
+            start_times = self.gti[:, 0]
+            stop_times = self.gti[:, 1]
+            start = np.searchsorted(self.time, start_times)
+            stop = np.searchsorted(self.time, stop_times)
+        elif self.dt > 0:
+            start, stop = bin_intervals_from_gtis(
+                self.gti, segment_size, self.time, fraction_step=fraction_step, dt=self.dt
+            )
+            start_times = self.time[start] - 0.5 * self.dt
+            # Remember that stop is one element above the last element, because
+            # it's defined to be used in intervals start:stop
+            stop_times = self.time[stop - 1] + self.dt * 1.5
+        else:
+            start_times, stop_times = time_intervals_from_gtis(
+                self.gti, segment_size, fraction_step=fraction_step
+            )
+            start = np.searchsorted(self.time, start_times)
+            stop = np.searchsorted(self.time, stop_times)
+
+        results = []
+        for i, (st, sp, tst, tsp) in enumerate(zip(start, stop, start_times, stop_times)):
+            if sp - st <= 1:
+                warnings.warn(
+                    f"Segment {i} ({tst}--{tsp}) has one data point or less. Skipping it."
+                )
+                res = np.nan
+            else:
+                lc_filt = self[st:sp]
+                lc_filt.gti = np.asarray([[tst, tsp]])
+
+                res = func(lc_filt, **kwargs)
+            results.append(res)
+
+        results = np.array(results)
+
+        if len(results.shape) == 2:
+            results = [results[:, i] for i in range(results.shape[1])]
+        return start_times, stop_times, results
+
+    def analyze_by_gti(self, func, fraction_step=1, **kwargs):
+        """Analyze the light curve with any function, on a GTI-by-GTI base.
+
+        Parameters
+        ----------
+        func : function
+            Function accepting a :class:`StingrayTimeseries` object as single argument, plus
+            possible additional keyword arguments, and returning a number or a
+            tuple - e.g., ``(result, error)`` where both ``result`` and ``error`` are
+            numbers.
+
+        Other parameters
+        ----------------
+        fraction_step : float
+            By default, segments do not overlap (``fraction_step`` = 1). If ``fraction_step`` < 1,
+            then the start points of consecutive segments are ``fraction_step * segment_size``
+            apart, and consecutive segments overlap. For example, for ``fraction_step`` = 0.5,
+            the window shifts one half of ``segment_size``)
+        kwargs : keyword arguments
+            These additional keyword arguments, if present, they will be passed
+            to ``func``
+
+        Returns
+        -------
+        start_times : array
+            Lower time boundaries of all time segments.
+        stop_times : array
+            upper time boundaries of all segments.
+        result : array of N elements
+            The result of ``func`` for each segment of the light curve
+        """
+        return self.analyze_segments(func, segment_size=None, fraction_step=fraction_step, **kwargs)
+
 
 def interpret_times(time: TTime, mjdref: float = 0) -> tuple[npt.ArrayLike, float]:
     """Understand the format of input times, and return seconds from MJDREF
diff --git a/stingray/events.py b/stingray/events.py
index 6214b2f7..a45d0bde 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -4,24 +4,19 @@ Definition of :class:`EventList`.
 :class:`EventList` is used to handle photon arrival times.
 """
 
-import copy
 import logging
-import pickle
 import warnings
-from collections.abc import Iterable
 
 import numpy as np
-import numpy.random as ra
-from astropy.table import Table
 
 from stingray.utils import _int_sum_non_zero
 
-from .base import StingrayObject, StingrayTimeseries
+from .base import StingrayTimeseries
 from .filters import get_deadtime_mask
-from .gti import append_gtis, check_separate, cross_gtis, generate_indices_of_boundaries
+from .gti import generate_indices_of_boundaries
 from .io import load_events_and_gtis
 from .lightcurve import Lightcurve
-from .utils import assign_value_if_none, simon, njit
+from .utils import simon, njit
 from .utils import histogram
 
 __all__ = ["EventList"]
@@ -627,6 +622,25 @@ class EventList(StingrayTimeseries):
 
         return super().read(filename=filename, fmt=fmt)
 
+    def get_energy_mask(self, energy_range, use_pi=False):
+        """Get a mask corresponding to events with a given energy range.
+
+        Parameters
+        ----------
+        energy_range: [float, float]
+            Energy range in keV, or in PI channel (if ``use_pi`` is True)
+
+        Other Parameters
+        ----------------
+        use_pi : bool, default False
+            Use PI channel instead of energy in keV
+        """
+        if use_pi:
+            energies = self.pi
+        else:
+            energies = self.energy
+        return (energies >= energy_range[0]) & (energies < energy_range[1])
+
     def filter_energy_range(self, energy_range, inplace=False, use_pi=False):
         """Filter the event list from a given energy range.
 
@@ -654,12 +668,7 @@ class EventList(StingrayTimeseries):
         >>> assert np.allclose(events.time, [0, 1])
 
         """
-        if use_pi:
-            energies = self.pi
-        else:
-            energies = self.energy
-        mask = (energies >= energy_range[0]) & (energies < energy_range[1])
-
+        mask = self.get_energy_mask(energy_range, use_pi=use_pi)
         return self.apply_mask(mask, inplace=inplace)
 
     def apply_deadtime(self, deadtime, inplace=False, **kwargs):
@@ -714,3 +723,80 @@ class EventList(StingrayTimeseries):
             new_ev = [new_ev, retall]
 
         return new_ev
+
+    def get_color_evolution(self, energy_ranges, segment_size=None, use_pi=False):
+        """Compute the color in equal-length segments of the event list.
+
+        Parameters
+        ----------
+        energy_ranges : 2x2 list
+            List of energy ranges to compute the color:
+            ``[[en1_min, en1_max], [en2_min, en2_max]]``
+        segment_size : float
+            Segment size in seconds. If None, the full GTIs are considered
+            instead as segments.
+
+        Other Parameters
+        ----------------
+        use_pi : bool, default False
+            Use PI channel instead of energy in keV
+
+        Returns
+        -------
+        color : array-like
+            Array of colors, computed in each segment as the ratio of the
+            counts in the second energy range to the counts in the first energy
+            range.
+        """
+        if energy_ranges is None or np.shape(energy_ranges) != (2, 2):
+            raise ValueError("Energy ranges must be specified as a 2x2 array")
+
+        def color(ev):
+            mask1 = ev.get_energy_mask(energy_ranges[0], use_pi=use_pi)
+            mask2 = ev.get_energy_mask(energy_ranges[1], use_pi=use_pi)
+            en1_ct = np.count_nonzero(mask1)
+            en2_ct = np.count_nonzero(mask2)
+
+            color = en2_ct / en1_ct
+            color_err = color * (np.sqrt(en1_ct) / en1_ct + np.sqrt(en2_ct) / en2_ct)
+            return color, color_err
+
+        starts, stops, (colors, color_errs) = self.analyze_segments(color, segment_size)
+
+        return starts, stops, colors, color_errs
+
+    def get_intensity_evolution(self, energy_range, segment_size=None, use_pi=False):
+        """Compute the intensity in equal-length segments (or full GTIs) of the event list.
+
+        Parameters
+        ----------
+        energy_range : ``[en1_min, en1_max]``
+            Energy range to compute the intensity
+        segment_size : float
+            Segment size in seconds. If None, the full GTIs are considered
+            instead as segments.
+
+        Other Parameters
+        ----------------
+        use_pi : bool, default False
+            Use PI channel instead of energy in keV
+
+        Returns
+        -------
+        intensity : array-like
+            Array of intensities (in counts/s), computed in each segment.
+        """
+        if energy_range is None or np.shape(energy_range) != (2,):
+            raise ValueError("Energy ranges must be specified as a 2-element list")
+
+        def intensity(ev):
+            mask1 = ev.get_energy_mask(energy_range, use_pi=use_pi)
+            en1_ct = np.count_nonzero(mask1)
+            segment_size = ev.gti[0, 1] - ev.gti[0, 0]
+            rate = en1_ct / segment_size
+            rate_err = np.sqrt(en1_ct) / segment_size
+            return rate, rate_err
+
+        starts, stops, (rate, rate_err) = self.analyze_segments(intensity, segment_size)
+
+        return starts, stops, rate, rate_err
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index f5393338..86be3734 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -5,14 +5,11 @@ Definition of :class::class:`Lightcurve`.
 or to save existing light curves in a class that's easy to use.
 """
 import os
-import copy
 import logging
 import warnings
-import pickle
 from collections.abc import Iterable
 
 import numpy as np
-import matplotlib.pyplot as plt
 from astropy.table import Table
 from astropy.time import TimeDelta, Time
 from astropy import units as u
@@ -21,11 +18,9 @@ from stingray.base import StingrayTimeseries, reduce_precision_if_extended
 import stingray.utils as utils
 from stingray.exceptions import StingrayError
 from stingray.gti import (
-    bin_intervals_from_gtis,
     check_gtis,
     create_gti_mask,
     cross_two_gtis,
-    gti_border_bins,
     join_gtis,
 )
 from stingray.utils import (
@@ -1292,21 +1287,25 @@ class Lightcurve(StingrayTimeseries):
         warnings.warn("This function was renamed to estimate_segment_size", DeprecationWarning)
         return self.estimate_segment_size(*args, **kwargs)
 
-    def estimate_segment_size(self, min_total_counts=100, min_time_bins=100):
-        """Estimate a reasonable segment length for chunk-by-chunk analysis.
+    def estimate_segment_size(self, min_counts=100, min_samples=100, even_sampling=None):
+        """Estimate a reasonable segment length for segment-by-segment analysis.
 
-        Choose a reasonable length for time segments, given a minimum number of total
-        counts in the segment, and a minimum number of time bins in the segment.
-
-        The user specifies a condition on the total counts in each segment and
-        the minimum number of time bins.
+        The user has to specify a criterion based on a minimum number of counts (if
+        the time series has a ``counts`` attribute) or a minimum number of time samples.
+        At least one between ``min_counts`` and ``min_samples`` must be specified.
 
         Other Parameters
         ----------------
-        min_total_counts : int
-            Minimum number of counts for each chunk
-        min_time_bins : int
-            Minimum number of time bins
+        min_counts : int
+            Minimum number of counts for each chunk. Optional (but needs ``min_samples``
+            if left unspecified). Only makes sense if the series has a ``counts`` attribute and
+            it is evenly sampled.
+        min_samples : int
+            Minimum number of time bins. Optional (but needs ``min_counts`` if left unspecified).
+        even_sampling : bool
+            Force the treatment of the data as evenly sampled or not. If None, the data are
+            considered evenly sampled if ``self.dt`` is larger than zero and the median
+            separation between subsequent times is within 1% of ``self.dt``.
 
         Returns
         -------
@@ -1320,11 +1319,11 @@ class Lightcurve(StingrayTimeseries):
         >>> count = np.zeros_like(time) + 3
         >>> lc = Lightcurve(time, count, dt=1)
         >>> assert np.isclose(
-        ...     lc.estimate_segment_size(min_total_counts=10, min_time_bins=3), 4)
-        >>> assert np.isclose(lc.estimate_segment_size(min_total_counts=10, min_time_bins=5), 5)
+        ...     lc.estimate_segment_size(min_counts=10, min_samples=3), 4)
+        >>> assert np.isclose(lc.estimate_segment_size(min_counts=10, min_samples=5), 5)
         >>> count[2:4] = 1
         >>> lc = Lightcurve(time, count, dt=1)
-        >>> assert np.isclose(lc.estimate_segment_size(min_total_counts=3, min_time_bins=1), 3)
+        >>> assert np.isclose(lc.estimate_segment_size(min_counts=3, min_samples=1), 3)
         >>> # A slightly more complex example
         >>> dt=0.2
         >>> time = np.arange(0, 1000, dt)
@@ -1334,25 +1333,14 @@ class Lightcurve(StingrayTimeseries):
         >>> min_total_bins = 40
         >>> assert np.isclose(lc.estimate_segment_size(100, 40), 8.0)
         """
-
-        rough_estimate = np.ceil(min_total_counts / self.meancounts) * self.dt
-
-        segment_size = np.max([rough_estimate, min_time_bins * self.dt])
-
-        keep_searching = True
-        while keep_searching:
-            start_times, stop_times, results = self.analyze_lc_chunks(segment_size, np.sum)
-            mincounts = np.min(results)
-            if mincounts >= min_total_counts:
-                keep_searching = False
-            else:
-                segment_size += self.dt
-
-        return segment_size
+        return super().estimate_segment_size(min_counts, min_samples, even_sampling=even_sampling)
 
     def analyze_lc_chunks(self, segment_size, func, fraction_step=1, **kwargs):
         """Analyze segments of the light curve with any function.
 
+        .. deprecated:: 2.0
+            Use :meth:`Lightcurve.analyze_segments(func, segment_size)` instead.
+
         Parameters
         ----------
         segment_size : float
@@ -1366,9 +1354,10 @@ class Lightcurve(StingrayTimeseries):
         Other parameters
         ----------------
         fraction_step : float
-            If the step is not a full ``segment_size`` but less (e.g. a moving window),
-            this indicates the ratio between step step and ``segment_size`` (e.g.
-            0.5 means that the window shifts of half ``segment_size``)
+            By default, segments do not overlap (``fraction_step`` = 1). If ``fraction_step`` < 1,
+            then the start points of consecutive segments are ``fraction_step * segment_size``
+            apart, and consecutive segments overlap. For example, for ``fraction_step`` = 0.5,
+            the window shifts one half of ``segment_size``)
         kwargs : keyword arguments
             These additional keyword arguments, if present, they will be passed
             to ``func``
@@ -1381,40 +1370,11 @@ class Lightcurve(StingrayTimeseries):
             upper time boundaries of all segments.
         result : array of N elements
             The result of ``func`` for each segment of the light curve
-
-        Examples
-        --------
-        >>> import numpy as np
-        >>> time = np.arange(0, 10, 0.1)
-        >>> counts = np.zeros_like(time) + 10
-        >>> lc = Lightcurve(time, counts, dt=0.1)
-        >>> # Define a function that calculates the mean
-        >>> mean_func = lambda x: np.mean(x)
-        >>> # Calculate the mean in segments of 5 seconds
-        >>> start, stop, res = lc.analyze_lc_chunks(5, mean_func)
-        >>> assert len(res) == 2
-        >>> assert np.allclose(res, 10)
         """
-        start, stop = bin_intervals_from_gtis(
-            self.gti, segment_size, self.time, fraction_step=fraction_step, dt=self.dt
+        warnings.warn(
+            "The analyze_lc_chunks method was superseded by analyze_segments", DeprecationWarning
         )
-        start_times = self.time[start] - 0.5 * self.dt
-
-        # Remember that stop is one element above the last element, because
-        # it's defined to be used in intervals start:stop
-        stop_times = self.time[stop - 1] + self.dt * 1.5
-
-        results = []
-        for i, (st, sp) in enumerate(zip(start, stop)):
-            lc_filt = self[st:sp]
-            res = func(lc_filt, **kwargs)
-            results.append(res)
-
-        results = np.array(results)
-
-        if len(results.shape) == 2:
-            results = [results[:, i] for i in range(results.shape[1])]
-        return start_times, stop_times, results
+        return super().analyze_segments(func, segment_size, fraction_step=fraction_step, **kwargs)
 
     def to_lightkurve(self):
         """


Test patch:
diff --git a/stingray/tests/test_base.py b/stingray/tests/test_base.py
index 56624dc5..ecdc8a18 100644
--- a/stingray/tests/test_base.py
+++ b/stingray/tests/test_base.py
@@ -1,33 +1,15 @@
 import os
+import importlib
 import copy
 import pytest
 import numpy as np
 import matplotlib.pyplot as plt
 from stingray.base import StingrayObject, StingrayTimeseries
 
-_HAS_XARRAY = _HAS_PANDAS = _HAS_H5PY = _HAS_YAML = True
-
-try:
-    import xarray
-    from xarray import Dataset
-except ImportError:
-    _HAS_XARRAY = False
-
-try:
-    import pandas
-    from pandas import DataFrame
-except ImportError:
-    _HAS_PANDAS = False
-
-try:
-    import h5py
-except ImportError:
-    _HAS_H5PY = False
-
-try:
-    import yaml
-except ImportError:
-    _HAS_YAML = False
+_HAS_XARRAY = importlib.util.find_spec("xarray") is not None
+_HAS_PANDAS = importlib.util.find_spec("pandas") is not None
+_HAS_H5PY = importlib.util.find_spec("h5py") is not None
+_HAS_YAML = importlib.util.find_spec("yaml") is not None
 
 
 class DummyStingrayObj(StingrayObject):
@@ -68,8 +50,10 @@ class TestStingrayObject:
         sting_obj.panesapa = [[41, 25], [98, 3]]
         cls.sting_obj = sting_obj
 
-    def test_print(self):
+    def test_print(self, capsys):
         print(self.sting_obj)
+        captured = capsys.readouterr()
+        assert "guefus" in captured.out
 
     def test_preliminary(self):
         assert np.allclose(self.sting_obj.guefus, self.arr)
@@ -406,8 +390,10 @@ class TestStingrayTimeseries:
         cls.sting_obj = sting_obj
         cls.sting_obj_highp = sting_obj_highp
 
-    def test_print(self):
+    def test_print(self, capsys):
         print(self.sting_obj)
+        captured = capsys.readouterr()
+        assert "59777" in captured.out
 
     def test_invalid_instantiation(self):
         with pytest.raises(ValueError, match="Lengths of time and guefus must be equal"):
@@ -415,6 +401,14 @@ class TestStingrayTimeseries:
         with pytest.raises(ValueError, match="Lengths of time and guefus must be equal"):
             StingrayTimeseries(time=np.arange(10), array_attrs=dict(guefus=np.zeros((5, 2))))
 
+    def test_mask_is_none_then_isnt_no_gti(self):
+        ts = copy.deepcopy(self.sting_obj)
+        assert ts._mask is None
+        # Unset GTIs
+        ts.gti = None
+        # But when I use the mask property, it's an array
+        assert np.array_equal(ts.mask, np.ones(len(ts.time), dtype=bool))
+
     def test_apply_mask(self):
         ts = copy.deepcopy(self.sting_obj)
         mask = [True, True] + 8 * [False]
@@ -963,7 +957,7 @@ class TestStingrayTimeseries:
         so.guefus = np.random.randint(0, 4, 3)
         so.panesapa = np.random.randint(5, 9, (6, 2))
         with pytest.warns(
-            UserWarning, match=f".* output does not serialize the metadata at the moment"
+            UserWarning, match=".* output does not serialize the metadata at the moment"
         ):
             so.write(f"dummy.{fmt}", fmt=fmt)
         new_so = StingrayTimeseries.read(f"dummy.{fmt}", fmt=fmt)
@@ -1048,8 +1042,10 @@ class TestStingrayTimeseriesSubclass:
         sting_obj.gti = np.asarray([[-0.5, 2.5]])
         cls.sting_obj = sting_obj
 
-    def test_print(self):
+    def test_print(self, capsys):
         print(self.sting_obj)
+        captured = capsys.readouterr()
+        assert "59777" in captured.out
 
     def test_astropy_roundtrip(self):
         so = copy.deepcopy(self.sting_obj)
@@ -1421,3 +1417,72 @@ class TestFillBTI(object):
         lc_like_filt.gti = np.asarray([[0, 900], [950, 999], [999.5, 1000]])
         lc_new = lc_like_filt.fill_bad_time_intervals()
         assert np.allclose(lc_new.gti, self.gti)
+
+
+class TestAnalyzeChunks(object):
+    @classmethod
+    def setup_class(cls):
+        cls.time = np.arange(150)
+        counts = np.zeros_like(cls.time) + 3
+        cls.ts = StingrayTimeseries(cls.time, counts=counts, dt=1)
+        cls.ts_no_dt = StingrayTimeseries(cls.time, counts=counts, dt=0)
+        cls.ts_no_counts = StingrayTimeseries(cls.time, dt=0)
+
+    def test_invalid_input(self):
+        with pytest.raises(ValueError, match="You have to specify at least one of"):
+            self.ts.estimate_segment_size()
+        with pytest.raises(ValueError, match="You have to specify at least one of"):
+            self.ts_no_dt.estimate_segment_size()
+
+    def test_no_total_counts(self):
+        assert self.ts.estimate_segment_size(min_samples=2) == 2
+        assert self.ts_no_dt.estimate_segment_size(min_samples=2) == 2
+        assert self.ts_no_counts.estimate_segment_size(min_samples=2) == 2
+
+    def test_estimate_segment_size(self):
+        # Here, the total counts dominate
+        assert self.ts.estimate_segment_size(min_counts=10, min_samples=3) == 4.0
+
+    def test_estimate_segment_size_more_bins(self):
+        # Here, the time bins dominate
+        assert self.ts.estimate_segment_size(min_counts=10, min_samples=5) == 5.0
+
+    def test_estimate_segment_size_lower_counts(self):
+        counts = np.zeros_like(self.time) + 3
+        counts[2:4] = 1
+        ts = StingrayTimeseries(self.time, counts=counts, dt=1)
+        assert ts.estimate_segment_size(min_counts=3, min_samples=1) == 3.0
+
+    def test_estimate_segment_size_lower_dt(self):
+        # A slightly more complex example
+        dt = 0.2
+        time = np.arange(0, 1000, dt)
+        counts = np.random.poisson(100, size=len(time))
+        ts = StingrayTimeseries(time, counts=counts, dt=dt)
+        assert ts.estimate_segment_size(min_counts=100, min_samples=2) == 0.4
+
+        assert ts.estimate_segment_size(100, min_samples=40) == 8.0
+
+    def test_analyze_segments_bad_intv(self):
+        ts = StingrayTimeseries(time=np.arange(10), dt=1, gti=[[-0.5, 0.5], [1.5, 10.5]])
+
+        def func(x):
+            return np.size(x.time)
+
+        # I do not specify the segment_size, which means results will be calculated per-GTI
+        with pytest.warns(UserWarning, match="has one data point or less."):
+            _, _, results = ts.analyze_segments(func, segment_size=None)
+        # the first GTI contains only one bin, the result will be invalid
+        assert np.isnan(results[0])
+
+    def test_analyze_segments_by_gti(self):
+        ts = StingrayTimeseries(time=np.arange(11), dt=1, gti=[[-0.5, 5.5], [6.5, 10.5]])
+
+        def func(x):
+            return np.size(x.time)
+
+        _, _, results_as = ts.analyze_segments(func, segment_size=None)
+        _, _, results_ag = ts.analyze_by_gti(func)
+
+        assert np.allclose(results_as, results_ag)
+        assert np.allclose(results_as, [6, 4])
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index f2bc04af..83367e22 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -1,8 +1,8 @@
-import warnings
-import numpy as np
+import importlib
 import copy
 import os
 import pytest
+import numpy as np
 from astropy.time import Time
 
 from ..events import EventList
@@ -11,34 +11,10 @@ from ..lightcurve import Lightcurve
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
 
-_H5PY_INSTALLED = True
-_HAS_YAML = True
-_HAS_XARRAY = _HAS_PANDAS = True
-
-try:
-    import h5py
-except ImportError:
-    _H5PY_INSTALLED = False
-
-import astropy.timeseries
-from astropy.timeseries import TimeSeries
-
-try:
-    import xarray
-    from xarray import Dataset
-except ImportError:
-    _HAS_XARRAY = False
-
-try:
-    import pandas
-    from pandas import DataFrame
-except ImportError:
-    _HAS_PANDAS = False
-
-try:
-    import yaml
-except ImportError:
-    _HAS_YAML = False
+_HAS_XARRAY = importlib.util.find_spec("xarray") is not None
+_HAS_PANDAS = importlib.util.find_spec("pandas") is not None
+_HAS_H5PY = importlib.util.find_spec("h5py") is not None
+_HAS_YAML = importlib.util.find_spec("yaml") is not None
 
 
 class TestEvents(object):
@@ -56,7 +32,7 @@ class TestEvents(object):
             _ = EventList(self.time, self.counts, gti=self.gti, bubu="settete")
         assert np.any(["Unrecognized keywords:" in r.message.args[0] for r in record])
 
-    def test_warn_wrong_keywords(self):
+    def test_warn_wrong_keywords_ncounts(self):
         with pytest.warns(DeprecationWarning, match="The ncounts keyword does nothing"):
             _ = EventList(self.time, self.counts, gti=self.gti, ncounts=10)
 
@@ -210,7 +186,7 @@ class TestEvents(object):
         assert np.allclose(ev.time, self.time)
         os.remove("ev.pickle")
 
-    @pytest.mark.skipif("not _H5PY_INSTALLED")
+    @pytest.mark.skipif("not _HAS_H5PY")
     def test_io_with_hdf5_auto(self):
         ev = EventList(time=self.time, mjdref=54000)
         ev.write("ev.hdf5")
@@ -219,7 +195,7 @@ class TestEvents(object):
         assert np.allclose(ev.time, self.time)
         os.remove("ev.hdf5")
 
-    @pytest.mark.skipif("not _H5PY_INSTALLED")
+    @pytest.mark.skipif("not _HAS_H5PY")
     def test_io_with_hdf5(self):
         ev = EventList(time=self.time, mjdref=54000)
         ev.write("ev.hdf5", fmt="hdf5")
@@ -598,3 +574,65 @@ class TestFilters(object):
         assert np.allclose(filt_events.time, expected)
         assert np.allclose(filt_events.pi, 1)
         assert np.allclose(filt_events.energy, 1)
+
+
+class TestColors(object):
+    @classmethod
+    def setup_class(cls):
+        cls.events = EventList(
+            time=np.arange(100000) + 0.5, energy=np.random.choice([2, 5], 100000), gti=[[0, 100000]]
+        )
+
+    def test_bad_interval_color(self):
+        with pytest.raises(ValueError, match=" 2x2 array"):
+            self.events.get_color_evolution([[0, 3], [4, 6], [7, 8]], 10000)
+        with pytest.raises(ValueError, match=" 2x2 array"):
+            self.events.get_color_evolution([[0, 3, 8]], 10000)
+        with pytest.raises(ValueError, match=" 2x2 array"):
+            self.events.get_color_evolution([0], 10000)
+        with pytest.raises(ValueError, match=" 2x2 array"):
+            self.events.get_color_evolution([[0, 1]], 10000)
+
+    def test_bad_interval_intensity(self):
+        with pytest.raises(ValueError, match="2-element list"):
+            self.events.get_intensity_evolution([[0, 3], [4, 6], [7, 8]], 10000)
+        with pytest.raises(ValueError, match="2-element list"):
+            self.events.get_intensity_evolution([[0, 3, 8]], 10000)
+        with pytest.raises(ValueError, match="2-element list"):
+            self.events.get_intensity_evolution([0], 10000)
+        with pytest.raises(ValueError, match="2-element list"):
+            self.events.get_intensity_evolution([[0, 1]], 10000)
+
+    def test_colors(self):
+        start, stop, colors, color_errs = self.events.get_color_evolution([[0, 3], [4, 6]], 10000)
+        # 5000 / 5000 = 1
+        # 2 x sqrt(5000) / 5000 = 0.0282
+        assert np.allclose(colors, 1, rtol=0.1)
+        assert np.allclose(color_errs, 0.0282, atol=0.003)
+        assert np.allclose(start, np.arange(10) * 10000)
+        assert np.allclose(stop, np.arange(1, 11) * 10000)
+
+    def test_colors_no_segment(self):
+        start, stop, colors, color_errs = self.events.get_color_evolution([[0, 3], [4, 6]])
+        # 50000 / 50000 = 1
+        # 2 x sqrt(50000) / 50000 = 0.0089
+        assert np.allclose(colors, 1, rtol=0.1)
+        assert np.allclose(color_errs, 0.0089, atol=0.001)
+        assert np.allclose(start, 0)
+        assert np.allclose(stop, 100000)
+
+    def test_intensity(self):
+        start, stop, rate, rate_errs = self.events.get_intensity_evolution([0, 6], 10000)
+
+        assert np.allclose(rate, 1, rtol=0.1)
+        assert np.allclose(rate_errs, 0.01, atol=0.003)
+        assert np.allclose(start, np.arange(10) * 10000)
+        assert np.allclose(stop, np.arange(1, 11) * 10000)
+
+    def test_intensity_no_segment(self):
+        start, stop, rate, rate_errs = self.events.get_intensity_evolution([0, 6])
+
+        assert np.allclose(rate, 1, rtol=0.1)
+        assert np.allclose(rate_errs, 0.003, atol=0.001)
+        assert np.allclose(start, 0)
+        assert np.allclose(stop, 100000)
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 6b770bbd..78d82e1f 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -1,14 +1,13 @@
 import os
+import importlib
 import copy
 import numpy as np
 import pytest
 import warnings
-import os
 import matplotlib.pyplot as plt
 from numpy.testing import assert_allclose
 import astropy.units as u
 from astropy.time import Time
-import astropy.timeseries
 from astropy.timeseries import TimeSeries
 from astropy.table import Table
 import scipy.stats
@@ -19,31 +18,12 @@ from stingray.gti import create_gti_mask
 
 np.random.seed(20150907)
 
-_H5PY_INSTALLED = True
-_HAS_LIGHTKURVE = True
-_HAS_YAML = True
-_IS_WINDOWS = os.name == "nt"
-_HAS_ULTRANEST = True
-
-try:
-    import h5py
-except ImportError:
-    _H5PY_INSTALLED = False
-
-try:
-    import lightkurve
-except ImportError:
-    _HAS_LIGHTKURVE = False
+_H5PY_INSTALLED = importlib.util.find_spec("h5py") is not None
+_HAS_LIGHTKURVE = importlib.util.find_spec("lightkurve") is not None
+_HAS_YAML = importlib.util.find_spec("yaml") is not None
+_HAS_ULTRANEST = importlib.util.find_spec("ultranest") is not None
 
-try:
-    import yaml
-except ImportError:
-    _HAS_YAML = False
-
-try:
-    import ultranest
-except ImportError:
-    _HAS_ULTRANEST = False
+_IS_WINDOWS = os.name == "nt"
 
 curdir = os.path.abspath(os.path.dirname(__file__))
 datadir = os.path.join(curdir, "data")
@@ -117,7 +97,7 @@ class TestProperties(object):
         assert np.any(["On Windows, the size of an integer" in r.message.args[0] for r in record])
 
     @pytest.mark.skipif("_IS_WINDOWS")
-    def test_warn_on_windows(self, monkeypatch):
+    def test_warn_on_windows_monkeypatching_elsewhere(self, monkeypatch):
         monkeypatch.setattr(os, "name", "nt")
         with pytest.warns(UserWarning) as record:
             _ = Lightcurve(self.lc.time, self.lc.counts, gti=self.lc.gti)
@@ -262,7 +242,7 @@ class TestProperties(object):
     def test_add_data_to_empty_lightcurve_wrong(self, attr):
         lc0 = Lightcurve()
         lc0.time = [1, 2, 3]
-        with pytest.raises(ValueError, match=f".*the same shape as the time array"):
+        with pytest.raises(ValueError, match=".*the same shape as the time array"):
             setattr(lc0, attr, [1, 2, 3, 4])
 
     @pytest.mark.parametrize("attr", "counts,countrate".split(","))
@@ -304,7 +284,8 @@ class TestChunks(object):
         cls.lc = Lightcurve(times, counts, gti=cls.gti)
 
     def test_analyze_lc_chunks_fvar_fracstep(self):
-        start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun, fraction_step=0.5)
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun, fraction_step=0.5)
         # excess_variance returns fvar and fvar_err
         fvar, fvar_err = res
 
@@ -314,20 +295,55 @@ class TestChunks(object):
         assert np.all(fvar > fvar_err)
 
     def test_analyze_lc_chunks_nvar_fracstep(self):
-        start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun, fraction_step=0.5)
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun, fraction_step=0.5)
         # excess_variance returns fvar and fvar_err
         fvar, fvar_err = res
-        start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun, fraction_step=0.5)
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun, fraction_step=0.5)
         # excess_variance returns fvar and fvar_err
         nevar, nevar_err = res
         assert np.allclose(nevar, fvar**2, rtol=0.01)
 
     def test_analyze_lc_chunks_nvar_fracstep_mean(self):
-        start, stop, mean = self.lc.analyze_lc_chunks(20, np.mean, fraction_step=0.5)
-        start, stop, res = self.lc.analyze_lc_chunks(20, evar_fun, fraction_step=0.5)
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, mean = self.lc.analyze_lc_chunks(20, np.mean, fraction_step=0.5)
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, res = self.lc.analyze_lc_chunks(20, evar_fun, fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        evar, evar_err = res
+        with pytest.warns(DeprecationWarning, match="The analyze_lc_chunks method was superseded"):
+            start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun, fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        nevar, nevar_err = res
+        assert np.allclose(nevar * mean**2, evar, rtol=0.01)
+        assert np.allclose(nevar_err * mean**2, evar_err, rtol=0.01)
+
+    def test_analyze_segments_fvar_fracstep(self):
+        start, stop, res = self.lc.analyze_segments(fvar_fun, 20, fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        fvar, fvar_err = res
+
+        assert np.allclose(start[0], self.gti[0, 0])
+        assert np.all(fvar > 0)
+        # This must be a clear measurement of fvar
+        assert np.all(fvar > fvar_err)
+
+    def test_analyze_segments_nvar_fracstep(self):
+        start, stop, res = self.lc.analyze_segments(fvar_fun, 20, fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        fvar, fvar_err = res
+        start, stop, res = self.lc.analyze_segments(nvar_fun, 20, fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        nevar, nevar_err = res
+        assert np.allclose(nevar, fvar**2, rtol=0.01)
+
+    def test_analyze_segments_nvar_fracstep_mean(self):
+        start, stop, mean = self.lc.analyze_segments(np.mean, 20, fraction_step=0.5)
+        start, stop, res = self.lc.analyze_segments(evar_fun, 20, fraction_step=0.5)
         # excess_variance returns fvar and fvar_err
         evar, evar_err = res
-        start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun, fraction_step=0.5)
+        start, stop, res = self.lc.analyze_segments(nvar_fun, 20, fraction_step=0.5)
         # excess_variance returns fvar and fvar_err
         nevar, nevar_err = res
         assert np.allclose(nevar * mean**2, evar, rtol=0.01)
@@ -426,13 +442,13 @@ class TestLightcurve(object):
         lc = Lightcurve(self.times, self.counts)
         assert lc.n == 4
 
-    def test_analyze_lc_chunks(self):
+    def test_analyze_segments(self):
         lc = Lightcurve(self.times, self.counts, gti=self.gti)
 
         def func(lc):
             return lc.time[0]
 
-        start, stop, res = lc.analyze_lc_chunks(2, func)
+        start, stop, res = lc.analyze_segments(func, 2)
         assert start[0] == 0.5
         assert np.allclose(start + lc.dt / 2, res)
 
@@ -556,7 +572,7 @@ class TestLightcurve(object):
         t = [0.5, 1.5, 2.5, 3.5, 4.5]
         lc = [5, 5, 0, 5, 5]
         gtis = [[0, 2], [3, 5]]
-        lc = Lightcurve(t, lc, gti=gtis)
+        lc = Lightcurve(t, lc, gti=gtis, dt=1)
 
         assert lc.meanrate == 5
         assert lc.meancounts == 5
@@ -1522,7 +1538,7 @@ class TestLightcurveRebin(object):
         lc = Lightcurve(time, count, gti=[[-0.5, 150.5]])
         lc.gti = [[-0.5, 2.5], [12.5, 14.5]]
         lc_new = lc.apply_gtis(inplace=inplace)
-        if inplace == True:
+        if inplace:
             assert lc_new is lc
         assert lc_new.n == 5
         for attr in lc_new.array_attrs():
@@ -1542,7 +1558,7 @@ class TestLightcurveRebin(object):
         lc_rate = Lightcurve(time, counts=countrate, input_counts=False, gti=[[-0.5, 10.5]])
         lc_rate.gti = [[-0.5, 2.5]]
         lc_rate_new = lc_rate.apply_gtis(inplace=inplace)
-        if inplace == True:
+        if inplace:
             assert lc_rate_new is lc_rate
         assert lc_rate_new.n == 2
         for attr in lc_rate_new.array_attrs():