From cb001a7c6947cb2bfa17895dbae5f71b96b78758 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Fri, 13 Oct 2023 13:53:32 -0400
Subject: [PATCH 1/6] Fixed bug in computation of initial bracket for inverse
 evaluation.

---
 MParT/Utilities/RootFinding.h | 82 ++++++++++++++++++++++++++---------
 tests/Test_RootFinding.cpp    | 26 +++--------
 2 files changed, 67 insertions(+), 41 deletions(-)

diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
index 6d8702ea..c290df4a 100644
--- a/MParT/Utilities/RootFinding.h
+++ b/MParT/Utilities/RootFinding.h
@@ -11,22 +11,62 @@ KOKKOS_INLINE_FUNCTION void swapPair(T& x1, T& x2, T& y1, T& y2) {
     simple_swap(y1, y2);
 }
 
+/** Finds a bracket [xlb, xub] such that f(xlb)<yd and f(xub)>yd. */
 template<typename MemorySpace, typename FunctorType>
-KOKKOS_INLINE_FUNCTION void FindBound(bool haveLowerBound, double yd,
-    FunctorType f, double& x_bound, double& y_bound,
-    double& x_unbound, double& y_unbound, const unsigned int maxIts) {
-    double boundSign = haveLowerBound ? 1.0 : -1.0;
+KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
+                                        double& xlb, double& ylb,
+                                        double& xub, double& yub,
+                                        const double yd)
+{
+    double xb, xf; // Bisection point and regula falsi point
+    double xc, yc;
+    const unsigned int maxIts = 1000;
     double stepSize = 1.0;
-    unsigned int iter;
-    for(iter = 0; iter < maxIts; iter++) {
-        x_unbound += boundSign*stepSize;
-        y_unbound = f(x_unbound);
-        if(boundSign*(y_unbound - yd) > 0) return;
-        swapPair(x_bound, x_unbound, y_bound, y_unbound);
-        stepSize *= 2;
+
+    ylb = f(xb);
+
+    // We actually found an upper bound...
+    if(ylb>yd){
+
+        mpart::simple_swap(ylb,yub);
+        mpart::simple_swap(xlb,xub);
+
+        // Now find a lower bound...
+        unsigned int i;
+        for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
+            xlb = xub-stepSize;
+            ylb = f(xlb);
+            if(ylb>yd){
+                mpart::simple_swap(ylb,yub);
+                mpart::simple_swap(xlb,xub);
+                stepSize *= 2.0;
+            }else{
+                break;
+            }
+        }
+
+        if(i>=maxIts)
+            ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
+
+    // We have a lower bound...
+    }else{
+        // Now find an upper bound...
+        unsigned int i;
+        for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
+            xub = xlb+stepSize;
+            yub = f(xub);
+            if(yub<yd){
+                mpart::simple_swap(ylb,yub);
+                mpart::simple_swap(xlb,xub);
+                stepSize *= 2.0;
+            }else{
+                break;
+            }
+        }
+
+        if(i>=maxIts)
+            ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
     }
-    if(iter>=maxIts)
-            ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: bound calculation exceeds maxIts");
 }
 
 KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, double ylb, double yub,
@@ -47,7 +87,8 @@ KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, doub
 
 template<typename MemorySpace, typename FunctorType>
 KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, double x0, const double xtol, const double ftol)
-{
+{   
+    std::cout << "Hereherehere" << std::endl;
     double stepSize=1.0;
     const unsigned int maxIts = 10000;
 
@@ -59,13 +100,10 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     xlb = xub = x0;
     ylb = yub = f(xlb);
 
-    bool haveLowerBound = ylb < yd; // if ylb is a lower bound for yd
-    if(haveLowerBound) {
-        FindBound<MemorySpace>(haveLowerBound, yd, f, xlb, ylb, xub, yub, maxIts);
-    } else {
-        FindBound<MemorySpace>(haveLowerBound, yd, f, xub, yub, xlb, ylb, maxIts);
-    }
-
+    // Compute bounds
+    std::cout << "About to call find bracket..." << std::endl;
+    FindBracket<MemorySpace>(f, xlb, ylb, xub, yub, yd);
+    std::cout << "done " << xlb << ", " << xub << ", " << ylb << ", " << yub << std::endl;
     assert(ylb<yub);
     assert(xlb<xub);
 
@@ -77,7 +115,9 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
 
     unsigned int it;
     for(it=0; it<maxIts; ++it){
+        std::cout << " about to call Find_x_ITP" << std::endl;
         xc = Find_x_ITP(xlb, xub, yd, ylb, yub, k1, k2, nhalf, n0, it, xtol);
+        std::cout << " done" << std::endl;
 
         yc = f(xc);
 
diff --git a/tests/Test_RootFinding.cpp b/tests/Test_RootFinding.cpp
index 5bfd4b0c..f098597c 100644
--- a/tests/Test_RootFinding.cpp
+++ b/tests/Test_RootFinding.cpp
@@ -23,36 +23,22 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
     auto identity = [](double x){return x;};
     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("FindBound lower linear") {
+    
+    
+    SECTION("FindBracket linear") {
         double xd = -1.1;
         double yd = identity(xd);
         double xub = 2., yub = 2.;
         double xlb = 0., ylb = 0.;
-        FindBound<HostSpace>(false, yd, identity, xub, yub, xlb, ylb, 10'000);
+        FindBracket<HostSpace>(identity, xlb, ylb, xub, ylb, 10'000);
         CheckFoundBounds(identity, xlb, xd, xub, ylb, yd, yub);
     }
-    SECTION("FindBound upper linear") {
-        double xd = 1.1;
-        double yd = identity(xd);
-        double xub = 0., yub = 0.;
-        double xlb = -2., ylb = -2.;
-        FindBound<HostSpace>(true, yd, identity, xlb, ylb, xub, yub, 10'000);
-        CheckFoundBounds(identity, xlb, xd, xub, ylb, yd, yub);
-    }
-    SECTION("FindBound lower sigmoid") {
+    SECTION("FindBracker sigmoid") {
         double xd = -0.5;
         double yd = sigmoid(xd);
         double xub =  2., yub = sigmoid(xub);
         double xlb =  0., ylb = sigmoid(xlb);
-        FindBound<HostSpace>(false, yd, sigmoid, xub, yub, xlb, ylb, 10'000);
-        CheckFoundBounds(sigmoid, xlb, xd, xub, ylb, yd, yub);
-    }
-    SECTION("FindBound upper sigmoid") {
-        double xd = 0.5;
-        double yd = sigmoid(xd);
-        double xub =  0., yub = sigmoid(xub);
-        double xlb = -2., ylb = sigmoid(xlb);
-        FindBound<HostSpace>(true, yd, sigmoid, xlb, ylb, xub, yub, 10'000);
+        FindBracket<HostSpace>(sigmoid, xlb, ylb, xub, ylb, 10'000);
         CheckFoundBounds(sigmoid, xlb, xd, xub, ylb, yd, yub);
     }
     SECTION("Test Inverse Linear, low x0") {

From 6b6adc1afb25be4b432f4626feaffc420a93c836 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Fri, 13 Oct 2023 14:29:19 -0400
Subject: [PATCH 2/6] More bug fixes related to root finding tests.

---
 MParT/MonotoneComponent.h     |  2 +-
 MParT/Utilities/RootFinding.h | 18 ++++++------------
 tests/Test_RootFinding.cpp    |  4 ++--
 3 files changed, 9 insertions(+), 15 deletions(-)

diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h
index 30b9a9d4..ed9c0c80 100644
--- a/MParT/MonotoneComponent.h
+++ b/MParT/MonotoneComponent.h
@@ -400,7 +400,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
 
                 // Compute the inverse
                 Kokkos::View<double*,MemorySpace> workspace(team_member.thread_scratch(1), workspaceSize);
-                auto eval = SingleEvaluator(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_);
+                auto eval = SingleEvaluator<decltype(pt),decltype(coeffs)>(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_);
                 output(ptInd) = RootFinding::InverseSingleBracket<MemorySpace>(ys(ptInd), eval, pt(pt.extent(0)-1), xtol, ytol);
             }
         };
diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
index c290df4a..5228174c 100644
--- a/MParT/Utilities/RootFinding.h
+++ b/MParT/Utilities/RootFinding.h
@@ -18,16 +18,15 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
                                         double& xub, double& yub,
                                         const double yd)
 {
-    double xb, xf; // Bisection point and regula falsi point
-    double xc, yc;
     const unsigned int maxIts = 1000;
     double stepSize = 1.0;
-
-    ylb = f(xb);
-
+    
+    ylb = f(xlb);
+    yub = f(xub);
+    
     // We actually found an upper bound...
     if(ylb>yd){
-
+    
         mpart::simple_swap(ylb,yub);
         mpart::simple_swap(xlb,xub);
 
@@ -67,7 +66,7 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
         if(i>=maxIts)
             ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
     }
-}
+ }
 
 KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, double ylb, double yub,
                            double k1, double k2, double nhalf, double n0, int it, double xtol) {
@@ -88,7 +87,6 @@ KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, doub
 template<typename MemorySpace, typename FunctorType>
 KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, double x0, const double xtol, const double ftol)
 {   
-    std::cout << "Hereherehere" << std::endl;
     double stepSize=1.0;
     const unsigned int maxIts = 10000;
 
@@ -101,9 +99,7 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     ylb = yub = f(xlb);
 
     // Compute bounds
-    std::cout << "About to call find bracket..." << std::endl;
     FindBracket<MemorySpace>(f, xlb, ylb, xub, yub, yd);
-    std::cout << "done " << xlb << ", " << xub << ", " << ylb << ", " << yub << std::endl;
     assert(ylb<yub);
     assert(xlb<xub);
 
@@ -115,9 +111,7 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
 
     unsigned int it;
     for(it=0; it<maxIts; ++it){
-        std::cout << " about to call Find_x_ITP" << std::endl;
         xc = Find_x_ITP(xlb, xub, yd, ylb, yub, k1, k2, nhalf, n0, it, xtol);
-        std::cout << " done" << std::endl;
 
         yc = f(xc);
 
diff --git a/tests/Test_RootFinding.cpp b/tests/Test_RootFinding.cpp
index f098597c..b630ee35 100644
--- a/tests/Test_RootFinding.cpp
+++ b/tests/Test_RootFinding.cpp
@@ -30,7 +30,7 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
         double yd = identity(xd);
         double xub = 2., yub = 2.;
         double xlb = 0., ylb = 0.;
-        FindBracket<HostSpace>(identity, xlb, ylb, xub, ylb, 10'000);
+        FindBracket<HostSpace>(identity, xlb, ylb, xub, yub, yd);
         CheckFoundBounds(identity, xlb, xd, xub, ylb, yd, yub);
     }
     SECTION("FindBracker sigmoid") {
@@ -38,7 +38,7 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
         double yd = sigmoid(xd);
         double xub =  2., yub = sigmoid(xub);
         double xlb =  0., ylb = sigmoid(xlb);
-        FindBracket<HostSpace>(sigmoid, xlb, ylb, xub, ylb, 10'000);
+        FindBracket<HostSpace>(sigmoid, xlb, ylb, xub, yub, yd);
         CheckFoundBounds(sigmoid, xlb, xd, xub, ylb, yd, yub);
     }
     SECTION("Test Inverse Linear, low x0") {

From e4eeb5d6562da6f20e0091123f924e4849793174 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Fri, 13 Oct 2023 21:42:27 -0400
Subject: [PATCH 3/6] Now throwing nan if inverse blows up.

---
 MParT/MonotoneComponent.h       | 270 +++++++++++++++++++-------------
 MParT/Utilities/Miscellaneous.h |   2 +-
 MParT/Utilities/RootFinding.h   |  50 ++++--
 3 files changed, 198 insertions(+), 124 deletions(-)

diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h
index ed9c0c80..6b6be68f 100644
--- a/MParT/MonotoneComponent.h
+++ b/MParT/MonotoneComponent.h
@@ -312,6 +312,16 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                      StridedVector<double, MemorySpace>              output,
                  std::map<std::string, std::string>                  options=std::map<std::string,std::string>())
     {
+        // std::cout << "xs.shape = " << xs.extent(0) << "," << xs.extent(1) << std::endl;
+        // std::cout << "xs = \n" << std::endl;
+        // for(unsigned int row=0; row<xs.extent(0); ++row){
+        //     for(unsigned int col=0; col<xs.extent(1); ++col){
+        //         std::cout << "  " << xs(row,col);
+        //     }
+        //     std::cout << std::endl;
+        // }
+        // std::cout << std::endl;
+
         // Extract the method from the options map
         std::string method;
         if(options.count("Method")){
@@ -393,15 +403,39 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
 
                 // Create a subview containing x_{1:d-1}
                 auto pt = Kokkos::subview(xs, Kokkos::ALL(), xInd);
+                
+                //std::cout << "\n\nPt0 = \n" << std::endl;
+                for(unsigned int ii=0; ii<pt.size(); ++ii){
+                    //std::cout << "  " << pt(ii);
+                    if(std::isnan(pt(ii))){
+                        std::cout << "Warning: Found nan in point, not computing inverse!" << std::endl;
+                        //throw std::runtime_error("Found nan in point!");
+                        return;
+                    }
+                }
+                //std::cout << std::endl << std::endl;
 
                 // Fill in the cache with everything that doesn't depend on x_d
                 Kokkos::View<double*,MemorySpace> cache(team_member.thread_scratch(1), cacheSize);
                 expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None);
 
+                // std::cout << "\n\nPt1 = \n" << std::endl;
+                // double max_val = 0.0;
+                // for(unsigned int ii=0; ii<pt.size(); ++ii){
+                //     std::cout << "  " << pt(ii);
+                //     max_val  = (abs(pt(ii))>max_val) ?  abs(pt(ii)) : max_val;
+                // }
+                // std::cout << std::endl << std::endl;
+                // if (max_val > 1e-4){
+                //     std::cout << "Max val = " << max_val;
+                //     assert(max_val < 1e4);
+                // }
                 // Compute the inverse
                 Kokkos::View<double*,MemorySpace> workspace(team_member.thread_scratch(1), workspaceSize);
                 auto eval = SingleEvaluator<decltype(pt),decltype(coeffs)>(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_);
+                //std::cout << "Here 0" << std::endl;
                 output(ptInd) = RootFinding::InverseSingleBracket<MemorySpace>(ys(ptInd), eval, pt(pt.extent(0)-1), xtol, ytol);
+                //std::cout << "Here 1" << std::endl;
             }
         };
 
@@ -967,6 +1001,16 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                                                  QuadratureType    const& quad,
                                                  ExpansionType     const& expansion)
     {
+        // std::cout << "pt = " << std::endl;
+        // for(unsigned int i=0; i<pt.size(); ++i)
+        //     std::cout << "  " << pt(i);
+        // std::cout << std::endl;
+        // std::cout << "xd = " << xd << std::endl;
+        // std::cout << "Coefficients = " << std::endl;
+        // for(unsigned int i=0; i<coeffs.size(); ++i)
+        //     std::cout << "  " << coeffs(i);
+        // std::cout << std::endl;
+        
         double output = 0.0;
         // Compute the integral \int_0^1 g( \partial_D f(x_1,...,x_{D-1},t*x_d)) dt
         MonotoneIntegrand<ExpansionType, PosFuncType, PointType, CoeffsType, MemorySpace> integrand(cache,
@@ -977,11 +1021,17 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                                                                                        DerivativeFlags::None);
 
         quad.Integrate(workspace, integrand, 0, 1, &output);
+        // std::cout << "output 1 = " << output << std::endl;
 
         // Finish filling in the cache for an evaluation of the expansion with x_d=0
+        // std::cout << "pt = " << std::endl;
+        // for(unsigned int i=0; i<pt.size(); ++i)
+        //     std::cout << "  " << pt(i);
+        // std::cout << std::endl;
+
         expansion.FillCache2(cache, pt, 0.0, DerivativeFlags::None);
         output += expansion.Evaluate(cache, coeffs);
-
+        // std::cout << "output 2 = " << output << std::endl;
         return output;
     }
 
@@ -1003,115 +1053,115 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
      @tparam PointType The type of the point.  Typically either Kokkos::View<double*> or some subview with a similar 1d signature.
      @tparam CoeffsType The type of the coefficients.  Typically Kokkos::View<double*> or a similarly structured subview.
      */
-    template<typename PointType, typename CoeffsType>
-    KOKKOS_FUNCTION static double InverseSingleBracket(double*                    cache,
-                                                       double*                    workspace,
-                                                       PointType           const& pt,
-                                                       double                     yd,
-                                                       CoeffsType          const& coeffs,
-                                                       const double               xtol,
-                                                       const double               ftol,
-                                                       QuadratureType      const& quad,
-                                                       ExpansionType       const& expansion)
-    {
-        double stepSize=1.0;
-        const unsigned int maxIts = 10000;
-
-        // First, we need to find two points that bound the solution.
-        double xlb, xub;
-        double ylb, yub;
-        double xb, xf; // Bisection point and regula falsi point
-        double xc, yc;
-
-        xlb = pt(pt.extent(0)-1);
-        ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
-
-        // We actually found an upper bound...
-        if(ylb>yd){
-
-            mpart::simple_swap(ylb,yub);
-            mpart::simple_swap(xlb,xub);
-
-            // Now find a lower bound...
-            unsigned int i;
-            for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
-                xlb = xub-stepSize;
-                ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
-                if(ylb>yd){
-                    mpart::simple_swap(ylb,yub);
-                    mpart::simple_swap(xlb,xub);
-                    stepSize *= 2.0;
-                }else{
-                    break;
-                }
-            }
-            if(i>maxIts)
-                ProcAgnosticError<MemorySpace, std::runtime_error>::error("InverseSingleBracket: lower bound iterations exceed maxIts");
-
-        // We have a lower bound...
-        }else{
-            // Now find an upper bound...
-            unsigned int i;
-            for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
-                xub = xlb+stepSize;
-                yub = EvaluateSingle(workspace, cache, pt, xub, coeffs, quad, expansion);
-                if(yub<yd){
-                    mpart::simple_swap(ylb,yub);
-                    mpart::simple_swap(xlb,xub);
-                    stepSize *= 2.0;
-                }else{
-                    break;
-                }
-            }
-            if(i>maxIts)
-                ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: upper bound calculation exceeds maxIts");
-        }
-
-        assert(ylb<yub);
-        assert(xlb<xub);
-
-        // Bracketed search
-        const double k1 = 0.1;
-        const double k2 = 2.0;
-        const double nhalf = ceil(log2(0.5*(xub-xlb)/xtol));
-        const double n0 = 1.0;
-
-        double sigma, delta, rho;
-        unsigned int it;
-        for(it=0; it<maxIts; ++it){
-
-            xb = 0.5*(xub+xlb); // bisection point
-            xf = xlb - (yd-ylb)*(xub-xlb) / (yub-ylb); // regula-falsi point
-
-            sigma = ((xb-xf)>0)?1.0:-1.0; // sign(xb-xf)
-            delta = fmin(k1*pow((xub-xlb), k2), fabs(xb-xf));
-
-            xf += delta*sigma;
-
-            rho = fmin(xtol*pow(2.0, nhalf + n0 - it) - 0.5*(xub-xlb), fabs(xf - xb));
-            xc = xb - sigma*rho;
-
-            yc = EvaluateSingle(workspace, cache, pt, xc, coeffs, quad, expansion);
-
-            if(abs(yc-yd)<ftol){
-                return xc;
-            }else if(yc>yd){
-                mpart::simple_swap(yc,yub);
-                mpart::simple_swap(xc,xub);
-            }else{
-                mpart::simple_swap(yc,ylb);
-                mpart::simple_swap(xc,xlb);
-            }
-
-            // Check for convergence
-            if(((xub-xlb)<xtol)||((yub-ylb)<ftol))
-                break;
-        };
-
-        if(it>maxIts)
-            ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
-        return 0.5*(xub+xlb);
-    }
+    // template<typename PointType, typename CoeffsType>
+    // KOKKOS_FUNCTION static double InverseSingleBracket(double*                    cache,
+    //                                                    double*                    workspace,
+    //                                                    PointType           const& pt,
+    //                                                    double                     yd,
+    //                                                    CoeffsType          const& coeffs,
+    //                                                    const double               xtol,
+    //                                                    const double               ftol,
+    //                                                    QuadratureType      const& quad,
+    //                                                    ExpansionType       const& expansion)
+    // {
+    //     double stepSize=1.0;
+    //     const unsigned int maxIts = 10000;
+
+    //     // First, we need to find two points that bound the solution.
+    //     double xlb, xub;
+    //     double ylb, yub;
+    //     double xb, xf; // Bisection point and regula falsi point
+    //     double xc, yc;
+
+    //     xlb = pt(pt.extent(0)-1);
+    //     ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
+
+    //     // We actually found an upper bound...
+    //     if(ylb>yd){
+
+    //         mpart::simple_swap(ylb,yub);
+    //         mpart::simple_swap(xlb,xub);
+
+    //         // Now find a lower bound...
+    //         unsigned int i;
+    //         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
+    //             xlb = xub-stepSize;
+    //             ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
+    //             if(ylb>yd){
+    //                 mpart::simple_swap(ylb,yub);
+    //                 mpart::simple_swap(xlb,xub);
+    //                 stepSize *= 2.0;
+    //             }else{
+    //                 break;
+    //             }
+    //         }
+    //         if(i>maxIts)
+    //             ProcAgnosticError<MemorySpace, std::runtime_error>::error("InverseSingleBracket: lower bound iterations exceed maxIts");
+
+    //     // We have a lower bound...
+    //     }else{
+    //         // Now find an upper bound...
+    //         unsigned int i;
+    //         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
+    //             xub = xlb+stepSize;
+    //             yub = EvaluateSingle(workspace, cache, pt, xub, coeffs, quad, expansion);
+    //             if(yub<yd){
+    //                 mpart::simple_swap(ylb,yub);
+    //                 mpart::simple_swap(xlb,xub);
+    //                 stepSize *= 2.0;
+    //             }else{
+    //                 break;
+    //             }
+    //         }
+    //         if(i>maxIts)
+    //             ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: upper bound calculation exceeds maxIts");
+    //     }
+
+    //     assert(ylb<yub);
+    //     assert(xlb<xub);
+
+    //     // Bracketed search
+    //     const double k1 = 0.1;
+    //     const double k2 = 2.0;
+    //     const double nhalf = ceil(log2(0.5*(xub-xlb)/xtol));
+    //     const double n0 = 1.0;
+
+    //     double sigma, delta, rho;
+    //     unsigned int it;
+    //     for(it=0; it<maxIts; ++it){
+
+    //         xb = 0.5*(xub+xlb); // bisection point
+    //         xf = xlb - (yd-ylb)*(xub-xlb) / (yub-ylb); // regula-falsi point
+
+    //         sigma = ((xb-xf)>0)?1.0:-1.0; // sign(xb-xf)
+    //         delta = fmin(k1*pow((xub-xlb), k2), fabs(xb-xf));
+
+    //         xf += delta*sigma;
+
+    //         rho = fmin(xtol*pow(2.0, nhalf + n0 - it) - 0.5*(xub-xlb), fabs(xf - xb));
+    //         xc = xb - sigma*rho;
+
+    //         yc = EvaluateSingle(workspace, cache, pt, xc, coeffs, quad, expansion);
+
+    //         if(abs(yc-yd)<ftol){
+    //             return xc;
+    //         }else if(yc>yd){
+    //             mpart::simple_swap(yc,yub);
+    //             mpart::simple_swap(xc,xub);
+    //         }else{
+    //             mpart::simple_swap(yc,ylb);
+    //             mpart::simple_swap(xc,xlb);
+    //         }
+
+    //         // Check for convergence
+    //         if(((xub-xlb)<xtol)||((yub-ylb)<ftol))
+    //             break;
+    //     };
+
+    //     if(it>maxIts)
+    //         ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
+    //     return 0.5*(xub+xlb);
+    // }
 
     /** Give access to the underlying FixedMultiIndexSet
      * @return The FixedMultiIndexSet
diff --git a/MParT/Utilities/Miscellaneous.h b/MParT/Utilities/Miscellaneous.h
index f1dfa790..6647bf26 100644
--- a/MParT/Utilities/Miscellaneous.h
+++ b/MParT/Utilities/Miscellaneous.h
@@ -22,7 +22,7 @@ namespace mpart{
     template<typename MemorySpace, typename ErrorType>
     struct ProcAgnosticError {
         static void error(const char*) {
-            assert(0);
+            assert(false);
         }
     };
 
diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
index 5228174c..03e5c65d 100644
--- a/MParT/Utilities/RootFinding.h
+++ b/MParT/Utilities/RootFinding.h
@@ -18,12 +18,14 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
                                         double& xub, double& yub,
                                         const double yd)
 {
-    const unsigned int maxIts = 1000;
+    const unsigned int maxIts = 40;
     double stepSize = 1.0;
     
     ylb = f(xlb);
     yub = f(xub);
-    
+    //std::cout << "-1" << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
+            
+
     // We actually found an upper bound...
     if(ylb>yd){
     
@@ -35,17 +37,26 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
             xlb = xub-stepSize;
             ylb = f(xlb);
+            //std::cout << "\n\n1 " << i << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
+            
+            if(abs((yub-ylb)/(xub-xlb))<1e-12){
+                break;
+                //std::cout << "slope = " << (yub-ylb)/(xub-xlb) << std::endl;
+                //assert(abs(yub-ylb)>1e-10);
+            }
+
             if(ylb>yd){
-                mpart::simple_swap(ylb,yub);
-                mpart::simple_swap(xlb,xub);
+                yub = ylb;
+                xub = xlb;
                 stepSize *= 2.0;
             }else{
                 break;
             }
         }
 
-        if(i>=maxIts)
-            ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
+
+        // if(i>=maxIts)
+        //     ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
 
     // We have a lower bound...
     }else{
@@ -54,17 +65,24 @@ 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);
+            //std::cout << "\n\n2 " << i << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
+
+            // if(abs((yub-ylb)/(xub-xlb))<1e-10){
+            //     std::cout << "Function seems to be flat!" << std::endl;
+            //     std::cout << "slope = " << (yub-ylb)/(xub-xlb) << std::endl;
+            //     assert(abs(yub-ylb)>1e-10);
+            // }
             if(yub<yd){
-                mpart::simple_swap(ylb,yub);
-                mpart::simple_swap(xlb,xub);
+                ylb = yub;
+                xlb = xub;
                 stepSize *= 2.0;
             }else{
                 break;
             }
         }
 
-        if(i>=maxIts)
-            ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
+        // if(i>=maxIts)
+        //     ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
     }
  }
 
@@ -95,13 +113,15 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     double ylb, yub;
     double xc, yc;
 
+    //std::cout << "xlb = " << xlb << std::endl;
     xlb = xub = x0;
     ylb = yub = f(xlb);
 
     // Compute bounds
     FindBracket<MemorySpace>(f, xlb, ylb, xub, yub, yd);
-    assert(ylb<yub);
-    assert(xlb<xub);
+    
+    if ((ylb>yd)||(yub<yd))
+        return std::numeric_limits<double>::quiet_NaN();
 
     // Bracketed search
     const double k1 = 0.1;
@@ -111,6 +131,8 @@ 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);
@@ -127,8 +149,10 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
         if(((xub-xlb)<xtol)||((yub-ylb)<ftol)) break;
     };
 
+
     if(it>maxIts)
-        ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
+        return std::numeric_limits<double>::quiet_NaN();
+        //ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
 
     return 0.5*(xub+xlb);
 }

From 23a76046fc87b588aa41055fed670f85fb756242 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Tue, 17 Oct 2023 16:55:30 -0400
Subject: [PATCH 4/6] Improved robustness of error handling in bracketed
 inverse.

---
 MParT/MonotoneComponent.h       | 23 ++---------
 MParT/Utilities/Miscellaneous.h |  3 ++
 MParT/Utilities/RootFinding.h   | 67 +++++++++++++++++++--------------
 tests/Test_RootFinding.cpp      | 51 +++++++++++++++++++++----
 4 files changed, 88 insertions(+), 56 deletions(-)

diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h
index 6b6be68f..aa7744fb 100644
--- a/MParT/MonotoneComponent.h
+++ b/MParT/MonotoneComponent.h
@@ -395,6 +395,7 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
         auto functor = KOKKOS_CLASS_LAMBDA (typename Kokkos::TeamPolicy<ExecutionSpace>::member_type team_member) {
 
             unsigned int ptInd = team_member.league_rank () * team_member.team_size () + team_member.team_rank ();
+            int info;
 
             if(ptInd<numPts){
                 unsigned int xInd = ptInd;
@@ -404,38 +405,22 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                 // Create a subview containing x_{1:d-1}
                 auto pt = Kokkos::subview(xs, Kokkos::ALL(), xInd);
                 
-                //std::cout << "\n\nPt0 = \n" << std::endl;
+                // Check for NaNs.  If found, set output to nan and return
                 for(unsigned int ii=0; ii<pt.size(); ++ii){
-                    //std::cout << "  " << pt(ii);
                     if(std::isnan(pt(ii))){
-                        std::cout << "Warning: Found nan in point, not computing inverse!" << std::endl;
-                        //throw std::runtime_error("Found nan in point!");
+                        output(ptInd) = std::numeric_limits<double>::quiet_NaN();
                         return;
                     }
                 }
-                //std::cout << std::endl << std::endl;
 
                 // Fill in the cache with everything that doesn't depend on x_d
                 Kokkos::View<double*,MemorySpace> cache(team_member.thread_scratch(1), cacheSize);
                 expansion_.FillCache1(cache.data(), pt, DerivativeFlags::None);
 
-                // std::cout << "\n\nPt1 = \n" << std::endl;
-                // double max_val = 0.0;
-                // for(unsigned int ii=0; ii<pt.size(); ++ii){
-                //     std::cout << "  " << pt(ii);
-                //     max_val  = (abs(pt(ii))>max_val) ?  abs(pt(ii)) : max_val;
-                // }
-                // std::cout << std::endl << std::endl;
-                // if (max_val > 1e-4){
-                //     std::cout << "Max val = " << max_val;
-                //     assert(max_val < 1e4);
-                // }
                 // Compute the inverse
                 Kokkos::View<double*,MemorySpace> workspace(team_member.thread_scratch(1), workspaceSize);
                 auto eval = SingleEvaluator<decltype(pt),decltype(coeffs)>(workspace.data(), cache.data(), pt, coeffs, quad_, expansion_);
-                //std::cout << "Here 0" << std::endl;
-                output(ptInd) = RootFinding::InverseSingleBracket<MemorySpace>(ys(ptInd), eval, pt(pt.extent(0)-1), xtol, ytol);
-                //std::cout << "Here 1" << std::endl;
+                output(ptInd) = RootFinding::InverseSingleBracket<MemorySpace>(ys(ptInd), eval, pt(pt.extent(0)-1), xtol, ytol, info);
             }
         };
 
diff --git a/MParT/Utilities/Miscellaneous.h b/MParT/Utilities/Miscellaneous.h
index 6647bf26..b4e3b290 100644
--- a/MParT/Utilities/Miscellaneous.h
+++ b/MParT/Utilities/Miscellaneous.h
@@ -19,6 +19,9 @@ namespace mpart{
                           std::string                                 const& key,
                           std::string                                 const& defaultValue);
 
+    /** Provides a mechanism for raising exceptions in CPU code where recovery is possible 
+        and assertions in GPU code where exceptions aren't alllowed.
+     */
     template<typename MemorySpace, typename ErrorType>
     struct ProcAgnosticError {
         static void error(const char*) {
diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
index 03e5c65d..971ef235 100644
--- a/MParT/Utilities/RootFinding.h
+++ b/MParT/Utilities/RootFinding.h
@@ -11,20 +11,25 @@ KOKKOS_INLINE_FUNCTION void swapPair(T& x1, T& x2, T& y1, T& y2) {
     simple_swap(y1, y2);
 }
 
-/** Finds a bracket [xlb, xub] such that f(xlb)<yd and f(xub)>yd. */
+/** Finds a bracket [xlb, xub] such that f(xlb)<yd and f(xub)>yd. 
+ * The info argument can be used to detect when a bracket cannot be found.  Upon exit, a value of info=0
+ * indicates success while a negative value indicates failure.  info=-1 indicates that the function 
+ * seems to be perfectly flat and a root might not exist.  info=-2 indicates that the maximum number of 
+ * iterations (128) was exceeded.
+*/
 template<typename MemorySpace, typename FunctorType>
 KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
                                         double& xlb, double& ylb,
                                         double& xub, double& yub,
-                                        const double yd)
+                                        const double yd,
+                                        int& info)
 {
-    const unsigned int maxIts = 40;
+    const unsigned int maxIts = 128;
     double stepSize = 1.0;
-    
+    info = 0;
+
     ylb = f(xlb);
     yub = f(xub);
-    //std::cout << "-1" << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
-            
 
     // We actually found an upper bound...
     if(ylb>yd){
@@ -37,12 +42,10 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
             xlb = xub-stepSize;
             ylb = f(xlb);
-            //std::cout << "\n\n1 " << i << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
             
             if(abs((yub-ylb)/(xub-xlb))<1e-12){
+                info = -1;
                 break;
-                //std::cout << "slope = " << (yub-ylb)/(xub-xlb) << std::endl;
-                //assert(abs(yub-ylb)>1e-10);
             }
 
             if(ylb>yd){
@@ -54,9 +57,8 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
             }
         }
 
-
-        // if(i>=maxIts)
-        //     ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
+        if(i>=maxIts)
+            info = -2;
 
     // We have a lower bound...
     }else{
@@ -65,13 +67,13 @@ 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);
-            //std::cout << "\n\n2 " << i << ": " << xlb << ", " << ylb << ", " << xub << ", " << yub << " vs " << yd << std::endl;
+          
+            // Check to see if function is perfectly flat
+            if(abs((yub-ylb)/(xub-xlb))<1e-12){
+                info = -1;
+                break;
+            }
 
-            // if(abs((yub-ylb)/(xub-xlb))<1e-10){
-            //     std::cout << "Function seems to be flat!" << std::endl;
-            //     std::cout << "slope = " << (yub-ylb)/(xub-xlb) << std::endl;
-            //     assert(abs(yub-ylb)>1e-10);
-            // }
             if(yub<yd){
                 ylb = yub;
                 xlb = xub;
@@ -81,13 +83,14 @@ KOKKOS_INLINE_FUNCTION void FindBracket(FunctorType f,
             }
         }
 
-        // if(i>=maxIts)
-        //     ProcAgnosticError<MemorySpace,std::runtime_error>::error("FindBracket: Could not find initial bracket such that f(xlb)<yd and f(xub)>yd.");
+        if(i>=maxIts){
+            info = -2;
+        }
     }
  }
 
 KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, double ylb, double yub,
-                           double k1, double k2, double nhalf, double n0, int it, double xtol) {
+                                         double k1, double k2, double nhalf, double n0, int it, double xtol) {
 
         double xb = 0.5*(xub+xlb); // bisection point
         double xf = (xub*ylb - xlb*yub)/(ylb-yub); // regula-falsi point
@@ -102,8 +105,13 @@ KOKKOS_INLINE_FUNCTION double Find_x_ITP(double xlb, double xub, double yd, doub
         return xc;
 }
 
+/** Computes the inverse of a function using the ITP method.  
+ *  The info argument will be 0 upon successful completion and negative for failed inversions.  
+ *  A value of info=-2 indicates a failure to find a bracket that contains the root. In this case, a nan will be returned.
+ *  A value of info=-1 indicates that the maximum number of iterations was exceeded. 
+*/
 template<typename MemorySpace, typename FunctorType>
-KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, double x0, const double xtol, const double ftol)
+KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, double x0, const double xtol, const double ftol, int& info)
 {   
     double stepSize=1.0;
     const unsigned int maxIts = 10000;
@@ -112,16 +120,18 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
     double xlb, xub;
     double ylb, yub;
     double xc, yc;
+    info = 0;
 
-    //std::cout << "xlb = " << xlb << std::endl;
     xlb = xub = x0;
     ylb = yub = f(xlb);
 
-    // Compute bounds
-    FindBracket<MemorySpace>(f, xlb, ylb, xub, yub, yd);
-    
-    if ((ylb>yd)||(yub<yd))
+    // 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();
+    }
 
     // Bracketed search
     const double k1 = 0.1;
@@ -151,8 +161,7 @@ KOKKOS_INLINE_FUNCTION double InverseSingleBracket(double yd, FunctorType f, dou
 
 
     if(it>maxIts)
-        return std::numeric_limits<double>::quiet_NaN();
-        //ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
+        info = -1;
 
     return 0.5*(xub+xlb);
 }
diff --git a/tests/Test_RootFinding.cpp b/tests/Test_RootFinding.cpp
index b630ee35..67462966 100644
--- a/tests/Test_RootFinding.cpp
+++ b/tests/Test_RootFinding.cpp
@@ -30,52 +30,87 @@ TEST_CASE( "RootFindingUtils", "[RootFindingUtils]") {
         double yd = identity(xd);
         double xub = 2., yub = 2.;
         double xlb = 0., ylb = 0.;
-        FindBracket<HostSpace>(identity, xlb, ylb, xub, yub, yd);
+        int info = 0;
+        FindBracket<HostSpace>(identity, xlb, ylb, xub, yub, yd, info);
         CheckFoundBounds(identity, xlb, xd, xub, ylb, yd, yub);
+        CHECK(info==0);
     }
     SECTION("FindBracker sigmoid") {
         double xd = -0.5;
         double yd = sigmoid(xd);
         double xub =  2., yub = sigmoid(xub);
         double xlb =  0., ylb = sigmoid(xlb);
-        FindBracket<HostSpace>(sigmoid, xlb, ylb, xub, yub, yd);
+        int info = 0;
+        FindBracket<HostSpace>(sigmoid, xlb, ylb, xub, yub, yd, info);
         CheckFoundBounds(sigmoid, xlb, xd, xub, ylb, yd, yub);
+        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;
+        int info = 0;
+        auto f = [](double x){return -1.0;};
+        FindBracket<HostSpace>(f, xlb, ylb, xub, yub, yd, info);
+        CHECK(info==-1);
     }
     SECTION("Test Inverse Linear, low x0") {
         double xd = 0.5, yd = identity(xd);
         double x0 = 0.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, identity, x0, xtol, ftol);
+        int info = 0;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, identity, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
     }
     SECTION("Test Inverse Linear, high x0") {
         double xd = 0.5, yd = identity(xd);
         double x0 = 1.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, identity, x0, xtol, ftol);
+        int info = 0;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, identity, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
+    }
+    SECTION("Test Inverse Flat") {
+        double xd = 0.5, yd = -1.0;
+        double x0 = 0.0, xtol = 1e-5, ftol = 1e-5;
+        auto f = [](double x){return -1.0;};
+        int info = 0;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, f, x0, xtol, ftol, info);
+        CHECK( std::isnan(xd_found));
+        CHECK(info==-2);
     }
     SECTION("Test Inverse Sigmoid, low x0") {
         double xd = 0.5, yd = sigmoid(xd);
         double x0 = 0.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid, x0, xtol, ftol);
+        int info = 0;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
     }
     SECTION("Test Inverse Sigmoid, high x0") {
         double xd = 0.5, yd = sigmoid(xd);
         double x0 = 1.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid, x0, xtol, ftol);
+        int info;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
     }
     SECTION("Test Inverse Sigmoid Combo, low x0") {
         double xd = 0.5, yd = sigmoid_combo(xd);
         double x0 = -5.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid_combo, x0, xtol, ftol);
+        int info;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid_combo, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
     }
     SECTION("Test Inverse Sigmoid Combo, high x0") {
         double xd = 0.5, yd = sigmoid_combo(xd);
         double x0 = 5.0, xtol = 1e-5, ftol = 1e-5;
-        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid_combo, x0, xtol, ftol);
+        int info;
+        double xd_found = InverseSingleBracket<HostSpace>(yd, sigmoid_combo, x0, xtol, ftol, info);
         CHECK( xd_found == Approx(xd).epsilon(2*xtol));
+        CHECK(info==0);
     }
 
 }
\ No newline at end of file

From 892c746342c218421b4c3213799148f007625718 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Tue, 17 Oct 2023 17:21:45 -0400
Subject: [PATCH 5/6] Added additional check on slope.

---
 MParT/Utilities/RootFinding.h | 14 +++++++-------
 tests/Test_RootFinding.cpp    |  4 +---
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/MParT/Utilities/RootFinding.h b/MParT/Utilities/RootFinding.h
index 971ef235..384d755f 100644
--- 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
index 67462966..e9a68cb1 100644
--- 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;

From 280987189c82f86b2fad93d00325e22d3fd75079 Mon Sep 17 00:00:00 2001
From: Matthew Parno <mparno@solea.com>
Date: Wed, 18 Oct 2023 14:00:25 -0400
Subject: [PATCH 6/6] Cleaning up unnecessary comments.

---
 MParT/MonotoneComponent.h | 157 +-------------------------------------
 1 file changed, 1 insertion(+), 156 deletions(-)

diff --git a/MParT/MonotoneComponent.h b/MParT/MonotoneComponent.h
index aa7744fb..d5fab1d5 100644
--- a/MParT/MonotoneComponent.h
+++ b/MParT/MonotoneComponent.h
@@ -312,16 +312,6 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                      StridedVector<double, MemorySpace>              output,
                  std::map<std::string, std::string>                  options=std::map<std::string,std::string>())
     {
-        // std::cout << "xs.shape = " << xs.extent(0) << "," << xs.extent(1) << std::endl;
-        // std::cout << "xs = \n" << std::endl;
-        // for(unsigned int row=0; row<xs.extent(0); ++row){
-        //     for(unsigned int col=0; col<xs.extent(1); ++col){
-        //         std::cout << "  " << xs(row,col);
-        //     }
-        //     std::cout << std::endl;
-        // }
-        // std::cout << std::endl;
-
         // Extract the method from the options map
         std::string method;
         if(options.count("Method")){
@@ -986,16 +976,6 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                                                  QuadratureType    const& quad,
                                                  ExpansionType     const& expansion)
     {
-        // std::cout << "pt = " << std::endl;
-        // for(unsigned int i=0; i<pt.size(); ++i)
-        //     std::cout << "  " << pt(i);
-        // std::cout << std::endl;
-        // std::cout << "xd = " << xd << std::endl;
-        // std::cout << "Coefficients = " << std::endl;
-        // for(unsigned int i=0; i<coeffs.size(); ++i)
-        //     std::cout << "  " << coeffs(i);
-        // std::cout << std::endl;
-        
         double output = 0.0;
         // Compute the integral \int_0^1 g( \partial_D f(x_1,...,x_{D-1},t*x_d)) dt
         MonotoneIntegrand<ExpansionType, PosFuncType, PointType, CoeffsType, MemorySpace> integrand(cache,
@@ -1006,148 +986,13 @@ class MonotoneComponent : public ConditionalMapBase<MemorySpace>
                                                                                        DerivativeFlags::None);
 
         quad.Integrate(workspace, integrand, 0, 1, &output);
-        // std::cout << "output 1 = " << output << std::endl;
-
-        // Finish filling in the cache for an evaluation of the expansion with x_d=0
-        // std::cout << "pt = " << std::endl;
-        // for(unsigned int i=0; i<pt.size(); ++i)
-        //     std::cout << "  " << pt(i);
-        // std::cout << std::endl;
 
         expansion.FillCache2(cache, pt, 0.0, DerivativeFlags::None);
         output += expansion.Evaluate(cache, coeffs);
-        // std::cout << "output 2 = " << output << std::endl;
+        
         return output;
     }
 
-
-    /**
-     @brief Solves \f$y_D = T(x_1,\ldots,x_D)\f$ for \f$x_D\f$ at a single point using the ITP bracketing method.
-     @details Uses the [Interpolate Truncate and Project (ITP)](https://en.wikipedia.org/wiki/ITP_method) method to solve
-              \f$y_D = T(x_1,\ldots,x_D)\f$ for \f$x_D\f$.  This method is a bracketing method similar to bisection, but
-              with faster convergence.   The worst case performance of ITP, in terms of required evaluations of \f$T\f$,
-              is the same as bisection.
-
-     @param cache Memory set up by Kokkos for this evaluation.  The `expansion_.FillCache1` function must be
-                  called before calling this function.
-     @param pt Array of length \f$D\f$ containing the fixed values of \f$x_{1:D-1}\f$ and an initial guess for \f$x_D\f$ in the last component.
-     @param coeffs An array of coefficients for the expansion.
-     @param options
-     @return The value of \f$x_D\f$ solving \f$y_D = T(x_1,\ldots,x_D)\f$
-
-     @tparam PointType The type of the point.  Typically either Kokkos::View<double*> or some subview with a similar 1d signature.
-     @tparam CoeffsType The type of the coefficients.  Typically Kokkos::View<double*> or a similarly structured subview.
-     */
-    // template<typename PointType, typename CoeffsType>
-    // KOKKOS_FUNCTION static double InverseSingleBracket(double*                    cache,
-    //                                                    double*                    workspace,
-    //                                                    PointType           const& pt,
-    //                                                    double                     yd,
-    //                                                    CoeffsType          const& coeffs,
-    //                                                    const double               xtol,
-    //                                                    const double               ftol,
-    //                                                    QuadratureType      const& quad,
-    //                                                    ExpansionType       const& expansion)
-    // {
-    //     double stepSize=1.0;
-    //     const unsigned int maxIts = 10000;
-
-    //     // First, we need to find two points that bound the solution.
-    //     double xlb, xub;
-    //     double ylb, yub;
-    //     double xb, xf; // Bisection point and regula falsi point
-    //     double xc, yc;
-
-    //     xlb = pt(pt.extent(0)-1);
-    //     ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
-
-    //     // We actually found an upper bound...
-    //     if(ylb>yd){
-
-    //         mpart::simple_swap(ylb,yub);
-    //         mpart::simple_swap(xlb,xub);
-
-    //         // Now find a lower bound...
-    //         unsigned int i;
-    //         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
-    //             xlb = xub-stepSize;
-    //             ylb = EvaluateSingle(workspace, cache, pt, xlb, coeffs, quad, expansion);
-    //             if(ylb>yd){
-    //                 mpart::simple_swap(ylb,yub);
-    //                 mpart::simple_swap(xlb,xub);
-    //                 stepSize *= 2.0;
-    //             }else{
-    //                 break;
-    //             }
-    //         }
-    //         if(i>maxIts)
-    //             ProcAgnosticError<MemorySpace, std::runtime_error>::error("InverseSingleBracket: lower bound iterations exceed maxIts");
-
-    //     // We have a lower bound...
-    //     }else{
-    //         // Now find an upper bound...
-    //         unsigned int i;
-    //         for(i=0; i<maxIts; ++i){ // Could just be while(true), but want to avoid infinite loop
-    //             xub = xlb+stepSize;
-    //             yub = EvaluateSingle(workspace, cache, pt, xub, coeffs, quad, expansion);
-    //             if(yub<yd){
-    //                 mpart::simple_swap(ylb,yub);
-    //                 mpart::simple_swap(xlb,xub);
-    //                 stepSize *= 2.0;
-    //             }else{
-    //                 break;
-    //             }
-    //         }
-    //         if(i>maxIts)
-    //             ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: upper bound calculation exceeds maxIts");
-    //     }
-
-    //     assert(ylb<yub);
-    //     assert(xlb<xub);
-
-    //     // Bracketed search
-    //     const double k1 = 0.1;
-    //     const double k2 = 2.0;
-    //     const double nhalf = ceil(log2(0.5*(xub-xlb)/xtol));
-    //     const double n0 = 1.0;
-
-    //     double sigma, delta, rho;
-    //     unsigned int it;
-    //     for(it=0; it<maxIts; ++it){
-
-    //         xb = 0.5*(xub+xlb); // bisection point
-    //         xf = xlb - (yd-ylb)*(xub-xlb) / (yub-ylb); // regula-falsi point
-
-    //         sigma = ((xb-xf)>0)?1.0:-1.0; // sign(xb-xf)
-    //         delta = fmin(k1*pow((xub-xlb), k2), fabs(xb-xf));
-
-    //         xf += delta*sigma;
-
-    //         rho = fmin(xtol*pow(2.0, nhalf + n0 - it) - 0.5*(xub-xlb), fabs(xf - xb));
-    //         xc = xb - sigma*rho;
-
-    //         yc = EvaluateSingle(workspace, cache, pt, xc, coeffs, quad, expansion);
-
-    //         if(abs(yc-yd)<ftol){
-    //             return xc;
-    //         }else if(yc>yd){
-    //             mpart::simple_swap(yc,yub);
-    //             mpart::simple_swap(xc,xub);
-    //         }else{
-    //             mpart::simple_swap(yc,ylb);
-    //             mpart::simple_swap(xc,xlb);
-    //         }
-
-    //         // Check for convergence
-    //         if(((xub-xlb)<xtol)||((yub-ylb)<ftol))
-    //             break;
-    //     };
-
-    //     if(it>maxIts)
-    //         ProcAgnosticError<MemorySpace,std::runtime_error>::error("InverseSingleBracket: Bracket search iterations exceeds maxIts");
-    //     return 0.5*(xub+xlb);
-    // }
-
     /** Give access to the underlying FixedMultiIndexSet
      * @return The FixedMultiIndexSet
      */