Merge pull request #252 from adonath/ts_map_performance_improvement

Implement TS map computation in Cython

docs/detect/index.rst

@@ -50,10 +50,9 @@ In the following the computation of a TS map for prepared Fermi survey data, whi
 	hdu_list = fits.open('all.fits.gz')
 	kernel = Gaussian2DKernel(5)
 	result = compute_ts_map(hdu_list['On'].data, hdu_list['Background'].data,
-							hdu_list['ExpGammaMap'].data, kernel, threshold=0)
+							hdu_list['ExpGammaMap'].data, kernel)
 
-The option ``threshold=0`` sets a minimal required TS value based on the initial flux estimate, that the pixel is
-processed at all. The function returns a `~gammapy.detect.TSMapResult` object, that bundles all relevant
+The function returns a `~gammapy.detect.TSMapResult` object, that bundles all relevant
 data. E.g. the time needed for the TS map computation can be checked by:
 
 .. code-block:: python
@@ -75,7 +74,7 @@ on the Fermi example dataset by:
 .. code-block:: bash
 
 	$ cd gammapy-extra/datasets/fermi_survey
-	$ gammapy-ts-image all.fits.gz --threshold 0 --scale 0 
+	$ gammapy-ts-image all.fits.gz ts_map_0.00.fits --scale 0 
 
 The command line tool additionally requires a psf json file, where the psf shape is defined by the parameters
 of a triple Gaussian model. See also `gammapy.irf.multi_gauss_psf_kernel`. By default the command line tool uses
@@ -88,7 +87,7 @@ Furthermore it is possible to compute residual TS maps. Using the following opti
 
 .. code-block:: bash
 
-	$ gammapy-ts-image all.fits.gz --threshold 0 --scale 0 --residual --model model.fits.gz 
+	$ gammapy-ts-image all.fits.gz residual_ts_map_0.00.fits --scale 0 --residual --model model.fits.gz 
 
 When ``--residual`` is set an excess model must be provided using the ``--model`` option.
 

gammapy/astro/source/pulsar.py

@@ -15,7 +15,7 @@
 class SimplePulsar(object):
     """Magnetic dipole spin-down model for a pulsar.
 
-    Reference: http://home.strw.leidenuniv.nl/~oberg/Pulsars/Verbunt_Heise_NeutronStars.pdf
+    Reference: http://www.cv.nrao.edu/course/astr534/Pulsars.html
 
     Parameters
     ----------
@@ -86,7 +86,7 @@ def magnetic_field(self):
 class Pulsar(SimplePulsar):
     """Magnetic dipole spin-down pulsar model.
 
-    Reference: http://home.strw.leidenuniv.nl/~oberg/Pulsars/Verbunt_Heise_NeutronStars.pdf
+    Reference: http://www.cv.nrao.edu/course/astr534/Pulsars.html
 
     Parameters
     ----------

gammapy/detect/_test_statistics_cython.pyx

@@ -0,0 +1,129 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+from __future__ import print_function, division
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+cdef np.float_t FLUX_FACTOR = 1E-12
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+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):
+    """
+    Function to find root of. Described in Appendix A, Stewart (2009).
+
+    Parameters
+    ----------
+    x : float
+        Model amplitude.
+    counts : `~numpy.ndarray`
+        Count map slice, where model is defined.
+    background : `~numpy.ndarray`
+        Background map slice, where model is defined.
+    model : `~numpy.ndarray`
+        Source template (multiplied with exposure).
+    """
+    cdef np.float_t sum
+    cdef unsigned int i, j, ni, nj
+    ni = counts.shape[1]
+    nj = counts.shape[0]
+    sum = 0
+    for j in range(nj):
+        for i in range(ni):
+            sum += (model[j, i] * (counts[j, i] / (x * FLUX_FACTOR * model[j, i]
+                                                   + background[j, i]) - 1))
+    return sum
+
+
+@cython.cdivision(True)
+@cython.boundscheck(False)
+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):
+    """
+    Compute bounds for the root of `_f_cash_root_cython`.
+
+    Parameters
+    ----------
+    counts : `~numpy.ndarray`
+        Count map.
+    background : `~numpy.ndarray`
+        Background map.
+    model : `~numpy.ndarray`
+        Source template (multiplied with exposure).
+    """
+
+    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
+    ni = counts.shape[1]
+    nj = counts.shape[0]
+    for j in range(ni):
+        for i in range(nj):
+            s_model += model[j, i]
+            if counts[j, i] > 0:
+                s_counts += counts[j, i]
+                if model[j, i] != 0:
+                    sn = background[j, i] / model[j, i]
+                    if sn < sn_min:
+                        sn_min = sn
+                        c_min = counts[j, i]
+    b_min = c_min / s_model - sn_min
+    b_max = s_counts / s_model - sn_min
+    return b_min / FLUX_FACTOR, b_max / FLUX_FACTOR
+
+
+cdef extern from "math.h":
+    float log(float x)
+
+
+def _cash_cython(np.ndarray[np.float_t, ndim=2] counts,
+                 np.ndarray[np.float_t, ndim=2] model):
+    """
+    Cash fit statistics.
+
+    Parameters
+    ----------
+    counts : `~numpy.ndarray`
+        Count map slice, where model is defined.
+    model : `~numpy.ndarray`
+        Source template (multiplied with exposure).
+    """
+
+    cdef unsigned int i, j, ni, nj
+    ni = counts.shape[1]
+    nj = counts.shape[0]
+    cdef np.ndarray[np.float_t, ndim=2] cash = np.zeros([nj, ni], dtype=float)
+
+    for j in range(nj):
+        for i in range(ni):
+            if model[j, i] > 0:
+                cash[j, i] = 2 * (model[j, i] - counts[j, i] * log(model[j, i]))
+    return cash
+
+
+def _cash_sum_cython(np.ndarray[np.float_t, ndim=2] counts,
+                     np.ndarray[np.float_t, ndim=2] model):
+    """
+    Summed cash fit statistics.
+
+    Parameters
+    ----------
+    counts : `~numpy.ndarray`
+        Count map slice, where model is defined.
+    model : `~numpy.ndarray`
+        Source template (multiplied with exposure).
+    """
+    cdef np.float_t sum = 0
+    cdef unsigned int i, j, ni, nj
+    ni = counts.shape[1]
+    nj = counts.shape[0]
+    for j in range(nj):
+        for i in range(ni):
+            if model[j, i] > 0:
+                sum += model[j, i] - counts[j, i] * log(model[j, i])
+    return 2 * sum