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 Implementations of miscellaneous tools
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/IncidenceMatrix.h"
35 #include "polymake/linalg.h"
36 #include "polymake/tropical/thomog.h"
37 #include "polymake/tropical/morphism_values.h"
38 #include "polymake/tropical/homogeneous_convex_hull.h"
39 #include "polymake/tropical/misc_tools.h"
40
41 namespace polymake { namespace tropical {
42
43 using matrix_pair = std::pair<Matrix<Rational>, Matrix<Rational>>;
44
all_cones_as_incidence(BigObject complex)45 IncidenceMatrix<> all_cones_as_incidence(BigObject complex)
46 {
47 Array<IncidenceMatrix<>> all_cones = complex.give("CONES");
48 if (all_cones.size() == 0) return IncidenceMatrix<>();
49 return IncidenceMatrix<>(rowwise(), all_cones);
50 }
51
binaryMatrix(Int n)52 Matrix<Rational> binaryMatrix(Int n)
53 {
54 ListMatrix<Vector<Rational>> result(0,n);
55 Vector<Rational> prev_row = -ones_vector<Rational>(n);
56 result /= prev_row;
57 // Now increase the last row of result by "one" in each iteration and append the new row to result
58 Integer iterations = Integer::pow(2,n)-1;
59 for (Int i = 1; i <= iterations; ++i) {
60 // Find the first -1-entry
61 Vector<Rational> newrow(prev_row);
62 auto neg_it = find_in_range_if(entire(newrow), operations::negative());
63 // Now toggle this and all previous entries
64 auto nr_it = entire(newrow);
65 while (nr_it != neg_it) {
66 *nr_it = -1;
67 ++nr_it;
68 }
69 *nr_it = 1;
70 result /= newrow;
71 prev_row = newrow;
72 }
73 return result;
74 }
75
76 /**
77 @brief Computes the labels for each cell of a function domain, given its ray and lineality values
78 @param BigObject domain A Cycle object, representing the domain of the function
79 @param Matrix<Rationl> ray_values Values of the function on the rays, given as row vectors in non-homog. coordinates
80 @param Matrix<Rational> lin_values Values of the function on the lineality space, given as row vectors in non-homog. coordinates.
81 @param bool values_are_homogeneous Whether vertex and lineality values are given in tropical projective
82 coordinates.
83 @return A list of std::strings
84 */
computeFunctionLabels(BigObject domain,Matrix<Rational> ray_values,Matrix<Rational> lin_values,const bool values_are_homogeneous)85 ListReturn computeFunctionLabels(BigObject domain, Matrix<Rational> ray_values, Matrix<Rational> lin_values, const bool values_are_homogeneous)
86 {
87 // Extract values
88 const Matrix<Rational>& rays_ref = domain.give("SEPARATED_VERTICES");
89 const Matrix<Rational> rays = tdehomog(rays_ref,0);
90 const IncidenceMatrix<>& cones = domain.give("SEPARATED_MAXIMAL_POLYTOPES");
91 const Matrix<Rational>& lineality_ref = domain.give("LINEALITY_SPACE");
92 const Matrix<Rational> lineality = tdehomog(lineality_ref,0);
93
94 if (values_are_homogeneous) {
95 ray_values = tdehomog(ray_values,0);
96 lin_values = tdehomog(lin_values,0);
97 }
98
99 ListReturn result;
100
101 for (auto mc = entire(rows(cones)); !mc.at_end(); ++mc) {
102 Matrix<Rational> matrix;
103 Vector<Rational> translate;
104 computeConeFunction(rays.minor(*mc,All), lineality, ray_values.minor(*mc,All), lin_values, translate, matrix);
105
106 matrix = thomog(matrix,0,false);
107 translate = thomog_vec(translate,0,false);
108
109 std::ostringstream rep;
110 if (matrix.rows() > 1) {
111 wrap(rep) << "(" << translate << ")" << " + ";
112 for (auto i = entire(rows(matrix)); !i.at_end(); ++i)
113 wrap(rep) << "[" << (*i) << "]";
114 }
115 // We have a special representation format for functions to R
116 else {
117 bool hadnonzeroterm = false;
118 for (Int i = 0; i < matrix.cols(); ++i) {
119 if (matrix(0,i) != 0) {
120 if (hadnonzeroterm) rep << " + ";
121 hadnonzeroterm = true;
122 if (matrix(0,i) < 0)
123 rep << "(" << matrix(0,i) << ")";
124 else
125 rep << matrix(0,i);
126 rep << "*x_" << (i+1);
127 }
128 }
129 if (translate[0] < 0 && hadnonzeroterm) rep << " - " << -translate[0];
130 if (translate[0] > 0 && hadnonzeroterm) rep << " + " << translate[0];
131 if (!hadnonzeroterm) rep << translate[0];
132 }
133 result << rep.str();
134 }
135
136 return result;
137 }
138
139 // Documentation see perl wrapper
contains_point(BigObject complex,const Vector<Rational> & point)140 bool contains_point(BigObject complex, const Vector<Rational>& point)
141 {
142 // Special case: Empty cycle
143 if (call_function("is_empty", complex))
144 return false;
145
146 // Extract values
147 Matrix<Rational> rays = complex.give("VERTICES");
148 Matrix<Rational> linspace = complex.give("LINEALITY_SPACE");
149 IncidenceMatrix<> cones = complex.give("MAXIMAL_POLYTOPES");
150
151 if (point.dim() != rays.cols() && point.dim() != linspace.cols()) {
152 throw std::runtime_error("Point does not have the same ambient dimension as the complex.");
153 }
154
155 for (Int mc = 0; mc < cones.rows(); ++mc) {
156 if (is_ray_in_cone(rays.minor(cones.row(mc),All), linspace, point, true))
157 return true;
158 }
159
160 return false;
161 }
162
163 UserFunction4perl("# @category Lattices"
164 "# Returns n random integers in the range 0.. (max_arg-1),inclusive"
165 "# Note that this algorithm is not optimal for real randomness:"
166 "# If you change the range parameter and then change it back, you will"
167 "# usually get the exact same sequence as the first time"
168 "# @param Int max_arg The upper bound for the random integers"
169 "# @param Int n The number of integers to be created"
170 "# @return Vector<Integer>",
171 &randomInteger,"randomInteger($, $)");
172
173 UserFunction4perl("# @category Basic polyhedral operations"
174 "# Takes a weighted complex and a point and computed whether that point lies in "
175 "# the complex"
176 "# @param Cycle A weighted complex"
177 "# @param Vector<Rational> point An arbitrary vector in the same ambient"
178 "# dimension as complex. Given in tropical projective coordinates with leading coordinate."
179 "# @return Bool Whether the point lies in the support of complex",
180 &contains_point,"contains_point(Cycle,$)");
181
182 Function4perl(&computeFunctionLabels, "computeFunctionLabels(Cycle, Matrix,Matrix,$)");
183
184 } }
185