# OK to post homework
# Lucy Martinez, 02-27-22, Assignment 10
read `C10.txt`:
#----------------------Problem 1-----------------------#
# Write a procedure CountInstabilities(M,W,pi)
# That enters ranking tables M,W, and a current matching pi, and counts how many of the n*(n-1) ordered pairs of
# husbands 1 d i1, i2 d 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 count,i1,i2:
count:=0:
for i1 from 1 to nops(pi) do
for i2 from 1 to nops(pi) do
 if i1<>i2 then
  if not IsStable1(M,W,pi,i1,i2) then
   count:=count+1;
  fi:
 fi:
od:
od:
count:
end:
#----------------------Problem 2-----------------------#
# Write a procedure GFstab(n,K,x)
# That 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. By taking the first and second 
# derivatives w.r.t. to x, estimate the expectation and variance of the "number of instabilities".

GFstab:=proc(n,K,x) local i, T,P,pi,E,V,Der:
P:=0:
E:=0:
V:=0:
Der:=0:
for i from 1 to K do
 T:=RT(n):
 pi:=randperm(n):
 P:=P+x^(CountInstabilities(T[1],T[2],pi)):
od:
Der:=diff(P,x):
E:=eval(Der/K,x=1):
V:=eval(diff(Der,x)/K,x=1)+E-E^2:
print(`Polynomial:`),print(P):
print(`Expectation of the Number of Instabilities:`),print(E):
print(`Variance`),print(V):
end:
# The following is just a small test:
GFstab(3,1000,x);
                          Polynomial:
          x^6 + 7*x^5 + 61*x^4 + 123*x^3 + 263*x^2 + 350*x + 195

             Expectation of the Number of Instabilities:

                           153/100

                           Variance:

                          13551/10000