#OK to post 
#Johnny Lenn Fonseca, 04/12/2020, Assignment 21


#These are all of the procedures I coded.
Help := proc() print(`BRmixedRow(G,q), BRmixedCol(G,p), MNE(G), BRplayerj(G,j,others), multIndex(bounds), PNEn(M), blankNestedList(bounds), nRandGame(bounds, K1, K2), nRandGameGenPayOff(bounds,payOffLimits), biMatSplitter(M), biMatFuser(R,C), dotProd(v,w), BRmixedRowG(G,colProbs), BRmixedColG(G,rowProbs), MNEG(G)`) end:

#1
# Write procedures BRmixedRow(G,q) and BRmixedCol(G,p) for the set of Best Responses (now possibly infinite subses of the unit inerval) of The Row player to the mixed strategy [q,1-q] of the Column Player, an Column Player to the mixed stragegy of the Row player for a two-player game with two strategies available for each player.

BRmixedRow := proc(G,q) local q1, a:

q1 := solve(G[1,1,1]*a + G[1,2,1]*(1-a) - G[2,1,1]*a  - G[2,2,1]*(1-a), a):

if evalb((q1 = ()) and (q <> 1) and (q <> 0))then
   print(`One of the column player's strategies must be played with probability zero.`):
   RETURN(q1)
else
   if evalb(whattype(q) = symbol) then
      RETURN([{1, q > q1}, {2, q < q1}])
   elif type(q,numeric) then
      if evalb((q < 0) or (q > 1)) then
         print(`Invalid input for q.`):
         RETURN(FAIL)
      elif q = 0 then
         RETURN(BRrow(G,2)):
      elif q = 1 then
         RETURN(BRrow(G,1)):
      elif q > q1 then
         RETURN({1})
      else
          RETURN({2})
      fi:
   else
      print(`Invalid input for q.`):
      RETURN(FAIL)
   fi:
fi:
end:






BRmixedCol := proc(G,p) local p1, b:

p1 := solve(G[1,1,2]*b + G[2,1,2]*(1-b) - G[1,2,2]*b  - G[2,2,2]*(1-b), b):

if evalb((p1 = ()) and (p <> 1) and (p <> 0))then
   print(`One of the row player's strategies must be played with probability zero.`):
   RETURN(p1)
else
   if evalb(whattype(p) = symbol) then
      RETURN([{1, p > p1}, {2, p < p1}])
   elif type(p,numeric) then
      if evalb((p < 0) or (p > 1)) then
         print(`Invalid input for p.`):
         RETURN(FAIL)
      elif p = 0 then
         RETURN(BRcol(G,2)):
      elif p = 1 then
         RETURN(BRcol(G,1)):
      elif p > p1 then
         RETURN({1})
      else
          RETURN({2})
      fi:
   else
      print(`Invalid input for p.`):
      RETURN(FAIL)
   fi:
fi:
end:









#2 Use these to write a procedure MNE(G) that inputs a 2 by 2 bimatrix G and outputs the set of mixed Nash Equilibria. 

MNE := proc(G) local p, eR, q, eC, equiliForRow, equiliForCol, ei:

eR := BRmixedCol(G,p):
eC := BRmixedRow(G,q):
if evalb((eR = ()) or (eC = ())) then
   print(`No Mixed Nash Equilibrium exists.`): #Recall that one only exists if it is assumed every strategy is assigned NONZERO probability. 
   RETURN(FAIL)
else
   equiliForRow := lhs(eR[1,2]):
   equiliForCol := lhs(eC[1,2]):
   if evalb((0 < equiliForCol) and (1 > equiliForCol) and (0 < equiliForRow) and (1 > equiliForRow)) then
      RETURN([[equiliForRow, 1 - equiliForRow], [equiliForCol, 1 - equiliForCol]])
   else
      print(`No Mixed Nash Equilibrium exists.`): #Recall that one only exists if it is assumed every strategy is assigned NONZERO probability. 
      RETURN(FAIL)
   fi:
fi:
end:









#3
# [Optional Challenge] Write a procedure PNE3(G) that inputs a three-dimensional array describing a game with THREE players, where player ONE has A strategies, Player TWO has B strategies and Player THREE has C strategies and G is a list-of-lists-of-lists such that G[i][j][k] is the triple [PayOffForOne,PayOffForTwo,PayOffForThree] if they adopt strategies [i,j,k] respectively. What is PNE3(G) if
#
# G:= [ [ [[3,3,3],[0,4,2]], [[4,0,1],[1,1,2]] ], [ [[1,5,3],[2,-1,3]], [[2,0,-1],[-1,1,-2]] ] ]:
#
# Also write a procedure RandGameThree(A,B,C,K1,K2): A random 3-player game with A available strategies for Player ONE and B available strategies for Player TWO and C available strategies for Player TWO with payoffs between K1 and K2. Then try out PNE3(G) on many such random games. 

# Let us actually consider a game with n players!

#Inputs a game with n players, with player i having Ai strategies, and determines the best response for player j if the other n - 1 players choose strategies given by others := [a1,...a_{j-1}, a_{j +1},...a_{n}].
BRplayerj := proc(G,j,others) local numAjStrats, currentBest, bestStrats, allGood, k:
numAjStrats := 0:
if j = 1 then
   numAjStrats := nops(G):
else
   numAjStrats := nops(G[1$(j-1)]):
fi:
bestStrats := {1}:
currentBest := 0:
if j <> 1 then
   currentBest := G[op(others[1..j-1]),1,op(others[j..nops(others)]),j]:
   allGood := {currentBest}:
   for k from 2 to numAjStrats do
      if G[op(others[1..j-1]),k,op(others[j..nops(others)]),j] > currentBest then   
         currentBest := G[op(others[1..j-1]),k,op(others[j..nops(others)]),j]:
         bestStrats := {k}:
      elif G[op(others[1..j-1]),k,op(others[j..nops(others)]),j] = currentBest then
         allGood := {op(allGood), G[op(others[1..j-1]),k,op(others[j..nops(others)]),j]}:
         bestStrats := {op(bestStrats),k}:
      fi:
   od:
else
   currentBest := G[1,op(others[j..nops(others)]),1]:
   allGood := {currentBest}:
   for k from 2 to numAjStrats do
      if G[k,op(others),1] > currentBest then
         currentBest := G[k,op(others),1]:
         bestStrats := {k}:
      elif G[k,op(others),1] = currentBest then
         allGood := {op(allGood), G[k,op(others),1]}:
         bestStrats := {op(bestStrats),k}:
      fi:
   od:
fi:

RETURN(bestStrats)
end:






#Outputs all multi-indices with upper limits the elements of bounds e.g. if bounds := [3,2] then the output is all 6 pairs of indices starting from 1 in each component and running to [3,2].
multIndex := proc(bounds) local innerIndices, i:
if nops(bounds) = 1 then
   RETURN([seq(i, i = 1..bounds[1])])
else
   innerIndices := multIndex(bounds[2..nops(bounds)]):
   RETURN([seq(seq([i,op(innerIndices[j])], j = 1..nops(innerIndices)), i = 1..bounds[1])])
fi:
end:






#This determines Pure Nash Equilibria for n players where the i-th player has strategies A_i = {1,...,N_i}. In particular, it outputs the strategies each of the n players should play, AND their corresponding payoffs in doing so.
#Try it out with M := [[[[3, 3, 3], [0, 4, 2]], [[4, 0, 1], [1, 1, 2]]], [[[1, 5, 3], [2, -1, 3]], [[2, 0, -1], [-1, 1, -2]]], [[[2, 1, 3], [-1, 0, 1]], [[5, 1, 2], [-2, 7, -3]]]].
PNEn:=proc(M) local sizes, playerStrats, multiIndices, n, S, currentIndex, good, i, j, equilStratAndPayOffs:
sizes := [nops(M), seq(nops(M[1$j]), j = 1..2)]:
playerStrats := [seq([seq(i, i = 1..sizes[j])], j = 1..nops(sizes))]:
multiIndices := multIndex(sizes):
n := nops(playerStrats):
S := {}:
for i from 1 to nops(multiIndices) do
   currentIndex := multiIndices[i]:
   good := 0:
   if currentIndex[1] in BRplayerj(M,1,currentIndex[2..n]) then
      good := good + 1:
   fi:
   for j from 2 to n-1 do
      if currentIndex[j] in BRplayerj(M,j,[op(currentIndex[1..j - 1]),op(currentIndex[j + 1..n])]) then
         good := good + 1:
      fi:
   od:
   if currentIndex[n] in BRplayerj(M,n,currentIndex[1..n-1]) then
      good := good + 1:
   fi:
   if good = n then
      S := {op(S), currentIndex}:
   fi:
od:

equilStratAndPayOffs := {seq([S[i],G[op(S[i])]], i = 1..nops(S))}:

RETURN(equilStratAndPayOffs)
end:






#Outputs a sequence of nested lists with the final "innermost" element being of the form [0$bounds[-1]]
blankNestedList := proc(bounds) local i:
if nops(bounds) = 1 then
   RETURN([0$bounds[1]])
else
   RETURN([blankNestedList(bounds[2..nops(bounds)])$bounds[1]])
fi:

end:






#Generates a random game for n players with player i having bounds[i] many strategies, with payoffs ranging from K1 to K2.
#
nRandGame:=proc(bounds,K1,K2) local sizes, ra, initPayOffs, inds, i:
sizes := [op(bounds), nops(bounds)]:
ra:=rand(K1..K2):
initPayOffs := blankNestedList(sizes):
inds := multIndex(sizes):
for i from 1 to nops(inds) do
   initPayOffs[op(inds[i])] := ra():
od:
RETURN(initPayOffs)
end:






#The following generalizes the previous code by allowing player i to have random payoffs particular to player i i.e. the input payOffLimits is of the form
#payOffLimits := [seq([K1(i), K2(i)], i = 1..numPlayers)]:
#
#The code above is then equivalent to the general one with
#payOffLimits := [[K1,K2]$nops(bounds)]:
#
nRandGameGenPayOff:=proc(bounds,payOffLimits) local n, sizes, multiRAS, initPayOffs, inds, i:
n := nops(payOffLimits):
sizes := [op(bounds), n]:
multiRAS := [seq(rand(payOffLimits[i,1]..payOffLimits[i,2]), i = 1..n)]:

initPayOffs := blankNestedList(sizes):
inds := multIndex(sizes):
for i from 1 to nops(inds) do
   initPayOffs[op(inds[i])] := multiRAS[inds[i,n+1]]():
od:
RETURN(initPayOffs)
end:









#4
# [Optional (possibly difficult) Challenge] Write MNE(G) for an arbitrary two-person game, or at least for the 3 by 3 case. 
#

#"Splits" the input bimatrix into the two underlying matrices corresponding to the payoffs of each player in each possible choice of game strategy, but returned as lists.
biMatSplitter := proc(M) local R, C:
R := [seq([seq(M[i,j,1], j = 1..nops(M[1]))], i = 1..nops(M))]:
C := [seq([seq(M[i,j,2], j = 1..nops(M[1]))], i = 1..nops(M))]:
RETURN([R,C])
end:






#Combines two matrices, ASSUMED TO BE OF THE SAME SIZE, into a bimatrix.
biMatFuser := proc(R,C)
RETURN([seq([seq([R[i,j],C[i,j]], j = 1..nops(R[1]))], i = 1..nops(R))])
end:






dotProd := proc(v,w)
RETURN(sum(v[i]*w[i], i = 1..nops(v)))
end:






#This code determines any Mixed Nash Equilibria in a two-player game where player 1 has N1 strategies and player 2 has N2 strategies i.e. they can have the same number of strategies or differing amounts. However, note that this code assumes a 'nondegenerate' mixed strategy i.e. at least two of the probabilities are nonzero. Indeed, if only one was nonzero, then it would necessarily be with probability one. But then one is really just considering the Best Response of the other player in response to that particular action being chosen. To this end, see my procedure 'BRplayerj(G,j,others)' above for this case.

MNEG := proc(G) local sepMats, rowsMat, colsMat, rowProbs, colProbs, results, n, allPossC, sys, currIndex, currColProbs, nonzeroProbs, j, i, k, l, sys2, killedRows, goodRows, x, z, tempRes:

sepMats := biMatSplitter(G):
rowsMat := sepMats[1]:
colsMat := Transpose(sepMats[2]):

rowProbs := [seq(p[i], i = 1..nops(G))]:
colProbs := [seq(q[i], i = 1..nops(G[1]))]:
results := {}:

n := nops(colProbs):
allPossC := [seq(op(choose(n, k)), k = 0..n - 2)]:
for i from 1 to 2^n - n - 1 do
   sys := {}:
   currIndex := allPossC[i]:
   currColProbs := colProbs:
   for j from 1 to nops(currIndex) do
      currColProbs[currIndex[j]] := 0:
   od:

   nonzeroProbs := {seq(i, i = 1..n)} minus convert(currIndex,set):
   k := nonzeroProbs[1]:
   for j in ({seq(x, x = 1 ..n)} minus {k}) do
      if j in (nonzeroProbs minus {k}) then
         sys := {op(sys), dotProd(colsMat[k],rowProbs) = dotProd(colsMat[j],rowProbs)}:
      else
         sys := {op(sys), dotProd(colsMat[k],rowProbs) > dotProd(colsMat[j],rowProbs)}:
      fi:
   od:


   killedRows := {}:
   for l from 1 to nops(rowProbs) do
      if not evalb(`true` in [seq(evalb(rowsMat[l,t] >= sepMats[2,l,t]), t in nonzeroProbs)]) then
         sys := {op(sys), rowProbs[l] = 0}:
         killedRows := {op(killedRows),l}:
      fi:
   od:
   goodRows := {seq(i, i = 1..nops(rowsMat))} minus killedRows:
   sys := {op(sys), add(p[goodRows[i]], i = 1..nops(goodRows)) = 1, seq(op({p[goodRows[i]] > 0, p[goodRows[i]] < 1}), i = 1..nops(goodRows))}:


   if nops(goodRows) <> 0 then
      sys2 := {}:
      z := goodRows[1]:
      for j in ({seq(y, y = 1 ..nops(rowsMat))} minus {z}) do
         if j in (goodRows minus {z}) then
            sys2 := {op(sys2), dotProd(rowsMat[z],colProbs) = dotProd(rowsMat[j],colProbs)}:
         fi:
      od:
      sys2 := {op(sys2), add(q[i], i = 1..nops(G[1])) = 1, seq(q[currIndex[j]] = 0, j = 1..nops(currIndex)), seq(op({q[k] > 0, q[k] < 1}), k in nonzeroProbs)}:

      tempRes := [solve(sys), solve(sys2)]:
      if nops(tempRes) = 2 then
         results := {op(results),tempRes}:
      fi:
   fi:
od:

RETURN(results)
end:






































#########################
### Code from lecture ###
#########################


#C21.txt: April 9, 2020
##Nash Equilibrium

Help21:=proc(): print(`DB(), BRrow(M,j), BRcol(M,i) , PNE(M) , RandGame(A,B,K1,K2), ExGain(M,p,q) , All2Games()`):end:

with(combinat):

#DB(): A set of 5 classical 2-player games taken from the bible of Game Theory:
#A Course in Game Theory by Martin J. Osborne and Ariel Rubinstein [ a true masterpiece]
#https://sites.math.rutgers.edu/~zeilberg/EM20/OsborneRubinsteinMasterpiece.pdf
#

DB:=proc():
[
#Bach Or Stravisky, Fig. 16.1 (p. 16)
 [
  [[2,1],[0,0]],
  [[0,0],[1,2]]
 ],
#Mozart or Mahler, Fig. 16.2 (p. 16)
 [
  [[2,2],[0,0]],
  [[0,0],[1,1]]
 ],
#Prisoner's dillema, Fig. 17.1 (p. 17)
 [
  [[3,3],[0,4]],
  [[4,0],[1,1]]
 ],
#Hawk-Dove, Fig. 17.2 (p. 17)
 [
  [[3,3],[1,4]],
  [[4,1],[0,0]]
 ],
#Matching Pennies, Fig. 17.3 (p. 17)
 [
  [[1,-1],[-1,1]],
  [[-1,1],[1,-1]]
 ]

]
:
end:



#BRrow(M,j): Inputs
#(i) a bimatrix M , the so-called NORMAL FORM of a finite two-player game, that
#is  A by B list of lists where M[a][b] is the pair [PayoffForRowPlayer,PayoffForColPlayer]
#if Row Player adopts stragey a amd Col player adopts strategy b
#(ii) an integer j from 1 to B (One of the Column player's avaiable strategy
#outputs
#the Best Response(s) of the Row Player to that strategy. The set is usually a singleton, but in case of ties
#may have more than one element. For example, for the Best Response of the Row Player
#to Column's strategy 2  "Confess" type:
#BRrow(DB()[3],2);
BRrow:=proc(M,j) local rec,champs,A,B,i:
A:=nops(M): B:=nops(M[1]):

if not (1<=j and j<=B) then
 print(j, `is not a valid strategy for Col player`):
 RETURN(FAIL):
fi:

#champs is the current set of best responses to Col's strategy j, as we examine
#Row's possible responses, rec is the current record
champs:={1}:
rec:=M[1][j][1]:

for i from 2 to A do
if M[i][j][1]>rec then
#if strategy i is better than the previous best response it is now the only champion
  champs:={i}:
  rec:=M[i][j][1]:
elif M[i][j][1]=rec then
#the record is the same but there is a co-winnder
 champs:=champs union {i}:
fi:
od:
#We return the set of Champions
champs:
end:


#BRcol(M,i): Inputs
#(i) a bimatrix M , the so-called NORMAL FORM of a finite two-player game, that
#is  A by B list of lists where M[a][b] is the pair [PayoffForRowPlayer,PayoffForColPlayer]
#if Row Player adopts stragey a amd Col player adopts strategy b
#(ii) an integer i from 1 to A (One of the Row    player's avaiable strategy
#outputs
#the Best Response(s) of the Col Player to that strategy. The set is usually a singleton, but in case of ties
#may have more than one element. For example, for the Best Response of the Col Player
#to Row's strategy 2  "Confess" type:
#BRcol(DB()[3],1);
BRcol:=proc(M,i) local rec,champs,A,B,j:
A:=nops(M): B:=nops(M[1]):

if not (1<=i and i<=A) then
 print(i, `is not a valid strategy for Row player`):
 RETURN(FAIL):
fi:

#champs is the current set of best responses to Row's strategy i, as we examine
#Col's possible responses, rec is the current record
champs:={1}:
rec:=M[i][1][2]:

for j from 2 to B do
if M[i][j][2]>rec then
#if strategy j is better than the previous best response it is now the only champion
  champs:={j}:
  rec:=M[i][j][2]:
elif M[i][j][2]=rec then
#the record is the same but there is a co-winnder
 champs:=champs union {j}:
fi:
od:
#We return the set of Champions
champs:
end:





#PNE(M): Inputs
#a bimatrix M , the so-called NORMAL FORM of a finite two-player game, that
#is  A by B list of lists where M[a][b] is the pair [PayoffForRowPlayer,PayoffForColPlayer]
#if Row Player adopts stragey a amd Col player adopts strategy b
#outputs
#the (possibly empty) set of PURE NASH EQUILBRIA. Try:
#PNE(DB()[3]);

PNE:=proc(M) local A,B,i,j,S:
A:=nops(M):
B:=nops(M[1]):

#Row Player has A strategies {1, ..., A}
#Col Player has B strategies {1, ..., B}

#Recall that [i,j] is a Pure Nash Equilibrium if i is a member of BRcol(M,j) and j is a member of BRrow(M,i)
#Let S be the set of Pure Nash Equilibria
S:={};

for i from 1 to A do
 for j from 1 to B do

   if member(i, BRrow(M,j)) and member(j, BRcol(M,i)) then
    S:=S union {[i,j]}:
   fi:
 od:
od:

S:

end:

#RandGame(A,B,K1,K2): A random 2-player game 
#with A available strategies for the Row Player
#and  B available strategies for the Column Player
#with random pay-offs from K1 to K2, Try:
#RandGame(2,2,-1,1);
RandGame:=proc(A,B,K1,K2) local i,j,ra:
ra:=rand(K1..K2):

[seq([seq([ra(),ra()],j=1..B)],i=1..A)]:

end:

#ExGain(M,p,q): inputs a bi-matrix of  2 by 2 bi-matrix describing the payoffs of Row and Column player
#outputs the par [ExpectedGainForRow,ExpectedGainForCol] if 
#Player Row adopts strategy 1 with prob. p (and strategy 2 with prob. 1-p)
#Player Col adopts strategy 1 with prob. q (and strategy 2 with prob. 1-q)
#Note that these are polynomials of degree 1 in p and in q separately
ExGain:=proc(M,p,q) local k:

if not (nops(M)=2 and nops(M[1])=2) then
 print(`Not 2 by 2`):
 RETURN(FAIL):
fi:

expand([seq(M[1][1][k]*p*q+M[1][2][k]*p*(1-q)+M[2][1][k]*(1-p)*q+M[2][2][k]*(1-p)*(1-q),k=1..2)]):

end:


##Added right before class, APril 9, 2020
#All2Games(): All 2-person games with different payoffs for each player from {1,2,3,4}. Try:
#All2Games();
All2Games:=proc() local P,pi1,pi2:

P:=permute(4):

{
seq(
  seq(
     [
       [ [pi1[1],pi2[1]],[pi1[2],pi2[2]] ],
       [ [pi1[3],pi2[3]],[pi1[4],pi2[4]] ]
     ],
pi1 in P),
pi2 in P)
}:

end: