diff --git a/ChangeLog b/ChangeLog index 15c785f8c0..b3528cba0f 100644 --- a/ChangeLog +++ b/ChangeLog @@ -8,6 +8,8 @@ - Added `-R` option to `gambit-liap`, `gambit-simpdiv`, `gambit-ipa`, and `gambit-gnm` to reproducibly seed random starting point generation. (#448) - Added `Game.relabel_players`, which simultaneously reassigns the labels of the game's players. (#1058) +- Added `max_rectangles` to `enumpoly_solve` (and `-r` to `gambit-enumpoly`), bounding the number + of cells examined when searching for equilibria on a single support. (#1055) ### Fixed - `Game.add_player`, on a game in strategic form, no longer leaves the new player's sole @@ -17,6 +19,7 @@ - `enumpoly_solve` correctly handles failed perturbations at singular Jacobians (#591) - `enumpoly_solve` iterates sequences of actions in a deterministic order, pre-empting the possibility of equilibria being found in orders that differ across runs. (#1016) +- Improved the numerical stability of `enumpoly_solve` handling less-well-conditioned systems. (#1055) ### Changed - All long-running calculations are now done in a thread in the GUI process, rather than calling the diff --git a/Makefile.am b/Makefile.am index 5e93fa0a27..c1a8af428a 100644 --- a/Makefile.am +++ b/Makefile.am @@ -352,6 +352,8 @@ enumpoly_SOURCES = \ ${nashsupport_SOURCES} \ src/solvers/enumpoly/indexproduct.h \ src/solvers/enumpoly/rectangle.h \ + src/solvers/enumpoly/lufactor.cc \ + src/solvers/enumpoly/lufactor.h \ src/solvers/enumpoly/poly.cc \ src/solvers/enumpoly/poly.h \ src/solvers/enumpoly/polysystem.h \ diff --git a/doc/algorithms.rst b/doc/algorithms.rst index 260ea688db..4238a15559 100644 --- a/doc/algorithms.rst +++ b/doc/algorithms.rst @@ -81,6 +81,17 @@ Subsequent analysis of unreached information sets can yield alternative profiles that specify different choices at unreached information sets while still satisfying the Nash equilibrium conditions. +The subdivision process for a given support is guaranteed to terminate, because +subdivision stops once a cell shrinks below a minimum size. However, this bound +is on the *depth* of subdivision, not on the *number* of cells examined to reach +that depth. As protection against pathological cases, the number of cells +examined per support is capped (via ``max_rectangles`` in +:py:func:`pygambit.nash.enumpoly_solve`, or :option:`gambit-enumpoly -r` for +:program:`gambit-enumpoly`); if this limit is reached before the search of a +support is complete, that support may be reported with incomplete equilibria, +and (in verbose mode) is identified with the label "budget-exceeded" -- see +:ref:`gambit-enumpoly ` for details. + .. _lp: lp diff --git a/doc/tools.enumpoly.rst b/doc/tools.enumpoly.rst index 78042c9706..91996b2a91 100644 --- a/doc/tools.enumpoly.rst +++ b/doc/tools.enumpoly.rst @@ -23,6 +23,15 @@ share the listed support. However, the converse is not true: not all supports labeled as "singular" will necessarily be the support of some set of equilibria. +Independently of degeneracy, the search on a single support can require +examining an impractically large number of subdivided cells before it +can either confirm or rule out an equilibrium. To bound the time spent +on any one support, the number of cells examined is capped (see +:option:`gambit-enumpoly -r`). When this limit is reached, the search for that support +is abandoned, any equilibria already confirmed on it are still reported, +and (in verbose mode) the support is identified with the label +"budget-exceeded,". + .. program:: gambit-enumpoly .. cmdoption:: -d @@ -55,6 +64,14 @@ support of some set of equilibria. By default, the program will search all support profiles. This switch instructs the program to terminate when EQA equilibria have been found. +.. cmdoption:: -r RECTANGLES + + .. versionadded:: 17.0.0 + + Specify the maximum number of cells to examine when searching for equilibria + on a single support, before giving up on that support (default 20000). See + the :ref:`algorithm description ` for why this is necessary. + .. cmdoption:: -q Suppresses printing of the banner at program launch. diff --git a/src/core/matrix.h b/src/core/matrix.h index 870bbf2120..89b95563c1 100644 --- a/src/core/matrix.h +++ b/src/core/matrix.h @@ -41,10 +41,6 @@ class SingularMatrixException final : public std::runtime_error { /// /// @tparam T Scalar element type. Must support arithmetic operations, /// comparison with zero, and abs(T). -/// -/// @note Inverse and Determinant use legacy Gaussian elimination algorithms. -/// They are known not to be numerically optimal for near-singular matrices. -/// Current behaviour is temporarily preserved for historical compatibility. template class Matrix { RectArray m_data; @@ -205,8 +201,6 @@ template class Matrix { /// @name Other operations ///@{ Matrix Transpose() const; - Matrix Inverse() const; - T Determinant() const; ///@} }; @@ -412,158 +406,6 @@ template Matrix Matrix::Transpose() const return tmp; } -// ---------------------------------------------------------------------------- -// Implementation of additional operations -// ---------------------------------------------------------------------------- - -template Matrix Matrix::Inverse() const -{ - if (!IsSquare()) { - throw DimensionException(); - } - const int rmin = MinRow(); - const int rmax = MaxRow(); - const int cmin = MinCol(); - const int cmax = MaxCol(); - using Gambit::abs; - - Matrix copy(*this); - Matrix inv(rmin, rmax, cmin, cmax); - - inv = T{0}; - - // initialize inverse matrix and prescale row vectors - for (int i = rmin; i <= rmax; ++i) { - auto copy_row = copy.m_data.GetRowView(i); - auto inv_row = inv.m_data.GetRowView(i); - - T max = maximize_function(copy_row, [](const T &v) { return abs(v); }); - if (max == T{0}) { - throw SingularMatrixException(); - } - const T scale = T{1} / max; - for (auto &v : copy_row) { - v *= scale; - } - inv_row[i] = scale; - } - - for (int i = cmin; i <= cmax; ++i) { - // find pivot row - auto col_i = copy.m_data.GetColumnView(i); - T max = abs(col_i[i]); - int row = i; - for (int j = i + 1; j <= rmax; ++j) { - const T v = abs(col_i[j]); - if (v > max) { - max = v; - row = j; - } - } - - if (max <= T{0}) { - throw SingularMatrixException(); - } - - copy.SwitchRows(i, row); - inv.SwitchRows(i, row); - // scale pivot row - T factor = T{1} / copy(i, i); - auto copy_row = copy.m_data.GetRowView(i); - auto inv_row = inv.m_data.GetRowView(i); - auto copy_it = copy_row.begin(); - auto inv_it = inv_row.begin(); - for (; copy_it != copy_row.end(); ++copy_it, ++inv_it) { - *copy_it *= factor; - *inv_it *= factor; - } - - // reduce other rows - auto pivot_copy_row = copy.m_data.GetRowView(i); - auto pivot_inv_row = inv.m_data.GetRowView(i); - - for (int j = rmin; j <= rmax; ++j) { - if (j == i) { - continue; - } - auto row_copy = copy.m_data.GetRowView(j); - auto row_inv = inv.m_data.GetRowView(j); - const T mult = row_copy[i]; - auto pivot_copy_it = pivot_copy_row.begin(); - auto pivot_inv_it = pivot_inv_row.begin(); - auto row_copy_it = row_copy.begin(); - auto row_inv_it = row_inv.begin(); - - for (; pivot_copy_it != pivot_copy_row.end(); - ++pivot_copy_it, ++pivot_inv_it, ++row_copy_it, ++row_inv_it) { - *row_copy_it -= (*pivot_copy_it) * mult; - *row_inv_it -= (*pivot_inv_it) * mult; - } - } - } - - return inv; -} - -template T Matrix::Determinant() const -{ - if (!IsSquare()) { - throw DimensionException(); - } - const int rmin = MinRow(); - const int rmax = MaxRow(); - using Gambit::abs; - - Matrix M(*this); - - for (int row = rmin; row <= rmax; ++row) { - // Experience (as of 3/22/99) suggests that, in the interest of - // numerical stability, it might be best to do Gaussian - // elimination with respect to the row (of those feasible) - // whose entry has the largest absolute value. - int swap_row = row; - T max = abs(M(row, row)); - for (int i = row + 1; i <= rmax; ++i) { - const T v = abs(M(i, row)); - if (v > max) { - max = v; - swap_row = i; - } - } - - if (swap_row != row) { - M.SwitchRows(row, swap_row); - auto pivot_row = M.m_data.GetRowView(row); - for (auto &v : pivot_row) { - v = -v; - } - } - - if (M(row, row) == T{0}) { - return T{0}; - } - // now do row operations to clear the row'th column - // below the diagonal - auto pivot_row = M.m_data.GetRowView(row); - for (int row1 = row + 1; row1 <= rmax; ++row1) { - auto elim_row = M.m_data.GetRowView(row1); - const T factor = -elim_row[row] / pivot_row[row]; - auto pivot_it = pivot_row.begin(); - auto elim_it = elim_row.begin(); - for (; pivot_it != pivot_row.end(); ++pivot_it, ++elim_it) { - *elim_it += (*pivot_it) * factor; - } - } - } - - // finally we multiply the diagonal elements - T det = T{1}; - for (int row = rmin; row <= rmax; ++row) { - det *= M(row, row); - } - return det; -} - // ---------------------------------------------------------------------------- // Implementation of operators // ---------------------------------------------------------------------------- diff --git a/src/gui/nashspec.cc b/src/gui/nashspec.cc index 36130495e6..4c8ee29a23 100644 --- a/src/gui/nashspec.cc +++ b/src/gui/nashspec.cc @@ -65,7 +65,7 @@ EnumPolyNashSpec::MakeSolver(NashRepresentation p_representation) const return [spec](const Game &p_game, const ProfileFoundCallback &p_callback, const CancelToken &p_cancel) { Nash::EnumPolyBehaviorSolve( - p_game, spec.stopAfter, spec.maxRegret, + p_game, spec.stopAfter, spec.maxRegret, Nash::kDefaultEnumPolyMaxRectangles, [&p_callback](const MixedBehaviorProfile &p) { p_callback(ComputedProfile(p)); }, Nash::NullEnumPolyEventCallback, p_cancel); }; @@ -73,7 +73,7 @@ EnumPolyNashSpec::MakeSolver(NashRepresentation p_representation) const return [spec](const Game &p_game, const ProfileFoundCallback &p_callback, const CancelToken &p_cancel) { Nash::EnumPolyStrategySolve( - p_game, spec.stopAfter, spec.maxRegret, + p_game, spec.stopAfter, spec.maxRegret, Nash::kDefaultEnumPolyMaxRectangles, [&p_callback](const MixedStrategyProfile &p) { p_callback(ComputedProfile(p)); }, Nash::NullEnumPolyEventCallback, p_cancel); }; diff --git a/src/pygambit/gambit.pxd b/src/pygambit/gambit.pxd index 785008659b..29ece38f0c 100644 --- a/src/pygambit/gambit.pxd +++ b/src/pygambit/gambit.pxd @@ -594,10 +594,10 @@ cdef extern from "solvers/nashsupport/nashsupport.h": cdef extern from "solvers/enumpoly/enumpoly.h": stdlist[c_MixedStrategyProfile[double]] EnumPolyStrategySolve( - c_Game, int, float + c_Game, int, float, size_t ) except +RuntimeError stdlist[c_MixedBehaviorProfile[double]] EnumPolyBehaviorSolve( - c_Game, int, float + c_Game, int, float, size_t ) except +RuntimeError cdef extern from "solvers/logit/logit.h": diff --git a/src/pygambit/nash.pxi b/src/pygambit/nash.pxi index 4d5a547221..64fffa696d 100644 --- a/src/pygambit/nash.pxi +++ b/src/pygambit/nash.pxi @@ -180,16 +180,18 @@ def _enumpoly_strategy_solve( game: Game, stop_after: int, maxregret: float, + max_rectangles: int, ) -> list[MixedStrategyProfileDouble]: - return _convert_mspd(EnumPolyStrategySolve(game.game, stop_after, maxregret)) + return _convert_mspd(EnumPolyStrategySolve(game.game, stop_after, maxregret, max_rectangles)) def _enumpoly_behavior_solve( game: Game, stop_after: int, maxregret: float, + max_rectangles: int, ) -> list[MixedBehaviorProfileDouble]: - return _convert_mbpd(EnumPolyBehaviorSolve(game.game, stop_after, maxregret)) + return _convert_mbpd(EnumPolyBehaviorSolve(game.game, stop_after, maxregret, max_rectangles)) def _logit_strategy_solve( diff --git a/src/pygambit/nash.py b/src/pygambit/nash.py index 21d8f52b37..11639dad6d 100644 --- a/src/pygambit/nash.py +++ b/src/pygambit/nash.py @@ -662,7 +662,8 @@ def enumpoly_solve( game: libgbt.Game, use_strategic: bool = False, stop_after: int | None = None, - maxregret: float = 1.0e-4, + maxregret: float = 1.0e-8, + max_rectangles: int = 20_000, phcpack_path: pathlib.Path | str | None = None ) -> NashComputationResult: """:ref:`Compute Nash equilibria by enumerating all support profiles @@ -682,11 +683,18 @@ def enumpoly_solve( Maximum number of equilibria to compute. If not specified, examines all support profiles of the game. - maxregret : float, default 1e-4 + 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 + max_rectangles : int, default 20000 + The maximum number of rectangles to examine when searching for equilibria on a + single support, before giving up on that support. See the + :ref:`algorithm description ` for why this is necessary. + + .. versionadded:: 17.0.0 + phcpack_path : str or pathlib.Path, optional If specified, use PHCpack to solve the systems of equations. This argument specifies the path to the PHCpack executable. @@ -712,6 +720,11 @@ def enumpoly_solve( raise ValueError( f"enumpoly_solve(): maxregret must be a positive number; got {maxregret}" ) + if max_rectangles <= 0: + raise ValueError( + f"enumpoly_solve(): " + f"max_rectangles argument must be a positive number; got {max_rectangles}" + ) if phcpack_path is not None: if game.is_tree and not use_strategic: raise ValueError( @@ -730,15 +743,16 @@ def enumpoly_solve( ) if not game.is_tree or use_strategic: - equilibria = libgbt._enumpoly_strategy_solve(game, stop_after, maxregret) + equilibria = libgbt._enumpoly_strategy_solve(game, stop_after, maxregret, max_rectangles) else: - equilibria = libgbt._enumpoly_behavior_solve(game, stop_after, maxregret) + equilibria = libgbt._enumpoly_behavior_solve(game, stop_after, maxregret, max_rectangles) return NashComputationResult( game=game, method="enumpoly", rational=False, use_strategic=not game.is_tree or use_strategic, - parameters={"stop_after": stop_after, "maxregret": maxregret}, + parameters={"stop_after": stop_after, "maxregret": maxregret, + "max_rectangles": max_rectangles}, equilibria=equilibria, ) diff --git a/src/solvers/enumpoly/efgpoly.cc b/src/solvers/enumpoly/efgpoly.cc index c7d81383a7..b8b51645c4 100644 --- a/src/solvers/enumpoly/efgpoly.cc +++ b/src/solvers/enumpoly/efgpoly.cc @@ -182,8 +182,9 @@ FindNashExtension(const MixedBehaviorProfile &p_baseProfile, double p_ma } std::list> SolveSupport(const BehaviorSupportProfile &p_support, - bool &p_isSingular, int p_stopAfter, - double p_maxRegret, + bool &p_isSingular, bool &p_budgetExceeded, + int p_stopAfter, double p_maxRegret, + size_t p_maxRectangles, const CancelToken &p_cancel) { ProblemData data(p_support); @@ -201,7 +202,7 @@ std::list> SolveSupport(const BehaviorSupportProfil try { roots = solver.FindRoots({bottoms, tops}, (p_stopAfter > 0) ? p_stopAfter : std::numeric_limits::max(), - p_cancel); + p_maxRectangles, p_budgetExceeded, p_cancel); } catch (const SingularMatrixException &) { p_isSingular = true; @@ -228,7 +229,7 @@ namespace Gambit::Nash { std::list> EnumPolyBehaviorSolve(const Game &p_game, int p_stopAfter, double p_maxregret, - BehaviorCallbackType p_onEquilibrium, + size_t p_maxRectangles, BehaviorCallbackType p_onEquilibrium, EnumPolyEventCallbackType p_onEvent, const CancelToken &p_cancel) { @@ -249,15 +250,20 @@ EnumPolyBehaviorSolve(const Game &p_game, int p_stopAfter, double p_maxregret, p_cancel.Check(); p_onEvent(EnumPolyCandidateSupportEvent{support}); bool isSingular = false; + bool budgetExceeded = false; for (const auto &solution : - SolveSupport(support, isSingular, std::max(p_stopAfter - static_cast(ret.size()), 0), - p_maxregret, p_cancel)) { + SolveSupport(support, isSingular, budgetExceeded, + std::max(p_stopAfter - static_cast(ret.size()), 0), p_maxregret, + p_maxRectangles, p_cancel)) { p_onEquilibrium(solution); ret.push_back(solution); } if (isSingular) { p_onEvent(EnumPolySingularSupportEvent{support}); } + if (budgetExceeded) { + p_onEvent(EnumPolyBudgetExceededSupportEvent{support}); + } if (p_stopAfter > 0 && static_cast(ret.size()) >= p_stopAfter) { break; } diff --git a/src/solvers/enumpoly/enumpoly.h b/src/solvers/enumpoly/enumpoly.h index 0f3aa75a2d..f148a1e3e5 100644 --- a/src/solvers/enumpoly/enumpoly.h +++ b/src/solvers/enumpoly/enumpoly.h @@ -39,17 +39,31 @@ template struct EnumPolySingularSupportEvent { const Support &support; }; +/// Reported when the search for roots on a support exhausted its rectangle +/// budget before completing. Distinct from EnumPolySingularSupportEvent: +/// this can happen even when the equilibria being sought are regular. +template struct EnumPolyBudgetExceededSupportEvent { + const Support &support; +}; + template using EnumPolyEvent = - std::variant, EnumPolySingularSupportEvent>; + std::variant, EnumPolySingularSupportEvent, + EnumPolyBudgetExceededSupportEvent>; template using EnumPolyEventCallbackType = std::function &)>; template void NullEnumPolyEventCallback(const EnumPolyEvent &) {} +/// Default budget on the number of rectangles examined per support before +/// giving up. Ordinary supports need only a handful, leaving ample headroom +/// for difficult-but-tractable ones while bounding the worst case. +constexpr size_t kDefaultEnumPolyMaxRectangles = 20'000; + std::list> EnumPolyStrategySolve(const Game &p_game, int p_stopAfter, double p_maxregret, + size_t p_maxRectangles = kDefaultEnumPolyMaxRectangles, StrategyCallbackType p_onEquilibrium = NullStrategyCallback, EnumPolyEventCallbackType p_onEvent = NullEnumPolyEventCallback, @@ -57,6 +71,7 @@ EnumPolyStrategySolve(const Game &p_game, int p_stopAfter, double p_maxregret, std::list> EnumPolyBehaviorSolve(const Game &, int p_stopAfter, double p_maxregret, + size_t p_maxRectangles = kDefaultEnumPolyMaxRectangles, BehaviorCallbackType p_onEquilibrium = NullBehaviorCallback, EnumPolyEventCallbackType p_onEvent = NullEnumPolyEventCallback, diff --git a/src/solvers/enumpoly/lufactor.cc b/src/solvers/enumpoly/lufactor.cc new file mode 100644 index 0000000000..551348f3a4 --- /dev/null +++ b/src/solvers/enumpoly/lufactor.cc @@ -0,0 +1,111 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (http://www.gambit-project.org) +// +// FILE: src/solvers/enumpoly/lufactor.cc +// LU factorization of a square matrix of doubles +// +// 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 "lufactor.h" + +namespace Gambit { + +namespace { +constexpr double kSingularityTolerance = 1.0e-9; +} // end anonymous namespace + +LUFactorization::LUFactorization(const Matrix &m) + : m_lu(m), m_pivot(m.MinRow(), m.MaxRow()) +{ + if (!m.IsSquare()) { + throw SingularMatrixException(); + } + const int lo = m_lu.MinRow(), hi = m_lu.MaxRow(); + for (int i = lo; i <= hi; i++) { + m_pivot[i] = i; + } + for (int k = lo; k <= hi; k++) { + int pivotRow = k; + double maxVal = std::abs(m_lu(k, k)); + for (int i = k + 1; i <= hi; i++) { + if (std::abs(m_lu(i, k)) > maxVal) { + maxVal = std::abs(m_lu(i, k)); + pivotRow = i; + } + } + if (maxVal <= kSingularityTolerance) { + throw SingularMatrixException(); + } + if (pivotRow != k) { + m_lu.SwitchRows(k, pivotRow); + std::swap(m_pivot[k], m_pivot[pivotRow]); + m_determinant = -m_determinant; + } + m_determinant *= m_lu(k, k); + for (int i = k + 1; i <= hi; i++) { + const double factor = m_lu(i, k) / m_lu(k, k); + m_lu(i, k) = factor; + for (int j = k + 1; j <= hi; j++) { + m_lu(i, j) -= factor * m_lu(k, j); + } + } + } +} + +Vector LUFactorization::Solve(const Vector &b) const +{ + if (!m_lu.ConformsToColumn(b)) { + throw DimensionException(); + } + const int lo = m_lu.MinRow(), hi = m_lu.MaxRow(); + Vector x(lo, hi); + for (int i = lo; i <= hi; i++) { + x[i] = b[m_pivot[i]]; + } + // Forward substitution: L y = Pb (L has unit diagonal, stored below it in m_lu) + for (int i = lo + 1; i <= hi; i++) { + for (int j = lo; j < i; j++) { + x[i] -= m_lu(i, j) * x[j]; + } + } + // Back substitution: U x = y + for (int i = hi; i >= lo; i--) { + for (int j = i + 1; j <= hi; j++) { + x[i] -= m_lu(i, j) * x[j]; + } + x[i] /= m_lu(i, i); + } + return x; +} + +Matrix LUFactorization::Inverse() const +{ + const int lo = m_lu.MinRow(), hi = m_lu.MaxRow(); + Matrix inv(lo, hi, lo, hi); + for (int j = lo; j <= hi; j++) { + Vector e(lo, hi); + e = 0.0; + e[j] = 1.0; + const Vector col = Solve(e); + for (int i = lo; i <= hi; i++) { + inv(i, j) = col[i]; + } + } + return inv; +} + +} // end namespace Gambit diff --git a/src/solvers/enumpoly/lufactor.h b/src/solvers/enumpoly/lufactor.h new file mode 100644 index 0000000000..a785254ba4 --- /dev/null +++ b/src/solvers/enumpoly/lufactor.h @@ -0,0 +1,52 @@ +// +// This file is part of Gambit +// Copyright (c) 1994-2026, The Gambit Project (http://www.gambit-project.org) +// +// FILE: src/solvers/enumpoly/lufactor.h +// LU factorization of a square matrix of doubles +// +// 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_ENUMPOLY_LUFACTOR_H +#define GAMBIT_SOLVERS_ENUMPOLY_LUFACTOR_H + +#include "gambit.h" + +namespace Gambit { + +/// LU factorization (with partial pivoting) of a square matrix, computed once +/// and reused for Solve()/Determinant()/Inverse() without repeating the +/// O(n^3) elimination. Throws SingularMatrixException if the matrix is +/// singular, or any pivot is too small to trust the factorization. +class LUFactorization { +public: + explicit LUFactorization(const Matrix &m); + + double Determinant() const { return m_determinant; } + + /// @throws DimensionException if `b` does not conform to the factored matrix + Vector Solve(const Vector &b) const; + Matrix Inverse() const; + +private: + Matrix m_lu; + Array m_pivot; + double m_determinant{1.0}; +}; + +} // end namespace Gambit + +#endif // GAMBIT_SOLVERS_ENUMPOLY_LUFACTOR_H diff --git a/src/solvers/enumpoly/nfgpoly.cc b/src/solvers/enumpoly/nfgpoly.cc index 361cee4280..d257e23edd 100644 --- a/src/solvers/enumpoly/nfgpoly.cc +++ b/src/solvers/enumpoly/nfgpoly.cc @@ -103,7 +103,8 @@ namespace Gambit::Nash { std::list> EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_singular, - int p_stopAfter, const CancelToken &p_cancel) + bool &p_budgetExceeded, int p_stopAfter, size_t p_maxRectangles, + const CancelToken &p_cancel) { auto space = std::make_shared(support.MixedProfileLength() - support.GetGame()->NumPlayers()); @@ -120,7 +121,7 @@ EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_sin try { roots = solver.FindRoots({bottoms, tops}, (p_stopAfter > 0) ? p_stopAfter : std::numeric_limits::max(), - p_cancel); + p_maxRectangles, p_budgetExceeded, p_cancel); } catch (const SingularMatrixException &) { is_singular = true; @@ -142,7 +143,7 @@ EnumPolyStrategySupportSolve(const StrategySupportProfile &support, bool &is_sin std::list> EnumPolyStrategySolve(const Game &p_game, int p_stopAfter, double p_maxregret, - StrategyCallbackType p_onEquilibrium, + size_t p_maxRectangles, StrategyCallbackType p_onEquilibrium, EnumPolyEventCallbackType p_onEvent, const CancelToken &p_cancel) { @@ -163,9 +164,10 @@ EnumPolyStrategySolve(const Game &p_game, int p_stopAfter, double p_maxregret, p_cancel.Check(); p_onEvent(EnumPolyCandidateSupportEvent{support}); bool is_singular; + bool budget_exceeded; for (auto solution : EnumPolyStrategySupportSolve( - support, is_singular, std::max(p_stopAfter - static_cast(ret.size()), 0), - p_cancel)) { + support, is_singular, budget_exceeded, + std::max(p_stopAfter - static_cast(ret.size()), 0), p_maxRectangles, p_cancel)) { const MixedStrategyProfile fullProfile = solution.ToFullSupport(); if (fullProfile.GetMaxRegret() < p_maxregret) { p_onEquilibrium(fullProfile); @@ -176,6 +178,9 @@ EnumPolyStrategySolve(const Game &p_game, int p_stopAfter, double p_maxregret, if (is_singular) { p_onEvent(EnumPolySingularSupportEvent{support}); } + if (budget_exceeded) { + p_onEvent(EnumPolyBudgetExceededSupportEvent{support}); + } if (p_stopAfter > 0 && static_cast(ret.size()) >= p_stopAfter) { break; } diff --git a/src/solvers/enumpoly/polysolver.cc b/src/solvers/enumpoly/polysolver.cc index bffd66929a..69983f44a0 100644 --- a/src/solvers/enumpoly/polysolver.cc +++ b/src/solvers/enumpoly/polysolver.cc @@ -22,6 +22,7 @@ #include "polysolver.h" #include "indexproduct.h" +#include "lufactor.h" namespace Gambit { @@ -149,7 +150,7 @@ double PolynomialSystemSolver::MaxDistanceFromPointToVertexAfterTransformation( bool PolynomialSystemSolver::HasNoOtherRootsIn(const Rectangle &r, const Vector &p) const { - auto M = m_derivatives.SquareDerivativeMatrix(p).Inverse(); + auto M = LUFactorization(m_derivatives.SquareDerivativeMatrix(p)).Inverse(); auto system2 = m_normalizedSystem.Translate(p).TransformCoords(M); double radius = MaxDistanceFromPointToVertexAfterTransformation(r, p, M); return (std::accumulate(system2.begin(), system2.end(), 0.0, @@ -162,51 +163,52 @@ Vector PolynomialSystemSolver::NewtonStep(const Vector &point) c { const Vector evals = m_derivatives.ValuesOfRootPolys(point, NumEquations()); const Matrix deriv = m_derivatives.DerivativeMatrix(point, NumEquations()); - auto transpose = deriv.Transpose(); - auto sqmat = Matrix(deriv * transpose); - if (std::abs(sqmat.Determinant()) <= 1.0e-9) { - throw SingularMatrixException(); - } - return point - (transpose * sqmat.Inverse()) * evals; + return point - LUFactorization(deriv).Solve(evals); } Vector PolynomialSystemSolver::ImprovingNewtonStep(const Vector &point) const { const Vector evals = m_derivatives.ValuesOfRootPolys(point, NumEquations()); const Matrix deriv = m_derivatives.DerivativeMatrix(point, NumEquations()); - auto transpose = deriv.Transpose(); - auto sqmat = Matrix(deriv * transpose); - if (std::abs(sqmat.Determinant()) <= 1.0e-9) { - throw SingularMatrixException(); - } - auto step = -(transpose * sqmat.Inverse()) * evals; + const Vector delta = LUFactorization(deriv).Solve(evals); - while (m_derivatives.ValuesOfRootPolys(point + step, NumEquations()).NormSquared() > + double scale = 1.0; + while (m_derivatives.ValuesOfRootPolys(point - delta * scale, NumEquations()).NormSquared() > evals.NormSquared()) { - step /= 2; + scale /= 2; } - return point + step; + return point - delta * scale; } std::list> PolynomialSystemSolver::FindRoots(const Rectangle &r, const int max_roots, + const size_t max_rectangles, + bool &p_budgetExceeded, const CancelToken &p_cancel) { + p_budgetExceeded = false; std::list> roots; if (NumEquations() == 0) { roots.emplace_back(); } else { - FindRoots(roots, r, max_roots, p_cancel); + size_t budget = max_rectangles; + FindRoots(roots, r, max_roots, budget, p_budgetExceeded, p_cancel); } return roots; } void PolynomialSystemSolver::FindRoots(std::list> &rootlist, const Rectangle &r, const size_t max_roots, + size_t &p_budget, bool &p_budgetExceeded, const CancelToken &p_cancel) const { p_cancel.Check(); + if (p_budget == 0) { + p_budgetExceeded = true; + return; + } + --p_budget; if (SystemHasNoRootsIn(r)) { return; } @@ -230,8 +232,8 @@ void PolynomialSystemSolver::FindRoots(std::list> &rootlist, return; } for (const auto &cell : r.Orthants()) { - FindRoots(rootlist, cell, max_roots, p_cancel); - if (rootlist.size() >= max_roots) { + FindRoots(rootlist, cell, max_roots, p_budget, p_budgetExceeded, p_cancel); + if (rootlist.size() >= max_roots || p_budgetExceeded) { return; } } diff --git a/src/solvers/enumpoly/polysolver.h b/src/solvers/enumpoly/polysolver.h index d074deaf14..7d1cda68b2 100644 --- a/src/solvers/enumpoly/polysolver.h +++ b/src/solvers/enumpoly/polysolver.h @@ -43,9 +43,11 @@ namespace Gambit { /// rectangle are not resolved by these queries, the rectangle is /// subdivided into 2^d subrectangles, and the process is repeated on /// each. This continues until it has been shown that all roots have been -/// found, or a predetermined search depth is reached. The bound on depth -/// is necessary because the procedure will not terminate if there are -/// singular roots. +/// found, a predetermined search depth is reached, or a budget on the +/// number of rectangles examined is exhausted. The depth bound alone is +/// not a practical guarantee of termination, since the number of +/// rectangles examined can still grow very quickly with depth; the +/// rectangle budget caps the total work directly. /// The list of polynomials must /// be at least as long as the dimension GetDimension() of the space of the @@ -71,8 +73,12 @@ class PolynomialSystemSolver { bool HasNoOtherRootsIn(const Rectangle &, const Vector &) const; + // Recursive worker for FindRoots(). `p_budget` counts down the number of + // rectangles still allowed to be examined; when it reaches zero, sets + // `p_budgetExceeded` and the recursion unwinds without exploring further. void FindRoots(std::list> &rootlist, const Rectangle &r, - const size_t max_roots, const CancelToken &p_cancel) const; + const size_t max_roots, size_t &p_budget, bool &p_budgetExceeded, + const CancelToken &p_cancel) const; public: explicit PolynomialSystemSolver(const PolynomialSystem &p_system) @@ -96,8 +102,11 @@ class PolynomialSystemSolver { /// is smaller than at the original point Vector ImprovingNewtonStep(const Vector &) const; - // Find up to `max_roots` roots inside rectangle `r` + // Find up to `max_roots` roots inside rectangle `r`, examining at most + // `max_rectangles` subdivisions. Sets `p_budgetExceeded` if that budget + // runs out first, in which case the roots returned may be incomplete. std::list> FindRoots(const Rectangle &r, int max_roots, + size_t max_rectangles, bool &p_budgetExceeded, const CancelToken &p_cancel = CancelToken()); }; diff --git a/src/tools/enumpoly/enumpoly.cc b/src/tools/enumpoly/enumpoly.cc index 0387ee4d08..e43e9d7923 100644 --- a/src/tools/enumpoly/enumpoly.cc +++ b/src/tools/enumpoly/enumpoly.cc @@ -20,6 +20,7 @@ // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. // +#include #include #include #include @@ -52,9 +53,13 @@ void PrintHelp(char *progname) std::cerr << " -h, --help print this help message\n"; std::cerr << " -S use strategic game\n"; std::cerr << " -m MAXREGRET maximum regret acceptable as a proportion of range of\n"; - std::cerr << " payoffs in the game\n"; + std::cerr << " payoffs in the game (default is 1e-8)\n"; std::cerr << " -e EQA terminate after finding EQA equilibria\n"; std::cerr << " (default is to search in all supports)\n"; + std::cerr << " -r RECTANGLES maximum number of rectangles to examine when searching\n"; + std::cerr << " for roots on a single support, before giving up on that\n"; + std::cerr << " support (default " << Nash::kDefaultEnumPolyMaxRectangles + << ")\n"; std::cerr << " -q quiet mode (suppresses banner)\n"; std::cerr << " -V, --verbose verbose mode (shows supports investigated)\n"; std::cerr << " (default is only to show equilibria)\n"; @@ -128,6 +133,9 @@ template void PrintEnumPolyEvent(const EnumPolyEvent &p else if constexpr (std::is_same_v>) { PrintSupport(std::cout, "singular", event.support); } + else if constexpr (std::is_same_v>) { + PrintSupport(std::cout, "budget-exceeded", event.support); + } }, p_event); } @@ -138,8 +146,9 @@ int main(int argc, char *argv[]) bool quiet = false; bool useStrategic = false; - double maxregret = 1.0e-4; + double maxregret = 1.0e-8; int stopAfter = 0; + size_t maxRectangles = Nash::kDefaultEnumPolyMaxRectangles; int long_opt_index = 0; option long_options[] = {{"help", 0, nullptr, 'h'}, @@ -147,7 +156,7 @@ int main(int argc, char *argv[]) {"verbose", 0, nullptr, 'V'}, {nullptr, 0, nullptr, 0}}; int c; - while ((c = getopt_long(argc, argv, "d:hSm:e:qvV", long_options, &long_opt_index)) != -1) { + while ((c = getopt_long(argc, argv, "d:hSm:e:r:qvV", long_options, &long_opt_index)) != -1) { switch (c) { case 'v': PrintBanner(std::cerr); @@ -167,6 +176,9 @@ int main(int argc, char *argv[]) case 'e': stopAfter = atoi(optarg); break; + case 'r': + maxRectangles = std::strtoull(optarg, nullptr, 10); + break; case 'q': quiet = true; break; @@ -212,13 +224,13 @@ int main(int argc, char *argv[]) if (!game->IsTree() || useStrategic) { EnumPolyStrategySolve( - game, stopAfter, maxregret, + game, stopAfter, maxregret, maxRectangles, [](const MixedStrategyProfile &eqm) { PrintProfile(std::cout, "NE", eqm); }, [](const EnumPolyEvent &event) { PrintEnumPolyEvent(event); }); } else { EnumPolyBehaviorSolve( - game, stopAfter, maxregret, + game, stopAfter, maxregret, maxRectangles, [](const MixedBehaviorProfile &eqm) { PrintProfile(std::cout, "NE", eqm); }, [](const EnumPolyEvent &event) { PrintEnumPolyEvent(event); }); } diff --git a/tests/games.py b/tests/games.py index 1a2dda3ed6..dbe5602f96 100644 --- a/tests/games.py +++ b/tests/games.py @@ -25,8 +25,20 @@ UNICODE_LABELS = ["é", "naïve", "日本語", "😀"] VALID_LABELS = ["x", "a b", "a b", "a b c", *UNICODE_LABELS] INVALID_LABELS = [ - " x", "x ", " ", "a b", "a\tb", "a\nb", "a\x01b", "a\x7fb", "a…b", - " x", "x ", " ", "a  b", "a
b", + " x", + "x ", + " ", + "a b", + "a\tb", + "a\nb", + "a\x01b", + "a\x7fb", + "a…b", + " x", + "x ", + " ", + "a  b", + "a
b", ] @@ -43,7 +55,7 @@ def read_from_file(fn: str) -> gbt.Game: raise ValueError(f"Unknown file extension in {fn}") -def create_efg_corresponding_to_bimatrix_game( +def create_efg_corresponding_to_bimatrix_game_arrays( A: np.ndarray, B: np.ndarray, title: str ) -> gbt.Game: """ @@ -61,11 +73,22 @@ def create_efg_corresponding_to_bimatrix_game( for i, j in itertools.product(range(m), range(n)): g.set_outcome( g.root.children[str(i)].children[str(j)], - g.add_outcome(f"({i},{j})", [A[i, j], B[i, j]]) + g.add_outcome(f"({i},{j})", [A[i, j], B[i, j]]), ) return g +def create_efg_corresponding_to_bimatrix_game(g: gbt.Game) -> gbt.Game: + """ + There is no direct pygambit method to create an EFG from a stategic-form game. + Here we create an EFG corresponding to a bimatrix pygambit.Game. + Player 1 moves first. + """ + assert len(g.players) == 2 + A, B = g.to_arrays() + return create_efg_corresponding_to_bimatrix_game_arrays(A, B, g.title) + + ################################################################################################ # Extensive-form games (efg) @@ -85,7 +108,7 @@ def create_2x2_zero_sum_efg(variant: None | str = None) -> gbt.Game: A = np.eye(2) B = -A - g = create_efg_corresponding_to_bimatrix_game(A, B, title) + g = create_efg_corresponding_to_bimatrix_game_arrays(A, B, title) if variant == "missing term outcome": g.delete_outcome(g.root.children["0"].children["1"].outcome) @@ -406,24 +429,16 @@ def create_one_shot_trust_efg(unique_NE_variant: bool = False) -> gbt.Game: ) g.append_move(g.root, "Buyer", ["Trust", "Not trust"]) g.append_move(g.root.children["Trust"], "Seller", ["Honor", "Abuse"]) - g.set_outcome( - g.root.children["Trust"].children["Honor"], - g.add_outcome("Trustworthy", [1, 1]) - ) + g.set_outcome(g.root.children["Trust"].children["Honor"], g.add_outcome("Trustworthy", [1, 1])) if unique_NE_variant: g.set_outcome( - g.root.children["Trust"].children["Abuse"], - g.add_outcome("Untrustworthy", ["1/2", 2]) + g.root.children["Trust"].children["Abuse"], g.add_outcome("Untrustworthy", ["1/2", 2]) ) else: g.set_outcome( - g.root.children["Trust"].children["Abuse"], - g.add_outcome("Untrustworthy", [-1, 2]) - ) - g.set_outcome( - g.root.children["Not trust"], - g.add_outcome("Opt-out", [0, 0]) + g.root.children["Trust"].children["Abuse"], g.add_outcome("Untrustworthy", [-1, 2]) ) + g.set_outcome(g.root.children["Not trust"], g.add_outcome("Opt-out", [0, 0])) return g @@ -431,7 +446,7 @@ def create_EFG_for_nxn_bimatrix_coordination_game(n: int) -> gbt.Game: A = np.eye(n, dtype=int) B = A title = f"{n}x{n} coordination game, {2**n - 1} equilibria" - return create_efg_corresponding_to_bimatrix_game(A, B, title) + return create_efg_corresponding_to_bimatrix_game_arrays(A, B, title) def create_EFG_for_6x6_bimatrix_with_long_LH_paths_and_unique_eq() -> gbt.Game: @@ -456,7 +471,7 @@ def create_EFG_for_6x6_bimatrix_with_long_LH_paths_and_unique_eq() -> gbt.Game: A = np.array(A) B = np.array(B) title = "6x6 Long Lemke-Howson Paths, unique eq" - return create_efg_corresponding_to_bimatrix_game(A, B, title) + return create_efg_corresponding_to_bimatrix_game_arrays(A, B, title) class EfgFamilyForReducedStrategicFormTests(ABC): @@ -686,9 +701,11 @@ def gbt_game(self): for n in g.nodes: if not n.is_terminal and not n.children["L"].is_terminal: left = n.children["L"] - g.make_infoset(list(left.infoset.members) + [n.children["R"]], - left.infoset.player.label, - left.infoset.label or None) + g.make_infoset( + list(left.infoset.members) + [n.children["R"]], + left.infoset.player.label, + left.infoset.label or None, + ) return g def reduced_strategic_form(self): diff --git a/tests/test_nash.py b/tests/test_nash.py index e715dca138..607bdce1e3 100644 --- a/tests/test_nash.py +++ b/tests/test_nash.py @@ -421,7 +421,7 @@ class QREquilibriumTestCase: d("165/347", 0, 0, "42/347", "70/347", "70/347"), d(0, 0, "162/2297", "28/2297", "1540/2297", "567/2297"), ], - [d(0, "1/1", 0, 0, 0, 0), d(0, 0, 0, 0, 0, "1/1")], + [d(0, 1, 0, 0, 0, 0), d(0, 0, 0, 0, 0, 1)], [d(0, 0, "3/5", "2/5", 0, 0), d("3/13", "10/13", 0, 0, 0, 0)], [d(0, 0, "7/12", "5/12", 0, 0), d(0, "25/32", "7/32", 0, 0, 0)], [d(0, 0, "4/7", "3/7", 0, 0), d(0, 0, "4/7", "3/7", 0, 0)], @@ -490,7 +490,7 @@ class QREquilibriumTestCase: d(0, "77/122", "50/183", 0, 0, "35/366"), d("160/629", "315/629", 0, 0, 0, "154/629"), ], - [d(0, 0, 0, 0, "1/1", 0), d("1/1", 0, 0, 0, 0, 0)], + [d(0, 0, 0, 0, 1, 0), d(1, 0, 0, 0, 0, 0)], [d(0, "11/15", 0, 0, "4/15", 0), d("4/15", 0, 0, 0, 0, "11/15")], [d(0, 0, 0, 0, "2/7", "5/7"), d("2/7", "5/7", 0, 0, 0, 0)], [d(0, 0, 0, 0, "1/4", "3/4"), d(0, "2/3", "1/3", 0, 0, 0)], @@ -756,7 +756,7 @@ class QREquilibriumTestCase: d("165/347", 0, 0, "42/347", "70/347", "70/347"), d(0, 0, "162/2297", "28/2297", "1540/2297", "567/2297"), ], - [d(0, "1/1", 0, 0, 0, 0), d(0, 0, 0, 0, 0, "1/1")], + [d(0, 1, 0, 0, 0, 0), d(0, 0, 0, 0, 0, 1)], [d(0, 0, "3/5", "2/5", 0, 0), d("3/13", "10/13", 0, 0, 0, 0)], [d(0, 0, "7/12", "5/12", 0, 0), d(0, "25/32", "7/32", 0, 0, 0)], [d(0, 0, "4/7", "3/7", 0, 0), d(0, 0, "4/7", "3/7", 0, 0)], @@ -825,7 +825,7 @@ class QREquilibriumTestCase: d(0, "77/122", "50/183", 0, 0, "35/366"), d("160/629", "315/629", 0, 0, 0, "154/629"), ], - [d(0, 0, 0, 0, "1/1", 0), d("1/1", 0, 0, 0, 0, 0)], + [d(0, 0, 0, 0, 1, 0), d(1, 0, 0, 0, 0, 0)], [d(0, "11/15", 0, 0, "4/15", 0), d("4/15", 0, 0, 0, 0, "11/15")], [d(0, 0, 0, 0, "2/7", "5/7"), d("2/7", "5/7", 0, 0, 0, 0)], [d(0, 0, 0, 0, "1/4", "3/4"), d(0, "2/3", "1/3", 0, 0, 0)], @@ -917,7 +917,34 @@ class QREquilibriumTestCase: marks=pytest.mark.nash_enumpoly_strategy, id="test_enumpoly_strategy_3", ), + # 6x6 bimatrix game with 75 equilibria + pytest.param( + EquilibriumTestCase( + factory=functools.partial( + gbt.catalog.load, "journals/dcg/vonstengel1999/6x6_game_with_75_eq" + ), + solver=functools.partial(gbt.nash.enumpoly_solve, stop_after=10), + expected=[ + [d(0, 1, 0, 0, 0, 0), d(0, 0, 0, 0, 0, 1)], + [d(0, 0, 0, 0, 1, 0), d(1, 0, 0, 0, 0, 0)], + [d("20/23", "3/23", 0, 0, 0, 0), d("4/19", "15/19", 0, 0, 0, 0)], + [d("45/52", "7/52", 0, 0, 0, 0), d(0, "14/17", "3/17", 0, 0, 0)], + [d("25/29", "4/29", 0, 0, 0, 0), d(0, 0, "5/9", "4/9", 0, 0)], + [d("55/64", "9/64", 0, 0, 0, 0), d(0, 0, 0, "3/58", "55/58", 0)], + [d("6/7", "1/7", 0, 0, 0, 0), d(0, 0, 0, 0, "6/7", "1/7")], + [d("28/29", 0, 0, "1/29", 0, 0), d("7/32", "25/32", 0, 0, 0, 0)], + [d("24/25", 0, 0, "1/25", 0, 0), d(0, "25/31", "6/31", 0, 0, 0)], + [d("150/157", 0, 0, "7/157", 0, 0), d(0, 0, "9/16", "7/16", 0, 0)], + ], + regret_tol=TOL_LARGE, + prob_tol=TOL_LARGE, + ), + marks=pytest.mark.nash_enumpoly_strategy, + id="test_enumpoly_strategy_4", + ), + ################################# # 3-player games + ################################# # 2x2x2 strategic form game based on local max cut -- 2 pure and 1 mixed pytest.param( EquilibriumTestCase( @@ -934,7 +961,7 @@ class QREquilibriumTestCase: regret_tol=TOL, ), marks=pytest.mark.nash_enumpoly_strategy, - id="test_enumpoly_strategy_4", + id="test_enumpoly_strategy_5", ), # A three-player game with a unique Nash equilibrium in irrational mixed strategies # (nau2004 sec4 catalog game) @@ -953,7 +980,7 @@ class QREquilibriumTestCase: regret_tol=TOL, ), marks=pytest.mark.nash_enumpoly_strategy, - id="test_enumpoly_strategy_5", + id="test_enumpoly_strategy_6", ), # A three-player 2x2x2 game with 3 pure, 2 incompletely mixed, and a # continuum of completely mixed Nash equilibria (nau2004 sec5 catalog game) @@ -970,7 +997,7 @@ class QREquilibriumTestCase: regret_tol=TOL, ), marks=pytest.mark.nash_enumpoly_strategy, - id="test_enumpoly_strategy_6", + id="test_enumpoly_strategy_7", ), ] @@ -2895,7 +2922,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s expected=[ [[d("1/2", "1/2")], [d("2/5", "3/5"), d(0, 1), d(1, 0)]], [[d("1/4", "3/4")], [d(0, 1), d(1, 0), d("1/3", "2/3")]], - [[d(1, 0)], [d(1, 0), d(1, 0), d(1, 0)]] + [[d(1, 0)], [d(1, 0), d(1, 0), d(1, 0)]], ], regret_tol=TOL, prob_tol=TOL, @@ -2903,6 +2930,31 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s marks=pytest.mark.nash_enumpoly_behavior, id="test_enumpoly_behavior_10", ), + # 6x6 bimatrix game with 75 equilibria -- runs slow so only do first 2 equilibria + pytest.param( + EquilibriumTestCase( + factory=functools.partial( + games.create_efg_corresponding_to_bimatrix_game, + gbt.catalog.load("journals/dcg/vonstengel1999/6x6_game_with_75_eq"), + ), + solver=functools.partial(gbt.nash.enumpoly_solve, stop_after=2), + expected=[ + [ + [d("132/245", "66/245", "6/49", "12/245", "4/245", "1/245")], + [d("4/245", "1/245", "6/49", "12/245", "132/245", "66/245")], + ], + [ + [d("5775/9871", "2640/9871", "1050/9871", "336/9871", "70/9871", 0)], + [d("70/9871", 0, "1050/9871", "336/9871", "5775/9871", "2640/9871")], + ], + ], + regret_tol=TOL_LARGE, + prob_tol=TOL_LARGE, + ), + marks=pytest.mark.nash_enumpoly_behavior, + id="test_enumpoly_behavior_11", + ), + ######################################################################## # 3-player games pytest.param( EquilibriumTestCase( @@ -2915,7 +2967,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_11", + id="test_enumpoly_behavior_12", ), pytest.param( EquilibriumTestCase( @@ -2928,7 +2980,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_12", + id="test_enumpoly_behavior_13", ), # 3-player perfect info game to test behavior two off equilibrium path pytest.param( @@ -2945,7 +2997,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_13", + id="test_enumpoly_behavior_14", ), pytest.param( EquilibriumTestCase( @@ -2963,7 +3015,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_14", + id="test_enumpoly_behavior_15", ), pytest.param( EquilibriumTestCase( @@ -2983,7 +3035,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_15", + id="test_enumpoly_behavior_16", ), pytest.param( EquilibriumTestCase( @@ -2999,7 +3051,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_16", + id="test_enumpoly_behavior_17", ), # 3-player game equivalent to a simultaneous-move 2x2x2 game (nine equilibria total, # including two totally-mixed ones); this previously required an unordered check due @@ -3023,7 +3075,7 @@ def test_nash_strategy_solver_w_start(test_case: EquilibriumTestCaseWithStart, s prob_tol=TOL, ), marks=pytest.mark.nash_enumpoly_behavior, - id="test_enumpoly_behavior_17", + id="test_enumpoly_behavior_18", ), ]