# OK to post homework
# Vikrant, Jan 30 2022, Assignment 3

# ================================================================================
# 0. Code that has been given.
# ================================================================================

with(combinat):

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

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

BR2d:=
    proc(G,a1)
        local a2:
        max[index]([seq(G[a1][a2][2],a2=1..nops(G[a1]))]):
    end:

BR1dv:=
    proc(G)
        local a2:
        [seq(BR1d(G,a2),a2=1..nops(G[1]))]:
    end:

BR2dv:=
    proc(G)
        local a1:
        [seq(BR2d(G,a1),a1=1..nops(G))]:
    end:

DynRC:=
    proc(G)
        local L,i,iBest,jBest:
        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:=
    proc(G)
        local L,j,iBest,jBest:
        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:

# ================================================================================
# 1. Samples from Garvan's Maple Primer.
# ================================================================================

(*

plot(cosh(x),x=-Pi..Pi);
plot([[cos(t),cosh(t),t=0..2*Pi],[sin(t)/2,t,t=0..cosh(2*Pi)]]);
# Looks like this package is deprecated.
with(linalg):
factor(det(matrix(3,3,(i,j)->x^(j*i))));

*)

# ================================================================================
# 2. Read and understand C3.txt.
# ================================================================================

Done.

# ================================================================================
# 3. RandDisGame(4,4) and (5,7) data.
# ================================================================================

randomize(1879406939):

# [[[4, 4], [8, 6], [15, 16], [6, 3]], [[9, 13], [1, 2], [7, 9], [12, 1]], [[10, 7], [16, 8], [13, 11], [3, 15]], [[11, 10], [14, 5], [5, 12], [2, 14]]]
G1:= RandDisGame(4,4):
(*
# 1->3->[15,16]
# 2->1->[9,13]
# 3->4->[3,15]
# 4->4->[2,14]
RowFirstByHand:= [1,3],[15,16]:
evalb(RowFirstByHand = DynRC(G1));
# 1->4->[11,10]
# 2->3->[16,8]
# 3->1->[15,16]
# 4->2->[12,1]
ColFirstByHand:= [1,3],[15,16]:
evalb(ColFirstByHand = DynCR(G1));
NashEquiByHand:= {[1,3]}:
evalb(NashEquiByHand = PureNashEqui(G1));
*)


# Omitted due to length.
G2:= RandDisGame(5,7):
(*
# 1->3->[27,35]
# 2->6->[34,33]
# 3->3->[9,34]
# 4->1->[18,26]
# 5->4->[26,30]
RowFirstByHand:= [2,6],[34,33]:
evalb(RowFirstByHand = DynRC(G2));
# 1->3->[24,31]
# 2->5->[17,27]
# 3->5->[29,8]
# 4->4->[35,11]
# 5->3->[31,2]
# 6->2->[34,33]
# 7->3->[32,31]
ColFirstByHand:=[2,6],[34,33]:
evalb(ColFirstByHand = DynCR(G2));
NashEquiByHand:= {[2,6]}:
evalb(NashEquiByHand = PureNashEqui(G2));
*)

# ================================================================================
# 3. All games.
# ================================================================================

AllGames:=
    proc(a,b)
        local S,pi1,pi2:
        S:= permute(a*b):
        {seq(seq(MakeGame(pi1,pi2,a,b),pi2 in S),pi1 in S)}: 
    end:

MakeGame:=
    proc(pi1,pi2,a,b)
        local i,j:
        [seq([seq([pi1[b*i+j],pi2[b*i+j]],j=1..b)],i=0..a-1)]:
    end:

# ================================================================================
# 4. Full statistics.
# ================================================================================

FullStat:=
    proc(a,b)
        local G:
        local Gs:= AllGames(a,b):
        [
            nops([seq(`if`(nops(PureNashEqui(G)) = 0,1,NULL),G in Gs)]),
            nops([seq(`if`(nops(PureNashEqui(G)) = 1,1,NULL),G in Gs)]),
            nops([seq(`if`(nops(PureNashEqui(G)) = 2,1,NULL),G in Gs)]),
            nops([seq(`if`(DynCR(G)[2][1] > DynRC(G)[2][1],1,NULL),G in Gs)]),
            nops([seq(`if`(DynRC(G)[2][2] > DynCR(G)[2][2],1,NULL),G in Gs)]),
            nops([seq(`if`(DynRC(G) = DynCR(G),1,NULL),G in Gs)])
        ]:
    end:

# FullStat(2,2);
# [72, 432, 72, 72, 72, 420]

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

# ================================================================================
# 5. Sample statistics.
# ================================================================================

SampleStat:=
    proc(a,b,K)
        local i,G:
        local Gs:= {seq(RandDisGame(a,b),i=1..K)}:
        evalf(
            [
                nops([seq(`if`(nops(PureNashEqui(G)) = 0,1,NULL),G in Gs)])/K,
                nops([seq(`if`(nops(PureNashEqui(G)) = 1,1,NULL),G in Gs)])/K,
                nops([seq(`if`(nops(PureNashEqui(G)) = 2,1,NULL),G in Gs)])/K,
                nops([seq(`if`(DynCR(G)[2][1] > DynRC(G)[2][1],1,NULL),G in Gs)])/K,
                nops([seq(`if`(DynRC(G)[2][2] > DynCR(G)[2][2],1,NULL),G in Gs)])/K,
                nops([seq(`if`(DynRC(G) = DynCR(G),1,NULL),G in Gs)])/K
            ]
        ):
    end:

(*

for i from 1 to 5 do:
    print(SampleStat(10,10,10000));
od:

# Experiment 1
# [0.3270000000, 0.4189000000, 0.1979000000, 0.2963000000, 0.2953000000, 0.3414000000]
# Experiment 2
# [0.3319000000, 0.4047000000, 0.2035000000, 0.2929000000, 0.2939000000, 0.3410000000]
# Experiment 3
# [0.3280000000, 0.4127000000, 0.2011000000, 0.2916000000, 0.2968000000, 0.3446000000]
# Experiment 4
# [0.3331000000, 0.4075000000, 0.2044000000, 0.2993000000, 0.2947000000, 0.3435000000]
# Experiment 5
# [0.3259000000, 0.4120000000, 0.2040000000, 0.2891000000, 0.2909000000, 0.3444000000]

*)

# ================================================================================
# oo. Helpers from hw2.
# ================================================================================

IsNash:=
    proc(G,a1,a2)
        a1 = BR1d(G,a2) and a2 = BR2d(G,a1):
    end:

PureNashEqui:=
    proc(G)
        local i,j:
        {seq(seq(`if`(IsNash(G,i,j),[i,j],NULL),j=1..nops(G[i])),i=1..nops(G))}:
    end: