with(combinat):

Help:=proc()

if nargs=0 then 
print("Welcome to bVATTER, a Maple package of Lara Pudwell"):
print("for counting permutations avoiding barred patterns."):
print(""):
print("The main functions are : SchemeSlow(Patterns,Depth),"):
print("SchemeFast(Patterns,Depth,GapDepth), and SeqS(Scheme,Integer)"):
print("You may also see Example1(); and Example2(); for the "):
print("best known examples of barred patterns in the literature."):
print(""):
print("notation: a barred permutation is represented as:"):
print("[[p1,a1],[p2,a2],...,[pn,an]]"):
print("where the permutation is p1...pn, and ai=1 if the entry is barred, ai=0 otherwise"):
print(""):
print("For example, try SchemeSlow({[[1,1],[3,0],[2,0]]},2);"):
print("versus SchemeFast({[[1,1],[3,0],[2,0]]},2,1);"):


elif args[1]=Append1 then
print("Append1(pi,i): Given a permutation pi, and an integer"):
print("i between 1 and nops(pi)+1, constructs the permutation"):
print("of length nops(pi)+1 whose last entry is i, and such that"):
print("reduced form of the first nops(pi) letters is pi"):
print("For example, try: Append1([2,1,3],2);"):

elif args[1]=redu then
print("redu(perm): Given a permutation on a set of positive integers"):
print("finds its reduced form"):

elif args[1]=IV1 then
print("IV1(n,k) all increasing vectors of length"):
print("k of the form 1<=i_1< ...<i_k <=n"):

elif args[1]=Wilf then
print("Wilf(n,P): inputs a set of barred patterns P and an integer n"):
print("outputs the set of permutations of length n avoiding all patterns in P."):
print(""):
print("notation: a barred permutation is represented as:"):
print("[[p1,a1],[p2,a2],...,[pn,an]]"):
print("where the permutation is p1...pn, and ai=1 if the entry is barred, ai=0 otherwise"):
print(""):
print("for example, try: Wilf(5,{[[3,0],[5,1],[2,0],[4,0],[1,0]],[[2,0],[3,0],[4,0],[1,0]]});"):

elif args[1]=f then
print("f(p,S): inputs a single barred permutation pattern p"):
print("and a set of permutations S and returns all elements of S"):
print("that avoid p"):
print("for example, try: f([[3,0],[5,1],[2,0],[4,0],[1,0]],{op(permute(5))});"):

elif args[1]=reduce then
print("reduce(w): returns the reduced form of permutation w"):
print("for example, try reduce([4,6,5,2]);"):

elif args[1]=collapse then
print("collapse(S,pos): removes all positions in pos from every list in set S"):
print("for example, try: collapse({[1,2,3,4],[5,6,7,8],[9,10,11,12]},{2,4});"):

elif args[1]=isgood then
print("isgood(p,q): returns true if p avoids q,"):
print("false if p contains q, and contains a set of all lists of positions where a bad pattern occurs"):
print("for example, try: isgood([1,2,3],[1,2]]); versus isgood([1,2,3],[2,1]]);"):

elif args[1]=WilfTable then
print("WilfTable(n,Patterns,k): Given a pos. integer n, a set of"):
print("patterns Patterns and a pos. integer k<=n, creates a"):
print("table T such that for every k-permutation  pi and"):
print("every increasing vector of integers v"):
print("T[pi,v] is the subset of Wilf(n,Patterns) of perms whose"):
print("first k entries reduce to pi and consist of the elements of v"):
print("For example, try: op(WilfTable(5,{[[1,1],[3,0],[2,0]]},2));"):

elif args[1]=WilfPre then
print("WilfPre(n,Patterns,pi,v): The set of members of Wilf(n,Patterns)"):
print("whose first k members reduce to pi and whose elements are"):
print("v (v is increasing). "):
print("For example, try: WilfPre(5,{[[1,1],[3,0],[2,0]]},[1,2],[1,5]);"):

elif args[1]=Delete1 then
print("Delete1(v,r): deleting the r^th entry of the permutation v"):
print("and reducing"):
print("For example, try: Delete1([1,4,2,3],2);"):

elif args[1]=DeleteSet1 then
print("DeleteSet1(S,r): deletes the r^th entry and reduces from"):
print("each of the permutations in S"):
print("For example, try: DeleteSet1({[1,4,2,3]},2);"):

elif args[1]=PD1 then
print("PD1(n,Patterns,pi,v,r): The `partial derivative` of"):
print("WilfPre(n,Patterns,pi,v) w.r.t. to the r^th entry"):
print("of pi, i.e. the symmetric difference of"):
print("WilfPre(n-1,Patterns,redu(pi w/o pi[r]), redu(v,w/o v[pi[r]]))"):
print("and the image of WilfPre(n,Patterns,pi,v) under the map"):
print("delete the r^th entry and reduce"):
print("For example, try: PD1(5,{[[1,1],[3,0],[2,0]]},[1,2],[1,5],2);"):

elif args[1]=IsRevDel11 then
print("IsRevDel11(n0,Patterns,pi,r): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(n0,Patterns)?"):
print("For example, try:"):
print("IsRevDel11(6,{[[1,1],[3,0],[2,0]]},[2,1],1);"):

elif args[1]=IsRevDel1 then
print("IsRevDel1(n0,Patterns,pi,r): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(i,Patterns), i<=n0?"):
print("For example, try:"):
print("IsRevDel1(6,{[[1,1],[3,0],[2,0]]},[2,1],1);"):

elif args[1]=IV1Gap then
print("IV1Gap(n,k,G): the subset of IV1(n,k) of vectors v"):
print("satisfying v[1]-0>=G[1]+1, v[2]-v[1]>=G[2]+1 , ..."):
print("v[k]-v[k-1]>=G[k]+1, n+1-v[k]>=G[k+1]+1"):
print("For example, try:IV1Gap(5,3,[1,0,0]);"):

elif args[1]=IV1GapC then
print("IV1GapC(n,k,S): all the vectors v of IV1(n,k) that"):
print("do NOT satisfy v[1]-0>=G[1]+1, v[2]-v[1]>=G[2]+1 , ..."):
print("v[k]-v[k-1]>=G[k]+1, n+1-v[k]>=G[k+1]+1"):
print("for each G in S"):

elif args[1]=IsRevDel11Gap then
print("IsRevDel11Gap(n0,Patterns,pi,r,G): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(n0,Patterns) satisying the gap condition"):
print("given by the set of gap vectors G?"):
print("For example, try:"):
print("IsRevDel11Gap(6,{[[1,1],[3,0],[2,0]]},[1,2],1,{[0,1]});"):

elif args[1]=IsRevDel11Gap then
print("IsRevDel1Gap(n0,Patterns,pi,r,G): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(i,Patterns), i<=n0?"):
print("For example, try:"):
print("IsRevDel1Gap(6,{[[1,1],[3,0],[2,0]]},[2,1],1,{[0,1]});"):

elif args[1]=hofkhi then
print("hofkhi(pi): inverse of pi"):

elif args[1]=Hofkhi then
print("Hofkhi(S): inverse of every permutation in set S"):

elif args[1]=SetRevDel then
print("SetRevDel(n0,Patterns,pi): the set of reversely deleteable"):
print("entries for the Wilf class Wilf(n0,Patterns) for the"):
print("prefix permutation pi. For example, try:"):
print("SetRevDel(6,{[[1,1],[3,0],[2,0]]},[1,2]);"):

elif args[1]=SetRevDelGap then
print("SetRevDelGap(n0,Patterns,pi,G): the set of reversely deleteable"):
print("entries for the Wilf class Wilf(n0,Patterns) for the"):
print("prefix permutation pi and gap vector G. For example, try:"):
print("SetRevDelGap(6,{[[1,1],[3,0],[2,0]]},[1,2],{[0,0,1]);"):

elif args[1]=EmptiesD1 then
print("EmptiesD1(n0,Patterns,pi): Given an integer n0 and"):
print("a set of patterns Patterns and a prefix permutation"):
print("pi  (of size k, say)"):
print("finds the differences"):
print("[v[1]-1,v[2]-v[1]-1, ..., v[k]-v[k-1]-1,n-v[k]]"):
print("for all vectors that for which WilfPre(n0,Patterns,pi,v)"):
print("is empty"):
print("For example, try: EmptiesD1(5,{[[1,1],[3,0],[2,0]]},[1,2]);"):

elif args[1]=NewEmptiesD then
print("NewEmptiesD(n0,Patterns,pi): Given an integer n0 and"):
print("a set of patterns Patterns and a prefix permutation"):
print("pi  (of size k, say)"):
print("finds the differences"):
print("[v[1]-1,v[2]-v[1]-1, ..., v[k]-v[k-1]-1,n-v[k]]"):
print("for all vectors that for which WilfPre(n0,Patterns,pi,v)"):
print("is empty and that are not implied by gap vectors for smaller n"):
print("For example, try: NewEmptiesD(5,{[[1,1],[3,0],[2,0]]},[1,2]);"):

elif args[1]=GapVectors then
print("GapVectors(n0,Patterns,pi): Given an integer n0 and"):
print("a set of patterns Patterns and a prefix permutation"):
print("pi  (of size k, say)"):
print("finds all the gap vectors by investigating permutation"):
print("up to length n0"):
print("For example, try: GapVectors(7,{[[1,1],[3,0],[2,0]]},[1,2]);"):

elif args[1]=Yeladim then
print("Yeladim(pi): all the children of pi"):

elif args[1]=GaBchildren then
print("GaBchildren(n0,Patterns,pi): the good and bad children"):
print("of a prefix permutation pi w.r.t. to the set of"):
print("patterns Patterns, using empirical investigations of"):
print("permutations up to size n0. For example, try:"):
print("GaBchildren(7,{[[1,1],[3,0],[2,0]]},[1]):"):

elif args[1]=DeleteEntries then
print("DeleteEntries(pi,S): deletes all the entries of"):
print("pi[i] for i in S and reduces"):

elif args[1]=Bdok then
print("Bdok(S): checks whether the tentative scheme S is OK"):

elif args[1]=SimplifyScheme then
print("SimplifyScheme(S): simplifies the scheme S"):
print("to avoid, if necessary, the third component with cardinality larger"):
print("than 1"):

elif args[1]=Scheme1 then
print("Scheme1(n0,Gvul,Patterns): Tries to find a scheme"):
print("of depth<=Gvul for the Wilf class avoiding"):
print("the patterns in Patterns using empirical investigations"):
print("of perms of length<=n0."):
print("if it fails it returns FAIL, followed by the partial"):
print("scheme and the leftovers. For example, try:"):
print("Scheme1(6,2,{[[1,1],[3,0],[2,0]]});"):

elif args[1]=Shrink1 then
print("Shrink1(v,S): Given a vector v and a list S from {1, ..., nops(v)-1}"):
print("returns the vector obtained by removing the barriers from the"):
print("fences of S"):
print("For example,try: Shrink1([4,1,5],[1]);"):

elif args[1]=DeleteEntriesV then
print("DeleteEntriesV(v,S): given an increasing vector v and"):
print("a set S of entries, deletes the entries indexed by S"):
print("and reduces"):

elif args[1]=PD1set then
print("PD1set(n,Patterns,pi,v,S): The `partial derivative` of"):
print("WilfPre(n,Patterns,pi,v) w.r.t. to the entries"):
print("of S, i.e. the symmetric difference of"):
print("WilfPre(n-nops(S),Patterns,redu(pi w/o S), redu(v,w/o v[pi[S]]))"):
print("and the image of WilfPre(n,Patterns,pi,v) under the map"):
print("delete the entries of S and reduce"):
print("For example, try:"):
print("PD1set(5,{[[1,1],[3,0],[2,0]]},[2,1],[4,5],{2});"):

elif args[1]=IsRevDel11set then
print("IsRevDel11set(n0,Patterns,pi,S): are the  entries  of S in the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(n0,Patterns)?"):
print("For example, try:"):
print("IsRevDel11set(6,{[[1,1],[3,0],[2,0]]},[2,1],{1});"):

elif args[1]=IsRevDel1set then
print("IsRevDel1set(n0,Patterns,pi,S): are the entries of S in the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(i,Patterns), i<=n0?"):
print("For example, try:"):
print("IsRevDel1set(6,{[[1,1],[3,0],[2,0]]},[2,1],{1});"):

elif args[1]=IsRevDel11GapSet then
print("IsRevDel11GapSet(n0,Patterns,pi,S,G): Are the entries of S"):
print("for the prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(n0,Patterns) satisying the gap condition"):
print("given by the set of gap vectors G?"):
print("For example, try:"):
print("IsRevDel11GapSet(6,{[[1,1],[3,0],[2,0]]},[1,2],{1},{[0,1]});"):

elif args[1]=IsRevDel1GapSet then
print("IsRevDel1GapSet(n0,Patterns,pi,S,G): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class Wilf(i,Patterns), i<=n0?"):
print("For example, try:"):
print("IsRevDel1GapSet(6,{[[1,1],[3,0],[2,0]]},[2,1],{1},{[0,1]});"):

elif args[1]=IsIndeedScheme then
print("IsIndeedScheme(n0,Patterns,S):checks whether the scheme S"):
print("is OK for example, try:"):
print("IsIndeedScheme(7,{[[1,1],[3,0],[2,0]]},Scheme(7,3,{[[1,1],[3,0],[2,0]]}));"):

elif args[1]=DeleteEntriesG then
print("DeleteEntriesG(g,S): shrinks the vector g by removing the"):
print("barriers in the places indicated by S"):
print("For example, try: DeleteEntries([2,3,2],{2});"):

elif args[1]=Miklos then
print("Miklos(S,pi,v): Given a scheme S, and a prefix permutation"):
print("pi (that is part of the scheme) and a vector v of non-negative"):
print("integers v (of size nops(pi)+1) outputs the number"):
print("of permutations whose prefix reduces to pi and has"):
print("the gap-vector v. For example try:"):
print("Miklos(Scheme(6,2,{[[1,1],[3,0],[2,0]]}),[],[5]);"):

elif args[1]=SeqS then
print("SeqS(S,K): the first K+1 terms computed by the"):
print("enumeration scheme S, for example try:"):
print("SeqS(Scheme(6,2,{[[1,1],[3,0],[2,0]]}),20);"):

elif args[1]=SchemeSlow then
print("SchemeSlow(Patterns,Gvul): Finds a scheme for the Wilf class"):
print("of restricted permutations avoiding the set of patterns"):
print("Patterns, of depth<=Gvul. If it fails it returns"):
print("FAIL followed by the partial scheme of depth Gvul"):
print("followed by those prefix permutations that are not"):
print("reducible (that have yet to be explored)."):
print("For example, try:"):
print("SchemeSlow({[[1,1],[3,0],[2,0]]},2);"):

elif args[1]=SchemeFast then
print("SchemeFast(Patterns,Gvul,GvulGap): Finds a scheme for the Wilf class"):
print("of restricted permutations avoiding the set of patterns"):
print("Patterns, of depth<=Gvul and with gap vectors"):
print("of sum <=GvulGap . If it fails it returns"):
print("FAIL followed by the partial scheme of depth Gvul"):
print("followed by those prefix permutations that are not"):
print("reducible (that have yet to be explored)."):
print("For example, try:"):
print("SchemeFast({[[1,1],[3,0],[2,0]]},2,1);"):

else 
print(args[1], " is not a valid function"):
fi:
end:

Example1:=proc()
print("The number of 2-stack sortable permutations of length n is known to be:"):
print("2*(3*i)!/(i+1)!/(2*i+1)!"):
print("Computing this for i=1..8 gives: 1, 2, 6, 22, 91, 408, 1938, 9614"):
print("The number of 2-stack sortable permutations of length n is also known to be"):
print("exactly those which avoid [[3,0],[5,1],[2,0],[4,0],[1,0]] and [[2,0],[3,0],[4,0],[1,0]]"):
print("Compute this with seq(nops(Wilf(i, {[[3,0],[5,1],[2,0],[4,0],[1,0]],[[2,0],[3,0],[4,0],[1,0]]})),i=1..8);"):
print("to get the same sequences a before"):
end:

Example2:=proc()
print("The number of forrest-like permutations of length n is known to have OGF:"):
print("F:=((1-x)*(1-4*x-2*x^2)-(1-5*x)*sqrt(1-4*x))/(2*(1-5*x+2*x^2-x^3))"):
print("Computing the first 8 coefficients gives: 1, 2, 6, 22, 89, 379, 1661, 7405"):
print("The number of forrest-like permutations of length n is also known to be"):
print("exactly those which avoid [[2,0],[1,0],[3,1],[5,0],[4,0]] and [[1,0],[3,0],[2,0],[4,0]]"):
print("Compute this with seq(nops(Wilf(i,{[[2,0],[1,0],[3,1],[5,0],[4,0]],[[1,0],[3,0],[2,0],[4,0]]})),i=1..8);"):
print("to get the same sequence as before"):
end:

#Append1(pi,i): Given a permutation pi, and an integer
#i between 1 and nops(pi)+1, constructs the permutation
#of length nops(pi)+1 whose last entry is i, and such that
#reduced form of the first nops(pi) letters is pi
#For example, try: Append1([2,1,3],2);
Append1:=proc(pi,i): redu([op(pi),i-1/2]): end:

#redu(perm): Given a permutation on a set of positive integers
#finds its reduced form
redu:=proc(perm) local gu,i,ra:
gu:=sort(perm):
for i from 1 to nops(gu) do ra[gu[i]]:=i: od:
[seq(ra[perm[i]],i=1..nops(perm))]:
end:

#IV1(n,k) all increasing vectors of length
#k of the form 1<=i_1< ...<i_k <=n
IV1:=proc(n,k) local i: convert(choose([seq(i,i=1..n)],k),set):end:

#Wilf(n,P): inputs a set of barred patterns P and an integer n
#outputs the set of permutations of length n avoiding all patterns in P.
#
#notation: a barred permutation is represented as:
#[[p1,a1],[p2,a2],...,[pn,an]]
#where the permutation is p1...pn, and ai=1 if the entry is barred, ai=0 otherwise
#
#for example, try: WestAv({[[3,0],[5,1],[2,0],[4,0],[1,0]],[[2,0],[3,0],[4,0],[1,0]]},5);
Wilf:=proc(n,P) local S, i, S2:
S:=permute(n): 

for i from 1 to nops(P) do
S2:=f(P[i],S):
S:=S2:
od:

S:
end:

#f(p,S): inputs a single barred permutation pattern p
#and a set of permutations S and returns all elements of S
#that avoid p
#
#for example, try: f([[3,0],[5,1],[2,0],[4,0],[1,0]],{op(permute(5))});
f:=proc(p,S) local p1, pb, ppos,i,Sav,Ssmall,Sbig:
p1:=[]: pb:=[]: ppos:={}:
#every copy of p1 must be part of a copy of pb
for i from 1 to nops(p) do:
if p[i][2]=0 then
p1:=[op(p1),p[i][1]]:
else
ppos:=ppos union {i}:
fi:
pb:=[op(pb),p[i][1]]:
od:
Sav:={}:

for i from 1 to nops(S) do
Ssmall:=isgood(S[i],reduce(p1)):
if evalb(Ssmall[1]=true) then Sav:=Sav union {S[i]}:
elif ppos<>{} then
Ssmall:=Ssmall[2]:
Sbig:=collapse(isgood(S[i],reduce(pb))[2],ppos):
if evalb(Ssmall subset Sbig) then Sav:=Sav union {S[i]}: fi:
fi:
od:

Sav:
end:

#reduce(w): returns the reduced form of permutation w
#for example, try reduce([4,6,5,2]);
reduce:=proc(w) local s:
s:=sort(w):
subs({seq(s[i]=i,i=1..nops(w))},w):
end:

#collapse(S,pos):
#removes all positions in pos from every list in set S
#for example, try: collapse({[1,2,3,4],[5,6,7,8],[9,10,11,12]},{2,4});
collapse:=proc(S,pos) local S2, i, w, b:
S2:={}:
for i from 1 to nops(S) do
 w:=[]:
 for b from 1 to nops(S[1]) do
 if not member(b,pos) then w:=[op(w),S[i][b]]: fi:
 od:
 S2:=S2 union {w}:
od:
S2:
end:

#isgood(p,q): returns true if p avoids q,
#false if p contains q, and contains a set of all lists of positions where a bad pattern occurs
#
#for example, try: isgood([1,2,3],[2,1]); versus isgood([1,2,3],[1,2]);
isgood:=proc(p,q) local S,i,j,w,places:
if nops(q)>nops(p) then return [true,{}]: fi:

S:={op(choose([seq(i,i=1..nops(p))],nops(q)))}:
places:={}:
for i from 1 to nops(S) do
w:=[seq(p[S[i][j]],j=1..nops(S[i]))]:
if evalb(reduce(w)=q) then places:=places union {S[i]}: fi:
od:
if places<>{} then return [false,places]: fi:

[true,{}]:
end:

#WilfTable(n,Patterns,k): Given a pos. integer n, a set of
#patterns Patterns and a pos. integer k<=n, creates a
#table T such that for every k-permutation  pi and
#every increasing vector of integers v
#T[pi,v] is the subset of Wilf(n,Patterns) of perms whose
#first k entries reduce to pi and consist of the elements of v
#For example, try:
#op(WilfTable(5,{[1,2,3]},2));

WilfTable:=proc(n,Patterns,k) local T,gu1,gu2,i1,i2,gu,lu,sig:
option remember:
gu1:=Wilf(k,{});
gu2:=IV1(n,k):

for i1 from 1 to nops(gu1) do
 for i2 from 1 to nops(gu2) do
  T[gu1[i1],gu2[i2]]:={}:
 od:
od:

gu:=Wilf(n,Patterns):
for i1 from 1 to nops(gu) do
 sig:=gu[i1]:
 lu:=[op(1..k,sig)]:
 T[redu(lu),sort(lu)]:= T[redu(lu),sort(lu)] union {sig}:
od:
T:
end:

#WilfPre(n,Patterns,pi,v): The set of members of Wilf(n,Patterns)
#whose first k members reduce to pi and whose elements are
#v (v is increasing). For example, try:
#WilfPre(5,{[1,2,3]},[2,1],[4,5]);
WilfPre:=proc(n,Patterns,pi,v) local T:

if nops(pi)<>nops(v) or sort(v)<>v then
 print(`Bad output`):
 RETURN(FAIL):
fi:

T:=WilfTable(n,Patterns,nops(pi)) :
T[pi,v]:
end:

#Delete1(v,r): deleting the r^th entry of the permutation v
#and reducing
#For example, try: Delete1([1,4,2,3],2);
Delete1:=proc(v,r): redu([op(1..r-1,v),op(r+1..nops(v),v)]):end:

#DeleteSet1(S,r): deletes the r^th entry and reduces from
#each of the permutations in S
#For example, try: DeleteSet1({[1,4,2,3]},2);
DeleteSet1:=proc(S,r) local i:
{seq(Delete1(S[i],r),i=1..nops(S))}:
end:

#PD1(n,Patterns,pi,v,r): The "partial derivative" of
#WilfPre(n,Patterns,pi,v) w.r.t. to the r^th entry
#of pi, i.e. the symmetric difference of
#WilfPre(n-1,Patterns,redu(pi w/o pi[r]), redu(v,w/o v[pi[r]]))
#and the image of WilfPre(n,Patterns,pi,v) under the map
#delete the r^th entry and reduce
#For example, try: PD1(5,{[[1,1],[3,0],[2,0]]},[1,2],[1,5],2);
PD1:=proc(n,Patterns,pi,v,r) local gu,gu1,pi1,v1,i,gua:
gu:=WilfPre(n,Patterns,pi,v):
gua:=gu:
gu:=DeleteSet1(gu,r):
 pi1:=redu([op(1..r-1,pi),op(r+1..nops(pi),pi)]):

 v1:=[op(1..pi[r]-1,v),seq(v[i]-1,i=pi[r]+1..nops(v))]:

 gu1:=WilfPre(n-1,Patterns,pi1,v1):
(gu1 minus gu) union (gu minus gu1):
end:


#IsRevDel11(n0,Patterns,pi,r): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(n0,Patterns)?
#For example, try:
#IsRevDel11(6,{[1,2,3]},[2,1],1);
IsRevDel11:=proc(n0,Patterns,pi,r) local gu,v:
gu:=IV1(n0,nops(pi)):
evalb({seq(evalb(PD1(n0,Patterns,pi,v,r)={}),v in gu)}={true}):
end:


#IsRevDel1(n0,Patterns,pi,r): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(i,Patterns), i<=n0?
#For example, try:
#IsRevDel1(6,{[1,2,3]},[2,1],1);
IsRevDel1:=proc(n0,Patterns,pi,r) local i:
evalb({seq(IsRevDel11(i,Patterns,pi,r),i=nops(pi)..n0)}={true}):
end:

#IV1Gap(n,k,G): the subset of IV1(n,k) of vectors v
#satisfying v[1]-0>=G[1]+1, v[2]-v[1]>=G[2]+1 , ...
#v[k]-v[k-1]>=G[k]+1, n+1-v[k]>=G[k+1]+1
#For example, try:IV1Gap(5,3,[1,0,0]);
IV1Gap:=proc(n,k,G) local mu,i,v:
mu:=IV1(n-convert(G,`+`),k):
{seq([seq(v[i]+convert([op(1..i,G)],`+`),i=1..nops(v))],v in mu)}:
end:

#IV1GapC(n,k,S): all the vectors v of IV1(n,k) that
#do NOT satisfy v[1]-0>=G[1]+1, v[2]-v[1]>=G[2]+1 , ...
#v[k]-v[k-1]>=G[k]+1, n+1-v[k]>=G[k+1]+1
#for each G in S
IV1GapC:=proc(n,k,S) local G:
IV1(n,k) minus {seq(op(IV1Gap(n,k,G)),G in S)}:end:

#IsRevDel11Gap(n0,Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(n0,Patterns) satisying the gap condition
#given by the set of gap vectors G?
#For example, try:
#IsRevDel11Gap(6,{[1,2,3]},[1,2],1,{[0,1]});
IsRevDel11Gap:=proc(n0,Patterns,pi,r,G) local gu,v,i:
if G<>{} and {seq(nops(G[i])-nops(pi),i=1..nops(G))}<>{1} then
 ERROR(`Bad input`):
fi:

gu:=IV1GapC(n0,nops(pi),G):

if gu={} then
 RETURN(true):
fi:
evalb({seq(evalb(PD1(n0,Patterns,pi,v,r)={}),v in gu)}={true}):
end:

#IsRevDel1Gap(n0,Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(i,Patterns), i<=n0?
#For example, try:
#IsRevDel1Gap(6,{[1,2,3]},[2,1],1,{[0,1]});
IsRevDel1Gap:=proc(n0,Patterns,pi,r,G) local i:
evalb({seq(IsRevDel11Gap(i,Patterns,pi,r,G),i=nops(pi)..n0)}={true}):
end:

#hofkhi=inverse
hofkhi:=proc(pi)  local i, T:
for  i from 1 to nops(pi) do
T[pi[i]]:=i:
od:
[seq(T[i],i=1..nops(pi))]:
end:

#inverse of every permutation in set S
Hofkhi:=proc(S)  local i:
 {seq(hofkhi(S[i]),i=1..nops(S))}:
end:

#SetRevDel(n0,Patterns,pi): the set of reversely deleteable
#entries for the Wilf class Wilf(n0,Patterns) for the
#prefix permutation pi. For example, try:
#SetRevDel(6,{[1,2,3]},[1,2]);
SetRevDel:=proc(n0,Patterns,pi) local gu,r:
gu:={}:
for r from 1 to nops(pi) do
 if IsRevDel1(n0,Patterns,pi,r) then
   gu:=gu union {r}:
 fi:
od:
gu:
end:

#SetRevDelGap(n0,Patterns,pi,G): the set of reversely deleteable
#entries for the Wilf class Wilf(n0,Patterns) for the
#prefix permutation pi and gap vector G. For example, try:
#SetRevDelGap(6,{[1,2,3]},[1,2],{[0,0,1]);
SetRevDelGap:=proc(n0,Patterns,pi,G) local gu,r:
gu:={}:
for r from 1 to nops(pi) do
 if IsRevDel1Gap(n0,Patterns,pi,r,G) then
   gu:=gu union {r}:
 fi:
od:
gu:
end:


#EmptiesD1(n0,Patterns,pi): Given an integer n0 and
#a set of patterns Patterns and a prefix permutation
#pi  (of size k, say)
#finds the differences
#[v[1]-1,v[2]-v[1]-1, ..., v[k]-v[k-1]-1,n-v[k]]
#for all vectors that for which WilfPre(n0,Patterns,pi,v)
#is empty
#For example, try: EmptiesD1(5,{[1,2,3]},[1,2]);
EmptiesD1:=proc(n0,Patterns,pi)
local mu,gu,k,v,i:
option remember:

if pi=[] then
  RETURN({}):
fi:

k:=nops(pi):
mu:=IV1(n0,k):

gu:={}:
for v in mu do
 if WilfPre(n0,Patterns,pi,v)={} then
  gu:=gu union {[v[1]-1,seq(v[i]-v[i-1]-1,i=2..k),n0-v[k]]}:
 fi:
od:

gu:

end:


#NewEmptiesD(n0,Patterns,pi): Given an integer n0 and
#a set of patterns Patterns and a prefix permutation
#pi  (of size k, say)
#finds the differences
#[v[1]-1,v[2]-v[1]-1, ..., v[k]-v[k-1]-1,n-v[k]]
#for all vectors that for which WilfPre(n0,Patterns,pi,v)
#is empty and that are not implied by gap vectors for smaller n
#For example, try: NewEmptiesD(5,{[1,2,3]},[1,2]);
NewEmptiesD:=proc(n0,Patterns,pi)
local gu1,gu2,v,i:
option remember:

gu1:=EmptiesD1(n0-1,Patterns,pi):

gu1:={seq(seq([op(1..i-1,v),v[i]+1,op(i+1..nops(v),v)],
i=1..nops(v)),v in gu1)}:

gu2:=EmptiesD1(n0,Patterns,pi):

if gu1 minus gu2<>{} then
# print(`something may be wrong`):
fi:

gu2 minus gu1:

end:

#GapVectors(n0,Patterns,pi): Given an integer n0 and
#a set of patterns Patterns and a prefix permutation
#pi  (of size k, say)
#finds all the gap vectors by investigating permutation
#up to length n0
#For example, try: GapVectors(7,{[1,2,3]},[1,2]);
GapVectors:=proc(n0,Patterns,pi) local i:
{seq(op(NewEmptiesD(i,Patterns,pi)),i=1..n0)}:
end:


#Yeladim(pi): all the children of pi
Yeladim:=proc(pi) local i:[seq(Append1(pi,i),i=1..nops(pi)+1)]:end:

#GaBchildren(n0,Patterns,pi): the good and bad children
#of a prefix permutation pi w.r.t. to the set of
#patterns Patterns, using empirical investigations of
#permutations up to size n0. For example, try:
#GaBchildren(7,{[1,2,3]},[1]):
GaBchildren:=proc(n0,Patterns,pi) local gu,G,B,pi1,Gaps1,rd1:
gu:=Yeladim(pi):
G:={}: B:={}:

for pi1 in gu do
 Gaps1:=GapVectors(n0,Patterns,pi1):

   if Gaps1={} then
    rd1:=SetRevDel(n0,Patterns,pi1):
   else
    rd1:=SetRevDelGap(n0,Patterns,pi1,Gaps1):
   fi:

    if rd1<>{} then
       G:=G union {[pi1,Gaps1,rd1]}:
    else
       B:=B union {[pi1,Gaps1,rd1]}:
    fi:
od:

G,B:

end:

#DeleteEntries(pi,S): deletes all the entries of
#pi[i] for i in S and reduces
DeleteEntries:=proc(pi,S) local k,i,S1:
k:=nops(pi):

if S minus {seq(i,i=1..k)}<>{} then
 ERROR(`bad input`):
fi:

S1:= {seq(i,i=1..k)} minus S:
S1:=sort(convert(S1,list)):

redu([seq(pi[S1[i]],i=1..nops(S1))]):
end:

#Bdok(S): checks whether the tentative scheme S is OK
Bdok:=proc(S) local Expa,Redu,i,pi,pi1:
Expa:={}; Redu:={}:




for i from 1 to nops(S) do
if S[i][3]={} then
  Expa:=Expa union {S[i][1]}:
else
  Redu:=Redu union {S[i][1]}:
fi:
od:


for  i from 1 to nops(S) do
 pi:=S[i][1]:
  if S[i][3]={} then
    if {op(Yeladim(pi))} minus (Redu union Expa)<>{} then
     print(pi, `does not have all its children in the scheme`):
     RETURN(false):
    fi:

 else

  pi1:=DeleteEntries(pi,S[i][3]):
   if not member(pi1, Expa union Redu) then
      print(`the reduction of`, pi, `which is`, pi1, `does not belong`):
      RETURN(false):
   fi:

 fi:
od:
true:

end:

#SimplifyScheme(S): simplifies the scheme S
#to avoid, if necessary, the third component with cardinality larger
#than 1
SimplifyScheme:=proc(S)
local Prefs,i,S1,ku,j:
 Prefs:={seq(S[i][1],i=1..nops(S))}:
S1:={}:
 for i from 1 to nops(S) do
   ku:=S[i][3]:
   if nops(ku)<2 then
    S1:=S1 union {S[i]}:
   else
    for j from 1 to nops(ku) do
     if member(DeleteEntries( S[i][1],{ku[j]}),Prefs) then
      S1:=S1 union {[S[i][1],S[i][2],{ku[j]}]} :
      break:
     fi:
    od:
     if j=nops(ku)+1 then
         S1:=S1 union {S[i]}:
      fi:
    fi:
 od:
  S1:
end:

#Scheme1(n0,Gvul,Patterns): Tries to find a scheme
#of depth<=Gvul for the Wilf class avoiding
#the patterns in Patterns using empirical investigations
#of perms of length<=n0.
#if it fails it returns FAIL, followed by the partial
#scheme and the leftovers. For example, try:
#Scheme1(6,2,{[1,2,3]});
Scheme1:=proc(n0,Gvul,Patterns) local LeftToDo,S,pi,j,i,gu:
LeftToDo:={[]}:

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

for i from 0 to Gvul while LeftToDo<>{} do

  for pi in LeftToDo do
     gu:=GaBchildren(n0,Patterns,pi) :
     S:=[op(S), op(gu[1]),op(gu[2])]:
     LeftToDo:=LeftToDo minus {pi} union
       {seq(gu[2][j][1],j=1..nops(gu[2]))}:
  od:

if LeftToDo={} then
    S:=SimplifyScheme(S):
if IsIndeedScheme(n0,Patterns,S) then
   RETURN(S):
else
#print(`Not closed under deletion`):
RETURN(FAIL,S):
fi:
fi:
od:

FAIL,S,LeftToDo:
end:

#Shrink1(v,S): Given a vector v and a list S from {1, ..., nops(v)-1}
#returns the vector obtained by removing the barriers from the
#fences of S
#For example,try: Shrink1([4,1,5],[1]);
Shrink1:=proc(v,S) local i,k,v1,S1,i1:
k:=nops(v)-1:
if not type(v,list) or not type(S,list) or
convert(S,set) minus {seq(i,i=1..k)}<>{} then
 ERROR(`bad input`):
fi:

i:=S[1]:
v1:=[op(1..i-1,v),v[i]+v[i+1],op(i+2..nops(v),v)]:
if nops(S)=1 then
 v1:
else
S1:=[seq(S[i1]-1,i1=2..nops(S))]:
 Shrink1(v1,S1):
fi:
end:

#KickOut(pi,S): given a permutation pi and a set of indices S
#deletes the entries in the places of S and reduces
KickOut:=proc(pi,S) local i,S1,k:
k:=nops(pi):

if convert(S, set) minus {seq(i,i=1..k)}<>{} then
 ERROR(`Bad input`):
fi:

S1:={seq(i,i=1..nops(pi))} minus convert(S,set):
S1:=sort(convert(S1,list)):
redu([seq(pi[S1[i]],i=1..nops(S1))]):

end:










#DeleteEntriesV(v,S): given an increasing vector v and
#a set S of entries, deletes the entries indexed by S
#and reduces
DeleteEntriesV:=proc(v,S) local r,v1,i,S1:

if nops(S)=0 then
   RETURN(v):
fi:

S1:=sort(convert(S,list)):
r:=S1[nops(S1)]:

S1:=[op(1..nops(S1)-1,S1)]:

v1:=[op(1..r-1,v),seq(v[i]-1,i=r+1..nops(v))]:
DeleteEntriesV(v1,S1):
end:

DeleteEntriesSet:=proc(Setpi,S) local pi:
 {seq(DeleteEntries(pi,S),pi in Setpi)}:
end:



#PD1set(n,Patterns,pi,v,S): The "partial derivative" of
#WilfPre(n,Patterns,pi,v) w.r.t. to the entries
#of S, i.e.
#WilfPre(n-nops(S),Patterns,redu(pi w/o S), redu(v,w/o v[pi[S]]))
#minus the image of WilfPre(n,Patterns,pi,v) under the map
#delete the entries of S and reduce
#For example, try:
#PD1set(5,{[1,2,3]},[2,1],[4,5],{2});
PD1set:=proc(n,Patterns,pi,v,S) local gu,gu1,pi1,v1,i,S1:
gu:=WilfPre(n,Patterns,pi,v):

gu:=DeleteEntriesSet(gu,S):

 pi1:=DeleteEntries(pi,S):
S1:=[seq(pi[S[i]],i=1..nops(S))]:
S1:=sort(S1):
v1:=DeleteEntriesV(v,S1):
 gu1:=WilfPre(n-nops(S),Patterns,pi1,v1):

 gu1 minus gu:
end:




#IsRevDel11set(n0,Patterns,pi,S): are the  entries  of S in the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(n0,Patterns)?
#For example, try:
#IsRevDel11set(6,{[1,2,3]},[2,1],{1});
IsRevDel11set:=proc(n0,Patterns,pi,S) local gu,v:
gu:=IV1(n0,nops(pi)):
evalb({seq(evalb(PD1set(n0,Patterns,pi,v,S)={}),v in gu)}={true}):
end:


#IsRevDel1set(n0,Patterns,pi,S): are the entries of S in the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(i,Patterns), i<=n0?
#For example, try:
#IsRevDel1set(6,{[1,2,3]},[2,1],{1});
IsRevDel1set:=proc(n0,Patterns,pi,S) local i:
evalb({seq(IsRevDel11set(i,Patterns,pi,S),i=nops(pi)..n0)}={true}):
end:


#IsRevDel11GapSet(n0,Patterns,pi,S,G): Are the entries of S
#for the prefix permutation pi  reversely deletable
#for the Wilf class Wilf(n0,Patterns) satisying the gap condition
#given by the set of gap vectors G?
#For example, try:
#IsRevDel11GapSet(6,{[1,2,3]},[1,2],{1},{[0,1]});
IsRevDel11GapSet:=proc(n0,Patterns,pi,S,G) local gu,v,i:
if {seq(nops(G[i])-nops(pi),i=1..nops(G))}<>{1} then
 ERROR(`Bad input`):
fi:

gu:=IV1GapC(n0,nops(pi),G):

if gu={} then
 RETURN(true):
fi:
evalb({seq(evalb(PD1set(n0,Patterns,pi,v,S)={}),v in gu)}={true}):
end:

#IsRevDel1GapSet(n0,Patterns,pi,S,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class Wilf(i,Patterns), i<=n0?
#For example, try:
#IsRevDel1GapSet(6,{[1,2,3]},[2,1],{1},{[0,1]});
IsRevDel1GapSet:=proc(n0,Patterns,pi,S,G) local i:
evalb({seq(IsRevDel11GapSet(i,Patterns,pi,S,G),i=nops(pi)..n0)}={true}):
end:

hofkhi:=proc(pi)  local i, T:
for  i from 1 to nops(pi) do
T[pi[i]]:=i:
od:
[seq(T[i],i=1..nops(pi))]:
end:

Hofkhi:=proc(S)  local i:
 {seq(hofkhi(S[i]),i=1..nops(S))}:
end:



#IsIndeedScheme(n0,Patterns,S):checks whether the scheme S
#is OK for example, try:
#IsIndeedScheme(7,{[1,2,3,4]},Scheme(7,3,{[1,2,3,4]}));
IsIndeedScheme:=proc(n0,Patterns,S) local i,S1,pi,G,RevD,G1:

if not Bdok(S)  then
 print(S, `is not a genuine scheme`):
 RETURN(false):
fi:

for i from 1 to nops(S) do
 S1:=S[i]:
 pi:=S1[1]:
 G:=S1[2]:
 RevD:=S1[3]:
 G1:=GapVectors(n0,Patterns,pi):
 if G<>G1 then
   print(G , `is not the set of gap-vectors  for`, pi):
   RETURN(false):
 fi:

if G={} then
 if not IsRevDel1set(n0,Patterns,pi,RevD) then
   print(RevD, ` is not  reversely deletable for `, pi):
   RETURN(false):
 fi:
else
 if not IsRevDel1GapSet(n0,Patterns,pi,RevD,G) then
   print(RevD, ` is not  reversely deletable for `, pi):
   RETURN(false):
  fi:
fi:
od:

if SeqS(S,n0)<>[seq(nops(Wilf(i,Patterns)),i=0..n0)] then
# print(`The predicted sequence and the actual sequence do not agree`):
 RETURN(false):
fi:


true:
end:



#DeleteEntriesG(g,S): shrinks the vector g by removing the
#barriers in the places indicated by S
#For example, try: DeleteEntries([2,3,2],{2});
DeleteEntriesG:=proc(g,S) local r,S1,g1,i:

if S={} then
 RETURN(g):
fi:

r:=min(op(S)):

g1:=[op(1..r-1,g),g[r]+g[r+1],op(r+2..nops(g),g)]:
if nops(S)=1 then
RETURN(g1):
 else
 S1:=S minus {r}:
 S1:={seq(S1[i]-1,i=1..nops(S1))}:
 RETURN(DeleteEntriesG(g1,S1)):
fi:

end:


IsBadUntil:=proc(pi,P) local n, S, i, Ch:
n:=nops(pi):
if IsBad(pi,P) then S:={n}: else return {}: fi:
Ch:={pi}:
for i from n+1 do
Ch:={seq(op(Yeladim(Ch[j])),j=1..nops(Ch))}:
if {seq(IsBad(Ch[j],P),j=1..nops(Ch))}={true} then
S:=S union {i}:
else
return S:
fi:
od:
S:
end:



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

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

len:=add(v[i],i=1..nops(v))+nops(pi):

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

if nops(lu)=4 then
if member(len, lu[4]) then return 0: fi:
fi:

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

GapVectors:=lu[2]:



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

if v=[0$nops(v)] then
  RETURN(1):
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(Miklos(S,pi1,v1)):
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+Miklos(S,pi1,v1):
 od:
od:
su:
end:






#SeqS(S,K): the first K+1 terms computed by the
#enumeration scheme S, for example try:
#SeqS(Scheme(6,2,{[1,2,3]}),20);
SeqS:=proc(S,K) local i: [seq(Miklos(S,[],[i]) ,i=0..K)]:end:



#SchemeSlow(Patterns,Gvul): Finds a scheme for the Wilf class
#of restricted permutations avoiding the set of patterns
#Patterns, of depth<=Gvul. If it fails it retunrs
#FAIL followed by the partial scheme of depth Gvul
#followed by those prefix permutations that are not
#reducible (that have yet to be explored).
#For example, try:
#SchemeSlow({[1,2,3]},2);
SchemeSlow:=proc(Patterns,Gvul) local gu,Gvul1,godel,i:

if Patterns<>{} then
godel:=max(seq(nops(Patterns[i]),i=1..nops(Patterns))):
else
 godel:=3:
fi:

for Gvul1 from 0 to Gvul-1 do
gu:=Scheme1(godel+Gvul1,Gvul1,Patterns):

 if nops([gu])=1  then
    RETURN(gu):
 fi:
od:

gu:

end:


##############################
# functions for faster scheme
##############################

#ProveGapVector(g,pi,Patterns):: Given a gap vector g
#a prefix permutation pi, and a set of patterns Patterns
#proves that it is indeed a gap vector
#For example try: ProveGapVector([0,0,1],[1,2],{[1,2,3]});
ProveGapVector1:=proc(g,pi,Patterns) local lu,k,i,j:
k:=nops(pi):

if nops(g)<>k+1 then
 ERROR(`bad input`):
fi:

lu:={seq(seq(i-1+j/(g[i]+1),j=1..g[i]),i=1..nops(g))}:
lu:=permute(lu):
#print(lu):
lu:=[seq([op(pi),op(lu[i])],i=1..nops(lu))]:


evalb({seq(IsBad(lu[i],Patterns),i=1..nops(lu))}={true}):
end:

#ProveGapVector(g,pi,Patterns):: Given a gap vector g
#a prefix permutation pi, and a set of patterns Patterns
#proves that it is indeed a gap vector
#For example try: ProveGapVector([0,0,1],[1,2],{[1,2,3]});
ProveGapVector:=proc(g,pi,Patterns) local lu,k,i,j,a,b,C,co, Pat, Pat2:
Pat2:=Patterns:
for Pat in Patterns do

if member(splitbar(Pat)[4],{seq({$s..nops(Pat)},s=1..nops(Pat))}) then Pat2:=Pat2 minus {Pat}: fi:
od:
if Pat2<>Patterns then return ProveGapVector(g,pi,Pat2): fi:

k:=nops(pi):

if nops(g)<>k+1 then
 ERROR(`bad input`):
fi:

lu:={seq(seq(i-1+j/(g[i]+1),j=1..g[i]),i=1..nops(g))}:
lu:=permute(lu):
lu:=[seq([op(pi),op(lu[i])],i=1..nops(lu))]:

a:=nops(g):
co:=maxbars(Patterns):


C:=[seq(op(Comps(b,a)),b=1..co)];
#print(C):
for i from 1 to nops(C) do
C[i]:=Khaber(C[i],g);
od:

#print(lu,C):

evalb({seq(IsBad(lu[i],Patterns),i=1..nops(lu)),seq(ProveGapVector1(C[i],pi,Patterns),i=1..nops(C))}={true}):
end:


#maxbars(P) inputs the maximum number of bars 
#on any pattern in pattern set P (actually now returns the total number of bars
maxbars:=proc(P) local m, co, i, j, su:
m:=0: su:=0:
for j from 1 to nops(P) do
co:=0:
for i from 1 to nops(P[j]) do
if P[j][i][2]=1 then co:=co+1: fi:
od:
su:=su + co:
if co>m then m:=co: fi:
od:
#m:
su:
end:

#IsBad(perm1,Patterns): Inputs a permutation perm1 and
#a set of Patterns, outputs true if perm1 contains
#(at least one) pattern in Patterns and false otherwise
#For example try:
#IsBad([1,2,3,4],{[[1,0],[2,0],[3,0]]});
IsBad:=proc(perm1,Patterns) local pat1:
for pat1 in Patterns do
 if f(pat1,{perm1})={} then
     RETURN(true):
 fi:
od:
false:
end:

#Comps(a,n): the set of vectors of non-negative integers
#of length n that add up to a. For example, try:
#Comps(4,3);
Comps:=proc(a,n)
local gu,i,j,mu:
option remember:

if n=0 then
 if a=0 then
   RETURN({[]}):
 else
   RETURN({}):
 fi:
fi:

gu:={}:

for i from 0 to a do
 mu:=Comps(a-i,n-1):
 gu:=gu union {seq([op(mu[j]),i],j=1..nops(mu))}:
od:

gu:
end:

#CompsAd(a,n): the set of vectors of non-negative integers
#of length n that add up to <=a. For example, try:
#CompsAd(4,3);
CompsAd:=proc(a,n) local i:
{seq(op(Comps(i,n)),i=0..a)}:
end:

#Khaber(u,v): the component-wise sum of vectors u and v
Khaber:=proc(u,v) local i:[seq(u[i]+v[i],i=1..nops(u))]:end:

#ImpliedGaps1(g,s): given a gap vector g, and a positive
#integer s, finds the set of all vectors of non-negative
#integers of the same size as g and that add up to
#to s that are implied by g
#For example, try:
#ImpliedGaps1([1,1],3);
ImpliedGaps1:=proc(g,s)
local gu,n,g1:
n:=nops(g):
gu:=Comps(s-convert(g,`+`) ,n):
{seq(Khaber(g1,g), g1 in gu)}:
end:

#ImpliedGaps(G,s): given a set of gap vectors G, and a positive
#integer s, finds the set of all vectors of non-negative
#integers of the same size as g and that add up to
#to s that are implied by g in G
#For example, try:
#ImpliedGaps({[1,1],[0,2]},4);
ImpliedGaps:=proc(G,s) local g:
{seq(op(ImpliedGaps1(g,s)),g in G)}:
end:

#NewGapVectorsF(Patterns,pi,s,G): Given 
#a set of patterns Patterns and a prefix permutation
#pi  (of size k, say) and a positive integer s
#and a set of previously found gap vectors G
#finds all new gap vectors of size s , in the Fast way.
#For example, try: NewGapVectorsF({[[1,1],[3,0],[2,0]]},[1,2],1,{});
NewGapVectorsF:=proc(Patterns,pi,s,G) local gu,mu,g:

mu:=Comps(s,nops(pi)+1) minus ImpliedGaps(G,s):
gu:={}:
for g in mu do
 if ProveGapVector(g,pi,Patterns) then
   gu:=gu union {g}:
 fi:
od:
gu:

end:

#GapVectorsF(Patterns,pi,s): Given 
#a set of patterns Patterns and a prefix permutation
#pi  (of size k, say) and a positive integer s
#finds all gap vectors of size <=s , in the Fast way.
#For example, try: GapVectorsF({[[1,1],[3,0],[2,0]]},[1,2],1);
GapVectorsF:=proc(Patterns,pi,s) local G,s1:
G:={}:

for s1 from 0 to s do
 G:=G union NewGapVectorsF(Patterns,pi,s1,G):
od:
G:
end:

########################
#rd stuff begins here
########################
#UnfortunateEvents(pi,p,r): Given a prefix permutation pi,
#a pattern p, and a place r (1 between 1 and nops(pi)) 
#outputs all the subests of {1, ..., nops(pi)}
#that contain r such that the subperm of pi in these
#places reduces to the reduction of the appropriate
#prefix of p (of the same length)
#For example, try:
#UnfortunateEvents([2,1],[1,2,3],1);
UnfortunateEvents:=proc(pi,p,r) local s,S,i,k,s1,gu,T:
if r<1 or r>nops(pi) then
 ERROR(`Bad input`):
fi:
k:=nops(pi):
T:={seq(i,i=1..k)} minus {r}:
S:={seq(op(choose(T,i)),i=0..nops(p)-1)}:

S:={seq(s union {r}, s in S)}:

gu:={}:
for s in S do
 s1:=sort(convert(s,list)):
 if redu([seq(pi[s1[i]],i=1..nops(s1))])=
   redu([op(1..nops(s1),p)]) then
   gu:=gu union {s1}:
 fi:
od:


gu:
end:

#Scenarios1(pi,p,r,Gaps): Given a prefix permutation pi,
#a pattern p, and a place r (1 between 1 and nops(pi)) 
#and a set of gap-vectors
#outputs all the subests of {1, ..., nops(pi)}
#that contain r such that the subperm of pi in these
#places reduces to the reduction of the appropriate
#prefix of p (of the same length)
#For example, try:
#Scenarios1([2,1],[1,2,3],1,{});
Scenarios1:=proc(pi,p,r,Gaps) local s1,gu,mu,mu1,m:

gu:=UnfortunateEvents(pi,p,r):

mu:={}:

for s1 in gu do
mu1:=SubScenarios1(pi,p,s1,Gaps):
mu:= mu union {seq([s1,m],m in mu1)}:
od:

mu:

end:

#Scenarios(pi,P,r,Gaps): Given a prefix permutation pi,
#a set of patterns P, and a place r (1 between 1 and nops(pi)) 
#and a set of gap-vectors Gaps,
#outputs all pairs [p,S] where p
#is in P and S is a set of assignements to the
#pattern p relative to the elements of pi arising
#from secnarios that contain r such that the subperm of pi in these
#places reduces to the reduction of the appropriate
#prefix of p (of the same length)
#For example, try:
#Scenarios([2,1],{[1,2,3]},1,{});
Scenarios:=proc(pi,P,r,Gaps) local gu,mu,m,p:
gu:={}:

for p in P do
 mu:=Scenarios1(pi,p,r,Gaps):
 gu:=gu union {seq([p,m],m in mu)}:
od:


gu:
end:



#RestrictedComps(k,S,B): Given a positive integer k,
#an increasing sequence of integers S from [1, ..., k]
#and a vector B of non-negative integers of length
#nops(S)+1, returns all vectors v of length k+1
#with the property that if S=[i_1,i_2, ..., i_s]
#then v[1]+...+v[i_1]=B[1]
#v[i_1+1]+...+v[i2]=B[2]
#...
#v[i_s+1]+ ...+ v[k]=B[s+1]
#For example, try:
#RestrictedComps(4,[2],[3,3]);

RestrictedComps:=proc(K,S,B) local S1,B1,b,gu1,gu2,K1,i1,i2:
if nops(B)-nops(S)<>1 then
 ERROR(`Bad input`):
fi:

if S<>[] and S[nops(S)]>K then
 ERROR(`Bad input`):
fi:
 
if nops(S)=0 then
 RETURN(Comps(B[1],K)):
fi:

K1:=S[nops(S)]:
S1:=[op(1..nops(S)-1,S)]:
B1:=[op(1..nops(B)-1,B)]:
gu1:=RestrictedComps(K1,S1,B1):
b:=B[nops(B)]:
gu2:=Comps(b,K-K1):

{seq(seq([op(gu1[i1]),op(gu2[i2])],i2=1..nops(gu2)),i1=1..nops(gu1))}:
end:



#IsBigger1(v,g): is vector v point-wise bigger than
#vector g. For example, try: IsBigger1([1,2],[2,4]);
IsBigger1:=proc(v,g) local i:
if nops(v)<>nops(g) then
 ERROR(`Bad input`):
fi:
evalb({seq(evalb(v[i]-g[i]>=0),i=1..nops(v))}={true});
end:


#IsBigger(v,G): is v bigger than at least one of the vectors in G?
IsBigger:=proc(v,G) local g:
evalb(G<>{} and {seq(IsBigger1(v,g),g in G)}<>{false}):
end:

#RestrictedCompsGaps(k,S,B,Gaps): Given a positive integer k,
#an increasing sequence of integers S from [1, ..., k]
#and a vector B of non-negative integers of length
#nops(S)+1, and a set of vectors Gaps
#returns all vectors v of length k+1
#not implied by any the vectors in Gaps,
#with the property that if S=[i_1,i_2, ..., i_s]
#then v[1]+...+v[i_1]=B[1]
#v[i_1+1]+...+v[i2]=B[2]
#...
#v[i_s+1]+ ...+ v[k]=B[s+1]
#For example, try:
#RestrictedComps(4,[2],[3,3],{});
RestrictedCompsGaps:=proc(K,S,B,Gaps) local mu,gu,v:
gu:=RestrictedComps(K,S,B):
mu:={}:

for v in gu do
 if not IsBigger(v,Gaps) then
   mu:=mu union {v}:
 fi:
od:
mu:
end:

#SubScenarios1(pi,p,s1,Gaps): Given a prefix permutation pi,
#a pattern p,  and a scenario s1,
#such that the subperm of pi in these
#places reduces to the reduction of the appropriate
#prefix of p (of the same length)
#and a set of gap vectors Gaps
#outputs the possible ways ALL the members of p
#can be placed relative to the members of the pi
#subject to the gap-conditions
#For example, try:
#SubScenarios1([2,1],[1,2,3],1,[1],{});
SubScenarios1:=proc(pi,p,s1,Gaps) 
local gu,pi1,p1,i,L1,T,gu1,lu,co,a,Gu,j:

pi1:=[seq(pi[s1[i]],i=1..nops(s1))]:
p1:=[op(1..nops(s1),p)]:

if redu(pi1)<>redu(p1) then
  print(`Bad input`):
  RETURN(FAIL):
fi:

pi1:=sort(pi1):
p1:=sort(p1):

L1:=[p1[1]-1,seq(p1[i]-p1[i-1]-1,i=2..nops(p1)),nops(p)-p1[nops(p1)]]:

lu:=[]:
for i from 1 to nops(pi) do
 if member(i,convert(pi1,set)) then
  lu:=[op(lu),1]:
  else
  lu:=[op(lu),0]:
 fi:
od:


gu:=RestrictedCompsGaps(nops(pi)+1,pi1,L1,Gaps):


Gu:={}:
for a from 1 to nops(gu) do
gu1:=gu[a]:

co:=0:

for i from 1 to nops(gu1) do
  for j from 1 to gu1[i] do

    co:=co+1:
    T[co]:=i-1+j/(gu1[i]+1):
  od:
  if i<=nops(lu) and lu[i]=1 then
     co:=co+1:
     T[co]:=i:
  fi:
od:

Gu:=Gu union {[seq(T[j],j=1..nops(p))]}:

od:

Gu:
  
end:

#SetRevDelGapF(Patterns,pi,G): 
#Using the Fast way, the set of reversely deleteable
#entries for the Wilf class avooding Patterns for the
#prefix permutation pi and set of gap vectors G. For example, try:
#SetRevDelGapF({[1,2,3]},[1,2],{[0,0,1]);
SetRevDelGapF:=proc(Patterns,pi,G) local gu,r,S,S2:
gu:={}:
for r from 1 to nops(pi) do
 if {seq(IsRevDelGapF(Patterns[p],Patterns,pi,r,G),p=1..nops(Patterns))}={true} then
   gu:=gu union {r}:
 fi:
od:
gu:
end:


#IsRevDel1GapF(Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class of forbidden patterns Patterns
#For example, try:
#IsRevDel1GapF({[1,2,3]},[2,1],1,{[0,1]});
IsRevDelGapF:=proc(P,Pats,pi,r,G)
if IsRevDel1GapF(P,Pats,pi,r,G) and IsRevDel2GapF(P,Pats,pi,r,G) then return true:
else return false:
fi:
end:

#############################
#extra functions to deal w/
#r.d. elts with barred pats
#############################

#BarredScenarios(Pat,pi,r,G): inputs a barred pattern Pat, a prefix pi, 
#a position r of pi and a set of gaps G
#returns all copies of the forbidden unbarred pattern which are NOT
#part of a copy of the barred pattern
#for example, try: BarredScenarios([[1,1],[3,0],[2,0]],[1],1,{});
BarredScenarios:=proc(Pat,pi,r,G) local p, Sl, Ss, i, Saved,temp,w1,w2,j:
p:=splitbar(Pat):

Sl:=Scenarios(pi,{p[2]},r,G):
Ss:=Scenarios(pi,{p[1]},r,G):

#if evalb(pi=p[1]) then Ss:=Ss union {[pi,[[seq(i,i=1..nops(pi))],sort(pi)]  ]}: fi:
Saved:=Ss:

#print(Sl, Ss):

if evalb(Sl=Ss) then return Saved: fi:


for i from 1 to nops(Sl) do
	w1:=[]:
	for j from 1 to nops(Sl[i][2][1]) do
	if j in p[3] then
	w1:=[op(w1),Sl[i][2][1][j]]:
	fi:
	od:
	
	w2:={op(Sl[i][2][2])}: #what to delete?
	for j from 1 to nops(pi) do
	if j in p[4] then
	w2:=w2 minus {pi[j]}:
	fi:
	od:
	w2:=sort([op(w2)]):

temp:=[p[1],[w1,w2]]:
#print(Sl[i], " becomes ", temp):

Saved:=Saved minus {temp}:
od:

Saved:
end:

#IsRevDel1GapF(Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class of forbidden patterns Patterns
#For example, try:
#IsRevDel1GapF({[1,2,3]},[2,1],1,{[0,1]});
IsRevDel1GapF:=proc(P,Pats,pi,r,Gaps) local gu,i,p,gu1,mu1,j1,
Overlap,beyakhad:
gu:=BarredScenarios(P,pi,r,Gaps):
#print(gu):

for i from 1 to nops(gu) do
 gu1:=gu[i]:
 p:=gu1[1]: Overlap:=nops(gu1[2][1]): mu1:=gu1[2][2]:

   beyakhad:=[op(1..r-1,pi),op(r+1..nops(pi),pi),
      seq(mu1[p[j1]],j1=Overlap+1..nops(p))]:
#print(beyakhad, Pats, IsBad(beyakhad,Pats)):
     if not IsBad(beyakhad,Pats) then
       RETURN(false):
     fi:

 od:
if not IsRevDel3GapF(P,Pats,pi,r,Gaps) then return false: fi:

true:
end:

BarredScenarios2:=proc(P,pi,r,Gaps) local p,Sl, i, j, S2,pos:
p:=splitbar(P):
#print(p,pi,{p[2]},r,Gaps):
Sl:=Scenarios(pi,{p[2]},r,Gaps):
S2:=Sl:
#print(S2):

if evalb(p[1]=p[2]) then return {}: fi:

for i from 1 to nops(Sl) do
pos:={}:
	for j from 1 to nops(Sl[i][2][1]) do
	if member(j,p[4]) then
        pos:=pos union {Sl[i][2][1][j]}:
	fi:
	od:
#print(Sl[i],pos):
if not member(r,pos) then
S2:=S2 minus {Sl[i]}:
fi:
od:
S2:
end:

#IsRevDel2GapF(Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely insertable
#for the Wilf class of forbidden patterns Patterns
#For example, try:
#IsRevDel2GapF({[1,2,3]},[2,1],1,{[0,1]});
IsRevDel2GapF:=proc(P,Pats,pi,r,Gaps) local gu,i,p,gu1,mu1,j1,
Overlap,beyakhad:
gu:=BarredScenarios2(P,pi,r,Gaps):
#print(gu):

for i from 1 to nops(gu) do
 gu1:=gu[i]:
 p:=gu1[1]: Overlap:=nops(gu1[2][1]): mu1:=gu1[2][2]:

   beyakhad:=[op(1..r-1,pi),op(r+1..nops(pi),pi),
      seq(mu1[p[j1]],j1=Overlap+1..nops(p))]:
#print(beyakhad, Pats, IsBad(beyakhad,Pats)):
     if IsBad(beyakhad,Pats) then
       RETURN(false):
     fi:

 od:

true:
end:



#inputs a barred permutation p and outputs
#a pair of permutations [p1,pb], where 
#p1 is the unbarred part of p, and pb is the entire
#permutation p, sans bars, for example try
#splitbar([[1,1],[3,0],[2,0]]);
splitbar:=proc(p) local i, p1, pb, pos:
p1:=[]: pb:=[]: pos:={}:
for i from 1 to nops(p) do
if p[i][2]=0 then
p1:=[op(p1),p[i][1]]: pos:=pos union {i}:
fi:
pb:=[op(pb),p[i][1]]:
od:

[redu(p1),pb, pos, {seq(i,i=1..nops(p))} minus pos]:
end:


BarredScenarios3:=proc(P,Pats,pi,r,G) local gu, i, gu1, p, Overlap, mu1, beyakhad, L, L2, j, gu2, k,z,perm, max, temp:
gu:=BarredScenarios(P,pi,r,G):
gu2:={}:
max:=nops(splitbar(Pats[1])[4]):
for i from 2 to nops(Pats) do
temp:=nops(splitbar(Pats[i])[4]):
if temp>max then max:=temp: fi:
od:

#print(gu):

perm:=splitbar(P):

for z from 1 to max do
for i from 1 to nops(gu) do
 gu1:=gu[i]:

 p:=gu1[1]: Overlap:=nops(gu1[2][1]): mu1:=gu1[2][2]:

   beyakhad:=[op(pi), seq(mu1[p[j1]],j1=Overlap+1..nops(p))]:
L:=sort(beyakhad):
	L2:=[L[1]/2]:
	for j from 1 to nops(L)-1 do
		L2:=[op(L2),(L[j]+L[j+1])/2]:
	od:
	L2:=[op(L2),L[nops(L)]+1/2]:

for k from 1 to nops(L2) do
 gu2:=gu2 union {[gu1[1],[gu1[2][1],[op(gu1[2][2])],[L2[k]]]]}:
od:

od:
gu:=gu2:
od:
gu:
end:

IsRevDel3GapF:=proc(P,Pats,pi,r,Gaps) local gu,i,p,gu1,mu1,j1,
Overlap,beyakhad, num,j,Inter,Seq2,Seq1,S,k,l,temp,q,m,q2,beyakhadlong:
gu:=BarredScenarios3(P,Pats,pi,r,Gaps):
#print("check 3GapF ",gu):

if nops(gu)>0 and nops(gu[1][2])=3 then
num:=nops(gu[1][2][2])+nops(gu[1][2][3]):
fi:

if nops(gu)>0 and nops(gu[1][2])<3 then return true: fi:

#print(num):

for i from 1 to nops(gu) do
 gu1:=gu[i]:
 p:=gu1[1]: Overlap:=nops(gu1[2][1]): mu1:=gu1[2][2]:
#print(num,Overlap,gu1,choose([seq(j,j=Overlap+1..num)],nops(gu1[2][3]))):
 Inter:=choose([seq(j,j=Overlap+1..num)],nops(gu1[2][3])):

 Seq2:=permute(gu1[2][3]):
 Seq1:=[seq(mu1[p[j1]],j1=Overlap+1..nops(p)),seq(gu1[2][2][j2],j2=nops(p)+1..nops(gu1[2][2]))]:

#print(Seq1,Seq2):

 S:={}:

 for k from 1 to nops(Seq2) do
  for l from 1 to nops(Inter) do
   temp:=[]: p:=1: q:=1:

   for m from Overlap+1 to num do
    if member(m, {op(Inter[l])}) then
    temp:=[op(temp),Seq2[k][p]]: p:=p+1:
    else temp:=[op(temp),Seq1[q]]: q:=q+1:
    fi:
   od:


   S:=S union {temp}:
  od:
 od:

#print(S):
for q2 from 1 to nops(S) do

   beyakhad:=[op(1..r-1,pi),op(r+1..nops(pi),pi), op(S[q2])]:
 beyakhadlong:=[op(1..nops(pi),pi), op(S[q2])]:
#print(beyakhad, Pats, IsBad(beyakhad,Pats)):
     if IsBad(beyakhadlong,Pats) and not IsBad(beyakhad,Pats) then
       RETURN(false):
     fi:
od:


od:

true:
end:



##############################
#schemefast proc's
##############################

#GaBchildrenF(Patterns,GvulGap,pi): the good and bad children
#of a prefix permutation pi w.r.t. to the set of
#patterns Patterns, with gap vectors of sum<=GvulGap
#For example, try:
#GaBchildrenF({[1,2,3]},1,[1]):
GaBchildrenF:=proc(Patterns,GvulGap,pi) local gu,G,B,pi1,Gaps1,rd1:
gu:=Yeladim(pi):
G:={}: B:={}:

for pi1 in gu do

 Gaps1:=GapVectorsF(Patterns,pi1,GvulGap):
    rd1:=SetRevDelGapF(Patterns,pi1,Gaps1):
    if rd1<>{} then
       G:=G union {[pi1,Gaps1,rd1]}:
    else
       B:=B union {[pi1,Gaps1,rd1]}:
    fi:
od:

G,B:

end:

#Scheme1F(Gvul,GvulGap,Patterns): Tries to find a scheme
#of depth<=Gvul and with gap-vectors of sum<=GvulGap
#for the Wilf class avoiding
#the patterns in Patterns using empirical investigations
#if it fails it returns FAIL, followed by the partial
#scheme and the leftovers. For example, try:
#Scheme1F(2,1,{[1,2,3]});
Scheme1F:=proc(Gvul,GvulGap,Patterns) local LeftToDo,S,pi,j,i,gu:
LeftToDo:={[]}:

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

for i from 0 to Gvul while LeftToDo<>{} do 

  for pi in LeftToDo do
     gu:=GaBchildrenF(Patterns,GvulGap,pi) :
     S:=[op(S), op(gu[1]),op(gu[2])]:
     LeftToDo:=LeftToDo minus {pi} union 
       {seq(gu[2][j][1],j=1..nops(gu[2]))}:
  od:

if LeftToDo={} then
    S:=SimplifyScheme(S):
 if Bdok(S) then
   RETURN(S):
 else
 print(`Not closed under deletion`):
 RETURN(FAIL,S):
 fi:
fi:
od:

FAIL,S,LeftToDo:
end:






#SchemeFast(Patterns,Gvul,GvulGap): Finds a scheme for the Wilf class
#of restricted permutations avoiding the set of patterns
#Patterns, of depth<=Gvul and with gap vectors
#of sum <=GvulGap . If it fails it retunrs
#FAIL followed by the partial scheme of depth Gvul
#followed by those prefix permutations that are not
#reducible (that have yet to be explored).
#For example, try:
#SchemeFast({[1,2,3]},2,1);
SchemeFast:=proc(Patterns,Gvul,GvulGap) local gu,Gvul1,i, S3:

for Gvul1 from 0 to Gvul-1 do
gu:=Scheme1F(Gvul1,GvulGap,Patterns):

 if nops([gu])=1  then

S3:={}:
for i from 1 to nops(gu) do
if evalb(gu[i][2]={[0$(nops(gu[i][1])+1)]}) then
S3:=S3 union {[gu[i][1],gu[i][2],gu[i][3],{}]}:
elif evalb(gu[i][1]=[]) then 
S3:=S3 union {[gu[i][1],gu[i][2],gu[i][3],{}]}:
else
S3:=S3 union {[gu[i][1],gu[i][2],gu[i][3],IsBadUntil(gu[i][1],Patterns)]}:
fi:
od:

    RETURN(S3):
 fi:
od:

gu:

end:

##############################
# functions for Sipur
##############################

rev:=proc(L) local L2, i:
L2:=[]:
for i from 1 to nops(L) do
L2:=[L[i],op(L2)]:
od:
L2:
end:

Rev:=proc(S) local S2, i:
S2:={}:
for i from 1 to nops(S) do
S2:=S2 union {rev(S[i])}:
od:
end:

comp:=proc(L) local n, w,i:
n:=nops(L):
w:=[]:
for i from 1 to n do
w:=[op(w),[n+1-L[i][1],L[i][2]]]:
od:
w:
end:

Comp:=proc(S) local S2, i:
S2:={}:
for i from 1 to nops(S) do
S2:=S2 union {comp(S[i])}:
od:
end:

inv:=proc(pi)  local i, T,S:
for  i from 1 to nops(pi) do
T[pi[i][1]]:=i:
 
S[pi[i][1]]:=pi[i][2]: 
od:
[seq([T[i],S[i]],i=1..nops(pi))]:
end:

Inv:=proc(S) local S2, i:
S2:={}:
for i from 1 to nops(S) do
S2:=S2 union {inv(S[i])}:
od:
end:

#build ALL barred patterns of length n
Pats:=proc(n) local i, j, k, P, Bt, B, S:
P:=permute(n):
Bt:=choose([0$n,1$n-1],n):
B:={seq(op(permute(Bt[i])),i=1..nops(Bt))}:
S:={}:

for i from 1 to nops(P) do
for j from 1 to nops(B) do
S:=S union {[seq([P[i][k],B[j][k]],k=1..n)]}:
od:
od:
S:
end:

#build ALL barred patterns of length n
PatsB:=proc(n,a) local i, j, k, P, Bt, B, S:
P:=permute(n):
Bt:=[0$n-a,1$a]:
B:={seq(op(permute(Bt)),i=1..nops(Bt))}:
S:={}:

for i from 1 to nops(P) do
for j from 1 to nops(B) do
S:=S union {[seq([P[i][k],B[j][k]],k=1..n)]}:
od:
od:
S:
end:

KHAVERIM:=proc(Perms):
{Perms} union {Rev(Perms)} union {Inv(Perms)} union
 {Comp(Perms)} union
{Rev(Inv(Perms))} union {Comp(Rev(Perms))} union {Comp(Inv(Perms))} union {Comp(Inv(Comp(Perms)))}:
end:

#AllPatternSets1(L): inputs a list of pairs and outputs all pattern
#sets obeying that list for example
#AllPatternSets1([[3,1],[3,1]]); returns all pairs of barred patterns of length 3
#with one bar on each pattern
AllPatternSets1:=proc(L) local S, i, Stemp,j,k,S2:
S:=PatsB(op(L[1])):
S2:={}:
for i from 1 to nops(S) do
S2:=S2 union {{S[i]}}:
od:
S:=S2:


#print("first perms are ", S):
for i from 2 to nops(L) do
Stemp:=PatsB(op(L[i])):
#print("S is ", S, " S2 is ",S2):
#print("Stemp is ", Stemp):
S2:={}:
for j from 1 to nops(S) do
for k from 1 to nops(Stemp) do
#print("union ",S[j], {Stemp[k]}):

S2:=S2 union {S[j] union {Stemp[k]}}:
od:
od:
#print("S is ", S, " S2 is ",S2):

S:=S2:
od:

S2:=S:
for i from 1 to nops(S2) do
if nops(S2[i])<>nops(L) then S:=S minus {S2[i]}: fi:
od:

S:
end:

#takes AllPatternSets1(L) and breaks them down into Wilf-equivalent groups
AllPatternSets:=proc(L) local gu,Gu,S,K,i:
S:=AllPatternSets1(L):
Gu:={}:

while S<>{} do
 gu:=S[1]:
K:=KHAVERIM(gu):
 Gu:=Gu union {K}:
 S:=S minus K:
od:
Gu:
end:

#SchemeImage(Patterns,Gvul): tries to find a scheme
#for an equivalent set of patterns
SchemeImage:=proc(Patterns,GVUL,GvulGap) local mu,S,i:
mu:=KHAVERIM(Patterns):

for i from 1 to nops(mu) do
 S:=SchemeFast(mu[i],GVUL,GvulGap):
  if S[1]<>FAIL then
   RETURN([mu[i],S]):
  fi:
od:

FAIL:

end:

#Sipur(L,Gvul,K): Everything you ever wanted to know
#about all patterns of type L, that schemes of depth <=Gvul
#and printing out, in the successful cases, the first K+1
#terms. For example, try:
#Sipur([3],2,20);
Sipur:=proc(L,Gvul,GvulGap,K)
local gu,S,F,sch,i,ka,mu, Se, Se2:
gu:=AllPatternSets(L);
print(`There all together`, nops(gu), `different equivalence classes `):
S:={}:
F:={}:

for i from 1 to nops(gu) do
sch:=SchemeImage(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[i][1]),i=1..nops(mu))):
   print(`the member `, sch[1], `has a scheme of depth `, ka):
   print(`here it is: `):
   print(mu):
if K>7 then
Se:=[seq(nops(Wilf(j,gu[i][1])),j=0..6)]:
Se2:=[op(1..7,SeqS(mu,K))]:
else Se:=[]: Se2:=[]:
fi:
 print(`Naively, we would expect the sequence to begin `, seq(nops(Wilf(j,gu[i][1])),j=0..6));
   print(`Using  the scheme, the first, `, K+1, `terms are `):
   print(SeqS(mu,K)):
if not evalb(Se=Se2) then print("CHECK FOR A BAD SCHEME!!!!!!!!!!!!"): fi:
   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):
F:
end: