#OK to post homework
#Yuxuan Yang, Feb 24th, Assignment 11
with(combinat):



#2
GaleShTw:=proc(M,W) 
inv(GaleShT(W,M)[1]),GaleShT(W,M)[2]:
end:

#3

CheckOpt:=proc(M,W) local pi1,pi2,i:
pi1:=GaleShT(M,W)[1]:
pi2:=inv(GaleShT(W,M)[1]):
for i from 1 to nops(pi1) do
if M[i][pi1[i]]>M[i][pi2[i]] then
 RETURN(FAIL):
fi:
if W[i][inv(pi1)[i]]<W[i][inv(pi2)[i]] then
 RETURN(FAIL):
fi:
od:
RETURN(TRUE):
end:

#4
#I have read and understood the algorithm that introduces addtional "dummy" man and woman.
#but I have no ideas about the input W,M for this case. 

#5
#a, b := RT(50);st := time();FindStable(a, b, 10000);time() - st;st2 := time();GaleShT(a, b);time() - st2;
#    6.391                0.
#the first one is always worse. Maybe one reason is it involves "print".

#6
EstStat:=proc(n,K) local f,x,M,W,i,ex,via:
f:=0:
for i from 1 to K do
M,W:=RT(n):
f:=f+x^(GaleShT(M,W)[2]):
od:
f:=f/K:
ex:=subs(x=1,diff(f,x)):
via:=subs(x=1,diff(diff(f,x),x))-ex^2:
evalf([ex,sqrt(via)]):
end:

#[23.96740000, 3.757464203], [62.33150000, 13.80445970], [107.2785000, 24.50480846], [156.7089000, 36.63427440], 
#[208.0608000, 48.40287908], [262.1291000, 60.87098926], [318.9593000, 74.14335941], [376.4519000, 86.15244504], 
#[435.5173000, 98.03289346], [497.6458000, 108.7833652], [555.9758000, 121.8682420], [619.2498000, 135.1021066], 
#[684.4617000, 146.9452001], [745.8299000, 159.1502116], [812.0498000, 174.4173418]
#Expectation is about degree two
#standard-deviation is about degree one


#7
MayBeNotSoStupidFindStable:=proc(M,W) local j,i,pi,pi1,n,ins,ins1,pi2,val:
n:=nops(M[1]):
pi:=randperm(n):
ins:=CountInstablities(M,W,pi):
for j from 1 to 10000 do
if ins=0 then RETURN([pi,j]): fi:
val:=1:
 for i from 1 to n-1 do
 pi1:=[op(1..i-1,pi),pi[i+1],pi[i],op(i+2..n,pi)]:
 print(pi1):
 ins1:=CountInstablities(M,W,pi1):
 if ins1<ins then
 ins:=ins1:
 val:=0:
 pi2:=pi1:
 fi:
 od:
if val=1 then RETURN(FAIL): fi:
pi:=pi2:
od:
end:
#FAIL rate is too high to be a good algorithm.


Help11:=proc(): print(` inv(pi), GaleSh(M,W), GaleShT(M,W) `): end:

with(combinat):

#inv(pi): The inverse of the permutation pi
inv:=proc(pi) local n,i,T:
n:=nops(pi):
for i from 1 to n do
 T[pi[i]]:=i:
od:
[seq(T[i],i=1..n)  ]:
end:

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


#Yuxuan Yang, Feb 24th, Assignment 10
with(combinat):

CountInstablities:=proc(M,W,pi) local n,i1,i2,num:
num:=0:
n:=nops(pi):
for i1 from 1 to n do
 for i2 from 1 to n do
  if i1<>i2 then
   if not IsStable1(M,W,pi,i1,i2) then
      num:=num+1:
   fi:
 fi:
od:
od:
num:
end:

GFstab:=proc(n,K,x) local i,pi,f,rt:
f:=0:
for i from 1 to K do
 pi:=randperm(n):
 rt:=RT(n):
 f:=f+x^(CountInstablities(rt[1],rt[2],pi)):
od:
f:
end:



#c10.txt
Help10:=proc():print(` RT(n), IsStable1(M,W,pi,i1,i2), IsStable(M,W,pi), FindStable(M,W,K) `):end:
with(combinat):

#RT(n): Generates a pair of random preferance M,W, two lists-of-lists
#M[i][j]:= The ranking of Ms. j in the eyes of Mr. i
#W[j][i]:= The ranking of Mr. i  in the eyes of Ms. j

RT:=proc(n) local i,j:
[seq(randperm(n),i=1..n)],[seq(randperm(n),j=1..n)]:
end:

# i1 is married to j1
#i2 is married to j2
#i1 would rather be married to j2 and j2 would rather be married to i1
IsStable1:=proc(M,W,pi,i1,i2) local j1,j2:
j1:=pi[i1]:
j2:=pi[i2]:

  not(M[i1][j2]<M[i1][j1] and W[j2][i1]<W[j2][i2]) :
  
end:


#IsStable(M,W,pi): inputs preferenace lists-of-lists M,W, and a proposed matching is it stable?
#If it is not it outputs false, and an unstable pair i1,i2:
#
IsStable:=proc(M,W,pi) local n,i1,i2:

n:=nops(pi):

for i1 from 1 to n do
 for i2 from 1 to n do
  if i1<>i2 then
   if not IsStable1(M,W,pi,i1,i2) then
      RETURN(false, [i1,i2]):
   fi:
 fi:
od:
od:
true,true:
end:

#FindStable(M,W,K): Given the Men's and Women's preferance tables, starts with a random matching
#and keeps "fixing it" until hopefully you get a stable matching, we limit the number of
#iterations to K, if none found it will return FAIL. Otherwise it will return a stable
#mathing, followed by the number if iterations it took
FindStable:=proc(M,W,K) local n,pi,co,a,i1,i2,j1,j2:
print(`Warning: This is a stupid program that usually FAILS, for larger n`):
n:=nops(M):
pi:=randperm(n):

for co from 1 to K do
 a:=IsStable(M,W,pi):
 if a[1] then
   RETURN(pi,co):
 else
 i1:=a[2][1]:
 i2:=a[2][2]:
  j1:=pi[i1]:
  j2:=pi[i2]:
 if i1<i2 then
  pi:=[op(1..i1-1,pi),j2,op(i1+1..i2-1,pi),j1,op(i2+1..n,pi)]:
else
pi:=[op(1..i2-1,pi),j1,op(i2+1..i1-1,pi),j2,op(i1+1..n,pi)]:
fi:
fi:
od:
FAIL:
end: