print(`ClusterGPP - A package by Andrew Baxter`):
print(`to accompany "Applying the Cluster Method to Generalized Permutation Patterns"`):
print(``):
print(`Use Help() for more information`):
print(`This version last updated 7/26/09`):
print(`Warning -- Des12, DesWt, Ris12,RisWt are now global variables.  Do not redefine them.`):

with(combinat):
with(ListTools):

#################
#
#This package contains the following:
#
#  1. Pattern functions, including some hard-coded patterns.
#     Also includes major index maj (inv is [2,0,1])
#
#
#  2. A brute-force method of computing the distribution of a pattern
#       This uses part 1.
#   
#
#  3. A recurrence-based method of computing the distribution of a pattern
#       Base-cases use part 2.
#
#    3a.  Distributions over words with a given alphabet vector.
#           Each type-(1,2) pattern has its own function, directed by Rdist.
#           These are all subsumed under RdistM in section 4.
#    3b.  Distributions over permutations of length n.
#
#
#  4.  A recurrence-based method of computing the distribution of multiple patterns
#          over words with given alphabet vector.  
#
#################


Help:=proc()
 print(`Patterns must be entered in the alphabet {0,1,2,3...},`):
 print(`where a 0 represents a dash.  Hence (a-cb) would be written`):
 print(`as [1,0,3,2]`):
 print(``):
 print(`PatternCount(pattern,word) will count for any user-entered pattern`):
 print(`e.g. PatternCount([1,3,0,2], [1,3,4,6,2,5])=5 for 13-2 and 46-5`):
 print(`Yes, you can have repeated letters in both the pattern and the word`):
 print(`e.g. PatternCount([1,2,0,2],[1,2,2,2,4,5,3,5]) returns 3.`):
 print(``):
 print(`BFdist(pattern, set, q) is the distribution of pattern over set`):
 print(`BFdist([2,0,1], permute(4), q) is the distribution of inv over S_4`):
 print(`BFdist can also handle multistatistics via lists of patterns`):
 print(`e.g. BFdist([[2,0,1],[2,1]],permute(4),[q,t]) is bistatistic (inv, des)`):
 print(`BFdist uses PatternCount and brute force computation, so beware large sets`):
 print(``):
 print(`Rdist(pattern,m,q) is the distribution of pattern over all words with`):
 print(`alphabet vector m.  The distribution is calculated recursively for`):
 print(`certain pre-programmed patterns: (a-cb), (b-ca), (c-ba), (a-ba), (b-ba),`):
 print(`their reverses and their complements.  Otherwise BFdist is used.`):
 print(`When m=[1,...,1], permutation-specific recurrences are used.`):
 print(``):
 print(`RdistM(patterns,m,q) is the distribution of the multi-statistic`):
 print(`patterns over all words with alphabet vector m.  The distribution`):
 print(`is calculated recursively for those multi-statistics which use`):
 print(`patterns which all lie within same category, as discussed in`):
 print(`the accompanying article.`):
 print(`e.g. RdistM([[2,1],[1,0,3,2]], [1,1,2,1], [t,q])`):
end:


#################
#
#  1. Pattern Functions
#
#################

#########
#
# 1a. A few 3-patterns have been pre-programmed
#
#########

a0cb := proc (p) local count, pair;
count := 0; 
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else 
  for pair in choose(nops(p)-1, 2) do 
   if p[pair[1]] < p[pair[2]+1] and p[pair[2]+1] < p[pair[2]] then count := count+1: fi: 
  od:
fi:
count: 
end proc:

b0ca := proc (p) local count, pair;
count := 0; 
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
  for pair in choose(nops(p)-1, 2) do
   if p[pair[2]+1] < p[pair[1]] and p[pair[1]] < p[pair[2]] then count := count+1 fi:
  end do:
 end if;
count:
end proc: 

c0ba := proc (p) local count, pair;
count := 0;
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else 
  for pair in choose(nops(p)-1, 2) do
   if p[pair[2]] < p[pair[1]] and p[pair[2]+1] < p[pair[2]] then count := count+1 end if:
  end do:
 end if; 
count: 
end proc:

a0bc := proc (p) local count, pair; 
count := 0;
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
 
  for pair in choose(nops(p)-1, 2) do
   if p[pair[1]] < p[pair[2]] and p[pair[2]] < p[pair[2]+1] then count := count+1 end if:
  end do:
 end if; 
count:
end proc: 

b0ac := proc (p) local count, pair;
count := 0; 
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
  for pair in choose(nops(p)-1, 2) do 
   if p[pair[1]] < p[pair[2]+1] and p[pair[2]] < p[pair[1]] then count := count+1 end if:
  end do 
 end if;
count:
end proc:

c0ab := proc (p) local count, pair;
count := 0; 
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
  for pair in choose(nops(p)-1, 2) do
   if p[pair[2]+1] < p[pair[1]] and p[pair[2]] < p[pair[2]+1] then count := count+1 end if:
  end do:
 end if;
count:
end proc:


cb0a := proc (p) local count, pair;
count := 0;
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
  for pair in choose([`$`(2 .. nops(p))], 2) do
   if p[pair[2]] < p[pair[1]] and p[pair[1]] < p[pair[1]-1] then count := count+1 end if 
  end do 
 end if;
count:
end proc:


maj := proc(p) local n, i, M; 
n := nops(p); M := 0; 
for i to n-1 do 
  if p[i+1] < p[i] then M := M+i fi:
 od:
M:
end:

b0ba := proc (p) local count, pair; 
count := 0; 
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else 
  for pair in choose(nops(p)-1, 2) do
   if p[pair[1]] = p[pair[2]] and p[pair[2]+1] < p[pair[2]] then count := count+1 end if 
  end do 
 end if;
count:
 end proc:

a0ba := proc (p) local count, pair; 
count := 0;
if nops(p) = 0 or nops(p) = 1 or nops(p) = 2 then return 0:
 else
  for pair in choose(nops(p)-1, 2) do 
   if p[pair[1]] = p[pair[2]+1] and p[pair[2]+1] < p[pair[2]] then count := count+1 end if 
  end do 
 end if;
count:
end proc:




#########
#
# 1b. PatternCount(pattern,word) will count for any user-entered pattern
#
#########


#PatternCount(pattern,p) is the number of occurences
#of pattern in word p
PatternCount:=proc(patt,p) 
local patt1, n,m, Dashes, BlockSizes, Tuples, tuple, indices, NumBlocks, M:

n:=nops(p):
if n=0 then return 0: fi:

#strip any leading or terminal 0's:
patt1:=[op({patt[1]} minus {0}), op(2..nops(patt)-1,patt), op({patt[nops(patt)]} minus {0})]:

#Find the indices of the dashes
Dashes:= [SearchAll(0,patt1)]:

m:=nops(patt1)-nops(Dashes):

if nops(Dashes)=0 then
 BlockSizes:=[nops(patt1)]:
else
 BlockSizes:=[Dashes[1]-1, seq(Dashes[i]-Dashes[i-1]-1,i=2..nops(Dashes)), nops(patt1)-Dashes[nops(Dashes)]]:
fi:

NumBlocks:=nops(BlockSizes):


Tuples:=choose(n-BlockSizes[NumBlocks]+1,NumBlocks):

M:=0:
for tuple in Tuples do
  indices:=[seq(j$j=tuple[i]..tuple[i]+BlockSizes[i]-1,i=1..nops(tuple))]:
  if nops(indices)=nops({op(indices)}) then
   if PatternCheck(patt1,p,indices) then M:=M+1: fi:
  fi:
od:
return M:
end:

######
#
# PatternCount requires PatternCheck, which requires ReduceWord
#
######

#Checks if the sub-word of word p at indices reduces to patt
PatternCheck:=proc(patt,p,indices) local subword, dashless:

subword:=[seq(p[i], i in indices)]:
dashless:=remove(n->evalb(n=0),patt):

return evalb(ReduceWord(subword)=dashless):
end:


ReduceWord := proc (L) local letters, L1, subitin;
if min(L) < 1 then L1 := L+[`$`(1-min(L), nops(L))] else L1 := L end if;
letters := sort(MakeUnique(L1)); 
subitin := [seq(letters[i] = i, i = 1 .. nops(letters))]; 
subs(subitin, L1) :
end proc:




#################
#
#  2.  Distributions
#
#################



#########
#
# 2. Brute-force Method: BFdist
#
#########


###
#Example: BFdist([2,0,1], permute(n), q) is the distribution of inv over S_n.
###


BFdist:=proc(pattern, set, q)
if type(pattern[1], list) then
 if type(q,list) and nops(q)<>nops(pattern) then ERROR(`q must be symbolic or have same length as first argument`): fi:
 add( mul(q[i]^PatternCount(pattern[i], word),i=1..nops(pattern)),   word in set):
else
 add( q^PatternCount(pattern, word),   word in set):
fi:
end:



#########
#
# 3. Single patterns via Cluster Method Recurrences
#
#########

######
# 3a.  Single patterns over words with given alphabet vector.
#
# Rdist(pattern,alphabetvector,q)
#  Checks if pattern fits one of the pre-programmed 
#  patterns and uses them if there is a match
#  Otherwise BFdist is used.
#####

Rdist:=proc(pattern,m,q)

  if pattern=[2,1]     or pattern=[1,2]     then Rdistba(m,q):
elif pattern=[1,0,3,2] or pattern=[2,3,0,1] then Rdista0cb(m,q):
elif pattern=[2,0,3,1] or pattern=[1,3,0,2] then Rdistb0ca(m,q):
elif pattern=[3,0,2,1] or pattern=[1,2,0,3] then Rdistc0ba(m,q):
elif pattern=[1,0,2,1] or pattern=[1,2,0,1] then Rdista0ba(m,q):
elif pattern=[2,0,2,1] or pattern=[1,2,0,2] then Rdistb0ba(m,q):
elif pattern=[1,0,2,3] or pattern=[3,2,0,1] then Rdista0bc(m,q):
elif pattern=[2,0,1,3] or pattern=[3,1,0,2] then Rdistb0ac(m,q):
elif pattern=[3,0,1,2] or pattern=[2,1,0,3] then Rdistc0ab(m,q):
elif pattern=[1,0,1,2] or pattern=[2,1,0,1] then Rdista0ab(m,q):
elif pattern=[2,0,1,2] or pattern=[2,1,0,2] then Rdistb0ab(m,q):
else
  print(`Using Brute Force....`):
  BFdist(pattern,permute([seq(i$m[i],i=1..nops(m))]),q):
fi:
end:

###
#All of the recursions will require the E shift operator
###

E:=proc (steps, L) local m, step;
 m := L:
 for step in steps do
   m[step] := m[step]-1
 od:
m:
end:

##
#Rchoose is a variant of choose, wherein the sets are written in descending order
#This is necessary to match the convention for T written in the paper
##
Rchoose:=(n,k)->map(Reverse,choose(n,k)):


##
# Each Rdistx0yz(m,q) computes the distribution for (x-yz) recursively
# using the recursion obtained via the cluster method.
# These are subsumed under RdistM(m,q), but are retained to demonstrate
# recurrences shown in the paper.
##


Rdistba:=proc(m1,q) local m, n, T, k: option remember:

 if {op(m)}={1} then return PRdistba(n,q): fi: 
 if member(0,m) then BFdist([2,1], permute([seq(i$m[i],i=1..nops(m))]), q):
  n:=nops(m):

  else add( add( mul( q-1 , j=1..k-1) * Rdistba(E(T,m),q) , T in Rchoose(n,k)), k=1..n):
 fi:
end:



Rdista0cb:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if {op(m)}={1} then return PRdista0cb(n,q): fi:
 if member(0, m) then BFdist([1,0,3,2], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = 1 .. T[j+1]-1)-k+j+1)-1, j = 1 .. k-1)*Rdista0cb(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
fi:
end:



Rdistb0ca:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if {op(m)}={1} then return PRdistb0ca(n,q): fi:
 if member(0, m) then BFdist([2,0,3,1], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = T[j+1]+1 .. T[j]-1)) -1, j = 1 .. k-1)*Rdistb0ca(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
fi:
end:



Rdistc0ba:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if {op(m)}={1} then return PRdistc0ba(n,q): fi:
 if member(0, m) then BFdist([3,0,2,1], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = T[j]+1 .. n))-1, j = 1 .. k-1)*Rdistc0ba(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
fi:
end:



Rdista0ba:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if max(m)=1 then return (convert(m,`+`))!: fi:
 if member(0, m) then BFdist([1,0,2,1],permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul(q^(m[T[j+1]]-1)-1, j = 1 .. k-1)*Rdista0ba(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) :
fi:
end:



Rdistb0ba:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if max(m)=1 then return (convert(m,`+`))!: fi:
 if member(0, m) then BFdist([2,0,2,1],permute([seq(i$m[i],i=1..nops(m))]), q):
  else add(add(mul(q^(m[T[j]]-1)-1, j = 1 .. k-1)*Rdistb0ba(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n): 
fi:
end:




Rdista0bc:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if {op(m)}={1} then return PRdistc0ba(n,q): fi: #This uses the complement's equivalence in permutations
 if member(0, m) then BFdist([1,0,2,3], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = 1 .. T[j+1]-1))-1, j = 1 .. k-1)*Rdista0bc(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
fi:
end:


Rdistb0ac:=proc(m,q) local n, T, k;
 n := nops(m);
 if {op(m)}={1} then return PRdistb0ca(n,q): fi: #This uses the complement's equivalence in permutations
 if member(0, m) then BFdist([2,0,1,3], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = T[j+1]+1 .. T[j]-1))-1, j = 1 .. k-1)*Rdistb0ac(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
fi:
end:




Rdistc0ab:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if {op(m)}={1} then return PRdista0cb(n,q): fi: #This uses the complement's equivalence in permutations
 if member(0, m) then BFdist([3,0,1,2], permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul( q^(add(m[i], i = T[j]+1 .. n)-j+1)-1, j = 1 .. k-1)*Rdistc0ab(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) 
 fi:
end:



Rdista0ab:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if max(m)=1 then return (convert(m,`+`))!: fi:
 if member(0, m) then BFdist([1,0,1,2],permute([seq(i$m[i],i=1..nops(m))]), q):
  else add( add( mul(q^(m[T[j+1]]-1)-1, j = 1 .. k-1)*Rdista0ab(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n) :
 fi:
end:



Rdistb0ab:=proc(m,q) local n, T, k; option remember:
 n := nops(m);
 if max(m)=1 then return (convert(m,`+`))!: fi:
 if member(0, m) then BFdist([2,0,1,2],permute([seq(i$m[i],i=1..nops(m))]), q):
  else add(add(mul(q^(m[T[j]]-1)-1, j = 1 .. k-1)*Rdistb0ab(E(T, m),q), T in Rchoose(n, k)), k = 1 .. n): 
 fi:
end:


######
#  
# 3b. Distributions of Single patterns over permutations of length n
#
######


###
#  PRdistx0yz(n,q) computes the distribution of pattern (x-yz) over symmetric group S_n
#  They use the permutation-specific recurrences, which
#  run faster than their word counterparts.
#  Each recurrence requires 2 or 3 sub-routines, arecX, brecX, and crecX for X=1,2, or 3.
###


PRdistba:=proc(n,q) local k: option remember:
 if n=0 or n=1 then 1: else add( binomial(n,k) * (q-1)^(k-1) * PRdistba(n-k,q), k=1..n): fi:
end:



PRdista0cb:=proc(n,q) local k: option remember:
 if n=0 or n=1 then 1: else add(arec1(n,k,q)*PRdista0cb(n-k,q),k=1..n): fi:
end:

   arec1:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then n:
    else arec1(n-1,k,q)+brec1(n-1,k-1,q): fi:
   end:

   brec1:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then add(q^(i-1)-1,i=1..n):
    else brec1(n-1,k,q)+(q^(n-k)-1)*brec1(n-1,k-1,q): fi:
   end:




PRdistb0ca:=proc(n,q) local k: option remember:
 if n=0 or n=1 then 1: else add(arec2(n,k,q)*PRdistb0ca(n-k,q),k=1..n): fi:
end:

   arec2:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then n:
    else arec2(n-1,k,q)+ q^(n-1)*brec2(n-1,k-1,q) - arec2(n-1,k-1,q): fi:
   end:

   brec2:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then add(q^(-i),i=1..n):
    else brec2(n-1,k,q)+ q^(-1)*brec2(n-1,k-1,q) - q^(-n)*arec2(n-1,k-1,q): fi:
   end:




PRdistc0ba:=proc(n,q) local k: option remember:
 if n=0 or n=1 then 1: else add(arec3(n,k,q)*PRdistc0ba(n-k,q),k=1..n): fi:
end:

   arec3:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then n:
    else brec3(n-1,k-1,q): fi:
   end:

   brec3:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then add( i*(q^(n-i)-1), i=1..n):
    else brec3(n-1,k,q)+ crec3(n-1,k,q) +(q^(n-1)-1)*crec3(n-1,k-1,q): fi:
   end:

   crec3:=proc(n,k,q) option remember:
    if k>=n then 0:
    elif k=1 then add(q^(n-i)-1,i=1..n):
    else crec3(n-1,k,q)+ (q^(n-1)-1) *crec3(n-1,k-1,q): fi:
   end:


######
#
#  4.  Distributions of multiple patterns over words
#
######


###
#RdistM(patterns, m, q)
#  Find the distribution of multistatistic patterns [pattern1, pattern2,...,patternsN] using the
#  indeterminates [q1, q2, ..., qN]
#  This works recursively and has been designed to be quicker than Rdist above.  
###



#Des12, DesWt, Ris12, and RisWt are global variables, necessary for RdistM:
 #Des12 are the descent-based patterns of type-(1,2), along with (ba)
 Des12:=[[2,1], [1,0,3,2], [2,0,3,1], [3,0,2,1], [1,0,2,1], [2,0,2,1]]:
 #DesWt is the descent-weights for each of the patterns in Des12 (in order)
 DesWt:=[1, 
 (T,j,k,n,m)->add(m[i], i = 1 .. T[j+1]-1)-k+j+1,
 (T,j,k,n,m)->add(m[i], i = T[j+1]+1 .. T[j]-1),
 (T,j,k,n,m)->add(m[i], i = T[j]+1 .. n),
 (T,j,k,n,m)->m[T[j+1]]-1,
 (T,j,k,n,m)->m[T[j]]-1
 ]:


 #Ris12 are the rise-based patterns of type-(1,2), along with (ab)
 Ris12:=[[1,2], [1,0,2,3], [2,0,1,3], [3,0,1,2], [1,0,1,2], [2,0,1,2]]:
 #RisWt is the rise-weights for each of the patterns in Ris12 (in order)
 RisWt:=[1, 
  (T,j,k,n,m)->add(m[i], i = 1 .. T[j+1]-1),
  (T,j,k,n,m)->add(m[i], i = T[j+1]+1 .. T[j]-1), 
  (T,j,k,n,m)->add(m[i], i = T[j]+1 .. n)-j+1, 
  (T,j,k,n,m)->m[T[j+1]]-1,
  (T,j,k,n,m)->m[T[j]]-1
 ]:


RdistM:=proc(patterns, m1, q, PattWt1:=0) local m, S,categ, patts,PattWt, T,n,k,i: option remember:

 
 if not type(patterns[1],list) then ERROR(`First argument must be a list of lists`): fi:
 if type(q,list) and nops(q)<>nops(patterns) then ERROR(`q must be symbolic or have same length as first argument`): fi:

 m:=remove(i->i=0,m1):
 if m=[] or m=[1] then return 1: fi:


if PattWt1=0 then 
   #First check if the patterns are category-consistent
    categ:=map(S->`subset`({op(patterns)},{op(S)}), [Des12, Ris12, map(Reverse, Des12), map(Reverse, Ris12)]):
    if not member(true, categ) then return FAIL, `Cross-category multistatistic.  Use BFDist`: fi:
   #Get the pattern into a multistatistic of rise-based or descent-based patterns of type-(1,2)
   #  along with a marker patts[2] which says whether des (1) or rise (2).
    patts:=[patterns, Search(true, categ)]:
    if patts[2]=3 or patts[2]=4 then patts[1]:=map(Reverse,patts[1]): patts[2]:=patts[2]-2:  fi:
   
   if patts[2]=1 then PattWt:=map(p-> DesWt[Search(p,Des12)] ,  patts[1]): 
   elif patts[2]=2 then PattWt:=map(p-> RisWt[Search(p,Ris12)] ,  patts[1]):
   else return FAIL:
   fi:
 else PattWt:=PattWt1:
 fi:

 n := nops(m):

 add( add( 
         mul(   mul(q[d]^( PattWt[d](T,j,k,n,m) ), d=1..nops(PattWt) ) -1  , j = 1 .. k-1)*RdistM(patterns, E(T, m), q, PattWt), 
       T in Rchoose(n, k)), k = 1 .. n): 

end: