diff --git a/ChangeLog b/ChangeLog index 336dc42..d5e91a9 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2021-05-11 +2021-05-19 * Version 2.0.0 released ======================== diff --git a/NEWS b/NEWS index a132cf2..a92a8c2 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in ctools 2.0.0 -11 May 2021 +19 May 2021 Introduction @@ -255,7 +255,7 @@ the width, height and position angle of the rectangular sky regions (#3248). cspull - Generate pull distributions ------------------------------------ -None +The ASCII output file was replaced by a FITS file (#1707). csresmap - Generate residual map @@ -361,8 +361,9 @@ regions (#3248). Examples -------- -The show_butterfly and show_sensitivity.py scripts were adapted to deal also -with FITS output files (#1707). +The show_butterfly, show_pull_evolution.py, show_pull_histogram.py and +show_sensitivity.py scripts were adapted to deal also with FITS output files +(#1707). caldb diff --git a/README.md b/README.md index 63f1a13..8e9f654 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ctools information ================== -* Version: 2.0.0.dev (11 May 2021) +* Version: 2.0.0.dev (19 May 2021) * GammaLib dependency: 2.0.0.dev [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/) diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am index f731dd5..7cf54b0 100755 --- a/cscripts/Makefile.am +++ b/cscripts/Makefile.am @@ -1,7 +1,7 @@ ############################################################################# # Makefile for cscripts scripts # # ------------------------------------------------------------------------- # -# Copyright (C) 2016-2020 Juergen Knoedlseder # +# Copyright (C) 2016-2021 Juergen Knoedlseder # # ------------------------------------------------------------------------- # # # # This program is free software: you can redistribute it and/or modify # @@ -203,6 +203,7 @@ testdata = $(top_srcdir)/test/data/crab.xml \ $(top_srcdir)/test/data/pha.fits \ $(top_srcdir)/test/data/pha_channels.fits \ $(top_srcdir)/test/data/pull.dat \ + $(top_srcdir)/test/data/pull.fits \ $(top_srcdir)/test/data/sensitivity.dat \ $(top_srcdir)/test/data/sensitivity_integral.dat \ $(top_srcdir)/test/data/sensitivity.fits \ diff --git a/cscripts/cspull.par b/cscripts/cspull.par index 54ab578..0c18f70 100755 --- a/cscripts/cspull.par +++ b/cscripts/cspull.par @@ -24,7 +24,7 @@ caldb, s, a, prod2,,, "Calibration database" irf, s, a, South_0.5h,,, "Instrument response function" edisp, b, h, no,,, "Apply energy dispersion?" deadc, r, h, 0.98,0,1, "Deadtime correction factor" -outfile, f, a, pull.dat,,, "Output pull distribution file" +outfile, f, a, pull.fits,,, "Output pull distribution FITS file" # # Script parameters diff --git a/cscripts/cspull.py b/cscripts/cspull.py index b7e92c6..cc86e97 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -2,7 +2,7 @@ # ========================================================================== # This script generates the pull distribution for all model parameters. # -# Copyright (C) 2011-2018 Juergen Knoedlseder +# Copyright (C) 2011-2021 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -41,7 +41,8 @@ class cspull(ctools.csobservation): # Initialise application by calling the appropriate class constructor self._init_csobservation(self.__class__.__name__, ctools.__version__, argv) - # Initialise some members + # Initialise class members + self._fits = None self._nthreads = 0 # Return @@ -59,6 +60,7 @@ class cspull(ctools.csobservation): """ # Set pickled dictionary state = {'base' : ctools.csobservation.__getstate__(self), + 'fits' : self._fits, 'nthreads' : self._nthreads} # Return pickled dictionary @@ -74,11 +76,13 @@ class cspull(ctools.csobservation): Pickled instance """ ctools.csobservation.__setstate__(self, state['base']) + self._fits = state['fits'] self._nthreads = state['nthreads'] # Return return + # Private methods def _get_parameters(self): """ @@ -121,9 +125,12 @@ class cspull(ctools.csobservation): self['chatter'].integer() # Query some parameters - self['outfile'].filename() self['profile'].boolean() + # Read ahead output parameters + if self._read_ahead(): + self['outfile'].filename() + # Write input parameters into logger self._log_parameters(gammalib.TERSE) @@ -297,13 +304,18 @@ class cspull(ctools.csobservation): # Set name as a combination of model name and parameter # name separated by an underscore. In that way each # parameter has a unique name. - name = model.name()+'_'+par.name() + name = model.name() + '_' + par.name() + + # Set "Parameter", "Parameter_error" and "Parameter_Pull" + # column names + col_par = name + col_par_error = name+'_error' + col_par_pull = name+'_pull' - # Append parameter, Pull_parameter and e_parameter column - # names - colnames.append(name) - colnames.append('Pull_'+name) - colnames.append('e_'+name) + # Append column names + colnames.append(col_par) + colnames.append(col_par_error) + colnames.append(col_par_pull) # Compute pull for this parameter as the difference # (fitted - true) / error @@ -317,9 +329,9 @@ class cspull(ctools.csobservation): pull = 99.0 # Store results in dictionary - values[name] = fitted_value - values['Pull_'+name] = pull - values['e_'+name] = error + values[col_par] = fitted_value + values[col_par_error] = error + values[col_par_pull] = pull # Write results into logger value = '%.4f (%e +/- %e)' % (pull, fitted_value, error) @@ -331,6 +343,100 @@ class cspull(ctools.csobservation): # Return return result + def _create_fits(self, results): + """ + Create FITS file from results + + Parameters + ---------- + results : list of dict + List of result dictionaries + """ + # Gather headers for parameter columns + headers = [] + for colname in results[0]['colnames']: + if colname != 'LogL' and colname != 'Sim_Events' and \ + colname != 'Npred_Events': + headers.append(colname) + + # Create FITS table columns + nrows = len(results) + logl = gammalib.GFitsTableDoubleCol('LOGL', nrows) + nevents = gammalib.GFitsTableDoubleCol('NEVENTS', nrows) + npred = gammalib.GFitsTableDoubleCol('NPRED', nrows) + logl.unit('') + nevents.unit('counts') + npred.unit('counts') + columns = [] + for header in headers: + name = gammalib.toupper(header) + column = gammalib.GFitsTableDoubleCol(name, nrows) + column.unit('') + columns.append(column) + + # Fill FITS table columns + for i, result in enumerate(results): + logl[i] = result['values']['LogL'] + nevents[i] = result['values']['Sim_Events'] + npred[i] = result['values']['Npred_Events'] + for k, column in enumerate(columns): + column[i] = result['values'][headers[k]] + + # Initialise FITS Table with extension "PULL_DISTRIBUTION" + table = gammalib.GFitsBinTable(nrows) + table.extname('PULL_DISTRIBUTION') + + # Add keywors for compatibility with gammalib.GMWLSpectrum + table.card('INSTRUME', 'CTA', 'Name of Instrument') + table.card('TELESCOP', 'CTA', 'Name of Telescope') + + # Stamp header + self._stamp(table) + + # Set optional keyword values + if gammalib.toupper(self['onsrc'].string()) != 'NONE': + onrad = self['onrad'].real() + else: + onrad = 0.0 + if self['enumbins'].integer() > 0: + npix = self['npix'].integer() + binsz = self['binsz'].real() + coordsys = self['coordsys'].string() + proj = self['proj'].string() + else: + npix = 0 + binsz = 0.0 + coordsys = '' + proj = '' + + # Add script keywords + table.card('NPULLS', self['ntrials'].integer(), 'Number of pulls') + table.card('SEED', self['seed'].integer(), 'Seed value for pulls') + table.card('ONSRC', self['onsrc'].string(), 'Name of On surce for On/Off analysis') + table.card('ONRAD', onrad, '[deg] Radius of On region') + table.card('ENUMBINS', self['enumbins'].integer(), 'Number of energy bins') + table.card('NPIX', npix, 'Number of pixels for binned analysis') + table.card('BINSZ', binsz, 'Pixel size for binned analysis') + table.card('COORDSYS', coordsys, 'Coordinate system for binned analysis') + table.card('PROJ', proj, 'Projection for binned analysis') + table.card('STAT', self['statistic'].string(), 'Optimization statistic') + table.card('EDISP', self['edisp'].boolean(), 'Use energy dispersion?') + table.card('PROFILE', self['profile'].boolean(), 'Use likelihood profile method for errors?') + + # Append filled columns to fits table + table.append(logl) + table.append(nevents) + table.append(npred) + for column in columns: + table.append(column) + + # Create the FITS file now + self._fits = gammalib.GFits() + self._fits.append(table) + + # Return + return + # Public methods def run(self): @@ -356,6 +462,9 @@ class cspull(ctools.csobservation): # Get seed value seed = self['seed'].integer() + # Initialise results + results = [] + # If more than a single thread is requested then use multiprocessing if self._nthreads > 1: args = [(self, '_trial', i + seed) for i in range(ntrials)] @@ -367,16 +476,21 @@ class cspull(ctools.csobservation): # If multiprocessing was used then recover results and put them # into the log file if self._nthreads > 1: - result = poolresults[i][0] + #result = poolresults[i][0] + results.append(poolresults[i][0]) self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False) # ... otherwise make a trial now else: + + # Run trial result = self._trial(i + seed) - # Write out trial result - ioutils.write_csv_row(self['outfile'].filename().url(), i, - result['colnames'], result['values']) + # Append results + results.append(result) + + # Create FITS file + self._create_fits(results) # Return return @@ -396,6 +510,38 @@ class cspull(ctools.csobservation): # Return return + def save(self): + """ + Save pull FITS file + """ + # Write header + self._log_header1(gammalib.TERSE, 'Save pull distribution') + + # Continue only if FITS file is valid + if self._fits != None: + + # Get outmap parameter + outfile = self['outfile'].filename() + + # Log file name + self._log_value(gammalib.NORMAL, 'Pull distribution file', outfile.url()) + + # Save pull distribution + self._fits.saveto(outfile, self['clobber'].boolean()) + + # Return + return + + def pull_distribution(self): + """ + Return pull distribution FITS file + + Returns: + FITS file containing pull distribution + """ + # Return + return self._fits + # ======================== # # Main routine entry point # diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 485550d..c85d03b 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,6 +1,6 @@ .. _2.0: -ctools 2.0 release (11 May 2021) +ctools 2.0 release (19 May 2021) ================================ Overview diff --git a/doc/source/users/reference_manual/cspull.rst b/doc/source/users/reference_manual/cspull.rst index d48b9cf..7ee15d1 100755 --- a/doc/source/users/reference_manual/cspull.rst +++ b/doc/source/users/reference_manual/cspull.rst @@ -34,12 +34,10 @@ be performed in parallel over as many processes as the number of CPUs available machine. The maximum number of parallel processes can be set by the user through the ``nthreads`` hidden parameter. -:ref:`cspull` will generate an ASCII file in comma-separated value (CSV) format, -containing one row per pull. The first row is a header row providing the -column names. The following rows give the pull results, one row per pull. -This includes for each parameter the parameter value, error and pull. The -maximum likelihood value and the observed and the estimated number of counts -are also given. +:ref:`cspull` will generate a FITS file containing a binary table with one row per +pull. FITS table columns include the maximum log-likelihood value of the pull, +the number of simulated and fitted events, and for each parameter the parameter +value, error and pull. From the output file, pull distribution plots can be generated using for example the ``show_pull_histogram.py`` script in the examples folder. The diff --git a/examples/show_pull_evolution.py b/examples/show_pull_evolution.py index ba9cf40..d012743 100755 --- a/examples/show_pull_evolution.py +++ b/examples/show_pull_evolution.py @@ -2,7 +2,7 @@ # ========================================================================== # Shows the evolution of the mean and rms pull # -# Copyright (C) 2015-2020 Juergen Knoedlseder +# Copyright (C) 2015-2021 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -31,9 +31,80 @@ try: except ImportError: print('This script needs the "numpy" module') sys.exit() +import gammalib import cscripts +# ====================== # +# Read pull distribution # +# ====================== # +def read_pull_distribution(filename, parname): + """ + Read pull distribution + + Parameters + ---------- + filename : str + Name of FITS or CSV file + parname : str + Parameter name + + Returns + ------- + pulls : dict + Dictionary with pull distribution + """ + # Read pull distribution dependent on file type + fname = gammalib.GFilename(filename) + if fname.is_fits(): + pulls = read_pull_distribution_fits(filename, parname) + else: + pulls = cscripts.ioutils.read_pull_values(filename, parname) + + # Return pulls + return pulls + + +# ===================================== # +# Read pull distribution from FITS file # +# ===================================== # +def read_pull_distribution_fits(filename, parname): + """ + Read pull distribution from FITS file + + Parameters + ---------- + filename : str + Name of FITS file + parname : str + Parameter name + + Returns + ------- + pulls : dict + Dictionary with pull distribution + """ + # Open FITS file + fits = gammalib.GFits(filename) + + # Get pull distribution table + table = fits.table('PULL_DISTRIBUTION') + + # Get relevant column + column = table[parname] + + # Initialise vector + pulls = [] + + # Fill vectors + nrows = table.nrows() + for row in range(nrows): + pulls.append(column[row]) + + # Return + return pulls + + # =================== # # Plot pull evolution # # =================== # @@ -50,8 +121,11 @@ def plot_pull_evolution(filename, parname, plotfile): plotfile : str Plot filename """ - # Read values from CSV file - values = np.array(cscripts.ioutils.read_pull_values(filename, parname)) + # Read pull distribution from file + pulls = read_pull_distribution(filename, parname) + + # Create Numpy array + values = np.array(pulls) # Compute mean and rms nsamples = len(values) @@ -97,7 +171,7 @@ def show_pull_evolution(): usage = 'show_pull_evolution.py [-p plotfile] file parameter' # Set default options - options = [{'option': '-p', 'value': ''}] + options = [{'option': '-p', 'value': ''}] # Get arguments and options from command line arguments args, options = cscripts.ioutils.get_args_options(options, usage) diff --git a/examples/show_pull_histogram.py b/examples/show_pull_histogram.py index e90fb82..2328672 100755 --- a/examples/show_pull_histogram.py +++ b/examples/show_pull_histogram.py @@ -2,7 +2,7 @@ # ========================================================================== # Shows the pull histogram # -# Copyright (C) 2011-2020 Juergen Knoedlseder +# Copyright (C) 2011-2021 Juergen Knoedlseder # # This program is free software: you can redistribute it and/or modify # it under the terms of the GNU General Public License as published by @@ -31,9 +31,80 @@ try: except ImportError: print('This script needs the "numpy" module') sys.exit() +import gammalib import cscripts +# ====================== # +# Read pull distribution # +# ====================== # +def read_pull_distribution(filename, parname): + """ + Read pull distribution + + Parameters + ---------- + filename : str + Name of FITS or CSV file + parname : str + Parameter name + + Returns + ------- + pulls : dict + Dictionary with pull distribution + """ + # Read pull distribution dependent on file type + fname = gammalib.GFilename(filename) + if fname.is_fits(): + pulls = read_pull_distribution_fits(filename, parname) + else: + pulls = cscripts.ioutils.read_pull_values(filename, parname) + + # Return pulls + return pulls + + +# ===================================== # +# Read pull distribution from FITS file # +# ===================================== # +def read_pull_distribution_fits(filename, parname): + """ + Read pull distribution from FITS file + + Parameters + ---------- + filename : str + Name of FITS file + parname : str + Parameter name + + Returns + ------- + pulls : dict + Dictionary with pull distribution + """ + # Open FITS file + fits = gammalib.GFits(filename) + + # Get pull distribution table + table = fits.table('PULL_DISTRIBUTION') + + # Get relevant column + column = table[parname] + + # Initialise vector + pulls = [] + + # Fill vectors + nrows = table.nrows() + for row in range(nrows): + pulls.append(column[row]) + + # Return + return pulls + + # =================== # # Plot pull histogram # # =================== # @@ -52,8 +123,11 @@ def plot_pull_histogram(filename, parname, nbins, plotfile): plotfile : str Plot filename """ - # Read values from CSV file - values = np.array(cscripts.ioutils.read_pull_values(filename, parname)) + # Read pull distribution from file + pulls = read_pull_distribution(filename, parname) + + # Create Numpy array + values = np.array(pulls) # Create histogram _, bins, _ = plt.hist(values, nbins, range=[-4.0,4.0], diff --git a/test/data/pull.fits b/test/data/pull.fits new file mode 100644 index 0000000..6299c5f --- /dev/null +++ b/test/data/pull.fits @@ -0,0 +1,2 @@ +SIMPLE = T / file does conform to FITS standard BITPIX = 8 / number of bits per data pixel NAXIS = 0 / number of data axes EXTEND = T / FITS dataset may contain extensions COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H COMMENT FITS (Flexible Image Transport System) format is defined in 'AstronomyCOMMENT and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / 8-bit bytes NAXIS = 2 / 2-dimensional binary table NAXIS1 = 120 / width of table in bytes NAXIS2 = 10 / number of rows in table PCOUNT = 0 / size of special data area GCOUNT = 1 / one data group (required keyword) TFIELDS = 15 / number of table fields TTYPE1 = 'LOGL ' / label for field 1 TFORM1 = '1D ' / data format of field 1 TTYPE2 = 'NEVENTS ' / label for field 2 TFORM2 = '1D ' / data format of field 2 TUNIT2 = 'counts ' / physical unit of field2 TTYPE3 = 'NPRED ' / label for field 3 TFORM3 = '1D ' / data format of field 3 TUNIT3 = 'counts ' / physical unit of field3 TTYPE4 = 'CRAB_PREFACTOR' / label for field 4 TFORM4 = '1D ' / data format of field 4 TTYPE5 = 'CRAB_PREFACTOR_ERROR' / label for field 5 TFORM5 = '1D ' / data format of field 5 TTYPE6 = 'CRAB_PREFACTOR_PULL' / label for field 6 TFORM6 = '1D ' / data format of field 6 TTYPE7 = 'CRAB_INDEX' / label for field 7 TFORM7 = '1D ' / data format of field 7 TTYPE8 = 'CRAB_INDEX_ERROR' / label for field 8 TFORM8 = '1D ' / data format of field 8 TTYPE9 = 'CRAB_INDEX_PULL' / label for field 9 TFORM9 = '1D ' / data format of field 9 TTYPE10 = 'CTABACKGROUNDMODEL_PREFACTOR' / label for field 10 TFORM10 = '1D ' / data format of field 10 TTYPE11 = 'CTABACKGROUNDMODEL_PREFACTOR_ERROR' / label for field 11 TFORM11 = '1D ' / data format of field 11 TTYPE12 = 'CTABACKGROUNDMODEL_PREFACTOR_PULL' / label for field 12 TFORM12 = '1D ' / data format of field 12 TTYPE13 = 'CTABACKGROUNDMODEL_INDEX' / label for field 13 TFORM13 = '1D ' / data format of field 13 TTYPE14 = 'CTABACKGROUNDMODEL_INDEX_ERROR' / label for field 14 TFORM14 = '1D ' / data format of field 14 TTYPE15 = 'CTABACKGROUNDMODEL_INDEX_PULL' / label for field 15 TFORM15 = '1D ' / data format of field 15 EXTNAME = 'PULL_DISTRIBUTION' / name of this extension INSTRUME= 'CTA ' / Name of Instrument TELESCOP= 'CTA ' / Name of Telescope CREATOR = 'cspull v2.0.0.dev' / Program which created the file NPULLS = 10 / Number of pulls SEED = 1 / Seed value for pulls ONSRC = 'NONE ' / Name of On surce for On/Off analysis ONRAD = 0.0000000000E+00 / [deg] Radius of On region ENUMBINS= 0 / Number of energy bins NPIX = 0 / Number of pixels for binned analysis BINSZ = 0.0000000000E+00 / Pixel size for binned analysis COORDSYS= ' ' / Coordinate system for binned analysis PROJ = ' ' / Projection for binned analysis STAT = 'DEFAULT ' / Optimization statistic EDISP = F / Use energy dispersion? PROFILE = F / Use likelihood profile method for errors? TUNIT1 = ' ' / physical unit of field1 TUNIT4 = ' ' / physical unit of field4 TUNIT5 = ' ' / physical unit of field5 TUNIT6 = ' ' / physical unit of field6 TUNIT7 = ' ' / physical unit of field7 TUNIT8 = ' ' / physical unit of field8 TUNIT9 = ' ' / physical unit of field9 TUNIT10 = ' ' / physical unit of field10 TUNIT11 = ' ' / physical unit of field11 TUNIT12 = ' ' / physical unit of field12 TUNIT13 = ' ' / physical unit of field13 TUNIT14 = ' ' / physical unit of field14 TUNIT15 = ' ' / physical unit of field15 END Av*@{@{n< ǩ`J?lmx3? F(M@?&ްQ?膠??wYB?}U[E?xh;Agu.VY@B@A<ĝm3J,vwDB?}aT$&4"Aפb@W@@W?ѩ?㑬Φ5?K1D? a5?%Ҕh?l+tT?}Dn[r,?A_R@׏@@׏?vH<$=IV