Compute_NE.py

from Instances import *
from Initial_str import *
# control time
from time import time
from copy import deepcopy
import csv

import numpy as np

import gurobipy as grb


# under unix: to limit time
# import signal
#
# def signal_handler(signum, frame):
#     raise Exception("Timed out!")

# m-SGM
def IterativeSG(G, max_iter, opt_solver=1, S=[]):
    r"""Create instances in a standard format.

    Parameters:
    ----------
    G: Game class (see Instances file)
    max_iter: maximum number of sampled games to be solved
    opt_solver: 0 if cplex is used and 1 otherwise (use gurobi); in the paper it is always 1.
    S: intial set of strategies (optinal)
    Returns:
    -------
    ne: list of vectors indicating the Nash equilibrium associated with the strategies in S
    Profits: List of profits for each player under ne
    S: final set of strategies
    count: number of iterations, i.e., sampled games solved
    numb_back: number of backtracking
    cpu_time: computational time
    """
    # create log file
    # f = open("log_ComputeOneNE.txt",'w')
    # f.close()
    # STEP 0 - INITIALIZATION
    # initialize set of strategies
    S, U_p, Best_m = InitialStrategies(G, opt_solver)
    S_new = [[] for p in range(G.m())]
    if [[]] in S:
        print("ERROR: There is a player without feasible strategies")
        return [], [], S, 0, 0, 0
    Numb_stra = [1 for p in range(G.m())]
    U_depend = [[{} for k in range(G.m())] for p in range(G.m())]
    # set mandatory strategy in the support
    Numb_stra_S = [0] * G.m()
    M_pre = [G.m() - 1, 0, deepcopy(Numb_stra_S), deepcopy(S_new), [[]] * G.m(), deepcopy(U_depend), deepcopy(S_new),
             [0, 0], 0]  # strategy S[0] must be in the supp of player m
    # set memory
    Memory = [deepcopy(M_pre)]
    Back = False  # start computation of equilibria from the start
    numb_back = 0
    # STEP 2 - COMPUTE EQUILIBRIA OF RESTRICTED (SAMPLED) GAME
    # compute Nash equilibrium taking account M and Back
    count = 1
    U_depend = Utilities_Poymatrix(G.m(), G.Q(), U_depend, S_new, S, Numb_stra_S)
    # M_pos = [player p, best response of p, Numb_stra of current node, strategies S of current node, U_p of S, U_depend of S, new best responses from S]
    M_pos = [None, None, deepcopy(Numb_stra), deepcopy(S), deepcopy(U_p), deepcopy(U_depend), deepcopy(S_new), None, 0]
    list_best = list(range(G.m()))
    time_aux = time()
    ne = []
    ne_previous = ne[:]
    while True and count <= max_iter and time() - time_aux <= 3600:
        # f = open("log_ComputeOneNE.txt",'a')
        # f.write("\n\n Processing node ... "+str(count))
        # f.write("\nComputing equilibria of game of size "+str(Numb_stra)+" .... \n")
        # f.close()
        print("\n\n Processing node ... ", count)
        print("Computing equilibria.... \n")
        ne_previous = ne[:]
        # under unix: to limit time
        # signal.signal(signal.SIGALRM, signal_handler)
        # signal.alarm(3600-int(time()-time_aux))   #  seconds
        try:
            ne, Profits, S0_j = ComputeNE(M_pre, Back, U_depend, U_p, G.m(), G.n_I(), G.n_C(), Numb_stra, opt_solver,
                                          M_pos[2], M_pos[-1])
        # except Exception, msg: python 2.7
        except Exception:
            print("Time limit exceeded")
            return ne_previous, [], S, count, numb_back, time() - time_aux
        if ne == []:  # fail to compute such equilibrium
            # f = open("log_ComputeOneNE.txt",'a')
            # f.write("\nFail to compute equilibrium\n#############################\n Backtracking ....\n#############################\n")
            # f.close()
            print(" Fail to compute equilibrium")
            print("#############################")
            print("\n Backtracking ....\n")
            print("#############################")
            S_new = deepcopy(M_pre[-3])
            U_p = deepcopy(M_pre[4])
            S = deepcopy(M_pre[3])
            U_depend = deepcopy(M_pre[5])
            Numb_stra_S = deepcopy(M_pre[2])
            U_depend = Utilities_Poymatrix(G.m(), G.Q(), U_depend, S, S_new, Numb_stra_S)
            U_p, S = IndUtilities(G.m(), G.c(), G.Q(), S, U_p, S_new)
            Numb_stra = [Numb_stra_S[p] + len(S_new[p]) for p in range(G.m())]  # update number of strategies
            M_pos = deepcopy(M_pre)
            Memory.pop()
            M_pre = deepcopy(Memory[-1])
            Back = True
            numb_back = numb_back + 1
            # count = count -1
        else:
            # f = open("log_ComputeOneNE.txt",'a')
            # f.write("\n Equilibrium computed sucessfully\n")
            # f.close()
            print(" Equilibrium computed sucessfully")
            # M_pre[2] = ne
            # Memory[-1] = M
            Back = False  # ne found
            aux = True  # no player has incentive to deviate
            S_new = [[] for p in range(G.m())]  # set of new strategies to be considered
            Profile = [np.array([sum(ne[int(k + sum(Numb_stra[:p]))] * S[p][k][i] for k in range(Numb_stra[p])) for i in
                                 range(G.n_I()[p] + G.n_C()[p])]) for p in range(G.m())]
            aux_p = 0
            while aux and aux_p < G.m():  # FEED BEST RESPONSES WITH NE solution
                p = list_best[aux_p]
                s_p, u_max, _ = BestReactionGurobiNASP(G.name(), G.m(), G.c()[p], G.Q()[p], Profile, p, False,
                                                       Best_m[p])
                if abs(Profits[p] - u_max) >= 1e-5:
                    # print("Player",p)
                    # print("\tBR",s_p)
                    # print("\tIncumbent", Profile[p])
                    aux = False
                    S_new[p].append(s_p)
                    Numb_stra_S = deepcopy(Numb_stra)
                    # update M_pos
                    M_pos[0] = p
                    M_pos[1] = Numb_stra[p]
                    M_pos[-2] = deepcopy(ne)
                    M_pos[-3][p].append(s_p)
                    M_pos[-1] = S0_j
                    Memory.append(deepcopy(M_pos))
                    Numb_stra[p] = Numb_stra[p] + 1
                    M_pre = deepcopy(M_pos)
                    U_depend = Utilities_Poymatrix(G.m(), G.Q(), U_depend, S, S_new, Numb_stra_S)
                    U_p, S = IndUtilities(G.m(), G.c(), G.Q(), S, U_p, S_new)
                    S_new = [[] for _ in range(G.m())]
                    M_pos = [None, None, deepcopy(Numb_stra), deepcopy(S), deepcopy(U_p), deepcopy(U_depend),
                             deepcopy(S_new), None, 0]
                    # print(list_best)
                    list_best.append(p)
                    list_best = list_best[:aux_p] + list_best[aux_p + 1:]
                aux_p = aux_p + 1
            # count = count+1
            if aux:
                # f = open("log_ComputeOneNE.txt",'a')
                # f.write("\n Total time: "+str(time()-cpu_time))
                # f.close()
                return ne, Profits, S, count, numb_back, time() - time_aux
        count = count + 1
    if time() - time_aux > 3600:
        print("Time Limit Exceeded")
    else:
        print(" Maximum number of iterations was attained")
    return ne_previous, [], S, count, numb_back, time() - time_aux


#######################################################################################################################

#######################################################
##        COMPUTE INDIVIDUAL PROFITS                 ##
#######################################################

# INPUT
# m = number of players
# c = linear objective function coefficients for each player (list of vectors)
# S = list of strategies for each player
# U_p = list of individual profits for each player
# S_new = new strategies to be added to S and to compute individual profit

# OUTPUT
# U_p = list of players individual profits
# S = new set of strategies
def IndUtilities(m, c, Q, S, U_p, S_new):
    for p in range(m):
        for s in S_new[p]:
            U_p[p].append(float(np.dot(c[p], s) - 0.5 * np.dot(s, np.dot(Q[p][p], s))))
            S[p].append(s)
    return U_p, S


#######################################################################################################################

#######################################################
##        POLYMATRIX PART OF THE PROFITS             ##
#######################################################

# INPUT
# m = number of players
# Q = bilinear coefficients in the objective function for each player (list of matrices)
# p = player for which we are fixing the strategy
# U_p = list of individual profits for each player
# U_depend = list of the players' profit
# S = strategies of each player (list)
# s = profile of strategies being fixed
# numb = last strategy fixed
# Numb_stra_S = number of strategies in S[p]
# OUTPUT
# U_depend = matrice of utilities (in fact it is a dictionary)
def Utilities_Poymatrix(m, Q, U_depend, S, S_new, Numb_stra_S):
    for p in range(m):
        for k in range(p + 1, m):
            for sp in enumerate(S_new[p]):
                for sk in enumerate(S[k] + S_new[k]):
                    U_depend[p][k][(Numb_stra_S[p] + sp[0], sk[0])] = float(np.dot(sk[1], np.dot(Q[p][k], sp[1])))
                    U_depend[k][p][(sk[0], Numb_stra_S[p] + sp[0])] = float(np.dot(sp[1], np.dot(Q[k][p], sk[1])))
        for k in range(p):
            for sp in enumerate(S_new[p]):
                for sk in enumerate(S[k]):
                    U_depend[p][k][(Numb_stra_S[p] + sp[0], sk[0])] = float(np.dot(sk[1], np.dot(Q[p][k], sp[1])))
                    U_depend[k][p][(sk[0], Numb_stra_S[p] + sp[0])] = float(np.dot(sp[1], np.dot(Q[k][p], sk[1])))
    return U_depend


#######################################################################################################################

#######################################################
##        COMPUTE Nash Equilibrium                   ##
#######################################################

# INPUT
# S = set of strategies for each player (list)
# M = (p, numb, sigma)
# Back = computation continues from previous computed equilibrium sigma (if Back = True)
# U_depend = polymatrix
# U_p = individual profits
# m = number of players
# n_I = number of binary variables for each player (list)
# n_C = number of continuous variables for each player (list)
# Numb_stra = size of S; number of strategies available for each player (list)
# opt_solver = 0 then use CLEP, = 1 then use Gurobi
# Supp_Stra = M_pos[2] strategies to consider in the support(new strategies should not be considered = S_new of M_pos)
# OUTPUT
# ne = a Nash equilibrium with strategy S[p][numb] of player p in the support

from itertools import combinations_with_replacement, combinations, product, chain


def ComputeNE(M, Back, U_depend, U_p, m, n_I, n_C, Numb_stra, opt_solver, Supp_Stra, start):
    ##### HEURISTIC 6 ####
    size_pre_ne = [sum(1 for j in M[-2][int(sum(M[2][:p])):int(sum(M[2][:p + 1]))] if j > 10 ** -5) for p in range(m)]
    S0 = Heuristic6(Supp_Stra, m, size_pre_ne, n_I, n_C)
    ##### HEURISTIC 6 ####
    for S0_j, s0 in enumerate(S0[start:]):
        # now we have the support sizes
        # if Back is true, how can we start s0?
        A_supp = [None for p in range(m)]
        # ATTENTION: THE SETS IN D ARE NOT SORTED
        D = []
        for p in range(m):
            # D[p] is the set of strategies in the support of player p with size s0[p]
            # from S[p] take s0[p] strategies
            if p != M[0]:
                # D.append([candidate for candidate in combinations(range(Supp_Stra[p]),s0[p])])
                # HEURISTIC II
                D.append([candidate for candidate in combinations(
                    HeuristicII(Supp_Stra[p], M[-2][int(sum(M[2][:p])):int(sum(M[2][:p + 1]))], -1, M[-3][p]), s0[p])])
            else:  # M[1] must be in the support of player M[0]
                # D.append([candidate+(M[1],) for candidate in combinations(range(M[1])+range(M[1]+1,Supp_Stra[p]),s0[p]-1)])
                # HEURISTIC II: give priority to strategies choosen in the previous ne
                D.append([candidate + (M[1],) for candidate in combinations(
                    HeuristicII(Supp_Stra[p], M[-2][int(sum(M[2][:p])):int(sum(M[2][:p + 1]))], M[1], M[-3][p]),
                    s0[p] - 1)])
        ne, Profits = Recursive_Backtracking(m, A_supp, D, 0, U_depend, U_p, opt_solver, Numb_stra)
        if ne != []:  # Nash equilibrium found!
            return ne, Profits, start + S0_j
    return [], [], start + S0_j


def ComputeNE_NOT_DFS(U_depend, U_p, m, n_I, n_C, Numb_stra, opt_solver, ne_previous, deviator):
    ##### HEURISTIC 6 ####
    size_pre_ne = [sum(1 for j in ne_previous[int(sum(Numb_stra[:p])):int(sum(Numb_stra[:p + 1]))] if j > 10 ** -5) for
                   p in range(deviator)] + [sum(
        1 for j in ne_previous[int(sum(Numb_stra[:deviator])):int(sum(Numb_stra[:deviator + 1]) - 1)] if
        j > 10 ** -5)] + [
                      sum(1 for j in ne_previous[int(sum(Numb_stra[:p])):int(sum(Numb_stra[:p + 1]))] if j > 10 ** -5)
                      for p in range(deviator + 1, m)]
    S0 = Heuristic6(Numb_stra, m, size_pre_ne, n_I, n_C)
    ##### HEURISTIC 6 ####
    Numb_stra_previous = deepcopy(Numb_stra)
    Numb_stra_previous[deviator] = Numb_stra_previous[deviator] - 1
    for S0_j, s0 in enumerate(S0):
        # now we have the support sizes
        # if Back is true, how can we start s0?
        A_supp = [None for p in range(m)]
        # ATTENTION: THE SETS IN D ARE NOT SORTED
        D = []
        for p in range(m):
            # D[p] is the set of strategies in the support of player p with size s0[p]
            # from S[p] take s0[p] strategies
            # HEURISTIC II
            if p != deviator:
                D.append([candidate for candidate in combinations(HeuristicII(Numb_stra[p], ne_previous[int(sum(
                    Numb_stra_previous[:p])):int(sum(Numb_stra_previous[:p + 1]))], -1, []), s0[p])])
            else:
                D.append([candidate for candidate in combinations(HeuristicII(Numb_stra[p], ne_previous[int(sum(
                    Numb_stra_previous[:p])):int(sum(Numb_stra_previous[:p + 1]))], -1, [1]), s0[p])])
        ne, Profits = Recursive_Backtracking(m, A_supp, D, 0, U_depend, U_p, opt_solver, Numb_stra)
        if ne != []:  # Nash equilibrium found!
            return ne, Profits
    return [], []


def HeuristicII(Supp_Stra_p, M_ne, M_1, S_new_p):
    order_index_str = chain(range(M_1), range(M_1 + 1, Supp_Stra_p))
    M_ne_aux = M_ne + [0 for _ in S_new_p]
    if M_ne_aux != []:
        return sorted(order_index_str, key=lambda x: M_ne_aux[x])
    else:
        return order_index_str


def Heuristic6(Supp_Stra, m, size_ne, n_I, n_C):
    S0 = []
    # it is n_I[p]+n_C[p]+1 in case the objectives are linear: an optimum is in a facet which has dimension n_I+n_C-1
    # and thus, any point of it can be written as a convex combinatiuon of n_I+n_C extreme points of that facet
    for s0 in product(*[range(1, min(Supp_Stra[p] + 1, n_I[p] + n_C[p] + 2)) for p in range(m)]):
        S0.append(list(s0))
    if m == 2:
        return sorted(S0, key=lambda x: (abs(x[0] - x[1]), max(abs(size_ne[0] - x[0]), abs(size_ne[1] - x[1])),
                                         max(abs(size_ne[0] + 1 - x[0]), abs(size_ne[1] + 1 - x[1])), x[0] + x[1]))
    else:
        return sorted(S0, key=lambda x: (
            max(abs(size_ne[p] - x[p]) for p in range(m)), max(abs(size_ne[p] + 1 - x[p]) for p in range(m)), sum(x),
            max(abs(x[i] - x[j]) for i in range(m) for j in range(i, m))))


#######################################################################################################################

#######################################################
##        RECURSIVE BACKTRACKING                     ##
#######################################################

# INPUT
# m = number of players
# A_supp = set of strategies in the support for each player (list)
# D = candidates to be a support for each player (list)
# i = player for whom we are fixing the strategy
# S = set of strategies for each player in the restricted game (list)
# U_depend = polymatrix of utilities
# U_p = indepedent utilities
# opt_solver = 0 then use CPLEX, = 1 then use Gurobi
# Numb_stra = number strategies for each player (list)
# OUPUT
# ne - Nash equilibrium (list)
# Profits - profits for each player in the equilibrium ne (list)

def Recursive_Backtracking(m, A_supp, D, i, U_depend, U_p, opt_solver, Numb_stra):
    if i == m:  # this is, we have fixed the support for each player
        # Solve Feasibility Problem
        return FeasibilityProblem(m, A_supp, U_depend, U_p, opt_solver, Numb_stra)  # ne, Profits
    else:
        while D[i] != []:
            d_i = D[i].pop()  # remove d_i from D[i]
            A_supp[i] = d_i
            D_new = RS([[A_supp[p]] for p in range(i + 1)] + deepcopy(D[i + 1:]), Numb_stra, U_depend, U_p, m)
            if D_new != None:
                ne, Profits = Recursive_Backtracking(m, deepcopy(A_supp), deepcopy(D_new), i + 1, U_depend, U_p,
                                                     opt_solver, Numb_stra)
                if ne != []:
                    return ne, Profits
    return [], []


#######################################################################################################################

#######################################################
##        FEASIBILITY PROBLEM                       ##
#######################################################

# INPUT
# m = number of players
# A_supp  = strategies to which each player associates positive probability (list)
# U_depend = polymatrix of utilities
# U_p = independent utilities
# opt_solver = 0 then use CPLEX, = 1 then use Gurobi
# Numb_stra = number of strategies for each player (list)
# OUTPUT
# ne = Nash equilibrium (list)
# Profits = profit of each player for the equilibrium ne (list)

def FeasibilityProblem(m, A_supp, U_depend, U_p, opt_solver, Numb_stra):
    return FeasibilityProblem_Gurobi(m, A_supp, U_depend, U_p, Numb_stra)


def FeasibilityProblem_Gurobi(m, A_supp, U_depend, U_p, Numb_stra, m_p=None):
    # print "\n\n Solving Problem with Supports: ", A_supp
    if m_p == None:
        # initiate model
        m_p = grb.Model("FeasibilityProblem")
        m_p.setParam("Threads", 2)
        # no pritting of the output
        m_p.setParam('OutputFlag', False)
        # set objective function direction
        m_p.ModelSense = -1  # maximize
        m_p.update()
        # probability variables
        sigma = [{sp: m_p.addVar(lb=0, vtype="C", name="sigma_" + str(p) + "_" + str(sp)) for sp in A_supp[p]} for p in
                 range(m)]
        m_p.update()

        ########################################################################################################
        ############# WHEN FEASIBILITY PROBLEM HAS MORE THAN ONE SOLUTION ######################################
        ###### MAXIMIZE THE NUMBER OF VARIABLES WITH POSITIVE PROBABILITY ######################################
        # aux = [m_p.addVar(obj = 1, lb=0,vtype="C",name="aux_"+str(p)) for p in range(m)] # aux <= sigma_p_sp
        # m_p.update()
        # for p, sp in enumerate(A_supp):
        #     for s in sp:
        #         m_p.addConstr(aux[p] <= sigma[p][s])
        #         m_p.update()
        ########################################################################################################
        ########################################################################################################

        # profit variables
        v = [m_p.addVar(lb=-1 * grb.GRB.INFINITY, vtype="C", name="v_" + str(p)) for p in range(m)]
        m_p.update()
        for p in range(m):
            m_p.addConstr(grb.quicksum(sigma[p].values()) == 1)
            m_p.update()
        for p, S_p in enumerate(Numb_stra):
            for sp in range(S_p):
                if sp in A_supp[p]:
                    m_p.addConstr(U_p[p][sp] + grb.quicksum(
                        sigma[k][sk] * U_depend[p][k][(sp, sk)] for k in range(m) if k != p for sk in A_supp[k]) == v[
                                      p])
                    m_p.update()
                else:
                    m_p.addConstr(U_p[p][sp] + grb.quicksum(
                        sigma[k][sk] * U_depend[p][k][(sp, sk)] for k in range(m) if k != p for sk in A_supp[k]) <= v[
                                      p])
                    m_p.update()
        # m_p.write("apagar.lp")
        m_p.optimize()
        ne = []
        Profits = []
        # print "Solution status for Feasibility Problem: ", m_p.status
        if m_p.status not in [3, 4]:
            for p, sp in enumerate(Numb_stra):
                for j in range(sp):
                    if j in A_supp[p]:
                        ne.append(sigma[p][j].x)
                    else:
                        ne.append(0)
                Profits.append(v[p].x)
        return ne, Profits


#######################################################################################################################

#######################################################
##    REMOVAL OF STRICTLY DOMINATED STRATEGIES       ##
#######################################################

# INPUT
# D = list of strategies to which each player is restricted to play (list; D[p] is a list that contains tuples = sets of strategies)
# Numb_stra = strategies available for each player (list)
# U_depend = polymatrix of utilities
# U_p = indepedent utilities
# OUTPUT
# D_new =  D_new[p] strategies that are not strictly dominated given D[-p]

# REMARK: conditionally dominated
# sp in S[p] is conditionally dominated given a profile of sets of available actions R[-p] contained in S[-p],
# if the following conditions holds:
# there is sp' in Sp[p], forall s[-p] in R[-p]: Profit[p](sp,s[-p]) < Profit[p](sp',s[-p])

def RS(D, Numb_stra, U_depend, U_p, m):
    changed = True
    while changed:
        changed = False
        for p in range(m):
            for a in set(dp for Ap in D[p] for dp in Ap):  # for all pure strategies in the possible supports
                for a_prime in chain(range(a), range(a + 1, Numb_stra[p])):
                    aux = True  # it is conditionally dominated
                    # all possible outcomes
                    Outcomes_minus_p = [set(d_minus_p for A_minus_p in D[k] for d_minus_p in A_minus_p) for k in
                                        range(p)] + [set([0])] + [
                                           set(d_minus_p for A_minus_p in D[k] for d_minus_p in A_minus_p) for k in
                                           range(p + 1, m)]
                    for s in product(*Outcomes_minus_p):
                        # if a is conditionally dominated by a_prime given D[-p]:
                        if U_p[p][a] + sum(U_depend[p][k][(a, s[k])] for k in range(m) if k != p) >= U_p[p][
                            a_prime] + sum(U_depend[p][k][(a_prime, s[k])] for k in range(m) if k != p):
                            aux = False
                            break
                    if aux:  # it is dominated
                        D_new_p = []
                        for Ap in D[p]:
                            if a not in Ap:
                                D_new_p.append(Ap)
                        D[p] = D_new_p
                        changed = True
                        if D[p] == []:
                            # print "### Not solving feasibility problem "
                            return None
    return D


if __name__ == "__main__":

    f = open('resultsNASP.csv', 'a')
    writer = csv.writer(f)
    writer.writerow(['Instance', 'iterations', 'cpu_time', 'status'])
    f.flush()

    for instance in range(149):
        np.random.seed(6)
        G = Game('NASP', NASPfile="./Instances/NASP/Instance_" + str(instance+1))

        # DFS: EXECUTE m-SGM; max numb of iterations 50
        ne, Profits_mSGM, S, numb_iter, numb_back, cpu_time = IterativeSG(G, 700)
        print("NE:", ne)
        print("Profits", Profits_mSGM)
        for p in range(G.m()):
            print("Strategy of", p, "is:", S[p])
        status = "YES" if len(Profits_mSGM) > 1 else "NO"
        writer.writerow([instance+1, numb_iter, cpu_time, status])
        f.flush()