Justine Devin / gammalib

Commit f949a672d322d867d589b8c9f3aeebdcfa9892e1

Authored by Jürgen Knödlseder
1 parent 969d7a4d

Adapt model and Npred gradient computation to cope correctly with parameter boun… 

…daries. Also streamlined the interface to avoid passing parameter indices. We now pass parameter references.
Inline Side-by-side
Showing 3 changed files with 134 additions and 182 deletions
include/GObservation.hpp
  View file @f949a67
... ... @@ -98,10 +98,9 @@ public:
98 98 GVector* gradient = NULL) const;
99 99 virtual double model_grad(const GModel& model,
100 100 const GModelPar& par,
101   - const int& ipar,
102 101 const GEvent& event) const;
103   - virtual double npred_grad(const GModel& model,
104   - const int& ipar) const;
  102 + virtual double npred_grad(const GModel& model,
  103 + const GModelPar& par) const;
105 104  
106 105 // Implemented methods
107 106 void name(const std::string& name);
... ... @@ -138,18 +137,18 @@ protected:
138 137 public:
139 138 model_func(const GObservation* parent,
140 139 const GModel& model,
141   - const GEvent& event,
142   - int ipar) :
  140 + GModelPar* par,
  141 + const GEvent& event) :
143 142 m_parent(parent),
144   - m_model(&model),
145   - m_event(&event),
146   - m_ipar(ipar) { }
  143 + m_model(model),
  144 + m_par(par),
  145 + m_event(event) { }
147 146 double eval(const double& x);
148 147 protected:
149 148 const GObservation* m_parent; //!< Pointer to parent
150   - const GModel* m_model; //!< Pointer to model
151   - const GEvent* m_event; //!< Pointer to event
152   - int m_ipar; //!< Parameter index
  149 + const GModel& m_model; //!< Reference to model
  150 + GModelPar* m_par; //!< Pointer to parameter
  151 + const GEvent& m_event; //!< Reference to event
153 152 };
154 153  
155 154 // Npred methods
... ... @@ -189,15 +188,15 @@ protected:
189 188 public:
190 189 npred_func(const GObservation* parent,
191 190 const GModel& model,
192   - int ipar) :
  191 + GModelPar* par) :
193 192 m_parent(parent),
194   - m_model(&model),
195   - m_ipar(ipar) { }
  193 + m_model(model),
  194 + m_par(par) { }
196 195 double eval(const double& x);
197 196 protected:
198 197 const GObservation* m_parent; //!< Pointer to parent
199   - const GModel* m_model; //!< Pointer to model
200   - int m_ipar; //!< Parameter index
  198 + const GModel& m_model; //!< Reference to model
  199 + GModelPar* m_par; //!< Pointer to parameter
201 200 };
202 201  
203 202 // Protected data area
... ...
pyext/GObservation.i
  View file @f949a67
... ... @@ -70,6 +70,22 @@ public:
70 70 virtual void read(const GXmlElement& xml) = 0;
71 71 virtual void write(GXmlElement& xml) const = 0;
72 72  
  73 + // Virtual methods
  74 + virtual double likelihood(const GModels& models,
  75 + GVector* gradient,
  76 + GMatrixSparse* curvature,
  77 + double* npred) const;
  78 + virtual double model(const GModels& models,
  79 + const GEvent& event,
  80 + GVector* gradient = NULL) const;
  81 + virtual double npred(const GModels& models,
  82 + GVector* gradient = NULL) const;
  83 + virtual double model_grad(const GModel& model,
  84 + const GModelPar& par,
  85 + const GEvent& event) const;
  86 + virtual double npred_grad(const GModel& model,
  87 + const GModelPar& par) const;
  88 +
73 89 // Implemented methods
74 90 void name(const std::string& name);
75 91 void id(const std::string& id);
... ... @@ -79,21 +95,6 @@ public:
79 95 const std::string& id(void) const;
80 96 const GEvents* events(void) const;
81 97 const std::string& statistics(void) const;
82   - virtual double likelihood(const GModels& models,
83   - GVector* gradient,
84   - GMatrixSparse* curvature,
85   - double* npred) const;
86   - virtual double model(const GModels& models,
87   - const GEvent& event,
88   - GVector* gradient = NULL) const;
89   - virtual double npred(const GModels& models,
90   - GVector* gradient = NULL) const;
91   - virtual double model_grad(const GModel& model,
92   - const GModelPar& par,
93   - const int& ipar,
94   - const GEvent& event) const;
95   - virtual double npred_grad(const GModel& model,
96   - const int& ipar) const;
97 98 };
98 99  
99 100  
... ...
src/obs/GObservation.cpp
  View file @f949a67
... ... @@ -67,8 +67,6 @@ const double minerr = 1.0e-100; //!&lt; Minimum statistical error
67 67  
68 68 /* __ Coding definitions _________________________________________________ */
69 69 #define G_LN_ENERGY_INT //!< ln(E) variable substitution for integration
70   -//#define G_GRAD_RIDDLER //!< Use Riddler's method for computing derivatives
71   -#define G_TEST_CODE
72 70  
73 71 /* __ Debug definitions __________________________________________________ */
74 72 //#define G_OPT_DEBUG
... ... @@ -335,7 +333,7 @@ double GObservation::model(const GModels&amp; models, const GEvent&amp; event,
335 333 (*gradient)[igrad+ipar] = par.factor_gradient();
336 334 }
337 335 else {
338   - (*gradient)[igrad+ipar] = model_grad(*mptr, par, ipar, event);
  336 + (*gradient)[igrad+ipar] = model_grad(*mptr, par, event);
339 337 }
340 338 }
341 339 else {
... ... @@ -415,7 +413,8 @@ double GObservation::npred(const GModels&amp; models, GVector* gradient) const
415 413 // Optionally determine Npred gradients
416 414 if (gradient != NULL) {
417 415 for (int k = 0; k < mptr->size(); ++k) {
418   - (*gradient)[igrad+k] = npred_grad(*mptr, k);
  416 + const GModelPar& par = (*mptr)[k];
  417 + (*gradient)[igrad+k] = npred_grad(*mptr, par);
419 418 }
420 419 }
421 420  
... ... @@ -437,36 +436,23 @@ double GObservation::npred(const GModels&amp; models, GVector* gradient) const
437 436 * @brief Returns parameter gradient of model for a given event
438 437 *
439 438 * @param[in] model Model.
  439 + * @param[in] par Model parameter.
440 440 * @param[in] event Event.
441   - * @param[in] ipar Parameter index for which gradient should be returned.
442 441 *
443   - * This method uses a robust but dumb method to estimate parameter
444   - * gradients that have not been provided by the model. We use here a dumb
445   - * method as this method is likely used for spatial model parameters, and
446   - * the spatial model may eventually be noisy due to numerical integration
  442 + * This method uses a robust but simple difference method to estimate
  443 + * parameter gradients that have not been provided by the model. We use here
  444 + * a simple method as this method is likely used for spatial model parameters,
  445 + * and the spatial model may eventually be noisy due to numerical integration
447 446 * limits.
448 447 *
449   - * The step size for the dumb method has been fixed to 0.0002, which
450   - * corresponds to abound 1 arcsec for parameters that are given in degrees.
  448 + * The step size for the simple method has been fixed to 0.0002, which
  449 + * corresponds to about 1 arcsec for parameters that are given in degrees.
451 450 * The reasoning behind this value is that parameters that use numerical
452 451 * gradients are typically angles, such as for example the position, and
453 452 * we want to achieve arcsec precision with this method.
454   - *
455   - * @todo Implement a more precise numerical derivation scheme. This needs
456   - * a deep investigation as the scheme needs to be precise but also
457   - * robust (not sensitive to noise in the function).
458   - * @todo For the Riddler method, we simply remove any parameter boundaries
459   - * here for the computation to avoid any out of boundary errors.
460   - * We may have models, however, for which out of bound parameters lead
461   - * to illegal computations, such as division by zero or taking the
462   - * square root of negative values.
463   - * I cannot see any elegant method to catch this at this level.
464   - * Eventually, the higher level method should avoid going in a
465   - * parameter domain that is not defined.
466 453 ***************************************************************************/
467 454 double GObservation::model_grad(const GModel& model,
468 455 const GModelPar& par,
469   - const int& ipar,
470 456 const GEvent& event) const
471 457 {
472 458 // Initialise gradient
... ... @@ -485,58 +471,47 @@ double GObservation::model_grad(const GModel&amp; model,
485 471 double x = par.factor_value();
486 472  
487 473 // Set fixed step size for computation of derivative.
488   - // By default, the step size is fixed to 0.0002, but if this would
489   - // violate a boundary, dx is reduced accordingly. In case that x
490   - // is right on the boundary, x is displaced slightly from the
491   - // boundary to allow evaluation of the derivative.
492   - #if !defined(G_GRAD_RIDDLER)
493   - //const double step_size = 0.0002; // ~1 arcsec
494   - const double step_size = 0.001; // ~4 arcsec
495   - //const double step_size = 1.0e-5;
496   - double dx = step_size;
497   - if (par.has_min()) {
498   - double dx_min = x - par.factor_min();
499   - if (dx_min == 0.0) {
500   - dx = step_size * x;
501   - if (dx == 0.0) {
502   - dx = step_size;
503   - }
504   - x += dx;
505   - }
506   - else if (dx_min < dx) {
507   - dx = dx_min;
508   - }
509   - }
510   - if (par.has_max()) {
511   - double dx_max = par.factor_max() - x;
512   - if (dx_max == 0.0) {
513   - dx = step_size * x;
514   - if (dx == 0.0) {
515   - dx = step_size;
516   - }
517   - x -= dx;
518   - }
519   - else if (dx_max < dx) {
520   - dx = dx_max;
  474 + // By default, the step size is fixed to 0.0002.
  475 + const double step_size = 0.0002; // ~1 arcsec
  476 + double h = step_size;
  477 +
  478 + // Re-adjust the step-size h in case that the initial step size is
  479 + // larger than the allowed parameter range
  480 + if (par.has_min() && par.has_max()) {
  481 + double par_h = par.factor_max() - par.factor_min();
  482 + if (par_h < h) {
  483 + h = par_h;
521 484 }
522 485 }
523   - #endif
524 486  
525   - // Remove any boundaries to avoid limitations
526   - ptr->remove_range();
  487 + // Continue only if step size is positive
  488 + if (h > 0.0) {
527 489  
528   - // Setup derivative function
529   - GObservation::model_func function(this, model, event, ipar);
  490 + // Remove any boundaries to avoid limitations
  491 + //ptr->remove_range(); // Not needed in principle
530 492  
531   - // Get derivative. We use a fixed step size here that has been
532   - // checked on spatial parameters of models
533   - GDerivative derivative(&function);
534   - #if defined(G_GRAD_RIDDLER)
535   - grad = derivative.value(x);
536   - #else
537   - //grad = derivative.difference(x, dx);
538   - grad = derivative.smooth_robust(x, dx, 2, 5);
539   - #endif
  493 + // Setup derivative function
  494 + GObservation::model_func function(this, model, ptr, event);
  495 + GDerivative derivative(&function);
  496 +
  497 + // If we are too close to the minimum boundary use a right sided
  498 + // difference ...
  499 + if (par.has_min() && ((x-par.factor_min()) < h)) {
  500 + grad = derivative.right_difference(x, h);
  501 + }
  502 +
  503 + // ... otherwise if we are too close to the maximum boundary use
  504 + // a left sided difference ...
  505 + else if (par.has_max() && ((par.factor_max()-x) < h)) {
  506 + grad = derivative.left_difference(x, h);
  507 + }
  508 +
  509 + // ... otherwise use a symmetric difference
  510 + else {
  511 + grad = derivative.difference(x, h);
  512 + }
  513 +
  514 + } // endif: step size was positive
540 515  
541 516 // Restore current model parameter
542 517 *ptr = current;
... ... @@ -552,7 +527,7 @@ double GObservation::model_grad(const GModel&amp; model,
552 527 * @brief Returns parameter gradient of Npred
553 528 *
554 529 * @param[in] model Gamma-ray source model.
555   - * @param[in] ipar Parameter index for which gradient should be returned.
  530 + * @param[in] par Model parameter.
556 531 *
557 532 * Computes
558 533 * \f[\frac{{\rm d} N_{\rm pred}}{{\rm d} a_i}\f]
... ... @@ -574,97 +549,80 @@ double GObservation::model_grad(const GModel&amp; model,
574 549 * \f$t\f$ is the true photon arrival time, and
575 550 * \f$d\f$ is the instrument pointing.
576 551 *
577   - * This method uses a robust but dumb method to estimate gradients. This
578   - * method has turned out more robust then the Riddler's method implement
579   - * by the GDerivative::value() method.
  552 + * This method uses a robust but simple difference method to estimate
  553 + * parameter gradients that have not been provided by the model. We use here
  554 + * a simple method as this method is likely used for spatial model parameters,
  555 + * and the spatial model may eventually be noisy due to numerical integration
  556 + * limits.
580 557 *
581   - * The step size for the dumb method has been fixed to 0.0002, which
582   - * corresponds to abound 1 arcsec for parameters that are given in degrees.
  558 + * The step size for the simple method has been fixed to 0.0002, which
  559 + * corresponds to about 1 arcsec for parameters that are given in degrees.
583 560 * The reasoning behind this value is that parameters that use numerical
584 561 * gradients are typically angles, such as for example the position, and
585 562 * we want to achieve arcsec precision with this method.
586   - *
587   - * @todo Implement a more precise numerical derivation scheme. This needs
588   - * a deep investigation as the scheme needs to be precise but also
589   - * robust (not sensitive to noise in the function).
590   - * @todo For Riddler's method we simply remove any parameter boundaries here
591   - * for the computation to avoid any out of boundary errors. We may have
592   - * models, however, for which out of bound parameters lead to illegal
593   - * computations, such as division by zero or taking the square root of
594   - * negative values.
595   - * I cannot see any elegant method to catch this at this level.
596   - * Eventually, the higher level method should avoid going in a
597   - * parameter domain that is not defined.
598 563 ***************************************************************************/
599   -double GObservation::npred_grad(const GModel& model, const int& ipar) const
  564 +double GObservation::npred_grad(const GModel& model, const GModelPar& par) const
600 565 {
601 566 // Initialise result
602 567 double grad = 0.0;
603 568  
604 569 // Compute gradient only if parameter is free
605   - if (model[ipar].is_free()) {
  570 + if (par.is_free()) {
606 571  
607   - // Get non-const model pointer (circumvent const correctness)
608   - GModel* ptr = const_cast<GModel*>(&model);
  572 + // Get non-const model pointer
  573 + GModelPar* ptr = const_cast<GModelPar*>(&par);
609 574  
610 575 // Save current model parameter
611   - GModelPar current = (*ptr)[ipar];
  576 + GModelPar current = par;
612 577  
613 578 // Get actual parameter value
614   - double x = model[ipar].factor_value();
615   -
616   - // Determine fixed step size for computation of derivative.
617   - // By default, the step size is fixed to 0.0002, but if this would
618   - // violate a boundary, dx is reduced accordingly. In case that x
619   - // is right on the boundary, x is displaced slightly from the
620   - // boundary to allow evaluation of the derivative.
621   - #if !defined(G_GRAD_RIDDLER)
622   - const double step_size = 0.0002;
623   - double dx = step_size;
624   - if (model[ipar].has_min()) {
625   - double dx_min = x - model[ipar].factor_min();
626   - if (dx_min == 0.0) {
627   - dx = step_size * x;
628   - if (dx == 0.0) {
629   - dx = step_size;
630   - }
631   - x += dx;
632   - }
633   - else if (dx_min < dx) {
634   - dx = dx_min;
635   - }
636   - }
637   - if (model[ipar].has_max()) {
638   - double dx_max = model[ipar].factor_max() - x;
639   - if (dx_max == 0.0) {
640   - dx = step_size * x;
641   - if (dx == 0.0) {
642   - dx = step_size;
643   - }
644   - x -= dx;
645   - }
646   - else if (dx_max < dx) {
647   - dx = dx_max;
  579 + double x = par.factor_value();
  580 +
  581 + // Set fixed step size for computation of derivative.
  582 + // By default, the step size is fixed to 0.0002.
  583 + const double step_size = 0.0002; // ~1 arcsec
  584 + double h = step_size;
  585 +
  586 + // Re-adjust the step-size h in case that the initial step size is
  587 + // larger than the allowed parameter range
  588 + if (par.has_min() && par.has_max()) {
  589 + double par_h = par.factor_max() - par.factor_min();
  590 + if (par_h < h) {
  591 + h = par_h;
648 592 }
649 593 }
650   - #endif
651 594  
652   - // Remove any boundaries to avoid limitations
653   - (*ptr)[ipar].remove_range();
  595 + // Continue only if step size is positive
  596 + if (h > 0.0) {
654 597  
655   - // Setup derivative function
656   - GObservation::npred_func function(this, model, ipar);
  598 + // Remove any boundaries to avoid limitations
  599 + //ptr->remove_range(); // Not needed in principle
657 600  
658   - // Get derivative.
659   - GDerivative derivative(&function);
660   - #if defined(G_GRAD_RIDDLER)
661   - grad = derivative.value(x);
662   - #else
663   - grad = derivative.difference(x, dx);
664   - #endif
  601 + // Setup derivative function
  602 + GObservation::npred_func function(this, model, ptr);
  603 + GDerivative derivative(&function);
  604 +
  605 + // If we are too close to the minimum boundary use a right sided
  606 + // difference ...
  607 + if (par.has_min() && ((x-par.factor_min()) < h)) {
  608 + grad = derivative.right_difference(x, h);
  609 + }
  610 +
  611 + // ... otherwise if we are too close to the maximum boundary use
  612 + // a left sided difference ...
  613 + else if (par.has_max() && ((par.factor_max()-x) < h)) {
  614 + grad = derivative.left_difference(x, h);
  615 + }
  616 +
  617 + // ... otherwise use a symmetric difference
  618 + else {
  619 + grad = derivative.difference(x, h);
  620 + }
  621 +
  622 + } // endif: step size was positive
665 623  
666 624 // Restore current model parameter
667   - (*ptr)[ipar] = current;
  625 + *ptr = current;
668 626  
669 627 } // endif: model parameter was free
670 628  
... ... @@ -1252,14 +1210,11 @@ double GObservation::likelihood_gaussian_binned(const GModels&amp; models,
1252 1210 ***************************************************************************/
1253 1211 double GObservation::model_func::eval(const double& x)
1254 1212 {
1255   - // Get non-const model pointer (circumvent const correctness)
1256   - GModel* model = const_cast<GModel*>(m_model);
1257   -
1258 1213 // Set value
1259   - (*model)[m_ipar].factor_value(x);
  1214 + m_par->factor_value(x);
1260 1215  
1261 1216 // Compute model value
1262   - double value = model->eval(*m_event, *m_parent);
  1217 + double value = m_model.eval(m_event, *m_parent);
1263 1218  
1264 1219 // Return value
1265 1220 return value;
... ... @@ -1279,14 +1234,11 @@ double GObservation::model_func::eval(const double&amp; x)
1279 1234 ***************************************************************************/
1280 1235 double GObservation::npred_func::eval(const double& x)
1281 1236 {
1282   - // Get non-const model pointer (circumvent const correctness)
1283   - GModel* model = const_cast<GModel*>(m_model);
1284   -
1285 1237 // Set value
1286   - (*model)[m_ipar].factor_value(x);
  1238 + m_par->factor_value(x);
1287 1239  
1288 1240 // Compute Npred value
1289   - double npred = m_parent->npred_temp(*model);
  1241 + double npred = m_parent->npred_temp(m_model);
1290 1242  
1291 1243 // Return value
1292 1244 return npred;
... ...