Repository: StingraySoftware/stingray
Pull request #241 (2017-08-11 22:37:46)

Fixed behaviour of sort
`sort()` currently sorts only by counts. This is great, but it's not what I expected for a time series (I expected it to sort by time). It also turns out there are some issues that appear when the time bins are out of order. 

This PR should fix those issues, and implements `sort()` as sorting by time bin, and `sort_counts()` to sort by counts.

dhuppenkothen: Added switch. Does that look better?
matteobachetti: Sure! At least a few operations on the full array are avoided.
But.. why not just making `unsorted` a key argument? False by default, and only True if the user knows what they do?
dhuppenkothen: Because it seems like we can check this automatically, so there is no strict requirement for user input. 
matteobachetti: This is cheap if the light curve is small, but if it's large l think we should try to cut down all operations performed on the full array
dhuppenkothen: My feeling is that the uncaught errors that an unsorted light curve will produce in some of the methods (quietly dropping data points etc) are more serious than the speedups by not doing this check.
I am assuming that people only ever read the documentation very cursory, and will therefore largely miss/ignore the presence of an optional keyword like that.
matteobachetti: @dhuppenkothen I'll leave to you the choice if merging it as it is or making the suggested modifications

Patch:
diff --git a/stingray/gti.py b/stingray/gti.py
index 3b9d784b..3fa3d340 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -225,7 +225,7 @@ def create_gti_mask_complete(time, gtis, safe_interval=0, min_length=0,
 
     dt = assign_value_if_none(dt,
                               np.zeros_like(time) +
-                              np.median(np.diff(time)))
+                              np.median(np.diff(np.sort(time)) / 2))
 
     mask = np.zeros(len(time), dtype=bool)
 
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index 0ccfc638..46effe4c 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -66,7 +66,6 @@ class Lightcurve(object):
         mjdref: float
             MJD reference (useful in most high-energy mission data)
 
-
         Attributes
         ----------
         time: numpy.ndarray
@@ -163,8 +162,19 @@ class Lightcurve(object):
 
         self.mjdref = mjdref
         self.time = np.asarray(time)
+        dt_array = np.diff(np.sort(self.time))
+        dt_array_unsorted = np.diff(self.time)
+        unsorted = np.any(dt_array_unsorted < 0)
+
         if dt is None:
-            self.dt = np.median(self.time[1:] - self.time[:-1])
+            if unsorted:
+                logging.warning("The light curve is unordered! This may cause "
+                                "unexpected behaviour in some methods! Use "
+                                "sort() to order the light curve in time and "
+                                "check that the time resolution `dt` is "
+                                "calculated correctly!")
+
+            self.dt = np.median(dt_array)
         else:
             self.dt = dt
 
@@ -173,8 +183,12 @@ class Lightcurve(object):
 
         self.err_dist = err_dist
 
-        self.tstart = self.time[0] - 0.5*self.dt
-        self.tseg = self.time[-1] - self.time[0] + self.dt
+        if unsorted:
+            self.tstart = np.min(self.time) - 0.5 * self.dt
+            self.tseg = np.max(self.time) - np.min(self.time) + self.dt
+        else:
+            self.tstart = self.time[0] - 0.5*self.dt
+            self.tseg = self.time[-1] - self.time[0] + self.dt
 
         self.gti = \
             np.asarray(assign_value_if_none(gti,
@@ -216,7 +230,6 @@ class Lightcurve(object):
                   "This could cause problems with Fourier transforms. "
                   "Please make the input time evenly sampled.")
 
-
     def change_mjdref(self, new_mjdref):
         """Change the MJDREF of the light curve.
 
@@ -277,8 +290,9 @@ class Lightcurve(object):
                   "We are setting the errors to zero to avoid complications.")
             new_counts_err = np.zeros_like(new_counts)
         elif self.err_dist.lower() in valid_statistics:
-                new_counts_err = np.sqrt(np.add(self.counts_err[mask_self]**2,
-                                                other.counts_err[mask_other]**2))
+                new_counts_err = \
+                    np.sqrt(np.add(self.counts_err[mask_self]**2,
+                                   other.counts_err[mask_other]**2))
             # More conditions can be implemented for other statistics
         else:
             raise StingrayError("Statistics not recognized."
@@ -436,7 +450,7 @@ class Lightcurve(object):
             raise IndexError("The index must be either an integer or a slice "
                              "object !")
 
-    def __eq__(self,other_lc):
+    def __eq__(self, other_lc):
         """
         Compares two :class:`Lightcurve` objects.
 
@@ -454,7 +468,8 @@ class Lightcurve(object):
         """
         if not isinstance(other_lc, Lightcurve):
             raise ValueError('Lightcurve can only be compared with a Lightcurve Object')
-        if (np.allclose(self.time, other_lc.time) and np.allclose(self.counts, other_lc.counts)):
+        if (np.allclose(self.time, other_lc.time) and
+                np.allclose(self.counts, other_lc.counts)):
             return True
         return False
 
@@ -605,7 +620,6 @@ class Lightcurve(object):
             raise ValueError("New time resolution must be larger than "
                              "old time resolution!")
 
-
         if self.gti is None:
 
             bin_time, bin_counts, bin_err, _ = \
@@ -615,7 +629,7 @@ class Lightcurve(object):
         else:
             bin_time, bin_counts, bin_err = [], [], []
             gti_new = []
-            for g  in self.gti:
+            for g in self.gti:
                 if g[1] - g[0] < dt_new:
                     continue
                 else:
@@ -683,7 +697,7 @@ class Lightcurve(object):
         if self.dt != other.dt:
             utils.simon("The two light curves have different bin widths.")
 
-        if( self.tstart < other.tstart ):
+        if(self.tstart < other.tstart):
             first_lc = self
             second_lc = other
         else:
@@ -703,7 +717,6 @@ class Lightcurve(object):
                 simon("Lightcurves have different statistics!"
                       "We are setting the errors to zero.")
 
-
             elif self.err_dist.lower() in valid_statistics:
                 valid_err = True
             # More conditions can be implemented for other statistics
@@ -712,7 +725,6 @@ class Lightcurve(object):
                                     " Please use one of these: "
                                     "{}".format(valid_statistics))
 
-
             from collections import Counter
             counts = Counter()
             counts_err = Counter()
@@ -723,10 +735,11 @@ class Lightcurve(object):
 
             for i, time in enumerate(second_lc.time):
 
-                if (counts.get(time) != None): #Common time
-
-                    counts[time] = (counts[time] + second_lc.counts[i]) / 2  #avg
-                    counts_err[time] = np.sqrt(( ((counts_err[time]**2) + (second_lc.counts_err[i] **2)) / 2))
+                if counts.get(time) is not None:  # Common time
+                    counts[time] = (counts[time] + second_lc.counts[i]) / 2
+                    counts_err[time] = \
+                        np.sqrt(((counts_err[time] **2 ) +
+                                 (second_lc.counts_err[i] ** 2)) / 2)
 
                 else:
                     counts[time] = second_lc.counts[i]
@@ -745,7 +758,8 @@ class Lightcurve(object):
 
             new_time = np.concatenate([first_lc.time, second_lc.time])
             new_counts = np.concatenate([first_lc.counts, second_lc.counts])
-            new_counts_err = np.concatenate([first_lc.counts_err, second_lc.counts_err])
+            new_counts_err = \
+                np.concatenate([first_lc.counts_err, second_lc.counts_err])
 
         new_time = np.asarray(new_time)
         new_counts = np.asarray(new_counts)
@@ -757,7 +771,6 @@ class Lightcurve(object):
 
         return lc_new
 
-
     def truncate(self, start=0, stop=None, method="index"):
         """
         Truncate a Lightcurve object from points on the time array.
@@ -847,7 +860,47 @@ class Lightcurve(object):
 
     def sort(self, reverse=False):
         """
-        Sort a Lightcurve object in accordance with its counts array.
+        Sort a Lightcurve object by time.
+
+        A Lightcurve can be sorted in either increasing or decreasing order
+        using this method. The time array gets sorted and the counts array is
+        changed accordingly.
+
+        Parameters
+        ----------
+        reverse : boolean, default False
+            If True then the object is sorted in reverse order.
+
+        Example
+        -------
+        >>> time = [2, 1, 3]
+        >>> count = [200, 100, 300]
+        >>> lc = Lightcurve(time, count)
+        >>> lc_new = lc.sort()
+        >>> lc_new.time
+        array([1, 2, 3])
+        >>> lc_new.counts
+        array([100, 200, 300])
+
+        Returns
+        -------
+        lc_new: :class:`Lightcurve` object
+            The :class:`Lightcurve` object with sorted time and counts
+            arrays.
+        """
+
+        new_time, new_counts, new_counts_err = \
+            zip(*sorted(zip(self.time, self.counts, self.counts_err),
+                        reverse=reverse))
+
+        new_lc = Lightcurve(new_time, new_counts, err=new_counts_err,
+                            gti=self.gti, dt=self.dt, mjdref=self.mjdref)
+
+        return new_lc
+
+    def sort_counts(self, reverse=False):
+        """
+        Sort a Lightcurve object by counts.
 
         A Lightcurve can be sorted in either increasing or decreasing order
         using this method. The counts array gets sorted and the time array is
@@ -863,24 +916,27 @@ class Lightcurve(object):
         >>> time = [1, 2, 3]
         >>> count = [200, 100, 300]
         >>> lc = Lightcurve(time, count)
-        >>> lc.sort()
-        >>> lc.counts
-        array([100, 200, 300])
-        >>> lc.time
+        >>> lc_new = lc.sort_counts()
+        >>> lc_new.time
         array([2, 1, 3])
+        >>> lc_new.counts
+        array([100, 200, 300])
 
         Returns
         -------
         lc_new: :class:`Lightcurve` object
-            The :class:`Lightcurve` object with truncated time and counts
+            The :class:`Lightcurve` object with sorted time and counts
             arrays.
         """
 
-        new_counts, new_time, new_counts_err = zip(*sorted(zip(self.counts, self.time, self.counts_err) , reverse=reverse))
+        new_counts, new_time, new_counts_err = \
+            zip(*sorted(zip(self.counts, self.time, self.counts_err),
+                        reverse=reverse))
 
-        self.time = np.asarray(new_time)
-        self.counts = np.asarray(new_counts)
-        self.counts_err = np.asarray(new_counts_err)
+        new_lc = Lightcurve(new_time, new_counts, err=new_counts_err,
+                            gti=self.gti, dt=self.dt, mjdref=self.mjdref)
+
+        return new_lc
 
     def estimate_chunk_length(self, min_total_counts=100, min_time_bins=100):
         """Choose a reasonable chunk length.
@@ -948,8 +1004,8 @@ class Lightcurve(object):
         Other parameters
         ----------------
         fraction_step : float
-            If the step is not a full chunk_length but less (e.g. a moving window),
-            this indicates the ratio between step step and `chunk_length` (e.g.
+            If the step is not a full chunk_length but less (e.g. a moving
+            window), this indicates the ratio between step step and `chunk_length` (e.g.
             0.5 means that the window shifts of half chunk_length)
         kwargs : keyword arguments
             These additional keyword arguments, if present, they will be passed


Test patch:
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index 7c10e6d1..09c84e02 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -564,22 +564,42 @@ class TestLightcurve(object):
         assert lc2.mjdref == lc.mjdref
 
     def test_sort(self):
+        _times = [2, 1, 3, 4]
+        _counts = [40, 10, 20, 5]
+        lc = Lightcurve(_times, _counts, mjdref=57000)
+        mjdref = lc.mjdref
+
+        lc_new = lc.sort()
+
+        assert np.all(lc_new.counts == np.array([10, 40, 20, 5]))
+        assert np.all(lc_new.time == np.array([1, 2, 3, 4]))
+        assert lc_new.mjdref == mjdref
+
+        lc_new = lc.sort(reverse=True)
+
+        assert np.all(lc_new.counts == np.array([5, 20, 40,  10]))
+        assert np.all(lc_new.time == np.array([4, 3, 2, 1]))
+        assert lc_new.mjdref == mjdref
+
+    def test_sort_counts(self):
         _times = [1, 2, 3, 4]
         _counts = [40, 10, 20, 5]
         lc = Lightcurve(_times, _counts, mjdref=57000)
         mjdref = lc.mjdref
 
-        lc.sort()
+        lc_new = lc.sort_counts()
+
+        assert np.all(lc_new.counts == np.array([5, 10, 20, 40]))
+        assert np.all(lc_new.time == np.array([4, 2, 3, 1]))
+        assert lc_new.mjdref == mjdref
+
+        lc_new = lc.sort_counts(reverse=True)
 
-        assert np.all(lc.counts == np.array([5, 10, 20, 40]))
-        assert np.all(lc.time == np.array([4, 2, 3, 1]))
-        assert lc.mjdref == mjdref
+        assert np.all(lc_new.counts == np.array([40, 20, 10,  5]))
+        assert np.all(lc_new.time == np.array([1, 3, 2, 4]))
+        assert lc_new.mjdref == mjdref
 
-        lc.sort(reverse=True)
 
-        assert np.all(lc.counts == np.array([40, 20, 10,  5]))
-        assert np.all(lc.time == np.array([1, 3, 2, 4]))
-        assert lc.mjdref == mjdref
 
     def test_sort_reverse(self):
         times = np.arange(1000)