Efficient calculation of payoff derivatives

src/games/game.h

@@ -402,6 +402,7 @@ class GameStrategyRep : public std::enable_shared_from_this<GameStrategyRep> {
   template <class T> friend class MixedStrategyProfile;
   template <class T> friend class TableMixedStrategyProfileRep;
   template <class T> friend class MixedBehaviorProfile;
+  friend struct TreeSequenceLayout;
 
   bool m_valid{true};
   GamePlayerRep *m_player;

src/games/gameagg.cc

@@ -20,6 +20,7 @@
 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 //
 
+#include <type_traits>
 #include <iostream>
 
 #include "games.h"
@@ -90,6 +91,7 @@ template <class T> class AGGMixedStrategyProfileRep : public MixedStrategyProfil
   T GetPayoff(int pl) const override;
   T GetPayoffDeriv(int pl, const GameStrategy &) const override;
   T GetPayoffDeriv(int pl, const GameStrategy &, const GameStrategy &) const override;
+  void GetPayoffDerivBlock(int pl, Vector<T> &p_values, Matrix<T> &p_derivs) const override;
 };
 
 template <class T> T AGGMixedStrategyProfileRep<T>::GetPayoff(int pl) const
@@ -168,6 +170,51 @@ T AGGMixedStrategyProfileRep<T>::GetPayoffDeriv(int pl, const GameStrategy &ps1,
   return g.aggPtr->getMixedPayoff(pl - 1, s);
 }
 
+// Each row of derivatives comes from a single leave-one-out convolution in the AGG, which is
+// available in floating point only; exact computation uses the default.
+template <class T>
+void AGGMixedStrategyProfileRep<T>::GetPayoffDerivBlock(int pl, Vector<T> &p_values,
+                                                        Matrix<T> &p_derivs) const
+{
+  if constexpr (!std::is_same_v<T, double>) {
+    MixedStrategyProfileRep<T>::GetPayoffDerivBlock(pl, p_values, p_derivs);
+  }
+  else {
+    const auto game = this->m_support.GetGame();
+    auto &g = dynamic_cast<GameAGGRep &>(*game);
+    const auto &agg = *g.aggPtr;
+    std::vector<double> s(agg.getNumActions());
+    for (int i = 0; i < agg.getNumPlayers(); ++i) {
+      for (int j = 0; j < agg.getNumActions(i); ++j) {
+        const int ind = this->m_profileIndex.at(game->GetPlayer(i + 1)->GetStrategy(j + 1));
+        s[agg.firstAction(i) + j] = (ind == -1) ? 0.0 : this->m_probs.GetFlattened()[ind];
+      }
+    }
+
+    const auto player = game->GetPlayer(pl);
+    const auto &strategies = this->m_support.GetStrategies(player);
+    p_values = Vector<double>(strategies.size());
+    p_derivs = Matrix<double>(strategies.size(), this->m_probs.GetFlattened().size());
+    p_derivs = 0.0;
+    std::vector<double> row(agg.getNumActions());
+    size_t r = 1;
+    for (const auto &strategy : strategies) {
+      p_values[r] = GetPayoffDeriv(pl, strategy);
+      g.aggPtr->getPayoffJacobianRow(pl - 1, strategy->GetNumber() - 1, s, row);
+      for (const auto &other : game->GetPlayers()) {
+        if (other == player) {
+          continue;
+        }
+        for (const auto &otherStrategy : this->m_support.GetStrategies(other)) {
+          p_derivs(r, this->m_profileIndex.at(otherStrategy)) =
+              row[agg.firstAction(other->GetNumber() - 1) + otherStrategy->GetNumber() - 1];
+        }
+      }
+      r++;
+    }
+  }
+}
+
 template class AGGMixedStrategyProfileRep<double>;
 template class AGGMixedStrategyProfileRep<Rational>;
 

src/games/gamebagg.cc

@@ -20,7 +20,9 @@
 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
 //
 
+#include <algorithm>
 #include <iostream>
+#include <type_traits>
 
 #include "games.h"
 #include "gamebagg.h"
@@ -93,6 +95,10 @@ template <class T> class BAGGMixedStrategyProfileRep : public MixedStrategyProfi
   T GetPayoff(int pl) const override;
   T GetPayoffDeriv(int pl, const GameStrategy &) const override;
   T GetPayoffDeriv(int pl, const GameStrategy &, const GameStrategy &) const override;
+  void GetPayoffDerivBlock(int pl, Vector<T> &p_values, Matrix<T> &p_derivs) const override;
+
+private:
+  T GetPartialDeriv(int pl, const std::vector<GameStrategy> &p_strategies) const;
 };
 
 template <class T> T BAGGMixedStrategyProfileRep<T>::GetPayoff(int pl) const
@@ -119,84 +125,168 @@ template <class T> T BAGGMixedStrategyProfileRep<T>::GetPayoff(int pl) const
   return g.baggPtr->getMixedPayoff(bplayer, btype, s);
 }
 
+// The payoff to agent pl is linear in its own mixed strategy and in the mixed strategy each
+// other BAGG player induces on its AGG actions, to which each of that player's types t
+// contributes with weight P(t); it does not depend on the strategies of the other types of pl's
+// own BAGG player.  The partial derivative with respect to the probabilities of p_strategies is
+// therefore the payoff when each is played for sure (for another BAGG player, by all of its
+// types), times the probabilities of the types of the other BAGG players involved.  It is zero
+// if two of the strategies belong to one BAGG player, or one belongs to another type of pl's.
 template <class T>
-T BAGGMixedStrategyProfileRep<T>::GetPayoffDeriv(int pl, const GameStrategy &ps) const
+T BAGGMixedStrategyProfileRep<T>::GetPartialDeriv(
+    int pl, const std::vector<GameStrategy> &p_strategies) const
 {
-  auto &g = dynamic_cast<GameBAGGRep &>(*(this->m_support.GetGame()));
-  std::vector<T> s(g.GetStrategies().size());
-  int bplayer = -1, btype = -1;
-  for (int i = 0; i < g.baggPtr->getNumPlayers(); ++i) {
-    for (int tp = 0; tp < g.baggPtr->getNumTypes(i); ++tp) {
-      if (pl == g.baggPtr->typeOffset[i] + tp + 1) {
-        bplayer = i;
-        btype = tp;
-      }
-      if (g.baggPtr->typeOffset[i] + tp + 1 == ps->GetPlayer()->GetNumber()) {
-        for (unsigned int j = 0; j < g.baggPtr->typeActionSets.at(i).at(tp).size(); ++j) {
-          s.at(g.baggPtr->firstAction(i, tp) + j) = T(0);
-        }
-        s.at(g.baggPtr->firstAction(i, tp) + ps->GetNumber() - 1) = T(1);
-      }
-      else {
-        for (int j = 0; j < g.baggPtr->getNumActions(i, tp); ++j) {
-          const GameStrategy strategy = this->m_support.GetGame()
-                                            ->GetPlayer(g.baggPtr->typeOffset[i] + tp + 1)
-                                            ->GetStrategy(j + 1);
-          const int ind = this->m_profileIndex.at(strategy);
-          s.at(g.baggPtr->firstAction(i, tp) + j) =
-              (ind == -1) ? T(0) : this->m_probs.GetFlattened()[ind];
-        }
+  const auto game = this->m_support.GetGame();
+  auto &g = dynamic_cast<GameBAGGRep &>(*game);
+  const auto &bagg = *g.baggPtr;
+  auto &agg = *bagg.aggPtr;
+  const int bplayer = g.agent2baggPlayer[pl] - 1;
+  const int btype = pl - 1 - bagg.typeOffset[bplayer];
+
+  int ownAction = -1;
+  // Each other BAGG player involved, with the AGG action it plays for sure
+  std::vector<std::pair<int, int>> fixed;
+  T weight(1);
+  for (const auto &strategy : p_strategies) {
+    const int agent = strategy->GetPlayer()->GetNumber();
+    const int player = g.agent2baggPlayer[agent] - 1;
+    const int type = agent - 1 - bagg.typeOffset[player];
+    if (agent == pl) {
+      if (ownAction != -1) {
+        return T(0);
       }
+      ownAction = strategy->GetNumber() - 1;
+      continue;
+    }
+    if (player == bplayer || std::any_of(fixed.begin(), fixed.end(),
+                                         [&](const auto &f) { return f.first == player; })) {
+      return T(0);
+    }
+    fixed.emplace_back(player,
+                       agg.firstAction(player) +
+                           bagg.typeAction2ActionIndex[player][type][strategy->GetNumber() - 1]);
+    if constexpr (std::is_same_v<T, Rational>) {
+      weight *= static_cast<Rational>(bagg.exactIndepTypeDist[player][type]);
+    }
+    else {
+      weight *= bagg.indepTypeDist[player][type];
     }
   }
-  return g.baggPtr->getMixedPayoff(bplayer, btype, s);
+
+  std::vector<T> s(bagg.strategyOffset[bagg.getNumTypes()]);
+  for (int k = 1, offset = 0; k <= bagg.getNumTypes(); k++) {
+    for (const auto &strategy : game->GetPlayer(k)->GetStrategies()) {
+      const int ind = this->m_profileIndex.at(strategy);
+      s[offset++] = (ind == -1) ? T(0) : this->m_probs.GetFlattened()[ind];
+    }
+  }
+
+  std::vector<T> as(agg.getNumActions());
+  T value(0);
+  for (int action = 0; action < bagg.getNumActions(bplayer, btype); action++) {
+    const T prob = (ownAction == -1) ? s[bagg.firstAction(bplayer, btype) + action]
+                                     : T((action == ownAction) ? 1 : 0);
+    if (prob <= T(0)) {
+      continue;
+    }
+    bagg.getAGGStrat(as, s, bplayer, btype, action);
+    for (const auto &[player, index] : fixed) {
+      std::fill(as.begin() + agg.firstAction(player),
+                as.begin() + agg.firstAction(player) + agg.getNumActions(player), T(0));
+      as[index] = T(1);
+    }
+    value += prob * agg.getV(bplayer, bagg.typeAction2ActionIndex[bplayer][btype][action], as);
+  }
+  return (fixed.empty()) ? value : weight * value;
+}
+
+template <class T>
+T BAGGMixedStrategyProfileRep<T>::GetPayoffDeriv(int pl, const GameStrategy &ps) const
+{
+  return GetPartialDeriv(pl, {ps});
 }
 
 template <class T>
 T BAGGMixedStrategyProfileRep<T>::GetPayoffDeriv(int pl, const GameStrategy &ps1,
                                                  const GameStrategy &ps2) const
 {
-  const auto player1 = ps1->GetPlayer().get();
-  const auto player2 = ps2->GetPlayer().get();
-  if (player1 == player2) {
-    return T(0);
-  }
+  return GetPartialDeriv(pl, {ps1, ps2});
+}
 
-  auto &g = dynamic_cast<GameBAGGRep &>(*(this->m_support.GetGame()));
-  std::vector<T> s(g.GetStrategies().size());
-  int bplayer = -1, btype = -1;
-  for (int i = 0; i < g.baggPtr->getNumPlayers(); ++i) {
-    for (int tp = 0; tp < g.baggPtr->getNumTypes(i); ++tp) {
-      if (pl == g.baggPtr->typeOffset[i] + tp + 1) {
-        bplayer = i;
-        btype = tp;
+// See GetPartialDeriv.  The derivative with respect to an action b of type tq of another BAGG
+// player j is P(tq) times the payoff when j plays b's AGG action for sure, which is one entry of
+// a row of the AGG payoff Jacobian.  The Jacobian row is available in floating point only; exact
+// computation uses the default.
+template <class T>
+void BAGGMixedStrategyProfileRep<T>::GetPayoffDerivBlock(int pl, Vector<T> &p_values,
+                                                         Matrix<T> &p_derivs) const
+{
+  if constexpr (!std::is_same_v<T, double>) {
+    MixedStrategyProfileRep<T>::GetPayoffDerivBlock(pl, p_values, p_derivs);
+  }
+  else {
+    const auto game = this->m_support.GetGame();
+    const auto &bagg = *dynamic_cast<GameBAGGRep &>(*game).baggPtr;
+    auto &agg = *bagg.aggPtr;
+    const int numAgents = bagg.getNumTypes();
+
+    // The BAGG strategy vector, and the profile position of each of its entries (-1 if the
+    // strategy is not in the support); agent k's strategies start at strategyOffset[k - 1]
+    std::vector<double> s(bagg.strategyOffset[numAgents]);
+    std::vector<int> column(s.size());
+    for (int k = 1; k <= numAgents; k++) {
+      int offset = bagg.strategyOffset[k - 1];
+      for (const auto &strategy : game->GetPlayer(k)->GetStrategies()) {
+        const int ind = this->m_profileIndex.at(strategy);
+        s[offset] = (ind == -1) ? 0.0 : this->m_probs.GetFlattened()[ind];
+        column[offset++] = ind;
       }
+    }
 
-      if (g.baggPtr->typeOffset[i] + tp + 1 == player1->GetNumber()) {
-        for (unsigned int j = 0; j < g.baggPtr->typeActionSets.at(i).at(tp).size(); ++j) {
-          s.at(g.baggPtr->firstAction(i, tp) + j) = T(0);
-        }
-        s.at(g.baggPtr->firstAction(i, tp) + ps1->GetNumber() - 1) = T(1);
-      }
-      else if (g.baggPtr->typeOffset[i] + tp + 1 == player2->GetNumber()) {
-        for (int j = 0; j < g.baggPtr->getNumActions(i, tp); ++j) {
-          s.at(g.baggPtr->firstAction(i, tp) + j) = T(0);
+    int bplayer = 0;
+    while (bagg.typeOffset[bplayer + 1] < pl) {
+      bplayer++;
+    }
+    const int btype = pl - 1 - bagg.typeOffset[bplayer];
+    std::vector<double> as(agg.getNumActions());
+    bagg.getAGGStrat(as, s, bplayer, btype, 0);
+    const int ownFirst = agg.firstAction(bplayer);
+    const auto &ownActions = bagg.typeAction2ActionIndex[bplayer][btype];
+    as[ownFirst + ownActions[0]] = 0.0;
+
+    const auto &strategies = this->m_support.GetStrategies(game->GetPlayer(pl));
+    p_values = Vector<double>(strategies.size());
+    p_derivs = Matrix<double>(strategies.size(), this->m_probs.GetFlattened().size());
+    p_derivs = 0.0;
+    std::vector<double> row(agg.getNumActions());
+    size_t r = 1;
+    for (const auto &strategy : strategies) {
+      const int action = ownActions[strategy->GetNumber() - 1];
+      as[ownFirst + action] = 1.0;
+      const double value = agg.getV(bplayer, action, as);
+      p_values[r] = value;
+      agg.getPayoffJacobianRow(bplayer, action, as, row);
+      as[ownFirst + action] = 0.0;
+
+      for (int j = 0; j < bagg.getNumPlayers(); j++) {
+        if (j == bplayer) {
+          continue;
         }
-        s.at(g.baggPtr->firstAction(i, tp) + ps2->GetNumber() - 1) = T(1);
-      }
-      else {
-        for (unsigned int j = 0; j < g.baggPtr->typeActionSets.at(i).at(tp).size(); ++j) {
-          const GameStrategy strategy = this->m_support.GetGame()
-                                            ->GetPlayer(g.baggPtr->typeOffset[i] + tp + 1)
-                                            ->GetStrategy(j + 1);
-          const int ind = this->m_profileIndex.at(strategy);
-          s.at(g.baggPtr->firstAction(i, tp) + j) =
-              (ind == -1) ? T(0) : this->m_probs.GetFlattened()[ind];
+        const double *payoffs = row.data() + agg.firstAction(j);
+        for (int tq = 0; tq < bagg.getNumTypes(j); tq++) {
+          const int first = bagg.strategyOffset[bagg.typeOffset[j] + tq];
+          const auto &actions = bagg.typeAction2ActionIndex[j][tq];
+          const double weight = bagg.indepTypeDist[j][tq];
+          for (int b = 0; b < bagg.getNumActions(j, tq); b++) {
+            if (column[first + b] != -1) {
+              p_derivs(r, column[first + b]) = weight * payoffs[actions[b]];
+            }
+          }
         }
       }
+      r++;
     }
   }
-  return g.baggPtr->getMixedPayoff(bplayer, btype, s);
 }
 
 template class BAGGMixedStrategyProfileRep<double>;