#Flexible_Vincular_Scheme
with(combinat):

IsLeaf := proc(bad_pats, downfix, max_gap) 
	local pi, Y, vector_len, poss_gap_restrictions, real_gap_restrictions, restrict,
	redundant, prev_restrict, reduction_rules, found_new_restrict, scenarios_with_i, 
	scenarios_without_i, deletable, scenario, i:
	pi := downfix[1]:
	Y := downfix[2]:
	vector_len := nops(pi)+1:
	poss_gap_restrictions := AllGapVectors(vector_len, max_gap):
	real_gap_restrictions := {}:
	reduction_rules := []:
	found_new_restrict := true:
	while (not member([0$nops(pi)+1],real_gap_restrictions)) and found_new_restrict do
		found_new_restrict := false:
		for restrict in poss_gap_restrictions do
			redundant := false:
			for prev_restrict in real_gap_restrictions do
				if Satisfies(restrict, prev_restrict) then
					redundant := true:
					break:
				fi:
			od:
			if redundant then
				continue:
			fi:

			if not redundant then
				if RestrictIsNonviable(bad_pats, downfix, restrict) then
					real_gap_restrictions := real_gap_restrictions union {restrict}:
					reduction_rules := [op(reduction_rules), [restrict, 0]]:
					poss_gap_restrictions := poss_gap_restrictions minus {restrict}:
					found_new_restrict := true:
				else
					for i from 1 to nops(pi) do
						scenarios_with_i := Sp(bad_pats, downfix, i, restrict, real_gap_restrictions):
						scenarios_without_i := map(curry(DeleteElt,i), scenarios_with_i):
						deletable := true:
						for scenario in scenarios_without_i do
							if AvoidsAll(scenario, bad_pats) then
								deletable := false:
								break:
							fi:
						od:
						if deletable then
							real_gap_restrictions := real_gap_restrictions union {restrict}:
							poss_gap_restrictions := poss_gap_restrictions minus {restrict}:
							reduction_rules := [op(reduction_rules),[restrict, i]]:
							found_new_restrict := true:
							break:
						fi:
					od:
				fi:
			fi:
			if member(restrict, real_gap_restrictions) then
				break:
			fi:
		od:
	od:
	if member([0$nops(pi)+1],real_gap_restrictions) then
		return(reduction_rules):
	else
		return(false):
	fi:
end:

BPPats2Ypats := proc(pre_bad_pats) local bad_pats, pat:
	bad_pats := {}:
		for pat in pre_bad_pats do
			bad_pats := bad_pats union {[pat[1], select(i -> pat[2][i]=0,{seq(i,i=1..nops(pat[2]))})]}:
		od:
	return(bad_pats):
end:

HasScheme := proc(pre_bad_pats, max_depth, max_gap)
	local need_reduction, has_reduction, partial_scheme, downfix, putative, deletable,
	rule, bad_pats, pat, i:
	#test to see if pattern is given in my format or Baxter & Pudwell's 
	if type(pre_bad_pats[1][2], `list`) then
		bad_pats := BPPats2Ypats(pre_bad_pats):
	else
		bad_pats := pre_bad_pats:
	fi:
	need_reduction := {[[],{}]}:
	has_reduction := {}:
	partial_scheme := {}:
	while need_reduction <> {} do
		downfix := need_reduction[1]:
		putative := IsLeaf(bad_pats, downfix, max_gap):
		if putative = false then
			if nops(downfix[1]) > max_depth then
				return(false):
			fi:
			has_reduction := has_reduction union {downfix}:
			need_reduction := need_reduction union Refinements(downfix):
			partial_scheme := partial_scheme union {[downfix, []]}
		else
			has_reduction := has_reduction union {downfix}:
			for rule in putative do
				deletable := rule[2]:
				if deletable <> 0 then
					need_reduction := need_reduction union {DeleteElt(deletable,downfix)}:
				fi:
			od:
			partial_scheme := partial_scheme union {[downfix, putative]}:
		fi:
		need_reduction := need_reduction minus has_reduction:
	od:
	return(partial_scheme):
end:

SeqS1 := proc(scheme, downfix, gap_vector)
	local pi, Y, k, downfix_rule, rule, output, i, new_pi, gap_split, new_gap_vector,
	gap_restrict_subrule, deletable, new_downfix, deletable_index:
	option remember:
	#print(gap_vector,downfix):
		pi := downfix[1]:
		Y := downfix[2]:
		k := nops(pi):
		#select the rule from the scheme corresponding to the downfix
		downfix_rule := select(rule -> rule[1] = downfix, scheme)[1]:
		if {op(gap_vector)} = {0} and (not member([gap_vector,0], {op(downfix_rule[2])})) then
			return(1):
		fi:
		if downfix_rule[2] = [] then
			output := 0:
			for i from 1 to k+1 do
				new_pi := [op(1..i-1, pi), k+1, op(i..k, pi)]:
				for gap_split from 0 to gap_vector[i]-1 do
					new_gap_vector := [op(1..i-1, gap_vector), gap_split, gap_vector[i]-gap_split-1, op(i+1..k+1, gap_vector)]:
					output := output + SeqS1(scheme, [new_pi, Y], new_gap_vector):
				od:
			od:
			return(output):
		else
			#select the subrule corresponding to the first satisfied gap restriction
			gap_restrict_subrule := select(subrule -> Satisfies(gap_vector,subrule[1]), downfix_rule[2])[1]:
			deletable := gap_restrict_subrule[2]:
			if deletable = 0 then
				return(0):
			else
				new_downfix := DeleteElt(deletable, downfix):
				deletable_index := select(i -> pi[i] = deletable, [seq(i,i=1..k)])[1]:
				new_gap_vector := [op(1..deletable_index-1,gap_vector), gap_vector[deletable_index]+gap_vector[deletable_index+1], op(deletable_index+2..k+1,gap_vector)]:
				return(SeqS1(scheme, new_downfix, new_gap_vector)):
			fi:
		fi:
end:

SeqS := proc(scheme, N) local n:
	return(seq(SeqS1(scheme, [[],{}], [n]), n=1..N)):
end:


Sp1 := proc(bad_pat, downfix, ind, gap_restrict, prev_gap_restrict) 
	option remember:
	local pi, Y, k, other_inds, sub_pi, sub_pi_sets, p, sig, X, vals, poss_downfix, 
		all_downfixes, output, downfix_overlap, bins_of_downfix_elts, bins_of_occurrence_elts,
		pi_overlap,i,raised_sig,sig_overlap,new_scenarios, partial_scenarios, complete_scenarios,
		ps, gap_vector, filler_needed, amt_filler, filler_cands, filler_choices, choice,
		raised_filler, raised_ps, ordered_filler, bins_of_filler, bins_of_content, content, 
		combined_content, combined_bins, scenario, raised_X, scenario_Y, all_occurrence_elts,
		partially_raised_X:
	pi := downfix[1]:
	Y := downfix[2]:
	sig := bad_pat[1]:
	X := bad_pat[2]:
	k := nops(pi):
	sub_pi_sets := AllSublistsWithI(pi, ind):
	all_downfixes := {}:
	for sub_pi in sub_pi_sets do
		vals := {op(sub_pi)}:
		#print(vals):
		poss_downfix := [sub_pi, ((vals intersect Y) union (vals minus {seq(i-1, i in vals)} minus {k}))]:
		if IsRepDownfix(poss_downfix, bad_pat) then
			all_downfixes := all_downfixes union {poss_downfix}:
		fi:
	od:
	#print(all_downfixes):
	partial_scenarios := {}:
	complete_scenarios := {}:
	for downfix_overlap in all_downfixes do
		pi_overlap := downfix_overlap[1]:
		raised_sig := [seq(s+nops(pi)-nops(pi_overlap),s in sig)]:
		sig_overlap := [seq(s+nops(pi) - nops(pi_overlap), s in redu(pi_overlap))]:
		#print(raised_sig, sig_overlap):
		bins_of_downfix_elts := [select(curry(OccursBetween,pi,0, pi_overlap[1]), pi)]:
		bins_of_occurrence_elts := [select(curry(OccursBetween, raised_sig, 0, sig_overlap[1]), raised_sig)]:
		
		for i from 2 to nops(pi_overlap) do
			bins_of_downfix_elts := [op(bins_of_downfix_elts), select(curry(OccursBetween, pi, pi_overlap[i-1], pi_overlap[i]), pi)]:
			bins_of_occurrence_elts := [op(bins_of_occurrence_elts), select(curry(OccursBetween, raised_sig, sig_overlap[i-1], sig_overlap[i]), raised_sig)]:
		od:
		bins_of_downfix_elts := [op(bins_of_downfix_elts),select(curry(OccursBetween,pi,pi_overlap[-1], 0), pi)]:
		bins_of_occurrence_elts := [op(bins_of_occurrence_elts),select(curry(OccursBetween,raised_sig,sig_overlap[-1], 0), raised_sig)]:
		#print(pi_overlap, bins_of_occurrence_elts):
		all_occurrence_elts := sort([op(pi_overlap),seq(op(bin),bin in bins_of_occurrence_elts)]):
		#print(all_occurrence_elts):
		#print(bins_of_occurrence_elts):
		partial_scenarios := AllMerges(bins_of_downfix_elts[1],bins_of_occurrence_elts[1]):
		for i from 2 to nops(pi_overlap)+1 do
			partial_scenarios := {seq(seq([op(p), pi_overlap[i-1], op(m)], p in partial_scenarios), m in AllMerges(bins_of_downfix_elts[i],bins_of_occurrence_elts[i]))}:
		od:
		#print(partial_scenarios):
		for ps in partial_scenarios do
			#print(ps):
			gap_vector := GetGapVector(ps, pi):
			if member(true, map(curry(Satisfies, gap_vector), prev_gap_restrict)) then
				next:
			fi:
			#print(gap_vector):
			filler_needed := [seq(max(0,gap_restrict[i]-gap_vector[i]), i=1..nops(gap_vector))]:
			amt_filler := convert(filler_needed,`+`):
			filler_cands := ({seq(i,i=1..nops(ps)+1)} minus {op(pi)}) minus {seq(x+nops(ps)-nops(sig)+1, x in X)}:
			filler_choices := ChooseWithReplacement(filler_cands, amt_filler):
			#print(filler_choices):
			for choice in filler_choices do
				raised_filler := [seq(choice[i]+i-1, i=1..nops(choice))]:
				raised_ps := [seq(ps[i] + nops(select(x->evalb(x<=ps[i]), choice)),i=1..nops(ps))]:
				#print(raised_ps):
				for ordered_filler in permute(raised_filler) do
					#print(ordered_filler):
					bins_of_filler := PutInBins(ordered_filler, filler_needed):
					content := select(x->evalb(x>nops(pi)), raised_ps):
					bins_of_content := PutInBins(content, gap_vector):
					#print(bins_of_filler, bins_of_content):
					for combined_content in MergeByBins(bins_of_filler, bins_of_content) do 
						combined_bins := PutInBins(combined_content, [seq(gap_vector[i]+filler_needed[i],i=1..nops(gap_vector))]):
						scenario := combined_bins[1]:
						for i from 1 to nops(pi) do
							scenario := [op(scenario), pi[i], op(combined_bins[i+1])]:
						od:
						partially_raised_X := [seq(all_occurrence_elts[i],i in X)]:
						raised_X := {seq(x + nops(select(y->evalb(y<=x), choice)),x in partially_raised_X)}:
						scenario_Y := {op(scenario)} minus ({seq(i,i=1..nops(pi)-1)} minus Y) minus raised_X minus {nops(scenario)}:
						#print(ps, scenario,scenario_Y):
						complete_scenarios := complete_scenarios union {[scenario,scenario_Y]}:
					od:
				od:
			od:
		od:
	od:
	return(complete_scenarios):
end:

Sp := proc(bad_pats, downfix, ind, gap_restrict, prev_gap_restrict) local all_scenarios, bad_pat:
	all_scenarios := {}:
	for bad_pat in bad_pats do
		all_scenarios := all_scenarios union Sp1(bad_pat, downfix, ind, gap_restrict, prev_gap_restrict):
	od:
	return(all_scenarios):
end:

DeleteElt := proc(elt, perm) local pi, Y, elt_ind, new_pi, new_Y, f:
	pi := perm[1]:
	Y := perm[2]:
	new_pi := redu(select(i -> i <> elt, pi)):
	f := proc(i)
		if i >= elt then
			i-1:
		else
			i:
		fi:
	end:
	new_Y := map(f, Y) union {elt-1} minus {0}:
	return([new_pi,new_Y]):
end:

Refinements := proc(downfix) local i, refinements, k, pi, Y, new_pi:
	pi := downfix[1]:
	Y := downfix[2]:
	k := nops(pi):
	refinements := {}:
	for i from 1 to k+1 do
		new_pi := [op(1..i-1, pi), k+1, op(i..k, pi)]:
		refinements := refinements union {[new_pi, Y]}:
	od:
	return(refinements):
end:


PermsByDownfixGaps := proc(downfix, gap_vector) 
	local pi, Y, gap_norm, n, new_elt_ords, new_Y, all_perms, ordered_elts, perm, i, prev_used:
	pi := downfix[1]:
	Y := downfix[2]:
	gap_norm := convert(gap_vector, `+`):
	n := nops(pi) + gap_norm:
	new_elt_ords := permute([seq(i,i=nops(pi)+1..n)]):
	new_Y := Y union {seq(i,i=nops(pi)..n-1)}:
	all_perms := {}:
	for ordered_elts in new_elt_ords do
		perm := [op(1..gap_vector[1], ordered_elts)]:
		for i from 1 to nops(pi) do
			prev_used := add(gap_vector[j],j=1..i):
			perm := [op(perm), pi[i], op(prev_used+1..prev_used+gap_vector[i+1], ordered_elts)]:
		od:
		all_perms := all_perms union {[perm, new_Y]}:
	od:
	return(all_perms):
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 := output union {[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 := output union GapVectorsBySum(len, i):
	od:
	return(output):
end:

RestrictIsNonviable := proc(bad_pats, downfix, gap_restrict) local perms, perm, pat, perm_allowed:
	perms := PermsByDownfixGaps(downfix, gap_restrict):
	for perm in perms do
		perm_allowed := true:
		for pat in bad_pats do
			if not Avoids(perm, pat) then
				perm_allowed := false:
				break:
			fi:
		od:
		if perm_allowed then
			return(false):
		fi:
	od:
	return(true):
end:


AllSublists := proc(li) local prev_lists, prev_list:
	if li = [] then
		return({[]}):
	else
		prev_lists := AllSublists([op(2..-1,li)]):
		return(prev_lists union {seq([li[1], op(prev_list)], prev_list in prev_lists)}):
	fi:
end:
	
AllSublistsWithI := proc(li, i) local tail, li1:
	if li = [] then
		return({}):
	elif li[1] = i then
		tail := [op(2..-1, li)]:
		return({seq([i, op(li1)], li1 in AllSublists(tail))}):
	else
		tail := [op(2..-1, li)]:
		return({seq([li[1], op(li1)], li1 in AllSublistsWithI(tail,i))} union {seq([op(li1)], li1 in AllSublistsWithI(tail,i))})
	fi:
end:

AllMerges := proc(li1, li2) local old_merges1, old_merges2:
	if li1 = [] then
		return({li2}):
	elif li2 = [] then
		return({li1}):
	else
		old_merges1 := AllMerges([op(2..-1,li1)],li2):
		old_merges2 := AllMerges(li1,[op(2..-1,li2)]):
		return({seq([li1[1], op(rest)], rest in old_merges1)} union {seq([li2[1], op(rest)], rest in old_merges2)}):
	fi:
end:

MergeByBins := proc(li1, li2) local first_merge, rest_merge:
	if nops(li1)<>nops(li2) then
		print(`Lists must have the same number of bins to merge`):
		return([]):
	fi:
	if nops(li1) = 1 then
		return(AllMerges(li1[1],li2[1])):
	else
		first_merge := AllMerges(li1[1],li2[1]):
		rest_merge := MergeByBins([op(2..-1,li1)], [op(2..-1,li2)]):
		return({seq(seq([op(fm), op(rm)], fm in first_merge), rm in rest_merge)}):
	fi:
end:

PutInBins := proc(li, bin_sizes) local new_li, new_sizes:
	if nops(bin_sizes) = 1 then
		return([li]):
	else
		new_li := [op(bin_sizes[1]+1..-1, li)]:
		new_sizes := [op(2..-1,bin_sizes)]:
		return([[op(1..bin_sizes[1],li)],op(PutInBins(new_li,new_sizes))]):
	fi:
end:

GetDownfix := proc(perm,level) local x,pi,Y:
	pi := perm[1]:
	Y := perm[2]:
	return([select(x->evalb(x<=level),pi), select(x->evalb(x<=level),Y)]):
end:

GetGapVector := proc(perm, downfix) local output, current_gap, i:
	output := [0$(nops(downfix)+1)]:
	current_gap := 1:
	for i from 1 to nops(perm) do
		if perm[i] <> downfix[current_gap] then
			output[current_gap] := output[current_gap] + 1:
		elif current_gap < nops(downfix) then
			current_gap := current_gap + 1:
		else
			output[-1] := nops(perm) - i:
			return(output):
		fi:
	od:
end:

Satisfies := proc(gap_vector, gap_restrict) local i:
	for i from 1 to nops(gap_vector) do
		if gap_vector[i] < gap_restrict[i] then
			return(false):
		fi:
	od:
	return(true):
end:

ChooseWithReplacement := proc(li,num)
	if num = 0 then
		return({[]}):
	else
		return({seq(seq(sort([li[i],op(li1)]), i=1..nops(li)), li1 in ChooseWithReplacement(li,num-1))}):
	fi:
end:

OccursBetween := proc(li,lower,upper,target) local reduced:
	if lower = target or upper = target then
		return(false):
	fi:
	if lower <> 0 and upper <> 0 then
		reduced := select(x->member(x,{lower,upper,target}), li):
		return(evalb(target = reduced[2] and lower=reduced[1])):
	elif upper <> 0 then
		reduced := select(x->member(x,{upper,target}), li):
		return(evalb(target = reduced[1])):
	elif lower <> 0 then
		reduced := select(x->member(x,{lower,target}),li):
		return(evalb(target = reduced[2])):
	else
		return(true):
	fi:
end:


IsRepDownfix := proc(downfix, pat) local small_pat:
	small_pat := GetDownfix(pat,nops(downfix[1])):
	return(evalb(not Avoids(downfix, small_pat))):
end:


Wilf:=proc(n, pre_bad_pats) local old, o, i, cand, output, bad_pats:
	if type(pre_bad_pats[1][2], `list`) then
		bad_pats := BPPats2Ypats(pre_bad_pats):
	else
		bad_pats := pre_bad_pats:
	fi:
	if n = 0 then
		return({[]}):
	fi:
	old := Wilf(n-1, bad_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, bad_pats) then
				output := [op(output), cand]:
			fi:
		od:
	od:
	return(output):
end:

Avoids := proc(perm, pat) local pi, Y, sig, X, n, inds, ind, new_pi, vals, new_Y, compare:
	pi := perm[1]:
	Y := perm[2]:
	sig := pat[1]:
	X := pat[2]:
	n := nops(pi):
	inds := choose([seq(i, i=1..n)], nops(sig)):
	for ind in inds do
		new_pi := pi[ind]:
		vals := {op(new_pi)}:
		new_Y := Y intersect vals:
		compare := redu_spaced([new_pi, new_Y]):
		if compare[1] = sig and (X intersect compare[2] = {}) then
			return(false):
		fi:
	od:
	return(true):
end:

AvoidsAll := proc(perm, pats) local pat:
	for pat in pats do
		if not Avoids(perm, pat) then
			return(false):
		fi:
	od:
	return(true):
end:

newCheckAvoid:=proc(perm, pats) local n, i, n_ind, inds_avail, pat, pre_inds, pre_ind, subword, k, inds, sig, X, covinc_sat,f,j:
	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
		sig := pat[1]:
		X := pat[2]:
		k := nops(sig):
		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 sig <> redu(subword) then
				next:
			fi:
			covinc_sat := true:
			for i from 1 to k do
				#assume that X can't contain k
				if member(sig[i], X) then
					f := x -> evalb(sig[x] = sig[i]+1):
 					j := select(f,[seq(l,l=1..k)])[1]:
 					if subword[j] <> subword[i]+1 then
 						covinc_sat:=false:
 					fi:
 				fi:
 			od:
 			if covinc_sat then
 				return(false):
 			fi:

		od:
	od:
	return(true):
end:

index_one_larger := proc(x,sig)
	if sig[x] = sig[i] + 1 then
		true:
	else
		false:
	fi:
end:

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

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

redu_spaced := proc(perm) local pi, Y, ordered_vals, new_pi, new_Y, i:
	pi := perm[1]:
	Y := perm[2]:
	ordered_vals := sort(pi):
	new_pi := sort(sort(pi, 'output=permutation'), 'output=permutation'):
	new_Y := {}:
	for i from 1 to nops(pi) do
		if member(ordered_vals[i], Y) or (i < nops(pi) and ordered_vals[i+1]<>ordered_vals[i]+1) then
			new_Y := new_Y union {i}:
		fi:
	od:
	return([new_pi, new_Y]):
end:

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

PatInvert:=proc(pat) local sig, X:
	sig := pat[1]:
	X := pat[2]:
	return(Invert(sig),X):
end:

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:

AllPerms := proc(downfix, gap_vector) 
	local pi, Y, k, gap_len, all_elts, all_orders, output, order, bins_of_elts, 
	i, new_perm:
	pi := downfix[1]:
	Y := downfix[2]:
	k := nops(pi):
	gap_len := convert(gap_vector, `+`):
	all_elts := [seq(i+k, i=1..gap_len)]:
	all_orders := permute(all_elts):
	output := {}:
	for order in all_orders do
		bins_of_elts := PutInBins(order, gap_vector):
		new_perm := bins_of_elts[1]:
		for i from 1 to k do
			new_perm := [op(new_perm), pi[i], op(bins_of_elts[i+1])]:
		od:
		output:= output union {[new_perm, Y]}:
	od:
	return(output):
end:

SeqS1BF := proc(pre_bad_pats, downfix, gap_vector) local universe, bad_pats:
	option remember:
	if type(pre_bad_pats[1][2], `list`) then
		bad_pats := BPPats2Ypats(pre_bad_pats):
	else
		bad_pats := pre_bad_pats:
	fi:
	universe := AllPerms(downfix, gap_vector):
	return(nops(select(rcurry(AvoidsAll,bad_pats), universe))):
end:

ContainsSubsequence := proc(big, little) local all_inds, ind:
	all_inds := choose(nops(big), nops(little)):
	for ind in all_inds do
		if big[ind] = little then
			return(true):
		fi:
	od:
	return(false):
end: