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 Contains some miscellaneous tools.
27
28 */
29
30
31 #pragma once
32
33 #include "polymake/Rational.h"
34 #include "polymake/Matrix.h"
35 #include "polymake/Set.h"
36 #include "polymake/IncidenceMatrix.h"
37 #include "polymake/Array.h"
38 #include "polymake/RandomGenerators.h"
39 #include "polymake/linalg.h"
40 #include "polymake/tropical/homogeneous_convex_hull.h"
41
42 namespace polymake { namespace tropical {
43
44 /**
45 @brief Takes a matrix and returns the row indices where the first coordinate is nonzero and where the first coordinate is zero in two different sets
46 @param Matrix<Rational> m The matrix whose rows we consider
47 @return std::pair<Set<Int>, Set<Int>> The first set contains the row indices of rows that start with a zero entry, the second set is the complement
48 */
49 template <typename MType>
far_and_nonfar_vertices(const GenericMatrix<MType> & m)50 std::pair<Set<Int>, Set<Int>> far_and_nonfar_vertices(const GenericMatrix<MType>& m)
51 {
52 const auto& first_col_supp = support(m.col(0));
53 return std::pair<Set<Int>, Set<Int>>( sequence(0, m.rows()) - first_col_supp, first_col_supp);
54 }
55
56 /**
57 * @brief Takes a polyhedral complex and returns [[CONES]] summarized into one single incidence matrix.
58 * @param PolyhedralComplex
59 * @return IncidenceMatrix<>
60 */
61 IncidenceMatrix<> all_cones_as_incidence(BigObject complex);
62
63 /**
64 * @brief Converts an incidence matrix to a Vector<Set<Int>>
65 */
66 template <typename MType>
incMatrixToVector(const GenericIncidenceMatrix<MType> & i)67 Vector<Set<Int>> incMatrixToVector(const GenericIncidenceMatrix<MType>& i)
68 {
69 return Vector<Set<Int>>(i.rows(), entire(rows(i)));
70 }
71
72 inline
randomInteger(const Int max_arg,const Int n)73 Vector<Integer> randomInteger(const Int max_arg, const Int n)
74 {
75 static UniformlyRandomRanged<Integer> rg(max_arg);
76 return Vector<Integer>(n, rg.begin());
77 }
78
79 /**
80 @brief Computes all vectors of dimension n with entries +1 and -1.
81 They are sorted such that each vector v has the row index determined by the sum:
82 sum_{i: v_i = 1} 2^i (where i runs from 0 to n-1)
83 @param Int n The column dimension of the matrix
84 @return Matrix<Rational> A 2^n by n matrix containing all +-1-vectors of dimension n
85 */
86 Matrix<Rational> binaryMatrix(Int n);
87
88 /**
89 @brief Assumes v is a vector with entries +1 and -1 only.
90 Returns sum_{i: v_i = 1} 2^i (where i runs from 0 to n-1
91 */
92 template <typename VType>
binaryIndex(const GenericVector<VType> & v)93 Int binaryIndex(const GenericVector<VType>& v)
94 {
95 Int result = 0;
96 for (const Int i : indices(attach_selector(v.top(), operations::positive())))
97 result += pow(2,i);
98 return result;
99 }
100
101 /**
102 @brief Helper function for the refinement function.
103 Given a polyhedral cell in terms of rays and lineality space, it computes, whether a given ray
104 is contained in this cell (possibly modulo (1,..,1)).
105 @param Matrix<Rational> rays The rays of the cell
106 @param Matrix<Rational> lineality The lineality space of the cell
107 @param Vector<Rational> ray The ray to be tested
108 @param bool is_projective Whether coordinates are given as tropical projective coordinates.
109 (False means they're affine).
110 @param solver A convex hull solver
111 @returns true, if and only if ray lies in the cone
112 */
113 inline
is_ray_in_cone(const Matrix<Rational> & rays,const Matrix<Rational> & lineality,const Vector<Rational> & ray,const bool is_projective)114 bool is_ray_in_cone(const Matrix<Rational>& rays, const Matrix<Rational>& lineality,
115 const Vector<Rational>& ray, const bool is_projective)
116 {
117 const auto facets = is_projective ? enumerate_homogeneous_facets(rays, lineality)
118 : polytope::enumerate_facets(rays, lineality, false);
119 // Check equations
120 for (auto l = entire(rows(facets.second)); !l.at_end(); ++l) {
121 if (*l * ray != 0) return false;
122 }
123 // Check facets
124 for (auto f = entire(rows(facets.first)); !f.at_end(); ++f) {
125 if (*f * ray < 0) return false;
126 }
127 return true;
128 }
129
130 } }
131
132