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("This version last updated March 21, 2008."):
print(""):
print("The main functions are : Sipur(L,Depth,K)"):
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 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]=Yeladim then
print("Yeladim(pi): all the children of pi"):

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]=IsBadUntil then
print("IsBadUntil(pi,P): the set of stop points for prefix pi and patterns P"):

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]=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]=ProveGapVector1 then
print("ProveGapVector1(g,pi,Patterns):: Given a gap vector g"):
print("a prefix permutation pi, and a set of patterns Patterns"):
print("proves that there are no permutations beginning with pi that avoid"):
print("Patterns and obey g exactly"):
print("For example try: ProveGapVector1([0,0,1],[1,2],{[[1,0],[2,0],[3,0]]});"):

elif args[1]=ProveGapVector then
print("ProveGapVector(g,pi,Patterns):: Given a potential gap vector g"):
print("a prefix permutation pi, and a set of patterns Patterns proves that it is "):
print("indeed a gap vector by checking that its children are also gap vectors"):
print("For example try: ProveGapVector([0,0,1],[1,2],{[[1,0],[2,0],[3,0]]});"):

elif args[1]=maxbars then
print("maxbars(P) inputs the total number of bars on all patterns in pattern set P"):

elif args[1]=IsBad then
print("IsBad(perm1,Patterns): Inputs a permutation perm1 and"):
print("a set of Patterns, outputs true if perm1 contains"):
print("(at least one) pattern in Patterns and false otherwise"):
print("For example try: IsBad([1,2,3,4],{[[1,0],[2,0],[3,0]]});"):

elif args[1]=Comps then
print("Comps(a,n): the set of vectors of non-negative integers"):
print("of length n that add up to a. For example, try: Comps(4,3);"):

elif args[1]=CompsAd then
print("CompsAd(a,n): the set of vectors of non-negative integers"):
print("of length n that add up to <=a. For example, try: CompsAd(4,3);"):

elif args[1]=Khaber then
print("Khaber(u,v): the component-wise sum of vectors u and v"):

elif args[1]=ImpliedGaps1 then
print("ImpliedGaps1(g,s): given a gap vector g, and a positive integer s, finds the "):
print("set of all vectors of non-negative integers of the same size as g and that add up to"):
print("s that are implied by g For example, try: ImpliedGaps1([1,1],3);"):

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

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

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

elif args[1]=UnfortunateEvents then
print("UnfortunateEvents(pi,p,r): Given a prefix permutation pi, a pattern p, "):
print("and a place r (1 between 1 and nops(pi)) outputs all the subests of {1, ..., nops(pi)}"):
print("that contain r such that the subperm of pi in these places reduces to the reduction "):
print("of the appropriate prefix of p (of the same length) "):
print("For example, try: UnfortunateEvents([2,1],[1,2,3],1);"):

elif args[1]=Scenarios1 then
print("Scenarios1(pi,p,r,Gaps): Given a prefix permutation pi, a pattern p, "):
print("and a place r (1 between 1 and nops(pi)) and a set of gap-vectors"):
print("outputs all the subsets of {1, ..., nops(pi)} that contain r such that "):
print("the subperm of pi in these places reduces to the reduction of the appropriate"):
print("prefix of p (of the same length) For example, try:"):
print("Scenarios1([2,1],[1,2,3],1,{});"):

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

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

elif args[1]=IsBigger1 then
print("IsBigger1(v,g): is vector v point-wise bigger than"):
print("vector g. For example, try: IsBigger1([1,2],[2,4]);"):

elif args[1]=IsBigger then
print("IsBigger(v,G): is v bigger than at least one of the vectors in G?"):

elif args[1]=RestrictedCompsGaps then
print("RestrictedCompsGaps(k,S,B,Gaps): Given a positive integer k,"):
print("an increasing sequence of integers S from [1, ..., k] and a vector B of non-negative "):
print("integers of length nops(S)+1, and a set of vectors Gaps returns all vectors v of length k+1"):
print("not implied by any the vectors in Gaps, with the property that if S=[i_1,i_2, ..., i_s]"):
print("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]."):
print("For example, try: RestrictedComps(4,[2],[3,3],{});"):

elif args[1]=SubScenarios1 then
print("SubScenarios1(pi,p,s1,Gaps): Given a prefix permutation pi, a pattern p,  and a scenario s1,"):
print("such that the subperm of pi in these places reduces to the reduction of the appropriate"):
print("prefix of p (of the same length) and a set of gap vectors Gaps"):
print("outputs the possible ways ALL the members of p can be placed relative to the "):
print("members of the pi subject to the gap-conditions"):

elif args[1]=SetRevDelGapF then
print("SetRevDelGapF(Patterns,pi,G): Using the Fast way, the set of reversely deleteable"):
print("entries for the Wilf class avooding Patterns for the prefix permutation pi and set of gap vectors G."):

elif args[1]=IsRevDel1GapF then
print("IsRevDel1GapF(Patterns,pi,r,G): Is the r^th entry of the"):
print("prefix permutation pi  reversely deletable"):
print("for the Wilf class of forbidden patterns Patterns"):

elif args[1]=BarredScenarios then
print("BarredScenarios(Pat,pi,r,G): inputs a barred pattern Pat, a prefix pi, "):
print("a position r of pi and a set of gaps G"):
print("returns all copies of the forbidden unbarred pattern which are NOT"):
print("part of a copy of the barred pattern"):
print("for example, try: BarredScenarios([[1,1],[3,0],[2,0]],[1],1,{});"):

elif args[1]=IsRevDel1GapF then
print("IsRevDel1GapF(P, Patterns,pi,r,G): checks if insertion of pi_r is a bijection"):
print("of permutations avoiding P (in Patterns) and obeying the gap conditions in Gaps"):
print("for scenarios using the minimum number of letters needed to create a forbidden pattern"):

elif args[1]=BarredScenarios2 then
print("BarredScenarios2(P,pi,r,Gaps): makes a list of all possible ways for pi to be involved"):
print("in a forbidden pattern P without using pi_r but while obeying Gaps"):

elif args[1]=IsRevDel2GapF then
print("IsRevDel2GapF(P, Patterns,pi,r,G): checks if deletion of pi_r is a bijection"):
print("of permutations avoiding P (in Patterns) and obeying the gap conditions in G"):

elif args[1]=splitbar then
print("splitbar(p): inputs a barred permutation p and outputs"):
print("a pair of permutations [p1,pb], followed by the positions occupied by p1"):
print("and the set of positions not occupied by p1 where "):
print("p1 is the unbarred part of p, and pb is the entire"):
print("permutation p, sans bars, for example try"):
print("splitbar([[1,1],[3,0],[2,0]]);"):

elif args[1]=BarredScenarios3 then
print("BarredScenarios3(P,Pats,pi,r,G): makes a list of all possible ways for pi_r to be"):
print("involved in a forbidden pattern P (in Pats), while still obeying gaps G"):

elif args[1]=IsRevDel3GapF then
print("IsRevDel3GapF(P,Pats,pi,r,Gaps): checks if insertion of pi_r is a bijection"):
print("of permutations avoiding P (in Pats) and obeying the gap conditions in Gaps"):
print("for scenarios with extra letters"):

elif args[1]=Scheme1F then
print("Scheme1F(Depth,GapDepth,Patterns): Tries to find a scheme"):
print("of depth<=Depth and with gap-vectors of sum<=GapDepth"):
print("without stop points"):
print("for the Wilf class avoiding the patterns in Patterns"):
print("if it fails it returns FAIL, followed by the partial"):
print("scheme and the leftovers. For example, try:"):
print("Scheme1F(2,1,{[[1,1],[3,0],[2,0]]});"):

elif args[1]=GaBchildrenF then
print("GaBchildrenF(Patterns,GapDepth,pi): the good and bad children"):
print("of a prefix permutation pi w.r.t. to the set of"):
print("patterns Patterns, with gap vectors of sum<=GapDepth"):
print("For example, try: GaBchildrenF({[[1,1],[3,0],[2,0]]},1,[1]):"):

elif args[1]=SchemeImage then
print("SchemeImage(Patterns,Depth, GapDepth): tries to find a scheme"):
print("for an equivalent set of patterns"):

elif args[1]=SchemeFast then
print("SchemeFast(Patterns,Depth,GapDepth): Finds a scheme for the Wilf class"):
print("of restricted permutations avoiding the set of patterns including stop points"):
print("Patterns, of depth<=Depth and with gap vectors"):
print("of sum <=GapDepth . If it fails it returns"):
print("FAIL followed by the partial scheme of depth Depth"):
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);"):

elif args[1]=KHAVERIM then
print("KHAVERIM(Perms): the set of patterns that are"):
print("Wilf-equivalent to Perms"):
print("for example, try: KHAVERIM({[[1,1],[3,0],[2,0]],[[1,1],[3,0],[2,0],[4,0]]});"):

elif args[1]=PatsB then
print("PatsB(n,a): builds ALL barred patterns of length n with a bars"):

elif args[1]=Pats then
print("Pats(n): builds ALL barred patterns of length n with <n bars"):

elif args[1]=Inv then
print("Inv(S): finds the inverse of each pattern in set S"):

elif args[1]=inv then
print("inv(pi): finds the inverse of pattern pi"):

elif args[1]=Comp then
print("Comp(S): finds the complement of each pattern in set S"):

elif args[1]=comp then
print("comp(pi): finds the complement of pattern pi"):

elif args[1]=Rev then
print("Rev(S): finds the reverse of each pattern in set S"):

elif args[1]=rev then
print("rev(pi): finds the reverse of pattern pi"):

elif args[1]=AllPatternSets1 then
print("AllPatternSets1(L): inputs a list of pairs and outputs all pattern"):
print("sets obeying that list for example"):
print("AllPatternSets1([[3,1],[3,1]]); returns all pairs of barred patterns of length 3"):
print("with one bar on each pattern"):

elif args[1]=AllPatternSets then
print("AllPatternSets(L): takes AllPatternSets1(L) and"):
print("breaks them down into Wilf-equivalent groups"):

elif args[1]=Sipur then
print("Sipur(L,Depth, GapDepth, K): Everything you ever wanted to know"):
print("about all patterns of type L, that schemes of depth <=Depth"):
print("with GapVectors of maximum weight GapDepth"):
print("and printing out, in the successful cases, the first K+1"):
print("terms. For example, try:"):
print("Sipur([[3,1]],2,1,20);"):

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:


#Yeladim(pi): all the children of pi
Yeladim:=proc(pi) local i:[seq(Append1(pi,i),i=1..nops(pi)+1)]: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:




#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:


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:


##############################
# 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):
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)];
for i from 1 to nops(C) do
C[i]:=Khaber(C[i],g);
od:

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
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.
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
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):

Saved:=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]]:

Saved:=Saved minus {temp}:
od:

Saved:
end:

#IsRevDel1GapF(P, Patterns,pi,r,G): Is the r^th entry of the
#prefix permutation pi  reversely deletable
#for the Wilf class of forbidden patterns Patterns
IsRevDel1GapF:=proc(P,Pats,pi,r,Gaps) local gu,i,p,gu1,mu1,j1,
Overlap,beyakhad:
gu:=BarredScenarios(P,pi,r,Gaps):
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))]:
     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(P,pi,r,Gaps)
BarredScenarios2:=proc(P,pi,r,Gaps) local p,Sl, i, j, S2,pos:
p:=splitbar(P):
Sl:=Scenarios(pi,{p[2]},r,Gaps):
S2:=Sl:

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:
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):

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))]:
     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(P,Pats,pi,r,G)
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:

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(P,Pats,pi,r,Gaps)
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):
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:

for i from 1 to nops(gu) do
 gu1:=gu[i]:
 p:=gu1[1]: Overlap:=nops(gu1[2][1]): mu1:=gu1[2][2]:
 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]))]:
 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:

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])]:
     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,1],[3,0],[2,0]]},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,1],[3,0],[2,0]]});
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,1],[3,0],[2,0]]},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:

if evalb(SeqS(S3,7)<>[seq(nops(Wilf(i,Patterns)),i=0..7)]) then return FAIL:
fi:
    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,Depth, GapDepth): 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,Depth, GapDepth, K): Everything you ever wanted to know
#about all patterns of type L, that schemes of depth <=Gvul
#and maximum gap vector weight GapDepth
#and printing out, in the successful cases, the first K+1
#terms. For example, try:
#Sipur([[3,1]],2,1,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: