#!/usr/local/bin/maple
# -*- maplev -*-

# Nathaniel Shar
# HW 20
# Experimental Mathematics

# It is okay to link to this assignment on the course webpage.

# From c20.txt:

#April 8, 2013, Nash Equilibria for general two player games

Help:=proc(): print(`MaxF(F,x,n) , BRrow(A,y) , BRcol(B,x) `):
print(`PNE(A,B) , GenRanMat(m,n,K) , tra(A) `):
end:

#Given two matrices A and B m by n
#PayOff for A: xAy Payoff for B: xBy (x is a row vector of
#size m and y is a column of size n)
#(x,y) is a Nash Equilibrium if x is in BestResponse(A,y)
#and y is in BestResponse(B,x)

#MaxF(F,x,n): inputs sumbol x, positive integer n
#and a linear combination of x[1], ..., x[n]
#outputs the set of i such that x[i]=1 makes F max
#For example,
#MaxF(3*x[1]+4*x[2]+4*x[3]+2*X[4],x,4);
MaxF:=proc(F,x,n) local i,champs,rec,nathan:

champs:={1}:
rec:=coeff(F,x[1],1):

for i from 2 to n do

nathan:=coeff(F,x[i],1):

 if nathan=rec then
   champs:=champs union  {i}:
 elif nathan>rec then
    champs:={i}:
   rec:=nathan:
 fi:

od:

champs:

end:

#BRrow(A,y): The best response by the ROW player
#(whose payoff matrix is the m by n matrix A
#to the strategy y of the column player
#output is a set of pure strategies
BRrow:=proc(A,y) local i,j,F,x,m,n:
m:=nops(A):
n:=nops(A[1]):
F:=add(add(x[i]*A[i][j]*y[j],i=1..m),j=1..n):

MaxF(F,x,m):

end:

#BRcol(B,x): The best response by the COL player
#(whose payoff matrix is the m by n matrix B
#to the strategy x of the row player
#output is a set of pure strategies
BRcol:=proc(B,x) local i,j,F,y,m,n:
m:=nops(B):
n:=nops(B[1]):
F:=add(add(x[i]*B[i][j]*y[j],i=1..m),j=1..n):

MaxF(F,y,n):

end:


#PNE(A,B)
PNE:=proc(A,B) local m,n,nash,i,j,john:

m:=nops(A):
n:=nops(A[1]):


nash:={}:

for j from 1 to n do
john:=BRrow(A,[0$(j-1),1,0$(n-j)]):

  for i in john do

     if j in BRcol(B,[0$(i-1),1,0$(m-i)]) then
       nash:=nash union {[i,j]}:
     fi:
  od:
od:

nash:

end:

#GenRanMat(m,n,K): a random m by n matrix with
#entries from 0 to K
GenRanMat:=proc(m,n,K) local i,j,tim:

tim:=rand(0..K):

[seq([seq(tim(),j=1..n)],i=1..m)]:

end:


#transpose of A
tra:=proc(A) local i,j,m,n:
m:=nops(A):
n:=nops(A[1]):
[seq([seq(A[j][i],j=1..m)],i=1..n)]:
end:

#############
# Problem 1 #
#############

GatherStat := proc(n,K,L,t) local i, j, k, r, matrices:
    r := rand(0..K):
    matrices := [seq([seq([seq(r(), i=1..n)], j=1..n)], k=1..L)]:
    return [seq(numboccur(i, map(x->nops(PNE(x, tra(x))), matrices)), i=0..n^2)]:
end:

#############
# Problem 2 #
#############

# The distribution appears to be roughly normal (or at least
# bell-shaped). I would conjecture that the mean approaches 2 as n
# approaches infinity, since each diagonal entry is a Nash equilibrium
# with probability roughly 1/n, and each off-diagonal entry is a Nash
# equilibrium with probability roughly 1/n^2.

#############
# Problem 5 #
#############

# Not sure if this actually works in general, but it worked in the
# cases I tested and it should definitely work for 2x2.

NE := proc(A, B) local S, i, m, n:
    m := nops(A):
    n := nops(A[1]):
    S := PNE(A,B):
    S := map(x->[[seq(`if`(x[1]=i, 1, 0), i=1..m)], [seq(`if`(x[2]=i, 1, 0), i=1..n)]], S):
    S := S union MNE(A,B):
    return S:
end:


# Find mixed Nash equilibria

MNE:=proc(A,B) local m,n,i, j, eqs, vars, solns, L, r:
    m:=nops(A):
    n:=nops(A[1]):
    L := NonDominated(A, B, {seq(i, i=1..m)}, {seq(j, j=1..n)}):
    print(L):
    vars := {seq(x[i], i=L[1]), seq(y[i], i=L[2])}:
    eqs := {seq(add(A[L[1][i]][j]*y[j], j=L[2]) = add(A[L[1][i+1]][j]*y[j], j=L[2]), i=1..nops(L[1])-1)}:
    eqs := eqs union {seq(add(B[i][L[2][j]]*x[i], i=L[1]) = add(B[i][L[2][j+1]]*x[i], i=L[1]), j=1..nops(L[2])-1)}:
    eqs := eqs union {add(x[i], i=L[1]) = 1, add(y[i], i=L[2]) = 1}:
    solns := solve(eqs, vars):
    print(eqs):
    r := {[subs(solns, [seq(x[i], i=1..m)]), subs(solns, [seq(y[i], i=1..n)])]}:
    return subs({seq(x[i] = 0, i=1..m), seq(y[i]=0, i=1..n)}, r):
end:

NonDominated := proc(A, B, row, col) local i, j, k, S, T:
    S := {}:
    T := {}:
    for i in row do:
        for j in row minus {i} do:
            if {true} = {true, seq(evalb(A[i][k] < A[j][k]), k in col)} then:
                S := S union {i}:
            fi:
        od:
    od:
    for i in col do:
        for j in col minus {i} do:
            if {true} = {true, seq(evalb(B[k][i] < B[k][j]), k in row)} then:
                T := T union {i}:
            fi:
        od:
    od:
    if S = {} and T = {} then:
        return [row, col]:
    else:
        return NonDominated(A, B, row minus S, col minus T):
    fi:
end:


#############
# Problem 6 #
#############

LetsCooperate := proc(A, B, i, j) local m, n, a, b, S:
    m := nops(A):
    n := nops(B):
    S := {}:
    for a from 1 to m do:
        for b from 1 to n do:
            if a <> i and b <> j and A[a][b] > A[i][j] and B[a][b] > B[i][j] then:
                S := S union {[a,b]}:
            fi:
        od:
    od:
    return S:
end: