Andreas Specovius / ctools

Commit 18cf2954cd9f03e873d1d6773b96f7f221f09b81

Authored by Andreas Specovius
1 parent 1c81f76b

First implementation of script.

Inline Side-by-side
Showing 1 changed file with 368 additions and 0 deletions
examples/show_significance_distribution.py 0 → 100755
  View file @18cf295
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# Shows the distribution of significances in a given significance map.
  4 +#
  5 +# Copyright (C) 2011-2017 Andreas Specovius
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import sys
  22 +import gammalib
  23 +import cscripts
  24 +try:
  25 + import matplotlib.pyplot as plt
  26 + import matplotlib.mlab as mlab
  27 + plt.figure()
  28 + plt.close()
  29 +except (ImportError, RuntimeError):
  30 + print('This script needs the "matplotlib" module')
  31 + sys.exit()
  32 +try:
  33 + import numpy as np
  34 +except ImportError:
  35 + print('This script needs the "numpy" module')
  36 + sys.exit()
  37 +
  38 +try:
  39 + import scipy.optimize as optimize
  40 +except ImportError:
  41 + print('This script needs the "scipy" module')
  42 + sys.exit()
  43 +
  44 +
  45 +#
  46 +#
  47 +#
  48 +# What do I need for plotting:
  49 +#
  50 +#- significance map itself [mappath]
  51 +# Provided as a file path guiding to a GSkyMap in FITS format
  52 +# No default
  53 +#
  54 +#- included regions (for shrinking data space i.e. for cutting certain edges or at a certain radius) [includedreg]
  55 +# Provided as a file path to a FITS WCS map or DS9 regions file
  56 +# Default should be None, which means all is allowed
  57 +#
  58 +#- excluded regions [excludedreg]
  59 +# Provided as a file path to a FITS WCS map or DS9 regions file
  60 +# Default should be None, which means there are no excluded regions. In this case no bkg will be plotted and fitted
  61 +#
  62 +#- binning parameters for the histogram: [nbins], [sign_min], [sign_max]
  63 +# Provided as floats
  64 +# Default should be -5 to 8 with 130 bins (131 edges)
  65 +#
  66 +#- title of the histogram [title]
  67 +# Provided as string
  68 +# Default should be an empty string
  69 +#
  70 +#- Path to plotted file [plotfile]
  71 +#
  72 +#
  73 +
  74 +
  75 +# ================ #
  76 +# Helper functions #
  77 +# ================ #
  78 +def gaussian(x, norm, xmean, sigma):
  79 + """
  80 + Compute function values for gaussian function.
  81 + """
  82 + return norm*np.exp(-(x-xmean)**2/float(2*sigma**2))
  83 +
  84 +
  85 +def skymap_to_numpy_ndarray(skymap):
  86 + """
  87 + Return numpy.ndarray of pixel values of a skymap.
  88 + """
  89 + # Get map size
  90 + npix = skymap.npix()
  91 +
  92 + # Initialise out array
  93 + array = np.empty(npix)
  94 +
  95 + # Loop over map pixels
  96 + for ipix in range(npix):
  97 + array[ipix] = skymap[ipix]
  98 +
  99 + # Return
  100 + return array
  101 +
  102 +
  103 +def read_regions(fpath, init_map):
  104 + """
  105 + Read a regions WCS FITS or ds9 file and return a mask WCS map.
  106 +
  107 + Parameters
  108 + ----------
  109 + fpath : str
  110 + Path to regions file
  111 + init_map : GSkyMap
  112 + Sky map representing fov of interest. Used for initialisation of
  113 + regions map.
  114 + """
  115 +
  116 + # Create starter map for user fov filled with zeros
  117 + regions_map = init_map.copy()
  118 + regions_map *= 0
  119 +
  120 + # Take care about ds9 region files
  121 + if fpath.lower().endswith('.reg'):
  122 +
  123 + # Read ds9 regions
  124 + ds9_regions = gammalib.GSkyRegions(fpath)
  125 +
  126 + # Loop over regions
  127 + for reg in ds9_regions:
  128 +
  129 + # Make map from region
  130 + ds9_map = gammalib.GSkyRegionMap(reg)
  131 +
  132 + # Add ds9 region map to global map
  133 + regions_map += ds9_map.map()
  134 +
  135 + # Take care about fits wcs files
  136 + elif fpath.lower().endswith('.fits'):
  137 +
  138 + # Read wcs regions map
  139 + wcs_regions = gammalib.GSkyMap(fpath)
  140 +
  141 + # Add wcs regions map to global map
  142 + regions_map += wcs_regions
  143 +
  144 + else:
  145 + raise RuntimeError("Invalid regions file detected. Please provide " +
  146 + "a valid ds9 or FITS WCS regions file.")
  147 +
  148 + # Return
  149 + return regions_map
  150 +
  151 +
  152 +# ============================== #
  153 +# Plot significance distribution #
  154 +# ============================== #
  155 +def plot_significance_distribution(mappath, nbins, sigma_min, sigma_max,
  156 + includedreg, excludedreg, title):
  157 + """
  158 + Plot the significance distribution and return instance of pyplot
  159 + figure and axis.
  160 +
  161 + Parameters
  162 + ----------
  163 + mappath : str
  164 + Path to significance map
  165 + nbins : int
  166 + Number of bins in histogram
  167 + sigma_min : float
  168 + Lower limit of the x axis in the plot
  169 + sigma_max : float
  170 + Upper limit of the x axis in the plot
  171 + includedreg : str
  172 + Path to global included region file
  173 + excludedreg : str
  174 + Path to global excluded region file
  175 + title : str
  176 + Title of the plot
  177 +
  178 + Figure and axis objects are returned.
  179 + """
  180 +
  181 + # Read significance map
  182 + # =====================
  183 +
  184 + # Read map
  185 + significance_map = gammalib.GSkyMap(mappath+'[SIGNIFICANCE]')
  186 +
  187 +
  188 + # Read included regions
  189 + # =====================
  190 +
  191 + # Set availability to false
  192 + inclusion_available = False
  193 +
  194 + # Read in case provided
  195 + if len(includedreg) > 0:
  196 + # Read regions map
  197 + regions_map = read_regions(includedreg, init_map=significance_map)
  198 +
  199 + # Make numpy boolean mask
  200 + mask_include = skymap_to_numpy_ndarray(regions_map).astype(bool)
  201 +
  202 + # Enable masked map
  203 + inclusion_available = True
  204 +
  205 +
  206 + # Read excluded regions
  207 + # =====================
  208 +
  209 + # Set availability to false
  210 + src_exclusion_available = False
  211 +
  212 + # Read in case provided
  213 + if len(excludedreg) > 0:
  214 + # Read regions file
  215 + regions_map = read_regions(excludedreg, init_map=significance_map)
  216 +
  217 + # Make numpy boolean mask
  218 + mask_exclude = skymap_to_numpy_ndarray(regions_map).astype(bool)
  219 +
  220 + # Enable masked map
  221 + src_exclusion_available = True
  222 +
  223 +
  224 + # Convet maps to flat arrays and apply spatial cuts
  225 + # =================================================
  226 +
  227 + # significance data to flat array
  228 + data = skymap_to_numpy_ndarray(significance_map)
  229 +
  230 + # Initialise dummy mask allowing all pixels
  231 + mask = np.ones(significance_map.npix()).astype(bool)
  232 +
  233 + # Apply inclusion region spatial cuts
  234 + if inclusion_available:
  235 + mask &= mask_include
  236 + data = data[mask]
  237 +
  238 + # Apply exclusion region spatial cuts
  239 + if src_exclusion_available:
  240 + data_without_sources = skymap_to_numpy_ndarray(significance_map)
  241 + mask &= ~mask_exclude
  242 + data_without_sources = data_without_sources[mask]
  243 +
  244 +
  245 + # Compute bin edges to user parameters
  246 + # ====================================
  247 +
  248 + # Compute bin edges, hence use nbins+1
  249 + bin_edges = np.linspace(sigma_min, sigma_max, nbins+1)
  250 +
  251 +
  252 + # Begin with plotting
  253 + # ===================
  254 +
  255 + # Create figure
  256 + fig = plt.figure()
  257 + ax = fig.gca()
  258 +
  259 + # Draw significance histogram for full fov
  260 + y, _, _ = ax.hist(data, bins=bin_edges, histtype='step', color='k', label="significance")
  261 +
  262 + # Draw the default normal distribution. Scale to fit maximum of histogram
  263 + ax.plot(bin_edges, gaussian(bin_edges, y.max(),0,1), 'b-')
  264 + msg = 'mean: $0$\nwidth: $1$'
  265 + ax.text(0.98, 0.80, msg, ha='right', va='top', bbox=dict(edgecolor='blue', facecolor='white'), transform=ax.transAxes)
  266 +
  267 + # If exclusion was provided also handle significanes with excluded src data
  268 + if src_exclusion_available:
  269 +
  270 + # Draw histogram of significances for data without excluded sources
  271 + y, _, _ = ax.hist(data_without_sources, bins=bin_edges, histtype='step',
  272 + color='k', alpha=0.5, label="significance without exclusions")
  273 +
  274 + # Fit a normal distribution to the binned masked significance data
  275 + pars = [1.0, 0.0, 1.0]
  276 + bin_centers = 0.5*(bin_edges[:-1] + bin_edges[1:])
  277 +
  278 + pars_opt, pars_cov = optimize.curve_fit(gaussian, xdata=bin_centers,
  279 + ydata=y, p0=pars)
  280 +
  281 + # Compute parameters fit errors
  282 + pars_err = np.sqrt(np.diag(pars_cov))
  283 +
  284 + # Log fit result
  285 + print('Fit of normal distribution to excluded source significance data ' +\
  286 + 'results in: norm=%f+=%f xmean=%f+-%f sigma=%f+-%f' %
  287 + (pars_opt[0], pars_err[0], pars_opt[1], pars_err[1], pars_opt[2], pars_err[2])
  288 + )
  289 +
  290 + # Draw text box
  291 + msg = 'mean: $%.3f\pm%.3f$\nwidth: $%.3f\pm%.3f$' % (pars_opt[1], pars_err[1], pars_opt[2], pars_err[2])
  292 + ax.text(0.98, 0.95, msg, ha='right', va='top', bbox=dict(edgecolor='red', facecolor='white'), transform=ax.transAxes)
  293 +
  294 + # Plot the normal (2)
  295 + ax.plot(bin_edges, gaussian(bin_edges, *(pars_opt)), 'r-')
  296 +
  297 + # Configure the plot
  298 + ax.set_ylim(0.5, ax.get_ylim()[1]*2.0)
  299 + ax.set_xlim(sigma_min, sigma_max)
  300 + ax.set_title(title)
  301 + ax.set_xlabel("Significance")
  302 + ax.set_ylabel("Entries")
  303 + ax.grid()
  304 + ax.set_yscale('log')
  305 +
  306 + # Return
  307 + return fig, ax
  308 +
  309 +
  310 +# ============================== #
  311 +# Show significance distribution #
  312 +# ============================== #
  313 +def show_significance_distribution():
  314 + """
  315 + Plot and show significance distribution
  316 + """
  317 + # Set usage string
  318 + usage = 'show_significance_distribution.py [--nbins bins] ' +\
  319 + '[--smin sigma_min] [--smax sigma_max] [--include regfile] ' +\
  320 + '[--exclude regfile] [--title title] [-p plotfile] mapfile'
  321 +
  322 + # Set default options
  323 + options = [{'option': '--nbins', 'value': '130'},
  324 + {'option': '--smin', 'value': '-5'},
  325 + {'option': '--smax', 'value': '8'},
  326 + {'option': '--include', 'value': ''},
  327 + {'option': '--exclude', 'value': ''},
  328 + {'option': '--title', 'value': ''},
  329 + {'option': '-p', 'value': ''}]
  330 +
  331 + # Get arguments and options from command line arguments
  332 + args, options = cscripts.ioutils.get_args_options(options, usage)
  333 +
  334 + # Print help message
  335 + if args[0]=='--help' or args[0]=='-h':
  336 + print(usage)
  337 + sys.exit()
  338 +
  339 + # Extract script parameters from options
  340 + nbins = int(options[0]['value'])
  341 + sigma_min = int(options[1]['value'])
  342 + sigma_max = int(options[2]['value'])
  343 + includedreg = options[3]['value']
  344 + excludedreg = options[4]['value']
  345 + title = options[5]['value']
  346 + plotfile = options[6]['value']
  347 +
  348 + # Plot significance distribution histogram
  349 + fig, ax = plot_significance_distribution(args[0], nbins, sigma_min, sigma_max,
  350 + includedreg, excludedreg, title)
  351 +
  352 + # Show or save the plot
  353 + if len(plotfile) > 0:
  354 + plt.savefig(plotfile)
  355 + else:
  356 + plt.show()
  357 +
  358 + # Return
  359 + return
  360 +
  361 +
  362 +# ======================== #
  363 +# Main routine entry point #
  364 +# ======================== #
  365 +if __name__ == '__main__':
  366 +
  367 + # Plot and show significance distribution histogram
  368 + show_significance_distribution()
... ...