with(combinat):

Help:=proc()
print("binarytrees(n), stricttrees(n), readtree(L), readlayers(L), red(L)"):
print("pcontainsq(p,q), pavoidsq(p,q), pavoidsQset(p,Q)"):
print("treesavoidQset(n,Q), numavoidersQset(n,Q)"):
print("stricttreesavoidQset(n,Q), strictnumavoidersQset(n,Q)"):
end:

#binary tree data structure is [r,A,B] where r is the root, A is the left subtree and B is the right subtree
#[r,A] is a tree where the root has only one subtree
#[r] is a one leaf tree
#[r,[a],[]] is a two leaf tree
#[r,[a],[b]] and [r,[a,[b]]] are the two possible two-leaf tree structures
#the number of such trees given increasing labels is given by OEIS A080635
binarytrees:=proc(n) local newlabels,S,L,oldlabels,i,l,j,leftlabels,rightlabels,lefttrees,righttrees,possleftlabels,temp, lt, rt,oldleftlabels,oldrightlabels: option remember:
if n=1 then return {[1]}: fi:

S:={}:
newlabels:=[$2..n]:

#case 1: the root has 1 subtree
L:=binarytrees(n-1):
oldlabels:=[$1..n-1]:


for l in L do
S:=S union {[1,subs({seq(oldlabels[i]=newlabels[i],i=1..n-1)},l)]}:
od:

#case 2: the root has 2 subtrees
#left subtree has i vertices (1<=i<=n-2)
#right subtree has n-1-i vertices


for i from 1 to n-2 do
lefttrees:=binarytrees(i):
righttrees:=binarytrees(n-1-i):

possleftlabels:=choose({$2..n},i):

for temp in possleftlabels do
rightlabels:={$2..n} minus temp:
leftlabels:=sort([op(temp)]):
rightlabels:=sort([op(rightlabels)]):
oldleftlabels:=[$1..i]:
oldrightlabels:=[$1..n-1-i]:


for lt in lefttrees do
for rt in righttrees do

S:=S union {[1,subs({seq(oldleftlabels[j]=leftlabels[j],j=1..i)},lt), subs({seq(oldrightlabels[j]=rightlabels[j],j=1..n-i-1)},rt)]}:

od: #undo lefttrees
od: #undo righttrees

od: #undo leftlabels
od: #undo i

S:
end:

#strict tree data structure is [r,A,B] where r is the root, A is the left subtree and B is the right subtree
#[r] is a one leaf tree
#[r,[a],[b]] and [r,[a,[b]]] are the two possible two-leaf tree structures
#the number of such trees given increasing labels is given by OEIS A080635
stricttrees:=proc(n) local newlabels,S,L,oldlabels,i,l,j,leftlabels,rightlabels,lefttrees,righttrees,possleftlabels,temp, lt, rt,oldleftlabels,oldrightlabels: option remember:
if n=1 then return {[1]}: fi:

S:={}:
newlabels:=[$2..n]:

#the root has 2 subtrees
#left subtree has i vertices (1<=i<=n-2)
#right subtree has n-1-i vertices


for i from 1 to n-2 do
lefttrees:=stricttrees(i):
righttrees:=stricttrees(n-1-i):

possleftlabels:=choose({$2..n},i):

for temp in possleftlabels do
rightlabels:={$2..n} minus temp:
leftlabels:=sort([op(temp)]):
rightlabels:=sort([op(rightlabels)]):
oldleftlabels:=[$1..i]:
oldrightlabels:=[$1..n-1-i]:


for lt in lefttrees do
for rt in righttrees do

S:=S union {[1,subs({seq(oldleftlabels[j]=leftlabels[j],j=1..i)},lt), subs({seq(oldrightlabels[j]=rightlabels[j],j=1..n-i-1)},rt)]}:

od: #undo lefttrees
od: #undo righttrees

od: #undo leftlabels
od: #undo i

S:
end:


#readtree: inputs a labeled tree and outputs and corresponding permutation read from
#taking each layer beginning with the root and reading from left to right.
readtree:=proc(L) local P:
P:=readlayers(L):

[seq(op(P[i]),i=1..nops(P))]:

end:

#same as readtree, but groups labels by their depth in the tree
#output looks like readtree but with extra braces
readlayers:=proc(L) local L2,leftlayers,rightlayers,i:
if nops(L)=1 then return [L]: fi:
if nops(L)=2 then return [[L[1]],op(readlayers(L[2]))]: fi:

leftlayers:=readlayers(L[2]):
rightlayers:=readlayers(L[3]):

L2:=[[L[1]]]:
for i from 1 to min(nops(leftlayers),nops(rightlayers)) do:
L2:=[op(L2),[op(leftlayers[i]),op(rightlayers[i])]]:
od:

if nops(leftlayers)>min(nops(leftlayers),nops(rightlayers)) then
L2:=[op(L2),seq(leftlayers[j],j=min(nops(leftlayers),nops(rightlayers))+1..nops(leftlayers))]:
fi:

if nops(rightlayers)>min(nops(leftlayers),nops(rightlayers)) then
L2:=[op(L2),seq(rightlayers[j],j=min(nops(leftlayers),nops(rightlayers))+1..nops(rightlayers))]:
fi:

L2:
end:

#given a list L without repeated digits, returns the reduction of L
red:=proc(L) local A:
A:=sort(L):
subs({seq(A[i]=i,i=1..nops(L))},L):
end:


pcontainsq:=proc(p,q) local poss, temp,pos1:
poss:=choose([$1..nops(p)],nops(q)):
for pos1 in poss do
temp:=[seq(p[pos1[j]],j=1..nops(q))]:
if evalb(red(temp)=red(q)) then return true: fi:
od:
false:
end:

pavoidsq:=proc(p,q)
not pcontainsq(p,q):
end:

pavoidsQset:=proc(p,Q)
evalb({seq(pavoidsq(p,q),q in Q)}={true}):
end:

treesavoidQset:=proc(n,Q) local B,S,b:
B:=binarytrees(n):
S:={}:

for b in B do
if pavoidsQset(readtree(b),Q) then S:=S union {b}: fi:
od:

S:
end:

numavoidersQset:=proc(n,Q)
nops(treesavoidQset(n,Q)):
end:



stricttreesavoidQset:=proc(n,Q) local B,S,b:
B:=stricttrees(n):
S:={}:

for b in B do
if pavoidsQset(readtree(b),Q) then S:=S union {b}: fi:
od:

S:
end:

strictnumavoidersQset:=proc(n,Q)
nops(stricttreesavoidQset(n,Q)):
end: