HP provisional termination

src/solvers/hp/hp.cc

@@ -39,11 +39,47 @@ HPStrategySolve(const MixedStrategyProfile<double> &p_prior,
   const PathTracer tracer;
   const PathTracer::TraceDirection direction = PathTracer::TraceDirection::Positive;
   const size_t tracking_index = 1; // Track the first variable (t) for orientation
-  auto termination_condition = [](const Vector<double> &point) { return point[1] >= 1.5; };
-  auto criterion_function = [](const Vector<double> &point,
-                               const Vector<double> &tangent) -> double { return point[1] - 1.0; };
+  const double t_target = 1.0;
+  const double tol = 1e-8;
+  double last_t = 0.0;
+  bool has_crossed = false;
 
-  tracer.TracePath(
+  auto termination_condition = [t_target, &last_t, &has_crossed, tol,
+                                &system](const Vector<double> &point) {
+    const double t = point[1];
+
+    // Path tracer reaches tol
+    if (system.ExtractEquilibrium(point).GetMaxRegret() <= tol && t >= t_target - tol) {
+      return true;
+    }
+
+    if (t > t_target) {
+      if (!has_crossed) {
+        has_crossed = true;
+      }
+      else if (t > last_t + tol) { // Criterion function is not working; polish will do the job
+        return true;
+      }
+    }
+    else {
+      // Criterion function might take t back to being less than t_target
+      has_crossed = false;
+    }
+
+    last_t = t;
+    return false;
+  };
+
+  auto criterion_function = [t_target](const Vector<double> &point,
+                                       const Vector<double> &tangent) -> double {
+    return point[1] - t_target;
+  };
+
+  auto polishing_termination_condition = [tol, &system](const Vector<double> &point) -> bool {
+    return system.ExtractEquilibrium(point).GetMaxRegret() <= tol;
+  };
+
+  const auto tracing_result = tracer.TracePath(
       [&system](const Vector<double> &point, Vector<double> &lhs) { system.GetValue(point, lhs); },
       [&system](const Vector<double> &point, Matrix<double> &jac) {
         system.GetJacobian(point, jac);
@@ -55,6 +91,20 @@ HPStrategySolve(const MixedStrategyProfile<double> &p_prior,
       },
       criterion_function, NullCriterionBracketFunction, p_cancel);
 
+  const PolishResult polishing_result = PolishPoint(
+      [&system](const Vector<double> &point, Vector<double> &lhs) { system.GetValue(point, lhs); },
+      [&system](const Vector<double> &point, Matrix<double> &jac) {
+        system.GetJacobian(point, jac);
+      },
+      x, t_target, 1, polishing_termination_condition, 100,
+      [&system, &p_onEvent](const Vector<double> &point) {
+        const MixedStrategyProfile<double> profile = system.ExtractEquilibrium(point);
+        p_onEvent(HPStepEvent{.profile = profile, .t = point[1]});
+      });
+
+  if (!polishing_result.status) {
+    return {};
+  }
   const MixedStrategyProfile<double> equilibrium = system.ExtractEquilibrium(x);
   p_onEquilibrium(equilibrium);
   equilibria.push_back(equilibrium);

src/solvers/path/path.cc

@@ -132,21 +132,25 @@ TracePathResult PathTracer::TracePath(
     CallbackFunctionType p_callback, CriterionFunctionType p_criterion,
     CriterionBracketFunctionType p_criterionBracket, const CancelToken &p_cancel) const
 {
-  const double c_tol = 1.0e-4;   // tolerance for corrector iteration
-  const double c_maxDist = 0.4;  // maximal distance to curve
-  const double c_maxContr = 0.6; // maximal contraction rate in corrector
-  const double c_eta = 0.1;      // perturbation to avoid cancellation
-                                 // in calculating contraction rate
-  double h = m_hStart;           // initial stepsize
-  const double c_hmin = 1.0e-8;  // minimal stepsize
-  const int c_maxIter = 100;     // maximum iterations in corrector
+  const double c_tol = 1.0e-4;       // tolerance for corrector iteration
+  const double c_maxDist = 0.4;      // maximal distance to curve
+  const double c_maxContr = 0.6;     // maximal contraction rate in corrector
+  const double c_eta = 0.1;          // perturbation to avoid cancellation
+                                     // in calculating contraction rate
+  double h = m_hStart;               // initial stepsize
+  const double c_hmin = 1.0e-8;      // minimal stepsize
+  const int c_maxIter = 100;         // maximum iterations in corrector
+  const double c_newtonTol = 1.0e-8; // tolerance for Newton convergence
 
   bool newton = false;             // using Newton steplength (for zero-finding)
   const double c_pert = 0.0000001; // The size of perturbation to apply to avoid bifurcation traps
   double pert = 0.0;               // The current version of the perturbation being applied
   double pert_countdown = 0.0;     // How much longer (in arclength) to apply perturbation
   const double c_orientTol = 1.0e-8; // tolerance for detecting change in orientation
 
+  const double b_tol = 1.0e-10; // Tolerance for perturbing the b matrix in case of singularity
+  const double b_pert = 1.0e-8; // Perturbation of the b matrix in case of singularity
+
   Vector<double> u(x.size());
   // t is current tangent at x; newT is tangent at u, which is the next point.
   Vector<double> t(x.size()), newT(x.size());
@@ -158,11 +162,22 @@ TracePathResult PathTracer::TracePath(
   QRDecomp(b, q);
   q.GetRow(q.NumRows(), t);
   p_callback(x);
+  int steps = 0;
+
+  auto stepsizeBelowMinimum = [&]() -> TracePathResult {
+    if (newton && std::abs(p_criterion(x, t)) < c_newtonTol) {
+      return {x, true,
+              "Path following terminated successfully at point satisfying criterion function.",
+              steps};
+    }
+    return {x, false, "Stepsize fell below minimum threshold.", steps};
+  };
+
   bool first_step = true;
   double omega = (p_direction == TraceDirection::Positive) ? 1.0 : -1.0;
 
   if (p_trackingIndex > x.size() || p_trackingIndex < 1) {
-    return {x, false, "Tracking index exceeds dimension of point vector."};
+    return {x, false, "Tracking index exceeds dimension of point vector.", steps};
   }
 
   while (!p_terminate(x)) {
@@ -171,16 +186,17 @@ TracePathResult PathTracer::TracePath(
     bool accept = true;
 
     if (std::abs(h) <= c_hmin) {
-      return {x, false, "Stepsize fell below minimum threshold."};
+      return stepsizeBelowMinimum();
     }
 
     if (first_step) {
       if (std::abs(t[p_trackingIndex]) <= c_orientTol) {
-        return {x, false, "Initial tangent vector is orthogonal to path-following direction."};
+        return {x, false, "Initial tangent vector is orthogonal to path-following direction.",
+                steps};
       }
       // Ensure that the tangent is oriented in the same direction as
       // the path-following direction.
-      else if (t[p_trackingIndex] < -c_orientTol) {
+      if (t[p_trackingIndex] < -c_orientTol) {
         omega *= -1.0;
       }
       first_step = false;
@@ -195,6 +211,18 @@ TracePathResult PathTracer::TracePath(
     p_jacobian(u, b);
     QRDecomp(b, q);
 
+    // Perturb the b matrix if it is singular or nearly singular
+    for (size_t i = 1; i < b.NumRows(); i++) {
+      if (std::abs(b(i, i)) < b_tol) {
+        if (b(i, i) < 0) {
+          b(i, i) -= b_pert;
+        }
+        else {
+          b(i, i) += b_pert;
+        }
+      }
+    }
+
     int iter = 1;
     double disto = 0.0;
     while (true) {
@@ -226,7 +254,7 @@ TracePathResult PathTracer::TracePath(
       disto = dist;
       iter++;
       if (iter > c_maxIter) {
-        return {x, false, "Maximum iterations exceeded."};
+        return {x, false, "Maximum iterations exceeded.", steps};
       }
     }
 
@@ -240,15 +268,15 @@ TracePathResult PathTracer::TracePath(
       // is oriented in the same direction as we were originally following
       if (pert_countdown == 0.0) {
         pert = c_pert;
-        pert_countdown = abs(2 * h);
+        pert_countdown = std::abs(2 * h);
       }
       accept = false;
     }
 
     if (!accept) {
       h /= m_maxDecel; // PC not accepted; change stepsize and retry
       if (std::abs(h) <= c_hmin) {
-        return {x, false, "Stepsize fell below minimum threshold."};
+        return stepsizeBelowMinimum();
       }
       continue;
     }
@@ -279,18 +307,98 @@ TracePathResult PathTracer::TracePath(
     x = u;
     t = newT;
     p_callback(x);
+    steps++;
 
     if (pert_countdown > 0.0) {
       // If we are currently perturbing in the neighborhood of a bifurcation, check to see
       // whether we think we are likely past it, and switch off if we are.
-      pert_countdown -= abs(h);
+      pert_countdown -= std::abs(h);
       if (pert_countdown < 0.0) {
         pert = 0.0;
         pert_countdown = 0.0;
       }
     }
   }
-  return {x, true, "Path tracing terminated successfully."};
+  return {x, true, "Path tracing terminated successfully.", steps};
+}
+
+PolishResult PolishPoint(std::function<void(const Vector<double> &, Vector<double> &)> p_function,
+                         std::function<void(const Vector<double> &, Matrix<double> &)> p_jacobian,
+                         Vector<double> &x, double fixed_value, size_t fixed_index,
+                         TerminationFunctionType p_terminate, int max_iter,
+                         CallbackFunctionType p_callback)
+{
+  x[fixed_index] = fixed_value;
+
+  const Vector<double> original_x = x;
+
+  const size_t N = x.size() - 1;
+  Vector<double> y(N);               // Equations results
+  Matrix<double> jac_full(N + 1, N); // Full Jacobian matrix (N+1 unknowns, N equations)
+  Matrix<double> jac_square(N, N);   // Jacobian matrix with fixed_index row removed
+  Matrix<double> Q(N, N);            // Orthogonal matrix from QR decomposition
+  Vector<double> x_reduced(N);       // Reduced x vector with fixed_index removed
+
+  double const eq_tol = 1e-2;
+
+  int steps = 0;
+  double dist = 0.0;
+
+  while (!p_terminate(x)) {
+    if (steps >= max_iter) {
+      return {x, false, "Polishing exceeded maximum iterations.", steps};
+    }
+
+    p_function(x, y);
+    p_jacobian(x, jac_full);
+
+    size_t row_index = 1;
+    for (size_t i = 1; i <= N + 1; ++i) { // Newton step expects the transposed Jacobian
+      if (i != fixed_index) {
+        for (size_t j = 1; j <= N; ++j) {
+          jac_square(row_index, j) = jac_full(i, j);
+        }
+        row_index++;
+      }
+    }
+
+    // Reduced x vector removing fixed_index
+    size_t temp_idx = 1;
+    for (size_t i = 1; i <= N + 1; ++i) {
+      if (i != fixed_index) {
+        x_reduced[temp_idx++] = x[i];
+      }
+    }
+
+    QRDecomp(jac_square, Q);
+
+    // Solve jac_square * x_reduced = -y
+    NewtonStep(Q, jac_square, x_reduced, y, dist);
+
+    // Update x, keeping fixed_index constant
+    temp_idx = 1;
+    for (size_t i = 1; i <= N + 1; ++i) {
+      if (i != fixed_index) {
+        x[i] = x_reduced[temp_idx++];
+      }
+    }
+
+    steps++;
+
+    if (p_callback) {
+      p_callback(x);
+    }
+  }
+
+  // Checking that the profile satisfies the system of equations
+  for (size_t i = 1; i <= N; ++i) {
+    if (std::abs(y[i]) > eq_tol) {
+      x = original_x;
+      return {x, false, "Polishing converged to an invalid mathematical state. Reverted.", steps};
+    }
+  }
+
+  return {x, true, "Polishing terminated successfully.", steps};
 }
 
 } // end namespace Gambit

src/solvers/path/path.h

@@ -63,8 +63,15 @@ struct TracePathResult {
   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
+  int steps;           // Step at which the tracing terminated
 };
 
+struct PolishResult {
+  Vector<double> final_point;
+  bool status; // true if polishing terminated successfully, false if it terminated due to error
+  std::string message; // error message if status is false
+  int steps;           // Step at which the polishing terminated
+};
 //
 // 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
@@ -96,6 +103,15 @@ class PathTracer {
   double m_maxDecel{1.1}, m_hStart{0.03};
 };
 
+// This function reduces the regret of a point that is close to an equilibrium that has been found
+// by the path-following algorithm.  Fixing the value of a component of the point, it uses a
+// Newton-type method to find a nearby point with lower regret.
+PolishResult PolishPoint(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 fixed_value, size_t fixed_index,
+                         TerminationFunctionType p_terminate, int max_iter = 100,
+                         CallbackFunctionType p_callback = NullCallbackFunction);
+
 } // end namespace Gambit
 
 #endif // GAMBIT_SOLVERS_LOGIT_PATH_H