Merge branch 'integration' into devel

Jürgen Knödlseder
2 parents 624999d9 36b7c210
Showing 5 changed files with 82 additions and 30 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.4.rst
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
-2017-07-27
+2017-07-28
         * Version 1.4.0 released
         ========================
         Fix ctmodel bug for stacked analysis with diffuse maps or cubes (#2171)
         Install ctools header files (#2145)
+        Parallelize event binning in ctbin using OpenMP (#1970)
 2017-07-07
 New Features and Important Changes in ctools 1.4.0
-27 July 2017
+28 July 2017
 Introduction
@@ -37,7 +37,9 @@ None
 ctbin - CTA event binning
 -------------------------
-None
+The loop over the observations was parallelized using OpenMP for faster
+execution. In addition, the sky directions were cached for faster computation
+of the weights (#1970).
 ctlike - Maximum likelihood fitting
@@ -25,6 +25,8 @@ Improvements
   Fix ctmodel bug for stacked analysis with diffuse maps or cubes
 * [`2145 <https://cta-redmine.irap.omp.eu/issues/2145>`_] -
   Install C++ header files
+* [`1970 <https://cta-redmine.irap.omp.eu/issues/1970>`_] -
+  Parallelize :ref:`ctbin` using OpenMP
 New features
@@ -233,11 +233,32 @@ void ctbin::run(void)
     log_observations(NORMAL, m_obs, "Input observation");
     // Write header into logger
+    log_header1(TERSE, gammalib::number("Find unbinned observation",
+                m_obs.size()));
+
+    // Find all unbinned CTA observations in m_obs
+    std::vector<GCTAObservation*> obs_list(0);
+    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
+            obs = next_unbinned_observation()) {
+
+        // Push observation into list
+        obs_list.push_back(obs);
+
+        // Write message
+        std::string msg = " Including unbinned "+obs->instrument()+
+                          " observation";
+        log_string(NORMAL, msg);
+    }
+
+    // Write header into logger
     log_header1(TERSE, gammalib::number("Bin observation", m_obs.size()));
     // Loop over all unbinned CTA observations in the container
-    for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
-         obs = next_unbinned_observation()) {
+    #pragma omp parallel for
+    for (int i = 0; i < obs_list.size(); ++i) {
+
+        // Get pointer to observation
+        GCTAObservation* obs = obs_list[i];
         // Fill the cube
         fill_cube(obs);
@@ -361,6 +382,9 @@ void ctbin::init_members(void)
     m_ontime   = 0.0;
     m_livetime = 0.0;
+    // Initialise cache members
+    m_dirs.clear();
+
     // Set CTA reference time
     m_gti.reference(m_cta_ref);
@@ -383,13 +407,16 @@ void ctbin::copy_members(const ctbin&amp; app)
     m_chatter = app.m_chatter;
     // Copy protected members
-    m_counts    = app.m_counts;
-    m_weights   = app.m_weights;
-    m_ebounds   = app.m_ebounds;
-    m_gti       = app.m_gti;
-    m_cube      = app.m_cube;
-    m_ontime    = app.m_ontime;
-    m_livetime  = app.m_livetime;
+    m_counts   = app.m_counts;
+    m_weights  = app.m_weights;
+    m_ebounds  = app.m_ebounds;
+    m_gti      = app.m_gti;
+    m_cube     = app.m_cube;
+    m_ontime   = app.m_ontime;
+    m_livetime = app.m_livetime;
+
+    // Copy cache members
+    m_dirs = app.m_dirs;
     // Return
     return;
@@ -443,6 +470,13 @@ void ctbin::get_parameters(void)
         m_outcube = (*this)["outcube"].filename();
     }
+    // Cache sky directions in m_counts
+    m_dirs.resize(m_counts.npix());
+    for (int pix = 0; pix < m_counts.npix(); ++pix) {
+        m_dirs[pix] = m_counts.inx2dir(pix);
+        m_dirs[pix].radec(m_dirs[pix].ra(), m_dirs[pix].dec());
+    }
+
     // Write parameters into logger
     log_parameters(TERSE);
@@ -503,8 +537,8 @@ void ctbin::fill_cube(GCTAObservation* obs)
         const GCTAEventAtom* event = (*events)[i];
         // Determine event sky direction
-        GCTAInstDir* inst  = (GCTAInstDir*)&(event->dir());
-        GSkyDir      dir   = inst->dir();
+        GCTAInstDir* inst = (GCTAInstDir*)&(event->dir());
+        GSkyDir      dir  = inst->dir();
         // Skip event if it is outside the RoI
         if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
@@ -542,25 +576,36 @@ void ctbin::fill_cube(GCTAObservation* obs)
         }
         // Fill event in skymap
+        #pragma omp critical(ctbin_fill_cube)
         m_counts(pixel, iebin) += 1.0;
+
+        // Increment number of maps
         num_in_map++;
     } // endfor: looped over all events
-    // Append GTIs
-    m_gti.extend(events->gti());
+    // Update time information
+    #pragma omp critical(ctbin_fill_cube)
+    {
+        // Append GTIs
+        m_gti.extend(events->gti());
-    // Update ontime and livetime
-    m_ontime   += obs->ontime();
-    m_livetime += obs->livetime();
+        // Update ontime and livetime
+        m_ontime   += obs->ontime();
+        m_livetime += obs->livetime();
+    }
     // Log filling results
-    log_value(NORMAL, "Events in list", obs->events()->size());
-    log_value(NORMAL, "Events in cube", num_in_map);
-    log_value(NORMAL, "Events outside RoI", num_outside_roi);
-    log_value(NORMAL, "Events with invalid WCS", num_invalid_wcs);
-    log_value(NORMAL, "Events outside cube area", num_outside_map);
-    log_value(NORMAL, "Events outside energy bins", num_outside_ebds);
+    #pragma omp critical(ctbin_fill_cube)
+    {
+        log_header3(TERSE, get_obs_header(obs));
+        log_value(NORMAL, "Events in list", obs->events()->size());
+        log_value(NORMAL, "Events in cube", num_in_map);
+        log_value(NORMAL, "Events outside RoI", num_outside_roi);
+        log_value(NORMAL, "Events with invalid WCS", num_invalid_wcs);
+        log_value(NORMAL, "Events outside cube area", num_outside_map);
+        log_value(NORMAL, "Events outside energy bins", num_outside_ebds);
+    }
     // Return
     return;
@@ -609,11 +654,8 @@ void ctbin::set_weights(GCTAObservation* obs)
     // Loop over all pixels in counts cube
     for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
-        // Get pixel sky direction
-        GSkyDir dir = m_counts.inx2dir(pixel);
-
         // Skip pixel if it is outside the RoI
-        if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
+        if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) {
             continue;
         }
@@ -624,7 +666,9 @@ void ctbin::set_weights(GCTAObservation* obs)
             if (!usage[iebin]) {
                 continue;
             }
+
             // Signal that bin was filled
+            #pragma omp critical(ctbin_set_weights)
             m_weights(pixel, iebin) = 1.0;
         } // endfor: looped over energy layers of counts cube
 /***************************************************************************
  *                        ctbin - Event binning tool                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -102,6 +102,9 @@ protected:
     GCTAEventCube m_cube;     //!< Events cube (for cube() method)
     double        m_ontime;   //!< Total ontime
     double        m_livetime; //!< Total livetime
+
+    // Cache members
+    std::vector<GSkyDir> m_dirs; //!< Cached GSkyDir for each pixel in m_counts
 };