Make gammapy.time better

docs/time/index.rst

@@ -11,38 +11,64 @@ Introduction
 
 `gammapy.time` contains methods for timing analysis.
 
-At the moment there's almost no functionality here.
-Any contributions welcome, e.g.
+At the moment there isn't a lot of functionality yet ... contributions welcome!
 
-* utility functions for time computations
-* lightcurve analysis functions
-* pulsar periodogram
-* event dead time calculations and checks for constant rate
 
-Although where possible we shouldn't duplicate timing analysis functionality
-that's already available elsewhere. In some cases it's enough to refer
-gamma-ray astronomers to other packages, sometimes a convenience wrapper for
-our data formats or units might make sense.
+Getting Started
+===============
 
-Some references (please add what you find useful
+.. _time-lc:
 
-* https://github.com/astroML/gatspy
-* https://github.com/matteobachetti/MaLTPyNT
-* https://github.com/nanograv/PINT
-* http://www.astroml.org/modules/classes.html#module-astroML.time_series
-* http://www.astroml.org/book_figures/chapter10/index.html
-* http://docs.astropy.org/en/latest/api/astropy.stats.bayesian_blocks.html
-* https://github.com/samconnolly/DELightcurveSimulation
-* https://github.com/cokelaer/spectrum
-* https://github.com/YSOVAR/YSOVAR
-* http://nbviewer.ipython.org/github/YSOVAR/Analysis/blob/master/TimeScalesinYSOVAR.ipynb
+Lightcurve
+----------
 
+The `~gammapy.time.LightCurve` class can be used to read a lightcurve,
+plot it and compute some summary statistics:
 
-Getting Started
-===============
+.. code-block:: python
+
+    >>> from gammapy.time import LightCurve
+    >>> lc = LightCurve.read('$GAMMAPY_EXTRA/todo/make/example-lightcurve.fits.gz')
+    >>> lc.plot()
+    >>> lc.info()
+
+
+.. _time-variability:
+
+Variability test
+----------------
+
+The `~gammapy.time.exptest` function can be used to compute the significance
+of variability (compared to the null hypothesis of constant rate)
+for a list of event time differences.
+
+Here's an example how to use the `~gammapy.time.make_random_times_poisson_process` helper
+function to simulate a `~astropy.time.TimeDelta` array for a given constant rate
+and use `~gammapy.time.exptest` to assess the level of variability (1.3 sigma in this case,
+not variable):
 
-TODO: document
+.. code-block:: python
+
+    >>> from gammapy.time import exptest
+    >>> from gammapy.time import make_random_times_poisson_process
+    >>> rate = Quantity(100, 's^-1')
+    >>> time = make_random_times_poisson_process(size=size, rate=rate, random_state=0)
+    >>>
+    >>> mr = exptest(time_delta)
+    >>> print(mr)
+    1.3790634240947202
+
+
+See ``examples/example_exptest.py`` for a longer example.
+
+TODO: apply this to the 2FHL events and check which sources are variable as a nice example.
+
+.. code-block:: python
 
+    from gammapy.data import EventList
+    from gammapy.time import exptest
+    events = EventList.read('$GAMMAPY_EXTRA/datasets/fermi_2fhl/2fhl_events.fits.gz ', hdu='EVENTS')
+    # TODO: cone select events for 2FHL catalog sources, compute mr for each and print 10 most variable sources
 
 .. _time_handling:
 
@@ -126,6 +152,28 @@ Time differences
 TODO: discuss when to use `~astropy.time.TimeDelta` or `~astropy.units.Quantity` or [MET]_ floats and
 where one needs to convert between those and what to watch out for.
 
+Other codes
+===========
+
+Where possible we shouldn't duplicate timing analysis functionality
+that's already available elsewhere. In some cases it's enough to refer
+gamma-ray astronomers to other packages, sometimes a convenience wrapper for
+our data formats or units might make sense.
+
+Some references (please add what you find useful
+
+* https://github.com/astroML/gatspy
+* https://github.com/matteobachetti/MaLTPyNT
+* https://github.com/nanograv/PINT
+* http://www.astroml.org/modules/classes.html#module-astroML.time_series
+* http://www.astroml.org/book_figures/chapter10/index.html
+* http://docs.astropy.org/en/latest/api/astropy.stats.bayesian_blocks.html
+* https://github.com/samconnolly/DELightcurveSimulation
+* https://github.com/cokelaer/spectrum
+* https://github.com/YSOVAR/YSOVAR
+* http://nbviewer.ipython.org/github/YSOVAR/Analysis/blob/master/TimeScalesinYSOVAR.ipynb
+
+
 Reference/API
 =============
 

examples/example_exptest.py

@@ -1,9 +1,14 @@
-#  An example how to compute exptest for multiple runs.
-from gammapy.time import exptest
+"""
+An example for the exptest variability test.
+
+- Simulate constant rate events for several observations.
+- Check that the ``mr`` distribution is a standard normal, as expected.
+"""
+import numpy as np
+from scipy.stats import norm
 import matplotlib.pyplot as plt
 from astropy.table import Table
-from scipy.stats import norm
-import numpy as np
+from gammapy.time import exptest
 
 
 def simulation(n_obs):
@@ -34,8 +39,7 @@ def simulation(n_obs):
 
 
 def exptest_multi(table_events):
-    """
-    Compute mr value for each run and whole dataset.
+    """Compute mr value for each run and whole dataset.
 
     Parameter
     ---------
@@ -60,19 +64,12 @@ def exptest_multi(table_events):
 
     for i in range(0, len(table_obs)):
         time_delta_each_run = []
+
         for j in range(0, len(table_events) - 1):
-            if table_obs['n_events'][i] > 20 and table_obs['obs_id'][i] == table_events[
-                'obs_id'][j] and table_events['obs_id'][j] == table_events['obs_id'][j + 1]:
-                time_delta_each_run.append(
-                    (table_events['mjd'][
-                         j + 1] - table_events['mjd'][j]) * 0.5 * (
-                        table_events['expCount'][
-                            j + 1] + table_events['expCount'][j]))
-                time_delta_all.append(
-                    (table_events['mjd'][
-                         j + 1] - table_events['mjd'][j]) * 0.5 * (
-                        table_events['expCount'][
-                            j + 1] + table_events['expCount'][j]))
+            if table_obs['n_events'][i] > 20 and table_obs['obs_id'][i] == table_events['obs_id'][j] and table_events['obs_id'][j] == table_events['obs_id'][j + 1]:
+                time_delta_each_run.append((table_events['mjd'][j + 1] - table_events['mjd'][j]) * 0.5 * (table_events['expCount'][j + 1] + table_events['expCount'][j]))
+                time_delta_all.append((table_events['mjd'][j + 1] - table_events['mjd'][j]) * 0.5 * (table_events['expCount'][j + 1] + table_events['expCount'][j]))
+
         if len(time_delta_each_run) == 0:
             continue
         mr = exptest(time_delta_each_run)
@@ -98,9 +95,8 @@ def plot(m_value):
     print("mu:{:10.3f}".format(mu), " sigma:{:10.4f}".format(sigma))
     plt.xlabel('Mr value')
     plt.ylabel('counts')
-    plt.title(
-        r'$\mathrm{Histogram\ of\ IQ:}\ \mu=%.3f,\ \sigma=%.3f$' %
-        (mu, sigma))
+    title = r'$\mathrm{Histogram\ of\ IQ:}\ \mu={:.3f},\ \sigma={:.3f}$'.format(mu, sigma)
+    plt.title(title)
     plt.grid(True)
     plt.show()
 

gammapy/time/exptest.py

@@ -6,24 +6,31 @@
     'exptest',
 ]
 
+
 def exptest(time_delta):
-    """Compute Mr value, the level of variability for a certain period of time.
+    """Compute the level of variability for a certain period of time.
 
-     A single Mr value can be calculated, which shows the level of variability
-     for the whole period, or the Mr value for each run can be shown.
+    The level of variability is quantified by ``mr``, as defined in Prahl (1999).
 
+    For constant rate events, it follows a standard normal distribution,
+    i.e. it can be used directly as the significance of variability.
 
-    Ref:Prah(1999),http://adsabs.harvard.edu/abs/1999astro.ph..9399P
+    For example usage, see :ref:`time-variability`.
 
     Parameters
     ----------
-    time_delta : list of float
-        the time difference between two adjacent events
+    time_delta : array-like
+        Time differences between consecutive events
 
     Returns
     -------
     mr : float
         Level of variability
+
+    References
+    ----------
+    .. [1] Prahl (1999), "A fast unbinned test on event clustering in Poisson processes",
+       `Link <http://adsabs.harvard.edu/abs/1999astro.ph..9399P>`_
     """
     mean_time = np.mean(time_delta)
     normalized_time_delta = time_delta / mean_time
@@ -32,5 +39,7 @@ def exptest(time_delta):
     sum_time_all = np.sum([sum_time])
     m_value = sum_time_all / len(time_delta)
     # the numbers are from Prahl(1999), derived from simulations
-    mr = (m_value - (1 / 2.71828 - 0.189 / len(time_delta))) / (0.2427 / np.sqrt(len(time_delta)))
+    term1 = m_value - (1 / 2.71828 - 0.189 / len(time_delta))
+    term2 = 0.2427 / np.sqrt(len(time_delta))
+    mr = term1 / term2
     return mr

gammapy/time/simulate.py

@@ -4,33 +4,42 @@
 from ..utils.random import get_random_state
 
 __all__ = [
+    'make_random_size_poisson_process',
     'make_random_times_poisson_process',
 ]
 
+def make_random_size_poisson_process(rate):
+    """
+
+    Parameters
+    ----------
+    rate : `~astropy.units.Quantity`
+        Event rate (dimension: 1 / TIME)
+    dead_time : `~astropy.units.Quantity` or `~astropy.time.TimeDelta`, optional
+        Dead time after event (dimension: TIME)
+    random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
+        Defines random number generator initialisation.
+        Passed to `~gammapy.utils.random.get_random_state`.
 
-def make_random_times_poisson_process(size, rate, dead_time=TimeDelta(0, format='sec'),
+    Returns
+    -------
+    size : int
+        Number of events
+    """
+    size = rate * (1obs_time
+
+
+def make_random_times_poisson_process(size, dead_time=TimeDelta(0, format='sec'),
                                       random_state='random-seed'):
     """Make random times assuming a Poisson process.
 
     This function can be used to test event time series,
     to have a comparison what completely random data looks like.
 
-    For the implementation see
-    `here <http://preshing.com/20111007/how-to-generate-random-timings-for-a-poisson-process/>`__ and
-    `here <http://stackoverflow.com/questions/1155539/how-do-i-generate-a-poisson-process>`__,
-    as well as `numpy.random.exponential`.
-
-    TODO: I think usually one has a given observation duration,
-    not a given number of events to generate.
-    Implementing this is more difficult because then the number
-    of samples to generate is variable.
-
     Parameters
     ----------
     size : int
         Number of samples
-    rate : `~astropy.units.Quantity`
-        Event rate (dimension: 1 / TIME)
     dead_time : `~astropy.units.Quantity` or `~astropy.time.TimeDelta`, optional
         Dead time after event (dimension: TIME)
     random_state : {int, 'random-seed', 'global-rng', `~numpy.random.RandomState`}
@@ -41,6 +50,9 @@ def make_random_times_poisson_process(size, rate, dead_time=TimeDelta(0, format=
     -------
     time : `~astropy.time.TimeDelta`
         Time differences (second) after time zero.
+
+    Examples
+    --------
     """
     # initialise random number generator
     random_state = get_random_state(random_state)