#OK to post homework
#Blair Seidler, 2/23/22, Assignment 10
with(combinat):
Help:=proc(): print(`CountInstablities(M,W,pi),GFstab(n,K,x)`): end:

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

#2.
#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 ≤ 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 T,i,f:
add(x^CountInstablities(RT(n),randperm(n)),i in 1..K):
end:

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

MeanVar:=proc(f,x) local f1,mu,sig2:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig2:=subs(x=1,f1):
evalf(mu),evalf(sig2):
end:

(* Sample outputs:
F := GFstab(5, 10000, x);
          15       14       13       12        11        10
  F := 4 x   + 14 x   + 38 x   + 55 x   + 121 x   + 232 x  

            9        8         7         6         5         4
     + 408 x  + 677 x  + 1035 x  + 1346 x  + 1575 x  + 1510 x 

             3         2              
     + 1334 x  + 1010 x  + 501 x + 140

MeanVar(F, x);
                    4.990800000, 6.237515360



F := GFstab(10, 10000, x);
      53      52      51      49      48      47      46       45
F := x   + 2 x   + 3 x   + 3 x   + 2 x   + 4 x   + 9 x   + 10 x  

         44       43       42       41       40       39       38
   + 10 x   + 16 x   + 15 x   + 24 x   + 36 x   + 43 x   + 49 x  

         37       36       35        34        33        32
   + 60 x   + 83 x   + 98 x   + 147 x   + 178 x   + 227 x  

          31        30        29        28        27        26
   + 238 x   + 266 x   + 342 x   + 390 x   + 385 x   + 481 x  

          25        24        23        22        21        20
   + 458 x   + 520 x   + 606 x   + 604 x   + 589 x   + 566 x  

          19        18        17        16        15        14
   + 561 x   + 529 x   + 511 x   + 432 x   + 374 x   + 272 x  

          13        12        11        10       9       8
   + 242 x   + 199 x   + 152 x   + 108 x   + 66 x  + 40 x 

         7       6      5      4      3
   + 18 x  + 21 x  + 6 x  + 2 x  + 2 x 

MeanVar(F, x);
                    22.47360000, 47.11110304

*)

#### From C10.txt ####

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