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 provides functionality to compute certain special tropical varieties
27 */
28
29 #pragma once
30
31 #include "polymake/client.h"
32 #include "polymake/PowerSet.h"
33 #include "polymake/Set.h"
34 #include "polymake/Matrix.h"
35 #include "polymake/Vector.h"
36 #include "polymake/Rational.h"
37 #include "polymake/Array.h"
38 #include "polymake/IncidenceMatrix.h"
39 #include "polymake/Graph.h"
40 #include "polymake/linalg.h"
41 #include "polymake/tropical/thomog.h"
42 #include "polymake/tropical/misc_tools.h"
43
44 namespace polymake { namespace tropical {
45
46 template <typename Addition>
empty_cycle(Int ambient_dim)47 BigObject empty_cycle(Int ambient_dim)
48 {
49 BigObject cycle("Cycle", mlist<Addition>(),
50 "VERTICES", Matrix<Rational>(0, ambient_dim+2),
51 "MAXIMAL_POLYTOPES", Array<Set<Int>>(),
52 "WEIGHTS", Vector<Integer>(),
53 "PROJECTIVE_AMBIENT_DIM", ambient_dim);
54 cycle.set_description() << "Empty cycle in dimension " << ambient_dim;
55 return cycle;
56 }
57
58 template <typename Addition>
point_collection(Matrix<Rational> m,const Vector<Integer> & weights)59 BigObject point_collection(Matrix<Rational> m, const Vector<Integer>& weights)
60 {
61 // Sanity check
62 if (m.rows() == 0)
63 throw std::runtime_error("No points given.");
64 if (m.rows() != weights.dim())
65 throw std::runtime_error("Number of points does not match number of weights");
66
67 // Create vertices
68 m = ones_vector<Rational>() | m;
69
70 // Create polytopes
71 Array<Set<Int>> polytopes(m.rows());
72 for (Int i = 0; i < polytopes.size(); ++i)
73 polytopes[i] = scalar2set(i);
74
75 return BigObject("Cycle", mlist<Addition>(),
76 "PROJECTIVE_VERTICES", m,
77 "MAXIMAL_POLYTOPES", polytopes,
78 "WEIGHTS", weights);
79 }
80
81 template <typename Addition>
82 BigObject uniform_linear_space(const Int n, const Int k, Integer weight = 1)
83 {
84 // Ensure that dimensions match
85 if (k > n)
86 throw std::runtime_error("Cannot create uniform linear space. Fan dimension is larger than ambient dimension.");
87 if (k < 0 || n < 0)
88 throw std::runtime_error("Cannot create uniform linear space. Negative dimension provided.");
89 if (k == 0) {
90 return point_collection<Addition>( Matrix<Rational>(1,n+1), ones_vector<Integer>(1));
91 }
92
93 // Create rays
94 Matrix<Rational> vertices(unit_matrix<Rational>(n+1));
95 vertices = zero_vector<Rational>(n+1) | vertices;
96 vertices *= Addition::orientation();
97 vertices = unit_vector<Rational>(n+2,0) / vertices;
98
99 // Create cones
100 Array<Set<Int>> polytopes{ all_subsets_of_k(sequence(1,n+1), k) };
101 for (Int i = 0; i < polytopes.size(); ++i)
102 polytopes[i] += 0;
103
104 // Create weights
105 Vector<Integer> weights = weight * ones_vector<Integer>(polytopes.size());
106
107 // Create final object
108 BigObject fan("Cycle", mlist<Addition>(),
109 "PROJECTIVE_VERTICES", vertices,
110 "MAXIMAL_POLYTOPES", polytopes,
111 "WEIGHTS", weights);
112 fan.set_description() << "Uniform linear space of dimension " << k << " in dimension " << n;
113 return fan;
114 }
115
116 template <typename Addition>
halfspace_subdivision(const Rational & a,const Vector<Rational> & g,const Integer & weight)117 BigObject halfspace_subdivision(const Rational& a, const Vector<Rational>& g, const Integer& weight)
118 {
119 // Sanity check
120 if (is_zero(g))
121 throw std::runtime_error("Zero vector does not define a hyperplane.");
122 if (!is_zero(accumulate(g, operations::add())))
123 throw std::runtime_error("Normal vector must be homogenous, i.e. sum of entries must be zero");
124
125 const Vector<Rational> apex = Rational(1) | (a / sqr(g))*g;
126 const Matrix<Rational> vertices = vector2row(apex) / (0 | g) / (0 | -g);
127 Matrix<Rational> lineality = zero_vector<Rational>() | null_space(g).minor(range_from(1), All);
128
129 const IncidenceMatrix<> polytopes({{0, 2}, {0, 1}});
130
131 return BigObject("Cycle", mlist<Addition>(),
132 "PROJECTIVE_VERTICES", vertices,
133 "MAXIMAL_POLYTOPES", polytopes,
134 "LINEALITY_SPACE", lineality,
135 "WEIGHTS", same_element_vector(weight, 2));
136 }
137
138 template <typename Addition>
projective_torus(Int n,Integer weight)139 BigObject projective_torus(Int n, Integer weight)
140 {
141 // Sanity check
142 if (n < 0) throw std::runtime_error("Negative ambient dimension is not allowed.");
143
144 const Matrix<Rational> vertex = vector2row(unit_vector<Rational>(n+2, 0));
145 const Matrix<Rational> lineality = zero_matrix<Rational>(n, 2) | unit_matrix<Rational>(n);
146
147 const IncidenceMatrix<> polytopes({{0}});
148
149 return BigObject("Cycle", mlist<Addition>(),
150 "PROJECTIVE_VERTICES", vertex,
151 "MAXIMAL_POLYTOPES", polytopes,
152 "LINEALITY_SPACE", lineality,
153 "WEIGHTS", same_element_vector(weight, 1));
154 }
155
156 template <typename Addition>
157 BigObject orthant_subdivision(Vector<Rational> point, Int chart = 0, Integer weight = 1)
158 {
159 if (point.dim() <= 2) {
160 throw std::runtime_error("Cannot create orthant subdivision. Vector dimension too small");
161 }
162
163 // Dehomogenize
164 point = tdehomog_vec(point,chart);
165 Int dim = point.dim() -1;
166 // Create ray matrix - first positive rays, then negative rays
167 Matrix<Rational> rays = unit_matrix<Rational>(dim);
168 rays /= (-unit_matrix<Rational>(dim));
169 // Prepend a zero and set the vertex as last ray
170 rays = zero_vector<Rational>() | rays;
171 rays /= point;
172
173 // Create cones
174 const auto seq = sequence(0,dim);
175 // All possible sign choices
176 RestrictedIncidenceMatrix<only_cols> cones_growing(0, rays.rows());
177 for (auto s = entire(all_subsets(seq)); !s.at_end(); ++s) {
178 Set<Int> complement = seq - *s;
179 Set<Int> rayset = *s;
180 // Add all rays from the current set with positive sign and all the others with negative sign
181 for (auto c = entire(complement); !c.at_end(); ++c) {
182 rayset += (*c + dim);
183 }
184 // Finally add the vertex
185 rayset += (rays.rows()-1);
186 cones_growing /= rayset;
187 } // END create cones
188
189 const IncidenceMatrix<> cones(std::move(cones_growing));
190
191 return BigObject("Cycle", mlist<Addition>(),
192 "PROJECTIVE_VERTICES", thomog(rays, chart),
193 "MAXIMAL_POLYTOPES", cones,
194 "WEIGHTS", same_element_vector(weight, cones.rows()));
195 }
196
197 template <typename Addition>
198 BigObject affine_linear_space(const Matrix<Rational> &generators, Vector<Rational> translate = Vector<Rational>(), Integer weight = 1)
199 {
200 // Sanity check
201 if (translate.dim() > 0 && translate.dim() != generators.cols()) {
202 throw std::runtime_error("affine_linear_space: Dimension mismatch.");
203 }
204 if (translate.dim() == 0)
205 translate = Vector<Rational>(generators.cols());
206
207 Matrix<Rational> vertices(1,generators.cols()+1);
208 vertices(0,0) = 1;
209 vertices.row(0).slice(range_from(1)) = translate;
210 Vector<Set<Int>> polytopes;
211 polytopes |= scalar2set(0);
212 Vector<Integer> weights(1);
213 weights[0] = weight;
214
215 return BigObject("Cycle", mlist<Addition>(),
216 "PROJECTIVE_VERTICES", vertices,
217 "MAXIMAL_POLYTOPES", polytopes,
218 "LINEALITY_SPACE", zero_vector<Rational>() | generators,
219 "WEIGHTS", weights);
220 }
221
222 ///////////////////////////////////////////////////////////////////////////////////////
223
224 template <typename Addition>
225 BigObject cross_variety(Int n, Int k, Rational h = 1, Integer weight = 1)
226 {
227 // Create the cube vertices
228 Matrix<Rational> rays = binaryMatrix(n);
229 Vector<Set<Int>> cones;
230
231 // Sanity check
232 if (n < k || k < 0 || h < 0) {
233 throw std::runtime_error("cross_variety: Invalid input parameters.");
234 }
235
236 // First we treat the special case of k = 0
237 if (k == 0) {
238 if (h == 0) {
239 rays = Matrix<Rational>(1, n+1);
240 rays(0, 0) = 1;
241 } else {
242 rays = ones_vector<Rational>() | (weight * rays);
243 }
244 const Int n_rays = rays.rows();
245 IncidenceMatrix<> polytopes(n_rays, n_rays);
246 for (Int r = 0; r < n_rays; ++r)
247 polytopes(r, r) = true;
248
249 return BigObject("Cycle", mlist<Addition>(),
250 "VERTICES", thomog(rays),
251 "MAXIMAL_POLYTOPES", polytopes,
252 "WEIGHTS", same_element_vector(weight, n_rays));
253 }
254
255 // Now create the k-skeleton of the n-cube: For each n-k-set S of 0,..,n-1 and for each vertex
256 // v of the n-k-dimensional cube: Insert the entries of v in S and then insert all possible
257 // vertices of the k-dimensional cube in S^c to obtain a k-dimensional face of the cube
258 Array<Set<Int>> nmkSets{ all_subsets_of_k(sequence(0,n), n-k) };
259 Matrix<Rational> nmkVertices = binaryMatrix(n-k);
260 Matrix<Rational> kVertices = binaryMatrix(k);
261
262 for (Int s = 0; s < nmkSets.size(); ++s) {
263 for (Int v = 0; v < nmkVertices.rows(); ++v) {
264 Set<Int> S = nmkSets[s];
265 Set<Int> newface;
266 Vector<Rational> vertex(n);
267 vertex.slice(S) = nmkVertices.row(v);
268 for (Int w = 0; w < kVertices.rows(); ++w) {
269 vertex.slice(~S) = kVertices.row(w);
270 newface += binaryIndex(vertex);
271 }
272 cones |= newface;
273 }
274 } //End create k-skeleton
275
276 Int vertexnumber = rays.rows();
277
278 // Now we also create the k-1-skeleton of the cube to compute the ray faces
279 Array<Set<Int>> nmlSets{ all_subsets_of_k(sequence(0,n), n-k+1) };
280 Matrix<Rational> nmlVertices = binaryMatrix(n-k+1);
281 Matrix<Rational> lVertices = binaryMatrix(k-1);
282 Vector<Set<Int>> raycones;
283
284 for (Int s = 0; s < nmlSets.size(); ++s) {
285 for (Int v = 0; v < nmlVertices.rows(); ++v) {
286 Set<Int> S = nmlSets[s];
287 Set<Int> newface;
288 Vector<Rational> vertex(n);
289 vertex.slice(S) = nmlVertices.row(v);
290 for (Int w = 0; w < lVertices.rows(); ++w) {
291 vertex.slice(~S) = lVertices.row(w);
292 newface += binaryIndex(vertex);
293 }
294 raycones |= newface;
295 }
296 } //End create k-1-skeleton
297
298 // We add a copy of each vertex and consider it a ray.
299 // Now, for each face S of the k-1-skeleton, we add a cone that contains for each i in S:
300 // The vertex i and its corresponding ray
301 if (k > 0) rays = rays / rays;
302 Int iter = raycones.size();
303 for (Int c = 0; c < iter; ++c) {
304 Set<Int> newface;
305 Set<Int> cubeface = raycones[c];
306 for (auto v = entire(cubeface); !v.at_end(); ++v) {
307 newface += *v;
308 Int rayindex = *v + vertexnumber;
309 newface += rayindex;
310 }
311 cones |= newface;
312 }
313
314 //In the degenerate case, where h = 0, we replace all nonfar vertices by a single one.
315 if (h == 0) {
316 rays = zero_vector<Rational>(rays.cols()) / rays.minor(~sequence(0,vertexnumber),All);
317 rays = unit_vector<Rational>(rays.rows(),0) | rays;
318 // Re-index cones and remove bounded ones
319 Set<Int> bad_indices = sequence(0,vertexnumber);
320 Set<Int> bounded_cones;
321 for (Int c = 0; c < cones.dim(); ++c) {
322 Set<Int> bad_of_c = cones[c] * bad_indices;
323 if (bad_of_c.size() == cones[c].size()) {
324 bounded_cones += c;
325 } else {
326 Set<Int> rest_of_c = cones[c] - bad_indices;
327 cones[c] = Set<Int>();
328 for (auto shiftindex = entire(rest_of_c); !shiftindex.at_end(); ++shiftindex) {
329 cones[c] += (*shiftindex - vertexnumber + 1);
330 }
331 cones[c] += 0;
332 }
333 }
334 cones = cones.slice(~bounded_cones);
335 }
336 // Otherwise scale vertices appropriately
337 else {
338 rays.minor(sequence(0, vertexnumber), All) *= h;
339 Vector<Rational> leading_coordinate = ones_vector<Rational>(vertexnumber) | zero_vector<Rational>(vertexnumber);
340 rays = leading_coordinate | rays;
341 }
342
343 return BigObject("Cycle", mlist<Addition>(),
344 "VERTICES", thomog(rays),
345 "MAXIMAL_POLYTOPES", cones,
346 "WEIGHTS", same_element_vector(weight, cones.dim()));
347 }
348
349 } }
350
351