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 geometric properties of a MatroidRingCycle.
27 */
28
29 #include "polymake/client.h"
30 #include "polymake/Rational.h"
31 #include "polymake/Matrix.h"
32 #include "polymake/Vector.h"
33 #include "polymake/linalg.h"
34 #include "polymake/IncidenceMatrix.h"
35 #include "polymake/tropical/misc_tools.h"
36 #include "polymake/tropical/thomog.h"
37 #include "polymake/tropical/specialcycles.h"
38
39 namespace polymake { namespace tropical {
40
41 // Finds a vertex in a list of row vectors.
find_index(const Vector<Rational> & v,const Matrix<Rational> & m)42 Int find_index(const Vector<Rational>& v, const Matrix<Rational>& m)
43 {
44 Int i = 0;
45 for (auto r = entire(rows(m)); !r.at_end(); ++r, ++i) {
46 if (*r == v) return i;
47 }
48 throw std::runtime_error("Vertex not found");
49 }
50
51 /*
52 * @brief Takes a list of tropical cycles which are all defined on (a subset of)
53 * the same polyhedral structure and computes their cycle sum
54 * @return Cycle<Addition>
55 */
56 template <typename Addition>
add_refined_cycles(Array<BigObject> cycles)57 BigObject add_refined_cycles(Array<BigObject> cycles)
58 {
59 Array<Matrix<Rational>> vertices(cycles.size());
60 Matrix<Rational> vertex_union;
61 Array<IncidenceMatrix<>> cones(cycles.size());
62 Array<Vector<Integer> > weights(cycles.size());
63 for (Int i =0; i < cycles.size(); ++i) {
64 cycles[i].give("VERTICES") >> vertices[i];
65 vertex_union /= vertices[i];
66 cycles[i].give("MAXIMAL_POLYTOPES") >> cones[i];
67 cycles[i].give("WEIGHTS") >> weights[i];
68 }
69 Set<Vector<Rational> > vset (rows(Matrix<Rational>(vertex_union)));
70 Matrix<Rational> vertices_total(vset);
71
72 // Renumber vertices and map cones
73 Vector<Set<Int>> new_cones;
74 Map<Set<Int>, Int> new_cone_indices;
75 Int next_index = 0;
76 Vector<Integer> new_weights;
77 for (Int i = 0; i < cycles.size(); ++i) {
78 Map<Int, Int> vmap;
79 for (Int r =0; r < vertices[i].rows(); ++r) {
80 vmap[r] = find_index(vertices[i].row(r), vertices_total);
81 }
82 Int cindex = 0;
83 for (auto c = entire(rows(cones[i])); !c.at_end(); ++c, ++cindex) {
84 Set<Int> mapped_cone{ vmap.map(*c) };
85 Integer cw = (weights[i])[cindex];
86 if (!new_cone_indices.contains(mapped_cone)) {
87 new_cone_indices[mapped_cone] = next_index;
88 new_cones |= mapped_cone;
89 new_weights |= cw;
90 ++next_index;
91 } else {
92 new_weights[ new_cone_indices[mapped_cone] ] += cw;
93 }
94 }
95 }
96
97 // Remove zero-weight cones
98 Set<Int> supp = support(new_weights);
99 Set<Int> used_vertices = accumulate( new_cones.slice(supp), operations::add());
100 IncidenceMatrix<> new_cones_matrix(new_cones);
101
102 return BigObject("Cycle", mlist<Addition>(),
103 "VERTICES", vertices_total.minor(used_vertices, All),
104 "MAXIMAL_POLYTOPES", new_cones_matrix.minor(supp, used_vertices),
105 "WEIGHTS", new_weights.slice(supp));
106 }
107
108 FunctionTemplate4perl("add_refined_cycles<Addition>(Cycle<Addition>+)");
109
110 } }
111