Repository: StingraySoftware/stingray
Pull request #253 (2017-09-04 19:44:29)

Fix excess variance error bar; allow normalized excess variance; fix GTI behavior during slicing
This has been triggered by Issue 282 in DAVE. That's probably not an issue, but I did find a mistake in the error bar calculation for the unnormalized excess variance. I took the chance to document the function better and to allow both normalized and unnormalized excess variance.

ADD: During tests of this new functionality, I discovered a bug in the slice operation with light curves, where GTIs were forgotten. I fixed this and tested it properly.


Patch:
diff --git a/.travis.yml b/.travis.yml
index 00626f2f..136836be 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -72,14 +72,15 @@ matrix:
         - python: 3.4
           env: NUMPY_VERSION=1.11
 
+        - python: 3.6
+          env: SETUP_CMD='build_docs -w'
+
         # Try older scipy versions
         - python: 2.7
           env: SCIPY_VERSION=0.16.0
         - python: 2.7
           env: SCIPY_VERSION=0.14.0
 
-        - python: 2.7
-          env: SETUP_CMD='build_sphinx -w'
         # Try Astropy development version
         - python: 2.7
           env: ASTROPY_VERSION=development
@@ -93,8 +94,6 @@ matrix:
         - python: 2.7
           env: SCIPY_VERSION=0.14.0
 
-        - python: 2.7
-          env: SETUP_CMD='build_sphinx -w'
         # Try Astropy development version
         - python: 2.7
           env: ASTROPY_VERSION=development
diff --git a/appveyor.yml b/appveyor.yml
index e6f20dfe..194171d2 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -10,8 +10,8 @@ environment:
         PYTHON_ARCH: "64" # needs to be set for CMD_IN_ENV to succeed. If a mix
                           # of 32 bit and 64 bit builds are needed, move this
                           # to the matrix section.
-        CONDA_DEPENDENCIES: "scipy numpy nose h5py astropy matplotlib"
-        PIP_DEPENDENCIES: "emcee statsmodels pytest-catchlog"
+        CONDA_DEPENDENCIES: "scipy numpy nose h5py astropy matplotlib statsmodels"
+        PIP_DEPENDENCIES: "emcee pytest-catchlog"
 
     matrix:
 
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index a4f1e517..d9ee3533 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -43,8 +43,8 @@ class Lightcurve(object):
             the bins defined in `time` (note: use `input_counts=False` to
             input the count rage, i.e. counts/second, otherwise use
             counts/bin). If None, we assume the data is poisson distributed
-            and calculate the error from the average of the lower and upper 
-            1-sigma confidence intervals for the Poissonian distribution with 
+            and calculate the error from the average of the lower and upper
+            1-sigma confidence intervals for the Poissonian distribution with
             mean equal to `counts`.
 
         input_counts: bool, optional, default True
@@ -388,7 +388,11 @@ class Lightcurve(object):
 
         This method adds functionality to retrieve the count value at
         a particular index. This also can be used for slicing and generating
-        a new Lightcurve object.
+        a new Lightcurve object. GTIs are recalculated based on the new light
+        curve segment
+
+        If the slice object is of kind start:stop:step, GTIs are also sliced,
+        and rewritten as zip(time - self.dt /2, time + self.dt / 2)
 
         Parameters
         ----------
@@ -408,9 +412,23 @@ class Lightcurve(object):
         if isinstance(index, int):
             return self.counts[index]
         elif isinstance(index, slice):
-            new_counts = self.counts[index.start:index.stop:index.step]
-            new_time = self.time[index.start:index.stop:index.step]
-            return Lightcurve(new_time, new_counts, mjdref=self.mjdref)
+            start = assign_value_if_none(index.start, 0)
+            stop = assign_value_if_none(index.stop, len(self.counts))
+            step = assign_value_if_none(index.step, 1)
+
+            new_counts = self.counts[start:stop:step]
+            new_time = self.time[start:stop:step]
+            new_gti = [[self.time[start] - 0.5 * self.dt,
+                        self.time[stop - 1] + 0.5 * self.dt]]
+            new_gti = np.asarray(new_gti)
+            if step > 1:
+                new_gt1 = np.array(list(zip(new_time - self.dt / 2,
+                                            new_time + self.dt / 2)))
+                new_gti = cross_two_gtis(new_gti, new_gt1)
+
+            new_gti = cross_two_gtis(self.gti, new_gti)
+            return Lightcurve(new_time, new_counts, mjdref=self.mjdref,
+                              gti=new_gti)
         else:
             raise IndexError("The index must be either an integer or a slice "
                              "object !")
@@ -424,7 +442,7 @@ class Lightcurve(object):
             "smoothness" parameter. Larger values make the baseline stiffer
             Typically 1e2 < lam < 1e9
         p : float
-            "asymmetry" parameter. Smaller values make the baseline more 
+            "asymmetry" parameter. Smaller values make the baseline more
             "horizontal". Typically 0.001 < p < 0.1, but not necessary.
         """
         baseline = np.zeros_like(self.time)
@@ -619,7 +637,7 @@ class Lightcurve(object):
 
             elif self.err_dist.lower() in valid_statistics:
                 valid_err = True
-            # More conditions can be implemented for other statistics 
+            # More conditions can be implemented for other statistics
             else:
                 raise StingrayError("Statistics not recognized."
                                     " Please use one of these: "
@@ -635,7 +653,7 @@ class Lightcurve(object):
                 counts_err[time] = first_lc.counts_err[i]
 
             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
@@ -651,7 +669,7 @@ class Lightcurve(object):
                 new_counts_err = list(counts_err.values())
             else:
                 new_counts_err = np.zeros_like(new_counts)
-            
+
             del[counts, counts_err]
 
         else:
diff --git a/stingray/utils.py b/stingray/utils.py
index 8eaebc1d..0ce02728 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -32,6 +32,12 @@ except ImportError:
             return wrapped_f
 
 
+__all__ = ['simon', 'rebin_data', 'rebin_data_log', 'look_for_array_in_array',
+           'is_string', 'is_iterable', 'order_list_of_arrays',
+           'optimal_bin_time', 'contiguous_regions', 'is_int',
+           'get_random_state', 'baseline_als', 'excess_variance',
+           'create_window']
+
 def _root_squared_mean(array):
     return np.sqrt(np.sum(array ** 2)) / len(array)
 
@@ -356,7 +362,7 @@ def baseline_als(y, lam, p, niter=10):
         "smoothness" parameter. Larger values make the baseline stiffer
         Typically 1e2 < lam < 1e9
     p : float
-        "asymmetry" parameter. Smaller values make the baseline more 
+        "asymmetry" parameter. Smaller values make the baseline more
         "horizontal". Typically 0.001 < p < 0.1, but not necessary.
     """
     from scipy import sparse
@@ -374,15 +380,31 @@ def baseline_als(y, lam, p, niter=10):
 def excess_variance(lc, normalization='fvar'):
     """Calculate the excess variance.
 
-    Vaughan+03
+    Vaughan et al. 2003, MNRAS 345, 1271 give three measurements of source
+    intrinsic variance: if a light curve has a total variance of :math:`S^2`,
+    and each point has an errorbar :math:`\sigma_{err}`, the *excess variance*
+    is defined as
+
+    .. math:: \sigma_{XS} = S^2 - \overline{\sigma_{err}}^2;
+
+    the *normalized excess variance* is the excess variance divided by the
+    square of the mean intensity:
+
+    .. math:: \sigma_{NXS} = \dfrac{\sigma_{XS}}{\overline{x}^2};
+
+    the *fractional mean square variability amplitude*, or
+    :math:`F_{var}`, is finally defined as
+
+    .. math:: F_{var} = \sqrt{\dfrac{\sigma_{XS}}{\overline{x}^2}}
 
     Parameters
     ----------
     lc : a :class:`Lightcurve` object
     normalization : str
-        if 'fvar', return normalized square-root excess variance. If 'none',
-        return the unnormalized variance
-
+        if 'fvar', return the fractional mean square variability :math:`F_{var}`.
+        If 'none', return the unnormalized excess variance variance
+        :math:`\sigma_{XS}`. If 'norm_xs', return the normalized excess variance
+        :math:`\sigma_{XS}`
     Returns
     -------
     var_xs : float
@@ -393,37 +415,42 @@ def excess_variance(lc, normalization='fvar'):
     var_xs = lc_actual_var - lc_mean_var
     mean_lc = np.mean(lc.counts)
     mean_ctvar = mean_lc ** 2
+    var_nxs = var_xs / mean_lc ** 2
 
     fvar = np.sqrt(var_xs / mean_ctvar)
 
     N = len(lc.counts)
-    var_xs_err_A = np.sqrt(2 / N) * lc_mean_var / mean_lc ** 2
-    var_xs_err_B = np.sqrt(mean_lc ** 2 / N) * 2 * fvar / mean_lc
-    var_xs_err = np.sqrt(var_xs_err_A ** 2 + var_xs_err_B ** 2)
+    var_nxs_err_A = np.sqrt(2 / N) * lc_mean_var / mean_lc ** 2
+    var_nxs_err_B = np.sqrt(mean_lc ** 2 / N) * 2 * fvar / mean_lc
+    var_nxs_err = np.sqrt(var_nxs_err_A ** 2 + var_nxs_err_B ** 2)
 
-    fvar_err = var_xs_err / (2 * fvar)
+    fvar_err = var_nxs_err / (2 * fvar)
 
     if normalization == 'fvar':
         return fvar, fvar_err
+    elif normalization == 'norm_xs':
+        return var_nxs, var_nxs_err
     elif normalization == 'none' or normalization is None:
-        return var_xs, var_xs_err
+        return var_xs, var_nxs_err * mean_lc **2
 
 
 def create_window(N, window_type='uniform'):
     """ A method to create window functions commonly used in signal processing.
-        Windows supported are:
-        Hamming, Hanning, uniform(rectangular window), triangular window, blackmann window among others.
-
-        Parameters
-        ----------
-        N : int
-            Total number of data points in window. If negative, abs is taken.
-        window_type : {'uniform', 'parzen', 'hamming', 'hanning', 'traingular', 'welch', 'blackmann', 'flat-top'}, optional, default 'uniform'
-            Type of window to create.
-         Returns
-        -------
-        window: numpy.ndarray
-            Window function of length N.
+
+    Windows supported are:
+    Hamming, Hanning, uniform(rectangular window), triangular window, blackmann window among others.
+
+    Parameters
+    ----------
+    N : int
+        Total number of data points in window. If negative, abs is taken.
+    window_type : {'uniform', 'parzen', 'hamming', 'hanning', 'traingular', 'welch', 'blackmann', 'flat-top'}, optional, default 'uniform'
+        Type of window to create.
+
+    Returns
+    -------
+    window: numpy.ndarray
+        Window function of length N.
     """
 
     if not isinstance(N, int):
@@ -481,7 +508,7 @@ def create_window(N, window_type='uniform'):
     if window_type == 'welch':
         N_minus_1_by_2 = N_minus_1 / 2
         window = 1 - np.square((n - N_minus_1_by_2) / N_minus_1_by_2)
-        
+
     if window_type == 'blackmann':
         a0 = 0.42659
         a1 = 0.49656


Test patch:
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index c493c433..15fb489f 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -18,6 +18,70 @@ try:
 except ImportError:
     _H5PY_INSTALLED = False
 
+def fvar_fun(lc):
+    from stingray.utils import excess_variance
+    return excess_variance(lc, normalization='fvar')
+
+def nvar_fun(lc):
+    from stingray.utils import excess_variance
+    return excess_variance(lc, normalization='norm_xs')
+
+def evar_fun(lc):
+    from stingray.utils import excess_variance
+    return excess_variance(lc, normalization='none')
+
+
+class TestChunks(object):
+    @classmethod
+    def setup_class(cls):
+        dt = 0.1
+        tstart = 0
+        tstop = 100
+        times = np.arange(tstart, tstop, dt)
+        cls.gti = np.array([[tstart - dt/2, tstop - dt/2]])
+        # Simulate something *clearly* non-constant
+        counts = np.random.poisson(
+            10000 + 2000 * np.sin(2 * np.pi * times))
+
+        cls.lc = Lightcurve(times, counts, gti=cls.gti)
+
+    def test_analyze_lc_chunks_fvar_fracstep(self):
+
+        start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun,
+                                                     fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        fvar, fvar_err = res
+
+        assert np.allclose(start[0], self.gti[0, 0])
+        assert np.all(fvar > 0)
+        # This must be a clear measurement of fvar
+        assert np.all(fvar > fvar_err)
+
+    def test_analyze_lc_chunks_nvar_fracstep(self):
+        start, stop, res = self.lc.analyze_lc_chunks(20, fvar_fun,
+                                                     fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        fvar, fvar_err = res
+        start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun,
+                                                     fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        nevar, nevar_err = res
+        assert np.allclose(nevar, fvar**2, rtol=0.01)
+
+    def test_analyze_lc_chunks_nvar_fracstep(self):
+        start, stop, mean = self.lc.analyze_lc_chunks(20, np.mean,
+                                                      fraction_step=0.5)
+        start, stop, res = self.lc.analyze_lc_chunks(20, evar_fun,
+                                                     fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        evar, evar_err = res
+        start, stop, res = self.lc.analyze_lc_chunks(20, nvar_fun,
+                                                     fraction_step=0.5)
+        # excess_variance returns fvar and fvar_err
+        nevar, nevar_err = res
+        assert np.allclose(nevar * mean ** 2, evar, rtol=0.01)
+        assert np.allclose(nevar_err * mean ** 2, evar_err, rtol=0.01)
+
 
 class TestLightcurve(object):
 
@@ -33,7 +97,7 @@ class TestLightcurve(object):
         Demonstrate that we can create a trivial Lightcurve object.
         """
         lc = Lightcurve(self.times, self.counts)
-    
+
     def test_irregular_time_warning(self):
         """
         Check if inputting an irregularly spaced time iterable throws out
@@ -98,31 +162,6 @@ class TestLightcurve(object):
         assert start[0] == 0.5
         assert np.all(start + lc.dt / 2 == res)
 
-    def test_analyze_lc_chunks_fvar_fracstep(self):
-        dt = 0.1
-        tstart = 0
-        tstop = 100
-        times = np.arange(tstart, tstop, dt)
-        gti = np.array([[tstart - dt/2, tstop - dt/2]])
-        # Simulate something *clearly* non-constant
-        counts = np.random.poisson(
-            10000 + 2000 * np.sin(2 * np.pi * times))
-
-        lc = Lightcurve(times, counts, gti=gti)
-
-        def excvar(lc):
-            from stingray.utils import excess_variance
-            return excess_variance(lc, normalization='fvar')
-
-        start, stop, res = lc.analyze_lc_chunks(20, excvar, fraction_step=0.5)
-        # excess_variance returns fvar and fvar_err
-        res, res_err = res
-
-        assert np.allclose(start[0], gti[0, 0])
-        assert np.all(res > 0)
-        # This must be a clear measurement of fvar
-        assert np.all(res > res_err)
-
     def test_bin_edges(self):
         bin_lo = [0.5,  1.5,  2.5,  3.5]
         bin_hi = [1.5,  2.5,  3.5,  4.5]
@@ -365,12 +404,16 @@ class TestLightcurve(object):
         assert lc[0] == lc[1] == lc[2] == lc[3] == 2
 
     def test_slicing(self):
-        lc = Lightcurve(self.times, self.counts)
+        lc = Lightcurve(self.times, self.counts, gti=self.gti)
 
         assert np.all(lc[1:3].counts == np.array([2, 2]))
         assert np.all(lc[:2].counts == np.array([2, 2]))
+        assert np.all(lc[:2].gti == [[0.5, 2.5]])
         assert np.all(lc[2:].counts == np.array([2, 2]))
+        assert np.all(lc[2:].gti == [[2.5, 4.5]])
         assert np.all(lc[:].counts == np.array([2, 2, 2, 2]))
+        assert np.all(lc[::2].gti == [[0.5, 1.5], [2.5, 3.5]])
+        assert np.all(lc[:].gti == lc.gti)
         assert lc[:].mjdref == lc.mjdref
 
 
@@ -442,7 +485,7 @@ class TestLightcurve(object):
 
         lc1 = Lightcurve(self.times, self.counts, err_dist = "poisson")
         lc2 = Lightcurve(_times, _counts, err_dist = "gauss")
-        
+
         with warnings.catch_warnings(record=True) as w:
             lc3 = lc1.join(lc2)
             assert "We are setting the errors to zero." in str(w[1].message)