From 9bd154fab718c3827a175925342410eb63449332 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Fri, 3 Jul 2026 10:32:24 +0200 Subject: [PATCH 1/9] Add shell infrastructure for Herings-Peeters (2001) solver (#967) --- Makefile.am | 4 ++++ doc/references.bib | 10 ++++++++++ setup.py | 3 ++- src/pygambit/gambit.pxd | 4 ++++ src/pygambit/nash.pxi | 4 ++++ src/pygambit/nash.py | 29 +++++++++++++++++++++++++++++ src/solvers/hp/hp.cc | 38 ++++++++++++++++++++++++++++++++++++++ src/solvers/hp/hp.h | 32 ++++++++++++++++++++++++++++++++ 8 files changed, 123 insertions(+), 1 deletion(-) create mode 100644 src/solvers/hp/hp.cc create mode 100644 src/solvers/hp/hp.h diff --git a/Makefile.am b/Makefile.am index 5ef18af9bd..f2ea3fda3d 100644 --- a/Makefile.am +++ b/Makefile.am @@ -313,6 +313,10 @@ gtracer_SOURCES = \ src/solvers/gtracer/gnm.cc \ src/solvers/gtracer/ipa.cc +hp_SOURCES = \ + src/solvers/hp/hp.h \ + src/solvers/hp/hp.cc + if IS_WIN32 AM_LDFLAGS = -static -static-libgcc -static-libstdc++ endif diff --git a/doc/references.bib b/doc/references.bib index f39c14834c..36ae6110e6 100644 --- a/doc/references.bib +++ b/doc/references.bib @@ -38,6 +38,16 @@ @article{GovWil04 category = {articles_equilibria} } +@article{HerPee01, + author = {Herings, P. J.-J. and Peeters, R. J. A. P.}, + title = {A differentiable homotopy to compute {N}ash equilibria of n-person games}, + journal = {Economic Theory}, + volume = {18}, + pages = {159--185}, + year = {2001}, + category = {articles_equilibria} +} + @article{HalPas21, author = {Halpern, J. Y. and Pass, R.}, title = {Sequential equilibrium in games of imperfect recall}, diff --git a/setup.py b/setup.py index 9b47eb5333..945b827147 100644 --- a/setup.py +++ b/setup.py @@ -100,6 +100,7 @@ def run(self) -> None: cppgambit_gtracer = solver_library_config("cppgambit_gtracer", ["gtracer", "ipa", "gnm"]) cppgambit_simpdiv = solver_library_config("cppgambit_simpdiv", ["simpdiv"]) cppgambit_enumpoly = solver_library_config("cppgambit_enumpoly", ["nashsupport", "enumpoly"]) +cppgambit_hp = solver_library_config("cppgambit_hp", ["hp"]) libgambit = setuptools.Extension( @@ -117,7 +118,7 @@ def run(self) -> None: setuptools.setup( cmdclass={"build_py": GambitBuildPy}, libraries=[cppgambit_bimatrix, cppgambit_liap, cppgambit_logit, cppgambit_simpdiv, - cppgambit_gtracer, cppgambit_enumpoly, + cppgambit_gtracer, cppgambit_enumpoly, cppgambit_hp, cppgambit_games, cppgambit_core], ext_modules=Cython.Build.cythonize(libgambit, language_level="3str", diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index 099e4aa2a6..6442287340 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -621,6 +621,10 @@ cdef extern from "solvers/logit/logit.h": double getitem "operator[]"(int) except +IndexError +cdef extern from "solvers/hp/hp.h": + stdlist[c_MixedStrategyProfile[double]] HPStrategySolve(c_Game) except +RuntimeError + + cdef extern from "nash.h": stdlist[c_MixedBehaviorProfile[double]] LogitBehaviorSolveWrapper( c_Game, double, double, double diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index 00c74cf3b9..881bd99e15 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -368,3 +368,7 @@ def _logit_behavior_branch(game: Game, p.thisptr = profile_ptr ret.append(p) return ret + + +def _hp_strategy_solve(game: Game) -> list[MixedStrategyProfileDouble]: + return _convert_mspd(HPStrategySolve(game.game)) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 96201df0b6..a233448dfe 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -804,3 +804,32 @@ def logit_solve( equilibria=equilibria, parameters={"first_step": first_step, "max_accel": max_accel}, ) + + +def hp_solve( + game: libgbt.Game, +) -> NashComputationResult: + """Compute Nash equilibria of a game using :cite:p:`HerPee01` + + Returns an approximation to the limiting point on the principal branch of + the homotopy path for the game. + + Parameters + ---------- + game : Game + The game to compute equilibria in. + + Returns + ------- + res : NashComputationResult + The result represented as a ``NashComputationResult`` object. + """ + equilibria = libgbt._hp_strategy_solve(game) + return NashComputationResult( + game=game, + method="hp", + rational=False, + use_strategic=True, + equilibria=equilibria, + parameters={}, + ) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc new file mode 100644 index 0000000000..15e57b6e33 --- /dev/null +++ b/src/solvers/hp/hp.cc @@ -0,0 +1,38 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (https://www.gambit-project.org) +// +// FILE: src/solvers/hp/hp.cc +// Computation of a Nash equilibria using a differentiable homotopy +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// + +#include +#include "gambit.h" +#include "solvers/hp/hp.h" + +namespace Gambit { +std::list> HPStrategySolve(const Game &p_game) +{ + std::list> result; + const StrategySupportProfile support(p_game); + + const MixedStrategyProfile trivial_profile = support.NewMixedStrategyProfile(); + + result.push_back(trivial_profile); + return result; +} +} // namespace Gambit diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h new file mode 100644 index 0000000000..d51a0dc6c4 --- /dev/null +++ b/src/solvers/hp/hp.h @@ -0,0 +1,32 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (http://www.gambit-project.org) +// +// FILE: src/solvers/hp/hp.h +// Computation of a Nash equilibria using a differentiable homotopy +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// + +#ifndef HP_H +#define HP_H + +#include + +namespace Gambit { +std::list> HPStrategySolve(const Game &p_game); +} // namespace Gambit + +#endif // HP_H From 76fd80983d55f3cc5f3346757f1540b9c56f987a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Mon, 20 Jul 2026 20:42:00 +0200 Subject: [PATCH 2/9] Structure of hp method, ComputeInitialPoint and ExtractEquilibrium methods (#973) Implementation of algorithm equation system and termination condition at t=1. Passes several simple tests; more robust test suite to come. --- Makefile.am | 4 +- src/pygambit/gambit.pxd | 4 +- src/pygambit/nash.pxi | 5 +- src/pygambit/nash.py | 10 +- src/solvers/hp/hp.cc | 46 ++++++- src/solvers/hp/hp.h | 3 +- src/solvers/hp/hpsystem.cc | 239 +++++++++++++++++++++++++++++++++++++ src/solvers/hp/hpsystem.h | 70 +++++++++++ tests/test_hp.py | 132 ++++++++++++++++++++ 9 files changed, 497 insertions(+), 16 deletions(-) create mode 100644 src/solvers/hp/hpsystem.cc create mode 100644 src/solvers/hp/hpsystem.h create mode 100644 tests/test_hp.py diff --git a/Makefile.am b/Makefile.am index f2ea3fda3d..fde0f3ca11 100644 --- a/Makefile.am +++ b/Makefile.am @@ -315,7 +315,9 @@ gtracer_SOURCES = \ hp_SOURCES = \ src/solvers/hp/hp.h \ - src/solvers/hp/hp.cc + src/solvers/hp/hp.cc \ + src/solvers/hp/hpsystem.h \ + src/solvers/hp/hpsystem.cc if IS_WIN32 AM_LDFLAGS = -static -static-libgcc -static-libstdc++ diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index 6442287340..f147ad8ec1 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -622,7 +622,9 @@ cdef extern from "solvers/logit/logit.h": cdef extern from "solvers/hp/hp.h": - stdlist[c_MixedStrategyProfile[double]] HPStrategySolve(c_Game) except +RuntimeError + stdlist[c_MixedStrategyProfile[double]] HPStrategySolve( + c_MixedStrategyProfile[double] + ) except +RuntimeError cdef extern from "nash.h": diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index 881bd99e15..9dc1019323 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -370,5 +370,6 @@ def _logit_behavior_branch(game: Game, return ret -def _hp_strategy_solve(game: Game) -> list[MixedStrategyProfileDouble]: - return _convert_mspd(HPStrategySolve(game.game)) +def _hp_strategy_solve( + prior: MixedStrategyProfileDouble) -> list[MixedStrategyProfileDouble]: + return _convert_mspd(HPStrategySolve(deref(prior.profile))) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index a233448dfe..e1fdd0b09f 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -807,7 +807,7 @@ def logit_solve( def hp_solve( - game: libgbt.Game, + prior: libgbt.MixedStrategyProfileDouble, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -816,17 +816,17 @@ def hp_solve( Parameters ---------- - game : Game - The game to compute equilibria in. + prior : MixedStrategyProfileDouble + The prior distribution over strategies. Returns ------- res : NashComputationResult The result represented as a ``NashComputationResult`` object. """ - equilibria = libgbt._hp_strategy_solve(game) + equilibria = libgbt._hp_strategy_solve(prior) return NashComputationResult( - game=game, + game=prior.game, method="hp", rational=False, use_strategic=True, diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 15e57b6e33..9b71fdd07e 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -23,16 +23,50 @@ #include #include "gambit.h" #include "solvers/hp/hp.h" +#include "solvers/hp/hpsystem.h" +#include "solvers/logit/path.h" namespace Gambit { -std::list> HPStrategySolve(const Game &p_game) +std::list> +HPStrategySolve(const MixedStrategyProfile &p_prior) { - std::list> result; - const StrategySupportProfile support(p_game); - const MixedStrategyProfile trivial_profile = support.NewMixedStrategyProfile(); + std::list> equilibria; - result.push_back(trivial_profile); - return result; + HPEquationSystem system(p_prior); + Vector x = system.ComputeInitialPoint(); + + const PathTracer tracer; + double omega = 1.0; + + 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 TracePathResult result = tracer.TracePath( + [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, + [&system](const Vector &point, Matrix &jac) { + system.GetJacobian(point, jac); + }, + x, omega, termination_condition, + [&system](const Vector &point) { + std::cout << "[Path Tracer Step] t = " << point[1]; + std::cout << " | Alfas: "; + for (size_t i = 2; i <= 5; ++i) { + std::cout << point[i] << " "; + } + std::cout << "| Mu: " << point[6] << " " << point[7] << std::endl; + + std::cout << "Full point vector in probabilities: "; + Vector prob_vector = system.ExtractEquilibrium(point).GetProbVector(); + for (size_t i = 1; i <= prob_vector.size(); ++i) { + std::cout << prob_vector[i] << " "; + } + std::cout << std::endl; + }, + criterion_function); + + equilibria.push_back(system.ExtractEquilibrium(x)); + return equilibria; } } // namespace Gambit diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index d51a0dc6c4..2ed16a4303 100644 --- a/src/solvers/hp/hp.h +++ b/src/solvers/hp/hp.h @@ -26,7 +26,8 @@ #include namespace Gambit { -std::list> HPStrategySolve(const Game &p_game); +std::list> +HPStrategySolve(const MixedStrategyProfile &p_prior); } // namespace Gambit #endif // HP_H diff --git a/src/solvers/hp/hpsystem.cc b/src/solvers/hp/hpsystem.cc new file mode 100644 index 0000000000..1a5898722b --- /dev/null +++ b/src/solvers/hp/hpsystem.cc @@ -0,0 +1,239 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (https://www.gambit-project.org) +// +// FILE: src/solvers/hp/hpsystem.cc +// Computation of a Nash equilibria using a differentiable homotopy +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// + +#include +#include "gambit.h" +#include "solvers/hp/hpsystem.h" + +namespace Gambit { +HPEquationSystem::HPEquationSystem(const MixedStrategyProfile &prior) + : m_game(prior.GetGame()), m_prior(prior), + m_current_sigma(prior.GetGame()->NewMixedStrategyProfile(0.0)), + m_star(prior.MixedProfileLength()) +{ + m_payoffs_against_prior.reserve(m_prior.MixedProfileLength()); + for (const auto &player : m_game->GetPlayers()) { + for (const auto &strategy : player->GetStrategies()) { + m_payoffs_against_prior.push_back(m_prior.GetPayoff(strategy)); + } + } +} + +void HPEquationSystem::GetValue(const Vector &point, Vector &lhs) const +{ + const double t = point[1]; + + int temp_alpha_idx = 2; + for (const auto &player : m_game->GetPlayers()) { + for (const auto &strategy : player->GetStrategies()) { + m_current_sigma[strategy] = AlphaToSigma(point[temp_alpha_idx++]); + } + } + + int alpha_idx = 2; + int eq_idx = 1; + int player_idx = 1; + int flat_strategy_idx = 0; + + for (const auto &player : m_game->GetPlayers()) { + + const double mu = point[1 + m_star + player_idx]; + double sum_sigma = 0.0; + + // (a) Best response equations + for (const auto &strategy : player->GetStrategies()) { + const double alpha = point[alpha_idx++]; + const double lambda = AlphaToLambda(alpha); + const double sigma = AlphaToSigma(alpha); + sum_sigma += sigma; + + const double v_i = CalculateDynamicPayoff(flat_strategy_idx++, strategy, m_current_sigma, t); + + lhs[eq_idx++] = v_i + lambda - mu; + } + + // (b) Probability sum equation + lhs[eq_idx++] = sum_sigma - 1.0; + + player_idx++; + } +} + +void HPEquationSystem::GetJacobian(const Vector &point, Matrix &p_jac) const +{ + const double t = point[1]; + + // Compute current sigma from alpha values + int temp_alpha_idx = 2; + for (const auto &player : m_game->GetPlayers()) { + for (const auto &strategy : player->GetStrategies()) { + m_current_sigma[strategy] = AlphaToSigma(point[temp_alpha_idx++]); + } + } + + // Initialize the Jacobian matrix to zero + p_jac = 0.0; + + int eq_idx = 1; + int player_idx = 1; + int flat_s1_idx = 0; + + for (const auto &player1 : m_game->GetPlayers()) { + + // (a) Best response equations + for (const auto &strat1 : player1->GetStrategies()) { + + // Column of t: + const double payoff_vs_sigma = m_current_sigma.GetPayoff(strat1); + const double payoff_vs_prior = m_payoffs_against_prior[flat_s1_idx]; + p_jac(1, eq_idx) = payoff_vs_sigma - payoff_vs_prior; + + // Column of mu_i: Derivative with respect to mu of this player (-1.0) + p_jac(1 + m_star + player_idx, eq_idx) = -1.0; + + // Alpha columns: + int alpha_col = 2; + for (const auto &player2 : m_game->GetPlayers()) { + for (const auto &strat2 : player2->GetStrategies()) { + const double alpha2 = point[alpha_col]; + + if (player1 == player2) { + // Same strategy, use the derivative of lambda with respect to alpha + if (strat1 == strat2) { + p_jac(alpha_col, eq_idx) = AlphaToLambdaDeriv(alpha2); + } + } + else { + // Using chain rule + const double deriv_u = + m_current_sigma.GetPayoffDeriv(player1->GetNumber(), strat1, strat2); + const double deriv_sigma = AlphaToSigmaDeriv(alpha2); + p_jac(alpha_col, eq_idx) = t * deriv_u * deriv_sigma; + } + alpha_col++; + } + } + eq_idx++; + flat_s1_idx++; + } + + // (b) Probability sum equation + // Derivative with respect to 't' and 'mu' is 0 + // Only derivatives with respect to the alphas of this player + int alpha_col = 2; + for (const auto &player2 : m_game->GetPlayers()) { + for (const auto &strat2 : player2->GetStrategies()) { + if (player1 == player2) { + p_jac(alpha_col, eq_idx) = AlphaToSigmaDeriv(point[alpha_col]); + } + alpha_col++; + } + } + eq_idx++; + player_idx++; + } +} + +Vector HPEquationSystem::ComputeInitialPoint() const +{ + const int n_players = m_game->GetPlayers().size(); + const double tol = 1e-9; // Tolerance for floating-point comparisons + + // Dimension: 1 (t) + m_star (total strategies) + n (number of players) + const int vector_size = 1 + m_star + n_players; + Vector start_point(vector_size); + + start_point[1] = 0.0; // t = 0 + + int alpha_idx = 2; + int player_idx = 1; + int flat_strategy_idx = 0; // Index for accessing m_payoffs_against_prior + + for (const auto &player : m_game->GetPlayers()) { + const int temp_idx = flat_strategy_idx; + // Finding mu^i (the maximum payoff for player i against the prior) + double max_payoff = -std::numeric_limits::infinity(); + for (const auto &strategy : player->GetStrategies()) { + const double payoff = m_payoffs_against_prior[flat_strategy_idx++]; + if (payoff > max_payoff) { + max_payoff = payoff; + } + } + + // Store mu^i + start_point[1 + m_star + player_idx] = max_payoff; + + // Compute alpha^i_s for each strategy s of player i + + bool found_br = false; // Flag to check if a best response has been found + int local_s_idx = temp_idx; + for (const auto &strategy : player->GetStrategies()) { + const double lambda = max_payoff - m_payoffs_against_prior[local_s_idx++]; + if (std::abs(lambda) < tol && !found_br) { + start_point[alpha_idx++] = 1.0; // Best response + found_br = true; + } + else if (std::abs(lambda) < tol) { + throw std::runtime_error("Multiple best responses found for player " + + std::to_string(player_idx) + + ". Only one best response is allowed."); + } + else { + // Avoid sqrt of negative numbers + start_point[alpha_idx++] = -std::sqrt(std::max(0.0, lambda)); + } + } + player_idx++; + } + + return start_point; +} + +MixedStrategyProfile +HPEquationSystem::ExtractEquilibrium(const Vector &final_point) const +{ + MixedStrategyProfile ret = m_game->NewMixedStrategyProfile(0.0); + int alpha_idx = 2; // First position is reserved to t + + for (const auto &player : m_game->GetPlayers()) { + for (const auto &strategy : player->GetStrategies()) { + const double alpha_val = final_point[alpha_idx++]; + const double prob = this->AlphaToSigma(alpha_val); + ret[strategy] = prob; + } + } + ret = ret.Normalize(); + + return ret; +} + +// v^i(t, s) +double HPEquationSystem::CalculateDynamicPayoff(int action_index, const GameStrategy &strategy, + const MixedStrategyProfile ¤t_sigma, + double t) const +{ + const double payoff_against_sigma = current_sigma.GetPayoff(strategy); + const double payoff_against_prior = m_payoffs_against_prior[action_index]; + return t * payoff_against_sigma + (1.0 - t) * payoff_against_prior; +} + +} // end namespace Gambit diff --git a/src/solvers/hp/hpsystem.h b/src/solvers/hp/hpsystem.h new file mode 100644 index 0000000000..fd8d2e9cf1 --- /dev/null +++ b/src/solvers/hp/hpsystem.h @@ -0,0 +1,70 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (http://www.gambit-project.org) +// +// FILE: src/solvers/hp/hpsystem.h +// Computation of a Nash equilibria using a differentiable homotopy +// +// This program is free software; you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation; either version 2 of the License, or +// (at your option) any later version. +// +// This program is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with this program; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +// + +#ifndef HPSYSTEM_H +#define HPSYSTEM_H + +#include + +namespace Gambit { + +class HPEquationSystem { +public: + HPEquationSystem(const MixedStrategyProfile &prior); + + // Evaluates H(t, alpha, mu) = 0 + void GetValue(const Vector &point, Vector &lhs) const; + + void GetJacobian(const Vector &point, Matrix &jac) const; + + // Computes the initial point for the homotopy path tracing (t=0) + Vector ComputeInitialPoint() const; + + // Transforms the final vector into an equilibrium mixed strategy profile + MixedStrategyProfile ExtractEquilibrium(const Vector &final_point) const; + +private: + const Game m_game; + MixedStrategyProfile m_prior; + std::vector m_payoffs_against_prior; + int m_star; + mutable MixedStrategyProfile m_current_sigma; + + // Transforms alpha to sigma and lambda + inline double AlphaToSigma(double alpha) const { return (alpha > 0.0) ? (alpha * alpha) : 0.0; } + inline double AlphaToLambda(double alpha) const { return (alpha < 0.0) ? (alpha * alpha) : 0.0; } + + // d(sigma)/d(alpha) + inline double AlphaToSigmaDeriv(double alpha) const { return (alpha > 0.0) ? 2.0 * alpha : 0.0; } + // d(lambda)/d(alpha) + inline double AlphaToLambdaDeriv(double alpha) const + { + return (alpha < 0.0) ? 2.0 * alpha : 0.0; + } + + // v^i(t, s) + double CalculateDynamicPayoff(int action_index, const GameStrategy &strategy, + const MixedStrategyProfile ¤t_sigma, double t) const; +}; + +} // namespace Gambit +#endif // HPSYSTEM_H diff --git a/tests/test_hp.py b/tests/test_hp.py new file mode 100644 index 0000000000..4b6af3cca6 --- /dev/null +++ b/tests/test_hp.py @@ -0,0 +1,132 @@ +"""Test of calls to the Herings & Peeters (2001) homotopy solver.""" + +import dataclasses +import typing + +import numpy as np +import pytest + +import pygambit as gbt + +TOL = 1e-6 + + +def d(*probs) -> tuple: + """Helper function to let us write d() to be suggestive of + "probability distribution on simplex" ("Delta") + """ + return tuple(probs) + + +@dataclasses.dataclass +class HPSolverTestCase: + """Summarising the data relevant for a test fixture of a call to the HP solver.""" + factory: typing.Callable[[], gbt.MixedStrategyProfileDouble] + expected: list + prob_tol: float = TOL + + +def create_hs_base_game() -> gbt.Game: + """Creates the base 2x2 game used in all examples from Harsanyi & Selten (1988) Section 4.11 + and also featured in Herings & Peeters (2001). + """ + p1_payoffs = np.array([[2, 0], [0, 1]]) + p2_payoffs = np.array([[1, 0], [0, 4]]) + return gbt.Game.from_arrays(p1_payoffs, p2_payoffs, title="HS 1988 Base Game") + + +def create_hp_paper_example() -> gbt.MixedStrategyProfileDouble: + """Creates the example from Herings & Peeters (2001) Figure 1. + Also used in Harsanyi & Selten (1988) Section 4.11. -Second Example.""" + game = create_hs_base_game() + prior = game.mixed_strategy_profile() + p1, p2 = list(game.players) + + prior[list(p1.strategies)[0]] = 0.5 + prior[list(p1.strategies)[1]] = 0.5 + prior[list(p2.strategies)[0]] = 2.0 / 3.0 + prior[list(p2.strategies)[1]] = 1.0 / 3.0 + + return prior + + +def create_hs_example_1() -> gbt.MixedStrategyProfileDouble: + """Harsanyi & Selten (1988) Section 4.11 - First Example.""" + game = create_hs_base_game() + prior = game.mixed_strategy_profile() + p1, p2 = list(game.players) + + prior[list(p1.strategies)[0]] = 1.0 / 3.0 + prior[list(p1.strategies)[1]] = 2.0 / 3.0 + prior[list(p2.strategies)[0]] = 1.0 / 6.0 + prior[list(p2.strategies)[1]] = 5.0 / 6.0 + + return prior + + +def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: + """A prior that causes multiple best responses exactly at t=0.""" + game = create_hs_base_game() + prior = game.mixed_strategy_profile() + p1, p2 = list(game.players) + + prior[list(p1.strategies)[0]] = 2.0 / 3.0 + prior[list(p1.strategies)[1]] = 1.0 / 3.0 + prior[list(p2.strategies)[0]] = 1.0 / 3.0 + prior[list(p2.strategies)[1]] = 2.0 / 3.0 + + return prior + + +HP_CASES = [ + pytest.param( + HPSolverTestCase( + factory=create_hp_paper_example, + expected=[d(0.0, 1.0), d(0.0, 1.0)], + ), + id="test_hp_herings_peeters_example", + ), + pytest.param( + HPSolverTestCase( + factory=create_hs_example_1, + expected=[d(0.0, 1.0), d(0.0, 1.0)], + ), + id="test_hp_hs_example_1", + ), +] + + +@pytest.mark.nash +@pytest.mark.parametrize("test_case", HP_CASES) +def test_hp_strategy_solver(test_case: HPSolverTestCase, subtests) -> None: + """Test calls of the HP solver with starting priors. + + Subtests: + - Number of equilibria found is exactly 1. + - Equilibrium profile matches the expected theoretical result. + """ + prior = test_case.factory() + game = prior.game + + result = gbt.nash.hp_solve(prior=prior) + + with subtests.test("number of equilibria found"): + # The HP method uniquely selects exactly 1 equilibrium. + assert len(result.equilibria) == 1 + + eq = result.equilibria[0] + expected = game.mixed_strategy_profile(rational=False, data=test_case.expected) + + with subtests.test("strategy_profile matches expected"): + for player in game.players: + for strategy in player.strategies: + assert abs(eq[strategy] - expected[strategy]) <= test_case.prob_tol + + +@pytest.mark.nash +def test_hp_degenerate_t0_prior_raises_error() -> None: + """Test that the HP solver correctly identifies when given a degenerate prior.""" + prior = create_t0_degenerate_example() + with pytest.raises(RuntimeError, match="Multiple best responses found for player 1. " + "Only one best response is allowed."): + gbt.nash.hp_solve(prior=prior) From 72c965093dfc553ebe57cbae0801a337327f19e8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Mon, 3 Aug 2026 15:44:01 +0200 Subject: [PATCH 3/9] Refactoring HPsystem (#1019) Refactors the HP system of equations to the class-based equation pattern. --- src/solvers/hp/hpsystem.cc | 279 +++++++++++++++++++++++-------------- src/solvers/hp/hpsystem.h | 22 +-- 2 files changed, 181 insertions(+), 120 deletions(-) diff --git a/src/solvers/hp/hpsystem.cc b/src/solvers/hp/hpsystem.cc index 1a5898722b..3f590b0b87 100644 --- a/src/solvers/hp/hpsystem.cc +++ b/src/solvers/hp/hpsystem.cc @@ -25,131 +25,202 @@ #include "solvers/hp/hpsystem.h" namespace Gambit { -HPEquationSystem::HPEquationSystem(const MixedStrategyProfile &prior) - : m_game(prior.GetGame()), m_prior(prior), - m_current_sigma(prior.GetGame()->NewMixedStrategyProfile(0.0)), - m_star(prior.MixedProfileLength()) -{ - m_payoffs_against_prior.reserve(m_prior.MixedProfileLength()); - for (const auto &player : m_game->GetPlayers()) { - for (const auto &strategy : player->GetStrategies()) { - m_payoffs_against_prior.push_back(m_prior.GetPayoff(strategy)); + +class HPEquation { +public: + virtual ~HPEquation() = default; + + virtual double Value(const Vector &point, + const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior) const = 0; + + virtual void Gradient(const Vector &point, + const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior, + Vector &gradient) const = 0; +}; + +namespace { + +// Transforms alpha to sigma and lambda +inline double AlphaToSigma(double alpha) { return (alpha > 0.0) ? (alpha * alpha) : 0.0; } +inline double AlphaToLambda(double alpha) { return (alpha < 0.0) ? (alpha * alpha) : 0.0; } + +// d(sigma)/d(alpha) +inline double AlphaToSigmaDeriv(double alpha) { return (alpha > 0.0) ? 2.0 * alpha : 0.0; } +// d(lambda)/d(alpha) +inline double AlphaToLambdaDeriv(double alpha) { return (alpha < 0.0) ? 2.0 * alpha : 0.0; } + +// Eq (a): Best Response Equation +class BestResponseEquation final : public HPEquation { + GameStrategy m_strategy; + int m_alpha_idx; + int m_mu_idx; + int m_flat_s_idx; + +public: + BestResponseEquation(const GameStrategy &strat, int alpha_idx, int mu_idx, int flat_s_idx) + : m_strategy(strat), m_alpha_idx(alpha_idx), m_mu_idx(mu_idx), m_flat_s_idx(flat_s_idx) + { + } + + ~BestResponseEquation() override = default; + + double Value(const Vector &point, const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior) const override + { + const double t = point[1]; + const double alpha = point[m_alpha_idx]; + const double lambda = AlphaToLambda(alpha); + const double mu = point[m_mu_idx]; + + // Calculate Dynamic Payoff + const double payoff_vs_sigma = current_sigma.GetPayoff(m_strategy); + const double payoff_vs_prior = payoffs_against_prior[m_flat_s_idx]; + + // v^i(t, s) + const double v_i = t * payoff_vs_sigma + (1.0 - t) * payoff_vs_prior; + + return v_i + lambda - mu; + } + + void Gradient(const Vector &point, const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior, + Vector &gradient) const override + { + gradient = 0.0; + const double t = point[1]; + const GamePlayer my_player = m_strategy->GetPlayer(); + const Game game = m_strategy->GetGame(); + + // Derivative wrt t: + const double payoff_vs_sigma = current_sigma.GetPayoff(m_strategy); + const double payoff_vs_prior = payoffs_against_prior[m_flat_s_idx]; + gradient[1] = payoff_vs_sigma - payoff_vs_prior; + + // Derivative wrt mu_i: + gradient[m_mu_idx] = -1.0; + + // Derivatives wrt all alphas: + int alpha_col = 2; + for (const auto &player2 : game->GetPlayers()) { + for (const auto &strat2 : player2->GetStrategies()) { + const double alpha2 = point[alpha_col]; + + if (my_player == player2) { + if (m_strategy == strat2) { + gradient[alpha_col] = AlphaToLambdaDeriv(alpha2); + } + } + else { + // Chain rule for opposing players' strategies + const double deriv_u = + current_sigma.GetPayoffDeriv(my_player->GetNumber(), m_strategy, strat2); + const double deriv_sigma = AlphaToSigmaDeriv(alpha2); + gradient[alpha_col] = t * deriv_u * deriv_sigma; + } + alpha_col++; + } } } -} +}; -void HPEquationSystem::GetValue(const Vector &point, Vector &lhs) const -{ - const double t = point[1]; +// Eq (b): Probability Sum Equation +class ProbabilitySumEquation final : public HPEquation { + int m_first_alpha_idx; + int m_last_alpha_idx; - int temp_alpha_idx = 2; +public: + ProbabilitySumEquation(int first_idx, int last_idx) + : m_first_alpha_idx(first_idx), m_last_alpha_idx(last_idx) + { + } + + ~ProbabilitySumEquation() override = default; + + double Value(const Vector &point, const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior) const override + { + double sum_sigma = 0.0; + for (int i = m_first_alpha_idx; i < m_last_alpha_idx; ++i) { + sum_sigma += AlphaToSigma(point[i]); + } + return sum_sigma - 1.0; + } + + void Gradient(const Vector &point, const MixedStrategyProfile ¤t_sigma, + const std::vector &payoffs_against_prior, + Vector &gradient) const override + { + gradient = 0.0; + // Only non-zero derivatives are those with respect to the player's own alphas + for (int i = m_first_alpha_idx; i < m_last_alpha_idx; ++i) { + gradient[i] = AlphaToSigmaDeriv(point[i]); + } + } +}; + +} // end namespace + +HPEquationSystem::HPEquationSystem(const MixedStrategyProfile &prior) + : m_game(prior.GetGame()), m_prior(prior), m_star(prior.MixedProfileLength()), + m_current_sigma(prior.GetGame()->NewMixedStrategyProfile(0.0)) +{ + m_payoffs_against_prior.reserve(m_prior.MixedProfileLength()); for (const auto &player : m_game->GetPlayers()) { for (const auto &strategy : player->GetStrategies()) { - m_current_sigma[strategy] = AlphaToSigma(point[temp_alpha_idx++]); + m_payoffs_against_prior.push_back(m_prior.GetPayoff(strategy)); } } + // Pre-allocate space for all equations: m_star (Best Responses) + n (Prob Sums) + m_equations.reserve(m_star + m_game->GetPlayers().size()); int alpha_idx = 2; - int eq_idx = 1; int player_idx = 1; int flat_strategy_idx = 0; for (const auto &player : m_game->GetPlayers()) { + const int first_alpha = alpha_idx; + const int mu_idx = 1 + m_star + player_idx; - const double mu = point[1 + m_star + player_idx]; - double sum_sigma = 0.0; - - // (a) Best response equations + // Instantiate Best Response Equations for (const auto &strategy : player->GetStrategies()) { - const double alpha = point[alpha_idx++]; - const double lambda = AlphaToLambda(alpha); - const double sigma = AlphaToSigma(alpha); - sum_sigma += sigma; - - const double v_i = CalculateDynamicPayoff(flat_strategy_idx++, strategy, m_current_sigma, t); - - lhs[eq_idx++] = v_i + lambda - mu; + m_equations.push_back( + std::make_shared(strategy, alpha_idx, mu_idx, flat_strategy_idx)); + alpha_idx++; + flat_strategy_idx++; } - // (b) Probability sum equation - lhs[eq_idx++] = sum_sigma - 1.0; + // Instantiate Probability Sum Equations + m_equations.push_back(std::make_shared(first_alpha, alpha_idx)); player_idx++; } } -void HPEquationSystem::GetJacobian(const Vector &point, Matrix &p_jac) const +void HPEquationSystem::GetValue(const Vector &point, Vector &lhs) const { - const double t = point[1]; + // Update internal mutable state + UpdateSigma(point); - // Compute current sigma from alpha values - int temp_alpha_idx = 2; - for (const auto &player : m_game->GetPlayers()) { - for (const auto &strategy : player->GetStrategies()) { - m_current_sigma[strategy] = AlphaToSigma(point[temp_alpha_idx++]); - } + // Evaluate all equations + for (size_t i = 1; i <= m_equations.size(); ++i) { + lhs[i] = m_equations[i - 1]->Value(point, m_current_sigma, m_payoffs_against_prior); } +} - // Initialize the Jacobian matrix to zero - p_jac = 0.0; - - int eq_idx = 1; - int player_idx = 1; - int flat_s1_idx = 0; - - for (const auto &player1 : m_game->GetPlayers()) { - - // (a) Best response equations - for (const auto &strat1 : player1->GetStrategies()) { - - // Column of t: - const double payoff_vs_sigma = m_current_sigma.GetPayoff(strat1); - const double payoff_vs_prior = m_payoffs_against_prior[flat_s1_idx]; - p_jac(1, eq_idx) = payoff_vs_sigma - payoff_vs_prior; - - // Column of mu_i: Derivative with respect to mu of this player (-1.0) - p_jac(1 + m_star + player_idx, eq_idx) = -1.0; - - // Alpha columns: - int alpha_col = 2; - for (const auto &player2 : m_game->GetPlayers()) { - for (const auto &strat2 : player2->GetStrategies()) { - const double alpha2 = point[alpha_col]; +void HPEquationSystem::GetJacobian(const Vector &point, Matrix &p_jac) const +{ + // Update internal mutable state + UpdateSigma(point); - if (player1 == player2) { - // Same strategy, use the derivative of lambda with respect to alpha - if (strat1 == strat2) { - p_jac(alpha_col, eq_idx) = AlphaToLambdaDeriv(alpha2); - } - } - else { - // Using chain rule - const double deriv_u = - m_current_sigma.GetPayoffDeriv(player1->GetNumber(), strat1, strat2); - const double deriv_sigma = AlphaToSigmaDeriv(alpha2); - p_jac(alpha_col, eq_idx) = t * deriv_u * deriv_sigma; - } - alpha_col++; - } - } - eq_idx++; - flat_s1_idx++; - } + p_jac = 0.0; + Vector column(point.size()); // Temp vector matching Jacobian column size - // (b) Probability sum equation - // Derivative with respect to 't' and 'mu' is 0 - // Only derivatives with respect to the alphas of this player - int alpha_col = 2; - for (const auto &player2 : m_game->GetPlayers()) { - for (const auto &strat2 : player2->GetStrategies()) { - if (player1 == player2) { - p_jac(alpha_col, eq_idx) = AlphaToSigmaDeriv(point[alpha_col]); - } - alpha_col++; - } - } - eq_idx++; - player_idx++; + // Compute the Jacobian + for (size_t i = 1; i <= m_equations.size(); ++i) { + m_equations[i - 1]->Gradient(point, m_current_sigma, m_payoffs_against_prior, column); + p_jac.SetColumn(i, column); } } @@ -217,7 +288,7 @@ HPEquationSystem::ExtractEquilibrium(const Vector &final_point) const for (const auto &player : m_game->GetPlayers()) { for (const auto &strategy : player->GetStrategies()) { const double alpha_val = final_point[alpha_idx++]; - const double prob = this->AlphaToSigma(alpha_val); + const double prob = AlphaToSigma(alpha_val); ret[strategy] = prob; } } @@ -226,14 +297,14 @@ HPEquationSystem::ExtractEquilibrium(const Vector &final_point) const return ret; } -// v^i(t, s) -double HPEquationSystem::CalculateDynamicPayoff(int action_index, const GameStrategy &strategy, - const MixedStrategyProfile ¤t_sigma, - double t) const +void HPEquationSystem::UpdateSigma(const Vector &point) const { - const double payoff_against_sigma = current_sigma.GetPayoff(strategy); - const double payoff_against_prior = m_payoffs_against_prior[action_index]; - return t * payoff_against_sigma + (1.0 - t) * payoff_against_prior; + int temp_alpha_idx = 2; + for (const auto &player : m_game->GetPlayers()) { + for (const auto &strategy : player->GetStrategies()) { + m_current_sigma[strategy] = AlphaToSigma(point[temp_alpha_idx++]); + } + } } } // end namespace Gambit diff --git a/src/solvers/hp/hpsystem.h b/src/solvers/hp/hpsystem.h index fd8d2e9cf1..3f3db28ba9 100644 --- a/src/solvers/hp/hpsystem.h +++ b/src/solvers/hp/hpsystem.h @@ -27,9 +27,12 @@ namespace Gambit { +class HPEquation; + class HPEquationSystem { public: - HPEquationSystem(const MixedStrategyProfile &prior); + explicit HPEquationSystem(const MixedStrategyProfile &prior); + ~HPEquationSystem() = default; // Evaluates H(t, alpha, mu) = 0 void GetValue(const Vector &point, Vector &lhs) const; @@ -48,22 +51,9 @@ class HPEquationSystem { std::vector m_payoffs_against_prior; int m_star; mutable MixedStrategyProfile m_current_sigma; + std::vector> m_equations; - // Transforms alpha to sigma and lambda - inline double AlphaToSigma(double alpha) const { return (alpha > 0.0) ? (alpha * alpha) : 0.0; } - inline double AlphaToLambda(double alpha) const { return (alpha < 0.0) ? (alpha * alpha) : 0.0; } - - // d(sigma)/d(alpha) - inline double AlphaToSigmaDeriv(double alpha) const { return (alpha > 0.0) ? 2.0 * alpha : 0.0; } - // d(lambda)/d(alpha) - inline double AlphaToLambdaDeriv(double alpha) const - { - return (alpha < 0.0) ? 2.0 * alpha : 0.0; - } - - // v^i(t, s) - double CalculateDynamicPayoff(int action_index, const GameStrategy &strategy, - const MixedStrategyProfile ¤t_sigma, double t) const; + void UpdateSigma(const Vector &point) const; }; } // namespace Gambit From a7d358c2e38e7fbf507ff814613335dea948da61 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= <153729439+AndresFerCervell@users.noreply.github.com> Date: Wed, 12 Aug 2026 17:44:05 +0200 Subject: [PATCH 4/9] Adding a controlled direction in path.cc (#1030) This updates `path.cc ` so that the tracer always follows the curve in the direction chosen by the programmer. It has been proven useful in the new `HP`algorithm to avoid a random direction in `t=0`. `logit` has also been affected by the modifications: the headers have been updated. The mechanism implemented is the following: -In `path.h` an enum type named **`TraceDirection`** has been declared. It can be `Positive` associated to the number 1, and `Negative`, (-1). Positive means that the tracer will move forward; negative backwards. -To control whether the tracer is headed into the right direction in the first iteration, it is required to know which is the index of the variable that we are following (`lambda`in `logit`, `t` in `HP`). Another parameter has been added to `TracePath` named **`tracking_index`** that stores the position of that variable in the vector `x`. Note that in `HP` it is the first position and in `logit` is the last one. The check that we are heading the right way is as discussed: if the first step goes as expected (equivalent to the tangent being positive), the internal omega does not change. Otherwise, it is multiplied by -1. Tracing returns with an error result if the tangent in the first step is 0 (up to numerical tolerance) --- src/pygambit/nash.h | 33 +++++++++++++++++------------ src/solvers/hp/hp.cc | 6 +++--- src/solvers/logit/efglogit.cc | 19 +++++++++-------- src/solvers/logit/logit.h | 26 ++++++++++++----------- src/solvers/logit/nfglogit.cc | 19 +++++++++-------- src/solvers/logit/path.cc | 40 ++++++++++++++++++++++++++--------- src/solvers/logit/path.h | 4 +++- src/tools/logit/logit.cc | 17 +++++++++------ 8 files changed, 99 insertions(+), 65 deletions(-) diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 4643e4e381..0707488efb 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -22,6 +22,7 @@ #include "gambit.h" #include "solvers/logit/logit.h" +#include "solvers/logit/path.h" using namespace std; using namespace Gambit; @@ -32,8 +33,8 @@ std::list> LogitBehaviorSolveWrapper(const Game &p_ double p_maxAccel) { std::list> ret; - ret.push_back(LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, 1.0, - p_firstStep, p_maxAccel) + ret.push_back(LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel) .back() .GetProfile()); return ret; @@ -43,16 +44,17 @@ inline std::list LogitBehaviorPrincipalBranchWrapper(const Game &p_game, double p_regret, double p_firstStep, double p_maxAccel) { - return LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, 1.0, p_firstStep, - p_maxAccel); + return LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel); } std::shared_ptr LogitBehaviorEstimateWrapper(std::shared_ptr> p_frequencies, bool p_stopAtLocal, double p_firstStep, double p_maxAccel) { - return make_shared(LogitBehaviorEstimate( - *p_frequencies, 1000000.0, 1.0, p_stopAtLocal, p_firstStep, p_maxAccel)); + return make_shared( + LogitBehaviorEstimate(*p_frequencies, 1000000.0, PathTracer::TraceDirection::Positive, + p_stopAtLocal, p_firstStep, p_maxAccel)); } std::list> @@ -61,7 +63,8 @@ LogitBehaviorAtLambdaWrapper(const Game &p_game, const std::list &p_targ { LogitQREMixedBehaviorProfile start(p_game); std::list> ret; - for (auto &qre : LogitBehaviorSolveLambda(start, p_targetLambda, 1.0, p_firstStep, p_maxAccel)) { + for (auto &qre : LogitBehaviorSolveLambda( + start, p_targetLambda, PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel)) { ret.push_back(std::make_shared(qre)); } return ret; @@ -73,8 +76,8 @@ std::list> LogitStrategySolveWrapper(const Game &p_ double p_maxAccel) { std::list> ret; - ret.push_back(LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, 1.0, - p_firstStep, p_maxAccel) + ret.push_back(LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel) .back() .GetProfile()); return ret; @@ -84,8 +87,8 @@ inline std::list LogitStrategyPrincipalBranchWrapper(const Game &p_game, double p_regret, double p_firstStep, double p_maxAccel) { - return LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, 1.0, p_firstStep, - p_maxAccel); + return LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel); } std::list> @@ -94,7 +97,8 @@ LogitStrategyAtLambdaWrapper(const Game &p_game, const std::list &p_targ { LogitQREMixedStrategyProfile start(p_game); std::list> ret; - for (auto &qre : LogitStrategySolveLambda(start, p_targetLambda, 1.0, p_firstStep, p_maxAccel)) { + for (auto &qre : LogitStrategySolveLambda( + start, p_targetLambda, PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel)) { ret.push_back(std::make_shared(qre)); } return ret; @@ -104,6 +108,7 @@ std::shared_ptr LogitStrategyEstimateWrapper(std::shared_ptr> p_frequencies, bool p_stopAtLocal, double p_firstStep, double p_maxAccel) { - return make_shared(LogitStrategyEstimate( - *p_frequencies, 1000000.0, 1.0, p_stopAtLocal, p_firstStep, p_maxAccel)); + return make_shared( + LogitStrategyEstimate(*p_frequencies, 1000000.0, PathTracer::TraceDirection::Positive, + p_stopAtLocal, p_firstStep, p_maxAccel)); } diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 9b71fdd07e..1e75d8dddb 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -37,8 +37,8 @@ HPStrategySolve(const MixedStrategyProfile &p_prior) Vector x = system.ComputeInitialPoint(); const PathTracer tracer; - double omega = 1.0; - + 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; }; @@ -48,7 +48,7 @@ HPStrategySolve(const MixedStrategyProfile &p_prior) [&system](const Vector &point, Matrix &jac) { system.GetJacobian(point, jac); }, - x, omega, termination_condition, + x, direction, tracking_index, termination_condition, [&system](const Vector &point) { std::cout << "[Path Tracer Step] t = " << point[1]; std::cout << " | Alfas: "; diff --git a/src/solvers/logit/efglogit.cc b/src/solvers/logit/efglogit.cc index cae2f769d6..280bafc89f 100644 --- a/src/solvers/logit/efglogit.cc +++ b/src/solvers/logit/efglogit.cc @@ -309,8 +309,8 @@ void EstimatorCallbackFunction::EvaluatePoint(const Vector &p_point) namespace Gambit { std::list -LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, double p_omega, - double p_firstStep, double p_maxAccel, +LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, MixedBehaviorObserverFunctionType p_observer) { if (p_start.size() == 0) { @@ -335,7 +335,7 @@ LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, + x, p_direction, x.size(), [game, p_regret](const Vector &p_point) { return RegretTerminationFunction(game, p_point, p_regret); }, @@ -345,9 +345,9 @@ LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, std::list LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - MixedBehaviorObserverFunctionType p_observer) + const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, + double p_maxAccel, MixedBehaviorObserverFunctionType p_observer) { if (p_start.size() == 0) { return {p_start}; @@ -369,7 +369,7 @@ LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, LambdaPositiveTerminationFunction, + x, p_direction, x.size(), LambdaPositiveTerminationFunction, [&callback](const Vector &p_point) -> void { callback.AppendPoint(p_point); }, [lam](const Vector &x, const Vector &) -> double { return x.back() - lam; @@ -381,7 +381,8 @@ LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, LogitQREMixedBehaviorProfile LogitBehaviorEstimate(const MixedBehaviorProfile &p_frequencies, double p_maxLambda, - double p_omega, double p_stopAtLocal, double p_firstStep, double p_maxAccel, + PathTracer::TraceDirection p_direction, double p_stopAtLocal, + double p_firstStep, double p_maxAccel, MixedBehaviorObserverFunctionType p_observer) { const LogitQREMixedBehaviorProfile start(p_frequencies.GetGame()); @@ -405,7 +406,7 @@ LogitBehaviorEstimate(const MixedBehaviorProfile &p_frequencies, double [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, + x, p_direction, x.size(), [p_maxLambda](const Vector &p_point) { return LambdaRangeTerminationFunction(p_point, 0, p_maxLambda); }, diff --git a/src/solvers/logit/logit.h b/src/solvers/logit/logit.h index 8d619d44e0..a2d19fcd45 100644 --- a/src/solvers/logit/logit.h +++ b/src/solvers/logit/logit.h @@ -24,6 +24,7 @@ #define SOLVERS_LOGIT_H #include +#include "solvers/logit/path.h" namespace Gambit { @@ -85,18 +86,19 @@ using MixedStrategyObserverFunctionType = inline void NullMixedStrategyObserver(const LogitQREMixedStrategyProfile &) {} std::list LogitStrategySolve( - const LogitQREMixedStrategyProfile &p_start, double p_regret, double p_omega, - double p_firstStep, double p_maxAccel, + const LogitQREMixedStrategyProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, const MixedStrategyObserverFunctionType &p_observer = NullMixedStrategyObserver); std::list LogitStrategySolveLambda( const LogitQREMixedStrategyProfile &p_start, const std::list &p_targetLambda, - double p_omega, double p_firstStep, double p_maxAccel, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, const MixedStrategyObserverFunctionType &p_observer = NullMixedStrategyObserver); LogitQREMixedStrategyProfile LogitStrategyEstimate(const MixedStrategyProfile &p_frequencies, double p_maxLambda, - double p_omega, double p_stopAtLocal, double p_firstStep, double p_maxAccel, + PathTracer::TraceDirection p_direction, double p_stopAtLocal, + double p_firstStep, double p_maxAccel, MixedStrategyObserverFunctionType p_observer = NullMixedStrategyObserver); using LogitQREMixedBehaviorProfile = LogitQRE>; @@ -107,19 +109,19 @@ using MixedBehaviorObserverFunctionType = inline void NullMixedBehaviorObserver(const LogitQREMixedBehaviorProfile &) {} std::list -LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, double p_omega, - double p_firstStep, double p_maxAccel, +LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, MixedBehaviorObserverFunctionType p_observer = NullMixedBehaviorObserver); -std::list -LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - MixedBehaviorObserverFunctionType p_observer = NullMixedBehaviorObserver); +std::list LogitBehaviorSolveLambda( + const LogitQREMixedBehaviorProfile &p_start, const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + MixedBehaviorObserverFunctionType p_observer = NullMixedBehaviorObserver); LogitQREMixedBehaviorProfile LogitBehaviorEstimate(const MixedBehaviorProfile &p_frequencies, double p_maxLambda, - double p_omega, double p_stopAtLocal, double p_firstStep, double p_maxAccel, + PathTracer::TraceDirection p_direction, double p_stopAtLocal, + double p_firstStep, double p_maxAccel, MixedBehaviorObserverFunctionType p_observer = NullMixedBehaviorObserver); } // namespace Gambit diff --git a/src/solvers/logit/nfglogit.cc b/src/solvers/logit/nfglogit.cc index e8e9f81679..224e9566de 100644 --- a/src/solvers/logit/nfglogit.cc +++ b/src/solvers/logit/nfglogit.cc @@ -347,8 +347,8 @@ void EstimatorCallbackFunction::EvaluatePoint(const Vector &p_point) } // namespace std::list -LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret, double p_omega, - double p_firstStep, double p_maxAccel, +LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, const MixedStrategyObserverFunctionType &p_observer) { if (p_start.size() == 0) { @@ -374,7 +374,7 @@ LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret, [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, + x, p_direction, x.size(), [p_start, p_regret](const Vector &p_point) { return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret); }, @@ -384,9 +384,9 @@ LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret, std::list LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - const MixedStrategyObserverFunctionType &p_observer) + const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, + double p_maxAccel, const MixedStrategyObserverFunctionType &p_observer) { if (p_start.size() == 0) { return {p_start}; @@ -408,7 +408,7 @@ LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start, [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, LambdaPositiveTerminationFunction, + x, p_direction, x.size(), LambdaPositiveTerminationFunction, [&callback](const Vector &p_point) -> void { callback.AppendPoint(p_point); }, [lam](const Vector &x, const Vector &) -> double { return x.back() - lam; @@ -420,7 +420,8 @@ LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start, LogitQREMixedStrategyProfile LogitStrategyEstimate(const MixedStrategyProfile &p_frequencies, double p_maxLambda, - double p_omega, double p_stopAtLocal, double p_firstStep, double p_maxAccel, + PathTracer::TraceDirection p_direction, double p_stopAtLocal, + double p_firstStep, double p_maxAccel, MixedStrategyObserverFunctionType p_observer) { const LogitQREMixedStrategyProfile start(p_frequencies.GetGame()); @@ -444,7 +445,7 @@ LogitStrategyEstimate(const MixedStrategyProfile &p_frequencies, double [&system](const Vector &p_point, Matrix &p_jac) { system.GetJacobian(p_point, p_jac); }, - x, p_omega, + x, p_direction, x.size(), [p_maxLambda](const Vector &p_point) { return LambdaRangeTerminationFunction(p_point, 0, p_maxLambda); }, diff --git a/src/solvers/logit/path.cc b/src/solvers/logit/path.cc index 267d8f9148..d881561668 100644 --- a/src/solvers/logit/path.cc +++ b/src/solvers/logit/path.cc @@ -38,16 +38,16 @@ inline double sqr(double x) { return x * x; } void Givens(Matrix &b, Matrix &q, double &c1, double &c2, int l1, int l2, int l3) { - if (fabs(c1) + fabs(c2) == 0.0) { + if (std::abs(c1) + std::abs(c2) == 0.0) { return; } double sn; - if (fabs(c2) >= fabs(c1)) { - sn = std::sqrt(1.0 + sqr(c1 / c2)) * fabs(c2); + if (std::abs(c2) >= std::abs(c1)) { + sn = std::sqrt(1.0 + sqr(c1 / c2)) * std::abs(c2); } else { - sn = std::sqrt(1.0 + sqr(c2 / c1)) * fabs(c1); + sn = std::sqrt(1.0 + sqr(c2 / c1)) * std::abs(c1); } const double s1 = c1 / sn; const double s2 = c2 / sn; @@ -128,8 +128,9 @@ void NewtonStep(Matrix &q, Matrix &b, Vector &u, Vector< TracePathResult PathTracer::TracePath(std::function &, Vector &)> p_function, std::function &, Matrix &)> p_jacobian, - Vector &x, double &p_omega, TerminationFunctionType p_terminate, - CallbackFunctionType p_callback, CriterionFunctionType p_criterion, + Vector &x, TraceDirection p_direction, size_t p_trackingIndex, + TerminationFunctionType p_terminate, CallbackFunctionType p_callback, + CriterionFunctionType p_criterion, CriterionBracketFunctionType p_criterionBracket) const { const double c_tol = 1.0e-4; // tolerance for corrector iteration @@ -145,6 +146,7 @@ PathTracer::TracePath(std::function &, Vector 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 Vector u(x.size()); // t is current tangent at x; newT is tangent at u, which is the next point. @@ -157,17 +159,35 @@ PathTracer::TracePath(std::function &, Vector QRDecomp(b, q); q.GetRow(q.NumRows(), t); p_callback(x); + 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."}; + } while (!p_terminate(x)) { bool accept = true; - if (fabs(h) <= c_hmin) { + if (std::abs(h) <= c_hmin) { return {x, false, "Stepsize fell below minimum threshold."}; } + if (first_step) { + if (std::abs(t[p_trackingIndex]) <= c_orientTol) { + return {x, false, "Initial tangent vector is orthogonal to path-following direction."}; + } + // Ensure that the tangent is oriented in the same direction as + // the path-following direction. + else if (t[p_trackingIndex] < -c_orientTol) { + omega *= -1.0; + } + first_step = false; + } + // Predictor step for (size_t k = 1; k <= x.size(); k++) { - u[k] = x[k] + h * p_omega * t[k]; + u[k] = x[k] + h * omega * t[k]; } double decel = 1.0 / m_maxDecel; // initialize deceleration factor @@ -226,7 +246,7 @@ PathTracer::TracePath(std::function &, Vector if (!accept) { h /= m_maxDecel; // PC not accepted; change stepsize and retry - if (fabs(h) <= c_hmin) { + if (std::abs(h) <= c_hmin) { return {x, false, "Stepsize fell below minimum threshold."}; } continue; @@ -251,7 +271,7 @@ PathTracer::TracePath(std::function &, Vector } else { // Standard steplength adaptation - h = fabs(h / decel); + h = std::abs(h / decel); } // PC step was successful; update and iterate diff --git a/src/solvers/logit/path.h b/src/solvers/logit/path.h index 655194e42a..8193c5764d 100644 --- a/src/solvers/logit/path.h +++ b/src/solvers/logit/path.h @@ -71,6 +71,7 @@ struct TracePathResult { // class PathTracer { public: + enum class TraceDirection { Positive = 1, Negative = -1 }; PathTracer() = default; virtual ~PathTracer() = default; @@ -83,7 +84,8 @@ class PathTracer { TracePathResult TracePath(std::function &, Vector &)> p_function, std::function &, Matrix &)> p_jacobian, - Vector &p_x, double &p_omega, TerminationFunctionType p_terminate, + Vector &p_x, TraceDirection p_direction, size_t p_trackingIndex, + TerminationFunctionType p_terminate, CallbackFunctionType p_callback = NullCallbackFunction, CriterionFunctionType p_criterion = NullCriterionFunction, CriterionBracketFunctionType p_criterionBracker = NullCriterionBracketFunction) const; diff --git a/src/tools/logit/logit.cc b/src/tools/logit/logit.cc index a10fbebd2d..01f1a93ef5 100644 --- a/src/tools/logit/logit.cc +++ b/src/tools/logit/logit.cc @@ -206,7 +206,8 @@ int main(int argc, char *argv[]) } }; auto result = - LogitStrategyEstimate(frequencies, maxLambda, 1.0, false, hStart, maxDecel, printer); + LogitStrategyEstimate(frequencies, maxLambda, PathTracer::TraceDirection::Positive, + false, hStart, maxDecel, printer); PrintProfile(std::cout, decimals, result); return 0; } @@ -219,14 +220,15 @@ int main(int argc, char *argv[]) }; const LogitQREMixedStrategyProfile start(game); if (!targetLambda.empty()) { - auto result = - LogitStrategySolveLambda(start, targetLambda, 1.0, hStart, maxDecel, printer); + auto result = LogitStrategySolveLambda( + start, targetLambda, PathTracer::TraceDirection::Positive, hStart, maxDecel, printer); for (auto &profile : result) { PrintProfile(std::cout, decimals, profile); } } else { - auto result = LogitStrategySolve(start, maxregret, 1.0, hStart, maxDecel, printer); + auto result = LogitStrategySolve(start, maxregret, PathTracer::TraceDirection::Positive, + hStart, maxDecel, printer); PrintProfile(std::cout, decimals, result.back(), true); } } @@ -238,14 +240,15 @@ int main(int argc, char *argv[]) }; const LogitQREMixedBehaviorProfile start(game); if (!targetLambda.empty()) { - auto result = - LogitBehaviorSolveLambda(start, targetLambda, 1.0, hStart, maxDecel, printer); + auto result = LogitBehaviorSolveLambda( + start, targetLambda, PathTracer::TraceDirection::Positive, hStart, maxDecel, printer); for (auto &profile : result) { PrintProfile(std::cout, decimals, profile); } } else { - auto result = LogitBehaviorSolve(start, maxregret, 1.0, hStart, maxDecel, printer); + auto result = LogitBehaviorSolve(start, maxregret, PathTracer::TraceDirection::Positive, + hStart, maxDecel, printer); PrintProfile(std::cout, decimals, result.back(), true); } } From 59525dde8db02cb24f34ed2429c247567923f9a8 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Tue, 1 Sep 2026 20:18:22 +0100 Subject: [PATCH 5/9] Instrument HP with callbacks and cooperative cancel --- src/pygambit/callback.h | 24 ++++++++++++++++++++++++ src/pygambit/gambit.pxd | 12 ++++++------ src/pygambit/nash.h | 8 ++++++++ src/pygambit/nash.pxi | 25 +++++++++++++++++++++++-- src/pygambit/nash.py | 10 +++++++++- src/solvers/hp/hp.cc | 35 +++++++++++++---------------------- src/solvers/hp/hp.h | 34 ++++++++++++++++++++++++++++------ 7 files changed, 111 insertions(+), 37 deletions(-) diff --git a/src/pygambit/callback.h b/src/pygambit/callback.h index 7982c86b34..bd99dbbaab 100644 --- a/src/pygambit/callback.h +++ b/src/pygambit/callback.h @@ -33,6 +33,7 @@ #include "core/rational.h" #include "solvers/enumpoly/enumpoly.h" #include "solvers/gnm/gnm.h" +#include "solvers/hp/hp.h" #include "solvers/ipa/ipa.h" #include "solvers/liap/liap.h" #include "solvers/logit/logit.h" @@ -63,6 +64,10 @@ InvokeLogitStrategyEventCallback(PyObject *p_callback, std::string InvokeLogitBehaviorEventCallback(PyObject *p_callback, std::shared_ptr p_qre); +std::string +InvokeHPStrategyEventCallback(PyObject *p_callback, + std::shared_ptr> p_profile, + double p_t); std::string InvokeGNMPerturbationEventCallback( PyObject *p_callback, std::shared_ptr> p_profile); std::string @@ -228,6 +233,25 @@ MakeLogitEventCallback(PyObject *p_callbac }; } +/// +/// Builds an HPEventCallbackType which, when invoked with a point traced +/// along the HP homotopy path, calls a Python callable with the mixed +/// strategy profile and homotopy parameter t, converted to the +/// corresponding pygambit types. A null callback (Python `None`) yields the +/// solver's own no-op default. +/// +inline Gambit::Nash::HPEventCallbackType MakeHPEventCallback(PyObject *p_callback) +{ + if (!p_callback || p_callback == Py_None) { + return Gambit::Nash::NullHPEventCallback; + } + return [p_callback](const Gambit::Nash::HPEvent &p_event) { + const auto &step = std::get(p_event); + Gambit::ThrowIfPythonError(InvokeHPStrategyEventCallback( + p_callback, std::make_shared>(step.profile), step.t)); + }; +} + /// /// Builds a Nash::GNMEventCallbackType which, when invoked, dispatches to /// whichever Invoke*EventCallback trampoline matches the alternative held by diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index c176d84d8b..a29b0f08e5 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -549,6 +549,8 @@ cdef extern from "callback.h": pass cppclass LogitEventCallbackType "Gambit::LogitEventCallbackType"[T]: pass + cppclass HPEventCallbackType "Gambit::Nash::HPEventCallbackType": + pass cppclass GNMEventCallbackType "Gambit::Nash::GNMEventCallbackType": pass cppclass LiapEventCallbackType "Gambit::Nash::LiapEventCallbackType"[T]: @@ -563,6 +565,7 @@ cdef extern from "callback.h": StrategyCallbackType[T] MakeStrategyCallback[T](object) BehaviorCallbackType[T] MakeBehaviorCallback[T](object) LogitEventCallbackType[T] MakeLogitEventCallback[T](object) + HPEventCallbackType MakeHPEventCallback(object) GNMEventCallbackType MakeGNMEventCallback(object) LiapEventCallbackType[T] MakeLiapEventCallback[T](object) SimpdivEventCallbackType MakeSimpdivEventCallback(object) @@ -663,12 +666,6 @@ cdef extern from "solvers/logit/logit.h": double getitem "operator[]"(int) except +IndexError -cdef extern from "solvers/hp/hp.h": - stdlist[c_MixedStrategyProfile[double]] HPStrategySolve( - c_MixedStrategyProfile[double] - ) except +RuntimeError - - cdef extern from "nash.h": pair[ stdlist[c_MixedStrategyProfile[T]], stdlist[stdlist[c_MixedStrategyProfile[T]]] @@ -703,3 +700,6 @@ cdef extern from "nash.h": shared_ptr[c_MixedStrategyProfile[double]], bool, double, double, LogitEventCallbackType[c_LogitQREMixedStrategyProfile] ) except + + stdlist[c_MixedStrategyProfile[double]] HPStrategySolveWrapper( + c_MixedStrategyProfile[double], HPEventCallbackType + ) except +RuntimeError diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 632f330f92..787fd60919 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -21,6 +21,7 @@ // #include "solvers/enummixed/enummixed.h" +#include "solvers/hp/hp.h" #include "solvers/logit/logit.h" #include "solvers/logit/path.h" @@ -140,3 +141,10 @@ LogitStrategyEstimateWrapper(std::shared_ptr> p_fre LogitStrategyEstimate(*p_frequencies, 1000000.0, PathTracer::TraceDirection::Positive, p_stopAtLocal, p_firstStep, p_maxAccel, p_onEvent)); } + +std::list> +HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, + Nash::HPEventCallbackType p_onEvent = Nash::NullHPEventCallback) +{ + return Nash::HPStrategySolve(p_prior, Nash::NullStrategyCallback, p_onEvent); +} diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index f8ebac266e..e319846a3e 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -63,6 +63,13 @@ class GNMTerminationEvent: message: str +@dataclasses.dataclass(frozen=True) +class HPStepEvent: + """Reports one point traced along the HP homotopy path.""" + profile: MixedStrategyProfileDouble + t: float + + @dataclasses.dataclass(frozen=True) class LiapStartEvent: """Reports the starting point of a :ref:`Lyapunov function minimization ` run.""" @@ -207,6 +214,16 @@ cdef public string InvokeLogitBehaviorEventCallback( return b"" +cdef public string InvokeHPStrategyEventCallback( + callback, profile: shared_ptr[c_MixedStrategyProfile[float]], t: float +): + try: + callback(HPStepEvent(profile=MixedStrategyProfileDouble.wrap(profile), t=t)) + except BaseException as e: + return f"{type(e).__name__}: {e}".encode("utf-8") + return b"" + + cdef public string InvokeGNMPerturbationEventCallback( callback, profile: shared_ptr[c_MixedStrategyProfile[float]] ): @@ -896,5 +913,9 @@ def _logit_behavior_branch(game: Game, def _hp_strategy_solve( - prior: MixedStrategyProfileDouble) -> list[MixedStrategyProfileDouble]: - return _convert_mspd(HPStrategySolve(deref(prior.profile))) + prior: MixedStrategyProfileDouble, + event_callback: object = None, +) -> list[MixedStrategyProfileDouble]: + return _convert_mspd(HPStrategySolveWrapper( + deref(prior.profile), MakeHPEventCallback(event_callback) + )) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 95b3be189f..32aa3aa07c 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -1113,6 +1113,7 @@ def logit_solve( def hp_solve( prior: libgbt.MixedStrategyProfileDouble, + event_callback: Callable[[libgbt.HPStepEvent], None] | None = None, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -1124,12 +1125,19 @@ def hp_solve( prior : MixedStrategyProfileDouble The prior distribution over strategies. + 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 + to the returned equilibrium. + + .. versionadded:: 17.0.0 + Returns ------- res : NashComputationResult The result represented as a ``NashComputationResult`` object. """ - equilibria = libgbt._hp_strategy_solve(prior) + equilibria = libgbt._hp_strategy_solve(prior, event_callback) return NashComputationResult( game=prior.game, method="hp", diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 1e75d8dddb..5f7c303e2e 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -20,17 +20,17 @@ // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // -#include #include "gambit.h" #include "solvers/hp/hp.h" #include "solvers/hp/hpsystem.h" #include "solvers/logit/path.h" -namespace Gambit { +namespace Gambit::Nash { std::list> -HPStrategySolve(const MixedStrategyProfile &p_prior) +HPStrategySolve(const MixedStrategyProfile &p_prior, + StrategyCallbackType p_onEquilibrium, HPEventCallbackType p_onEvent, + const CancelToken &p_cancel) { - std::list> equilibria; HPEquationSystem system(p_prior); @@ -43,30 +43,21 @@ HPStrategySolve(const MixedStrategyProfile &p_prior) auto criterion_function = [](const Vector &point, const Vector &tangent) -> double { return point[1] - 1.0; }; - const TracePathResult result = tracer.TracePath( + tracer.TracePath( [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, [&system](const Vector &point, Matrix &jac) { system.GetJacobian(point, jac); }, x, direction, tracking_index, termination_condition, - [&system](const Vector &point) { - std::cout << "[Path Tracer Step] t = " << point[1]; - std::cout << " | Alfas: "; - for (size_t i = 2; i <= 5; ++i) { - std::cout << point[i] << " "; - } - std::cout << "| Mu: " << point[6] << " " << point[7] << std::endl; - - std::cout << "Full point vector in probabilities: "; - Vector prob_vector = system.ExtractEquilibrium(point).GetProbVector(); - for (size_t i = 1; i <= prob_vector.size(); ++i) { - std::cout << prob_vector[i] << " "; - } - std::cout << std::endl; + [&system, &p_onEvent](const Vector &point) { + const MixedStrategyProfile profile = system.ExtractEquilibrium(point); + p_onEvent(HPStepEvent{.profile = profile, .t = point[1]}); }, - criterion_function); + criterion_function, NullCriterionBracketFunction, p_cancel); - equilibria.push_back(system.ExtractEquilibrium(x)); + const MixedStrategyProfile equilibrium = system.ExtractEquilibrium(x); + p_onEquilibrium(equilibrium); + equilibria.push_back(equilibrium); return equilibria; } -} // namespace Gambit +} // namespace Gambit::Nash diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index 2ed16a4303..4bc5a7280f 100644 --- a/src/solvers/hp/hp.h +++ b/src/solvers/hp/hp.h @@ -20,14 +20,36 @@ // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // -#ifndef HP_H -#define HP_H +#ifndef GAMBIT_SOLVERS_HP_HP_H +#define GAMBIT_SOLVERS_HP_HP_H +#include #include +#include -namespace Gambit { +#include "solvers/nash.h" + +namespace Gambit::Nash { + +/// @brief Reports a point traced along the HP homotopy path, at homotopy parameter \p t +struct HPStepEvent { + const MixedStrategyProfile &profile; + double t; +}; + +using HPEvent = std::variant; +using HPEventCallbackType = std::function; + +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); -} // namespace Gambit +HPStrategySolve(const MixedStrategyProfile &p_prior, + StrategyCallbackType p_onEquilibrium = NullStrategyCallback, + HPEventCallbackType p_onEvent = NullHPEventCallback, + const CancelToken &p_cancel = CancelToken()); + +} // namespace Gambit::Nash -#endif // HP_H +#endif // GAMBIT_SOLVERS_HP_HP_H From 4ba0e642926743d5fdf0502f79b1bc3a71eaf0d5 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 2 Sep 2026 07:18:36 +0100 Subject: [PATCH 6/9] Add command-line wrapper for HP --- doc/pygambit.api.rst | 1 + doc/tools.rst | 1 + pyproject.toml | 1 + src/pygambit/cli/hp.py | 118 +++++++++++++++++++++++++++++++++++++ tests/cli/test_contract.py | 2 + tests/cli/test_hp.py | 96 ++++++++++++++++++++++++++++++ 6 files changed, 219 insertions(+) create mode 100644 src/pygambit/cli/hp.py create mode 100644 tests/cli/test_hp.py diff --git a/doc/pygambit.api.rst b/doc/pygambit.api.rst index e8ef2c8a65..67b286597d 100644 --- a/doc/pygambit.api.rst +++ b/doc/pygambit.api.rst @@ -403,6 +403,7 @@ Computation of Nash equilibria simpdiv_solve ipa_solve gnm_solve + hp_solve Computation of quantal response equilibria diff --git a/doc/tools.rst b/doc/tools.rst index 2efab1a955..99445adf14 100644 --- a/doc/tools.rst +++ b/doc/tools.rst @@ -53,3 +53,4 @@ documentation. tools.logit tools.gnm tools.ipa + tools.hp diff --git a/pyproject.toml b/pyproject.toml index 9dd6042676..78ce2dedde 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -43,6 +43,7 @@ gambit-simpdiv = "pygambit.cli.simpdiv:main" gambit-gnm = "pygambit.cli.gnm:main" gambit-ipa = "pygambit.cli.ipa:main" gambit-enumpoly = "pygambit.cli.enumpoly:main" +gambit-hp = "pygambit.cli.hp:main" [project.optional-dependencies] test = ["pytest", "pytest-subtests", "nbformat", "nbclient", "ipykernel"] diff --git a/src/pygambit/cli/hp.py b/src/pygambit/cli/hp.py new file mode 100644 index 0000000000..0daa664b61 --- /dev/null +++ b/src/pygambit/cli/hp.py @@ -0,0 +1,118 @@ +# +# This file is part of Gambit +# Copyright (c) 1994-2026, The Gambit Project (https://www.gambit-project.org) +# +# FILE: src/pygambit/cli/hp.py +# Command-line driver program for Nash equilibrium computation via the +# Herings-Peeters (2001) homotopy method +# +# This program is free software; you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation; either version 2 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program; if not, write to the Free Software +# Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. +# +"""Command-line driver program for Nash equilibrium computation via the +Herings-Peeters (2001) homotopy method. +""" + +from __future__ import annotations + +import click + +import pygambit as gbt + +from .common import ( + handle_errors, + load_game, + render_profile_csv, + resolve_strategy_starts, + version_option, +) + +DESCRIPTION = "Compute a Nash equilibrium using the Herings-Peeters (2001) homotopy method" +PROG_NAME = "gambit-hp" + + +@click.command( + context_settings={"help_option_names": ["-h", "--help"]}, + help=( + f"{DESCRIPTION}.\n\n" + "Reads a game from FILE, or from standard input if FILE is not specified." + ), +) +@click.argument("file", required=False, default=None) +@click.option( + "-d", + "decimals", + default=6, + show_default=True, + type=int, + help="show equilibria as floating point with DECIMALS digits", +) +@click.option( + "-n", + "n_priors", + type=int, + default=None, + help="number of prior distributions to generate randomly (mutually exclusive with -s)", +) +@click.option( + "-R", + "seed", + type=int, + default=None, + help="seed the random number generator used to generate prior distributions " + "(default is to seed from system entropy); requires -n", +) +@click.option( + "-s", + "start_file", + type=str, + default=None, + help="file containing prior distributions (mutually exclusive with -n)", +) +@click.option("-q", "--quiet", is_flag=True, help="quiet mode (suppresses banner)") +@click.option( + "-V", + "--verbose", + is_flag=True, + help="verbose mode (shows each point traced along the homotopy path)", +) +@version_option(DESCRIPTION) +@handle_errors +def main( + file: str | None, + decimals: int, + n_priors: int | None, + seed: int | None, + start_file: str | None, + quiet: bool, + verbose: bool, +) -> None: + game = load_game(quiet, DESCRIPTION, file, PROG_NAME) + priors = resolve_strategy_starts(game, n_priors, seed, start_file) + + def render_event(event: gbt.HPStepEvent) -> None: + if verbose: + click.echo(render_profile_csv(event.profile, f"{event.t:.6g}", decimals)) + + for prior in priors: + prior = prior.as_float() + if verbose: + click.echo(render_profile_csv(prior, "prior", decimals)) + result = gbt.nash.hp_solve(prior, event_callback=render_event) + for eq in result.equilibria: + click.echo(render_profile_csv(eq, "NE", decimals)) + + +if __name__ == "__main__": + main() diff --git a/tests/cli/test_contract.py b/tests/cli/test_contract.py index f836ebc50c..d1d9b4f68f 100644 --- a/tests/cli/test_contract.py +++ b/tests/cli/test_contract.py @@ -12,6 +12,7 @@ enumpoly, enumpure, gnm, + hp, ipa, lcp, liap, @@ -31,6 +32,7 @@ gnm, ipa, enumpoly, + hp, ] diff --git a/tests/cli/test_hp.py b/tests/cli/test_hp.py new file mode 100644 index 0000000000..5fa4991024 --- /dev/null +++ b/tests/cli/test_hp.py @@ -0,0 +1,96 @@ +"""Tests that gambit-hp's switches produce the behavior they document.""" + +from pygambit.cli import hp + + +def _start_file(tmp_path, text="0.9,0.1,0.9,0.1\n"): + path = tmp_path / "prior.csv" + path.write_text(text) + return path + + +def test_default_reports_one_equilibrium_per_prior( + cli_runner, nfg_asymmetric_table_text, tmp_path +): + start_file = _start_file(tmp_path) + result = cli_runner.invoke( + hp.main, ["-q", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + assert result.exit_code == 0 + assert result.stdout.strip().splitlines() == ["NE,1.000000,0.000000,1.000000,0.000000"] + + +def test_starting_file_accepts_exact_fractions(cli_runner, nfg_asymmetric_table_text, tmp_path): + start_file = _start_file(tmp_path, "9/10,1/10,9/10,1/10\n") + result = cli_runner.invoke( + hp.main, ["-q", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + assert result.exit_code == 0 + assert result.stdout.strip().splitlines() == ["NE,1.000000,0.000000,1.000000,0.000000"] + + +def test_decimals_flag_changes_precision(cli_runner, nfg_asymmetric_table_text, tmp_path): + start_file = _start_file(tmp_path) + result = cli_runner.invoke( + hp.main, ["-q", "-d", "2", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + assert result.exit_code == 0 + assert result.stdout.strip().splitlines() == ["NE,1.00,0.00,1.00,0.00"] + + +def test_verbose_flag_adds_prior_and_step_lines(cli_runner, nfg_asymmetric_table_text, tmp_path): + start_file = _start_file(tmp_path) + plain = cli_runner.invoke( + hp.main, ["-q", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + verbose = cli_runner.invoke( + hp.main, ["-q", "-V", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + assert verbose.exit_code == 0 + lines = verbose.stdout.strip().splitlines() + # The first line reports the (generally mixed) prior itself; the second is the + # t=0 step, the best response to that prior, which is pure by construction. + assert lines[0] == "prior,0.900000,0.100000,0.900000,0.100000" + assert lines[1].startswith("0,") + assert len(lines[1].split(",")) == 5 + assert len(lines) > len(plain.stdout.strip().splitlines()) + for line in plain.stdout.strip().splitlines(): + assert line in verbose.stdout + + +def test_n_and_s_are_mutually_exclusive(cli_runner, nfg_asymmetric_table_text, tmp_path): + start_file = _start_file(tmp_path) + result = cli_runner.invoke( + hp.main, ["-q", "-n", "2", "-s", str(start_file)], input=nfg_asymmetric_table_text + ) + assert result.exit_code == 1 + assert result.stderr == "Error: The -n and -s options are mutually exclusive.\n" + + +def test_seed_requires_n(cli_runner, nfg_asymmetric_table_text): + result = cli_runner.invoke(hp.main, ["-q", "-R", "5"], input=nfg_asymmetric_table_text) + assert result.exit_code == 1 + assert result.stderr == "Error: The -R option requires -n.\n" + + +def test_n_controls_the_number_of_priors_reported(cli_runner, nfg_asymmetric_table_text): + one = cli_runner.invoke(hp.main, ["-q", "-n", "1", "-R", "1"], input=nfg_asymmetric_table_text) + three = cli_runner.invoke( + hp.main, ["-q", "-n", "3", "-R", "1"], input=nfg_asymmetric_table_text + ) + assert one.exit_code == 0 + assert three.exit_code == 0 + # Each prior contributes exactly one reported equilibrium. + assert len(one.stdout.strip().splitlines()) == 1 + assert len(three.stdout.strip().splitlines()) == 3 + + +def test_degenerate_prior_is_a_clean_error(cli_runner, nfg_coordination_text, tmp_path): + # A prior exactly tied between two best responses has no unique best response + # at t=0, which HPStrategySolve rejects rather than picking one arbitrarily. + start_file = _start_file(tmp_path, "0.5,0.5,0.5,0.5\n") + result = cli_runner.invoke( + hp.main, ["-q", "-s", str(start_file)], input=nfg_coordination_text + ) + assert result.exit_code == 1 + assert "Multiple best responses" in result.stderr From c1b99b2dcdae7faa010bddc8d54985f6110eb9cc Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 2 Sep 2026 07:41:56 +0100 Subject: [PATCH 7/9] Expose HP as an option in the GUI --- Makefile.am | 5 +++-- src/gui/dllogit.h | 8 ++++---- src/gui/dlnash.cc | 16 +++++++++++++++- src/gui/nashspec.cc | 20 ++++++++++++++++++-- src/gui/nashspec.h | 22 +++++++++++++++++----- 5 files changed, 57 insertions(+), 14 deletions(-) diff --git a/Makefile.am b/Makefile.am index 9abfdf7a55..ef1c4a6536 100644 --- a/Makefile.am +++ b/Makefile.am @@ -363,13 +363,14 @@ AM_CXXFLAGS = ${LLVM_CXXFLAGS} ## recompiling the same core/games/solver sources from scratch. noinst_LIBRARIES = libcore.a libgames.a libbimatrix.a libgtracer.a \ - libliap.a liblogit.a libsimpdiv.a libenumpoly.a + libliap.a liblogit.a libsimpdiv.a libenumpoly.a libhp.a libcore_a_SOURCES = ${core_SOURCES} libgames_a_SOURCES = ${game_SOURCES} libbimatrix_a_SOURCES = ${bimatrix_SOURCES} libgtracer_a_SOURCES = ${gtracerlib_SOURCES} libliap_a_SOURCES = ${liap_SOURCES} +libhp_a_SOURCES = ${hp_SOURCES} liblogit_a_SOURCES = ${logit_SOURCES} libsimpdiv_a_SOURCES = ${simpdiv_SOURCES} libenumpoly_a_SOURCES = ${enumpoly_SOURCES} @@ -450,7 +451,7 @@ gambit_CXXFLAGS = $(AM_CXXFLAGS) $(WX_CXXFLAGS) gambit_CPPFLAGS = $(AM_CPPFLAGS) $(WX_CXXFLAGS) gambit_LDADD_LIBS = libbimatrix.a libliap.a liblogit.a libgtracer.a \ - libsimpdiv.a libenumpoly.a libgames.a libcore.a + libsimpdiv.a libenumpoly.a libhp.a libgames.a libcore.a gambit_DEPENDENCIES = $(RC_OBJECT_PATH) $(gambit_LDADD_LIBS) diff --git a/src/gui/dllogit.h b/src/gui/dllogit.h index d6462b2cc0..c946e64198 100644 --- a/src/gui/dllogit.h +++ b/src/gui/dllogit.h @@ -87,8 +87,8 @@ struct BehavLogitTraits { const CancelToken &p_cancel) { const QREType start(p_game); - LogitBehaviorSolve(start, 1.0e-8, 1.0, 0.03, 1.1, Nash::NullBehaviorCallback, - p_onEvent, p_cancel); + LogitBehaviorSolve(start, 1.0e-8, PathTracer::TraceDirection::Positive, 0.03, 1.1, + Nash::NullBehaviorCallback, p_onEvent, p_cancel); } }; @@ -123,8 +123,8 @@ struct MixedLogitTraits { const CancelToken &p_cancel) { const QREType start(p_game); - LogitStrategySolve(start, 1.0e-8, 1.0, 0.03, 1.1, Nash::NullStrategyCallback, - p_onEvent, p_cancel); + LogitStrategySolve(start, 1.0e-8, PathTracer::TraceDirection::Positive, 0.03, 1.1, + Nash::NullStrategyCallback, p_onEvent, p_cancel); } }; diff --git a/src/gui/dlnash.cc b/src/gui/dlnash.cc index 47c2a47a0b..702b3b63d2 100644 --- a/src/gui/dlnash.cc +++ b/src/gui/dlnash.cc @@ -38,6 +38,7 @@ static wxString s_enumpure(wxT("by looking for pure strategy equilibria")); static wxString s_enummixed(wxT("by enumerating extreme points")); static wxString s_enumpoly(wxT("by solving systems of polynomial equations")); static wxString s_gnm(wxT("by global Newton tracing")); +static wxString s_hp(wxT("by the Herings-Peeters homotopy")); static wxString s_ipa(wxT("by iterated polymatrix approximation")); static wxString s_lp(wxT("by solving a linear program")); static wxString s_lcp(wxT("by solving a linear complementarity program")); @@ -91,6 +92,9 @@ NashMethodSpec ResolveMethod(const wxString &p_method, NashEquilibriumTarget p_t if (p_method == s_gnm) { return GNMNashSpec{}; } + if (p_method == s_hp) { + return HPNashSpec{}; + } if (p_method == s_ipa) { return IPANashSpec{}; } @@ -118,7 +122,7 @@ NashMethodSpec ResolveMethod(const wxString &p_method, NashEquilibriumTarget p_t template concept StrategicMethod = std::same_as || std::same_as || - std::same_as || std::same_as || + std::same_as || std::same_as || std::same_as || std::same_as || std::same_as; template @@ -170,6 +174,9 @@ wxString ExternalCommand(const NashComputationSpec &p_spec) method.localNewtonInterval, method.localNewtonMaxIterations); } + else if constexpr (std::is_same_v) { + return prefix + wxString::Format("hp -d 10 -n %d", method.priors); + } else if constexpr (std::is_same_v) { return prefix + wxString::Format("ipa -d 10 -n %d", method.perturbations); } @@ -223,6 +230,9 @@ wxString MethodDescription(const NashMethodSpec &p_method) else if constexpr (std::is_same_v) { return wxT("by global Newton tracing"); } + else if constexpr (std::is_same_v) { + return wxT("by the Herings-Peeters homotopy"); + } else if constexpr (std::is_same_v) { return wxT("by iterated polymatrix approximation"); } @@ -263,6 +273,9 @@ wxString ParameterDescription(const NashMethodSpec &p_method) method.perturbations, method.lambdaEnd, method.steps, method.localNewtonInterval, method.localNewtonMaxIterations); } + else if constexpr (std::is_same_v) { + return wxString::Format(" (%d random priors)", method.priors); + } else if constexpr (std::is_same_v) { return wxString::Format(" (%d perturbation)", method.perturbations); } @@ -387,6 +400,7 @@ void NashChoiceDialog::OnCount(wxCommandEvent &p_event) m_methodChoice->Append(s_liap); m_methodChoice->Append(s_gnm); m_methodChoice->Append(s_ipa); + m_methodChoice->Append(s_hp); m_methodChoice->Append(s_enumpoly); } else { diff --git a/src/gui/nashspec.cc b/src/gui/nashspec.cc index 4c8ee29a23..77e5d17f31 100644 --- a/src/gui/nashspec.cc +++ b/src/gui/nashspec.cc @@ -26,6 +26,7 @@ #include "solvers/enumpoly/enumpoly.h" #include "solvers/enumpure/enumpure.h" #include "solvers/gnm/gnm.h" +#include "solvers/hp/hp.h" #include "solvers/ipa/ipa.h" #include "solvers/lcp/lcp.h" #include "solvers/liap/liap.h" @@ -95,6 +96,21 @@ std::optional GNMNashSpec::MakeSolver(NashRepresentation) const }; } +std::optional HPNashSpec::MakeSolver(NashRepresentation) const +{ + const HPNashSpec spec = *this; + return [spec](const Game &p_game, const ProfileFoundCallback &p_callback, + const CancelToken &p_cancel) { + for (const auto &prior : NewRandomStrategyProfiles(p_game, spec.priors)) { + p_cancel.Check(); + Nash::HPStrategySolve( + prior, + [&p_callback](const MixedStrategyProfile &p) { p_callback(ComputedProfile(p)); }, + Nash::NullHPEventCallback, p_cancel); + } + }; +} + std::optional IPANashSpec::MakeSolver(NashRepresentation) const { const IPANashSpec spec = *this; @@ -178,7 +194,7 @@ std::optional LogitNashSpec::MakeSolver(NashRepresentation p_rep const CancelToken &p_cancel) { const LogitQREMixedBehaviorProfile start(p_game); LogitBehaviorSolve( - start, spec.maxRegret, spec.omega, spec.firstStep, spec.maxAcceleration, + start, spec.maxRegret, spec.direction, spec.firstStep, spec.maxAcceleration, [&p_callback](const MixedBehaviorProfile &p) { p_callback(ComputedProfile(p)); }, NullLogitEventCallback, p_cancel); }; @@ -187,7 +203,7 @@ std::optional LogitNashSpec::MakeSolver(NashRepresentation p_rep const CancelToken &p_cancel) { const LogitQREMixedStrategyProfile start(p_game); LogitStrategySolve( - start, spec.maxRegret, spec.omega, spec.firstStep, spec.maxAcceleration, + start, spec.maxRegret, spec.direction, spec.firstStep, spec.maxAcceleration, [&p_callback](const MixedStrategyProfile &p) { p_callback(ComputedProfile(p)); }, NullLogitEventCallback, p_cancel); }; diff --git a/src/gui/nashspec.h b/src/gui/nashspec.h index 50af9a87d3..51698a8ed2 100644 --- a/src/gui/nashspec.h +++ b/src/gui/nashspec.h @@ -27,8 +27,11 @@ #include #include "core/cancel.h" -#include "solvers/nash.h" +#include "core/matrix.h" #include "core/rational.h" +#include "core/vector.h" +#include "solvers/logit/path.h" +#include "solvers/nash.h" namespace Gambit::GUI { @@ -94,6 +97,15 @@ struct IPANashSpec { std::optional MakeSolver(NashRepresentation) const; }; +struct HPNashSpec { + // Unlike GNM/IPA, a single prior yields at most one equilibrium (no internal + // path-tracing can surface more), so this defaults high, like LiapNashSpec's + // startingPoints, rather than to 1. + int priors{10}; + + std::optional MakeSolver(NashRepresentation) const; +}; + struct LPNashSpec { std::optional MakeSolver(NashRepresentation p_representation) const; }; @@ -115,7 +127,7 @@ struct LiapNashSpec { struct LogitNashSpec { double maxRegret{1.0e-8}; - double omega{1.0}; + PathTracer::TraceDirection direction{PathTracer::TraceDirection::Positive}; double firstStep{0.03}; double maxAcceleration{1.1}; @@ -132,9 +144,9 @@ struct SimpdivNashSpec { std::optional MakeSolver(NashRepresentation) const; }; -using NashMethodSpec = - std::variant; +using NashMethodSpec = std::variant; struct NashComputationSpec { NashRepresentation representation; From ea8a55e525b330bee93bc0e7eedf1dcd93bbb5e9 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 2 Sep 2026 08:02:48 +0100 Subject: [PATCH 8/9] Reorganise path-following methods internally. --- Makefile.am | 20 +++++++++++++------- setup.py | 7 +++---- src/gui/nashspec.h | 2 +- src/pygambit/nash.h | 2 +- src/solvers/hp/hp.cc | 2 +- src/solvers/logit/efglogit.cc | 2 +- src/solvers/logit/logit.h | 2 +- src/solvers/logit/nfglogit.cc | 2 +- src/solvers/{logit => path}/path.cc | 2 +- src/solvers/{logit => path}/path.h | 2 +- 10 files changed, 24 insertions(+), 19 deletions(-) rename src/solvers/{logit => path}/path.cc (99%) rename src/solvers/{logit => path}/path.h (99%) diff --git a/Makefile.am b/Makefile.am index ef1c4a6536..44fe6c2d6f 100644 --- a/Makefile.am +++ b/Makefile.am @@ -312,15 +312,22 @@ liap_SOURCES = \ src/solvers/liap/nfgliap.cc \ src/solvers/liap/liap.h +path_SOURCES = \ + src/solvers/path/path.cc \ + src/solvers/path/path.h + logit_SOURCES = \ src/solvers/logit/logbehav.h \ src/solvers/logit/logbehav.imp \ - src/solvers/logit/path.cc \ - src/solvers/logit/path.h \ src/solvers/logit/logit.h \ src/solvers/logit/efglogit.cc \ src/solvers/logit/nfglogit.cc +homotopy_SOURCES = \ + ${path_SOURCES} \ + ${logit_SOURCES} \ + ${hp_SOURCES} + simpdiv_SOURCES = \ src/solvers/simpdiv/simpdiv.cc \ src/solvers/simpdiv/simpdiv.h @@ -363,15 +370,14 @@ AM_CXXFLAGS = ${LLVM_CXXFLAGS} ## recompiling the same core/games/solver sources from scratch. noinst_LIBRARIES = libcore.a libgames.a libbimatrix.a libgtracer.a \ - libliap.a liblogit.a libsimpdiv.a libenumpoly.a libhp.a + libliap.a libhomotopy.a libsimpdiv.a libenumpoly.a libcore_a_SOURCES = ${core_SOURCES} libgames_a_SOURCES = ${game_SOURCES} libbimatrix_a_SOURCES = ${bimatrix_SOURCES} libgtracer_a_SOURCES = ${gtracerlib_SOURCES} libliap_a_SOURCES = ${liap_SOURCES} -libhp_a_SOURCES = ${hp_SOURCES} -liblogit_a_SOURCES = ${logit_SOURCES} +libhomotopy_a_SOURCES = ${homotopy_SOURCES} libsimpdiv_a_SOURCES = ${simpdiv_SOURCES} libenumpoly_a_SOURCES = ${enumpoly_SOURCES} @@ -450,8 +456,8 @@ gambit_SOURCES = \ gambit_CXXFLAGS = $(AM_CXXFLAGS) $(WX_CXXFLAGS) gambit_CPPFLAGS = $(AM_CPPFLAGS) $(WX_CXXFLAGS) -gambit_LDADD_LIBS = libbimatrix.a libliap.a liblogit.a libgtracer.a \ - libsimpdiv.a libenumpoly.a libhp.a libgames.a libcore.a +gambit_LDADD_LIBS = libbimatrix.a libliap.a libhomotopy.a libgtracer.a \ + libsimpdiv.a libenumpoly.a libgames.a libcore.a gambit_DEPENDENCIES = $(RC_OBJECT_PATH) $(gambit_LDADD_LIBS) diff --git a/setup.py b/setup.py index 4a9e213469..35a70b8410 100644 --- a/setup.py +++ b/setup.py @@ -95,11 +95,10 @@ def run(self) -> None: cppgambit_bimatrix = solver_library_config("cppgambit_bimatrix", ["linalg", "lp", "lcp", "enummixed"]) cppgambit_liap = solver_library_config("cppgambit_liap", ["liap"]) -cppgambit_logit = solver_library_config("cppgambit_logit", ["logit"]) +cppgambit_homotopy = solver_library_config("cppgambit_homotopy", ["path", "logit", "hp"]) cppgambit_gtracer = solver_library_config("cppgambit_gtracer", ["gtracer", "ipa", "gnm"]) cppgambit_simpdiv = solver_library_config("cppgambit_simpdiv", ["simpdiv"]) cppgambit_enumpoly = solver_library_config("cppgambit_enumpoly", ["nashsupport", "enumpoly"]) -cppgambit_hp = solver_library_config("cppgambit_hp", ["hp"]) libgambit = setuptools.Extension( @@ -112,8 +111,8 @@ def run(self) -> None: setuptools.setup( cmdclass={"build_py": GambitBuildPy}, - libraries=[cppgambit_bimatrix, cppgambit_liap, cppgambit_logit, cppgambit_simpdiv, - cppgambit_gtracer, cppgambit_enumpoly, cppgambit_hp, + libraries=[cppgambit_bimatrix, cppgambit_liap, cppgambit_homotopy, cppgambit_simpdiv, + cppgambit_gtracer, cppgambit_enumpoly, cppgambit_games, cppgambit_core], ext_modules=Cython.Build.cythonize(libgambit, language_level="3str", diff --git a/src/gui/nashspec.h b/src/gui/nashspec.h index 51698a8ed2..de07b4d933 100644 --- a/src/gui/nashspec.h +++ b/src/gui/nashspec.h @@ -30,8 +30,8 @@ #include "core/matrix.h" #include "core/rational.h" #include "core/vector.h" -#include "solvers/logit/path.h" #include "solvers/nash.h" +#include "solvers/path/path.h" namespace Gambit::GUI { diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 787fd60919..75dbfb2e20 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -23,7 +23,7 @@ #include "solvers/enummixed/enummixed.h" #include "solvers/hp/hp.h" #include "solvers/logit/logit.h" -#include "solvers/logit/path.h" +#include "solvers/path/path.h" using namespace std; using namespace Gambit; diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 5f7c303e2e..25b06325c2 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -23,7 +23,7 @@ #include "gambit.h" #include "solvers/hp/hp.h" #include "solvers/hp/hpsystem.h" -#include "solvers/logit/path.h" +#include "solvers/path/path.h" namespace Gambit::Nash { std::list> diff --git a/src/solvers/logit/efglogit.cc b/src/solvers/logit/efglogit.cc index f82ee94899..093460e0f4 100644 --- a/src/solvers/logit/efglogit.cc +++ b/src/solvers/logit/efglogit.cc @@ -26,7 +26,7 @@ #include "games.h" #include "logit.h" #include "logbehav.imp" -#include "path.h" +#include "solvers/path/path.h" namespace { diff --git a/src/solvers/logit/logit.h b/src/solvers/logit/logit.h index eb24bce6ca..a8f0faded2 100644 --- a/src/solvers/logit/logit.h +++ b/src/solvers/logit/logit.h @@ -26,8 +26,8 @@ #include #include -#include "solvers/logit/path.h" #include "solvers/nash.h" +#include "solvers/path/path.h" namespace Gambit { diff --git a/src/solvers/logit/nfglogit.cc b/src/solvers/logit/nfglogit.cc index e473aa0b1b..56fc3cdf08 100644 --- a/src/solvers/logit/nfglogit.cc +++ b/src/solvers/logit/nfglogit.cc @@ -25,7 +25,7 @@ #include "games.h" #include "logit.h" -#include "path.h" +#include "solvers/path/path.h" namespace Gambit { diff --git a/src/solvers/logit/path.cc b/src/solvers/path/path.cc similarity index 99% rename from src/solvers/logit/path.cc rename to src/solvers/path/path.cc index 01d0cf9390..0a6fcba301 100644 --- a/src/solvers/logit/path.cc +++ b/src/solvers/path/path.cc @@ -2,7 +2,7 @@ // This file is part of Gambit // Copyright (c) 1994-2026, The Gambit Project (https://www.gambit-project.org) // -// FILE: src/solvers/logit/path.cc +// FILE: src/solvers/path/path.cc // Implementation of generic smooth path-following algorithm. // // This program is free software; you can redistribute it and/or modify diff --git a/src/solvers/logit/path.h b/src/solvers/path/path.h similarity index 99% rename from src/solvers/logit/path.h rename to src/solvers/path/path.h index 8a7a935267..d05b1c9adf 100644 --- a/src/solvers/logit/path.h +++ b/src/solvers/path/path.h @@ -2,7 +2,7 @@ // This file is part of Gambit // Copyright (c) 1994-2026, The Gambit Project (https://www.gambit-project.org) // -// FILE: src/solvers/logit/path.h +// FILE: src/solvers/path/path.h // Interface to generic smooth path-following algorithm. // // This program is free software; you can redistribute it and/or modify From 7f1b71ea3ba57632e3d2e600ac7e7db9aa4d4ec2 Mon Sep 17 00:00:00 2001 From: Theodore Turocy Date: Wed, 2 Sep 2026 09:04:31 +0100 Subject: [PATCH 9/9] Integrate HP tests into main Nash tests --- pyproject.toml | 1 + tests/games.py | 9 ++++ tests/test_hp.py | 129 --------------------------------------------- tests/test_nash.py | 48 +++++++++++++++++ 4 files changed, 58 insertions(+), 129 deletions(-) delete mode 100644 tests/test_hp.py diff --git a/pyproject.toml b/pyproject.toml index 78ce2dedde..b3c9e30c48 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -125,6 +125,7 @@ markers = [ "nash_logit_behavior: tests of logit_solve in behavior strategies", "nash_gnm_strategy: tests of gnm_solve in mixed strategies", "nash_ipa_strategy: tests of lpa_solve in mixed strategies", + "nash_hp_strategy: tests of hp_solve in mixed strategies", "nash_simpdiv: tests of simpdiv_solve (in mixed strategies)", "nash_liap_strategy: tests of liap_solve (in mixed strategies)", "nash_liap_agent: tests of liap_agent_solve (in mixed behaviors)", diff --git a/tests/games.py b/tests/games.py index 07fc521c73..7e2ee54f7f 100644 --- a/tests/games.py +++ b/tests/games.py @@ -112,6 +112,15 @@ def create_efg_corresponding_to_bimatrix_game(g: gbt.Game) -> gbt.Game: return create_efg_corresponding_to_bimatrix_game_arrays(A, B, g.title) +def create_hs1988_base_game() -> gbt.Game: + """The base 2x2 game used in all examples from Harsanyi & Selten (1988) Section 4.11, + also featured as Figure 1 of Herings & Peeters (2001). + """ + p1_payoffs = np.array([[2, 0], [0, 1]]) + p2_payoffs = np.array([[1, 0], [0, 4]]) + return gbt.Game.from_arrays(p1_payoffs, p2_payoffs, title="HS 1988 Base Game") + + ################################################################################################ # Extensive-form games (efg) diff --git a/tests/test_hp.py b/tests/test_hp.py deleted file mode 100644 index 9c5bcf1508..0000000000 --- a/tests/test_hp.py +++ /dev/null @@ -1,129 +0,0 @@ -"""Test of calls to the Herings & Peeters (2001) homotopy solver.""" - -import dataclasses -import typing - -import numpy as np -import pytest - -import pygambit as gbt - -TOL = 1e-6 - - -def d(*probs) -> tuple: - """Helper function to let us write d() to be suggestive of - "probability distribution on simplex" ("Delta") - """ - return tuple(probs) - - -@dataclasses.dataclass -class HPSolverTestCase: - """Summarising the data relevant for a test fixture of a call to the HP solver.""" - factory: typing.Callable[[], gbt.MixedStrategyProfileDouble] - expected: list - prob_tol: float = TOL - - -def create_hs_base_game() -> gbt.Game: - """Creates the base 2x2 game used in all examples from Harsanyi & Selten (1988) Section 4.11 - and also featured in Herings & Peeters (2001). - """ - p1_payoffs = np.array([[2, 0], [0, 1]]) - p2_payoffs = np.array([[1, 0], [0, 4]]) - return gbt.Game.from_arrays(p1_payoffs, p2_payoffs, title="HS 1988 Base Game") - - -def create_hp_paper_example() -> gbt.MixedStrategyProfileDouble: - """Creates the example from Herings & Peeters (2001) Figure 1. - Also used in Harsanyi & Selten (1988) Section 4.11. -Second Example.""" - game = create_hs_base_game() - prior = game.mixed_strategy_profile() - p1, p2 = list(game.players) - s1, s2 = list(game.get_strategies(p1)), list(game.get_strategies(p2)) - - prior[p1] = {s1[0]: 0.5, s1[1]: 0.5} - prior[p2] = {s2[0]: 2.0 / 3.0, s2[1]: 1.0 / 3.0} - - return prior - - -def create_hs_example_1() -> gbt.MixedStrategyProfileDouble: - """Harsanyi & Selten (1988) Section 4.11 - First Example.""" - game = create_hs_base_game() - prior = game.mixed_strategy_profile() - p1, p2 = list(game.players) - s1, s2 = list(game.get_strategies(p1)), list(game.get_strategies(p2)) - - prior[p1] = {s1[0]: 1.0 / 3.0, s1[1]: 2.0 / 3.0} - prior[p2] = {s2[0]: 1.0 / 6.0, s2[1]: 5.0 / 6.0} - - return prior - - -def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: - """A prior that causes multiple best responses exactly at t=0.""" - game = create_hs_base_game() - prior = game.mixed_strategy_profile() - p1, p2 = list(game.players) - s1, s2 = list(game.get_strategies(p1)), list(game.get_strategies(p2)) - - prior[p1] = {s1[0]: 2.0 / 3.0, s1[1]: 1.0 / 3.0} - prior[p2] = {s2[0]: 1.0 / 3.0, s2[1]: 2.0 / 3.0} - - return prior - - -HP_CASES = [ - pytest.param( - HPSolverTestCase( - factory=create_hp_paper_example, - expected=[d(0.0, 1.0), d(0.0, 1.0)], - ), - id="test_hp_herings_peeters_example", - ), - pytest.param( - HPSolverTestCase( - factory=create_hs_example_1, - expected=[d(0.0, 1.0), d(0.0, 1.0)], - ), - id="test_hp_hs_example_1", - ), -] - - -@pytest.mark.nash -@pytest.mark.parametrize("test_case", HP_CASES) -def test_hp_strategy_solver(test_case: HPSolverTestCase, subtests) -> None: - """Test calls of the HP solver with starting priors. - - Subtests: - - Number of equilibria found is exactly 1. - - Equilibrium profile matches the expected theoretical result. - """ - prior = test_case.factory() - game = prior.game - - result = gbt.nash.hp_solve(prior=prior) - - with subtests.test("number of equilibria found"): - # The HP method uniquely selects exactly 1 equilibrium. - assert len(result.equilibria) == 1 - - eq = result.equilibria[0] - expected = game.mixed_strategy_profile(rational=False, data=test_case.expected) - - with subtests.test("strategy_profile matches expected"): - for player in game.players: - for strategy in game.get_strategies(player): - assert abs(eq[player][strategy] - expected[player][strategy]) <= test_case.prob_tol - - -@pytest.mark.nash -def test_hp_degenerate_t0_prior_raises_error() -> None: - """Test that the HP solver correctly identifies when given a degenerate prior.""" - prior = create_t0_degenerate_example() - with pytest.raises(RuntimeError, match="Multiple best responses found for player 1. " - "Only one best response is allowed."): - gbt.nash.hp_solve(prior=prior) diff --git a/tests/test_nash.py b/tests/test_nash.py index d49b599686..80030858e4 100644 --- a/tests/test_nash.py +++ b/tests/test_nash.py @@ -1696,6 +1696,36 @@ def _one_hot_perturbation(rational: bool) -> gbt.MixedStrategyProfile: ] +HP_STRATEGY_CASES = [ + pytest.param( + EquilibriumTestCaseWithStart( + factory=games.create_hs1988_base_game, + solver=gbt.nash.hp_solve, + start_data=dict(data=[[0.5, 0.5], [2.0 / 3.0, 1.0 / 3.0]], rational=False), + expected=[[d(0.0, 1.0), d(0.0, 1.0)]], + regret_tol=TOL_LARGE, + prob_tol=TOL_LARGE, + ), + marks=pytest.mark.nash_hp_strategy, + id="test_hp_herings_peeters_example", + ), + pytest.param( + EquilibriumTestCaseWithStart( + factory=games.create_hs1988_base_game, + solver=gbt.nash.hp_solve, + start_data=dict( + data=[[1.0 / 3.0, 2.0 / 3.0], [1.0 / 6.0, 5.0 / 6.0]], rational=False + ), + expected=[[d(0.0, 1.0), d(0.0, 1.0)]], + regret_tol=TOL_LARGE, + prob_tol=TOL_LARGE, + ), + marks=pytest.mark.nash_hp_strategy, + id="test_hp_hs_example_1", + ), +] + + SIMPDIV_CASES = [ pytest.param( EquilibriumTestCaseWithStart( @@ -1718,6 +1748,7 @@ def _one_hot_perturbation(rational: bool) -> gbt.MixedStrategyProfile: CASES = [] CASES += LIAP_STRATEGY_CASES +CASES += HP_STRATEGY_CASES CASES += SIMPDIV_CASES @@ -1748,6 +1779,23 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s assert abs(eq_prob - exp_prob) <= test_case.prob_tol +@pytest.mark.nash +@pytest.mark.nash_hp_strategy +def test_hp_degenerate_t0_prior_raises_error() -> None: + """hp_solve() rejects a prior without a unique best response for some player at t=0, + rather than picking one of the tied best responses arbitrarily. + """ + game = games.create_hs1988_base_game() + prior = game.mixed_strategy_profile( + data=[[2.0 / 3.0, 1.0 / 3.0], [1.0 / 3.0, 2.0 / 3.0]], rational=False + ) + with pytest.raises( + RuntimeError, + match="Multiple best responses found for player 1. Only one best response is allowed.", + ): + gbt.nash.hp_solve(prior) + + ################################################################################################## # NASH SOLVER IN MIXED BEHAVIORS ##################################################################################################