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/Matrix.h"
20 #include "polymake/Vector.h"
21 #include "polymake/Graph.h"
22 
23 namespace polymake { namespace polytope {
24 
25 template <typename Scalar>
26 Graph<Directed> dgraph(BigObject p, BigObject lp, OptionSet options)
27 {
28    Graph<Directed> DG=p.give("GRAPH.ADJACENCY");
29    const Matrix<Scalar> V=p.give("VERTICES");
30 
31    const bool inverse=options["inverse"], generic=options["generic"];
32    bool upper_bound=true, lower_bound=true, check_for_rays=false;
33    Vector<Scalar> obj;
34    Set<Int> rays;
35 
36    if (lp.lookup("LINEAR_OBJECTIVE") >> obj) {
37       obj=V*obj;
38       p.give("FAR_FACE") >> rays;
39       check_for_rays=!rays.empty();
40    } else {
41       lp.give("ABSTRACT_OBJECTIVE") >> obj;
42    }
43    if (inverse) obj.negate();
44 
45    auto obj_value=obj.begin();
46    for (auto n=entire(nodes(DG));  !n.at_end();  ++n, ++obj_value) {
47       if (check_for_rays && rays.contains(*n)) {
48          const Int ray_orientation = sign(*obj_value);
49          if (ray_orientation > 0) upper_bound = false;
50          if (ray_orientation < 0) lower_bound = false;
51          if (ray_orientation >= 0) n.out_edges().clear();
52          if (ray_orientation <= 0) n.in_edges().clear();
53       } else {
54          // affine vertex
55          for (auto e=n.out_edges().find_nearest(*n, operations::gt()); !e.at_end(); ) {
56             if (check_for_rays && rays.contains(e.to_node())) {
57                ++e;
58                continue;
59             }
60             Int idiff = sign(*obj_value - obj[e.to_node()]);
61             if (idiff == 0 && generic)
62                idiff = lex_compare(V.row(*n), V.row(e.to_node())) == (inverse ? pm::cmp_lt : pm::cmp_gt) ? 1 : -1;
63 
64             if (idiff <= 0)
65                n.in_edges().erase(e.to_node());
66             if (idiff >= 0)
67                n.out_edges().erase(e++);
68             else
69                ++e;
70          }
71       }
72    }
73 
74    if (!inverse && !generic) {
75       Set<Int> minface, maxface;
76       for (auto n=entire(nodes(DG));  !n.at_end();  ++n) {
77          if (!n.in_degree()) minface.push_back(n.index());
78          if (!n.out_degree()) maxface.push_back(n.index());
79       }
80       lp.take("MINIMAL_FACE") << minface;
81       lp.take("MAXIMAL_FACE")  << maxface;
82       if (lower_bound)
83          lp.take("MINIMAL_VALUE") << obj[minface.front()];
84       else
85          lp.take("MINIMAL_VALUE") << -std::numeric_limits<Scalar>::infinity();
86       if (upper_bound)
87          lp.take("MAXIMAL_VALUE") << obj[maxface.front()];
88       else
89          lp.take("MAXIMAL_VALUE") << std::numeric_limits<Scalar>::infinity();
90    }
91    return DG;
92 }
93 
94 template <typename Scalar>
95 Vector<Scalar> objective_values_for_embedding(BigObject p, BigObject lp)
96 {
97    const Matrix<Scalar> V = p.give("VERTICES");
98    const Vector<Scalar> Obj = lp.give("LINEAR_OBJECTIVE");
99    Vector<Scalar> val = V*Obj;
100    const Set<Int> rays = p.give("FAR_FACE");
101    if (!rays.empty()) {
102       const Scalar max_obj=accumulate(val.slice(~rays), operations::max()),
103          min_obj=accumulate(val.slice(~rays), operations::min());
104       for (auto r=entire(rays); !r.at_end(); ++r)
105          if (val[*r]>0)
106             val[*r]=2*max_obj-min_obj;
107          else
108             val[*r]=2*min_obj-max_obj;
109    }
110    return val;
111 }
112 
113 /* @category Optimization
114  Direct the graph of a polytope //P// according to a linear or abstract objective function.
115  The maximal and minimal values, which are attained by the objective function, as well
116  as the minimal and the maximal face are written into separate sections.
117  The option //inverse// directs the graph with decreasing instead of increasing edges.
118  If the option //generic// is set, ties will be broken by lexicographical ordering.
119  @param Polytope P
120  @param LinearProgram LP
121  @option Bool inverse inverts the direction
122  @option Bool generic breaks ties
123  @return Graph<Directed>
124 */
125 
126 FunctionTemplate4perl("dgraph<Scalar>(Polytope<Scalar>, LinearProgram<Scalar> { inverse => undef, generic => undef })");
127 
128 FunctionTemplate4perl("objective_values_for_embedding<Scalar>(Polytope<Scalar> LinearProgram<Scalar>)");
129 
130 } }
131 
132 // Local Variables:
133 // mode:C++
134 // c-basic-offset:3
135 // indent-tabs-mode:nil
136 // End:
137