Fix case when root is at the interval boundary

Fix also root finding in roostats in hypotestinverter that edges are returned when interval does not contain a root.
root-project · Feb 15, 2019 · 7aae501 · 7aae501
1 parent 6e00fcc
commit 7aae501
Show file tree

Hide file tree

Showing 3 changed files with 50 additions and 16 deletions.
diff --git a/math/mathcore/src/BrentMethods.cxx b/math/mathcore/src/BrentMethods.cxx
@@ -51,9 +51,16 @@ namespace BrentMethods {
       yymin = (*function)(xxmin);
    else if (type < 4)
       yymin = -(*function)(xxmin);
-   else
+   else { 
       yymin = (*function)(xxmin)-fy;
 
+      // case root is at the interval boundaries
+      if (yymin==0) {
+         xmin = xxmin; 
+         xmax = xxmin; 
+         return xxmin; 
+      }
+   }
    for (int i=1; i<=npx-1; i++) {
       double x = xmin + i*dx;
       if (logStep) x = std::exp(x);
@@ -71,8 +78,16 @@ namespace BrentMethods {
       }
       // when looking for root break at first instance
       if (type == 4 ) {
-         // found good interval if sign product is negative
-         if ( std::copysign(1.,y)*std::copysign(1.,yymin) < 0 ) {
+         // if root is at interval boundaries
+         if (y == 0) {
+            xmin = x;
+            xmax = x;
+            xmiddle = x;
+            foundInterval = true;
+            break; 
+         }
+         // found good interval if sign product is negative or equal zero to
+         if (std::copysign(1.,y)*std::copysign(1.,yymin) < 0 ) {
             xmin = xxmin; // previous value
             xmax = x;  // current value
             xmiddle = 0.5*(xmax+xmin);
@@ -81,7 +96,7 @@ namespace BrentMethods {
          }
          // continue bracketing
          xxmin = x;
-         yymin = y; 
+         yymin = y;
       }
    }
 
@@ -102,6 +117,8 @@ namespace BrentMethods {
       xmin = 1;
       xmax = 0; 
    }
+   // else 
+   //    std::cout << " root found !!! " << std::endl;
 
    return xmiddle;
 }

diff --git a/math/mathcore/test/stressTF1.cxx b/math/mathcore/test/stressTF1.cxx
@@ -10,7 +10,7 @@
 
 using namespace std;
 
-const double XMIN = 0, XMAX = 2*TMath::Pi();
+const double XMIN = 0, XMAX = 2*TMath::Pi() + 1.E-15;  // add an eps to avoid failing the search
 const Int_t NB = 100;
 const Int_t REP = 100000;
 const double TNORM = REP / 1000000.;
@@ -80,13 +80,15 @@ int TestRoot(TF1* f1)
    TStopwatch w;
    double totalTime = 0;
 
+   double eps = 0;
+
    cout << "ROOT TEST\n"
         << "---------------------------------------------------------"
         << endl;
 
    w.Start(kTRUE);
    for ( int j = 0; j < REP; ++j )
-      x = f1->GetX(0, XMIN, 1);
+      x = f1->GetX(0, XMIN-eps, 1);
    w.Stop();
    root = 0;
    status += PrintStatus("Root", x, root, w.RealTime()/TNORM );
@@ -102,7 +104,7 @@ int TestRoot(TF1* f1)
 
    w.Start(kTRUE);
    for ( int j = 0; j < REP; ++j )
-      x = f1->GetX(0, 6, XMAX);
+      x = f1->GetX(0, 6, XMAX+eps);
    w.Stop();
    root = TMath::Pi() * 2;
    status += PrintStatus("Root", x, root, w.RealTime()/ TNORM );

diff --git a/roofit/roostats/src/HypoTestInverterResult.cxx b/roofit/roostats/src/HypoTestInverterResult.cxx
@@ -610,6 +610,18 @@ double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool l
        varmax = var->getMax();
    }
 
+
+   // find ymin and ymax  and corresponding values
+   double ymin = TMath::MinElement(n,y);
+   double ymax = TMath::MaxElement(n,y);
+   // cannot find intercept in the full range - return min or max valie
+   if (ymax < y0) {
+      return (lowSearch) ? varmax : varmin; 
+   }
+   if (ymin > y0) {
+      return (lowSearch) ? varmin : varmax; 
+   }
+
    double xmin = axmin;
    double xmax = axmax;
 
@@ -626,13 +638,10 @@ double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool l
       double yfirst = graph.GetY()[0];
       double ylast = graph.GetY()[n-1];
 
-      // case we have lower/upper limits
-
-      // find ymin and ymax  and corresponding values
-      //double ymin = TMath::MinElement(n,y);
-      double ymax = TMath::MaxElement(n,y);
+      // distinguish the case we have lower /upper limits
+      // check if a possible crossing exists otherwise return variable min/max
+
       // do lower extrapolation
-
       if ( (ymax < y0 && !lowSearch) || ( yfirst > y0 && lowSearch) ) {
          xmin = varmin;
       }
@@ -650,9 +659,16 @@ double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool l
    ROOT::Math::BrentRootFinder brf;
    brf.SetFunction(f1d,xmin,xmax);
    brf.SetNpx(TMath::Max(graph.GetN()*2,100) );
+#ifdef DO_DEBUG
+   std::cout << "findind root for " << xmin << " ,  "<< xmax << "f(x) : " << graph.Eval(xmin) << " , " << graph.Eval(0.5*(xmax+xmin))
+             << " , " << graph.Eval(xmax) << " target " << y0 << std::endl;
+#endif
+
    bool ret = brf.Solve(100, 1.E-16, 1.E-6);
    if (!ret) {
-      ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed - return inf" << std::endl;
+      ooccoutE(this,Eval) << "HypoTestInverterResult - interpolation failed for interval [" << xmin << "," << xmax
+                          << " ]  g(xmin,xmax) =" << graph.Eval(xmin) << "," << graph.Eval(xmax)
+                          << " target=" << y0 << " return inf" << std::endl;
          return TMath::Infinity();
    }
    double limit =  brf.Root();
@@ -674,7 +690,6 @@ double HypoTestInverterResult::GetGraphX(const TGraph & graph, double y0, bool l
    auto file = TFile::Open(fname,"RECREATE");
    graph.Write("graph");
    file->Close();
-
 #endif
 
    // look in case if a new intersection exists
@@ -834,7 +849,7 @@ double HypoTestInverterResult::FindInterpolatedLimit(double target, bool lowSear
       << "the computed " << limitType << " limit is " << limit << " +/- " << error << std::endl;
 
 #ifdef DO_DEBUG
-   std::cout << "Found limit is " << limit << " +/- " << error << astd::endl;
+   std::cout << "Found limit is " << limit << " +/- " << error << std::endl;
 #endif