1package $StatTrees; 2 3export StagedTrees; 4export PrintTree; 5export PrintTrees; 6 7-- export LabelledEventTrees; 8-- export florets; 9-- export PrintTree; 10-- export PrintTrees; 11-- export TreeToPoly; 12-- export AreGoodFlorets; 13-- export OneTree; 14 15---------------------------------------------------------------------- 16 17-- define OneTree(c_T) 18-- return OneTreeRec(support(c_T)); 19-- enddefine; -- OneTree 20 21-- define OneTreeRec(mons) 22-- PD := PrimaryDecomposition(ideal(mons)); 23-- if len(PD) = 0 then return []; endif; 24-- if len(PD) = 1 then return gens(PD[1]); endif; 25-- Root := gens(PD[1]); -- finds the shortest one 26-- foreach C in PD do 27-- if len(gens(C))<len(Root) then Comp:=gens(Root); endif; 28-- endforeach; 29-- L := [ [X, OneTreeRec(diff([T/X | T in mons and IsDivisible(T,X)],[X]))] 30-- | X In Root ]; 31-- return L; 32-- enddefine; -- OneTree 33 34---------------------------------------------------------------------- 35 36define StagedTrees(c_T) 37 AllT := LabelledEventTrees(c_T); 38 return [T in AllT | AreGoodFlorets(florets(T))]; 39enddefine; -- StagedTrees 40 41 42define LabelledEventTrees(c_T) 43 S := support(c_T); 44 foreach pp in S do 45 if max(exponents(pp))>1 then error("not square-free"); endif; 46 endforeach; 47 try 48 return AllTreesRec(support(c_T)); 49 UponError E do /*nothing*/; 50 endtry; 51 return []; 52enddefine; -- AllTrees 53 54 55define AllTreesRec(mons) 56 if mons=[1] then return [[]]; endif; 57 AllTrees := []; 58 foreach X in PrimaryDecomposition(ideal(mons)) do 59 if NumGens(X)=1 then error("Single branch"); endif; 60 try 61 T := AllTreesFromX(gens(X), mons); 62 AllTrees := concat(AllTrees, T); 63 UponError E do 64 if GetErrMesg(E) isin ["Single branch", 65 "More than one divisor", 66 "No subtrees"] then 67 continue; 68 endif; 69 error("unexpected error: "+GetErrMesg(E)); 70 endtry; 71 endforeach; 72 if AllTrees = [] then error("No subtrees"); endif; 73 return AllTrees; 74enddefine; -- AllTreesRec 75 76 77define AllTreesFromX(X, mons) 78 r := len(X); 79 Ix := [ [ pp in mons | IsDivisible(pp,X[i])] | i in 1..r ]; 80 if not(ArePairwiseDisjoint(Ix)) then error("More than one divisor"); endif; 81 T := [ AllTreesRec(Ix[i]/X[i]) | i in 1..r ]; ---> might give error 82 xT := [ [ [X[i],tree] | tree in T[i] ] | i in 1..r ]; 83 return CartesianProductList(xT); 84enddefine; -- AllTreesFromX 85 86 87define ArePairwiseDisjoint(Ix) 88 LL := [ L in Ix | L<>[1] ]; 89 for i := 1 to len(LL) do 90 for j := i+1 to len(LL) do 91 if intersection(LL[i], LL[j]) <> [] then return false; endif; 92 endfor; 93 endfor; 94 return true; 95enddefine; -- ArePairwiseDisjoint 96 97 98---------------------------------------------------------------------- 99 100define florets(tree) 101 if tree = [] then return []; endif; 102 if type(tree[1])=RINGELEM then return [tree]; endif; 103 F := [[node[1] | node in tree]]; 104 RecFlorets := flatten([ florets(node[2]) | node in tree and len(node)=2], 1); 105 return MakeSet(concat(F, RecFlorets)); 106enddefine; -- florets 107 108 109define AreGoodFlorets(florets) 110 for i := 1 to len(florets) do 111 for j := i+1 to len(florets) do 112 if [ petal in florets[i] | petal isin florets[j] ]<>[] then 113 return false; 114 endif; 115 endfor; 116 endfor; 117 return true; 118enddefine; -- AreGoodFlorets 119 120---------------------------------------------------------------------- 121 122-- define PrintTrees(trees) 123-- println "-- number of trees = ", len(trees), " -----------------------"; 124-- for i := 1 to len(trees) do 125-- println "-------- tree ",i, " ----------------------------------"; 126-- println "-- florets: ", florets(trees[i]); 127-- PrintTreeRec(trees[i], ""); 128-- println; 129-- endfor; 130-- println "-------- end trees -------------------------------"; 131-- enddefine; -- PrintTrees 132 133 134-- Define PrintNodeAndSubtree(L, spaces) 135-- print spaces, "[", L[1], ", "; 136-- PrintTreeRec(last(L), spaces+" "); 137-- print "]"; 138-- EndDefine; -- PrintNodeAndSubtree 139 140 141-- Define PrintTreeRec(L, spaces) 142-- if L=[] then print L; 143-- else 144-- print "[\n"; 145-- for i := 1 To len(L)-1 do 146-- PrintNodeAndSubtree(L[i], spaces+" "); 147-- print ",\n"; 148-- EndFor; 149-- PrintNodeAndSubtree(last(L), spaces+" "); 150-- print " ]"; 151-- EndIf; 152-- Return; 153-- EndDefine; -- PrintTreeRec 154 155---------------------------------------------------------------------- 156define PrintTrees(trees) 157 println "-- number of trees = ", len(trees), " -----------------------"; 158 for i := 1 to len(trees) do 159 println "-------- tree ",i, " ----------------------------------"; 160 println "-- florets: ", florets(trees[i]); 161 PrintTree(trees[i], ""); 162 println; 163 endfor; 164 println "-------- end trees -------------------------------"; 165enddefine; -- PrintTrees 166 167 168define PrintNodeAndSubtree(L, spaces) 169 print spaces, L[1], " "; 170 PrintTree(last(L), spaces+" "); 171enddefine; -- PrintNodeAndSubtree 172 173 174define PrintTree(L, opt spaces) 175 if not(IsDefined(spaces)) then spaces := ""; endif; 176 if L=[] then print ""; 177 else 178 print "< \n"; 179 for i := 1 To len(L)-1 do 180 PrintNodeAndSubtree(L[i], spaces+" "); 181 print "\n"; 182 endfor; 183 PrintNodeAndSubtree(last(L), spaces+" "); 184 endif; 185 return; 186enddefine; -- PrintTree 187 188---------------------------------------------------------------------- 189 190define TreeToPoly(T) 191 if type(T[2])=RINGELEM then return sum(T); endif; 192 return sum([T[i,1]*TreeToPoly(T[i,2]) | i in 1..len(T)]); 193enddefine; -- TreeToPoly 194 195 196endpackage; 197---------------------------------------------------------------------- 198