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 all [[LATTICE...]- related properties
27 */
28
29 #include "polymake/client.h"
30 #include "polymake/Set.h"
31 #include "polymake/Array.h"
32 #include "polymake/Matrix.h"
33 #include "polymake/SparseMatrix.h"
34 #include "polymake/IncidenceMatrix.h"
35 #include "polymake/Vector.h"
36 #include "polymake/Rational.h"
37 #include "polymake/Map.h"
38 #include "polymake/linalg.h"
39 #include "polymake/common/lattice_tools.h"
40 #include "polymake/integer_linalg.h"
41 #include "polymake/tropical/thomog.h"
42 #include "polymake/tropical/lattice.h"
43 #include "polymake/tropical/linear_algebra_tools.h"
44 #include "polymake/tropical/misc_tools.h"
45
46
47 namespace polymake { namespace tropical {
48
49 using LatticeMap = Map< std::pair<Int, Int>, Vector<Integer>>;
50
51 /*
52 * @brief Computes [[LATTICE_NORMAL_SUM]]
53 */
computeLatticeNormalSum(BigObject cycle)54 void computeLatticeNormalSum(BigObject cycle)
55 {
56 const LatticeMap& latticeNormals = cycle.give("LATTICE_NORMALS");
57 const Int ambient_dim = cycle.give("FAN_AMBIENT_DIM");
58 const Vector<Integer> &weights = cycle.give("WEIGHTS");
59 const IncidenceMatrix<> &codimInc = cycle.give("MAXIMAL_AT_CODIM_ONE");
60
61 // This will contain the result
62 ListMatrix< Vector<Integer> > summatrix(0, ambient_dim);
63
64 // Iterate over all codim one faces
65 for (auto facet = entire<indexed>(rows(codimInc)); !facet.at_end(); ++facet) {
66 // This will contain the weighted sum of the lattice normals
67 Vector<Integer> result = zero_vector<Integer>(ambient_dim);
68 // Go through all adjacent cones
69 for (auto e = entire(facet->top()); !e.at_end(); ++e) {
70 result = result +latticeNormals[std::make_pair(facet.index(),*e)] * weights[*e];
71 }
72 summatrix /= result;
73 }
74
75 cycle.take("LATTICE_NORMAL_SUM") << summatrix;
76
77 } // END computeLatticeNormalSum
78
79 /*
80 * @brief Computes properties [[LATTICE_NORMAL_FCT_VECTOR]], [[LATTICE_NORMAL_SUM_FCT_VECTOR]]
81 */
computeLatticeFunctionData(BigObject cycle)82 void computeLatticeFunctionData(BigObject cycle)
83 {
84 // Extract properties from the cycle
85 const Matrix<Rational> &linealitySpace_ref = cycle.give("LINEALITY_SPACE");
86 const Matrix<Rational> linealitySpace = tdehomog(linealitySpace_ref);
87 const Int lineality_dim = linealitySpace.rows();
88
89 const Matrix<Rational>& rays_ref = cycle.give("SEPARATED_VERTICES");
90 const Matrix<Rational> rays = tdehomog(rays_ref);
91 const LatticeMap &latticeNormals = cycle.give("LATTICE_NORMALS");
92 const Matrix<Rational>& normalsums_ref = cycle.give("LATTICE_NORMAL_SUM");
93 const Matrix<Rational> normalsums = tdehomog(normalsums_ref);
94 const IncidenceMatrix<>& codimOneCones = cycle.give("SEPARATED_CODIMENSION_ONE_POLYTOPES");
95 const IncidenceMatrix<>& maximalCones = cycle.give("SEPARATED_MAXIMAL_POLYTOPES");
96 const IncidenceMatrix<>& coneIncidences = cycle.give("MAXIMAL_AT_CODIM_ONE");
97
98 // Result variables
99 Map<std::pair<Int, Int>, Vector<Rational>> summap;
100 ListMatrix<Vector<Rational>> summatrix;
101
102 // Iterate over all codim 1 faces
103 auto coInc = entire(rows(coneIncidences));
104 for (auto fct = entire<indexed>(rows(codimOneCones)); !fct.at_end(); ++fct, ++coInc) {
105 for (auto mc_it = entire(*coInc); !mc_it.at_end(); ++mc_it) {
106 const Int mc = *mc_it;
107 Vector<Rational> normalvector(latticeNormals[std::make_pair(fct.index(), mc)]);
108 normalvector = tdehomog_vec(normalvector);
109 // Compute the representation of the normal vector
110 summap[std::make_pair(fct.index(), mc)] =
111 functionRepresentationVector(maximalCones.row(mc),
112 normalvector,
113 rays, linealitySpace);
114 }
115
116 // Now compute the representation of the sum of the normals
117 try {
118 summatrix /= functionRepresentationVector(*fct,
119 normalsums.row(fct.index()),
120 rays, linealitySpace);
121 }
122 catch (const std::runtime_error& e) {
123 // This goes wrong, if X is not balanced at a given codim 1 face
124 summatrix /= zero_vector<Rational>(rays.rows() + lineality_dim);
125 }
126 }
127
128 // Set fan properties
129 cycle.take("LATTICE_NORMAL_FCT_VECTOR") << summap;
130 cycle.take("LATTICE_NORMAL_SUM_FCT_VECTOR") << summatrix;
131 } // END computeLatticeFunctionData
132
133
lattice_basis_of_cone(const Matrix<Rational> & rays,const Matrix<Rational> & lineality,Int dim,bool has_leading_coordinate)134 Matrix<Integer> lattice_basis_of_cone(const Matrix<Rational>& rays,
135 const Matrix<Rational> &lineality,
136 Int dim, bool has_leading_coordinate)
137 {
138 // Special case: If the cone is full-dimensional, return the standard basis
139 Int ambient_dim = std::max(rays.cols(), lineality.cols()) - (has_leading_coordinate ? 1 : 0);
140 if (dim == ambient_dim)
141 return unit_matrix<Integer>(ambient_dim);
142
143 // Compute span matrix if there is a leading coordinate
144 Matrix<Rational> span_matrix(0, ambient_dim);
145 if (has_leading_coordinate) {
146 std::pair<Set<Int>, Set<Int>> sortedVertices = far_and_nonfar_vertices(rays);
147 span_matrix = rays.minor(sortedVertices.first, range_from(1));
148 if(lineality.rows() > 0) span_matrix /= lineality.minor(All, range_from(1));
149 Int first = *(sortedVertices.second.begin());
150 sortedVertices.second -= first;
151 for (auto vtx = entire(rows(rays.minor(sortedVertices.second,All))); !vtx.at_end(); ++vtx) {
152 span_matrix /= (*vtx - rays.row(first)).slice(range_from(1));
153 }
154 } else {
155 span_matrix = rays / lineality;
156 }
157
158 // Compute span of cone
159 Matrix<Rational> linspan = null_space(span_matrix);
160 SparseMatrix<Integer> transformation =
161 pm::hermite_normal_form( common::eliminate_denominators_in_rows(linspan),false).companion;
162
163 // The last dim columns are a Z-basis for the cone
164 return Matrix<Integer>(T(transformation.minor(All,sequence(transformation.cols()-dim,dim))));
165 } // END lattice_basis_of_cone
166
167 /*
168 * @brief Computes properties [[LATTICE_BASES]] and [[LATTICE_GENERATORS]]
169 */
computeLatticeBases(BigObject cycle)170 void computeLatticeBases(BigObject cycle)
171 {
172 // Extract properties
173 const Matrix<Rational> &rays_ref = cycle.give("VERTICES");
174 const Matrix<Rational> rays = tdehomog(rays_ref).minor(All, range_from(1));
175 const Matrix<Rational> &linspace_ref = cycle.give("LINEALITY_SPACE");
176 const Matrix<Rational> linspace = tdehomog(linspace_ref).minor(All, range_from(1));
177 const IncidenceMatrix<> &cones = cycle.give("MAXIMAL_POLYTOPES");
178 // FIXME Use unimodularity?
179 const Set<Int> &directional = cycle.give("FAR_VERTICES");
180 const Set<Int> vertices = sequence(0,rays.rows()) - directional;
181 const Int dim = cycle.give("PROJECTIVE_DIM");
182
183 Matrix<Integer> generators;
184 std::vector<Set<Int>> bases;
185
186 // Iterate over all cones
187 for (auto mc = entire(rows(cones)); !mc.at_end(); ++mc) {
188 // Compute a lattice basis for the cone:
189 // Construct ray matrix of cone:
190 Matrix<Rational> mc_rays = rays.minor((*mc) * directional, All);
191 Matrix<Rational> mc_vert = rays.minor((*mc) * vertices,All);
192 for (Int v = 1; v < mc_vert.rows(); v++) {
193 mc_rays /= (mc_vert.row(v) - mc_vert.row(0));
194 }
195
196 Matrix<Integer> basis;
197 // FIXME Use unimodularity?
198 basis = lattice_basis_of_cone(mc_rays, linspace,dim,false);
199 basis = zero_vector<Integer>(basis.rows()) | basis;
200
201 Set<Int> basis_set;
202 // Add rays, looking for doubles
203 for (auto b = entire(rows(basis)); !b.at_end(); ++b) {
204 // We normalize s.t. the first non-zero entry is > 0
205 auto first_non_zero = find_in_range_if(entire(b->top()), operations::non_zero());
206 if (*first_non_zero < 0) b->negate();
207 Int ray_index = -1;
208 for (auto gr = entire<indexed>(rows(generators)); !gr.at_end(); ++gr) {
209 if (*gr == *b) {
210 ray_index = gr.index(); break;
211 }
212 }
213 if (ray_index == -1) {
214 generators /= *b;
215 ray_index = generators.rows()-1;
216 }
217 basis_set += ray_index;
218 } // END go through basis elements
219 bases.push_back(basis_set);
220 } // END iterate over all maximal cones
221
222 // Set properties
223 cycle.take("LATTICE_GENERATORS") << thomog(generators);
224 cycle.take("LATTICE_BASES") << bases;
225 }
226
227 Function4perl(&computeLatticeNormalSum,"computeLatticeNormalSum(Cycle)");
228 Function4perl(&computeLatticeFunctionData,"computeLatticeFunctionData(Cycle)");
229 Function4perl(&computeLatticeBases, "computeLatticeBases(Cycle)");
230 Function4perl(&lattice_basis_of_cone,"lattice_basis_of_cone(Matrix,Matrix,$,$)");
231
232 } }
233