Refactor TracePath to return a struct with relevant data

src/solvers/logit/efglogit.cc

@@ -338,6 +338,7 @@ LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret,
         return RegretTerminationFunction(game, p_point, p_regret);
       },
       [&callback](const Vector<double> &p_point) -> void { callback.AppendPoint(p_point); });
+
   return callback.GetProfiles();
 }
 
@@ -369,6 +370,7 @@ LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start,
         [lam](const Vector<double> &x, const Vector<double> &) -> double {
           return x.back() - lam;
         });
+
     ret.push_back(callback.GetProfiles().back());
   }
   return ret;

src/solvers/logit/path.cc

@@ -125,11 +125,12 @@ void NewtonStep(Matrix<double> &q, Matrix<double> &b, Vector<double> &u, Vector<
 // bifurcation point that the tracing gets stuck there as it is not possible
 // to find a small enough step size to avoid stepping over the bifurcation
 // point.
-void PathTracer::TracePath(
-    std::function<void(const Vector<double> &, Vector<double> &)> p_function,
-    std::function<void(const Vector<double> &, Matrix<double> &)> p_jacobian, Vector<double> &x,
-    double &p_omega, TerminationFunctionType p_terminate, CallbackFunctionType p_callback,
-    CriterionFunctionType p_criterion, CriterionBracketFunctionType p_criterionBracket) const
+TracePathResult
+PathTracer::TracePath(std::function<void(const Vector<double> &, Vector<double> &)> p_function,
+                      std::function<void(const Vector<double> &, Matrix<double> &)> p_jacobian,
+                      Vector<double> &x, double &p_omega, TerminationFunctionType p_terminate,
+                      CallbackFunctionType p_callback, CriterionFunctionType p_criterion,
+                      CriterionBracketFunctionType p_criterionBracket) const
 {
   const double c_tol = 1.0e-4;   // tolerance for corrector iteration
   const double c_maxDist = 0.4;  // maximal distance to curve
@@ -152,6 +153,7 @@ void PathTracer::TracePath(
   Matrix<double> b(x.size(), x.size() - 1);
   Matrix<double> q(x.size(), x.size());
 
+  // Obtain the tangent at the initial point
   p_jacobian(x, b);
   QRDecomp(b, q);
   q.GetRow(q.NumRows(), t);
@@ -161,7 +163,7 @@ void PathTracer::TracePath(
     bool accept = true;
 
     if (fabs(h) <= c_hmin) {
-      return;
+      return {x.back(), x, false, "Stepsize fell below minimum threshold."};
     }
 
     // Predictor step
@@ -204,7 +206,7 @@ void PathTracer::TracePath(
       disto = dist;
       iter++;
       if (iter > c_maxIter) {
-        return;
+        return {x.back(), x, false, "Maximum iterations exceeded."};
       }
     }
 
@@ -226,7 +228,7 @@ void PathTracer::TracePath(
     if (!accept) {
       h /= m_maxDecel; // PC not accepted; change stepsize and retry
       if (fabs(h) <= c_hmin) {
-        return;
+        return {x.back(), x, false, "Stepsize fell below minimum threshold."};
       }
       continue;
     }
@@ -268,6 +270,7 @@ void PathTracer::TracePath(
       }
     }
   }
+  return {x.back(), x, true, "Path tracing terminated successfully."};
 }
 
 } // end namespace Gambit

src/solvers/logit/path.h

@@ -58,6 +58,13 @@ using CallbackFunctionType = std::function<void(const Vector<double> &)>;
 
 inline void NullCallbackFunction(const Vector<double> &) {}
 
+struct TracePathResult {
+  double final_parameter_value;
+  Vector<double> final_point;
+  bool status; // true if path tracing terminated successfully, false if it terminated due to error
+  std::string message; // error message if status is false
+};
+
 //
 // This class implements a generic path-following algorithm for smooth curves.
 // It is based on the ideas and codes presented in Allgower and Georg's
@@ -74,7 +81,7 @@ class PathTracer {
   void SetStepsize(double p_hStart) { m_hStart = p_hStart; }
   double GetStepsize() const { return m_hStart; }
 
-  void
+  TracePathResult
   TracePath(std::function<void(const Vector<double> &, Vector<double> &)> p_function,
             std::function<void(const Vector<double> &, Matrix<double> &)> p_jacobian,
             Vector<double> &p_x, double &p_omega, TerminationFunctionType p_terminate,

src/tools/logit/logit.cc

@@ -255,4 +255,4 @@ int main(int argc, char *argv[])
     std::cerr << "Error: " << e.what() << std::endl;
     return 1;
   }
-}
+}