Repository: StingraySoftware/stingray
Pull request #386 (2019-02-04 23:27:24)

Add n_summed_spectra to Z2n stats
A missing bit of information in the Z2 statistics. We have the number of summed spectra in the PDS detection level, why not adding it to Z2 :)


Patch:
diff --git a/.travis.yml b/.travis.yml
index c39b6267..d21a558a 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -34,7 +34,7 @@ env:
 matrix:
     fast_finish: true
     include:
-        - env: PYTHON_VERSION=3.7 SETUP_CMD='egg_info' SETUP_XVFB=False
+        - env: PYTHON_VERSION=3.7 SETUP_CMD='egg_info'
 
         - env: PYTHON_VERSION=3.7
 
@@ -64,13 +64,18 @@ matrix:
         # Try Astropy development version
         - env: PYTHON_VERSION=3.7 ASTROPY_VERSION=development
 
+        - env: PYTHON_VERSION=2.7
+
     allow_failures:
+        - env: PYTHON_VERSION=3.7 SETUP_CMD='egg_info'
         # Try older scipy versions
         - env: PYTHON_VERSION=3.4 SCIPY_VERSION=0.16.0 NUMPY_VERSION=1.13 ASTROPY_VERSION=lts
 
         # Try Astropy development version
         - env: PYTHON_VERSION=3.7 ASTROPY_VERSION=development
 
+        - env: PYTHON_VERSION=2.7
+
 install:
 
     # We now use the ci-helpers package to set up our testing environment.
diff --git a/astropy_helpers b/astropy_helpers
index 231c409a..54d4ca98 160000
--- a/astropy_helpers
+++ b/astropy_helpers
@@ -1,1 +1,1 @@
-Subproject commit 231c409a632dcbf2beae1c2dea5b843d81ede511
+Subproject commit 54d4ca984ec81321001f72c888076376381ec0a6
diff --git a/stingray/_astropy_init.py b/stingray/_astropy_init.py
index 1e229e19..6d94a38b 100644
--- a/stingray/_astropy_init.py
+++ b/stingray/_astropy_init.py
@@ -1,6 +1,6 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
 
-__all__ = ['__version__', '__githash__']
+__all__ = ['__version__', '__githash__', 'test']
 
 # this indicates whether or not we are in the package's setup.py
 try:
@@ -22,20 +22,98 @@ try:
 except ImportError:
     __githash__ = ''
 
+# set up the test command
+def _get_test_runner():
+    import os
+    from astropy.tests.helper import TestRunner
+    return TestRunner(os.path.dirname(__file__))
+
+def test(package=None, test_path=None, args=None, plugins=None,
+         verbose=False, pastebin=None, remote_data=False, pep8=False,
+         pdb=False, coverage=False, open_files=False, **kwargs):
+    """
+    Run the tests using `py.test <http://pytest.org/latest>`__. A proper set
+    of arguments is constructed and passed to `pytest.main`_.
+
+    .. _py.test: http://pytest.org/latest/
+    .. _pytest.main: http://pytest.org/latest/builtin.html#pytest.main
+
+    Parameters
+    ----------
+    package : str, optional
+        The name of a specific package to test, e.g. 'io.fits' or 'utils'.
+        If nothing is specified all default tests are run.
+
+    test_path : str, optional
+        Specify location to test by path. May be a single file or
+        directory. Must be specified absolutely or relative to the
+        calling directory.
+
+    args : str, optional
+        Additional arguments to be passed to pytest.main_ in the ``args``
+        keyword argument.
+
+    plugins : list, optional
+        Plugins to be passed to pytest.main_ in the ``plugins`` keyword
+        argument.
+
+    verbose : bool, optional
+        Convenience option to turn on verbose output from py.test_. Passing
+        True is the same as specifying ``'-v'`` in ``args``.
+
+    pastebin : {'failed','all',None}, optional
+        Convenience option for turning on py.test_ pastebin output. Set to
+        ``'failed'`` to upload info for failed tests, or ``'all'`` to upload
+        info for all tests.
 
-if not _ASTROPY_SETUP_:  # noqa
+    remote_data : bool, optional
+        Controls whether to run tests marked with @remote_data. These
+        tests use online data and are not run by default. Set to True to
+        run these tests.
+
+    pep8 : bool, optional
+        Turn on PEP8 checking via the `pytest-pep8 plugin
+        <http://pypi.python.org/pypi/pytest-pep8>`_ and disable normal
+        tests. Same as specifying ``'--pep8 -k pep8'`` in ``args``.
+
+    pdb : bool, optional
+        Turn on PDB post-mortem analysis for failing tests. Same as
+        specifying ``'--pdb'`` in ``args``.
+
+    coverage : bool, optional
+        Generate a test coverage report.  The result will be placed in
+        the directory htmlcov.
+
+    open_files : bool, optional
+        Fail when any tests leave files open.  Off by default, because
+        this adds extra run time to the test suite.  Requires the
+        `psutil <https://pypi.python.org/pypi/psutil>`_ package.
+
+    parallel : int, optional
+        When provided, run the tests in parallel on the specified
+        number of CPUs.  If parallel is negative, it will use the all
+        the cores on the machine.  Requires the
+        `pytest-xdist <https://pypi.python.org/pypi/pytest-xdist>`_ plugin
+        installed. Only available when using Astropy 0.3 or later.
+
+    kwargs
+        Any additional keywords passed into this function will be passed
+        on to the astropy test runner.  This allows use of test-related
+        functionality implemented in later versions of astropy without
+        explicitly updating the package template.
+
+    """
+    test_runner = _get_test_runner()
+    return test_runner.run_tests(
+        package=package, test_path=test_path, args=args,
+        plugins=plugins, verbose=verbose, pastebin=pastebin,
+        remote_data=remote_data, pep8=pep8, pdb=pdb,
+        coverage=coverage, open_files=open_files, **kwargs)
+
+if not _ASTROPY_SETUP_:
     import os
     from warnings import warn
-    from astropy.config.configuration import (
-        update_default_config,
-        ConfigurationDefaultMissingError,
-        ConfigurationDefaultMissingWarning)
-
-    # Create the test function for self test
-    from astropy.tests.runner import TestRunner
-    test = TestRunner.make_test_runner_in(os.path.dirname(__file__))
-    test.__test__ = False
-    __all__ += ['test']
+    from astropy import config
 
     # add these here so we only need to cleanup the namespace at the end
     config_dir = None
@@ -45,16 +123,16 @@ if not _ASTROPY_SETUP_:  # noqa
         config_template = os.path.join(config_dir, __package__ + ".cfg")
         if os.path.isfile(config_template):
             try:
-                update_default_config(
+                config.configuration.update_default_config(
                     __package__, config_dir, version=__version__)
             except TypeError as orig_error:
                 try:
-                    update_default_config(__package__, config_dir)
-                except ConfigurationDefaultMissingError as e:
-                    wmsg = (e.args[0] +
-                            " Cannot install default profile. If you are "
+                    config.configuration.update_default_config(
+                        __package__, config_dir)
+                except config.configuration.ConfigurationDefaultMissingError as e:
+                    wmsg = (e.args[0] + " Cannot install default profile. If you are "
                             "importing from source, this is expected.")
-                    warn(ConfigurationDefaultMissingWarning(wmsg))
+                    warn(config.configuration.ConfigurationDefaultMissingWarning(wmsg))
                     del e
-                except Exception:
+                except:
                     raise orig_error
diff --git a/stingray/gti.py b/stingray/gti.py
index def76e49..95382c40 100644
--- a/stingray/gti.py
+++ b/stingray/gti.py
@@ -100,6 +100,9 @@ def check_gtis(gti):
         If GTIs have overlapping or displaced values
     """
     gti = np.asarray(gti)
+    if len(gti) < 1:
+        raise ValueError("Empty GTIs")
+
     if len(gti) != gti.shape[0] or len(gti.shape) != 2 or \
             len(gti) != gti.shape[0]:
         raise TypeError("Please check formatting of GTIs. They need to be"
@@ -158,7 +161,6 @@ def create_gti_mask_jit(time, gtis, mask, gti_mask, min_length=0):  # pragma: no
             if length < min_length:
                 next_gti = True
                 continue
-
             next_gti = False
             gti_mask[gti_el] = True
 
@@ -229,10 +231,15 @@ def create_gti_mask(time, gtis, safe_interval=0, min_length=0,
                                         dt=dt, epsilon=epsilon)
 
     gtis = np.array(gtis, dtype=np.longdouble)
+
     check_gtis(gtis)
 
     dt = assign_value_if_none(dt, np.median(np.diff(time)))
 
+    lengths = gtis[:, 1] - gtis[:, 0]
+    good = lengths > max(min_length, dt)
+    gtis = gtis[good]
+
     mask = np.zeros(len(time), dtype=bool)
 
     if not isinstance(safe_interval, collections.Iterable):
diff --git a/stingray/pulse/pulsar.py b/stingray/pulse/pulsar.py
index 24f0f268..f7c11127 100644
--- a/stingray/pulse/pulsar.py
+++ b/stingray/pulse/pulsar.py
@@ -7,6 +7,8 @@ import collections
 import warnings
 from ..utils import simon, jit, mad
 from scipy.optimize import minimize, basinhopping, curve_fit
+from scipy import stats
+
 try:
     import pint.toa as toa
     import pint
@@ -259,7 +261,6 @@ def fold_events(times, *frequency_derivatives, **opts):
     start_phase, stop_phase = pulse_phase(np.array([start_time, stop_time]),
                                           *frequency_derivatives,
                                           to_1=False)
-
     raw_profile, bins = np.histogram(phases,
                                      bins=np.linspace(0, 1, nbin + 1),
                                      weights=weights)
@@ -355,6 +356,7 @@ def fold_detection_level(nbin, epsilon=0.01, ntrial=1):
     if ntrial > 1:
         simon("fold: The treatment of ntrial is very rough. Use with caution")
     from scipy import stats
+
     return stats.chi2.isf(epsilon/ntrial, nbin - 1)
 
 
@@ -366,7 +368,7 @@ def z_n(phase, n=2, norm=1):
     phase : array of floats
         The phases of the events
     n : int, default 2
-        The ``n`` in $Z^2_n$.
+        Number of harmonics, including the fundamental
 
     Other Parameters
     ----------------
@@ -395,7 +397,7 @@ def z_n(phase, n=2, norm=1):
                 for k in range(1, n + 1)])
 
 
-def z2_n_detection_level(n=2, epsilon=0.01, ntrial=1):
+def z2_n_detection_level(n=2, epsilon=0.01, ntrial=1, n_summed_spectra=1):
     """Return the detection level for the Z^2_n statistics.
 
     See Buccheri et al. (1983), Bendat and Piersol (1971).
@@ -403,7 +405,7 @@ def z2_n_detection_level(n=2, epsilon=0.01, ntrial=1):
     Parameters
     ----------
     n : int, default 2
-        The ``n`` in $Z^2_n$
+        The ``n`` in $Z^2_n$ (number of harmonics, including the fundamental)
     epsilon : float, default 0.01
         The fractional probability that the signal has been produced by noise
 
@@ -411,6 +413,8 @@ def z2_n_detection_level(n=2, epsilon=0.01, ntrial=1):
     ----------------
     ntrial : int
         The number of trials executed to find this profile
+    n_summed_spectra : int
+        Number of Z_2^n periodograms that are being averaged
 
     Returns
     -------
@@ -418,13 +422,13 @@ def z2_n_detection_level(n=2, epsilon=0.01, ntrial=1):
         The epoch folding statistics corresponding to a probability
         epsilon * 100 % that the signal has been produced by noise
     """
-    if ntrial > 1:
-        simon("Z2_n: The treatment of ntrial is very rough. Use with caution")
-    from scipy import stats
-    return stats.chi2.isf(epsilon / ntrial, 2 * n)
+    retlev = stats.chi2.isf(epsilon / ntrial, 2 * n_summed_spectra * n) \
+        / (n_summed_spectra)
+
+    return retlev
 
 
-def z2_n_probability(z2, n=2, ntrial=1):
+def z2_n_probability(z2, n=2, ntrial=1, n_summed_spectra=1):
     """Calculate the probability of a certain folded profile, due to noise.
 
     Parameters
@@ -432,12 +436,14 @@ def z2_n_probability(z2, n=2, ntrial=1):
     z2 : float
         A Z^2_n statistics value
     n : int, default 2
-        The ``n`` in $Z^2_n$
+        The ``n`` in $Z^2_n$ (number of harmonics, including the fundamental)
 
     Other Parameters
     ----------------
     ntrial : int
         The number of trials executed to find this profile
+    n_summed_spectra : int
+        Number of Z_2^n periodograms that were averaged to obtain z2
 
     Returns
     -------
@@ -446,8 +452,10 @@ def z2_n_probability(z2, n=2, ntrial=1):
     """
     if ntrial > 1:
         simon("Z2_n: The treatment of ntrial is very rough. Use with caution")
-    from scipy import stats
-    return stats.chi2.sf(z2, 2 * n) * ntrial
+
+    epsilon = ntrial * stats.chi2.sf(z2 * n_summed_spectra,
+                                     2 * n * n_summed_spectra)
+    return epsilon
 
 
 def fftfit_fun(profile, template, amplitude, phase):
@@ -657,16 +665,19 @@ def _plot_TOA_fit(profile, template, toa, mod=None, toaerr=None,
     import matplotlib.pyplot as plt
     from scipy.interpolate import interp1d
     import time
-    phases = np.arange(0, 1, 1 / len(profile))
+    phases = np.arange(0, 2, 1 / len(profile))
+    profile = np.concatenate((profile, profile))
+    template = np.concatenate((template, template))
     if mod is None:
         mod = interp1d(phases, template, fill_value='extrapolate')
 
     fig = plt.figure()
-    plt.plot(phases - np.floor(phases), profile, drawstyle='steps-mid')
+    plt.plot(phases, profile, drawstyle='steps-mid')
     fine_phases = np.linspace(0, 1, 1000)
     fine_phases_shifted = fine_phases - toa / period + additional_phase
-    plt.plot(fine_phases,
-             mod(fine_phases_shifted - np.floor(fine_phases_shifted)))
+    model = mod(fine_phases_shifted - np.floor(fine_phases_shifted))
+    model = np.concatenate((model, model))
+    plt.plot(np.linspace(0, 2, 2000), model)
     if toaerr is not None:
         plt.axvline((toa - toaerr) / period)
         plt.axvline((toa + toaerr) / period)
diff --git a/stingray/utils.py b/stingray/utils.py
index 5f667f4a..0547e119 100644
--- a/stingray/utils.py
+++ b/stingray/utils.py
@@ -17,6 +17,7 @@ try:
     from numba import jit
 
     HAS_NUMBA = True
+    from numba import njit, prange
 except ImportError:
     warnings.warn("Numba not installed. Faking it")
 
@@ -30,6 +31,19 @@ except ImportError:
 
             return wrapped_f
 
+    class njit(object):
+        def __init__(self, *args, **kwargs):
+            pass
+
+        def __call__(self, func):
+            def wrapped_f(*args, **kwargs):
+                return func(*args, **kwargs)
+
+            return wrapped_f
+
+    def prange(x):
+        return range(x)
+
 try:
     from statsmodels.robust import mad as mad  # pylint: disable=unused-import
 except ImportError:


Test patch:
diff --git a/stingray/tests/test_gti.py b/stingray/tests/test_gti.py
index e7d2007d..e1f55380 100644
--- a/stingray/tests/test_gti.py
+++ b/stingray/tests/test_gti.py
@@ -237,3 +237,9 @@ class TestGTI(object):
 
         with pytest.raises(ValueError):
             check_gtis([[1, 0]])
+
+    def test_check_gti_fails_empty(self):
+        with pytest.raises(ValueError) as excinfo:
+            check_gtis([])
+        assert 'Empty' in str(excinfo)
+