#OK to post homework
#Blair Seidler, 1/30/22, Assignment 3
with(combinat):
with(ListTools):

Help:=proc(): print(`AllGames(a,b), FullStat(a,b), SampleStat(a,b,K)`): end:

#2. Done

#3.
(*
      G := [[[13, 8], [14, 6], [11, 4], [6, 5]], 

            [[9, 10], [1, 11], [16, 3], [7, 9]], 

            [[3, 2], [2, 7], [12, 15], [5, 1]], 

            [[4, 13], [10, 16], [8, 12], [15, 14]]]

Row goes first:

Row selects R1, COL vector [8,6,4,5], Col selects C1, outcome [1,1] result [13,8]
Row selects R2, COL vector [10,11,3,9], Col selects C2, outcome [2,2] result [1,11]
Row selects R3, COL vector [2,7,15,1], Col selects C3, outcome [3,3] result [12,15]
Row selects R4, COL vector [13,16,12,14], Col selects C2, outcome [4,2] result [10,16]

Best option of row is R1 result, so outcome of game is [1,1] with payout [13,8]

Col goes first:

Col selects C1, ROW vector [13,9,3,4], Row selects R1, outcome [1,1] result [13,8]
Col selects C2, ROW vector [14,1,2,10], Row selects R1, outcome [1,2] result [14,6]
Col selects C3, ROW vector [11,16,12,8], Row selects R2, outcome [2,3] result [16,3]
Col selects C4, ROW vector [6,7,5,15], Row selects R4, outcome [4,4] result [15,14]

Best option of col is C4 result, so outcome of game is [4,4] with payout [15,14]

Pure Nash Equilibrium: [1,1] result [13,8]

Confirmed:
> DynRC(G);
                        [1, 1], [13, 8]
> DynCR(G);
                        [4, 4], [15, 14]
> PureNashEqui(G);
                            {[1, 1]}

*)

(*
   G := [[[30, 28], [16, 13], [15, 21], [11, 2],  [17, 23], [29, 15], [28, 33]],
         [[32, 27], [1, 35],  [26, 17], [31, 16], [12, 11], [9, 19],  [22, 1]],
         [[18, 6],  [13, 9],  [24, 29], [8, 26],  [25, 32], [33, 34], [10, 10]],
         [[34, 30], [27, 25], [2, 12],  [23, 20], [20, 22], [4, 8],   [14, 24]],
         [[21, 3],  [6, 5],   [3, 4],   [19, 31], [5, 7],   [7, 14],  [35, 18]]]

Row goes first:

Row selects R1, COL vector [28,13,21,2,23,15,33], Col selects C7, outcome [1,7] result [28,33]
Row selects R2, COL vector [27,35,17,16,11,19,1], Col selects C2, outcome [2,2] result [1,35]
Row selects R3, COL vector [6,9,29,26,32,34,10], Col selects C6, outcome [3,6] result [33,34]
Row selects R4, COL vector [30,25,12,20,22,8,24], Col selects C1, outcome [4,1] result [34,30]
Row selects R5, COL vector [3,5,4,31,7,14,18], Col selects C4, outcome [5,4] result [19,31]

Best option of row is R4 result, so outcome of game is [4,1] with payout [34,30]

Col goes first:

Col selects C1, ROW vector [30,32,18,34,21], Row selects R4, outcome [4,1] result [34,30]
Col selects C2, ROW vector [16,1,13,27,25], Row selects R4, outcome [4,2] result [27,25]
Col selects C3, ROW vector [15,26,24,2,3], Row selects R2, outcome [2,3] result [26,17]
Col selects C4, ROW vector [11,31,8,23,19], Row selects R2, outcome [2,4] result [31,16]
Col selects C5, ROW vector [17,12,25,20,5], Row selects R3, outcome [3,5] result [25,32]
Col selects C6, ROW vector [29,9,33,4,7], Row selects R3, outcome [3,6] result [33,34]
Col selects C7, ROW vector [28,22,10,14,35], Row selects R5, outcome [5,7] result [35,18]

Best option of col is C6 result, so outcome of game is [3,6] with payout [33,34]

Pure Nash Equilibrium: [3,6] result [33,34] AND [4,1] result [34,30]

Confirmed:
>DynRC(G);
                        [4, 1], [34, 30]
>DynCR(G);
                        [3, 6], [33, 34]
>PureNashEqui(G);
                        {[3, 6], [4, 1]}

*)

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

#5.
#FullStat(a,b): inputs positive integes a and b (not too big!) and outputs the following list
#[Number of Games with No Nash Equilibria, Number of Games with One Nash Equilibria,
# Number of Games with Two Nash Equilibria,
# Number of Games where the Row Player was better off playing second,
# Number of Games where the Column Player was better off playing second,
# Number of Games where DynRC(G)=DynCR(G)]
FullStat:=proc(a,b) local AG,G,nash,res,r2bet,c2bet,rceq:
AG:=AllGames(a,b):
nash:=[seq(nops(PureNashEqui(G)),G in AG)]:
res:=[seq([DynRC(G)[2],DynCR(G)[2]],G in AG)]:
r2bet:=[seq(evalb(i[2,1]>i[1,1]),i in res)]:
c2bet:=[seq(evalb(i[1,2]>i[2,2]),i in res)]:
rceq:=[seq(evalb(i[1]=i[2]),i in res)]:
[Occurrences(0,nash),Occurrences(1,nash),Occurrences(2,nash),
Occurrences(true,r2bet),Occurrences(true,c2bet),Occurrences(true,rceq)]:
end:

(* Results:
> FullStat(2,2);
          [72, 432, 72, 72, 72, 420]

> FullStat(3, 2);
          [86400, 345600, 86400, 86760, 83952, 331560]

*)

#6.
#SampleStat(a,b,K): inputs positive integers a and b, and a fairly large positive integer K,
#and generates K random games where Row has a strategies and Col has b strategies
#evalf( [(Number of Games with No Nash Equilibria)/K, (Number of Games with One Nash Equilibria)/K,
#(Number of Games with Two Nash Equilibria)/K,
#(Number of Games where the Row Player was better off playing second)/K,
#(Number of Games where the Column Player was better off playing second)/K,
#(Number of Games where DynRC(G)=DynCR(G))/K]):
SampleStat:=proc(a,b,K) local AG,G,nash,res,r2bet,c2bet,rceq:
AG:=[seq(RandDisGame(a,b),i in 1..K)]:
nash:=[seq(nops(PureNashEqui(G)),G in AG)]:
res:=[seq([DynRC(G)[2],DynCR(G)[2]],G in AG)]:
r2bet:=[seq(evalb(i[2,1]>i[1,1]),i in res)]:
c2bet:=[seq(evalb(i[1,2]>i[2,2]),i in res)]:
rceq:=[seq(evalb(i[1]=i[2]),i in res)]:
evalf([Occurrences(0,nash),Occurrences(1,nash),Occurrences(2,nash),
Occurrences(true,r2bet),Occurrences(true,c2bet),Occurrences(true,rceq)]/K):
end:

(* Results:
> for i to 10 do SampleStat(10, 10, 10000); od;
   [0.3225000000, 0.4148000000, 0.2027000000, 0.2887000000, 0.2916000000, 0.3384000000]
   [0.3219000000, 0.4154000000, 0.2024000000, 0.2841000000, 0.2879000000, 0.3447000000]
   [0.3283000000, 0.4132000000, 0.1967000000, 0.2887000000, 0.2957000000, 0.3474000000]
   [0.3272000000, 0.4103000000, 0.2011000000, 0.2962000000, 0.2913000000, 0.3377000000]
   [0.3323000000, 0.4117000000, 0.1996000000, 0.2995000000, 0.3010000000, 0.3368000000]
   [0.3243000000, 0.4124000000, 0.2081000000, 0.2962000000, 0.2916000000, 0.3369000000]
   [0.3316000000, 0.4082000000, 0.2014000000, 0.2954000000, 0.2964000000, 0.3437000000]
   [0.3150000000, 0.4238000000, 0.2019000000, 0.2884000000, 0.2817000000, 0.3494000000]
   [0.3268000000, 0.4061000000, 0.2075000000, 0.2927000000, 0.2928000000, 0.3389000000]
   [0.3255000000, 0.4110000000, 0.2054000000, 0.2925000000, 0.2868000000, 0.3432000000]

Yes, these results all look similar.
*)

#### Code included from hw2 ####
IsNash:= proc(G,a1,a2): evalb((a2 in BR2(G,a1)) and (a1 in BR1(G,a2))): end:

PureNashEqui:=proc(G) local row,col,nash:
nash:={}:
for row from 1 to nops(G) do
  for col from 1 to nops(G[1]) do
    if IsNash(G,row,col) then
      nash:=nash union {[row,col]}:
    fi:
  od:
od:
nash:
end:

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



#BR1(G,a2): Inputs a 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 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[1]))]):
end:

#BR1v(G): 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:=proc(G) local a1: [seq(BR2(G,a1),a1=1..nops(G))]:end:


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

#### Code included from C3.txt ####

#C3.txt: Maple Code for Lecture 3 of Math640 (Spring 2022)

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


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