#OK to post homework 
#Natalya Ter-Saakov, Feb 27, Assignment 11
read "C11.txt":
read "hw10NatalyaTer-Saakov.txt":

#Problem 1:

#The Gale-Shapley algorithm: Verbose version
GaleSh:=proc(M,W) local n,co,Mrank,Wrank,i,j,CurW,CurH, SinM, SinW,i1:
n:=nops(M):
if not (type(M,list) and type(W,list) and nops(W)=n) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
for i from 1 to n do 
if not (type(M[i], list) and nops(M[i])=n and nops(convert(M[i],set))=n and type(W[i], list) and nops(W[i])=n and nops(convert(W[i],set))=n) then  
print(`Bad input`):
 RETURN(FAIL):
fi:
od:

co:=0:
#Mrank[i]: Ranking of Mr. i of CURRENTLY AVAILABLE WOMEN
#Wrank[j]:=Ranking of Ms. j of CURRENTLY AVAILABLE MEN
for i from 1 to n do
 Mrank[i]:=inv(M[i]):
od:
for j from 1 to n do
 Wrank[j]:=inv(W[j]):
od:
#CurW[i]:=Current wife of Mr. i (0 if he is single)
#CurH[j]:=Current husband of Ms. j (0 if she is single)
for i from 1 to n do
 CurW[i]:=0:
od:
for j from 1 to n do
 CurH[j]:=0:
od:

SinM:={seq(i,i=1..n)}:
SinW:={seq(j,j=1..n)}:

while SinM<>{} do
#i is the current proposer
i:=SinM[1]:
#Mr. i proposes to the best remaining woman
j:=Mrank[i][1]:
#Mr. i1 is the current husband of Ms. j
i1:=CurH[j]:
co:=co+1:
print(`I am Mr.`, i, `Ms. `, j, `would you marry me?`):
if i1=0 then
 print(`Sure, you are my`, W[j][i], `-th choice, but since I am single it does not matter, so sure`):
SinM:=SinM minus {i}:
SinW:=SinW minus {j}:
CurW[i]:=j:
CurH[j]:=i:

   elif W[j][i1]<W[j][i] then
    print(`How dare you!, I am currently married to Mr.`, i1, `whom I rank `, W[j][i1], `-th best `):
    print(`while you rank only`, W[j][i], `-th in my list, so go to to hell`):
    Mrank[i]:=[op(2..nops(Mrank[i]),Mrank[i])]:

   else

   print(`Sure, I am currently married to Mr.`, i1, `whom I rank `, W[j][i1], `-th best `):
    print(`while you rank even better`, W[j][i], `-th in my list, so yes! I am all yours`):

   SinM:=SinM minus {i} union {i1}:
   
    Mrank[i1]:=[op(2..nops(Mrank[i1]),Mrank[i1])]:
   CurW[i]:=j:
    CurH[j]:=i:
   CurW[i1]:=0:
fi:
od:

print(`The stable matching gotten from Gale Shapley is as follows`):
for i from 1 to n do
 print(`Mr.`, i , `is married to Ms. `, CurW[i]):
od:


print(`It took`, co, `proposals `):
[seq(CurW[i],i=1..n)],co:

end:

#GaleShT(M,W):terse version of GaleSh(M,W):
GaleShT:=proc(M,W) local n,co,Mrank,Wrank,i,j,CurW,CurH, SinM, SinW,i1:
n:=nops(M):
if not (type(M,list) and type(W,list) and nops(W)=n) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
for i from 1 to n do 
if not (type(M[i], list) and nops(M[i])=n and nops(convert(M[i],set))=n and type(W[i], list) and nops(W[i])=n and nops(convert(W[i],set))=n) then 
print(`Bad input`):
 RETURN(FAIL):
fi:
od:
co:=0:
#Mrank[i]: Ranking of Mr. i of CURRENTLY AVAILABLE WOMEN
#Wrank[j]:=Ranking of Ms. j of CURRENTLY AVAILABLE MEN
for i from 1 to n do
 Mrank[i]:=inv(M[i]):
od:
for j from 1 to n do
 Wrank[j]:=inv(W[j]):
od:
#CurW[i]:=Current wife of Mr. i (0 if he is single)
#CurH[j]:=Current husband of Ms. j (0 if she is single)
for i from 1 to n do
 CurW[i]:=0:
od:
for j from 1 to n do
 CurH[j]:=0:
od:

SinM:={seq(i,i=1..n)}:
SinW:={seq(j,j=1..n)}:

while SinM<>{} do
#i is the current proposer
i:=SinM[1]:
#Mr. i proposes to the best remaining woman
j:=Mrank[i][1]:
#Mr. i1 is the current husband of Ms. j
i1:=CurH[j]:
co:=co+1:

if i1=0 then

SinM:=SinM minus {i}:
SinW:=SinW minus {j}:
CurW[i]:=j:
CurH[j]:=i:

   elif W[j][i1]<W[j][i] then
   
    Mrank[i]:=[op(2..nops(Mrank[i]),Mrank[i])]:

   else

  
   SinM:=SinM minus {i} union {i1}:
   
    Mrank[i1]:=[op(2..nops(Mrank[i1]),Mrank[i1])]:
   CurW[i]:=j:
    CurH[j]:=i:
   CurW[i1]:=0:
fi:
od:
[seq(CurW[i],i=1..n)],co:

end:

#Problem 2: women proposing
GaleShTw:=proc(M,W) local temp:
temp:=GaleShT(W,M):
inv(temp[1]), temp[2]:
end:

#Problem 3
#CheckOpt(M,W): inputs M and W, finds the matchings where the men do the proposing and the one where the women do it,
# and verify that each man got a better (or equal) wife (in his eyes) and each woman get a worse (or equal) husband in the first case.
CheckOpt:=proc(M,W) local HMprop, HWprop, WMprop, WWprop, i:
HMprop:=GaleShT[1]:
HWprop:=GaleShTw[1]:
WMprop:=inv(HMprop):
WWprop:=inv(HWprop):
for i from 1 to nops(Mprop) do
if M[i][MMprop[i]]>M[i][MWprop[i]] then
print("Mr.",i,"prefers the wife he would have had if the women proposed"):
return(false):
fi:
if W[i][WMprop[i]]<W[i][WWprop[i]] then
print("Ms.",i,"prefers the husband she would have had if the men proposed"):
return(false):
fi:
od:
end:

#Problem 4
#CompGaleSh(M,W):inputs partial preferences and outputs a matching or FAILS if no complete matching exists
CompGaleSh:=proc(M,W) local n, S, nM, nW, res,i:
n:=nops(M):
if not (type(M,list) and type(W,list) and nops(W)=n) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
S:={seq(i,1..n+1)}:
for i from 1 to n do
nM[i]:=[ops(M[i]), n+1, ops(S-convert(M[i],set))]:
nW[i]:=[ops(W[i]), n+1, ops(S-convert(W[i],set))]:
od:
nM[n+1]:=[seq(i,1..n+1)]:
nW[n+1]:=[seq(i,1..n+1)]:
res:=GaleShT(nM,nW):
if res[n+1]<>n+1 then return(FAIL) fi:
[ops(GaleShT[1], 1..n)], GaleShT[2]:
end:

#Problem 5
Compare:=proc(n) local K, M,W, fails, stupid, notStupid, st:
stupid:=0: notStupid:=0:fails:=0:
for K from 1 to 1000 do
M,W:=RT(n):
st:=FindStable(M,W,10000):
if st=FAIL then fails+=1:
elif st[2]<GaleShT(M,W)[2] then stupid+=1: 
else notStupid+=1: fi:
od:
evalf([n,fails/1000,stupid/1000,notStupid/1000]):
end:

#for n from 1 to 20 do Compare(n); od;
#                        [1., 0., 0., 1.]
#              [2., 0., 0.8000000000, 0.2000000000]
#             [3., 0., 0.9080000000, 0.09200000000]
#              [4., 0., 0.8950000000, 0.1050000000]
#        [5., 0.002000000000, 0.8730000000, 0.1250000000]
#        [6., 0.007000000000, 0.8410000000, 0.1520000000]
#        [7., 0.02000000000, 0.7790000000, 0.2010000000]
#        [8., 0.02700000000, 0.7360000000, 0.2370000000]
#        [9., 0.04000000000, 0.6810000000, 0.2790000000]
#        [10., 0.06800000000, 0.5970000000, 0.3350000000]
#        [11., 0.08600000000, 0.5170000000, 0.3970000000]
#        [12., 0.1160000000, 0.4590000000, 0.4250000000]
#       [13., 0.1850000000, 0.3430000000, 0.4720000000]
#       [14., 0.1940000000, 0.2950000000, 0.5110000000]
#       [15., 0.2670000000, 0.2410000000, 0.4920000000]
#       [16., 0.3200000000, 0.1900000000, 0.4900000000]
#       [17., 0.3610000000, 0.1320000000, 0.5070000000]
#       [18., 0.4170000000, 0.09900000000, 0.4840000000]
#       [19., 0.5230000000, 0.07800000000, 0.3990000000]
#       [20., 0.5810000000, 0.03900000000, 0.3800000000]

# For small n (until about 10), the stupid way is faster most of the time. 
# After this, around n=18, the stupid way starts failing more than half of the time.

#Problem 6
#EstStat(n,K): inputs a positive integer n , and a large positive integer K (at least 1000),
# generates K exaples of M,W (using RT(n)) and outputs the average number of proposals, followed by the standard-deviation.
EstStat:=proc(n,K) local i, avg, st,co, M,W:
avg:=0: st:=0:
for i from 1 to K do
M,W:=RT(n):
co:=GaleShT(M,W)[2]:
avg+=co: st+=co^2: 
od:
evalf([avg/K, sqrt(st/K-(avg/K)^2)]):
end:

#for i from 1 to 15 do EstStat(10*i,10000); od;
#                   [23.85580000, 6.089072702]
#                   [62.47710000, 15.75528088]
#                   [107.0135000, 27.14670362]
#                   [156.6369000, 38.66851766]
#                   [208.7178000, 50.46475764]
#                   [262.3259000, 62.36792035]
#                   [317.9663000, 75.03471839]
#                   [376.0865000, 87.53446760]
#                   [435.9536000, 100.0851450]
#                   [495.9903000, 112.2017416]
#                   [557.6206000, 125.1234776]
#                   [620.2300000, 135.6476874]
#                   [683.3427000, 151.7878205]
#                   [747.2442000, 163.5359189]
#                   [813.6449000, 177.0906288]

#Problem 7
#FindNeighbors(pi): inputs a permutation pi, and outputs the n-1 permutations formed by swapping two adjacent entries
FindNeighbors:=proc(pi) local neighbors, n,i:
neighbors:={}:
n:=nops(pi):
for i from 1 to n-1 do
neighbors:=neighbors union {[op(1..i-1, pi), pi[i+1],pi[i], op(i+2..n,pi)]}:
od:
neighbors:
end:

MayBeNotSoStupidFindStable:=proc(M,W) local n, min, pi, altPi, best, t, co:
n:=nops(M): co:=0:
pi:=randperm(n):
best:=pi:
min:=CountInstabilities(M,W,pi):
while CountInstabilities(M,W,pi)<>0 do
co+=1:
for altPi in FindNeighbors(pi) do
t:=CountInstabilities(M,W,altPi):
if t<min then 
min:=t: 
best:=altPi:
fi:
od:
if best = pi then return(FAIL): fi:
pi:=best:
od:
pi,co:
end:

Compare2:=proc(n) local K, M,W, fails, semiStupid, notStupid, st:
semiStupid:=0: notStupid:=0:fails:=0:
for K from 1 to 1000 do
M,W:=RT(n):
st:=MayBeNotSoStupidFindStable(M,W):
if st=FAIL then fails+=1:
elif st[2]<GaleShT(M,W)[2] then semiStupid+=1: 
else notStupid+=1: fi:
od:
evalf([n,fails/1000,semiStupid/1000,notStupid/1000]):
end:

for n from 1 to 15 do Compare2(n); od;
#                        [1., 0., 1., 0.]
#                        [2., 0., 1., 0.]
#        [3., 0.2940000000, 0.7030000000, 0.003000000000]
#        [4., 0.6140000000, 0.3850000000, 0.001000000000]
#              [5., 0.8240000000, 0.1760000000, 0.]
#             [6., 0.9390000000, 0.06100000000, 0.]
#             [7., 0.9830000000, 0.01700000000, 0.]
#             [8., 0.9930000000, 0.007000000000, 0.]
#             [9., 0.9950000000, 0.005000000000, 0.]
#            [10., 0.9990000000, 0.001000000000, 0.]
#                       [11., 1., 0., 0.]
#                       [12., 1., 0., 0.]
#                       [13., 1., 0., 0.]
#                       [14., 1., 0., 0.]
#                       [15., 1., 0., 0.]

#Conclusion: When this works, it beats Gale-Shapley almost always.
# However, it rarely works.