Repository: StingraySoftware/stingray

Issue:
Improve performance of Lightcurve for astropy.time.Time objects
Right now, it looks like we *can* use `Lightcurve.make_lightcurve` on `astropy.time.Time` objects, but for reasons I don't understand that's vastly slower than sticking in a numpy array. It would be useful to check why that is, and speed up that operation.

Discussion:

Patch:
diff --git a/stingray/events.py b/stingray/events.py
index d8ebf954..af92100e 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -15,7 +15,7 @@ from .filters import get_deadtime_mask
 from .gti import append_gtis, check_separate, cross_gtis
 from .io import read, write
 from .lightcurve import Lightcurve
-from .utils import assign_value_if_none, simon
+from .utils import assign_value_if_none, simon, interpret_times
 
 __all__ = ['EventList']
 
@@ -97,6 +97,7 @@ class EventList(object):
         self.ncounts = ncounts
 
         if time is not None:
+            time, mjdref = interpret_times(time, mjdref)
             if not high_precision:
                 self.time = np.asarray(time)
             else:
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 196e7e2b..9cc9aa1a 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -15,7 +15,8 @@ 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 (assign_value_if_none, baseline_als,
-                            poisson_symmetrical_errors, simon)
+                            poisson_symmetrical_errors, simon, interpret_times)
+
 
 __all__ = ["Lightcurve"]
 
@@ -138,8 +139,8 @@ class Lightcurve(object):
                  gti=None, err_dist='poisson', mjdref=0, dt=None,
                  skip_checks=False, low_memory=False):
 
-        if hasattr(time, 'value'):
-            time = time.value
+        time, mjdref = interpret_times(time, mjdref=mjdref)
+
         time = np.asarray(time)
         counts = np.asarray(counts)
         if err is not None:
@@ -793,6 +794,7 @@ class Lightcurve(object):
         lc: :class:`Lightcurve` object
             A :class:`Lightcurve` object with the binned light curve
         """
+        toa, mjdref = interpret_times(toa, mjdref=mjdref)
 
         toa = np.sort(np.asarray(toa))
         # tstart is an optional parameter to set a starting time for
diff --git a/stingray/utils.py b/stingray/utils.py
index 8565d60d..cfd18a9d 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -8,6 +8,9 @@ from collections.abc import Iterable
 
 import numpy as np
 import scipy
+from astropy.time import Time, TimeDelta
+import astropy.units as u
+from astropy.units import Quantity
 
 # If numba is installed, import jit. Otherwise, define an empty decorator with
 # the same name.
@@ -130,7 +133,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
     Parameters
     ----------
     x: iterable
-        The dependent variable with some resolution, which can vary throughout 
+        The dependent variable with some resolution, which can vary throughout
         the time series.
 
     y: iterable
@@ -149,7 +152,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         each new bin of ``x``, or take the arithmetic mean.
 
     dx: float
-        The old resolution (otherwise, calculated from difference between 
+        The old resolution (otherwise, calculated from difference between
         time bins)
 
     Returns
@@ -200,7 +203,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
                          "old frequency resolution.")
 
     # left and right bin edges
-    # assumes that the points given in `x` correspond to 
+    # assumes that the points given in `x` correspond to
     # the left bin edges
     xedges = np.hstack([x, x[-1]+dx_old[-1]])
 
@@ -217,7 +220,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         xmax = xbin[i+1]
         min_ind = xedges.searchsorted(xmin)
         max_ind = xedges.searchsorted(xmax)
-        
+
         output[i] = np.sum(y[min_ind:max_ind-1])
         outputerr[i] = np.sum(yerr[min_ind:max_ind-1])
         step_size[i] = len(y[min_ind:max_ind-1])
@@ -234,7 +237,7 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
             output[i] += y[max_ind-1]*post_frac
             outputerr[i] += yerr[max_ind-1]*post_frac
             step_size[i] += post_frac
-        
+
     if method in ['mean', 'avg', 'average', 'arithmetic mean']:
         ybin = output / step_size
         ybinerr = np.sqrt(outputerr) / step_size
@@ -253,11 +256,11 @@ def rebin_data(x, y, dx_new, yerr=None, method='sum', dx=None):
         ybin = ybin[:-1]
         ybinerr = ybinerr[:-1]
         step_size = step_size[:-1]
-        
+
     dx_var = np.var(dx_old) / np.mean(dx_old)
 
     if np.size(dx_old) == 1 or dx_var < 1e-6:
-        step_size = step_size[0] 
+        step_size = step_size[0]
 
     new_x0 = (x[0] - (0.5 * dx_old[0])) + (0.5 * dx_new)
     xbin = np.arange(ybin.shape[0]) * dx_new + new_x0
@@ -1054,3 +1057,80 @@ def genDataPath(dir_path):
 
     else:
         raise IOError(("Directory does not exist."))
+
+
+def interpret_times(time, mjdref=0):
+    """Get time interval in seconds from an astropy Time object
+
+    Examples
+    --------
+    >>> time = Time(57483, format='mjd')
+    >>> newt, mjdref = interpret_times(time)
+    >>> newt == 0
+    True
+    >>> mjdref == 57483
+    True
+    >>> time = Time([57483], format='mjd')
+    >>> newt, mjdref = interpret_times(time)
+    >>> np.all(newt == 0)
+    True
+    >>> mjdref == 57483
+    True
+    >>> time = TimeDelta([3, 4, 5] * u.s)
+    >>> newt, mjdref = interpret_times(time)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> time = np.array([3, 4, 5])
+    >>> newt, mjdref = interpret_times(time, mjdref=45000)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> mjdref == 45000
+    True
+    >>> time = np.array([3, 4, 5] * u.s)
+    >>> newt, mjdref = interpret_times(time, mjdref=45000)
+    >>> np.allclose(newt, [3, 4, 5])
+    True
+    >>> mjdref == 45000
+    True
+    >>> newt, mjdref = interpret_times(1, mjdref=45000)
+    >>> newt == 1
+    True
+    >>> newt, mjdref = interpret_times(list, mjdref=45000)
+    Traceback (most recent call last):
+    ...
+    ValueError: Unknown time format: ...
+    >>> newt, mjdref = interpret_times("guadfkljfd", mjdref=45000)
+    Traceback (most recent call last):
+    ...
+    ValueError: Unknown time format: ...
+    """
+    if isinstance(time, TimeDelta):
+        out_times = time.to('s').value
+        return out_times, mjdref
+
+    if isinstance(time, Time):
+        mjds = time.mjd
+        if mjdref == 0:
+            if isinstance(mjds, Iterable):
+                mjdref = mjds[0]
+            else:
+                mjdref = mjds
+
+        out_times = (mjds - mjdref) * 86400
+        return out_times, mjdref
+
+    if isinstance(time, Quantity):
+        out_times = time.to('s').value
+        return out_times, mjdref
+
+    if isinstance(time, (tuple, list, np.ndarray)):
+        return time, mjdref
+
+    if not isinstance(time, Iterable):
+        try:
+            float(time)
+            return time, mjdref
+        except (ValueError, TypeError):
+            pass
+
+    raise ValueError(f"Unknown time format: {type(time)}")


Test patch:
diff --git a/stingray/tests/test_events.py b/stingray/tests/test_events.py
index 846d4135..1b1280ab 100644
--- a/stingray/tests/test_events.py
+++ b/stingray/tests/test_events.py
@@ -2,6 +2,7 @@
 import numpy as np
 import os
 import pytest
+from astropy.time import Time
 
 from ..events import EventList
 from ..lightcurve import Lightcurve
@@ -28,6 +29,22 @@ class TestEvents(object):
                          [1000, 2040, 1000, 3000, 4020, 2070]]
         self.gti = np.asarray([[0, 4]])
 
+    def test_initiate_from_ndarray(self):
+        times = np.sort(
+            np.random.uniform(1e8, 1e8 + 1000, 101).astype(np.longdouble))
+        ev = EventList(times, mjdref=54600)
+        assert np.allclose(ev.time, times, atol=1e-15)
+        assert np.allclose(ev.mjdref, 54600)
+
+    def test_initiate_from_astropy_time(self):
+        times = np.sort(
+            np.random.uniform(1e8, 1e8 + 1000, 101).astype(np.longdouble))
+        mjdref = 54600
+        mjds = Time(mjdref + times / 86400, format='mjd')
+        ev = EventList(mjds, mjdref=mjdref)
+        assert np.allclose(ev.time, times, atol=1e-15)
+        assert np.allclose(ev.mjdref, mjdref)
+
     def test_create_high_precision_object(self):
         times = np.sort(
             np.random.uniform(1e8, 1e8 + 1000, 101).astype(np.longdouble))
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 24a7be3d..ce8e8384 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -67,19 +67,23 @@ class TestProperties(object):
         _ = lc.bin_lo
         assert lc._bin_lo is not None
 
+    def test_lightcurve_from_astropy_time(self):
+        time = Time([57483, 57484], format='mjd')
+        counts = np.array([2, 2])
+        lc = Lightcurve(time, counts)
+        assert lc.dt == 86400
+        assert np.all(lc.counts == counts)
+
     def test_time_is_quantity_or_astropy_time(self):
         counts = [34, 21.425]
-        times = [57000, 58000]
+        times = np.array([57000, 58000])
 
-        times_q = times * u.d
+        times_q = (times - times[0]) * u.d
         times_t = Time(times, format='mjd')
 
-        lc = Lightcurve(time=times, counts=counts)
-        lc_q = Lightcurve(time=times_q, counts=counts)
+        lc_q = Lightcurve(time=times_q, counts=counts, mjdref=times[0])
         lc_t = Lightcurve(time=times_t, counts=counts)
-        assert_allclose(lc.time, lc_q.time)
-        assert_allclose(lc.time, lc_t.time)
-
+        assert_allclose(lc_q.time, lc_t.time)
 
     def test_gti(self):
         lc = copy.deepcopy(self.lc)
@@ -321,6 +325,25 @@ class TestLightcurve(object):
         assert np.allclose(lc.time, lc2.time)
         assert np.all(lc.counts == lc2.counts)
 
+    def test_lightcurve_from_toa_quantity(self):
+        lc = Lightcurve.make_lightcurve(self.times * u.s, self.dt,
+                                        use_hist=True, tstart=0.5)
+        lc2 = Lightcurve.make_lightcurve(self.times, self.dt, use_hist=False,
+                                         tstart=0.5)
+        assert np.allclose(lc.time, lc2.time)
+        assert np.all(lc.counts == lc2.counts)
+
+    def test_lightcurve_from_toa_Time(self):
+        mjdref = 56789
+        mjds = Time(self.times / 86400 + mjdref, format='mjd')
+
+        lc = Lightcurve.make_lightcurve(mjds, self.dt, mjdref=mjdref,
+                                        use_hist=True, tstart=0.5)
+        lc2 = Lightcurve.make_lightcurve(self.times, self.dt, use_hist=False,
+                                         tstart=0.5, mjdref=mjdref)
+        assert np.allclose(lc.time, lc2.time)
+        assert np.all(lc.counts == lc2.counts)
+
     def test_lightcurve_from_toa_halfbin(self):
         lc = Lightcurve.make_lightcurve(self.times + 0.5, self.dt,
                                         use_hist=True,