#Tim Naumovitz
#hw20.txt
#I give permission to post this.

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

GatherStat:=proc(n,K,L,t) local M,i,j,k,A,N,count:

M:=[seq(GenRanMat(n,n,K), j=1..L)]:
N:=[]:
for i from 0 to t do
 count:=0:
 for A in M do
  if nops(PNE(A,tra(A)))=i then
   count:=count+1:
  fi:
 od:
N:=[op(N),count]:
od:
N:
end:

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

#> GatherStat(3, 5, 1000, 3);
#                      [29, 204, 369, 270]
#> GatherStat(2, 5, 1000, 3);
#                       [0, 356, 463, 146]
#> GatherStat(2, 10, 1000, 3);
#                       [0, 389, 531, 67]
#> GatherStat(3, 10, 1000, 3);
#                      [57, 270, 379, 217]
#> GatherStat(3, 5, 1000, 3);
#                      [28, 177, 387, 264]
#> GatherStat(3, 10, 1000, 3);
#                      [52, 264, 406, 222]
#> GatherStat(2, 5, 10000, 3);
#                     [0, 3508, 4840, 1378]

#It appears that for 2x2 matrices we never get
#0 pure nash equilibria, but for 3x3 matrices, we
#get a small number of 0 NE matrices(<1%).
#As K gets larger, we're more likely to have fewer
#NE.  We get 1 NE somewhere around 30-35% of the time
# and 2 NE around  40-45% of the time.

###################Problem 3###################

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

m:=2:
n:=2:


nash:={}:

for i from 0 to K do
john:=BRrow(A,[i/K,1-i/K]):

  for j in john do

     P:=B[j][1]*i/K+B[j][2]*(1-i/K):
     if B[j][BRcol(B,[0$(j-1),1,0$(m-j)])]=P then
       nash:=nash union {[[0$(j-1),1,0$(m-j)],[i/K,1-i/K]]}:
     fi:
  od:
od:

nash:

end:

###################Problem 4###################

#NE3(A,B,K)
NE3:=proc(A,B,K) local m,n,nash,i,j,k,r,john,P:

m:=3:
n:=3:


nash:={}:

for i from 0 to K do
 for j from 0 to K-i do
  k:=(K-i-j)/K:
john:=BRrow(A,[i/K,j/K,k/K]):

  for r in john do

     P:=B[r][1]*i/K+B[r][2]*(j/K)+B[r][3]*k/K:
     if B[r][BRcol(B,[0$(r-1),1,0$(m-r)])]=P then
       nash:=nash union {[[0$(r-1),1,0$(m-r)],[i/K,j/K,k/K]]}:
     fi:
  od:
od:

nash:

end:

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

LetsCooporate:=proc(A,B,i,j) local m,n,r,s,S:

m:=nops(A):
n:=nops(A[1]):
S:={}:

for r from 1 to m do
 for s from 1 to n do
  if A[r][s]>A[i][j] and B[r][s]>B[i][j] then
   S:=S union {[r,s]}:
  fi:
 od:
od:
S:
end:

###################Old Stuff###################

#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: