print(`qVATTER Maple package by Andrew Baxter`):
print(`This package requires the VATTER package`):
print(`Last updated: 7/11/09`):
print(`Use the Help() command for more details`):

Help:=proc()
print(`Main procedures are:`):
print(`qWilf(n,Patterns,q): Brute-force computation of inv over S_n(Patterns)`):
print(`qMiklos(Scheme, prefix, gap vector, q): Uses Scheme to compute inv over permutations starting with prefix with given gap vector.  To all B-avoiding permutations of length n, use prefix=[] and gapvector=[n] `):
print(`qMiklos(SchemeFast(Patterns,?,?), [], [n],q) = qWilf(n,Patterns,q)`):
print(`qSeqS(S,K,q): the terms n=1..K for scheme S, as computed by qMiklos`):
end:



with(combinat):

#inv(pi) counts the number of inversions of pi
inv:=proc(pi) local pair, M:
 M:=0:
 for pair in choose(nops(pi), 2) do
  if pi[pair[1]]>pi[pair[2]] then M:=M+1: fi:
 od:
M:
end:

#qnorm(S,q) finds the distribution of inv over set S
qnorm:=(S,q)->add(q^inv(p), p in S):


#qWilf(n, Patterns,q) finds the distribution of inv
#over Wilf(n,Patterns)
qWilf:=(n,Patterns,q)->qnorm(Wilf(n,Patterns), q):


#qMiklos(S,pi,v,q): Given a scheme S, and a prefix permutation
#pi (that is part of the scheme) and a vector v of non-negative
#integers v (of size nops(pi)+1) outputs the polynomial in q
#enumerating (according to inv) the number of permutations whose 
#prefix reduces to pi and has the gap-vector v. For example try:
#qMiklos(SchemeFast({[1,2,3]},2,1),[],[5]);
qMiklos:=proc(S,pi,v,q) local i,lu,GapVectors,g,pi1,v1,su,j,L1:
option remember:

if nops(pi)+1<>nops(v) then
 ERROR(`Bad input: the third argument should be one longer than the sec.`):
fi:



for i from 1 to nops(S) do
 if S[i][1]=pi then
  lu:=S[i]:
  break:
 fi:
od:

if i=nops(S)+1 then
 RETURN(FAIL):
fi:

GapVectors:=lu[2]:


#Check gap vector criterion
for g in GapVectors do
 if {seq(evalb(v[j]>=g[j]),j=1..nops(v))}={true} then
  RETURN(0):
 fi:
od:

#Check if prefix is exhaustive, i.e. all letters appear in the prefix
if v=[0$nops(v)] then
  RETURN(q^inv(pi)):
fi:


if lu[3]<>{} then
 pi1:=DeleteEntries(pi,lu[3]):
  L1:=sort(convert(lu[3],list)):
  L1:={seq(pi[L1[i]],i=1..nops(L1))}:
  v1:=DeleteEntriesG(v,L1):
  RETURN( q^(qinveff(pi, v, lu[3])) *qMiklos(S,pi1,v1,q)):

fi:

su:=0:

for i from 1 to nops(pi)+1 do
 pi1:=Append1(pi,i):

 for j from 0 to v[i]-1 do
   v1:=[op(1..i-1,v),j,v[i]-1-j,op(i+1..nops(v),v)]:
   su:=su+qMiklos(S,pi1,v1,q):
 od:
od:
su:
end:



#qSeqS(S,K,q): the first K terms (starting at 1)
#computed by the enumeration scheme S, 
#for example try: SeqS(SchemeFast(6,2,{[1,2,3]}),20,q);
qSeqS:=proc(S,K,q,nice:=true) local i: 
if nice and type(q,symbol) then
   return [seq(sort(expand(qMiklos(S,[],[i],q)),q,ascending) ,i=1..K)]:
else
   return [seq(qMiklos(S,[],[i],q) ,i=1..K)]:
fi:
end:



#####
#Tools for qMiklos
####

#qinveff(p, v, T) counts the number of inversions lost
#when simultaneously removing p[t] for each t in T 
#from a permutation with prefix p and gap vector v.
qinveff := proc (p, v, T) local cast, t, M, contrib; 

cast := IV1Gap(convert(v, `+`)+nops(v)-1, nops(v)-1, v)[1];
M := 0; 

for t in T do 
 contrib := map(i-> `if`(p[t] < i, 1, -1), p[1 .. t-1]); 
 M := M+cast[p[t]]-1+convert(contrib, `+`):
od:

#any inversions occuring in the sub-permutation of rev del indices
#will be counted twice.
M - inv( [seq(p[t], t in convert(T,list) )] ):
end:


####
#End of Tools for qMiklos
####



#qSipur(L,Gvul,GvulGAP,K,q): Everything you ever wanted to know
#about all patterns of type L, that schemes of depth <=Gvul
#with gap vectors with sum <=GvulGAP
#and printing out, in the successful cases, the first K+1
#terms. For example, try:
#Sipur([3],2,20);
qSipur:=proc(L,Gvul, GvulGAP,K,q)
local gu,S,F,sch,i,j,ka,mu,sqn1, sqn2, n, meansqn, centsqn, moments:
gu:=AllPatternSets(L);
print(`There all together`, nops(gu), `different equivalence classes `):
S:={}:
F:={}:

for i from 1 to nops(gu) do
 sch:=SchemeImageFast(gu[i][1],Gvul,GvulGAP):
  if sch<>FAIL then
   print(`For the equivalence class of patterns`, gu[i]):
    mu:=sch[2]:
    ka:=max(seq(nops(mu[j][1]),j=1..nops(mu))):
   print(`the member `, sch[1], `has a scheme of depth `, ka):
   print(`here it is: `):
   print(mu):
   print(`Using  the scheme, the first, `, K, `terms are `):
   sqn1:=qSeqS(mu,K,q):
   lprint(sqn1):
   print(`with the reverse patterns and complement patterns having distributions`):
   sqn2:=[seq( normal(q^binomial(n,2) * subs(q=q^(-1), sqn1[n])), n=1..K)]:
   lprint(sqn2):
   print(``):

   print(`The number of permutations avoiding`, sch[1], `is given by`):
   print( subs(q=1, sqn1)):
   print(`The number of EVEN permutations avoiding`, sch[1], ` is given by`):
   print( (subs(q=1, sqn1) + subs(q=-1,sqn1))/2):
   print(`The number of ODD permutations avoiding`, sch[1], ` is given by`):
   print( (subs(q=1, sqn1) - subs(q=-1,sqn1))/2):
   print(`For the reverse patterns and complement patterns, we get`):
   print(`EVEN:`, (subs(q=1, sqn2) + subs(q=-1,sqn2))/2):
   print(` ODD:`, (subs(q=1, sqn2) - subs(q=-1,sqn2))/2):   
   print(``):
   print(`The average number of inversions for each n is given by`):
   meansqn:=map(p-> `if`(p<>0, subs(q=1, evalf(diff(p,q)/p)), FAIL), sqn1):   
   print(meansqn):
   print(`The standard deviation for each n is given by`):
   print(map(p-> `if`(p<>0, evalf(sqrt(subs(q=1, diff(p,q,q)/p + diff(p,q)/p - (diff(p,q)/p)^2))), FAIL), sqn1)):
   print(`The centralized moments are`):
   #centsqn is the centralized sequences
   centsqn:=[seq( normal(sqn1[n]*q^(-meansqn[n])), n=1..K)]:

   moments:=map(p->rMoments(p,q,4), centsqn):
   print(`Second: `, map(L->L[2], moments)):
   print(`Skewness: `, map(L->L[3]/(L[2]^(3/2)), moments)):
   print(`Kurtosis: `, map(L->L[4]/(L[2]^2), moments)):
   print(`end of this data`): print(``):
   
   S:=S union {gu[i]}:
 else
   F:=F union {gu[i]}:
  fi:
od:

print(`Out of a total of `, nops(gu), `cases `):
print(nops(S), `were successful and `, nops(F) , `failed `):
print(`Success Rate: `, evalf(nops(S)/nops(gu),3) ):
print(`Here are the failures `):
print(F):
end:


#####
#Tools for qSipur
#####


#rMoments(f,t,M) returns the first M centralized moments of distribution f(t)
rMoments:=proc(f,t,M) local mu, cf, factorial, r, k:
  if f=0 then return FAIL: fi:
  mu:=subs(t=1,diff(f,t)/f):
  cf:=expand(normal(t^(-mu)*f)):
  factorial:=subs(t=1,[seq( diff(cf, t$r)/cf, r=1..M)]):
  evalf[6]([seq( add(stirling2(r,k)*factorial[k] ,k=1..r), r=1..M)]):
end: