Addresses two sources of inconsistency in the output of enumpoly:

ChangeLog

@@ -8,6 +8,9 @@
   strategy with an empty label; it is now labeled `"1"`, as `Game.new_table` already labels
   every other player's initial strategies.  (#1047)
 - Numerical data in action graph games is now read with exact precision, to match other game formats.  (#879)
+- `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)
 
 ### Removed
 - Assigning to `Strategy.label` has been removed; use `Game.relabel_strategies`, which enforces

src/games/behavspt.cc

@@ -223,39 +223,39 @@ size_t BehaviorSupportProfile::Sequences::size() const
 
 BehaviorSupportProfile::Sequences::iterator BehaviorSupportProfile::Sequences::begin() const
 {
-  return {m_support->GetSequenceMap(), false};
+  return {m_support->GetSequenceMap(), m_support->GetPlayers(), false};
 }
 BehaviorSupportProfile::Sequences::iterator BehaviorSupportProfile::Sequences::end() const
 {
-  return {m_support->GetSequenceMap(), true};
+  return {m_support->GetSequenceMap(), m_support->GetPlayers(), true};
 }
 
 BehaviorSupportProfile::Sequences::iterator::iterator(
-    const std::shared_ptr<SequenceMap> p_sequences, bool p_end)
-  : m_sequences(p_sequences)
+    const std::shared_ptr<SequenceMap> p_sequences, const GameRep::Players &p_players, bool p_end)
+  : m_sequences(p_sequences), m_players(p_players)
 {
   if (p_end) {
-    m_currentPlayer = m_sequences->cend();
+    m_currentPlayer = m_players.end();
   }
   else {
-    m_currentPlayer = m_sequences->cbegin();
-    m_currentSequence = m_currentPlayer->second.cbegin();
+    m_currentPlayer = m_players.begin();
+    m_currentSequence = m_sequences->at(*m_currentPlayer).cbegin();
   }
 }
 
 BehaviorSupportProfile::Sequences::iterator &
 BehaviorSupportProfile::Sequences::iterator::operator++()
 {
-  if (m_currentPlayer == m_sequences->cend()) {
+  if (m_currentPlayer == m_players.end()) {
     return *this;
   }
   m_currentSequence++;
-  if (m_currentSequence != m_currentPlayer->second.cend()) {
+  if (m_currentSequence != m_sequences->at(*m_currentPlayer).cend()) {
     return *this;
   }
-  m_currentPlayer++;
-  if (m_currentPlayer != m_sequences->cend()) {
-    m_currentSequence = m_currentPlayer->second.cbegin();
+  ++m_currentPlayer;
+  if (m_currentPlayer != m_players.end()) {
+    m_currentSequence = m_sequences->at(*m_currentPlayer).cbegin();
   }
   return *this;
 }
@@ -265,7 +265,7 @@ bool BehaviorSupportProfile::Sequences::iterator::operator==(const iterator &it)
   if (m_sequences != it.m_sequences || m_currentPlayer != it.m_currentPlayer) {
     return false;
   }
-  if (m_currentPlayer == m_sequences->end()) {
+  if (m_currentPlayer == m_players.end()) {
     return true;
   }
   return (m_currentSequence == it.m_currentSequence);

src/games/behavspt.h

@@ -156,12 +156,20 @@ class BehaviorSupportProfile {
 
   public:
     class iterator {
+      // Sequences are visited grouped by player, in the game's canonical player order
+      // (GameRep::Players), rather than in SequenceMap's own key order.  SequenceMap is
+      // keyed by GamePlayer, whose ordering is based on object addresses; iterating the
+      // map directly would therefore make the visitation order (and hence, downstream,
+      // any polynomial variable numbering derived from it) depend on heap layout, which
+      // can differ from run to run.
       const std::shared_ptr<SequenceMap> m_sequences;
-      SequenceMap::const_iterator m_currentPlayer;
+      GameRep::Players m_players;
+      GameRep::Players::iterator m_currentPlayer;
       std::vector<GameSequence>::const_iterator m_currentSequence;
 
     public:
-      iterator(const std::shared_ptr<SequenceMap> p_sequences, bool p_end);
+      iterator(const std::shared_ptr<SequenceMap> p_sequences, const GameRep::Players &p_players,
+               bool p_end);
 
       GameSequence operator*() const { return *m_currentSequence; }
       GameSequence operator->() const { return *m_currentSequence; }

src/solvers/enumpoly/polysolver.cc

@@ -72,6 +72,7 @@ bool PolynomialSystemSolver::NewtonRootInRectangle(const Rectangle<double> &r,
       newpoint = NewtonStep(point);
     }
     catch (SingularMatrixException &) {
+      bool perturbed = false;
       for (int i = 1; i <= point.size(); i++) {
         Vector<double> perturbed_point(point);
         if (r.Side(i).UpperBound() > point[i]) {
@@ -82,11 +83,21 @@ bool PolynomialSystemSolver::NewtonRootInRectangle(const Rectangle<double> &r,
         }
         try {
           newpoint = point + (NewtonStep(perturbed_point) - perturbed_point);
+          perturbed = true;
           break;
         }
         catch (SingularMatrixException &) {
         }
       }
+      if (!perturbed) {
+        // The Jacobian is singular at this point, and remains singular under every
+        // single-coordinate perturbation tried.  This is a local degeneracy of the
+        // search at this point/rectangle, not evidence that the system of equations
+        // for the support as a whole is degenerate -- so we treat it the same as any
+        // other failure to confirm a root here, leaving it to the caller to subdivide
+        // the rectangle further.
+        return false;
+      }
     }
 
     if (!r.SameCenterDoubleLengths().Contains(newpoint)) {