diff --git a/src/games/agg/agg.cc b/src/games/agg/agg.cc index 05be8d684e..1767d38041 100644 --- a/src/games/agg/agg.cc +++ b/src/games/agg/agg.cc @@ -36,16 +36,15 @@ inline int select2nd(const pair &x) { return x.second; } AGG::AGG(int numPlayers, std::vector &_actions, int numANodes, int _numPNodes, vector> &_actionSets, vector> &neighb, - vector &projTypes, vector>> &projS, - vector>> &proj, vector> &projF, - vector>> &Po, vector &_payoffs, - vector &_exactPayoffs) + vector &projTypes, vector>> &proj, + vector> &projF, vector>> &Po, + vector &_payoffs, vector &_exactPayoffs) : numPlayers(numPlayers), numActionNodes(numANodes), numPNodes(_numPNodes), actionSets(_actionSets), neighbors(neighb), projectionTypes(projTypes), payoffs(_payoffs), exactPayoffs(_exactPayoffs), projection(proj), m_state(numANodes, numPlayers), - m_exactState(numANodes, numPlayers), fullProjectedStrat(projS), projFunctions(projF), - Porder(Po), isPure(numANodes, true), node2Action(numANodes, vector(numPlayers)), - cache(numPlayers + 1), player2Class(numPlayers), kSymStrategyOffset(1, 0) + m_exactState(numANodes, numPlayers), projFunctions(projF), Porder(Po), isPure(numANodes, true), + node2Action(numANodes, vector(numPlayers)), player2Class(numPlayers), + kSymStrategyOffset(1, 0) { // actions @@ -293,8 +292,8 @@ std::shared_ptr AGG::makeAGG(istream &in) throw std::runtime_error("Unknown payoff type " + std::to_string(t)); } } - return std::make_shared(n, size, S, P, ASets, neighb, projTypes, projS, proj, projF, Po, - pays, exactPays); + return std::make_shared(n, size, S, P, ASets, neighb, projTypes, proj, projF, Po, pays, + exactPays); } void AGG::setProjections(vector>> &projS, @@ -392,8 +391,10 @@ void AGG::initPorder(vector &Po, int i, int N, vector()) varies with V, the -// action graph structure itself is shared. +// compute the induced distribution. The action graph structure (projection, Porder, +// projFunctions, actionSets, neighbors) is purely combinatorial -- built once from integer data, +// independent of payoff/probability values -- so it's shared by both V; only the numeric working +// state, reached via state(), differs by V. template void AGG::computeP(int player, int act, int player2, int act2) { auto &s = state(); @@ -518,6 +519,114 @@ V AGG::getJ(int player1, int act1, int player2, int act2, const StrategyProfile< return payoffInnerProd(actionSets[player1][act1], state().Pr[numPlayers - 1]); } +namespace { + +// Fills allBut[indices[i]] with the convolution of factors[indices[k]] for all k in [start,end) +// except i, for every i in [start,end); temp receives the convolution of factors[indices[k]] +// for ALL k in [start,end) (the "full product" of this range), for the caller to combine with a +// sibling range. +void computeAllButOne(const std::vector &indices, size_t start, size_t end, + const std::vector> &factors, int numNei, + std::vector &projFuncs, + std::vector> &allBut, + ConfigDistribution &temp) +{ + if (end - start == 1) { + allBut[indices[start]].reset(); + return; + } + const size_t mid = start + (end - start) / 2; + + computeAllButOne(indices, start, mid, factors, numNei, projFuncs, allBut, temp); + computeAllButOne(indices, mid, end, factors, numNei, projFuncs, allBut, temp); + + // temp := full product of [start, mid) + temp = factors[indices[start]]; + if (mid - start > 1) { + temp.multiply(allBut[indices[start]], numNei, projFuncs); + } + + // extend each i in [start, mid)'s "all but i, within [start,mid)" by the full product of + // [mid, end), giving "all but i, within [start,end)" + if (mid - start == 1) { + allBut[indices[start]] = factors[indices[mid]]; + if (end - mid > 1) { + allBut[indices[start]].multiply(allBut[indices[mid]], numNei, projFuncs); + } + } + else { + for (size_t i = start; i < mid; ++i) { + allBut[indices[i]].multiply(factors[indices[mid]], numNei, projFuncs); + if (end - mid > 1) { + allBut[indices[i]].multiply(allBut[indices[mid]], numNei, projFuncs); + } + } + } + + // symmetrically, extend each i in [mid, end) by the full product of [start, mid), i.e. temp + if (end - mid == 1) { + allBut[indices[mid]] = temp; + } + else { + for (size_t i = mid; i < end; ++i) { + allBut[indices[i]].multiply(temp, numNei, projFuncs); + } + } +} + +} // namespace + +void AGG::getPayoffJacobianRow(int player1, int act1, const StrategyProfile &s, + StrategyProfile &dest) +{ + const int Node = actionSets[player1][act1]; + const int numNei = neighbors[Node].size(); + + // base[player1] is pinned to the one-hot at act1 (fixed for the whole row); base[k] for + // every other player k is k's actual projected mixed strategy at Node, computed fresh from s. + std::vector> base(numPlayers); + base[player1].insert(std::make_pair(projection[Node][player1][act1], 1.0)); + + std::vector others; + others.reserve(numPlayers - 1); + for (int k = 0; k < numPlayers; ++k) { + if (k == player1) { + continue; + } + others.push_back(k); + for (int j = 0; j < actions[k]; ++j) { + if (s[j + firstAction(k)] > 0.0) { + base[k] += std::make_pair(projection[Node][k][j], s[j + firstAction(k)]); + } + } + } + if (others.empty()) { + return; // only one player in the game; no off-diagonal entries to fill + } + + std::vector> allBut(numPlayers); + if (others.size() == 1) { + // With only one other player, "everyone except player1 and that player" is nobody, so the + // induced distribution is just player1's own pinned action. + allBut[others[0]] = base[player1]; + } + else { + ConfigDistribution temp; + computeAllButOne(others, size_t{0}, others.size(), base, numNei, projFunctions[Node], allBut, + temp); + for (const int j : others) { + allBut[j].multiply(base[player1], numNei, projFunctions[Node]); + } + } + + for (const int j : others) { + for (int act2 = 0; act2 < actions[j]; ++act2) { + dest[act2 + firstAction(j)] = allBut[j].inner_prod(projection[Node][j][act2], numNei, + projFunctions[Node], payoffs[Node]); + } + } +} + template double AGG::payoffInnerProd(int node, const ConfigDistribution &dist); template Rational AGG::payoffInnerProd(int node, const ConfigDistribution &dist); diff --git a/src/games/agg/agg.h b/src/games/agg/agg.h index 8d0d42e65b..3dfa5b2b92 100644 --- a/src/games/agg/agg.h +++ b/src/games/agg/agg.h @@ -33,39 +33,50 @@ #include "trie_map.h" #include "games/number.h" -namespace Gambit { +namespace Gambit::agg { -namespace gametracer { -class aggame; -} - -namespace agg { - -// mixed strategy profile: one probability per action; V is double or Rational. +// data structure for a mixed strategy profile: one probability per action, using numeric type V +// (double for the floating-point engine, Rational for the exact-arithmetic counterpart). template using StrategyProfile = std::vector; -// payoff function for one action node (neighbor configuration -> payoff). Not a template like -// ConfigDistribution below: it's consumed via trie_map's own (floating-point) inner_prod. +// the local payoff function for one action node: maps a neighbor configuration to the payoff of +// playing that node under it. Populated at parse time; consumed by the floating-point +// convolution engine (getMixedPayoff/getV/getJ) via trie_map's own inner_prod, so it has to stay +// a trie_map -- unlike StrategyProfile/ConfigDistribution above, this is deliberately not +// a template (see ExactPayoffTable below for why a Rational instantiation would be unsound). using PayoffTable = trie_map; -// exact counterpart of PayoffTable. A std::map, not a trie_map: it stores Number -// (not Rational) to preserve the original source text, and trie_map's inner_prod has an -// epsilon-discard threshold that is unsound for exact arithmetic. +// exact (arbitrary-precision) counterpart of PayoffTable, populated alongside it at parse time +// from the same source text. Deliberately not a trie_map, for two independent reasons: +// (1) trie_map's own inner_prod has an epsilon-tolerance discard/throw policy +// (trie_map::THRESH = 1e-12) that is meaningful for double but unsound for exact arithmetic -- +// there is no notion of "too small to matter" once payoffs are exact; and (2) this stores Number, +// not Rational, to preserve the original entered text for exact serialization (printPayoffs), +// which a bare Rational cannot do. As a result this is only ever used for storage/lookup +// (find/insert/iterate) and a hand-written inner product (AGG::exactInnerProd), never trie_map's +// own arithmetic. using ExactPayoffTable = std::map, Number>; -// probability distribution over configurations; V is double or Rational. Templated (unlike -// PayoffTable/ExactPayoffTable) because it participates in convolution arithmetic. +// data struct for a probability distribution over configurations, using numeric type V (double +// for the floating-point engine, Rational for the exact-arithmetic counterpart). Unlike +// PayoffTable/ExactPayoffTable above, this participates in the convolution algorithm +// (multiply/inner_prod), so V must support arithmetic -- which Rational, unlike Number, +// provides -- letting both instantiations reuse trie_map directly. template using ConfigDistribution = trie_map; -// AGG's mutable working state for the convolution algorithm. AGG keeps one instance per V (see -// AGG::state()) since double and Rational profiles may be in use simultaneously. +// AGG's mutable working state for the convolution algorithm, at numeric type V. AGG holds one +// instance per V (see AGG::state()) rather than a single instance of one or the other, since both +// may be in use simultaneously -- e.g. from different MixedStrategyProfile instantiations over +// the same AGG -- and each is rebuilt fresh per call regardless of its previous contents. template struct ConvolutionState { - // for each action node and player, the projected mixed strategy (a distribution over - // 'contributions' -- see AGG::projection). + // foreach s \in S, foreach i \in N, the projected mixed strat: a prob distribution over the + // set of 'contributions' (see AGG::projection). std::vector>> projectedStrat; - // used by AGG::computeP(): the distribution induced by a partial strategy profile, or (for - // the payoff jacobian) by all agents except one. + // when computing the induced distribution via AGG::computeP(): foreach k <= n-1, the prob. + // distribution P_k induced by the partial strat profile of agents o_1..o_k. Also reused, when + // computing the partial distributions for the payoff jacobian, as the partial distribution + // induced by all agents except j. std::vector> Pr; ConvolutionState(int numActionNodes, int numPlayers) @@ -74,8 +85,8 @@ template struct ConvolutionState { { } - // set projectedStrat[Node][player] to the weighted sum of contributions[j] over actions with - // s[j] > 0. + // reset projectedStrat[Node][player] to the weighted sum, over actions j with s[j] > 0, of + // contributions[j] weighted by s[j]. void project(int Node, int player, const std::vector> &contributions, int numActions, const V *s) { @@ -102,18 +113,14 @@ class AGG { static const char LBRACKET = '['; static const char RBRACKET = ']'; - friend class gametracer::aggame; // wrapper class for gametracer - // read an AGG from input stream static std::shared_ptr makeAGG(std::istream &in); - // constructor; projS seeds fullProjectedStrat only, not the convolution engine's own working - // state, which is sized independently. + // constructor. The convolution engine's own working state (ConvolutionState) is sized from + // numANodes/numPlayers directly, not seeded from a caller-supplied array. AGG(int numPlayers, std::vector &actions, int numANodes, int numPNodes, std::vector> &actionSets, std::vector> &neighbors, - std::vector &projTypes, - std::vector>> &projS, - std::vector>> &proj, + std::vector &projTypes, std::vector>> &proj, std::vector> &projF, std::vector>> &Po, std::vector &payoffs, std::vector &exactPayoffs); @@ -122,6 +129,8 @@ class AGG { int getNumPlayers() const { return numPlayers; } int getNumActions() const { return totalActions; } int getNumActions(int i) const { return actions[i]; } + // per-player action counts, e.g. for constructing a gnmgame from an AGG. + const std::vector &getActionCounts() const { return actions; } int getMaxActions() const { return maxActions; } int firstAction(int i) const { return strategyOffset[i]; } int lastAction(int i) const { return strategyOffset[i + 1]; } @@ -135,15 +144,28 @@ class AGG { int firstKSymAction(int i) const { return kSymStrategyOffset[i]; } int lastKSymAction(int i) const { return kSymStrategyOffset[i + 1]; } - // expected payoff under a mixed profile, via the convolution algorithm. + // exp. payoff under mixed strat profile, via the convolution algorithm (doProjection/computeP, + // templated on V; inner_prod is trie_map's own method for V=double, exactInnerProd for + // V=Rational -- see exactInnerProd's comment for why that one isn't merged too). template V getMixedPayoff(int player, const StrategyProfile &s); void getPayoffVector(std::vector &dest, int player, const StrategyProfile &s); template V getV(int player, int action, const StrategyProfile &s); template V getJ(int player, int action, int player2, int action2, const StrategyProfile &s); - // payoff of the pure profile s: a direct table lookup. V is never deducible from s, so callers - // must specify it explicitly, e.g. getPurePayoff(...). + // Computes an entire row of the payoff Jacobian for player1 playing action1: for every other + // player j and every one of their actions, the payoff to player1 of playing action1 while j + // deviates to that action and everyone else plays according to s -- i.e. getJ(player1, + // action1, j, ., s) for every j and every one of j's actions, but computed for the whole row + // via a single leave-one-out convolution. double-only: its only caller (gametracer, for + // GNM/IPA) never needs Rational, and the implementation looks up PayoffTable directly rather + // than dispatching through payoffInnerProd, so it has no Rational specialization to offer. + void getPayoffJacobianRow(int player1, int action1, const StrategyProfile &s, + StrategyProfile &dest); + + // payoff of the pure profile s: a direct lookup, no floating-point computation involved. + // V=double looks up PayoffTable; V=Rational looks up ExactPayoffTable and casts the stored + // Number to Rational. template V getPurePayoff(int player, const std::vector &s) const; bool isSymmetric() const @@ -185,7 +207,8 @@ class AGG { const std::vector &getActionSet(int player) { return actionSets.at(player); } const PayoffTable &getPayoffMap(int node) { return payoffs.at(node); } - // largest/smallest payoff in any pure profile; callers write getMaxPayoff(), etc. + // largest/smallest payoff in any pure profile. Same V dispatch as getPurePayoff above; likewise + // never deducible (no arguments at all), so callers write getMaxPayoff(), etc. template V getMaxPayoff() const; template V getMinPayoff() const; @@ -247,13 +270,11 @@ class AGG { // the 'contribution' of s_i to D^(s) std::vector>> projection; + // working state for the convolution algorithm (see ConvolutionState above), one instance per + // numeric type V, since both may be in use simultaneously. ConvolutionState m_state; ConvolutionState m_exactState; - // full set of projected actions, filled once from the constructor's projS argument; read by - // gametracer::aggame for jacobian bookkeeping, not used by AGG's own convolution algorithm. - std::vector>> fullProjectedStrat; - // foreach s in S, foreach neighbor of s, its projection function std::vector> projFunctions; @@ -267,9 +288,6 @@ class AGG { // foreach s in S, j in N, the index of s in j's action set, or -1 if N/A std::vector> node2Action; - // cache of jacobian entries. - trie_map cache; - // the unique action sets std::vector uniqueActionSets; @@ -321,7 +339,9 @@ class AGG { static void initPorder(std::vector &Po, int i, int N, std::vector> &projS); - // dispatches to the V-specific working state, so computeP/doProjection don't need to branch. + // private methods. computeP/doProjection are the same algorithm for both V; they're made + // non-branching by going through state() for the V-specific working arrays (m_state vs + // m_exactState), rather than an if constexpr per call. template ConvolutionState &state() { if constexpr (std::is_same_v) { @@ -343,19 +363,22 @@ class AGG { template void doProjection(int Node, V *s); template void doProjection(int Node, int player, V *s); - // exactPayoffs is a std::map, not a trie_map, so this does the inner product by hand. + // exactPayoffs is a std::map<..., Number>, not a trie_map, so trie_map's own + // inner_prod can't be used directly against it (unlike the double path, which calls + // dist.inner_prod(payoffs[node]) directly) -- this does the equivalent sum manually. Rational exactInnerProd(int node, const ConfigDistribution &dist) const; - // dispatch point used by getV/getJ: V=double uses trie_map's own inner_prod, V=Rational uses - // exactInnerProd. + // dispatch point used by getV/getJ: V=double calls trie_map's own inner_prod against + // PayoffTable, V=Rational calls exactInnerProd against ExactPayoffTable. Not merged further + // than this one function -- same reason PayoffTable/ExactPayoffTable themselves aren't a + // template pair (see their comments): the two sides genuinely differ in container and value + // type, not just in V. template V payoffInnerProd(int node, const ConfigDistribution &dist); void getSymConfigProb(int plClass, StrategyProfile &s, int ownPlClass, int act, ConfigDistribution &dest, int plClass2 = -1, int act2 = -1); }; -} // namespace agg - -} // end namespace Gambit +} // namespace Gambit::agg #endif // GAMBIT_AGG_AGG_H diff --git a/src/games/agg/bagg.cc b/src/games/agg/bagg.cc index 80b07a824b..c938c4f668 100644 --- a/src/games/agg/bagg.cc +++ b/src/games/agg/bagg.cc @@ -269,7 +269,11 @@ template Rational BAGG::getMixedPayoff(int player, int tp, const StrategyProfile &s) const; // payoff of the pure profile ps, via the degenerate mixed profile and the convolution engine -// (mirrors AGG::getPurePayoff). +// (mirrors AGG::getPurePayoff). Deliberately does NOT enumerate the Cartesian product of the +// other players' types directly: that product grows with the number of OTHER PLAYERS, not just +// their type counts, defeating the whole point of the action-graph representation (efficient +// computation for many-player games via the graph structure). getMixedPayoff's convolution +// algorithm computes the same expectation in time polynomial in the number of players. template V BAGG::getPurePayoff(int player, int tp, const std::vector &ps) const { StrategyProfile s(strategyOffset[typeOffset[numPlayers]]); diff --git a/src/games/agg/bagg.h b/src/games/agg/bagg.h index 642947dd03..6cdc08ac08 100644 --- a/src/games/agg/bagg.h +++ b/src/games/agg/bagg.h @@ -71,13 +71,16 @@ class BAGG { double getMixedPayoff(int player, StrategyProfile &s); // exp payoff for player, conditioned on her receiving type tp, via the convolution algorithm. + // V=double runs the floating-point engine; V=Rational runs the exact-arithmetic counterpart + // throughout -- see agg.h::getMixedPayoff for how the dispatch works. template V getMixedPayoff(int player, int tp, const StrategyProfile &s) const; void getPayoffVector(std::vector &dest, int player, int tp, const StrategyProfile &s); template V getV(int player, int tp, int action, const StrategyProfile &s) const; - // payoff for player, conditioned on receiving type tp, under the pure profile s. + // payoff for player, conditioned on receiving type tp, under the pure profile s, via the + // degenerate mixed profile and the convolution algorithm (mirrors AGG::getPurePayoff). template V getPurePayoff(int player, int tp, const std::vector &s) const; double getPurePayoff(int player, std::vector &s) { diff --git a/src/games/gameagg.cc b/src/games/gameagg.cc index 25fc2c88a6..1d5050fa73 100644 --- a/src/games/gameagg.cc +++ b/src/games/gameagg.cc @@ -215,6 +215,9 @@ MixedStrategyProfile GameAGGRep::NewMixedStrategyProfile(double) const MixedStrategyProfile GameAGGRep::NewMixedStrategyProfile(const Rational &) const { + // AGG supports exact payoff computation on Rational profiles throughout, via the same + // convolution algorithm as the double engine (getMixedPayoff et al.); see + // AGGMixedStrategyProfileRep::GetPayoff. return MixedStrategyProfile(std::make_unique>( StrategySupportProfile(std::const_pointer_cast(shared_from_this())))); } diff --git a/src/games/gamebagg.cc b/src/games/gamebagg.cc index 53041a73e1..039b913cb9 100644 --- a/src/games/gamebagg.cc +++ b/src/games/gamebagg.cc @@ -287,6 +287,9 @@ MixedStrategyProfile GameBAGGRep::NewMixedStrategyProfile(double) const MixedStrategyProfile GameBAGGRep::NewMixedStrategyProfile(const Rational &) const { + // BAGG supports exact payoff computation on Rational profiles throughout, via the same + // convolution algorithm as the double engine (getMixedPayoff et al.); see + // BAGGMixedStrategyProfileRep::GetPayoff. return MixedStrategyProfile(std::make_unique>( StrategySupportProfile(std::const_pointer_cast(shared_from_this())))); } diff --git a/src/solvers/gtracer/aggame.cc b/src/solvers/gtracer/aggame.cc index 0ef1f43a14..006bc44451 100644 --- a/src/solvers/gtracer/aggame.cc +++ b/src/solvers/gtracer/aggame.cc @@ -25,380 +25,40 @@ namespace Gambit::gametracer { -void aggame::computePartialP_PureNode(int player1, int act1, std::vector &tasks) const -{ - int i, j; - const int Node = aggPtr->actionSets[player1][act1]; - const int numNei = aggPtr->neighbors[Node].size(); - - // assert(aggPtr->isPure[Node]||tasks.size()==0); - std::vector strat(numNei); - agg::AGG::Config a(numNei, 0); - // compute the full distrib - aggPtr->computeP(player1, act1); - - // store the full distrib in Pr[player1] - aggPtr->m_state.Pr[player1].swap(aggPtr->m_state.Pr[numPlayers - 1]); - for (i = 0; i < (int)tasks.size(); i++) { - // assert(tasks[i]!=player1); - agg::ConfigDistribution &P = aggPtr->m_state.Pr[tasks[i]]; - // P.clear(); // to get ready for division, we need clear() - P = aggPtr->m_state.Pr[player1]; - - bool NullOnly = true; - for (j = 0; j < numNei; ++j) { - a[j]++; - auto pp = aggPtr->m_state.projectedStrat[Node][tasks[i]].find(a); - if (pp == aggPtr->m_state.projectedStrat[Node][tasks[i]].end()) { - strat[j] = 0; - } - else { - strat[j] = pp->second; - if (strat[j] > (double)0) { - NullOnly = false; - } - } - - a[j]--; - } -#ifdef AGGDEBUG - cout << "dividing " << endl; - P.print_in_order(); - cout << endl << "by ["; - copy(strat.begin(), strat.end(), ostream_iterator(cout, " ")); - cout << "]\n"; -#endif - if (!NullOnly) { - P /= strat; - } -#ifdef AGGDEBUG - cout << "result is: " << endl; - P.print_in_order(); -#endif - } // end for(i -} - -void aggame::computePartialP_bisect(int player1, int act1, std::vector::iterator start, - std::vector::iterator endp, - agg::ConfigDistribution &temp) const -{ -#ifdef AGGDEBUG - cout << "calling computePartialP_bisect with player1=" << player1 << ", act1=" << act1 - << " *start=" << *start << " *(endp-1)=" << *(endp - 1) << ", (endp-start)=" << endp - start - << endl; -#endif - if (endp - start == 1) { - aggPtr->m_state.Pr[*start].reset(); - return; - } - const int Node = aggPtr->actionSets[player1][act1]; - const int numNei = aggPtr->neighbors[Node].size(); - - int player2; - std::vector::iterator ptr, mid = start + (endp - start) / 2; -#ifdef AGGDEBUG - cout << "*mid=" << *mid << " mid-start=" << mid - start << " endp-mid=" << endp - mid << endl; -#endif - computePartialP_bisect(player1, act1, start, mid, temp); - computePartialP_bisect(player1, act1, mid, endp, temp); - - temp.reset(); - temp = aggPtr->m_state.projectedStrat[Node][*start]; - if (mid - start > 1) { - temp.multiply(aggPtr->m_state.Pr[*start], numNei, aggPtr->projFunctions[Node]); - } - - if (mid - start == 1) { - // assert(aggPtr->m_state.Pr[*start].empty()); - aggPtr->m_state.Pr[*start] = aggPtr->m_state.projectedStrat[Node][*mid]; - if (endp - mid > 1) { - aggPtr->m_state.Pr[*start].multiply(aggPtr->m_state.Pr[*mid], numNei, - aggPtr->projFunctions[Node]); - } - } - else { - for (ptr = start; ptr != mid; ++ptr) { - player2 = *ptr; - aggPtr->m_state.Pr[player2].multiply(aggPtr->m_state.projectedStrat[Node][*mid], numNei, - aggPtr->projFunctions[Node]); - if (endp - mid > 1) { - aggPtr->m_state.Pr[player2].multiply(aggPtr->m_state.Pr[*mid], numNei, - aggPtr->projFunctions[Node]); - } - } - } - - if (endp - mid == 1) { - // assert(aggPtr->m_state.Pr[*mid].empty()); - aggPtr->m_state.Pr[*mid] = temp; - } - else { - for (ptr = mid; ptr != endp; ++ptr) { - player2 = *ptr; - aggPtr->m_state.Pr[player2].multiply(temp, numNei, aggPtr->projFunctions[Node]); - } - } -} - +// The payoff Jacobian: dest(i, j) is the payoff to the owner of action i if it deviates to i +// while the owner of action j simultaneously deviates to j, everyone else playing s. Computed a +// full row at a time via AGG::getPayoffJacobianRow -- see its comment for the algorithm -- so +// this only touches AGG's public API and avoids redoing the projection/convolution work +// independently for every (row, column) pair. void aggame::payoffMatrix(cmatrix &dest, const cvector &s, double fuzz) const { - // compute jacobian - // s: mixed strat - -#ifdef AGGDEBUG - cout << "calling payoffMatrix with stratety s=" << endl << s << endl; -#endif - double fuzzcount; - int rown, coln, rowi, coli, act1, act2, currNode, numNei; - static std::vector::iterator p; - static std::vector tasks, spares, nontasks; - tasks.reserve(aggPtr->numPlayers); - spares.reserve(aggPtr->numPlayers); - nontasks.reserve(aggPtr->numPlayers); - aggPtr->cache.reset(); - - // do projection - for (int Node = 0; Node < aggPtr->numActionNodes; Node++) { - aggPtr->doProjection(Node, s.values()); - } + const std::vector sp(s.values(), s.values() + s.getm()); // deal with the diagonal - for (rown = 0; rown < aggPtr->numPlayers; ++rown) { - fuzzcount = fuzz; - for (rowi = firstAction(rown); rowi < lastAction(rown); rowi++) { - for (coli = firstAction(rown); coli < lastAction(rown); coli++) { + for (int rown = 0; rown < aggPtr->getNumPlayers(); ++rown) { + double fuzzcount = fuzz; + for (int rowi = firstAction(rown); rowi < lastAction(rown); rowi++) { + for (int coli = firstAction(rown); coli < lastAction(rown); coli++) { dest(rowi, coli) = fuzzcount; fuzzcount += fuzz; } } } - for (rown = 0; rown < aggPtr->numPlayers; ++rown) { // rown: the row player - for (act1 = 0; act1 < aggPtr->actions[rown]; act1++) { // act1: player rown's action - - currNode = aggPtr->actionSets[rown][act1]; - numNei = aggPtr->neighbors[currNode].size(); -#ifdef AGGDEBUG - cout << "for player " << rown << ", action " << act1 << ", action node " << currNode << endl; - cout << "cache is: " << endl << aggPtr->cache << endl; -#endif - tasks.clear(); // for these col players, we need to compute the distribution induced by their - // complements. input of the bisection alg - spares.clear(); // these col players have only one projected action - nontasks.clear(); // complement of tasks. includes spares, and completely cached col players - std::vector key(numNei + 3); - key[numNei] = rown; - key[numNei + 1] = act1; - - // first, populate tasks, spares and nontasks - for (coln = 0; coln < aggPtr->numPlayers; ++coln) { - if (rown != coln) { // coln: col player - - key[numNei + 2] = coln; - bool allCached = true; - for (act2 = 0; act2 < aggPtr->actions[coln]; ++act2) { - - copy(aggPtr->projection[currNode][coln][act2].begin(), - aggPtr->projection[currNode][coln][act2].end(), key.begin()); -#ifdef AGGDEBUG - cout << "for player2=" << coln << " act2=" << act2 << endl; - cout << "checking cache for: ["; - copy(key.begin(), key.end(), ostream_iterator(cout, " ")); - cout << "]\n"; -#endif - auto r = aggPtr->cache.findExact(key); - if (r != aggPtr->cache.end()) { - dest(act1 + firstAction(rown), act2 + firstAction(coln)) = r->second; - } - else { - allCached = false; - } - } - if (allCached) { // if all coln's actions are already cached: - nontasks.push_back(coln); - } - else { - if (aggPtr->fullProjectedStrat[currNode][coln].size() == - 1) { // if coln has only one projected action - spares.push_back(coln); - nontasks.push_back(coln); - } - else { - tasks.push_back(coln); - } - } - } - } -#ifdef AGGDEBUG - cout << "spares are: ["; - copy(spares.begin(), spares.end(), ostream_iterator(cout, " ")); - cout << "]\ntasks are :["; - copy(tasks.begin(), tasks.end(), ostream_iterator(cout, " ")); - cout << "]\n"; -#endif - - // compute partial prob distributions - if (tasks.empty() && spares.empty()) { - continue; - } // nothing to be done for this row - - if (aggPtr->isPure[currNode] || tasks.empty()) { - computePartialP_PureNode(rown, act1, tasks); - } - else { // do bisection - computePartialP_bisect(rown, act1, tasks.begin(), tasks.end(), aggPtr->m_state.Pr[rown]); -#ifdef AGGDEBUG - cout << "after calling computePartialP_bisect:" << endl; - for (int tt = 0; tt < tasks.size(); tt++) { - cout << "for player " << tasks[tt] << endl; - cout << aggPtr->m_state.Pr[tasks[tt]] << endl; - } -#endif - // now apply rown's action (act1), and the strategies of - // players in nontasks - aggPtr->m_state.Pr[rown].reset(); - aggPtr->m_state.Pr[rown].insert(make_pair(aggPtr->projection[currNode][rown][act1], 1.0)); - for (p = nontasks.begin(); p != nontasks.end(); ++p) { - aggPtr->m_state.Pr[rown].multiply(aggPtr->m_state.projectedStrat[currNode][*p], numNei, - aggPtr->projFunctions[currNode]); + // off-diagonal entries + std::vector row(aggPtr->getNumActions()); + for (int rown = 0; rown < aggPtr->getNumPlayers(); ++rown) { + for (int act1 = 0; act1 < aggPtr->getNumActions(rown); act1++) { + aggPtr->getPayoffJacobianRow(rown, act1, sp, row); + for (int coln = 0; coln < aggPtr->getNumPlayers(); ++coln) { + if (coln == rown) { + continue; } -#ifdef AGGDEBUG - cout << "the polynomial product of strats of player " << rown - << " and players in the vector nontasks is:" << endl; - cout << aggPtr->m_state.Pr[rown] << endl; -#endif - if (tasks.size() == 1) { - aggPtr->m_state.Pr[tasks[0]] = aggPtr->m_state.Pr[rown]; + for (int act2 = 0; act2 < aggPtr->getNumActions(coln); act2++) { + dest(act1 + firstAction(rown), act2 + firstAction(coln)) = row[act2 + firstAction(coln)]; } - else { - for (p = tasks.begin(); p != tasks.end(); ++p) { - if (aggPtr->m_state.Pr[*p].empty()) { - std::cerr << "AGG::payoffMatrix() ERROR for rown=" << rown << " act1=" << act1 - << " *p=" << *p << ": the distribution should not be empty!" << std::endl; -#ifdef AGGDEBUG - std::cerr << "strategy profile is: " << s << std::endl; -#endif - } - aggPtr->m_state.Pr[*p].multiply(aggPtr->m_state.Pr[rown], numNei, - aggPtr->projFunctions[currNode]); - } // end for(p=tasks.begin... - } - - // if spares not empty, we need to compute nondisturbed payoffs - // which requires the distrib induced by everyone (except rown). - // we store this distrib in Pr[rown][act1][rown] - if (!spares.empty()) { - // assert(tasks.size()>0); - aggPtr->m_state.Pr[rown].reset(); - aggPtr->m_state.Pr[rown].multiply(aggPtr->m_state.Pr[tasks[0]], - aggPtr->m_state.projectedStrat[currNode][tasks[0]], - numNei, aggPtr->projFunctions[currNode]); - } - } // end else -#ifdef AGGDEBUG - cout << "after computing parital distributions, the distributions are" << endl; - for (int tt = 0; tt < numPlayers; tt++) { - cout << "for player " << tt << endl; - cout << aggPtr->m_state.Pr[tt]; - cout << endl; } -#endif - // compute entries - double undisturbedPayoff; - bool hasUndisturbed = false; - - if (!spares.empty()) { // for players in spares, we compute one undisturbed payoff - computeUndisturbedPayoff(undisturbedPayoff, hasUndisturbed, rown, act1, rown); - for (p = spares.begin(); p != spares.end(); ++p) { - for (act2 = 0; act2 < aggPtr->actions[*p]; act2++) { - savePayoff(dest, rown, act1, *p, act2, undisturbedPayoff, aggPtr->cache); - } - } - } - for (p = tasks.begin(); p != tasks.end(); ++p) { - for (act2 = 0; act2 < aggPtr->actions[*p]; act2++) { // act2: col action - - if (aggPtr->m_state.projectedStrat[currNode][*p].size() == 1 && - aggPtr->m_state.projectedStrat[currNode][*p].begin()->first == - aggPtr->projection[currNode][*p][act2]) { - computeUndisturbedPayoff(undisturbedPayoff, hasUndisturbed, rown, act1, *p); - savePayoff(dest, rown, act1, *p, act2, undisturbedPayoff, aggPtr->cache); - } - computePayoff(dest, rown, act1, *p, act2, aggPtr->cache); - } // end for(act2 - } // end for(p - } // end for(act1 - } // end for(rown -} - -void aggame::computeUndisturbedPayoff(double &undisturbedPayoff, bool &has, int player1, int act1, - int player2) const -{ - if (has) { - return; - } - const int Node = aggPtr->actionSets[player1][act1]; - const int numNei = aggPtr->neighbors[Node].size(); - if (player2 == player1) { - undisturbedPayoff = aggPtr->m_state.Pr[player2].inner_prod(aggPtr->payoffs[Node]); - } - else { - // assert(aggPtr->m_state.projectedStrat[Node][player2].size()==1); - undisturbedPayoff = aggPtr->m_state.Pr[player2].inner_prod( - aggPtr->m_state.projectedStrat[Node][player2].begin()->first, numNei, - aggPtr->projFunctions[Node], aggPtr->payoffs[Node]); - } - has = true; -} - -void aggame::savePayoff(cmatrix &dest, int player1, int act1, int player2, int act2, double result, - agg::trie_map &cache, bool partial) const -{ - - const int Node = aggPtr->actionSets[player1][act1]; - const int numNei = aggPtr->neighbors[Node].size(); - - if (!partial) { - std::pair, double> pair1(aggPtr->projection[Node][player2][act2], result); - pair1.first.reserve(numNei + 3); - pair1.first.push_back(player1); - pair1.first.push_back(act1); - pair1.first.push_back(player2); - cache.insert(pair1); - } - - if (aggPtr->node2Action[Node][player2] != -1 && - aggPtr->fullProjectedStrat[Node][player1].count(aggPtr->projection[Node][player2][act2])) { - std::pair, double> pair2(aggPtr->projection[Node][player2][act2], result); - pair2.first.reserve(numNei + 3); - pair2.first.push_back(player2); - pair2.first.push_back(aggPtr->node2Action[Node][player2]); - pair2.first.push_back(player1); - cache.insert(pair2); - } - dest(act1 + firstAction(player1), act2 + firstAction(player2)) = result; -} - -void aggame::computePayoff(cmatrix &dest, int player1, int act1, int player2, int act2, - agg::trie_map &cache) const -{ - const int Node = aggPtr->actionSets[player1][act1]; - const int numNei = aggPtr->neighbors[Node].size(); - - std::pair, double> insPair(aggPtr->projection[Node][player2][act2], 0); - insPair.first.reserve(numNei + 3); - insPair.first.push_back(player1); - insPair.first.push_back(act1); - insPair.first.push_back(player2); - const std::pair::iterator, bool> r = cache.insert(insPair); - if (!r.second) { - dest(act1 + firstAction(player1), act2 + firstAction(player2)) = r.first->second; - } - else { - r.first->second = - aggPtr->m_state.Pr[player2].inner_prod(aggPtr->projection[Node][player2][act2], numNei, - aggPtr->projFunctions[Node], aggPtr->payoffs[Node]); - savePayoff(dest, player1, act1, player2, act2, r.first->second, cache, r.second); + } } } @@ -409,7 +69,7 @@ void aggame::KSymPayoffMatrix(cmatrix &dest, const cvector &s, double fuzz) cons for (int rowcls = 0; rowcls < getNumPlayerClasses(); ++rowcls) { for (int rowa = 0; rowa < getNumKSymActions(rowcls); ++rowa) { for (int colcls = 0; colcls < getNumPlayerClasses(); ++colcls) { - int multiplier = aggPtr->playerClasses[colcls].size(); + int multiplier = aggPtr->getPlayerClass(colcls).size(); if (rowcls == colcls) { multiplier -= 1; } diff --git a/src/solvers/gtracer/aggame.h b/src/solvers/gtracer/aggame.h index d1b4b8346d..1516b2b6cd 100644 --- a/src/solvers/gtracer/aggame.h +++ b/src/solvers/gtracer/aggame.h @@ -35,7 +35,7 @@ namespace Gambit::gametracer { class aggame : public gnmgame { public: explicit aggame(const Gambit::GameAGGRep &g) - : gnmgame(g.GetUnderlyingAGG()->actions), aggPtr(g.GetUnderlyingAGG()) + : gnmgame(g.GetUnderlyingAGG()->getActionCounts()), aggPtr(g.GetUnderlyingAGG()) { } @@ -51,7 +51,7 @@ class aggame : public gnmgame { double getKSymMixedPayoff(int cls, cvector &s) { std::vector sp(s.values(), s.values() + s.getm()); - return (double)aggPtr->getKSymMixedPayoff(cls, sp); + return aggPtr->getKSymMixedPayoff(cls, sp); } void payoffMatrix(cmatrix &dest, const cvector &s, double fuzz) const override; @@ -89,18 +89,6 @@ class aggame : public gnmgame { private: std::shared_ptr aggPtr; - - // helper functions for computing jacobian - void computePartialP_PureNode(int player, int act, std::vector &tasks) const; - void computePartialP_bisect(int player, int act, std::vector::iterator f, - std::vector::iterator l, - Gambit::agg::ConfigDistribution &temp) const; - void computePayoff(cmatrix &dest, int player1, int act1, int player2, int act2, - Gambit::agg::trie_map &cache) const; - void savePayoff(cmatrix &dest, int player1, int act1, int player2, int act2, double result, - Gambit::agg::trie_map &cache, bool partial = false) const; - void computeUndisturbedPayoff(double &undisturbedPayoff, bool &has, int player1, int act1, - int player2) const; }; } // end namespace Gambit::gametracer diff --git a/src/tools/gt/nfgipa.cc b/src/tools/gt/nfgipa.cc index 6848b25d2c..207eceefa1 100644 --- a/src/tools/gt/nfgipa.cc +++ b/src/tools/gt/nfgipa.cc @@ -29,7 +29,6 @@ using namespace Gambit; using namespace Gambit::Nash; -using namespace Gambit::gametracer; extern std::vector> ReadStrategyPerturbations(const Game &p_game, std::istream &p_stream); diff --git a/tests/test_agg.py b/tests/test_agg.py index c417a6f3af..c22f656cd4 100644 --- a/tests/test_agg.py +++ b/tests/test_agg.py @@ -45,8 +45,10 @@ def test_agg_bagg_to_nfg(game_path): def test_agg_fraction_and_long_decimal_payoffs_parsed_exactly(): - """Payoffs written as a fraction ("1/3") or a decimal too precise for a double to round-trip - are parsed into exact Rationals, not silently rounded through a bare double read. + """Payoffs written in an .agg file as a fraction ("1/3") or a decimal needing more than + double's ~15-17 significant digits of precision to round-trip are parsed into exact + Rationals (via Number, agg::AGG::exactPayoffs), not silently rounded through a bare + `istream >> double` read. 2x2_fraction_payoffs.agg is 2x2.agg with two of its four payoffs replaced: 35 -> 1/3 and 95 -> 0.123456789012345. @@ -58,22 +60,33 @@ def test_agg_fraction_and_long_decimal_payoffs_parsed_exactly(): assert player.max_payoff == gbt.Rational(1, 3) assert player.min_payoff == -10 - # pure-strategy payoff lookup also reports the exact value, not a rounded approximation - profile = game.mixed_strategy_profile(rational=True) + # pure-strategy payoff lookup (AGGPureStrategyProfileRep::GetPayoff) also reports the + # exact value, not a double-rounded approximation -- for both the fraction and the long + # decimal payoff p0, p1 = game.players s0 = list(p0.strategies) s1 = list(p1.strategies) + + profile = game.mixed_strategy_profile(rational=True) profile[s0[0]] = gbt.Rational(1) profile[s0[1]] = gbt.Rational(0) profile[s1[0]] = gbt.Rational(0) profile[s1[1]] = gbt.Rational(1) assert profile.payoff(p0) == gbt.Rational(1, 3) + profile = game.mixed_strategy_profile(rational=True) + profile[s0[0]] = gbt.Rational(0) + profile[s0[1]] = gbt.Rational(1) + profile[s1[0]] = gbt.Rational(0) + profile[s1[1]] = gbt.Rational(1) + assert profile.payoff(p0) == gbt.Rational("0.123456789012345") + assert profile.payoff(p1) == gbt.Rational("0.123456789012345") + def test_bagg_fraction_type_distribution_parsed_exactly(): - """BAGG type-distribution probabilities written as a fraction ("1/3") are parsed exactly, - and the exact payoff computation over them agrees with the double engine on the same - profile. + """BAGG type-distribution probabilities written as a fraction ("1/3") are parsed exactly + (agg::BAGG::exactIndepTypeDist) and the exact weighted-sum payoff computation over them + (agg::BAGG::getMixedPayoff) agrees with the double engine on the same profile. 2x2_fraction_types.bagg: player 1 has two types (weights 1/3, 2/3), player 2 has one. """ @@ -92,10 +105,13 @@ def test_bagg_fraction_type_distribution_parsed_exactly(): @pytest.mark.parametrize("game_path", ["2x2.agg", "2x2.bagg"]) def test_agg_bagg_mixed_strategy_profile_rational_exact_payoff(game_path): - """AGG/BAGG mixed-strategy payoffs support exact (rational) computation, not just double. + """AGG/BAGG mixed-strategy payoffs support exact (rational) computation: the convolution + algorithm (agg::AGG::getMixedPayoff et al.) is generic in its numeric type, so it also runs + with Rational arithmetic throughout (agg::AGG::getMixedPayoff), not just double. 2x2.agg has a known mixed equilibrium at (10/11, 1/11) for both players, with payoff - exactly -5/11 to each -- computed here directly (not via a solver), confirming zero regret. + exactly -5/11 to each -- computed here directly (not via a solver), confirming the exact + engine reproduces it with zero regret. """ game = games.read_from_file(game_path) profile = game.mixed_strategy_profile(rational=True)