Domenico Tiziani / gammalib

Commit cd76bca36751ca26ff6901c3212ca35e127b0f64

Authored by knodlseder
1 parent f47ed897

add difference() method; minuit2 not necessarily operational ... (still needs some fix of eps)

Inline Side-by-side
Showing 2 changed files with 167 additions and 23 deletions
include/GDerivative.hpp
  View file @cd76bca
... ... @@ -68,19 +68,33 @@ public:
68 68 const GFunction* function(void) const { return m_func; }
69 69 double value(const double& x, double step = 0.0);
70 70 double ridder(const double& x, const double& h, double& err);
71   - double minuit(const double& x, double& err);
  71 + double minuit2(const double& x, double& err);
  72 + double difference(const double& x, const double& h);
72 73 std::string print(void) const;
73 74  
74 75 protected:
75 76 // Protected methods
76   - void init_members(void);
77   - void copy_members(const GDerivative& dx);
78   - void free_members(void);
  77 + void init_members(void);
  78 + void copy_members(const GDerivative& dx);
  79 + void free_members(void);
  80 + void set_tiny(void);
  81 +
  82 + // Tiny number computation
  83 + class tiny {
  84 + public:
  85 + tiny(void) : m_one(1.0) {}
  86 + ~tiny(void) {}
  87 + double one(void) const;
  88 + double operator()(double eps) const;
  89 + private:
  90 + double m_one;
  91 + };
79 92  
80 93 // Protected members
81 94 GFunction* m_func; //!< Pointer to function
82 95 double m_eps; //!< Derivative precision
83 96 double m_step_frac; //!< Value fraction to use for initial step
  97 + double m_tiny; //!< Tiny number for minuit2
84 98 int m_max_iter; //!< Maximum number of iterations
85 99 int m_iter; //!< Number of iterations used
86 100 bool m_silent; //!< Suppress warnings
... ...
src/numerics/GDerivative.cpp
  View file @cd76bca
... ... @@ -33,7 +33,7 @@
33 33  
34 34 /* __ Debug definitions __________________________________________________ */
35 35 #define G_VALUE_DEBUG 0 //!< Debug value method
36   -#define G_MINUIT_DEBUG 1 //!< Debug minuit method
  36 +#define G_MINUIT_DEBUG 0 //!< Debug minuit method
37 37  
38 38  
39 39 /*==========================================================================
... ... @@ -333,7 +333,7 @@ double GDerivative::ridder(const double&amp; x, const double&amp; h, double&amp; err)
333 333  
334 334 } // endfor: looped over tableau
335 335  
336   - // If higher order is worse by a significant factor, the quit early
  336 + // If higher order is worse by a significant factor, then quit early
337 337 if (std::abs( a(i,i) - a(i-1,i-1) ) >= safe*err)
338 338 break;
339 339  
... ... @@ -347,57 +347,90 @@ double GDerivative::ridder(const double&amp; x, const double&amp; h, double&amp; err)
347 347  
348 348  
349 349 /***********************************************************************//**
350   - * @brief Returns derivative using minuit2 algorithm
  350 + * @brief Returns derivative using Minuit2 algorithm
351 351 *
352 352 * @param[in] x Function value.
353 353 * @param[out] err Estimated error in the derivative.
  354 + *
  355 + * This method has been inspired by corresponding code in the Minuit2
  356 + * package. Please check the following files in Minuit2 for more information:
  357 + * Numerical2PGradientCalculator.cxx for the main code
  358 + * MnStrategy.cxx for the definition of algorithm parameters
  359 + * MnMachinePrecision.h for the definition of eps and eps2
  360 + * InitialGradientCalculator.cxx for the Minuit2 way of estimating initial
  361 + * gradients.
  362 + *
  363 + * The exact value of fcnup has little impact on the results.
354 364 ***************************************************************************/
355   -double GDerivative::minuit(const double& x, double& err)
  365 +double GDerivative::minuit2(const double& x, double& err)
356 366 {
357 367 // Evaluate function at x
358 368 double fcnmin = m_func->eval(x);
359 369  
360 370 // Set precision
361   - double eps = m_eps; // Minuit2: smallest number that gives 1.+eps>1
  371 +// double eps = m_eps; // Minuit2: smallest number that gives 1.+eps>1
  372 +// double eps = m_tiny; // Bad
  373 +// double eps = 1.0e-30; // Bad
  374 +// double eps = 1.0e-10; // Poor
  375 +// double eps = 1.0e-6; // Quite good
  376 +// double eps = 1.0e-8; // Worse than 1.0e-6
  377 + double eps = 1.0e-4; // Pretty good !
  378 +// double eps = 1.0e-2; // Wrong
362 379 double eps2 = 2.0 * std::sqrt(eps);
363 380  
364 381 // Set ...
365   - double fcnup = 1.0; // Minuit2: Fcn().Up()
366   - double dfmin = 8.0 * eps2 * (std::abs(fcnmin) + fcnup);
367   - double vrysml = 8.0 * eps * eps;
  382 + //double fcnup = 1.0; // Minuit2: Fcn().Up() - function delta for 1 sigma errors
  383 + double fcnup = 0.01 * fcnmin;
  384 + double dfmin = 8.0 * eps2 * (std::abs(fcnmin) + fcnup);
  385 + double vrysml = 8.0 * eps * eps;
368 386  
369 387 // Set high strategy (see MnStrategy)
370 388 int ncycles = 5;
371 389 double step_tol = 0.1;
372 390 double grad_tol = 0.02;
373 391  
374   - // Initialise values
375   - double epspri = eps2; // Minuit2: eps2 + fabs(grd(i)*eps2)
376   - double stepb4 = 0.0; //
377   - err = 0.0; // Initial error estimate
378   -
379   - // Initial gradient estimation
380   - double gstep = 0.1; // Initial gradient step
  392 + // Initial gradient estimation. Note that this is brute force method
  393 + // that does not correspond to the Minuit2 approach. For the Minuit2
  394 + // approach, check InitialGradientCalculator.cxx
  395 + double gstep = 0.001; // Initial gradient step
381 396 double fs1 = m_func->eval(x + gstep);
382 397 double fs2 = m_func->eval(x - gstep);
383 398 double grd = 0.5 * (fs1 - fs2) / gstep;
384 399 double g2 = (fs1 + fs2 - 2.0*fcnmin) / gstep / gstep;
385   -
  400 +
  401 + // Debug option: Show actual results
  402 + #if G_MINUIT_DEBUG
  403 + std::cout << "GDerivative::minuit(";
  404 + std::cout << "iter=0";
  405 + std::cout << ", x=" << x;
  406 + std::cout << ", grd=" << grd;
  407 + std::cout << ", step=" << gstep;
  408 + std::cout << ")" << std::endl;
  409 + #endif
  410 +
  411 + // Initialise values
  412 + double epspri = eps2 + std::abs(grd*eps2);
  413 + double step = 0.0;
  414 + double stepb4 = 0.0;
  415 + double step_change = 0.0;
  416 + err = 0.0;
  417 +
386 418 // Loop over cycles
387   - for (m_iter = 0; m_iter < ncycles; ++m_iter) {
  419 + for (m_iter = 1; m_iter <= ncycles; ++m_iter) {
388 420  
389 421 // Compute optimum step size
390 422 double optstp = std::sqrt(dfmin/(std::abs(g2)+epspri));
391 423  
392 424 // Compute step size
393   - double step = std::max(optstp, std::abs(0.1*gstep));
  425 + step = std::max(optstp, std::abs(0.1*gstep));
394 426 double stpmax = 10.0 * std::abs(gstep);
395 427 double stpmin = std::max(vrysml, 8.0 * std::abs(eps2*x));
396 428 if (step > stpmax) step = stpmax;
397 429 if (step < stpmin) step = stpmin;
398 430  
399 431 // Break is we are below the step tolerance
400   - if (std::abs((step-stepb4)/step) < step_tol)
  432 + step_change = std::abs((step-stepb4)/step);
  433 + if (step_change < step_tol)
401 434 break;
402 435  
403 436 // Store step size information
... ... @@ -433,6 +466,43 @@ double GDerivative::minuit(const double&amp; x, double&amp; err)
433 466  
434 467 } // endfor: looped over cycles
435 468  
  469 + // Compile option: signal if we exceed the tolerance
  470 + if (!silent()) {
  471 + if (err >= grad_tol) {
  472 + std::cout << "WARNING: GDerivative::minuit(";
  473 + std::cout << "iter=" << m_iter;
  474 + std::cout << ", x=" << x;
  475 + std::cout << ", grd=" << grd;
  476 + std::cout << ", step=" << step;
  477 + std::cout << ", step_change=" << step_change;
  478 + std::cout << ", err=" << err;
  479 + std::cout << "): error exceeds tolerance of ";
  480 + std::cout << grad_tol << std::endl;
  481 + }
  482 + }
  483 +
  484 + // Return gradient
  485 + return grd;
  486 +}
  487 +
  488 +
  489 +/***********************************************************************//**
  490 + * @brief Returns gradient computed from function difference
  491 + *
  492 + * @param[in] x Function value.
  493 + * @param[in] h Step size.
  494 + *
  495 + * This is the most simple and dumb gradient computation method we can think
  496 + * of. It does a reasonable good job if the problem is well controlled, and
  497 + * in particular, if the step size is known in advance.
  498 + ***************************************************************************/
  499 +double GDerivative::difference(const double& x, const double& h)
  500 +{
  501 + // Compute gradient
  502 + double fs1 = m_func->eval(x + h);
  503 + double fs2 = m_func->eval(x - h);
  504 + double grd = 0.5 * (fs1 - fs2) / h;
  505 +
436 506 // Return gradient
437 507 return grd;
438 508 }
... ... @@ -478,10 +548,14 @@ void GDerivative::init_members(void)
478 548 m_func = NULL;
479 549 m_eps = 1.0e-6;
480 550 m_step_frac = 0.02;
  551 + m_tiny = 1.0e-7;
481 552 m_max_iter = 5;
482 553 m_iter = 0;
483 554 m_silent = false;
484 555  
  556 + // Compute machine precision
  557 + //set_tiny();
  558 +
485 559 // Return
486 560 return;
487 561 }
... ... @@ -498,6 +572,7 @@ void GDerivative::copy_members(const GDerivative&amp; dx)
498 572 m_func = dx.m_func;
499 573 m_eps = dx.m_eps;
500 574 m_step_frac = dx.m_step_frac;
  575 + m_tiny = dx.m_tiny;
501 576 m_max_iter = dx.m_max_iter;
502 577 m_iter = dx.m_iter;
503 578 m_silent = dx.m_silent;
... ... @@ -517,6 +592,61 @@ void GDerivative::free_members(void)
517 592 }
518 593  
519 594  
  595 +/***********************************************************************//**
  596 + * @brief Compute tiny number for Minuit2
  597 + *
  598 + * On fermi gives m_tiny = 8.88178e-16 for i=51.
  599 + ***************************************************************************/
  600 +void GDerivative::set_tiny(void)
  601 +{
  602 + // Allocate tiny instance
  603 + GDerivative::tiny tiny;
  604 +
  605 + // Calculate machine precision
  606 + double epstry = 0.5;
  607 + double epsbak = 0.0;
  608 + double epsp1 = 0.0;
  609 + double one = 1.0;
  610 +
  611 + // Loop until we found
  612 + for (int i = 0; i < 100; ++i) {
  613 + epstry *= 0.5;
  614 + epsp1 = one + epstry;
  615 + epsbak = tiny(epsp1);
  616 + if (epsbak < epstry) {
  617 + m_tiny = 8.0 * epstry;
  618 + break;
  619 + }
  620 + }
  621 +
  622 + // Return
  623 + return;
  624 +}
  625 +
  626 +
  627 +/***********************************************************************//**
  628 + * @brief Return one
  629 + ***************************************************************************/
  630 +double GDerivative::tiny::one(void) const
  631 +{
  632 + // Return
  633 + return m_one;
  634 +}
  635 +
  636 +
  637 +/***********************************************************************//**
  638 + * @brief Evaluate minimal difference between two floating points
  639 + ***************************************************************************/
  640 +double GDerivative::tiny::operator()(double eps) const
  641 +{
  642 + // Compute difference
  643 + double result = eps - one();
  644 +
  645 + // Return result
  646 + return result;
  647 +}
  648 +
  649 +
520 650 /*==========================================================================
521 651 = =
522 652 = Friends =
... ...