[RF] Avoid using out-of-range values in RooDerivative

Avoid using out-of-range values for the numerical differentiation in
RooDerivative.

This adds three new special cases:

  * derivative is first order and x-value is close to limit: just use
    one-directional derivative already implemented in the
    RichardsonDerivator

  * higher-order derivative but function is constant within floating
    point precision: just return zero

  * higher-order derivative for non-constant function: throw error,
    because we don't want to risk wrong derivatives by evaluating
    outside the range, which will result in x-value-clipping and
    therefore invalid y value differences

Closes #18615.

Author: guitargeek

roofit/roofitcore/src/RooDerivative.cxx

@@ -90,22 +90,60 @@ RooDerivative::~RooDerivative() = default;
 
 double RooDerivative::evaluate() const
 {
-  if (!_ftor) {
-    _ftor = std::unique_ptr<RooFunctor>{_func.arg().functor(_x.arg(),RooArgSet(),_nset)};
-    ROOT::Math::WrappedFunction<RooFunctor&> wf(*_ftor);
-    _rd = std::make_unique<ROOT::Math::RichardsonDerivator>(wf,_eps,true);
-  }
-
-  switch (_order) {
-  case 1: return _rd->Derivative1(_x);
-  case 2: return _rd->Derivative2(_x);
-  case 3: return _rd->Derivative3(_x);
-  }
-  return 0 ;
+   if (!_ftor) {
+      _ftor = std::unique_ptr<RooFunctor>{_func.arg().functor(_x.arg(), RooArgSet(), _nset)};
+      ROOT::Math::WrappedFunction<RooFunctor &> wf(*_ftor);
+      _rd = std::make_unique<ROOT::Math::RichardsonDerivator>(wf, _eps, true);
+   }
+
+   // Figure out if we are close to the variable boundaries
+   double val = _x;
+   auto &xVar = static_cast<RooRealVar &>(*_x);
+   double valMin = xVar.getMin();
+   double valMax = xVar.getMax();
+   bool isCloseLo = val - valMin < _eps;
+   bool isCloseHi = valMax - val < _eps;
+
+   // If we hit the boundary left and right, there is obviously a mistake when setting epsilon
+   if (isCloseLo && isCloseHi) {
+      std::stringstream errMsg;
+      errMsg << "error in numerical derivator: 2 * epsilon is larger than the variable range!";
+      coutE(Eval) << errMsg.str() << std::endl;
+      throw std::runtime_error(errMsg.str());
+   }
+
+   // If we are close to the variable boundary on either side
+   if (isCloseLo || isCloseHi) {
+      // For first-order derivatives we are good: the RichardsonDerivator can
+      // also calculate the derivative using only forward or backward
+      // variations:
+      if (_order == 1) {
+         return isCloseLo ? _rd->DerivativeForward(val) : _rd->DerivativeBackward(val);
+      }
+      // If the function is constant within floating point precision anyway,
+      // we don't have a problem.
+      const double eps = std::numeric_limits<double>::epsilon();
+      const double yval1 = _ftor->eval(val);
+      const double yval2 = isCloseLo ? _ftor->eval(val + _eps) : _ftor->eval(val - _eps);
+      if (std::abs(yval2 - yval1) <= eps) {
+         return 0.0;
+      }
+
+      // Give up
+      std::stringstream errMsg;
+      errMsg << "error in numerical derivator: variable value is to close to limits to compute finite differences";
+      coutE(Eval) << errMsg.str() << std::endl;
+      throw std::runtime_error(errMsg.str());
+   }
+
+   switch (_order) {
+   case 1: return _rd->Derivative1(val);
+   case 2: return _rd->Derivative2(val);
+   case 3: return _rd->Derivative3(val);
+   }
+   return 0;
 }
 
-
-
 ////////////////////////////////////////////////////////////////////////////////
 /// Zap functor and derivator ;
 

tutorials/roofit/roofit/rf111_derivatives.C

@@ -35,12 +35,6 @@ void rf111_derivatives()
    RooRealVar mean("mean", "mean of gaussian", 1, -10, 10);
    RooRealVar sigma("sigma", "width of gaussian", 1, 0.1, 10);
 
-   // Ranges for plotting derivatives. They can't be the full range, because
-   // the derivatives is calculated by variation around the central point. We
-   // need to ensure that the varied values are still in the full range.
-   x.setRange("plotrange", -9.95, 9.95);
-   sigma.setRange("plotrange", 0.15, 1.95);
-
    // Build gaussian pdf in terms of x,mean and sigma
    RooGaussian gauss("gauss", "gaussian PDF", x, mean, sigma);
 
@@ -61,9 +55,9 @@ void rf111_derivatives()
    gauss.plotOn(xframe);
 
    // Plot derivatives in same frame
-   dgdx->plotOn(xframe, LineColor(kMagenta), Range("plotrange"));
-   d2gdx2->plotOn(xframe, LineColor(kRed), Range("plotrange"));
-   d3gdx3->plotOn(xframe, LineColor(kOrange), Range("plotrange"));
+   dgdx->plotOn(xframe, LineColor(kMagenta));
+   d2gdx2->plotOn(xframe, LineColor(kRed));
+   d3gdx3->plotOn(xframe, LineColor(kOrange));
 
    // C r e a t e   a n d   p l o t  d e r i v a t i v e s   w . r . t .   s i g m a
    // ------------------------------------------------------------------------------
@@ -82,9 +76,9 @@ void rf111_derivatives()
    gauss.plotOn(sframe);
 
    // Plot derivatives in same frame
-   dgds->plotOn(sframe, LineColor(kMagenta), Range("plotrange"));
-   d2gds2->plotOn(sframe, LineColor(kRed), Range("plotrange"));
-   d3gds3->plotOn(sframe, LineColor(kOrange), Range("plotrange"));
+   dgds->plotOn(sframe, LineColor(kMagenta));
+   d2gds2->plotOn(sframe, LineColor(kRed));
+   d3gds3->plotOn(sframe, LineColor(kOrange));
 
    // Draw all frames on a canvas
    TCanvas *c = new TCanvas("rf111_derivatives", "rf111_derivatives", 800, 400);

tutorials/roofit/roofit/rf111_derivatives.py

@@ -25,12 +25,6 @@
 mean = ROOT.RooRealVar("mean", "mean of gaussian", 1, -10, 10)
 sigma = ROOT.RooRealVar("sigma", "width of gaussian", 1, 0.1, 10)
 
-# Ranges for plotting derivatives. They can't be the full range, because
-# the derivatives is calculated by variation around the central point. We
-# need to ensure that the varied values are still in the full range.
-x.setRange("plotrange", -9.95, 9.95)
-sigma.setRange("plotrange", 0.15, 1.95)
-
 # Build gaussian pdf in terms of x, and sigma
 gauss = ROOT.RooGaussian("gauss", "gaussian PDF", x, mean, sigma)
 
@@ -51,9 +45,9 @@
 gauss.plotOn(xframe)
 
 # Plot derivatives in same frame
-dgdx.plotOn(xframe, LineColor="m", Range="plotrange")
-d2gdx2.plotOn(xframe, LineColor="r", Range="plotrange")
-d3gdx3.plotOn(xframe, LineColor="kOrange", Range="plotrange")
+dgdx.plotOn(xframe, LineColor="m")
+d2gdx2.plotOn(xframe, LineColor="r")
+d3gdx3.plotOn(xframe, LineColor="kOrange")
 
 # Create and plot derivatives w.r.t. sigma
 # ------------------------------------------------------------------------------
@@ -72,9 +66,9 @@
 gauss.plotOn(sframe)
 
 # Plot derivatives in same frame
-dgds.plotOn(sframe, LineColor="m", Range="plotrange")
-d2gds2.plotOn(sframe, LineColor="r", Range="plotrange")
-d3gds3.plotOn(sframe, LineColor="kOrange", Range="plotrange")
+dgds.plotOn(sframe, LineColor="m")
+d2gds2.plotOn(sframe, LineColor="r")
+d3gds3.plotOn(sframe, LineColor="kOrange")
 
 # Draw all frames on a canvas
 c = ROOT.TCanvas("rf111_derivatives", "rf111_derivatives", 800, 400)