GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GDerivative.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GDerivative.cpp - Derivative class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2011-2014 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GDerivative.cpp
23  * @brief GDerivative class implementation.
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #include <cstdlib> // Definition of NULL
29 #include <cmath> // for abs()
30 #include <cfloat> // for DBL_MAX
31 #include "GDerivative.hpp"
32 #include "GMatrix.hpp"
33 #include "GTools.hpp"
34 
35 /* __ Method name definitions ____________________________________________ */
36 #define G_VALUE "GDerivative::value(double&, double&)"
37 #define G_RIDDER "GDerivative::ridder(double&, double&, double*)"
38 #define G_MINUIT2 "GDerivative::minuit2(double&, double*)"
39 #define G_SMOOTH_ROBUST "GDerivative::smooth_robust(double&, double&, int&,"\
40  " int&)"
41 
42 /* __ Macros _____________________________________________________________ */
43 
44 /* __ Coding definitions _________________________________________________ */
45 
46 /* __ Debug definitions __________________________________________________ */
47 #define G_VALUE_DEBUG 0 //!< Debug value method
48 #define G_MINUIT_DEBUG 0 //!< Debug minuit method
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  ***************************************************************************/
61 {
62  // Initialise private members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief Function constructor
72  *
73  * @param[in] func Function pointer.
74  *
75  * Construct class instance from single parameter function. Note that the
76  * function is not copied but just a pointer to it is stored. The function
77  * should thus not be destroyed before computation of derivatives is
78  * finished.
79  ***************************************************************************/
81 {
82  // Initialise private members
83  init_members();
84 
85  // Set function
86  this->function(func);
87 
88  // Return
89  return;
90 }
91 
92 
93 /***********************************************************************//**
94  * @brief Copy constructor
95  *
96  * @param[in] dx Derivative.
97  ***************************************************************************/
99 {
100  // Initialise members
101  init_members();
102 
103  // Copy members
104  copy_members(dx);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Destructor
113  ***************************************************************************/
115 {
116  // Free members
117  free_members();
118 
119  // Return
120  return;
121 }
122 
123 
124 /*==========================================================================
125  = =
126  = Operators =
127  = =
128  ==========================================================================*/
129 
130 /***********************************************************************//**
131  * @brief Assignment operator
132  *
133  * @param[in] dx Derivative.
134  * @return Derivative.
135  ***************************************************************************/
137 {
138  // Execute only if object is not identical
139  if (this != &dx) {
140 
141  // Free members
142  free_members();
143 
144  // Initialise members
145  init_members();
146 
147  // Copy members
148  copy_members(dx);
149 
150  } // endif: object was not identical
151 
152  // Return
153  return *this;
154 }
155 
156 
157 /*==========================================================================
158  = =
159  = Public methods =
160  = =
161  ==========================================================================*/
162 
163 /***********************************************************************//**
164  * @brief Clear derivative
165  ***************************************************************************/
167 {
168  // Free members
169  free_members();
170 
171  // Initialise private members
172  init_members();
173 
174  // Return
175  return;
176 }
177 
178 
179 /***********************************************************************//**
180  * @brief Clone derivative
181  *
182  * @return Pointer to deep copy of derivative.
183  ***************************************************************************/
185 {
186  return new GDerivative(*this);
187 }
188 
189 
190 /***********************************************************************//**
191  * @brief Returns derivative
192  *
193  * @param[in] x Function value.
194  * @param[in] step Initial step size (default: automatic)
195  *
196  * Computes derivative using the Ridders' method. If the specified initial
197  * step size is non-zero, the method starts with this step size and
198  * iteratively decreases the step size until the error becomes less than
199  * a threshold defined with the eps() method (by default the threshold is
200  * set to 1e-6). If the initial step size is zero, the method takes
201  * \f$h={\tt m\_step\_frac} |x|\f$ as initial step size.
202  * The maximum number of iterations is controlled by the max_iter()
203  * method (by default, the maximum is set to 5).
204  ***************************************************************************/
205 double GDerivative::value(const double& x, const double& step)
206 {
207  // Set constants
208  const double min_h = 1.0e-6;
209 
210  // Initialise result
211  double result = 0.0;
212 
213  // Set initial step size. Use either the specified initial step size,
214  // or if the initial step size is 0, use a fixed fraction of the function
215  // value
216  double h = step;
217  if (h == 0.0) {
218  h = m_step_frac * std::abs(x);
219  if (h < min_h) {
220  h = min_h;
221  }
222  }
223 
224  // Initialise tolerance
225  double err = DBL_MAX;
226 
227  // Loop over Ridder's method until we have an acceptable tolerance
228  m_iter = 0;
229  for (; m_iter < m_max_iter; ++m_iter) {
230 
231  // Compute derivative using Ridder's method
232  result = ridder(x, h, &err);
233 
234  // Debug option: Show actual results
235  #if G_VALUE_DEBUG
236  std::cout << "GDerivative::value(";
237  std::cout << "iter=" << m_iter;
238  std::cout << ", x=" << x;
239  std::cout << ", dy/dx=" << result;
240  std::cout << ", h=" << h;
241  std::cout << ", err=" << err;
242  std::cout << ")" << std::endl;
243  #endif
244 
245  // If uncertainty is below tolerance then exit now
246  if (err < m_eps) {
247  break;
248  }
249 
250  // Reduce the step size
251  h /= 5.0;
252 
253  } // endfor: tolerance matching loop
254 
255  // Compile option: signal if we exceed the tolerance
256  if (!silent()) {
257  if (err >= m_eps) {
258  std::string msg = "Derivative uncertainty "+gammalib::str(err)+
259  " exceeds tolerance of "+gammalib::str(m_eps)+
260  " at function value "+gammalib::str(x)+
261  " (df/dx="+gammalib::str(result)+
262  ", iter="+gammalib::str(m_iter)+
263  ", h="+gammalib::str(h)+").";
265  }
266  }
267 
268  // Return derivative
269  return result;
270 }
271 
272 /***********************************************************************//**
273  * @brief Returns derivative by Ridders' method
274  *
275  * @param[in] x Function value.
276  * @param[in] h Estimated initial step size.
277  * @param[out] err Estimated error in the derivative.
278  *
279  * @exception GException::invalid_argument
280  * Step size of 0 has been specified.
281  *
282  * Returns the derivative of a function at point x by Ridders's method of
283  * polynomial extrapolation. The value h is input as an estimated initial
284  * stepsize; it needs not to be small, but rather should be an increment
285  * in x over which the function changes substantially.
286  *
287  * This method has been adopted from Numerical Recipes, 3rd edition, page
288  * 321f.
289  ***************************************************************************/
290 double GDerivative::ridder(const double& x, const double& h, double* err)
291 {
292  // Set constants
293  const int ntab = 10; // Maximum table size
294  const double con = 1.4; // Stepsize decreases by con at each iteration
295  const double con2 = con * con;
296  const double big = DBL_MAX;
297  const double safe = 2.0; // Return when error is safe worse than best so far
298 
299  // Initialise result
300  double dx = 0.0;
301  *err = big;
302 
303  // Continue only if we have a valid function
304  if (m_func != NULL) {
305 
306  // Allocate Neville table
307  GMatrix a(ntab, ntab);
308 
309  // Throw exception if step size is zero
310  if (h <= 0) {
311  throw GException::invalid_argument(G_RIDDER, "Step size must be positive.");
312  }
313 
314  // Initialise adaptive step size
315  double hh = h;
316 
317  // Perform initial function evaluation
318  a(0,0) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh);
319 
320  // Perform successive computations
321  for (int i = 1; i < ntab; ++i) {
322 
323  // Reduce step size
324  hh /= con;
325 
326  // Try new smaller step size
327  a(0,i) = (m_func->eval(x+hh) - m_func->eval(x-hh))/(2.0*hh);
328 
329  // Compute extrapolations of various orders, requiring no new
330  // function evaluations
331  double fac = con2;
332  for (int j = 1; j <= i; ++j) {
333 
334  // Compute extrapolation
335  a(j,i) = ( a(j-1,i) * fac - a(j-1,i-1) ) / (fac - 1.0);
336  fac *= con2;
337 
338  // If the error is reduced then save the improved answer
339  double err1 = std::abs(a(j,i) - a(j-1,i));
340  double err2 = std::abs(a(j,i) - a(j-1,i-1));
341  double errt = (err1 > err2) ? err1 : err2;
342  if (errt <= *err) {
343  *err = errt;
344  dx = a(j,i);
345  }
346 
347  } // endfor: looped over tableau
348 
349  // If higher order is worse by a significant factor, then quit early
350  if (std::abs( a(i,i) - a(i-1,i-1) ) >= safe*(*err)) {
351  break;
352  }
353 
354  } // endfor: looped over computations
355 
356  } // endif: we had a valid function
357 
358  // Return derivative
359  return dx;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Returns derivative using Minuit2 algorithm
365  *
366  * @param[in] x Function value.
367  * @param[out] err Estimated error in the derivative.
368  *
369  * This method has been inspired by corresponding code in the Minuit2
370  * package. Please check the following files in Minuit2 for more information:
371  *
372  * Numerical2PGradientCalculator.cxx for the main code
373  * MnStrategy.cxx for the definition of algorithm parameters
374  * MnMachinePrecision.h for the definition of eps and eps2
375  * InitialGradientCalculator.cxx for the Minuit2 way of estimating initial
376  * gradients.
377  *
378  * The exact value of fcnup has little impact on the results.
379  ***************************************************************************/
380 double GDerivative::minuit2(const double& x, double* err)
381 {
382  // Evaluate function at x
383  double fcnmin = m_func->eval(x);
384 
385  // Set precision
386 // double eps = m_eps; // Minuit2: smallest number that gives 1.+eps>1
387 // double eps = m_tiny; // Bad
388 // double eps = 1.0e-30; // Bad
389 // double eps = 1.0e-10; // Poor
390 // double eps = 1.0e-6; // Quite good
391 // double eps = 1.0e-8; // Worse than 1.0e-6
392  double eps = 1.0e-4; // Pretty good !
393 // double eps = 1.0e-2; // Wrong
394  double eps2 = 2.0 * std::sqrt(eps);
395 
396  // Set ...
397  //double fcnup = 1.0; // Minuit2: Fcn().Up() - function delta for 1 sigma errors
398  double fcnup = 0.01 * fcnmin;
399  double dfmin = 8.0 * eps2 * (std::abs(fcnmin) + fcnup);
400  double vrysml = 8.0 * eps * eps;
401 
402  // Set high strategy (see MnStrategy)
403  int ncycles = 5;
404  double step_tol = 0.1;
405  double grad_tol = 0.02;
406 
407  // Initial gradient estimation. Note that this is brute force method
408  // that does not correspond to the Minuit2 approach. For the Minuit2
409  // approach, check InitialGradientCalculator.cxx
410  double gstep = 0.001; // Initial gradient step
411  double fs1 = m_func->eval(x + gstep);
412  double fs2 = m_func->eval(x - gstep);
413  double grd = 0.5 * (fs1 - fs2) / gstep;
414  double g2 = (fs1 + fs2 - 2.0*fcnmin) / gstep / gstep;
415 
416  // Debug option: Show actual results
417  #if G_MINUIT_DEBUG
418  std::cout << "GDerivative::minuit(";
419  std::cout << "iter=0";
420  std::cout << ", x=" << x;
421  std::cout << ", grd=" << grd;
422  std::cout << ", step=" << gstep;
423  std::cout << ")" << std::endl;
424  #endif
425 
426  // Initialise values
427  double epspri = eps2 + std::abs(grd*eps2);
428  double step = 0.0;
429  double stepb4 = 0.0;
430  double step_change = 0.0;
431  *err = 0.0;
432 
433  // Loop over cycles
434  for (m_iter = 1; m_iter <= ncycles; ++m_iter) {
435 
436  // Compute optimum step size
437  double optstp = std::sqrt(dfmin/(std::abs(g2)+epspri));
438 
439  // Compute step size
440  step = std::max(optstp, std::abs(0.1*gstep));
441  double stpmax = 10.0 * std::abs(gstep);
442  double stpmin = std::max(vrysml, 8.0 * std::abs(eps2*x));
443  if (step > stpmax) step = stpmax;
444  if (step < stpmin) step = stpmin;
445 
446  // Break is we are below the step tolerance
447  step_change = std::abs((step-stepb4)/step);
448  if (step_change < step_tol) {
449  break;
450  }
451 
452  // Store step size information
453  gstep = step;
454  stepb4 = step;
455 
456  // Bookkeeping of last gradient
457  double grdb4 = grd;
458 
459  // Evaluate gradient
460  double fs1 = m_func->eval(x + step);
461  double fs2 = m_func->eval(x - step);
462  grd = 0.5 * (fs1 - fs2) / step;
463  g2 = (fs1 + fs2 - 2.0*fcnmin) / step / step;
464 
465  // Compute error
466  *err = std::abs(grdb4-grd) / (std::abs(grd)+dfmin/step);
467 
468  // Debug option: Show actual results
469  #if G_MINUIT_DEBUG
470  std::cout << "GDerivative::minuit(";
471  std::cout << "iter=" << m_iter;
472  std::cout << ", x=" << x;
473  std::cout << ", grd=" << grd;
474  std::cout << ", step=" << step;
475  std::cout << ", err=" << *err;
476  std::cout << ")" << std::endl;
477  #endif
478 
479  // Break if gradient is accurate enough
480  if (*err < grad_tol) {
481  break;
482  }
483 
484  } // endfor: looped over cycles
485 
486  // Compile option: signal if we exceed the tolerance
487  if (!silent()) {
488  if (*err >= grad_tol) {
489  std::string msg = "Derivative uncertainty "+gammalib::str(*err)+
490  " exceeds tolerance of "+gammalib::str(grad_tol)+
491  " at function value "+gammalib::str(x)+
492  " (df/dx="+gammalib::str(grd)+
493  ", step="+gammalib::str(step)+
494  ", step_change="+gammalib::str(step_change)+").";
496  }
497  }
498 
499  // Return gradient
500  return grd;
501 }
502 
503 
504 /***********************************************************************//**
505  * @brief Returns gradient computed from symmetric function difference
506  *
507  * @param[in] x Function value.
508  * @param[in] h Step size.
509  *
510  * This is the most simple and dumb gradient computation method we can think
511  * of. It does a reasonable good job if the problem is well controlled, and
512  * in particular, if the step size is known in advance.
513  ***************************************************************************/
514 double GDerivative::difference(const double& x, const double& h)
515 {
516  // Compute derivative
517  double fs1 = m_func->eval(x + h);
518  double fs2 = m_func->eval(x - h);
519  double derivative = 0.5 * (fs1 - fs2) / h;
520 
521  // Return derivative
522  return derivative;
523 }
524 
525 
526 /***********************************************************************//**
527  * @brief Returns gradient computed from left-sided function difference
528  *
529  * @param[in] x Function value.
530  * @param[in] h Step size.
531  *
532  * This is the most simple and dumb gradient computation method we can think
533  * of. It does a reasonable good job if the problem is well controlled, and
534  * in particular, if the step size is known in advance.
535  ***************************************************************************/
536 double GDerivative::left_difference(const double& x, const double& h)
537 {
538  // Compute derivative
539  double fs1 = m_func->eval(x);
540  double fs2 = m_func->eval(x - h);
541  double derivative = (fs1 - fs2) / h;
542 
543  // Return derivative
544  return derivative;
545 }
546 
547 
548 /***********************************************************************//**
549  * @brief Returns gradient computed from right-sided function difference
550  *
551  * @param[in] x Function value.
552  * @param[in] h Step size.
553  *
554  * This is the most simple and dumb gradient computation method we can think
555  * of. It does a reasonable good job if the problem is well controlled, and
556  * in particular, if the step size is known in advance.
557  ***************************************************************************/
558 double GDerivative::right_difference(const double& x, const double& h)
559 {
560  // Compute derivative
561  double fs1 = m_func->eval(x + h);
562  double fs2 = m_func->eval(x);
563  double derivative = (fs1 - fs2) / h;
564 
565  // Return derivative
566  return derivative;
567 }
568 
569 
570 /***********************************************************************//**
571  * @brief Returns smooth noise-robust derivative
572  *
573  * @param[in] x Function value.
574  * @param[in] h Step size.
575  * @param[in] degree Differentiator degree (default: 2).
576  * @param[in] length Filter length (default: 5).
577  *
578  * Compute a smooth noise-robust numerical derivative using the method
579  * proposed by Pavel Holoborodko
580  * (http://www.holoborodko.com/pavel/numerical-methods/)
581  * Formulae for the following filter lengths and degrees have been
582  * implemented:
583  *
584  * degree=2 length=5,7,9,11
585  * degree=4 length=7,9,11
586  ***************************************************************************/
587 double GDerivative::smooth_robust(const double& x, const double& h,
588  const int& degree, const int& length)
589 {
590  // Initialise derivative
591  double derivative = 0.0;
592 
593  // Use method dependent on length and degree
594  switch (degree) {
595 
596  // Degree 2
597  case 2:
598  switch (length) {
599  case 5:
600  {
601  double h2 = 2.0 * h;
602  double fp1 = m_func->eval(x + h);
603  double fm1 = m_func->eval(x - h);
604  double fp2 = m_func->eval(x + h2);
605  double fm2 = m_func->eval(x - h2);
606  derivative = (2.0*(fp1-fm1) + (fp2-fm2))/(8.0*h);
607  } // block needed against cross-initialization
608  break;
609  case 7:
610  {
611  double h2 = 2.0 * h;
612  double h3 = 3.0 * h;
613  double fp1 = m_func->eval(x + h);
614  double fm1 = m_func->eval(x - h);
615  double fp2 = m_func->eval(x + h2);
616  double fm2 = m_func->eval(x - h2);
617  double fp3 = m_func->eval(x + h3);
618  double fm3 = m_func->eval(x - h3);
619  derivative = (5.0*(fp1-fm1) + 4.0*(fp2-fm2) + (fp3-fm3))/(32.0*h);
620  } // block needed against cross-initialization
621  break;
622  case 9:
623  {
624  double h2 = 2.0 * h;
625  double h3 = 3.0 * h;
626  double h4 = 4.0 * h;
627  double fp1 = m_func->eval(x + h);
628  double fm1 = m_func->eval(x - h);
629  double fp2 = m_func->eval(x + h2);
630  double fm2 = m_func->eval(x - h2);
631  double fp3 = m_func->eval(x + h3);
632  double fm3 = m_func->eval(x - h3);
633  double fp4 = m_func->eval(x + h4);
634  double fm4 = m_func->eval(x - h4);
635  derivative = (14.0*(fp1-fm1) + 14.0*(fp2-fm2) + 6.0*(fp3-fm3) +
636  (fp4-fm4))/(128.0*h);
637  } // block needed against cross-initialization
638  break;
639  case 11:
640  {
641  double h2 = 2.0 * h;
642  double h3 = 3.0 * h;
643  double h4 = 4.0 * h;
644  double h5 = 5.0 * h;
645  double fp1 = m_func->eval(x + h);
646  double fm1 = m_func->eval(x - h);
647  double fp2 = m_func->eval(x + h2);
648  double fm2 = m_func->eval(x - h2);
649  double fp3 = m_func->eval(x + h3);
650  double fm3 = m_func->eval(x - h3);
651  double fp4 = m_func->eval(x + h4);
652  double fm4 = m_func->eval(x - h4);
653  double fp5 = m_func->eval(x + h5);
654  double fm5 = m_func->eval(x - h5);
655  derivative = (42.0*(fp1-fm1) + 48.0*(fp2-fm2) + 27.0*(fp3-fm3) +
656  8.0*(fp4-fm4) + (fp5-fm5))/(512.0*h);
657  } // block needed against cross-initialization
658  break;
659  default:
660  std::string msg = "Specified filter length "+
661  gammalib::str(length)+" not implemented. Please "
662  "specify 5,7,9 or 11 as filter length.";
664  break;
665  }
666  break;
667 
668  // Degree 4
669  case 4:
670  switch (length) {
671  case 7:
672  {
673  double h2 = 2.0 * h;
674  double h3 = 3.0 * h;
675  double fp1 = m_func->eval(x + h);
676  double fm1 = m_func->eval(x - h);
677  double fp2 = m_func->eval(x + h2);
678  double fm2 = m_func->eval(x - h2);
679  double fp3 = m_func->eval(x + h3);
680  double fm3 = m_func->eval(x - h3);
681  derivative = (39.0*(fp1-fm1) + 12.0*(fp2-fm2) - 5.0*(fp3-fm3)) /
682  (96.0*h);
683  } // block needed against cross-initialization
684  break;
685  case 9:
686  {
687  double h2 = 2.0 * h;
688  double h3 = 3.0 * h;
689  double h4 = 4.0 * h;
690  double fp1 = m_func->eval(x + h);
691  double fm1 = m_func->eval(x - h);
692  double fp2 = m_func->eval(x + h2);
693  double fm2 = m_func->eval(x - h2);
694  double fp3 = m_func->eval(x + h3);
695  double fm3 = m_func->eval(x - h3);
696  double fp4 = m_func->eval(x + h4);
697  double fm4 = m_func->eval(x - h4);
698  derivative = (27.0*(fp1-fm1) + 16.0*(fp2-fm2) - (fp3-fm3) -
699  2.0*(fp4-fm4)) / (96.0*h);
700  } // block needed against cross-initialization
701  break;
702  case 11:
703  {
704  double h2 = 2.0 * h;
705  double h3 = 3.0 * h;
706  double h4 = 4.0 * h;
707  double h5 = 4.0 * h;
708  double fp1 = m_func->eval(x + h);
709  double fm1 = m_func->eval(x - h);
710  double fp2 = m_func->eval(x + h2);
711  double fm2 = m_func->eval(x - h2);
712  double fp3 = m_func->eval(x + h3);
713  double fm3 = m_func->eval(x - h3);
714  double fp4 = m_func->eval(x + h4);
715  double fm4 = m_func->eval(x - h4);
716  double fp5 = m_func->eval(x + h5);
717  double fm5 = m_func->eval(x - h5);
718  derivative = (322.0*(fp1-fm1) + 256.0*(fp2-fm2) + 39.0*(fp3-fm3) -
719  32.0*(fp4-fm4) - 11.0*(fp5-fm5)) / (1536.0*h);
720  } // block needed against cross-initialization
721  break;
722  default:
723  std::string msg = "Specified filter length "+
724  gammalib::str(length)+" not implemented. Please "
725  "specify 7,9 or 11 as filter length.";
727  break;
728  }
729  break;
730 
731  // Invalid degree
732  default:
733  std::string msg = "Specified degree "+gammalib::str(degree)+
734  " not implemented. Please specify 2 as degree.";
736  break;
737  }
738 
739  // Return derivative
740  return derivative;
741 }
742 
743 
744 /***********************************************************************//**
745  * @brief Print derivative information
746  *
747  * @param[in] chatter Chattiness (defaults to NORMAL).
748  * @return String containing derivative information.
749  ***************************************************************************/
750 std::string GDerivative::print(const GChatter& chatter) const
751 {
752  // Initialise result string
753  std::string result;
754 
755  // Continue only if chatter is not silent
756  if (chatter != SILENT) {
757 
758  // Append header
759  result.append("=== GDerivative ===");
760 
761  // Append information
762  result.append("\n"+gammalib::parformat("Relative precision"));
763  result.append(gammalib::str(eps()));
764  result.append("\n"+gammalib::parformat("Max. number of iterations"));
765  result.append(gammalib::str(max_iter()));
766  result.append("\n"+gammalib::parformat("Initial step fraction"));
767  result.append(gammalib::str(step_frac()));
768  if (silent()) {
769  result.append("\n"+gammalib::parformat("Warnings")+"suppressed");
770  }
771  else {
772  result.append("\n"+gammalib::parformat("Warnings"));
773  result.append("in standard output");
774  }
775 
776  } // endif: chatter was not silent
777 
778  // Return result
779  return result;
780 }
781 
782 
783 /*==========================================================================
784  = =
785  = Private methods =
786  = =
787  ==========================================================================*/
788 
789 /***********************************************************************//**
790  * @brief Initialise class members
791  ***************************************************************************/
793 {
794  // Initialise members
795  m_func = NULL;
796  m_eps = 1.0e-6;
797  m_step_frac = 0.02;
798  m_tiny = 1.0e-7;
799  m_max_iter = 5;
800  m_iter = 0;
801  m_silent = false;
802 
803  // Compute machine precision
804  //set_tiny();
805 
806  // Return
807  return;
808 }
809 
810 
811 /***********************************************************************//**
812  * @brief Copy class members
813  *
814  * @param[in] dx Derivative.
815  ***************************************************************************/
817 {
818  // Copy attributes
819  m_func = dx.m_func;
820  m_eps = dx.m_eps;
822  m_tiny = dx.m_tiny;
823  m_max_iter = dx.m_max_iter;
824  m_iter = dx.m_iter;
825  m_silent = dx.m_silent;
826 
827  // Return
828  return;
829 }
830 
831 
832 /***********************************************************************//**
833  * @brief Delete class members
834  ***************************************************************************/
836 {
837  // Return
838  return;
839 }
840 
841 
842 /***********************************************************************//**
843  * @brief Compute tiny number for Minuit2
844  *
845  * On fermi gives m_tiny = 8.88178e-16 for i=51.
846  ***************************************************************************/
848 {
849  // Allocate tiny instance
851 
852  // Calculate machine precision
853  double epstry = 0.5;
854  double epsbak = 0.0;
855  double epsp1 = 0.0;
856  double one = 1.0;
857 
858  // Loop until we found
859  for (int i = 0; i < 100; ++i) {
860  epstry *= 0.5;
861  epsp1 = one + epstry;
862  epsbak = tiny(epsp1);
863  if (epsbak < epstry) {
864  m_tiny = 8.0 * epstry;
865  break;
866  }
867  }
868 
869  // Return
870  return;
871 }
872 
873 
874 /***********************************************************************//**
875  * @brief Return one
876  ***************************************************************************/
877 double GDerivative::tiny::one(void) const
878 {
879  // Return
880  return m_one;
881 }
882 
883 
884 /***********************************************************************//**
885  * @brief Evaluate minimal difference between two floating points
886  ***************************************************************************/
887 double GDerivative::tiny::operator()(double eps) const
888 {
889  // Compute difference
890  double result = eps - one();
891 
892  // Return result
893  return result;
894 }
double right_difference(const double &x, const double &h)
Returns gradient computed from right-sided function difference.
double m_step_frac
Value fraction to use for initial step.
int m_max_iter
Maximum number of iterations.
double value(const double &x, const double &step=0.0)
Returns derivative.
const double & step_frac(void) const
Get step fraction.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1265
double one(void) const
Return one.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
Generic matrix class definition.
double operator()(double eps) const
Evaluate minimal difference between two floating points.
double m_eps
Derivative precision.
double difference(const double &x, const double &h)
Returns gradient computed from symmetric function difference.
Gammalib tools definition.
double ridder(const double &x, const double &h, double *err)
Returns derivative by Ridders&#39; method.
GDerivative * clone(void) const
Clone derivative.
GDerivative & operator=(const GDerivative &dx)
Assignment operator.
GFunction * m_func
Pointer to function.
const int & max_iter(void) const
Return maximum number of iterations.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
void init_members(void)
Initialise class members.
double m_tiny
Tiny number for minuit2.
const double & eps(void) const
Get precision.
double smooth_robust(const double &x, const double &h, const int &degree=2, const int &length=5)
Returns smooth noise-robust derivative.
#define G_VALUE
Definition: GDerivative.cpp:36
GChatter
Definition: GTypemaps.hpp:33
virtual double eval(const double &x)=0
const bool & silent(void) const
Get silence flag.
void set_tiny(void)
Compute tiny number for Minuit2.
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:872
GDerivative(void)
Void constructor.
Definition: GDerivative.cpp:60
std::string print(const GChatter &chatter=NORMAL) const
Print derivative information.
double left_difference(const double &x, const double &h)
Returns gradient computed from left-sided function difference.
void clear(void)
Clear derivative.
Numerical derivatives class.
Definition: GDerivative.hpp:45
GDerivative class interface definition.
Single parameter function abstract base class.
Definition: GFunction.hpp:44
virtual ~GDerivative(void)
Destructor.
#define G_RIDDER
Definition: GDerivative.cpp:37
bool m_silent
Suppress warnings.
void copy_members(const GDerivative &dx)
Copy class members.
Generic matrix class definition.
Definition: GMatrix.hpp:79
#define G_SMOOTH_ROBUST
Definition: GDerivative.cpp:39
#define G_MINUIT2
Definition: GDerivative.cpp:38
int m_iter
Number of iterations used.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void free_members(void)
Delete class members.
double minuit2(const double &x, double *err)
Returns derivative using Minuit2 algorithm.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413