From e607c99590a90554a6e1f7d3e9b0ab97975a8d49 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 20 Jul 2020 16:55:36 +0200 Subject: [PATCH] Speed-up stacked event binning, in particular when parallel processing is enabled. --- src/ctbin/ctbin.cpp | 99 ++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------- src/ctbin/ctbin.hpp | 10 +++++----- 2 files changed, 45 insertions(+), 64 deletions(-) diff --git a/src/ctbin/ctbin.cpp b/src/ctbin/ctbin.cpp index b703f04..a41108e 100755 --- a/src/ctbin/ctbin.cpp +++ b/src/ctbin/ctbin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2018 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -89,7 +89,6 @@ ctbin::ctbin(const GObservations& obs) : } - /***********************************************************************//** * @brief Command line constructor * @@ -221,9 +220,9 @@ void ctbin::clear(void) * @brief Run the event binning tool * * Gets the user parameters and loops over all CTA observations in the - * observation container to bin the events into a single counts cube. All - * observations in the observation container that do not contain CTA event - * lists will be skipped. + * observation container to bin the events into either a single or a set of + * counts cube. All observations in the observation container that do not + * contain CTA event lists will be skipped. ***************************************************************************/ void ctbin::run(void) { @@ -281,9 +280,15 @@ void ctbin::run(void) // Initialise arrays if (nobs > 0) { - m_counts = std::vector(nobs); - m_weights = std::vector(nobs); - m_cubes = std::vector(nobs); + + // Determine number of cubes + int ncubes = (m_stack) ? 1 : nobs; + + // Allocate cubes + m_counts = std::vector(ncubes); + m_weights = std::vector(ncubes); + m_cubes = std::vector(ncubes); + } else { m_counts.clear(); @@ -316,46 +321,25 @@ void ctbin::run(void) // Get pointer to observation GCTAObservation* obs = obs_list[i]; - // Fill the cube - GSkyMap counts = fill_cube(obs); + // Determine counts cube slot + int islot = (m_stack) ? 0 : i; + + // If slot is empty then allocate cubes + #pragma omp critical(ctbin_run) + if (m_counts[islot].is_empty()) { + m_counts[islot] = create_cube(obs); + m_weights[islot] = create_cube(obs); + } + + // Fill counts cube + fill_cube(obs, m_counts[islot]); // Set the counts cube weights - GSkyMap weights = set_weights(obs); + set_weights(obs, m_weights[islot]); // Dispose events to free memory obs->dispose_events(); - // Post process cubes - #pragma omp critical(ctbin_fill_cube) - { - // If stacking of cubes is requested then add up counts cube - // and set weights in the first slot of the m_counts and - // m_weights arrays - if (m_stack) { - if (m_counts[0].is_empty()) { - m_counts[0] = counts; - m_weights[0] = weights; - } - else { - for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) { - for (int pixel = 0; pixel < counts.npix(); ++pixel) { - m_counts[0](pixel, iebin) += counts(pixel, iebin); - if (weights(pixel, iebin) == 1.0) { - m_weights[0](pixel, iebin) = 1.0; - } - } - } - } - } - - // ... otherwise store counts and weights - else { - m_counts[i] = counts; - m_weights[i] = weights; - } - - } // end: omp critical section - } // endfor: looped over observations // Combine GTIs of all observations and compute total ontime and livetime @@ -369,16 +353,16 @@ void ctbin::run(void) // If resulting GTI is not filled then copy over the reference from // the observation and use it for the resulting GTI - if (m_gti.is_empty()) { + if (m_gti.is_empty()) { m_gti.reference(gti.reference()); - } + } - // Append GTIs - m_gti.extend(gti); + // Append GTIs + m_gti.extend(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(); } // endfor: looped over all observations @@ -842,13 +826,14 @@ GSkyMap ctbin::create_cube(const GCTAObservation* obs) * @brief Fill events into counts cube * * @param[in] obs CTA observation. + * @param[in,out] counts Counts cube * * @exception GException::invalid_value * No event list or valid RoI found in observation. * * Fills the events from an event list into the counts cube. ***************************************************************************/ -GSkyMap ctbin::fill_cube(const GCTAObservation* obs) +void ctbin::fill_cube(const GCTAObservation* obs, GSkyMap& counts) { // Make sure that the observation holds a CTA event list. If this is // not the case then throw an exception. @@ -883,9 +868,6 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) int num_outside_ebds = 0; int num_in_map = 0; - // Create counts cube - GSkyMap counts = create_cube(obs); - // Extract counts cube boundaries double pixel_x_max = counts.nx() - 0.5; double pixel_y_max = counts.ny() - 0.5; @@ -936,6 +918,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) } // Fill event in skymap + #pragma omp atomic counts(pixel, iebin) += 1.0; // Increment number of maps @@ -956,7 +939,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) } // Return - return counts; + return; } @@ -964,6 +947,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) * @brief Set counts cube weights for a given observation * * @param[in] obs CTA observation. + * @param[in,out] weights Counts cube weights. * * @exception GException::invalid_value * No event list or valid RoI found in observation. @@ -971,7 +955,7 @@ GSkyMap ctbin::fill_cube(const GCTAObservation* obs) * Sets the counts cube weights for all bins that are considered for the * specific observation to unity. ***************************************************************************/ -GSkyMap ctbin::set_weights(const GCTAObservation* obs) +void ctbin::set_weights(const GCTAObservation* obs, GSkyMap& weights) { // Get the RoI const GCTARoi& roi = obs->roi(); @@ -992,9 +976,6 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs) // Get counts cube usage flags std::vector usage = cube_layer_usage(m_ebounds, obs->ebounds()); - // Create weights cube - GSkyMap weights = create_cube(obs); - // Set sky direction cache availability flag bool cache = (m_stack && m_dirs.size() == weights.npix()); @@ -1029,7 +1010,7 @@ GSkyMap ctbin::set_weights(const GCTAObservation* obs) } // endfor: looped over pixels of counts cube // Return - return weights; + return; } diff --git a/src/ctbin/ctbin.hpp b/src/ctbin/ctbin.hpp index e60f01a..063c32e 100755 --- a/src/ctbin/ctbin.hpp +++ b/src/ctbin/ctbin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * ctbin - Event binning tool * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2018 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -44,8 +44,8 @@ * FITS files, and on observation definition XML files. * * If multiple event lists are specified in the observation container or the - * XML definition file, the class will merge these events into a single - * counts cube. + * XML definition file, the class will upon request stack these events into a + * single counts cube. * * If the hidden parameter stack=no is used, one counts cube is generated for * each individual observation. These cubes are saved to a path that can be @@ -94,8 +94,8 @@ protected: void get_parameters(void); void init_sky_dir_cache(void); GSkyMap create_cube(const GCTAObservation* obs); - GSkyMap fill_cube(const GCTAObservation* obs); - GSkyMap set_weights(const GCTAObservation* obs); + void fill_cube(const GCTAObservation* obs, GSkyMap& counts); + void set_weights(const GCTAObservation* obs, GSkyMap& weights); void obs_cube_stacked(void); // User parameters -- libgit2 0.21.4