Pierrick Martin / gammalib

Commit 82c53c6fa15eba93ee1967bab649afa6d901f71e

Authored by Jürgen Knödlseder
1 parent ca8209d4

Add and use GModel::has_eval_indices() and GModel::eval_indices() methods (#745) 

The methods GModel::has_eval_indices() and GModel::eval_indices() allow for a more efficient parameter handling in GObservation::model() that is specifically relevant for a model with many parameters, such as the SPI data space models. These methods allow the GModel::eval() method to signal which of the parameters were actually updated, so that the gradients for only these parameters will be considered in GObservation::model(), while all other gradients will be ignored.
Inline Side-by-side
Showing 4 changed files with 170 additions and 91 deletions
include/GModel.hpp
  View file @82c53c6
1 1 /***************************************************************************
2 2 * GModel.hpp - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -126,29 +126,31 @@ public:
126 126 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
127 127  
128 128 // Implemented methods
129   - int size(void) const;
130   - int scales(void) const;
131   - GModelPar& at(const int& index);
132   - const GModelPar& at(const int& index) const;
133   - bool has_par(const std::string& name) const;
134   - bool has_scales(void) const;
135   - const std::string& name(void) const;
136   - void name(const std::string& name);
137   - const double& ts(void) const;
138   - void ts(const double& ts);
139   - const bool& tscalc(void) const;
140   - void tscalc(const bool& tscalc);
141   - const bool& has_ts(void) const;
142   - std::string instruments(void) const;
143   - void instruments(const std::string& instruments);
144   - GModelPar& scale(const int& index);
145   - const GModelPar& scale(const int& index) const;
146   - GModelPar scale(const std::string& instrument) const;
147   - void scale(const GModelPar& par);
148   - std::string ids(void) const;
149   - void ids(const std::string& ids);
150   - bool is_valid(const std::string& instruments,
151   - const std::string& ids) const;
  129 + int size(void) const;
  130 + int scales(void) const;
  131 + GModelPar& at(const int& index);
  132 + const GModelPar& at(const int& index) const;
  133 + bool has_par(const std::string& name) const;
  134 + bool has_scales(void) const;
  135 + const std::string& name(void) const;
  136 + void name(const std::string& name);
  137 + const double& ts(void) const;
  138 + void ts(const double& ts);
  139 + const bool& tscalc(void) const;
  140 + void tscalc(const bool& tscalc);
  141 + const bool& has_ts(void) const;
  142 + std::string instruments(void) const;
  143 + void instruments(const std::string& instruments);
  144 + GModelPar& scale(const int& index);
  145 + const GModelPar& scale(const int& index) const;
  146 + GModelPar scale(const std::string& instrument) const;
  147 + void scale(const GModelPar& par);
  148 + std::string ids(void) const;
  149 + void ids(const std::string& ids);
  150 + bool is_valid(const std::string& instruments,
  151 + const std::string& ids) const;
  152 + const bool& has_eval_indices(void) const;
  153 + const std::vector<int>& eval_indices(void) const;
152 154  
153 155 protected:
154 156 // Protected methods
... ... @@ -171,6 +173,12 @@ protected:
171 173 bool m_has_tscalc; //!< Signals if tscalc attribute is available
172 174 bool m_tscalc; //!< Signals if TS should be computed
173 175 double m_ts; //!< Test Statistic of the model
  176 +
  177 + // Protected members for efficient gradient access
  178 + mutable bool m_has_eval_inx; //!< Signal that parameter indices
  179 + //!< updated by eval() method exist
  180 + mutable std::vector<int> m_eval_inx; //!< List of parameter indices updated
  181 + //!< by eval() method
174 182 };
175 183  
176 184  
... ... @@ -352,4 +360,28 @@ const bool&amp; GModel::has_ts(void) const
352 360 return m_has_ts;
353 361 }
354 362  
  363 +
  364 +/***********************************************************************//**
  365 + * @brief Signals that parameter indices updated by eval() method were set
  366 + *
  367 + * @return True if parameter indices updated by eval() method were set.
  368 + ***************************************************************************/
  369 +inline
  370 +const bool& GModel::has_eval_indices(void) const
  371 +{
  372 + return m_has_eval_inx;
  373 +}
  374 +
  375 +
  376 +/***********************************************************************//**
  377 + * @brief Return vector of parameter indices updated by eval() method
  378 + *
  379 + * @return Vector of parameter indices updated by eval() method.
  380 + ***************************************************************************/
  381 +inline
  382 +const std::vector<int>& GModel::eval_indices(void) const
  383 +{
  384 + return m_eval_inx;
  385 +}
  386 +
355 387 #endif /* GMODEL_HPP */
... ...
pyext/GModel.i
  View file @82c53c6
1 1 /***************************************************************************
2 2 * GModel.i - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -59,26 +59,28 @@ public:
59 59 virtual void write(GXmlElement& xml) const = 0;
60 60  
61 61 // Implemented methods
62   - int size(void) const;
63   - int scales(void) const;
64   - bool has_par(const std::string& name) const;
65   - bool has_scales(void) const;
66   - const std::string& name(void) const;
67   - void name(const std::string& name);
68   - const double& ts(void) const;
69   - void ts(const double& ts);
70   - const bool& tscalc(void) const;
71   - void tscalc(const bool& tscalc);
72   - const bool& has_ts(void) const;
73   - std::string instruments(void) const;
74   - void instruments(const std::string& instruments);
75   - GModelPar& scale(const int& index);
76   - GModelPar scale(const std::string& instrument) const;
77   - void scale(const GModelPar& par);
78   - std::string ids(void) const;
79   - void ids(const std::string& ids);
80   - bool is_valid(const std::string& instruments,
81   - const std::string& ids) const;
  62 + int size(void) const;
  63 + int scales(void) const;
  64 + bool has_par(const std::string& name) const;
  65 + bool has_scales(void) const;
  66 + const std::string& name(void) const;
  67 + void name(const std::string& name);
  68 + const double& ts(void) const;
  69 + void ts(const double& ts);
  70 + const bool& tscalc(void) const;
  71 + void tscalc(const bool& tscalc);
  72 + const bool& has_ts(void) const;
  73 + std::string instruments(void) const;
  74 + void instruments(const std::string& instruments);
  75 + GModelPar& scale(const int& index);
  76 + GModelPar scale(const std::string& instrument) const;
  77 + void scale(const GModelPar& par);
  78 + std::string ids(void) const;
  79 + void ids(const std::string& ids);
  80 + bool is_valid(const std::string& instruments,
  81 + const std::string& ids) const;
  82 + const bool& has_eval_indices(void) const;
  83 + const std::vector<int>& eval_indices(void) const;
82 84 };
83 85  
84 86  
... ...
src/model/GModel.cpp
  View file @82c53c6
1 1 /***************************************************************************
2 2 * GModel.cpp - Abstract virtual model base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2009-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2009-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -617,11 +617,13 @@ void GModel::init_members(void)
617 617 m_scales.clear();
618 618 m_ids.clear();
619 619 m_pars.clear();
620   - m_ts = 0.0;
621   - m_has_ts = false;
622   - m_has_tscalc = false;
623   - m_tscalc = false;
624   -
  620 + m_ts = 0.0;
  621 + m_has_ts = false;
  622 + m_has_tscalc = false;
  623 + m_tscalc = false;
  624 + m_has_eval_inx = false;
  625 + m_eval_inx.clear();
  626 +
625 627 // Return
626 628 return;
627 629 }
... ... @@ -635,15 +637,17 @@ void GModel::init_members(void)
635 637 void GModel::copy_members(const GModel& model)
636 638 {
637 639 // Copy members
638   - m_name = model.m_name;
639   - m_instruments = model.m_instruments;
640   - m_scales = model.m_scales;
641   - m_ids = model.m_ids;
642   - m_pars = model.m_pars;
643   - m_ts = model.m_ts;
644   - m_has_ts = model.m_has_ts;
645   - m_has_tscalc = model.m_has_tscalc;
646   - m_tscalc = model.m_tscalc;
  640 + m_name = model.m_name;
  641 + m_instruments = model.m_instruments;
  642 + m_scales = model.m_scales;
  643 + m_ids = model.m_ids;
  644 + m_pars = model.m_pars;
  645 + m_ts = model.m_ts;
  646 + m_has_ts = model.m_has_ts;
  647 + m_has_tscalc = model.m_has_tscalc;
  648 + m_tscalc = model.m_tscalc;
  649 + m_has_eval_inx = model.m_has_eval_inx;
  650 + m_eval_inx = model.m_eval_inx;
647 651  
648 652 // Return
649 653 return;
... ...
src/obs/GObservation.cpp
  View file @82c53c6
1 1 /***************************************************************************
2 2 * GObservation.cpp - Abstract observation base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -302,39 +302,80 @@ double GObservation::model(const GModels&amp; models,
302 302 // over the energy dispersion. Here it's simpler to just use
303 303 // numerical gradients.
304 304 if (gradient != NULL) {
305   - for (int ipar = 0; ipar < mptr->size(); ++ipar) {
306   -
307   - // Get reference to model parameter
308   - const GModelPar& par = (*mptr)[ipar];
309   -
310   - // Make sure that we have a slot for the gradient
311   - #if defined(G_RANGE_CHECK)
312   - if (igrad >= grad_size) {
313   - std::string msg = "Vector has not enough elements "
314   - "to store the model parameter "
315   - "gradients. "+
316   - gammalib::str(models.npars())+
317   - " elements requested while vector "
318   - "only contains "+
319   - gammalib::str(gradient->size())+
320   - " elements.";
321   - throw GException::invalid_value(G_MODEL, msg);
322   - }
323   - #endif
324   -
325   - if (par.is_free()) {
326   - if (par.has_grad() && !response()->use_edisp()) {
327   - (*gradient)[igrad+ipar] = par.factor_gradient();
  305 +
  306 + // Make sure that we have a slot for the gradient
  307 + #if defined(G_RANGE_CHECK)
  308 + if (igrad+mptr->size() > grad_size) {
  309 + std::string msg = "Vector has not enough elements to "
  310 + "store the model parameter gradients. "+
  311 + gammalib::str(models.npars())+" elements "
  312 + "requested while vector only contains "+
  313 + gammalib::str(gradient->size())+" elements.";
  314 + throw GException::invalid_value(G_MODEL, msg);
  315 + }
  316 + #endif
  317 +
  318 + // If model provides the parameter indices that were
  319 + // updated by the eval() method then use them as this is
  320 + // less time consuming than looping over all indices
  321 + if (mptr->has_eval_indices()) {
  322 +
  323 + // Get number of relevant parameters
  324 + int npars = mptr->eval_indices().size();
  325 +
  326 + // Loop over all relevant parameters
  327 + for (int index = 0; index < npars; ++index) {
  328 +
  329 + // Retrieve parameter index
  330 + int ipar = mptr->eval_indices()[index];
  331 +
  332 + // Get reference to model parameter
  333 + const GModelPar& par = (*mptr)[ipar];
  334 +
  335 + // Set gradient
  336 + if (par.is_free()) {
  337 + if (par.has_grad() && !response()->use_edisp()) {
  338 + (*gradient)[igrad+ipar] = par.factor_gradient();
  339 + }
  340 + else {
  341 + (*gradient)[igrad+ipar] = model_grad(*mptr, par, event);
  342 + }
328 343 }
329 344 else {
330   - (*gradient)[igrad+ipar] = model_grad(*mptr, par, event);
  345 + (*gradient)[igrad+ipar] = 0.0;
331 346 }
332   - }
333   - else {
334   - (*gradient)[igrad+ipar] = 0.0;
335   - }
336   - }
337   - }
  347 +
  348 + } // endfor: looped over all parameter indices
  349 +
  350 + } // endif: parameter indices were available
  351 +
  352 + // ... otherwise loop over all parameter indices
  353 + else {
  354 +
  355 + // Loop over all parameters
  356 + for (int ipar = 0; ipar < mptr->size(); ++ipar) {
  357 +
  358 + // Get reference to model parameter
  359 + const GModelPar& par = (*mptr)[ipar];
  360 +
  361 + // Set gradient
  362 + if (par.is_free()) {
  363 + if (par.has_grad() && !response()->use_edisp()) {
  364 + (*gradient)[igrad+ipar] = par.factor_gradient();
  365 + }
  366 + else {
  367 + (*gradient)[igrad+ipar] = model_grad(*mptr, par, event);
  368 + }
  369 + }
  370 + else {
  371 + (*gradient)[igrad+ipar] = 0.0;
  372 + }
  373 +
  374 + } // endfor: looped over all parameter indices
  375 +
  376 + } // endelse: no parameter indices were available
  377 +
  378 + } // endif: use model gradients
338 379  
339 380 } // endif: model component was valid for instrument
340 381  
... ...