From 42c737c13279f5fa71d51575b2a6ab58a8e98bb9 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 24 Oct 2014 16:27:14 +0200 Subject: [PATCH] Reduce memory requirements for ctobssim. --- ChangeLog | 3 ++- NEWS | 8 +++++++- src/ctobssim/ctobssim.cpp | 432 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------------------------------------------------------------------------------------------------------------------------------------- src/ctobssim/ctobssim.hpp | 27 ++++++++------------------- 4 files changed, 284 insertions(+), 186 deletions(-) diff --git a/ChangeLog b/ChangeLog index 2a23821..53cb6d5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-10-17 Juergen Knoedlseder +2014-10-24 Juergen Knoedlseder * Version 0.8.0 released ======================== @@ -17,6 +17,7 @@ Support energy dispersion in ctobssim, ctlike and ctmodel Correct ctobssim for thread safety Correctly apply time references in ctselect + Reduce memory usage Move obsutils to ctools namespace Add several functions to obsutils Enhance documentation diff --git a/NEWS b/NEWS index dac6f72..726c6ec 100644 --- a/NEWS +++ b/NEWS @@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00 Author(s) : Juergen Knoedlseder -17 October 2014 +24 October 2014 1. Introduction @@ -48,6 +48,12 @@ avoiding concurrent memory access (#1258). ctobssim now also supports the handling of energy dispersion. Energy dispersion can be switched on or off using the hidden edisp parameter. +The memory needs for ctobssim have been reduced by first working on a +local copy of the CTA observation (this de-allocates properly and memory +cache associated to the response functions), and by saving and disposing +the event lists immediately if called using the execute() method (this +is the method called when the tool is run from the command line). + 6. ctselect - CTA event selection --------------------------------- diff --git a/src/ctobssim/ctobssim.cpp b/src/ctobssim/ctobssim.cpp index d70e096..3507537 100644 --- a/src/ctobssim/ctobssim.cpp +++ b/src/ctobssim/ctobssim.cpp @@ -92,7 +92,7 @@ ctobssim::ctobssim(const GObservations& obs) : * * Constructs application using command line arguments for parameter setting. ***************************************************************************/ -ctobssim::ctobssim(int argc, char *argv[]) : +ctobssim::ctobssim(int argc, char *argv[]) : ctool(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv) { // Initialise members @@ -218,18 +218,16 @@ void ctobssim::run(void) log << std::endl; } - // Write seed values into logger - if (logTerse()) { - log << std::endl; - log.header1("Seed values"); - for (int i = 0; i < m_rans.size(); ++i) { - log << gammalib::parformat("Seed "+gammalib::str(i)); - log << gammalib::str(m_rans[i].seed()) << std::endl; - } + // Special mode: if read ahead is specified we know that we called + // the execute() method, hence files are saved immediately and event + // lists are disposed afterwards. + if (read_ahead()) { + m_save_and_dispose = true; } - // Set energy dispersion flag for all CTA observations and save old - // values in save_edisp vector + // Determine the number of valid CTA observations, set energy dispersion flag + // for all CTA observations and save old values in save_edisp vector + int n_observations = 0; std::vector save_edisp; save_edisp.assign(m_obs.size(), false); for (int i = 0; i < m_obs.size(); ++i) { @@ -237,6 +235,43 @@ void ctobssim::run(void) if (obs != NULL) { save_edisp[i] = obs->response()->apply_edisp(); obs->response()->apply_edisp(m_apply_edisp); + n_observations++; + } + } + + // If more than a single observation has been handled then make sure that + // an XML file will be used for storage + if (n_observations > 1) { + m_use_xml = true; + } + + // Write execution mode into logger + if (logTerse()) { + log << std::endl; + log.header1("Execution mode"); + log << gammalib::parformat("Event list management"); + if (m_save_and_dispose) { + log << "Save and dispose (reduces memory needs)" << std::endl; + } + else { + log << "Keep events in memory" << std::endl; + } + log << gammalib::parformat("Output format"); + if (m_use_xml) { + log << "Write Observation Definition XML file" << std::endl; + } + else { + log << "Write single event list FITS file" << std::endl; + } + } + + // Write seed values into logger + if (logTerse()) { + log << std::endl; + log.header1("Seed values"); + for (int i = 0; i < m_rans.size(); ++i) { + log << gammalib::parformat("Seed "+gammalib::str(i)); + log << gammalib::str(m_rans[i].seed()) << std::endl; } } @@ -263,9 +298,6 @@ void ctobssim::run(void) } } - // Initialise counters - int n_observations = 0; - // From here on the code can be parallelized if OpenMP support // is enabled. The code in the following block corresponds to the // code that will be executed in each thread @@ -288,11 +320,11 @@ void ctobssim::run(void) wrklog.max_size(10000000); // Loop over all observation in the container. If OpenMP support - // is enabled, this looped will be parallelized. + // is enabled, this loop will be parallelized. #pragma omp for for (int i = 0; i < m_obs.size(); ++i) { - // Get CTA observation + // Get pointer on CTA observation GCTAObservation* obs = dynamic_cast(m_obs[i]); // Continue only if observation is a CTA observation @@ -308,26 +340,55 @@ void ctobssim::run(void) } } - // Increment counter - n_observations++; + // Work on a clone of the CTA observation. This makes sure that + // any memory allocated for computing (for example a response + // cache) is properly de-allocated on exit of this run + GCTAObservation obs_clone = *obs; // Save number of events before entering simulation - int events_before = obs->events()->size(); + int events_before = obs_clone.events()->size(); // Simulate source events - simulate_source(obs, models, m_rans[i], &wrklog); + simulate_source(&obs_clone, models, m_rans[i], &wrklog); // Simulate source events - simulate_background(obs, models, m_rans[i], &wrklog); + simulate_background(&obs_clone, models, m_rans[i], &wrklog); // Dump simulation results if (logNormal()) { wrklog << gammalib::parformat("MC events"); - wrklog << obs->events()->size() - events_before; + wrklog << obs_clone.events()->size() - events_before; wrklog << " (all models)"; wrklog << std::endl; } + // If requested, event lists are saved immediately + if (m_save_and_dispose) { + + // Set event output file name. If multiple observations are + // handled, build the filename from prefix and observation + // index. Otherwise use the outfile parameter. + std::string outfile; + if (m_use_xml) { + m_prefix = (*this)["prefix"].string(); + outfile = m_prefix + gammalib::str(i) + ".fits"; + } + else { + outfile = (*this)["outfile"].filename(); + } + + // Store output file name in original observation + obs->eventfile(outfile); + + // Save observation into FITS file + obs_clone.save(outfile, clobber()); + + } + + // ... otherwise append the event list to the original observation + else { + obs->events(*(obs_clone.events())); + } } // endif: CTA observation found @@ -341,12 +402,6 @@ void ctobssim::run(void) } } // end pragma omp parallel - - // If more than a single observation has been handled then make sure that - // an XML file will be used for storage - if (n_observations > 1) { - m_use_xml = true; - } // Restore energy dispersion flag for all CTA observations for (int i = 0; i < m_obs.size(); ++i) { @@ -409,6 +464,142 @@ void ctobssim::save(void) /***********************************************************************//** + * @brief Return observation container + * + * @return Reference to observation container + ***************************************************************************/ +const GObservations& ctobssim::obs(void) const +{ + // If event lists have been disposed, load them now into the observation + // container + if (m_save_and_dispose) { + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); + + // Handle only CTA observations + if (obs != NULL) { + + // Load observation from FITS file + obs->load(obs->eventfile()); + + } // endif: observation was a CTA observations + + } // endfor: looped over CTA observations + + // Now we have the event lists hence we can reset the disposal flag + m_save_and_dispose = false; + + } // endif: event lists have been disposed + + // Return observation container + return m_obs; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + * + * The thrown area is fixed to pi*(2500^2) m2, which is the same value that + * is used in the Monte Carlo simulations (information from Konrad Bernloehr) + ***************************************************************************/ +void ctobssim::init_members(void) +{ + // Initialise user parameters + m_infile.clear(); + m_outfile.clear(); + m_prefix.clear(); + m_seed = 1; + m_ra = 0.0; + m_dec = 0.0; + m_rad = 0.0; + m_tmin = 0.0; + m_tmax = 0.0; + m_emin = 0.0; + m_emax = 0.0; + m_apply_edisp = false; + m_deadc = 1.0; + + // Initialise protected members + m_rans.clear(); + m_obs.clear(); + m_use_xml = false; + m_save_and_dispose = false; + + // Set fixed parameters + m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2 + m_max_photons = 1000000; //!< Maximum number of photons / time slice + + // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD, + // which is defined in GCTALib.hpp. This is somehow a kluge. We need + // a better mechanism to implement the CTA reference MJD. + m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL"); + + // Initialise first event identifier + m_event_id = 1; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctobssim::copy_members(const ctobssim& app) +{ + // Copy user parameters + m_infile = app.m_infile; + m_outfile = app.m_outfile; + m_prefix = app.m_prefix; + m_seed = app.m_seed; + m_ra = app.m_ra; + m_dec = app.m_dec; + m_rad = app.m_rad; + m_tmin = app.m_tmin; + m_tmax = app.m_tmax; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_apply_edisp = app.m_apply_edisp; + m_deadc = app.m_deadc; + + // Copy protected members + m_area = app.m_area; + m_max_photons = app.m_max_photons; + m_rans = app.m_rans; + m_obs = app.m_obs; + m_use_xml = app.m_use_xml; + m_save_and_dispose = app.m_save_and_dispose; + m_cta_ref = app.m_cta_ref; + m_event_id = app.m_event_id; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctobssim::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** * @brief Get application parameters * * Get all task parameters from parameter file or (if required) by querying @@ -473,6 +664,16 @@ void ctobssim::get_parameters(void) m_prefix = (*this)["prefix"].string(); } + // Make sure that all observations hold a CTA event list + for (int i = 0; i < m_obs.size(); ++i) { + GCTAObservation* obs = dynamic_cast(m_obs[i]); + if (obs != NULL) { + if (dynamic_cast(obs->events()) == NULL) { + set_list(obs); + } + } + } + // Initialise random number generators. We initialise here one random // number generator per observation so that each observation will // get it's own random number generator. This will lead to identical @@ -624,21 +825,17 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, } // Get CTA response - const GCTAResponseIrf* rsp = dynamic_cast(obs->response()); + const GCTAResponseIrf* rsp = + dynamic_cast(obs->response()); if (rsp == NULL) { std::string msg = "Response is not an IRF response.\n" + obs->response()->print(); throw GException::invalid_value(G_SIMULATE_SOURCE, msg); } - // Make sure that the observation holds a CTA event list. If this - // is not the case then allocate and attach a CTA event list now. - if (dynamic_cast(obs->events()) == NULL) { - set_list(obs); - } - // Get pointer on event list (circumvent const correctness) - GCTAEventList* events = static_cast(const_cast(obs->events())); + GCTAEventList* events = + static_cast(const_cast(obs->events())); // Extract simulation region. GSkyDir dir = events->roi().centre().dir(); @@ -898,7 +1095,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels& models, ***************************************************************************/ void ctobssim::simulate_background(GCTAObservation* obs, const GModels& models, - GRan& ran, + GRan& ran, GLog* wrklog) { // Continue only if observation pointer is valid @@ -909,12 +1106,6 @@ void ctobssim::simulate_background(GCTAObservation* obs, wrklog = &log; } - // Make sure that the observation holds a CTA event list. If this - // is not the case then allocate and attach a CTA event list now. - if (dynamic_cast(obs->events()) == NULL) { - set_list(obs); - } - // Get pointer on event list (circumvent const correctness) GCTAEventList* events = static_cast(const_cast(obs->events())); @@ -934,7 +1125,7 @@ void ctobssim::simulate_background(GCTAObservation* obs, // Get simulated CTA event list. Note that this method // includes the deadtime correction. GCTAEventList* list = - dynamic_cast(model->mc(*obs, ran)); + dynamic_cast(model->mc(*obs, ran)); // Continue only if we got a CTA event list if (list != NULL) { @@ -984,103 +1175,6 @@ void ctobssim::simulate_background(GCTAObservation* obs, } -/*========================================================================== - = = - = Private methods = - = = - ==========================================================================*/ - -/***********************************************************************//** - * @brief Initialise class members - * - * The thrown area is fixed to pi*(2500^2) m2, which is the same value that - * is used in the Monte Carlo simulations (information from Konrad Bernloehr) - ***************************************************************************/ -void ctobssim::init_members(void) -{ - // Initialise user parameters - m_infile.clear(); - m_outfile.clear(); - m_prefix.clear(); - m_seed = 1; - m_ra = 0.0; - m_dec = 0.0; - m_rad = 0.0; - m_tmin = 0.0; - m_tmax = 0.0; - m_emin = 0.0; - m_emax = 0.0; - m_apply_edisp = false; - m_deadc = 1.0; - - // Initialise protected members - m_rans.clear(); - m_obs.clear(); - m_use_xml = false; - - // Set fixed parameters - m_area = 19634954.0 * 1.0e4; //!< pi*(2500^2) m^2 - m_max_photons = 1000000; //!< Maximum number of photons / time slice - - // Set CTA time reference. G_CTA_MJDREF is the CTA reference MJD, - // which is defined in GCTALib.hpp. This is somehow a kluge. We need - // a better mechanism to implement the CTA reference MJD. - m_cta_ref.set(G_CTA_MJDREF, "s", "TT", "LOCAL"); - - // Initialise first event identifier - m_event_id = 1; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Copy class members - * - * @param[in] app Application. - ***************************************************************************/ -void ctobssim::copy_members(const ctobssim& app) -{ - // Copy user parameters - m_infile = app.m_infile; - m_outfile = app.m_outfile; - m_prefix = app.m_prefix; - m_seed = app.m_seed; - m_ra = app.m_ra; - m_dec = app.m_dec; - m_rad = app.m_rad; - m_tmin = app.m_tmin; - m_tmax = app.m_tmax; - m_emin = app.m_emin; - m_emax = app.m_emax; - m_apply_edisp = app.m_apply_edisp; - m_deadc = app.m_deadc; - - // Copy protected members - m_area = app.m_area; - m_max_photons = app.m_max_photons; - m_rans = app.m_rans; - m_obs = app.m_obs; - m_use_xml = app.m_use_xml; - m_cta_ref = app.m_cta_ref; - m_event_id = app.m_event_id; - - // Return - return; -} - - -/***********************************************************************//** - * @brief Delete class members - ***************************************************************************/ -void ctobssim::free_members(void) -{ - // Return - return; -} - - /***********************************************************************//** * @brief Save event list in FITS format. * @@ -1089,14 +1183,19 @@ void ctobssim::free_members(void) ***************************************************************************/ void ctobssim::save_fits(void) { - // Get output filename - m_outfile = (*this)["outfile"].filename(); + // Save only if event list has not yet been saved and disposed + if (!m_save_and_dispose) { + + // Get output filename + m_outfile = (*this)["outfile"].filename(); - // Get CTA observation from observation container - GCTAObservation* obs = dynamic_cast(m_obs[0]); + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(m_obs[0]); - // Save observation into FITS file - obs->save(m_outfile, clobber()); + // Save observation into FITS file + obs->save(m_outfile, clobber()); + + } // endif: event list has not yet been saved and disposed // Return return; @@ -1131,34 +1230,37 @@ void ctobssim::save_xml(void) " definition files." << std::endl; } - // Initialise file number - int file_num = 0; + // Save only if event lists have not yet been saved and disposed + if (!m_save_and_dispose) { - // Loop over all observation in the container - for (int i = 0; i < m_obs.size(); ++i) { + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { - // Get CTA observation - GCTAObservation* obs = dynamic_cast(m_obs[i]); + // Get CTA observation + GCTAObservation* obs = dynamic_cast(m_obs[i]); - // Handle only CTA observations - if (obs != NULL) { + // Handle only CTA observations + if (obs != NULL) { - // Set event output file name - std::string outfile = m_prefix + gammalib::str(file_num) + ".fits"; + // Continue only if there is an event list (it may have been disposed) + if (obs->events()->size() != 0) { - // Store output file name in observation - obs->eventfile(outfile); + // Set event output file name + std::string outfile = m_prefix + gammalib::str(i) + ".fits"; - // Save observation into FITS file - obs->save(outfile, clobber()); + // Store output file name in observation + obs->eventfile(outfile); - // Increment file number - file_num++; + // Save observation into FITS file + obs->save(outfile, clobber()); + } - } // endif: observation was a CTA observations + } // endif: observation was a CTA observations - } // endfor: looped over observations + } // endfor: looped over observations + + } // endif: event list has not yet been saved and disposed // Save observations in XML file m_obs.save(m_outfile); diff --git a/src/ctobssim/ctobssim.hpp b/src/ctobssim/ctobssim.hpp index 07a9de4..3e8636d 100644 --- a/src/ctobssim/ctobssim.hpp +++ b/src/ctobssim/ctobssim.hpp @@ -101,25 +101,14 @@ protected: double m_deadc; //!< Average deadtime correction // Protected members - double m_area; //!< Surface area for simulation (cm2) - int m_max_photons; //!< Maximum number of photons/slice - std::vector m_rans; //!< Random number generators - GObservations m_obs; //!< Observation container - bool m_use_xml; //!< Use XML file instead of FITS file - GTimeReference m_cta_ref; //!< CTA time reference - int m_event_id; //!< Event identifier + mutable GObservations m_obs; //!< Observation container + mutable bool m_save_and_dispose; //!< Save and dispose immediately + double m_area; //!< Surface area for simulation (cm2) + int m_max_photons; //!< Maximum number of photons/slice + std::vector m_rans; //!< Random number generators + bool m_use_xml; //!< Use XML file instead of FITS file + GTimeReference m_cta_ref; //!< CTA time reference + int m_event_id; //!< Event identifier }; - -/***********************************************************************//** - * @brief Return observation container - * - * @return Reference to observation container - ***************************************************************************/ -inline -const GObservations& ctobssim::obs(void) const -{ - return m_obs; -} - #endif /* CTOBSSIM_HPP */ -- libgit2 0.21.4