1 /* PIP_Tree related class implementation: non-inline functions.
2 Copyright (C) 2001-2010 Roberto Bagnara <bagnara@cs.unipr.it>
3 Copyright (C) 2010-2016 BUGSENG srl (http://bugseng.com)
4
5 This file is part of the Parma Polyhedra Library (PPL).
6
7 The PPL is free software; you can redistribute it and/or modify it
8 under the terms of the GNU General Public License as published by the
9 Free Software Foundation; either version 3 of the License, or (at your
10 option) any later version.
11
12 The PPL is distributed in the hope that it will be useful, but WITHOUT
13 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
15 for more details.
16
17 You should have received a copy of the GNU General Public License
18 along with this program; if not, write to the Free Software Foundation,
19 Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02111-1307, USA.
20
21 For the most up-to-date information see the Parma Polyhedra Library
22 site: http://bugseng.com/products/ppl/ . */
23
24 #include "ppl-config.h"
25 #include "PIP_Tree_defs.hh"
26 #include "PIP_Problem_defs.hh"
27 #include <algorithm>
28 #include <memory>
29 #include <map>
30
31 #if 0
32 #define NOISY_PIP_TREE_STRUCTURE 1
33 #define NOISY_PIP 1
34 #define VERY_NOISY_PIP 1
35 #endif
36
37 namespace Parma_Polyhedra_Library {
38
39 namespace {
40
41 //! Assigns to \p x the positive remainder of the division of \p y by \p z.
42 inline void
pos_rem_assign(Coefficient & x,Coefficient_traits::const_reference y,Coefficient_traits::const_reference z)43 pos_rem_assign(Coefficient& x,
44 Coefficient_traits::const_reference y,
45 Coefficient_traits::const_reference z) {
46 rem_assign(x, y, z);
47 if (x < 0) {
48 x += z;
49 }
50 }
51
52 class Add_Mul_Assign_Row_Helper1 {
53 public:
Add_Mul_Assign_Row_Helper1(Coefficient_traits::const_reference coeff)54 Add_Mul_Assign_Row_Helper1(Coefficient_traits::const_reference coeff)
55 : c(coeff) {
56 }
57
58 void
operator ()(Coefficient & x,Coefficient_traits::const_reference y) const59 operator()(Coefficient& x, Coefficient_traits::const_reference y) const {
60 add_mul_assign(x, c, y);
61 }
62
63 private:
64 Coefficient c;
65 }; // class Add_Mul_Assign_Row_Helper1
66
67
68 class Add_Mul_Assign_Row_Helper2 {
69 public:
Add_Mul_Assign_Row_Helper2(Coefficient_traits::const_reference coeff)70 Add_Mul_Assign_Row_Helper2(Coefficient_traits::const_reference coeff)
71 : c(coeff) {
72 }
73
74 void
operator ()(Coefficient & x,Coefficient_traits::const_reference y) const75 operator()(Coefficient& x, Coefficient_traits::const_reference y) const {
76 x = y;
77 x *= c;
78 }
79
80 private:
81 Coefficient c;
82 }; // class Add_Mul_Assign_Row_Helper2
83
84 // Compute x += c * y
85 inline void
add_mul_assign_row(PIP_Tree_Node::Row & x,Coefficient_traits::const_reference c,const PIP_Tree_Node::Row & y)86 add_mul_assign_row(PIP_Tree_Node::Row& x,
87 Coefficient_traits::const_reference c,
88 const PIP_Tree_Node::Row& y) {
89 WEIGHT_BEGIN();
90 x.combine_needs_second(y,
91 Add_Mul_Assign_Row_Helper1(c),
92 Add_Mul_Assign_Row_Helper2(c));
93 WEIGHT_ADD_MUL(108, x.size());
94 }
95
96
97 struct Sub_Assign_Helper1 {
98 void
operator ()Parma_Polyhedra_Library::__anon8d897b600111::Sub_Assign_Helper199 operator()(Coefficient& x, Coefficient_traits::const_reference y) const {
100 x -= y;
101 }
102 }; // struct Sub_Assign_Helper1
103
104 struct Sub_Assign_Helper2 {
105 void
operator ()Parma_Polyhedra_Library::__anon8d897b600111::Sub_Assign_Helper2106 operator()(Coefficient& x, Coefficient_traits::const_reference y) const {
107 x = y;
108 neg_assign(x);
109 }
110 }; // struct Sub_Assign_Helper2
111
112 // Compute x -= y
113 inline void
sub_assign(PIP_Tree_Node::Row & x,const PIP_Tree_Node::Row & y)114 sub_assign(PIP_Tree_Node::Row& x, const PIP_Tree_Node::Row& y) {
115 WEIGHT_BEGIN();
116 x.combine_needs_second(y, Sub_Assign_Helper1(), Sub_Assign_Helper2());
117 WEIGHT_ADD_MUL(10, x.size());
118 }
119
120 // Merge constraint system to a matrix-form context such as x = x U y
121 void
merge_assign(Matrix<PIP_Tree_Node::Row> & x,const Constraint_System & y,const Variables_Set & parameters)122 merge_assign(Matrix<PIP_Tree_Node::Row>& x, const Constraint_System& y,
123 const Variables_Set& parameters) {
124 const dimension_type params_size = parameters.size();
125 PPL_ASSERT(params_size == x.num_columns() - 1);
126 const dimension_type new_rows = Implementation::num_constraints(y);
127 if (new_rows == 0) {
128 return;
129 }
130 const dimension_type old_num_rows = x.num_rows();
131 x.add_zero_rows(new_rows);
132
133 // Compute once for all.
134 const dimension_type cs_space_dim = y.space_dimension();
135 const Variables_Set::const_iterator param_begin = parameters.begin();
136 const Variables_Set::const_iterator param_end = parameters.end();
137
138 dimension_type i = old_num_rows;
139 for (Constraint_System::const_iterator y_i = y.begin(),
140 y_end = y.end(); y_i != y_end; ++y_i, ++i) {
141 WEIGHT_BEGIN();
142 PPL_ASSERT(y_i->is_nonstrict_inequality());
143 PIP_Tree_Node::Row& x_i = x[i];
144 Coefficient_traits::const_reference inhomogeneous_term
145 = y_i->inhomogeneous_term();
146 Variables_Set::const_iterator pj = param_begin;
147 PIP_Tree_Node::Row::iterator itr = x_i.end();
148 if (inhomogeneous_term != 0) {
149 itr = x_i.insert(0, inhomogeneous_term);
150 }
151 // itr may still be end() but it can still be used as a hint.
152 // TODO: This code could be optimized more (if it's expected that the
153 // size of `parameters' will be greater than the number of nonzero
154 // coefficients in y_i).
155 for (dimension_type j = 1; pj != param_end; ++pj, ++j) {
156 const Variable vj(*pj);
157 if (vj.space_dimension() > cs_space_dim) {
158 break;
159 }
160 Coefficient_traits::const_reference c = y_i->coefficient(vj);
161 if (c != 0) {
162 itr = x_i.insert(itr, j, c);
163 }
164 }
165 WEIGHT_ADD_MUL(102, params_size);
166 }
167 }
168
169 #if PPL_USE_SPARSE_MATRIX
170
171 // Assigns to row x the negation of row y.
172 inline void
neg_assign_row(PIP_Tree_Node::Row & x,const PIP_Tree_Node::Row & y)173 neg_assign_row(PIP_Tree_Node::Row& x, const PIP_Tree_Node::Row& y) {
174 WEIGHT_BEGIN();
175 x = y;
176 for (PIP_Tree_Node::Row::iterator i = x.begin(),
177 i_end = x.end(); i != i_end; ++i) {
178 neg_assign(*i);
179 WEIGHT_ADD(1);
180 }
181 }
182
183 #else // !PPL_USE_SPARSE_MATRIX
184
185 inline void
neg_assign_row(PIP_Tree_Node::Row & x,const PIP_Tree_Node::Row & y)186 neg_assign_row(PIP_Tree_Node::Row& x, const PIP_Tree_Node::Row& y) {
187 WEIGHT_BEGIN();
188 const dimension_type x_size = x.size();
189 for (dimension_type i = x_size; i-- > 0; )
190 neg_assign(x[i], y[i]);
191 WEIGHT_ADD_MUL(14, x_size);
192 }
193
194 #endif // !PPL_USE_SPARSE_MATRIX
195
196 // Given context row \p y and denominator \p denom,
197 // to be interpreted as expression expr = y / denom,
198 // assigns to context row \p x a new value such that
199 // x / denom == - expr - 1.
200 inline void
complement_assign(PIP_Tree_Node::Row & x,const PIP_Tree_Node::Row & y,Coefficient_traits::const_reference denom)201 complement_assign(PIP_Tree_Node::Row& x,
202 const PIP_Tree_Node::Row& y,
203 Coefficient_traits::const_reference denom) {
204 PPL_ASSERT(denom > 0);
205 neg_assign_row(x, y);
206 PIP_Tree_Node::Row::iterator itr = x.insert(0);
207 Coefficient& x_0 = *itr;
208 if (denom == 1) {
209 --x_0;
210 }
211 else {
212 PPL_DIRTY_TEMP_COEFFICIENT(mod);
213 pos_rem_assign(mod, x_0, denom);
214 x_0 -= (mod == 0) ? denom : mod;
215 }
216 if (x_0 == 0) {
217 x.reset(itr);
218 }
219 }
220
221 // Add to `context' the columns for new artificial parameters.
222 inline void
add_artificial_parameters(Matrix<PIP_Tree_Node::Row> & context,const dimension_type num_art_params)223 add_artificial_parameters(Matrix<PIP_Tree_Node::Row>& context,
224 const dimension_type num_art_params) {
225 if (num_art_params > 0) {
226 context.add_zero_columns(num_art_params);
227 }
228 }
229
230 // Add to `params' the indices of new artificial parameters.
231 inline void
add_artificial_parameters(Variables_Set & params,const dimension_type space_dim,const dimension_type num_art_params)232 add_artificial_parameters(Variables_Set& params,
233 const dimension_type space_dim,
234 const dimension_type num_art_params) {
235 for (dimension_type i = 0; i < num_art_params; ++i) {
236 params.insert(space_dim + i);
237 }
238 }
239
240 // Update `context', `params' and `space_dim' to account for
241 // the addition of the new artificial parameters.
242 inline void
add_artificial_parameters(Matrix<PIP_Tree_Node::Row> & context,Variables_Set & params,dimension_type & space_dim,const dimension_type num_art_params)243 add_artificial_parameters(Matrix<PIP_Tree_Node::Row>& context,
244 Variables_Set& params,
245 dimension_type& space_dim,
246 const dimension_type num_art_params) {
247 add_artificial_parameters(context, num_art_params);
248 add_artificial_parameters(params, space_dim, num_art_params);
249 space_dim += num_art_params;
250 }
251
252 /* Compares two columns lexicographically in a revised simplex tableau:
253 - returns true if
254 <CODE>
255 (column ja)*(-cst_a)/pivot_a[ja] < (column jb)*(-cst_b)/pivot_b[jb];
256 </CODE>
257 - returns false otherwise.
258 */
259 bool
column_lower(const Matrix<PIP_Tree_Node::Row> & tableau,const std::vector<dimension_type> & mapping,const std::vector<bool> & basis,const PIP_Tree_Node::Row & pivot_a,const dimension_type ja,const PIP_Tree_Node::Row & pivot_b,const dimension_type jb,Coefficient_traits::const_reference cst_a=-1,Coefficient_traits::const_reference cst_b=-1)260 column_lower(const Matrix<PIP_Tree_Node::Row>& tableau,
261 const std::vector<dimension_type>& mapping,
262 const std::vector<bool>& basis,
263 const PIP_Tree_Node::Row& pivot_a, const dimension_type ja,
264 const PIP_Tree_Node::Row& pivot_b, const dimension_type jb,
265 Coefficient_traits::const_reference cst_a = -1,
266 Coefficient_traits::const_reference cst_b = -1) {
267 Coefficient_traits::const_reference sij_a = pivot_a.get(ja);
268 Coefficient_traits::const_reference sij_b = pivot_b.get(jb);
269 PPL_ASSERT(sij_a > 0);
270 PPL_ASSERT(sij_b > 0);
271
272 PPL_DIRTY_TEMP_COEFFICIENT(lhs_coeff);
273 PPL_DIRTY_TEMP_COEFFICIENT(rhs_coeff);
274 lhs_coeff = cst_a * sij_b;
275 rhs_coeff = cst_b * sij_a;
276
277 const int lhs_coeff_sign = sgn(lhs_coeff);
278 const int rhs_coeff_sign = sgn(rhs_coeff);
279
280 if (ja == jb) {
281 // Same column: just compare the ratios.
282 // This works since all columns are lexico-positive.
283 return lhs_coeff > rhs_coeff;
284 }
285
286 PPL_DIRTY_TEMP_COEFFICIENT(lhs);
287 PPL_DIRTY_TEMP_COEFFICIENT(rhs);
288 const dimension_type num_vars = mapping.size();
289 dimension_type k = 0;
290 // While loop guard is: (k < num_rows && lhs == rhs).
291 // Return value is false, if k >= num_rows; it is equivalent to
292 // lhs < rhs, otherwise.
293 // Try to optimize the computation of lhs and rhs.
294 WEIGHT_BEGIN();
295 while (true) {
296 const dimension_type mk = mapping[k];
297 const bool in_base = basis[k];
298 if (++k >= num_vars) {
299 return false;
300 }
301 if (in_base) {
302 // Reconstitute the identity submatrix part of tableau.
303 if (mk == ja) {
304 // Optimizing for: lhs == lhs_coeff && rhs == 0;
305 if (lhs_coeff == 0) {
306 continue;
307 }
308 else {
309 return lhs_coeff > 0;
310 }
311 }
312 if (mk == jb) {
313 // Optimizing for: lhs == 0 && rhs == rhs_coeff;
314 if (rhs_coeff == 0) {
315 continue;
316 }
317 else {
318 return 0 > rhs_coeff;
319 }
320 }
321 // Optimizing for: lhs == 0 && rhs == 0;
322 continue;
323 }
324 else {
325 // Not in base.
326 const PIP_Tree_Node::Row& t_mk = tableau[mk];
327 Coefficient_traits::const_reference t_mk_ja = t_mk.get(ja);
328 Coefficient_traits::const_reference t_mk_jb = t_mk.get(jb);
329 if (t_mk_ja == 0) {
330 if (t_mk_jb == 0) {
331 continue;
332 }
333 else {
334 const int rhs_sign = rhs_coeff_sign * sgn(t_mk_jb);
335 if (0 == rhs_sign) {
336 continue;
337 }
338 else {
339 return 0 > rhs_sign;
340 }
341 }
342 }
343 else {
344 if (t_mk_jb == 0) {
345 const int lhs_sign = lhs_coeff_sign * sgn(t_mk_ja);
346 if (lhs_sign == 0) {
347 continue;
348 }
349 else {
350 return lhs_sign > 0;
351 }
352 }
353 else {
354 WEIGHT_ADD(2);
355 lhs = lhs_coeff * t_mk_ja;
356 rhs = rhs_coeff * t_mk_jb;
357 if (lhs == rhs) {
358 continue;
359 }
360 else {
361 return lhs > rhs;
362 }
363 }
364 }
365 }
366 }
367 // This point should be unreachable.
368 PPL_UNREACHABLE;
369 return false;
370 }
371
372 /* Find the column j in revised simplex tableau such that
373 - j is in candidates
374 - (column j) / pivot_row[j] is lexico-minimal
375 When this function returns, candidates contains the minimum(s) column(s)
376 index(es).
377 */
378 void
find_lexico_minimal_column_in_set(std::vector<dimension_type> & candidates,const Matrix<PIP_Tree_Node::Row> & tableau,const std::vector<dimension_type> & mapping,const std::vector<bool> & basis,const PIP_Tree_Node::Row & pivot_row)379 find_lexico_minimal_column_in_set(std::vector<dimension_type>& candidates,
380 const Matrix<PIP_Tree_Node::Row>& tableau,
381 const std::vector<dimension_type>& mapping,
382 const std::vector<bool>& basis,
383 const PIP_Tree_Node::Row& pivot_row) {
384 WEIGHT_BEGIN();
385 const dimension_type num_vars = mapping.size();
386
387 PPL_ASSERT(!candidates.empty());
388 // This is used as a set, it is always sorted.
389 std::vector<dimension_type> new_candidates;
390 for (dimension_type var_index = 0; var_index < num_vars; ++var_index) {
391 new_candidates.clear();
392 std::vector<dimension_type>::const_iterator i = candidates.begin();
393 std::vector<dimension_type>::const_iterator i_end = candidates.end();
394 PPL_ASSERT(!candidates.empty());
395 new_candidates.push_back(*i);
396 dimension_type min_column = *i;
397 ++i;
398 if (i == i_end) {
399 // Only one candidate left, so it is the minimum.
400 break;
401 }
402 PIP_Tree_Node::Row::const_iterator pivot_itr;
403 pivot_itr = pivot_row.find(min_column);
404 PPL_ASSERT(pivot_itr != pivot_row.end());
405 Coefficient sij_b = *pivot_itr;
406 ++pivot_itr;
407 const dimension_type row_index = mapping[var_index];
408 const bool in_base = basis[var_index];
409 if (in_base) {
410 for ( ; i != i_end; ++i) {
411 pivot_itr = pivot_row.find(pivot_itr, *i);
412 PPL_ASSERT(pivot_itr != pivot_row.end());
413 Coefficient_traits::const_reference sij_a = *pivot_itr;
414 ++pivot_itr;
415 PPL_ASSERT(sij_a > 0);
416 PPL_ASSERT(sij_b > 0);
417
418 // Reconstitute the identity submatrix part of tableau.
419 if (row_index != *i) {
420 if (row_index == min_column) {
421 // Optimizing for: lhs == 0 && rhs == rhs_coeff;
422 new_candidates.clear();
423 min_column = *i;
424 sij_b = sij_a;
425 new_candidates.push_back(min_column);
426 }
427 else {
428 // Optimizing for: lhs == 0 && rhs == 0;
429 new_candidates.push_back(*i);
430 }
431 }
432 }
433 WEIGHT_ADD_MUL(44, candidates.size());
434 }
435 else {
436 // Not in base.
437 const PIP_Tree_Node::Row& row = tableau[row_index];
438 PIP_Tree_Node::Row::const_iterator row_itr = row.lower_bound(min_column);
439 PIP_Tree_Node::Row::const_iterator row_end = row.end();
440 PPL_DIRTY_TEMP_COEFFICIENT(row_jb);
441 if (row_itr == row_end || row_itr.index() > min_column) {
442 row_jb = 0;
443 }
444 else {
445 PPL_ASSERT(row_itr.index() == min_column);
446 row_jb = *row_itr;
447 ++row_itr;
448 }
449 for ( ; i != i_end; ++i) {
450 pivot_itr = pivot_row.find(pivot_itr, *i);
451 PPL_ASSERT(pivot_itr != pivot_row.end());
452 Coefficient_traits::const_reference sij_a = *pivot_itr;
453 PPL_ASSERT(sij_a > 0);
454 PPL_ASSERT(sij_b > 0);
455
456 PPL_DIRTY_TEMP_COEFFICIENT(lhs);
457 PPL_DIRTY_TEMP_COEFFICIENT(rhs);
458 if (row_itr != row_end && row_itr.index() < *i) {
459 row_itr = row.lower_bound(row_itr, *i);
460 }
461 PPL_DIRTY_TEMP_COEFFICIENT(row_ja);
462 if (row_itr == row_end || row_itr.index() > *i) {
463 row_ja = 0;
464 }
465 else {
466 PPL_ASSERT(row_itr.index() == *i);
467 row_ja = *row_itr;
468 ++row_itr;
469 }
470
471 // lhs is actually the left-hand side with toggled sign.
472 // rhs is actually the right-hand side with toggled sign.
473 lhs = sij_b * row_ja;
474 rhs = sij_a * row_jb;
475 if (lhs == rhs) {
476 new_candidates.push_back(*i);
477 }
478 else {
479 if (lhs < rhs) {
480 new_candidates.clear();
481 min_column = *i;
482 row_jb = row_ja;
483 sij_b = sij_a;
484 new_candidates.push_back(min_column);
485 }
486 }
487 }
488 WEIGHT_ADD_MUL(68, candidates.size());
489 }
490 using std::swap;
491 swap(candidates, new_candidates);
492 }
493 }
494
495 /* Find the column j in revised simplex tableau such that
496 - pivot_row[j] is positive;
497 - (column j) / pivot_row[j] is lexico-minimal.
498 */
499 bool
find_lexico_minimal_column(const Matrix<PIP_Tree_Node::Row> & tableau,const std::vector<dimension_type> & mapping,const std::vector<bool> & basis,const PIP_Tree_Node::Row & pivot_row,const dimension_type start_j,dimension_type & j_out)500 find_lexico_minimal_column(const Matrix<PIP_Tree_Node::Row>& tableau,
501 const std::vector<dimension_type>& mapping,
502 const std::vector<bool>& basis,
503 const PIP_Tree_Node::Row& pivot_row,
504 const dimension_type start_j,
505 dimension_type& j_out) {
506 WEIGHT_BEGIN();
507 const dimension_type num_columns = tableau.num_columns();
508
509 PPL_ASSERT(start_j <= pivot_row.size());
510 if (start_j == pivot_row.size()) {
511 // There are no candidates, so there is no minimum.
512 return false;
513 }
514
515 // This is used as a set, it is always sorted.
516 std::vector<dimension_type> candidates;
517 for (PIP_Tree_Node::Row::const_iterator i = pivot_row.lower_bound(start_j),
518 i_end = pivot_row.end(); i != i_end; ++i) {
519 if (*i > 0) {
520 candidates.push_back(i.index());
521 }
522 }
523 WEIGHT_ADD_MUL(201, candidates.size());
524
525 if (candidates.empty()) {
526 j_out = num_columns;
527 return false;
528 }
529
530 find_lexico_minimal_column_in_set(candidates, tableau,
531 mapping, basis, pivot_row);
532 PPL_ASSERT(!candidates.empty());
533 j_out = *(candidates.begin());
534
535 return true;
536 }
537
538 // Computes into gcd the GCD of gcd and all coefficients in [first, last).
539 template <typename Iter>
540 void
gcd_assign_iter(Coefficient & gcd,Iter first,Iter last)541 gcd_assign_iter(Coefficient& gcd, Iter first, Iter last) {
542 PPL_ASSERT(gcd != 0);
543 if (gcd < 0) {
544 neg_assign(gcd);
545 }
546 if (gcd == 1) {
547 return;
548 }
549 WEIGHT_BEGIN();
550 for ( ; first != last; ++first) {
551 Coefficient_traits::const_reference coeff = *first;
552 if (coeff != 0) {
553 WEIGHT_ADD(24);
554 gcd_assign(gcd, coeff, gcd);
555 if (gcd == 1) {
556 return;
557 }
558 }
559 }
560 }
561
562 // Simplify row by exploiting variable integrality.
563 void
integral_simplification(PIP_Tree_Node::Row & row)564 integral_simplification(PIP_Tree_Node::Row& row) {
565 if (row[0] != 0) {
566 PIP_Tree_Node::Row::const_iterator j_begin = row.begin();
567 PIP_Tree_Node::Row::const_iterator j_end = row.end();
568 PPL_ASSERT(j_begin != j_end && j_begin.index() == 0 && *j_begin != 0);
569 /* Find next column with a non-zero value (there should be one). */
570 ++j_begin;
571 PPL_ASSERT(j_begin != j_end);
572 for ( ; *j_begin == 0; ++j_begin) {
573 PPL_ASSERT(j_begin != j_end);
574 }
575 /* Use it to initialize gcd. */
576 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
577 gcd = *j_begin;
578 ++j_begin;
579 gcd_assign_iter(gcd, j_begin, j_end);
580 if (gcd != 1) {
581 PPL_DIRTY_TEMP_COEFFICIENT(mod);
582 pos_rem_assign(mod, row[0], gcd);
583 row[0] -= mod;
584 }
585 }
586 /* Final normalization. */
587 row.normalize();
588 }
589
590 // Divide all coefficients in row x and denominator y by their GCD.
591 void
row_normalize(PIP_Tree_Node::Row & x,Coefficient & denom)592 row_normalize(PIP_Tree_Node::Row& x, Coefficient& denom) {
593 if (denom == 1) {
594 return;
595 }
596 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
597 gcd = denom;
598 gcd_assign_iter(gcd, x.begin(), x.end());
599
600 // Divide the coefficients by the GCD.
601 WEIGHT_BEGIN();
602 for (PIP_Tree_Node::Row::iterator i = x.begin(),
603 i_end = x.end(); i != i_end; ++i) {
604 Coefficient& x_i = *i;
605 exact_div_assign(x_i, x_i, gcd);
606 WEIGHT_ADD(22);
607 }
608 // Divide the denominator by the GCD.
609 exact_div_assign(denom, denom, gcd);
610 }
611
612 // This is here because it is used as a template argument in
613 // compatibility_check_find_pivot, so it must be a global declaration.
614 struct compatibility_check_find_pivot_in_set_data {
615 dimension_type row_index;
616 // We cache cost and value to avoid calling get() multiple times.
617 Coefficient cost;
618 Coefficient value;
operator ==Parma_Polyhedra_Library::__anon8d897b600111::compatibility_check_find_pivot_in_set_data619 bool operator==(const compatibility_check_find_pivot_in_set_data& x) const {
620 return row_index == x.row_index;
621 }
622 // Needed by std::vector to sort the values.
operator <Parma_Polyhedra_Library::__anon8d897b600111::compatibility_check_find_pivot_in_set_data623 bool operator<(const compatibility_check_find_pivot_in_set_data& x) const {
624 return row_index < x.row_index;
625 }
626 };
627
628 void
compatibility_check_find_pivot_in_set(std::vector<std::pair<dimension_type,compatibility_check_find_pivot_in_set_data>> & candidates,const Matrix<PIP_Tree_Node::Row> & s,const std::vector<dimension_type> & mapping,const std::vector<bool> & basis)629 compatibility_check_find_pivot_in_set(
630 std::vector<std::pair<dimension_type,
631 compatibility_check_find_pivot_in_set_data> >&
632 candidates,
633 const Matrix<PIP_Tree_Node::Row>& s,
634 const std::vector<dimension_type>& mapping,
635 const std::vector<bool>& basis) {
636
637 typedef compatibility_check_find_pivot_in_set_data data_struct;
638 typedef std::vector<std::pair<dimension_type, data_struct> > candidates_t;
639 // This is used as a set, it is always sorted.
640 candidates_t new_candidates;
641 const dimension_type num_vars = mapping.size();
642 for (dimension_type var_index = 0; var_index < num_vars; ++var_index) {
643 const dimension_type row_index = mapping[var_index];
644 const bool in_base = basis[var_index];
645 candidates_t::const_iterator i = candidates.begin();
646 candidates_t::const_iterator i_end = candidates.end();
647 PPL_ASSERT(i != i_end);
648 dimension_type pj = i->first;
649 Coefficient cost = i->second.cost;
650 Coefficient value = i->second.value;
651 new_candidates.clear();
652 new_candidates.push_back(*i);
653 if (in_base) {
654 for (++i; i != i_end; ++i) {
655 bool found_better_pivot = false;
656
657 const dimension_type challenger_j = i->first;
658 Coefficient_traits::const_reference challenger_cost = i->second.cost;
659 PPL_ASSERT(value > 0);
660 PPL_ASSERT(i->second.value > 0);
661 PPL_ASSERT(pj < challenger_j);
662
663 const int lhs_coeff_sgn = sgn(cost);
664 const int rhs_coeff_sgn = sgn(challenger_cost);
665
666 PPL_ASSERT(pj != challenger_j);
667
668 // Reconstitute the identity submatrix part of tableau.
669 if (row_index == pj) {
670 // Optimizing for: lhs == lhs_coeff && rhs == 0;
671 if (lhs_coeff_sgn == 0) {
672 new_candidates.push_back(*i);
673 }
674 else {
675 found_better_pivot = (lhs_coeff_sgn > 0);
676 }
677 }
678 else {
679 if (row_index == challenger_j) {
680 // Optimizing for: lhs == 0 && rhs == rhs_coeff;
681 if (rhs_coeff_sgn == 0) {
682 new_candidates.push_back(*i);
683 }
684 else {
685 found_better_pivot = (0 > rhs_coeff_sgn);
686 }
687 }
688 else {
689 // Optimizing for: lhs == 0 && rhs == 0;
690 new_candidates.push_back(*i);
691 }
692 }
693 if (found_better_pivot) {
694 pj = challenger_j;
695 cost = challenger_cost;
696 value = i->second.value;
697 new_candidates.clear();
698 new_candidates.push_back(*i);
699 }
700 }
701 }
702 else {
703 // Not in base.
704 const PIP_Tree_Node::Row& row = s[row_index];
705 PIP_Tree_Node::Row::const_iterator row_itr = row.lower_bound(pj);
706 PIP_Tree_Node::Row::const_iterator new_row_itr;
707 PIP_Tree_Node::Row::const_iterator row_end = row.end();
708 PPL_DIRTY_TEMP_COEFFICIENT(row_value);
709 if (row_itr != row_end && row_itr.index() == pj) {
710 row_value = *row_itr;
711 ++row_itr;
712 }
713 else {
714 row_value = 0;
715 }
716 PPL_DIRTY_TEMP_COEFFICIENT(row_challenger_value);
717 for (++i; i != i_end; ++i) {
718 const dimension_type challenger_j = i->first;
719 Coefficient_traits::const_reference challenger_cost = i->second.cost;
720 Coefficient_traits::const_reference challenger_value
721 = i->second.value;
722 PPL_ASSERT(value > 0);
723 PPL_ASSERT(challenger_value > 0);
724 PPL_ASSERT(pj < challenger_j);
725
726 new_row_itr = row.find(row_itr, challenger_j);
727 if (new_row_itr != row.end()) {
728 row_challenger_value = *new_row_itr;
729 // Use new_row_itr as a hint in next iterations
730 row_itr = new_row_itr;
731 }
732 else {
733 row_challenger_value = 0;
734 // Using end() as a hint is not useful, keep the current hint.
735 }
736 PPL_ASSERT(row_challenger_value == row.get(challenger_j));
737
738 // Before computing and comparing the actual values, the signs are
739 // compared. This speeds up the code, because the values' computation
740 // is a bit expensive.
741
742 const int lhs_sign = sgn(cost) * sgn(row_value);
743 const int rhs_sign = sgn(challenger_cost) * sgn(row_challenger_value);
744
745 if (lhs_sign != rhs_sign) {
746 if (lhs_sign > rhs_sign) {
747 #ifndef NDEBUG
748 pj = challenger_j;
749 #endif
750 cost = challenger_cost;
751 value = challenger_value;
752 row_value = row_challenger_value;
753 new_candidates.clear();
754 new_candidates.push_back(*i);
755 }
756 }
757 else {
758
759 // Sign comparison is not enough this time.
760 // Do the full computation.
761
762 PPL_DIRTY_TEMP_COEFFICIENT(lhs);
763 lhs = cost;
764 lhs *= challenger_value;
765 PPL_DIRTY_TEMP_COEFFICIENT(rhs);
766 rhs = challenger_cost;
767 rhs *= value;
768
769 lhs *= row_value;
770 rhs *= row_challenger_value;
771
772 if (lhs == rhs) {
773 new_candidates.push_back(*i);
774 }
775 else {
776 if (lhs > rhs) {
777 #ifndef NDEBUG
778 pj = challenger_j;
779 #endif
780 cost = challenger_cost;
781 value = challenger_value;
782 row_value = row_challenger_value;
783 new_candidates.clear();
784 new_candidates.push_back(*i);
785 }
786 }
787 }
788 }
789 }
790 using std::swap;
791 swap(candidates, new_candidates);
792 }
793 }
794
795 // Returns false if there is not a positive pivot candidate.
796 // Otherwise, it sets pi, pj to the coordinates of the pivot in s.
797 bool
compatibility_check_find_pivot(const Matrix<PIP_Tree_Node::Row> & s,const std::vector<dimension_type> & mapping,const std::vector<bool> & basis,dimension_type & pi,dimension_type & pj)798 compatibility_check_find_pivot(const Matrix<PIP_Tree_Node::Row>& s,
799 const std::vector<dimension_type>& mapping,
800 const std::vector<bool>& basis,
801 dimension_type& pi, dimension_type& pj) {
802 // Look for a negative RHS (i.e., constant term, stored in column 0),
803 // maximizing pivot column.
804 const dimension_type num_rows = s.num_rows();
805 typedef compatibility_check_find_pivot_in_set_data data_struct;
806 // This is used as a set, it is always sorted.
807 typedef std::vector<std::pair<dimension_type, data_struct> > candidates_t;
808 typedef std::map<dimension_type,data_struct> candidates_map_t;
809 candidates_map_t candidates_map;
810 for (dimension_type i = 0; i < num_rows; ++i) {
811 const PIP_Tree_Node::Row& s_i = s[i];
812 Coefficient_traits::const_reference s_i0 = s_i.get(0);
813 if (s_i0 < 0) {
814 dimension_type j;
815 if (!find_lexico_minimal_column(s, mapping, basis, s_i, 1, j)) {
816 // No positive pivot candidate: unfeasible problem.
817 return false;
818 }
819 Coefficient_traits::const_reference s_ij = s_i.get(j);
820 candidates_map_t::iterator itr = candidates_map.find(j);
821 if (itr == candidates_map.end()) {
822 data_struct& current_data = candidates_map[j];
823 current_data.row_index = i;
824 current_data.cost = s_i0;
825 current_data.value = s_ij;
826 }
827 else {
828 data_struct& current_data = candidates_map[j];
829 PPL_ASSERT(current_data.value > 0);
830
831 // Before computing and comparing the actual values, the signs are
832 // compared. This speeds up the code, because the values' computation
833 // is a bit expensive.
834 const int lhs_coeff_sgn = sgn(current_data.cost);
835 const int rhs_coeff_sgn = sgn(s_i0);
836
837 if (lhs_coeff_sgn != rhs_coeff_sgn) {
838 // Same column: just compare the ratios.
839 // This works since all columns are lexico-positive:
840 // return cst_a * sij_b > cst_b * sij_a.
841 if (lhs_coeff_sgn > rhs_coeff_sgn) {
842 // Found better pivot
843 current_data.row_index = i;
844 current_data.cost = s_i0;
845 current_data.value = s_ij;
846 }
847 // Otherwise, keep current pivot for this column.
848 }
849 else {
850 // Sign comparison is not enough this time.
851 // Do the full computation.
852 Coefficient_traits::const_reference value_b = s_i.get(j);
853 PPL_ASSERT(value_b > 0);
854
855 PPL_DIRTY_TEMP_COEFFICIENT(lhs_coeff);
856 lhs_coeff = current_data.cost;
857 lhs_coeff *= value_b;
858
859 PPL_DIRTY_TEMP_COEFFICIENT(rhs_coeff);
860 rhs_coeff = s_i0;
861 rhs_coeff *= current_data.value;
862
863 // Same column: just compare the ratios.
864 // This works since all columns are lexico-positive:
865 // return cst_a * sij_b > cst_b * sij_a.
866 if (lhs_coeff > rhs_coeff) {
867 // Found better pivot
868 current_data.row_index = i;
869 current_data.cost = s_i0;
870 current_data.value = s_ij;
871 }
872 // Otherwise, keep current pivot for this column.
873 }
874 }
875 }
876 }
877 candidates_t candidates;
878 for (candidates_map_t::iterator
879 i = candidates_map.begin(), i_end = candidates_map.end();
880 i != i_end; ++i) {
881 candidates.push_back(*i);
882 }
883 if (!candidates.empty()) {
884 compatibility_check_find_pivot_in_set(candidates, s, mapping, basis);
885 PPL_ASSERT(!candidates.empty());
886 pi = candidates.begin()->second.row_index;
887 pj = candidates.begin()->first;
888 }
889 else {
890 pi = s.num_rows();
891 pj = 0;
892 }
893
894 return true;
895 }
896
897 } // namespace
898
899 namespace IO_Operators {
900
901 std::ostream&
operator <<(std::ostream & os,const PIP_Tree_Node & x)902 operator<<(std::ostream& os, const PIP_Tree_Node& x) {
903 x.print(os);
904 return os;
905 }
906
907 std::ostream&
operator <<(std::ostream & os,const PIP_Tree_Node::Artificial_Parameter & x)908 operator<<(std::ostream& os, const PIP_Tree_Node::Artificial_Parameter& x) {
909 const Linear_Expression& expr = static_cast<const Linear_Expression&>(x);
910 os << "(" << expr << ") div " << x.denominator();
911 return os;
912 }
913
914 } // namespace IO_Operators
915
PIP_Tree_Node(const PIP_Problem * owner)916 PIP_Tree_Node::PIP_Tree_Node(const PIP_Problem* owner)
917 : owner_(owner),
918 parent_(0),
919 constraints_(),
920 artificial_parameters() {
921 }
922
PIP_Tree_Node(const PIP_Tree_Node & y)923 PIP_Tree_Node::PIP_Tree_Node(const PIP_Tree_Node& y)
924 : owner_(y.owner_),
925 parent_(0), // NOTE: parent is not copied.
926 constraints_(y.constraints_),
927 artificial_parameters(y.artificial_parameters) {
928 }
929
~PIP_Tree_Node()930 PIP_Tree_Node::~PIP_Tree_Node() {
931 }
932
933 PIP_Tree_Node::Artificial_Parameter
Artificial_Parameter(const Linear_Expression & expr,Coefficient_traits::const_reference d)934 ::Artificial_Parameter(const Linear_Expression& expr,
935 Coefficient_traits::const_reference d)
936 : Linear_Expression(expr), denom(d) {
937 if (denom == 0) {
938 throw std::invalid_argument("PIP_Tree_Node::Artificial_Parameter(e, d): "
939 "denominator d is zero.");
940 }
941
942 // Normalize if needed.
943 // FIXME: Provide a proper normalization helper.
944 Linear_Expression& param_expr = *this;
945 if (denom < 0) {
946 neg_assign(denom);
947 neg_assign(param_expr);
948 }
949
950 // Compute GCD of parameter expression and denominator.
951
952 Coefficient gcd = param_expr.gcd(0, space_dimension() + 1);
953
954 if (gcd == 1) {
955 return;
956 }
957
958 if (gcd == 0) {
959 gcd = denom;
960 }
961 else {
962 gcd_assign(gcd, denom, gcd);
963 }
964
965 if (gcd == 1) {
966 return;
967 }
968
969 // Divide coefficients and denominator by their (non-trivial) GCD.
970 PPL_ASSERT(gcd > 1);
971 param_expr.exact_div_assign(gcd, 0, space_dimension() + 1);
972 Parma_Polyhedra_Library::exact_div_assign(denom, denom, gcd);
973
974 PPL_ASSERT(OK());
975 }
976
977 bool
978 PIP_Tree_Node::Artificial_Parameter
operator ==(const PIP_Tree_Node::Artificial_Parameter & y) const979 ::operator==(const PIP_Tree_Node::Artificial_Parameter& y) const {
980 const Artificial_Parameter& x = *this;
981 if (x.space_dimension() != y.space_dimension()) {
982 return false;
983 }
984 if (x.denom != y.denom) {
985 return false;
986 }
987 return x.is_equal_to(y);
988 }
989
990 bool
991 PIP_Tree_Node::Artificial_Parameter
operator !=(const PIP_Tree_Node::Artificial_Parameter & y) const992 ::operator!=(const PIP_Tree_Node::Artificial_Parameter& y) const {
993 return !operator==(y);
994 }
995
996 bool
OK() const997 PIP_Tree_Node::Artificial_Parameter::OK() const {
998 if (denom <= 0) {
999 #ifndef NDEBUG
1000 std::cerr << "PIP_Tree_Node::Artificial_Parameter "
1001 << "has a non-positive denominator.\n";
1002 #endif
1003 return false;
1004 }
1005 return true;
1006 }
1007
1008 void
ascii_dump(std::ostream & s) const1009 PIP_Tree_Node::Artificial_Parameter::ascii_dump(std::ostream& s) const {
1010 s << "artificial_parameter ";
1011 Linear_Expression::ascii_dump(s);
1012 s << " / " << denom << "\n";
1013 }
1014
1015 bool
ascii_load(std::istream & s)1016 PIP_Tree_Node::Artificial_Parameter::ascii_load(std::istream& s) {
1017 std::string str;
1018 if (!(s >> str) || str != "artificial_parameter") {
1019 return false;
1020 }
1021 if (!Linear_Expression::ascii_load(s)) {
1022 return false;
1023 }
1024 if (!(s >> str) || str != "/") {
1025 return false;
1026 }
1027 if (!(s >> denom)) {
1028 return false;
1029 }
1030 PPL_ASSERT(OK());
1031 return true;
1032 }
1033
PPL_OUTPUT_DEFINITIONS(PIP_Tree_Node::Artificial_Parameter)1034 PPL_OUTPUT_DEFINITIONS(PIP_Tree_Node::Artificial_Parameter)
1035
1036 PIP_Solution_Node::PIP_Solution_Node(const PIP_Problem* owner)
1037 : PIP_Tree_Node(owner),
1038 tableau(),
1039 basis(),
1040 mapping(),
1041 var_row(),
1042 var_column(),
1043 special_equality_row(0),
1044 big_dimension(not_a_dimension()),
1045 sign(),
1046 solution(),
1047 solution_valid(false) {
1048 }
1049
PIP_Solution_Node(const PIP_Solution_Node & y)1050 PIP_Solution_Node::PIP_Solution_Node(const PIP_Solution_Node& y)
1051 : PIP_Tree_Node(y),
1052 tableau(y.tableau),
1053 basis(y.basis),
1054 mapping(y.mapping),
1055 var_row(y.var_row),
1056 var_column(y.var_column),
1057 special_equality_row(y.special_equality_row),
1058 big_dimension(y.big_dimension),
1059 sign(y.sign),
1060 solution(y.solution),
1061 solution_valid(y.solution_valid) {
1062 }
1063
PIP_Solution_Node(const PIP_Solution_Node & y,No_Constraints)1064 PIP_Solution_Node::PIP_Solution_Node(const PIP_Solution_Node& y,
1065 No_Constraints)
1066 : PIP_Tree_Node(y.owner_), // NOTE: only copy owner.
1067 tableau(y.tableau),
1068 basis(y.basis),
1069 mapping(y.mapping),
1070 var_row(y.var_row),
1071 var_column(y.var_column),
1072 special_equality_row(y.special_equality_row),
1073 big_dimension(y.big_dimension),
1074 sign(y.sign),
1075 solution(y.solution),
1076 solution_valid(y.solution_valid) {
1077 }
1078
~PIP_Solution_Node()1079 PIP_Solution_Node::~PIP_Solution_Node() {
1080 }
1081
PIP_Decision_Node(const PIP_Problem * owner,PIP_Tree_Node * fcp,PIP_Tree_Node * tcp)1082 PIP_Decision_Node::PIP_Decision_Node(const PIP_Problem* owner,
1083 PIP_Tree_Node* fcp,
1084 PIP_Tree_Node* tcp)
1085 : PIP_Tree_Node(owner),
1086 false_child(fcp),
1087 true_child(tcp) {
1088 if (false_child != 0) {
1089 false_child->set_parent(this);
1090 }
1091 if (true_child != 0) {
1092 true_child->set_parent(this);
1093 }
1094 }
1095
PIP_Decision_Node(const PIP_Decision_Node & y)1096 PIP_Decision_Node::PIP_Decision_Node(const PIP_Decision_Node& y)
1097 : PIP_Tree_Node(y),
1098 false_child(0),
1099 true_child(0) {
1100 if (y.false_child != 0) {
1101 false_child = y.false_child->clone();
1102 false_child->set_parent(this);
1103 }
1104 // Protect false_child from exception safety issues via std::auto_ptr.
1105 std::auto_ptr<PIP_Tree_Node> wrapped_node(false_child);
1106 if (y.true_child != 0) {
1107 true_child = y.true_child->clone();
1108 true_child->set_parent(this);
1109 }
1110 // It is now safe to release false_child.
1111 wrapped_node.release();
1112 }
1113
~PIP_Decision_Node()1114 PIP_Decision_Node::~PIP_Decision_Node() {
1115 delete false_child;
1116 delete true_child;
1117 }
1118
1119 void
set_owner(const PIP_Problem * owner)1120 PIP_Solution_Node::set_owner(const PIP_Problem* owner) {
1121 owner_ = owner;
1122 }
1123
1124 void
set_owner(const PIP_Problem * owner)1125 PIP_Decision_Node::set_owner(const PIP_Problem* owner) {
1126 owner_ = owner;
1127 if (false_child != 0) {
1128 false_child->set_owner(owner);
1129 }
1130 if (true_child != 0) {
1131 true_child->set_owner(owner);
1132 }
1133 }
1134
1135 bool
check_ownership(const PIP_Problem * owner) const1136 PIP_Solution_Node::check_ownership(const PIP_Problem* owner) const {
1137 return get_owner() == owner;
1138 }
1139
1140 bool
check_ownership(const PIP_Problem * owner) const1141 PIP_Decision_Node::check_ownership(const PIP_Problem* owner) const {
1142 return get_owner() == owner
1143 && (false_child == 0 || false_child->check_ownership(owner))
1144 && (true_child == 0 || true_child->check_ownership(owner));
1145 }
1146
1147 const PIP_Decision_Node*
as_decision() const1148 PIP_Decision_Node::as_decision() const {
1149 return this;
1150 }
1151
1152 const PIP_Decision_Node*
as_decision() const1153 PIP_Solution_Node::as_decision() const {
1154 return 0;
1155 }
1156
1157 const PIP_Solution_Node*
as_solution() const1158 PIP_Decision_Node::as_solution() const {
1159 return 0;
1160 }
1161
1162 const PIP_Solution_Node*
as_solution() const1163 PIP_Solution_Node::as_solution() const {
1164 return this;
1165 }
1166
1167 bool
OK() const1168 PIP_Solution_Node::Tableau::OK() const {
1169 if (s.num_rows() != t.num_rows()) {
1170 #ifndef NDEBUG
1171 std::cerr << "PIP_Solution_Node::Tableau matrices "
1172 << "have a different number of rows.\n";
1173 #endif
1174 return false;
1175 }
1176
1177 if (!s.OK() || !t.OK()) {
1178 #ifndef NDEBUG
1179 std::cerr << "A PIP_Solution_Node::Tableau matrix is broken.\n";
1180 #endif
1181 return false;
1182 }
1183
1184 if (denom <= 0) {
1185 #ifndef NDEBUG
1186 std::cerr << "PIP_Solution_Node::Tableau with non-positive "
1187 << "denominator.\n";
1188 #endif
1189 return false;
1190 }
1191
1192 // All tests passed.
1193 return true;
1194 }
1195
1196 bool
OK() const1197 PIP_Tree_Node::OK() const {
1198 #ifndef NDEBUG
1199 using std::endl;
1200 using std::cerr;
1201 #endif
1202
1203 // Parameter constraint system should contain no strict inequalities.
1204 for (Constraint_System::const_iterator
1205 i = constraints_.begin(), i_end = constraints_.end(); i != i_end; ++i) {
1206 if (i->is_strict_inequality()) {
1207 #ifndef NDEBUG
1208 cerr << "The feasible region of the PIP_Problem parameter context"
1209 << "is defined by a constraint system containing strict "
1210 << "inequalities."
1211 << endl;
1212 ascii_dump(cerr);
1213 #endif
1214 return false;
1215 }
1216 }
1217 return true;
1218 }
1219
1220 void
1221 PIP_Tree_Node
add_constraint(const Row & row,const Variables_Set & parameters)1222 ::add_constraint(const Row& row, const Variables_Set& parameters) {
1223 // Compute the expression for the parameter constraint.
1224 Linear_Expression expr = Linear_Expression(row.get(0));
1225 Variables_Set::const_iterator j = parameters.begin();
1226 if (!parameters.empty()) {
1227 // Needed to avoid reallocations in expr when iterating upward.
1228 add_mul_assign(expr, 0, Variable(*(parameters.rbegin())));
1229 // The number of increments of j plus one.
1230 dimension_type j_index = 1;
1231 Row::const_iterator i = row.begin();
1232 Row::const_iterator i_end = row.end();
1233 if (i != i_end && i.index() == 0) {
1234 ++i;
1235 }
1236 // NOTE: iterating in [1..num_params].
1237 WEIGHT_BEGIN();
1238 for ( ; i != i_end; ++i) {
1239 PPL_ASSERT(i.index() <= parameters.size());
1240 std::advance(j, i.index() - j_index);
1241 j_index = i.index();
1242 WEIGHT_ADD(74);
1243 add_mul_assign(expr, *i, Variable(*j));
1244 }
1245 }
1246 // Add the parameter constraint.
1247 constraints_.insert(expr >= 0);
1248 }
1249
1250 void
parent_merge()1251 PIP_Tree_Node::parent_merge() {
1252 const PIP_Decision_Node& parent = *parent_;
1253
1254 // Merge the parent's artificial parameters.
1255 artificial_parameters.insert(artificial_parameters.begin(),
1256 parent.art_parameter_begin(),
1257 parent.art_parameter_end());
1258
1259 PPL_ASSERT(OK());
1260 }
1261
1262 bool
OK() const1263 PIP_Solution_Node::OK() const {
1264 #ifndef NDEBUG
1265 using std::cerr;
1266 #endif
1267 if (!PIP_Tree_Node::OK()) {
1268 return false;
1269 }
1270
1271 // Check that every member used is OK.
1272
1273 if (!tableau.OK()) {
1274 return false;
1275 }
1276
1277 // Check coherency of basis, mapping, var_row and var_column
1278 if (basis.size() != mapping.size()) {
1279 #ifndef NDEBUG
1280 cerr << "The PIP_Solution_Node::basis and PIP_Solution_Node::mapping "
1281 << "vectors do not have the same number of elements.\n";
1282 #endif
1283 return false;
1284 }
1285 if (basis.size() != var_row.size() + var_column.size()) {
1286 #ifndef NDEBUG
1287 cerr << "The sum of number of elements in the PIP_Solution_Node::var_row "
1288 << "and PIP_Solution_Node::var_column vectors is different from the "
1289 << "number of elements in the PIP_Solution_Node::basis vector.\n";
1290 #endif
1291 return false;
1292 }
1293 if (var_column.size() != tableau.s.num_columns()) {
1294 #ifndef NDEBUG
1295 cerr << "The number of elements in the PIP_Solution_Node::var_column "
1296 << "vector is different from the number of columns in the "
1297 << "PIP_Solution_Node::tableau.s matrix.\n";
1298 #endif
1299 return false;
1300 }
1301 if (var_row.size() != tableau.s.num_rows()) {
1302 #ifndef NDEBUG
1303 cerr << "The number of elements in the PIP_Solution_Node::var_row "
1304 << "vector is different from the number of rows in the "
1305 << "PIP_Solution_Node::tableau.s matrix.\n";
1306 #endif
1307 return false;
1308 }
1309 for (dimension_type i = mapping.size(); i-- > 0; ) {
1310 const dimension_type row_column = mapping[i];
1311 if (basis[i] && var_column[row_column] != i) {
1312 #ifndef NDEBUG
1313 cerr << "Variable " << i << " is basic and corresponds to column "
1314 << row_column << " but PIP_Solution_Node::var_column["
1315 << row_column << "] does not correspond to variable " << i
1316 << ".\n";
1317 #endif
1318 return false;
1319 }
1320 if (!basis[i] && var_row[row_column] != i) {
1321 #ifndef NDEBUG
1322 cerr << "Variable " << i << " is nonbasic and corresponds to row "
1323 << row_column << " but PIP_Solution_Node::var_row["
1324 << row_column << "] does not correspond to variable " << i
1325 << ".\n";
1326 #endif
1327 return false;
1328 }
1329 }
1330 // All checks passed.
1331 return true;
1332 }
1333
1334 bool
OK() const1335 PIP_Decision_Node::OK() const {
1336 // Perform base class well-formedness check on this node.
1337 if (!PIP_Tree_Node::OK()) {
1338 return false;
1339 }
1340
1341 // Recursively check if child nodes are well-formed.
1342 if (false_child != 0 && !false_child->OK()) {
1343 return false;
1344 }
1345 if (true_child != 0 && !true_child->OK()) {
1346 return false;
1347 }
1348
1349 // Decision nodes should always have a true child.
1350 if (true_child == 0) {
1351 #ifndef NDEBUG
1352 std::cerr << "PIP_Decision_Node with no 'true' child.\n";
1353 #endif
1354 return false;
1355 }
1356
1357 // Decision nodes with a false child must have exactly one constraint.
1358 if (false_child != 0) {
1359 const dimension_type dist = Implementation::num_constraints(constraints_);
1360 if (dist != 1) {
1361 #ifndef NDEBUG
1362 std::cerr << "PIP_Decision_Node with a 'false' child has "
1363 << dist << " parametric constraints (should be 1).\n";
1364 #endif
1365 return false;
1366 }
1367 }
1368
1369 // All checks passed.
1370 return true;
1371 }
1372
1373 void
update_tableau(const PIP_Problem & pip,const dimension_type external_space_dim,const dimension_type first_pending_constraint,const Constraint_Sequence & input_cs,const Variables_Set & parameters)1374 PIP_Decision_Node::update_tableau(
1375 const PIP_Problem& pip,
1376 const dimension_type external_space_dim,
1377 const dimension_type first_pending_constraint,
1378 const Constraint_Sequence& input_cs,
1379 const Variables_Set& parameters) {
1380
1381 true_child->update_tableau(pip,
1382 external_space_dim,
1383 first_pending_constraint,
1384 input_cs,
1385 parameters);
1386 if (false_child != 0) {
1387 false_child->update_tableau(pip,
1388 external_space_dim,
1389 first_pending_constraint,
1390 input_cs,
1391 parameters);
1392 }
1393 PPL_ASSERT(OK());
1394 }
1395
1396 PIP_Tree_Node*
solve(const PIP_Problem & pip,const bool check_feasible_context,const Matrix<Row> & context,const Variables_Set & params,dimension_type space_dim,const int indent_level)1397 PIP_Decision_Node::solve(const PIP_Problem& pip,
1398 const bool check_feasible_context,
1399 const Matrix<Row>& context,
1400 const Variables_Set& params,
1401 dimension_type space_dim,
1402 const int indent_level) {
1403 PPL_ASSERT(indent_level >= 0);
1404 #ifdef NOISY_PIP_TREE_STRUCTURE
1405 indent_and_print(std::cerr, indent_level, "=== SOLVING DECISION NODE\n");
1406 #else
1407 PPL_USED(indent_level);
1408 #endif
1409 PPL_ASSERT(true_child != 0);
1410 Matrix<Row> context_true(context);
1411 Variables_Set all_params(params);
1412 const dimension_type num_art_params = artificial_parameters.size();
1413 add_artificial_parameters(context_true, all_params, space_dim,
1414 num_art_params);
1415 merge_assign(context_true, constraints_, all_params);
1416 const bool has_false_child = (false_child != 0);
1417 const bool has_true_child = (true_child != 0);
1418 #ifdef NOISY_PIP_TREE_STRUCTURE
1419 indent_and_print(std::cerr, indent_level,
1420 "=== DECISION: SOLVING THEN CHILD\n");
1421 #endif
1422 true_child = true_child->solve(pip, check_feasible_context,
1423 context_true, all_params, space_dim,
1424 indent_level + 1);
1425
1426 if (has_false_child) {
1427 // Decision nodes with false child must have exactly one constraint
1428 PPL_ASSERT(1 == Implementation::num_constraints(constraints_));
1429 // NOTE: modify context_true in place, complementing its last constraint.
1430 Matrix<Row>& context_false = context_true;
1431 Row& last = context_false[context_false.num_rows() - 1];
1432 complement_assign(last, last, 1);
1433 #ifdef NOISY_PIP_TREE_STRUCTURE
1434 indent_and_print(std::cerr, indent_level,
1435 "=== DECISION: SOLVING ELSE CHILD\n");
1436 #endif
1437 false_child = false_child->solve(pip, check_feasible_context,
1438 context_false, all_params, space_dim,
1439 indent_level + 1);
1440 }
1441
1442 if (true_child == 0 && false_child == 0) {
1443 // No children: the whole subtree is unfeasible.
1444 #ifdef NOISY_PIP_TREE_STRUCTURE
1445 indent_and_print(std::cerr, indent_level,
1446 "=== DECISION: BOTH BRANCHES NOW UNFEASIBLE: _|_\n");
1447 #endif
1448 delete this;
1449 return 0;
1450 }
1451
1452 if (has_false_child && false_child == 0) {
1453 // False child has become unfeasible: merge this node's artificials with
1454 // the true child, while removing the local parameter constraints, which
1455 // are no longer discriminative.
1456 #ifdef NOISY_PIP_TREE_STRUCTURE
1457 indent_and_print(std::cerr, indent_level,
1458 "=== DECISION: ELSE BRANCH NOW UNFEASIBLE\n");
1459 indent_and_print(std::cerr, indent_level,
1460 "==> merge then branch with parent.\n");
1461 #endif
1462 PIP_Tree_Node* const node = true_child;
1463 node->parent_merge();
1464 node->set_parent(parent());
1465 true_child = 0;
1466 delete this;
1467 PPL_ASSERT(node->OK());
1468 return node;
1469 }
1470 else if (has_true_child && true_child == 0) {
1471 // True child has become unfeasible: merge this node's artificials
1472 // with the false child.
1473 #ifdef NOISY_PIP_TREE_STRUCTURE
1474 indent_and_print(std::cerr, indent_level,
1475 "=== DECISION: THEN BRANCH NOW UNFEASIBLE\n");
1476 indent_and_print(std::cerr, indent_level,
1477 "==> merge else branch with parent.\n");
1478 #endif
1479 PIP_Tree_Node* const node = false_child;
1480 node->parent_merge();
1481 node->set_parent(parent());
1482 false_child = 0;
1483 delete this;
1484 PPL_ASSERT(node->OK());
1485 return node;
1486 }
1487 else if (check_feasible_context) {
1488 // Test all constraints for redundancy with the context, and eliminate
1489 // them if not necessary.
1490 Constraint_System cs;
1491 swap(cs, constraints_);
1492 for (Constraint_System::const_iterator ci = cs.begin(),
1493 ci_end = cs.end(); ci != ci_end; ++ci) {
1494 Matrix<Row> ctx_copy(context);
1495 merge_assign(ctx_copy, Constraint_System(*ci), all_params);
1496 Row& last = ctx_copy[ctx_copy.num_rows()-1];
1497 complement_assign(last, last, 1);
1498 if (compatibility_check(ctx_copy)) {
1499 // The constraint is not redundant with the context: keep it.
1500 constraints_.insert(*ci);
1501 }
1502 }
1503 // If the constraints set has become empty, only keep the true child.
1504 if (constraints_.empty()) {
1505 #ifdef NOISY_PIP_TREE_STRUCTURE
1506 indent_and_print(std::cerr, indent_level,
1507 "=== DECISION: NO BRANCHING CONSTRAINTS LEFT\n");
1508 indent_and_print(std::cerr, indent_level,
1509 "==> merge then branch with parent.\n");
1510 #endif
1511 PIP_Tree_Node* const node = true_child;
1512 node->parent_merge();
1513 node->set_parent(parent());
1514 true_child = 0;
1515 delete this;
1516 PPL_ASSERT(node->OK());
1517 return node;
1518 }
1519 }
1520 PPL_ASSERT(OK());
1521 return this;
1522 }
1523
1524 void
ascii_dump(std::ostream & s) const1525 PIP_Decision_Node::ascii_dump(std::ostream& s) const {
1526 // Dump base class info.
1527 PIP_Tree_Node::ascii_dump(s);
1528
1529 // Dump true child (if any).
1530 s << "\ntrue_child: ";
1531 if (true_child == 0) {
1532 // Note: this branch should normally be unreachable code, since a
1533 // well-formed decision node always has a true child. We keep this code
1534 // for debugging purposes (since we want to dump broken nodes).
1535 s << "BOTTOM\n";
1536 }
1537 else if (const PIP_Decision_Node* const dec = true_child->as_decision()) {
1538 s << "DECISION\n";
1539 dec->ascii_dump(s);
1540 }
1541 else {
1542 const PIP_Solution_Node* const sol = true_child->as_solution();
1543 PPL_ASSERT(sol != 0);
1544 s << "SOLUTION\n";
1545 sol->ascii_dump(s);
1546 }
1547
1548 // Dump false child (if any).
1549 s << "\nfalse_child: ";
1550 if (false_child == 0) {
1551 s << "BOTTOM\n";
1552 }
1553 else if (const PIP_Decision_Node* const dec = false_child->as_decision()) {
1554 // Note: this branch should normally be unreachable code.
1555 // Since a well-formed decision node having a false child should have
1556 // a single context constraint, its false child will have no context
1557 // constraints at all, so that no further branch is possible.
1558 // We keep this code for debugging purposes.
1559 s << "DECISION\n";
1560 dec->ascii_dump(s);
1561 }
1562 else {
1563 const PIP_Solution_Node* const sol = false_child->as_solution();
1564 PPL_ASSERT(sol != 0);
1565 s << "SOLUTION\n";
1566 sol->ascii_dump(s);
1567 }
1568 }
1569
1570 bool
ascii_load(std::istream & s)1571 PIP_Decision_Node::ascii_load(std::istream& s) {
1572 std::string str;
1573
1574 // Load base class info.
1575 if (!PIP_Tree_Node::ascii_load(s)) {
1576 return false;
1577 }
1578
1579 // Release the "true" subtree (if any).
1580 delete true_child;
1581 true_child = 0;
1582
1583 // Load true child (if any).
1584 if (!(s >> str) || str != "true_child:") {
1585 return false;
1586 }
1587 if (!(s >> str)) {
1588 return false;
1589 }
1590 if (str == "BOTTOM") {
1591 // Note: normally unreachable code (see comment on ascii_dump).
1592 true_child = 0;
1593 }
1594 else if (str == "DECISION") {
1595 PIP_Decision_Node* const dec = new PIP_Decision_Node(0, 0, 0);
1596 true_child = dec;
1597 if (!dec->ascii_load(s)) {
1598 return false;
1599 }
1600 }
1601 else if (str == "SOLUTION") {
1602 PIP_Solution_Node* const sol = new PIP_Solution_Node(0);
1603 true_child = sol;
1604 if (!sol->ascii_load(s)) {
1605 return false;
1606 }
1607 }
1608 else {
1609 // Unknown node kind.
1610 return false;
1611 }
1612
1613 // Release the "false" subtree (if any).
1614 delete false_child;
1615 false_child = 0;
1616
1617 // Load false child (if any).
1618 if (!(s >> str) || str != "false_child:") {
1619 return false;
1620 }
1621 if (!(s >> str)) {
1622 return false;
1623 }
1624 if (str == "BOTTOM") {
1625 false_child = 0;
1626 }
1627 else if (str == "DECISION") {
1628 // Note: normally unreachable code (see comment on ascii_dump).
1629 PIP_Decision_Node* const dec = new PIP_Decision_Node(0, 0, 0);
1630 false_child = dec;
1631 if (!dec->ascii_load(s)) {
1632 return false;
1633 }
1634 }
1635 else if (str == "SOLUTION") {
1636 PIP_Solution_Node* const sol = new PIP_Solution_Node(0);
1637 false_child = sol;
1638 if (!sol->ascii_load(s)) {
1639 return false;
1640 }
1641 }
1642 else {
1643 // Unknown node kind.
1644 return false;
1645 }
1646
1647 // Loaded all info.
1648 PPL_ASSERT(OK());
1649 return true;
1650 }
1651
1652
1653 void
normalize()1654 PIP_Solution_Node::Tableau::normalize() {
1655 if (denom == 1) {
1656 return;
1657 }
1658
1659 const dimension_type num_rows = s.num_rows();
1660
1661 // Compute global gcd.
1662 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
1663 gcd = denom;
1664 for (dimension_type i = num_rows; i-- > 0; ) {
1665 WEIGHT_BEGIN();
1666 const Row& s_i = s[i];
1667 for (Row::const_iterator
1668 j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
1669 Coefficient_traits::const_reference s_ij = *j;
1670 if (s_ij != 0) {
1671 WEIGHT_ADD(30);
1672 gcd_assign(gcd, s_ij, gcd);
1673 if (gcd == 1) {
1674 return;
1675 }
1676 }
1677 }
1678 WEIGHT_BEGIN();
1679 const Row& t_i = t[i];
1680 for (Row::const_iterator
1681 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
1682 Coefficient_traits::const_reference t_ij = *j;
1683 if (t_ij != 0) {
1684 WEIGHT_ADD(14);
1685 gcd_assign(gcd, t_ij, gcd);
1686 if (gcd == 1) {
1687 return;
1688 }
1689 }
1690 }
1691 }
1692 PPL_ASSERT(gcd > 1);
1693 // Normalize all coefficients.
1694 WEIGHT_BEGIN();
1695 for (dimension_type i = num_rows; i-- > 0; ) {
1696 Row& s_i = s[i];
1697 for (Row::iterator j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
1698 Coefficient& s_ij = *j;
1699 WEIGHT_ADD(19);
1700 exact_div_assign(s_ij, s_ij, gcd);
1701 }
1702 Row& t_i = t[i];
1703 for (Row::iterator j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
1704 Coefficient& t_ij = *j;
1705 WEIGHT_ADD(27);
1706 exact_div_assign(t_ij, t_ij, gcd);
1707 }
1708 }
1709 // Normalize denominator.
1710 exact_div_assign(denom, denom, gcd);
1711 }
1712
1713 void
scale(Coefficient_traits::const_reference ratio)1714 PIP_Solution_Node::Tableau::scale(Coefficient_traits::const_reference ratio) {
1715 WEIGHT_BEGIN();
1716 for (dimension_type i = s.num_rows(); i-- > 0; ) {
1717 Row& s_i = s[i];
1718 for (Row::iterator j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
1719 WEIGHT_ADD(19);
1720 *j *= ratio;
1721 }
1722 Row& t_i = t[i];
1723 for (Row::iterator j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
1724 WEIGHT_ADD(25);
1725 *j *= ratio;
1726 }
1727 }
1728 denom *= ratio;
1729 }
1730
1731 bool
1732 PIP_Solution_Node::Tableau
is_better_pivot(const std::vector<dimension_type> & mapping,const std::vector<bool> & basis,const dimension_type row_0,const dimension_type col_0,const dimension_type row_1,const dimension_type col_1) const1733 ::is_better_pivot(const std::vector<dimension_type>& mapping,
1734 const std::vector<bool>& basis,
1735 const dimension_type row_0,
1736 const dimension_type col_0,
1737 const dimension_type row_1,
1738 const dimension_type col_1) const {
1739 const dimension_type num_params = t.num_columns();
1740 const dimension_type num_rows = s.num_rows();
1741 const Row& s_0 = s[row_0];
1742 const Row& s_1 = s[row_1];
1743 Coefficient_traits::const_reference s_0_0 = s_0.get(col_0);
1744 Coefficient_traits::const_reference s_1_1 = s_1.get(col_1);
1745 const Row& t_0 = t[row_0];
1746 const Row& t_1 = t[row_1];
1747 PPL_DIRTY_TEMP_COEFFICIENT(product_0);
1748 PPL_DIRTY_TEMP_COEFFICIENT(product_1);
1749 WEIGHT_BEGIN();
1750 // On exit from the loop, if j_mismatch == num_params then
1751 // no column mismatch was found.
1752 dimension_type j_mismatch = num_params;
1753 Row::const_iterator j0 = t_0.end();
1754 Row::const_iterator j0_end = t_0.end();
1755 Row::const_iterator j1 = t_1.end();
1756 Row::const_iterator j1_end = t_1.end();
1757 for (dimension_type i = 0; i < num_rows; ++i) {
1758 const Row& s_i = s[i];
1759 Coefficient_traits::const_reference s_i_col_0 = s_i.get(col_0);
1760 Coefficient_traits::const_reference s_i_col_1 = s_i.get(col_1);
1761 j0 = t_0.begin();
1762 j1 = t_1.begin();
1763 while (j0 != j0_end && j1 != j1_end) {
1764 if (j0.index() == j1.index()) {
1765 WEIGHT_ADD(1361);
1766 product_0 = (*j0) * s_1_1 * s_i_col_0;
1767 product_1 = (*j1) * s_0_0 * s_i_col_1;
1768 if (product_0 != product_1) {
1769 // Mismatch found: exit from both loops.
1770 j_mismatch = j0.index();
1771 goto end_loop;
1772 }
1773 ++j0;
1774 ++j1;
1775 }
1776 else
1777 if (j0.index() < j1.index()) {
1778 if (*j0 != 0 && s_1_1 != 0 && s_i_col_0 != 0) {
1779 // Mismatch found: exit from both loops.
1780 j_mismatch = j0.index();
1781 goto end_loop;
1782 }
1783 ++j0;
1784 }
1785 else {
1786 PPL_ASSERT(j0.index() > j1.index());
1787 if (*j1 != 0 && s_0_0 != 0 && s_i_col_1 != 0) {
1788 // Mismatch found: exit from both loops.
1789 j_mismatch = j1.index();
1790 goto end_loop;
1791 }
1792 ++j1;
1793 }
1794 }
1795 while (j0 != j0_end) {
1796 if (*j0 != 0 && s_1_1 != 0 && s_i_col_0 != 0) {
1797 // Mismatch found: exit from both loops.
1798 j_mismatch = j0.index();
1799 goto end_loop;
1800 }
1801 ++j0;
1802 }
1803 while (j1 != j1_end) {
1804 if (*j1 != 0 && s_0_0 != 0 && s_i_col_1 != 0) {
1805 // Mismatch found: exit from both loops.
1806 j_mismatch = j1.index();
1807 goto end_loop;
1808 }
1809 }
1810 }
1811
1812 end_loop:
1813 return (j_mismatch != num_params)
1814 && column_lower(s, mapping, basis, s_0, col_0, s_1, col_1, *j0, *j1);
1815 }
1816
1817 void
ascii_dump(std::ostream & s) const1818 PIP_Tree_Node::ascii_dump(std::ostream& s) const {
1819 s << "constraints_\n";
1820 constraints_.ascii_dump(s);
1821 const dimension_type artificial_parameters_size
1822 = artificial_parameters.size();
1823 s << "\nartificial_parameters( " << artificial_parameters_size << " )\n";
1824 for (dimension_type i = 0; i < artificial_parameters_size; ++i) {
1825 artificial_parameters[i].ascii_dump(s);
1826 }
1827 }
1828
1829 bool
ascii_load(std::istream & s)1830 PIP_Tree_Node::ascii_load(std::istream& s) {
1831 std::string str;
1832 if (!(s >> str) || str != "constraints_") {
1833 return false;
1834 }
1835 constraints_.ascii_load(s);
1836
1837 if (!(s >> str) || str != "artificial_parameters(") {
1838 return false;
1839 }
1840
1841 dimension_type artificial_parameters_size;
1842 if (!(s >> artificial_parameters_size)) {
1843 return false;
1844 }
1845
1846 if (!(s >> str) || str != ")") {
1847 return false;
1848 }
1849 Artificial_Parameter ap;
1850 for (dimension_type i = 0; i < artificial_parameters_size; ++i) {
1851 if (!ap.ascii_load(s)) {
1852 return false;
1853 }
1854 artificial_parameters.push_back(ap);
1855 }
1856
1857 // Note: do not assert OK() here.
1858 // The node invariants should be checked on derived nodes.
1859 return true;
1860 }
1861
1862 PIP_Tree_Node*
clone() const1863 PIP_Solution_Node::clone() const {
1864 return new PIP_Solution_Node(*this);
1865 }
1866
1867 PIP_Tree_Node*
clone() const1868 PIP_Decision_Node::clone() const {
1869 return new PIP_Decision_Node(*this);
1870 }
1871
1872 void
ascii_dump(std::ostream & os) const1873 PIP_Solution_Node::Tableau::ascii_dump(std::ostream& os) const {
1874 os << "denominator " << denom << "\n";
1875 os << "variables ";
1876 s.ascii_dump(os);
1877 os << "parameters ";
1878 t.ascii_dump(os);
1879 }
1880
1881 bool
ascii_load(std::istream & is)1882 PIP_Solution_Node::Tableau::ascii_load(std::istream& is) {
1883 std::string str;
1884 if (!(is >> str) || str != "denominator") {
1885 return false;
1886 }
1887 Coefficient d;
1888 if (!(is >> d)) {
1889 return false;
1890 }
1891 denom = d;
1892 if (!(is >> str) || str != "variables") {
1893 return false;
1894 }
1895 if (!s.ascii_load(is)) {
1896 return false;
1897 }
1898 if (!(is >> str) || str != "parameters") {
1899 return false;
1900 }
1901 if (!t.ascii_load(is)) {
1902 return false;
1903 }
1904 PPL_ASSERT(OK());
1905 return true;
1906 }
1907
1908 void
ascii_dump(std::ostream & os) const1909 PIP_Solution_Node::ascii_dump(std::ostream& os) const {
1910 PIP_Tree_Node::ascii_dump(os);
1911
1912 os << "\ntableau\n";
1913 tableau.ascii_dump(os);
1914
1915 os << "\nbasis ";
1916 const dimension_type basis_size = basis.size();
1917 os << basis_size;
1918 for (dimension_type i = 0; i < basis_size; ++i) {
1919 os << (basis[i] ? " true" : " false");
1920 }
1921
1922 os << "\nmapping ";
1923 const dimension_type mapping_size = mapping.size();
1924 os << mapping_size;
1925 for (dimension_type i = 0; i < mapping_size; ++i) {
1926 os << " " << mapping[i];
1927 }
1928
1929 os << "\nvar_row ";
1930 const dimension_type var_row_size = var_row.size();
1931 os << var_row_size;
1932 for (dimension_type i = 0; i < var_row_size; ++i) {
1933 os << " " << var_row[i];
1934 }
1935
1936 os << "\nvar_column ";
1937 const dimension_type var_column_size = var_column.size();
1938 os << var_column_size;
1939 for (dimension_type i = 0; i < var_column_size; ++i) {
1940 os << " " << var_column[i];
1941 }
1942
1943 os << "\n";
1944
1945 os << "special_equality_row " << special_equality_row << "\n";
1946 os << "big_dimension " << big_dimension << "\n";
1947
1948 os << "sign ";
1949 const dimension_type sign_size = sign.size();
1950 os << sign_size;
1951 for (dimension_type i = 0; i < sign_size; ++i) {
1952 os << " ";
1953 switch (sign[i]) {
1954 case UNKNOWN:
1955 os << "UNKNOWN";
1956 break;
1957 case ZERO:
1958 os << "ZERO";
1959 break;
1960 case POSITIVE:
1961 os << "POSITIVE";
1962 break;
1963 case NEGATIVE:
1964 os << "NEGATIVE";
1965 break;
1966 case MIXED:
1967 os << "MIXED";
1968 break;
1969 }
1970 }
1971 os << "\n";
1972
1973 const dimension_type solution_size = solution.size();
1974 os << "solution " << solution_size << "\n";
1975 for (dimension_type i = 0; i < solution_size; ++i) {
1976 solution[i].ascii_dump(os);
1977 }
1978 os << "\n";
1979
1980 os << "solution_valid " << (solution_valid ? "true" : "false") << "\n";
1981 }
1982
1983 bool
ascii_load(std::istream & is)1984 PIP_Solution_Node::ascii_load(std::istream& is) {
1985 if (!PIP_Tree_Node::ascii_load(is)) {
1986 return false;
1987 }
1988
1989 std::string str;
1990 if (!(is >> str) || str != "tableau") {
1991 return false;
1992 }
1993 if (!tableau.ascii_load(is)) {
1994 return false;
1995 }
1996 if (!(is >> str) || str != "basis") {
1997 return false;
1998 }
1999 dimension_type basis_size;
2000 if (!(is >> basis_size)) {
2001 return false;
2002 }
2003 basis.clear();
2004 for (dimension_type i = 0; i < basis_size; ++i) {
2005 if (!(is >> str)) {
2006 return false;
2007 }
2008 bool val = false;
2009 if (str == "true") {
2010 val = true;
2011 }
2012 else if (str != "false") {
2013 return false;
2014 }
2015 basis.push_back(val);
2016 }
2017
2018 if (!(is >> str) || str != "mapping") {
2019 return false;
2020 }
2021 dimension_type mapping_size;
2022 if (!(is >> mapping_size)) {
2023 return false;
2024 }
2025 mapping.clear();
2026 for (dimension_type i = 0; i < mapping_size; ++i) {
2027 dimension_type val;
2028 if (!(is >> val)) {
2029 return false;
2030 }
2031 mapping.push_back(val);
2032 }
2033
2034 if (!(is >> str) || str != "var_row") {
2035 return false;
2036 }
2037 dimension_type var_row_size;
2038 if (!(is >> var_row_size)) {
2039 return false;
2040 }
2041 var_row.clear();
2042 for (dimension_type i = 0; i < var_row_size; ++i) {
2043 dimension_type val;
2044 if (!(is >> val)) {
2045 return false;
2046 }
2047 var_row.push_back(val);
2048 }
2049
2050 if (!(is >> str) || str != "var_column") {
2051 return false;
2052 }
2053 dimension_type var_column_size;
2054 if (!(is >> var_column_size)) {
2055 return false;
2056 }
2057 var_column.clear();
2058 for (dimension_type i = 0; i < var_column_size; ++i) {
2059 dimension_type val;
2060 if (!(is >> val)) {
2061 return false;
2062 }
2063 var_column.push_back(val);
2064 }
2065
2066 if (!(is >> str) || str != "special_equality_row") {
2067 return false;
2068 }
2069 if (!(is >> special_equality_row)) {
2070 return false;
2071 }
2072
2073 if (!(is >> str) || str != "big_dimension") {
2074 return false;
2075 }
2076 if (!(is >> big_dimension)) {
2077 return false;
2078 }
2079
2080 if (!(is >> str) || str != "sign") {
2081 return false;
2082 }
2083 dimension_type sign_size;
2084 if (!(is >> sign_size)) {
2085 return false;
2086 }
2087
2088 sign.clear();
2089 for (dimension_type i = 0; i < sign_size; ++i) {
2090 if (!(is >> str)) {
2091 return false;
2092 }
2093 Row_Sign val;
2094 if (str == "UNKNOWN") {
2095 val = UNKNOWN;
2096 }
2097 else if (str == "ZERO") {
2098 val = ZERO;
2099 }
2100 else if (str == "POSITIVE") {
2101 val = POSITIVE;
2102 }
2103 else if (str == "NEGATIVE") {
2104 val = NEGATIVE;
2105 }
2106 else if (str == "MIXED") {
2107 val = MIXED;
2108 }
2109 else {
2110 return false;
2111 }
2112 sign.push_back(val);
2113 }
2114
2115 if (!(is >> str) || str != "solution") {
2116 return false;
2117 }
2118 dimension_type solution_size;
2119 if (!(is >> solution_size)) {
2120 return false;
2121 }
2122 solution.clear();
2123 for (dimension_type i = 0; i < solution_size; ++i) {
2124 Linear_Expression val;
2125 if (!val.ascii_load(is)) {
2126 return false;
2127 }
2128 solution.push_back(val);
2129 }
2130
2131 if (!(is >> str) || str != "solution_valid") {
2132 return false;
2133 }
2134 if (!(is >> str)) {
2135 return false;
2136 }
2137 if (str == "true") {
2138 solution_valid = true;
2139 }
2140 else if (str == "false") {
2141 solution_valid = false;
2142 }
2143 else {
2144 return false;
2145 }
2146
2147 PPL_ASSERT(OK());
2148 return true;
2149 }
2150
2151 PIP_Solution_Node::Row_Sign
row_sign(const Row & x,const dimension_type big_dimension)2152 PIP_Solution_Node::row_sign(const Row& x,
2153 const dimension_type big_dimension) {
2154 if (big_dimension != not_a_dimension()) {
2155 // If a big parameter has been set and its coefficient is not zero,
2156 // then return the sign of the coefficient.
2157 Coefficient_traits::const_reference x_big = x.get(big_dimension);
2158 if (x_big > 0) {
2159 return POSITIVE;
2160 }
2161 if (x_big < 0) {
2162 return NEGATIVE;
2163 }
2164 // Otherwise x_big == 0, then no big parameter involved.
2165 }
2166
2167 PIP_Solution_Node::Row_Sign sign = ZERO;
2168 for (Row::const_iterator i = x.begin(), i_end = x.end(); i != i_end; ++i) {
2169 Coefficient_traits::const_reference x_i = *i;
2170 if (x_i > 0) {
2171 if (sign == NEGATIVE) {
2172 return MIXED;
2173 }
2174 sign = POSITIVE;
2175 }
2176 else if (x_i < 0) {
2177 if (sign == POSITIVE) {
2178 return MIXED;
2179 }
2180 sign = NEGATIVE;
2181 }
2182 }
2183 return sign;
2184 }
2185
2186 bool
compatibility_check(const Matrix<Row> & context,const Row & row)2187 PIP_Tree_Node::compatibility_check(const Matrix<Row>& context, const Row& row) {
2188 // CHECKME: do `context' and `row' have compatible (row) capacity?
2189 Matrix<Row> s(context);
2190 s.add_row(row);
2191 return compatibility_check(s);
2192 }
2193
2194 bool
compatibility_check(Matrix<Row> & s)2195 PIP_Tree_Node::compatibility_check(Matrix<Row>& s) {
2196 PPL_ASSERT(s.OK());
2197 // Note: num_rows may increase.
2198 dimension_type num_rows = s.num_rows();
2199 const dimension_type num_columns = s.num_columns();
2200 const dimension_type num_vars = num_columns - 1;
2201
2202 std::vector<Coefficient> scaling(num_rows, 1);
2203 std::vector<bool> basis;
2204 basis.reserve(num_vars + num_rows);
2205 std::vector<dimension_type> mapping;
2206 mapping.reserve(num_vars + num_rows);
2207 std::vector<dimension_type> var_row;
2208 var_row.reserve(num_rows);
2209 std::vector<dimension_type> var_column;
2210 var_column.reserve(num_columns);
2211
2212 // Column 0 is the constant term, not a variable
2213 var_column.push_back(not_a_dimension());
2214 for (dimension_type j = 1; j <= num_vars; ++j) {
2215 basis.push_back(true);
2216 mapping.push_back(j);
2217 var_column.push_back(j - 1);
2218 }
2219 for (dimension_type i = 0; i < num_rows; ++i) {
2220 basis.push_back(false);
2221 mapping.push_back(i);
2222 var_row.push_back(i + num_vars);
2223 }
2224
2225 // Scaling factor (i.e., denominator) for pivot coefficients.
2226 PPL_DIRTY_TEMP_COEFFICIENT(pivot_denom);
2227 // Allocate once and for all: short life temporaries.
2228 PPL_DIRTY_TEMP_COEFFICIENT(product);
2229 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
2230 PPL_DIRTY_TEMP_COEFFICIENT(scale_factor);
2231
2232 // Perform simplex pivots on the context
2233 // until we find an empty solution or an optimum.
2234 while (true) {
2235 // Check if the client has requested abandoning all expensive
2236 // computations. If so, the exception specified by the client
2237 // is thrown now.
2238 maybe_abandon();
2239
2240 // pi is the pivot's row index.
2241 dimension_type pi = num_rows;
2242 // pj is the pivot's column index.
2243 dimension_type pj = 0;
2244
2245 const bool found_positive_pivot_candidate
2246 = compatibility_check_find_pivot(s, mapping, basis, pi, pj);
2247
2248 if (!found_positive_pivot_candidate) {
2249 return false;
2250 }
2251
2252 if (pj == 0) {
2253 // No negative RHS: fractional optimum found.
2254 // If it is integer, then the test is successful.
2255 // Otherwise, generate a new cut.
2256 bool all_integer_vars = true;
2257 // NOTE: iterating downwards would be correct, but it would change
2258 // the ordering of cut generation.
2259 WEIGHT_BEGIN();
2260 for (dimension_type i = 0; i < num_vars; ++i) {
2261 if (basis[i]) {
2262 // Basic variable = 0, hence integer.
2263 continue;
2264 }
2265 // Not a basic variable.
2266 WEIGHT_ADD(70);
2267 const dimension_type mi = mapping[i];
2268 Coefficient_traits::const_reference denom = scaling[mi];
2269 if (s[mi].get(0) % denom == 0) {
2270 continue;
2271 }
2272 // Here constant term is not integer.
2273 all_integer_vars = false;
2274 // Generate a new cut.
2275 var_row.push_back(mapping.size());
2276 basis.push_back(false);
2277 mapping.push_back(num_rows);
2278 s.add_zero_rows(1);
2279 Row& cut = s[num_rows];
2280 ++num_rows;
2281 const Row& s_mi = s[mi];
2282 cut = s_mi;
2283 for (Row::iterator
2284 j = cut.begin(), j_end = cut.end(); j != j_end; ++j) {
2285 WEIGHT_ADD(32);
2286 pos_rem_assign(*j, *j, denom);
2287 }
2288 cut[0] -= denom;
2289 scaling.push_back(denom);
2290 }
2291 // Check if an integer solution was found.
2292 if (all_integer_vars) {
2293 return true;
2294 }
2295 else {
2296 continue;
2297 }
2298 }
2299
2300 // Here we have a positive s[pi][pj] pivot.
2301
2302 // Normalize the tableau before pivoting.
2303 for (dimension_type i = num_rows; i-- > 0; ) {
2304 row_normalize(s[i], scaling[i]);
2305 }
2306
2307 // Update basis.
2308 {
2309 const dimension_type var_pi = var_row[pi];
2310 const dimension_type var_pj = var_column[pj];
2311 var_row[pi] = var_pj;
2312 var_column[pj] = var_pi;
2313 basis[var_pi] = true;
2314 basis[var_pj] = false;
2315 mapping[var_pi] = pj;
2316 mapping[var_pj] = pi;
2317 }
2318
2319 // Create an identity row corresponding to basic variable pj.
2320 s.add_zero_rows(1);
2321 Row& pivot = s[num_rows];
2322 pivot[pj] = 1;
2323
2324 // Swap identity row with the pivot row previously found.
2325 using std::swap;
2326 swap(pivot, s[pi]);
2327 // Save original pivot scaling factor in a temporary,
2328 // then reset scaling factor for identity row.
2329 pivot_denom = scaling[pi];
2330 scaling[pi] = 1;
2331
2332 // Perform a pivot operation on the matrix.
2333 Coefficient_traits::const_reference pivot_pj = pivot.get(pj);
2334 {
2335 for (Row::const_iterator
2336 j = pivot.begin(), j_end = pivot.end(); j != j_end; ++j) {
2337 if (j.index() == pj) {
2338 continue;
2339 }
2340 Coefficient_traits::const_reference pivot_j = *j;
2341 // Do nothing if the j-th pivot element is zero.
2342 if (pivot_j == 0) {
2343 continue;
2344 }
2345
2346 WEIGHT_BEGIN();
2347 for (dimension_type i = num_rows; i-- > 0; ) {
2348 Row& s_i = s[i];
2349 product = s_i.get(pj) * pivot_j;
2350 if (product % pivot_pj != 0) {
2351 WEIGHT_ADD(103);
2352 // Must scale row s_i to stay in integer case.
2353 gcd_assign(gcd, product, pivot_pj);
2354 exact_div_assign(scale_factor, pivot_pj, gcd);
2355 for (Row::iterator
2356 k = s_i.begin(), k_end = s_i.end(); k != k_end; ++k) {
2357 WEIGHT_ADD(30);
2358 *k *= scale_factor;
2359 }
2360 product *= scale_factor;
2361 scaling[i] *= scale_factor;
2362 }
2363 PPL_ASSERT(product % pivot_pj == 0);
2364 exact_div_assign(product, product, pivot_pj);
2365 s_i[j.index()] -= product;
2366 WEIGHT_ADD(134);
2367 }
2368 }
2369 }
2370 // Update column only if pivot coordinate != 1.
2371 if (pivot_pj != pivot_denom) {
2372 WEIGHT_BEGIN();
2373 for (dimension_type i = num_rows; i-- > 0; ) {
2374 Row& s_i = s[i];
2375 Coefficient& s_i_pj = s_i[pj];
2376 product = s_i_pj * pivot_denom;
2377 if (product % pivot_pj != 0) {
2378 WEIGHT_ADD(98);
2379 // As above, perform row scaling.
2380 gcd_assign(gcd, product, pivot_pj);
2381 exact_div_assign(scale_factor, pivot_pj, gcd);
2382 for (Row::iterator
2383 k = s_i.begin(), k_end = s_i.end(); k != k_end; ++k) {
2384 WEIGHT_ADD(26);
2385 *k *= scale_factor;
2386 }
2387 product *= scale_factor;
2388 scaling[i] *= scale_factor;
2389 }
2390 PPL_ASSERT(product % pivot_pj == 0);
2391 exact_div_assign(s_i_pj, product, pivot_pj);
2392 WEIGHT_ADD(97);
2393 }
2394 }
2395 // Drop pivot to restore proper matrix size.
2396 s.remove_trailing_rows(1);
2397 }
2398
2399 // This point should be unreachable.
2400 PPL_UNREACHABLE;
2401 return false;
2402 }
2403
2404 void
2405 PIP_Solution_Node
update_tableau(const PIP_Problem & pip,const dimension_type external_space_dim,const dimension_type first_pending_constraint,const Constraint_Sequence & input_cs,const Variables_Set & parameters)2406 ::update_tableau(const PIP_Problem& pip,
2407 const dimension_type external_space_dim,
2408 const dimension_type first_pending_constraint,
2409 const Constraint_Sequence& input_cs,
2410 const Variables_Set& parameters) {
2411
2412 // Make sure a parameter column exists, for the inhomogeneous term.
2413 if (tableau.t.num_columns() == 0) {
2414 tableau.t.add_zero_columns(1);
2415 }
2416
2417 // NOTE: here 'params' stands for problem (i.e., non artificial) parameters.
2418 const dimension_type old_num_vars = tableau.s.num_columns();
2419 const dimension_type old_num_params
2420 = pip.internal_space_dim - old_num_vars;
2421 const dimension_type num_added_dims
2422 = pip.external_space_dim - pip.internal_space_dim;
2423 const dimension_type new_num_params = parameters.size();
2424 const dimension_type num_added_params = new_num_params - old_num_params;
2425 const dimension_type num_added_vars = num_added_dims - num_added_params;
2426
2427 // Resize the two tableau matrices.
2428 if (num_added_vars > 0) {
2429 tableau.s.add_zero_columns(num_added_vars);
2430 }
2431
2432 if (num_added_params > 0) {
2433 tableau.t.add_zero_columns(num_added_params, old_num_params + 1);
2434 }
2435
2436 dimension_type new_var_column = old_num_vars;
2437 const dimension_type initial_space_dim = old_num_vars + old_num_params;
2438 for (dimension_type i = initial_space_dim; i < external_space_dim; ++i) {
2439 if (parameters.count(i) == 0) {
2440 // A new problem variable.
2441 if (tableau.s.num_rows() == 0) {
2442 // No rows have been added yet
2443 basis.push_back(true);
2444 mapping.push_back(new_var_column);
2445 }
2446 else {
2447 /*
2448 Need to insert the original variable id
2449 before the slack variable id's to respect variable ordering.
2450 */
2451 basis.insert(nth_iter(basis, new_var_column), true);
2452 mapping.insert(nth_iter(mapping, new_var_column), new_var_column);
2453 // Update variable id's of slack variables.
2454 for (dimension_type j = var_row.size(); j-- > 0; ) {
2455 if (var_row[j] >= new_var_column) {
2456 ++var_row[j];
2457 }
2458 }
2459 for (dimension_type j = var_column.size(); j-- > 0; ) {
2460 if (var_column[j] >= new_var_column) {
2461 ++var_column[j];
2462 }
2463 }
2464 if (special_equality_row > 0) {
2465 ++special_equality_row;
2466 }
2467 }
2468 var_column.push_back(new_var_column);
2469 ++new_var_column;
2470 }
2471 }
2472
2473 if (big_dimension == not_a_dimension()
2474 && pip.big_parameter_dimension != not_a_dimension()) {
2475 // Compute the column number of big parameter in tableau.t matrix.
2476 Variables_Set::const_iterator pos
2477 = parameters.find(pip.big_parameter_dimension);
2478 big_dimension = 1U
2479 + static_cast<dimension_type>(std::distance(parameters.begin(), pos));
2480 }
2481
2482 Coefficient_traits::const_reference denom = tableau.denominator();
2483
2484 for (Constraint_Sequence::const_iterator
2485 c_iter = nth_iter(input_cs, first_pending_constraint),
2486 c_end = input_cs.end(); c_iter != c_end; ++c_iter) {
2487 const Constraint& constraint = *c_iter;
2488 // (Tentatively) Add new rows to s and t matrices.
2489 // These will be removed at the end if they turn out to be useless.
2490 const dimension_type row_id = tableau.s.num_rows();
2491 tableau.s.add_zero_rows(1);
2492 tableau.t.add_zero_rows(1);
2493 Row& v_row = tableau.s[row_id];
2494 Row& p_row = tableau.t[row_id];
2495
2496 {
2497 dimension_type p_index = 1;
2498 dimension_type v_index = 0;
2499 // Setting the inhomogeneous term.
2500 if (constraint.inhomogeneous_term() != 0) {
2501 Coefficient& p_row0 = p_row[0];
2502 p_row0 = constraint.inhomogeneous_term();
2503 if (constraint.is_strict_inequality()) {
2504 // Transform (expr > 0) into (expr - 1 >= 0).
2505 --p_row0;
2506 }
2507 p_row0 *= denom;
2508 }
2509 else {
2510 if (constraint.is_strict_inequality()) {
2511 // Transform (expr > 0) into (expr - 1 >= 0).
2512 neg_assign(p_row[0], denom);
2513 }
2514 }
2515 WEIGHT_BEGIN();
2516 dimension_type last_dim = 0;
2517 const Constraint::expr_type e = constraint.expression();
2518 for (Constraint::expr_type::const_iterator
2519 i = e.begin(), i_end = e.end(); i != i_end; ++i) {
2520 const dimension_type dim = i.variable().space_dimension();
2521 if (dim != last_dim + 1) {
2522 // We have skipped some zero coefficients.
2523 // Update p_index and v_index accordingly.
2524 const dimension_type n
2525 = std::distance(parameters.lower_bound(last_dim),
2526 parameters.lower_bound(dim - 1));
2527 const dimension_type num_skipped = dim - last_dim - 1;
2528 p_index += n;
2529 v_index += (num_skipped - n);
2530 }
2531 PPL_ASSERT(p_index + v_index == i.variable().id() + 1);
2532 const bool is_parameter = (1 == parameters.count(dim - 1));
2533 Coefficient_traits::const_reference coeff_i = *i;
2534
2535 WEIGHT_ADD(140);
2536 if (is_parameter) {
2537 p_row.insert(p_index, coeff_i * denom);
2538 ++p_index;
2539 }
2540 else {
2541 const dimension_type mv = mapping[v_index];
2542 if (basis[v_index]) {
2543 // Basic variable: add coeff_i * x_i
2544 add_mul_assign(v_row[mv], coeff_i, denom);
2545 }
2546 else {
2547 // Non-basic variable: add coeff_i * row_i
2548 add_mul_assign_row(v_row, coeff_i, tableau.s[mv]);
2549 add_mul_assign_row(p_row, coeff_i, tableau.t[mv]);
2550 }
2551 ++v_index;
2552 }
2553
2554 last_dim = dim;
2555 }
2556 }
2557
2558 if (row_sign(v_row, not_a_dimension()) == ZERO) {
2559 // Parametric-only constraints have already been inserted in
2560 // initial context, so no need to insert them in the tableau.
2561 tableau.s.remove_trailing_rows(1);
2562 tableau.t.remove_trailing_rows(1);
2563 }
2564 else {
2565 const dimension_type var_id = mapping.size();
2566 sign.push_back(row_sign(p_row, big_dimension));
2567 basis.push_back(false);
2568 mapping.push_back(row_id);
2569 var_row.push_back(var_id);
2570 if (constraint.is_equality()) {
2571 // Handle equality constraints.
2572 // After having added the f_i(x,p) >= 0 constraint,
2573 // we must add -f_i(x,p) to the special equality row.
2574 if (special_equality_row == 0 || basis[special_equality_row]) {
2575 // The special constraint has not been created yet
2576 // FIXME: for now, we do not handle the case where the variable
2577 // is basic, and we just create a new row.
2578 // This might be faster however.
2579 tableau.s.add_zero_rows(1);
2580 tableau.t.add_zero_rows(1);
2581 // NOTE: addition of rows invalidates references v_row and p_row
2582 // due to possible matrix reallocations: recompute them.
2583 neg_assign_row(tableau.s[1 + row_id], tableau.s[row_id]);
2584 neg_assign_row(tableau.t[1 + row_id], tableau.t[row_id]);
2585 sign.push_back(row_sign(tableau.t[1 + row_id], big_dimension));
2586 special_equality_row = mapping.size();
2587 basis.push_back(false);
2588 mapping.push_back(1 + row_id);
2589 var_row.push_back(1 + var_id);
2590 }
2591 else {
2592 // The special constraint already exists and is nonbasic.
2593 const dimension_type m_eq = mapping[special_equality_row];
2594 sub_assign(tableau.s[m_eq], v_row);
2595 sub_assign(tableau.t[m_eq], p_row);
2596 sign[m_eq] = row_sign(tableau.t[m_eq], big_dimension);
2597 }
2598 }
2599 }
2600 }
2601 PPL_ASSERT(OK());
2602 }
2603
2604 PIP_Tree_Node*
solve(const PIP_Problem & pip,const bool check_feasible_context,const Matrix<Row> & context,const Variables_Set & params,dimension_type space_dim,const int indent_level)2605 PIP_Solution_Node::solve(const PIP_Problem& pip,
2606 const bool check_feasible_context,
2607 const Matrix<Row>& context,
2608 const Variables_Set& params,
2609 dimension_type space_dim,
2610 const int indent_level) {
2611 PPL_ASSERT(indent_level >= 0);
2612 #ifdef NOISY_PIP_TREE_STRUCTURE
2613 indent_and_print(std::cerr, indent_level, "=== SOLVING NODE\n");
2614 #else
2615 PPL_USED(indent_level);
2616 #endif
2617 // Reset current solution as invalid.
2618 solution_valid = false;
2619
2620 Matrix<Row> ctx(context);
2621 Variables_Set all_params(params);
2622 const dimension_type num_art_params = artificial_parameters.size();
2623 add_artificial_parameters(ctx, all_params, space_dim, num_art_params);
2624 merge_assign(ctx, constraints_, all_params);
2625
2626 // If needed, (re-)check feasibility of context.
2627 if (check_feasible_context) {
2628 Matrix<Row> ctx_copy(ctx);
2629 if (!compatibility_check(ctx_copy)) {
2630 delete this;
2631 return 0;
2632 }
2633 }
2634
2635 const dimension_type not_a_dim = not_a_dimension();
2636
2637 // Main loop of the simplex algorithm.
2638 while (true) {
2639 // Check if the client has requested abandoning all expensive
2640 // computations. If so, the exception specified by the client
2641 // is thrown now.
2642 maybe_abandon();
2643
2644 PPL_ASSERT(OK());
2645
2646 const dimension_type num_rows = tableau.t.num_rows();
2647 const dimension_type num_vars = tableau.s.num_columns();
2648 const dimension_type num_params = tableau.t.num_columns();
2649 Coefficient_traits::const_reference tableau_denom = tableau.denominator();
2650
2651 #ifdef VERY_NOISY_PIP
2652 tableau.ascii_dump(std::cerr);
2653 std::cerr << "context ";
2654 ctx.ascii_dump(std::cerr);
2655 #endif // #ifdef VERY_NOISY_PIP
2656
2657 // (Re-) Compute parameter row signs.
2658 // While at it, keep track of the first parameter rows
2659 // having negative and mixed sign.
2660 dimension_type first_negative = not_a_dim;
2661 dimension_type first_mixed = not_a_dim;
2662 for (dimension_type i = 0; i < num_rows; ++i) {
2663 Row_Sign& sign_i = sign[i];
2664 if (sign_i == UNKNOWN || sign_i == MIXED) {
2665 sign_i = row_sign(tableau.t[i], big_dimension);
2666 }
2667
2668 if (sign_i == NEGATIVE && first_negative == not_a_dim) {
2669 first_negative = i;
2670 }
2671 else if (sign_i == MIXED && first_mixed == not_a_dim) {
2672 first_mixed = i;
2673 }
2674 }
2675
2676 // If no negative parameter row was found, try to refine the sign of
2677 // mixed rows using compatibility checks with the current context.
2678 if (first_negative == not_a_dim && first_mixed != not_a_dim) {
2679 for (dimension_type i = first_mixed; i < num_rows; ++i) {
2680 // Consider mixed sign parameter rows only.
2681 if (sign[i] != MIXED) {
2682 continue;
2683 }
2684 const Row& t_i = tableau.t[i];
2685 Row_Sign new_sign = ZERO;
2686 // Check compatibility for constraint t_i(z) >= 0.
2687 if (compatibility_check(ctx, t_i)) {
2688 new_sign = POSITIVE;
2689 }
2690 // Check compatibility for constraint t_i(z) < 0,
2691 // i.e., -t_i(z) - 1 >= 0.
2692 Row t_i_complement(num_params);
2693 complement_assign(t_i_complement, t_i, tableau_denom);
2694 if (compatibility_check(ctx, t_i_complement)) {
2695 new_sign = (new_sign == POSITIVE) ? MIXED : NEGATIVE;
2696 }
2697 // Update sign for parameter row i.
2698 sign[i] = new_sign;
2699 // Maybe update first_negative and first_mixed.
2700 if (new_sign == NEGATIVE && first_negative == not_a_dim) {
2701 first_negative = i;
2702 if (i == first_mixed) {
2703 first_mixed = not_a_dim;
2704 }
2705 }
2706 else if (new_sign == MIXED) {
2707 if (first_mixed == not_a_dim) {
2708 first_mixed = i;
2709 }
2710 }
2711 else if (i == first_mixed) {
2712 first_mixed = not_a_dim;
2713 }
2714 }
2715 }
2716
2717 // If there still is no negative parameter row and a mixed sign
2718 // parameter row (first_mixed) such that:
2719 // - it has at least one positive variable coefficient;
2720 // - constraint t_i(z) > 0 is not compatible with the context;
2721 // then this parameter row can be considered negative.
2722 if (first_negative == not_a_dim && first_mixed != not_a_dim) {
2723 WEIGHT_BEGIN();
2724 for (dimension_type i = first_mixed; i < num_rows; ++i) {
2725 // Consider mixed sign parameter rows only.
2726 if (sign[i] != MIXED) {
2727 continue;
2728 }
2729 // Check for a positive variable coefficient.
2730 const Row& s_i = tableau.s[i];
2731 bool has_positive = false;
2732 {
2733 for (Row::const_iterator
2734 j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
2735 if (*j > 0) {
2736 has_positive = true;
2737 break;
2738 }
2739 }
2740 }
2741 if (!has_positive) {
2742 continue;
2743 }
2744 // Check compatibility of constraint t_i(z) > 0.
2745 Row row(tableau.t[i]);
2746 PPL_DIRTY_TEMP_COEFFICIENT(mod);
2747 Coefficient& row0 = row[0];
2748 pos_rem_assign(mod, row0, tableau_denom);
2749 row0 -= (mod == 0) ? tableau_denom : mod;
2750 WEIGHT_ADD(210);
2751 const bool compatible = compatibility_check(ctx, row);
2752 // Maybe update sign (and first_* indices).
2753 if (compatible) {
2754 // Sign is still mixed.
2755 if (first_mixed == not_a_dim) {
2756 first_mixed = i;
2757 }
2758 }
2759 else {
2760 // Sign becomes negative (i.e., no longer mixed).
2761 sign[i] = NEGATIVE;
2762 if (first_negative == not_a_dim) {
2763 first_negative = i;
2764 }
2765 if (first_mixed == i) {
2766 first_mixed = not_a_dim;
2767 }
2768 }
2769 }
2770 }
2771
2772 #ifdef VERY_NOISY_PIP
2773 std::cerr << "sign =";
2774 for (dimension_type i = 0; i < sign.size(); ++i)
2775 std::cerr << " " << "?0+-*"[sign[i]];
2776 std::cerr << std::endl;
2777 #endif // #ifdef VERY_NOISY_PIP
2778
2779 // If we have found a negative parameter row, then
2780 // either the problem is unfeasible, or a pivoting step is required.
2781 if (first_negative != not_a_dim) {
2782
2783 // Search for the best pivot row.
2784 dimension_type pi = not_a_dim;
2785 dimension_type pj = not_a_dim;
2786 for (dimension_type i = first_negative; i < num_rows; ++i) {
2787 if (sign[i] != NEGATIVE) {
2788 continue;
2789 }
2790 dimension_type j;
2791 if (!find_lexico_minimal_column(tableau.s, mapping, basis,
2792 tableau.s[i], 0, j)) {
2793 // No positive s_ij was found: problem is unfeasible.
2794 #ifdef NOISY_PIP_TREE_STRUCTURE
2795 indent_and_print(std::cerr, indent_level,
2796 "No positive pivot: Solution = _|_\n");
2797 #endif // #ifdef NOISY_PIP_TREE_STRUCTURE
2798 delete this;
2799 return 0;
2800 }
2801 if (pj == not_a_dim
2802 || tableau.is_better_pivot(mapping, basis, i, j, pi, pj)) {
2803 // Update pivot indices.
2804 pi = i;
2805 pj = j;
2806 if (pip.control_parameters[PIP_Problem::PIVOT_ROW_STRATEGY]
2807 == PIP_Problem::PIVOT_ROW_STRATEGY_FIRST) {
2808 // Stop at first valid row.
2809 break;
2810 }
2811 }
2812 }
2813
2814 #ifdef VERY_NOISY_PIP
2815 std::cerr << "Pivot (pi, pj) = (" << pi << ", " << pj << ")\n";
2816 #endif // #ifdef VERY_NOISY_PIP
2817
2818 // Normalize the tableau before pivoting.
2819 tableau.normalize();
2820
2821 // Perform pivot operation.
2822
2823 // Update basis.
2824 {
2825 const dimension_type var_pi = var_row[pi];
2826 const dimension_type var_pj = var_column[pj];
2827 var_row[pi] = var_pj;
2828 var_column[pj] = var_pi;
2829 basis[var_pi] = true;
2830 basis[var_pj] = false;
2831 mapping[var_pi] = pj;
2832 mapping[var_pj] = pi;
2833 }
2834
2835 PPL_DIRTY_TEMP_COEFFICIENT(product);
2836 PPL_DIRTY_TEMP_COEFFICIENT(gcd);
2837 PPL_DIRTY_TEMP_COEFFICIENT(scale_factor);
2838
2839 // Creating identity rows corresponding to basic variable pj:
2840 // 1. add them to tableau so as to have proper size and capacity;
2841 tableau.s.add_zero_rows(1);
2842 tableau.t.add_zero_rows(1);
2843 // 2. swap the rows just added with empty ones.
2844 Row s_pivot(0);
2845 Row t_pivot(0);
2846 swap(s_pivot, tableau.s[num_rows]);
2847 swap(t_pivot, tableau.t[num_rows]);
2848 // 3. drop rows previously added at end of tableau.
2849 tableau.s.remove_trailing_rows(1);
2850 tableau.t.remove_trailing_rows(1);
2851
2852 // Save current pivot denominator.
2853 PPL_DIRTY_TEMP_COEFFICIENT(pivot_denom);
2854 pivot_denom = tableau.denominator();
2855 // Let the (scaled) pivot coordinate be 1.
2856 s_pivot[pj] = pivot_denom;
2857
2858 // Swap identity row with the pivot row previously found.
2859 s_pivot.m_swap(tableau.s[pi]);
2860 t_pivot.m_swap(tableau.t[pi]);
2861 sign[pi] = ZERO;
2862
2863 PPL_DIRTY_TEMP_COEFFICIENT(s_pivot_pj);
2864 s_pivot_pj = s_pivot.get(pj);
2865
2866 // Compute columns s[*][j]:
2867 //
2868 // <CODE>
2869 // s[i][j] -= s[i][pj] * s_pivot[j] / s_pivot_pj;
2870 // </CODE>
2871 for (dimension_type i = num_rows; i-- > 0; ) {
2872 Row& s_i = tableau.s[i];
2873 PPL_DIRTY_TEMP_COEFFICIENT(s_i_pj);
2874 s_i_pj = s_i.get(pj);
2875
2876 if (s_i_pj == 0) {
2877 continue;
2878 }
2879
2880 WEIGHT_BEGIN();
2881 Row::iterator itr = s_i.end();
2882 for (Row::const_iterator
2883 j = s_pivot.begin(), j_end = s_pivot.end(); j != j_end; ++j) {
2884 if (j.index() != pj) {
2885 Coefficient_traits::const_reference s_pivot_j = *j;
2886 // Do nothing if the j-th pivot element is zero.
2887 if (s_pivot_j != 0) {
2888 product = s_pivot_j * s_i_pj;
2889 if (product % s_pivot_pj != 0) {
2890 // Must scale matrix to stay in integer case.
2891 gcd_assign(gcd, product, s_pivot_pj);
2892 exact_div_assign(scale_factor, s_pivot_pj, gcd);
2893 tableau.scale(scale_factor);
2894 s_i_pj *= scale_factor;
2895 product *= scale_factor;
2896 WEIGHT_ADD(102);
2897 }
2898 PPL_ASSERT(product % s_pivot_pj == 0);
2899 exact_div_assign(product, product, s_pivot_pj);
2900 WEIGHT_ADD(130);
2901 if (product != 0) {
2902 itr = s_i.insert(itr, j.index());
2903 *itr -= product;
2904 WEIGHT_ADD(34);
2905 }
2906 }
2907 }
2908 }
2909 }
2910
2911 // Compute columns t[*][j]:
2912 //
2913 // <CODE>
2914 // t[i][j] -= s[i][pj] * t_pivot[j] / s_pivot_pj;
2915 // </CODE>
2916 for (dimension_type i = num_rows; i-- > 0; ) {
2917 Row& s_i = tableau.s[i];
2918 Row& t_i = tableau.t[i];
2919
2920 Row::iterator s_i_pj_itr = s_i.find(pj);
2921
2922 if (s_i_pj_itr == s_i.end()) {
2923 continue;
2924 }
2925
2926 // FIXME: the following comment does not make sense.
2927 // NOTE: This is a Coefficient& instead of a
2928 // Coefficient_traits::const_reference, because scale() may silently
2929 // modify it.
2930 const Coefficient& s_i_pj = *s_i_pj_itr;
2931
2932 if (s_i_pj == 0) {
2933 continue;
2934 }
2935
2936 WEIGHT_BEGIN();
2937 Row::iterator k = t_i.end();
2938 for (Row::const_iterator
2939 j = t_pivot.begin(), j_end = t_pivot.end(); j != j_end; ++j) {
2940 Coefficient_traits::const_reference t_pivot_j = *j;
2941 // Do nothing if the j-th pivot element is zero.
2942 if (t_pivot_j != 0) {
2943 product = t_pivot_j * s_i_pj;
2944 if (product % s_pivot_pj != 0) {
2945 // Must scale matrix to stay in integer case.
2946 gcd_assign(gcd, product, s_pivot_pj);
2947 exact_div_assign(scale_factor, s_pivot_pj, gcd);
2948 tableau.scale(scale_factor);
2949 product *= scale_factor;
2950 WEIGHT_ADD(261);
2951 }
2952 PPL_ASSERT(product % s_pivot_pj == 0);
2953 exact_div_assign(product, product, s_pivot_pj);
2954 WEIGHT_ADD(115);
2955 if (product != 0) {
2956 k = t_i.insert(k, j.index());
2957 *k -= product;
2958 WEIGHT_ADD(41);
2959 }
2960
2961 // Update row sign.
2962 Row_Sign& sign_i = sign[i];
2963 switch (sign_i) {
2964 case ZERO:
2965 if (product > 0) {
2966 sign_i = NEGATIVE;
2967 }
2968 else if (product < 0) {
2969 sign_i = POSITIVE;
2970 }
2971 break;
2972 case POSITIVE:
2973 if (product > 0) {
2974 sign_i = MIXED;
2975 }
2976 break;
2977 case NEGATIVE:
2978 if (product < 0) {
2979 sign_i = MIXED;
2980 }
2981 break;
2982 default:
2983 break;
2984 }
2985 }
2986 }
2987 }
2988
2989 // Compute column s[*][pj]: s[i][pj] /= s_pivot_pj;
2990 // Update column only if pivot coordinate != 1.
2991 if (s_pivot_pj != pivot_denom) {
2992 WEIGHT_BEGIN();
2993 Row::iterator itr;
2994 for (dimension_type i = num_rows; i-- > 0; ) {
2995 Row& s_i = tableau.s[i];
2996 itr = s_i.find(pj);
2997 if (itr == s_i.end()) {
2998 continue;
2999 }
3000 WEIGHT_ADD(43);
3001 product = *itr * pivot_denom;
3002 if (product % s_pivot_pj != 0) {
3003 // As above, perform matrix scaling.
3004 gcd_assign(gcd, product, s_pivot_pj);
3005 exact_div_assign(scale_factor, s_pivot_pj, gcd);
3006 tableau.scale(scale_factor);
3007 product *= scale_factor;
3008 WEIGHT_ADD(177);
3009 }
3010 PPL_ASSERT(product % s_pivot_pj == 0);
3011 if (product != 0 || *itr != 0) {
3012 WEIGHT_ADD(106);
3013 exact_div_assign(*itr, product, s_pivot_pj);
3014 }
3015 }
3016 }
3017
3018 // Pivoting process ended: jump to next iteration.
3019 continue;
3020 } // if (first_negative != not_a_dim)
3021
3022
3023 PPL_ASSERT(first_negative == not_a_dim);
3024 // If no negative parameter row was found,
3025 // but a mixed parameter row was found ...
3026 if (first_mixed != not_a_dim) {
3027 // Look for a constraint (i_neg):
3028 // - having mixed parameter sign;
3029 // - having no positive variable coefficient;
3030 // - minimizing the score (sum of parameter coefficients).
3031 dimension_type i_neg = not_a_dim;
3032 PPL_DIRTY_TEMP_COEFFICIENT(best_score);
3033 PPL_DIRTY_TEMP_COEFFICIENT(score);
3034 for (dimension_type i = first_mixed; i < num_rows; ++i) {
3035 // Mixed parameter sign.
3036 if (sign[i] != MIXED) {
3037 continue;
3038 }
3039 // No positive variable coefficient.
3040 bool has_positive = false;
3041 {
3042 const Row& s_i = tableau.s[i];
3043 for (Row::const_iterator
3044 j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
3045 if (*j > 0) {
3046 has_positive = true;
3047 break;
3048 }
3049 }
3050 }
3051 if (has_positive) {
3052 continue;
3053 }
3054 // Minimize parameter coefficient score,
3055 // eliminating implicated tautologies (if any).
3056 score = 0;
3057 {
3058 WEIGHT_BEGIN();
3059 const Row& t_i = tableau.t[i];
3060 for (Row::const_iterator
3061 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
3062 WEIGHT_ADD(55);
3063 score += *j;
3064 }
3065 }
3066 if (i_neg == not_a_dim || score < best_score) {
3067 i_neg = i;
3068 best_score = score;
3069 }
3070 }
3071
3072 if (i_neg != not_a_dim) {
3073 Row tautology = tableau.t[i_neg];
3074 /* Simplify tautology by exploiting integrality. */
3075 integral_simplification(tautology);
3076 ctx.add_row(tautology);
3077 add_constraint(tautology, all_params);
3078 sign[i_neg] = POSITIVE;
3079 #ifdef NOISY_PIP
3080 {
3081 Linear_Expression expr = Linear_Expression(tautology.get(0));
3082 dimension_type j = 1;
3083 for (Variables_Set::const_iterator p = all_params.begin(),
3084 p_end = all_params.end(); p != p_end; ++p, ++j)
3085 add_mul_assign(expr, tautology.get(j), Variable(*p));
3086 using namespace IO_Operators;
3087 std::cerr << std::setw(2 * indent_level) << ""
3088 << "Row " << i_neg
3089 << ": mixed param sign, negative var coeffs\n";
3090 std::cerr << std::setw(2 * indent_level) << ""
3091 << "==> adding tautology: "
3092 << Constraint(expr >= 0) << ".\n";
3093 }
3094 #endif // #ifdef NOISY_PIP
3095 // Jump to next iteration.
3096 continue;
3097 }
3098
3099 PPL_ASSERT(i_neg == not_a_dim);
3100 // Heuristically choose "best" (mixed) pivoting row.
3101 dimension_type best_i = not_a_dim;
3102 for (dimension_type i = first_mixed; i < num_rows; ++i) {
3103 if (sign[i] != MIXED) {
3104 continue;
3105 }
3106 score = 0;
3107 {
3108 WEIGHT_BEGIN();
3109 const Row& t_i = tableau.t[i];
3110 for (Row::const_iterator
3111 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
3112 WEIGHT_ADD(51);
3113 score += *j;
3114 }
3115 }
3116 if (best_i == not_a_dim || score < best_score) {
3117 best_score = score;
3118 best_i = i;
3119 }
3120 }
3121
3122 Row t_test(tableau.t[best_i]);
3123 /* Simplify t_test by exploiting integrality. */
3124 integral_simplification(t_test);
3125
3126 #ifdef NOISY_PIP
3127 {
3128 Linear_Expression expr = Linear_Expression(t_test.get(0));
3129 dimension_type j = 1;
3130 for (Variables_Set::const_iterator p = all_params.begin(),
3131 p_end = all_params.end(); p != p_end; ++p, ++j)
3132 add_mul_assign(expr, t_test.get(j), Variable(*p));
3133 using namespace IO_Operators;
3134 std::cerr << std::setw(2 * indent_level) << ""
3135 << "Row " << best_i << ": mixed param sign\n";
3136 std::cerr << std::setw(2 * indent_level) << ""
3137 << "==> depends on sign of " << expr << ".\n";
3138 }
3139 #endif // #ifdef NOISY_PIP
3140
3141 // Create a solution node for the "true" version of current node.
3142 PIP_Tree_Node* t_node = new PIP_Solution_Node(*this, No_Constraints());
3143 // Protect it from exception safety issues via std::auto_ptr.
3144 std::auto_ptr<PIP_Tree_Node> wrapped_node(t_node);
3145
3146 // Add parametric constraint to context.
3147 ctx.add_row(t_test);
3148 // Recursively solve true node with respect to updated context.
3149 #ifdef NOISY_PIP_TREE_STRUCTURE
3150 indent_and_print(std::cerr, indent_level, "=== SOLVING THEN CHILD\n");
3151 #endif
3152 t_node = t_node->solve(pip, check_feasible_context,
3153 ctx, all_params, space_dim,
3154 indent_level + 1);
3155 // Resolution may have changed t_node: in case, re-wrap it.
3156 if (t_node != wrapped_node.get()) {
3157 wrapped_node.release();
3158 wrapped_node.reset(t_node);
3159 }
3160
3161 // Modify *this in place to become the "false" version of current node.
3162 PIP_Tree_Node* f_node = this;
3163 // Swap aside constraints and artificial parameters
3164 // (these will be later restored if needed).
3165 Constraint_System cs;
3166 Artificial_Parameter_Sequence aps;
3167 swap(cs, f_node->constraints_);
3168 swap(aps, f_node->artificial_parameters);
3169 // Compute the complement of the constraint used for the "true" node.
3170 Row& f_test = ctx[ctx.num_rows() - 1];
3171 complement_assign(f_test, t_test, 1);
3172
3173 // Recursively solve false node with respect to updated context.
3174 #ifdef NOISY_PIP_TREE_STRUCTURE
3175 indent_and_print(std::cerr, indent_level, "=== SOLVING ELSE CHILD\n");
3176 #endif
3177 f_node = f_node->solve(pip, check_feasible_context,
3178 ctx, all_params, space_dim,
3179 indent_level + 1);
3180
3181 // Case analysis on recursive resolution calls outcome.
3182 if (t_node == 0) {
3183 if (f_node == 0) {
3184 // Both t_node and f_node unfeasible.
3185 #ifdef NOISY_PIP_TREE_STRUCTURE
3186 indent_and_print(std::cerr, indent_level,
3187 "=== EXIT: BOTH BRANCHES UNFEASIBLE: _|_\n");
3188 #endif
3189 return 0;
3190 }
3191 else {
3192 // t_node unfeasible, f_node feasible:
3193 // restore cs and aps into f_node (i.e., this).
3194 PPL_ASSERT(f_node == this);
3195 swap(f_node->constraints_, cs);
3196 swap(f_node->artificial_parameters, aps);
3197 // Add f_test to constraints.
3198 f_node->add_constraint(f_test, all_params);
3199 #ifdef NOISY_PIP_TREE_STRUCTURE
3200 indent_and_print(std::cerr, indent_level,
3201 "=== EXIT: THEN BRANCH UNFEASIBLE: SWAP BRANCHES\n");
3202 #endif
3203 return f_node;
3204 }
3205 }
3206 else if (f_node == 0) {
3207 // t_node feasible, f_node unfeasible.
3208 #ifdef NOISY_PIP_TREE_STRUCTURE
3209 indent_and_print(std::cerr, indent_level,
3210 "=== EXIT: THEN BRANCH FEASIBLE\n");
3211 #endif
3212 // NOTE: in principle, we could merge t_node into its parent.
3213 // However, if t_node is a decision node having both children,
3214 // then we would obtain a node violating the PIP_Decision_Node
3215 // invariant saying that t_node should have a single constraint:
3216 // it will have, at least, the two splitting constraints.
3217 const PIP_Decision_Node* const decision_node_p
3218 = dynamic_cast<PIP_Decision_Node*>(t_node);
3219 if (decision_node_p != 0 && decision_node_p->false_child != 0) {
3220 // Do NOT merge: create a new decision node.
3221 PIP_Tree_Node* const parent
3222 = new PIP_Decision_Node(t_node->get_owner(), 0, t_node);
3223 // Previously wrapped 't_node' is now safe: release it
3224 // and protect new 'parent' node from exception safety issues.
3225 wrapped_node.release();
3226 wrapped_node.reset(parent);
3227 // Restore into parent `cs' and `aps'.
3228 swap(parent->constraints_, cs);
3229 swap(parent->artificial_parameters, aps);
3230 // Add t_test to parent's constraints.
3231 parent->add_constraint(t_test, all_params);
3232 // It is now safe to release previously wrapped parent pointer
3233 // and return it to caller.
3234 return wrapped_node.release();
3235 }
3236 else {
3237 // Merge t_node with its parent:
3238 // a) append into `cs' the constraints of t_node;
3239 for (Constraint_System::const_iterator
3240 i = t_node->constraints_.begin(),
3241 i_end = t_node->constraints_.end(); i != i_end; ++i) {
3242 cs.insert(*i);
3243 }
3244 // b) append into `aps' the parameters of t_node;
3245 aps.insert(aps.end(),
3246 t_node->artificial_parameters.begin(),
3247 t_node->artificial_parameters.end());
3248 // c) swap the updated `cs' and `aps' into t_node.
3249 swap(cs, t_node->constraints_);
3250 swap(aps, t_node->artificial_parameters);
3251 // d) add t_test to t_nodes's constraints.
3252 t_node->add_constraint(t_test, all_params);
3253 // It is now safe to release previously wrapped t_node pointer
3254 // and return it to caller.
3255 return wrapped_node.release();
3256 }
3257 }
3258
3259 // Here both t_node and f_node are feasible:
3260 // create a new decision node.
3261 #ifdef NOISY_PIP_TREE_STRUCTURE
3262 indent_and_print(std::cerr, indent_level,
3263 "=== EXIT: BOTH BRANCHES FEASIBLE: NEW DECISION NODE\n");
3264 #endif
3265 PIP_Tree_Node* parent
3266 = new PIP_Decision_Node(f_node->get_owner(), f_node, t_node);
3267 // Previously wrapped 't_node' is now safe: release it
3268 // and protect new 'parent' node from exception safety issues.
3269 wrapped_node.release();
3270 wrapped_node.reset(parent);
3271
3272 // Add t_test to the constraints of the new decision node.
3273 parent->add_constraint(t_test, all_params);
3274
3275 if (!cs.empty()) {
3276 #ifdef NOISY_PIP_TREE_STRUCTURE
3277 indent_and_print(std::cerr, indent_level,
3278 "=== NODE HAS BOTH BRANCHES AND TAUTOLOGIES:\n");
3279 indent_and_print(std::cerr, indent_level,
3280 "=== CREATE NEW PARENT FOR TAUTOLOGIES\n");
3281 #endif
3282 // If node to be solved had tautologies,
3283 // store them in a new decision node.
3284 parent = new PIP_Decision_Node(parent->get_owner(), 0, parent);
3285 // Previously wrapped 'parent' node is now safe: release it
3286 // and protect new 'parent' node from exception safety issues.
3287 wrapped_node.release();
3288 wrapped_node.reset(parent);
3289 swap(parent->constraints_, cs);
3290 }
3291 swap(parent->artificial_parameters, aps);
3292 // It is now safe to release previously wrapped decision node
3293 // and return it to the caller.
3294 return wrapped_node.release();
3295 } // if (first_mixed != not_a_dim)
3296
3297
3298 PPL_ASSERT(first_negative == not_a_dim);
3299 PPL_ASSERT(first_mixed == not_a_dim);
3300 // Here all parameters are positive: we have found a continuous
3301 // solution. If the solution happens to be integer, then it is the
3302 // solution of the integer problem. Otherwise, we may need to generate
3303 // a new cut to try and get back into the integer case.
3304 #ifdef NOISY_PIP
3305 indent_and_print(std::cerr, indent_level,
3306 "All parameters are positive.\n");
3307 #endif // #ifdef NOISY_PIP
3308 tableau.normalize();
3309
3310 // Look for any row having non integer parameter coefficients.
3311 Coefficient_traits::const_reference denom = tableau.denominator();
3312 for (dimension_type k = 0; k < num_vars; ++k) {
3313 if (basis[k]) {
3314 // Basic variable = 0, hence integer.
3315 continue;
3316 }
3317 const dimension_type i = mapping[k];
3318 const Row& t_i = tableau.t[i];
3319 WEIGHT_BEGIN();
3320 for (Row::const_iterator
3321 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
3322 WEIGHT_ADD(27);
3323 if (*j % denom != 0) {
3324 goto non_integer;
3325 }
3326 }
3327 }
3328 // The goto was not taken, the solution is integer.
3329 #ifdef NOISY_PIP_TREE_STRUCTURE
3330 indent_and_print(std::cerr, indent_level,
3331 "EXIT: solution found.\n");
3332 #endif // #ifdef NOISY_PIP
3333 return this;
3334
3335 non_integer:
3336 // The solution is non-integer: generate a cut.
3337 PPL_DIRTY_TEMP_COEFFICIENT(mod);
3338 dimension_type best_i = not_a_dim;
3339 dimension_type best_pcount = not_a_dim;
3340
3341 const PIP_Problem::Control_Parameter_Value cutting_strategy
3342 = pip.control_parameters[PIP_Problem::CUTTING_STRATEGY];
3343
3344 if (cutting_strategy == PIP_Problem::CUTTING_STRATEGY_FIRST) {
3345 // Find the first row with simplest parametric part.
3346 for (dimension_type k = 0; k < num_vars; ++k) {
3347 if (basis[k]) {
3348 continue;
3349 }
3350 const dimension_type i = mapping[k];
3351 // Count the number of non-integer parameter coefficients.
3352 WEIGHT_BEGIN();
3353 dimension_type pcount = 0;
3354 const Row& t_i = tableau.t[i];
3355 for (Row::const_iterator
3356 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
3357 WEIGHT_ADD(18);
3358 pos_rem_assign(mod, *j, denom);
3359 if (mod != 0) {
3360 ++pcount;
3361 }
3362 }
3363 if (pcount > 0 && (best_i == not_a_dim || pcount < best_pcount)) {
3364 best_pcount = pcount;
3365 best_i = i;
3366 }
3367 }
3368 // Generate cut using 'best_i'.
3369 generate_cut(best_i, all_params, ctx, space_dim, indent_level);
3370 }
3371 else {
3372 PPL_ASSERT(cutting_strategy == PIP_Problem::CUTTING_STRATEGY_DEEPEST
3373 || cutting_strategy == PIP_Problem::CUTTING_STRATEGY_ALL);
3374 // Find the row with simplest parametric part
3375 // which will generate the "deepest" cut.
3376 PPL_DIRTY_TEMP_COEFFICIENT(best_score);
3377 best_score = 0;
3378 PPL_DIRTY_TEMP_COEFFICIENT(score);
3379 PPL_DIRTY_TEMP_COEFFICIENT(s_score);
3380 std::vector<dimension_type> all_best_is;
3381
3382 for (dimension_type k = 0; k < num_vars; ++k) {
3383 if (basis[k]) {
3384 continue;
3385 }
3386 const dimension_type i = mapping[k];
3387 // Compute score and pcount.
3388 score = 0;
3389 dimension_type pcount = 0;
3390 {
3391 WEIGHT_BEGIN();
3392 const Row& t_i = tableau.t[i];
3393 for (Row::const_iterator
3394 j = t_i.begin(), j_end = t_i.end(); j != j_end; ++j) {
3395 WEIGHT_ADD(46);
3396 pos_rem_assign(mod, *j, denom);
3397 if (mod != 0) {
3398 WEIGHT_ADD(94);
3399 score += denom;
3400 score -= mod;
3401 ++pcount;
3402 }
3403 }
3404 }
3405
3406 // Compute s_score.
3407 s_score = 0;
3408 {
3409 WEIGHT_BEGIN();
3410 const Row& s_i = tableau.s[i];
3411 for (Row::const_iterator
3412 j = s_i.begin(), j_end = s_i.end(); j != j_end; ++j) {
3413 WEIGHT_ADD(94);
3414 pos_rem_assign(mod, *j, denom);
3415 s_score += denom;
3416 s_score -= mod;
3417 }
3418 }
3419 // Combine 'score' and 's_score'.
3420 score *= s_score;
3421 /*
3422 Select row i if it is non integer AND
3423 - no row has been chosen yet; OR
3424 - it has fewer non-integer parameter coefficients; OR
3425 - it has the same number of non-integer parameter coefficients,
3426 but its score is greater.
3427 */
3428 if (pcount != 0
3429 && (best_i == not_a_dim
3430 || pcount < best_pcount
3431 || (pcount == best_pcount && score > best_score))) {
3432 if (pcount < best_pcount) {
3433 all_best_is.clear();
3434 }
3435 best_i = i;
3436 best_pcount = pcount;
3437 best_score = score;
3438 }
3439 if (pcount > 0) {
3440 all_best_is.push_back(i);
3441 }
3442 }
3443 if (cutting_strategy == PIP_Problem::CUTTING_STRATEGY_DEEPEST) {
3444 generate_cut(best_i, all_params, ctx, space_dim, indent_level);
3445 }
3446 else {
3447 PPL_ASSERT(cutting_strategy == PIP_Problem::CUTTING_STRATEGY_ALL);
3448 for (dimension_type k = all_best_is.size(); k-- > 0; ) {
3449 generate_cut(all_best_is[k], all_params, ctx,
3450 space_dim, indent_level);
3451 }
3452 }
3453 } // End of processing for non-integer solutions.
3454
3455 } // Main loop of the simplex algorithm
3456
3457 // This point should be unreachable.
3458 PPL_UNREACHABLE;
3459 return 0;
3460 }
3461
3462 void
generate_cut(const dimension_type index,Variables_Set & parameters,Matrix<Row> & context,dimension_type & space_dimension,const int indent_level)3463 PIP_Solution_Node::generate_cut(const dimension_type index,
3464 Variables_Set& parameters,
3465 Matrix<Row>& context,
3466 dimension_type& space_dimension,
3467 const int indent_level) {
3468 PPL_ASSERT(indent_level >= 0);
3469 #ifdef NOISY_PIP
3470 std::cerr << std::setw(2 * indent_level) << ""
3471 << "Row " << index << " requires cut generation.\n";
3472 #else
3473 PPL_USED(indent_level);
3474 #endif // #ifdef NOISY_PIP
3475
3476 const dimension_type num_rows = tableau.t.num_rows();
3477 PPL_ASSERT(index < num_rows);
3478 const dimension_type num_vars = tableau.s.num_columns();
3479 const dimension_type num_params = tableau.t.num_columns();
3480 PPL_ASSERT(num_params == 1 + parameters.size());
3481 Coefficient_traits::const_reference denom = tableau.denominator();
3482
3483 PPL_DIRTY_TEMP_COEFFICIENT(mod);
3484 PPL_DIRTY_TEMP_COEFFICIENT(coeff);
3485
3486 // Test if cut to be generated must be parametric or not.
3487 bool generate_parametric_cut = false;
3488 {
3489 // Limiting the scope of reference row_t (may be later invalidated).
3490 const Row& row_t = tableau.t[index];
3491 Row::const_iterator j = row_t.begin();
3492 Row::const_iterator j_end = row_t.end();
3493 // Skip the element with index 0.
3494 if (j != j_end && j.index() == 0) {
3495 ++j;
3496 }
3497 WEIGHT_BEGIN();
3498 for ( ; j != j_end; ++j) {
3499 WEIGHT_ADD(7);
3500 if (*j % denom != 0) {
3501 generate_parametric_cut = true;
3502 break;
3503 }
3504 }
3505 }
3506
3507 // Column index of already existing Artificial_Parameter.
3508 dimension_type ap_column = not_a_dimension();
3509
3510 if (generate_parametric_cut) {
3511 // Fractional parameter coefficient found: generate parametric cut.
3512 bool reuse_ap = false;
3513 Linear_Expression expr;
3514
3515 // Limiting the scope of reference row_t (may be later invalidated).
3516 {
3517 const Row& row_t = tableau.t[index];
3518 Row::const_iterator j = row_t.begin();
3519 Row::const_iterator j_end = row_t.end();
3520 if (j != j_end && j.index() == 0) {
3521 pos_rem_assign(mod, *j, denom);
3522 ++j;
3523 if (mod != 0) {
3524 // Optimizing computation: expr += (denom - mod);
3525 expr += denom;
3526 expr -= mod;
3527 }
3528 }
3529 if (!parameters.empty()) {
3530 // To avoid reallocations of expr.
3531 add_mul_assign(expr, 0, Variable(*(parameters.rbegin())));
3532 Variables_Set::const_iterator p_j = parameters.begin();
3533 dimension_type last_index = 1;
3534 WEIGHT_BEGIN();
3535 for ( ; j != j_end; ++j) {
3536 WEIGHT_ADD(69);
3537 pos_rem_assign(mod, *j, denom);
3538 if (mod != 0) {
3539 // Optimizing computation: expr += (denom - mod) * Variable(*p_j);
3540 WEIGHT_ADD(164);
3541 coeff = denom - mod;
3542 PPL_ASSERT(last_index <= j.index());
3543 std::advance(p_j, j.index() - last_index);
3544 last_index = j.index();
3545 add_mul_assign(expr, coeff, Variable(*p_j));
3546 }
3547 }
3548 }
3549 }
3550 // Generate new artificial parameter.
3551 const Artificial_Parameter ap(expr, denom);
3552
3553 // Search if the Artificial_Parameter has already been generated.
3554 ap_column = space_dimension;
3555 const PIP_Tree_Node* node = this;
3556 do {
3557 for (dimension_type j = node->artificial_parameters.size(); j-- > 0; ) {
3558 --ap_column;
3559 if (node->artificial_parameters[j] == ap) {
3560 reuse_ap = true;
3561 break;
3562 }
3563 }
3564 node = node->parent();
3565 } while (!reuse_ap && node != 0);
3566
3567 if (reuse_ap) {
3568 // We can re-use an existing Artificial_Parameter.
3569 #ifdef NOISY_PIP
3570 using namespace IO_Operators;
3571 std::cerr << std::setw(2 * indent_level) << ""
3572 << "Re-using parameter " << Variable(ap_column)
3573 << " = " << ap << std::endl;
3574 #endif // #ifdef NOISY_PIP
3575 ap_column = ap_column - num_vars + 1;
3576 }
3577 else {
3578 // Here reuse_ap == false: the Artificial_Parameter does not exist yet.
3579 // Beware: possible reallocation invalidates row references.
3580 tableau.t.add_zero_columns(1);
3581 context.add_zero_columns(1);
3582 artificial_parameters.push_back(ap);
3583 parameters.insert(space_dimension);
3584 #ifdef NOISY_PIP
3585 using namespace IO_Operators;
3586 std::cerr << std::setw(2 * indent_level) << ""
3587 << "New parameter " << Variable(space_dimension)
3588 << " = " << ap << std::endl;
3589 #endif // #ifdef NOISY_PIP
3590 ++space_dimension;
3591 ap_column = num_params;
3592
3593 // Update current context with constraints on the new parameter.
3594 const dimension_type ctx_num_rows = context.num_rows();
3595 context.add_zero_rows(2);
3596 Row& ctx1 = context[ctx_num_rows];
3597 Row& ctx2 = context[ctx_num_rows+1];
3598 // Recompute row reference after possible reallocation.
3599 const Row& row_t = tableau.t[index];
3600 {
3601 Row::const_iterator j = row_t.begin();
3602 Row::const_iterator j_end = row_t.end();
3603 Row::iterator itr1 = ctx1.end();
3604 Row::iterator itr2 = ctx2.end();
3605 if (j != j_end && j.index() == 0) {
3606 pos_rem_assign(mod, *j, denom);
3607 if (mod != 0) {
3608 itr1 = ctx1.insert(0, denom);
3609 *itr1 -= mod;
3610 itr2 = ctx2.insert(0, *itr1);
3611 neg_assign(*itr2);
3612 // Compute <CODE> ctx2[0] += denom-1; </CODE>
3613 *itr2 += denom;
3614 --(*itr2);
3615 }
3616 else {
3617 // Compute <CODE> ctx2[0] += denom-1; </CODE>
3618 itr2 = ctx2.insert(0, denom);
3619 --(*itr2);
3620 }
3621 ++j;
3622 }
3623 else {
3624 // Compute <CODE> ctx2[0] += denom-1; </CODE>
3625 itr2 = ctx2.insert(0, denom);
3626 --(*itr2);
3627 }
3628 WEIGHT_BEGIN();
3629 for ( ; j != j_end; ++j) {
3630 pos_rem_assign(mod, *j, denom);
3631 if (mod != 0) {
3632 const dimension_type j_index = j.index();
3633 itr1 = ctx1.insert(itr1, j_index, denom);
3634 *itr1 -= mod;
3635 itr2 = ctx2.insert(itr2, j_index, *itr1);
3636 neg_assign(*itr2);
3637 WEIGHT_ADD(294);
3638 }
3639 }
3640 itr1 = ctx1.insert(itr1, num_params, denom);
3641 neg_assign(*itr1);
3642 itr2 = ctx2.insert(itr2, num_params, denom);
3643 WEIGHT_ADD(122);
3644 }
3645
3646 #ifdef NOISY_PIP
3647 {
3648 using namespace IO_Operators;
3649 Variables_Set::const_iterator p = parameters.begin();
3650 Linear_Expression expr1(ctx1.get(0));
3651 Linear_Expression expr2(ctx2.get(0));
3652 for (dimension_type j = 1; j <= num_params; ++j, ++p) {
3653 add_mul_assign(expr1, ctx1.get(j), Variable(*p));
3654 add_mul_assign(expr2, ctx2.get(j), Variable(*p));
3655 }
3656 std::cerr << std::setw(2 * indent_level) << ""
3657 << "Adding to context: "
3658 << Constraint(expr1 >= 0) << " ; "
3659 << Constraint(expr2 >= 0) << std::endl;
3660 }
3661 #endif // #ifdef NOISY_PIP
3662 }
3663 }
3664
3665 // Generate new cut.
3666 tableau.s.add_zero_rows(1);
3667 tableau.t.add_zero_rows(1);
3668 Row& cut_s = tableau.s[num_rows];
3669 Row& cut_t = tableau.t[num_rows];
3670 // Recompute references after possible reallocation.
3671 const Row& row_s = tableau.s[index];
3672 const Row& row_t = tableau.t[index];
3673 WEIGHT_BEGIN();
3674 {
3675 Row::iterator itr = cut_s.end();
3676 for (Row::const_iterator
3677 j = row_s.begin(), j_end = row_s.end(); j != j_end; ++j) {
3678 WEIGHT_ADD(55);
3679 itr = cut_s.insert(itr, j.index(), *j);
3680 pos_rem_assign(*itr, *itr, denom);
3681 }
3682 }
3683 {
3684 Row::iterator cut_t_itr = cut_t.end();
3685 for (Row::const_iterator
3686 j = row_t.begin(), j_end = row_t.end(); j!=j_end; ++j) {
3687 WEIGHT_ADD(37);
3688 pos_rem_assign(mod, *j, denom);
3689 if (mod != 0) {
3690 WEIGHT_ADD(108);
3691 cut_t_itr = cut_t.insert(cut_t_itr, j.index(), mod);
3692 *cut_t_itr -= denom;
3693 }
3694 }
3695 }
3696 if (ap_column != not_a_dimension()) {
3697 // If we re-use an existing Artificial_Parameter
3698 cut_t[ap_column] = denom;
3699 }
3700 #ifdef NOISY_PIP
3701 {
3702 using namespace IO_Operators;
3703 Linear_Expression expr;
3704 dimension_type ti = 1;
3705 dimension_type si = 0;
3706 for (dimension_type j = 0; j < space_dimension; ++j) {
3707 if (parameters.count(j) == 1) {
3708 add_mul_assign(expr, cut_t.get(ti), Variable(j));
3709 ++ti;
3710 }
3711 else {
3712 add_mul_assign(expr, cut_s.get(si), Variable(j));
3713 ++si;
3714 }
3715 }
3716 std::cerr << std::setw(2 * indent_level) << ""
3717 << "Adding cut: "
3718 << Constraint(expr + cut_t.get(0) >= 0)
3719 << std::endl;
3720 }
3721 #endif // #ifdef NOISY_PIP
3722 var_row.push_back(num_rows + num_vars);
3723 basis.push_back(false);
3724 mapping.push_back(num_rows);
3725 sign.push_back(NEGATIVE);
3726 }
3727
3728
3729 memory_size_type
external_memory_in_bytes() const3730 PIP_Tree_Node::Artificial_Parameter::external_memory_in_bytes() const {
3731 return Linear_Expression::external_memory_in_bytes()
3732 + Parma_Polyhedra_Library::external_memory_in_bytes(denom);
3733 }
3734
3735 memory_size_type
total_memory_in_bytes() const3736 PIP_Tree_Node::Artificial_Parameter::total_memory_in_bytes() const {
3737 return sizeof(*this) + external_memory_in_bytes();
3738 }
3739
3740 memory_size_type
external_memory_in_bytes() const3741 PIP_Tree_Node::external_memory_in_bytes() const {
3742 memory_size_type n = constraints_.external_memory_in_bytes();
3743 // Adding the external memory for `artificial_parameters'.
3744 n += artificial_parameters.capacity() * sizeof(Artificial_Parameter);
3745 for (Artificial_Parameter_Sequence::const_iterator
3746 ap = art_parameter_begin(),
3747 ap_end = art_parameter_end(); ap != ap_end; ++ap) {
3748 n += (ap->external_memory_in_bytes());
3749 }
3750
3751 return n;
3752 }
3753
3754 memory_size_type
external_memory_in_bytes() const3755 PIP_Decision_Node::external_memory_in_bytes() const {
3756 memory_size_type n = PIP_Tree_Node::external_memory_in_bytes();
3757 PPL_ASSERT(true_child != 0);
3758 n += true_child->total_memory_in_bytes();
3759 if (false_child != 0) {
3760 n += false_child->total_memory_in_bytes();
3761 }
3762 return n;
3763 }
3764
3765 memory_size_type
total_memory_in_bytes() const3766 PIP_Decision_Node::total_memory_in_bytes() const {
3767 return sizeof(*this) + external_memory_in_bytes();
3768 }
3769
3770 memory_size_type
external_memory_in_bytes() const3771 PIP_Solution_Node::Tableau::external_memory_in_bytes() const {
3772 return Parma_Polyhedra_Library::external_memory_in_bytes(denom)
3773 + s.external_memory_in_bytes()
3774 + t.external_memory_in_bytes();
3775 }
3776
3777 memory_size_type
external_memory_in_bytes() const3778 PIP_Solution_Node::external_memory_in_bytes() const {
3779 memory_size_type n = PIP_Tree_Node::external_memory_in_bytes();
3780 n += tableau.external_memory_in_bytes();
3781 // FIXME: size of std::vector<bool> ?
3782 n += basis.capacity() * sizeof(bool);
3783 n += sizeof(dimension_type)
3784 * (mapping.capacity() + var_row.capacity() + var_column.capacity());
3785 n += sign.capacity() * sizeof(Row_Sign);
3786 // FIXME: Adding the external memory for `solution'.
3787 n += solution.capacity() * sizeof(Linear_Expression);
3788 for (std::vector<Linear_Expression>::const_iterator
3789 i = solution.begin(), i_end = solution.end();
3790 i != i_end; ++i) {
3791 n += (i->external_memory_in_bytes());
3792 }
3793
3794 return n;
3795 }
3796
3797 memory_size_type
total_memory_in_bytes() const3798 PIP_Solution_Node::total_memory_in_bytes() const {
3799 return sizeof(*this) + external_memory_in_bytes();
3800 }
3801
3802 void
indent_and_print(std::ostream & s,const int indent,const char * str)3803 PIP_Tree_Node::indent_and_print(std::ostream& s,
3804 const int indent,
3805 const char* str) {
3806 PPL_ASSERT(indent >= 0);
3807 s << std::setw(2 * indent) << "" << str;
3808 }
3809
3810 void
print(std::ostream & s,const int indent) const3811 PIP_Tree_Node::print(std::ostream& s, const int indent) const {
3812 const dimension_type pip_space_dim = get_owner()->space_dimension();
3813 const Variables_Set& pip_params = get_owner()->parameter_space_dimensions();
3814
3815 std::vector<bool> pip_dim_is_param(pip_space_dim);
3816 for (Variables_Set::const_iterator p = pip_params.begin(),
3817 p_end = pip_params.end(); p != p_end; ++p) {
3818 pip_dim_is_param[*p] = true;
3819 }
3820
3821 dimension_type first_art_dim = pip_space_dim;
3822 for (const PIP_Tree_Node* node = parent();
3823 node != 0; node = node->parent()) {
3824 first_art_dim += node->art_parameter_count();
3825 }
3826
3827 print_tree(s, indent, pip_dim_is_param, first_art_dim);
3828 }
3829
3830 void
print_tree(std::ostream & s,const int indent,const std::vector<bool> &,dimension_type first_art_dim) const3831 PIP_Tree_Node::print_tree(std::ostream& s, const int indent,
3832 const std::vector<bool>&,
3833 dimension_type first_art_dim) const {
3834 using namespace IO_Operators;
3835
3836 // Print artificial parameters.
3837 for (Artificial_Parameter_Sequence::const_iterator
3838 api = art_parameter_begin(),
3839 api_end = art_parameter_end(); api != api_end; ++api) {
3840 indent_and_print(s, indent, "Parameter ");
3841 s << Variable(first_art_dim) << " = " << *api << "\n";
3842 ++first_art_dim;
3843 }
3844
3845 // Print constraints, if any.
3846 if (!constraints_.empty()) {
3847 indent_and_print(s, indent, "if ");
3848
3849 Constraint_System::const_iterator ci = constraints_.begin();
3850 Constraint_System::const_iterator ci_end = constraints_.end();
3851 PPL_ASSERT(ci != ci_end);
3852 s << *ci;
3853 for (++ci; ci != ci_end; ++ci) {
3854 s << " and " << *ci;
3855 }
3856
3857 s << " then\n";
3858 }
3859 }
3860
3861 void
print_tree(std::ostream & s,const int indent,const std::vector<bool> & pip_dim_is_param,const dimension_type first_art_dim) const3862 PIP_Decision_Node::print_tree(std::ostream& s, const int indent,
3863 const std::vector<bool>& pip_dim_is_param,
3864 const dimension_type first_art_dim) const {
3865 // First print info common to decision and solution nodes.
3866 PIP_Tree_Node::print_tree(s, indent, pip_dim_is_param, first_art_dim);
3867
3868 // Then print info specific of decision nodes.
3869 const dimension_type child_first_art_dim
3870 = first_art_dim + art_parameter_count();
3871
3872 PPL_ASSERT(true_child != 0);
3873 true_child->print_tree(s, indent+1, pip_dim_is_param, child_first_art_dim);
3874
3875 indent_and_print(s, indent, "else\n");
3876
3877 if (false_child != 0) {
3878 false_child->print_tree(s, indent+1, pip_dim_is_param,
3879 child_first_art_dim);
3880 }
3881 else {
3882 indent_and_print(s, indent+1, "_|_\n");
3883 }
3884 }
3885
3886 void
print_tree(std::ostream & s,const int indent,const std::vector<bool> & pip_dim_is_param,const dimension_type first_art_dim) const3887 PIP_Solution_Node::print_tree(std::ostream& s, const int indent,
3888 const std::vector<bool>& pip_dim_is_param,
3889 const dimension_type first_art_dim) const {
3890 // Print info common to decision and solution nodes.
3891 PIP_Tree_Node::print_tree(s, indent, pip_dim_is_param, first_art_dim);
3892
3893 // Print info specific of solution nodes:
3894 // first update solution if needed ...
3895 update_solution(pip_dim_is_param);
3896 // ... and then actually print it.
3897 const bool no_constraints = constraints_.empty();
3898 indent_and_print(s, indent + (no_constraints ? 0 : 1), "{");
3899 const dimension_type pip_space_dim = pip_dim_is_param.size();
3900 for (dimension_type i = 0, num_var = 0; i < pip_space_dim; ++i) {
3901 if (pip_dim_is_param[i]) {
3902 continue;
3903 }
3904 if (num_var > 0) {
3905 s << " ; ";
3906 }
3907 using namespace IO_Operators;
3908 s << solution[num_var];
3909 ++num_var;
3910 }
3911 s << "}\n";
3912
3913 if (!no_constraints) {
3914 indent_and_print(s, indent, "else\n");
3915 indent_and_print(s, indent+1, "_|_\n");
3916 }
3917 }
3918
3919 const Linear_Expression&
parametric_values(const Variable var) const3920 PIP_Solution_Node::parametric_values(const Variable var) const {
3921 const PIP_Problem* const pip = get_owner();
3922 PPL_ASSERT(pip != 0);
3923
3924 const dimension_type space_dim = pip->space_dimension();
3925 if (var.space_dimension() > space_dim) {
3926 std::ostringstream s;
3927 s << "PPL::PIP_Solution_Node::parametric_values(v):\n"
3928 << "v.space_dimension() == " << var.space_dimension()
3929 << " is incompatible with the owning PIP_Problem "
3930 << " (space dim == " << space_dim << ").";
3931 throw std::invalid_argument(s.str());
3932 }
3933
3934 dimension_type solution_index = var.id();
3935 const Variables_Set& params = pip->parameter_space_dimensions();
3936 for (Variables_Set::const_iterator p = params.begin(),
3937 p_end = params.end(); p != p_end; ++p) {
3938 const dimension_type param_index = *p;
3939 if (param_index < var.id()) {
3940 --solution_index;
3941 }
3942 else if (param_index == var.id()) {
3943 throw std::invalid_argument("PPL::PIP_Solution_Node"
3944 "::parametric_values(v):\n"
3945 "v is a problem parameter.");
3946 }
3947 else {
3948 break;
3949 }
3950 }
3951
3952 update_solution();
3953 return solution[solution_index];
3954 }
3955
3956
3957 void
update_solution() const3958 PIP_Solution_Node::update_solution() const {
3959 // Avoid doing useless work.
3960 if (solution_valid) {
3961 return;
3962 }
3963 const PIP_Problem* const pip = get_owner();
3964 PPL_ASSERT(pip != 0);
3965 std::vector<bool> pip_dim_is_param(pip->space_dimension());
3966 const Variables_Set& params = pip->parameter_space_dimensions();
3967 for (Variables_Set::const_iterator p = params.begin(),
3968 p_end = params.end(); p != p_end; ++p) {
3969 pip_dim_is_param[*p] = true;
3970 }
3971
3972 update_solution(pip_dim_is_param);
3973 }
3974
3975 void
3976 PIP_Solution_Node
update_solution(const std::vector<bool> & pip_dim_is_param) const3977 ::update_solution(const std::vector<bool>& pip_dim_is_param) const {
3978 // Avoid doing useless work.
3979 if (solution_valid) {
3980 return;
3981 }
3982
3983 // const_cast required so as to refresh the solution cache.
3984 PIP_Solution_Node& x = const_cast<PIP_Solution_Node&>(*this);
3985
3986 const dimension_type num_pip_dims = pip_dim_is_param.size();
3987 const dimension_type num_pip_vars = tableau.s.num_columns();
3988 const dimension_type num_pip_params = num_pip_dims - num_pip_vars;
3989 const dimension_type num_all_params = tableau.t.num_columns() - 1;
3990 const dimension_type num_art_params = num_all_params - num_pip_params;
3991
3992 if (solution.size() != num_pip_vars) {
3993 x.solution.resize(num_pip_vars);
3994 }
3995
3996 // Compute external "names" (i.e., indices) for all parameters.
3997 std::vector<dimension_type> all_param_names(num_all_params);
3998
3999 // External indices for problem parameters.
4000 for (dimension_type i = 0, p_index = 0; i < num_pip_dims; ++i) {
4001 if (pip_dim_is_param[i]) {
4002 all_param_names[p_index] = i;
4003 ++p_index;
4004 }
4005 }
4006 // External indices for artificial parameters.
4007 for (dimension_type i = 0; i < num_art_params; ++i) {
4008 all_param_names[num_pip_params + i] = num_pip_dims + i;
4009 }
4010
4011
4012 PPL_DIRTY_TEMP_COEFFICIENT(norm_coeff);
4013 Coefficient_traits::const_reference denom = tableau.denominator();
4014 for (dimension_type i = num_pip_vars; i-- > 0; ) {
4015 Linear_Expression& sol_i = x.solution[i];
4016 sol_i = Linear_Expression(0);
4017 if (basis[i]) {
4018 continue;
4019 }
4020 const Row& row = tableau.t[mapping[i]];
4021
4022 // Start from index 1 to skip the inhomogeneous term.
4023 Row::const_iterator j = row.begin();
4024 Row::const_iterator j_end = row.end();
4025 // Skip the element with index 0.
4026 if (j != j_end && j.index() == 0) {
4027 ++j;
4028 }
4029 for ( ; j != j_end; ++j) {
4030 Coefficient_traits::const_reference coeff = *j;
4031 if (coeff == 0) {
4032 continue;
4033 }
4034 norm_coeff = coeff / denom;
4035 if (norm_coeff != 0) {
4036 add_mul_assign(sol_i, norm_coeff,
4037 Variable(all_param_names[j.index() - 1]));
4038 }
4039 }
4040 norm_coeff = row.get(0) / denom;
4041 sol_i += norm_coeff;
4042 }
4043
4044 // Mark solution as valid.
4045 x.solution_valid = true;
4046 }
4047
4048 } // namespace Parma_Polyhedra_Library
4049