#OK to post homework
#Blair Seidler, 2/6/22, Assignment 5
with(combinat):
with(ListTools):
with(plots):

Help:=proc(): print(`BeatNE(a,b), EstBeatNE(a,b,K,M)`): end:

#1.
(*
> P1 := (p1, p2) -> (a - p1 + b*p2)*(p1 - c);
P1 := proc (p1, p2) options operator, arrow; (a-p1+b*p2)*(p1-c) 

   end proc
> P2 := (p1, p2) -> (a - p2 + b*p1)*(p2 - c);
P2 := proc (p1, p2) options operator, arrow; (a-p2+b*p1)*(p2-c) 

   end proc
> dP1 := diff(P1(p1, p2), p1);
                   dP1 := b p2 + a + c - 2 p1
> solve(dP1 = 0, p1);
                       1        1     1  
                       - b p2 + - a + - c
                       2        2     2  
> dP2 := diff(P2(p1, p2), p2);
                   dP2 := b p1 + a + c - 2 p2
> solve(dP2 = 0, p2);
                       1        1     1  
                       - b p1 + - a + - c
                       2        2     2  
> solve({dP1 = 0, dP2 = 0}, {p1, p2});
                  /       a + c         a + c\ 
                 { p1 = - -----, p2 = - ----- }
                  \       b - 2         b - 2/ 

Confirming the results in the book.

Running the command

> implicitplot([dP1 = 0, dP2 = 0], p1 = 0 .. 100, p2 = 0 .. 100, color = ["NavyBlue", "Teal"], legend = [P1, P2])

with a=100, b=0.8, c=2 produced the graph in the attached file hw5p1BlairSeidler.png, which confirms
the solution as (85,85). Note that (100+2)/(2-0.8)=85.

*)

#2. 
(*
> G := RG(30, 30, 10000);
> BestTot(G);
                        {[5, 14]}, 19661
> NE(G);
                               {}


> G := RG(30, 30, 10000);
> BestTot(G);
                        {[6, 18]}, 18956
> NE(G);
                               {}



> G := RG(30, 30, 10000);
> BestTot(G);
                       {[29, 26]}, 19310
> NE(G);
                           {[7, 28]}
> BetterForBoth(G, 7, 28);
                               {}
*)

#3.
#BeatNE(a,b): inputs positive integers a and b and outputs the set of all games 
#(given as bi-matrices) for which there is a unique Nash Equilibrium, AND there exists
#a strategy choice that is better for BOTH players.
BeatNE:=proc(a,b) local AG,G,neset,OUT:
AG:=AllGames(a,b):
OUT:={}:
for G in AG do
  neset:=NE(G):
  if nops(neset)=1 then
    if nops(BetterForBoth(G,op(neset[1])))>0 then
      OUT:=OUT union {G}:
    fi:
  fi:
od:
OUT:
end:

(* Output

> BeatNE(2, 2);
           {[[[1, 1], [3, 2]], [[4, 3], [2, 4]]], 

            [[[1, 1], [3, 4]], [[2, 3], [4, 2]]], 

            [[[1, 2], [3, 4]], [[2, 3], [4, 1]]], 

            [[[1, 3], [3, 4]], [[2, 2], [4, 1]]], 

            [[[1, 4], [3, 3]], [[2, 2], [4, 1]]], 

            [[[1, 4], [4, 3]], [[2, 2], [3, 1]]], 

            [[[1, 4], [4, 3]], [[3, 2], [2, 1]]], 

            [[[2, 1], [3, 2]], [[4, 3], [1, 4]]], 

            [[[2, 2], [3, 1]], [[1, 4], [4, 3]]], 

            [[[2, 2], [4, 1]], [[1, 3], [3, 4]]], 

            [[[2, 2], [4, 1]], [[1, 4], [3, 3]]], 

            [[[2, 3], [4, 1]], [[1, 2], [3, 4]]], 

            [[[2, 3], [4, 2]], [[1, 1], [3, 4]]], 

            [[[2, 4], [4, 3]], [[3, 2], [1, 1]]], 

            [[[3, 1], [2, 2]], [[4, 3], [1, 4]]], 

            [[[3, 2], [1, 1]], [[2, 4], [4, 3]]], 

            [[[3, 2], [2, 1]], [[1, 4], [4, 3]]], 

            [[[3, 3], [1, 4]], [[4, 1], [2, 2]]], 

            [[[3, 4], [1, 1]], [[4, 2], [2, 3]]], 

            [[[3, 4], [1, 2]], [[4, 1], [2, 3]]], 

            [[[3, 4], [1, 3]], [[4, 1], [2, 2]]], 

            [[[4, 1], [2, 2]], [[3, 3], [1, 4]]], 

            [[[4, 1], [2, 2]], [[3, 4], [1, 3]]], 

            [[[4, 1], [2, 3]], [[3, 4], [1, 2]]], 

            [[[4, 2], [2, 3]], [[3, 4], [1, 1]]], 

            [[[4, 3], [1, 4]], [[2, 1], [3, 2]]], 

            [[[4, 3], [1, 4]], [[3, 1], [2, 2]]], 

            [[[4, 3], [2, 4]], [[1, 1], [3, 2]]]}

There are 28 of these games (out of 576 total games). For each of them, there is a unique NE,
but the opposite corner of the bimatrix produces a better result for both. So if the NE is
[2,1], then [1,2] produces a better result for both players. This is the only location where
the better for both result can be, otherwise the NE wouldn't be a NE. This is also true in
Prisoner's Dilemma, where the NE is [2,2] but [1,1] is better for both players.


I did run BeatNE(3, 2), but there are 32400 games which I will not list here. Of those games,
there are 5400 with each of the six combined strategies as the unique NE. Again, none of the
cells in the same row or column of the bimatrix can be better for both because otherwise it
would not be a NE.

For the 5400 games with unique NE at [1,1], 2400 are better for both at [2,2], 2400 are better
for both at [3,2], and 600 are better for both at [2,2] and [3,2].

*)

#4. EstBeatNE(a,b,K,M): inputs a,b,K (as for RG(a,b,K) and a large integer M, generates M random 
#games, and counts, out of these M games, how many have the property that they have exactly one 
#Nash Equilibria, followed by the number of those, among the former, that have the property that 
#the total pay-offs is better than the total pay-off of the Nash Equilibria, followed by the 
#number of those (games with exactly one NE) that have the property that there exists a strategy 
#choice that is NOT a NE but nevertheless is better for both of them (like (Silent, Silent) in the 
#Prisoner's Dilemma). Divide these numbers by M, and take evalf
EstBeatNE:=proc(a,b,K,M) local AG,G,bt,neset,has1ne,hasbettot,hasbetboth:
AG:={seq(RG(a,b,K),i=1..M)}:
has1ne:=0:
hasbettot:=0:
hasbetboth:=0:
for G in AG do
  neset:=NE(G):
  if nops(neset)=1 then
    has1ne:=has1ne+1:
    if nops(BetterForBoth(G,op(neset[1])))>0 then
      hasbetboth:=hasbetboth+1:
    fi:
    bt:=BestTot(G):
    if not neset[1] in bt[1] then
      hasbettot:=hasbettot+1:
    fi:
  fi:
od:
evalf([has1ne,hasbettot,hasbetboth]/M):
end:

(* Output:

The numbers I was seeing with M=1000 were all over the place, so I boosted M to 10000. This
produced much more consistent results:

> for i to 5 do EstBeatNE(10, 10, 1000, 10000); end do;

           [0.4115, 0.2100, 0.1053]
           [0.4068, 0.2134, 0.1043]
           [0.4117, 0.2121, 0.1052]
           [0.4143, 0.2137, 0.1073]
           [0.4041, 0.2040, 0.0966]

*)

#5.
(*
> u1 := q1*(1 - q1 - q2)^d1; u2 := q2*(1 - q1 - q2)^d2;
                                          d1
                    u1 := q1 (1 - q1 - q2)  
                                          d2
                    u2 := q2 (1 - q1 - q2)  
> u1p := diff(u1, q1);
                                                  d1   
                             d1   q1 (1 - q1 - q2)   d1
         u1p := (1 - q1 - q2)   - ---------------------
                                       1 - q1 - q2     
> u2p := diff(u2, q2);
                                                  d2   
                             d2   q2 (1 - q1 - q2)   d2
         u2p := (1 - q1 - q2)   - ---------------------
                                       1 - q1 - q2     
> solve({u1p = 0, u2p = 0}, {q1, q2});
          /           d2                    d1       \ 
         { q1 = ---------------, q2 = --------------- }
          \     d1 d2 + d1 + d2       d1 d2 + d1 + d2/ 


When d1=d2=1, this yields the known NE (1/3,1/3).

*)

#### Code included from hw3 ####
#AllGames(a,b): inputs positive integers a and b and outputs the set of (a*b)!2
#Games where each player's payoffs are distinct and drawn from {1, ..., a*b}
AllGames:=proc(a,b) local p1:
p1 := permute(a*b):
{seq(seq(
#the below line of code was copied from RandDisGame in C3.txt
[seq([seq([pi1[b*i1+j1],pi2[b*i1+j1]],j1=1..b)],i1=0..a-1)] ,
pi2 in p1), pi1 in p1)}:
end:


#### Code included from C5.txt ####
#C5.txt Maple Code for Lecture 5

Help5:=proc(): print(` IsNE(G,a1,a2), NE(G), BestTot(G) , BetterForBoth(G,a1,a2), RG(a,b,K) `):end:

#BetterForBoth(G,a1,a2): Given a game G and a strategy choice (a1,a2) finds all the strategy choices that are better for BOTH players
BetterForBoth:=proc(G,a1,a2) local b1,b2,S:
S:={}:

for b1 from 1 to nops(G) do
 for b2 from 1 to nops(G[1]) do
  if G[b1][b2][1]>G[a1][a2][1] and G[b1][b2][2]>G[a1][a2][2] then
   S:=S union {[b1,b2]}:
  fi:
 od:
od:
S:

end:

#BestTot(G): Given a 2-player game G, outputs the places where the total pay-off is the best, together with the pay-off
BestTot:=proc(G) local a1,a2,S,rec:
S:={[1,1]}:
rec:=G[1][1][1]+G[1][1][2]:

for a1 from 1 to nops(G) do
 for a2 from 1 to nops(G[1]) do
   if G[a1][a2][1]+G[a1][a2][2]>rec then
     S:={[a1,a2]}:
  rec:= G[a1][a2][1]+G[a1][a2][2]:
   elif   G[a1][a2][1]+G[a1][a2][2]=rec then
     S:=S union {[a1,a2]}:
   fi:
 od:
od:
S,rec:

end:

#IsNE(G,a1,a2): Is [a1,a2] a Nash Equilibrium of the game G?
IsNE:=proc(G,a1,a2)
member(a1,BR1(G,a2)) and member(a2,BR2(G,a1)):
end:

#NE(G): Given a 2-player game G given by a bi-matrix G finds the set of Nash Equilibria
NE:=proc(G) local a1,a2,S:
S:={}:

for a1 from 1 to nops(G) do
 for a2 from 1 to nops(G[1]) do
  if IsNE(G,a1,a2) then
   S:=S union {[a1,a2]}:
  fi:
 od:
od:
S:

end:


#RG(a,b,K): A random 2-player (static) game where the Row player has a strategies (Row 1, .., Row a), the Column player has b strategies (Col. 1, ..., Crol. b)
#and the pay-offs are from 0 to K. For example, to see a random game where Player Row has four strategy choices, and Player Column has five strategy choices
#and the pay-offs are integers from 0 to 20 type:
RG:=proc(a,b,K)  local ra,i,j:
ra:=rand(0..K):

[seq([seq([ra(),ra()],j=1..b)],i=1..a)]:

end:

#Start FROM C4.txt

#C4.txt: Maple Code for Lecture 4 of Math640 (Spring 2022)

Help4:=proc(): print(` FP(F), BR12(G), BR21(G) `): end:

#FP(f): Given a discrete function f:= {1,...,n} ->{1,...,n} described by a list of length n, find the set of fixed points
FP:=proc(L) local S,i:
S:={}:
for i from 1 to nops(L) do
 if L[i]=i then
  S:=S union {i}:
 fi:
od:
S:
end:

#BR12(G): Given a game G outputs the discrete function f: A1->A1 such that f(a1)= Best Response of Player 1 to the Best Response of Player 2 to Player 1's strategy a1.
#Try: G:=RandDisGame(G,10,20): BR12(G);
BR12:=proc(G) local L1,L2,a1:
L1:=BR1dv(G):
L2:=BR2dv(G):
[seq(L1[L2[a1]],a1=1..nops(L2))]:
end:


#BR21(G): Given a game G outputs the discrete function g: A2->A2 such that g(a2)= Best Response of Player 2 to the Best Response of Player 1 to Player 2's strategy a2.
#Try: G:=RandDisGame(G,10,20): BR21(G);
BR21:=proc(G) local L1,L2,a2:
L1:=BR1dv(G):
L2:=BR2dv(G):
[seq(L2[L1[a2]],a2=1..nops(L1))]:
end:



#FROM C3.txt

Help3:=proc(): print(`RandDisGame(a,b), BR1d(G,a2), BR2d(G,a1),  BR1dv(G), BR2dv(G), DynRC(G), DynCR(G)  `):end:
with(combinat):

#RandDisGame(a,b): A random game where Player Row has a strategies R1, R2, ..., Ra; and Player Col. has b strategies: C1, C2, ..., Cb and all pay-offs are DISCTINCT
#and consist of the set {1,2,...,ab}. Try:
#RandDisGame(4,6);
RandDisGame:=proc(a,b) local pi1,pi2,i1,j1:
pi1:=randperm(a*b):
pi2:=randperm(a*b):
[seq([seq([pi1[b*i1+j1],pi2[b*i1+j1]],j1=1..b)],i1=0..a-1)]:
end:



#BR1d(G,a2): Inputs a bimatrix of a game G and a member a2 of A2 outputs the UNIQUE member of A1 that consists of the best response  
BR1d:=proc(G,a2) local a1:
max[index]([seq(G[a1][a2][1],a1=1..nops(G))]):
end:



#BR2d(G,a1): Inputs a bimatrix of a game G and a member a1 of A1 outputs the UNIQUE member of A2 that consists of the best response  
BR2d:=proc(G,a1) local a2:
max[index]([seq(G[a1][a2][2],a2=1..nops(G[a1]))]):
end:

#BR1dv(G): Inputs a bimatrix of a game G, outputs The list of size A2 with all the  best responses
#
BR1dv:=proc(G) local a2: [seq(BR1d(G,a2),a2=1..nops(G[1]))]:end:

#BR2dv(G): Inputs a bimatrix of a game G, outputs The list of size A1 with all the  best responses
BR2dv:=proc(G) local a1: [seq(BR2d(G,a1),a1=1..nops(G))]:end:


#DynRC(G): The outcome of the Dynamical version of the game G if Row goes first and Column goes next
DynRC:=proc(G) local L,i,iBest, jBest:
#L is the list of size A1 where L[i] is the best response of Player Column to Strategy i and Row gets G[i][L[i]][1]
L:=BR2dv(G):
iBest:=max[index]([seq(G[i][L[i]][1],i=1..nops(G))]):
jBest:=L[iBest]:
[iBest,jBest], G[iBest][jBest]:
end:


#DynCR(G): The outcome of the Dynamical version of the game G if Column goes first and Row goes next
DynCR:=proc(G) local L,j,iBest, jBest:
#L is the list of size A2 where L[j] is the best response of Player Row to Strategy j and Col gets G[L[j]][j]][2]
L:=BR1dv(G):
jBest:=max[index]([seq(G[L[j]][j][2],j=1..nops(L))]):
iBest:=L[jBest]:
[iBest,jBest], G[iBest][jBest]:
end:

#End C3.txt



#FROM C2.txt

Help2:=proc(): print(` GameDB(), Rand2PlayerGame(a,b,K), IsStictDom(v1,v2), FindR(G), FindC(G), ShrinkGame(G), ReducedGame(G) , MyMaxIndex(L), LoadedCoin(p), BR1(G,a2), BR2(G,a1), BR1v(G), BR2v(G)  `):end:



###PREPARED BEFORE CLASS

#GameDB(): A list of length 5 consisiting of a a "data base"  of famous gamess (and less famous ones): 
#Prisoner's dillema, Boattle of the Sexes,  Matching pennies, Figure 1.1.1. in Gibbons, Figure 1.1.4 in Gibbons
#For example, to see the Prisoner's dillema, type
#GameDB()[1])
GameDB:=proc():
[

[ [ [[-1,-1],[-9,0]], [[0,-9],[-6,-6]]], [`Mum`, `Fink`], [`Mum`, `Fink`]],

[ [ [[2,1],[0,0]], [[0,0],[1,2]]], [`Box`, `Opera`],[`Box`, `Opera`]],

[ [ [[-1,1],[1,-1]], [[1,-1],[-1,1]]], [`Odd`, `Even`], [`Odd`, `Even`]],

[ [ [ [1,0],[1,2],[0,1]], [ [0,3],[0,1],[2,0]]],[`Up`, `Down`], [`Left`, `Middle`, `Right`]],

[[ [ [0,4],[4,0],[5,3]],[ [4,0],[0,4],[5,3]], [ [3,5],[3,5],[6,6]]], [`T`, `M`, `B`],[`L`,`C`,`R`]],
[
[
[[0,0],[-1,1],[1,-1]],
[[1,-1],[0,0],[-1,1]],
[[-1,1],[1,-1],[0,0]]
],
[`Scissors`, `Rock`, `Paper`],
[`Scissors`, `Rock`, `Paper`]
]

]:
end:



PrintGame:=proc(G) local i:
matrix(
[
[" ", op(G[3])],
seq([G[2][i],op(G[1][i])],i=1..nops(G[2]))
]):

end:



#Rand2PlayerGame(a,b,K): A random 2-player (static) game where the Row player has a strategies (Row 1, .., Row a), the Column player has b strategies (Col. 1, ..., Crol. b)
#and the pay-offs are from 0 to K. For example, to see a random game where Player Row has four strategy choices, and Player Column has five strategy choices
#and the pay-offs are integers from 0 to 20 type:
#matrix(Rand2PlayerGame(4,5,20));
Rand2PlayerGame:=proc(a,b,K)  local ra,i,j:
ra:=rand(0..K):

[ [seq([seq([ra(),ra()],j=1..b)],i=1..a)],[seq(i,i=1..a)],[seq(j,j=1..b)]]:

end:

#IsStictDom(v1,v2): Given two lists of numbers is v1[i]<=v2[i] for all i. For example
#IsStrictDom([1,3,2],[2,4,3]); should return true but IsDom([1,3,2],[2,4,1]); should return false
IsStrictDom:=proc(v1,v2) local i:

for i from 1 to nops(v1) do
 if v1[i]>=v2[i] then
  RETURN(false):
 fi:
od:
true:
end:


#FindR(G): Finds a strictly dominated row stategy
FindR:=proc(G) local G1,RowS,i1,i2,j,v1,v2:
G1:=G[1]:
RowS:=G[2]:

for i1 from 1 to nops(RowS) do
for i2 from i1+1 to nops(RowS) do
  v1:=[seq(G1[i1][j][1],j=1..nops(G1[i1]))]:
  v2:=[seq(G1[i2][j][1],j=1..nops(G1[i2]))]:
  if IsStrictDom(v1,v2) then
   RETURN(i1):
  elif IsStrictDom(v2,v1) then
   RETURN(i2):
  fi:
od:
od:
FAIL:
end:


#FindC(G): Finds a strictly dominated row stategy
FindC:=proc(G) local G1,ColS,j1,j2,i,v1,v2:
G1:=G[1]:
ColS:=G[3]:

for j1 from 1 to nops(ColS) do
for j2 from j1+1 to nops(ColS) do
  v1:=[seq(G1[i][j1][2],i=1..nops(G1))]:
  v2:=[seq(G1[i][j2][2],i=1..nops(G1))]:
  if IsStrictDom(v1,v2) then
   RETURN(j1):
  elif IsStrictDom(v2,v1) then
   RETURN(j2):
  fi:
od:
od:
FAIL:
end:


#ShrinkGame(G): Given a game G, tries to shrink it, if it can't it returs it
ShrinkGame:=proc(G) local i,j,G1,RowS,ColS,i1:

G1:=G[1]:
RowS:=G[2];
ColS:=G[3]:

i:=FindR(G):

if i<>FAIL then
 RETURN([
[op(1..i-1,G1),op(i+1..nops(G1),G1)], 
[op(1..i-1,RowS),op(i+1..nops(RowS),RowS)],ColS
]):
fi:

j:=FindC(G):

if j<>FAIL then

 RETURN([
[seq( [op(1.. j-1,G1[i1]),op(j+1..nops(G1[i1]),G1[i1])],i1=1..nops(G1))],RowS, [op(1..j-1,ColS),op(j+1..nops(ColS),ColS)]]):
fi:

FAIL:

end:

#ReducedGame(G): The reduced game of G after all the possible Elimination of dominated  strategy
ReducedGame:=proc(G) local G1,G2:
G1:=G:
G2:=ShrinkGame(G1):

while G2<>FAIL do
G1:=G2:
G2:=ShrinkGame(G1):
od:
G1:
end:

####END PREPARED BEFORE CLASS

#DONE DURING CLASS WITH STUDENTS' HELP

#MyMaxIndex(L): [Suggested by Victoria Chayes]. 
#Inputs a list L of numbers, output the SET of places (indices) where it is maximum. For example
#MyMaxIndex([5,6,7,7,1,7]); should give {3,4,6}
MyMaxIndex:=proc(L) local i,m,S:
m:=max(L):
S:={}:
for i from 1 to nops(L) do
 if L[i]=m then
    S:=S union {i}:
 fi:
od:
S:
end:



#Clarification (Jan. 26, 2022) thanks to Blair Seidler: Note that G is the actual BI-MATRIX, so to experiment with it use, for example
#G:=Rand2PlayerGame(10,20,100)[1] 
#BR1(G,a2): Inputs a bi-matrix of a 2-player game G and a member a2 of A2 outputs the SUBSET of A1 that consists of the best response  
BR1:=proc(G,a2) local a1:
MyMaxIndex([seq(G[a1][a2][1],a1=1..nops(G))]):
end:



#BR2(G,a1): Inputs a bimatrix of a game G and a member a1 of A1 outputs the SUBSET of A2 that consists of the best response  
BR2:=proc(G,a1) local a2:
MyMaxIndex([seq(G[a1][a2][2],a2=1..nops(G[a1]))]):
end:

#BR1v(G): Inputs a bimatrix of a game G, outputs The list of size A2 with all the set of best responses
#
BR1v:=proc(G) local a2: [seq(BR1(G,a2),a2=1..nops(G[1]))]:end:

#BR2v(G): Inputs a bimatrix of a game G, outputs The list of size A1 with all the set of best responses
BR2v:=proc(G) local a1: [seq(BR2(G,a1),a1=1..nops(G))]:end:

#end DURING CLASS WITH STUDENTS' HELP

#Dne after class

#LoadedCoin(p): outputs 1 or 2 with probability p and 1-p respecively. Try:
#[seq(LoadedCon(1/3)),i=1..300)];
LoadedCoin:=proc(p) local m,n,ra:
m:=numer(p):
n:=denom(p):

ra:=rand(1..n)():

if ra<=m then
 1:
else
 2:
fi:

end:

#End FROM C2.txt