External variables as parameters; evaluate NonlinearExpr through recursive replace

Parameters

Author: brunompacheco

src/Game/Assignment.jl

@@ -0,0 +1,73 @@
+
+const AssignmentDict = Dict{VariableRef,Float64}
+
+"""
+Create a dictionary of variable assignments (JuMP-style) from a pure strategy.
+"""
+Assignment(player::Player, x::PureStrategy) = AssignmentDict(zip(all_variables(player), x))
+Assignment(x::Profile{PureStrategy})::AssignmentDict = merge(collect(Assignment(p, x_p) for (p, x_p) in x)...)
+export Assignment
+
+"Replace the variables in an expression with their assigned values."
+function replace(expr::AbstractJuMPScalar, assignment::AssignmentDict)
+    _recursive_replace(expr::Number) = expr
+    function _recursive_replace(expr::AbstractJuMPScalar)
+        if expr isa VariableRef
+            return get(assignment, expr, expr)
+        elseif expr isa AffExpr
+            # TODO: this behavior of value is unintended and may lead to problems in the future
+            return value(v -> get(assignment, v, v), expr)
+        elseif expr isa QuadExpr
+            return value(v -> get(assignment, v, v), expr)
+        elseif expr isa NonlinearExpr
+            replaced_expr = NonlinearExpr(expr.head, Vector{Any}(map(_recursive_replace, expr.args)))
+            # If there are no nonlinear arguments, we can try to simplify the resulting expression
+            if !any(arg isa NonlinearExpr for arg in replaced_expr.args)
+                # TODO: the following is to be replaced by JuMP.simplify once https://github.com/jump-dev/JuMP.jl/pull/4047 gets merged
+                g = MOI.Nonlinear.SymbolicAD.simplify(moi_function(replaced_expr))
+
+                # TODO: this owner model assignment is really ugly. maybe the owner model should be an argument
+                terms_in_replaced_expr = filter(v -> v isa JuMP.AbstractJuMPScalar, replaced_expr.args)
+                owner = nothing
+                for term in terms_in_replaced_expr
+                    term_owner = owner_model(term)
+                    if ~isnothing(term_owner)
+                        owner = term_owner
+                        break
+                    end
+                end
+
+                if isnothing(owner)
+                    # Explicitly handle the case where g is not a JuMP expression
+                    if g isa Number
+                        replaced_expr = g
+                    else
+                        error("replace: Unable to convert simplified expression to a number or JuMP expression. Got: $(typeof(g))")
+                    end
+                else
+                    # update owner model
+                    replaced_expr = jump_function(owner, g)
+                end
+            end
+            return replaced_expr
+        else
+            return expr
+        end
+    end
+
+    return _recursive_replace(expr)
+end
+
+"Translate variable references of the assignment to internal references."
+function _internalize_assignment(player::Player, assignment::AssignmentDict)
+    internal_assignment = AssignmentDict()
+    for (v_ref, v_val) in assignment
+        if v_ref ∈ all_variables(player.X)
+            internal_assignment[v_ref] = v_val
+        elseif v_ref ∈ keys(player._param_dict)
+            internal_assignment[player._param_dict[v_ref]] = v_val
+        end
+    end
+
+    return internal_assignment
+end

src/Game/BestResponse.jl

@@ -1,9 +1,9 @@
 
 "Compute `player`'s best response to the mixed strategy profile `σ_others`."
 function best_response(player::Player, σ_others::Profile{DiscreteMixedStrategy})
-    vars_player = all_variables(player.X)
+    @assert player ∉ keys(σ_others) "Player must not be in the profile of others."
 
-    obj = payoff(player, vars_player, σ_others)
+    obj = expected_value(x_others -> replace_in_payoff(player, Assignment(x_others)), σ_others)
 
     # I don't know why, but it was raising an error without changing the sense to feasibility first
     set_objective_sense(player.X, JuMP.MOI.FEASIBILITY_SENSE)
@@ -12,32 +12,8 @@ function best_response(player::Player, σ_others::Profile{DiscreteMixedStrategy}
     set_silent(player.X)
     optimize!(player.X)
 
-    return value.(vars_player)
+    return value.(all_variables(player))
 end
+
 "Compute `player`'s best response to the pure strategy profile `x_others`."
 best_response(player::Player, x_others::Profile{PureStrategy}) = best_response(player, convert(Profile{DiscreteMixedStrategy}, x_others))
-# function best_response(player::Player{<:AbstractBilateralPayoff}, σ::Vector{DiscreteMixedStrategy})
-    # error("best_response for player with bilateral payoff not implemented yet") # TODO
-
-    # xp = all_variables(player.X)
-
-    # # 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 = AffExpr()
-    # for k in eachindex(σ)
-    #     if k == player.p
-    #         obj += IPG.bilateral_payoff(player.Π, xp)
-    #     else
-    #         obj += IPG.bilateral_payoff(player.Π, xp, σ[k], k)
-    #     end
-    # end
-    # # I don't know why, but it was raising an error without changing the sense to feasibility first
-    # set_objective_sense(player.X, JuMP.MOI.FEASIBILITY_SENSE)
-    # @objective(player.X, JuMP.MOI.MAX_SENSE, obj)
-
-    # set_silent(player.X)
-    # optimize!(player.X)
-
-    # return value.(xp)
-# end

src/Game/Game.jl

@@ -2,5 +2,6 @@
 include("Strategies.jl")
 include("Player.jl")
 include("Profile.jl")
+include("Assignment.jl")
 include("Payoff.jl")
 include("BestResponse.jl")

src/Game/Payoff.jl

@@ -1,40 +1,32 @@
 
-"""
-Create a dictionary of variable assignments (JuMP-style) from a pure strategy.
-"""
-function build_var_assignments(player::Player, x::Vector{<:Any})
-    variable_assignments = Dict{VariableRef, Any}()
-    for (v, val) in zip(all_variables(player.X), x)
-        variable_assignments[v] = val
-    end
-    return variable_assignments
+function replace_in_payoff(player::Player, assignment::AssignmentDict)
+    internal_assignment = _internalize_assignment(player, assignment)
+    return replace(player.Π, internal_assignment)
 end
 
 """
 Get the payoff map for `player` given the pure strategy profile `x_others`.
 The payoff map is a function that takes the player's strategy and returns the payoff.
 """
 function get_payoff_map(player::Player, x_others::Profile{PureStrategy})
-    others_var_assignments = [build_var_assignments(other, x_other)
-                              for (other, x_other) in x_others]
-    var_assignments = merge(others_var_assignments...)
-
-    function payoff_map(x_player::Vector{<:Any})
-        complete_var_assignments = merge(var_assignments, build_var_assignments(player, x_player))
+    assignment_others = Assignment(x_others)
+    internal_assignment_others = _internalize_assignment(player, assignment_others)
 
-        return value(v -> complete_var_assignments[v], player.Π)
+    function payoff_map(x_player::Vector{Float64})
+        complete_assignment = merge(internal_assignment_others, Assignment(player, x_player))
+        return value(v -> complete_assignment[v], player.Π)
     end
 
     return payoff_map
 end
 
 "Evaluate the player's payoff when she plays `x_player` and the others play `x_others`."
-function payoff(player::Player, x_player::Vector{<:Any}, x_others::Profile{PureStrategy})
+function payoff(player::Player, x_player::PureStrategy, x_others::Profile{PureStrategy})
     return get_payoff_map(player, x_others)(x_player)
 end
 
 "Expected payoff of a pure strategy (`x_player`) against a mixed profile (`σ_others`)."
-function payoff(player::Player, x_player::Vector{<:Any}, σ_others::Profile{DiscreteMixedStrategy})
+function payoff(player::Player, x_player::PureStrategy, σ_others::Profile{DiscreteMixedStrategy})
     return expected_value(x_others -> payoff(player, x_player, x_others), σ_others)
 end
 

src/Game/Player.jl

@@ -1,20 +1,70 @@
 using JuMP, JSON3
 
+const VarToParamDict = Dict{VariableRef,VariableRef}
+
 "A player in an IPG."
 mutable struct Player
     "Strategy space."
     X::Model
     # TODO: using value(...) to manipulate expressions does not work for NonlinearExpr, see https://github.com/jump-dev/JuMP.jl/issues/4044 for an appropriate solution (another huge refactor))
     "Payoff expression."
     Π::AbstractJuMPScalar
+    _param_dict::VarToParamDict
 end
-function Player()
-    return Player(Model(), AffExpr(0))
+Player() = Player(Model(), AffExpr(NaN), VarToParamDict())
+Player() = Player(Model(), AffExpr(0.0), VarToParamDict())
+Player(X::Model) = Player(X, AffExpr(0.0), VarToParamDict())
+function Player(X::Model, Π::AbstractJuMPScalar) 
+    player = Player(X)
+    set_payoff!(player, Π)
+    return player
 end
 export Player
 
+JuMP.all_variables(p::Player) = filter(v -> !is_parameter(v), all_variables(p.X))
+
+"Maps external variables to internal parameters. Creates a new parameter if it does not exist."
+function _maybe_create_parameter_for_external_var(player::Player, var::VariableRef)::VariableRef
+    var ∈ all_variables(player.X) && return var
+
+    if !haskey(player._param_dict, var)
+        # create anonymous parameter with the same name as the variable
+        param = @variable(player.X, base_name=name(var), set=Parameter(1))
+
+        player._param_dict[var] = param
+    end
+
+    return player._param_dict[var]
+end
+
 function set_payoff!(player::Player, payoff::AbstractJuMPScalar)
-    player.Π = payoff
+    _recursive_internalize_expr(expr::Number) = expr
+    function _recursive_internalize_expr(expr::AbstractJuMPScalar)::AbstractJuMPScalar
+        if expr isa VariableRef
+            return _maybe_create_parameter_for_external_var(player, expr)
+        elseif expr isa AffExpr
+            internal_terms = typeof(expr.terms)(
+                _maybe_create_parameter_for_external_var(player, var) => coeff
+                for (var, coeff) in expr.terms
+            )
+            return AffExpr(expr.constant, internal_terms)
+        elseif expr isa QuadExpr
+            internal_terms = typeof(expr.terms)(
+                UnorderedPair{VariableRef}(
+                    _maybe_create_parameter_for_external_var(player, vars.a),
+                    _maybe_create_parameter_for_external_var(player, vars.b)
+                ) => coeff
+                for (vars, coeff) in expr.terms
+            )
+            return QuadExpr(_recursive_internalize_expr(expr.aff), internal_terms)
+        elseif expr isa NonlinearExpr
+            return NonlinearExpr(expr.head, Vector{Any}(map(_recursive_internalize_expr, expr.args)))
+        else
+            error("Unsupported expression type: $(typeof(expr))")
+        end
+    end
+
+    player.Π = _recursive_internalize_expr(payoff)
 end
 function set_payoff!(player::Player, payoff::Real)
     player.Π = AffExpr(payoff)
@@ -39,7 +89,7 @@ function find_feasible_pure_strategy(player::Player)::PureStrategy
     set_silent(player.X)
     optimize!(player.X)
 
-    return value.(all_variables(player.X))
+    return value.(all_variables(player))
 end
 
 "Solve the feasibility problem of all players, returning a feasible profile."

src/SGM/DeviationReaction.jl

@@ -5,9 +5,10 @@ function find_deviation_best_response(players::Vector{Player}, σ::Profile{Discr
 
     for p in player_order
         player = players[p]
-        new_x_p = best_response(player, σ)
         σ_others = others(σ, player)
 
+        new_x_p = best_response(player, σ_others)
+
         payoff_improvement = payoff(player, new_x_p, σ_others) - payoff(player, σ)
         if payoff_improvement > dev_tol
             return payoff_improvement, player, new_x_p

src/SGM/Initialization.jl

@@ -5,7 +5,7 @@ empty_S_X(players::Vector{Player}) = Dict{Player, Vector{PureStrategy}}(p => Vec
 function initialize_strategies_feasibility(players::Vector{Player})
     S_X = empty_S_X(players)
     for player in players
-        xp_init = start_value.(all_variables(player.X))
+        xp_init = start_value.(all_variables(player))
 
         if nothing in xp_init
             # TODO: if `initial_sol` is just a partial solution, I could fix its values
@@ -24,13 +24,13 @@ function initialize_strategies_player_alone(players::Vector{Player})
     S_X = empty_S_X(players)
 
     # profile that simulates players being alone (all others play 0)
-    x_others_dummy = Profile{PureStrategy}(player => zeros(length(all_variables(player.X))) for player in players)
+    x_dummy = Profile{PureStrategy}(player => zeros(length(all_variables(player))) for player in players)
 
     for player in players
-        xp_init = start_value.(all_variables(player.X))
+        xp_init = start_value.(all_variables(player))
 
         if nothing in xp_init
-            xp_init = best_response(player, x_others_dummy)
+            xp_init = best_response(player, others(x_dummy, player))
         end
 
         push!(S_X[player], xp_init)

src/SGM/PolymatrixGame/Polymatrix.jl

@@ -3,28 +3,60 @@ using LinearAlgebra
 
 const Polymatrix = Dict{Tuple{Player, Player}, Matrix{Float64}}
 
-"Compute the component of the payoff that doesn't depend on other players."
-function compute_self_payoff(p::Player, x_p::PureStrategy)
-    var_assignments_p = build_var_assignments(p, x_p)
 
-    self_linear_payoff = sum(get(p.Π.aff.terms, v, 0) * var_assignments_p[v] for v in all_variables(p.X))
+function compute_self_payoff(Π::AffExpr, v_bar::AssignmentDict)::Float64
     # note the get() may be necessary as there may not be terms for all variables
-    self_affine_payoff = p.Π.aff.constant + self_linear_payoff
+    self_linear_payoff = sum(get(Π.terms, ref, 0) * val for (ref, val) in v_bar)
+
+    # constant is here by convention (inherited from NormalGames.jl)
+    return Π.constant + self_linear_payoff
+end
+
+function compute_self_payoff(Π::QuadExpr, v_bar::AssignmentDict)::Float64
     # TODO: maybe Dict{VariableRef, Number} should be the standard for assignments
-    self_quad_payoff = sum(get(p.Π.terms, UnorderedPair(v,v), 0) * var_assignments_p[v]^2 for v in all_variables(p.X))
+    self_quad_payoff = sum(get(Π.terms, UnorderedPair(ref,ref), 0) * val^2 for (ref, val) in v_bar)
+
+    return self_quad_payoff + compute_self_payoff(Π.aff, v_bar)
+end
+
+"Compute the component of the payoff that doesn't depend on other players."
+function compute_self_payoff(p::Player, x_p::PureStrategy)
+    v_bar = Assignment(p, x_p)
 
-    return self_affine_payoff + self_quad_payoff
+    return compute_self_payoff(p.Π, v_bar)
+end
+
+"Compute player p's payoff component from some *other* player playing v_bar_k."
+function compute_others_payoff(Π::AffExpr, v_bar_k::AssignmentDict)::Float64
+    # TODO: identical to compute_self_payoff(::AffExpr, ::AssignmentDict), except for the constant. I could refactor
+    return sum(  # terms of the form q_j * xk_j, where xk_j belongs to player k
+        get(Π.terms, ref, 0) * val for (ref, val) in v_bar_k
+    )
+end
+
+"Compute player p's payoff component from her playing v_bar_p and some *other* player playing v_bar_k."
+function compute_bilateral_payoff(Π::QuadExpr, v_bar_p::AssignmentDict, v_bar_k::AssignmentDict)::Float64
+    mixed_components = sum(  # terms of the form q_ij * xp_i * xk_j, where xp_i belongs to player p and xk_j belongs to player k
+        get(Π.terms, UnorderedPair(ref_i,ref_j), 0) * val_i * val_j
+        for (ref_i, val_i) in v_bar_p
+        for (ref_j, val_j) in v_bar_k
+    )
+    other_components = sum(  # terms of the form q_ij * xk_i * xk_j, where xk_i,xk_j belong to player k
+        get(Π.terms, UnorderedPair(ref_i,ref_j), 0) * val_i * val_j
+        for (ref_i, val_i) in v_bar_k
+        for (ref_j, val_j) in v_bar_k
+    )
+    other_components = other_components / 2  # I'm iterating over al possible unordered pairs twice!
+
+    return mixed_components + other_components + compute_others_payoff(Π.aff, v_bar_k)
 end
 
 function compute_bilateral_payoff(p::Player, x_p::PureStrategy, k::Player, x_k::PureStrategy)
-    var_assignments_k = build_var_assignments(k, x_k)
-    var_assignments_p = build_var_assignments(p, x_p)  # TODO: this could be cached.
     # In fact, +1 for having Dict{VariableRef, Number} as the standard for assignments
+    v_bar_p = Assignment(p, x_p)  # TODO: this could be cached.
+    v_bar_k = _internalize_assignment(p, Assignment(k, x_k))
 
-    return sum(
-        get(p.Π.terms, UnorderedPair(vp,vk), 0) * var_assignments_p[vp] * var_assignments_k[vk]
-        for vp in all_variables(p.X), vk in all_variables(k.X)
-    )
+    return compute_bilateral_payoff(p.Π, v_bar_p, v_bar_k)
 end
 
 "Compute polymatrix for normal form game from sample of strategies."

src/SGM/SGM.jl

@@ -20,7 +20,6 @@ function SGM(players::Vector{Player}, optimizer_factory; max_iter=100, dev_tol=1
     # the strategy space of each player or, in case there is none, a feasibility problem is
     # solved
 
-    # TODO: I should use S_X for the sampled game, rather than the strategies
     S_X = initialize_strategies(players)
     sampled_game = PolymatrixSampledGame(players, S_X)
     verbose && println("Game initialized with strategies: ", S_X)