Updated model trimming and parallelization in ctmodel (#2189)

Updated the model trimming so that the models that are outside an observation are offloaded to another model container and then re-appended after the evaluation. Parallelization scheduling set to "static" with a chunk size of 1. This will slightly help with the issue in which multiple observations in a row do not contribute to the model cube and thus dividing the observations into "n-observations/n-cores" is not ideal.

Updated model trimming and parallelization in ctmodel (#2189)
Updated the model trimming so that the models that are outside an observation are offloaded to another model container and then re-appended after the evaluation. Parallelization scheduling set to "static" with a chunk size of 1. This will slightly help with the issue in which multiple observations in a row do not contribute to the model cube and thus dividing the observations into "n-observations/n-cores" is not ideal.
Josh Cardenzana
1 parent 1c4c3571
Showing 2 changed files with 71 additions and 61 deletions
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
@@ -235,48 +235,55 @@ void ctmodel::run(void)
     log_header1(TERSE, "Generate model cube");
  
     // Loop over all observations in the container
-    #pragma omp parallel for
-    for (int i = 0; i < m_obs.size(); ++i) {
+    #pragma omp parallel
+    {
+        // Create a models object for this core
+        GModels models(m_obs.models());
  
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Skip observation if it's not CTA
-        if (obs == NULL) {
-            log_header3(TERSE, get_obs_header(m_obs[i]));
-            std::string msg = " Skipping "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-        }
-
-        // Fill cube and leave loop if we are binned mode (meaning we
-        // only have one binned observation)
-        else if (m_binned) {
-            fill_cube(obs);
-            i = m_obs.size();
-        }
-
-        // Skip observation if we have a binned observation
-        else if (obs->eventtype() == "CountsCube") {
-            log_header3(TERSE, get_obs_header(m_obs[i]));
-            std::string msg = " Skipping binned "+m_obs[i]->instrument()+
-                              " observation";
-            log_string(NORMAL, msg);
-        }
-
-        // Otherwise, everything seems to be fine, so fill the cube from obs
-        else {
-            // Fill the cube
-            fill_cube(obs);
+        #pragma omp for schedule(static,1)
+        for (int i = 0; i < m_obs.size(); ++i) {
+            
+            // Get CTA observation
+            GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
-            // Dispose events to free memory if event file exists on disk
-            if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
-                obs->dispose_events();
+            // Skip observation if it's not CTA
+            if (obs == NULL) {
+                log_header3(TERSE, get_obs_header(m_obs[i]));
+                std::string msg = " Skipping "+m_obs[i]->instrument()+
+                " observation";
+                log_string(NORMAL, msg);
             }
-        }
-
-    } // endfor: looped over observations
+            
+            // Fill cube and leave loop if we are binned mode (meaning we
+            // only have one binned observation)
+            else if (m_binned) {
+                fill_cube(obs, models);
+                i = m_obs.size();
+            }
+            
+            // Skip observation if we have a binned observation
+            else if (obs->eventtype() == "CountsCube") {
+                log_header3(TERSE, get_obs_header(m_obs[i]));
+                std::string msg = " Skipping binned "+m_obs[i]->instrument()+
+                " observation";
+                log_string(NORMAL, msg);
+            }
+            
+            // Otherwise, everything seems to be fine, so fill the cube from obs
+            else {
+                // Fill the cube
+                fill_cube(obs, models);
+                
+                // Dispose events to free memory if event file exists on disk
+                if (obs->eventfile().length() > 0 && obs->eventfile().exists()) {
+                    obs->dispose_events();
+                }
+            }
+            
+        } // endfor: looped over observations
  
+    } // end openmp parallel section
+    
     // Write model cube into header
     log_header1(NORMAL, "Model cube");
     log_string(NORMAL, m_cube.print());
@@ -789,7 +796,7 @@ void ctmodel::extract_cube_properties(void)
  * GTI of the model cube so that cube GTI is a list of the GTIs of all
  * observations that were used to generate the model cube.
  ***************************************************************************/
-void ctmodel::fill_cube(const GCTAObservation* obs)
+void ctmodel::fill_cube(const GCTAObservation* obs, GModels& models)
 {
     // Get number of spatial and spectral bins in counts cube
     int npix   = m_cube.npix();
@@ -819,7 +826,7 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     int    num_outside_roi  = 0;
  
     // Extract copy of models from observation container
-    GModels models = trim_models(m_obs.models(), roi);
+    GModels ex_models = trim_models(models, roi);
  
     // Get pointer to event cube pixels
     double* pixels = const_cast<double*>(m_cube.counts().pixels());
@@ -877,6 +884,9 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
  
     } // endfor: looped over all cube bins
  
+    // Re-append the models to the main model
+    models.extend(ex_models);
+    
     // Update GTIs
     #pragma omp critical(ctmodel_fill_cube)
     {
@@ -918,27 +928,24 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
  * point sources (which have radius=0) are appropriately included if they
  * fall near the edge of the observation.
  ***************************************************************************/
-GModels ctmodel::trim_models(const GModels& all_models, const GCTARoi& roi)
+GModels ctmodel::trim_models(GModels& all_models, const GCTARoi& roi)
 {
+    // Create the model object to be returned
+    GModels excluded_models;
+    
     // Do nothing if roi is not valid. This will be the case for binned data
     // sets as they have the default roi radius of 0. In this case we would
     // remove ALL models, which is obviously incorrect. So we keep all models
     // when filling based on a counts cube.
     if (!roi.is_valid()) {
-        return all_models;
+        return excluded_models;
     }
  
     // Remove all models that dont overlap with the region of interest. Note
     // that an extra factor is used since point sources have regions of
     // radius 0, which is a problem when they fall just barely outside the
     // ROI. This ensures point sources are appropriately included.
-    GSkyRegionCircle obsreg(
-        roi.centre().dir().ra_deg(),
-        roi.centre().dir().dec_deg(),
-        roi.radius()+0.5);                  // Note the buffer
-
-    // Create the model object to be returned
-    GModels contained_models;
+    GSkyRegionCircle obsreg(roi.centre().dir(), roi.radius()+0.5);
  
     // Loop over the models in the passed container
     for (int mod=0; mod<all_models.size(); mod++) {
@@ -947,18 +954,21 @@ GModels ctmodel::trim_models(const GModels&amp; all_models, const GCTARoi&amp; roi)
  
         if (model == NULL) {
             // Model is most likely a background model and should be kept
-            contained_models.append(*all_models.at(mod)->clone());
-            
-        } else if (model->spatial()->region()->overlaps(obsreg)) {
-            // Model represents a source, that overlaps the observation
-            // and should be kept.
-            contained_models.append(*model);
-            
-            // Delete model since it isnt needed anymore
-            delete model;
+            continue;
+        }
+        // Otherwise, exclude the model if it isn't a diffuse cube and it
+        // doesn't overlap with the observation.
+        else if (!(model->spatial()->code() == GMODEL_SPATIAL_DIFFUSE) &&
+            (!model->spatial()->region()->overlaps(obsreg))) {
+            // Model represents a source, that doesn't overlap the observation
+            // and should be removed for now.
+            excluded_models.append(*model);
+            all_models.remove(mod);
+            mod--;  // decrement to prevent skipping models
         }
  
+        delete model;
     }
  
-    return contained_models;
+    return excluded_models;
 }
@@ -80,9 +80,9 @@ protected:
     void get_parameters(void);
     void get_obs(void);
     void extract_cube_properties(void);
-    void fill_cube(const GCTAObservation* obs);
+    void fill_cube(const GCTAObservation* obs, GModels& models);
     bool has_cube(void) const;
-    GModels trim_models(const GModels& all_models, const GCTARoi& roi);
+    GModels trim_models(GModels& all_models, const GCTARoi& roi);
  
     // User parameters
     GFilename m_outcube;      //!< Output model cube