Enhance coverage and comment density (#1828)

Jürgen Knödlseder
1 parent d9f9f2bd
Showing 9 changed files with 163 additions and 119 deletions
cscripts/csiactobs.py
cscripts/csinfo.py
cscripts/csobsdef.py
cscripts/cspull.py
cscripts/cssens.py
pyext/ctool.i
pyext/ctools/tools.i
test/test_cscript.py
test/test_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)
@@ -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 .......... ')
 #! /usr/bin/env python
 # ==========================================================================
-# Generation of an observation definition file.
+# Generation of an observation definition file
 #
 # Copyright (C) 2015-2016 Juergen Knoedlseder
 #
@@ -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):
@@ -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 #
 # ======================== #
@@ -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<bool> 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
  ***************************************************************************/
@@ -36,6 +36,12 @@
 /* __ Include standard typemaps for vectors and strings __________________ */
 %include stl.i
+%include "std_vector.i"
+namespace std {
+   %template(vector_int)    vector<int>;
+   %template(vector_double) vector<double>;
+   %template(vector_bool)   vector<bool>;
+};
 /* __ Make sure that exceptions are catched ______________________________ */
 %import(module="gammalib.support") "GException.i";
@@ -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
@@ -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)