Squashed commit of the following:

commit a25d065a0618b8be6c78c6303a7318b5557a78c5
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Thu Feb 27 22:06:39 2025 +0000

    refactor NormalGames.jl dependency into SearchNE.jl

commit 1c24319e18b41ae5e5b20a1b5f5a56d6964fb06b
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Thu Feb 27 21:57:11 2025 +0000

    restrict polymatrix approach to players with bilateral payoffs

commit e4aa3c7ac1e6f7d78fd5a234159f871eb70b02da
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Thu Feb 27 21:45:55 2025 +0000

    simplify payoff functions

commit 44ba2d5c8d972aace19516a6fc97e59951dfaa09
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Thu Feb 27 21:45:37 2025 +0000

    fix expected total payoff for a player

commit 1ea9395277f35d9eb7f320cf616f59f45edc28ca
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Thu Feb 27 21:38:22 2025 +0000

    define type for bilateral payoff functions and add documentation

commit fc5171a
Author: Bruno Machado Pacheco <mpacheco.bruno@gmail.com>
Date:   Mon Feb 24 23:26:04 2025 +0000

    create payoff directory; correct abstractness in function arguments

Author: brunompacheco

src/Game/Payoff.jl

@@ -0,0 +1,62 @@
+
+"""
+A bilateral payoff can be computed as the sum of bilateral interactions between the players.
+
+If all players have bilateral payoffs, an equilibrium for the sampled game can be solved
+using its polymatrix representation, which yields an efficient formulation for the
+optimization problem.
+
+It is expected that each bilateral payoff type implements the `bilateral_payoff` function.
+The function must have the following signature:
+```julia
+function bilateral_payoff(Πp::MyBilateralPayoff, p::Integer, xp::Vector{<:Union{Real,VariableRef}}, k::Integer, xk::Vector{<:Real})
+    [...]
+end
+```
+Note that the `xp` argument can also be a vector of JuMP variable references. This is
+because the payoff functions are also used to build JuMP expressions, e.g., in
+`best_response` (see `src/Game/Player.jl`).
+
+"""
+abstract type AbstractBilateralPayoff <: AbstractPayoff end
+
+"Payoff function of player `p` with quadratic bilateral (pairwise) interactions."
+struct QuadraticPayoff <: AbstractBilateralPayoff
+    cp::Vector{Float64}
+    "`x[k]' * Qp[k] * x[p]` is the payoff component for player `p` with respect to the strategy of player `k`."
+    Qp::Vector{Matrix{Float64}}
+    function QuadraticPayoff(cp::Vector{<:Real}, Qp::Vector{<:Matrix{<:Real}})
+        if !all(Qpk->size(Qpk,2)==length(cp), Qp)
+            error("All Qp matrices must have the same number of columns as elements in cp (= dimension of p's strategy space).")
+        end
+        return new(cp, Qp)
+    end
+end
+function QuadraticPayoff(cp::Real, Qp::Vector{<:Real})
+    # constructor for simpler, unidimensional cases
+    return QuadraticPayoff([cp], [qpk * ones(1,1) for qpk in Qp])
+end
+
+"Compute each component of the payoff of player `p` with respect to player `k`."
+function bilateral_payoff(Πp::QuadraticPayoff, p::Integer, xp::Vector{<:Union{Real,VariableRef}}, k::Integer, xk::Vector{<:Real})
+    if p == k
+        return Πp.cp' * xp - 0.5 * xp' * Πp.Qp[p] * xp
+    else
+        return xk' * Πp.Qp[k] * xp
+    end
+end
+
+
+# Utils
+
+function bilateral_payoff(Πp::AbstractBilateralPayoff, p::Integer, xp::Vector{<:Union{Real,VariableRef}}, k::Integer, σk::DiscreteMixedStrategy)
+    return expected_value(xk -> bilateral_payoff(Πp, p, xp, k, xk), σk)
+end
+function bilateral_payoff(Πp::AbstractBilateralPayoff, p::Integer, σp::DiscreteMixedStrategy, k::Integer, σk::DiscreteMixedStrategy)
+    return expected_value(xp -> bilateral_payoff(Πp, p, xp, k, σk), σp)
+end
+
+"Compute the payoff of player `p` given strategies x."
+function payoff(Πp::AbstractBilateralPayoff, x::Vector{<:Vector{<:Union{Real,VariableRef}}}, p::Integer)
+    return sum([bilateral_payoff(Πp, p, x[p], k, x[k]) for k in 1:length(x)])
+end

src/Game/Payoff/BilateralPayoff.jl

@@ -0,0 +1,13 @@
+using JuMP
+
+abstract type AbstractPayoff end
+
+include("BilateralPayoff.jl")
+
+
+# Utils
+
+"Compute expected payoff for player `p` given mixed strategy profile `σ`."
+function payoff(Πp::AbstractPayoff, σ::Vector{DiscreteMixedStrategy}, p::Integer)
+    return expected_value(x -> payoff(Πp, x, p), σ)
+end

src/Game/Payoff/Payoff.jl

@@ -27,33 +27,14 @@ function set_optimizer(player::AbstractPlayer, optimizer_factory)
     JuMP.set_optimizer(player.Xp, optimizer_factory)
 end
 
-"Compute the utility that `player_p` receives from `player_k` when they play, resp., `xp` and `xk`."
-function bilateral_payoff(player_p::Player{QuadraticPayoff}, xp::Vector{<:Union{Real,VariableRef}}, player_k::Player{QuadraticPayoff}, xk::Vector{<:Real})
-    return bilateral_payoff(player_p.Πp, player_p.p, xp, player_k.p, xk)
-end
-"Compute the utility that `player_p` receives from `player_k` when they play, resp., `xp` and `σk`."
-function bilateral_payoff(player_p::Player{QuadraticPayoff}, xp::Vector{<:Union{Real,VariableRef}}, player_k::Player{QuadraticPayoff}, σk::DiscreteMixedStrategy)
-    return expected_value(xk -> bilateral_payoff(player_p.Πp, player_p.p, xp, player_k.p, xk), σk)
-end
-
-"Compute the payoff of player `player` given pure strategy profile `x`."
-function payoff(player::AbstractPlayer, x::Vector{<:Vector{<:Real}})
-    return payoff(player.Πp, x, player.p)
-end
-"Compute the payoff of player `player` given mixed strategy profile `σ`."
-function payoff(player::AbstractPlayer, σ::Vector{DiscreteMixedStrategy})
-    _payoff = x -> payoff(player, x)
-    return expected_value(_payoff, σ)
-end
-
 "Compute `player`'s best response to the mixed strategy profile `σp`."
-function best_response(player::Player{QuadraticPayoff}, σ::Vector{DiscreteMixedStrategy})
+function best_response(player::Player{<:AbstractPayoff}, σ::Vector{DiscreteMixedStrategy})
     xp = all_variables(player.Xp)
 
     # TODO: No idea why this doesn't work
     # @objective(model, Max, sum([IPG.bilateral_payoff(Πp, p, xp, k, σ[k]) for k in 1:m]))
 
-    obj = QuadExpr()
+    obj = AffExpr()
     for k in 1:length(σ)
         obj += IPG.bilateral_payoff(player.Πp, player.p, xp, k, σ[k])
     end

src/Game/Player.jl

@@ -1,9 +1,9 @@
 module IPG
 
-using NormalGames, JuMP
+using JuMP
 
 include("Game/Strategies.jl")
-include("Game/Payoff.jl")
+include("Game/Payoff/Payoff.jl")
 include("Game/Player.jl")
 include("SGM/SGM.jl")
 

src/IPG.jl

@@ -11,7 +11,7 @@ function find_deviation_best_response(players::Vector{<:AbstractPlayer}, σ::Vec
 
         new_σ = copy(σ)
         new_σ[p] = DiscreteMixedStrategy([1], [new_x_p])
-        payoff_improvement = payoff(player, new_σ) - payoff(player, σ)
+        payoff_improvement = payoff(player.Πp, new_σ, player.p) - payoff(player.Πp, σ, player.p)
         if payoff_improvement > dev_tol
             return payoff_improvement, p, new_x_p
         end

src/SGM/DeviationReaction.jl

@@ -1,6 +1,6 @@
 
 "Compute polymatrix for normal form game from sample of strategies."
-function get_polymatrix(players::Vector{<:AbstractPlayer}, S_X::Vector{<:Vector{<:Vector{<:Real}}})
+function get_polymatrix(players::Vector{<:Player{<:AbstractBilateralPayoff}}, S_X::Vector{<:Vector{<:Vector{<:Real}}})
     polymatrix = Dict{Tuple{Integer, Integer}, Matrix{Float64}}()
 
     # compute utility of each player `p` using strategy `i_p` against player `k` using strategy `i_k`
@@ -15,8 +15,8 @@ function get_polymatrix(players::Vector{<:AbstractPlayer}, S_X::Vector{<:Vector{
                 for i_p in 1:length(S_X[p])
                     for i_k in 1:length(S_X[k])
                         polymatrix[p,k][i_p,i_k] = (
-                            IPG.bilateral_payoff(players[p], S_X[p][i_p], players[k], S_X[k][i_k])
-                            + IPG.bilateral_payoff(players[p], S_X[p][i_p], players[p], S_X[p][i_p])
+                            IPG.bilateral_payoff(players[p].Πp, p, S_X[p][i_p], k, S_X[k][i_k])
+                            + IPG.bilateral_payoff(players[p].Πp, p, S_X[p][i_p], p, S_X[p][i_p])
                         )
                     end
                 end
@@ -27,59 +27,49 @@ function get_polymatrix(players::Vector{<:AbstractPlayer}, S_X::Vector{<:Vector{
     return polymatrix
 end
 
-"Wrapper for NormalGames.NormalGame that includes the sample of strategies."
+"Normal-form representation of the sampled game."
 mutable struct SampledGame
     S_X::Vector{Vector{Vector{Float64}}}  # sample of strategies (finite subset of the strategy space X)
-    normal_game::NormalGames.NormalGame
+    polymatrix::Dict{Tuple{Int, Int}, Matrix{Float64}}
 end
-function SampledGame(players::Vector{<:AbstractPlayer}, S_X::Vector{<:Vector{<:Vector{<:Real}}})
-    payoff_polymatrix = get_polymatrix(players, S_X)
-    normal_game = NormalGames.NormalGame(length(players), length.(S_X), payoff_polymatrix)
-
-    return SampledGame(S_X, normal_game)
+function SampledGame(players::Vector{<:Player{<:AbstractBilateralPayoff}}, S_X::Vector{<:Vector{<:Vector{<:Real}}})
+    return SampledGame(S_X, get_polymatrix(players, S_X))
 end
 
-function add_new_strategy!(sampled_game::SampledGame, players::Vector{<:AbstractPlayer}, new_xp::Vector{<:Real}, p::Integer)
+function add_new_strategy!(sg::SampledGame, players::Vector{<:Player{<:AbstractBilateralPayoff}}, new_xp::Vector{<:Real}, p::Integer)
     # first part is easy, just add the new strategy to the set
-    push!(sampled_game.S_X[p], new_xp)
-
-    n = sampled_game.normal_game.n
-    strat = sampled_game.normal_game.strat
-    polymatrix = sampled_game.normal_game.polymatrix
-
-    strat[p] += 1
+    push!(sg.S_X[p], new_xp)
+    strat = length.(sg.S_X)
 
     # now we need to update the normal game (polymatrix)
-    for (p1, p2) in keys(polymatrix)
+    for (p1, p2) in keys(sg.polymatrix)
         if p1 == p2 == p
             # add new row to polymatrix to store the utilities wrt the new strategy
-            polymatrix[p,p] = zeros(strat[p], strat[p])
+            sg.polymatrix[p,p] = zeros(strat[p], strat[p])
         elseif (p1 != p) & (p2 == p)
             # add new column to polymatrix to store the utilities wrt the new strategy
-            polymatrix[p1,p] = hcat(polymatrix[p1,p], zeros(strat[p1], 1))
+            sg.polymatrix[p1,p] = hcat(sg.polymatrix[p1,p], zeros(strat[p1], 1))
 
             for i in 1:strat[p1]
                 # compute utility of player `p1` using strategy `i` against the new strategy of player `p`
-                polymatrix[p1,p][i,end] = (
-                    IPG.bilateral_payoff(players[p1], sampled_game.S_X[p1][i], players[p], new_xp)
-                    + IPG.bilateral_payoff(players[p1], sampled_game.S_X[p1][i], players[p1], sampled_game.S_X[p1][i])
+                sg.polymatrix[p1,p][i,end] = (
+                    IPG.bilateral_payoff(players[p1].Πp, p1, sg.S_X[p1][i], p, new_xp)
+                    + IPG.bilateral_payoff(players[p1].Πp, p1, sg.S_X[p1][i], p1, sg.S_X[p1][i])
                 )
             end
         elseif (p1 == p) & (p2 != p)
             # add new row to polymatrix to store the utilities wrt the new strategy
-            polymatrix[p,p2] = vcat(polymatrix[p,p2], zeros(1, strat[p2]))
+            sg.polymatrix[p,p2] = vcat(sg.polymatrix[p,p2], zeros(1, strat[p2]))
 
             for i in 1:strat[p2]
                 # compute utility of player `p1` using strategy `i` against the new strategy of player `p`
-                polymatrix[p,p2][end,i] = (
-                    IPG.bilateral_payoff(players[p], new_xp, players[p2], sampled_game.S_X[p2][i])
-                    + IPG.bilateral_payoff(players[p], new_xp, players[p], new_xp)
+                sg.polymatrix[p,p2][end,i] = (
+                    IPG.bilateral_payoff(players[p].Πp, p, new_xp, p2, sg.S_X[p2][i])
+                    + IPG.bilateral_payoff(players[p].Πp, p, new_xp, p, new_xp)
                 )
             end
         end
     end
-
-    sampled_game.normal_game = NormalGames.NormalGame(n, strat, polymatrix)
 end
 
 include("SearchNE.jl")

src/SGM/SampledGame/SampledGame.jl

@@ -1,7 +1,14 @@
+using NormalGames
 
 "Compute a (mixed) nash equilibrium for the sampled game using PNS."
 function solve_PNS(sampled_game::SampledGame, optimizer_factory)::Vector{DiscreteMixedStrategy}
-    _, _, NE_mixed = NormalGames.NashEquilibriaPNS(sampled_game.normal_game, optimizer_factory, false, false, false)
+    normal_game = NormalGames.NormalGame(
+        length(sampled_game.S_X),
+        length.(sampled_game.S_X),
+        sampled_game.polymatrix
+    )
+
+    _, _, NE_mixed = NormalGames.NashEquilibriaPNS(normal_game, optimizer_factory, false, false, false)
 
     # each element in NE_mixed is a mixed NE, represented as a vector of probabilities in 
     # the same shape as S_X
@@ -17,15 +24,19 @@ end
 function solve_Sandholm1(sampled_game::SampledGame, optimizer_factory)::Vector{DiscreteMixedStrategy}
     # the method doesn't support polymatrices with negative entries, so a quick
     # preprocessing is performed
-    polymatrix = copy(sampled_game.normal_game.polymatrix)
+    polymatrix = copy(sampled_game.polymatrix)
     offset = 0
     for k in keys(polymatrix)
         offset = min(offset, minimum(polymatrix[k]))
     end
     for k in keys(polymatrix)
         polymatrix[k] = polymatrix[k] .- offset
     end
-    normal_game = NormalGames.NormalGame(sampled_game.normal_game.n, sampled_game.normal_game.strat, polymatrix)
+    normal_game = NormalGames.NormalGame(
+        length(sampled_game.S_X),
+        length.(sampled_game.S_X),
+        polymatrix
+    )
 
     _, _, NE_mixed, _, _, _, _, _ = NormalGames.NashEquilibria2(normal_game, optimizer_factory)