1 2 /* 3 This program is free software; you can redistribute it and/or 4 modify it under the terms of the GNU General Public License 5 as published by the Free Software Foundation; either version 2 6 of the License, or (at your option) any later version. 7 8 This program is distributed in the hope that it will be useful, 9 but WITHOUT ANY WARRANTY; without even the implied warranty of 10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 11 GNU General Public License for more details. 12 13 You should have received a copy of the GNU General Public License 14 along with this program; if not, write to the Free Software 15 Foundation, Inc., 51 Franklin Street, Fifth Floor, 16 Boston, MA 02110-1301, USA. 17 18 --- 19 Copyright (C) 2011 - 2015, Simon Hampe <simon.hampe@googlemail.com> 20 21 --- 22 Copyright (c) 2016-2021 23 Ewgenij Gawrilow, Michael Joswig, and the polymake team 24 Technische Universität Berlin, Germany 25 https://polymake.org 26 27 Computes the minimal interior cells of a subdivision of a polyhedron 28 */ 29 30 #include "polymake/client.h" 31 #include "polymake/Set.h" 32 #include "polymake/Array.h" 33 #include "polymake/Matrix.h" 34 #include "polymake/SparseMatrix.h" 35 #include "polymake/IncidenceMatrix.h" 36 #include "polymake/Vector.h" 37 #include "polymake/Rational.h" 38 #include "polymake/Map.h" 39 #include "polymake/linalg.h" 40 #include "polymake/integer_linalg.h" 41 #include "polymake/tropical/thomog.h" 42 #include "polymake/tropical/misc_tools.h" 43 #include "polymake/tropical/lattice.h" 44 #include "polymake/tropical/linear_algebra_tools.h" 45 #include "polymake/polytope/convex_hull.h" 46 47 48 namespace polymake { namespace tropical { 49 50 /* 51 * @brief Computes the set of indices such that the corresponding entries 52 * in a given vector are either all zero or nonzero. 53 * @param Vector v The vector 54 * @param bool search_zeros: If true, find all zero entries, if false all non-zero 55 * entries 56 * @return Set<Int> 57 */ 58 template <typename T> 59 Set<Int> binaryFinder(const GenericVector<T>& v, bool search_zeros) 60 { 61 if (search_zeros) 62 return indices(attach_selector(v.top(), operations::is_zero())); 63 else 64 return indices(attach_selector(v.top(), operations::non_zero())); 65 } 66 67 /** 68 * @brief Finds all the maximal sets and returns the row indices 69 */ 70 Set<Int> find_maximal_faces(const IncidenceMatrix<>& mat) 71 { 72 Set<Int> smaller_sets; 73 for (Int i = 0; i < mat.rows(); ++i) { 74 for (Int j = 0; j < mat.rows(); ++j) { 75 if (i != j) { 76 if (mat.row(i) * mat.row(j) == mat.row(i)) smaller_sets += i; 77 } 78 } 79 } 80 return sequence(0,mat.rows()) - smaller_sets; 81 } 82 83 inline Set<Int> pairset(Int i, Int j) 84 { 85 return Set<Int>{ i, j }; 86 } 87 88 IncidenceMatrix<> minimal_interior(const Matrix<Rational>& vertices, const IncidenceMatrix<>& polytopes) 89 { 90 // If it only has one cone, it is minimal 91 if (polytopes.rows() == 1) { 92 return polytopes; 93 } 94 95 IncidenceMatrix<> result(0, vertices.rows()); 96 97 // Compute facet normals of the original polyhedron 98 const Matrix<Rational> facet_normals = polytope::enumerate_facets(vertices, false).first; 99 100 // The first step is computing interior and exterior codimension one cells 101 102 // Exterior one are the maximal intersections of the facet normals with the polytopes 103 // FIXME: growing IncidenceMatrix 104 IncidenceMatrix<> intersect_fn_with_polytopes(0, vertices.rows()); 105 for (Int f = 0; f < facet_normals.rows(); f++) { 106 // Find the vertices of the polytopes lying in the facet 107 Set<Int> facet_vertices = binaryFinder(vertices * facet_normals.row(f), true); 108 for (Int p = 0; p < polytopes.rows(); ++p) { 109 intersect_fn_with_polytopes /= facet_vertices * polytopes.row(p); 110 } 111 } 112 113 IncidenceMatrix<> exterior_codim = 114 intersect_fn_with_polytopes.minor(find_maximal_faces(intersect_fn_with_polytopes), All); 115 116 // Interior ones are the maximal intersections of the polytopes 117 IncidenceMatrix<> pairwise_intersections(0,vertices.rows()); 118 IncidenceMatrix<> associated_pairs(0,polytopes.rows()); 119 for (Int p1 = 0; p1 < polytopes.rows(); ++p1) { 120 for (Int p2 = p1+1; p2 < polytopes.rows(); ++p2) { 121 pairwise_intersections /= (polytopes.row(p1)*polytopes.row(p2)); 122 associated_pairs /= pairset(p1,p2); 123 } 124 } 125 Set<Int> pairwise_maximal = find_maximal_faces(pairwise_intersections); 126 IncidenceMatrix<> interior_codim = pairwise_intersections.minor(pairwise_maximal,All); 127 // Tells for all polytopes, which interior codim 1 cells they contain. 128 IncidenceMatrix<> interior_in_max = T(associated_pairs.minor(pairwise_maximal,All)); 129 130 // If there are no codimension one cones, return the maximal cells 131 if (exterior_codim.rows() + interior_codim.rows() == 0) { 132 return polytopes; 133 } 134 135 // For each maximal cone, we compute all its minimal interior faces as maximal intersections 136 // of interior codimension one faces. However, we only use codim-1-faces, that we haven't used 137 // in another maximal cone so far (Since any minimal face in such a face is hence also a minimal 138 // face of this other maximal cone, so we already computed it). 139 140 Set<Int> markedFaces; 141 142 for (Int mc = 0; mc < polytopes.rows(); ++mc) { 143 // Compute all non-marked codim-1-cells of mc. If there are none left, go to the next cone 144 Vector<Int> nonmarked(interior_in_max.row(mc) - markedFaces); 145 if (nonmarked.dim() == 0) continue; 146 Int k = nonmarked.dim(); 147 // ordered list of indices of interior codim-1-cells (in nonmarked) 148 // indices != -1 correspond to codim-1-cells that we intersect to obtain a minimal face 149 Vector<Int> currentSet(k,-1); 150 currentSet[0] = 0; 151 // Indicates the index below the next cone index (in nonmarked) we should try to add. Is always 152 // larger equal than the last element != -1 in currentSet 153 Int lowerBound = 0; 154 // Indicates the current position in currentSet we're trying to fill 155 Int currentPosition = 1; 156 // Indicates at position i < currentPosition the intersection of the cones specified by the 157 // elements currentSet[0] .. currentSet[i] 158 Vector<Set<Int>> currentIntersections(k); 159 currentIntersections[0] = interior_codim.row(nonmarked[0]); 160 // Now iterate all intersections in a backtrack algorithm: 161 // If an intersection is maximal, we don't need to go any further 162 // We stop when we have tried all possibilities at the first position 163 while (!(currentPosition == 0 && lowerBound == k-1)) { 164 // Try the next posssible index 165 Int j = lowerBound+1; 166 // If we're already beyond k-1, we have found a maximal intersection 167 // Check if it is a minimal face, then go back one step 168 if (j == k) { 169 // We test whether the set is not contained in any border face and does not contain 170 // any existing minimal face 171 Set<Int> potentialMinimal(currentIntersections[currentPosition-1]); 172 bool invalid = false; 173 for (Int ec = 0; ec < exterior_codim.rows(); ++ec) { 174 if (potentialMinimal.size() == (potentialMinimal * exterior_codim.row(ec)).size()) { 175 invalid = true; 176 break; 177 } 178 } 179 for (Int mf = 0; mf < result.rows() && !invalid; ++mf) { 180 if (result.row(mf).size() == (result.row(mf) * potentialMinimal).size()) { 181 invalid = true; 182 break; 183 } 184 } 185 if (!invalid) 186 result /= potentialMinimal; 187 // Go back one step 188 lowerBound = currentSet[currentPosition-1]; 189 --currentPosition; 190 } else { 191 // Compute the intersection with the next codim 1 cell 192 // (If we're at the first position, we just insert a cell) 193 Set<Int> intersection; 194 if (currentPosition == 0) 195 intersection = interior_codim.row(nonmarked[j]); 196 else 197 intersection = currentIntersections[currentPosition-1] * interior_codim.row(nonmarked[j]); 198 // If its still not empty, go forward one step 199 if (!intersection.empty()) { 200 currentSet[currentPosition] = j; 201 currentIntersections[currentPosition] = intersection; 202 ++currentPosition; 203 lowerBound = j; 204 } else { 205 // Otherwise go up one step 206 ++lowerBound; 207 } 208 } //END if j == k 209 } //END while(...) 210 211 // Now mark all codim -1-faces of mc 212 markedFaces += interior_in_max.row(mc); 213 214 } //END iterate maximal cones 215 216 return result; 217 } //END minimal_interior 218 219 220 IncidenceMatrix<> refined_local_cones(BigObject localized_cycle, BigObject refining_cycle) 221 { 222 IncidenceMatrix<> local_restriction = localized_cycle.give("LOCAL_RESTRICTION"); 223 Matrix<Rational> local_vertices = localized_cycle.give("VERTICES"); 224 local_vertices = tdehomog(local_vertices); 225 Matrix<Rational> local_lineality = localized_cycle.give("LINEALITY_SPACE"); 226 local_lineality = tdehomog(local_lineality); 227 Int local_lineality_dim = localized_cycle.give("LINEALITY_DIM"); 228 229 Matrix<Rational> refining_vertices = refining_cycle.give("VERTICES"); 230 refining_vertices = tdehomog(refining_vertices); 231 IncidenceMatrix<> refining_polytopes = refining_cycle.give("MAXIMAL_POLYTOPES"); 232 Int refining_lineality_dim = refining_cycle.give("LINEALITY_DIM"); 233 234 Set<Set<Int>> result; 235 236 // Iterate all local cones of the localized cycle 237 for (Int lc = 0; lc < local_restriction.rows(); ++lc) { 238 // Go through maximal cones of refining complex and check which ones intersect in the right dimension 239 Matrix<Rational> lc_vertices = local_vertices.minor(local_restriction.row(lc),All); 240 Int local_dim = rank(lc_vertices) + local_lineality_dim; 241 RestrictedIncidenceMatrix<> lc_refiners; 242 for (Int mc = 0; mc < refining_polytopes.rows(); ++mc) { 243 Set<Int> contained_vertices; 244 Set<Int> mcset = refining_polytopes.row(mc); 245 for (auto vert = entire(mcset); !vert.at_end(); ++vert) { 246 if (is_ray_in_cone(lc_vertices, local_lineality, refining_vertices.row(*vert), false)) 247 contained_vertices += (*vert); 248 249 if (rank( refining_vertices.minor(contained_vertices,All) ) + refining_lineality_dim == local_dim) 250 lc_refiners /= contained_vertices; 251 } 252 } //END iterate maximal cones 253 254 // Now that we have the subdividing cones, compute the minimal interior ones 255 IncidenceMatrix<> mi_res = minimal_interior(refining_vertices, IncidenceMatrix<>(std::move(lc_refiners))); 256 for (auto mi_cell = entire(rows(mi_res)); !mi_cell.at_end(); mi_cell++) 257 result += *mi_cell; 258 259 } //END iterate local cones 260 return IncidenceMatrix<>(result); 261 262 } //END refined_local_cones 263 264 } } 265