diff --git a/src/core/segment.h b/src/core/segment.h index 84f8f7abc..00edd74b6 100644 --- a/src/core/segment.h +++ b/src/core/segment.h @@ -99,6 +99,7 @@ template class Segmented { void SetFlattened(const Storage &v) { m_values = v; } const Storage &GetFlattened() const { return m_values; } + Storage &GetFlattened() { return m_values; } }; template using SegmentedArray = Segmented>; diff --git a/src/games/game.cc b/src/games/game.cc index c33e1840e..e3298dd1d 100644 --- a/src/games/game.cc +++ b/src/games/game.cc @@ -249,14 +249,16 @@ void GameRep::WriteNfgFile(std::ostream &p_file) const template MixedStrategyProfileRep::MixedStrategyProfileRep(const StrategySupportProfile &p_support) - : m_probs(p_support.MixedProfileLength()), m_support(p_support), + : m_probs(p_support.GetShape()), m_offsets(p_support.GetShape()), m_support(p_support), m_gameversion(p_support.GetGame()->GetVersion()) { int index = 1; for (auto player : p_support.GetGame()->GetPlayers()) { for (auto strategy : player->GetStrategies()) { if (p_support.Contains(strategy)) { - m_profileIndex[strategy] = index++; + m_offsets.GetFlattened()[index] = StrategyOffset(strategy); + m_profileIndex[strategy] = index; + index++; } else { m_profileIndex[strategy] = -1; @@ -367,8 +369,20 @@ template void MixedStrategyProfile::ComputePayoffs() const Cache newCache; for (const auto &player : m_rep->GetSupport().GetPlayers()) { newCache.m_payoffs[player] = GetPayoff(player); - for (const auto &strategy : m_rep->GetSupport().GetStrategies(player)) { - newCache.m_strategyValues[player][strategy] = GetPayoff(strategy); + const auto &strategies = m_rep->GetSupport().GetStrategies(player); + Vector values(strategies.size()); + if (m_rep->GetPayoffDerivs(player->GetNumber(), values)) { + auto value_it = values.begin(); + for (const auto &strategy : strategies) { + newCache.m_strategyValues[player][strategy] = *value_it; + ++value_it; + ; + } + } + else { + for (const auto &strategy : m_rep->GetSupport().GetStrategies(player)) { + newCache.m_strategyValues[player][strategy] = GetPayoff(strategy); + } } } newCache.m_valid = true; diff --git a/src/games/game.h b/src/games/game.h index 420d072bd..63ea083ac 100644 --- a/src/games/game.h +++ b/src/games/game.h @@ -624,6 +624,7 @@ class GameRep : public std::enable_shared_from_this { friend class StrategySupportProfile; template friend class MixedBehaviorProfile; template friend class MixedStrategyProfile; + template friend class MixedStrategyProfileRep; template friend class TableMixedStrategyProfileRep; protected: diff --git a/src/games/gameagg.cc b/src/games/gameagg.cc index ccdda7d17..6341459c4 100644 --- a/src/games/gameagg.cc +++ b/src/games/gameagg.cc @@ -102,7 +102,7 @@ template T AGGMixedStrategyProfileRep::GetPayoff(int pl) const const GameStrategy strategy = this->m_support.GetGame()->GetPlayer(i + 1)->GetStrategy(j + 1); const int ind = this->m_profileIndex.at(strategy); - s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs[ind]; + s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs.GetFlattened()[ind]; } } return (T)g.aggPtr->getMixedPayoff(pl - 1, s); @@ -125,7 +125,7 @@ T AGGMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &ps) const GameStrategy strategy = this->m_support.GetGame()->GetPlayer(i + 1)->GetStrategy(j + 1); const int &ind = this->m_profileIndex.at(strategy); - s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs[ind]; + s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs.GetFlattened()[ind]; } } } @@ -162,7 +162,7 @@ T AGGMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &ps1, const GameStrategy strategy = this->m_support.GetGame()->GetPlayer(i + 1)->GetStrategy(j + 1); const int ind = this->m_profileIndex.at(strategy); - s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs[ind]; + s[g.aggPtr->firstAction(i) + j] = (ind == -1) ? (T)0 : this->m_probs.GetFlattened()[ind]; } } } diff --git a/src/games/gamebagg.cc b/src/games/gamebagg.cc index 30ee5848f..fcc669357 100644 --- a/src/games/gamebagg.cc +++ b/src/games/gamebagg.cc @@ -113,7 +113,7 @@ template T BAGGMixedStrategyProfileRep::GetPayoff(int pl) const ->GetPlayer(g.baggPtr->typeOffset[i] + tp + 1) ->GetStrategy(j + 1); const int ind = this->m_profileIndex.at(strategy); - s.at(offs) = (ind == -1) ? (T)0 : this->m_probs[ind]; + s.at(offs) = (ind == -1) ? (T)0 : this->m_probs.GetFlattened()[ind]; } } } @@ -144,7 +144,8 @@ T BAGGMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &ps) ->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) ? Rational(0) : this->m_probs[ind]; + s.at(g.baggPtr->firstAction(i, tp) + j) = + (ind == -1) ? Rational(0) : this->m_probs.GetFlattened()[ind]; } } } @@ -191,7 +192,7 @@ T BAGGMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &ps1 ->GetStrategy(j + 1); const int ind = this->m_profileIndex.at(strategy); s.at(g.baggPtr->firstAction(i, tp) + j) = - static_cast((ind == -1) ? T(0) : this->m_probs[ind]); + static_cast((ind == -1) ? T(0) : this->m_probs.GetFlattened()[ind]); } } } diff --git a/src/games/gametable.cc b/src/games/gametable.cc index 2a2d11952..bb810d950 100644 --- a/src/games/gametable.cc +++ b/src/games/gametable.cc @@ -102,30 +102,179 @@ PureStrategyProfile GameTableRep::NewPureStrategyProfile() const // TableMixedStrategyProfileRep //======================================================================== -template class TableMixedStrategyProfileRep : public MixedStrategyProfileRep { -private: - /// @name Private recursive payoff functions - //@{ - /// Recursive computation of payoff to player pl - T GetPayoff(int pl, int index, int i) const; - /// Recursive computation of payoff derivative - void GetPayoffDeriv(int pl, int const_pl, int cur_pl, long index, const T &prob, T &value) const; - /// Recursive computation of payoff second derivative - void GetPayoffDeriv(int pl, int const_pl1, int const_pl2, int cur_pl, long index, const T &prob, - T &value) const; - //@} - - const CartesianProductSpace &m_pureStrategies; - - long StrategyOffset(const GameStrategy &s) const +template class ProductDistribution { +public: + using index_type = long; + using prob_type = T; + using value_type = std::pair; + + class iterator { + public: + using iterator_category = std::forward_iterator_tag; + using difference_type = std::ptrdiff_t; + using value_type = std::pair; + using reference = value_type; + using pointer = void; + + iterator() = default; + + iterator(const SegmentedVector &probs, const SegmentedArray &offsets, size_t skip1, + size_t skip2, bool end) + : m_probs(&probs), m_offsets(&offsets), m_done(end) + { + if (m_done) { + return; + } + + const size_t P = m_probs->GetShape().size(); + + // Build active dimension list + for (size_t p = 1; p <= P; ++p) { + if (p != skip1 && p != skip2) { + m_dims.push_back(p); + } + } + + m_K = m_dims.size(); + + m_digit.assign(m_K, 0); + m_radix.assign(m_K, 0); + m_cum_prob.assign(m_K + 1, T{}); + + // initialise radices + for (size_t j = 0; j < m_K; ++j) { + const size_t p = m_dims[j]; + m_radix[j] = m_probs->segment(p).size(); + if (m_radix[j] == 0) { + m_done = true; + return; + } + } + + // initial recompute + recompute_from(0); + advance_to_next_nonzero(); + } + + reference operator*() const { return {m_index, m_cum_prob[m_K]}; } + + iterator &operator++() + { + if (m_done) { + return *this; + } + + // increment odometer + size_t j = 0; + for (; j < m_K; ++j) { + if (++m_digit[j] < m_radix[j]) { + break; + } + m_digit[j] = 0; + } + + if (j == m_K) { + m_done = true; + return *this; + } + + recompute_from(j); + advance_to_next_nonzero(); + return *this; + } + + iterator operator++(int) + { + iterator tmp = *this; + ++(*this); + return tmp; + } + + bool operator==(const iterator &other) const { return m_done == other.m_done; } + + bool operator!=(const iterator &other) const { return !(*this == other); } + + private: + void recompute_from(size_t j0) + { + m_index = 0; + m_cum_prob[0] = T{1}; + + for (size_t j = 0; j < m_K; ++j) { + const size_t p = m_dims[j]; + const size_t d = m_digit[j] + 1; + + const T pi = m_probs->segment(p)[d]; + m_cum_prob[j + 1] = m_cum_prob[j] * pi; + m_index += m_offsets->segment(p)[d]; + } + } + + void advance_to_next_nonzero() + { + while (!m_done && m_cum_prob[m_K] == T{0}) { + size_t j = 0; + for (; j < m_K; ++j) { + if (++m_digit[j] < m_radix[j]) { + break; + } + m_digit[j] = 0; + } + + if (j == m_K) { + m_done = true; + return; + } + + recompute_from(j); + } + } + + const SegmentedVector *m_probs{nullptr}; + const SegmentedArray *m_offsets{nullptr}; + + std::vector m_dims; // active player numbers + size_t m_K{0}; + bool m_done{true}; + + std::vector m_digit; + std::vector m_radix; + std::vector m_cum_prob; + index_type m_index{0}; + }; + + ProductDistribution(const SegmentedVector &probs, const SegmentedArray &offsets) + : m_probs(probs), m_offsets(offsets), m_skip1(0), m_skip2(0) { - return (s->GetNumber() - 1) * m_pureStrategies.m_strides[s->m_player->GetNumber() - 1]; } + ProductDistribution(const SegmentedVector &probs, const SegmentedArray &offsets, + size_t skip1) + : m_probs(probs), m_offsets(offsets), m_skip1(skip1), m_skip2(0) + { + } + + ProductDistribution(const SegmentedVector &probs, const SegmentedArray &offsets, + size_t skip1, size_t skip2) + : m_probs(probs), m_offsets(offsets), m_skip1(skip1), m_skip2(skip2) + { + } + + iterator begin() const { return iterator(m_probs, m_offsets, m_skip1, m_skip2, false); } + + iterator end() const { return iterator(m_probs, m_offsets, m_skip1, m_skip2, true); } + +private: + const SegmentedVector &m_probs; + const SegmentedArray &m_offsets; + size_t m_skip1; + size_t m_skip2; +}; + +template class TableMixedStrategyProfileRep : public MixedStrategyProfileRep { public: explicit TableMixedStrategyProfileRep(const StrategySupportProfile &p_support) - : MixedStrategyProfileRep(p_support), - m_pureStrategies(p_support.GetGame()->m_pureStrategies) + : MixedStrategyProfileRep(p_support) { } ~TableMixedStrategyProfileRep() override = default; @@ -133,6 +282,7 @@ template class TableMixedStrategyProfileRep : public MixedStrategyProf std::unique_ptr> Copy() const override; T GetPayoff(int pl) const override; T GetPayoffDeriv(int pl, const GameStrategy &) const override; + bool GetPayoffDerivs(int pl, Vector &p_derivs) const override; T GetPayoffDeriv(int pl, const GameStrategy &, const GameStrategy &) const override; }; @@ -142,105 +292,76 @@ std::unique_ptr> TableMixedStrategyProfileRep::Cop return std::make_unique(*this); } -template -T TableMixedStrategyProfileRep::GetPayoff(int pl, int index, int current) const -{ - if (current > static_cast(this->GetSupport().GetGame()->NumPlayers())) { - const Game game = this->GetSupport().GetGame(); - auto &g = dynamic_cast(*game); - if (const auto outcome = g.m_results[index]) { - return outcome->GetPayoff(this->GetSupport().GetGame()->GetPlayer(pl)); - } - return static_cast(0); - } - - T sum = static_cast(0); - for (auto s : - this->GetSupport().GetStrategies(this->GetSupport().GetGame()->GetPlayer(current))) { - if ((*this)[s] != T(0)) { - sum += ((*this)[s] * GetPayoff(pl, index + StrategyOffset(s), current + 1)); - } - } - return sum; -} - template T TableMixedStrategyProfileRep::GetPayoff(int pl) const { - return GetPayoff(pl, 0, 1); -} - -template -void TableMixedStrategyProfileRep::GetPayoffDeriv(int pl, int const_pl, int cur_pl, long index, - const T &prob, T &value) const -{ - if (cur_pl == const_pl) { - cur_pl++; - } - if (cur_pl > static_cast(this->GetSupport().GetGame()->NumPlayers())) { - const Game game = this->GetSupport().GetGame(); - auto &g = dynamic_cast(*game); + const auto game = this->GetSupport().GetGame(); + auto &g = dynamic_cast(*game); + const auto player = game->GetPlayer(pl); + T value{0}; + for (auto [index, prob] : ProductDistribution(this->m_probs, this->m_offsets)) { if (const auto outcome = g.m_results[index]) { - value += prob * outcome->GetPayoff(this->GetSupport().GetGame()->GetPlayer(pl)); - } - } - else { - for (auto s : - this->GetSupport().GetStrategies(this->GetSupport().GetGame()->GetPlayer(cur_pl))) { - if ((*this)[s] > T(0)) { - GetPayoffDeriv(pl, const_pl, cur_pl + 1, index + StrategyOffset(s), prob * (*this)[s], - value); - } + value += prob * outcome->template GetPayoff(player); } } + return value; } template T TableMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &strategy) const { - T value = T(0); - GetPayoffDeriv(pl, strategy->GetPlayer()->GetNumber(), 1, StrategyOffset(strategy), T(1), value); + const auto game = this->GetSupport().GetGame(); + auto &g = dynamic_cast(*game); + auto base_index = this->StrategyOffset(strategy); + const auto player = game->GetPlayer(pl); + T value{0}; + for (auto [index, prob] : ProductDistribution(this->m_probs, this->m_offsets, + strategy->GetPlayer()->GetNumber())) { + if (const auto outcome = g.m_results[base_index + index]) { + value += prob * outcome->template GetPayoff(player); + } + } return value; } template -void TableMixedStrategyProfileRep::GetPayoffDeriv(int pl, int const_pl1, int const_pl2, - int cur_pl, long index, const T &prob, - T &value) const +bool TableMixedStrategyProfileRep::GetPayoffDerivs(int pl, Vector &p_derivs) const { - while (cur_pl == const_pl1 || cur_pl == const_pl2) { - cur_pl++; - } - if (cur_pl > static_cast(this->GetSupport().GetGame()->NumPlayers())) { - const Game game = this->GetSupport().GetGame(); - auto &g = dynamic_cast(*game); - if (const auto outcome = g.m_results[index]) { - value += prob * outcome->GetPayoff(this->GetSupport().GetGame()->GetPlayer(pl)); - } - } - else { - for (auto s : - this->GetSupport().GetStrategies(this->GetSupport().GetGame()->GetPlayer(cur_pl))) { - if ((*this)[s] > static_cast(0)) { - GetPayoffDeriv(pl, const_pl1, const_pl2, cur_pl + 1, index + StrategyOffset(s), - prob * (*this)[s], value); + const auto game = this->GetSupport().GetGame(); + auto &g = dynamic_cast(*game); + const auto player = game->GetPlayer(pl); + p_derivs = T{0}; + auto segment = this->m_offsets.segment(pl); + for (auto [index, prob] : ProductDistribution(this->m_probs, this->m_offsets, pl)) { + auto deriv_it = p_derivs.begin(); + for (const auto base_index : segment) { + if (const auto outcome = g.m_results[base_index + index]) { + *deriv_it += prob * outcome->template GetPayoff(player); + ++deriv_it; } } } + return true; } template T TableMixedStrategyProfileRep::GetPayoffDeriv(int pl, const GameStrategy &strategy1, const GameStrategy &strategy2) const { - const auto player1 = strategy1->GetPlayer().get(); - const auto player2 = strategy2->GetPlayer().get(); - if (player1 == player2) { - return T(0); + if (strategy1->GetPlayer() == strategy2->GetPlayer()) { + return T{0}; + } + const auto game = this->GetSupport().GetGame(); + auto &g = dynamic_cast(*game); + auto base_index = this->StrategyOffset(strategy1) + this->StrategyOffset(strategy2); + const auto player = game->GetPlayer(pl); + T value{0}; + for (auto [index, prob] : + ProductDistribution(this->m_probs, this->m_offsets, strategy1->GetPlayer()->GetNumber(), + strategy2->GetPlayer()->GetNumber())) { + if (const auto outcome = g.m_results[base_index + index]) { + value += prob * outcome->template GetPayoff(player); + } } - - T value = T(0); - GetPayoffDeriv(pl, player1->GetNumber(), player2->GetNumber(), 1, - StrategyOffset(strategy1) + StrategyOffset(strategy2), T(1), value); return value; } diff --git a/src/games/stratmixed.h b/src/games/stratmixed.h index 3a5da5b59..bce9f0e16 100644 --- a/src/games/stratmixed.h +++ b/src/games/stratmixed.h @@ -30,12 +30,21 @@ namespace Gambit { template class MixedStrategyProfileRep { protected: - Vector m_probs; + SegmentedVector m_probs; + SegmentedArray m_offsets; StrategySupportProfile m_support; /// The index into the strategy profile for a strategy (-1 if not in support) std::map m_profileIndex; unsigned int m_gameversion; + long StrategyOffset(const GameStrategy &s) const + { + const auto &space = this->GetSupport().GetGame()->m_pureStrategies; + const auto &player = s->GetPlayer(); + const size_t i = player->GetNumber() - 1; + return (s->GetNumber() - 1) * space.m_strides[i]; + } + public: explicit MixedStrategyProfileRep(const StrategySupportProfile &); virtual ~MixedStrategyProfileRep() = default; @@ -46,22 +55,22 @@ template class MixedStrategyProfileRep { void SetCentroid(); std::unique_ptr Normalize() const; - const T &operator[](int i) const { return m_probs[i]; } + const T &operator[](int i) const { return m_probs.GetFlattened()[i]; } T &operator[](int i) { OnProfileChanged(); - return m_probs[i]; + return m_probs.GetFlattened()[i]; } /// Returns the probability the strategy is played const T &operator[](const GameStrategy &p_strategy) const { - return m_probs[m_profileIndex.at(p_strategy)]; + return m_probs.GetFlattened()[m_profileIndex.at(p_strategy)]; } /// Returns the probability the strategy is played T &operator[](const GameStrategy &p_strategy) { OnProfileChanged(); - return m_probs[m_profileIndex.at(p_strategy)]; + return m_probs.GetFlattened()[m_profileIndex.at(p_strategy)]; } /// Set the strategy of the corresponding player to a pure strategy void SetStrategy(const GameStrategy &p_strategy) @@ -72,15 +81,15 @@ template class MixedStrategyProfileRep { (*this)[p_strategy] = static_cast(1); OnProfileChanged(); } - const Vector &GetProbVector() const { return m_probs; } + const Vector &GetProbVector() const { return m_probs.GetFlattened(); } void SetProbVector(const Vector &p_vector) { - m_probs = p_vector; + m_probs.SetFlattened(p_vector); OnProfileChanged(); } void SetProbConstant(const T &c) { - m_probs = c; + m_probs.GetFlattened() = c; OnProfileChanged(); } unsigned int GetProfileIndex(const GameStrategy &p_strategy) const @@ -89,6 +98,7 @@ template class MixedStrategyProfileRep { } virtual T GetPayoff(int pl) const = 0; virtual T GetPayoffDeriv(int pl, const GameStrategy &) const = 0; + virtual bool GetPayoffDerivs(int pl, Vector &p_derivs) const { return false; } virtual T GetPayoffDeriv(int pl, const GameStrategy &, const GameStrategy &) const = 0; T GetPayoff(const GamePlayer &p_player) const { return GetPayoff(p_player->GetNumber()); } diff --git a/src/games/stratspt.h b/src/games/stratspt.h index 0ceef7ac3..22990aeb5 100644 --- a/src/games/stratspt.h +++ b/src/games/stratspt.h @@ -98,6 +98,14 @@ class StrategySupportProfile { /// Returns the total number of strategies in the support. int MixedProfileLength() const; + Array GetShape() const + { + Array shape(NumPlayers()); + std::transform(m_strategyDigits.m_allowedDigits.begin(), + m_strategyDigits.m_allowedDigits.end(), shape.begin(), + [](const auto &c) { return c.size(); }); + return shape; + } template MixedStrategyProfile NewMixedStrategyProfile() const; /// Returns the number of players in the game