Grid refinement and leash length in simplicial subdivision.

This adds parameters to simpdiv_solve in pygambit to modify the grid refinement rate and leash length.

Internally, the leash length implementation was refactored in C++ for better clarity.

This closes #191.

ChangeLog

@@ -1,6 +1,13 @@
 # Changelog
 
-## [16.1.0a2] - unreleased
+## [16.1.0a3] - 2023-09-29
+
+### Changed
+- The `refine` and `leash` parameters to simpdiv have been made available and documented in
+  pygambit.
+
+
+## [16.1.0a2] - 2023-09-22
 
 ### Changed
 - Most operations which modify games have been moved to being operations on `Game` instead of

src/pygambit/gambit.pxd

@@ -399,7 +399,9 @@ cdef extern from "solvers/liap/liap.h":
     c_List[c_MixedBehaviorProfileDouble] LiapBehaviorSolve(c_Game, int p_maxitsN) except +RuntimeError
 
 cdef extern from "solvers/simpdiv/simpdiv.h":
-    c_List[c_MixedStrategyProfileRational] SimpdivStrategySolve(c_Game) except +RuntimeError
+    c_List[c_MixedStrategyProfileRational] SimpdivStrategySolve(c_Game,
+                                                                int p_gridResize,
+                                                                int p_leashLength) except +RuntimeError
 
 cdef extern from "solvers/ipa/ipa.h":
     c_List[c_MixedStrategyProfileDouble] IPAStrategySolve(c_Game) except +RuntimeError

src/pygambit/nash.pxi

@@ -123,8 +123,10 @@ def _liap_strategy_solve(game: Game, maxiter: int) -> typing.List[MixedStrategyP
 def _liap_behavior_solve(game: Game, maxiter: int) -> typing.List[MixedBehaviorProfileDouble]:
     return _convert_mbpd(LiapBehaviorSolve(game.game, maxiter))
 
-def _simpdiv_strategy_solve(game: Game) -> typing.List[MixedStrategyProfileRational]:
-    return _convert_mspr(SimpdivStrategySolve(game.game))
+def _simpdiv_strategy_solve(
+        game: Game, gridstep: int, leash: int
+) -> typing.List[MixedStrategyProfileRational]:
+    return _convert_mspr(SimpdivStrategySolve(game.game, gridstep, leash))
 
 def _ipa_strategy_solve(game: Game) -> typing.List[MixedStrategyProfileDouble]:
     return _convert_mspd(IPAStrategySolve(game.game))

src/pygambit/nash.py

@@ -500,6 +500,8 @@ def liap_solve(
 
 def simpdiv_solve(
         game: libgbt.Game,
+        refine: int = 2,
+        leash: typing.Optional[int] = None,
         external: bool = False
 ) -> typing.List[libgbt.MixedStrategyProfile]:
     """Compute Nash equilibria of a game using :ref:`simplicial
@@ -509,6 +511,15 @@ def simpdiv_solve(
     ----------
     game : Game
         The game to compute equilibria in.
+    refine : int, default 2
+        This controls the rate at which the triangulation of the space of mixed strategy
+        profiles is made more fine at each iteration.
+    leash : int, optional
+        Simplicial subdivision is guaranteed to converge to an (approximate) Nash equilibrium.
+        The method may take arbitrarily long paths through the space of mixed strategies in
+        doing so.  If specified, `leash` sets a maximum number of grid steps the method
+        may explore.  This trades off the possibility of finding an equilibrium more
+        quickly by giving up the guarantee than an equilibrium will necessarily be found.
     external : bool, default False
         Call the external command-line solver instead of the internally-linked
         implementation.  Requires the command-line solvers to be installed somewhere
@@ -521,7 +532,12 @@ def simpdiv_solve(
     """
     if external:
         return ExternalSimpdivSolver().solve(game)
-    return libgbt._simpdiv_strategy_solve(game)
+    if not isinstance(refine, int) or refine < 2:
+        raise ValueError("simpdiv_solve(): refine must be an integer no less than 2")
+    if leash is not None:
+        if not isinstance(leash, int) or leash <= 0:
+            raise ValueError("simpdiv_solve(): leash must be a non-negative integer")
+    return libgbt._simpdiv_strategy_solve(game, refine, leash or 0)
 
 
 def ipa_solve(

src/solvers/simpdiv/simpdiv.cc

@@ -35,12 +35,31 @@ inline GameStrategy GetStrategy(const Game &game, int pl, int st)
   return game->GetPlayer(pl)->Strategies()[st];
 }
 
+class NashSimpdivStrategySolver::State {
+  public:
+    int m_leashLength;
+    int t, ibar;
+    Rational d, pay, maxz, bestz;
+
+    State(int p_leashLength) : m_leashLength(p_leashLength), t(0), ibar(1), bestz(1.0e30) { }
+    Rational getlabel(MixedStrategyProfile<Rational> &yy, Array<int> &,
+		      PVector<Rational> &);
+
+    /* Check whether the distance p_dist is "too far" given the leash length, if set. */
+    bool CheckLeashOK(const Rational &p_dist) const {
+      if (m_leashLength == 0) {
+        return true;
+      }
+      return p_dist < m_leashLength * d;
+    }
+  };
+
 Rational 
 NashSimpdivStrategySolver::Simplex(MixedStrategyProfile<Rational> &y,
 				   const Rational &d) const
 {
   Game game = y.GetGame();
-  State state;
+  State state(m_leashLength);
   state.d = d;
   Array<int> nstrats(game->NumStrategies());
   Array<int> ylabel(2);
@@ -156,55 +175,55 @@ NashSimpdivStrategySolver::Simplex(MixedStrategyProfile<Rational> &y,
   getY(state, y, v, U, TT, ab, pi, ii);
 
   /* case3a */
-  if (i==1 && 
-      (y[GetStrategy(game, j, k)]<=Rational(0) || 
-       (v(j,k)-y[GetStrategy(game, j, k)]) >= Rational(m_leashLength)*state.d)) {
+  if (i == 1 &&
+      (y[GetStrategy(game, j, k)] <= Rational(0) ||
+       !state.CheckLeashOK(v(j, k) - y[GetStrategy(game, j, k)]))) {
     for (hh = 1, tot = 0; hh <= nstrats[j]; hh++) {
-      if (TT(j,hh)==1 || U(j,hh)==1)  {
-	tot++;
+      if (TT(j, hh) == 1 || U(j, hh) == 1) {
+        tot++;
       }
     }
     if (tot == nstrats[j] - 1) {
-      U(j,k)=1;
+      U(j, k) = 1;
       goto end;
     }
     else {
       update(state, pi, labels, ab, U, j, i);
-      U(j,k) = 1;
+      U(j, k) = 1;
       getnexty(state, y, pi, U, state.t);
       goto step1;
     }
   }
   /* case3b */
-  else if (i>=2 && i<=state.t &&
+  else if (i >= 2 && i <= state.t &&
 	   (y[GetStrategy(game, j, k)] <= Rational(0) || 
-	    (v(j,k)-y[GetStrategy(game, j, k)]) >= Rational(m_leashLength)*state.d)) {
+	    !state.CheckLeashOK(v(j,k)-y[GetStrategy(game, j, k)]))) {
     goto step4;
   }
   /* case3c */
-  else if (i==state.t+1 && ab(j,kk) == Rational(0)) {
-    if (y[GetStrategy(game, j, h)] <= Rational(0) || 
-	(v(j,h)-y[GetStrategy(game, j, h)]) >= Rational(m_leashLength)*state.d) {
+  else if (i == state.t + 1 && ab(j, kk) == Rational(0)) {
+    if (y[GetStrategy(game, j, h)] <= Rational(0) ||
+        !state.CheckLeashOK(v(j, h) - y[GetStrategy(game, j, h)])) {
       goto step4;
     }
     else {
-      k=0;
-      while (ab(j,kk) == Rational(0) && k==0) {
-	if(kk==h)k=1;
-	kk++;
-	if (kk > nstrats[j]) {
-	  kk=1; 
-	}
+      k = 0;
+      while (ab(j, kk) == Rational(0) && k == 0) {
+        if (kk == h)k = 1;
+        kk++;
+        if (kk > nstrats[j]) {
+          kk = 1;
+        }
       }
       kk--;
       if (kk == 0) {
-	kk = nstrats[j];
+        kk = nstrats[j];
       }
       if (kk == h) {
-	goto step4;
+        goto step4;
       }
       else {
-	goto step5;
+        goto step5;
       }
     }
   }
@@ -232,7 +251,7 @@ NashSimpdivStrategySolver::Simplex(MixedStrategyProfile<Rational> &y,
   h = pi(i-1,2);
   TT(j,h) = 0;
   if (y[GetStrategy(game, j, h)] <= Rational(0) || 
-      (v(j,h)-y[GetStrategy(game, j, h)]) >= Rational(m_leashLength)*state.d) {
+      !state.CheckLeashOK(v(j,h)-y[GetStrategy(game, j, h)])) {
     U(j,h) = 1;
   }
   labels.RotateUp(i,state.t+1);

src/solvers/simpdiv/simpdiv.h

@@ -53,15 +53,7 @@ class NashSimpdivStrategySolver : public StrategySolver<Rational> {
   int m_gridResize, m_leashLength;
   bool m_verbose;
 
-  class State {
-  public:
-    int t, ibar;
-    Rational d, pay, maxz, bestz;
-
-    State() : t(0), ibar(1), bestz(1.0e30) { }
-    Rational getlabel(MixedStrategyProfile<Rational> &yy, Array<int> &, 
-		      PVector<Rational> &);
-  };
+  class State;
 
   Rational Simplex(MixedStrategyProfile<Rational> &, const Rational &d) const;
   void update(State &, RectArray<int> &, RectArray<int> &, PVector<Rational> &,
@@ -76,9 +68,10 @@ class NashSimpdivStrategySolver : public StrategySolver<Rational> {
   int get_b(int j, int h, int nstrats, const PVector<int> &) const;
 };
 
-inline List<MixedStrategyProfile<Rational> > SimpdivStrategySolve(const Game &p_game)
+inline List<MixedStrategyProfile<Rational> >
+SimpdivStrategySolve(const Game &p_game, int p_gridResize = 2, int p_leashLength = 0)
 {
-  return NashSimpdivStrategySolver().Solve(p_game);
+  return NashSimpdivStrategySolver(p_gridResize, p_leashLength).Solve(p_game);
 }
 
 }  // end namespace Gambit::Nash