with(combinat):

Help:=proc()
if args = NULL then
	print(`HasScheme(pats,max_depth,max_gap_size,quick_fail), SeqS(scheme, n)`):
elif args = 'HasScheme' then
	print(`HasScheme(pats,max_depth, max_gap_size,quick_fail:=false) attempts to find a scheme to count the permutations avoiding all the patterns of pats with depth <= max_depth and gap restrictions of size <=max_gap_size. If it succeeds, it returns the scheme, and if it fails it returns FAIL along with the set of prefixes it was unable to reduce. By default, quick_fail is set to false. If it is set to true instead, then HasScheme does a depth-first search of all prefixes of length <= max_depth, and stops as soon as it finds one of length max_depth that doesn't reduce. The advantage is that it will return FAIL more quickly when no scheme exists, but the disadvantage is that it will only list one prefix that it is unable to reduce instead of all of them.`):
	print(`For example, try`):
	print(`HasScheme({[1, 3, 4, 2], [3, 1, 2, 4]}, 4, 2);`):
elif args = 'SeqS' then
	print(`SeqS(scheme, n) gives the first n terms of the enumeration sequence defined by scheme`):
	print(`For example, try`):
	print(`SeqS({[[1, 2], [[[0, 0, 1], 1, []]], [2, [1]]], [[2, 1], [[[0, 0, 0], 2, [1]]], [1, [1]]]},15);`)
fi:
end:

VatterByPrefix := proc(pats, pref, gaps) local one_ind, old_pref, old_gaps, i, output, prev, old, new_ind, new_out:
	option remember:
	if pref = [] then
		return(Wilf(gaps[1],pats)):
	fi:
	
	one_ind:=1:
	while pref[one_ind] <> 1 do
		one_ind:=one_ind+1:
	od:
	
	old_pref := [op(1..one_ind-1, pref), op(one_ind+1..nops(pref), pref)]:
	old_pref := redu(old_pref):
	old_gaps := [op(1..one_ind, gaps), op(one_ind+2..nops(gaps),gaps)]:
	old_gaps[one_ind] := old_gaps[one_ind] + gaps[one_ind+1]:
	
	old := VatterByPrefix(pats, old_pref, old_gaps):
	
	output := []:
	for prev in old do
		prev := prev + [1$nops(prev)]:
		new_ind := add(gaps[i],i=1..one_ind)+one_ind:
		new_out := [op(1..new_ind-1, prev), 1, op(new_ind..nops(prev),prev)]:
		if OneCheckAvoid(new_out, pats) then
			output:=[op(output), new_out]:
		fi:
	od:
	return(output):
end:
	

##Inspired by Dr. Z's procedure Wilf
Wilf:=proc(n, pats) local old, o, i, cand, output:
	if n = 0 then
		return({[]}):
	fi:
	old := Wilf(n-1, pats):
	output := []:
	for o in old do
		for i from 1 to n do
			cand := [op(1.. i-1, o), n, op(i..(n-1), o)]:
			if newCheckAvoid(cand, pats) then
				output := [op(output), cand]:
			fi:
		od:
	od:
	return(output):
end:
			

##tests to see if perm contains some pattern in pats using the largest element of perm
newCheckAvoid:=proc(perm, pats) local n, i, n_ind, inds_avail, pat, pre_inds, pre_ind, subword, k, inds:
	n := max(perm):
	for i from 1 to n do
		if perm[i] = n then
			n_ind := i:
			break:
		fi:
	od:
	inds_avail:=[op({seq(i, i=1..n)} minus {n_ind})]:
	for pat in pats do
		k := nops(pat):
		pre_inds := choose(inds_avail, k-1):
		for pre_ind in pre_inds do
			inds := sort([op(pre_ind), n_ind]):
			subword := [seq(perm[i], i in inds)]:
			if pat = redu(subword) then
				return(false):
			fi:
		od:
	od:
	return(true):
end:

##tests to see if perm contains some pattern in pats using the first element of perm
InitCheckAvoid := proc(perm, pats) local inds, ind, pat, subword:
	for pat in pats do
		inds := choose(nops(perm)-1, nops(pat) - 1):
		for ind in inds do
			subword := [perm[1], seq(perm[i+1], i in ind)]:
			if redu(subword) = pat then
				return(false):
			fi:
		od:
	od:
	return(true):
end:

OneCheckAvoid := proc(perm, pats) local new_pats, pat:
	new_pats:= {seq(Invert(pat), pat in pats)}:
	return(InitCheckAvoid(Invert(perm),new_pats)):
end:
			

redu:=proc(perm)
	sort(sort(perm, 'output=permutation'), 'output=permutation'):
end:	

Invert:=proc(perm)
	sort(perm, output = permutation):
end:

		
	
GapVectorsBySum := proc(len, tot) local output, i, old, o:
	option remember:
	if tot = 0 then
		return([[0$len]]):
	fi:
	if len = 1 then
		return([tot]):
	fi:
	output := []:
	for i from 0 to tot do
		old := GapVectorsBySum(len-1, tot-i):
		for o in old do
			output := [op(output), [i, op(o)]]:
		od:
	od:
	return(output):
end:
	
AllGapVectors := proc(len, max_tot) local output, i:
	option remember:
	output := []:
	for i from 0 to max_tot do
		output := [op(output), op(GapVectorsBySum(len, i))]:
	od:
	return(output):
end:

RestrictedGapVectors := proc(len, max_tot, restrict) local restrict_sum, i, pre_vectors, pv, output:
	option remember:
	restrict_sum := add(restrict[i], i=1..nops(restrict)):
	if restrict_sum > max_tot then
		return([])
	fi:
	pre_vectors := AllGapVectors(len, max_tot - restrict_sum):
	output := []:
	for pv in pre_vectors do
		output := [op(output), pv+restrict]:
	od:
	return(output):
end:

RestrictedGapVectorsBySum := proc(len, tot, restrict) local restrict_sum, pre_vectors, i, pv, output:
	option remember:
	restrict_sum := add(restrict[i], i=1..nops(restrict)):
	if restrict_sum > tot then
		return([]):
	fi:
	pre_vectors := GapVectorsBySum(len, tot-restrict_sum):
	output:=[]:
	for pv in pre_vectors do
		output := [op(output), pv+restrict]:
	od:
	return(output):
end:

ReverseRestrictedGapVectors := proc(len, max_tot, restricts) local cands, cand, output, forbid, restrict:
	option remember:
	cands := AllGapVectors(len, max_tot):
	forbid := {seq(op(RestrictedGapVectors(len, max_tot, restrict)), restrict in restricts)}:
	output := []:
	for cand in cands do
		if not member(cand, forbid) then
			output := [op(output), cand]:
		fi:
	od:
	return(output):
end:

Refine := proc(perm) local output, new_perm, i,j:
	output := []:
	for i from 1 to nops(perm)+1 do
		new_perm := perm:
		for j from 1 to nops(new_perm) do
			if new_perm[j] >= i then
				new_perm[j] := new_perm[j] +1:
			fi:
		od:
		output := [op(output), [op(new_perm),i]]:
	od:
	return(output):
end:	

Remove:=proc(i,perm) local j, output:
	output:=[]:
	for j from 1 to nops(perm) do
		if perm[j] <> i then
			output:=[op(output), perm[j]]:
		fi:
	od:
	output:
end:
		


#naive_scheme.txt
#same as naive_scheme.jl, but written in Maple

#read `utils.txt`:

CheckEquivalence := proc(pats, pref1, pref2, collapsed_ind, restrict, restricts_found) local vecs2check, pat, v, filler1, v2, filler2, i, cands, forbid,r:
	cands := RestrictedGapVectors(nops(pref1)+1, max([seq(nops(pat), pat in pats)])+add(restrict[i], i=1..nops(restrict))-1, restrict):
	forbid := {seq(op(RestrictedGapVectors(nops(pref1)+1, max([seq(nops(pat), pat in pats)])+add(restrict[i],i=1..nops(restrict)) -1, r)), r in restricts_found)}:
	vecs2check:=[]:
	for v in cands do
		if not member(v,forbid) then
			vecs2check:=[op(vecs2check),v]:
		fi:
	od:
		
	for v in vecs2check do
		filler1 := nops(VatterByPrefix(pats, pref1,v)):
		v2 := [op(1..collapsed_ind, v), op(collapsed_ind+2..nops(v), v)]:
		v2[collapsed_ind] := v2[collapsed_ind]+v[collapsed_ind+1]:
		filler2 := nops(VatterByPrefix(pats, pref2,v2)):
		if filler1 <> filler2 then
			return(false):
		fi:
	od:
	return(true):
end:

IsImpossible := proc(pats, pref, r, restricts_found) local r1:
	for r1 in restricts_found do
		if ArrayLess(r1,r) then
			return(true):
		fi:
	od:
	if nops(VatterByPrefix(pats, pref, r)) <> 0 then
		return(false):
	fi:
	return(true):
end:

CheckEquivalenceDown := proc(pats, pref1, pref2, collapsed_ind, restrict) local vecs2check, pat, v, filler1, v2, filler2, i:
	vecs2check := RestrictedGapVectors(nops(pref1)+1, max([seq(nops(pat), pat in pats)])+add(restrict[i],i=1..nops(restrict)) -1, restrict):
	for v in vecs2check do
		filler1:=nops(VatterByPrefix(pats, pref1,v)):
		v2 := [op(1..collapsed_ind, v), op(collapsed_ind+2..nops(v), v)]:
		v2[collapsed_ind] := v2[collapsed_ind]+v[collapsed_ind+1]:
		filler2 := nops(VatterByPrefix(pats, pref2,v2)):
		if filler1 <> filler2 then
			return(false):
		fi:
	od:
	return(true):
end:

CheckEquivalenceUp := proc(pats, pref1, pref2, collapsed_ind, restricts) local vecs2check, pat, v, filler1, v2, filler2, i:
	vecs2check := ReverseRestrictedGapVectors(nops(pref1)+1, max([seq(nops(pat), pat in pats)])-1, restricts):
	for v in vecs2check do
		filler1 := nops(VatterByPrefix(pats, pref1,v)):
		v2 := [op(1..collapsed_ind, v), op(collapsed_ind+2..nops(v), v)]:
		v2[collapsed_ind] := v2[collapsed_ind]+v[collapsed_ind+1]:
		filler2 := nops(VatterByPrefix(pats, pref2,v2)):
		if filler1 <> filler2 then
			return(false):
		fi:
	od:
	return(true):
end:

IsImpossibleDown := proc(pats, pref, r)
	if nops(VatterByPrefix(pats, pref, r)) <> 0 then
		return(false):
	fi:
	return(true):
end:

IsImpossibleUp := proc(pats, pref, r)
	if {op(r)} = {0} then
		return(true):
	fi:
	if nops(VatterByPrefix(pats,pref,[0$(nops(pref)+1)])) <> 0 then
		return(false):
	fi:
	return(true):
end:

IsLeaf1 := proc(pats, pref, max_gap_size) local pref_cands, restrict_cands, r, p, i, output, restricts_found, unreduced, redundant, r1:
	if pref = [] then
		return([false,false]):
	fi:
	output:=[pref, [], []]:
	restrict_cands := RestrictedGapVectors(nops(pref)+1, max_gap_size, [0$(nops(pref)+1)]):
	restricts_found := []:
	for r in restrict_cands do
		redundant:=false:
		for r1 in restricts_found do
			if ArrayLess(r1,r) then
				redundant := true:
				break:
			fi:
		od:
		if redundant then
			next:
		fi: 
		if IsImpossibleDown(pats,pref,r) then
			output[2]:=[op(output[2]), [r,0]]:
			restricts_found := [op(restricts_found), r]:
		else
			for i from 1 to nops(pref) do
				p:=redu([op(1..i-1,pref),op(i+1..nops(pref), pref)]):
				if CheckEquivalenceDown(pats, pref, p, i, r) then
					output[2] := [op(output[2]), [r,i]]:
					restricts_found := [op(restricts_found), r]:
					break:
				fi:
			od:
		fi:
	od:
	#print(restricts_found):
	if IsImpossibleUp(pats, pref, restricts_found) then
		output[3]:=0:
		return(output):
	fi:
	for i from 1 to nops(pref) do
		p:=redu([op(1..i-1,pref),op(i+1..nops(pref), pref)]):
		if CheckEquivalenceUp(pats, pref, p, i, restricts_found) then
			output[3] := i:
			return(output):
		fi:
	od:
	if nops(restricts_found) = 0 then
		return([false, false]):
	else
		return([false, true]):
	fi:
end:

IsLeaf2 := proc(pats, pref, max_gap_size) local pref_cands, restrict_cands, r, p, i, output, restricts_found, unreduced, redundant, r1:
	if pref = [] then
		return(false):
	fi:
	output:=[pref, []]:
	restrict_cands := Rev(RestrictedGapVectors(nops(pref)+1, max_gap_size, [0$(nops(pref)+1)])):
	restricts_found := []:
	for r in restrict_cands do
		if IsImpossible(pats,pref,r,restricts_found) then
			output[2]:=[op(output[2]), [r,0]]:
			restricts_found := [op(restricts_found), r]:
		else
			for i from 1 to nops(pref) do
				p:=redu([op(1..i-1,pref),op(i+1..nops(pref), pref)]):
				if CheckEquivalence(pats, pref, p, i, r,restricts_found) then
					output[2] := [op(output[2]), [r,i]]:
					restricts_found := [op(restricts_found), r]:
					break:
				fi:
			od:
		fi:
		if ReverseRestrictedGapVectors(nops(pref)+1, max_gap_size, restricts_found) = [] then
			return(output):
		fi:
	od:
	return(false):
end:
			
HasScheme1 := proc(pats, max_depth, max_gap_size,  current_rules,  current_need_reduction, current_found_reduction) local rules, need_reduction, found_reduction, pref, putative, opt, nr, unreduced, r:
	#print(current_rules, current_need_reduction, current_found_reduction):
	if current_need_reduction = {} then
		return(current_rules):
	fi:
	rules := current_rules:
	need_reduction := current_need_reduction:
	found_reduction := current_found_reduction:
	for pref in need_reduction do
		#print(pref):
		putative := IsLeaf1(pats, pref, max_gap_size):
		if putative[1] <> false then
			rules := rules union {putative}:
			need_reduction := need_reduction minus {pref}:
			found_reduction := found_reduction union {pref}:
			for r in putative[2] do
				if r[2] > 0 then
					need_reduction := need_reduction union {redu([op(1..r[2]-1, pref), op(r[2]+1..nops(pref),pref)])}:
				fi:
			od:
			if putative[3] > 0 then
				need_reduction := need_reduction union {redu([op(1..putative[3]-1, pref), op(putative[3]+1..nops(pref),pref)])}:
			fi:
		else
			if putative[2] <> false then
				putative:=IsLeaf2(pats,pref,max_gap_size):
				if putative <> false then
					rules := rules union {putative}:
					need_reduction := need_reduction minus {pref}:
					found_reduction := found_reduction union {pref}:
					for r in putative[2] do
						if r[2] > 0 then
							need_reduction := need_reduction union {redu([op(1..r[2]-1, pref), op(r[2]+1..nops(pref),pref)])}:
						fi:
					od:
					#if putative[3] > 0 then
					#	need_reduction := need_reduction union {redu([op(1..putative[3]-1, pref), op(putative[3]+1..nops(pref),pref)])}:
					#fi:
				else
					rules := rules union {[pref, []]}:
					need_reduction := need_reduction union {op(IRefine(pref))}:
					found_reduction := found_reduction union {pref}:
				fi:
			else
				rules := rules union {[pref, []]}:
				need_reduction := need_reduction union {op(IRefine(pref))}:
				found_reduction := found_reduction union {pref}:
			fi:
		fi:
	od:
	need_reduction := need_reduction minus found_reduction:
	if max({seq(nops(nr), nr in need_reduction)}) > max_depth then
		return([FAIL, rules, need_reduction]):
	fi:
	return(HasScheme1(pats,max_depth,max_gap_size, rules, need_reduction, found_reduction)):
end:

#Depth-first search yields quicker failures
HasScheme2 := proc(pats, max_depth, max_gap_size, current_rules, current_need_reduction, current_found_reduction) local rules, need_reduction, found_reduction, max_len, max_pref, pref, putative, opt, r:
	rules:=current_rules:
	need_reduction := current_need_reduction:
	found_reduction:= current_found_reduction:
	while need_reduction <> {} do
		max_len := 0:
		max_pref:=[]:
		for pref in need_reduction do
			if nops(pref) > max_len then
				max_len:=nops(pref):
				max_pref:=pref:
			fi:
		od:
		if max_len > max_depth then
			return([FAIL,{max_pref}]):
		fi:
		putative := IsLeaf1(pats, pref, max_gap_size):
		if putative[1] <> false then
			rules := rules union {putative}:
			need_reduction := need_reduction minus {pref}:
			found_reduction := found_reduction union {pref}:
			for r in putative[2] do
				if r[2] > 0 then
					need_reduction := need_reduction union {redu([op(1..r[2]-1, pref), op(r[2]+1..nops(pref),pref)])}:
				fi:
			od:
			if putative[3] > 0 then
				need_reduction := need_reduction union {redu([op(1..putative[3]-1, pref), op(putative[3]+1..nops(pref),pref)])}:
			fi:
		else
			if putative[2] <> false then
				putative:=IsLeaf2(pats,pref,max_gap_size):
				if putative <> false then
					rules := rules union {putative}:
					need_reduction := need_reduction minus {pref}:
					found_reduction := found_reduction union {pref}:
					for r in putative[2] do
						if r[2] > 0 then
							need_reduction := need_reduction union {redu([op(1..r[2]-1, pref), op(r[2]+1..nops(pref),pref)])}:
						fi:
					od:
					#if putative[3] > 0 then
					#	need_reduction := need_reduction union {redu([op(1..putative[3]-1, pref), op(putative[3]+1..nops(pref),pref)])}:
					#fi:
				else
					rules := rules union {[pref,[]]}:
					need_reduction := need_reduction union {op(IRefine(pref))}:
					found_reduction := found_reduction union {pref}:
				fi:
			else
				rules := rules union {[pref,[]]}:
				need_reduction := need_reduction union {op(IRefine(pref))}:
				found_reduction := found_reduction union {pref}:
			fi:
		fi:
		need_reduction := need_reduction minus found_reduction:
	od:
	return(rules):
end:
		
HasScheme := proc(pats, max_depth, max_gap_size, quick_fail:=false)
	if quick_fail then
		HasScheme2(pats, max_depth, max_gap_size, {}, {[]}, {}):
	else
		HasScheme1(pats, max_depth, max_gap_size, {}, {[]}, {}):
	fi:
end:

SeqS := proc(scheme, n) local n1:
	if type(scheme, `table`) then
		seq(FindTermTable(scheme,[],[n1]), n1=1..n):
	else
		seq(FindTermScheme(scheme,[],[n1]), n1=1..n):
	fi:
end:

FindTermTable:=proc(T, pref, gap) local rule, gap_restricts, i, refi, r, new_ind, new_gap, j, output, g, down_collapsed_ind, down_pref, down_gap, up_collapsed_ind, up_pref, up_gap:
	option remember:
	 
	rule := T[pref]:
	gap_restricts := rule[1]:
	if gap_restricts = [] then
		if add(gap[i], i=1..nops(gap)) = 0 then
			return(1):
		fi:
		output := 0:
		refi := IRefine(pref):
		for r in refi do
			new_ind := max[index](r):
			new_gap := [0$(nops(gap)+1)]:
			for j from 1 to gap[new_ind] do
				new_gap := [op(1..new_ind-1, gap), j-1, gap[new_ind]-j, op(new_ind+1..nops(gap), gap)]:
				output := output + FindTermTable(T, r, new_gap):
			od:
		od:
		return(output):
	fi:
	for i from 1 to nops(gap_restricts) do
		g:=gap_restricts[i]:
		if ArrayLess(g[1], gap) then
			down_collapsed_ind := g[2]:
			#this is a convention to indicate nonviable gap vectors
			if down_collapsed_ind = 0 then
				return(0):
			fi:
					
			down_pref := redu([op(1..down_collapsed_ind-1, pref), op(down_collapsed_ind+1..nops(pref),pref)]):
					
			down_gap := [op(1..down_collapsed_ind, gap),op(down_collapsed_ind+2..nops(gap), gap)]:
			down_gap[down_collapsed_ind]:= down_gap[down_collapsed_ind] + gap[down_collapsed_ind+1]:
			return(FindTermTable(T, down_pref, down_gap)):
		fi:
	od:
	
	up_collapsed_ind := rule[2]:
	if up_collapsed_ind = 0 then
		return(0):
	fi:
	up_pref := redu([op(1..up_collapsed_ind-1, pref), op(up_collapsed_ind+1..nops(pref), pref)]):
	up_gap := [op(1..up_collapsed_ind, gap), op(up_collapsed_ind+2..nops(gap), gap)]:
	up_gap[up_collapsed_ind] := up_gap[up_collapsed_ind] + gap[up_collapsed_ind+1]:
	return(FindTermTable(T, up_pref, up_gap)):
end:
	
	

FindTermScheme:=proc(scheme, pref, gap) local i,j,rule,gap_restricts, down_collapsed_ind, down_pref, down_gap, up_collapsed_ind, up_pref, up_gap, output, refi, r, new_ind, new_gap,g:
	option remember:
	#print(pref,gap):
	
	
	for rule in scheme do
		if pref = rule[1] then
			gap_restricts := rule[2]:
			
			if gap_restricts = [] then
				if add(gap[i], i=1..nops(gap)) = 0 then
					return(1):
				fi:
				output:=0:
				refi := IRefine(pref):
				for r in refi do
					new_ind := max[index](r):
					new_gap := [0$(nops(gap)+1)]:
					for j from 1 to gap[new_ind] do
						new_gap := [op(1..new_ind-1, gap), j-1, gap[new_ind]-j, op(new_ind+1..nops(gap), gap)]:
						output := output + FindTermScheme(scheme, r, new_gap):
					od:
				od:
				return(output):
			fi:

			for i from 1 to nops(gap_restricts) do
				g:=gap_restricts[i]:
				if ArrayLess(g[1], gap) then
					down_collapsed_ind := g[2]:
					
					#this is a convention to indicate nonviable gap vectors
					if down_collapsed_ind = 0 then
						return(0):
					fi:
					
					down_pref := redu([op(1..down_collapsed_ind-1, pref), op(down_collapsed_ind+1..nops(pref),pref)]):
					
					down_gap := [op(1..down_collapsed_ind, gap),op(down_collapsed_ind+2..nops(gap), gap)]:
					down_gap[down_collapsed_ind]:= down_gap[down_collapsed_ind] + gap[down_collapsed_ind+1]:
					return(FindTermScheme(scheme, down_pref, down_gap)):
				fi:
			od:
			
			up_collapsed_ind := rule[3]:
			if up_collapsed_ind = 0 then
				return(0):
			fi:
			up_pref := redu([op(1..up_collapsed_ind-1, pref), op(up_collapsed_ind+1..nops(pref), pref)]):
			up_gap := [op(1..up_collapsed_ind, gap), op(up_collapsed_ind+2..nops(gap), gap)]:
			up_gap[up_collapsed_ind] := up_gap[up_collapsed_ind] + gap[up_collapsed_ind+1]:
			return(FindTermScheme(scheme, up_pref, up_gap)):
		fi:
	od:
end:
	
ArrayLess:=proc(arr1, arr2) local i:
	if nops(arr1) <> nops(arr2) then
		 print(`Two arrays must have the same size`):
		 return(FAIL):
	fi:
	for i from 1 to nops(arr1) do
		if arr1[i] > arr2[i] then
			return(false):
		fi:
	od:
	return(true):
end:

Scheme2Table:=proc(scheme) local T,rule:
	T:=table():
	for rule in scheme do
		T[rule[1]]:=[op(2..-1,rule)]:
	od:
	return(op(T)):
end:
	

IRefine:=proc(perm) local i,n:
	n:=nops(perm):
	return([seq([op(1..i,perm),n+1,op(i+1..n, perm)], i=0..n)]):
end:

##Example Scheme for 123: {[[1, 2], [[[0, 0, 1], 1, []]], [2, [1]]], [[2, 1], [[[0, 0, 0], 2, [1]]], [1, [1]]]}

DecreaseToOne:=proc(perm) local i:
	for i from 1 to nops(perm) do
		if perm[i] = 1 then
			return(true):
		elif i >=2 and perm[i] > perm[i-1] then
			return(false):
		fi:
	od:
	6:
end:

Rev:=proc(L) local i:
	[seq(L[-i],i=1..nops(L))]:
end:

downfix := proc(len,pi) local i, output:
	output := []:
	for i from 1 to nops(pi) do
		if pi[i] <= len then
			output := [op(output), pi[i]]:
		fi:
	od:
end: