diff --git a/doc/algorithms.rst b/doc/algorithms.rst index 317306b029..db71bc87d2 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -22,6 +22,42 @@ Algorithm Description :ref:`hp` Compute a specific Nash equilibrium using a homotopy path-following method :py:func:`pygambit.nash.hp_solve` :ref:`gambit-hp ` ================ =========================================================================== ======================================== ========================================== +.. _pygambit-nash-maxregret: + +Acceptance criteria for approximate Nash equilibria +---------------------------------------------------- + +Some methods for computing Nash equilibria operate using floating-point arithmetic and/or +generate candidate equilibrium profiles using methods which involve some form of successive +approximation. The outputs of these methods therefore are in general +:math:`\varepsilon`-equilibria, for some positive :math:`\varepsilon`: a strategy profile at +which no player can gain more than :math:`\varepsilon` in expected payoff by unilaterally +deviating. Every Nash equilibrium is an :math:`\varepsilon`-equilibrium with +:math:`\varepsilon = 0`. + +To provide a uniform interface across methods, where relevant Gambit provides a parameter +`maxregret`, which specifies the acceptance criterion for labeling the output of the +algorithm as an equilibrium. This parameter is interpreted *proportionally* to the range of +payoffs in the game: any profile returned as an equilibrium is guaranteed to be an +:math:`\varepsilon`-equilibrium, for :math:`\varepsilon` no more than `maxregret` times the +difference of the game's maximum and minimum payoffs. For example, with the default +`maxregret` of :math:`10^{-8}` in a game whose payoffs range over 4 units, any equilibrium +returned is guaranteed to have a regret, measured directly in the game's own payoffs, of no +more than :math:`4 \times 10^{-8}`. + +Expressing `maxregret` scaled by the game's payoffs in this way standardises the behavior of +methods across games: for instance, doubling all the payoffs in a game does not change the +`maxregret` value needed to obtain equilibria of comparable quality. + +Methods differ in the guarantees they offer once a `maxregret` criterion is specified. +Globally-convergent methods, such as :ref:`logit` and :ref:`gnm`, are guaranteed to eventually +satisfy any `maxregret` criterion, though a tighter criterion generally requires more +computation. Other methods, such as :ref:`liap`, are not globally convergent, and may fail to +find any equilibrium satisfying the criterion from a given starting point. + +See the :doc:`stripped-down poker tutorial ` for a worked +example comparing `maxregret` across several methods. + .. _enumpure: enumpure diff --git a/doc/tools.hp.rst b/doc/tools.hp.rst index b2eac367c4..86770fa328 100644 --- a/doc/tools.hp.rst +++ b/doc/tools.hp.rst @@ -24,6 +24,11 @@ different equilibria being found. Prints a help message listing the available options. +.. cmdoption:: -m + + Specify the maximum regret criterion for acceptance as an approximate Nash equilibrium + (default is 1e-8). See :ref:`pygambit-nash-maxregret` for interpretation and guidance. + .. cmdoption:: -n Randomly generate the specified number of prior distributions. diff --git a/doc/tutorials/03_stripped_down_poker.ipynb b/doc/tutorials/03_stripped_down_poker.ipynb index fab307ac26..cc890eec5a 100644 --- a/doc/tutorials/03_stripped_down_poker.ipynb +++ b/doc/tutorials/03_stripped_down_poker.ipynb @@ -675,29 +675,7 @@ "cell_type": "markdown", "id": "b2867dca", "metadata": {}, - "source": [ - "Acceptance criteria for Nash equilibria\n", - "---------------------------------------\n", - "\n", - "Some methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations.\n", - "The outputs of these methods therefore are in general $\\varepsilon$-equilibria, for some positive $\\varepsilon$.\n", - "\n", - "$\\varepsilon$-equilibria (from [Wikipedia](https://en.wikipedia.org/wiki/Epsilon-equilibrium)):\n", - "\n", - "> In game theory, an epsilon-equilibrium, or near-Nash equilibrium, is a strategy profile that approximately satisfies the condition of Nash equilibrium. In a Nash equilibrium, no player has an incentive to change his behavior. In an approximate Nash equilibrium, this requirement is weakened to allow the possibility that a player may have a small incentive to do something different.\n", - "\n", - "> Given a game and a real non-negative parameter $\\varepsilon$, a strategy profile is said to be an $\\varepsilon$-equilibrium if it is not possible for any player to gain more than $\\varepsilon$ in expected payoff by unilaterally deviating from his strategy. Every Nash Equilibrium is an $\\varepsilon$-equilibrium where $\\varepsilon = 0$.\n", - "\n", - "\n", - "To provide a uniform interface across methods, where relevant Gambit provides a parameter\n", - "`maxregret`, which specifies the acceptance criterion for labeling the output of the\n", - "algorithm as an equilibrium.\n", - "This parameter is interpreted *proportionally* to the range of payoffs in the game.\n", - "Any profile returned as an equilibrium is guaranteed to be an $\\varepsilon$-equilibrium, for $\\varepsilon$ no more than `maxregret`\n", - "times the difference of the game's maximum and minimum payoffs.\n", - "\n", - "As an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" - ] + "source": "Acceptance criteria for Nash equilibria\n---------------------------------------\n\nSome methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations, and so return only approximate equilibria.\nGambit expresses how close an approximation must be, before it is reported as an equilibrium, via a `maxregret` parameter, interpreted *proportionally* to the range of payoffs in the game.\nSee [Acceptance criteria for approximate Nash equilibria](../algorithms.html#pygambit-nash-maxregret) in the algorithms documentation for the full explanation of `maxregret` and the guarantees it provides.\n\nAs an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" }, { "cell_type": "code", diff --git a/src/gui/dlnash.cc b/src/gui/dlnash.cc index 702b3b63d2..5572e7b1f8 100644 --- a/src/gui/dlnash.cc +++ b/src/gui/dlnash.cc @@ -175,7 +175,8 @@ wxString ExternalCommand(const NashComputationSpec &p_spec) method.localNewtonMaxIterations); } else if constexpr (std::is_same_v) { - return prefix + wxString::Format("hp -d 10 -n %d", method.priors); + return prefix + + wxString::Format("hp -d 10 -n %d -m %.17g", method.priors, method.maxRegret); } else if constexpr (std::is_same_v) { return prefix + wxString::Format("ipa -d 10 -n %d", method.perturbations); @@ -274,7 +275,8 @@ wxString ParameterDescription(const NashMethodSpec &p_method) method.localNewtonMaxIterations); } else if constexpr (std::is_same_v) { - return wxString::Format(" (%d random priors)", method.priors); + return wxString::Format(" (%d random priors; maximum regret %.4g)", method.priors, + method.maxRegret); } else if constexpr (std::is_same_v) { return wxString::Format(" (%d perturbation)", method.perturbations); diff --git a/src/gui/nashspec.cc b/src/gui/nashspec.cc index 77e5d17f31..54be81996f 100644 --- a/src/gui/nashspec.cc +++ b/src/gui/nashspec.cc @@ -104,7 +104,7 @@ std::optional HPNashSpec::MakeSolver(NashRepresentation) const for (const auto &prior : NewRandomStrategyProfiles(p_game, spec.priors)) { p_cancel.Check(); Nash::HPStrategySolve( - prior, + prior, spec.maxRegret, [&p_callback](const MixedStrategyProfile &p) { p_callback(ComputedProfile(p)); }, Nash::NullHPEventCallback, p_cancel); } diff --git a/src/gui/nashspec.h b/src/gui/nashspec.h index de07b4d933..d138cf9f8b 100644 --- a/src/gui/nashspec.h +++ b/src/gui/nashspec.h @@ -102,6 +102,7 @@ struct HPNashSpec { // path-tracing can surface more), so this defaults high, like LiapNashSpec's // startingPoints, rather than to 1. int priors{10}; + double maxRegret{1.0e-8}; std::optional MakeSolver(NashRepresentation) const; }; diff --git a/src/pygambit/cli/hp.py b/src/pygambit/cli/hp.py index 0daa664b61..41acc7890c 100644 --- a/src/pygambit/cli/hp.py +++ b/src/pygambit/cli/hp.py @@ -80,6 +80,14 @@ default=None, help="file containing prior distributions (mutually exclusive with -n)", ) +@click.option( + "-m", + "maxregret", + default=1.0e-8, + show_default=True, + type=float, + help="maximum regret acceptable as a proportion of the range of payoffs in the game", +) @click.option("-q", "--quiet", is_flag=True, help="quiet mode (suppresses banner)") @click.option( "-V", @@ -95,6 +103,7 @@ def main( n_priors: int | None, seed: int | None, start_file: str | None, + maxregret: float, quiet: bool, verbose: bool, ) -> None: @@ -109,7 +118,7 @@ def render_event(event: gbt.HPStepEvent) -> None: prior = prior.as_float() if verbose: click.echo(render_profile_csv(prior, "prior", decimals)) - result = gbt.nash.hp_solve(prior, event_callback=render_event) + result = gbt.nash.hp_solve(prior, maxregret=maxregret, event_callback=render_event) for eq in result.equilibria: click.echo(render_profile_csv(eq, "NE", decimals)) diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index 8dce490284..1c8c16f73f 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -701,5 +701,5 @@ cdef extern from "nash.h": LogitEventCallbackType[c_LogitQREMixedStrategyProfile] ) except + stdlist[c_MixedStrategyProfile[double]] HPStrategySolveWrapper( - c_MixedStrategyProfile[double], HPEventCallbackType + c_MixedStrategyProfile[double], double, HPEventCallbackType ) except +RuntimeError diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 75dbfb2e20..815a27d786 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -143,8 +143,9 @@ LogitStrategyEstimateWrapper(std::shared_ptr> p_fre } std::list> -HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, +HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, double p_maxRegret, Nash::HPEventCallbackType p_onEvent = Nash::NullHPEventCallback) { - return Nash::HPStrategySolve(p_prior, Nash::NullStrategyCallback, p_onEvent); + return Nash::HPStrategySolve(p_prior, p_maxRegret, Nash::NullStrategyCallback, + p_onEvent); } diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index e319846a3e..595a422235 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -914,8 +914,9 @@ def _logit_behavior_branch(game: Game, def _hp_strategy_solve( prior: MixedStrategyProfileDouble, + maxregret: float, event_callback: object = None, ) -> list[MixedStrategyProfileDouble]: return _convert_mspd(HPStrategySolveWrapper( - deref(prior.profile), MakeHPEventCallback(event_callback) + deref(prior.profile), maxregret, MakeHPEventCallback(event_callback) )) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 32aa3aa07c..8789079022 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -1113,6 +1113,7 @@ def logit_solve( def hp_solve( prior: libgbt.MixedStrategyProfileDouble, + maxregret: float = 1.0e-8, event_callback: Callable[[libgbt.HPStepEvent], None] | None = None, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -1125,6 +1126,11 @@ def hp_solve( prior : MixedStrategyProfileDouble The prior distribution over strategies. + maxregret : float, default 1e-8 + The acceptance criterion for approximate Nash equilibrium; the maximum + regret of any player must be no more than `maxregret` times the + difference of the maximum and minimum payoffs of the game + event_callback : Callable[[HPStepEvent], None], optional If specified, called with each point traced along the homotopy path, and the homotopy parameter ``t`` at which it was reached, on the way @@ -1137,12 +1143,14 @@ def hp_solve( res : NashComputationResult The result represented as a ``NashComputationResult`` object. """ - equilibria = libgbt._hp_strategy_solve(prior, event_callback) + if maxregret <= 0.0: + raise ValueError("hp_solve(): maxregret argument must be positive") + equilibria = libgbt._hp_strategy_solve(prior, maxregret, event_callback) return NashComputationResult( game=prior.game, method="hp", rational=False, use_strategic=True, equilibria=equilibria, - parameters={}, + parameters={"maxregret": maxregret}, ) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 25b06325c2..95e640b875 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -27,7 +27,7 @@ namespace Gambit::Nash { std::list> -HPStrategySolve(const MixedStrategyProfile &p_prior, +HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret, StrategyCallbackType p_onEquilibrium, HPEventCallbackType p_onEvent, const CancelToken &p_cancel) { @@ -39,11 +39,49 @@ HPStrategySolve(const MixedStrategyProfile &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 &point) { return point[1] >= 1.5; }; - auto criterion_function = [](const Vector &point, - const Vector &tangent) -> double { return point[1] - 1.0; }; + const double t_target = 1.0; + double last_t = 0.0; + bool has_crossed = false; - tracer.TracePath( + auto termination_condition = [t_target, &last_t, &has_crossed, p_maxRegret, + &system](const Vector &point) { + const double t = point[1]; + + // Path tracer reaches maximum acceptable regret + if (system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret && + t >= t_target - p_maxRegret) { + return true; + } + + if (t > t_target) { + if (!has_crossed) { + has_crossed = true; + } + else if (t > + last_t + p_maxRegret) { // 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 &point, + const Vector &tangent) -> double { + return point[1] - t_target; + }; + + auto polishing_termination_condition = [p_maxRegret, + &system](const Vector &point) -> bool { + return system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret; + }; + + const auto tracing_result = tracer.TracePath( [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, [&system](const Vector &point, Matrix &jac) { system.GetJacobian(point, jac); @@ -55,7 +93,24 @@ HPStrategySolve(const MixedStrategyProfile &p_prior, }, criterion_function, NullCriterionBracketFunction, p_cancel); + const PolishResult polishing_result = PolishPoint( + [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, + [&system](const Vector &point, Matrix &jac) { + system.GetJacobian(point, jac); + }, + x, t_target, 1, polishing_termination_condition, 100, + [&system, &p_onEvent](const Vector &point) { + const MixedStrategyProfile profile = system.ExtractEquilibrium(point); + p_onEvent(HPStepEvent{.profile = profile, .t = point[1]}); + }); + + if (!polishing_result.status) { + return {}; + } const MixedStrategyProfile equilibrium = system.ExtractEquilibrium(x); + if (equilibrium.GetMaxRegret() > p_maxRegret) { + return {}; + } p_onEquilibrium(equilibrium); equilibria.push_back(equilibrium); return equilibria; diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index 4bc5a7280f..b600a1e020 100644 --- a/src/solvers/hp/hp.h +++ b/src/solvers/hp/hp.h @@ -45,7 +45,7 @@ inline void NullHPEventCallback(const HPEvent &) {} /// @brief Compute a Nash equilibrium of a game using the homotopy method of /// Herings and Peeters (2001) std::list> -HPStrategySolve(const MixedStrategyProfile &p_prior, +HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret = 1.0e-8, StrategyCallbackType p_onEquilibrium = NullStrategyCallback, HPEventCallbackType p_onEvent = NullHPEventCallback, const CancelToken &p_cancel = CancelToken()); diff --git a/src/solvers/path/path.cc b/src/solvers/path/path.cc index 0a6fcba301..7db29c9bbd 100644 --- a/src/solvers/path/path.cc +++ b/src/solvers/path/path.cc @@ -132,14 +132,15 @@ 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 @@ -147,6 +148,9 @@ TracePathResult PathTracer::TracePath( 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 u(x.size()); // t is current tangent at x; newT is tangent at u, which is the next point. Vector 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,7 +268,7 @@ 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; } @@ -248,7 +276,7 @@ TracePathResult PathTracer::TracePath( 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; } @@ -266,9 +294,10 @@ TracePathResult PathTracer::TracePath( p_criterionBracket(x, u); } - if (newton) { + const double diff = p_criterion(u, newT) - p_criterion(x, t); + if (newton && std::abs(diff) > c_newtonTol) { // Newton-type steplength adaptation, secant method - h *= -p_criterion(u, newT) / (p_criterion(u, newT) - p_criterion(x, t)); + h *= -p_criterion(u, newT) / diff; } else { // Standard steplength adaptation @@ -279,18 +308,88 @@ 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 &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &x, double fixed_value, size_t fixed_index, + TerminationFunctionType p_terminate, int max_iter, + CallbackFunctionType p_callback) +{ + x[fixed_index] = fixed_value; + + const size_t N = x.size() - 1; + Vector y(N); // Equations results + Matrix jac_full(N + 1, N); // Full Jacobian matrix (N+1 unknowns, N equations) + Matrix jac_square(N, N); // Jacobian matrix with fixed_index row removed + Matrix Q(N, N); // Orthogonal matrix from QR decomposition + Vector x_reduced(N); // Reduced x vector with fixed_index removed + + 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); + } + } + + // The loop only exits here once p_terminate(x) holds for the actual, current x; + // there is nothing further to validate against a separate tolerance. + return {x, true, "Polishing terminated successfully.", steps}; } } // end namespace Gambit diff --git a/src/solvers/path/path.h b/src/solvers/path/path.h index d05b1c9adf..43cf9aff12 100644 --- a/src/solvers/path/path.h +++ b/src/solvers/path/path.h @@ -63,8 +63,15 @@ struct TracePathResult { Vector 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 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 &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &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