1 /* Copyright (c) 1997-2021
2 Ewgenij Gawrilow, Michael Joswig, and the polymake team
3 Technische Universität Berlin, Germany
4 https://polymake.org
5
6 This program is free software; you can redistribute it and/or modify it
7 under the terms of the GNU General Public License as published by the
8 Free Software Foundation; either version 2, or (at your option) any
9 later version: http://www.gnu.org/licenses/gpl.txt.
10
11 This program is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15 --------------------------------------------------------------------------------
16 */
17
18 #include "polymake/client.h"
19 #include "polymake/Vector.h"
20 #include "polymake/Matrix.h"
21 #include "polymake/ListMatrix.h"
22 #include "polymake/tropical/dual_addition_version.h"
23 #include "polymake/TropicalNumber.h"
24
25 namespace polymake { namespace tropical {
26
27 // Constructs the dual polynomial of hyperplanes defined by the [[POINTS]] of a cone.
28 template <typename Addition, typename Scalar>
29 Polynomial<TropicalNumber<typename Addition::dual,Scalar>>
cone_polynomial(const Matrix<TropicalNumber<Addition,Scalar>> & points)30 cone_polynomial(const Matrix<TropicalNumber<Addition,Scalar> > &points)
31 {
32 using TNumber = TropicalNumber<typename Addition::dual, Scalar>;
33
34 // Invert the points
35 const Matrix<TNumber> dual_points = dual_addition_version(points,1);
36
37 // Construct the linear polynomials
38 Polynomial<TNumber> h (TNumber::one(), points.cols());
39 for (auto vec=entire(rows(dual_points)); !vec.at_end(); ++vec) {
40 h *= Polynomial<TNumber>(*vec, unit_matrix<Int>(points.cols()));
41 }
42
43 return h;
44 }
45
46 // Constructs the dome of the above Polynomial.
47 template <typename Addition, typename Scalar>
dome_hyperplane_arrangement(const Matrix<TropicalNumber<Addition,Scalar>> & points)48 BigObject dome_hyperplane_arrangement(const Matrix<TropicalNumber<Addition,Scalar>>& points)
49 {
50 using TNumber = TropicalNumber<typename Addition::dual, Scalar>;
51
52 Polynomial<TNumber> h = cone_polynomial(points);
53
54 const Matrix<Int> monoms_int = h.monomials_as_matrix();
55 Matrix<Rational> monoms(monoms_int); // cast coefficients Integer -> Rational
56 const Vector< TNumber > coefs=h.coefficients_as_vector();
57 const Int d = monoms.cols();
58 const Int n = monoms.rows();
59
60 // We have to make all exponents positive, otherwise the below equations produce
61 // a wrong result. We multiply the polynomial with a single monomial, which
62 // does not change the hypersurface.
63 Vector<Rational> min_degrees(monoms.cols());
64 for (Int v = 0; v < monoms.cols(); ++v) {
65 min_degrees[v] = accumulate(monoms.col(v),operations::min());
66 // If the minimal degree is positive, we're good
67 min_degrees[v] = std::min(min_degrees[v],Rational(0));
68 }
69 for (Int m = 0; m < monoms.rows(); ++m) {
70 monoms.row(m) -= min_degrees;
71 }
72
73 // dual to extended Newton polyhedron
74 ListMatrix< Vector<Rational> > ineq;
75 const TNumber zero=TNumber::zero();
76 for (Int i = 0; i < n; ++i) {
77 if (coefs[i]==zero)
78 ineq /= unit_vector<Rational>(d+1,0);
79 else
80 ineq /= (-1)*Addition::orientation()*(Rational(coefs[i])|monoms[i]);
81 }
82
83 return BigObject("polytope::Polytope", mlist<Scalar>(),
84 "INEQUALITIES", ineq,
85 "FEASIBLE", true,
86 "BOUNDED", false);
87 }
88
89 FunctionTemplate4perl("cone_polynomial<Addition,Scalar>(Matrix<TropicalNumber<Addition, Scalar>>)");
90 FunctionTemplate4perl("dome_hyperplane_arrangement<Addition,Scalar>(Matrix<TropicalNumber<Addition, Scalar>>)");
91
92 } }
93