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 <fcntl.h>
19 #include <stdio.h>
20
21 #include "polymake/polytope/lrs_interface.h"
22 #include "polymake/hash_set"
23 #include "polymake/ListMatrix.h"
24
25 #define MA
26 #define GMP
27 extern "C" {
28 #ifdef HAVE_LRSDRIVER
29 #include <lrsdriver.h>
30 #endif
31 #include <lrslib.h>
32 }
33 #undef positive
34 #undef negative
35 #undef zero
36 #undef one
37 #undef gcd
38 #undef copy
39 #undef sign
40 #undef GMP
41 #undef MA
42
43
44 namespace polymake { namespace polytope { namespace lrs_interface {
45
46 #ifdef POLYMAKE_LRS_STANDALONE_GLOBAL_INIT
47
48 namespace {
49
global_construct_standalone()50 void global_construct_standalone()
51 {
52 FILE* dummy_out = nullptr;
53 #ifdef POLYMAKE_LRS_SUPPRESS_OUTPUT
54 dummy_out = fopen("/dev/null", "w");
55 #endif
56 lrs_mp_init(0, nullptr, dummy_out);
57 }
58
global_destroy_standalone()59 void global_destroy_standalone()
60 {
61 #ifdef POLYMAKE_LRS_SUPPRESS_OUTPUT
62 fclose(lrs_ofp);
63 #endif
64 }
65
66 }
67
68 void (* const LrsInstance::Initializer::global_construct)() = &global_construct_standalone;
69 void (* const LrsInstance::Initializer::global_destroy)() = &global_destroy_standalone;
70
71 #endif
72
Initializer()73 LrsInstance::Initializer::Initializer()
74 {
75 global_construct();
76 }
77
~Initializer()78 LrsInstance::Initializer::~Initializer()
79 {
80 global_destroy();
81 }
82
83 class lrs_mp_vector_output {
84 public:
lrs_mp_vector_output(Int n)85 explicit lrs_mp_vector_output(Int n)
86 : d(n-1)
87 , ptr(lrs_alloc_mp_vector(d))
88 {
89 if (!ptr) throw std::bad_alloc();
90 }
91
~lrs_mp_vector_output()92 ~lrs_mp_vector_output() { lrs_clear_mp_vector(ptr, d); }
93
operator lrs_mp_vector()94 operator lrs_mp_vector() { return ptr; }
95
front() const96 const mpz_t& front() const { return ptr[0]; }
97
98 class iterator {
99 public:
100 typedef std::input_iterator_tag iterator_category;
101 typedef Rational value_type;
102 typedef Rational reference;
103 typedef void pointer;
104 typedef ptrdiff_t difference_type;
105
iterator(lrs_mp_vector_output & vec,bool or_arg)106 iterator(lrs_mp_vector_output& vec, bool or_arg)
107 : leading(vec.ptr)
108 , cur(leading)
109 , last(leading+vec.d)
110 , oriented(or_arg)
111 {}
112
operator *()113 Rational operator* ()
114 {
115 if (cur==leading) {
116 // looking for the leading non-zero
117 if (int sgn = mpz_sgn(*leading)) {
118 if (!oriented)
119 // all following elements will be divided through the leading non-zero
120 sgn = 1;
121 else if (sgn < 0)
122 // all following elements will be divided through abs(leading non-zero)
123 mpz_neg(*leading, *leading);
124 ++cur;
125 return Rational(sgn);
126 } else {
127 ++leading;
128 return Rational(std::move(*cur++));
129 }
130 } else if (cur < last) {
131 return Rational(std::move(*cur++), *leading);
132 } else {
133 return Rational(std::move(*cur++), std::move(*leading));
134 }
135 }
136
operator ++()137 iterator& operator++ () { return *this; }
138
139 private:
140 mpz_t* leading;
141 mpz_t* cur;
142 mpz_t* const last;
143 const bool oriented;
144 };
145
make_Vector(const bool oriented,const bool repair=true)146 Vector<Rational> make_Vector(const bool oriented, const bool repair = true)
147 {
148 Vector<Rational> result(d+1, iterator(*this, oriented));
149 if (repair) {
150 mpz_t* last=ptr+d;
151 if ((**last)._mp_alloc) --last;
152 for (mpz_t* cur=ptr; cur<=last; ++cur)
153 mpz_init(*cur);
154 }
155 return result;
156 }
157
158 private:
159 Int d;
160 lrs_mp_vector ptr;
161 };
162
163 // lrs stores numerators and denominators in separate integer vectors
164 class lrs_mp_vector_input {
165 public:
lrs_mp_vector_input(Int n_arg)166 explicit lrs_mp_vector_input(Int n_arg)
167 : n(n_arg)
168 , nums(new mpz_t[n])
169 , dens(new mpz_t[n]) {}
170
~lrs_mp_vector_input()171 ~lrs_mp_vector_input()
172 {
173 delete[] nums;
174 delete[] dens;
175 }
176
get_nums() const177 lrs_mp_vector get_nums() const { return nums; }
get_dens() const178 lrs_mp_vector get_dens() const { return dens; }
179
180 template <typename Iterator>
consume(Iterator && src)181 void consume(Iterator&& src)
182 {
183 for (Int i = 0; i < n; ++i, ++src) {
184 *nums[i] = *mpq_numref(src->get_rep());
185 *dens[i] = *mpq_denref(src->get_rep());
186 }
187 }
188
189 private:
190 Int n;
191 lrs_mp_vector nums;
192 lrs_mp_vector dens;
193 };
194
195 class lrs_mp_matrix_output {
196 public:
lrs_mp_matrix_output(lrs_mp_matrix A,Int m_arg,Int n_arg)197 lrs_mp_matrix_output(lrs_mp_matrix A, Int m_arg, Int n_arg)
198 : ptr(A)
199 , m(m_arg)
200 , n(n_arg) {}
201
~lrs_mp_matrix_output()202 ~lrs_mp_matrix_output()
203 {
204 if (ptr) lrs_clear_mp_matrix(ptr, m, n);
205 }
206
207 class iterator {
208 public:
209 typedef std::input_iterator_tag iterator_category;
210 typedef Rational value_type;
211 typedef Rational reference;
212 typedef void pointer;
213 typedef ptrdiff_t difference_type;
214
iterator(lrs_mp_matrix_output & mat)215 explicit iterator(lrs_mp_matrix_output& mat)
216 : vec(mat.ptr)
217 , i(0)
218 , n(mat.n)
219 {}
220
operator *()221 Rational operator* ()
222 {
223 return Rational(std::move((*vec)[i]));
224 }
225
operator ++()226 iterator& operator++ ()
227 {
228 if (++i==n) {
229 i=0;
230 ++vec;
231 }
232 return *this;
233 }
234
235 private:
236 lrs_mp_vector* vec;
237 Int i;
238 const Int n;
239 };
240
make_Matrix()241 Matrix<Rational> make_Matrix()
242 {
243 return Matrix<Rational>(m, n, iterator(*this));
244 }
245
246 private:
247 lrs_mp_matrix ptr;
248 Int m, n;
249 };
250
251 struct dictionary {
252 lrs_dat *Q;
253 lrs_dic_struct *P;
254 lrs_mp_matrix Lin;
255 FILE* save_lrs_ofp = nullptr;
256 #if defined(POLYMAKE_LRS_SUPPRESS_OUTPUT) && POLYMAKE_LRS_SUPPRESS_OUTPUT == 2
257 int save_stdout = -1;
258 #endif
259
260 // stream cleanup and restore stdout
restore_ofppolymake::polytope::lrs_interface::dictionary261 void restore_ofp() {
262 if (lrs_ofp == stderr) {
263 fflush(lrs_ofp);
264 lrs_ofp = save_lrs_ofp;
265 }
266 #if defined(POLYMAKE_LRS_SUPPRESS_OUTPUT) && POLYMAKE_LRS_SUPPRESS_OUTPUT == 2
267 else if (save_stdout != -1) {
268 if (stdout != nullptr)
269 fflush(stdout);
270 dup2(save_stdout, 1);
271 close(save_stdout);
272 }
273 #endif
274 }
275
~dictionarypolymake::polytope::lrs_interface::dictionary276 ~dictionary()
277 {
278 if (Lin) lrs_clear_mp_matrix(Lin, Q->nredundcol, Q->n);
279 lrs_free_dic(P,Q);
280 lrs_free_dat(Q);
281 restore_ofp();
282 }
283
284 // parameter ge: primal case: true for vertex, false for linearity
285 // dual case: true for inequality, false for equation
set_matrixpolymake::polytope::lrs_interface::dictionary286 void set_matrix(const Matrix<Rational>& A, Int start_row = 0, bool ge = true)
287 {
288 lrs_mp_vector_input vec(A.cols());
289 auto x=concat_rows(A).begin();
290
291 // lrs enumerates rows starting with 1
292 for (Int r = start_row+1, r_end = r+A.rows(); r != r_end; ++r) {
293 vec.consume(x);
294 lrs_set_row_mp(P, Q, r, vec.get_nums(), vec.get_dens(), ge);
295 }
296 }
297
set_obj_vectorpolymake::polytope::lrs_interface::dictionary298 void set_obj_vector(const Vector<Rational>& V, bool maximize)
299 {
300 const Int n = V.size();
301 if (n != Q->n)
302 throw std::runtime_error("lrs_interface - inequalities/objective dimension mismatch");
303 lrs_mp_vector_input vec(n);
304 vec.consume(V.begin());
305 lrs_set_obj_mp(P, Q, vec.get_nums(), vec.get_dens(), maximize);
306 Q->lponly=1;
307 }
308
309
dictionarypolymake::polytope::lrs_interface::dictionary310 dictionary(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations, const bool dual, const bool verbose=false)
311 {
312 // this avoids a segfault in some lrs versions
313 if (dual && Inequalities.cols() == 0 && Equations.cols() == 0)
314 throw std::runtime_error("lrs_interface - cannot handle ambient dimension 0 in dual mode");
315 // initialize static lrs data
316 Lin = nullptr;
317
318 if (verbose) {
319 save_lrs_ofp = lrs_ofp;
320 lrs_ofp = stderr;
321 }
322 #if defined(POLYMAKE_LRS_SUPPRESS_OUTPUT) && POLYMAKE_LRS_SUPPRESS_OUTPUT == 2
323 else {
324 save_stdout = dup(1);
325 dup2(fileno(lrs_ofp), 1);
326 }
327 #endif
328
329 char name[] = "polymake";
330 Q=lrs_alloc_dat(name);
331 if (!Q) {
332 restore_ofp();
333 throw std::bad_alloc();
334 }
335 if (verbose)
336 Q->debug=1;
337 Q->m=Inequalities.rows()+Equations.rows();
338 Q->n=Inequalities.cols();
339 if (!Q->n) Q->n=Equations.cols();
340 Q->hull=dual?0:1;
341
342 // initialize dynamic lrs data
343 P=lrs_alloc_dic(Q);
344 if (!P) {
345 restore_ofp();
346 lrs_free_dat(Q);
347 throw std::bad_alloc();
348 }
349 // store inequalities/points and equations/lineality in Q
350 if (Inequalities.rows()) set_matrix(Inequalities);
351 if (Equations.rows()) set_matrix(Equations, Inequalities.rows(), false);
352 }
353
354
355 // the following functions "run" lrs
356 enum filter_nothing_t { filter_nothing };
357 enum filter_bounded_t { filter_bounded };
358 enum filter_rays_t { filter_rays };
359 enum filter_facets_t { filter_facets };
360
get_solution_matrixpolymake::polytope::lrs_interface::dictionary361 Matrix<Rational> get_solution_matrix(filter_nothing_t)
362 {
363 ListMatrix<Vector<Rational>> facets(0,Q->n);
364
365 lrs_mp_vector_output output(Q->n);
366 do {
367 for (Int col = 0; col <= P->d; ++col)
368 if (lrs_getsolution(P, Q, output, col))
369 facets /= output.make_Vector(true);
370 } while (lrs_getnextbasis (&P, Q, 0));
371
372 return Matrix<Rational>(facets.rows(), facets.cols(), operations::move(), entire(rows(facets)));
373 }
374
count_solutionspolymake::polytope::lrs_interface::dictionary375 long count_solutions(filter_nothing_t)
376 {
377 long facets=0;
378
379 lrs_mp_vector_output output(Q->n);
380 do {
381 for (Int col = 0; col <= P->d; ++col)
382 if (lrs_getsolution(P, Q, output, col))
383 ++facets;
384 } while (lrs_getnextbasis (&P, Q, 0));
385
386 return facets;
387 }
388
get_solution_matrixpolymake::polytope::lrs_interface::dictionary389 Matrix<Rational> get_solution_matrix(filter_rays_t, bool isCone)
390 {
391 // each vertex is computed only once, but rays can appear multiple times.
392 ListMatrix<Vector<Rational>> vertices(0,Q->n);
393 hash_set<Vector<Rational>> rays;
394
395 lrs_mp_vector_output output(Q->n);
396 do {
397 for (Int col = 0; col <= P->d; ++col) {
398 if (lrs_getsolution(P, Q, output, col)) {
399 if (!mpz_sgn(output.front())) { // a ray starts with 0
400 rays.insert(output.make_Vector(true));
401 } else if (!isCone) {
402 // lrs returns the origin as a vertex for cones
403 // we have to remove this in our interpretation
404 vertices /= output.make_Vector(false);
405 }
406 }
407 }
408 } while (lrs_getnextbasis (&P, Q, 0));
409
410 if (isCone)
411 return Matrix<Rational>(rays.size(), Q->n, operations::move(), entire(rays));
412 else
413 return Matrix<Rational>(rays.size()+vertices.rows(), Q->n, operations::move(), entire(rays), entire(rows(vertices)));
414 }
415
get_linearitiespolymake::polytope::lrs_interface::dictionary416 Matrix<Rational> get_linearities()
417 {
418 const Int m = Q->nredundcol, n = Q->n;
419 lrs_mp_matrix_output output(Lin, m, n);
420 Lin=nullptr;
421 return output.make_Matrix();
422 }
423
count_solutionspolymake::polytope::lrs_interface::dictionary424 std::pair<long, long> count_solutions(filter_rays_t)
425 {
426 std::pair<long, long> vertices(0, 0);
427 hash_set<Vector<Rational>> rays;
428
429 lrs_mp_vector_output output(Q->n);
430 do {
431 for (Int col = 0; col <= P->d; ++col)
432 if (lrs_getsolution(P, Q, output, col)) {
433 if (mpz_sgn(output.front()))
434 ++vertices.second;
435 else
436 rays.insert(output.make_Vector(true));
437 }
438
439 } while (lrs_getnextbasis (&P, Q, 0));
440
441 vertices.first = vertices.second + rays.size();
442 return vertices;
443 }
444
count_solutionspolymake::polytope::lrs_interface::dictionary445 long count_solutions(filter_bounded_t)
446 {
447 long vertices=0;
448
449 lrs_mp_vector_output output(Q->n);
450 do {
451 for (Int col = 0; col <= P->d; ++col)
452 if (lrs_getsolution(P, Q, output, col)) {
453 if (mpz_sgn(output.front()))
454 ++vertices;
455 }
456 } while (lrs_getnextbasis (&P, Q, 0));
457
458 return vertices;
459 }
460
get_solution_matrixpolymake::polytope::lrs_interface::dictionary461 Matrix<Rational> get_solution_matrix(filter_facets_t)
462 {
463 hash_set<Vector<Rational>> facets(Q->m * Q->n);
464
465 lrs_mp_vector_output output(Q->n);
466 do {
467 for (Int col = 0; col <= P->d; ++col)
468 if (lrs_getsolution(P, Q, output, col))
469 facets.insert(output.make_Vector(true));
470 } while (lrs_getnextbasis (&P, Q, 0));
471
472 return Matrix<Rational>(facets.size(), Q->n, operations::move(), entire(facets));
473 }
474 };
475
476 convex_hull_result<Rational>
enumerate_facets(const Matrix<Rational> & Points,const Matrix<Rational> & Lineality,const bool isCone) const477 ConvexHullSolver::enumerate_facets(const Matrix<Rational>& Points, const Matrix<Rational>& Lineality, const bool isCone) const
478 {
479 dictionary D(Points, Lineality, false, verbose);
480 // we have a polytope if and only if all first coordinates are !=0
481 // FIXME find better name, vertex enumeration is the same for cones and polytopes
482 D.Q->polytope= isCone || attach_selector(Points.col(0), operations::is_zero()).empty();
483
484 if (!lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1) && !D.Q->nredundcol) throw infeasible();
485
486 Matrix<Rational> AH = isCone ? D.get_linearities().minor(range_from(1), All) : D.get_linearities(); // always lrs returns the functional [1,0,0,0,...]
487 Matrix<Rational> F = D.Q->polytope
488 ? D.get_solution_matrix(dictionary::filter_nothing) // lrs computes facets only once if input is a polytope
489 : D.get_solution_matrix(dictionary::filter_facets); // FIXME can facets appear several times for unbounded polyhedra?
490 // TODO: std::move
491 return { F, AH };
492 }
493
494 // FIXME check: why are unbounded polyhedra not allowed?
495 long
count_facets(const Matrix<Rational> & Points,const Matrix<Rational> & Lineality,const bool isCone) const496 ConvexHullSolver::count_facets(const Matrix<Rational>& Points, const Matrix<Rational>& Lineality, const bool isCone) const
497 {
498 dictionary D(Points, Lineality, verbose);
499
500 // CHECK: the restriction to bounded polyhedra has been added prior to the implementation
501 // of convex hull rules for unbounded polyhedra, and apparently the reason is not known anymore
502 // there is no documentation of this function in the lrs manual, so if we want to remove this, then someone has to check the code
503 if ( !isCone && !attach_selector(Points.col(0), operations::is_zero()).empty())
504 throw std::runtime_error("count_facets is not applicable to unbounded polyhedra");
505
506 if (!lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) throw infeasible();
507
508 return D.Q->nredundcol+1==D.Q->n
509 ? 0 // lrs does not treat the special case of a single point correctly
510 : D.count_solutions(dictionary::filter_nothing);
511 }
512
513 convex_hull_result<Rational>
enumerate_vertices(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations,const bool isCone) const514 ConvexHullSolver::enumerate_vertices(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations, const bool isCone) const
515 {
516 dictionary D(Inequalities, Equations, true, verbose);
517
518 if (!lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) throw infeasible();
519
520 Matrix<Rational> Lineality = D.get_linearities();
521 Matrix<Rational> Vertices = D.get_solution_matrix(dictionary::filter_rays, isCone);
522
523 // TODO: std::move
524 return { Vertices, Lineality };
525 }
526
527 ConvexHullSolver::vertex_count
count_vertices(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations,const bool only_bounded) const528 ConvexHullSolver::count_vertices(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations, const bool only_bounded) const
529 {
530 dictionary D(Inequalities, Equations, true, verbose);
531
532 if (!lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) throw infeasible();
533
534 vertex_count count;
535 count.lineality_dim = D.Q->nredundcol;
536 if (only_bounded) {
537 count.n_vertices = 0;
538 count.n_bounded_vertices = D.count_solutions(dictionary::filter_bounded);
539 } else {
540 std::tie(count.n_vertices, count.n_bounded_vertices) = D.count_solutions(dictionary::filter_rays);
541 }
542
543 return count;
544 }
545
546 //primal or dual
547 std::pair< Bitset, Matrix<Rational> >
find_irredundant_representation(const Matrix<Rational> & Points,const Matrix<Rational> & Lineality,const bool dual) const548 ConvexHullSolver::find_irredundant_representation(const Matrix<Rational>& Points, const Matrix<Rational>& Lineality, const bool dual) const
549 {
550 dictionary D(Points, Lineality, dual, verbose);
551
552 if (!lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) throw infeasible();
553 const Matrix<Rational> AH=D.get_linearities();
554
555 Bitset V(Points.rows());
556 for (Int index = D.Q->lastdv+1, end = D.P->m_A+D.P->d; index <= end; ++index)
557 if ( !checkindex(D.P,D.Q,index) )
558 V += D.Q->inequality[index - D.Q->lastdv]-1;
559
560 return std::pair< Bitset, Matrix<Rational> >(V,AH);
561 }
562
check_feasibility(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations) const563 bool LP_Solver::check_feasibility(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations) const
564 {
565 dictionary D(Inequalities, Equations,true);
566 return lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1);
567 }
568
check_feasibility(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations,Vector<Rational> & ValidPoint) const569 bool LP_Solver::check_feasibility(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations, Vector<Rational>& ValidPoint) const
570 {
571 dictionary D(Inequalities, Equations,true);
572
573 if (lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) {
574 lrs_mp_vector_output output(D.Q->n);
575 for (Int col = 0; col <= D.P->d; ++col)
576 if (lrs_getsolution(D.P, D.Q, output, col)) break;
577 ValidPoint = output.make_Vector(false, false);
578 return true;
579 }
580 return false;
581 }
582
583 LP_Solution<Rational>
solve(const Matrix<Rational> & Inequalities,const Matrix<Rational> & Equations,const Vector<Rational> & Objective,bool maximize,bool) const584 LP_Solver::solve(const Matrix<Rational>& Inequalities, const Matrix<Rational>& Equations,
585 const Vector<Rational>& Objective, bool maximize, bool) const
586 {
587 dictionary D(Inequalities, Equations, true);
588 D.set_obj_vector(Objective, maximize);
589
590 LP_Solution<Rational> result;
591 if (lrs_getfirstbasis(&D.P, D.Q, &D.Lin, 1)) {
592 result.lineality_dim = D.Q->nredundcol;
593 if (D.Q->unbounded) {
594 result.status = LP_status::unbounded;
595 } else {
596 result.status = LP_status::valid;
597
598 // sometimes there is lineality that fits the objective function but unbounded is not set
599 // hence we check all lineality rows that lrs computed manually
600 if (result.lineality_dim) {
601 Matrix<Rational> lin = D.get_linearities();
602 for (auto r = entire(rows(lin)); !r.at_end(); ++r) {
603 if (Objective * (*r) != 0) {
604 result.status = LP_status::unbounded;
605 break;
606 }
607 }
608 }
609
610 if (result.status == LP_status::valid) {
611 lrs_mp_vector_output output(D.Q->n);
612 for (Int col = 0; col <= D.P->d; ++col)
613 if (lrs_getsolution(D.P, D.Q, output, col)) break;
614
615 result.objective_value.set(std::move(D.P->objnum), std::move(D.P->objden));
616 result.solution = output.make_Vector(false, false);
617 }
618 }
619 } else {
620 result.status = LP_status::infeasible;
621 result.lineality_dim = 0;
622 }
623
624 return result;
625 }
626
627 } } }
628
629 // Local Variables:
630 // mode:C++
631 // c-basic-offset:3
632 // indent-tabs-mode:nil
633 // End:
634