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