Add Herings-Peeters homotopy method

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}
 ## 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)

doc/pygambit.api.rst

@@ -403,6 +403,7 @@ Computation of Nash equilibria
    simpdiv_solve
    ipa_solve
    gnm_solve
+   hp_solve
 
 
 Computation of quantal response equilibria

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},

doc/tools.rst

@@ -53,3 +53,4 @@ documentation.
    tools.logit
    tools.gnm
    tools.ipa
+   tools.hp

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)",

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,

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<double>,
-                       p_onEvent, p_cancel);
+    LogitBehaviorSolve(start, 1.0e-8, PathTracer::TraceDirection::Positive, 0.03, 1.1,
+                       Nash::NullBehaviorCallback<double>, 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<double>,
-                       p_onEvent, p_cancel);
+    LogitStrategySolve(start, 1.0e-8, PathTracer::TraceDirection::Positive, 0.03, 1.1,
+                       Nash::NullStrategyCallback<double>, p_onEvent, p_cancel);
   }
 };
 

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 <class M>
 concept StrategicMethod =
     std::same_as<M, EnumPureNashSpec> || std::same_as<M, EnumMixedNashSpec> ||
-    std::same_as<M, GNMNashSpec> || std::same_as<M, IPANashSpec> ||
+    std::same_as<M, GNMNashSpec> || std::same_as<M, HPNashSpec> || std::same_as<M, IPANashSpec> ||
     std::same_as<M, LiapNashSpec> || std::same_as<M, SimpdivNashSpec>;
 
 template <class M>
@@ -170,6 +174,9 @@ wxString ExternalCommand(const NashComputationSpec &p_spec)
                                            method.localNewtonInterval,
                                            method.localNewtonMaxIterations);
         }
+        else if constexpr (std::is_same_v<Method, HPNashSpec>) {
+          return prefix + wxString::Format("hp -d 10 -n %d", method.priors);
+        }
         else if constexpr (std::is_same_v<Method, IPANashSpec>) {
           return prefix + wxString::Format("ipa -d 10 -n %d", method.perturbations);
         }
@@ -223,6 +230,9 @@ wxString MethodDescription(const NashMethodSpec &p_method)
         else if constexpr (std::is_same_v<Method, GNMNashSpec>) {
           return wxT("by global Newton tracing");
         }
+        else if constexpr (std::is_same_v<Method, HPNashSpec>) {
+          return wxT("by the Herings-Peeters homotopy");
+        }
         else if constexpr (std::is_same_v<Method, IPANashSpec>) {
           return wxT("by iterated polymatrix approximation");
         }
@@ -263,6 +273,9 @@ wxString ParameterDescription(const NashMethodSpec &p_method)
               method.perturbations, method.lambdaEnd, method.steps, method.localNewtonInterval,
               method.localNewtonMaxIterations);
         }
+        else if constexpr (std::is_same_v<Method, HPNashSpec>) {
+          return wxString::Format(" (%d random priors)", method.priors);
+        }
         else if constexpr (std::is_same_v<Method, IPANashSpec>) {
           return wxString::Format(" (%d perturbation)", method.perturbations);
         }
@@ -387,6 +400,7 @@ void NashChoiceDialog::OnCount(wxCommandEvent &p_event)
     m_methodChoice->Append(s_liap);
     m_methodChoice->Append(s_gnm);
     m_methodChoice->Append(s_ipa);
+    m_methodChoice->Append(s_hp);
     m_methodChoice->Append(s_enumpoly);
   }
   else {

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<SolverFunction> GNMNashSpec::MakeSolver(NashRepresentation) const
   };
 }
 
+std::optional<SolverFunction> HPNashSpec::MakeSolver(NashRepresentation) const
+{
+  const HPNashSpec spec = *this;
+  return [spec](const Game &p_game, const ProfileFoundCallback &p_callback,
+                const CancelToken &p_cancel) {
+    for (const auto &prior : NewRandomStrategyProfiles(p_game, spec.priors)) {
+      p_cancel.Check();
+      Nash::HPStrategySolve(
+          prior,
+          [&p_callback](const MixedStrategyProfile<double> &p) { p_callback(ComputedProfile(p)); },
+          Nash::NullHPEventCallback, p_cancel);
+    }
+  };
+}
+
 std::optional<SolverFunction> IPANashSpec::MakeSolver(NashRepresentation) const
 {
   const IPANashSpec spec = *this;
@@ -178,7 +194,7 @@ std::optional<SolverFunction> 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<double> &p) { p_callback(ComputedProfile(p)); },
           NullLogitEventCallback<LogitQREMixedBehaviorProfile>, p_cancel);
     };
@@ -187,7 +203,7 @@ std::optional<SolverFunction> 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<double> &p) { p_callback(ComputedProfile(p)); },
         NullLogitEventCallback<LogitQREMixedStrategyProfile>, p_cancel);
   };

src/gui/nashspec.h

@@ -27,8 +27,11 @@
 #include <variant>
 
 #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,15 @@ struct IPANashSpec {
   std::optional<SolverFunction> MakeSolver(NashRepresentation) const;
 };
 
+struct HPNashSpec {
+  // Unlike GNM/IPA, a single prior yields at most one equilibrium (no internal
+  // path-tracing can surface more), so this defaults high, like LiapNashSpec's
+  // startingPoints, rather than to 1.
+  int priors{10};
+
+  std::optional<SolverFunction> MakeSolver(NashRepresentation) const;
+};
+
 struct LPNashSpec {
   std::optional<SolverFunction> MakeSolver(NashRepresentation p_representation) const;
 };
@@ -115,7 +127,7 @@ struct LiapNashSpec {
 
 struct LogitNashSpec {
   double maxRegret{1.0e-8};
-  double omega{1.0};
+  PathTracer::TraceDirection direction{PathTracer::TraceDirection::Positive};
   double firstStep{0.03};
   double maxAcceleration{1.1};
 
@@ -132,9 +144,9 @@ struct SimpdivNashSpec {
   std::optional<SolverFunction> MakeSolver(NashRepresentation) const;
 };
 
-using NashMethodSpec =
-    std::variant<EnumPureNashSpec, EnumMixedNashSpec, EnumPolyNashSpec, GNMNashSpec, IPANashSpec,
-                 LPNashSpec, LCPNashSpec, LiapNashSpec, LogitNashSpec, SimpdivNashSpec>;
+using NashMethodSpec = std::variant<EnumPureNashSpec, EnumMixedNashSpec, EnumPolyNashSpec,
+                                    GNMNashSpec, HPNashSpec, IPANashSpec, LPNashSpec, LCPNashSpec,
+                                    LiapNashSpec, LogitNashSpec, SimpdivNashSpec>;
 
 struct NashComputationSpec {
   NashRepresentation representation;

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<Gambit::LogitQREMixedBehaviorProfile> p_qre);
+std::string
+InvokeHPStrategyEventCallback(PyObject *p_callback,
+                              std::shared_ptr<Gambit::MixedStrategyProfile<double>> p_profile,
+                              double p_t);
 std::string InvokeGNMPerturbationEventCallback(
     PyObject *p_callback, std::shared_ptr<Gambit::MixedStrategyProfile<double>> p_profile);
 std::string
@@ -228,6 +233,25 @@ MakeLogitEventCallback<Gambit::LogitQREMixedBehaviorProfile>(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<Gambit::Nash::HPStepEvent>(p_event);
+    Gambit::ThrowIfPythonError(InvokeHPStrategyEventCallback(
+        p_callback, std::make_shared<Gambit::MixedStrategyProfile<double>>(step.profile), step.t));
+  };
+}
+
 ///
 /// Builds a Nash::GNMEventCallbackType which, when invoked, dispatches to
 /// whichever Invoke*EventCallback trampoline matches the alternative held by