From 348ad289aa7cb1d33d63ed1ea38b82533bf48abf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sat, 4 Jul 2026 22:45:29 +0200 Subject: [PATCH 1/6] Structure ofhp method, ComputeInitialPoint and ExtractEquilibrium methods --- Makefile.am | 4 +- src/pygambit/gambit.pxd | 4 +- src/pygambit/nash.pxi | 5 +- src/pygambit/nash.py | 5 +- src/solvers/hp/hp.cc | 24 ++++++-- src/solvers/hp/hp.h | 3 +- src/solvers/hp/hpsystem.cc | 109 +++++++++++++++++++++++++++++++++++++ src/solvers/hp/hpsystem.h | 60 ++++++++++++++++++++ 8 files changed, 203 insertions(+), 11 deletions(-) create mode 100644 src/solvers/hp/hpsystem.cc create mode 100644 src/solvers/hp/hpsystem.h diff --git a/Makefile.am b/Makefile.am index f2ea3fda3..fde0f3ca1 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 644228734..9701eecdc 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_Game, c_MixedStrategyProfile[double] + ) except +RuntimeError cdef extern from "nash.h": diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index 881bd99e1..c45e45dfb 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(game: Game, + prior: MixedStrategyProfileDouble) -> list[MixedStrategyProfileDouble]: + return _convert_mspd(HPStrategySolve(game.game, deref(prior.profile))) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index a233448df..27fa34355 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -808,6 +808,7 @@ def logit_solve( def hp_solve( game: libgbt.Game, + prior: libgbt.MixedStrategyProfileDouble, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -818,13 +819,15 @@ def hp_solve( ---------- 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(game, prior) return NashComputationResult( game=game, method="hp", diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 15e57b6e3..f50290869 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -23,16 +23,30 @@ #include #include "gambit.h" #include "solvers/hp/hp.h" +#include "solvers/hp/hpsystem.h" namespace Gambit { -std::list> HPStrategySolve(const Game &p_game) +std::list> +HPStrategySolve(const Game &p_game, const MixedStrategyProfile &p_prior) { - std::list> result; + const std::list> result; const StrategySupportProfile support(p_game); - const MixedStrategyProfile trivial_profile = support.NewMixedStrategyProfile(); + const MixedStrategyProfile prior(p_prior); - result.push_back(trivial_profile); - return result; + const HPEquationSystem system(p_game, p_prior); + Vector init_pt = system.ComputeInitialPoint(); + + std::cout << "Computed starting vector (alpha form): "; + for (size_t i = 1; i <= init_pt.size(); ++i) { + std::cout << init_pt[i] << " "; + } + std::cout << std::endl; + + std::list> ret; + const MixedStrategyProfile T_0 = system.ExtractEquilibrium(init_pt); + ret.push_back(T_0); + + return ret; } } // namespace Gambit diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index d51a0dc6c..bda3fa71f 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 Game &p_game, 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 000000000..8eeb7d4b0 --- /dev/null +++ b/src/solvers/hp/hpsystem.cc @@ -0,0 +1,109 @@ +// +// 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 Game &game, const MixedStrategyProfile &prior) + : m_game(game), m_prior(prior) +{ +} + +Vector HPEquationSystem::ComputeInitialPoint() const +{ + const int n_players = m_game->GetPlayers().size(); + int m_star = 0; + const double tol = 1e-9; // Tolerance for floating-point comparisons + for (const auto &player : m_game->GetPlayers()) { + m_star += player->GetStrategies().size(); + } + + // 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; + + for (const auto &player : m_game->GetPlayers()) { + + double max_payoff = -std::numeric_limits::infinity(); + std::vector payoffs_against_prior; + payoffs_against_prior.reserve(player->GetStrategies().size()); + + // Finding mu^i (the maximum payoff for player i against the prior) + for (const auto &strategy : player->GetStrategies()) { + const double payoff = m_prior.GetPayoff(strategy); + payoffs_against_prior.push_back(payoff); + + 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 + int local_s_idx = 0; + bool found_br = false; // Flag to check if a best response has been found + + for (const auto &strategy : player->GetStrategies()) { + const double lambda = max_payoff - 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 { + // 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; +} + +} // end namespace Gambit diff --git a/src/solvers/hp/hpsystem.h b/src/solvers/hp/hpsystem.h new file mode 100644 index 000000000..177ca7847 --- /dev/null +++ b/src/solvers/hp/hpsystem.h @@ -0,0 +1,60 @@ +// +// 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 Game &game, 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; + + // 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; } + + // Dynamic payoffs + double CalculateMarginalPayoff(int player, int strategy_idx, + const MixedStrategyProfile ¤t_sigma, + double t) const; +}; + +} // namespace Gambit +#endif // HPSYSTEM_H From 93f72b1008e67cedd83bec8ab96c31608ea7b2e1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sat, 11 Jul 2026 22:18:27 +0200 Subject: [PATCH 2/6] Exception and elimination of game parameter --- src/pygambit/gambit.pxd | 2 +- src/pygambit/nash.pxi | 4 ++-- src/pygambit/nash.py | 7 ++----- src/solvers/hp/hp.cc | 7 ++----- src/solvers/hp/hp.h | 2 +- src/solvers/hp/hpsystem.cc | 9 +++++++-- src/solvers/hp/hpsystem.h | 4 ++-- 7 files changed, 17 insertions(+), 18 deletions(-) diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index 9701eecdc..f147ad8ec 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -623,7 +623,7 @@ cdef extern from "solvers/logit/logit.h": cdef extern from "solvers/hp/hp.h": stdlist[c_MixedStrategyProfile[double]] HPStrategySolve( - c_Game, c_MixedStrategyProfile[double] + c_MixedStrategyProfile[double] ) except +RuntimeError diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index c45e45dfb..9dc101932 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -370,6 +370,6 @@ def _logit_behavior_branch(game: Game, return ret -def _hp_strategy_solve(game: Game, +def _hp_strategy_solve( prior: MixedStrategyProfileDouble) -> list[MixedStrategyProfileDouble]: - return _convert_mspd(HPStrategySolve(game.game, deref(prior.profile))) + return _convert_mspd(HPStrategySolve(deref(prior.profile))) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 27fa34355..e1fdd0b09 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -807,7 +807,6 @@ def logit_solve( def hp_solve( - game: libgbt.Game, prior: libgbt.MixedStrategyProfileDouble, ) -> NashComputationResult: """Compute Nash equilibria of a game using :cite:p:`HerPee01` @@ -817,8 +816,6 @@ def hp_solve( Parameters ---------- - game : Game - The game to compute equilibria in. prior : MixedStrategyProfileDouble The prior distribution over strategies. @@ -827,9 +824,9 @@ def hp_solve( res : NashComputationResult The result represented as a ``NashComputationResult`` object. """ - equilibria = libgbt._hp_strategy_solve(game, prior) + 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 f50290869..02521f58c 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -27,14 +27,11 @@ namespace Gambit { std::list> -HPStrategySolve(const Game &p_game, const MixedStrategyProfile &p_prior) +HPStrategySolve(const MixedStrategyProfile &p_prior) { const std::list> result; - const StrategySupportProfile support(p_game); - const MixedStrategyProfile prior(p_prior); - - const HPEquationSystem system(p_game, p_prior); + const HPEquationSystem system(p_prior); Vector init_pt = system.ComputeInitialPoint(); std::cout << "Computed starting vector (alpha form): "; diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h index bda3fa71f..2ed16a430 100644 --- a/src/solvers/hp/hp.h +++ b/src/solvers/hp/hp.h @@ -27,7 +27,7 @@ namespace Gambit { std::list> -HPStrategySolve(const Game &p_game, const MixedStrategyProfile &p_prior); +HPStrategySolve(const MixedStrategyProfile &p_prior); } // namespace Gambit #endif // HP_H diff --git a/src/solvers/hp/hpsystem.cc b/src/solvers/hp/hpsystem.cc index 8eeb7d4b0..1d0b269b8 100644 --- a/src/solvers/hp/hpsystem.cc +++ b/src/solvers/hp/hpsystem.cc @@ -25,8 +25,8 @@ #include "solvers/hp/hpsystem.h" namespace Gambit { -HPEquationSystem::HPEquationSystem(const Game &game, const MixedStrategyProfile &prior) - : m_game(game), m_prior(prior) +HPEquationSystem::HPEquationSystem(const MixedStrategyProfile &prior) + : m_game(prior.GetGame()), m_prior(prior) { } @@ -77,6 +77,11 @@ Vector HPEquationSystem::ComputeInitialPoint() const 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)); diff --git a/src/solvers/hp/hpsystem.h b/src/solvers/hp/hpsystem.h index 177ca7847..6bbe71bf0 100644 --- a/src/solvers/hp/hpsystem.h +++ b/src/solvers/hp/hpsystem.h @@ -29,7 +29,7 @@ namespace Gambit { class HPEquationSystem { public: - HPEquationSystem(const Game &game, const MixedStrategyProfile &prior); + HPEquationSystem(const MixedStrategyProfile &prior); // Evaluates H(t, alpha, mu) = 0 void GetValue(const Vector &point, Vector &lhs) const; @@ -43,7 +43,7 @@ class HPEquationSystem { MixedStrategyProfile ExtractEquilibrium(const Vector &final_point) const; private: - const Game &m_game; + const Game m_game; MixedStrategyProfile m_prior; // Transforms alpha to sigma and lambda From 25e707ba02ea650d383ab85960c83179655328c2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sun, 12 Jul 2026 22:15:22 +0200 Subject: [PATCH 3/6] Tests --- tests/test_hp.py | 135 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 135 insertions(+) create mode 100644 tests/test_hp.py diff --git a/tests/test_hp.py b/tests/test_hp.py new file mode 100644 index 000000000..52e5387b3 --- /dev/null +++ b/tests/test_hp.py @@ -0,0 +1,135 @@ +"""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) + + # El prior que descubriste accidentalmente + 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.8, 0.2), d(1.0 / 3.0, 2.0 / 3.0)], + ), + marks=pytest.mark.xfail(reason="Mathematical curve tracking not yet fully implemented"), + 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)], + ), + marks=pytest.mark.xfail(reason="Mathematical curve tracking not yet fully implemented"), + 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 d21e976e93f91055b032078aa5039d41d11c105b 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: Tue, 14 Jul 2026 15:02:54 +0200 Subject: [PATCH 4/6] Tests correction Updated expected results for HPSolver test case. --- tests/test_hp.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/test_hp.py b/tests/test_hp.py index 52e5387b3..4bfa01ec1 100644 --- a/tests/test_hp.py +++ b/tests/test_hp.py @@ -70,7 +70,6 @@ def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: prior = game.mixed_strategy_profile() p1, p2 = list(game.players) - # El prior que descubriste accidentalmente 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 @@ -83,7 +82,7 @@ def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: pytest.param( HPSolverTestCase( factory=create_hp_paper_example, - expected=[d(0.8, 0.2), d(1.0 / 3.0, 2.0 / 3.0)], + expected=[d(0.0, 1.0), d(0.0, 1.0)], ), marks=pytest.mark.xfail(reason="Mathematical curve tracking not yet fully implemented"), id="test_hp_herings_peeters_example", From 8503f5c880674bb4b7e36e726cf6618aa10e67ff Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Tue, 14 Jul 2026 15:40:10 +0200 Subject: [PATCH 5/6] HP Algorithm --- src/solvers/hp/hp.cc | 43 +++++++--- src/solvers/hp/hpsystem.cc | 157 +++++++++++++++++++++++++++++++++---- src/solvers/hp/hpsystem.h | 18 ++++- tests/test_hp.py | 2 - 4 files changed, 186 insertions(+), 34 deletions(-) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 02521f58c..5dfea6993 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -24,26 +24,45 @@ #include "gambit.h" #include "solvers/hp/hp.h" #include "solvers/hp/hpsystem.h" +#include "solvers/logit/path.h" namespace Gambit { std::list> HPStrategySolve(const MixedStrategyProfile &p_prior) { - const std::list> result; - const HPEquationSystem system(p_prior); - Vector init_pt = system.ComputeInitialPoint(); + std::list> equilibria; - std::cout << "Computed starting vector (alpha form): "; - for (size_t i = 1; i <= init_pt.size(); ++i) { - std::cout << init_pt[i] << " "; - } - std::cout << std::endl; + HPEquationSystem system(p_prior); + Vector x = system.ComputeInitialPoint(); - std::list> ret; - const MixedStrategyProfile T_0 = system.ExtractEquilibrium(init_pt); - ret.push_back(T_0); + const PathTracer tracer; + double omega = 1.0; - return ret; + auto termination_condition = [](const Vector &point) { 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; + }); + + equilibria.push_back(system.ExtractEquilibrium(x)); + return equilibria; } } // namespace Gambit diff --git a/src/solvers/hp/hpsystem.cc b/src/solvers/hp/hpsystem.cc index 1d0b269b8..1a5898722 100644 --- a/src/solvers/hp/hpsystem.cc +++ b/src/solvers/hp/hpsystem.cc @@ -26,18 +26,137 @@ namespace Gambit { HPEquationSystem::HPEquationSystem(const MixedStrategyProfile &prior) - : m_game(prior.GetGame()), m_prior(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(); - int m_star = 0; const double tol = 1e-9; // Tolerance for floating-point comparisons - for (const auto &player : m_game->GetPlayers()) { - m_star += player->GetStrategies().size(); - } // Dimension: 1 (t) + m_star (total strategies) + n (number of players) const int vector_size = 1 + m_star + n_players; @@ -47,18 +166,14 @@ Vector HPEquationSystem::ComputeInitialPoint() const 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()) { - - double max_payoff = -std::numeric_limits::infinity(); - std::vector payoffs_against_prior; - payoffs_against_prior.reserve(player->GetStrategies().size()); - + 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_prior.GetPayoff(strategy); - payoffs_against_prior.push_back(payoff); - + const double payoff = m_payoffs_against_prior[flat_strategy_idx++]; if (payoff > max_payoff) { max_payoff = payoff; } @@ -68,11 +183,11 @@ Vector HPEquationSystem::ComputeInitialPoint() const start_point[1 + m_star + player_idx] = max_payoff; // Compute alpha^i_s for each strategy s of player i - int local_s_idx = 0; - bool found_br = false; // Flag to check if a best response has been found + 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 - payoffs_against_prior[local_s_idx++]; + 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; @@ -111,4 +226,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 +{ + 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 index 6bbe71bf0..fd8d2e9cf 100644 --- a/src/solvers/hp/hpsystem.h +++ b/src/solvers/hp/hpsystem.h @@ -45,15 +45,25 @@ class HPEquationSystem { 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; } - // Dynamic payoffs - double CalculateMarginalPayoff(int player, int strategy_idx, - const MixedStrategyProfile ¤t_sigma, - double t) const; + // 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 diff --git a/tests/test_hp.py b/tests/test_hp.py index 4bfa01ec1..4b6af3cca 100644 --- a/tests/test_hp.py +++ b/tests/test_hp.py @@ -84,7 +84,6 @@ def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: factory=create_hp_paper_example, expected=[d(0.0, 1.0), d(0.0, 1.0)], ), - marks=pytest.mark.xfail(reason="Mathematical curve tracking not yet fully implemented"), id="test_hp_herings_peeters_example", ), pytest.param( @@ -92,7 +91,6 @@ def create_t0_degenerate_example() -> gbt.MixedStrategyProfileDouble: factory=create_hs_example_1, expected=[d(0.0, 1.0), d(0.0, 1.0)], ), - marks=pytest.mark.xfail(reason="Mathematical curve tracking not yet fully implemented"), id="test_hp_hs_example_1", ), ] From fb530db2a8461127f925d99dcf9f6eca431d624c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Andr=C3=A9s=20Fern=C3=A1ndez=20Cervell?= Date: Sun, 19 Jul 2026 18:54:42 +0200 Subject: [PATCH 6/6] termination in t=1 --- src/solvers/hp/hp.cc | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc index 5dfea6993..9b71fdd07 100644 --- a/src/solvers/hp/hp.cc +++ b/src/solvers/hp/hp.cc @@ -39,7 +39,10 @@ HPStrategySolve(const MixedStrategyProfile &p_prior) const PathTracer tracer; double omega = 1.0; - auto termination_condition = [](const Vector &point) { return point[1] >= 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) { @@ -60,7 +63,8 @@ HPStrategySolve(const MixedStrategyProfile &p_prior) std::cout << prob_vector[i] << " "; } std::cout << std::endl; - }); + }, + criterion_function); equilibria.push_back(system.ExtractEquilibrium(x)); return equilibria;