Change logit numerical continuation to terminate on regret-based crit…

…erion.

This changes the numerical continuation in logit to terminate when it reaches a
profile which has less than a specified level of regret (epsilon), which is
interpreted relative to the scale of payoffs of the game.

Closes #365.

doc/tools.logit.rst

@@ -33,6 +33,13 @@ if an information set is not reached due to being the successor of chance
 moves with zero probability.  In such games, the implementation treats
 the beliefs at such information sets as being uniform across all member nodes.
 
+.. versionchanged:: 16.2.0
+
+   The criterion for accepting whether a point is sufficiently close to a
+   Nash equilibrium to terminate the path-following is specified
+   in terms of the maximum regret.  This regret is interpreted as a fraction
+   of the difference between the maximum and minimum payoffs in the game.
+
 .. program:: gambit-logit
 
 .. cmdoption:: -d
@@ -61,14 +68,9 @@ the beliefs at such information sets as being uniform across all member nodes.
 
 .. cmdoption:: -m
 
-   Stop when reaching the specified value of the
-   parameter lambda. By default, the tracing stops when lambda reaches
-   1,000,000, which is usually suitable for computing a good
-   approximation to a Nash equilibrium. For applications, such as to
-   laboratory experiments, where the behavior of the correspondence for
-   small values of lambda is of interest and the asymptotic behavior is
-   not relevant, setting MAXLAMBDA to a much smaller value may be
-   indicated.
+   .. versionadded:: 16.2.0
+
+   Specify the maximum regret criterion for acceptance as an approximate Nash equilibrium.
 
 .. cmdoption:: -l
 

src/pygambit/gambit.pxd

@@ -455,10 +455,10 @@ cdef extern from "solvers/gnm/gnm.h":
     ) except +RuntimeError
 
 cdef extern from "solvers/logit/nfglogit.h":
-    c_List[c_MixedStrategyProfileDouble] LogitStrategySolve(c_Game) except +RuntimeError
+    c_List[c_MixedStrategyProfileDouble] LogitStrategySolve(c_Game, double) except +RuntimeError
 
 cdef extern from "solvers/logit/efglogit.h":
-    c_List[c_MixedBehaviorProfileDouble] LogitBehaviorSolve(c_Game) except +RuntimeError
+    c_List[c_MixedBehaviorProfileDouble] LogitBehaviorSolve(c_Game, double) except +RuntimeError
 
 cdef extern from "solvers/logit/nfglogit.h":
     cdef cppclass c_LogitQREMixedStrategyProfile "LogitQREMixedStrategyProfile":

src/pygambit/nash.h

@@ -53,12 +53,11 @@ std::shared_ptr<LogitQREMixedStrategyProfile> logit_atlambda(const Game &p_game,
       alg.SolveAtLambda(start, null_stream, p_lambda, 1.0));
 }
 
-List<LogitQREMixedStrategyProfile> logit_principal_branch(const Game &p_game,
-                                                          double p_maxLambda = 1000000.0)
+List<LogitQREMixedStrategyProfile> logit_principal_branch(const Game &p_game, double p_maxregret)
 {
   LogitQREMixedStrategyProfile start(p_game);
   StrategicQREPathTracer alg;
   NullBuffer null_buffer;
   std::ostream null_stream(&null_buffer);
-  return alg.TraceStrategicPath(start, null_stream, p_maxLambda, 1.0);
+  return alg.TraceStrategicPath(start, null_stream, p_maxregret, 1.0);
 }

src/pygambit/nash.pxi

@@ -33,6 +33,7 @@ def _convert_mspd(
 ) -> typing.List[MixedStrategyProfileDouble]:
     ret = []
     for i in range(inlist.Length()):
+        print(i)
         p = MixedStrategyProfileDouble()
         p.profile = copyitem_list_mspd(inlist, i+1)
         ret.append(p)
@@ -182,12 +183,12 @@ def _gnm_strategy_solve(
         raise
 
 
-def _logit_strategy_solve(game: Game) -> typing.List[MixedStrategyProfileDouble]:
-    return _convert_mspd(LogitStrategySolve(game.game))
+def _logit_strategy_solve(game: Game, maxregret: float) -> typing.List[MixedStrategyProfileDouble]:
+    return _convert_mspd(LogitStrategySolve(game.game, maxregret))
 
 
-def _logit_behavior_solve(game: Game) -> typing.List[MixedBehaviorProfileDouble]:
-    return _convert_mbpd(LogitBehaviorSolve(game.game))
+def _logit_behavior_solve(game: Game, maxregret: float) -> typing.List[MixedBehaviorProfileDouble]:
+    return _convert_mbpd(LogitBehaviorSolve(game.game, maxregret))
 
 
 @cython.cclass
@@ -254,8 +255,8 @@ def logit_atlambda(game: Game, lam: float) -> LogitQREMixedStrategyProfile:
     return ret
 
 
-def logit_principal_branch(game: Game, maxlam: float = 100000.0):
-    solns = _logit_principal_branch(game.game, maxlam)
+def logit_principal_branch(game: Game):
+    solns = _logit_principal_branch(game.game, 1.0e-6)
     ret = []
     for i in range(solns.Length()):
         p = LogitQREMixedStrategyProfile()

src/pygambit/nash.py

@@ -516,7 +516,9 @@ def gnm_solve(
 
 
 def logit_solve(
-        game: libgbt.Game, use_strategic: bool = False
+        game: libgbt.Game,
+        use_strategic: bool = False,
+        maxregret: float = 0.0001,
 ) -> NashComputationResult:
     """Compute Nash equilibria of a game using :ref:`the logit quantal response
     equilibrium correspondence <gambit-logit>`.
@@ -528,19 +530,29 @@ def logit_solve(
     ----------
     game : Game
         The game to compute equilibria in.
+
     use_strategic : bool, default False
         Whether to use the strategic form.  If True, always uses the strategic
         representation even if the game's native representation is extensive.
 
+    maxregret : float, default 0.0001
+        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
+
+        .. versionadded: 16.2.0
+
     Returns
     -------
     res : NashComputationResult
         The result represented as a ``NashComputationResult`` object.
     """
+    if maxregret <= 0.0:
+        raise ValueError("logit_solve(): maxregret argument must be positive")
     if not game.is_tree or use_strategic:
-        equilibria = libgbt._logit_strategy_solve(game)
+        equilibria = libgbt._logit_strategy_solve(game, maxregret)
     else:
-        equilibria = libgbt._logit_behavior_solve(game)
+        equilibria = libgbt._logit_behavior_solve(game, maxregret)
     return NashComputationResult(
         game=game,
         method="logit",

src/solvers/logit/efglogit.cc

@@ -36,9 +36,12 @@ namespace Gambit {
 class AgentQREPathTracer::EquationSystem : public PathTracer::EquationSystem {
 public:
   explicit EquationSystem(const Game &p_game);
+
   ~EquationSystem() override;
+
   // Compute the value of the system of equations at the specified point.
   void GetValue(const Vector<double> &p_point, Vector<double> &p_lhs) const override;
+
   // Compute the Jacobian matrix at the specified point.
   void GetJacobian(const Vector<double> &p_point, Matrix<double> &p_matrix) const override;
 
@@ -49,7 +52,9 @@ class AgentQREPathTracer::EquationSystem : public PathTracer::EquationSystem {
   class Equation {
   public:
     virtual ~Equation() = default;
+
     virtual double Value(const LogBehavProfile<double> &p_point, double p_lambda) const = 0;
+
     virtual void Gradient(const LogBehavProfile<double> &p_point, double p_lambda,
                           Vector<double> &p_gradient) const = 0;
   };
@@ -72,6 +77,7 @@ class AgentQREPathTracer::EquationSystem : public PathTracer::EquationSystem {
     }
 
     double Value(const LogBehavProfile<double> &p_profile, double p_lambda) const override;
+
     void Gradient(const LogBehavProfile<double> &p_profile, double p_lambda,
                   Vector<double> &p_gradient) const override;
   };
@@ -95,6 +101,7 @@ class AgentQREPathTracer::EquationSystem : public PathTracer::EquationSystem {
     }
 
     double Value(const LogBehavProfile<double> &p_profile, double p_lambda) const override;
+
     void Gradient(const LogBehavProfile<double> &p_profile, double p_lambda,
                   Vector<double> &p_gradient) const override;
   };
@@ -238,9 +245,11 @@ class AgentQREPathTracer::CallbackFunction : public PathTracer::CallbackFunction
     : m_stream(p_stream), m_game(p_game), m_fullGraph(p_fullGraph), m_decimals(p_decimals)
   {
   }
+
   ~CallbackFunction() override = default;
 
   void operator()(const Vector<double> &p_point, bool p_isTerminal) const override;
+
   const List<LogitQREMixedBehaviorProfile> &GetProfiles() const { return m_profiles; }
 
 private:
@@ -302,27 +311,62 @@ class AgentQREPathTracer::LambdaCriterion : public PathTracer::CriterionFunction
 //              AgentQREPathTracer: Wrapper to the tracing engine
 //------------------------------------------------------------------------------
 
+namespace {
+
+bool RegretTerminationFunction(const Game &p_game, const Vector<double> &p_point, double p_regret)
+{
+  if (p_point.back() < 0.0) {
+    return false;
+  }
+  MixedBehaviorProfile<double> profile(p_game);
+  for (int i = 1; i < p_point.Length(); i++) {
+    profile[i] = exp(p_point[i]);
+  }
+  return profile.GetMaxRegret() < p_regret;
+}
+
+} // namespace
+
 List<LogitQREMixedBehaviorProfile>
 AgentQREPathTracer::TraceAgentPath(const LogitQREMixedBehaviorProfile &p_start,
-                                   std::ostream &p_stream, double p_maxLambda, double p_omega,
-                                   double p_targetLambda)
+                                   std::ostream &p_stream, double p_regret, double p_omega) const
 {
+  double scale = p_start.GetGame()->GetMaxPayoff() - p_start.GetGame()->GetMinPayoff();
+  if (scale != 0.0) {
+    p_regret *= scale;
+  }
+
   List<LogitQREMixedBehaviorProfile> ret;
   Vector<double> x(p_start.BehaviorProfileLength() + 1);
   for (size_t i = 1; i <= p_start.BehaviorProfileLength(); i++) {
     x[i] = log(p_start[i]);
   }
-  x[x.Length()] = p_start.GetLambda();
+  x.back() = p_start.GetLambda();
 
   CallbackFunction func(p_stream, p_start.GetGame(), m_fullGraph, m_decimals);
-  if (p_targetLambda > 0.0) {
-    TracePath(EquationSystem(p_start.GetGame()), x, p_maxLambda, p_omega, func,
-              LambdaCriterion(p_targetLambda));
-  }
-  else {
-    TracePath(EquationSystem(p_start.GetGame()), x, p_maxLambda, p_omega, func);
-  }
+  TracePath(
+      EquationSystem(p_start.GetGame()), x, p_omega,
+      [p_start, p_regret](const Vector<double> &p_point) {
+        return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret);
+      },
+      func);
   return func.GetProfiles();
 }
 
+LogitQREMixedBehaviorProfile
+AgentQREPathTracer::SolveAtLambda(const LogitQREMixedBehaviorProfile &p_start,
+                                  std::ostream &p_stream, double p_targetLambda,
+                                  double p_omega) const
+{
+  Vector<double> x(p_start.BehaviorProfileLength() + 1);
+  for (int i = 1; i <= p_start.BehaviorProfileLength(); i++) {
+    x[i] = log(p_start[i]);
+  }
+  x.back() = p_start.GetLambda();
+  CallbackFunction func(p_stream, p_start.GetGame(), m_fullGraph, m_decimals);
+  TracePath(EquationSystem(p_start.GetGame()), x, p_omega, LambdaPositiveTerminationFunction, func,
+            LambdaCriterion(p_targetLambda));
+  return func.GetProfiles().back();
+}
+
 } // end namespace Gambit

src/solvers/logit/efglogit.h

@@ -53,8 +53,11 @@ class AgentQREPathTracer : public PathTracer {
   ~AgentQREPathTracer() override = default;
 
   List<LogitQREMixedBehaviorProfile> TraceAgentPath(const LogitQREMixedBehaviorProfile &p_start,
-                                                    std::ostream &p_stream, double p_maxLambda,
-                                                    double p_omega, double p_targetLambda = -1.0);
+                                                    std::ostream &p_stream, double p_regret,
+                                                    double p_omega) const;
+  LogitQREMixedBehaviorProfile SolveAtLambda(const LogitQREMixedBehaviorProfile &p_start,
+                                             std::ostream &p_logStream, double p_targetLambda,
+                                             double p_omega) const;
 
   void SetFullGraph(bool p_fullGraph) { m_fullGraph = p_fullGraph; }
   bool GetFullGraph() const { return m_fullGraph; }
@@ -71,13 +74,13 @@ class AgentQREPathTracer : public PathTracer {
   class LambdaCriterion;
 };
 
-inline List<MixedBehaviorProfile<double>> LogitBehaviorSolve(const Game &p_game)
+inline List<MixedBehaviorProfile<double>> LogitBehaviorSolve(const Game &p_game, double p_epsilon)
 {
   AgentQREPathTracer tracer;
   tracer.SetFullGraph(false);
   std::ostringstream ostream;
   auto result =
-      tracer.TraceAgentPath(LogitQREMixedBehaviorProfile(p_game), ostream, 1000000.0, 1.0);
+      tracer.TraceAgentPath(LogitQREMixedBehaviorProfile(p_game), ostream, p_epsilon, 1.0);
   auto ret = List<MixedBehaviorProfile<double>>();
   ret.push_back(result[1].GetProfile());
   return ret;

src/solvers/logit/nfglogit.cc

@@ -209,18 +209,44 @@ void StrategicQREPathTracer::CallbackFunction::operator()(const Vector<double> &
 //               StrategicQREPathTracer: Main driver routines
 //----------------------------------------------------------------------------
 
+namespace {
+
+bool RegretTerminationFunction(const Game &p_game, const Vector<double> &p_point, double p_regret)
+{
+  if (p_point.back() < 0.0) {
+    return false;
+  }
+  MixedStrategyProfile<double> profile(p_game->NewMixedStrategyProfile(0.0));
+  for (int i = 1; i < p_point.Length(); i++) {
+    profile[i] = exp(p_point[i]);
+  }
+  return profile.GetMaxRegret() < p_regret;
+}
+
+} // namespace
+
 List<LogitQREMixedStrategyProfile>
 StrategicQREPathTracer::TraceStrategicPath(const LogitQREMixedStrategyProfile &p_start,
-                                           std::ostream &p_stream, double p_maxLambda,
+                                           std::ostream &p_stream, double p_regret,
                                            double p_omega) const
 {
+  double scale = p_start.GetGame()->GetMaxPayoff() - p_start.GetGame()->GetMinPayoff();
+  if (scale != 0.0) {
+    p_regret *= scale;
+  }
+
   Vector<double> x(p_start.MixedProfileLength() + 1);
   for (int i = 1; i <= p_start.MixedProfileLength(); i++) {
     x[i] = log(p_start[i]);
   }
-  x[x.Length()] = p_start.GetLambda();
+  x.back() = p_start.GetLambda();
   CallbackFunction func(p_stream, p_start.GetGame(), m_fullGraph, m_decimals);
-  TracePath(EquationSystem(p_start.GetGame()), x, p_maxLambda, p_omega, func);
+  TracePath(
+      EquationSystem(p_start.GetGame()), x, p_omega,
+      [p_start, p_regret](const Vector<double> &p_point) {
+        return RegretTerminationFunction(p_start.GetGame(), p_point, p_regret);
+      },
+      func);
   return func.GetProfiles();
 }
 
@@ -233,10 +259,10 @@ StrategicQREPathTracer::SolveAtLambda(const LogitQREMixedStrategyProfile &p_star
   for (int i = 1; i <= p_start.MixedProfileLength(); i++) {
     x[i] = log(p_start[i]);
   }
-  x[x.Length()] = p_start.GetLambda();
+  x.back() = p_start.GetLambda();
   CallbackFunction func(p_stream, p_start.GetGame(), m_fullGraph, m_decimals);
-  TracePath(EquationSystem(p_start.GetGame()), x, std::max(1.0, 3.0 * p_targetLambda), p_omega,
-            func, LambdaCriterion(p_targetLambda));
+  TracePath(EquationSystem(p_start.GetGame()), x, p_omega, LambdaPositiveTerminationFunction, func,
+            LambdaCriterion(p_targetLambda));
   return func.GetProfiles().back();
 }
 
@@ -380,19 +406,18 @@ StrategicQREEstimator::Estimate(const LogitQREMixedStrategyProfile &p_start,
   for (int i = 1; i <= p_start.MixedProfileLength(); i++) {
     x[i] = log(p_start[i]);
   }
-  x[x.Length()] = p_start.GetLambda();
+  x.back() = p_start.GetLambda();
 
   CallbackFunction callback(p_stream, p_start.GetGame(),
                             static_cast<const Vector<double> &>(p_frequencies), m_fullGraph,
                             m_decimals);
-  while (x[x.Length()] < p_maxLambda) {
-    TracePath(EquationSystem(p_start.GetGame()), x, p_maxLambda, p_omega, callback,
-              CriterionFunction(static_cast<const Vector<double> &>(p_frequencies)));
-    if (x[x.Length()] < p_maxLambda) {
-      // Found an extremum of the likelihood function
-      // start iterating again from the same point in case of
-      // local optima.
-    }
+  while (x.back() < p_maxLambda) {
+    TracePath(
+        EquationSystem(p_start.GetGame()), x, p_omega,
+        [p_maxLambda](const Vector<double> &p_point) {
+          return LambdaRangeTerminationFunction(p_point, 0, p_maxLambda);
+        },
+        callback, CriterionFunction(static_cast<const Vector<double> &>(p_frequencies)));
   }
   callback.PrintMaximizer();
   return callback.GetMaximizer();