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..f147ad8ec 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 881bd99e1..9dc101932 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 a233448df..e1fdd0b09 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 15e57b6e3..9b71fdd07 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 d51a0dc6c..2ed16a430 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 000000000..1a5898722 --- /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 000000000..fd8d2e9cf --- /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 000000000..4b6af3cca --- /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)