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 License
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 This file contains a function to compute the coarsest polyhedral structure of a
27 tropical variety
28 */
29
30 #include "polymake/client.h"
31 #include "polymake/Matrix.h"
32 #include "polymake/Rational.h"
33 #include "polymake/Vector.h"
34 #include "polymake/IncidenceMatrix.h"
35 #include "polymake/linalg.h"
36 #include "polymake/tropical/thomog.h"
37 #include "polymake/polytope/convex_hull.h"
38 #include "polymake/vector"
39
40 namespace polymake { namespace tropical {
41
42 /**
43 * @brief Takes a polyhedron in its V-description and tests whether it fulfills a given inequality
44 * @param Matrix<Rational> rays The rays/vertices of the polyhedron
45 * @param Matrix<Rational> linspace The lineality space generators of the polyhedron
46 * @param Vector<Rational> facet The inequality, in standard polymake format. Note that all three should be given in the same format, i.e. in homog. or non-homog. coordinates.
47 * @return True, if the polyhedron fulfills the inequality; false otherwise.
48 */
coneInHalfspace(const Matrix<Rational> & rays,const Matrix<Rational> & linspace,const Vector<Rational> & facet)49 bool coneInHalfspace(const Matrix<Rational>& rays, const Matrix<Rational>& linspace,
50 const Vector<Rational>& facet)
51 {
52 Matrix<Rational> allgenerators = rays / linspace;
53 Vector<Rational> product = allgenerators * facet;
54 for (Int i = 0; i < product.dim(); ++i) {
55 if (product[i] < 0) return false;
56 }
57 return true;
58 }
59
60 // Documentation see perl wrapper
61 template <typename Addition>
coarsen(BigObject complex,bool testFan=false)62 BigObject coarsen(BigObject complex, bool testFan = false)
63 {
64 // Extract values
65 Matrix<Rational> rays = complex.give("VERTICES");
66 rays = tdehomog(rays);
67 Matrix<Rational> linspace = complex.give("LINEALITY_SPACE");
68 linspace = tdehomog(linspace);
69 Int cmplx_dim = complex.give("PROJECTIVE_DIM");
70 Int cmplx_ambient_dim = complex.give("PROJECTIVE_AMBIENT_DIM");
71 bool weights_exist = false;
72 Vector<Integer> weights;
73 if (complex.lookup("WEIGHTS") >> weights) {
74 weights_exist = true;
75 }
76 IncidenceMatrix<> maximalCones = complex.give("MAXIMAL_POLYTOPES");
77 IncidenceMatrix<> codimOneCones = complex.give("CODIMENSION_ONE_POLYTOPES");
78 IncidenceMatrix<> codimInMaximal = complex.give("MAXIMAL_AT_CODIM_ONE");
79 IncidenceMatrix<> maximalOverCodim = T(codimInMaximal);
80 Int noOfCones = maximalCones.rows();
81
82 // For fan testing, we need the equations of all non-twovalent codimension one faces
83 Vector<Matrix<Rational>> codimOneEquations(codimOneCones.rows());
84 Set<Int> highervalentCodim;
85 for (Int cd = 0; cd < codimOneCones.rows(); ++cd) {
86 if (codimInMaximal.row(cd).size() > 2) {
87 codimOneEquations[cd] = null_space( rays.minor(codimOneCones.row(cd),All) / linspace);
88 highervalentCodim += cd;
89 } //END only for >2-valent
90 } //END compute codim one cone equations.
91
92 // If the fan has no rays, it has only lineality space and we return the complex itself
93 if (rays.rows() == 0) return complex;
94
95 // If the fan has no codim 1 faces, it must be 0- dimensional or a lineality space and we return
96 // the complex itself
97 if (codimOneCones.rows() == 0) return complex;
98
99 // Compute equivalence classes of maximal cones
100 std::vector<Set<Int>> equivalenceClasses;
101 std::vector<bool> hasBeenAdded(noOfCones); //contains whether a cone has been added to an equivalence class
102
103 for (Int mc = 0; mc < noOfCones; ++mc) {
104 if (!hasBeenAdded[mc]) {
105 equivalenceClasses.push_back(Set<Int>{ mc });
106 // Do a breadth-first search for all other cones in the component
107 std::list<Int> queue;
108 queue.push_back(mc);
109 // Semantics: Elements in that queue have been added but their neighbours might not
110 while (!queue.empty()) {
111 // Take the first element and find its neighbours
112 Int node = queue.front();
113 queue.pop_front();
114 Set<Int> cdset = maximalOverCodim.row(node);
115 for (auto cd = entire(cdset); !cd.at_end(); ++cd) {
116 Set<Int> otherMaximals = codimInMaximal.row(*cd) - node;
117 // We are only interested in codim-one-faces that have exactly two adjacent maximal cones
118 if (otherMaximals.size() == 1) {
119 Int othermc = *(otherMaximals.begin());
120 if (!hasBeenAdded[othermc]) {
121 // Now add the cone to the equivalence class of mc
122 equivalenceClasses.back() += othermc;
123 queue.push_back(othermc);
124 hasBeenAdded[othermc] = true;
125 }
126 }
127 }
128 } //End iterate queue
129 }
130 } //END iterate maximal cones
131
132 /* Test equivalence classes for convexit:
133 An equivalence class S = (s1,...,sr) is convex, iff it fulfills the following criterion: Let T be the set of all outer (i.e. not two-valent) codim one facets and write f_t >= a_t for the corresponding facet inequality. Then for each pair t,t' in T the polyhedron t must fulfill the inequality f_t' >= a_t'.
134 */
135 if (testFan) {
136 // Check each equivalence class
137 for (const auto& conesInClass : equivalenceClasses) {
138 // First, we compute the equation of the subdivision class
139 Matrix<Rational> classEquation =
140 null_space(rays.minor(maximalCones.row(*(conesInClass.begin())), All) / linspace);
141 // Compute all outer codim one faces and their facet inequality;
142 Set<Int> outerFacets;
143 Map<Int, Vector<Rational>> facetInequality;
144 for (auto mc = entire(conesInClass); !mc.at_end(); ++mc) {
145 Set<Int> outerCodimInMc = maximalOverCodim.row(*mc) * highervalentCodim;
146 outerFacets += outerCodimInMc;
147 for (auto fct = entire(outerCodimInMc); !fct.at_end(); ++fct) {
148 // Find the one equation of fct that is not an equation of the class
149 for (Int r = 0; r < codimOneEquations[*fct].rows(); ++r) {
150 if (rank(classEquation / codimOneEquations[*fct].row(r)) > cmplx_ambient_dim - cmplx_dim) {
151 // Find the right sign and save the inequality
152 Int additionalRay = *( (maximalCones.row(*mc) - codimOneCones.row(*fct)).begin());
153 facetInequality[*fct] = codimOneEquations[*fct].row(r);
154 if (facetInequality[*fct] * rays.row(additionalRay) < 0) {
155 facetInequality[*fct] = - facetInequality[*fct];
156 }
157 break;
158 }
159 } //END iterate equations of fct
160 } //END iterate all outer facets of the cone
161 } //END iterate all maximal cones in the class
162
163 // Now go through all pairs of outer facets
164 for (const Int out1 : outerFacets) {
165 for (const Int out2 : outerFacets - out1) {
166 // Check if out2 fulfills the facet inequality of out1
167 if (!coneInHalfspace(rays.minor(codimOneCones.row(out2),All), linspace, facetInequality[out1])) {
168 throw std::runtime_error("The equivalence classes are not convex. There is no coarsest structure.");
169 }
170 } //END iterate second facet.
171 } //END iterate first facet
172 } //END iterate classes
173 } //END test convexity of equivalence classes
174
175 // Now compute the new cones as unions of the cones in each equivalence class
176 Matrix<Rational> newlin;
177 Int newlindim = 0;
178 bool newlin_computed = false;
179 Vector<Set<Int>> newcones;
180 Vector<Integer> newweights;
181 Matrix<Rational> complete_matrix = rays / linspace;
182 Set<Int> used_rays;
183
184 for (const auto& conesInClass : equivalenceClasses) {
185 Matrix<Rational> union_rays(0, rays.cols());
186 Vector<Int> union_ray_list;
187 for (auto mc = entire(conesInClass); !mc.at_end(); ++mc) {
188 union_rays /= rays.minor(maximalCones.row(*mc),All);
189 union_ray_list |= Vector<Int>(maximalCones.row(*mc));
190 }
191
192 const auto union_cone = polytope::get_non_redundant_points(union_rays, linspace, true);
193
194 // If we need to test fan-ness, we need to check if different classes have different lin. spaces
195 if (testFan && newlin_computed) {
196 Matrix<Rational> otherlin = (union_rays / linspace).minor(union_cone.second, All);
197 if (newlindim != rank(otherlin)) {
198 throw std::runtime_error("Different classes have different lineality spaces. There is no coarsest structure.");
199 } else if (rank(newlin / otherlin) > rank(newlin)) {
200 throw std::runtime_error("Different classes have different lineality spaces. There is no coarsest structure.");
201 }
202 } //END check equality of lineality spaces.
203
204 if (!newlin_computed) {
205 newlin = (union_rays / linspace).minor(union_cone.second, All);
206 if (testFan) newlindim = rank(newlin);
207 newlin_computed = true;
208 }
209
210 // Convert indices of rays in union_rays to indices in rays
211 Set<Int> ray_set(union_ray_list.slice(union_cone.first));
212
213 newcones |= ray_set;
214 used_rays += ray_set;
215 if (weights_exist) newweights |= weights[*(conesInClass.begin())];
216 }
217
218 // Some rays might become equal (modulo lineality space) when coarsening, so we have to clean up
219 Map<Int, Int> ray_index_conversion;
220 Int next_index = 0;
221 Matrix<Rational> final_rays(0, complete_matrix.cols());
222 Int linrank = rank(newlin);
223
224 for (const Int r1 : used_rays) {
225 if (!keys(ray_index_conversion).contains(r1)) {
226 ray_index_conversion[r1] = next_index;
227 final_rays /= complete_matrix.row(r1);
228
229 for (const Int r2 : used_rays) {
230 if (!keys(ray_index_conversion).contains(r2)) {
231 // Check if both rays are equal mod lineality space
232 Vector<Rational> diff_vector = complete_matrix.row(r1) - complete_matrix.row(r2);
233 if (rank(newlin / diff_vector) == linrank) {
234 ray_index_conversion[r2] = next_index;
235 }
236 } //END if second not mapped yet
237 } //END iterate second ray index
238 ++next_index;
239 } //END if first not mapped yet
240 } //END iterate first ray index
241
242 // Now convert the cones
243 for (Int nc = 0; nc < newcones.dim(); ++nc) {
244 newcones[nc] = Set<Int>{ ray_index_conversion.map(newcones[nc]) };
245 }
246
247 // Produce final result
248 BigObject result("Cycle", mlist<Addition>());
249 result.take("VERTICES") << thomog(final_rays);
250 result.take("MAXIMAL_POLYTOPES") << newcones;
251 result.take("LINEALITY_SPACE") << thomog(newlin);
252 if (weights_exist) {
253 result.take("WEIGHTS") << newweights;
254 }
255
256 return result;
257 }
258
259 UserFunctionTemplate4perl("# @category Basic polyhedral operations"
260 "# Takes a tropical variety on which a coarsest polyhedral structure exists"
261 "# and computes this structure."
262 "# @param Cycle<Addition> complex A tropical variety which has a unique "
263 "# coarsest polyhedral structure "
264 "# @param Bool testFan (Optional, FALSE by default). Whether the algorithm should perform some consistency "
265 "# checks on the result. If true, it will check the following: "
266 "# - That equivalence classes of cones have convex support"
267 "# - That all equivalence classes have the same lineality space"
268 "# If any condition is violated, the algorithm throws an exception"
269 "# Note that it does not check whether equivalence classes form a fan"
270 "# This can be done via [[fan::check_fan]] afterwards, but it is potentially slow."
271 "# @return Cycle<Addition> The corresponding coarse complex. Throws an "
272 "# exception if testFan = True and consistency checks fail.",
273 "coarsen<Addition>(Cycle<Addition>; $=0)");
274 } }
275