Parallelise cslightcrv (#2421)

Jürgen Knödlseder
1 parent dcc6376a
Showing 8 changed files with 127 additions and 57 deletions
ChangeLog
NEWS
cscripts/cslightcrv.par
cscripts/cslightcrv.py
cscripts/mputils.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/cslightcrv.rst
pyext/ctool.i
-2018-06-29
+2018-07-03
         * Version 1.6.0 released
         ========================
@@ -7,6 +7,7 @@
         Change formula for number of true energy bins in csphagen (#2490)
         Remove ctools::is_valid_filename() method
         Fix bug when using precomputed model in csresspec (#2448)
+        Parallelise cslightcrv (#2421)
         Keep existing library path in installation script on Mac OS X (#2420)
         csphagen produces Xspec compliant data (#2404)
         Add show_rmf.py example script
 New Features and Important Changes in ctools 1.6.0
-29 June 2018
+3 July 2018
 Introduction
@@ -193,6 +193,9 @@ None
 cslightcrv - Generate light curve
 ---------------------------------
+Add hidden "nthreads" parameter and implement parallelisation of light
+curve computation (#2421).
+
 Add hidden "mjdref" parameter for the specification of the MET reference
 (#2054).
@@ -62,13 +62,13 @@ fix_bkg,   b, h, no,,, &quot;Fix background model parameters?&quot;
 #
 # Standard parameters
 #====================
-publish, b, h, no,,, "Publish light curve on VO Hub?"
-chatter, i, h, 2,0,4, "Chattiness of output"
-clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
-debug,   b, h, no,,, "Activate debugging mode?"
-mode,    s, h, ql,,, "Mode of automatic parameters"
-logfile, f, h, cslightcrv.log,,, "Log filename"
-nthreads,i, h, 0,0,1000, "Number of parallel processes. 0 = number of available CPUs."
+nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
+publish,  b, h, no,,, "Publish light curve on VO Hub?"
+chatter,  i, h, 2,0,4, "Chattiness of output"
+clobber,  b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,    b, h, no,,, "Activate debugging mode?"
+mode,     s, h, ql,,, "Mode of automatic parameters"
+logfile,  f, h, cslightcrv.log,,, "Log filename"
@@ -26,6 +26,33 @@ from cscripts import ioutils
 from cscripts import mputils
+# ============================================ #
+# Global functions for multiprocessing support #
+# ============================================ #
+def _multiprocessing_func_wrapper(args):
+   return _multiprocessing_func(*args)
+def _multiprocessing_func(cls, i):
+
+    # Initialise thread logger
+    cls._log.clear()
+    cls._log.buffer_size(100000)
+
+    # Compute light curve bin
+    cstart  = cls.celapse()
+    result  = cls._timebin(i)
+    celapse = cls.celapse() - cstart
+    buffer  = cls._log.buffer()
+
+    # Close logger
+    cls._log.close()
+
+    # Collect thread information
+    info = {'celapse': celapse, 'log': buffer}
+
+    # Return light curve bin result and thread information
+    return result, info
+
+
 # ================ #
 # cslightcrv class #
 # ================ #
@@ -73,11 +100,12 @@ class cslightcrv(ctools.csobservation):
         self._init_csobservation(self.__class__.__name__, ctools.__version__, argv)
         # Initialise some members
-        self._srcname = ''
-        self._tbins   = gammalib.GGti()
-        self._onoff   = False
-        self._stacked = False
-        self._fits    = gammalib.GFits()
+        self._srcname  = ''
+        self._tbins    = gammalib.GGti()
+        self._onoff    = False
+        self._stacked  = False
+        self._fits     = gammalib.GFits()
+        self._nthreads = 0
         # Return
         return
@@ -87,31 +115,42 @@ class cslightcrv(ctools.csobservation):
         """
         Extend ctools.csobservation getstate method to include some members
-        :return: dictionary
-            state
+        Returns
+        -------
+        state : dict
+            Pickled instance
         """
-
+        # Set pickled dictionary
         state = {'base'     : ctools.csobservation.__getstate__(self),
                  'srcname'  : self._srcname,
                  'tbins'    : self._tbins,
                  'stacked'  : self._stacked,
                  'onoff'    : self._onoff,
-                 'fits'     : self._fits}
+                 'fits'     : self._fits,
+                 'nthreads' : self._nthreads}
+
+        # Return pickled dictionary
         return state
     def __setstate__(self, state):
         """
         Extend ctools.csobservation setstate method to include some members
-        :param state: dictionary of class state
+        Parameters
+        ----------
+        state : dict
+            Pickled instance
         """
-
         ctools.csobservation.__setstate__(self, state['base'])
-        self._srcname   = state['srcname']
-        self._tbins     = state['tbins']
-        self._onoff     = state['onoff']
-        self._stacked   = state['stacked']
-        self._fits      = state['fits']
+        self._srcname  = state['srcname']
+        self._tbins    = state['tbins']
+        self._onoff    = state['onoff']
+        self._stacked  = state['stacked']
+        self._fits     = state['fits']
+        self._nthreads = state['nthreads']
+
+        # Return
+        return
     # Private methods
     def _get_parameters(self):
@@ -399,16 +438,20 @@ class cslightcrv(ctools.csobservation):
         # Return upper limit tuple
         return ul_diff, ul_flux, ul_eflux
-
-    def _timebin(self,i):
+    def _timebin(self, i):
         """
-        Run likelihood analysis in one time bin.
+        Run likelihood analysis in one time bin
-        :param i: integer, time bin number
-        :return: dictionary,
-                results of the likelihood analysis
-        """
+        Parameters
+        ----------
+        i : int
+            time bin number
+        Returns
+        -------
+        result : dict
+            Results of the likelihood analysis
+        """
         # Get names of free parameters
         pars = self._get_free_par_names()
@@ -611,8 +654,19 @@ class cslightcrv(ctools.csobservation):
             pool = Pool(processes = self._nthreads)
             # Run time bin analysis in parallel with map
-            results = pool.map(self._timebin,range(self._tbins.size()))
-            
+            args        = [(self, i) for i in range(self._tbins.size())]
+            poolresults = pool.map(_multiprocessing_func_wrapper, args)
+
+            # Close pool and join
+            pool.close()
+            pool.join()
+
+            # Construct results
+            results = []
+            for i in range(self._tbins.size()):
+                results.append(poolresults[i][0])
+                self._log_string(gammalib.TERSE, poolresults[i][1]['log'], False)
+
         # Otherwise loop over time bins and run time bin analysis
         else:
             results = []
 # ==========================================================================
 # Utility functions for multiprocessing handling
 #
-# Copyright (C) 2018- Luigi Tibaldo
+# Copyright (C) 2018 Luigi Tibaldo
 #
 # 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
@@ -23,44 +23,50 @@ import gammalib
 # ==================== #
 # Number of processes  #
 # ==================== #
-
 def nthreads(cls):
     """
-    Determines the number of parallel processes to use,
-    based on user parameters and availability of multiprocessing module
+    Determines the number of parallel processes to use.
+
+    The number is based on the user parameter "nthreads" and the availability
+    of the multiprocessing module.
-    :param cls: `~ctools.cscript`
+    Parameters
+    ----------
+    cls : `~ctools.cscript`
         cscript class
-    :return: int
-        number of parallel processes
-    """
+    Returns
+    -------
+    nthreads : int
+        Number of parallel processes
+    """
     # Log multiprocessing configuration
-    cls._log_header1(gammalib.NORMAL, "Multiprocessing")
+    cls._log_header1(gammalib.NORMAL, 'Multiprocessing')
+    # Try getting multiprocessing support
     try:
         from multiprocessing import cpu_count
         ncpus = cpu_count()
+        cls._log_value(gammalib.NORMAL, 'Multiprocessing', 'available')
         cls._log_value(gammalib.EXPLICIT, 'Number of CPUs available', ncpus)
         # Set processes to number of CPUs if requested by the user
         if cls['nthreads'].integer() == 0:
-            nthr = ncpus
+            nthreads = ncpus
-        # Otherwise use the number requested
+        # ... otherwise use the number requested
         else:
-            cls._log_value(gammalib.EXPLICIT, 'Number of processes requested by the user',
+            cls._log_value(gammalib.EXPLICIT,
+                           'Number of processes requested',
                            cls['nthreads'].integer())
-            nthr = cls['nthreads'].integer()
+            nthreads = cls['nthreads'].integer()
     except:
-        nthr = 1
-        cls._log_value(gammalib.EXPLICIT,
-                       'Multiprocessing not available. Number of processes', nthr)
-
+        nthreads = 1
+        cls._log_value(gammalib.NORMAL, 'Multiprocessing', 'not available')
     # Log number of processes
-    cls._log_value(gammalib.NORMAL, 'Number of parallel processes in use', nthr)
+    cls._log_value(gammalib.NORMAL, 'Number of processes in use', nthreads)
-    # return number of processes to be used
-    return nthr
+    # Return number of processes to be used
+    return nthreads
@@ -12,6 +12,7 @@ In particular, this release provides:
 * Support for Xspec spectral analysis
 * Improved computation speed for :ref:`ctskymap` and :ref:`csphagen`
+* Parallelisation of :ref:`cslightcrv`
 * Support for Python pickeling
@@ -31,6 +32,8 @@ Improvements
   Transparently update parameter files
 * [`2490 <https://cta-redmine.irap.omp.eu/issues/2490>`_] -
   Change formula for number of true energy bins in :ref:`csphagen`
+* [`2421 <https://cta-redmine.irap.omp.eu/issues/2421>`_] -
+  Parallelise :ref:`cslightcrv`
 * [`2404 <https://cta-redmine.irap.omp.eu/issues/2404>`_] -
   :ref:`csphagen` produces Xspec compliant data
 * Add ``show_rmf.py`` example script
@@ -177,6 +177,9 @@ General parameters
 Standard parameters
 -------------------
+``(nthreads = 0) [integer]``
+    Number of parallel processes (0=use all available CPUs).
+
 ``(publish = no) [boolean]``
     Specifies whether the light curve should be published on VO Hub.
@@ -222,9 +222,9 @@ public:
                            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") {
+    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) {