######################################################################
## FindingMatchings.txt Save this file as  FindingMatchings.txt  to  #
# use it, stay in the                                                #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  FindingMatchings.txt<Enter>                 #
## Then follow the instructions given there                          #
##                                                                   #
## Written by AJ Bu, Kayla Gibson, Lucy Martinez,  Natalya           #
##  Ter-Saakov, and Doron Zeilberger, Rutgers University             # 
######################################################################


print(`First Written : May 2022 : tested for Maple 2021 `):
print(`Version: May 2022  `):
print():
print(`This is FindingMatchings.txt, A Maple package`):
print(`accompanying AJ Bu, Kayla Gibson, Lucy Martinez, and Natalya Ter-Saakov's  article: `):
print(`Analysis of the Gale-Shapley Algorithm for the Stable Marriage Problem `):

print():
print(`The most current version is available on WWW at:`):
print(` INSERT WEBSITE HERE .`):
print(`Please report all bugs to: ab1854 at math dot rutgers dot edu .`):
print(``):

print(`----------------------------------------------------`):

print():
print(`To see the major procedures in this package, type "Help()"`):

print():
print(`For  a list of the SUPPORTING functions type "Help1();" `):
print(`For explanation of each procedure, type "WhatIs(procedure_name)"`):
print():
 print(`----------------------------------------------------`):


print(`----------------------------------------------------`):

with(combinat):

Help:=proc():
	print(`EstHappMenProp(n,K), EstHappWomenProp(n,K), AvgHapOverAllMatch(r)`):
	print(`AvgHapMenPropOverAllRank(n),AvgHapWomenPropOverAllRank(n)`):
	print(`AvgHapOverAllRankMatch(n),EstHappAllMatch(n,K)`):
	print(`GaleShCollegeAdmissions(A,C,Q), TerseGaleShCollegeAdmissions(A,C,Q)`):
end:



Help1:=proc():
	print(` RT(n), IsStable1(M,W,pi,i1,i2), IsStable(M,W,pi), FindStable(M,W,K) `): 
	print(`CountInstablities(M,W,pi), GFstab(n,K,x) `):
	print(` inv(pi), GaleSh(M,W), GaleShT(M,W), GaleShTw(M,W) `):
	print(`CheckOpt(M,W),CompGaleSh(M,W),EstStat(n,K),FindNeighbors(pi)`):
	print(`FindHappiness(M,W,matching), AllRank(n), AllMatchings(M,W), SlowAvgHapOverAllMatch(r)`):
	print(`inverse(pi,N)`):
end:

WhatIs:=proc(procedure):
	if procedure=FindHappiness then
		print(`FindHappiness(M,W,matching): inputs ranking tables M,W, and a matching`):
		print(`outputs the average happiness of the men and the average happiness of the women,`):
		print(`where a person's happiness is the ranking they gave their spouse they are matched with`):
	elif procedure=EstHappMenProp then
		print(`EstHappMenProp(n,K): inputs an integer n and a large integer K,`):
		print(`generates K examples of M,W (using RT(n)), applies the Gale-Shapley algorithm where the men propose,`):
		print(`and outputs the average happiness of the men and the average happiness of the women`):
	elif procedure=EstHappWomenProp then
		print(`EstHappWomenProp(n,K): inputs an integer n and a large integer K,`):
		print(`generates K examples of M,W (using RT(n)), applies the Gale-Shapley algorithm where the women propose,`):
		print(`and outputs the average happiness of the men and the average happiness of the women`):
	elif procedure=AllRank then
		print(`AllRank(n): The list of all ways n men (respectively, women) can rank n women (respectively, men)`):
	elif procedure=AllMatchings then 
		print(`AllMatchings(M,W): inputs ranking tables M,W and outputs all stable matchings`):
	elif procedure=SlowAvgHapOverAllMatch then
		print(`SlowAvgHapOverAllMatch(M,W):  inputs ranking tables M,W`):
		print(`and finds the average happiness of the men and the average happiness of the women over all stable matchings`):
	elif procedure=AvgHapOverAllMatch then
		print(`AvgHapOverAllMatch(M,W):  inputs ranking tables M,W`):
		print(`and finds the average happiness of the men and the average happiness of the women over all stable matchings`):
	elif procedure=AvgHapMenPropOverAllRank then
		print(`AvgHapMenPropOverAllRank(n): Inputs an integer n and calculates, over all possible rankings of n men and n women,`):
		print(`the average happiness of the men and the average happiness of the women when men are proposing in the Gale-Shapley algorithm`):
	elif procedure=AvgHapWomenPropOverAllRank then
		print(`AvgHapWomenPropOverAllRank(n): Inputs an integer n and calculates, over all possible rankings of n men and n women,`):
		print(`the average happiness of the men and the average happiness of the women when women are proposing in the Gale-Shapley algorithm`):
	elif procedure=AvgHapOverAllRankMatch then
		print(`AvgHapOverAllRankMatch(n): Inputs an integer n and calculates the average happiness of men and the average happiness of women`):
		print(`over all possible rankings of n men and women and all stable matchings of each ranking`):
	elif procedure=EstHappAllMatch then
		print(`EstHappAllMatch(n,K): Inputs an integer n and a large integer K,`):
		print(`generates K examples of M,W (using RT(n)),`):
		print(`and outputs the average happiness of the men and the average happiness of the women`):
		print(`over all stable matchings of the K rankings`):
	elif procedure=RT then
		print(`RT(n): Generates a pair of random preference M,W, two lists-of-lists`):
	elif procedure=GaleShCollegeAdmissions then
		print(`GaleShCollegeAdmissions(A,C,Q): The Gale-Shapley algorithm: admissions version`):
	elif procedure=TerseGaleShCollegeAdmissions then
		print(`TerseGaleShCollegeAdmissions(A,C,Q): The Gale-Shapley algorithm: admissions version (TERSE)`):
	elif procedure=inverse then
		print(`inverse(pi,N): The inverse of an incomplete permutation pi`):
	fi:
	
end:


#RT(n): Generates a pair of random preference 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 preference 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:



#CountInstablities(M,W,pi): inputs ranking tables M,W, and a current matching pi,
# and counts how many of the n*(n-1) ordered pairs of husbands 1<=i1,i2<=n are such that
# i1 is married to j1, i2 is married to j2,
# but i1 would rather be married to j2 (instead of his current wife j1)
# and j2 would rather be married to i1 (rather than her current husband i2)

CountInstabilities:=proc(M,W,pi) local i1, i2, leaving:
 leaving:=0:
 for i1 from 1 to nops(pi) do
  for i2 from 1 to nops(pi) do 
   if not(IsStable1(M,W,pi,i1,i2)) then leaving+=1: fi:
  od:
 od:
 leaving:
end:
#Note that IsItStable1 returns true if i1=i2 so we do not need to account for this case.

#GFstab(n,K,x): inputs positive integers n and a large positive integer K,
# generates K times, random ranking tables using RT(n), and a random matching, pi,
# and outputs the polynomial in x (of degree d n(n-1)),
# whose coeff. of xi is the number of those random choices that lead to exactly i instablities.

GFstab:=proc(n,K,x) local R,pi,k,poly:
 poly:=0:
 for k from 1 to K do
  R:=RT(n):
  pi:=randperm(n):
  poly+=x^(CountInstabilities(R[1],R[2],pi)):
 od:
 poly:
end:

#By taking the first and second derivatives w.r.t. to x,
# estimate the expectation and variance of the "number of instabilities".

ExpVarNumInstab:=proc(n,K) local p,dp:
	p:=GFstab(n,K,x):
 	dp:=diff(p,x):
 	evalf(subs(x=1,dp)/K),evalf(subs(x=1,diff(dp,x)/K+dp/K-(dp/K)^2)) :
end:



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


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


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

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


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:



#EstStat(n,K): inputs a positive integer n , and a large positive integer K (at least 1000),
# generates K examples 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:


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







#FindHappiness(M,W,matching): inputs ranking tables M,W, and a matching,
#outputs the average happiness of the men and the average happiness of the women,
#where a person's happiness is the ranking they gave their spouse they are matched with
FindHappiness:=proc(M,W,matching) local i,j, men, women,n:
	n:=nops(matching):
	men:=0:women:=0:
	for i from 1 to n do
		men+=M[i][matching[i]]:
		women+=W[matching[i]][i]:
	od:
[men/n, women/n]:
end:




#EstHappMenProp(n,K): inputs an integer n and a large integer K,
#generates K examples of M,W (using RT(n)), applies the Gale-Shapley algorithm where the men propose,
#and outputs the average happiness of the men and the average happiness of the women

EstHappMenProp:=proc(n,K) local i,R,M,W, avgM, avgW, FH:
	avgM:=0:avgW:=0:
	for i from 1 to K do
		R:=RT(n):
		M:=R[1]: W:=R[2]:
		FH:=FindHappiness(M,W, GaleShT(M,W)[1]):
		avgM+=FH[1]: avgW+=FH[2]:
	od:
	[avgM/K, avgW/K]:
end:


#EstHappWomenProp(n,K): inputs an integer n and a large integer K,
#generates K examples of M,W (using RT(n)), applies the Gale-Shapley algorithm where the women propose,
#and outputs the average happiness of the men and the average happiness of the women

EstHappWomenProp:=proc(n,K) local i,R,M,W, avgM, avgW, FH:
	avgM:=0:avgW:=0:
	for i from 1 to K do
		R:=RT(n):
		M:=R[1]: W:=R[2]:
		FH:=FindHappiness(M,W, GaleShTw(M,W)[1]):
		avgM+=FH[1]: avgW+=FH[2]:
	od:
	[avgM/K, avgW/K]:
end:


#AllMatchings(M,W): inputs ranking tables M,W and outputs all stable matchings
AllMatchings:=proc(M,W)local n,Match,P,m:
	Match:={}:
	n:=nops(M):
	P:=permute(n):
	for m in P do
		if IsStable(M,W,m)[1] then Match:=Match union {m}: fi:
	od:
	Match:
end:

#SlowAvgHapOverAllMatch(M,W):  inputs ranking tables M,W
#and finds the average happiness of the men and the average happiness of the women
#over all stable matchings 
SlowAvgHapOverAllMatch:=proc(M,W) local Match, avgM, avgW, m,FH,K:
	avgM:=0:avgW:=0:
	Match:=AllMatchings(M,W):
	K:=nops(Match):
	for m in Match do
		FH:=FindHappiness(M,W,m):
		avgM+=FH[1]: avgW+=FH[2]:
	od:
	[avgM/K, avgW/K]:
end:

#AvgHapOverAllMatch(M,W):  inputs ranking tables M,W
#and finds the average happiness of the men and the average happiness of the women
#over all stable matchings 
AvgHapOverAllMatch:=proc(M,W) local n,K, avgM, avgW, FH,m, P:
	n:=nops(M):
	K:=0: avgM:=0:avgW:=0:
	P:=permute(n):
	for m in P do
		if IsStable(M,W,m)[1] then
			FH:=FindHappiness(M,W,m):
			avgM+=FH[1]: avgW+=FH[2]:
			K+=1:
		fi:
	od:
	[avgM/K, avgW/K]:
end:



#AllRank(n): The list of all ways n men (respectively, women) can rank n women (respectively, men)
AllRank:=proc(n) local p,L,i:
	p:=permute(n):
	L:=seq([r],r in p): 
	for i from 1 to n-1 do
 		L:=seq(seq([op(r),r1], r in L), r1 in p   ):
	od:
 	[L]:
end:


#AvgHapMenPropOverAllRank(n): Inputs an integer n and calculates, over all possible rankings of n men and n women,
#the average happiness of the men and the average happiness of the women 
#when men are proposing in the Gale-Shapley algorithm
AvgHapMenPropOverAllRank:=proc(n) local allRank, avgM, avgW, M,W, FH:
	avgM:=0:avgW:=0:
	allRank:=AllRank(n):

	for M in allRank do
		for W in allRank do
			FH:=FindHappiness(M,W, GaleShT(M,W)[1]):
			avgM+=FH[1]: avgW+=FH[2]:
		od:
	od:
	[avgM/(n!)^(2*n), avgW/(n!)^(2*n)]:
end:

#AvgHapWomenPropOverAllRank(n):  Inputs an integer n and calculates, over all possible rankings of n men and n women,
#the average happiness of the men and the average happiness of the women 
#when women are proposing in the Gale-Shapley algorithm

AvgHapWomenPropOverAllRank:=proc(n) local allRank, avgM, avgW, M,W, FH:
	avgM:=0:avgW:=0:
	allRank:=AllRank(n):
	for M in allRank do
		for W in allRank do
			FH:=FindHappiness(M,W, GaleShTw(M,W)[1]):
			avgM+=FH[1]: avgW+=FH[2]:
		od:
	od:
	[avgM/(n!)^(2*n), avgW/(n!)^(2*n)]:
end:


#AvgHapOverAllRankMatch(n): Inputs an integer n and calculates the average happiness of men and the average happiness of women
#over all possible rankings of n men and women and all stable matchings of each ranking 
AvgHapOverAllRankMatch:=proc(n) local allRank, avgM, avgW, M,W, AH:
	avgM:=0:avgW:=0:
	allRank:=AllRank(n):
	for M in allRank do
		for W in allRank do
			AH:=AvgHapOverAllMatch(M,W):
			avgM+=AH[1]: avgW+=AH[2]:
		od:
	od:
	[avgM/(n!)^(2*n), avgW/(n!)^(2*n)]:
end:


#EstHappAllMatch(n,K): Inputs an integer n and a large integer K,
#generates K examples of M,W (using RT(n)),
#and outputs the average happiness of the men and the average happiness of the women
#over all stable matchings of the K rankings
EstHappAllMatch:=proc(n,K) local i,R,M,W, avgM, avgW, AH:
	avgM:=0:avgW:=0:
	for i from 1 to K do
		R:=RT(n):
		AH:=AvgHapOverAllMatch(R):
		avgM+=AH[1]: avgW+=AH[2]:
	od:
	[avgM/K, avgW/K]:
end:
 


with(combinat):
#inverse(pi,N): given a list of integers length m<=N (with no repeated integers and integers are between 1 and N) 
# inverse(pi,N) adds the remaining N-m integers to the end of the list to turn it into a permutation on N letters;
# it then inverts the permutation. 

inverse:=proc(pi,N) local i,j,T,p1:
p1:=pi:
if not nops(pi)=N then

for i from 1 to N do
 if not member(i,pi) then
 p1:=[op(p1),i]: 
 fi:
od:
fi: 
for j from 1 to N do
 T[p1[j]]:=j:
od:
[seq(T[j],j=1..N)  ]:
end:


#GaleShCollegeAdmissions(A,C,Q): The Gale-Shapley algorithm: admissions version (VERBOSE Version)
#  Given $A$ (the set of $N$ student preference lists), $C$ (the set of $M$ college preference lists),
# and $Q$ (the list of quotas for the $M$ colleges), such that M=|C|=|Q|, and college i has quota Q[i],
# this algorithm will output  a list of students accepted at each college in the form f
#This will also print out what is happening during each round of applications 

GaleShCollegeAdmissions:=proc(A,C,Q) local M,N,i,j,k,co,c,CRank,waitlist,rejected,applicants,worst,Acopy,left,finished:
N:=nops(A):
M:=nops(C):
co:=0:

if not (type(C,list) and type(A,list)) or not nops(C)=nops(Q) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

Acopy:=A:
CRank:=[seq(inverse(C[i],N),i=1..M)]:

#nobody has been accepted yet.
rejected:={seq(i,i=1..N)}:

#initializing waitlists
waitlist:=[seq({},i=1..M)]:

#student i is applying to college _  (0 if not applying)
applicants:=[seq(0,i=1..N)]:

finished:=false:

while not finished do

 co:=co+1:
  
 #each applicant applies to its preferred college
 printf(`\n beginning of application round %d \n\n`,co):
 for i from 1 to N do
  printf(`Is student %d applying this round? %a \n`,i, member(i,rejected)):   
  if member(i,rejected) then

   while not evalb(Acopy[i]=[]) and not member(i,C[Acopy[i][1]]) do
    printf(`%d tried to apply to college %d but they are not allowed to \n`,i, Acopy[i][1] ):
    Acopy[i]:=subsop(1=NULL, Acopy[i]):
   od:
   
   if not evalb(Acopy[i]=[]) then
    applicants[i]:=Acopy[i][1]:
    printf(`preference list for student %d is %a \n`,i, Acopy[i]):
    Acopy[i]:=subsop(1=NULL,Acopy[i]):   
    printf(`Applicant %d is applying to college %d \n`,i, applicants[i]):
   else
    applicants[i]:=0:
   fi:

  else 
   applicants[i]:=0:
  fi:
 od:
  
 printf(`\n beginning of rejecting/waitlisting students \n\n`):
 
 for j from 1 to N do
   
  if applicants[j]<>0 then

   if Q[applicants[j]]>nops(waitlist[applicants[j]]) then
    printf(`Quota is non empty so %d is added to the waitlist and taken off of the rejected list. \n`,j):
    waitlist[applicants[j]]:=waitlist[applicants[j]] union {j}:
    rejected:=rejected minus {j}:   
   else
    printf(`Quota is  empty so we need to make a decision about %d \n`,j):
    #We need to make a decision.
    #Who is the highest (worst) ranked person in my waitlist?
    worst:=0:

    for k in waitlist[applicants[j]] do    
     if CRank[applicants[j]][k]>worst then
      worst:=k:
     fi:
    od:

    printf(`The worst applicant in the waitlist is %d \n`, worst):

    if CRank[applicants[j]][j]<CRank[applicants[j]][worst] then
     printf(`%d is better ranked than the worst, so we replace the worst in the waitlist with %d \n`,j,j):

     rejected:=rejected union {worst}: 
     waitlist[applicants[j]]:=waitlist[applicants[j]] minus {worst}:
     waitlist[applicants[j]]:= waitlist[applicants[j]] union {j}:
     rejected:= rejected minus {j}:

    else
     printf(`%d is not better ranked than the worst, so we don't change anything, they are still rejected. \n`,j):
    fi:

   fi: 
    
  fi:

 od:

 #are we done yet?
 if evalb(rejected={}) then
  printf(`\n All students are on a waitlist! \n`):
  finished:=true:
 elif add(i, i in [seq(nops(Acopy[j]),j in rejected)])=0 then
  printf(`\n The rejected students have nowhere else to apply. \n`):
  finished:=true:
 fi:

od:

printf(`\nThe stable matching gotten from Gale Shapley is as follows: \n`):

for c from 1 to M do
 left:=Q[c]-nops(waitlist[c]):
 if left=0 then
  printf(`College %d accepted applicants %s and filled it's quota \n`,c,convert(waitlist[c],string)):
 else
  printf(`College %d accepted applicants %s and wished to accept %d more students \n`,c,convert(waitlist[c],string),left):
 fi:

od:

printf(`\nThe following students were not accepted to any colleges:   `):
printf(` %s \n`, convert(rejected, string) ):
printf(`The algorithm took %d rounds of applications \n`,co):

waitlist:

end:



#TerseGaleShCollegeAdmissions(A,C,Q): The Gale-Shapley algorithm: admissions version (TERSE)
#  Given $A$ (the set of $N$ student preference lists), $C$ (the set of $M$ college preference lists),
# and $Q$ (the list of quotas for the $M$ colleges), such that M=|C|=|Q|, and college i has quota Q[i],
# this algorithm will output  a list of students accepted at each college in the form f

TerseGaleShCollegeAdmissions:=proc(A,C,Q) local M,N,i,j,k,co,c,CRank,waitlist,rejected,applicants,worst,Acopy,left,finished:
N:=nops(A):
M:=nops(C):
co:=0:

if not (type(C,list) and type(A,list)) or not nops(C)=nops(Q)  then
 print(`Bad input`):
 RETURN(FAIL):
fi:

Acopy:=A:
CRank:=[seq(inverse(C[i],N),i=1..M)]:

#nobody has been accepted yet.
rejected:={seq(i,i=1..N)}:

#initializing waitlists
waitlist:=[seq({},i=1..M)]:

#student i is applying to college _  (0 if not applying)
applicants:=[seq(0,i=1..N)]:

finished:=false:

while not finished do

 co:=co+1:
  
 #each applicant applies to its preferred college
 for i from 1 to N do

  if member(i,rejected) then
  
   while not evalb(Acopy[i]=[]) and not member(i,C[Acopy[i][1]]) do
    Acopy[i]:=subsop(1=NULL, Acopy[i]):
   od:


   if not evalb(Acopy[i]=[]) then
    applicants[i]:=Acopy[i][1]: 
    Acopy[i]:=subsop(1=NULL,Acopy[i]):
   else
    applicants[i]:=0:
   fi:
  else 
   applicants[i]:=0:    
  fi:
 od:
  
 for j from 1 to N do  
  if applicants[j]<>0 then
   if Q[applicants[j]]>nops(waitlist[applicants[j]]) then
    waitlist[applicants[j]]:=waitlist[applicants[j]] union {j}:
    rejected:=rejected minus {j}:
   else
    #We need to make a decision.
    #Who is the highest (worst) ranked person in my waitlist?
    worst:=0:
    for k in waitlist[applicants[j]] do
     if CRank[applicants[j]][k]>worst then
      worst:=k:
     fi:
    od:

    if CRank[applicants[j]][j]<CRank[applicants[j]][worst] then

     rejected:=rejected union {worst}: 
     waitlist[applicants[j]]:=waitlist[applicants[j]] minus {worst}:
     waitlist[applicants[j]]:= waitlist[applicants[j]] union {j}:
     rejected:= rejected minus {j}:
    fi:
   fi: 
  fi:
 od:

 #are we done yet?
 if evalb(rejected={}) then
  finished:=true:
 elif add(i, i in [seq(nops(Acopy[j]),j in rejected)])=0 then
  finished:=true:
 fi:

od:

waitlist:

end: