#OK to post homework 
#Natalya Ter-Saakov, Feb 27, Assignment 10 
read "C10.txt":

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


for n from 1 to 100 do n,ExpVarNumInstab(n, 10000); od;
#                           1, 0., 0.
#                 2, 0.4913000000, 0.3769243100
#                  3, 1.495000000, 1.370375000
#                  4, 2.977000000, 3.219471000
#                  5, 5.000900000, 6.112099190
#                  6, 7.462800000, 10.54341616
#                  7, 10.59080000, 16.59875536
#                  8, 14.06830000, 24.00983511
#                  9, 18.00420000, 34.93798236
#                  10, 22.57610000, 46.89540879
#                  11, 27.49270000, 62.91014671
#                  12, 33.07160000, 79.12227344
#                  13, 38.97790000, 103.5942116
#                  14, 45.59930000, 127.1187395
#                  15, 52.52900000, 151.2877590
#                  16, 59.83020000, 180.5517680
#                  17, 68.00530000, 217.1290719
#                  18, 76.36310000, 266.4026584
#                  19, 85.45200000, 314.5442960
#                  20, 94.55920000, 357.6048954
#                  21, 105.0486000, 412.1516380
#                  22, 115.5683000, 472.7559351
#                  23, 126.4387000, 526.5080423
#                  24, 137.9064000, 605.3532390
#                  25, 150.2319000, 700.0603224
#                  26, 162.7155000, 748.4375598
#                  27, 175.4229000, 869.1018556
#                  28, 189.1055000, 970.4501698
#                  29, 203.1577000, 1085.534431
#                  30, 217.6044000, 1175.402701
#                  31, 232.5406000, 1294.507952
#                  32, 247.9026000, 1418.830113
#                  33, 263.6921000, 1517.888698
#                  34, 280.9677000, 1721.880057
#                  35, 298.4970000, 1857.673191
#                  36, 314.8797000, 2024.217228
#                  37, 333.1149000, 2185.022898
#                  38, 352.1577000, 2364.929631
#                  39, 370.7772000, 2513.108160
#                  40, 390.3736000, 2777.585423
#                  41, 409.8001000, 2912.591940
#                  42, 430.3268000, 3141.185602
#                  43, 450.2561000, 3454.975313
#                  44, 472.7923000, 3562.627961
#                  45, 495.5295000, 4005.073530
#                  46, 517.3841000, 4124.613567
#                  47, 539.7064000, 4429.450399
#                  48, 562.7271000, 4792.343626
#                  49, 588.1834000, 5020.392564
#                  50, 612.8343000, 5347.940044
#                  51, 637.7541000, 5669.550233
#                  52, 662.3917000, 6146.434671
#                  53, 688.7645000, 6323.493040
#                  54, 716.5615000, 6774.731018
#                  55, 742.5893000, 7147.935226
#                  56, 769.4146000, 7499.946107
#                  57, 797.8646000, 7923.934067
#                  58, 826.9617000, 8561.904633
#                  59, 857.7900000, 8638.830900
#                  60, 884.3801000, 9144.622624
#                  61, 915.5807000, 9550.630288
#                  62, 944.7394000, 10274.68129
#                  63, 976.0413000, 10560.03679
#                  64, 1006.151600, 11350.94342
#                  65, 1037.736400, 11772.64112
#                  66, 1072.500400, 12227.62040
#                  67, 1106.285800, 12804.25692
#                  68, 1139.548700, 13752.43363
#                  69, 1172.162600, 13993.79416
#                  70, 1206.698400, 14438.83024
#                  71, 1242.440700, 14886.82388
#                  72, 1275.644100, 16130.50524
#                  73, 1314.205000, 16438.73738
#                  74, 1350.003400, 17427.87819
#                  75, 1385.163000, 17718.17723
#                  76, 1424.226800, 18642.02316
#                  77, 1465.645300, 19203.26489
#                  78, 1502.766500, 20474.13938
#                  79, 1540.422600, 20894.87521
#                  80, 1579.444900, 21802.63916
#                  81, 1620.721800, 22885.10680
#                  82, 1657.929900, 22831.48239
#                  83, 1701.154400, 24890.56196
#                  84, 1744.185800, 25175.20128
#                  85, 1783.283000, 25700.94751
#                  86, 1830.310800, 26752.30940
#                  87, 1870.213400, 27926.71406
#                  88, 1916.337100, 28970.18366
#                  89, 1956.230700, 28973.56948
#                  90, 2003.904300, 30851.87914
#                  91, 2049.320100, 32133.18544
#                  92, 2094.633400, 33049.95360
#                  93, 2141.501600, 34408.47320
#                  94, 2185.149800, 35216.51496
#                  95, 2232.426200, 36307.50195
#                  96, 2275.155700, 37578.98826
#                  97, 2327.824200, 38671.50629
#                  98, 2377.315400, 40792.88172
#                  99, 2428.059700, 39408.71514
#                 100, 2476.401500, 43050.90870