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