From e8978f2ad7c05526b4e0c2537d8aad3397102efc Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 14 Aug 2016 02:35:25 +0200 Subject: [PATCH] Enhance coverage and comment density (#1828) --- cscripts/csiactobs.py | 4 ++++ cscripts/csinfo.py | 11 ++++------- cscripts/csobsdef.py | 2 +- cscripts/cspull.py | 58 +++++++++++++++++++++++++++------------------------------- cscripts/cssens.py | 109 +++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------ pyext/ctool.i | 21 +++++++++++---------- pyext/ctools/tools.i | 6 ++++++ test/test_cscript.py | 69 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------- test/test_csiactobs.py | 2 +- 9 files changed, 163 insertions(+), 119 deletions(-) diff --git a/cscripts/csiactobs.py b/cscripts/csiactobs.py index 1feb9b1..0972e35 100755 --- a/cscripts/csiactobs.py +++ b/cscripts/csiactobs.py @@ -78,6 +78,10 @@ class csiactobs(ctools.cscript): self._subdir = '' self._debug = False + # Initialise empty observation definition XML file + self._xml.append(gammalib.GXmlElement('observation_list ' + 'title="observation list"')) + # Initialise application by calling the appropriate class constructor self._init_cscript(argv) diff --git a/cscripts/csinfo.py b/cscripts/csinfo.py index 5c8852e..6e0281f 100755 --- a/cscripts/csinfo.py +++ b/cscripts/csinfo.py @@ -129,22 +129,19 @@ def csinfo_setup_check(): # Check if the "GAMMALIB" environment variable is set sys.stdout.write(' GAMMALIB environment variable ... ') - if 'GAMMALIB' in os.environ: + gammalib_environ = 'GAMMALIB' in os.environ + if gammalib_environ: print('ok') - gammalib_environ = True else: print('NOT OK') - gammalib_environ = False # Check if the "CTOOLS" environment variable is set sys.stdout.write(' CTOOLS environment variable ... ') - if 'CTOOLS' in os.environ: + ctools_environ = 'CTOOLS' in os.environ + if ctools_environ: print('ok') - ctools_environ = True else: print('NOT OK') - ctools_environ = False - # Check if gammalib Python module import works sys.stdout.write(' gammalib Python import .......... ') diff --git a/cscripts/csobsdef.py b/cscripts/csobsdef.py index 271fde1..29ae8c1 100755 --- a/cscripts/csobsdef.py +++ b/cscripts/csobsdef.py @@ -1,6 +1,6 @@ #! /usr/bin/env python # ========================================================================== -# Generation of an observation definition file. +# Generation of an observation definition file # # Copyright (C) 2015-2016 Juergen Knoedlseder # diff --git a/cscripts/cspull.py b/cscripts/cspull.py index 4739701..e77857f 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -236,31 +236,37 @@ class cspull(ctools.cscript): if self._logNormal(): self._log.header3('Pulls') - # Gather results - colnames = [] - values = {} - colnames.append('LogL') - colnames.append('Sim_Events') - colnames.append('Npred_Events') - values['LogL'] = logL - values['Sim_Events'] = nevents - values['Npred_Events'] = npred + # Gather results in form of a list of result columns and a + # dictionary containing the results. The result contains the + # log-likelihood, the number of simulated events, the number of + # predicted events and for each fitted parameter the fitted value, + # the pull and the fit error. + # + # Note that we do not use the model and parameter iterators + # because we need the indices to get the true (or real) parameter + # values from the input models. + colnames = ['LogL', 'Sim_Events', 'Npred_Events'] + values = {'LogL': logL, 'Sim_Events': nevents, 'Npred_Events': npred} for i in range(models.size()): - model = models[i] - model_name = model.name() + model = models[i] for k in range(model.size()): par = model[k] if par.is_free(): - # Set parameter name - name = model_name+'_'+par.name() + # 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() - # Append parameter, Pull_parameter and e_parameter + # Append parameter, Pull_parameter and e_parameter column + # names colnames.append(name) colnames.append('Pull_'+name) colnames.append('e_'+name) - # Compute pull + # Compute pull for this parameter as the difference + # (fitted - true) / error + # In case that the error is 0 the pull is set to 99 fitted_value = par.value() real_value = self._obs.models()[i][k].value() error = par.error() @@ -269,16 +275,16 @@ class cspull(ctools.cscript): else: pull = 99.0 - # Store results - values[name] = fitted_value + # Store results in dictionary + values[name] = fitted_value values['Pull_'+name] = pull values['e_'+name] = error - # Write result + # Write results into logger value = '%.4f (%e +/- %e)' % (pull, fitted_value, error) self._log_value(gammalib.NORMAL, name, value) - # Bundle together results + # Bundle together results in a dictionary result = {'colnames': colnames, 'values': values} # Return @@ -298,20 +304,10 @@ class cspull(ctools.cscript): self._get_parameters() # Write observation into logger - if self._logTerse(): - self._log('\n') - self._log.header1(gammalib.number('Observation',len(self._obs))) - self._log(str(self._obs)) - self._log('\n') - if self._logExplicit(): - for obs in self._obs: - self._log(str(obs)) - self._log('\n') + self._log_observations(gammalib.NORMAL, self._obs, 'Input observation') # Write header - if self._logTerse(): - self._log('\n') - self._log.header1('Generate pull distribution') + self._log_header1(gammalib.TERSE, 'Generate pull distribution') # Loop over trials for seed in range(self._ntrials): diff --git a/cscripts/cssens.py b/cscripts/cssens.py index b84bb6c..9d99909 100755 --- a/cscripts/cssens.py +++ b/cscripts/cssens.py @@ -73,8 +73,8 @@ class cssens(ctools.cscript): self._num_avg = 3 self._log_clients = False - # Initialise application by calling the appropriate class - # constructor. + # Initialise application by calling the appropriate base class + # constructor self._init_cscript(argv) # Return @@ -294,13 +294,13 @@ class cssens(ctools.cscript): e_mean = math.pow(10.0, loge) erg_mean = e_mean * tev2erg - # Compute Crab unit (this is the factor with which the Prefactor needs - # to be multiplied to get 1 Crab + # Compute Crab unit. This is the factor with which the Prefactor needs + # to be multiplied to get 1 Crab. crab_flux = self._get_crab_flux(emin, emax) src_flux = test_model[self._srcname].spectral().flux(emin, emax) crab_unit = crab_flux/src_flux - # Write header + # Write header and some introductory information into logger self._log_header2(gammalib.TERSE, 'Energies: '+str(emin)+' - '+str(emax)) self._log_value(gammalib.TERSE, 'Crab flux', str(crab_flux)+' ph/cm2/s') self._log_value(gammalib.TERSE, 'Source model flux', str(src_flux)+ @@ -321,37 +321,38 @@ class cssens(ctools.cscript): # Update iteration counter iterations += 1 - # Write header - if self._logExplicit(): - self._log.header2('Iteration '+str(iterations)) + # Write header for iteration into logger + self._log_header2(gammalib.EXPLICIT, 'Iteration '+str(iterations)) - # Set test source model. crab_prefactor is the Prefactor that - # corresponds to 1 Crab + # Create a copy of the test models, set the prefactor of the test + # source in the models, and append the models to the observation. + # "crab_prefactor" is the Prefactor that corresponds to a flux of + # 1 Crab. models = test_model.copy() - crab_prefactor = models[self._srcname]['Prefactor'].value() * \ - crab_unit - models[self._srcname]['Prefactor'].value(crab_prefactor * \ - test_crab_flux) + crab_prefactor = models[self._srcname]['Prefactor'].value() * crab_unit + models[self._srcname]['Prefactor'].value(crab_prefactor * test_crab_flux) self._obs.models(models) - # Simulate events + # Simulate events for the models. "sim" holds an observation + # container with observations containing the simulated events. sim = obsutils.sim(self._obs, nbins=self._enumbins, seed=iterations, binsz=self._binsz, npix=self._npix, log=self._log_clients, debug=self['debug'].boolean(), edisp=self._edisp) - # Determine number of events in simulation + # Determine number of events in simulation by summing the events + # over all observations in the observation container nevents = 0.0 for run in sim: nevents += run.events().number() - # Write simulation results + # Write simulation results into logger self._log_header3(gammalib.EXPLICIT, 'Simulation') - self._log_value(gammalib.EXPLICIT, 'Number of simulated events', - nevents) + self._log_value(gammalib.EXPLICIT, 'Number of simulated events', nevents) - # Fit test source + # Fit test source to the simulated events in the observation + # container fit = ctools.ctlike(sim) fit['edisp'] = edisp=self._edisp fit['debug'] = self['debug'].boolean() @@ -370,24 +371,28 @@ class cssens(ctools.cscript): photon_flux = source.spectral().flux(emin, emax) energy_flux = source.spectral().eflux(emin, emax) - # Compute differential sensitivity in unit erg/cm2/s + # Compute differential sensitivity in unit erg/cm2/s by evaluating + # the spectral model at the "e_mean" energy and by multipling the + # result with the energy squared. Since the "eval()" method returns + # an intensity in units of ph/cm2/s/MeV we multiply by 1.0e6 to + # convert into ph/cm2/s/TeV, by "e_mean" to convert into ph/cm2/s, + # and finally by "erg_mean" to convert to erg/cm2/s. energy = gammalib.GEnergy(e_mean, 'TeV') - sensitivity = source.spectral().eval(energy) * e_mean*erg_mean * 1.0e6 + sensitivity = source.spectral().eval(energy) * e_mean*erg_mean*1.0e6 # Assess quality based on a comparison between Npred and Nevents quality = npred - nevents - # Write test source fit results + # Write test source fit results into logger self._log_header3(gammalib.EXPLICIT, 'Test source model fit') self._log_value(gammalib.EXPLICIT, 'Test statistics', ts) self._log_value(gammalib.EXPLICIT, 'log likelihood', logL) - self._log_value(gammalib.EXPLICIT, 'Number of predicted events', - quality) + self._log_value(gammalib.EXPLICIT, 'Number of predicted events', quality) self._log_value(gammalib.EXPLICIT, 'Fit quality', npred) for model in models: self._log_value(gammalib.EXPLICIT, 'Model', model.name()) for par in model: - self._log(str(par)+'\n') + self._log_string(gammalib.EXPLICIT, str(par)) # If TS was non-positive then increase the test flux and start # over @@ -395,26 +400,27 @@ class cssens(ctools.cscript): # If the number of iterations was exceeded then stop if (iterations >= self._max_iter): - if self._logTerse(): - self._log(' Test ended after %d iterations.\n' % - self._max_iter) + self._log_string(gammalib.TERSE, + ' Test ended after %d iterations.' % self._max_iter) break # Increase test flux by a factor of 2 test_crab_flux = test_crab_flux * 2.0 # Signal start we start over - if self._logExplicit(): - self._log('Non positive TS, increase test flux and start ' - 'over.\n') + self._log_string(gammalib.EXPLICIT, + 'Non positive TS, increase test flux and start over.') # ... and start over continue - # Compute flux correction factor based on average TS + # Compute flux correction factor based on average TS. We + # assume here that the significance is proportional to the source + # flux, and that the significance scales as the square root of + # the Test Statistic correct = math.sqrt(self._ts_thres/ts) - # Compute extrapolated fluxes + # Compute extrapolated fluxes based on the flux correction factor crab_flux = correct * crab_flux photon_flux = correct * photon_flux energy_flux = correct * energy_flux @@ -437,7 +443,7 @@ class cssens(ctools.cscript): for model in models: self._log_value(gammalib.EXPLICIT, 'Model', model.name()) for par in model: - self._log(str(par)+'\n') + self._log_string(gammalib.EXPLICIT, str(par)) elif self._logTerse(): name = 'Iteration %d' % iterations value = ('TS=%.4f corr=%e %e ph/cm2/s = %e erg/cm2/s = %.2f ' @@ -447,7 +453,7 @@ class cssens(ctools.cscript): # Compute sliding average of extrapolated fitted prefactor, # photon and energy flux. This damps out fluctuations and - # improves convergence + # improves convergence. crab_flux = 0.0 photon_flux = 0.0 energy_flux = 0.0 @@ -494,7 +500,7 @@ class cssens(ctools.cscript): ' iterations.\n') break - # Write fit results + # Write fit results into logger self._log_header3(gammalib.TERSE, 'Fit results') self._log_value(gammalib.TERSE, 'Test statistics', ts) self._log_value(gammalib.TERSE, 'Photon flux', @@ -512,7 +518,7 @@ class cssens(ctools.cscript): for model in models: self._log_value(gammalib.TERSE, 'Model', model.name()) for par in model: - self._log_string(gammalib.TERSE, str(par)+'\n') + self._log_string(gammalib.TERSE, str(par)) # Restore energy boundaries of observation container for i, obs in enumerate(self._obs): @@ -553,31 +559,15 @@ class cssens(ctools.cscript): models = modutils.test_source(self._obs.models(), self._srcname, ra=self._ra, dec=self._dec) - # Write observations into logger - if self._logTerse(): - self._log('\n') - self._log.header1(gammalib.number('Observation',len(self._obs))) - self._log(str(self._obs)) - self._log('\n') - if self._logExplicit(): - for obs in self._obs: - self._log(str(obs)) - self._log('\n') + # Write observation into logger + self._log_observations(gammalib.NORMAL, self._obs, 'Observation') # Write models into logger - if self._logTerse(): - self._log('\n') - self._log.header1('Models') - self._log(str(models)) - self._log('\n') + self._log_models(gammalib.NORMAL, models, 'Model') # Write header - if self._logTerse(): - self._log('\n') - self._log.header1('Sensitivity determination') - self._log.parformat('Type') - self._log(self._type) - self._log('\n') + self._log_header1(gammalib.TERSE, 'Sensitivity determination') + self._log_value(gammalib.TERSE, 'Type', str(self._type)) # Loop over energy bins for ieng in range(self._ebounds.size()): @@ -620,7 +610,6 @@ class cssens(ctools.cscript): return - # ======================== # # Main routine entry point # # ======================== # diff --git a/pyext/ctool.i b/pyext/ctool.i index 076dc42..0d80e84 100644 --- a/pyext/ctool.i +++ b/pyext/ctool.i @@ -93,8 +93,6 @@ public: %rename(_create_cta_obs) create_cta_obs; %rename(_require_inobs) require_inobs; %rename(_require_inobs_nocube) require_inobs_nocube; - %rename(_log_observations) log_observations; - %rename(_log_models) log_models; %rename(_get_roi) get_roi; %rename(_set_response) set_response; %rename(_set_edisp) set_edisp; @@ -133,14 +131,6 @@ public: // Protected methods that extract user parameters GCTARoi get_roi(void); - // Protected methods for logging - void log_observations(const GChatter& chatter, - const GObservations& obs, - const std::string& what = "Observation"); - void log_models(const GChatter& chatter, - const GModels& models, - const std::string& what = "Model"); - // Protected support methods void set_response(GObservations& obs); std::vector set_edisp(GObservations& obs, const bool& edisp) const; @@ -193,12 +183,23 @@ public: * @brief ctool base class Python extensions ***************************************************************************/ %extend ctool { + void _log_observations(const int& chatter, + const GObservations& obs, + const std::string& what = "Observation") { + self->log_observations(GChatter(chatter), obs, what); + } + void _log_models(const int& chatter, + const GModels& models, + const std::string& what = "Model") { + self->log_models(GChatter(chatter), models, what); + } void _read_ahead(const bool& flag) { self->m_read_ahead = flag; return; } } + /***********************************************************************//** * @brief ctool and cscript base class Python extensions ***************************************************************************/ diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i index 56b7375..d5350be 100755 --- a/pyext/ctools/tools.i +++ b/pyext/ctools/tools.i @@ -36,6 +36,12 @@ /* __ Include standard typemaps for vectors and strings __________________ */ %include stl.i +%include "std_vector.i" +namespace std { + %template(vector_int) vector; + %template(vector_double) vector; + %template(vector_bool) vector; +}; /* __ Make sure that exceptions are catched ______________________________ */ %import(module="gammalib.support") "GException.i"; diff --git a/test/test_cscript.py b/test/test_cscript.py index 2bba959..08d0ae9 100755 --- a/test/test_cscript.py +++ b/test/test_cscript.py @@ -108,6 +108,23 @@ class cscript_test(ctools.cscript): obs = gammalib.GObservations() return self._get_mean_pointing(obs) + # Check set_obs_response() method + def check_set_obs_response(self, expcube, psfcube, bkgcube, edispcube, edisp): + obs = gammalib.GCTAObservation() + counts = gammalib.GCTAEventCube() + obs.events(counts) + self.pars()['expcube'].value(expcube) + self.pars()['psfcube'].value(psfcube) + self.pars()['bkgcube'].value(bkgcube) + self.pars()['edispcube'].value(edispcube) + self.pars()['edisp'].value(edisp) + self._set_obs_response(obs) + return + + # Check get_current_rss() method + def check_get_current_rss(self): + return self._get_current_rss() + # ======================== # # ctobservation_test class # @@ -196,6 +213,13 @@ class Test(test): pars.append(gammalib.GApplicationPar('rad','r','h','1.0','','','')) pars.append(gammalib.GApplicationPar('emin','r','h','1.0','','','')) pars.append(gammalib.GApplicationPar('emax','r','h','10.0','','','')) + pars.append(gammalib.GApplicationPar('expcube','f','h','NONE','','','')) + pars.append(gammalib.GApplicationPar('psfcube','f','h','NONE','','','')) + pars.append(gammalib.GApplicationPar('bkgcube','f','h','NONE','','','')) + pars.append(gammalib.GApplicationPar('edispcube','f','h','NONE','','','')) + pars.append(gammalib.GApplicationPar('edisp','b','h','yes','','','')) + pars.append(gammalib.GApplicationPar('caldb','s','h','prod2','','','')) + pars.append(gammalib.GApplicationPar('irf','s','h','South_0.5h','','','')) pars.append_standard() pars.save('cscript_test.par') @@ -333,6 +357,14 @@ class Test(test): except ValueError: self.test_try_success() + # Test require_inobs_nocube() method for valid file name + self.test_try('Check require_inobs_nocube() method for event list') + try: + empty.check_require_inobs_nocube(self._events) + self.test_try_success() + except ValueError: + self.test_try_failure('Exception thrown') + # Test require_inobs_nocube() method for invalid file name self.test_try('Check require_inobs_nocube() method for counts cube') try: @@ -346,15 +378,12 @@ class Test(test): self.test_assert(not roi.is_valid(), 'Check get_roi() method') # Test restore_edisp() method for invalid vector size - #self.test_try('Check restore_edisp() method for invalid vector size') - #try: - # empty.check_restore_edisp() - # self.test_try_failure('Exception not thrown') - #except ValueError: - # self.test_try_success() - # - # NOTE: This code does not work because of a missing list to - # integer vector conversion. + self.test_try('Check restore_edisp() method for invalid vector size') + try: + empty.check_restore_edisp() + self.test_try_failure('Exception not thrown') + except ValueError: + self.test_try_success() # Test set_obs_bounds() method obs = empty.check_set_obs_bounds() @@ -377,6 +406,28 @@ class Test(test): except ValueError: self.test_try_success() + # Test set_obs_response() method with 'NONE' files. If energy dispersion + # is requested then an exception will be thrown. + empty.check_set_obs_response('NONE', 'NONE', 'NONE', 'NONE', 'no') + empty.check_set_obs_response('NONE', 'NONE', 'NONE', 'NONE', 'yes') + self.test_try('Check set_obs_response() method for file names') + try: + empty.check_set_obs_response(self._datadir+'/crab_expcube.fits', + self._datadir+'/crab_psfcube.fits', + self._datadir+'/crab_bkgcube.fits', + 'NONE', 'yes') + self.test_try_failure('Exception not thrown') + except ValueError: + self.test_try_success() + empty.check_set_obs_response(self._datadir+'/crab_expcube.fits', + self._datadir+'/crab_psfcube.fits', + self._datadir+'/crab_bkgcube.fits', + self._datadir+'/crab_edispcube.fits', + 'yes') + + # Test get_current_rss() method + empty.check_get_current_rss() + # Return return diff --git a/test/test_csiactobs.py b/test/test_csiactobs.py index a4363eb..a55bfde 100755 --- a/test/test_csiactobs.py +++ b/test/test_csiactobs.py @@ -147,7 +147,7 @@ class Test(test): iactobs.save() # Check empty observation definition XML file - #self._check_obsdef('csiactobs_obs_py0.xml', 0) # See #1841 + self._check_obsdef('csiactobs_obs_py0.xml', 0) # Check empty model definition XML file self._check_moddef('csiactobs_bgd_py0.xml', 0) -- libgit2 0.21.4