Repository: StingraySoftware/stingray
Pull request #774 (2023-10-18 23:19:38)

Gti fixes
Fixes some problems in `cross_two_gtis` and `join_gtis `. Notably:

+ `cross_two_gtis([[1,2]], [[2, 3]])` would give `[[2,3]]` instead of an empty list
+ `join_gtis([[1,2], [3, 4]], [[2, 3]])` would give `[[1,3], [3,4]]` instead of `[[1,4]]`

I checked that both things now work

Also, I add a new function, `merge_gtis`, that can do any operation on multiple GTIs (append, join, cross), while also checking them for consistency. This is useful for operations such as `events.join` and `lightcurve.join` that duplicate a lot of code trying to be smart about what to do when merging different datasets. Option `gti_treatment="infer"` makes the same kind of choices once.
 


Patch:
diff --git a/docs/changes/774.bugfix.rst b/docs/changes/774.bugfix.rst
new file mode 100644
index 00000000..64179168
--- /dev/null
+++ b/docs/changes/774.bugfix.rst
@@ -0,0 +1,1 @@
+Various bugfixes in `gti.py`, and a new function to interpret the mix of multiple GTIs.
\ No newline at end of file
diff --git a/stingray/gti.py b/stingray/gti.py
index 6efdca40..7c9e97bf 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -674,6 +674,11 @@ def cross_two_gtis(gti0, gti1):
     >>> newgti = cross_two_gtis(gti1, gti2)
     >>> np.allclose(newgti, [[1, 4]])
     True
+    >>> gti1 = np.array([[1, 2]])
+    >>> gti2 = np.array([[2, 3]])
+    >>> newgti = cross_two_gtis(gti1, gti2)
+    >>> len(newgti)
+    0
     """
     gti0 = join_equal_gti_boundaries(np.asarray(gti0))
     gti1 = join_equal_gti_boundaries(np.asarray(gti1))
@@ -704,7 +709,11 @@ def cross_two_gtis(gti0, gti1):
     conc_end = conc_end[order]
     conc_tag = conc_tag[order]
 
-    last_end = conc_start[0] - 1
+    last_end = conc_start[0] - 1.0
+
+    # The maximum end must not be larger than the second last end!
+    max_end = conc_end[-2]
+
     final_gti = []
     for ie, e in enumerate(conc_end):
         # Is this ending in series 0 or 1?
@@ -733,6 +742,9 @@ def cross_two_gtis(gti0, gti1):
         cond1 = (gti_end[other_series] > s) * (gti_end[other_series] < e)
         cond2 = gti_end[other_series][so_pos] < s
         condition = np.any(np.logical_or(cond1, cond2))
+        if e > max_end:
+            condition = True
+        # Also, the last closed interval in the other series must be before e
         # Well, if none of the conditions at point 2 apply, then you can
         # create the new gti!
         if not condition:
@@ -786,6 +798,8 @@ def cross_gtis(gti_list):
 
     for gti in gti_list[1:]:
         gti0 = cross_two_gtis(gti0, gti)
+        if len(gti0) == 0:
+            return []
 
     return gti0
 
@@ -954,6 +968,86 @@ def join_equal_gti_boundaries(gti, threshold=0.0):
     return np.asarray(ng)
 
 
+def merge_gtis(gti_list, gti_treatment):
+    """Merge a list of GTIs using the specified method.
+
+    Invalid GTI lists (None or empty) are ignored.
+
+    Parameters
+    ----------
+    gti_list : list of 2-d float arrays
+        List of GTIs.
+
+    Other parameters
+    ----------------
+    gti_treatment : {"intersection", "union", "append", "infer", "none"}
+        Method to use to merge the GTIs. If "intersection", the GTIs are merged
+        using the intersection of the GTIs. If "union", the GTIs are merged
+        using the union of the GTIs. If "none", a single GTI with the minimum and
+        the maximum time stamps of all GTIs is returned. If "infer", the strategy
+        is decided based on the GTIs. If there are no overlaps, "union" is used,
+        otherwise "intersection" is used. If "append", the GTIs are simply appended
+        but they must be mutually exclusive.
+
+    Examples
+    --------
+    >>> gti1 = np.array([[1, 2], [3, 4], [5, 6]])
+    >>> gti2 = np.array([[1, 2]])
+    >>> gti3 = np.array([[4, 5]])
+    >>> gti = merge_gtis([gti1, gti2], "intersection")
+    >>> np.array_equal(gti, [[1, 2]])
+    True
+    >>> merge_gtis([gti1, gti2, gti3], "intersection") is None
+    True
+    >>> merge_gtis([gti2, gti3], "intersection") is None
+    True
+    >>> gti = merge_gtis([gti1, gti2], "infer")
+    >>> np.array_equal(gti, [[1, 2]])
+    True
+    >>> gti = merge_gtis([gti2, gti3], "infer")
+    >>> np.array_equal(gti, [[1, 2], [4, 5]])
+    True
+    """
+    all_gti_lists = []
+    global_min = np.inf
+    global_max = -np.inf
+    for gti in gti_list:
+        if gti is None or len(gti) == 0:
+            continue
+        all_gti_lists.append(gti)
+        global_min = min(global_min, np.min(gti))
+        global_max = max(global_max, np.max(gti))
+
+    if len(all_gti_lists) == 0:
+        return None
+
+    if gti_treatment == "none":
+        return np.asarray([[global_min, global_max]])
+
+    if len(all_gti_lists) == 1:
+        return all_gti_lists[0]
+
+    cross = cross_gtis(all_gti_lists)
+    if len(cross) == 0:
+        cross = None
+    if gti_treatment == "infer":
+        if cross is None:
+            gti_treatment = "union"
+        else:
+            gti_treatment = "intersection"
+
+    if gti_treatment == "intersection":
+        return cross
+
+    gti0 = all_gti_lists[0]
+    for gti in all_gti_lists[1:]:
+        if gti_treatment == "union":
+            gti0 = join_gtis(gti0, gti)
+        elif gti_treatment == "append":
+            gti0 = append_gtis(gti0, gti)
+    return gti0
+
+
 def append_gtis(gti0, gti1):
     """
     Union of two non-overlapping GTIs.
@@ -993,7 +1087,7 @@ def append_gtis(gti0, gti1):
 
     # Check if GTIs are mutually exclusive.
     if not check_separate(gti0, gti1):
-        raise ValueError("In order to append, GTIs must be mutually" "exclusive.")
+        raise ValueError("In order to append, GTIs must be mutually exclusive.")
 
     new_gtis = np.concatenate([gti0, gti1])
     order = np.argsort(new_gtis[:, 0])
@@ -1016,10 +1110,23 @@ def join_gtis(gti0, gti1):
 
     ::
 
-        (cumsum)   -1   -2         -1   0   -1 -2           -1  -2  -1        0
-        GTI A      |-----:----------|   :    |--:------------|   |---:--------|
-        FINAL GTI  |-----:--------------|    |--:--------------------:--------|
-        GTI B            |--------------|       |--------------------|
+        (g_all)    0     1     2     3     4     5     6     7     8     9
+        (cumsum)   -1   -2    -1     0    -1    -2     -1   -2    -1     0
+        GTI A      |-----:-----|     :     |-----:-----|     |-----:-----|
+        FINAL GTI  |-----:-----------|     |-----:-----------------:-----|
+        GTI B            |-----------|           |-----------------|
+
+    In case one GTI ends exactly where another one starts, the cumulative sum is 0
+    but we do not want to close. In this case, we make a check that the next element
+    of the sequence is not equal to the one where we would close.
+
+    ::
+
+        (g_all)    0    1,1         3,3          5
+        (cumsum)   -1   0,-1       -1,-2         0
+        GTI A      |-----|           |-----------|
+        FINAL GTI  |-----------------------------|
+        GTI B            |-----------|
 
     Parameters
     ----------
@@ -1047,9 +1154,13 @@ def join_gtis(gti0, gti1):
 
     g0 = gti0.flatten()
     # Opening GTI: type = 1; Closing: type = -1
-    g0_type = np.asarray(list(zip(-np.ones(int(len(g0) / 2)), np.ones(int(len(g0) / 2)))))
+    g0_type = np.asarray(
+        list(zip(-np.ones(int(len(g0) / 2), dtype=int), np.ones(int(len(g0) / 2), dtype=int)))
+    )
     g1 = gti1.flatten()
-    g1_type = np.asarray(list(zip(-np.ones(int(len(g1) / 2)), np.ones(int(len(g1) / 2)))))
+    g1_type = np.asarray(
+        list(zip(-np.ones(int(len(g1) / 2), dtype=int), np.ones(int(len(g1) / 2), dtype=int)))
+    )
 
     g_all = np.append(g0, g1)
     g_type_all = np.append(g0_type, g1_type)
@@ -1059,8 +1170,12 @@ def join_gtis(gti0, gti1):
 
     sums = np.cumsum(g_type_all)
 
-    # Where the cumulative sum is zero, we close the GTI
-    closing_bins = sums == 0
+    # Where the cumulative sum is zero, we close the GTI.
+    # But pay attention! If one GTI ends exactly where another one starts,
+    # the cumulative sum is zero, but we do not want to close the GTI.
+    # So we check that the next element of g_all is not equal to the one where
+    # we would close.
+    closing_bins = (sums == 0) & (g_all != np.roll(g_all, -1))
     # The next element in the sequence is the start of the new GTI. In the case
     # of the last element, the next is the first. Numpy.roll gives this for
     # free.


Test patch:
diff --git a/stingray/tests/test_gti.py b/stingray/tests/test_gti.py
index c4af5514..b3cf4647 100644
--- a/stingray/tests/test_gti.py
+++ b/stingray/tests/test_gti.py
@@ -3,7 +3,7 @@ import pytest
 import os
 
 from stingray.gti import cross_gtis, append_gtis, load_gtis, get_btis, join_gtis
-from stingray.gti import check_separate, create_gti_mask, check_gtis
+from stingray.gti import check_separate, create_gti_mask, check_gtis, merge_gtis
 from stingray.gti import create_gti_from_condition, gti_len, gti_border_bins
 from stingray.gti import time_intervals_from_gtis, bin_intervals_from_gtis
 from stingray.gti import create_gti_mask_complete, join_equal_gti_boundaries
@@ -40,6 +40,53 @@ class TestGTI(object):
 
         assert np.allclose(newgti, gti1), "GTIs do not coincide!"
 
+    def test_crossgti4(self):
+        """A more complicated example of intersection of GTIs."""
+        gti1 = np.array([[2, 3]])
+        gti2 = np.array([[3, 4]])
+        newgti = cross_gtis([gti1, gti2])
+        gti3 = np.array([[3, 5]])
+        assert len(newgti) == 0
+
+        newgti = cross_gtis([gti1, gti2, gti3])
+        assert len(newgti) == 0
+
+    def test_crossgti5(self):
+        """A more complicated example of intersection of GTIs."""
+        gti1 = np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]])
+        gti2 = np.array([[0.5, 14]])
+        newgti0 = cross_gtis([gti1, gti2])
+        newgti1 = cross_gtis([gti2, gti1])
+
+        assert np.allclose(gti1, np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]]))
+        assert np.allclose(gti2, np.array([[0.5, 14]]))
+        for newgti in [newgti0, newgti1]:
+            assert np.allclose(newgti, gti1)
+
+    def test_crossgti6(self):
+        """A more complicated example of intersection of GTIs."""
+        gti1 = np.array([[1.5, 12.5]])
+        gti2 = np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]])
+        newgti0 = cross_gtis([gti1, gti2])
+        newgti1 = cross_gtis([gti2, gti1])
+
+        for newgti in [newgti0, newgti1]:
+            assert np.allclose(
+                newgti, np.array([[1.5, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 12.5]])
+            )
+
+    def test_crossgti7(self):
+        """A more complicated example of intersection of GTIs."""
+        gti1 = np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]])
+        gti2 = np.array([[0.5, 3], [4.5, 4.7], [10, 14]])
+        newgti0 = cross_gtis([gti1, gti2])
+        newgti1 = cross_gtis([gti2, gti1])
+
+        assert np.allclose(gti1, np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]]))
+        assert np.allclose(gti2, np.array([[0.5, 3], [4.5, 4.7], [10, 14]]))
+        for newgti in [newgti0, newgti1]:
+            assert np.allclose(newgti, np.array([[1, 2], [4.5, 4.7], [11, 11.2], [12.2, 13.2]]))
+
     def test_bti(self):
         """Test the inversion of GTIs."""
         gti = np.array([[1, 2], [4, 5], [7, 10], [11, 11.2], [12.2, 13.2]])
@@ -204,6 +251,11 @@ class TestGTI(object):
             join_gtis(gti0, gti1) == np.array([[0, 1], [2, 3], [4, 8], [10, 11], [12, 13]])
         )
 
+    def test_join_gtis_in_middle(self):
+        gti0 = [[0, 1], [2, 3], [4, 8]]
+        gti1 = [[1, 2], [3, 4]]
+        assert np.all(join_gtis(gti0, gti1) == np.array([[0, 8]]))
+
     def test_time_intervals_from_gtis(self):
         """Test the division of start and end times to calculate spectra."""
         start_times, stop_times = time_intervals_from_gtis(
@@ -311,3 +363,54 @@ class TestGTI(object):
         )
         newg = join_equal_gti_boundaries(gti)
         assert np.allclose(newg, np.array([[1.16703354e08, 1.16703514e08]]))
+
+
+_ALL_METHODS = ["intersection", "union", "infer", "append"]
+
+
+class TestMergeGTIs(object):
+    @classmethod
+    def setup_class(cls):
+        cls.gti1 = np.array([[1, 2], [3, 4], [5, 6]])
+        cls.gti2 = np.array([[1, 2]])
+        cls.gti3 = np.array([[2, 3]])
+        cls.gti4 = np.array([[4, 5]])
+
+    @pytest.mark.parametrize("method", _ALL_METHODS)
+    def test_merge_gti_empty(self, method):
+        assert merge_gtis([], method) is None
+        assert merge_gtis([None], method) is None
+        assert merge_gtis([None, []], method) is None
+        assert merge_gtis([[]], method) is None
+
+    @pytest.mark.parametrize("method", _ALL_METHODS)
+    def test_merge_gti_single(self, method):
+        # all methods but `none` should just return the unaltered GTI
+        assert np.array_equal(merge_gtis([self.gti1], method), self.gti1)
+
+    def test_merge_gti_none(self):
+        assert np.array_equal(merge_gtis([self.gti1], "none"), [[1, 6]])
+        assert np.array_equal(merge_gtis([self.gti1, self.gti2], "none"), [[1, 6]])
+
+    def test_merge_gti_intersection(self):
+        gti = merge_gtis([self.gti1, self.gti2], "intersection")
+        assert np.array_equal(gti, [[1, 2]])
+        assert merge_gtis([self.gti1, self.gti2, self.gti3], "intersection") is None
+        assert merge_gtis([self.gti2, self.gti3], "intersection") is None
+
+    def test_merge_gti_union(self):
+        assert np.array_equal(merge_gtis([self.gti1, self.gti2], "union"), self.gti1)
+        assert np.array_equal(merge_gtis([self.gti1, self.gti3], "union"), [[1, 4], [5, 6]])
+
+    def test_merge_gti_append(self):
+        assert np.array_equal(merge_gtis([self.gti2, self.gti4], "append"), [[1, 2], [4, 5]])
+        assert np.array_equal(merge_gtis([self.gti2, self.gti3], "append"), [[1, 3]])
+        with pytest.raises(ValueError, match="In order to append, GTIs must be mutually"):
+            merge_gtis([self.gti1, self.gti3], "append")
+
+    def test_merge_gti_infer(self):
+        gti = merge_gtis([self.gti1, self.gti2], "infer")
+        assert np.array_equal(gti, [[1, 2]])
+
+        gti = merge_gtis([self.gti2, self.gti4], "infer")
+        assert np.array_equal(gti, [[1, 2], [4, 5]])