Johan Bregeon / gammalib

Commit 73ea5ed7381ffded95c0d4c01b4f066ae1d39dad

Authored by Jürgen Knödlseder
1 parent 6bae1947

Add smooth_robust() method to compute smooth noise-robust derivatives.

Inline Side-by-side
Showing 3 changed files with 191 additions and 9 deletions
include/GDerivative.hpp
  View file @73ea5ed
1 1 /***************************************************************************
2 2 * GDerivative.hpp - Derivative class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2014 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -73,6 +73,9 @@ public:
73 73 double ridder(const double& x, const double& h, double* err);
74 74 double minuit2(const double& x, double* err);
75 75 double difference(const double& x, const double& h);
  76 + double smooth_robust(const double& x, const double& h,
  77 + const int& degree = 2,
  78 + const int& length = 5);
76 79 std::string print(const GChatter& chatter = NORMAL) const;
77 80  
78 81 protected:
... ...
pyext/GDerivative.i
  View file @73ea5ed
1 1 /***************************************************************************
2 2 * GDerivative.i - Derivative class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2014 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -60,6 +60,9 @@ public:
60 60 double value(const double& x, const double& step = 0.0);
61 61 double ridder(const double& x, const double& h, double* err);
62 62 double minuit2(const double& x, double* err);
  63 + double smooth_robust(const double& x, const double& h,
  64 + const int& degree = 2,
  65 + const int& length = 5);
63 66 double difference(const double& x, const double& h);
64 67 };
65 68  
... ...
src/numerics/GDerivative.cpp
  View file @73ea5ed
1 1 /***************************************************************************
2 2 * GDerivative.cpp - Derivative class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2011-2013 by Juergen Knoedlseder *
  4 + * copyright (C) 2011-2014 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -36,6 +36,8 @@
36 36 #define G_VALUE "GDerivative::value(double&, double&)"
37 37 #define G_RIDDER "GDerivative::ridder(double&, double&, double*)"
38 38 #define G_MINUIT2 "GDerivative::minuit2(double&, double*)"
  39 +#define G_SMOOTH_ROBUST "GDerivative::smooth_robust(double&, double&, int&,"\
  40 + " int&)"
39 41  
40 42 /* __ Macros _____________________________________________________________ */
41 43  
... ... @@ -511,13 +513,187 @@ double GDerivative::minuit2(const double& x, double* err)
511 513 ***************************************************************************/
512 514 double GDerivative::difference(const double& x, const double& h)
513 515 {
514   - // Compute gradient
515   - double fs1 = m_func->eval(x + h);
516   - double fs2 = m_func->eval(x - h);
517   - double grd = 0.5 * (fs1 - fs2) / h;
  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;
518 520  
519   - // Return gradient
520   - return grd;
  521 + // Return derivative
  522 + return derivative;
  523 +}
  524 +
  525 +
  526 +/***********************************************************************//**
  527 + * @brief Returns smooth noise-robust derivative
  528 + *
  529 + * @param[in] x Function value.
  530 + * @param[in] h Step size.
  531 + * @param[in] degree Differentiator degree (default: 2).
  532 + * @param[in] length Filter length (default: 5).
  533 + *
  534 + * Compute a smooth noise-robust numerical derivative using the method
  535 + * proposed by Pavel Holoborodko
  536 + * (http://www.holoborodko.com/pavel/numerical-methods/)
  537 + * Formulae for the following filter lengths and degrees have been
  538 + * implemented:
  539 + *
  540 + * degree=2 length=5,7,9,11
  541 + * degree=4 length=7,9,11
  542 + ***************************************************************************/
  543 +double GDerivative::smooth_robust(const double& x, const double& h,
  544 + const int& degree, const int& length)
  545 +{
  546 + // Initialise derivative
  547 + double derivative = 0.0;
  548 +
  549 + // Use method dependent on length and degree
  550 + switch (degree) {
  551 +
  552 + // Degree 2
  553 + case 2:
  554 + switch (length) {
  555 + case 5:
  556 + {
  557 + double h2 = 2.0 * h;
  558 + double fp1 = m_func->eval(x + h);
  559 + double fm1 = m_func->eval(x - h);
  560 + double fp2 = m_func->eval(x + h2);
  561 + double fm2 = m_func->eval(x - h2);
  562 + derivative = (2.0*(fp1-fm1) + (fp2-fm2))/(8.0*h);
  563 + } // block needed against cross-initialization
  564 + break;
  565 + case 7:
  566 + {
  567 + double h2 = 2.0 * h;
  568 + double h3 = 3.0 * h;
  569 + double fp1 = m_func->eval(x + h);
  570 + double fm1 = m_func->eval(x - h);
  571 + double fp2 = m_func->eval(x + h2);
  572 + double fm2 = m_func->eval(x - h2);
  573 + double fp3 = m_func->eval(x + h3);
  574 + double fm3 = m_func->eval(x - h3);
  575 + derivative = (5.0*(fp1-fm1) + 4.0*(fp2-fm2) + (fp3-fm3))/(32.0*h);
  576 + } // block needed against cross-initialization
  577 + break;
  578 + case 9:
  579 + {
  580 + double h2 = 2.0 * h;
  581 + double h3 = 3.0 * h;
  582 + double h4 = 4.0 * h;
  583 + double fp1 = m_func->eval(x + h);
  584 + double fm1 = m_func->eval(x - h);
  585 + double fp2 = m_func->eval(x + h2);
  586 + double fm2 = m_func->eval(x - h2);
  587 + double fp3 = m_func->eval(x + h3);
  588 + double fm3 = m_func->eval(x - h3);
  589 + double fp4 = m_func->eval(x + h4);
  590 + double fm4 = m_func->eval(x - h4);
  591 + derivative = (14.0*(fp1-fm1) + 14.0*(fp2-fm2) + 6.0*(fp3-fm3) +
  592 + (fp4-fm4))/(128.0*h);
  593 + } // block needed against cross-initialization
  594 + break;
  595 + case 11:
  596 + {
  597 + double h2 = 2.0 * h;
  598 + double h3 = 3.0 * h;
  599 + double h4 = 4.0 * h;
  600 + double h5 = 5.0 * h;
  601 + double fp1 = m_func->eval(x + h);
  602 + double fm1 = m_func->eval(x - h);
  603 + double fp2 = m_func->eval(x + h2);
  604 + double fm2 = m_func->eval(x - h2);
  605 + double fp3 = m_func->eval(x + h3);
  606 + double fm3 = m_func->eval(x - h3);
  607 + double fp4 = m_func->eval(x + h4);
  608 + double fm4 = m_func->eval(x - h4);
  609 + double fp5 = m_func->eval(x + h5);
  610 + double fm5 = m_func->eval(x - h5);
  611 + derivative = (42.0*(fp1-fm1) + 48.0*(fp2-fm2) + 27.0*(fp3-fm3) +
  612 + 8.0*(fp4-fm4) + (fp5-fm5))/(512.0*h);
  613 + } // block needed against cross-initialization
  614 + break;
  615 + default:
  616 + std::string msg = "Specified filter length "+
  617 + gammalib::str(length)+" not implemented. Please "
  618 + "specify 5,7,9 or 11 as filter length.";
  619 + throw GException::invalid_argument(G_SMOOTH_ROBUST, msg);
  620 + break;
  621 + }
  622 + break;
  623 +
  624 + // Degree 4
  625 + case 4:
  626 + switch (length) {
  627 + case 7:
  628 + {
  629 + double h2 = 2.0 * h;
  630 + double h3 = 3.0 * h;
  631 + double fp1 = m_func->eval(x + h);
  632 + double fm1 = m_func->eval(x - h);
  633 + double fp2 = m_func->eval(x + h2);
  634 + double fm2 = m_func->eval(x - h2);
  635 + double fp3 = m_func->eval(x + h3);
  636 + double fm3 = m_func->eval(x - h3);
  637 + derivative = (39.0*(fp1-fm1) + 12.0*(fp2-fm2) - 5.0*(fp3-fm3)) /
  638 + (96.0*h);
  639 + } // block needed against cross-initialization
  640 + break;
  641 + case 9:
  642 + {
  643 + double h2 = 2.0 * h;
  644 + double h3 = 3.0 * h;
  645 + double h4 = 4.0 * h;
  646 + double fp1 = m_func->eval(x + h);
  647 + double fm1 = m_func->eval(x - h);
  648 + double fp2 = m_func->eval(x + h2);
  649 + double fm2 = m_func->eval(x - h2);
  650 + double fp3 = m_func->eval(x + h3);
  651 + double fm3 = m_func->eval(x - h3);
  652 + double fp4 = m_func->eval(x + h4);
  653 + double fm4 = m_func->eval(x - h4);
  654 + derivative = (27.0*(fp1-fm1) + 16.0*(fp2-fm2) - (fp3-fm3) -
  655 + 2.0*(fp4-fm4)) / (96.0*h);
  656 + } // block needed against cross-initialization
  657 + break;
  658 + case 11:
  659 + {
  660 + double h2 = 2.0 * h;
  661 + double h3 = 3.0 * h;
  662 + double h4 = 4.0 * h;
  663 + double h5 = 4.0 * h;
  664 + double fp1 = m_func->eval(x + h);
  665 + double fm1 = m_func->eval(x - h);
  666 + double fp2 = m_func->eval(x + h2);
  667 + double fm2 = m_func->eval(x - h2);
  668 + double fp3 = m_func->eval(x + h3);
  669 + double fm3 = m_func->eval(x - h3);
  670 + double fp4 = m_func->eval(x + h4);
  671 + double fm4 = m_func->eval(x - h4);
  672 + double fp5 = m_func->eval(x + h5);
  673 + double fm5 = m_func->eval(x - h5);
  674 + derivative = (322.0*(fp1-fm1) + 256.0*(fp2-fm2) + 39.0*(fp3-fm3) -
  675 + 32.0*(fp4-fm4) - 11.0*(fp5-fm5)) / (1536.0*h);
  676 + } // block needed against cross-initialization
  677 + break;
  678 + default:
  679 + std::string msg = "Specified filter length "+
  680 + gammalib::str(length)+" not implemented. Please "
  681 + "specify 7,9 or 11 as filter length.";
  682 + throw GException::invalid_argument(G_SMOOTH_ROBUST, msg);
  683 + break;
  684 + }
  685 + break;
  686 +
  687 + // Invalid degree
  688 + default:
  689 + std::string msg = "Specified degree "+gammalib::str(degree)+
  690 + " not implemented. Please specify 2 as degree.";
  691 + throw GException::invalid_argument(G_SMOOTH_ROBUST, msg);
  692 + break;
  693 + }
  694 +
  695 + // Return derivative
  696 + return derivative;
521 697 }
522 698  
523 699  
... ...