Merge pull request #468 from StingraySoftware/allow_transparent_mjdref

Allow transparent MJDREF change

stingray/events.py

@@ -297,6 +297,10 @@ def join(self, other):
             simon("One of the event lists you are concatenating is empty.")
             other.time = np.asarray([])
 
+        # Tolerance for MJDREF:1 microsecond
+        if not np.isclose(self.mjdref, other.mjdref, atol=1e-6 / 86400):
+            other = other.change_mjdref(self.mjdref)
+
         ev_new.time = np.concatenate([self.time, other.time])
         order = np.argsort(ev_new.time)
         ev_new.time = ev_new.time[order]
@@ -511,3 +515,46 @@ def apply_deadtime(self, deadtime, inplace=False, **kwargs):
             new_ev = [new_ev, retall]
 
         return new_ev
+
+    def change_mjdref(self, new_mjdref):
+        """Change the MJD reference time (MJDREF) of the light curve.
+
+        Times will be now referred to this new MJDREF
+
+        Parameters
+        ----------
+        new_mjdref : float
+            New MJDREF
+
+        Returns
+        -------
+        new_lc : :class:`EventList` object
+            The new LC shifted by MJDREF
+        """
+        time_shift = (self.mjdref - new_mjdref) * 86400
+
+        new_ev = self.shift(time_shift)
+        new_ev.mjdref = new_mjdref
+        return new_ev
+
+    def shift(self, time_shift):
+        """
+        Shift the events and the GTIs in time.
+
+        Parameters
+        ----------
+        time_shift: float
+            The time interval by which the light curve will be shifted (in
+            the same units as the time array in :class:`Lightcurve`
+
+        Returns
+        -------
+        new_ev : lightcurve.Lightcurve object
+            The new event list shifted by ``time_shift``
+
+        """
+        new_ev = copy.deepcopy(self)
+        new_ev.time = new_ev.time + time_shift
+        new_ev.gti = new_ev.gti + time_shift
+
+        return new_ev

stingray/gti.py

@@ -620,6 +620,28 @@ def gti_len(gti):
     return np.sum(gti[:, 1] - gti[:, 0])
 
 
+@jit(nopython=True)
+def _check_separate(gti0, gti1):
+    """Numba-compiled core of ``check_separate``."""
+    gti0_start = gti0[:, 0]
+    gti0_end = gti0[:, 1]
+    gti1_start = gti1[:, 0]
+    gti1_end = gti1[:, 1]
+
+    if (gti0_end[-1] <= gti1_start[0]) or (gti1_end[-1] <= gti0_start[0]):
+        return True
+
+    for g in gti1.flatten():
+        for g0, g1 in zip(gti0[:, 0], gti0[:, 1]):
+            if (g <= g1) and (g >= g0):
+                return False
+    for g in gti0.flatten():
+        for g0, g1 in zip(gti1[:, 0], gti1[:, 1]):
+            if (g <= g1) and (g >= g0):
+                return False
+    return True
+
+
 def check_separate(gti0, gti1):
     """
     Check if two GTIs do not overlap.
@@ -636,6 +658,33 @@ def check_separate(gti0, gti1):
     -------
     separate: bool
         ``True`` if GTIs are mutually exclusive, ``False`` if not
+
+    Examples
+    --------
+    >>> gti0 = [[0, 10]]
+    >>> gti1 = [[20, 30]]
+    >>> check_separate(gti0, gti1)
+    True
+    >>> gti0 = [[0, 10]]
+    >>> gti1 = [[0, 10]]
+    >>> check_separate(gti0, gti1)
+    False
+    >>> gti0 = [[0, 10]]
+    >>> gti1 = [[10, 20]]
+    >>> check_separate(gti0, gti1)
+    True
+    >>> gti0 = [[0, 11]]
+    >>> gti1 = [[10, 20]]
+    >>> check_separate(gti0, gti1)
+    False
+    >>> gti0 = [[0, 11]]
+    >>> gti1 = [[10, 20]]
+    >>> check_separate(gti1, gti0)
+    False
+    >>> gti0 = [[0, 10], [30, 40]]
+    >>> gti1 = [[11, 28]]
+    >>> check_separate(gti0, gti1)
+    True
     """
 
     gti0 = np.asarray(gti0)
@@ -646,16 +695,9 @@ def check_separate(gti0, gti1):
     # Check if independently GTIs are well behaved
     check_gtis(gti0)
     check_gtis(gti1)
-
-    gti0_start = gti0[:, 0][0]
-    gti0_end = gti0[:, 1][-1]
-    gti1_start = gti1[:, 0][0]
-    gti1_end = gti1[:, 1][-1]
-
-    if (gti0_end <= gti1_start) or (gti1_end <= gti0_start):
-        return True
-    else:
-        return False
+    t0 = min(gti0[0, 0], gti1[0, 0])
+    return _check_separate((gti0 - t0).astype(np.double),
+                           (gti1 - t0).astype(np.double))
 
 
 def join_equal_gti_boundaries(gti):
@@ -701,13 +743,15 @@ def append_gtis(gti0, gti1):
     --------
     >>> np.all(append_gtis([[0, 1]], [[2, 3]]) == [[0, 1], [2, 3]])
     True
+    >>> np.allclose(append_gtis([[0, 1], [4, 5]], [[2, 3]]),
+    ...             [[0, 1], [2, 3], [4, 5]])
+    True
     >>> np.all(append_gtis([[0, 1]], [[1, 3]]) == [[0, 3]])
     True
     """
 
     gti0 = np.asarray(gti0)
     gti1 = np.asarray(gti1)
-
     # Check if independently GTIs are well behaved.
     check_gtis(gti0)
     check_gtis(gti1)
@@ -717,9 +761,9 @@ def append_gtis(gti0, gti1):
         raise ValueError('In order to append, GTIs must be mutually'
                          'exclusive.')
 
-    new_gtis = np.sort(np.concatenate([gti0, gti1]))
-
-    return join_equal_gti_boundaries(new_gtis)
+    new_gtis = np.concatenate([gti0, gti1])
+    order = np.argsort(new_gtis[:, 0])
+    return join_equal_gti_boundaries(new_gtis[order])
 
 
 def join_gtis(gti0, gti1):

stingray/lightcurve.py

@@ -4,7 +4,7 @@
 :class::class:`Lightcurve` is used to create light curves out of photon counting data
 or to save existing light curves in a class that's easy to use.
 """
-
+import warnings
 import logging
 import numpy as np
 import stingray.io as io
@@ -435,7 +435,7 @@ def change_mjdref(self, new_mjdref):
         new_lc : lightcurve.Lightcurve object
             The new LC shifted by MJDREF
         """
-        time_shift = (new_mjdref - self.mjdref) * 86400
+        time_shift = -(new_mjdref - self.mjdref) * 86400
 
         new_lc = self.shift(time_shift)
         new_lc.mjdref = new_mjdref
@@ -485,7 +485,8 @@ def _operation_with_other_lc(self, other, operation):
             The new light curve calculated in ``operation``
         """
         if self.mjdref != other.mjdref:
-            raise ValueError("MJDref is different in the two light curves")
+            warnings.warn("MJDref is different in the two light curves")
+            other = other.change_mjdref(self.mjdref)
 
         common_gti = cross_two_gtis(self.gti, other.gti)
         mask_self = create_gti_mask(self.time, common_gti, dt=self.dt)
@@ -924,7 +925,8 @@ def join(self, other):
         array([ 300,  100,  400,  600, 1200,  800])
         """
         if self.mjdref != other.mjdref:
-            raise ValueError("MJDref is different in the two light curves")
+            warnings.warn("MJDref is different in the two light curves")
+            other = other.change_mjdref(self.mjdref)
 
         if self.dt != other.dt:
             utils.simon("The two light curves have different bin widths.")

stingray/tests/test_events.py

@@ -226,6 +226,24 @@ def test_overlapping_join(self):
                 np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
         assert (ev_new.gti == np.array([[5, 6]])).all()
 
+    def test_overlapping_join_change_mjdref(self):
+        """Join two non-overlapping event lists.
+        """
+        ev = EventList(time=[1, 1, 10, 6, 5],
+                       energy=[10, 6, 3, 11, 2], gti=[[1, 3],[5, 6]],
+                       mjdref=57001)
+        ev_other = EventList(time=np.asarray([5, 7, 6, 6, 10]) + 86400,
+                             energy=[2, 3, 8, 1, 2],
+                             gti=np.asarray([[5, 7],[8, 10]]) + 86400,
+                             mjdref=57000)
+        ev_new = ev.join(ev_other)
+
+        assert np.allclose(ev_new.time,
+                           np.array([1, 1, 5, 5, 6, 6, 6, 7, 10, 10]))
+        assert (ev_new.energy ==
+                np.array([10, 6, 2, 2, 11, 8, 1, 3, 3, 2])).all()
+        assert np.allclose(ev_new.gti, np.array([[5, 6]]))
+
     def test_io_with_ascii(self):
         ev = EventList(self.time)
         ev.write('ascii_ev.txt',format_='ascii')

stingray/tests/test_lightcurve.py

@@ -214,7 +214,7 @@ def setup_class(cls):
         cls.times = np.array([1, 2, 3, 4])
         cls.counts = np.array([2, 2, 2, 2])
         cls.dt = 1.0
-        cls.gti = [[0.5, 4.5]]
+        cls.gti = np.array([[0.5, 4.5]])
 
     def test_create(self):
         """
@@ -462,12 +462,6 @@ def test_add_with_same_gtis(self):
         lc = lc1 + lc2
         np.testing.assert_almost_equal(lc.gti, self.gti)
 
-    def test_add_with_different_mjdref(self):
-        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
-        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
-        with pytest.raises(ValueError):
-            lc1 + lc2
-
     def test_add_with_different_gtis(self):
         gti = [[0., 3.5]]
         lc1 = Lightcurve(self.times, self.counts, gti=self.gti)
@@ -515,12 +509,6 @@ def test_sub_with_different_err_dist(self):
             assert np.any(["ightcurves have different statistics"
                            in str(wi.message) for wi in w])
 
-    def test_sub_with_different_mjdref(self):
-        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
-        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
-        with pytest.raises(ValueError):
-            lc1 - lc2
-
     def test_subtraction(self):
         _counts = [3, 4, 5, 6]
 
@@ -597,10 +585,37 @@ def test_join_with_different_dt(self):
                            in str(wi.message) for wi in w])
 
     def test_join_with_different_mjdref(self):
-        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
+        shift = 86400.  # day
+        lc1 = Lightcurve(self.times + shift, self.counts, gti=self.gti + shift,
+                         mjdref=57000)
         lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
-        with pytest.raises(ValueError):
-            lc1.join(lc2)
+        newlc = lc1.join(lc2)
+        # The join operation *averages* the overlapping arrays
+        assert np.allclose(newlc.counts, lc1.counts)
+
+    def test_sum_with_different_mjdref(self):
+        shift = 86400.  # day
+        lc1 = Lightcurve(self.times + shift, self.counts, gti=self.gti + shift,
+                         mjdref=57000)
+        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
+        with pytest.warns(UserWarning) as record:
+            newlc = lc1 + lc2
+        assert np.any(["MJDref"
+                       in r.message.args[0] for r in record])
+
+        assert np.allclose(newlc.counts, lc1.counts * 2)
+
+    def test_subtract_with_different_mjdref(self):
+        shift = 86400.  # day
+        lc1 = Lightcurve(self.times + shift, self.counts, gti=self.gti + shift,
+                         mjdref=57000)
+        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
+        with pytest.warns(UserWarning) as record:
+            newlc = lc1 - lc2
+        assert np.any(["MJDref"
+                       in r.message.args[0] for r in record])
+
+        assert np.allclose(newlc.counts, 0)
 
     def test_join_disjoint_time_arrays(self):
         _times = [5, 6, 7, 8]