Added Cython cash and model function

gammapy/detect/_test_statistics.pyx

@@ -0,0 +1,83 @@
+# 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
+
+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):
+    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(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
+    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(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]
+    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(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]
+    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

gammapy/detect/test_statistics.py

@@ -16,7 +16,7 @@
 from astropy.nddata.utils import extract_array
 from astropy.io import fits
 
-from ..stats import cash
+from _test_statistics import _cash, __amplitude_bounds, _cash_sum, _f_cash_root
 from ..irf import multi_gauss_psf_kernel
 from ..morphology import Shell2D
 from ..extern.zeros import newton
@@ -41,6 +41,8 @@ class TSMapResult(Bunch):
     ----------
     ts : `~numpy.ndarray`
         Estimated TS map
+    ts : `~numpy.ndarray`
+        Estimated sqrt(TS) map
     amplitude : `~numpy.ndarray`
         Estimated best fit flux amplitude map
     niter : `~numpy.ndarray`
@@ -60,13 +62,14 @@ def read(filename):
         """
         hdu_list = fits.open(filename)
         ts = hdu_list['ts'].data
+        sqrt_ts = hdu_list['sqrt_ts'].data
         amplitude = hdu_list['amplitude'].data
         niter = hdu_list['niter'].data
         scale = hdu_list[0].header['SCALE']
         if scale == 'max':
             scale = hdu_list['scale'].data
         morphology = hdu_list[0].header['MORPH']
-        return TSMapResult(ts=ts, amplitude=amplitude, niter=niter,
+        return TSMapResult(ts=ts, sqrt_ts=sqrt_ts, amplitude=amplitude, niter=niter,
                            scale=scale, morphology=morphology)
 
     def write(self, filename, header, overwrite=False):
@@ -83,7 +86,7 @@ def write(self, filename, header, overwrite=False):
             hdu_list.append(fits.ImageHDU(self.scale.astype('float32'), header))
         else:
             header['SCALE'] = self.scale,  'Source morphology scale parameter.'
-        for key in ['ts', 'amplitude', 'niter']:
+        for key in ['ts', 'sqrt_ts', 'amplitude', 'niter']:
             header['EXTNAME'] = key
             header['HDUNAME'] = key
             hdu_list.append(fits.ImageHDU(self[key].astype('float32'), header))
@@ -105,7 +108,7 @@ def f_cash_root(x, counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    return (model * (counts / (x * FLUX_FACTOR * model + background) - 1)).sum()
+    return _f_cash_root(x, counts, background, model)
 
 
 def f_cash(x, counts, background, model):
@@ -123,8 +126,7 @@ def f_cash(x, counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    with np.errstate(invalid='ignore', divide='ignore'):
-        return cash(counts, background + x * FLUX_FACTOR * model).sum()
+    return _cash_sum(counts, background + x * FLUX_FACTOR * model)
 
 
 def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto',
@@ -169,7 +171,7 @@ def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto
 
     for scale in scales:
         logging.info('Computing {0}TS map for scale {1:.3f} deg and {2}'
-                     ' morphology.'.format('residual' if residual else '',
+                     ' morphology.'.format('residual ' if residual else '',
                                            scale, morphology))
 
         # Sample down and require that scale parameters is at least 5 pix
@@ -219,16 +221,16 @@ def compute_ts_map_multiscale(maps, psf_parameters, scales=[0], downsample='auto
 
         # Compute TS map
         if residual:
-            background = (maps_['background'] + maps_['diffuse'] + maps_['OnModel'])
+            background = (maps_['background'] + maps_['diffuse'] + maps_['onmodel'])
         else:
-            background = maps_['background'] + maps_['diffuse']
+            background = maps_['background']  # + maps_['diffuse']
         ts_results = compute_ts_map(maps_['on'], background, maps_['expgammamap'],
                                     kernel, *args, **kwargs)
         logging.info('TS map computation took {0:.1f} s \n'.format(ts_results.runtime))
         ts_results['scale'] = scale
         ts_results['morphology'] = morphology
         if downsampled:
-            for name, order in zip(['ts', 'amplitude', 'niter'], [1, 1, 0]):
+            for name, order in zip(['ts', 'sqrt_ts', 'amplitude', 'niter'], [1, 1, 1, 0]):
                 ts_results[name] = upsample_2N(ts_results[name], factor,
                                                order=order, shape=shape)
         multiscale_result.append(ts_results)
@@ -307,6 +309,10 @@ 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
     -------
@@ -334,19 +340,20 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
             flux = (counts - background) / exposure / FLUX_FACTOR
         flux = convolve(flux, tophat.array) / CONTAINMENT
 
-    TS = np.zeros(counts.shape)
+    # Compute null statistics for the whole map
+    C_0_map = _cash(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
     positions = product(range(y_min, y_max), range(x_min, x_max))
 
     # Positions where exposure == 0 are not processed
     if mask is None:
-        mask = binary_erosion(exposure > 0, np.ones(np.array(kernel.shape) + 1))
+        mask = binary_erosion(exposure > 0, np.ones(np.array(kernel.shape) + 2))
     positions = [(j, i) for j, i in positions if mask[j][i]]
 
     wrap = partial(_ts_value, counts=counts, exposure=exposure,
-                   background=background, kernel=kernel, flux=flux,
+                   background=background, C_0_map=C_0_map, kernel=kernel, flux=flux,
                    method=method, optimizer=optimizer, threshold=threshold)
 
     if parallel:
@@ -360,16 +367,20 @@ def compute_ts_map(counts, background, exposure, kernel, mask=None, flux=None,
 
     # Set TS values at given positions
     j, i = zip(*positions)
+    TS = np.ones(counts.shape) * np.nan
     amplitudes = np.zeros(TS.shape)
     niter = np.zeros(TS.shape)
     TS[j, i] = [_[0] for _ in results]
     amplitudes[j, i] = [_[1] for _ in results]
     niter[j, i] = [_[2] for _ in results]
-    return TSMapResult(ts=TS, amplitude=amplitudes,
+
+    # Compute edge cleaned sqrt TS map
+    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))
 
 
-def _ts_value(position, counts, exposure, background, kernel, flux,
+def _ts_value(position, counts, exposure, background, C_0_map, kernel, flux,
               method, optimizer, threshold):
     """
     Compute TS value at a given pixel position i, j using the approach described
@@ -400,10 +411,10 @@ def _ts_value(position, counts, exposure, background, kernel, flux,
     counts_ = extract_array(counts, kernel.shape, position).astype(float)
     background_ = extract_array(background, kernel.shape, position).astype(float)
     exposure_ = extract_array(exposure, kernel.shape, position).astype(float)
+    C_0_ = extract_array(C_0_map, kernel.shape, position)
     model = (exposure_ * kernel._array).astype(float)
 
-    with np.errstate(invalid='ignore', divide='ignore'):
-        C_0 = cash(counts_, background_).sum()
+    C_0 = C_0_.sum()
 
     if threshold is not None:
         with np.errstate(invalid='ignore', divide='ignore'):
@@ -433,13 +444,13 @@ def _ts_value(position, counts, exposure, background, kernel, flux,
         C_1 = f_cash(amplitude, counts_, background_, model)
 
     # Compute and return TS value
-    return C_0 - C_1, amplitude * FLUX_FACTOR, niter
+    return (C_0 - C_1) * np.sign(amplitude), amplitude * FLUX_FACTOR, niter
 
 
 def _amplitude_bounds(counts, background, model):
     """
     Compute bounds for the root of `f_cash_root`.
- 
+
     Parameters
     ----------
     counts : `~numpy.ndarray`
@@ -449,14 +460,7 @@ def _amplitude_bounds(counts, background, model):
     model : `~numpy.ndarray`
         Source template (multiplied with exposure).
     """
-    # Check that counts slice contains at least one count
-    # if not, set flux estimate to -min(B/S).
-    with np.errstate(invalid='ignore', divide='ignore'):
-        sn = background / model
-    sn[counts == 0] = np.inf
-    index = np.unravel_index(np.nanargmin(sn), sn.shape)
-    bounds = np.array([counts[index], counts.sum()])
-    return (bounds / model.sum() - np.nanmin(sn)) / FLUX_FACTOR
+    return __amplitude_bounds(counts, background, model)
 
 
 def _root_amplitude(counts, background, model, flux):
@@ -620,4 +624,3 @@ def _flux_correlation_radius(kernel, containment=CONTAINMENT):
     # Containment radius of Tophat kernel is given by r_c_tophat = r_0 * sqrt(C)
     # by setting r_c = r_c_tophat we can estimate the equivalent containment radius r_0
     return r_c / np.sqrt(containment)
-

gammapy/detect/tests/test_test_statistics.py

@@ -30,10 +30,10 @@ 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]):
+    for name, order in zip(['ts', 'amplitude', 'niter'], [1, 1, 0]):
         result[name] = upsample_2N(result[name], 2, order=order)
 
-    assert_allclose([[99], [99]], np.where(result.ts == result.ts.max()))
-    assert_allclose(1705.840212274973, result.ts[99, 99], rtol=1e-3)
+    assert_allclose(1675.174869581638, result.ts[99, 99], rtol=1e-3)
+    assert_allclose([[99], [99]], np.where(result.ts == np.nanmax(result.ts)))
     assert_allclose(6, result.niter[99, 99])
-    assert_allclose(1.0227934338735763e-09, result.amplitude[99, 99], rtol=1e-3)
+    assert_allclose(1.0142865975138443e-09, result.amplitude[99, 99], rtol=1e-3)

gammapy/image/plotting.py

@@ -252,7 +252,7 @@ class GalacticPlaneSurveyPanelPlot(object):
     """
 
     def __init__(self, npanels=4, center=(0, 0), fov=(10, 1),
-                 xsize=10, xborder=0.5, yborder=0.5, yspacing=0.5):
+                 xsize=10, ysize=None, xborder=0.5, yborder=0.5, yspacing=0.5):
         """Compute panel parameters and make a matplotlib Figure.
         """
         import matplotlib.pyplot as plt
@@ -261,6 +261,7 @@ def __init__(self, npanels=4, center=(0, 0), fov=(10, 1),
                                                   center=center,
                                                   fov=fov,
                                                   xsize=xsize,
+                                                  ysize=ysize,
                                                   xborder=xborder,
                                                   yborder=yborder,
                                                   yspacing=yspacing)
@@ -529,7 +530,7 @@ def fitsfigure_add_colorbar_inset(ff, box, linewidth=1, color='w', normalize=Non
     ticks_pos = np.linspace(0, 1, n_ticks)
     if normalize is not None:
         ticks_pos = normalize.inverse(ticks_pos)
-    tick_labels = ['{0:' + ticklabel_format + '}'.format(_) for _ in ticks_pos]
+    tick_labels = [('{0:' + ticklabel_format + '}').format(_) for _ in ticks_pos]
     cbar.set_ticks(np.linspace(0, 1, n_ticks))
     cbar_axes.set_yticklabels(tick_labels, color=color)
     if label_position == 'bottom':

gammapy/scripts/ts_image.py

@@ -10,10 +10,10 @@
 
 def main(args=None):
     parser = get_parser(ts_image)
-    parser.add_argument('input', type=str,
+    parser.add_argument('input_file', type=str,
+                        help='Input data FITS file name')
+    parser.add_argument('output_file', type=str,
                         help='Input data FITS file name')
-    parser.add_argument('--folder', type=str, default='gaussian2d',
-                        help='Output folder name.')
     parser.add_argument('--psf', type=str, default='psf.json',
                         help='JSON file containing PSF information. ')
     parser.add_argument('--morphology', type=str, default='Gaussian2D',
@@ -44,8 +44,8 @@ def main(args=None):
     ts_image(**vars(args))
 
 
-def ts_image(input_file, psf, model, scales, downsample, residual, morphology,
-             width, overwrite, folder, threshold):
+def ts_image(input_file, output_file, psf, model, scales, downsample, residual,
+             morphology, width, overwrite, threshold):
     """
     Compute source model residual images.
 
@@ -77,14 +77,18 @@ def ts_image(input_file, psf, model, scales, downsample, residual, morphology,
         maps.append(fits.ImageHDU(data, header, 'OnModel'))
     results = compute_ts_map_multiscale(maps, psf_parameters, scales, downsample,
                                         residual, morphology, width)
-    folder = morphology.lower()
-    if not os.path.exists(folder):
+
+    folder, filename = os.path.split(output_file)
+    if not os.path.exists(folder) and folder != '':
         os.mkdir(folder)
 
     # Write results to file
     header = maps[0].header
-    for scale, result in zip(scales, results):
-        filename = os.path.join(folder, 'ts_{0:.3f}.fits'.format(scale))
-        logging.info('Writing {0}'.format(filename))
-        result.write(filename, header, overwrite=overwrite)
-
+    if len(results) > 1:
+        for scale, result in zip(scales, results):
+            filename_ = filename.replace('.fits', '_{0:.3f}.fits'.format(scale))
+            logging.info('Writing {0}'.format(os.path.join(folder, filename_)))
+            result.write(os.path.join(folder, filename_), header, overwrite=overwrite)
+    else:
+        logging.info('Writing {0}'.format(output_file))
+        results[0].write(output_file, header, overwrite=overwrite)