Added additional check on slope.

MeasureTransport · Oct 17, 2023 · 892c746 · 892c746
1 parent 23a7604
commit 892c746
Show file tree

Hide file tree

Showing 2 changed files with 8 additions and 10 deletions.
diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
@@ -43,9 +43,9 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
             xlb = xub-stepSize;
             ylb = f(xlb);
 
-            if(abs((yub-ylb)/(xub-xlb))<1e-12){
+            if((fabs((yub-ylb)/(xub-xlb))<1e-12)&&((xub-xlb)>10)){
                 info = -1;
-                break;
+                return;
             }
 
             if(ylb>yd){
@@ -67,11 +67,11 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
             xub = xlb+stepSize;
             yub = f(xub);
-          
-            // Check to see if function is perfectly flat
-            if(abs((yub-ylb)/(xub-xlb))<1e-12){
+
+            // Check to see if function is perfectly flat over a wide region
+            if((fabs((yub-ylb)/(xub-xlb))<1e-12)&&((xub-xlb)>10)){
                 info = -1;
-                break;
+                return;
             }
 
             if(yub<yd){
@@ -128,6 +128,7 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     // Compute initial bracket containing the root
     int bracket_info = 0;
     FindBracket<MemorySpace>(f, xlb, ylb, xub, yub, yd, bracket_info);
+
     if((bracket_info<0)||(((ylb>yd)||(yub<yd)))){
         info = -2;
         return std::numeric_limits<double>::quiet_NaN();
@@ -142,7 +143,6 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     unsigned int it;
     for(it=0; it<maxIts; ++it){
 
-        //std::cout << "Iteration " << it << std::endl;
         xc = Find_x_ITP(xlb, xub, yd, ylb, yub, k1, k2, nhalf, n0, it, xtol);
 
         yc = f(xc);

diff --git a/tests/Test_RootFinding.cpp b/tests/Test_RootFinding.cpp
@@ -24,7 +24,6 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
     auto sigmoid = [](double x){return 1./(std::exp(-x)+1);};
     auto sigmoid_combo = [sigmoid](double x){return sigmoid(2*(x-1)) + 3*sigmoid(x*0.5) + 0.5*sigmoid(1.5*(x+1));};
 
-
     SECTION("FindBracket linear") {
         double xd = -1.1;
         double yd = identity(xd);
@@ -46,7 +45,6 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
         CHECK(info==0);
     }
     SECTION("FindBracket flat") {
-        double xd = -1.1;
         double yd = 0.0;
         double xub = 2., yub = -1.0;
         double xlb = 0., ylb = -1.0;
@@ -72,7 +70,7 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
         CHECK(info==0);
     }
     SECTION("Test Inverse Flat") {
-        double xd = 0.5, yd = -1.0;
+        double xd = 0.5, yd = 0.0;
         double x0 = 0.0, xtol = 1e-5, ftol = 1e-5;
         auto f = [](double x){return -1.0;};
         int info = 0;