1 /*
2 This program is free software; you can redistribute it and/or
3 modify it under the terms of the GNU General Public License
4 as published by the Free Software Foundation; either version 2
5 of the License, or (at your option) any later version.
6
7 This program is distributed in the hope that it will be useful,
8 but WITHOUT ANY WARRANTY; without even the implied warranty of
9 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10 GNU General Public License for more details.
11
12 You should have received a copy of the GNU General Public LicenseRestrictedIncidenceMatrix
13 along with this program; if not, write to the Free Software
14 Foundation, Inc., 51 Franklin Street, Fifth Floor,
15 Boston, MA 02110-1301, USA.
16
17 ---
18 Copyright (C) 2011 - 2015, Simon Hampe <simon.hampe@googlemail.com>
19
20 ---
21 Copyright (c) 2016-2021
22 Ewgenij Gawrilow, Michael Joswig, and the polymake team
23 Technische Universität Berlin, Germany
24 https://polymake.org
25
26 Computes an augmented bergman matroid fan by extending the chains of the
27 lattice of flats with their respective compatible independent sets.
28 */
29
30
31 #include "polymake/client.h"
32 #include "polymake/Matrix.h"
33 #include "polymake/IncidenceMatrix.h"
34 #include "polymake/Rational.h"
35 #include "polymake/Vector.h"
36 #include "polymake/Set.h"
37 #include "polymake/graph/Lattice.h"
38 #include "polymake/graph/Decoration.h"
39 #include "polymake/graph/maximal_chains.h"
40 #include "polymake/tropical/specialcycles.h"
41
42 namespace polymake { namespace tropical {
43
44 using graph::Lattice;
45 using graph::lattice::Sequential;
46 using graph::lattice::BasicDecoration;
47
48 template <typename Addition>
augmented_matroid_fan(BigObject matroid)49 BigObject augmented_matroid_fan(BigObject matroid)
50 {
51 const Int n = matroid.give("N_ELEMENTS");
52 const Array<Set<Int>> bases = matroid.give("BASES");
53 const Int full_rank = matroid.give("RANK");
54 const Set<Int> loops = matroid.give("LOOPS");
55
56 const BigObject lattice_flats = matroid.give("LATTICE_OF_FLATS");
57 const int top_index = lattice_flats.give("TOP_NODE");
58 const NodeMap<Directed,Set<Int>> flats = lattice_flats.give("FACES");
59
60 const bool ascending = (top_index != 0); // Is FACES ordered rank-ascending, or rank-descending?
61
62 if (loops.size() != 0) {
63 throw std::runtime_error("augmented_matroid_fan: The passed matroid is not loopless!");
64 }
65
66 // Create the rays.
67
68 // We have two homogenizing coordinates (in the columns).
69 // First 0 says "this is a ray"
70 // Second 0 comes from homogenization onto the tropical torus
71 Matrix<Rational> rays(n + flats.size(), n + 2);
72
73 // First type of ray corresponds to elements of the matroid's ground set.
74 // These are the rays generated by all the standard basis vectors.
75 // Hence we start the identity matrix (with two homogenizing coordinates).
76 for (int i = 0; i < n; i++) {
77 rays(i, i + 2) = Addition::orientation();
78 }
79
80 // Second type of ray corresponds to nonempty proper flats.
81 // These rays are generated by the negative complement of the
82 // indicator vector.
83 // E.g. if the flat is {0, 2, 3} out of 5 Elements
84 // Its indicator vector is ( 1, 0, 1, 1, 0)
85 // And the form we require is ( 0, -1, 0, 0, -1)
86
87 for (int i = 0; i < flats.size(); i++) {
88 if (i == top_index) { // Skip the ground set, it's not a proper flat.
89 continue;
90 }
91
92 Set<Int> flat = flats[i];
93 for (int j = 0; j < n; j++) {
94 if(!flat.exists(j)) {
95 // We already wrote first n rows previously. Hence + n
96 // If we have a desecending order, we skip the first flat because its the ground set.
97 const Int row_index = ascending ? i + n : i + n - 1;
98
99 const Int col_index = j + 2; // + 2 due to the two homogenizing coords.
100 rays(row_index, col_index) = -1 * Addition::orientation();
101 }
102 }
103 }
104
105 // Lastly, we add an artificial "ray" that corresponds to the origin in a fan.
106 // This is required because a tropical cycle is in general a polyhedral complex.
107 rays(rays.rows() - 1, 0) = 1;
108
109
110 // Next we must compute the maximal cones of the fan.
111 // These maximal cones correspond to tuples (I, F), where
112 // I is an independent set in the matroid, and F is a maximal flag of proper
113 // flats with the condition that every proper flat in F must contain I.
114
115 // We generate all the maximal flags of flats without considering I
116 // We want the empty face in our flags, but not the ground set - false, true.
117 const Lattice<BasicDecoration, Sequential> dec_flats(lattice_flats);
118 const Array<Set<Int>> maximal_flags = maximal_chains(dec_flats, false, true);
119
120 // Output maximal cones of the fan.
121 RestrictedIncidenceMatrix<> maximal_cones(0);
122 Int maximal_cones_rows = 0;
123
124 // All the maximal flags of flats give us the maximal cones of the form ({}, F)
125 // These are the first maximal cones we create.
126 for (int i = 0; i < maximal_flags.size(); i++) {
127 const Set<Int> flag = maximal_flags[i];
128 Set<Int> new_row;
129 for (Int face_index : flag) {
130 // + n Since we wish to adress the rays which correspond to faces.
131 // - 1 for descending order, since our 0 corresponds to the "first"
132 // face which is not the ground set, but 0 corresponds to the
133 // ground set in a descending order.
134
135 Int new_elem = ascending ? face_index + n : face_index + n - 1;
136 new_row += new_elem;
137 }
138 maximal_cones /= new_row;
139 maximal_cones_rows++;
140 }
141
142
143 // We now iteratively create the maximal cones for the intermediate ranks.
144 // We iterate on the rank 0 < r < full_rank of the independent set I in (I, F).
145
146 // Assume we are in the situtation that we wish to construct all (I, F)
147 // with rank(I) = r, and we know how all the (I', F') look like with rank(I') < r.
148
149 // Take a (I', F'), with rank(I') = r - 1, we write
150 // F' = { F_(r-1), F_(r), F_(r+1), ..., F_(full_rank - 1) }
151 // with F_(r-1) ⊂ F_(r) ⊂ F_(r+1) ⊂ ... ⊂ F_(full_rank - 1)
152
153 // Then these (I', F') generate new (I, F) with rank(I) = r.
154 // We take F = F' - F_(r-1).
155 // And I = I' + e, for any e in F_(r) \ F_(r-1).
156
157 // Notice that we get all possible pairs (I, F) in this way if we take all
158 // previous (I', F'). Let I = (x_1, ..., x_r). Then we may arrive at
159 // (I,F) by considering the sequence of flats cl(x_1), cl(x_1, x_2), ..., cl(x_1, x_2, x_r)
160 // By permuting these x_1, ..., x_r we arrive at dupliacates of (I,F).
161 // In fact, these are exactly all the duplicates.
162 // So, in order to prevent duplicates, we demand that x_1 > x_2 > ... > x_r. (viewed by their indices)
163 // Exactly one permutation will satisify this property.
164
165 // These two variables indicate the rows of the matrix we have
166 // written in the previous outermost iteration.
167 // I.e. they correspond to the (I', F') from which we will generate
168 // the next iteration's (I, F).
169 Int last_rows_from = 0;
170 Int last_rows_to = maximal_cones_rows; // Notice we begin with the already generated rank 0 cones.
171
172 for (int rank = 1; rank < full_rank; rank++) {
173 // Iterate through all (I', F')
174 for (Int row_index = last_rows_from; row_index < last_rows_to; row_index++) {
175
176 const Set<Int> the_row = maximal_cones.row(row_index);
177 const Set<Int> i_prime = the_row * range(0, n - 1); // Corresponds to I'.
178 const Set<Int> f_prime = the_row - i_prime; // Corresponds to F' (in indices).
179
180 Set<Int> f = f_prime; // Will correspond to F after if branch.
181 Set<Int> smallest_face; // Will correspond to F_(r-1) after if branch.
182 Set<Int> second_smallest_face; // Will correspond to F_(r) after if branch.
183 if (ascending) { // Smallest face is indexed by smallest integer.
184 Int index = f.front(); f.pop_front();
185 // Need to correct index.
186 // Must subtract by n, since the first n elements correspond to matroid elements in the_row.
187 index -= n;
188 smallest_face = flats[index];
189 index = f.front();
190 index -= n; // Ditto with the correction.
191 second_smallest_face = flats[index];
192 } else { // Smallest face is indexed by largest integer.
193 Int index = f.back(); f.pop_back();
194 // Need to correct the index.
195 // Must subtract by n, since the first n elements correspond to matroid elements in the_row.
196 // 0 refers to ground set in flats;
197 // But 0 refers to "first" set which is not ground set in maximal_cones after - n.
198 // So we add 1 back.
199 index = index - n + 1;
200 smallest_face = flats[index];
201 index = f.back();
202 index = index - n + 1; // Ditto with the correction.
203 second_smallest_face = flats[index];
204 }
205
206 // We now have constructed F, F_(r) and F_(r-1).
207 // We can now easily itearte over all e with e in F_(r) \ F_(r-1)
208 // These e allow us to construct all possible I = I' + e.
209 // This is all we need to generate another row.
210 Set<Int> all_es = second_smallest_face - smallest_face;
211 // The smallest value of I'. If I' is empty, then effectively no min value.
212 const Int min_value = i_prime.size() == 0 ? n + 1 : i_prime.front();
213
214 for (const Int e : all_es) {
215 if (e > min_value) {
216 // To create (I, F) with I = (x_1, ..., x_r)
217 // Where x_1 was added first, then x_2, and so on.
218 // We only allow the creation via x_1 > ... > x_r.
219 // This prevents duplicates.
220 break;
221 }
222
223 Set<Int> new_row = i_prime + scalar2set(e) + f;
224 maximal_cones /= new_row;
225 maximal_cones_rows++;
226 }
227 }
228
229 last_rows_from = last_rows_to;
230 last_rows_to = maximal_cones_rows;
231 }
232
233 // Now add the maximal cones which correspond to bases
234 // I.e. they are of the form (B, {}) where B is a basis.
235 IncidenceMatrix<> bases_matrix(bases);
236 maximal_cones /= bases_matrix;
237 maximal_cones_rows += bases_matrix.rows();
238
239 // Since a cycle is in general a polyhedral complex we add
240 // the zero as a vertex of all the maximal cones.
241 for (int i = 0; i < maximal_cones_rows; i++) {
242 maximal_cones(i, rays.rows() - 1) = 1;
243 }
244
245 return BigObject("Cycle", mlist<Min>(),
246 "PROJECTIVE_VERTICES", rays,
247 "MAXIMAL_POLYTOPES", IncidenceMatrix<>(std::move(maximal_cones)),
248 "WEIGHTS", ones_vector<Integer>(maximal_cones_rows));
249 }
250
251
252 UserFunctionTemplate4perl("# @category Matroids"
253 "# Computes the augmented Bergman fan of a matroid."
254 "# Note that this is potentially very slow for large matroids."
255 "# A definition of the augmented Bergman fan can be found in arXiv:2002.03341. See also the follow up paper arXiv:2010.06088."
256 "# The algorithim used to construct the augemented Bergman fan closely follows its description in the first paper."
257 "# @param matroid::Matroid A matroid. Should be loopfree."
258 "# @tparam Addition Min or max, determines the matroid fan coordinates."
259 "# @example [application matroid]"
260 "# > $m = matroid::fano_matroid;"
261 "# > $f = tropical::augmented_matroid_fan<Min>($m);"
262 "# @return tropical::Cycle<Addition>",
263 "augmented_matroid_fan<Addition>(matroid::Matroid)");
264
265 } }
266