445 enh add starting point for simpdiv solve

ChangeLog

@@ -33,6 +33,7 @@
   to always take an explicit starting point (in `liap_solve`), and to specify a regret-based
   acceptance criterion (#330)
 - `simpdiv_solve`/`gambit-simpdiv` now accepts a regret-based acceptance criterion (#439)
+- `simpdiv_solve` now takes an explicit starting point (#445)
 - Converted test suite for mixed behavior profiles to pytest style; added parametrizations for
   test_realiz_prob; added test_martingale_property_of_node_value (#375)
 - Improved test suite for mixed strategy profiles (#374)

doc/pygambit.user.rst

@@ -601,6 +601,90 @@ If, instead, we double all payoffs, the output of the method is unchanged.
    result.equilibria[0].max_regret() / (g.max_payoff - g.min_payoff)
 
 
+Generating starting points for algorithms
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+Some methods for computation of Nash equilibria take as an initial condition a
+:py:class:`.MixedStrategyProfile` or :py:class:`MixedBehaviorProfile` which is used
+as a starting point.  The equilibria found will depend on which starting point is
+selected.  To facilitate generating starting points, :py:class:`.Game` provides
+methods :py:meth:`.Game.random_strategy_profile` and :py:meth:`.Game.random_behavior_profile`,
+to generate profiles which are drawn from the uniform distribution on the product
+of simplices.
+
+As an example, we consider a three-player game from McKelvey and McLennan (1997),
+in which each player has two strategies.  This game has nine equilibria in total, and
+in particular has two totally mixed Nash equilibria, which is the maximum possible number
+of regular totally mixed equilbria in games of this size.
+
+We first consider finding Nash equilibria in this game using :py:func:`.liap_solve`.
+If we run this method starting from the centroid (uniform randomization across all
+strategies for each player), :py:func:`.liap_solve` finds one of the totally-mixed equilibria.
+
+.. ipython:: python
+
+   g = gbt.Game.read_game("2x2x2.nfg")
+   gbt.nash.liap_solve(g.mixed_strategy_profile())
+
+Which equilibrium is found depends on the starting point.  With a different starting point,
+we can find, for example, one of the pure-strategy equilibria.
+
+.. ipython:: python
+
+   gbt.nash.liap_solve(g.mixed_strategy_profile([[.9, .1], [.9, .1], [.9, .1]]))
+
+To search for more equilibria, we can instead generate strategy profiles at random.
+
+.. ipython:: python
+
+   gbt.nash.liap_solve(g.random_strategy_profile())
+
+Note that methods which take starting points do record the starting points used in the
+result object returned.  However, the random profiles which are generated will differ
+in different runs of a program.  To support making the generation of random strategy
+profiles reproducible, and for finer-grained control of the generation of these profiles
+if desired, :py:meth:`.Game.random_strategy_profile` and :py:meth:`.Game.random_behavior_profile`
+optionally take a :py:class:`numpy.random.Generator` object, which is used as the source
+of randomness for creating the profile.
+
+.. ipython:: python
+
+   import numpy as np
+   gen = np.random.default_rng(seed=1234567890)
+   p1 = g.random_strategy_profile(gen=gen)
+   p1
+   gen = np.random.default_rng(seed=1234567890)
+   p2 = g.random_strategy_profile(gen=gen)
+   p2
+   p1 == p2
+
+When creating profiles in which probabilities are represented as floating-point numbers,
+:py:meth:`.Game.random_strategy_profile` and :py:meth:`.Game.random_behavior_profile`
+internally use the Dirichlet distribution for each simplex to generate correctly uniform
+sampling over probabilities.  However, in some applications generation of random profiles
+with probabilities as rational numbers is desired.  For example, :py:func:`.simpdiv_solve`
+takes such a starting point, because it operates by successively refining a triangulation
+over the space of mixed strategy profiles.
+:py:meth:`.Game.random_strategy_profile` and :py:meth:`.Game.random_behavior_profile`
+both take an optional parameter `denom` which, if specified, generates a profile in which
+probabilities are generated uniformly from the grid in each simplex in which all probabilities
+have denominator `denom`.
+
+.. ipython:: python
+
+   gen = np.random.default_rng(seed=1234567890)
+   g.random_strategy_profile(denom=10, gen=gen)
+   g.random_strategy_profile(denom=10, gen=gen)
+
+These can then be used in conjunction with :py:func:`.simpdiv_solve` to search for equilibria
+from different starting points.
+
+.. ipython:: python
+
+   gbt.nash.simpdiv_solve(g.random_strategy_profile(denom=10, gen=gen))
+   gbt.nash.simpdiv_solve(g.random_strategy_profile(denom=10, gen=gen))
+   gbt.nash.simpdiv_solve(g.random_strategy_profile(denom=10, gen=gen))
+
 
 Estimating quantal response equilibria
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

src/pygambit/gambit.pxd

@@ -440,7 +440,8 @@ cdef extern from "solvers/liap/liap.h":
 
 cdef extern from "solvers/simpdiv/simpdiv.h":
     c_List[c_MixedStrategyProfileRational] SimpdivStrategySolve(
-            c_Game, c_Rational p_maxregret, int p_gridResize, int p_leashLength
+            c_MixedStrategyProfileRational start, c_Rational p_maxregret, int p_gridResize,
+            int p_leashLength
     ) except +RuntimeError
 
 cdef extern from "solvers/ipa/ipa.h":

src/pygambit/nash.pxi

@@ -149,9 +149,9 @@ def _liap_behavior_solve(start: MixedBehaviorProfileDouble,
 
 
 def _simpdiv_strategy_solve(
-        game: Game, maxregret: Rational, gridstep: int, leash: int
+        start: MixedStrategyProfileRational, maxregret: Rational, gridstep: int, leash: int
 ) -> typing.List[MixedStrategyProfileRational]:
-    return _convert_mspr(SimpdivStrategySolve(game.game,
+    return _convert_mspr(SimpdivStrategySolve(deref(start.profile),
                                               to_rational(str(maxregret).encode("ascii")),
                                               gridstep, leash))
 

src/pygambit/nash.py

@@ -264,7 +264,9 @@ def liap_solve(
     Parameters
     ----------
     start : MixedStrategyProfileDouble or MixedBehaviorProfileDouble
-        The starting point for function minimization.
+        The starting profile for function minimization.  Up to one equilibrium will be found
+        from any starting profile, and the equilibrium found may (and generally will)
+        depend on the initial profile chosen.
 
     maxregret : float, default 1e-4
         The acceptance criterion for approximate Nash equilibrium; the maximum
@@ -307,18 +309,24 @@ def liap_solve(
 
 
 def simpdiv_solve(
-        game: libgbt.Game,
+        start: libgbt.MixedStrategyProfileRational,
         maxregret: libgbt.Rational = None,
         refine: int = 2,
         leash: typing.Optional[int] = None
 ) -> NashComputationResult:
     """Compute Nash equilibria of a game using :ref:`simplicial
     subdivision <gambit-simpdiv>`.
 
+    .. versionchanged:: 16.2.0
+
+       Method now takes a starting point, as a mixed strategy profile, instead of a game.
+
     Parameters
     ----------
-    game : Game
-        The game to compute equilibria in.
+    start: MixedStrategyProfileRational
+        The starting profile for the algorithm.  Up to one equilibrium will be found
+        from any starting profile, and the equilibrium found may (and generally will)
+        depend on the initial profile chosen.
 
     maxregret : Rational, default 1e-8
         The acceptance criterion for approximate Nash equilibrium; the maximum
@@ -351,14 +359,14 @@ def simpdiv_solve(
         maxregret = libgbt.Rational(1, 10000000)
     elif maxregret < libgbt.Rational(0):
         raise ValueError("simpdiv_solve(): maxregret must be positive")
-    equilibria = libgbt._simpdiv_strategy_solve(game, maxregret, refine, leash or 0)
+    equilibria = libgbt._simpdiv_strategy_solve(start, maxregret, refine, leash or 0)
     return NashComputationResult(
-        game=game,
+        game=start.game,
         method="simpdiv",
         rational=True,
         use_strategic=True,
         equilibria=equilibria,
-        parameters={"maxregret": maxregret, "refine": refine, "leash": leash}
+        parameters={"start": start, "maxregret": maxregret, "refine": refine, "leash": leash}
     )
 
 

src/solvers/simpdiv/simpdiv.h

@@ -31,7 +31,7 @@ namespace Nash {
 ///
 /// This is a simplicial subdivision algorithm with restart, for finding
 /// mixed strategy solutions to general finite n-person games.  It is based on
-/// van Der Laan, Talman and van Der Heyden, Math in Oper Res, 1987.
+/// van Der Laan, Talman and van Der Heyden, Math of Oper Res, 1987.
 ///
 class NashSimpdivStrategySolver : public StrategySolver<Rational> {
 public:
@@ -68,10 +68,11 @@ class NashSimpdivStrategySolver : public StrategySolver<Rational> {
 };
 
 inline List<MixedStrategyProfile<Rational>>
-SimpdivStrategySolve(const Game &p_game, const Rational &p_maxregret = Rational(1, 1000000),
-                     int p_gridResize = 2, int p_leashLength = 0)
+SimpdivStrategySolve(const MixedStrategyProfile<Rational> &p_start,
+                     const Rational &p_maxregret = Rational(1, 1000000), int p_gridResize = 2,
+                     int p_leashLength = 0)
 {
-  return NashSimpdivStrategySolver(p_gridResize, p_leashLength, p_maxregret).Solve(p_game);
+  return NashSimpdivStrategySolver(p_gridResize, p_leashLength, p_maxregret).Solve(p_start);
 }
 
 } // namespace Nash

tests/test_nash.py

@@ -108,7 +108,7 @@ def test_liap_behavior(self):
 
     def test_simpdiv_strategy(self):
         """Test calls of simplicial subdivision for mixed strategy equilibria."""
-        result = gbt.nash.simpdiv_solve(self.poker)
+        result = gbt.nash.simpdiv_solve(self.poker.mixed_strategy_profile(rational=True))
         assert len(result.equilibria) == 1
 
     def test_ipa_strategy(self):