Jürgen Knödlseder / gammalib

Commit 0125a7f64f495f5fbfba93f0739d4c3d29d11d08

Authored by knodlseder
1 parent 7e8be82b

Redcuce fixed step size of numerical derivative computation in GObservation::mod… 

…el_grad and GObservation::npred_grad to 0.0002.

This value is equivalent to 1 arcsec for angular parameters, and improves the precision of the position determination.
Inline Side-by-side
Showing 1 changed file with 33 additions and 13 deletions
src/obs/GObservation.cpp
  View file @0125a7f
... ... @@ -459,6 +459,15 @@ void GObservation::free_members(void)
459 459 * the spatial model may eventually be noisy due to numerical integration
460 460 * limits.
461 461 *
  462 + * The step size for the dumb method has been fixed to 0.0002, which
  463 + * corresponds to abound 1 arcsec for parameters that are given in degrees.
  464 + * The reasoning behind this value is that parameters that use numerical
  465 + * gradients are typically angles, such as for example the position, and
  466 + * we want to achieve arcsec precision with this method.
  467 + *
  468 + * @todo Implement a more precise numerical derivation scheme. This needs
  469 + * a deep investigation as the scheme needs to be precise but also
  470 + * robust (not sensitive to noise in the function).
462 471 * @todo For the Riddler method, we simply remove any parameter boundaries
463 472 * here for the computation to avoid any out of boundary errors.
464 473 * We may have models, however, for which out of bound parameters lead
... ... @@ -493,19 +502,20 @@ double GObservation::model_grad(const GModel& model, const GEvent& event,
493 502 // Get actual parameter value
494 503 double x = model[ipar].value();
495 504  
496   - // Determine fixed step size for computation of derivative.
497   - // By default, the step size is fixed to 0.05, but if this would
  505 + // Set fixed step size for computation of derivative.
  506 + // By default, the step size is fixed to 0.0002, but if this would
498 507 // violate a boundary, dx is reduced accordingly. In case that x
499 508 // is right on the boundary, x is displaced slightly from the
500 509 // boundary to allow evaluation of the derivative.
501 510 #if !defined(G_GRAD_RIDDLER)
502   - double dx = 0.05;
  511 + const double step_size = 0.0002;
  512 + double dx = step_size;
503 513 if (model[ipar].hasmin()) {
504 514 double dx_min = x - model[ipar].min();
505 515 if (dx_min == 0.0) {
506   - dx = 0.05 * x;
  516 + dx = step_size * x;
507 517 if (dx == 0.0) {
508   - dx = 0.05;
  518 + dx = step_size;
509 519 }
510 520 x += dx;
511 521 }
... ... @@ -516,9 +526,9 @@ double GObservation::model_grad(const GModel& model, const GEvent& event,
516 526 if (model[ipar].hasmax()) {
517 527 double dx_max = model[ipar].max() - x;
518 528 if (dx_max == 0.0) {
519   - dx = 0.05 * x;
  529 + dx = step_size * x;
520 530 if (dx == 0.0) {
521   - dx = 0.05;
  531 + dx = step_size;
522 532 }
523 533 x -= dx;
524 534 }
... ... @@ -612,6 +622,15 @@ double GObservation::model_func::eval(double x)
612 622 * method has turned out more robust then the Riddler's method implement
613 623 * by the GDerivative::value() method.
614 624 *
  625 + * The step size for the dumb method has been fixed to 0.0002, which
  626 + * corresponds to abound 1 arcsec for parameters that are given in degrees.
  627 + * The reasoning behind this value is that parameters that use numerical
  628 + * gradients are typically angles, such as for example the position, and
  629 + * we want to achieve arcsec precision with this method.
  630 + *
  631 + * @todo Implement a more precise numerical derivation scheme. This needs
  632 + * a deep investigation as the scheme needs to be precise but also
  633 + * robust (not sensitive to noise in the function).
615 634 * @todo For Riddler's method we simply remove any parameter boundaries here
616 635 * for the computation to avoid any out of boundary errors. We may have
617 636 * models, however, for which out of bound parameters lead to illegal
... ... @@ -639,18 +658,19 @@ double GObservation::npred_grad(const GModel& model, int ipar) const
639 658 double x = model[ipar].value();
640 659  
641 660 // Determine fixed step size for computation of derivative.
642   - // By default, the step size is fixed to 0.05, but if this would
  661 + // By default, the step size is fixed to 0.0002, but if this would
643 662 // violate a boundary, dx is reduced accordingly. In case that x
644 663 // is right on the boundary, x is displaced slightly from the
645 664 // boundary to allow evaluation of the derivative.
646 665 #if !defined(G_GRAD_RIDDLER)
647   - double dx = 0.05;
  666 + const double step_size = 0.0002;
  667 + double dx = step_size;
648 668 if (model[ipar].hasmin()) {
649 669 double dx_min = x - model[ipar].min();
650 670 if (dx_min == 0.0) {
651   - dx = 0.05 * x;
  671 + dx = step_size * x;
652 672 if (dx == 0.0) {
653   - dx = 0.05;
  673 + dx = step_size;
654 674 }
655 675 x += dx;
656 676 }
... ... @@ -661,9 +681,9 @@ double GObservation::npred_grad(const GModel& model, int ipar) const
661 681 if (model[ipar].hasmax()) {
662 682 double dx_max = model[ipar].max() - x;
663 683 if (dx_max == 0.0) {
664   - dx = 0.05 * x;
  684 + dx = step_size * x;
665 685 if (dx == 0.0) {
666   - dx = 0.05;
  686 + dx = step_size;
667 687 }
668 688 x -= dx;
669 689 }
... ...