Renamed cython functions, removed overlap.py

gammapy/detect/_test_statistics.pyx → gammapy/detect/_test_statistics_cython.pyx

@@ -12,9 +12,9 @@ FLUX_FACTOR = 1E-12
 
 @cython.cdivision(True)
 @cython.boundscheck(False)
-def _f_cash_root(np.float_t x, np.ndarray[np.float_t, ndim=2] counts,
-                 np.ndarray[np.float_t, ndim=2] background,
-                 np.ndarray[np.float_t, ndim=2] model):
+def _f_cash_root_cython(np.float_t x, np.ndarray[np.float_t, ndim=2] counts,
+                        np.ndarray[np.float_t, ndim=2] background,
+                        np.ndarray[np.float_t, ndim=2] model):
     cdef np.float_t sum
     cdef unsigned int i, j, ni, nj
     ni = counts.shape[1]
@@ -29,9 +29,9 @@ def _f_cash_root(np.float_t x, np.ndarray[np.float_t, ndim=2] counts,
 
 @cython.cdivision(True)
 @cython.boundscheck(False)
-def __amplitude_bounds(np.ndarray[np.float_t, ndim=2] counts,
-                       np.ndarray[np.float_t, ndim=2] background,
-                       np.ndarray[np.float_t, ndim=2] model):
+def _amplitude_bounds_cython(np.ndarray[np.float_t, ndim=2] counts,
+                             np.ndarray[np.float_t, ndim=2] background,
+                             np.ndarray[np.float_t, ndim=2] model):
     cdef float s_model = 0, s_counts = 0, sn, sn_min = 1E14, c_min
     cdef float b_min, b_max
     cdef unsigned int i, j, ni, nj
@@ -56,8 +56,8 @@ cdef extern from "math.h":
     float log(float x)
 
 
-def _cash(np.ndarray[np.float_t, ndim=2] counts,
-          np.ndarray[np.float_t, ndim=2] model):
+def _cash_cython(np.ndarray[np.float_t, ndim=2] counts,
+                 np.ndarray[np.float_t, ndim=2] model):
     cdef unsigned int i, j, ni, nj
     ni = counts.shape[1]
     nj = counts.shape[0]
@@ -70,8 +70,8 @@ def _cash(np.ndarray[np.float_t, ndim=2] counts,
     return cash
 
 
-def _cash_sum(np.ndarray[np.float_t, ndim=2] counts,
-              np.ndarray[np.float_t, ndim=2] model):
+def _cash_sum_cython(np.ndarray[np.float_t, ndim=2] counts,
+                     np.ndarray[np.float_t, ndim=2] model):
     cdef np.float_t sum = 0
     cdef unsigned int i, j, ni, nj
     ni = counts.shape[1]

gammapy/detect/test_statistics.py

@@ -16,7 +16,8 @@
 from astropy.nddata.utils import extract_array
 from astropy.io import fits
 
-from _test_statistics import _cash, __amplitude_bounds, _cash_sum, _f_cash_root
+from ._test_statistics_cython import (_cash_cython, _amplitude_bounds_cython,
+                                      _cash_sum_cython, _f_cash_root_cython)
 from ..irf import multi_gauss_psf_kernel
 from ..morphology import Shell2D
 from ..extern.zeros import newton
@@ -108,7 +109,7 @@ def f_cash_root(x, counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    return _f_cash_root(x, counts, background, model)
+    return _f_cash_root_cython(x, counts, background, model)
 
 
 def f_cash(x, counts, background, model):
@@ -126,7 +127,7 @@ def f_cash(x, counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    return _cash_sum(counts, background + x * FLUX_FACTOR * model)
+    return _cash_sum_cython(counts, background + x * FLUX_FACTOR * model)
 
 
 def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto',
@@ -152,7 +153,7 @@ def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto
         See `~gammapy.irf.multi_gauss_psf` for details.
     scales : list ([0])
         List of scales to use for TS map computation.
-    downsample : int (auto)
+    downsample : int ('auto')
         Down sampling factor. Can be set to 'auto' if the down sampling
         factor should be chosen automatically.
     residual : bool (False)
@@ -309,16 +310,28 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
     threshold : float (None)
         If the TS value corresponding to the initial flux estimate is not above
         this threshold, the optimizing step is omitted to save computing time.
-    sqrt : bool
-        Whether to return a sqrt(TS) map.
-    clean : bool
-        Whether to apply an edge cleaning algorithm to the TS map.
 
     Returns
     -------
     TS : `TSMapResult`
         `TSMapResult` object.
 
+
+    Notes
+    -----
+    Negative :math:`TS` values are defined as following:
+
+        .. math::
+
+            TS = \\left \\{
+                     \\begin{array}{ll}
+                       -TS & : \\textnormal{if} \\ F < 0 \\\\
+                       \\ \\ TS & : \\textnormal{else}
+                     \\end{array}
+                   \\right.
+
+    Where :math:`F` is the fitted flux amplitude.
+
     References
     ----------
     [Stewart2009]_
@@ -341,7 +354,7 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
         flux = convolve(flux, tophat.array) / CONTAINMENT
 
     # Compute null statistics for the whole map
-    C_0_map = _cash(counts.astype(float), background.astype(float))
+    C_0_map = _cash_cython(counts.astype(float), background.astype(float))
 
     x_min, x_max = kernel.shape[1] // 2, counts.shape[1] - kernel.shape[1] // 2
     y_min, y_max = kernel.shape[0] // 2, counts.shape[0] - kernel.shape[0] // 2
@@ -374,7 +387,7 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
     amplitudes[j, i] = [_[1] for _ in results]
     niter[j, i] = [_[2] for _ in results]
 
-    # Compute edge cleaned sqrt TS map
+    # Compute negative TS values
     sqrt_TS = np.where(TS > 0, np.sqrt(TS), -np.sqrt(-TS))
     return TSMapResult(ts=TS, sqrt_ts=sqrt_TS, amplitude=amplitudes,
                        niter=niter, runtime=np.round(time() - t_0, 2))
@@ -460,7 +473,7 @@ def _amplitude_bounds(counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    return __amplitude_bounds(counts, background, model)
+    return _amplitude_bounds_cython(counts, background, model)
 
 
 def _root_amplitude(counts, background, model, flux):

gammapy/detect/tests/test_test_statistics.py

@@ -31,6 +31,7 @@ def test_compute_ts_map():
     result = compute_ts_map(data['counts'], data['background'], data['exposure'],
                             kernel)
     for name, order in zip(['ts', 'amplitude', 'niter'], [2, 5, 0]):
+        result[name] = np.nan_to_num(result[name])
         result[name] = upsample_2N(result[name], 2, order=order)
 
     assert_allclose(1705.840212274973, result.ts[99, 99], rtol=1e-3)