Integrate code

Jürgen Knödlseder
1 parent 75bde4e8
Showing 5 changed files with 71 additions and 45 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.4.rst
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
-2017-07-03
+2017-07-28
         * Version 1.4.0 released
         ========================
         Install ctools header files (#2145)
+        Parallelize event binning in ctbin using OpenMP (#1970)
 2017-06-23
 New Features and Important Changes in ctools 1.4.0
-3 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
@@ -23,6 +23,8 @@ Improvements
 * [`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,25 +233,33 @@ 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()));
+    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()) {
-        std::string msg = " Including unbinned "+obs->instrument()+" observation";
-        log_string(NORMAL, msg);
+
+        // 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
-    #pragma omp parallel for
-    for (int oindx=0; oindx<obs_list.size(); ++oindx) {
-        GCTAObservation* obs = obs_list[oindx];
-        
+    #pragma omp parallel for(ctbin_run)
+    for (int i = 0; i < obs_list.size(); ++i) {
+
+        // Get pointer to observation
+        GCTAObservation* obs = obs_list[i];
+
         // Fill the cube
         fill_cube(obs);
@@ -371,10 +379,12 @@ void ctbin::init_members(void)
     m_ebounds.clear();
     m_gti.clear();
     m_cube.clear();
-    m_dirs.clear();
     m_ontime   = 0.0;
     m_livetime = 0.0;
+    // Initialise cache members
+    m_dirs.clear();
+
     // Set CTA reference time
     m_gti.reference(m_cta_ref);
@@ -397,14 +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_dirs      = app.m_dirs;
-    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;
@@ -445,13 +457,6 @@ void ctbin::get_parameters(void)
     m_counts  = 0.0;
     m_weights = m_counts;
-    // 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());
-    }
-    
     // Get energy boundaries
     m_ebounds = cube.ebounds();
@@ -465,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);
@@ -525,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()) {
@@ -564,32 +576,37 @@ void ctbin::fill_cube(GCTAObservation* obs)
         }
         // Fill event in skymap
-        #pragma omp critical
+        #pragma omp critical(ctbin_fill_cube)
         m_counts(pixel, iebin) += 1.0;
-        
+
+        // Increment number of maps
         num_in_map++;
     } // endfor: looped over all events
-    #pragma omp critical
+    // 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();
     }
-    
+
     // Log filling results
-    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);
-    
+    #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;
 }
@@ -636,7 +653,7 @@ void ctbin::set_weights(GCTAObservation* obs)
     // Loop over all pixels in counts cube
     for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
-        
+
         // Skip pixel if it is outside the RoI
         if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) {
             continue;
@@ -649,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,7 +102,9 @@ protected:
     GCTAEventCube m_cube;     //!< Events cube (for cube() method)
     double        m_ontime;   //!< Total ontime
     double        m_livetime; //!< Total livetime
-    std::vector<GSkyDir> m_dirs;  //!< Cached GSkyDir for each pixel in m_counts
+
+    // Cache members
+    std::vector<GSkyDir> m_dirs; //!< Cached GSkyDir for each pixel in m_counts
 };