00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include <config.h>
00024 #include "LP_Problem.defs.hh"
00025 #include "globals.types.hh"
00026 #include "globals.defs.hh"
00027 #include "Row.defs.hh"
00028 #include "Matrix.defs.hh"
00029 #include "Linear_Row.defs.hh"
00030 #include "Linear_System.defs.hh"
00031 #include "Linear_Expression.defs.hh"
00032 #include "Constraint_System.defs.hh"
00033 #include "Constraint_System.inlines.hh"
00034 #include "Generator.defs.hh"
00035 #include <stdexcept>
00036 #include <sstream>
00037 #include <map>
00038 #include <deque>
00039 #include <set>
00040 #include <algorithm>
00041
00042 #ifndef PPL_NOISY_SIMPLEX
00043 #define PPL_NOISY_SIMPLEX 0
00044 #endif
00045
00046 #ifdef PPL_NOISY_SIMPLEX
00047 #include <iostream>
00048 #endif
00049
00050 #ifndef PPL_SIMPLEX_ENABLE_STEEPEST_EDGE
00051 #define PPL_SIMPLEX_ENABLE_STEEPEST_EDGE 1
00052 #endif
00053
00054 namespace PPL = Parma_Polyhedra_Library;
00055
00056 #if PPL_NOISY_SIMPLEX
00057 namespace {
00058
00059 unsigned long num_iterations = 0;
00060
00061 }
00062 #endif // PPL_NOISY_SIMPLEX
00063
00064 PPL::dimension_type
00065 PPL::LP_Problem::steepest_edge() const {
00066 const dimension_type tableau_num_rows = tableau.num_rows();
00067 assert(tableau_num_rows == base.size());
00068
00069 TEMP_INTEGER(squared_lcm_basis);
00070
00071 std::vector<Coefficient> norm_factor(tableau_num_rows);
00072 {
00073
00074 TEMP_INTEGER(lcm_basis);
00075 lcm_basis = 1;
00076 for (dimension_type i = tableau_num_rows; i-- > 0; )
00077 lcm_assign(lcm_basis, lcm_basis, tableau[i][base[i]]);
00078
00079 for (dimension_type i = tableau_num_rows; i-- > 0; )
00080 exact_div_assign(norm_factor[i], lcm_basis, tableau[i][base[i]]);
00081
00082
00083 lcm_basis *= lcm_basis;
00084 std::swap(squared_lcm_basis, lcm_basis);
00085 }
00086
00087
00088 TEMP_INTEGER(challenger_num);
00089 TEMP_INTEGER(scalar_value);
00090 TEMP_INTEGER(challenger_den);
00091 TEMP_INTEGER(challenger_value);
00092 TEMP_INTEGER(current_value);
00093
00094 TEMP_INTEGER(current_num);
00095 TEMP_INTEGER(current_den);
00096 dimension_type entering_index = 0;
00097 const int cost_sign = sgn(working_cost[working_cost.size() - 1]);
00098 for (dimension_type j = tableau.num_columns() - 1; j-- > 1; ) {
00099 const Coefficient& cost_j = working_cost[j];
00100 if (sgn(cost_j) == cost_sign) {
00101
00102
00103 challenger_num = cost_j * cost_j;
00104
00105
00106 challenger_den = squared_lcm_basis;
00107 for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00108 const Coefficient& tableau_ij = tableau[i][j];
00109
00110 if (tableau_ij != 0) {
00111 scalar_value = tableau_ij * norm_factor[i];
00112 add_mul_assign(challenger_den, scalar_value, scalar_value);
00113 }
00114 }
00115
00116 if (entering_index == 0) {
00117 std::swap(current_num, challenger_num);
00118 std::swap(current_den, challenger_den);
00119 entering_index = j;
00120 continue;
00121 }
00122 challenger_value = challenger_num * current_den;
00123 current_value = current_num * challenger_den;
00124
00125 if (challenger_value > current_value) {
00126 std::swap(current_num, challenger_num);
00127 std::swap(current_den, challenger_den);
00128 entering_index = j;
00129 }
00130 }
00131 }
00132 return entering_index;
00133 }
00134
00135
00136
00137 PPL::dimension_type
00138 PPL::LP_Problem::get_entering_var_index() const {
00139
00140
00141
00142
00143
00144
00145 const dimension_type cost_sign_index = working_cost.size() - 1;
00146 const int cost_sign = sgn(working_cost[cost_sign_index]);
00147 assert(cost_sign != 0);
00148 for (dimension_type i = 1; i < cost_sign_index; ++i)
00149 if (sgn(working_cost[i]) == cost_sign)
00150 return i;
00151
00152
00153 return 0;
00154 }
00155
00156 void
00157 PPL::LP_Problem::linear_combine(Row& x,
00158 const Row& y,
00159 const dimension_type k) {
00160 assert(x.size() == y.size());
00161 assert(y[k] != 0 && x[k] != 0);
00162
00163
00164
00165 TEMP_INTEGER(normalized_x_k);
00166 TEMP_INTEGER(normalized_y_k);
00167 normalize2(x[k], y[k], normalized_x_k, normalized_y_k);
00168 for (dimension_type i = x.size(); i-- > 0; )
00169 if (i != k) {
00170 Coefficient& x_i = x[i];
00171 x_i *= normalized_y_k;
00172
00173 const Coefficient& y_i = y[i];
00174 if (y_i != 0)
00175 sub_mul_assign(x_i, y_i, normalized_x_k);
00176 }
00177 x[k] = 0;
00178 x.normalize();
00179 }
00180
00181
00182
00183 void
00184 PPL::LP_Problem::swap_base(const dimension_type entering_var_index,
00185 const dimension_type exiting_base_index) {
00186 const Row& tableau_out = tableau[exiting_base_index];
00187
00188 for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
00189 Row& tableau_i = tableau[i];
00190 if (i != exiting_base_index && tableau_i[entering_var_index] != 0)
00191 linear_combine(tableau_i, tableau_out, entering_var_index);
00192 }
00193
00194 if (working_cost[entering_var_index] != 0)
00195 linear_combine(working_cost, tableau_out, entering_var_index);
00196
00197 base[exiting_base_index] = entering_var_index;
00198 }
00199
00200
00201
00202 PPL::dimension_type
00203 PPL::LP_Problem
00204 ::get_exiting_base_index(const dimension_type entering_var_index) const {
00205
00206
00207
00208
00209
00210
00211
00212 const dimension_type tableau_num_rows = tableau.num_rows();
00213 dimension_type exiting_base_index = tableau_num_rows;
00214 for (dimension_type i = 0; i < tableau_num_rows; ++i) {
00215 const Row& t_i = tableau[i];
00216 const int num_sign = sgn(t_i[entering_var_index]);
00217 if (num_sign != 0 && num_sign == sgn(t_i[base[i]])) {
00218 exiting_base_index = i;
00219 break;
00220 }
00221 }
00222
00223 if (exiting_base_index == tableau_num_rows)
00224 return tableau_num_rows;
00225
00226
00227 TEMP_INTEGER(lcm);
00228 TEMP_INTEGER(current_min);
00229 TEMP_INTEGER(challenger);
00230 for (dimension_type i = exiting_base_index + 1; i < tableau_num_rows; ++i) {
00231 const Row& t_i = tableau[i];
00232 const Coefficient& t_ie = t_i[entering_var_index];
00233 const Coefficient& t_ib = t_i[base[i]];
00234 const int t_ie_sign = sgn(t_ie);
00235 if (t_ie_sign != 0 && t_ie_sign == sgn(t_ib)) {
00236 const Row& t_e = tableau[exiting_base_index];
00237 const Coefficient& t_ee = t_e[entering_var_index];
00238 lcm_assign(lcm, t_ee, t_ie);
00239 exact_div_assign(current_min, lcm, t_ee);
00240 current_min *= t_e[0];
00241 current_min = abs(current_min);
00242 exact_div_assign(challenger, lcm, t_ie);
00243 challenger *= t_i[0];
00244 challenger = abs(challenger);
00245 current_min -= challenger;
00246 const int sign = sgn(current_min);
00247 if (sign > 0
00248 || (sign == 0 && base[i] < base[exiting_base_index]))
00249 exiting_base_index = i;
00250 }
00251 }
00252 return exiting_base_index;
00253 }
00254
00255
00256
00257 bool
00258 PPL::LP_Problem::compute_simplex() {
00259 assert(tableau.num_columns() == working_cost.size());
00260 const dimension_type tableau_num_rows = tableau.num_rows();
00261 while (true) {
00262
00263 const dimension_type entering_var_index
00264 #if PPL_SIMPLEX_ENABLE_STEEPEST_EDGE
00265 = steepest_edge();
00266 #else
00267 = get_entering_var_index();
00268 #endif
00269
00270 if (entering_var_index == 0)
00271 return true;
00272
00273
00274 const dimension_type exiting_base_index
00275 = get_exiting_base_index(entering_var_index);
00276
00277 if (exiting_base_index == tableau_num_rows)
00278 return false;
00279
00280
00281
00282
00283 swap_base(entering_var_index, exiting_base_index);
00284 #if PPL_NOISY_SIMPLEX
00285 ++num_iterations;
00286 if (num_iterations % 200 == 0)
00287 std::cout << "Primal Simplex: iteration "
00288 << num_iterations << "." << std::endl;
00289 #endif
00290 }
00291 }
00292
00293
00294
00295 void
00296 PPL::LP_Problem::prepare_first_phase() {
00297
00298
00299
00300 const dimension_type tableau_old_n_cols = tableau.num_columns();
00301 for (dimension_type i = tableau.num_rows(); i-- > 0 ; ) {
00302 Row& tableau_i = tableau[i];
00303 if (tableau_i[0] > 0)
00304 for (dimension_type j = tableau_old_n_cols; j-- > 0; )
00305 neg_assign(tableau_i[j]);
00306 }
00307
00308
00309
00310
00311 if (tableau.max_num_columns() - tableau_old_n_cols <= tableau.num_rows())
00312 throw std::length_error("PPL::LP_Problem:\nthe maximum size of an "
00313 "internal data structure has been exceeded "
00314 "while solving the LP_Problem.");
00315 tableau.add_zero_columns(tableau.num_rows() + 1);
00316
00317 working_cost = Row(tableau.num_columns(), Row::Flags());
00318
00319
00320
00321
00322
00323 for (dimension_type i = 0; i < tableau.num_rows(); ++i) {
00324 const dimension_type j = tableau_old_n_cols + i;
00325 tableau[i][j] = 1;
00326 working_cost[j] = -1;
00327 base[i] = j;
00328 }
00329
00330
00331
00332 const dimension_type last_obj_index = working_cost.size() - 1;
00333 working_cost[last_obj_index] = 1;
00334
00335
00336 for (dimension_type i = tableau.num_rows(); i-- > 0; )
00337 linear_combine(working_cost, tableau[i], base[i]);
00338 }
00339
00340
00341
00342 void
00343 PPL::LP_Problem::erase_slacks() {
00344 const dimension_type tableau_last_index = tableau.num_columns() - 1;
00345 dimension_type tableau_n_rows = tableau.num_rows();
00346 const dimension_type first_slack_index = tableau_last_index - tableau_n_rows;
00347
00348
00349 for (dimension_type i = 0; i < tableau_n_rows; ++i)
00350 if (base[i] >= first_slack_index) {
00351
00352 Row& tableau_i = tableau[i];
00353 bool redundant = true;
00354 for (dimension_type j = first_slack_index; j-- > 1; )
00355 if (tableau_i[j] != 0) {
00356 swap_base(j, i);
00357 redundant = false;
00358 break;
00359 }
00360 if (redundant) {
00361
00362
00363 --tableau_n_rows;
00364 if (i < tableau_n_rows) {
00365
00366
00367 tableau_i.swap(tableau[tableau_n_rows]);
00368 base[i] = base[tableau_n_rows];
00369 --i;
00370 }
00371 tableau.erase_to_end(tableau_n_rows);
00372 base.pop_back();
00373 }
00374 }
00375
00376
00377
00378
00379 const dimension_type new_tableau_n_cols = first_slack_index + 1;
00380 const dimension_type new_tableau_last_index = first_slack_index;
00381
00382
00383 tableau.remove_trailing_columns(tableau.num_columns() - new_tableau_n_cols);
00384
00385 for (dimension_type i = tableau_n_rows; i-- > 0; )
00386 tableau[i][new_tableau_last_index] = 0;
00387
00388
00389
00390 working_cost[new_tableau_last_index] = working_cost[tableau_last_index];
00391
00392 const dimension_type working_cost_new_size = working_cost.size() -
00393 (tableau_last_index - new_tableau_last_index);
00394 working_cost.shrink(working_cost_new_size);
00395 }
00396
00397
00398
00399 PPL::LP_Problem_Status
00400 PPL::LP_Problem::compute_tableau() {
00401 assert(tableau.num_rows() == 0);
00402 assert(dim_map.size() == 0);
00403
00404
00405 Linear_System& cs = input_cs;
00406 const dimension_type cs_num_rows = cs.num_rows();
00407 const dimension_type cs_num_cols = cs.num_columns();
00408
00409
00410
00411
00412
00413
00414
00415
00416 dimension_type tableau_num_rows = cs_num_rows;
00417 dimension_type tableau_num_cols = 2*cs_num_cols - 1;
00418 dimension_type num_slack_variables = 0;
00419
00420
00421
00422
00423 std::deque<bool> is_tableau_constraint(cs_num_rows, true);
00424
00425
00426
00427 std::deque<bool> nonnegative_variable(cs_num_cols - 1, false);
00428
00429
00430 for (dimension_type i = cs_num_rows; i-- > 0; ) {
00431 const Linear_Row& cs_i = cs[i];
00432 bool found_a_nonzero_coeff = false;
00433 bool found_many_nonzero_coeffs = false;
00434 dimension_type nonzero_coeff_column_index = 0;
00435 for (dimension_type j = cs_num_cols; j-- > 1; ) {
00436 if (cs_i[j] != 0)
00437 if (found_a_nonzero_coeff) {
00438 found_many_nonzero_coeffs = true;
00439 if (cs_i.is_ray_or_point_or_inequality())
00440 ++num_slack_variables;
00441 break;
00442 }
00443 else {
00444 nonzero_coeff_column_index = j;
00445 found_a_nonzero_coeff = true;
00446 }
00447 }
00448
00449
00450 if (found_many_nonzero_coeffs)
00451 continue;
00452
00453 if (!found_a_nonzero_coeff) {
00454
00455
00456 if (cs_i.is_ray_or_point_or_inequality()) {
00457 if (cs_i[0] < 0)
00458
00459 return UNFEASIBLE_LP_PROBLEM;
00460 }
00461 else
00462
00463 if (cs_i[0] != 0)
00464
00465 return UNFEASIBLE_LP_PROBLEM;
00466
00467 is_tableau_constraint[i] = false;
00468 --tableau_num_rows;
00469 continue;
00470 }
00471 else {
00472
00473
00474
00475
00476
00477
00478
00479
00480
00481
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501 const dimension_type nonzero_var_index = nonzero_coeff_column_index - 1;
00502
00503 const int sgn_a = sgn(cs_i[nonzero_coeff_column_index]);
00504 const int sgn_b = sgn(cs_i[0]);
00505
00506 if (sgn_a == sgn_b) {
00507 if (cs_i.is_ray_or_point_or_inequality())
00508 ++num_slack_variables;
00509 }
00510
00511 else if (cs_i.is_line_or_equality()) {
00512 if (!nonnegative_variable[nonzero_var_index]) {
00513 nonnegative_variable[nonzero_var_index] = true;
00514 --tableau_num_cols;
00515 }
00516 }
00517
00518 else if (sgn_b < 0) {
00519 if (!nonnegative_variable[nonzero_var_index]) {
00520 nonnegative_variable[nonzero_var_index] = true;
00521 --tableau_num_cols;
00522 }
00523 ++num_slack_variables;
00524 }
00525
00526 else if (sgn_a > 0) {
00527 if (!nonnegative_variable[nonzero_var_index]) {
00528 nonnegative_variable[nonzero_var_index] = true;
00529 --tableau_num_cols;
00530 }
00531 is_tableau_constraint[i] = false;
00532 --tableau_num_rows;
00533 }
00534
00535 else
00536 ++num_slack_variables;
00537 }
00538 }
00539
00540
00541 tableau_num_cols += num_slack_variables;
00542
00543
00544 for (dimension_type i = 0, j = nonnegative_variable.size(),
00545 nnv_size = j; i < nnv_size; ++i)
00546 if (!nonnegative_variable[i]) {
00547 dim_map.insert(std::make_pair(i, j));
00548 ++j;
00549 }
00550
00551
00552
00553 if (tableau_num_rows > 0) {
00554 if (tableau_num_cols > tableau.max_num_columns())
00555 throw std::length_error("PPL::LP_Problem:\nthe maximum size of an "
00556 "internal data structure has been exceeded "
00557 "while solving the LP_Problem.");
00558 tableau.add_zero_rows_and_columns(tableau_num_rows,
00559 tableau_num_cols,
00560 Row::Flags());
00561 }
00562
00563
00564
00565
00566
00567
00568
00569
00570 for (dimension_type k = tableau_num_rows, slack_index = tableau_num_cols,
00571 i = cs_num_rows; i-- > 0; )
00572 if (is_tableau_constraint[i]) {
00573
00574 Row& tableau_k = tableau[--k];
00575 const Linear_Row& cs_i = cs[i];
00576 for (dimension_type j = cs_num_cols; j-- > 0; )
00577 tableau_k[j] = cs_i[j];
00578
00579 if (cs_i.is_ray_or_point_or_inequality())
00580 tableau_k[--slack_index] = -1;
00581 }
00582
00583
00584 typedef std::map<dimension_type, dimension_type>::const_iterator iter;
00585 for (iter map_itr = dim_map.begin(),
00586 map_end = dim_map.end(); map_itr != map_end; ++map_itr) {
00587 const dimension_type original_var = (map_itr->first) + 1;
00588 const dimension_type split_var = (map_itr->second) + 1;
00589 for (dimension_type i = tableau_num_rows; i-- > 0; ) {
00590 Row& tableau_i = tableau[i];
00591 tableau_i[split_var] = -tableau_i[original_var];
00592 }
00593 }
00594
00595
00596
00597
00598
00599 if (tableau_num_rows == 0)
00600 for (dimension_type i = tableau_num_cols; i-- > 1; )
00601 if (input_obj_function[i] > 0){
00602 status = UNBOUNDED;
00603 return UNBOUNDED_LP_PROBLEM;
00604 }
00605
00606
00607
00608 status = OPTIMIZED;
00609 return OPTIMIZED_LP_PROBLEM;
00610 }
00611
00612 bool
00613 PPL::LP_Problem::is_in_base(const dimension_type var_index,
00614 dimension_type& row_index) const {
00615 for (row_index = base.size(); row_index-- > 0; )
00616 if (base[row_index] == var_index)
00617 return true;
00618 return false;
00619 }
00620
00621 PPL::Generator
00622 PPL::LP_Problem::compute_generator() const {
00623
00624
00625 dimension_type original_space_dim = input_cs.space_dimension();
00626 std::vector<Coefficient> num(original_space_dim);
00627 std::vector<Coefficient> den(original_space_dim);
00628 dimension_type row = 0;
00629
00630
00631 typedef std::map<dimension_type, dimension_type>::const_iterator iter;
00632 iter map_end = dim_map.end();
00633
00634 for (dimension_type i = original_space_dim; i-- > 0; ) {
00635 Coefficient& num_i = num[i];
00636 Coefficient& den_i = den[i];
00637
00638
00639 if (is_in_base(i+1, row)) {
00640 const Row& t_row = tableau[row];
00641 if (t_row[i+1] > 0) {
00642 num_i= -t_row[0];
00643 den_i= t_row[i+1];
00644 }
00645 else {
00646 num_i= t_row[0];
00647 den_i= -t_row[i+1];
00648 }
00649 }
00650 else {
00651 num_i = 0;
00652 den_i = 1;
00653 }
00654
00655 iter map_iter = dim_map.find(i);
00656 if (map_iter != map_end) {
00657
00658
00659 const dimension_type split_i = map_iter->second;
00660
00661 if (is_in_base(split_i+1, row)) {
00662 const Row& t_row = tableau[row];
00663 TEMP_INTEGER(split_num);
00664 TEMP_INTEGER(split_den);
00665 if (t_row[split_i+1] > 0) {
00666 split_num = -t_row[0];
00667 split_den = t_row[split_i+1];
00668 }
00669 else {
00670 split_num = t_row[0];
00671 split_den = -t_row[split_i+1];
00672 }
00673
00674
00675 TEMP_INTEGER(lcm);
00676 lcm_assign(lcm, den_i, split_den);
00677 exact_div_assign(den_i, lcm, den_i);
00678 exact_div_assign(split_den, lcm, split_den);
00679 num_i *= den_i;
00680 sub_mul_assign(num_i, split_num, split_den);
00681 if (num_i == 0)
00682 den_i = 1;
00683 else
00684 den_i = lcm;
00685 }
00686
00687
00688 }
00689 }
00690
00691
00692 TEMP_INTEGER(lcm);
00693 lcm = den[0];
00694 for (dimension_type i = 1; i < original_space_dim; ++i)
00695 lcm_assign(lcm, lcm, den[i]);
00696
00697
00698 for (dimension_type i = original_space_dim; i-- > 0; ) {
00699 exact_div_assign(den[i], lcm, den[i]);
00700 num[i] *= den[i];
00701 }
00702
00703
00704 Linear_Expression expr;
00705 for (dimension_type i = original_space_dim; i-- > 0; )
00706 expr += num[i] * Variable(i);
00707 return point(expr, lcm);
00708 }
00709
00710 void
00711 PPL::LP_Problem::second_phase() {
00712
00713 assert(status == SATISFIABLE || status == UNBOUNDED || status == OPTIMIZED);
00714
00715 if (status == UNBOUNDED || status == OPTIMIZED)
00716 return;
00717
00718
00719 Row new_cost = input_obj_function;
00720 if (opt_mode == MINIMIZATION)
00721 for (dimension_type i = new_cost.size(); i-- > 0; )
00722 neg_assign(new_cost[i]);
00723
00724
00725 const dimension_type cost_zero_size = working_cost.size();
00726 Row tmp_cost = Row(new_cost, cost_zero_size, cost_zero_size);
00727 tmp_cost.swap(working_cost);
00728 working_cost[cost_zero_size-1] = 1;
00729
00730
00731 typedef std::map<dimension_type, dimension_type>::const_iterator iter;
00732 for (iter map_itr = dim_map.begin(),
00733 map_end = dim_map.end(); map_itr != map_end; ++map_itr){
00734 const dimension_type original_var = (map_itr->first) + 1;
00735 const dimension_type split_var = (map_itr->second) + 1;
00736 working_cost[split_var] = -working_cost[original_var];
00737 }
00738
00739
00740
00741 for (dimension_type i = tableau.num_rows(); i-- > 0; ) {
00742 const dimension_type base_i = base[i];
00743 if (working_cost[base_i] != 0)
00744 linear_combine(working_cost, tableau[i], base_i);
00745 }
00746
00747 bool second_phase_successful = compute_simplex();
00748
00749 #if PPL_NOISY_SIMPLEX
00750 std::cout << "LP_Problem::solve: 2nd phase ended at iteration "
00751 << num_iterations << "." << std::endl;
00752 #endif
00753 if (second_phase_successful) {
00754 last_generator = compute_generator();
00755 status = OPTIMIZED;
00756 }
00757 else
00758 status = UNBOUNDED;
00759 assert(OK());
00760 }
00761
00762 void
00763 PPL::LP_Problem::evaluate_objective_function(const Generator& evaluating_point,
00764 Coefficient& ext_n,
00765 Coefficient& ext_d) const {
00766 const dimension_type ep_space_dim = evaluating_point.space_dimension();
00767 if (space_dimension() < ep_space_dim)
00768 throw std::invalid_argument("PPL::LP_Problem::"
00769 "evaluate_objective_function(p, n, d):\n"
00770 "*this and p are dimension incompatible.");
00771 if (!evaluating_point.is_point())
00772 throw std::invalid_argument("PPL::LP_Problem::"
00773 "evaluate_objective_function(p, n, d):\n"
00774 "p is not a point.");
00775
00776
00777
00778 const dimension_type space_dim
00779 = std::min(ep_space_dim, input_obj_function.space_dimension());
00780
00781 ext_n = input_obj_function.inhomogeneous_term();
00782 for (dimension_type i = space_dim; i-- > 0; )
00783 ext_n += evaluating_point.coefficient(Variable(i))
00784 * input_obj_function.coefficient(Variable(i));
00785
00786 normalize2(ext_n, evaluating_point.divisor(), ext_n, ext_d);
00787 }
00788
00789 bool
00790 PPL::LP_Problem::is_satisfiable() const {
00791 #if PPL_NOISY_SIMPLEX
00792 num_iterations = 0;
00793 #endif
00794
00795 switch (status) {
00796 case UNSATISFIABLE:
00797 return false;
00798 case SATISFIABLE:
00799 return true;
00800 case UNBOUNDED:
00801 return true;
00802 case OPTIMIZED:
00803 return true;
00804 case PARTIALLY_SATISFIABLE:
00805 return false;
00806 case UNSOLVED:
00807 break;
00808 }
00809
00810 LP_Problem& x = const_cast<LP_Problem&>(*this);
00811
00812
00813
00814
00815
00816 const dimension_type space_dim = x.input_cs.num_columns() - 1;
00817
00818
00819 x.tableau.clear();
00820 x.dim_map.clear();
00821
00822 LP_Problem_Status s_status = x.compute_tableau();
00823
00824
00825 switch (s_status) {
00826 case UNFEASIBLE_LP_PROBLEM:
00827 return false;
00828 case UNBOUNDED_LP_PROBLEM:
00829
00830
00831 x.last_generator = point(0*Variable(space_dim-1));
00832 return true;
00833 case OPTIMIZED_LP_PROBLEM:
00834
00835
00836 if (x.tableau.num_rows() == 0) {
00837
00838 x.last_generator = point(0*Variable(space_dim-1));
00839 return true;
00840 }
00841 break;
00842 }
00843
00844 #if PPL_NOISY_SIMPLEX
00845 num_iterations = 0;
00846 #endif
00847
00848
00849 x.base = std::vector<dimension_type> (x.tableau.num_rows());
00850
00851
00852
00853 x.prepare_first_phase();
00854
00855 bool first_phase_successful = x.compute_simplex();
00856
00857 #if PPL_NOISY_SIMPLEX
00858 std::cout << "LP_Problem::solve: 1st phase ended at iteration "
00859 << num_iterations << "." << std::endl;
00860 #endif
00861
00862
00863 if (!first_phase_successful || x.working_cost[0] != 0){
00864 x.status = UNSATISFIABLE;
00865 return false;
00866 }
00867
00868
00869
00870
00871 x.last_generator = compute_generator();
00872 x.status = SATISFIABLE;
00873
00874 x.erase_slacks();
00875 return true;
00876 }
00877
00878 bool
00879 PPL::LP_Problem::OK() const {
00880 #ifndef NDEBUG
00881 using std::endl;
00882 using std::cerr;
00883 #endif
00884
00885
00886 if (input_cs.has_strict_inequalities()) {
00887 #ifndef NDEBUG
00888 cerr << "The feasible region of the LP_Problem is defined by "
00889 << "a constraint system containing strict inequalities."
00890 << endl;
00891 #endif
00892 return false;
00893 }
00894
00895
00896 const dimension_type space_dim = input_cs.space_dimension();
00897 if (space_dim < input_obj_function.space_dimension()) {
00898 #ifndef NDEBUG
00899 cerr << "The LP_Problem and the objective function have "
00900 << "incompatible space dimensions ("
00901 << space_dim << " < " << input_obj_function.space_dimension() << ")."
00902 << endl;
00903 #endif
00904 return false;
00905 }
00906
00907 if (status == SATISFIABLE || status == UNBOUNDED || status == OPTIMIZED) {
00908
00909
00910 if (space_dim != last_generator.space_dimension()) {
00911 #ifndef NDEBUG
00912 cerr << "The LP_Problem and the cached feasible point have "
00913 << "incompatible space dimensions ("
00914 << space_dim << " != " << last_generator.space_dimension() << ")."
00915 << endl;
00916 #endif
00917 return false;
00918 }
00919 if (!input_cs.satisfies_all_constraints(last_generator)) {
00920 #ifndef NDEBUG
00921 cerr << "The cached feasible point does not belong to "
00922 << "the feasible region of the LP_Problem."
00923 << endl;
00924 #endif
00925 return false;
00926 }
00927
00928 const dimension_type tableau_nrows = tableau.num_rows();
00929 const dimension_type tableau_ncols = tableau.num_columns();
00930
00931
00932 if (tableau_nrows != base.size()) {
00933 #ifndef NDEBUG
00934 cerr << "tableau and base have incompatible sizes" << endl;
00935 #endif
00936 return false;
00937 }
00938
00939
00940 if (tableau_ncols != working_cost.size()) {
00941 #ifndef NDEBUG
00942 cerr << "tableau and working_cost have incompatible sizes" << endl;
00943 #endif
00944 return false;
00945 }
00946
00947
00948 for (dimension_type i = base.size(); i-- > 0; )
00949 if (base[i] > tableau_ncols) {
00950 #ifndef NDEBUG
00951 cerr << "base contains an invalid column index" << endl;
00952 #endif
00953 return false;
00954 }
00955
00956
00957
00958 std::set<dimension_type> domain;
00959 std::set<dimension_type> range;
00960 typedef std::map<dimension_type, dimension_type>::const_iterator Iter;
00961 for (Iter i = dim_map.begin(), iend = dim_map.end(); i != iend; ++i) {
00962 domain.insert(i->first);
00963 range.insert(i->second);
00964 }
00965 if (domain.size() != range.size()
00966 || domain.end() != std::find_first_of(domain.begin(), domain.end(),
00967 range.begin(), range.end())) {
00968 #ifndef NDEBUG
00969 cerr << "dim_map encodes an invalid map" << endl;
00970 #endif
00971 return false;
00972 }
00973 }
00974
00975
00976
00977
00978 return true;
00979 }
00980
00981 void
00982 PPL::LP_Problem::ascii_dump(std::ostream& s) const {
00983 using namespace IO_Operators;
00984
00985 s << "input_cs\n";
00986 input_cs.ascii_dump(s);
00987 s << "\ninput_obj_function\n";
00988 input_obj_function.ascii_dump(s);
00989 s << "\nopt_mode " << (opt_mode == MAXIMIZATION ? "MAX" : "MIN") << "\n";
00990
00991 s << "\nstatus: ";
00992 switch (status) {
00993 case UNSATISFIABLE:
00994 s << "UNSAT";
00995 break;
00996 case SATISFIABLE:
00997 s << "SATIS";
00998 break;
00999 case UNBOUNDED:
01000 s << "UNBOU";
01001 break;
01002 case OPTIMIZED:
01003 s << "OPTIM";
01004 break;
01005 case PARTIALLY_SATISFIABLE:
01006 s << "P_SAT";
01007 break;
01008 case UNSOLVED:
01009 s << "UNSOL";
01010 break;
01011 }
01012 s << "\n";
01013
01014 s << "\ntableau\n";
01015 tableau.ascii_dump(s);
01016 s << "\nworking_cost\n";
01017 working_cost.ascii_dump(s);
01018
01019 const dimension_type base_size = base.size();
01020 s << "\nbase (" << base_size << ")\n";
01021 for (dimension_type i = 0; i != base_size; ++i)
01022 s << base[i] << ' ';
01023
01024 const dimension_type dim_map_size = dim_map.size();
01025 s << "\ndim_map (" << dim_map_size << ")\n";
01026 for (std::map<dimension_type, dimension_type>::const_iterator
01027 i = dim_map.begin(), iend = dim_map.end(); i != iend; ++i)
01028 s << Variable(i->first) << "->" << Variable(i->second) << ' ';
01029
01030 s << "\nlast_generator\n";
01031 last_generator.ascii_dump(s);
01032 s << "\n";
01033 }
01034
01035 PPL_OUTPUT_DEFINITIONS_ASCII_ONLY(LP_Problem);