#OK to post homework
#Blair Seidler, 2/27/22, Assignment 11
with(combinat):
with(Statistics):
Help:=proc(): print(`GaleShTw(M,W), CheckOpt(M,W), GSPartial(M,W)`): end:

#1. Not entirely sure what was being asked for here.

#2. GaleShTw(M,W), where Women do the proposing and men respond. Hint: it should be a one-liner.
# Make sure to invert the first output of GaleShTw(W,M).
GaleShTw:=proc(M,W) local m:
m:=GaleShT(W,M): inv(m[1]),m[2]:
end:

#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 m,w,mi,wi,i,j,n:
 n:=nops(M):
 m:=GaleShT(M,W)[1]:
 mi:=inv(m):
 w:=GaleShTw(M,W)[1]:
 wi:=inv(w):
 for i from 1 to n do
  if M[i][m[i]]>M[i][w[i]] then
    printf("Man %d got choice %d when men proposed, but choice %d when women proposed!",
           i,M[i][m[i]],M[i][w[i]]):
    return(false):
  fi:
  if W[i][wi[i]]>W[i][mi[i]] then
    printf("Woman %d got choice %d when men proposed, but choice %d when women proposed!",
           i,W[i][mi[i]],W[i][wi[i]]):
    return(false):
  fi:
 end:
 return(true):
end:

(* Output:

This seems to have worked based on the following (one of several runs I tried):

> {seq(CheckOpt(RT(30)), i = 1 .. 1000)};
                             {true}
*)

#4. GSPartial(M,W): note that the inputs here cannot be rankings, because not every person ranks
# every person from the other set. These have to be akin to M[1]=[4,2], where Mr. 1 prefers Ms. 4,
# would marry Ms.2, and would rather remain single than marry Ms. 1, Ms. 3, or Ms. 5+ (if they exist).
GSPartial:=proc(M,W) local full,i,j,n,m,w,tm,tw,match:
n:=nops(M):
full:={seq(i,i=1..n)}:
m:=[seq([],i=1..n),[op(full),n+1]]:
w:=[seq([],i=1..n),[op(full),n+1]]:
for i from 1 to n do
  m[i]:=inv([op(M[i]),n+1,op(full minus convert(M[i],set))]):
  w[i]:=inv([op(W[i]),n+1,op(full minus convert(W[i],set))]):
od:

match:=GaleShT(m,w)[1]:
if match[n+1]=n+1 then
  return([op(1..n,match)]):
else
  print("No complete stable matching exists"):
  return(FAIL):
fi:
end:

(* Output:

> M:=[[4,1],[1,3,4],[3],[4,1,3,2]]:
> W:=[[1,3],[4,1,2],[3],[2,4,3,1]]:
> GSPartial(M, W);
                          [1, 4, 3, 2]

*)

#5. Compare the "stupid" way, FindStable(M,W,K), with K=10000, with GaleShT(M,W) for many random
# choices given by RT(n). For small n, does the stupid way fare better sometimes?
# What is the "cutoff" of n where FindStable(M,W,K) starts to be useless?
TestFSGS:=proc(T,K) local n,R,i,smart,stupid,fail,FSout:
for n from 4 to 19 do
  fail:=0:
  stupid:=0:
  smart:=0:
  for i from 1 to T do
    R := RT(n):
    FSout:=[FindStable(R, K)]:
    if nops(FSout)>1 then
      stupid:=stupid+FSout[2]:
    else
      fail++:
    fi:
    smart:=smart+GaleShT(R)[2]:
  od:
  printf("n=%d, smart=%.1f, stupid=%.1f, stupid failures=%d\n",
          n,evalf(smart/T),evalf(stupid/(max(1,T-fail))),fail);
od:
for n from 20 to 100 by 10 do
  fail:=0:
  stupid:=0:
  smart:=0:
  for i from 1 to T do
    R := RT(n):
    FSout:=[FindStable(R, K)]:
    if nops(FSout)>1 then
      stupid:=stupid+FSout[2]:
    else
      fail++:
    fi:
    smart:=smart+GaleShT(R)[2]:
  od:
  printf("n=%d, smart=%.1f, stupid=%.1f, stupid failures=%d\n",
          n,evalf(smart/T),evalf(stupid/(max(1,T-fail))),fail);
od:
end:

(* Output:

> TestFSGS(10, 10000);
n=4, smart=6.0, stupid=3.3, stupid failures=0
n=5, smart=9.3, stupid=3.7, stupid failures=0
n=6, smart=12.5, stupid=7.0, stupid failures=0
n=7, smart=14.9, stupid=10.5, stupid failures=0
n=8, smart=19.6, stupid=10.7, stupid failures=0
n=9, smart=20.2, stupid=15.2, stupid failures=1
n=10, smart=24.5, stupid=22.8, stupid failures=0
n=11, smart=24.9, stupid=25.9, stupid failures=0
n=12, smart=27.0, stupid=29.9, stupid failures=1
n=13, smart=39.2, stupid=32.0, stupid failures=2
n=14, smart=37.5, stupid=37.6, stupid failures=3
n=15, smart=43.9, stupid=46.9, stupid failures=2
n=16, smart=51.1, stupid=59.6, stupid failures=3
n=17, smart=46.8, stupid=73.5, stupid failures=4
n=18, smart=62.4, stupid=88.2, stupid failures=6
n=19, smart=57.3, stupid=76.5, stupid failures=6
n=20, smart=55.0, stupid=148.8, stupid failures=6
n=30, smart=88.4, stupid=156.0, stupid failures=9
n=40, smart=167.4, stupid=0.0, stupid failures=10
n=50, smart=204.1, stupid=0.0, stupid failures=10
n=60, smart=279.8, stupid=0.0, stupid failures=10
n=70, smart=333.6, stupid=0.0, stupid failures=10
n=80, smart=341.9, stupid=0.0, stupid failures=10
n=90, smart=478.4, stupid=0.0, stupid failures=10
n=100, smart=499.6, stupid=0.0, stupid failures=10

The smart and stupid numbers are the average number of steps required to find a solution.
For small n, the stupid algorithm is faster. When n is in the low teens, the two algorithms
are about the same, but the stupid one fails occasionally. By the upper teens, GaleSh is
faster and FindStable is failing regularly. When n is 30 or higher, FindStable fails nearly
every time at 10000 iterations.
*)

#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 A,i:
A:=[seq(GaleShT(RT(n))[2],i=1..K)]:
printf("n=%d, Mean=%.2f, StdDev=%.2f\n",
        n,Mean(A),StandardDeviation(A)):
end:

(* Output:

Run EstStat(n,10000) for n=10..., 150, (with increments of 10)

> for n from 10 by 10 to 150 do
>     EstStat(n, 1000);
> end do;
n=10, Mean=23.83, StdDev=6.25
n=20, Mean=62.76, StdDev=15.61
n=30, Mean=106.37, StdDev=25.92
n=40, Mean=155.66, StdDev=40.75
n=50, Mean=208.03, StdDev=50.98
n=60, Mean=265.00, StdDev=60.52
n=70, Mean=319.91, StdDev=74.57
n=80, Mean=374.44, StdDev=86.80
n=90, Mean=436.56, StdDev=97.78
n=100, Mean=492.98, StdDev=113.48
n=110, Mean=552.05, StdDev=117.11
n=120, Mean=621.38, StdDev=138.27
n=130, Mean=682.71, StdDev=151.59
n=140, Mean=753.49, StdDev=167.67
n=150, Mean=804.92, StdDev=171.06

Can you conjecture an approximate "trend"?

1.08*n*ln(n) is a pretty good fit. I'm not sure why that should be the case. Is "wife collector"
related to coupon collector somehow?

*)

#7. Let the "neighbors" of a permutation pi of length n be the n-1 permutations obtained by swapping
# ADJACENT entries. For example, the three neighbors of [3,1,2,4] are [1,3,2,4], [3,2,1,4], [3,1,4,2]
# By using CountInstablities(M,W,pi) that you wrote for hw10, write a procedure
# MayBeNotSoStupidFindStable(M,W) that starts out with a random permutation, looks at all its n-1
# neighbors, and moves to the one with the least number of instablities. It keeps doing it until you
# either get a stable matching, or get into a "valley", where there are still instabilities, but none
# of the neighbors has less instabilities (then return FAIL). It should also record the number of
# iterations. How often does it give you a stable matching? Is it sometimes better than Gale-Shapley?

MayBeNotSoStupidFindStable:=proc(M,W) local n,pi,co,curi,neighbors,i,j,pi1,pi1i,bestpi:
K:=10000: #Just to avoid an infinite loop while debugging
n:=nops(M):
pi:=randperm(n):

for co from 1 to K do
 curi:=CountInstablities(M,W,pi):
 if curi=0 then
   RETURN(pi,co):
 fi:
 neighbors:={seq(seq([op(1..i-1,pi),pi[j],op(i+1..j-1,pi),pi[i],op(j+1..n,pi)],j=i+1..n),i=1..n-1)}:
 bestpi:=[]:
 for pi1 in neighbors do
   pi1i:=CountInstablities(M,W,pi1):
   if pi1i<curi then
     bestpi:=pi1:
     curi:=pi1i:
   fi:
  od:
 if bestpi=[] then return(FAIL): fi:
 pi:=bestpi:   
od:
FAIL:
end:

(* Output:

So this one is still pretty stupid. The assigned problem stated that we should only consider the n-1
permutations where we are swapping adjacent positions in the matching. That fails almost every time
even for small n. I expanded the neighbors to include all possible swaps, and it still fails almost
every time for n=15.

*)



#### From homework 10 ####
#CountInstablities(M,W,pi): enters 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)
CountInstablities:=proc(M,W,pi) local n,i1,i2,inst:
n:=nops(pi):
inst:=0:
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
      inst:=inst+1:
   fi:
 fi:
od:
od:
inst:

end:

#### From C11.txt ####

#Algorithm:
#Each person starts out with no people "cancelled" 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: