with(combinat):

print("Welcome to ROOKSCHEMES, a Maple package to complement the paper"):
print("Pattern avoidance in the rook monoid by Dan Daly and Lara Pudwell"):
print("available electronically at: (insert web address)"):
print(""):
print("The main brute force procedures are:"):
print("numavoidersS(Q,n),setavoidersS(Q,n),numavoidersS2d(Q,n),setavoidersS2d(Q,n)"):
print("The main recursive procedures are:"):
print("BuildScheme(Pats,maxw,depth),SchemeTo0Triangle(S,N),SchemeToSeq(S,N),classify(L,maxw,depth,N)"):
print("To learn more about one of these procedures, type Help(functionname);"):
print("To see a list of all procedures in this package type Help2();"):


Help:=proc()
if args=NULL then 
print("The main brute force procedures are:"):
print("numavoidersS(Q,n),setavoidersS(Q,n),numavoidersS2d(Q,n),setavoidersS2d(Q,n)"):
print("The main recursive procedures are:"):
print("BuildScheme(Pats,maxw,depth),SchemeTo0Triangle(S,N),SchemeToSeq(S,N),classify(L,maxw,depth,N)"):
print("To learn more about one of these procedures, type Help(functionname);"):

elif args[1]=numavoidersS then
print("numavoidersS(Q,n): inputs a pattern set Q and a positive integer n"):
print("uses brute force computation to compute the number of rook placements on an nxn board"):
print("avoiding all patterns in set Q"):
print("For example, try, numavoidersS({[0,0,1]},4)"):
print("Because this procedure is brute force, values of n larger than 7 are not recommended."):

elif args[1]=setavoidersS then
print("setavoidersS(Q,n): inputs a pattern set Q and a positive integer n"):
print("uses brute force computation to compute the set of rook placements on an nxn board"):
print("avoiding all patterns in set Q"):
print("For example, try, setavoidersS({[0,0,1]},4)"):
print("Because this procedure is brute force, values of n larger than 7 are not recommended."):

elif args[1]=numavoidersS2d then
print("numavoidersS2d(Q,n): inputs a pattern set Q and a positive integer n"):
print("uses brute force computation to compute the number of rook placements on an nxn board"):
print("avoiding all patterns in set Q in the 2d sense"):
print("For example, try, numavoidersS2d({[0,0,1]},4)"):
print("Because this procedure is brute force, values of n larger than 5 are not recommended."):

elif args[1]=setavoidersS2d then
print("setavoidersS2d(Q,n): inputs a pattern set Q and a positive integer n"):
print("uses brute force computation to compute the set of rook placements on an nxn board"):
print("avoiding all patterns in set Q in the 2d sense"):
print("For example, try, setavoidersS2d({[0,0,1]},4)"):
print("Because this procedure is brute force, values of n larger than 5 are not recommended."):

elif args[1]=BuildScheme then
print("BuildScheme(Pats,maxw,depth): inputs a pattern set Pats and a positive integers maxw and depth"):
print("uses a divide-and-conquer algorithm to find a recurrence for rook placements avoiding Pats with max"):
print("gap vector weight maxw and depth at most depth"):
print("For example, try, BuildScheme({[1,0],[2,1]},2,2)"):

elif args[1]=SchemeTo0Triangle then
print("SchemeTo0Triangle(S,N): inputs an enumeration scheme S and a positive integer N"):
print("returns a sequence of lists of length N where the j+1st entry of the ith list is the"):
print("number of rook placements on an ixi board with exactly j 0s that obey scheme S"):
print("For example, try, SchemeTo0Triangle(BuildScheme({[1,0],[2,1]},2,2),10)"):


elif args[1]=SchemeToSeq then
print("SchemeToSeq(S,N): inputs an enumeration scheme S and a positive integer N"):
print("returns a sequence of integers of length N where the ith entry is the total number"):
print("of rook placements on an ixi board that obey scheme S"):
print("For example, try, SchemeToSeq(BuildScheme({[1,0],[2,1]},2,2),10)"):



elif args[1]=classify then
print("classify(L,maxw,depth,N): inputs a list of positive integers L, and positive integers maxw, depth, and N"):
print("and considers all pattern sets with one pattern of each length in L"):
print("If a pattern set (or one of its symmetries) has an enumeration scheme with max gap weight maxw and"):
print("depth at most depth, then it reports the scheme along with the number of rook placements obeying the scheme"):
print("of length at most N.  These enumerations are reported both as total number of enumerations and divided"):
print("according to number of 0s"):
print("For example, try, classify([2],2,2,5)"):

else
print("Sorry, there is Help for a procedure called ", args[1],"."):
print("Type Help(); to see a list of available major procedures."):

fi:

end:

Help2:=proc()
print("Functions for brute force counting:"):
print("R(n), makewords(zeros,nums), makeword(p,zeros), red(p), pcontainsq(p,q), pavoidsq(p,q)"):
print("pavoidsQset(p,Q), split0s(L), match0s(L1,L2), numavoidersS(Q,n), setavoidersS(Q,n)"):
print("vertical0s(L), horizontal0s(L), redvert(L,zeros), zeropartition(L,zeros), listleq(L1,L2)"):
print("pcontainsq2d(p,q), pavoidsq2d(p,q), pavoidsQ2dset(p,Q), numavoidersS2d(Q,n), setavoidersS2d(Q,n)"):
print(""):
print("Functions for schemes:"):
print("redwith0s(pre), getChildren(pre), GoodAndBadChildren(currpres,pre,Pats,d), BuildScheme(Pats,maxw,depth)"):
print("forb(lets,gap,rest,z), makesets(L), numnonzero(L,S), SchemeWZerosToNum(S,pre,lets, n,z,rest)"):
print("SchemeAndPreToNum(S,pre,lets, n), SchemeTo0Triangle(S,N), SchemeToSeq(S,N)"):
print(""):
print("Functions for reversibly deletable elements:"):
print("prefs(pat), occ(pre,pat), endings(alph,size), fracletters(minl,maxl,n), scenarios(pre,pat)"):
print("posscenarios(pre,pat,pos), setposscenarios(pre,Pats,pos), del(L,S), isrevdel1(pre,Pats,l,gaps)"):
print("obeysgap(scen,g), obeysgaps(scen,gaps), GapScenarios(pre,Pats,pos,gaps)"):
print("isrevdel(pre,Pats,S,gaps), findRevDel(currpres,pre,Pats,gaps)"):
print(""):
print("Functions for gap vectors:"):
print("IsGapVector(pre,Pats,v), parts(w,n), weight(v), AddGapPair(v,L), ImpliedGaps(v,w) "):
print("AllGapCands(w,n), FindGaps(pre,Pats,w)"):
print(""):
print("Functions for generating data files:"):
print("allpats(n), allpatsSet(L), redundant(S), rev(L), comp(L), allsym(L), redSet(S)"):
print("revSet(S), compSet(S), allsymSet(S), allClasses(L), classify(L,maxw,depth,N)"):
end:

#R(n): the set of rook monoid elements of length n
R:=proc(n) local S,k,zeros,nums,z,nu:
S:={op(permute(n))}:
for k from 1 to n do
zeros:=choose([$1..n],k):
nums:=choose([$1..n],n-k):

for z in zeros do
for nu in nums do
S:=S union makewords({op(z)},nu):
od:
od:
od:

S:

end:


#makewords(zeros,nums): inputs a set of positive integers "zeros"
#and a set of positive integers "nums" and outputs the set of all rook monoid elements
#that use the values in nums as entries but have zeros in the positions in "zeros"
makewords:=proc(zeros,nums) local P,S,p:
P:=permute(nums):
S:={}:
for p in P do
S:=S union {makeword(p,zeros)}:
od:

S:

end:


#makeword(p,zeros): inputs a permutation p and a set of positive integers "zeros" and
#outputs the word formed by inserting 0s in the locations "zeros" into permutation p.
makeword:=proc(p,zeros) local L,i,p2:
L:=[]:
p2:=p:
for i from 1 to nops(p)+nops(zeros) do
if member(i,zeros) then
L:=[op(L),0]:
else
L:=[op(L),p2[1]]:
p2:=[op(2..nops(p2),p2)]:
fi:
od:
L:
end:


#red(p): the reduction of word p.
red:=proc(p) local L,p2:
L:=[$1..nops(p)]:
p2:=sort(p):
subs({seq(p2[i]=L[i],i=1..nops(p))},p):
end:

#pcontainsq(p,q): returns true if p contains q in the 1d sense.
pcontainsq:=proc(p,q) local pos,pos1,ptemp:
pos:=choose([$1..nops(p)],nops(q)):

for pos1 in pos do
ptemp:=[seq(p[pos1[i]],i=1..nops(pos1))]:


if evalb(red(split0s(ptemp))=split0s(q)) and match0s(ptemp,q) then return true: fi:
od:
false:
end:

#pavoidsq(p,q): returns true if p avoids q in the 1d sense.
pavoidsq:=proc(p,q)
not pcontainsq(p,q):
end:

#pavoidsQset(p,Q): returns true if p avoids Q in the 1d sense.
pavoidsQset:=proc(p,Q)
evalb({seq(pavoidsq(p,q),q in Q)}={true}):
end:

#split0s(L): inputs a string L and returns L with all 0s deleted.
split0s:=proc(L) local i,L2:
L2:=[]:
for i from 1 to nops(L) do
if L[i]<>0 then L2:=[op(L2),L[i]]: fi:
od:
L2:
end:

#match0s(L1,L2): inputs lists L1 and L2 of equal length and returns true if L1
#and L2 have 0s in identical positions
match0s:=proc(L1,L2) local Z1,Z2,i:
if nops(L1)<>nops(L2) then return FAIL: fi:

Z1:={}: Z2:={}:

for i from 1 to nops(L1) do
if L1[i]=0 then Z1:=Z1 union {i}: fi:
if L2[i]=0 then Z2:=Z2 union {i}: fi:
od:

evalb(Z1=Z2):
end:

#numavoidersS(Q,n): inputs a set of reduced rook patterns and a positive integer n
#and returns the number of members of R(n) that avoid all patterns
#in Q in the 1d sense
numavoidersS:=proc(Q,n) local S,c,s:
S:=R(n):
c:=0:
for s in S do
if evalb({seq(pavoidsq(s,q), q in Q)}={true}) then c:=c+1: fi:
od:
c:
end:

#setavoidersS(Q,n): inputs a set of reduced rook patterns and a positive integer n
#and returns the set of members of R(n) that avoid all patterns
#in Q in the 1d sense
setavoidersS:=proc(Q,n) local S,R2,s:
S:=R(n):
R2:={}:
for s in S do
if evalb({seq(pavoidsq(s,q), q in Q)}={true}) then R2:=R2 union {s}: fi:
od:
R2:
end:

#vertical0s(L): inputs a rook monoid element and outputs the set of rows 
#that have no rook
vertical0s:=proc(L) local S,n:
n:=nops(L):
S:={$1..n} minus {op(split0s(L))}:
end:

#horizontal0s(L): inputs a rook monoid element and outputs the set of columns 
#that have no rook
horizontal0s:=proc(L) local S,i:
S:={}:
for i from 1 to nops(L) do
if L[i]=0 then S:=S union {i}: fi:
od:
S:
end:

#redvert(L,zeros): inputs a word and a set of vertical 0s and reduces the word, respecting the
#vertical 0s.  For example, try redvert([2,4],{3}) which returns [[1,3],{2}]
redvert:=proc(L,zeros) local p:
p:=red([op(L),op(zeros)]):
return [[seq(p[i],i=1..nops(L))], {seq(p[j],j=nops(L)+1..nops(p))}]:
end:


#zeropartition(L,zeros): inputs a word and a set of vertical 0s and returns a vector of length
#nops(L)+1 where entry i of the vector tells the number of vertical 0s between element i-1 and element i
#in the word.
zeropartition:=proc(L,zeros) local L2,v,z,zeros2,curr,curr2,j, flag,a:
L2:=[op(sort(L)),infinity]:

v:=[0$(nops(L)+1)]:
zeros2:=zeros:
j:=1:
curr:=0:
curr2:=L2[1]:

while nops(zeros2)>0 do
a:=min(op(zeros2)):

flag:=false:

while flag=false do
if a>curr and a<curr2 then 
v[j]:=v[j]+1: flag:=true:

else 
j:=j+1: 

curr:=L2[j-1]: curr2:=L2[j]:

fi:
od:

zeros2:=zeros2 minus {a}:


od:

v:

end:


#listleq(L1,L2): inputs two lists of equal length and returns true of component-wise L1 is less
#than or equal to L2
listleq:=proc(L1,L2)
if nops(L1)<>nops(L2) then return false:

else return evalb({seq(evalb(L1[i]<=L2[i]),i=1..nops(L1))}={true}): fi:

end:

#pcontainsq2d(p,q): returns true if p contains q in the 2d sense
pcontainsq2d:=proc(p,q) local pos,pos1,ptemp:
pos:=choose([$1..nops(p)],nops(q)):

for pos1 in pos do
ptemp:=[seq(p[pos1[i]],i=1..nops(pos1))]:


#need to check (1) horizontal 0s match, and (2) horizontal permutation matches
#(3) vertical 0s match



if evalb(red(split0s(ptemp))=red(split0s(q))) and match0s(ptemp,q) and 
listleq(zeropartition(split0s(q),vertical0s(q)),zeropartition(op(redvert(split0s(ptemp),vertical0s(p))))) then return true: fi:
od:
false:
end:

#pavoidsq2d(p,q): returns true if p avoids q in the 2d sense
pavoidsq2d:=proc(p,q)
not pcontainsq2d(p,q):
end:

#pavoidsQ2dset(p,Q): returns true if p avoids Q in the 2d sense
pavoidsQ2dset:=proc(p,Q)
evalb({seq(pavoidsQ2d(p,q),q in Q)}={true}):
end:

#numavoidersS2d(Q,n): inputs a set of reduced rook patterns and a positive integer n
#and returns the number of members of R(n) that avoid all patterns
#in Q in the 2d sense
numavoidersS2d:=proc(Q,n) local S,c,s:
S:=R(n):
c:=0:
for s in S do
if evalb({seq(pavoidsq2d(s,q), q in Q)}={true}) then c:=c+1: fi:
od:
c:
end:

#setavoidersS2d(Q,n): inputs a set of reduced rook patterns and a positive integer n
#and returns the set of members of R(n) that avoid all patterns
#in Q in the 2d sense
setavoidersS2d:=proc(Q,n) local S,R2,s:
S:=R(n):
R2:={}:
for s in S do
if evalb({seq(pavoidsq2d(s,q), q in Q)}={true}) then R2:=R2 union {s}: fi:
od:
R2:
end:

#####################################################
#
# functions for making schemes
#
#####################################################

#redwith0s(pre): inputs a prefix and reduces the non-zero entries.
redwith0s:=proc(pre) local H, ptemp,v,i:
H:=horizontal0s(pre):
ptemp:=red(split0s(pre)):
v:=[]:
for i from 1 to nops(pre) do:
if member(i,H) then v:=[op(v),0]:
else v:=[op(v),ptemp[1]]: ptemp:=[seq(ptemp[j],j=2..nops(ptemp))]:
fi:
od:
v:
end:

#getChildren(pre): inputs a rook prefix and outputs the set of children prefixes of length 
#nops(pre)+1
getChildren:=proc(pre) local n:
n:=nops(split0s(pre)):
{[op(pre),0],seq(redwith0s([op(pre),i+0.5]),i=0..n)}:
end:

#GoodAndBadChildren(currpres,pre,Pats,d): inputs a prefix pre and a set of patterns Pats and returns
#an ordered triple [G,B,C] of sets
#G is a set of pairs of [prefix,reversibly deletable positions]
#B is a set of prefixes with no reversible deletable positions
#C is a set of prefixes containing a pattern from Pats
#all children of pre appear as prefixes in G or in B or in C.
GoodAndBadChildren:=proc(currpres,pre,Pats,d) local G, B,temp,t, C,g:
G:={}:
B:={}:
C:={}:
temp:=getChildren(pre):

for t in temp do
#if not pavoidsQset(t,Pats) then C:=C union {t}:
#elif findRevDel(t,Pats)={} then B:=B union {t}: 
#else G:=G union {[t,findRevDel(t,Pats), FindGaps(t,Pats,d)]}:
#fi:

g:=FindGaps(t,Pats,d):

if findRevDel(currpres,t,Pats,g)={} then B:=B union {t}: 
G:=G union {[t,{}, {}]}:
else G:=G union {[t,findRevDel(currpres,t,Pats,g), g]}:
fi:
od:

[G,B]:
end:

#BuildScheme(Pats,maxw,depth): inputs a set of forbidden patterns Pats, a 
#max gap vector pair weight, and a scheme depth
#and returns an enumeration scheme if possible, or FAIL and the partial scheme if
#no scheme was found with depth less than or equal to the given input
BuildScheme:=proc(Pats,maxw, depth) local S, todo, d, t,temptodo, dead,par,currpres:

S:={[[],{},{}]}:
todo:={[]}:
dead:={}:

for d from 0 to depth do

if nops(todo)=0 then return [op(S)]: fi:
for t in todo do
temptodo:={}:
currpres:={seq(s[1],s in S)}:
par:=GoodAndBadChildren(currpres,t,Pats,maxw):
S:=S union par[1]:
temptodo:=temptodo union par[2]:
#dead:=dead union par[3]:

od:
todo:=todo minus {t} union temptodo:

od:

return [FAIL,S,todo]:
end:


#forb(lets,gap,rest,z) given the current literal prefix lets, the desired
#set of nonzero digits rest, and the number of zeros z, returns all sets of digits
#that may not appear according to the gap vectors in gap
#the set returned is a set of pairs of the form [set of nonzero digits, number of zeros]
forb:=proc(lets,gap,rest,z) local letssort,gapallow,a,i,forbset,nonzero,g,temp,t:
letssort:=sort(split0s(lets)):
#sorted version of lets without 0s

if nops(letssort)>0 then
gapallow:=[{}$(nops(letssort)+1)]:

for a in rest do
if a<letssort[1] then gapallow[1]:=gapallow[1] union {a}:
elif a>letssort[nops(letssort)] then gapallow[nops(letssort)+1]:=gapallow[nops(letssort)+1] union {a}:
else
for i from 1 to nops(letssort)-1 do
if a>letssort[i] and a<letssort[i+1] then gapallow[i+1]:=gapallow[i+1] union {a}: fi:
od:
fi:
od:
#gapallow is now the allowed digits sorted by which prefix digits they fit between
else
gapallow:=[rest]:
fi:

forbset:={}:

for g in gap do
if g[2]<=z then
nonzero:=g[1]:

temp:=makesets([seq([gapallow[i],nonzero[i]],i=1..nops(nonzero))]):

for t in temp do
forbset:=forbset union {[t,g[2]]}:
od:

fi:
od:
forbset:
end:

#makesets(L) is a subprocedure of forb
makesets:=proc(L) local S,Stemp,S2,i,s,s2:
Stemp:=[seq(choose(L[i][1],L[i][2]),i=1..nops(L))]:
S:=Stemp[1]:
S2:={}:
for i from 2 to nops(L) do
for s in S do
for s2 in Stemp[i] do
S2:=S2 union {s union s2}:
od:
od:
S:=S2:
S2:={}:
od:

S:
end:


#numnonzero(L,S): given a list L and a set S returns the number of non-zero entries of L in positions given by S
numnonzero:=proc(L,S) local c,i:
c:=0:
for i from 1 to nops(S) do:
if L[S[i]]<>0 then c:=c+1: fi:
od:
c:
end:

#SchemeWZerosToNum(S,pre,lets, n,z,rest): n is width, z in number of zeros after prefix, allow is the set of allowed digits
SchemeWZerosToNum:=proc(S,pre,lets,n,z,rest) local schemepres,schemeduo,s,gap, rd,allow2,g,F,f,tot: option remember:

#print("start with ", pre, lets, n, z, rest):

if z+nops(rest)+nops(pre)<>n then ERROR("bad input: nops(pre)+z + nops(rest) should equal n.  We have pre=", pre, "z=", z, "rest=",rest, "and n=", n): fi:

if z>n-nops(pre) then return 0: fi:
#if we are supposed to have more zeros than remaining digits
#in the rook element, then return 0

schemepres:={seq(s[1],s in S)}:
schemeduo:=[]:

if not member(pre,schemepres) then ERROR("The prefix ", pre, " is not in the scheme!"): fi:

for s in S do
if pre=s[1] then schemeduo:=s: fi:
od:

#schemeduo is the tuple of the scheme corresponding to prefix pre


gap:=schemeduo[3]:


F:=forb(lets,gap,rest,z):

#print("forb is ", F):

for f in F do
if f[1] subset rest and z>=f[2] then return 0: fi:
od:
#this loop checks the digits and number of zeros
#we're expecting against gap vectors.  If [rest,z] violates
#a gap vector, there are no rook elements to report


if nops(gap)>0 and weight(gap[1])=0 then 
return 0:
fi: 
if nops(pre)=n then return 1:
fi:
#two base cases: if we get the all 0s vector, we know the prefix contains a forbidden pattern
#also if the prefix length equals the length of the entire word, then lets is
#the only possible rook placement for the given input


#we've made it past all the base cases!  now, two options
#(a) there are reversibly deletable elements
#(b) there are no rd elements and we must further partition by prefixes

rd:=schemeduo[2]:

if rd<>{} then


#need prefix with letters deleted
#lets with letters deleted
#n reduced by nops(rd)
#rest and z remain unchanged
#print("recurse with ", redwith0s(del(pre,rd)),del(lets,rd), n-nops(rd),z,rest):

return SchemeWZerosToNum(S,redwith0s(del(pre,rd)),del(lets,rd), n-nops(rd),z,rest):

else

#in this case, prefix and lets get one digit longer
#using digits from rest union 0
#n stays the same
#either rest or z reduce depending on which digit we used


if z>0 then
tot:=SchemeWZerosToNum(S,[op(pre),0],[op(lets),0],n,z-1,rest):
else
tot:=0:
fi:

return tot + add(SchemeWZerosToNum(S,redwith0s([op(lets),r]),[op(lets),r],n,z,rest minus {r}),r in rest):

fi:

end:


#SchemeAndPreToNum(S,pre,lets, n): n is width, z is number of zeros
#allow is the set of digits that appear after the prefix
SchemeAndPreToNum:=proc(S,pre,lets,n) local z: option remember:
add(binomial(n,n-z)*SchemeWZerosToNum(S,pre,lets, n,z,{$1..(n-z)}),z=0..n):
end:

#SchemeTo0Triangle(S,N): inputs a scheme and a positive integer n 
#and outputs a list of lists where list i has length i rook elements with 0,... i 0s 
#obeying the scheme
SchemeTo0Triangle:=proc(S,N) local i,j,r:
[seq([seq(binomial(i,i-j)*SchemeWZerosToNum(S,[],[],i,j,{$1..(i-j)}),j=0..i)],i=1..N)]:
end:


#SchemeToSeq(S,N): inputs a scheme
#and outputs the number of elements of R(n)
#that obey the scheme for 1<=n<=N
SchemeToSeq:=proc(S,N)
[seq(SchemeAndPreToNum(S,[],[],n,n),n=1..N)]:
end:



########################################################
#
# functions for finding reversibly deletable elements
#
########################################################

#prefs(pat): inputs a pattern pat and outputs all prefixes of pat
prefs:=proc(pat) local S, i:
S:={}:
for i from 1 to nops(pat) do
S:=S union {redwith0s([seq(pat[j],j=1..i)])}:
od:
S:
end:

#occ(pre,pat): inputs a prefix pre, and a pattern pat and outputs all partial occurences of
#pat in pre
occ:=proc(pre,pat) local S, precand, p, pos, p2:
S:={}:
precand:=prefs(pat):

for p in precand do
pos:=choose({$1..nops(pre)},nops(p)):

for p2 in pos do
if evalb(redwith0s([seq(pre[p2[i]],i=1..nops(p2))])=p) then S:=S union {[[seq(pre[p2[i]],i=1..nops(p2))],p2]}:fi:
od:
od:

S:
end:

#endings(alph,size): inputs a set of letters alph and a positive integer size
#and outputs all permutations of letters chosen from alph of length size
endings:=proc(alph,size) local sets, s:
sets:=choose(alph,size):

{seq(op(permute(s)),s in sets)}:

end:

#fracletters(minl,maxl,n)
fracletters:=proc(minl,maxl,n) local f:
f:=(maxl-minl)/(n+1):
{seq(minl+f*i,i=1..n)}:
end:

#scenarios(pre,pat): inputs a prefix pre, and a pattern pat and outputs all possible ways 
#pat can be partially contained in pre
scenarios:=proc(pre,pat) local alph,occurences,S,o,newalph,e,E,n:
alph:=sort([op(split0s(pre))]):



occurences:=occ(pre,pat):

S:={}:
for o in occurences do

n:=nops(pat)-nops(o[1]):
if nops(alph)>=1 then
newalph:=[op(fracletters(0,alph[1],n)),seq(op(fracletters(alph[j],alph[j+1],n)),j=1..nops(alph)-1),op(fracletters(alph[nops(alph)],alph[nops(alph)]+1,n)),0$nops(horizontal0s(pat))]:
else
newalph:=[0$nops(horizontal0s(pat)),$1..nops(pat)-1]:
fi:



E:=endings(newalph,nops(pat)-nops(o[1])):

for e in E do
if redwith0s([op(o[1]),op(e)])=pat then

S:=S union {[[op(pre),op(e)],op(o),e]}:
fi:
od:
od:
S:

end:

#posscenarios(pre,pat,pos): inputs a prefix pre and a pattern pat, and a set of positions pos
#returns all ways for pre to partially contain pat using all positions in pos
posscenarios:=proc(pre,pat,pos) local S, s, S2:
S:=scenarios(pre,pat):

S2:={}:
for s in S do
if pos subset s[3] then S2:=S2 union {s}: fi:
od:
S2:
end:

#setposscenarios(pre,Pats,pos): same as posscenarios, but can handle a *set* of patterns instead of just one
setposscenarios:=proc(pre,Pats,pos)
{seq(op(posscenarios(pre,p,pos)),p in Pats)}:
end:

#del(L,S): inputs a list L and a set of indices S, and removes L[s], s in S from L
del:=proc(L,S) local i, L2:
L2:=[]:
for i from 1 to nops(L) do
if not member(i,S) then L2:=[op(L2),L[i]]: fi:
od:
L2:
end:

#isrevdel1(pre,Pats,l,gaps): inputs a prefix pre, a set of patterns Pats, and a position l that is a
#reversibly deletable candidate, returns true if l is reversibly deletable
isrevdel1:=proc(pre,Pats,l,gaps) local scen,predelwords,delwords:
scen:=GapScenarios(pre,Pats,{l},gaps):
predelwords:={seq(s[1],s in scen)}:

delwords:={seq(del(s,{l}),s in predelwords)}:

({seq(pavoidsQset(s,Pats),s in delwords)} subset {false}):


end:

#obeysgap(scen,g)
obeysgap:=proc(scen,g) local L,i,v,pre,j:

L:=sort(scen[4]):
pre:=sort([seq(scen[1][i],i=1..nops(scen[1])-nops(scen[4]))]):

while nops(pre)>0 and pre[1]=0 do
pre:=[seq(pre[i],i=2..nops(pre))]:
od:

pre:=[op(pre),infinity]:
v:=[[0$nops(g[1])],0]:


while nops(L)>0 and L[1]=0 do
v[2]:=v[2]+1:
L:=[seq(L[i],i=2..nops(L))]:
od:


#print(L,pre, v):

for i from 1 to nops(L) do

if L[i]<pre[1] then v[1][1]:=v[1][1]+1:
else
for j from 1 to nops(pre)-1 do
if L[i]>pre[j] and L[i]<pre[j+1] then
v[1][j+1]:=v[1][j+1]+1:
fi:
od:
fi:
od:
v:

#print(v):

not evalb({seq(evalb(v[1][i]>=g[1][i]),i=1..nops(v[1])), evalb(v[2]>=g[2])}={true}):

end:

#obeysgaps(scen,gaps)
obeysgaps:=proc(scen,gaps)
if gaps={} then return true: fi:

evalb({seq(obeysgap(scen,g),g in gaps)}={true}):
end:

#GapScenarios(pre,Pats,pos,gaps)
GapScenarios:=proc(pre,Pats,pos,gaps) local S,R,s,g:
S:=setposscenarios(pre,Pats,pos):
R:={}:

#for g in gaps do
#if max(pre)>0 and weight(g)>0 and weight(g)=g[2] then return S: fi:
#od:

#print(S):

for s in S do
#print(s,obeysgaps(s,gaps)):
if obeysgaps(s,gaps) then R:=R union {s}: fi:
od:

R:
end:

#isrevdel(pre,Pats,S,gaps): inputs a prefix pre, a set of patterns Pats, and a set S of reversibly
#deletable candidates, returns true if the set is reversibly deletable
isrevdel:=proc(pre,Pats,S,gaps) local scen,predelwords,delwords:

if nops(S)=1 then return isrevdel1(pre,Pats,op(S),gaps): fi:

if member(false, {seq(isrevdel1(pre,Pats,S[i], gaps),i=1..nops(S))}) then return false: fi:

if nops(S)>1 and nops(GapScenarios(pre,Pats,S,gaps))=0 and {seq(nops(GapScenarios(pre,Pats,{s},gaps)),s in S)}<>{0} then return false: fi:

scen:=GapScenarios(pre,Pats,S,gaps):



predelwords:={seq(s[1],s in scen)}:

delwords:={seq(del(s,S),s in predelwords)}:

({seq(pavoidsQset(s,Pats),s in delwords)} subset {false}):


end:

#findRevDel(currpres,pre,Pats,gaps): inputs a prefix pre and a set of patterns Pats and outputs 
#the largest possible set of reversibly deletable positions of pre
findRevDel:=proc(currpres,pre,Pats,gaps) local pos, curr,p:
pos:=powerset({$1..nops(pre)}) minus {{}}:

curr:={}:
for p in pos do
if isrevdel(pre,Pats,p,gaps) and nops(p)>nops(curr) and member(redwith0s(del(pre,p)),currpres) then curr:=p: 
fi:
od:

curr:
end:


########################################################
#
# functions for finding gap pairs
# (here, a gap pair for a prefix with n non-zero elements is a vector of length n+1 (like a normal
# gap vector) and an integer number of 0s to consider)
#
########################################################



#IsGapVector(pre,Pats,v)
IsGapVector:=proc(pre,Pats,v) local n,newlets,i,poss,p:
n:=numnonzero(pre,{$1..nops(pre)}):

if nops(v[1])<>n+1 then return false: fi:

newlets:=fracletters(0,1,v[1][1]):
for i from 1 to nops(v[1])-1 do
newlets:=newlets union fracletters(i,i+1,v[1][i+1]):
od:
newlets:=[op(newlets),0$v[2]]:

poss:=permute(newlets):

for p in poss do
if pavoidsQset([op(pre),op(p)],Pats) then return false: fi:
od:

true:

end:

#parts(w,n): all ways to partition w into n parts
parts:=proc(w,n) local S,i,R,s:
R:={}:
if n=0 and w>0 then return {}:
elif n=0 and w=0 then return {[]}:
elif n=1 then return {[w]}:
else
for i from 0 to w do
S:=parts(w-i,n-1):
R:=R union {seq([i,op(s)], s in S)}:
od:
fi:
R:
end:

#weight(v): weight of gap pair v
weight:=proc(v)
v[2]+add(v[1][i],i=1..nops(v[1])):
end:


#AddGapPair(v,L)
AddGapPair:=proc(v,L)

[[seq(v[1][i]+L[i],i=1..nops(v[1]))],v[2]+L[nops(L)]]:

end:

#ImpliedGaps(v,w): all gap pairs of weight w that are implied by the
#pair v
ImpliedGaps:=proc(v,w) local P,p:

P:=parts(w-weight(v),nops(v[1])+1):

{seq(AddGapPair(v,p), p in P)}:

end:

#AllGapCands(w,n)
AllGapCands:=proc(w,n) local poss:
poss:=parts(w,n+1):

{seq([[seq(p[i],i=1..n)],p[n+1]],p in poss)}:

end:


#FindGaps(pre,Pats,w)
FindGaps:=proc(pre,Pats,w) local w1,n,Cand,c,B:
n:=numnonzero(pre,{$1..nops(pre)}):
B:={}:

for w1 from 0 to w do
Cand:=AllGapCands(w1,n+1) minus {seq(op(ImpliedGaps(b,w1)),b in B)}:

for c in Cand do
if IsGapVector(pre,Pats,c) then 

#if not (weight(c)>0 and weight(c)=c[2]) then
B:=B union {c}: 
#fi:
fi:
od:
od:


B:
end:

#####################################################
#
# functions for generating data files
#
#####################################################

#allpats(n): generates the set of all
#one-dimensional rook patterns of length n
allpats:=proc(n) local S,S2,s,z,w:
S:=R(n):
S2:={}:
for s in S do
z:=horizontal0s(s):
w:=split0s(s):
S2:=S2 union {makeword(red(w),z)}:
od:
S2:
end:

#allpatsSet(L): generates the set of all
#sets of one-dimensional rook patterns
#with one pattern of each length in list L
allpatsSet:=proc(L) local s, S1, S2,R,s1,s2,L2,R2,r:
if nops(L)=1 then return {seq({s},s in allpats(L[1]))}: fi:
S1:=allpats(L[1]):
S2:=allpatsSet([seq(L[i],i=2..nops(L))]):
R:={}:
for s1 in S1 do
for s2 in S2 do
if not member(s1,s2) then
R:=R union {{s1,op(s2)}}:
fi:
od:
od:
R2:=R:

for r in R2 do
if redundant(r) then
R:=R minus {r}:
fi:
od:
R:
end:

redundant:=proc(S) local s,s2:
for s in S do
for s2 in S minus {s} do
if pcontainsq(s,s2) then return true: fi:
od:
od:
false:
end:

#rev(L): the reverse of list L
rev:=proc(L)
[seq(L[nops(L)+1-i],i=1..nops(L))]:
end:

#comp(L): the complement of rook pattern L
comp:=proc(L) local L2,i:
L2:=[]:
for i from 1 to nops(L) do
if L[i]=0 then L2:=[op(L2),0]:
else L2:=[op(L2),nops(L)+1-L[i]]:
fi:
od:
L2:
end:


#allsym(L): all trivially Wilf-equivalent patterns to one-D pattern L
allsym:=proc(L)
{L,comp(rev(L)),comp(L), rev(L)}:
end:

#redSet(S): the reverse of every pattern in set S
redSet:=proc(S) local L:
{seq(redwith0s(L),L in S)}:
end:

#revSet(S): the reverse of every pattern in set S
revSet:=proc(S) local L:
{seq(rev(L),L in S)}:
end:

#compSet(S): the complement of every pattern in set S
compSet:=proc(S) local L:
{seq(comp(L),L in S)}:
end:


#allsymSet(S): all trivially Wilf-equivalent pattern sets 
#to one-D pattern set S
allsymSet:=proc(S)
{S,redSet(compSet(revSet(S))), redSet(compSet(S)), redSet(revSet(S))}:
end:

#allClasses(L): same as allpatsSet(L), but divided into trivial symmetry classes
allClasses:=proc(L) local S,R,T:
S:=allpatsSet(L):
R:={}:
while S<>{} do
T:=allsymSet(S[1]) intersect S:
S:=S minus T:
R:=R union {T}:
od:
R:
end:

#classify(L,maxw,depth,N): inputs a list of pattern lengths,
#positive integers maxw, depth, and N, and looks for one
#pattern set in each trivial equivalence class of allpatsSet(L)
#that has an enumeration scheme with parameters maxw and depth
#If successful, returns the first N terms of total number of 
#rook placements, as well as the first N rows of the 
#triangle of rook placements of length i, counted according to number
#of empty rows
classify:=proc(L,maxw,depth,N) local pats,p,p2,flag,F,S,T,i:
pats:=allClasses(L):
print("There are ", nops(allpatsSet(L)), " pattern sets of type ", L):
print("They fall into ", nops(pats), " trivial pattern classes."):
F:={}:

for p in pats do
#p is a trivial equivalence class, so p[i]
#is an actual pattern set.
flag:=false:
for p2 in p do
if not flag then
 S:=BuildScheme(p2,maxw,depth):
 if S[1]<>FAIL then
  flag:=true:
  print("---------------------------------------------"):
  print("For the equivalence class ", p, " the pattern set ", p2, " has the scheme ", S): 
  print("First ", N, " terms of the counting sequence are ", SchemeToSeq(S,N)):
  print("First ", N, " rows of the 0s triangle are: "):
T:=SchemeTo0Triangle(S,N):
for i from 1 to nops(T) do
print(T[i]):
od:

if SchemeToSeq(S,5)<>[seq(numavoidersS(p2,n),n=1..5)] then
print("DOES NOT MATCH BRUTE FORCE DATA.  NEED TO DEBUG!"):
fi:

 fi:
fi:
od:
if flag=false then F:=F union {p}: fi:
od:

if nops(F)>0 then
  print("---------------------------------------------"):
print("Here are the failures: ", F):
fi:

  print("---------------------------------------------"):
print("Overall, we found schemes for ", nops(pats)-nops(F), " out of ", nops(pats), " trivial classes, for a success rate of ", evalf((nops(pats)-nops(F))/nops(pats))):

end: