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 The functions here deal with all the [[SEPARATED..]] properties
27 */
28
29 #include "polymake/client.h"
30 #include "polymake/Set.h"
31 #include "polymake/Matrix.h"
32 #include "polymake/Vector.h"
33 #include "polymake/Rational.h"
34 #include "polymake/Array.h"
35 #include "polymake/IncidenceMatrix.h"
36 #include "polymake/tropical/separated_data.h"
37
38
39 namespace polymake { namespace tropical {
40
computeSeparatedData(BigObject X)41 void computeSeparatedData(BigObject X)
42 {
43 // Extract properties of X
44 Matrix<Rational> rays = X.give("VERTICES");
45
46 IncidenceMatrix<> codimOneCones = X.give("CODIMENSION_ONE_POLYTOPES");
47 IncidenceMatrix<> maximalCones = X.give("MAXIMAL_POLYTOPES");
48 IncidenceMatrix<> facet_incidences = X.give("MAXIMAL_AT_CODIM_ONE");
49 facet_incidences = T(facet_incidences);
50
51 // Result variables
52 Matrix<Rational> cmplxrays(0,rays.cols());
53 Vector<Set<Int>> maxcones(maximalCones.rows());
54 Vector<Set<Int>> codimone(codimOneCones.rows());
55 Vector<Int> conversion;
56
57 // Divide the set of rays into those with x0 != 0 and those with x0 = 0
58 Set<Int> affineRays;
59 Set<Int> directionalRays;
60 Map<Int, Int> newAffineIndices; //This maps the old ray indices to the new ones in cmplxrays
61 for (Int r = 0; r < rays.rows(); ++r) {
62 if (is_zero(rays.row(r)[0])) {
63 directionalRays = directionalRays + r;
64 } else {
65 affineRays = affineRays + r;
66 cmplxrays = cmplxrays / rays.row(r);
67 conversion |= r;
68 newAffineIndices[r] = cmplxrays.rows()-1;
69 }
70 }
71
72 // Insert the indices of the new affine rays for each cone
73 for (Int co = 0; co < codimOneCones.rows(); ++co) {
74 Set<Int> corays = codimOneCones.row(co) * affineRays;
75 codimone[co] = Set<Int>();
76 for (auto e = entire( corays); !e.at_end(); ++e) {
77 codimone[co] = codimone[co] + newAffineIndices[*e];
78 }
79 }
80 for (Int mc = 0; mc < maximalCones.rows(); ++mc) {
81 Set<Int> mcrays = maximalCones.row(mc) * affineRays;
82 maxcones[mc] = Set<Int>();
83 for (auto e = entire( mcrays); !e.at_end(); ++e) {
84 maxcones[mc] = maxcones[mc] + newAffineIndices[*e];
85 }
86 }
87
88 // Now we go through the directional rays and compute the connected component for each one
89 for (auto r = entire(directionalRays); !r.at_end(); ++r) {
90
91 // List of connected components of this ray, each element is a component
92 // containing the indices of the maximal cones
93 Vector<Set<Int>> connectedComponents;
94 // The inverse of the component matrix, i.e. maps cone indices to row indices of connectedComponents
95 Map<Int, Int> inverseMap;
96
97 // Compute the set of maximal cones containing r
98 Set<Int> rcones;
99 for (Int mc = 0; mc < maximalCones.rows(); ++mc) {
100 if (maximalCones.row(mc).contains(*r)) {
101 rcones += mc;
102 }
103 }
104
105 // For each such maximal cone, compute its component (if it hasnt been computed yet).
106 for (auto mc = entire(rcones); !mc.at_end(); ++mc) {
107 if (!inverseMap.exists(*mc)) {
108 // Create new component
109 Set<Int> newset{ *mc };
110 connectedComponents |= newset;
111 inverseMap[*mc] = connectedComponents.dim()-1;
112
113 // Do a breadth-first search for all other cones in the component
114 std::list<Int> queue;
115 queue.push_back(*mc);
116 // Semantics: Elements in that queue have been added but their neighbours might not
117 while (!queue.empty()) {
118 Int node = queue.front(); //Take the first element and find its neighbours
119 queue.pop_front();
120 for (auto othercone = entire(rcones); !othercone.at_end(); ++othercone) {
121 // We only want 'homeless' cones
122 if (!inverseMap.exists(*othercone)) {
123 // This checks whether both cones share a ray with x0=1
124 if (!(maximalCones.row(node) * maximalCones.row(*othercone) * affineRays).empty()) {
125 // Add this cone to the component
126 connectedComponents[connectedComponents.dim()-1] += *othercone;
127 inverseMap[*othercone] = connectedComponents.dim()-1;
128 queue.push_back(*othercone);
129 }
130 }
131 }
132 }
133
134 }
135 } //END computation of connected components
136
137 // Now add r once for each connected component to the appropriate cones
138 for (Int cc = 0; cc < connectedComponents.dim(); ++cc) {
139 cmplxrays = cmplxrays / rays.row(*r);
140 conversion |= (*r);
141 Int rowindex = cmplxrays.rows()-1;
142 Set<Int> ccset = connectedComponents[cc];
143 for (auto mc = entire(ccset); !mc.at_end(); ++mc) {
144 maxcones[*mc] = maxcones[*mc] + rowindex;
145 // For each facet of mc that contains r, add rowindex
146 Set<Int> fcset;
147 // If there are maximal cones not intersecting x0 = 1, they have no facets
148 // in facet_incidences, hence the following check
149 if (*mc < facet_incidences.rows()) {
150 fcset = facet_incidences.row(*mc);
151 }
152 for (auto fct = entire(fcset); !fct.at_end(); ++fct) {
153 if (codimOneCones.row(*fct).contains(*r)) {
154 codimone[*fct] += rowindex;
155 }
156 }
157 }
158 }
159
160 } // END iterate over all rays
161
162 X.take("SEPARATED_VERTICES") << cmplxrays;
163 X.take("SEPARATED_MAXIMAL_POLYTOPES") << maxcones;
164 X.take("SEPARATED_CODIMENSION_ONE_POLYTOPES") << codimone;
165 X.take("SEPARATED_CONVERSION_VECTOR") << conversion;
166
167 } // END computeSeparatedData
168
169 Function4perl(&computeSeparatedData,"computeSeparatedData(Cycle)");
170
171 } }
172