Repository: StingraySoftware/stingray

Issue:
Allow merging event lists from different instruments
Need to have some `change_mjdref` method in `EventList`, similar to the one in `Lightcurve`. 

Discussion:
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/tests/test_events.py`](https://github.com/StingraySoftware/stingray/blob/6e393fd2a564269b902491ac81e1c96ab5e35f47/stingray/tests/test_events.py):

> [Line 233:61](https://github.com/StingraySoftware/stingray/blob/6e393fd2a564269b902491ac81e1c96ab5e35f47/stingray/tests/test_events.py#L233): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','
> [Line 237:52](https://github.com/StingraySoftware/stingray/blob/6e393fd2a564269b902491ac81e1c96ab5e35f47/stingray/tests/test_events.py#L237): [E231](https://duckduckgo.com/?q=pep8%20E231) missing whitespace after ','




matteobachetti: @dhuppenkothen I rebased to the current master
codecov-commenter: # [Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=h1) Report
> Merging [#468](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=desc) into [master](https://codecov.io/gh/StingraySoftware/stingray/commit/870f122354d1763b4542640c75d2280bfc2a10c3&el=desc) will **increase** coverage by `0.04%`.
> The diff coverage is `97.50%`.

[![Impacted file tree graph](https://codecov.io/gh/StingraySoftware/stingray/pull/468/graphs/tree.svg?width=650&height=150&src=pr&token=FjWeFfhU9F)](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master     #468      +/-   ##
==========================================
+ Coverage   95.17%   95.21%   +0.04%     
==========================================
  Files          35       35              
  Lines        5176     5204      +28     
==========================================
+ Hits         4926     4955      +29     
+ Misses        250      249       -1     
```


| [Impacted Files](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [stingray/gti.py](https://codecov.io/gh/StingraySoftware/stingray/pull/468/diff?src=pr&el=tree#diff-c3RpbmdyYXkvZ3RpLnB5) | `97.89% <95.45%> (-0.23%)` | :arrow_down: |
| [stingray/events.py](https://codecov.io/gh/StingraySoftware/stingray/pull/468/diff?src=pr&el=tree#diff-c3RpbmdyYXkvZXZlbnRzLnB5) | `100.00% <100.00%> (ø)` | |
| [stingray/lightcurve.py](https://codecov.io/gh/StingraySoftware/stingray/pull/468/diff?src=pr&el=tree#diff-c3RpbmdyYXkvbGlnaHRjdXJ2ZS5weQ==) | `97.02% <100.00%> (+0.39%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=footer). Last update [870f122...f2794c9](https://codecov.io/gh/StingraySoftware/stingray/pull/468?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).


Patch:
diff --git a/stingray/events.py b/stingray/events.py
index 8c998d4a..e5df3659 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -297,6 +297,10 @@ class EventList(object):
             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 @@ class EventList(object):
             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
diff --git a/stingray/gti.py b/stingray/gti.py
index 90d6de0f..dfedd4ef 100644
--- a/stingray/gti.py
+++ b/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):
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 9ebbca56..238d456e 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -4,7 +4,7 @@ Definition of :class::class:`Lightcurve`.
 :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 @@ class Lightcurve(object):
         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 @@ class Lightcurve(object):
             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 @@ class Lightcurve(object):
         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.")


Test patch:
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index a7e42ca6..0b84f870 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -226,6 +226,24 @@ class TestEvents(object):
                 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')
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 40b7e62b..52bfde72 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -214,7 +214,7 @@ class TestLightcurve(object):
         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 @@ class TestLightcurve(object):
         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 @@ class TestLightcurve(object):
             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 @@ class TestLightcurve(object):
                            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]