#OK to post homework
#Victoria Chayes, 2/28/21, Assignment 10

Help:=proc(): print(`ExactGFinv(n,x), ExactGFmaj(n,x), ExactGFdes(n,x), EstGFinv(n,x,K) ,EstGFmaj(n,x,K) , EstGFdes(n,x,K), PnSeq(N,t)`): end:
with(combinat):


#1: Given a permutation of {1,...,n} the number of inversion, denoted by inv(pi), is the number of pairs 1 <= i < j <= n such pi[i] > pi[j]. [For example inv(631425)=8]
# The number of descents, denoted by des(pi), is the number of places i, such that pi[i] > pi[i+1] [For example inv(631425)=3]
# The major index, denoted by maj(pi), is the SUM of the places i, such that pi[i] > pi[i+1] [For example maj(631425)=1+2+4=7].
# Write procedures ExactGFinv(n,x) ,ExactGFmaj(n,x) , ExactGFdes(n,x) that inputs a positive integer n and a variable name n, and outputs the generating function whose coefficient of x^k is the exact number of permutations with inv=k, maj=k, des=k, respectively.

ExactGFinv:=proc(n,x) local n1, K,S,i,j,k,count,s:
n1:=add(i,i=1..(n-1)):
K:=[0$(n1+1)]:
S:=permute(n):
for i from 1 to nops(S) do
count:=1:
s:=S[i]:
 for j from 1 to (n-1) do
  for k from (j+1) to n do
    if s[j]>s[k] then
      count:=count+1
    fi:
  od:
 od:
K[count]:=K[count]+1:
od:
add(K[i+1]*x^i, i=0..n1):
end:

ExactGFmaj:=proc(n,x) local K,S,s,i,j,count:
K:=[0$n]:
S:=permute(n):
for i from 1 to nops(S) do
count:=1:
s:=S[i]:
  for j from 1 to (n-1) do
    if s[j]>s[j+1] then
      count:=count+1
    fi:
  od:
K[count]:=K[count]+1
od:
add(K[i+1]*x^i, i=0..(n-1)):
end:

ExactGFdes:=proc(n,x) local n1, K,S,s,i,j,count:
n1:=add(i, i=1..(n-1)):
K:=[0$(n1+1)]:
S:=permute(n):
for i from 1 to nops(S) do
count:=1:
s:=S[i]:
  for j from 1 to (n-1) do
    if s[j]>s[j+1] then
      count:=count+j
    fi:
  od:
K[count]:=K[count]+1
od:
add(K[i+1]*x^i, i=0..n1):
end:

# Try to conjecture (or derive mathematically) closed form or efficient recurrences for these generating functions that will enable you to easily compute them for say n=100. [Warning: I don't think that the g.f. according to des has a nice closed form, but it has a recurrence]

# First, we note that ExactGFinv and ExactGFdes produce the same polynomials, so we only have search for the pattern with one of them. Looking at the progression, we see:
#ExactGFinv(2,x)
#                             1 + x
#
#ExactGFinv(3,x)
#                       3      2          
#                      x  + 2 x  + 2 x + 1
#
#ExactGFinv(4,x)
#             6      5      4      3      2          
#            x  + 3 x  + 5 x  + 6 x  + 5 x  + 3 x + 1
#
#ExactGFinv(5,x)
# 10      9      8       7       6       5       4       3      2
#x   + 4 x  + 9 x  + 15 x  + 20 x  + 22 x  + 20 x  + 15 x  + 9 x  + 4 x + 1
#
#ExactGFinv(6,x)
# 15      14       13       12       11       10       9        8
#x   + 5 x   + 14 x   + 29 x   + 49 x   + 71 x   + 90 x  + 101 x 
#
#          7       6       5       4       3       2          
#   + 101 x  + 90 x  + 71 x  + 49 x  + 29 x  + 14 x  + 5 x + 1
# or coefficient sequences:
# 1, 1
# 1, 2, 2, 1
# 1, 3, 5, 6, 5, 3, 1
# 1, 4, 9, 15, 20, 22, 20, 15, 9, 4, 1
# 1, 5, 14, 29, 49, 71, 90, 101, 101, 90, 71, 49, 29, 14, 5, 1
#
# There is some sort of triangle-like setup going on, with the following properties:
# --the number of terms is sequence A000124, n(n+1)/2 + 1
# --the first coefficient is always 1
# --the second coefficient seems to be (n-1)
# and from there it is harder for me to see an obvious pattern
#
# For ExactGFmaj, 
# ExactGFmaj(2,x)
#                             1 + x
#
# ExactGFmaj(3,x)
#                           2          
#                          x  + 4 x + 1
#
# ExactGFmaj(4,x)
#                      3       2           
#                     x  + 11 x  + 11 x + 1
#
# ExactGFmaj(5,x)
#                  4       3       2           
#                 x  + 26 x  + 66 x  + 26 x + 1
#
# ExactGFmaj(6,x)
#             5       4        3        2           
#            x  + 57 x  + 302 x  + 302 x  + 57 x + 1
#
# ExactGFmaj(7,x)
#      6        5         4         3         2            
#     x  + 120 x  + 1191 x  + 2416 x  + 1191 x  + 120 x + 1
#
# The patterns that immediately emerge are:
# --n terms 
# --same triangle-like formation 
# --the first coefficient seems to be sequence A000295, the Eulerian numbers 
# --the second coefficient seems to be sequence A000460, the Eulerian numbers second column, which redirects to sequence A008292 that I believe will give all the coefficients. 



#2 By using randperm(n), write procedures EstGFinv(n,x,K) ,EstGFmaj(n,x,K) , EstGFdes(n,x,K) that samples K random permutations and finds the empircal prob. gen. function. Compare plotDist and SAc for up to the 8th moment, with n=100, and K=2000.

EstGFinv:=proc(n,x,K) local n1, S,i,j,k,count,s:
n1:=add(i,i=1..(n-1)):
S:=Array([0$(n1+1)]):
for i from 1 to K do
count:=1:
s:=randperm(n):
 for j from 1 to (n-1) do
  for k from (j+1) to n do
    if s[j]>s[k] then
      count:=count+1
    fi:
  od:
 od:
S[count]:=S[count]+1:
od:
(n!/K)*add(S[i+1]*x^i, i=0..n1):
end:

EstGFmaj:=proc(n,x,K) local S,s,i,j,count:
S:=[0$n]:
for i from 1 to K do
count:=1:
s:=randperm(n):
  for j from 1 to (n-1) do
    if s[j]>s[j+1] then
      count:=count+1
    fi:
  od:
S[count]:=S[count]+1
od:
(n!/K)*add(S[i+1]*x^i, i=0..(n-1)):
end:

EstGFdes:=proc(n,x,K) local n1, S,s,i,j,count:
n1:=add(i, i=1..(n-1)):
S:=Array([0$(n1+1)]):
for i from 1 to K do
count:=1:
s:=randperm(n):
  for j from 1 to (n-1) do
    if s[j]>s[j+1] then
      count:=count+j
    fi:
  od:
S[count]:=S[count]+1
od:
(n!/K)*add(S[i+1]*x^i, i=0..n1):
end:

# Pictures of PlotDist are given as hw10VictoriaChayesfig1, 2, 3. 
# evalf(SAc(EstGFinv(100,x,2000),x,8))
#    [2472.001000, 27667.03500, -0.07307826647, 2.914031516, -1.047584078, 14.41900431, -13.98421632, 108.5064884]
# evalf(SAc(EstGFmaj(100,x,2000),x,8))
#    [49.58150000, 8.258357750, -0.04802876048, 2.800543896, -0.4695491758, 12.72958960, -5.462371712, 80.74816403]
# evalf(SAc(EstGFdes(100,x,2000),x,8))
#    [2472.989000, 29073.07188, -0.06945474171, 2.907227463, -0.7919565967, 13.33799650, -7.994123298, 80.06813228]



###########################
### Programs From Class ###
###########################

SAc:=proc(f,x,K) local i,mu,sig2,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:
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):
M:=[mu,sig2]:
if sig2=0 then RETURN(M): fi:
for i from 3 to K do
  f1:=x*diff(f1,x):
  M:=[op(M), subs(x=1,f1)/sig2^(i/2)]:
od:
M:
end:

PlotDist:=proc(f,x,K) local i,mu,sig,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:
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):
sig:=sqrt(subs(x=1,f1)):
if sig=0 then RETURN(FAIL): fi:
plot([seq([(i-mu)/sig,coeff(f,x,i)/subs(x=1,f)*sig],i=trunc(mu-K*sig)..trunc(mu+K*sig))]):
end:

RandSPn:=proc(n) local i,k,die1:
die1:=[seq(NuSPnk(n,i),i=1..n)]:
k:=LD(die1):
RandSPnk(n,k):
end:

LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:
for i from 1 to n do
su:=su+L[i]:
if r<=su then
RETURN(i):
fi:
od:
end:

NuSPnk:=proc(n,k) option remember:
if n=1 then
if k=1 then
RETURN(1):
else
RETURN(0):
fi:
fi:
NuSPnk(n-1,k-1)+ k*NuSPnk(n-1,k):
end: