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 functions to compute the star of a cycle at a point.
27 */
28
29 #include "polymake/client.h"
30 #include "polymake/Rational.h"
31 #include "polymake/Matrix.h"
32 #include "polymake/Vector.h"
33 #include "polymake/linalg.h"
34 #include "polymake/IncidenceMatrix.h"
35 #include "polymake/polytope/convex_hull.h"
36 #include "polymake/tropical/convex_hull_tools.h"
37 #include "polymake/tropical/thomog.h"
38 #include "polymake/tropical/star.h"
39
40 namespace polymake { namespace tropical {
41
42 using StarResult = std::pair< Matrix<Rational>, std::vector<Set<Int>>>;
43
44 template <typename Addition>
normalized_star_data(BigObject local_cycle,const Vector<Rational> & point)45 BigObject normalized_star_data(BigObject local_cycle, const Vector<Rational>& point)
46 {
47 const Matrix<Rational>& local_vertices = local_cycle.give("VERTICES");
48 const IncidenceMatrix<>& local_cones = local_cycle.give("MAXIMAL_POLYTOPES");
49 Matrix<Rational> lineality = local_cycle.give("LINEALITY_SPACE");
50 Vector<Integer> weights = local_cycle.give("WEIGHTS");
51 StarResult r = computeStar(point, local_vertices, local_cones);
52
53 // Find all zero rays and make the first coordinate a one
54 Set<Int> apices = indices( attach_selector( rows(r.first), operations::is_zero()));
55 r.first.col(0).slice(apices).fill(1);
56 Matrix<Rational> vertices(0, r.first.cols());
57 Vector<Int> vertex_map_vector = insert_rays(vertices, r.first, true);
58 Map<Int, Int> vertex_map;
59 for (auto vm = entire<indexed>(vertex_map_vector); !vm.at_end(); ++vm) {
60 vertex_map[ vm.index() ] = *vm;
61 }
62
63 // Map cones to new ray indices
64 RestrictedIncidenceMatrix<> polytopes;
65 for (auto cone : r.second) {
66 Set<Int> new_cone ( attach_operation( cone, pm::operations::associative_access<Map<Int, Int>, Int>(&vertex_map)));
67 polytopes /= new_cone;
68 }
69 IncidenceMatrix<> final_polytopes(std::move(polytopes));
70
71 // Now find redundant rays
72 Set<Int> used_rays;
73
74 Matrix<Rational> lineality_dehom = tdehomog(lineality);
75 for (auto p = entire(rows(final_polytopes)); !p.at_end(); ++p) {
76 const Set<Int> can_rays = polytope::get_non_redundant_points(vertices.minor(*p, All), lineality_dehom, true).first;
77 used_rays += select(*p, can_rays);
78 }
79
80 return BigObject("Cycle", mlist<Addition>(),
81 "VERTICES", thomog(vertices.minor(used_rays, All)),
82 "MAXIMAL_POLYTOPES", final_polytopes.minor(All, used_rays),
83 "LINEALITY_SPACE", lineality,
84 "WEIGHTS", weights);
85 }
86
87 template <typename Addition>
star_at_vertex(BigObject cycle,Int vertex_index)88 BigObject star_at_vertex(BigObject cycle, Int vertex_index)
89 {
90 BigObject local_cycle = call_function("local_vertex", cycle, vertex_index);
91 const Matrix<Rational>& vertices = cycle.give("VERTICES");
92 return normalized_star_data<Addition>(local_cycle, vertices.row(vertex_index));
93 }
94
95 template <typename Addition>
star_at_point(BigObject cycle,const Vector<Rational> & point)96 BigObject star_at_point(BigObject cycle, const Vector<Rational>& point)
97 {
98 BigObject local_cycle = call_function("local_point",cycle, point);
99 return normalized_star_data<Addition>(local_cycle, point);
100 }
101
102
103 UserFunctionTemplate4perl("# @category Local computations"
104 "# Computes the Star of a tropical cycle at one of its vertices."
105 "# @param Cycle<Addition> C a tropical cycle"
106 "# @param Int i The index of a vertex in [[VERTICES]], which should not be a ray"
107 "# @return Cycle<Addition> The Star of C at the vertex",
108 "star_at_vertex<Addition>(Cycle<Addition>,$)");
109
110 UserFunctionTemplate4perl("# @category Local computations"
111 "# Computes the Star of a tropical cycle at an arbitrary point in its support"
112 "# @param Cycle<Addition> C a tropical cycle "
113 "# @param Vector<Rational> v A point, given in tropical projective coordinates with"
114 "# leading coordinate and which should lie in the support of C"
115 "# @return Cycle<Addition> The Star of C at v (Note that the subdivision may be finer than"
116 "# a potential coarsest structure",
117 "star_at_point<Addition>(Cycle<Addition>,Vector<Rational>)");
118 } }
119