Add nthreads parameter to ctlike and cterror (#2421)

Jürgen Knödlseder
1 parent 0ac251c1
Showing 17 changed files with 83 additions and 33 deletions
ChangeLog
NEWS
cscripts/cslightcrv.py
cscripts/csphasecrv.py
cscripts/cspull.py
cscripts/cssens.py
cscripts/csspec.py
cscripts/cstsdist.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/cslightcrv.rst
doc/source/users/reference_manual/csphasecrv.rst
doc/source/users/reference_manual/cterror.rst
doc/source/users/reference_manual/ctlike.rst
src/cterror/cterror.cpp
src/cterror/cterror.par
src/ctlike/ctlike.cpp
src/ctlike/ctlike.par
@@ -10,7 +10,7 @@
         Add Exponential spectral model to science verification (#2454)
         Fix bug when using precomputed model in csresspec (#2448)
         Parallelise cslightcrv, csphasecrv, csspec, cssens, csphagen, cspull and
-          cstsdist (#2421)
+          cstsdist and add nthreads parameter to ctlike and cterror (#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
@@ -74,7 +74,7 @@ None
 cterror - Likelihood profile error computation
 ----------------------------------------------
-None
+Add hidden "nthreads" parameter to control parallel computations (#2421).
 ctexpcube - CTA exposure cube generation
@@ -84,7 +84,7 @@ None
 ctlike - Maximum likelihood fitting
 -----------------------------------
-None
+Add hidden "nthreads" parameter to control parallel computations (#2421).
 ctmapcube - CTA map cube generation
@@ -530,7 +530,8 @@ class cslightcrv(ctools.csobservation):
         # Do maximum likelihood model fitting
         if obs.size() > 0:
             like = ctools.ctlike(obs)
-            like['edisp'] = self['edisp'].boolean()
+            like['edisp']    = self['edisp'].boolean()
+            like['nthreads'] = 1  # Avoids OpenMP conflict
             like.run()
             # Skip bin if no event was present
@@ -443,7 +443,8 @@ class csphasecrv(ctools.csobservation):
             # Do maximum likelihood model fitting
             like = ctools.ctlike(obs)
-            like['edisp'] = self['edisp'].boolean()
+            like['edisp']    = self['edisp'].boolean()
+            like['nthreads'] = 1  # Avoids OpenMP conflict
             like.run()
             # Renormalize models to phase selection
@@ -464,11 +465,13 @@ class csphasecrv(ctools.csobservation):
             # Set empty models container
             fitmodels = gammalib.GModels()
-        result = {'phstr'       :phstr,
-                  'fitmodels'   :fitmodels}
+        # Set results
+        result = {'phstr': phstr, 'fitmodels': fitmodels}
+        # Return results
         return result
+
     # Public methods
     def run(self):
         """
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This script generates the pull distribution for all model parameters.
 #
-# Copyright (C) 2011-2017 Juergen Knoedlseder
+# Copyright (C) 2011-2018 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
@@ -25,6 +25,7 @@ from cscripts import obsutils
 from cscripts import ioutils
 from cscripts import mputils
+
 # ============================================ #
 # Global functions for multiprocessing support #
 # ============================================ #
@@ -282,6 +283,7 @@ class cspull(ctools.csobservation):
                 like['srcname']   = model.name()
                 like['edisp']     = edisp
                 like['statistic'] = statistic
+                like['nthreads']  = 1  # Avoids OpenMP conflict
                 like['debug']     = self._logDebug()
                 like['chatter']   = self['chatter'].integer()
                 like.run()
@@ -289,6 +291,7 @@ class cspull(ctools.csobservation):
             like = ctools.ctlike(obs)
             like['edisp']     = edisp
             like['statistic'] = statistic
+            like['nthreads']  = 1  # Avoids OpenMP conflict
             like['debug']     = self._logDebug()
             like['chatter']   = self['chatter'].integer()
             like.run()
@@ -27,6 +27,7 @@ from cscripts import modutils
 from cscripts import ioutils
 from cscripts import mputils
+
 # ============================================ #
 # Global functions for multiprocessing support #
 # ============================================ #
@@ -434,9 +435,10 @@ class cssens(ctools.csobservation):
             # Fit test source to the simulated events in the observation
             # container
             fit = ctools.ctlike(sim)
-            fit['edisp']   = self['edisp'].boolean()
-            fit['debug']   = self['debug'].boolean()
-            fit['chatter'] = self['chatter'].integer()
+            fit['edisp']    = self['edisp'].boolean()
+            fit['nthreads'] = 1  # Avoids OpenMP conflict
+            fit['debug']    = self['debug'].boolean()
+            fit['chatter']  = self['chatter'].integer()
             fit.run()
             # Get model fitting results
@@ -23,6 +23,7 @@ import gammalib
 import ctools
 from cscripts import mputils
+
 # ============================================ #
 # Global functions for multiprocessing support #
 # ============================================ #
@@ -674,8 +675,9 @@ class csspec(ctools.csobservation):
         self._log_header3(gammalib.EXPLICIT, 'Performing fit')
         # Perform maximum likelihood fit. Skip bins that have no events.
-        like          = ctools.ctlike(obs)
-        like['edisp'] = self['edisp'].boolean()
+        like = ctools.ctlike(obs)
+        like['edisp']    = self['edisp'].boolean()
+        like['nthreads'] = 1  # Avoids OpenMP conflict
         if like.obs().nobserved() > 0:
             like.run()
@@ -26,6 +26,7 @@ from cscripts import modutils
 from cscripts import ioutils
 from cscripts import mputils
+
 # ============================================ #
 # Global functions for multiprocessing support #
 # ============================================ #
@@ -235,8 +236,9 @@ class cstsdist(ctools.csobservation):
         # Fit model
         fit = ctools.ctlike(sim)
-        fit['debug']   = self['debug'].boolean()
-        fit['chatter'] = self['chatter'].integer()
+        fit['nthreads'] = 1  # Avoids OpenMP conflict
+        fit['debug']    = self['debug'].boolean()
+        fit['chatter']  = self['chatter'].integer()
         fit.run()
         # Get model fitting results
@@ -34,6 +34,9 @@ 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>`_] -
+  Add `nthreads` parameter to :ref:`ctlike` and :ref:`cterror` to control
+  parallel computations
 * [`2404 <https://cta-redmine.irap.omp.eu/issues/2404>`_] -
   :ref:`csphagen` produces Xspec compliant data
 * Add ``show_rmf.py`` example script
@@ -31,10 +31,10 @@ are necessary to run internally the :ref:`csphagen` script (``inexclusion``, ``c
 ``xref``, ``yref``, ``srcshape``, ``rad``, ``bkgmethod``, ``bkgregmin``, ``maxoffset``,
 ``etruemin``, ``etruemax``, ``etruebins``).
-:ref:`cslightcrv` supports multiprocessing. By default the analysis in each time bin will
-be performed in parallel over as many processes as the number of CPUs available on your
-machine. The maximum number of parallel processes can be set by the user through the
-``nthreads`` hidden parameter.
+:ref:`cslightcrv` supports multiprocessing. By default the analysis in each time bin
+will be performed in parallel over as many processes as the number of CPUs
+available on your machine. The maximum number of parallel processes can be set
+by the user through the ``nthreads`` hidden parameter.
 :ref:`cslightcrv` writes the fitted model parameters and their statistical errors
 in a FITS file. In addition, it computes for each time bin the statistical 
@@ -31,10 +31,10 @@ are necessary to run internally the :ref:`csphagen` script (``inexclusion``, ``c
 ``xref``, ``yref``, ``srcshape``, ``rad``, ``bkgmethod``, ``bkgregmin``, ``maxoffset``,
 ``etruemin``, ``etruemax``, ``etruebins``).
-:ref:`csphasecrv` supports multiprocessing. By default the analysis in each phase bin will
-be performed in parallel over as many processes as the number of CPUs available on your
-machine. The maximum number of parallel processes can be set by the user through the
-``nthreads`` hidden parameter.
+:ref:`csphasecrv` supports multiprocessing. By default the analysis in each phase bin
+will be performed in parallel over as many processes as the number of CPUs
+available on your machine. The maximum number of parallel processes can be set
+by the user through the ``nthreads`` hidden parameter.
 The script writes the fit results into a FITS file. The script also produces one
 XML file per phase bin that contains the best-fit model for that bin.
@@ -83,6 +83,9 @@ General parameters
 Standard parameters
 -------------------
+``(nthreads = 0) [integer]``
+    Number of parallel processes (0=use all available CPUs).
+
 ``(chatter = 2) [integer]``
     Verbosity of the executable:
      ``chatter = 0``: no information will be logged
@@ -112,6 +112,9 @@ General parameters
 Standard parameters
 -------------------
+``(nthreads = 0) [integer]``
+    Number of parallel processes (0=use all available CPUs).
+
 ``(chatter = 2) [integer]``
     Verbosity of the executable:
      ``chatter = 0``: no information will be logged
@@ -33,6 +33,11 @@
 #include "GTools.hpp"
 #include "GOptimizer.hpp"
+/* __ OpenMP section _____________________________________________________ */
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 /* __ Method name definitions ____________________________________________ */
 #define G_ERR_BISECTION         "cterror::error_bisection(double&, double&)"
@@ -477,6 +482,14 @@ void cterror::get_parameters(void)
         (*this)["outmodel"].query();
     }
+    // Set number of OpenMP threads
+    #ifdef _OPENMP
+    int nthreads = (*this)["nthreads"].integer();
+    if (nthreads > 0) {
+        omp_set_num_threads(nthreads);
+    }
+    #endif
+
     // Write parameters into logger
     log_parameters(TERSE);
@@ -40,8 +40,9 @@ max_iter,   i, h, 50,1,1000, &quot;Maximum number of iterations&quot;
 #
 # Standard parameters
 #====================
-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, cterror.log,,, "Log filename"
+nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
+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, cterror.log,,, "Log filename"
@@ -32,6 +32,11 @@
 #include "ctlike.hpp"
 #include "GTools.hpp"
+/* __ OpenMP section _____________________________________________________ */
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
 /* __ Method name definitions ____________________________________________ */
 /* __ Debug definitions __________________________________________________ */
@@ -422,6 +427,14 @@ void ctlike::get_parameters(void)
         static_cast<GOptimizerLM*>(&m_opt)->logger(NULL);
     }
+    // Set number of OpenMP threads
+    #ifdef _OPENMP
+    int nthreads = (*this)["nthreads"].integer();
+    if (nthreads > 0) {
+        omp_set_num_threads(nthreads);
+    }
+    #endif
+
     // Write parameters into logger
     log_parameters(TERSE);
@@ -38,8 +38,9 @@ fix_spat_for_ts, b, h, no,,, &quot;Fix spatial parameters for TS computation?&quot;
 #
 # Standard parameters
 #====================
-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, ctlike.log,,, "Log filename"
+nthreads, i, h, 0,0,1000, "Number of parallel processes (0=use all available CPUs)"
+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, ctlike.log,,, "Log filename"