#Algorithm:
#Each person starts out with no people “canceiled” from his or her list. People will be
#cancelled from lists as the algorithm progresses.
#For each man m, do propose(m), as defined below.
#propose(m):
#Let W be the first uncancelled woman on m’s preference list.
#Do: refuse(PV,m), as defined below.
#refuse(W,m):
#Let m’ be W’s current partner (if any).
#If W prefers m’ to m, then she rejects m, in which case:
#(1) cancel m off W’s list and W off m’s list;
#(2) do: propose(m). (Now m must propose to someone else.)
#Otherwise, W accepts m as her new partner, in which case:
#(1) cancel m’ off W’s list and W off m”s list;
#(2) do: propose( (Now m’ must propose to someone else.)

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:

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


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