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 Computes the moduli space of rational n-marked curves.
27 */
28
29
30 #include "polymake/client.h"
31 #include "polymake/Vector.h"
32 #include "polymake/Matrix.h"
33 #include "polymake/Set.h"
34 #include "polymake/Rational.h"
35 #include "polymake/Integer.h"
36 #include "polymake/PowerSet.h"
37 #include "polymake/Array.h"
38 #include "polymake/tropical/thomog.h"
39 #include "polymake/tropical/specialcycles.h"
40 //#include "polymake/atint/moduli.h"
41
42 namespace polymake { namespace tropical {
43
44
45 //Counts maximal cones by a simple formula
count_maximal_mn_cones(Int n)46 Integer count_maximal_mn_cones(Int n)
47 {
48 if (n == 3) {
49 return 1;
50 }
51 Integer result = 1;
52 Integer nint(n);
53 for (Int i = 0; i <= n-4; ++i) {
54 result *= 2*(nint-i) - 5;
55 }
56 return result;
57 }
58
59
60 // Documentation see perl wrapper
count_mn_cones(Int n,Int k)61 Integer count_mn_cones(Int n, Int k)
62 {
63 if (n == 3) {
64 return Integer(1);
65 }
66 if (k == n-3)
67 return count_maximal_mn_cones(n);
68
69 Int vertex_count = k+1;
70 Int seq_length = n + k-1;
71
72 // We compute the number of ways that a Prüfer sequence of appropriate length and
73 // order can be created:
74 // We first compute the number of distributions of total valences (i.e. the distribution of numbers
75 // of free spaces in the sequence) on the interior vertices p_0,...,p_k as Integer points of a
76 // polytope.
77 // For each such valence distribution we compute the number of ways to realize it:
78 // p_0 has to fill the first position to create an ordered Prüfer sequence and then
79 // we have (# of remaining space CHOOSE valence of p_0) other possibilities. Multiplying over all
80 // p_i gives all realizations of the valence distribution and summing over all valence distributions
81 // gives all possibilities.
82
83 Matrix<Rational> eq(0,vertex_count+1);
84 Vector<Rational> eqvec = ones_vector<Rational>(vertex_count);
85 eqvec = Rational(-seq_length) | eqvec;
86 eq /= eqvec;
87
88 Matrix<Rational> ineq = unit_matrix<Rational>(vertex_count);
89 ineq = ( -2 * ones_vector<Rational>(vertex_count)) | ineq;
90
91 BigObject p("polytope::Polytope",
92 "INEQUALITIES", ineq,
93 "EQUATIONS", eq);
94 Matrix<Integer> latt = p.call_method("LATTICE_POINTS");
95 latt = latt.minor(All, range_from(1));
96
97 Integer total(0);
98 for (Int l = 0; l < latt.rows(); ++l) {
99 Integer prod(1);
100 Int sum_vi = 0;
101 for (Int v = 0; v < vertex_count-1; ++v) {
102 Int vi(latt(l, v));
103 prod *= Integer::binom(seq_length-sum_vi-1, vi-1);
104 sum_vi += vi;
105 }
106 total += prod;
107 }
108 return total;
109 }
110
111 ///////////////////////////////////////////////////////////////////////////////////////
112
113 //Documentation see perl wrapper
count_mn_rays(Int n)114 Integer count_mn_rays(Int n)
115 {
116 if (n == 3) {
117 return Integer(0);
118 }
119 Integer result(0);
120 Integer nint(n);
121 for (long i = 1; i <= n-3; ++i) {
122 result += Integer::binom(nint-1, i);
123 }
124 return result;
125 }
126
127 ///////////////////////////////////////////////////////////////////////////////////////
128
129 /**
130 @brief Does exactly the same as count_mn_rays, but returns an Int. Only works for n<=12, since larger values produce too large integers
131 */
count_mn_rays_int(Int n)132 Int count_mn_rays_int(Int n)
133 {
134 if (n == 3) {
135 return 0;
136 }
137 Int result = 0;
138 Int nint = n;
139 for (Int i = 1; i <= n-3; ++i) {
140 result += Int(Integer::binom(nint-1,i));
141 }
142 return result;
143 }
144
145 ///////////////////////////////////////////////////////////////////////////////////////
146
147 // Documentation see header
pair_index_map(Int n)148 Matrix<Int> pair_index_map(Int n)
149 {
150 Matrix<Int> E(n,n);
151 Int nextindex = 0;
152 for (Int i = 0; i < n-1; ++i) {
153 for (Int j = i+1; j < n; ++j) {
154 E(i,j) = E(j,i) = nextindex;
155 ++nextindex;
156 }
157 }
158 return E;
159 }
160
161 ///////////////////////////////////////////////////////////////////////////////////////
162
163 // Documentation see header
decodePrueferSequence(const Vector<Int> & pseq,Int n)164 Vector<Set<Int>> decodePrueferSequence(const Vector<Int>& pseq, Int n)
165 {
166 // Construct vertex set
167 if (n < 0) n = pseq[0]; // The first element is always the number of leaves
168 // Compute number of bounded edges
169 Int no_of_edges = pseq.dim()-n+1;
170 Set<Int> V = sequence(0, n+no_of_edges+1);
171 Vector<Set<Int>> adjacencies(no_of_edges+1); // Which leaves lie "behind" which interior vertex?
172 Vector<Set<Int>> result;
173 const auto allLeafs = sequence(0,n);
174
175 Int firstindex = 0; //We pretend that pseq starts at this index
176 // Connect leaves
177 for (Int i = 0; i < n; ++i) {
178 adjacencies[pseq[firstindex]-n] += i;
179 V -= i;
180 ++firstindex;
181 } //END add leaves
182
183
184 // Now create edges
185 for (Int i = 1; i <= no_of_edges; ++i) {
186 Set<Int> rayset;
187 // If there are only two vertices left, connect them
188 if (i == no_of_edges) {
189 Vector<Int> lasttwo(V);
190 rayset = adjacencies[lasttwo[0]-n];
191 } else {
192 // Find the minimal element in V that is not in the sequence (starting at firstindex)
193 Set<Int> pset(pseq.slice(range_from(firstindex)));
194 Int smallest = -1;
195 for (auto vit = entire(V); !vit.at_end(); ++vit) {
196 if (!pset.contains(*vit)) {
197 smallest = *vit;
198 break;
199 }
200 } //END look for smallest in V\P
201 Set<Int> Av = adjacencies[smallest-n];
202 rayset = Av;
203 adjacencies[pseq[firstindex]-n] += Av;
204 V -= smallest;
205 ++firstindex;
206 }
207
208 // If rayset contains the last leaf, take the complement
209 if (rayset.contains(n-1)) {
210 rayset = allLeafs - rayset;
211 }
212
213 result |= rayset;
214 } //END create edges
215
216 return result;
217 } //END decodePrueferSequence
218
219 ///////////////////////////////////////////////////////////////////////////////////////
220
221 //Documentation see perl wrapper
222 template <typename Addition>
m0n(Int n)223 BigObject m0n(Int n)
224 {
225 if (n == 3) {
226 return projective_torus<Addition>(0,1);
227 }
228 if (n < 3) {
229 throw std::runtime_error("Number of leaves should be at least 3 for M_0,n computation");
230 }
231
232 // First we create the edge index matrix E(i,j) that contains at element i,j the edge index of edge (i,j)
233 // in the complete graph on n-1 nodes
234 Int nextindex = 0;
235 Matrix<Int> E(n-1, n-1);
236 for (Int i = 0; i < n-2; ++i) {
237 for (Int j = i+1; j < n-1; ++j) {
238 E(i,j) = nextindex;
239 E(j,i) = nextindex;
240 ++nextindex;
241 }
242 }
243
244 // We compute the set of all ordered paired Prüfer sequences on n,...,2n-3
245 // (i.e. all sequences of length 2n-4 where each element from n,...,2n-3 occurs twice
246 // and where removing every second occurence of an element yields an ascending sequence)
247 // From each such Prüfer sequence we then construct a maximal cone
248
249 // Will contain the rays of the moduli space in matroid coordinates
250 Int raydim = (n*(n-3))/2 + 1;
251 Int raycount = count_mn_rays_int(n);
252 Matrix<Rational> rays(raycount, raydim);
253
254 // Will contain value 'true' for each ray that has been computed
255 std::vector<bool> raysComputed(count_mn_rays_int(n));
256 //Will contain the set of maximal cones
257 RestrictedIncidenceMatrix<> cones;
258
259 // Compute the number of sequences = number of maximal cones
260 Int noOfMax(count_mn_cones(n, n-3));
261
262 // Things we will need:
263 const auto allLeafs = sequence(0,n); //The complete sequence of leaves (for taking complements)
264 Vector<Int> rayIndices(n-2); //Entry k contains the sum from i = 1 to k of binomial(n-1,i)
265 rayIndices[0] = 0;
266 for (Int i = 1; i < rayIndices.dim(); ++i) {
267 rayIndices[i] = rayIndices[i-1] + Int(Integer::binom(n-1,i));
268 }
269
270 // Iterate through all Prüfer sequences -------------------------------------------------
271
272 Vector<Int> indices = ones_vector<Int>(n-2);
273 Vector<Int> baseSequence(2*n-4);
274 Vector<Set<Int>> adjacent(n-2); // These will be the partitions of the edges
275 Vector<Rational> newray(raydim); // Container for new rays
276 for (Int iteration = 0; iteration < noOfMax; ++iteration) {
277
278 // Create the sequence currently represented by indices and append it------------------
279 baseSequence.resize(2*n-4);
280 baseSequence.fill(0);
281 for (Int i = 0; i < n-1; ++i) {
282 // Go through the non-zero entries of baseSequence. If it is the first or the indices[i]+1-th,
283 // insert an n+i
284 Int nonzero_count = -1;
285 for (Int entry = 0; entry < baseSequence.dim(); ++entry) {
286 if (baseSequence[entry] == 0) {
287 ++nonzero_count;
288 if (nonzero_count == 0) {
289 baseSequence[entry] = n+i;
290 }
291 if (nonzero_count == indices[i]) {
292 baseSequence[entry] = n+i;
293 break;
294 }
295 }
296 }
297 }
298
299 // We now decode the Prüfer sequence to obtain the corresponding cone---------------------
300 Set<Int> newcone;
301
302 Set<Int> V = sequence(0,2*n-2);
303 adjacent.fill(Set<Int>());
304 // First: Connect the leaves
305 for (Int i = 0; i < n; ++i) {
306 adjacent[baseSequence[0]-n] += i;
307 V -= i;
308 baseSequence = baseSequence.slice(range_from(1));
309 }
310 // Now create edges:
311 Int enumber = n-3;
312 for (Int i = 1; i <= enumber; ++i) {
313 // Construct the leaf partition represented by the curve corresponding to the sequence
314 Set<Int> rayset;
315 if (i == enumber) { // If V only has two elements left, simply connect these
316 rayset = adjacent[V.front() - n];
317 } else {
318 Set<Int> pset(baseSequence);
319 Int smallest = -1;
320 // Find the smallest element in V that is not in P
321 for (auto vit = entire(V); !vit.at_end(); ++vit) {
322 if (!pset.contains(*vit)) {
323 smallest = *vit;
324 break;
325 }
326 }
327 rayset = adjacent[smallest-n];
328 // Add the leaves of this partition to the adjacency of the newly connected p_i
329 adjacent[baseSequence[0]-n] += adjacent[smallest-n];
330 // Remove v and p_i
331 V -= smallest;
332 baseSequence = baseSequence.slice(range_from(1));
333 }
334 // The new edge is: v_{adjacent[smallest]}. If it containst the last leaf, take the complement
335 if (rayset.contains(n-1)) {
336 rayset = allLeafs - rayset;
337 }
338
339 // Now we compute the index of the ray -----------------------------------------
340 // Consider ray as a vector of length n filled with a's and b's where entry i is a iff i is in I
341 Int k = n - rayset.size();
342 Int bsleft = k-1;
343 Int l = 1;
344 Int rIndex = rayIndices[k-2];
345 while (bsleft > 1) {
346 if (rayset.contains(n-l-1)) {
347 rIndex += Int(Integer::binom(n-l-1,bsleft-1));
348 } else {
349 --bsleft;
350 }
351 ++l;
352 }
353 Int m = 0;
354 while (rayset.contains(m)) { ++m; }
355 // at last we add the difference of the indices of the second b' and the first b (-1)
356 rIndex += (n-1-l)-m;
357 newcone += rIndex;
358
359
360 // If not, create the corresponding matroid coordinates
361 if (!raysComputed[rIndex]) {
362 raysComputed[rIndex] = true;
363 newray.fill(0);
364 for (auto raypair = entire(all_subsets_of_k(rayset,2)); !raypair.at_end(); ++raypair) {
365 Int newrayindex = E((*raypair).front(),(*raypair).back());
366 // If the newrayindex is one higher than the ray dimension,
367 // this means it is the last pair. Also, we don't
368 // add -e_n but e_1 + ... + e_{n-1} (as we mod out lineality)
369 newray[newrayindex] = Addition::orientation();
370 }
371 rays.row(rIndex) = newray;
372 }
373 } //END iterate edges
374 cones /= newcone;
375
376 // Increase the indices vector by "1"---------------------------------------------------
377 if (iteration < noOfMax-1) {
378 Int counterindex = n-3;
379 while (indices[counterindex] == 2*(n-counterindex)-5) {
380 indices[counterindex] = 1;
381 --counterindex;
382 }
383 ++indices[counterindex];
384 }
385 } //END iterate cones
386
387 // Add the vertex at the origin
388 rays = zero_vector<Rational>(rays.rows()) | rays;
389
390 rays /= unit_vector<Rational>(rays.cols(), 0);
391
392 // Add the vertex to all cones
393 const Int vertex = rays.rows() - 1;
394 for (auto mc = entire(rows(cones)); !mc.at_end(); ++mc) {
395 *mc += vertex;
396 }
397
398 const IncidenceMatrix<> result_cones{std::move(cones)};
399
400 BigObject result("Cycle", mlist<Addition>(),
401 "WEIGHTS", ones_vector<Int>(result_cones.rows()),
402 "PROJECTIVE_VERTICES", rays,
403 "MAXIMAL_POLYTOPES", result_cones);
404 result.set_description() << "Moduli space M_0," << n;
405 return result;
406 }
407
408
409 template <typename Addition>
space_of_stable_maps(Int n,Int d,Int r)410 BigObject space_of_stable_maps(Int n, Int d, Int r)
411 {
412 BigObject moduli = m0n<Addition>(n+d);
413 BigObject torus = projective_torus<Addition>(r,1);
414 BigObject result = call_function("cartesian_product", moduli, torus);
415 result.set_description() << "Moduli space of stable rational maps with " << n << " contracted ends, " << d << " non-contracted ends into the torus of dimension " << d;
416 return result;
417 }
418
419 // ------------------------- PERL WRAPPERS ---------------------------------------------------
420
421 UserFunction4perl("# @category Moduli of rational curves"
422 "# Computes the number of k-dimensional cones of the tropical moduli space M_0,n"
423 "# @param Int n The number of leaves. Should be >= 3"
424 "# @param Int k The number of bounded edges. This argument is optional and n-3 by default"
425 "# @return Integer The number of k-dimensional cones of M_0,n",
426 &count_mn_cones,"count_mn_cones($;$=$_[0]-3)");
427
428
429 UserFunction4perl("# @category Moduli of rational curves"
430 "# Computes the number of rays of the tropical moduli space M_0,n"
431 "# @param Int n The number of leaves. Should be >= 3"
432 "# @return Integer The number of rays",
433 &count_mn_rays,"count_mn_rays($)");
434
435 UserFunctionTemplate4perl("# @category Moduli of rational curves"
436 "# Creates the moduli space of abstract rational n-marked curves. Its coordinates are"
437 "# given as the coordinates of the bergman fan of the matroid of the complete graph on "
438 "# n-1 nodes (but not computed as such)"
439 "# The isomorphism to the space of curve metrics is obtained by choosing"
440 "# the last leaf as special leaf"
441 "# @param Int n The number of leaves. Should be at least 3"
442 "# @tparam Addition Min or Max"
443 "# @return Cycle The tropical moduli space M_0,n",
444 "m0n<Addition>($)");
445
446 UserFunctionTemplate4perl("# @category Moduli of rational curves"
447 "# Creates the moduli space of stable maps of rational n-marked curves into a "
448 "# projective torus. It is given as the cartesian product of M_{0,n+d} and R^r,"
449 "# where n is the number of contracted leaves, d the number of non-contracted leaves"
450 "# and r is the dimension of the target torus. The R^r - coordinate is interpreted as "
451 "# the image of the last (n-th) contracted leaf."
452 "# Due to the implementation of [[cartesian_product]], the projective coordinates are"
453 "# non-canonical: Both M_{0,n+d} and R^r are dehomogenized after the first coordinate, then"
454 "# the product is taken and homogenized after the first coordinate again."
455 "# Note that functions in a-tint will usually treat this space in such a way that the"
456 "# first d leaves are the non-contracted ones and the remaining n leaves are the "
457 "# contracted ones."
458 "# @param Int n The number of contracted leaves"
459 "# @param Int d The number of non-contracted leaves"
460 "# @param Int r The dimension of the target space for the stable maps."
461 "# @tparam Addition Min or Max. Determines the coordinates."
462 "# @return Cycle The moduli space of rational stable maps.",
463 "space_of_stable_maps<Addition>($,$,$)");
464
465 //Function4perl(&decodePrueferSequence,"dcp(Vector<Int>;$=-1)");
466 // UserFunction4perl("",&adjacentRays,"adjacentRays(RationalCurve)");
467
468 } }
469