diff --git a/stingray/events.py b/stingray/events.py index 8c998d4a5..e5df36595 100644 --- a/stingray/events.py +++ b/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 diff --git a/stingray/gti.py b/stingray/gti.py index 90d6de0f8..dfedd4ef9 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 9ebbca56b..238d456ef 100644 --- a/stingray/lightcurve.py +++ b/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.") diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py index a7e42ca60..0b84f8701 100644 --- a/stingray/tests/test_events.py +++ b/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') diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py index 40b7e62b5..52bfde72a 100644 --- a/stingray/tests/test_lightcurve.py +++ b/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]