Johan Bregeon / ctools

Commit 23dd86c9dace430ff5ac4bc2253e67ebe7ec4657

Authored by Jürgen Knödlseder
2 parents 1e0ec9e5 2c3a446b

Merge branch '1479-ctulimit-should-reoptimize-other-parameters' into integration

Inline Side-by-side
Showing 6 changed files with 83 additions and 35 deletions
pyext/ctulimit.i
  View file @23dd86c
... ... @@ -48,6 +48,7 @@ public:
48 48 void run(void);
49 49 void save(void);
50 50 const GObservations& obs(void) const;
  51 + const GOptimizer* opt(void) const;
51 52 const double& diff_ulimit(void) const;
52 53 const double& flux_ulimit(void) const;
53 54 const double& eflux_ulimit(void) const;
... ...
scripts/csspec.py
  View file @23dd86c
... ... @@ -21,7 +21,7 @@
21 21 import gammalib
22 22 import ctools
23 23 import sys
24   -import tempfile
  24 +
25 25  
26 26 # ============ #
27 27 # csspec class #
... ...
src/cterror/cterror.cpp
  View file @23dd86c
... ... @@ -232,20 +232,17 @@ void cterror::run(void)
232 232 log.header1("Compute best-fit likelihood");
233 233 }
234 234  
235   - // Allocate optimizer
236   - m_opt = new GOptimizerLM();
237   -
238 235 // Optimize and save best log-likelihood
239   - m_obs.optimize(*m_opt);
240   - m_obs.errors(*m_opt);
  236 + m_obs.optimize(m_opt);
  237 + m_obs.errors(m_opt);
241 238 m_best_logL = m_obs.logL();
242 239  
243 240 // Store optimizer for later recovery
244   - GOptimizerLM best_opt = *m_opt;
  241 + GOptimizerLM best_opt = m_opt;
245 242  
246 243 // Write optimised model into logger
247 244 if (logTerse()) {
248   - log << *m_opt << std::endl;
  245 + log << m_opt << std::endl;
249 246 log << gammalib::parformat("Maximum log likelihood");
250 247 log << gammalib::str(m_best_logL,3) << std::endl;
251 248 log << m_obs.models() << std::endl;
... ... @@ -356,7 +353,7 @@ void cterror::run(void)
356 353 } // endif: source model exists
357 354  
358 355 // Recover optimizer
359   - *m_opt = best_opt;
  356 + m_opt = best_opt;
360 357  
361 358 // Return
362 359 return;
... ... @@ -409,10 +406,15 @@ void cterror::init_members(void)
409 406  
410 407 // Initialise protected members
411 408 m_obs.clear();
  409 + m_opt.clear();
412 410 m_dlogL = 0.0;
413 411 m_best_logL = 0.0;
414 412 m_model_par = NULL;
415 413  
  414 + // Set optimizer parameters
  415 + m_opt.max_iter(m_max_iter);
  416 + m_opt.max_stalls(10);
  417 +
416 418 // Return
417 419 return;
418 420 }
... ... @@ -435,6 +437,7 @@ void cterror::copy_members(const cterror&amp; app)
435 437 m_obs = app.m_obs;
436 438 m_dlogL = app.m_dlogL;
437 439 m_best_logL = app.m_best_logL;
  440 + m_opt = app.m_opt;
438 441 m_model_par = NULL;
439 442  
440 443 // Return
... ... @@ -657,7 +660,7 @@ double cterror::evaluate(const double&amp; value)
657 660 m_model_par->fix();
658 661  
659 662 // Re-optimize
660   - m_obs.optimize(*m_opt);
  663 + m_obs.optimize(m_opt);
661 664  
662 665 // Retrieve likelihood
663 666 logL = m_obs.logL();
... ...
src/cterror/cterror.hpp
  View file @23dd86c
... ... @@ -84,7 +84,7 @@ protected:
84 84 double m_dlogL; //!< Likelihood difference for upper limit computation
85 85 GModelPar* m_model_par; //!< Pointer to model parameter
86 86 double m_best_logL; //!< Best fit log likelihood of given model
87   - GOptimizerLM* m_opt; //!< Optimizer
  87 + GOptimizerLM m_opt; //!< Optimizer
88 88  
89 89 };
90 90  
... ... @@ -92,7 +92,7 @@ protected:
92 92 /***********************************************************************//**
93 93 * @brief Return observation container
94 94 *
95   - * @return Reference to observation container
  95 + * @return Reference to observation container.
96 96 ***************************************************************************/
97 97 inline
98 98 const GObservations& cterror::obs(void) const
... ... @@ -100,15 +100,16 @@ const GObservations&amp; cterror::obs(void) const
100 100 return m_obs;
101 101 }
102 102  
  103 +
103 104 /***********************************************************************//**
104 105 * @brief Return optimizer
105 106 *
106   - * @return Pointer to optimizer
  107 + * @return Pointer to optimizer.
107 108 ***************************************************************************/
108 109 inline
109 110 const GOptimizer* cterror::opt(void) const
110 111 {
111   - return m_opt;
  112 + return &m_opt;
112 113 }
113 114  
114 115 #endif /* CTERROR_HPP */
... ...
src/ctulimit/ctulimit.cpp
  View file @23dd86c
... ... @@ -232,6 +232,9 @@ void ctulimit::run(void)
232 232 // Save original models
233 233 GModels models_orig = m_obs.models();
234 234  
  235 + // Initialise best fit optimizer
  236 + GOptimizerLM best_opt;
  237 +
235 238 // Save original log-likelihood. If the value is zero it has never been
236 239 // computed hence we compute it now.
237 240 m_best_logL = m_obs.logL();
... ... @@ -243,24 +246,32 @@ void ctulimit::run(void)
243 246 log.header1("Compute best-fit likelihood");
244 247 }
245 248  
246   - // Reoptimize if likelihood was not given before
247   - GOptimizerLM* opt = new GOptimizerLM();
248   - m_obs.optimize(*opt);
249   - m_obs.errors(*opt);
250   - m_best_logL = m_obs.logL();
  249 + // Optimize and save best log-likelihood
  250 + m_obs.optimize(m_opt);
  251 + m_obs.errors(m_opt);
  252 + m_best_logL = m_obs.logL();
251 253  
252   - // Write optimised model into logger
253   - if (logTerse()) {
254   - log << m_obs.models() << std::endl;
255   - }
  254 + // Store optimizer for later recovery
  255 + best_opt = m_opt;
  256 +
  257 + // Write optimised model into logger
  258 + if (logTerse()) {
  259 + log << m_opt << std::endl;
  260 + log << gammalib::parformat("Maximum log likelihood");
  261 + log << gammalib::str(m_best_logL,3) << std::endl;
  262 + log << m_obs.models() << std::endl;
  263 + }
256 264  
257 265 } // endif: likelihood was zero
258 266  
  267 + // Extract current value
  268 + double value = m_model_par->factor_value();
  269 +
259 270 // Compute parameter bracketing
260   - double value = m_model_par->value();
261   - double error = m_model_par->error();
262   - double parmin = value + m_sigma_min * error;
263   - double parmax = value + m_sigma_max * error;
  271 + double parmin = std::max(m_model_par->factor_min(),
  272 + value - m_sigma_min * m_model_par->factor_error());
  273 + double parmax = std::min(m_model_par->factor_max(),
  274 + value + m_sigma_max * m_model_par->factor_error());
264 275  
265 276 // Write header
266 277 if (logTerse()) {
... ... @@ -329,6 +340,9 @@ void ctulimit::run(void)
329 340 // Recover original models
330 341 m_obs.models(models_orig);
331 342  
  343 + // Recover optimizer
  344 + m_opt = best_opt;
  345 +
332 346 // Return
333 347 return;
334 348 }
... ... @@ -394,6 +408,7 @@ void ctulimit::init_members(void)
394 408  
395 409 // Initialise protected members
396 410 m_obs.clear();
  411 + m_opt.clear();
397 412 m_dlogL = 0.0;
398 413 m_skymodel = NULL;
399 414 m_model_par = NULL;
... ... @@ -402,6 +417,10 @@ void ctulimit::init_members(void)
402 417 m_diff_ulimit = 0.0;
403 418 m_eflux_ulimit = 0.0;
404 419  
  420 + // Set optimizer parameters
  421 + m_opt.max_iter(m_max_iter);
  422 + m_opt.max_stalls(10);
  423 +
405 424 // Return
406 425 return;
407 426 }
... ... @@ -433,6 +452,7 @@ void ctulimit::copy_members(const ctulimit&amp; app)
433 452 m_diff_ulimit = app.m_diff_ulimit;
434 453 m_flux_ulimit = app.m_flux_ulimit;
435 454 m_eflux_ulimit = app.m_eflux_ulimit;
  455 + m_opt = app.m_opt;
436 456  
437 457 // Extract model parameter
438 458 get_model_parameter();
... ... @@ -662,17 +682,23 @@ double ctulimit::evaluate(const double&amp; value)
662 682 double logL = 0.0;
663 683  
664 684 // Check if given parameter is within boundaries
665   - if (value > m_model_par->min() && value < m_model_par->max()) {
  685 + if (value > m_model_par->factor_min() && value < m_model_par->factor_max()) {
666 686  
667 687 // Change parameter factor
668   - m_model_par->value(value);
  688 + m_model_par->factor_value(value);
669 689  
670   - // Evaluate likelihood for new model container
671   - m_obs.eval();
  690 + // Fix parameter
  691 + m_model_par->fix();
  692 +
  693 + // Re-optimize
  694 + m_obs.optimize(m_opt);
672 695  
673 696 // Retrieve likelihood
674 697 logL = m_obs.logL();
675 698  
  699 + // Free parameter
  700 + m_model_par->free();
  701 +
676 702 } // endif: value was inside allowed range
677 703  
678 704 // ... otherwise signal that the parameter went outside the boundaries
... ...
src/ctulimit/ctulimit.hpp
  View file @23dd86c
... ... @@ -68,6 +68,7 @@ public:
68 68 void run(void);
69 69 void save(void);
70 70 const GObservations& obs(void) const;
  71 + const GOptimizer* opt(void) const;
71 72 const double& diff_ulimit(void) const;
72 73 const double& flux_ulimit(void) const;
73 74 const double& eflux_ulimit(void) const;
... ... @@ -102,13 +103,14 @@ protected:
102 103 double m_diff_ulimit; //!< Differential upper limit value
103 104 double m_flux_ulimit; //!< Flux upper limit value
104 105 double m_eflux_ulimit; //!< Energy flux upper limits
  106 + GOptimizerLM m_opt; //!< Optimizer
105 107 };
106 108  
107 109  
108 110 /***********************************************************************//**
109 111 * @brief Return observation container
110 112 *
111   - * @return Reference to observation container
  113 + * @return Reference to observation container.
112 114 ***************************************************************************/
113 115 inline
114 116 const GObservations& ctulimit::obs(void) const
... ... @@ -116,10 +118,23 @@ const GObservations&amp; ctulimit::obs(void) const
116 118 return m_obs;
117 119 }
118 120  
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Return optimizer
  124 + *
  125 + * @return Pointer to optimizer.
  126 + ***************************************************************************/
  127 +inline
  128 +const GOptimizer* ctulimit::opt(void) const
  129 +{
  130 + return &m_opt;
  131 +}
  132 +
  133 +
119 134 /***********************************************************************//**
120 135 * @brief Return differential upper limit
121 136 *
122   - * @return differential upper limit
  137 + * @return Differential upper flux limit.
123 138 ***************************************************************************/
124 139 inline
125 140 const double& ctulimit::diff_ulimit(void) const
... ... @@ -127,10 +142,11 @@ const double&amp; ctulimit::diff_ulimit(void) const
127 142 return m_diff_ulimit;
128 143 }
129 144  
  145 +
130 146 /***********************************************************************//**
131 147 * @brief Return flux upper limit
132 148 *
133   - * @return flux upper limit
  149 + * @return Upper flux limit.
134 150 ***************************************************************************/
135 151 inline
136 152 const double& ctulimit::flux_ulimit(void) const
... ... @@ -138,10 +154,11 @@ const double&amp; ctulimit::flux_ulimit(void) const
138 154 return m_flux_ulimit;
139 155 }
140 156  
  157 +
141 158 /***********************************************************************//**
142 159 * @brief return energy flux upper limit
143 160 *
144   - * @return energy flux upper limits
  161 + * @return Upper energy flux limit.
145 162 ***************************************************************************/
146 163 inline
147 164 const double& ctulimit::eflux_ulimit(void) const
... ...