diff --git a/Makefile.am b/Makefile.am index 1e984d47e..232ed6e38 100644 --- a/Makefile.am +++ b/Makefile.am @@ -277,6 +277,12 @@ gtracer_SOURCES = \ src/solvers/gtracer/gnm.cc \ src/solvers/gtracer/ipa.cc +hp_SOURCES = \ + src/solvers/hp/hp.h \ + src/solvers/hp/hp.cc \ + src/solvers/hp/hpsystem.h \ + src/solvers/hp/hpsystem.cc + nashsupport_SOURCES = \ src/solvers/nashsupport/efgsupport.cc \ src/solvers/nashsupport/nfgsupport.cc \ @@ -306,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 @@ -357,14 +370,14 @@ AM_CXXFLAGS = ${LLVM_CXXFLAGS} -Wall -Wsign-compare -Wunreachable-code ## 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 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} -liblogit_a_SOURCES = ${logit_SOURCES} +libhomotopy_a_SOURCES = ${homotopy_SOURCES} libsimpdiv_a_SOURCES = ${simpdiv_SOURCES} libenumpoly_a_SOURCES = ${enumpoly_SOURCES} @@ -443,7 +456,7 @@ 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 \ +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/doc/algorithms.rst b/doc/algorithms.rst index 4238a1555..db71bc87d 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -15,12 +15,49 @@ Algorithm Description :ref:`lp` Compute equilibria in a two-player constant-sum game via linear programming :py:func:`pygambit.nash.lp_solve` :ref:`gambit-lp ` :ref:`lcp` Compute equilibria in a two-player game via linear complementarity :py:func:`pygambit.nash.lcp_solve` :ref:`gambit-lcp ` :ref:`liap` Compute equilibria using function minimization :py:func:`pygambit.nash.liap_solve` :ref:`gambit-liap ` -:ref:`logit` Trace logit QRE and approximate a Nash equilibrium at high precision :py:func:`pygambit.nash.logit_solve` :ref:`gambit-logit ` +:ref:`logit` Trace logit QRE and approximate a Nash equilibrium at high precision :py:func:`pygambit.nash.logit_solve` :ref:`gambit-logit ` :ref:`simpdiv` Compute equilibria via simplicial subdivision :py:func:`pygambit.nash.simpdiv_solve` :ref:`gambit-simpdiv ` :ref:`ipa` Compute equilibria using iterated polymatrix approximation :py:func:`pygambit.nash.ipa_solve` :ref:`gambit-ipa ` :ref:`gnm` Compute equilibria using a global Newton method :py:func:`pygambit.nash.gnm_solve` :ref:`gambit-gnm ` +:ref:`hp` Compute a specific Nash equilibrium using a homotopy path-following method :py:func:`pygambit.nash.hp_solve` :ref:`gambit-hp ` ================ =========================================================================== ======================================== ========================================== +.. _pygambit-nash-maxregret: + +Acceptance criteria for approximate Nash equilibria +---------------------------------------------------- + +Some methods for computing Nash equilibria operate using floating-point arithmetic and/or +generate candidate equilibrium profiles using methods which involve some form of successive +approximation. The outputs of these methods therefore are in general +:math:`\varepsilon`-equilibria, for some positive :math:`\varepsilon`: a strategy profile at +which no player can gain more than :math:`\varepsilon` in expected payoff by unilaterally +deviating. Every Nash equilibrium is an :math:`\varepsilon`-equilibrium with +:math:`\varepsilon = 0`. + +To provide a uniform interface across methods, where relevant Gambit provides a parameter +`maxregret`, which specifies the acceptance criterion for labeling the output of the +algorithm as an equilibrium. This parameter is interpreted *proportionally* to the range of +payoffs in the game: any profile returned as an equilibrium is guaranteed to be an +:math:`\varepsilon`-equilibrium, for :math:`\varepsilon` no more than `maxregret` times the +difference of the game's maximum and minimum payoffs. For example, with the default +`maxregret` of :math:`10^{-8}` in a game whose payoffs range over 4 units, any equilibrium +returned is guaranteed to have a regret, measured directly in the game's own payoffs, of no +more than :math:`4 \times 10^{-8}`. + +Expressing `maxregret` scaled by the game's payoffs in this way standardises the behavior of +methods across games: for instance, doubling all the payoffs in a game does not change the +`maxregret` value needed to obtain equilibria of comparable quality. + +Methods differ in the guarantees they offer once a `maxregret` criterion is specified. +Globally-convergent methods, such as :ref:`logit` and :ref:`gnm`, are guaranteed to eventually +satisfy any `maxregret` criterion, though a tighter criterion generally requires more +computation. Other methods, such as :ref:`liap`, are not globally convergent, and may fail to +find any equilibrium satisfying the criterion from a given starting point. + +See the :doc:`stripped-down poker tutorial ` for a worked +example comparing `maxregret` across several methods. + .. _enumpure: enumpure @@ -234,3 +271,20 @@ The algorithm takes as a parameter a mixed strategy profile. This profile is interpreted as defining a ray in the space of games. The profile must have the property that, for each player, the most frequently played strategy must be unique. + +.. _hp: + +hp +--- +Computes the Nash equilibrium selected by the tracing procedure +of Harsanyi and Selten using a homotopy path-following method. The algorithm +was first described by Herings and Peeters :cite:p:`HerPee01`. + +The algorithm takes as a parameter a mixed strategy profile, which acts as +the subjective prior beliefs of the players. +The profile must have the property that, for each player, +there must only exist one best response. + +For generic games, the algorithm converges to the unique Nash equilibrium selected by +the tracing procedure of Harsanyi and Selten. For non-generic games, the algorithm may +converge to a Nash equilibrium that is not selected by the tracing procedure. diff --git a/doc/gui.nash.rst b/doc/gui.nash.rst index 2df8238e4..5a07164af 100644 --- a/doc/gui.nash.rst +++ b/doc/gui.nash.rst @@ -87,6 +87,7 @@ Method Parameters used by the graphical interface ``ipa`` One random perturbation. ``gnm`` One random perturbation; ending lambda ``-10``; 100 steps per support cell; local Newton refinement every 3 steps, with at most 10 iterations. +``hp`` No method-specific parameters. ================ ============================================================================ For extensive games, there is an option of whether to use the diff --git a/doc/pygambit.api.rst b/doc/pygambit.api.rst index 6a66ed49e..cd7da0de3 100644 --- a/doc/pygambit.api.rst +++ b/doc/pygambit.api.rst @@ -355,6 +355,7 @@ Computation of Nash equilibria simpdiv_solve ipa_solve gnm_solve + hp_solve Computation of quantal response equilibria diff --git a/doc/references.bib b/doc/references.bib index 1998c5ef4..a655b2f14 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/doc/tools.hp.rst b/doc/tools.hp.rst new file mode 100644 index 000000000..86770fa32 --- /dev/null +++ b/doc/tools.hp.rst @@ -0,0 +1,81 @@ +.. _gambit-hp: + +:program:`gambit-hp` +===================== + +Compute a Nash equilibrium in a strategic game using the homotopy method of +:cite:t:`HerPee01`. + +The algorithm finds one equilibrium starting from any given prior +distribution over strategies, which must have a unique best response for +each player. Multiple prior distributions may be generated via the `-n` +option or specified via the `-s` option; different priors may result in +different equilibria being found. + + +.. program:: gambit-hp + +.. cmdoption:: -d + + Express all output using decimal representations + with the specified number of digits. + +.. cmdoption:: -h + + Prints a help message listing the available options. + +.. cmdoption:: -m + + Specify the maximum regret criterion for acceptance as an approximate Nash equilibrium + (default is 1e-8). See :ref:`pygambit-nash-maxregret` for interpretation and guidance. + +.. cmdoption:: -n + + Randomly generate the specified number of prior distributions. + Mutually exclusive with :option:`-s`. + +.. cmdoption:: -R + + Seeds the random number generator used to generate prior + distributions with the specified value, so that the sequence of priors + generated by :option:`-n` can be reproduced across runs. If not + specified, the generator is seeded from system entropy. Requires + :option:`-n`. + +.. cmdoption:: -q + + Suppresses printing of the banner at program launch. + +.. cmdoption:: -s + + Specifies a file containing a list of prior distributions over + strategies. The format of the file is comma-separated values, + one mixed strategy profile per line, in the same format used for + output of equilibria (excluding the initial NE tag). + Mutually exclusive with :option:`-n`. + +.. cmdoption:: -V, --verbose + + Show the prior distribution itself, tagged `prior`, followed by each + point traced along the homotopy path, tagged with the homotopy + parameter t in place of the NE tag. Note that the point at t=0 is the + best response to the prior, and so is generally a pure strategy profile + even when the prior itself is not. If this option is not specified, + only the equilibrium found is reported. + +.. cmdoption:: -v, --version + + Prints version information and exits. + + +Computing an equilibrium of the reduced strategic form of the example in +Figure 2 of :cite:p:`Sel75`, starting from the prior in which player 1 +plays (0.5, 0.3, 0.2) and player 2 plays (0.6, 0.4):: + + $ echo "0.5,0.3,0.2,0.6,0.4" > prior.csv + $ gambit-hp -s prior.csv catalog/journals/ijgt/selten1975/fig2.efg + Compute a Nash equilibrium using the Herings-Peeters (2001) homotopy method + Gambit version |release|, Copyright (C) 1994-2026, The Gambit Project + This is free software, distributed under the GNU GPL + + NE,1.000000,0.000000,0.000000,0.999976,0.000024 diff --git a/doc/tools.rst b/doc/tools.rst index 2efab1a95..99445adf1 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/doc/tutorials/03_stripped_down_poker.ipynb b/doc/tutorials/03_stripped_down_poker.ipynb index fab307ac2..cc890eec5 100644 --- a/doc/tutorials/03_stripped_down_poker.ipynb +++ b/doc/tutorials/03_stripped_down_poker.ipynb @@ -675,29 +675,7 @@ "cell_type": "markdown", "id": "b2867dca", "metadata": {}, - "source": [ - "Acceptance criteria for Nash equilibria\n", - "---------------------------------------\n", - "\n", - "Some methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations.\n", - "The outputs of these methods therefore are in general $\\varepsilon$-equilibria, for some positive $\\varepsilon$.\n", - "\n", - "$\\varepsilon$-equilibria (from [Wikipedia](https://en.wikipedia.org/wiki/Epsilon-equilibrium)):\n", - "\n", - "> In game theory, an epsilon-equilibrium, or near-Nash equilibrium, is a strategy profile that approximately satisfies the condition of Nash equilibrium. In a Nash equilibrium, no player has an incentive to change his behavior. In an approximate Nash equilibrium, this requirement is weakened to allow the possibility that a player may have a small incentive to do something different.\n", - "\n", - "> Given a game and a real non-negative parameter $\\varepsilon$, a strategy profile is said to be an $\\varepsilon$-equilibrium if it is not possible for any player to gain more than $\\varepsilon$ in expected payoff by unilaterally deviating from his strategy. Every Nash Equilibrium is an $\\varepsilon$-equilibrium where $\\varepsilon = 0$.\n", - "\n", - "\n", - "To provide a uniform interface across methods, where relevant Gambit provides a parameter\n", - "`maxregret`, which specifies the acceptance criterion for labeling the output of the\n", - "algorithm as an equilibrium.\n", - "This parameter is interpreted *proportionally* to the range of payoffs in the game.\n", - "Any profile returned as an equilibrium is guaranteed to be an $\\varepsilon$-equilibrium, for $\\varepsilon$ no more than `maxregret`\n", - "times the difference of the game's maximum and minimum payoffs.\n", - "\n", - "As an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" - ] + "source": "Acceptance criteria for Nash equilibria\n---------------------------------------\n\nSome methods for computing Nash equilibria operate using floating-point arithmetic and/or generate candidate equilibrium profiles using methods which involve some form of successive approximations, and so return only approximate equilibria.\nGambit expresses how close an approximation must be, before it is reported as an equilibrium, via a `maxregret` parameter, interpreted *proportionally* to the range of payoffs in the game.\nSee [Acceptance criteria for approximate Nash equilibria](../algorithms.html#pygambit-nash-maxregret) in the algorithms documentation for the full explanation of `maxregret` and the guarantees it provides.\n\nAs an example, consider solving our one-card poker game using `logit_solve`. The range of the payoffs in this game is 4 (from +2 to -2):\n" }, { "cell_type": "code", diff --git a/pyproject.toml b/pyproject.toml index 9dd604267..b3c9e30c4 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"] @@ -124,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/setup.py b/setup.py index 4ebb4d8a3..35a70b841 100644 --- a/setup.py +++ b/setup.py @@ -95,7 +95,7 @@ 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"]) @@ -111,7 +111,7 @@ def run(self) -> None: setuptools.setup( cmdclass={"build_py": GambitBuildPy}, - libraries=[cppgambit_bimatrix, cppgambit_liap, cppgambit_logit, cppgambit_simpdiv, + libraries=[cppgambit_bimatrix, cppgambit_liap, cppgambit_homotopy, cppgambit_simpdiv, cppgambit_gtracer, cppgambit_enumpoly, cppgambit_games, cppgambit_core], ext_modules=Cython.Build.cythonize(libgambit, diff --git a/src/gui/dllogit.h b/src/gui/dllogit.h index d6462b2cc..c946e6419 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 47c2a47a0..5572e7b1f 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,10 @@ 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 -m %.17g", method.priors, method.maxRegret); + } else if constexpr (std::is_same_v) { return prefix + wxString::Format("ipa -d 10 -n %d", method.perturbations); } @@ -223,6 +231,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 +274,10 @@ 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; maximum regret %.4g)", method.priors, + method.maxRegret); + } else if constexpr (std::is_same_v) { return wxString::Format(" (%d perturbation)", method.perturbations); } @@ -387,6 +402,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 4c8ee29a2..54be81996 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, spec.maxRegret, + [&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 50af9a87d..d138cf9f8 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/nash.h" +#include "solvers/path/path.h" namespace Gambit::GUI { @@ -94,6 +97,16 @@ 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}; + double maxRegret{1.0e-8}; + + std::optional MakeSolver(NashRepresentation) const; +}; + struct LPNashSpec { std::optional MakeSolver(NashRepresentation p_representation) const; }; @@ -115,7 +128,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 +145,9 @@ struct SimpdivNashSpec { std::optional MakeSolver(NashRepresentation) const; }; -using NashMethodSpec = - std::variant; +using NashMethodSpec = std::variant; struct NashComputationSpec { NashRepresentation representation; diff --git a/src/pygambit/callback.h b/src/pygambit/callback.h index 7982c86b3..bd99dbbaa 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/cli/hp.py b/src/pygambit/cli/hp.py new file mode 100644 index 000000000..41acc7890 --- /dev/null +++ b/src/pygambit/cli/hp.py @@ -0,0 +1,127 @@ +# +# 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( + "-m", + "maxregret", + default=1.0e-8, + show_default=True, + type=float, + help="maximum regret acceptable as a proportion of the range of payoffs in the game", +) +@click.option("-q", "--quiet", is_flag=True, help="quiet mode (suppresses banner)") +@click.option( + "-V", + "--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, + maxregret: float, + 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, maxregret=maxregret, event_callback=render_event) + for eq in result.equilibria: + click.echo(render_profile_csv(eq, "NE", decimals)) + + +if __name__ == "__main__": + main() diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index e0140f340..1c8c16f73 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) @@ -697,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], double, HPEventCallbackType + ) except +RuntimeError diff --git a/src/pygambit/nash.h b/src/pygambit/nash.h index 96091b85b..815a27d78 100644 --- a/src/pygambit/nash.h +++ b/src/pygambit/nash.h @@ -21,7 +21,9 @@ // #include "solvers/enummixed/enummixed.h" +#include "solvers/hp/hp.h" #include "solvers/logit/logit.h" +#include "solvers/path/path.h" using namespace std; using namespace Gambit; @@ -47,9 +49,9 @@ LogitBehaviorSolveWrapper(const Game &p_game, double p_regret, double p_firstSte NullLogitEventCallback) { std::list> ret; - ret.push_back(LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, 1.0, - p_firstStep, p_maxAccel, Nash::NullBehaviorCallback, - p_onEvent) + ret.push_back(LogitBehaviorSolve(LogitQREMixedBehaviorProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel, + Nash::NullBehaviorCallback, p_onEvent) .back() .GetProfile()); return ret; @@ -59,8 +61,8 @@ 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 @@ -69,8 +71,9 @@ LogitBehaviorEstimateWrapper(std::shared_ptr> p_fre LogitEventCallbackType p_onEvent = NullLogitEventCallback) { - return make_shared(LogitBehaviorEstimate( - *p_frequencies, 1000000.0, 1.0, p_stopAtLocal, p_firstStep, p_maxAccel, p_onEvent)); + return make_shared( + LogitBehaviorEstimate(*p_frequencies, 1000000.0, PathTracer::TraceDirection::Positive, + p_stopAtLocal, p_firstStep, p_maxAccel, p_onEvent)); } std::list> @@ -82,7 +85,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, p_onEvent)) { + LogitBehaviorSolveLambda(start, p_targetLambda, PathTracer::TraceDirection::Positive, + p_firstStep, p_maxAccel, p_onEvent)) { ret.push_back(std::make_shared(qre)); } return ret; @@ -95,9 +99,9 @@ LogitStrategySolveWrapper(const Game &p_game, double p_regret, double p_firstSte NullLogitEventCallback) { std::list> ret; - ret.push_back(LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, 1.0, - p_firstStep, p_maxAccel, Nash::NullStrategyCallback, - p_onEvent) + ret.push_back(LogitStrategySolve(LogitQREMixedStrategyProfile(p_game), p_regret, + PathTracer::TraceDirection::Positive, p_firstStep, p_maxAccel, + Nash::NullStrategyCallback, p_onEvent) .back() .GetProfile()); return ret; @@ -107,8 +111,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> @@ -120,7 +124,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, p_onEvent)) { + LogitStrategySolveLambda(start, p_targetLambda, PathTracer::TraceDirection::Positive, + p_firstStep, p_maxAccel, p_onEvent)) { ret.push_back(std::make_shared(qre)); } return ret; @@ -132,6 +137,15 @@ LogitStrategyEstimateWrapper(std::shared_ptr> p_fre LogitEventCallbackType p_onEvent = NullLogitEventCallback) { - return make_shared(LogitStrategyEstimate( - *p_frequencies, 1000000.0, 1.0, p_stopAtLocal, p_firstStep, p_maxAccel, p_onEvent)); + return make_shared( + LogitStrategyEstimate(*p_frequencies, 1000000.0, PathTracer::TraceDirection::Positive, + p_stopAtLocal, p_firstStep, p_maxAccel, p_onEvent)); +} + +std::list> +HPStrategySolveWrapper(const MixedStrategyProfile &p_prior, double p_maxRegret, + Nash::HPEventCallbackType p_onEvent = Nash::NullHPEventCallback) +{ + return Nash::HPStrategySolve(p_prior, p_maxRegret, Nash::NullStrategyCallback, + p_onEvent); } diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index b4e8d6677..595a42223 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]] ): @@ -893,3 +910,13 @@ def _logit_behavior_branch(game: Game, max_accel: float): solns = LogitBehaviorPrincipalBranchWrapper(game.game, maxregret, first_step, max_accel) return [LogitQREMixedBehaviorProfile.wrap(profile) for profile in make_list_of_pointer(solns)] + + +def _hp_strategy_solve( + prior: MixedStrategyProfileDouble, + maxregret: float, + event_callback: object = None, +) -> list[MixedStrategyProfileDouble]: + return _convert_mspd(HPStrategySolveWrapper( + deref(prior.profile), maxregret, MakeHPEventCallback(event_callback) + )) diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 5c330d6cf..878907902 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -1109,3 +1109,48 @@ def logit_solve( equilibria=equilibria, parameters={"first_step": first_step, "max_accel": max_accel}, ) + + +def hp_solve( + prior: libgbt.MixedStrategyProfileDouble, + maxregret: float = 1.0e-8, + event_callback: Callable[[libgbt.HPStepEvent], None] | None = None, +) -> NashComputationResult: + """Compute Nash equilibria of a game using :cite:p:`HerPee01` + + Returns an approximation to the limiting point on the principal branch of + the homotopy path for the game. + + Parameters + ---------- + prior : MixedStrategyProfileDouble + The prior distribution over strategies. + + maxregret : float, default 1e-8 + The acceptance criterion for approximate Nash equilibrium; the maximum + regret of any player must be no more than `maxregret` times the + difference of the maximum and minimum payoffs of the game + + event_callback : Callable[[HPStepEvent], None], optional + If specified, called with each point traced along the homotopy path, + and the homotopy parameter ``t`` at which it was reached, on the way + to the returned equilibrium. + + .. versionadded:: 17.0.0 + + Returns + ------- + res : NashComputationResult + The result represented as a ``NashComputationResult`` object. + """ + if maxregret <= 0.0: + raise ValueError("hp_solve(): maxregret argument must be positive") + equilibria = libgbt._hp_strategy_solve(prior, maxregret, event_callback) + return NashComputationResult( + game=prior.game, + method="hp", + rational=False, + use_strategic=True, + equilibria=equilibria, + parameters={"maxregret": maxregret}, + ) diff --git a/src/solvers/hp/hp.cc b/src/solvers/hp/hp.cc new file mode 100644 index 000000000..95e640b87 --- /dev/null +++ b/src/solvers/hp/hp.cc @@ -0,0 +1,118 @@ +// +// 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 "gambit.h" +#include "solvers/hp/hp.h" +#include "solvers/hp/hpsystem.h" +#include "solvers/path/path.h" + +namespace Gambit::Nash { +std::list> +HPStrategySolve(const MixedStrategyProfile &p_prior, double p_maxRegret, + StrategyCallbackType p_onEquilibrium, HPEventCallbackType p_onEvent, + const CancelToken &p_cancel) +{ + std::list> equilibria; + + HPEquationSystem system(p_prior); + Vector x = system.ComputeInitialPoint(); + + const PathTracer tracer; + const PathTracer::TraceDirection direction = PathTracer::TraceDirection::Positive; + const size_t tracking_index = 1; // Track the first variable (t) for orientation + const double t_target = 1.0; + double last_t = 0.0; + bool has_crossed = false; + + auto termination_condition = [t_target, &last_t, &has_crossed, p_maxRegret, + &system](const Vector &point) { + const double t = point[1]; + + // Path tracer reaches maximum acceptable regret + if (system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret && + t >= t_target - p_maxRegret) { + return true; + } + + if (t > t_target) { + if (!has_crossed) { + has_crossed = true; + } + else if (t > + last_t + p_maxRegret) { // Criterion function is not working; polish will do the job + return true; + } + } + else { + // Criterion function might take t back to being less than t_target + has_crossed = false; + } + + last_t = t; + return false; + }; + + auto criterion_function = [t_target](const Vector &point, + const Vector &tangent) -> double { + return point[1] - t_target; + }; + + auto polishing_termination_condition = [p_maxRegret, + &system](const Vector &point) -> bool { + return system.ExtractEquilibrium(point).GetMaxRegret() <= p_maxRegret; + }; + + const auto tracing_result = tracer.TracePath( + [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, + [&system](const Vector &point, Matrix &jac) { + system.GetJacobian(point, jac); + }, + x, direction, tracking_index, termination_condition, + [&system, &p_onEvent](const Vector &point) { + const MixedStrategyProfile profile = system.ExtractEquilibrium(point); + p_onEvent(HPStepEvent{.profile = profile, .t = point[1]}); + }, + criterion_function, NullCriterionBracketFunction, p_cancel); + + const PolishResult polishing_result = PolishPoint( + [&system](const Vector &point, Vector &lhs) { system.GetValue(point, lhs); }, + [&system](const Vector &point, Matrix &jac) { + system.GetJacobian(point, jac); + }, + x, t_target, 1, polishing_termination_condition, 100, + [&system, &p_onEvent](const Vector &point) { + const MixedStrategyProfile profile = system.ExtractEquilibrium(point); + p_onEvent(HPStepEvent{.profile = profile, .t = point[1]}); + }); + + if (!polishing_result.status) { + return {}; + } + const MixedStrategyProfile equilibrium = system.ExtractEquilibrium(x); + if (equilibrium.GetMaxRegret() > p_maxRegret) { + return {}; + } + p_onEquilibrium(equilibrium); + equilibria.push_back(equilibrium); + return equilibria; +} +} // namespace Gambit::Nash diff --git a/src/solvers/hp/hp.h b/src/solvers/hp/hp.h new file mode 100644 index 000000000..b600a1e02 --- /dev/null +++ b/src/solvers/hp/hp.h @@ -0,0 +1,55 @@ +// +// 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 GAMBIT_SOLVERS_HP_HP_H +#define GAMBIT_SOLVERS_HP_HP_H + +#include +#include +#include + +#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, double p_maxRegret = 1.0e-8, + StrategyCallbackType p_onEquilibrium = NullStrategyCallback, + HPEventCallbackType p_onEvent = NullHPEventCallback, + const CancelToken &p_cancel = CancelToken()); + +} // namespace Gambit::Nash + +#endif // GAMBIT_SOLVERS_HP_HP_H diff --git a/src/solvers/hp/hpsystem.cc b/src/solvers/hp/hpsystem.cc new file mode 100644 index 000000000..3f590b0b8 --- /dev/null +++ b/src/solvers/hp/hpsystem.cc @@ -0,0 +1,310 @@ +// +// 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 { + +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++; + } + } + } +}; + +// Eq (b): Probability Sum Equation +class ProbabilitySumEquation final : public HPEquation { + int m_first_alpha_idx; + int m_last_alpha_idx; + +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_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 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; + + // Instantiate Best Response Equations + for (const auto &strategy : player->GetStrategies()) { + m_equations.push_back( + std::make_shared(strategy, alpha_idx, mu_idx, flat_strategy_idx)); + alpha_idx++; + flat_strategy_idx++; + } + + // Instantiate Probability Sum Equations + m_equations.push_back(std::make_shared(first_alpha, alpha_idx)); + + player_idx++; + } +} + +void HPEquationSystem::GetValue(const Vector &point, Vector &lhs) const +{ + // Update internal mutable state + UpdateSigma(point); + + // 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); + } +} + +void HPEquationSystem::GetJacobian(const Vector &point, Matrix &p_jac) const +{ + // Update internal mutable state + UpdateSigma(point); + + p_jac = 0.0; + Vector column(point.size()); // Temp vector matching Jacobian column size + + // 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); + } +} + +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 = AlphaToSigma(alpha_val); + ret[strategy] = prob; + } + } + ret = ret.Normalize(); + + return ret; +} + +void HPEquationSystem::UpdateSigma(const Vector &point) const +{ + 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 new file mode 100644 index 000000000..3f3db28ba --- /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 HPEquation; + +class HPEquationSystem { +public: + explicit HPEquationSystem(const MixedStrategyProfile &prior); + ~HPEquationSystem() = default; + + // 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; + std::vector> m_equations; + + void UpdateSigma(const Vector &point) const; +}; + +} // namespace Gambit +#endif // HPSYSTEM_H diff --git a/src/solvers/logit/efglogit.cc b/src/solvers/logit/efglogit.cc index ba9e1b5f1..093460e0f 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 { @@ -309,10 +309,12 @@ 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, Nash::BehaviorCallbackType p_onEquilibrium, - LogitEventCallbackType p_onEvent, const CancelToken &p_cancel) +std::list +LogitBehaviorSolve(const LogitQREMixedBehaviorProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + Nash::BehaviorCallbackType p_onEquilibrium, + LogitEventCallbackType p_onEvent, + const CancelToken &p_cancel) { if (p_start.size() == 0) { return {p_start}; @@ -336,7 +338,7 @@ std::list LogitBehaviorSolve( [&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); }, @@ -347,11 +349,10 @@ std::list LogitBehaviorSolve( return profiles; } -std::list -LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - LogitEventCallbackType p_onEvent) +std::list LogitBehaviorSolveLambda( + const LogitQREMixedBehaviorProfile &p_start, const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + LogitEventCallbackType p_onEvent) { if (p_start.size() == 0) { return {p_start}; @@ -373,7 +374,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; @@ -385,7 +386,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, LogitEventCallbackType p_onEvent) { const LogitQREMixedBehaviorProfile start(p_frequencies.GetGame()); @@ -409,7 +411,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 81e9f6537..a8f0faded 100644 --- a/src/solvers/logit/logit.h +++ b/src/solvers/logit/logit.h @@ -27,6 +27,7 @@ #include #include "solvers/nash.h" +#include "solvers/path/path.h" namespace Gambit { @@ -92,46 +93,46 @@ template using LogitEventCallbackType = std::function void NullLogitEventCallback(const LogitEvent &) {} 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, Nash::StrategyCallbackType p_onEquilibrium = Nash::NullStrategyCallback, LogitEventCallbackType p_onEvent = NullLogitEventCallback, const CancelToken &p_cancel = CancelToken()); -std::list -LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - LogitEventCallbackType p_onEvent = - NullLogitEventCallback); +std::list LogitStrategySolveLambda( + const LogitQREMixedStrategyProfile &p_start, const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + LogitEventCallbackType p_onEvent = + NullLogitEventCallback); 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, LogitEventCallbackType p_onEvent = NullLogitEventCallback); using LogitQREMixedBehaviorProfile = LogitQRE>; std::list LogitBehaviorSolve( - const LogitQREMixedBehaviorProfile &p_start, double p_regret, double p_omega, - double p_firstStep, double p_maxAccel, + const LogitQREMixedBehaviorProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, Nash::BehaviorCallbackType p_onEquilibrium = Nash::NullBehaviorCallback, LogitEventCallbackType p_onEvent = NullLogitEventCallback, const CancelToken &p_cancel = CancelToken()); -std::list -LogitBehaviorSolveLambda(const LogitQREMixedBehaviorProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - LogitEventCallbackType p_onEvent = - NullLogitEventCallback); +std::list LogitBehaviorSolveLambda( + const LogitQREMixedBehaviorProfile &p_start, const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + LogitEventCallbackType p_onEvent = + NullLogitEventCallback); 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, LogitEventCallbackType p_onEvent = NullLogitEventCallback); diff --git a/src/solvers/logit/nfglogit.cc b/src/solvers/logit/nfglogit.cc index d838a2b44..56fc3cdf0 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 { @@ -347,10 +347,12 @@ 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, Nash::StrategyCallbackType p_onEquilibrium, - LogitEventCallbackType p_onEvent, const CancelToken &p_cancel) +std::list +LogitStrategySolve(const LogitQREMixedStrategyProfile &p_start, double p_regret, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + Nash::StrategyCallbackType p_onEquilibrium, + LogitEventCallbackType p_onEvent, + const CancelToken &p_cancel) { if (p_start.size() == 0) { return {p_start}; @@ -375,7 +377,7 @@ std::list LogitStrategySolve( [&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); }, @@ -386,11 +388,10 @@ std::list LogitStrategySolve( return profiles; } -std::list -LogitStrategySolveLambda(const LogitQREMixedStrategyProfile &p_start, - const std::list &p_targetLambda, double p_omega, - double p_firstStep, double p_maxAccel, - LogitEventCallbackType p_onEvent) +std::list LogitStrategySolveLambda( + const LogitQREMixedStrategyProfile &p_start, const std::list &p_targetLambda, + PathTracer::TraceDirection p_direction, double p_firstStep, double p_maxAccel, + LogitEventCallbackType p_onEvent) { if (p_start.size() == 0) { return {p_start}; @@ -412,7 +413,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; @@ -424,7 +425,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, LogitEventCallbackType p_onEvent) { const LogitQREMixedStrategyProfile start(p_frequencies.GetGame()); @@ -448,7 +450,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/path/path.cc similarity index 55% rename from src/solvers/logit/path.cc rename to src/solvers/path/path.cc index 073213044..7db29c9bb 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 @@ -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; @@ -125,27 +125,31 @@ void NewtonStep(Matrix &q, Matrix &b, Vector &u, Vector< // bifurcation point that the tracing gets stuck there as it is not possible // to find a small enough step size to avoid stepping over the bifurcation // point. -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, - CriterionBracketFunctionType p_criterionBracket, - const CancelToken &p_cancel) const +TracePathResult PathTracer::TracePath( + std::function &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, Vector &x, + TraceDirection p_direction, size_t p_trackingIndex, TerminationFunctionType p_terminate, + CallbackFunctionType p_callback, CriterionFunctionType p_criterion, + CriterionBracketFunctionType p_criterionBracket, const CancelToken &p_cancel) const { - const double c_tol = 1.0e-4; // tolerance for corrector iteration - const double c_maxDist = 0.4; // maximal distance to curve - const double c_maxContr = 0.6; // maximal contraction rate in corrector - const double c_eta = 0.1; // perturbation to avoid cancellation - // in calculating contraction rate - double h = m_hStart; // initial stepsize - const double c_hmin = 1.0e-8; // minimal stepsize - const int c_maxIter = 100; // maximum iterations in corrector + const double c_tol = 1.0e-4; // tolerance for corrector iteration + const double c_maxDist = 0.4; // maximal distance to curve + const double c_maxContr = 0.6; // maximal contraction rate in corrector + const double c_eta = 0.1; // perturbation to avoid cancellation + // in calculating contraction rate + double h = m_hStart; // initial stepsize + const double c_hmin = 1.0e-8; // minimal stepsize + const int c_maxIter = 100; // maximum iterations in corrector + const double c_newtonTol = 1.0e-8; // tolerance for Newton convergence bool newton = false; // using Newton steplength (for zero-finding) const double c_pert = 0.0000001; // The size of perturbation to apply to avoid bifurcation traps double pert = 0.0; // The current version of the perturbation being applied double pert_countdown = 0.0; // How much longer (in arclength) to apply perturbation + const double c_orientTol = 1.0e-8; // tolerance for detecting change in orientation + + const double b_tol = 1.0e-10; // Tolerance for perturbing the b matrix in case of singularity + const double b_pert = 1.0e-8; // Perturbation of the b matrix in case of singularity Vector u(x.size()); // t is current tangent at x; newT is tangent at u, which is the next point. @@ -158,25 +162,67 @@ PathTracer::TracePath(std::function &, Vector QRDecomp(b, q); q.GetRow(q.NumRows(), t); p_callback(x); + int steps = 0; + + auto stepsizeBelowMinimum = [&]() -> TracePathResult { + if (newton && std::abs(p_criterion(x, t)) < c_newtonTol) { + return {x, true, + "Path following terminated successfully at point satisfying criterion function.", + steps}; + } + return {x, false, "Stepsize fell below minimum threshold.", steps}; + }; + + bool first_step = true; + double omega = (p_direction == TraceDirection::Positive) ? 1.0 : -1.0; + + if (p_trackingIndex > x.size() || p_trackingIndex < 1) { + return {x, false, "Tracking index exceeds dimension of point vector.", steps}; + } while (!p_terminate(x)) { p_cancel.Check(); bool accept = true; - if (fabs(h) <= c_hmin) { - return {x, false, "Stepsize fell below minimum threshold."}; + if (std::abs(h) <= c_hmin) { + return stepsizeBelowMinimum(); + } + + if (first_step) { + if (std::abs(t[p_trackingIndex]) <= c_orientTol) { + return {x, false, "Initial tangent vector is orthogonal to path-following direction.", + steps}; + } + // Ensure that the tangent is oriented in the same direction as + // the path-following direction. + 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 p_jacobian(u, b); QRDecomp(b, q); + // Perturb the b matrix if it is singular or nearly singular + for (size_t i = 1; i < b.NumRows(); i++) { + if (std::abs(b(i, i)) < b_tol) { + if (b(i, i) < 0) { + b(i, i) -= b_pert; + } + else { + b(i, i) += b_pert; + } + } + } + int iter = 1; double disto = 0.0; while (true) { @@ -208,7 +254,7 @@ PathTracer::TracePath(std::function &, Vector disto = dist; iter++; if (iter > c_maxIter) { - return {x, false, "Maximum iterations exceeded."}; + return {x, false, "Maximum iterations exceeded.", steps}; } } @@ -222,15 +268,15 @@ PathTracer::TracePath(std::function &, Vector // is oriented in the same direction as we were originally following if (pert_countdown == 0.0) { pert = c_pert; - pert_countdown = abs(2 * h); + pert_countdown = std::abs(2 * h); } accept = false; } if (!accept) { h /= m_maxDecel; // PC not accepted; change stepsize and retry - if (fabs(h) <= c_hmin) { - return {x, false, "Stepsize fell below minimum threshold."}; + if (std::abs(h) <= c_hmin) { + return stepsizeBelowMinimum(); } continue; } @@ -248,31 +294,102 @@ PathTracer::TracePath(std::function &, Vector p_criterionBracket(x, u); } - if (newton) { + const double diff = p_criterion(u, newT) - p_criterion(x, t); + if (newton && std::abs(diff) > c_newtonTol) { // Newton-type steplength adaptation, secant method - h *= -p_criterion(u, newT) / (p_criterion(u, newT) - p_criterion(x, t)); + h *= -p_criterion(u, newT) / diff; } else { // Standard steplength adaptation - h = fabs(h / decel); + h = std::abs(h / decel); } // PC step was successful; update and iterate x = u; t = newT; p_callback(x); + steps++; if (pert_countdown > 0.0) { // If we are currently perturbing in the neighborhood of a bifurcation, check to see // whether we think we are likely past it, and switch off if we are. - pert_countdown -= abs(h); + pert_countdown -= std::abs(h); if (pert_countdown < 0.0) { pert = 0.0; pert_countdown = 0.0; } } } - return {x, true, "Path tracing terminated successfully."}; + return {x, true, "Path tracing terminated successfully.", steps}; +} + +PolishResult PolishPoint(std::function &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &x, double fixed_value, size_t fixed_index, + TerminationFunctionType p_terminate, int max_iter, + CallbackFunctionType p_callback) +{ + x[fixed_index] = fixed_value; + + const size_t N = x.size() - 1; + Vector y(N); // Equations results + Matrix jac_full(N + 1, N); // Full Jacobian matrix (N+1 unknowns, N equations) + Matrix jac_square(N, N); // Jacobian matrix with fixed_index row removed + Matrix Q(N, N); // Orthogonal matrix from QR decomposition + Vector x_reduced(N); // Reduced x vector with fixed_index removed + + int steps = 0; + double dist = 0.0; + + while (!p_terminate(x)) { + if (steps >= max_iter) { + return {x, false, "Polishing exceeded maximum iterations.", steps}; + } + + p_function(x, y); + p_jacobian(x, jac_full); + + size_t row_index = 1; + for (size_t i = 1; i <= N + 1; ++i) { // Newton step expects the transposed Jacobian + if (i != fixed_index) { + for (size_t j = 1; j <= N; ++j) { + jac_square(row_index, j) = jac_full(i, j); + } + row_index++; + } + } + + // Reduced x vector removing fixed_index + size_t temp_idx = 1; + for (size_t i = 1; i <= N + 1; ++i) { + if (i != fixed_index) { + x_reduced[temp_idx++] = x[i]; + } + } + + QRDecomp(jac_square, Q); + + // Solve jac_square * x_reduced = -y + NewtonStep(Q, jac_square, x_reduced, y, dist); + + // Update x, keeping fixed_index constant + temp_idx = 1; + for (size_t i = 1; i <= N + 1; ++i) { + if (i != fixed_index) { + x[i] = x_reduced[temp_idx++]; + } + } + + steps++; + + if (p_callback) { + p_callback(x); + } + } + + // The loop only exits here once p_terminate(x) holds for the actual, current x; + // there is nothing further to validate against a separate tolerance. + return {x, true, "Polishing terminated successfully.", steps}; } } // end namespace Gambit diff --git a/src/solvers/logit/path.h b/src/solvers/path/path.h similarity index 72% rename from src/solvers/logit/path.h rename to src/solvers/path/path.h index 5fa8e83a8..43cf9aff1 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 @@ -63,8 +63,15 @@ struct TracePathResult { Vector final_point; bool status; // true if path tracing terminated successfully, false if it terminated due to error std::string message; // error message if status is false + int steps; // Step at which the tracing terminated }; +struct PolishResult { + Vector final_point; + bool status; // true if polishing terminated successfully, false if it terminated due to error + std::string message; // error message if status is false + int steps; // Step at which the polishing terminated +}; // // This class implements a generic path-following algorithm for smooth curves. // It is based on the ideas and codes presented in Allgower and Georg's @@ -72,6 +79,7 @@ struct TracePathResult { // class PathTracer { public: + enum class TraceDirection { Positive = 1, Negative = -1 }; PathTracer() = default; virtual ~PathTracer() = default; @@ -84,7 +92,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, @@ -94,6 +103,15 @@ class PathTracer { double m_maxDecel{1.1}, m_hStart{0.03}; }; +// This function reduces the regret of a point that is close to an equilibrium that has been found +// by the path-following algorithm. Fixing the value of a component of the point, it uses a +// Newton-type method to find a nearby point with lower regret. +PolishResult PolishPoint(std::function &, Vector &)> p_function, + std::function &, Matrix &)> p_jacobian, + Vector &p_x, double fixed_value, size_t fixed_index, + TerminationFunctionType p_terminate, int max_iter = 100, + CallbackFunctionType p_callback = NullCallbackFunction); + } // end namespace Gambit #endif // GAMBIT_SOLVERS_LOGIT_PATH_H diff --git a/tests/cli/test_contract.py b/tests/cli/test_contract.py index f836ebc50..d1d9b4f68 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 000000000..5fa499102 --- /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 diff --git a/tests/games.py b/tests/games.py index 9b49eb2d3..657deae3b 100644 --- a/tests/games.py +++ b/tests/games.py @@ -227,6 +227,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_nash.py b/tests/test_nash.py index f8be93430..265affebf 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 ##################################################################################################